带软时间窗的冷链电动汽车路径问题研究

2022-07-14刘志硕李秋雨董子琦

北京交通大学学报 2022年3期

刘志硕，李秋雨，董子琦，陈哲

（北京交通大学交通运输学院，北京 100044）

在食品生鲜配送领域，传统的冷链配送通常采用燃油汽车，但是燃油燃烧产生的二氧化碳和不完全燃烧产生的有害气体会加重环境的负担，生鲜产品的制冷运输也会消耗更多能源.与燃油车相比，电动车以电能为动力源，使用过程零排放，对环境和社会资源更加友好.有关研究表明，相比燃油车，采用电动车进行物流配送能够有效降低能源消耗和运营成本［1-2］.因此，将电动车应用于冷链物流配送领域是未来的发展趋势.目前我国已有一些物流企业开始尝试将电动冷藏车投入到冷链配送环节，但电动冷藏车的应用仍然存在一些问题.电动冷藏车一般采用磷酸铁锂动力电池，虽然其比三元电池成本更低且更加安全，但存在能量密度低、续驶里程短、充电时间长等问题［3］.冷链物流配送过程中电动车需要制冷，因此需要消耗更多电能，而大部分城市仍存在社会充电桩密度小、覆盖不全面等问题.因此，如何利用现有资源合理规划路径，在提高配送效率的同时降低配送成本显得尤为重要.

针对冷链物流配送优化问题，国内外学者从冷链物流中多目标、动态需求、多温共配等角度研究了车辆路径问题［4-7］.Osvald 等［5］考虑到新鲜蔬菜易变质的特点，研究了在时间窗和行程时间约束下的冷链物流车辆路径问题.Chen 等［6］考虑了冷链产品配送中客户需求的不确定性，提出了带有生产调度的车辆路径问题.Tsang 等［7］提出了一种基于多温区共同配送的路径规划系统，该方案可降低食品变质率，缩短交付时间，同时提高客户满意度.

针对电动汽车路径优化问题，国内外学者主要研究了里程约束、充电方式、充电策略等问题，并开发了相关的高效智能算法［8-11］.Goeke 等［9］提出了传统汽车和电动汽车的混合车队下带时间窗的路径优化问题.Keskin 等［10］在研究电动汽车路径优化问题时放宽了完全充电的限制，提出了部分充电的策略.Verma 等［11］通过结合充电站充电和交换站更换电池两种方式来提高客户交付的准时率.

然而，目前鲜有针对电动车进行冷链物流配送路径问题的优化研究.主要原因可能是电动车在物流领域的应用尚处于起步阶段，社会充电桩等基础设施的建设还有待完善.但是随着电动汽车不断发展和其应用范围的不断扩大，电动车冷链物流配送路径优化将是一个非常值得研究的问题.因此，本文提出一种带软时间窗的冷链电动汽车路径问题（Cold-Chain Electric Vehicle Routing Problem with Soft Time Windows，CEVRPTW），构建以配送总成本最低为优化目标的数学模型，设计自适应大邻域搜索算法并进行求解.

1 模型建立

1.1 问题描述及基本假设

一个生鲜配送中心派出同质的电动车，向多个客户配送同一温区的生鲜商品，配送中心只在固定时间段内运营，车辆只能在该时段内出发或者返回.客户有软时间窗要求，即给客户配送商品，允许车辆早到或晚到，但会产生相应的等待成本或延误成本.车辆容量固定，续驶里程有限，中途可能需要在社会公共充电站充电，充电方式为快充.要求制定最优的配送方案，使配送总成本最低.

基本假设如下：

1）一个配送中心，运营时段固定，所有车辆出发前在配送中心充满电，中途可多次充电，且每次充满；

2）每辆车只允许配送一次且访问同一个充电站最多一次；

3）各客户有软时间窗，需求量已知，且不超过车辆最大容量；

4）每个客户只允许由一辆车配送一次；

5）温区、车型单一，充电速度恒定；

6）为提供车辆动力所消耗的电能与里程成正比；

7）制冷所消耗的电能与制冷的时间成正比.

1.2 模型变量

为解决CEVRPTW 问题，建立了一系列模型参数与决策变量，其设置及具体含义分别见表1 与表2.

表2 模型决策变量Tab.2 Decision variables of model

1.3 目标函数

CEVRPTW 问题以配送总成本最小为目标，它由固定成本和可变成本构成，其中可变成本包含运输成本、制冷成本、等待成本和延误成本.

1）固定成本包括车辆购置、支付员工工资等，被分配到每辆配送车辆上，计算式为

2）运输成本为车辆行驶过程中为提供车辆动力所消耗的电能成本，计算式为

3）制冷成本包括电动车通过制冷装备在行驶途中(不包括车辆完成配送任务回到配送中心的这段路程)、卸货期间和充电期间消耗的电能，计算式为

4）等待成本是车辆在客户最早接受服务时间之前到达目的地所产生的，它通过该时间段的制冷成本来衡量，计算式为

5）车辆在客户要求的最晚配货时间之后到达，则产生延误成本，计算式为

1.4 数学模型

目标函数：

约束条件：

式（7）表示车辆给客户配送后必须离开；式（8）表示车辆从配送中心出发，最后再返回配送中心；式（9）表示每个客户只允许访问一次；式（10）表示车辆容量约束；式（11）表示每辆车在配送过程中允许多次充电；式（12）表示车辆的行驶时间；式（13）表示给客户的服务时间，即卸货时间；式（14）表示车辆在充电站的充电量；式（15）表示车辆在充电站的充电时间；式（16）表示车辆在节点的等待时间；式（17）表示车辆在节点的延误时间；式（18）表示车辆在节点开始接受服务的时刻；式（19）表示车辆到达任意节点的时刻；式（20）表示车辆到达任意节点的装载量为非负值；式（21）表示车辆在开始配送任务之前的装载量不超过车辆的最大容量；式（22）表示客户i与其他节点j的时间关系；式（23）表示充电站i与其他节点j的时间关系；式（24）表示客户或充电站与配送中心的时间关系；式（25）表示客户及配送中心满足的时间窗要求；式（26）～式（28）表示车辆从某个节点出发并访问客户后有足够电量到达其他节点；式（29）表示车辆在充电站充完电后有足够电量回到配送中心；式（30）表示两个变量之间的关系；式（31）表示决策变量取值范围.

2 求解算法

上述模型是一个0-1 整数规划模型，属于车辆路径问题（Vehicle Routing Problem，VRP）的范畴.VRP 问题是一个典型的NP-Hard 难题，当节点规模很小时，可以采用数学规划方法求得精确解.但随着节点规模增加，解空间呈指数级增长，采用数学规划方法无法求解，因此，一般采用启发式方法或元启发式方法求解.Elshaer 等［12］对求解VRP 问题的元启发式算法进行了统计分析，结果显示，最受学者青睐的算法有禁忌搜索算法、自适应大规模邻域搜索算法（Adaptive Large Neighborhood Search，ALNS）和蚁群算法等.因此，本文中针对CEVRPTW 问题特点设计了ALNS 算法求解.

2.1 ALNS 算法

ALNS 算法是大邻域搜索算法的延伸［13］，可用于解决各种车辆路径问题.它使用多种破坏和修复算子，每次迭代时采用随机方式来选择要执行的算子，而选择概率由该算子之前的表现来决定，历史表现越好，被选中的概率越高.

本文在Shaw 等人［13-15］对ALNS 算法研究的基础上，增加了针对CEVRPTW 问题的相关移除算子和插入算子.算法总体思路是首先构造初始解，然后尝试采用不同邻域算子对其进行改进，不断迭代，直到满足算法终止条件.在每次迭代中，先随机选择某个移除算子从线路中移除一些节点（客户、充电站），以破坏当前解，然后选择插入算子来修复该解.插入算子尝试性地将移除的客户或充电站节点再插回到线路中，以获得比之前更好的解.

2.2 构造初始解

采用贪婪插入法构造初始解，步骤如下：

1）初始化一条新线路，不含任何客户节点和充电站节点.

2）从当前所有未插入到线路的客户节点中选取地理位置上离配送中心最远的客户节点，作为第一个节点插入到新线路中.

3）判断当前是否还有未插入的客户节点，如果仍有未插入到线路中的客户节点，转4）；否则算法终止.

4）判断当前所有未插入到线路的客户节点中，是否存在满足插入到当前线路条件的客户节点.若存在，则从这些客户节点中选择一个使线路成本增加最少的节点插入到当前线路，并跳转到3）；否则转1）.

关于判断节点是否满足插入条件，具体指的是要满足车辆容量约束以及配送中心的硬时间窗约束.如果客户节点满足以上两个条件，却因电池电量低无法插入，则可以同时插入客户节点和一个充电站节点.节点插入到线路的成本，指的是分别计算该节点插入到当前线路所有各个位置时的成本，以最小成本所在的位置作为该节点的插入位置.

2.3 移除算子

2.3.1 客户节点移除

本文中采用Shaw、Avci 和Masson 等设计的随机移除、相关移除、最差成本移除、最差时间移除、区域移除和随机回路移除算子［13-15］. 此外，针对CEVRPTW 问题中的软时间窗约束设计了最长延误客户节点移除算子.

1）随机移除（CRR）：从当前解中随机移除若干客户节点.

2）相关移除（CSRR）：根据客户节点之间的相似度进行移除，客户节点i和j的相似度R(i，j)根据时间窗要求、距离和需求量计算，权重分别用α、β、γ来表示，计算式为

式中：di0和dj0分别表示客户节点i和j距离配送中心的距离；ei和ej分别表示客户节点i和j的最早服务时间；li和lj分别表示客户节点i和j的最晚服务时间；qi和qj分别表示客户节点i和j的需求量.R(i，j)越小，客户节点i和j的相似度越高.

3）最差成本移除（CWCR）：移除使当前解增加成本最多的若干客户节点，即若移除这些节点，能够使当前解的成本降低幅度最大.

4）最差时间移除（CWTR）：本文中考虑了客户节点的软时间窗要求，时间代价是指车辆服务该客户节点的等待或延误所造成的成本.因此移除等待或延误时间最长的客户节点.

5）区域移除（CZR）：节点所在的笛卡尔坐标系被划分为许多较小区域.随机选择一个区域，将该区域中所有客户节点从当前解中移除.

6）随机回路移除（CRRR）：随机选择若干条回路并删除回路中所有客户节点.回路数量的取值范围为当前解回路总数的10%～12%之间.

7）最大延误客户节点移除（CBLR）：客户节点有软时间窗要求，晚到会产生相应的延误成本，且这部分成本较高，因此将延误时间最长的客户节点移除.

由于CEVRPTW 问题针对的是电动车，行驶过程中必须考虑电动车电量，所以在移除客户节点后，一些充电可能不再需要.因此在客户节点移除过程中，存在客户和充电站节点关联移除.即如果某个客户节点被移除，恰好访问该客户节点的电动车是从充电站节点再过来的，则同时移除该客户节点和充电站节点，另一种情况是恰好访问完某个客户节点后需要去下一个充电站节点充电，如果该客户节点被移除，那么接下来要访问的充电站节点也被移除.

2.3.2 充电站节点移除

充电站节点是EVRP 问题的独特之处，既要满足电动车行驶过程中的电量要求，又要避免多余的充电时间.因此考虑充电站节点的移除也可以在一定程度上优化初始解.

1）随机充电站节点移除（SRR）：与客户节点随机移除机制类似.

2）最差充电站节点使用移除（SWR）：移除线路中电动车充电量最少的充电站节点.该算子的目的是在需要充电之前尽可能多地利用电动车电池，提高充电站效率.

2.4 插入算子

2.4.1 客户节点插入

客户节点插入算子包括随机插入、贪婪插入、后悔插入算子［13-15］.

1）随机插入（CI）：在满足初始解生成过程中插入条件的前提下随机插入.

2）贪婪插入（CGI）：参照初始解生成的算法.

3）后悔插入（CRI）：用xik表示第k极插入成本的路径，当k≤k′时(k′∈K)，满足Δzxik≤Δzxik′，用后悔值ci=Δzxi2≤Δzxi1表示把花费成本最高客户节点i插入当前最优路径和次优路径之间的成本差值.插入最佳位置为在满足插入条件的前提下对应成本最小.

2.4.2 充电站节点插入

在移除充电站节点后，某些线路会因电动车电量不足而导致不可行，因而需要在合适的位置插入充电站节点，充电站节点插入如下：

1）最佳充电站节点插入（SBI）：若车辆到达客户节点的电池电量为负，则在该客户节点和前一个客户节点之间的弧上插入“最佳”（增加距离最小）充电站节点，如果此插入不可行，则以相同方式尝试之前的弧.

2）比较插入（SCI）：由于在选择充电站节点时遵循就近原则，可能会导致充电完成后，访问下一个客户节点时行驶更远的距离，花费更多的成本，因此，在一条回路内部，可以考虑让车辆提前充电.若提前充电成本较之前有所降低且电量能够支持访问回路里的剩余节点回到配送中心，则插入新的充电站节点.具体插入示意图见图1.

图1 比较插入算子示意图Fig.1 Schematic diagram of comparison insertion operator

ALNS 算法伪代码如下：

3 算例分析

针对CEVRPTW 问题的特点，对Goeke 构建的算例进行改造［9］，制作出CEVRPSTW 问题的算例集，以验证ALNS 算法的有效性.Goeke 的算例改编自Solomon 基准算例集［16］.所设计的算例集按节点数目分为大规模和小规模两个子集.其中，小规模算例有15 个，包含5、10 和15 个客户节点的算例各5 个；大规模算例有30 个，每个算例包含100 个客户节点、20 个充电站节点.采用Dev-C++编程实现本文中所设计的ALNS 算法，测试环境为CPU：3.4 GHz、内存：16 G.

3.1 调参试验

本文中的模型和算法涉及很多参数，包括CEVRPTW 问题自身的参数和算法的参数.对于问题自身的参数，一部分参考相关文献来设置，如车辆的固定成本、单位耗电成本等；一部分由笔者根据现实情况酌情设置，表示为⊙.对于算法的参数，一部分参考相关文献来设置，比如时间窗要求权重、距离权重和需求量权重、模拟退火的初始温度等，其来源直接标记为参考文献；一部分通过试验来优化设置，如ALNS 算法的迭代次数、因子得分增量等，其来源表示为⊕，其余由笔者自行设置，表示为⊙，具体如表3 所示.

表3 参数设置Tab.3 Parameter setting

对于上表中通过试验来优化设置的参数，选取小算例进行调参试验，设计如下：5、10 和15 个客户节点的算例各5 个.针对15 个算例分别进行ALNS算法迭代次数的参数试验、因子得分增量参数试验，得到最佳的参数配置，以便于后续的试验.

3.1.1 ALNS 算法迭代次数参数试验

本节对ALNS 算法的迭代次数进行试验.迭代次数分别为100 次、200 次和300 次.试验结果如表4所示，f表示运行10 次程序后得到的最优值，t表示求解的时长.

从表4 可以看出，ALNS 算法的迭代次数为200 次时最佳，仅迭代100 次无法找到最优解，迭代300 次花费时间太长.

表4 ALNS 算法参数试验结果Tab.4 Experimental results of ALNS algorithm parameter

3.1.2 因子得分增量参数试验

本节针对因子得分增量参数进行试验，通过试验效果对因子得分增量(θ1，θ2，θ3)进行调整控制.理论上来讲，算子的性能越好，越有可能搜索到更好的解，越有可能出现“优化后的解优于全局最优解”“优化后的解优于当前解，劣于全局最优解”这两种情形，该算子的得分就应该越高，越容易在下一次调用时被选中.也就是说，θ1＞θ2＞θ3才是最合适的参数方案.

为了验证这一假想，(θ1，θ2，θ3)共选择（40，30，20）、（40，20，30）、（30，20，40）、（30，40，20）、（20，30，40）、（20，40，30）共6 种参数组合方案进行试验，试验结果如表5 所示.

从表5 可以看出，各参数组合结果相差不大，其中（30，40，20）表现最好，（40，30，20）次之，（20，30，40）效果最差.

表5 因子得分增量参数试验结果Tab.5 Experimental results of increment factor parameter

试验结果并未完全支持θ1＞θ2＞θ3这一假想.比较θ2与θ3，方案（40，30，20）对于“优化后的解优于当前解，劣于全局最优解”给30 分，而对于“优化后的解劣于当前解，但被接受”给20 分.然而，该参数方案的平均最优值仅为488.79，优于方案（40，20，30）的491.89，即θ2＞θ3的方案要优于θ2＜θ3的方案.另外两组对比（30，40，20）vs（30，20，40），（20，40，30）vs（20，30，40），均得出了类似结论.显然，3 组对比均支持θ2＞θ3这一假想.

从θ1与θ2相比较来看，方案（40，30，20）对于“优化后的解优于全局最优解”给40 分，而对于“优化后的解优于当前解，劣于全局最优解”给30 分.然而，该方案的平均最优解高于方案（30，40，20），即θ1＞θ2的方案要劣于θ1＜θ2的方案.而另外两组对比（30，20，40）vs（20，30，40），（40，20，30）vs（20，40，30），却支持θ1＞θ2的方案.

因此，仅从试验结果来看，θ1＞θ2＞θ3这一假想不一定合理.由于在所有方案中，方案（30，40，20）的优化效果最好，因此本文后续试验均采用（30，40，20）的参数组合.

3.2 小规模算例试验结果与分析

为验证所建数学模型的正确性和ALNS 算法的性能，使用CPLEX 和ALNS 算法对小规模算例进行求解，其中ALNS 算法对每个算例执行10 次，取最优结果.试验结果如表6 所示，gap 表示ALNS所求得的解与CPLEX 解的偏差.

从表6 中可以看出，针对所有算例，ALNS 均能找到CPLEX 的精确解，且花费的时间远小于CPLEX.其中15 个客户节点的最优解对应的配送路线为1-3-13-6-19-10-9-16-5-2-1，1-4-17-11-14-1，1-8-7-18-15-12-1（节点1 表示车厂，节点17-21表示充电站，其余为客户节点），如图2 所示.

图2 15 个客户节点的最优解对应的配送路线图Fig.2 Distribution route map corresponding to the optimal solution of 15 customer nodes

表6 小规模算例实验结果Tab.6 Small-scale experimental results of algorithm example

3.3 大规模算例试验结果与分析

针对大规模算例，将节点的分布类型分为C 类、R 类、RC 类，C 类算例中节点的分布形态为集群式分布，R 类算例中节点的分布形态为随机分布，RC算例中节点的分布形态则结合C 类分布和R 类分布的特征.每类算例又根据时间窗要求的不同划分成两种.其中C 类算例有10 个，R 类算例有10 个，RC类算例有10 个，每类分布形态包括宽和窄时间窗的算例各5 个.针对每个算例，运行10 次，得到初始值f1、平均值f2、最优值f3和程序的运行时间t，试验结果如表7 和图3 所示.

从表7 可以看出，所设计的ALNS 算法解决大规模的CEVRPTW 问题是可行的，且求解效果较好，但是代码运行时间较长.

表7 大规模算例试验结果Tab.7 Large-scale experimental results of algorithm example

从图3 可以得出以下两个结论：一是不同节点分布形态的求解结果有较大差异，特别是在C 类算例的路径方案中车辆单次配送能够访问更多的客户节点，车辆装载率较高，致使C 类算例的平均值少于R类和RC 类算例，且C 类算例运行时间更短.二是每类算例中不同种算例之间的目标值存在较大的差异，这是由于C2 系列、R2 系列、RC2 系列算例的时间窗要求比C1 系列、R1 系列、RC1 系列的更加宽松，车辆在一次配送任务中访问的客户节点较多，车辆装载率较高，使用车辆数较少，其配送成本也随之减少.

图3 大规模算例试验结果比较图Fig.3 Comparison chart of large-scale algorithm examples’experimental results

3.4 ALNS 算子表现分析

研究了所设计的ALNS 算子的性能，如图4～图6 所示.每个图的横轴表示ALNS 算法的进度：算法从0%开始，到100%结束.图4 展示了每个算子的累积贡献率，用算子发现新的可接受解的次数表示，并用找到新的可接受解的总次数（百分比）来归一化.实线表示移除算子，虚线表示插入算子.从图中可以看出，在算法搜索的早期就发现了较大数量的改进解，搜索快结束时图形逐渐变平.其中，最差充电站节点使用移除算子（SWR）、客户节点相关移除算子（CSRR）和充电站节点比较插入算子（SCI）、客户节点后悔插入算子（CRI）发现了最大数量的改进解，性能较好.

图4 每个算子的完成度Fig.4 Degree of completion for each operator

图6 插入算子的权重变化Fig.6 Weight change of insertion operator

图5和图6 分别展现了移除、插入算子在算法执行过程中的平均权重.尽管从百分比上讲，最大延误客户节点移除算子（CBLR）和客户节点相关移除算子（CSRR）发现了数量相近的改进解（图4），但从图5 可以看出随着时间的推移，该算子的性能会下降，并被客户节点相关移除算子（CSRR）所超越.从图6 可以看出，充电站节点比较插入算子（SCI）和客户节点后悔插入算子（CRI）的性能始终优于其他插入算子.

图5 移除算子的权重变化Fig.5 Weight change of remove operator

综合图4～图6 可知，针对CEVRPTW 问题所设计的ALNS 邻域算子中，移除算子表现较好的是客户节点相关移除（CSRR），最差充电站节点使用移除（SWR）次之；表现较差的是客户节点随机移除（CRR）、客户节点随机回路移除（CRRR）和客户节点区域移除（CZR）.插入算子表现较好的是充电站节点比较插入（SCI）和客户节点后悔插入（CRI）；表现较差的是客户节点随机插入（CI）和客户节点贪婪插入（CGI）.