答家瑞，鄭瀾波

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

帶時(shí)間窗的車(chē)輛路徑問(wèn)題的精確算法研究

答家瑞，鄭瀾波

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

將CVRP（Capacitated Vehicle Routing Problem）中的二維車(chē)流模型擴(kuò)展至VRPTW中，用它來(lái)替代列生成算法中的分支-切割過(guò)程，為解決VRPTW提供了一種新思路。同時(shí)對(duì)最少車(chē)輛數(shù)量的理論上界進(jìn)行了猜想，并用Solomon基準(zhǔn)測(cè)試包進(jìn)行了實(shí)驗(yàn)，求解出的算例均肯定了這一猜想。

時(shí)間窗；車(chē)輛路徑問(wèn)題；運(yùn)籌學(xué)；整數(shù)線(xiàn)性規(guī)劃；列生成；精確算法

1 引言

帶時(shí)間窗的車(chē)輛路徑問(wèn)題（Vehicle Routing Problem with Time Windows，VRPTW）是對(duì)車(chē)輛路徑問(wèn)題的擴(kuò)展，最早由Pullen在1967年提出[1]，其主要不同點(diǎn)在于添加了一個(gè)時(shí)間窗來(lái)限制每個(gè)客戶(hù)地點(diǎn)的服務(wù)時(shí)間，該限制可以針對(duì)于客戶(hù)服務(wù)的最早開(kāi)始時(shí)間、車(chē)輛停留在客戶(hù)點(diǎn)的最短時(shí)間以及客戶(hù)服務(wù)的最晚開(kāi)始時(shí)間等，這些限制統(tǒng)稱(chēng)為時(shí)間窗約束。

VRPTW可以定義為一組車(chē)輛(V)，一些客戶(hù)(C)和一個(gè)有向圖(G)。該有向圖中包括|C|+2個(gè)點(diǎn)，其中客戶(hù)點(diǎn)用1，2，…，n表示，倉(cāng)庫(kù)點(diǎn)用點(diǎn)0（開(kāi)始點(diǎn)）和n+1（結(jié)束點(diǎn)）表示。用N來(lái)表示所有點(diǎn)的集合，A表示所有邊的集合，即G(N,A) 。對(duì)于每條邊（i，j），cij表示該邊上的行駛成本，tij表示途經(jīng)該邊的行駛時(shí)間。對(duì)于每個(gè)客戶(hù)點(diǎn)i都有一個(gè)時(shí)間窗[ai,bi]，即到達(dá)i點(diǎn)的車(chē)輛最早開(kāi)始服務(wù)時(shí)間為ai，最晚開(kāi)始服務(wù)時(shí)間為bi，每個(gè)點(diǎn)的服務(wù)時(shí)間長(zhǎng)度為sti，需要裝載的貨物重量（容量）為di。對(duì)于每輛車(chē)都有一個(gè)相同的載重量（最大容量）q。

2 列生成算法

列生成（Column generation）[2]是一種求解大規(guī)模整數(shù)線(xiàn)性規(guī)劃的方法。使用它進(jìn)行求解時(shí)，不需要用到原始問(wèn)題中所有的變量，而是求解一個(gè)子集（subset）。具體來(lái)說(shuō)列生成通常會(huì)產(chǎn)生一個(gè)主問(wèn)題，這個(gè)主問(wèn)題的解依賴(lài)于一個(gè)或多個(gè)子問(wèn)題的多重解（repeated solution）。

接下來(lái)用數(shù)學(xué)的方法來(lái)描述，我們用線(xiàn)性規(guī)劃的方法定義主問(wèn)題（master problem，MP），數(shù)學(xué)模型如下：

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

約束條件：

這是一個(gè)通用的線(xiàn)性規(guī)劃模型，其中集合J代表所有的列（column）。MP的對(duì)偶變量用非負(fù)向量π來(lái)表示，根據(jù)對(duì)偶理論，我們希望找到一個(gè) j∈J能夠使得：=cj-πtaj的值最小。當(dāng) ||J的值很大時(shí)，這種直接定價(jià)（pricing）的方法運(yùn)算成本很高。所以我們考慮取一個(gè)合理的子集 j'∈J來(lái)操作，用間接枚舉（implicit enumeration）的方法來(lái)計(jì)算檢驗(yàn)成本（reduced cost），它被稱(chēng)為限制主問(wèn)題（restricted master problem，RMP）。假設(shè)λ和π分別是當(dāng)前RMP的原始最優(yōu)解和對(duì)偶最優(yōu)解，給定一個(gè)集合A，它由列（column）aj,j∈J組成，目標(biāo)函數(shù)的成本系數(shù)cj可以通過(guò)關(guān)于aj的函數(shù)式c計(jì)算得到，由此得到子問(wèn)題的修改成本：

帶時(shí)間窗的車(chē)輛路徑問(wèn)題（VRPTW）可以用集合劃分問(wèn)題（Set-partitioning problem）的模型來(lái)構(gòu)造[3]。用P來(lái)表示所有可能的路徑集合（其中的所有路徑均滿(mǎn)足時(shí)間窗與容量約束），對(duì)于每條路徑p∈P，cp表示路徑p的行駛成本，對(duì)于每個(gè)客戶(hù)點(diǎn)i，當(dāng)路徑p訪(fǎng)問(wèn)了客戶(hù)點(diǎn)i則aip取值為1，否則取值為0。定義路徑變量Yp，代表該路徑是否在最優(yōu)路徑上，如果路徑p在最優(yōu)路徑中則Yp=1，若不在則取0。數(shù)學(xué)模型如下：

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

約束條件：

目標(biāo)函數(shù)（5）所有路徑的總成本最小。集合劃分約束（6）保證了每個(gè)客戶(hù)點(diǎn)剛好被一輛車(chē)訪(fǎng)問(wèn)一次。

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

將約束條件（8）替代（7）加入到主問(wèn)題模型中，就得到了線(xiàn)性主問(wèn)題的線(xiàn)性模型。此時(shí)，變量Yp代表每條路徑p的權(quán)重，約束條件（6）等式左側(cè)用矩陣的方式描述可以表示為：每一列代表一條可行路徑乘以權(quán)重Yp，每一行代表一個(gè)客戶(hù)點(diǎn)乘以權(quán)重Yp。

基于列生成的理論，可以得到限制主問(wèn)題的數(shù)學(xué)模型：

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

約束條件：

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

算法主要思路如下：首先對(duì)主問(wèn)題進(jìn)行線(xiàn)性松弛，得到原始的限制主問(wèn)題。接著將初始解放入限制主問(wèn)題后，對(duì)其用單純形法求解，同時(shí)生成每個(gè)客戶(hù)點(diǎn)的對(duì)偶變量（π）并求得其值。接著對(duì)子問(wèn)題目標(biāo)函數(shù)中的路線(xiàn)行駛成本進(jìn)行修改是邊（i，j）上修改后的成本，=cij-(πi+πj)/2，并求解子問(wèn)題，子問(wèn)題的解值稱(chēng)為檢驗(yàn)成本（reduced cost），解（solution）可以構(gòu)成一條新的路徑。如果檢驗(yàn)成本為負(fù)數(shù)（negative reduced cost），則將新的路徑作為一列（column）加入到限制主問(wèn)題（restricted master problem）當(dāng)中，重新求解，生成新的對(duì)偶變量，如此循環(huán)直到檢驗(yàn)成本非負(fù)。最后對(duì)主問(wèn)題進(jìn)行求解，若可以求得整數(shù)解，則該解為最優(yōu)解；若求得不為整數(shù)解，則該值為此問(wèn)題的下界，接著用分支-切割法進(jìn)行求解，重復(fù)上述循環(huán)直到求得最優(yōu)解。

3 基于二維車(chē)流的數(shù)學(xué)模型

傳統(tǒng)的帶容量約束的車(chē)輛路徑問(wèn)題（Capacitated VRP，CVRP）的數(shù)學(xué)模型十分豐富，但由于時(shí)間窗的約束，能夠應(yīng)用到VRPTW的數(shù)學(xué)模型很少，主要是因?yàn)闀r(shí)間窗的約束難以和傳統(tǒng)的容量-子回路消除約束相結(jié)合。我們嘗試著將CVRP的經(jīng)典二維車(chē)流模型加以改進(jìn)，演變?yōu)檫m用于VRPTW的數(shù)學(xué)模型。

首先介紹傳統(tǒng)的容量-子回路消除約束，以下描述中出現(xiàn)的變量與參數(shù)均與前文敘述一致。給定一個(gè)集合W?C，定義r(W)代表W集合中能夠服務(wù)所有客戶(hù)點(diǎn)并滿(mǎn)足車(chē)輛容量約束的最小車(chē)輛數(shù)，定義，代表W集合中的所有貨物總重量。根據(jù)裝箱問(wèn)題（Bin Packing Problem，BPP）的理論研究，r(W)的上界（upper bound）為d(W)/q（向上取整），憑借對(duì)CVRP整數(shù)規(guī)劃模型的多面體結(jié)構(gòu)分析，學(xué)者們將該上界提高到2d(W)/q（向上取整）?；诖?，引出兩個(gè)重要的約束條件，容量切割約束（Capacity-Cut Constraints，CCCs）：

與通用子回路消除約束（Generalized Subtour Elimination Constraints,GSECs）：

可以看出這類(lèi)約束在模型中的個(gè)數(shù)是指數(shù)級(jí)的，需要用分離（separation）的算法來(lái)尋找違反約束條件的集合，加入時(shí)間窗約束后，分離問(wèn)題變得更加復(fù)雜，難以快速得到所需子集。由于理論上界研究的突破，應(yīng)用這類(lèi)約束處理CVRP十分有效，但該下界無(wú)法應(yīng)用在VRPTW中。

在VRPTW中得出類(lèi)似r(W)的下界，即W集合中能夠服務(wù)所有客戶(hù)點(diǎn)并滿(mǎn)足車(chē)輛容量約束和時(shí)間窗約束的最小車(chē)輛數(shù)目十分困難，所以我們考慮使用1960年Miller，Tucker和Zemlin對(duì)非對(duì)稱(chēng)旅行商問(wèn)題（ATSP）進(jìn)行子回路消除的思路進(jìn)行處理（簡(jiǎn)稱(chēng)MTZ）。引入一個(gè)額外變量ui∈R,i∈C，MTZ的約束可表述為：

不難看出經(jīng)典的三維多商品網(wǎng)絡(luò)流模型中的時(shí)間窗約束就是由MTZ演變而來(lái)。1990年Desrochers和Laporte對(duì)MTZ進(jìn)行了增強(qiáng)，提出了新的約束條件（簡(jiǎn)稱(chēng)DL）：

這種構(gòu)造保證了當(dāng)xij=1時(shí)，ui+1=uj，從而使得其多面體結(jié)構(gòu)更小。

我們將這種思路應(yīng)用在容量-子回路消除約束中，提出新的數(shù)學(xué)模型，其中定義K是所需車(chē)輛的數(shù)量，即路徑數(shù)量，模型如下：

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

約束條件：

目標(biāo)函數(shù)（5）同樣是求總行駛路程最短。入度（indegree）和出度（outdegree）約束條件（17）和（18）保證每個(gè)客戶(hù)點(diǎn)進(jìn)入的邊數(shù)與離開(kāi)的邊數(shù)均為1；（19）與（20）保證了從倉(cāng)庫(kù)點(diǎn)0出發(fā)的車(chē)輛數(shù)與到達(dá)倉(cāng)庫(kù)點(diǎn)n+1的車(chē)輛數(shù)就是圖中的路徑數(shù)量；（21）和（22）確保了容量約束的滿(mǎn)足；（23）與（24）確保了時(shí)間窗約束的滿(mǎn)足，其中Mij通常代表一個(gè)非常大的常數(shù)，在該問(wèn)題中它的值可以減少為

該模型相比于經(jīng)典的三維多商品網(wǎng)絡(luò)流模型[4]，整數(shù)變量x減少了一維，但增加了一個(gè)線(xiàn)性變量u。提升該模型計(jì)算性能的關(guān)鍵在于如何降低K的上界，即至少需要多少輛車(chē)可以形成最優(yōu)路徑。在利用列生成算法解決基于集合劃分的數(shù)學(xué)模型時(shí)，會(huì)得到線(xiàn)性主問(wèn)題的下界（檢驗(yàn)成本為0時(shí)），同時(shí)產(chǎn)生線(xiàn)性解中的路徑數(shù)量值K'（該數(shù)量值可能不是整數(shù)），我們認(rèn)為K'（向上取整）或K'+1（向上取整）很有可能就是K的理論上界，所以用該模型來(lái)替代列生成算法中得到主問(wèn)題線(xiàn)性最優(yōu)解后的分支-切割過(guò)程，以測(cè)試模型質(zhì)量與最優(yōu)性。

4 帶資源約束的基本最短路徑問(wèn)題

帶資源約束的基本最短路徑問(wèn)題（Elementary Shortest Path Problem with Resource Constraints,ESPPRC）是VRPTW基于列生成算法的子問(wèn)題。當(dāng)用集合劃分模型作為主問(wèn)題，利用列生成算法來(lái)解決VRPTW時(shí)，子問(wèn)題被分解成多個(gè)性質(zhì)相同的ESPPRC。更加明確地來(lái)說(shuō)，子問(wèn)題是一個(gè)帶時(shí)間窗與容量約束的基本最短路徑問(wèn)題（Elementary Shortest Path Problem with Time Windows and Capacity Constraints，ESPPTWCC），其中“基本”（elementary）的意思是每個(gè)客戶(hù)點(diǎn)只能出現(xiàn)一次，即最優(yōu)解中不存在回路（cycle）[5]。

本文基于整數(shù)線(xiàn)性規(guī)劃，加入以下三組有效不等式對(duì)子問(wèn)題進(jìn)行了建模求解。

4.1 點(diǎn)邊不等式

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

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

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

如果i點(diǎn)在時(shí)間窗約束下不能到達(dá)j點(diǎn)，即ai-bj+sti+tij＞0時(shí)，則i點(diǎn)只能在j點(diǎn)訪(fǎng)問(wèn)后再進(jìn)行訪(fǎng)問(wèn)或者不訪(fǎng)問(wèn)。當(dāng)ai-bj-stj-tji≤0時(shí)，時(shí)間前后不等式可以表述為：

4.3 順序前后不等式

對(duì)任意客戶(hù)點(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的路徑，如果點(diǎn)i在最短路徑上，則該路徑一定會(huì)經(jīng)過(guò)i點(diǎn)，故順序前后不等式如下：

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

實(shí)驗(yàn)數(shù)據(jù)來(lái)自于經(jīng)典的Solomon基準(zhǔn)測(cè)試包[3]，測(cè)試包中的數(shù)據(jù)集合分為三類(lèi)，

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

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

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

每個(gè)算例都包括100個(gè)客戶(hù)點(diǎn)、50個(gè)客戶(hù)點(diǎn)和25個(gè)客戶(hù)點(diǎn)三種規(guī)模，后兩種分別是取包括100個(gè)客戶(hù)點(diǎn)的完整算例中的前50個(gè)和25個(gè)客戶(hù)點(diǎn)。

計(jì)算實(shí)驗(yàn)都是在含有2.5GHz的Intel（i7-4710MQ）處理器、12GB內(nèi)存以及配有64位Windows操作系統(tǒng)的戴爾一體機(jī)上運(yùn)行。使用Java（編譯版本1.7）在Eclipse上進(jìn)行算法編譯，借助IBM Ilog Cplex/Concert 12.6進(jìn)行建模優(yōu)化。

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

子問(wèn)題：加入三種有效不等式后的整數(shù)線(xiàn)性規(guī)劃模型。

得到主問(wèn)題的下界后用基于二維車(chē)流的數(shù)學(xué)模型求最優(yōu)解，該方法簡(jiǎn)稱(chēng)為ILP算法。

ILP算法是通過(guò)在Ilog Cplex中的Cplex優(yōu)化器上建立列生成結(jié)構(gòu)與數(shù)學(xué)模型進(jìn)行求解，Cplex參數(shù)設(shè)置為默認(rèn)。

實(shí)驗(yàn)在求解完所有子問(wèn)題得到主問(wèn)題的下界后，會(huì)得到一個(gè)線(xiàn)性解中的路徑數(shù)量值K'，我們將二維車(chē)流模型中的K（所需的車(chē)輛數(shù)量或路徑數(shù)量）定義為整數(shù)變量，其取值范圍設(shè)置為[K ',K'+1]，直接求解該模型，輸出最優(yōu)解。為了測(cè)試修改后二維車(chē)流模型的性能，我們不會(huì)將下界加入到模型中去。

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

表中具體說(shuō)明如下：

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

Opt：最優(yōu)解（Opimal solution），空白格代表500秒內(nèi)未求解出最優(yōu)解；

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

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

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

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

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

表1 r組25個(gè)點(diǎn)的VRPTW測(cè)試結(jié)果

表2 rc組25個(gè)點(diǎn)的VRPTW測(cè)試結(jié)果

在表1、表2和表3中，可以看出在求解出下界后，用修改后的二維車(chē)流模型來(lái)求解得最優(yōu)解的時(shí)間（OptT）非常短，可以認(rèn)為它能夠替代列生成算法后期的分支-定界部分。在目前求解出的所有算例中，二維車(chē)流模型的解與前人提供的最優(yōu)解相同，但由于尚未求解出更多的算例，無(wú)法肯定K'+1（向上取整）就是K（最優(yōu)解中的路徑數(shù)量）的理論上界，目前的結(jié)果都傾向于肯定這一猜想。

表3 c組25個(gè)點(diǎn)的VRPTW測(cè)試結(jié)果

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

帶時(shí)間窗的車(chē)輛路徑問(wèn)題是經(jīng)典的組合優(yōu)化問(wèn)題，也是目前應(yīng)用最廣泛的運(yùn)輸問(wèn)題之一，具有很高的研究?jī)r(jià)值。國(guó)內(nèi)的研究往往傾向于用啟發(fā)式方法解決該問(wèn)題，精確算法的研究非常匱乏。對(duì)于啟發(fā)式方法，我們認(rèn)為它最大的優(yōu)勢(shì)在于高效性和通用性，而這也在一定程度上導(dǎo)致了它精確性的不足。針對(duì)一個(gè)NP-hard問(wèn)題，如果沒(méi)有透徹的理論分析和精確算法做支撐，就很難建立適合的啟發(fā)式算法，更不能保證其求解質(zhì)量。

我們將CVRP中的二維車(chē)流模型擴(kuò)展至VRPTW中，用它來(lái)替代列生成算法中得到主問(wèn)題線(xiàn)性最優(yōu)解后的分支-切割過(guò)程，為解決VRPTW提供了一種新思路。

[1]Pullen H G M,Webb M H J.A Computer Application to a Transport Scheduling Problem[J].Computer Journal,1967,10(1).

[2]Danzig G B,Wolfe P.The decomposition algorithm for linear programming[J].Econometrica,1961,(4)：767-778.

[3]Desrochers M,Desrosiers J,Solomon M.A new optimization algorithm for the vehicle routing problem with time windows[J].Operations research,1992,40(2)：342-354.

[4]Toth P,Vigo D.The vehicle routing problem[M].Beijing：Tsinghua University Press,2011.

[5]Jepsen M K,Petersen B,Spoorendonk S,et al.A branch-andcut algorithm for the capacitated profitable tour problem[J].Discrete Optimization,2014,14(C)：78-96.

Study on Exact Algorithm for Vehicle Routing Problem with Time Window

Da Jiarui,Zheng Lanbo
(School of Logistics Engineering,Wuhan University of Technology,Wuhan 430063,China)

In this paper,we extended the 2D car flow model from the capacitated vehicle routing problem(CVRP)to the vehicle routing problem with time window(VRPTW)to replace the branch-and-cut process of the column generation algorithm and provide a new line of thinking for the solution of the VRPTW.At the same time,we hypothesized the theoretical upper bound of the minimal vehicle quantity and tested and verified it using the Solomon benchmark test package.

time window;vehicle routing problem;operations research;integer linear programming;column generation;exact algorithm

F224.0

A

1005-152X(2017)06-0095-05

10.3969/j.issn.1005-152X.2017.06.023

2017-05-01

國(guó)家自然科學(xué)基金（71501152）

答家瑞（1992-），男，湖北武漢人，碩士研究生，研究方向：運(yùn)籌學(xué)與供應(yīng)鏈網(wǎng)絡(luò)；鄭瀾波（1980-），女，武漢理工大學(xué)副教授，博士，研究方向：運(yùn)籌學(xué)與供應(yīng)鏈網(wǎng)絡(luò)。