针对突发公共事件后应急车辆调度问题, 国内外学者做了较多研究^[1–5], 文献[1]针对多地区同时发生灾害时, 需及时展开救援工作问题, 设计了改进遗传算法对所建模型求解, 并做了多组仿真实验, 结果表明, 改进算法在满足救援物资时间上取得了较大提高. 文献[2]考虑灾害发生后受损道路对车辆调度的影响, 根据道路实际状况建立了运输车辆调度模型; 文献[3]针对大规模突发事件下应急车辆调度问题,用物资未满足度的形式对灾民的损失进行量化,构建最小化灾民损失和应急车辆调度费双目标的混合整数规划模型. 文献[4]针对灾后需求不确定问题, 假设需求服从正态分布, 构建了带时间窗的应急车辆调度模型, 设计了两阶段遗传算法, 最后通过实例验证了模型与算法的有效性. 然而, 学者在建模时默认受灾点对于期待被救援的时间要求不存在差异, 本文考虑到不同程度受灾点对救援时间的紧迫性不同, 引入Sigmoid时间满意函数, 构建应急车辆调度模型. 车辆调度问题作为NP问题的一种, 一些经典的算法(动态规划、分支定界等)在寻优速度上效果欠佳, 因此需要寻找一种在可接受的时间范围内能够求得问题最优解或近似最优解的方法. 有研究表明^[6], 鲸鱼群算法(Whale Swarm algorithm, WSA)在求解高维函数优化问题上具有一定优势. 鉴于鲸鱼群算法局部搜索部分是一个随机搜索的过程, 局部寻优能力欠佳, 本文设计了基于混沌序列的搜索算子, 为了提高算法搜索效率, 又对搜索区域进行动态化处理, 进而提出了本文改进的混沌鲸鱼群算法.

1 应急情景描述与数学模型 1.1 情景描述

假定某地突发大规模公共事件, 现有一个应急救援中心,该救援中心有同种车型的救援车辆若干, 负责对多个受灾点进行救援物资的运输工作, 且受灾点的受灾程度不同. 为了确保救援任务的顺利进行, 现需要对救援车辆的路线进行合理安排, 要求使得总配送任务的平均时间满意度较高以及行驶距离较低(成本较低).

本文通过借鉴文献[7]关于地震发时被困人员生存概率函数以及文献[8]关于突发大规模疫情后城市救援有效度函数, 引入时间满意度函数来评价救援效果. 综合考虑, 采用了降指数Sigmoid时间满意度函数(见图1)对救援效果评价, 这里假设降指数Sigmoid时间满意度函数是通过Tangent Sigmoid函数截取和翻转得到, 见公式(1).

$C(t) = \left\{ \begin{gathered} 1{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t \leqslant EL \hfill \\ \frac{{2{e^{ - \beta (t - EL)}}}}{{1 + {e^{ - \beta (t - EL)}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t > EL \hfill \\ \end{gathered} \right.$

(1)

式中, $\beta $ 是正的时间敏感系数, 参数 $\beta $ 越大, 函数的敏感度越高, 即对应于受灾点对救援时间的敏感度高.

1.2 模型建立

建模之前作几方面假设: 各受灾点与救援中心的位置已知; 单个受灾点的需求量不大于单车最大容量, 且不允许出现车辆超载现象; 各受灾点的物资需求量已知且为同一类物品; 计算中只考虑单程配送问题, 即只研究救援中心到受灾点这一配送过程, 但救援车辆完成配送后需要返回救援中心; 救援中心同时对多个受灾点进行救援, 每辆车只进行一次救援, 无缺货发生.

建立的数学模型如下:

$\max \;\;\;{z_1} = \frac{{\sum\limits_{i = 1}^N {C({t_i})} }}{N}$

(2)

$\min \;\;\;{z_2} = \sum\limits_{i = 0}^N {\sum\limits_{i = 0}^N {\sum\limits_{k = 1}^M {{d_{ij}}{x_{ijk}}} } } $

(3)

s.t.

$\sum\limits_{k = 1}^M {{y_{ik}}} = 1,\;\;i = 1,2, \cdots ,N$

(4)

$\sum\limits_{j = 0,i \ne j}^N {{x_{ijk}}} = {y_{ik}},\;\;i = 0,1,2, \cdots ,N,k = 1,2, \cdots ,M$

(5)

$\sum\limits_{i = 0,i \ne j}^N {{x_{ijk}}} = {y_{jk}},\;\;j = 1,2, \cdots ,N,k = 1,2, \cdots ,M$

(6)

${s_j} = {s_i} + {t_i} + {t_{ij}},\;\;i,j = 0,1, \cdots N,i \ne j$

(7)

$\begin{split}& {x_{ijk}}({x_{ijk}} - 1) = 0,\;\; i = 0,1, \cdots ,N, \\& j = 1,2, \cdots ,N,k = 1,2, \cdots ,M \end{split} $

(8)

${x_{ik}}({x_{ik}} - 1) = 0,\;\; i = 1, \cdots ,N,k = 1,2, \cdots ,M$

(9)

${s_0} = 0, $

(10)

${C_{}}({s_i}) = \left\{ \begin{gathered} 1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;{s_i}_{}\; \leqslant E{L_{}} \hfill \\ \frac{{2{e^{ - {\beta _i}({s_i} - E{L_{}})}}}}{{1 + {e^{^{ - {\beta _i}({s_i} - EL)}}}}}\;\;\;\;{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \;{s_i}_{}\; > E{L_{}}\; \hfill \\ \end{gathered} \right.\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$

(11)

式(2)表示最大化所有受灾点的平均时间满意度; 式(3)表示最小化救援活动中车辆行驶的距离; 式(4)表示对于救灾点 $i$ 的救援任务由车辆 $k$ 来完成, 不允许分车分批运输; 式(5)和式(6)为了确保到达受灾点的救援车辆必须去下一地点; 式(7)表示车辆到达救灾点 $j$ 处的时间 ${s_j}$ 等于到达上一救灾点的时间 ${s_i}$ 加上逗留时间 ${t_i}$ 再加上 $i$ 到 $j$ 运输时间 ${t_{ij}}$ ; 式(8)和式(9)表示两个变量间的关系, 其中式(8)表示给受灾点 $j$ 实施救援的车辆 $k$ 来自 $i$ 处, 即前序节点的唯一性, 式(9)表示到达受灾点 $i$ 处的救援车辆为车辆 $k$ , 即到达某个受灾点的车辆唯一性; 式(10)表示规定每次救援任务均从救援中心出发且出发时刻记为0, 即 ${s_0} = 0$ ; 式(11)定义了受灾点 $i$ 对救援时间的满意度.

2 基于鲸鱼群算法的应急调度算法 2.1 标准鲸鱼群算法

鲸鱼群算法是通过模拟自然界中鲸鱼进行捕食等群体活动时, 通过超声波与同伴进行交流的行为特性而发展来的一种新型的元启发式算法. 当鲸鱼发现了食物源, 它会发出声音(超声波)通知附近的其它鲸鱼关于食物质量好坏和数量多少的信息. 因此, 每只鲸鱼将收到大量来自附近鲸鱼的通知信息, 然后根据这些信息移动到适当的地方寻找食物.

介绍鲸鱼群算法产生候选种群之前, 首先介绍鲸鱼发出超声波的相对强度. 并从数学角度对鲸鱼群算法的优化机理作如下描述.

定义1. 鲸鱼发出超声波的相对强度为

$\rho = {\rho _0} \cdot {e^{ - \eta \cdot d}}$

(12)

其中, ${\rho _0}$ 指超声波源的强度, 根据大量实验的经验, 对几乎所有的案例, ${\rho _0}$ 都可以设置为2. $\eta $ 为衰减系数, 它取决于介质的物理化学性质和超声波本身的属性(如超声波频率), 对于函数优化问题, 影响 $\eta $ 的的因素与目标函数的特征有关, 包括函数的维度、定义域和峰值分布. 因此在对不同的函数进行优化时, 需要设置适当的 $\eta $ 值.

鲸鱼在捕食的过程中, 如果距离“最近且较优”的鲸鱼较近, 鲸鱼将积极地向它随机移动; 反之, 鲸鱼会消极地随机移动, 因此经过一段时间的迭代就会形成一些子种群. 每只鲸鱼都随机游向“最近且较优”的鲸鱼来寻找更好的食物, 由于随机运动是鲸鱼行为的一个重要特征, 故本文采用超声波衰减的随机移动规则来探索获得新位置的迭代公式, 如公式(13)所示.

定义2. 鲸鱼x被吸引向鲸鱼y移动的位置更新公式由以公式(13)决定.

$x_i^{t + 1} = x_i^t + rand(0,{\rho _0} \cdot {e^{ - \eta \cdot {d_{x,y}}}}) * \left( {y_i^t - x_i^t} \right)$

(13)

其中, ${x_i}^t$ 、 $x_i^{t + 1}$ 分别表示鲸鱼x的第 $i$ 个元素在公式迭代第 $t$ 步与第 $t + 1$ 步时的空间位置, $y_i^t$ 表示鲸鱼y的第 $i$ 个元素在公式迭代第 $t$ 步所处的空间位置; ${d_{x,y}}$ 表示鲸鱼x与鲸鱼y之间距离, $rand(0,{\rho _0} \cdot {e^{ - \eta \cdot {d_{x,y}}}})$ 表示0到 ${\rho _0} \cdot {e^{ - \eta \cdot {d_{x,y}}}}$ 之间均匀分布的随机数.

算法实现优化的过程是: 先将鲸鱼群体随机散布在解空间中, 每一头鲸鱼所处的位置不同发出的超声波强度也不一样, 对鲸鱼进行初始化; 通过公式(12)比较, 鲸鱼积极地向超声波强度(“最近且较优”)的鲸鱼移动; 根据公式(13)来计算跟新后的位置, 这样通过多次迭代之后, 所有个体将会向超声波强度最高的鲸鱼聚集, 从而实现寻优的目的. 其中, 在每一次迭代前, 判断算法是否达到预先设计的最大迭代次数, 若达到, 算法终止, 否则继续迭代. 鲸鱼群算法流程如图2所示.

2.2 混沌鲸鱼群优化算法 2.2.1 算法原理

本文采用逻辑自映射函数对鲸鱼群算法的寻优化过程进行扰动, 可提高算法的效率. 逻辑自映射函数的数学表达式为式(14):

${L_{n{\rm{ + 1}},d}} = 1 - 2{L^2}_{n,d},\;\;n = 1,2, \cdots ,\infty ;{L_{n,d}} \in ( - 1,1)$

(14)

从上式中可以看出, 当迭代初始值不为0时, 就会发生混沌, 映射的定义域为(–1,1), 且不为0和0.5. $d$ 表示搜索的空间维度. 混沌优化的过程为: 在搜索过程的某个时刻, 鲸鱼个体 $i$ 位于D维空间内的某一位置. 根据逻辑自映射函数的性质, 首先,按照式(15)将个体空间的每一维映射到(–1,1)上. 其次, 按照式(14)进行载波操作, 从而获得了新的新混沌变量序列. 最后, 将获得的混沌变量序列按照式(16)变换到原解空间, 在此寻优过程中, 如果发现更优的解, 则将这个最新位置代替鲸鱼 $i$ 原先的位置. 否则, 进入新一轮混沌搜索, 直到搜索次数达到预先设定的次数.

${L_{id}} = 2 \times \frac{{{y_{id}} - {a_{id}}}}{{{b_{id}} - {a_{id}}}} - 1,d = 1,2, \cdots D$

(15)

$x{'_{id}} = \frac{1}{2} \times ({b_{id}} - {a_{id}}) \times {L_{id}} + \frac{1}{2} \times ({b_{id}} + {a_{id}}),d = 1,2, \cdots D$

(16)

其中, 上述两式中的 ${a_{id}}$ 和 ${b_{id}}$ 分别表示鲸鱼个体 $i$ 第d维变量的搜索上下界.

为了平衡运算时间与求解精度, 本文没有将全部的鲸鱼个体进行混沌化, 选取了部分鲸鱼作为精英个体进行混沌运算. 与此同时, 为了使得搜索效率提高, 本文使搜索区域进行动态化处理, 使算法在初期搜索范围广一些, 以免过早陷入局部最优; 在后期搜索范围小一些, 使得收敛速度加快, 用式(17)来动态收缩搜索区域.

${x_{\min ,j}} = \max \{ {x_{\min ,j}},{x_{g,j}} - \varphi \times ({x_{\max ,j}} - {x_{\min ,j}})\} $

${x_{\max ,j}} = \min \{ {x_{\max ,j}},{x_{g,j}} + \varphi \times ({x_{\max ,j}} - {x_{\min ,j}})\} $

(17)

其中, $\left[ {{x_{\min ,j}},{x_{\max ,j}}} \right]$ 表示第 $j$ 维变量的搜索范围, ${x_{g,j}}$ 表示当前处于最好空间位置鲸鱼个体第 $j$ 维对应的位置, 收缩因子 $\varphi $ 采用式(18)进行计算, $t$ 是当前算法的搜索次数.

$\varphi \left( t \right) = 1 - \frac{1}{{1 + \exp (4 - 0.04t)}}$

(18)

2.2.2 基于问题的鲸鱼编码与译码

鲸鱼群算法是一种基于种群的全局智能寻优方法, 搜索空间中的每一头鲸鱼都对应寻优空间中的一个解, 且对应一个目标函数值. 每头鲸鱼的位置向量X是一个2N维的向量, 前半部分N维向量用来定义受灾点所使用的车辆, 每一维随机数取(0,M–1)之间整数(向上取整), M表示车辆数, 后半部分N维向量用来定义车辆的救援顺序, 每一维随机数取(0,1)之间的实数, 对于同一辆车救援多个受灾点, 救援顺序由这N维基因位置上的数值决定, 数值大的先救援, 数值相等时, 近左优先. 假设有8个受灾点, 3辆救援车辆, 某一鲸鱼的向量表示如图3.

前面8个基因位使得受灾点对应分配一辆救援车辆, 8个受灾点对应的救援车辆分别为1、2、1、2、3、1、2、3, 从而可以得出车辆1救援受灾点1、3、6, 救援顺序则根据后面L维向量基因位的数值大小, 跟据上文所述, 车辆1的救援安排为0-1-6-3-0. 同理, 车辆2的救援安排为0-4-2-7-0, 车辆3的救援安排为0-5-8-0.

2.3 适应度值计算

为了计算鲸鱼个体的适应度值, 综合考虑 ${z_1}$ 和 ${z_2}$ 两个目标, 给出计算公式如下:

$fitness{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} z = \lambda \cdot \frac{1}{{{z_2} + \varepsilon }} + (1 - \lambda ) \cdot \frac{1}{{{z_2} + \varepsilon }} \cdot {z_1}$

(19)

其中, $\varepsilon $ 为一极小正数, 防止分母为零; $\lambda $ 为决策偏好因子, 取值(0,0.5), 这样取一方面考虑了应急物流时间紧迫性的特点, 使得对时间满意度赋予了较高的权重, 另一方面又兼顾了应急物流的弱经济性. 适应度函数前半部分为与成本有关的函数, 后半部分为与时间效用有关的函数, 分别表示救援成本和救援有效度对适应度函数的影响. 从上述公式可以看出, 在距离函数 ${z_2}$ 较小、救援满意度函数 ${z_1}$ 较大时适应度值越大. 与一般的多目标加权求和适应度函数不同, 本文适应度函数考虑到不同量级的目标函数对适应度值影响程度的差异性, 因此在后项中引入调节因子 $\displaystyle\frac{1}{{{z_2} + \varepsilon }}$ , 以均衡 ${z_1}$ 和 ${z_2}$ 对适应度函数的影响.

2.4 算法步骤

综上所述, 混沌鲸鱼群算法的基本步骤如下:

Step 1. 初始化算法基本参数: 设置鲸鱼数目Q, 衰减系数 $\eta $ , 精英群体比例 $n\% $ , 混沌搜索迭代次数H, 最大搜索次数Iter_MAX或搜索精度.

Step 2. 在搜索区域内随机初始化鲸鱼的空间位置, 根据所处位置计算鲸鱼的目标函数值, 将其作为各自的最大超声波强度 ${\rho _0}$ .

Step 3. 根据式(12)计算群体中鲸鱼的相对超声波 $\rho $ , 根据式(13)更新鲸鱼的空间位置.

Step 4. 对鲸鱼个体进行评估, 选取性能最好的 $n\% $ 的个体作为精英, 采用混沌优化策略进行优化; 选取性能最差的 $n\% $ 的个体, 重新随机产生新的鲸鱼个体予以替代.

Step 5. 根据式(17)动态收缩搜索区域.

Step 6. 根据移动后的鲸鱼位置, 重新计算各鲸鱼的最大超声波强度.

Step 7. 当满足最大搜索次数或达到搜索精度时, 转至Step 8, 否则转Step 3, 进行下一次搜索.

Step 8. 输出全局最优个体与全局极值点, 算法结束.

3 实验仿真及分析

本文实验环境为: PC计算机处理器主频2.3 GHz, 酷睿i7 3610QM处理器, 内存8 GB, 在Win10系统下的Matlab 2010a 的编程软件. 针对本文所设计的突发公共事件下的应急车辆调度问题进行对比实验. 鲸鱼群算法中所涉及的各种参数设置目前均没有严格的理论依据, 故本文所设置的参数值都是经过反复实验最终确定如下: 设置鲸鱼数目Q=100, 衰减系数 $\eta $ =1.5, 超声波源的强度 ${\rho _0}$ =2, 精英群体比例20%, 混沌搜索迭代次数50, 最大搜索次数Iter_MAX=200. 模拟退火算法是一种发展较成熟的启发式算法, 学者们利用此方法求解车辆调度问题取得了一定成果^[9,10], 本文拟采用模拟退火算法进行对比实验, 相关参数设置为: 初始温度T₀=1000, 下降比率p=0.95, 终止温度T_end=10, 内循环次数N_nei=200.

3.1 算例设计

现假设某地区发生地震, 政府部门需从一个应急救援中心向多个应急救援待救点进行救援物资的配送, 救援中心有多辆配送车辆可以利用. 由于缺少真实的实验数据, 本文构建了三组实验数据(见表1、表2、表3), 针对三组实验分别采用载货量为80、130和160单位的车辆配送, 每组实验仅采用一种车型且均以75 KM/H匀速行驶. 本文考虑四个受灾等级(灾情1-灾情4), 数字越大表示受灾点对救援时间的要求越高, sigmoid时间满意度函数需要设置四个不同的 $\beta $ 值,分别设置为0.1、0.2、0.3、0.5, EL设为1H,即救援车辆在1小时之内到达受灾点的满意度为100%, 超过1小时后会随着时间的推移, 不同等级的受灾点的满意度会相应降低.

3.2 实验结果及分析

对三个算例分别采用模拟退火SA、鲸鱼群算法WSA和混沌鲸鱼群算法CWSA各独立求解20次,对三组方案的求解效果进行比较, 求解结果如表4所示. 由于篇幅有限, 文章仅列出了三种算法针对25个受灾点最优实验结果, 如最优行驶方案对比(表5)、最优行驶路径(图4–图6)以及三种规模迭代对比(图7–图9).

从表4可以看出, 在处理较小规模车辆调度情况下, 三种算法处理能力差距较小, SA求解的稳定性上相对有优势、但不明显. 随着求解规模的扩大, 改进的鲸鱼群算法效果较明显. 相比原始鲸鱼群算法, 无论在满意度与行驶距离上都得到了的改进, 且求解结果较为稳定.

针对25个受灾点的救援问题, 表5给出了三种算法最优寻优结果, 可以看出: 三种算法给出解的救援满意度均较好; CWSA算法相比另外两种算法, 救援成本较低. 综上, 针对基本鲸鱼群算法的改进策略是有效的. 此外, 对于14个受灾点的应急救援问题, 三种算法求得的最优解满意度较低, 这可能是由案例本身决定的. 救援车辆数目的选取对受灾点的满意度存在影响, 当增加救援车辆的数目时, 可以提高救援满意度.

4 结束语

本文通过分析大规模突发公共事件后特征, 引入时间满意度函数, 建立了平均满意度与行驶距离为目标的数学模型, 针对鲸鱼群算法局部搜索能力较弱, 设计了改进的鲸鱼群算法对模型进行求解. 通过实验仿真, 验证了模型及改进算法的有效性. 应急救援具有复杂性、动态性、难解性等特点, 建立不同的应急车辆调度模型应对不同救援场景, 并设计更为有效的算法将是我们今后研究的重点.