吳坤霖，趙佳虹

(1.廣東工業(yè)大學(xué) 土木與交通工程學(xué)院，廣東 廣州 511006；2. 武漢理工大學(xué) 交通與物流工程學(xué)院，湖北 武漢 430063)

0 引言

2020年新冠疫情爆發(fā)，感染性醫(yī)療廢物產(chǎn)量激增，城市醫(yī)療廢物運輸安全面臨極大挑戰(zhàn)[1]。合理評估運輸風(fēng)險、優(yōu)化運輸?shù)脑O(shè)施選址和車輛路徑?jīng)Q策是疫情下感染性醫(yī)療廢物運輸急待解決的安全問題。

醫(yī)療廢物屬于危險廢物范疇，近年來國內(nèi)外學(xué)者在危險廢物運輸風(fēng)險度量和選址-路徑優(yōu)化問題等領(lǐng)域進(jìn)行了大量研究[2]。Yilmaz等[3]設(shè)計了考慮人口和脆弱環(huán)境區(qū)域影響的風(fēng)險模型；Taslimi等[4]以最小化最大區(qū)域風(fēng)險來實現(xiàn)風(fēng)險公平性；Zhao等[5]設(shè)計了不同水環(huán)境中的重油擴散風(fēng)險模型；Yu等[6]設(shè)計了人口感染風(fēng)險模型；趙佳虹等[7]考慮新冠病毒的氣溶膠傳播特性，構(gòu)建了時變“環(huán)境-人口”風(fēng)險模型；Taslimi等[8]建立了醫(yī)療廢物運輸-存放風(fēng)險模型；趙佳虹等[9]基于環(huán)境風(fēng)險控制建立了環(huán)境風(fēng)險和成本最小化的選址-路徑模型；陳泊梨等[10]構(gòu)建了不確定條件下的多式聯(lián)運路徑優(yōu)化模型；唐燕等[11]設(shè)計了響應(yīng)動態(tài)環(huán)境溫度的大規(guī)模鄰域搜索車輛調(diào)度模型；Saeidi等[12]考慮多方?jīng)Q策利益，設(shè)計了雙層規(guī)劃模型；Shi等[13]提出了成本最小的醫(yī)療廢物回收模型；Nikzamir等[14]考慮隨機傳遞時間和車輛廢氣排放，構(gòu)建了醫(yī)療廢物選址-路徑模型；Tirkolaee等[15]基于時間窗約束構(gòu)建了醫(yī)療廢物選址-路徑模型；Kargar等[16]考慮了醫(yī)療廢物生產(chǎn)源點的差異性和多樣性，建立了成本、風(fēng)險和未收集廢物量最小化的逆向物流系統(tǒng)優(yōu)化模型。

綜上所述，以上研究成果在疫情下感染性醫(yī)療廢物運輸管理應(yīng)用中存在如下不足：(1)風(fēng)險度量模型未考慮病毒的環(huán)境傳播特性和風(fēng)險源的結(jié)構(gòu)特征；(2)傳統(tǒng)選址-路徑模型忽略了路段遞增運量對優(yōu)化目標(biāo)的直接影響；(3)多目標(biāo)組合優(yōu)化問題的求解缺少穩(wěn)定性和敏感性分析。因此，本研究提出了一類考慮風(fēng)險擴散的多目標(biāo)選址-路徑優(yōu)化模型，根據(jù)病毒的環(huán)境傳播特性構(gòu)建立體式風(fēng)險度量模型。引入路段運量遞增約束，構(gòu)建總成本和總風(fēng)險最小的選址-路徑模型，并基于最小包絡(luò)聚類分析法和NSGA-Ⅱ算法設(shè)計求解方法。最后以武漢實例和測試來驗證模型及算法的有效性。

1 立體式風(fēng)險度量

1.1 風(fēng)險定義

感染性醫(yī)療廢物的運輸風(fēng)險源自病毒泄露，風(fēng)險影響對象是周邊居民，居民吸入的病毒數(shù)量反映了風(fēng)險危害程度。本研究重點考慮病毒的環(huán)境傳播特性[17]，整合環(huán)境擴散模型和暴露人口模型，將感染性醫(yī)療廢物的立體式運輸風(fēng)險定義為運輸發(fā)生泄漏事故后，在應(yīng)急救援時間范圍內(nèi)暴露于有毒環(huán)境中的人口數(shù)量與吸入有毒物質(zhì)濃度的乘積。

1.2 風(fēng)險度量模型

如圖1所示，將運輸感染性醫(yī)療廢物的車輛設(shè)為風(fēng)險源，考慮風(fēng)速影響，風(fēng)險擴散過程可描述為：風(fēng)險以立體車廂為初始狀態(tài)，基于水平地面以一定擴散半徑進(jìn)行立體式傳播，病毒在有限擴散時間內(nèi)逐漸形成一個半圓柱體，則擴散半徑為：

圖1 運輸過程中風(fēng)險擴散示意圖Fig.1 Schematic diagram of risk diffusion during transport

R=u×SPT，

(1)

式中，R為擴散半徑；u為風(fēng)速；SPT為擴散時間。

擴散時間可表示為：

(2)

式中，veh為救援車輛平均行駛速度；Qij為感染性醫(yī)療廢物運量；DDij為應(yīng)急中心到弧(i，j)的距離；F為應(yīng)急中心救援能力。若H為車輛高度，則擴散角度值α計算為：

(3)

半圓柱體橫截面積S表示為：

(4)

設(shè)定Lij為弧(i，j)的距離，圖1的不規(guī)則圓柱體體積為：

(5)

設(shè)定θ為廢物的陽性率，ρij為弧(i，j)上的人口密度，暴露人口數(shù)量Pij的計算公式為：

Pij=πu2SPT2ρijθLij。

(6)

結(jié)合式(5)～(6)，構(gòu)建運輸風(fēng)險Bij的度量模型為：

(7)

2 數(shù)學(xué)模型

2.1 問題描述

如圖2所示，感染性醫(yī)療廢物運輸系統(tǒng)架構(gòu)于含有應(yīng)急系統(tǒng)的城市路網(wǎng)之上，包含了生產(chǎn)源點、回收站和臨時回收站3類節(jié)點。運輸車輛由一個回收站或臨時回收站出發(fā)，依次訪問各生產(chǎn)源點，然后返回起始的回收站或臨時回收站。若運輸時突發(fā)事故，應(yīng)急救援會迅速啟動，控制感染風(fēng)險的立體式擴散。因此，本研究的感染性醫(yī)療廢物運輸選址-路徑多目標(biāo)優(yōu)化問題是在應(yīng)急時間內(nèi)協(xié)同優(yōu)化多個回收站選址和多車輛路徑設(shè)計決策。

圖2 感染性醫(yī)療廢物運輸系統(tǒng)Fig.2 Infectious medical wastes transport system

2.2 基本假設(shè)

在構(gòu)建模型之前，設(shè)定如下假設(shè)條件：(1)不考慮環(huán)境因素動態(tài)變化；(2)設(shè)施和車輛都滿足處理感染性醫(yī)療廢物的安全要求；(3)運輸網(wǎng)絡(luò)沒有容量限制；(4)泄漏事故發(fā)生時，外界環(huán)境總處于有風(fēng)的狀態(tài)，導(dǎo)致病毒在空氣中擴散。

2.3 參數(shù)說明

根據(jù)設(shè)定的假設(shè)條件和建模背景，在最小化成本和風(fēng)險的基礎(chǔ)上構(gòu)建感染性醫(yī)療廢物運輸?shù)倪x址-路徑優(yōu)化模型，設(shè)置模型的集合如下：N(V，E)為運輸網(wǎng)絡(luò)，V=G∪T∪S為網(wǎng)絡(luò)節(jié)點，E為網(wǎng)絡(luò)弧。G∈{1，…，g}為生產(chǎn)源點集合；T∈{1，…，t}為臨時回收站建設(shè)候選點集合；S∈{1，…，s}為回收站集合；K∈{1，…，k}為運輸車輛集合。

2.4 多目標(biāo)選址-路徑模型

本研究以總成本和總風(fēng)險最小為優(yōu)化目標(biāo)，建立多目標(biāo)選址-路徑模型，下式為2個目標(biāo)函數(shù)：

(8)

(9)

式中，f1為最小化回收站的固定和變動成本、車輛購置和運輸成本；f2為總運輸風(fēng)險最小；oi為0～1變量，若建設(shè)回收站i∈T∪S，則為1，反之為0；TFCi為回收站i∈T∪S的固定成本；qijtk為連續(xù)變量，代表車輛k∈K在弧(i，j)∈E上運量，且運輸終點是回收站t∈T∪S；C為單位運輸成本；τk為0～1變量，若使用車輛k∈K，則為1，反之為0；VFCk為車輛k∈K的運營成本；ρρij為弧(i，j)∈E附近的人口密度；Hk為車輛k∈K的高度。

式(10)是唯一性約束，表示每個生產(chǎn)源點只屬于一條運輸路線，且這條運輸路線只有一輛車通過，而該車最終返回一個回收站。同時，每個生產(chǎn)源點只能由一個回收站服務(wù)，只能被一輛車訪問。式(11)表示每條運輸路線至多只有一輛車通過，且該車輛至多只屬于一個回收站。式(12)確保每條運輸路線的車輛從某一個回收站出發(fā)，且最終回到同一個回收站。式(13)表示運輸路線上任意兩節(jié)點的連接性。

(10)

(11)

?j∈G，?k∈K，

(12)

?k∈K，

(13)

式中，xijtk為0-1變量，若車輛k∈K經(jīng)過弧(i，j)∈E，且運輸終點為回收站t∈T∪S，則為1，反之為0；yitk為0-1變量，若點i∈G的廢物被車輛k∈K收集且運往回收站t∈T∪S，則為1，反之為0；ztk為0-1變量，若車輛k∈K服務(wù)回收站t∈T∪S，則為1，反之為0。

式(14)為流量守恒約束，表示該路段流量與上一路段流量之差等于連接這2個路段的節(jié)點廢物產(chǎn)量：

?t∈T∪S，?k∈K，

(14)

式中g(shù)i為源點i∈G的廢物產(chǎn)量。

式(15)為車輛的最大載重約束，式(16)為回收站的最大能力約束：

?i，j∈V，?t∈T∪S，?k∈K，

(15)

(16)

式中，VCAPk為車輛k∈K的最大載重量；CAPi為回收站i∈T∪S的最大回收能力。

式(17)和(18)為支路消除約束：

wjtk≥witk+1-(1-xijtk)M，?i，j∈G，

?t∈T∪S，?k∈K，

(17)

oi+oj+xijtk≤2，?i，j∈T，?t∈T∪S，

?k∈K，

(18)

式中，witk為整數(shù)決策變量，表示點i∈V在回收站t∈T∪S派出車輛k∈K的訪問次序；M為一個無窮大的正整數(shù)。

式(19)為決策變量的邏輯約束：

(19)

3 算法設(shè)計

本研究構(gòu)建了一個雙目標(biāo)0-1混合整數(shù)非線性模型，采用分階段求解思路，將選址-路徑組合優(yōu)化問題分解為選址-分配和路徑優(yōu)化問題，基于最小包絡(luò)聚類分析法和NSGA-II算法設(shè)計求解步驟[18-19]。

3.1 基于最小包絡(luò)聚類分析法的選址-分配算法

將相異度[20]設(shè)定為特征空間的距離函數(shù)，設(shè)定x*和y*的第i個變量值為xi*和yi*，變量數(shù)為q*，設(shè)計距離函數(shù)為：

(20)

3.1.1 帶價值屬性的歐式距離

考慮各生產(chǎn)源點的產(chǎn)量和單位運價及回收站固定成本的影響，設(shè)計單位歐式距離價值為C*，并構(gòu)建帶有價值屬性的歐式距離計算公式：

(21)

3.1.2 求解步驟

首先，引入如下定義：

(1)將還未分配給任何1個回收站t∈T∪S的生產(chǎn)源點i∈G構(gòu)成待分配集合A。

其次，設(shè)計具體求解步驟如下：

(2)聚類分析。對集合A的生產(chǎn)源點進(jìn)行聚類分析，使用最小包絡(luò)法確定回收站t∈T的選址位置，建立源點i與回收站t的分配關(guān)系。

3.2 基于NSGA-II的路徑優(yōu)化算法

3.2.1 步驟與流程

借鑒多目標(biāo)優(yōu)化算法[21]的求解思路，設(shè)計基于NSGA-Ⅱ的求解步驟如下：

(1)染色體編碼?；蛐痛硎具\輸路線，生產(chǎn)源點被某條路線服務(wù)，則基因為1，否則為0。

(2)隨機生成初始種群。初始種群集合N，設(shè)置進(jìn)化代數(shù)t=0，初始種群記為P(0)，第t代種群記為P(t)。

(3)適應(yīng)度計算。以個體非支配排序以及擁擠度綜合計算。

快速非支配排序：種群中的每個個體都有支配它的個體數(shù)量n(i)和被它支配的個體數(shù)量S(i)。當(dāng)n(i)=0時，則表示其沒有被其他任何個體所支配，將其存入第1個非支配層F(1)。將之前被它所支配的個體集合K中的所有n(K)減去1，解除被第1層支配的數(shù)量。以此類推，直至所有個體都被分級。

擁擠度計算：在同一等級的非支配排序下，計算某點在空間中與相鄰2點在多個目標(biāo)函數(shù)值下的相對距離之和。擁擠度能保證算法收斂于一個比較均勻的Pareto域，達(dá)到維護(hù)種群多樣性的目的。

(4)選擇。采用競標(biāo)賽選擇法隨機選擇2個個體，比較個體的非支配排序和擁擠度算子來評價其優(yōu)劣性。其中，非支配排序越低，越接近Pareto前沿，其性狀更優(yōu)；在同一排序級別下，擁擠度算子越小，其勝出概率越大。

(5)交叉。隨機交換2個染色體之間的信息，產(chǎn)生新個體。

(6)變異。通過變異算子作用于種群，以一定概率隨機改變某個體的一個基因，產(chǎn)生新個體。

(7)終止條件。設(shè)定最大遺傳代數(shù)MaxGen，當(dāng)t=MaxGen時，迭代結(jié)束。

3.2.2 Pareto最優(yōu)解方案選擇策略

首先引出如下定義：

(1)Pareto最優(yōu)解集構(gòu)成集合U*，其中方案i∈U*的成本和風(fēng)險分別記為cost(i)和risk(i)；

其次，方案選擇策略如下：

4 算例

4.1 基本信息

以武漢市新冠疫情下醫(yī)療廢物運輸管理為背景，設(shè)置網(wǎng)絡(luò)節(jié)點36個，其中醫(yī)療廢物生產(chǎn)源點為各大醫(yī)院，共有27個(G1～G27)，各醫(yī)院病床數(shù)占用量為100～1 200張，每張病床的醫(yī)療廢物產(chǎn)量為2.5～2.85 kg/d[22]。已有回收站設(shè)置2個(S1和S2)，臨時回收站候選點6個(T1～T6)，回收站最大能力均為12 t/d，回收站信息如表1所示。設(shè)定廢物陽性率為76.5%，平均風(fēng)速為10.8 km/h，運輸車輛最大載重量為4 t，高度為2.5 m，購置成本為15萬元/輛，單位運輸成本為30 元/km/t，應(yīng)急中心救援能力為4 t/h，所有車輛平均行駛速度為80 km/h。

表1 回收站信息Tab.1 Recycling stations

4.2 計算結(jié)果

初始種群數(shù)100個，交叉和變異概率為0.8和0.2，最大迭代次數(shù)為500次，在Intel(R)/CPU2.2 GHz/2.5 GHz 環(huán)境下，采用JAVA編程調(diào)用IBM ILOG CPLEX12.6.3求解選址-分配子問題，使用Matlab R2015a編程求解路徑優(yōu)化子問題。

如表2所示，經(jīng)過0.61 s求得總成本為6.01×106元的選址-分配方案，將該方案引入路徑優(yōu)化問題求解，經(jīng)過208.17 s共求得優(yōu)化方案40 320個。

表2 選址-分配方案Tab.2 Schemes of location-allocation

根據(jù)最優(yōu)方案選擇策略，推薦總運輸成本為0.83×104元、總運輸風(fēng)險為4.43×104(t·人)/km3的路徑優(yōu)化方案。如表3所示，推薦的選址-路徑方案共設(shè)立4個回收站和15輛車?？梢?，新模型和算法能夠求解含有778 527 個約束條件和400 833個決策變量的優(yōu)化問題，并在208.78 s內(nèi)求得優(yōu)化方案。

表3 選址-路徑方案Tab.3 Schemes of location-routing

4.3 結(jié)果分析

如表4所示，相較于武漢市的先行方案，新的優(yōu)化方案具有如下優(yōu)點：(1)降低醫(yī)療廢物運輸網(wǎng)絡(luò)的工作負(fù)荷，疫情期間，現(xiàn)有回收站都超負(fù)荷運轉(zhuǎn)，新方案增加2個臨時回收站，擴大總體回收能力，回收站負(fù)荷降低50%；(2)減少運輸車輛，新方案采用環(huán)式回收路線，減少車輛使用量44.44%；(3)減少運輸成本和風(fēng)險，新方案縮短了車輛行駛距離，運輸成本和風(fēng)險分別降低了57.87%和75.46%。

表4 方案對比結(jié)果Tab.4 Comparison of schemes

4.4 測試對比

為進(jìn)一步驗證模型和算法的有效性，基于武漢實例的運輸網(wǎng)絡(luò)，首先驗證不同風(fēng)險度量模型和計算方法的有效性。然后調(diào)整參數(shù)取值，測試模型敏感性。最后驗證不同計算規(guī)模下新算法的求解穩(wěn)定性。

4.4.1 風(fēng)險度量對比

分別以人口暴露噸數(shù)和新建的立體式風(fēng)險度量模型為風(fēng)險目標(biāo)函數(shù)，求解風(fēng)險最小的單目標(biāo)問題。如表5所示，相較于傳統(tǒng)模型，新建的風(fēng)險度量模型僅需增加3.75%的求解時間就可以降低10.08%的運輸成本。

表5 風(fēng)險度量模型對比結(jié)果Tab.5 Comparison of risk measurement models

4.4.2 求解方法對比

對比分析常規(guī)的線性加權(quán)多目標(biāo)優(yōu)化方法和新算法的求解效果。其中，線性加權(quán)的求解過程采用CPLEX編程計算，設(shè)定治療新冠確診患者的人均費用為1.7萬元來等價轉(zhuǎn)換風(fēng)險值，2個子目標(biāo)的權(quán)重系數(shù)均為0.5。如表6所示，相較于常規(guī)的多目標(biāo)求解方法，新算法能夠縮短79.53%的求解時間，使得優(yōu)化方案的最優(yōu)值平均差異率降低了2.93%。

表6 求解方法對比結(jié)果Tab.6 Comparison of different solution methods

4.4.3 參數(shù)敏感性

以武漢實例為基礎(chǔ)算例，根據(jù)以下情景進(jìn)行敏感性分析。情景1：車輛最大載重量增至5 t；情景2：車輛平均行駛速度降為40 km/h；情景3：廢物陽性率降至30%；情景4：平均風(fēng)速降為20 km/h。模型和算法敏感性分析結(jié)果如表7所示，增加車輛最大載重量，總風(fēng)險和成本分別降低了20%和1.82%；提高車輛平均行駛速度，總風(fēng)險降低了7.67%；降低陽性率，總風(fēng)險降低了67.55%；提升風(fēng)速，總風(fēng)險降低了36.27%。

表7 各情景下的計算結(jié)果對比Tab.7 Comparison of computation results in different scenarios

4.4.4 計算穩(wěn)定性

為驗證新算法的求解穩(wěn)定性，隨機生成4種不同規(guī)模的測試算例。除各類設(shè)施數(shù)量差異外，算例其他參數(shù)都相同。如表8所示，本研究設(shè)計的求解算法具有較強的計算穩(wěn)定性，能在307 s內(nèi)求得不同計算規(guī)模問題的最優(yōu)解。

表8 計算規(guī)模測試結(jié)果Tab.8 Testing result of computation scales

5 結(jié)論

為減少疫情下感染性醫(yī)療廢物的風(fēng)險，提出了多目標(biāo)選址-路徑模型。首先，根據(jù)病毒的環(huán)境傳播特性構(gòu)建了立體式的風(fēng)險度量模型；其次，設(shè)計了運量遞增約束，構(gòu)建了總成本和總風(fēng)險最小化的選址-路徑模型；然后，根據(jù)模型的復(fù)雜度設(shè)計了兩階段求解方法；最后，通過武漢實例和多組測試算例驗證了模型和算法的有效性。本研究結(jié)果可為相關(guān)管理部門制定疫情下的醫(yī)療廢物設(shè)施選址和路徑優(yōu)化提供決策支持。根據(jù)計算結(jié)果，本研究主要結(jié)論如下：

(1)新模型和算法能夠在209 s內(nèi)求得多個有效方案，且具有一定的參數(shù)敏感性。

(2)相較于現(xiàn)行方案，新方案可以分別降低50%的回收站工作負(fù)荷、57.87%的運輸成本和75.46%的運輸風(fēng)險。

(3)相較于傳統(tǒng)風(fēng)險模型，新建的風(fēng)險度量模型能夠降低10.08%的運輸成本。

(4)相較于常規(guī)多目標(biāo)求解方法，新算法能縮短79.53%的求解時間，并能在307 s內(nèi)求解不同計算規(guī)模的優(yōu)化問題。