基于龙卷风优化算法（ Tornado Optimizer with Coriolis force ，TOC) 的多个无人机协同路径规划（可以自定义无人机数量及起始点），MATLAB代码-优快云博客

一、龙卷风优化算法

龙卷风优化算法（ Tornado Optimizer with Coriolis force ，TOC) 是2025年提出的一种新型的基于自然启发的元启发式算法，其灵感来源于自然界中龙卷风的形成和演化过程。龙卷风的形成是一个复杂的自然现象，通常从雷暴或风暴中发展而来，并在地面上形成强烈的旋转气流。TOC 算法通过模拟龙卷风的形成、发展和消散过程，将其转化为优化问题中的搜索过程。
在这里插入图片描述

龙卷风的形成和科里奥利力

龙卷风的形成过程可以分为以下几个阶段：

初始阶段：在风暴系统中，风速和风向的变化会形成旋转效应，这种旋转效应在上升气流的作用下被拉直，形成旋转的气柱。
发展阶段：当风暴增强时，会形成超级单体雷暴（supercell thunderstorm），这种雷暴具有旋转的上升气流，是龙卷风形成的关键条件。
成熟阶段：当旋转的气柱与超级单体雷暴结合时，龙卷风开始形成并逐渐延伸至地面。
消散阶段：随着条件的变化，龙卷风逐渐变窄并逐渐升回云层，最终消散。

科里奥利力是地球自转引起的惯性力，它对大气和海洋中的大规模环流系统有重要影响。在北半球，Coriolis 力会使气流向右偏转，而在南半球则向左偏转。科里奥利力的大小与纬度和风速有关，其公式为：
$a_{COR}| = fV$
其中， $f$ 是 Coriolis 参数，定义为：
$\cdot \Omega \cdot \sin(\phi)$
其中， $\Omega$ 是地球的角速度， $\phi$ 是纬度。

气旋流和梯度风

气旋流（Cyclonic flow）是指气流呈逆时针方向旋转（北半球）或顺时针方向旋转（南半球），通常与低气压系统相关。反气旋流（Anticyclonic flow）则相反，与高气压系统相关。在气旋流中，Coriolis 力和离心力（Centrifugal force, CeF）的方向相同，而压力梯度力（Pressure Gradient Force, PGF）则方向相反。其平衡关系可以表示为：
$PGF = CF + C e F$
在反气旋流中，PGF 和 CeF 方向相同，CF 方向相反，平衡关系为：
$PGF + C e F - CF = 0$

当气流的曲率半径 $R$ 很小时，Coriolis 力可以忽略不计，此时的风速称为气旋风速（Cyclostrophic wind speed），其公式为：
$\sqrt{\frac{-R \partial \phi}{\partial n}}$
其中， $\partial \phi / \partial n$ 是压力梯度力的法向分量。

梯度风近似

在实际应用中，梯度风近似（Gradient wind approximation）用于描述气流在气旋和反气旋中的运动。其公式为：
$V^2 / R + fV = -\frac{\partial \phi}{\partial n}$
解此方程可得梯度风速 $V$ ：
$\frac{-fR \pm \sqrt{f^2 R^2 + 4f R \partial \phi / \partial n}}{2}$

基于科里奥利力的龙卷风优化器 (TOC)

TOC 算法是一种基于群体的优化算法，其核心思想是模拟龙卷风的形成和演化过程。算法中包含风暴（Windstorms）、雷暴（Thunderstorms）和龙卷风（Tornadoes）三种角色，它们在搜索空间中移动并相互作用，最终找到全局最优解。
TOC算法的实现步骤如下：
初始化种群：随机生成风暴、雷暴和龙卷风的位置。
适应度评估：计算每个个体的适应度值。
风速更新：根据科里奥利力和气流速度公式更新风速。
位置更新：根据风速更新个体的位置。
随机形成风暴：在搜索空间中随机生成新的风暴。
迭代终止：重复上述步骤，直到达到最大迭代次数或满足终止条件

初始化种群

TOC 算法的初始化过程是随机生成一组设计变量（风暴、雷暴和龙卷风）的位置。这些位置在搜索空间的上下界之间均匀分布。初始化公式为：
$y_{i,j} = l_j + \text{rand} \times (u_j - l_j)$
其中， $y_{i,j}$ 是第 $i$ 个个体在第 $j$ 维的初始位置， $l_j$ 和 $u_j$ 分别是第 $j$ 维的下界和上界， $\text{rand}$ 是 $[0, 1]$ 之间的随机数。

适应度评估

每个风暴、雷暴和龙卷风的适应度值通过目标函数计算得出。适应度值越小，表示该个体越接近全局最优解。适应度计算公式为：
$\text{fit}_i = \text{fit}(y_{i,1}, y_{i,2}, \dots, y_{i,d})$
其中， $\text{fit}_i$ 是第 $i$ 个个体的适应度值， $d$ 是问题的维度。

风暴的演化

风暴会根据其强度和演化速度向雷暴和龙卷风移动。风暴的演化过程通过适应度值来衡量，适应度值越低的风暴越有可能演化为雷暴或龙卷风。演化过程的数学模型如下：

风暴分配到雷暴和龙卷风

根据雷暴和龙卷风的适应度值，计算每个雷暴或龙卷风分配到的风暴数量。假设 $n t$ 为雷暴数量， $n o$ 为龙卷风数量， $n w$ 为风暴数量，则有：
$n t o = n t + n o$
每个雷暴或龙卷风分配到的风暴数量为：
$n_{\dot{w}k} = \left\lfloor \frac{f_k}{\sum_{k=1}^{nto} f_k} \times nw \right\rfloor$
其中， $f_k$ 是第 $k$ 个雷暴或龙卷风的适应度值， $\lfloor \cdot \rfloor$ 表示向下取整。

风暴的速度更新

风暴的速度更新考虑了 Coriolis 力的影响。更新公式为：
$\vec{v}_{t+1}^{i} = \begin{cases} \eta (\mu \vec{v}_{t}^{i} - c \frac{f \times R_l}{2} + \sqrt{CF_l}) & \text{if } \text{rand} \geq 0.5 \\ \eta (\mu \vec{v}_{t}^{i} - c \frac{f \times R_r}{2} + \sqrt{CF_r}) & \text{if } \text{rand} < 0.5 \end{cases}$
其中：

$\vec{v}_{t+1}^{i}$ 是第 $i$ 个风暴在第 $t + 1$ 次迭代的速度；
$\vec{v}_{t}^{i}$ 是第 $i$ 个风暴在第 $t$ 次迭代的速度；
$\eta$ 是收缩因子，用于控制风暴的收敛速度，其值为：
$\eta = \frac{2}{\sqrt{2 - \chi - \frac{\chi^2}{4}}}$
其中， $\chi$ 是加速因子，取值为 4.10；
$\mu$ 是模糊自适应动能因子，取值为：
$\mu = 0.5 + \frac{\text{rand}}{2}$
$R_l$ 和 $R_r$ 分别是北半球和南半球的曲率半径，其值为：
$R_l = \frac{2}{1 + e^{-(t-T)/2}}, \quad R_r = -\frac{2}{1 + e^{-(t-T)/2}}$
其中， $t$ 是当前迭代次数， $T$ 是最大迭代次数；
$c$ 是随机因子，取值范围为：
$\times \delta_1 \times wr$
其中， $b r$ 是常数，取值为 100000； $\delta_1$ 是符号变化因子，

参考文献
[1]Braik, M., Al-Hiary, H., Alzoubi, H. et al. Tornado optimizer with Coriolis force: a novel bio-inspired meta-heuristic algorithm for solving engineering problems. Artif Intell Rev 58, 123 (2025). https://doi.org/10.1007/s10462-025-11118-9

二、无人机（UAV）三维路径规划

单个无人机三维路径规划数学模型参考如下文献：

Phung M D , Ha Q P . Safety-enhanced UAV Path Planning with Spherical Vector-based Particle Swarm Optimization[J]. arXiv e-prints, 2021.

每个无人机的目标函数由路径长度成本，安全性与可行性成本、飞行高度成本和路径平滑成本共同组成：

2.1路径长度成本

路径长度成本由相邻两个节点之间的欧氏距离和构成，其计算公式如下：
在这里插入图片描述

2.2路径安全性与可行性成本

在这里插入图片描述

路径安全性与可行性成本通过下式计算：

在这里插入图片描述

2.3路径飞行高度成本

在这里插入图片描述

飞行高度成本通过如下公式计算所得：
在这里插入图片描述

2.4路径平滑成本

在这里插入图片描述

投影向量通过如下公式计算：

在这里插入图片描述

转弯角度的计算公式为：
在这里插入图片描述

爬坡角度的计算公式为：

在这里插入图片描述

平滑成本的计算公式为：
在这里插入图片描述

2.5总成本（目标函数）

在这里插入图片描述

总成本由最优路径成本，安全性与可行性成本、飞行高度成本和路径平滑成本的线性加权所得。其中，b为加权系数。

三、实验结果

在三维无人机路径规划中，无人机的路径由起点，终点以及起始点间的点共同连接而成。因此，自变量为无人机起始点间的各点坐标，每个无人机的目标函数为总成本（公式9）。本文研究3个无人机协同路径规划，总的目标函数为3个无人机的总成本之和。

%% 第一个无人机 起始点
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;

figure
plot(Convergence_curve,'LineWidth',2)
xlabel('Iteration');
ylabel('Best Cost');
grid on;