多目标低碳车辆路径优化模型及求解算法*

2020-06-17李文霞张春民马昌喜

交通信息与安全 2020年1期

李文霞张春民马昌喜

（兰州交通大学交通运输学院兰州 730070）

0 引言

自2009年12月哥本哈根世界气候大会以来，低碳环保充分引起了全球节能减排意识，中国政府也提出了控制企业排放、建立碳税政策和企业碳排放权利交易等相关措施以减少运输业碳排放。运输业作为温室气体排放量最大的产业之一，其全球温室气体排放占总排放量的24%，公路运输作为物流配送的主要运输方式，排放量占到总量的17%[1]。为进一步减少物流配送带来的大量CO2排放，低碳条件下的车辆路径问题（vehicle routing problem,VRP）逐渐成为物流配送行业发展的必然趋势，对运输业减少碳排放具有十分重要的意义。

随着全球能源问题的日益严峻，低碳环保问题引起了各行业的重视，低碳车辆路径问题（green vehicle routing problem,GVRP）也引起了交通领域学者的广泛关注[2-4]。低碳车辆路径问题最早由Palmer[5]提出，该车辆路径优化模型中考虑车辆行驶速度及道路行驶条件对车辆碳排放的影响，但未研究车辆负载对油耗量的影响。近年来，有更多学者结合物流配送的实际情况，对低碳条件下的物流配送进行了深入研究。张得志[6]构建了基于油耗模型的多车型车辆路径优化模型，并基于遗传算法启发式进行求解，实例分析结果表明考虑碳排放的车辆路径行驶里程较长，但其综合成本较低。李进等[7]以车辆行驶速度作为决策变量，考虑油耗成本、碳排放成本和旅行时间费用最小化，建立基于碳排放及速度优化的车辆路径问题，并深入分析了碳排放与旅行时间之间的替换关系。葛显龙等[8]针对不同碳交易机制分别建立低碳配送模型，深入研究了碳限额、碳税率、碳交易3种机制对配送总成本、配送路径及碳排放量的影响。E.Jabir等[9]建立了多车场车辆路径问题模型（multi-depotgreen vehiclerouting problem,MDGVRP），所提出的基于蚂蚁系统的元启发式算法在合理的范围内为大规模MDGVRP实例提供了近似最优解。由于GVRP是经典VRP的衍生问题之一，启发式算法是解决该类问题的有效手段，尤其是面对大规模优化问题，启发式算法表现出良好的求解效率，但在求解精度、收敛性等方面仍具有很大的改进空间，因此，不少学者对一些算法有针对性地进行了改进[10-12]。

通过横向对比可知，既有文献目前存在以下问题：①油耗及碳排放量的计算未能综合考虑行驶距离、载重、道路因素、车辆因素等多个影响因素。②通过均值法或直接标定法对辆行驶速度进行估算，所损失的交通信息量较大，进而对配送时间及燃油消耗量的计算产生影响。③既有研究中的优化目标多以单目标为主[13-15]，考虑总运输成本最低，没有兼顾多方的利益。在实际配送过程中，运输企业的决策受多个因素影响，根据实际需求考虑多方面的多目标优化更具指导意义。因此，笔者首先考虑了行驶速度、车辆负载、路面属性及车辆属性对油耗及碳排放的影响，依据车辆行驶速度与燃油消耗量间的相互作用关系进一步改进综合油耗模型，其中车辆行驶速度采用三角分布函数在合理速度区间进行估算；其次为了提高企业的配送效率，以车辆周转时间作为设备周转速度的量化依据。基于以上分析，以总成本最低和车辆周转时间最少建立多目标低碳车辆路径优化模型，并针对该模型设计了增强型非支配排序遗传算法（enhanced non-dominated sorting genetic algorithm-ii,ENSGA-II），提高了模型的求解质量。

1 多目标低碳车辆路径模型建立

1.1 问题描述及相关参数定义

多目标低碳车辆路径问题（multi-objective low carbon vehicle routing problem，MO-LCVRP）可描述为：单个配送中心派遣异质车队为多个客户点配送货物，以总成本和车辆周转时间最小建立多目标优化模型，合理规划运输路径。为进一步明确本文研究范围，提出以下假设。

1）单个配送中心向多个客户点配送货物，地理坐标位置均已知，完成配送任务后，配送车辆全部返回配送中心。

2）各客户点需求量已知，且货物需求不可拆分。

3）每个客户点仅由1辆车服务，配送过程中仅考虑送货服务。

4）单个客户点的货物需求量不超过配送车辆的最大容量。

5）采用多车型进行配送，每种车型的载重量和固定启用成本均不相同。

6）车辆必须在客户规定的时间窗内提供服务，违背时间窗会带来额外的惩罚成本。

7）不考虑拥堵，车辆行驶速度通过三角分布函数进行估算。

模型相关变量及参数定义见表1。

1.2 基于三角分布的车辆速度量化

由于路网中不同路段的路况具有差异性，导致车辆在各路段的旅行速度（以下简称弧速度）不尽相同。而车辆行驶速度又与车辆油耗密切相关，因此建立优化目标前需要对速度进行合理的量化。既有研究中普遍认为弧速度在已知的最大值与最小值之间是均匀分布的，采用均值法或在给定的区间随机生成来标定弧速度[16]，但该方法仅考虑2个边界值，所包含的信息较少，所得的估算值与实际值差距较大，过于理想化，不足以反映实际路况，为此，本文拟采用适用于随机建模的三角形概率分布模型估算弧速度。

表1 模型相关变量及参数定义Tab.1 Model related variables and parameter definitions

三角形分布除了有上界a和下界b的限定外，还包含了最可能出现的值c，是一种连续概率分布模型，与正态分布相似，但比正态分布更灵活、更直观，即如果c的值更接近a或b，则三角形分布也会偏斜，见图1。三角分布虽然是一种对于实际情况的简化，但这种简化所损失的信息量较小，由此所得到的结果与真实情况相差较小。尤其适用于样本数据不全但已知变量之间关系的情景，应用于在最大值与最小值之间出现最可能值，即已知模型中样本最可能出现结果的概率预测模型。

图1 三角分布示意Fig.1 Illustration of triangle distribution

在现实中，车辆不可能以恒定速度行驶。根据先验信息，有可能在最大速度vmax和最小速度vmin内行驶。得益于目前交通大数据的支持，可知车辆在最小速度和最大速度之间的最可能行驶速度，因此，使用三角形分布能更好地量化弧速度，能够使速度在服从概率分布的情况下具有一定的随机性，相较于均匀分布和正态分布，其应用更为灵活，适应性更强。因此，引入三角分布来量化弧速度约束，可以有效的结合自由流速度vu、拥堵速度vd和可能速度vm这3个速度信息，对期望速度进行估算，为进一步分析车辆速度与燃料消耗之间的关系奠定基础。三角分布概率密度函数图像见图1，表达式见式（1）。

可以看出，三角分布概率密度函数是分段函数，通过式（2）对速度在不同的取值区间进行积分，求得速度的期望值。

式中：p和q为速度区间的上/下限。

1.3 优化模型建立

1.3.1 优化目标

根据相关研究，车辆的碳排放与油耗呈正相关关系。因此，以燃油消耗量作为衡量车辆碳排放的标准。车辆油耗量的计算有瞬时排放模型、载重碳排放模型、综合排放模型、MEET模型等[17]，其中综合排放模型全面考虑了车辆在行驶过程中运行速度、货物负载重量、行驶道路路面因素等的影响。但在实际的应用过程中，综合模型的计算过于复杂，需要进行测算的数据较多，根据文献[17]的相关研究，车辆在行驶过程中，影响油耗量大小的决定因素随速度的变化而改变，一般货物运输车辆行驶速度与燃油消耗关系见图2。

图2 车辆行驶速度与燃油消耗关系Fig.2 Relationshipbetweenvehiclespeedandfuelconsumption

由图2可知，发动机系统做功和牵引力做功共同决定油耗量的变化。当车辆行驶速度小于40 km/h时，燃油消耗量以发动机做功为主；当车速高于40 km/h时，燃油消耗以牵引做功为主，并且随着车速的不断增加，牵引力做功比重不断上升，发动机系统做功比重不断下降。当车辆行驶速度大于60 km/h时，牵引力做功引起的燃油消耗可占到总燃油消耗量的70%以上，因此，不考虑拥堵情况，假定车辆行驶速度大于40 km/h，以牵引力做功引起的燃油消耗近似代替总的燃油消耗量，进一步简化计算，并求得碳排放量。

当车辆由运行弧(i，j)到运行弧(j，k)时，其运行道路的倾斜角度，车辆迎风面积等均为定值，在此，仅考虑车辆旅行速度及负载量的变化。车辆在弧(i，j)上的期望速度为vij，总重量为W（包括车辆自重ωm和负载量ωlij），路面倾斜角度θ，弧上里程为dij，则弧上耗油量计算见式（3）。

式中：p1为车辆上的总牵引力需求；αij=a+gsinθ+gCrcosθ为弧相关参数，与车辆加速度a、坡度θ、重力加速度g和滚动阻力系数Cr有关；β=0.5CdAρ为车辆属性相关的参数，与迎风面积A、空气密度 ρ和空气阻力系数Cd有关；其中，表征了由负载变化引起的燃油量消表征了由速度变化引起的燃油量消耗。

基于以上油耗模型的分析，综合考虑配送过程中产生的油耗成本、碳排放成本、固定成本及违背时间窗惩罚成本，建立总经济成本最小的优化目标。

为加快设备的周转速度，以车辆周转时间来有效量化车辆配送效率，模型的第2个优化目标为周转时间最小化。由于运输任务采用的是多个子路径同时配送模式，因此，车辆周转时间为耗时最长的子路径完成配送任务的时间为

1.3.2 约束条件

式（6）表示每个客户点仅被1辆车服务且仅服务1次；式（7）表示车辆从仓库出发最终返回仓库，且启用车辆数不超过固有车辆数；式（8）表示车辆配送过程中服务的连续性；式（9）为车辆容量约束；式（10）为车辆的弧段运量约束；式（11）为配送过程中时间限制条件；式（12）表示消除子回路；式（13）为决策变量约束；式（14）为非负变量约束。

2 求解算法设计

NSGA-II算法[18]是求解多目标优化问题一种有效算法，但经典NSGA-II算法在求解VRP及其衍生问题时表现出搜索能力弱，收敛性不好等，因此，本文对NSGA-II算法进行改进。以Ong等[19]于2016年提出的多因子优化算法（multifactorial evolutionary algorithm，MFEA）为指导思想，并针对不同的优化目标设计了相应的邻域搜索算子，MFEA被应用于NSGA-II父代种群的优化过程，基于此，提出了求解多目标的增强型非支配排序遗传算法（enhanced non-dominated sorting genetic algorithm,ENSGA-II）。

2.1 改进ENSGA-II算法

笔者提出的ENSGA-II算法的改进重点是：以NSGA-II算法为架构，对第l代种群Pl执行精英保留策略得到第l+1代种群Pl+1的父代集合parent_set()，对parent_set()进行多因子优化，改善算法的局部搜索能力，并利用多个目标之间的协同进化提高算法的收敛性。

2.1.1 改进多因子优化方法

多因子优化打破了传统多目标算法仅在统一的解空间搜索折衷解的限制，提出基于不同的技能因子（优化目标）在不同的解空间中搜索可行解，即同时进行多个优化任务，最后具有不同技能因子的个体间执行交叉操作，为跨任务的个体间的信息交换提供机会[20]。

尽管多因子优化表现出较好的算法性能，但还存在一定的制约，尤其当技能因子间不满足协同性和相关性时，可能会出现个体间信息的负迁移效应。本文模型中考虑的碳排放成本和周转时间在旅行距离上具有相关性，保证了信息交换的正迁移效应，因此，以多因子优化为指导思想设计了求解GVRPSTW问题的MFEA算法，其中多优化任务的划分和选型交配阶段的变异、交叉算子具体如下。

步骤1。多优化任务。针对1条完整的可行路径R，将其划分为多条子路径ri，每条子路径作为独立的优化任务被分配1个技能因子τh,h=(1,2,…,H)，见图3。

图3 多优化任务划分Fig.3 Dividing multiple optimization tasks

步骤2。变异操作。

1）Reverse算子。碳排放与车辆在旅行弧上的载重直接相关，在此定义了全局弧载重，来阐明Reverse算子在绿色车辆路径中的适用性。全局弧载重Gd是指对全局每个旅行弧上的载重量di与旅行弧里程Li的乘积的求和，即Reverse算子的优势在于路径r与路径rˉ具有相同的旅行距离成本，但具有更低的全局弧载重，见图4。

Reverse算子利用了路径r与其反向路径rˉ全局弧载重不同的事实，比较图4(a)和图4(b)可看出，全局弧载重Gdr＞Gdrˉ，因此，Reverse算子在优化路径碳排放成本上表现出较好的效果。

2）Relocate算子。针对路径r，随机产生移出节点u1和插入节点u2，将u1对应的元素插入到u2位置，见图5。

图4 路径与反向路径的Gd值的比较Fig.4 Comparison ofGdvalues of path and reverse path

图5 插入操作及对应路径Fig.5 Insertion operation and corresponding path

步骤3。交叉操作。

1）Exchange算子。分别对路径r1和路径r2随机产生交换节点u1和u2，将路径r1在节点u1对应的元素与路径r2在节点u2对应的元素交换，见图6。

图6 单点交换操作及对应路径Fig.6 Single-point exchange operation and corresponding path

2）Cross-exchange算子。对路径r1随机产生2个交换节点u1和u2，对路径r2随机产生2个交换节点u3和u4，将节点u1和u2之间的序列sq1与节点u3和u4之间的序列sq2交换，见图7。

图7 匹配交换操作及对应路径Fig.7 Matching exchange operation and corresponding path

具有相同技能因子的个体执行交叉操作或单个体执行变异操作是在统一的优化方向上的进化，具有不同技能因子的个体执行交叉操作是在不同优化方向上的信息交互，以期达到不同优化目标的协同进化效果，增强算法搜索能力。所设计的多因子优化过程不破坏解的可行性，不需要对解进行修复。

2.2 算法流程

图8 ENSGA-II算法流程Fig.8 Flowchart of ENSGA-II algorithm

3 算例分析

3.1 算例数据及参数设置

为了模拟需求点位置的随机性，测试算例配送中心及需求点位置通过matlab随机生成，需求点需求量是区间[0,10]均匀分布的随机数，需求点时间窗是区间[0,50]均匀分布的随机数；每单位货物的卸货时间为2 min/t，需求点信息见表2；车辆由配送中心出发，向19个需求点提供服务，车辆其他相关信息见表3；为了表征路网属性的随机性，最可能的车速在区间[40,100]范围内随机生成，得到服从三角分布的速度矩阵见表4；针对模型中特定弧（i,j）参数αij的取值，参考文献[15]在区间[0.09,0.15]随机生成，单位油耗费用为7.47元/L，单位油耗碳排放费用为0.64元/L。

表2 配送中心和需求点信息Tab.2 Distribution center and demand point information

表3 配送车辆信息Tab.3 Distribution vehicle information

表4 服从三角分布的速度矩阵Tab.4 Speed matrix obeying triangular distribution km/h

3.2 模型计算结果及分析

3.2.1 结果分析

采用所提出的ENSGA-II算法在Matlab R2016a软件上编程对模型进行求解，设置算法循环迭代Max_Cycle=200，多因子优化过程各优化任务最大循环迭代Max_subCycle=20，初始种群规模pop=200；交叉概率为0.9，变异概率为0.1。算法在主频为1.6 GHz的Intel Core i5处理器、操作系统为 Win 10的环境下进行，最优车辆路径规划见图9，目标f1(x)和目标f2(x)对应的极值解及车辆路径见表5。

图9 最优车辆路径规划Tab.9 Optimal vehicle path planning

表5 目标极值及对应路径集合Tab.5 Target extreme value and corresponding path set

表5中碳排放占比能够有效反映节能减排成本在总运输费用中的占比情况，碳排放占比=碳排放成本/（碳排放成本+燃油成本+惩罚成本+车辆启动成本），以表5中成本最优时的子路径1为例，42.01/（42.01+490.28+486.37+930）=2.16%。

由于本文针对多目标进行优化，因此，所求得的结果为Pareto解集，即使得2个目标函数同时达到最优的解集不存在。在实际决策过程中，决策者在选择路径时一般会兼顾多个方面，即在不同的优化目标之间取折衷解，这更符合兼顾各方经济效益的原则和路径优化的实际背景，因此，由表5可知，当运输商更偏好总经济成本最优时，其对应的车辆周转时间取得非劣解，该非劣解相较于最小周转时间增加了26.8%，此时车辆平均满载率、平均碳排放成本、平均燃油成本、平均旅行时间分别为83%，20.43元，238.52元，87.27min；当运输商更偏好车辆周转时间最优时，对应的总经济成本取得非劣解，该非劣解相较于最低经济成本增加了14.8%，此时车辆平均满载率、平均碳排放成本、平均燃油成本、平均旅行时间分别为83.67%,31.94元，372.89元，73.33 min。由此可见，在平均燃油和碳排放成本方面，由于碳排放与燃油消耗成正比例关系，其优化结果相统一，总经济成本最优时比车辆周转时间最优时降低了36.04%；在平均旅行时间方面，车辆周转时间最优时比总经济成本最优时降低了15.97%。

由表5的结果，可以得到以下结论。

1）车辆在运输过程中应优先考虑重型车辆。其具有单位重量的运行成本小的出行优势，当重型车辆数无法满足运输需求时，考虑使用轻型车辆。

最佳适应算法BF(Best Fit):在装入货品时装入到最合适这个货品的箱子里，这个箱子不是第一个可装的箱子，而是最合适的。当没有适合该物体的箱子时，打开一个空箱子。

2）在配送过程中应当尽量减少车辆的启用。尽管多启用1辆车进行配送会使得时间窗满足率提高以及车辆周转时间减少从而带来成本和时间上的节省，但同时也会带来车辆成本的增加，由于车辆固定启用成本在经济成本中占有较大比重（50%～60%），因此，在成本和时间上的节省不足以弥补车辆成本的增加。

3）适度提高碳排放价格有利于节能减排。无论是经济成本达最优还是车辆周转时间达最优，碳排放成本在总成本的占比最高不超过3%，这是由于目前中国正处于工业快速发展的时期，而单位碳排放的成本价格仅占单位油耗价格的8.57%，因此企业在制定物流配送路径时对碳排放的考量较少，为进一步落实中国向环境友好型社会的转变，政府应当在一定合理范围内提高碳排放成本的价格，使得配送商能够充分考虑碳排放成本并合理规划配送路径。

基于以上分析，一方面，由于政府强调环境保护，使得配送商更倾向于选择在满足政府节能减排政策下配送成本最小的运输方案；另一方面，为了提高企业经济效益，配送商期望选择车辆周转时间最短的运输方案，因此，在配送方案的选择上，可根据不同的情景要求，优化一个方面的同时兼顾另一方面取得非劣解，选择合理的运输路径。

3.2.2 算法性能分析

为测试ENSGA-II算法的性能，以本文算例为基础试验数据，算法参数设置同3.2节，分别对ENSGA-II算法和NSGA-II算法运行10次，2种算法Pareto前沿对比见图10，统计结果及算法改进百分比见表6。

图10 2种算法Pareto前沿对比Fig.10 Comparison of two algorithms Pareto frontier

表6 不同算法求解统计结果及改进百分比Tab.6 Statistical results and improvement percentages of different algorithms

经过10次算法测试后，将ENSGA-II算法与经典NSGA-II算法在求解精度和收敛性2个方面进行对比。由图10中2种算法Pareto前沿对比图可以知，相较于经典NSGA-II算法，ENSGA-II算法的Pareto解集中最优解数量更多，并且Pareto前沿具有更好分布性，在求解质量和收敛速度上表现出较强的优势。由表6可知，在经济成本方面，最优路径、最差路径和平均路径的改进均超过3%，其中最差路径的改进程度达到4.45%，并且在平均收敛次数上提前了44个迭代循环，提高了28.03%；在车辆周转时间方面，平均路径的改进值最为明显，超过5%，最优路径及最差路径也有一定程度上的改进，并且在平均收敛次数上提前了25个迭代循环，提高了15.24%。基于以上分析可知，本文提出的ENSGA-II算法充分发挥了多因子优化的协同进化和信息交互的优势，表现出更好的搜索性能和求解效率。