基于分布鲁棒优化的无人机配送设施选址分配问题

Distributionally robust optimization for drone delivery facility location and allocation problem

摘要

本文围绕无人机服务模式下的城市末端配送过程，通过两阶段数学规划模型，优化无人机配送设施选址及服务匹配策略。利用分布鲁棒优化方法刻画了需求不确定性对选址分配决策的影响，并将原始模型等价转化成混合整数二次锥规划问题，提出外逼近算法提升求解效率。本文依据上海市松江区的药品配送数据实例进行仿真验证，结果表明：本文提出的外逼近算法能将商业求解的运算时间缩短32.07%；与确定性模型对比，本文提出的分布鲁棒优化模型能够在需求具有较大波动时将系统总利润提高4.21倍。

Abstract

This paper focuses on the urban last-mile delivery process under the drone service model， optimizing drone delivery facility location and service allocation strategies through a two-stage mathematical programming model. The impact of demand uncertainty on location-allocation decisions is characterized using a distributionally robust optimization method. The original model is equivalently transformed into a mixed-integer second-order cone programming problem， and an outer approximation algorithm is proposed to improve solving efficiency. Simulation verification is conducted based on the pharmaceutical delivery data in Songjiang District， Shanghai. The results indicate that the proposed outer approximation algorithm can reduce the computational time of commercial solvers by 32.07%； compared with the deterministic model， the proposed distributionally robust optimization model can increase the system's total profit by 4.21 times when demand is highly volatile.

Graphical abstract

关键词

交通运输规划与管理 / 无人机配送 / 设施选址分配 / 需求不确定性 / 分布鲁棒优化 / 外逼近法 / 低空经济

Key words

transportation planning and management / drone delivery / facility location-allocation / demand uncertainty / distributionally robust optimization / outer approximation / low-altitude economy

0 引 言

当前，传统通用航空、无人驾驶航空等低空经济各业态的应用场景逐渐增多，低空经济的发展已成为必然趋势^［1］。由于无人机配送可以不受地理条件限制地将货物快速送达目的地^［2］，现已在全球范围内展开了商用实践。2019年，UPS将药品通过无人机直接送到急需的顾客手中或配送到往来不便的村庄^［3］；2022年8月，美团无人机配送已在深圳部分区域落地，可为近2万户居民服务。无人机配送设施选址是布局无人机配送网络的基础条件，如京东将在陕西省建设超过100个无人机机场^［4］。无人机配送服务设施主要涵盖货物存储、装载、调度管理、维护与保养等。无人机配送设施主要包括两层网络结构：设施层和需求层。设施层主要指无人机配送系统中固定的基础设施；需求层则涉及最终用户的需求管理。

既有文献对无人机配送设施选址及配送问题展开了大量研究。Golabi等^［5］提出了无人机救灾配送中心选址模型；Hong等^［6］研究了考虑飞行范围限制的无人机充电站选址问题；钱欣悦等^［7］研究了不同空域设置背景下末端配送中无人机起降点选址及分配问题；Chauhan等^［8］基于最大覆盖理论构建了无人机设施选址及分配模型。

目前已有大量研究考虑不确定性对选址决策的影响。李艳波等^［9］综合考虑路径上电动汽车流量和初始剩余电量的随机特性；赖志柱等^［10］研究了考虑物资需求不确定性的应急物流中心的选址及应急物资的调运问题。鲁棒优化将不确定性限制在一个可能的集合中，鲜斌等^［11］提出了一种新型非线性鲁棒控制策略；何勇等^［12］研究需求不确定情形下的载机平台调度问题；Zhu等^［13］提出面向人道主义救援过程的无人机选址-分配-指派鲁棒优化模型。然而鲁棒优化模型由于其过于保守导致日常运营成本过高。

分布鲁棒优化方法相比传统鲁棒优化可更好地利用实际运营中的数据描述随机参数分布的模糊集合。鉴于此，本文采用的基于矩信息的分布鲁棒优化模型，能够充分利用历史数据中的有效信息，在理论上保障了数据量与模糊集大小的统计关系^［14］，能够得到给定数据量前提下真实分布在模糊集中的概率，克服了传统鲁棒优化过分保守的缺陷^［15］。本文针对需求不确定场景，基于分布鲁棒优化方法建立无人机配送设施选址及分配模型；为提高求解效率，依据线性松弛模型的凸性，提出外逼近算法。

1 问题描述与建模

1.1　问题描述

本文基于城市无人机商业配送场景，研究无人机配送设施选址及分配优化问题。考虑到无人机的有效载荷和飞行里程双重限制，末端配送中，无人机一般执行直送服务。即无人机每次从配送设施出发，都只飞往一个需求点进行配送服务，之后直接返回配送设施。本文研究的配送服务设施是无人机枢纽站，即无人机的起降点。这些枢纽站是无人机配送网络的终端节点，负责执行具体的配送任务，包括货物的装载、起飞、配送以及返回^［16］。

研究问题考虑一个包含备选配送设施、需求点和无人机的系统，已知无人机的总量为定值，需求点的个数和位置确定，每个需求点的需求量都不超过无人机的最大载重。配送设施的位置在给定的备选配送设施集合中选取。为模拟实际情况，首先在当前的周期内，从备选配送设施集合中选择若干套配送设施，作为开放使用的无人机配送设施，即选址决策。其次，对每个选定的配送设施分配一定数量的无人机，即无人机分配决策。由于无人机配送需要规定航线，提前审批，需要划定配送设施的服务范围，每个需求点只能由一套指定的配送设施来提供服务，即划定服务范围决策。最后，在分配无人机和划定服务范围的条件下，对每架无人机具体服务哪些需求点进行调度，即无人机服务匹配决策。配送设施选址、无人机分配和划定服务范围决策如图1所示。

1.2　符号含义

模型涉及的符号定义如表1所示。加粗符号代表矩阵或者向量。

1.3　确定性模型建立

模型的优化目标是使配送活动的日均总利润最大化。配送设施运营成本设为c_{\mathrm{f}}，每个配送设施的日均运营成本为c_i，所有配送设施的运行成本表示为：

c_{\mathrm{f}}={\sum }_{i\in I}{c}_i{x}_i（1）

无人机使用成本为c_{\mathrm{d}}，每架无人机的日均成本c_{\mathrm{0}}是相同的，因此，无人机的使用成本可表示为：

c_{\mathrm{d}}=c_{\mathrm{0}}{\sum }_{i\in I}{\sum }_{k\in K}{z}_{ik}（2）

无人机配送收益为\vartheta，对单位重量的货物收费为\theta，且人机配送收益与无人机配送所满足的需求量成正比，因此，无人机的配送收益表示为：

\vartheta =\theta {\sum }_{i\in I}{\sum }_{j\in J}{\sum }_{k\in K}{w}_{ijk}{q}_j（3）

无人机配送的总利润为-c_{\mathrm{f}}-c_{\mathrm{d}}+\vartheta，需求覆盖率为\eta，且无人机至少为⌈\eta m⌉个需求点进行配送，因此，无人机配送设施选址和分配联合优化模型可表示为：

\begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}-{\sum }_{i\in I}{c}_i{x}_i-c_{\mathrm{0}}{\sum }_{i\in I}{\sum }_{k\in K}{z}_{ik}+\\ \theta {\sum }_{i\in I}{\sum }_{j\in J}{\sum }_{k\in K}{w}_{ijk}{q}_j\end{array}（4）

\mathrm{s}.\mathrm{t}.{\sum }_{i\in I}{x}_i\le P（5）

y_{ij}\le x_i,\forall i\in I,j\in J（6）

{\sum }_{i\in I}{y}_{ij}=\mathrm{1},\forall j\in J（7）

z_{ik}\le x_i,\forall i\in I,k\in K（8）

{\sum }_{i\in I}{z}_{ik}\le \mathrm{1},\forall k\in K（9）

w_{ijk}\le y_{ij},\forall i\in I,j\in J,k\in K（10）

w_{ijk}\le z_{ik},\forall i\in I,j\in J,k\in K（11）

{\sum }_{i\in I}{\sum }_{j\in J}{w}_{ijk}{l}_{ij}\le L,\forall k\in K（12）

{\sum }_{i\in I}{\sum }_{k\in K}{w}_{ijk}\le \mathrm{1},\forall j\in J（13）

{\sum }_{k\in K}{z}_{ik}\le R,\forall i\in I（14）

{\sum }_{i\in I}{\sum }_{j\in J}{\sum }_{k\in K}{w}_{ijk}\ge ⌈\eta m⌉（15）

x_i,y_{ij},z_{ik},w_{ijk}\in \{\mathrm{0,1}\},\forall i\in I,j\in J,k\in K（16）

其中，目标函数式（4）表示最大化日均总利润，约束（5）表示选定的配送设施不超过P个；约束（6）表示只有选定的配送设施才能服务需求点；约束（7）表示每个需求点必须指定给一个配送设施进行服务；约束（8）表示只有选定的配送设施才能分配无人机；约束（9）表示每架无人机只能分配给一个配送设施；约束（10）表示无人机只能服务配送设施覆盖范围内的需求点；约束（11）表示被分配给配送设施的无人机，才能为对应的需求点完成配送任务；约束（12）表示无人机的里程限制；约束（13）表示需求至多被配送一次；约束（14）表示每个选定的配送设施配备的无人机不能超过R架；约束（15）表示配送服务覆盖率不得低于\eta；约束（16）说明了决策变量x_i\mathrm{、}{y}_{ij}\mathrm{、}{z}_{ik}\mathrm{、}{w}_{ijk}均为0-1变量。

1.4　分布鲁棒优化模型

相关优化决策可以划分成两个阶段：首先确定无人机配送设施的选址，对选定的配送设施需要划定配送范围并配置无人机，即第一阶段决策变量x_i\mathrm{、}{y}_{ij}\mathrm{、}{z}_{ik}；然后第二阶段确定无人机的分配情况，即第二阶段决策变量w_{ijk}。

实际情况中，需求点的日均需求量\mathit{q}是一个具有不确定性的参数。在进行选址和分配决策前，无法知道\mathit{q}的真实取值，只能凭已有的统计量信息来进行建模和决策，在原模型中，需求量\mathit{q}仅出现在目标函数中，即需求量\mathit{q}只会直接影响到无人机配送的收益。研究假设不论实际需求如何波动，都不会超过无人机的有效载荷，因此，不需要再额外添加约束条件。不妨将收益函数简单记作e(\mathit{q})，为了使配送系统的总利润最大化，在做出调度决策时必须使配送收益尽可能大，即寻找\mathrm{m}\mathrm{a}\mathrm{x}e(\mathit{q})。令\tilde{\mathit{q}}表示需求对应的随机变量。当需求为随机变量且采用分布鲁棒优化方法处理不确定参数时，随机变量的分布函数P不确定且属于模糊集合\mathrm{𝒫}，为使模型具有较好的鲁棒性，即在任何可能发生的情况下都能保证可行，需要在\mathrm{𝒫}集合所有的分布函数中，找到使e(\tilde{\mathit{q}})的期望值最小（目标收益最差）的分布，并使在取最差分布的条件下，仍能保证配送系统的总利润最大化。因此，原问题的目标函数可以改为：

\mathrm{m}\mathrm{a}\mathrm{x}-c_{\mathrm{f}}-c_{\mathrm{d}}+\underset{P\in \mathrm{𝒫}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathbb{E}}_P[e(\mathit{q})]（17）

由原始的目标函数可知，配送的收益e(\mathit{q})是关于\mathit{q}的线性函数，因此，{\mathbb{E}}_P[e(\tilde{\mathit{q}})]只取决于\tilde{\mathit{q}}的期望值{\mathbb{E}}_P(\tilde{\mathit{q}})。\tilde{\mathit{q}}为一随机向量，在本文中，假设需求量\tilde{\mathit{q}}的分布函数为P，该函数的具体形式未知。但可以建立一个关于分布函数P的不确定集\mathrm{𝒫}：

\mathrm{𝒫}=\left\{P:{\left({\mathbb{E}}_P(\tilde{\mathit{q}})-\mathit{\mu }\right)}^{\mathrm{T}}{\mathit{\Sigma }}^{-\mathrm{1}}\left({\mathbb{E}}_P(\tilde{\mathit{q}})-\mathit{\mu }\right)\le {ϵ}^{\mathrm{2}}\right\}（18）

集合（18）表示，随机变量的一阶矩{\mathbb{E}}_P(\tilde{\mathit{q}})在一个由历史数据观测得到的均值\mathit{\mu }、协方差矩阵\Sigma \succcurlyeq \mathrm{0}和半径ϵ定义的椭球集内。所有满足这一条件的分布函数P共同构成了集合^［15］。

由原模型可知，e(\tilde{\mathit{q}})是一个关于第一阶段决策变量\mathit{x}\mathbf{、}\mathit{y}\mathbf{、}\mathit{z}的函数，因此，将其改为g(\mathit{x}\mathbf{、}\mathit{y}\mathbf{、}\mathit{z}\mathbf{、}\tilde{\mathit{q}})。根据上述分析，可以建立商业配送情景下考虑需求不确定性的无人机配送设施选址和分配两阶段分布式鲁棒模型：

第一阶段表示为：

\begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}-{\sum }_{i\in I}{c}_i{x}_i-c_{\mathrm{0}}{\sum }_{i\in I}{\sum }_{k\in K}{z}_{ik}+\\ \theta \underset{P\in \mathrm{𝒫}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathbb{E}}_p[g(\mathit{x},\mathit{y},\mathit{z},\tilde{\mathit{q}})]\end{array}（19）

\mathrm{s}.\mathrm{t}.x_i,y_{ij},z_{ik}\in \{\mathrm{0,1}\},\forall i\in I,j\in J,k\in K

(\mathrm{5})~(\mathrm{9})\mathrm{、}(\mathrm{14})（20）

式中：目标函数式（19）表示即使在需求的分布为最差情况时，仍然使总利润最大化；约束（5）~（9）、（14）的形式及含义同上文所述；约束（20）表示决策变量x_i\mathrm{、}{y}_{ij}\mathrm{、}{z}_{ik}均为0-1变量。

第二阶段表示为： g(\mathit{x}\mathbf{、}\mathit{y}\mathbf{、}\mathit{z}\mathbf{、}\tilde{\mathit{q}})=\mathrm{m}\mathrm{a}\mathrm{x}{\sum }_{i\in I}{\sum }_{j\in J}{\sum }_{k\in K}{w}_{ijk}{q}_j（21）

\mathrm{s}.\mathrm{t}.w_{ijk}\in \{\mathrm{0,1}\},\forall i\in I,j\in J,k\in K(\mathrm{14})~(\mathrm{16})\mathrm{、}(\mathrm{19})（22）

其中，目标函数式（21）表示最大化覆盖需求；约束（14）~（16）（19）的形式及含义同上文所述；约束（22）表示决策变量w_{ijk}为0-1变量。

下面，使用一种分布鲁棒等价转化方法^［14］来对上述两阶段鲁棒模型进行处理，将该模型转化为一个可解的二阶锥规划（Second order cone program， SOCP）模型。

首先将原问题的目标函数进行转化，由式（19）（21）变为：

\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{\sum }_{i\in I}{c}_i{x}_i+c_{\mathrm{0}}{\sum }_{i\in I}{\sum }_{k\in K}{z}_{ik}+\\ \theta \underset{p\in \mathrm{𝒫}}{\mathrm{m}\mathrm{a}\mathrm{x}}{\mathbb{E}}_p[g(\mathit{x},\mathit{y},\mathit{z},\tilde{\mathit{q}})]\end{array}

g(\mathit{x},\mathit{y},\mathit{z},\tilde{\mathit{q}})=\mathrm{m}\mathrm{i}\mathrm{n}-{\sum }_{i\in I}{\sum }_{j\in J}{\sum }_{k\in K}{w}_{ijk}{q}_j（23）

由于原第二阶段的决策变量\mathit{w}与\mathit{q}维数不同，为方便推导，引入一个新的决策变量\mathit{u}，定义：

u_j=\sum _{i\in I}\sum _{k\in K}{w}_{ijk},j\in J（24）

根据约束（24）可知，u_j也为一个0-1变量，表示需求点j是否被无人机执行派送服务，若是则取1，否则取0。以\mathit{u}代换原来目标函数中的\mathit{w}，则式（23）又可以变为：

\begin{array}{c}g\left(\mathit{x},\mathit{y},\mathit{z},\tilde{\mathit{q}}\right)=\mathrm{m}\mathrm{i}\mathrm{n}-\sum _{j\in J}{u}_j{q}_j\end{array}（25）

通过引入一个辅助决策变量a，该两阶段鲁棒优化问题可以被转化为如下的形式：\underset{\mathit{x},\mathit{y},\mathit{z},\mathit{w},\mathit{u},a}{\mathrm{m}\mathrm{i}\mathrm{n}}\sum _{i\in I}{c}_i{x}_i+c_{\mathrm{0}}\sum _{i\in I}\sum _{k\in K}{z}_{ik}-\theta \left({\mathit{\mu }}^{\mathrm{T}}\mathit{u}-ϵa\right)

\mathrm{s}.\mathrm{t}.{‖\mathit{\Sigma }}^{\frac{\mathrm{1}}{\mathrm{2}}}{\mathit{u}‖}_{\mathrm{2}}\le a（26）

u_j\in \left\{\mathrm{0,1}\right\},\forall j\in J

\begin{array}{c}a\ge \mathrm{0}\end{array}

\left(\mathrm{5}\right)~\left(\mathrm{16}\right)\mathrm{、}\left(\mathrm{24}\right)

由于最初对模型的目标函数做了取负的变换，从最大化变为最小化，因此，最终需要对得到的最优值取负，得到配送利润的最大值。

2 求解方法

针对所建立的混合整数二次锥规划（Mixed-integer second order cone programming， MISOCP）问题，本文采用外逼近算法（Outer approximation method， OAM）进行优化。OAM是一种精确求解混合整数非线性规划问题的全局优化方法，它将原始问题分解成较简单的混合整数线性规划（Mixed integer linear programming，MILP）主问题，和一个非线性规划（Nonlinear program，NLP）子问题，通过不断交替迭代求解可以得到一系列原问题的上界和下界^［15］。在原始问题满足一定条件时，OAM能够通过有限次迭代获得原始问题的最优解，相较于演化算法在算法精确度上具有显著优势。

首先，在命题中证明该MISOCP问题的线性松弛问题是凸的。

命题1 定义函数\delta (\mathit{u},a)=\sqrt[]{{\mathit{u}}^{\mathrm{T}}\mathit{\Sigma }\mathit{u}}-a，则\delta (\mathit{u},a)的线性松弛项为凸函数。

证明：

根据平方根法将随机变量的协方差矩阵\mathit{\Sigma }表示为\mathit{\Sigma }=\mathit{L}{\mathit{L}}^{\mathrm{T}}，\sqrt[]{{\mathit{u}}^{\mathrm{T}}\mathit{\Sigma }\mathit{u}}=\sqrt[]{{\mathit{u}}^{\mathrm{T}}\mathit{L}{\mathit{L}}^{\mathrm{T}}\mathit{u}}，令{\mathit{A}}^{\mathrm{T}}={\mathit{u}}^{\mathrm{T}}\mathit{L}，{\mathit{A}}^{\mathrm{T}\mathrm{*}}=\left({\mathit{A}}_{\mathrm{1}},{\mathit{A}}_{\mathrm{2}},\cdots ,{\mathit{A}}_I\right)可得\sqrt[]{{\mathit{u}}^{\mathrm{T}}\mathit{\Sigma }\mathit{u}}={\left|\left|{\mathit{A}}^{\mathrm{T}\mathrm{*}}\right|\right|}_{\mathrm{2}}，对于任意向量{\mathit{A}}_{\mathrm{1}}^{\mathrm{T}\mathrm{*}}和{\mathit{A}}_{\mathrm{2}}^{\mathrm{T}\mathrm{*}}而言，三角不等式成立。\begin{array}{c}\left|\left|\lambda {\mathit{A}}_{\mathrm{1}}^{\mathrm{T}\mathrm{*}}+\left(\mathrm{1}-\lambda \right){\mathit{A}}_{\mathrm{2}}^{\mathrm{T}\mathrm{*}}\right|\right|\le \left|\left|\lambda {\mathit{A}}_{\mathrm{1}}^{\mathrm{T}\mathrm{*}}\right|\right|+\left|\left|\left(\mathrm{1}-\lambda \right){\mathit{A}}_{\mathrm{2}}^{\mathrm{T}\mathrm{*}}\right|\right|\end{array}，即证明函数\delta (\mathit{u},a)为凸函数，命题成立。

2.1　初始解构建

由于模型的任何可行解均可看作OAM的初始解，给定一个足够保守的可行解，确保满足所有约束条件，作为原问题的初始解：

\left\{\begin{array}{l}{x}_i^{\mathrm{0}}=\left\{\begin{array}{l}\mathrm{1},i=\mathrm{1}\\ \mathrm{0},\mathrm{其他}\end{array}\right.\\ y_{ij}^{\mathrm{0}}=\left\{\begin{array}{l}\mathrm{1},i=\mathrm{1},j=\mathrm{1,2},\cdots ,J\\ \mathrm{0},\mathrm{其他}\end{array}\right.\\ z_{ik}^{\mathrm{0}}=\left\{\begin{array}{l}\mathrm{1},i=\mathrm{1},k=\mathrm{1,2},\cdots ,K\\ \mathrm{0},\mathrm{其他}\end{array}\right.\\ w_{ijk}^{\mathrm{0}}=\left\{\begin{array}{l}\mathrm{1},i=\mathrm{1},j=\mathrm{1,2},\cdots ,J,k=\mathrm{1,2},\cdots ,K\\ \mathrm{0},\mathrm{其他}\end{array}\right.\end{array}\right.（27）

2.2　NLP子问题

当固定模型中的整数为变量时，其退化为一个经典非线性规划问题，仅与连续变量有关。令{\tilde{\mathit{u}}}^h为在第h次迭代时的最优解。在转化后的模型中，非线性项可以由整数变量的值为：{\tilde{a}}^h=\sqrt[]{{\left({\tilde{\mathit{u}}}^h\right)}^{\mathrm{T}}\mathit{\Sigma }{\tilde{\mathit{u}}}^h}。

因而，本文提出的OAM不用计算较难求解的非线性子问题，可以直接代入取值。代入{\tilde{\mathit{u}}}^h之后，因为非线性项总是满足的，所以当前返回值为可行解，可以得到原问题的一个上界。

2.3　MILP主问题

模型的主问题通过不断加入外逼近切更新而得到原问题的下界，外逼近切通过MILP主问题的最优解{\tilde{\mathit{u}}}^h\mathrm{、}{\tilde{a}}^h得到，命题2给出了非线性约束（26）的线性有效不等式。

命题2 约束（26）在第h次迭代中的外逼近切为：

\begin{array}{c}{\left({\tilde{\mathit{u}}}^h\right)}^{\mathrm{T}}\mathit{\Sigma }{\mathit{u}}^h-{\tilde{a}}^{h}a\le \mathrm{0}\end{array}（28）

证明：由命题1有\delta (\mathit{u},a)为凸函数，对其做泰勒展开可得，\delta ({\tilde{\mathit{u}}}^h,{\tilde{a}}^h)+\nabla \delta {\left({\tilde{\mathit{u}}}^h,{\tilde{a}}^h\right)}^{\mathrm{T}}\cdot \left[\begin{array}{c}\mathit{u}-{\tilde{\mathit{u}}}^h\\ a-{\tilde{a}}^h\end{array}\right]\le \delta (\mathit{u},a)\le \mathrm{0}。

其中，\nabla \delta {\left({\tilde{\mathit{u}}}^h,{\tilde{a}}^h\right)}^{\mathrm{T}}=\left[\frac{{\left({\tilde{\mathit{u}}}^h\right)}^{\mathrm{T}}\mathit{\Sigma }}{\sqrt[]{{\left({\tilde{\mathit{u}}}^h\right)}^{\mathrm{T}}\mathit{\Sigma }{\tilde{\mathit{u}}}^h}}-\mathrm{1}\right]，又因为{\tilde{a}}^h=\sqrt[]{{\left({\tilde{\mathit{u}}}^h\right)}^{\mathrm{T}}\mathit{\Sigma }{\tilde{\mathit{u}}}^h}，上式可简化为{\left({\tilde{\mathit{u}}}^{\mathrm{h}}\right)}^{\mathrm{T}}\mathit{\Sigma }{\mathit{u}}^h-{\tilde{a}}^{h}a\le \mathrm{0}，与式（28）相同，证毕。

引入辅助变量\Omega，将目标函数放到约束条件中，综上所述，\mathrm{M}\mathrm{I}\mathrm{L}\mathrm{P}汇总为：

\begin{array}{c}\mathrm{m}\mathrm{i}\mathrm{n}\end{array}\Omega

\mathrm{s}.\mathrm{t}.\sum _{i\in I}{c}_i{x}_i+c_{\mathrm{0}}\sum _{i\in I}\sum _{k\in K}{z}_{ik}-\theta \left({\mathit{\mu }}^{\mathrm{T}}{\tilde{\mathit{u}}}^h-\mathrm{ϵ}a\right)\le \Omega

\begin{array}{c}{\left({\tilde{\mathit{u}}}^h\right)}^{\mathrm{T}}\mathit{\Sigma }{\tilde{\mathit{u}}}^h-{\tilde{a}}^{h}a\le \mathrm{0}\end{array}

\begin{array}{c}\Omega \in R,\mathit{\mu }\in {\mathbb{R}}^n,a\ge \mathrm{0}\end{array}

\left(\mathrm{5}\right)~\left(\mathrm{16}\right)\mathrm{、}\left(\mathrm{24}\right)（29）

2.4　算法主要流程

基于以上对主问题和子问题的介绍，本文模型的外逼近法具体算法步骤如下：

步骤1 给定初始解为x_i^{\mathrm{0}}\mathrm{、}{y}_{ij}^{\mathrm{0}}\mathrm{、}{z}_{ik}^{\mathrm{0}}\mathrm{、}{w}_{ijk}^{\mathrm{0}}，选取\mathrm{U}{\mathrm{B}}^{\mathrm{0}}=\mathrm{\infty }，\mathrm{L}{\mathrm{B}}^{\mathrm{0}}-\mathrm{\infty }，h=\mathrm{0}，c={\mathrm{10}}^{-\mathrm{6}}。

步骤2 若\mathrm{U}{\mathrm{B}}^{\mathrm{0}}-\mathrm{L}{\mathrm{B}}^{\mathrm{0}}\le \mathrm{c}，则停止计算，并返回当前解。

步骤3 依{\tilde{x}}_i^h\mathrm{、}{\tilde{y}}_{ij}^h\mathrm{、}{\tilde{z}}_{ik}^h\mathrm{、}{\tilde{w}}_{ijk}^h计算{\tilde{\mathit{u}}}^h、{\tilde{a}}^h，其中{\tilde{a}}^h=\sqrt[]{{\left({\tilde{\mathit{u}}}^{\mathrm{h}}\right)}^{\mathrm{T}}\mathit{\Sigma }{\tilde{\mathit{u}}}^{\mathrm{h}}}，以此更新下界。

步骤4 更新上界\mathrm{U}{\mathrm{B}}^{\mathrm{r}}=\sum _{i\in I}{c}_i{x}_i+c_{\mathrm{0}}\sum _{i\in I}\sum _{k\in K}{z}_{ik}-\theta \left({\mathit{\mu }}^{\mathrm{T}}{\tilde{\mathit{u}}}^h-ϵ{\tilde{a}}^h\right)，将OA切（22）并入MILP并求解。

步骤5 更新整数变量{\tilde{x}}_i^h\mathrm{、}{\tilde{y}}_{ij}^h\mathrm{、}{\tilde{z}}_{ik}^h\mathrm{、}{\tilde{w}}_{ijk}^h，h=h+\mathrm{1}，返回步骤2。

3 算例分析

本文以上海市松江区的实际数据集为例，研究无人机在药品配送过程中的相关优化决策。在案例中，以区域内药店为无人机配送设施备选点，以各乡镇中心点位置为需求点。模型中其他主要参数包括无人机的有效飞行里程、最大载重、配送设施运作成本、无人机使用成本以及单位配送收益。本配送场景所使用的无人机拟定为迅蚁多旋翼无人机RA3，其最大有效载荷为6 kg，货箱容积为25 L，最大航程为25 km。每架无人机每天的使用成本为无人机充电成本和电池损耗成本之和，约为11元。根据无人机枢纽站售价和使用年限折算，每台每天的成本约为27元。对于药店放置枢纽站和无人机所需要的相关改造工作，可能产生额外的一部分成本，这部分成本会因各个设施的具体条件而不同，不妨假设一次性改造成本为10 000~30 000元不等，以使用期限5年进行折算，每天的成本为6~16元。人力成本为每天100元，设施的每日运营成本为133~143元。迅蚁每架无人机的最大载荷为6 kg，每千克产品的配送价格为5~6元，取中值计算即配送价格为5.5元/kg。根据松江区人口及乡镇数量设定需求点。根据配送设施的间距及无人机配送里程限制筛选出10家门店组成备选配送设施，如图2所示。

所有计算皆在Lenovo XiaoXin Air 15IKBR电脑上进行，CPU型号为Intel（R） Core（TM） i5-8250U CPU @ 1.60 GHz，内存大小为8 GB，内存频率为2 400 MHz。所有代码均由Python编程实现。

3.1　算法性能分析

在上述数据处理前提下，进行算法性能分析。

首先，根据历史观测数据的样本均值和样本方差，分别记作\widehat{\mathit{\mu }}和{\mathit{\Sigma }}_{\mathrm{0}}，ϵ=\mathrm{0.01}，以计算估计随机变量均值{\mathbb{E}}_P[e(\tilde{\mathit{q}})]的分布集合。为了较为直观地进行算法性能分析，分别用求解器、考虑强制约束的遗传算法（规模=100，迭代次数上限为500次，交叉率为0.8，变异率为0.01）和OAM对不同算例进行多次计算，求解时间取平均值，在多个算例下的求解结果汇总于表2，

可知，大部分算例下，OAM相比求解器有更快的求解速度，尤其是在当求解器的求解时间超过500 s后，对应算例下的OAM求解速度平均提高32.07%；遗传算法与OAM相比求解速度更快，但是求解结果相比OAM求得的最优解差距过大，仅能在有限迭代步骤中得到可行解，无法保障解的质量。

3.2　模型效果对比

初始模型的决策变量分为选址决策（x_i，z_{ik}）和分配决策（y_j,w_{ijk}）。本文对确定性模型和分布鲁棒模型进行对比，以比较二者在不同实证数据集下的表现。

首先，利用蒙特卡洛仿真方法在随机参数不确定集合中生成S（=100）个可能的取值（M^s），其中s=\mathrm{1},\cdots ,S。M^s中的每一项{\widehat{M}}_{ij}服从某正态分布，即M_{ij}~N\left({\mu }_{ij},{\sigma }_{ij}^{\mathrm{2}}\right)，其中{\mu }_{ij}为配送设施到需求点之间的货物需求量历史数据方差。各算例在本节分析中表现基本一致，本节\left(\widehat{m},\widehat{n},\widehat{r}\right)采用(\mathrm{20,20,10})举例说明，各项分析均基于多次实验平均取值数据。

表3展示了确定性模型和分布鲁棒优化模型的配送设施选址位置、配置的无人机数量以及单个无人机服务的需求点范围，为了更为直观地表示两类优化模型的服务匹配关系，绘制图3及图4，展示松江区范围内的无人机设施选址、配置及分配方案。

模型运行结果对比见表4，分布鲁棒模型的选址成本等于确定性成本，总体来看分布鲁棒模型的配送阶段利润得到提高，总利润提高至确定性模型的4.21倍。将随机生成的需求量数据代入模型得到的第一阶段决策中，鲁棒模型得到的配送利润总是优于确定性模型。这是因为鲁棒模型能使需求分布较差时仍然保持总体利润的最大化，而确定性模型不能预见需求量的波动，得到的最优决策容易受不确定性的干扰。由实验结果可知：考虑需求不确定性的鲁棒模型具有较好的效果，模型的鲁棒性得到了验证。

3.3　灵敏度分析

本节主要对无人机安全系数、需求的变异系数进行灵敏度分析。在算例\left(\widehat{m},\widehat{n},\widehat{r}\right)=(\mathrm{10,33,12})和\left(\widehat{m},\widehat{n},\widehat{r}\right)=(\mathrm{40,20,5})下进行灵敏度分析，分析时控制其他参数一致。

3.3.1　对安全系数\mathrm{\lambda }的灵敏度分析

各算例在本节分析中表现基本一致，令\left(\widehat{m},\widehat{n},\widehat{r}\right)=(\mathrm{10,33,12})。无人机在飞行过程中，会受到天气状况、起降行为、实际载重以及电池本身特性等因素的影响，其有效飞行里程往往不能达到额定的最大值。为保证无人机在飞行过程中的安全性，在实际情况中，往往不以无人机最大飞行里程或最长飞行时间来进行计算，而是要考虑无人机电池的安全系数，即无人机使用的电量不能超过其电池容量的一定比例。若仍认为电量消耗和飞行距离成正比，可以将这一系数加在本模型中的里程约束上，即每架无人机总飞行里程不超过最大里程乘以安全系数。

如表5所示，随着安全系数\mathrm{\lambda }的提高，无人机有效最大飞行里程不断缩短。在使用同样数量无人机的情况下，安全系数越高，则无人机服务半径越小，其需求覆盖率也越低，配送总利润也随之减少。

3.3.2　对变异系数的灵敏度分析

各算例在本节分析中表现基本一致，令(\widehat{m},\widehat{n},\widehat{r})=(\mathrm{40,20,5})。假设各需求点之间需求量互不相关，则矩阵\mathit{\Sigma }中的对角线元素为各点需求量期望的方差，其余元素为0。对于需求量期望的波动，可以用需求的变异系数来计算，其含义是标准差与均值之比。图5展示了总利润值随鲁棒椭球集半径和需求变异系数CV变化的曲线。

由图5可知，随着椭球集半径参数\mathrm{ϵ}和需求变异系数的增大，需求量期望的变动范围随之变大，需求量期望的最差分布情况逐步变差，总利润值的期望不断下降。当椭球集半径为0时，需求的波动为0，分布鲁棒模型退化为确定性模型，因此，3条曲线交汇于点（0，107）；当椭球集半径从0提高至1.5时，在\mathrm{C}\mathrm{V}=\mathrm{0.1}和0.3两种情况下的总利润分别降低了11.2%和34.9%，表明鲁棒程度增加会使求解结果更为保守。尽管鲁棒模型和确定性模型得到的决策可能是相同的，但鲁棒模型总能给出该种决策在最差情况下的总利润期望值，使得系统即使面临较大波动时的效果依然良好。

4 结 论

（1）外逼近算法在求解本文模型的大规模算例时，相较于商业求解器Gurobi，其求解速度可提高32.07%。

（2）当面临需求波动时，分布鲁棒优化模型的系统总成本与确定性模型持平，分配调度阶段利润更高，稳定性提升，总利润提高至确定性模型的4.21倍。

（3）通过灵敏度分析，当分布鲁棒集半径ϵ=\mathrm{0.3}时，若需求变异系数由0增加到1.5，将使得总利润降低34.9%，需求不确定性会显著降低系统总利润。

参考文献

[1] 张世昌, 许芳菲. 低空经济生成与发展的伦理考量——基于马克思空间思想的视角[J]. 北京航空航天大学学报:社会科学版, 2024, 37(6): 1-9.

[2] Zhang Shi-chang, Xu Fang-fei. Ethical considerations on the generation and development of low-altitude economy: a perspective based on marx's spatial[J]. Journal of Beihang University (Social Sciences Edition), 2024, 37(6): 1-9.

[3] 何颖, 樊俊松, 王巍, 等. 无人机空地安全通信与航迹规划的多目标联合优化方法[J]. 吉林大学学报:工学版, 2023, 53(3): 913-922.

[4] He Ying, Fan Jun-song, Wang Wei, et al. Joint optimization of secure communication and trajectory planning in unmanned aerial vehicle air⁃to⁃ground[J]. Journal of Jilin University (Engineering and Technology Edition), 2023, 53(3): 913-922.

[5] 陈明舟. 无人机物流配送服务动向[J]. 物流工程与管理, 2020, 42(10): 90-92.

[6] Chen Ming-zhou. Trends in logistics delivery services using UAV[J]. Logistics Engineering and Management, 2020, 42(10): 90-92.

[7] 陈刚, 付江月. 军民融合背景下无人机配送中心选址问题研究[J]. 计算机工程与应用, 2019, 55(8): 226-231.

[8] Chen Gang, Fu Jiang-yue. Research on the location problem of drone delivery centers under the background of military-civilian integration[J]. Computer Engineering and Applications, 2019, 55(8): 226-231.

[9] Golabi M, Shavarani S M, Izbirak G. An edge-based stochastic facility location problem in UAV supported humanitarian relief logistics: a case study of Tehran earthquake[J]. Natural Hazards, 2017, 87(3):1545-1565.

[10] Hong I, Kuby M, Murray A T. A range-restricted recharging station coverage model for drone delivery service planning[J]. Transportation Research Part C: Emerging Technologies, 2018, 90:198-212.

[11] 钱欣悦, 张洪海, 张芳,等. 末端配送物流无人机起降点选址分配问题研究[J].武汉理工大学学报:交通科学与工程版, 2021, 45(4): 682-687.

[12] Qian Xin-yue, Zhang Hong-hai, Zhang Fang, et al. Research on location allocation of UAV landing points for terminal distribution logistics[J]. Journal of Wuhan University of Technology (Transportation Science & Engineering), 2021, 45(4): 682-687.

[13] Chauhan D, Unnikrishnan A, Figliozzi M. Maximum coverage capacitated facility location problem with range constrained drones[J]. Transportation Research Part C: Emerging Technologies, 2019, 99:1-18.

[14] 李艳波, 柳柏松, 姚博彬, 等. 考虑路网随机特性的高速公路换电站选址[J]. 吉林大学学报:工学版, 2023, 53(5): 1364-1371.

[15] Li Yan-bo, Liu Bai-song, Yao Bo-bin, et al. Location of electrical changing station of expressway considering stochastic characteristics of road network[J]. Journal of Jilin University (Engineering and Technology Edition), 2023, 53(5): 1364-1371.

[16] 赖志柱, 王铮, 戈冬梅, 等. 多目标应急物流中心选址的鲁棒优化模型[J]. 运筹与管理, 2020, 29(5): 74-83.

[17] Lai Zhi-zhu, Wang Zheng, Ge Dong-mei, et al. A multi-objective robust optimization model for emergency logistics center location[J]. Operations Research and Management Science, 2020, 29(5): 74-83.

[18] 鲜斌, 王光怡, 蔡佳明. 多无人机吊挂负载运输系统的非线性鲁棒控制设计[J]. 吉林大学学报: 工学版, 2024, 54(6): 1788-1795.

[19] Xian Bin, Wang Guang-yi, Cai Jia-ming. Nonlinear robust control design for multi unmanned aerial vehicles suspended payload transportation system[J]. Journal of Jilin University (Engineering and Technology Edition), 2024, 54(6): 1788-1795.

[20] 何勇, 张成义, 李姗姗. 基于两阶段鲁棒优化的无人机载机平台调度问题[J].系统工程学报, 2020, 35(6): 838-848.

[21] He Yong, Zhang Cheng-yi, Li Shan-shan. Unmanned aerial vehicle carriers scheduling problem based on two-stage robust optimization[J]. Journal of Systems Engineering, 2020, 35(6): 838-848.

[22] Zhu T K, Boyles S D, Unnikrisshnan A. Two-stage robust facility location problem with drones[J]. Transportation Research Part C: Emerging Technologies, 2022, 137: 103563.

[23] Zhang Y L, Shen Z J M, Song S J. Distributionally robust optimization of two-Stage lot-Sizing problems[J]. Production and Operations Management, 2016, 25(12): 2116-2131.

[24] Delage E, Ye Y Y. Distributionally robust optimization under moment uncertainty with application to data-driven problems[J]. Operations Research, 2010, 58 (3):595-612.

[25] 郑雨欣. 城市物流无人机运行中若干关键问题研究[D]. 天津: 中国民航大学空中交通管理学院, 2022.

[26] Zheng Yu-xin. Research on several key issues in urban logistics drone operations[D]. Tianjin: College of Air Traffic Management, Civil Aviation University of China, 2022.

基金资助

中央高校基本科研业务费专项项目(2023JBMC006)

国家自然科学基金项目(72471025)

国家自然科学基金项目(62433005)

国家自然科学基金项目(72101021)