改进的ALNS算法求解MDVRPTW问题

Improved ALNS algorithm for solving multi-depots vehicle routing problem with time windows

摘要

研究了带时间窗多车场车辆路径问题（multi-depots vehicle routing problem with time windows，MDVRPTW），建立MDVRPTW模型，设计了结合混合高斯模型（Gaussian mixture model，GMM）聚类算法的自适应大邻域搜索（adaptive large neighborhood search，ALNS）算法。通过在邻域变换前将客户集进行分类，优化初始解，提高算法运算效率。算法使用6种不同变换因子，采用得分系统对变换因子进行评价，使算法能够在迭代的不同阶段自适应地选择合适的变换因子。分析了参数设置值的合理性，设计了3组仿真实验，实验结果验证了算法的高效性。

关键词

多车场车辆路径问题 / 时间窗 / 改进的ALNS算法 / GMM聚类算法 / 邻域变换

Key words

multi-depots vehicle routing problem / time windows / improved ALNS algorithm / GMM clustering algorithm / neighborhood transformation

多车场车辆路径问题（multi-depots vehicle routing problem，MDVRP）作为车辆路径问题（vehicle routing problem，VRP）衍生类的一种，相较VRP更接近于现实情形。Afshar等^［1］研究了在异构车辆及时间依赖下的MDVRP问题。Kubil等^［2］使用多目标蚁群优化算法求解了MDVRP问题。Reza^［3］则同时讨论了在MDVRP问题下的封闭式车辆路径问题及开放式车辆路径问题。在求解方式上，MDVRP存在两种基本求解模式，一种是先设立虚拟中心，将多车场问题转化成单车场问题，而后在求解过程中对每条子路径匹配最近的车场中心^［4］；而另一种是将客户进行分类，而后对每个客户子群体匹配合适的车场中心，将单个多车场问题转化成多个单车场问题进行求解^［5］。

在实际运送过程中，客户时间窗（time windows，TW）又分为硬时间窗^［6］和软时间窗^［7］。硬时间窗指每个客户的服务时间不能超出客户提供的可服务时间范围；软时间窗则可以超出其范围，但模型中会设置客户满意度惩罚项进行约束，在尽量满足客户需求的同时尽可能减少运输成本。Goeke等^［8］研究包含时间窗及充电站的VRP模型，通过改进的ALNS算法得到较好结果。

目前文献对MDVRPTW问题的研究较少。Fan等^［9］设计了混合遗传大邻域搜索算法（hybrid genetic algorithm with variable neighborhood search considering the temporal-spatial distance，HGAVNS_TS）来求解MDVRPTW模型，在初始化阶段设计了时空间距离计算方法用以将客户集合进行分类。吴凡等^［10］将MDVRPTW模型应用在应急物资运输上，提出一种变长基因型混合小生境遗传算法，运用模拟退火（simulated annealing，SA）算法的Metropolis准则和路径编码的海明距离，在每次迭代时保留较为优秀且结构互不相同的个体，保持较好的种群多样性，以避免遗传算法（genetic algorithm，GA）的早熟现象。

ALNS作为一种非常便捷的离散模型处理算法框架，易与其他算法融合，如将大邻域搜索（large neighborhood search，LNS）算法与GA算法融合^［9］，或将其与禁忌搜索（tabu search，TS）算法结合^［11］。Wu等^［12］用迭代自适应大邻域搜索算法求解异构车辆模型（heterogeneous fixed fleet vehicle routing problem，HFFVRP），并利用360个客户数据的算例得到了与BKs相近的解。

本文用先聚类后求解的方法对MDVRPTW进行计算，使用混合高斯模型将客户分类，再将MDVRPTW模型转换成多个VRPTW模型求解初始种群，优化种群结构，达到优化算法求解效率的目的。

1 MDVRPTW数学模型

设无向图 G = \left(V,E\right)， V = V_C \bigcup V_D为客户集合 V_C = \left\{\mathrm{1},\dots ,N\right\}和车场集合 V_D = \left\{\mathrm{1},\dots ,M\right\}的总和，N为客户总数，M为车场总数。 E = \left\{x_{ij}^k|i\ne j,i,j\in V,k\in K\right\}为可行边集，其中 x_{ij}^k等于1，表示车辆k从客户i到客户j，若为0则表示车辆k未从客户i到客户j，令 \begin{array}{ccc}{c}_{ij}& (i\ne j,& i,j\in V\end{array}）为点i到点j的距离。本文考虑的是同种车辆进行运输，车辆总数为K，车辆载荷为Q，车辆速度为v，单位行驶成本为P。

设客户i要求货运公司尽量在时间段［l_i， u_i ］内服务，为节省运输成本，存在部分客户无法完全按照客户时间段配送的情况，因此需设置客户满意度惩罚函数F_i。令 {\theta }_i为客户服务时间的惩罚因子，设车辆到达客户i的时间点为t_i，则存在以下3种可能：

（1） ，则令车辆等待，直到客户可服务时进行服务， \begin{array}{cc}{t}_i=l_i,F_i=\mathrm{0};& \end{array}

（2） l_i \le t_i \le u_i，则车辆可以直接服务客户，无须等待， F_i = \mathrm{0}；

（3） ，则车辆可以服务客户，但需要承担惩罚， F_i = {\theta }_i \left(t_i-u_i\right)。

MDVRPTW模型的基本假设为：客户完全得到服务；所有客户都在时间窗范围内被服务；车辆为同种车辆；车辆的各项属性已知；车辆从车场出发，在服务完所有安排好的客户后返回原车场；车辆在服务客户时的任意时刻，所载货物质量不能超过其载荷；客户仅被一辆运输车服务一次；车辆不允许重复经过某一个客户。

令q_i 为客户需求，v为车辆速度，c_ij 为从客户i到客户j的距离，t_i 为车辆到达客户i的时间点。tp_ij 为从客户i到j的时间，h_i 是车辆为客户i的服务时间， x_{ij}^k∈［0，1］，当车辆k从客户i到j时为1，否则为0。 y_i^k∈［0，1］，当车辆k经过客户i时为1，否则为0。F_i 是客户的满意度惩罚函数，L是混合高斯聚类每个高斯模型的个数。

MDVRPTW模型为

\mathrm{m} \mathrm{i} \mathrm{n} : F = p {\displaystyle \sum _{k=\mathrm{1}}^K} \sum _{i\in V}\sum _{j\in V}{c}_{ij}{x}_{ij}^k + \sum _{i\in V_C} F_i （1）

s.t.

\begin{array}{cc}\sum _{k\in K}\sum _{j\in V}{x}_{ij}^k=\sum _{k\in K}\sum _{j\in V}{x}_{ji}^k=\mathrm{1},& \forall i\in V_C\end{array} （2）

\begin{array}{ccc}\sum _{j\in V_C}{x}_{dj}^k=\sum _{i\in V_C}{x}_{id}^k\le \mathrm{1},& \forall k\in K,& \forall d\in V_D\end{array} （3）

\begin{array}{cc}\sum _{d\in V_D}\sum _{j\in V_C}{x}_{dj}^k\le \mathrm{1},& \forall k\in K\end{array} （4）

\begin{array}{cc}\begin{array}{cc}\sum _{i\in S}\sum _{j\in S}{x}_{ij}^k\le \left|S\right|-\mathrm{1},& \forall S\subseteq V\end{array},& \forall k\in K\end{array} （5）

\begin{array}{cc}\sum _{i\in V_D}\sum _{j\in V_D}{x}_{ij}^k=\mathrm{0},& \forall k\in K\end{array} （6）

\begin{array}{ccc}{t}_i^k+h_i+t{p}_{ij}^k\le t_j^k,& \forall k\in K,& \forall i,j\in V\end{array} （7）

\begin{array}{cc}\sum _{i\in V_C}{q}_i{y}_i^k\le Q,& \forall k\in K\end{array} （8）

\begin{array}{ccc}{x}_{ij}^k,y_i^k\in \left\{\mathrm{0,1}\right\},& \forall k\in K,& \forall i,j\in V\end{array} （9）

（10）

F_i = \left\{\begin{array}{ll}{\theta }_i\left(t_i-u_i\right),& t_i>u_i\\ \mathrm{0},& \mathrm{其他}\end{array}\right. （11）

式（1）为目标函数，即最小化车辆行驶成本与客户满意度成本之和；式（2）表示各客户被服务次数有且仅有一次；式（3）表示车辆k被使用时必须从车场出发并最终返回原车场；式（4）表示每辆车最多只能使用一次；式（5）表示车辆路径中不能产生子循环；式（6）表示车辆不能直接从一个车场行驶到另一个车场；式（7）表示车辆开始服务客户j的时间必须比到达客户j的时间要晚；式（8）表示车辆不能超载；式（9）中 x_{ij}^k \mathrm{、} y_i^k是布尔变量，分别表示当车辆经过路径ij或经过客户i时为1，否则为0；式（10）、式（11）指代客户i的时间及客户i的客户满意度惩罚函数。

2 算法设计

在多车场问题中有两种初始化方式，第一种是使用合适的聚类算法，将客户依照单个或多个属性^［9，13］分类，再依据类别初始化；另一种是生成虚拟中心^［4］，依照单车场模型初始化后再对每条路径分配最近的车场。本文使用第一种处理方式进行求解。

2.1 初始化

Wang等^［13］在研究MDVRPTW模型初始化时，将客户时间窗转换成空间属性，构建了三维k-means算法来分类客户，证明经过聚类后的客户集合生成的解较未聚类的客户集合生成的解更好。鲁建厦等^［14］设计了一种类似三标准^［15］的聚类算法，将供给成本与客户到每个聚类的距离相结合，以此计算每个客户计入不同子类的成本。Fernando^［5］依据各客户到各中心的距离对客户群体进行分类，证明了用聚类算法的初始解比直接初始化的初始解更优秀。

为寻找客户密度属性对车辆路径成本的影响，本文使用基于混合高斯模型的聚类算法进行求解。首先，假设所有客户点呈高斯分布，这使客户分类具有更多的可能性。使用均值和协方差来描述聚类的形状，所以这些类可以呈不同形状的椭圆形。采用最大期望（expectation-maximization，EM）优化算法找到数据集的均值和协方差。分为E步和M步两个主要的迭代步骤，其中混合高斯模型如式（12）所示。

p ( x ) = {\displaystyle \sum _{i=\mathrm{1}}^L} {\alpha }_i\cdot p(x|{\mu }_i,Co{v}_i) （12）

L L ( D ) = \mathrm{l} \mathrm{n} ( {\displaystyle \prod _{j=\mathrm{1}}^N} p(x_j)) = {\displaystyle \sum _{j=\mathrm{1}}^N} ln({\displaystyle \sum _{i=\mathrm{1}}^L}{\alpha }_i\cdot p(x_j|{\mu }_i,Co{v}_i)) （13）

p（x| μ_i， \mathrm{C} \mathrm{o} {\mathrm{v}}_{\mathrm{i}}）代表假设的各聚类高斯子模型，μ_i 、 \mathrm{C} \mathrm{o} {\mathrm{v}}_{\mathrm{i}}分别是对应的均值和协方差。式（13）是极大似然函数的对数似然函数，估计h_ji。

E ( h_{ji} | x_j ) = \frac{{\alpha }_{i}p(x_j|{\mu }_i,\mathrm{C}\mathrm{o}{\mathrm{v}}_{\mathrm{i}})}{{\displaystyle \sum _{l=\mathrm{1}}^L}{\alpha }_{l}p(x_j|{\mu }_l,\mathrm{C}\mathrm{o}{\mathrm{v}}_{\mathrm{i}})} （14）

式中：h_ji 为x_j 属于第i类的概率；E（h_ji | x_j ）是其期望。EM算法的具体步骤如下：

E步：将n-1次迭代中求得的各参数值带入到期望函数中，求得第n次迭代中的期望值。

M步：用极大似然估计对式（13）求各参数偏导，代入E步中求第n次迭代期望值，求解第n次迭代中的参数值，循环求解，最终得到样本个体属于每个集群的概率。

2.2 改进的ALNS

元启发式算法主要依靠大量种群个体不断按一定机制寻找可行域中的不同点，逐步逼近局部最优解或全局最优解。本文在优化车辆路径过程中使用ALNS算法进行求解，并改进算法中的删除因子及修复因子，加入得分系统，使算法能够在迭代的不同阶段自适应选择合适的变换因子，达到加快求解速度的目的。

2.3 邻域变换

本节先介绍在算法迭代过程中客户及子路径优劣的比较方法，再介绍删除、修复因子。

2.3.1 客户分摊成本方法

从式（1）中可得，因P在假设中为固定值，与客户i相关的路径使用成本只与 c_{ij} \mathrm{、} c_{j^{\text{'}}i} ( j \ne j^{\text{'}} )及F_i 有关。所以，可通过 \mathrm{W} \mathrm{e} \mathrm{i} \mathrm{g} \mathrm{h} {\mathrm{t}}_i = ( ( {c_i}_j + c_{j^{\text{'}}i} ) / \mathrm{2} ) + F_i计算每次迭代时每个粒子中客户分担的成本。

2.3.2 子路径效率比较方法

以客户分摊成本为基础，通过计算每条子路径客户分摊成本的平均值来比较不同子路径间的优劣程度，如式（15）所示。

\begin{array}{cc}\mathrm{A}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{W}\mathrm{e}{\mathrm{i}}_S=\sum _{i\in S}\mathrm{W}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}{\mathrm{t}}_i|S|,& S\subseteq V\end{array} （15）

2.3.3 客户相似性比较方法

客户共需以下3个属性：坐标、需求、服务时间窗。坐标要素计算两客户间距离并除以点集内两点间的最大距离；需求要素计算目标客户需求与客户最大需求的商；服务时间窗要素计算两客户时间窗最早服务时间和最晚服务时间绝对差值的和，并除以最大时间窗。此3个属性的单位不一致，若直接使用会使某些属性值过大或过小，无法体现客户间所有属性的内在关系。因此，在比较客户间关系时，需使用相对值进行比较，如式（16）所示。

\mathrm{A} \mathrm{t} \mathrm{t} \mathrm{r} \mathrm{i} \mathrm{b} \mathrm{u} \mathrm{t} {\mathrm{e}}_{ij} = \frac{c_{ij}}{\mathrm{M}\mathrm{a}\mathrm{x}(c_{ij})} + \frac{\left|d_i-d_j\right|}{\mathrm{M}\mathrm{a}\mathrm{x}(d_i)} + \frac{\left|l_i-l_j\right|+\left|u_i-u_j\right|}{\mathrm{M}\mathrm{a}\mathrm{x}(l_i-u_i)} （16）

2.3.4 删除因子

采用5种删除因子，分为单客户删除方法及多客户删除方法两类。单客户方法为随机单客户及非随机客户两类距离删除因子。随机单客户删除因子由swap^［16］及switch变换转变而来。客户距离删除因子计算每个客户到最近车场的距离与到当前子路径车场的距离的差值，排序后选择。某客户差值越大，则表示其越需要变换现有路径位置。

多客户删除方法有随机子路径删除因子、最差子路径删除因子及Shaw删除因子。随机子路径删除因子^［17］是随机从被选择的某个迭代粒子中选择一条子路径删除，再将删除的客户重新插入的方法。最差子路径删除因子是比较每条子路径的AverWei _S，从大到小排序，运用轮盘赌方式选择需要删除的子路径，并将此子路径中的客户序号输出到待删除列表中。Shaw删除因子由Shaw提出，是通过寻找与初始删除客户属性相似的其他多个客户，并重新插入的方法。

2.3.5 修复因子

修复因子使用约束内修复因子、无约束修复因子和最差起始修复因子。3种修复因子都以插入后增加的子路径目标函数值为标准，按照轮盘赌方式进行选择重排位置。约束内修复因子在重排客户时需要考虑被插入子路径时车辆能否承载此客户，若不能则选择其他子路径。无约束修复因子指客户点插入时不用考虑是否会超重，若有超重则在此子路径中使用删除因子选定待删除客户并按照约束修复因子重排。最差起始修复因子主要针对多客户删除后的重排顺序。与随机选择不同，其主要按照删除客户的分摊成本Weight _i 进行排序，再按从大到小的顺序用轮盘赌进行选择。

结合随机单客户删除因子与修复因子，生成两种变换因子：单客户无约束变换因子 {\Omega }_{\mathrm{1}}及单客户约束变换因子 {\Omega }_{\mathrm{2}}。客户距离删除因子 {\Omega }_{\mathrm{3}}采用最差起始修复因子选择需要重排的客户。多客户删除方法统一先与最差起始修复因子结合，对需要重新插入的客户集合排序，而后依靠无约束修复因子进行插入，得到随机子路径变换因子 {\Omega }_{\mathrm{4}}、最差子路径变换因子 {\Omega }_{\mathrm{5}}及Shaw变换因子 {\Omega }_{\mathrm{6}}。

2.4 评分规则及接受准则

2.4.1 评分规则

种群设置了当前总体最优解Gbest和每个粒子的当前个体最优解Pbest _i。假设当前解为S_i，变换后解为 S_i^{\text{'}}，且目标函数值求最小，则经过变换后的解存在以下3种情况。

（1） λ ₁： ，变换后的解比当前个体最优解更优；

（2） λ ₂： ，变换后的解优于当前解，但差于当前个体最优解；

（3） λ ₃： ，变换后的解差于当前解。

为评价每次变换的结果，对以上3种情况设置3个不同的分数累加统计。鉴于模型求目标函数值最小，令λ ₁得分为 {\omega }_{\mathrm{1}}、λ ₂得分为 {\omega }_{\mathrm{2}}、λ ₃得分为 {\omega }_{\mathrm{3}}，且 {\omega }_{\mathrm{1}} > {\omega }_{\mathrm{2}} > {\omega }_{\mathrm{3}}。

设所有变换方法初始权重都为1。每次迭代时，各粒子按方法权重大小选取需要的变换方法。依据比较规则判定此次变换方法的得分，并累加到对应方法的权重中。

2.4.2 接受准则

为保证迭代过程中粒子的不断优化，且使得粒子具有能够逃逸局部最优解的能力，本文使用SA算法的Metropolis接收准则判断变换后解 S_i^{\text{'}}是否能替换当前解S_i。设T ₀为初始温度，令当前温度T=T ₀，β为温度变化率。若变换后情况为λ ₁，则替换当前解和当前个体最优解；若为λ ₂，则用变换后解替换当前解；若为λ ₃则生成随机数α，并依照如式（17）所示的接受规则判断。

p_i = \mathrm{e} \mathrm{x} \mathrm{p} \left(\left(F\left(S_i\right)-F\left(S_i^{\text{'}}\right)\right)/T\right) （17）

若 ，则满足接受准则，替换当前解，并令 T = \beta T；否则，保持当前解不变。

2.5 禁忌链表及特赦规则

对每个粒子单独设置禁忌链表TB，定义其容量为size，定义remove及push两种操作。remove为删除TB中现有部分客户序号；push则将部分客户序号链接到TB末端。每次执行push操作时需判断TB何时达到容量上限，若达到上限则删除与末端剩余未添加的客户子集相同大小的顶端客户子集，再将剩余客户添加进TB中。

记录每次变换时重排的客户点集tb_i，若在某次迭代时需要删除的客户子集 t b_i \bigcap T B = \mathrm{\varnothing }，则执行ALNS算法的后续步骤，并令 T B = T B \bigcup t b_i；若 t b_i \bigcap T B \ne \mathrm{\varnothing }，则先执行插入因子，并令 \mathrm{c} \mathrm{o} \mathrm{m} \mathrm{m} \mathrm{o} {\mathrm{n}}_i = t b_i \bigcap T B，得到变换解 S_i^{\text{'}}后依靠评分规则进行判断，即若为λ ₁或λ ₂则执行特赦规则，先执行remove操作，即 T B = T B - \mathrm{c} \mathrm{o} \mathrm{m} \mathrm{m} \mathrm{o} {\mathrm{n}}_i，后执行push操作，即 T B = T B \bigcup S_i；若为λ ₃则重新选择需删除客户子集tb_i。使 t b_i \bigcap T B = \mathrm{\varnothing }，进行变换后按照接受准则，若变换解 S_i^{\text{'}}替换当前解S_i 则执行push操作，否则不执行操作。

3 仿真实验

本文共设计3组实验：实验一，使用汉明距离比较聚类初始方法与非聚类初始方法对算法初始解的影响，验证聚类初始化对算法的优化；实验二，构建灵敏性分析，比较不同的变换因子对算法的影响效果，验证算法鲁棒性；实验三，将本文算法与3种不同算法的结果进行比较，验证算法的有效性。

本文使用Java编译算法，在Windows10平台上运行，机器配置为i5-4200H CPU @ 2.80 GHz，16 G内存。算例使用MDVRPTW标准算例中的pr01~pr10进行计算。实验种群个数为150、速度为1 km/h，单位行驶成本P为每公里一元，两点间距离c_ij 单位为千米，由于本文讨论硬时间窗问题，将惩罚因子θ_i 设为每小时500元。为保证计算便捷，分别用相同的向量将以上算例的所有客户点坐标平移到坐标轴的第二象限，得到新的算例坐标。为使算法能够顺利执行，需要先确定算法的各项系数，即评分规则中的得分 {\omega }_{\mathrm{1}}、 {\omega }_{\mathrm{2}}、 {\omega }_{\mathrm{3}}。

3.1 参数设计

为使粒子在不同迭代阶段能自适应地选择合适的变换因子，需要运用合理的计分方式对每个变换因子进行评价。设初始状态下 {\omega }_{\mathrm{1}} = {\omega }_{\mathrm{2}} = {\omega }_{\mathrm{3}} = \mathrm{1}，种群个数为150，这样每次迭代时就有足够多的样本进行评价，且保证初始选择的概率一致。由于 {\omega }_{\mathrm{2}}及 {\omega }_{\mathrm{3}}都代表变换后解为 S_i^{\text{'}}优于当前解为S_i，可先设置 {\omega }_{\mathrm{2}} = {\omega }_{\mathrm{3}}，验证 {\omega }_{\mathrm{1}}与 {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}的关系。对 {\omega }_{\mathrm{1}}、 {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}分别设置一组数据，计算在不同系数组合下算法计算10次得出的结果。使用MDVRPTW模型的算例pr01进行验证。目标函数平均值变化情况见图1。

图1a、1b分别为不同系数组下的10次计算中目标函数的平均值变化及最小值变化，方框灰度越深表示值越小。由图1a可看到图右下部分颜色较左上部分更深，说明 {\omega }_{\mathrm{1}}小于 {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}时的值普遍要大于 {\omega }_{\mathrm{1}}大于 {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}时的值；当 {\omega }_{\mathrm{1}}值较小时，图1a中灰度深的区块更多，说明模型取得较优解的次数更多。由图1b可看到，当 {\omega }_{\mathrm{1}}在［0.5，1.5］， {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}在［1.5，3］时模型能够得到较优解。综合两图比较稳定性和优化能力后可确定 {\omega }_{\mathrm{1}} = \mathrm{0.5}较合适。在确定 {\omega }_{\mathrm{1}} = \mathrm{0.5}后还需确定 {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}比较合适的值。为合理分区，在区间［1.4， 3］设5组系数组，并得到如图2的数据。

如图2所示， {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}的值在变化时，模型目标函数值并无明显变化，当 {\omega }_{\mathrm{2}} \mathrm{、} {\omega }_{\mathrm{3}}的值较大时，图中深色色块个数更多，表明取较大值时得到更优解的概率更高。

结合图2可发现，当 {\omega }_{\mathrm{2}} = \mathrm{2.2}、 {\omega }_{\mathrm{3}} = \mathrm{2.6}时，模型10次计算中不论是平均结果还是最优结果都有较好表现。

3.2 实验一

为验证聚类算法对模型求解的优化，引入汉明距离对比聚类初始化和非聚类初始解集的多样性差异，对初始种群中所有不同个体进行比较，由实验种群个数可知样本个数为11 175。其结果如图3所示，图3b~3g分别为算例pr01到pr06的实验结果。

由图3中的6组实验可看出聚类初始化在6种不同算例中生成的初始种群，在多样性程度上较非聚类初始化有不同程度减弱，但从图3中可知这种削弱随算例客户数的增多而减弱，在算例pr04及pr06实验中聚类初始种群与非聚类初始种群的多样性均值及方差的差值比分别为0.035 8、0.112及0.0117、0.191。

在目标函数值的影响方面，由图3a可看出，随客户量增加聚类初始种群目标函数值优势逐渐增大，在算例pr05及pr06实验中平均值、最小值、最大值的差值比分别为0.114、0.105、0.102及0.094、0.105、0.11。

综上，在初始化种群效果方面，聚类初始化方法较非聚类初始化方法有较大改进。

3.3 实验二

本文共使用6种变换因子，为检测不同变换因子在算法中的作用，将此实验分为两个阶段。第一阶段将6种变换因子分为7种不同组合，混合表示运用所有变换因子，缺Ω ₁表示只使用其他5种变换因子，而不使用Ω ₁，其余同理，保证组合间各不相同。分别对算例pr01~pr04进行计算，并比较结果。

从表1的计算数据可以看出，每组变换因子对应4个算例的结果大致相同。其中混合6种变换因子的结果最优，而缺少 {\Omega }_{\mathrm{3}}及 {\Omega }_{\mathrm{4}}会对算法结果产生较大影响，算法结果平均增加0.891%和0.942%。

分析数据可知，缺少 {\Omega }_{\mathrm{1}}、 {\Omega }_{\mathrm{2}}、 {\Omega }_{\mathrm{3}}时结果增加0.024 7%。就变换因子分类而言，单客户变换因子对模型结果的影响比多客户变换因子的影响更大。在单客户变换因子种类中，缺 {\Omega }_{\mathrm{1}}比缺 {\Omega }_{\mathrm{2}}的影响更大，说明在路径变换中先不考虑车辆载荷约束的情况下，能使算法得到更好结果。但单从结果上看并不能完全解释每个变换因子对算法运行期间的作用，需比较每种组合在算法不同迭代时期的当前全局最优解的变化情况。

算例pr01在不同组别下前400次迭代时的目标函数值变化情况如图4所示。

图4中，pr02、pr03、pr04这3个算例迭代趋势与pr01相似。从图4中前200次迭代趋势可看出，在缺少 {\Omega }_{\mathrm{4}}、 {\Omega }_{\mathrm{5}}时算法收敛速度会减慢，尤其是缺少 {\Omega }_{\mathrm{4}}的情况更为明显。在缺少 {\Omega }_{\mathrm{6}}时前100次迭代的收敛速度比其他组别更快，但到50次迭代后其效率差于其他组别，且在170~220次迭代期间，其当前全局最优解已弱于混合组当前全局最优解，这说明 {\Omega }_{\mathrm{6}}对算法后期迭代收敛有一定影响。与其类似，缺少 {\Omega }_{\mathrm{3}}时算法迭代效率在50~100次迭代中开始减弱。

第二阶段分别单独测量每一种变换因子对算法的影响。记录不同迭代时期的当前全局最优解和当前得分。

表3　单变换因子得分比较

从表2和表3可以很明显看出单独使用 {\Omega }_{\mathrm{5}}时算法无法收敛。结合图4可知 {\Omega }_{\mathrm{5}}的主要作用是前期优化种群现有解集，就结果而言， {\Omega }_{\mathrm{2}}的结果最优，而得分情况则是 {\Omega }_{\mathrm{1}}最优，对比两类数据可知在单次优化能力上 {\Omega }_{\mathrm{2}}要强于 {\Omega }_{\mathrm{1}}。就变换类型而论，单独使用单客户变换类型要优于多客户变换类型，这也验证了表1的数据结果。

图5为单变换因子比较情况，其分别为6组变换因子计算10次的当前最优解变化趋势图及因子得分趋势比较图。此图以pr01数据为例，算例pr02、pr03、pr04数据结果与其类似。图5a中，由于 {\Omega }_{\mathrm{5}}结果太差，所以不再显示。

从图5a中可以看出迭代前期收敛速度最快的是 {\Omega }_{\mathrm{4}}， {\Omega }_{\mathrm{6}}收敛能力最差，在340~380次迭代期间 {\Omega }_{\mathrm{4}}的当前最优解值已弱于 {\Omega }_{\mathrm{2}}、 {\Omega }_{\mathrm{3}}，随后更弱于 {\Omega }_{\mathrm{1}}。由此也可以说明 {\Omega }_{\mathrm{4}}后期跳出局部最优解的能力比 {\Omega }_{\mathrm{1}}、 {\Omega }_{\mathrm{2}}、 {\Omega }_{\mathrm{3}}要差。 {\Omega }_{\mathrm{3}}的情况与 {\Omega }_{\mathrm{4}}类似，其平均值在507次迭代后便弱于 {\Omega }_{\mathrm{2}}。这一点从图5b也可以明显看出， {\Omega }_{\mathrm{3}}和 {\Omega }_{\mathrm{4}}的得分比在短暂上升后便呈现下跌趋势，直到后期趋于平稳。虽然结果并不理想，但从图5b中得分比变化趋势可以看出， {\Omega }_{\mathrm{6}}在后期依然能够对解产生比其他变换因子更好的优化作用，证明了图4中描述的 {\Omega }_{\mathrm{6}}作用的观点。两阶段实验证明了本文算法使用变换因子的有效性及鲁棒性。

3.4 实验三

将本文算法与TS算法、变邻域搜索算法（variable neighbourhood search algorithm，VNS）及混合遗传算法（hybrid genetic algorithm，HGA）^［18］进行比较，求解最短路径的长度，目标函数如式（18）所示。

\mathrm{m} \mathrm{i} \mathrm{n} : F = {\displaystyle \sum _{k=\mathrm{1}}^K} \sum _{i\in V}\sum _{j\in V}{c}_{ij}{x}_{ij}^k （18）

分别计算pr01~pr06共6组算例10次，在表4中列出计算结果。%符号表示本文算法与表中最优解的差值百分比。表4中所有时间的单位为min。从表4中可看出，本文算法在pr01中的结果超过其他3种算法的结果，本文算法在pr06和pr10两个拥有288个客户的算例上有一定差距。在多数算例中本算法在时效上有一定优势，如pr01求解时间为0.44 min，较HGA算法节省0.22 min。pr05求解结果所花时间较HGA算法节省9.83 min，较VNS节省84.78 min。

4 结论

本文结合客户的空间属性，使用变换因子得分系统，提出结合GMM聚类的改进自适应大邻域搜索算法。该算法能够在迭代的不同阶段自适应地选择合适的变换因子以此提高算法收敛速度。本文设计了3组实验，实验一通过比较聚类初始化与非聚类初始化验证了初始化方法对算法初始解的优化；实验二针对不同变换因子设计了灵敏度实验，验证了不同变换因子的作用，说明算法的有效性及鲁棒性；实验三将本文算法与其他3种算法进行了比较，验证了本文算法高效的计算效率。

今后可对本文算法继续改进，提高算法求解能力。本文仅假设车辆速度恒定不变，但实际生活中车辆速度不一定恒定，下一步可研究车辆变速情况下的MDVRPTW问题。

参考文献

[1] Afshar N B， Afshar-Nadjafi A.A constructive heuristic for time-dependent multi-depot vehicle rou-ting problem with time-windows and heterogeneous fleet［J］.Journal of King Saud University-Engineering Sciences，2017，29（1）：29-34.

[2] Kubil V N， Mokhov V A， Grinchenkov D V.Multi-objective ant colony optimization for multi-depot heterogenous vehicle routing problem［C］//2018 International Conference on Industrial Engineering，Applications and Manufacturing.Moscow：IEEE，2018：1-6.

[3] Reza， Tavakkoli M， Meskini M，et al.A multi-depot close and open vehicle routing problem with heterogeneous vehicles［C］//Industrial Engineering and Systems Management.Shanghai：IEEE， 2019， 154-159.

[4] 凌海峰，谷俊辉.带软时间窗的多车场开放式车辆调度［J］.计算机工程与应用，2017，53（14）：232-239.

[5] Fernando B D， Rasul E， Hossein J S，et al. A cooperative coevolutionary algorithm for the multi-depot vehicle routing problem［J］. Expert Systems with Applications， 2016， 43（1）： 117-130.

[6] 李宁，邹彤，孙德宝.带时间窗车辆路径问题的粒子群算法［J］.系统工程理论与实践，2004，24（4）：130-135.

[7] 宾松，符卓.求解带软时间窗的车辆路径问题的改进遗传算法［J］.系统工程，2003，21（6）：12-15.

[8] Goeke D， Schneider M.Routing a mixed fleet of electric and conventional vehicles［J］.European Journal of Operational Research，2015，245（1）：81-99.

[9] Fan H M， Zhang Y G， Tian P J，et al.Time-dependent multi-depot green vehicle routing problem with time windows considering temporal-spatial distance［J］.Computers & Operations Research，2021，129（6）：268-281.

[10] 吴凡，杨冰，洪思.基于变长基因型遗传算法的多供应点应急物资调度优化［J］.计算机应用研究，2022，39（4）：1148-1154.

[11] Hesam Sadati M E， Çatay B， Aksen D.An efficient variable neighborhood search with tabu shaking for a class of multi-depot vehicle routing problems［J］.Computers & Operations Research，2021，133（4）：1-22.

[12] Wu Y， Yang W， He G C，et al.An improved adaptive large neighborhood search algorithm for the heterogeneous fixed fleet vehicle routing problem［C］//2017 8th IEEE International Conference on Software Engineering and Service Science.Beijing：IEEE，2017：657-663.

[13] Wang Y， Li Q， Guan X Y，et al.Collaborative multi-depot pickup and delivery vehicle routing problem with split loads and time windows［J］.Knowledge-Based Systems，2021，231（6）：1-24.

[14] 鲁建厦，洪欢蕾，陈青丰.考虑供给商品价格的多车场车辆路径问题［J］.浙江工业大学学报，2016，44（5）：553-558.

[15] 王征，张俊，王旭坪.多车场带时间窗车辆路径问题的变邻域搜索算法［J］.中国管理科学，2011，19（2）：99-109.

[16] Deng Y F， Xiang J L， Ou Z L.Improvement of genetic algorithm for vehicle routing problems with time windows［C］//2013 Third International Conference on Intelligent System Design and Engineering Applications.Hong Kong：IEEE，2013：866-869.

[17] Mühlbauer F， Fontaine P.A parallelised large neighbourhood search heuristic for the asymmetric two-echelon vehicle routing problem with swap containers for cargo-bicycles［J］.European Journal of Operational Research，2021，289（2）：742-757.

[18] 洪联系，董绍华.MDVRPTW问题多阶段迭代启发式算法［J］.计算机工程与应用，2007，43（26）：217-222.

基金资助

国家自然科学基金(61972266)

辽宁省自然科学基金(2020-MS-233)

辽宁省兴辽英才计划(XLYC2002017)