車輛路徑規(guī)劃問題的逆向優(yōu)化方法

陳禹伊， 陳 璐

(上海交通大學 機械與動力工程學院， 上海 200240)

隨著信息技術的不斷發(fā)展，電子商務(簡稱電商)已成為最受歡迎的消費方式之一.而在線交易量的急劇增加為物流網(wǎng)絡帶來了巨大壓力，因此高效的商品交付是電商公司提高服務水平的關鍵因素之一.電商物流中的“最后一公里”交付問題，即商品從配送站到客戶的配送過程是一個典型的有容量約束的車輛路徑規(guī)劃問題(Capacitated Vehicle Routing Problem，CVRP)，可以描述為由一個配送中心向n個客戶派遣K輛車配送貨物，每個客戶由一輛車服務且僅被服務一次.每輛車設有運送容量上限，而每個客戶有特定的貨物量需求.其目的為得到成本最低的客戶配送服務路線.

目前，利用CVRP求解算法[1-2]估算現(xiàn)實生活中各條道路的行駛成本較為困難，得到的優(yōu)化結果存在較大誤差，實際應用價值往往較低.Bertsimas等[3]認為，各條道路的行駛成本之間存在相互依賴性，因此對基于最短路徑或最短行駛時間成本矩陣的CVRP模型進行優(yōu)化難以得到符合實際需求的路徑.此外，經(jīng)驗豐富的駕駛員(專家)并不總是遵循最短路徑成本矩陣選擇配送路徑，如優(yōu)先避免易產(chǎn)生擁堵或交通信號燈多的道路.因此，在路徑規(guī)劃中還需要融入駕駛員的經(jīng)驗或偏好.Kovacs等[4]研究考慮一致性的車輛路徑規(guī)劃問題.其中，一致性包括人員一致性，即客戶希望被所熟悉的駕駛員服務[5]；區(qū)域一致性，即讓駕駛員在更為熟悉的服務區(qū)域進行配送[6-7].然而，駕駛員的偏好并不總需要完全遵循一致性原則.Liu等[8]通過分析駕駛員的連續(xù)數(shù)據(jù)軌跡來揭示駕駛員的操作模式.韋增欣等[9]建立基于多屬性決策方法的路徑優(yōu)化模型，并根據(jù)駕駛員的偏好調整模型參數(shù).Pahlavani等[10]利用局部線性神經(jīng)模糊模型學習駕駛員偏好來預測路線.He等[11]根據(jù)市區(qū)出租車司機的經(jīng)驗，利用協(xié)作偏好發(fā)現(xiàn)算法和智能司機網(wǎng)絡生成算法計算在旅行時間變化的情況下更可靠的行駛路線.黃敏等[12]根據(jù)出租車司機的經(jīng)驗, 提出約束深度強化學習算法，并在線計算不同時間段內(nèi)的最快行駛路線.以上研究多針對特定應用場景進行數(shù)據(jù)或實驗分析，因此方法的應用價值有限.

逆向優(yōu)化是指已知一個優(yōu)化問題的專家解(最優(yōu)解)，從而逆向推斷優(yōu)化模型的參數(shù)[13].Ahuja等[13-14]研究了基于最優(yōu)性條件和對偶理論的逆線性優(yōu)化問題，并將其擴展到逆向網(wǎng)絡流問題.Zhang等[15-16]利用牛頓法解決了在l∞范數(shù)下的逆向組合問題.Heuberger[17]綜述了不同逆向組合優(yōu)化問題的算法.隨著大數(shù)據(jù)應用的發(fā)展，研究者提出在逆向優(yōu)化中結合數(shù)據(jù)來設計開發(fā)算法.You 等[18]對城市卡車貨運量進行建模和預測，并提出利用觀測數(shù)據(jù)對車輛路徑規(guī)劃問題的模型參數(shù)進行校準，使得優(yōu)化結果更符合實際需求.牟健慧等[19]將逆向優(yōu)化理論和方法引入車間調度領域，建立考慮調度效率和穩(wěn)定性的數(shù)學模型，以解決多目標流水車間的逆調度問題.B?rmann等[20]結合逆向優(yōu)化和在線學習算法的思想，將模型解決方案與專家決策之間的差異設為成本矩陣的梯度，并利用迭代法令梯度逐漸降低并收斂.Aswani等[21-22]考慮了在凸優(yōu)化問題最優(yōu)解的測量結果被噪聲破壞情況下的逆向優(yōu)化問題.Saez-gallego等[23]利用逆向優(yōu)化法預測對價格敏感的消費者的電力總需求，建立電力總需求的價格響應模型，并利用廣義逆向優(yōu)化法估算模型參數(shù).

綜上分析，提出一種逆向優(yōu)化方法.通過學習配送過程中專家的過往路徑?jīng)Q策，逆向優(yōu)化得到能夠表示專家對道路選擇偏好的成本矩陣，并將該矩陣應用于CVRP模型，令配送路徑能夠融入專家經(jīng)驗.

1 逆向優(yōu)化模型

對“最后一公里”的配送問題進行建模，并根據(jù)該模型的對偶模型建立逆向優(yōu)化數(shù)學模型.

1.1 CVRP模型

配送問題描述如下：共有K輛車為n個客戶進行配送，每個客戶僅由一輛車負責服務，每輛車的容量約束為Q，第i個客戶的需求量為qi(i=1, 2, …,n), 0≤qi≤Q.定義集合V={0, 1, …,n}，其中0表示配送中心.C={cij,i,j∈V}為成本矩陣，其中cij為從客戶i到客戶j的行駛成本.優(yōu)化目標為在滿足車輛容量約束的前提下，得到總行駛成本最小的車輛配送路徑.

模型的決策變量包括xij(i,j∈V)和車輛到達客戶i時的負載ui(i∈V/{0}).如果車輛直接從客戶i行駛到客戶j，則xij=1，否則xij=0.

CVRP的數(shù)學模型(模型L)表示如下：

min ∑i∈V∑j∈Vcijxij

(1)

s.t. ∑j∈Vxij=1， ?i∈V/{0}

(2)

∑i∈Vxij=1， ?j∈V/{0}

(3)

∑j∈Vx0j=K

(4)

ui-uj+Q(1-xij)≥qi

(5)

?i,j∈V/{0}

ui≤Q， ?i∈V/{0}

(6)

ui≥qi， ?i∈V/{0}

(7)

xij∈{0, 1}， ?i,j∈V

(8)

目標函數(shù)式(1)使得總行駛成本最小化；式(2)和(3)確保每一位客戶僅被服務一次并控制網(wǎng)絡流守恒；式(4)確保配送車輛數(shù)目為K；式(5)為子回路消除約束和車輛容量約束；式(6)和(7)用于定義ui；式(8)為決策變量的屬性定義.

1.2 逆向優(yōu)化模型

min ∑i∈V∑j∈V(ηij+σij)

(9)

∑i∈V∑j∈Vφij-∑j∈V/{0}∑j∈V/{0}

[αi+βj+(Q-qi)γij+Qλi+qiμi]=0

(10)

δ≤c00+η00

(11)

(12)

βj+δ+φ0j=c0j+η0j-σ0j

(13)

?(0,j)∈J

(14)

αi+φi0=ci0+ηi0-σi0， ?(i, 0)∈J

(15)

αi+βj-Qγij≤cij+ηij

(16)

αi+βj-Qγij+φij=cij+ηij-σij

(17)

?(i,j)∈J,i,j∈V/{0}

∑j∈V/{0}γij-∑j∈V/{0}γji+λi+μi≤0

(18)

?i∈V/{0}

αi,βj,δ∈R， ?i,j∈V/{0}

(19)

λi,φij≤0， ?i,j∈V/{0}

(20)

γij,μi,ηij,σij≥0， ?i,j∈V/{0}

(21)

式中：αi,βj,δ,γij,λi,μi和φij分別為模型L中式(2)～(8)的對偶變量.目標函數(shù)式(9)使C的變動量最小化；式(10)確保x*成為模型L的最優(yōu)解；式(11)～(18)為對偶約束；式(19)～(21)為對偶變量的屬性定義.

利用模型Inv-L的迭代算法(Inv-L Algorithm，ILA)，輸入過往專家解，迭代多次求解Inv-L模型，每次求解時對C進行變動，最終得到一個能夠代表專家經(jīng)驗的成本矩陣.然而，該算法在每次迭代時均需要對模型Inv-L求最優(yōu)解，因此求解時間較長.對此，開發(fā)一種在線學習算法以實現(xiàn)對過往專家解的學習.

2 乘性權重更新算法

乘性權重更新(Multiplicative Weights Updates，MWU)算法[24]采用迭代更新的方式，在每一次迭代中，通過定義每輪的損失成本，并利用乘法更新規(guī)則對決策進行更新，能夠使得總損失成本的期望最小化.采用改進的MWU算法，將原MWU算法中基于列向量的運算擴展至基于二維矩陣的運算，并針對問題屬性確定算法學習率，使得改進后的算法具備收斂性.

2.1 算法描述

算法中的相關參數(shù)如表1所示.

表1 參數(shù)定義Tab.1 Definition of parameters

wt+1=wt-μ(wt⊙mt)

(22)

式中：a⊙b表示維數(shù)相同的兩個矩陣對應分量相乘.

MWU算法偽代碼

fort=1, 2, …,Tdo

xt=arg min (Ct·xt|xt∈X(pt))

mt=0

else

wt+1=wt-μ(wt⊙mt)

end if

end for

returnC1,C2, …,CT

對于第t輪，算法首先通過標準化wt得到Ct，并基于Ct和X(pt)正向求得xt.然后計算mt.最后基于mt和式(22)將wt更新為wt+1.經(jīng)過反復迭代，得到一個成本矩陣序列C1,C2, …,CT，其中CT為最終學習得到的代表專家經(jīng)驗的成本矩陣.由于MWU算法具有收斂性，所以可以保證經(jīng)過足夠多次學習迭代后得到的CT無限接近于Ctrue.

在 MWU算法中, 僅存在一個循環(huán)體，在每個循環(huán)中均需要調用數(shù)學優(yōu)化技術CPLEX執(zhí)行一次CVRP的求解.由于CPLEX求解線性規(guī)劃的內(nèi)置算法為分支定界算法，所以MWU算法的時間復雜度為O(n′2n′)，其中n′為輸入?yún)?shù).

2.2 算法收斂性

(23)

式中：abs(·)為求數(shù)據(jù)絕對值的函數(shù).

(24)

即MWU算法具有收斂性.

3 算例分析

利用不同維度的算法指標評價MWU算法的有效性，然后分別利用VRPLIB算例庫中的隨機算例和電商真實算例進行實驗分析，所需計算機配置為Intel?CoreTMi7-8500@3.60 GHz，內(nèi)存為8 GB.程序由MATLAB R2019b進行編譯，模型求解調用CPLEX IBM 12.9.0.

VRPLIB算例庫中的隨機算例包括成本矩陣、車輛容量和客戶需求等信息.基于給定的成本矩陣，通過更改客戶需求和車輛容量能夠生成多個不同算例，求得各算例的最優(yōu)解并作為專家解.將同一規(guī)模下的專家解分成兩部分，其中90%的專家解為訓練集，10%的專家解為驗證集.

3.1 算法指標

從算法的有效性、可實現(xiàn)性和收斂性等方面提出以下評價指標.

平均目標函數(shù)值偏差：

其中：M為驗證集的數(shù)量.該指標用于衡量由逆向優(yōu)化學習算法得到的最優(yōu)解與真實專家解目標函數(shù)值之間的平均偏差.

平均收斂誤差：

該指標用于衡量算法的收斂誤差.

平均編輯距離：

其中：Edit為求解兩個變量編輯距離的函數(shù).編輯距離為將一個字符串轉換為另一個字符串所需要的最小操作數(shù).該指標用于衡量驗證集中的專家解和由逆向優(yōu)化學習法得到的最優(yōu)解之間的平均編輯距離.

平均相似度：

其中：|H|為算例規(guī)模.判斷兩個解的相似度需要同時考慮E和|H|.此外，評價指標還包括算法的學習時間(CPU).

3.2 基于隨機算例的實驗分析

3.2.1學習效果分析 分別利用ILA和MWU算法求解具有10個客戶點的算例.當訓練集數(shù)(N1)和驗證集數(shù)(N2)逐漸增大時，兩種算法的性能分析如圖1所示.可知，在指標ε1，ε2，S和CPU方面，MWU算法均優(yōu)于ILA.隨著訓練集數(shù)量的增多，兩種算法的學習效果均有所提升，提升趨勢均為先激增后趨于平緩，且MWU算法的提升幅度更顯著.因此，在求解大規(guī)模算例時，可以考慮權衡訓練時間和訓練效果.雖然MWU算法需要花費較長時間學習CT，但在后續(xù)過程中，可以直接利用CT求解新問題.而其在求解新問題時的效率等同于求解一個CVRP問題時的效率.因此，前期訓練完成后，MWU算法求解新問題的效率較高.

圖1 兩種算法的效果對比Fig.1 Comparison of two learning algorithms

根據(jù)上述結論，選擇MWU算法求解后續(xù)算例，結果如表2所示.可知，與驗證集中已知專家解相比，MWU算法可以提供S=88.05%和ε1=0.16%的配送路徑，表明了算法在學習專家經(jīng)驗方面的有效性.此外，MWU算法的ε2=1.40，表明算法的收斂性較好.

表2 MWU算法求解不同規(guī)模算例的結果

3.2.2網(wǎng)絡圖結構的敏感性分析 上述隨機算例均為非對稱網(wǎng)絡，而在較多配送場景中，網(wǎng)絡結構呈對稱性.因此，需要分析MWU算法應用于兩種不同網(wǎng)絡結構時的性能，結果如圖2所示.保留非對稱網(wǎng)絡成本矩陣的上半部分不變，可以將其轉化為對稱成本矩陣.

圖2中，與求解非對稱網(wǎng)絡的算例相比，MWU算法在求解具有對稱網(wǎng)絡結構算例時的ε1值略低，而ε2值相對一致.表明算法的收斂性不受網(wǎng)絡對稱性影響.由于對稱網(wǎng)絡中存在多個最優(yōu)解，會對算法學習造成干擾，所以對其進行學習后所得求解結果的相似度較差.此外，隨著算例規(guī)模增大，對稱網(wǎng)絡的訓練時間呈逐漸增長趨勢.

圖2 網(wǎng)絡結構對稱性敏感性分析Fig.2 Sensitivity analysis of two networks

圖3 研究區(qū)域的客戶聚類分布Fig.3 Clustered customer network of studied district

3.3 基于電商真實數(shù)據(jù)的算例分析

采用某電商在2018年9月至11月的日常家用電器配送路線進行算例分析.服務區(qū)域為上海市徐匯區(qū)，收集到包括422輛廂型貨車在內(nèi)的20萬條訂單的配送路線.數(shù)據(jù)集包括配送的車次、時間和地址等.根據(jù)企業(yè)反饋，選擇工作效率較高的駕駛員，收集其配送路徑，并剔除異常解和存在客戶時間窗約束等不屬于研究范圍的解，但此時所得專家解并不等同于基于最短路徑矩陣的CVRP模型最優(yōu)解.將駕駛員的配送路徑作為專家解進行學習，以得到能夠代表駕駛員經(jīng)驗的成本矩陣，并用于求解新問題.

3.3.1數(shù)據(jù)預處理 數(shù)據(jù)預處理包括缺省值、異常值和地址的經(jīng)緯度轉換處理.為了增加客戶點的重復度，減少模型訓練難度，根據(jù)客戶的地理位置采用K-means算法進行聚類.從算法的可實現(xiàn)性、聚類合理性及其聚類指標等方面綜合分析，選擇聚類數(shù)為55，客戶聚類情況如圖3所示.數(shù)據(jù)預處理后共得到185個真實算例.

3.3.2結果分析 各聚類質心的歐氏距離成本矩陣為C1，算法學習過程為從185個算例中隨機選擇184個并利用MWU算法進行訓練，將剩余的一個算例作為驗證算例.將學習得到的CT代入路徑規(guī)劃模型L，求解驗證算例，并對比所得解與已知專家解.按上述過程訓練185次，求每一次訓練效果的平均值.由于真實場景下無法得到真實的成本矩陣，所以本節(jié)不考慮指標ε1和ε2.此外，利用C1求解算例，并對比采用算法前后的優(yōu)化效果，如表3所示.

表3 真實算例結果對比Tab.3 Comparison of real instances

可知，當基于CT求解算例時，求解結果明顯優(yōu)于基于C1的求解結果，在指標E和S上分別優(yōu)化了37.38%和24.67%.同時，所得配送路徑與專家解之間的平均相似度為71.84%，表明當真實算例的數(shù)量增加時，算法的學習效應也將隨之加強.

4 結語與展望

對電商物流中的“最后一公里”交付問題進行研究.在現(xiàn)實生活中，經(jīng)驗豐富的駕駛員(專家)并不總是遵循最短路徑成本矩陣選擇配送路徑.因此，采用逆向優(yōu)化法學習專家的過往路徑?jīng)Q策方案，訓練得到代表專家對路徑選擇偏好的成本矩陣，將其應用到路徑規(guī)劃模型中求解.其中，利用MWU算法實現(xiàn)對專家經(jīng)驗的學習，并設計隨機算例和電商真實物流算例進行驗證.結果表明：① MWU算法在收斂性、有效性和可實現(xiàn)性方面表現(xiàn)優(yōu)異，利用學習所得成本矩陣可以得到近似專家的配送路徑；② 算法收斂性不受網(wǎng)絡結構非對稱性和對稱性的影響；③ 在求解真實算例時，MWU算法將所得解與專家解的平均相似度提升至71.84%，表明該算法能夠有效利用過往經(jīng)驗，提供與專家解相似度極高的配送路徑.

在后續(xù)研究中，可以將逆向優(yōu)化法應用于具有不確定性車輛路徑規(guī)劃等更復雜的優(yōu)化問題.同時，現(xiàn)實生活中的專家解集往往來源于多個專家，如何對專家解進行分類并提供全局逆向優(yōu)化方案，同樣值得探討.