软件工具 | Python调用运筹优化求解器(一):以CVRPVRPTW为例

news/2025/1/11 7:37:10/

目录

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


http://www.ppmy.cn/news/64522.html

相关文章

软件测试面试宝典,最常见的7个高频面试题(附答案,建议收藏)

收集了2022年所有黑马学员的面试题后&#xff0c;负责就业的黑马讲师们整理出了7个高频出现的面试题&#xff0c;一起来看看。 高频问题1&#xff1a;请自我介绍下&#xff1f; 高频问题2&#xff1a;请介绍下最近做过的项目&#xff1f; 高频问题3&#xff1a;请介绍下你印象…

【Ubuntu18.04】Docker配置镜像源

作者主页&#xff1a;爱笑的男孩。的博客_CSDN博客-深度学习,活动,YOLO领域博主爱笑的男孩。擅长深度学习,活动,YOLO,等方面的知识,爱笑的男孩。关注算法,python,计算机视觉,图像处理,深度学习,pytorch,神经网络,opencv领域.https://blog.csdn.net/Code_and516?typeblog个人简…

ADS-B教学实验方案

ADS-B教学系统是为了让学生学习ADS-B原理、ADS-B系统组成、ADS-B信号处理技术。可以通过ADS-B教学系统进一步研究分析ADS-B位置的精度、准确性、稳定性、实时性&#xff0c;设计基于ADS-B的空中碰撞告警系统&#xff0c;混合空域的空中交通管理系统(UTM)设计。也可以研究ADS-B报…

iOS与Android应用开发的对比:如何选择最佳开发平台?

第一章&#xff1a;引言 在移动应用开发领域&#xff0c;iOS和Android是最为流行的操作系统。选择最佳的开发平台可以使开发人员更有效地开发和发布应用程序。本文将分析iOS和Android应用开发的优缺点&#xff0c;并提供一些有关如何选择最佳开发平台的建议。 第二章&#xf…

【工具】使用VS Code调试Docker Container中的代码

目录 使用VS Code调试Docker Container中的Autoware.ai代码Part 1 -- 在VS Code中设置并进行DebugStep 1Step 2Step 3Step 4launch.jsontasks.jsonc_cpp_properties.jsonsettings.json Step 5Step 6Step 7Error Solutions参考链接 Part 2 -- cmake重新编译cmake使用方法&#x…

icevision环境安装

Installation - IceVision # 1. git clone 代码# pip 换源&#xff1a; ~/.pip/pip.conf 隐藏文件[global] index-url https://pypi.tuna.tsinghua.edu.cn/simple [install] trusted-hostmirrors.aliyun.compip install -e .[all,dev]ImportError: cannot import name Multi…

MySQL数据库——MySQL SELECT:数据表查询语句

在 MySQL 中&#xff0c;可以使用 SELECT 语句来查询数据。查询数据是指从数据库中根据需求&#xff0c;使用不同的查询方式来获取不同的数据&#xff0c;是使用频率最高、最重要的操作。 SELECT 的语法格式如下&#xff1a; SELECT {* | <字段列名>} [ FROM <表 1&…

Git 从入门到精通

在软件开发和版本控制领域&#xff0c;Git已经成为了最流行和广泛使用的工具之一。它不仅仅是一个版本控制系统&#xff0c;还是一个强大而灵活的工具&#xff0c;可以帮助开发者更好地管理代码、协作团队以及追踪项目的历史。本文将从Git的基础概念开始&#xff0c;逐步介绍Gi…