提出了一种名为动麦优化算法（AOO）的新型元启发式算法，该算法受到环境中动麦自然行为的启发。AOO模拟了动画燕麦的三种独特行为：（i）通过风、水和动物等自然元素进行种子传播；（ii）在吸湿运动的影响下，动画燕麦种子的主芒发生扭曲和旋转，使整个种子能够滚动和传播；（iii）在滚动传播过程中，遇到障碍物时储存能量，在特定条件下触发推进机制以进一步传播种子。

为了模拟动画燕麦种子的传播过程，借鉴了Rivka Elbaum、Dennis Evangelista和其他研究人员对滚动和喷射机制以及种子传播过程中涉及的力学和运动模型进行的研究，提出了一种包括探索和开发的两阶段方法来模拟动画燕麦种子的传播过程，如图3所示。这个过程涉及三种传播模式：

图1 勘探开发阶段。
1.种子传播是通过风、水或动物的影响或作用实现的。
2.在传播过程中没有障碍的情况下，种子仅通过吸湿运动诱导的滚动来经历位置变化。
3.在散布过程中，当遇到障碍物时，我们假设主芒在储存能量的驱动下经历突然弹射。整个种子散布过程可以认为是主要由弹射驱动的弹丸运动。

3.2 初始化

AOO 从一组随机解开始，如公式(1)所示。

$$X=\left[\begin{array}{cccccc}{x}_{1,1}& \dots & x_{1,j}& \dots & x_{1,Dim-1}& x_{1,Dim}\\ x_{2,1}& \dots & x_{2,j}& \dots & x_{2,Dim-1}& x_{2,Dim}\\ ⋮& & ⋮& & ⋮& \\ x_{i,j}& \dots & \dots & & \dots & x_{i,j}\\ x_{N-1,1}& \dots & x_{N-1,j}& \dots & x_{N-1,Dim-1}& x_{N-1,Dim}\\ x_{N,1}& \dots & x_{N,j}& \dots & x_{N,Dim-1}& x_{N,Dim}\end{array}\right]\text{\hspace{1em}(1)}$$

其中 $x_i$ 表示第 $i$ 个子群中的个体及其对应的位置信息。$Dim$ 表示问题的维度大小，$N$ 表示子群中的个体数量。矩阵 $X$ 中的每个 $x_{i,j}$ 可以使用公式(2)计算：

$$x_{i,j}=r\times (U{B}_j-L{B}_j)+L{B}_j,i=1,2,\dots ,N,j=1,2,\dots ,Dim\text{\hspace{1em}(2)}$$

其中 $r$ 是区间 [0, 1] 内的随机数，$U{B}_j$ 和 $L{B}_j$ 是给定问题第 $j$ 维的上界和下界。

3.3 参数计算

在我们开发的数学模型中，动画燕麦种子传播过程的运动特征与种子主芒的长度、其质量以及滚动过程中偏心率系数相关。这些可以通过以下公式计算：

$$m=0.5\times \frac{r}{dim}\text{\hspace{1em}(3)}$$

$$L=N\times \frac{r}{dim}\text{\hspace{1em}(4)}$$

$$e=0.5\times \frac{r}{dim}\text{\hspace{1em}(5)}$$

$$c=1-{\left(\frac{t}{T}\right)}^3\text{\hspace{1em}(6)}$$

其中 $m$ 表示动画燕麦种子的质量，$r$ 是 0 到 1 之间的随机数，$dim$ 是目标问题的维度，$N$ 是种群大小，$L$ 是动画燕麦种子主芒的长度，$e$ 是种子滚动过程中的偏心率系数，$t$ 是当前迭代次数，$T$ 是最大迭代次数，$c$ 是动态调整因子。

3.4 探索阶段

在部分动画燕麦种子从植物上掉落后，它们的传播主要通过风、水或动物的影响或作用实现。这种传播模式表现出显著的随机性，允许在广泛的解空间中进行探索。位置更新如下：

$$W=\frac{e}{r}\times (2\times r_{dim}-1)\otimes UB\text{\hspace{1em}(7)}$$

$$\{\begin{array}{l}{X}_{i+1}(i)=\frac{1}{N}\times {\sum }_{i=1}^N{X}_i(i)+W,\text{if }\text{mod}(i,N/10)=0,\\ X_{i+1}(i)=X_{best}+W,\text{if }\text{mod}(i,N/10)=1,\\ X_{i+1}(i)=X_i(i)+W,\text{else.}\end{array}\text{\hspace{1em}(8)}$$

其中 $UB$ 是目标问题的上界，$X_{i+1}(i)$ 是第 $i+1$ 代中第 $i$ 个个体的位置，$X_{best}$ 是种群中最佳个体的位置。

3.5 开发阶段

在这个阶段，剩余的动画燕麦种子将根据是否遇到障碍物分为两种传播方法。此外，我们假设两种情况的概率相等。

在没有障碍物的情况下，种子位置的变化通过吸湿滚动实现，这是由湿度引起的应力梯度驱动的。受 Lindner 等人 [2] 的启发，他们证明了细胞微纤维取向控制各向异性膨胀，我们使用临界曲率进行快照屈曲来模拟这种运动。滚动机制通过扭矩公式的偏心旋转数学表达：

$$A=UB-\lfloor \frac{UB\times t\times \sin (2\times \pi \times r)}{T}\rfloor \text{\hspace{1em}(9)}$$

$$R=\left(m\times e+L^2\right)\times \frac{r_{dim}(-A,A)}{dim}\text{\hspace{1em}(10)}$$

$$\text{Levy}(\text{dim})=0.01\times \frac{\mu \times \sigma }{|v{|}^{1/\beta }}\text{\hspace{1em}(11)}$$

$$\sigma ={\left(\frac{\Gamma (1+\beta )\times \sin \left(\frac{\pi \beta }{2}\right)}{\Gamma \left(\frac{1}{2}\right)\times \beta \times 2^{(\beta -1)/2}}\right)}^{1/\beta }\text{\hspace{1em}(12)}$$

$$\{\begin{array}{l}{X}_{i+1}(i)=X_{best}+R+c\times \text{Levy}(\text{dim})\otimes X_{best}\\ B=UB-\lfloor \frac{UB\times t\times \cos (2\times \pi \times r)}{T}\rfloor \end{array}\text{\hspace{1em}(13)}$$

其中 $r_{dim}(-A,A)$ 是从 $-A$ 到 $A$ 的随机矩阵，其维度与目标问题的维度匹配。在 Levy 飞行中，$\mu$ 表示期望位置或平均值，通常用于调整步长大小，$r$ 是 0 到 1 之间的随机值。这种随机性有助于控制飞行过程中运动的方向和距离。$\sigma$ 是步长分布的尺度参数，控制步长分布的宽度并确定步长变化的范围。$v$ 表示当前速度向量，反映粒子或个体的运动状态。$\beta$ 是稳定分布的参数，决定步长分布的形状并影响步长变化的随机性和多样性。最后，$\Gamma$ 表示伽马函数，是阶乘函数的连续扩展，用于计算非整数参数的阶乘值。在使用 Levy 飞行时，我们只需要预定义问题的维度，本文中 $\beta$ 设置为 1.5。

当种子在传播过程中遇到障碍物时，我们假设主芒在储存能量后经历由喷射驱动的喷射运动。整个种子传播过程以投射运动为特征。为了模拟这一过程，我们使用简化的投射运动模型进行位置更新，位置更新公式如下：

$$k=0.5\times r\text{\hspace{1em}(14)}$$

$$x=3\times \frac{r}{dim}\text{\hspace{1em}(15)}$$

$$\theta =\pi \times r\text{\hspace{1em}(16)}$$

$$\alpha =\frac{2}{\pi }\times e^{\frac{r}{T}}\text{\hspace{1em}(17)}$$

$$J=\frac{2\times k\times x^2\times \sin (2\theta )}{mg}\times \frac{r_{dim}(-B,B)}{dim}\times (1-\alpha )\text{\hspace{1em}(18)}$$

$$X_i(i)=X_{best}+J+c\times \text{Levy}(\text{dim})\otimes X_{best}\text{\hspace{1em}(19)}$$

其中 $k$ 是种子主芒的弹性系数，$x$ 是主芒在能量储存期间用于喷射的长度，$\theta$ 是投射轨迹与地面之间的角度，$\alpha$ 是投射运动期间的空气阻力系数，$r$ 是 0 到 1 之间的随机数，$r^{\prime }$ 是 0 到 1 之间的随机数。

3.6 AOO 分析

图1 AOO算法流程图

4.1 二元感知模型

在当前无线传感器网络（WSN）覆盖场景中，最常见且理想的感知模型是二元感知模型[27]。如图2所示，$Q_1$和$Q_2$代表平面区域内的两个随机点。$Q_1$位于传感器的检测区域内，其感知概率为1。$Q_2$的感知概率为0，因为它位于传感器的检测区域之外。二元感知模型的公式可以表示为：

$$P(M_i,Q)=\{\begin{array}{ll}1,& d(M_i,Q)\le r\\ 0,& d(M_i,Q)>r\end{array}\text{\hspace{1em}(1)}$$

其中，$M$表示监控区域内的传感器节点，所有节点具有相同的感知半径$r$，节点集表示为$(M_1,M_2,M_3\cdots M_n)$。$d(M,Q)$是监控点$Q$与传感器节点$M$之间的欧几里得距离。$P(M_i,Q)$是传感器节点$M$检测点$Q$的概率。

4.1.2 传感器节点的二维和三维覆盖模型

无线传感器网络的监控区域被划分为二维矩形区域$E\times F$和三维监控区域$E\times F\times L$，其中在监控区域内部署了$n$个传感器节点。

在二维区域中，传感器节点$M_i$的位置坐标为$(X_i,Y_i)$，监控节点$Q$的位置坐标为$(X_j,Y_j)$。传感器节点$M_i$和目标监控点$Q$之间的欧几里得距离定义为：

$$d(M_i,Q)=\sqrt{(X_i-X_j{)}^2+(Y_i-Y_j{)}^2}\text{\hspace{1em}(2)}$$

在三维区域中，传感器节点$M_i$的位置坐标为$(X_i,Y_i,Z_i)$，监控节点$Q$的位置坐标为$(X_j,Y_j,Z_j)$。传感器节点$M_i$和目标监控点$Q$之间的欧几里得距离定义为：

$$d(M_i,Q)=\sqrt{(X_i-X_j{)}^2+(Y_i-Y_j{)}^2+(Z_i-Z_j{)}^2}\text{\hspace{1em}(3)}$$

此外，多个节点可以覆盖监控区域内的同一目标点。节点集$M$检测监控点$Q$的概率由下式给出：

$$P(M,Q)=1-\prod _{i=1}^n(1-P(M_i,Q))\text{\hspace{1em}(4)}$$

覆盖率可以定义为监控区域总面积中所有传感器节点覆盖的比例[28]。因此，在二维WSN环境中的覆盖率表示为：

$${\text{Cov}}_{2D}=\frac{{\sum }_{Q\in E\times F}P(M,Q)}{E\times F}\text{\hspace{1em}(5)}$$

在三维WSN环境中的覆盖率表示为：

$${\text{Cov}}_{3D}=\frac{{\sum }_{Q\in E\times F\times L}P(M,Q)}{E\times F\times L}\text{\hspace{1em}(6)}$$

完整代码

%%Animated Oat Optimization Algorithm
function [Best_Score,Best_X,convergence_curve]= AOO(pop_size,iter_max,lb,ub,dim,fobj)
    %% initialize
    % search_history,ave_fit,x_1st
    Best_Score=inf; 
    convergence_curve=zeros(iter_max,1);
    search_history = zeros(iter_max*pop_size, dim); % search history
    ave_fit = zeros(iter_max, 1); % average fitness
    x_1st= zeros(iter_max, dim); % 1st trajectory
    
    x = 3 * rand(1, pop_size) / dim;
    m = 0.5 * rand(1, pop_size) / dim;
    L = pop_size / dim * rand(1, pop_size);
    e = m;
    g = 9.8/dim;
    X=initialization(pop_size,dim,ub,lb);
    
    
    for i = 1:pop_size
        Pop_Fit(i) = fobj(X(i,:));
        if Pop_Fit(i)<Best_Score
            Best_Score=Pop_Fit(i);
            Best_X=X(i,:);
        end
    end
    convergence_curve(1) = Best_Score;
    
    %% iteration part
    
    for t = 1 :1:iter_max
        theta = pi * rand(1,pop_size);
        c = (1 -  t / iter_max)^3;
        P = levy(pop_size,dim,1.5);
    
        for i = 1 :pop_size
            if rand > 0.5
                W = c/pi *(2 * rand(1,dim) - 1) .* ub  ;
                if mod(i,pop_size/10) == 0
                    X(i,:) = mean(X) + W;
                elseif mod(i,pop_size/10) == 1
                    X(i,:) = Best_X + W ;
                else
                    X(i,:) = X(i,:) + W ;
                end
            else
                if rand > 0.5
                    A = ub - abs(ub * t * sin(2 * pi * rand) / iter_max);
                    R = (m(i) * e(i) + L(i) ^2) /dim * unifrnd( -A , A, 1, dim);
                    X(i,:) = Best_X + R  + c * P(i,:) .* Best_X  ;
                else
                    k = 0.5 + 0.5 * rand;
                    B = ub - abs(ub * t * cos(2 * pi * rand) / iter_max);
                    alpha = 1 / pi * exp((randi([0,t]) / iter_max));
                    J = 2 * k * x(i)^2 * sin (2 * theta(i)) / m(i) / g * (1 - alpha) /dim  * unifrnd( -B , B, 1, dim);
                    X(i,:) = Best_X + J + c * P(i,:) .* Best_X ;
                end
            end
        end
    
        X = boundaryCheck(X, lb, ub);
        Pop_Fit = fobj(X);
        for i=1:pop_size
            Pop_Fit(i) = fobj(X(i,:));
            if Pop_Fit(i)<Best_Score
                Best_Score=Pop_Fit(i);
                Best_X=X(i,:);
            end
        end
        search_history ((t-1)*pop_size+1:t*pop_size,:)=X; % search history
        ave_fit(t) = mean(Pop_Fit); % average fitness
        x_1st(t,:)= X(1,:); % 1st trajectory
        convergence_curve(t)=Best_Score;
    end
end
function [ X ] = boundaryCheck(X, lb, ub)
    for i=1:size(X,1)
        FU=X(i,:)>ub;
        FL=X(i,:)<lb;
        X(i,:)=(X(i,:).*(~(FU+FL)))+ub.*FU+lb.*FL;
    end
end
function X=initialization(SearchAgents_no,dim,ub,lb)
    Boundary_no= size(ub,2); % numnber of boundaries
    if Boundary_no==1
        X=rand(SearchAgents_no,dim).*(ub-lb)+lb;
    end
    % If each variable has a different lb and ub
    if Boundary_no>1
        for i=1:dim
            ub_i=ub(i);
            lb_i=lb(i);
            X(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
        end
    end
end

function [z] = levy(n,m,beta)
    num = gamma(1+beta)*sin(pi*beta/2); % used for Numerator 
    den = gamma((1+beta)/2)*beta*2^((beta-1)/2); % used for Denominator
    sigma_u = (num/den)^(1/beta);% Standard deviation
    u = random('Normal',0,sigma_u,n,m); 
    v = random('Normal',0,1,n,m);
    z =u./(abs(v).^(1/beta));
  end

Wang R B, Hu R B, Geng F D, et al. The Animated Oat Optimization Algorithm: A Nature-Inspired Metaheuristic for Engineering Optimization and a Case Study on Wireless Sensor Networks[J]. Knowledge-Based Systems, 2025: 113589.