2025年9月青蒿素黏菌算法，附完整MATLAB代码

1、简介

这项研究提出了一种新型的仿生元启发式，即青蒿素黏液霉菌算法（ASMA），它将黏液霉菌的自适应觅食行为与青蒿素的趋化特性相结合，青蒿素是一种以其抗疟药特性而闻名的化合物。ASMA通过动态平衡探索和开发的混合机制增强了全局优化性能。代理通过类趋化行为遍历搜索空间，而青蒿素启发的组件引导它们走向有希望的解决方案，从而加速收敛并提高解决方案的准确性。该算法在CEC 2017基准集上进行了评估，并展示了优于几种领先算法的性能。此外，ASMA用于动态频率约束下桁架结构的尺寸优化，这需要大量的特征值计算。为了应对这一挑战，采用稀疏矩阵的最大特征值（LESM）技术仅提取与动态约束相关的主要特征值，避免了全谱特征值分析。这大大降低了计算成本，同时保持了准确性，特别是对于具有固有稀疏刚度矩阵的大规模桁架系统。ASMA-LESM框架在一系列圆顶结构上得到了验证，展示了动态约束下结构优化的高效率、可扩展性和鲁棒性。

3. 优化算法

3.1. 蛞蝓模具算法的概念

蛞蝓模具算法(SMA)是一种基于Physarum polycephalum的觅食行为的生物学启发式度量，Physarum polycephalum是一种以动态重组其网络结构以响应食物刺激而闻名的蛞蝓物种。SMA通过正负反馈机制模拟这种自适应行为，通过调节搜索代理的位置来实现。

3.1.1. 觅食和包裹食物的方法

蛞蝓模具(SM)模拟了Physarum polycephalum中静脉状结构的收缩行为。这种收缩由生物振荡器调节，当气味浓度（即适应度）高时产生更强的信号。这导致细胞质流动增加和静脉增厚。这种自适应行为提高了生物体探索和利用搜索空间的能力。代理位置的更新结合了正负反馈机制，并在数学上如下公式化：

其中 是蛞蝓模具的新位置， 是高阶浓度位置， 和 是来自SM的随机参数。Z = 0.03， 是参数，从1到0，t是当前迭代，W是SM的权重因子。Both ub和lb分别表示搜索空间的上下界。

概率参数 的产生如下：

其中 ， 是第i个解的适应度， 是迄今为止获得的全局最佳适应度。

参数 ，控制搜索强度，定义为：

权重因子 根据解决方案的排名和气味指数更新为：

其中Smell Index (t) = sort(S)，基于其适应度值对种群进行排名，r [0, 1]是随机值，  = 最大迭代次数，bF和wF分别表示种群中的最佳和最差适应度值。

3.1.2. 振荡

为了确定最佳食物浓度，SMA算法使用生物振荡器机制。该振荡器通过参数 ， 和 调节静脉宽度的变化，反映观察到的自适应增厚行为。振荡器提供了动态反馈循环，在搜索空间的有前景区域中放大信号。这个过程模拟了一种正强化机制，增强了蛞蝓寻找食物的能力，无论环境条件如何，都能增强光线稀缺或周围环境[39]。

3.2. 青蒿素

青蒿素的灵感来自青蒿素的药理行为，这是一种著名的抗疟疾动物化合物[40]。青蒿素在体内的生物机制通过模拟寄生虫在体内扩散并消灭寄生虫，从而在AO算法中探索和细化复杂优化问题的解决方案。

3.2.1. 广泛的药物分布和消除

在疟疾治疗的初始阶段，高剂量的青蒿素被施用以快速缓解症状并减少寄生虫种群[41]。这个阶段模拟了青蒿素在体内的控制扩散，允许化合物高效地扩散到广泛的区域。类似地，在AO算法中，搜索代理最初分散在解决方案空间中，以实现全局探索并获得对景观的整体理解。随着搜索进程的推进，代理的运动越来越多地受到青蒿素启发的扩散机制的影响。代理开始专注于搜索空间中更精细的区域，减少探索范围并加强开发。这种行为反映了青蒿素随时间递减效果的逐渐过渡，反映了从探索到开发的逐渐过渡。随着函数评估次数的增加，代理越来越多地优先考虑局部搜索策略，将其努力引向最有希望的区域。

探索和开发之间的平衡由控制参数K控制，定义为[42]：

其中 是探索衰减率（通常  = 6），t是当前迭代， 是最大评估次数。

在开发阶段，代理的位置根据以下公式更新：

其中 是代理i的第j维， 是全局最佳位置。药物浓度衰减由指数衰减函数建模：

其中 ，药物衰减率设置为 ，F的指数衰减调节了代理运动强度。随着算法的进展，搜索的影响逐渐减少探索并增加开发，使搜索集中在更有希望的区域。

3.2.2. 选择性药物靶向

在青蒿素治疗的第二阶段，药物选择性地靶向残留寄生虫。类似地，AO通过采用适应度归一化突变机制，加强了围绕有希望区域的搜索。归一化适应度在时间i计算为：

这种机制确保了具有较高相对适应度的代理对搜索方向有更大的贡献，从而增强了对有希望解决方案的收敛。

3.2.3. 治疗后愈合阶段

为了避免过早收敛并增强搜索的鲁棒性，引入了最终愈合阶段。该机制允许代理利用以下更新规则从局部最优中概率性地逃脱：

这种机制模拟了治疗后阶段，其中休眠的寄生虫被消灭。它允许代理从局部陷阱中逃脱，从而随着时间的推移提高收敛的可靠性和整体解决方案质量。

3.3. 青蒿素蛞蝓模具算法(ASMA)

ASMA是一种生物学启发式度量，融合了青蒿素的药理学和SM的自适应觅食行为。该算法由药物扩散和浓度衰减的原理引导，动态控制其搜索策略。如图2所示，ASMA引入了关键动态参数：探索因子K、药物浓度E和突变概率M。最初，K和F都设置为高值以鼓励广泛的探索。随着搜索的进行，这些参数逐渐衰减，将算法转向局部开发。相比之下，突变率M逐渐降低，有助于保持多样性并避免过早收敛。其他参数的同步衰减促进了从探索到开发的平稳过渡。ASMA的另一个组成部分引入了自适应加权机制，其中代理根据候选解决方案的适应度质量获得更大的权重。高适应度区域被赋予更大的权重以加强局部搜索，而低适应度区域则被赋予较小的权重以促进探索并从局部最优中逃脱。这种加权机制通过引导代理朝向有希望的区域，同时保持种群多样性，提高了搜索效率。

3.3.1. 使用改进的 范围的搜索过程和更新规则

靠近已知最佳解的代理通过调整其位置来参与局部开发，以增强在有希望区域内的搜索。更新规则定义如下：

这种公式允许每个代理通过围绕其自身位置或全局最佳位置进行扰动来利用已知的高质量区域。

为了促进全局探索，代理根据邻近代理的加权影响调整其位置（见式(6)）。这种调整通过随机排列策略实现。增强的全局探索建模如下：

在原始SMA中，变量边界 在式(5)中定义。在ASMA中，这已改进为 ，其中E是一个指数衰减参数（见式(10)）。这种动态范围提供了自适应步长控制。在早期迭代中，E的大值允许广泛的探索，而在后期迭代中，其衰减导致更细粒度的搜索行为。搜索范围的逐渐缩小如图3所示。从原始范围[-α, α]到改进范围[-E, E]的过渡增强了ASMA在全局和局部搜索之间保持平衡的能力，促进了高质量解决方案的发现，同时避免了在局部最优中的停滞。

3.3.2. 后整合阶段

ASMA中的后整合机制旨在解决潜在的搜索停滞问题，类似于疟疾复发，其中休眠的寄生虫可能在治疗过早停止时重新活化。为了缓解过早收敛，ASMA引入了一种时间衰减的突变概率，由以下公式控制：

其中 促进初始多样性，而 确保在探索阶段的突变基线水平。衰减率由时间常数 控制，使从探索到开发的平滑过渡。在突变过程中，如果随机值小于M，则使用以下公式选择并更新一个随机解ID：

自适应方法通过式(17)得到加强，该公式控制个体解的选择用于突变（见式(13)）。本质上，这种频率自适应突变机制使ASMA能够平滑地从广泛探索转向集中开发，从而增强收敛稳定性并保持种群多样性。

3.3.3. 收敛标准

为了确保优化过程的有效终止，采用基于种群中惩罚目标值稳定性的收敛标准。当平均和最佳惩罚适应度值之间的差异变得足够小，表明几代中进展有限时，算法停止。收敛指标 定义为[43]：

其中 表示当前代中最低的惩罚目标值，而 表示种群中的平均值。其中 是用户指定的控制收敛敏感性的阈值。 的选择取决于问题特征。对于复杂或高维问题，建议使用较小的阈值 （例如， 到 ），以确保有意义的收敛检测[44]，[45]。该策略允许早期识别停滞并减少计算成本，同时保留解决方案质量。

3.3.4. ASMA的伪代码

ASMA的完整过程步骤在算法1中呈现。伪代码概述了核心阶段。

%Artemisinin slime mold algorithm (ASMA)
% Author and programmer: 
% Mr. Arnut Sutha
% Center of Excellence in Applied Mechanics and Structures, Department of Civil Engineering, Chulalongkorn University, 10330 Bangkok, Thailand
% e-mail:       mynut2009@gmail.com
% Researchgate: https://www.researchgate.net/profile/Arnut_Sutha
%_____________________________________________________________________________________________________    
function [bestfitness,Leader_pos,Convergence_curve]=ASMA(N,MaxFEs,lb,ub,dim,fobj) %SMA(N,Max_iter,lb,ub,dim,fobj)
% Initialization parameters
FEs=0;
it=1;
tFitness = inf*ones(1,N);%record the fitness of all slime mold
%% Initialization of the solution set
pop=initialization(N,dim,ub,lb);
%Calculate the fitness value of the initial solution set
for i=1:N
    Fitness(i)=fobj(pop(i,:));
    FEs=FEs+1;
end
[fmin,x]=min(Fitness);

%Container
New_pop=zeros(N,dim);
Fitnorm=zeros(1,N);
Convergence_curve=[];
%Record the current optimal solution
best=pop(x,:);
bestfitness=fmin;
time_per_iteration = zeros(1, MaxFEs); % Preallocate array for time per iteration
delta_threshold = 1e-4; % δ^t
%% Main loop
while FEs<=MaxFEs
    tic; % Start timer for iteration  
    [SmellOrder,SmellIndex] = sort(Fitness);  %Eq.(2.6)
    worstFitness = SmellOrder(N);
    bestFitness = SmellOrder(1);
    S=bestFitness-worstFitness+eps;  % plus eps to avoid denominator zero
    %calculate the fitness weight of each slime mold
    for i=1:N
        for j=1:dim
            if i<=(N/2)  %Eq.(2.5)
                weight(SmellIndex(i),j) = 1+rand()*log10((bestfitness-SmellOrder(i))/(S)+1);
            else
                weight(SmellIndex(i),j) = 1-rand()*log10((bestfitness-SmellOrder(i))/(S)+1);
            end
        end
    end
    K= 1-((FEs)^(1/6)/(MaxFEs)^(1/6));
    E=1*exp(-4*(FEs/MaxFEs)^(1/2));
    a = 1;
    b = 0.003;
    % c = 0.05;
    d = 0.01; %0.3;
     % IMF inertia weight calculation
    % t = FEs/MaxFEs; % normalized iteration number (0 to 1)
    %E1 = a * exp(-b * FEs) + d; % IMF inertia weight equation
    vb = unifrnd(-E,E,1,dim);  %Eq.(2.3)
    Fitness_threshold = mean(Fitness);  % ใช้ค่าเฉลี่ยของ Fitness เป็นเกณฑ์
    for i=1: N
        Fitnorm(i)= (Fitness(i)-min(Fitness))/(max(Fitness)-min(Fitness));
        for j=1:dim
            if rand<K
                if rand<0.5
                    New_pop(i,j) = pop(j,j) + E .* sin((2*pi)*rand()) * pop(i,j) * (-1)^FEs;
                else
                    New_pop(i,j) = pop(i,j) + E .* cos((2*pi)*rand()) * best(j) * (-1)^FEs;
                end
            else
               New_pop(i,j)=pop(i,j);
                % New_pop(i,j) = (ub-lb)*rand+lb;
            end
            if rand<Fitnorm(i)                
                AA=randperm(N);   %randperm(N) สุ่มเรียงลำดับตัวเลขทั้งหมดจาก 1 ถึง N โดยไม่ซ้ำ.
                New_pop(i,j)=pop(AA(3),j)+ vb(j)*(weight(i,j)*pop(AA(1),j)-pop(AA(2),j));
            else
                    New_pop(i,j) = vb(j)*pop(i,j);
            end
        end
      
        % Exponential Decrease Adaptive Mutation
        P_max = 0.5; % Initial mutation probability
        P_min = 0.1; % Minimum mutation probability
        tau = MaxFEs / log(P_max / P_min); % Time constant for exponential decrease
        % Exponential Decrease Adaptive Mutation
        mutation_prob = P_max * exp(-FEs / tau);% exp(-FEs/ tau)% Exponential decrease formula
        if rand() < mutation_prob
            mutation_index = randi([1 N]);
            New_pop(mutation_index, :) = lb + (ub - lb) .* rand(1, dim);
        end
        New_pop(i,:)=Mutation(New_pop(i,:),pop(i,:),best,dim);
        New_pop(i,:)=Transborder_reset(New_pop(i,:),ub,lb,dim,best);
        
        tFitness=fobj(New_pop(i,:));
        FEs=FEs+1;
        if tFitness<Fitness(i)
            pop(i,:)= New_pop(i,:);
            Fitness(i)=tFitness;
        end
    end
    [fmin,x]=min(Fitness);
    if fmin<bestfitness
        best=pop(x,:);
        bestfitness=fmin;
    end
    Convergence_curve(it)=bestfitness;
    Leader_pos=best;
    time_per_iteration(it) = toc; % End timer and record elapsed time
    disp(['Iteration ' num2str(it) ':   Best Cost(ASMA) = ' num2str(Convergence_curve(it)) ': Time =' num2str(time_per_iteration(it)) ]);
    
    %% Relative difference convergence check
    P_best = min(Fitness);        % Best penalized fitness
    P_mean = mean(Fitness);       % Mean penalized fitness
    delta = abs((P_mean - P_best)/P_best);  % Eq. (21)
    if delta < delta_threshold
        disp(['Terminated due to convergence: delta = ' num2str(delta) ', threshold = ' num2str(delta_threshold)]);
        break;
    end

    it=it+1;
end
disp(['Time(s) = ' num2str(sum(time_per_iteration))]);
end

function z=Mutation(z,x,b,dim)
for j=1:dim
    if rand<0.05
        z(j)=x(j);
    end
    if rand<0.2
        z(j)=b(j);
    end
end
end

% 
% function z=Transborder_reset(z,ub,lb,dim,best)
% for j=1:dim
%     if z(j)>ub || z(j)<lb
% 
%         z(j)=best(j);
% 
%     end
% end
% end

function z=Transborder_reset(z,ub,lb,dim,best)
    for j=1:dim
        if z(j) > ub || z(j) < lb
            z(j) = best(j);
        end
    end
end

%% This function initialize the first population of search agents
function Positions=initialization(SearchAgents_no,dim,ub,lb)
%Boundary_no = numel(ub); % number of boundaries
Boundary_no= size(ub,2); % numnber of boundaries

% If the boundaries of all variables are equal and user enter a signle
% number for both ub and lb
if Boundary_no==1
    Positions=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);
        Positions(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
    end
end
end

Sutha A, Tangaramvong S, Pyone E C, et al. Artemisinin slime mold algorithm for large-scale truss optimization under frequency constraints[C]//Structures. Elsevier, 2025, 80: 110045.

https://github.com/suthels/-/blob/main/README.md