?；肥鹿蔬B續(xù)消耗型應(yīng)急物資調(diào)度模型

盧建鋒，牟瑞芳，趙佳虹，王小霞

(1.西南交通大學 交通運輸與物流學院，四川 成都 610031；2.廣東工業(yè)大學 土木與交通工程學院，廣東 廣州 510006)

危化品在其生產(chǎn)、儲存以及使用過程中極容易發(fā)生事故。在過去的幾十年里，危險品事故頻發(fā)，給自然環(huán)境和人類生活帶來嚴重危害。在?；肥鹿实膽?yīng)急救援中，除了需要一次性消耗的應(yīng)急物資，如衣服、帳篷外，還需要連續(xù)消耗的應(yīng)急物資，即應(yīng)急活動中需持續(xù)不斷供應(yīng)的資源類型，比如燃氣、食品、醫(yī)藥等。一直以來，連續(xù)消耗型物資調(diào)度問題得到許多學者的關(guān)注[1-2]。劉春林等[3]分析了應(yīng)急物資調(diào)度系統(tǒng)多點出救的特點，引入了連續(xù)可行方案的概念，并提出以應(yīng)急開始時間最早為目標的數(shù)學模型。潘郁等[4]研究了連續(xù)性消耗應(yīng)急物資的調(diào)度，從應(yīng)急活動成本和受災(zāi)損失2方面考慮，構(gòu)建了總成本最低的應(yīng)急物資調(diào)度模型。魏國強等[5]建立規(guī)劃模型，解決了多應(yīng)急點連續(xù)消耗應(yīng)急物資調(diào)度預(yù)案的制定問題。鄭昊等[6]假設(shè)應(yīng)急資源的消耗速率為非負可積函數(shù)，構(gòu)建了應(yīng)急時間最早為目標的連續(xù)型應(yīng)急物資調(diào)度模型。張力丹等[7]將災(zāi)民安置與物資運輸一體化規(guī)劃，以施救不及時損失和應(yīng)急系統(tǒng)施救成本為優(yōu)化目標，構(gòu)建連續(xù)消耗應(yīng)急物資調(diào)配模型，并用遺傳算法和序列線性規(guī)劃算法相結(jié)合來求解模型。在應(yīng)急物資調(diào)度優(yōu)化目標方面，目前國內(nèi)外研究較多的包括: 出救點最少、應(yīng)急時間開始最早、最短運輸時間和最小運輸成本等[8-11]。

以上關(guān)于連續(xù)消耗應(yīng)急物資調(diào)度的研究主要考慮調(diào)度效率、成本以及受災(zāi)損失等，并不完全適用于?；肥鹿蕬?yīng)急物資調(diào)度。因為對于?；肥鹿识裕l(fā)生事故會對周邊污染環(huán)境，有毒氣體將擴散到道路網(wǎng)中形成毒性負荷，在調(diào)度時應(yīng)該考慮毒路徑毒性負荷，本文在現(xiàn)有研究基礎(chǔ)上，針對危化品事故應(yīng)急物資調(diào)度特征，構(gòu)建以缺貨損失最小、調(diào)運時間最短和調(diào)運毒性負荷最小為優(yōu)化目標的連續(xù)消耗應(yīng)急物資調(diào)度模型，并采用NSGAⅡ遺傳算法求解模型，得出多種可供參考的應(yīng)急物資調(diào)度方案，采用逼近理想解的排序方法選擇推薦方案。

1 問題描述及模型建立

1.1 問題描述

設(shè)?；肥鹿手挥幸粋€受災(zāi)點S，共需要m種連續(xù)消耗型應(yīng)急物資，周邊有n個應(yīng)急中心點能向受災(zāi)點S供應(yīng)應(yīng)急物資。每個應(yīng)急中心到受災(zāi)點之間有多條推薦路徑，每條路徑時間和毒性負荷均已提前預(yù)測得到。要求給出一個應(yīng)急中心應(yīng)急資源調(diào)度方案，包括確定應(yīng)急中心點發(fā)出的應(yīng)急物資數(shù)量及其所選擇路徑。要求調(diào)度模型盡量實現(xiàn)以下目標：1) 受災(zāi)點缺貨損失最小，即使得受災(zāi)點因物資缺乏導(dǎo)致的損失最??；2) 完成應(yīng)急調(diào)度運輸?shù)乃臅r間最短；3) 調(diào)度毒性負荷最小。

1.1.1 假設(shè)條件

模型建立的假設(shè)條件為

1) 各種類型物資都可以混裝在一輛貨車；

2) 受災(zāi)點的多種物資的消耗速度是確定不變的；

3) 每個出救點運力充足；

4) 運輸過程中不同路徑的運輸時間和運輸風險是己知的。

1.1.2 符號說明

I為出救點集，i∈I；

J為資源種類集，j∈J；

K為出救點到受災(zāi)點的路徑標號集，k∈K；

Ai為 第i個出救點(i=1,2,···,n)；

S為受災(zāi)點；

hik為從出救點i出發(fā)的第k條路徑；

wik為出救點i到受災(zāi)區(qū)第k條路徑的環(huán)境風險值，用路徑毒性負荷來表示，計算方法見文獻[12]；

Tik為出救點i到受災(zāi)區(qū)第k條路徑的行程時間；

yik表示出救點是否選擇第k條路徑，yik=1時為選擇，yik=0時不選擇；

G0為標準貨車載重量；

Ni為出救點i派出的貨車數(shù)量；

vj為應(yīng)急活動開始后j類物資消耗速率；

ti為從出救點i到達受災(zāi)點花費的時間；

t0為救援的機會時間；

l為各出救點車輛到達受災(zāi)點的順序，l=1，2，…，n；

t(l)為第l個到達受災(zāi)點所用的時間；

T為應(yīng)急活動結(jié)束時間；

dmaxij為 出救點Ai處第j類物資最大可供應(yīng)量；

dj(l)為t(l)時刻到達受災(zāi)點的j類應(yīng)急物資數(shù)量；

lj(t)為t時刻受災(zāi)點第j類應(yīng)急物資存量；

Fj為受災(zāi)點因j類應(yīng)急物資缺失而導(dǎo)致的單位時間損失費用；

Ej0為 從機會時間t0到第1個出救點物資到達時，因j類應(yīng)急物資缺失導(dǎo)致的損失費用；

Ejl為第1個出救點物資到達以后，因j類應(yīng)急物資缺失而導(dǎo)致的損失費用；

Ej為受災(zāi)點因j類應(yīng)急物資缺失而導(dǎo)致的總損失費用；

xij為 出救點Ai提供的j類物資量。

1.2 應(yīng)急物資缺貨損失計算方法

在受災(zāi)點應(yīng)急物資連續(xù)消耗的情況下，受災(zāi)點的物資庫存量如圖1所示。

引入一個判斷變量Ril，設(shè)定Ril=1時為應(yīng)急中心i在t(l)時 刻為受災(zāi)點供貨，Ril=0時為其 他情況，則即n個應(yīng)急中心有且只有一個是第l個到達受災(zāi)點。

t(l)第l個到達受災(zāi)點的時間只能是n個應(yīng)急中心中一個。

圖1 第j類應(yīng)急物資連續(xù)消耗下庫存量變化圖Figure 1 Change chart of inventory under continuous consumption of category j emergency material

x(l)第l個到達受災(zāi)點的物資只能是n個應(yīng)急中心中的一個。

從事故發(fā)生時刻t=0到機會時間t0的損失不可避免。應(yīng)急物資缺失導(dǎo)致受災(zāi)點的損失費用計算方法如下[13]。

事故發(fā)生時，從機會時間t0時到第1個應(yīng)急中心物資到達時刻t(1)，這段時間內(nèi)受災(zāi)點物資庫存量必然缺失，因此得出第j類物資缺失的損失費用為

為不失一般性，設(shè)t(k)≤T，且t(n+1)=T，則 在t(k)時刻受災(zāi)點j類物資庫存量為

如圖1所示，存在t(F)∈[t(k),t(k+1)]，且t(k))=0時，受災(zāi)點第j類物資有部分存余，但卻無法滿足整個周期 [t(k),t(k+1)]的需求，物資存在部分缺失。此時，有且有0 ≤Ij(t(k))≤vj[t(k+1)?t(k)]，得第j類物資缺失的損失費用為

其他情況下，Bjk=0。

因此，應(yīng)急物資缺貨損失費用按以下2種情況處理。

1)Bjk=0，第i時刻受災(zāi)點各類應(yīng)急物資都充足，即缺失量為0。則，因第j類應(yīng)急物資缺失導(dǎo)致的損失費用只是從機會時間t0時到第一個應(yīng)急中心物資到達時刻t(1)的時間段內(nèi)的物資缺失。第j類物資缺貨損失費用為

2) 當?shù)趈類物資應(yīng)急物資在周期 [t(k),t(k+1)]內(nèi)缺失，即存在t(F)∈[t(k),t(k+1)]，使得

第j類物資缺貨損失費用為

1.3 模型建立

要求給出應(yīng)急物資調(diào)度方案{dij}，具體包括參與應(yīng)急的應(yīng)急中心選擇及其提供的物資數(shù)量。模型建立如下。

Ni表示某事故階段某應(yīng)急中心要派的貨車數(shù)量，其值為該點派出的總貨物質(zhì)量除以單位貨車載重量向上取整數(shù)。

式(7)表示受災(zāi)點缺貨損失最少；式(8)表示使物資調(diào)運過程中環(huán)境風險最??；式(9)表示應(yīng)急物資運送到受災(zāi)區(qū)的總調(diào)運時間最短；式(10)表示應(yīng)急中心提供的物資數(shù)量介于d0和dmaxij之間，或能提供的物資數(shù)量為0；式(11)表示各類物資總供應(yīng)量大于需求量；式(12)表示t(k)不 小于機會時間t0，也不大于應(yīng)急結(jié)束時間T；式(13)表示ti大于機會時間t0，小于應(yīng)急結(jié)束時間T；式(14)表示第k個到達受災(zāi)點的時間只能是n個應(yīng)急中心中的一個；式(15)表示第k個到達受災(zāi)點的物資只能是n個應(yīng)急中心中的一個；式(16)表示每個應(yīng)急中心到達受災(zāi)點選擇的路徑有且只有一條。

2 算法設(shè)計

由式(7)、式(8)、式(9)及其約束條件可知，模型屬于非線性多目標優(yōu)化問題，存在多個Pareto解，且需要結(jié)合決策者的偏好對解進行選擇?？紤]求解效率，可采用非劣分層遺傳算法(NSGAⅡ)求出模型Pareto最優(yōu)解集，構(gòu)成決策矩陣。在選擇最優(yōu)解時，用信息熵法[14]確定3個目標屬性的權(quán)重，最后采用逼近理想解的排序方法(TOPSIS)對Pareto最優(yōu)解排序，依據(jù)決策者偏好選出最優(yōu)解。

2.1 非劣分層遺傳算法(NSGAⅡ)

非劣分層遺傳算法(NSGAⅡ)是目前求解多目標問題應(yīng)用較好的一種算法。NSGAⅡ在NSGA算法基礎(chǔ)上加入精英策略、快速非支配排序策略以及密度值估計策略。借鑒NSGA-Ⅱ算法中的選擇機制，對交叉和變異算子進行適當修改。算法步驟如下[15]。

步驟1初始化種群規(guī)模 Popsize，對最大迭代次數(shù)、交叉概率和變異概率進行設(shè)置，計數(shù)器gen=1。

步驟2生成滿足約束的個體Popsize個，計算每個個體的目標函數(shù)值，對所有個體進行非支配排序，生成父代種群。

步驟3利用錦標賽選擇、多點交叉和實值變異方法，對父代進行遺傳操作，生成臨時種群。

步驟4判斷臨時種群的個體滿足約束條件否，找出不滿足這些約束的個體，產(chǎn)生滿足這些約束的相同數(shù)量的個體替代，得出子代種群，計算個體后代的目標函數(shù)值。

步驟5將父群體和子群體合并，并進行非支配排序，將Popsize值最高的個體作為父代種群，gen=gen+1。

步驟6如果 gen≤Maxgen，轉(zhuǎn)入第3步，否則結(jié)束。

以上算法的關(guān)鍵部分操作說明如下。

1) 染色體的編碼。用1個染色體來表示從應(yīng)急中心Ai至受災(zāi)點S的物資調(diào)度方案。染色體的每個基因?qū)?yīng)每個應(yīng)急中心各類物資配送到受災(zāi)點的數(shù)量及該應(yīng)急中心所選路徑。

2) 適應(yīng)度函數(shù)。適應(yīng)度由目標函數(shù)值決定，每條染色體都具有3個適應(yīng)度值f1、f2和f3。

3) 遺傳操作。選擇: 根據(jù)快速非支配的排序以及擁擠距離大小，采用錦標賽法做選擇。交叉和變異: 采用遺傳算法的多點交叉算子，加快了求解空間的搜索速度，避免了早熟收斂。利用實值變異和區(qū)域掃描限制變異范圍，使新個體無法跳出解空間中相應(yīng)解。

2.2 逼近理想解的排序方法(TOPSIS)

TOPSIS方法通過比較和理想解的相近程度，計算與理想解和負理想解的距離來對方案排序，計算流程如下[16]。

1) 用向量規(guī)范法得出規(guī)范決策矩陣。假設(shè)有模型有n個多目標屬性，m個決策方案。原決策矩陣Y={yij}, 規(guī)范化后的決策矩陣為Z={zij}, 則

2) 構(gòu)成加權(quán)規(guī)范矩陣。

設(shè)多屬性的權(quán)向量w={w1,w2,···,wn}T，則

3) 求理想解x+和負理想解x?。設(shè)理想解x+的第j個屬性值為設(shè)負理想解x?的第j個屬性值為則正理想解為為效益型屬性)，或(j為成本型屬性)。負理想解為為成本型屬性)，或(j為效益型屬性)。

4) 計算各備選方案到理想點和負理想點距離。備選方案到理想點距離為

備選方案到負理想點距離為

5)計算各方案排隊指示值(綜合評價指數(shù))。

6)根據(jù)Ci值由大到小排列方案的優(yōu)劣次序。

3 數(shù)值與算例分析

3.1 算例數(shù)據(jù)

設(shè)某?；穫}庫發(fā)生泄露并引發(fā)爆炸等事故，倉庫周邊一定范圍成為受災(zāi)點?，F(xiàn)急需4種應(yīng)急物資，可從5個應(yīng)急中心點安排物資調(diào)度。受災(zāi)點和各應(yīng)急中心點的分布如圖2所示，各應(yīng)急中心到受災(zāi)點的路徑的環(huán)境風險值和時間如表1所示，各應(yīng)急中心物資存有量如表2所示。各類應(yīng)急物資的需求總量為X=(40，30，40，45) t，各類應(yīng)急物資的消耗速率V=(3，5，4，1) t/min，各類物資缺失后的單位時間損失費為F=(7，2，9，5) 元/min，應(yīng)急救援活動的終止時間T為80 min，標準貨車載重設(shè)定為18 t。假設(shè)到達受災(zāi)點的不同路徑時間和毒性負荷提前計算得到，要求在滿足受災(zāi)點應(yīng)急物資需求的條件下，求出最優(yōu)調(diào)度方案。

3.2 求解結(jié)果

圖2 ?；肥鹿适転?zāi)點和救援點位置示意圖Figure 2 Schematic map of disaster site and rescue site of dangerous chemicals accident

表1 各應(yīng)急中心點到受災(zāi)點的路徑相關(guān)數(shù)據(jù)Table 1 Path-related data from emergency centers to the disaster site

表2 各應(yīng)急中心物資數(shù)量Table 2 Quantity of material at material emergency center t

在Intel Core Intel Core i5-337U＠1.8 GHz處 理器，8 GB 內(nèi)存的計算機上用Matlab語言編寫NSGAⅡ算法程序?qū)λ憷抡?。多次測試后，參數(shù)確定為：初始種群數(shù)量為500，交叉概率為0.8，變異概率為0.04，最大進化代數(shù)為1000。模型目標函數(shù)共求得52個Pareto解，如圖3所示。

圖3 NSGAⅡ算法帕累托最優(yōu)解示意圖Figure 3 Schematic diagram of Pareto optimal solution for NSGA II algorithm

用這些解構(gòu)造52×3的決策矩陣。先利用信息熵法計算得缺貨損失、環(huán)境風險、調(diào)運時間屬性的權(quán)值向量為w=[0.04，0.86，0.01]。

再運用TOPSIS法，對52個Pareto解進行排序，排序第1的解為滿意解。該解作為應(yīng)急物資調(diào)度推薦方案，如表3所示。方案的缺貨損失、環(huán)境風險、調(diào)運時間屬性值向量為P=[107150，850，40]。

表3 應(yīng)急物資調(diào)度推薦方案Table 3 Recommended schemes for emergency material dispatching

此外，缺貨損失最小方案如表4所示。該方案的缺貨損失、環(huán)境風險、調(diào)運時間屬性值向量為P=[104510，2140，40]。

表4 缺貨損失最小方案Table 4 Shortage loss minimum scheme

將推薦方案和受災(zāi)損失最小方案對比，結(jié)果如表5所示?？砂l(fā)現(xiàn)推薦方案的環(huán)境風險值比受災(zāi)損失最小方案節(jié)約近一倍，而受災(zāi)損失相差很小，行程時間相同?？梢娡扑]方案有較大的優(yōu)勢，也說明了算法的合理性。

表5 兩方案目標值對比表Table 5 Target value contrast table of two schemes

4 結(jié)論

?；肥鹿食３０殡S毒氣泄露，污染周邊道路網(wǎng)。進行應(yīng)急物資調(diào)度時考慮路徑毒性負荷，能減少調(diào)度工作人員二次傷害。本文針對?；肥鹿蕬?yīng)急物資調(diào)度的特點，建立了危化品事故連續(xù)性消耗的應(yīng)急物資調(diào)度優(yōu)化模型，以缺貨損失最小、調(diào)運時間最短和調(diào)運總環(huán)境風險最小為優(yōu)化目標，采用NSGAⅡ遺傳算法求解，得出多種可供參考的應(yīng)急物資調(diào)度方案，采用逼近理想解的排序方法選擇推薦方案。通過算例分析，將推薦方案和受災(zāi)損失最小方案對比，發(fā)現(xiàn)推薦方案的毒性負荷比受災(zāi)損失最小方案小近一倍，而受災(zāi)損失相差很小，調(diào)運時間相同。說明推薦方案有較大的優(yōu)勢，也說明算法合理可行。