随机时变下带时间窗的取送货车辆路径问题优化研究

靳 鹏，张歆悦 （合肥工业大学 管理学院，安徽 合肥 230009）

0 引 言

车辆路径问题（Vehicle Routing Problem，VRP） 自提出以来，一直是物流组织优化领域中的研究难点。在实际物流系统中，为满足取货或配送及服务时间窗限制的客户需求，带时间窗的取送货车辆路径问题（Pickup and Delivery with Time Windows，PDPTW） 逐渐受到众多学者的关注。然而，近年来城市交通日趋拥堵，受交通管制、交通事故等不确定因素的影响，城市交通路网存在一定的时变性和随机性，车辆的行驶速度是时变的。在此背景下，研究随机时变下带时间窗的取送货车辆路径问题（Stochastic Time-dependent Pickup and Delivery with Time Windows, STDPDPTW） 具有重要意义。

带时间窗的取送货车辆路径问题（PDPTW） 是经典车辆路径问题（Vehicle Routing Problem, VRP） 的扩展，属于NP-hard问题，最早由Wilson提出。Bent 等人首次提出两阶段混合算法，第一阶段选用模拟退火算法优化路径数量，第二阶段选用大邻域搜索算法降低路径成本，有效求解PDPTW 问题。Al Chami 等人使用贪婪随机的自适应搜索方法生成初始解，后使用遗传算法优化解决方案，有效提高求解PDPTW 问题的效率。目前用于求解带时间窗的取送货车辆路径问题的启发式算法主要包括遗传算法、模拟退火算法、自适应大邻域搜索算法、蚁群算法等。

然而，城市化进程逐步加快，道路交通超负荷运转，加之交通管制和交通事故等因素的影响，在物流运输过程中，考虑车辆行驶时间的时变性和随机性更贴近实际。随机时变下车辆路径问题（Stochastic Time-dependent Vehicle Routing Problem,STDVRP） 最早由随机时变最优路径问题（Stochastic Time-dependent Optimal Path Problem, STDOPP） 发展而来，Hang 等以行驶时间负效用函数计算随机时变路径行驶时间，随机路网中的最佳路径受随机依赖性影响。张春苗指出以车辆行驶时间的概率分布为基础，网络拓扑结构复杂，算法求解时间复杂度随路网规模增大呈指数增长。由Sim提出鲁棒优化模型，将随机时变路网转化为确定型时变，曹慧等将该方法应用于随机时变路网的最优路径问题，结果证明满足随机一致性条件，车辆行驶时间可依据确定型时变路网计算。Duan也指出鲁棒优化时间模型适用于实际的道路网络。段征宇等人建立STDVRP 鲁棒优化模型，考虑客户的时间窗约束，增大了求解的复杂度，并设计改进蚁群算法进行有效求解。

综上，多数文献独立研究STDVRP 和PDPTW，同时考虑随机时变、时间窗约束和取送货的车辆路径问题的研究较少，本文建立了随机时变下带时间窗的取送货车辆路径问题（STDPDPTW） 模型，基于鲁棒优化方法计算随机时变路网下车辆行驶时间，为快速且按时满足客户需求，以车辆总行驶时间最小化为目标函数建立混合整数规划模型，并设计一种混合遗传模拟退火算法进行求解。

1 问题描述与数学模型

1.1 问题描述。STDPDPTW 研究一定数量的车辆从地理位置不同的车场出发在客户能接受的时间范围内为客户提供取货或配送服务，受道路类型和行驶时段的影响，同一路段可能具有不同的行驶速度，制定总行驶时间最小的路径规划方案。对研究问题做如下描述：（1） 配送网络中有多个车场，每个车场有一定数量相同类型的运输车辆，每辆车的起止点为同一车场；（2）每个客户仅由一辆车进行服务；（3） 客户点的需求量不大于车辆载重量；（4） 每个客户点有服务时间窗的要求，车辆可早于最早开始服务时间到达客户点，等待直至最早开始服务时间进行服务，但不能晚于最迟开始时间到达。

1.3 问题建模。以服务完所有客户所需的总行驶时间最小化为目标函数，则随机时变下带时间窗的取送货车辆路径问题的数学模型建立如下：

约束条件如下：

式（1） 表示最小化车辆总行驶时间；式（2） 表示每个客户点仅由一辆车进行访问；式（3） 表示每个客户点的流量平衡约束；式（4） 表示每辆车仅使用一次；式（5）、式（6） 表示车辆到达客户点的时间和客户点的开始服务时间之间的关系；式（7） 表示一个订单包含的取货和送货请求均由同一辆车提供服务；式（8） 表示车辆从车场出发完成任务后并返回该车场；式（9） 表示完成客户点i 的取货时间早于完成客户点n+m+i 的送货时间；式（10） 表示车辆从车场出发以及返回车场承载货物量均为0；式（11） 表示车辆在执行任务过程中承载货物容量不能超过车辆容量限制；式（12）、式（13） 表示保证车辆离开上一节点到下一节点之间容量的一致性，车辆在客户点之间承载容量的变化；式（14） 表示决策变量取值约束。

1.4 时变路网下车辆行驶时间计算方法。在实际交通中，车辆在不同时间段内道路的拥堵情况不同，因而车辆的行驶速度存在差异，考虑实际生活中城市路段早晚高峰的情况，车辆的行驶速度不是阶跃变化，而是更近似平滑变化。将一天的工作时间等分为η 个分段，改进Figliozzi的时间依赖函数，使用平滑变化的速度时间函数，如图1 表示车辆的行驶速度时间函数。

图1 速度时间函数

输入：车辆k 离开节点i 的出发时刻λ，时间分段τ，时段范围(ts,te]，节点i 与节点j 的距离d。

1.5 随机时变车辆行驶时间的鲁棒优化。随机时变路网的车辆行驶时间可以通过路段概率分布函数以及期望估算模型进行优化，然而车辆必须在客户规定的时间窗内进行服务，完成所有客户的需求，而不是以一定概率进行满足。因此，选用最大最小准则作为车辆行驶时间的鲁棒优化方法，对于车辆k 以时刻λ从节点i 到节点j 的最坏情况发生，从最坏的情况中找到最优路径的行驶时间为：

2 遗传模拟退火算法设计

遗传算法通过适应度函数对搜索空间中的多个解进行评估，能够快速搜索到新解，适用于本文组合优化问题的求解，但后期搜索能力减弱致使该算法易过早收敛。本文针对问题特点和遗传算法的不足，结合模拟退火算法能有效跳出局部循环的特点，提出了一种高效求解STDPDPTW 模型的混合遗传模拟退火算法（Hybrid Genetic Simulated Annealing Algorithm, HGSA）。

2.1 染色体编码方式。本文设计了一种三行染色体编码方式，染色体的第一行为车辆访问的客户点，第二行为与第一行对应的车辆编号，第三行为客户点的需求类型，1 为取货点，-1 为送货点，0 为车场。每条染色体的编码方式如图2 所示，以8 个客户点为例，车辆1 的路径为D-8-3-1-5-D，车辆2 的路径为D-2-4-7-6-D。

图2 染色体编码示意图

2.2 初始种群的生成。根据客户点的需求类型和车辆的容量约束，按照以下5 个步骤生成初始种群：

Step1：从客户点集合P 中随机选择一个取货点和一个送货点进行两两组合形成待插入客户点序列P；

Step2：根据车辆容量约束，按待插入客户点序列P顺序为车辆划分客户点，得到车辆所访问的客户点序列H，形成染色体的第一行编码；

Step3：在每辆车所访问的客户点序列H的最前面和最后面加入该车辆对应的车场编号，用来表示车辆出发的起点和返回的终点，得到每辆车辆的初始路径访问的客户点序列P，形成染色体的第一行编码；

Step4：在每辆车的初始路径的编码的第二行添加对应的车辆编号，在第三行添加客户点对应的任务类型编号，形成一条完整的染色体，如图3 所示；

图3 种群初始化过程示意图

Step5：根据预设的种群规模N重复Step1～Step4，得到初始种群。

2.3 多段多点交叉算子。多点交叉操作不仅增加染色体的多样性，而且能够有效避免个体过早收敛，根据染色体编码方式的特点，染色体中的每一段代表一辆车访问的客户点序列，本文设计了多段多点交叉算子，算法步骤如下：

Step1：采用轮盘赌方式从种群中选择两条染色体作为父代染色体F和F；

Step3：分别从父代染色体F和F中相同车辆编号的分段F和F中随机选择连续的多个基因进行交换；

图4 染色体交叉操作示意图

2.4 修复算子。使用多段多点交叉方式后的染色体可能存在部分客户点缺失和重复的问题，则需要将未访问的客户点插入到路径规划方案中，且为缩短车辆的行驶时间，加快全局寻优过程，设计了修复算子，算法步骤如下：

Step1：从交叉操作后的染色体中标记出所有重复访问的客户点，利用公式（16） 依次计算染色体中重复点位上一位客户点n到达该客户点n的行驶时间与等待时间之和st；

Step2：重复访问客户点中相同编号的客户点两两比较st的大小，保留st较小的客户点，并去除标记；

Step3：利用公式（17） 从未访问客户点集合N中选择客户点依次替换标记出的重复访问客户点，生成修复后的新子代个体。

式（16） 表示染色体中重复点位上一位客户点n到达该客户点n的行驶时间与等待时间之和。式（17） 表示从未访问客户点集合N中选择待插入客户点n，以依次待替换的染色体中重复点位的上一位客户点n到达未插入的客户点n之间的行驶时间和等待时间最小化为替换客户点的选择准则。

2.5 变异算子。变异操作有助于提升该算法的全局搜索能力，避免陷入局部最优的困境。本文选择两点交换变异方式，随机从染色体中选择两个客户点进行位置交换，得到新的子代染色体，具体操作如图5 所示。

图5 染色体变异操作示意图

2.6 种群更新。通过选择部分父代种群中的精英染色体替代新产生的子代种群中相同数量的染色体来完成种群更新操作，选择替代染色体的数量为：N=N× 1-Ga（ ）p ，其中N为种群规模，Gap 为代沟。具体操作如下：一个新的种群由父代种群中适应度值前N位的染色体和子代种群中前（N-N）位的染色体组成，作为遗传算法下一轮迭代的新种群。

2.7 邻域结构。邻域搜索操作可以扩大搜索范围增加解的多样性，避免算法陷入局部最优，增大算法找到全局最优解的可能性。本文采用了四种邻域搜索算子，包括逆序搜索算子、单点交换算子、合并搜索算子和两点交换算子，邻域算子操作如下：（1） 逆序搜索算子：随机选择两个客户点进行交换，并将两个客户点中的客户点反转逆序排列，如图6（a） 所示。（2） 单点交换算子：从原路径中随机选择客户点移除，插入到其它路径中，如图6（b） 所示。（3） 合并搜索算子：随机选择一条路径，将该路径上的客户点插入到其它路径中，如图6（c） 所示。（4） 两点交换算子：从两条不同的路径中分别各选择一个客户点，并交换他们的位置，如图6（d） 所示。

图6 领域算子示意图

2.8 混合遗传模拟退火算法求解步骤。混合遗传模拟退火算法将遗传算法与模拟退火算法结合，通过遗传算法得到较优解，以此作为模拟退火算法的初始解，再通过邻域搜索最终得到近似最优路径方案。具体算法步骤如下：

Step1：初始化混合遗传模拟退火算法的参数，种群规模N、代沟Gap、交叉概率P、变异概率P、初始温度T、终止温度T、步长L、冷却速率R；

Step2：种群初始化，获得车辆初始路径序列方案集合；

Step3：通过式（1） 计算染色体的适应度；

Step4：根据交叉概率P对选择的两条父代染色体进行交叉操作生成交叉后的子代染色体；

Step5：对交叉后的染色体使用修复算子进行修复，补全因交叉过程缺失的客户点；

Step6：根据变异概率P对交叉后的子代染色体进行变异操作生成变异后的子代染色体；

Step7：种群更新操作；

Step8：判断种群中最优解的适应度值连续10 代是否下降，若没有下降，则输出遗传算法的当前最优解S，转Step9，否则，转Step4；

Step9：初始化当前温度T=T，当前温度下迭代次数L=0，将Step8 输出的当前最优解S 作为当前解S；

Step10：如果温度T＞T，转Step11，否则，输出算法终止输出最优解；

Step11：如果迭代次数L＜L，则L=L+1，转12，否则，L=0，T=T×R，转Step10；

3 实验分析

本文所有实验使用Java 语言进行编写，eclipse2018-12 的编程环境，Windows 10 操作系统，Intel Core i7-7700 CPU@3.60 GHz 16.0 GB 的运行环境。

3.1 参数设置。通过合理的参数设置能够有效提高算法的求解质量，考虑到实际应用需求，选取多个不同算例，在合理的时间范围内，采用控制变量法，最终确定一个最佳的参数组合。首先根据经验确定参数的测试区间，再通过更改参数组合进行多组数值实验，记录实验结果和CPU 运行时间，最终确定HGSA 参数设置最佳组合如表1 所示。

表1 算法的最佳参数值

3.2 PDPTW 算例。目前关于STDPDPTW 研究所进行的实验并没有公认的数据集，且STDPDPTW 是在经典的PDPTW 问题中加部分约束形成的衍生问题。因此，选取Li & Lim PDPTW 标准数据集pdp_100的56 个算例进行数值实验，验证本文算法有效性。

本实验选用Li & Lim PDPTW 标准数据集中的56 个算例实验，算例包括LC 数据集（地点规则分布）、LR 数据集（地随机分布） 和LRC 数据集（地点混合分布）。所有车辆行驶道路类型相同，且以速度为1 的恒定速度行驶，以车辆总行驶距为优化目标。表2 记录了PDPTW 的已知最好解R、HGSA 算法的总距离S及其运算时间，Gap 表示R与S之间的差百分比，Gap=（R-S）/S。增的 点离值

表2 PDPTW 数据集实验结果

由表2 中数据可知，HGSA 算法更新了5 个PDPTW 算例分别为lc109、lc203、lr108、lrc106 和lrc204 的已知最好解。在其它的51 个PDPTW 算例中，HGSA 算法的求解结果劣于算例已知最好解，但与已知最好解的差距均值在1%之内，实验数据说明了HGSA 算法求解PDPTW 的有效性。

3.3 STDPDPTW 算例。基于Li & Lim 基准测试实例，根据实例规模将单一车场扩充为5 个车场，每个车场位置随机生成，且每个车场分配5 辆车用于完成任务，并增加路段中车辆行驶速度的波动范围，构建STDPDPTW 测试实例。依据城市道路通行情况，将道路分为五种类型，总服务时域分为9 个等分的时间分段，依据图1 将车辆的行驶速度时间函数设置为连续函数，不同的时段车辆的最高行驶速度和最低行驶速度不同，路段的道路类型由式（18） 判断：

式（18） 中，mod 为取余计算，例如从7 号节点到9 号节点的道路类型为2。

表3 表示不同道路类型在9 个等分的时间分段中车辆行驶速度的变化范围。对各道路类型和时间分段设置高、中、低三种波动范围，对于时间分段τ、τ和τ中行驶速度增加±15%的波动范围，对于时间分段τ、τ和τ中行驶速度增加±20%的波动范围，对于时间分段τ、τ和τ中行驶速度增加±10%的波动范围。根据行驶速度波动范围通过随机时变车辆行驶时间的鲁棒优化方法将随机时变路网转化为确定型路网。

表3 各道路类型车辆行驶速度变化范围

表4 记录了随机选取的10 个STDPDPTW 算例的实验结果，在相同实验条件下每个算例运行20 次。其中，20 次运行结果的最大值记录在C列中，平均值记录在C列中，最小值记录在C列中，20 次运行的平均运行时间记录在CPU列中，Gap表示C与C之间的差值百分比，Gap表示C与C之间的差值百分比。其中，Gap=（C-C）/C，Gap=（C-C）/C。

表4 STDPDPTW 数据集实验结果

由表4 数据可知，考虑车辆行驶时间的时变性和随机性，在车辆数量有限的情况下，通过优化车辆的出发时间和出发的车场位置可以缩短完成配送任务的时间。10 个STPDPTW 算例均在较短的时间内获得近似最优解，将引入多段多点交叉算子和修复算子的遗传算法与模拟退火算法结合，有效提高了局部寻优能力和收敛速度。另外，表4 的C与C之间的差值百分比的平均值仅相差2.01%，C与C之间的差值百分比的平均值仅相差1.68%，表明HGSA 具有较强的稳定性和寻优能力。

4 结 论

本文考虑城市交通拥堵、交通管制等不确定因素的影响，建立了随机时变下带时间窗的取送货车辆路径问题模型，提出考虑车辆加速度的行驶时间计算方式，运用随机时变车辆行驶时间的鲁棒优化方法将随机性时变转化为确定型时变，根据STDPDPTW 模型特点，设计了两阶段的混合遗传模拟退火算法，引入多段多点交叉算子和修复算子，保证了前期迭代的收敛速度后使用多领域搜索算子，提高了算法的搜索效率。通过标准数据集进行验证，HGSA 算法更新了56 个PDPTW 算例中5 个已知最好解，其余算例与已知最好解的差距均值在1%之内，同时STDPDPTW 算例表明，HGSA 算法求解STDPDPTW 具有较强的稳定性，能有效缩短车辆的总行驶时间，降低物流配送成本。