問題描述

以一個(gè)經(jīng)典的選址問題為例，建立數(shù)學(xué)模型，假設(shè)代表4s店的集合，代表配送中心的侯選地點(diǎn)的集合，而弧長則代表可能的配送中心地點(diǎn)與需求節(jié)點(diǎn)之間的物流量。假設(shè)有個(gè)起始點(diǎn)（4s店）和個(gè)配送中心。
變量
:表示4s店i運(yùn)送到配送中心j的數(shù)量；
:表示j點(diǎn)是否設(shè)為配送中心；
:表示i到j(luò)的運(yùn)輸費(fèi)用，通過運(yùn)輸距離乘以單價(jià)得到；
:表示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')
#['序號(hào)', '起始地-目的地', '緯度', '經(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:])

# 計(jì)算起點(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了，對(duì)于約束和變量的寫法都比較輕松，類似于matlab里的Yalmip工具。