混合蚁群算法求解双目标时间窗VRP

邓丽娟，张纪会

(青岛大学a.复杂性科学研究所；b.山东省工业控制技术重点实验室，山东 青岛 266071)

0 引言

车辆路径问题[1]，又称车辆调度问题，由Dantzig等学者于1959年首次提出。该问题的一般提法是：给定配送中心、客户的位置以及各客户的需求量，车辆从配送中心出发，完成配送任务后回到配送中心，确定合理的配送方案，在满足车辆装载能力和行程限制等约束条件下使得某些目标达到最优(如总行程最短、使用车辆数最少、费用最小等)。该问题具有广泛的实际应用背景，也是著名的NP难问题。根据应用背景、目标函数和约束条件的不同，该问题有很多变种。1985年Savelsbergh提出带时间窗的车辆路径问题[2],这类问题更加关注服务的时效性。Solomon将带有时间窗的车辆路径问题定义为：车辆从配送中心出发，在满足各个客户的服务时间约束下，将货物运送到客户手中，最后返回配送中心，使行驶的总行程最短[3]。在实际配送问题中，除了运输成本外，还需要考虑许多其它的要求，因此，现实中的车辆调度问题往往是多目标优化问题，求解难度更大。研究带有时间窗的多目标车辆路径问题具有重要意义。

众所周知，VRPTW是典型的NP难问题，还没有一种有效的算法能完全解决它。目前的主要求解方法可分为：精确算法和启发式算法。精确算法有动态规划法、分支定界法、两阶段法、最短路径算法等，Christofides[4]用动态规划求解VRP，Laporte和Nobert[5]提出了分支定界法，Fisher[6]提出两阶段法求解VRPTW，Azi[7]提出最短路径算法。只有小规模问题才能用精确算法在较短的时间内寻求其精确解。因此，设计有效的启发式算法来寻求大规模车辆路径问题的近似解是一个重要的研究方向。启发式算法包括经典启发式算法和智能启发式算法。经典启发式算法采用局部搜索寻求满意解，算法易实现但容易陷入局部最优。如Clarke和Wright[8]使用节约法求解VRP，适用于小规模的VRP问题，Glover[9]将禁忌搜索法用于VRP的求解，但其求解速度较慢。智能启发式算法采用全局搜索来寻求满意解，能避免陷入局部最优，但算法设计较为复杂，参数设置对算法影响较大。如近年来发展起来的现代优化算法：遗传算法[10-14]、粒子群算法[15-19]、蝙蝠算法[20]、蚁群算法[21-26]等越来越多地应用于VRP领域。然而在实际应用中，粒子群优化算法和人工神经网络在迭代后期易陷入局部最优，遗传算法初始种群敏感易早熟，蚁群算法虽然寻优能力较强却依赖于参数设置，故启发式优化算法的实用性、寻优效率及应用价值均有待提高。

大多数VRP问题都将成本作为主要优化目标，但许多基于成本优化的解决方案被证明是不可能在实践中应用的。因为如果没有其他性能度量进行平衡，只考虑成本最小可能产生反效果[27]。因此，本文设置的目标是以对物流企业来说达到总成本最小，对客户来说达到服务的满意度最高，以此实现企业与客户的双赢，有利于长期合作。

本文提出一种混合蚁群算法(HACO)来求解VRPTW。设计了更符合实际情况的宽松配送时间窗。利用精英蚂蚁策略分别探索两个目标函数，获取更好的非支配解，通过重新定义自适应挥发因子及加入变邻域搜索算法来平衡算法的局部和全局搜索能力。使用NSGA-Ⅱ获得Pareto解集。通过Solomon标准算例对算法性能进行测试,并对车速进行了灵敏度分析。仿真结果表明，HACO能找到质量更好的解，能有效减少物流配送成本，提高服务质量。

1 VRPTW的数学模型

本文所考虑的优化目标为：最小化总成本和最大化客户满意度。总成本包括车辆离开仓库的固定成本、车辆在行驶过程中产生的运输成本以及车辆未在时间窗内服务的惩罚成本。

1.1 基本假设

1)已知配送中心及客户的位置以及各客户的需求量。

2)已知各客户的配送时间窗。

3)配送车辆数充足，所有车辆的属性(油耗、载重、速度)一致。

4)每个客户的需求量不超过车辆运载能力，每个客户只能由一辆车配送。

5)每辆车所服务的客户需求之和不超过车辆的载重量。

1.2 满意度函数

在带时间窗的车辆路径问题中，时间窗往往设置为单时间窗，但由于物流企业需要降低成本，客户需要及时服务。为有效提高服务质量，本文设计了宽松时间窗，使之更符合实际情况。假定客户i的服务时间窗为[EETi,LLTi]，若车辆为客户i服务的时间在此时间窗外，则客户i的满意度为0，且会对物流企业带来不良影响，产生惩罚成本。客户i的最满意的服务时间窗为[ETi,LTi](其中，[ETi,LTi]⊆[EETi,LLTi])，若车辆在该时间窗内为客户i服务，则客户i的满意度为1；车辆为客户i服务的时间离该时间窗越远，客户i的满意度就越低。假设μi为配送过程中客户i的服务满意度，则客户满意度函数表示为:

(1)

其中，ti表示车辆到达客户i的时刻。若车辆早于EETi到达客户i，则需等待，并且产生等待成本。若晚于LLTi到达客户i，会产生惩罚成本。

1.3 惩罚函数

客户的最优配送时间窗由客户指定。配送车辆在配送时间窗[EETi,LLTi]内完成任务，不予处罚。若配送车辆早于或晚于时间窗配送，则影响服务质量，给予一定惩罚。惩罚成本函数设计如式(2)：

(2)

其中,EETi为客户i的最早服务时间，LLTi为客户i的最晚服务时间，ti为车辆到达客户i的实际时刻；p为单位时间等待成本，q为单位时间惩罚成本。

1.4 建模

1.4.1 参数和变量

节点集记为V={0,1,2，…,n}，其中节点0表示配送中心，是配送车辆的起点和终点，客户节点集为{1,2，…,n}。K={1,2,…,k}为配送车辆集合。其它参数定义为:dij：客户i和j之间的距离；qi：客户i的需求量；W：车辆的最大承载能力；C：车辆行驶单位距离产生的成本；tij：车辆从客户i行驶到客户j所需要的时间；fij：车辆从客户i到客户j的行驶成本；c：车辆离开配送中心的固定成本；Ti：客户i的服务时长；ti：车辆到达客户i的时刻；v：车辆的速度；Qk:车辆k的服务客户集；[EETi,LLTi]：客户i的服务时间窗；[ETi,LTi]：客户i的最满意服务时间窗；决策变量：

(3)

(4)

1.4.2 目标函数

基于以上假设，所研究的VRPTW的目标函数为

(5)

(6)

目标函数(5)表示最小化总成本，包括车辆离开仓库的固定成本、运输成本、车辆未在时间窗内送达的惩罚成本；目标函数(6)表示最大化客户的平均满意度。

1.4.3 约束条件

(7)

(8)

fij=Cdij∀i,j∈V

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

约束条件(7)表示满足车辆容量限制，每辆车不能超过其最大载重量；约束条件(8)表示车辆行驶时间与距离成正比；约束条件(9)表示车辆行驶成本与距离成正比；约束条件(10)表示每个客户只能由一辆车配送；约束条件(11)表示满足车辆数量限制；约束条件(12)表示服务时间限制，其中M是充分大的正数；约束条件(13)～(15)表示车辆从配送中心出发又回到配送中心，形成一个回路。约束(16)和(17)定义了决策变量的性质。

2 混合蚁群算法

蚁群算法由意大利学者Marco Dorigo提出，是用来探索优化路径的概率型算法。传统蚁群算法的鲁棒性较好，且具有优良的分布式计算能力，易与其它算法结合使用，但易陷入局部最优。为增强搜索能力，本文引入变邻域搜索算法，将NSGA-Ⅱ算法用于双目标的择优过程，同时将精英蚂蚁分为两部分，分别探索两个目标，达到同时优化双目标的目的；并将挥发因子设置为自适应变化。

2.1 Pareto解集

定义1：(最小化问题的Pareto最优)对于两个决策变量a,b∈X，称a支配b，记为a≻b，当且仅当∀i∈{1,2,3…,n},fi(a)≤fi(b),∃j∈{1,2,3…,n},fj(a)

2.1.1 非支配排序

为了解决双目标优化问题并保证Pareto解集的多样性，用NSGA-Ⅱ来处理双目标的择优过程。它具有以下两个特点:非支配排序和拥挤度设计。对每次迭代后得到的解进行非支配排序，解集中的每个不被其他解支配的解，记为第一层F1，其非支配序值irank为1；然后从所有解集中删除该层解，在其余解中找到非支配解集，记为第二层F2，其非支配序值irank为2；以此类推，直至每个解都分层完毕，其中，同一分层内的解的非支配序值irank相同。

2.1.2 拥挤度设计

拥挤度是指在解集中每个解周围所包含的解的密度，解的拥挤度是相同分层的两个相邻解的目标函数的距离差。拥挤度能衡量解的多样性。其算法描述如下：

1)令拥挤度nd=0,n∈1,2,…,n。

2)对每个目标函数fm：

当蚂蚁选择路径及进行局部搜索时，根据以下规则判断两个解的优劣：

1)当两个解不在同一非支配等级时，非支配序值irank较低的解更优。

2)当两个解在同一非支配等级时，拥挤度较大的解更优。即此解周围的其他解较少，所在区域解的分布较为分散。

2.2 算法流程

2.2.1 伪随机概率选择规则

为避免搜索陷入停滞，蚂蚁在选择下一个客户时采用确定性与随机性相结合的选择策略，引入伪随机因子q0，位于客户i的蚂蚁k按式(18)选择下一个客户。

(18)

其中，q为区间(0,1)内的均匀分布随机数,q0为可调参数且q0∈[0,1]。每当蚂蚁选择下一个客户时，就随机产生一个在[0,1]内的数，根据随机数q的值按式(18)确定蚂蚁选择下一个客户。当q

每只蚂蚁考虑以下两个因素来确定选择下一个客户的概率：信息素及时间窗。蚂蚁本身根据信息素浓度的高低来选择下一个客户，同时考虑到达下一个客户的时间以及时间窗的大小，优先选择时间窗较小且到达时间较合适的客户。

(19)

其中，Jk表示t时刻蚂蚁k下一步允许服务的客户集(即还没有服务的客户)；ξ和ψ为权重系数，满足0≤ξ,ψ≤1，且ξ+ψ=1。α为信息启发式因子，反映了残留的信息量在蚁群的探索过程中的相对重要程度；β为期望启发式因子，反映了路径期望值在蚁群的探索过程中的相对重要程度；ηij表示由客户i转移到客户j的期望程度，其计算公式为：

(20)

2.2.2 自适应挥发因子

在实际情况中，挥发因子ρ的值往往是随机设定的，且挥发因子ρ的大小会影响算法的全局搜索能力以及收敛速度；若ρ过大，则路径上的信息素总量快速减少，会降低了算法的收敛速度，但利于寻找全局最优解；若ρ过小，则路径上的信息素总量增加，能提高算法的收敛速度，但易陷入局部最优解，全局搜索能力较差。为平衡算法的局部和全局搜索能力，避免陷入早熟，本文将挥发因子ρ设置为如式(21)的自适应调整方式：

(21)

2.2.3 信息素的设置

为避免由于某条路径上的信息素浓度过高而导致搜索陷入局部最优。将每条路径上的信息素τij(t)的界限设置为[τmin,τmax]，若τij(t)>τmax,则τij(t)=τmax;若τij(t)<τmin,则τij(t)=τmin。其中，τmax和τmin的计算公式为

(22)

(23)

其中，f为当前解中的最小成本，u为当前解中的最大满意度，σ为精英蚂蚁的个数。

2.2.4 信息素的更新

信息素的更新方法设置为两种：局部信息素更新和全局信息素更新。

局部信息素更新规则为：每只蚂蚁访问完一个客户后对各条路径上的信息素浓度进行更新，其计算公式为

(24)

其中，

(25)

全局信息素更新规则为：在每次迭代结束后，对精英蚂蚁经过的路径进行信息素更新。由于本文为双目标优化问题，因此在仿真时将精英蚂蚁分为两部分，一部分对应成本最低，另一部分对应满意度最高。通过调整两部分精英蚂蚁对信息素的影响来达到同时降低成本和提高客户满意度的目的。信息素计算公式为:

(26)

(27)

其中，将成本按升序排序，f为所有路径中的最小成本，u为所有路径中的最大满意度，fμ为排名第μ的路径花费的成本，该路径上的信息素增加量为△τij(t)。

2.2.5 局部搜索

变邻域搜索(VNS)由Hansen等[28]学者于1997年提出，基本思想是通过不同的邻域结构进行搜索，不断提高解的质量，具有较强的局部搜索能力。当所有蚂蚁搜索完之后，对当前最优解进行局部搜索。由于本文考虑的是双目标优化问题，因此当所有蚂蚁搜索完之后，找出当前解中成本最小及满意度最大的两条路径，分别对其进行变邻域搜索，扩大搜索空间。

结合所研究问题的特点，本文设计了3种邻域结构，分别为：1)2-opt算子，随机将路径中两点间的客户反序；2)Swap算子，随机交换路径中两个客户的位置；3)Insert算子，随机选取路径中两个客户，将后面位置的客户插入到前面客户的前面。

基于上述3种邻域结构的变邻域搜索的步骤如下：

1)确定邻域结构LYk,k∈{1,2,3}，设置最大循环次数itermax，令当前循环次数iter=1。

2)若达到最大循环次数，则输出最优解；否则令k=1。

3)对当前最优解进行局部搜索，若搜索后的新解比当前最优解更优，则更新当前最优解；否则k=k+1。

4)若k≤3则转入第3步；否则iter=iter+1，转入第2步。

需要说明的是，子路径间的邻域结构有可能会破坏解的可行性,即可能会违反车辆容量限制，因此,在进行变邻域搜索前，去掉路径中的配送中心，将路径转换为只有客户的路径，再进行变邻域搜索，搜索完的路径再按照车辆容量来划分子路径，分配车辆。

2.2.6 算法步骤

综上所述，利用HACO算法求解带时间窗的车辆路径问题的流程如图1所示。

图1 HACO算法流程图

3 仿真实验与结果分析

为了验证HACO算法求解带时间窗的车辆路径问题上的有效性，利用Matlab R2015b软件进行两组仿真实验，操作系统为Windows 10，处理器为Intel Core i7-8550U，主频为1.8GHz，内存为8GB。

3.1 参数设定

表1 正交组合及测试结果表

表2 参数水平取值表

表3 参数平均响应值表

由表3可知，参数ξ对算法的影响最为显著，其次是参数α和β，由图2可知，算法中合适的参数组合为α=4，β=1，ξ=0.9。

图2 各参数对算法性能的影响趋势

3.2 评价标准

为衡量算法的性能，选用以下常用的性能标准：

1)非支配解的数量(NPS)：算法所得非支配解越多，说明算法能保持较好的解的多样性。

2)解集间的覆盖率(C)[29]：C(E,F)表示集合F中至少被集合E中的一个解支配的解所占比例，其计算公式如式(28)：

(28)

C(E,F)的取值范围为[0,1]，值越大，则说明F中有更多的解被E中的解支配，而E中解的质量比F更好，需要注意的是：C(E,F)与C(F,E)需要分别计算，C(E,F)≠1-C(F,E)。

3.3 仿真实验1

仿真实验1以文献[20]中的算例进行仿真实验，在算例中增加了宽松时间窗。算例参数设置如下：由一个配送中心为20个客户服务，配送中心位于各客户的中间位置，坐标为[70,70]，客户分布为辐射型。配送车辆最大载重为8 t，所有配送车辆规格一致，平均车速为40 km/h，车辆使用成本为100元，行驶成本为2元/km。等待成本为10元/h，惩罚成本为40元/h。由于如今加油站随处可见，将配送车辆的最大行驶路程设为无限制。客户信息如表4所示(其中编号0为配送中心，编号1～20为需要服务的客户)。

表4 客户信息表

将本文所讲算法与传统蚁群算法及NSGA-Ⅱ算法进行对比。种群大小等于客户数，最大迭代次数NCmax=500。每个算法独立运行40次，结果如表5所示，图4为HACO算法求解的最优路径。

图4 最优路径

表5 与改进前算法比较结果

表5为改进的HACO算法与改进前的ACO算法、NSGA-Ⅱ算法比较结果，从表中可知：HACO的平均值、最优值、最差值均优于ACO算法及NSGA-Ⅱ算法。HACO算法获得的非支配解的数量最多，NSGA-Ⅱ算法获得的非支配解的数量高于ACO算法，但是解的质量低于ACO算法及HACO算法。这主要是由于HACO算法相比于NSGA-Ⅱ算法增加了对当前最优解进行局部探索，HACO算法相比于ACO算法不仅增强了其局部探索能力、再加上挥发因子ρ在陷入局部最优解的时候自动衰减，有助于跳出局部最优，能找到更多高质量的解，且精英蚂蚁策略针对两个目标分别进行搜索的方式可以保证解的多样性。表6为3种算法两两配对计算得到的覆盖率，表格对应C(E,F)，其中E为行对应的算法，F为列对应的算法。如C(HACO,ACO)=0.826 1,说明ACO获得的非支配解中有82.61%的解被HACO的非支配解所支配；C(HACO,NSGA-Ⅱ)=0.866 7,说明NSGA-Ⅱ的非支配解中有86.67%的解被HACO的非支配解所支配，而C(ACO,NSGA-Ⅱ)=0.666 7，说明NSGA-Ⅱ的非支配解中有66.67%的解被ACO的非支配解所支配。表6的实验结果也证实了本文所提的HACO算法的有效性。从图3可以直观地看出HACO的性能显著优于初始的ACO算法及NSGA-Ⅱ算法，其Pareto解占据明显的支配地位。ACO及NSGA-Ⅱ的解呈交叉式，即成本较高的NSGA-Ⅱ的非支配解优于ACO的非支配解，成本较低的则呈现相反的状态。实验结果显示：三种算法均没有找到满意度为1的路径，这并不代表不存在满意度为1的配送方案，而是因为本文研究的为双目标车辆路径问题，算法在运行过程中会同时兼顾两个目标探索方向，不存在不顾成本地提高满意度以及单纯降低成本而不顾客户满意度的情况。

表6 与改进前算法的覆盖率(C)比较结果

图3 Pareto前沿

表7为本文算法与文献[30]的文化狼群算法、文献[31]的改进离散花朵授粉算法、文献[32]的单点单亲遗传混合蚁群算法以及文献[20]所提算法中性能最优的单点重组精英遗传混合蝙蝠算法就仿真实验1的结果进行对比。从最少车辆数来说，只有单点单亲遗传混合蚁群算法的最少车辆数为7，其余的均为6。从最短路径来说，HACO算法的最短距离为999.26 km，比文化狼群算法减少了0.64%，比改进离散花朵授粉算法减少了5.27%，比单点单亲遗传混合蚁群算法减少了13.66%，比单点重组精英遗传混合蝙蝠算法减少了8.34%。HACO算法求得的平均最短路径长度也优于其余4种算法，说明HACO具有更好的稳定性。在5种算法的比较中，HACO算法求得的解质量更好，而且更稳定，说明HACO算法能有效减少物流配送成本，提高服务质量。

表7 与其他算法比较结果

灵敏度分析。带时间窗的车辆配送受车辆行驶速度的影响。从平均20 km/h的车速开始，以5 km/h的速度增加，不断改变车辆平均行驶速度并分别运行算法20次，记录最小总成本及最大客户满意度。比较不同车速下的总成本及满意度最优值的变化如图5所示。

由图5可知，车速的提高在一定程度上能降低物流企业的成本并保持较高的客户满意度，但当车速增长过大时，会导致成本增加并降低满意度，其原因是由于车速的增加，车辆在路上行驶的时间较短，过早到达客户处，导致在客户处的等待时间增加，增加了车辆的等待成本并降低了满意度，且过快的行驶速度易导致交通事故，造成重大的财产和人员损失。因此，控制车辆合适的行驶速度对物流企业至关重要。

图5 平均车速对最优解的影响

3.4 仿真实验2

为进一步验证本文所提HACO算法求解不同类型和规模问题的性能，本文选用Solomon标准测试集中的R101、C101、RC101三个算例作为测试问题，其中R101中的客户坐标为随机分布，C101中的客户坐标为聚类分布，RC101中的客户坐标为混合随机聚类分布，每个算例的客户数量为100。为满足本文时间窗设定，对算例中的时间窗做出调整。将HACO算法与ACO算法、NSGA-Ⅱ算法进行性能比较，不同算法的比较结果如表8所示，图6为C101算例的优化结果。

表8 算例对比结果

从表7可以看出，在最小总成本、最大满意度方面，HACO明显优于初始的ACO与NSGA-Ⅱ算法。HACO算法在ACO中加入了局部搜索算法，综合提高了HACO算法的局部搜索和全局搜索能力，并通过NSGA-Ⅱ算法更快地选择出更好的解。总体来说，本文提出的HACO算法能稳定地获得更优的解，能有效解决带时间窗的车辆路径问题，求解性能明显提高。

4 结论

在带时间窗的车辆路径问题中，影响配送成本的因素主要是客户的位置及时间窗，客户位置不可变动，客户时间窗分布较为复杂但可变。当客户根据自身需求提出服务时间窗时，物流企业可根据路径规划对客户反馈，对能及时服务的客户协调好服务时间，对不能及时服务的客户调整至下一时间段优先配送。从而达到物流企业与客户的双赢。本文对更大客户规模的算例实验不足,今后可以通过改进算法、提高程序的运行效率、提升硬件环境来解决更复杂的大规模求解问题。