狄衛民,杜慧莉,張鵬閣
(鄭州大學 管理工程學院,河南 鄭州 450001)
在規劃配送車輛路線時,既要追求經濟目標,又要注重環境影響[1],因此,綠色車輛路徑問題(green vehicle routing problem,GVRP)引起了學者們的關注。Li等[2]、張亞明等[3]和段鳳華等[4]研究了多車型GVRP問題,結果表明多車型配送有利于碳減排,但是這些學者均未考慮配送過程中的道路交通狀況。近些年來,城市道路擁堵問題日益突出,Xiao等[5]、徐梅等[6]和周鮮成等[7]在GVRP問題中考慮了道路擁堵因素,但尚未涉及多車型配送的情形。趙志學等[8]雖然綜合考慮了道路擁堵、碳排放、多車型與時間窗,但僅以靜態區域來描述擁堵狀況,車輛行駛速度僅與所處的區域有關,未涉及考慮時變速度的動態GVRP。
本文在已有研究的基礎上,將配送服務時間劃分為若干時段,并以道路擁堵系數反映不同時段的道路擁堵狀況;分析速度、距離、載重和車輛類型對碳排放的影響,建立了以系統總成本最小化為目標的非線性規劃模型。根據模型特點,設計了頭腦風暴優化算法(brain storm optimization,BSO),通過算例驗證了模型和算法的有效性。
考慮動態擁堵的多車型GVRP可描述為:配送中心派遣多種車輛為若干客戶送貨,已知每種車輛的負載能力、啟動成本和單位運輸成本。車輛從配送中心出發,為客戶服務后返回配送中心。客戶具有時間窗要求且服務時間已知,每個客戶僅由一臺車服務。不同時段的道路擁堵程度不同,道路擁堵影響車輛行駛速度和行駛時間。配送車輛的碳排放量與車輛自身狀況和道路擁堵程度有關。本文要解決的主要問題是:為使配送系統的總成本最小,應如何確定配送車輛的類型和數量?如何安排各車輛的配送路徑?
為方便模型構建,提出以下假設:①同一時段內車輛行駛速度恒定,不同時段的車輛行駛速度不同;②僅考慮單一產品的配送;③配送中心和客戶的位置已知,各節點間距離已知;④客戶必須全部被服務,并且每個客戶只允許訪問一次;⑤客戶的需求量已知,且小于車輛的最大負載能力;⑥車輛需在客戶規定的時間窗內完成配送任務,車輛提前或者延遲到達均要承擔相應的懲罰。
G=(N,R)為配送網絡,N為節點集合,N={0,1,2,…,n},其中,0代表配送中心,N0={1,2,…n} 代表客戶;R為節點間的弧集,R={(i,j)|i,j∈N,i≠j};L為車型集合,L={1,2,…,l};K為某類型車輛的編號集合,K={1,2,…,k};T為配送服務時間的長度,T={T1,T2,…,Th,…,TM} 表示將T劃分為M個時段;Zl為l類型車輛的數量;fl為l類型車輛的啟動成本;cl為l類型車輛的單位運輸成本;qi為客戶i的需求量;Ql為l類型車輛的最大負載能力;dij為弧(i,j)的長度;rh為Th時段的道路擁堵系數,1≤rh≤10,取值越大表示擁堵狀況越嚴重;v為車輛在道路暢通狀況下的行駛速度;λ為消耗每千克CO2的環境成本;Ai為違反客戶i規定時間窗的懲罰成本;si為客戶i的服務時間;[ETi,LTi]為客戶i規定的時間窗。

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


(3)
當1≤p≤M-h時
(4)

(5)
(6)

(7)
(8)
當g∈[0,p-1]時
(9)
否則,當g∈[p+1,M-h]時
(10)
式(5)、式(7)、式(9)和式(10)計算了不同情形下車輛在每個時段內的實際行駛距離;式(6)、式(8)計算了不同情形下車輛的行駛時間。
一般情況下,運輸車輛在行駛過程中必然會產生碳排放,碳排放量與車輛行駛速度、距離、載重及車型有關。本文采用歐盟委員會在MEET報告中給出的碳排放計算函數[10]來刻畫碳排放量。
假設車輛在時段Th從節點i出發駛向節點j,途徑時段Th+g內產生的碳排放量可用式(11)表示


(11)
此公式僅適用于計算空載車輛在平緩道路上行駛時產生的碳排放量(單位:克)。
然而,車輛在配送過程中的載重不可能完全為零,因此,Mansoureh等[11]在MEET模型的基礎上考慮了載重約束進行修正。載重約束的計算公式為


(12)

因此,車輛在弧(i,j)上產生的碳排放量為
(13)
由上述分析,構建的考慮動態擁堵的多車型GVRP的非線性規劃模型如下



(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
xlk={0,1}, ?l∈L, ?k∈K
(24)
(25)
(26)
式(14)為目標函數,表示由車輛啟動成本、運輸成本、碳排放成本和違反時間窗的懲罰成本構成的系統總成本最小。式(15)為車輛載重約束,表示該車輛服務的客戶的總需求量不得超過車輛負載能力;式(16)保證每個客戶僅能由一輛車提供服務;式(17)表示使用某類型的車輛總數不得超過該類型車輛的原有數量;式(18)確保每個客戶只服務一次;式(19)表示車輛為客戶服務后必須離開;式(20)保證車輛從配送中心駛出并在完成配送任務后返回配送中心;式(21)表示車輛在零時刻從配送中心駛出;式(22)、式(23)分別為到達和離開客戶的時間約束;式(24)~式(26)為決策變量取值約束。
考慮動態擁堵的多車型GVRP是VRP的延伸,屬于NP-hard問題,精確算法無法避開指數爆炸的問題,只適合求解小規模問題,對大規模問題難以求得最優解,而元啟發式算法在求解此類問題上效果較好[12]。BSO算法是在模擬人類頭腦風暴過程的基礎上形成的一種群體智能算法,頭腦風暴過程的每個想法代表解集合中的一個個體,通過聚類方法分析解集合構成,基于解的分布生成新個體,經過不斷迭代,得出滿意解[13]。BSO算法對初始值沒有要求,具有極強的全局搜索和快速收斂能力,因此,本文使用BSO算法求解考慮動態擁堵的多車型GVRP模型。BSO算法的具體步驟如下:
步驟1 生成初始種群

步驟2 計算適應度值
對初始種群中的每個個體以車輛載重和客戶規定的時間窗進行檢驗。若該車輛服務的客戶需求總量不超過車輛載重量,則此路線成立;若該車輛服務的客戶需求總量大于車輛載重量,則按客戶的服務順序依次判斷。例如,對于一條路線[1 2 3 4 5],若該車輛服務客戶3時符合載重要求,一旦服務客戶4則車輛超載,則規定該車輛在服務完客戶3后返回配送中心,放棄服務客戶4和5。對于放棄服務的每個客戶,產生一個極大的缺貨成本,表示未服務該客戶產生的懲罰。然后依次計算相關聯節點間的行駛時間及路段載重率,判斷是否符合客戶的時間窗要求,不符合要求的計算違反時間窗的懲罰成本。令個體的適應度值為車輛啟動成本、運輸成本、碳排放成本、違反時間窗的懲罰成本和缺貨成本之和。
步驟3 聚類操作
用K-means聚類方法將種群中的個體聚成E類,并選擇該類中適應度值最小的個體作為聚類中心。
步驟4 判斷聚類中心是否被取代
隨機產生(0,1)之間的數,比較該值與給定概率P1的大小(0 步驟5 個體更新 在頭腦風暴過程中,需要不斷提出新的“想法”以期找出更優的決策方案,同理,BSO算法中使用在原個體上添加“隨機擾動”的方式進行個體更新。BSO算法中添加“隨機擾動”的方式為在原個體上加入高斯隨機數,然而本文采用的編碼方式為整數編碼,不適用此方式。因此,本文采用遺傳算法的交叉和變異操作作為個體更新的“隨機擾動”。 首先,產生(0,1)之間的隨機數Pb,并與給定的概率P2比較(0 (1)變異操作:若Pb (2)交叉操作:若Pb≥P2,隨機選擇兩個類,并產生(0,1)之間的隨機數Pd。對于給定的概率P4(0 步驟6 判斷是否完成個體更新 若個體更新達到SIZE次,則完成個體更新,轉入步驟7否則返回步驟5。 步驟7 判斷是否滿足停止條件 若算法已達到最大迭代次數,則算法終止,得到滿意解;否則,轉入步驟3。 BSO算法的流程如圖1所示。 圖1 BSO算法流程 本文使用一個隨機生成的仿真算例來檢驗構建的模型和BSO算法。假設某地區有1個編號為0的配送中心和30個客戶。客戶的需求量、時間窗及服務時間見表1;節點之間的距離見表2;車輛的相關參數見表3。 表1 客戶信息 表2 配送網絡中部分節點之間的距離/km 表3 車輛參數 其余數據設置如下:消耗每千克CO2的環境成本λ=0.5元;違反客戶規定時間窗的懲罰成本Ai=600元/h;未服務客戶的缺貨懲罰為10 000元/個;車輛在暢通路況下的行駛速度為40 km/h;配送服務時間長度為6∶00~13∶00,規定6∶00為配送服務開始的零時刻。根據城市道路擁堵情況,將配送中心的服務時間劃分為6∶00~7∶00、7∶00~9∶00、9∶00~13∶00這3個時段,定義6∶00~7∶00、9∶00~13∶00為正常時段,道路擁堵系數r1=r3=1;7∶00~9∶00為早高峰時段,道路擁堵系數r2=2。 程序采用MATLAB R2014a編程實現,其中,在實驗測試的基礎上,BSO算法的參數設置如下:種群大小為100;最大迭代次數為200;聚類數目為5;概率P1、P2、P3、P4分別為0.2、0.8、0.4、0.5。 4.2.1 算例結果分析 將BSO算法在操作系統為Win 10,主頻為2.7 Ghz的Intel Core i5 處理器上運行10次,平均運行時間為82.3 s。其中,最優運算結果的迭代趨勢如圖2所示。由圖2可知,在第156代求得當前最優值5167.3,種群的最小成本和平均成本隨著迭代次數增加均有明顯的下降趨勢,兩者之間的差值逐漸減小,說明BSO算法能夠有效求解該模型并且具有良好的收斂性。 圖2 BSO算法迭代趨勢 表4列舉了圖2對應的最優車輛行駛路徑。結果表明,該配送中心共有7條配送路徑,其中,使用2臺類型一車輛、3臺類型二車輛、2臺類型三車輛。 表4 最優車輛行駛路徑 4.2.2 算法有效性分析 遺傳算法(genetic algorithm,GA)已被許多學者證明可以有效求解VRP問題,且應用較廣[12],并且BSO算法中的個體更新方式使用了遺傳算法的交叉和變異操作,于是本文使用GA與BSO算法進行對比。在實驗測試的基礎上,GA的設置如下:編碼與解碼方式與BSO算法相同,適應度函數設置為總成本的倒數,采用輪盤賭法進行選擇,初始種群為100,最大迭代次數為200,交叉概率為0.7,變異概率為0.2。將GA程序運行10次,運行結果與BSO算法比較,見表5。 表5 BSO算法與GA的結果對比 由表5可知,BSO算法的最小目標值和平均目標值均小于GA,并且BSO算法的平均目標值較GA降低了19.19%,兩者的運行時間差距不大,BSO算法比GA耗時多2.11%。通過對比可以看出,BSO算法的運算效果較好,且運算結果優于GA,驗證BSO算法性能較好,可以有效求解本文研究的問題。 4.2.3 多車型與單車型配送對比 為驗證多車型配送有利于降低成本,進行多車型與單車型配送的對比分析。針對此算例,假設配送中心分別采用4 t、6 t和8 t的單一車型進行配送服務,除車輛相關參數外,其余設置保持不變,使用BSO算法求解。將計算結果與圖2最優值進行比較,見表6。其中,TC表示總成本;C1表示車輛啟動成本;C2表示運輸成本;C3表示碳排放成本;C4表示違反時間窗的懲罰成本(單位:元)。 表6 單車型與多車型配送成本對比 由表6結果可知:當僅使用4 t、6 t和8 t的單一車型進行配送時,總成本比多車型配送分別高出0.97%、6.48%和15.41%。因此,使用多車型進行配送服務是合理的。 4.2.4 考慮動態擁堵的GVRP與傳統GVRP的對比 為驗證在配送決策中考慮道路擁堵狀況有利于降低配送成本,將本文研究的問題與傳統GVRP進行對比。在本文的算例中去除道路擁堵系數這一參數,其余參數保持不變,于是,原問題變成車速確定的GVRP,使用BSO算法予以求解。將求得的最優配送路徑帶入考慮動態擁堵的GVRP,計算在未提前考慮交通狀況時,一旦發生擁堵的配送成本。表7顯示了考慮動態擁堵的GVRP與傳統GVRP的結果對比。其中,各參數的含義與表6相同。 表7 考慮動態擁堵的GVRP與傳統GVRP的結果對比 表7數據表明,傳統GVRP求得的配送路徑在擁堵狀況下的各項成本都高于圖2結果,且總成本比考慮動態擁堵的GVRP高出7.10%。因此,在優化配送路徑時考慮道路擁堵可以有效減少配送費用,更具經濟效益。 配送車輛的碳排放會對環境造成影響,物流企業應該引起高度重視。本文研究了考慮動態擁堵的多車型GVRP,將道路交通狀況、碳排放、多車型及時間窗融入配送車輛路徑優化中,使用時變速度描述道路擁堵狀況,確定了道路擁堵狀況下的車輛行駛時間,并引入碳排放計算公式,建立了以系統總成本最小化為目標的數學模型。由于該模型是非線性規劃模型,根據其特點,設計BSO算法,通過算例及與遺傳算法的對比驗證了模型及算法的有效性。研究結果表明,在優化配送車輛路徑時考慮動態擁堵及車型因素有利于降低配送成本。期望本文能為企業綠色配送和碳減排提供決策參考。
4 算例分析
4.1 仿真數據



4.2 結果分析





5 結束語