文章参考了ZhaoDongqiang（C++调用cplex求解VRPTW模型），补充了在Jyputer编辑器中用Python调用Cplex求解VRPTW问题，修正了模型中的小错误，具体如下：

目录

1 VRPTW数学模型

2 Python调用Cplex求解

2.1 调用所需的库 

2.2  初始化参数

2.3 调用docplex建模

2.4 求解结果

 3 结语

1 VRPTW数学模型

2 Python调用Cplex求解

2.1 调用所需的库 

import pandas as pd
import numpy as np
from docplex.mp.model import Model
import matplotlib.pyplot as plt

2.2  初始化参数

 data.txt 内容如下：

40	50	0	0	1236	0
45	68	10	912	967	90
45	70	30	825	870	90
42	66	10	65	146	90
42	68	10	727	782	90
42	65	10	15	67	90
40	69	20	621	702	90
40	66	20	170	225	90
38	68	20	255	324	90
38	70	10	534	605	90
35	66	10	357	410	90
35	69	10	448	505	90
25	85	20	652	721	90
22	75	30	30	92	90
22	85	10	567	620	90
20	80	40	384	429	90
20	85	40	475	528	90
18	75	20	99	148	90
15	75	20	179	254	90
15	80	10	278	345	90
30	50	10	10	73	90
30	52	20	914	965	90
28	52	20	812	883	90
28	55	10	732	777	90
25	50	10	65	144	90
25	52	40	169	224	90
40	50	0	0	1236	0

df = pd.read_csv('data.txt',sep='\t',header=None)
M=10000
c_num = 25 #顾客数
v_num = 3  #卡车数
q = 500    #卡车的容量
CM = [i for i in range (1,c_num+1)]  #顾客列表
V = [i for i in range(1,v_num+1)]   #卡车列表
N = [0]+CM+[26]   #结点列表

#读取txt文件
x_loc = list(df.iloc[:,0])
y_loc = list(df.iloc[:,1])
d = list(df.iloc[:,2])
lbtw = list(df.iloc[:,3])
ubtw = list(df.iloc[:,4])
t_service = list(df.iloc[:,5])
t_service[0]

#各个节点之间的弧
A = [(i,j) for i in N for j in N if j!=i]

#各节点之间的距离
C={(i,j):np.floor(np.hypot((x_loc[i]-x_loc[j]),(y_loc[i]-y_loc[j]))*10)/10 for i,j in A}

#将两点距离作为时间
t={(i,j):0 for i,j in A}
for i,j in C:
    t[i,j]=C[i,j]+t_service[i]

#画出坐标散点图
import matplotlib.pyplot as plt
plt.scatter(x_loc[1:25],y_loc[1:25],c='b')
plt.plot(x_loc[0],y_loc[0],c='r',marker='s')
for i in CM:
    plt.annotate('%d'%i,(x_loc[i]+0.5,y_loc[i]))    
plt.annotate("%d"%0,(x_loc[0]+0.5,y_loc[0]))


#初始化A_car 和 S_time，作为决策变量s和x字典的键值 
A_car=[(k,i,j) for k in V for i,j in A ]
S_time=[(k,i) for k in V for i in N]

 图1 各点的坐标散点图

2.3 调用docplex建模

mdl=Model("VRPTW")

#决策变量
x=mdl.binary_var_dict(A_car,name='x')
s=mdl.continuous_var_dict(S_time,name='s')

#目标函数
mdl.minimize(mdl.sum(C[i,j]*x[k,i,j] for k,i,j in A_car))

#约束1-8：
#约束1
mdl.add_constraints(mdl.sum(x[k,i,j] for k in V for j in N if j!=i\
                           )==1 for i in CM)
#约束2
mdl.add_constraints(mdl.sum(d[i]*mdl.sum(x[k,i,j] for j in N if j!=i)\
                            for i in CM)<=q for k in V)
#约束3
mdl.add_constraints(mdl.sum(x[k,0,j] for j in N if j!=0)==1 for k in V)

#约束4
mdl.add_constraints(mdl.sum(x[k,i,h] for i in N if i!=h)==mdl.sum(x[k,h,j] for j in N if j!=h) for h in CM for k in V)

#约束5
mdl.add_constraints(mdl.sum(x[k,i,c_num+1] for i in N if i!=(c_num+1))==1 for k in V)
#约束6
mdl.add_constraint(mdl.sum(x[k,0,j] for k in V for j in N if j!=0)<=v_num)
#约束7

mdl.add_constraints(s[k,i]+t[i,j]-M*(1-x[k,i,j])<=s[k,j] for k in V for i,j in A)

#约束8
mdl.add_constraints(s[k,i]<=ubtw[i] for k in V for i in N)
mdl.add_constraints(s[k,i]>=lbtw[i] for k in V for i in N);

#求解
solution=mdl.solve()
print(solution)

solution for: VRPTW
objective: 191.3
x_1_0_20=1
x_1_20_24=1
x_1_21_26=1
x_1_22_21=1
x_1_23_22=1
x_1_24_25=1
x_1_25_23=1
x_2_0_13=1
x_2_12_26=1
x_2_13_17=1
x_2_14_12=1
x_2_15_16=1
x_2_16_14=1
x_2_17_18=1
x_2_18_19=1
x_2_19_15=1
x_3_0_5=1
x_3_1_26=1
x_3_2_1=1
x_3_3_7=1
x_3_4_2=1
x_3_5_3=1
x_3_6_4=1
x_3_7_8=1
x_3_8_10=1
x_3_9_6=1
x_3_10_11=1
x_3_11_9=1
s_1_1=967.000
s_1_2=825.000
s_1_3=65.000
s_1_4=727.000
s_1_5=15.000
s_1_6=621.000
s_1_7=170.000
s_1_8=255.000
s_1_9=534.000
s_1_10=357.000
s_1_11=448.000
s_1_12=652.000
s_1_13=30.000
s_1_14=567.000
s_1_15=384.000
s_1_16=475.000
s_1_17=99.000
s_1_18=179.000
s_1_19=278.000
s_1_20=10.000
s_1_21=917.000
s_1_22=825.000
s_1_23=732.000
s_1_24=132.000
s_1_25=224.000
s_1_26=1075.600
s_2_1=912.000
s_2_2=825.000
s_2_3=65.000
s_2_4=727.000
s_2_5=15.000
s_2_6=621.000
s_2_7=170.000
s_2_8=255.000
s_2_9=534.000
s_2_10=357.000
s_2_11=448.000
s_2_12=687.800
s_2_13=30.800
s_2_14=594.800
s_2_15=407.800
s_2_16=502.800
s_2_17=124.800
s_2_18=217.800
s_2_19=312.800
s_2_20=10.000
s_2_21=914.000
s_2_22=812.000
s_2_23=732.000
s_2_24=65.000
s_2_25=169.000
s_2_26=815.800
s_3_1=940.600
s_3_2=848.600
s_3_3=106.100
s_3_4=755.000
s_3_5=15.100
s_3_6=662.800
s_3_7=198.100
s_3_8=290.900
s_3_9=570.600
s_3_10=384.500
s_3_11=477.500
s_3_12=652.000
s_3_13=30.000
s_3_14=567.000
s_3_15=384.000
s_3_16=475.000
s_3_17=99.000
s_3_18=179.000
s_3_19=278.000
s_3_20=10.000
s_3_21=914.000
s_3_22=812.000
s_3_23=732.000
s_3_24=65.000
s_3_25=169.000
s_3_26=1049.200

2.4 求解结果

#选出那些决策变量为1时所经过的弧
active_arcs=[ a for a in A_car if x[a].solution_value > 0.9]

#画出VRPTW的结果
import matplotlib.pyplot as plt
plt.scatter(x_loc[1:25],y_loc[1:25])
plt.plot(x_loc[0],y_loc[0],c='r',marker='s')
for i in CM:
    plt.annotate('%d'%i,(x_loc[i]+0.5,y_loc[i]))
    
    
plt.annotate("%d"%0,(x_loc[0]+0.5,y_loc[0]))
for k,i,j in active_arcs:
    if k==1:
        plt.plot((x_loc[i],x_loc[j]),(y_loc[i],y_loc[j]),c='r')
    if k==2:
        plt.plot((x_loc[i],x_loc[j]),(y_loc[i],y_loc[j]),c='g')
    if k==3:
        plt.plot((x_loc[i],x_loc[j]),(y_loc[i],y_loc[j]),c='b')
        

 3 结语

目前Python调用Cplex的相关文章较少，本文在ZhaoDongqiang（C++调用cplex求解VRPTW模型）的基础之上，利用Python调用Cplex解决了VRPTW问题。