超启发式分布估计算法求解带软时间窗的同时取送货车辆路径问题

张烜荧,胡 蓉,钱 斌

(昆明理工大学信息工程与自动化学院自动化系,云南昆明 650500;昆明理工大学云南省人工智能重点实验室,云南昆明 650500)

1 引言

随着“制造强国战略”的实施,在推动先进制造业迅速发展的同时,现代物流业的高质量发展也成为大势所趋.在物流系统中,配送是核心部分,而车辆路径问题(vehicle routing problem,VRP)是其中的重要一环.传统的VRP主要描述为组织一个车队,安排适当的行车路线为一定数量的客户提供送货服务,并能在满足车载量、行车里程等约束条件下,达到诸如车辆行驶总里程最短,成本最小,耗时最少等目的.随着资源需求的增加与供给双向流通方式的快速发展,以及物流企业对降低成本、保护环境的注重,考虑逆向物流的车辆路径问题开始受到重视.此外,日益激烈的市场竞争逼迫企业在满足客户需求的同时要考虑提高客户满意度.在此背景下,研究带软时间窗的同时取送货车辆路径问题(vehicle routing problem with simultaneous pickup and delivery and soft time windows,VRPSPDSTW)具有重要意义.由于带时间窗的同时取送货车辆路径问题(vehicle routing problem with simultaneous pickup and delivery and time windows,VRPSPDTW)为NP-hard问题[1],而VRPSPDTW可归约为VRPSPDSTW,故VRPSPDSTW也属于NPhard问题,对其展开研究具有较大理论价值.

VRPSPDSTW的求解算法主要分为两类.一类是运筹学算法,包括分支定价、动态规划和拉格朗日松弛等[2–8],基本都是用于求解单目标问题.此类算法以线性代数和几何分析为基本工具,利用问题优化目标函数和约束式的结构信息构造搜索,可在几分钟至几十分钟内获取较小规模问题(客户数小于等于20)的最优解.但是,VRPSPDSTW具有非凸特性,其几何结构和最优解间的关系仍属开放问题,不存在多项式时间获取最优解的算法,故对应的运筹学算法均是靠遍历或部分遍历解空间来确保求解质量,这使其需要数小时至数天才可获得较大规模问题(客户数在80～200之间)的最优解或满意解.另一类是智能算法,包括遗传算法、禁忌搜索算法、蚁群算法等[9–29],既用于求解单目标问题,也用于求解多目标问题.此类算法不依赖于问题结构,而是将问题的约束条件隐式地包含在所设计的解的编码、解码规则中来处理,并利用某种拟人、拟物的机制不断生成新的可行个体或解,从而引导算法执行搜索,往往可在几秒或十几秒内就能获得不同规模问题的满意解.因此,设计智能算法以实现对VRPSPDSTW的快速有效求解,是合理且必要的.

近年来,智能算法求解VRPSPDSTW逐渐受到重视.针对单目标VRPSPDTW,Lai等[9]提出改进差分进化算法(improved differential evolution algorithm,IDE)求解问题,取得较好的效果.Wang等[1]提出一种基于最小代价插入法的协同进化遗传算法(genetic algorithm,GA)进行求解,实验结果表明所提算法的有效性.王等[10]提出模拟退火算法(simulated annealing,SA),并在模拟退火过程中融合局部搜索方法,在可接受的计算时间内取得了较好的优化结果.王等[11]提出一种离散布谷鸟算法(discrete cuckoo search,DCS),其L´evy飞行机制协调局部搜索和全局搜索,实验结果更新了国际最优解.韩等[12]针对带折线型软时间窗的车辆路径问题(vehicle routing problem with soft time windows,VRPSTW),以最小化运输总成本为目标,提出一种超启发式遗传算法(hyper-heuristic genetic algorithm,HHGA),以GA 作为上层搜索算法,CW 节约法、MJ插入法、Kilby插入法作为底层搜索规则,并通过加入排序规则和解内部调整算法以进一步提高解的质量.实验结果验证了HHGA 的可行性和有效性.王等[13]提出一种回溯搜索优化算法(backtracking search optimization algorithm,BSA),使用启发式规则产生初始种群,通过多种搜索算子更新局部最优解,实验验证了BSA的有效性.Qin等[14]进一步考虑了低碳需求,提出一种自适应遗传爬山算法(adaptive genetic hill-climbing algorithm,AGHCA)进行求解,同时研究了碳税对总成本和碳排放的影响.

对于多目标VRPSPDSTW(multi-objective VRPSPDSTW,MOVRPSPDSTW)及相关问题,李等[15]针对VRPSPDTW,以最小化总车辆数、运输总时间、车辆行驶总里程的线性组合为优化目标,提出一种基于路线集合划分的分解迭代算法,将路线集合分解为路线子集合,再用记录更新法求解子集合,不断迭代逐渐逼近最优解,实验结果表明了算法的有效性.Banos等[16]针对考虑负载均衡的带硬时间窗的车辆路径问题(vehicle routing problem with time windows,VRPTW),以同时最小化各车辆间最大行驶距离之差和车辆间最大负载之差为优化目标,提出一种结合模拟退火算法的多起点多目标进化算法(multi-start multiobjective evolutionary algorithm with simulated annealing,MMOEASA).仿真实验和对比结果验证了算法的有效性.Yang等[17]针对带软时间窗的车辆路径问题(VRPSTW),以最小化总运营成本、时间窗偏离总量、碳排放为优化目标,提出一种混合遗传算法(hybrid genetic algorithm,HGA),通过实验验证了算法的有效性,并通过帕累托分析,说明了目标之间的关系.Sumaiya等[18]针对VRPSTW,以最小化总行驶距离、违反时间窗的服务次数、车辆数为优化目标,提出一种人工蜂群算法(artificial bee colony,ABC),结合车辆内和车辆间的两步局部搜索进行求解,实验结果表明ABC能在合理的计算时间内获得高质量的解.Zhang等[19]针对时间依赖型同时取送货车辆路径问题(time dependent vehicle routing problem with simultaneous delivery and pickup,TDVRPSDP),以最小化碳排放和旅行时间为优化目标,提出一种多目标协同进化量子算法(multi-objective cooperative quantum evolutionary algorithm,MCQEA),将种群分为两个亚种群,使用非劣解集和亚种群中的解共同产生新种群,实验结果验证了MCQEA的有效性.Li等[20]针对MOVRPSPDSTW,以最小化车辆数、车辆固定成本、行驶总距离、旅行时间、服务成本为优化目标,使用基于分解的多目标化学反应优化算法(chemical reaction optimization,CRO),将复杂的MOVRPSPDSTW转换为多个单目标子问题,并同时在CRO中对其进行优化,实验结果验证了算法的有效性.陈等[21]针对多目标同时取送货车辆路径问题(multi-objective vehicle routing problem with simultaneous pickup and delivery,MOVRPSPD),以最小化总路径成本、各车辆最大行驶里程差为优化目标,提出一种嵌入禁忌表、且具有贪婪转移准则的多目标改进蚁群算法(improved ant colony optimization,IACO),实验结果表明IACO可求得权衡各优化目标且使单一优化目标近似最优的Pareto解.柴等[22]考虑路径具有双重属性,建立基于双属性的多目标VRPTW,进而提出一种蚁群算法(ant colony algorithm,ACO),在算法中定义了一种综合考虑各优化目标、时间窗和信息素等启发信息的状态转移概率公式.实验结果表明,所提算法在运行时间、收敛性和群体多样性方面具有优势.Wang等[23]针对带时间窗的周期性车辆路径问题(periodic vehicle routing problem with time windows,PVRPTW),以同时最小化车辆行驶总距离、最小化车辆总数、最大化客户访问频率为优化目标,提出一种多目标模拟退火蚁群优化算法(multi-objective simululate annealing-ant colony optimization,MOSAACO),在算法中设计了基于欧氏距离的优质解接受方式和可生成不同初始解的参数控制策略,并加入了有效局部搜索策略.实验表明,该算法在求解PVRPTW问题上性能良好.由文献调研可知,虽然VRPSPDSTW已得到一定研究,但尚无人考虑客户满意度(服务准时率)这一重要优化目标.此外,上述算法在全局搜索时基本都直接采用算法原有搜索机制执行,未考虑结合问题特性做合理更改来增强算法发现优质解的能力,且局部搜索时大多迭代执行较为单一的插入、变换操作,难以实现对复杂问题解空间中优质解区域的深入搜索.显然,针对VRPSPDSTW这类复杂问题,如何设计能获取高质量解的有效算法是需要进一步研究的课题.

分布估计算法(estimation of distribution algorithm,EDA)是一种新兴的基于统计学原理的随机优化算法.EDA采用统计学习手段从群体宏观角度建立概率模型,用来学习和积累优质解在搜索空间的分布信息,然后对此概率模型采样产生新种群,如此反复实现种群进化[24].此进化方式可一定程度上避免传统进化算法中交叉、变异等操作所带来的模式破坏问题,从而使其搜索具有更好的引导性[25].目前,使用EDA求解VRP的研究较为有限.Wan等[26]针对VRPSTW,提出一种混合粒子群算法的混合EDA,通过粒子群算法的快速更新策略修改了EDA的概率矩阵,实验结果验证了算法的有效性.孟等[27]针对带硬时间窗的车辆路径问题(VRPTW),提出一种混合种群增量学习算法(hybrid population-based incremental learning algorithm,HPBIL),为每辆车建立2维EDA概率模型,同时设计了基于插入和两点邻域交换的两阶段局部搜索来增强算法的局部开发能力,实验结果表明混合后的算法求解效果优良.Wang等[28]针对客户需求和运输成本不确定的车辆路径问题(uncertain capacitated vehicle routing problems,UCVRP),提出一种改进EDA,不仅利用概率模型学习优质解信息并引导全局搜索,而且引入插入操作和2-opt操作来提高局部搜索能力,同时采用多场景技术,增强算法求解不确定问题的鲁棒性.Ji等[29]针对考虑碳排放的多目标联合运输VRP,提出一种混合EDA(hybrid estimation of distribution algorithm,HEDA),采用扩展的2维概率矩阵学习各阶段运输模式对应的优质非支配解信息,同时设计了基于部分交换操作的局部搜索,对问题进行了有效求解.从EDA在VRP相关问题上的研究现状来看,采用基于问题的概率模型可以合理引导全局搜索的方向,引入有效的局部搜索策略能进一步提高算法的性能,两者的协同搜索是算法设计的关键.据文献调研可知,针对VRPSPDSTW,尚无EDA求解的报道,故十分必要开展相关研究.

综上,本文研究多目标VRPSPDSTW的建模与求解.建模方面,建立同时以最小化总行驶路径和最大化服务准时率为优化目标的VRPSPDSTW.求解方面,在对问题模型特点分析的基础上,提出一种超启发式分布估计算法(hyper-heuristic estimation of distribution algorithm,HHEDA).HHEDA在全局搜索阶段执行分布估计算法EDA,在局部搜索阶段执行学习型超启发式局部搜索(learning hyper-heuristic local search,LHHLS).其中,EDA利用结合问题特性所设计的3类概率矩阵协同学习和积累优质解(即可更新当前非劣解集的解)信息并引导搜索尽快到达优质解区域;LHHLS采用学习型评分机制和11种有效邻域操作构造一系列启发式算法或搜索,然后分别对当前种群的部分个体和当前非劣解集中的解逐一执行相应的启发式搜索,可进一步加大搜索深度,从而增大算法在较短时间内获取高质量解的机率.最后,通过对不同测试问题的仿真实验和算法比较,验证了HHEDA的有效性.

2 问题描述

关于本文所涉及的数学符号及定义如表1所示,VRPSPDSTW示例如图1所示.定义一个无向图G(V,E),其中,V为节点集且VVc ∪{0},E为路径集.车辆从配送中心出发对若干客户进行服务,在约束条件下同时满足客户取送货需求.要求合理安排路径,最小化车辆行驶总里程和最大化服务准时率.

表1 符号表Table 1 Symbol table

图1 VRPSPDSTW问题示例Fig.1 VRPSPDSTW problem example

根据上述描述,建立如下VRPSPDSTW的问题模型:

其中:式(1)–(2)为优化目标,式(1)表示最小化车辆行驶总里程,式(2)表示最大化服务准时率;式(3)–(12)为约束条件,式(3)表示每个客户只被一辆车服务一次,式(4)表示客户点车辆流平衡约束,即车辆到达客户点h后必须离开该客户点;式(5)–(6)表示每辆车起始于配送中心并且完成服务后最终回到配送中心;式(7)–(8)表示车辆从配送中心出发时的载重,且该载重不能超过车辆最大载重;式(9)–(10)表示车辆服务完客户j的载重,且该载重不允许超过车辆最大载重;式(11)表示车辆从i出发服务j,到达j的时间.若车辆仅服务客户i,且满足式(12),那么这个客户可以被准时服务.

由上述问题模型可知,相对于传统VRP,本文的VRPSPDSTW进一步考虑软时间窗和同时取送货,使得模型约束增多,问题的解空间更加复杂;而且,增加了优化目标(即最大化服务准时率),使得解空间中优质解(即非劣解)的分布更加分散;此外,随着客户数量的增加,VRPSPDSTW的解空间规模呈指数增长[11].这些表明求解VRPSPDSTW的难度很大.

3 HHEDA

HHEDA的框架如图2所示.从图2可知,在HHEDA的初始化阶段,为提高初始种群的质量并加快算法收敛速度,提出3条规则用于生成初始解.在HHEDA 的全局搜索阶段,执行分布估计算法EDA.在EDA中,结合问题特性,设计排序概率矩阵、距离概率矩阵、捆绑关系概率矩阵共同学习和积累优质解(即可更新当前非劣解集的解)信息,然后通过采样这些概率矩阵以生成较高质量的新个体或解,并更新非劣解集和概率矩阵,可较快引导算法到达解空间中存在优质解的区域.在HHEDA的局部搜索阶段,执行学习型超启发式局部搜索LHHLS,以增强算法的局部寻优能力.在LHHLS中,低层问题域采用11种有效邻域操作组成备选集合,高层控制域设计改进的学习型评分机制动态评价底层不同邻域操作在搜索过程中的贡献(即对两个优化目标的改进程度),进而依据各邻域操作的贡献率,对全局搜索阶段EDA种群中部分个体和非劣解集中的解,逐一执行由4种邻域操作(用轮盘赌方法确定)构成的启发式算法,以实现较深入的局部搜索并更新当前非劣解集和概率矩阵.

图2 HHEDA流程图Fig.2 HHEDA’s flow chart

本节后续对HHEDA各环节进行介绍.

3.1 解的表示

令解π的长度为L,π(π[1],π[2],···,π[L]),其中Nc+2 ≤L≤2×Nc+1.配送中心编号为0,客户编号为大于0的整数.车辆从配送中心开出前往客户处,服务一定数量的客户后返回配送中心.如图1所示,编码为π(0,1,2,3,4,0,5,6,7,0,8,9,10,0)的解表示车辆1从配送中心开出,依次服务客户1、客户2、客户3、客户4后返回配送中心.车辆2、车辆3的服务顺序释义以此类推.

3.2 全局搜索

3.2.1 种群初始化

为了保证种群的多样性和分散性,使用规则和随机方法产生规模为popsize的第一代种群.使用规则1和规则2产生2个解,剩下popsize-2个解使用规则3 产生.令未被服务的客户集合为C.具体规则如下:

规则1最近邻规则.

步骤如果车辆从配送中心出发,计算配送中心与C中每个客户的距离,选择最近的一个客户进行服务,如果有多个最近客户,随机选择其中一个进行服务,并将这个客户从C中去除;如果车辆从客户i出发,计算i与C中客户的距离,选择最近的一个客户进行服务,如果有多个最近客户,随机选择其中一个进行服务,如果车辆继续服务该客户不会违反约束,则将这个城市从C中去除,否则,车辆放弃服务该客户并返回配送中心.

规则2准时率最高规则.

步骤如果车辆从配送中心出发,计算C中每个客户的到达时间,选择在时间窗内到达的客户进行服务,如果有多个可在时间窗内到达的客户,选择其中最近的进行服务,如果没有可在时间窗内到达的客户,在C中选择距离最近的客户进行服务,将选中的客户从C中去除;如果车辆从客户i出发,计算车辆从i到C中每个客户的时间,选择在时间窗内到达的客户进行服务,如果有多个可在时间窗内到达的客户,选择其中最近的进行服务,如果车辆继续服务该客户不会违反约束,则将这个城市从C中去除,否则,车辆放弃服务该客户并返回配送中心,如果没有可在时间窗内到达的客户,车辆返回配送中心.

规则3随机规则.

步骤在C中随机选取客户进行服务,如果车辆服务该客户不会违反约束条件,那么将这个客户从C中去除,否则,车辆放弃服务该客户并返回配送中心.

3.2.2 概率矩阵更新

在设计概率矩阵时,通过研究问题的特性,引入问题的特定信息和知识,是提高算法搜索能力的有效途径.对于VRPSPDSTW,为了最小化车辆行驶总里程,车辆总是倾向于先服务那些离当前位置比较近的客户,所以在车辆行驶总里程较小的解中,存在一些由客户间距离决定的位置信息;为了最大化服务准时率,车辆总是倾向于按照时间窗先后顺序服务客户,若客户i,j的时间窗和距离满足式(13)的关系,那么当车辆连续服务客户i,j时,只要车辆在时间窗内到达i,那么j一定能在时间窗内被服务,即ij满足捆绑关系.所以,在服务准时率较高的解中,客户间由于时间窗和距离较为匹配,存在一些捆绑关系信息.在HHEDA中,用ρ,λ,τ三个矩阵分别学习优质解的排序信息、客户间的距离信息和捆绑关系.

概率矩阵初始化阶段,ρ和λ为Nc×Nc的全零矩阵,τ为Nc×Nc的矩阵,用式(14)对τ中元素值进行初始化.在概率矩阵更新时,选取非劣解集中新增加的个体集合Ω按式(16)–(20)更新概率矩阵.其中,在每张车辆服务顺序中,ρj′j表示客户j出现在位置j′的次数;Kij为决策变量,当客户i与j被连续服务时为1,否则为0;Ij′j为决策变量,当客户j出现在位置j′时为1,否则为0;γij为学习速率,其值为客户ij之间距离的倒数;Lij为决策变量,当客户i与j被连续服务并且车辆到达i,j的时间都在时间窗内时其值为1,否则为0.

3.2.3 概率矩阵遗忘操作

在迭代过程中,新产生的非劣解会支配掉部分原非劣解集中的解,为了使从新非劣解中学习到的信息在概率矩阵中占较大比重,在算法每一次迭代结束后对概率矩阵按式(21)–(23)执行遗忘操作,遗忘因子为α(α<1)

3.2.4 对概率模型采样产生种群

算法的第iteration(iteration ≥2)代,新种群由精英保留机制和采样概率矩阵产生的解共同构成.随机选择m个非劣解(mmin(ps,「popsize×15%」),ps为非劣解集中解的个数)作为精英解集,剩余popsize−m个解由概率矩阵采样生成.

采样过程中,根据式(24)–(26),考虑客户排序信息、距离信息和捆绑关系,产生较优个体.其中C为未被服务的客户集合,(i −1)′表示在车辆中被第i −1个被服务的客户,pij为客户j在车辆中被第i个服务的概率,每次采样均使用轮盘赌操作,在一定程度上保证了种群的多样性.

当产生位置i1的客户时,由于是车辆第1个服务的客户,所以仅考虑客户在解中的排序位置,由式(24)计算C中客户被选中的概率;当产生位置i(i≥2)的客户时,若客户(i −1)′被准时服务,则在计算C中客户被选中的概率时,着重考虑客户之间的捆绑关系,以最大化服务准时率为侧重优化目标,由客户排序信息和捆绑关系综合计算概率,计算公式如式(25)所示;若客户(i −1)′没有被准时服务,着重考虑客户之间的距离关系,以最小化总距离为侧重优化目标,由客户排序信息和距离关系综合计算概率,计算公式如式(26)所示:

综上,算法任意代数的种群初始化的伪代码如表2所示.

表2 种群初始化Table 2 Population initialization

3.3 局部搜索

对于多目标组合优化问题,全局最优非劣解集中的任何解(即全局最优非劣解)在所有邻域结构下均为局部最优非劣解,而接近全局最优非劣解的优质解往往在多种邻域结构下均为局部最优非劣解.采用单一邻域操作(对应一种邻域结构)的迭代搜索容易较早达到并陷入该邻域结构的局部最优非劣解,但该解的质量大都一般;而采用多种有效邻域操作(对应多种邻域结构)的迭代搜索可在到达多种邻域结构共同的局部最优非劣解前一直持续进行搜索,从而容易获取优质非劣解.此外,多目标组合优化问题的解空间十分复杂,大量优质非劣解分散在庞大的“峡谷型”解空间中靠近底部的多个区域[30].因此,为能对HHEDA全局搜索发现的优质解区域(即非劣解集)进行有效搜索,本文设计学习型超启发式局部搜索(LHHLS),用于执行丰富的多种变邻域搜索,以增强算法局部搜索的深度.

超启发算法为两层结构,底层问题域包含一组直接搜索问题解空间的启发式算法(low-level heuristics,LLHs),其中每个LLH可为启发式算法,也可以是规则或邻域操作;高层策略域包含确定底层不同LLH组合的某种算法.高层策略域的算法(high-level strategy,HLS)通过动态控制底层问题域的LLHs来不断组合生产新的启发式算法,进而使用这些新的启发式算法在底层搜索问题解空间.对于本文的LHHLS,底层问题域设计11种LLH,每种LLH为一类有效邻域操作;高层策略域的算法由改进的学习型评分机制和轮盘赌方法共同构成.前者用于定量地评估和统计每个LLH在搜索过程中的执行效果;后者依据所有LLHs的评分(即贡献)对其进行轮盘赌,以生成一系列新的启发式算法(即LLH的组合).

显然,搜索能力较强的LLH能以较大概率被选中,其他LLH也能以一定概率被选中,这有利于生成多种新的有效启发式算法.

下面对LHHLS各环节进行介绍.

3.3.1 LLH设计

LLH要根据具体问题进行设计.本文设计了2类共11种邻域变换操作作为低层启发式算子,LLH1至LLH5为车辆内路径局部搜索操作,LLH6–LLH11为车辆间路径局部搜索操作.

1) LLH1:车辆内客户顺序逆转操作.选择被同一张车连续服务的几个客户,进行逆转.

2) LLH2:车辆内不相邻客户交换操作.选择被同一张车服务的两个不连续的客户,进行交换.

3) LLH3:车辆内相邻客户交换操作.选择被同一张车服务的两个相邻的客户,进行交换.

4) LLH4:车辆内客户插入操作.选择一张车辆的服务路线ri,在ri中随机选择一个城市c,将其插入到ri中的其他位置.

5) LLH5:车辆内客户基于准时性的交换操作.选择一张车辆的服务路线ri,在ri中选择车辆晚于时间窗下限到达的城市c1和车辆早于时间窗上限到达的城市c2,交换c1和c2.

6) LLH6:车辆间客户插入操作.选中两张车辆的服务路线ri和rj,在ri中随机选择一个城市c,将其插入到rj中.

7) LLH7:车辆间客户交换操作.选中两张车辆的服务路线ri和rj,在ri中随机选择一个城市c1,在rj中随机选择一个城市c2,交换c1和c2.

8) LLH8:车辆间连续客户插入操作.选中两张车辆的服务路线ri和rj,在ri中随机选择两个被连续服务的城市c1和c2,将c1和c2插入到rj中.

9) LLH9:车辆间连续客户交换操作.选中两张车辆的服务路线ri和rj,在ri中随机选择两个被连续服务的城市c1和c2,在rj中随机选择两个被连续服务的城市c3和c4,交换c1,c2和c3,c4.

10) LLH10:车辆间客户交换操作.选中两张车辆的服务路线ri和rj,在ri中随机选择两个连续服务的城市c1和c2,在rj中随机选择一个城市c3,交换c1,c2和c3.

11) LLH11:车辆间客户基于准时性的交换操作.选择两张车辆服务路线ri和rj,在ri中选择车辆晚于时间窗下限到达的城市c1,在rj中选择车辆早于时间窗上限到达的城市c2,交换c1和c2.

在对解执行任意一种LLHi(1 ≤i≤11)操作后,若产生的新解为可行解且支配原解,则用新解替代原解,否则,不接受该新解.

3.3.2 基于改进程度的学习型评分机制

对于双目标优化问题,当前解进行LLHi(1 ≤i≤11)操作后得到的可行解存在3种可能:1)新解所有的目标都得到改进或一个目标得到改进,另一个目标保持不变;2)新解的一个目标得到改进,另一个目标变差了;3)新解的两个目标都变差了.每次使用LLHi对解π进行局部搜索后,采用所设计的评分机制(式(27)–(29))定量学习和积累LLHi的贡献率或得分信息,并依此采用轮盘赌方法构造新的启发式算法HLS(即局部搜索策略),用于对下一代的各解进行局部搜索.LLHi初始得分为Sstart,最高得分为Supper,最低得分为Slower.l(π)为解π的行驶总里程,otr(π)为解π的时间窗准时率.若otr(π)0,则LLHi得分更新公式如式(27)所示,否则用式(28)更新LLHi得分.

3.3.3 算法高层选择策略

为增强算法的局部搜索能力,保证非劣解集周围解空间得到充分搜索,需要对其进行细致的局部搜索以保证搜索的深度;同时,为了保证局部搜索的广度,需要对一部分种群做局部搜索,充分利用种群来寻找更多的非劣解.为了在可接受时间内平衡算法局部搜索的深度与广度,通常对非劣解集和1/5到1/4的种群做局部搜索[30].

当算法运行第1代,对非劣解集和部分种群做局部搜索时,为了在同样使用次数下全面评价各个LLH的搜索性能,并且搜索出更多的非劣解供概率矩阵学习优质信息,对非劣解集和部分种群使用所有LLH构成的启发式算法HLS进行局部搜索,并在搜索过程中计算LLH得分.

当算法运行至第i(i>1)代,对非劣解集和部分种群做局部搜索时,按照LLH得分,采用轮盘赌法在LLH1至LLH11中选择4种组成启发式算法HLS,对非劣解集和种群中1/θ的个体做局部搜索,并更新LLHi得分.每产生一个可行解,就更新一次LLH的得分.

综上,表3为对一个可行解做局部搜索操作的伪代码.

表3 对一个解的局部搜索策略Table 3 Local search strategy for a solution

本节根据局部搜索操作的范围,将LLH分为2类共11种局部搜索操作(车辆内的局部搜索操作和车辆间的局部搜索操作),采用LHHLS,用对解的改进程度来评价LLH的搜索性能,并对其评分;在算法高层选择LLH阶段,依据LLH得分使用轮盘赌法选择部分算子,在选择LLH时,轮盘赌操作使得搜索能力较强的LLH能以较大概率被选中,其他LLH也能以一定概率被选中.组合选中的部分算子产生启发式算法HLS,对可行解进行高效多样的局部搜索.

综上所述,HHEDA整体伪代码如表4所示.

表4 HHEDA伪代码Table 4 HHEDA pseudo code

4 实验设计与分析

4.1 实验设置

本文测试算例采用文献[1]中测试算例集的部分算例.该算例集由在VRPTW测试算例集(即著名的所罗门测试算例集)[31]基础上生成的6种VRPPDTW测试算例集组成,是目前唯一通用的VRPPDTW测试算例集.各算法采用Python3.7编程实现,操作系统为windows10,电脑内存为16G,CPU为Intel i7–1065G7,主频1.30 GHz.

4.2 参数设置

HHEDA的关键参数为:种群的局部搜索比例1/θ、概率矩阵遗忘因子α、局部搜索策略中每种LLH重复的次数n.为考察参数设置对算法效果的影响,本文针对中等规模的测试算例rdp101(Nc100),采用实验设计方法(design of experiment,DOE)[32]进行实验分析.各关键参数设置水平表如表4所示.HHEDA的其他参数设置为:popsize100,Sstart50,Slower1,Supper200.

针对各参数的水平设置,选择规模为L16(21×42)的正交实验,采用Ishibuchi等[33]提出的方法对非支配解集进行评价,评价指标如式(30)所示:

其中:S表示在不同算法获得的非劣解集合并后的总非支配解;y ≺x表示解x被y支配;|Sj|表示第j种算法获得的非支配解集中解的个数.以R-NDS的平均值AVG作为平均响应值,AVG越大意味着这组参数下的算法性能越强,参数设置的正交表如表5所示,参数响应值如表6所示,各参数的响应趋势如图3所示.

表5 参数水平表Table 5 Parameter level table

表6 正交表和AVG统计值Table 6 Orthogonal table and AVG statistics

图3 各参数响应趋势Fig.3 The influence trend of each parameter

由表4–6和图3可知,1/θ,α和n的不同组合方式对算法性能会产生较大影响.为了发挥算法的最佳性能,对算法进行了参数实验.基于实验结果,HHEDA的参数设置为:1/θ1/5,α0.8和n3,此时,算法可表现出良好性能.因此,后文将基于此参数设置开展进一步的算法比较研究.

4.3 仿真结果与比较

为评价算法性能,采用式(30)–(31)所示的两个指标对非支配解集评价.第1个指标R-N如式(30)所示,表示Sj中的解不被S中的任意解支配的解所占的比例,衡量解的质量,第2个指标N-N如式(31)所示,评价的是Sj中的解不被S中的解支配的数量,衡量算法获得非支配解的能力.两个指标的值越大,算法的性能越好.

表7 各参数平均响应值Table 7 Average response value of each parameter

4.3.1 验证LHHLS的有效性

为验证HHEDA中局部搜索LHHLS的有效性,将其修改,形成2种算法HHEDA-V1和HHEDA-V2,并进行比较.HHEDA-V1仅采用算法运行第1代后11种LLH中评分最高的4个作为后续生成启发式算法HLS的备选集合,HHEDA-V2采用全部11种LLH的随机排序生成启发式算法HLS.对于每个算例,各算法均在相同时间(1.5×Ncs)下独立运行20次.测试结果如表8所示,其中最好指标值用粗体表示.

表8 LHHLS性能验证Table 8 LHHLS performance verification

由表8可知,HHEDA在各问题上均优于HHEDAV1和HHEDA-V2.这表明选择较大的备选集合可确保启发式算法HLS的多样性(对比HHEDA-V1),同时验证了利用各LLH在搜索过程中的贡献率来动态选择LLH构造启发式算法HLS的必要性(HHEDA-V2).因此,对于本文VRPSPDSTW这类解空间十分复杂的问题,采用LHHLS执行局部搜索是合理有效的.

4.3.2 验证HHEDA的有效性

为验证HHEDA的有效性,本节将HHEDA与现有的5 种有效算法进行比较,包括ACO[22],CRO[20],HEDA[29],MMOEASA[16],MOSA-ACO[23].对于每个测试算例,所有算法均在相同时间(秒)下独立运行20次.测试结果如表9所示,其中最好指标值用粗体表示.此外,图4–9给出了各算法在部分算例上获得的非劣解集.由于优化目标为最小化车辆行驶总里程(对应X轴)和最大化服务准时率(对应Y轴),故帕累托前沿(即最优非劣解集)趋向左上角(即坐标为(0,1)的点).

图4 算例rcdp2501非劣解集Fig.4 Example rcdp2501 non-inferior solution set

图5 算例rcdp5001非劣解集Fig.5 Example rcdp5001 non-inferior solution set

由表9可知,HHEDA在绝大多数问题上占优.除算例rcdp5001之外,HHEDA获得的非劣解集中解的个数均多于其他算法,同时HHEDA获得的非劣解集中的解不被总非劣解集中的解支配的比例明显较高.从图4–9可看出,HHEDA的分散性整体较好,覆盖区域较广,但HHEDA与HEDA,MMOEASA这3个整体性能最好的算法获得的非劣解集分布不均匀.造成这一现象的原因在于本文问题的两个优化目标(最小化车辆行驶总里程和最大化服务准时率)间存在较明显冲突,导致问题解空间中最优和优质非劣解的分布本身就不均匀,从而使得性能较好算法的非劣解集中会出现“断裂”现象.譬如,图6–7中服务准时率约在0.3～0.4间且行驶总距离约在2000～3500间的区域为各算法均未能发现非劣解的“断裂”区域.虽然问题解空间中最优和优质非劣解分布的不规则性增加了算法的搜索难度,但HHEDA利用结合问题特点设计的初始解生成规则、全局搜索和局部搜索,可实现对问题解空间更为有效的搜索,进而能发现更多趋向于图中左上角的优质非劣解(非劣解越接近坐标为(0,1)的左上角则其质量越好).仍以图6–7为例,服务准时率在0.3以下时HHEDA的非劣解基本支配了其他算法的非劣解,同时服务准时率在0.4以上时HHEDA的非劣解个数和质量整体占优.类似地,从其他图中也可看出HHEDA的非劣解在总非劣解集S中所占比例最大,覆盖区域最广.

表9 算法比较结果Table 9 Algorithm comparison results

图6 算例cdp101非劣解集Fig.6 Example cdp101 non-inferior solution set

图7 算例rcdp101非劣解集Fig.7 Example rcdp101 non-inferior solution set

具体到各算法的设计,ACO,CRO,HEDA,MMOEASA和MOSA-ACO侧重于简单采用算法自身进化机制执行全局搜索,同时采用交叉、变异等邻域操作或者较为单一的客户交换操作执行局部搜索.HHEDA进一步考虑问题特点及其解空间的复杂性,设计3种概率矩阵较全面有效地学习和积累优质解信息,并通过对其采样生成新个体以实现全局搜索,有助于尽快发现优质解区域,同时在局部搜索过程中动态构造丰富的变邻域搜索(即启发式算法HLS),有利于增加算法在优质解区域的搜索深度.因此,HHEDA对搜索方向的引导更为合理,能在相同时间内更有效地发现高质量的非劣解集.

图8 算例rdp101非劣解集Fig.8 Example rdp101 non-inferior solution set

5 结语

本文针对带软时间窗的同时取送货车辆路径问题,在考虑车辆载重等约束条件的同时,构建了最小化辆行驶总里程和最大化服务准时率的双目标模型,并提出一种超启发式分布估计算法HHEDA 进行求解.HHEDA具有如下优点:1)采用多种启发式规则共同生成初始解,可确保初始种群的质量和分散性;2)利用多个概率矩阵充分学习优质解中不同的有用信息,以推动全局搜索较快到达不同的优质解区域进行搜索;3)采用学习型超启发式局部搜索执行丰富的变邻域搜索,能有效增加搜索的深度,从而提升算法的性能.仿真实验和算法比较验证了HHEDA是求解VRPSPDSTW的有效算法.后续研究将进一步考虑道路拥堵导致速度变化的问题,并设计有效求解算法.