徐林浩，錢 斌，胡 蓉，于乃康

（1. 昆明理工大學(xué) 信息工程與自動(dòng)化學(xué)院，云南 昆明 650500；2. 昆明理工大學(xué) 機(jī)電工程學(xué)院，云南 昆明 650500）

車輛路徑問題(Vehicle Routing Problem,VRP)由Dantzig和Ramser于1959年首次提出[1]，VRP作為經(jīng)典的組合優(yōu)化問題，得到了廣大學(xué)者的研究。帶容量的車輛路徑問題(Capacitated Vehicle Routing Problem,CVRP)是VRP的一類，廣泛存在交通運(yùn)輸、物流配送等方面，同時(shí)，當(dāng)今社會(huì)推行綠色低碳發(fā)展，綠色車輛路徑問題受到了重視。在上述背景下，研究綠色帶容量的車輛路徑問題(Green Capacitated Vehicle Routing Problem, GCVRP)具有重要的現(xiàn)實(shí)意義。由于VRP為NP-hard問題，而VRP可歸約為GCVRP，故GCVRP也屬于NP-hard問題，對(duì)其展開研究具有較大的理論價(jià)值。

VRP的求解方法可以分為兩類，一類為智能優(yōu)化算法，如遺傳算法、蟻群算法和禁忌搜索算法等[2-4]，該類算法通過建立問題的排序模型進(jìn)行求解，一般不能保證得到全局最優(yōu)解以及無(wú)法判斷所求解的優(yōu)劣程度；另一類為精確算法，如分支定界、列生成和拉格朗日松弛算法等[5-10]，該類算法通過建立問題嚴(yán)格的數(shù)學(xué)模型，能在合理時(shí)間內(nèi)求得小規(guī)模問題最優(yōu)解，針對(duì)大規(guī)模問題，該類算法能在合理時(shí)間內(nèi)求得問題的精確下界，由于該類算法所求解為問題最優(yōu)解或逼近最優(yōu)解的下界，能較大程度保證解的求解質(zhì)量，同時(shí)求解時(shí)間也基本合理。

近年來(lái)，相關(guān)學(xué)者對(duì)VRP及其拓展問題的精確算法或基于精確算法的混合算法進(jìn)行了研究。Bouzid等[11]針對(duì)CVRP，以最小化路徑成本為優(yōu)化目標(biāo)，建立了整數(shù)規(guī)劃模型，并提出一種拉格朗日分裂(Lagrangian Split, LS)與變鄰域搜索(Variable Neighbourhood Search, VNS)結(jié)合的混合算法進(jìn)行求解。Singh等[5]針對(duì)具有模糊需求的CVRP，建立了該問題的數(shù)學(xué)模型，并以最小化成本為優(yōu)化目標(biāo)，提出一種分支定界算法(Branch and Bound, B&B)進(jìn)行求解。Lysgaard等[6]針對(duì)累積容量車輛路徑問題(Cumulative Capacitated Vehicle Routing Problem,CCVRP)，建立了問題的集劃分模型，以最小化路徑總成本為優(yōu)化目標(biāo)，采用(Branch-and-Cut-and-Price,BCP)進(jìn)行求解。Majid等[7]針對(duì)隨機(jī)需求的車輛路徑問題(Vehicle Routing Problem with Stochastic Demands, VRPSD)，建立了該問題的數(shù)學(xué)模型，以最小化總費(fèi)用為優(yōu)化目標(biāo)，采用整數(shù)L形(L-shaped)算法在分支切割框架下求解VRPSD，實(shí)驗(yàn)結(jié)果驗(yàn)證了所提算法的有效性。Christiansen等[8]針對(duì)VRPSD，以最小化成本為求解目標(biāo)，提出了一種分支定價(jià)算法(Branch-and-price, B&P)進(jìn)行求解。

當(dāng)今社會(huì)持續(xù)推進(jìn)節(jié)能減排，堅(jiān)持綠色發(fā)展，考慮燃油消耗和碳排放等因素的綠色車輛路徑問題(Green Vehicle Routing Problem, GVRP)受到了重視。Andelmin等[9]針對(duì)求解目標(biāo)為最小化行駛總距離的GVRP，并將該問題建模為集合劃分問題，提出了一種結(jié)合K-路(K-path)切割的列生成(Column Generation, CG)算法進(jìn)行求解。Dabia等[10]研究污染路徑問題(Pollution-routing Problem, PRP)，以最小化運(yùn)送成本為求解目標(biāo)，運(yùn)用CG算法對(duì)主問題進(jìn)行了求解，并提出一種定制的標(biāo)記算法解決了定價(jià)問題，實(shí)驗(yàn)結(jié)果驗(yàn)證了所提算法的有效性。?a?r?等[12]提出了GCVRP的混合整數(shù)規(guī)劃模型(Mixed Integer Programming, MIP)，并以最小化行駛距離為求解目標(biāo)，同時(shí)提出了一種基于模擬退火(Simulated Annealing, SA)的精確解方法進(jìn)行求解，作者運(yùn)用改進(jìn)B&P算法改進(jìn)下界，并引用了一種基于SA的啟發(fā)式算法來(lái)獲得上界，實(shí)驗(yàn)結(jié)果驗(yàn)證了所提算法的有效性。Qiu等[13]研究碳排放限額交易下的污染生產(chǎn)路徑問題(Pollution Production-routing Problem, PPRP)，建立了該問題的數(shù)學(xué)模型并以最小化運(yùn)營(yíng)成本為求解目標(biāo)，提出一種分支定價(jià)啟發(fā)式算法進(jìn)行求解，仿真實(shí)驗(yàn)證明了該模型具有降低二氧化碳排放水平和運(yùn)營(yíng)成本的潛力。由文獻(xiàn)調(diào)研可知，目前求解GCVRP的精確算法或與精確算法相結(jié)合的算法相對(duì)較少。

拉格朗日松弛算法(Lagrange Relaxation, LR)是一種求解復(fù)雜組合優(yōu)化問題的有效方法，LR主要思想是將問題的難約束引入目標(biāo)函數(shù)得到松弛問題，再通過不斷更新迭代乘子，逐步求得松弛問題的解。運(yùn)用傳統(tǒng)LR求解GCVRP，往往只能求得問題的下界，該下界一般為不可行解。因此，本文綜合LR和啟發(fā)式算法優(yōu)勢(shì)，采用兩者相結(jié)合的算法求解GCVRP，首先運(yùn)用LR求得問題下界，再根據(jù)問題性質(zhì)，設(shè)計(jì)相應(yīng)啟發(fā)式算法對(duì)下界進(jìn)行修復(fù)及優(yōu)化，得到問題的上界，通過上界和下界的間隙(Gap)，可以衡量算法性能，是一種求解GCVRP的有效方法。

綜上，本文針對(duì)實(shí)際生活中廣泛存在的GCVRP，建立以最小化運(yùn)輸成本為優(yōu)化目標(biāo)的混合整數(shù)規(guī)劃模型，并提出一種改進(jìn)拉格朗日算法進(jìn)行求解，具體為運(yùn)用拉格朗日松弛算法得到GCVRP的松弛問題，進(jìn)而得到原問題的對(duì)偶模型，并采用一種次梯度法求解該對(duì)偶模型，得到原問題的一個(gè)下界；然后，采用修復(fù)算法及鄰域搜索算法對(duì)下界進(jìn)行修復(fù)和優(yōu)化，得到原問題的一個(gè)上界。最后，通過上下界的間隙大小以及對(duì)比Gurobi求解器所求解來(lái)衡量所提算法的有效性。

1 GCVRP問題描述與數(shù)學(xué)模型

1.1 GCVRP問題描述

GCVRP是在CVRP的基礎(chǔ)上進(jìn)一步考慮了綠色因素，該問題建立在已知倉(cāng)庫(kù)和客戶點(diǎn)坐標(biāo)(坐標(biāo)以km為單位)、客戶需求及車輛載重條件下，合理調(diào)度倉(cāng)庫(kù)中的車輛為客戶送貨，使運(yùn)輸費(fèi)用最小。圖1為問題示意圖。

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

問題相關(guān)假設(shè)如下：

(1) 倉(cāng)庫(kù)內(nèi)所有車輛均參與配送；

(2) 每輛車僅執(zhí)行一次配送任務(wù)；

(3) 車輛從倉(cāng)庫(kù)出發(fā)，完成服務(wù)后返回原倉(cāng)庫(kù)；

(4) 每個(gè)客戶只能被一輛車服務(wù)一次；

(5) 每輛車配送的貨物總量不能超過其載重；

(6) 車輛在配送過程中保持勻速行駛，忽略道路等因素對(duì)車輛的影響。

GCVRP相關(guān)符號(hào)定義見表1。

表1 符號(hào)及定義Table 1 Symbols and definitions

1.2 GCVRP混合整數(shù)規(guī)劃模型

其中：式(1)為目標(biāo)函數(shù)；式(2)要求從倉(cāng)庫(kù)出入車輛數(shù)目相等；式(3)要求倉(cāng)庫(kù)內(nèi)所有車均參與配送，且配送完返回原倉(cāng)庫(kù)；式(4)～(6)要求每個(gè)客戶只能被一輛車服務(wù)一次；式(7)為流平衡約束；式(8)要求倉(cāng)庫(kù)的每輛車的載貨量不能超過其額定載重；式(9)為MTZ子環(huán)消除約束；式(10)表示車輛從倉(cāng)庫(kù)出發(fā)時(shí)的載貨量等于所服務(wù)客戶的總需求量；式(11)表示車輛完成服務(wù)任務(wù)，返回倉(cāng)庫(kù)的載貨量為0；式(12)表示某輛車服務(wù)完當(dāng)前客戶點(diǎn)后的貨物量等于從該客戶點(diǎn)離開的貨物量；式(13)和(14)為決策變量。

1.3 碳排放量計(jì)算方法

車輛行駛過程中的碳排放與燃油的消耗直接相關(guān)，影響燃油消耗的因素較多，計(jì)算也相對(duì)復(fù)雜。本文在計(jì)算碳排放量上，采用文獻(xiàn)[14-16]中的計(jì)算方法，相關(guān)參數(shù)定義見表2。計(jì)算公式為

表2 參數(shù)定義與取值Table 2 Parameter definition and value

2 改進(jìn)拉格朗日松弛算法

本節(jié)針對(duì)上節(jié)所提問題模型，提出一種結(jié)合拉格朗日松弛算法與啟發(fā)式算法的改進(jìn)拉格朗日松弛算法進(jìn)行求解。GCVRP求解流程如圖2所示。

圖2 GCVRP求解流程圖Fig.2 GCVRP solution flow chart

2.1 松弛載重約束

由于上節(jié)所提GCVRP具有NP難屬性，直接求解其最優(yōu)解非常困難，但求取該問題的下界較為簡(jiǎn)單。LR是一種求解問題下界的有效方法，通過松弛問題中的難約束，可以大大減小問題的求解難度。通過對(duì)GCVRP進(jìn)行分析，車輛的載重約束(8)影響了車輛的行駛距離，會(huì)影響目標(biāo)值的求取，松弛這條約束會(huì)使原問題容易求解。本節(jié)將約束條件(8)松弛到目標(biāo)函數(shù)中，得到拉格朗日松弛模型為

2.2 構(gòu)造可行解

求解對(duì)偶問題所得解一般是不可行解，因?yàn)榭赡苓`反車輛的載重約束條件，因此需要先對(duì)不可行解進(jìn)行修復(fù)，得到原問題的可行解，再運(yùn)用局部搜索操作對(duì)可行解進(jìn)行優(yōu)化，得到原問題高質(zhì)量解。本小節(jié)采用修復(fù)算法和鄰域搜索算法對(duì)不可行解進(jìn)行修復(fù)及優(yōu)化。

2.2.1 修復(fù)算法

2) 修復(fù)不可行解。

針對(duì)上述獲得的不可行解，運(yùn)用以下算法進(jìn)行修復(fù)，設(shè)車輛數(shù)為M，車輛編號(hào)為m，算法步驟如下：

(1) 將二進(jìn)制編碼解轉(zhuǎn)換成排序解。

(2) 依次判斷車輛m(m=1,2,···,M)是否滿足載重約束，若滿足，將車輛m放入可行車輛集合X；否則，將車輛m中超過其服務(wù)能力的客戶點(diǎn)按服務(wù)先后順序摘出，放入剩余客戶點(diǎn)集合Y中，操作后，車輛m滿足載重約束，將車輛m放入可行車輛集合X中。直至判斷完所有車輛。

(3) 以最小化行駛距離為目標(biāo)，并以車輛載重為約束條件，將剩余客戶點(diǎn)插入所有車輛(M輛)的所有位置，找出最佳插入點(diǎn)，將客戶插入，對(duì)剩余客戶按順序逐個(gè)執(zhí)行步驟3。若剩余客戶均分配完，輸出修復(fù)后的可行解，執(zhí)行步驟7；否則，執(zhí)行步驟4。

(4) 選出M輛車中剩余容量最大的車m1，將剩余客戶插入車m1的 最后一個(gè)位置，同時(shí)，從車m1中選出某一客戶點(diǎn)插入車m2(m2=1,2,···,M,m1≠m2)的最后一個(gè)位置，使兩車均滿足容量約束。對(duì)剩余客戶按順序逐個(gè)執(zhí)行步驟4。

(5) 若剩余客戶點(diǎn)分配完，輸出修復(fù)后的可行解，執(zhí)行步驟(7)；否則，執(zhí)行步驟(6)。

(6) 將剩余客戶點(diǎn)插入車輛m(m=1,2,···,M)的最后一個(gè)位置，同時(shí)，從車輛m選出某一客戶插入除車輛m的其他車輛的最后一個(gè)位置，使兩車均滿足容量約束。如果滿足，則插入；否則，直至判斷完所有車輛。對(duì)剩余客戶點(diǎn)按順序逐個(gè)執(zhí)行步驟(6)。

(7) 程序結(jié)束。

2.2.2 鄰域搜索算法

本小節(jié)針對(duì)上述獲得的可行解，設(shè)計(jì)了車輛間和車輛內(nèi)兩類共4種鄰域搜索操作，以進(jìn)一步提高可行解的質(zhì)量。算法步驟如下：

(1) 設(shè)可行解為 Π，令當(dāng)前最優(yōu)解Π?=Π，局部搜索次數(shù)為m axs， 車輛間搜索次數(shù)為vout，車輛內(nèi)搜索次數(shù)為vin， 車輛數(shù)為M，循環(huán)參數(shù)：l oop=1、Nout=1、Nin=1、Nv=1， 標(biāo)志位 S I=1。

(2) 在 Π中隨機(jī)選取兩輛車，在兩輛車中分別隨機(jī)選取一個(gè)客戶，得到客戶r和t，將t客戶插到r客戶之前，若不滿足載重約束條件，令 Π =Π?，執(zhí)行步驟(3)；反之，若獲得更優(yōu)解 Π′， 則令Π?=Π′，Π =Π′，執(zhí)行步驟(4)，若沒獲得更優(yōu)解，執(zhí)行步驟(3)。

(3) 在 Π中隨機(jī)選取兩輛車，在兩輛車中分別隨機(jī)選取一個(gè)客戶進(jìn)行交換，若不滿足載重約束條件，令Π =Π?，執(zhí)行步驟(4)；反之，若獲得更優(yōu)解Π′，則令Π?=Π′， Π =Π′，執(zhí)行步驟(4)。

2.3 次梯度算法

針對(duì)上述所提的拉格朗日對(duì)偶問題，本文采用次梯度算法進(jìn)行求解。其求解基本思想是：根據(jù)已知條件，求解對(duì)偶問題，獲得問題的下界和決策變量值，進(jìn)而計(jì)算更新次梯度以及拉格朗日乘子，通過不斷更新迭代求解對(duì)偶問題，直到找到滿足條件的對(duì)偶值。用次梯度算法求解對(duì)偶問題時(shí)，拉格朗日乘子λm(m∈C)更新公式為

2.4 改進(jìn)拉格朗日松弛算法步驟

整個(gè)求解算法步驟如下：

(1) 初始化拉格朗日乘子 λ0=0 ，迭代次數(shù)n=0，設(shè)初始上界U B為原問題較大的可行目標(biāo)值。

(2) 求解對(duì)偶問題(17)，獲得決策變量xn及下界LBn。

(3) 運(yùn)用2.2.1中的修復(fù)算法對(duì)下界進(jìn)行修復(fù)，若能修復(fù)得可行解 Π，執(zhí)行步驟(4)；否則上界U B保持不變，執(zhí)行步驟(5)。

(4) 運(yùn)用2.2.2中的二階段算法對(duì)可行解 Π進(jìn)行優(yōu)化，獲得優(yōu)化后的解 Π?，如果優(yōu)化后的目標(biāo)值Z(Π?)

(5) 運(yùn)用公式(18)、(19)和(20)更新次梯度、步長(zhǎng)和拉格朗日乘子。

(6) 判斷是否滿足終止條件，若滿足，輸出最佳上界 UB 和下界L Bn，程序結(jié)束；否則繼續(xù)執(zhí)行步驟(2)～(5)。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)環(huán)境及數(shù)據(jù)

3.2 結(jié)果分析

本文提出一種改進(jìn)拉格朗日松弛算法求解GCVRP，與商業(yè)求解器Gurobi所求解進(jìn)行了對(duì)比，求解結(jié)果見表3，從表可以看出：針對(duì)編號(hào)1～7的小規(guī)模算例，Gurobi求解器在1 800 s內(nèi)能求得4個(gè)算例的最優(yōu)解。小規(guī)模算例所求解平均G ap為5.78%，求解平均時(shí)間為825.277 s。改進(jìn)拉格朗日松弛算法在合理時(shí)間內(nèi)求得了問題的滿意解，求解平均G ap為7.17%，求解平均時(shí)間為219.896 s。這也體現(xiàn)Gurobi求解器在求解小規(guī)模問題的優(yōu)勢(shì)，能大概率求解問題的最優(yōu)解，但求解時(shí)間較長(zhǎng)。改進(jìn)拉格朗日松弛算法求解效率較高，在較短時(shí)間內(nèi)能求得問題的上下界，能求得問題的滿意解。

表3 不同規(guī)模問題算法對(duì)比1)Table 3 Comparison of algorithms for different scale problems

針對(duì)編號(hào)為8～19的中大規(guī)模算例，隨著問題規(guī)模的變大，解空間更加復(fù)雜，Gurobi求解器弊端凸顯，12個(gè)算例的平均G ap為21.12%，平均求解時(shí)間為5 100 s，求解質(zhì)量和求解效率不高。而改進(jìn)拉格朗日松弛算法通過松弛難約束，使問題變得容易求解，求解效率更高，求解質(zhì)量令人滿意，所求12個(gè)算例的平均G ap為7.87%，平均求解時(shí)間為2 422.016 s，求解質(zhì)量遠(yuǎn)優(yōu)于Gurobi求解器。

綜合來(lái)看，針對(duì)19個(gè)不同規(guī)模的算例，改進(jìn)拉格朗日松弛算法能較好地求解GCVRP，在合理時(shí)間內(nèi)，所有算例的平均G ap為7.61%。其在求得問題上界的同時(shí)，也能為原問題提供一個(gè)下界，可以衡量所求解的質(zhì)量，綜合了傳統(tǒng)拉格朗日松弛算法與啟發(fā)式算法的優(yōu)勢(shì)，是求解GCVRP的有效算法。

4 結(jié)論

本文在VRP的基礎(chǔ)上，進(jìn)一步考慮車輛載重和綠色因素，建立了GCVRP的混合整數(shù)規(guī)劃模型，并以總費(fèi)用最小為求解目標(biāo)，提出一種改進(jìn)拉格朗日松弛算法進(jìn)行求解。首先，采用次梯度法求解對(duì)偶問題，得到原問題的下界；其次利用修復(fù)算法和鄰域操作對(duì)下界進(jìn)行修復(fù)及優(yōu)化，得到原問題高質(zhì)量下界，實(shí)現(xiàn)上界下界并行求解。實(shí)驗(yàn)結(jié)果表明，本文算法的綜合求解性能優(yōu)于Gurobi，是一種求解GCVRP的有效方法。后續(xù)將進(jìn)一步考慮多車場(chǎng)多車型因素，并設(shè)計(jì)有效求解算法。