張聆曄, 呂 靖, 班 豪, 范瀚文

(大連海事大學(xué) 交通運(yùn)輸工程學(xué)院,遼寧 大連 116026)

0 引言

隨著海上油氣勘探的深入和原油運(yùn)輸量的遞增，海上發(fā)生重大溢油事故的風(fēng)險(xiǎn)逐年上升[1]。一旦出現(xiàn)海上重大溢油必將導(dǎo)致嚴(yán)重的財(cái)產(chǎn)損失與難以估計(jì)的海洋生態(tài)環(huán)境破壞。此時(shí)，科學(xué)高效的物資調(diào)度對(duì)溢油應(yīng)急工作的開展具有十分重要的意義。

目前，學(xué)者們已對(duì)陸上應(yīng)急物資的調(diào)度進(jìn)行了廣泛研究。為解決應(yīng)急響應(yīng)過程中存在的物資需求動(dòng)態(tài)變化以及運(yùn)輸網(wǎng)絡(luò)的不確定性問題，現(xiàn)有研究常用的方法主要有三種，分別為周期規(guī)劃法[2](Time Period Planning Horizon)、滾動(dòng)規(guī)劃法[3](Rolling Planning Horizon)以及隨機(jī)規(guī)劃法[4](Stochastic Programming)。

現(xiàn)有研究重在探討陸上災(zāi)害或事故的應(yīng)急響應(yīng)，針對(duì)海上應(yīng)急物資調(diào)度的研究則較為有限[5]。事實(shí)上，由于海陸環(huán)境存在較大差異，適用于陸上災(zāi)害應(yīng)急響應(yīng)的常規(guī)方法不足以應(yīng)對(duì)海上溢油應(yīng)急響應(yīng)。這些方法的應(yīng)用是基于假設(shè)受災(zāi)區(qū)域的物資需求以及運(yùn)輸網(wǎng)絡(luò)的狀態(tài)在不同的周期或情景內(nèi)能維持相對(duì)穩(wěn)定，直至下一個(gè)周期開始或轉(zhuǎn)換到另一個(gè)情景中才對(duì)應(yīng)急系統(tǒng)內(nèi)關(guān)于需求量與運(yùn)輸路徑相關(guān)的參數(shù)進(jìn)行更新以反映應(yīng)急響應(yīng)過程中存在的動(dòng)態(tài)特性。在溢油應(yīng)急響應(yīng)過程中，溢油油膜具有擴(kuò)散和漂移特征，這些動(dòng)態(tài)過程不存在停頓與周期性變更。這樣的特點(diǎn)決定了各需求點(diǎn)對(duì)應(yīng)急物資的需求量以及各需求點(diǎn)之間的距離具有時(shí)變特征。溢油應(yīng)急物資調(diào)度優(yōu)化問題勢(shì)必需要通過構(gòu)建具有時(shí)變特征的動(dòng)態(tài)優(yōu)化模型進(jìn)行探討。

為此，本文在綜合考慮受污染區(qū)域的時(shí)變物資需求、相關(guān)運(yùn)輸網(wǎng)絡(luò)的不確定性以及應(yīng)急物資調(diào)度決策與外部決策環(huán)境之間的相互作用關(guān)系后，以生態(tài)環(huán)境損失值最小和應(yīng)急響應(yīng)總成本最小為目標(biāo)構(gòu)建了一個(gè)動(dòng)態(tài)多目標(biāo)選址-路徑規(guī)劃模型，并基于鯨魚算法提出了相應(yīng)的求解方法。

1 問題描述以及假設(shè)提出

本文主要以圍油欄調(diào)度為例探討海上溢油應(yīng)急物資的調(diào)度問題。由于圍油欄的需求量與溢油油膜周長(zhǎng)具有正相關(guān)關(guān)系，溢油油膜擴(kuò)散時(shí)對(duì)圍油欄的需求量也將相應(yīng)增加。此外，油膜的漂移運(yùn)動(dòng)也將導(dǎo)致應(yīng)急船的運(yùn)輸網(wǎng)絡(luò)也具有動(dòng)態(tài)特征。更為重要的是，在進(jìn)行海上溢油應(yīng)急物資調(diào)度時(shí)除了需要考慮油膜動(dòng)態(tài)變化對(duì)應(yīng)急物資調(diào)度決策產(chǎn)生的影響外，也需要考慮物資調(diào)度本身對(duì)油膜變化產(chǎn)生的反作用(如迅速地利用圍油欄對(duì)溢油進(jìn)行圍控能及時(shí)抑制油膜擴(kuò)散，從而減少需求點(diǎn)對(duì)圍油欄的需求量)。此時(shí)外部決策環(huán)境與應(yīng)急物資調(diào)度決策間存在相互作用。

模型構(gòu)建基于以下假設(shè)：

(1)各受污染點(diǎn)的分布情況以及溢油數(shù)量可通過地理信息系統(tǒng)(GIS)獲得；

(2)油膜擴(kuò)散速率和漂移速度能維持穩(wěn)定，應(yīng)急船也可保持勻速行駛；

(3)在發(fā)生物資或運(yùn)力短缺的情況時(shí)可及時(shí)獲得物資供應(yīng)方和鄰近港口的支援(無等待時(shí)長(zhǎng))；

(4)同一個(gè)船隊(duì)中的所有應(yīng)急船同屬于一個(gè)應(yīng)急基地，不存在異地應(yīng)急船混編現(xiàn)象。

2 模型建立

2.1 相關(guān)符號(hào)說明

集合:

A：需求點(diǎn)集合，其中a為需求點(diǎn)點(diǎn)引索，a∈A={1,2,…,n}；

B：應(yīng)急基地集合，其中b為應(yīng)急基地引索，b∈B={n+1,n+2,…,n+m}；

V：應(yīng)急系統(tǒng)中所有節(jié)點(diǎn)集合，V=A∪B；

R：應(yīng)急船隊(duì)集合，r為船隊(duì)引索，r∈R。

決策變量:

urb：應(yīng)急船隊(duì)r是否屬于應(yīng)急基地b，若是urb=1，否則，urb=0，r∈R,b∈B；

yrij：應(yīng)急船隊(duì)r是否從節(jié)點(diǎn)i出發(fā)駛向節(jié)點(diǎn)j，若是yrij=1，否則，yrij=0，r∈R,i,j∈V，；

zra：需求點(diǎn)a是否由應(yīng)急船隊(duì)r完成物資的運(yùn)輸與投放，若是zra=1，否則，zra=0，r∈R,a∈A，；

kr：應(yīng)急船隊(duì)r中所包含的應(yīng)急船數(shù)量，r∈R；

qr：應(yīng)急船隊(duì)r所裝載的物資總量，r∈R。

2.2 溢油應(yīng)急物資調(diào)度優(yōu)化模型

考慮到成功的突發(fā)事件應(yīng)急響應(yīng)應(yīng)該兼顧效率與成本[6]，溢油應(yīng)急物資調(diào)度的優(yōu)化目標(biāo)也據(jù)此設(shè)定。具體模型構(gòu)建如下：

目標(biāo)函數(shù)

(1)

為對(duì)生態(tài)環(huán)境損失進(jìn)行估值，本文根據(jù)佛羅里達(dá)公式[7]構(gòu)建了如下等式：

(2)

此外，本文主要采用Liu和Leendertes[8]提出的一種溢油擴(kuò)散模型來反映油膜的變化過程，具體模型如下:

(3)

其中，f(t,Voa)為需求點(diǎn)a處體積為Voa的溢油在t時(shí)刻的油膜直徑；ρw為海水密度；ρ0為溢油密度，β=1-ρ0/ρw；g為重力加速度；Vo為此處的溢油體積；vw為粘滯系數(shù)；t為時(shí)間；δ為凈表面張力系數(shù)。

minZ2=OVOC+SVOC+ORCC+SRCC+TNC+TDC

(4)

式(4)為成本目標(biāo)函數(shù)，表示最小化應(yīng)急響應(yīng)成本，整個(gè)應(yīng)急響應(yīng)過程中產(chǎn)生的成本主要包括如下幾個(gè)部分：應(yīng)急基地自有應(yīng)急船的固定運(yùn)營成本(OVOC)，臨時(shí)征用的應(yīng)急船的固定運(yùn)營成本(SVOC)，基地自有應(yīng)急物資的消耗成本(ORCC)，臨時(shí)增購的應(yīng)急物資的消耗成本(SRCC)，應(yīng)急船隊(duì)在節(jié)點(diǎn)間的運(yùn)輸成本(TNC)，應(yīng)急船隊(duì)隨需求點(diǎn)漂移的運(yùn)輸成本(TDC)。式中，kb應(yīng)急基地b所配備的應(yīng)急船總數(shù)，b∈B；qb表示應(yīng)急基地b所配備的應(yīng)急物資數(shù)量，b∈B；vk表示應(yīng)急船的航速；trij表示應(yīng)急船隊(duì)r在節(jié)點(diǎn)i和j之間的運(yùn)輸時(shí)長(zhǎng)，i,j∈V；ct表示應(yīng)急船的單位距離運(yùn)輸成本；trj表示應(yīng)急船隊(duì)r到達(dá)節(jié)點(diǎn)j的時(shí)刻，j∈V;Tri表示應(yīng)急船隊(duì)r正要離開節(jié)點(diǎn)i的時(shí)刻，i∈V；ooc表示應(yīng)急基地自有應(yīng)急船的固定運(yùn)營成本；soc表示向港口臨時(shí)征調(diào)補(bǔ)充的應(yīng)急；oqc表示應(yīng)急基地自有物資單位使用成本;sqc表示臨時(shí)補(bǔ)給物資單位使用成本。

約束條件

(20)

式(5)表示應(yīng)急船隊(duì)r到達(dá)節(jié)點(diǎn)j的時(shí)刻；式(6)表示運(yùn)輸船隊(duì)r在節(jié)點(diǎn)間的運(yùn)輸時(shí)間，F(xiàn)(Tri)的具體展開內(nèi)容見式(21)；式(7)、(8)共同表示應(yīng)急船隊(duì)r離開不同種類節(jié)點(diǎn)的時(shí)刻，式中er表示應(yīng)急船隊(duì)r的物資投放效率，r∈R；式(9)為基于式(3)所構(gòu)建的需求量時(shí)變函數(shù)，當(dāng)應(yīng)急物資的種類確定為圍油欄時(shí)，需求點(diǎn)對(duì)圍油欄的需求數(shù)量等于油膜周長(zhǎng)，式中da(t)表示在t時(shí)刻需求點(diǎn)a對(duì)應(yīng)急物資的需求量，而gi(t,Voi)表示在t時(shí)刻需求點(diǎn)a的需求量增長(zhǎng)率；式(10)表示應(yīng)急船隊(duì)r的物資投放效率(圍油欄布置效率)，式中ek表示應(yīng)急船的物資投放效率；式(11)表示每個(gè)應(yīng)急船隊(duì)所配備的運(yùn)力恰好滿足由其負(fù)責(zé)的所有需求點(diǎn)的應(yīng)急物資需求總量，式中qk表示應(yīng)急船的最大載貨量；式(12)表示由應(yīng)急基地i組建的應(yīng)急船隊(duì)只能從該應(yīng)急基地出發(fā)；式(13)表示應(yīng)急船隊(duì)不會(huì)在應(yīng)急基地間往返；式(14)表示每個(gè)節(jié)點(diǎn)有應(yīng)急船隊(duì)進(jìn)，就一定有應(yīng)急船隊(duì)出，保持流量平衡；式(15)表示每個(gè)節(jié)點(diǎn)只通過一次應(yīng)急船隊(duì)，避免形成未經(jīng)過應(yīng)急基地的循環(huán)路線；式(16)表示每一個(gè)需求點(diǎn)有且僅有一個(gè)應(yīng)急船隊(duì)對(duì)其負(fù)責(zé)；式(17)表示應(yīng)急船隊(duì)途經(jīng)路徑上的節(jié)點(diǎn)時(shí)存在先后順序；式(18)～(20)是決策變量的0-1約束。

此外，由于油膜具有漂移特征，當(dāng)應(yīng)急船隊(duì)r完成當(dāng)前需求點(diǎn)i的運(yùn)輸任務(wù)后向下一個(gè)需求點(diǎn)航行時(shí)，航行距離也因此會(huì)處于時(shí)變狀態(tài)。海上各節(jié)點(diǎn)間的位置關(guān)系可表示如下:

圖1中ai、aj為需求點(diǎn)，此時(shí)(t=Tri)應(yīng)急船隊(duì)r剛好完成對(duì)點(diǎn)ai的物資投放并準(zhǔn)備前往點(diǎn)aj，應(yīng)急船隊(duì)r的航速為vk，點(diǎn)aj的漂移速度為vj，船隊(duì)從點(diǎn)ai行進(jìn)到點(diǎn)aj的航行時(shí)長(zhǎng)為trij，γ為方位夾角，disij(Tri)表示點(diǎn)ai與點(diǎn)aj的當(dāng)前距離，其中γ、disij(Tri)可通過已知初始條件求得。利用余弦定理構(gòu)建三邊關(guān)系等式后，經(jīng)數(shù)學(xué)推導(dǎo)可將應(yīng)急船隊(duì)的航行時(shí)間表示如下:

trij=F(Tri)=[-vj·disij(Tri)·cosγ+disij(Tri)·

(21)

3 算法設(shè)計(jì)

本文基于鯨魚算法[9]提出了一種改進(jìn)的多目標(biāo)鯨魚算法(Improved Multi-objective Whale Optimization Algorithm, IMWOA)對(duì)構(gòu)建的優(yōu)化模型進(jìn)行求解。

IMWOA的具體計(jì)算步驟如下:

Step1初始化鯨魚種群。

Step2計(jì)算鯨魚個(gè)體所在位置的適應(yīng)度值。

Step3非支配解篩選。

Step4外部檔案的更新與維護(hù)。

Step5確定領(lǐng)頭鯨魚X*。

現(xiàn)有研究大多對(duì)領(lǐng)頭鯨魚進(jìn)行隨機(jī)選定[11]，這樣做容易導(dǎo)致計(jì)算結(jié)果局部收斂與早熟。在小生境共享機(jī)制中相似個(gè)體數(shù)越多，則此類個(gè)體間的共享度越大，其對(duì)應(yīng)的小生境適應(yīng)度就越小，反之亦然[12]。鑒于小生境共享機(jī)制這一特點(diǎn)，本文決定據(jù)此確定外部檔案中各個(gè)體的小生境適應(yīng)度值，并在此基礎(chǔ)上通過輪盤賭篩選作為領(lǐng)頭鯨魚X*的個(gè)體，從而確保解的多樣性，以避免局部收斂和早熟情況的出現(xiàn)。

Step6更新種群中鯨魚個(gè)體的位置。

A=2·a·r-a

(22)

(23)

式中，a的取值介于區(qū)間[0,2]且隨迭代次數(shù)τ的遞增而線性遞減；r為介于區(qū)間[0,1]的均勻分布隨機(jī)數(shù)。具體更新方式見文獻(xiàn)[9]。

在標(biāo)準(zhǔn)WOA中，全局搜索和局部搜索之間的取舍主要由收斂因子A決定。隨著a的線性遞減，A也隨之遞減，易導(dǎo)致領(lǐng)頭鯨魚陷入局部最優(yōu)。為此，本文采用一種先緩慢遞減后急劇遞增的函數(shù)(見式(24))代替標(biāo)準(zhǔn)WOA中關(guān)于收斂參數(shù)a的線性遞減函數(shù)(見式(23))。

(24)

Step7重復(fù)Step 2～4，并獲得更新后的外部存檔，隨后轉(zhuǎn)到Step 8，判斷是否終止迭代計(jì)算。

Step8判斷終止條件。若滿足最大迭代次數(shù)(τ=T)，則直接輸出計(jì)算結(jié)果；否則，τ←τ+1，轉(zhuǎn)到Step 5，繼續(xù)搜索。

4 仿真案例分析

4.1 背景及數(shù)據(jù)

從事故發(fā)生到對(duì)事故海域完成溢油監(jiān)測(cè)并采取實(shí)際的應(yīng)急行動(dòng)，時(shí)間間隔有7小時(shí)。此時(shí)，海面已出現(xiàn)數(shù)量眾多且分散的油膜區(qū)域，其中共有20處油膜較厚的需求點(diǎn)急需圍油欄對(duì)其進(jìn)行圍控。這些油膜在兩股不同海流作用下主要出現(xiàn)兩種不同的漂移方向與速度，部分油膜以1.8kn速度向南偏東32°方向漂移，另一部分油膜以1.2kn速度向北偏西28°方向漂移，油膜數(shù)量、位置及其漂移速度等相關(guān)信息見圖2。事故海域近鄰的應(yīng)急基地共有3處，分別是大連基地、煙臺(tái)基地以及威海基地，各應(yīng)急基地的物資庫存、船舶配備、物資消耗成本等與應(yīng)急物資相關(guān)的信息見表1。應(yīng)急船的航速、運(yùn)力以及使用成本等信息見表2。

圖2 需求點(diǎn)分布情況

表1 應(yīng)急基地物資儲(chǔ)備情況

表2 應(yīng)急船性能及其使用成本

4.2 算法比較

為求解模型，本文利用MATLAB軟件對(duì)IMWOA算法進(jìn)行程序編寫，并將IMWOA算法相關(guān)參數(shù)經(jīng)多次試驗(yàn)后設(shè)置如下：種群規(guī)模N=50，最大迭代次數(shù)T=300，共享距離σshare=2.8。由于NSGA-II被廣泛應(yīng)用于多目標(biāo)問題求解且該算法的有效性也獲得了研究領(lǐng)域的普遍認(rèn)可，故選擇其作為對(duì)比算法以檢驗(yàn)利用IMWOA所求得的Pareto解集的質(zhì)量。為保證分析結(jié)果對(duì)比的有效性，兩種算法的參數(shù)盡量一致，關(guān)于NSGA-II的相關(guān)參數(shù)設(shè)置如下：種群規(guī)模為50，最大迭代次數(shù)為300，交叉概率為0.8，變異概率為0.05。為求得精確Pareto前沿的近似解，兩種算法各運(yùn)行10次并取整合后的結(jié)果作為計(jì)算結(jié)果，如圖3所示。

圖3 IMWOA和NSGA-II的Pareto解分布

經(jīng)對(duì)比可發(fā)現(xiàn)，NSGA-II求得的Pareto解均受IMWOA求得的解的支配，IMWOA在精確度方面明顯優(yōu)于NSGA-II。就多樣性而言，IMWOA求得的解相較于NSGA-II而言數(shù)量更多，分布更為均勻。而NSGA-II求得的解則明顯分布不均，魯棒性較弱，主要是由于算法運(yùn)行到中后期時(shí)難以逃脫局部最優(yōu)而陷于局部區(qū)域搜索。由此可見，IMWOA具有更高的求解精度與更強(qiáng)的搜索能力，不僅能有效避免陷入局部最優(yōu)，同時(shí)還保證了解的多樣性，使Pareto解集分布更加廣泛，能夠?yàn)闆Q策者提供更優(yōu)質(zhì)、更充足的應(yīng)急物資調(diào)度決策參考。

4.3 計(jì)算結(jié)果

利用IMWOA求出Pareto解集后根據(jù)理想解選取模型[4]求得評(píng)分最高的解的得分為0.6268，其調(diào)度方案如表3所示。在表3中，每個(gè)應(yīng)急船隊(duì)運(yùn)輸應(yīng)急物資的路徑由其抵達(dá)的節(jié)點(diǎn)順序表示。以應(yīng)急船隊(duì)1為例，其運(yùn)輸路徑為21- 6- 8-13-3-21，表示該船隊(duì)從其所屬的大連基地(節(jié)點(diǎn)21)出發(fā)，依次向需求點(diǎn)6、8、13以及3運(yùn)輸應(yīng)急物資，最終回到大連基地(節(jié)點(diǎn)21)，完成其應(yīng)急物資運(yùn)輸任務(wù)。每個(gè)船隊(duì)的配船情況以及應(yīng)急物資(圍油欄)的分配情況也可由表3獲知。此時(shí)生態(tài)環(huán)境損失為8.902×108元，應(yīng)急響應(yīng)成本為1.323×106元。

表3 動(dòng)態(tài)環(huán)境下的理想調(diào)度方案

5 結(jié)論

本文通過對(duì)油膜動(dòng)態(tài)過程進(jìn)行分析，構(gòu)建了動(dòng)態(tài)多目標(biāo)選址-路徑規(guī)劃模型以解決溢油應(yīng)急物資調(diào)度優(yōu)化問題。然后根據(jù)模型特點(diǎn)，設(shè)計(jì)了IMWOA算法對(duì)其進(jìn)行求解。仿真案例分析結(jié)果表明，該動(dòng)態(tài)模型能良好地反映油膜運(yùn)動(dòng)與應(yīng)急決策之間的相互作用關(guān)系；在算法方面，與NSGA-II相比，IMWOA搜索能力更強(qiáng)，能保障模型的求解效率。

在進(jìn)一步研究中，將充分考慮不同應(yīng)急物資的使用條件與儲(chǔ)運(yùn)特點(diǎn)以探討多種物資協(xié)同調(diào)度情景下的溢油應(yīng)急物資調(diào)度問題。