考虑时空相关随机行驶时间的车辆路径问题模型与算法

张冬青， 郭钊侠， 张殷杰

(四川大学商学院， 成都 610065)

1 引 言

城市化加快了电子商务产业的发展，也催生了以及时、快捷、高效为目标的电商快递、生鲜配送、冷链物流等行业. 在城市物流体系中，最后一公里配送是对接顾客的重要环节，其质量直接影响顾客的体验和满意度. 因此，如何制定最优的车辆配送方案和行驶路径，高效地完成配送任务，一直是学术界研究的重要课题. 学术界通常将最后一公里配送问题建模为车辆路径问题[1].

城市物流配送过程中存在着各种随机因素. 根据所考虑的随机元素的不同，随机车辆路径问题(Stochastic Vehicle Routing Problem, SVRP)主要可以分为考虑随机顾客需求的车辆路径问题[2]、考虑随机顾客的车辆路径问题[3]和考虑随机时间(包含随机的顾客服务时间、随机的车辆行驶时间)的车辆路径问题[4-5].

根据所考虑的道路网络上车辆行驶时间是否具有时间依赖性，即车辆通过一条路段的时间是否会随着从该路段端点出发时间的不同而发生改变，考虑随机车辆行驶时间的SVRP又可以分为非时间依赖的SVRP[5]和时间依赖的SVRP[1]. 绝大多数考虑随机车辆行驶时间的SVRP文献都假设各路段的车辆行驶时间是服从于某种概率分布的独立随机变量，如正态分布、gamma分布和log-gamma分布等. 但是这些文献并没有考虑真实道路网络中车辆行驶时间所存在的相关性.

基于真实的车辆行驶时间或速度数据，有学者[6-8]发现交通流在道路网络上的演变会导致各路段的车辆行驶时间(或速度)之间存在时间和空间上的相关性(简称时空相关性)，即某一路段某一时段与另一时段另一路段的车辆行驶时间之间存在相关性. 虽然，道路网络车辆行驶时间的时空相关性已经在随机最短路径问题[9]和交通预测[10]等领域中被广泛考虑和研究，并且，Gendreau等[4]也指出考虑随机量之间的相关性是SVRP未来重要的研究方向，然而，考虑随机相关行驶时间的SVRP研究在国内还未见报道. 在国外，现有SVRP研究中考虑相关性的文献也非常有限. Rostami等[11]研究了考虑随机相关行驶时间的SVRP，以最大化行驶时间的可靠性为目标，提出了一种分支-定价-分割算法进行求解. Rajabi-Bahaabadi等[12]进一步考虑了软时间窗约束，并以最小化车辆的行驶成本和违反时间窗的惩罚的总和为目标，提出一种最大-最小蚁群系统算法来求解. 然而，这两篇文献都没有考虑真实的道路网络和行驶时间的时间依赖性以及在时间维度上的相关性.Avraham和Raviv[13]进一步考虑了时间依赖的行驶时间和时空相关性，并研究了单个车辆的随机相关SVRP，但是他们同样没有考虑真实的道路网络.

基于真实道路网络，本文考虑车辆行驶时间存在的时空相关性和时间依赖性，对相应的SVRP进行了研究. 该问题旨在将顾客所需要的产品分配到合适的车辆上，并安排每一辆车的行驶路线，来达到最小化司机加班费用与车辆CO2排放成本的总和的目标.

一方面，SVRP集成了整数规划和随机规划的性质，是典型的NP-hard类组合优化问题. 因此，精确算法往往无法在合理的计算时间内进行求解，而启发式方法和智能优化算法是广泛应用于求解规模较大路径优化和分配调度等问题的主流方法. 另一方面，随机行驶时间之间的时空相关性会导致高维相关随机量. 例如，一个包含1 000条路段和10个时间周期的道路网络共有10 000个相关随机量. 当处理如此高维度的相关随机量时，传统的蒙特卡洛仿真方法往往需要大量的场景才能反应出随机变量之间的相关性[4]，进而导致过长的求解时间. 而基于copula的情景生成(Scenario Generation, SG)技术[14]已经被证明能以很少数量的场景有效地反映高维随机变量之间的相关性[15]. 因此，本文旨在结合智能优化算法良好的寻优能力和SG技术处理高维相关随机变量的优异性能，提出一类适用于不同规模SVRP的智能随机优化方法，解决考虑时空相关随机行驶时间的SVRP.

2.1 问题描述

在目标函数中，本文综合考虑了以司机加班成本为代表的经济成本和以二氧化碳排放成本为代表的环境成本. 针对经济成本，物流企业的运作中，除了人力成本，还存在油耗成本、维修成本、折旧成本等. 然而，由于汽车的油耗量与二氧化碳排放量近似成正比[16]，本文优化二氧化碳排放成本可以间接地降低油耗成本. 其次，车辆的维修成本和折旧成本受多种因素影响，如货物的类别和司机的驾驶技术等，而难以在数学模型中准确地进行定量化表达. 因此，本文选取人力成本(司机加班成本)作为经济成本的代表. 另一方面，针对环境成本，虽然汽车尾气中包含二氧化碳、碳氢化合物、一氧化碳等多种气体，但与二氧化碳排放量(对于轻型车，与油耗量的比值约为2.3 kg/L)相比，尾气中的碳氢化合物和一氧化碳等气体的排放水平很低(与油耗量的比值约0.01 kg/L)[16]. 并且，由于碳交易市场的建立，二氧化碳排放量可以直接转化为企业的一项经济成本来进行优化. 因此，本文选取二氧化碳排放成本来代表环境成本.

本文利用运输排放和能源消耗计算方法(MEET)[17]来计算车辆行驶过程中CO2的排放量. 对于以速度v(单位： km/h)行驶的7.5～16 T的卡车，式(1)计算了其行驶单位距离的CO2排放量e(单位： g/km).

e=871-16v+0.143v2+32031/v2

(1)

2.2 数学模型

所研究的SVRP的两阶段随机规划模型如下式.

(2)

s. t.

(3)

(4)

(5)

(6)

(7)

(8)

2≤|S|≤n-1, k=1,...,K

(9)

xijk,yik∈{0,1},∀i,j=0,...,n,i≠j,

k=1,...,V

(10)

对于每一个场景ξm：

(11)

[Lp,Up], m=1,...,M

(12)

(13)

sj, ∀j=1,...,c,m=1,...,M

(14)

∀j=c+1,...,n, (i′,j′)

∈E, m=1,...,M

(15)

(16)

∀k=1,...,V, m=1,...,M

(17)

∀m=1,...,M

(18)

m=1,...,M, (i′,j′)∈E

(19)

3 智能随机优化方法

本文旨在针对考虑时空相关随机行驶时间的SVRP提出一个高效的算法求解框架. 为了适用于现实中问题规模较大的SVRP，本文利用智能优化算法进行解的寻优，并创新性地将其与SG技术相结合，提出一种智能随机优化方法. 该方法首先利用适用于车辆路径问题的智能优化算法在解空间中搜索候选解；然后，利用SG技术产生的代表时空相关随机行驶速度的场景集合对候选解进行评价. 混合粒子群优化(Hybrid Particle Swarm Optimization, HPSO)算法具有局部优化解与全局优化解相互结合的寻优机制，能够成功地求解确定环境中的车辆路径问题等离散优化问题. 因此，本文设计了一种HPSO算法来进行解的寻优. 在HPSO算法中，我们设计了一种贪婪的启发式方法来得到真实道路网络中两个顾客点之间的具体行驶路径；引入了可变邻域下降算法来对粒子群优化算法搜寻到的解进行局部提升；并设计退火策略来更新自身最优解以避免陷入局部最优.

图1 基于HPSO的智能随机优化方法流程图Fig.1 Flow chart of the HPSO-based intelligent stochastic optimization method

图1展示了本文提出的基于HPSO的智能随机优化方法的流程. 其中，我们首先利用SG技术产生能够代表真实道路网络中时空相关随机速度的M个随机场景；然后，利用HPSO算法搜寻解空间，找到候选解；再基于SG技术产生的随机场景计算并返回给HPSO算法该候选解的适应度值，以帮助HPSO算法继续搜寻解空间找到最优解. 算法的主要部分介绍如下.

3.1 随机速度场景的产生

在随机优化的框架下，首先我们使用基于copula的SG技术[16]产生能代表时空相关的随机行驶速度的随机场景集合. 在给定相同输入数据的情况下，每次运行此SG技术会产生相同的随机场景集. 因此，当所研究的时空道路网络不发生改变时，随机场景只需要产生一次. 利用SG技术产生随机场景集的具体操作步骤如下.

Step1： 输入所有随机速度变量的边际分布函数；

Step2： 输入能代表随机速度变量之间时空相关性的相关系数矩阵；

Step3： 使用基于copula的SG技术产生随机场景集合ξ={ξ1,ξ2,...,ξM}，其中每个场景的发生概率为1/M.

3.2 粒子的编码与解码

本文采用没有分隔符的所有顾客节点的序列来编码所有粒子. 然后，利用传统的随机分割算子进行解码，即将粒子随机分割为与车辆数相等的V段. 如图2所示，需要配送的10个顾客节点被依次划分到3辆车上进行配送.

图2 粒子表示方式Fig.2 Representation of a particle

3.3 种群初始化

为了更好地引导解的搜寻方向，我们首先利用最近邻点启发式方法得到一个初始粒子. 在这个启发式方法中，车辆总会选择距离当前节点最近的顾客点作为下一个配送目标，若配送下一个顾客点的货物会导致车辆超载，则安排下一辆车从车场出发以下一个顾客点为第一个配送点继续进行配送，直至所有顾客点的货物全部配送完毕. 然后，种群中剩余的其他粒子均随机产生.

3.4 粒子的评价

为了增强算法搜索的多样性，本文允许非可行解以一定的惩罚费用进入寻优过程，并利用式(21)对候选解进行评价. 其中，F为SVRP的原始目标函数(式(2))的值；cl为对超载货物的惩罚系数；OL为当前粒子超载的货物量；t为算法迭代的次数.

Z=F+cl·OL·t

(21)

基于SG技术产生的随机场景集ξ，分别对每一辆车的配送路线进行评价. 具体过程如下.

Step1： 令m=1；

Step2： 若m≤M，重复step3和step4，否则转到step5；

Step3： 在场景ξm中，对所有相邻的顾客点对(i,j)，利用3.2节所述方法计算顾客点i到顾客点j之间的行驶路径，并根据式(12)计算hξm的值；

Step4： 置m←m+1；

Step5： 根据式(2)计算F的值，然后根据式(21)计算Z的值.

由于在随机环境中对一个粒子进行评价需要计算其在M个场景下的适应度值，所以随机问题中解的评价时间会比其对应确定问题中解的评价时间长很多. 为了节约计算时间，在用于全局搜索的粒子群优化算法中，我们使用M个场景对解进行评价，而在用于局部提升的可变邻域下降算法中，我们仅使用一个确定场景. 在此确定场景中，各路段各时间周期的车辆行驶速度设置为相应随机速度变量的均值.

3.5 粒子位置的更新

用于全局搜索的粒子群优化算法根据一定的执行概率利用三个作用于粒子的随机变化操作，即突变、与自身最优解交叉和与种群最优解交叉，来更新粒子的位置. 具体使用的算子分别是swapmutation算子、erXover算子和erXover算子. swapmutation算子是经典的变异算子，即随机选择粒子中的两个顾客点交换其位置；erXover算子基于所选父辈粒子中每一个节点的邻居节点数量，通过循环地选出当前节点的相邻节点中具有最少剩余邻居节点数量的节点，来进行边的重组以得到交叉后的粒子位置.

用于局部提升的可变邻域下降算法依次使用三种路径间变化算子crossover、swap(2, 2)和shift(1, 0)，以及两种路径内变化算子2-opt和relocate.对于以一定概率从种群中随机选择出的一个粒子，首先对其执行crossover算子，如果找到更好的解就更新当前粒子位置，并依次执行2-opt算子和relocate算子；如果没有找到更好的解就按此规则依次执行路径间变化算子swap(2, 2)算子和shift(1, 0)算子. 以此类推，直至执行完shift(1, 0)算子且没有找到更好的解，则结束该可变邻域下降算法，返回粒子的当前位置.

其中，路径间变化算子crossover是随机选择分段点，将两个父辈粒子都分割为前后两部分，再将父辈粒子1的前一部分和父辈粒子2的后一部分组合成为改进的粒子，并将剩下的两部分组合成为另一个改进的粒子.swap(2, 2)算子是从两个父辈粒子中分别选择两个相连的节点片段，并交换这两个所选片段的位置以组成改进的粒子.shift(1, 0)算子是将一个父辈粒子中的一个顾客节点插入另一个父辈粒子中，以得到改进的粒子. 路径内变化算子2-opt算子先将当前粒子的两条边断开，并组成另外两条新的边来代替断开的边，得到改进的粒子.relocate算子是将一个顾客节点从粒子中移除，再重新插入到当前粒子的其他相邻节点之间以改进粒子性能.

3.6 更新自身最优解和种群最优解

每次迭代后，检查是否达到了给定算法的最大迭代次数或达到种群最优解连续未更新的最大迭代次数. 如果满足任一条件，算法运行终止，返回种群最优解作为此智能随机优化方法找到的最佳解.

4 数值实验

4.1 算例描述

为了验证设计的智能随机优化方法的有效性，本文基于Augerat等[18]提出的带有装载限制的车辆路径问题算例A-n36-k5，构建了车场和顾客数量的总和分别为19、27和36的三个算例集. 这3个算例集的顾客需求分别为A-n36-k5中前19、27和36个顾客的需求，其中第一个节点为需求等于0的车场. 车场中可利用的车辆数分别为3、4、5，且每辆车的最大载重与算例A-n36-k5中的设置一致. 所有测试算例均定义在包含327个节点和1 207条边的北京市区道路网络上. 在每一个算例集中，利用Huang等[1]使用的方法得到顾客随机分布(R)、集中分布(C)、一部分随机分布另一部分集中分布(RC)的3个算例.图3展示了在简化的北京道路网络中算例R36的顾客分布情况.

图3 算例R36的顾客分布

司机的单位时间加班费用c1设置为25元/h. 根据报告[19]，单位CO2排放成本c2设置为248.15元/吨. 以早班为例，设定车辆在早上8∶30离开车场，中午12∶00返回，则司机的工作时长WT为3.5 h. 每位顾客的服务时间设为20 min. 道路交通情况变化的时间周期设置为30 min，则整个计划期有8个时间周期(如表1所示). 因此，本文需要处理含有9 656(1 207×8)个相关随机量的时空网络.

与文献[15]一致，本文设定每条路段每个时间周期的车辆行驶速度服从有界beta分布.其中,参数α和β等于Dmax/(l·Kp)，Dmax表示道路网络中距离最长的路段长度，l为当前路段的长度，Kp为当前时间周期道路拥堵程度的常数，且Kp越小代表道路越拥堵. 每一个时间周期道路网络的Kp值、最大、最小和平均的车辆行驶速度如表1所示. 为了表示道路网络行驶时间的时空相关性，本文设定两个连续时间周期之间的时间相关系数为0.4；两个相邻路段之间的空间相关系数为0.4；相邻路段在相邻时间周期的时空相关系数为其时间相关系数与空间相关系数的乘积；其他情况下的相关系数为0. 利用样本内外测试方法，发现对于所构建的测试算例，10个随机场景即可保证对解评价误差在2%以内，因此本文设定随机场景数量M为10.

通过多次实验比较，设定智能随机优化方法的参数如下：种群中的粒子总数为20，超载货物的惩罚系数cl为0.5，用于更新粒子位置的swapmutation算子的执行概率为0.2，erXover算子的执行概率为0.6，执行可变邻域下降算法时从种群中选择粒子的概率为0.3，更新自身最优解时的初始温度T(0)为30，程序的最大迭代次数为500，种群最优解连续未更新的最大迭代次数为100. 所有程序使用MATLAB R2016b编写，程序运行环境为Inter Core i5-6200U CPU @ 2.30 GHz和12 GB RAM的笔记本电脑.

表1 各时间周期速度分布的设定

4.2 实验结果

4.2.1 算法比较 为了检验设计的HPSO算法的求解性能，本文将其与单纯的粒子群优化算法和混合遗传算法进行对比. 其中，提出的HPSO算法中若不使用用于局部提升的可变邻域下降算法，即为粒子群优化算法；传统的遗传算法结合上述可变邻域下降算法，即构成混合遗传算法.其中所用到的选择算子为锦标赛选择算子，交叉和变异算子分别为erXover算子和swapmutation算子. 所有算法运行10次，结果的平均值如表2所示. 由于对于同一个时空网络，利用SG技术产生的随机场景可以重复使用，表格中的程序运行时间不包含产生随机场景的时间.

表2中的结果表明,混合粒子群优化算法在求解质量和计算效率上均具有一定的优势. 虽然，粒子群优化算法所用时间很短，但由于没有对解进行局部提升，其求解性能明显差于另外两种混合算法. 混合粒子群优化算法和混合遗传算法的求解时间差异不大，但是很明显，对于所有算例，混合粒子群算法的求解性能总是更优.

4.2.2 时空相关性分析 为了进一步分析在SVRP中考虑道路网络车辆行驶时间的相关性的意义，本文以算例R27为例比较了不同的时间和空间相关系数组合下的求解结果. 取10次运行的最佳解，结果如表3所示. 在相关系数不同的随机场景中，同一个解的评价值会存在差异；为了研究不同相关系数组合对解的影响，本文利用时间和空间相关系数均为0.4的随机场景重新计算了所有解的目标值，并将其作为比较基准值.

表3中的结果表明，考虑与不考虑时空相关性所得到的解会存在一定差异；考虑不同的时间和空间相关系数组合往往会得到不同的解，进而导致不同的总成本；并且，若考虑的时间和空间相关系数与实际情况不符，往往会得到在实际中表现更差的解. 因此，准确地考虑真实道路网络中车辆行驶时间所存在的时空相关性才能得到适用于实际交通环境的最优路径决策方案.

本文进一步分析了不同空间和时间相关系数分别对各项成本的影响. 图4展示了时间相关系数(或空间相关系数)为0.4时，不同空间相关系数(或时间相关系数)下各项成本的数值. 由图4可以看出，总成本的波动主要来源于加班费用的波动，这说明相比于CO2排放成本，司机加班时间(费用)受不同空间和时间相关系数的影响较大. 因此，考虑道路网络车辆行驶时间的时空相关性对企业减少司机的加班时间和提高物流配送效率有重要的现实意义.

表2 不同算法结果比较

表3 不同时空相关系数组合下的总成本(单位：元)

本文以真实的城市道路网络为背景，以最小化司机加班费用和车辆CO2排放成本的总和为目标，研究了考虑时空相关的随机行驶时间的车辆路径问题. 首先，建立了该问题的两阶段随机规划模型. 然后，结合智能优化和随机优化技术，提出了一种智能随机优化方法. 该方法利用一种结合了可变邻域下降算法的HPSO算法进行解的寻优，并利用基于copula的SG技术产生随机场景来对候选解进行评价.

基于北京市区道路网络，本文进行了一系列的数值实验. 实验结果验证了所提出的智能随机优化方法的有效性. 同时，结果表明对于所研究的问题，提出的混合粒子群优化算法的寻优性能优于混合遗传算法以及传统的粒子群优化算法；考虑真实道路网络中车辆行驶时间之间的时空相关性能够得到更加符合实际交通环境需求的车辆路径方案，从而帮助企业在综合考虑经济和环境成本的情况下做出更高效的车辆路径决策.

本文的数值实验是基于假设的时空相关车辆行驶速度数据，未来的研究将通过获取真实的交通数据来进行实证研究. 另一方面，近年来，一些精确算法被提出求解具有特定特征的VRP. 未来的研究可以在本文提出的算法框架下，利用这些精确算法替换本文用于寻优的HPSO算法，来求解考虑时空相关的随机行驶时间的SVRP，并对比研究其与本文提出的智能随机优化方法的性能差异以及各自的适用场景. 未来的研究还可以集中在考虑更多的实际应用场景，比如多个优化目标、异质型的车辆、顾客服务时间窗的限制等.