海上溢油應(yīng)急物資調(diào)度雙目標(biāo)規(guī)劃模型*

李 晶，徐文輝，李亞東

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

0 引言

自20世紀(jì)70年代以來，我國(guó)共發(fā)生溢油事故約2 200起，總溢油量約22 000 t，溢油事故給海洋生態(tài)環(huán)境造成巨大破壞,同時(shí)帶來巨大經(jīng)濟(jì)損失。因此，在面對(duì)突發(fā)溢油事故時(shí)，政府應(yīng)合理選擇出救點(diǎn)，恰當(dāng)規(guī)劃應(yīng)急物資調(diào)度方案，降低溢油事故造成的不利影響。

針對(duì)海上應(yīng)急事故特點(diǎn)進(jìn)行物資調(diào)度的研究尚處于起步階段：王軍等[1]考慮海上環(huán)境等因素對(duì)海上突發(fā)公共安全事故的影響，構(gòu)建以滿足物資供應(yīng)為目標(biāo)的應(yīng)急物資調(diào)度模型；文獻(xiàn)[2-3] 考慮到需求點(diǎn)對(duì)應(yīng)急物資的需求不確定性，用三角模糊數(shù)描述應(yīng)急物資需求量等不確定因素。

相對(duì)海上事故，陸上突發(fā)公共安全事件物資調(diào)度研究比較成熟：文獻(xiàn)[4-8]根據(jù)應(yīng)急事故特點(diǎn)提出“多出救點(diǎn)和單一需求點(diǎn)”的應(yīng)急物資調(diào)度模型；針對(duì)部分突發(fā)事故有多個(gè)需求點(diǎn)，文獻(xiàn)[9-11]提出“多出救點(diǎn)、多需求點(diǎn)及需求點(diǎn)對(duì)應(yīng)急物資的需求量確定”；文獻(xiàn)[12-15]提出“多出救點(diǎn)、多需求點(diǎn)及需求點(diǎn)對(duì)應(yīng)急物資的需求量不確定”的應(yīng)急物資調(diào)度優(yōu)化模型。

綜上，在海上突發(fā)溢油事故情況下，現(xiàn)有研究多采用模糊規(guī)劃或隨機(jī)規(guī)劃描述需求不確定性，而模糊規(guī)劃則需要根據(jù)決策者經(jīng)驗(yàn)得到不確定參數(shù)的模糊隸屬函數(shù)，具有主觀性；隨機(jī)規(guī)劃需要假設(shè)不確定因素服從某種概率分布，結(jié)果不夠準(zhǔn)確。

突發(fā)溢油事故下，應(yīng)急物資應(yīng)以最短時(shí)間抵達(dá)需求點(diǎn)，但在實(shí)際應(yīng)急調(diào)度過程中決策者通常需要考慮應(yīng)急成本，從而選擇合適出救點(diǎn)和救援路徑，2個(gè)目標(biāo)相互沖突，決策者需要在2者之間平衡利弊，以尋求最佳溢油應(yīng)急方案。因此，本文充分考慮海上溢油事故特點(diǎn)，采用魯棒優(yōu)化方法處理受災(zāi)點(diǎn)對(duì)物資需求量的不確定，將溢油點(diǎn)漂移特性納入應(yīng)急物資調(diào)度模型，構(gòu)建應(yīng)急成本最小、應(yīng)急時(shí)間最短的物資調(diào)度雙目標(biāo)規(guī)劃模型，并用NSGA-Ⅱ算法進(jìn)行求解。

1 問題描述

溢油事故發(fā)生在風(fēng)速較大、水流湍急的海面上時(shí)，油膜發(fā)生漂移，形成新的“需求點(diǎn)”，即漂移點(diǎn)。應(yīng)急基地救援示意如圖1所示，其中實(shí)線為救援路徑，虛線為油膜漂移路徑。

圖1 應(yīng)急基地救援示意Fig.1 Schematic diagram of emergency base rescue

為將溢油事故造成的損失降至最低，考慮溢油量不確定，首先選擇到達(dá)“漂移點(diǎn)”時(shí)間最短的應(yīng)急基地，以阻斷油膜漂移路徑，使“漂移點(diǎn)”地理位置不發(fā)生變化，同時(shí)對(duì)溢油點(diǎn)和漂移點(diǎn)展開后續(xù)救援。本文主要針對(duì)需求不確定情況下，多出救點(diǎn)、多需求點(diǎn)的應(yīng)急物資調(diào)度問題展開研究。

針對(duì)事故特點(diǎn)，做出以下4點(diǎn)假設(shè)：

1)溢油事故發(fā)生后，海風(fēng)風(fēng)向和水流方向不變，風(fēng)速等級(jí)和水流速度大小不變。

2)應(yīng)急基地可以供應(yīng)多個(gè)需求點(diǎn)物資，但基地之間不可以進(jìn)行物資轉(zhuǎn)運(yùn)。

3)漂移點(diǎn)一旦受到首批應(yīng)急物資的救助，其地理位置將不會(huì)發(fā)生變動(dòng)。

4)引用“救助資源包”這一概念，將溢油所需救援物資進(jìn)行整合，組成1個(gè)救助資源包，且單個(gè)“救助資源包”可處理1 kg的油污。

2 模型構(gòu)建

2.1 “漂移點(diǎn)”確定

以假設(shè)3)為依據(jù)，在首批應(yīng)急物資到達(dá)漂移點(diǎn)之前，對(duì)油膜漂移動(dòng)態(tài)過程進(jìn)行求解，確定漂移點(diǎn)位置。

根據(jù)油膜漂移速度、船舶速度及應(yīng)急基地至溢油點(diǎn)距離，利用式(1)求得溢油點(diǎn)漂移時(shí)間，以確定溢油點(diǎn)漂移位置：

(D3)2=(D1)2+(D2)2-2D1D2cos(θ)

(1)

式中：D1為油膜漂移距離，m；D2為應(yīng)急基地與溢油點(diǎn)之間距離，m；D3為應(yīng)急基地與漂移點(diǎn)間距離，m；θ為基地至溢油點(diǎn)路徑與漂移路徑兩邊夾角，(°)。

船舶在水中實(shí)際移動(dòng)速度如式(2)所示：

Vz=V-(0.745h-0.257qh)(1.0-1.35×10-6DV)

(2)

式中：Vz為水中實(shí)際行駛速度，kn；V為船舶在靜水中的速度，kn；h為波高，m；q為船首方向與波浪來向之夾角(°)；D為船舶排水量，t。

油膜漂移速度如式(3)所示：

Vo=K×Vw+Vf

(3)

式中：K為風(fēng)力系數(shù)，取0.035；VW為風(fēng)速，kn；Vf為水的表面流速，kn。

油膜漂移距離如式(4)所示：

(4)

式中：Vo為油膜漂移速度，kn；Δt為油膜漂移時(shí)間，h。

2.2 考慮漂移點(diǎn)的物資調(diào)度模型構(gòu)建

溢油事故發(fā)生后，若油污不能及時(shí)清理，將在風(fēng)力和水流共同作用下發(fā)生漂移，擴(kuò)散至更廣闊的海域，以應(yīng)急時(shí)間最少為目標(biāo)函數(shù)，同時(shí)考慮救援成本，構(gòu)建應(yīng)急成本最小和應(yīng)急時(shí)間最少的雙目標(biāo)需求不確定應(yīng)急物資調(diào)度模型如式(5)～(15)所示。其中成本函數(shù)如式(5)所示：

(5)

式中：C為應(yīng)急成本，元；Dij為應(yīng)急基地i與需求點(diǎn)j之間的距離，m；xijs為是否使用救助船舶s將應(yīng)急物資從i點(diǎn)運(yùn)至j點(diǎn)；f為船舶每海里行駛的油耗成本,元；z為隨應(yīng)急救援時(shí)間延長(zhǎng)，溢油導(dǎo)致的環(huán)境污染治理成本，元,本文取值1 000 000。時(shí)間函數(shù)應(yīng)急時(shí)間T如式(6)所示：

(6)

模型約束條件如式(7)～(15)所示：

ri≥0,?i∈I

(7)

Dij≥0,?i∈I,j∈J

(8)

(9)

(10)

(11)

xijs=0，i∈I,j∈I,s∈S

(12)

xi={0,1},?i∈I

(13)

xijs={0,1},?i∈I,j∈J,s∈S

(14)

yijs≥0,?i∈I,j∈J,s∈S

(15)

模型中式(5)～(6)為目標(biāo)函數(shù)，式(5)為成本函數(shù)，表示各出救點(diǎn)向需求點(diǎn)運(yùn)送物資的總運(yùn)輸費(fèi)用與油膜給環(huán)境帶來的損失之和最??；式(6)為時(shí)間函數(shù)，即各船從出救點(diǎn)向需求點(diǎn)的總救援時(shí)間最短。

式(7)～(15)為約束條件：約束(7)指應(yīng)急基地i的救援物資庫存量不小于0；約束(8)指應(yīng)急基地i與需求點(diǎn)j的距離不小于0；約束(9)指從各應(yīng)急基地運(yùn)至需求點(diǎn)的應(yīng)急物資量不超過其庫存；約束(10)指船舶的容量約束；約束(11)指充分滿足需求點(diǎn)的物資需求量；約束(12)指各應(yīng)急基地間的救援物資不可以相互轉(zhuǎn)運(yùn)；約束(13)是0～1變量，如果選擇應(yīng)急基地i為出救點(diǎn)，則xi為1，否則為0；約束(14)是0～1變量，如果選擇使用船舶s將物資從基地i運(yùn)至需求點(diǎn)j，則xijs為1，否則為0；約束(15)指船舶s運(yùn)輸?shù)奈镔Y量不小于0。

2.3 基于魯棒優(yōu)化的需求量不確定性處理

海上突發(fā)溢油事故的溢油量具有不確定性，使用機(jī)會(huì)約束規(guī)劃模型和模糊規(guī)劃描述應(yīng)急物資需求的不確定性具有很大的局限性，但魯棒優(yōu)化具有以下2個(gè)優(yōu)點(diǎn)：魯棒優(yōu)化在建模過程中充分考慮了不確定性，且不需要不確定參數(shù)的分布模型或其模糊隸屬函數(shù)；魯棒優(yōu)化的約束條件是嚴(yán)格成立的，即只要不確定參數(shù)屬于不確定集合，所求出的解就滿足約束條件。因此，本文使用魯棒優(yōu)化方法處理需求不確定性。

(16)

式中：dj為需求點(diǎn)j的名義需求量；aj為需求點(diǎn)j的需求擾動(dòng)量，aj=εjdj；εj為擾動(dòng)比例；uj為不確定性因子。

不確定需求集合如式(17)所示：

(17)

式中：Γu表示需求不確定集合的不確定水平參數(shù)，客觀衡量約束條件的保守程度，體現(xiàn)決策者的風(fēng)險(xiǎn)偏好程度，Γu越小，表明模型越不保守，Γu∈[0,n]；n為模型中不確定需求點(diǎn)的個(gè)數(shù)(n可取小數(shù))[15]。將約束(11)進(jìn)行魯棒優(yōu)化得到約束(18)～(19)：

(18)

(19)

因dj是名義需求，不受不確定性影響，所以將約束(18)轉(zhuǎn)化成約束(20)：

(20)

因約束(20)是非線性規(guī)劃，不易求解，引入對(duì)偶變量yj，zj進(jìn)行對(duì)偶轉(zhuǎn)化，將約束(19)～(20)轉(zhuǎn)化為約束(21)～(23)。

(21)

yj+zj≥εjdj,?j∈J

(22)

yj≥0,zj≥0,?j∈J

(23)

約束條件(11)經(jīng)魯棒優(yōu)化處理后，轉(zhuǎn)化為約束條件如(21)～(23)所示。

3 NSGA-Ⅱ算法設(shè)計(jì)

3.1 NSGA-Ⅱ

針對(duì)雙目標(biāo)應(yīng)急物資調(diào)度模型，本文用NSGA-Ⅱ算法求解，由GA改進(jìn)而來，需經(jīng)歷種群初始化、選擇、交叉、變異、迭代等過程。與GA最大區(qū)別在于在算法中加入非支配排序和擁擠度比較；引入精英策略，保證個(gè)體的多樣性。

3.2 編碼

本文使用二進(jìn)制編碼，根據(jù)問題特點(diǎn)及規(guī)模，生成的染色體有22個(gè)基因，如圖2所示。將染色體分成11部分，每部分代表1個(gè)出救點(diǎn)，且有2個(gè)基因，基因位1代表是否將物資運(yùn)至溢油點(diǎn)；基因位2代表是否將物資運(yùn)至漂移點(diǎn)；若2個(gè)基因均為1，表示既將物資運(yùn)至溢油點(diǎn)，也將物資運(yùn)至漂移點(diǎn)。圖2中染色體表示將應(yīng)急物資從基地2運(yùn)至溢油點(diǎn)，從基地5運(yùn)至漂移點(diǎn)，并且從基地10運(yùn)往溢油點(diǎn)和漂移點(diǎn)。

圖2 染色體編碼示意Fig.2 Schematic diagram of chromosome coding

3.3 種群初始化、選擇、交叉及變異

1)種群初始化。本文種群中染色體為隨機(jī)生成，每個(gè)應(yīng)急基地都有幾率將物資運(yùn)往需求點(diǎn)，但隨機(jī)產(chǎn)生的染色體應(yīng)滿足模型中的約束條件。

2)選擇。采用錦標(biāo)賽選擇策略，基于擁擠比較算子及非支配排序等級(jí)。擁擠比較算子在算法中可以將選擇過程導(dǎo)向均勻分布的Pareto最優(yōu)前沿面，在2個(gè)不同支配等級(jí)排名的解中，選取排名較低的解，若2個(gè)解屬于同一Pareto前沿面，選取處于非擁擠區(qū)域的解。

3)交叉。選取鄰近的2個(gè)父代染色體尋找交叉點(diǎn)，進(jìn)行單點(diǎn)交叉，單點(diǎn)交叉能較好地保存父親的優(yōu)良基因，假如取4號(hào)基因位為交叉位置，染色體A和染色體B交叉所得染色體為A′和B′，如式(24)所示：

(24)

4)變異。本文采取的變異方式為隨機(jī)變異，在眾多染色體中，以一定概率選擇需要變異的染色體，隨機(jī)選擇該染色體中的1個(gè)或多個(gè)基因進(jìn)行變異。例如，選擇染色體C的第2個(gè)和第7個(gè)基因進(jìn)行變異，形成新的染色體為C′，如式(25)所示：

(25)

4 算例分析

4.1 相關(guān)數(shù)據(jù)

本文假設(shè)蓬萊溢油事件名義溢油量為250 t；根據(jù)第4)條假設(shè)可知“救助資源包”的名義需求量為250 000個(gè)；假設(shè)溢油點(diǎn)對(duì)“救助資源包”的名義需求量為150 000個(gè)，“漂移點(diǎn)”的名義需求量為100 000個(gè)；根據(jù)商務(wù)部數(shù)據(jù)顯示，救助船舶行駛每海里的成本f為315元。各備選應(yīng)急基地基本信息見表1。

表1 備選應(yīng)急基地基本信息Table 1 Basic information of alternative emergency bases

2011年蓬萊石油泄漏事故發(fā)生時(shí)，渤海海域風(fēng)向?yàn)槠巷L(fēng)，平均風(fēng)速為2.214 kn，水流流速較小，取船首與浪向夾角q為0°，水流流速Vf取0 kn[16]，浪高0.6 m；應(yīng)急救助船舶平均航行速度為21 kn，總噸位D=4 000 t(噸)，托運(yùn)力100 t。利用式(2)計(jì)算船舶實(shí)際行駛速度為20.164 kn，式(3)計(jì)算出油膜漂移速度為0.279 kn。

4.2 模型結(jié)果

1)“漂移點(diǎn)”位置確定

“漂移點(diǎn)”位置確定主要包括以下3個(gè)步驟：

①利用應(yīng)急基地和溢油點(diǎn)經(jīng)緯度(38.373N，120.118E)求出其之間的距離。

②由于溢油事故發(fā)生時(shí)風(fēng)向?yàn)槟巷L(fēng)，且水流對(duì)油膜漂移的影響可忽略不計(jì)，可知漂移點(diǎn)的經(jīng)度與溢油點(diǎn)的經(jīng)度一致。首先假設(shè) “擬定漂移點(diǎn)”的經(jīng)緯度為(39.373N，120.118E)，然后利用ArcGIS軟件將溢油點(diǎn)、擬定漂移點(diǎn)和應(yīng)急基地作為三角形3個(gè)頂點(diǎn)，最后求出以溢油點(diǎn)位頂點(diǎn)的夾角θ。

③將船舶實(shí)際速度、油膜漂移速度、夾角θ及應(yīng)急基地與溢油點(diǎn)之間的距離帶入式(1)，得到應(yīng)急基地至漂移點(diǎn)的時(shí)間，結(jié)果顯示東營(yíng)至漂移點(diǎn)的時(shí)間最短為2.838 h。

因此，首先選擇將存放在東營(yíng)應(yīng)急基地的物資運(yùn)往漂移點(diǎn)，以阻礙油膜繼續(xù)漂移，減少對(duì)生態(tài)環(huán)境的破壞，利用式(4)求得溢油點(diǎn)漂移的最長(zhǎng)距離為1 466.784 m，求得漂移點(diǎn)的經(jīng)緯度為(38.387145 N，120.118486 E)。

2)調(diào)度結(jié)果

本文用Matlab2016a進(jìn)行計(jì)算，而需求點(diǎn)對(duì)“救助資源包”的實(shí)際需求量取決于決策者對(duì)風(fēng)險(xiǎn)水平的偏好程度，其在魯棒模型中可以從Γu和εj的取值大小體現(xiàn)出來。由于本文的需求點(diǎn)只有2個(gè)，因此，Γu和εj在本文中的取值集合分別為{0，0.5，1，1.5，2}和{3%，5%，10%}；取種群規(guī)模為200，迭代次數(shù)為1 000，得到的調(diào)度結(jié)果見表2。當(dāng)Γ為0時(shí)，無論擾動(dòng)比例取3%、5%、10%，其所代表的方案均為：3-a,7-b,8-a,9-b,11-a。表示將應(yīng)急物資從營(yíng)口運(yùn)至溢油點(diǎn)、萊州運(yùn)至漂移點(diǎn)、濰坊運(yùn)至溢油點(diǎn)、東營(yíng)運(yùn)至漂移點(diǎn)和葫蘆島運(yùn)至溢油點(diǎn)。

表2 物資調(diào)度方案Table 2 Schemes of materials dispatching

注：表中數(shù)字所代表的基地序號(hào)見表1；a，b分別指溢油點(diǎn)和漂移點(diǎn)。

應(yīng)急調(diào)度具有弱經(jīng)濟(jì)性的特點(diǎn)，勢(shì)必以效率優(yōu)先，但實(shí)際過程中會(huì)受嚴(yán)格資金約束限制，所以帕累托解可為決策者提供參考，決策者可根據(jù)實(shí)際應(yīng)急事故，選擇相應(yīng)應(yīng)急調(diào)度方案，名義模型結(jié)果示意如圖3所示。由圖3可知，得到的解具有前沿面，若決策者判斷此事故為重大事故，則選擇應(yīng)急時(shí)間較短的調(diào)度方案，否則選擇應(yīng)急成本和應(yīng)急時(shí)間綜合結(jié)果較優(yōu)的應(yīng)急物資調(diào)度方案。

圖3 名義模型結(jié)果Fig.3 Results of nominal model

3)結(jié)果分析

雖然不同擾動(dòng)比例及不確定水平下的應(yīng)急成本和時(shí)間不同，但應(yīng)急出救點(diǎn)的選擇方案共有8個(gè)，其中某個(gè)方案被使用了8次，說明該模型的魯棒性較好，優(yōu)化方案對(duì)參數(shù)擾動(dòng)不敏感。

根據(jù)模型，不同風(fēng)險(xiǎn)類型偏好的決策者可以得到相應(yīng)的救援方案，方案見表2。如果決策者為風(fēng)險(xiǎn)追求偏好，則可選取較小的Γu和εj，得到相應(yīng)的應(yīng)急方案，該方案應(yīng)急成本較小并且應(yīng)急時(shí)間較少；如果決策者為風(fēng)險(xiǎn)規(guī)避型，則可選取較大的Γu和εj，但是會(huì)以較高的應(yīng)急成本和較長(zhǎng)的應(yīng)急時(shí)間為代價(jià)；如果決策者可接受的風(fēng)險(xiǎn)水平介于2者之間，則選擇折中方案。在不同擾動(dòng)比例εj和不確定水平Γu下，2個(gè)函數(shù)的目標(biāo)值變化如圖4所示。當(dāng)Γu=0時(shí)，魯棒物資調(diào)度模型等價(jià)于名義模型，應(yīng)急成本C=742 770.8元，應(yīng)急時(shí)間T=24.31 h。Γu相同時(shí)，隨εj變大，2個(gè)目標(biāo)函數(shù)值上升幅度呈現(xiàn)從小到大的趨勢(shì)；當(dāng)εj為3%時(shí)，隨Γu增大，目標(biāo)函數(shù)上升趨勢(shì)并不明顯；當(dāng)εj為5%時(shí)，隨Γu增大，目標(biāo)函數(shù)值呈現(xiàn)較慢的上升趨勢(shì)；當(dāng)εj為10%時(shí)，隨Γu的增大，目標(biāo)函數(shù)上升趨勢(shì)比較明顯，應(yīng)急成本相對(duì)名義模型的增加率分別為0.6%、1.7%、2.3%和5.1%，應(yīng)急時(shí)間相對(duì)名義模型的增加率分別為1.3%、5.6%、10.2%和21.4%。Γu和εj在一定程度上可以衡量決策者的風(fēng)險(xiǎn)偏好，因此決策者可根據(jù)自己對(duì)不確定風(fēng)險(xiǎn)水平的偏好程度選擇相應(yīng)的Γu和εj。

圖4 需求不確定下的目標(biāo)函數(shù)值Fig.4 Objective function values under uncertainty of demand

在不同Γu和εj情況下獲得的物資調(diào)度方案見表2。隨Γu和εj的變大，為滿足溢油點(diǎn)和漂移點(diǎn)對(duì)應(yīng)急物資的需求，采取以下3種措施：①當(dāng)原有方案能滿足現(xiàn)在的需求，則原來的方案保持不變。②當(dāng)原有出救點(diǎn)能滿足需求，但物資調(diào)度分布的數(shù)量不能滿足受災(zāi)點(diǎn)對(duì)應(yīng)急物資的需求量，則應(yīng)急出救點(diǎn)的選擇保持不變，將調(diào)度方案重新規(guī)劃。③當(dāng)原有應(yīng)急物資調(diào)度方案不能滿足需求點(diǎn)對(duì)應(yīng)急物資的需求時(shí)，系統(tǒng)會(huì)重新進(jìn)行計(jì)算評(píng)估獲得新的救援方案。

5 結(jié)論

1)充分考慮油膜漂移特點(diǎn)和需求不確定性因素，基于魯棒優(yōu)化思想建立以救援應(yīng)急成本最小和應(yīng)急時(shí)間最短為目標(biāo)的規(guī)劃模型，并運(yùn)用NSGA-Ⅱ算法進(jìn)行求解。

2)以蓬萊溢油事件為背景，設(shè)計(jì)具體算例，運(yùn)用所構(gòu)建模型確定應(yīng)急物資調(diào)度方案，驗(yàn)證模型有效性。

3)針對(duì)需求點(diǎn)對(duì)應(yīng)急物資需求不確定性，使用魯棒優(yōu)化方法為不同風(fēng)險(xiǎn)類型的決策者提供相應(yīng)應(yīng)急物資調(diào)度方案。

4)本文僅考慮油膜的1次性漂移，在后續(xù)研究中，將會(huì)把油膜漂移的動(dòng)態(tài)過程納入研究范圍，以期得到的應(yīng)急物資調(diào)度方案更加精準(zhǔn)有效。