蔣海青，趙燕偉，冷龍龍

(1.浙江工業(yè)大學 特種裝備制造與先進加工技術(shù)教育部重點實驗室，浙江 杭州 310014; 2.中國計量大學 現(xiàn)代科技學院，浙江 杭州 310014)

0 引言

半個多世紀的研究使車輛路徑問題(Vehicle Routing Problem, VRP)從問題描述到求解方法都取得了豐碩的成果，但如何獲得滿意的求解方法依然是VRP中具有挑戰(zhàn)性的研究領(lǐng)域[1]。VRP的求解方法大體分為精確式、啟發(fā)式和智能優(yōu)化3種，精確求解方法如動態(tài)規(guī)劃、分支定界等[2]，啟發(fā)式算法如節(jié)約法[3]、掃描法[4]及文獻[5-6]對這兩種算法的改進。由于啟發(fā)式算法一般求到的是局部最優(yōu)解，VRP目前的求解方法主要以智能優(yōu)化算法為主，如遺傳算法(Genetic Algorithm, GA)[7-8]、禁忌搜索(Tabu Search, TS)算法[9-10]、模擬退火(Simulated Annealing, SA)算法[11-13]、蟻群優(yōu)化(Ant Colony Optimization, ACO)算法[14-17]。

石兆[18]總結(jié)智能優(yōu)化算法在VRP的應用認為，禁忌TS算法搜索時間短，但禁忌規(guī)則難以確定;SA算法在不考慮計算時間時求得的解最好;ACO算法的搜索速度較慢，且較易出現(xiàn)早熟問題。鑒于智能優(yōu)化算法的不足，近年來人們嘗試從以下兩方面進行算法改進：

(1)將多種算法進行集成 如Boon等[19]提出的集成差分算法，通過擴展TS算法求解車輛路徑問題；Teymourian等[20]將改進智能水流算法(Enhanced Intelligent Water Drops, EIWD)和布谷鳥搜索(Cuckoo Search, CS)算法結(jié)合，并應用局部搜索混合算法(Local Search Hybrid Algorithm, LSHA)和后優(yōu)化混合算法(Post-Optimization Hybrid Algorithm, POHA)局部搜索解決帶有容量約束的車輛路徑問題(Capacity Vehicle Routing Problem, CVRP)；Akpinar[21]將大規(guī)模鄰域搜索(Large Neighbourhood Search, LNS)與蟻群算法ACO結(jié)合組成LNS-ACO算法求解CVRP問題，該算法應用ACO構(gòu)造新的解，應用LNS進行解的改進，結(jié)合基于插入的局部優(yōu)化方法在大量CVRP案例中獲得了滿意的結(jié)果；張景玲等[22]將量子進化算法與免疫算法結(jié)合，提出自適應免疫量子進化算法求取動態(tài)VRP；曹高立等[23]設(shè)計了混合量子算法求解CVRP，并采用二維量子位觀測模型和可見度生成解改善求解效果。

(2)將新算法引入VRP的求解 例如Szeto等[24]應用人工蜂群(Artificial Bee Colony, ABC)算法求解CVRP,并結(jié)合insertion，swap多種局部優(yōu)化方法改善解；Korayem等[25]提出運用灰狼算法求解CVRP；Tahayassene等[26]提出一種和聲搜索算法meta-HSA(meta-harmony search algorithm)，通過設(shè)置求解器和優(yōu)化器求解帶時間窗的車輛路徑問題(Vehicle Routing Problems with Time Windows, VRPTW)，并采用交換、2-opt的方法進行局部優(yōu)化搜索,以提高解的性能；Yurtkuran等[27]提出類電磁場算法(Electromagnetic Field Algorithms, EFA)求解CVRP，并應用swap局部優(yōu)化方法改善解的性能;Dorronsoro等[28]提出元胞遺傳算法(Cellular Genetic Algorithm, CGA)求解CVRP，并構(gòu)建基于交換和2-opt的局部優(yōu)化方法提高解的性能。

上述改進算法雖然在解決VRP的有效性和優(yōu)越性都得到了驗證，但還是存在一些問題：①在子代解的構(gòu)造中，當父代和母代解不合法時，就不能由其再產(chǎn)生子代解，從而限制了全局搜索范圍；②兩種算法集成時一般是由一種算法生成一個解，另外一種算法進行解的改善，較難實現(xiàn)兩種算法的徹底融合。

化學反應優(yōu)化(Chemical Reaction Optimization, CRO)算法是2010年由Lam等[29]提出的用于解決NP-hard問題的新型算法，算法根據(jù)分子在其勢能面(Potential Energy Surface, PES)中能量越低狀態(tài)越穩(wěn)定的特性，通過模擬化學反應中分子間的相互作用來尋求分子最小勢能。在求解時，通過其基本組成單元——分子結(jié)構(gòu)來描述問題特性，即用分子表示問題的解，通過分子內(nèi)部原子之間的反應(如碰撞、交換等)改變分子結(jié)構(gòu)，尋求分子勢能降低的途徑。隨著解的改進，分子勢能逐步下降，當其勢能最低時，分子達到穩(wěn)定狀態(tài)，這時的分子結(jié)構(gòu)即為問題的最優(yōu)解。與ACO,SA,GA等智能算法的不同之處是,CRO的參數(shù)意義簡單且易于設(shè)置，只需要設(shè)置一個參數(shù)就能發(fā)生化學反應，每次反應都可以獲得新的混合解，使其每代都有廣泛的搜索空間。除此之外，其4種化學反應極易與其他智能算法結(jié)合，實現(xiàn)算法集成。因此，CRO已經(jīng)被成功應用于連續(xù)性問題求解[30]，如電流優(yōu)化問題[31]、網(wǎng)絡(luò)編碼優(yōu)化[32]、自行車定位問題[33]、集合覆蓋問題(Set CoveringProblem, SCP)[34]、背包問題[35]、任務調(diào)度[36]、頻譜分配問題[37]等各類組合優(yōu)化問題。盡管如此，文獻[38]在對CRO的分析中認為其適合于目標問題的求解，而VRP與背包問題、分配問題及定位問題一樣屬于離散型NP-hard問題，且具有諸多相似性，因此本文探討應用CRO算法解決CVRP，以期豐富CVRP的有效求解方法。為適應CVRP的特征，本文設(shè)計了包括配送中心、顧客整數(shù)序列的分子結(jié)構(gòu)，采用節(jié)約法進行了分子的初始化，以提高初始解的質(zhì)量；在分子的交換、撞墻反應中借鑒遺傳算法的交叉、變異操作，應用分解反應擴大求解范圍，將單個分子分解成保留父代、母代一定特征的兩個分子，在合成反應中充分考慮了距離對成本的影響，將兩點間的最短距離作為子代分子的選擇依據(jù)；通過設(shè)置動態(tài)變化參數(shù)控制4種化學反應發(fā)生的頻率，來提高算法的執(zhí)行效率。為測試算法的有效性，本文對經(jīng)典的benchmark基準實例的P，A問題進行了求解，并與當前先進算法進行了比較，實驗結(jié)果顯示，所提CRO算法具有一定的優(yōu)勢。

1 問題描述

CVRP可以描述為：由一個配送中心為多個顧客服務，車輛從配送中心出發(fā)，配送完成后返回出發(fā)的配送中心，問題的目標是獲得總距離最小的最佳配送方案。將顧客與配送中心組成的網(wǎng)絡(luò)圖設(shè)為G=(N,A),其中：N=Nm∪Nc,Nm=0為配送中心集合，Nc=(1,2,…,c)為顧客點集合；A={(i,j):i,j∈N,i≠j}為點i,j的弧集合。K=(1,2,…,k)為車輛集合，dij為i，j的距離，車輛的載重限制為CQ,qi為i點需求。為便于建模，本文對問題進行如下假設(shè)：①車輛勻速行駛，路上無擁堵現(xiàn)象，路況平穩(wěn)；②配送車輛數(shù)量不限，車輛類型只有1種；③客戶需求已知并確定，每個客戶必須且只能被服務一次。用模型表達如下：

(1)

s.t.

(2)

(3)

(4)

(5)

(6)

xijk∈{0,1},?i,j∈N,k∈K。

(7)

其中：式(1)表示距離最短的目標函數(shù)；式(2)和式(3)保證每個顧客只能被服務一次；式(4)表示顧客需求限制；式(5)表示消除子回路；式(6)保證每輛車配送完后回配送中心；式(7)表示變量的取值。

2 化學反應優(yōu)化算法

設(shè)分子ω的能量由勢能PEω和動能KEω組成，在反應過程中，PEω逐漸減小， CVRP的目標是獲得最短距離，因此本算法將PEω作為目標函數(shù)。但如果在算法設(shè)計時，單純追求勢能最小值，則極有可能只能求取到局部最優(yōu)解。為了求得全局最優(yōu)解，在CRO算法中通過調(diào)節(jié)動能，使新產(chǎn)生的分子ω′的PEω′即使大于PEω，ω′也能以一定的概率保留。CRO算法尋優(yōu)的基本反應包括單分子反應和雙分子反應，單分子反應主要有撞墻反應和分解反應，雙分子反應主要有交換反應和合成反應。當分子之間的碰撞強度不大時，發(fā)生撞墻反應或交換反應，這兩個反應一般獲得原分子的鄰域解，分解反應和合成反應在當分子發(fā)生劇烈碰撞時用于打破局域解。下面簡單介紹這些反應：

(1)撞墻反應 指分子與容器壁發(fā)生碰撞后形成新分子的過程。若原分子為ω,則新分子ω′=neibour(ω)，即新分子為原分子的鄰域函數(shù)，撞墻反應發(fā)生的條件為

PEω+KEω≥PEω′。

(8)

反應后，KEω′=(PEω+KEω-PEω′)×q，其中q∈[KELossRate,1]，KELossRate為能量損失最大百分比，損失的能量值保存在能量中心，用于提供分解反應所需能量，將能量中心的能量值表示為buffer。

(2)分解反應 指分子在容器中發(fā)生猛烈碰撞后，由1個分子分解為2個分子的過程。若原分子為ω，新分子為ω1和ω2，則分解反應發(fā)生的條件為

PEω+KEω≥PEω1+PEω2。

(9)

因在撞墻反應中，KEω逐漸減小，當式(9)不易滿足時，若能滿足

PEω+KEω+buffer≥PEω1+PEω2，

(10)

則能發(fā)生分解反應。

若式(9)成立，則KEω 1=temp×k，KEω 2=temp×(1-k)，其中temp=PEω+KEω-(PEω1+PEω2)，k為[0,1]的隨機數(shù)；若式(10)成立，則KEω 1=(temp+buffer)×m1×m2，KEω2=(temp+buffer-KEω1)×m3×m4，其中m1，m2，m3，m4為[0,1]的隨機數(shù)。

(3)交換反應 指容器中兩個分子發(fā)生碰撞，各自被彈開后產(chǎn)生的變化。若原分子為ω1和ω2，新分子為ω′1和ω′2，則交換反應發(fā)生的條件為

KEω1+KEω2+PEω1+PEω2≥PEω′1+PEω′2。

(11)

反應后，temp=KEω1+KEω2+PEω1+PEω2-(PEω′1+PEω′2)，KEω′1=temp×k，KEω′2=temp×(1-k)，其中k為[0,1]的隨機數(shù)。

(4)合成反應 指容器中兩個分子發(fā)生劇烈碰撞后合成一個分子的過程。若原分子為ω1和ω2，新分子為ω′，則合成反應發(fā)生的條件為

KEω1+KEω2+PEω1+PEω2≥PEω′。

(12)

反應后，KEω′=KEω1+KEω2+PEω1+PEω2-PEω′。

3 CRO應用于CVRP問題

3.1 初始分子

將分子結(jié)構(gòu)表示為由配送中心和顧客組成的整數(shù)序列，對于有n個顧客的配送，將配送中心編號設(shè)置為1，顧客編號從2～n開始。首先對每個顧客i形成1-i-1的線路，按Sij=Ci1+Cj1-Cij計算插入顧客j后的Sij值,并將其由大到小排序，若最大Sij值對應的j顧客插入能滿足載重條件，則將該j顧客插入1-i-1路線，直到所有顧客都能安排為止。如圖1所示為構(gòu)建的某初始分子，表示2輛車服務8個顧客，第1輛車服務2，3，7顧客，第2輛車服務1，4，5，6，8顧客，將此時的距離設(shè)置為初始勢能值。

3.2 撞墻反應

借鑒GA的互換變異操作，撞墻反應采用位置交換法改變分子結(jié)構(gòu)的流程如下：

(1)在容器庫中隨機選擇一個分子ω。

(2)隨機產(chǎn)生兩個位置，交換兩個位置的顧客，產(chǎn)生新分子ω′，其執(zhí)行流程如圖2所示。

3.3 分解反應

CRO算法分解形成的兩條新路線不僅應保留原路線的部分特點，還要能擴大搜索范圍。借鑒GA的循環(huán)交叉設(shè)計以下分解反應：從種群中隨機選一個分子ω，再隨機選一個位置pos，分解后的分子ω1保留pos前部分，將后部分反轉(zhuǎn)，ω2保留pos后部分，將前部分反轉(zhuǎn)，若原分子為34825679，pos=4,則分解后的兩個分子分別為34829765和28435679，分解反應的執(zhí)行流程如圖3所示。

3.4 合成反應

合成反應是將兩條路線合并生成一條路線的過程?？紤]到問題的目標是總路徑最短，在合成過程中參考文獻[39]，將相鄰點距離最短的顧客點作為新路徑節(jié)點。具體實施過程為：首先在種群中隨機選擇兩個不相同的分子A1和A2；在分子A1中，隨機選擇一個顧客作為合成路徑的第1個顧客點，分別計算A1，A2中該顧客點與其后相鄰點的距離，選擇最短距離所對應的顧客點為合成路徑的下1個顧客，重復上述過程直到合成后的分子包含所有顧客。如圖4所示，A1中pos=4，若A22～9點的距離大于A1中2～5點的距離，則將顧客5放入新分子的位置2。合成反應的具體流程如圖5所示。

3.5 碰撞反應

參考文獻[40]的交叉方法隨機選擇兩個不同的分子A和B，隨機選擇分子A的顧客點ai和分子B的顧客點bj(i=j)，交換ai與bj后的顧客點，產(chǎn)生兩個新分子，對新分子中重復的顧客點進行映射。當選擇交換位置為4時，若交換前A，B分子分別為34825679和67534829，則交換后的分子A1和B1分別為34827569和67538249，碰撞反應的具體程序如圖6所示。

3.6 算法流程

CRO求解CVRP問題的算法流程如圖7所示。首先生成初始分子并計算其勢能值，確定初始最小值，然后隨機產(chǎn)生1個0～1均勻分布的t值，如果t大于預先設(shè)定的多分子反應概率MoleColl值，則執(zhí)行雙分子反應，否則執(zhí)行單分子反應。對于單分子，將HitNumber與α進行比較，然后選擇分解反應或撞墻反應；對于雙分子，通過比較KE與β來選擇相應的反應，具體執(zhí)行過程如圖8所示。

4 數(shù)據(jù)實驗及分析

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

CRO算法的主要參數(shù)有：初始分子的個數(shù)Popsize，影響單雙分子反應選擇的MoleColl,初始動能InitialKE，合成、交換兩種雙分子反應控制參數(shù)β1，分解、撞墻單分子反應控制參數(shù)α1，最大迭代次數(shù)num。當MoleColl<0.5時，雙分子反應權(quán)重增大，反之單分子反應權(quán)重較大。圖9所示為對P-n19-k2案例MoleColl取值的10次模擬平均結(jié)果，由圖可知，該案例中MoleColl=0.3，0.8時取得了相同的目標函數(shù)，即單雙分子反應權(quán)重對解的影響沒有明顯差異，因此本文取MoleColl=0.5。

InitialKE是保證反應過程中能量守恒的參數(shù)，其設(shè)置過小易導致局域解[33];β控制合成、交換雙分子反應的概率，其值越小，合成反應執(zhí)行的幾率越小，影響全局搜索，但如果過大，則碰撞反應的幾率太?。沪林涤绊懛纸?、撞墻反應的發(fā)生概率，若α設(shè)置過大，則影響分解反應的發(fā)生概率，降低全局搜索，但過小又會影響碰撞反應的發(fā)生。本文按式(13)和式(14)設(shè)定α，β的值:

α=i/N×α1;

(13)

β=N/i×β1。

(14)

式中：N為總循環(huán)次數(shù)，i為當前循環(huán)次數(shù)，α1和β1為常數(shù)。

參考文獻[33]的方法，設(shè)置初始參數(shù)α1=0，β1=0，InitialKE=179，Popsize=50，然后依次改變參數(shù)InitialKE，β1，Popsize，α1的取值。以P-n19-k2為測試案例，圖10所示為參數(shù)取值與目標函數(shù)和CPU用時的關(guān)系。由圖10a和圖10b可知，α1，β1設(shè)置太大或太小均不能獲得最佳解，說明幾種反應的權(quán)重系數(shù)對求解結(jié)果的影響具有隨機性，當α1=40，β1=80時解最小，其CPU用時的變化極小(不超過4 s)；由圖10c確定InitialKE=716；由圖10 c確定InitialKE=716;由圖10d種群與解及CPU的關(guān)系確定Popsize=150。

迭代次數(shù)越多，分子越能充分地在解空間進行搜索，分子間的信息交換越充分，越能搜索到更好的解。圖11及表1所示為迭代次數(shù)與解的關(guān)系，表1第2列是P-n19-k2取上述參數(shù)隨機運行10次的平均值，第3列是其平均CPU用時，圖11所示為與表1對應的圖。由圖11可知，隨著迭代次數(shù)的增大，解得到了改善，但當其增大到一定程度后(3 000代后)，對優(yōu)化結(jié)果的提升效果趨小，但其CPU用時卻大幅提高。因此，本文在優(yōu)化結(jié)果和計算時間之間尋求一個平衡，將最大代數(shù)設(shè)為3 000。

結(jié)合圖10和圖11的結(jié)果及后期大量的實驗，本文參數(shù)設(shè)置如下：MoleColl=0.5，InitialKE=716，Popsize=150，α1=40，β1=80，num=3 000。

4.2 結(jié)果分析與比較

本次計算采用MATLAB 7.0編程,計算機的配置為Pentium(R)Dual-Core E5300 2.60 GHz CPU，8 GB RAM，每個案例計算10次，以最大循環(huán)次數(shù)為結(jié)束條件。為了測試算法的效果，將本文設(shè)計的CRO算法與標準庫案例結(jié)果進行對比，標準庫案例結(jié)果來自網(wǎng)址為http://www.bernabe.dorronsoro.es/vrp/的網(wǎng)站，數(shù)據(jù)采集的時間為2017年2月10日。表2所示為CRO算法的計算結(jié)果和標準庫最優(yōu)解的比較。其中：BKR為該標準庫的最優(yōu)解，Best為本算法求得的最佳解，Avg為平均值，第五列表示平均CPU用時，Dev為偏差。偏差的計算公式[25]如下：

由表2可知，在對標準庫9個代表性的P，A案例求解中，CRO有4個案例取得了與標準庫相同的最優(yōu)解，占總案例(9個)45%，其他未獲得最優(yōu)解案例的偏差值絕大多數(shù)不超過5%。為了進一步驗證CRO算法的求解效果，將其與文獻[21]3種策略的大規(guī)模領(lǐng)域搜索(LNS)算法進行比較，比較結(jié)果如表3所示，表中：fbLNS，fwLNS分別表示不同策略LNS算法的最好解值和最差解值，fCRO表示本算法的最好解值；bRPD，wRPD表示本算法與LNS的最好解和最差解的相對差異百分比：

(16)

(17)

兩種算法的CPU用時比較如表4所示。

表3 CRO與LNS結(jié)果比較

由表3可知，CRO算法在絕大多數(shù)案例中都取得了與LNS相同的最優(yōu)解，部分結(jié)果優(yōu)于LNSi，沒有結(jié)果優(yōu)于LNSa和LNS-ACO，但平均解偏差與LNSa和LNS-ACOCRO相比均不超過1%，與LNSi的平均結(jié)果相近，說明CRO算法具有較強的搜索能力。由表4可知，相比LNS算法，CRO算法的CPU用時在多個算例中減少近30%，因此CRO算法更適合于實時調(diào)度，將CRO與基于擴展節(jié)約法的集合覆蓋算法(SetCovering Based Extended Savings Algorithm, SC-ESA)、粒子群優(yōu)化(Particle Swarm Optimization, PSO)算法、帶局部搜索的差分進化算法(Differential Evolution algorithm with Local Search, DELS)[19]幾種常用算法進行對比，比較結(jié)果如表5所示，表中RPD為相對差異百分比，

(18)

式中：fCRO為本算法的最優(yōu)解，fmin為其他算法的最優(yōu)解。

表5 CRO與其他算法結(jié)果比較

由表5可知，CRO算法的9個案例中有2個案例的求解結(jié)果優(yōu)于SC-ESA，總體結(jié)果與SC-ESA相近，略差于經(jīng)典的PSO，DELS算法，但也有近50%的案例結(jié)果與PSO，DELS算法相同，平均解偏差不超過1.5%，因此證明CRO算法具有較強的搜索能力。

表6 CRO算子的實驗結(jié)果

圖12和圖13所示分別為應用CRO算法得到的表2中案例2的線路圖和收斂圖。由圖13的收斂圖可知，算法初始階段，由于α較小、β值較大，合成、分解反應發(fā)生的概率較高，目標函數(shù)值減小的幅度較大；隨著迭代次數(shù)的增加，α逐步增大，β值逐步減小，合成、分解反應發(fā)生的概率降低，目標函數(shù)值降低的幅度變小。

5 結(jié)束語

本文將CRO算法應用于CVRP，主要進行了以下研究：

(1)將CRO算法的分子設(shè)計為反映顧客與配送中心序列的實數(shù)碼，借鑒GA的交叉、變異和局部優(yōu)化方法對CRO算法的撞墻、分解、合并和交換4個反應進行設(shè)計。具體為在撞墻反應中采用互換變異操作，分解反應中采用循環(huán)交叉操作，將單條路徑分解為帶有母體部分特征的兩條路徑，以擴大解的搜索范圍；交換反應采用次序交叉操作，借鑒局部優(yōu)化思想設(shè)計合成反應，在合成反應過程中獲得局部路徑最小解。

(2)采用動態(tài)方式設(shè)置CRO算法的兩個重要參數(shù)α和β，將全局搜索與局部搜索有效結(jié)合，獲得了較強的搜索能力。將求取的解與標準benchmark的6個P類和3個A類案例進行對比，結(jié)果顯示CRO算法能夠很好地求解CVRP，絕大多數(shù)結(jié)果與標準庫的數(shù)據(jù)誤差不超過5%，但其CPU用時減少幅度較大。隨著實時配送對顧客滿意度影響的增大，在不增加或少增加配送成本的基礎(chǔ)上，增大配送實時性有助于提高顧客滿意度，由此證明CRO算法是一種求解CVRP的有效算法。

下一步主要研究將CRO算法應用于帶時間窗的車輛路徑問題、動態(tài)車輛路徑問題(Dynamic Vehicle Routing Problem, DVRP)等其他類型的VRP中，進一步拓展CRO算法的應用范圍。