狄衛(wèi)民，杜慧莉，張鵬閣

(鄭州大學(xué) 管理工程學(xué)院，河南 鄭州 450001)

0 引 言

在規(guī)劃配送車輛路線時，既要追求經(jīng)濟(jì)目標(biāo)，又要注重環(huán)境影響[1]，因此，綠色車輛路徑問題(green vehicle routing problem，GVRP)引起了學(xué)者們的關(guān)注。Li等[2]、張亞明等[3]和段鳳華等[4]研究了多車型GVRP問題，結(jié)果表明多車型配送有利于碳減排，但是這些學(xué)者均未考慮配送過程中的道路交通狀況。近些年來，城市道路擁堵問題日益突出，Xiao等[5]、徐梅等[6]和周鮮成等[7]在GVRP問題中考慮了道路擁堵因素，但尚未涉及多車型配送的情形。趙志學(xué)等[8]雖然綜合考慮了道路擁堵、碳排放、多車型與時間窗，但僅以靜態(tài)區(qū)域來描述擁堵狀況，車輛行駛速度僅與所處的區(qū)域有關(guān)，未涉及考慮時變速度的動態(tài)GVRP。

本文在已有研究的基礎(chǔ)上，將配送服務(wù)時間劃分為若干時段，并以道路擁堵系數(shù)反映不同時段的道路擁堵狀況；分析速度、距離、載重和車輛類型對碳排放的影響，建立了以系統(tǒng)總成本最小化為目標(biāo)的非線性規(guī)劃模型。根據(jù)模型特點，設(shè)計了頭腦風(fēng)暴優(yōu)化算法(brain storm optimization，BSO)，通過算例驗證了模型和算法的有效性。

1 問題描述與假設(shè)

1.1 問題描述

考慮動態(tài)擁堵的多車型GVRP可描述為：配送中心派遣多種車輛為若干客戶送貨，已知每種車輛的負(fù)載能力、啟動成本和單位運輸成本。車輛從配送中心出發(fā)，為客戶服務(wù)后返回配送中心。客戶具有時間窗要求且服務(wù)時間已知，每個客戶僅由一臺車服務(wù)。不同時段的道路擁堵程度不同，道路擁堵影響車輛行駛速度和行駛時間。配送車輛的碳排放量與車輛自身狀況和道路擁堵程度有關(guān)。本文要解決的主要問題是：為使配送系統(tǒng)的總成本最小，應(yīng)如何確定配送車輛的類型和數(shù)量？如何安排各車輛的配送路徑？

1.2 假設(shè)條件

為方便模型構(gòu)建，提出以下假設(shè)：①同一時段內(nèi)車輛行駛速度恒定，不同時段的車輛行駛速度不同；②僅考慮單一產(chǎn)品的配送；③配送中心和客戶的位置已知，各節(jié)點間距離已知；④客戶必須全部被服務(wù)，并且每個客戶只允許訪問一次；⑤客戶的需求量已知，且小于車輛的最大負(fù)載能力；⑥車輛需在客戶規(guī)定的時間窗內(nèi)完成配送任務(wù)，車輛提前或者延遲到達(dá)均要承擔(dān)相應(yīng)的懲罰。

2 模型構(gòu)建

2.1 符號與變量

G=(N,R)為配送網(wǎng)絡(luò)，N為節(jié)點集合，N={0,1,2,…,n}，其中，0代表配送中心，N0={1,2,…n} 代表客戶；R為節(jié)點間的弧集，R={(i,j)|i,j∈N,i≠j}；L為車型集合，L={1,2,…,l}；K為某類型車輛的編號集合，K={1,2,…,k}；T為配送服務(wù)時間的長度，T={T1,T2,…,Th,…,TM} 表示將T劃分為M個時段；Zl為l類型車輛的數(shù)量；fl為l類型車輛的啟動成本；cl為l類型車輛的單位運輸成本；qi為客戶i的需求量；Ql為l類型車輛的最大負(fù)載能力；dij為弧(i,j)的長度；rh為Th時段的道路擁堵系數(shù)，1≤rh≤10，取值越大表示擁堵狀況越嚴(yán)重；v為車輛在道路暢通狀況下的行駛速度；λ為消耗每千克CO2的環(huán)境成本；Ai為違反客戶i規(guī)定時間窗的懲罰成本；si為客戶i的服務(wù)時間；[ETi,LTi]為客戶i規(guī)定的時間窗。

2.2 擁堵狀況下車輛行駛時間分析

道路擁堵系數(shù)是指平均一次出行實際旅行時間與自由流狀態(tài)下旅行時間的比值，可通過百度地圖獲取。本文使用的道路擁堵系數(shù)以近七日同時段道路擁堵系數(shù)的平均值表示。由于一天內(nèi)城市交通狀況呈規(guī)律性變化，因此車輛行駛速度具有明顯的時間相關(guān)性[9]。定義車輛在時段Th=[th,t′h]內(nèi)的行駛速度為vh，則有

(1)

t′h=th+1

(2)

(3)

當(dāng)1≤p≤M-h時

(4)

(5)

(6)

(7)

(8)

當(dāng)g∈[0,p-1]時

(9)

否則，當(dāng)g∈[p+1,M-h]時

(10)

式(5)、式(7)、式(9)和式(10)計算了不同情形下車輛在每個時段內(nèi)的實際行駛距離；式(6)、式(8)計算了不同情形下車輛的行駛時間。

2.3 碳排放量的計算

一般情況下，運輸車輛在行駛過程中必然會產(chǎn)生碳排放，碳排放量與車輛行駛速度、距離、載重及車型有關(guān)。本文采用歐盟委員會在MEET報告中給出的碳排放計算函數(shù)[10]來刻畫碳排放量。

假設(shè)車輛在時段Th從節(jié)點i出發(fā)駛向節(jié)點j，途徑時段Th+g內(nèi)產(chǎn)生的碳排放量可用式(11)表示

(11)

此公式僅適用于計算空載車輛在平緩道路上行駛時產(chǎn)生的碳排放量(單位：克)。

然而，車輛在配送過程中的載重不可能完全為零，因此，Mansoureh等[11]在MEET模型的基礎(chǔ)上考慮了載重約束進(jìn)行修正。載重約束的計算公式為

(12)

因此，車輛在弧(i,j)上產(chǎn)生的碳排放量為

(13)

2.4 數(shù)學(xué)模型

由上述分析，構(gòu)建的考慮動態(tài)擁堵的多車型GVRP的非線性規(guī)劃模型如下

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

xlk={0,1}, ?l∈L, ?k∈K

(24)

(25)

(26)

式(14)為目標(biāo)函數(shù)，表示由車輛啟動成本、運輸成本、碳排放成本和違反時間窗的懲罰成本構(gòu)成的系統(tǒng)總成本最小。式(15)為車輛載重約束，表示該車輛服務(wù)的客戶的總需求量不得超過車輛負(fù)載能力；式(16)保證每個客戶僅能由一輛車提供服務(wù)；式(17)表示使用某類型的車輛總數(shù)不得超過該類型車輛的原有數(shù)量；式(18)確保每個客戶只服務(wù)一次；式(19)表示車輛為客戶服務(wù)后必須離開；式(20)保證車輛從配送中心駛出并在完成配送任務(wù)后返回配送中心；式(21)表示車輛在零時刻從配送中心駛出；式(22)、式(23)分別為到達(dá)和離開客戶的時間約束；式(24)～式(26)為決策變量取值約束。

3 算法設(shè)計

考慮動態(tài)擁堵的多車型GVRP是VRP的延伸，屬于NP-hard問題，精確算法無法避開指數(shù)爆炸的問題，只適合求解小規(guī)模問題，對大規(guī)模問題難以求得最優(yōu)解，而元啟發(fā)式算法在求解此類問題上效果較好[12]。BSO算法是在模擬人類頭腦風(fēng)暴過程的基礎(chǔ)上形成的一種群體智能算法，頭腦風(fēng)暴過程的每個想法代表解集合中的一個個體，通過聚類方法分析解集合構(gòu)成，基于解的分布生成新個體，經(jīng)過不斷迭代，得出滿意解[13]。BSO算法對初始值沒有要求，具有極強(qiáng)的全局搜索和快速收斂能力，因此，本文使用BSO算法求解考慮動態(tài)擁堵的多車型GVRP模型。BSO算法的具體步驟如下：

步驟1 生成初始種群

步驟2 計算適應(yīng)度值

對初始種群中的每個個體以車輛載重和客戶規(guī)定的時間窗進(jìn)行檢驗。若該車輛服務(wù)的客戶需求總量不超過車輛載重量，則此路線成立；若該車輛服務(wù)的客戶需求總量大于車輛載重量，則按客戶的服務(wù)順序依次判斷。例如，對于一條路線[1 2 3 4 5]，若該車輛服務(wù)客戶3時符合載重要求，一旦服務(wù)客戶4則車輛超載，則規(guī)定該車輛在服務(wù)完客戶3后返回配送中心，放棄服務(wù)客戶4和5。對于放棄服務(wù)的每個客戶，產(chǎn)生一個極大的缺貨成本，表示未服務(wù)該客戶產(chǎn)生的懲罰。然后依次計算相關(guān)聯(lián)節(jié)點間的行駛時間及路段載重率，判斷是否符合客戶的時間窗要求，不符合要求的計算違反時間窗的懲罰成本。令個體的適應(yīng)度值為車輛啟動成本、運輸成本、碳排放成本、違反時間窗的懲罰成本和缺貨成本之和。

步驟3 聚類操作

用K-means聚類方法將種群中的個體聚成E類，并選擇該類中適應(yīng)度值最小的個體作為聚類中心。

步驟4 判斷聚類中心是否被取代

隨機(jī)產(chǎn)生(0,1)之間的數(shù)，比較該值與給定概率P1的大小(0

步驟5 個體更新

在頭腦風(fēng)暴過程中，需要不斷提出新的“想法”以期找出更優(yōu)的決策方案，同理，BSO算法中使用在原個體上添加“隨機(jī)擾動”的方式進(jìn)行個體更新。BSO算法中添加“隨機(jī)擾動”的方式為在原個體上加入高斯隨機(jī)數(shù)，然而本文采用的編碼方式為整數(shù)編碼，不適用此方式。因此，本文采用遺傳算法的交叉和變異操作作為個體更新的“隨機(jī)擾動”。

首先，產(chǎn)生(0,1)之間的隨機(jī)數(shù)Pb，并與給定的概率P2比較(0

(1)變異操作：若Pb

(2)交叉操作：若Pb≥P2，隨機(jī)選擇兩個類，并產(chǎn)生(0,1)之間的隨機(jī)數(shù)Pd。對于給定的概率P4(0

步驟6 判斷是否完成個體更新

若個體更新達(dá)到SIZE次，則完成個體更新，轉(zhuǎn)入步驟7否則返回步驟5。

步驟7 判斷是否滿足停止條件

若算法已達(dá)到最大迭代次數(shù)，則算法終止，得到滿意解；否則，轉(zhuǎn)入步驟3。

BSO算法的流程如圖1所示。

圖1 BSO算法流程

4 算例分析

4.1 仿真數(shù)據(jù)

本文使用一個隨機(jī)生成的仿真算例來檢驗構(gòu)建的模型和BSO算法。假設(shè)某地區(qū)有1個編號為0的配送中心和30個客戶。客戶的需求量、時間窗及服務(wù)時間見表1；節(jié)點之間的距離見表2；車輛的相關(guān)參數(shù)見表3。

表1 客戶信息

表2 配送網(wǎng)絡(luò)中部分節(jié)點之間的距離/km

表3 車輛參數(shù)

其余數(shù)據(jù)設(shè)置如下：消耗每千克CO2的環(huán)境成本λ=0.5元；違反客戶規(guī)定時間窗的懲罰成本Ai=600元/h；未服務(wù)客戶的缺貨懲罰為10 000元/個；車輛在暢通路況下的行駛速度為40 km/h；配送服務(wù)時間長度為6∶00～13∶00，規(guī)定6∶00為配送服務(wù)開始的零時刻。根據(jù)城市道路擁堵情況，將配送中心的服務(wù)時間劃分為6∶00～7∶00、7∶00～9∶00、9∶00～13∶00這3個時段，定義6∶00～7∶00、9∶00～13∶00為正常時段，道路擁堵系數(shù)r1=r3=1；7∶00～9∶00為早高峰時段，道路擁堵系數(shù)r2=2。

程序采用MATLAB R2014a編程實現(xiàn)，其中，在實驗測試的基礎(chǔ)上，BSO算法的參數(shù)設(shè)置如下：種群大小為100；最大迭代次數(shù)為200；聚類數(shù)目為5；概率P1、P2、P3、P4分別為0.2、0.8、0.4、0.5。

4.2 結(jié)果分析

4.2.1 算例結(jié)果分析

將BSO算法在操作系統(tǒng)為Win 10，主頻為2.7 Ghz的Intel Core i5 處理器上運行10次，平均運行時間為82.3 s。其中，最優(yōu)運算結(jié)果的迭代趨勢如圖2所示。由圖2可知，在第156代求得當(dāng)前最優(yōu)值5167.3，種群的最小成本和平均成本隨著迭代次數(shù)增加均有明顯的下降趨勢，兩者之間的差值逐漸減小，說明BSO算法能夠有效求解該模型并且具有良好的收斂性。

圖2 BSO算法迭代趨勢

表4列舉了圖2對應(yīng)的最優(yōu)車輛行駛路徑。結(jié)果表明，該配送中心共有7條配送路徑，其中，使用2臺類型一車輛、3臺類型二車輛、2臺類型三車輛。

表4 最優(yōu)車輛行駛路徑

4.2.2 算法有效性分析

遺傳算法(genetic algorithm，GA)已被許多學(xué)者證明可以有效求解VRP問題，且應(yīng)用較廣[12]，并且BSO算法中的個體更新方式使用了遺傳算法的交叉和變異操作，于是本文使用GA與BSO算法進(jìn)行對比。在實驗測試的基礎(chǔ)上，GA的設(shè)置如下：編碼與解碼方式與BSO算法相同，適應(yīng)度函數(shù)設(shè)置為總成本的倒數(shù)，采用輪盤賭法進(jìn)行選擇，初始種群為100，最大迭代次數(shù)為200，交叉概率為0.7，變異概率為0.2。將GA程序運行10次，運行結(jié)果與BSO算法比較，見表5。

表5 BSO算法與GA的結(jié)果對比

由表5可知，BSO算法的最小目標(biāo)值和平均目標(biāo)值均小于GA，并且BSO算法的平均目標(biāo)值較GA降低了19.19%，兩者的運行時間差距不大，BSO算法比GA耗時多2.11%。通過對比可以看出，BSO算法的運算效果較好，且運算結(jié)果優(yōu)于GA，驗證BSO算法性能較好，可以有效求解本文研究的問題。

4.2.3 多車型與單車型配送對比

為驗證多車型配送有利于降低成本，進(jìn)行多車型與單車型配送的對比分析。針對此算例，假設(shè)配送中心分別采用4 t、6 t和8 t的單一車型進(jìn)行配送服務(wù)，除車輛相關(guān)參數(shù)外，其余設(shè)置保持不變，使用BSO算法求解。將計算結(jié)果與圖2最優(yōu)值進(jìn)行比較，見表6。其中，TC表示總成本；C1表示車輛啟動成本；C2表示運輸成本；C3表示碳排放成本；C4表示違反時間窗的懲罰成本(單位：元)。

表6 單車型與多車型配送成本對比

由表6結(jié)果可知：當(dāng)僅使用4 t、6 t和8 t的單一車型進(jìn)行配送時，總成本比多車型配送分別高出0.97%、6.48%和15.41%。因此，使用多車型進(jìn)行配送服務(wù)是合理的。

4.2.4 考慮動態(tài)擁堵的GVRP與傳統(tǒng)GVRP的對比

為驗證在配送決策中考慮道路擁堵狀況有利于降低配送成本，將本文研究的問題與傳統(tǒng)GVRP進(jìn)行對比。在本文的算例中去除道路擁堵系數(shù)這一參數(shù)，其余參數(shù)保持不變，于是，原問題變成車速確定的GVRP，使用BSO算法予以求解。將求得的最優(yōu)配送路徑帶入考慮動態(tài)擁堵的GVRP，計算在未提前考慮交通狀況時，一旦發(fā)生擁堵的配送成本。表7顯示了考慮動態(tài)擁堵的GVRP與傳統(tǒng)GVRP的結(jié)果對比。其中，各參數(shù)的含義與表6相同。

表7 考慮動態(tài)擁堵的GVRP與傳統(tǒng)GVRP的結(jié)果對比

表7數(shù)據(jù)表明，傳統(tǒng)GVRP求得的配送路徑在擁堵狀況下的各項成本都高于圖2結(jié)果，且總成本比考慮動態(tài)擁堵的GVRP高出7.10%。因此，在優(yōu)化配送路徑時考慮道路擁堵可以有效減少配送費用，更具經(jīng)濟(jì)效益。

5 結(jié)束語

配送車輛的碳排放會對環(huán)境造成影響，物流企業(yè)應(yīng)該引起高度重視。本文研究了考慮動態(tài)擁堵的多車型GVRP，將道路交通狀況、碳排放、多車型及時間窗融入配送車輛路徑優(yōu)化中，使用時變速度描述道路擁堵狀況，確定了道路擁堵狀況下的車輛行駛時間，并引入碳排放計算公式，建立了以系統(tǒng)總成本最小化為目標(biāo)的數(shù)學(xué)模型。由于該模型是非線性規(guī)劃模型，根據(jù)其特點，設(shè)計BSO算法，通過算例及與遺傳算法的對比驗證了模型及算法的有效性。研究結(jié)果表明，在優(yōu)化配送車輛路徑時考慮動態(tài)擁堵及車型因素有利于降低配送成本。期望本文能為企業(yè)綠色配送和碳減排提供決策參考。