The stochastic location-routing problem with parallel truck–drone operations for humanitarian aid delivery

随机选址-路由问题：基于卡车-无人机并行操作的救援物资配送。

信息：
作者：Hannan Tureci-Isik, Melih Çelik , Ece Sanci
单位：University of Bath, Amazon Logistics
期刊：European Journal of Operational Research
日期：2025-09-02

---

1. 概述

1.1. 背景：

不确定性 + depot 选址 + truck–drone 并行：

灾后救援场景下，物流系统的核心目标并非传统意义上的“运输成本最小化（路径最短）”，而是尽可能早地将救援物资送达所有需求点（到达延迟）。但在地震、洪水等突发灾害发生后，道路网络往往遭受不同程度的破坏，使得卡车的行驶时间呈现出显著的不确定性。在这种条件下，仅依赖地面车辆进行配送往往难以保证救援时效。

无人机由于能够绕开受损道路，在灾害救援中展现出独特优势，但又受到续航能力、飞行距离和载重等因素的严格限制。因此，更好的救援模式是卡车与无人机协同、并行执行任务：卡车负责在地面网络中服务部分需求点，无人机则从救援中心起飞，对部分可达节点进行快速点对点配送。

在这一背景下，灾前规划阶段就需要确定哪些地点作为救援中心（depots），并配置相应的卡车与无人机资源，但这一决策必须在灾害尚未发生、道路状态未知的情况下做出，因此不可避免地面临灾后道路条件的不确定性。

1.2. 贡献：模型与算法各占一半

1. 新问题定义 + 新两阶段随机规划模型：提出了一类新的随机 location-routing 问题，将“灾前 depot 选址决策（第一阶段）”与“灾后卡车路径规划与无人机服务顺序决策（第二阶段）”统一纳入一个两阶段随机规划框架中，并以平均到达时间（latency）而非运输成本作为唯一优化目标（与传统的LRP是多目标优化不同）。

2. 面向稀疏路网与重复访点的流式建模方法：针对灾后道路网络稀疏、非三角不等式成立以及卡车可能多次经过同一节点的特点，作者放弃了传统 LRP 中基于路径变量与子回路消除的建模方式，转而采用 commodity flow（商品流） 表达卡车服务进度，并通过一组辅助变量精确刻画“首次离开、最后一次离开与非首次访问”等路径结构特征。

3. 可扩展求解：VNS 元启发式：由于所提出的随机混合整数模型在实际规模下难以由商业求解器直接求解，本文设计了一种面向两阶段结构的变邻域搜索（VNS）算法。该算法在外层对 depot 选址进行搜索，在内层对给定 depot 集合下的卡车–无人机协同路由进行跨场景评估与改进，从而在解质量与计算效率之间取得平衡。

---

2. 研究方法

2.1. 问题定义

2.1.1. 两阶段随机 location-routing with parallel truck–drone

• 节点集合：候选 depot 节点集合 $W$ 与需求节点集合 $D$，两者统称为图上的节点 $N$。论文用无向图 $G=(N,E)$ 表示路网, $E$ 是真实存在的道路边
• 不确定性建模：用离散场景 $\xi \in \Omega$ 描述灾后道路状态，在每个场景下，卡车在边 $(i,j)$ 上的旅行时间表示为 $t_{ij}(\xi )$。用“行驶时间膨胀”来隐式表达道路受损程度
• 并行作业假设：每个被选中的 depot 配置一辆卡车和一架无人机。卡车执行一条地面服务路径（tour），无人机可从 depot 多次起飞，每次服务一个需求点并返回。卡车与无人机在时间上并行执行、相互独立。（即每个需求点要么被卡车服务，要么被无人机服务）。
• 无人机可达性（续航约束）：$W_i$为无人机能在续航内从 depot $k$ 往返服务需求点 $i$ 的 depot 集合；同时，对于每个给定的 depot $k$，将其可服务的需求点按往返飞行时间非降序排序形成集合 $D_k$。(作者指出：要最小化无人机到达时间的平均值，访问顺序应按往返时间从小到大) 。这一排序结构在后续用于刻画无人机的到达时间，本质上等价于一个单机调度中的累计完工时间模型。

注：唯一被建模为随机变量的是：卡车在道路上的行驶时间。路网结构本身（哪些边存在）是确定的、需求点位置、需求规模是确定的、无人机的飞行时间是确定的、只有卡车在边 $(i,j)$ 上的行驶时间 $t_{ij}(\xi )$ 随灾害场景 $\xi$ 变化。

2.1.2. 目标函数：期望平均到达时间（latency）

二阶段目标 $g(\xi )$ 由两部分组成：$$\min \frac{{\sum }_{(i,j)\in E}{f}_{ij}(\xi )t_{ij}(\xi )-{\sum }_{k\in W}{\sum }_{(i,j)\in E}{y}_{ijk}(\xi )t_{ij}(\xi )+{\sum }_{i\in D}{p}_i(\xi )}{|D|}$$

1. 卡车服务的到达时间总和
   $$\sum _{(i,j)\in E}{f}_{ij}(\xi )t_{ij}(\xi )-\sum _{k\in W}\sum _{(i,j)\in E}{y}_{ijk}(\xi )t_{ij}(\xi )$$

   其中 $f_{ij}(\xi )$ 是商品流，$y_{ijk}(\xi )$ 是“卡车（从 depot k 出发）是否走边(i,j)”的二进制变量。$\text{商品流}\times \text{边时间}$ - 校准项$\sum yt$
   校准项：由于商品流 $f_{ij}$ 在“离开 depot 的第一条边”上，会多算一次（因为 depot 本身不是需求点）同时，卡车在完成最后一个服务点后，仍可能需要走一段路到虚拟汇点 $s$ ，但这段路不应该再贡献任何到达时间。所以要 $ -\sum y_{ijk} t_{ij}$

2. 无人机服务的到达时间总和

   无人机部分则通过“从 depot 出发访问按飞行时间排序的点”来表达其到达时间 $p_i(\xi )$，并直接以$\sum p_i(\xi )$ 纳入目标函数。

2.1.3. 二阶段随机 MIP

第一阶段（灾前）：depot 选址

• $x_k\in 0,1$：表示候选点 $k$ 是否被选为 depot，并满足 depot 数量上限约束 ${\sum }_k{x}_k\le n$。

第二阶段（灾后，按场景 $\xi$）：分配、路径与服务顺序决策

• 卡车边变量：$y_{ijk}(\xi )\in 0,1$: 场景 $\xi$ 下，来自 depot k 的卡车是否走边 $(i,j)$
• 服务分配变量：
  • c i k ( ξ ) ∈ { 0 , 1 } c_{ik}(\xi)\in\{0,1\} cik​(ξ)∈{0,1} 表示需求点 $i$ 是否由 depot $k$ 的卡车服务
  • z i k ( ξ ) ∈ { 0 , 1 } z_{ik}(\xi)\in\{0,1\} zik​(ξ)∈{0,1} 表示需求点 $i$ 是否由 depot $k$ 的无人机服务（仅对 $k\in W_i$定义）
  • 每个需求点必须且仅能由一种方式服务：${\sum }_k{c}_{ik}(\xi )+{\sum }_{k\in W_i}{z}_{ik}(\xi )=1$。
• 无人机到达时间：$p_i(\xi )\ge 0$
• 商品流：$f_{ij}(\xi )$ 在扩展网络 $E^{\prime }$ 上流动。作者引入一个虚拟汇点 $s$，使卡车路径在完成最后一次服务后直接终止于 $s$，从而避免强制返回 depot，因为回程不会对到达时间目标产生正面作用。

2.1.3.1. 随机场景的产生

1. 在每个场景 $\xi$ 中随机选一个点作为灾害中心（在实验里是均匀分布在区域内）。
2. 然后对图上的每一个节点（以及每一条边），计算它到灾害中心的距离，然后按距离落入三个同心区域之一：破坏区、损伤区、强影响区
3. 对每一条边，根据它“最严重的端点”决定影响等级：对每一条边$(i,j)$，看节点 $i$ 和 $j$ 分别落在哪个区域，取影响更严重的那个区域
4. 用“随机膨胀系数”生成行驶时间：对每条边$(i,j)$，在场景 $\xi$ 下 $t_{ij}(\xi )=t_{ij}^{\text{base}}\times {\lambda }_{ij}(\xi )$，其中 ${\lambda }_{ij}(\xi )$ 是膨胀系数，来自与影响区对应的区间均匀分布。

也就是说这是“连续随机性”，不是“结构性不确定性”，所以让模型始终“可行”，没有考虑“断路”这种极端情况。

2.1.3.2. 商品流 commodity flow

在本文中，优化目标是所有需求点的平均到达时间。直观上，最自然的建模方式是：给每个需求点一个“到达时刻变量”，然后把它们加起来。但这样做会立刻遇到两个致命问题：

• 到达时刻取决于路径顺序，而路径顺序在稀疏网络中可能包含回退和重复访问；
• 在随机场景下，每个需求点、每个场景都要一个时间变量，模型规模会爆炸。

所以作者采用了一种不显式引入到达时间变量的建模技巧：用商品流（commodity flow）来“隐式累计”到达时间。 commodity flow 代表“当前卡车上仍然对应多少尚未完成服务的需求点”。也就是说：刚从 depot 出发时，车上“背着”所有还没被服务的需求点；每当卡车第一次到达并服务一个需求点，这个“剩余数量”就减 1；如果之后只是“路过”该节点（重复访问），这个数不减少。

本文的问题并不强制每个节点只能访问一次： “访问（visit）”与“服务（service）”被明确区分。卡车在执行救援任务过程中，可能由于灾后路网稀疏或受损，需要多次经过同一节点；但每一个需求点只在其第一次被满足需求时才被视为完成服务，并对目标函数产生贡献。因此，模型只要求每个需求点恰好被服务一次，而不要求其在路径中仅被访问一次。该结构的作用可以通过图1中的示例直观理解：

有 3 个需求点：2, 3, 4；1为 depot
卡车路径（允许重复访问）为：1 → 2 → 3 → 2 → 4 → s（虚拟终点）

1 → 2 初始流量 $f_{1,2}(\xi )=4$ (还剩3个尚未完成服务的需求,再加上一个虚拟汇点s)，所以这条边的行驶时间，会被乘以 4 计入目标；以此类推…第一次到达节点 3，剩余需求点数 = 2 (4, s) ；3 → 2 这条边的行驶时间，会被乘以 2。此时再次到达节点 2（只是路过，不是服务），剩余需求点数仍然是 2，所以2 → 4 这条边的时间，仍然乘以 2。

---

2.2. 求解方法：两阶段结构上的 VNS 元启发式

作者指出，即使在较小规模实例下，所提出的随机混合整数模型也难以由商业求解器在合理时间内求得最优解，因此设计了一种面向问题结构的可变邻域搜索（VNS）算法。

2.2.1. 算法1：StochasticLocationRouting

本质上是外层只做“换 depot”的局部搜索

总体算法从一个随机 depot 集合开始。定义 $N_{\text{depot}}$ 为swap 邻域：关闭一个已开 depot，同时打开一个未开的 depot。并以期望 latency 作为评价指标，逐步迭代改进 depot 选址方案。

算法1自己不会碰“路线细节”，它只会反复调用算法2去评估“这个 depot 方案好不好”。

2.2.2. 算法2：DepotEvaluation

评估一个 depot 集合在所有场景的期望 latency

输入：固定开设的 depot 集合 $X$ 以及所有场景集合 $\Omega$;
对每个场景 $\xi \in \Omega$:

1. 先构造一个初始解 $R(X,\xi )$ （算法3）
2. 再对这个初始解做改进 $Improvement(\cdot )$ （算法5）最后返回期望目标 $E_{\xi |X}[h(R(X,\xi ))]$

相当于对同一个 depot 方案，要在很多种灾后路况下都跑一遍路线，并把平均表现当作这个 depot 方案的分数。

2.2.3. 算法3：InitialConstruction

算法3先做一件很粗但稳定的事：把每个需求点分配给“地理距离最近”的开设 depot

truck tours：然后对每个 depot $k\in X$ 单独造卡车 tour：从depot k 出发，贪心地加最近的可达节点（在场景 $\xi$ 下，用旅行时间 $t_{ij}(\xi )$ 做nearest-neighbour 插入）

drone schedules：算法3把“被分给无人机的点”放进 $D{L}_k$ ，然后按照某个顺序形成无人机访问序列 $D{S}_k$ ，并且它会把“可行的（续航内）”优先放，若出现不可行分配也会先保留，交给后面 算法5 VNS 处理

难点：如果“当前点”到任何未访问点都没有直接边怎么办？ – 引入 truck backtracking
当发现没有直接路可走时，卡车必须“沿着走过的路往回退”，找一个历史节点作为“回退点 BackNode”，从那里能接上一条边去到未访问点。如果一路回退都找不到这样的点（dead-end），算法3干脆承认卡车走不通了：把剩余所有未访问点都扔给相同 depot 的无人机去服务

回退点例子：

当前卡车 tour 是1 → 2 → 3 → 4 → 3 → 5 → 6，未访问的只有节点 7 ：

• 从当前节点 6 出发，发现 6 和 7 没有边（走不过去）
• 于是沿 tour 倒着检查：6 前面是 5、3、4、3、2……
• 最终发现 节点2 与 7 有边，所以回退点是 2
• 回退路径不是傻傻按原路倒回，而是用算法4找“最短回退路”，最终形成 tour：
  1 → 2 → 3 → 4 → 3 → 5 → 6 → 3 → 2 → 7

2.2.4. 算法4：Find shortest backtracking route

算法4在需要回退时被调用：

• Step 1：截取 partial tour：从回退点（2）到当前点（6），把这段原 tour 片段抽出来 （2 → 3 → 4 → 3 → 5 → 6）
• Step 2：简化 partial tour：如果这段里有重复访问的节点，就删掉“从第一次到最后一次访问之间”的冗余段 （4），得到一个更短、更干净的“骨架路径” （ 2→3→5 + 当前点6）
• Step 3：找共享边/更短连接：简化后是 2→3→5，然后从当前点 6 去找能接上的点：6 与 2 不连，但 6 与 3 连，所以先从 6 走到 3，再退到 2，再去 7.
• 最终形成 tour：1 → 2 → 3 → 4 → 3 → 5 → 6 → 3 → 2 → 7

2.2.5. 算法5：VNS

当算法3给出初始的 truck tours + drone schedules 之后，就进入算法5的 VNS 改进阶段

VNS 的邻域（N1~N7）：

• N1：truck→truck 重新分配：把一个需求点从一条卡车 tour 挪到另一条卡车 tour
• N2：truck→drone（同 depot）：从卡车 tour 拿掉一个点，改由同 depot 的无人机去送
• N3：truck→drone（不同 depot）：从卡车 tour 拿掉一个点，改由另一个 depot 的无人机送
• N4：drone→truck（同 depot）：把无人机送的点改回同 depot 的卡车来跑
• N5：drone→truck（不同 depot）：把无人机送的点改给别的 depot 的卡车
• N6：drone→drone：在不同无人机之间转移一个需求点
• N7：重建一条卡车 tour：随机选一辆卡车，改变其“第一个访问的需求点”并重建整条路线；原来那条 tour 上的点先临时分配给最近的无人机（一个很强的 shaking，N7 只用于 shaking，不用于 local search）

VNS 的流程：先 shake 再 local search（并且 best-improvement）

• 每次迭代先做 shake：随机从当前解的某个邻域里跳到一个邻居解（用N1-N7中的一个）
• 然后做 local search：用 best-improvement，在同一个邻域里反复找改进，找不到才换下一个邻域；如果找到改进就回到第一个邻域重新开始
• 如果 local search 没带来更好解，就换一个更大的 shaking 邻域；直到所有 shaking 邻域都试过

也就是说整个流程是：

• 算法1：
  随机选 $n$ 个候选点作为初始 depot 集合 $X_0$（灾前决策的起点）。

• 算法2–4：
  在固定 depot 集合 $X$ 和固定场景 $\xi$ 下，构造一个“能跑”的初始解（卡车路径 + 无人机分配 + 回退）。

• 算法5（VNS）：在不改变 depot 的前提下，反复改进

  • 卡车 tour
  • 卡车 / 无人机分配
  • 需求点归属哪个 depot，以降低该场景下的 latency。

• 算法1 的外层搜索：用算法2–5算出来的“期望 latency”作为评分，判断是否要替换 depot 集合。

---

3. 实验

3.1. 实验设置

• 数学模型用 Gurobi 11.0.3 + Python API；
• 启发式用 C++（Visual Studio 2022）；
• 机器：Windows 10, i7-10510U, 16GB RAM。

3.2. 小规模随机实例：和 Gurobi 的对照实验

作者生成了 10 个小实例：15 个需求点、5 个候选 depot，开 2 个 depot，考虑 10 个灾害场景。节点随机放在 $5000\times 5000$ 平面上；道路网络是稀疏的随机连接。卡车/无人机速度相同；为了简化，先忽略无人机航程限制。

对Gurobi施加了3小时的计算时间限制。且由于求解两阶段随机混合整数规划模型的计算量较大，Gurobi无法在时间限制内缩小MIP差距。为了得到最优解，通过列举所有可能的开放仓库集来运行Gurobi，并对每个组合求解，其中期望目标值最小的开放仓库集作为最优解。

表3显示了作者提出的启发式算法在所有场景下都能找到接近最优的卡车路线和无人机调度，且速度快。

邻域（N1–N7）谁更关键（消融实验）

表4说明了 N7 和 N4 最关键：拿掉 N7，目标平均上升 5.66%；拿掉 N4，上升 3.24%。

3.3. 真实案例：2011 年土耳其 Van 地震

• 96个需求点，15个 depot 候选点。

• 卡车运行时间通过谷歌地图得到。

• 使用25 mph(约40 km/h)作为无人机速度，30 min作为无人机续航时间。

• 使用Python中的“geodesic”库计算地点之间的测地线距离，并将距离除以无人机速度，得到无人机飞行时间。

• 输入是需求点坐标与灾害场景，输出包括：

  • 开哪些 depot（以及每个 depot 配置 1 truck + 1 drone 的并行系统）；
  • 每个场景下卡车路线、无人机服务顺序；
  • 最终用“跨场景的期望平均到达时间”作为指标，同时也报一些“最大到达时间”等辅助指标。

3.3.1. 基线策略对比

三组基准对照：

• Truck-only（只有卡车）：平均 gap 高达 116.8%；作者把它解释成：只用卡车导致所有点都得靠地面不确定路网送达，延迟大幅上升。（图5b）
• No road damage assessment（不做路损评估）：只用“灾前旅行时间”的单场景来选址，然后再把选址拿去真实多场景里评估；平均会更差，文中说这会带来大约 7.4% 的期望目标下降/改进（换句话说：不评估路损会明显变差），并用图5c展示 instance 10 下 depot 覆盖差异。（图5d）
• K-means 选址（先按地理聚类定 depot）：平均目标恶化 19.1%，并且在某些实例下会直接不可行（例如无人机续航约束下 infeasible）。图5d展示了 instance 10 下 K-means 导致 depot 分布更“均匀”，但对“灾害断裂带附近”的需求反而覆盖不对。

3.3.2. 敏感性分析

1. 开 depot 数量的边际收益递减：当开放 depot 从 3 增到 15，期望平均到达时间从约 3.4h 降到 0.4h，且每加一个 depot 都能改善但收益递减（隐含“投资成本 vs 响应时间”的权衡）。见图6a。

2. 无人机续航：从 30min 到 1h 收益大，之后趋于饱和：把续航从 30min 提到 1h，平均到达时间减少约 18min（约 20%）；从 1h 再增加 1h 仅约 1% 改善；甚至“无限续航”也不比 2h 更好（饱和）。这实际上给了一个很强的工程结论：现有技术的无人机已经能带来主要收益，继续堆续航不是最优投资方向。见图6b。

3. 无人机速度：提升显著但仍有递减收益：例如速度到 56km/h 相比基线 40km/h，平均到达时间减少约 24min（约 26%）。见图6c。

3.4. 不确定性的“价值量化”：VSS 与 EVPI

作者用两个经典指标把“考虑不确定性到底值不值”讲清楚：

• VSS（Value of Stochastic Solution）：比较“真正的随机模型（考虑多场景）”与“用期望旅行时间替代不确定性、当作单场景来做”的差距。他们的计算方式是：先用期望 travel time 求 depot 位置，再在所有场景下生成路由评估其期望目标；与随机模型的期望目标差值就是 VSS。文中给出：平均 %VSS 为 3.4%，但最大可达 17.4%；作者把它换算成总响应时间，最大情形可改善 24.5 小时量级，强调在救灾语境下非常关键。

• EVPI（Expected Value of Perfect Information）：“如果你提前知道真实场景，会最多好多少”。做法是对每个场景当作已知，分别跑启发式求最优决策并再取期望，得到 $z^{PI}$；与随机模型解 $z^h$ 的差就是 EVPI。

这两段的逻辑价值在于：它把“随机建模不是学术花活”变成一个可量化的管理结论：哪怕平均看起来只有几个百分点，换成‘小时级生命线’的救援时效就是硬收益。

---

4. 总结

4.1. 结论

1. 卡车+无人机并行能显著改善灾后到达时效，尤其在道路受损不确定的情形下，无人机对“被道路隔离”的节点提供了关键补充。
2. 灾前 depot 选址必须把灾害不确定性纳入决策：VSS 的量化说明“只用期望路况做确定性规划”会在一些实例里导致显著的时效损失。
3. 工程参数上，提升无人机续航/速度与增设 depot 的收益存在明确的递减区间，从而能为资源投入提供更实用的边界判断。

4.2. 限制

这篇论文的限制主要来自建模的简化：

• 每个 depot 固定 1 truck + 1 drone，未考虑多车多机规模配置、或车辆在 depots 间调度。
• 无人机每次起飞只服务一个点，且忽略装载/换电等待时间；这使得模型更聚焦“到达时效”，但弱化了更细的操作约束。
• 需求被假设为“每个点一趟即可满足”，从而把重点放在时效与可达性，而不是多趟补给或容量约束。
• 不确定性只体现在“道路 travel time 膨胀”，未显式建模“道路彻底不可达（binary failure）”、天气导致无人机不可飞等更复杂联动。

4.3. 未来方向

基于上述限制，合理的未来方向通常会落在：

• 更丰富的 fleet 维度：每个 depot 多车多机、或者车辆共享与跨 depot 调度；
• 更精细的无人机作业模型：允许一次飞行访问多个点、显式建模换电/充电/排队；
• 更广的随机性/鲁棒性：把“不可达边/道路失效概率”“风雨导致无人机停飞”等纳入场景，或引入风险度量。