龙格库塔算法 Runge Kutta Method及其Matlab代码

最新推荐文章于 2024-07-17 20:43:45 发布

Lzh_023016

最新推荐文章于 2024-07-17 20:43:45 发布

阅读量660

点赞数

分类专栏： Matlab 文章标签：算法 matlab 开发语言

本文链接：https://blog.csdn.net/Lzh_023016/article/details/132309101

版权

Matlab 专栏收录该内容

20 篇文章 2 订阅

订阅专栏

龙格库塔算法 Runge Kutta Method (RUN)

优化领域受到基于隐喻的“伪创新”或“花哨”优化器的影响。这些老套的方法大多模仿动物的搜索趋势，对优化过程本身的贡献很小。这些方法大多存在局部高效的性能、对简单问题的验证方法存在偏差、组件之间的交互高度相似等问题。在数学基础和数学中广为人知的Runge Kutta 方法的思想基础上，引入一种新的无隐喻的基于种群的优化方法。提出了RUNge - Kutta优化器(RUN)来处理未来各种类型的优化问题。
在这里插入图片描述
参考文献：Ahmadianfar I, Heidari A A, Gandomi A H, et al. RUN beyond the metaphor: An efficient optimization algorithm based on Runge Kutta method[J]. Expert Systems with Applications, 2021, 181: 115079.

龙格库塔算法的数学模型

RUN这种方法的好处是隐藏了优化器中使用的方程的真实性质。因此，RUN解释了Runge Kutta 技术的主要逻辑以及智能体群体基于种群的进化。包括：初始化，搜索机制，更新解，提高解的质量。
位置更新图：
在这里插入图片描述
搜索机制：

优化过程：

流程图：

龙格库塔算法的Matlab代码

%--------------------------------------------------------------------------
% RUN beyond the metaphor(RUN)
% Refernece:RUN beyond the metaphor: An efficient optimization algorithm based on Runge Kutta method
% Code by Luzhenhui Yangtze University
%--------------------------------------------------------------------------
% 输入
% nP：种群个体数目
% MaxIt：算法的最大迭代次数
% lb：变量的取值下界
% ub：变量的取值上界
% dim：变量的维度
% fob：优化的目标函数
%--------------------------------------------------------------------------
% 输出
% Best_Cost：最优适应度值
% Best_X：最优适应度值所对应的个体位置信息
% Convergence_curve：每次迭代的最优适应度值
%--------------------------------------------------------------------------
function [Best_Cost,Best_X,Convergence_curve]=RUN(nP,MaxIt,lb,ub,dim,fobj)
Cost=zeros(nP,1);                    % 记录种群个体的适应度值
X=initialization(nP,dim,ub,lb);      % 初始化种群个体   
Xnew2=zeros(1,dim);
Convergence_curve = zeros(1,MaxIt);  % 记录每一次迭代的最优适应度值
for i=1:nP
    Cost(i) = fobj(X(i,:));          % 计算种群个体的适应度值
end
[Best_Cost,ind] = min(Cost);         % 确定最优适应度值及其对应的最优个体
Best_X = X(ind,:);
Convergence_curve(1) = Best_Cost;
%% 算法主要步骤 
it=1;
while it<MaxIt
    it=it+1;
    f=20.*exp(-(12.*(it/MaxIt))); % 公式（17.6） 
    Xavg = mean(X);               % 种群个体的平均解
    SF=2.*(0.5-rand(1,nP)).*f;    % 通过公式（17.5）确定自适应系数SF   
    for i=1:nP
            [~,ind_l] = min(Cost);
            lBest = X(ind_l,:);   
            [A,B,C]=RndX(nP,i);   % 确定三个随机个体
            [~,ind1] = min(Cost([A B C]));
            % 通过公式（11.1）至（11.3）确定detla X
            gama = rand.*(X(i,:)-rand(1,dim).*(ub-lb)).*exp(-4*it/MaxIt);  
            Stp=rand(1,dim).*((Best_X-rand.*Xavg)+gama);
            DelX = 2*rand(1,dim).*(abs(Stp));
            % Determine Xb and Xw for using in Runge Kutta method
            if Cost(i)<Cost(ind1)                
                Xb = X(i,:);
                Xw = X(ind1,:);
            else
                Xb = X(ind1,:);
                Xw = X(i,:);
            end
            SM = RungeKutta(Xb,Xw,DelX);   % Search Mechanism (SM) of RUN based on Runge Kutta Method
            L=rand(1,dim)<0.5;
            Xc = L.*X(i,:)+(1-L).*X(A,:);  % 公式（17.3）
            Xm = L.*Best_X+(1-L).*lBest;   % 公式（17.4）
            vec=[1,-1];
            flag = floor(2*rand(1,dim)+1);
            r=vec(flag);                   
            g = 2*rand;
            mu = 0.5+.1*randn(1,dim);
            % 公式（18）
            if rand<0.5
                Xnew = (Xc+r.*SF(i).*g.*Xc) + SF(i).*(SM) + mu.*(Xm-Xc);
            else
                Xnew = (Xm+r.*SF(i).*g.*Xm) + SF(i).*(SM)+ mu.*(X(A,:)-X(B,:));
            end  
        % 边界处理
        FU=Xnew>ub;FL=Xnew<lb;Xnew=(Xnew.*(~(FU+FL)))+ub.*FU+lb.*FL; 
        CostNew=fobj(Xnew);
        if CostNew<Cost(i)
            X(i,:)=Xnew;
            Cost(i)=CostNew;
        end
%% Enhanced solution quality (ESQ)  (Eq. 19)      
        if rand<0.5
            EXP=exp(-5*rand*it/MaxIt);
            r = floor(Unifrnd(-1,2,1,1));
            u=2*rand(1,dim); 
            w=Unifrnd(0,2,1,dim).*EXP;               % 公式（19-1）
            [A,B,C]=RndX(nP,i);
            Xavg=(X(A,:)+X(B,:)+X(C,:))/3;           % 公式（19-2）         
            beta=rand(1,dim);
            Xnew1 = beta.*(Best_X)+(1-beta).*(Xavg); % 公式（19-3）
            for j=1:dim
                if w(j)<1 
                    Xnew2(j) = Xnew1(j)+r*w(j)*abs((Xnew1(j)-Xavg(j))+randn);
                else
                    Xnew2(j) = (Xnew1(j)-Xavg(j))+r*w(j)*abs((u(j).*Xnew1(j)-Xavg(j))+randn);
                end
            end
            FU=Xnew2>ub;FL=Xnew2<lb;Xnew2=(Xnew2.*(~(FU+FL)))+ub.*FU+lb.*FL;
            CostNew=fobj(Xnew2);            
            if CostNew<Cost(i)
                X(i,:)=Xnew2;
                Cost(i)=CostNew;
            else
                if rand<w(randi(dim)) 
                    SM = RungeKutta(X(i,:),Xnew2,DelX);
                    Xnew = (Xnew2-rand.*Xnew2)+ SF(i)*(SM+(2*rand(1,dim).*Best_X-Xnew2));  % (Eq. 20)
                    
                    FU=Xnew>ub;FL=Xnew<lb;Xnew=(Xnew.*(~(FU+FL)))+ub.*FU+lb.*FL;
                    CostNew=fobj(Xnew);
                    
                    if CostNew<Cost(i)
                        X(i,:)=Xnew;
                        Cost(i)=CostNew;
                    end
                end
            end
        end
%% 确定最优适应度值及对应的个体
        if Cost(i)<Best_Cost
            Best_X=X(i,:);
            Best_Cost=Cost(i);
        end

    end
% 保存每一次迭代的最优适应度值  
Convergence_curve(it) = Best_Cost;
disp(['it : ' num2str(it) ', Best Cost = ' num2str(Convergence_curve(it) )]);
end
end

function z=Unifrnd(a,b,c,dim)
a2 = a/2;
b2 = b/2;
mu = a2+b2;
sig = b2-a2;
z = mu + sig .* (2*rand(c,dim)-1);
end

% A function to determine thress random indices of solutions
function [A,B,C]=RndX(nP,i)
Qi=randperm(nP);Qi(Qi==i)=[];
A=Qi(1);B=Qi(2);C=Qi(3);
end

function SM=RungeKutta(XB,XW,DelX)
dim=size(XB,2);
C=randi([1 2])*(1-rand);
r1=rand(1,dim);
r2=rand(1,dim);
K1 = 0.5*(rand*XW-C.*XB);
K2 = 0.5*(rand*(XW+r2.*K1.*DelX/2)-(C*XB+r1.*K1.*DelX/2));
K3 = 0.5*(rand*(XW+r2.*K2.*DelX/2)-(C*XB+r1.*K2.*DelX/2));
K4 = 0.5*(rand*(XW+r2.*K3.*DelX)-(C*XB+r1.*K3.*DelX));
XRK = (K1+2.*K2+2.*K3+K4);
SM=1/6*XRK;
end


function Positions=initialization(SearchAgents_no,dim,ub,lb)
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