基于改进变邻域搜索的多隔室车辆路径优化算法

姚冠新，范雪茹，张冬梅

(1.江苏大学 管理学院，江苏 镇江 212013；2.扬州大学 江苏现代物流研究基地，江苏 扬州 225009)

1 问题的提出

随着人们生活品质的提高，多品种、小批量、定制化的产品需求逐渐增加，越来越多的物流企业和配送中心采用多隔室车辆(Multi-Compartment Vehicles, MCVs)代替单隔室车辆(Single-Compartment Vehicles, SCVs)来实现不同种类(或特性)产品的联合交付[1]，既降低了配送总成本[1-2]，又满足了消费者需求。SCVs的整个车厢是一个独立的空间[3]，仅能运输一种特性的产品，如冷冻产品；MCVs的车厢被分割成多个相对独立的空间[3]，能够装载并共同运输不同种类的、不相容的、具有严格分装要求的产品[4-5]，如冷藏和冷冻产品。多隔室车辆路径优化问题(Multi-Compartment Vehicle Routing Problem, MCVRP)是具有容量限制的车辆路径优化问题(Capacitated Vehicle Routing Problem, CVRP)的扩展。CVRP即SCVs从配送中心出发，为需求节点配送一种特性产品，当所装载产品配送完毕后返回配送中心，其中每个节点仅能被访问一次，最终实现配送距离最短或配送成本最低[6]；MCVRP则是采用MCVs，设计运输线路以满足一组节点对多种特性产品的需求[4]，产品需装载在同一车辆的不同隔室内进行共同运输[2,7]，从配送中心出发依次抵达需求节点进行配送，配送完毕后返回配送中心，最小化配送距离或配送成本。以1个配送中心和3个需求节点片段为例，在总负载上限和节点总需求相同、车辆类型和产品细分需求不同、任意整车和隔室不超载且分装要求严格的前提下，一辆SCV即可满足CVRP中节点ai-1、ai、ai+1的配送需求(如图1a)，但却需要两辆MCVs才能满足MCVRP中节点bj-1、bj、bj+1的配送需求(如图1b)。由此可见，CVRP和MCVRP存在差异性，且由于条件严苛，MCVs的路径规划变得更加困难。

实际生产和生活中MCVs的应用逐渐增加。零售商的杂货分销中[1,5]，不同产品需储存在不同温度环境下运送[1]；废物和垃圾收集中[8-9]，不同颜色的玻璃、不同种类的垃圾需被严格分装在不同的隔室内进行联合运输；橄榄油收集[10]、成品油配送[4]、奶制品收集[11]、燃料分配或补给[11-12]、冷链配送[11,13]、动物饲料及牲畜集配[14]均会用到MCVs，其应用广泛且前景良好。此外，基于MCVs应用背景的MCVRP的理论研究热度也逐渐提升。近年来，关于MCVRP的相关研究主要集中在以下3方面：①基于不同情境的MCVRP模型创新和扩展，丰富了相关研究理论体系，如基于MCVRP的车辆选择[5]、多车库的MCVRP[11]、带时间窗的MCVRP[13]等；②求解MCVRP的元启发式算法的创新设计及改进，降低了配送距离(成本)，提供了有效的路径规划手段，如混合蚁群算法[2]、改进粒子群算法[3]等；③具有实际约束的MCVRP的案例研究，为生产、生活实践提供了科学的理论依据，如不同型号的成品油配送、分类垃圾收集等(如表1)。

表1 多隔室车辆路径优化相关问题及算法设计

续表1

续表1

综合来看，MCVRP优化研究以及相应的算法创新设计兼具理论价值和实际意义。然而，MCVRP中多种特性产品不能混装的特质以及相应隔室的增加、隔室容量的限制造成了研究的差异性及困难性，且MCVRP的平均解决方案质量不如CVRP[6]。因此，为优化MCVRP解决方案，本文基于变邻域搜索算法(Variable Neighborhood Search, VNS)框架，针对MCVRP的特点设计了一种改进变邻域搜索算法(Improved Variable Neighborhood Search, IVNS)。首先，设计多起点寻优机制，采用扫描法构造历次迭代的初始解，增加解的多样性；其次，细分主程序下子程序1和子程序2，子程序1中设计Shaking过程，以实现小范围解空间的探索，子程序2中设计全局扰动过程，以实现大范围解空间的探索，并设计一种还原及再分配策略处理不可行解，探寻解空间中的不可行区域；再次，结合使用贪婪算法和多种混合算子来设计Local Search过程，以实现局部搜索和优化，提升求解质量和计算效率；然后，应用最大迭代次数停止准则结束主循环并保留最优解；最后，采用改编算例进行实验及结果对比分析，验证了IVNS算法的有效性、优越性及收敛性。

2 问题描述及模型构建

2.1 问题描述

多隔室车辆路径优化问题(MCVRP)可以描述为一个配送中心服务于多个需求节点，每个需求节点具有不同种类产品需求，任意一类产品需求不可拆分配送，不同种类产品需严格分装、联合配送交付，且仅接受一次车辆访问；需采用车型相同的、每个隔室的容量固定不变、具有容量限制的多隔室车辆进行配送，其中每个隔室仅能装载一种产品，每辆车负责线路中所有需求节点对某种产品的需求已知且不超过对应隔室的负载上限、对某种产品的需求总量不超过对应隔室的负载上限、对所有种类产品的需求总量不超过整车的负载上限；配送中心的库存最大容量、多隔室车辆的台数均能够满足所覆盖需求节点的所有种类产品需求；多隔室车辆从配送中心出发，沿规划路径配送完毕后，返回配送中心，最终实现总的配送距离最小化。

(1)相关符号G=(N,E)，N={0}∪N′={0,1,2,…,n}。其中：G表示配送系统;N表示物流节点集合;{0}表示配送中心;N′表示需求节点集合;E={(i,j)|i,j∈N}，E表示无向边集合，i,j表示物流节点，i,j∈N。dij表示物流节点i到j的距离，其中dij=dji，即两节点往返距离相同。V={1,2,3,…,k}，V表示多隔室车辆集合，k表示多隔室车辆，k∈V。M={1,2,3,…,m}，M表示车辆隔室集合。P={1,2,3,…,p}，P表示产品种类集合，其中m=p，即产品种类数与隔室数量相同。Qp表示某隔室所能够负载对应种类产品p的最大载重量，Qmax表示多隔室车辆的总的最大载重量。qjp表示需求节点j对产品p的需求量。

2.2 模型构建

MCVRP的数学模型如下：

(1)

s.t.

(2)

(3)

(4)

∀k∈V,∀p∈P；

(5)

∀k∈V,∀p∈P,i≠j；

(6)

∀j∈N,∀k∈V,∀p∈P,i≠j。

(7)

其中：式(1)为目标函数，表示最小化多隔室车辆总配送距离；式(2)表示一个需求节点只接受一辆多隔室车的服务且仅能被访问一次；式(3)表示路径具有连续性；式(4)表示多隔室车辆配送过程中的任意种类产品需求之和不超过对应的隔室负载上限；式(5)表示任意需求节点的任意种类产品需求量不超过相应隔室的负载上限；式(6)表示多隔室车辆返回配送中心时，某种产品的荷载量等于该多隔室车辆本次配送路线中各需求节点的该种产品的需求量之和；式(7)表示配送车辆在任意节点的某种产品的荷载量满足在该点需求量的递推关系。

3 算法设计

VNS算法是MLADENOVIC等[29]和HANSEN等[30]提出的一种优化算法，被广泛用于解决组合优化问题，其基本思想是在搜索过程中系统地改变邻域结构集以拓展搜索范围，通过局部搜索流程以获得局部最优解，再基于此局部最优解重复上述过程获得另一个局部最优解，如此不断迭代最终获得收敛，达到优化目的[21,31-32]。VNS算法的基本环节包括初始解构造、Shaking过程、Local Search过程及停止准则，具有环节少、参数少、易于操作和改进等优势。为了在合适的求解时间内获得MCVRP的高质量解决方案，本文基于VNS算法，针对MCVRP的特点，提出改进变邻域搜索算法(IVNS)，具体流程如图2所示。下面将详细介绍IVNS算法的初始解构造、Shaking过程、全局扰动过程、还原及再分配策略、Local Search过程及停止准则等重要组成部分。

3.1 初始解构造

VNS算法是基于单个解向量展开的寻优过程[31]，即从单点出发进行寻优，遍历所有邻域结构后保留最优解，求解的时效性较好，但是解的多样性较差。鉴于此，设计一种多起点寻优机制[8]，采用经典的扫描法构造每次迭代过程的初始解，基于MCVRP不同种类产品需求以及不同隔室的负载上限进行容量及路径可行性判定，具体步骤如下：

步骤1坐标系转换。将配送中心{0}和所有需求节点{1,2,3,…,n}所在的笛卡尔坐标系转换成以配送中心为极点、以任选1个需求节点和配送中心的连线为极轴的极坐标系，其他需求节点基于极轴转换为极坐标角度。

步骤2将需求节点分组。最小角度的需求节点作为起始点，按照逆时针方向，将需求节点逐一纳入组Group1中，直到该组中所有需求节点的某特性产品的需求总量超过该隔室负载上限，然后以本组最后1个需求节点作为下一组起始点建立新组Group2，继续扫描剩余需求节点并逐一添加到新组内。

步骤3重复步骤2，直至所有需求节点都被分组为止。此时，任意一组内所有需求节点的任意种类产品的需求总量未超过当前隔室的负载上限，且组内所有需求节点的所有种类产品需求总量未超过多隔室车辆整车负载上限；

步骤4形成初始解。分别将各组内的第1个需求节点、最后1个需求节点与配送中心连接，组内的第i个点与第i+1个点连接(i=1,…,j-1，其中j为当前组内需求节点总数)，每组各形成一条子路径，所有子路径共同构成初始解。

以1个配送中心和10个需求节点为例，初始解构造流程如图3所示。需注意，IVNS算法中可能出现较少点构成的子路径或较少条子路径构成的解，甚至出现极端情况，如1个点构成的子路径或1条子路径构成的解。若某段程序中需具备2个及以上需求节点构成的子路径或2条及以上子路径构成的解才能继续执行时，通过设计并执行特定程序来重新选择子路径或跳出循环，确保算法运行畅通。

3.2 Shaking过程

Shaking过程即系统地改变邻域结构集以拓展当前解的搜索范围，避免陷入局部最优解，从而获得更优解。鉴于此，基于MCVRP的子路径形式初始解设计子程序1中的Shaking过程，从子路径转换的角度构造邻域结构集合，小范围扩展当前解搜索空间，具体步骤如下：

步骤1定义邻域结构集合Nk(k=1,…,kmax)，其中kmax=9，即设计了9个邻域结构。

步骤2选取邻域结构。按照预先设定顺序，选取邻域结构集合中的一个邻域结构Nk。

步骤3选取两条子路径。确定当前解s包含的子路径的数量m，先选取当前解s的一条随机子路径Subpathi，其中i=1,2,…,m，若i不等于m，则再选取子路径Subpathi+1；若i等于m，则再选取子路径Subpath1。

步骤4获得两条(或一条)新子路径。根据Nk的定义，采用当前定义下的特定算子对所选择的两条子路径进行节点以及节点的不同种类产品需求转换操作，获得两条新子路径，若因节点变化导致仅剩一条子路径，则删除空子路径以不影响后续算法执行。

步骤5形成新解。通过保留当前解s的大部分特征，改变小部分特征，形成新解s′，加快算法收敛速度[32]。

单一的邻域结构可能造成当前解的较大波动[21]，而多个邻域结构可广泛探索解空间[8]，进而提高算法的求解效率和稳定性[32]。因此，基于插入、交换、2-opt*等经典算子来设计邻域结构集合，并按照固定顺序逐一实现邻域拓展动作。

首先，基于插入算子设计邻域结构N1～N4：1-插入、2-插入、3-插入、4-插入，即分别基于1～4个(连续)节点的插入动作。插入算子(如图4)的基本原理即指定两条子路径SubpathK(k1,…,ki-1,ki,ki+1,…,km)和SubpathT(t1,…,tj-1,tj,tj+1,…,tn)，提取SubpathK中ki插入到SubpathT中tj的后置位，形成两条新的子路径newSubpathK(k1,…,ki-1,ki+1,…,km)和newSubpathT(t1,…,tj-1,tj,ki,tj+1,…,tn)，然后判断是否demandQ(newSubpathK)≤vehiclecompartmentcapacityQ且demandQ(newSubpathT)≤vehiclecompartmentcapacityQ(Q=1,2,…,q，q为MCVs隔室的数量)，即判断每个隔室的负载是否合规，若否，则重新进行选择；若是，则进一步判断是否distance(newSubpathK)+distance(newSubpathT)

其次，基于交换算子设计邻域结构N5～N8：1-交换、2-交换、3-交换、4-交换，即分别基于1～4个(连续)节点的交换动作。交换算子的基本原理即指定两条子路径，提取并交换SubpathK中的ki和SubpathT中的tj，形成两条新的子路径newSubpathK(k1,…,ki-1,tj,ki+1,…,km)和newSubpathT(t1,…,tj-1,ki,tj+1,…,tn)，之后的判断、选择动作也需要结合还原及再分配策略并进行每个隔室的负载判断和每条子路径的可行性判断。此外，引入Cross-Exchange和iCross-Exchange两种算子，假设交换SubpathK中的(ki-1,ki,ki+1)与SubpathT中的(tj-1,tj,tj+1)，Cross-Exchange算子(如图5a)即将两个连续点片段分别在对方子路径中进行正向放置，形成newSubpathK(k1,…,ki-2,tj-1,tj,tj+1,ki+2,…,km)和newSubpathT(t1,…,ti-2,ki-1,ki,ki+1,ti+2,…,tn)；iCross-Exchange算子(如图5b)即将两个连续点片段在对方子路径中进行反向放置，形成newSubpathK(k1,…,ki-2,tj+1,tj,tj-1,ki+2,…,km)和newSubpathT(t1,…,ti-2,ki+1,ki,ki-1,ti+2,…,tn)。因此，出于对连续点片段的方向性以及解的可行性考虑[21]，设置基准概率P=0.8，随机生成0-1的随机数p，当pP时，应用iCross-Exchange算子。

最后，基于2-opt*算子设计邻域结构N9：2-opt*。2-opt*算子的基本原理即指定两条子路径，分别截断任意连续两点间线路(ki-1,ki)和(tj,tj+1)，获得SubpathK1(k1,…,ki-1)、SubpathK2(ki,ki+1,…,km)、SubpathT1(t1,…,tj-1,tj)以及SubpathT2(tj+1,…,tn)4个片段，将两条子路径的片段重新连接形成两种组合情形：情形1是newSubpathK(k1,…,ki-1,tj+1,…,tn)(K1和T2)和newSubpathT(t1,…,tj-1,tj,ki,ki+1,…,km)(K2和T1)(如图6a)；情形2是SubpathK(k1,…,ki-1,tj,tj-1,…,t1)(K1和T1)和newSubpathT(km,…,ki+1,ki,tj+1,…,tn)(K2和T2)(如图6b)，之后的判断、选择动作也需要结合还原及再分配策略并进行每个隔室的负载判断和每条子路径的可行性判断。此外，可设计一定概率获得不同组合情形以增加解的多样性，因此，设置基准概率T=0.5，随机生成0～1的随机数t，当tT时，采用情形2。

3.3 全局扰动过程

上述Shaking过程是基于子路径展开的小范围解空间的探索，为了进一步扩大解空间的搜索范围，通过有效的邻域变换保证多样化的探索，并最终实现全局的优化[20]，借鉴扰动方法设计思想[8]，基于总路径形式解设计子程序2中的全局扰动过程，从总路径转换的角度构造邻域结构集合，大范围扩展当前解搜索空间。具体地，任选总路径形式当前解的两点执行置换、转置、插入动作：置换即将原总路径Rout(r1,…,ri-1,ri,ri+1,…,rj-1,rj,rj+1,…,rn)中的ri和rj交换，获得新的总路径newRout(r1,…,ri-1,rj,ri+1,…,rj-1,ri,rj+1,…,rn)；转置即将ri和rj节点及之间节点顺序颠倒，获得新的总路径newRout(r1,…,ri-1,rj,rj-1,…,ri+1,ri,rj+1,…,rn)；插入即将ri插入到rj的后置位，获得新的总路径newRout(r1,…,ri-1,ri+1,…,rj-1,rj,ri,rj+1,…,rn)。需要注意的是，全局扰动过程也要结合还原及再分配策略进行之后的判断和选择，而且，MCVs每个隔室的负载判断及每条子路径的可行性判断均不可或缺，具体步骤如下：

步骤1输入当前解w的总路径形式解x，设置最大循环次数Maxl=9，i=0。

步骤2令i=i+1，结合判断、选择动作重复以下步骤，直至i>Maxl。

步骤3从x中任选两点xi和xj，产生1～3的伪随机整数m，进行以下操作以形成新排列顺序，记为x′。

(1)当m=1时，对xi和xj两节点执行置换操作(如图7a)；

(2)当m=2时，对xi和xj两节点及之间执行转置操作(如图7b)；

(3)当m=3时，对xi和xj两节点执行插入操作(如图7c)。

步骤4将x′的子路径形式解记为w′，实现了解空间较大范围拓展搜索。

3.4 还原及再分配策略

一般算法设计中，相应动作的执行以新解负载合规、配送距离变小为前提[6]，导致部分解被过于严苛的约束条件剔除，此时可采用一些特殊的策略和机制保留并处理不可行解以探索解空间中的不可行区域[6,8,21,23]，从而对解空间进行更大范围的探索，增加解的多样性的同时，提升获得更优解的可能性[26]。IVNS算法中，经过Shaking过程和全局扰动过程可能获得不可行新解，即子路径中某种产品总需求超过当前隔室的负载上限(隔室超载)或所有产品的总需求超过MCVs整车的负载上限(整车超载)，实际运输中不可行。因此，设计一种还原策略和再分配策略，保留并处理不可行解，将其还原成总路径形式解并再度分配成子路径形式的可行新解。

还原策略的具体步骤如下：

步骤1输入子路径形式的不可行当前解r，计算其子路径的数量o，即为分组数量o。

步骤2对第1组到第o组执行以下操作：去掉节点与节点、节点与配送中心之间的所有路径，并保存所有分组内的未去掉路径之前的节点顺序，即获得总路径形式的解y(如图8a)。

再分配策略的具体步骤如下：

步骤1基于总路径形式的当前解y，从第1组的第1个需求节点开始按照原顺序进行重新分组，假设分成j组，任意组内的负载判断均合规，隔室和整车均未超载。

步骤2第1组到第j组重新形成子路径，形成子路径形式的可行新解r′(如图8b)。

此外，由于Shaking过程和全局扰动过程基于的MCVRP解形式不同，还原及再分配策略在子程序1和2中的顺序也不同：子程序1中，基于子路径形式解执行Shaking过程，随后采用还原策略和再分配策略，最后进行局部搜索和优化(如图9a)；子程序2中，先执行还原策略获得总路径形式解，为全局扰动过程做准备，随后执行再分配策略，最后进行局部搜索和优化(如图9b)。可见，还原及再分配策略的使用使得不止两条子路径上的节点数量、顺序、总需求等发生变化。

3.5 Local Search过程

Local Search过程即通过局部搜索流程求得局部最优解，决定着最终的求解质量和计算效率[32]。为了对子程序1和子程序2产生的邻域解进行局部搜索以获得优化，借鉴文献[6,17-18]的思想，设计基于子路径内和子路径间搜索优化机制的Local Search过程。

基于经典2-opt算子和贪婪算法设计子路径内搜索优化机制。2-opt算子(如图10)即针对当前子路径SubpathK(k1,…,ki-1,ki,ki+1,…,kj-1,kj,kj+1,…,km)，选择两个需求节点ki和kj，将两需求节点及之间节点顺序颠倒，获得新子路径newSubpathK(k1,…,ki-1,kj,kj-1,…,ki+1,ki,kj+1,…,km)；贪婪算法即针对当前子路径，每次只选择距离当前节点最近的节点作为下一个配送节点，直至遍历当前子路径中所有节点。需注意，两者均是基于子路径内部节点的动作，每个隔室及整车负载均未变化，仅需要判断是否distance(newSubpathK)

基于2-opt*、插入、交换等经典算子(基本原理见3.2节)设计子路径间搜索优化机制。需注意，算例分析及实际应用中，互为近邻的子路径更容易产生无效交叉部分，消除部分交叉能够减少配送距离，而任意选择两条子路径进行优化时，获得较优解的概率较小，且会大幅增加计算时间。因此，扩展文献[6,22]的思想，以加快搜索速度、减少计算量为目的，设计一种子路径近邻优先选取原则，即先后选取相邻和相隔子路径进行局部搜索，相邻子路径即Subpathi和Subpathi+1(如图11a)，相隔子路径即Subpathi和Subpathi+2(如图11b)。此外，还采用best改进策略来实现求解质量和运行时间的最佳平衡[33]，即遍历并计算当前解的所有邻居解并保留最优解。特别地，子路径间搜索优化机制伴随着节点及节点需求转换，因此，MCVs的每个隔室都要重新进行负载判断。

据此，设计Local Search过程包含两个组成部分：①基于2-opt算子、贪婪算法交替使用的子路径内搜索优化机制;②基于近邻优先选取原则的、多种混合算子的子路径间搜索优化机制。其中，子路径间搜索优化机制包含相邻子路径2-opt*、相隔子路径2-opt*、相邻子路径1-插入、相隔子路径1-插入、相邻子路径1-交换、相隔子路径1-交换6个搜索流程。需注意，IVNS算法部分过程基于相同的经典算子原理，但是都进行了一定程度的结合创新和改进设计，因而具有差异性。

3.6 停止准则

与IVNS算法主程序的多起点寻优机制相对应，采用最大迭代次数停止准则，迭代一定次数后结束循环，在历次迭代的较优解中保留最优解。其中，子程序1停止准则即遍历Shaking过程的9种邻域结构后，结束循环并保留最优解；子程序2停止准则即对全局扰动过程循环9次后，结束循环并保留最优解。

4 算例分析

MCVRP缺乏国际标准算例，本文基于REED等[9]的改编办法对14个CVRP国际标准算例(http://neo.lcc.uma.es/vrp/vrp-instances/capacitated-vrp-instances/)进行改编，获得28个MCVRP算例。具体改编办法：①坐标：配送中心和需求节点的坐标不变；②车辆：原单隔室车辆容量按照3∶1的比例拆分成两个隔室容量；③节点需求：当需求节点的位置,x,y∈(0,35)时，其需求量按照3∶1的比例分割成两种特性产品需求量，其他需求节点的需求量按照2∶1的比例(S2改编方式)或4∶1的比例(S3改编方式)分割成两种特性产品的需求量。本文采用MATLAB R2016a进行计算，运行环境是Windows10专业版、Intel(R)Core(TM)i5-4210H CPU @ 2.90 GHz处理器、4 G内存、64位操作系统。通过改编算例的实验以及结果对比，分析并验证算法的有效性、优越性及收敛性。

4.1 算法的有效性分析

基于vrpnc2-S2算例，在其他程序不变的情况下，仅测试当前动作对结果的影响，每个变换下均运行5次、迭代500次，获得最好解、最差解、平均解及标准方差，所有表格中加粗数据为对比之下的最优解；此外，由于最好解、最差解、平均解等指标能够反映算法的寻优能力，极差(最差解—最好解)、标准方差、极差百分比(极差/最好解)等指标能够反映算法的求解稳定性[34]，故再计算极差和极差百分比来绘制寻优能力和求解稳定性统计分析图并进行统计学分析。

4.1.1 初始解策略的有效性

分别采用贪婪算法、随机策略、扫描算法产生初始解。结果显示，基于扫描法的IVNS的最好解、最差解、平均解均最小，随机策略次之，贪婪算法最大。相对于贪婪算法和随机策略，基于扫描算法的IVNS的平均解分别减少了10.902 8和7.482 2，减少百分比分别是1.23%和0.85%。基于扫描算法的IVNS的结果标准方差最小，求解稳定性最好(如表2)。此外，观察不同构造策略的初始解和最好解的折线图(其中：G为贪婪算法，R为随机策略，S为扫描算法；I为初始解，O为最好解)，5次运行中随机策略的初始解的波动最大，极差达到84.351 4，扫描法的初始解的波动最小，极差仅有22.824 7，说明扫描法能够更稳定地获得更优初始解(如图12)。可见，以扫描法构造初始解的IVNS算法更有效。

表2 不同初始解构造策略对结果的影响

4.1.2 Shaking过程算子的有效性

Shaking过程中依次叠加使用插入、交换、2-opt*算子。结果显示，相比于未采用任何算子，基于插入、交换、2-opt*算子的平均解减少了8.895 4，减少百分比为1.01%，结果质量更好，标准方差减少了0.90，求解稳定性更强。其中插入算子对解质量的提升最大，平均解减少值为6.566 5，为总减少值的73.82%(如表3)。此外，观察Shaking过程不同算子设计的寻优能力和稳定性相关指标的变化趋势，最好解、最差解及平均解均逐渐降低，而极差、标准方差和极差百分比变化有波动，但整体呈下降趋势，说明所设计Shaking过程算子的寻优能力更强、稳定性更高(如图13a和图13b)。可见，Shaking过程使用插入、交换、2-opt*算子使IVNS算法更有效。

表3 Shaking过程算子对结果的影响

4.1.3 全局扰动过程的有效性

在Shaking过程后增加全局扰动过程。结果显示，相比于仅有Shaking过程，结合了全局扰动过程的IVNS算法减小了最好解、最差解、平均解，略增加了标准方差。其中，平均解减少了0.525 4，提升了求解质量，跳出了局部最优，标准方差增加了0.36，降低了稳定性(如表4)。此外，观察全局扰动过程的寻优能力和求解稳定性的相关指标变化趋势，最好解、最差解、平均解均呈下降趋势，但是极差、标准方差、极差百分比均呈上升趋势，说明全局扰动过程的设计使得IVNS算法的寻优能力变强，虽然导致求解稳定性变差，但与全局扰动过程的大范围解空间探索的设计初衷吻合(如图13c和图13d)。可见，全局扰动过程的设计和应用使得IVNS算法更有效。

表4 全局扰动过程对结果的影响

4.1.4 Local Search过程优化机制及算子的有效性

Local Search过程中依次叠加使用子路径内和子路径间搜索优化机制。结果显示，子路径内及子路径间搜索优化机制均有效降低了最好解、最差解、平均解及标准方差，两者导致平均解分别减少了250.775 6和47.143 5，减少百分比分别是21.37%和4.02%，标准方差减少了19.41和1.27，提升了求解质量及稳定性(如表5)。此外，观察不同搜索优化机制的寻优能力及求解稳定性相关指标变化趋势，最好解、最差解、平均解以及极差、标准方差、极差百分比均呈较大幅度下降趋势，说明所设计搜索优化机制的寻优能力和求解稳定性均实现较大幅度提升(如图14a和图14b)。可见，Local Search过程中子路径内及子路径间搜索优化机制的设计提升了IVNS算法的有效性。

表5 Local Search过程子路径内及子路径间搜索优化机制对结果的影响

进一步地，子路径间搜索优化机制中依次叠加使用2-opt*、插入、交换算子。结果显示，2-opt*、插入、交换算子导致平均解分别减少了33.990 8、7.391 8及5.760 9，优化效果2-opt*最大，插入算子次之，交换算子最小(如表6)。此外，观察子路径间搜索优化机制中不同算子设计的寻优能力和求解稳定性相关指标变化趋势，最好解、最差解、平均解均呈下降趋势，且最好解和最差解的差距始终较小，而极差、极差百分比呈“先上升、后下降”趋势，但整体上极差、极差百分比和标准方差是下降的，说明所设计的子路径间搜索优化机制中的算子寻优能力和求解稳定性均加强(如图14c和图14d)。可见，Local Search过程中子路径间搜索优化机制2-opt*、插入、交换算子的采用使得IVNS算法更有效。

表6 Local Search过程子路径间搜索优化机制中算子对结果的影响

4.2 与已有文献结果对比

4.2.1 配送路径规划合理性比较

梳理现有MCVRP研究文献，REED等[9]和陈久梅等[3]的研究中绘制了vrpnc1-S2和vrpnc1-S3算例的结果网络结构图(配送路径规划图)，其中REED、陈久梅以及本文的最好结果如表7所示，相应算例目前最好的，以及本文最好结果的网络结构如图15～图18所示。相比于IPSO结果，IVNS算法所获得的vrpnc1-S2、vrpnc1-S3算例的最好解分别减少了1.494 2和1.729 2，减少百分比分别为0.27%和0.31%，结果质量有所提升，且IVNS算法所获得的网络结构图减少了无效交叉路径，降低了配送距离，均优于ACO和IPSO求解结果。

表7 与已有文献的总配送距离的比较

4.2.2 求解稳定性比较

梳理相关文献，陈久梅等[3]的研究中展示了多次运行结果的最好解、最差解、平均值、标准方差等具体数据，结果新鲜、全面且优质，故将本文结果与之进行求解稳定性的比较。陈久梅等[3]采用IPSO，运行10次，迭代次数1 000次。对比结果显示，在有限但不相同的运行、迭代次数内，IVNS算法能够用更少的运行次数和迭代次数获得质量和稳定性更高的结果，获得的最好解、最差解、平均解和标准方差均优于IPSO算法。平均解最大差距是272.944 4，减少百分比最大为17.83%，标准方差最大差距是38.71(如表8)。此外，从两个算法的不同算例的平均解和标准方差对比中可以看出(其中：AveS为平均解，StaD为标准方差)，IVNS算法的平均解和标准方差始终低于IPSO算法，且IVNS算法的标准方差的波动幅度更小，说明所设计的IVNS算法的求解质量和稳定性更高(如图19a和图19b)。

表8 与IPSO算法的求解稳定性比较

4.2.3 求解优越性比较

梳理现有相关文献，REED[9]、ABDULKADER[2]、KAABACHI[26]、陈久梅[3]等获得了目前最优解，将IVNS算法的求解结果与之对比分析发现，IVNS算法更新了28个算例中的23个目前最优解，优化值最高达到175.126 3，缩小了其余5个结果与目前最优解的差距，差距百分比仅在0.05%～8.98%之间，已经趋近目前最优解。进一步对比28个最好解的平均值，其中HABC算法的结果平均值为1 038.62，为历史最优，IVNS算法的结果平均值为980.225 6，相对于前者减少58.394 4，优化百分比为5.62%(如表9)。此外，观察不同算法的不同算例的历史最优解对比和IVNS算法相对于历史最优解的优化程度(其中：O为最优解)，仅有vrpnc1-S2、vrpnc1-S3、vrpnc5-S2、vrpnc6-S3、vrpnc12-S3算例未实现优化，其中vrpnc1-S3和vrpnc6-S3算例的未优化程度较大，为-8.98%和-8.76%，而优化算例中vrpnc13-S2的优化程度最大，为13.68%，且所有算例整体的优化程度为正，即整体上实现了优化效果，说明IVNS算法的求解质量更高(如图20所示，其中HVANS和HABC算法的vrpnc14-S3算例缺失数据按照HAC算法结果填充)。可见，相对于已获得历史较优解的IPSO、ACS、HAC、HVANS、HABC等算法，IVNS算法具有优越性。

表9 与历史最优解的求解质量比较

续表9

4.3 算法收敛性分析

以vrpnc2-S2算例为基准，采用其5次运行的迭代结果来展示所设计IVNS算法的收敛性(如图21)。结果显示，随着迭代次数的增加，求解结果逐渐收敛，并且在求解该较大规模算例的情况下，算法并不会过早或过晚收敛，收敛性较好。

5 结束语

本文通过设计针对MCVRP的IVNS算法并展开算例分析，验证了IVNS算法的有效性、优越性和收敛性。首先，采用扫描法产生初始解使得结果质量更优；运用插入、交换、2-opt*算子设计Shaking过程，在Shaking过程之后进行全局扰动，设计子路径内和子路径间搜索优化机制以实现Local Search过程，以及运用2-opt*、插入、交换算子设计子路径间搜索优化机制，均提高了IVNS算法的有效性。其次，IVNS算法具有优越性，能够更合理地规划网络结构、更稳定地获得高质量解，且更新了23个(28个算例)目前最优解、缩小了其余5个结果与目前最优解的差距。最后，IVNS算法的收敛性较好，不会过早、过晚收敛。由于篇幅限制，本文未将所设计IVNS算法用于实际MCVRP案例研究中，未来的研究应针对实际的MCVRP进行深入探讨，为实际应用提供科学指导。