2025最新智能优化算法：龙卷风优化算法（ Tornado Optimizer with Coriolis force ，TOC) 求解23个经典函数测试集，MATLAB

一、龙卷风优化算法

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

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

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

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

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

气旋流和梯度风

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

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

梯度风近似

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

风暴的演化

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

风暴分配到雷暴和龙卷风

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

风暴的速度更新

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

• ${\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}},R_r=-\frac{2}{1+e^{-(t-T)/2}}$$
  其中，$t$ 是当前迭代次数，$T$ 是最大迭代次数；
• $c$ 是随机因子，取值范围为：
  $$c=br\times {\delta }_1\times wr$$
  其中，$br$ 是常数，取值为 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

二、23个函数介绍

参考文献：

[1] Yao X, Liu Y, Lin G M. Evolutionary programming made faster[J]. IEEE transactions on evolutionary computation, 1999, 3(2):82-102.

三、部分代码及结果

clear;
clc;
close all;
warning off all;

SearchAgents_no=50;    %Number of search solutions
Max_iteration=500;    %Maximum number of iterations

Func_name='F1'; % Name of the test function

% Load details of the selected benchmark function
[lb,ub,dim,fobj]=Get_F(Func_name); 

tic;
[Best_score,Best_pos,cg_curve]=SGA(SearchAgents_no,Max_iteration,lb,ub,dim,fobj); 
tend=toc;

% figure('Position',[500 500 901 345])
%Draw search space
subplot(1,2,1);
func_plot(Func_name);
title('Parameter space')
xlabel('x_1');
ylabel('x_2');
zlabel([Func_name,'( x_1 , x_2 )'])

%Draw objective space
subplot(1,2,2);
semilogy(cg_curve,'Color','m',LineWidth=2.5)
title(Func_name)

% title('Objective space')
xlabel('Iteration');
ylabel('Best score obtained so far');

axis tight
grid on
box on
legend('SGA')

display(['The running time is:', num2str(tend)]);
display(['The best fitness is:', num2str(Best_score)]);
display(['The best position is: ', num2str(Best_pos)]);

四、完整MATLAB代码见下方名片