SCI一区梯度优化算法，原理详解，MATLAB代码免费获取

基于梯度的优化算法（Gradient-based optimizer，GBO）是一种元启发式的数学优化算法。该算法受基于梯度的牛顿方法的启发，使用两个主要算子：梯度搜索规则（GSR）和局部逃逸算子（LEO）以及一组向量来探索搜索空间。GSR采用基于梯度的方法来增强搜索趋势，加快收敛速度，以获得更好的搜索空间中的位置。LEO使得所提出的GBO能够逃离局部最优。

该成果于2020年发表在计算机领域一区SCI期刊“Information Sciences”上，目前被引高达51960，是计算机领域十分热门的优化技术。

1、算法原理 

（1）种群初始化： 优化问题涉及一组决策变量、约束和目标函数。GBO的控制参数包括用于从勘探到开采的过渡的参数（a）和概率。迭代次数和种群大小取决于问题的复杂性。在该算法中，种群中的每个成员都被称为“向量”。因此，GBO包括D维搜索空间中的N个向量。因此，向量可以表示为：

通常，GBO的初始向量在D维搜索域中随机生成，其可以定义为：

其中Xmin和Xmax是决策变量X的边界，rand(0,1)是[0，1]中的随机数。

（2）梯度搜索规则

在梯度搜索规则中，控制向量的移动以更好地在可行域中搜索并获得更好的位置。为了增强GBO的勘探趋势和加速收敛，在GB方法的基础上提出了GSR方法。然而，这个规则是从牛顿的基于梯度的方法中提取的。考虑到许多优化问题是不可微的，数值梯度法被用来代替函数的直接求导。通常，GB方法从猜测的初始解开始，并沿梯度指定的方向沿着向下一个位置移动。根据下式，一阶导数必须通过利用泰勒级数来计算。函数f(x+∆x)和f(x-∆x)的泰勒级数可以分别表示为：

一阶导数由以下中心差分公式给出：

基于以上公式，则新位置（xn+1）被定义为：

由于GSR被认为是所提出的算法的主要核心，一些修改是必不可少的处理基于人口的搜索。xn的相邻位置是如图所示的xn+∆x和xn-∆x。在GBO算法中，这些相邻位置被替换为群体中的两个其他位置（向量）。由于f(x)是一个最小化问题，如图所示，位置xn+∆x的适应度比xn-∆x差。因此，GBO算法用在位置xn的邻域中具有较好位置的xbest替换位置xn-∆x，而用在xn的邻域中的较差位置的xworst替换xn+∆x。此外，所提出的算法采用位置（xn），而不是其适应度f(xn)，因为使用位置的适应度在计算中更耗时。然后，GSR的公式如下：

（3）寻找食物

通常，优化算法应该能够平衡全局探索和局部开发，以探索搜索域中有希望的区域，并最终收敛到全局最优解。为了实现这一目标，可以通过利用自适应系数来改变GSR。在这项研究中，下式被引入作为GBO中最重要的参数，以平衡勘探和开采搜索过程，并且它可以表示为：

其中βmin和βmax分别为0.2和1.2，m是迭代次数，M是迭代总数。下图描绘了参数a如何随着迭代次数而变化。最大迭代次数为1000。此参数可以在每次迭代时更改。它在迭代初期具有较大的值以增强种群多样性，然后随着迭代次数的增加而减小以加速收敛。生成的解决方案应该能够探索其相应的最佳解决方案周围的搜索空间。在这方面，参数值在迭代次数从550到750的范围内增加，这有助于所提出的算法摆脱任何局部最优，因为它可以增加种群的多样性，以围绕所获得的最佳解进行搜索。因此，可得到公式：

提出的GSR有助于GBO考虑优化过程中的随机行为，促进探索和逃避局部最优。在式中，∆x是根据最佳解(xbest)与随机选择位置(xm r1)之间的差值确定的。为了保证∆x在每次迭代中都是变化的，参数d由下式定义。

其中rand（1:N）为N维随机数，r1; r2; r3;和r4是从[1，N]中随机选取的不同整数，step是一个步长，由xbest和xm r1决定。GSR可以重写为：

为了更好地利用xn附近的区域，还添加了移动方向(DM)。这一项使用最佳矢量，并将当前矢量(xn)向(xbest-xn)方向移动。因此，该过程产生了合适的局部搜索倾向，提高了GBO算法的收敛速度。建议的DM公式如下:

其中rand是[0,1]中的一个随机数，q2是一个随机参数，这有助于每个向量具有不同的步长。

下图显示当前位置是如何更新的。如图所示，在搜索空间中由GSR和DM指定的随机点上创建位置X1m n。

这种搜索方向方法强调开发过程。下式表示的搜索方法适合局部搜索，但对全局搜索有限制。因此，GBO利用了这两种搜索方法加强勘探开发。因此，根据位置X1m n, X2m n和当前位置(Xm n)，下一次迭代(xm+1 n)的新解可以定义为:

其中ra和rb为[0,1]中的两个随机数。

下图描述了向量如何在二维搜索空间中更新其相对于X1m n、X2m n和X3m n的位置。

（4）局部转义算子

为了提高GBO算法求解复杂问题的效率，引入了LEO算法。该运算符可以显著地改变解xm+1 n的位置。LEO通过使用若干解来生成具有上级性能的解（Xm LEO），所述若干解包括最佳位置（xbest）、解X1m n和X2m n、两个随机解xm r1和xm r2以及新的随机生成的解（xm k）。f1是[-1，1]范围内的均匀随机数，f2是来自正态分布的随机数，平均值为0，标准差为1，pr是概率，u1，u2和u3是三个随机数，定义为：

其中rand为[0,1]范围内的随机数。以上方程可以化简为

其中L1是一个值为0或1的二进制参数。为确定式中的解xm k，采用以下公式:

其中Xrand是新解，Xmp是群体，u2是一个在[0,1]范围内的随机数。对应的简化公式为：

其中L2是一个值为0或1的二进制参数。如果μ2小于0.5，则μ2的值为1，否则为0。这种选择参数u1、u2和u3值的随机行为有助于增加种群的多样性，避免局部最优解。GBO算法伪代码如表所示。

2、结果展示

3、MATLAB核心代码

%% 淘个代码 %%
% 微信公众号搜索：淘个代码，获取更多代码
% 基于梯度的优化算法（Gradient-based optimizer，GBO）
function [Best_Cost,Best_X,Convergence_curve]=GBO(nP,MaxIt,lb,ub,dim,fobj)
%% Initialization
nV = dim;                        % Number f Variables
pr = 0.5;                        % Probability Parameter
lb = ones(1,dim).*lb;            % lower boundary 
ub = ones(1,dim).*ub;            % upper boundary
Cost = zeros(nP,1);
X = initialization(nP,nV,ub,lb); %Initialize the set of random solutions
Convergence_curve = zeros(1,MaxIt);
for i=1:nP
    Cost(i) = fobj(X(i,:));      % Calculate the Value of Objective Function
end
[~,Ind] = sort(Cost);     
Best_Cost = Cost(Ind(1));        % Determine the vale of Best Fitness
Best_X = X(Ind(1),:);
Worst_Cost=Cost(Ind(end));       % Determine the vale of Worst Fitness
Worst_X=X(Ind(end),:);
%% Main Loop
for it=1:MaxIt
    beta = 0.2+(1.2-0.2)*(1-(it/MaxIt)^3)^2;                        % Eq.(14.2)
    alpha = abs(beta.*sin((3*pi/2+sin(3*pi/2*beta))));              % Eq.(14.1)
    for i=1:nP                
        A1 = fix(rand(1,nP)*nP)+1;                                  % Four positions randomly selected from population        
        r1 = A1(1);r2 = A1(2);   
        r3 = A1(3);r4 = A1(4);        
        Xm = (X(r1,:)+X(r2,:)+X(r3,:)+X(r4,:))/4;                   % Average of Four positions randomly selected from population        
        ro = alpha.*(2*rand-1);ro1 = alpha.*(2*rand-1);        
        eps = 5e-3*rand;                                            % Randomization Epsilon
        DM = rand.*ro.*(Best_X-X(r1,:));Flag = 1;                   % Direction of Movement Eq.(18)
        GSR=GradientSearchRule(ro1,Best_X,Worst_X,X(i,:),X(r1,:),DM,eps,Xm,Flag);      
        DM = rand.*ro.*(Best_X-X(r1,:));
        X1 = X(i,:) - GSR + DM;                                     % Eq.(25)
        DM = rand.*ro.*(X(r1,:)-X(r2,:));Flag = 2;
        GSR=GradientSearchRule(ro1,Best_X,Worst_X,X(i,:),X(r1,:),DM,eps,Xm,Flag); 
        DM = rand.*ro.*(X(r1,:)-X(r2,:));
        X2 = Best_X - GSR + DM;                                     % Eq.(26)            
        Xnew=zeros(1,nV);
        for j=1:nV                                                  
            ro=alpha.*(2*rand-1);                       
            X3=X(i,j)-ro.*(X2(j)-X1(j));           
            ra=rand;rb=rand;
            Xnew(j) = ra.*(rb.*X1(j)+(1-rb).*X2(j))+(1-ra).*X3;     % Eq.(27)          
        end
        % Local escaping operator(LEO)                              % Eq.(28)
        if rand<pr           
            k = fix(rand*nP)+1;
            f1 = -1+(1-(-1)).*rand();f2 = -1+(1-(-1)).*rand();         
            ro = alpha.*(2*rand-1);
            Xk = unifrnd(lb,ub,1,nV);%lb+(ub-lb).*rand(1,nV);       % Eq.(28.8)
            L1=rand<0.5;u1 = L1.*2*rand+(1-L1).*1;u2 = L1.*rand+(1-L1).*1;
            u3 = L1.*rand+(1-L1).*1;                                    
            L2=rand<0.5;            
            Xp = (1-L2).*X(k,:)+(L2).*Xk;                           % Eq.(28.7)
            if u1<0.5
                Xnew = Xnew + f1.*(u1.*Best_X-u2.*Xp)+f2.*ro.*(u3.*(X2-X1)+u2.*(X(r1,:)-X(r2,:)))/2;     
            else
                Xnew = Best_X + f1.*(u1.*Best_X-u2.*Xp)+f2.*ro.*(u3.*(X2-X1)+u2.*(X(r1,:)-X(r2,:)))/2;   
            end
        end
        % Check if solutions go outside the search space and bring them back
        Flag4ub=Xnew>ub;
        Flag4lb=Xnew<lb;
        Xnew=(Xnew.*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;                
        Xnew_Cost = fobj(Xnew);
        % Update the Best Position        
        if Xnew_Cost<Cost(i)
            X(i,:)=Xnew;
            Cost(i)=Xnew_Cost;
            if Cost(i)<Best_Cost
                Best_X=X(i,:);
                Best_Cost=Cost(i);
            end            
        end
       % Update the Worst Position 
        if Cost(i)>Worst_Cost
            Worst_X= X(i,:);
            Worst_Cost= Cost(i);
        end               
    end
    
    Convergence_curve(it) = Best_Cost;
    % Show Iteration Information
    disp(['Iteration ' num2str(it) ': Best Fitness = ' num2str(Convergence_curve(it))]);
end
end
% _________________________________________________
% Gradient Search Rule
function GSR=GradientSearchRule(ro1,Best_X,Worst_X,X,Xr1,DM,eps,Xm,Flag)
    nV = size(X,2);
    Delta = 2.*rand.*abs(Xm-X);                            % Eq.(16.2)
    Step = ((Best_X-Xr1)+Delta)/2;                         % Eq.(16.1)
    DelX = rand(1,nV).*(abs(Step));                        % Eq.(16)
    GSR = randn.*ro1.*(2*DelX.*X)./(Best_X-Worst_X+eps);   % Gradient search rule  Eq.(15)
    if Flag == 1
      Xs = X - GSR + DM;                                   % Eq.(21)
    else
      Xs = Best_X - GSR + DM;
    end    
    yp = rand.*(0.5*(Xs+X)+rand.*DelX);                    % Eq.(22.6)
    yq = rand.*(0.5*(Xs+X)-rand.*DelX);                    % Eq.(22.7)
    GSR = randn.*ro1.*(2*DelX.*X)./(yp-yq+eps);            % Eq.(23)   
end

参考文献

[1]Ahmadianfar I, Bozorg-Haddad O, Chu X. Gradient-based optimizer: A new metaheuristic optimization algorithm[J]. Information Sciences, 2020, 540: 131-159.

完整代码获取

后台回复关键词：

TGDM834

获取更多代码：

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