不确定环境下的危险品运输车辆路径优化

王 琴，杨信丰，李 楠，李 平，方晟浩

1.兰州交通大学 交通运输学院，兰州 730070

2.中国铁路兰州局集团有限公司 兰州通信段，兰州 730030

危险品运输车辆路径优化问题（vehicle routing problem for hazardous material，HVRP）自21世纪以来一直是路径优化领域的热点问题。由于危险品的危险特性，一旦发生危险，不仅会对生产、运输企业造成利益损失，更会危及周边人员的人身安全，还会对环境造成污染，因此车辆路径安排合理与否对危险品运输的安全至关重要。

近年来，HVRP的研究已从风险最小的单目标模型趋向于更贴近实际的多目标优化模型。柴获等[1]考虑经过人口密集区域行驶距离最小目标，建立危险品车辆运输路径问题优化模型，并设计基于概率模型的多目标进化算法进行求解。代存杰等[2]采用行程时间和运输风险的随机属性值为优化准则建立了双目标优化模型，设计了两阶段多维标号修正算法求解模型。Liu等[3]在分析路径行程时间和可靠性的基础上，构建了基于行程时间可靠性的多目标调度模型。Bula等[4]以最小化两个冲突目标总路由风险和总运输成本建立优化模型，提出了基于邻域搜索的多目标邻域优势算法和ε约束元启发式算法进行求解。以上研究主要针对静态确定网络环境下的多目标优化，但实际运输过程中往往存在许多动态不确定信息，因此研究动态不确定环境下的危险品运输车辆路径优化问题更具有现实意义。李立等[5]对构建时空矩阵模拟道路信息，建立了动态多目标危险品运输路径优化模型。薛翔等[6]研究不确定需求下的危险品运输车辆问题，建立了随机优化模型。Ghannadpour等[7]将顾客的偏好信息表示为关于服务时间满意度的凸模糊数，基于此提出了带模糊时间窗的多目标动态车辆路径模型。上述文献着力研究单一不确定条件下的危险品运输车辆路径优化问题，尚未考虑实际运输网络中客户需求量、受影响的人口数量以及运输速度等多重不确定因素下的路径选择问题。

目前，关于风险评估模型的研究尚不充分，成果主要集中在如何度量危险品的运输风险上，风险度量方式的合理性直接影响车辆配送方案风险规避的有效性，提出更具鲁棒性和实用性的风险评估模型一直是学者们关注的重点。Kang等[8]提出了一种新的风险价值（valueat-risk，VaR）模型，定义了一种求解多行程危险品运输路径规划问题的广义方法。Men等[9]在风险评估时考虑到人口数量的不确定性，采用IT2-FV型模糊变量（trapezoid interval type-2 fuzzy variables，IT2-FV）表示人口数量建立模型并求解。魏福禄等[10]提出了基于风险价值的危险品运输路径优化设计方法。Zhang等[11]考虑空气污染等构建了综合风险评价模型，将风险因素考虑到维护运营成本中建立了多目标优化模型。葛显龙等[12]提出利用泊松分布生成动态客户建立模型，设计组合优化算法进行求解。张翔等[13]建立了基于模糊综合评价的风险指标体系，利用距离加权法计算危险品运输路线的风险值。Holeczek[14]比较了不同的风险模型，说明了车辆负载动态变化对风险的重要影响。然而现有研究较少考虑载货量对风险评估的影响，实际上，不同载货量的车辆在同一处发生事故所造成的后果截然不同。

鉴于以上分析，本文视动态载货量为影响风险的重要因素，建立了基于动态载货量的风险评估模型，并引入区间数和三角模糊数表征运输过程中的人口密度、行驶速度与运输时间以及客户需求量等不确定变量，以运输总风险、车辆总行程、车辆使用数最小为优化目标，同时考虑时间窗、载货量、事故概率等约束建立了不确定环境下的危险品车辆路径多目标优化模型，并对模型进行清晰化处理，设计了NSGA-II算法（non-dominated sorting genetic algorithm-II，NSGA-II）与LNS算法（large neighborhood search，LNS）相结合的混合NSGAII算法求解模型，并利用solomon算例[15]验证模型和算法的有效性。

1 问题描述

不确定环境下危险品车辆路径优化问题可以描述为：设有向无环道路网络G=(V,E)，V={0,1,2,…,n}为该网络的节点集合0为配送中心，E为网络中边的集合。所有配送车辆须从配送中心出发，在客户时间窗内为每个需求点配送满足需求的危险品，最后返回配送中心。在人口密度、行驶速度与运输时间以及客户需求量不确定的条件下，以运输总风险、车辆总行程及车辆使用数最小为优化目标，求配送车辆的最优路径。

为便于研究，本文作如下假设：

（1）各需求节点需求不可拆分。

（2）所有配送车辆额定载重量相同。

（3）配送中心有足够的库存为客户点配送。

1.1 模型参数说明

定义的参数和变量如表1所示。

表1 参数及变量说明Table 1 Parameter description

1.2 不确定因素分析

1.2.1 人口密度

人口密度是运输风险的评估的重要考量因素，传统风险模型中的人口密度一般为依据经验给定一个确定值。考虑到人口密度的不确定性会影响风险评估的准确性，采用区间数来表示人口密度，即人口密度表示为

1.2.2 行驶速度与运输时间

由于路段的风险、距离、拥堵等情况的不同，危险品运输车辆行驶速度往往无法准确预知，因此将危险品车辆行驶速度以区间数表示为，所以危险品车辆在节点(i,j)之间的运输时间为：

1.2.3 客户需求量

由于配送中心和客户之间的信息不完全共享，客户需求因为时间、空间的差异导致需求结构的不确定性。需求不确定性会造成决策的障碍，致使服务水平降低，为描述不确定需求下的供需关系，引入三角模糊变量，具体形式为͂=(ql,qe,qn)，其中，ql为最小客户需求量，qe为隶属度为1的客户需求量，q n为最大客户需求量。͂的隶属度函数表示为：

1.3 动态风险评估模型

传统风险模型定义事故概率与受影响人数的乘积为运输风险值。危险品运输网络中，设发生危险时受影响区域半径为λ，则以λ为半径范围内的人都可能受到危险的潜在影响，因此，定义受影响范围的面积如图1所示，其面积公式为：

图1 影响面积示意图Fig.1 Schematic diagram of affected area

而发生危险事故的概率与运输距离有关，表示为：

因此运输风险为：

除事故概率和受影响人数以外，载货量的动态变化亦是影响风险评估的关键因素，不同载货量的车辆在同一处发生事故所造成的后果截然不同。本文借鉴文献[17]，认为风险影响后果与载货量存在线性关系，即：

由于客户点i的货物需求量q i为不确定变量，因此车辆k在配送路径(j-1,j)上的动态载货量为：

车辆k途径m路段的风险为：

车辆k的配送路径累计风险为：

运输总风险为：

综上，运输总风险͂表示为：

2 模型构建

2.1 模型建立

综上，建立不确定环境下的危险品运输车辆路径优化多目标模型如下：

式（13）～（15）为目标函数，分别表示总风险最小、总行程最小及车辆使用数最小；式（16）表示某一配送路径风险约束；式（17）表示事故概率阈值；式（18）载货量约束；式（19）表示车辆k到达节点j的时间；式（20）表示客户点时间窗约束；式（21）表示每个需求点仅有一辆车服务；式（22）表示每辆车出发最多访问一个客户；式（23）表示车辆到达需求点i完成配送后须从此需求点出发；式（24）～（25）为决策变量。

2.2 模型清晰化处理

为了将模型中的模糊变量清晰表达，通过区间数排序法及模糊机会约束规划将含有模糊变量的模型转化为确定的数学表达式以便设计算法进行求解。

2.2.1 区间数排序法

对于区间数，将区间数的最大值和最小值代入参数，即可得到包含区间数的参数的最大值和最小值：

因此运输风险为：

同理运输时间为：

其中，ηij为决策者承担运输风险的意愿，规定可接受偏差为：

表示客观变量值与最小变量值之间的偏差，其中φ1max为决策者可接受的最大偏差值。

同理ϑij为决策者承担运输时间的意愿，规定可接受偏差为：

2.2.2 模糊机会约束规划

对采用三角模糊数表示的不确定客户需求量，根据Iwamura等[18]提出的模糊机会约束规划模型，给定置信水平，解决其不确定性，模型形式如下：

式中，x为模糊决策变量；ξ为模糊向量；Pos{·}表示{·}里事件发生的可能性；f(x,ξ)为目标函数；g i(x,ξ)为该模糊规划的约束条件；∂为目标函数的置信水平；βi为约束条件的置信水平，是事先给定的值。

3 算法设计

多目标优化问题的求解主要以精确算法和启发式算法两种解法为主。带时间窗的危险品运输车辆路径优化问题（vehicle routing problem with time windows for hazardous material，HVRPTW）是VRP的衍生问题，属于NP-Hard问题。该类问题的计算量通常伴随客户规模增大呈指数级增长，大规模车辆路径问题只能通过设计有效的启发式算法求解。非支配排序遗传算法（NSGA-II）求解车辆路径问题时在求解精度和计算时间等方面性能优越。本文将全局搜索能力强的非支配排序遗传算法与局部搜索能力突出的大规模邻域算法结合，在遗传算法交叉、变异等寻优过程的基础上引入删除算子和修复算子，设计了求解本文模型的混合NSGA-II算法，该算法主要流程如下：编码方式确定、种群初始化、快速非支配排序、拥挤度计算、选择、交叉、变异、算子的删除及修复过程。

3.1 染色体的编码和解码

染色体编码方式影响算法求解效率，本文采用自然数编码方式来产生染色体，客户数量为n，最大车辆数为k时，染色体长度为n+k-1，染色体为(1,2,…,n,n+1,…,n+k-1)的随机排列，1～n表示客户编号，n+1～(n+k-1)表示不同车辆配送方案的分割点，该种编码方式可以直观高效地反映车辆的配送路径，图2以10个客户(n=(1,2,…,10))和3辆车(k=(1,2,3))的配送方案为例，染色体为（1，7，8，11，2，3，6，5，12，4，9，10）时，11、12为分割点将客户群体划分为三段，每段表示一个车辆的配送方案。个体解码后用0表示配送中心，车辆1从配送中心出发，访问1、7、8后回到该中心，同理可知车辆2、3配送路径。

图2 染色体解码示意图Fig.2 Schematic diagram of chromosome decoding

3.2 种群初始化

遗传算法的收敛性对初始解的质量依赖性较强，为产生质量较高的初始解，本文采用贪婪插入法生成初始解。在满足时间窗约束和载货量约束的前提下，在一条路径上尽可能多插入客户点，直至生成初始种群。该方法旨在尽量减少配送车数，具体步骤如下：

步骤1对于客户i(i=1,2,…,n)，随机生成一个遍历序号为j(j=1,2,…,n)的遍历序列，并初始化车辆数k=1。

步骤2按照遍历序列在车辆1的路径依次插入客户点，直至新插入的客户点不满足时间窗约束或载货量约束。

步骤3不满足约束时保存本条路径上的客户，更新车辆数开辟新路径依次插入剩余客户。

步骤4循环步骤2和步骤3，直至完全遍历客户并产生满足约束条件的初始种群。

3.3 交叉和变异

在交叉操作时，为尽可能保存优秀的父代染色体，随机选取子路径进行交叉，交叉过程中会产生基因冲突，即子代个体中有二次访问的客户点或者遗漏客户点，子路径交叉结束和删除重复的客户点，依次补上遗漏客户点即可。交叉过程如图3所示。

图3 交叉操作图Fig.3 Cross operation diagram

变异过程采用逆转变异法，即在染色体上随机选取两个基因片段，将片段上的个体顺序逆转产生新的个体。

3.4 改进策略

为避免算法早熟，在遗传算法寻优操作的基础上引入邻域搜索的思想改进算法，改进策略根据本文模型设计了相应的删除、修复算子，以此来提高算法的寻优搜索效率。

文献[19]首次提出了相似度移除算子的概念，本文在考虑风险值、运输成本及约束条件的基础上，对相似度Hij作如下定义：

其中，r ij为客户(i,j)间的路段风险，d ij为客户(i,j)间的距离，|Eai-Ea j|为客户左时间窗的差值，为客户间需求量的差值。由定义可知H ij越小，两个客户越相关。首先随机选取一个客户，计算与其他客户相似度，依次删除n个删除相似度最大的客户。α1、α2、α3、α4、为各个影响因素的权重，本文取6、5、4、3[20]。

修复算子：在修复解时采用贪婪启发式插入方法，依次将删除的客户插入使总风险增量最小的位置，执行n次迭代，直至全部插入被删除的客户形成新的解。在修复操作时，解的成本包含总风险和惩罚成本两部分，当产生的解不满足时间窗和载货量约束时则在目标函数添加相应的惩罚成本，两种约束惩罚成本取相同值。

4 数值实验及分析

选用solomon标准测试算例的c101中50个客户点进行数值实验，用混合NSGA-II算法进行500次迭代实验，算法在操作系统为Windows10家庭版的计算机上基于Matlab 2016a软件实现。参数设置为：种群规模N=20，交叉概率Pc=0.9；变异概率P m=0.05，危险品运输车辆为25辆，人口密度上下界为区间[0，800]均匀分布的随机数，行驶速度上下界为[40，80]，客户需求量为区间[0，50]均匀分布的随机数，所有随机数均使用Matlab随机生成。其余实验参数取值见表2。

表2 参数取值Table 2 Parameter values

4.1 结果分析

用混合NSGA-II算法求解上述模型，得到的各目标最优时车辆路径如图4所示，各目标最优值及路径安排见表3。

表3 各目标最优值及路径集合Table 3 Optimal value and path set of each target

图4 混合NSGA-II算法各目标路径优化示意图Fig.4 Schematic diagram of each target path optimization of hybrid NSGA-II algorithm

多目标优化问题的各个目标函数相互影响，无法同时取得最优解，一个目标函数最优时，其余目标函数只能取得对应的非劣解。利用混合NSGA-II算法可以获得危险品运输车辆路径选择的一系列Pareto解，由表3可知，当Z1取最优时，对应的Z2和Z3目标函数值分别比其各自最优值增加了0.9%和28%，配送路径如图4（a），说明追求总风险最小需要决策者根据路段风险值合理安排车辆路径和车辆数目；当Z2取最优时，对应的Z1和Z3分别比其各自最优解增加了30.7%和28%，配送路径如图4（b），说明行程并不是影响风险的最重要因素，路段上的风险值大小对风险的影响更加显著，因此决策者可以考虑绕行部分风险小的路段以降低运输风险；当Z3取最优时，对应的Z1、Z2分别比其各自最优解增加了15.4%和8%，配送路径如图4（c），说明合理的路径安排可以在减少车辆使用数的同时不使车辆行驶距离显著增加，从而有效降低运输成本。

对求得的Pareto解集进行横向分析，考虑不同偏好决策者的需求，可以为多方面的运输参与者提供不同因素考量的危险品运输路径决策。分析表3中数据可知：（1）从运输部门安全的角度考虑，J1总风险最小的路径安排是运输部门的最佳选择，此时车辆总行程和车辆使用数取得非劣解，该非劣解相对于最小总行程和最小车辆使用数分别增加了0.9%和28%；（2）从承运人追求总成本尽可能小的角度考虑，则J2车辆使用数目最小的路径集合是最优方案，因为每多一辆车，就会增加司机工资、车辆固定使用成本及车辆油耗等多种成本，车辆使用数越多的方案，成本也将越高[21]。虽然J2方案的总行程比J3方案的多8%，但J3方案增加了2辆车的成本将远远高于增加的8%行程的成本；（3）从物流外包方的角度考虑，则J3总行程最小的路径方案是最好决策，是因为对于物流外包方来说希望租用出去的车总行驶距离最小且车辆周转时间最小，以便服务下一个需求客户从而获得更大利益。

4.2 算法性能分析

为对比NSGA-II算法和本文设计的混合NSGA-II算法的算法性能，以本文算例为基础算例，其他实验条件及参数设置不变，分别执行两种算法各10次，生成的Pareto前沿分布如图5所示，按车辆数分类的Pareto前沿分布如图6所示，算法性能统计见表4。

表4 不同算法求解结果统计表Table 4 Statistics of solution results of different algorithms

图6 不同算法按车辆数分类的Pareto前沿分布Fig.6 Pareto frontier distribution of different algorithms classified by number of vehicles

通过上述图表，从Pareto解数目、解的分布、求解精度、算法复杂度四个方面对算法性能进行对比分析：

（1）比较图5（a）和图5（b）可得，NSGA-II算法得到了8个Pareto解，混合NSGA-II算法得到了18个Pareto解，混合NSGA-II算法得到的Pareto解更多，解集的空间分布均匀且具有收敛性和鲁棒性，说明该算法可以搜索到较大的解空间，算法搜索性能得到了提高。

图5 不同算法的Pareto前沿分布Fig.5 Pareto frontier distribution of different algorithms

（2）由图6可得，混合NSGA-II算法所得的解中车辆使用数为7、8的解较多，相比于NSGA-II算法，混合NSGA-II算法得到的解分布更均匀，差异性更大，该算法求得的解质量更高。

（3）由表4可知，混合NSGA-II算法得到的最优总风险、总行程及车辆使用数目分别比NSGA-II算法优化了11.5%、1.0%和14.3%，平均总风险、总行程及车辆使用数目分别比NSGA-II算法优化了1.3%、5.4%和-13.6%，可见该算法求解精度明显更高。

（4）算法复杂度方面，NSGA-II算法采用快速非支配排序法，算法总体复杂度为O(mN2)，其中m为目标函数个数，N为种群大小，本文中m=3，N=20。混合NSGA-II算法依旧沿用NSGA-II算法的排序方法，算法复杂度上没有改变，但在寻优过程中增加了修复和删除算子进行改进，使得算法没有增加算法复杂度的同时提高了搜索性能。

综上本文提出的混合NSGA-II算法性能上总体优于NSGA-II算法。

5 结语

为提高危险品运输车辆路径决策的安全性、准确性和可靠性，本文以模糊变量来刻画运输过程中的不确定影响因素，引入基于动态载货量的风险评估模型，构建了不确定环境下的危险品运输车辆多目标优化模型，设计了NSGA-II算法与LNS算法相结合的混合NSGA-II算法来求解模型。最后，采用solomon算例验证了模型和算法的正确性和有效性，得出如下结论：

（1）模型综合考虑了安全性、经济性、时效性以及运输过程中的不确定因素，所提出的多目标优化模型更加贴近实际运输过程，能为决策者提供更加准确的决策支持，具有较强的实用性。

（2）多目标优化问题往往不能求得满足所有目标的最优解，现实生活中的各目标间甚至存在相悖的关系。因此在考虑单个决策目标最优时，其他目标只能取得相应的非劣解，本文以三方运输参与者的不同决策方案说明了模型结果可以为决策者提供不同偏好的决策支持。但在实际决策情境中，决策者在路径选择时一般会兼顾多个方面，即在不同优化目标间取得兼顾各方利益的折衷解更符合实际应用背景。

（3）设计了NSGA-II算法与LNS相结合的混合NSGA-II算法求解提出的多目标优化模型，并得到了Pareto最优解集。与NSGA-II算法相比，混合NSGA-II算法复杂度仍为O(mN2)，该算法在复杂度没有增加的同时在Pareto解数目、求解精度等方面的性能均表现得更好。

本文仅研究了不确定环境下单车型的危险品运输车辆路径优化问题，具有一定局限性，更具有实际研究背景的多车型危险品运输车辆路径优化问题将是下一步的研究方向。