一、杜鵑鳥鯰魚優化算法

杜鵑鳥鯰魚優化（Cuckoo Catfish Optimizer，CCO）算法模擬了杜鵑鳥鯰魚的搜索、捕食和寄生慈鯛行為。該算法的早期迭代側重于執行多維包絡搜索策略和壓縮空間策略，并結合輔助搜索策略來有效限制慈鰾的逃逸空間。此階段確保對解決方案空間進行廣泛探索。在迭代的中間階段，該算法采用過渡策略促進從勘探到開發的平滑過渡，賦予了算法一定的勘探能力和開發能力。在后期階段，該算法使用混沌捕食機制在慈鳷周圍制造干擾，以提高對最優解的利用。在整個優化過程中，整合了個體的引導、寄生和死亡機制，讓個體能夠實時調整位置，提高整體收斂精度。

參考文獻
[1]Tian-Lei Wang, Shao-Wei Gu, Ren-Ju Liu, Le-Qing Chen, Zhu Wang*, Zhi-Qiang Zeng*. Cuckoo Catfish Optimizer: A New Meta-Heuristic Optimization Algorithm. Artificial Intelligence Review,2024

二. 無人機路徑規劃數學模型

2.1 路徑最優性

為了提高無人機的操作效率，規劃的路徑需要在特定的應用標準下達到最優。在我們的研究中，主要關注空中攝影、測繪和表面檢查，因此選擇最小化路徑長度作為優化目標。由于無人機通過地面控制站（GCS）進行控制，飛行路徑 $X_i$ 被表示為無人機需要飛越的一系列 $n$ 個航路點的列表。每個航路點對應于搜索地圖中的一個路徑節點，其坐標為 $P_{ij}=(x_{ij},y_{ij},z_{ij})$。通過表示兩個節點之間的歐幾里得距離為 $| \overrightarrow{P_{ij}P_{i,j+1}} |，與路徑長度相關的成本 $F_1$ 可以計算為：

$$F_1(X)=\sum _{j=1}^{n?1}\Vert \overrightarrow{P_{ij}{P}_{i,j+1}}\Vert$$

2.2 安全性和可行性約束

除了最優性之外，規劃的路徑還需要確保無人機的安全操作，引導其避開操作空間中可能出現的威脅，這些威脅通常由障礙物引起。設 $K$ 為所有威脅的集合，每個威脅被假設為一個圓柱體，其投影的中心坐標為 $C_k$，半徑為 $R_k$，如下圖 所示。

對于給定的路徑段 $\Vert \overrightarrow{P_{ij}{P}_{i,j+1}}\Vert$，其相關的威脅成本與它到 $C_k$ 的距離 $d_k$ 成正比。考慮到無人機的直徑 $D$ 和到碰撞區域的危險距離 $S$，威脅成本 $F_2$ 在障礙物集合 $K$ 上計算如下：

$$F_2(X_i)=\sum _{j=1}^{n?1}\sum _{k=1}^K{T}_k(\overrightarrow{P_{ij}{P}_{i,j+1}}),$$

其中

$$T_k(\overrightarrow{P_{ij}{P}_{i,j+1}})=?\begin{array}{ll}0,& \text{if?}{d}_k>S+D+R_k\\ (S+D+R_k)?d_k,& \text{if?}D+R_k<d_k\le S+D+R_k\\ \infty ,& \text{if?}{d}_k\le D+R_k\end{array}$$

在操作過程中，飛行高度通常被限制在給定的最小和最大高度之間，例如在調查和搜索應用中，需要相機以特定的分辨率和視場收集視覺數據，從而限制飛行高度。設最小和最大高度分別為 $h_{\text{min}}$ 和 $h_{\text{max}}$。與航路點 $P_{ij}$ 相關的高度成本計算為：

$$H_{ij}=\{\begin{array}{ll}|h_{ij}?\frac{h_{\text{max}}+h_{\text{min}}}{2}|,& \text{if?}{h}_{\text{min}}\le h_{ij}\le h_{\text{max}}\\ \infty ,& \text{otherwise}\end{array}$$

其中 $h_{ij}$ 表示相對于地面的飛行高度，如下圖所示。

可以看出，$H_{ij}$ 保持平均高度并懲罰超出范圍的值。對所有航路點求和得到高度成本：

$$F_3(X)=\sum _{j=1}^n{H}_{ij}$$

平滑成本評估轉彎率和爬升率，這對于生成可行路徑至關重要。如下圖 所示。

轉彎角 ${?}_{ij}$ 是兩個連續路徑段 $\overrightarrow{P_{ij}^{\prime }{P}_{i,j+1}^{\prime }}$ 和 $\overrightarrow{P_{i,j+1}^{\prime }{P}_{i,j+2}^{\prime }}$ 在水平面 Oxy 上的投影之間的角度。設 $\overrightarrow{k}$ 是 z 軸方向的單位向量，投影向量可以計算為：

$$\overrightarrow{P_{ij}^{\prime }{P}_{i,j+1}^{\prime }}=\overrightarrow{k}\times (\overrightarrow{P_{ij}{P}_{i,j+1}}\times \overrightarrow{k})$$

因此，轉彎角計算為：

$${?}_{ij}=\arctan \left(\frac{\Vert \overrightarrow{P_{ij}^{\prime }{P}_{i,j+1}^{\prime }}\times \overrightarrow{P_{i,j+1}^{\prime }{P}_{i,j+2}^{\prime }}\Vert }{\overrightarrow{P_{ij}{P}_{i,j+1}^{\prime }}?\overrightarrow{P_{i,j+1}^{\prime }{P}_{i,j+2}^{\prime }}}\right)$$

爬升角 ${\psi }_{ij}$ 是路徑段 $\overrightarrow{P_{ij}{P}_{i,j+1}}$ 與其在水平面上的投影 $\overrightarrow{P_{ij}^{\prime }{P}_{i,j+1}^{\prime }}$ 之間的角度，由下式給出：

$${\psi }_{ij}=\arctan \left(\frac{z_{i,j+1}?z_{ij}}{\Vert \overrightarrow{P_{ij}^{\prime }{P}_{i,j+1}^{\prime }}\Vert }\right)$$

然后，平滑成本計算為：

$$F_4(X)=a_1\sum _{j=1}^{n?2}{?}_{ij}+a_2\sum _{j=1}^{n?1}|{\psi }_{ij}?{\psi }_{j?1}|$$

其中 $a_1$ 和 $a_2$ 分別是轉彎角和爬升角的懲罰系數。

2.3 總體成本函數

2.3.1 單個無人成本計算

考慮到路徑 $X$ 的最優性、安全性和可行性約束，第$i$ 個無人機總體成本函數可以定義為以下形式：

$$f_i(X)=\sum _{k=1}^4{b}_k{F}_k(X_i)$$

其中 $b_k$ 是權重系數，$F_1(X_i)$ 到 $F_4(X_i)$ 分別是路徑長度、威脅、平滑度和飛行高度相關的成本。決策變量是 $X$，包括 $n$ 個航路點 $P_{ij}=(x_{ij},y_{ij},z_{ij})$ 的列表，使得 $P_{ij}\in O$，其中 $O$ 是無人機的操作空間。根據這些定義，成本函數 $F$ 是完全確定的，可以作為路徑規劃過程的輸入。

2.3.2 多無人機總成本計算

若共有$m$ 個無人機，其總成本為單個無人機成本和，計算公式如下：
$$fitness(X)=\sum _{i=1}^m{f}_i(X)$$
參考文獻：
[1] Phung M D , Ha Q P .Safety-enhanced UAV path planning with spherical vector-based particle swarm optimization[J].Applied Soft Computing, 2021(2):107376.DOI:10.1016/j.asoc.2021.107376.

三、部分代碼及結果

close all
clear
clc
% dbstop if all error
pop=100;%種群大小（可以修改）
maxgen=200;%最大迭代（可以修改）%% 模型建立
model=Create_Model();
UAVnum=4;%無人機數量(可以修改)  必須與無人機的起始點保持一致%% 初始化每個無人機的模型
for i=1:UAVnumModelUAV(i).model=model;
end%% 第一個無人機 起始點
start_location = [120;200;100];
end_location = [800;800;150];
ModelUAV(1).model.start=start_location;
ModelUAV(1).model.end=end_location;
%% 第二個無人機 起始點
start_location = [400;100;100];
end_location = [900;600;150];
ModelUAV(2).model.start=start_location;
ModelUAV(2).model.end=end_location;
%% 第三個無人機 起始點
start_location = [200;150;150];
end_location =[850;750;150];
ModelUAV(3).model.start=start_location;
ModelUAV(3).model.end=end_location;
%% 第四個無人機 起始點
start_location = [100;100;150];
end_location = [800;730;150];
ModelUAV(4).model.start=start_location;
ModelUAV(4).model.end=end_location;
%% 第5個無人機 起始點
% start_location = [500;100;130];
% end_location = [850;650;150];
% ModelUAV(5).model.start=start_location;
% ModelUAV(5).model.end=end_location;
% %% 第6個無人機 起始點
% start_location = [100;100;150];
% end_location =   [800;800;150];
% ModelUAV(6).model.start=start_location;
% ModelUAV(6).model.end=end_location;

五個無人機：

四、完整MATLAB代碼見下方名片