結合聚類分解的增強蟻群算法求解復雜綠色車輛路徑問題

胡蓉 李洋 錢斌 金懷平 向鳳紅

傳統的車輛路徑問題(vehicle routing problem,VRP)由Dantzig 和Ramser 于1959 年首次提出[1].該問題主要描述為在滿足車輛載重、容積、行駛里程及客戶服務要求等條件的同時,合理調度車輛出行數量、行車路線、出行時間,使得總運費用最優化.隨著社會經濟的高速發展,跨區域多車隊聯合車輛配送需求明顯增加[2?4].同時,當今環保要求越來越嚴格,考慮燃油消耗和碳排放等因素的低能耗車輛配送已開始受到重視[5].此外,日益激烈的市場競爭逼迫企業注重降低物流配送成本和提高客戶滿意度[6].在此背景下,研究帶時間窗的低能耗多車場多車型車輛路徑問題(Low-energy-consumption multi-depots heterogeneous-fleet vehicle routing problem with time windows,LMHFVPR_TW),

具有十分重要的現實意義.由于VRP 為NP-hard問題,而VRP 可歸約為LMHFVPR_TW,故LMHFVPR_TW 也屬于NP-hard 問題,對其開展研究亦具有較大理論價值.

低能耗車輛路徑問題近年已開始受到重視,但相關研究仍較為有限.在低能耗多車場車輛路徑問題 (Low-energy-consumption multi-depots VRP,LMVRP)方面,Jabir等[7]在問題中考慮碳排放等因素,利用蟻群算法與變鄰域搜索(Variable neighborhood search,VNS)相結合進行求解.仿真實驗表明,該算法可有效求解小規模和大規模問題.Kaabachi等[8]在問題中考慮距離、碳排放和燃油消耗等因素,并在蟻群算法中加入插入、交換和2-opt等鄰域搜索策略以提高算法局部搜索能力,取得良好的效果.在低能耗多車型車輛路徑問題 (Low-energy-consumption heterogeneous-fleet VRP,LHFVRP)方面,Xiao等[9]在問題中考慮交通擁堵和碳排放等因素,并設計融合變鄰域搜索和整數規劃的迭代搜索算法進行求解.Kwon等[10]在問題中考慮碳排放交易機制(即各國之間通過碳排放交易市場相互買賣碳排放量,以保證該國的碳排放量在其規定的配額內),然后采用混合禁忌搜索算法進行求解.目前尚無LMHFVPR_TW 的相關研究.

對于多車場車輛路徑問題 (Multi-depot VRP,MVRP)、多車型車輛路徑問題(Heterogeneousfleet VRP,HFVRP)和多車場多車型車輛路徑問題(Multi-depot heterogeneous-fleet VRP,MHFVRP)這幾類復雜VRP,智能算法已有一定的研究.現有研究大多對問題進行整體編碼和求解,但由于這些問題比單車場單車型VRP (即傳統VRP)具有更多的決策變量和約束條件,使得解空間擴大且編碼解碼較為繁瑣,這時僅靠智能算法本身已難以將搜索快速引導至目標值整體較優的不同區域(即較優區域)執行,容易導致算法實際搜索效率下降.因此,近年來部分學者在所提算法中,先采用某種分解策略將問題合理分解為多個相對簡單的子問題,從而將算法搜索盡量限定于解空間中部分存在優質解的較優區域之內,然后再用智能算法等求解各子問題并獲得原問題的解,取得良好效果.在MVRP方面,主要考慮地理條件、客戶需求和交貨時間等因素,并采用聚類算法把MVRP 分解為一系列單個車場VRP,然后采用粒子群優化(Particle swarm optimization,PSO)算法[11?12]、遺傳算法[13]等求解.在HFVRP方面,主要考慮車輛載重量、車輛運輸費用和客戶需求等因素,并采用路徑分割法把HFVRP 分解為一系列單車輛VRP (即旅行商問題(Traveling salesman problem,TSP)),進而采用禁忌搜索算法[14]等求解.在MHFVRP 方面,由于該問題復雜度高,為HFVRP 和MVRP 的綜合,故大都先采用分解策略將MHFVRP 分解為一系列HFVRP (即單車場異構車輛的VRP),然后采用路徑分割法將各HFVRP 進一步分解為多個TSP.譬如,Dondo等[15]先采用一種啟發式聚類算法將全部客戶分成若干類,然后利用路徑分割算法分配車輛,從而將原問題轉化為多個小規模TSP,然后再利用混合整數規劃求解器進行求解.Tang等[16]先利用最近鄰方法把客戶聚類到相應的車場,然后再用掃描法對每個車場中的客戶進行路徑分割,進而將原問題轉化為多個TSP,最后采用蟻群算法求解.上述文獻中的仿真實驗和算法對比驗證了分解策略和智能算法結合的必要性.然而,對于現實物流配送中大量存在且更為復雜的LMHFVPR_TW,尚無結合分解策略的智能求解算法,故開展相關研究意義重大.

蟻群算法(Ant colony optimization,ACO)是一種模擬螞蟻尋找食物的群智能優化算法,最早由Dorigo等[17]提出并首次用于求解TSP.ACO 借鑒螞蟻覓食的群體智能行為,利用信息素濃度矩陣(以下簡稱信息素矩陣)學習和保留螞蟻在各自行駛路徑(即問題解)上留下的信息素濃度信息,并基于該矩陣構建路徑轉移概率矩陣(以下簡稱概率矩陣),然后通過對概率矩陣采樣生成新種群來實現路徑搜素(即解空間搜索)并引導搜索方向.ACO 采用正反饋機制使得較短行駛路徑(較優解)的優良信息得以更多積累,可引導算法在較短時間內到達解空間中存在優質解的區域.這使得算法具有較強的全局搜索能力,從而在多種VRP 上得到成功應用[7?8, 18?24].從ACO 在VRP 系列問題上的研究現狀來看,采用合理的信息素濃度計算策略可控制全局搜索的方向,引入基于有效鄰域操作的局部搜索能進一步提高算法性能,兩者是設計有效ACO 的關鍵.文獻調研表明,ACO 求解LMHFVPR_TW 的研究目前處于空白狀態.

本文研究LMHFVPR_TW 的建模與求解.在建模方面,給出包含經濟費用、燃油消耗費用和客戶滿意度費用的運輸總費用計算模型,并建立以最小化運輸總費用為優化目標的LMHFVPR_TW.在求解方面,考慮到LMHFVPR_TW 這類問題的解空間龐大且整體編碼解碼復雜,直接采用智能算法難以在較短時間內實現有效搜索,故提出一種結合聚類分解的增強蟻群算法(Enhanced ant colony optimization based on clustering decomposition,EACO_CD)進行求解.

EACO_CD 的結構如圖1 所示(以2 個車場2 類車型為例).由圖1 可知,EACO_CD 由問題分解階段和子問題求解階段組成.1)在問題分解階段,為確保分解后的各子問題區域盡量覆蓋解空間中的優質解區域,設計兩層基于K-means 的聚類算法對問題進行逐步分解.第1 層聚類算法為所提的一種改進平衡K-means 聚類算法(Improved balanced K-means algorithm,IBKA),用于給每個車場分配一定數量的客戶,從而形成一系列帶時間窗的低能耗單車場多車型VRP (LHFVRP with time windows,LHFVRP_TW);第2 層聚類算法為所提的一種帶粒子群優化的混合K-means 聚類算法(Hybrid K-means clustering algorithm with PSO,HKMA),用于給每個LHFVRP_TW 中的每類車型合理分配一定數量的客戶,進而得到一系列帶時間窗的低能耗VRP (Low-energy-consumption VRP with time windows,LVRP_TW).兩層聚類分解算法將原問題最終分解為子問題LVRP_TW,而不是分解為比LVRP_TW 規模更小的帶時間窗的低能耗TSP,這使其子問題的解空間區域相對較大,有利于算法對原問題解空間中更多區域進行搜索,可望獲得更好的解.2)在子問題求解階段,采用所提的增強蟻群算法(Enhanced ant colony optimization,EACO)對每個分解后子問題(即LVRP_TW)的解空間區域進行搜索,進而將子問題的解合并后得到原問題的解.在EACO 中,采用動態的信息素濃度揮發系數(以下簡稱信息素揮發系數),并進一步加入信息素揮發系數控制因子以調節其取值,可避免算法過早收斂,并進一步引導算法全局搜索到達更多的不同區域;同時,設計基于4 種變鄰域操作的兩階段局部搜索策略,用于對全局搜索發現的優質解區域進行細致且高效的搜索,以平衡全局和局部搜索并增強算法性能.最后,通過仿真實驗和算法比較驗證所提算法的有效性.

圖1 EACO_CD (EACO_IBKA_HKMA)結構Fig.1 Framework of EACO_CD (EACO_IBKA_HKMA)

1 LMHFVRP_TW 建模與分析

本節建立LMHFVPR_TW 模型并對問題特點和求解算法進行分析.LMHFVPR_TW 為帶軟時間窗約束的多車場多車型車輛配送問題,優化目標為最小化運輸總費用(即經濟費用、燃油消耗費用和客戶滿意度費用之和).該問題在已知客戶位置坐標、車場位置坐標、客戶貨物需求量、客戶時間窗、各車場車型特征(表1 中給出相應的符號定義)的情況下,希望找到各車場中車輛的最佳配送路線(由第1.3 節問題模型中的決策變量具體取值確定),從而使總運輸費用達到最優.LMHFVPR_TW的應用場景很多,譬如電商和連鎖超市的商品配送、戰區和疫區的緊急物資配送、石化企業成品油的一次和二次配送等.

1.1 問題假設和符號定義

LMHFVPR_TW 滿足如下假設:

1) 每個車場當中的每一類車型都參與配送任務;

2) 每個車場所含的車型種類相同;

3) 車輛在一次任務分配中僅進行一次配送;

5) 車輛從車場出發,完成服務后需回到原車場;

6) 每個客戶都被服務,被服務的次數有且僅有一次;

7) 車輛配送中的貨物需盡量在每個客戶要求的到達時間段內送達.

LMHFVPR_TW 的相關符號定義如表1 所示.

表1 符號及定義Table 1 Symbols and definitions

1.2 問題優化目標計算模型

LMHFVPR_TW 的優化目標為運輸總費用ZTotal,該優化目標由3 部分組成: 經濟費用ZEmission-cost、排放費用ZEmission-cost及客戶滿意度費用ZCustomer-Satisfaction.其中,經濟費用由運輸距離費用F1和車輛固定成本F2組成,排放費用設定為燃油消耗費用F3,客戶滿意度費用由軟時間窗懲罰費用F4構成.

其中,式(3)中燃油消耗量FUMij的計算基于文獻[10]中的綜合燃油消耗模型和參數設置,并假設各類車型的油耗參數相同,可得到FUMij的計算式為

1.3 問題模型

問題優化目標為

其中,式(9)要求從某車場出發的車輛數目和回到該車場的車輛數目相等;式(10)和式(11)要求每個客戶只能被一輛車服務1 次;式(12)要求所有車輛均需參與配送,且配送完后返回原車場;式(13)要求生成的配送方案或問題解不出現不含車場的子回路;式(14)要求車輛配送過程中不能超過其對應車型的額定載重量.

1.4 問題特點分析

由第1.3 節可知,本文的LMHFVPR_TW 為具有非線性和大規模約束的0-1 整數規劃問題.相對于傳統VRP,式(1)～ (4)使優化目標不同且更加復雜,同時式(4)使優化目標具有非線性;多車場和多車型的引入,使得決策變量由傳統VRP 的3維變量xijk變為5 維變量xP Mijk,導致決策變量和約束的數量增多,這使得問題解空間更加龐大和復雜.此外,相對于常規0-1 整數規劃問題,VRP 系列問題均含子回路消除約束(式(13)),該類約束的數量為 2N+1?2,這種指數級的約束數量限制了問題的求解規模.雖然可以通過增加連續決策變量,使式(13)的約束個數變為客戶數N的多項式,從而擴大問題的求解規模,但此時原問題近似轉化為同時含離散和連續決策變量的混合0-1 整數規劃問題,該近似問題的求解難度并未降低[26].

1.5 問題求解算法分析

VRP 系列問題(包括本文的LMHFVPR_TW)均可建模為0-1 整數或混合整數規劃模型,其求解算法分為兩類.一類為運籌學算法,包括分支定界、分支切割、動態規劃等[27].該類算法主要以線性代數和幾何分析為基本工具,利用問題優化目標函數和約束條件式的結構信息構造遍歷或部分遍歷解空間的搜索,可在幾分鐘至幾十分鐘內求解較小規模問題(客戶數小于等于20),同時經過合適設計,也可用于求解較大規模問題(客戶數大于等于80),但往往求解時間很長[27?29].不同于線性規劃問題,0-1規劃問題具有非凸特性,其內在幾何結構與最優解間的關系仍屬開放問題,目前僅在無約束0-1二次規劃等少數具有特殊結構的簡單問題上存在有效多項式時間求解算法[6,28,30],尚無通用的最優解多項式時間求解算法.對于LMHFVPR_TW 這類復雜的0-1 規劃問題,采用已有的運籌學算法框架雖然也可設計相應的求解算法,但算法難以在較短時間獲取較大規模問題的滿意解.另一類是智能優化算法(以下簡稱智能算法),包括蟻群優化算法、遺傳算法、禁忌搜索算法等[27].該類算法不依賴于問題結構,而是將問題的約束條件隱式地包含在所設計的解的編碼、解碼規則中處理,并利用某種擬人、擬物的機制不斷生成新的可行個體或解,從而引導算法執行搜索,往往可在幾秒或十幾秒內就能獲得各類VRP的滿意解[6, 27, 31].

目前,各類文獻對所提智能算法使用少量個體(一般種群規模為 20～100)運行少量代數(一般進化代數為 100～1000)后,為何都可獲得問題的滿意解,基本都只從算法本身的機制進行解釋.這導致現有研究過于追求新方法層面的創新.實際上,智能算法之所以有效,是由其編碼解碼規則和智能算法自身機制共同決定的.智能算法并不直接對0-1變量編碼,而是對客戶序號排列進行編碼解碼規則的設計.這樣的排序編碼解碼,容易將解中元素的取值限定在滿足約束的可行范圍,從而避免繁雜的約束處理,可明顯提升算法搜索效率.同時,排序編碼解碼規則所確定的解空間極為扁平,問題目標值的變化范圍遠遠小于排序模型解空間的規模,這非常有利于算法短時間內獲取滿意解.具體來說,譬如優化目標為最小化總行駛距離的單車輛VRP(即TSP),假如有100 個客戶,任意兩個客戶間的距離為在區間[1,1000]上均勻分布的隨機數,且要求車輛從某個客戶駐地(車場)出發,服務完客戶后返回該駐地,則其目標值變化范圍在(1,1000000)之內(1000000 為各客戶間的距離之和的上限,實際的最大目標值小于此值),而解空間規模為99!(車輛出發位置的客戶需扣除),平均每一個具體目標值對應約 9×10149個不同排列或解,這表明數量巨大的不同解具有相同的目標值.對于其他更為復雜的VRP (包括本文的LMHFVPR_TW),也存在這一情況.另外,離散問題本身解之間沒有梯度,相鄰解之間的目標值差別可能很小,也可能較大.排序編碼解碼規則所確定解空間的 “極為扁平”性和“無梯度” 性,使得各種智能算法即使只用幾秒(CPU為 2.6 GHz 的主流PC 在 1 s 內可搜索并評價含100個客戶的傳統VRP 問題的數萬個解)搜索解空間中極小的區域(數萬的數量級相對于99! 類似于1根針或更小的面積相對于體育場),也可能到達較廣的目標值區域,同時其內在的尋優機制可驅動算法到達目標值較優的不同區域搜索,從而獲得滿意解.智能算法這種通過對解空間極小區域的搜索來實現對目標值較廣且較優區域的搜索,是其有效的本質原因,也是現有運籌學算法這類偏遍歷或部分遍歷解空間的算法難以做到的.因此,設計智能算法以實現對LMHFVPR_TW 快速有效求解,是合理且必要的.

由第1.4 節和第1.5 節的分析可知,對于LMHFVPR_TW 這類非凸、非線性問題,直接采用運籌學方法在短時間內獲取滿意解難度很大,故本文設計一種智能算法(即結合聚類分解的增強蟻群算法EACO_CD)進行求解.

2 問題求解算法EACO_CD

本文所提的求解算法EACO_CD 包含問題分解階段和子問題求解階段(見圖1).通過兩層聚類分解,原問題LMHFVPR_TW 轉換成為Pt×Mt個單車場單車型子問題LVRP_TWs,然后采用本文設計的子問題求解算法EACO 依次對每個LVRP_TW 進行求解,進而可獲得原問題的解和優化目標值.第2.1 節介紹EACO_CD 中問題分解階段的細節,包含兩層聚類算法的細節和復雜度分析;第2.2 節介紹EACO_CD 中子問題求解階段的細節,包含EACO 的細節和復雜度分析;第2.3 節介紹EACO_CD 的總體框架并分析整體復雜度.

2.1 EACO_CD 的問題分解階段

EACO_CD 的問題分解階段執行兩層聚類算法,第1 層聚類算法為IBKA,用于將原問題LMHFVPR_TW 分解為Pt個單車場子問題LHFVRP_TWs;第2 層聚類算法為HKMA,用于將每個單車場問題LHFVRP_TW 進一步分解為Mt個單車場單車型子問題LVRP_TWs (見圖1).在對數值型數據進行聚類時,大都采用K-means、模糊聚類(Fuzzy cluster)等算法[14, 32?34].而針對混合型數據,可采用K-prototypes 算法或先將混合型數據轉換為數值型數據[35].由于本文問題中的客戶屬性數據均為值型,不需要對數據類型進行轉換,同時,Kmeans 具有線性計算復雜度,已成功應用于VRP系列問題[12, 14, 33, 36?37],故本文在IBKA 和HKMA中采用K-means 作為基本分解策略.

2.1.1 第1 層聚類算法(IBKA)

本文問題的優化目標(見式(8)) 由4 部分組成,其中有兩部分與車輛運輸距離有關.當車輛距客戶較遠時,配送所產生的燃油消耗費用和距離費用都隨之增加(見第1.2 節模型),故IBKA 采用客戶之間的歐氏距離作為評價標準,將客戶聚成與車場數目相同的幾類,使得每個車場服務其中一類客戶.但是,隨著客戶數量逐漸增大,聚類后會出現每個車場服務的客戶數量相差較大.為確保每個車場的客戶資源能均衡分配,設計IBKA 進行聚類.由于文獻[37]提出的K-means 平衡算法僅適用于2 類客戶之間的平衡(Pt=2),本節將文獻[37]提出的Kmeans 平衡算法推廣客戶群,同時設計車場分配規則,可更為合理地確定各車場服務的客戶群.進而,獲得Pt個單車場子問題LHFVRP_TWs.IBKA 的步驟如下.

步驟1.對全體客戶進行K-means 聚類.初始化聚類重心為每個車場的坐標位置,初始化聚類數量為車場數量,以歐氏距離作為評價指標,以迭代次數作為終止條件.

步驟2.平衡客戶群.完成上述步驟后進行k類客戶群之間的平衡計算,具體如圖2 所示,按每類客戶群數目由大到小排序,根據兩類客戶間的邊緣算法[37],將客戶較多的客戶群依次向其余客戶較少的客戶群移動,循環操作直至各類客戶群中客戶數目達到平衡.

圖2 4 類客戶平衡移動示意圖Fig.2 Diagram of balanced movement for four customer groups

步驟3.計算每個車場到各聚類重心的距離.列出距離矩陣A,設A為Pt×Pt的二維矩陣,矩陣中的每一個元素表示某一車場坐標到某聚類重心的距離,即

其中,d(i,j) 為車場坐標i到聚類重心坐標j的歐氏距離.

步驟4.按設計的分配規則對車場所服務的客戶群進行分配.選取矩陣A中不同行(車場不同)不同列(客戶群重心不同)的Pt個元素并相加,共有Pt! 種不同的組合及對應的累加和選出其中最小的Smin(見式(17)) 所對應的組合中每個元素d(i,j) 的行、列標號即為車場i所服務的客戶群j. 采用Smin確定車場服務的客戶群有利于整體減小車輛運輸距離.

其中,Sl表示矩陣A中對不同行不同列元素的第l種累加和.

對于K-means 迭代次數為gen1K、平衡迭代次數為genb、問題規模為N_Pt(客戶數量_車場數量) 的問題,IBKA 的計算復雜度(以下簡稱復雜度)TIBKA分別由K-means 聚類復雜度O(gen1K×Pt×N)、客戶平衡復雜度 O (genb×(Pt×N)2) (最壞情況)、獲取方陣A的復雜度 O () 和獲取Smin的復雜度 O (Pt×) 決定.因此,

針對3 車場和144 客戶的本文問題LMHFVPR_TW,采用IBKA 進行車場劃分.根據列出距離矩陣A(3×3),分別給出矩陣A(3×3)中不同行不同列的Pt個元素累加和的全部組合方式,并計算不同的累加和共3! 次(見式(18)).若Smin=S1,即車場1 服務A類客戶,車場2 服務B類客戶,車場3 服務C類客戶,劃分結果如圖3 所示.圖3 中上、下子圖分別為未采用和采用平衡聚類的結果.從圖中可以看出,采用平衡聚類可使各車場服務的客戶數量更加接近,有利于提高各車場的整體效率和控制子問題規模.

圖3 3 車場K-means 未平衡聚類與平衡聚類比較Fig.3 Comparison of unbalanced K-means cluster and balanced K-means cluster of three depots

2.1.2 第2 層聚類算法(HKMA)

客戶的時間窗、客戶的貨物需求重量、客戶與車場間的距離是影響本文問題優化目標(見式(8))的重要因素,故HKMA 將其作為客戶的3 類屬性,采用K-means 對客戶聚類后再為各類客戶分配相應的車型.K-means 聚類時是將每個客戶作為一個對象(質點).在多屬性客戶進行K-means 聚類時,需將各屬性做加權和后才可進行.然而,上述3 類屬性對優化目標(見式(8))的影響程度無法直接確定,故需對各屬性權重系數進行優化,從而提高子問題分解的合理性.

HKMA 為內外兩層結構,外層執行微粒子群優化(PSO) 算法,用于優化屬性權重系數向量λ(λ=(λ1,λ2,λ3)且λ1+λ2+λ3=1);內層對基于每個λi的聚類分解效果進行評價,并將評價值fi(λi)作為PSO 種群個體的適應度值.HKMA的工作機制如圖4 所示,其中PSO 種群PPSO的每個λi(i= 1,···,M1)表示一組權重系數取值,基于λi對單車場問題LHFVRP_TW 做K-means 聚類可得到Mt類客戶群,將Mt類車型和Mt類客戶群逐一配對,得到單車場單車型子問題LVRP_TW(記為Ai,j,i和j分別為車型和客戶群編號,i,j=1,···,Mt),f(Ai,j) 表示對問題Ai,j采用第3.1 節的編碼規則隨機生成10 個解的評價值(即子問題總運輸費用)的平均值,Fi為f(Ai,j) 組成的方陣,Sum_fk(k=1,···,Mt! ) 表示第k個由Fi中Mt個相互不同行(車型不同)不同列(客戶群不同)元素累加所得的值,fi_min 取所有Sum_fk中的最小值并賦值給fi(λi). 由Mt個車型和客戶群均不同的Ai,j可構成一個LHFVRP_TW,故Sum_fk實際上是對第k種子問題分解方式合理性的評估,同時fi_min用于標識相對最為合理的子問題分解方式.

圖4 HKMA 工作機制Fig.4 Running mechanism of HKMA

令ZTotal(Ai,j) (見式(8))為個體Ai,j的目標函數值,λopt為PPSO的歷史最優個體,subP(λopt)為f(λopt) 對應的Mt個單車場單車型問題LVRP_TWs 的集合(具體見前述),Maxgen(PSO) 為HKMA 中PSO 算法的運行代數,gen2K為HKMA中每次執行K-means 算法的運行代數.HKMA 的具體步驟如下.

步驟1.令gen=1,隨機生成PPSO中的每個λi.

步驟2.(內層): 確定PPSO中所有λi的適應度值(具體見前述).

步驟3.(外層): 若gen>1,則將當代PPSO和父代PPSO的個體進行冒泡排序,并擇優保留前M1 個組成PPSO. 更新λopt和subP(λopt).

步驟4.(外層):PPSO通過PSO 算法進化后得到新的種群PPSO,也即得到M1 組新的權重系數向量.

步驟5.令gen=gen+1,若gen≤Maxgen(PSO),則返回步驟2,否則輸出subP(λopt).

在本文算法EACO_CD 中,HKMA 的參數設置為M1=10,M2=10,Maxgen(PSO) =10和gen2K=20,PSO 中的慣性權重和加速因子的設置與文獻[38]相同,具體設置為慣性權重ω=0.9 ,速度因子c1=1.2,c2=1.5.針對第2.1 節中示例問題的車場1 和其服務的A類客戶(即IBKA 劃分得到一個單車場子問題LHFVRP_TW),采用HKMA 進行車型劃分,可得圖5 和圖6.由圖5 和圖6 可見,A類客戶依據車型和客戶的3 類屬性得以進一步劃分,從而將1 個LHFVRP_TW 分解為規模更小的2 個單車場單車型子問題LVRP_TWs.

圖5 HKMA 三維聚類效果Fig.5 The 3D clustering results of HKMA

圖6 HKMA 二維結果Fig.6 The 2D results of HKMA

2.2 EACO_CD 的子問題求解階段

通過上一階段的兩層聚類分解后,原問題LMHFVPR_TW 轉換成為Pt×Mt個單車場單車型子問題LVRP_TWs,從而本階段可采用本文設計的子問題求解算法EACO 依次對每個LVRP_TW進行求解,進而可獲得原問題的解和優化目標值(見圖1).第2.2.1～ 2.2.4 節介紹EACO 的細節,第2.2.5節給出EACO 的復雜度分析.

2.2.1 編碼與解碼規則

EACO 對各車輛服務的所有客戶進行整體編碼.每只螞蟻代表多輛車,每只螞蟻的行駛路徑(即子問題LVRP_TW 的l個解)為對應各車輛的總行駛路徑.譬如對于第l只螞蟻服務10 個客戶,可記為:πant(l)=[1,4,9,5,10,2,3,6,7,8].

EACO 采用文獻[39]中的方法對每只螞蟻的行駛路徑進行解碼,即逐一選擇車場中的每輛車,將πant(l)中客戶從左向右依次加入當前車輛,如果當前車輛的載重量大于其服務客戶的總需求量,則繼續往該車輛中加入客戶,否則選擇下一輛車繼續依次加入剩余客戶.譬如,對于上一段的πant(l),假設需兩輛車進行服務,解碼后可得πant(l)=[πvehicle(1),πvehicle(2)]=[1,4,9,5],[10,2,3,6,7,8].

顯然,當采用上述編碼解碼規則時,式(13)的子回路消除約束自然得到滿足,同時也不違反其他約束.

2.2.2 全局搜索

2.2.2.1 初始化

首先,利用掃描法(Sweep algorithm,SWA)構造1 只螞蟻的行駛路徑[40],并采用該路徑對信息素濃度τij進行初始化

其中,i,j為客戶編號,Pm為信息素濃度值的初始值(設置為大于1 的值).為防止算法過早陷入局部最優[41],設置信息素濃度τij的最大值和最小值為

然后,設置信息素濃度增量 ?τij=0 以及種群大小popsize=(2/3)×N.

2.2.2.2 螞蟻行駛路徑搜索

EACO 每代通過對當前種群中各螞蟻行駛路徑的確定來實現對問題解空間的全局搜索.每只螞蟻均利用式(22)確定其行駛路徑.

其中,i為當前客戶,u為下一客戶,為第l只螞蟻從i到u的轉移概率,s為螞蟻l還未服務過的客戶,τiu(τis) 為i和u(s) 間的信息素濃度,ηiu(ηis)為i和u(s) 之間距離的倒數,α和β為重要程度權重,tabu為螞蟻已經搜索過的客戶地點集合.對tabu}進行輪盤賭選擇,即可得到u的取值,從而確定螞蟻l從客戶i出發需到達的下一客戶.

2.2.2.3 信息素濃度更新

EACO 每代對每只螞蟻的信息素濃度采用式(23)～ (26)依次進行更新

其中,γ為信息素揮發系數控制因子(初始設置為1),Ma和Mb為迭代次數系數(均設置為5),T為基準值(設置為50),Lb為當前最好解對應的目標函數值,t表示算法進化的代數(即算法第t次迭代),ρt為第t代的信息素揮發系數(初始設置為0.9),ρmin和ρmax分別為ρt的最小值和最大值(分別設置為0 和1),為第t代的信息素濃度增量,W為信息素增量常數(設置為500),為第t代的信息素濃度.γ根據算法每代的實際運行效果動態調整其取值,可在運行效果好(差)時增加(減小)信息素濃度,有利于增強算法合理引導搜索的能力.

2.2.3 局部搜索

本節在基本變鄰域搜索(VNS)[42]基礎上提出一種兩階段變鄰域局部搜索策略(Two-stage variable neighborhood search,TVNS),對EACO 每代發現的當前最好解進一步執行高效且細致的鄰域搜索,以提高解的質量.

TVNS 為在車輛間和車輛內部分別執行多種鄰域操作(見圖7)的兩階段局部搜索,兩個階段中各鄰域操作均在當前最好解上執行且每次執行的最大次數均設置為genlocal=20.首先,第1 階段隨機選擇兩輛車,依次對車輛間的客戶執行 “Insert”和“Exchange”鄰域操作,如果在執行操作期間發現更好解,則更新當前最好解并立刻進入第2 階段,否則各鄰域操作均執行至最大次數genlocal后進入第2 階段.然后,第2 階段逐一選擇每輛車,對每輛車的內部客戶依次執行 “2-opt”、“Insert”、“Exchange”和 “Swap”鄰域操作,如果在執行操作期間發現更好解,則更新當前最好解并立刻對下一輛車執行操作,否則各鄰域操作均執行至最大次數genlocal后再對下一輛車執行操作.第2 階段結束后完成本次局部搜索.

圖7 局部搜索策略Fig.7 Local search strategy

2.2.4 EACO 終止條件

設定終止條件為算法的運行時間,如果滿足運行時間要求,則輸出服務所有客戶的每輛車的行駛路徑以及相應的行駛總費用.

2.2.5 EACO 復雜度分析

其中,l ～ogN表示低于N的亞線性復雜度.

2.3 EACO_CD 整體結構及復雜度分析

2.3.1 EACO_CD 整體結構

算法EACO_CD 的整體結構如圖1 所示(以2 個車場2 類車型為例).由圖1 可知,EACO_CD對Pt×Mt個子問題LVRP_TWs 依次執行EACO進行求解,最終將獲得的各子問題解合并,得到原問題的解,同時將各子問題解的優化目標值累加,得到原問題的優化目標值,從而實現對原問題的求解.

2.3.2 EACO_CD 復雜度分析

EACO_CD 是對本文問題LMHFVPR_TW的整體求解,其復雜度TEACO_CD由分解算法IBKA的復雜度TIBKA(見第2.1.1 節)、分解算法HKMA總的復雜度total(THKMA) (見第2.1.2 節)、求解算法EACO 總的復雜度total(TEACO) 組成(見第2.2.5節).即

雖然TIBKA,total(THKMA)和total(TEACO) 中的變量較多,但除N(客戶數量)、Pt(車場數量)、Mt(車型數量)與問題規模相關,Z(局部搜索每次迭代搜索的鄰域或個體)與具體問題相關,popsizeE設置為2/3×N(見第2.2.2.1節)以外,其余變量均與算法相關且實際使用中都設置為常數.因此,可得簡化式如下:

顯然,所有變量的多項式次數不超過2,影響TEACO_CD的是TIBKA中的 O(Pt!) 和total(THKMA) 中 的O(Pt×Mt!). 然而,實際問題中Pt和Mt一般都不大于10,而若僅計算10!,用CPU 為2.6 GHz 的主流PC 可在1 s左右完成計算,故EACO_CD 對大多數實際問題可在較短時間內結束運算.另外,如果Pt(Mt)很大(譬如為50),可將IBKA 的步驟4 (HKMA 的步驟2)中窮舉不同行不同列的組合改為采樣Pt(Mt)的多項式次數的組合,從而擴大EACO_CD 求解問題的規模.

3 實驗比較與分析

本節所有算法的測試均采用以下軟硬件配置:英特爾I7 處理器(3.2 GHz),8 GB 內存,Win10 操作系統,MATLAB2018a 編程環境.

3.1 EACO 參數設定

在第1 節所提出的目標函數中(見式(7)),距離費用系數設定參考文獻[7],車輛固定費用參考文獻[9],燃油費用系數設定參考文獻[43],時間窗懲罰費用系數參考文獻[16],具體取值見表2.

表2 目標函數中的相關系數Table 2 Coefficients in the object function

EACO 的殘留信息素重要程度α,啟發信息素重要程度β,初始化信息素濃度參數Pm,信息素增量常數W為4 個主要的參數,為確定合適的參數組合,對這4 個參數進行正交實驗,各參數設置水平如表3 所示.由于聚類分解后原問題LMHFVPR_TW 轉變為Pt×Mt個LVRP_TWs,故采用Solomon 標準VRPTW 數據集中的問題c101 (100 個客戶) 進行測試.每組參數組合下的EACO 在c101 上獨立運行20 次,取20 次的平均最小運輸總費用值作為平均響應值ARV,參數設置的正交表見表4,參數的平均響應值和影響力見表5.由表4和表5 可得,EACO 的參數設置為:α=1.25,β=2.5,Pm=1.1,W=500.

表3 主要參數與水平Table 3 Main parameters and level

表4 參數設置的正交表Table 4 Orthogonal table of parameter settings

表5 各參數不同水平下的平均響應值和影響力Table 5 Average response values and influences table at different levels of each parameter

3.2 實驗設計及測試問題選取

3.2.1 實驗設計

EACO_CD 由IBKA、HKMA 和EACO 組成.為驗證EACO_CD 的有效性,先將這3 種算法分別與國際期刊上的相近算法進行比較,最后再進行總體比較.測試實驗車輛參數見表6.

表6 4 種不同車型相關參數設置Table 6 Related parameter settings for four different vehicle types

首先,為驗證EACO 的有效性,在多車場單車型問題LMVRP_TW 上,將EACO1與國際期刊中的有效算法DHACO[37]進行對比.EACO1 為EACO_CD 不含兩層聚類的算法,而DHACO 為不含分解且直接求解整個問題的一類蟻群算法.EACO1 采用DHACO 的編碼.測試結果見表7.

其次,為單獨驗證IBKA 的有效性,在多車場單車型問題LMVRP_TW 上,將EACO_IBKA與EACO_KM、EACO_NNA 進行對比.這3 種算法為在EACO1 中分別加入3 類第1 層聚類算法(即IBKA、K-means[12, 14]、NNA[11])后得到.測試結果見表7.

表7 EACO_IBKA與其他算法對比結果Table 7 Comparison results of EACO_IBKA with other algorithms

然后,為單獨驗證HKMA 的有效性,在單車場多車型問題LHFVPR_TW 上,將EACO_ HKMA 和EACO_RDA、EACO_RAA、EACO_ KEW、EACO2、TSA_RDA[10]進行對比.前4 種算法為在EACO1中分別加入4 類第2 層聚類算法(即HKMA、RDA[10]、RAA、KEW)后得到,其中RDA根據客戶需求量優先選擇最大載重量的車型,進而將問題分解為一系列單車輛LVRPs (即TSPs),RAA隨機將客戶劃分給某種車型,KEW 給每類客戶屬性分配相同權重后用K-means 劃分車型(即簡化的HKMA),RAA 和KEW 均將問題分解為一系列單車型LVRP_TWs;EACO2 為EACO1 改用多車型編碼后得到,用于直接求解整個問題;TSA_RDA 為國際期刊中的有效算法,該算法在禁忌搜索算法中加入RDA.測試結果見表8.

表8 HKMA與其他劃分算法的對比結果Table 8 Comparison results of HKMA and the other dividing algorithms

最后,為驗證EACO_CD 的有效性,在本文問題LMHFVPR_TW 上,將EACO_CD 和國際期刊中的有效算法IACO_CD[16]和IHGA[13]進行對比.IACO_CD 為兩層分解算法,該算法第2 層將問題分解為單車輛LVRP (即TSPs),而IHGA 通過加入車型編碼后對問題進行整體求解.測試結果見表9.

3.2.2 測試問題選取

本文的測試數據來源于網站(http://neo.lcc.uma.es/vrp/)中的MDVRP 數據集(數據中不存在明顯聚類特征的客戶位置分布).其中包括問題pr01 (48 個客戶),pr02 (96 個客戶),pr03 (144 個客戶),pr04 (192 個客戶),pr05 (240 個客戶),pr06(288 個客戶),p22 (360 個客戶).由于該數據中未含有車型數據,故在HKMA 和EACO_CD 的有效性驗證時加入車型數據(見表6),各比較算法在每個問題上獨立運行20 次,每個算法設置相同運行時間T(s).性能指標為各算法20 次運行中的最優值(Best)、平均值(Average)、最差值(Worst)和標準差(SD).各種帶分解策略的比較算法每次運行時間為問題分解算法與問題求解算法運行時間的總和,譬如,EACO_CD 每次運行的時間為EACO、IBKA和HKMA 的運行時間總和,表7～ 9 中粗體為算法運行結果的較小值.

3.3 實驗結果比較和分析

3.3.1 驗證EACO 有效性

由表7 可知,EACO1 (即不含問題分解策略的EACO_CD)的最優值和平均值均優于DHACO.在算法全局搜索部分,DHACO 和EACO1 的復雜度無明顯差異,但EACO1 在全局階段引入信息素揮發系數控制因子,進一步動態調節信息素揮發系數,可有效控制信息素的揮發,從而提高了算法的全局搜索能力.在算法局部搜索部分,相比DHACO的VNS 策略,EACO1 的TVNS 能有效搜索更多不同區域,具有更強的局部搜索能力.

3.3.2 驗證IBKA 有效性

由表7 可知,EACO_IBKA 總體優于EACO_KM 和EACO_NNA.這表明IBKA 通過聚類和平衡客戶群,可以更合理地分配車場并較好實現子問題解耦.同時,EACO_IBKA 也總體優于不含問題分解策略的整體求解算法EACO1 和DHACO.這說明采用有效的問題分解策略為算法提前確定需搜索的較優區域,可在一定程度上避免過多的低效搜索,有利于增強算法性能.此外,測試結果顯示,帶分解策略的算法更適用于較大規模問題,原因在于這些問題的解空間更加龐大,僅依靠智能算法自身進化機制難以快速引導算法發現優質解區域.

3.3.3 驗證HKMA 有效性

由表8 可知,EACO_HKMA 總體優于其他比較算法.這表明HKMA 綜合考慮客戶的3 類屬性并為各屬性設定合適的權重,可以更合理地劃分子問題.具體而言,EACO_HKMA 明顯優于EACO_RDA 和TSA_RDA,這說明RDA 將問題分解為更小規模的單車輛LVRPs (即TSPs)會將算法實際搜索空間限定得過小,從而遺漏較多存在優質解的區域;EACO_HKMA 明顯優于EACO_RAA(除N=96的1 個較小規模問題)和EACO_KEW,這說明RAA 隨機劃分子問題的方式和KEW 給各屬性分配相同權重后再劃分子問題的方式,均無法獲取真正優質的搜索區域;EACO_HKMA 也明顯優于整體求解算法EACO2,這再次驗證了將問題先合理分解可把算法搜索直接限定在較優搜索區域進行,有利于提高算法效率.此外,與第3.3.2 節類似,本節的分解策略同樣對較大規模問題更有效.

3.3.4 驗證EACO_CD 有效性

由表9 可知,EACO_CD 總體優于IACO_CD和IHGA.具體而言,EACO_CD 的最優值、平均值、最差值均明顯優于IACO_CD,且標準差整體較小,這說明IACO_CD 將問題分解為更小規模的單車輛LVRPs (即TSPs)會遺漏較多存在優質解的區域,同時IACO_CD 自身算法的搜索能力有限,從而導致該算法的性能和魯棒性較差;EACO_CD 的最優值、平均值、最差值、標準差明顯優于IHGA (除N=96 的2 個較小規模問題),這表明隨著車場、車型和客戶的增多,本文問題LMHFVPR_TW 的可行解空間呈指數倍增加,僅采用智能算法(如IHGA)難以在較短時間內完成大范圍搜索以獲得原問題的優質解,而EACO_CD 依靠合理的問題分解能將EACO 的搜索限定在較優區域之內,進而利用EACO 較強的搜索能力可短時間內從較優區域中獲得優質解,從而使其性能和魯棒性均較強.

表9 EACO_CD 性能驗證Table 9 Performance verification of EACO_CD

綜上可知,結合IBKA 和KHMA 的EACO_CD是求解本文問題LMHFVPR_TW 的有效魯棒算法.

4 結束語

針對多車場多車型綠色車輛路徑問題LMHFVRP_TW,本文提出一種結合聚類分解策略的增強蟻群算法(EACO_CD) 進行求解.這是首次采用智能優化算法求解該問題.EACO_CD 包括兩個階段.第1 階段為問題分解,該階段設計兩種聚類分解策略IBKA 和HKMA,用于對LMHFVRP_TW 進行分解,從而有效控制問題規模并快速引導算法在較優區域內搜索;第2 階段為子問題求解,該階段采用EACO 對每個分解后的子問題LVRP_TW 進行求解.在EACO 的全局搜索中,通過加入信息素揮發系數控制因子動態調節信息素揮發系數的取值,可較好地避免算法出現過早收斂,進而引導算法搜索達到解空間中更多的不同區域.同時,為增強EACO 的局部搜索能力,設計一種兩階段變鄰域搜索對全局搜索獲得的優質解區域進行細致且高效的搜索.仿真實驗和算法比較驗證了所提EACO_CD 是求解LMHFVRP_TW 的有效算法.對于HFVRP 和MHFVRP,已有的問題分解方法大都將其最終分解為一系列更為簡單的單車輛VRP(即標準TSP),這一分解方式是否合理且唯一,目前尚無人探討.本文所提的IBKA 和HKMA 把LMHFVRP_TW 有效分解為一系列單車型VRP(即多車輛TSP),使其子問題具有較大的可行解區域,進而適當擴大了算法的實際搜索空間,取得了更好的效果.這表明合理選取子問題可在一定程度上提高算法求解性能,也對復雜車輛路徑問題如何進行分解優化設計具有借鑒意義.后續研究將把EACO_CD 和隨機處理技術結合,用于求解動態復雜車輛路徑問題.