带时间窗的多车型需求可拆分揽收配送问题

张 凯，靳 鹏，崔 勇

合肥工业大学 管理学院，合肥230009

传统物流企业利用自有车辆为客户提供揽收/配送服务。由于自有车辆具有使用成本高、管理难度大等缺点，在共享经济大背景下，物流企业开始采用一种基于外部车辆的新型揽收配送模式，利用外部车辆为客户提供揽收/配送服务，外部车辆主要是社会闲置车辆。新型揽收配送模式可以减少企业自有车辆的使用，降低物流成本，并对提高社会闲置车辆利用率和环境保护方面具有重要意义。

本文研究了基于新型揽收配送模式下的一种带时间窗的多车型需求可拆分揽收配送问题（Multi-Vehicle Split Pickup and Delivery Problem with Time Windows，MVSPDPTW），使用社会闲置车辆作为揽收/配送任务的资源，社会闲置车辆拥有自己的起点、终点以及最早出发时间和最晚返回时间，车辆必须在最早出发时间后从起点出发为客户提供揽收/配送服务，并在最晚返回时间前返回终点。在服务客户的过程中，车辆需要在客户时间窗内提供服务，并且允许需求拆分，同一客户需求可由多辆车完成，以便充分利用车辆的装载能力。

车辆路径问题（Vehicle Routing Problems，VRP）[1]自1959 年被提出以来，一直是物流研究领域的热点问题。目前，在VRP 的研究中一般限制客户的需求只能由一辆车来完成，属于需求不可拆分的VRP。文献[2]在1989年首次研究了需求可拆分的VRP（Split Delivery Vehicle Routing Problems，SDVRP），指出对客户需求进行合理的拆分有利于提高车辆装载率和节约总成本，并在1994 年证明SDVRP 是NP-hard[3]问题。目前已有不少学者从不同角度对SDVRP展开了研究。由于客户通常要求在给定的时间内得到服务，因此一种考虑时间窗约束的SDVRP[4-7]得到了研究，文献[4]首次将时间限制引入到SDVRP，构造了考虑时间窗约束的SDVRP模型，提出了一种两阶段启发式算法进行求解。文献[5]提出一种禁忌搜索算法求解考虑时间窗约束的SDVRP，设计了多种新型邻域搜索算子。文献[6-7]研究了考虑软时间窗约束的SDVRP，软时间窗约束通过添加惩罚函数允许客户在给定的时间窗外被服务。在SDVRP中只考虑了客户的配送需求，近些年随着逆向物流的发展，客户往往同时含有揽收/配送需求，因此一种考虑揽收需求的SDVRP[8-9]得到了研究，文献[8]首次将需求可拆分引入揽收配送问题，揭示了需求可拆分带来的好处，证明对于给定的一组揽收和配送位置，当客户需求量略大于车辆载重量一半时可以带来最大的收益。文献[9]提出了一个改进蚁群算法求解需求可拆分的揽收配送问题，实验结果表明该算法可以提高配送中心的运作效率。

以往对SDVRP 的研究基本建立在单一车型上，但是在实际配送中，车辆往往拥有多种车型[10-12]。文献[10]研究了多车型的SDVRP，车辆的装载能力与最大行驶里程各不相同，提出了一种模拟退火算法进行求解。文献[11]研究需求可拆分的多车型开放式车辆路径问题，车辆从不同的配送中心出发为客户服务，在服务完最后一个客户后结束服务而不需要返回配送中心，提出了一种禁忌搜索算法进行求解。文献[12]分别研究了单车型和多车型的SDVRP，通过实验结果说明了不同装载能力的车辆更符合实际需求并且可以有效提高提高车辆的装载率。

上述多数文献大都将揽收需求、配送需求、多车型、客户时间窗等单独研究，但是随着物流技术的快速发展，用户需求更加复杂，这种单一模式研究已经难以满足实际应用需求。因此，本文综合考虑上述多种要素，研究了一个带时间窗的多车型需求可拆分揽收配送问题，针对这个问题本文以执行任务车辆行驶路径总长度最小为目标函数，建立了一个混合整数线性规划模型，提出了一种高效禁忌模拟退火（Tabu Simulated Annealing，TSA）算法进行求解。

1 问题描述

在MVSPDPTW中存在给定的一个配送中心，一组地理位置不同的客户和一组多车型车辆，MVSPDPTW的目标是选择合适的车辆执行完所有客户任务，使得执行任务车辆行驶路径总长度最小。

0表示配送中心，T={1,2,…,n}表示客户点集合，每个客户点对应一个揽收/配送任务，将对应揽收任务的客户点称为pickup 点，对应配送任务的客户点称为delivery点。对于客户点i∈T，[ei,li] 表示i的时间窗，ei表示最早开始服务时间，li表示最迟开始服务时间。qi表示i需要揽收/配送的货物量。每个客户点可以被多个车辆访问，但只允许被同一车辆至多访问一次。K={1,2,…,m}表示车辆集合，对于k∈K，sk、nk、Qk、ok、fk分别表示车辆k的起点、终点、最大容量、最早出发时间以及最晚返回时间。如图1所示，车辆需要在最早出发时间后从起点出发，在pickup点揽收货物送回配送中心，在配送中心装载货物配送至delivery点，最终在最晚返回时间前返回终点。车辆执行任务的过程中允许多次访问配送中心。

图1 车辆执行任务示意图

2 数学模型

2.1 虚拟化处理

由于车辆可能多次访问配送中心，如果模型中以一个变量去区分车辆第几次访问配送中心，则建模难度和求解难度都会显著增加。为了降低难度，采用了以下处理方式：T={1,2,…,n}表示客户点集合，n为实际客户点的数量，对于编号为i的pickup 点，构建一个编号为n+i的虚拟delivery点与之对应；对于编号为j的delivery点，则重新编号为n+j，并构建一个编号为j的虚拟pickup 点与之对应。经过上述处理后，T=TP⋃TD，TP={1,2,…,n} 表示pickup 点集合，TD={n+1,n+2,…,2n}表示delivery点集合，车辆需要从点i(i∈TP)揽收货物配送至点n+i。该种处理方式将配送中心替换成了n个虚拟客户点，所有虚拟客户点的位置与配送中心位置相同且无服务时间窗限制；每个虚拟客户点揽收/配送的货物量与对应实际客户点相等。如图2 所示，图1中的配送中心经过虚拟化处理后替换成了方框内的虚拟客户点，其中i*=n+i。

图2 虚拟化示意图

2.2 模型符号

经过上述处理后，将MVSPDPTW在全有向图上图上定义G=(V,A),V={S⋃N⋃T⋃0}表示顶点集，A={(i,j):i,j∈V,i≠j}表示弧集。

S：车辆的起点集合，S={s1,s2,…,sm}；

N：车辆的终点集合，N={n1,n2,…,nm}；

T：客户点集合，T={1,2,…,2n}；

dij：点i和点j之间的距离，i,j∈V；

tij：点i和点j之间的时间，i,j∈V；

bik：车辆k到达点i的时间，k∈K,i∈T⋃0；

aik：车辆k在点i的开始服务时间，k∈K,i∈T⋃0；

rik：车辆k在点i的服务时长，k∈K,i∈T⋃0；

qik：车辆k在点i的装载/卸载的货物量，k∈K,i∈T⋃0；

wik：车辆k经过点i后承载货物容量，k∈K,i∈T⋃0；

xijk：表示车辆k是否直接从点i到点j，如果是则为1，否则为0；

M：表示非常大的整数。

2.3 线性模型

其中，公式（1）为目标函数，使得执行任务车辆行驶路径总长度最小；公式（2）、（3）表示车辆从起点出执行任务并最终返回终点；公式（4）、（5）表示每个客户点可以被多个车辆访问，但只允许被同一车辆至多访问一次；公式（6）、（7）、（8）表示车辆需要从pickup点揽收货物配送至对应的delivery点；公式（9）表示车辆从起点出发以及返回终点的承载货物容量均为0；公式（10）表示容量约束；公式（11）、（12）表示车辆在客户点承载货物容量变化；公式（13）表示每个客户点的分割货物量之和等于该客户点的总货物量；公式（14）、（15）表示车辆需要在最早出发时间后从起点出发为客户提供服务，并在最晚返回时间前返回终点；公式（16）、（17）表示车辆到达客户点时间和在客户点开始服务时间之间的关系；公式（18）表示车辆只能在客户时间窗内提供服务。

3 算法设计

为了解决MVSPDPTW，提出了一种禁忌模拟退火算法，设计了两种新的邻域搜索算子，分别用于修复违反车辆容量约束以及换车操作。在算法内加入禁忌机制以及违反约束惩罚机制，实现了搜索空间的有效裁剪，提高了算法的全局寻优能力。首先利用动态贪婪算法生成初始可行解，并将初始可行解作为禁忌模拟退火算法的输入，最终输出全局最好解。

3.1 初始解生成

模拟退火算法本质上是一种迭代求解算法，初始解决定了迭代的起点，因此一个好的初始解有助于算法找到较好的最终解，本文参考文献[5]中初始解的生成方法设计了一种简单有效的动态贪婪算法生成初始解。

公式（20）、（21）分别表示车辆以及客户点的选择准则，车辆的选择以服务时间最大化为原则。客户点的选择以当前规划路径中最后的访问对象i和待选择客户点j之间总的旅行时间和等待时间最小化为原则。利用公式（21）为当前车辆依次选择服务的客户点，检验所选客户点是否满足插入约束，包括时间窗约束、容量约束、客户点访问次数约束以及需求约束，其中需求约束是指所有配送的货物来源于配送中心，所有揽收的货物送回配送中心。如果满足约束则将客户点插入到当前车辆规划路径中，重复这个过程直到所有的客户需求被满足。动态贪婪算法的具体流程如下：

输入：车辆集合K、客户点集合T、配送中心0。

输出：初始可行解Sinitial。

步骤1 如果集合T不为空集，转步骤2，否则，算法结束输出结果。

步骤2 利用公式（20）从K中选出当前车辆k。

步骤3 构造待选择客户点集合C=T。

步骤4 如果C不为空集，转步骤5，否则，转步骤9。

步骤5 利用公式（21）从C中选出客户j。

步骤6 检验j是否满足插入约束，如果满足则将j插入到当前车辆规划路径Rk中。

步骤7 检验j的需求是否全部被满足，如果满足则从T中删除j。

步骤8 将j从C中删除，转步骤4。

步骤9 检验Rk最后一个访问对象是否为配送中心0，如果是，转步骤11，否则，转步骤10。

步骤10 检验0是否满足时间窗约束，如果满足将0插入到Rk中，转步骤3，否则，转步骤11。

步骤11 从K中删除k，转步骤1。

3.2 违反约束惩罚机制

在迭代过程中接受不可行解有助于增大邻域搜索范围以及提高算法跳出局部最优的能力，帮助算法在迭代步骤[13]更好的可行解，因此本文设计了一种违反约束惩罚机制。整体约束惩罚函数设计如公式（22），B表示路径规划解的适应度值，β表示参数惩罚权重，η(i,k)、φ(i,k)、λ(i,k)分别表示Rk中点i容量约束评价值、时间窗约束评价值以及需求约束评价值，它们的计算方式由公式（23）、（24）、（25）给出。Rk={ok,…,nk}表示车辆k的执行任务路径，它按照车辆的访问顺序排列而成。

3.3 禁忌机制

禁忌搜索是一种模拟人类思维过程的元启发式算法，它的基本思想是避免对已搜索过的空间进行重复搜索，迫使算法去探索其他新的解决方案。

将每次迭代的当前解[14]作为禁忌对象存储在禁忌表中，在以后的l（禁忌表的长度）次迭代中，该对象被设置成禁忌状态，这与在禁忌结构中常用的边移动禁忌[15-16]不同。禁忌表的长度指禁忌表储存对象的个数，当储存对象个数达到上限，再有对象加入时，则采用先进先出的原则将前面的对象从禁忌表中移除，再进行加入。

3.4 多邻域结构

多邻域结构可以增加解的多样性，提高算法找到全局最好解的可能性。本文采用了五种邻域搜索算子，其中迁移算子、迁移分割算子、交换算子来源于经典邻域搜索算子[17]的引用，修复算子和换车算子针对本文在迭代过程中接受不可行解以及多车型的特点而设计，分别用于修复违反容量约束和换车操作。

为了更好地解释算子操作，定义了以下符号，Rk(k∈K)和Rh(h∈K)分别表示车辆k和车辆h的执行任务路径，对于客户点分别表示在算子操作前后车辆k在客户i装载/卸载的货物量。公式（26）给出了B(Rk)的计算方式，B(Rk)表示路径Rk的适应度值。

任务集：在已规划好的车辆路径中，以车辆的起点、终点以及配送中心等非客户点为界，两个相邻的非客户点之间连续的客户点组成一个任务集。如果两个相邻的非客户点之间不包含客户点，则不构成任务集。如图3所示，客户点1和2组成一个任务集，终点与配送中心之间不含客户点则不构成任务集。

图3 任务集示意图

迁移算子：迁移算子用于将客户点从原路径中删除插入到其他路径中，如图4所示，对于客户点i∈Rk，将i从Rk中删除，以B(Rh)最小化为原则将i插入到Rh中。

图4 迁移算子操作示意图

迁移分割算子：如图5所示，客户i由车辆k和h共同服务，客户j由车辆h服务，i和j属于路径Rh中同一个任务集。将i从路径Rk中删除，然后以B(Rk)最小化原则将j插入到Rk中，此时i由车辆h服务，j由车辆k和h共同服务，k和h在i和j装载/卸载的货物量调整为。

图5 迁移分割算子操作示意图

交换算子：交换算子用于交换两条路径中的客户点，如图6 所示，客户点i∈Rk，客户点j∈Rh，将i和j进行交换，交换后，i∈Rh,j∈Rk。

图6 交换算子操作示意图

修复算子：修复算子用于修复违反车辆容量约束，如图7 所示，假定车辆k在服务客户i时承载货物容量超过了车辆的最大容量Qk，此时在i所属任务集中任选一个客户点j,j可以为i，任选一条非空路径Rh，将j以B(Rh)最小化原则插入到Rh中，k和h在j装载/卸载的货物量调整为：的计算方式由公式（23）给出。

图7 修复算子操作示意图

换车算子：对于非空路径Rk以及空路径Rh，将Rk中的所有客户点构成客户点集合C，然后按照初始解中路径的生成方式，将C中的客户点插入到Rh中。如果C中所有的客户点都成功插入到Rh中，那么路径Rk清空，否则的话将插入到Rh的客户点从Rk中删除。

图8 换车算子操作示意图

由于单个客户点只允许被同一车辆至多访问一次，所以在算子操作后，如果出现如图9所示情况，Rk中含有两个客户点i，则以B(Rk)最小化原则选择一个i从Rk中删除，将k两次服务i装载/卸载的货物量合并。

图9 合并客户点示意图

3.5 禁忌模拟退火算法描述

首先利用动态贪婪算法生成初始可行解，并将其设置为第0代的当前解Snow和全局最好解Sbest，然后开始迭代，在每次迭代中，利用设计的准则更新Snow和Sbest。最后当当前温度小于设置的最低温度时，算法终止输出Sbest，具体算法流程如下：

输入：初始温度T0、终止温度Tend、冷却速率RT、迭代步长L、邻域解数量Nu、初始可行解Sinitial。

输出：全局最好解Sbest。

步骤1 将Sinitial设置为Snow以及Sbest。

步骤2 初始化当前温度T'=T0，初始化当前温度下迭代次数L'=0。

步骤3 初始化禁忌表，将Snow加入到禁忌表中。

步骤4 如果T'>Tend，转步骤5，否则，算法结束输出结果。

步骤5 如果L'

步骤6 从五个邻域算子中随机选择一个邻域算子转换当前解生成邻域解。

步骤7 重复步骤6直到邻域解的数量达到Nu。

步骤8 从Nu个邻域解中选出最优的没有被禁忌的邻域解Scan。

步骤9 如果Scan

步骤10 如果>Random(0,1),Scan将成为新的Snow,Random(0,1)用于生成（0，1）上均匀分布的随机数据，否则保留Snow。

步骤11 如果Scan为可行解并且满足Scan

步骤12 如果Snow发生变化，则将Snow加入禁忌表。

步骤13 转步骤5。

4 仿真算例

本文利用文献[18]使用的考虑时间窗约束的VRP（Vehicle Routing Problem with Time Windows，VRPTW）数据集、文献[19]使用的考虑时间窗约束的VPRSD（Split Delivery Vehicle Routing Problem with Time Windows，SDVRPTW）数据集以及构造的MVSPDPTW 仿真数据集对TSA 算法进行了测试并对实验结果进行对比分析，每个算例都运行20 次。所有实验均在i5-8500 CPU 3.0 GHz、内存16 GB 的台式计算机上、myeclipse2017的环境下进行。

4.1 参数设置方法

通过设置合理的参数可以有效提高算法的性能[20]，首先参考文献[21-22]中的参数设置确定本文各个参数的测试区间，然后选取不同规模算例，通过更改参数组合进行多组测试，最后确定参数的最佳组合。

各参数取值范围以及步长具体参考表1，为简化计算复杂度，本文采用控制变量法依次确定各参数的最佳值。TSA 参数设置的最终结果如下所示：初始温度T0=100、终止温度Tend=0.1、迭代长度L=30、降温速率RT=0.93、惩罚权重β=4、邻域解数量Nu=2n+150（n为客户数量）、禁忌表的长度l=8。

表1 算法参数组合

4.2 基于VRPTW、SDVRPTW数据集的数值实验

由于目前在MVSPDPTW 领域没有公认的数据集可供使用，导致没有办法对本文算法直接进行测试，MVSPDPTW 是经典VRPTW 以及SDVRPTW 的拓展，问题中很多约束条件相似，因此如果本文算法可以很好地解决VRPTW 以及SDVRPTW，那么有理由相信本文算法也可以有效地求解MVSPDPTW。

4.2.1 基于VRPTW数据集的数值实验

从文献[18]使用的VRPTW数据集中有针对性地选择了27 个VRPTW 算例进行测试，VRPTW 算例中车辆容量皆为200。表2记录了TSA算法求解每一个算例结果Rbest和算例已知最好解Kbest对比，Gap表示Rbest与Kbest之间差距，Gap=(Rbest-Kbest)/Kbest。Kbest数据来源于http：//web.cba.neu.edu/～msolomon/problems.htm。

表2 TSA算法求解结果与已知最好解对比

根据表2中数据可知，TSA算法更新了4个VRPTW算例（算例编号C101，客户数量25；算例编号C102，客户数量50；算例编号C107，客户数量50；算例编号C109，客户数量50）的已知最好解。在剩余23个VRPTW算例中，TSA 算法的求解结果虽然劣于已知最好解，但是它们之间的差距均保持在1%之内。实验结果说明了TSA算法求解VRPTW效果好。

4.2.2 基于SDVRPTW数据集的数值实验

从文献[19]使用的SDVRPTW 数据集中有针对性地选择了12 个SDVRPTW 算例进行测试，SDVRPTW算例中车辆容量皆为100。表3记录了TSA算法和分支与价格和削减[19（]Branch and Price and Cut algorithm，BPC）算法求解12 个SDVRPTW 算例结果对比。Gap表示TSA 算法求解每一个算例结果(Rbest)和BPC 算法求解每一个算例结果(Cbest)之间的差距，Gap=(Rbest-Cbest)/Cbest。

从表3中数据可知，TSA算法在客户数量为25的小规模算例下的求解结果与BPC算法相同。随着客户数量的增大，虽然TSA算法的求解结果劣于BPC算法，但是它们之间的差距均保持1%以内。在求解时间方面，在小规模算例下，BPC算法求解时间小于TSA算法，但是随着客户数量的增大，BPC算法求解时间远超TSA算法，并且在相同客户数量下不同算例的求解时间上，TSA 算法基本保持不变而BPC 算法则差异明显。实验结果充分说明了TSA 算法求解SDVRPTW 效果接近BPC算法。

表3 TSA算法与BPC算法求解结果对比

4.3 基于仿真数据集的数值实验

基于4.2.2 节的12 个SDVRPTW 算例构建12 个MVSPSPTW 仿真算例，不改变SDVRPTW 算例中配送中心位置、客户位置、客户货物量以及客户时间窗，仅随机选择半数客户将其配送任务改为揽收任务，对于车辆的起点和终点不再设置为配送中心，而是在客户所属100×100的正方形区域内随机生成；车辆的开始服务时间为0～800间随机生成的整数，结束服务时间等于500～900间随机生成的整数加上车辆的开始服务时间；车辆最大容量为70～120 间随机生成的整数。此外，当客户数量为25时随机生成15辆车，当客户数量为50时随机生成30辆车，保证车辆资源充足。

表4 记录了TSA 算法、传统模拟退火（Simulated Annealing，SA）算法、粒子群算法（Particles Swarm Optimization，PSO）求解MVSPDPTW 仿真算例结果对比，PSO算法的设计参考了文献[23-24]中的处理方式，传统SA 算法与PSO 算法的参数通过4.1 节参数设置方法来确定。Gap1、Gap2 分别表示TSA算法求解每一个算例结果(Rbest)与SA算法求解每一个算例结果(Tbest)、PSO算法求解每一个算例结果(Pbest)之间的差距，Gap1=(Tbest-Rbest)/Rbest，Gap2=(Pbest-Rbest)/Rbest。

根据表4中数据可知，本文所提出的TSA算法对比传统SA 算法和PSO 算法求解性能更优。主要体现在：（1）TSA算法的求解结果最优，相比传统SA算法最高提升了2.40%，平均提升了1.14%；相比于PSO算法最高提升了5.19%，平均提升了3.52%。（2）在求解时间方面，SA 算法最长，TSA 算法次之，PSO 算法最短，但TSA 算法与PSO算法耗时总体差距不大。此外，TSA算法求解相同客户数量算例的时间基本保持不变，具有较好的稳定性。TSA 算法求解性能更优的原因主要在于：（1）由于本文多车型的特点，因此选车结果会对最终解的优劣产生至关重要的影响，换车算子有利于寻找到最合适的车辆执行任务。除此之外多种算子配合的方式有利于扩大邻域搜索范围，提高算法的全局寻优能力。（2）TSA算法内加入了禁忌机制以及违反约束惩罚机制，实现了搜索空间的有效裁剪，避免算法陷入局部最优的同时节约了求解时间，相比之下传统SA 算法和PSO 算法易陷入局部最优。

表4 TSA、SA、PSO算法求解结果对比

5 总结展望

本文研究了基于新型揽收配送模式下的一种带时间窗的多车型需求可拆分揽收配送问题。针对这个问题建立一个混合整数线性规划模型。在构造模型的过程中，通过为每个客户点构建一个对应虚拟揽收/配送点的方式减少模型中部分变量的增加，并使用大M法对非线性约束做了线性化处理。设计了一种高效TSA 算法，TSA 算法在传统SA 算法框架内加入了禁忌机制以及违反约束惩罚机制，实现了搜索空间的有效裁剪，提高了算法的全局寻优能力。TSA 算法更新了27 个VRPTW算例中4个算例的已知最好解；12个SDVRPTW算例的求解结果中，TSA算法与BPC算法的差距均不超过1%，且求解时间更短；就12个MVSPDPTW仿真算例而言，TSA 求解结果均不差于传统SA 算法以及PSO 算法，且分别最高提升了2.40%和5.19%。实验结果充分证明了本文所提出的TSA 算法求解MVSPDPTW 的可行性和有效性。

本文所研究的内容依旧有许多亟待解决的问题，将来的研究工作可以从以下两个方面展开：（1）考虑需求依订单拆分而不是任意连续拆分；（2）考虑使用外部车辆和自有车辆共同配送方式。