​牛顿-拉夫逊优化算法算法，原理详解，MATLAB代码免费带走

牛顿-拉夫逊优化算法算法（Newton-Raphson-based optimizer，NRBO）是一种全新的元启发式优化方法，其灵感来源主要基于两个关键原理：Newton-Raphson搜索规则(NRSR)和陷阱避免算子(TAO)。NRSR使用Newton-Raphson方法来提高NRBO的探索能力，并提高收敛速度以达到改进的搜索空间位置。TAO有助于NRBO避免局部最优陷阱。NRBO具有进化能力强、搜索速度快、寻优能力强的特点。这一成果由Sowmya等人于2024年2月发表在中科院2区顶级SCI期刊《Engineering Applications of Artificial Intelligence》上。

1、算法原理 

（1）初始化 

想象一下以下场景：对具有以下定义的无约束单目标优化问题进行优化。f(x)的二次插值函数L(x)可表示为:

式中f(x)为待最小化适应度函数，xj为决策向量，dim为问题维数，lb为下界，u b为上界。与其他MH算法一样，NRBO通过在候选解的边界内产生初始随机种群来开始搜索最优解。基于有Np个种群，每个种群由模糊决策变量/向量组成的事实。因此，使用下式生成随机总体。

式中，xi j表示第n个总体的第j维位置，rand表示(0,1)之间的随机数。下式给出了可以描述所有维度的种群的种群矩阵。

（2）牛顿-拉夫逊搜索规则 

向量由NRSR控制，允许更准确地探索可行区域并获得更好的位置。NRM的提出是为了促进勘探趋势和加快收敛速度。由于许多优化技术是不可微的，因此在这种情况下，使用数学NRM来代替函数的显式公式。NRM从一个假定的初始解开始，沿着一个确定的方向前进到下一个位置。为了从下式中得到NRSR，必须使用TS来确定二阶导数。f(x−Δx)和f(x+Δx)的TS如下所示。

将上式相减或相加，可得f ' (x)和f ' (x)的表达式:

更新后的根位置重写如下:

考虑到NRSR应该是NRBO的主要组成部分，为了管理基于人口的搜索，有必要进行一定的调整。由于f(xn)是一个最小化问题，如图所示，位置xn+Δx的适应度值比位置x差，而位置xn - Δx的适应度值比位置xn高。

因此，NRBO将位置xn - Δx替换为位置Xb，因为Xb在其邻域的位置比位置xn好。而位置xn+Δx被位置Xw取代，Xw在其邻域的位置比位置xn差。建议的方法的另一个优点是它使用位置(xn)而不是其适应度(f(xn))，这节省了计算时间。然后，建议的NRSR表示为:

式中，randn为均值为0，方差为1的正态分布随机数，Xw为最差位置，Xb为最佳位置。该式可以通过帮助当前解更新其位置来增强当前解，通过随机参数增加了NRBO的搜索能力，更好地平衡了开采和勘探能力。

作为经验法则，所提出的算法必须能够在多样化和集约化之间达到平衡，以便在搜索空间中发现最优解并最终收敛到全局解。可以通过应用一个称为δ的自适应系数来增强算法。δ的表达式见下式。

其中IT表示当前迭代，而Max IT表示迭代的最大次数。为了保持勘探阶段和开发阶段之间的平衡，参数δ在迭代过程中自适应。根据上式，δ的值在1到−1之间变化。 

结合自适应参数δ，该算法在优化过程中考虑了随机行为，增加了多样性，避免了局部最优，同时显著减少了迭代次数，证明了NRBO算法的有效性，如下式所示。

通过加入另一个参数ρ来改进所提出的NRBO的利用，该参数将种群引向正确的方向。ρ的表达式如下。

其中a和b是(0,1)之间的随机数，r1和r2是从总体中随机选择的不同整数。然而，r1和r2的值并不相等。向量(Xn IT)的当前位置已由下式更新。

其中X1IT n是通过更新X1IT n得到的新向量位置。新提出的NRM进一步完善了NRSR。 

（3）陷阱避免操作员（TAO）

TAO是采用(ahmadanfar等人，2020)改进和增强的作业器。使用TAO可以显著改变xIT+ 1n的位置。通过结合最佳位置Xb和当前矢量位置xitn，它产生了一个具有增强质量的解决方案。使用下式，当rand的值小于DF时，得到解XIT TAO。

式中，rand为(0,1)之间的均匀随机数，θ1和θ2分别为(- 1,1)和(- 0.5,0.5)之间的均匀随机数，DF为控制NRBO性能的决定因子，μ1和μ2为随机数，分别由下式生成。

其中rand表示(0,1)之间的随机数，Δ表示(0,1)之间的随机数。将上式进一步简化如下:

其中β表示二进制数，1或0,rand表示随机数。如果Δ的值大于等于0.5，则β的值为0;否则，该值为1。由于参数μ1和μ2选择的随机性，使得种群更加多样化，从而摆脱了局部最优解的约束，从而提高了种群的多样性。与NRBO类似，GBO也受到牛顿方法的启发。因此，NRBO的概念可能看起来与GBO相似，但由于NRBO的独特特性，其性能比GBO要好得多。算法中给出了NRBO的伪代码。所提出的NRBO流程图如下图所示。所提出的NRBO伪代码如下所示。

所提出的NRBO伪代码如下所示。

2.结果展示

3.MATLAB核心代码

function [Best_Score, Best_Pos, CG_curve] = NRBO(N, MaxIt, LB, UB, dim, fobj)
    % Input arguments:
    % N     - Number of particles in the population
    % MaxIt - Maximum number of iterations
    % LB    - Lower bound of the search space
    % UB    - Upper bound of the search space
    % dim   - Dimensionality of the search space
    % fobj  - Objective function to minimize/maximize
    % Deciding Factor for Trap Avoidance Operator
    DF = 0.6;
    % Initialize the bounds for each dimension
    LB = ones(1, dim) * LB;           
    UB = ones(1, dim) * UB;
    % Initialization of the population
    Position = initialization(N, dim, UB, LB);
    Fitness = zeros(N, 1); % Vector to store individual costs
    % Calculate the initial fitness for each particle
    for i = 1:N
        Fitness(i) = fobj(Position(i,:));      
    end
    % Determine the best and worst fitness in the initial population
    [~, Ind] = sort(Fitness);     
    Best_Score = Fitness(Ind(1));
    Best_Pos = Position(Ind(1),:);
    Worst_Cost = Fitness(Ind(end));
    Worst_Pos = Position(Ind(end),:);
    % Initialize convergence curve
    CG_curve = zeros(1, MaxIt);
    % Main optimization loop
    for it = 1:MaxIt
        % Dynamic parameter delta, decreases over iterations
        delta = (1 - ((2 * it) / MaxIt)) .^ 5;
        % Loop over all particles in the population
        for i = 1:N                
            % Randomly select two different indices for differential evolution
            P1 = randperm(N, 2);                                       
            a1 = P1(1); a2 = P1(2);
            % Calculate the step size rho
            rho = rand * (Best_Pos - Position(i,:)) + rand * (Position(a1,:) - Position(a2,:));
            % Apply Newton-Raphson Search Rule
            Flag = 1;                   
            NRSR = SearchRule(Best_Pos, Worst_Pos, Position(i,:), rho, Flag);      
            X1 = Position(i,:) - NRSR + rho;                                  
            X2 = Best_Pos - NRSR + rho;                                            
            % Update position of particle
            Xupdate = zeros(1, dim);
            for j = 1:dim                                                                       
                X3 = Position(i,j) - delta * (X2(j) - X1(j));           
                a1 = rand; a2 = rand;
                Xupdate(j) = a1 * (a1 * X1(j) + (1 - a2) * X2(j)) + (1 - a2) * X3;             
            end
            % Trap Avoidance Operator to prevent local optima
            if rand < DF
                theta1 = -1 + 2 * rand(); theta2 = -0.5 + rand();      
                beta = rand < 0.5;
                u1 = beta * 3 * rand + (1 - beta); u2 = beta * rand + (1 - beta);          
                if u1 < 0.5
                    X_TAO = Xupdate +  theta1 * (u1 * Best_Pos - u2 * Position(i,:)) + theta2 * delta * (u1 * mean(Position) - u2 * Position(i,:));
                else
                    X_TAO = Best_Pos + theta1 * (u1 * Best_Pos - u2 * Position(i,:)) + theta2 * delta * (u1 * mean(Position) - u2 * Position(i,:));  
                end
                Xnew = X_TAO;
            else
                Xnew = Xupdate;
            end
            % Enforce boundary conditions
            Xnew = min(max(Xnew, LB), UB);
            % Evaluate new solution
            Xnew_Cost = fobj(Xnew);
            % Update the best and worst positions
            if Xnew_Cost < Fitness(i)
                Position(i,:) = Xnew;
                Fitness(i) = Xnew_Cost;
                % Update the global best solution
                if Fitness(i) < Best_Score
                    Best_Pos = Position(i,:);
                    Best_Score = Fitness(i);
                end
            end
            % Update the global worst solution
            if Fitness(i) > Worst_Cost
                Worst_Pos = Position(i,:);
                Worst_Cost = Fitness(i);
            end
        end
        % Update convergence curve
        CG_curve(it) = Best_Score;
        % Display iteration information
        disp(['Iteration ' num2str(it) ': Best Fitness = ' num2str(CG_curve(it))]);
    end
end

参考文献

[1] Sowmya R, Premkumar M, Jangir P. Newton-Raphson-based optimizer: A new population-based metaheuristic algorithm for continuous optimization problems[J]. Engineering Applications of Artificial Intelligence, 2024, 128: 107532.

完整代码获取方式：后台回复关键字：

TGDM877