改進(jìn)混合螢火蟲算法求解CVRP

白雪媛,張 磊,李 琳,武文喆

(1.沈陽航空航天大學(xué) 理學(xué)院,遼寧 沈陽 110136;2.沈陽航空航天大學(xué) 電子信息工程學(xué)院,遼寧 沈陽 110136)

0 引 言

有容量約束車輛路徑問題(CVRP)是NP-hard問題,其目標(biāo)函數(shù)為車輛服務(wù)完全部客戶后,所有車輛所走路程最短或所用費用最少[1]。目前,該問題求解方法有精確算法、啟發(fā)式算法和元啟發(fā)式算法。

因使用精確算法在可接受時間內(nèi)無法求得大規(guī)模CVRP最優(yōu)解,國內(nèi)外學(xué)者多采用啟發(fā)式或元啟發(fā)式算法進(jìn)行求解[2]。黃戈文等[3]利用整數(shù)編碼和先路由后分組實現(xiàn)解的轉(zhuǎn)換,提出自適應(yīng)遺傳灰狼優(yōu)化算法(AGGWOA)以研究CVRP。邵可南等[4]提出模擬退火-遺傳混合算法(GA-SA)解決帶硬時間窗和容量約束的單中心車輛路徑問題。Francesco等[5]用K-Means聚類提出CC-CVRP,將客戶點分為不同集群,將CVRP轉(zhuǎn)化為多個TSP。

螢火蟲算法(Firefly Algorithm,FA)是Yang[6]開發(fā)的啟發(fā)式算法。FA和不同方法結(jié)合,衍生出許多解決NP-Hard問題的算法[7]。標(biāo)準(zhǔn)螢火蟲算法因每只螢火蟲迭代尋優(yōu)時都會隨機(jī)移動,尋優(yōu)搜索過程復(fù)雜,算法收斂速度較慢,易陷入局部最優(yōu)。為提高算法性能,許多學(xué)者對FA做出多方面改進(jìn)。

Mohammadi等[8]用FA在考慮訂單的交貨時間和成本基礎(chǔ)上最小化總完成時間。李媛媛等[9]以K-Means聚類中心作為尋優(yōu)螢火蟲,結(jié)合鄰域與隨機(jī)吸引(N-R吸引)提出了KNRFA模型。閆軍等[10]將K-Means聚類和Q-leaning自啟發(fā)式引入蟻群算法,求解帶時間窗同時取送貨車輛路徑問題。

Li等[11]在FA中加入自適應(yīng)對數(shù)慣性權(quán)重,用步長減小因子動態(tài)調(diào)整隨機(jī)步長,提出LWFA。Mohsen等[12]為解決TSP將FA與局部搜索算法(2-opt和3-opt),遺傳算法中的變異交叉因子相結(jié)合,提出離散混合螢火蟲算法(HDFA)。Matthopoulos等[13]結(jié)合貪婪算法提出混合螢火蟲算法,解決異構(gòu)固定車隊VRP。Jaradat等[14]首先利用K-Means聚類將TSP分為子問題,再利用FA在每一類中尋找最優(yōu)路徑。

Altabeeb等[15-16]將FA與局部搜索(改進(jìn)的2-opt和2-h-opt)、變異交叉算子結(jié)合提出CVRP-FA,用田口方法統(tǒng)計確定參數(shù)的最佳值。但因CVRP-FA在2-h-opt處理后只接受有改進(jìn)方案,算法易陷入局部最優(yōu),無法產(chǎn)生更好的解,故將2-h-opt作為一種變異,而非局部搜索,提出CVRP-CHFA。該方式使算法可接受2-h-opt產(chǎn)生的所有解。該算法還引入合作島模型概念,保持解的多樣性,避免算法陷入局部最小值,達(dá)到加快算法收斂速度的目的。

該文提出一新的混合螢火蟲算法(KM-HFA)求解CVRP。引入有約束條件的K-Means聚類,以聚類結(jié)果為基礎(chǔ)構(gòu)建較好可行解作為螢火蟲尋優(yōu)初始解。結(jié)合部分匹配交叉和Intra-Swap,Inter-Swap變異算子,防止算法陷入局部最優(yōu)。引入2-opt和2H-opt局部搜索算法提高算法收斂速度,將算法用于求解標(biāo)準(zhǔn)算例。

對A,B,P,E組數(shù)據(jù)集72組算例,算法KM-HFA所得可行解比文獻(xiàn)[5]提出的CC-CVRP、文獻(xiàn)[13]提出的HFA結(jié)果好。當(dāng)客戶點逐漸增多時,KM-HFA的計算穩(wěn)定性顯著優(yōu)于HFA,最優(yōu)解與經(jīng)典解的百分比差異要比另兩種算法求解的百分比差異小得多,算例驗證了文中算法的有效性。

1 CVRP描述和數(shù)學(xué)模型

CVRP是經(jīng)典VRP,其定義在圖G(N,E)上,節(jié)點集為N={0,1,…,n}。節(jié)點0為倉庫(配送中心),節(jié)點{1,2,…,n}表示需服務(wù)客戶點,i∈N{0}的需求為qi,從節(jié)點i到j(luò)的運(yùn)輸成本為Cij。K為倉庫所擁有的車輛集合,Q為車輛最大載重。

CVRP數(shù)學(xué)模型描述如下:

(1)

s.t.

(2)

(3)

(4)

(5)

(6)

目標(biāo)函數(shù)(1)為最小化車輛總運(yùn)輸成本,約束(2)和(3)保證節(jié)點僅被訪問一次,約束(4)保證車輛運(yùn)輸過程中載重不超過車輛最大載重,約束(5)保證車輛從倉庫出發(fā)并最終回到倉庫。

2 改進(jìn)混合螢火蟲算法

FA作為元啟發(fā)式算法,隨機(jī)生成初始解,但尋優(yōu)結(jié)果在很大程度上依賴初始解,初始解較好時,其尋優(yōu)結(jié)果會更好。該文提出的KM-HFA先引入帶約束條件K-Means聚類,以聚類結(jié)果為基礎(chǔ),構(gòu)建較好可行解作為螢火蟲尋優(yōu)初始解,克服尋優(yōu)結(jié)果依賴初始解的問題。為增加螢火蟲種群的多樣性,將GA的部分匹配交叉和2H-opt交換算子與FA結(jié)合,進(jìn)行可行解間操作。最后,引入局部搜索算法2-opt和兩種變異算子提高算法收斂速度,防止其陷入局部最優(yōu)。KM-HFA的偽代碼如下,各步驟將在后續(xù)部分中具體說明。

KM-HFA算法

Begin

Input: 目標(biāo)函數(shù)f(S);

Output:全局最優(yōu)解S*;

Initialize:MI,C-R,M-R;

t=0; /*t代表迭代次數(shù)*/

初始種群P:S=(S1,S2,…,Sn);

Fori=1-P-Sdo

計算螢火蟲Si的目標(biāo)函數(shù)f(Si);

找出初始種群內(nèi)最優(yōu)解S*;

End for

計算漢明距離和亮度

whilet≤MI do

fori=1 toP-Sdo

forj=1 toido

IfIi>Ijthen

2H-opt procedureSi;

If rand

PMX procedureSi;

End if

2-opt procedureSi;

mutation procedureSi;

End if

End forj

End fori

排序,找出當(dāng)前最優(yōu)螢火蟲;

t=t+1;

End while

返回全局最優(yōu)解S*;

End

2.1 路徑編碼轉(zhuǎn)換

在計算過程中,可行解由K條(K為所用車輛數(shù)目)路徑組成,每條路都由n個點構(gòu)成,表示該路徑所分配到的客戶點個數(shù)。該文在計算螢火蟲亮度和可行解之間的操作時,需將多路徑CVRP轉(zhuǎn)換為單路徑TSP,進(jìn)行路徑編碼轉(zhuǎn)換(見圖1)。在完成可行解間操作后,編碼方式轉(zhuǎn)換回多路徑編碼,再對各可行解內(nèi)的路徑進(jìn)行局部搜索和變異操作。

圖1 路徑編碼轉(zhuǎn)換

2.2 有約束的K-Means聚類構(gòu)建初始解

標(biāo)準(zhǔn)FA的尋優(yōu)隨機(jī)生成初始解,尋優(yōu)結(jié)果對初始解有較強(qiáng)依賴性,當(dāng)初始解較好時,最終可找到更好解。該文先引入K-Means聚類,聚類后的客戶點間距離較接近,再進(jìn)行路線構(gòu)建。聚類使距離較接近的客戶點歸為一類,利用該特性進(jìn)行路線構(gòu)建得到的初始解比隨機(jī)生成的初始解好,將聚類后構(gòu)建的解為啟動螢火蟲尋優(yōu)算法的初始解。

2.2.1 K-Means聚類

K-Means聚類采用距離作為對象相似性的評價指標(biāo),兩個點離得越近就越相似,相似度和距離成反比。兩點間距離采用歐氏距離來計算:

(7)

每類有且只有一個聚類中心,該中心是類內(nèi)所有數(shù)據(jù)的平均值,計算點與聚類中心的距離,判斷出該點屬于哪一類。最終所有數(shù)據(jù)點被分為獨立的,由距離相近的點所組成的K類。K-Means聚類的K值需預(yù)先設(shè)定。為滿足車輛載重約束,K值由客戶點總需求和車輛最大載重計算,結(jié)果向上取整,即:

(8)

K-Means聚類的步驟為:

(1)計算兩點間歐氏距離;

(2)由公式8設(shè)定聚類數(shù)K,并隨機(jī)選取K個客戶點作為聚類中心Ci;

(3)將客戶點分配給距離其最近的聚類中心,直至所有客戶被分為K類;

(4)更新類的中心點,并重復(fù)步驟3直至簇不再發(fā)生變化。

2.2.2 構(gòu)建初始種群

初始種群由多個可行解構(gòu)成,每個可行解都按如下步驟進(jìn)行構(gòu)建:

(1)派出一輛車從倉庫出發(fā)隨機(jī)選取一類CLi進(jìn)行服務(wù);

(2)在該類內(nèi)隨機(jī)選擇一個客戶點并將該客戶點從列表刪除,以免重復(fù)服務(wù);

(3)將選中的客戶點加入到該路徑;

(4)判斷加入該客戶點后的車輛是否超載,若車上的載重未超過最大載貨量則重復(fù)步驟2和步驟3;否則將車輛返回倉庫,并派出下一輛車?yán)^續(xù)服務(wù);

(5)當(dāng)該類客戶都服務(wù)完,找到和該類中心最近的下一類中心Cj,對CLj類內(nèi)的客戶進(jìn)行服務(wù);

(6)重復(fù)操作至數(shù)據(jù)集內(nèi)客戶點都服務(wù)完。

上述步驟每進(jìn)行一次則得到一個可行解,多次重復(fù)后得到初始種群,即后續(xù)螢火蟲尋優(yōu)算法的初始解。對每個可行解計算目標(biāo)函數(shù),即車輛行駛的總距離。初始種群中,螢火蟲的亮度和目標(biāo)函數(shù)值成反比,目標(biāo)函數(shù)值越大亮度越小,目標(biāo)函數(shù)值較低的可行解就是更具吸引力的螢火蟲。根據(jù)亮度對初始種群內(nèi)的可行解進(jìn)行排序,目標(biāo)函數(shù)值最小的螢火蟲作為當(dāng)前的最優(yōu)解Sbest。

2.3 亮度計算

該文使用漢明距離(Hamming Distance,HD)計算尋優(yōu)過程中的螢火蟲亮度。HD為兩相等長度字符串對應(yīng)位元素不同的數(shù)量,即比較兩種可行解對應(yīng)位置的客戶是否相同,若客戶相同則兩者間的漢明距離不增加,若不同則漢明距離加1。

圖2展示了計算漢明距離的例子,Si和Sbest都有12個客戶點,在客戶點中,對應(yīng)位不同的客戶共7個,兩者之間的漢明距離為7。

圖2 漢明距離計算

設(shè)螢火蟲i的亮度為1到HD之間的隨機(jī)數(shù):

Ii=Random(1,HDi,best)

(9)

比較兩只螢火蟲的亮度,對亮度大的螢火蟲進(jìn)行后續(xù)操作。

2.4 可行解間的交換

2H-opt交換和部分匹配交叉都是針對整個可行解進(jìn)行處理的變換,涉及到路徑轉(zhuǎn)換,所以將2H-opt和部分匹配交叉放在一起進(jìn)行操作,可避免多次路徑編碼的轉(zhuǎn)換增加算法復(fù)雜度。在可行解間進(jìn)行2H-opt交換和部分匹配交叉可增加螢火蟲種群的多樣性,擴(kuò)大尋優(yōu)范圍,增加找到更優(yōu)解的可能性。

2.4.1 2H-opt交換

2.4.2 部分匹配交叉

交叉是遺傳算法中的操作算子,它可以挖掘種群的多樣性,克服啟發(fā)式算法容易陷入局部解的問題。該文將Si和Sbest作為父代,對兩者進(jìn)行部分匹配交叉(Partially Matching Crossover,PMX)操作,具體步驟通過圖3中的例子展示。

圖3 部分匹配交叉

在兩個父代中隨機(jī)選取兩者中相同位置的基因片段進(jìn)行交換,圖中灰色片段被交換;如圖3所示,由于交換后染色體中可能會出現(xiàn)如基因5,8,3等與交換片段所重復(fù)的基因,所以需要染色體進(jìn)行重復(fù)檢測并進(jìn)行去重處理;交換的片段建立映射關(guān)系,如基因5的映射點為基因3;重復(fù)的點將會由其映射點替換,注意此時替換的是未被交換的基因,被交換的基因不做任何處理;將重復(fù)的點全部替換完成后,最終形成兩個子代。交叉后得到兩子代,計算兩子代目標(biāo)函數(shù)值,與Si目標(biāo)值比較,若子代中有更小的解則更新Si,否則不更新。

2.5 可行解內(nèi)部交換

可行解內(nèi)變換利用局部搜索和變異算子對可行解內(nèi)路徑進(jìn)行處理,不僅可以提高算法的收斂速度,維持螢火蟲種群的多樣性,從而獲得更優(yōu)可行解,還能避免螢火蟲算法尋優(yōu)過程陷入局部最優(yōu)。

2.5.1 局部搜索

用局部搜索2-opt算子在種群內(nèi)產(chǎn)生新可行解。在可行解內(nèi)隨機(jī)選取某一條路徑,在該路徑中隨機(jī)選取兩個客戶點,將這兩點間的子片段進(jìn)行反轉(zhuǎn)操作后重新插入到該片段中,即可產(chǎn)生新的可行解。

2.5.2 變異算子

變異能維持螢火蟲種群的多樣性,一定程度克服FA易陷入局部最優(yōu)的缺陷。給定變異率,在區(qū)間(0,1)內(nèi)產(chǎn)生隨機(jī)數(shù),當(dāng)隨機(jī)數(shù)小于設(shè)定變異率,利用變異算子產(chǎn)生新的可行解。該文采用兩種變異算子:Intra-Swap是在可行解內(nèi)隨機(jī)選取一條路徑r1,隨機(jī)選取該路徑的兩個點c1,c2進(jìn)行交換得到新的可行解;Inter-Swap是在可行解內(nèi)選取兩條路徑r1,r2,在r1中隨機(jī)選擇一客戶c1,在r2中隨機(jī)選擇一客戶c2,將c1,c2進(jìn)行交換,得到新的可行解。

3 實驗及結(jié)果分析

文中算法用Python在Intel i5 1.8 GHz和4 GB內(nèi)存計算機(jī)上實現(xiàn)。測試實例包括以下在http://vrp.galgos.inf.puc-rio.br上公開提供的實例:客戶數(shù)量在100以下的小規(guī)模數(shù)據(jù)集A,B,P,E,共79組實例;客戶數(shù)量為101-1 001的中規(guī)模數(shù)據(jù)集X,選取其中15組實例。KM-HFA對每個實例獨立運(yùn)行10次,最大迭代次數(shù)為1 000。該算法參考了Altabeeb等[17]通過田口統(tǒng)計實驗所提出的參數(shù)值,種群數(shù)為20,變異率為0.1,交叉率為0.95。

表1列舉了文中所用數(shù)據(jù)的相關(guān)信息,其中Set為數(shù)據(jù)集名稱,Dimension為數(shù)據(jù)維度,Q表示該數(shù)據(jù)所用車輛的容量,Instance Number為從數(shù)據(jù)集中選取實驗實例的個數(shù),Time為該數(shù)據(jù)集創(chuàng)建時間。

表1 數(shù)據(jù)相關(guān)信息

文獻(xiàn)[5]利用K-Means聚類將客戶分為不同集群,提出CC-CVRP算法,把CVRP轉(zhuǎn)為TSP,對Set A, X進(jìn)行求解;文獻(xiàn)[13]結(jié)合貪婪算法提出HFA對Set A, B, P, E進(jìn)行求解。該文將KM-HFA同上述兩種算法進(jìn)行比較,比較各算法得出的最優(yōu)解所需的路徑距離(Best)及與網(wǎng)站公布的已知最優(yōu)解(B-K)之間的百分比差異(Gap/%)。

表2～表5為KM-HFA同CC-CVRP和HFA計算79組小規(guī)模實例的結(jié)果比較,表6為KM-HFA與CC-CVRP計算中規(guī)模數(shù)據(jù)集X中15組實例結(jié)果比較。其中,第一列為CVRP實例數(shù)據(jù),第二列為公布的已知最優(yōu)解,第三、四列為HFA最佳結(jié)果和它與最優(yōu)解的百分比差異,最后兩列為KM-HFA求得的最好解和它與最優(yōu)解的百分比差異。表中加粗部分為與其他算法相比該算法結(jié)果更好,“*”代表和經(jīng)典解有相同路徑數(shù)的更優(yōu)新解,“**”代表比經(jīng)典解增加一條路徑,車輛行駛總距離更小的新解。

表2 數(shù)據(jù)集A的KM-HFA與其他算法結(jié)果比較

由表2知,對A組27個實例,KM-HFA的結(jié)果26個優(yōu)于HFA和CC-CVRP的解,且A-n33-k5的解幾乎接近經(jīng)典解。值得注意的是,對實例A-n33-k6和A-n37-k6,在同路徑數(shù)條件下,KM-HFA找到了比經(jīng)典解更好的新可行解。路徑展示如圖4所示。

(a)A-n33-k6

對表3展示的B組23個實例,KM-HFA可行解有21個優(yōu)于HFA得到的可行解,且計算B-n50-k8的解幾乎接近經(jīng)典解。對B-n31-k5和B-n34-k5兩個實例,文中算法得到了與經(jīng)典解相同的結(jié)果。

表3 數(shù)據(jù)集B的KM-HFA與混合螢火蟲結(jié)果比較

KM-HFA和HFA對P組23個實例計算結(jié)果如表4所示,KM-HFA可行解有21個顯著優(yōu)于HFA的可行解。對實例P-n21-k2和P-n22-k2,兩種算法都搜尋到經(jīng)典解。KM-HFA算法在計算P-n16-k8,P-n19-k2,P-n20-k2時,在不增加路徑數(shù)的同時得到比經(jīng)典解更好的方案;在計算P-n22-k8,P-n23-k8兩組實例時,KM-HFA算法在比經(jīng)典解增加一條路徑的前提下,找到了車輛總行駛距離更小的解。

表4 數(shù)據(jù)集P的KM-HFA與混合螢火蟲結(jié)果比較

各實例對應(yīng)解為:P-n16-k8共8條路線:0-5-14-0; 0-1-0;0-4-11-0; 0-3-8-0; 0-10-12-15-0; 0-7-9-13-0; 0-2-0; 0-6-0; P-n19-k2共2條路線: 0-18-5-13-15-9-7-2-6-1-0;0-4-11-14-12-3-17-16-8-10-0; P-n20-k2共2條路線: 0-19-5-14-16-9-13-2-7-6-0;0-1-10-8-17-18-3-12-15-11-4-0; P-n22-k8共9條路線0-9-2-1-10-0;0-18-17-0;0-14-0;0-16-0;0-5-7-0;0-19-21-20-0;0-4-3-6-8-0;0-11-13-0;0-12-15-0; P-n23-k8共9條路線: 0-5-14-22-0;0-6-16-0;0-7-20-0;0-8-2-0;0-4-11-0;0-21-0;0-3-19-18-0;0-17-9-13-0;0-1-10-12-15-0。

表5呈現(xiàn)E組6個實例結(jié)果,其中KM-HFA的可行解有4個明顯優(yōu)于HFA的解。但在計算E-n33-k4實例時,HFA搜尋到經(jīng)典解,但是KM-HFA得到的最優(yōu)可行解也與經(jīng)典解非常接近。

表5 數(shù)據(jù)集E的KM-HFA與混合螢火蟲結(jié)果比較

表6是X中15個實例結(jié)果,其中KM-HFA的可行解有7個都優(yōu)于CC-CVRP的解;但CC-CVRP有8組實例結(jié)果更好;表明KM-HFA在計算中規(guī)模數(shù)據(jù)實例時,計算質(zhì)量有所降低。分析數(shù)據(jù)發(fā)現(xiàn),是因為中規(guī)模數(shù)據(jù)相比于小規(guī)模數(shù)據(jù),客戶數(shù)量較多、位置較分散。客戶點較分散時,K-Means聚類不能將客戶進(jìn)行更有效的分類,導(dǎo)致螢火蟲算法不能以較好的初始解開始尋優(yōu),算法計算性能下降。

表6 數(shù)據(jù)集X的KM-HFA與CC-CVRP結(jié)果比較

4 結(jié)束語

該文提出KM-HFA求解CVRP。通過引入K-Means聚類算法獲得更好初始解,將所得解作為螢火蟲尋優(yōu)初始解。結(jié)合遺傳算法中的PMX交叉在種群內(nèi)產(chǎn)生新的子代,以及兩種變異算子Intra-Swap、Inter-Swap,防止算法陷入局部最優(yōu)。為提升螢火蟲種群的多樣性,增強(qiáng)算法的局部搜索能力,加入2H-opt交換算子和2-opt局部搜索算子。

采用小規(guī)模數(shù)據(jù)79組,中規(guī)模數(shù)據(jù)15組,共94組數(shù)據(jù)檢驗KM-HFA的計算效果。實驗結(jié)果表明:在72組小規(guī)模實例中,KM-HFA得到的可行解比文獻(xiàn)[5]提出的CC-CVRP、文獻(xiàn)[13]提出的HFA得到的結(jié)果好。且當(dāng)客戶點增多時,KM-HFA計算穩(wěn)定性顯著優(yōu)于HFA的,得到的最優(yōu)解與經(jīng)典解的百分比差異要比另外兩種算法求得的解的百分比差異小得多。在15組中規(guī)模實例中,由于客戶數(shù)量較多、位置較分散,K-Means聚類不能將客戶進(jìn)行更有效的分類,導(dǎo)致螢火蟲算法不能以較好的初始解開始尋優(yōu),最終使KM-HFA的計算質(zhì)量有所降低。

KM-HFA在計算5組小規(guī)模實例A-n33-k6,A-n37-k6,P-n16-k8,P-n19-k2,P-n20-k2時,不增加路徑數(shù)的同時得到比經(jīng)典解更好的方案。在計算實例P-n22-k8和P-n23-k8時,比經(jīng)典解路徑數(shù)增加了一條的前提下,找到了車輛行駛總距離更小的解。

提出的KM-HFA還可應(yīng)用于解決其他VRP,例如VRPTW、綠色車輛路徑問題等。未來將研究如何把該算法與其它精確求解算法相結(jié)合,得到運(yùn)算效果更好的混合算法。