0 引言

近年来，中国持续推行“清废行动”和“无废城市”建设，城市危险废物运输网络优化受到了社会公众的普遍关注。城市危险废物大多来源于工业化建设和医疗健康活动，具有毒性、易燃易爆性或感染性^[1]。在运输网络中优化设施建设位置，合理配置运输车辆并设计运输路线，是降低城市危险废物运输的公众感知风险、提升网络运营能力的关键问题。

国内外的城市危险废物运输网络优化研究主要涉及运输风险度量、优化模型构建与求解方法设计3个领域^[2]。首先，在运输风险度量方面，Erkut等^[3]考虑了社会公众对事故概率和后果的认知差异，提出了一类感知风险模型；Ke等^[4]将模糊事件率、人口暴露数和暴露时间的乘积作为风险衡量指标来构建度量模型；Govindan等^[5]根据危险废物的感染性，将暴露的人口数量设定为风险度量的建模基础。其次，在优化建模与求解方面，赵佳虹等^[6]建立了双目标的数学模型，协同优化设施选址和车辆路径设计问题；黄丽霞等^[7]研究了区域型的危险废物运输网络规划问题，改进了多目标优化模型，并设计了基于增广加权切比雪夫的求解方法；Rabbani等^[8]研究了多周期多种类危险废物选址-路径-库存问题，并采用NSGA-Ⅱ算法设计求解步骤；由于疫情产生的医疗废物同样具有危险性，因此Aydemir^[9]建立了医疗废物的多周期选址-路径双目标模型，并设计了两阶段的多目标自适应大领域搜索算法。另外，由于危险废物属于特定类型的危险品，因此Hao等^[10]在考虑城际道路交通限制的约束下，研究了危险品运输的选址-路径优化建模与求解；Hassanpour等^[11]研究了具有随机中断路段、时变参数和配送时间窗的危险物品选址-路径问题，并根据增广ε-约束方法设计了求解步骤; Jiao等^[12]将基于分解的多目标进化算法(MOEA/D)运用于求解高度约束的多目标优化问题，验证了算法求解多组基准优化问题的有效性；Wang等^[13]也采用了MOEA/D求解电动汽车复合制动系统的能量消耗优化问题，并取得了更好的优化效果。

上述研究成果能够为城市危险废物运输网络的优化建模与求解提供借鉴，但是仍然存在不足：(1)没有考虑城市系统对危险废物运输网络设计的影响。(2)忽视了社会公众的风险感知差异。(3)忽略了运输车辆配置和路线选择决策的协同优化。(4)求解组合优化问题时，建模和算法设计有待进一步提升。(5)对模型和方法的求解稳定性和敏感性分析不足。因此，本研究充分考虑城市系统对危险废物运输网络的影响，基于公众风险感知差异，提出危险废物运输网络多目标优化建模与求解方法，解决多品类危险废物的运输设施选址、异构型运输车辆配置及运输网络设计的组合优化问题。考虑应急响应时间对公众感知的影响，设计感知风险度量模型。采用双商品流建模方法，构建总成本和总感知风险最小的运输网络优化模型，并基于MOEA/D-分支切割设计三阶段求解算法。最后，以昆明市实例和多个测试分析模型和算法的有效性、敏感性和稳定性。

1 感知风险度量 1.1 风险定义

在城市路网中运输危险废物时，若突发泄露事故，周边的社会公众将受到有害影响。因此，基于暴露人口吨数模型^[1]，综合考虑社会公众的风险感知差异^[14]，将危险废物运输网络的风险定义为：在城市路网中突发危险废物泄漏事故时，一定范围内的社会公众所感知的风险量。其具体表示为危险废物泄露吨数与带有风险感知系数的暴露人口数量的乘积，其中，风险感知系数与城市路网的应急响应时间密切相关。

1.2 风险度量模型

城市路网中设有固定的应急设施，面对突发事故时，应急系统迅速响应，应急响应时间的长短，对已有风险的消除和潜在风险的预防，具有至关重要的作用，这直接影响着社会公众的风险感知态度。因此，本研究以应急响应时间阈值RTTV为标准^[10]，RT为城市路网的实际应急响应时间，则风险感知系数PC为：

(1)

当实际应急时间小于阈值时(PC＜1)，应急效果最好，此时社会公众的风险感知态度为乐观；当实际应急时间等于阈值时(PC=1)，社会公众的风险感知为中观；当实际应急时间大于阈值时(PC>1)，社会公众的风险感知态度为悲观。节点和路段的感知风险PR_i^node和PR_i^edge为：

(2)

(3)

式中，Q_i^node和Q_ij^edge分别为节点i和路段(i, j)泄露的危险废物量；POP_i^node和POP_ij^edge分别为节点和路段的人口密度；PC_i^node和PC_ij^edge分别为节点和路段发生事故时公众的风险感知系数。

2 数学模型 2.1 问题描述

城市危险废物运输网络的构建基础是既有的城市应急系统，如图 1所示。该网络包含了生产源点、储存中心、集成处理中心和应急中心4类节点。不同类型的危险废物从生产源点收集至储存中心进行短期存储，然后转运至集成处理中心进行集中处置。在废物收集阶段，不同类型的车辆从任意一个储存中心出发，依次访问多个生产源点收集对应类型的危险废物，满载后回到原储存中心，收集路线的起点和终点是同一个储存中心。在废物转运阶段，车辆在储存中心和集成处理中心之间进行往复式的运输。在以上2个阶段，若突发泄露事故，城市既有的应急求援中心会迅速启动应急响应。本研究基于既有的城市应急体统，以总成本和总感知风险最小化为目标，协同优化危险品运输网络的多类设施选址、异构型车辆配置和运输路线设计问题。其中，为高效协调优化收集路线和转运路线，本研究采用双商品流建模方法设计运输车辆路径优化内容。

2.2 基本假设

为构建城市危险废物运输网络优化模型，设定基本假设条件为：(1)运输网络的外部环境相对稳定，各路段是双向道路，不考虑道路运输容量限制。(2)各类危险废物具有运输和存储相容性，即不相容的危险废物不能同时在同一车辆上运输，在存储中心需分区存储。(3)突发事故时，储存中心和车辆的危险废物量视为泄露量。

2.3 参数说明

参数说明如表 1所示。

2.4 多种类多车辆危险废物选址-路径模型

式(4)为总成本最小化，包含储存中心固定建设成本、废物处理和储存成本、车辆购置成本、废物运输和转运成本。式(5)为总感知风险的最小化，包括储存中心选址风险、废物收集和转运风险。本研究采用双商品流建模方法，每条收集路线包括2个方向的废物流量边，选择收集路线时，每1条流量边都参与计算。为避免重复统计，总成本最小化的目标函数如下，其中，收集路线的运输成本和感知风险的计算项都乘以1/2。

(4)

(5)

式(6)~(24)为约束条件。其中，式(6)为每个生产源点的净增流量等于该点产生的各类危险废物产量的2倍：

(6)

式(7)为确保复制储存中心的总流出量等于剩余的车辆容量：

(7)

式(8)为收集车辆的额定载重量约束：

(8)

式(9)规定了每个复制存储中心的危险废物总量等于危险废物的总产量：

(9)

式(10)规定了任意2点之间的收集路线含有2个方向的流量边：

(10)

式(11)为确保每条收集路线配备1个收集车辆：

(11)

式(12)为在收集路线的任意2个节点之间，2个相反方向的流量之和等于该收集车辆的额定载重量：

(12)

式(13)~(16)为决策变量逻辑约束：

(13)

(14)

(15)

(16)

式(17) 为每个储存中心的各类危险废物储存量等于流入其对应复制储存中心的危险废物总量：

(17)

式(18) 为各储存中心的危险废物储存量都将被转运至集成处理中心，以及储存中心能力约束：

(18)

式(19) 为每个生产源点的废物都会被1个储存中心使用特定类型的车辆收集：

(19)

式(20)~(21)为支路消除约束：

(20)

(21)

式(22)为收集车辆的数量限制：

(22)

3 算法设计

本研究构建的是一个双目标0-1混合整数线性网络规划模型，含有成本和风险2个不同量纲的目标函数。针对模型特征，本研究借鉴基于分解的多目标进化算法(Multi-objective Evolutionary Algorithm Based on Decomposition, MOEA/D)求解思路，设计了一种基于MOEA/D-分支切割三阶段混合多目标求解算法(Three-stage Hybrid Multi-objective Algorithm, THMOA)。MOEA/D具有较强的计算能力，能够采用聚合函数将多目标优化问题转化为多个子问题，以协作方式对多个子问题进行协同优化，算法设计流程如图 2所示。本研究在第1阶段采用MOEA/D，输出对应的选址-分配Pareto解集；在第2阶段输入上一阶段的Pareto解集，并沿用权重分配方案，通过分支-切割算法快速求解车辆路径子问题；最后，整合前2个阶段的求解方案，形成网络优化方案的Pareto解集，通过设计非支配解选择策略，选取最终推荐方案。

3.1 基于MOEA/D的求解步骤(选址-分配子问题)

本研究改进了MOEA/D，在分解策略上选择切比雪夫方法对作为聚合函数，并采用实数编码的交叉变异规则来进化种群。具体计算步骤如下。

步骤1：输入参数，初始化种群。该步骤主要包括随机生成规定数量的初始解，生成均匀权重组合，为个体分配规定数量的邻居，计算总成本和总风险目标函数值，使用切比雪夫法计算适应度值，筛选当前种群中的非支配解，建立精英种群。本研究采用实数向量对个体进行编码。假设当前求解的问题包含2类危险废物，共有10个节点，分别为2个集成处理中心，3个储存中心和5个生产源点，则构造的代表个体的向量大小为2行10列。每一行包含3个部分：第1部分为表示集成中心的开放状态；第2部分表示每个储存中心的选址及分配给的集成处理中心；第3部分表示生产源点-储存中心的分配关系。个体编码如图 3所示。

如图 3所示，每一行特定类型的危险废物对应集成处理中心、存储中心、生产源点编号3列内容。其中，“集成处理中心”列表示集成处理中心1和2都可以处理危险废物1；由“储存中心”列可知，储存中心3的危险废物1由处理中心2负责，储存中心5的危险废物1由处理中心1负责；根据“生产源点”列可知，生产源点6和10分配给储存中心3，生产源点7，8，9分配给储存中心5。

步骤2：交叉变异操作，生成新解。使用随机次序，依次在种群中的个体中随机抽取出邻居个体，进行交叉变异产生新解。比较交叉变异的子代解和父代解的适应度值ψ^*与ψ，若ψ^*≥ψ，则使用子代解替换原解。同时，抽取父代解的指定数量的邻居，比较在邻居个体的权重组合，计算其适应度值，若子代解的适应度值更小，则将邻居个体的解也替换为新的子代解。本研究采用部分映射交叉方法，通过交叉个体与邻居个体来更新解。为保障种群的多样性和扩展搜索空间，本研究采用互换和逆序2类交叉变异操作，如图 4所示，图中黑色底区域的基因位进行了设定的交叉和变异操作。互换变异是对个体内部2个位置互换数值，逆序变异是对个体中随机数量的一部分进行位置逆转。

步骤3：精英种群更新。循环所有个体，筛选出种群中新的非支配解，加入精英种群，再进行非支配解的筛选，重新组成精英种群，进入下一轮迭代。

步骤4：判断终止条件。若达到终止条件，则停止，输出结果。否则转至步骤2。

3.2 基于分支切割算法的求解步骤(车辆路径子问题)

导出第一阶段求得的选址-分配非支配解集，将集合中的每一个解代入车辆路径子问题的数学模型中，利用分支切割算法直接编程求解。具体步骤如下。

步骤1：按废物种类和储存中心的顺序，在每一个选址-分配非支配解中，依次找到储存中心分配的生产源点，并且按照储存中心-生产源点-复制储存中心的顺序排列编号，生成对应的距离矩阵、路段人口密度和风险感知系数等数据。

步骤2：输入数据，求解模型，并输出最优方案，获得成本和风险理想点。

步骤3：输入基础数据、理想点和该个体的权重方案，采用切比雪夫法聚合模型的2个优化子目标，求解转化后的单目标模型，得到每一个选址-分配非支配解对应的车辆路径最优方案，将相关结果记入个体中。

步骤4：重复步骤1~3，直至所有个体的车辆路径方案均求解完毕，再对选址分配和车辆路径的成本风险进行整合计算，得到所有个体的总成本和总风险。根据所有个体的总成本和总风险进行非支配排序，得到选址-路径的非支配解集。

3.3 基于灰色关联分析法的非支配解选择策略

本研究采用灰色关联分析法^[15]，以灰色相关系数描述每个非支配解与理想点之间的相似度。在规避决策者主观偏好的前提下，灰色关联分析法能够以客观数据为导向，设计了从非支配解集中推荐最优方案的选择策略。具体步骤如下。

步骤1，假设具有n个目标的优化问题，共求得m个非支配解，f_ij为第i个非支配解的第j个目标值，对每个非支配解的目标值都进行归一化处理后，得到目标值F_ij为：

(23)

步骤2，在归一化处理后的目标值F_ij中找到理想目标值F_j^*：

(24)

步骤3，求出每个目标值F_ij与理想目标值F_j^*之间的差异值∂_ij：

(25)

步骤4，求出第i个非支配解的灰色关联系数GRC_i：

(26)

式中，∂_max和∂_min分别为∂_ij (i∈m, j∈n)中的最大值和最小值。

步骤5，找到最大的GRC_i，对应的第i个非支配解就是推荐的最优方案。

4 算例 4.1 实例信息

以昆明市危险废物运输管理为背景，选取20个市内生产废酸和废矿物油的工业园区和企业，并设定为生产源点G1~G20，废酸和废矿物油产量的取值区间为0.49~4.50 t/月和0.41~5.40 t/月。依据文献[16-17]的参数设置，建立4个储存中心候选点S1~S4，建设废酸和废矿物油的固定建设成本，如表 2所示。各存储中心候选点的最大容量均为30 t，储存废酸和废矿物油的单位成本分别为500元/t和550元/t。设立1个集成处理中心T1，处理废酸和废矿物油的单位成本分别为6 800元/t和4 500元/t。废酸收集车辆2类，最大载重量为7 t和10 t，购置成本为12万元和18万元，单位收运成本为330元/km和380元/km。废矿物油收集车辆2类，最大载重量为10 t和12 t，购置成本为20万元和26万元，单位收运成本为450元/km和540元/km。2种废物的单位转运成本分别为30元/ (km·t)和35元/ (km·t)，其影响半径分别为0.8 km和0.5 km，最佳应急响应时间阈值为10 min^[18]。

4.2 计算结果与分析

本研究结合文献[19]的参数设置方法和测试，设置初始种群数=200，最大迭代数=300，精英种群大小=400，权重向量组数为200，邻域大小为30，交叉、变异概率分别为0.8和0.2，在AMD/Ryzen/75800H/CPU3.2GHz的计算环境下，采用MATLAB的R2018b版本和Gurobi的9.5.0版本软件进行编程计算。

MOEA/D经过124.47 s求得38组有非支配关系的优化方案。将所有方案代入第2阶段进行求解，用时1 053.33 s，通过非支配排序最终得到11个优化方案。根据设计的灰色关联分析法，推荐总成本为4.583×10⁶元、总风险为11.651×10⁶ t·人的优化方案。该网络优化方案的设施选址、收集车辆配置和收集路径规划的决策内容如表 3所示。该方案为不同的危险废物建立了2个储存中心，购置车型1和2的收集车辆共13辆。新模型和算法能够求解含有28 363个约束条件和14 471个决策变量的优化问题，并在1 177.8 s内求得推荐的优化方案。

将本研究在实际案例中求得的推荐方案与只考虑直达运输形式的原始方案进行对比分析，如表 4所示。相比于原始方案，新方案采用环式回收路线，仅需13辆车，减少车辆使用量为60%，极大缩短了车辆行驶距离，提高运转效率，分别降低了53.58%的运输成本和30.56%的运输风险。

4.3 测试对比

为进一步验证模型和算法有效性，基于昆明市实例的运输网络，分别验证风险度量模型、数学模型和计算方法的有效性，测试模型参数敏感性和算法的求解稳定性。

4.3.1 风险度量对比

分别以人口暴露风险模型和新建的感知风险度量模型为风险目标函数，求解风险最小的单目标问题。风险度量模型对比结果如表 5所示。相较于传统风险度量模型，新模型在增加11.0%的求解时间情况下，降低了10.2%的成本。

4.3.2 优化模型对比

采用求解软件Gurobi 9.5.0，限定求解时间为3 600 s，比较传统模型与本研究新建模型的求解性能，结果如表 6所示。相较于传统优化模型，新建模型分别减少了6.08%的决策变量数量和48.98%的约束条件数量，具备了更优的求解性能，能在限定时间内求得总成本和总风险更小的优化方案，其最优解的差距值优化了5.22%。

4.3.3 求解方法对比

采用Gurobi软件的版本9.5.0编程增广加权切比雪夫方法，并求解昆明市实例，将其计算结果与本研究新算法的计算结果进行对比分析。其中，增广加权切比雪夫方法的目标权重系数沿用新算法的设置参数。根据求解方法对比结果(见表 7)显示，新算法能缩短67.28%的求解时间，其最终方案的差距值优化了3.20%。

4.3.4 参数敏感性

以昆明市实例为基础算例，设定不同情景分析参数敏感性。情景1，收集车辆额定载重量增加50%；情景2，储存中心容量增加50%；情景3，应急响应时间阈值增加为15 min。根据不同情景下的计算结果(见表 8)可知：新模型和算法对收集车辆额定载重量、储存中心容量和应急响应时间阈值参数具有一定敏感性，其中，增大车辆载重能力，总成本降低10.26%，总风险提高28.15%，求解时间缩短9.47%；提高储存中心容量，总成本降低2.62%，总风险增加8.24%；增加应急响应时间阈值，总风险降低12.10%。

4.3.5 计算稳定性

随机生成3种不同规模的测试算例，测试新算法的求解稳定性。根据不同计算规模的测试结果(见表 9)可知，新算法能在2 772 s以内求解不同计算规模问题，并且在第1阶段和第2阶段的计算过程中，新算法能够求得稳定数量的非支配解。

5 结论

为降低危险废物的运输风险，提高网络运转效率，本研究考虑了公众感知差异，并提出了一类城市危险废物运输网络的多目标优化建模与求解方法，协同优化了多品类危险废物的运输设施选址、异构型车辆配置及运输网络设计问题。首先，考虑城市应急响应时间对公众感知的影响，设计了感知风险度量模型；其次，采用双商品流建模方法，构建了总成本和总风险最小的运输网络优化模型；然后，针对模型计算复杂度，设计了基于MOEA/D-分支切割的三阶段多目标求解算法。最后，结合昆明市实例和测试验证了模型和算法有效性。

(1) 新模型和算法能够在1 177.8 s内求得11个有效方案；新模型和算法对试验参数具有一定敏感性。

(2) 相较于直达往复式的运输模式，考虑环式路径的最优方案可分别降低60%的车辆数量，53.58%的运输成本，30.56%的运输风险。

(3) 相较于传统风险模型，使用新风险度量模型求得的方案能降低运输成本10.20%。

(4) 相较于传统优化模型，在相同的求解时间内，新的运输网络优化模型可分别降低3.4%的总成本和8.78%的总风险，优化5.22%的差距值，并减少6.08%的决策变量和48.98%的约束条件数量。

(5) 相较于增广加权切比雪夫的常规多目标方案，新设计的算法能够缩短67.28%的求解时间，并求得差距值更优的方案，同时具有较高的计算稳定性，能在2 772 s以内求解不同计算规模的优化问题。