薛星群，王旭坪,2，韩 涛，阮俊虎

(1.大连理工大学商学院，辽宁 盘锦 124221；2.大连理工大学系统工程研究所，辽宁 大连 116024； 3.西北农林科技大学经济管理学院，陕西 杨凌 712100)

1 引言

近年来环境问题日益严峻，自然灾害频发，比如2005年南亚大地震、2008年艾奥瓦洪灾、2010年智利地震以及由地震引发的海啸等，这些自然灾害不仅使人们的财产受到损失，还危及人们的生命安全。在灾后救援实施工作中，及时有效地向受灾地区运输救援物资有利于降低灾区的受损程度。但是通常情况下，灾害发生会使道路呈现出不同程度的损害致使救援工具无法到达被救援地区，这不利于应急工作的开展。而救援任务时间紧急，在此背景下，相对于之前单一救援工具，运用车辆-直升机联合运输救援物资逐渐成为对受灾地点救助的主要应急方式[1-2]。

目前已有部分学者从单一救援工具的角度出发研究应急物资的运送问题。代颖等[3]考虑到车辆运载量的限制，为实现运送救援物资至各受灾点的总时间最少，构建物资运输点选址和车辆路径安排集成优化模型，设计一种两阶段启发式算法予以解决；詹沙磊和刘南[4]在研究多车型、多受灾点应急物资的运送问题中，将贝叶斯分析和最优停止决策规则纳入同一研究框架，通过建立贝叶斯风险函数解决实际问题；Mete和Zabinsky[5]根据地震后受灾地区实况建立一种随机优化方法来解决救援工具配送救援的医疗物资问题；Garrido等[6]为了在灾涝后帮助决策者合理调度救援物资而设计出一种量化应急物资的库存水平以及车辆可利用性的优化模型；Moreno等[7]研究出一种混合整数规划启发式方法来解决阶段灾害救援中车辆运力的重复使用问题；朱建明等[8]针对应急医疗物资的运输问题，考虑到医疗物资到达医院时间最小以及最大限度的满足医疗物资所需，构造出改进的随机算法进行模型求解。以上研究在应急物资运送模型构建的过程中，大部分只是针对车辆等陆地救援工具进行研究，关于多类型救援工具的研究相对较少。现实生活中，考虑到地震、海啸等大规模自然灾害实况，道路的流通性较差且余灾时有发生，仅靠车辆无法对灾区实施救援工作。因此，从多种类型救援工具的角度出发，Fikar等[9]研究了如何选取陆运和空运协调运输中转点等问题，根据仿真和优化方法设计出一种决策支持系统以应对最后一公里的救援物资的运输问题；Zhou Chao等[10]针对突发情况设计出多智能体方法来予以应对。但是上述多类型救援工具下的决策系统和方法仍不能解决大规模受灾的问题。为此，祁明亮等[11]以及阮俊虎等[12]分别采用直升机或者建立救援中转点方式，分两阶段运输救援物资以此处理道路受损情况下对难以进入救援区配送救援物资的问题，但是两阶段运输物资会对救灾时间产生影响，而且直升机在协同车辆帮助普通受灾地区运输物资的过程中剩余运力的使用率也不高。

多种运输工具的联合运输物资在应对大规模灾害的过程中发挥重要的作用。李双琳等[13]设计出一种NSGAII来解决模糊需求条件下救援物资配送选址及多式联运问题；Hu Zhihua[14]构建整数线性规划模型进而改进了免疫算法,优化了救援网络下的集装箱多式联运问题；郑斌等[15]基于上下层的双层规划构建模型并设计出混合遗传算法来解决物流选址及联运问题。以上文献大部分针对多种运输工具能够同时使用条件下的救援物资配送，但由于灾后情况的多变性，应急救援决策也需要适时调整，此时运输方式的选择受到限制，上述多种运输工具同时使用的条件也不能满足。王旭坪等[16]在运力受限的条件下通过混合整数规划模型并改进编码方式来配置运力；Yi和Özdamar[17]通过设计的集成位置分布模型来解决撤离人员问题和应急物资的运输问题。然而以有限运力为基础的紧急决策是分阶段的，而且运力的多次利用也具备约束性特点。Erdemir等[18]为了解决空运工具同陆地车辆联合实施救助等问题，设计出一种位置覆盖模型和贪婪启发式算法；Ruan等[19]基于模糊C-均值聚类方法，考虑医疗救助点数量和容量限制约束，构建平衡的车辆联运网络和直升机来解决大型自然灾害中的医疗物资多式联运问题。但是上述研究没有考虑运力的可重复利用，因此尚有可优化空间。

为了弥补上述研究的不足，我们将通行约束以及运力受限这两个条件考虑在内，在灾害造成通行受约束和救援工具数量有限的场景下，采用“车辆+直升机”的联合运送方式，其中孤立的受灾区只允许直升机到达，其他受灾区车辆和直升机均能到达，考虑救援的经济成本，对完成运送任务的运力重新进行救灾任务分配，构建满足灾后救援物资联合运送平均等待时间最短及决策体系总成本最低的优化模型，且基于带精英策略的非支配排序遗传算法(NSGA-II)改进了随机邻域搜索以及分级交叉，从而提出新的智能求解算法(NSHEA-II)，通过算例实验表明，本文设计的算法在有限的时间内得到更优的Pareto前沿，能够为不同偏好情况下的决策者提供更有效的调度方案。

2 模型构建

2.1 问题界定与模型假设

一般而言，灾后救助使用直升机和车辆运输救援物资。但是通常发生大规模灾害时，通向受灾地点的主要路径受到阻碍或是中断无法通行，由于道路疏通具有一定的时滞性或者短期内无法疏通，从而使一部分受困灾点的救援物资运输仅依赖高性能、运输路线灵活的直升机来完成，即特殊受灾点由直升机完成物资运输，一般受灾点依靠直升机的剩余运力和陆面救援车辆的的联合运输完成，具体运送路径如图1所示。同已有研究中直升机只进行孤立受灾区救援物资运送不同的是，本文考虑了直升机可以利用其剩余运力来帮助地面车辆运输救援物资至普通受灾区，减少了救援时间，提高了所有运力的利用率。本文综合考虑受灾地区获得救援的平均时间及应急决策体系下整体费用最小化这两个决策目标，合理安排直升机和地面救援车辆运力以实现运输应急物资至受灾点的高效化。

图1 “直升机+车辆”联合运送示意图

考虑到能及时调动的应急救助交通工具数量的限制性特点，一些直升机或是车辆救援工具可能会被反复派出进行多次物资的运送，假定一条子路径代表救援时一次任务，那么一个救援工具的所有运输路径是由其全部运输子路径构成。灾后，以集散中心作为起点，直升机以及救援车在完成一条子路径的物资运输后重返集散中心，然后进行下一次的物资配送直至任务完成。主要假设如下所述：

(1)假定决策系统中包含待通行约束的多个受灾点，受灾点的位置和两两受灾点之间的距离已知；还包含一个物资集散中心。

(2)本文对物资分配和多样化资源的配送不作细致性研究，所以认为物资集散中心的物资量充足而且种类单一，故救援工具仅运输一种物资；

(3)遥感飞机能够观测到受灾以后运输物资的路径状态，因此受灾点的种类以及所需救援资源的数量已知，并且各受灾点仅能被访问一次；

(4)假定能够应急调度救援的直升机、车辆的装载量和数量均已知，平均速度不变，运行时长不受限制，不考虑救援工具维护时间的限制。

2.2 符号定义

Di：节点i对救援物资的需求量；

dij：节点i到节点j的距离，且满足dij=dji，i,j∈V；

Nk：配送工具k的子路径集合；

Ck：配送工具k的单位距离运送成本；

Qk：配送工具k的容量上限；

si：配送工具在受灾点i工作的时间，s0=0；

ck：配送工具k的固定成本费用

vk：配送工具k的平均速度；

其中(i,j)∈A,k∈K,n∈Nk；

2.3 数学模型

通过对已有的应急救援进行研究，我们发现灾后救援物资的及时运输以及救灾的效果会受到救援效率的影响。救援物资(如医疗救护物资、电力物资等)的长时间匮乏不仅会错失抢救伤亡人员的最佳时间，而且会给灾民造成一定的心理恐慌，不利于应急救援工作的开展，所以本文将最小化受灾地区等待物资运送的平均时间设定为第一个目标。

(1)

(2)

救援物资协同调度优化模型如下：

Minf1

(3)

Minf2

(4)

(5)

(6)

(7)

(8)

(9)

∀i,j∈V,k∈K,n∈Nk

(10)

(11)

(12)

(13)

3 算法设计

割平面法和分支定界法是解决小规模整数规划问题的主流方法，然而由于实际情况越来越复杂，一定的时间限制下难以求出精准算法的解，通常利用遗传算法、模拟退火、粒子群算法等解决大规模调度问题以此提升运算效率。对于多目标优化问题，基于将多目标问题转变为单目标问题的思路，通常采用极小极大法、约束法等方法，然而权重值以及所给目标的次序会对其产生影响。对此，根据上述设立的模型的复杂度，我们运用带精英策略的非支配排序遗传算法(NSGA-II)计算其解，另外对随机邻域搜索和分级交叉进行改进，构造出带精英策略的非支配排序混合进化算法(NSHEA-II)，算法流程图如下所示。

图2 NSHEA-II算法流程图

3.1 染色体编码

我们采用二维矩阵实数编码方式，行表示救援工具，列表示受灾点的编号，为了满足式(7)的条件，染色体在初始化以及邻域操作的过程中始终确保带通行约束的受灾点仅分配给直升机的规则进行。假设有2架直升机，3辆卡车，20个受灾点，且编号为1、2、3的受灾点带通行约束(矩阵中0在编码和解码时忽略)，则其可能的编码方式如图3所示。由于每个救援工具有装载量限制，解码时按照应急工具的载容上限划分子路径，如卡车2的子路径可以是0-4-6-0和0-6-10-0。

图3 染色体编码图

3.2 遗传算子

这部分运用锦标赛选择法，从种群中选取50%的染色体作为父代种群进行交叉和变异。

通过对模型分析，我们不但要对救援工具的种类进行确定，还要对救援路径予以确定，以此获得应急物资的联合调度优化。而传统带精英策略的非支配排序遗传算法中变异算子形式单一，对当前解易产生较大的不稳定性，对解的局部优化产生了负面影响，为了弥补传统算法的不足，我们设立如下的邻域结构：

变换1：将某一应急工具路径中的一个节点插入到另一个应急工具路径中。具体做法：首先在一个应急工具的路线中任意抽一个节点a，然后在该路径中将节点a去除，接着查看该节点是否为带通行约束，如果满足条件那么将其放入直升机路径，如果不满足条件，随意将其插入某个路径。

变换2：交换某两个应急工具路径的节点。如果应急方式选择的是一架直升机、一辆车，并且从直升机路径中选带通行约束的受灾点，那么需要从直升机路径中剩余点中选择另外一个互换点，如果该直升机路经节点都是带通行约束的节点，那么需要重新选择路经。

变换3：优化某应急工具路径内的两路径。在满足总节点数大于3的应急工具路径中任意选择两个节点a和b，a之前和b之后的路径不做任何改变地加入至新路线，将a到b两节点间的路线翻转编号后加入至新的路线，从而优化现有路径[20]。

变换4：交换一条应急路径中的任意两个受灾点。在受灾点数多于3的应急路径中任意抽取两点a和b，将a与b互换。

3.3 非支配排序

第一步：将种群内全部ni=0染色体存到集合z1中；

第二步：考虑目前集合z1中的任意染色体i所支配的染色体集合si，然后将si中的每条染色体k的nk变成nk=nk-1(处于支配地位的染色体i已经存到集合z1中)，如果更新后的nk=0，那么把k保存在另一个集合z2。

第三步：集合z1中的染色体处于支配地位不受其他染色体支配，把它们认定为第一级最优的非支配染色体集合，给集合z1中的染色体一个相等的非支配序ranki=1，接着对z2重复以上分级操作并给予相应的非支配序，直至全部个体都被分级。

3.4 拥挤度比较

为了得到均匀的Pareto最优解集，需要计算染色体的拥挤度。在非支配排序和拥挤度计算之后，每条染色体i都拥有两个属性：拥挤度disi和非支配序ranki。因此可得到拥挤度比较算子：任选两条染色体i和j，如果染色体i所处的非支配层相对于染色体j所处的非支配层是占优的，即rankidisj，那么染色体i占优。染色体拥挤度的计算步骤如下：

第一步：初始化每条染色体的拥挤度disi=0；

第二步：将染色体按照f1和f2值排序，取最大的目标函数值fmax和最小的目标函数值fmin，把f1和f2的拥挤度设成无穷大；

第三步：运用式(13)计算其它染色体的拥挤度：

(14)

4 算例分析

我们将物资集散地的坐标设定成(50,50)，救灾点的坐标信息在100100的平面内随机生成，我们用欧氏距离描绘任意两受灾点间的距离。考虑通行以及运力受约束的情况，我们任选三个受灾点作为带通行约束的受灾点，并设定存在2架直升机、3辆卡车，表1为应援工具参数，表2 为受灾点信息。测试环境为Windows7操作系统，4GB内存，intel四核2.5GHz处理器，采用MATLAB 2015b编程算法。

4.1 算例测试结果

将NSHEA-II算法的参数设置为：染色体数目N=100，最大迭代次数Gmax=1000，最大邻域搜索次数Cmax=10，交叉概率Pc=0.8，变异概率Pm=0.2。图4和图5描绘了算法的收敛效果，由图4可知，初始种群的解分布有多个等级，具有许多可支配解，最少的平均等待时间是62分钟，最低的经济成本2719万元。当达到最大迭代次数后，全部解均为非支配解，从而形成Pareto最优前沿。而图5中平均等待时间最小值是39分钟，最小的经济成本是1770万元，优化比例分别达到37.1%和34.9%。

表1 应援工具参数

表2 受灾点信息

图4 初始解分布图

图5 最终解分布图

表3和表4为两种目标下的Pareto最优方案。表3为决策者偏好系统的总经济成本最小的最优方案，表4为受灾点平均等待时间最短的最优方案。其中，路径中的数字0代表初始点，括号内的数字表示从0时刻开始直至救援工具到达该受灾点的等待时间，数值均向上取整。以表3为例，我们对调度路径方案的含义进行简要说明，在决策者偏好系统总经济成本最小时，直升机1在第2分钟到达受灾点17完成运送任务，直升机2在第18分钟到达受灾点9，在第22分钟到达受灾点15，第35分钟到达受灾点2完成运送任务。卡车1有两条运送路径，第一条为0-16-12-5-1-6-0，第二条为0-3-0，卡车1完成第一条子路径的运输任务后才能向受灾点3运输救援物资，卡车1从0时刻开始直至到达受灾点3需要242分钟。剩余救援工具的具体调度方案解释同上，不再赘述。带通行约束的受灾点2、9、17都由直升机负责运输，符合式(7)。

表3 总成本最小

表4 平均等待时间最短

由图5中Pareto前沿的两端可以观察到，两个目标函数值并非呈等比例的反比关系，当牺牲很小一部分目标函数2时，可以大幅度优化目标函数1，反之亦然。当耗费的总经济成本由1770万元提高到1792万元时，提升比例为1.24%，但受灾点平均等待时间由86分钟降低为58分钟，下降比例为32.56%。由此可见，利用多目标进化算法得到Pareto最优解集是有效的，能够为决策者提供科学的应急物资调度方案。

为了检验算法的优越性，我们将用NSGA-II算法对算例求出的结果同用本文提出的算法求出的结果进行比较，图6显示了采用两种算法迭代到最大次数后得出的Pareto最优解集。根据图6显示的结果，我们可以看出，NSHEA-II的对应解集支配了NSGA-II的对应解集，NSHEA-II的寻优效果较大。运用NSGA-II求出的平均等待时间的最小值为46分钟，经济费用的最小值为1851万元，而用NSHEA-II求出的最短平均等待时间为39分钟，最小经济成本为1770万元，后者相对于前者的优化比例分别是15.2%、4.4%，因此，我们可以通过使用NSHEA-II制定更优的应急决策。

图6 NSGA-II与NSHEA-II寻优效果对比

本文分别以受灾点的平均等待时间最小化和应急体系下总经济成本最小化为目标，将实验重复进行10次，NSHEA-II和NSGA-II迭代到一定次数后得到的结果如图7 所示。根据图7显示，在一定时间约束的情况下运用NSGA-II得出的两个目标解不稳定，而运用NSHEA-II获得的目标解的相对来说较稳定，主要因为改进后的算法确保了染色体的多样性，并且多种邻域搜索算子使NSHEA-II规避了算法容易收敛到局部最优的风险。

图7 NSGA-II与NSHEA-II稳定性对比

本文构造出15组不同参数的算例以此检验算法的可行性，其中，num表示受灾点数量、p表示带通行约束的受灾点数量、k1表示直升机数量、k2表示车辆数量。由表5可以看出，通过结果对比我们发现，相对于NSGA-II，NSHEA-II取得的结果较好。通过对15组算例的目标一和目标二变量的计算，我们发现相对于比NSGA-II，NSHEA-II提高了3.83%和4.41%的优化。另外，我们还得到如下结果：通过对比只有受灾点量不同的算例我们发现，受灾点的个数与最小平均等待时间和最低总成本费用存在正比关系；通过对比结果中满足只有带约束受灾点数量不同的算例可以发现，在两种算法的结果中，带通行约束的受灾点数量同应急网络中的最小总成本成正相关，而与最短平均等待时间无明确的关系；通过对比结果中仅直升机数量不同的算例可知，直升机的数量同受灾点最短平均等待时间成负向关，同最低成本费用呈正向关系；根据仅车辆数量不同的算例，我们认为车辆的数量与最短平均等待时间呈现微弱的负向关系，而与最低成本不具有明确的关系。

表5 NSGA-II与NSHEA-II计算结果对比

5 结语

基于多目标的大规模灾害场景，我们分析了空运和陆运两种救援方式联合运输的调度问题。由于救援工具运力受限，在带通行约束的受灾地区仅准许直升机运送物资，在这样的条件下，本文将受灾点平均等待时间最短以及工具调度总成本最小作为主要目标，制订了决策者不同偏好情况下每架直升机和每辆卡车的救援路线，能够在较短时间内生成“直升机+车辆”联合调度方案，为特定灾害应对情景提供物资运送方案。本文在构建模型求解的过程中，按照参考分级交叉构造出新的交叉算子，然后综合使用多种邻域算子设计出随机邻域搜索的变异算子，构造新的多目标求解算法。通过实验检验，我们发现相对NSGA-II算法，NSHEA-II算法的结果更好，这对决策者制定相关应急决策具有一定的借鉴意义。

在后续的工作中，我们会考虑将受灾点服务的公平性这一条件纳入研究之中，也就是说，以平均等待时间最短为条件，在此基础上考虑满足各受灾点等待服务的时间差异性最小条件，以此规避受灾点等待物质运送的时间过长等风险，从而满足实际要求。此外，随着迭代次数的增多，本文设计的变异方法会减弱算法运行的速度，需要将算法的自适应性考虑在内。