基于概率围栏的无人机轨迹优化

无人机轨迹规划中基于概率地理围栏的迭代机会约束优化

摘要

机会约束优化因其能够刻画实际应用中的随机性，为解决具有不确定性的控制与规划问题提供了一个颇具前景的框架。本文基于机会约束优化方法，研究了带有概率地理围栏的无人机轨迹规划问题。在所考虑的问题中，模型的随机性（如地理围栏边界的不确定性）被纳入到问题建模中。通过采用一种新颖的基于采样的求解方法求解所建立的机会约束优化问题，获得最优无人机轨迹，同时将与地理围栏发生碰撞的概率限制在预设阈值内。此外，为了获得完全无碰撞轨迹（即不仅在离散时间步上，而且在整个时间范围内均避免碰撞），我们采用了一种迭代方案的思想：反复求解机会约束优化问题，直至在时间范围内的任意时刻均避免与概率地理围栏发生碰撞。最后，我们通过数值仿真验证了所提出方法的有效性。

索引词 —无人机，轨迹规划，电子围栏，机会约束优化。

一、引言

RECENTLY ，轨迹规划对于自动驾驶车辆[1]–[4],尤其是无人机（UAVs）[5]–[7]的应用而言，始终是一个根本性的关键问题。无人机轨迹规划通常指在特定空间内生成一条几何路径，同时确保所生成路径满足无人机运行的若干条件，例如安全（障碍物/碰撞规避）、最短时间或最低燃料消耗等。在这些基本条件中，安全性的保障尤为重要，因为无法保证安全可能导致严重后果。同时，地理空间限制（也称为geo-fence）有助于建立无人机不可穿越的受限区域的虚拟边界。电子围栏技术无疑有助于加强无人机操作的安全性。

事实上，过去几十年中，带地理围栏约束的无人机轨迹规划问题已受到越来越多的关注。为解决此类轨迹规划问题，文献中提出了大量方法，包括混合整数线性规划（MILP）方法[8],模型预测控制（MPC）[9],势场法（PFM）[10],等。在这些众多方法中，主要关注的典型挑战主要包括以下三个方面。

1) 非凸性 ：由于所考虑的规划空间中存在受限区域，无人机的安全区域本质上是非凸的，由此产生的优化模型也是非凸的。

2) 不确定性 ：考虑到无人机轨迹规划问题中存在多种随机因素，如运动扰动和环境不确定性。为了获得更可靠且鲁棒的无人机轨迹，必须将这些问题纳入考虑范围。

3) 离散化问题 ：当采用离散时间方法求解轨迹规划问题时，障碍物/碰撞规避通常仅在离散时间步上得到保证。因此，一个合理的问题是如何确保无人机在这些离散时间步之间的安全，并获得所谓的完全无碰撞轨迹，即在时间范围内的任意时刻均避开障碍物。

为了处理优化问题中的非凸性，一种常见的方法是引入一组二元变量，并利用大‐M约束[11]–[13]技术将非凸可行域重新表述为一组混合整数线性约束。尽管由于存在二元变量，所得的优化模型仍然是非凸的，但此类非凸性可通过一些现成的商业求解器（如CPLEX[14]和GUROBI[15]）轻松处理。

关于不确定性，大多数现有研究集中在无人机动力学建模中的随机性[16]–[21]。通过假设系统是线性的且噪声为白高斯噪声，无人机状态的分布可以在任意时刻表示为高斯统计量的传播。随后，线性机会约束被重新表述为关于无人机随机状态的均值和协方差的完全确定性形式[20],[21]。因此，可以应用标准优化方法来求解该问题。然而，由于这种重构会在确定性约束中引入高度非线性，因此计算复杂度仍是此类方法需要关注的问题。为进一步降低复杂度，[22]特别提出将球形碰撞区域扩大为半空间；但这种方法引入的保守性是其潜在缺点。此外，在少数研究环境不确定性的文献[23]–[25],中，研究了基于机会约束的快速探索随机树方法来解决轨迹规划问题。特别是，文献[23]将定位误差与运动扰动结合起来，并假设二者均为高斯分布。本文特别关注环境不确定性，即电子围栏边界的不确定性，但考虑其服从一般的概率分布，而非限制为白高斯分布。

为了获得完全无碰撞轨迹，文献中已采用多种方法来解决离散化问题。首先，应注意该问题仅在采用离散时间方法时出现，一种自然的思路是通过增加离散时间步数来实现更精细的离散化。然而，并没有理论保证需要多少时间步长才能确保无碰撞轨迹。此外，考虑到优化问题的规模随时间步长数量的增加而迅速增长，因此计算效率在此情况下也是一个关键问题。其次，另一种解决方案是扩大受限区域[8],[26],的实际尺寸，以确保无人机在这些离散时间步之间的安全性。类似地，仍然缺乏理论依据来确定受限区域应扩大多少；甚至更糟的是，如果引入过多保守性，优化模型可能变得不可行。第三，如近期研究[20],[27],所建议，当仅将多边形视为受限区域时，可通过确保同侧条件来生成无碰撞轨迹。然而，如[20]本身所述，在某些基本情况下，所得解的最优性可能会受到影响。

本文中，我们研究了考虑上述所有挑战的无人机轨迹规划问题。具体而言，我们的目标是在如下条件下生成最优无人机轨迹：假设无人机需要规划一条从初始位置到目的地的轨迹，同时该轨迹在控制能耗最小的意义上是最优的。假设环境中存在电子围栏规划空间。为了保证安全，在规划轨迹时应确保无人机能够避免与任何电子围栏发生碰撞。此外，考虑到电子围栏边界存在不确定性，碰撞规避应在概率方式下予以考虑。

值得指出的是，所考虑的无人机轨迹规划问题涉及上述三个挑战，本质上难以求解。针对这些困难，我们的思路是：i）引入一组二元决策变量以处理可行集的非凸性；ii）使用机会约束来刻画电子围栏的不确定性；iii）通过迭代求解机会约束优化，直到获得完全无碰撞轨迹。更具体地，我们将电子围栏的不确定性建模为边界距离的随机性，并确保避障概率必须高于预设阈值。对于非凸可行集，我们首先将电子围栏表示为一个由一组线性不等式定义的通用凸多边形，然后采用大‐M约束的思想，强制至少一个不等式被违反。这将得到一个带机会约束的混合整数二次规划（CC‐MIQP）问题。为求解该CC‐MIQP问题，我们开发了一种新的基于采样的求解方法。与现有的基于采样或基于场景的方法不同，我们的方法利用了所考虑的机会约束问题的特定结构，并建立在一种新颖的样本选择过程之上。通过设计适当的选择准则，仅选取部分采样场景，我们保证：i）通过求解近似确定性问题得到的解必须满足原始机会约束优化的可行性；ii）随着样本规模的增加，近似解可以任意接近真实最优解。此外，为了保证无人机在离散时间步之间的安全性，我们采用了迭代优化方案[8]，即求解CC‐MIQP问题并检查所得轨迹是否无碰撞。若否，则迭代执行求解与检查过程，直至获得期望的轨迹。因此，这种基于机会约束优化的方法比某些鲁棒方法[28],[29],更为宽松，因为后者通常仅考虑不确定性的最坏情况实现。同时，通过迭代方案，避障在整个时间范围内而非仅限于离散时间步上得以保证。

本文其余部分组织如下。第II‐A节首先介绍基本的CC‐MIQP公式，即仅在离散时间步上保证无人机的安全性。第II‐B节提出一种新颖的基于采样的求解方法来解决基本的CC‐MIQP问题。第II‐C节讨论所获得解的可行性与最优性。第III节给出了获得无碰撞轨迹的迭代方案。第IV节通过一组仿真结果验证了所提出方法的有效性。第V节对全文进行总结。为方便读者阅读，两个命题的证明见附录。

II. 基本的机会约束混合整数二次规划：在离散时间步的避障

为了便于阐述，我们首先介绍基本CC‐MIQP模型，其中避障仅在离散时间步上得到保证。无人机在离散时间步之间的行为将在后续进行研究。

A. 基本的机会约束混合整数二次规划

为简化问题表述，假设无人机位于二维平面中。在每个离散时间步t，无人机的位置表示为x(t) ∈ R²。我们研究有限时间范围T内的无人机轨迹规划问题，即总共考虑T个时间步长，因此每个离散时间步t满足t ∈ T:={0, 1,···, T}。此外，我们将无人机的控制输入记为u(t) ∈ R²。无人机轨迹规划问题的约束条件和目标函数如下所述。

1) 无人机动力学 ：除了无人机位置x(t)，我们还进一步考虑其速度，记为v(t) ∈ R²。无人机的动力学建模为以下线性时不变（LTI）状态空间表示，
$$
[x(t+ 1) v(t+ 1)]= A[x(t) v(t)]+ Bu(t), \tag{1}
$$
其中A ∈ R⁴×⁴和B ∈ R⁴ײ表示离散时间系统矩阵。更具体地，让我们在二维情况下定义无人机的位置和速度为
$$
x(t)=[px(t), py(t)]^T \text{ and } v(t)=[vx(t), vy(t)]^T, \tag{2}
$$
其中，下标x和y分别表示沿x轴和y轴的方向上的位置或速度。本文假设无人机的动力学遵循简单的双积分器模型，即 ˙v(t)= ¨x(t)= u(t)，其中控制输入u(t)表示无人机的加速度。因此，系统矩阵可被离散化为
$$
A= \begin{bmatrix}
1 & 0 & \Delta t & 0 \
0 & 1 & 0 & \Delta t \
0 & 0 & 1 & 0 \
0 & 0 & 0 & 1
\end{bmatrix}, \text{ and } B= \begin{bmatrix}
\Delta t^2/2 & 0 \
0 & \Delta t^2/2 \
\Delta t & 0 \
0 & \Delta t
\end{bmatrix},\tag{3}
$$
注意，这里Δt ∈ R⁺表示两个连续离散时间步之间的采样时间间隔。事实上，这种双积分无人机动力学模型已在许多与无人机轨迹规划相关的现有研究中被广泛采用，例如[20],[30],[31]。

2) 起讫点分配 ：考虑到无人机需要从特定的初始位置xorig ∈ R²出发，并在时间范围结束时到达目的地xdest ∈ R²，因此起讫点分配约束条件如下
$$
x(0)= x_{orig}, \text{ and } x(T)= x_{dest}. \tag{4}
$$

3) 有界规划空间 ：此外，假设无人机的整个轨迹被限制在有界规划空间内。我们将该空间定义为集合K ⊂ R²，通过一组线性不等式约束表示，即K:={x | Lx ≤b}。因此，有界规划空间约束表示为
$$
x(t) ∈ K, ∀t ∈ T. \tag{5}
$$

4) 地理围栏限制 ：假设总共有I个地理围栏，每个地理围栏由一个凸多边形Si表示，其索引为i ∈ I:={1, 2,···, I}。我们对每个凸电子围栏Si建模如下：设Si的中心为一个固定点pi ∈ R²，则每条边（索引为k ∈{1, 2,···, Ki}）对应于地理围栏的一个边界；参见图1中的实线。注意，每条边k可由其到点pi的垂距rik以及相应的角度θki唯一确定，如图1所示。借助垂距rik和角度θki，多边形的每条边界k可被解析表示为，
$$
\begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x - p_i) - r_k^i = 0, ∀k ∈{1, 2,···, Ki}.\tag{6}
$$
因此，由具有Ki边界的凸电子围栏Si定义的限制区域是，
$$
S_i:=\left{x \,\middle|\, \begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x - p_i) < r_k^i, ∀k ∈{1, 2,···, Ki}, x ∈ R^2\right}. \tag{7}
$$
为了保证安全，无人机在规划轨迹时的任何时刻t都必须避开这些限制区域，即
$$
x(t)\notin S_i, ∀i ∈I, t ∈ T. \tag{8}
$$
注释1 : 请注意，尽管大多数现有研究（例如[11],[20],[21],[23],）通过一组普通的线性不等式来建模凸多边形电子围栏/障碍物，即{(ak i)Tx ≤bk i}，但本文采用了一种新表示法(7)，利用了信息(pi, θki, rki)。这样做的一个显著优点是，可以明确地将不确定性距离rik提取到不等式的右侧。这一特性大大简化了后续开发我们的基于采样的求解方法以及后续章节中的理论性能分析。此外，我们还指出，这种新表示法(7)可以对二维平面上的任意凸多边形进行建模。

接下来，我们将不确定性引入电子围栏的模型中。此处，距离rik不再被假定为确定性常数，而是被视为一组遵循已知概率分布的随机变量。如图1所示，我们使用阴影区域表示边界的全部可能情况。根据该定义，(6)式中的实线代表随机距离的均值。此外需要注意的是，图1中仅展示了第k个边界的不确定性，而在我们的问题设定中，所有边界均存在不确定性。随后，通过以下机会约束来保证地理围栏限制，其中包含预设概率阈值 α ∈(0,1)。
$$
Pr(x(t)\notin S_r^i, ∀i ∈I) ≥ α, ∀t ∈ T. \tag{9}
$$
请注意，这里我们使用符号Sri，其定义类似于(7)，用于表示由于随机距离rik’引起的电子围栏的不确定性。Pr(·)是这些随机变量的概率测度。这些机会约束(9)背后的直观含义非常直接：假设电子围栏Sri受到某些不确定性的影响，要求无人机位置在任何随机变量实现情况下都位于Sri之外，这种做法过于保守（甚至不可行）。与这种保守方法不同的是，可以允许发生约束违反，但要求违反约束的概率不超过预设的阈值 1 − α。

值得指出的是，机会约束(9)中的限制条件，即x(t)∉Sri，其决策变量的可行域本质上是非凸的，因为对于rik的任意实现，集合Sri本身是凸的。为了处理这一非凸性问题，我们进一步采用大‐M约束的思想，并基于混合整数规划技术。回顾每个电子围栏由多边形Sri建模，如式(7)所定义。那么，无人机相对于电子围栏Sri的可行域可表示为该多边形的外部，即
$$
C_r^i:=\left{x \,\middle|\, \text{there exists } k ∈{1, 2,···, K_i} \text{ such that } \begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x - p_i) ≥ r_k^i, x ∈ R^2\right}. \tag{10}
$$
为了进一步处理定义Cri中k的存在性条件，我们引入一组二元变量zki ∈{0, 1}，其中k ∈{1, 2,···, Ki}。因此，集合Cri是等价的
$$
C_r^i=\left{x \,\middle|\, \begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x - p_i)+ M(1 - z_k^i) ≥ r_k^i, \sum_{k=1}^{K_i} z_k^i ≥ 1, z_k^i ∈{0, 1}, x ∈ R^2\right}. \tag{11}
$$
其中M是一个足够大的常数。因此，机会约束(9)现在被重写为，
$$
Pr(x(t) ∈ C_r^i, ∀i ∈I) ≥ α, ∀t ∈ T. \tag{12}
$$

5) 控制努力最小化 ：在本文中，我们考虑最小化控制努力的总成本，同时确保无人机的安全。因此，整体优化问题被表述为，
$$
\text{minimize} {{u(t)} {t∈T}} \sum_{t=0}^{T-1} |u(t)|^2, \text{ subject to }(1),(4),(5), \text{ and }(12). \tag{13}
$$
注意，问题（13）是一个包含二次目标函数和机会约束的混合整数规划，我们将其称为带机会约束的混合整数二次规划（CC‐MIQP）问题。接下来，我们将开发一种基于采样的求解方法来解决此类CC‐MIQP问题。

B. 基于采样的求解方法

在进一步开发我们的方法之前，让我们对基本的CC‐MIQP模型（13）做几点说明。首先，应注意机会约束的概率是联合应用于所有电子围栏Si,i ∈I的，如式（12）所示。作为一种更简单的方法，也可以分别考虑每个地理围栏，从而机会约束可表示为，
$$
Pr(x(t) ∈ C_r^i) ≥ α_i, ∀i ∈I, t ∈ T. \tag{14}
$$
在这种情况下，所得到的问题变得更容易求解，因为(14)中的每个单独的机会约束都可以等价地重写为以下确定性形式，
$$
\begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x(t) - p_i)+ M(1 -z_k^i(t)) ≥ q_{i,k}(α_i),\
\sum_{k=1}^{K_i} z_k^i(t) ≥ 1, z_k^i(t) ∈{0, 1}, ∀k ∈{1, 2,···, Ki}, t ∈ T. \tag{15}
$$
这里，qi,k(αi)是单个随机变量rik的αi‐分位数，即Pr(rk ≤ qi,k(αi)) ≥ αi。然而，本文中我们考虑的是联合概率，即需要确保联合概率不低于某个特定阈值α。因此，将存在无穷多种αi的组合方式使其乘积为α，从而考虑任何特定组合仅是一种可能存在保守性的情况。在本文中，我们采用基于采样的方法来处理联合概率，如下所述。其次，请注意另一类方法是利用可行域的凸性，当随机变量的分布具有某些特定性质时[32]–[34]。但在我们的设定中，由于模型中涉及二元变量zki’s，难以保持凸性属性，甚至无法构建非凸机会约束的凸近似。此处，我们考虑对基本的CC‐MIQP问题进行近似，其中随机变量的连续分布被基于蒙特卡洛采样所选定的一些特定样本替代。所得到的近似问题被表述为一个确定性的混合整数二次规划问题，因此可以很容易地由CPLEX[14]和GUROBI[15]等商业求解器处理。更重要的是，如后文所示，该方案还能为原始机会约束问题提供具有可行性保证的有希望的候选解。

为简化符号表示，此处将所有不确定距离rik堆叠在一起，并将其表示为一个维度为K= ∑Ii=1Ki的随机向量r=[r¹₁, r²₁ ···, rᴷ¹₁,···, r¹_I, r²_I,···, rᴷᴵ_I] ∈ RK。假设该随机向量r服从已知的概率分布，其联合累积分布函数（CDF）由下式给出：
$$
F(y):= Pr(r_k^i ≤ y_k^i, ∀k ∈{1, 2,···, Ki},i ∈I).\tag{16}
$$
假设根据均匀分布N个独立样本（也称为场景）被抽取，其分布为U(0,rmax)。记每个样本为r(n)，因此这N个样本的集合为 R={r(1),r(2),···,r(N)}。注意，此处要求每个样本n逐元素满足0 ≤ r(n) ≤ rmax，因为过大的样本对于求解该问题没有意义，这一点我们将在后文看到。实际上，根据电子围栏的中心以及（5）中定义的有界规划空间，可以预先计算出样本上界rmax。

从独立同分布（i.i.d.）的样本出发，文献中提出了多种基于场景的方法[35]–[38],，其中将r的真实分布替换为N个采样场景，从而将机会约束转化为N个独立的确定性约束。这些基于场景的方法背后的思想非常直观：当样本数量足够大，并且所获得的解对所有这些样本都可行时，该解极有可能（概率大于 α）也适用于原始的机会约束问题。然而，如[38],所指出的，由于在原始问题设定中允许以概率 1−α发生约束违反，要求解满足所有采样场景可能会带来较高的保守性。受此启发，我们接下来提出一种相较于现有基于场景方法具有更低保守性的求解方法。我们的核心思想是：不再强制考虑来自集合R的全部N个采样场景，而是有选择地筛选出其中一部分样本，使得满足以下两个条件：第一，每个选定样本单独都能产生对原始机会约束问题具有可行性保证的解；第二，所有选定样本共同能够反映约束违反容忍概率 1 − α。通过这种方式，我们期望构建一种方法，仅需选取一个选定样本即可得到可行解，同时实现最小的保守性。需要注意的是，解的可行性可由前述第一个关于选定样本的条件自动保证。当前的问题转化为如何通过引入第二个条件来实现最不保守的解。

考虑随机变量rik在机会约束的不等式右侧单独出现；参见定义(11)以及注释1。这意味着当距离rik的值较大时，获得可行轨迹（即位于不确定的地理围栏之外）的概率将更高。这一观察促使我们选择具有较大rik(n)值的采样场景r(n)。具体而言，我们通过以下准则确定是否应选择某个采样场景：若采样场景r(n)的联合累积分布函数满足F(r(n)) ≥ α，则称该采样场景被选定。注意，此处我们假设如(16)所示的联合累积分布函数可作为先验知识获得，因为随机向量r的分布是已知的。即使该随机向量可能服从任意的通用分布，我们仍指出其累积分布函数可通过蒙特卡洛仿真[39]提前近似得到。通过遵循此类样本选择准则，将在下一小节中证明每个选定样本均可产生原始机会约束优化的一个可行解，即满足第一个条件。为了考虑第二个条件，我们利用最小化问题的特性。

假设根据准则总共选定了S个样本，则我们将选定的样本子集表示为 R˜S={˜r(1), ˜r(2),···, ˜r(S)} ⊆ R。随后，我们用选定样本R˜S替代机会约束中的原始随机向量r，并在求解问题时要求该样本集合中至少有一个样本被满足。因此，机会约束被以下确定性约束所替代，
$$
x(t) ∈ C_{\tilde{r}(s)}^i, ∃s ∈{1, 2,···, S}, ∀i ∈I, t ∈ T. \tag{17}
$$
为了处理(17)中的存在性条件，我们沿用之前的方法，即引入一组二元变量w(s) ∈{0, 1}，s ∈{1, 2,···, S}。因此，对于t ∈ T，地理围栏约束的显式形式变为，
$$
\begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x(t) - p_i)+ M(1 -z_k^i(t)) ≥ \sum_{s=1}^{S} w(s)\tilde{r} k^i(s),\
\sum {k=1}^{K_i} z_k^i(t) ≥ 1, \sum_{s=1}^{S} w(s)= 1, z_k^i(t), w(s) ∈{0, 1}, ∀i ∈I. \tag{18}
$$
这里，我们使用˜rik(s)表示所选场景˜r(s)的对应分量。接下来，我们通过交替求解以下确定性优化问题来逼近原始CC‐MIQP(13)的解，
$$
\text{minimize} {{u(t)} {t ∈ T}} \sum_{t=0}^{T-1} |u(t)|^2, \text{ subject to }(1),(4),(5), \text{ and }(18). \tag{19}
$$
我们注意到，由于考虑的是一个最小化问题，并且只有一个选定样本会被激活，可以验证该优化倾向于激活具有相对较小˜rik(s)值的选定样本。在此基础上，前述条件中的第二个条件也有望得到满足，因为当距离rik的值较小时，所施加的保守性将更少。

C. 性能分析

值得注意的是，在选定样本子集R˜S的帮助下，原始的机会约束(12)已被一组确定性混合整数约束(18)进行了适当近似。为了验证该近似的有效性，本文通过以下两个命题提供其性能的理论保证。请注意，这两个命题也分别对应于前述选定样本的两个条件。在进入命题之前，我们首先引入两个额外的集合，它们分别由机会约束(12)及其确定性近似(18)定义。我们将由机会约束(12)诱导出的可行集X记为
$$
X:={x | Pr(x ∈ C_r^i, ∀i ∈I) ≥ α}, \tag{20}
$$
并将其由(18)诱导的近似集X˜(R˜S)表示为
$$
\tilde{X}(R_S):=\left{x \,\middle|\, \exists s \in {1,2,\dots,S} \text{ such that } \begin{bmatrix}
\cos(\theta_k^i) & \sin(\theta_k^i)
\end{bmatrix}(x - p_i)+ M(1 -z_k^i) ≥ \tilde{r} k^i(s), \sum {k=1}^{K_i} z_k^i ≥ 1, z_k^i \in {0, 1}, \forall i \in I \right}. \tag{21}
$$

命题1 ：对于任何非空的选定样本集R˜S，即S ≥ 1，有X˜(R˜S) ⊆ X成立。

证明 ：见附录A。

需要指出的是，上述命题1保证了通过求解(19)得到的解对于原始的随机约束优化问题(13)也是可行的。为了进一步研究所获得解的最优性，我们给出下一个命题。

命题2 ：假设样本R是独立同分布的，并设选定样本的数量为S → ∞，则有
$$
\lim_{S→∞} Pr(X ⊆ X˜(R˜S))= 1. \tag{22}
$$

证明 ：见附录B。

总之，考虑到处理概率地理围栏中所涉及的不确定性带来的挑战，本文提出求解近似确定性优化(19)，而非直接求解机会约束问题(13)。根据命题1，理论上可保证所得解对原始CC‐MIQP(13)是可行的。此外，命题2还验证了随着样本规模的增加，近似解可以任意接近真实最优解。

进一步地，在求解确定性优化(19)后，若某个样本对应的二进制变量满足s的w(s)= 1，则将该样本称为活跃场景。该活跃场景的确定有助于识别电子围栏中在所得解中反映出来的实际限制区域。

在明确了机会约束混合整数二次规划（CC‐MIQP）的公式化表述并给出了求解方法后，我们现在提供一个通过初步仿真获得的示例。考虑从原点xorig=[0, 0]^T到目的地xdest=[10, 10]^T的轨迹规划问题。有限时间范围被均匀离散为T= 20个时间步长，时间间隔为Δt= 1。假设在二维平面上存在两个电子围栏，分别以点p1=[3, 4]^T和p2=[7, 6]^T为中心。其受限区域被建模为两个圆盘，具有不确定性半径(r1,r2)。

如图2(a)所示。在我们的仿真中，假设随机变量对(r1,r2)服从均值为μ=[1, 1]^T、协方差为Σ=[0.01, 0; 0,0.01]的高斯分布N(μ,Σ)。机会约束的预设概率阈值设为α= 0.9。此外，无人机动力学遵循式(1)和(3)中规定的线性时不变系统。注意，我们使用GUROBI 9.0软件包求解得到的确定性混合整数规划。

为了适应第II‐A节中建立的CC‐MIQP模型，我们首先用一个16边形来近似类圆盘地理围栏，如图2(b)所示。注意，每个多边形近似可表示为式(7)中定义的集合Si，其角度为θki=(2k −1)π/16, k ∈{1, 2,···, 16}。接下来，为了应用基于采样的求解方法，我们在随机变量(r1,r2)上抽取N= 1000i.i.d.样本，并根据联合累积分布函数选取其中的S= 614个样本。所得到的无人机最优轨迹如图2(a)所示。我们注意到，以实线表示的边界半径为r1= r2= 1，即分布的均值；而以虚线表示的边界代表活跃场景，其半径分别为r1= 1.21和r2= 1.14。由于机会约束要求满足的概率至少为α= 0.9，因此可以容易地验证地理围栏的实际半径预期大于均值。这一点与图2(a)中的观察结果一致。

第三节 迭代式CC‐MIQP：实现无碰撞轨迹

回想一下，在基本CC‐MIQP模型中，无论是在原始机会约束(12)还是其确定性近似(18)中，无人机的安全性仅在离散时间步t ∈ T处得到保证。然而，在这种情况下，由于潜在的离散化问题，无人机在时间步长之间仍可能面临风险。在本节中，我们的目标是解决此类离散化问题，并获得完全无碰撞轨迹，即不仅在离散时间步上避开电子围栏，而且在整个时间范围内均保持无碰撞。

在继续之前，让我们首先通过两个简单的示例来说明无人机轨迹潜在的离散化问题。为了便于说明，此处我们假设两个示例中的每个电子围栏Si都具有确定性边界。如图3所示，电子围栏是两个圆盘，分别以[3, 4]^T和[7, 6]^T为中心，并具有固定半径ri= 1。有限时间范围被均匀离散化为T1= 10和T2= 20个时间步长，时间间隔为Δt= 1。从图3(a)和图3(b)中可以观察到，尽管无人机在每个离散时间步的位置都在圆盘之外，但所得到的轨迹仍然是不安全的，因为它们均穿越了电子围栏的受限区域。为了进一步解决上述离散化问题，我们在此采用文献[8]中建议的迭代方案思想。该方法的指导原则是将控制时间步长和避障时间步长分开考虑，从而可以独立地迭代添加避障时间步长，直到获得无碰撞的轨迹。

为了区分无人机控制和避障的离散化，我们用Tu表示控制时间步长的数量，用Ts表示避障时间步长的数量。注意，虽然控制时间步长可以如前所述进行指定，但避障时间步不再假设为均匀离散化。相反，我们用ts[k]（其中k ∈{1, 2,···, Ts}）表示避障时间步，并用Ts:={ts[k],k ∈{1, 2,···, s}}表示这些时间步的集合。

考虑这两种离散时间步后，基本CC‐MIQP模型变为，
$$
\text{minimize} {u(t)} \sum {t=0}^{T_u-1} |u(t)|^2
$$
subject to(1),(4),(5), and
$$
Pr(x(t_s[k]) ∈ C_r^i) ≥ α, ∀t_s[k] ∈ T_s.\tag{23}
$$
在公式(23)中，我们通过确定控制时间步长t ∈{0, 1,···, Tu−1}处的输入u(t)来控制无人机，同时确保无人机在避障时间步ts[k] ∈ Ts的安全。需要注意的是，当无人机的动力建模为线性时不变系统(1)时，其在任意时刻τ（包括位置x(τ)和速度v(τ)）的状态均可精确计算如下。不失一般性，我们假设时间τ落在区间[t,t+ 1]内，并且已知无人机在时刻t的状态，则根据无人机的动学特性，可得到x(τ)和v(τ)
$$
x(τ)= x(t)+(τ − t)v(t)+ 0.5(τ − t)^2 u(t); \tag{24}
$$
$$
v(τ)= v(t)+(τ − t)u(t). \tag{25}
$$
因此，一旦避障时间步长集合Ts确定，即可通过应用第II‐B节中提出的求解方法来解决优化问题(23)。

我们现在可以着手解决上述离散化问题。其思路是：在开始时仅使用少量避障时间步长来求解优化问题(23)，然后检查所得到的轨迹是否为无碰撞轨迹。如果不是，则通过安全检查过程识别出不安全穿越，并针对每次不安全穿越i计算相应的时间间隔[tli, tui]。注意，此处tli和tui分别表示无人机进入和离开电子围栏的时间。接下来，我们将新的避障时间步长tnews=(til+tui)/2加入集合Ts中，从而增强先前的优化问题，并重新求解该增强后的优化问题。此过程重复进行，直到获得完全无碰撞轨迹为止。综上所述，我们将该迭代方案总结为以下算法1，并对该算法提供若干说明。

算法1 基于迭代CC‐MIQP的轨迹规划
数据：初始化避障时间步长集合如Ts= ∅所示。求解CC‐MIQP问题(23)，并获得轨迹。

while所获得的轨迹不是无碰撞的do
(S.1) 对于每次不安全穿越i，计算穿越时间区间[til, tui];
(S.2) 通过向基本CC‐MIQP模型中添加新的时间步长到Ts中，进行扩展，
      Ts ←− Ts ⋃{ (til + tui) / 2}, ∀i; \tag{26}
(S.3) 使用基于采样的求解方法求解扩展后的CC‐MIQP问题，以及获得新的轨迹；
(S.4) 检测不安全穿越，并继续。
end

备注2 ：需要注意的是，不安全穿越的检测和穿越时间区间的计算可以通过多种方式高效完成。首先，如[8],所建议的，可以使用简单的二分法程序来识别碰撞时刻，进而计算穿越时间区间。尽管该过程仍然依赖于无人机连续时间动力学的离散化，但我们认为，由于此处不涉及优化问题，因此通过增加采样时刻的数量来实现更精细的离散化是可行的。此外，在不对连续时间动力学进行离散化的情况下，也可以通过解析求解涉及电子围栏边界模型和无人机动力学的非线性方程组来确定碰撞时刻。我们注意到，与简单的二分法检测相比，这种解析方法的计算成本可能显著更高；然而，如果通过求解非线性方程组来找到碰撞，则完全不存在离散化问题。因此，算法1生成的无人机轨迹可保证为无碰撞轨迹。

注释3 : 我们强调，在我们的概率设置中，不安全穿越的检测与计算也受到电子围栏不确定性的影响，因为它们高度依赖于边界的随机距离。回顾用于求解CC‐MIQP问题的基于采样的求解方法，仅存在一个决定所获得轨迹的活跃场景，该场景由优化求解后二元变量w(s)’s的值指示。因此，我们可以通过仅考虑所获得的活跃场景来实现对不安全穿越的检测以及穿越时间区间的计算。

IV. 数值示例

在本节中，我们首先通过展示两种不同的无碰撞轨迹规划案例，验证了所提出的迭代CC‐MIQP方法。此外，我们评估了概率阈值α和样本规模N的选择对方法性能的影响。

A. 无碰撞轨迹的演示

我们采用与第II‐C节类似的仿真设置。假设无人机需要根据线性时不变动力学(1)从原点xorig=[0, 0]^T规划轨迹至目的地xdest=[10, 10]^T，系统参数如(3)所示。控制时间步长被均匀离散化为Tu= 20，时间间隔为Δt= 1。在本仿真中，考虑两种不同的地理围栏场景：i) 在简单情况下，如图4(a)所示，两个矩形地理围栏的中心分别位于二维平面上的点p1=[2.0, 5.5]^T和p2=[7.0, 6.0]^T；ii) 在更复杂的情况下，如图4(b)所示，八个随机生成位置pi’s和角度θki’s的多边形被模拟为地理围栏。更具体地，在简单情况i)中，我们使用hi和wi表示点pi到水平和垂直边界之间的距离，并假设每对随机变量(hi, wi)服从联合高斯分布N(μi,Σi)，其均值为μ1=[3.5, 0.5]^T, μ2=[2.0, 0.4]^T，协方差矩阵为Σ1=[0.01, 0; 0; 0.01], Σ2=[0.01, 0; 0; 0.01]。在情况ii)中，为了便于演示，我们假设每个地理围栏i的所有Ki个边界具有相同的距离ri，即rik=ri,∀1 ≤k ≤Ki；此外，每个ri服从高斯分布N(μi, σ²i)，其中μi是随机生成的，σ²i= 0.01,∀i ∈ I。对于两种情况，机会约束的概率阈值均设为α= 0.9。为了应用基于采样的求解方法来解决机会约束优化问题，我们对两种情况中的每种随机变量分别抽取N= 1000个独立同分布样本。根据其联合累积分布函数的值，最终选择S= 332个样本用于情况i)，S= 170个样本用于情况ii)。注意，根据第II‐B节的分析，每一个样本都必须能够提供可行解。

两种情况下的轨迹规划结果分别如图4(a)和图4(b)所示。在这两幅图中，黑线表示将避障步长时间与控制时间步长视为相同的情况，即Ts={1, 2,···, Tu − 1}。在这种情况下可以观察到，尽管无人机在每个离散的控制时间步长（以黑色三角形标记）都是安全的，但其轨迹仍存在潜在的碰撞风险。

穿越一些电子围栏，因此无人机实际上处于风险之中。这一观察结果也与第III节开头介绍的两个示例一致。接下来，为了实现完全无碰撞轨迹，我们执行所提出的算法1，即通过添加避障时间步数来迭代增强优化模型。红线绘制了获得的轨迹，蓝色加号标记了每个避障时间步数时的无人机位置。从两种情况下的放大图可以看出，无人机刚好在电子围栏周围通过而未发生穿越。可以确认，这两条轨迹均为无碰撞轨迹。还应指出的是，对于每个电子围栏，实线表示由均值距离确定的边界，虚线表示主动场景中的边界。回顾不安全穿越检测仅考虑主动场景中的边界；参见注释3，我们得出结论：仿真结果符合预期。

最后，为了评估在机会约束下所获得轨迹的有效性，我们进一步通过对随机变量的采样数据集进行解的验证。具体而言，对于简单情况i)，我们通过在hi和wi上独立同分布地采样10⁵个场景来生成数据集，这些场景服从高斯分布N(μ1,Σ1)和N(μ2,Σ2)。接着，通过检查所获得的轨迹在每个采样场景下是否无碰撞，来测试该轨迹。为了确保在机会约束α= 0.9下的可行性，期望该轨迹在超过90%的采样场景中均无碰撞。事实上，我们观察到，在10⁵个样本中有95227个样本通过了安全验证。同样的过程也应用于情况ii)，观察到所获得的轨迹在10⁵个样本中有98961个样本是无碰撞的。这些结果证实了两个解在机会约束下的可行性。此外，考虑到情况ii)涉及的选定样本数（S= 170）少于情况i)中的选定样本数S= 332，因此在验证中所得轨迹更为保守是合理的（98.96%对比95.23%）。

B. 概率阈值选择对性能的影响

在本小节中，我们比较了算法在不同概率阈值α选择下的性能表现。对于第IV‐A节中的情况i）和ii），我们采用相同的仿真设置，仅将α的值从0.5变化到0.9。具体而言，我们关注所获得解的以下性能表现：i）最优目标函数值；以及ii）选定样本的数量S。考虑到较小的概率阈值α本质上意味着对约束违反的容忍度更高，因此理论上，当选择较小的α时，预期会观察到：i）获得更小的最优目标函数值，因为无人机具有更大的可行域来规划其轨迹；以及ii）更多的样本被选中，因为选择准则F(r(n)) ≥ α变得不那么严格。事实上，如图5所示的仿真结果证实了上述两种推断。请注意，此处报告的数据是通过每种情况运行了20次独立仿真，实线表示平均值。此外，从图5(a)和图5(c)中也可以看出，随着概率阈值α的增加，最优值的方差总体上有所减小。这主要是因为当指定较大的α时，所得解的变化可能性较小。然而，应注意在图5(c)中观察到了两个例外情况。首先，当α= 0.5时，最优解的方差异常小，因为无人机总能找到一条轨迹，无碰撞地穿过右上方区域的电子围栏（见图4(b)）。在这种情况下，所得到的轨迹非常接近于不考虑任何电子围栏时生成的轨迹。这解释了为何均值和方差都很小。其次，还可以观察到当α= 0.7时，方差相当大。事实上，这可以用与上述相同的原因来解释。当α= 0.7时，无人机有时可以穿越右上方的地理围栏，从而获得较小的最优值；但更多时候，它必须绕过这些地理围栏的聚集区域，导致相对较优的最优值较大。

C. 样本规模选择的性能

此外，我们评估了样本规模N的选择对轨迹规划结果的影响。根据第II‐B节的分析，随着样本规模的增加，基于采样的方法生成的解将更接近于直接求解随机问题所得到的真实最优解。为了验证这一推论，本文在相同的概率α= 0.8下，但采用不同的样本规模N，重新运行了前述仿真。与第IV‐B节的分析类似，我们关注所获得解的以下性能指标：i) 最优目标函数值；ii) 选定样本数S；此外还有iii) 算法的执行时间；iv) 为获得无碰撞轨迹而迭代添加的避障时间步数。需要注意的是，所需避障时间步数并不一定与迭代次数成正比，因为在算法1的同一次迭代中可能会添加多个时间步长。同样地，表I和表II提供了通过蒙特卡洛仿真运行所得结果的平均值，每个单元格中的第一个数字表示均值，而第二个是标准差。 )
)

从表中可以观察到，当样本规模N增加时：i）获得的目标函数值更小；ii）选定样本更多；iii）算法的执行时间更长；iv）所需的避障时间步数保持在同一水平。关于ii），其解释实际上是显而易见的，因为事先抽取了更多的样本。对于iii），由于计算过程中涉及大量选定样本，算法执行时间无疑会增加。观察结果i）本质上有助于验证命题2的正确性：当样本规模N增加时，机会约束的tighter bound更有可能被采样到。此外，表I和表II中的数据也反映了本方法在效率与精度之间的权衡：更大的样本规模有助于找到更接近真实最优解的解，但另一方面也增加了计算复杂度。有趣的是，在两种情况下均观察到，当样本规模从500增加到1000时，最优目标函数值的下降幅度明显大于从1000增加到3000或5000时的情况。这一现象表明，过大的样本规模可能带来的改善有限，尤其是在考虑计算复杂度的情况下。然而需要指出的是，上述结论仅基于两个相对简单的电子围栏情况进行分析得出。原则上，随着地理围栏数量的增加，为了获得机会约束的紧致边界，需要更多的样本以充分揭示不确定性环境。

五、结论

本文研究了带有概率地理围栏的无人机轨迹规划问题。提出了一种基本的CC‐MIQP问题，以描述地理围栏的不确定性边界，并处理可行集的非凸性。提出了一种新的基于采样的求解方法来解决该基本的CC‐MIQP问题。证明了通过求解近似确定性优化得到的解对原始CC‐MIQP问题是可行的。此外，随着样本规模的增加，近似解可以任意接近真实最优解。此外，为了解决潜在的离散化问题，进一步发展了迭代CC‐MIQP方法，有助于获得完全无碰撞轨迹。数值仿真最终验证了所提方法的有效性。

未来的研究方向可以从以下三个方面考虑。首先，由于本文提出的轨迹规划方法本质上是离线的，下一步是通过结合无人机在不确定性环境中的实时测量以及对预规划轨迹的实时调整，将我们的算法适应于在线框架。第二，考虑到概率阈值α对轨迹规划结果具有关键影响，正如仿真所示，进一步的问题是如何将α引入优化模型本身，从而确定α的最佳选择。第三，由于本文在规划无人机轨迹时仅考虑最小化总控制effort，其他目标函数也可能具有研究价值，例如最短飞行时间或最短轨迹长度。