目录
- 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问题系列求解算法,并将其求解效果与上述求解器进行比较。