禁忌搜索算法求解带时间窗的车辆路径规划问题详解(附Java代码)

VRPTW简介

VRPTW问题可描述为:假设一个配送中心为周围若干个位于不同地理位置、且对货物送达时间有不相同要求的客户点提供配送服务。其中,配送中心用于运行的车辆都是同一型号的(即拥有相同的容量、速度);配送中心对车辆出入的时间有限制。我们的任务是找出使所有车辆行使路径总和最小的路线。

VRPTW的更多详细介绍可以参考之前的推文:

干货|十分钟快速掌握CPLEX求解VRPTW数学模型(附JAVA代码及CPLEX安装流程

为了保持文章的独立型,同时方便后续讲解,这里给出建模实例(参考文献在文末标注):

所求的所有车辆路线需满足以下要求:

在此基础上求出每辆车辆的总时间最短(由于车辆速度相同,时间最短相当于路程最短)的路线。(允许不使用某些车辆)

TabuSearch简介

禁忌搜索算法(TabuSearchAlgorithm,简称TS)起源于对于人类记忆功能的模仿,是一种亚启发式算法(meta-heuristics)。它从一个初始可行解(initialfeasiblesolution)出发,试探一系列的特定搜索方向(移动),选择让特定的目标函数值提升最多的移动。为了避免陷入局部最优解,禁忌搜索对已经历过的搜索过程信息进行记录,从而指导下一步的搜索方向。

禁忌搜索是人工智能的一种体现,是局部搜索的一种扩展。禁忌搜索是在邻域搜索(localsearch)的基础上,通过设置禁忌表(tabulist)来禁忌一些曾经执行过的操作,并利用藐视准则来解禁一些优秀的解。

TS+VRPTW

对邻域搜索类算法而言,采取的搜索算子和评价函数至关重要。下面详细介绍代码中针对VRPTW的插入算子和评价函数。

插入算子:

评价函数:

算法概述

Java代码详解

代码主要分为以下几个类:

Main,主函数;

CustomerType,存放客户节点的信息;

RouteType,存放车辆路线信息;

Parameter,存放全局变量;

EvaluateRoute,处理路线方法;

InitAndPrint,初始化与输出对应方法;

TS,禁忌搜索方法。

接下来分别介绍。

Main

publicclassMain{publicstaticvoidmain(Stringarg[]){longbegintime=System.nanoTime();ReadIn();Construction();TabuSearchnewnew();Output();longendtime=System.nanoTime();doubleusedTime=(endtime-begintime)/(1e9);System.out.println();System.out.println(“程序耗时:”+usedTime+”s”);}}

程序的入口。

CustomerType

publicclassCustomerType{intNumber;//节点自身编号intR;//节点所属车辆路径编号doubleX,Y;//节点横纵坐标doubleBegin,End,Service;//节点被访问的最早时间,最晚时间以及服务时长doubleDemand;//节点的需求容量publicCustomerType(){this.Number=0;this.R=0;this.Begin=0;this.End=0;this.Service=0;this.X=0;this.Y=0;this.Demand=0;}publicCustomerType(CustomerTypec1){this.Number=c1.Number;this.R=c1.R;this.Begin=c1.Begin;this.End=c1.End;this.Service=c1.Service;this.X=c1.X;this.Y=c1.Y;this.Demand=c1.Demand;}}

客户类,对图中的每一个客户,分别构建客户类,存放自身编号,所属车辆路线,坐标位置,访问时间窗,服务所需时长、需求。

RouteType

importjava.util.ArrayList;publicclassRouteType{publicdoubleLoad;//单条路径承载量publicdoubleSubT;//单条路径违反各节点时间窗约束时长总和publicdoubleDis;//单条路径总长度publicArrayList

V=newArrayList<>();//单条路径上顾客节点序列。在route中,第0个、最后一个都为depot,第k个为第k位。}

路线类,记录该路线的总承载量,总长度,对时间窗约束的总违反量,以及单条路径上的客户节点序列。

Parameter

publicclassParameter{publicstaticdoubleINF=Double.MAX_VALUE;publicstaticintCustomerNumber=100;//算例中除仓库以外的顾客节点个数publicstaticintVehicleNumber=25;publicstaticintCapacity=200;//车辆最大容量publicstaticintIterMax=2000;//搜索迭代次数publicstaticintTabuTenure=20;//禁忌步长publicstaticint[][]Tabu=newint[CustomerNumber+10][VehicleNumber+10];//禁忌表用于禁忌节点插入操作:[i][j]将节点i插入路径j中publicstaticint[]TabuCreate=newint[CustomerNumber+10];//禁忌表用于紧急拓展新路径或使用新车辆publicstaticdoubleAns;//最优解距离publicstaticdoubleAlpha=1,Beta=1,Sita=0.5;//Alpha,Beta为系数,计算目标函数值;Sita控制系数改变的速度publicstaticdouble[][]Graph=newdouble[CustomerNumber+10][CustomerNumber+10];//记录图publicstaticCustomerType[]customers=newCustomerType[CustomerNumber+10];//存储客户数据publicstaticRouteType[]routes=newRouteType[CustomerNumber+10];//存储当前解路线数据publicstaticRouteType[]Route_Ans=newRouteType[CustomerNumber+10];//存储最优解路线数据}

参数类,有关VRPTW和TS的变量都存储在这里,在这里修改数据。

EvaluateRoute

publicstaticbooleanCheck(RouteTypeR[]){//检验解R是否满足所有约束doubleQ=0;doubleT=0;//检查是否满足容量约束for(inti=1;i<=VehicleNumber;++i)if(R[i].V.size()>2&&R[i].Load>Capacity)//对有客户且超过容量约束的路径,记录超过部分Q=Q+R[i].Load-Capacity;//检查是否满足时间窗约束for(inti=1;i<=VehicleNumber;++i)T+=R[i].SubT;//分别根据约束满足的情况和控制系数Sita更新Alpha和Beta值//新路径满足条件,惩罚系数减小,//新路径违反条件,惩罚系数加大。if(Q==0&&Alpha>=0.001)Alpha/=(1+Sita);elseif(Q!=0&&Alpha<=2000)Alpha*=(1+Sita);if(T==0&&Beta>=0.001)Beta/=(1+Sita);elseif(T!=0&&Beta<=2000)Beta*=(1+Sita);if(T==0&&Q==0)returntrue;elsereturnfalse;}

check函数是对产生解的检验。

由于插入算子产生的解并不都满足所有约束条件,对局部搜索产生的较优解需要判断是否满足时间窗约束和容量约束后,再决定是否为可行解。

在check局部最优解的过程中,修改惩罚系数Alpha、Beta的值。

publicstaticvoidUpdateSubT(RouteTyper){//更新路径r对时间窗的违反量doubleArriveTime=0;for(intj=1;j<r.V.size();++j){//对每一个节点分别计算超出时间窗的部分ArriveTime=ArriveTime+r.V.get(j-1).Service//服务时间+Graph[r.V.get(j-1).Number][r.V.get(j).Number];//路途经过时间if(ArriveTime>r.V.get(j).End)//超过,记录r.SubT=r.SubT+ArriveTime-r.V.get(j).End;elseif(ArriveTime<r.V.get(j).Begin)//未达到,等待ArriveTime=r.V.get(j).Begin;}}

UpdateSubT函数更新一条车辆路线中在每一个客户点的时间窗违反量。通过遍历整条路线累加得到结果。

//计算路径规划R的目标函数值,通过该目标函数判断解是否较优publicstaticdoubleCalculation(RouteTypeR[],intCus,intNewR){//目标函数主要由三个部分组成:D路径总长度(优化目标),Q超出容量约束总量,T超出时间窗约束总量//目标函数结构为f(R)=D+Alpha*Q+Beta*T,第一项为问题最小化目标,后两项为惩罚部分//其中Alpha与Beta为变量,分别根据当前解是否满足两个约束进行变化,根据每轮迭代得到的解在Check函数中更新doubleQ=0;doubleT=0;doubleD=0;//计算单条路径超出容量约束的总量for(inti=1;i<=VehicleNumber;++i)if(R[i].V.size()>2&&R[i].Load>Capacity)Q=Q+R[i].Load-Capacity;//计算总超出时间for(inti=1;i<=VehicleNumber;++i)T+=R[i].SubT;//计算路径总长度for(inti=1;i<=VehicleNumber;++i)D+=R[i].Dis;return(D+Alpha*Q+Beta*T);//返回目标函数值}

Calculate函数计算目标函数值,惩罚部分累加后乘以惩罚系数。

InitAndPrint

//计算图上各节点间的距离privatestaticdoubleDistance(CustomerTypeC1,CustomerTypeC2){returnsqrt((C1.X-C2.X)*(C1.X-C2.X)+(C1.Y-C2.Y)*(C1.Y-C2.Y));}

根据计算距离。

//读取数据publicstaticvoidReadIn(){for(inti=0;i<customernumber+10;i++){customers[i]=newcustomertype();routes[i]=newroutetype();route_ans[i]=newroutetype();}try{scannerin=newscanner(newfilereader(“c101.txt”));for(inti=1;i<=customernumber=””i=””.number=”in.nextint()+1;”.x=”in.nextDouble();customers[i].Y=in.nextDouble();customers[i].Demand=in.nextDouble();customers[i].Begin=in.nextDouble();customers[i].End=in.nextDouble();customers[i].Service=in.nextDouble();}in.close();}catch(FileNotFoundExceptione){//FilenotfoundSystem.out.println(“filenotfound!”);system.exit(-1);}for(inti=”1;i<=VehicleNumber;++i){if(routes[i].V.size()!=0)routes[i].V.clear();routes[i].V.add(newCustomerType(customers[1]));//尝试往这里加入一个复制,后面也都要改。routes[i].V.add(newCustomerType(customers[1]));routes[i].V.get(0).End=routes[i].V.get(0).Begin;//起点routes[i].V.get(1).Begin=routes[i].V.get(1).End;//终点//算例中给出节点0有起始时间0和终止时间,所以如上赋值。routes[i].Load=0;}Ans=INF;for(inti=1;i<=CustomerNumber+1;++i)for(intj=1;j<=CustomerNumber+1;++j)Graph[i][j]=Distance(customers[i],customers[j]);}

从文件中读取算例(在此修改算例,记得同时修改Parameter类中的参数),并对当前解routes[]的每条路线进行初始化,起终点都为配送中心。

记录客户间距离,存储在Graph数组中。

//构造初始解publicstaticvoidConstruction(){int[]Customer_Set=newint[CustomerNumber+10];for(inti=1;i<=CustomerNumber;++i)Customer_Set[i]=i+1;intSizeof_Customer_Set=CustomerNumber;intCurrent_Route=1;//以满足容量约束为目的的随机初始化//即随机挑选一个节点插入到第m条路径中,若超过容量约束,则插入第m+1条路径//且插入路径的位置由该路径上已存在的各节点的最早时间决定while(Sizeof_Customer_Set>0){intK=(int)(random()*Sizeof_Customer_Set+1);intC=Customer_Set[K];Customer_Set[K]=Customer_Set[Sizeof_Customer_Set];Sizeof_Customer_Set–;//将当前服务过的节点赋值为最末节点值,数组容量减1//随机提取出一个节点,类似产生乱序随机序列的代码if(routes[Current_Route].Load+customers[C].Demand>Capacity)Current_Route++;//不满足容量约束,下一条车辆路线for(inti=0;i<routes[Current_Route].V.size()-1;i++)//对路径中每一对节点查找,看是否能插入新节点if((routes[Current_Route].V.get(i).Begin<=customers[C].Begin)&&(customers[C].Begin<=routes[Current_Route].V.get(i+1).Begin)){routes[Current_Route].V.add(i+1,newCustomerType(customers[C]));//判断时间窗开始部分,满足,则加入该节点。routes[Current_Route].Load+=customers[C].Demand;customers[C].R=Current_Route;//更新路径容量,节点类。break;}}//初始化计算超过时间窗约束的总量for(inti=1;i<=VehicleNumber;++i){routes[i].SubT=0;routes[i].Dis=0;for(intj=1;j<routes[i].V.size();++j){routes[i].Dis+=Graph[routes[i].V.get(j-1).Number][routes[i].V.get(j).Number];}UpdateSubT(routes[i]);}}

Construction构造初始解。根据前文伪代码构造初始解,每次随机选择节点(类似打乱有序数列)。

针对该节点找到符合容量约束,同时时间窗开启时间符合要求的位置,插入该节点。记得在插入节点时同时更新该节点所属的路径。

对时间窗违反量进行初始化。

publicstaticvoidOutput(){//结果输出System.out.println(“************************************************************”);System.out.println(“TheMinimumTotalDistance=”+Ans);System.out.println(“ConcreteScheduleofEachRouteasFollowing:”);intM=0;for(inti=1;i<=VehicleNumber;++i)if(Route_Ans[i].V.size()>2){M++;System.out.print(“No.”+M+”:”);for(intj=0;j<Route_Ans[i].V.size()-1;++j)System.out.print(Route_Ans[i].V.get(j).Number+”->”);System.out.println(Route_Ans[i].V.get(Route_Ans[i].V.size()-1).Number);}System.out.println(“************************************************************”);}

输出结果,不再赘述。

publicstaticvoidCheckAns(){//检验距离计算是否正确doubleCheck_Ans=0;for(inti=1;i<=VehicleNumber;++i)for(intj=1;j<Route_Ans[i].V.size();++j)Check_Ans+=Graph[Route_Ans[i].V.get(j-1).Number][Route_Ans[i].V.get(j).Number];System.out.println(“Check_Ans=”+Check_Ans);//检验是否满足时间窗约束booleanflag=true;for(inti=1;i<=VehicleNumber;i++){UpdateSubT(Route_Ans[i]);if(Route_Ans[i].SubT>0)flag=false;}if(flag)System.out.println(“Solutionsatisfiestimewindowsconstruction”);elseSystem.out.println(“Solutionnotsatisfiestimewindowsconstruction”);}

最后加入一个CheckAns函数,检验一下输出解是否满足时间窗约束,计算的距离是否正确,有备无患~

TS

privatestaticvoidaddnode(intr,intpos,intCus){//节点加入的路径routes[r],节点customer[Cus],节点加入路径的位置pos//更新在路径r中加上节点Cus的需求routes[r].Load+=customers[Cus].Demand;//更新在路径r中插入节点Cus后所组成的路径距离routes[r].Dis=routes[r].Dis-Graph[routes[r].V.get(pos-1).Number][routes[r].V.get(pos).Number]+Graph[routes[r].V.get(pos-1).Number][customers[Cus].Number]+Graph[routes[r].V.get(pos).Number][customers[Cus].Number];//在路径r中插入节点Cusroutes[r].V.add(pos,newCustomerType(customers[Cus]));//插入i到下标为l处}privatestaticvoidremovenode(intr,intpos,intCus){//节点去除的路径routes[r],节点customer[cus],节点所在路径的位置pos//更新在路径r中去除节点Cus的需求routes[r].Load-=customers[Cus].Demand;//更新在路径r中去除节点Cus后所组成的路径的距离routes[r].Dis=routes[r].Dis-Graph[routes[r].V.get(pos-1).Number][routes[r].V.get(pos).Number]-Graph[routes[r].V.get(pos).Number][routes[r].V.get(pos+1).Number]+Graph[routes[r].V.get(pos-1).Number][routes[r].V.get(pos+1).Number];//从路径r中去除节点Cusroutes[r].V.remove(pos);}

先是两个辅助函数,addnode和removenode,他们是插入算子的执行部分。

publicstaticvoidTabuSearch(){//禁忌搜索//采取插入算子,即从一条路径中选择一点插入到另一条路径中//在该操作下形成的邻域中选取使目标函数最小化的解doubleTemp1;doubleTemp2;//初始化禁忌表for(inti=2;i<=CustomerNumber+1;++i){for(intj=1;j<=VehicleNumber;++j)Tabu[i][j]=0;TabuCreate[i]=0;}intIteration=0;while(Iteration<IterMax){intBestC=0;intBestR=0;intBestP=0;intP=0;doubleBestV=INF;for(inti=2;i<=CustomerNumber+1;++i){//对每一个客户节点for(intj=1;j<routes[customers[i].R].V.size();++j){//对其所在路径中的每一个节点if(routes[customers[i].R].V.get(j).Number==i){//找到节点i在其路径中所处的位置jP=j;//标记位置break;}}removenode(customers[i].R,P,i);//将客户i从原路径的第P个位置中移除//找到一条路径插入删去的节点for(intj=1;j<=VehicleNumber;++j)for(intl=1;l<routes[j].V.size();++l)//分别枚举每一个节点所在位置if(customers[i].R!=j){addnode(j,l,i);//将客户l插入路径j的第i个位置Temp1=routes[customers[i].R].SubT;//记录原先所在路径的时间窗违反总和Temp2=routes[j].SubT;//记录插入的路径时间窗违反总和//更新i节点移出的路径:routes[customers[i].R].SubT=0;UpdateSubT(routes[customers[i].R]);//更新i节点移入的路径j:routes[j].SubT=0;UpdateSubT(routes[j]);doubleTempV=Calculation(routes,i,j);//计算目标函数值if((TempV<Ans)||//藐视准则,如果优于全局最优解(TempV<BestV&&//或者为局部最优解,且未被禁忌(routes[j].V.size()>2&&Tabu[i][j]<=Iteration)||(routes[j].V.size()==2&&TabuCreate[i]<=Iteration)))//禁忌插入操作,前者为常规禁忌表,禁忌插入算子;后者为特殊禁忌表,禁忌使用新的车辆//路径中节点数超过2,判断是否禁忌插入算子;路径中只有起点、终点,判断是否禁忌使用新车辆。if(TempV<BestV){//记录局部最优情况BestV=TempV;//bestvehicle所属车辆BestC=i;//bestcustomer客户BestR=j;//bestroute所属路径BestP=l;//bestposition所在位置}//节点新路径复原routes[customers[i].R].SubT=Temp1;routes[j].SubT=Temp2;removenode(j,l,i);}//节点原路径复原addnode(customers[i].R,P,i);}//更新车辆禁忌表if(routes[BestR].V.size()==2)TabuCreate[BestC]=Iteration+2*TabuTenure+(int)(random()*10);//更新禁忌表Tabu[BestC][customers[BestC].R]=Iteration+TabuTenure+(int)(random()*10);//如果全局最优的节点正好属于当前路径,过for(inti=1;i<routes[customers[BestC].R].V.size();++i)if(routes[customers[BestC].R].V.get(i).Number==BestC){P=i;break;}//依据上述循环中挑选的结果,生成新的总体路径规划//依次更新改变过的路径的:载重,距离长度,超出时间窗重量//更新原路径removenode(customers[BestC].R,P,BestC);//更新新路径addnode(BestR,BestP,BestC);//更新超出时间routes[BestR].SubT=0;UpdateSubT(routes[BestR]);routes[customers[BestC].R].SubT=0;UpdateSubT(routes[customers[BestC].R]);//更新被操作的节点所属路径编号customers[BestC].R=BestR;//如果当前解合法且较优则更新存储结果if((Check(routes)==true)&&(Ans>BestV)){for(inti=1;i<=VehicleNumber;++i){Route_Ans[i].Load=routes[i].Load;Route_Ans[i].V.clear();for(intj=0;j<routes[i].V.size();++j)Route_Ans[i].V.add(routes[i].V.get(j));}Ans=BestV;}Iteration++;}}

终!于!到了最后一步了!

现在万事俱备,只欠东风,只需要按照禁忌搜索的套路,将所有工具整合在一起,搭建出代码框架就ok啦。

由于我们采用routes[]数组存储当前解,因此在进行插入操作之前要存储部分数据,在计算完目标函数之后要进行复原操作。

在更新禁忌表时,对禁忌步长的计算公式可以灵活改变。

记得对局部最优解进行判断,再选取为可行的全局最优解。

算例展示

我们采用标准solomon测试数据c101.txt进行测试。(算例可在留言区下载获取)

VehicleNumber=25;

Capacity=200;

分别测试节点数25,50,100的情况。精确解分别为:

CustomerNumber=25:

CustomerNumber=50:

CustomerNumber=100:

可见我们的代码精确度还是很可以滴~~

当然不排除运气不好,得出很差解的情况。不相信自己的人品可以手动调整迭代次数IterMax。

本期的内容到这里就差不多结束了!开心!

在这里提醒大家一下,在针对启发式算法的学习过程中,编写代码的能力是很重要的。VRPTW是一个很好的载体,建议有时间的读者尽量将学到的算法知识运用到实践中去。小编会和你们一起学习进步的!

参考文献:

Cordeau,J.F.,Laporte,G.,&Mercier,A..(2001).Aunifiedtabusearchheuristicforvehicleroutingproblemswithtimewindows.JournaloftheOperationalResearchSociety,52(8),928-936.

暂无评论

您必须登录才能参与评论!
立即登录
暂无评论...