无人机路径规划问题（Unmanned Aerial Vehicle Path Planning Problem, UAV-PP）旨在为多架无人机规划最优飞行路线，使其从各自起点出发，遍历所有指定目标点后返回，最小化总飞行距离或时间。该问题可建模为多旅行商问题（Multiple Traveling Salesman Problem, MTSP），属于 NP-难问题。本文提出混合离散浣熊优化算法（Hybrid Discrete Coati Optimization Algorithm, HDCOA）对该问题进行求解。

1. 问题定义与符号说明

考虑一个完全无向图 $G=(V,A)$，其中 V={1,2,…,n}V = \{1, 2, \dots, n\}V={1,2,…,n} 为目标点（或节点）集合，$A$ 为可行路径边集。每条边 $(i,j)$ 的权重为欧氏距离：

$$d_{ij}=\sqrt{(x_i-x_j{)}^2+(y_i-y_j{)}^2}\text{\hspace{1em}(1)}$$

设共有 $m$ 架无人机，每架无人机从各自起点（可相同或不同）出发，访问若干目标点后返回起点。路径规划的目标是将 $V$ 划分为 $m$ 个非空子集 {S1,S2,…,Sm}\{S_1, S_2, \dots, S_m\}{S1​,S2​,…,Sm​}，并为每个子集找到一条哈密顿回路，使得所有回路的总长度最小：

$$f=\sum _{t=1}^m\left(d_{{\text{depot}}_t,{\text{first}}_t}^{(t)}+\sum _{k=1}^{n_t-1}{d}_{k,k+1}^{(t)}+d_{{\text{last}}_t,{\text{depot}}_t}^{(t)}\right)\text{\hspace{1em}(2)}$$

其中 $n_t$ 为第 $t$ 架无人机需要访问的目标点数量，${\text{depot}}_t$ 为其起点。为平衡各无人机工作量，还需最小化路径长度的差异：

$$f_d=\sum _{i=1}^m\sum _{j=i+1}^m|R_i-R_j|\text{\hspace{1em}(3)}$$

式中 $R_i$ 为第 $i$ 架无人机的总路径长度。

2. 目标点分配与问题分解策略

2.1 扇形划分方法

为将目标点均匀分配给 $m$ 架无人机，采用扇形划分方法。以所有目标点的中心为原点 $O$，选取一条参考线（如水平线），计算每个目标点与原点连线的夹角 ${\theta }_i$：

$${\theta }_i=\arctan \left(\frac{y_i-O_y}{x_i-O_x}\right)\text{\hspace{1em}(4)}$$

按夹角大小对目标点排序后，依次划分给各无人机，确保每架无人机所辖目标点在空间上均匀分布，从而间接实现路径长度均衡。

2.2 K-近邻（KNN）子问题划分

当某架无人机需访问的目标点数量超过阈值 $K_{\text{th}}=150$ 时，采用 KNN 算法进一步划分。对于目标点 $p$，其 KNN 分类规则为：

$$\text{class}(p)=\arg \underset{c}{\max }\sum _{q\in \text{KNN}(p)}1(q\in c)\text{\hspace{1em}(5)}$$

其中 $\text{KNN}(p)$ 表示与 $p$ 距离最近的 $k$ 个目标点（本文取 $k=150$），$1(\cdot )$ 为指示函数。通过该划分，将大规模问题降维为多个小规模子问题分别求解，降低计算复杂度。

3. 离散化浣熊优化行为建模

原始 COA 中的捕猎（探索）与逃逸（开发）行为被重新设计为离散路径操作，用于更新无人机的飞行路线。

3.1 捕猎阶段（探索阶段）

原始 COA 中，浣熊通过攀爬树木恐吓鬣蜥，使其随机坠落。对应的连续更新公式为：

$$X_i^{p1}:x_{i,j}^{p1}=x_{i,j}+r\cdot (Iguan{a}_j-I\cdot x_{i,j}),i=1,2,\dots ,\lfloor \frac{N}{2}\rfloor \text{\hspace{1em}(6)}$$

在 HDCOA 中，该公式被重写为离散形式：

$$X_i^{p1}=x_i\oplus \left(Q\otimes Iguana\right),i=1,2,\dots ,\frac{N}{2}\text{\hspace{1em}(7)}$$

其中：

• $Iguana$ 为当前全局最优路线（目标点序列）；
• $\otimes$ 表示从 $Iguana$ 中随机选取一个子序列，选取比例为 $Q$；
• $\oplus$ 表示将选取的子序列插入到当前解 $x_i$ 中，替换相应长度的子段；
• $Q$ 随迭代次数单调递减：

$$Q=1-e^{\frac{1-t}{T}},t=1,2,\dots ,T\text{\hspace{1em}(8)}$$

式中 $t$ 为当前迭代次数，$T$ 为最大迭代次数。该机制保证了算法在早期更多保留全局最优路线的结构，后期则逐步增强局部扰动。

3.2 Swap 与 Shift 算子

为实现离散空间中的位置更新，定义以下两个算子：

• Swap 算子：随机选择两个不同位置 $u$ 和 $v$，交换其目标点：

$$\text{swap}(x,u,v):x[u]\leftrightarrow x[v]\text{\hspace{1em}(9)}$$

• Shift 算子：随机选择一个子序列 $[a,b]$ 和一个插入位置 $c$（$c\notin [a,b]$），将该子序列移动到 $c$ 之后：

$$\text{shift}(x,[a,b],c):x=x[1:a-1]\oplus x[b+1:c]\oplus x[a:b]\oplus x[c+1:\text{end}]\text{\hspace{1em}(10)}$$

在捕猎阶段，前 $N_1$ 个个体使用式 (7) 进行更新，其余个体使用 Swap 算子进行更新。

3.3 逃逸阶段（开发阶段）

对应浣熊被捕食者追赶时的随机逃逸行为。该阶段所有个体均使用 Shift 算子进行更新：

$$X_i^{p2}=\text{shift}(x_i,[a,b],c)\text{\hspace{1em}(11)}$$

若新路线的适应度（总路径长度）优于原路线，则接受更新，否则保留原解。该操作增强了算法的局部搜索能力，有助于精细化调整无人机路径。

4. 局部最优逃逸机制：Destroy & Rebuild (D&R) 算子

为防止算法陷入局部最优，本文提出了 D&R 算子。其核心思想是破坏当前劣质路线中的长边结构，并重建更优的子路径。

4.1 破坏阶段

选择当前种群中适应度最差的路线 $x_{\text{worst}}$。计算该路线中所有边的长度，并找出最长的 $d$ 条边。其中 $d$ 由下式确定：

$$d=\lfloor \frac{\text{point\_number}}{k}\rfloor \text{\hspace{1em}(12)}$$

这里 $k$ 为 KNN 算法中的邻居数（本文取 $k=150$）。将这 $d$ 条边从路线中移除，从而将原路线分割为 $d$ 个独立的子段。

4.2 重建阶段

对每个子段分别求解其最优哈密顿回路（可采用贪心或局部搜索方法），然后将这 $d$ 个回路按顺序连接，形成新的完整路线 $x_{\text{worst}}^{\prime }$。最后用 2-opt 算子对新路线进行局部优化。

5. 2-opt 路径优化算子

2-opt 算子用于消除路线中的交叉边，这在无人机实际飞行中可避免不必要的转弯和绕路。对于路线中的两条边 $(i,i+1)$ 和 $(j,j+1)$（$i<j$），若它们相交，则反转子路径 $[i+1,j]$：

$$x=x[1:i]\oplus \text{reverse}(x[i+1:j])\oplus x[j+1:\text{end}]\text{\hspace{1em}(13)}$$

该操作可保证路线不自交，从而缩短总飞行距离。

6. 算法整体流程伪代码

Algorithm: HDCOA for Unmanned Aerial Vehicle Path Planning
Input: 目标点集合 V，无人机数量 m，种群规模 N，最大迭代次数 T
Output: 每架无人机的最优飞行路线 {route_1, ..., route_m}

1: 使用扇形划分方法将 V 划分为 m 个子集，分配给各无人机
2: for each 无人机 t = 1 to m:
3:    使用 KNN 生成 NI 个初始路线，其余 N-NI 个随机初始化
4:    for iter = 1 to T/m:
5:        更新当前全局最优路线 Iguana
6:        // 捕猎阶段（探索）
7:        for i = 1 to N1:
8:            X_i_new = X_i ⊕ (Q ⊗ Iguana)   // 式(7)
9:            if fitness(X_i_new) < fitness(X_i) then X_i = X_i_new
10:       for i = N1+1 to N:
11:           X_i_new = swap(X_i)            // 式(9)
12:           if fitness(X_i_new) < fitness(X_i) then X_i = X_i_new
13:       // 逃逸阶段（开发）
14:       for i = 1 to N:
15:           X_i_new = shift(X_i)           // 式(10)
16:           if fitness(X_i_new) < fitness(X_i) then X_i = X_i_new
17:       // 局部最优检测与逃逸
18:       if 最优路线连续 50 代未改善 then
19:           对最差路线 x_worst 执行 D&R 算子  // 式(12)
20:           对 x_worst 执行 2-opt 优化     // 式(13)
21:       end if
22:    end for
23: end for
24: 输出所有无人机的最优飞行路线

7. 时间复杂度分析

HDCOA 求解无人机路径规划的时间复杂度由以下部分构成：

• 扇形划分：$O(n)$；
• 距离矩阵计算：$O(n^2)$；
• 种群初始化：$O(N)$；
• 每代操作（swap / shift）：$O(N\cdot n)$；
• D&R 与 2-opt 算子：$O(n^2)$。

因此，总时间复杂度为：

$$O\left(T\cdot N\cdot n+n^2\right)\text{\hspace{1em}(14)}$$

其中 $T$ 为总迭代次数，$N$ 为种群规模，$n$ 为目标点数量。该复杂度与主流离散元启发式算法相当，能够处理大规模无人机路径规划问题（目标点数量可达数万个）。

8. 部分MATLAB 及结果

clear all
close all
clc

%% ========================================== Initialize ==============================================================

N = 30;%种群
T = 500;%最大迭代数
m = 8;  % 机器人个数
[Dimension,data,NodeWeight,Name] = GetTSPData('./tsp/a280.tsp');%导入数据集

city_coord = data(:,2:3); % City coordinates
city_num = size(city_coord,1); % Number of cities

k = ceil(city_num/m) ;  % Divide the number of cities equally

9. 完整MATLAB见下方名片