基于客户满意度的车辆路径问题的混合蝙蝠算法

孙 奇，马 良

（上海理工大学 管理学院，上海 200093）

车辆路径问题（vehicle routing problem, VRP）作为运筹学中一个典型的组合优化问题，属于NP难题。该问题起源于车辆运输领域，最早由Dantzig等[1]于1959年首次提出并研究。由于车辆路径问题有着良好的实际应用价值，该问题自提出以来，便引起了管理学、运筹学、图论、组合数学、计算机应用等领域学者的关注，研究者们对其进行了大量的研究，并取得了很多成果。

带时间窗的车辆路径问题[2]（vehicle routing problem with time windows, VRPTW）最 早 由Savelsbergh提出，是在VRP的基础上考虑了客户对接受服务时间的需求。相比于传统的VRP，VRPTW更贴近实际情况。VRPTW也属于NP难题，当问题规模不断增大，其计算时间呈现爆炸式增长，精确算法难以甚至无法求出最优解。智能优化算法逐渐被大量使用，如遗传算法[3-4]、蚁群算法[5-8]、模拟退火算法[9]、粒子群算法[10-12]等,并取得了一定的成果。

蝙蝠算法（bat algorithm, BA）是由剑桥大学Yang[13]于2010年最早提出。随着对该算法的进一步研究，国内外学者提出了众多的改进蝙蝠算法[14-17]用于解决各领域优化难题。同大多数智能优化算法一样，蝙蝠算法也具有易陷入局部最优、后期收敛速度慢、局部搜索能力弱等缺点。针对蝙蝠算法的这些劣势，本文利用GRASP（greedy randomized adaptive search procedure）启发式算法为蝙蝠算法构建初始解，以提高蝙蝠算法的寻优精度。同时，在蝙蝠算法的主体中引入病毒进化（virus evolution, VE）机制，使病毒群体和蝙蝠群体以一种协同进化的方式进行演化，以增强蝙蝠算法跳出局部最优的能力。最后，通过实例验证了该算法的有效性。

1 基于客户满意度的车辆路径问题

1.1 问题描述

基于客户满意度的车辆路径问题是多目标VRPTW，该问题不仅要满足配送车辆的行驶里程最短，还需要考虑客户满意度，使得总客户满意度最大。为此，该问题需要满足以下假设条件：a. 所有车辆从配送中心出发，当前路径配送完毕后，所有客户必须返回配送中心；b. 各条路径上的客户需求总量不能超过车辆的最大装载量；c. 所有客户点都需要被服务，一辆车可以服务多个客户点，但一个客户只能由一辆车提供服务；d. 车辆到达客户的时间可以提前，车辆提前到达客户点需要等待至可接受服务的时间节点，但不能迟于客户可接受的最晚服务时间。

1.2 模型建立

1.2.1 引入客户满意度

现有客户i可接受服务的时间窗为 [E Ti,LTi]，现假设车辆到达客户i的时间为ti，在硬时间窗车辆路径问题的要求下，超出时间窗以外的客户将不会接受服务。而在实际物流配送过程中，虽然货物不在时间窗内到达的情况会引起客户不满，但客户对于这种情况都有一定的容忍度，只不过如果在时间窗外，客户的满意度会降低。本文在原有时间窗的基础上构建模糊时间窗，设客户能接受服务的最早时间为 E，能接受服务的最晚时间为 L T。货物如果在 [E Ti,LTi]内到达，客户的满意度最高为100，在时间窗之外，但在[ETL内到达，虽然客户满意度呈现下降趋势，但客户仍然会接受服务。构建的模糊时间窗的客户满意度如图1所示，对应的客户满意度函数如式（1）所示。

客户满意度函数：

1.2.2 数学模型

目标函数：

约束条件：

式中：M=[0，1，2···N]表示所有节点集合，其中0表示配送中心； 为客户点；

为配送车辆集合；Q表示车辆的最大载重量；di表示客户点i的需求量；ti为车辆到达客户i的时间；cij为客户i到客户j的距离；wi为车辆在客户i的等待服务时间；S为某辆车的配送客户集，|S|为集合S中所包含顶点的个数； 和

为决策变量。M={0,1,2,···,n}N={1,2,3,···,n}K={1,2,3,···,m}

xijk yik

式（2）和式（3）分别表示要求的最大化客户满意度和最小化总车辆行驶里程的目标函数；式（4）为车辆载重限制；式（5）、式（6）和式（7）保证每个客户点被且仅被服务一次，每个客户的需求只由一辆车来完成；式（8）为消除子回路；式（9）为时间窗约束。本文采用线性加权法将两个目标转化为单目标优化，即

2 求解基于客户满意度的车辆路径问题的混合蝙蝠算法设计

式中， 表示目标权重，ω∈(0,1)。

2.1 基本蝙蝠算法

蝙蝠算法是受到自然界蝙蝠回声定位的启发演变而来的新型智能优化算法，蝙蝠通过改变声波的频率、声音的响度、脉冲发生率，从而锁定目标进行捕食。

a. 蝙蝠速度位置更新。

基本蝙蝠算法是通过频率f实现对蝙蝠移动步长和范围的控制，从而达到更新蝙蝠的速度和位置。设m只蝙蝠于时刻t散布在D维搜索空间内，则第i只蝙蝠在t+1时刻位置和速度更新公式如下：

式中：fi表示第i只蝙蝠发出的脉冲频率；fmax和fmin分别表示蝙蝠发出脉冲频率的最大值和最小值；是一个随机变量，其取值范围为 ；v表示蝙蝠速度；X*表示当前蝙蝠种群中的最优位置。

同时基本蝙蝠算法内部可设定局部搜索行为，产生一个随机数rand1，通过比较rand1和脉冲发射率 ri，若 rand1＞ri，根据式（14）在当前最优解集中选取一个最优蝙蝠，进行随机扰动，从而产生一个新位置。β

式中：Xnew表示更新后的蝙蝠位置；Xold表示更新前的蝙蝠位置； 为一个随机数；At表示当前时刻下所有蝙蝠的平均响度。

b. 蝙蝠响度和脉冲发生率更新。

蝙蝠在找到猎物时，会降低自身的声波响度并且增加脉冲发生率。响度和脉冲发生率按照以下公式进行更新。

式中： 表示第i只蝙蝠在t时刻的脉冲发生率；

表示响度衰减系数，表示脉冲发生率增加系数。 与 均为常量，通过 与 的调节，从rti α∈(0,1)γ＞0

αγαγ而达到平衡整个算法的全局搜索和局部搜索。

2.2 混合蝙蝠算法

2.2.1 基本蝙蝠算法速度的重新定义与编码

基本蝙蝠算法最初被用来解决的是连续函数优化问题，在面对离散的组合优化问题时，基本蝙蝠算法很难直接应用，需要根据实际问题对蝙蝠位置的编码方式以及速度进行重新定义。本文借鉴文献[18], 利用汉明距离（HD）来对速度进行定义，并采取二维编码方式。现以一个有7个客户点的车辆路径问题为例，蝙蝠位置编码如下：

上述编码表示使用3辆车，其中，车辆1路径：0—1—3—0；车辆2路径：0—5—2—4—0；车辆3路径：0—7—6—0。

蝙蝠速度定义位置更新操作如式（17）和（18）所示：

汉明距离表示蝙蝠两个位置之间的差距，举例如下：

则在此例中，蝙蝠当前位置和最优位置之间的汉明距离为4。在本文中具体操作为：针对蝙蝠当前位置随机选取两对对位点进行2-opt操作则在蝙蝠当前位置中随机选取一对对位点进行2-opt操作。

2.2.2 GRASP启发式算法构造初始解

GRASP启发式算法[19]已经被证明是一种求解组合优化问题的有效方法，它由两个阶段构成，分别是构造阶段和局部搜索阶段。在构造阶段中通过带有随机特性的贪心算法进行初始解的选择，该算法具有贪心和随机两种特性。而在局部搜索阶段中，通过局部搜索策略对构造阶段所得到的解进行优化，以提高解的质量。本文局部搜索阶段所采取的局部搜索策略为2-opt方法，产生数量为0.8N（N为种群总数量）的初始解，剩余的种群初始解则由随机策略产生。

2.2.3 引入病毒进化机制

病毒是自然界一种神奇的生物，它具有特殊的感染能力和进化能力。在病毒进化算法中，病毒一方面通过感染能力，将自身的遗传信息传递给其他个体；另一方面通过自身的进化能力，从宿主个体中获取部分遗传信息来改变自身的遗传信息，从而达到进化的目的。病毒进化算法的具体设计如下所示。

a. 病毒个体的编码方式：病毒个体（virus）是从宿主个体（host）中抽取一定片段，长度与宿主个体保持一致均为n，在病毒个体中加入通配符“*”。

图2 病毒个体编码Fig.2 Virus individual coding

图3 病毒感染操作Fig.3 Virus infect operation

c. 病毒个体适应度：一个病毒个体是否优劣是由它的适应度函数值来表现，适应度值（fit）用病毒在感染宿主个体后，宿主个体适应度值的变化来表达。假设宿主个体集合为N，则第i个病毒的适应度计算为

d. 病毒的复制操作：病毒通过复制操作将宿主个体的基因以概率复制到自身基因中，以产生拥有更多有效字符串的病毒个体。在初始化病毒种群的时候随机选取宿主个体经复制操作产生病毒个体的初始种群，具体操作如图4所示。

图4 病毒复制操作Fig.4 Virus copy operation

图5 病毒剪切操作Fig.5 Virus cut operation

f. 病毒的生命力：每一个病毒都有自己的生命力（Life）, 当生命力小于0时，意味着病毒已经死亡，此时需要产生新的病毒来维持病毒种群，病毒生命力定义为

2.3 混合蝙蝠算法设计

本文针对基本蝙蝠算法的局限性以及车辆路径问题的具体情况，引入GRASP启发式算法和病毒进化机制，对基本蝙蝠算法中速度以及位置更新方式进行重新定义，提出了GRASP-VREBA算法，该算法具体步骤如下：

Step 1 初始化GRASP各参数，定义贪婪函数，通过构造阶段和局部搜索两个阶段得到数量为0.8 N的初始解，剩余初始解则采取随机生成的方法产生；

Step 2 初始化蝙蝠算法和病毒进化相应的参数，包括蝙蝠速度vi、响度Ai、脉冲频率ri、 响度衰减系数、脉冲频率增加系数，以及病毒感染、复制和剪切概率病毒生命力衰减系数等，并通过病毒复制操作产生病毒的初始种群；

Step 3 根据式（10）计算各蝙蝠的适应度值，排列蝙蝠，确定最优蝙蝠的位置；

Step 4 进行蝙蝠主种群的进化操作，根据重新定义的速度和位置更新公式，进行蝙蝠位置和速度的更新，并产生新解；

Step 5 将每一个病毒对所有蝙蝠进行感染操作，并计算蝙蝠的适应度值和生命力，若生命力小于0，则随机选择一个宿主个体进行病毒复制操作，否则对病毒个体进行剪切操作，当病毒生命力小于0时，通过病毒复制产生新的病毒；

Step 6 判断是否达到最大迭代次数或者已经达到最优解，否则返回Step 3继续运行。

3 仿真实验与分析

3.1 算例描述与数据

为了验证算法的有效性，本文对文献[8]算例的时间窗进行扩展，扩展方式为客户可接受服务的最早时间扩展为在原时间窗的基础上提前20%，客户可接受服务的最晚时间扩展为在原时间窗的基础上延长10%。

算例中，各参数如下：车辆最大载重量为8 t，为各客户点提供货物配送服务，各客户信息如表1所示；车辆行驶速度为40 km/h, 如果车辆提前到达必须等待至可接受服务时间节点后才可进行服务。在所有客户被服务完毕后，所选路径需要达到的目标为：车辆总行驶里程最少，客户满意度最大。本文权重取=0.2。

表1 配送中心与客户信息Tab.1 Distribution center and customer data

3.2 结果分析

本文应用Matlab R2014a对基本蝙蝠算法、GRASP-BA算法、GRASP-VEBA算法进行编程实现，这3种算法均在Intel Core i7-7700 3.6 GHz（8.00 GB RAM）、操作系统为Win10的环境下运行。3种算法的参数设置如下：=0.95，=0.95，初始响度；初始脉冲频率；蝙蝠种群数n=100，最大迭代次数为600。在基本蝙蝠算法中，蝙蝠种群为300。GRASP-VEBA算法中，病毒种群数目=0.3n；= 0.2；= 0.2，= 0.1。 GRASP-BA和 GRASPVEBA算法中，阀值=0.3。

3种算法均随机运行15次，优化结果如表2所示。基本蝙蝠算法求出的最优路径中，车辆总行驶里程为1 207.2 km，最大客户满意度为1 524.3；GRASP-BA算法求出的最优路径中，车辆总行驶里程为1 148.8 km，最大客户满意度为1 524.3；GRASP-VEBA算法求出的最优路径中，车辆总行驶里程为1 064.8 km，最大客户满意度为1 524.3。可见改进的蝙蝠算法均优于基本蝙蝠算法。尤其在GRASP-VEBA算法中，相比于BA算法和GRASPBA算法平均长度均减少超过10%，分别为16.9%和12.2%，平均客户满意度提高分别为3.1%和1.9%，验证了算法的有效性和稳定性。

表2 3种算法计算结果比较Tab.2 Comparison of the calculation results of three algorithms

最优路径图如图6所示，最优路径各车配送路线为

车辆 1：0—18—5—11—14—0

车辆 2：0—15—4—9—0

车辆 3：0—12—19—13—6—0

车辆 4：0—3—7—0

车辆 5：0—20—10—17—2—0

车辆 6：0—1—16—8—0

图6 最优方案路径图Fig.6 Optimal scheme routes

4 结 论

本文在传统单目标车辆路径问题的基础上增加客户满意度最大化的目标，使车辆路径问题更贴合生活实际。建立基于客户满意度的车辆路径问题模型，针对该问题，在基本蝙蝠算法的基础上引入GRASP启发式算法构造初始解和病毒进化机制进行协同进化，设计出求解多目标车辆路径问题的混合蝙蝠算法, 并用其求解测试算例。求解结果表明算法的有效性，能够较好地弥补基本蝙蝠算法后期收敛速度慢、缺乏有效局部搜索机制等的不足。