当前位置：首页 > news >正文

软件工具 | Python调用运筹优化求解器（一）：以CVRPVRPTW为例

news 文章来源：https://blog.csdn.net/python_n/article/details/130551994 2024/11/16 10:31:25

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] = distreturn N,cost,demand,x_coord,y_coorddef 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,Udef 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 = 2for 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 += 1wb.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 = 80cvrp_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 = 1200cvrp_model.optimize()route_list = []for k in K:route = [depot]cur_node = depotcur_k = Nonefor j in N[1:]:if X[depot, j, k].x > 0:cur_node = jcur_k = kroute.append(j)N.remove(j)breakif cur_k is None:continuewhile cur_node != depot:for j in N:if X[cur_node, j, cur_k].x > 0:cur_node = jroute.append(j)if j != depot:N.remove(j)breakroute_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,tuplelistdef 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] = distTT[f_n,t_n] = dist/1 # 假设速度为1return N,Q,TT,ET,LT,ST,Cost,x_coord,y_coorddef 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**5depot = 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 = 2for 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 += 1wb.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 = 80vrptw_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 = 1500vrptw_model.optimize()route_list = []route_timetable = []# 提取车辆路径for k in K:route = [depot]cur_node = depotcur_k = Nonefor j in N[1:]:if X[depot, j, k].x > 0:cur_node = jcur_k = kroute.append(j)N.remove(j)breakif cur_k is None:continuewhile cur_node != depot:for j in N:if X[cur_node, j, cur_k].x > 0:cur_node = jroute.append(j)if j != depot:N.remove(j)breakroute.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] = distreturn 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 = 2for 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 += 1wb.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 = 80cvrp_model, X, Y, U = build_model(N, K, depot, CAP, cost, demand)cvrp_model.param.timelimit = 1200cvrp_model.setLogFile('./copt_r101-31.log')cvrp_model.solve()route_list = []for k in K:route = [depot]cur_node = depotcur_k = Nonefor j in N[1:]:if X[depot, j, k].x > 0:cur_node = jcur_k = kroute.append(j)N.remove(j)breakif cur_k is None:continuewhile cur_node != depot:for j in N:if X[cur_node, j, cur_k].x > 0:cur_node = jroute.append(j)if j != depot:N.remove(j)breakroute_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,quicksumdef 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] = distTT[f_n,t_n] = dist/1 # 假设速度为1return 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**5depot = 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 = 2for 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 += 1wb.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 = 80vrptw_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 = 1500vrptw_model.solve()route_list = []route_timetable = []# 提取车辆路径for k in K:route = [depot]cur_node = depotcur_k = Nonefor j in N[1:]:if X[depot, j, k].x > 0:cur_node = jcur_k = kroute.append(j)N.remove(j)breakif cur_k is None:continuewhile cur_node != depot:for j in N:if X[cur_node, j, cur_k].x > 0:cur_node = jroute.append(j)if j != depot:N.remove(j)breakroute.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, multidictdef 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] = distreturn N,cost,demand,x_coord,y_coorddef 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, Udef 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 = 2for 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 += 1wb.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 = 80cvrp_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)  # 求解时间为1800scvrp_model.setLogfile('./scip_r101-31.log')cvrp_model.optimize()route_list = []for k in K:route = [depot]cur_node = depotcur_k = Nonefor j in N[1:]:if cvrp_model.getVal(X[depot, j, k]) > 0:cur_node = jcur_k = kroute.append(j)N.remove(j)breakif cur_k is None:continuewhile cur_node != depot:for j in N:if cvrp_model.getVal(X[cur_node, j, cur_k]) > 0:cur_node = jroute.append(j)if j != depot:N.remove(j)breakroute_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, quicksumdef 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] = distTT[f_n,t_n] = dist/1 # 假设速度为1return 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**5depot = 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 = 2for 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 += 1wb.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 = 80vrptw_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)  # 求解时间为1800svrptw_model.optimize()route_list = []route_timetable = []# 提取车辆路径for k in K:route = [depot]cur_node = depotcur_k = Nonefor j in N[1:]:if vrptw_model.getVal(X[depot, j, k]) > 0:cur_node = jcur_k = kroute.append(j)N.remove(j)breakif cur_k is None:continuewhile cur_node != depot:for j in N:if vrptw_model.getVal(X[cur_node, j, cur_k]) > 0:cur_node = jroute.append(j)if j != depot:N.remove(j)breakroute.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问题系列求解算法，并将其求解效果与上述求解器进行比较。

目录

1. 引言

2. 求解器介绍

3. 基础语言

3.1 创建模型

3.2 添加变量

3.3 添加目标函数

3.4 添加约束

3.5 设置参数

3.6 求解

4. 数学模型

4.1 CVRP数学模型

4.2 VRPTW数学模型

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. 小节

相关文章：