基于集合劃分的車輛路徑優(yōu)化精確算法研究

王維杰

（武漢理工大學(xué) 物流工程學(xué)院，湖北 武漢 430063）

1 引言

車輛路徑優(yōu)化（Vehicle Routing Problem）是運(yùn)輸問題的基礎(chǔ)模型，絕大多數(shù)實(shí)際應(yīng)用中的場(chǎng)景都是基于該模型的某種變形。Dantzig和Ramser在1959年首次研究了汽油配送卡車的最優(yōu)車輛路線問題[l]。Clarke和Wright于1964年提出了一個(gè)有效的貪婪啟發(fā)式算法[2]。這很快引起運(yùn)籌學(xué)、管理學(xué)、計(jì)算機(jī)科學(xué)、圖論等學(xué)科專家學(xué)者的高度重視。該問題的求解效率直接影響實(shí)際作業(yè)的效率，車輛配送路徑的優(yōu)化對(duì)于降低物流成本、減少污染和緩解交通壓力等方面均具有重要的社會(huì)意義。

在VRP中，按已知信息的特征將VRP分為確定性VRP和非確定性VRP[3]。在VRP剛提出時(shí)，研究的結(jié)點(diǎn)數(shù)的規(guī)模比較小，故研究學(xué)者對(duì)精確算法的研究比較多，各種精確算法紛紛涌現(xiàn)，整數(shù)規(guī)劃的應(yīng)用則更為廣泛，其主要的求解方法是割平面法、分支定界法以及列生成法。George Dantzig和Phil Wolfe使列生成法適合于具有可分解結(jié)構(gòu)的線性規(guī)劃問題并提出了 Dantzig-Wolfe分解。Desrochers，Desrosiers，Solomon則在D-W分解的基礎(chǔ)上最早應(yīng)用列生成算法并提出了Solomon基準(zhǔn)測(cè)試包[3]。本文主要通過研究VRPTW的子問題—帶資源約束的基本最短路徑問題（Elementary Shortest Path Problem with Resource Constraints，ESPPRC）來對(duì)算法進(jìn)行改善，該問題的松弛版本就是帶資源約束的最短路徑問題（Shortest Path Problem with Resource Constraints，SPPRC），兩者唯一不同的地方就在于是否允許重復(fù)的點(diǎn)出現(xiàn)在路徑中。大多數(shù)子問題的修改后行駛成本都是負(fù)數(shù)，允許最短路徑中包含回路（重復(fù)的點(diǎn)形成的循環(huán)路徑）能夠使其求得更小成本的路徑。同時(shí)在動(dòng)態(tài)規(guī)劃算法的框架下，利用帕累托最優(yōu)的原則能夠刪除更多的冗余路徑。因此，SPPRC的應(yīng)用非常成功，越來越多的學(xué)者選擇研究SPPRC而避開ESPPRC。但是該方法有個(gè)顯著的缺點(diǎn)，它產(chǎn)生的下界質(zhì)量難以保證，用ESPPRC來描述子問題則能夠保證更好的下界。

2 基于集合劃分的數(shù)學(xué)模型

帶時(shí)間窗的車輛路徑問題（Vehicle Routing Problem with Time Window）可以通過D-W分解劃分為主問題為集合劃分問題（Set-partitioning problem）及子問題為帶資源約束的基本最短路徑問題的模型。用P表示所有可能的路徑集合（其中的所有路徑均滿足時(shí)間窗與容量約束），對(duì)于每條路徑p∈P，cp表示路徑p的行駛成本，對(duì)于每個(gè)客戶點(diǎn)i，當(dāng)路徑p訪問了客戶點(diǎn)i則aip取值為1，否則取值為0。定義路徑變量yp，代表該路徑是否在最優(yōu)路徑上，如果路徑p在最優(yōu)路徑中則yp=1，若不在則取0。VRPTW的主問題數(shù)學(xué)模型如下：

目標(biāo)函數(shù)：

約束條件：

對(duì)約束條件（3）進(jìn)行線性松弛，使yp∈R，得到：

將約束條件（4）替代（3）加入到主問題模型中，就得到了線性主問題的線性模型。此時(shí)，變量yp代表每條路徑p的權(quán)重，可以得到限制主問題（restricted master problem，RMP）的數(shù)學(xué)模型如下：

目標(biāo)函數(shù)：

約束條件：

集合P′?P，僅包含所有已經(jīng)生成的路徑。限制主問題（RMP）的初始解選擇所有只含有一個(gè)客戶點(diǎn)的路徑，即路徑0→i→n+1（i∈C），此時(shí)約束條件（2）等式左側(cè)的矩陣為單位矩陣。

子問題為帶資源約束的基本最短路徑問題（ESPPRC），如何產(chǎn)生合適的路徑p是解決該問題的關(guān)鍵，其常見的求解模型則是基于整數(shù)線性規(guī)劃。定義是邊（i，j）上修改后的成本代表客戶點(diǎn)i上的服務(wù)開始時(shí)間。

子問題數(shù)學(xué)模型：

目標(biāo)函數(shù)（8）表示總行駛路程最短，即總成本（修改后）最小。約束條件（9）是容量約束；約束條件（13）和（14）是時(shí)間窗約束，Mij通常代表一個(gè)非常大的常數(shù) max{bi-aj+sti+tij，0},(i,j)∈A 。約束條件（10）、（11）和（12）是流（flow）約束，保證了一條不含循環(huán)的路徑從倉(cāng)庫(kù)點(diǎn)0到倉(cāng)庫(kù)點(diǎn)n+1。約束條件（15）定義了一個(gè)0-1變量。作為VRPTW的子問題，cij是一個(gè)非負(fù)實(shí)數(shù)，而可以是任意實(shí)數(shù)（在大多數(shù)情況下都為負(fù)數(shù)）。假設(shè)λ和π分別是當(dāng)前RMP的原始最優(yōu)解和對(duì)偶最優(yōu)解，給定一個(gè)集合A，它由列aj,j∈J組成，目標(biāo)函數(shù)的成本系數(shù)cj可以通過關(guān)于aj的函數(shù)式c計(jì)算得到，由此得到子問題的修改成本：

如果-c*≥0，且不存在負(fù)的-cj,j∈J，那么就說限制主問題的解λ就是主問題的最優(yōu)解。另一方面，向RMP添加的列都來自子問題的最優(yōu)解，并且不斷對(duì)RMP進(jìn)行重新優(yōu)化。為了提升子問題的求解效率，將以下三組不等式加入子問題模型用以減少子回路。

2.1 點(diǎn)邊不等式

定義集合M?N，給出不等式：

M集合僅包含同時(shí)滿足以下兩個(gè)條件的點(diǎn)k。i→j→k是可行路徑， 即在時(shí)間窗和容量約束下車輛能夠按順序依次訪問點(diǎn)i，j，k；邊（i，j）、（j，k）與（i，k）的修改后成本滿足+≤。

2.2 時(shí)間前后不等式

2.3 順序前后不等式

對(duì)任意客戶點(diǎn)i，定義兩組集合Pre,Suc?N，描述如下：

Pre代表一些只能出現(xiàn)在i點(diǎn)之前的點(diǎn)，Suc代表一些只能出現(xiàn)在i點(diǎn)之后的點(diǎn)。考慮所有從集合Pre中出發(fā)到集合Suc的路徑，故順序前后不等式如下：

在利用列生成算法解決基于集合劃分的數(shù)學(xué)模型時(shí)，會(huì)得到線性主問題的下界（檢驗(yàn)成本為0時(shí)），二維車流模型相比于經(jīng)典的三維多商品網(wǎng)絡(luò)流模型，整數(shù)變量x減少了一維，我們用該模型替代列生成算法中得到主問題線性最優(yōu)解后的分支-切割過程。其中定義K是所需車輛的數(shù)量，即路徑數(shù)量，模型如下：

目標(biāo)函數(shù)（20）同樣是求總行駛路程最短。入度（indegree）和出度（outdegree）約束條件（21）和（22）保證每個(gè)客戶點(diǎn)進(jìn)入的邊數(shù)與離開的邊數(shù)均為1；（23）與（24）保證了從倉(cāng)庫(kù)點(diǎn)0出發(fā)的車輛數(shù)與到達(dá)倉(cāng)庫(kù)點(diǎn)n+1的車輛數(shù)就是圖中的路徑數(shù)量；（25）和（26）確保了容量約束的滿足；（27）與（28）確保了時(shí)間窗約束的滿足；（29）定義了一個(gè)0-1變量。

3 兩點(diǎn)割集不等式

定義一個(gè)多面體（polyhedron）P?Rn，代表一組滿足有限個(gè)線性不等式的點(diǎn)的集合，即P={x?Rn:Ax≤b}，其中(A,b)是一個(gè)m×(n+1)階的矩陣，R表示所有實(shí)數(shù)集合。為了避免額外的時(shí)間變量s對(duì)模型多面體結(jié)構(gòu)的影響，Ascheuer在2000年提出了不可行路徑約束條件，并將它應(yīng)用在帶時(shí)間窗的非對(duì)稱旅行商問題（Asymmetric Traveling Salesman Problem with Time Windows，ATSPTW）中，用來替代傳統(tǒng)的時(shí)間窗約束與容量約束，效果很好[4]。2007年，Kallehauge把這項(xiàng)技術(shù)應(yīng)用在VRPTW上，首次探討VRPTW的多面體結(jié)構(gòu)[5]。接下來我們對(duì)Kallehauge的數(shù)學(xué)模型進(jìn)行描述，定義路徑P是邊的集合{(vi,vi+1)|i=1,…,k-1}，也可以表示為P=(v1,v2,…vk)，其中|P|=k-1，若i≠j則有vi≠vj。最小不可行路徑 P={(vi,vi+1)∈A|i=1,…,k-1}，即P{(v1,v2)}和P{(vp-1,vp)}均為可行路徑。定義在圖G中的所有最小不可行路徑的集合為PG。定義xij是0-1變量，如果邊(i,j)∈A在最優(yōu)路徑中則xij=1，若不在則取0。數(shù)學(xué)模型如下：

目標(biāo)函數(shù)：

約束條件：

目標(biāo)函數(shù)（30）是求圖中總路線長(zhǎng)度最短。約束條件（31）與（32）保證了每個(gè)客戶點(diǎn)剛好訪問一次；約束條件（33）是子回路不等式，確保了圖中的路線不包含子回路；約束條件（34）是路徑不等式，保證了所有路徑都滿足時(shí)間窗和容量約束。為了得到更易進(jìn)行多面體理論分析的模型，選取約束條件（10）、（11）、（12）與基于路徑不等式模型中的約束條件（33）、（34）進(jìn)行混合，得到新的模型：定義路徑P是邊（arc）的集合{(vi,vi+1)|i=1,…,k-1}，也可以表示為P=(v1,v2,…vk)，其中|P|=k-1，若i≠j則有vi≠vj。

目標(biāo)函數(shù)：

目標(biāo)函數(shù)與約束條件的含義與前文描述的一致，由此可以定義ESPPRC的多胞形（polytope）為：

經(jīng)典的Dantzig-Fulkerson-Johnson子回路刪除（subtour elimination）約束可以等價(jià)地用一組廣義割集不等式（generalized cutset inequalities，GCS）來表示：

GCS（43）對(duì)ESPPRC的多胞形（polytope）是有效的。

定義滿足以下形式的可行路徑集合P={(0,i,j,k)i,j∈C,k∈N}，對(duì)于任意一條可行路徑p∈P，兩點(diǎn)加強(qiáng)GCS不等式:

兩點(diǎn)情況是加強(qiáng)割集不等式的最小應(yīng)用示例，如果將它們?nèi)考尤胱訂栴}的數(shù)學(xué)模型中，會(huì)增加N（N-1）個(gè)不等式，數(shù)量并不多，但由于已直接將三組有效不等式直接加入子問題數(shù)學(xué)模型中，本次添加更會(huì)導(dǎo)致部分已滿足該條件的節(jié)點(diǎn)求解時(shí)間變長(zhǎng)，降低整體模型及算法的效率。本文將采用割平面回調(diào)形式加入兩點(diǎn)加強(qiáng)割集不等式，即不直接將不等式加入至子問題的數(shù)學(xué)模型，而是將兩點(diǎn)加強(qiáng)割集不等式加在根節(jié)點(diǎn)上，通過在每個(gè)根節(jié)點(diǎn)處檢查每個(gè)兩點(diǎn)加強(qiáng)割集不等式，查找調(diào)用該回調(diào)的根節(jié)點(diǎn)當(dāng)前線性松弛解所違反的不等式，根據(jù)該節(jié)點(diǎn)檢查的順序逐次將所違反的不等式作為割平面加入到當(dāng)前根節(jié)點(diǎn)中。

4 實(shí)驗(yàn)結(jié)果

將CVRP中的二維車流模型擴(kuò)展至VRPTW中，用它來替代列生成算法中得到主問題線性最優(yōu)解后的分支-切割過程以及加入子問題中的三組不等式已被證明有效[6]。為了進(jìn)一步檢驗(yàn)兩點(diǎn)加強(qiáng)割集不等式的性能，本文將以未加入割集不等式的實(shí)驗(yàn)組作為對(duì)照組對(duì)VRPTW進(jìn)行實(shí)驗(yàn)：

主問題：基于集合劃分的數(shù)學(xué)模型；

子問題：加入兩點(diǎn)加強(qiáng)割集不等式的整數(shù)線性規(guī)劃模型。

實(shí)驗(yàn)數(shù)據(jù)來自于經(jīng)典的Solomon基準(zhǔn)測(cè)試包，它是目前求解VRPTW各類算法最常用的標(biāo)準(zhǔn)測(cè)試集[3]。測(cè)試包中的數(shù)據(jù)集合分為三類：

r組：客戶點(diǎn)坐標(biāo)呈隨機(jī)分布；

c組：客戶點(diǎn)坐標(biāo)呈現(xiàn)分塊聚集分布；

rc組：客戶點(diǎn)坐標(biāo)既有隨機(jī)分布又有分塊聚集分布。

每個(gè)算例都包括100個(gè)客戶點(diǎn)、50個(gè)客戶點(diǎn)和25個(gè)客戶點(diǎn)三種規(guī)模，目前Solomon基準(zhǔn)測(cè)試包對(duì)于我們的算法過于困難，僅選擇其中25個(gè)客戶點(diǎn)的數(shù)據(jù)進(jìn)行測(cè)試，以驗(yàn)證修改后的二維車流模型的性能。

算法主要思路如下：通過D-W分解將VRPTW劃分為主問題為集合劃分以及子問題為帶資源約束的基本最短路徑問題。針對(duì)主問題：通過對(duì)其進(jìn)行線性松弛得到原始的限制主問題。接著將僅含有一個(gè)客戶點(diǎn)的初始解帶入限制主問題后，對(duì)其用單純形法求解，同時(shí)生成每個(gè)客戶點(diǎn)的對(duì)偶變量（π）并求得其值。針對(duì)子問題：根據(jù)主問題所生成的對(duì)偶變量對(duì)目標(biāo)函數(shù)中的路線行駛成本進(jìn)行修改，通過割平面回調(diào)形式加入兩點(diǎn)割集不等式，求解子問題。子問題的解值稱為檢驗(yàn)成本，解可以構(gòu)成一條新的路徑。如果檢驗(yàn)成本為負(fù)數(shù)，則將新的路徑作為一列加入到限制主問題當(dāng)中重新求解，生成新的對(duì)偶變量，如此循環(huán)直到檢驗(yàn)成本非負(fù)。最后對(duì)主問題進(jìn)行求解，若可以求得整數(shù)解，則該解為最優(yōu)解；若求得為非整數(shù)解，則該值為此問題的下界，不帶入下界用二維車流模型來進(jìn)行求解。

得到主問題的下界后用基于二維車流的數(shù)學(xué)模型求最優(yōu)解，該方法簡(jiǎn)稱為ILP算法。ILP算法是通過在Ilog Cplex中的Cplex優(yōu)化器上建立列生成結(jié)構(gòu)與數(shù)學(xué)模型進(jìn)行求解，Cplex參數(shù)設(shè)置為默認(rèn)。計(jì)算實(shí)驗(yàn)都是在含有2.5GHz的Intel（i7-4710MQ）處理器、12GB內(nèi)存以及配有64位Windows操作系統(tǒng)的戴爾一體機(jī)上運(yùn)行。使用Python（編譯版本2.7）進(jìn)行算法編譯，借助IBM Ilog Cplex/Concert 12.6進(jìn)行建模優(yōu)化。具體指標(biāo)定義如下：

LB：下界（Lower Bound），空白格代表500s內(nèi)未求解出下界；

Num：子問題的迭代次數(shù)，即求解了Num個(gè)子問題；

Paths：產(chǎn)生的路徑數(shù)量之和，即加入到主問題中列的數(shù)量；

LBT：求解到下界的時(shí)間，單位是s；

OptT：修改后二維車流模型的求解時(shí)間，單位是s；

TT：總時(shí)間（Total Time）。

表1 原VRPTW實(shí)驗(yàn)結(jié)果（25個(gè)點(diǎn)窄時(shí)間窗算例）

從表1、表2對(duì)比可以看出以割平面回調(diào)形式加入兩點(diǎn)割集不等式后模型的求解時(shí)間全面優(yōu)于原始模型，在均能求出最優(yōu)解的22個(gè)算例中，改進(jìn)模型的平均運(yùn)行時(shí)間均小于未改進(jìn)模型，可以認(rèn)為兩點(diǎn)加強(qiáng)割集不等式的加入有效提升了模型的運(yùn)算效率，具體數(shù)據(jù)見表3。

表2 加入兩點(diǎn)GCS后VRPTW實(shí)驗(yàn)結(jié)果（25個(gè)點(diǎn)窄時(shí)間窗算例）表3 VRPTW 平均時(shí)間對(duì)比表

表2 加入兩點(diǎn)GCS后VRPTW實(shí)驗(yàn)結(jié)果（25個(gè)點(diǎn)窄時(shí)間窗算例）表3 VRPTW 平均時(shí)間對(duì)比表

5 結(jié)束語(yǔ)

帶時(shí)間窗的車輛路徑優(yōu)化是現(xiàn)實(shí)生活中最常見的運(yùn)輸問題之一，具有很高的研究?jī)r(jià)值。但國(guó)內(nèi)外學(xué)者大都將研究方向集中于啟發(fā)式算法，因其通用性及較強(qiáng)的可擴(kuò)展性，研究方向更加豐富。但某種程度上也導(dǎo)致了啟發(fā)式算法的精確性不足。而針對(duì)組合優(yōu)化問題，若無完整的理論分析及精確算法作為基礎(chǔ)，就很難建立適合的啟發(fā)式算法，更不能保證其求解質(zhì)量。我們建立了基于集合劃分的數(shù)學(xué)模型，證明了兩點(diǎn)加強(qiáng)割集不等式對(duì)于子問題的多胞形是有效的，并通過割平面回調(diào)的形式對(duì)子問題加入兩點(diǎn)加強(qiáng)割集不等式，有效提升了算法的求解速度。