问题描述

以一个经典的选址问题为例，建立数学模型，假设代表4s店的集合，代表配送中心的侯选地点的集合，而弧长则代表可能的配送中心地点与需求节点之间的物流量。假设有个起始点（4s店）和个配送中心。
变量
:表示4s店i运送到配送中心j的数量；
:表示j点是否设为配送中心；
:表示i到j的运输费用，通过运输距离乘以单价得到；
:表示j点建设配送中心的固定成本;
:表示j点建造配送中心的容量；
：表示配送中心的建设数量；

模型

公式蛮难打的，我就截图了hh

image.png

程序

数据部分

from docplex.mp.model import Model  #调用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

# 读取和处理起始点（客户）数据
filename = '/Users/chenyanbo/Documents/python/cplextest/docplex/test.xlsx'
book = xlrd.open_workbook(filename)
sheet = book.sheet_by_name('Sheet1')
#['序号', '起始地-目的地', '纬度', '经度', '参考距离\n（km）', '汇总']
city_o = city(sheet.col_values(1),sheet.col_values(2),sheet.col_values(3))
demand = sheet.col_values(5)

# 读取和处理终点（候选配送中心）数据
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:])

# 计算起点到终点的距离矩阵 这是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)
            )       
  # 目标函数：最大化利润

# 添加约束条件
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("总费用为: %g" %total_cost)
print("被选为配送中心的城市为:")
for f_loc in range(0,len(open_facility)): 
   print(city_d.name[open_facility[f_loc]])

结果输出

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
总费用为: 1.64679e+08
被选为配送中心的城市为:
重庆
天津
深圳
上海
长沙

输出数据处理和可视化

# 输出运输量
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

选址结果.png

小结

docplex真好用，现在基本上更新到cplex129以后就支持python3.7了，对于约束和变量的写法都比较轻松，类似于matlab里的Yalmip工具。