馬檢, 靳文舟

( 華南理工大學(xué) 土木與交通學(xué)院, 廣東 廣州 510641)

0 引言

近年來,我國的危險品物流產(chǎn)業(yè)維持了高速增長的態(tài)勢,但無論是在理論研究還是生產(chǎn)實踐上,都仍處于初級階段,尚有很多待解決的問題。目前學(xué)界關(guān)于危險品物流方面的研究,如配送中心選址、配送路線優(yōu)化等,都還停留在理論階段，沒有與生產(chǎn)實際進(jìn)行進(jìn)一步的結(jié)合。而生產(chǎn)中的各個危險品供應(yīng)鏈上的企業(yè)各自為政,沿用著傳統(tǒng)普通貨物的物流系統(tǒng)進(jìn)行危險品的運輸。這種方式不僅不能有效降低物流的成本,還因缺乏安全管理意識和能力而引發(fā)危險品泄露、燃燒甚至爆炸的事故。如2018年發(fā)生的致22人死亡的石家莊危險品運輸車爆炸、2019年發(fā)生的江蘇響水化工廠爆炸等事故,就是物流中不同環(huán)節(jié)安全管理不善的直接后果，為群眾生命財產(chǎn)安全、政治穩(wěn)定、社會經(jīng)濟(jì)、環(huán)境保護(hù)等方面都造成了很大的負(fù)面影響，因此,危險品物流理論和產(chǎn)業(yè)的相關(guān)研究日益為管理部門、企業(yè)界及學(xué)術(shù)界所關(guān)注。

早期對危險品物流的研究主要集中在配送運輸路徑優(yōu)化上,Dantzig等[1]于1959年首先提出,作為車輛調(diào)度問題的延伸。目前,對危險品運輸?shù)难芯科毡橐褟脑缒甑膯文繕?biāo)問題轉(zhuǎn)變?yōu)殡p目標(biāo)乃至多目標(biāo)問題,成本、風(fēng)險、風(fēng)險公平性[2]、客戶滿意度[3]等均成為研究中的目標(biāo)。精確算法[4]、禁忌搜索算法[5]、粒子群算法[6]、多目標(biāo)遺傳算法[7-9]等求解方法也被應(yīng)用于求解。其中,由于運輸路徑問題的NP-hard問題特性和遺傳算法良好的適用性,啟發(fā)式算法中的多目標(biāo)遺傳算法應(yīng)用最為廣泛。

在模型上,越來越多的約束條件被加入到研究中，如時間窗[10]、風(fēng)險公平性[11]、多種類[12]等。朱偉等[3]根據(jù)隨機(jī)需求的概率密度函數(shù),將隨機(jī)需求轉(zhuǎn)化為一個約束,從而消除了模型的不確定性。風(fēng)險的表征方面,Toumazis等[13]將在險價值理論應(yīng)用在危險品風(fēng)險模型上(CVaR)。代存杰等[14]通過事故發(fā)生概率、影響半徑和事故發(fā)生點的人口計算風(fēng)險,而考慮了風(fēng)險分布特征。柴獲等[15]將事故發(fā)生點和人口聚集點分離開來,并定義沖擊系數(shù),使風(fēng)險大小與兩點間歐氏距離有關(guān)。上述過往研究中的風(fēng)險模型距離現(xiàn)實仍有一定距離。

本文將在產(chǎn)地-配送中心-客戶三級供應(yīng)鏈系統(tǒng)中,從最優(yōu)庫存模型出發(fā),對需求的不確定性進(jìn)行處理,提出應(yīng)對一條路線上多客戶需求不確定的配送和倉儲量公式，并研究危險品風(fēng)險表征中影響半徑和衰減系數(shù)的關(guān)系,提出計算風(fēng)險的更為一般化、可實踐的方法。最后開發(fā)并測試種群結(jié)構(gòu)控制選擇算子(PSC),提出適用于其他多目標(biāo)整數(shù)規(guī)劃問題的改進(jìn)型算法。

1 模型建立

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

本文以不確定需求的多種類危險品物流管理的配送中心定位-配送中心及需求點的庫存策略制定-配送路線設(shè)計為研究對象。每種危險品對應(yīng)一個產(chǎn)地和若干備選配送中心、客戶,而每個產(chǎn)地、配送中心和客戶也對應(yīng)一種危險品,危險品從產(chǎn)地運輸?shù)脚渌椭行牡倪^程稱補(bǔ)貨,從配送中心運輸?shù)娇蛻舻倪^程稱配送。因為產(chǎn)地和危險品一一對應(yīng),配送中心和補(bǔ)貨車輛一一對應(yīng),所以各用相同符號表示。

優(yōu)化的雙目標(biāo)為總成本及風(fēng)險值的最小化。總成本包括配送中心開放成本、危險品的儲存成本、訂貨成本和運輸成本,與危險品種類、庫存計劃、運輸路徑等有關(guān)?？傦L(fēng)險考慮的是危險品的爆炸風(fēng)險,根據(jù)危險品儲存和運輸時發(fā)生爆炸的可能性、爆炸強(qiáng)度和衰減特性等來計算,與危險品特性、庫存計劃、運輸路徑、車輛速度、人口分布等有關(guān)。

決策變量Wm=1時配送中心m開放,否則為0。Xi,k=1時車輛k服務(wù)于節(jié)點i,否則為0,i,j∈M∪N。Yi,j,k=1時車輛k先后在節(jié)點i,j作業(yè),否則為0,i,j∈M∪N。

1.2 配送策略分析

(1)

(2)

令總成本最小,求得客戶的訂貨周期和最優(yōu)訂貨批量分別為

(3)

(4)

(5)

取其期望值為車輛k實際的配送周期Pk，即

(6)

(7)

設(shè)置信度α的分位數(shù)為zα,上式可寫作

(8)

(9)

(10)

(11)

(12)

則總補(bǔ)貨運輸CR為

(13)

配送中心和客戶的平均庫存成本CI為

(14)

另外,設(shè)每個配送中心的開放成本為cm,則總開放成本為

(15)

1.3 風(fēng)險模型

定義危險品對平面任意坐標(biāo)的造成的風(fēng)險與種類、距離、受影響人口有關(guān)。根據(jù)(Carotenuto等, 2007),給定危險品坐標(biāo)(x0,y0),則其影響衰減公式為e-τl[(x-x0)2+(y-y0)2],其中τl為種類l危險品的影響衰減系數(shù)。對上式在整個平面上進(jìn)行二重積分,可求得危險品的全部影響為

(16)

庫存風(fēng)險空間分布如圖1所示。若給定一個以險品坐標(biāo)(x0,y0)為圓心的圓,分布在該圓形區(qū)域內(nèi)的影響大小與全部影響之比大于一定比例β(如99%),則認(rèn)為影響全部分布在該區(qū)域中。此時可求得該圓形區(qū)域的半徑,即為危險品l的影響半徑λl,與影響衰減系數(shù)有關(guān)。

(17)

圖1 庫存風(fēng)險空間分布示意圖Fig.1 Impact distribution of hazmat’s inventory

與庫存成本公式相似,配送中心和客戶庫存總風(fēng)險RI為

(18)

在運輸中,車輛就是圖1中鐘型風(fēng)險場的中心,這個風(fēng)險場隨車輛移動,擦過道路旁的某個點(x,y),那么點(x,y)的風(fēng)險為隨時間變化的曲線,與風(fēng)險空間分布的鐘型的豎直截面具有相同形狀。運輸風(fēng)險模型的構(gòu)建如圖2所示，運輸風(fēng)險的形成如圖2(a)所示。

e-τl(d2+z2)= e-τl(d2+v2t2)。

(19)

(20)

(21)

補(bǔ)貨運輸風(fēng)險為:(l既表示危險品種類,也表示產(chǎn)地節(jié)點)

(22)

模型為

min(CD+CR+CI+CL),

(23)

min(RD+RR+RI),

(24)

s.t.

(37)

式(22)表示總成本最小化,式(24)表示總風(fēng)險最小化,式(25)表示每種危險品至少開放一個配送中心,式(26)表示未開放的配送中心沒有車輛,開放的配送中心至少有1臺車輛，式(27)、(28)表示出發(fā)的車輛至少經(jīng)過一個配送中心和一個客戶，式(29)表示每個客戶都有車輛服務(wù)，式(30)表示每個客戶只被服務(wù)一次,且車輛形成閉合回路，式(31)表示同一車輛的配送中心和客戶危險品種類相同，式(32)、(33)、(34)分別是配送中心、客戶和車輛容量約束，式(35)、(36)、(37)為決策變量。

2 算法與算例分析

2.1 編碼、雜交和變異方法

首先采用NSGA-Ⅱ算法,即帶有精英選擇策略的非支配排序多目標(biāo)遺傳算法(non-dominated sorting genetic algorithm-Ⅱ，NSGA-Ⅱ),嘗試求解該問題。

編碼方式:將所有配送車輛路線排列成2行矩陣,第1行為車輛經(jīng)過的節(jié)點編號,第2行為該節(jié)點所屬的危險品種類。每輛車的起終點均為配送中心,在個體矩陣中出現(xiàn)的配送中心則說明已啟用,未出現(xiàn)則說明不啟用。

雜交方式:

步驟①:將個體1所有客戶序列提取出來,每種危險品保留前一半序列,得到半成品序列。

步驟②:剔除個體2中半成品序列已出現(xiàn)過的客戶,分種類將個體2剩余客戶依次填入半成品的后半段。

步驟③:將合成的新序列回填個體1空位,得到新的子代1。

編碼及雜交方式演示如圖3所示。

圖3 編碼及雜交方式演示Fig.3 Coding and crossover operator

變異方式1:隨機(jī)選擇同一車輛配送的2個客戶,互換位置；

變異方式2:隨機(jī)選擇一個客戶,移到另一輛同類車輛路線上；

變異方式3:隨機(jī)取消一輛車,所服務(wù)的客戶移到其他同類車輛路線上(如果配送中心上只有這一輛車,則關(guān)閉這個配送中心)；

變異方式4:隨機(jī)選擇一個配送中心生成一輛車,從其他同類車輛中隨機(jī)抽取客戶。

2.2 求解效果與分析

創(chuàng)建一個3類危險品-每種危險品3個候選配送中心-共30客戶的算例,其路網(wǎng)及節(jié)點性質(zhì)如圖4所示:

圖4 算例路網(wǎng)Fig.4 Road network of the example

運用NSGA-Ⅱ算法求解該算例,種群規(guī)模為200,迭代200次,迭代過程及結(jié)果如圖5所示。

其中,一個非支配等級中互不相同的個體數(shù)量稱為該等級的有效規(guī)模?？梢钥闯?隨著迭代的進(jìn)行,最高等級規(guī)模迅速增大,到70次迭代時已基本占據(jù)整個種群,同時最高等級的有效規(guī)模的增長也驟然減慢,風(fēng)險和成本的降低趨于平緩,可以認(rèn)為算法已達(dá)到收斂。此后雖有優(yōu)化,但非常緩慢。最終等級1的規(guī)模為200,而Pareto解集的規(guī)模僅為30。

在傳統(tǒng)NSGA-Ⅱ中,父代和雜交變異產(chǎn)生的子代會進(jìn)行合并,合并種群根據(jù)非支配等級和擁擠度進(jìn)行排序,排名在預(yù)設(shè)種群規(guī)模以內(nèi)的個體都會被選中,組成新一代種群。假設(shè)父代中的一個最高等級個體(稱個體1)經(jīng)過雜交變異后,仍存在于子代中,那么個體1就會在合并種群中出現(xiàn)2次,而精英選擇會將這2個相同的個體1都保留下來。這種情況的連續(xù)出現(xiàn)會讓個體1數(shù)量指數(shù)增大,而更多的相同個體1意味著下一代中出現(xiàn)個體1的概率越大——2個個體1雜交后仍是個體1。這種數(shù)量和出現(xiàn)概率的互相促進(jìn),會讓個體1在短時間內(nèi)大量復(fù)制,占據(jù)大量空間,排擠掉低等級的個體。顯然,在圖5(a)中,這種最高等級規(guī)模急速增長的現(xiàn)象出現(xiàn)了數(shù)次。

這種現(xiàn)象有利有弊：一方面,加速了算法的收斂；另一方面,僅靠少量重復(fù)的優(yōu)秀個體進(jìn)行雜交變異,此后算法對解空間的探索能力也會大大下降。

這種現(xiàn)象在用整數(shù)規(guī)劃模型研究的多目標(biāo)實際問題中影響更大。非整數(shù)規(guī)劃通常有無窮多解,雜交變異后的子代中,出現(xiàn)與父代相同個體的概率很低。而整數(shù)規(guī)劃中可行解的總數(shù)可能是很有限的。并且,在對實際問題的研究中,研究者要根據(jù)問題來設(shè)計個體的編碼、解碼及種群的初始化、變異和遺傳,特別是在種群初始化時,常常不能保證對潛在解空間的覆蓋程度,進(jìn)一步減少了算法實際運行時可行解的總數(shù),導(dǎo)致算法在合并種群和精英選擇的作用下提前收斂。

對此,很有必要將種群中各等級個體的數(shù)量控制在合理的范圍內(nèi),控制該現(xiàn)象,才能平衡算法收斂速度和對解空間的覆蓋度。基于此,本文提出帶種群結(jié)構(gòu)控制的非支配排序遺傳算法(non-dominated sorting genetic algorithm-Ⅱ with population structure control， NSGA-Ⅱ-PSC),在第三節(jié)中闡述。

3 基于種群結(jié)構(gòu)控制的選擇策略改進(jìn)

3.1 選擇算子模型

(38)

(39)

通過對各等級的復(fù)現(xiàn)倍數(shù)μi進(jìn)行控制,就可以合理地規(guī)劃種群等級結(jié)構(gòu),將各等級規(guī)?？刂圃诤侠淼姆秶鷥?nèi)。如圖6所示,本文設(shè)置復(fù)現(xiàn)倍數(shù)均等、線性遞減和負(fù)指數(shù)遞減3種帶種群結(jié)構(gòu)控制的選擇模型,以進(jìn)行種群結(jié)構(gòu)的控制,并觀察算法的表現(xiàn):

3.2 模型表達(dá)式及其參數(shù)

① 復(fù)現(xiàn)倍數(shù)均等(NSGA-Ⅱ-PSC-Ⅰ)

在復(fù)現(xiàn)倍數(shù)均等的種群結(jié)構(gòu)中,所有非支配等級的復(fù)現(xiàn)倍數(shù)相同,所以μi表達(dá)式中僅有一個常數(shù)a作為參數(shù)，因此,均等復(fù)現(xiàn)倍數(shù)模型為

μi=a,?i∈I。

(40)

該常數(shù)的表達(dá)式為

(41)

② 復(fù)現(xiàn)倍數(shù)線性遞減(NSGA-Ⅱ-PSC-Ⅱ)

設(shè)復(fù)現(xiàn)倍數(shù)線性遞減的模型為

μi=ai+b,?i∈I,a<0。

(42)

(43)

③ 復(fù)現(xiàn)倍數(shù)負(fù)指數(shù)遞減(NSGA-Ⅱ-PSC-Ⅲ)

設(shè)復(fù)現(xiàn)倍數(shù)負(fù)指數(shù)遞減的模型為

μi=ba-i,a>1,b>0。

(44)

I中元素最小值為1,設(shè)其最大值Imax=j,前四分之一分位數(shù)l=(1+j)/4,中值k=(1+j)/2。

(45)

(46)

3種模型的計算流程如圖7所示。

圖7 3種模型的計算流程Fig.7 Flow chart of three models

3.3 PSC選擇流程

在實際計算中,根據(jù)2.2中模型求得的μi可能仍不滿足約束(39)，而且s和ui都應(yīng)該是整數(shù),但μi往往不是整數(shù)。為解決這2個問題,采用以下步驟進(jìn)行進(jìn)一步的修整。

圖8 NSGA-Ⅱ-PSC流程圖 Fig.8 Flow chart of NSGA-Ⅱ-PSC

步驟1：i=1,根據(jù)上述流程計算μi的值。

步驟2：對μi向下取整得到[μi],使reffi復(fù)制[μi]次得到集合rseqi；

步驟3：對di(μi-[μi])向上取整得到[di(μi-[μi])],取qeffi中擁擠度最大的[di(μi-[μi])]個體,得到補(bǔ)充集合rsupi；

步驟4：步驟2和步驟3的集合相加得到新的等級i個體集合rnewi=rseqi+rsupi；

步驟6：新種群的規(guī)模為s′=|pnew|,此時有s′≥s。根據(jù)非支配等級和擁擠度,在pnew中選出最優(yōu)的s個個體組成最終種群psel。

NSGA-Ⅱ-PSC的總體流程及其中的PSC選擇流程如圖8所示。

3.4 模型效果對比

設(shè)計隨機(jī)路網(wǎng)并給定各節(jié)點的供應(yīng)鏈相關(guān)參數(shù),分別應(yīng)用帶傳統(tǒng)快速非支配選擇算子和3種種群結(jié)構(gòu)控制選擇算子的NSGA-Ⅱ進(jìn)行模型求解。種群規(guī)模設(shè)定為100,迭代次數(shù)設(shè)定為150次。

路網(wǎng)規(guī)模分3種,具體如下:

3*:3類危險品-每種危險品3個候選配送中心共30個客戶；

4*:4類危險品-每種危險品4個候選配送中心共40個客戶；

5*:5類危險品-每種危險品5個候選配送中心共50個客戶。

4種算法在3* 規(guī)模下,種群規(guī)模200和150次迭代的某次求解的迭代過程如圖9所示。

由于每個路網(wǎng)參數(shù)不一,因此為了方便比較,將所有結(jié)果進(jìn)行歸一化處理，將經(jīng)典NSGA-Ⅱ求得的Pareto最優(yōu)解集中各指標(biāo)值調(diào)整為設(shè)定值,其他值調(diào)整為其實際值與NSGA-Ⅱ?qū)嶋H值之比乘以NSGA-Ⅱ設(shè)定值。各設(shè)定值見表1。

表1 歸一化處理中的設(shè)定值Tab.1 Set values in normalization processing

將種群規(guī)模設(shè)定為100,迭代次數(shù)設(shè)定為150次,并分別記錄其第50、100、150次迭代時的目標(biāo)函數(shù)最小值。每種路網(wǎng)規(guī)模隨機(jī)生成10個路網(wǎng),每個路網(wǎng)重復(fù)求解10次,詳細(xì)的求解綜合對比見表2。

表2 不同迭代次數(shù)下4種算法性能比較Tab.2 Comparison of four algorithms under different iterations

可以看出,帶PSC算子的NSGA-Ⅱ比傳統(tǒng)NSGA-Ⅱ初期收斂較慢,在50次迭代時,PSC算子的最小成本和風(fēng)險平均比傳統(tǒng)NSGA-Ⅱ高10.56%和3.36%,等級1有效規(guī)模低38.25%。而從相對傳統(tǒng)NSGA-Ⅱ的平均增幅看,在PSC算法中,收斂最慢的是PSC-Ⅰ,最快的是PSC-Ⅲ。150次迭代時,PSC算子的表現(xiàn)則明顯反超傳統(tǒng)NSGA-Ⅱ選擇算子,成本和風(fēng)險平均比傳統(tǒng)NSGA-Ⅱ低3.45%和1.26%,Pareto解集規(guī)模則高出27.24%。

其中,在目標(biāo)函數(shù)最小化上表現(xiàn)最優(yōu)的是PSC-Ⅲ算子,2個目標(biāo)函數(shù)平均比傳統(tǒng)NSGA-Ⅱ低4.34%和1.92%。而在Pareto解集規(guī)模上表現(xiàn)最好的是PSC-Ⅰ,平均比傳統(tǒng)NSGA-Ⅱ高了35.59%,PSC-Ⅱ和PSC-Ⅲ則相差不多。

另外可以看出,隨著算例規(guī)模的增大,也就是問題的復(fù)雜化,PSC算子在Pareto解集規(guī)模上的優(yōu)勢越明顯,相比傳統(tǒng)NSGA-Ⅱ能得到更多的Pareto最優(yōu)解。

這是由于,在種群結(jié)構(gòu)上,高等級解占比由大到小分別是傳統(tǒng)NSGA-Ⅱ、PSC-Ⅲ、PSC-Ⅱ、PSC-Ⅰ。高等級占比大的種群結(jié)構(gòu)能夠加速算法的收斂,而等級占比小的種群結(jié)構(gòu)能讓算法探索更大的解空間。在問題太復(fù)雜時,收斂速度快的算法更有效率,而隨著問題變復(fù)雜,對解空間探索更為充分的算法往往能獲得更好的解集，因此,對種群結(jié)構(gòu)施加控制是提高遺傳算法質(zhì)量的一種思路。

最后,在運行時間上,PSC算法平均只多出4.05%,在可接受的范圍內(nèi)。

根據(jù)上述分析,平衡收斂速度和解集質(zhì)量,采用NSGA-Ⅱ-PSC-Ⅲ算法重新求解2.2節(jié)圖4的3*算例,圖10為求得的Pareto最優(yōu)解集。Pareto解集規(guī)模為44,比NSGA-Ⅱ提高46.7%,最低成本解的成本和風(fēng)險值分別為1 975.08和4 015.19,比NSGA-Ⅱ分別降低了4.64%和1.20%；最低風(fēng)險值的成本和風(fēng)險值分別為3 373.04和2 434.34,比NSGA-Ⅱ分別降低了2.31%和0.89%。圖11為Pareto解集中最低成本解、最低風(fēng)險解分別的運輸路線。

最低成本解共有3臺配送車輛,平均行駛里程為508.39 km,總里程為1525.16 km；最低風(fēng)險解共有5臺配送車輛,平均行駛里程為363.41 km,總里程為1 817.08 km；運輸車輛的增加不僅減少了單臺車輛的裝載量,降低了駕駛員的風(fēng)險,而且使得風(fēng)險得以分散化分布,也提高了整個區(qū)域的風(fēng)險公平性。

圖10 Pareto最優(yōu)解集Fig.10 Pareto optimal set

4 結(jié)語

本文結(jié)合需求不確定性與庫存理論,構(gòu)建了危險品定位庫存路徑模型,并提出了基于風(fēng)險分布特性和人口分布特性的風(fēng)險模型。此外,本文設(shè)計了能一次性求解LRP問題的算法,并在傳統(tǒng)NSGA-Ⅱ的基礎(chǔ)上提出了NSGA-Ⅱ-PSC及其3種類型。算例表明,對種群結(jié)構(gòu)進(jìn)行控制能夠影響NSGA-Ⅱ的性能,NSGA-Ⅱ-PSC在解集質(zhì)量上明顯優(yōu)于傳統(tǒng)NSGA-Ⅱ,且在復(fù)雜問題上優(yōu)勢更明顯，其中,PSC-Ⅰ能獲得最好的解集質(zhì)量,而PSC-Ⅲ的收斂速度最快。