危險品車輛路徑規(guī)劃的雙目標(biāo)模型與算法研究

沈 良,吳 婷,覃朝暉,馮 雨,邵 虎,邵 楓

(1.徐州醫(yī)科大學(xué)管理學(xué)院,江蘇 徐州 221116)(2.南京大學(xué)數(shù)學(xué)系,江蘇 南京 210093)(3.中國煤炭運銷協(xié)會,北京 100013)(4.中國礦業(yè)大學(xué)數(shù)學(xué)學(xué)院,江蘇 徐州 221116)

隨著城市化進(jìn)程的不斷發(fā)展,我國對于危險品的需求也逐漸增加,這就導(dǎo)致危險品運輸?shù)氖袌鲂枨笞兊猛? 需要指出的是危險物品具有易爆性、有毒性、傳染性、腐蝕性和其他危險的特點,一旦在運輸過程中出現(xiàn)交通事故導(dǎo)致危險品的泄露,往往容易造成人員傷亡、財產(chǎn)損失、環(huán)境污染等. 據(jù)估計,每年全世界有40億t危險材料被運輸. 中國每年通過公路運輸?shù)奈ｋU品超過10億t,全國約有1.16萬家企業(yè),約 3 200 萬人參與,從事道路危險品運輸?shù)能囕v超過36萬輛[1]. 對于承運企業(yè)而言,行程時間和運輸風(fēng)險是兩個非常重要的考慮因素. 從行程時間的角度看,承運企業(yè)希望盡可能地在較短的時間內(nèi)完成運輸任務(wù),降低企業(yè)的成本,這種情況下可能就會選擇經(jīng)過城區(qū)中心的短距離的運輸路線,此時就會使得較多的人口面臨潛在的風(fēng)險. 從運輸風(fēng)險的角度看,承運企業(yè)希望盡可能地降低運輸風(fēng)險,這種情況下可能就會選擇距離較長的郊區(qū)道路進(jìn)行運輸,此時就會使得行程時間變長,增加企業(yè)的運輸成本. 因此,如何更好地平衡這兩個因素之間的關(guān)系是一個非常值得研究的問題.

對于危險品運輸車輛的行程時間,可以看成是危險品運輸車輛的路徑優(yōu)化問題,此問題已經(jīng)引起了各國政府和主管部門的廣泛關(guān)注,許多學(xué)者從不同的方面進(jìn)行了研究[2-3]. 如Joy[4]在危險品運輸中運用最短路徑算法進(jìn)行了一些實證研究. Hall[5]研究了在出行時間具有隨機(jī)性和時變性的網(wǎng)絡(luò)中的最優(yōu)路徑問題. Toumazis等[6]提出了一種基于條件風(fēng)險值測度的車輛交通網(wǎng)絡(luò)的危險品路徑優(yōu)化方法. Wei等[7]在運輸風(fēng)險為時變模糊隨機(jī)變量的假設(shè)下,研究了危險品運輸?shù)奈恢谜{(diào)度問題,并給出了車輛的最優(yōu)出發(fā)時間. Fu等[8]在隨機(jī)的動態(tài)交通網(wǎng)絡(luò)條件下,提出了一種基于K最短路徑算法的啟發(fā)式算法用來求解動態(tài)隨機(jī)最短路徑問題. Chang等[9]在此基礎(chǔ)上,將此問題推廣到多準(zhǔn)則路徑優(yōu)化問題上,并給出詳細(xì)的應(yīng)用. 但是上述的研究都沒有考慮到危險品車輛運輸?shù)男谐虝r間可靠性[10-11]. 代存杰等以行程時間和運輸風(fēng)險的隨機(jī)屬性值為優(yōu)化準(zhǔn)則,建立了整數(shù)規(guī)劃模型,考慮了不確定條件下的危險品運輸車輛的可靠路徑和運輸風(fēng)險[12]. 然而,該研究中并沒有考慮到危險品運輸車輛在交通信號交叉口的等待時間以及路口和路段之間的時間相關(guān)性.

對于危險品運輸車輛的運輸風(fēng)險,其目的是評估事故的可能性和事故后果的嚴(yán)重程度,以尋求最低事故率和最少的損失[13]. 在現(xiàn)實的交通運輸網(wǎng)絡(luò)中,關(guān)于危險品運輸車輛路徑屬性值的相關(guān)數(shù)據(jù)的收集和采集的難度是非常大的,如發(fā)生意外時受影響的人數(shù)和受影響的車輛數(shù)、危險品運輸車輛的行駛時間等都具有很強(qiáng)的隨機(jī)性(并不是一個確定性的值,即數(shù)據(jù)的不確定性). 在這種情況下,更加需要進(jìn)行不確定條件下的危險品車輛運輸風(fēng)險的研究,使得結(jié)果更加符合實際情況. 需要注意的是,目前關(guān)于危險品車輛運輸風(fēng)險的計算方式有很多種,其中最為常用的就是危險品運輸車輛的傳統(tǒng)風(fēng)險度量模型[14]. 代存杰等[12]在傳統(tǒng)風(fēng)險模型的基礎(chǔ)上,將不同路段附近的人口數(shù)利用對數(shù)正態(tài)分布函數(shù)進(jìn)行分段擬合,可得不同路段附近人口的均值和方差,并建立多準(zhǔn)則路徑優(yōu)化模型. 接著,在此基礎(chǔ)上,代存杰等[15]建立了基于運輸風(fēng)險的隨機(jī)機(jī)會約束規(guī)劃模型,并利用分段線性逼近方法將非線性函數(shù)轉(zhuǎn)化為線性函數(shù)進(jìn)行求解. 需要注意的是,以上的研究中都只是考慮到了路段附近的人口數(shù),并沒有考慮路段上行駛的車輛中的人數(shù),當(dāng)危險品運輸車輛發(fā)生交通事故后,該路段上車輛中的駕駛員和乘客受到的影響會更嚴(yán)重.

本文擬對隨機(jī)道路交通網(wǎng)絡(luò)中危險品運輸車輛的可靠路徑和安全路徑進(jìn)行研究,建立一個雙目標(biāo)優(yōu)化模型,考慮到目標(biāo)函數(shù)的不可加性和非線性,采用基于K短路算法的啟發(fā)式算法進(jìn)行模型的求解,利用經(jīng)典的交通網(wǎng)絡(luò)算例對模型的準(zhǔn)確性和有效性進(jìn)行驗證.

1 問題描述與建模

1.1 基本符號定義

G=(N,A) 表示有向路網(wǎng),其中N表示節(jié)點的集合,A表示路段的集合;

Ξ表示路段和路口之間的相關(guān)系數(shù)矩陣;

ρa(bǔ)(i,j),a(s,t)表示路段a(i,j)和a(s,t)上的受影響人數(shù)的相關(guān)系數(shù);

Mk和Dk表示路徑k的近似對數(shù)正態(tài)分布的一階原點矩和二階中心矩;

Rk表示危險品運輸車輛在路徑k上總的運輸風(fēng)險;

Pa(i,j)表示兩個節(jié)點i和j之間的位置關(guān)系的變量;

γa(i,j)表示事故發(fā)生后路段a(i,j)附近受到影響的居民數(shù)量;

va(i,j)表示事故發(fā)生后路段a(i,j)上的車輛中受影響人數(shù);

Φ-1(θ) 表示置信水平θ下路徑通行時間的標(biāo)準(zhǔn)正態(tài)分布的逆累積分布函數(shù)(CDF);

pa(i,j)表示危險品運輸車輛在路段a(i,j)上發(fā)生事故的概率;

la(i,j)表示路段a(i,j)的長度;

βa(i,j)表示路段a(i,j)上的出行時間;

ψijl表示危險品車輛在節(jié)點j處的轉(zhuǎn)向;

fRk表示路徑k上的有效風(fēng)險值;

κa(i,j),a(j,l)表示運輸車輛在路口j的等待時間;

1.2 可靠路徑定義

(1)

(2)

(3)

(4)

需要注意的是,根據(jù)Sun等[21]和Yu等[22-23]關(guān)于隨機(jī)道路交通網(wǎng)絡(luò)中的交通信號燈的研究,可以發(fā)現(xiàn)若不考慮車輛在紅綠燈路口的等待時間,會使得求解出的結(jié)果與實際情況存在較大偏差.因此,需要將出行時間的不確定性、車輛在交通信號燈路口的等待時間以及不同路段之間的相關(guān)性等因素考慮到路徑優(yōu)化模型中,并進(jìn)行深入的研究.那么,運輸車輛在路口的等待時間參考Shen等[20]的模型,可以得到路徑k上的所有路口等待時間的均值和方差,分別為式(5)、(6)、(7)和(8).

(5)

(6)

(7)

(8)

式中,Pa(i,j)表示的是關(guān)于兩個節(jié)點i和j之間的位置關(guān)系的變量.Pa(i,j)=1,-1,2,-2分別表示點i在點j的上方、下方、左方和右方(點i和點j之間存在可行路段).本文采用ψijl=Pa(i,j)?Pa(j,l)表示危險品運輸車輛在節(jié)點j處的轉(zhuǎn)向(即通過路段a(i,j)和a(j,l),包含左轉(zhuǎn)、直行和右轉(zhuǎn)),具體的計算過程見式(9),可以發(fā)現(xiàn)ψijl=1表示車輛在節(jié)點j處左轉(zhuǎn),ψijl=2表示車輛在節(jié)點j處直行,ψijl=3表示車輛在節(jié)點j處右轉(zhuǎn).

(9)

圖1 節(jié)點2處的右轉(zhuǎn)示意圖Fig.1 Turn right at intersection node 2

為了說明車輛在路口轉(zhuǎn)彎的方向,本文只給出一個示意圖來描述右轉(zhuǎn),同理可得左轉(zhuǎn)和直行.如圖1所示,共有3個節(jié)點,分別為始點1、路口節(jié)點2和終點3.可以發(fā)現(xiàn),始點1在路口節(jié)點的左邊,故Pa(1,2)=2,中間節(jié)點2 在終點3的上方,故Pa(2,3)=1,即ψ123=2?1.此時根據(jù)等式(9)中的描述,可以發(fā)現(xiàn)車輛在節(jié)點2處的轉(zhuǎn)彎方向為右轉(zhuǎn)(見圖1).

(10)

(11)

式中,θ表示的是給定的置信水平.接著,可以將等式(11)進(jìn)行變形,得到有效出行時間模型:

(12)

上述給出的有效出行時間模型可以保證在最小化出行時間的同時使得出行時間的可靠性不小于給定的置信水平θ.接著,給出最優(yōu)可靠路徑的定義:

1.3 安全路徑定義

關(guān)于危險品運輸車輛的風(fēng)險度量模型有很多,其中最為常用的是傳統(tǒng)的風(fēng)險度量模型[13-15],具體形式如下:

ra(i,j)=pa(i,j)γa(i,j).

(13)

通常取pa(i,j)=la(i,j)×10-6.需要注意的是,當(dāng)事故發(fā)生后,不僅路段附近的居民會受到影響,路段上的車輛中的人也會受到影響,并且車輛中的人受到的影響會比附近的居民受到的影響更大.考慮到受影響的車輛中的人數(shù)與路段上的車流量是相關(guān)的,因此可以通過交通配流模型[25-26]得到不同路段上的車流量,進(jìn)而進(jìn)行運算.據(jù)此,給出如下新的風(fēng)險分析模型:

ra(i,j)=pa(i,j)(γa(i,j)+va(i,j)).

(14)

由于交通網(wǎng)絡(luò)中不同路段上的人口密度具有差異性和不確定性,本文中假設(shè)變量γa(i,j)是服從對數(shù)正態(tài)分布的隨機(jī)數(shù)[8,12],可得受影響的總?cè)藬?shù)的計算公式如下:

(15)

據(jù)此,可得危險品運輸車輛在路徑k上總的運輸風(fēng)險Rk為:

(16)

通過以上的討論,發(fā)現(xiàn)γa(i,j)是一個隨機(jī)變量,那么可得Rk也是一個隨機(jī)變量,為了控制運輸風(fēng)險的大小,建立了一個隨機(jī)機(jī)會約束規(guī)劃模型[12],給出有效風(fēng)險值fRk的定義,具體如下.

有效風(fēng)險值由概率約束模型確定:

(17)

根據(jù)概率約束模型確定的風(fēng)險值fR(x)即為滿足機(jī)會約束的有效風(fēng)險值.通過簡單的變形,可得對數(shù)正態(tài)分布下的有效風(fēng)險值為:

(18)

(19)

(20)

(21)

(22)

(23)

1.4 雙目標(biāo)優(yōu)化模型

根據(jù)上述關(guān)于可靠路徑和安全路徑的定義,給出了不確定條件下的危險品運輸車輛的雙目標(biāo)優(yōu)化模型,具體形式如下:

(24)

(25)

ψijl=Pa(i,j)?Pa(j,l),

(26)

(27)

Pa(i,j)∈{-2,-1,1,2}, ?a(i,j)∈A.

(28)

在該模型中,式(24)為該優(yōu)化模型的目標(biāo)函數(shù),即最優(yōu)可靠路徑和最優(yōu)安全路徑;式(25)保證了路徑的可行性;式(26)為運輸車輛在路口拐彎方向的確定;式(27)為路段路徑關(guān)聯(lián)變量;式(28)表示兩個相鄰節(jié)點的位置關(guān)系.

2 求解算法設(shè)計

根據(jù)雙目標(biāo)優(yōu)化模型中的目標(biāo)函數(shù)(24),可以發(fā)現(xiàn)兩個目標(biāo)都是具有不可加性和非線性的特點,是經(jīng)典的NP難問題. 傳統(tǒng)的路徑優(yōu)化算法(如Dijkstra算法、Floyd算法)并不能直接求解本文中提出的優(yōu)化模型.因此,為了克服該模型的不可加性和非線性的特點,本文擬通過采用K短路算法和不等式技巧(放縮法)進(jìn)行模型的求解.首先利用放縮法推導(dǎo)出每條可行路徑的有效出行時間和有效風(fēng)險值的下界.此時,此時放縮后的下界需要滿足可加性.接著,利用K短路算法生成一個備選可靠路徑集.需要注意的是,該備選可靠路徑集合中的路徑數(shù)量取決于選取的參數(shù)“K”的大小(小于等于參數(shù)“K”),其中該集合是有界的(因為K是有界的).最后,通過計算備選可靠路徑集合中的路徑的有效出行時間和有效風(fēng)險值,可以分別得到最優(yōu)可靠路徑和最優(yōu)安全路徑.通過以上算法的簡單描述,本文給出了如下的數(shù)學(xué)原理和證明.

2.1 求解算法的數(shù)學(xué)原理和證明

(29)

(30)

(31)

定理2最優(yōu)可靠路徑k*一定包含在備選路徑集合Q中,如果路徑k*滿足如下條件:

(32)

證明采用反證法.首先假設(shè)存在一條路徑k*,該路徑是OD對rs中所有可行路徑中有效出行時間最小的,并且假設(shè)該路徑k*不存在于備選可靠路徑集合Q中(即k*?Q).接下來,可以發(fā)現(xiàn)定理2中描述了若k*?Q,則:

(33)

根據(jù)定理2中所描述的內(nèi)容,發(fā)現(xiàn)路徑ki是所有路徑中有效出行時間下界中最小的.鑒于此,可以得到在所有的可行路徑中,一定存在某條路徑kj,滿足如下的不等式:

(34)

即路徑kj的有效出行時間小于或者等于路徑ki的有效出行時間.

接著,根據(jù)等式(33)和(34),可得:

(35)

接著,根據(jù)不等式(35),可得:

(36)

此時可以發(fā)現(xiàn)路徑kj的有效出行時間小于路徑k*的有效出行時間.那么可得路徑k*并不是當(dāng)前的最優(yōu)可靠路徑.據(jù)此可發(fā)現(xiàn),得到的結(jié)論與前面的假設(shè)是相互矛盾的.因此,可得路徑k*一定包含在備選可靠路徑集合Q中,定理2證畢.

定理3最優(yōu)安全路徑k**一定包含在備選路徑集合Q中,如果路徑k**滿足如下條件:

(37)

式中,fRki表示具有有效風(fēng)險值的下界中的最小值所對應(yīng)的路徑ki的有效風(fēng)險值.

證明根據(jù)定理2中的證明過程,定理3也是成立的.

2.2 算法的步驟

步驟1 進(jìn)行初始化

步驟2 確定備選路徑集合Q中的路徑

根據(jù)有效出行時間的下界,利用K短路算法[28,30]得到一組以遞增順序排列的路徑.

01:對于第ith短路路徑的有效出行時間的下界

03:如果i

04:然后,跳到步驟2開始的地方,否則

05:執(zhí)行步驟3,結(jié)束

06:否則,執(zhí)行步驟3,結(jié)束

07:結(jié)束

步驟3 確定備選路徑集合P中的路徑

根據(jù)有效風(fēng)險值的下界,利用K短路算法得到一組以遞增順序排列的路徑.

01:對于第jth短路路徑的有效風(fēng)險值的下界

02:如果fRkj,min≤fRkn

03:如果j

04:跳到步驟3開始的地方,否則

05:執(zhí)行步驟4,結(jié)束

06:否則,執(zhí)行步驟4,結(jié)束

07:結(jié)束

步驟4 尋找非支配有效解

對于備選路徑集合Q和P中的所有路徑,執(zhí)行如下操作:

注本文中的算法利用了K短路算法.由于K短路算法具有多項式時間計算復(fù)雜度[30],因此本節(jié)提出的算法也具有多項式時間復(fù)雜度.

圖2 路網(wǎng)結(jié)構(gòu)圖(km)Fig.2 The transportation network(km)

3 數(shù)值算例

為了驗證模型和算法的準(zhǔn)確性和有效性,本文中的數(shù)值算例采用代存杰[12]文章中的路網(wǎng)結(jié)構(gòu)圖,具體圖形見圖2.各路段上的行程時間的參數(shù)值以及區(qū)域內(nèi)人口數(shù)量的分布參數(shù)可以參考論文代存杰[12]中的數(shù)值算例.需要注意的是本文中還考慮了事故發(fā)生后受影響的路段上車輛中的人數(shù),該人數(shù)是根據(jù)路段上的車流量來確定的.根據(jù)經(jīng)典的BPR(bureau of public roads)函數(shù)模型,可以發(fā)現(xiàn)路段上的車流量和該路段上的通行時間具有正相關(guān)性,因此為了方便計算,在數(shù)值算例中通過參數(shù)和路段通行時間相乘的方法來近似模擬路段上受影響的人數(shù),函數(shù)模型中的參數(shù)設(shè)置為θ=0.95,該數(shù)值算例是在內(nèi)存8GB,CPU 3.00GHz的臺式計算機(jī)中,采用MATLAB2020a對算法進(jìn)行編程.

3.1 數(shù)值算例的解的分析

為說明本文提出的模型和算法的有效性,我們將本文中的數(shù)值算例的結(jié)果和代存杰[12]中的數(shù)值算例的結(jié)果進(jìn)行了比較,結(jié)果見表1. 通過表1中的數(shù)據(jù),可看出不同路徑的有效出行時間和有效風(fēng)險值均與文獻(xiàn)中的結(jié)果不同,這是由于本文中的模型考慮了更多的因素,使得結(jié)果更加符合實際情況.

表1 不同路徑上的有效出行時間和有效風(fēng)險值的比較Table 1 Comparisons of effective travel times and effective risk values for different paths

為了說明路口等待時間、路段和路口之間的相關(guān)性以及路段上的車流量、不同路段附近居民數(shù)量之間相關(guān)性是如何影響最優(yōu)可靠路徑和最優(yōu)安全路徑,我們給出了表2中所示的4種情況.

表2 不同情況下的有效出行時間和有效風(fēng)險值Table 2 The effective travel times and effective risk values under different scenarios

從表2中可以發(fā)現(xiàn),最優(yōu)可靠路徑和最優(yōu)安全路徑并不是同一條,這說明了這兩條路徑都是該模型的非支配有效解. 對于最優(yōu)可靠路徑,在不考慮路口等待時間和相關(guān)性的情況下,路徑1-2-6-10-13的有效出行時間為2.88 h,若考慮這些因素,路徑1-2-6-10-13的有效出行時間為3.15 h,可以發(fā)現(xiàn)路口等待時間和相關(guān)性對于有效出行時間影響還是比較大的. 對于最優(yōu)安全路徑,在不考慮路段上的車流量以及相關(guān)性的情況下,路徑1-5-8-10-11-13的有效風(fēng)險值為47.9×10-3人,若考慮上述兩個因素的話,路徑 1-5-8-10-11-13的有效風(fēng)險值變?yōu)?2.78×10-3人,說明路段上的車流量和相關(guān)性對于有效風(fēng)險值的計算有一定的影響. 至此,我們給出的數(shù)值算例說明了路口等待時間,路段和路口之間的相關(guān)性以及路段上的車流量,不同路段附近居民數(shù)量之間相關(guān)性對于提出的雙目標(biāo)優(yōu)化模型的求解結(jié)果具有重要的影響.

3.2 數(shù)值算例的非支配有效解

根據(jù)本文中提出的模型和算法,給出了OD對1和13之間的非支配有效解的具體結(jié)果,見圖3和表3. 根據(jù)圖3中的結(jié)果可以發(fā)現(xiàn)此時共有6個非支配的有效解,具體的非支配有效解對應(yīng)的有效出行時間和有效風(fēng)險值見表3.

圖3 OD對1和13之間的Pareto前沿面Fig.3 The Pareto solutions between OD pair 1 and 13

從表3中可以發(fā)現(xiàn)具有最小有效出行時間的路徑是1-2-6-10-13,具有最小有效風(fēng)險值的路徑是1-5-8-10-11-13. 那么,對于危險品運輸車輛來說,最優(yōu)路徑的選擇完全取決于司機(jī)或者運輸公司領(lǐng)導(dǎo)的個人偏好.

表3 OD對1和13之間的非支配有效解Table 3 The non-dominated solutions between OD pair 1 and 13

4 結(jié)論

本文討論了不確定條件下危險品車輛運輸?shù)目煽柯窂胶桶踩窂降碾p目標(biāo)優(yōu)化模型. 最后通過數(shù)值算例說明了模型和算法的準(zhǔn)確性和有效性,并得到了一系列非支配有效解. 通過本文的一些研究數(shù)值算例,我們論證了如下內(nèi)容:

(1)可靠路徑的有效出行時間取決于不確定網(wǎng)絡(luò)中的路段通行時間的均值和方差、路口等待時間的均值和方差以及準(zhǔn)時到達(dá)的概率和其相關(guān)系數(shù)矩陣;

(2)安全路徑的有效風(fēng)險值取決于不確定網(wǎng)絡(luò)中的路段通行時間(路段上的車流量)、路段的距離、路段附近人口數(shù)的均值和方差以及相關(guān)系數(shù)矩陣;

(3)此雙目標(biāo)優(yōu)化模型的最優(yōu)路徑選擇取決于通行者的個人偏好;

(4)提出的雙目標(biāo)模型及其求解算法可以顯著地幫助危險品車輛運輸者在可靠路徑和安全路徑之間進(jìn)行選擇和預(yù)估,即使在具有不確定性的大型實際交通網(wǎng)絡(luò)中也是適應(yīng)的.

未來研究的方向可以考慮以下幾個方面:本文中提出的基于可靠性的路徑搜索問題是在隨機(jī)道路網(wǎng)絡(luò)中計算有效出行時間的靜態(tài)模型. 因此,如何將本文提出的可靠路徑搜索模型拓展到一個隨時間變化的動態(tài)模型中是一個值得研究的問題. 此外,如何將該論文提出的算法應(yīng)用到更大規(guī)模的交通網(wǎng)絡(luò)中,從而驗證提出模型和算法的有效性是后續(xù)需要研究的內(nèi)容.