2025年中科院三区全新算法，恒星振荡优化器：受自然启发的元启发式优化，完整MATLAB代码免费获取

文章提出了一种受自然启发的元启发式优化算法，称为恒星振荡优化器（SOO），SOO的灵感来自星震学领域，该领域检查恒星的振荡行为，以了解其内部结构、物理特性和进化阶段。SOO模拟在恒星脉动中观察到的动态膨胀和收缩阶段，以平衡勘探和开发。

3 恒星振荡优化器（SOO）

本节首先描述了所提出的自然启发式元启发式优化算法（SOO）背后的灵感，然后是相应的数学模型和算法描述。

3.1 SOO灵感

恒星振荡优化器（SOO）受到恒星中观测到的振荡现象的启发，这是一个被称为天体物理学的领域，如图1所示，

图1 恒星振荡现象的插图
特别是恒星脉动中的正弦和余弦模式如图2所示。

图2 恒星振荡现象中的正弦和余弦模式图解

天体物理学专注于研究恒星的振荡或脉动，以了解其内部结构、物理属性和演化状态[34]。此外，在天体物理学中，恒星的角频率描述了恒星完成一个振荡周期的速度，而周期表示一个完整的膨胀和收缩周期所需的时间[62]。这些概念对于模拟恒星的脉动至关重要。在SOO中，角频率和周期被纳入控制搜索振荡器的振荡运动。角频率，与周期相关，控制振荡器在解空间中扩展和收缩的速度。周期随时间演变，允许搜索振荡器从早期迭代的广泛探索逐渐过渡到搜索进展的集中探索。这种动态调整最小化了恒星振荡的时间变化，为SOO提供了在复杂优化景观中导航的灵活性。

在SOO中，正弦函数表示恒星振荡的外向扩展阶段，这象征着优化过程中的探索阶段。此阶段使SOO能够探索解空间的广泛范围，类似于恒星在其脉动期间的扩展。受正弦运动控制的振荡器旨在探索解空间的新区域，确保覆盖不同的区域[63]。

另一方面，余弦函数表示恒星振荡的收缩阶段，对应于优化的利用阶段。在此阶段，SOO专注于通过缩小搜索范围来细化探索期间识别的解决方案，围绕最有希望的区域。余弦运动模拟的收缩效应使振荡器朝向全局最优移动，类似于恒星在膨胀后的收缩[64]。这种正弦基础探索和余弦基础利用的平衡方法确保了SOO能够有效地导航优化问题的全局和局部景观。

天体物理学中的一个关键概念是光曲线，它表示恒星亮度随时间的变化，由于其振荡。在SOO中，光曲线的概念反映在优化过程中对每个振荡器适应度值（或亮度）的连续跟踪中。振荡器的适应度值类似于恒星的亮度，随着恒星的膨胀和收缩而变化[65]。

SOO中的适应度评估对于确定表现最佳的振荡器至关重要，这些振荡器指导后续迭代中其他振荡器的运动。前3个表现最佳的振荡器通过对其适应度值进行排序来识别，并用于影响种群的其余部分，确保搜索保持集中在最有希望的区域[66]。这种自适应机制，受到恒星光曲线连续监测的启发，确保了SOO朝着最优解收敛。

SOO最创新的方面之一是它结合了正弦和余弦运动来生成新解决方案。通过平均正弦和余弦函数的运动向量，SOO确保振荡器在探索和开发之间保持平衡。这种平均机制类似于恒星脉动中发现的平衡，其中膨胀和收缩阶段之间的平衡是其稳定性的关键[67]。

3.2 SOO数学模型

SOO的数学模型由一系列方程描述，这些方程模拟了解空间中解的振荡运动，类似于恒星的脉动。

恒星振荡器的位置在$[b,ub]$的边界内初始化，在$D$维空间中。设振荡器的数量表示为$N_{\text{osc}}$。振荡器的位置由矩阵$\mathbf{X}\in {\mathbb{R}}^{N_{\text{osc}}\times D}$表示，其中每一行对应于一个振荡器的位置。每个振荡器的亮度通过适应度函数$f(\mathbf{x})$评估，其中$\mathbf{x}\in {\mathbb{R}}^D$ 。

恒星位置的初始化表示为公式(1)：

$${\mathbf{X}}_{\text{init}}=\text{initialization}(N_{\text{osc}},D,ub,lb)\text{\hspace{1em}(1)}$$

其中，${\mathbf{X}}_{\text{init}}$表示恒星振荡器的初始位置。$N_{\text{osc}}$是振荡器的数量，$D$是搜索空间的维度，向量$ub$和$lb$分别表示搜索空间的上界和下界[69]。

初始亮度和相应的最佳相位位置初始化如下公式(2)所示：

$$L_{\text{best}}=\infty ,{\mathbf{x}}_{\text{best}}=0_D\text{\hspace{1em}(2)}$$

其中，$L_{\text{best}}$表示迄今为止找到的最佳亮度，最初设置为无穷大。${\mathbf{x}}_{\text{best}}$表示振荡器在相空间中的当前最佳位置，最初设置为维度$D$的零向量，对于每次迭代，每个振荡器的亮度根据其位置[70]进行更新。

振荡器$P(t)$的周期在每次迭代$t$中更新，如公式(3)所示：

$$P(t)=P_0+\Delta P\cdot t\text{\hspace{1em}(3)}$$

其中，$P(t)$表示迭代$t$的振荡周期，$P_0$是初始周期，$\Delta P$是每次迭代的周期增量。该公式考虑了振荡周期的逐渐调整[71]。

振荡器$\omega (t)$的角频率如公式(4)所示：

$$\omega (t)=\frac{2\pi }{P(t)}\text{\hspace{1em}(4)}$$

其中，$\omega (t)$表示迭代$t$的振荡器的角频率，计算为周期$P(t)$的倒数。因子$2\pi$将周期转换为角频率，这是SOO振荡器行为中的一个参数[72]。

缩放因子$S(t)$如公式(5)所示更新：

$$S(t)=2-\frac{2\cdot t}{T}\text{\hspace{1em}(5)}$$

其中，$S(t)$表示迭代$t$的缩放因子，调整振荡器随时间的强度。$T$是总迭代次数。该因子随着迭代次数的增加而减小，允许SOO探索搜索空间[73]。

对于每个振荡器$i$，其位置根据天体物理振荡属性[74]进行更新。生成两个候选位置${\mathbf{x}}_{\text{osc1}}$和${\mathbf{x}}_{\text{osc2}}$，如公式(6)和(7)所示：

$${\mathbf{x}}_{\text{osc1}}(j)={\mathbf{x}}_{\text{best}}(j)-r_1{r}_3\left(\omega (t)S(t)(r_1-S(t))({\mathbf{x}}_i(j)-{\mathbf{x}}_{\text{best}}(j))\right)-|r_1\sin (r_2)\cdot |r_3\cdot {\mathbf{x}}_{\text{best}}(j)||\text{\hspace{1em}(6)}$$

其中，${\mathbf{x}}_{\text{osc1}}(j)$是振荡器$i$的第一个候选位置的第$j$个分量。术语$r_1,r_2,r_3$是均匀采样的随机值，来自区间$[0,1]$，并引入了振荡运动的随机性[75]。

第二个候选位置生成如公式(7)所示：

$${\mathbf{x}}_{\text{osc2}}(j)={\mathbf{x}}_{\text{best}}(j)-r_1{r}_3\left(\omega (t)S(t)(r_1-S(t))({\mathbf{x}}_i(j)-{\mathbf{x}}_{\text{best}}(j))\right)-|r_1\sin (r_2)\cdot |r_3\cdot {\mathbf{x}}_{\text{best}}(j)||\text{\hspace{1em}(7)}$$

其中，${\mathbf{x}}_{\text{osc2}}(j)$表示振荡器$i$的第二个候选位置。该方程类似于公式(6)，但余弦项取代了正弦项。这为振荡器提供了替代轨迹，引入了更多可能的运动[76]。

最终更新位置计算为两个候选位置的平均值，如公式(8)所示：

$${\mathbf{x}}_{\text{new}}(j)=r_3\cdot \frac{{\mathbf{x}}_{\text{osc1}}(j)+{\mathbf{x}}_{\text{osc2}}(j)}{2}\text{\hspace{1em}(8)}$$

其中，${\mathbf{x}}_{\text{new}}(j)$表示振荡器$i$的第$j$维的新位置。最终位置是通过平均两个候选位置获得的，加权由随机因子$r_3$，这确保了搜索过程的进一步随机性[77]。

在优化过程中，振荡器的位置受到前3个最佳恒星的平均位置的影响。这个平均位置由公式(9)给出：

$${\mathbf{x}}_{\text{avg}}=\frac{1}{3}\sum _{k=1}^3{\mathbf{x}}_{\text{top},k}\text{\hspace{1em}(9)}$$

其中，${\mathbf{x}}_{\text{avg}}$表示前3个最佳恒星的平均位置。${\mathbf{x}}_{\text{top},k}$表示第$k$个最佳振荡器的位置，求和平均化这些位置[78]。这个平均位置用作生成新位置的参考点[78]。

基于振荡运动的新位置由公式(10)给出：

$${\mathbf{x}}_{\text{osc}}={\mathbf{x}}_{\text{avg}}+0.5\left(\sin (r_7)({\mathbf{x}}_{r_1}-{\mathbf{x}}_{r_3})+\cos ((1-r_7)\pi )({\mathbf{x}}_{r_1}-{\mathbf{x}}_{r_3})\right)\text{\hspace{1em}(10)}$$

其中，${\mathbf{x}}_{\text{osc}}$表示基于振荡运动生成的新位置。随机值$r$控制正弦和余弦函数的影响，这些函数控制更新的方向和幅度。术语${\mathbf{x}}_{r_1},{\mathbf{x}}_{r_2},{\mathbf{x}}_{r_3}$是随机选择的振荡器的位置。此方程确保新位置受到其他振荡器多样性的影响[79]。

如果$r_j\ge 0.5$，则第$j$个位置以0.5的概率更新，如公式(11)所示：

$$x_{\text{new}}(j)=\{\begin{array}{ll}{\mathbf{x}}_{\text{osc}}(j),& \text{if }{r}_j\ge 0.5\\ {\mathbf{x}}_{\text{old}}(j),& \text{otherwise}\end{array}\text{\hspace{1em}(11)}$$

其中，$x_{\text{new}}(j)$表示第$j$维的更新位置。以0.5的概率，位置更新为振荡器位置${\mathbf{x}}_{\text{osc}}(j)$或保持为其先前值${\mathbf{x}}_{\text{old}}(j)$。这种概率更新引入了搜索过程中的额外随机性[80]。

在更新位置处评估适应度（亮度）。如果新适应度优于当前最佳，则位置更新如公式(12)所示：

$$L(x_{\text{new}})<L_{\text{best}}⟹L_{\text{best}}=L(x_{\text{new}}),{\mathbf{x}}_{\text{best}}={\mathbf{x}}_{\text{new}}\text{\hspace{1em}(12)}$$

其中，$L(x_{\text{new}})$表示新位置的适应度（亮度）。如果这个新适应度优于记录的最佳适应度$L_{\text{best}}$，则$L_{\text{best}}$更新为新位置[81]。

3.2 SOO算法描述

SOO算法（见算法3.2）基于天体物理学的原理。在天体物理学中，通过观察恒星的振荡行为，恒星的内部结构可以揭示关于恒星内部属性的重要见解。SOO通过将优化问题中的潜在解决方案视为恒星振荡器，每个振荡器在多维搜索空间中经历时间变化的周期运动。在每次迭代中，振荡器根据演变的振荡周期、角频率和缩放因子调整其位置（如公式(3)、(4)和(5)所示）。这些振荡运动受到随机因素的影响，因为SOO基于正弦和余弦函数生成候选位置（公式(6)和(7)），模拟了恒星振荡的自然性质。SOO还根据适应度值的概率更新振荡器的位置，反映了天体物理过程中固有的随机性。最佳解决方案，类似于恒星中最亮的振荡模式，根据适应度值不断更新（见公式(12)），确保SOO随着时间的推移朝着最优解移动。前3个最佳振荡器，代表最佳候选者，影响搜索空间中其他振荡器的运动。SOO的流程图如图3所示。

图3. SOO的资源采集阶段。

算法通过描述地位驱动的社会动态来复制，资源丰富的个体自然吸引更多机会，系统地引导种群朝着更好的解决方案发展，从而显著提高开发和增强整体优化性能。

算法1
确保高效的资源使用和详细的文档记录用于分析/应用。

%______________________________________________________________________________________________%
%  Stellar Oscillation Optimizer (SOO) source codes (version 1.1)                              %
%                                                                                              %
%  Developed in MATLAB R2018a (7.13)                                                           %
%  Author and programmer: Ali Rodan                                                            %
%                         e-Mail: alirodan@gmail.com                                           %
%                         Homepages:                                                           %
%                         1- https://scholar.google.co.uk/citations?user=n8Z3RMwAAAAJ&hl=en    %
%                         2- https://www.researchgate.net/profile/Ali-Rodan                    %
%                                                                                              %
%   Paper Title:Stellar Oscillation Optimizer: A Nature-inspired metaheuristic optimization algorithm.   %
%                   A. Rodan, A. Al-Tamimi, L. Alnemer, S. Mirjalili (2024)                                  %
%          Stellar Oscillation Optimizer: A Nature-inspired metaheuristic optimization algorithm.       %
%                                                                                 %
%______________________________________________________________________________________________%
% SOO Optimizer
function [best_luminosity, best_phase_position, luminosity_curve] = SOO(StarOsc_no, m_iter, lb, ub, dim, fobj)  
    % SOO initialization 
    best_phase_position = zeros(1, dim);  % Best position in the star's oscillation phase
    best_luminosity = inf;  % Best fitness value (luminosity)
    star_positions = initialization(StarOsc_no, dim, ub, lb);  % Initialize the positions of the star oscillators
    initial_period = 3.0;  % Initial period of oscillation
    lightcurve_iter = 1;  
    luminosity_curve = zeros(1, m_iter);  % Luminosity tracking for each iteration
    top_star_positions = zeros(3, dim);  % Top 3 best star positions
    top_luminosities = inf(1, 3);  % Top 3 best fitness values (luminosities)
    current_luminosities = zeros(1, StarOsc_no);  % Fitness values of the star oscillators
    updated_star_positions = star_positions;  % Initialize updated star positions
    % Initial evaluation of fitness (luminosity) for each star
    for i = 1:StarOsc_no
        current_luminosities(i) = fobj(star_positions(i, :));
        combined_luminosities = [top_luminosities, current_luminosities(i)];
        combined_star_positions = [top_star_positions; star_positions(i, :)];
        [sorted_luminosities, idx] = sort(combined_luminosities);
        sorted_star_positions = combined_star_positions(idx, :);
        top_luminosities = sorted_luminosities(1:3);
        top_star_positions = sorted_star_positions(1:3, :);
    end
    % Main loop 
    while lightcurve_iter <= m_iter 
    
        for i = 1:size(updated_star_positions, 1)
            current_period = initial_period + 0.001 * lightcurve_iter;  % Update period over time
            current_angular_freq = 2 * pi / current_period;  % Update angular frequency over time
            % Bound checks for star position 
            F_ub = updated_star_positions(i, :) > ub;
            F_lb = updated_star_positions(i, :) < lb;
            updated_star_positions(i, :) = (updated_star_positions(i, :) .* (~(F_ub + F_lb))) + ub .* F_ub + lb .* F_lb;
            current_luminosity = fobj(updated_star_positions(i, :));  % Evaluate new luminosity based on updated position
            if current_luminosity < best_luminosity
                best_luminosity = current_luminosity;
                best_phase_position = updated_star_positions(i, :);  % Update best position in the oscillation phase
            end
        end
        % Update the scaling factor for the oscillation
        scaling_factor = 2 - lightcurve_iter * ((2) / m_iter);
        % Update the location of the star oscillators 
        for i = 1:size(updated_star_positions, 1)
            for j = 1:size(updated_star_positions, 2)
                rand_factor1 = rand();
                rand_factor2 = rand();
                rand_factor3 = rand();
               osc_position1 = best_phase_position(j) - (rand_factor1*rand_factor3)* ((current_angular_freq * scaling_factor * rand_factor1 - scaling_factor) * ...
                    (updated_star_positions(i, j) - abs(rand_factor1 * sin(rand_factor2) * abs(rand_factor3  * best_phase_position(j)))));
                osc_position2 = best_phase_position(j) -  (rand_factor2*rand_factor3)* ((current_angular_freq * scaling_factor * rand_factor1 - scaling_factor) * ...
                    (updated_star_positions(i, j) - abs(rand_factor1 * cos(rand_factor2) * abs(rand_factor3  * best_phase_position(j)))));
                updated_star_positions(i, j) =rand_factor3*(osc_position1 + osc_position2/ 2);  % Update position based on averaged oscillation  % Update position based on averaged oscillation
            end
        end
        % Perform the update for oscillatory movement
        for i = 1:StarOsc_no
            avg_best_star_position = mean(top_star_positions, 1);  % Average of top star positions
            random_indices = randperm(StarOsc_no, 3);
            while any(random_indices == i)
                random_indices = randperm(StarOsc_no, 3);
            end
            % Generate new position based on oscillatory movement influenced by sine and cosine
            rand_factor = rand();
            oscillation_position = avg_best_star_position + 0.5 * (sin(rand_factor * pi) * (star_positions(random_indices(1), :) - star_positions(random_indices(2), :)) + ...
                cos((1 - rand_factor) * pi) * (star_positions(random_indices(1), :) - star_positions(random_indices(3), :)));
    
            star_update_position = star_positions(i, :);  % Initialize star update position
    
            % update the star's position based on oscillation_position
            for j = 1:dim
                if rand() <= 0.5
                    star_update_position(j) = oscillation_position(j); 
                end
            end
    
            % Ensure the updated position remains within bounds
            for j = 1:dim
                if isscalar(ub)
                    star_update_position(j) = min(max(star_update_position(j), lb), ub);
                else
                    star_update_position(j) = min(max(star_update_position(j), lb(j)), ub(j));
                end
            end
            % Update position if the new luminosity improves
            new_luminosity = fobj(star_update_position);
            if new_luminosity < current_luminosities(i)
                star_positions(i, :) = star_update_position;
                current_luminosities(i) = new_luminosity;
                % Update top star positions and luminosity values
                combined_luminosities = [top_luminosities, new_luminosity];
                combined_star_positions = [top_star_positions; star_update_position];
                [sorted_luminosities, idx] = sort(combined_luminosities);
                sorted_star_positions = combined_star_positions(idx, :);
                top_luminosities = sorted_luminosities(1:3);
                top_star_positions = sorted_star_positions(1:3, :);
            end
        end
        % Determine the best solution 
        best_primary_luminosity = best_luminosity;
        best_secondary_luminosity = min(top_luminosities);
        if best_primary_luminosity < best_secondary_luminosity
            best_luminosity = best_primary_luminosity;
        else
            best_luminosity = best_secondary_luminosity;
            best_phase_position = top_star_positions(1, :);
        end
        luminosity_curve(lightcurve_iter) = best_luminosity;  % Update luminosity curve based on best solution
        lightcurve_iter = lightcurve_iter + 1;
    end
end

Rodan, A., Al-Tamimi, AK., Al-Alnemer, L. et al. Stellar oscillation optimizer: a nature-inspired metaheuristic optimization algorithm. Cluster Comput 28, 362 (2025). https://doi.org/10.1007/s10586-024-04976-5