SCI一区龙格库塔算法，原理详解，MATLAB代码免费获取

龙格库塔算法（Runge Kutta Method，RUN）是一种基于龙格库塔法的数学基础而提出的一种基于种群的无隐喻优化方法。提出的Runge Kutta优化器（RUN）被开发用于处理未来各种类型的优化问题。RUN利用由龙格库塔方法计算的斜率变化的逻辑作为全局优化的有前途的逻辑搜索机制。这种搜索机制得益于两个积极的探索和开发阶段，探索有前途的地区在特征空间和建设性的运动，向全球最佳解决方案。通过对50个数学测试函数和4个实际工程问题的测试，比较了RUN算法和其他元启发式算法的效率。RUN提供了非常有前途和有竞争力的结果，显示出上级探索和开发趋势，快速收敛速度和避免局部最优。

该成果于2021年发表在知名一区SCI期刊“Expert Systems with Applications”上。

RUN使用计算出的斜率作为搜索逻辑，以探索搜索空间中的有希望区域，并根据基于群的优化算法的逻辑构建一组用于种群集合进化的规则。

1、算法原理

（1）种群初始化：

在这一步中，逻辑是在允许的迭代次数内设置要进化的初始群。在RUN中，为一个大小为N的种群随机生成N个位置。种群中的每个成员xn (N = 1,2，…，N)是一个优化问题的维数为D的解。一般来说，初始位置是通过以下思路随机创建的:

式中，Ll和Ul分别为问题第l个变量(l = 1,2，…，D)的下界和上界，rand为[0,1]范围内的随机数。这个规则只是在一定范围内生成一些解。

（2）搜索机制的根源

采用RK方法确定了该RUN中的搜索机制。利用一阶导数定义系数k1，由下式计算。并且，由于应用位置目标函数需要相当长的计算时间，因此本文提出的优化算法使用位置xn而不是其适应度y(xn)。根据下式， xn +Δx和xn−Δx是xn的两个相邻位置，将y(x)视为最小化问题，xn−Δx和xn +Δx位置分别有最佳和最差位置。因此，要创建基于种群的算法，位置xn−Δx等于xb(即xb是xn周围的最佳位置)，而位置xn +Δx等于xw(即xw是xk周围的最差位置)。因此，定义k1为:

其中xw和xb是每次迭代得到的最坏解和最优解，它们是根据从总体(xr1, xr2, xr3)中随机选择的三个解确定的，并且r1∕r2∕r3∕= n。

为了增强探索搜索，产生随机行为，可以将上式改写为:

其中rand为[0,1]范围内的随机数。总的来说，最佳解决方案(xb)在寻找有前途的领域和向全球最佳解决方案迈进方面发挥着至关重要的作用。因此，在本研究中，使用一个随机参数(u)来增加优化过程中最优解(xb)的重要性。在下式中，Δx可表示为:

其中rand为[0,1]范围内的随机数。总的来说，最佳解决方案(xb)在寻找有前途的领域和向全球最佳解决方案迈进方面发挥着至关重要的作用。因此，在本研究中，使用一个随机参数(u)来增加优化过程中最优解(xb)的重要性。参数γ是由解空间大小决定的尺度因子，在优化过程中呈指数递减。Xavg是每次迭代中所有解决方案的平均值。

据此，其他三个系数(即k2、k3、k4)可分别写成:

其中，rand1和rand2为[0,1]范围内的两个随机数。

（3）更新解决方案

RUN算法从一组随机个体(解)开始优化过程。在每次迭代中，解使用RK方法更新它们的位置。为此，RUN使用由RK方法获得的解决方案和搜索机制。图2描述了位置如何使用RK方法更新其位置。在本研究中，全局探索的公式如下

局部开发的公式如下

其中μ为随机数，randn为正态分布的随机数。

（4）提高解决方案质量

在RUN算法中，增强的解决方案质量（ESQ）被用来提高解决方案的质量，并避免在每次迭代局部最优。通过应用ESQ，RUN算法确保每个解决方案都朝着更好的位置移动。在所提出的ESQ中，计算三个随机解的平均值（xavg），并将其与最佳位置（xb）组合以生成新的解（xnew1）。执行以下方案以使用ESQ创建解决方案（xnew2）：

其中 ， ， 

β是[0，1]范围内的随机数。c是随机数，w是随机数，其随着迭代次数的增加而减小。r是一个整数，可以是1、0或-1。xbest是迄今为止探索的最佳解决方案。根据上述方案，对于w < 1（即，后面的迭代），解Xnew 2趋向于创建开发搜索，而对于w ≥ 1（i.即，早期迭代），解Xnew 2趋向于进行探索搜索。

如下图所示，在RUN中考虑三条路径进行优化。该算法首先使用RK搜索机制生成位置xn+1，然后使用ESQ机制在搜索空间中探索有希望的区域。

RUN优化算法对应的伪代码如下

RUN优化算法对应的流程图如下

2、结果展示

3、MATLAB核心代码

%% 淘个代码 %%
% 微信公众号搜索：淘个代码，获取更多代码
% 龙格库塔算法（RUN）
function [Best_Cost,Best_X,Convergence_curve]=RUN(nP,MaxIt,lb,ub,dim,fobj)
Cost=zeros(nP,1);                % Record the Fitness of all Solutions
X=initialization(nP,dim,ub,lb);  % Initialize the set of random solutions   
Xnew2=zeros(1,dim);
Convergence_curve = zeros(1,MaxIt);
for i=1:nP
    Cost(i) = fobj(X(i,:));      % Calculate the Value of Objective Function
end
[Best_Cost,ind] = min(Cost);     % Determine the Best Solution
Best_X = X(ind,:);
Convergence_curve(1) = Best_Cost;
%% Main Loop of RUN 
it=1;%Number of iterations
while it<MaxIt
    it=it+1;
    f=20.*exp(-(12.*(it/MaxIt))); % (Eq.17.6) 
    Xavg = mean(X);               % Determine the Average of Solutions
    SF=2.*(0.5-rand(1,nP)).*f;    % Determine the Adaptive Factor (Eq.17.5)
    
    for i=1:nP
            [~,ind_l] = min(Cost);
            lBest = X(ind_l,:);   
            [A,B,C]=RndX(nP,i);   % Determine Three Random Indices of Solutions
            [~,ind1] = min(Cost([A B C]));
            % Determine Delta X (Eqs. 11.1 to 11.3)
            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,:);  % (Eq. 17.3)
            Xm = L.*Best_X+(1-L).*lBest;   % (Eq. 17.4)
              
            vec=[1,-1];
            flag = floor(2*rand(1,dim)+1);
            r=vec(flag);                   % An Interger number 
            
            g = 2*rand;
            mu = 0.5+.1*randn(1,dim);
            
            % Determine New Solution Based on Runge Kutta Method (Eq.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  
            
        % Check if solutions go outside the search space and bring them back
        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;               %(Eq.19-1)
            
            [A,B,C]=RndX(nP,i);
            Xavg=(X(A,:)+X(B,:)+X(C,:))/3;           %(Eq.19-2)         
            
            beta=rand(1,dim);
            Xnew1 = beta.*(Best_X)+(1-beta).*(Xavg); %(Eq.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
% End of ESQ         
%% Determine the Best Solution
        if Cost(i)<Best_Cost
            Best_X=X(i,:);
            Best_Cost=Cost(i);
        end
    end
% Save Best Solution at each iteration    
Convergence_curve(it) = Best_Cost;
disp(['it : ' num2str(it) ', Best Cost = ' num2str(Convergence_curve(it) )]);
end
end
% A funtion to determine a random number 
%with uniform distribution (unifrnd function in Matlab) 
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

参考文献

[1]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.

完整代码获取

后台回复关键词：

TGDM834

获取更多代码：

或者复制链接跳转：
https://docs.qq.com/sheet/DU3NjYkF5TWdFUnpu