运筹优化带时间窗约束的车辆路径规划问题（VRPTW）详解+Python调用Gurobi建模求解

见贤思齐

文章目录一、概述1.1VRP问题1.2CVRP问题1.3VRPTW问题二、VRPTW的一般模型三、Python调用Gurobi建模求解3.1Solomn数据集3.2完整代码3.3运行结果展示3.3.1测试案例：c101.txt3.3.2测试案例：r101.txt一、概述1.1VRP问题车辆路径规划问题(VehicleRoutingProblem，VRP)一般指的是：对一系列发货点和收货点，组织调用一定的车辆，安排适当的行车路线，使车辆有序地通过它们，在满足指定的约束条件下（例如：货物的需求量与发货量，交发货时间，车辆容量限制，行驶里程限制，行驶时间限制等），力争实现一定的目标（如车辆空驶总里程最短，运输总费用最低，车辆按一定时间到达，使用的车辆数最小等）。下图给出了一个简单的VRP的例子1.2CVRP问题最基本的VRP问题叫做带容量约束的车辆路径规划问题（CapacitatedVehicleRoutingProblem，CVRP）。在CVRP中，需要考虑每辆车的容量约束、车辆的路径约束和装载量约束1.3VRPTW问题为了考虑配送时间要求，带时间窗的车辆路径规划问题（VehicleRoutingProblemwithTimeWindow，VRPTW）应运而生。VRPTW不仅考虑CVRP的所有约束，还需要考虑时间窗约束，也就是每个顾客对应一个时间窗[ei,li][e_i,l_i][ei​,li​]，其中eie_iei​和lil_ili​分别代表该点的最早到达时间和最晚到达时间。顾客点i∈Vi\inVi∈V的需求必须要在其时间窗内被送达VRPTW已经被证明是NP-hard问题，其求解复杂度随着问题规模的增加而急剧增加，求解较为困难。到目前为止，求解VRPTW的比较高效的精确算法是分支定价算法和分支定价切割算法。二、VRPTW的一般模型VRPTW可以建模为一个混合整数规划问题，在给出完整数学模型之前，先引入下面的决策变量：xijk={1，在最优解中，弧(i,j)被车辆k选中0，其他sik=车辆k到达i的时间模型中涉及的其他参数为:tij表示车辆在弧(i,j)上的行驶时间M为一个足够大的正数{x_{ijk}}={1，在最优解中，弧(i,j)被车辆k选中0，其他{1，在最优解中，弧(i,j)被车辆k选中0，其他\\{s_{ik}}=\text{车辆}k\text{到达}i\text{的时间}\\\text{模型中涉及的其他参数为}:\\{t_{ij}}\text{表示车辆在弧}\left(i,j\right)\text{上的行驶时间}\\M\text{为一个足够大的正数}xijk​={1，在最优解中，弧(i,j)被车辆k选中0，其他​sik​=车辆k到达i的时间模型中涉及的其他参数为:tij​表示车辆在弧(i,j)上的行驶时间M为一个足够大的正数关于MMM的取值，实际上可以直接取非常大的正数，但是为了提高求解效率，拉紧约束。我们可以采用下面的取值方法：M=max{bi+tij−aj},∀(i,j)∈AM=max\{b_i+t_{ij}-a_j\},\forall(i,j)\inAM=max{bi​+tij​−aj​},∀(i,j)∈A综合上面引出的决策变量，并参考文献（Desaulniersetal.，2006），给出的VRPTW的标准模型如下：min⁡∑k∈K∑i∈V∑j∈Vcijxijks.t.∑k∈K∑j∈Vxijk=1,∀i∈C  ∑j∈Vx0jk=1,∀k∈K  ∑i∈Vxihk−∑j∈Vxhjk=0,∀h∈C,∀k∈K  ∑i∈Vxi,n+1,k=1,∀k∈K  ∑i∈Cqi∑j∈Vxijk⩽Q,∀k∈K  sik+tij−M(1−xijk)⩽sjk  ,∀(i,j)∈A,∀k∈K  ei⩽sik⩽li  ,∀i∈V,∀k∈K  xijk∈{0,1}  ,∀(i,j)∈A,∀k∈K\min\sum_{k\inK}{\sum_{i\inV}{\sum_{j\inV}{{c_{ij}}{x_{ijk}}}}}\\s.t.\sum_{k\inK}{\sum_{j\inV}{{x_{ijk}}=1,\foralli\inC}}\\\,\,\sum_{j\inV}{{x_{0jk}}=1,\forallk\inK}\\\,\,\sum_{i\inV}{{x_{ihk}}-\sum_{j\inV}{{x_{hjk}}=0,\forallh\inC,\forallk\inK}}\\\,\,\sum_{i\inV}{x_{i,n+1,k}=1,\forallk\inK}\\\,\,\sum_{i\inC}{q_i\sum_{j\inV}{{x_{ijk}}\leqslantQ,\forallk\inK}}\\\,\,{s_{ik}}+{t_{ij}}-M\left(1-{x_{ijk}}\right)\leqslant{s_{jk}}\,\,,\forall\left(i,j\right)\inA,\forallk\inK\\\,\,e_i\leqslant{s_{ik}}\leqslantl_i\,\,,\foralli\inV,\forallk\inK\\\,\,{x_{ijk}}\in\left\{0,1\right\}\,\,,\forall\left(i,j\right)\inA,\forallk\inKmink∈K∑​i∈V∑​j∈V∑​cij​xijk​s.t.k∈K∑​j∈V∑​xijk​=1,∀i∈Cj∈V∑​x0jk​=1,∀k∈Ki∈V∑​xihk​−j∈V∑​xhjk​=0,∀h∈C,∀k∈Ki∈V∑​xi,n+1,k​=1,∀k∈Ki∈C∑​qi​j∈V∑​xijk​⩽Q,∀k∈Ksik​+tij​−M(1−xijk​)⩽sjk​,∀(i,j)∈A,∀k∈Kei​⩽sik​⩽li​,∀i∈V,∀k∈Kxijk​∈{0,1},∀(i,j)∈A,∀k∈K其中，QQQ为车容量，qiq_iqi​为第iii个顾客的需求：目标函数是为了最小化所有车辆的总行驶成本（距离）约束1~4保证了每辆车必须从仓库出发，经过一个点就离开那个点，最终返回仓库约束5为车辆的容量约束约束6~7是时间窗约束，保证了车辆到达每个顾客点的时间均在时间窗内，点n+1是点o的一个备份，是为了方便实现。三、Python调用Gurobi建模求解3.1Solomn数据集Solomn数据集下载地址3.2完整代码注意，在下面代码中，将弧iii到弧jjj所需的时间tijt_{ij}tij​和成本cijc_{ij}cij​都当作了弧iii到弧jjj所需的距离来看待#-*-coding:utf-8-*-##Author:WSKH#Blog:wskh0929.blog.csdn.net#Time:2023/2/811:14#Descriptionython调用Gurobi建模求解VRPTW问题importtimeimportmatplotlib.pyplotaspltimportnumpyasnpfromgurobipyimport*classData:customerNum=0nodeNum=0vehicleNum=0capacity=0corX=[]corY=[]demand=[]serviceTime=[]readyTime=[]dueTime=[]distanceMatrix=[[]]defreadData(path,customerNum):data=Data()data.customerNum=customerNumifcustomerNumisnotNone:data.nodeNum=customerNum+2withopen(path,'r')asf:lines=f.readlines()count=0forlineinlines:count+=1ifcount==5:line=line[:-1]s=re.split(r"+",line)data.vehicleNum=int(s[1])data.capacity=float(s[2])elifcount>=10and(customerNumisNoneorcount0.5:solution.X[int(split_arr[1])][int(split_arr[2])][int(split_arr[3])]=m.xelifsplit_arr[0]=='S'andm.x>0.5:solution.S[int(split_arr[1])][int(split_arr[2])]=m.xforkinrange(data.vehicleNum):i=0subRoute=[]subRoute.append(i)finish=Falsewhilenotfinish:forjinrange(data.nodeNum):ifsolution.X[i][j][k]>0.5:subRoute.append(j)i=jifj==data.nodeNum-1:finish=Trueiflen(subRoute)>=3:subRoute[-1]=0solution.routes.append(subRoute)solution.routeNum+=1returnsolutiondefplot_solution(solution,customer_num):plt.xlabel("x")plt.ylabel("y")plt.title(f"{data_type}:{customer_num}Customers")plt.scatter(data.corX[0],data.corY[0],c='blue',alpha=1,marker=',',linewidths=3,label='depot')#起点plt.scatter(data.corX[1:-1],data.corY[1:-1],c='black',alpha=1,marker='o',linewidths=3,label='customer')#普通站点forkinrange(solution.routeNum):foriinrange(len(solution.routes[k])-1):a=solution.routes[k][i]b=solution.routes[k][i+1]x=[data.corX[a],data.corX[b]]y=[data.corY[a],data.corY[b]]plt.plot(x,y,'k',linewidth=1)plt.grid(False)plt.legend(loc='best')plt.show()defprint_solution(solution,data):forindex,subRouteinenumerate(solution.routes):distance=0load=0foriinrange(len(subRoute)-1):distance+=data.distanceMatrix[subRoute[i]][subRoute[i+1]]load+=data.demand[subRoute[i]]print(f"Route-{index+1}:{subRoute},distance:{distance},load:{load}")defsolve(data):#声明模型model=Model("VRPTW")#模型设置#关闭输出model.setParam('OutputFlag',0)#定义变量X=[[[[]for_inrange(data.vehicleNum)]for_inrange(data.nodeNum)]for_inrange(data.nodeNum)]S=[[[]for_inrange(data.vehicleNum)]for_inrange(data.nodeNum)]foriinrange(data.nodeNum):forkinrange(data.vehicleNum):S[i][k]=model.addVar(data.readyTime[i],data.dueTime[i],vtype=GRB.CONTINUOUS,name=f'S_{i}_{k}')forjinrange(data.nodeNum):X[i][j][k]=model.addVar(vtype=GRB.BINARY,name=f"X_{i}_{j}_{k}")#目标函数obj=LinExpr(0)foriinrange(data.nodeNum):forjinrange(data.nodeNum):ifi!=j:forkinrange(data.vehicleNum)bj.addTerms(data.distanceMatrix[i][j],X[i][j][k])model.setObjective(obj,GRB.MINIMIZE)#约束1：车辆只能从一个点到另一个点foriinrange(1,data.nodeNum-1):expr=LinExpr(0)forjinrange(data.nodeNum):ifi!=j:forkinrange(data.vehicleNum):ifi!=0andi!=data.nodeNum-1:expr.addTerms(1,X[i][j][k])model.addConstr(expr==1)#约束2：车辆必须从仓库出发forkinrange(data.vehicleNum):expr=LinExpr(0)forjinrange(1,data.nodeNum):expr.addTerms(1,X[0][j][k])model.addConstr(expr==1)#约束3：车辆经过一个点就必须离开一个点forkinrange(data.vehicleNum):forhinrange(1,data.nodeNum-1):expr1=LinExpr(0)expr2=LinExpr(0)foriinrange(data.nodeNum):ifh!=i:expr1.addTerms(1,X[i][h][k])forjinrange(data.nodeNum):ifh!=j:expr2.addTerms(1,X[h][j][k])model.addConstr(expr1==expr2)#约束4：车辆最终返回仓库forkinrange(data.vehicleNum):expr=LinExpr(0)foriinrange(data.nodeNum-1):expr.addTerms(1,X[i][data.nodeNum-1][k])model.addConstr(expr==1)#约束5：车辆容量约束forkinrange(data.vehicleNum):expr=LinExpr(0)foriinrange(1,data.nodeNum-1):forjinrange(data.nodeNum):ifi!=0andi!=data.nodeNum-1andi!=j:expr.addTerms(data.demand[i],X[i][j][k])model.addConstr(expr