大规模公共卫生事件下城市即时配送网络优化模型与算法

摘要

大规模公共卫生事件下城市即时配送存在路网数据失真、供需侧信息不确定及网络中断问题，亟需考虑信息不确定性及路网中断可能性进行城市即时配送网络优化。首先，考虑封控导致道路限行下的路网构建问题，建立城市底层路网并提出改进的Floyd算法；其次，针对开放式多配送点的城市即时配送问题，考虑供需不确定性及设施服务中断问题，使用蒙特卡洛模拟方法构造情景树，建立多目标随机规划模型并设计混合进化算法求解；最后，以2022年上海新冠肺炎疫情事件为例，发现大规模公共卫生事件导致配送设施服务能力、路网容量及客户需求突变，配送系统容易因供需不匹配而发生“爆单”“爆仓”，但一方面设施服务中断未必导致配送成本增加，而是通过降低客户满意度来增加总成本，另一方面更多的车辆使用数目未必导致总成本增加。

关键词

大规模公共卫生事件; 中断; 即时配送; 多目标随机规划模型; 蒙特卡洛模拟; 混合进化算法

大规模公共卫生事件下，政府多根据应急预案启用重大突发公共卫生事件一级响应，限制人员和车辆等的流动，给城市即时配送带来严峻挑战^［1］。一方面，临时管控使得道路通行能力降低甚至中断^［2］，然而地图软件或商业配送软件的路网状态更新不及时，严重影响实际路网数据获取，导致配送服务难以满足客户的个性化服务需求，据此得到的配送策略会产生较大偏差。另一方面，民众恐慌心理导致应急物资需求存在较大不确定性，往往发生某类商品需求量超出系统能力而产生“爆单”“爆仓”现象，导致配送中心陷入失效困境，且配送中心的运行状态及配送能力也因区域管控存在不确定性、供需不匹配，导致既定配送方案无法正常完成^［3］。因此，系统考虑大规模公共卫生事件下应急设施中断及供需不确定性，优化城市即时配送网络，成为当前亟待解决的重要理论和现实问题。

传统的应急物资配送问题具有动态性和不确定性特征，已有研究聚焦于灾前预测、灾时疏散^［4］、灾后救援^［5-7］三个维度。针对灾前预测，部分研究从约束方法、参数及目标不确定性角度进行考虑^［8］，也有研究应用博弈论从应急准备角度考虑^［9］。针对灾时疏散，多采用设立损失函数、建立随机模型的方法进行规划^［10-11］。针对灾后救援，多考虑路网状态等现实因素，建立多阶段规划模型进行求解^［12］。但由于大规模公共卫生事件具有持续周期长、传播范围大、爆发性强等特点^［13］，即时配送网络优化问题对时效性和确定性要求更高，无法按照传统应急物资配送问题进行处理，已有研究多从事件整体角度考虑^［14］，或从灾后救援角度处理此类问题^［15-16］。

可见，突发事件下的应急物资配送网络优化已有较多研究，但对于大规模公共卫生事件下的城市即时配送网络优化研究仍然有限。本文研究大规模公共卫生事件下城市即时配送问题，主要内容包括：①考虑大规模公共卫生事件道路交通网络的动态性与不确定性，结合路网建立最短路径求解模型，突破直接利用坐标求取网络任意两点间距离的不足，使得求解结果更加准确实用；②综合考虑需求侧与供给侧路网状态、配送中心状态、配送车型、客户需求等的不确定性，针对开放式多配送点设施服务中断问题，构建多目标随机规划模型；③基于蒙特卡洛模拟方法生成情景树，设计混合进化算法，对不确定情况进行再分配，有效提高求解效率。

1 模型构建

在疫情缺乏准确预测导致交通设施破坏程度和受灾人数规模不确定性下，系统考虑路网数据的不确定性、供应能力及客户需求的随机性、设施中断的可能性，建立随机规划模型，研究框架见图1。

图1  研究框架

Fig. 1  Framework of research

首先，考虑封控导致的路段管控或中断问题，构建城市配送道路交通网络进行数据测算；其次，考虑疫情预测困难导致的配送中心储量、客户需求量、客户时间窗等的不确定性，引入情景树通过有限概率空间$\{\omega ,P\}$来表示固有不确定性，使用蒙特卡洛模拟方法生成一系列情景，情景集合$\omega$包括$n$个情景，其中第$s$个情景${\omega }_s$发生的概率为$P({\omega }_s)=P_s$，每个情景包括不同的输入参数组合（配送中心储量、客户需求量、客户时间窗）；最后，考虑封控导致的设施中断问题，在情景中融入配送中心中断概率，针对中断情景重新规划配送方案。

1.1　模型假设

（1）订单不可拆分，所有客户都需要被服务，每个客户仅接受一个配送中心的一辆车的一次配送服务，且所有客户的地理位置和需求量已知，顾客的服务时间窗为软时间窗，当配送车辆到达客户点的时间早于该客户点的最早时间或晚于该客户点的最晚时间时，需支付一定的惩罚成本。

（2）考虑单一类型应急物资（如预包装食品、水、医疗物资），以独立包装供应。

（3）车辆为常规货运车辆，运输能力有限，一辆车可服务多个客户点，配送服务过程中客户需求量要满足车辆容量限制，配送车辆只在行驶时间内产生油耗，车辆在完成配送任务后无需返回，选择就近的配送中心停靠即可。

（4）客户满意度仅与时间窗及需求满足有关。

（5）开放数量固定的多个配送中心，如果配送中心失效，服务的客户将分配给其他配送中心。

（6）通过情景树表示未来可能的情景，每一情景包括不同参数组合（配送中心储量、客户需求量、客户时间窗、配送中心中断概率）。

1.2　道路交通网络构建

大规模公共卫生事件下的封控管理易导致某些场所、路段的限行或中断^［17］。传统地图软件或商业配送软件不仅会产生不准确的路网数据，还可能满足不了配送企业实时变化的目标，有必要构建城市道路交通网络并对路网数据进行测算。路网搭建过程包括：提取开源地图软件OpenStreetMap某一区域交通路段，利用ArcGIS获取路网数据，结合当前道路交通政策，以道路交叉口作为节点、路段作为连接节点的边，经过映射得到道路交通网络。

一方面，城市道路的特殊性使得两点间的直线距离无法匹配具体路径，直接使用坐标距离会导致数据明显失真，极大降低模型的精度和应用性；另一方面，目前深度学习用于处理路网问题还存在较大误差^［18］。因此，本文针对构建的道路交通网络，提出改进的Floyd算法以求解最短配送路径，克服传统方法求解效率低、准确性低的问题。由于城市即时配送是在满足客户需求的条件下追求成本最低，因此最短距离如式（1）所示。

$$\mathrm{m}\mathrm{i}\mathrm{n}{d}_{ij}=\sum _{i,k\in I\bigcup O}{d}_{i,k}+\sum _{k,j\in I\bigcup O}{d}_{k,j}$$

（1）

式中：$I$为客户集合；$O$为配送中心集合；$i,j,k$分别表示出发节点、终止节点、中间节点；$d_{i,j}$为节点$i$到节点$j$的距离。出发点与目标点之间的距离由其关系来决定，如式（2）所示。

$$d_{ij}=\left\{\begin{array}{l}{d}_{i,j},\hspace{1em}\hspace{1em}\hspace{1em}i=j\\ d_{i,j},\hspace{1em}\hspace{1em}\hspace{1em}{d}_{i,k}>d_{i,j}\mathrm{或}{d}_{k,j}>d_{i,j}\\ d_{i,j},\hspace{1em}\hspace{1em}\hspace{1em}k=i\mathrm{或}k=j\\ d_{i,k}+d_{k,j}, d_{i,k}+d_{k,j}<d_{i,j}\end{array}\right.$$

（2）

1.3　多目标随机规划模型

1.3.1　目标函数

本文综合考虑配送成本最低、客户满意度最高两个目标，目标函数如式（3）所示。针对每一可能的情景，配送成本表示配送过程中的所有花费，包括配送中心开放成本$c_{\mathrm{1}}$、车辆派遣成本$c_{\mathrm{2}}$、配送成本$c_{\mathrm{3}}$，分别如式（4）~（6）所示；客户满意度$c_{\mathrm{4}}$表示客户对配送结果的评价，取决于需求物资是否全部得到满足以及车辆是否在规定的时间窗内到达，其中时间惩罚成本为违反时间与时间惩罚系数$a$的线性乘积，需求惩罚成本为未满足需求量与需求惩罚系数$b$的线性乘积，如式（7）所示。

$$\mathrm{m}\mathrm{i}\mathrm{n}z=\sum _{s\in S}{P}_s(\mathrm{1}-\alpha )(c_{\mathrm{1}}+c_{\mathrm{2}}+c_{\mathrm{3}})+P_s\alpha c_{\mathrm{4}}$$

（3）

$$c_{\mathrm{1}}=\sum _{o\in O}{c}_o{\beta }_{s,o}$$

（4）

$$c_{\mathrm{2}}=\sum _{v\in V}\sum _{i,j\in I}\sum _{o\in O}(c_v^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}{x}_{s,o,j,v}+c_v^{\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}{x}_{s,i,j,v})$$

（5）

$$c_{\mathrm{3}}=\sum _{i,j\in I}\sum _{v\in V}{d}_{i,j}{f}_v{c}_v{x}_{s,i,j,v}$$

（6）

$$\begin{array}{l}{c}_{\mathrm{4}}=\sum _{i\in I}\sum _{v\in V}\left(a\mathrm{m}\mathrm{a}\mathrm{x}\left\{\frac{t_{i,v}-t_{s,i}^{\mathrm{w}\mathrm{i}\mathrm{d}\mathrm{o}\mathrm{w}}}{t_{s,i}^{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{w}\mathrm{i}\mathrm{d}\mathrm{o}\mathrm{w}}-t_{s,i}^{\mathrm{w}\mathrm{i}\mathrm{d}\mathrm{o}\mathrm{w}}},\mathrm{0}\right\}\right)+\\ \hspace{1em}\sum _{i\in I}\left(b(q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}-q_{s,i}^{\mathrm{d}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}})\right)\end{array}$$

（7）

式中：$s\in S$为情景集合；$i,j\in I$为客户集合；$v\in V$为车辆集合；$P_s$为情景发生概率；$\alpha$为量化客户满意度的风险权重；$c_o$为配送中心的固定开放成本；${\beta }_{s,o}$表示配送中心$o$在情景$s$下是否正常运行；$c_v^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}$为车辆的发车费用；$c_v^{\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}}$为车辆的节点费用；$x_{s,o,j,v}$表示车辆$v$在情景$s$下是否经过$o-j$弧，同理$x_{s,i,j,v}$表示车辆$v$在情景$s$下是否经过$i-j$弧；$d_{i,j}$为从客户$i$到客户$j$的距离；$f_v$为车辆$v$的燃油消耗量；$c_v$为单位燃油价格；$a$、$b$分别为时间惩罚系数及需求惩罚系数；$t_{i,v}$为车辆开始服务客户的时间；$t_{s,i}^{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}}$、$t_{s,i}^{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}}$分别为情景$s$下客户的单侧时间窗及可接受的最迟服务时间窗；$q_{s,i}^{\mathrm{d}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}}$为情景$s$下客户的需求量；$q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$为情景$s$下客户的实际送货量。

1.3.2　约束条件

相比传统多目标优化模型，使用蒙特卡洛模拟的方法构造的概率统计模型具有较高的精度^［19］，原因在于蒙特卡洛通过大量的随机样本进行模拟，使结果更精确地反映系统的期望效果，并且可以通过增加随机样本数量的方式提高系统评估效果的准确性。式（8）表示情景发生概率，之后在每个情景中有如下约束：式（9）表示配送中心流量平衡约束；式（10）表示配送中心容量约束；式（11）表示客户节点流量平衡约束；式（12）表示每个客户至少由一个配送中心进行服务；式（13）表示每个客户至少由一辆车进行服务；式（14）表示配送中心之间不能互通；式（15）表示每辆车只能从唯一配送中心出发且只有一条配送路径；式（16）表示时间窗约束，配送遵守最晚服务时间窗口；式（17）表示车辆容量限制约束；式（18）表示订单任务点成对约束，即配送员不能在同一地点徘徊；式（19）表示时间连续性约束；式（20）、式（21）表示载重连续性约束；式（22）表示车辆作业结束后必须返回到系统中的一个配送中心；式（23）、式（24）、式（25）为决策变量约束。

$$P_s=\prod _{o\in O,{\beta }_{s,o}=\mathrm{1}}{p}_{s,o}\prod _{o\text{'}\in O,{\beta }_{s,o\text{'}}=\mathrm{0}}(\mathrm{1}-p_{s,o\mathrm{\text{'}}})$$

（8）

$$\sum _{o,o^{\mathrm{\text{'}}}\in O}{x}_{s,o,o\mathrm{\text{'}},v}=\sum _{o,o^{\mathrm{\text{'}}}\in O}{x}_{s,o\mathrm{\text{'}},o,v},\forall s\in S,\forall v\in V,o\ne o^{\mathrm{\text{'}}}$$

（9）

$$\sum _{v\in V}\sum _{i,j\in I}{q}_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}{x}_{s,i,j,v}{y}_{s,i,o}\le Q_{s,o}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}},\forall s\in S,\forall o\in O$$

（10）

$$\sum _{i,j\in I,i\ne j}(x_{s,i,j,v}-x_{s,j,i,v})=\mathrm{0},\forall s\in S,\forall v\in V$$

（11）

$$\sum _{o\in O}{y}_{s,i,o}\ge \mathrm{1},\forall s\in S,\forall i\in I$$

（12）

$$\sum _{i,j\in I,i\ne j}{x}_{s,i,j,v}\ge \mathrm{1},\forall s\in S,\forall v\in V$$

（13）

$$\sum _{o,o\mathrm{\text{'}}\in O,o\ne o^{\mathrm{\text{'}}}}{x}_{s,o,o\mathrm{\text{'}},v}=\mathrm{0},\forall s\in S,\forall v\in V$$

（14）

$$\sum _{o\in O}\sum _{j\in I}{x}_{s,o,j,v}\le \mathrm{1},\forall s\in S,\forall v\in V$$

（15）

$$t_{i,v}\le t_{s,i}^{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}},\forall s\in S,\forall i\in I,\forall v\in V$$

（16）

$$q_{s,i,v}^{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}+q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}\le Q_v^{\mathrm{v}\mathrm{e}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{l}\mathrm{e}},\forall s\in S,\forall i\in I,\forall v\in V$$

（17）

$$\sum _{i\in I}{x}_{s,i,i,v}=\mathrm{0},\forall s\in S,\forall v\in V$$

（18）

$$\begin{array}{l}{t}_{i,v}+t_{i,j,v}^{\mathrm{d}\mathrm{u}\mathrm{r}\mathrm{a}}-M(\mathrm{1}-x_{s,i,j,v})\le {\mathrm{t}}_{j,v},\\ \hspace{1em}\forall s\in S,\forall i,j\in I,\forall v\in V\end{array}$$

（19）

$$\begin{array}{l}{q}_{i,v}^{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}+q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}-M(\mathrm{1}-x_{s,i,j,v})\le q_{j,v}^{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}},\\ \hspace{1em}\forall s\in S,\forall i,j\in I,\forall v\in V\end{array}$$

（20）

$$\begin{array}{l}{q}_{i,v}^{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}+q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}+M(\mathrm{1}-x_{s,i,j,v})\ge q_{j,v}^{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}},\\ \hspace{1em}\forall s\in S,\forall i,j\in I,\forall v\in V\end{array}$$

（21）

$$\sum _{o\in O}{x}_{s,j,o,v}=\mathrm{1},\forall s\in S,\forall j\in I,\forall v\in V$$

（22）

$$x_{s,i,j,v}\in \{\mathrm{0,1}\},\forall s\in S,\forall i,j\in I,\forall v\in V$$

（23）

$$y_{s,i,o}^{}\in \left\{\mathrm{0,1}\right\},\forall s\in S,\forall o\in O,\forall i\in I$$

（24）

$${\beta }_{s,o}^{}\in \{\mathrm{0,1}\},\forall s\in S,\forall \mathrm{o}\in O$$

（25）

式（8）~（25）中：$p_{s,o}$为在情景$s$下配送中心$o$的中断概率；$y_{s,i,o}$表示在情景$s$下客户$i$的需求是否由配送中心$o$供应；$Q_{s,o}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}$为情景$s$下配送中心$o$的储量；$Q_v^{\mathrm{v}\mathrm{e}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{l}\mathrm{e}}$为车辆$v$的最大载重量；$M$为一个比较大的正数；$q_{i,v}^{\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}$为车辆$v$到达客户$i$时的载重量；$t_{i,j,v}^{\mathrm{d}\mathrm{u}\mathrm{r}\mathrm{a}}$为车辆$v$从客户$i$到客户$j$的旅行时间。

式（10）是非线性方程，无法直接使用CPLEX对模型进行精确求解，因此本文引入辅助变量${\Omega }_{s,i,o}$，并将式（10）分解成3个线性化公式进而转换为线性化模型。其中，式（26）、式（27）对${\Omega }_{s,i,o}=q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}{y}_{s,i,o}$进行约束，式（28）是对式（10）的线性化转换。

$${\Omega }_{s,i,o}\ge q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}-M<(\mathrm{1}-y_{s,i,o})$$

（26）

$${\Omega }_{s,i,o}\le q_{s,i}^{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}$$

（27）

$$\sum _{i\in I}{\Omega }_{s,i,o}\le Q_{s,o}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}},\forall s\in S,\forall o\in O$$

（28）

2 算法设计

本文模型为典型的NP-hard问题，因此设计基于GA-SA的混合进化算法求解。传统遗传算法（genetic algorithm， GA）可有效求解现实优化问题，但存在陷入局部最优的风险^［20］；因此一方面混合嵌入大邻域搜索算法（large neighborhood search algorithm， LNSA）的破坏、修复操作，以提高算法搜索效率；另一方面嵌入模拟退火算法（simulated annealing algorithm， SA），以避免遗传算法早熟收敛和算法稳定性差等缺点^［21］，提高算法的全局搜索性能。算法流程见图2，算法步骤详细操作见图3。

图2  混合进化算法流程图

Fig. 2  Flowchart of hybrid evolutionary algorithm

图3  混合进化算法示意图

Fig. 3  Schematic diagram of hybrid evolutionary algorithm

Step 1：构造概率统计模型，利用蒙特卡洛模拟方法构造情景树进行仿真模拟。在t=1阶段，引入确定性参数；在t=2阶段，确定随机变量；在t=3阶段，确定中断概率。以分解形式为例，其中共有$n$个情景，$b(n,\mathrm{3})=\mathrm{3}n/\mathrm{2}+\mathrm{1}$表示第n个情景第3个阶段为$\mathrm{3}n/\mathrm{2}+\mathrm{1}$。

Step 2：对染色体采用二维实数编码，每两个配送中心之间的路径即为一条行驶路线，针对每一条行驶路线指派一辆车对其进行配送。

Step 3：构造初始种群，从客户集合中随机选择一个客户，按照遍历顺序遍历客户集，将客户按照装载量约束、访问顺序依次添加到路径。

Step 4：模型的目标函数为配送成本最低、客户满意度最高，而客户满意度及约束条件违背情况可以用时间及需求惩罚函数表示，因此模型可以转换为求解极小化问题，适应度函数公式为目标函数的倒数。

Step 5：对种群$\mathrm{C}\mathrm{h}\mathrm{r}\mathrm{o}\mathrm{m}$中染色体进行轮盘赌选择、交叉、变异、重插入操作得到新种群。

Step 6：对种群$\mathrm{C}\mathrm{h}\mathrm{r}\mathrm{o}{\mathrm{m}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$进行局部搜索操作，首先，使用破坏算子按照相似性计算公式从当前解移除若干个客户。在式（29）、（30）中，$R_{i,j}$越大，客户$i$与客户$j$之间的相关性越大；其次，使用修复算子按照最小插入成本将被移除的客户重新插回到破坏的解中。

$$R_{i,j}=\frac{\mathrm{1}}{c_{i,j}^{\mathrm{\text{'}}}+V_{i,j}}$$

（29）

$$c_{i,j}^{\mathrm{\text{'}}}=\frac{c_{i,j}}{\mathrm{m}\mathrm{a}\mathrm{x}{c}_{i,j}}$$

（30）

式（29）、（30）中：$c_{i,j}^{\mathrm{\text{'}}}$为将$c_{i,j}$标准化后的值，$c_{i,j}^{\mathrm{\text{'}}}\in [\mathrm{0,1}]$；$c_{i,j}$为$i$与$j$之间的距离；$V_{i,j}$判断$i$与$j$是否在同一条路径上，即是否由同一辆车服务，如果$i$与$j$在同一条路径上则$V_{i,j}=\mathrm{0}$，否则为1。

Step 7：设计模拟退火接受准则，假设$\mathrm{C}\mathrm{h}\mathrm{r}\mathrm{o}\mathrm{m}$的解为$S^{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}$，$\mathrm{C}\mathrm{h}\mathrm{r}\mathrm{o}{\mathrm{m}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$的解为$S^{\mathrm{n}\mathrm{e}\mathrm{w}}$，如果$S^{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}$优于$S^{\mathrm{n}\mathrm{e}\mathrm{w}}$，则接受概率为1，如果$S^{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}}$劣于$S^{new}$，则接受概率为式（31），其中$T$表示当前温度，最终产生新种群$\mathrm{C}\mathrm{h}\mathrm{r}\mathrm{o}{\mathrm{m}}^{\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}$。

$$P=\left\{\begin{array}{l}\mathrm{1},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}f(S^{\mathrm{n}\mathrm{e}\mathrm{w}})<f(S^{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}})\\ {\mathrm{e}}^{-[f(S^{new})-f(S^{current})]/T},f(S^{\mathrm{n}\mathrm{e}\mathrm{w}})\ge f(S^{\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}})\end{array}\right.$$

（31）

Step 8：进行退火操作，随着迭代次数的增加，接受概率需要随之降低并更新迭代次数gen=gen+1，即温度T随着迭代次数的增加而衰减，退火公式如下，其中gen为当前迭代次数。

$$T^{\mathrm{g}\mathrm{e}\mathrm{n}+\mathrm{1}}=\alpha T^{\mathrm{g}\mathrm{e}\mathrm{n}}$$

（32）

Step 9：如果$\mathrm{g}\mathrm{e}\mathrm{n}\le \mathrm{m}\mathrm{a}\mathrm{x}\mathrm{g}\mathrm{e}\mathrm{n}$，则继续迭代，跳至Step 4；否则输出最优解，算法结束。

3 算例分析

3.1　算例描述

在模型算例分析中，相关参数数值的选择和估算是其中关键一环。客户需求数据包括需求量、配送时间窗、服务时间要求等，将直接影响配送效率和客户满意度；车辆系数和配送中心系数等数据可描述配送资源的限制条件和费用。数据的合理选择可以更准确地反映配送成本与客户满意度，因此本文参考历史订单记录、市场调研情况、地理信息等，根据实际进行相应的等比例修正和调整。例如，以药品配送为例，药店作为配送中心，储量一般相差不大，可满足10~20个客户的总需求。以2022年上海市新冠肺炎疫情期间管控为例，假设存在某一开放式配送系统，有4个配送中心、20个客户需求，配送中心编号为1、2、3、4，客户编号为5、6、7、…、24。其中，客户需求量、时间窗、配送中心物资储备量是不确定的，配送中心信息、客户需求信息、车辆信息如表1~3所示。

表1  配送中心信息

Tab. 1  Information of distribution center

编号	坐标	储量	时间窗/min	服务时间/min
1	（41，70）	210	［0，1 000］	0
2	（73，65）	220	［0，1 000］	0
3	（28，29）	200	［0，1 000］	0
4	（69，33）	200	［0，1 000］	0

表2  客户需求信息

Tab. 2  Information of customer demand

编号	坐标	需求量	时间窗/min	服务时间/min
5	（42，66）	15	［12，100］	10
6	（35，45）	10	［35，67］	20
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$

表3  车辆系数

Tab. 3  Coefficient of vehicles

车型编号	载货量/veh-1	租赁成本/（元·veh-1）	点位费/元
1	50	90	6
2	100	150	4
3	150	210	3

情景树的概率设计如图4所示，情景集S共包括25个情景：在t=1阶段，输入确定性参数（例如位置信息、车辆系数、成本信息等）；在t=2阶段，输入确定性情景的参数，其对应的情景发生概率0.3在同阶段最高，基于确定性情景按照比例生成离散的不确定性参数组合$(Q_{s,o}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}},q_{s,i}^{\mathrm{d}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}},t_{s,i}^{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}},t_{s,i}^{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}})$，作为不确定性情景输入参数，并赋予不同的情景发生概率；在t=3阶段，根据参数组合确定配送中心中断概率，例如需求量大而时间窗要求严格的客户处于疫区的概率较高，因此其对应情景的配送中心失效概率相对较高，需求量小而时间窗要求宽松的客户则往往具备较为良好的配送环境，对应的配送中心不易失效。在情景集中，情景1‒2‒7代表确定情况下的城市即时配送环境，对应的情景发生概率为$P_{\mathrm{1}}=\mathrm{0.3}\times \mathrm{0.3}=\mathrm{0.09}$；情景1‒6‒31代表极端封控情况下的城市即时配送环境，对应的情景发生概率为$P_{\mathrm{25}}=\mathrm{0.1}\times \mathrm{0.3}=\mathrm{0.03}$。其余介于二者之间。

图4  案例情景树概率设置

Fig. 4  Probability setting of case scenario tree

3.2　算例结果分析

在随机变量$(Q_{s,o}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}},q_{s,i}^{\mathrm{d}\mathrm{e}\mathrm{m}\mathrm{a}\mathrm{n}\mathrm{d}},t_{s,i}^{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}},t_{s,i}^{\mathrm{m}\mathrm{a}\mathrm{x}\_\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{w}})$确定且设施未发生中断的情况，即确定性情景下，总成本为1 857.931 8元；一旦设施发生中断，总成本达到 2 552.431 3元；在考虑参数随机性及中断可能性的再分配情景下，总成本为2 049.262 0元。3种情景下的案例求解结果见表4，配送路线见表5，配送方案示意见图5。

表4  案例求解结果

Tab. 4  Determination of solution results for cases under given scenarios

情景	可用配送中心编号	总成本 /元	车辆数/veh	配送成本 /元	总行驶距离/102 m	需求惩罚成本 /元	时间惩罚 成本/元	求解时间/s
正常	1、2、3、4	1 857.931 8	10	1 857.9318	515.931 8	0	0	19.993 8
中断	1、2、4	2 552.431 3	8	1 531.601 4	469.601 4	940.000 0	0
再分配	1、2、4	2 049.264 0	9	1 855.261 0	617.260 8	184.003 0	10.000 0	22.2347

表5  不确定情景下行驶路线

Tab. 5  Driving route under uncertain situations

编号	车型	路径	预计完成时间/min	成本/元
1	2	1→11→1	53.266 6	151.266 6
2	2	1→6→9→5→1	121.949 2	246.949 2
3	2	1→16→15→17→4	167.578 1	237.578 1
4	2	1→10→12→13→18→4	97.937 0	231.937 0
5	3	2→7→8→2	76.068 0	170.068 0
6	2	2→14→2	46.832 8	134.832 8
7	1	4→21→19→20→4	148.157 4	296.629 6
8	1	4→22→23→4	91.682 2	222.682 2
9	2	4→24→4	65.317 3	163.317 3

图5  配送方案示意

Fig. 5  Distribution scheme

根据表5行驶路线，并结合图6配送方案可以发现，一旦3号配送中心失效，一方面导致客户点15、16、19、20无设施提供服务，虽减少$\mathrm{157.084}\mathrm{7}\mathrm{元}+\mathrm{223.586}\mathrm{7}\mathrm{元}=\mathrm{380.671}\mathrm{4}$元的配送成本，但产生940.000 0元的需求惩罚成本，导致总成本增加559.328 6元；另一方面部分车辆完成服务后不再能选择距离最近的3号配送中心。最终，3号配送设施的中断导致总成本增加37.38%。

图6  客户需求量与配送中心储量三维曲面图

Fig. 6  3D surface of customer demand and distribution center reserves

故不论是出于人道主义物流还是经济效益，均有必要对3号配送中心的客户进行再分配，将15、16、17号客户分配给1号配送中心，将19、20、21、22、23号客户分配给4号配送中心，整合路线1‒10‒12‒（3）与1‒10‒13‒18‒4为新路线1‒10‒12‒13‒18‒4。中断后的再分配使得总成本增加10.30%、配送成本降低0.14%、配送距离增加19.64%、惩罚成本从0元增加到194.003 0元。可发现在大规模公共卫生事件背景下，对中断设施服务的再分配可有效降低成本27.08%，经再分配后的惩罚成本增加幅度较大，而配送成本反而降低。原因在于3号配送中心失效后，其服务的客户需要再分配给其他正常运行的配送中心，导致配送距离增加幅度较大，在同客户满意度的权衡下，再分配后的配送方案路线可能具有较低的配送成本（包括车辆租赁成本、车辆服务成本、燃油成本）和较高的惩罚成本，以达到尽可能低的总成本。

3.3　灵敏度分析

为进一步验证城市即时配送模型的有效性和适用性，从配送中心储量、客户需求量、客户单侧时间窗、客户最大可接受服务时间、客户满意度权重及中断概率六个方面分析总成本的波动情况。综合考虑研究目的、实际情况以及模型的稳定性、可靠性和合理性等因素，选择最具代表性的参数涨跌幅度设置特定的步长，模拟实际的业务变化情况。例如，配送中心储量的增加幅度为25%，因为配送中心通常会有一定的储备量，以应对业务峰值时的需求；客户需求量的增加幅度为20%，因为大规模突发事件下，客户需求的增加量一般会比平时更大。不同参数对优化结果的影响，见图6~11。

图7  客户需求量与配送中心储量等高线图

Fig. 7  Contour plot of customer demand and distribution center reserves

图8  单侧时间窗与最晚可接受时间三维曲面图

Fig. 8  3D surface diagram of unilateral time window and last acceptable time

图9  单侧时间窗与最晚可接受时间等高线图

Fig. 9  Contour plot of unilateral time window and latest acceptable time

图10  经济成本与惩罚成本的Tradeoff曲线

Fig. 10  Tradeoff curve of economic and penalty cost

图11  中断概率分析

Fig. 11  Analysis of interrupt probability

分析图6、7，随着配送中心储量不断增加，总成本呈现先减后增的变化趋势。当配送中心储量由100增加至250时，配送方案得到更优安排，总成本降低；当配送中心储量从250增加到450时，配送方案的优化空间有限，而配送中心规模的增大使配送中心运营成本增加。随着客户需求量的增加^［22］，总成本在初期可通过灵活选择车型来应对需求的波动，增速相对平缓；当需求量增幅从100%增加到200%时，供需失衡导致惩罚成本增加。此外，配送中心储量与客户需求量按照同趋势变化时，总成本在同等水平下相对较低，此时既满足了客户需求又选择了合适的配送中心储量，避免惩罚成本与过高的配送中心运行成本。

分析图8、9，随着客户单侧时间窗的增加，总成本先减少后趋于稳定。当单侧时间窗从50%提高到100%时，配送方案可进行更优的安排而不违背时间窗要求；当单侧时间窗从100%提高到200%时，此时配送方案已达到相对最优，客户时间窗对配送方案影响较小。最晚可接受服务时间对总成本影响不大，原因在于过高的惩罚系数导致配送方案最大程度避免延误。此外，图中存在部分区域总成本较高，此时左时间窗、单侧时间窗、最晚可接受服务时间接近，对客户满意度的鲁棒性较低，导致总成本较高；反之，部分区域总成本较低的主要原因在于多个时间窗口间存在较大的可调整范围，配送系统有更多可以选择的方案且不会违背客户的时间窗要求，使得总成本降低。

分析图10，分别设置客户满意度权重$\alpha$=0，0.05，…，1.00，对模型进行求解并将配送总成本、客户满意度分别归一化处理，使用多项式拟合曲线表达经济成本与惩罚成本（反映客户满意度水平）的Tradeoff关系。当客户满意度权重$\alpha$从0增至1，经济成本从0增至1，惩罚成本减少66.240 8%，95%置信水平条件下惩罚成本与经济成本的曲线为$f(x)=\mathrm{0.662}\mathrm{4}x+\mathrm{0.673}\mathrm{3}$，即客户满意度会在牺牲部分经济成本的条件下逐渐增加。

分析图11，随着中断概率的增加，正常运行的配送中心数量减少，原先相对最优的配送方案无法继续执行，总成本、配送成本、未满足需求量呈现增加的趋势；当失效配送中心数量达到一定程度时，此时配送方案仅对部分客户进行配送任务，配送系统需要服务的需求数量减少，无法满足的需求数量增加，导致配送成本降低但总成本增加。

3.4　算法有效性分析

为验证算法有效性，基于配备Intel（R） Core（TM） i5‒8300H处理器、搭载8GbCPU运行内存的Win10系统，以确定性情景作为求解算例，与商用求解软件CPLEX 22.1.1以及传统遗传算法进行最优解性能、求解时间的对比，结果见表6。

表6  基准算例对比

Tab. 6  Comparison of reference examples

算例规模	算法	总成本/元	求解时间/s
小	CPLEX	1 490.278 9	4.480 0
	遗传算法	1 644.291 9	14.674 2
	混合进化算法	1 526.8015	19.993 8
大	CPLEX
	遗传算法	12 903.925 0	25.351 2
	混合进化算法	10 556.887 0	31.491 3

针对小规模算例，混合进化算法相较于CPLEX的解性能差2.45%，求解时间多15.51 s；相较于遗传算法，混合进化算法的解性能高4.70%，求解时间多5.32 s；而针对大规模算例，CPLEX难以在可接受时间（7 200 s）内得到可行解，而混合进化算法仅比遗传算法多6.14 s的求解时间，解性能差异却达22.23%。显然，混合进化算法的解性能完全可以接受，且具有较快的求解速度。考虑到实际问题的复杂性，混合进化算法可以在求解效率和求解质量之间取得较好的平衡，为实际应用提供可行的解决方案，这也是美团等企业惯用的做法。

此外，在具体方案上，车辆使用数目与总成本可能呈现反比关系，原因在于，虽然增加车辆数目可能会导致更高的租赁成本、更长的行驶距离，但该解决方案一方面可以更好地考虑车型、车辆的访问顺序等因素，从而更好地满足客户需求，提高客户满意度；另一方面，针对不同客户需求的车辆车型使用情况也存在区别，若运往同一批客户的总需求量为400，所需访问的客户点数量为15，则情景a使用4辆2号车的配送成本为$\mathrm{4}\times \mathrm{150}+\mathrm{15}\times \mathrm{4}=\mathrm{660}$，情景b使用3辆3号车的配送成本为$\mathrm{3}\times \mathrm{210}+\mathrm{15}\times \mathrm{3}=\mathrm{665}$，显然，使用更多车辆的情景a反而具有更低的配送成本，而类似情况在现实情景是真实存在的。故在某些配送情景下，即使使用较多的车辆，也可获得更高的经济效益。

4 结语

本文针对大规模公共卫生事件背景下的城市即时配送问题，构建考虑道路限行的城市道路交通网络，建立多目标随机规划模型，使用蒙特卡洛模拟方法生成情景树，设计基于GA‒SA的混合进化算法进行求解，以2022年上海新冠肺炎疫情事件为例进行可行性验证，发现：大规模公共卫生事件一方面导致配送设施中断概率增加及储量降低，另一方面导致客户需求具有更多不确定性，配送系统容易因供需不匹配而发生“爆单”“爆仓”，虽然配送成本未必增加，但客户满意度降低导致总成本增加。因此一方面要提高城市即时配送系统的鲁棒性，提前做好设施服务中断的准备，储备冗余配送中心与运力，避免“爆单”“爆仓”问题发生；另一方面要加强对居民恐慌心理的疏导，提高居民在疫情期间的安全感，减少物资需求的不确定性。

作者贡献声明

孟令鹏：模型构建，算法设计，结果分析，论文撰写。

王旭东：文献梳理，算法实现，结果分析，论文撰写。

韩传峰：结果分析，论文撰写和修改。