0 引言

2020年新冠疫情爆发，感染性医疗废物产量激增，城市医疗废物运输安全面临极大挑战^[1]。合理评估运输风险、优化运输的设施选址和车辆路径决策是疫情下感染性医疗废物运输急待解决的安全问题。

医疗废物属于危险废物范畴，近年来国内外学者在危险废物运输风险度量和选址-路径优化问题等领域进行了大量研究^[2]。Yilmaz等^[3]设计了考虑人口和脆弱环境区域影响的风险模型；Taslimi等^[4]以最小化最大区域风险来实现风险公平性；Zhao等^[5]设计了不同水环境中的重油扩散风险模型；Yu等^[6]设计了人口感染风险模型；赵佳虹等^[7]考虑新冠病毒的气溶胶传播特性，构建了时变“环境-人口”风险模型；Taslimi等^[8]建立了医疗废物运输-存放风险模型；赵佳虹等^[9]基于环境风险控制建立了环境风险和成本最小化的选址-路径模型；陈泊梨等^[10]构建了不确定条件下的多式联运路径优化模型；唐燕等^[11]设计了响应动态环境温度的大规模邻域搜索车辆调度模型；Saeidi等^[12]考虑多方决策利益，设计了双层规划模型；Shi等^[13]提出了成本最小的医疗废物回收模型；Nikzamir等^[14]考虑随机传递时间和车辆废气排放，构建了医疗废物选址-路径模型；Tirkolaee等^[15]基于时间窗约束构建了医疗废物选址-路径模型；Kargar等^[16]考虑了医疗废物生产源点的差异性和多样性，建立了成本、风险和未收集废物量最小化的逆向物流系统优化模型。

综上所述，以上研究成果在疫情下感染性医疗废物运输管理应用中存在如下不足：(1) 风险度量模型未考虑病毒的环境传播特性和风险源的结构特征；(2) 传统选址-路径模型忽略了路段递增运量对优化目标的直接影响；(3) 多目标组合优化问题的求解缺少稳定性和敏感性分析。因此，本研究提出了一类考虑风险扩散的多目标选址-路径优化模型，根据病毒的环境传播特性构建立体式风险度量模型。引入路段运量递增约束，构建总成本和总风险最小的选址-路径模型，并基于最小包络聚类分析法和NSGA-Ⅱ算法设计求解方法。最后以武汉实例和测试来验证模型及算法的有效性。

1 立体式风险度量 1.1 风险定义

感染性医疗废物的运输风险源自病毒泄露，风险影响对象是周边居民，居民吸入的病毒数量反映了风险危害程度。本研究重点考虑病毒的环境传播特性^[17]，整合环境扩散模型和暴露人口模型，将感染性医疗废物的立体式运输风险定义为运输发生泄漏事故后，在应急救援时间范围内暴露于有毒环境中的人口数量与吸入有毒物质浓度的乘积。

1.2 风险度量模型

如图 1所示，将运输感染性医疗废物的车辆设为风险源，考虑风速影响，风险扩散过程可描述为：风险以立体车厢为初始状态，基于水平地面以一定扩散半径进行立体式传播，病毒在有限扩散时间内逐渐形成一个半圆柱体，则扩散半径为：

(1)

式中，R为扩散半径；u为风速；SPT为扩散时间。

扩散时间可表示为：

(2)

式中，veh为救援车辆平均行驶速度；Q_ij为感染性医疗废物运量；DD_ij为应急中心到弧(i，j)的距离；F为应急中心救援能力。若H为车辆高度，则扩散角度值α计算为：

(3)

半圆柱体横截面积S表示为：

(4)

设定L_ij为弧(i, j)的距离，图 1的不规则圆柱体体积为：

(5)

设定θ为废物的阳性率，ρ_ij为弧(i, j)上的人口密度，暴露人口数量Pij的计算公式为：

(6)

结合式(5)~(6)，构建运输风险B_ij的度量模型为：

(7)

2 数学模型 2.1 问题描述

如图 2所示，感染性医疗废物运输系统架构于含有应急系统的城市路网之上，包含了生产源点、回收站和临时回收站3类节点。运输车辆由一个回收站或临时回收站出发，依次访问各生产源点，然后返回起始的回收站或临时回收站。若运输时突发事故，应急救援会迅速启动，控制感染风险的立体式扩散。因此，本研究的感染性医疗废物运输选址-路径多目标优化问题是在应急时间内协同优化多个回收站选址和多车辆路径设计决策。

2.2 基本假设

在构建模型之前，设定如下假设条件：(1)不考虑环境因素动态变化；(2)设施和车辆都满足处理感染性医疗废物的安全要求；(3)运输网络没有容量限制；(4)泄漏事故发生时，外界环境总处于有风的状态，导致病毒在空气中扩散。

2.3 参数说明

根据设定的假设条件和建模背景，在最小化成本和风险的基础上构建感染性医疗废物运输的选址-路径优化模型，设置模型的集合如下：N(V, E)为运输网络，V=G∪T∪S为网络节点，E为网络弧。G∈ {1, ⋯, g}为生产源点集合；T∈ {1, ⋯, t}为临时回收站建设候选点集合；S∈ {1, ⋯, s}为回收站集合；K∈ {1, ⋯, k}为运输车辆集合。

2.4 多目标选址-路径模型

本研究以总成本和总风险最小为优化目标，建立多目标选址-路径模型，下式为2个目标函数：

(8)

(9)

式中，f₁为最小化回收站的固定和变动成本、车辆购置和运输成本；f₂为总运输风险最小；o_i为0~1变量，若建设回收站i∈T∪S，则为1，反之为0；TFC_i为回收站i∈T∪S的固定成本；q_ijtk为连续变量，代表车辆k∈K在弧(i, j)∈E上运量，且运输终点是回收站t∈T∪S；C为单位运输成本；τ_k为0~1变量，若使用车辆k∈K，则为1，反之为0；VFC_k为车辆k∈K的运营成本；ρρ_ij为弧(i, j)∈E附近的人口密度；H_k为车辆k∈K的高度。

式(10)是唯一性约束，表示每个生产源点只属于一条运输路线，且这条运输路线只有一辆车通过，而该车最终返回一个回收站。同时，每个生产源点只能由一个回收站服务，只能被一辆车访问。式(11)表示每条运输路线至多只有一辆车通过，且该车辆至多只属于一个回收站。式(12)确保每条运输路线的车辆从某一个回收站出发，且最终回到同一个回收站。式(13)表示运输路线上任意两节点的连接性。

(10)

(11)

(12)

(13)

式中，x_ijtk为0-1变量，若车辆k∈K经过弧(i, j)∈E，且运输终点为回收站t∈T∪S，则为1，反之为0；y_itk为0-1变量，若点i∈G的废物被车辆k∈K收集且运往回收站t∈T∪S，则为1，反之为0；z_tk为0-1变量，若车辆k∈K服务回收站t∈T∪S，则为1，反之为0。

式(14)为流量守恒约束，表示该路段流量与上一路段流量之差等于连接这2个路段的节点废物产量：

(14)

式中g_i为源点i∈G的废物产量。

式(15)为车辆的最大载重约束，式(16)为回收站的最大能力约束：

(15)

(16)

式中，VCAP_k为车辆k∈K的最大载重量；CAP_i为回收站i∈T∪S的最大回收能力。

式(17)和(18)为支路消除约束：

(17)

(18)

式中，w_itk为整数决策变量，表示点i∈V在回收站t∈T∪S派出车辆k∈K的访问次序；M为一个无穷大的正整数。

式(19)为决策变量的逻辑约束：

(19)

3 算法设计

本研究构建了一个双目标0-1混合整数非线性模型，采用分阶段求解思路，将选址-路径组合优化问题分解为选址-分配和路径优化问题，基于最小包络聚类分析法和NSGA-Ⅱ算法设计求解步骤^[18-19]。

3.1 基于最小包络聚类分析法的选址-分配算法

将相异度^[20]设定为特征空间的距离函数，设定x^*和y^*的第i个变量值为xi^*和yi^*，变量数为q^*，设计距离函数为：

(20)

3.1.1 带价值属性的欧式距离

考虑各生产源点的产量和单位运价及回收站固定成本的影响，设计单位欧式距离价值为C^*，并构建带有价值属性的欧式距离计算公式：

(21)

3.1.2 求解步骤

首先，引入如下定义：

(1) 将还未分配给任何1个回收站t∈T∪S的生产源点i∈G构成待分配集合A。

(2) 若源点i∈G与回收站t∈T∪S解除分配关系，则将源点计入此回收站的禁入列表R_t，该源点无法再次与回收站建立新的分配关系。

其次，设计具体求解步骤如下：

(1) 初始化。回收站初始容量为0，清空回收站t∈T禁入列表R_t，所有源点放入待分配集合A；

(2) 聚类分析。对集合A的生产源点进行聚类分析，使用最小包络法确定回收站t∈T的选址位置，建立源点i与回收站t的分配关系。

(3) 解除分配关系。当所有生产源点的废物总量超过回收站的最大能力时，则按源点与回收站的欧式距离降序顺序依次解除源点与回收站的分配关系，直至满足回收站最大能力为止。同时，将已解除分配关系的源点重新存入待分配集合A和相应的禁入列表R_t。

3.2 基于NSGA-Ⅱ的路径优化算法 3.2.1 步骤与流程

借鉴多目标优化算法^[21]的求解思路，设计基于NSGA-Ⅱ的求解步骤如下：

(1) 染色体编码。基因型串表示运输路线，生产源点被某条路线服务，则基因为1，否则为0。

(2) 随机生成初始种群。初始种群集合N，设置进化代数t=0，初始种群记为P (0)，第t代种群记为P(t)。

(3) 适应度计算。以个体非支配排序以及拥挤度综合计算。

快速非支配排序：种群中的每个个体都有支配它的个体数量n(i)和被它支配的个体数量S(i)。当n(i)=0时，则表示其没有被其他任何个体所支配，将其存入第1个非支配层F(1)。将之前被它所支配的个体集合K中的所有n(K)减去1，解除被第1层支配的数量。以此类推，直至所有个体都被分级。

拥挤度计算：在同一等级的非支配排序下，计算某点在空间中与相邻2点在多个目标函数值下的相对距离之和。拥挤度能保证算法收敛于一个比较均匀的Pareto域，达到维护种群多样性的目的。

(4) 选择。采用竞标赛选择法随机选择2个个体，比较个体的非支配排序和拥挤度算子来评价其优劣性。其中，非支配排序越低，越接近Pareto前沿，其性状更优；在同一排序级别下，拥挤度算子越小，其胜出概率越大。

(5) 交叉。随机交换2个染色体之间的信息，产生新个体。

(6) 变异。通过变异算子作用于种群，以一定概率随机改变某个体的一个基因，产生新个体。

(7) 终止条件。设定最大遗传代数MaxGen，当t=MaxGen时，迭代结束。

3.2.2 Pareto最优解方案选择策略

首先引出如下定义：

(1) Pareto最优解集构成集合U^*，其中方案i∈U^*的成本和风险分别记为cost(i)和risk(i)；

(2) 设所有车辆运输路线方案的成本和风险的等权平均值分别为cost和risk；

其次，方案选择策略如下：

(1) 若存在唯一方案i∈U^*，使得|risk(i)-risk|最小，则选择方案i；

(2) 若存在方案i, j∈U^*，使得|risk(i)-risk|=|risk(j)-risk|且最小，|cost(i)-cost| < |cost(j)-cost|，则选择方案i；

(3) 若存在方案i, j∈U^*，使得|risk(i)-risk|=|risk(j)-risk|且最小，|cost(i)-cost|=|cost(j)-cost|，则当risk(i) < risk(j)时，选择方案i。

4 算例 4.1 基本信息

以武汉市新冠疫情下医疗废物运输管理为背景，设置网络节点36个，其中医疗废物生产源点为各大医院，共有27个(G1~G27)，各医院病床数占用量为100~1 200张，每张病床的医疗废物产量为2.5~2.85 kg/d^[22]。已有回收站设置2个(S1和S2)，临时回收站候选点6个(T1~T6)，回收站最大能力均为12 t/d，回收站信息如表 1所示。设定废物阳性率为76.5%，平均风速为10.8 km/h，运输车辆最大载重量为4 t，高度为2.5 m，购置成本为15万元/辆，单位运输成本为30元/km/t，应急中心救援能力为4 t/h，所有车辆平均行驶速度为80 km/h。

4.2 计算结果

初始种群数100个，交叉和变异概率为0.8和0.2，最大迭代次数为500次，在Intel (R)/CPU2.2 GHz/2.5 GHz环境下，采用JAVA编程调用IBM ILOG CPLEX12.6.3求解选址-分配子问题，使用Matlab R2015a编程求解路径优化子问题。

如表 2所示，经过0.61 s求得总成本为6.01×10⁶元的选址-分配方案，将该方案引入路径优化问题求解，经过208.17 s共求得优化方案40 320个。

根据最优方案选择策略，推荐总运输成本为0.83×10⁴元、总运输风险为4.43×10⁴ (t·人)/km³的路径优化方案。如表 3所示，推荐的选址-路径方案共设立4个回收站和15辆车。可见，新模型和算法能够求解含有778 527个约束条件和400 833个决策变量的优化问题，并在208.78 s内求得优化方案。

4.3 结果分析

如表 4所示，相较于武汉市的先行方案，新的优化方案具有如下优点：(1) 降低医疗废物运输网络的工作负荷，疫情期间，现有回收站都超负荷运转，新方案增加2个临时回收站，扩大总体回收能力，回收站负荷降低50%；(2) 减少运输车辆，新方案采用环式回收路线，减少车辆使用量44.44%；(3) 减少运输成本和风险，新方案缩短了车辆行驶距离，运输成本和风险分别降低了57.87%和75.46%。

4.4 测试对比

为进一步验证模型和算法的有效性，基于武汉实例的运输网络，首先验证不同风险度量模型和计算方法的有效性。然后调整参数取值，测试模型敏感性。最后验证不同计算规模下新算法的求解稳定性。

4.4.1 风险度量对比

分别以人口暴露吨数和新建的立体式风险度量模型为风险目标函数，求解风险最小的单目标问题。如表 5所示，相较于传统模型，新建的风险度量模型仅需增加3.75%的求解时间就可以降低10.08%的运输成本。

4.4.2 求解方法对比

对比分析常规的线性加权多目标优化方法和新算法的求解效果。其中，线性加权的求解过程采用CPLEX编程计算，设定治疗新冠确诊患者的人均费用为1.7万元来等价转换风险值，2个子目标的权重系数均为0.5。如表 6所示，相较于常规的多目标求解方法，新算法能够缩短79.53%的求解时间，使得优化方案的最优值平均差异率降低了2.93%。

4.4.3 参数敏感性

以武汉实例为基础算例，根据以下情景进行敏感性分析。情景1：车辆最大载重量增至5 t；情景2：车辆平均行驶速度降为40 km/h；情景3：废物阳性率降至30%；情景4：平均风速降为20 km/h。模型和算法敏感性分析结果如表 7所示，增加车辆最大载重量，总风险和成本分别降低了20%和1.82%；提高车辆平均行驶速度，总风险降低了7.67%；降低阳性率，总风险降低了67.55%；提升风速，总风险降低了36.27%。

4.4.4 计算稳定性

为验证新算法的求解稳定性，随机生成4种不同规模的测试算例。除各类设施数量差异外，算例其他参数都相同。如表 8所示，本研究设计的求解算法具有较强的计算稳定性，能在307 s内求得不同计算规模问题的最优解。

5 结论

为减少疫情下感染性医疗废物的风险，提出了多目标选址-路径模型。首先，根据病毒的环境传播特性构建了立体式的风险度量模型；其次，设计了运量递增约束，构建了总成本和总风险最小化的选址-路径模型；然后，根据模型的复杂度设计了两阶段求解方法；最后，通过武汉实例和多组测试算例验证了模型和算法的有效性。本研究结果可为相关管理部门制定疫情下的医疗废物设施选址和路径优化提供决策支持。根据计算结果，本研究主要结论如下：

(1) 新模型和算法能够在209 s内求得多个有效方案，且具有一定的参数敏感性。

(2) 相较于现行方案，新方案可以分别降低50%的回收站工作负荷、57.87%的运输成本和75.46%的运输风险。

(3) 相较于传统风险模型，新建的风险度量模型能够降低10.08%的运输成本。

(4) 相较于常规多目标求解方法，新算法能缩短79.53%的求解时间，并能在307 s内求解不同计算规模的优化问题。