隨機時變車輛路徑問題的多目標魯棒優(yōu)化方法

段征宇 ，雷曾翔 ，孫 碩 ，楊東援

（1.同濟大學道路與交通工程教育部重點實驗室，上海 201804；2.上海市城市規(guī)劃設計研究院，上海 200040）

車輛路徑問題（vehicle routing problem，VRP）是物流配送的核心問題之一，研究如何安排車輛配送路線，使得配送成本最低.VRP 也是NP-hard（nondeterministic polynomial-time hard）問題，吸引眾多學者進行了大量的研究工作.以往研究大多假定路段行程時間是常數(shù)或確定型時變函數(shù).然而，隨著城市交通擁堵的日趨嚴重，交通需求的短期變化、交通事故和惡劣天氣等，導致了路網交通狀態(tài)的不確定性.因此，在物流配送中需要同時考慮路網交通狀態(tài)的時變性和隨機性，以滿足客戶的時效性和時間窗要求，提高客戶滿意度.

隨機時變車輛路徑問題 （stochastic time-dependent vehicle routing problem，STDVRP）是研究在隨機時變路網中，如何制定滿足多個客戶需求的配送車輛路徑，使得總配送成本最小.STDVRP 是時間依賴型車輛路徑問題 （time-dependent vehicle routing problem，TDVRP）和隨機車輛路徑問題 （stochastic vehicle routing problem，SVRP）發(fā)展而來，同時考慮了路段行程時間的時變性和隨機性.在實際路網中，STDVRP 需要大量計算客戶節(jié)點之間的最優(yōu)路徑，因此，隨機時變最優(yōu)路徑問題 （stochastic time-dependent optimal path problem，STDOPP）是STDVRP 的基礎.

早期STDOPP 研究以最小期望時間作為優(yōu)化標準.比如：Hall[1]提出了分支限界方法與k最短路徑相結合的求解方法，并發(fā)現(xiàn)由于不滿足Bellman 優(yōu)化原理，傳統(tǒng)最短路徑算法不能直接用于STDOPP；Miller-hooks 等[2]設計了求解最小期望時間路徑的標號擴展算法.20 世紀90年代后，效用理論被引入最優(yōu)路徑的評價.比如：Wellman 等[3]提出了隨機一致性條件，根據效用函數(shù)計算期望最短路徑；Huang等[4]采用行程時間的負效用函數(shù)評價隨機時變路徑，設計了標號修改算法.以上研究以路段行程時間的概率分布為基礎，求解算法的時間復雜度為指數(shù)形式，不適用于大規(guī)模路網.另一個思路是考察路段行程時間“最壞情況”下的“最優(yōu)路徑”.Sim[5]最早提出了魯棒優(yōu)化模型，將隨機路網的最優(yōu)路徑問題 （stochastic optimal path problem，SOPP）轉變?yōu)榇_定型路網的最短路徑問題.Sun 等[6-7]將該方法擴展到STDOPP，將其轉化為確定型時變網絡的最優(yōu)路徑問題，求解算法的時間復雜度為多項式形式.

目前STDVRP 的研究相對較少，主要根據路段行程時間的概率分布，采用機會約束模型進行建模，使得期望行程時間最小[8-9].無時間窗STDVRP 的模型相對簡單，但計算復雜度仍然較高，例如Nahum等[10]設計的遺傳算法，對于100～200 個客戶節(jié)點的STDVRP，平均計算時間為3.8 h.對于有時間窗STDVRP，需要考慮各客戶的時間窗約束，求解更為困難，相關研究還很少見.

目前的STDOPP 和STDVRP 研究，通常需要知道路段行程時間的概率分布，求解的計算復雜度較高.對于實際路網的應用，要標定每個路段的行程時間概率分布是比較困難的.魯棒優(yōu)化方法提供了一種新的思路，已被成功應用于SOPP[5]、STDOPP[6-7]和SVRP[11].此外，目前大多STDVRP 研究針對無時間窗問題，而對于有時間窗問題，需要考慮各客戶的時間窗約束，求解更為復雜.

本文建立了帶硬時間窗的STDVRP 的多目標魯棒優(yōu)化模型，設計了一種蟻群算法進行求解.與傳統(tǒng)的STDVRP 模型相比，該模型不需要知道路段行程時間的概率分布函數(shù)，只要根據歷史交通數(shù)據標定其波動范圍，在實際路網中使用更為方便.由于基于魯棒優(yōu)化模型的STDOPP，計算復雜度為多項式形式，因此，相比于基于路段行程時間概率分布的STDVRP 模型，魯棒STDVRP 模型的計算復雜度更低.此外，該模型考慮了各客戶的時間窗約束，與實際物流配送應用的情況更為一致.

1 STDVRP 建模

1.1 隨機時變路網的表示和標定

對于隨機時變路網G={V,A,T,Cab(t)}，其中：V為節(jié)點集合；A為路段集合；Cab(t)為路段lab在時段t（t∈T）的行程時間，如式（1），lab是從節(jié)點a到b的有向路段，lab∈A；T是時間分段的集合（1，2，· ··,M），M是路段行程時間函數(shù)Cab(t)的分段數(shù).

式中：Rab(t)在 時段t是一個確定值；δab(t)是一個取值范圍為[0,dab(t)]的 隨機變量，dab(t)在時段t也是一個確定值.

因此，Cab(t) 的取值范圍為 [Rab(t),Rab(t)+dab(t)].

假定隨機時變路網G滿足隨機一致性條件[3].也就是，對于任意路段lab∈A，若t≤t′（t'為晚于t的另一個時段），對于給定時間ξ，式（2）成立，也就是出發(fā)時間早的車輛早到終點的概率比出發(fā)時間晚的車輛大.

隨著浮動車技術的發(fā)展，獲取城市路網的路段行程車速也越來越容易.例如，上海2007年建成的浮動車采集系統(tǒng)，接入了2.5 萬輛出租車GPS 數(shù)據，可以計算得到各路段每5 min 時段的行程車速[12].采用歷史數(shù)據，可以標定各路段每5 min 時段的行程時間的波動范圍，從而得到路段行程時間函數(shù)Cab(t).

1.2 STDOPP 和STDVRP 模型的符號和變量

（1）常量

[0，te0]為配送中心的工作時間，開始時刻為0，結束時刻為te0；[tbi，tei]為客戶節(jié)點i要求的時間窗，開始時刻為tbi，結束時刻為tei；tsi為客戶i的服務時間；qi為客戶i的配送需求量；Qk為車輛k的載重量；B為無窮大的數(shù).

（2）集合

NS為始發(fā)節(jié)點集合；NE為終到節(jié)點集合；ND為需求節(jié)點集合；NSD為始發(fā)節(jié)點與需求節(jié)點的集合；NDE為需求節(jié)點與終到節(jié)點集合；NSDE為所有節(jié)點的集合；NSC為客戶點集合的任一子集合，|NSC|為集合NSC的客戶節(jié)點數(shù)；Λij為從節(jié)點i到j的可行路徑的集合.

（3）決策變量和參數(shù)

K為使用車輛數(shù)；xij(t)為若車輛在時段t內出發(fā)，從i駛向j，值為1，否則為0；ti為配送車輛到達節(jié)點i的時刻；τi為配送車輛從節(jié)點i的出發(fā)時刻；vik為若點i由車輛k訪問，值為1，否則為0；tij（τi）為時刻τi從節(jié)點i出發(fā)到j的最優(yōu)路徑的行程時間；λij為從節(jié)點i出發(fā)到j的一條可行路徑（λij∈Λij）；f~(λi j,τi)為時刻τi從節(jié)點i出發(fā)到j，路徑λij在最壞情況下的行程時間；lab為可行路徑λij中的路段，即lab∈λij；yab(t)為若路徑λij在時段t內經過路段lab，值為1，否則為0；t0k為車輛k從配送中心出發(fā)的時刻；tdk為車輛k返回配送中心的時刻.

1.3 STDOPP 的魯棒優(yōu)化模型

在STDVRP 中需要計算各客戶節(jié)點之間以及客戶節(jié)點與配送中心之間的最優(yōu)路徑，因此STDOPP是求解STDVRP 的前提.對于時刻τi從節(jié)點i到j的路徑λij，最壞情況下的行程時間為

所以，根據最小最大準則，時刻τi從節(jié)點i到j的最優(yōu)路徑的行程時間為

根據文獻[6]，結合式（3），可得

因此，時刻τi從節(jié)點i到j的最優(yōu)路徑的行程時間為

在式（6）中，Rab(t)+dab(t)是路段lab在時段t的行程時間上界，是一個確定值.若令Rab(t)+dab(t)為路段lab在時段t的行程時間，那么，可以得到一個確定型時變路網G′.可以證明G′是一個先入先出（first in first out，F(xiàn)IFO）網絡[6].這樣就將節(jié)點i到j的最優(yōu)路徑問題.轉換為FIFO 時變路網的行程時間最小路徑問題.傳統(tǒng)最短路徑的標號法可以用于求解該問題，并且時間復雜度為多項式形式.

1.4 STDVRP 的多目標魯棒優(yōu)化模型

下面給出帶硬時間窗的STDVRP 多目標優(yōu)化模型.硬時間窗約束要求在客戶指定的時間窗內配送，先到必須等待，且不允許遲到.VRP 常用的優(yōu)化目標有：使用車輛數(shù)最少、總行程時間最少、總等待時間最少、總行駛距離最小等.傳統(tǒng)方法通常采用加權求和方法，將多個目標加權組合為單個目標進行求解.但由于量綱不同，多個目標之間往往不能直接比較，權重也很難確定.

對于STDVRP，總的配送成本主要由車輛固定成本、路徑行程成本、客戶懲罰成本組成.車輛固定成本與使用車輛數(shù)有關；路徑行程成本與總行程時間有關；對于帶硬時間窗的STDVRP，不允許違反客戶時間窗要求，所以，不需要考慮客戶懲罰成本.因此，在本文中選取使用車輛數(shù)最少、總行程時間最少兩個優(yōu)化目標，采用最小最大準則的魯棒優(yōu)化方法對STDVRP 進行建模.

（1）目標函數(shù)

使用車輛數(shù)最少：

總行程時間最少：

（2）約束條件

式（9）、（10）表示每個客戶需要服務一次；式（11）、（12）表示同一條路線上的客戶由同一輛車服務；式（13）表示每個客戶只能由一輛車服務；式（14）表示每輛車的裝載量不超過該車的容量；式（15）表示所有的車輛必須從配送中心出發(fā)，并返回配送中心；式（16）防止在配送路徑中形成子回路；式（17）、（18）和式（19）分別是客戶節(jié)點和配送中心的時間窗約束.

在上述STDVRP 模型中，客戶點的時間窗和服務時間是固定的，客戶節(jié)點i到j的最優(yōu)行程時間tij(τi)根據1.3 節(jié)的式（6），采用路段在對應時段的行程時間上界計算，因此，該模型等價于求確定型時變路網G′的TDVRP.對于FIFO 網絡的TDVRP，可以采用傳統(tǒng)VRP 算法求解[13]，從而降低了STDVRP的求解困難.

2 STDVRP 的算法

Pareto 最優(yōu)解也稱為非支配解，是指不能再找到另外一些解，他們的所有目標都不差于Pareto 解的相應目標，且至少有一個目標優(yōu)于Pareto 解的相應目標.Pareto 最優(yōu)解通常是一個集合.設xA、xB是多目標優(yōu)化問題 minY=[f1(x)f2(x) ···fm(x)]T的可行 解，若 ?i1∈{1,2,···,m}，fi1(xA)≤fi1(xB)，且?j1∈{1,2,···,m}，fj1(xA)＜fj1(xB)，則稱xA與xB相比是Pareto 占優(yōu)的，記作xA?xB，也稱xA支配xB[14].

STDVRP 屬于NP-hard 問題，精確求解十分困難，通常采用亞啟發(fā)式算法進行求解[10].

改進型非支配排序遺傳算法（non-dominated sorting genetic algorithm II，NSGA-II）是Deb 等[15]于2002年提出的，在多目標VRP 問題求解中得到了成功的應用[16-17].由于本文已將STDVRP 轉換為TDVRP，因此，可以用NSGA-II 算法進行求解.

蟻 群 算 法（ant colony optimisation，ACO）在TDVRP 中得到了成功的應用，可用于求解客戶節(jié)點數(shù)量為100～1 000 的大規(guī)模TDVRP，具有較高的計算效率[18-20].本文將NSGA-II 算法中的基于Pareto占優(yōu)的非支配排序引入ACO 算法，設計了求解一種多目標STDVRP 的非支配排序蟻群算法（nondominated sorting ant colony optimisation，NSACO）.

2.1 基于Pareto 占優(yōu)的非支配排序

設解集P的規(guī)模為N，將P按照某種策略進行分級排序為m1個子集（P1、P2、…、Pm1），且滿足：（1）?i2,j2∈{1,2,···,m1}，且i2≠j2，Pi2∩Pj2=φ；（2）Pk1+1中 的可行解受Pk1的 可行解的支配（k1=1，2，…，m1-1）.快速非支配排序算法依據可行解之間的支配關系對P進行分級排序，將P劃分為互不相交的子集.解集P的每個可行解設定兩個參數(shù)ni2和Si2，ni2記 錄支配可行解i2的 可行解數(shù)量，Si2是記錄被可行解i2支 配的可行解的集合，令k1=1.算法步驟如下[21-22]：

步驟1計算每個可行解的ni2和Si2，找到解集P中ni2=0 的所有可行解，將其存入當前集合Pk1；

步驟2對當前集合Pk1中 每個可行解j2，考察它所支配的可行解集合S j2，將S j2中的每個可行解r2的nr2減去1（因為支配可行解r2的 可行解j2已經存入Pk1），若nr2-1=0則將可行解存入集合Q；

步驟3將Pk1作 為第k1級 非支配解集合，Pk1所有可行解的非支配序irank=k1，然后，令k1=k1+1，Pk1=Q；

步驟4若Pk1不為空，則轉入步驟2，否則算法停止.

2.2 NSACO 算法

（1）初始化

采用文獻[23]提出的初始解構造方法，生成質量較高的初始解，由此設定路段lab的信息素初始值ηab（0）.然后，采用“蟻周模型”更新規(guī)則[24]，對各路段的信息素進行更新.

（2）通過螞蟻移動，構造配送路徑，并進行局部信息素更新

每只螞蟻從配送中心出發(fā)，根據“偽隨機比例規(guī)則”[24]確定的轉移概率，選擇下一個節(jié)點：若滿足約束條件，則選擇該節(jié)點，并記錄到禁忌表中，然后，對行經路段進行信息素更新；否則，返回配送中心，重新開始一條新的路徑，直至所有客戶節(jié)點都被訪問，則該螞蟻完成一次周游.

（3）全局信息素更新

在所有螞蟻都完成周游以后，將得到的可行解與上次迭代保留的可行解合并，得到解集P；采用快速非支配排序算法，對P進行分級排序，得到第1 級非支配解集合P1；保留P1可行解，并對其行經路段進行信息素更新.路段lab的信息素ηab的更新方法如下：

式中：ηab（u）為第u次循環(huán)的當前信息素值；ηab（u+ 1）為更新后的信息素值；Δηab（u，u+ 1）為信息素增量；LP（u）為第u次循環(huán)，P1集合的可行解對應的總行程時間的平均值；ρ為信息素揮發(fā)系數(shù).

為了避免出現(xiàn)“過早收斂”，采用“最大-最小螞蟻系統(tǒng)”[25]的信息素更新策略，將路段信息素限制在[ηmin，ηmax]范圍內.

式中：n為客戶節(jié)點數(shù)；θb為選擇最優(yōu)解的概率，取0.05[25]；h（u）為第u次循環(huán)P1的可行解中最小的總行程時間.

（4）判斷是否達到預定最大迭代次數(shù)，若符合，則輸出最優(yōu)解并結束計算；否則，返回步驟2.

3 算例分析

3.1 算例構建

Solomon 基準算例[26]在傳統(tǒng)VRP 中被廣泛使用，包括6 組共56 個算例.每個算例包括100 個客戶節(jié)點，根據客戶的空間分布，可分為3 類：C 類算例（聚集分布）、R 類算例（隨機分布）和RC 類算例（混合分布）.參照文獻[19]的方法，基于Solomon基準算例構建TDVRP 算例，然后，增加路段行程車速的波動范圍，得到STDVRP 算例.具體如下：

將算例的路段分為5 類，路段lij的類型由g（i，j）表示，由式（24）確定.

表1給出了各類路段在4 個等分時段的行程車速.對時段1、3 的行程車速，增加±20%的波動范圍；對時段2、4 的行程車速，增加±10%的波動范圍，由此得到STDVRP 算例.為了方便仿真測試時計算期望行程時間和期望等待時間，假定路段行程車速的波動服從均勻分布.事實上本文的魯棒優(yōu)化模型和算法只與行程車速的波動范圍有關，而與行程車速的概率分布形式無關.

表1 各類路段的行程車速Tab.1 Travel speed of various links

3.2 Pareto 最優(yōu)解集合及邊界解

以擴展Solomon 基準算例R201 為例，采用NSACO 算法求解，得到的Pareto 最優(yōu)解集合的分布情況，如圖1所示.

Pareto 最優(yōu)解集包含多個非支配解，為便于比較NSACO 算法和NSGA-II 算法的結果，僅選取邊界解進行比較.為表述方便，將Pareto 最優(yōu)解集中“最壞行程時間”最小的邊界解稱為“邊界解A”，將“車輛數(shù)”最小的邊界解稱為“邊界解B”.

3.3 算例求解

采用試算法確定NSACO 算法和NSGA-II 算法的參數(shù)，即每次給出某參數(shù)的一組取值，固定其它參數(shù)，根據算例的計算結果確定該參數(shù)的最優(yōu)取值.得到NSACO 算法的參數(shù)：信息強度參數(shù)α= 1；能見度參數(shù)β= 2；信息素揮發(fā)系數(shù)ρ= 0.2；偽隨機比例參數(shù)ω= 0.9；螞蟻數(shù)量γ= 10.NSGA-II 算法的參數(shù)：交叉率取0.8；變異率取 0.2.兩種算法的初始解均采用文獻[23]提出的初始解構造方法生成.

圖1 R201 算例的Pareto 最優(yōu)解的分布Fig.1 Distribution of R201’s Pareto optimal solutions

選取擴展Solomon 基準算例的C101、C201、R101、R201、RC101 和RC201，分別采用NSACO 算法和NSGA-II 算法進行求解，兩種算法均迭代200次，得到各算例的Pareto 最優(yōu)解集合.由于Pareto 最優(yōu)解集合包含多個非支配解，為便于比較，這里僅給出“邊界解A”和“邊界解B”.

采用Monte Carlo 仿真的方法，根據表1的路段行程車速及波動范圍，對最優(yōu)解的配送執(zhí)行過程進行模擬仿真，計算100 次，得到最優(yōu)解的行程時間和等待時間的期望值.計算結果如表2所示.

表2 STDVRP 算例的計算結果Tab.2 Computational results of STDVRP instances

為了方便比較，將NSACO 算法最優(yōu)解的結果除以對應的NSGA-II 算法最優(yōu)解的結果，得到NSACO算法最優(yōu)解的相對值，如表3所示.

根據表2、3，NSACO 算法的優(yōu)化結果總體比NSGA-II 算法好.對于“邊界解B”，盡管NSACO算法的最壞行程時間和期望行程時間比NSGA-II算法大4%左右；但NSACO 算法的車輛數(shù)比NSGAII 算法小3.33%.對于“邊界解A”，NSACO 算法的車輛數(shù)、最壞行程時間和期望行程時間分別比NSGA-II 算法小11.98%、17.49%和23.03%.特別是對于C 類算例（C101 和C201），NSACO 算法的優(yōu)化結果明顯優(yōu)于NSGA-II 算法.這是由于C 類算例的客戶節(jié)點呈現(xiàn)聚集分布特征，能夠更好地發(fā)揮NSACO 算法的信息素正反饋作用，相比于NSGAII 算法的隨機搜索，NSACO 算法能更快搜索到最近的客戶點.

表3 NSACO 算法最優(yōu)解的相對值Tab.3 Relative values of NSACO algorithm’s optimal solutions

4 結束語

本文針對帶硬時間窗的STDVRP，提出了一種多目標魯棒優(yōu)化模型，并基于Pareto 占優(yōu)的非支配排序方法設計了一種NSACO 算法進行求解.主要結論如下：

（1）多目標STDVRP 魯棒優(yōu)化模型可以轉換為TDVRP，與基于路段行程時間概率分布的STDVRP模型相比，求解難度更低，適用于大規(guī)模實際網絡的應用.

（2）STDVRP 測試算例表明，與常用的多目標NSGA-II 算法相比，NSACO 算法的優(yōu)化結果更好.特別是對客戶節(jié)點聚集分布的情況，NSACO 算法的收斂速度更快.

本文的STDVRP 魯棒優(yōu)化模型，只需要標定路段行程時間在各時段的波動范圍，可以保證客戶的時間窗要求，適用于大規(guī)模實際網絡的物流配送應用.進一步的研究和改進包括：

（1）魯棒優(yōu)化模型主要考慮路段行程時間最壞的情況，得到的優(yōu)化解偏保守，下一步研究將適度放寬時間窗約束和路段行程時間波動的上下界，綜合考慮配送成本和準點送達率進行建模和優(yōu)化.

（2）粒子群算法、魚群算法等智能搜索算法在多目標優(yōu)化領域也具有良好的性能，下一步工作將探索更好的多目標STDVRP 優(yōu)化算法.