張明偉，李 波，屈曉龍

(1.天津仁愛學院 經濟與管理學院，天津 301636；2.天津大學 管理與經濟學部，天津 300072)

0 引 言

車輛的碳排放量與行駛速度密切相關，受早晚交通高峰的影響，城市內車輛平均行駛速度在不同時段變化明顯，形成時變網絡，從而影響車輛的碳排放量[1]。城市車輛配送路徑問題(vehicle routing problem，VRP)已經越來越多考慮了時變網絡影響[2,3]。但是還需考慮新到達訂單對調度方案造成的擾動[4]。目前已有學者對客戶訂單的調度周期等實際配送問題進行研究[5]，并對隨機到達的訂單擾動進行重調度來獲取更優(yōu)方案[6]。車輛選擇不同路徑，平均行駛速度及行駛距離會有差異，碳排放量也會受到顯著影響[7,8]。蟻群算法在路徑優(yōu)化過程中存在效率較低、易早熟等缺點，可通過改進算法流程和改進信息素更新方式，來提高算法效率，避免陷入局部最優(yōu)[9]。

本文將以碳排放量最小作為主要優(yōu)化目標，根據隨機訂單到達的時間安排調度，研究時變網絡、柔性路徑下配送貨車動態(tài)重調度的問題。即RTVVRP-PF(rescheduling of time-varying vehicle routing problem with path flexible)。采用模擬退火與蟻群混合的算法求解模型，以避免陷入局部最優(yōu)，提高算法收斂速度；引入碳排放因子，改進信息素濃度更新方式；使用自適應精英個體繁殖策略提高算法效率。

1 問題描述與假設條件

1.1 問題描述

本文以城市內貨物配送為研究對象，貨運車輛從一個配送中心出發(fā)，為多個不同時間窗口的客戶需求進行配送。車輛可以在多條柔性路徑中，選擇到達客戶節(jié)點的最佳路徑，并且在時變網絡下，車輛速度隨道路擁堵程度、交通高峰期動態(tài)變化，車輛速度的變化能夠直接影響車輛的碳排放量大小。

(1)動態(tài)重調度對解空間的影響

調度池中新訂單的增加，將減少優(yōu)化方案中的變量約束，擴大全局優(yōu)化的解空間，從而找到更優(yōu)的可行解，如圖1所示。

圖1中，原調度方案，車輛配送訂單2時，部分行駛時間處于交通擁堵期，車輛的速度降到25 km/h，由于速度較慢，車輛在交通擁堵期將產生大量碳排放。

假如在原調度方案中的車輛還沒有出發(fā)時，加入新訂單，一起優(yōu)化調度，新訂單路徑平均速度較快(35 km/h)，部分服務時間在交通擁堵期內不產生碳排放量，并且訂單2由于重調度，避開了高峰，速度變?yōu)?5 km/h，將減少大量碳排放。因此，以一定周期加入新訂單進行多次動態(tài)重調度，比一次性的靜態(tài)調度具有更大的解空間。

(2)時變網絡下柔性路徑對碳排放量的影響

在現實的城市配送道路網絡中，網絡節(jié)點之間往往有多條路徑連通，即存在柔性路徑。在早晚交通高峰期影響下，不同路徑(路段)上，車輛平均速度變化較大。車輛在時變網絡的不同時段，通過選擇碳排放量較小的路段組成的路徑，可以減少碳排放量。

RTVVRP-PF模型以動態(tài)重調度的方式為基礎，考慮了在車輛平均行駛速度受交通影響的時變網絡下，將一次性的靜態(tài)車輛調度分成多個階段，即隨著客戶訂單的不斷到達，在積累一定時間或一定數量的訂單后，及時安排重新調度，通過新舊訂單的混合來重構解空間，以使得車輛配送盡量避開交通擁堵期，提高車輛平均行駛速度來減少碳排放量。并且為了貼合實際工況，在模型中考慮柔性路徑、動態(tài)載重、時間窗口約束等因素，并在模型中這些因素整體優(yōu)化，以達到在VRP中減少碳排放的目的。

1.2 符號說明

模型中的變量定義如下：

V：節(jié)點集合，且V=1表示配送原點，V=(2, 3, …，N) 表示客戶節(jié)點集合；

E：網絡圖中的弧集，i，j為任意兩個客戶節(jié)點，E={(i,j)|(i,j∈V且i≠j)}；

r：路徑上的路段，r∈{1, 2, …，R}；

k：重調度次數，k∈{1, 2, …，K}；

m：配送車輛， ?m∈M；

g，G：車輛的實時載貨量和最大載重；

qj：節(jié)點j的貨物需求量；

tj：車輛到達客戶節(jié)點j的時刻；

tsj：節(jié)點j的服務時間；

tsm：車輛m的裝載時間；

1.3 條件假設

(1)配送原點擁有足夠多的載重量相同的車輛。

(2)每次調度中，每個客戶只接受一輛車輛配送服務。

(3)車輛到達客戶的時間具有時間窗口限制。

2 模型構建

2.1 碳排放量計算

車輛碳排放主要由車輛行駛過程中所消耗的燃油產生，車輛燃油消耗估算的方法分為實驗法(MEET)與模型法(CMEM)兩類。模型法多采用燃油消耗量為目標建立模型，來估算碳排放量。但是模型法受到燃油不能完全燃燒的影響，車輛在不同速度下，燃油消耗量與碳排放量之間關系并不是恒定的，會發(fā)生變化，因此會影響時變網絡下計算碳排放量的準確度[10]。MEET法應用范圍較廣[11]，本文采用MEET方法，即通過對貨運車輛的尾氣排放污染進行大量實驗分析，直接給出不同類型空載的貨車在不同速度v(km/h)下的碳排放量函數θu(g/km)和系數a，并可根據配送車輛類型和性能，對系數進行調整[12]，計算公式表示為

(1)

為適合在城市內范圍進行配送，本文選用載重量為4.5噸的貨車，根據文獻[12]，a0,a1,…，a6的取值依次為110,0,0,0.000 375,8702,0,0。

貨運車輛的碳排放量與載貨量成正比關系，文獻[12]中不同類型車輛在不同速度下，車輛滿載的碳排放量θl(v)的計算公式為

(2)

根據貨運車輛的載重量，b0,b1,…，b5取值不同，裝載量4.5噸車輛的取值依次為1.27,-0.002 35,0,0,-1.33。裝載量4.5噸車輛滿載與空載相比較，不同速度下，碳排放量系數的變化如圖2所示。

可以看出，4.5噸貨車滿載下，速度v高于23 km/h后，速度越快，碳排放系數越低，單位碳排放量越少。

綜上，某貨運車輛m在速度v下，實時載貨量gm，其動態(tài)實時碳排放量θ可以表示為

(3)

2.2 服務時間

服務時間包含貨運車輛裝車、卸車和與客戶進行手續(xù)交接的時間，在服務時間內車輛不產生碳排放量。車輛裝載在配送原點完成裝載，因此其裝載時間tsm不計入車輛行駛調度模型。客戶j的服務時間tsj只需要考慮卸載時間tsj1和手續(xù)交接時間tsj2，卸載時間與客戶j的需求量qj和單位時間卸載量qt(噸/小時)相關，即：tsj1=qj/qt。 如果tsj1=0， 即沒有貨物卸載，則tsj=tsj2=0， 表示在沒有到達客戶節(jié)點時，服務時間為0?？蛻鬸的服務時間可以表示為

(4)

2.3 模型建立

RTVVRP-PF重調度時變網絡下柔性路徑低碳城市配送模型中，選取3個優(yōu)化目標：所有車輛一天配送過程產生的碳排放量總和C最小，車輛行駛的總里程D最小，車輛總行駛時間T最小。C為主要優(yōu)化目標，車輛碳排放量和D密切相關，且D也為VRP主要優(yōu)化目標，另外，由于存在配送時間窗口約束，且T與車輛司機、車輛日常管理成本密切相關，因此也非常重要。本文選取C、D、T作為模型的3個優(yōu)化目標。

T、D、C優(yōu)化目標函數分別為

(5)

(6)

(7)

目標函數約束條件S.T.

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式(8)表示車輛最大負載限制；式(9)表示客戶時間窗口約束；式(10)表示車輛完成任一客戶e的配送任務后需要離開；式(11)表示車輛從配送原點出發(fā)一次，完成全部配送任務后需要返回原點；式(12)表示路徑路段的決策變量；式(13)表示選擇客戶節(jié)點的決策變量；式(14)表示重調度決策變量。

3 算法設計

VRP屬于旅行商問題，早被驗證為NP完全(non-deterministic polynomial complete)問題，難以使用數學解析的方法進行求解[13]。智能優(yōu)化算法已經公認對此類NP完全問題有較好的優(yōu)化效果。蟻群算法(ant colony optimization，ACO)和模擬退火算法(simulated annealing，SA)在VRP優(yōu)化求解中效果較好[14,15]。ACO是一種群智能優(yōu)化算法，具有記憶性，利用螞蟻的信息素濃度變化來選擇和構造可行解。ACO每次都需要對整個種群重新構造，因此效率較低，而且容易陷入停滯。SA以某種概率接受更新解，容易跳出局部最優(yōu)和停滯，且控制簡單，編碼易實現，但是計算時間較長，尋優(yōu)效果一般，優(yōu)化效率不高。

本文將上述兩種算法的優(yōu)點結合，提出混合蟻群算法SA-ACO，算法通過對SA中退火溫度的控制，實現快速收斂，并利用SA的Metropolis準則防止陷入局部最優(yōu)和過早停滯現象。通過ACO中的多參數信息素濃度計算來更新種群，增加種群的記憶性，并設計了碳排放因子，以增強種群優(yōu)化的方向性。SA-ACO的算法流程如圖3所示。

跳出內循環(huán)的判斷標準為達到設定的某一個循環(huán)次數的常數值Cr1，跳出外循環(huán)標準為設定的某一個循環(huán)次數常數值Cr2內最優(yōu)解的適應度不再發(fā)生變化。

SA-ACO算法中按照Metropolis抽樣準則來接受個體si(t)是否更新為s′i(t)， 避免早熟或停滯。接受概率函數Ω為

(15)

其中，Г為退火為溫度控制參數，其迭代計算見式(16)，L為算法的迭代次數

(16)

3.1 算法編碼

SA-ACO算法采用實數方式編碼，每個可行解為一個兩行的矩陣，矩陣的第一行為配送客戶節(jié)點順序的信息，數字“1”表示配送原點，車輛從原點出發(fā)，完成配送任務后返回原點?？尚薪饩仃嚨诙袨槿嵝月窂竭x擇信息，并且部分節(jié)點之間的路徑包含多個順序排列的路段，如圖4所示。

圖4第一行編碼中 [1,4,3,11,1] 表示一個配送車輛從原點“1”出發(fā)，順序到達客戶節(jié)點4、3、11進行配送服務，最后返回原點。第二行路徑選擇信息，即車輛在1、4節(jié)點之間選擇了第2條路徑，4、3節(jié)點之間選擇了第1條路徑，3、11之間選擇了第3條路徑，從客戶節(jié)點11返回原點選擇了第2條路徑。其它配送車輛行走信息依次類推，直到所有客戶配送服務均完成。

3.2 適應度計算

SA-ACO算法中，種群個體的適應度與優(yōu)化目標相關，模型包含碳排放量C、總行駛路徑D和總行駛時間T這3個優(yōu)化目標，存在Pareto解集。3個優(yōu)化目標單位不同，且數量級相差較大，在個體適應度計算過程中需對優(yōu)化目標的數量級進行調整，以便適合其權重系數。本文采用自適應的適應度函數f(s)來消除多個目標數量級的差別，f(s)表示為

(17)

其中，Cmax(S)、Dmax(S)、Tmax(S) 表示某一次迭代種群中，個體s的碳排放量、總行駛路程和總行駛時間的最大值，λ1、λ2和λ3分別表示三者權重系數，需要根據企業(yè)對三者的關注程度決定權重系數。

3.3 個體更新

本文提出自適應精英個體繁殖策略來更新個體，來保證能夠搜索到全部解空間，動態(tài)保留精英個體的優(yōu)良基因，以提高算法進化效率。

SA-ACO的種群個體采用兩點局部更新的方式，如圖5所示。

在客戶配送節(jié)點信息編碼中，將每臺車輛出發(fā)的原點設為編碼1，找到所有的U個車輛配送原點“1”，隨機選取第u1和u2兩個配送原點，將兩個配送原點中間的編碼按照ACO中信息素濃度進行重新構造，兩個隨機點之外的其它部分的信息編碼直接復制，構成新的種群個體。其中這兩個點之間的長度u12的計算公式為

(18)

式(18)中int[ ]向表示下取整，fmax(S) 為這一代群體中的最大適應度值，并且如果u12>U-u1， 則u2為個體的最后一個配送原點，即u1之后的部分全部更新。

自適應精英個體繁殖策略，能夠快速搜索到全部解空間，按照信息素濃度構造個體局部，能夠使其更新具有方向性的同時，保留其它部分的優(yōu)良信息。

為增加模型中的碳排放量優(yōu)化目標效果，并兼顧總行駛里程和總行駛時間的優(yōu)化目標，本文在構造ACO的信息素能見度和信息素增量更新時，設計了包含碳排放量因子的多因子計算方式。

(19)

其中，allowedj′表示還沒有配送的客戶節(jié)點集合。

(1)包含碳排放因子的多因子能見度更新

(20)

(2)包含碳排放的多因子信息素增量計算

(21)

其中，0<σ≤1，為信息素蒸發(fā)速率。

(22)

其中，ε1、ε2、ε3為權重系數。

4 仿真算例

4.1 算例數據

為驗證RTVVRP-PF模型與算法的有效性，盡量貼近現實城市配送狀況且不失一般性，本文依據我國城市的配送數據的統計分布，利用計算機生成符合統計分布的隨機算例數據。

城市內配送，客戶的需求量一般較小，許多城市限制配送貨車的載重量，一般不使用較大型貨車。由于5噸以下的異質貨車(車型不同、載重量不同的貨車)，碳排放量差異不大。所以這里使用最大載重量G=4.5噸的同質貨車(車型一致的貨車)進行配送運輸的模擬計算。

表1 隨機節(jié)點的坐標、需求量(噸)和時間窗口

由于相同時刻，城市內各條道路的車流密度不同，車輛的平均行駛速度也不盡相同，有些城市的快速路和擁堵路段與其它道路的平均速度差異較大。隨機選取10%的路徑，按照較慢速度[20～40] km/h和較快速度[55～70] km/h，產生以5為步長的整數平均分布作為車輛行駛速度。兩節(jié)點間不同路段，速度差異較大的路徑，則按照均勻分布隨機產生速度和路段長度百分比，見表2。

表2 柔性路徑的速度/(km/h)和路段長度/%

利用交通出行網以及地圖軟件，統計整理一般城市各個時段道路的平均行駛速度，見表3。表2中的路徑為柔性路徑，所有柔性路徑的平均速度為60 km/h。按照表2中柔性路徑不同時段的平均速度與柔性路徑平均速度60 km/h的比值，來確定柔性路徑在不同時段的速度系數(見表3)。則柔性路徑在不同時段的速度為表2中的速度乘以不同時段下的柔性路徑速度系數。

柔性路徑下，兩節(jié)點之間多條路徑的長度不同，這里將第二、三條路徑按照第一條路徑的長度以分布函數按照一定比例進行放大，放大的比例服從U(1，1.3)平均分布，見表4。

表3 車輛不同時間段的平均行駛速度

表4 柔性路徑的第二、三條路徑長度/km

4.2 仿真結果

本文按照4.1節(jié)中的數據進行仿真運算。使用Matlab2019a版本進行建模優(yōu)化，種群規(guī)模取500。參數設置直接影響算法求解的效率與效果，本文采用均勻設計試驗法來確定參數。算法參數中退火初始溫度Г0=1500，降溫系數0.99，內循環(huán)次數Cr1=1000，算法終止(跳出外循環(huán))條件為Cr2=1000次迭代內，算法最優(yōu)解適應度不再變化。算法的路徑選取概率函數中α和β的數值取值為1，信息素的蒸發(fā)率σ=0.58。優(yōu)化目標權重的系數λ1、λ2和λ3取值均為1，即企業(yè)對3種優(yōu)化目標的重視程度相等。

重調度過程中，需要統計訂單的到達時間，假設7點半之前，前一日未配送訂單和新增的訂單為客戶節(jié)點2～15的訂單，并且7點半之后，每1.5小時，新增5個客戶訂單，單位貨物的卸載時間qt=0.35噸/小時，每個客戶節(jié)點的手續(xù)交接時間tsj2=0.2小時，重調度后，揀貨和裝車的時間tsm=1，訂單到達時間見表5。

表5 訂單到達時間

分別采用動態(tài)重調度和靜態(tài)調度方式進行仿真對比運算。首先采用動態(tài)重調度方式，即按照表5的訂單截至時間進行4次調度，并且每次重調度，都將上一次沒有出發(fā)的貨車和訂單重新加入新的調度方案，調度完成后，經過1小時揀貨和裝車，貨車可出發(fā)。如果重調度時貨車已開始裝車，則該貨車和訂單不能加入新的調度方案。然后采用靜態(tài)調度方式，即只進行兩次調度，在7點半進行2～15客戶節(jié)點訂單調度，在12點進行16～30客戶節(jié)點訂單調度，并且第2次調度時，不加入第1次調度中沒有出發(fā)的貨車和訂單。分別采用兩種調度方式，各進行30次優(yōu)化運算，取平均值，結果見表6。

表6 動態(tài)重調度與靜態(tài)調度結果對比

由表6可以看出，與靜態(tài)調度方式相比較，采用動態(tài)重調度方式，在碳排放、總行駛時間、總行駛里程和使用車輛數方面均有所降低，尤其碳排放量C(kg)降低了15.3%，這主要是因為重調度使得優(yōu)化空間加大，更多的車輛可以避開交通擁堵期，從而使得碳排放量等指標能夠明顯減少。

以30次優(yōu)化計算中的一個重調度方案為例，重調度優(yōu)化方案見表7，優(yōu)化后的碳排放與總行駛時間和里程見表8。

表7 重調度優(yōu)化方案

由表7可以看出，在當日7點半時，對客戶節(jié)點第2～15的訂單進行調度，共安排4輛貨車，最早貨車8點半裝載貨物后出發(fā)。截至到當日上午9點時，接到客戶節(jié)點16～20的新訂單，這些訂單與第一次調度時尚未開始裝車的訂單(分別安排在9點半、10點半裝完車)，即客戶節(jié)點2、4、5、6、9、12、15的訂單一起進行第2次重調度，然后在上午10點半與上午12點，以相同方式分別進行了第3次和第4次重調度，共使用貨車8輛。

圖6為車輛配送路徑，可以看出，部分路徑有交叉，主要原因是訂單配送是經過多次重調度產生的，以及優(yōu)化目標中增加了減少碳排放目標。

多目標函數的Pareto解與各優(yōu)化目標的權重系數相關，RTVVRP-PF的3個優(yōu)化目標可以根據企業(yè)對優(yōu)化目標的需求程度來確定。本文對3個優(yōu)化目標給出參考范圍。創(chuàng)建Pareto非支配解集合，根據支配關系進行Pareto分級，將級別最高的解放入Pareto非支配集合中，若非支配解集規(guī)模超過設定值，則進行修剪。利用上述方法求解得到10個Pareto非支配解，見表8，并分別求出了3個優(yōu)化目標的極小值分別為碳排放252.45 kg，總行駛時間33.28小時和總行駛里程631.36 km。

表9為這次調度的碳排放、總行駛時間和里程。并且，為驗證柔性路徑和動態(tài)負載對配送優(yōu)化目標的影響，在這個調度方案中，分別統計了相同調度方案下，如果使用最短單一路徑和貨車空載、滿載的情況下，車輛的碳排放量、總行駛時間和里程變化。

由表9可以得出，如果不使用柔性路徑，而使用單一路徑，并且是最短路徑，優(yōu)化后，雖然總行駛里程稍有降低，但是碳排放量增加了11.2%，總行駛時間增加了9.2%，因為在柔性路徑狀況下，車輛在道路擁堵時可以選擇速度更快、碳排放量更低的較長路徑，說明柔性路徑更加符合實際配送工況。由表9還可以得出，不考慮動態(tài)負載的影響，分別直接使用固定的貨車空載和滿載質量來計算碳排放量，與考慮動態(tài)負載相比，碳排放量分別相差7.5%和6.5%，表明動態(tài)負載對貨車碳排放量的影響是比較顯著的。

表8 Pareto非支配解集

表9 柔性路徑及動態(tài)負載對優(yōu)化目標的影響

為驗證分路段對城市配送碳排放量的影響，計算表7優(yōu)化方案中的3個包含分路段的路徑碳排放量。當分路段時，分別計算每個路段的碳排放量并求和，當不分路段時，按照整條路徑的平均速度計算碳排放量，計算結果見表10。可以看出，城市配送中，分路段對碳排放量影響較大，尤其是各路段速度差異較大時，而同一時間，城市道路擁堵程度往往差異較大，所以在城市配送調度優(yōu)化過程中，不能忽略分路段對優(yōu)化結果的影響。

表10 分路段對優(yōu)化目標的影響

為驗證SA-ACO算法的效果，分別使用粒子群(particle swarm optimization，PSO)、SA和ACO算法，按照4.1節(jié)隨機算例數據，采用CPU3.3GHz，內存8 G，64位Windows10操作系統的計算機，SA中采用實數編碼，降溫系數取0.99，初始溫度1500，每個溫度迭代次數1000；PSO中慣性權重0.73，兩個加速因子均為1.21；ACO算法的參數取值與SA-ACO相同。各計算30次，結果取平均值。均采用均勻設計方法設置算法參數，結果見表11。

由表11可知，SA-ACO在碳排放量、總行駛時間、總行駛里程和使用運輸車輛數量方面實驗結果的數值均為最小。在計算耗時方面，由于SA-ACO混合蟻群算法在構造路徑方面計算量較大，相比較于PSO算法中尋找粒子方向的更新方式，SA-ACO計算耗時略高于PSO算法。

表11 算法效果比較

由圖7可知，與其它算法相比較，SA-ACO的收斂速度在100次迭代后，明顯優(yōu)于其它算法。

為進一步分析SA-ACO算法效率，本文采用TSPLIB中的3個經典的算例進行優(yōu)化分析，算法參數的取值與上述RTVVRP-PF模型相同，優(yōu)化結果見表12。

表12 TSP模型算法效果比較

由表12可以看出，SA-ACO對TSP以及VRP問題的優(yōu)化計算具有較好效果。

5 結束語

在城市配送RTVVRP-PF模型中，首先驗證了動態(tài)重調度和靜態(tài)調度相比較，對減少碳排放量的效果明顯。其次，提出SA-ACO算法，提高了算法效率和跳出局部最優(yōu)能力；設計了自適應精英個體繁殖策略，提高了種群優(yōu)良基因的個體數量；引入包含碳排放的多因子算子，增強信息素更新的方向性。最后，按照城市配送實際工況，隨機生成算例，通過對算例的運算結果進行分析，驗證了模型在動態(tài)重調度、柔性路徑、動態(tài)負載和分路段方面，對減少碳排放量效果顯著，同時驗證了SA-ACO算法的高效性。