目录
- 1. 引言
- 2. 求解器介绍
- 3. 基础语言
- 3.1 创建模型
- 3.2 添加变量
- 3.3 添加目标函数
- 3.4 添加约束
- 3.5 设置参数
- 3.6 求解
- 4. 数学模型
- 4.1 [CVRP数学模型](https://mp.weixin.qq.com/s/DYh-5WkrYxk1gCKo8ZjvAw)
- 4.2 [VRPTW数学模型](https://mp.weixin.qq.com/s/tF-ayzjpZfuZvelvItuecw)
- 5. 完整代码
- 5.1 Python调用Gurobi求解CVRP
- 5.2 Python调用Gurobi求解VRPTW
- 5.3 Python调用COPT求解CVRP
- 5.4 Python调用COPT求解VRPTW
- 5.5 Python调用SCIP求解CVRP
- 5.6 Python调用SCIP求解VRPTW
- 7. 测试案例
- 8. 测试参数
- 9. 测试结果
- 9.1 CVRP求解结果汇总
- 9.2 VRPTW求解结果汇总
- 9.3 上界下降曲线对比（以CVRP为例）
- 9.5 车辆路径可视化（以CVRP为例）
- c101-31（Gurobi）
- c201-31（Gurobi）
- r101-31（Gurobi）
- 10. 小节
欢迎关注个人微信公众号：Python助力交通
1. 引言
优化求解器是解决复杂工程问题不可或缺的工具，可以帮助我们验证模型的正确性、理解决策变量的耦合关系、获取最优决策方案（合适规模条件下）。小编搜罗了网上关于各类常见（其实并不常见）的优化求解器介绍的帖子：
- 优化求解器资源盘点
- 干货 | 运筹学、数学规划、离散优化求解器大PK，总有一款适合你
- 【学界】运筹学数学规划|离散优化求解器大搜罗
- 除了Gurobi /SCIP，国内外还有哪些优化求解器？
- 开源线性规划求解器（Linear Programming solver）LP_Solve和CLP的PK
- 有哪些简便好用的解凸优化的工具箱或者包？
- 开源建模框架+开源求解器 | 使用pyomo建模框架实现交通物流优化问题的求解
除了以上求解器外，还有一些针对特定问题而量身定制高效率求解器：
- VRP问题求解器
- 基于求解器的路径规划算法实现及性能分析
- 智能优化算法工具包scikit-opt介绍及vrp问题求解评测（一）
- Python主要智能优化算法库汇总
今天的主题是以CVRP和VRPTW问题为例，分享Python调用运筹优化求解器（Gurobi、COPT、SCIP）的教程。
2. 求解器介绍
（1）Gurobi
Gurobi是由美国 Gurobi Optimization 公司开发新一代大规模优化器，提供 C++, Java, Python, .Net, Matlab 和R等多种接口，支持求解以下模型：
（1）连续和混合整数线性问题
（2）凸目标或约束连续和混合整数二次问题
（3）非凸目标或约束连续和混合整数二次问题
（4）含有对数、指数、三角函数、高阶多项式目标或约束，以及任何形式的分段约束的非线性问题
（5）含有绝对值、最大值、最小值、逻辑与或非目标或约束的非线性问题
（2）COPT
杉数求解器COPT（Cardinal Optimizer）是杉数自主研发的针对大规模优化问题的高效数学规划求解器套件，也是支撑杉数端到端供应链平台的核心组件，是目前同时具备大规模混合整数规划、线性规划（单纯形法和内点法）、半定规划、（混合整数）二阶锥规划以及（混合整数）凸二次规划和（混合整数）凸二次约束规划问题求解能力的综合性能数学规划求解器，为企业应对高性能求解的需求提供了更多选择。COPT支持所有主流编程接口：C、C++、C#、Python、Julia、Java、AMPL、GAMS、Pyomo、PuLP、CVXPY。
（3）SCIP
SCIP起源于ZIB(Zuse Institute Berlin)，由Tobias Achterberg奠定整个框架，是目前用于混合整数规划（MIP）和混合整数非线性规划（MINLP）的最快的非商业解算器之一，它允许用户对求解过程进行全面控制，支持自定义搜索树中的各个模块，在分支限界(Branch and Bound)过程中添加变量等功能。SCIP支持Python Java AMPL GAMS MATLAB等编程语言。
3. 基础语言
3.1 创建模型
# gurobi
cvrp_model = Model('cvrp')
# copt
env = Envr()
cvrp_model = env.createModel('cvrp')
# scip
cvrp_model = Model('cvrp')
3.2 添加变量
# gurobi
x = cvrp_model.addVar(vtype=GRB.BINARY, name='x') # 单个变量
X = cvrp_model.addVars(N, N, K, vtype=GRB.BINARY, name='X[i,j,k]') # 多个变量
# copt
x = cvrp_model.addVars(vtype=COPT.BINARY, name='x') # 单个变量
X = cvrp_model.addVars(N, N, K, vtype=COPT.BINARY, nameprefix='X[i,j,k]') # 多个变量
# scip
X[i,j,k] = cvrp_model.addVar(vtype="B", name=f"x({i},{j},{k})") # 单个变量(貌似只能添加单个变量)
3.3 添加目标函数
# gurobi
cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K), GRB.MINIMIZE)
# copt
cvrp_model.setObjective(quicksum(X[i, j, k] * cost[i, j] for i in N for j in N for k in K), sense=COPT.MINIMIZE)
# scip
cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K),'minimize' )
3.4 添加约束
# gurobi
cvrp_model.addConstr() # 单个约束
cvrp_model.addConstrs() # 多个约束
# copt
cvrp_model.addConstr() # 单个约束
cvrp_model.addConstrs() # 多个约束
# scip
cvrp_model.addCons() # 单个约束
cvrp_model.addConss() # 多个约束
3.5 设置参数
# gurobi
cvrp_model.setParam(GRB.Param.LogFile, './gurobi_r101-31.log') # 保存日志
cvrp_model.Params.TimeLimit = 1200 # 设置求解时间
# copt
cvrp_model.setLogFile('./copt_r101-31.log') # 保存日志
cvrp_model.param.timelimit = 1200 # 设置求解时间
# scip
cvrp_model.setLogfile('./scip_r101-31.log') # 保存日志
cvrp_model.setRealParam('limits/time', 1200) # 设置求解时间
3.6 求解
# gurobi
cvrp_model.optimize()
# copt
cvrp_model.solve()
# scip
cvrp_model.optimize()
4. 数学模型
4.1 CVRP数学模型
4.2 VRPTW数学模型
5. 完整代码
5.1 Python调用Gurobi求解CVRP
import xlsxwriter
import math
import pandas as pd
import matplotlib.pyplot as plt
from gurobipy import Model,quicksum,GRB
def read_input(filename):
"""
:param filename: 数据文件路径
:return:
"""
df = pd.read_csv(filename)
x_coord = { df['id'][i]:df['x_coord'][i] for i in range(df.shape[0]) } # 节点横坐标
y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) } # 节点纵坐标
demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) } # 节点需求
cost = {}
N = df['id'].tolist()
for f_n in N:
for t_n in N:
dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
cost[f_n,t_n] = dist
return N,cost,demand,x_coord,y_coord
def build_model(N,K,depot,CAP,cost,demand):
"""
:param N: 网络节点集合
:param K: 车队集合
:param depot: 配送中心id
:param CAP: 车辆容量
:param cost: 网络弧费用集合
:param demand: 网络节点需求集合
:return:
"""
cvrp_model = Model('cvrp')
# 添加变量
X = cvrp_model.addVars(N, N, K, vtype=GRB.BINARY, name='X[i,j,k]')
Y = cvrp_model.addVars(N, K, vtype=GRB.BINARY, name='Y[i,k]')
U = cvrp_model.addVars(N, K, vtype=GRB.INTEGER, name='U[i,k]')
# 设置目标函数
cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K), GRB.MINIMIZE)
# 添加约束
# 需求覆盖约束
cvrp_model.addConstrs( quicksum(Y[i,k] for k in K) == 1 for i in N[1:] )
# 车辆启用约束
cvrp_model.addConstr( quicksum(Y[depot,k] for k in K) == len(K) )
# 车辆流平衡约束
cvrp_model.addConstrs( quicksum(X[i,j,k] for j in N ) == quicksum(X[j,i,k] for j in N ) for i in N for k in K )
# 车辆路径限制
cvrp_model.addConstrs( quicksum(X[i,j,k] for j in N) == Y[i,k] for i in N for k in K )
# 车辆容量约束
cvrp_model.addConstrs( quicksum(Y[i,k] * demand[i] for i in N) <= CAP for k in K )
# 破圈约束
cvrp_model.addConstrs( U[i,k] - U[j,k] + CAP * X[i,j,k] <= CAP - demand[i] for i in N[1:] for j in N[1:] for k in K )
cvrp_model.addConstrs( U[i,k] <= CAP for i in N[1:] for k in K)
cvrp_model.addConstrs( U[i,k] >= demand[i] for i in N[1:] for k in K )
return cvrp_model,X,Y,U
def draw_routes(route_list,x_coord,y_coord):
for route in route_list:
path_x = []
path_y = []
for n in route:
path_x.append(x_coord[n])
path_y.append(y_coord[n])
plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
plt.show()
def save_file(route_list,total_cost):
wb = xlsxwriter.Workbook('路径方案.xlsx')
ws = wb.add_worksheet()
ws.write(0,0,'总费用')
ws.write(0,1,total_cost)
ws.write(1,0,'车辆')
ws.write(1,1,'路径')
row = 2
for route in route_list:
ws.write(row,0,route[0])
route_str = [str(i) for i in route[1:]]
ws.write(row,1,'-'.join(route_str))
row += 1
wb.close()
if __name__ == '__main__':
filename = './datasets/CVRP/r101-31.csv'
N, cost, demand, x_coord, y_coord = read_input(filename)
depot = N[0]
K = list(range(1,10))
CAP = 80
cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)
cvrp_model.setParam(GRB.Param.LogFile, './gurobi_r101-31.log')
cvrp_model.Params.TimeLimit = 1200
cvrp_model.optimize()
route_list = []
for k in K:
route = [depot]
cur_node = depot
cur_k = None
for j in N[1:]:
if X[depot, j, k].x > 0:
cur_node = j
cur_k = k
route.append(j)
N.remove(j)
break
if cur_k is None:
continue
while cur_node != depot:
for j in N:
if X[cur_node, j, cur_k].x > 0:
cur_node = j
route.append(j)
if j != depot:
N.remove(j)
break
route_list.append(route)
print("最优路径距离:", cvrp_model.objVal)
print("最优路径使用车辆数:", len(route_list))
draw_routes(route_list, x_coord, y_coord)
save_file(route_list, cvrp_model.objVal)
5.2 Python调用Gurobi求解VRPTW
import math
import pandas as pd
import matplotlib.pyplot as plt
import xlsxwriter
from gurobipy import GRB,Model,quicksum,tupledict,tuplelist
def read_input(filename):
"""
:param filename: 数据文件路径
:return:
"""
N = [] #所有节点
Q = {} #节点需求
TT = {} #节点旅行时间
ET = {} #节点最早开始服务时间
LT = {} #节点最晚结束服务时间
ST = {} #节点服务时间
x_coord = {} # 节点横坐标
y_coord = {} # 节点纵坐标
Cost={}
df = pd.read_csv(filename)
for i in range(df.shape[0]):
id = df['id'][i]
N.append(id)
x_coord[id] = df['x_coord'][i]
y_coord[id] = df['y_coord'][i]
Q[id] = df['demand'][i]
ET[id] = df['start_time'][i]
LT[id] = df['end_time'][i]
ST[id] = df['service_time'][i]
for f_n in N:
for t_n in N:
dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
Cost[f_n,t_n] = dist
TT[f_n,t_n] = dist/1 # 假设速度为1
return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coord
def build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K):
"""
:param N: 所有节点集合，其中N[0]为车场
:param Q: 节点需求集合
:param TT: 旅行时间
:param ET: 节点最早开始服务时间
:param LT：节点最晚结束服务时间
:param ST: 节点服务时间
:param CAP: 车辆容量
:param Cost: 旅行费用
:param K: 车队
:return:
"""
C = N[1:] #需求节点
M=10**5
depot = N[0]
# 创建模型
vrptw_model=Model()
# 添加变量
X = vrptw_model.addVars(N,N,K,vtype=GRB.BINARY,name='X(i,j,k)')
T = vrptw_model.addVars( N,K,vtype=GRB.CONTINUOUS,lb=0,name='T[i,k]')
# 设置目标函数
z1 = quicksum( Cost[i,j]*X[i,j,k] for i in N for j in N for k in K if i!=j)
vrptw_model.setObjective(z1,GRB.MINIMIZE)
# 车辆起点约束
vrptw_model.addConstrs(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
# 车辆路径连续约束
vrptw_model.addConstrs( quicksum(X[i,j,k] for j in N if j!=i)==quicksum(X[j,i,k] for j in N if j!=i) for i in C for k in K)
# 车辆终点约束
vrptw_model.addConstrs(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
# 需求服务约束
vrptw_model.addConstrs( quicksum(X[i,j,k] for k in K for j in N if j!=i)==1 for i in C)
# 车辆容量约束
vrptw_model.addConstrs( quicksum(Q[i]*X[i,j,k] for i in C for j in N if i!=j)<=CAP for k in K )
# 时间窗约束
vrptw_model.addConstrs( T[i,k]+ST[i]+TT[i,j]-(1-X[i,j,k])*M<=T[j,k] for i in C for j in C for k in K if i!=j )
vrptw_model.addConstrs( T[i,k] >= ET[i] for i in N for k in K)
vrptw_model.addConstrs( T[i,k] <= LT[i] for i in N for k in K)
return vrptw_model,X,T
def draw_routes(route_list,x_coord,y_coord):
for route in route_list:
path_x = []
path_y = []
for n in route:
path_x.append(x_coord[n])
path_y.append(y_coord[n])
plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
plt.show()
def save_file(route_list,route_timetable,total_cost):
wb = xlsxwriter.Workbook('路径方案.xlsx')
ws = wb.add_worksheet()
ws.write(0,0,'总费用')
ws.write(0,1,total_cost)
ws.write(1,0,'车辆')
ws.write(1,1,'路径')
ws.write(1,2,'时刻')
row = 2
for id,route in enumerate(route_list):
ws.write(row,0,route[0])
route_str = [str(i) for i in route[1:]]
ws.write(row,1,'-'.join(route_str))
timetable_str = [str(i) for i in route_timetable[id]]
ws.write(row,2,'-'.join(timetable_str))
row += 1
wb.close()
if __name__=='__main__':
filename = './data/r101.csv'
N, Q, TT, ET, LT, ST, Cost,x_coord,y_coord = read_input(filename)
depot = N[0]
K = list(range(1, 30))
CAP = 80
vrptw_model,X,T = build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K)
vrptw_model.setParam(GRB.Param.LogFile, './log/gurobi_r101.log')
vrptw_model.Params.TimeLimit = 1500
vrptw_model.optimize()
route_list = []
route_timetable = []
# 提取车辆路径
for k in K:
route = [depot]
cur_node = depot
cur_k = None
for j in N[1:]:
if X[depot, j, k].x > 0:
cur_node = j
cur_k = k
route.append(j)
N.remove(j)
break
if cur_k is None:
continue
while cur_node != depot:
for j in N:
if X[cur_node, j, cur_k].x > 0:
cur_node = j
route.append(j)
if j != depot:
N.remove(j)
break
route.insert(0,k)
route_list.append(route)
# 提取车辆的时刻点
for route in route_list:
k = route[0]
timetable = []
for i in route[1:]:
timetable.append(T[i,k].x)
route_timetable.append(timetable)
print("最优路径距离:", vrptw_model.objVal)
print("最优路径使用车辆数:", len(route_list))
draw_routes(route_list, x_coord, y_coord)
save_file(route_list, route_timetable, vrptw_model.objVal)
5.3 Python调用COPT求解CVRP
import xlsxwriter
import math
import pandas as pd
import matplotlib.pyplot as plt
from coptpy import *
def read_input(filename):
"""
:param filename: 数据文件路径
:return:
"""
df = pd.read_csv(filename)
x_coord = { df['id'][i]:df['x_coord'][i] for i in range(df.shape[0]) } # 节点横坐标
y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) } # 节点纵坐标
demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) } # 节点需求
cost = {}
N = df['id'].tolist()
for f_n in N:
for t_n in N:
dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
cost[f_n,t_n] = dist
return N,cost,demand,x_coord,y_coord
def build_model(N,K,depot,CAP,cost,demand):
# 创建求解环境
env = Envr()
# 创建求解模型
cvrp_model = env.createModel('cvrp')
# 添加决策变量
X = cvrp_model.addVars(N, N, K, vtype=COPT.BINARY, nameprefix='X[i,j,k]')
Y = cvrp_model.addVars(N, K, vtype=COPT.BINARY, nameprefix='Y[i,k]')
U = cvrp_model.addVars(N, K, vtype=COPT.INTEGER, nameprefix='U[i,k]')
# 设置目标函数
cvrp_model.setObjective(quicksum(X[i, j, k] * cost[i, j] for i in N for j in N for k in K), sense=COPT.MINIMIZE)
# 添加约束
# 需求覆盖约束
cvrp_model.addConstrs(quicksum(Y[i, k] for k in K) == 1 for i in N[1:])
# 车辆启用约束
cvrp_model.addConstr(quicksum(Y[depot, k] for k in K) == len(K))
# 车辆流平衡约束
cvrp_model.addConstrs(quicksum(X[i, j, k] for j in N) == quicksum(X[j, i, k] for j in N) for i in N for k in K)
# 车辆路径限制
cvrp_model.addConstrs(quicksum(X[i, j, k] for j in N) == Y[i, k] for i in N for k in K)
# 车辆容量约束
cvrp_model.addConstrs(quicksum(Y[i, k] * demand[i] for i in N) <= CAP for k in K)
# 破圈约束
cvrp_model.addConstrs(
U[i, k] - U[j, k] + CAP * X[i, j, k] <= CAP - demand[i] for i in N[1:] for j in N[1:] for k in K)
cvrp_model.addConstrs(U[i, k] <= CAP for i in N[1:] for k in K)
cvrp_model.addConstrs(U[i, k] >= demand[i] for i in N[1:] for k in K)
return cvrp_model, X, Y, U
def draw_routes(route_list,x_coord,y_coord):
for route in route_list:
path_x = []
path_y = []
for n in route:
path_x.append(x_coord[n])
path_y.append(y_coord[n])
plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
plt.show()
def save_file(route_list,total_cost):
wb = xlsxwriter.Workbook('路径方案.xlsx')
ws = wb.add_worksheet()
ws.write(0,0,'总费用')
ws.write(0,1,total_cost)
ws.write(1,0,'车辆')
ws.write(1,1,'路径')
row = 2
for route in route_list:
ws.write(row,0,route[0])
route_str = [str(i) for i in route[1:]]
ws.write(row,1,'-'.join(route_str))
row += 1
wb.close()
if __name__ == '__main__':
filename = './datasets/CVRP/r101-31.csv'
N, cost, demand, x_coord,y_coord = read_input(filename)
depot = N[0]
K = list(range(1,10))
CAP = 80
cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)
cvrp_model.param.timelimit = 1200
cvrp_model.setLogFile('./copt_r101-31.log')
cvrp_model.solve()
route_list = []
for k in K:
route = [depot]
cur_node = depot
cur_k = None
for j in N[1:]:
if X[depot, j, k].x > 0:
cur_node = j
cur_k = k
route.append(j)
N.remove(j)
break
if cur_k is None:
continue
while cur_node != depot:
for j in N:
if X[cur_node, j, cur_k].x > 0:
cur_node = j
route.append(j)
if j != depot:
N.remove(j)
break
route_list.append(route)
print("最优路径距离:", cvrp_model.objval)
print("最优路径使用车辆数:", len(route_list))
draw_routes(route_list, x_coord, y_coord)
save_file(route_list, cvrp_model.objval)
5.4 Python调用COPT求解VRPTW
import math
import pandas as pd
import matplotlib.pyplot as plt
import xlsxwriter
from coptpy import Envr,COPT,quicksum
def read_input(filename):
"""
:param filename: 数据文件路径
:return:
"""
N = [] #所有节点
Q = {} #节点需求
TT = {} #节点旅行时间
ET = {} #节点最早开始服务时间
LT = {} #节点最晚结束服务时间
ST = {} #节点服务时间
x_coord = {} # 节点横坐标
y_coord = {} # 节点纵坐标
Cost={}
df = pd.read_csv(filename)
for i in range(df.shape[0]):
id = df['id'][i]
N.append(id)
x_coord[id] = df['x_coord'][i]
y_coord[id] = df['y_coord'][i]
Q[id] = df['demand'][i]
ET[id] = df['start_time'][i]
LT[id] = df['end_time'][i]
ST[id] = df['service_time'][i]
for f_n in N:
for t_n in N:
dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
Cost[f_n,t_n] = dist
TT[f_n,t_n] = dist/1 # 假设速度为1
return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coord
def build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K):
"""
:param N: 所有节点集合，其中N[0]为车场
:param Q: 节点需求集合
:param TT: 旅行时间
:param ET: 节点最早开始服务时间
:param LT：节点最晚结束服务时间
:param ST: 节点服务时间
:param CAP: 车辆容量
:param Cost: 旅行费用
:param K: 车队
:return:
"""
C = N[1:] #需求节点
M=10**5
depot = N[0]
# 创建求解环境
env = Envr()
# 创建求解模型
vrptw_model = env.createModel('vrptw')
# 添加决策变量
X = vrptw_model.addVars(N, N, K, vtype=COPT.BINARY, nameprefix='X[i,j,k]')
T = vrptw_model.addVars( N,K,vtype=COPT.CONTINUOUS,lb=0,nameprefix='T[i,k]')
# 设置目标函数
vrptw_model.setObjective(quicksum(X[i, j, k] * Cost[i, j] for i in N for j in N for k in K), sense=COPT.MINIMIZE)
# 车辆起点约束
vrptw_model.addConstrs(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
# 车辆路径连续约束
vrptw_model.addConstrs(
quicksum(X[i, j, k] for j in N if j != i) == quicksum(X[j, i, k] for j in N if j != i) for i in C for k in K)
# 车辆终点约束
vrptw_model.addConstrs(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
# 需求服务约束
vrptw_model.addConstrs(quicksum(X[i, j, k] for k in K for j in N if j != i) == 1 for i in C)
# 车辆容量约束
vrptw_model.addConstrs(quicksum(Q[i] * X[i, j, k] for i in C for j in N if i != j) <= CAP for k in K)
# 时间窗约束
vrptw_model.addConstrs(
T[i, k] + ST[i] + TT[i, j] - (1 - X[i, j, k]) * M <= T[j, k] for i in C for j in C for k in K if i != j)
vrptw_model.addConstrs(T[i, k] >= ET[i] for i in N for k in K)
vrptw_model.addConstrs(T[i, k] <= LT[i] for i in N for k in K)
return vrptw_model, X, T
def draw_routes(route_list,x_coord,y_coord):
for route in route_list:
path_x = []
path_y = []
for n in route:
path_x.append(x_coord[n])
path_y.append(y_coord[n])
plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
plt.show()
def save_file(route_list,route_timetable,total_cost):
wb = xlsxwriter.Workbook('路径方案.xlsx')
ws = wb.add_worksheet()
ws.write(0,0,'总费用')
ws.write(0,1,total_cost)
ws.write(1,0,'车辆')
ws.write(1,1,'路径')
ws.write(1,2,'时刻')
row = 2
for id,route in enumerate(route_list):
ws.write(row,0,route[0])
route_str = [str(i) for i in route[1:]]
ws.write(row,1,'-'.join(route_str))
timetable_str = [str(i) for i in route_timetable[id]]
ws.write(row,2,'-'.join(timetable_str))
row += 1
wb.close()
if __name__=='__main__':
filename = './data/r101.csv'
N, Q, TT, ET, LT, ST, Cost,x_coord,y_coord = read_input(filename)
depot = N[0]
K = list(range(1, 30))
CAP = 80
vrptw_model,X,T = build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K)
vrptw_model.setLogFile('./log/copt_c201.log')
vrptw_model.param.timelimit = 1500
vrptw_model.solve()
route_list = []
route_timetable = []
# 提取车辆路径
for k in K:
route = [depot]
cur_node = depot
cur_k = None
for j in N[1:]:
if X[depot, j, k].x > 0:
cur_node = j
cur_k = k
route.append(j)
N.remove(j)
break
if cur_k is None:
continue
while cur_node != depot:
for j in N:
if X[cur_node, j, cur_k].x > 0:
cur_node = j
route.append(j)
if j != depot:
N.remove(j)
break
route.insert(0,k)
route_list.append(route)
# 提取车辆的时刻点
for route in route_list:
k = route[0]
timetable = []
for i in route[1:]:
timetable.append(T[i,k].x)
route_timetable.append(timetable)
print("最优路径距离:", vrptw_model.objval)
print("最优路径使用车辆数:", len(route_list))
draw_routes(route_list, x_coord, y_coord)
save_file(route_list, route_timetable, vrptw_model.objval)
5.5 Python调用SCIP求解CVRP
import xlsxwriter
import math
import pandas as pd
import matplotlib.pyplot as plt
from pyscipopt import Model, quicksum, multidict
def read_input(filename):
"""
:param filename: 数据文件路径
:return:
"""
df = pd.read_csv(filename)
x_coord = { df['id'][i]:df['x_coord'][i] for i in range(df.shape[0]) } # 节点横坐标
y_coord = { df['id'][i]:df['y_coord'][i] for i in range(df.shape[0]) } # 节点纵坐标
demand = { df['id'][i]:df['demand'][i] for i in range(df.shape[0]) } # 节点需求
cost = {}
N = df['id'].tolist()
for f_n in N:
for t_n in N:
dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
cost[f_n,t_n] = dist
return N,cost,demand,x_coord,y_coord
def build_model(N,K,depot,CAP,cost,demand):
"""
:param N: 网络节点集合
:param K: 车队集合
:param depot: 配送中心id
:param CAP: 车辆容量
:param cost: 网络弧费用集合
:param demand: 网络节点需求集合
:return:
"""
cvrp_model = Model('cvrp')
# 添加变量
X = {}
Y = {}
U = {}
# 创建变量 X[i,j,k]
for i in N:
for j in N:
for k in K:
X[i,j,k] = cvrp_model.addVar(vtype="B", name=f"x({i},{j},{k})")
# 创建变量 Y[i,k]
for i in N:
for k in K:
Y[i,k] = cvrp_model.addVar(vtype="B", name=f"y({i},{k})")
# 创建辅助变量U[i,k]
for i in N:
for k in K:
U[i,k] = cvrp_model.addVar(vtype="I", name=f"u({i},{k})")
# 设置目标函数
cvrp_model.setObjective( quicksum(X[i,j,k] * cost[i,j] for i in N for j in N for k in K),'minimize' )
# 添加约束
# 需求覆盖约束
cvrp_model.addConss( quicksum(Y[i,k] for k in K) == 1 for i in N[1:] )
# 车辆启用约束
cvrp_model.addCons( quicksum(Y[depot,k] for k in K) == len(K) )
# 车辆流平衡约束
cvrp_model.addConss( quicksum(X[i,j,k] for j in N ) == quicksum(X[j,i,k] for j in N ) for i in N for k in K )
# 车辆路径限制
cvrp_model.addConss( quicksum(X[i,j,k] for j in N) == Y[i,k] for i in N for k in K )
# 车辆容量约束
cvrp_model.addConss( quicksum(Y[i,k] * demand[i] for i in N) <= CAP for k in K )
# 破圈约束
cvrp_model.addConss( U[i,k] - U[j,k] + CAP * X[i,j,k] <= CAP - demand[i] for i in N[1:] for j in N[1:] for k in K )
cvrp_model.addConss( U[i,k] <= CAP for i in N[1:] for k in K)
cvrp_model.addConss( U[i,k] >= demand[i] for i in N[1:] for k in K )
return cvrp_model, X, Y, U
def draw_routes(route_list,x_coord,y_coord):
for route in route_list:
path_x = []
path_y = []
for n in route:
path_x.append(x_coord[n])
path_y.append(y_coord[n])
plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
plt.show()
def save_file(route_list,total_cost):
wb = xlsxwriter.Workbook('路径方案.xlsx')
ws = wb.add_worksheet()
ws.write(0,0,'总费用')
ws.write(0,1,total_cost)
ws.write(1,0,'车辆')
ws.write(1,1,'路径')
row = 2
for route in route_list:
ws.write(row,0,route[0])
route_str = [str(i) for i in route[1:]]
ws.write(row,1,'-'.join(route_str))
row += 1
wb.close()
if __name__ == '__main__':
filename = './datasets/CVRP/r101-31.csv'
N, cost, demand, x_coord, y_coord = read_input(filename)
depot = N[0]
K = list(range(1,10))
CAP = 80
cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)
# cvrp_model.setRealParam('limits/gap', 0.1) # 求解gap为10%
cvrp_model.setRealParam('limits/time', 1200) # 求解时间为1800s
cvrp_model.setLogfile('./scip_r101-31.log')
cvrp_model.optimize()
route_list = []
for k in K:
route = [depot]
cur_node = depot
cur_k = None
for j in N[1:]:
if cvrp_model.getVal(X[depot, j, k]) > 0:
cur_node = j
cur_k = k
route.append(j)
N.remove(j)
break
if cur_k is None:
continue
while cur_node != depot:
for j in N:
if cvrp_model.getVal(X[cur_node, j, cur_k]) > 0:
cur_node = j
route.append(j)
if j != depot:
N.remove(j)
break
route_list.append(route)
print("最优路径距离:", cvrp_model.getObjVal())
print("最优路径使用车辆数:", len(route_list))
draw_routes(route_list, x_coord, y_coord)
save_file(route_list, cvrp_model.getObjVal())
5.6 Python调用SCIP求解VRPTW
import math
import pandas as pd
import matplotlib.pyplot as plt
import xlsxwriter
from pyscipopt import Model, quicksum
def read_input(filename):
"""
:param filename: 数据文件路径
:return:
"""
N = [] #所有节点
Q = {} #节点需求
TT = {} #节点旅行时间
ET = {} #节点最早开始服务时间
LT = {} #节点最晚结束服务时间
ST = {} #节点服务时间
x_coord = {} # 节点横坐标
y_coord = {} # 节点纵坐标
Cost={}
df = pd.read_csv(filename)
for i in range(df.shape[0]):
id = df['id'][i]
N.append(id)
x_coord[id] = df['x_coord'][i]
y_coord[id] = df['y_coord'][i]
Q[id] = df['demand'][i]
ET[id] = df['start_time'][i]
LT[id] = df['end_time'][i]
ST[id] = df['service_time'][i]
for f_n in N:
for t_n in N:
dist = math.sqrt( (x_coord[f_n]-x_coord[t_n])**2 + (y_coord[f_n] - y_coord[t_n])**2 )
Cost[f_n,t_n] = dist
TT[f_n,t_n] = dist/1 # 假设速度为1
return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coord
def build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K):
"""
:param N: 所有节点集合，其中N[0]为车场
:param Q: 节点需求集合
:param TT: 旅行时间
:param ET: 节点最早开始服务时间
:param LT：节点最晚结束服务时间
:param ST: 节点服务时间
:param CAP: 车辆容量
:param Cost: 旅行费用
:param K: 车队
:return:
"""
C = N[1:] #需求节点
M=10**5
depot = N[0]
# 创建模型
vrptw_model=Model()
X = {}
T = {}
# 创建变量 X[i,j,k]
for i in N:
for j in N:
for k in K:
X[i, j, k] = vrptw_model.addVar(vtype="B", name=f"x({i},{j},{k})")
# 创建变量 T[i,k]
for i in N:
for k in K:
T[i, k] = vrptw_model.addVar(vtype="I", name=f"t({i},{k})")
# 设置目标函数
vrptw_model.setObjective(quicksum(X[i, j, k] * Cost[i, j] for i in N for j in N for k in K), 'minimize')
# 车辆起点约束
# 车辆起点约束
vrptw_model.addConss(quicksum(X[depot, j, k] for j in N) == 1 for k in K)
# 车辆路径连续约束
vrptw_model.addConss(
quicksum(X[i, j, k] for j in N if j != i) == quicksum(X[j, i, k] for j in N if j != i) for i in C for k in K)
# 车辆终点约束
vrptw_model.addConss(quicksum(X[j, depot, k] for j in N) == 1 for k in K)
# 需求服务约束
vrptw_model.addConss(quicksum(X[i, j, k] for k in K for j in N if j != i) == 1 for i in C)
# 车辆容量约束
vrptw_model.addConss(quicksum(Q[i] * X[i, j, k] for i in C for j in N if i != j) <= CAP for k in K)
# 时间窗约束
vrptw_model.addConss(
T[i, k] + ST[i] + TT[i, j] - (1 - X[i, j, k]) * M <= T[j, k] for i in C for j in C for k in K if i != j)
vrptw_model.addConss(T[i, k] >= ET[i] for i in N for k in K)
vrptw_model.addConss(T[i, k] <= LT[i] for i in N for k in K)
return vrptw_model, X, T
def draw_routes(route_list,x_coord,y_coord):
for route in route_list:
path_x = []
path_y = []
for n in route:
path_x.append(x_coord[n])
path_y.append(y_coord[n])
plt.plot(path_x, path_y,linewidth=0.5, marker='s',ms=5)
plt.show()
def save_file(route_list,route_timetable,total_cost):
wb = xlsxwriter.Workbook('路径方案.xlsx')
ws = wb.add_worksheet()
ws.write(0,0,'总费用')
ws.write(0,1,total_cost)
ws.write(1,0,'车辆')
ws.write(1,1,'路径')
ws.write(1,2,'时刻')
row = 2
for id,route in enumerate(route_list):
ws.write(row,0,route[0])
route_str = [str(i) for i in route[1:]]
ws.write(row,1,'-'.join(route_str))
timetable_str = [str(i) for i in route_timetable[id]]
ws.write(row,2,'-'.join(timetable_str))
row += 1
wb.close()
if __name__=='__main__':
filename = './data/c201.csv'
N, Q, TT, ET, LT, ST, Cost,x_coord,y_coord = read_input(filename)
depot = N[0]
K = list(range(1, 30))
CAP = 80
vrptw_model,X,T = build_model(N,Q,TT,ET,LT,ST,Cost,CAP,K)
vrptw_model.setLogfile('./log/scip_c101.log')
vrptw_model.setRealParam('limits/time', 2000) # 求解时间为1800s
vrptw_model.optimize()
route_list = []
route_timetable = []
# 提取车辆路径
for k in K:
route = [depot]
cur_node = depot
cur_k = None
for j in N[1:]:
if vrptw_model.getVal(X[depot, j, k]) > 0:
cur_node = j
cur_k = k
route.append(j)
N.remove(j)
break
if cur_k is None:
continue
while cur_node != depot:
for j in N:
if vrptw_model.getVal(X[cur_node, j, cur_k]) > 0:
cur_node = j
route.append(j)
if j != depot:
N.remove(j)
break
route.insert(0,k)
route_list.append(route)
# 提取车辆的时刻点
for route in route_list:
k = route[0]
timetable = []
for i in route[1:]:
timetable.append(vrptw_model.getVal(T[i,k]))
route_timetable.append(timetable)
print("最优路径距离:", vrptw_model.getObjVal())
print("最优路径使用车辆数:", len(route_list))
draw_routes(route_list, x_coord, y_coord)
save_file(route_list, route_timetable, vrptw_model.getObjVal())
7. 测试案例
本期从solomon数据集中选择了c101，c201，cr101三个算例用于测评。
CVRP问题数据结构如下（保留30个节点）：
VRPTW问题数据结构如下（保留100个节点）：
8. 测试参数
软件版本号：
Gurbi：9.5.1
SCIP：8.0.0
COPT：5.0.1
CVRP相关参数设置：
车辆容量：80
车队规模：10
客户节点：30
求解时间：1200s
VRPTW相关参数设置：
车辆容量：80
车队规模：30
客户节点：100
求解时间：1500s
9. 测试结果
9.1 CVRP求解结果汇总
对于CVRP问题，在相同的求解时间下，Gurobi在求解质量方面有最好的表现，COPT比SCIP略差。
9.2 VRPTW求解结果汇总
对于VRPTW问题，在相同的求解时间下，Gurobi可以找到不错的优化解，而COPT和SCIP并未找到初始解。
9.3 上界下降曲线对比（以CVRP为例）
从上界下降曲线可以看出，Gurobi找到的初始可行解的质量最好，且收敛性最好，而COPT和SCIP的初始解质量相差不大，收敛性不如Gurobi。
9.5 车辆路径可视化（以CVRP为例）
c101-31（Gurobi）
c201-31（Gurobi）
r101-31（Gurobi）
10. 小节
根据以上测试结果，Gurobi的求解性能远优于SCIP和COPT，作为后起之秀COPT已经接近SCIP水平（最新版的COPT尚未测试，也许性能有更大提升）。
虽然求解器可以直接对数学模型进行求解，但仅在小规模案例上具有较好的性能，对于大规模问题仍需要设计更高效的算法。
小编后续将继续分享VRP问题系列求解算法，并将其求解效果与上述求解器进行比较。
