多目标多时间窗车辆路径问题的鸽群-水滴算法

马 龙，王春嬉，张正义，董 睿

西安航空学院 经济管理学院，西安710077

车辆路径规划问题是在时间窗约束下实现运输路径和作业成本控制的组合优化问题，该问题是由学者Dantzig与Ramser于1959年提出，并在物流运筹学领域获得广泛深入的研究。目前对于车辆路径规划主要以最小运输总成本、最短运输距离和最少配送车辆数目为独立研究目标，或以这三个独立研究目标进行两两组合，作为优化计算的目标，这样会使解算目标与实际问题脱节，且问题的求解结果差异较大，特别是针对规模较大的求解问题，传统的运筹学方法在求解车辆路径规划问题最优解时耗时费力[1-2]，因此，国内外学者开始使用基本的人工智能算法或两种基本算法的混合方法，包括下等双向搜索算法[3]、禁忌搜索算法[4-5]、量子烟花进化算法[6]、蚁群算法[7-15]、遗传算法[16-21]、蝙蝠算法[22-23]、狼群算法[24]、混合进化算法[25-27]等。但这些基本算法或混合算法主要以改进基本算法的参数或两种算法的混合形式，对带时间窗的车辆路径问题进行求解。然而，改进后的算法依然存在早熟或无法完全收敛到全局最优解的问题。因此，探索使用互补性较强的两种智能优化算法，对两种算法的性能进行互补改进，形成一种新型改进算法是求解带时间窗车辆路径问题的热点研究方向。

智能水滴算法（Intelligent Water Drops，IWD）是学者Hosseini在2007年提出的一种智能算法，该算法通过水滴算子和携带泥土量算子来模拟自然界的水系统与河道影响环境而形成的河水迭代流动计算过程。该算法起初在旅行商（Traveling Salesman Problem，TSP）问题[28]、多维背包问题[29]等方面得到应用，随后，Kamkar等学者率先将该算法应用于车辆路径问题[30]，且部分学者对该算法进一步改进应用[31-34]，然而，面对多目标多时间窗车辆路径问题，智能水滴算法及其改进算法依然存在计算速度慢、全局收敛性能差等问题，因此，能否利用一种具备精准搜索方向的智能优化算法对智能水滴算法进行彻底改善，形成一种具有优势互补的新型智能优化算法解决车辆路径问题还鲜有研究。

鸽群算法（Pigeon Inspired Optimization，PIO）是受到鸽子的归巢行为启发，由学者段海滨等于2014 年提出的新型精准搜索算法，该算法是由两个精准搜索算子（地图罗盘算子和地标算子）分阶段独立运行，弥补具有单一算子或复合算子的智能算法搜索不足问题，该算法具有明确的搜索方向和快速的搜索速度、计算精度，且在理论探索和实际应用方面的成果较多。如背包问题[35]、无人机路径规划[36]、基本智能算法改进与函数优化求解[37]等。但将基本鸽群算法或单独利用鸽群算子完全改进智能水滴算法，解决带有多目标多时间窗的车辆路径问题的成果还未见报道。

综上可知，虽然现有的人工智能算法在改进和车辆路径问题应用方面出现大量的成果，但能使算法在求解多时间窗车辆路径问题时，能够满足收敛速度、计算精度、均衡开发和探索能力的还不够完善。因此，综合考虑具有硬时间窗约束的车辆路径规划问题的特点，构建多目标多时间窗车辆路径规划问题（Vehicle Routing Problem of Multiple-Objectives and Multiple Time Windows，VRPMOMTW）模型，并针对基本智能水滴算法的不足之处，提出求解VRPMOMTW 问题的鸽群-智能水滴互补改进算法（Complementary Pigeon-Inspired optimization and Intelligent Water Drops，CPIIWD），并利用benchmark 函数对CPIIWD 算法进行测试，且将仿真结果与人工免疫系统算法（Artificial Immune System，AIS）、粒子群算法（Particle Swarm Optimization，PSO）、遗传算法（Genetic Algorithm，GA）的结果进行比较，表明CPIIWD 算法收敛性能要比其他三种算法的性能更好，也是求解VRPMOMTW问题的最佳方法。同时，将VRPMOMTW模型求解结果与遗传算法[38]、智能水滴算法[33]等结果进行比较，有效弥补了智能水滴算法在多时间窗车辆路径问题的应用不足，从而强化了人工智能优化算法在车辆路径规划方面的应用效果。

1 多目标多时间窗车辆路径规划模型

1.1 问题描述

多目标多时间窗路径规划问题是通过给一个物流配送中心分配若干辆配送车，为多个服务点提供配送服务，根据不同服务点上的供货需求量、车辆的最大载重量、服务点与配送中心之间的路径长度等制约因素，所有服务点上均有相互独立的服务时间窗。配送车辆为多个服务点提供货物供应服务，从配送中心出发再回到配送中心形成完整的闭合回路。假设一个服务点由一辆车在规定的时间窗内完成配送服务，且为服务点提供装卸货服务的时间窗固定，使用每一辆配送车辆的成本和单位行驶距离成本均固定，单辆车的行驶总距离和行驶时间要求小于给定路线的最长距离和行驶时间。如何使配送总费用最少、配送路径最短和配送车辆数最少，从而使得车辆配送总成本达到整体最小。为了建模需要，作出如下假设：

（1）车辆由配送中心出发后返回到配送中心。

（2）采用相同型号的配送车辆，且车辆与服务点之间具有一对一和一对多的关系，其服务点需求量小于车辆最大载重量，服务时间恒定。

（3）运输道路状况和车流量忽略不计。

（4）配送中心分配的车辆数小于车辆总数，不得重复调度。

（5）车辆到达服务点的时间应在允许的时间窗范围内波动。

（6）模型的适应度值是以配送车辆的行驶距离和成本为衡量标准。

1.2 变量符号定义

为了准确理解数学模型描述的抽象问题，需要对模型中涉及的变量符号进行定义，并说明符号的含义，如表1所示。

1.3 模型建立

1.3.1 多目标函数

多目标多时间窗车辆路径规划问题以最小运输总成本、最短运输距离和最少配送车辆数目为优化目标，构建的数学表达式为：

式中，z1中第一项和第二项表示m 辆车的固定运输成本和行驶距离的成本，第三项表示车辆为服务点提供服务的惩罚成本，即由服务点预订时间段a 与开始服务时间sik之差构成，M 表示车辆提前与滞后到达的惩罚系数；z2表示车辆k 从第i 服务点出发到达第j 个服务点提供服务的最短运输距离；z3表示车辆k 从第i 到第j个服务点配备的最少车辆数。

1.3.2 约束条件

（1）车辆出发地点。为了安排最少的车辆数配送货物，要求每一辆货车必须从物流配送中心出发，建立的数学表达式为：

（2）车辆返回配送中心。为了充分利用安排的配送车辆，要求完成配送任务的车辆必须返回配送中心，建立的数学表达式为：

（3）车辆出发的初始地点。配送车辆的出发地点从0点开始，且不作为车辆的返回地点，建立的数学表达式为：

（4）车辆返回地点。配送车辆的最终返回地点需要额外增加返回地点，而不能以原来的出发地点作为返回地点，建立的数学表达式为：

（5）车辆服务点。配送车辆中必须存在唯一车辆与之需求点提供对应的服务，建立的数学表达式为：

（6）车辆到达服务点。如果恰好存在第k 辆到达第j 个服务点，则该车辆必须为该需求点提供服务，建立的数学表达式为：

表1 变量符号与含义说明

（7）车辆服务后离开服务点。如果为第i 个需求点配备了第k 辆车，则该车辆服务完成后必须离开第i 个需求点，建立的数学表达式为：

（8）需求总量低于最大装载量。在物流配送过程中，任何需求点的需求总量必须小于其车辆的最大装载量，建立的数学表达式为：

（9）车辆行驶距离限制。在物流配送过程中，任何服务车辆的总里程小于等于最大行驶里程，建立的数学表达式为：

（10）车辆出发时刻限制。在物流配送过程中，要求配送车辆从配送中心点出发的时刻为零，建立的数学表达式为：

（11）车辆到达服务点的时间关系。在物流配送路径上，任何一辆配送车相继到达两个需求点的时刻存在一定的关联关系，建立的数学表达式为：

（12）车辆服务时间窗限制。在物流配送过程中，为第i 个需求点在给定的时间窗α 内配备了第k 辆车，则该车辆到达第i 个需求点的时刻需要在时间窗α 的范围内，建立的数学表达式为：

（13）车辆服务完成时间限制。在物流配送过程中，给定第i 个需求点的第k 辆车需要在给定的时间完成服务，建立的数学表达式为：

（14）车辆不到达服务点。如果第i 个需求点没有安排第k 辆车为其提供服务，则该车辆无需到达第i 个需求点，建立的数学表达式为：

（15）决策变量的取值范围为：

2 多目标鸽群-水滴算法

鸽群算法作为一种典型的仿生进化算法，通过使用两类独立搜索算子模拟鸽子的不同搜索方向。在地图罗盘算子中，如果鸽群偏离目标距离时，经过地磁场形成的脑地图感知目标方向，靠近目的地后，地磁场与太阳系的作用减弱。如果鸽群临近目标距离时，熟悉地标的鸽子主要依靠地标算子指引飞行方向，其他鸽子跟随飞行达到目标地。关于基本鸽群算法的原理可参考文献[39]。智能水滴算法是在2007 年，学者Hosseini 根据自然界河流中的水流形态和环境变化模拟而提出的智能优化算法，该算法具有水流速度和携带泥土量两类算子，关于水滴算法的原理可参考文献[28-29]。下面只对利用鸽群算子改进智能水滴算法中，涉及的最优解选择与存储机制、适应度函数、改进的水滴流动速度、流动方向和携带泥土量以及改进后的智能水滴算法的步骤等进行详细介绍。

2.1 最优解选择与存储机制

利用鸽群-水滴算法求解多目标车辆问题时，算法依然是局部和全局搜索的组合过程。因此，首要解决如何平衡局部与全局搜索获得的Pareto非最优解的关系，如图1 所示。根据鸽群经过地图罗盘算子的方向导引和地标算子的鸽群折半择优方法获得每个鸽子的局部最优搜索，搜索到的局部最优解作为水滴流动方向和携带泥土量的子代个体。利用存储池结构来存储局部搜索到的Pareto非最优解。通过比较存储池内的解集，如果算法在搜索过程中得到了新的Pareto非最优解，则用其更新存储池内的一个劣势解，局部搜索过程完成后，对于智能水滴算法和鸽群算法以及存储池内所有的解使用Pareto排序和拥挤度选择机制[40]进行选择操作，获得新的迭代搜索群体。

图1 最优解选择与存储机制示意图

2.2 适应度函数

根据鸽群-水滴互补优化算法的局部和全局搜索机制的差异，在同时考虑三个目标函数时，采用三种不同的适应度函数，分别用于全局最优解和局部搜索最优解的选择操作。

算法的选择操作采用Pareto 排序和拥挤度的度量方法[40]，fiti( x )=( rank(x),Crowding-distnce(x) )，其中，rank(x)为Pareto 占优的层级关系，rank(x)＜rank(y)为解x 占优y 。Crowding-distnce(x)为拥挤度距离，距离越大越好。在选择操作时，对于鸽子群体、水滴群体和存储池中的所有解按照第一个关键字rank(x)升序、第二个关键字Crowding-distnce(x)降序排序，然后选择前N 种群个体迭代进入下一代计算。

在鸽群算法的地图罗盘算子和地标算子方向的导引下，进行局部搜索，每次搜索均是从候选解的邻域解列表中选择一个最好解更新当前解，然后继续搜索其邻域，因此可引入Pareto 强度来定义局部搜索的适应度函数[41]：

式中，S( x )表示在存储池中被x 占优的解的个数，W( x)表示在存储池中占优x 的解的个数，局部适应度函数反映了解x 在Pareto最优解集中占优强弱程度，如果解x强度大，则适应度高，说明解x 在Pareto 最优解集中处于优势地位，可在局部搜索过程中优先搜索具有优势地位解的邻域。

2.3 利用鸽群搜索算子改进智能水滴算法

智能水滴算法中的水滴通常被抽象为离散化形式，这是因为每一滴水在运动过程携带的泥土量不同，且会以携带的泥土量大小来选择流动的方向，这种状态不利于河道内的水滴快速向最优解靠拢，也不利于携带泥土量多的水滴向最优解位置聚集。然而，离散二进制粒子群算法[42]以及利用改进鸽群搜索算子优化的粒子群算法[37]，对于求解种群个体的收敛速度和位置具有良好的效果，本文充分利用离散二进制的位置变化概率和鸽群中的两类导航算子对智能水滴算法中的水滴流动速度、流动方向和携带泥土量算子进行改进，下面分别对智能水滴算法的改进过程进行阐释。

2.3.1 改进的水滴流速算子

Kennedy与Eberhart提出的离散二进制粒子群算法是对粒子个体进行二进制编码后，将粒子飞行速度转化为变换概率[43]，并取得了良好的收敛效果。为了表示水滴速度值为二进制位取1的概率，本文将水滴流动速度值映射到区间[0，1]，映射的方法一般利用Sigmoid 函数，其数学表达式为：

式中，S( velocityij(IWD) )为水滴移动位置x( i,j )取1 的概率，水滴通过式（19）改变对应位值：

式中，rand( )表示在[0 ,1] 区间内产生的随机数，为了避免S( velocityij(IWD) )过度趋近0 或1 的边缘，使用参数velocitymax表达为水滴流动的最大速度值，用于限制水滴流动速度velocityij(IWD) 的取值范围，即velocityij(IWD)∈[- velocitymax,velocitymax] ，水滴流动速度最终限制了x( i,j )取0或1的概率。

智能水滴算法中的水滴经过离散化处理后，水滴流速快慢与携带的泥土量多少直接相关，特别是水滴流动过程中席卷的泥土量越大，则水滴的流速必会减慢，因此，利用鸽群算法中的地图罗盘算子对水滴的流速进行改进，改进后的数学表达式为：

2.3.2 改进的水滴流动方向

为了说明离散化水滴携带大量泥土而选择最优的流动方向，水滴位置x( )i,j 变化概率是关键，根据文献[43]中提出的位变化概率，对此进行分析，根据式（19），水滴流动轨迹是一种根据携带泥土重量交替选择概率，设为水滴在第t 次从位置i 流动到位置j 的速度，第t 次位为1的概率是，而且位为0的概率是1，此外，如果第t-1 次位值已成为0，则第t 次变化的概率为，同样，如果第t-1 次位值已成为1，则第t 次变化的概率为，而第t-1 次位值是0 的概率为，第t-1 次位值是1 的概率为。则第t 次位的变化概率p( t )为：

根据式（18）的结果，可得到如下的数学表达式：

根据水滴携带泥土重量的大小，选择流动概率，但流动方向还与地标指引方向有关，因此，利用鸽群算法中的地标算子改进水滴流动方向，其数学表达式为：

式中，xcenter为引导水滴流动的中心位置；其余变量的符合含义与式（20）一致。

2.3.3 自适应变邻域扰动携带泥土量算子

在以上改进的IWD 算子中，只是对水滴流经河道内的水滴速度和方向进行改进，并未对整个河道内的泥土量变化情况进行考虑，为了提高算法开发和探索能力，提出自适应变邻域扰动策略对水滴携带泥土量进行干扰，增加算法的全局收敛能力。自适应变邻域扰过程，如图2所示。

图2 自适应变邻域扰动过程

根据图2可知，在IWD水滴算法中，水滴位置S1 和S2 向目的地汇聚过程中，河道中央的水滴会不断对邻域位置上的水滴进行扰动，且会根据河道中的泥土量Soil(IWD)取值实时更新，并根据河道邻域内的水滴位置k 与汇聚点D 的间距大小来更新，如果间距小，则无需更新河岸edge( )k,j 的泥土量，否则位置k 上的水滴将产生新携带的泥土量Soil(IWD′)，该水滴个体携带的泥土量Soil(IWD)初始值设定为0，水滴流速velocitykD(IWD′)的初始值设定为Initvelocity，其携带新泥土量的水滴速度和位置的更新表达式为：

由此可推算出河岸edge( k,j )上的泥土量Soil(IWD)参数变化表达式为：

2.4 改进的智能水滴算法步骤

标准智能水滴算法的两类算子在迭代计算时，通常根据水滴流动过程中携带的泥土量来确定水流速度和移动方向，且携带的泥土量也会随着水流速度发生变化，这样算法在局部空间搜索计算时，由于某段河道内水滴携带泥土量较多，算法会陷入早熟收敛状态，然而鸽群算法中的两类独立运行算子能够根据地图罗盘和地标的导引，可较好地指明个体鸽子的飞行方向和飞行速度，因此，本文利用鸽群算法中的地图罗盘算子和地标算子以及自适应扰动策略分别对智能水滴算法中的流动速度、流动方向以及水滴携带泥土量进行改进，设计出引入鸽群搜索算子的智能水滴算法，较好地克服标准智能水滴算法在车辆路径问题求解方法的不足，算法的具体计算步骤如下：

步骤1 设置算法参数；最大迭代计算次数tmax；鸽子个数m；水滴个数N ；水滴流速参数av、bv、cv；泥土更新参数as、bs、cs；局部泥土更新系数ρ0、ρn；起止位置的初始泥土量

步骤2 初始化所有水滴的参数；设置算法的动态参数；水滴初始状态集合S( )IWD =∅，水滴初始速度、水滴初始携带泥土量Soil(IWD)以及初始携带重量。

步骤3 选择水滴局部更新信息；水滴根据流动路径上的泥土量大小，按照式（21）概率大小选择流经下一节点，根据被选择节点更新水滴的局部信息，按照式（24）～式（26）更新水滴的局部信息。

步骤4 根据式（24）、式（26）分别对水滴的流动速度和携带的泥土量进行更新。

步骤5 条件判断；如果河道内所有的水滴都会收敛至同一路径或达到最大迭代计算次数，转到步骤6，否则转到步骤3。

步骤6 保存最优解与最优值，求解计算结束。

3 多目标多时间窗车辆路径模型处理方法

多目标优化问题通常采用惩罚函数与帕累托最优解集进行约束条件处理，且这两种方法有着不同的特点。下面从惩罚函数和帕累托最优解集两个方面，对本文的多目标车辆路径模型的处理方法进行探讨分析，并提出多目标优化与惩罚函数的混合策略。

3.1 模型处理策略

多时间窗车辆路径规划问题可使用有向有环图来表达，假设图G=(V ,E )，V={0 ,1,2,…,n} 表示车辆配送中心0 与服务点i=1,2,…,n，E 表示车辆节点之间的路径集合。同时，使用水滴抽象表示每一辆车，水滴多次经过配送中心0，一次经过服务点，形成完整的流动路径过程恰好为VRPMOMTW的解。求解一次迭代结束后，可在每个水滴流经的完整路径中搜索出最优路径，利用该最优路径不断更新各个水滴节点之间的泥土量与全局最优路径，然后迭代到下一步迭代过程，直到达到算法的最大次数或期望的优化路径。

3.1.1 多目标函数的惩罚处理策略

多目标车辆路径规划问题涉及三个以上的目标函数和多类复杂的约束条件，该类问题的求解较为复杂，为了简化问题的求解难度，利用罚函数法[44-45]和理想点方法[46]分别对模型中的约束条件与多目标函数进行简化处理。根据理想点法，将最小运输总成本、最短运输距离和最少配送车辆数目的多目标转化为单目标问题，多目标车辆路径规划模型可转化为以下形式：

式中，Q1为A( x )的效用函数，取值为利用线性规划方法求解该问题的上界值；Q2为B( x )的效用函数；Q3为C( x )的效用函数，由约束违反度函数的特性，且取值为零；λi为第i 目标函数的权重，且∑λi=1，这种方法只能找到部分Pareto 最优解，为了找到更多的Pareto 最优解，可利用Pareto排序法确定各目标函数的权重系数。

3.1.2 多目标优化与惩罚函数的混合策略

通常采用惩罚函数法和Pareto 最优解集对多目标函数问题进行处理，但根据上述惩罚函数处理方法可知，虽然惩罚函数的选取具有一定的主观性，但根据Pareto 优超的定义和关系[47]，Pareto 最优解集并非适用于所有个体，根据需要引入偏好程度才更为有效，如图3所示，在Pareto 最优解集中包含3 个依次排列的非劣个体A、B、C，假设只有个体A 距离可行域较近，而个体B与C 距离可行域较远，所以A 是需要提取的重要信息。另外，不可行解是无法超越Pareto 可行解，这是因为全局最优点落在边界点上，且算法进入搜索后期会丢弃不可行解。因此，基于Pareto最优解集的多目标优化方法往往忽视了偏好程度，且引入偏好的简便方法是惩罚函数法。因此，本文提出多目标优化与惩罚函数混合方法，首先为了权衡目标函数与约束条件方程G(x)的标准变化，将从Pareto 最优解集中选择个体的搜索偏好，找出群体中的最大和最小目标函数值，然后对每个目标函数值和约束条件的程度进行归一化处理，这样处理的目标函数具有降维特征，且可以比较Pareto最优解集中的个体，根据偏好选择个体。同时，在算法搜索初期存在少量可行解，通过引入惩罚函数，可以引导算法尽早进入可行域。另外，优化过程中，惩罚系数会依据群体状态自动调整，当群体的可行解比例p ＞0.5 与惩罚系数ρ ＜1 时，可开发出目标函数与约束程度较小的不可行解。

图3 Pareto最优解集中的非劣解示意图

3.1.3 权重系数处理策略

在问题可行域上对各个目标函数zi( x )进行极小化处理，计算极小值点xi，得到：

利用算术平均法计算i 个分目标函数的极小点的相对离差、为了避免各个分目标函数值离差过大而影响权重系数的选取，对求解的相对离差做无量纲化处理。

同时，利用算术平均法计算出第i 个分目标函数的极小点的平均相对离差，得到：

显而易见，为了使权重系数规范化，当Δi＞0 时，假设Δi1≥Δi2≥…≥Δis，且Δis∈{ Δ1,Δ2,…,Δm} ，选取如下式作为各分目标对应的权系数：

根据上述分目标函数与权重系数的确定方法，构建的多目标车辆路径规划模型的数学表达式为：

式中，变量的含义与式（1）、式（27）的一致。

3.1.4 约束条件处理策略

多目标窗车辆路径规划问题是带有多种约束条件的复杂问题，采用罚函数法[44-45]将模型中的部分约束优化问题处理为无约束优化问题，由此对式（2）～式（15）中的约束条件处理如下：

式中，φi( x )为不等式约束条件；θi( x )为等式约束条件；σ 表示等式约束的惩罚因子；表明优化后的车辆配送成本应允许在[- σ,σ ]范围内波动，通常σ 取很小的值。

3.2 模型求解步骤

步骤1 输入初始静态变量；服务点与配送中心集合D；服务点个数n；配送中心点0，服务点i 的需求量为qi,多时间窗t=1,2,…,p；服务时间sti,i=1,2,…,n；服务中心至服务点的距离矩阵为d；车辆的单位运行成本c；车辆固定成本g；车辆行驶速度v；车辆的最大载重量Qmax；车辆在服务点与配送中心之间的行驶时间矩阵T ；最大迭代次数tmax；水滴数N ，水滴流动更新速度参数av、bv、cv，泥土携带量更新参数as、bs、cs，局部泥土量更新权值a，全局泥土量更新权值β，服务点上的初始泥土量InitSoil,水滴的初始流动速度InitVelocity,携带泥土量初始矩阵

步骤2 输入初始动态变量；初始化已访问节点集Visitnode 为空集；未访问节点集初始化为Novisitnode={0 ,1,2,…,n} ；水滴（车辆）从配送中心出发携带泥土量初始化为InitSoil=0；水滴（车辆）从配送中心出发的流动速度初始化为InitVelocity=0；水滴（车辆）从配送中心出发的载重量初始化为Qmax=0；水滴从配送中心出发的开始时间为t=0。

步骤3 根据路径规划模型的约束条件，随机选取水滴的初始访问节点，以时间窗动态更新水滴访问节点列表Visitnode( IWD )和未访问节点列表Novisitnode( IWD )。

步骤4 根据时间窗上未访问节点Novisitnodei的需求量，确定水滴移动访问节点集合DV ；如果未访问节点集合NDV 不符合容量和时间窗约束，水滴（车辆）返回配送中心0，初始化水滴对应的动态变量，形成一辆配送车辆的路径；如果未访问节点集合NDV=∅，水滴（车辆）返回配送中心0，形成水滴（车辆）的完整访问路径，该路径对应VRPMOMTW问题的可行解TIWD。

步骤5 利用多种节点选择方法，在水滴到达节点的服务概率与距离以及泥土量中加入重要度因子，计算水滴访问列表的概率p( t )，确定从当前节点i 到NDV 中的下一个服务点j 的概率，如下式：

εs为较小的整数，且εs=0.01。

步骤6 根据访问节点集合DV 中每个节点对应的概率，采用式（22）选择下一个被访问节点。并将被访问服务点j 列入已访问节点集合中，并修改当前水滴对应的已访问节点集合、访问顺序和未访问节点集合。

步骤7 根据式（24）、式（25）、式（26），分别对水滴（车辆）的流动速度velocityij(IWD) 、携带泥土量Soil(IWD)和路径滚动泥土量进行更新；其中，为了避免算法早熟问题，路径泥土携带量需要在Soil( i,j )∈[min Soil( i,j ),max Soil( i,j )]之间变动，且min Soil( i,j )与max Soil( i,j )的计算表达式如下[34]：

式中，Pbest=0.5，n 表示服务点数；表示当前局部最优解。

步骤8 根据步骤3 中计算的完整访问路径TIWD，并假设运输总成本cost(TIWD)，则搜索出目标函数中最小平均可行解作为模型的最优解TTB。

步骤9 根据当前迭代计算得到的最优解TTB，并将更新后的泥土量矩阵作为下次迭代计算的初始泥土量矩阵，其计算表达式为：

其中，nIB表示当前迭代计算路径上的服务节点数，表示最优路径上水滴（车辆）携带的泥土量。

步骤10 更新模型计算的全局最优解：

步骤11 更新迭代计算次数ti=ti+1，若ti＜tmax，则转入步骤2，否则，如果ti=tmax，则将最优解fmin作为VRPMOMTW的最优解，并输出计算结果。

4 实例仿真与结果分析

4.1 算法性能测试函数

为了验证所提CPIIWD 算法的性能，以文献[37]中的6 个经典函数为例，对算法性能进行仿真测试，其函数的具体形式如下：

（1）Sphere"s function

（2）Rosenbrock"s function

（3）Griewank"s function

（4）Ackley"s function

（5）Penalized function

（6）Shifted Griewank

4.2 实例数据来源

为了进一步验证本文提出的VRPMOMTW 的稳定性与算法的可行性，选用文献[33，38]中的案例，利用CPIIWD 算法对VRPMOMTW 进行优化求解，并与文献[33，38]中的求解结果进行比较。

实例1 某物流配送中心在不同时间窗内和服务时间长度，为10个需求量不同的客户提供物流配送服务，其配送中心与需求服务点的距离和服务时间如表2、表3所示。

表2 配送中心与服务点之间的距离 km

表3 多时间窗的服务点需求量与服务时间

实例2 某配送中心的车型相同，配送服务需求点为20，配送中心位置为（0，0），每个需求点的位置坐标、需求服务量和时间窗等信息如表4所示。

4.3 模型与算法参数设置

实验环境：Inter CoreTMi5-2450MCPU@ 3.2 GHz，内存为4 GB，Window7 操作系统，MatlabR2015a 版本，四种算法经过反复迭代计算后，选取如下参数后，计算结果较为理想。

表4 服务点位置、需求量、服务时间、时间窗数据

算法参数：（1）人工免疫系统算法[48]。抗体数n=30,交叉概率pc=0.6，变异程度pm=0.3，淘汰率pe=0.3，浓度阈值Tc=0.3，抗体亲和力阈值Tac1=0.75，记忆细胞亲和力阈值Tac2=0.8，记忆细胞数量Nm=5，最大迭代次数tmax=500。

（2）粒子群算法[37]。种群规模S=300,K=3,维度D=100,最大迭代次数tmax=500。

（3）遗传算法。种群规模N=300 ，交叉概率pc=0.7，变异概率pm=0.3，最大迭代次数tmax=500。

（4）鸽群-水滴算法。水滴数NIWD=100，水滴流动更新速度av=1,bv=0.1,cv=1，泥土携带量更新参数as=1,bs=0.1,cs=1，局部泥土量更新参数a=1，全局泥土量更新参数β=1，服务点上的初始泥土量InitSoil=2 000，初始流动速度InitVelocity=100，最大迭代次数tmax=500，水滴流动方向的罗盘导航因子R=0.9，局部泥土更新系数ρ0=ρn=0.5。

模型参数：车辆最大载重量Qmax=90 t ，车辆固定成本g=50 元，车辆行驶单位成本ck=20 元，车辆的平均行驶速度v=50 km/h，车辆在配送中心与服务点之间的行驶时间tij=dij/50 h。

实例参数：实例1中，水滴数NIWD=100，水滴流动更新速度av=1,bv=0.1,cv=1，泥土携带量更新参数as=1,bs=0.1,cs=1，局部泥土量更新参数a=1，全局泥土量更新参数β=1 ，服务点上的初始泥土量InitSoil=2 000 ，水滴的初始流动速度InitVelocity=100，最大迭代次数tmax=100，水滴流动方向的罗盘导航因子R=0.9，局部泥土量更新系数ρ0=ρn=0.5。

实例2 中，水滴数NIWD=200，水滴流动更新速度av=1 000,bv=0.1,cv=1 ，泥土携带量更新参数as=1 000,bs=0.1,cs=1，局部泥土量更新参数a=0.9，全局泥土量更新参数β=0.9 ，服务点上的初始泥土量InitSoil=2 000 ，水滴的初始流动速度InitVelocity=100，最大迭代次数tmax=800。

4.4 仿真结果与分析

4.4.1 算法性能测试结果与分析

根据表5 的仿真结果可知，四种算法分别经过500次的迭代测试后，发现除了Rosebrock 函数和Penalized函数无法收敛到理想值，其余4种函数几乎可收敛到理想的函数值。虽然AIS算法的收敛计算时间明显减少，但AIS算法的收敛能力较弱，这是因为AIS算法中保存在记忆库中的抗体数较少时，降低了记忆细胞的搜索能力。另外，经过Penalized 函数的仿真结果发现，虽然GA和CPIIWD收敛计算的结果几乎接近理想值，但GA算法的求解速度慢于AIS算法的求解速度，这是因为该函数是分段函数，算法在收敛计算时的参数设置的不同而造成的差异。同时，使用PSO算法仿真测试的多维单峰和多峰函数收敛计算时间较长，而采用CPIIWD算法的收敛计算时间明显减少，收敛效果要优于其他3种不同的算法，这说明鸽群-水滴互补优化算法在多维函数优化方面的改进效果要比单个人工智能算法的效果更好。

表5 四种算法的计算结果

图4 四种算法的仿真收敛曲线

根据图4的收敛曲线清晰看出，提出的鸽群-水滴互补优化算法由于引入精准的鸽群算子和多目标种群初始化方法，使得智能水滴算法在计算复杂多维函数时，可逃逸局部最优搜索空间，并在短时间内快速向全局最优解方向收敛。特别是AIS算法在6种不同函数上的收敛结果都不理想，且计算时间要比CPIIWD算法的计算时间长，但从部分函数也可看出，CPIIWD 算法在某些测试函数上也不能完全表现出最优性能，这是因为对水滴算法的初始参数选取和固有的弊端总会导致部分搜索计算效果较差，但从总体上来看，CPIIWD 算法要明显优于其他3种算法，这是因为利用鸽群搜索算子对水滴算法中的算子互补改进后产生的明显效果。

4.4.2 模型求解结果与分析

通过算例1 中的数据源与处理后的约束条件（33）和目标函数（1）的编程，利用CPIIWD 算法求解出本文最小运输总成本为6 705元、最少配送车辆数为3辆，最短配送路径如图5所示。

图5 3辆车的最短配送路径示意图

根据图5 的车辆配送路径分布结果可知，实例1 中的车辆从物流配送中心出发，经过不同时间点和服务点后，3辆配送车行经最短运输距离，最终返回配送中心，符合物流车辆配送作业的基本要求。

为了进一步验证本文CPIIWD算法求解VRPMOMTW问题的可行性，分别利用遗传算法、基本智能水滴算法和CPIIWD算法对车辆路径问题进行解算，在实例1中，三种不同的人工智能算法经过100次的迭代计算后，其收敛计算结果的曲线如图6所示。

根据图6中的收敛曲线可知，经过三个不同的人工智能算法仿真计算后，遗传算法经过23 次迭代后获得了问题的全局最优解，基本智能水滴算法经过19 次得到了全局最优解，而CPIIWD算法经过16次获得了问题最优解，这说明经过鸽群搜索算子对水滴算法的互补改进后，模型的求解效果要明显优于其他两种算法，且使用该算法求解的最小运输总成本要比其他两种算法少，为物流企业节约不少的运输费用。

表6 三种不同算法的求解计算时间比较

图6 三种算法的收敛计算过程

为了再次验证本文算法求解效果，通过对3种不同算法的寻优次数、最少迭代次数和计算时间等进行比较，如表6所示。

通过表6的计算结果可知，三种算法在相同的数据源和模型参数设置下，基本智能水滴算法的寻优次数为96，遗传算法的寻优次数为30，而CPIIWD算法为28，这是因为鸽群算子具备精准的罗盘算子和地标算子，使改进后的水滴算法的计算时间和迭代次数明显减少，但该算法得到的最差值要逊色于其他两种算法，说明根据文献[37]中的原理，任何算法不是万能的，不可能在各个方面都能表现出良好的效果。

通过实例2 中的数据来源与模型参数设置，利用CPIIWD 算法对多目标多时间窗车辆路径问题进行仿真计算，其20个服务点的最优配送路径示意图如图7所示，最小成本的收敛计算结果如图8 所示，并对20 个服务点的车辆路径距离进行优化计算，如图9所示。

图7 20个服务点的最优配送路径示意图

根据图7中的配送路径分布结构图可知，物流服务车辆经过不同的需求服务点后返回服务中心，且经过CPIIWD算法的求解计算后，20个服务点上需要最少分配4辆相同型号的车辆。

图8 鸽群-智能水滴算法的收敛计算过程

图9 20个服务点的车辆路径最优化距离

通过表7的计算结果可知，经过使用CPIIWD算法的优化计算后，20 个需求服务点上，使用4 辆运输车辆的最小运输总成本为1 765.836元，这与文献[38]中使用的基本智能水滴算法求解计算的最小运输总成本2 168.9元、最短运输距离353.804 km相比，每条运输路径的拓扑结构不同，路径距离分别减少20 km 左右，这说明使用精准搜索的鸽群算法改进水滴算法，求解VRPMOMTW问题的结果更为精确，优化效果更为明显。

表7 不同配送路径上车辆的载重量

根据图8 的收敛计算结果可知，为了增强CPIIWD算法求解VRPMOMTW问题的可信性，通过该算法800次迭代后，求解计算的最小运输成本达到最优，且该算法经过260 次的计算后，求解结果基本达到平稳状态，说明该算法能够快速收敛到问题的最优解。

由图9 的计算结果可知，使用CPIIWD 算法求解VRPMOMTW 中的路径距离后，车辆从配送中心出发后，经过20个需求服务点形成不同运输路径拓扑结构，最终返回配送中心，且计算获得的车辆路径最优化距离为273.191 1 km，这与文献[33]中使用基本智能水滴算法求解的车辆路径距离减少80 km左右，降低了车辆运输总成本。

5 结语

（1）针对带时间窗车辆路径组合规划目标的差异化问题，提出多目标多时间窗车辆路径规划模型，实现了物流车辆运输成本、运输距离和配置数量的抽象建模，拓展了物流车辆路径规划问题的研究范围。

（2）针对基本智能水滴算法求解车辆路径规划模型速度慢、结果精度低等问题，利用地图罗盘算子和地标算子分别改进了水滴的流动速度和方向，同时，利用自适应变邻域扰动策略干扰水滴携带的泥土量，提高水滴算法的搜索计算速度和求解精度。

（3）针对多目标模型的求解难度大、处理过程复杂等问题，采用理想点法和多目标优化与罚函数混合方法，分别处理多目标函数与复杂约束条件间的关系，简化了人工智能算法求解模型处理能力。

（4）通过利用经典标准测试函数与物流车辆路径规划问题作为仿真实例，测试了鸽群-水滴算法的优化性能以及验证了模型求解的可行性，并从寻优次数、收敛计算时间和平均值等指标考核了4 种不同算法的计算结果，证明了CPIIWD 算法求解车辆路径问题的优越性，虽然改进的智能水滴算法在车辆路径规划问题中获得较好的求解结果，但该算法能否在不同工程优化问题中表现出良好性能，这需要进一步对本文算法的深入应用展开研究。