狄衛(wèi)民,杜慧莉,張鵬閣
(鄭州大學(xué) 管理工程學(xué)院,河南 鄭州 450001)
在規(guī)劃配送車輛路線時(shí),既要追求經(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]雖然綜合考慮了道路擁堵、碳排放、多車型與時(shí)間窗,但僅以靜態(tài)區(qū)域來描述擁堵狀況,車輛行駛速度僅與所處的區(qū)域有關(guān),未涉及考慮時(shí)變速度的動(dòng)態(tài)GVRP。
本文在已有研究的基礎(chǔ)上,將配送服務(wù)時(shí)間劃分為若干時(shí)段,并以道路擁堵系數(shù)反映不同時(shí)段的道路擁堵狀況;分析速度、距離、載重和車輛類型對(duì)碳排放的影響,建立了以系統(tǒng)總成本最小化為目標(biāo)的非線性規(guī)劃模型。根據(jù)模型特點(diǎn),設(shè)計(jì)了頭腦風(fēng)暴優(yōu)化算法(brain storm optimization,BSO),通過算例驗(yàn)證了模型和算法的有效性。
考慮動(dòng)態(tài)擁堵的多車型GVRP可描述為:配送中心派遣多種車輛為若干客戶送貨,已知每種車輛的負(fù)載能力、啟動(dòng)成本和單位運(yùn)輸成本。車輛從配送中心出發(fā),為客戶服務(wù)后返回配送中心??蛻艟哂袝r(shí)間窗要求且服務(wù)時(shí)間已知,每個(gè)客戶僅由一臺(tái)車服務(wù)。不同時(shí)段的道路擁堵程度不同,道路擁堵影響車輛行駛速度和行駛時(shí)間。配送車輛的碳排放量與車輛自身狀況和道路擁堵程度有關(guān)。本文要解決的主要問題是:為使配送系統(tǒng)的總成本最小,應(yīng)如何確定配送車輛的類型和數(shù)量?如何安排各車輛的配送路徑?
為方便模型構(gòu)建,提出以下假設(shè):①同一時(shí)段內(nèi)車輛行駛速度恒定,不同時(shí)段的車輛行駛速度不同;②僅考慮單一產(chǎn)品的配送;③配送中心和客戶的位置已知,各節(jié)點(diǎn)間距離已知;④客戶必須全部被服務(wù),并且每個(gè)客戶只允許訪問一次;⑤客戶的需求量已知,且小于車輛的最大負(fù)載能力;⑥車輛需在客戶規(guī)定的時(shí)間窗內(nèi)完成配送任務(wù),車輛提前或者延遲到達(dá)均要承擔(dān)相應(yīng)的懲罰。
G=(N,R)為配送網(wǎng)絡(luò),N為節(jié)點(diǎn)集合,N={0,1,2,…,n},其中,0代表配送中心,N0={1,2,…n} 代表客戶;R為節(jié)點(diǎn)間的弧集,R={(i,j)|i,j∈N,i≠j};L為車型集合,L={1,2,…,l};K為某類型車輛的編號(hào)集合,K={1,2,…,k};T為配送服務(wù)時(shí)間的長(zhǎng)度,T={T1,T2,…,Th,…,TM} 表示將T劃分為M個(gè)時(shí)段;Zl為l類型車輛的數(shù)量;fl為l類型車輛的啟動(dòng)成本;cl為l類型車輛的單位運(yùn)輸成本;qi為客戶i的需求量;Ql為l類型車輛的最大負(fù)載能力;dij為弧(i,j)的長(zhǎng)度;rh為Th時(shí)段的道路擁堵系數(shù),1≤rh≤10,取值越大表示擁堵狀況越嚴(yán)重;v為車輛在道路暢通狀況下的行駛速度;λ為消耗每千克CO2的環(huán)境成本;Ai為違反客戶i規(guī)定時(shí)間窗的懲罰成本;si為客戶i的服務(wù)時(shí)間;[ETi,LTi]為客戶i規(guī)定的時(shí)間窗。
道路擁堵系數(shù)是指平均一次出行實(shí)際旅行時(shí)間與自由流狀態(tài)下旅行時(shí)間的比值,可通過百度地圖獲取。本文使用的道路擁堵系數(shù)以近七日同時(shí)段道路擁堵系數(shù)的平均值表示。由于一天內(nèi)城市交通狀況呈規(guī)律性變化,因此車輛行駛速度具有明顯的時(shí)間相關(guān)性[9]。定義車輛在時(shí)段Th=[th,t′h]內(nèi)的行駛速度為vh,則有
(1)
t′h=th+1
(2)
(3)
當(dāng)1≤p≤M-h時(shí)
(4)
(5)
(6)
(7)
(8)
當(dāng)g∈[0,p-1]時(shí)
(9)
否則,當(dāng)g∈[p+1,M-h]時(shí)
(10)
式(5)、式(7)、式(9)和式(10)計(jì)算了不同情形下車輛在每個(gè)時(shí)段內(nèi)的實(shí)際行駛距離;式(6)、式(8)計(jì)算了不同情形下車輛的行駛時(shí)間。
一般情況下,運(yùn)輸車輛在行駛過程中必然會(huì)產(chǎn)生碳排放,碳排放量與車輛行駛速度、距離、載重及車型有關(guān)。本文采用歐盟委員會(huì)在MEET報(bào)告中給出的碳排放計(jì)算函數(shù)[10]來刻畫碳排放量。
假設(shè)車輛在時(shí)段Th從節(jié)點(diǎn)i出發(fā)駛向節(jié)點(diǎn)j,途徑時(shí)段Th+g內(nèi)產(chǎn)生的碳排放量可用式(11)表示
(11)
此公式僅適用于計(jì)算空載車輛在平緩道路上行駛時(shí)產(chǎn)生的碳排放量(單位:克)。
然而,車輛在配送過程中的載重不可能完全為零,因此,Mansoureh等[11]在MEET模型的基礎(chǔ)上考慮了載重約束進(jìn)行修正。載重約束的計(jì)算公式為
(12)
因此,車輛在弧(i,j)上產(chǎn)生的碳排放量為
(13)
由上述分析,構(gòu)建的考慮動(dòng)態(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ù),表示由車輛啟動(dòng)成本、運(yùn)輸成本、碳排放成本和違反時(shí)間窗的懲罰成本構(gòu)成的系統(tǒng)總成本最小。式(15)為車輛載重約束,表示該車輛服務(wù)的客戶的總需求量不得超過車輛負(fù)載能力;式(16)保證每個(gè)客戶僅能由一輛車提供服務(wù);式(17)表示使用某類型的車輛總數(shù)不得超過該類型車輛的原有數(shù)量;式(18)確保每個(gè)客戶只服務(wù)一次;式(19)表示車輛為客戶服務(wù)后必須離開;式(20)保證車輛從配送中心駛出并在完成配送任務(wù)后返回配送中心;式(21)表示車輛在零時(shí)刻從配送中心駛出;式(22)、式(23)分別為到達(dá)和離開客戶的時(shí)間約束;式(24)~式(26)為決策變量取值約束。
考慮動(dòng)態(tài)擁堵的多車型GVRP是VRP的延伸,屬于NP-hard問題,精確算法無法避開指數(shù)爆炸的問題,只適合求解小規(guī)模問題,對(duì)大規(guī)模問題難以求得最優(yōu)解,而元啟發(fā)式算法在求解此類問題上效果較好[12]。BSO算法是在模擬人類頭腦風(fēng)暴過程的基礎(chǔ)上形成的一種群體智能算法,頭腦風(fēng)暴過程的每個(gè)想法代表解集合中的一個(gè)個(gè)體,通過聚類方法分析解集合構(gòu)成,基于解的分布生成新個(gè)體,經(jīng)過不斷迭代,得出滿意解[13]。BSO算法對(duì)初始值沒有要求,具有極強(qiáng)的全局搜索和快速收斂能力,因此,本文使用BSO算法求解考慮動(dòng)態(tài)擁堵的多車型GVRP模型。BSO算法的具體步驟如下:
步驟1 生成初始種群
步驟2 計(jì)算適應(yīng)度值
對(duì)初始種群中的每個(gè)個(gè)體以車輛載重和客戶規(guī)定的時(shí)間窗進(jìn)行檢驗(yàn)。若該車輛服務(wù)的客戶需求總量不超過車輛載重量,則此路線成立;若該車輛服務(wù)的客戶需求總量大于車輛載重量,則按客戶的服務(wù)順序依次判斷。例如,對(duì)于一條路線[1 2 3 4 5],若該車輛服務(wù)客戶3時(shí)符合載重要求,一旦服務(wù)客戶4則車輛超載,則規(guī)定該車輛在服務(wù)完客戶3后返回配送中心,放棄服務(wù)客戶4和5。對(duì)于放棄服務(wù)的每個(gè)客戶,產(chǎn)生一個(gè)極大的缺貨成本,表示未服務(wù)該客戶產(chǎn)生的懲罰。然后依次計(jì)算相關(guān)聯(lián)節(jié)點(diǎn)間的行駛時(shí)間及路段載重率,判斷是否符合客戶的時(shí)間窗要求,不符合要求的計(jì)算違反時(shí)間窗的懲罰成本。令個(gè)體的適應(yīng)度值為車輛啟動(dòng)成本、運(yùn)輸成本、碳排放成本、違反時(shí)間窗的懲罰成本和缺貨成本之和。
步驟3 聚類操作
用K-means聚類方法將種群中的個(gè)體聚成E類,并選擇該類中適應(yīng)度值最小的個(gè)體作為聚類中心。
步驟4 判斷聚類中心是否被取代
隨機(jī)產(chǎn)生(0,1)之間的數(shù),比較該值與給定概率P1的大小(0 步驟5 個(gè)體更新 在頭腦風(fēng)暴過程中,需要不斷提出新的“想法”以期找出更優(yōu)的決策方案,同理,BSO算法中使用在原個(gè)體上添加“隨機(jī)擾動(dòng)”的方式進(jìn)行個(gè)體更新。BSO算法中添加“隨機(jī)擾動(dòng)”的方式為在原個(gè)體上加入高斯隨機(jī)數(shù),然而本文采用的編碼方式為整數(shù)編碼,不適用此方式。因此,本文采用遺傳算法的交叉和變異操作作為個(gè)體更新的“隨機(jī)擾動(dòng)”。 首先,產(chǎn)生(0,1)之間的隨機(jī)數(shù)Pb,并與給定的概率P2比較(0 (1)變異操作:若Pb (2)交叉操作:若Pb≥P2,隨機(jī)選擇兩個(gè)類,并產(chǎn)生(0,1)之間的隨機(jī)數(shù)Pd。對(duì)于給定的概率P4(0 步驟6 判斷是否完成個(gè)體更新 若個(gè)體更新達(dá)到SIZE次,則完成個(gè)體更新,轉(zhuǎn)入步驟7否則返回步驟5。 步驟7 判斷是否滿足停止條件 若算法已達(dá)到最大迭代次數(shù),則算法終止,得到滿意解;否則,轉(zhuǎn)入步驟3。 BSO算法的流程如圖1所示。 圖1 BSO算法流程 本文使用一個(gè)隨機(jī)生成的仿真算例來檢驗(yàn)構(gòu)建的模型和BSO算法。假設(shè)某地區(qū)有1個(gè)編號(hào)為0的配送中心和30個(gè)客戶??蛻舻男枨罅俊r(shí)間窗及服務(wù)時(shí)間見表1;節(jié)點(diǎn)之間的距離見表2;車輛的相關(guān)參數(shù)見表3。 表1 客戶信息 表2 配送網(wǎng)絡(luò)中部分節(jié)點(diǎn)之間的距離/km 表3 車輛參數(shù) 其余數(shù)據(jù)設(shè)置如下:消耗每千克CO2的環(huán)境成本λ=0.5元;違反客戶規(guī)定時(shí)間窗的懲罰成本Ai=600元/h;未服務(wù)客戶的缺貨懲罰為10 000元/個(gè);車輛在暢通路況下的行駛速度為40 km/h;配送服務(wù)時(shí)間長(zhǎng)度為6∶00~13∶00,規(guī)定6∶00為配送服務(wù)開始的零時(shí)刻。根據(jù)城市道路擁堵情況,將配送中心的服務(wù)時(shí)間劃分為6∶00~7∶00、7∶00~9∶00、9∶00~13∶00這3個(gè)時(shí)段,定義6∶00~7∶00、9∶00~13∶00為正常時(shí)段,道路擁堵系數(shù)r1=r3=1;7∶00~9∶00為早高峰時(shí)段,道路擁堵系數(shù)r2=2。 程序采用MATLAB R2014a編程實(shí)現(xiàn),其中,在實(shí)驗(yàn)測(cè)試的基礎(chǔ)上,BSO算法的參數(shù)設(shè)置如下:種群大小為100;最大迭代次數(shù)為200;聚類數(shù)目為5;概率P1、P2、P3、P4分別為0.2、0.8、0.4、0.5。 4.2.1 算例結(jié)果分析 將BSO算法在操作系統(tǒng)為Win 10,主頻為2.7 Ghz的Intel Core i5 處理器上運(yùn)行10次,平均運(yùn)行時(shí)間為82.3 s。其中,最優(yōu)運(yùn)算結(jié)果的迭代趨勢(shì)如圖2所示。由圖2可知,在第156代求得當(dāng)前最優(yōu)值5167.3,種群的最小成本和平均成本隨著迭代次數(shù)增加均有明顯的下降趨勢(shì),兩者之間的差值逐漸減小,說明BSO算法能夠有效求解該模型并且具有良好的收斂性。 圖2 BSO算法迭代趨勢(shì) 表4列舉了圖2對(duì)應(yīng)的最優(yōu)車輛行駛路徑。結(jié)果表明,該配送中心共有7條配送路徑,其中,使用2臺(tái)類型一車輛、3臺(tái)類型二車輛、2臺(tái)類型三車輛。 表4 最優(yōu)車輛行駛路徑 4.2.2 算法有效性分析 遺傳算法(genetic algorithm,GA)已被許多學(xué)者證明可以有效求解VRP問題,且應(yīng)用較廣[12],并且BSO算法中的個(gè)體更新方式使用了遺傳算法的交叉和變異操作,于是本文使用GA與BSO算法進(jìn)行對(duì)比。在實(shí)驗(yàn)測(cè)試的基礎(chǔ)上,GA的設(shè)置如下:編碼與解碼方式與BSO算法相同,適應(yīng)度函數(shù)設(shè)置為總成本的倒數(shù),采用輪盤賭法進(jìn)行選擇,初始種群為100,最大迭代次數(shù)為200,交叉概率為0.7,變異概率為0.2。將GA程序運(yùn)行10次,運(yùn)行結(jié)果與BSO算法比較,見表5。 表5 BSO算法與GA的結(jié)果對(duì)比 由表5可知,BSO算法的最小目標(biāo)值和平均目標(biāo)值均小于GA,并且BSO算法的平均目標(biāo)值較GA降低了19.19%,兩者的運(yùn)行時(shí)間差距不大,BSO算法比GA耗時(shí)多2.11%。通過對(duì)比可以看出,BSO算法的運(yùn)算效果較好,且運(yùn)算結(jié)果優(yōu)于GA,驗(yàn)證BSO算法性能較好,可以有效求解本文研究的問題。 4.2.3 多車型與單車型配送對(duì)比 為驗(yàn)證多車型配送有利于降低成本,進(jìn)行多車型與單車型配送的對(duì)比分析。針對(duì)此算例,假設(shè)配送中心分別采用4 t、6 t和8 t的單一車型進(jìn)行配送服務(wù),除車輛相關(guān)參數(shù)外,其余設(shè)置保持不變,使用BSO算法求解。將計(jì)算結(jié)果與圖2最優(yōu)值進(jìn)行比較,見表6。其中,TC表示總成本;C1表示車輛啟動(dòng)成本;C2表示運(yùn)輸成本;C3表示碳排放成本;C4表示違反時(shí)間窗的懲罰成本(單位:元)。 表6 單車型與多車型配送成本對(duì)比 由表6結(jié)果可知:當(dāng)僅使用4 t、6 t和8 t的單一車型進(jìn)行配送時(shí),總成本比多車型配送分別高出0.97%、6.48%和15.41%。因此,使用多車型進(jìn)行配送服務(wù)是合理的。 4.2.4 考慮動(dòng)態(tài)擁堵的GVRP與傳統(tǒng)GVRP的對(duì)比 為驗(yàn)證在配送決策中考慮道路擁堵狀況有利于降低配送成本,將本文研究的問題與傳統(tǒng)GVRP進(jìn)行對(duì)比。在本文的算例中去除道路擁堵系數(shù)這一參數(shù),其余參數(shù)保持不變,于是,原問題變成車速確定的GVRP,使用BSO算法予以求解。將求得的最優(yōu)配送路徑帶入考慮動(dòng)態(tài)擁堵的GVRP,計(jì)算在未提前考慮交通狀況時(shí),一旦發(fā)生擁堵的配送成本。表7顯示了考慮動(dòng)態(tài)擁堵的GVRP與傳統(tǒng)GVRP的結(jié)果對(duì)比。其中,各參數(shù)的含義與表6相同。 表7 考慮動(dòng)態(tài)擁堵的GVRP與傳統(tǒng)GVRP的結(jié)果對(duì)比 表7數(shù)據(jù)表明,傳統(tǒng)GVRP求得的配送路徑在擁堵狀況下的各項(xiàng)成本都高于圖2結(jié)果,且總成本比考慮動(dòng)態(tài)擁堵的GVRP高出7.10%。因此,在優(yōu)化配送路徑時(shí)考慮道路擁堵可以有效減少配送費(fèi)用,更具經(jīng)濟(jì)效益。 配送車輛的碳排放會(huì)對(duì)環(huán)境造成影響,物流企業(yè)應(yīng)該引起高度重視。本文研究了考慮動(dòng)態(tài)擁堵的多車型GVRP,將道路交通狀況、碳排放、多車型及時(shí)間窗融入配送車輛路徑優(yōu)化中,使用時(shí)變速度描述道路擁堵狀況,確定了道路擁堵狀況下的車輛行駛時(shí)間,并引入碳排放計(jì)算公式,建立了以系統(tǒng)總成本最小化為目標(biāo)的數(shù)學(xué)模型。由于該模型是非線性規(guī)劃模型,根據(jù)其特點(diǎn),設(shè)計(jì)BSO算法,通過算例及與遺傳算法的對(duì)比驗(yàn)證了模型及算法的有效性。研究結(jié)果表明,在優(yōu)化配送車輛路徑時(shí)考慮動(dòng)態(tài)擁堵及車型因素有利于降低配送成本。期望本文能為企業(yè)綠色配送和碳減排提供決策參考。4 算例分析
4.1 仿真數(shù)據(jù)
4.2 結(jié)果分析
5 結(jié)束語