带时间窗和调货特性的多品类取送货优化算法

温 昆 ,郭 鹏,2 ,裴 霞 ,吴 晓,2

(1.西南交通大学 机械工程学院,成都 610031;2.轨道交通运维技术与装备四川省重点实验室,成都 610031)

快时尚品(譬如女鞋、女装等)款式更新快且上架展示周期短,需要零售管理方主动控制门店库存,尽可能保证门店销售需求得到满足。管理人员基于门店以往的销售数据对时下流行货品进行销量预测,估算出每种货品未来时段的销售量。以此为基础并结合每个门店现有的库存,对每个门店的每种货品进行补货和调货。当预测的销量小于现有的存量时,为了防止门店此类货品积压,需及时将多余的货品运走;当预测的销量大于现有的存量时,需对门店的此类货品进行补充。提倡门店之间货品互补,以降低调配成本。此外,门店取送货服务时段存在差异。为保证取送货的顺利进行,对配送车辆到达门店的时间也有限制。为此,本文提出带时间窗和调货特性的多品类取送货问题(Multi-commodity Pickup and Delivery Problem with Time Window and Transshipment,M-PDPTWT),以帮助快时尚零售企业降低成本、提高作业效率。M-PDPTWT是带时间窗车辆路径问题(Vehicle Routing Problem with Time Windows,VRPTW)更一般的描述。VRPTW 旨在考虑各个门店可配送的时间窗,在满足仓库辐射范围内各个门店货品需求和不超过车辆最大载容量的前提下,以实现配送车辆数最小化、车辆行驶距离最短、车辆装载均衡等目标。

Braekers等[1]对VRP进行了综述,讨论了精确算法与启发式算法应用情况。针对VRPTW,精确算法有基于集划分的分支定价与切割算法[2]、基于集合覆盖的列生成算法[3]、改进分支切割算法[4]等,启发式算法有遗传算法[5]、大规模邻域搜索算法[6]、进化学习算法[7]、融合禁忌搜索和人工蜂群算法的混合算法[8]、混合蚁群算法[9]以及混合遗传算法[10]等。VRPTW 作为经典的车辆路径问题,具有时间窗和车容量约束,但仅考虑送货需求。同时,带有取送货和时间窗的车辆路径问题(Vehicle Routing Problem with Simultaneous Pickup-Delivery and Time Windows,VRPSPDTW)和带时间窗的同时取送货问题(Pickup and Delivery Problem with Time Windows,PDPTW)均考虑了取送货需求。VRPSPDTW 和PDPTW 都是针对配送时间有要求的顾客,在满足车载容量、顾客时间窗以及每个顾客对货物的配送请求,规划出车辆数最少、运输距离最小以及仓库处理成本最小等目标的运输方案。VRPSPDTW 规定顾客需求的货物必须从仓库运输,顾客不需要的货物也必须带回仓库,如快递员上门取送包裹。Koç等[11]对同时取送货问题进行了综述,分析了取送货问题的研究趋势。模拟退火算法[12]、化学反应优化算法[13]、基于学习的两阶段搜索算法[14]等被用来求解VRPSPDTW,且面临最小化车辆数、最小化车辆行驶距离、最小化运营成本以及最大化工作效率等目标。PDPTW 则是将指定节点(包括仓库和顾客)的货物运送到另一指定节点,是一对一的运输问题,如经典的电话叫车问题(Dial A Ride Problem,DARP)。自适应大规模邻域搜索算法[15]、禁忌搜索算法[16]、基于数学规划的混合算法[17]、改进遗传算法[18]以及列生成算法与启发式规则相结合的CGA 混合算法[19]等用于求解不同场景下的PDPTW。

无论是VRPSPDTW 或PDPTW,配送请求都是针对一种货物,货物的配送和交付也不具有可重复使用性。而在快时尚服装销售邻域,需要考虑运送的货物多达几十、上百种,甚至上千种,这势必增加问题的求解难度。文献[20-21]中针对多货品调配和门店互补的车辆路径问题提出了精确与启发式求解算法。Li等[20]提出了一种基于强集合划分模型的分支定价切割算法。Zhang等[21]提出了一种基于自适应内存规划算法。对于快时尚品的配送货问题,栾玉麟等[22]基于顶点拆分策略,设计了遗传算法求解大规模实际算例。

上述文献中,与本文最为相关的文献[20-22]未考虑各个门店存在不同的取送货服务时间窗,而考虑VRPTW 的文献则未考虑多品类与门店间相互调货的实际情况。基于此,文本提出以车辆数最小化为第一优化目标,最小化转运成本(运输成本和仓库处理成本)为第二优化目标的M-PDPTWT,构建了混合整数规划模型。基于优化目标设计了两阶段搜索算法:第1阶段,利用最短路径搜索规则构建初始解;第2阶段,以变邻域搜索算法为框架。不同搜索时期用不同的适应度函数进行计算,以增强对最小化车辆数的搜索能力,并定义多种邻域解属性和优先保留策略。在搜索过程中不仅接受可行解,还允许接受具有潜在搜索能力的不可行解,以增强算法对全局解空间的搜索能力。在计算实验部分,采用标准算例和基于实际运营数据生成的算例进行了算法性能分析。

1 数学模型

1.1 问题描述与假设

配送车辆从仓库出发,携带不超过车辆最大载重量Q的多种货品用以补给缺少货品的门店。每个门店有规定的到达时间窗,配送车辆如果早于时刻ei到达则必须等待,但不能晚于li到达。车辆在达到某个门店后耗费停留时间sei用于装卸货品。每个门店除了有补货需求,也须将门店其他多余货品取走,避免货品在门店长期积压。被运走的货品可以用来补给同一条配送路径上后续门店,也可将其运回仓库。通过提倡货品在门店之间相互调配的方式,可以降低搭配成本和仓库对货品的处理成本。对整个配送网络体系的路径优化中,需要考虑的目标函数有两个:配送车辆数最少及在此基础上配送成本(距离)和仓库处理成本最小化。

建立模型时假设如下:

(1)每辆车载重量相同,行驶速度相同。

(2)货品装车时不考虑体积的影响。

(3)每个门店必须被访问,且只能访问一次。

1.2 符号说明及参数

假设有n个门店,每个门店取送货在同一地点、同一时段内完成。门店用i表示,仓库用0 和n+1表示。集合N为所有门店的集合,包含仓库的集合N0=N∪{0,n+1}。配送车辆集合V={1,2,…,h},其中,h为车场提供的车数量。K={1,2,…,m}为货品种类的集合,其中,m为配送货品的最大种类数目。

其余符号

cij——从门店i到门店j的运输成本,i,j∈N0

tij——从门店i到门店j的运输时间,i,j∈N0

sei——门店i的服务时间,i∈N0,其中,se0=sen+1=0

ei——门店i时间窗的最早达到时间,i∈N0,其中,e0=en+1表示仓库的最早工作时间

li——门店i时间窗的最晚到达时间,i∈N0,其中,l0=ln+1表示仓库的最晚工作时间

ω——仓库的处理成本系数

Q——车辆的最大载容量

决策变量

1.3 模型建立

目标函数式(1)表示最小化车辆数;式(2)表示最小化运输成本和仓库处理成本;约束条件式(3)、(4)表示每个门店必须访问且只能访问一次;式(5)为车流量守恒约束,表示第v辆车在访问p点之后必须从p点离开;式(6)表示车辆必须从仓库0点出发;式(7)表示车辆必须返回仓库n+1点;式(8)表示到达门店i的时间加上在门店i的服务时间再加上从门店i到达门店j的车辆行驶时间不能晚于门店j的时间窗最大值;式(9)为时间窗约束,表示车辆必须在规定的时间窗范围内访问门店,早于时间窗的最小值到达门店必须等待,不可晚于最晚的到达时间;式(10)为车辆的载重约束,表示车辆的实时载重不能超过车辆的最大载重;式(11)表示当第v辆车从i点到达j点时,离开i点时每类货品的载重量大于j点对该种货品的需求量;式(12)为载重平衡约束,表示当第v辆车从i点到达j点时,离开i点时每类货品的载量加上j点的需求量/供应量等于离开j点时该类货品的载量;式(13)为载重量守恒约束,表示车辆从车场带走的k类货品加上整条路径上所有门店对该类货品的需求/供应之和,等于返回车场时车辆装载k类货品的载量;式(14)～(16)为变量的取值范围约束。

2 两阶段搜索算法

本文提出的两阶段搜索算法(Two-Stage Search Algorithm,TSSA)采用自然数编码,“0”表示仓库,“1,2,…,i,…,n”分别表示不同的门店,单条车辆路径序列为“[0,…,i-1,i,i+1,…,0]”。第1阶段利用最短路径插入规则构建初始解,第2阶段以变邻域搜索算法为框架,通过两个不同适应度函数进行反复迭代搜索;定义了多种形式的更优解,引入多种邻域解的评价方式,接受各种形式的更优解;通过改进的扰动机制,保证算法能够跳出局部最优解的同时,又增强了算法对解空间的搜索能力。

2.1 构建初始解

初始解采用一种近似贪心算法的最短路径插入规则,在满足车容量和时间窗约束的前提下,尽可能多地在一条路径中插入门店。算法步骤如下:

(1)确定未被插入的门店集合A(A=N),生成一条车辆路径R={[0]};

(2)插入距离仓库最远的门店i(i∈A)作为空路径r(r∈R)的第1个点,A=A{i};

(3)将剩余的门店u(u∈A)依次插入路径r(r∈R)的每个位置a(不包括0点)前面,计算每个门店u插入路径r每个位置a的适应度值构成集合B1,集合B2记录了门店u每一次插入的位置,如果插入之后使得路径变成不可行解,则设置该次插入的适应度值为无穷大;

(4)若集合B1中记录所有适应度值均为无穷大,则表示当前所有路径没有可行的插入点,新增一条空路径R=R∪{[0]},当前插入路径r更新为新的空路径,循环步骤(2)～(3),否则执行步骤(5);

(5)B1中最小数值对应的门店即为最佳的插入门店u*,B2中门店u*对应的插入位置即为最佳的插入位置a*;

(6)将门店u*插入到路径r的位置a*,A=A{u*};

(7)循环步骤(3)～(6),直到A=φ时停止循环;

(8)在所有的路径最后插入0。

集合B1中每个元素的适应度值按下式计算:

式中:ra表示路径r的第a个门店;cra,u+cu,ra+1-cra,ra+1表示将u点插入路径r第a和第a+1个门店后该条路径运输成本增加值;frau表示将u点插入当前所有路径的成本增加值与将u点插入一条新的空路径的成本增加值的差值。

2.2 改进变邻域搜索算法

2.2.1 操作算子 常用的操作算子包括relocate、shift、cross、inter-exchange 以 及intro-exchange。Relocate是路径间的单点插入算子,是将一个点插入到其他路径中;shift是指路径内将一个点插入到其他位置;cross 是指交换两条路径的后半段;exchange是指交换两个点,包括路径间interexchange和路径内intro-exchange两种情况。本文除了使用这5个操作算子,还使用了两个扩展操作算子和k-opt算子。

图1 所示为路径间的随机长度片段交换算子λ-inter-exchange,对任意两条路径,取路径上的任意一个点后面的随机长度片段(片段可以只有一个点,不包括0)进行交换。图2所示为路径内的随机片段交换算子λ-intro-exchange,对长度大于等于4的路径,路径内任意两点将路径分割为3个片段,从第2和第3个片段最左边开始取随机长度的片段(片段可以只有一个点,不包括0)进行交换。两个操作算子是对exchange算子的扩展,增加了每一次邻域操作片段的长度,采用每次随机取每条路径中的某个片段进行交换,扩大了邻域解的搜索空间,在多次迭代之后保证了对最优解的搜索能力。

图1 λ-inter-exchange

图2 λ-intro-exchange

由Helsgaun[23-24]提出的k-opt算子是一种对单条路径寻优效果好且效率高的操作算子。k-opt总是在节点数较少的情况下进行交换,其效率可以相当于2.2-opt。k-opt首先使用不同长度的2-opt优化路径长度,在2-opt无法找到更优解时则使用3-opt操作,如果3-opt能找到更优解,则回到2-opt继续寻找;如果3-opt不能找到更优解,就使用4个节点的交换来寻找更优解,以此循环。

2.2.2 适应度函数 由于服务时间窗的导入,使得从仓库/门店离开到达下一个门店/仓库的时间以及在门店开始装卸的作业时间存在差异。根据时间窗早到等待的规则,一是车辆实际开始服务门店的时间,另一是车辆实际达到门店的时间,两者的计算方式为:

式(18)同式(9),车辆开始服务门店的时间须满足每个门店的时间窗限制。式(19)中,车辆v路径上第a个门店的实际到达时间等于路径上第a-1个门店的服务开始时间加上在该门店的服务时间再加上从该门店到第a个门店的运输时间,

式(20)计算了上述两个时间的累积差值。最小化时间差累计和能引导算法朝着时间窗更加紧凑的方向进行搜索。

下式给出了算法的第1 个适应度函数,依次为路径的运输成本、仓库的处理成本、时间窗的违背值、载重量的违背值及邻域解的紧时间窗引导函数T,计算中后面三部分通过权重系数进行调节,

由式(21)可知,对于邻域搜索出现的可行解,时间窗和载重量的违背值都等于0,并不影响转运成本;而不可行解中总有其中之一或两个违背值都不等于0。因此,在最小化的搜索过程中,算法总是朝着可行解的方向进行搜索。引导函数T的加入,使得适应度函数值总是大于实际的转运成本,故引导函数T不能一直存在于适应度函数中。

本算法设计通过适应度函数f1在经过连续Z次迭代后,就不再用引导函数进行评价,故算法提出了第2个适应度函数f2,即

2.2.3 邻域解评价及优先保留算子 邻域搜索中出现的邻域解分为可行解和不可行解两种。其中:①当邻域解的解序列中车辆的实时载量超过车辆的最大载量;②当邻域解的解序列中存在某个点的时间窗违背了该点前后两点的时间窗。只要这两种情况出现一种,就构成了不可行解。可行解中出现车辆数减少的情况,这种解称为优先可行解,记为sol1;可行解中车辆数未减少,但是总的距离缩短,这种解称为距离更优的可行解,记为sol2。同样,不可行解出现车辆数减少时,称为车辆数更优的不可行解,记为sol3;不可行解中车辆数未减少,总的距离缩短,称为距离更优的不可行解,记为sol4。sol1与sol2是算法前进的方向,并且车辆数减少的sol1是优先考虑接受的,但是产生sol3与sol4时也不能被即刻淘汰,也许基于当前不可行解进行搜索会产生sol1与sol2。因此,算法给出了两种不同情况下的优先保留算子:

(1)当产生sol1 时,适应度函数为:f1=f1/ξ1,f2=f2/ξ1。

(2)当产生sol3 时,适应度函数为:f1=f1/ξ2,f2=f2/ξ2。

其中,ξ1～ξ2 ＞1,∀ξ1,ξ2 ∈N*。在寻找最小值的邻域搜索中,这两种情况会被优先保留下来,两者的区别在于,前者的邻域解在向下传递邻域解时会将空路径删除,而后者不能删除空路径,除非在接下来的搜索中出现sol1或者sol2。

2.2.4 扰动机制 传统VNS算法在每次邻域搜索结束后,将邻域解传入新的操作算子中,通过随机指定操作节点的方式,在有限的时间内(通常时间小于1 s),将搜索过程中目标值最小的可行邻域解传入下一个操作算子进行搜索,如果在有限的时间内未找到可行的邻域解,则将原邻域解传入下一个操作算子进行搜索。可见,传统的邻域搜索在每次邻域搜索完成后都对邻域解进行扰动,并且扰动采用随机的方式,保证了VNS算法能够跳出局部最优解。然而,由于搜索过程中扰动频率太高,使得VNS算法对当前邻域解的搜索不够彻底,并且扰动采用随机的方式,扰动后的邻域解已大幅偏离原邻域解,再一次增大了对当前邻域解搜索不够彻底的可能。

在最小化的寻优算法中,操作算子在遍历当前所有可能的邻域解后,找到一个比当前邻域解更优的邻域解替换当前邻域解并向下传递,如果搜索出现的邻域解没有当前的邻域解好,则不做替换。针对这一特点,本算法设计在连续迭代Y次后都未找到优化解,则将当前邻域解的目标值设置为无穷大,在下一次操作算子搜索过程中,目标值为无穷大的邻域解很容易就会被新的邻域解替换。通过不断地寻优替换,在操作算子一次完整的搜索之后,又能重新得到一个质量较好的局部最优解。

2.3 算法流程

如图3所示,首先产生一个可行的初始解,初始化算法参数并将neibor_obj设置为无穷大,然后将初始解传入第二阶段搜索算法。邻域搜索过程用8个操作算子中neibor_num 对应的算子进行搜索。对搜索出现的每一个邻域解按照2.2.2和2.2.3节进行计算和评价。如果获得了优化解(sol1 或sol2),则将当前邻域解的目标值更新给neibor_obj,同时将当前优化解传递给全局优化解,然后迭代次数iteration加1,操作算子序号neibor_num 加1,count_break重置为0;如果本次操作算子搜索没有找到优化解,而搜索到更优的不可行解(sol3 或sol4)时,同样将当前邻域解的适应度函数值更新给neibor_obj,然 后iteration 加1,neibor_num 加1,count_break也加1;当上述4种邻域解的情况均未出现,即说明本次邻域搜索未找到任何更新解,只需将iteration加1,neibor_num 加1,count_break加1。在迭代过程中,邻域搜索算子循环使用。判断邻域解是否需要扰动,如果count_break是参数Y的整数倍,则重置邻域解目标值为无穷大。直到迭代次数达到最大迭代次数或连续迭代结果未改善次数超过给定阈值,则终止算法,输出结果。本文经过计算后将最大迭代次数设置为2 000,连续迭代结果未改善次数阈值设置为400。

图3 算法流程

3 算法验证及结果分析

TSSA 算法由Python 语言实现,运行环境为Intel(R)Core i3-8100@3.60 GHz CPU 及12 GB内存的个人电脑,模型求解采用Gurobi优化器,设置每个算例的最大计算时间为3 600 s。

3.1 VRPTW 计算分析

本算法设计针对M-PDPTWT,M-PDPTWT作为VRPTW 变体,在货品的取送方式上更为复杂,但因其具有同样的时间窗属性,使得针对MPDPTWT 的两阶段搜索算法在修改适应度函数和取送货约束后,算法同样适用于VRPTW。

算法中的两个适应度函数修改如下式所示,两个适应度函数中依次为运输成本、时间窗违背值、载重量违背值以及引导函数,

式中,Tlv为车辆v总的载重量。

邻域解中可能会出现Tlv＞Q的情况,即

式中,pi为每个顾客/门店需求的货物量。

2.2.4 节中的参数Z用以界定两个适应度函数,其取值会影响算法的收敛速度和可行解的求解质量。表1给出了不同取值的Z所求得的算例结果,其中3个数据依次表示车辆数/距离/求解时间。

当Z=0时,表示算法没有用适应度函数f1(或f3)进行评价。由表1 结果可见,当Z=0 时,算法将以最快的速度收敛,但是不论是在车辆数的优化还是距离的优化上,求解结果都不如加入参数Z的求解结果好,证明了本算法设计先用f1(或f3)进行评价再用f2(或f4)进行评价的搜索方式是有效的。当Z=10时,算法的综合性能保持得最好,在求解质量和求解速度上都能找到一个平衡,故算法采用Z=10作为最优参数。

表1 参数Z 的调试(100个门店部分数据)

附表1给出了TSSA 算法求解100 个门店的算例结果,同时给出了算例的最优解和近年发表的一些启发式算法的求解结果。计算结果表明,在大多数算例中,TSSA 算法都能找到最优的求解路径。

附表1 Solomon算例结果对比(100个顾客)

3.2 M-PDPTWT计算分析

采用某女鞋企业的销售数据进行算例设计,公司提供了近两年某城市包括仓库中心在内的31个门店的销售数据,数据包括4个季节中不同品类的女鞋样式18种。算例按照每个月每种鞋的平均销量作为该种鞋下一周的预期销售量,低于平均值就需要补货,高于平均值就需要运走货品。同时,为了应对未来城市规模的扩张,将门店的数量增加到50个,增加的门店位置在已有门店组成的地理位置网络上随机取点,新增门店的销售数据也通过随机的方式产生。因此,采用门店数S={5,12,18,30,50},货品数目C={2,3,6,10}组成18组,每组4个共计72个算例进行算例测试。其中,算例时间窗为8:00～20:00,共计720 min,每个门店的时间窗按照间隔大于30 min(服务时间为30 min)进行随机生成。

表2中S为除去中心仓库外门店的数量,C为货品的数量,NO 为算例序号,每组4个。每组算例在两阶段搜索算法中计算10次,Best-TSSA 为10次计算中最好的值,Aver-TSSA 为10次计算的平均值,Gurobi为通过Gurobi优化器求得的精确解,每一栏都包括车辆数(Veh)、距离(Dis)和计算时间(Time)。最后一栏为两阶段搜索算法的最好结果与Gurobi结果的差值。在用Gurobi优化器求解模型时,对两个目标函数设置不同的优先级顺序,具体的设计参数为:priorityobj1=0.9,priorityobj2=0.1。

表2 M-PDPTWT算例计算结果

续表2

由表2 计算结果可以看出,Gurobi在规定时间3 600 s内计算出精确解的22个算例中,TSSA算法都能求解到相同的结果,证明了TSSA 算法的准确性和数学模型的正确性。另外,有19个算例Gurobi能够找到与TSSA 算法相同的车辆数,但是在规定的时间内无法求解到距离的最优值;在剩下的31 个算例中,Gurobi无法找到与TSSA算法相同的车辆数和相同的距离值。在时间方面,Gurobi除了在小规模算例上(5个门店)花费的时间优于TSSA 算法,其余算例的计算时间都长于TSSA 算法。TSSA 算法的平均计算时间随着算例规模的增长成比例增长,均保持在合理的计算时间范围内。

4 结语

本文针对某快时尚品调配货品的实际问题,结合问题特性,将其归类为带时间窗和调货特性的多品类取送货问题(M-PDPTWT)。以最小化车辆数为第一优化目标,最小化配送成本为第二优化目标,建立了完整的混合整数规划模型。针对局部邻域搜索算法很容易陷入局部最优解这一缺点,提出了利用不同的适应度函数和不同的邻域解评价对邻域解进行定量和定性的处理方式,在利于减少车辆数的第2阶段和利于减少总距离的第3阶段反复交替迭代中,实现了对最优解的最大可能的寻找。将两阶段搜索算法用于求解Solomon算例的结果表明,两阶段搜索算法能求解到大部分与最优解相同的最优解。最后,根据企业提供的销售数据设计了用于求解M-PDPTWT 问题的测试算例,通过两阶段搜索算法与Gurobi求解的精确解对比发现,两阶段搜索算法在小规模算例上都能找到与精确解相同的最优解,在无法求解出精确解的较大规模算例上,两阶段搜索算法也能在较短时间求解出满意的结果。因此,本文设计的两阶段搜索算法能够满足快时尚品企业对城市各个网点的调配需求,为企业提供满足条件且运输成本最低的车辆行驶路径。