带时间窗的多重运输调度问题的自适应交易策略

带时间窗的多重运输调度问题的自适应交易策略

1mvrp体系化方法在大规模生产过程中，多因素运输（mvr）的应用有着广泛的前景，关于文献进行了深入的研究[5、6、7、8、9、10、11和12]。由于场地限制、成本控制等原因,很多实际应用(如汽车装配自动线上的零部件供应和JIT生产方式的原材料、在制品供应)对货物的送达时间也有着严格的要求,这便提出了带时间窗的多重运输调度问题(MultipleDemandsVehicleRoutingProblemwithTimeWindows,MVRPTW)。MVRPTW的网络模型如下:设G=(V,E)是一个边赋正权的连通无向简单图,V={v1,v2,…,vn}是顶点集,E为边集,Kij(i,j=1,2,…,n)是一组指标值(自然数)。运输任务定义如下:设tkij、Tkij(i,j=1,2,…,n;k=1,2,…,Kij;0≤tkij≤Tkij)是常数,以[tkij,Tkij]wkij(wkij≥0整数)表示在时间区间[tkij,Tkij],(即不早于tkij时间、不迟于Tkij时间)内供应顶点vi需要向需求顶点vj提供wkij个运输单位的货物。要求:完成全部运输任务,每辆车在运输完毕后回到出发点。问题:确定每辆车的行车路线,使总的行驶路程最小的方案。这里1个运输单位是指一整车货物,一个运输任务是指车辆需把1个运输单位的货物从供应顶点运送到需求顶点的一件任务。一辆车从车队驶出、完成若干个运输任务、最终回到车队的过程构成了车辆的一次出行。为了方便,将[tkij,Tkij]1(即wkij=1)简记为[tkij,Tkij]。对于wkij>1,把[tkij,Tkij]分解为wkij个[tkij,Tkij],这样的分解不但使MVRPTW在形式上变得简单、也是实际运输调度工作的含义所在:调度员不但为每一辆车分配运输任务、给出行车路线,也要计算出每个运输任务的开始时间、完成时间。因此,下文假定wkij≤1(i,j=1,2,…,n;k=1,2,…,Kij)。本文首先给出了MVRPTW的转换模型和相关的概念,基于该模型提出一台车辆的MVRPTW的自适应TabuSearch算法,该算法具有较多的智能特征,例如对搜索过程进行综合记忆、自动确定施行强化和多样化策略的时机。通过对大量随机产生的MVRPTW进行上机计算,结果表明算法的性能令人满意。文中未定义的图论术语和记号与文献一致。2模型重建2.1任务的时间限制tkij为了处理方便,将MVRPTW网络模型作如下转换:从G=(V,E)出发构造一个新的混合型网络G1=(V,E,A),其中A是有向边(弧)的集合,在G中若从顶点vi到顶点vj有运输任务[tkij,Tkij],则在G1中作一条从顶点vi到顶点vj的弧,弧的长度为vi到vj的距离。为了表示该弧与运输任务[tkij,Tkij]的联系,在此弧上方靠近箭头方向标出货物应到达的时间区间[tkij,Tkij]。将两个顶点之间的距离用时间表示。例如,如果运输任务的时间区间用小时作量纲,则顶点间的距离也用车辆的行驶小时数来表示。此外,不考虑车辆的装货时间和卸货时间,这样做并没有失去普遍性。根据以上转换,设A={a1,a2,…,am}是运输任务的全体(不失一般性,设有m个运输任务),以[tLi,tRi]表示任务ai的时间限制。并假设如果车辆提前到达,则需要等待,且等待无需支付费用。2.2顶和顶出位置定义,v共3.定义1设G=(V,E,A)是一个混合型网络,ai,aj∈A是两条有向边,则ai到aj的距离为ai的箭头顶点到aj的箭尾顶点的距离,记为d(ai,aj)。这里,若ai=(x,y),则称x为ai的箭尾顶点,y为ai的箭头顶点。定义2[3.5]设G=(V,E,A)是一个混合型网络,s={a1,a2,…,ak}是一组有向边的有序集合,则称k-1∑i=1[d(ai,ai+1)+|ai|]+|ak|为s的长度,记作len(s)或len(a1,a2,…,ak)。这里|ai|为弧ai的长度。定义3[3.5]设G=(V,E,A)是一个混合型网络,s={a1,a2,…,ak}是一组有向边的有序集合,v是V中的任意顶点,定义v到s的距离为v到a1的箭尾顶点的距离,记为d(v,s)或d(v;a1,a2,…,ak);s到v的距离定义为ak箭头顶点到v的距离。记为d(s,v)或d(a1,a2,…,ak;v)。2.3求解mvrp体的可行性解定义4车辆将所有货物从供应顶点运送到需求顶点并返回车队的任何方案称为MVRPTW的可行解,在不引起混淆时简称为解。然而,车辆必须在指定的时间区间内将货物运到需求顶点,因此应有足够的车辆可供使用。一般地可供使用的车辆越多(不一定全部投入运行),则总的行驶路程可望越短。从这个意义上讲,MVRPTW需要考虑两个目标,其一是使用尽量少的车辆,其二是使总的行驶路程最短。但是MVRPTW这种两目标问题可以分解为多项式个固定车辆数的MVRPTW。因此本文考虑固定车辆数的MVRPTW。由于车辆数目有限,故可能出现某些解违反时间限制的情况。因提前可以等待,因此只需考虑任务拖期的情形。定义5设s是MVRPTW的可行解、ai的实际完成时间记为ti(s)。以g(s)=∑1≤i≤m,ti(s)>tRiCi(ti(s)-tRi)(1)表示对应s任务拖期的总惩罚,其中Ci是ai拖期时的单位时间惩罚(损失)。设有λ台车辆可供使用,车队的位置在顶点v1(为了醒目,亦用p表示车队位置)。又设A=A1∪A2∪…∪Aλ,其中Ai(i=1,2,…,m)为车辆i的任务集且两两不交,si为Ai的一个排序(即车辆i任务的完成顺序,从而确定了每个任务的开始、完成时间,即si隐含着车辆i的任务集及其完成顺序。任务完成时间的计算见3.1节),故求MVRPTW的可行解s等价于确定s1,s2,…,sλ,形式上表示为s=s1+s2+…+sλ。定义6设s是MVRPTW的可行解,s=s1+s2+…+sλ,si是车辆i的任务集,以f(x)=λ∑i=1d(p,si)+len(si)+len(si,p)(2)表示全部车辆总的行驶路程。以后为了方便,简记lenp(si)=d(p,si)+len(si)+d(si,p)。定义71)设s′、s″是MVRPTW的可行解,以(g,f)表示MVRPTW的目标函数向量。如果(g(s′),f(s′))≤(g(s″),f(s″)),则称s′的质量不比s″劣,记作s′⪯s″或s″⪰s′;如果(g(s′),f(s′))<(g(s″),f(s″)),则称s′的质量比s″好,记作s′≺s″或s″≻s′。2)设s*是MVRPTW的可行解,如果对一切可行解s有(g(s*),f(s*))≤(g(s),f(s)),则称s*为MVRPTW的最优解。据上述分析,MVRPTW最终转化为MVRPTW转换模型使用λ台车辆,求可行解s,在(1)达到最小值的前提下,使(2)达到最小值。本文的余下讨论是对该转换模型进行的。3mvrp体的可解本节研究一台车辆的情形,它是解决多车辆MVRPTW的基础。不难看出一台车辆MVRPTW的可行解与A中全体元素的排序形成一一对应关系,故下文将A中的全体元素的一个排序作为一个可行解。因此一台车辆MVRPTW就是求A中全体元素的一个排序s,在(1)达到最小值的前提下使得lenp(s)最小化。引入虚任务a0=p、am+1=p,它们的箭头、箭尾顶点均为p。任务a0、am+1无时间限制。设s是可行解,置车辆出发时间为t0(s)=0,返回车队时间记为tm+1(s)。3.1ts与模拟退火算法TabuSearch算法(简称TS,直译为禁忌搜索算法),由Glover提出,最近在组合优化方面获得了成功应用,不断得到完善[3,7,9,17,18,19,20]。TS是一个启发式算法,为了避免在局部最优点产生搜索终止,TS允许解的质量变坏,从这个意义上讲TS与模拟退火算法(SimulatedAnnealing,简称SA)的思想是相同的,也与中国古代兵家的“以退为进”、“将欲取之,必姑与之”战略思想不谋而合。然而,为了实现这样的战略,TS和SA采取了不同的手法。TS主要使用禁忌表(Tabu-list,简称表。因此,TS意译为表搜索算法)。为了提高效率,TS使用强化(intensification)策略和多样化(diversification)策略。TS通过设置控制参数对计算过程进行控制,如表结构、表长度、表条件(禁用条件)、表刷新、算法终止条件等。下面分别构造这些控制参数。3.1.1改进向量dtij,dki设s是当前解,s=(a1′,a2′,…,am′),第i位置的元素为ai′。因为ti′(s)≥max{t(i-1)′(s)+d(a(i-1)′,ai′)+|ai′|,tLi′},而由(1)知ti′(s)越小越好,于是取ai′的完成时间为ti′(s):=max{t(i-1)′(s)+d(a(i-1)′,ai′)+|ai′|,tLi′}(i=1,2,…,m),其中t0′(s):=0系车辆出发时间。对于一切i,j(i=1,2,…,m-1;j=i+1,…,m),计算改进向量(DTi′j′,DPi′j′):1)时间惩罚改进值DΤi′j′:=g(s)-g(s′)=∑i≤k≤m,tk′(s)>tRk′Ck′(tk′(s)-tRk′)-∑i≤k≤m,tk′(s′)>tRk′Ck′)tk′(s′)-tRk′)其中s′是在s中将ai′、aj′的位置交换所得。2)距离改进值如果j=i+1则DPi′j′:=[d(a(i-1)′,ai′)+d(ai′,aj′)+d(aj′,a(j+1)′)]-[d(a(i-1)′,aj′)+d(aj′,ai′)+d(ai′,a(j+1)′)]如果j≠i+1则DPi′j′:=[d(a(i-1)′,ai′)+d(ai′,a(i+1)′)+d(a(j-1)′,aj′)+d(aj′,a(j+1)′)]-[d(a(i-1)′,aj′)+d(aj′,a(i+1)′)+d(a(j-1)′,ai′)+d(ai′,a(j+1)′)];迭代方法:寻求不违反表条件的使(DTi′j′,DPi′j′)最大的一对ai′、aj′,在s中交换ai′、aj′的位置得到新的s;此时也称ai′、aj′进行了交换活动。3.1.2搜索方向的确定由三个表T、T′、T″组成,分述如下。1)表T表T是对ai、aj近期交换活动的记忆,定义为ai、aj禁止交换的时间(次数),T为m×m矩阵(tij)m×m,TS实际上只关心T的主对角线以上的部分。构造如下:设T的表长度为L(L≥0,整数),则0≤tij≤L,初始tij:=0(i,j=1,2,…,m)。2)表T′设s,s′,s″是迭代过程中相继出现的三个解,考虑下述四种情况:①s′≺s,s′≺s″:s′是局部极小。②s′≻s,s′≻s″:s′是局部极大。③s′≺s,s′≻s″:s′比s好,但比s″差。搜索进入″下山″(解的质量连续变好)步调。④s′≻s,s′≺s″:s′比s差,但比s″好,搜索进入″上山″(解的质量连续变差)步调。表T′是对局部极小、极大点的记忆,即对情形①、②进行记忆。这两种情形为改变搜索方向提供机会,而情形③、④则不能。对于情形①、②,设s″是由s′交换ai、aj的位置得到,则T′的元素为三元对(s′,ai,aj)(由于交换的对称性,或者说(s′,ai,aj)和三元对(s′,aj,ai)是等价的,构造三元对(s′,ai,aj)时应满足i<j)。直观看来,此时ai、aj的交换决定着以后“上山”、“下山”的方向。事实上,在全局最优化中,确定搜索的大方向是最重要的,因为方向性错误将导致搜索无功而返。初始T′:=\～φ,T′的表长为L′(L′≥0,整数)。3)表T″表T″是对已经搜索和尚未搜索地区的记忆,T″为m×m矩阵T=(t″ij)m×m。设初始解为s,则T″的初值定义如次:对于一切i,j=1,2,…,m,如果ai在s中的位置为j,则t″ij:=1,否则t″ij:=0。3.1.3ai、aj被禁止的交换活动设s是当前解,如果以下两个条件之一成立,则ai、aj被禁止作交换(即迭代)活动:1)tij>0;2)tij=0但(s,ai,aj)在T′中。3.1.4t元素的更新如果当前迭代中,选中了ai、aj作交换,得到的解记作s′。1)T的表刷新tij:=L;对于一切i′,j′=1,2,…,m,(i′,j′)≠(i,j),ti′j′:=max{0,ti′j′-1}。2)T′的表刷新如果s′是局部极小点或局部极大点(必须等到下一次迭代完成后才能确定s′是否是极大、极小点),则若|T′|<L′,则将(s′,ai,aj)加到T′的尾部(即(s′,ai,aj)为T′的最后一个元素);如果|T′|=L′,则去掉T′的第1个元素,再将(s′,ai,aj)加到T′的尾部。3)T″的表刷新对于一切i′,j′=1,2,…,m,如果ai′在s′中的位置为j′,则t″i′j′:=t″i′j′+1,否则t″i′j′:=t″i′j′。某个地区未被充分搜索是指某个t″i′j′过小,某个地区被过分搜索是指某个t″i′j′过大,不充分或过分搜索都可能导致不能找到全局最优解的危险。注意到若迭代进行了k-1次,则包括初始解在内有k个可行解出现了。此时T″全部元素的和为km,全部元素的平均值为k/m。定义AverT″:=k/m,每次迭代之后需对AverT″进行更新。AverT″为度量T″元素过大、过小的尺度。AverT″亦可考虑其它的取值方式。3.1.53.1.6t的添加位置1)T的相伴表刷新对于一切i、j(i,j=1,2,…,m),如果tij>0,则tij:=max{0,tij-L+L1};其中L、L1分别为T变换前后的表长度。2)T′的相伴表刷新如果|T′|>L′1,则T′:={将T′的前|T′|-L′1)个元素去掉后的剩余部分},即新的T′只保留后L′1个元素;其中L′1为T′变换后的表长度。3.1.7把ai导入位置l如果解连续不能改进,且下列情形之一发生,则施行多样化策略:1)某个t″ij过小,则把ai插入到位置j。2)某个t″ij过大,则把ai插入到位置j′,其中t″ij′=min{t″ij|j=1,2,…m}。过大、过小的尺度如前述。3.1.8算法结束了如果在迭代过程中,解连续不能改进的次数达到阀值,则算法终止。3.2强化、多样性策略算法1(MVRPTW的TabuSearch算法)初始化:s*:={a1,a2,…,am};//s*为目前所获得的最好解s:=s*;//s是当前解iter_improved_num:=0;//解连续改进时迭代计数输入非负整数iter_improved_max;//施行强化策略时机iter_bad_num:=0;//解连续不能改进时迭代计数输入非负整数iter_bad_max_int_dev;//施行强化、多样化策略时机输入非负整数iter_bad_max_end;/*算法结束的时机。此数应比iter_bad_max_int_dev大一些*/输入非负整数L、L′;//表T、T′的长度初始化表T,T′,T″;Best_exchange_list:=\～φ;//最好的前L+L′+2个交换iter_int_dev:=0;//施行强化、多样化策略的等待计数输入非负整数iter_int_dev_max;//施行强化、多样化策略等待周期int_mode:=0;//选用强化策略或多样化策略的判别参数Whileiter_bad_num<iter_bad_max_enddo{i:=1;whilei≤m-1do{j:=i+1;whilej≤mdo{在s中分别找出第i,j位置的元素,不妨设为ai′,aj′;计算改进向量(DTi′j′,DPi′j′);//计算方法见3.1如果|Best_exchange_list|<L+L′+2则Best_exchange_list:=Best_exchange_list∪{(ai′,aj′)};否则如果Best_exchange_list中有Pi*j*比现行的Pi′j′小则Best_exchange_list:=Best_exchange_list∪{(ai′,aj′)}-{(ai*,aj*)};j:=j+1;}i:=i+1;}在Best_exchange_list中寻求不违反表条件的最好(ai′,aj′);交换ai′,aj′的位置得到新的s;刷新表T、T′、T″;Best_exchange_list:=ue001φ;如果(g(s),f(s))<(g(s*),f(s*))则{s*:=s;iter_bad_num:=0;iter_improved_num:=iter_improved_num+1;}否则{iter_bad_num:=iter_bad_num+1;iter_improved_num:=0;}如果iter_improved_num=iter_improved_max//施行强化策略{将L、L′的值缩小;//例如L:=L-1、L′:=L′-1相伴表刷新;iter_improved_num:=0;}iter_int_dev:=max(0,iter_int_dev-1);如果iter_int_dev=0且iter_bad_num≥iter_bad_max_int_dev{如果int_mode=0则//施行多样化策略{如果存在t″ij<AverT″则把ai*插入到位置j*;/*这里t″i*j*=min{t″ij|t″ij<AverT″,i,j=1,2,…,m}*/否则{如果存在t″ij>AverT″则把ai*插入到位置j′;/*这里t″i*j*=max{t″ij|t″ij>AverT″,i,j=1,2,…,m},t″i*j′=min{t″i*j|j=1,2,…,m}*/否则//若不能施行多样化策略,则施行强化策略int_mode:=1;}}如果int_mode=1//施行强化策略{将L、L′的值放大;//例如L:=L+1、L′:=L′+1相伴表刷新;}int_mode:=mod(int_mode+1,2);//取模2。强化与多样化策略交替进行iter_int_dev:=iter_int_dev_max;}}输出最优解s*;结束。4ti,b表征例1已知图1(a)所示网络G=(V,E)(边旁边的数字表示两个顶点之间以小时计的边长度),运输任务如表1。Ci=0(i=1,2,…,6)。转换的G1=(V,E,A)如图1(b)所示。今用算法1求解。控制参数取值为:iter_improved_max=4,iter_bad_max_int_dev=4,iter_bad_max_end=10,L=2,L′=2,iter_int_dev_max=4。控制变量初始化:iter_improved_num:=0,iter_bad_num:=0,Best_exchange_list:=ue001φ,iter_int_dev:=0,int_mode:=0。初始解s:={a1,a2,a3,a4,a5,a6},s*:=s。表初始化:tij:=0(i,j=1,2,…,6);T′:=ue001φ;t″ii:=1(i=1,2,…,6),t″ij:=0(i≠j;i,j=1,2,…,6)。迭代过程如下:1)对s={a1,a2,a3,a4,a5,a6},有f(s)=lenp(s)=30t0(s):=0,t1(s):=max{t0(s)+d(a0,a1)+|a1|,tL1}=max{0+4+2,9}=9t2(s):=max{t1(s)+d(a1,a2)+|a2|,tL2}=max{9+3+1,12}=13,……g(s)=64。DTij、DPij的计算结果见表2。Best_exchang_list:={(a2,a4),(a1,a3),(a2,a3),(a2,a5),(a1,a6),(a3,a5)}。在s中交换a2、a4的位置得到新的s:={a1,a4,a3,a2,a5,a6}。控制参数更新:iter_improved_num:=1,iter_b