基于時空狀態(tài)網(wǎng)絡的電動物流車輛路徑優(yōu)化方法

楊森炎，寧連舉，商攀

(1.北京郵電大學，a.現(xiàn)代郵政學院(自動化學院)，b.經(jīng)濟管理學院，北京100876；2.北京交通大學，交通運輸學院，北京100044)

0 引言

電動車輛具備低能耗、低排放、低噪聲等優(yōu)勢，已廣泛用于城市“最后一公里”物流配送過程。由于電池續(xù)航里程有限，電動車輛需要在途中補充電量。充電設施不足導致車輛繞路或排隊充電現(xiàn)象[1]。為提高電動物流車輛服務效率，學者們基于經(jīng)典的車輛路徑問題，考慮充耗電過程、電池容量等限制條件，提出電動車輛路徑優(yōu)化問題(EVRP)。Schneider 等[2]假設充電量與充電時間呈線性關系，研究帶時間窗的EVRP問題(EVRPTW)，設計了可變鄰域搜索算法與禁忌搜索結(jié)合的混合啟發(fā)式求解算法。Montoya等[3]提出分段線性逼近方法，研究非線性充電過程的EVRPTW 問題。Goeke 等[4]研究電動車輛和內(nèi)燃機車混合下的EVRPTW問題。Felipe等[5]研究部分充電策略和不同充電技術(shù)下的電動車輛路徑規(guī)劃和充電決策優(yōu)化方法。Keskin等[6]對比了部分充電和完全充電策略下的路徑優(yōu)化方案。Hiermann 等[7]研究考慮不同車輛配置的EVRPTW 問題，提出基于自適應大鄰域搜索的混合啟發(fā)式算法，通過智能局部搜索改進路線。

既有關于EVRPTW 問題的研究在車輛行駛、充耗電過程、充電站能力等方面考慮的約束條件較為復雜。傳統(tǒng)的基于物理運輸網(wǎng)絡的建模方法難以充分考慮客戶需求和電動車輛運行狀態(tài)的時空特性，需要設置時間、空間、載重狀態(tài)、電量狀態(tài)等多種類變量，建立大量的復雜約束條件以描述變量之間的相互關系，因此模型結(jié)構(gòu)較為復雜。常采用的智能算法雖然通用性強，時效性較高，但難以從理論上保證解的質(zhì)量。

近年來學者們提出采用時空網(wǎng)絡的建模方法求解車輛路徑優(yōu)化問題，通過系統(tǒng)離散化方法處理復雜約束，可以有效簡化模型，降低問題的求解難度。Mahmoudi等[8]提出時空網(wǎng)絡的建?？蚣?，求解城市交通場景下接送乘客的車輛路徑問題。Yang等[9]基于離散的時空狀態(tài)網(wǎng)絡表示方法，對考慮混合取送的物流車輛路徑優(yōu)化問題進行建模和求解。

本文運用離散的時空狀態(tài)網(wǎng)絡方法，研究考慮充電策略的電動物流車輛路徑優(yōu)化問題，對電動物流車輛路徑規(guī)劃和充耗電過程進行離散化表示，建立多商品網(wǎng)絡流優(yōu)化模型，在時間、空間和狀態(tài)維度上對電動物流車輛路徑以及充電策略實現(xiàn)同步優(yōu)化，并設計基于增廣拉格朗日松弛技術(shù)的求解算法，將原問題分解為一系列規(guī)模較小的最短路徑子問題，可以快速找到可行解。

1 模型構(gòu)建

1.1 問題描述

本文采用部分充電策略，即車輛中途只需補充支撐其完成剩余配送任務的部分電量，可以減少不必要的充電等待時間。考慮電池容量、車輛承載能力、充電站能力、客戶服務時間窗、路網(wǎng)空間結(jié)構(gòu)等約束，優(yōu)化客戶的服務順序、充電站點選擇及單次充電時間。模型假設如下：

(1)配送中心車輛充足，車輛出發(fā)時處于滿電狀態(tài)，完成配送后返回到原配送中心；

(2)車輛的耗電量與行駛時間成正比，當剩余電池容量低于最小閾值時需補充電量；

(3)充電速率固定，充電量與充電時間滿足線性關系；

(4)不考慮路網(wǎng)交通狀態(tài)的影響，車輛保持勻速行駛。

1.2 符號定義

(1)集合和元素

N——物理網(wǎng)絡節(jié)點集合；

T——時間點集合；

V——電動車輛集合；

P——客戶集合；

S——充電站點集合；

W——車輛剩余載貨量狀態(tài)值集合；

L——空間路段集合，L={(i,j)|i∈N,j∈N}

Av——車輛v對應的時空狀態(tài)弧集合，Av={(i,j,t,t′,w,w′,e,e′)|i∈N,j∈N}；

ψp,v——車輛v服務客戶配送需求的弧集合，ψp,v?Av；

ψs,v——車輛v在充電站充電的弧集合，ψs,v?Av；

i,j,j′,j″——節(jié)點編號，i,j,j′,j″∈N；

t,t′,t″——時間編號，t,t′,t″∈T；

v,v′——車輛編號，v,v′∈V；

p——客戶編號，p∈P；

s——充電站編號，s∈S；

w,w′,w″——剩余載貨量狀態(tài)值，w,w′,w″∈W；

e,e′,e″——剩余電量狀態(tài)；

(i,t,w,e)——時空狀態(tài)點；

(i,j,t,t′,w,w′,e,e′)——時空狀態(tài)??；

o——配送中心節(jié)點；

to,v——車輛v離開配送中心o的時間；

t′o,v——車輛v返回配送中心o的時間；

wo,v——車輛v離開配送中心o的剩余載貨量；

w′o,v——車輛v返回配送中心o的剩余載貨量；

evo——車輛v離開配送中心o的剩余電量；

e′o,v——車輛v返回配送中心o的剩余電量。

(2)參數(shù)

Emin——最小剩余電量閾值；

Emax——電池容量；

ξs——充電站s的充電能力，即最多可允許同時充電的車輛數(shù)量；

γs——充電站s的充電速率；

κv——車輛v的耗電速率；

χv——車輛v的最大承載能力；

[te,p,tl,p]——客戶p的服務時間窗，其中，te,p為可服務p的最早時間，tl,p為可服務p的最晚時刻；

c(i,j,t,t′,w,w′,e,e′)——時空狀態(tài)弧(i,j,t,t′,w,w′,e,e′)的配送成本。

(3)決策變量

y(i,j,t,t′,w,w′,e,e′),v——若車輛v經(jīng)過弧(i,j,t,t′,w,w′,e,e′)則為1，否則為0。

1.3 構(gòu)建離散時空狀態(tài)網(wǎng)絡

時空狀態(tài)網(wǎng)絡具有時間、空間和狀態(tài)3 個維度?？紤]到配送過程同時受車輛承載能力和電池續(xù)航能力限制，對狀態(tài)維度進行擴展，用狀態(tài)向量[w,e] 表示車輛的載貨量狀態(tài)和剩余電量狀態(tài)。車輛經(jīng)過節(jié)點i的時空狀態(tài)可以用[i,t,w,e] 表示，即車輛在t時刻行駛到節(jié)點i，載貨量w，剩余可用電量e。時空狀態(tài)網(wǎng)絡中任意兩個相鄰節(jié)點[i,t,w,e]和 [j,t′,w′,e′]之間存在一條時空狀態(tài)弧(i,j,t,t′,w,w′,e,e′)，共有4種不同類型的時空狀態(tài)弧，分別是運輸弧、充電弧、配送弧及等待弧。需要根據(jù)以下規(guī)則構(gòu)建時空狀態(tài)網(wǎng)絡。

(1)針對EVRPTW 問題，對時間和車輛載貨狀態(tài)維度進行離散化，構(gòu)建時空狀態(tài)點集合。

(2)基于物理空間網(wǎng)絡節(jié)點之間的鄰接關系，增加離散的時間和狀態(tài)維度，建立時空狀態(tài)弧集合A，包括運輸弧集合ψT、充電弧集合ψS、配送弧集合ψP和等待弧集合ψW，即A=ψT?ψS?ψP?ψW，具體如下：

添加運輸弧，對于任意i,j∈N，t∈T，Ti,j,t為t時刻路段(i,j)的運輸時間，車輛耗電速率是κv，添加 (i,j,t,t′,w,w′,e,e′)到ψT，且滿足關系：t′=t+Ti,j,t,w′=w,e′=e-κv·Ti,j,t。

添加充電弧，對于任意s∈S,t∈T，充電站s的充電速率為γs，添加(s,s,t,t+1,w,w′,e,e′)到ψS，且滿足關系：w′=w,e′=e+γs×1。

添加配送弧，對于任意t∈T，p∈P，客戶p位于弧(ip,jp)上，其配送需求為np，Tp,t為t時刻服務客戶p所需的時間，車輛耗電速率為κ′v，添加(ip,jp,t,t′,w,w′,e,e′)到ψP，且滿足關系 ：t′=t+Tp,t,w′=w-np,e′=e-κ′v·Tp,t。

添加等待弧，對于任意i∈N，t∈T，添加(i,i,t,t+1,w,w′,e,e′)到ψW，且滿足關系：w′=w,e′=e。

(3)對于任意一條時空狀態(tài)弧a∈A，設置對應的成本ca。

以8點配送網(wǎng)絡為例，客戶位于空間網(wǎng)絡上相鄰兩個節(jié)點之間，如圖1所示，假設電池最大容量Emax=10，最小電量閾值Emin=2。電動車輛從配送中心出發(fā)時處于滿載和滿電狀態(tài)，為5位客戶配送貨物，當e低于Emin時，需要選擇充電站進行充電。最優(yōu)路徑節(jié)點順序是“1→2→4→5→7→3→1”，剩余載貨量和電量的時空變化過程如圖2所示。車輛在節(jié)點7 完成充電的過程可以用兩個充電弧(7,7,4,5,3,3,2,6)和(7,7,5,6,3,3,6,10)表示。車輛完成對客戶1 的配送服務過程可以用配送弧(2,4,1,2,8,7,8,6)表示。車輛在完成客戶的配送服務之后載貨量和剩余電量均降低。

圖1 配送網(wǎng)絡示意圖Fig.1 Diagram of distribution network

圖2 電動物流車行駛過程中的時空狀態(tài)變化Fig.2 Spatial-temporal state changes during electric logistics vehicle driving

客戶服務時間窗口、車輛承載能力和電池容量等實際約束條件以及時間、載貨量和電量狀態(tài)的更新直接嵌入時空狀態(tài)網(wǎng)絡的構(gòu)建過程，縮減算法的搜索空間，提高計算效率。具體如下：

①客戶配送服務時間窗約束，車輛需要在客戶p的服務時間窗[te,p,tl,p] 內(nèi)配送貨物，即配送時間范圍是te,p≤t≤tl,p。

②車輛承載能力約束，任意車輛v在任意時刻t的剩余載貨量wv,t不能超過其承載能力χv，即0 ≤wv,t≤χv。

③電池容量約束，任意車輛v在任意時刻t的剩余可用電量狀態(tài)ev,t不小于閾值Emin，且不超過電池容量Emax，即Emin≤ev,t≤Emax。

④時間更新規(guī)則，對于任意時空狀態(tài)弧(i,j,t,t′,w,w′,e,e′)，車輛經(jīng)過相鄰節(jié)點的時間關系滿足t′=t+Ti,j,t。

⑤電量消耗過程，對于任意運輸弧(i,j,t,t′,w,w′,e,e′)，車輛經(jīng)過相鄰節(jié)點的能量更新規(guī)則為e′=e-κv·Ti,j,t。

⑥充電過程，對于任意充電弧(i,i,t,t+1w,w′,e,e′)，車輛的能量更新規(guī)則為e′=e+γs×1。

1.4 多商品網(wǎng)絡流優(yōu)化模型

以最小化總配送成本為目標，針對EVRPTW問題，構(gòu)建多商品網(wǎng)絡流優(yōu)化模型為

式(1)為目標函數(shù)，最小化總配送成本；式(2)～式(4)為流量平衡約束，式(2)保證車輛從配送中心出發(fā)，式(3)保證車輛最后返回配送中心，式(4)保證每條路徑任意中間節(jié)點的輸入流和輸出流平衡；式(5)確保每位客戶的配送需求被服務；式(6)保證每個充電站同時充電的車輛數(shù)不超過其充電站能力；式(7)為決策變量。

為保證使用盡可能少的電動車輛，設置從配送中心直接返回到配送中心的虛擬弧，虛擬弧的成本為0。對于無需完成配送任務的車輛，直接通過虛擬弧返回配送中心，不會進入正常的服務網(wǎng)絡。

時空狀態(tài)網(wǎng)絡建模方法是在傳統(tǒng)物理網(wǎng)絡上增加了時間和狀態(tài)維度，可以清晰地刻畫車輛運行過程中時間變化、空間變化及狀態(tài)轉(zhuǎn)移之間的耦合關系，構(gòu)建復雜場景下的網(wǎng)絡流優(yōu)化模型，適用于解決帶有時間窗約束或者時變運輸網(wǎng)絡優(yōu)化問題。雖然該方法可以顯著減少EVRPTW 問題的變量種類和約束數(shù)量，但是求解上述0-1 整數(shù)規(guī)劃模型仍然屬于NP-hard，時間和狀態(tài)維度的離散化導致搜索空間增大，在計算復雜度上存在較大的挑戰(zhàn)。

2 算法設計

2.1 模型分解

為使模型簡潔，用a表示時空狀態(tài)弧(i,j,t,t′,w,w′,e,e′)。采用拉格朗日松弛方法，把上述整數(shù)規(guī)劃模型的式(5)和式(6)松弛到目標函數(shù)式(1)中，得到拉格朗日松弛模型(LR)為

保留約束條件式(2)～式(4)和式(7)。

分解可得單輛車的LR子模型為

化簡為

成本為

該問題是標準最短路徑子問題，具有相同的拉格朗日乘子λp和λs,t。每次迭代，相同車輛會傾向于選擇服務同一客戶，導致解的對稱性問題，影響求解效率。為克服對稱性問題，在LR 模型的基礎上增加一個增廣二次項，構(gòu)建增廣拉格朗日松弛模型(ALR)為

式中：θp和θs為懲罰參數(shù)。二次懲罰項的引入導致ALR模型難以分解。為把ALR模型分解為易求解的子模型，對模型進行線性化處理。定義αp,v和βs,t,v為

式中：αp,v為客戶p除了車輛v之外被其他車輛服務的總次數(shù)；βs,t,v為除了車輛v之外在充電站s充電的總車輛數(shù)。車輛v對應的子問題目標函數(shù)為

由于ya,v是0-1 變量，二次項等價于，進一步化簡為

根據(jù)式(16)和式(19)，式(15)可進一步化簡為

式中：?為常數(shù)項。

2.2 求解算法

采用塊坐標下降法(BCD)的分解優(yōu)化框架，嵌入前向動態(tài)規(guī)劃(DP)算法，對ALR 模型分解得到的最短路徑子模型進行循環(huán)迭代求解。在BCD框架中，變量和約束條件按車輛被分為若干個模塊，針對每個模塊依次最小化目標函數(shù)[10]。假設配送中心有M輛電動車，編號m=1,2,…,M，任意車輛vm對應的時空狀態(tài)弧變量集合設為Ym，M輛電動車的時空狀態(tài)弧變量總集合設為Y。ALR 模型的目標函數(shù)為

式中：y∈Y=Y1×Y2×…×Ym?Rn，Ym?Rnm，變量y1,y2,…,yM依次更新公式為

在第k+1 次迭代，為了優(yōu)化第m輛電動車的路徑，完成前i-1 輛電動車的時空狀態(tài)路徑優(yōu)化后，更新車輛路徑方案，后M-m輛車的路徑仍然保持第k次迭代的優(yōu)化方案懲罰參數(shù)θp和θs分別是拉格朗日乘子和的更新步長，即

基于ALR 的方法可以通過增設二次懲罰項，克服LR 方法存在的對稱性問題，提高算法的收斂速率；另一方面，每次迭代計算最優(yōu)上界和下界的間隙GGAP，可以從理論上評估可行解的質(zhì)量，具體步驟如下。

Step 1 初始化

Step 2 計算ALR模型的下界解

Step 2.1 求解子模型

按式(21)更新成本c^a,v，不考慮充電站的充電能力約束，使用DP 算法，按式(23)的迭代方式，依次優(yōu)化單輛車的配送路徑。

Step 2.2 計算下界解，更新最優(yōu)下界值

根據(jù)Step 2.1 得到的車貨分配方案，計算下界解，規(guī)則為：若某些客戶被多輛車服務，則僅指派一輛車配送即可；若存在部分服務需求沒有被滿足，則重新指派一輛車完成服務。

計算下界值，并按更新最優(yōu)下界B*LB。

Step 3 計算ALR模型的上界解

Step 3.1 求解子模型

Step 3.2 計算上界解，更新最優(yōu)上界值

根據(jù)Step 3.1 得到的車貨分配方案，計算上界解，規(guī)則同Step 2.2。

計算上界值，并按更新最優(yōu)上界B*UB。

Step 4 評估解的質(zhì)量

計算B*UB和B*LB之間的最優(yōu)間隙GGAP為

GGAP值越小，解的質(zhì)量越好。

Step 5 更新拉格朗日乘子

采用次梯度方法，按式(24)和式(25)更新拉格朗日乘子λp,λs,t。

Step 6 結(jié)束終止條件

若k到達最大迭代次數(shù)K，終止算法，并輸出，B*UB，GGAP和最優(yōu)上界解；否則，返回Step 2，k=k+1，繼續(xù)迭代。

Step 2.1和Step 3.1的DP算法采用了隱枚舉的思想，盡早去除不滿足約束條件的變量取值組合，直至找到一個滿足條件的可行解，并記為當前最優(yōu)的可行解，計算其目標函數(shù)值。后續(xù)迭代得到的可行解的目標函數(shù)值與當前最優(yōu)值進行比較，對可行解的質(zhì)量進行改進。

求解最短路徑子問題的過程是上述ALR 算法最耗時的環(huán)節(jié)。 整個算法的復雜度是O(K×nV×nL×nT×nW)，其中，K為迭代次數(shù)，nV為電動車輛數(shù)，nL為空間路段數(shù)，nT為離散的時間點數(shù)，nW為車輛承載能力。采用時空狀態(tài)網(wǎng)絡的建模方法需要確定合適的離散精度，若離散的時間、空間及狀態(tài)精度過高，時空狀態(tài)網(wǎng)絡規(guī)模過大，會影響算法的計算效率；若離散精度設置過小，會影響模型的精確性。

3 算例分析

算法執(zhí)行的硬件環(huán)境為Intel(R)CPU(TM)i7-7700 CPU @ 3.60 GHz，內(nèi)存為16 GB。軟件環(huán)境為Windows 10系統(tǒng)，使用Python編寫求解算法。

3.1 算法分析

基于Sioux Falls 網(wǎng)絡構(gòu)建測試算例，如圖3所示，共有24個節(jié)點和76條路段，其中，節(jié)點11設為配送中心，節(jié)點1、2、8、13、15、18、23分別設為充電站s1,s2,s8,s13,s15,s18,s23，每條路段上標記的是行程時間。電動車輛從配送中心出發(fā)時均處于滿載(w0=8)和滿電狀態(tài)(e0=30)。充耗電過程的相關參數(shù)設置為：最小剩余電量閾值Emin=3，最大電池容量Emax=30，充電速率γ=10，耗電速率κ=1，充電站能力ξ=2。

圖3 Sioux Falls網(wǎng)絡Fig.3 Sioux Falls network

為測試算法求解不同規(guī)模算例的性能，設置兩組實驗，分別為24位客戶和50位客戶，得到電動車輛和充電站使用方案如表1所示。24 個客戶的算例需使用6 輛車和4 個充電站(s8,s13,s15,s23)，其余充電站均未有車輛充電。50 個客戶的算例需使用10 輛車和6 個充電站(s2,s8,s13,s15,s18,s23)?？蛻魯?shù)量增大需要調(diào)用更多的電動車輛完成配送任務，也會提高充電站的利用率。

表1 優(yōu)化方案Table 1 Optimization scheme

兩組不同規(guī)模算例的實驗結(jié)果如表2所示，可以看出，算法能夠在有限的迭代次數(shù)內(nèi)收斂，計算效率較高。最優(yōu)上界值與最優(yōu)下界值之間的GGAP較小，表明計算得到的可行解已接近最優(yōu)解，質(zhì)量較高?？蛻魯?shù)量越多，服務客戶需要的電動車輛越多，算法收斂需要迭代的次數(shù)越多，計算效率有所降低，但仍然可以較快地找到可行解。在24 個客戶算例中，算法迭代次數(shù)K設為100，ALR 算法每次迭代得到的最優(yōu)上界和下界變化如圖4所示，迭代43次后，最優(yōu)上界值、最優(yōu)下界值以及兩者之間的GGAP均保持不變，即算法收斂到較優(yōu)的解。

表2 不同規(guī)模算例的實驗結(jié)果Table 2 Experimental results of different-scale instances

圖4 ALR最優(yōu)上界和下界變化Fig.4 Variations of ALR best upper and lower bounds

3.2 參數(shù)靈敏度分析

(1)充電速率

把24 個客戶算例設為基準算例，充電速率γ調(diào)為5，其余參數(shù)保持不變。算法迭代63 次后收斂，CPU計算時間為79.51 s。充電速率降低使電動車輛在充電站的充電時間變長，由于客戶服務時間窗是固定的，導致車輛路徑和充電優(yōu)化方案整體上發(fā)生較大變化。如表3所示，與基準算例(γ=10)相比，總時間成本增大7.9%，總充電量提高14.3%，總充電時長顯著增大。

表3 充電速率的影響Table 3 Influence of recharge rate

(2)充電站能力

在基準算例的基礎上，充電站能力ξ調(diào)為1。算法迭代62 次后收斂，CPU 計算時間為74.93 s。如表4所示，與基準算例(ξ=2)相比，總時間成本有所降低，總剩余電量增大，故充電站能力參數(shù)也會影響車輛整體的路徑規(guī)劃和充電決策。

表4 充電站能力的影響Table 4 Influence of recharging station capacity

4 結(jié)論

本文綜合考慮了電動物流車輛路徑規(guī)劃、客戶配送服務需求及充耗電過程，基于高維的離散時空狀態(tài)網(wǎng)絡，構(gòu)建了多商品網(wǎng)絡流優(yōu)化模型，設計了基于ALR 和BCD 技術(shù)的分布式優(yōu)化算法，基于Sioux Falls 網(wǎng)絡測試算例，驗證了算法的時效性和可靠性，并得出以下結(jié)論：

(1)時空狀態(tài)網(wǎng)絡建模框架可以有效減少變量類型和耦合約束，簡化模型結(jié)構(gòu)；

(2)ALR 模型中懲罰項的引入可以克服LR 解的對稱性問題，加快算法的收斂速率，且最優(yōu)上下界的計算可以評估可行解的質(zhì)量。

(3)考慮部分充電策略的電動車輛路徑規(guī)劃可以降低能源消耗和運營成本，在時間、空間和狀態(tài)維度上實現(xiàn)電動物流資源的優(yōu)化配置。