張烜熒,胡 蓉,錢 斌

(昆明理工大學(xué)信息工程與自動化學(xué)院自動化系,云南昆明 650500;昆明理工大學(xué)云南省人工智能重點實驗室,云南昆明 650500)

1 引言

隨著“制造強國戰(zhàn)略”的實施,在推動先進制造業(yè)迅速發(fā)展的同時,現(xiàn)代物流業(yè)的高質(zhì)量發(fā)展也成為大勢所趨.在物流系統(tǒng)中,配送是核心部分,而車輛路徑問題(vehicle routing problem,VRP)是其中的重要一環(huán).傳統(tǒng)的VRP主要描述為組織一個車隊,安排適當(dāng)?shù)男熊嚶肪€為一定數(shù)量的客戶提供送貨服務(wù),并能在滿足車載量、行車里程等約束條件下,達到諸如車輛行駛總里程最短,成本最小,耗時最少等目的.隨著資源需求的增加與供給雙向流通方式的快速發(fā)展,以及物流企業(yè)對降低成本、保護環(huán)境的注重,考慮逆向物流的車輛路徑問題開始受到重視.此外,日益激烈的市場競爭逼迫企業(yè)在滿足客戶需求的同時要考慮提高客戶滿意度.在此背景下,研究帶軟時間窗的同時取送貨車輛路徑問題(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的求解算法主要分為兩類.一類是運籌學(xué)算法,包括分支定價、動態(tài)規(guī)劃和拉格朗日松弛等[2–8],基本都是用于求解單目標(biāo)問題.此類算法以線性代數(shù)和幾何分析為基本工具,利用問題優(yōu)化目標(biāo)函數(shù)和約束式的結(jié)構(gòu)信息構(gòu)造搜索,可在幾分鐘至幾十分鐘內(nèi)獲取較小規(guī)模問題(客戶數(shù)小于等于20)的最優(yōu)解.但是,VRPSPDSTW具有非凸特性,其幾何結(jié)構(gòu)和最優(yōu)解間的關(guān)系仍屬開放問題,不存在多項式時間獲取最優(yōu)解的算法,故對應(yīng)的運籌學(xué)算法均是靠遍歷或部分遍歷解空間來確保求解質(zhì)量,這使其需要數(shù)小時至數(shù)天才可獲得較大規(guī)模問題(客戶數(shù)在80～200之間)的最優(yōu)解或滿意解.另一類是智能算法,包括遺傳算法、禁忌搜索算法、蟻群算法等[9–29],既用于求解單目標(biāo)問題,也用于求解多目標(biāo)問題.此類算法不依賴于問題結(jié)構(gòu),而是將問題的約束條件隱式地包含在所設(shè)計的解的編碼、解碼規(guī)則中來處理,并利用某種擬人、擬物的機制不斷生成新的可行個體或解,從而引導(dǎo)算法執(zhí)行搜索,往往可在幾秒或十幾秒內(nèi)就能獲得不同規(guī)模問題的滿意解.因此,設(shè)計智能算法以實現(xiàn)對VRPSPDSTW的快速有效求解,是合理且必要的.

近年來,智能算法求解VRPSPDSTW逐漸受到重視.針對單目標(biāo)VRPSPDTW,Lai等[9]提出改進差分進化算法(improved differential evolution algorithm,IDE)求解問題,取得較好的效果.Wang等[1]提出一種基于最小代價插入法的協(xié)同進化遺傳算法(genetic algorithm,GA)進行求解,實驗結(jié)果表明所提算法的有效性.王等[10]提出模擬退火算法(simulated annealing,SA),并在模擬退火過程中融合局部搜索方法,在可接受的計算時間內(nèi)取得了較好的優(yōu)化結(jié)果.王等[11]提出一種離散布谷鳥算法(discrete cuckoo search,DCS),其L′evy飛行機制協(xié)調(diào)局部搜索和全局搜索,實驗結(jié)果更新了國際最優(yōu)解.韓等[12]針對帶折線型軟時間窗的車輛路徑問題(vehicle routing problem with soft time windows,VRPSTW),以最小化運輸總成本為目標(biāo),提出一種超啟發(fā)式遺傳算法(hyper-heuristic genetic algorithm,HHGA),以GA 作為上層搜索算法,CW 節(jié)約法、MJ插入法、Kilby插入法作為底層搜索規(guī)則,并通過加入排序規(guī)則和解內(nèi)部調(diào)整算法以進一步提高解的質(zhì)量.實驗結(jié)果驗證了HHGA 的可行性和有效性.王等[13]提出一種回溯搜索優(yōu)化算法(backtracking search optimization algorithm,BSA),使用啟發(fā)式規(guī)則產(chǎn)生初始種群,通過多種搜索算子更新局部最優(yōu)解,實驗驗證了BSA的有效性.Qin等[14]進一步考慮了低碳需求,提出一種自適應(yīng)遺傳爬山算法(adaptive genetic hill-climbing algorithm,AGHCA)進行求解,同時研究了碳稅對總成本和碳排放的影響.

對于多目標(biāo)VRPSPDSTW(multi-objective VRPSPDSTW,MOVRPSPDSTW)及相關(guān)問題,李等[15]針對VRPSPDTW,以最小化總車輛數(shù)、運輸總時間、車輛行駛總里程的線性組合為優(yōu)化目標(biāo),提出一種基于路線集合劃分的分解迭代算法,將路線集合分解為路線子集合,再用記錄更新法求解子集合,不斷迭代逐漸逼近最優(yōu)解,實驗結(jié)果表明了算法的有效性.Banos等[16]針對考慮負載均衡的帶硬時間窗的車輛路徑問題(vehicle routing problem with time windows,VRPTW),以同時最小化各車輛間最大行駛距離之差和車輛間最大負載之差為優(yōu)化目標(biāo),提出一種結(jié)合模擬退火算法的多起點多目標(biāo)進化算法(multi-start multiobjective evolutionary algorithm with simulated annealing,MMOEASA).仿真實驗和對比結(jié)果驗證了算法的有效性.Yang等[17]針對帶軟時間窗的車輛路徑問題(VRPSTW),以最小化總運營成本、時間窗偏離總量、碳排放為優(yōu)化目標(biāo),提出一種混合遺傳算法(hybrid genetic algorithm,HGA),通過實驗驗證了算法的有效性,并通過帕累托分析,說明了目標(biāo)之間的關(guān)系.Sumaiya等[18]針對VRPSTW,以最小化總行駛距離、違反時間窗的服務(wù)次數(shù)、車輛數(shù)為優(yōu)化目標(biāo),提出一種人工蜂群算法(artificial bee colony,ABC),結(jié)合車輛內(nèi)和車輛間的兩步局部搜索進行求解,實驗結(jié)果表明ABC能在合理的計算時間內(nèi)獲得高質(zhì)量的解.Zhang等[19]針對時間依賴型同時取送貨車輛路徑問題(time dependent vehicle routing problem with simultaneous delivery and pickup,TDVRPSDP),以最小化碳排放和旅行時間為優(yōu)化目標(biāo),提出一種多目標(biāo)協(xié)同進化量子算法(multi-objective cooperative quantum evolutionary algorithm,MCQEA),將種群分為兩個亞種群,使用非劣解集和亞種群中的解共同產(chǎn)生新種群,實驗結(jié)果驗證了MCQEA的有效性.Li等[20]針對MOVRPSPDSTW,以最小化車輛數(shù)、車輛固定成本、行駛總距離、旅行時間、服務(wù)成本為優(yōu)化目標(biāo),使用基于分解的多目標(biāo)化學(xué)反應(yīng)優(yōu)化算法(chemical reaction optimization,CRO),將復(fù)雜的MOVRPSPDSTW轉(zhuǎn)換為多個單目標(biāo)子問題,并同時在CRO中對其進行優(yōu)化,實驗結(jié)果驗證了算法的有效性.陳等[21]針對多目標(biāo)同時取送貨車輛路徑問題(multi-objective vehicle routing problem with simultaneous pickup and delivery,MOVRPSPD),以最小化總路徑成本、各車輛最大行駛里程差為優(yōu)化目標(biāo),提出一種嵌入禁忌表、且具有貪婪轉(zhuǎn)移準(zhǔn)則的多目標(biāo)改進蟻群算法(improved ant colony optimization,IACO),實驗結(jié)果表明IACO可求得權(quán)衡各優(yōu)化目標(biāo)且使單一優(yōu)化目標(biāo)近似最優(yōu)的Pareto解.柴等[22]考慮路徑具有雙重屬性,建立基于雙屬性的多目標(biāo)VRPTW,進而提出一種蟻群算法(ant colony algorithm,ACO),在算法中定義了一種綜合考慮各優(yōu)化目標(biāo)、時間窗和信息素等啟發(fā)信息的狀態(tài)轉(zhuǎn)移概率公式.實驗結(jié)果表明,所提算法在運行時間、收斂性和群體多樣性方面具有優(yōu)勢.Wang等[23]針對帶時間窗的周期性車輛路徑問題(periodic vehicle routing problem with time windows,PVRPTW),以同時最小化車輛行駛總距離、最小化車輛總數(shù)、最大化客戶訪問頻率為優(yōu)化目標(biāo),提出一種多目標(biāo)模擬退火蟻群優(yōu)化算法(multi-objective simululate annealing-ant colony optimization,MOSAACO),在算法中設(shè)計了基于歐氏距離的優(yōu)質(zhì)解接受方式和可生成不同初始解的參數(shù)控制策略,并加入了有效局部搜索策略.實驗表明,該算法在求解PVRPTW問題上性能良好.由文獻調(diào)研可知,雖然VRPSPDSTW已得到一定研究,但尚無人考慮客戶滿意度(服務(wù)準(zhǔn)時率)這一重要優(yōu)化目標(biāo).此外,上述算法在全局搜索時基本都直接采用算法原有搜索機制執(zhí)行,未考慮結(jié)合問題特性做合理更改來增強算法發(fā)現(xiàn)優(yōu)質(zhì)解的能力,且局部搜索時大多迭代執(zhí)行較為單一的插入、變換操作,難以實現(xiàn)對復(fù)雜問題解空間中優(yōu)質(zhì)解區(qū)域的深入搜索.顯然,針對VRPSPDSTW這類復(fù)雜問題,如何設(shè)計能獲取高質(zhì)量解的有效算法是需要進一步研究的課題.

分布估計算法(estimation of distribution algorithm,EDA)是一種新興的基于統(tǒng)計學(xué)原理的隨機優(yōu)化算法.EDA采用統(tǒng)計學(xué)習(xí)手段從群體宏觀角度建立概率模型,用來學(xué)習(xí)和積累優(yōu)質(zhì)解在搜索空間的分布信息,然后對此概率模型采樣產(chǎn)生新種群,如此反復(fù)實現(xiàn)種群進化[24].此進化方式可一定程度上避免傳統(tǒng)進化算法中交叉、變異等操作所帶來的模式破壞問題,從而使其搜索具有更好的引導(dǎo)性[25].目前,使用EDA求解VRP的研究較為有限.Wan等[26]針對VRPSTW,提出一種混合粒子群算法的混合EDA,通過粒子群算法的快速更新策略修改了EDA的概率矩陣,實驗結(jié)果驗證了算法的有效性.孟等[27]針對帶硬時間窗的車輛路徑問題(VRPTW),提出一種混合種群增量學(xué)習(xí)算法(hybrid population-based incremental learning algorithm,HPBIL),為每輛車建立2維EDA概率模型,同時設(shè)計了基于插入和兩點鄰域交換的兩階段局部搜索來增強算法的局部開發(fā)能力,實驗結(jié)果表明混合后的算法求解效果優(yōu)良.Wang等[28]針對客戶需求和運輸成本不確定的車輛路徑問題(uncertain capacitated vehicle routing problems,UCVRP),提出一種改進EDA,不僅利用概率模型學(xué)習(xí)優(yōu)質(zhì)解信息并引導(dǎo)全局搜索,而且引入插入操作和2-opt操作來提高局部搜索能力,同時采用多場景技術(shù),增強算法求解不確定問題的魯棒性.Ji等[29]針對考慮碳排放的多目標(biāo)聯(lián)合運輸VRP,提出一種混合EDA(hybrid estimation of distribution algorithm,HEDA),采用擴展的2維概率矩陣學(xué)習(xí)各階段運輸模式對應(yīng)的優(yōu)質(zhì)非支配解信息,同時設(shè)計了基于部分交換操作的局部搜索,對問題進行了有效求解.從EDA在VRP相關(guān)問題上的研究現(xiàn)狀來看,采用基于問題的概率模型可以合理引導(dǎo)全局搜索的方向,引入有效的局部搜索策略能進一步提高算法的性能,兩者的協(xié)同搜索是算法設(shè)計的關(guān)鍵.據(jù)文獻調(diào)研可知,針對VRPSPDSTW,尚無EDA求解的報道,故十分必要開展相關(guān)研究.

綜上,本文研究多目標(biāo)VRPSPDSTW的建模與求解.建模方面,建立同時以最小化總行駛路徑和最大化服務(wù)準(zhǔn)時率為優(yōu)化目標(biāo)的VRPSPDSTW.求解方面,在對問題模型特點分析的基礎(chǔ)上,提出一種超啟發(fā)式分布估計算法(hyper-heuristic estimation of distribution algorithm,HHEDA).HHEDA在全局搜索階段執(zhí)行分布估計算法EDA,在局部搜索階段執(zhí)行學(xué)習(xí)型超啟發(fā)式局部搜索(learning hyper-heuristic local search,LHHLS).其中,EDA利用結(jié)合問題特性所設(shè)計的3類概率矩陣協(xié)同學(xué)習(xí)和積累優(yōu)質(zhì)解(即可更新當(dāng)前非劣解集的解)信息并引導(dǎo)搜索盡快到達優(yōu)質(zhì)解區(qū)域;LHHLS采用學(xué)習(xí)型評分機制和11種有效鄰域操作構(gòu)造一系列啟發(fā)式算法或搜索,然后分別對當(dāng)前種群的部分個體和當(dāng)前非劣解集中的解逐一執(zhí)行相應(yīng)的啟發(fā)式搜索,可進一步加大搜索深度,從而增大算法在較短時間內(nèi)獲取高質(zhì)量解的機率.最后,通過對不同測試問題的仿真實驗和算法比較,驗證了HHEDA的有效性.

2 問題描述

關(guān)于本文所涉及的數(shù)學(xué)符號及定義如表1所示,VRPSPDSTW示例如圖1所示.定義一個無向圖G(V,E),其中,V為節(jié)點集且VVc ∪{0},E為路徑集.車輛從配送中心出發(fā)對若干客戶進行服務(wù),在約束條件下同時滿足客戶取送貨需求.要求合理安排路徑,最小化車輛行駛總里程和最大化服務(wù)準(zhǔn)時率.

表1 符號表Table 1 Symbol table

圖1 VRPSPDSTW問題示例Fig.1 VRPSPDSTW problem example

根據(jù)上述描述,建立如下VRPSPDSTW的問題模型:

其中:式(1)–(2)為優(yōu)化目標(biāo),式(1)表示最小化車輛行駛總里程,式(2)表示最大化服務(wù)準(zhǔn)時率;式(3)–(12)為約束條件,式(3)表示每個客戶只被一輛車服務(wù)一次,式(4)表示客戶點車輛流平衡約束,即車輛到達客戶點h后必須離開該客戶點;式(5)–(6)表示每輛車起始于配送中心并且完成服務(wù)后最終回到配送中心;式(7)–(8)表示車輛從配送中心出發(fā)時的載重,且該載重不能超過車輛最大載重;式(9)–(10)表示車輛服務(wù)完客戶j的載重,且該載重不允許超過車輛最大載重;式(11)表示車輛從i出發(fā)服務(wù)j,到達j的時間.若車輛僅服務(wù)客戶i,且滿足式(12),那么這個客戶可以被準(zhǔn)時服務(wù).

由上述問題模型可知,相對于傳統(tǒng)VRP,本文的VRPSPDSTW進一步考慮軟時間窗和同時取送貨,使得模型約束增多,問題的解空間更加復(fù)雜;而且,增加了優(yōu)化目標(biāo)(即最大化服務(wù)準(zhǔn)時率),使得解空間中優(yōu)質(zhì)解(即非劣解)的分布更加分散;此外,隨著客戶數(shù)量的增加,VRPSPDSTW的解空間規(guī)模呈指數(shù)增長[11].這些表明求解VRPSPDSTW的難度很大.

3 HHEDA

HHEDA的框架如圖2所示.從圖2可知,在HHEDA的初始化階段,為提高初始種群的質(zhì)量并加快算法收斂速度,提出3條規(guī)則用于生成初始解.在HHEDA 的全局搜索階段,執(zhí)行分布估計算法EDA.在EDA中,結(jié)合問題特性,設(shè)計排序概率矩陣、距離概率矩陣、捆綁關(guān)系概率矩陣共同學(xué)習(xí)和積累優(yōu)質(zhì)解(即可更新當(dāng)前非劣解集的解)信息,然后通過采樣這些概率矩陣以生成較高質(zhì)量的新個體或解,并更新非劣解集和概率矩陣,可較快引導(dǎo)算法到達解空間中存在優(yōu)質(zhì)解的區(qū)域.在HHEDA的局部搜索階段,執(zhí)行學(xué)習(xí)型超啟發(fā)式局部搜索LHHLS,以增強算法的局部尋優(yōu)能力.在LHHLS中,低層問題域采用11種有效鄰域操作組成備選集合,高層控制域設(shè)計改進的學(xué)習(xí)型評分機制動態(tài)評價底層不同鄰域操作在搜索過程中的貢獻(即對兩個優(yōu)化目標(biāo)的改進程度),進而依據(jù)各鄰域操作的貢獻率,對全局搜索階段EDA種群中部分個體和非劣解集中的解,逐一執(zhí)行由4種鄰域操作(用輪盤賭方法確定)構(gòu)成的啟發(fā)式算法,以實現(xiàn)較深入的局部搜索并更新當(dāng)前非劣解集和概率矩陣.

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

本節(jié)后續(xù)對HHEDA各環(huán)節(jié)進行介紹.

3.1 解的表示

令解π的長度為L,π(π[1],π[2],···,π[L]),其中Nc+2 ≤L≤2×Nc+1.配送中心編號為0,客戶編號為大于0的整數(shù).車輛從配送中心開出前往客戶處,服務(wù)一定數(shù)量的客戶后返回配送中心.如圖1所示,編碼為π(0,1,2,3,4,0,5,6,7,0,8,9,10,0)的解表示車輛1從配送中心開出,依次服務(wù)客戶1、客戶2、客戶3、客戶4后返回配送中心.車輛2、車輛3的服務(wù)順序釋義以此類推.

3.2 全局搜索

3.2.1 種群初始化

為了保證種群的多樣性和分散性,使用規(guī)則和隨機方法產(chǎn)生規(guī)模為popsize的第一代種群.使用規(guī)則1和規(guī)則2產(chǎn)生2個解,剩下popsize-2個解使用規(guī)則3 產(chǎn)生.令未被服務(wù)的客戶集合為C.具體規(guī)則如下:

規(guī)則1最近鄰規(guī)則.

步驟如果車輛從配送中心出發(fā),計算配送中心與C中每個客戶的距離,選擇最近的一個客戶進行服務(wù),如果有多個最近客戶,隨機選擇其中一個進行服務(wù),并將這個客戶從C中去除;如果車輛從客戶i出發(fā),計算i與C中客戶的距離,選擇最近的一個客戶進行服務(wù),如果有多個最近客戶,隨機選擇其中一個進行服務(wù),如果車輛繼續(xù)服務(wù)該客戶不會違反約束,則將這個城市從C中去除,否則,車輛放棄服務(wù)該客戶并返回配送中心.

規(guī)則2準(zhǔn)時率最高規(guī)則.

步驟如果車輛從配送中心出發(fā),計算C中每個客戶的到達時間,選擇在時間窗內(nèi)到達的客戶進行服務(wù),如果有多個可在時間窗內(nèi)到達的客戶,選擇其中最近的進行服務(wù),如果沒有可在時間窗內(nèi)到達的客戶,在C中選擇距離最近的客戶進行服務(wù),將選中的客戶從C中去除;如果車輛從客戶i出發(fā),計算車輛從i到C中每個客戶的時間,選擇在時間窗內(nèi)到達的客戶進行服務(wù),如果有多個可在時間窗內(nèi)到達的客戶,選擇其中最近的進行服務(wù),如果車輛繼續(xù)服務(wù)該客戶不會違反約束,則將這個城市從C中去除,否則,車輛放棄服務(wù)該客戶并返回配送中心,如果沒有可在時間窗內(nèi)到達的客戶,車輛返回配送中心.

規(guī)則3隨機規(guī)則.

步驟在C中隨機選取客戶進行服務(wù),如果車輛服務(wù)該客戶不會違反約束條件,那么將這個客戶從C中去除,否則,車輛放棄服務(wù)該客戶并返回配送中心.

3.2.2 概率矩陣更新

在設(shè)計概率矩陣時,通過研究問題的特性,引入問題的特定信息和知識,是提高算法搜索能力的有效途徑.對于VRPSPDSTW,為了最小化車輛行駛總里程,車輛總是傾向于先服務(wù)那些離當(dāng)前位置比較近的客戶,所以在車輛行駛總里程較小的解中,存在一些由客戶間距離決定的位置信息;為了最大化服務(wù)準(zhǔn)時率,車輛總是傾向于按照時間窗先后順序服務(wù)客戶,若客戶i,j的時間窗和距離滿足式(13)的關(guān)系,那么當(dāng)車輛連續(xù)服務(wù)客戶i,j時,只要車輛在時間窗內(nèi)到達i,那么j一定能在時間窗內(nèi)被服務(wù),即ij滿足捆綁關(guān)系.所以,在服務(wù)準(zhǔn)時率較高的解中,客戶間由于時間窗和距離較為匹配,存在一些捆綁關(guān)系信息.在HHEDA中,用ρ,λ,τ三個矩陣分別學(xué)習(xí)優(yōu)質(zhì)解的排序信息、客戶間的距離信息和捆綁關(guān)系.

概率矩陣初始化階段,ρ和λ為Nc×Nc的全零矩陣,τ為Nc×Nc的矩陣,用式(14)對τ中元素值進行初始化.在概率矩陣更新時,選取非劣解集中新增加的個體集合?按式(16)–(20)更新概率矩陣.其中,在每張車輛服務(wù)順序中,ρj′j表示客戶j出現(xiàn)在位置j′的次數(shù);Kij為決策變量,當(dāng)客戶i與j被連續(xù)服務(wù)時為1,否則為0;Ij′j為決策變量,當(dāng)客戶j出現(xiàn)在位置j′時為1,否則為0;γij為學(xué)習(xí)速率,其值為客戶ij之間距離的倒數(shù);Lij為決策變量,當(dāng)客戶i與j被連續(xù)服務(wù)并且車輛到達i,j的時間都在時間窗內(nèi)時其值為1,否則為0.

3.2.3 概率矩陣遺忘操作

在迭代過程中,新產(chǎn)生的非劣解會支配掉部分原非劣解集中的解,為了使從新非劣解中學(xué)習(xí)到的信息在概率矩陣中占較大比重,在算法每一次迭代結(jié)束后對概率矩陣按式(21)–(23)執(zhí)行遺忘操作,遺忘因子為α(α<1)

3.2.4 對概率模型采樣產(chǎn)生種群

算法的第iteration(iteration ≥2)代,新種群由精英保留機制和采樣概率矩陣產(chǎn)生的解共同構(gòu)成.隨機選擇m個非劣解(mmin(ps,「popsize×15%」),ps為非劣解集中解的個數(shù))作為精英解集,剩余popsize?m個解由概率矩陣采樣生成.

采樣過程中,根據(jù)式(24)–(26),考慮客戶排序信息、距離信息和捆綁關(guān)系,產(chǎn)生較優(yōu)個體.其中C為未被服務(wù)的客戶集合,(i ?1)′表示在車輛中被第i ?1個被服務(wù)的客戶,pij為客戶j在車輛中被第i個服務(wù)的概率,每次采樣均使用輪盤賭操作,在一定程度上保證了種群的多樣性.

當(dāng)產(chǎn)生位置i1的客戶時,由于是車輛第1個服務(wù)的客戶,所以僅考慮客戶在解中的排序位置,由式(24)計算C中客戶被選中的概率;當(dāng)產(chǎn)生位置i(i≥2)的客戶時,若客戶(i ?1)′被準(zhǔn)時服務(wù),則在計算C中客戶被選中的概率時,著重考慮客戶之間的捆綁關(guān)系,以最大化服務(wù)準(zhǔn)時率為側(cè)重優(yōu)化目標(biāo),由客戶排序信息和捆綁關(guān)系綜合計算概率,計算公式如式(25)所示;若客戶(i ?1)′沒有被準(zhǔn)時服務(wù),著重考慮客戶之間的距離關(guān)系,以最小化總距離為側(cè)重優(yōu)化目標(biāo),由客戶排序信息和距離關(guān)系綜合計算概率,計算公式如式(26)所示:

綜上,算法任意代數(shù)的種群初始化的偽代碼如表2所示.

表2 種群初始化Table 2 Population initialization

3.3 局部搜索

對于多目標(biāo)組合優(yōu)化問題,全局最優(yōu)非劣解集中的任何解(即全局最優(yōu)非劣解)在所有鄰域結(jié)構(gòu)下均為局部最優(yōu)非劣解,而接近全局最優(yōu)非劣解的優(yōu)質(zhì)解往往在多種鄰域結(jié)構(gòu)下均為局部最優(yōu)非劣解.采用單一鄰域操作(對應(yīng)一種鄰域結(jié)構(gòu))的迭代搜索容易較早達到并陷入該鄰域結(jié)構(gòu)的局部最優(yōu)非劣解,但該解的質(zhì)量大都一般;而采用多種有效鄰域操作(對應(yīng)多種鄰域結(jié)構(gòu))的迭代搜索可在到達多種鄰域結(jié)構(gòu)共同的局部最優(yōu)非劣解前一直持續(xù)進行搜索,從而容易獲取優(yōu)質(zhì)非劣解.此外,多目標(biāo)組合優(yōu)化問題的解空間十分復(fù)雜,大量優(yōu)質(zhì)非劣解分散在龐大的“峽谷型”解空間中靠近底部的多個區(qū)域[30].因此,為能對HHEDA全局搜索發(fā)現(xiàn)的優(yōu)質(zhì)解區(qū)域(即非劣解集)進行有效搜索,本文設(shè)計學(xué)習(xí)型超啟發(fā)式局部搜索(LHHLS),用于執(zhí)行豐富的多種變鄰域搜索,以增強算法局部搜索的深度.

超啟發(fā)算法為兩層結(jié)構(gòu),底層問題域包含一組直接搜索問題解空間的啟發(fā)式算法(low-level heuristics,LLHs),其中每個LLH可為啟發(fā)式算法,也可以是規(guī)則或鄰域操作;高層策略域包含確定底層不同LLH組合的某種算法.高層策略域的算法(high-level strategy,HLS)通過動態(tài)控制底層問題域的LLHs來不斷組合生產(chǎn)新的啟發(fā)式算法,進而使用這些新的啟發(fā)式算法在底層搜索問題解空間.對于本文的LHHLS,底層問題域設(shè)計11種LLH,每種LLH為一類有效鄰域操作;高層策略域的算法由改進的學(xué)習(xí)型評分機制和輪盤賭方法共同構(gòu)成.前者用于定量地評估和統(tǒng)計每個LLH在搜索過程中的執(zhí)行效果;后者依據(jù)所有LLHs的評分(即貢獻)對其進行輪盤賭,以生成一系列新的啟發(fā)式算法(即LLH的組合).

顯然,搜索能力較強的LLH能以較大概率被選中,其他LLH也能以一定概率被選中,這有利于生成多種新的有效啟發(fā)式算法.

下面對LHHLS各環(huán)節(jié)進行介紹.

3.3.1 LLH設(shè)計

LLH要根據(jù)具體問題進行設(shè)計.本文設(shè)計了2類共11種鄰域變換操作作為低層啟發(fā)式算子,LLH1至LLH5為車輛內(nèi)路徑局部搜索操作,LLH6–LLH11為車輛間路徑局部搜索操作.

1) LLH1:車輛內(nèi)客戶順序逆轉(zhuǎn)操作.選擇被同一張車連續(xù)服務(wù)的幾個客戶,進行逆轉(zhuǎn).

2) LLH2:車輛內(nèi)不相鄰客戶交換操作.選擇被同一張車服務(wù)的兩個不連續(xù)的客戶,進行交換.

3) LLH3:車輛內(nèi)相鄰客戶交換操作.選擇被同一張車服務(wù)的兩個相鄰的客戶,進行交換.

4) LLH4:車輛內(nèi)客戶插入操作.選擇一張車輛的服務(wù)路線ri,在ri中隨機選擇一個城市c,將其插入到ri中的其他位置.

5) LLH5:車輛內(nèi)客戶基于準(zhǔn)時性的交換操作.選擇一張車輛的服務(wù)路線ri,在ri中選擇車輛晚于時間窗下限到達的城市c1和車輛早于時間窗上限到達的城市c2,交換c1和c2.

6) LLH6:車輛間客戶插入操作.選中兩張車輛的服務(wù)路線ri和rj,在ri中隨機選擇一個城市c,將其插入到rj中.

7) LLH7:車輛間客戶交換操作.選中兩張車輛的服務(wù)路線ri和rj,在ri中隨機選擇一個城市c1,在rj中隨機選擇一個城市c2,交換c1和c2.

8) LLH8:車輛間連續(xù)客戶插入操作.選中兩張車輛的服務(wù)路線ri和rj,在ri中隨機選擇兩個被連續(xù)服務(wù)的城市c1和c2,將c1和c2插入到rj中.

9) LLH9:車輛間連續(xù)客戶交換操作.選中兩張車輛的服務(wù)路線ri和rj,在ri中隨機選擇兩個被連續(xù)服務(wù)的城市c1和c2,在rj中隨機選擇兩個被連續(xù)服務(wù)的城市c3和c4,交換c1,c2和c3,c4.

10) LLH10:車輛間客戶交換操作.選中兩張車輛的服務(wù)路線ri和rj,在ri中隨機選擇兩個連續(xù)服務(wù)的城市c1和c2,在rj中隨機選擇一個城市c3,交換c1,c2和c3.

11) LLH11:車輛間客戶基于準(zhǔn)時性的交換操作.選擇兩張車輛服務(wù)路線ri和rj,在ri中選擇車輛晚于時間窗下限到達的城市c1,在rj中選擇車輛早于時間窗上限到達的城市c2,交換c1和c2.

在對解執(zhí)行任意一種LLHi(1 ≤i≤11)操作后,若產(chǎn)生的新解為可行解且支配原解,則用新解替代原解,否則,不接受該新解.

3.3.2 基于改進程度的學(xué)習(xí)型評分機制

對于雙目標(biāo)優(yōu)化問題,當(dāng)前解進行LLHi(1 ≤i≤11)操作后得到的可行解存在3種可能:1)新解所有的目標(biāo)都得到改進或一個目標(biāo)得到改進,另一個目標(biāo)保持不變;2)新解的一個目標(biāo)得到改進,另一個目標(biāo)變差了;3)新解的兩個目標(biāo)都變差了.每次使用LLHi對解π進行局部搜索后,采用所設(shè)計的評分機制(式(27)–(29))定量學(xué)習(xí)和積累LLHi的貢獻率或得分信息,并依此采用輪盤賭方法構(gòu)造新的啟發(fā)式算法HLS(即局部搜索策略),用于對下一代的各解進行局部搜索.LLHi初始得分為Sstart,最高得分為Supper,最低得分為Slower.l(π)為解π的行駛總里程,otr(π)為解π的時間窗準(zhǔn)時率.若otr(π)0,則LLHi得分更新公式如式(27)所示,否則用式(28)更新LLHi得分.

3.3.3 算法高層選擇策略

為增強算法的局部搜索能力,保證非劣解集周圍解空間得到充分搜索,需要對其進行細致的局部搜索以保證搜索的深度;同時,為了保證局部搜索的廣度,需要對一部分種群做局部搜索,充分利用種群來尋找更多的非劣解.為了在可接受時間內(nèi)平衡算法局部搜索的深度與廣度,通常對非劣解集和1/5到1/4的種群做局部搜索[30].

當(dāng)算法運行第1代,對非劣解集和部分種群做局部搜索時,為了在同樣使用次數(shù)下全面評價各個LLH的搜索性能,并且搜索出更多的非劣解供概率矩陣學(xué)習(xí)優(yōu)質(zhì)信息,對非劣解集和部分種群使用所有LLH構(gòu)成的啟發(fā)式算法HLS進行局部搜索,并在搜索過程中計算LLH得分.

當(dāng)算法運行至第i(i>1)代,對非劣解集和部分種群做局部搜索時,按照LLH得分,采用輪盤賭法在LLH1至LLH11中選擇4種組成啟發(fā)式算法HLS,對非劣解集和種群中1/θ的個體做局部搜索,并更新LLHi得分.每產(chǎn)生一個可行解,就更新一次LLH的得分.

綜上,表3為對一個可行解做局部搜索操作的偽代碼.

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

本節(jié)根據(jù)局部搜索操作的范圍,將LLH分為2類共11種局部搜索操作(車輛內(nèi)的局部搜索操作和車輛間的局部搜索操作),采用LHHLS,用對解的改進程度來評價LLH的搜索性能,并對其評分;在算法高層選擇LLH階段,依據(jù)LLH得分使用輪盤賭法選擇部分算子,在選擇LLH時,輪盤賭操作使得搜索能力較強的LLH能以較大概率被選中,其他LLH也能以一定概率被選中.組合選中的部分算子產(chǎn)生啟發(fā)式算法HLS,對可行解進行高效多樣的局部搜索.

綜上所述,HHEDA整體偽代碼如表4所示.

表4 HHEDA偽代碼Table 4 HHEDA pseudo code

4 實驗設(shè)計與分析

4.1 實驗設(shè)置

本文測試算例采用文獻[1]中測試算例集的部分算例.該算例集由在VRPTW測試算例集(即著名的所羅門測試算例集)[31]基礎(chǔ)上生成的6種VRPPDTW測試算例集組成,是目前唯一通用的VRPPDTW測試算例集.各算法采用Python3.7編程實現(xiàn),操作系統(tǒng)為windows10,電腦內(nèi)存為16G,CPU為Intel i7–1065G7,主頻1.30 GHz.

4.2 參數(shù)設(shè)置

HHEDA的關(guān)鍵參數(shù)為:種群的局部搜索比例1/θ、概率矩陣遺忘因子α、局部搜索策略中每種LLH重復(fù)的次數(shù)n.為考察參數(shù)設(shè)置對算法效果的影響,本文針對中等規(guī)模的測試算例rdp101(Nc100),采用實驗設(shè)計方法(design of experiment,DOE)[32]進行實驗分析.各關(guān)鍵參數(shù)設(shè)置水平表如表4所示.HHEDA的其他參數(shù)設(shè)置為:popsize100,Sstart50,Slower1,Supper200.

針對各參數(shù)的水平設(shè)置,選擇規(guī)模為L16(21×42)的正交實驗,采用Ishibuchi等[33]提出的方法對非支配解集進行評價,評價指標(biāo)如式(30)所示:

其中:S表示在不同算法獲得的非劣解集合并后的總非支配解;y ?x表示解x被y支配;|Sj|表示第j種算法獲得的非支配解集中解的個數(shù).以R-NDS的平均值A(chǔ)VG作為平均響應(yīng)值,AVG越大意味著這組參數(shù)下的算法性能越強,參數(shù)設(shè)置的正交表如表5所示,參數(shù)響應(yīng)值如表6所示,各參數(shù)的響應(yīng)趨勢如圖3所示.

表5 參數(shù)水平表Table 5 Parameter level table

表6 正交表和AVG統(tǒng)計值Table 6 Orthogonal table and AVG statistics

圖3 各參數(shù)響應(yīng)趨勢Fig.3 The influence trend of each parameter

由表4–6和圖3可知,1/θ,α和n的不同組合方式對算法性能會產(chǎn)生較大影響.為了發(fā)揮算法的最佳性能,對算法進行了參數(shù)實驗.基于實驗結(jié)果,HHEDA的參數(shù)設(shè)置為:1/θ1/5,α0.8和n3,此時,算法可表現(xiàn)出良好性能.因此,后文將基于此參數(shù)設(shè)置開展進一步的算法比較研究.

4.3 仿真結(jié)果與比較

為評價算法性能,采用式(30)–(31)所示的兩個指標(biāo)對非支配解集評價.第1個指標(biāo)R-N如式(30)所示,表示Sj中的解不被S中的任意解支配的解所占的比例,衡量解的質(zhì)量,第2個指標(biāo)N-N如式(31)所示,評價的是Sj中的解不被S中的解支配的數(shù)量,衡量算法獲得非支配解的能力.兩個指標(biāo)的值越大,算法的性能越好.

表7 各參數(shù)平均響應(yīng)值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個作為后續(xù)生成啟發(fā)式算法HLS的備選集合,HHEDA-V2采用全部11種LLH的隨機排序生成啟發(fā)式算法HLS.對于每個算例,各算法均在相同時間(1.5×Ncs)下獨立運行20次.測試結(jié)果如表8所示,其中最好指標(biāo)值用粗體表示.

表8 LHHLS性能驗證Table 8 LHHLS performance verification

由表8可知,HHEDA在各問題上均優(yōu)于HHEDAV1和HHEDA-V2.這表明選擇較大的備選集合可確保啟發(fā)式算法HLS的多樣性(對比HHEDA-V1),同時驗證了利用各LLH在搜索過程中的貢獻率來動態(tài)選擇LLH構(gòu)造啟發(fā)式算法HLS的必要性(HHEDA-V2).因此,對于本文VRPSPDSTW這類解空間十分復(fù)雜的問題,采用LHHLS執(zhí)行局部搜索是合理有效的.

4.3.2 驗證HHEDA的有效性

為驗證HHEDA的有效性,本節(jié)將HHEDA與現(xiàn)有的5 種有效算法進行比較,包括ACO[22],CRO[20],HEDA[29],MMOEASA[16],MOSA-ACO[23].對于每個測試算例,所有算法均在相同時間(秒)下獨立運行20次.測試結(jié)果如表9所示,其中最好指標(biāo)值用粗體表示.此外,圖4–9給出了各算法在部分算例上獲得的非劣解集.由于優(yōu)化目標(biāo)為最小化車輛行駛總里程(對應(yīng)X軸)和最大化服務(wù)準(zhǔn)時率(對應(yīng)Y軸),故帕累托前沿(即最優(yōu)非劣解集)趨向左上角(即坐標(biāo)為(0,1)的點).

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

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

由表9可知,HHEDA在絕大多數(shù)問題上占優(yōu).除算例rcdp5001之外,HHEDA獲得的非劣解集中解的個數(shù)均多于其他算法,同時HHEDA獲得的非劣解集中的解不被總非劣解集中的解支配的比例明顯較高.從圖4–9可看出,HHEDA的分散性整體較好,覆蓋區(qū)域較廣,但HHEDA與HEDA,MMOEASA這3個整體性能最好的算法獲得的非劣解集分布不均勻.造成這一現(xiàn)象的原因在于本文問題的兩個優(yōu)化目標(biāo)(最小化車輛行駛總里程和最大化服務(wù)準(zhǔn)時率)間存在較明顯沖突,導(dǎo)致問題解空間中最優(yōu)和優(yōu)質(zhì)非劣解的分布本身就不均勻,從而使得性能較好算法的非劣解集中會出現(xiàn)“斷裂”現(xiàn)象.譬如,圖6–7中服務(wù)準(zhǔn)時率約在0.3～0.4間且行駛總距離約在2000～3500間的區(qū)域為各算法均未能發(fā)現(xiàn)非劣解的“斷裂”區(qū)域.雖然問題解空間中最優(yōu)和優(yōu)質(zhì)非劣解分布的不規(guī)則性增加了算法的搜索難度,但HHEDA利用結(jié)合問題特點設(shè)計的初始解生成規(guī)則、全局搜索和局部搜索,可實現(xiàn)對問題解空間更為有效的搜索,進而能發(fā)現(xiàn)更多趨向于圖中左上角的優(yōu)質(zhì)非劣解(非劣解越接近坐標(biāo)為(0,1)的左上角則其質(zhì)量越好).仍以圖6–7為例,服務(wù)準(zhǔn)時率在0.3以下時HHEDA的非劣解基本支配了其他算法的非劣解,同時服務(wù)準(zhǔn)時率在0.4以上時HHEDA的非劣解個數(shù)和質(zhì)量整體占優(yōu).類似地,從其他圖中也可看出HHEDA的非劣解在總非劣解集S中所占比例最大,覆蓋區(qū)域最廣.

表9 算法比較結(jié)果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

具體到各算法的設(shè)計,ACO,CRO,HEDA,MMOEASA和MOSA-ACO側(cè)重于簡單采用算法自身進化機制執(zhí)行全局搜索,同時采用交叉、變異等鄰域操作或者較為單一的客戶交換操作執(zhí)行局部搜索.HHEDA進一步考慮問題特點及其解空間的復(fù)雜性,設(shè)計3種概率矩陣較全面有效地學(xué)習(xí)和積累優(yōu)質(zhì)解信息,并通過對其采樣生成新個體以實現(xiàn)全局搜索,有助于盡快發(fā)現(xiàn)優(yōu)質(zhì)解區(qū)域,同時在局部搜索過程中動態(tài)構(gòu)造豐富的變鄰域搜索(即啟發(fā)式算法HLS),有利于增加算法在優(yōu)質(zhì)解區(qū)域的搜索深度.因此,HHEDA對搜索方向的引導(dǎo)更為合理,能在相同時間內(nèi)更有效地發(fā)現(xiàn)高質(zhì)量的非劣解集.

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

5 結(jié)語

本文針對帶軟時間窗的同時取送貨車輛路徑問題,在考慮車輛載重等約束條件的同時,構(gòu)建了最小化輛行駛總里程和最大化服務(wù)準(zhǔn)時率的雙目標(biāo)模型,并提出一種超啟發(fā)式分布估計算法HHEDA 進行求解.HHEDA具有如下優(yōu)點:1)采用多種啟發(fā)式規(guī)則共同生成初始解,可確保初始種群的質(zhì)量和分散性;2)利用多個概率矩陣充分學(xué)習(xí)優(yōu)質(zhì)解中不同的有用信息,以推動全局搜索較快到達不同的優(yōu)質(zhì)解區(qū)域進行搜索;3)采用學(xué)習(xí)型超啟發(fā)式局部搜索執(zhí)行豐富的變鄰域搜索,能有效增加搜索的深度,從而提升算法的性能.仿真實驗和算法比較驗證了HHEDA是求解VRPSPDSTW的有效算法.后續(xù)研究將進一步考慮道路擁堵導(dǎo)致速度變化的問題,并設(shè)計有效求解算法.