袁晓建, 张岐山,吴 伶,江义火
(1. 福州大学经济与管理学院, 福建 福州 350108; 2. 福州外语外贸学院理工学院,福建 福州 350202; 3. 福州大学数学与计算机科学学院, 福建 福州 350108)
带时间窗和同时送取货车辆路径优化问题(vehicle routing problem with simultaneous pickup-delivery and time windows,VRPSPDTW)作为车辆路径优化问题(vehicle routing problem,VRP)的扩展,该问题整合了企业正向物流和逆向物流,具有较强的应用价值,自提出以来得到了学者们的持续关注[1].
由于VRPSPDTW问题增加了同时送取货和时间窗两个约束,当配送客户规模的增大,求解的规模成指数级增长,传统的精确求解算法难以在可接受的时间范围内求得最优解. 随着启发式算法的提出,学者们开始尝试结合该算法求解VRPSPDTW问题. 文献[2]使用共同进化的遗传算法(genetic algorithm,GA)求解该问题,并在Solomon数据集基础上设计了65个数据集. 文献[3]使用遗传算法和禁忌搜索算法对家庭医疗物流问题进行求解. 文献[4]设计了并行模拟退火算法(parallel simulated annealing,p-SA)求解该问题,并在超大客户规模的测试集上验证了其算法的有效性. 文献[5]又提出了基于离散布谷鸟算法(discrete cuckoo algorithm,DCA)的仿生算法对该问题进行求解,并将部分算例的最优解进行更新.
文献[6-7]提出了量子进化算法,该算法将量子计算理论与进化算法进行巧妙结合,是一种概率搜索优化算法. 由于量子进化算法具有种群规模小、 收敛速度快和全局搜索能力强等特点,迅速得到了国内外研究人员的关注[8-11]. 近年来,也有部分学者尝试使用量子进化算法解决VRP相关问题[10, 12-15],但是利用量子算法求解VRP问题存在量子进化方向引导、 算法域到问题域映射等难题,国内外鲜见有学者进行深入研究.
本研究基于车辆装载率和客户-中心点方向的权重距离公式,构建量子元胞体完成算法域到问题域的映射,设计双层的改进量子算法,最后在VRPSPDTW问题通用测试集上对算法进行有效性检验.
该问题的图论描述如下:G(N,A)表示一个有向图,N={0, 1, …,n}表示有向图的节点集合,其中节点0表示中心车场,也是货物配送中心,是车辆出发和返回地点,要求所有的车辆从中心车场出发,最后回到中心车场; 节点1, …,n表示客户节点,每个客户i∈N同时有货物配送需求量Di>0和货物取货需求量Pi. 每两个节点之间有两条弧线相连,A={(i,j)∣i,j∈N;i≠j}表示所有弧的集合;cij{i,j∈N;i≠j}表示节点之间的配送成本,这里用两点间的欧式距离表示,且cij=cji. 假设配送中心有m辆同种类型的配送车辆,每辆车的核定载重为Q.
配送中心和顾客均有时间窗限制,配送中心的时间窗为[a0,b0],即车辆只能在a0之后(包含a0)离开,且只能在b0之前(包含b0)返回. 对于编号为i顾客也有预先设定的时间窗[ai,bi],车辆服务该顾客的时间必须在此区间内,即ai车辆提供配送服务的时间不能早于ai且不能晚于bi,同时设车辆为顾客i提供服务的时长为ti. 为研究方便,假设车辆的行驶速度为1,则车辆从客户i到客户j的行驶时间等于他们之间的欧氏距离,即tij=dij.
每辆车的启动成本设置为cd,考虑到市场实际的资源配置,本研究的第一优化目标为车辆数最少; 车辆数相等的条件下,第二优化目标为配送距离最短.
VRPSPDTW问题的目标函数是以最小的成本(最少配送车辆和最短配送距离)完成配送任务,并满足以下假设:① 所有顾客的配送需求和时间窗要求均要满足,且需求不可拆分,即每个顾客只能由同一辆车同时完成取货和送货; ② 每辆车从中心仓库出发并返回中心仓库; ③ 车辆在服务顾客时,先卸货后装货; ④ 车辆的最大载重不允许超过车辆的核定载重.
为方便建模,现将使用到的集合、 参数和变量定义如下:
集合.N为节点集(顾客集和中心仓库),N={0, 1, …,n}表示有向图的节点集合,其中0表示配送中心;N0为顾客集,N0={1, 2, …,n};K为车辆集,K={1, 2, …,m}.
参数.Q为车辆的额定装载量,即最大装载限制;cd为每启动一辆车的成本;ct为单个车辆单位距离的行驶成本;dij为顾客i和j之间的欧式距离,i∈N0,j∈N0,i≠j;Pi为顾客i的取货需求量,i∈N0;Di为顾客i的送货需求量,i∈N0;ai为顾客i允许最早开始服务的时间,i∈N0;bi为顾客i截止服务的时间,i∈N0;ti为车辆服务顾客i所需的时间,i∈N0;λ为权衡分派成本和行驶成本的参数,λ∈[0, 1];M为一个任意大的正数.
决策变量.xijk为车辆k(k∈K)的行驶变量,xijk∈{0, 1}, 如果车辆k从点i∈N行驶到点j∈N,则xijk=1,否则xijk=0;L0k为车辆k(k∈K)离开配送中心时的载货量;Li为车辆服务顾客i∈N0后的载货量;Sik为车辆开始服务顾客i∈N0的时间,如果车辆k没有服务顾客i,则Sik=0.
基于以上变量,参考文献[2]给出该问题MIP模型:
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
ai≤sik≤bi;xijk∈{0, 1} (∀i∈N, ∀j∈N, ∀k∈K)
(9)
其中,式(1)为该问题目标函数; 约束(2)表示每个顾客只能被一辆车服务,即服务不可拆分; 约束(3)避免形成不包含配送中心的回路; 约束(4)保证车辆的初始出发点和最终点均为配送中心; 约束(5)~(7)表示车辆载重约束条件,即车辆在配送全程任何时候的载重均不可超过其核定载重.
量子进化算法(quantum evolution algorithm,QEA)是由量子计算理论与进化算法相结合的算法,是一种概率搜索优化算法,其基本原理是:用量子位编码表示染色体,用量子门作用和量子旋转完成进化搜索[16-17].
量子进化算法的工作机理为:使用量子染色体表征问题解空间中任意解的叠加态,随机观测量子染色体,量子染色体以一定的概率观测到一个以二进制形式表示的解; 算法采用量子旋转门作为进化策略,原理是将需要的结果以概率的形式增加,不需要的结果以概率的形式减弱,保证解向着最优的方向进化.
(10)
其中:θi=S(αi,βi)Δθi,Δθi为量子旋转门旋转角,其大小控制算法的收敛速度;S(αi,βi)控制旋转角的方向,使算法向最优解方向搜素.S(αi,βi)和Δθi的取值可通过查表法[18]获得.
为避免由量子染色体坍塌为二进制解时出现不变解的情况,文献[17]对式(10)中(α′i,β′i)T,作如下处理:
(11)
式中:ε是一个较小的正数(本研究ε取值为0.01).
3.1.1初始解参数设计
由于初始解对算法的收敛速度影响较大,仅仅使用贪心算法又容易减少初始解的多样性. 如图1,车辆配送货物至客户i时,车辆已经接近满载,按照常理,此时车辆应该准备回场,最优的做法是“顺路配送”,即下一个客户尽量是在配送中心的方向上. 图1中,黑色原点表示中心车场,此时车辆在顾客i处,顾客i+1和顾客j+1相比较,明显顾客i+1距离车辆最近,如果单纯使用贪心算法则顾客i+1被加入到配送序列,以此类推,车辆距离中心车场越来越远,直至顾客i+n车辆满载,车辆必须回程,此时车辆距离配送中心比较远,最后车辆不得不长距离的回场. 理想状况是,车辆配送至顾客i时,车辆已经接近满载,即使顾客j与顾客i+1相比,顾客j比较远,但是由于j在车辆回场的方向上,优先选择顾客j作为下一个配送目标,依次类推,直至顾客j+k车辆满载,车辆距离配送中心越来越近. 为了实现此功能,本研究定义带回场权重的权重距离,具体如下:
图1 权重距离示意图Fig.1 Weight distance diagram
图2 向心角示意图Fig.2 Centripetal angle diagram
3.1.2初始解算法设计
初始解生成算法设计如下:① 启动一辆新车,随机选择一个客户作为路径起点; ②i为当前路径最后一个客户,j为下一个配送候选客户, 如果ωi≤ω0,车辆需要行驶的距离取i与j之间的欧式距离dij=dij, 否则,车辆需要行驶的距离取值带有权重的距离dij=wdij; ③ 遍历所有未配送的客户,建立车辆行驶dij由小到大的临时列表; ④ 遍历临时列表,校验时间窗与装载约束. 如果时间窗和车辆装载约束校验均通过,则当前候选客户入到当前路径尾部,算法返回步骤②, 不满足则继续校验临时列表的下一个客户, 直至临时列表中不再有未配送的客户,则当前车辆的配送方案生成完成; ⑤ 重复步骤①,直到所有客户均被添加进配送方案,返回整体配送方案.
3.1.3初始解转化成二进制解
上述步骤生成十进制的整体配送方案后,按照下述3.2编码设计生成二进制的初始解.
传统量子算法在解决VRP问题中,每次迭代本地搜索后,生成的全局最优解转化成对应的二进制解,但是量子的每次进化,无论是启用车辆数减少还是总距离减少,对应的量子的长度和位置均已发生改变,二进制的最优解失去了引导进化的意义. 针对上述问题,这里提出量子元胞体、 特定条件初始化和互换α和β位的解决方法.
3.2.1量子元胞体设计
取量子染色体的α位为观测位,设计长度一定的量子染色体,其中:初始的量子长度=量子胞体的个数×客户个数.
量子元胞体定义为:固定长度的量子染色体作为一个独立的观测单元,观测单元在染色体中的顺序代表客户编号.
元胞体具体形式为(车辆编号|配送顺序),其中该固定长度的前半段观测得到车辆编号的二进制数,后半段观测得到该配送车辆服务的顺序的二进制数. 这里,车辆编号位数的长度LV= log2NV,NV为当前最优解使用的车辆数; 配送顺序位数长度LC=log2n,n表示客户数. 定义量子胞体的长度为L0,L0= LV+ LC; 定义当前最优解下的量子染色体长度为Lq,Lq=nL0. 量子元胞体的具体设计如图3所示.
图3 量子胞体结构图Fig.3 Quantum cell structure
3.2.2量子染色体对齐操作设计
将现有量子染色体观测到的二进制,与由本地搜索到的最优解反向生成的二进制进行逐位对比,一致则对应量子位不变,不一致则交换对应量子位的α位和β位.
3.2.3量子染色体进化机制设计
量子染色体进化机制设计如下:
ⅱ) 由当前最优解(第一代则为初始解)引导方向,根据改进量子算法进行进化,并进行观测、 本地搜索操作;
ⅲ) 当车辆数NV变小,产生新的最优解,返回i),重新生成并初始化量子; 当 NV不变且配送总路程变小,产生新的最优解,进行量子对齐操作,并返回ii);
ⅳ) 有限回溯的全局优化操作;
ⅴ) 迭代达最大设定,输出最优配送方案.
3.3.1初始解多样化
经验证,算法的初始解对收敛速度和解的质量有较大的影响,为扩大搜索空间,这里使用3.1定义的满载率和权重距离生成多种初始解. 如ω0=1时,则满载率和权重距离设置失效,ω0取值不同,得到的初始解不同,灵活实现初始解的多样化.
3.3.2有限回溯操作
以搜索中间状态的所需车辆数划分阶段,当a阶段(a辆车)连续若干代无最优解更新认为陷入局部最优,回到a+1阶段的关键点,重新进行迭代,灵活实现大范围搜索.
为便于衡量算法的有效性,采用文献[2]构造的通用标准测试集. 该测试集包含65个VRPSPDTW测试算例,其中包括9个中小型顾客规模算例和56个大型顾客规模的算例(客户规模为100).
实验环境如下:Intel i5-5200U 2.2 GHz CPU,4 GB内存,操作系统为Windows 7,编程环境为Matlab R2014a.
从10、 25和50个客户规模中随机选取一个算例进行测试,考虑收敛速度和平均最优值,参数设置为:量子种群数目=3,最大迭代次数5 000,θ0=0.15和ω0=0.7对全部的9个算例进行有效性分析.
选取中小规模的9个算例在每个算例上独立运行10次,记录运行结果的最优值进行算法有效性验证,运行结果与GA算法[2]、 p-SA算法[4]和离散布谷鸟(DCS)算法[5]的运算结果进行对比,如表1所示.
表1 四种算法的运算效果对比Tab.1 Comparison of the effects of the four algorithms
为了衡量不同算法在同一个数据集上表现出来的优劣,大多数车辆路径优化算法的分析多采用求解质量和求解速度进行对比. 文献[19]提出了使用非参数检验的方法,使用Friedman检验和Wilcoxon检验分别将不同算法计算结果进行比较排序(Rank),检验算法在置信区间内是否存在显著差异.
在上面的算例实验中,由于求解的是测试集中的中小规模算例,因此求解的车辆数均一样,不具备检验条件,这里选取最短路径进行检验,检验统计量如表2所示.
表2 检验统计量Tab.2 Test statistics
由表2可以看出,在0.01显著性水平条件下,四种算法无显著性差异. 鉴于在Friedman检验中,渐近显著性为0.012,如果在常用的0.05显著性水平条件下,四种算法在Friedman检验中结论为有显著性差异,而算法两两之间显著性不明显.
综上所述,NIQA算法与目前几个著名优化算法不存在显著性差异.
本研究针对带时间窗和同时取送货的车辆路径问题,提出一种新的改进量子算法(NIQA). 通过Wang和Chen测试集的测试以及与现有最优解的比较,表明NIQA是有效的. 虽然所提的量子胞体、 特定条件初始化及互换α和β位等方法解决了量子算法在求解此类问题中的部分困难,给出了些许的思路,但是受限于算法设计能力,在求解大规模客户的问题上,并没有完全展现出量子算法的优越性. 如何设计从量子域到二进制域再到问题域完整的映射,尤其是量子域到二进制的映射,做到尽量少地丢失有效信息是今后的研究方向.