問題描述
以一個經(jīng)典的選址問題為例芭届，建立數(shù)學(xué)模型抹蚀，假設(shè)代表4s店的集合植影，
代表配送中心的侯選地點(diǎn)的集合丰涉，而弧長
則代表可能的配送中心地點(diǎn)與需求節(jié)點(diǎn)之間的物流量拓巧。假設(shè)有
個起始點(diǎn)（4s店）和
個配送中心。
變量
:表示4s店i運(yùn)送到配送中心j的數(shù)量一死；
:表示j點(diǎn)是否設(shè)為配送中心肛度；
:表示i到j(luò)的運(yùn)輸費(fèi)用，通過運(yùn)輸距離乘以單價得到摘符；
:表示j點(diǎn)建設(shè)配送中心的固定成本;
:表示j點(diǎn)建造配送中心的容量贤斜；
：表示配送中心的建設(shè)數(shù)量；
模型
公式蠻難打的逛裤，我就截圖了hh
image.png
程序
數(shù)據(jù)部分
from docplex.mp.model import Model #調(diào)用cplex solver
import numpy as np
import folium
import xlrd
from geopy.distance import great_circle
def get_distance(p1, p2):
return great_circle(p1,p2).meters
class city:
def __init__(self, name, x, y):
self.x = x
self.y = y
self.name = name
# 讀取和處理起始點(diǎn)（客戶）數(shù)據(jù)
filename = '/Users/chenyanbo/Documents/python/cplextest/docplex/test.xlsx'
book = xlrd.open_workbook(filename)
sheet = book.sheet_by_name('Sheet1')
#['序號', '起始地-目的地', '緯度', '經(jīng)度', '參考距離\n（km）', '匯總']
city_o = city(sheet.col_values(1),sheet.col_values(2),sheet.col_values(3))
demand = sheet.col_values(5)
# 讀取和處理終點(diǎn)（候選配送中心）數(shù)據(jù)
filename = '/Users/chenyanbo/Documents/python/cplextest/docplex/test2.xlsx'
book = xlrd.open_workbook(filename)
sheet = book.sheet_by_name('Sheet1')
city_d = city(sheet.col_values(0)[1:],sheet.col_values(1)[1:],sheet.col_values(2)[1:])
# 計算起點(diǎn)到終點(diǎn)的距離矩陣 這是20行35列
num_city_o = len(city_o.x)
num_city_d = len(city_d.x)
distance = np.zeros((num_city_o,num_city_d))
for i in range(0,num_city_o):
for j in range(0,num_city_d):
distance[i][j] = get_distance([city_o.x[i],city_o.y[i]],[city_d.x[j],city_d.y[j]])
模型部分
# 模型的建立和求解
mdl = Model(name='test model') # 初始化模型
customer = [i for i in range(0,num_city_o)]
facility = [j for j in range(0,num_city_d)]
y = mdl.binary_var_list(facility,name='is_open')
x = mdl.continuous_var_matrix(customer,facility,name = 'serviced_amount')
d = 4 * demand
f = 500000 * np.ones(num_city_d)
M = 1000 * np.ones(num_city_d)
# c = [[4,5,6,8,10],[6,4,3,5,8],[9,7,4,3,4]]
# c = np.transpose(np.array([cost,cost]))
price = 10
c = price * distance
mdl.minimize(
mdl.sum(f[j]*y[j] for j in facility)
+
mdl.sum(c[i,j]*x[i,j] for i in customer for j in facility)
)
# 目標(biāo)函數(shù)：最大化利潤
# 添加約束條件
mdl.add_constraints(mdl.sum(x[i,j] for j in facility) == d[i] for i in customer)
mdl.add_constraints(mdl.sum(x[i,j] for i in customer) <= M[j]*y[j] for j in facility)
mdl.add_constraints(x[i,j] >= 0 for i in customer for j in facility)
mdl.add_constraint(mdl.sum(y[j] for j in facility) == 5)
# mdl.print_information()
sol = mdl.solve() # 求解模型
mdl.print_solution()
total_cost = mdl.objective_value
open_facility = [f_loc for f_loc in range(0,num_city_d) if y[f_loc].solution_value == 1]
print("總費(fèi)用為: %g" %total_cost)
print("被選為配送中心的城市為:")
for f_loc in range(0,len(open_facility)):
print(city_d.name[open_facility[f_loc]])
結(jié)果輸出
objective: 164679104.574
is_open_0=1
is_open_9=1
is_open_13=1
is_open_15=1
is_open_33=1
serviced_amount_0_15=34.000
serviced_amount_1_15=3.000
serviced_amount_2_9=3.000
serviced_amount_3_9=1.000
serviced_amount_4_9=1.000
serviced_amount_5_15=1.000
serviced_amount_6_15=1.000
serviced_amount_7_15=1.000
serviced_amount_9_9=1.000
serviced_amount_10_15=1.000
serviced_amount_14_9=1.000
serviced_amount_15_9=1.000
serviced_amount_16_9=1.000
serviced_amount_17_9=1.000
serviced_amount_20_13=13.000
serviced_amount_21_13=2.000
serviced_amount_22_13=1.000
serviced_amount_23_13=2.000
serviced_amount_24_13=1.000
serviced_amount_25_13=1.000
serviced_amount_26_9=2.000
serviced_amount_27_9=1.000
serviced_amount_28_33=2.000
serviced_amount_29_0=2.000
serviced_amount_30_0=1.000
serviced_amount_31_33=2.000
serviced_amount_32_0=2.000
serviced_amount_33_0=1.000
serviced_amount_34_33=1.000
serviced_amount_35_9=1.000
serviced_amount_36_0=1.000
serviced_amount_38_0=1.000
serviced_amount_40_9=1.000
serviced_amount_41_0=1.000
總費(fèi)用為: 1.64679e+08
被選為配送中心的城市為:
重慶
天津
深圳
上海
長沙
輸出數(shù)據(jù)處理和可視化
# 輸出運(yùn)輸量
cf = [[i,j] for i in range(0,num_city_o) for j in range(0,num_city_d) if x[i,j].solution_value != 0]
cf_values = list(sol.get_value_dict(x).values())
cf_values = [cf_values[i] for i in range(0,len(cf_values)) if cf_values[i] != 0]
cf = np.column_stack((cf,cf_values))
# 繪圖部分
map_osm = folium.Map(location=[34,108], zoom_start=4.25)
for i in range(0,num_city_o):
folium.Marker([city_o.x[i],city_o.y[i]]).add_to(map_osm)
for j in range(0,len(open_facility)):
folium.Marker([city_d.x[open_facility[j]],city_d.y[open_facility[j]]]).add_to(map_osm)
for k in range(0,len(cf)):
cf = cf.astype(int)
map_osm.add_child(folium.PolyLine([[city_o.x[cf[k,0]],city_o.y[cf[k,0]]],[city_d.x[cf[k,1]],city_d.y[cf[k,1]]]], color=' #FF6A6A', weight=6))
map_osm
選址結(jié)果.png
小結(jié)
docplex真好用瘩绒，現(xiàn)在基本上更新到cplex129以后就支持python3.7了，對于約束和變量的寫法都比較輕松带族，類似于matlab里的Yalmip工具锁荔。
