求解帶容量和時間窗約束車輛路徑問題的改進蝙蝠算法*

張 瑾,洪 莉,戴二壯

(河南大學計算機與信息工程學院，河南 開封 475004)

1 引言

車輛路徑問題VRP(Vehicle Routing Problem)是物流理論體系中一類經典的組合優化問題，隨著物流業的不斷發展，越來越多的VRP問題模型得以產生和發展，如帶容量約束的VRP CVRP(Capacitated VRP)、帶時間窗約束的VRP VRPTW(VRP with Time Windows)以及帶容量和時間窗約束的VRP CVRPTW(Capacitated VRP with Time Windows)[1]等。

由于上述問題都屬于NP難題，對于較大規模問題，精確算法難以在有限時間內給出最優解，因此普遍采用智能算法進行求解?，F有關于CVRPTW問題的算法大致有遺傳算法GA(Genetic Algorithm)、蟻群算法、禁忌搜索算法、粒子群優化PSO(Particle Swarm Optimization)算法、蝙蝠算法和狼群算法等，問題的優化目標大多為運輸距離最短或時間最短，部分文獻考慮了違反時間窗約束的懲罰成本。采用遺傳算法GA求解的文獻中，Jiang等[2]通過合理構造染色體，利用遺傳算法使最優解的搜索在可行空間內進行;趙振華等[3]提出一種先考慮客戶服務的先后次序，然后在此基礎上構建初始種群的方法，二者都以運距最短為優化目標；李加玲等[4]在遺傳算法中引入輪盤賭策略尋找可行的配送路徑；黃務蘭等[5]引入精英保留策略，采用客戶點交叉和路段交叉算子相結合的方式以保證算法中種群的多樣性，二者都以配送時間最短作為優化目標。采用蟻群算法求解的文獻中，李琳等[6]在算法的不同階段采用不同的信息素揮發策略防止算法陷入局部最優；辜勇等[7]以改進蟻群算法為主體，插入遺傳操作算子作為局部優化方法，二者都以運距最短為優化目標；黃震等[8]在節點選擇概率公式中引入時間窗因素來初始化種群，然后引入遺傳操作的交叉和變異算子對路徑進行優化；李奕穎等[9]采用改進的狀態轉移規則和輪盤賭選擇機制構建初始解，并結合k元素優化k-opt(k-optimization)鄰域搜索進行優化，二者都以配送車輛最少和運距最短為優化目標。采用禁忌搜索算法求解的文獻中，鐘石泉等[10]引入多初始解和全局禁忌表等措施，旨在擴大搜索范圍并減少解的不穩定性，以運距和違反時間窗約束的懲罰成本之和最小為優化目標；Chen等[11]在禁忌搜索算法中引入遺傳操作算子進行初始化，并采用2-opt操作進行局部優化；李明燏等[12]引入一條虛擬的路徑作為保留表，已選擇路徑上的用戶點可以與保留表中的用戶點進行交換和移動，以便擴大搜索空間，二者都以運距最短為優化目標。采用粒子群優化算法PSO求解的文獻中，李寧等[13]在初始化時將粒子群劃分為若干個兩兩重疊的相鄰子群，使算法跳出局部最優；馬炫等[14]提出一種基于粒子交換原理的整數粒子更新方法；羅耀[15]引入微生物行為機制中的趨化、繁殖和遷移算子對算法進行優化，三者都以總運距和違反時間窗約束的懲罰成本之和最小為優化目標；王飛[16]在慣性權重遞減的基礎上通過群體極值進行t分布變異，使算法跳出局部最優；Marinakis等[17]提出一種多自適應粒子群優化算法，二者都以最短運距為優化目標。采用蝙蝠算法求解的文獻中，馬祥麗等[18]針對帶時間窗車輛路徑問題的具體特征對蝙蝠算法的操作算子進行重新設計，以運距最短為優化目標；戚遠航等[19]以硬時間窗為約束條件(即車輛在客戶時間窗之外到達不進行配送，反之則為軟時間窗約束)，引入隨機插入搜索、最少客戶車輛插入搜索、普通插入搜索和交換搜索等策略來擴大搜索空間，以配送車輛最少和運距最短為優化目標；孫奇等[20]加入貪婪隨機自適應啟發式策略提高求解精度，引入病毒進化機制使蝙蝠算法跳出局部最優，以客戶滿意度最大和運距最短為優化目標。采用狼群算法求解的文獻中，葉勇等[21]利用近鄰初始化方式構建初始解，結合狼群算法覓食行為中的游走、召喚和圍攻3種行為重新定義其智能行為，以運距最短為優化目標。

上述文獻在求解CVRPTW問題時均取得了一定的成果，但其多考慮硬時間窗約束，目標函數大多以運輸距離或時間最短為優化目標，且采用的測試案例基本都是小規模的，未涉及較大規模問題的求解。由于在實際物流運輸中軟時間窗約束更符合客戶需求，本文將綜合考慮包含車輛使用成本、運輸成本和違反時間窗約束的懲罰成本在內的帶軟時間窗和容量約束的車輛路徑問題。蝙蝠算法是一種較為新穎的智能優化算法，在離散優化領域應用較少，為了提高CVRPTW問題的求解精度，本文設計了一種改進的離散蝙蝠算法DBA(Discrete Bat Algorithm)，并考慮了較大規模問題的有效求解。

2 求解問題描述

本文研究的問題為：在同時滿足客戶軟時間窗約束和車載容量約束的前提下，通過合理安排車輛配送路線，使得包含車輛使用成本、運輸成本和違反時間窗約束的懲罰成本最小。問題研究的前提條件如下所示：

(1)配送起點和終點唯一；

(2)所有車輛容量相同，每個客戶點僅由1輛車提供服務；

(3)每個客戶點的需求量已知，時間窗約束唯一，無優先服務約束。

模型中相關重要參數定義如下所示：

Z：總成本；

N：客戶點總量；

Q：車輛的最大載重量；

c：單位車輛的固定費用；

b：單位距離油耗成本；

dij：客戶點i到客戶點j之間的距離；

Ai：車輛到達客戶點i的時間點；

ti：為客戶點i進行服務所需時間；

qi：客戶點i的貨物重量；

Sj：客戶點j的時間窗開始點；

Ej：客戶點j的時間窗結束點；

e：早于客戶最早時間窗到達的懲罰系數；

l：晚于客戶最晚時間窗到達的懲罰系數；

tij：車輛從客戶點i到客戶點j的運行時間；

S：客戶點的編號集合，S?{1,2,…,N}；

x0jk：表示車輛k從配送中心駛向客戶點j；

xi0k：表示車輛k從客戶點i返回配送中心。

決策變量為：

K：所需車輛總數；

xijk：0-1變量，當車輛k由客戶點i到達客戶點j時為1，否則為0；

yik：0-1變量，當車輛k為客戶點i提供服務時為1，否則為0。

構建模型如式(1)～式(11)所示，其中式(1)為目標函數，表示包含車輛租賃成本、違反時間窗約束的懲罰成本以及與行駛距離相關的油耗成本在內的總配送成本最小；式(2)表示每輛車裝載的貨物總量不超過車輛的最大容量；式(3)表示每個客戶點由且僅由1輛車提供服務；式(4)保證客戶點j之前的臨近節點只有1個；式(5)保證客戶點i之后的臨近節點只有1個；式(6)表示從配送中心出發的車輛在完成配送任務后要返回配送中心；式(7)表示消除子回路即消除車輛不是從車場出發的現象，式(4)～式(7)共同保證了可行回路；式(8)和式(9)表示客戶時間窗約束；式(10)和式(11)表示決策變量的取值范圍。

(1)

(2)

(3)

?k∈{1,2,…,K}

(4)

(5)

(6)

2≤|S|≤N-1

(7)

Sj

(8)

Aj=Ai+ti+tij,i,j∈{0,1,2,…,N},i≠j

(9)

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

(10)

yki∈{0,1},i∈{1,2,…,N},?k∈{1,2,…,K}

(11)

3 模型求解的改進蝙蝠算法

蝙蝠算法BA是由Yang[22]提出的一種新型啟發式算法，算法利用蝙蝠通過回聲定位行為進行捕食的原理進行問題求解，具有參數少、穩定性高、求解速度快、尋優能力強等優點，已在工程優化、特征選擇、故障診斷和數據挖掘等多個領域展現出較好的應用效果[23 - 27]。蝙蝠回聲定位原理為：蝙蝠以脈沖的形式發射一定頻率和響度的超聲波，當超聲波在傳播的過程中遇到物體時會返回回聲，通過對接收到的回聲進行處理，蝙蝠可以檢測到物體相對于自身的距離和方向，以及物體的大小和運動速度，從而捕食或者避開障礙物。為便于模擬蝙蝠的回聲定位行為，Yang給出2個理想化規則：

(1)搜索規則：蝙蝠隨機飛行，同時以固定的頻率、可變的波長和音量的超聲波來搜索獵物。

(2)參數變化規則：蝙蝠根據自身與獵物的距離來自動調整脈沖波長和脈沖發射率，并限定聲音響度在指定范圍內依照給定方式由大到小變化。

算法求解過程由若干獨立搜索的蝙蝠完成，算法首先將每只蝙蝠視為當前可行域內的一個解，每個解對應一個由所優化問題確定的適應值，每只蝙蝠通過調整其脈沖頻率、聲音響度、脈沖發射率3項參數來追隨當前最優蝙蝠，使得整個種群在問題求解空間中產生從無序到有序的衍化，進而獲取最優解。脈沖頻率、速度和位置的更新方式如式 (12)～式(14)所示：

fa=fmin+(fmax-fmin)β

(12)

(13)

(14)

Begin

初始化算法相關參數及蝙蝠的位置、速度和脈沖頻率；

根據蝙蝠的初始位置計算適應度值，得出初始解，找出最優蝙蝠位置；

While(t<最大迭代數)

根據式(12)～式(14)分別調整蝙蝠的脈沖頻率、速度和位置；

If(當前隨機數>當前脈沖發射率)

在當前最優解附近根據式(15)進行局部搜索；

xnew=x*+ε*μ

(15)

其中，xnew表示隨機擾動得到的新解；ε表示服從標準正態分布的隨機向量，且ε～N(0,1)；μ表示常量，且0<μ<1；x*表示當前最優解。

Endif

If(當前隨機數<當前聲音響度,且xnew的適應度值優于x*)

接受新解，分別根據式(16)和式(17)更新脈沖發射率和聲音響度；

(16)

(17)

Endif

輸出全局最優解；

Endwhile

End

3.1 編碼策略

車輛路徑問題屬于離散的組合優化問題，首先對所求問題中的各個變量設計編碼策略，然后將其轉變為蝙蝠算法可以優化的變量。假設求解問題的規模為N，將客戶點編碼為從1到N的不同整數。

(1)解編碼策略。

針對CVRPTW問題的特性，解的編碼以待訪問客戶點編號序列表示，解中每個分量對應一個客戶點編號，因而得到解X的編碼形式為：X=(m1,m2,…,mN)。其中，N表示編碼長度，mi表示第i個要訪問的客戶點編號(mi∈[1,N]，且任意mi≠mj)。

(2)速度編碼策略。

vi表示蝙蝠訪問到第i個客戶點時的速度值，則速度的編碼形式為：V=(v1,v2,…,vN)。其中，V表示速度，編碼中每一個速度的值可以為正值或負值，vi∈[-(N-1),N-1]。

3.2 局部搜索策略

由于標準蝙蝠算法缺少擾動機制，存在算法后期收斂速度慢、易陷入局部最優等缺陷[28]，為此引入變步長搜索策略和2元素優化2-opt操作以增加解的多樣性，跳出局部最優。

(1)變步長搜索。

在搜索過程中，蝙蝠的移動步長應伴隨搜索過程的進行發生改變，在算法運行初期，步長應保持一個較大值，以避免算法過早陷入局部最優；隨著迭代次數的增加，步長應自適應地減小，最后保持一個較小值，使得算法后期加快收斂，得到更精確的值。標準蝙蝠算法中，蝙蝠位置的更新受ε和μ2個因子的影響，二者不能保證步長值遵循上述規則發生改變。若常量u取值較大，易使算法收斂速度變慢，若μ取值較小，則算法易陷入局部最優。

本文引入變步長搜索策略，以增加擾動機制和解的多樣性。定義步長擾動因子，其計算方法如式(18)所示：

τj=(Nitermax-j)/Nitermax

(18)

其中，τj表示在第j次迭代時的步長擾動因子，Nitermax表示最大迭代次數，j≤Nitermax。

改進的蝙蝠算法中，蝙蝠以變步長的方式進行局部搜索，搜索方式如式(19)所示:

xnew=x*+ε*μ*τj

(19)

(2) 2-opt優化操作。

2-opt操作即選定路徑上的任意2個節點并將節點間的路徑進行翻轉，得到新的路徑，以增加路徑搜索的多樣性，提高算法局部搜索能力。以7個客戶點的配送為例，假設客戶點為A、B、C、D、E、F、G，s為當前最優解，s={A,B,C,D,E,F,G}。以一定的概率進行2-opt操作，首先在s中隨機選擇不相鄰的2個節點B和E，然后將2個節點之間的路徑翻轉獲得新路徑，節點B之前的路徑保持不變添加到新路徑中，將節點B到節點E之間的路徑逆序排號后添加到新路徑中，節點E之后的路徑不變添加到新路徑中，則得到的新路徑s′={A,E,D,C,B,F,G}。

3.3 算法實現

為降低算法在搜索過程中的隨機性，縮小搜索范圍，加快搜索速度，本文首先引入聚類操作，對所有客戶點按其所在位置采用K-means算法進行聚類，使得每個客戶點最終歸屬于若干個不同分區。本文DBA算法的適應度值根據式(1)計算，具體求解方法如下：

Step1根據每個待配送客戶點所需配送的貨物重量，并結合車輛的載重限制把滿足要求的客戶點放入車輛運行路線中。

Step2如果該客戶點的貨物重量超出了車輛的限載量Q，則再申請1輛車，并將當前所需車輛總數加1；當遍歷完所有客戶點時，即可確定所需車輛總數K。

Step3每個客戶點都對應唯一的時間窗，若車輛不在客戶點相應的時間窗內到達，則給予一定的懲罰。

Step4計算油耗成本、違反時間窗的懲罰成本和車輛的租賃費用之和。

DBA算法實現步驟如下所示：

Step1對蝙蝠位置、脈沖頻率和速度進行初始化操作。

Step2將Step 1中所有初始蝙蝠個體的位置(即所有客戶點)按其所處位置利用K-means算法進行分區處理。

Step3蝙蝠位置更新，即對于任意蝙蝠a，根據式(12)～式(14)分別調整其脈沖頻率、速度和位置，從而產生后代蝙蝠。

Step4若當前隨機數小于當前脈沖頻率，則按Step 3進行全局搜索，并計算新的適應度值；否則在當前解附近根據式(19)以變步長搜索方式進行局部搜索，并以一定的概率進行2-opt操作，計算新的適應度值。

Step5若Step 4得到的適應度值小于當前最小適應度值，并且當前聲音響度大于當前隨機數，則利用Step 4得到的適應度值更新最小適應度值，并記錄各個蝙蝠的位置，得到全局最優解，否則不更新。

Step6根據式(16)和式(17)分別更新脈沖發射率及聲音響度。

Step7重復Step 4～Step 6，直到滿足結束條件，輸出全局最優解。

4 仿真實驗

由于遺傳算法、蟻群算法和粒子群優化算法是當前在CVRPTW問題求解領域應用較為廣泛的群智能優化算法，而蟻群算法在求解較大規模問題時速度相對較慢，因此，本文選取GA和PSO算法進行對比實驗。實驗測試環境為Windows 10 64位操作系統，CPU為3.4 GHz，內存為4.0 GB；仿真實驗平臺為Matlab R2017b。

Solomon數據集是目前帶時間窗車輛路徑問題最常用的標準測試庫，按照節點間的位置關系可以將Solomon數據集中的測試數據分為C、R和RC 3大類，每類問題又劃分為1和2 2個子類。3大類問題客戶點坐標及時間窗設置方式不同，其中C類數據中節點呈集簇式分布，節點分布于若干中心位置附近；R類數據中節點呈隨機分布，節點位置間無明顯集簇關系；RC類數據介于兩者之間，部分節點呈隨機分布，部分節點呈集簇式分布。子類1中問題的時間窗設置相對集中，子類2中問題的時間窗設置相對分散，且時間跨度較大。為了更全面地驗證本文算法求解效果，本文分別選取C1、C2、R1、R2、RC1、RC2問題中12個經典測試實例(見表1，其中，測試實例C101、C102、C104都屬于C1類問題。C101、C102、C104是Solomon測試集中的測試實例名稱，其問題規模均為101，不同測試實例之間的客戶點分布不同)，分別使用DBA、PSO和GA算法對每個實例進行10次計算，所有測試均采用全浮點數運算，算法所得最優解為10次計算所得最小值。

算法中各參數取值如下：DBA中脈沖頻率的最大值fmax=1，脈沖頻率的最小值fmin=0，聲音響度的衰減系數p=0.9，搜索頻率的增強系數γ=0.9，聲音響度A∈(0,1)，脈沖發射率r∈(0,1)，早于和晚于時間窗到達的懲罰系數e=l=0.35，單位距離油耗成本b=0.35，單位車輛的費用c=100，車輛的最大載重量Q=200；GA算法中交叉概率pc=0.3，變異概率pm=0.2；PSO算法中慣性權重因子w=0.2，加速系數C1=C2=2。

表1～表3給出了種群規模等于城市規模數，迭代次數分別為200,600和1 000時3種算法所得最優值及平均耗費時間。

Table 1 Optimal values and average time consumption of three algorithms when iteration number is 200表1 迭代次數為200時3種算法最優值和平均耗費時間

Table 2 Optimal values and average time consumption of three algorithms when iteration number is 600表2 迭代次數為600時3種算法最優值和平均耗費時間

從表1～表3可以看出，種群規模相同的情況下，迭代次數分別為200,600和1 000時，DBA算法對于12個實例所得最優值均在很大程度上優于GA和PSO算法的，GA算法獲得的最優值多數優于PSO算法的。在平均時間耗費方面，DBA算法稍長于GA和PSO算法。

為進一步驗證DBA算法性能，分別針對各問題計算GA或PSO算法最優值與DBA算法最優值之差與 DBA算法最優值的百分比，得到如下結果：表1中，C101問題改進效果最好，DBA算法相對GA和PSO算法分別改進44.71%和66.41%；對于RC201和R201問題，DBA算法相對GA和PSO算法改進效果較差，前者分別改進了5.73%和4.46%，后者分別改進了6.37% 和4.37%；對于表1中12個實例，DBA算法相對GA和PSO算法平均改進了17.28%和24.11%。表2中，C101問題改進效果最好，DBA算法相對GA和PSO算法分別改進了59.47%和79.38%；RC201問題改進效果最差，DBA算法相對GA和PSO算法分別改進了5.49%和4.12%；對于表2中12個實例，DBA算法相對GA和PSO算法平均改進了17.64%和22.61%。表3中，C101問題改進效果最好，DBA算法相對GA和PSO算法分別改進了59.38%和73.69%；對于C201和RC201問題，DBA算法相對GA和PSO算法改進效果較差，前者分別改進了3.71%和5.54%，后者分別改進了6.26%和3.43%；對于表3中12個實例，DBA算法相對GA和PSO算法平均改進了18.02%和21.24%?？傮w來看，3種問題中C類問題改進效果最好，R類問題其次，RC類問題最差；每種問題內1類問題相對于2類問題的改進效果更好。

Table 3 Optimal values and average time consumption of three algorithms when iteration number is 1 000表3 迭代次數為1 000時3種算法最優值和平均耗費時間

圖1～圖3分別給出了與表1～表3相同條件下3種算法分別對12個測試實例進行10次運算所得解的平均值隨問題規模變化的曲線圖。圖4給出了對于表1～表3中12個實例，DBA算法分別相對于GA和PSO算法的平均改進百分比條形圖。

Figure 1 Change curves of average values of three algorithms when population size equals to the number of cities and iteration number is 200圖1 種群規模等于城市數時迭代 200次3種算法平均值變化曲線

Figure 2 Change curves of average values of three algorithms when population size equals to the number of cities and iteration number is 600圖2 種群規模等于城市規模時迭代 600次3種算法平均值變化曲線

Figure 3 Change curves of average values of three algorithms when population size equals to the number of cities and iteration number is 1 000圖3 種群規模等于城市規模時迭代 1 000次3種算法平均值變化曲線

Figure 4 Average improvement percentage of DBA relative to GA and PSO for the 12 problems圖4 對于12個問題DBA分別相對于 GA和PSO的平均改進百分比

從圖1～圖3及前述數據分析結果可以發現，DBA算法求解所得最優值及平均值均優于GA和PSO算法的，GA算法相對于PSO算法獲得較好最優解的比例更高，且在迭代1 000次時的DBA算法相對于GA算法改進效果最好。

從圖4可以看出，DBA算法相對于GA算法的平均改進百分比最大為18.02%，最小為17.28%；DBA算法相對于PSO算法的平均改進百分比最大為24.11%，最小為21.24%。上述圖形及數據表明，本文所設計的DBA算法相對GA和PSO算法具有較好的求解效果。

表4給出了迭代次數為1 000時，種群規模為30的DBA算法、種群規模等于城市規模數的GA和PSO算法所得平均值、最優值和平均耗費時間。

從表4可以看出，迭代次數相同的情況下，種群規模為30的DBA算法所得最優值、平均值和平均時間耗費均優于采用較大規模種群的GA和PSO算法。在上述12個測試實例中，種群規模為30的DBA算法的平均求解時間相對于同等規模的GA和PSO算法分別加快了158.16%和168.98%；GA和PSO算法的最優值相對于DBA算法的平均值的平均改進百分比分別為6.73%和9.27%，最好改進百分比分別為19.06%和33.97%。數據分析結果表明，DBA算法在問題求解精度和求解效率方面均有較為明顯的優勢。

5 結束語

帶時間窗和容量約束的車輛路徑問題在現實中有著廣泛的應用，對該問題進行優化的研究從未停歇。本文研究了求解該問題的離散蝙蝠算法，利用K-means算法對客戶點按其所在位置進行聚類，在DBA算法中引入了變步長搜索策略和2-opt操作進行局部搜索。

實驗結果表明：所設計的離散蝙蝠算法具有較快的收斂速度和較強的尋優能力，能夠有效降低配送成本。本文設計的算法只是對標準蝙蝠算法的初步改進，對測試庫內子類2中問題改進效果相對較差，后續研究中將考慮對子類2問題的時間窗進行聚類處理，并將蝙蝠算法與其它智能算法思想相結合，以進一步提高算法的求解性能。

Table 4 Optimal values and average time consumption of three algorithms with different population sizes when iteration number is 1 000表4 迭代次數為1 000時3種算法在不同種群規模下的最優值和平均耗費時間