基于MATLAB的麻雀搜索算法SSA代码解释（对照论文公式）

%% SSA:麻雀搜索算法 (doi:10.1080/21642583.2019.1708830)
%   Encoding format：utf-8
%   :param N:   种群数目
%   :param dim: 求解维度
%   :param x_min:   各维度搜索下限
%   :param x_max:   各维度搜索上限
%   :param iterate_max: 最大迭代次数
%   :param fitness: 适应度评价函数
%   :return:
%         optimal_value:  对应评价函数最优适应度
%         optimal_site: 最优位置
%% -----------------------------------------------------------------------------%%
function [optimal_value, optimal_site] = SSA(N, dim, x_min, x_max, iterate_max, fitness)
    %   初始化麻雀位置
    x = x_min + (x_max - x_min) .* rand(N, dim);
    %   最小常数
    varepsilon = exp(-16);
    %   安全阈值
    ST = 0.8;
    %   生产者数量
    PD = int64(0.2 * N);
    %   侦察麻雀数量
    SD = int64(0.1 * N);
    %   存储各个体适应度值
    Fx = ones(N, 1);
    %   求解对应位置的适应度值
    for i = 1:N
        Fx(i) = fitness(x(i, :));
    end
    %   初始化迭代起点
    iterate = 1;
    while iterate < iterate_max+1
        %   对适应度值进行排列并获取按照适应度从大到小的下标排列
        [~, Fx_range_index] = sort(Fx);
        %   获取当前全局最优位置
        x_best = x(Fx_range_index(1), :);
        f_g = Fx((Fx_range_index(1)));
        %   获取当前全局最差位置
        x_worst = x(Fx_range_index(N), :);
        f_w = Fx(Fx_range_index(N));
        %   用于记录新位置
        x_new = x;
        %   按照论文算法框架图，所有生产者是统一执行位置更新，见公式(3)
        if rand() < ST
            for i =1:PD
                x_new(Fx_range_index(i), :) = x(Fx_range_index(i), :) .* exp(-Fx_range_index(i) ./ (rand(1, dim) .* iterate_max));
            end
        else
            %   论文中的公式描述有点问题，前半部分描述的是指定了对应麻雀对应维度，即为标量，而后面乘的却又是向量，赋给的又是标量，因此在这里实现了两种写法
            %   在对Sphere函数测试中，发现第一种方式效果更好
            %   第一种：按照向量写法，即每个位置只有一个alpha
            x_new(Fx_range_index(1: PD), :) = x(Fx_range_index(1: PD), :) .* rand(PD, 1);
            %   第二种：按照标量写法，即每个位置的对应维度都有一个随机alpha
            %x_new(Fx_range_index(: PD), :) = x(Fx_range_index(: PD), :) .* np.random.random((PD, dim));
        end
        %   拾取者位置更新，公式(4)
        for i =PD+1:N
            if Fx_range_index(i) > N/2
                %   公式(4)和公式(3)一样，Q应该是对应维度都不同
                x_new(Fx_range_index(i), :) = rand(1, dim) .* exp((x_worst - x(Fx_range_index(i), :))/(Fx_range_index(i)^2));
            else
                %   A表示一个1×d的矩阵，其中每个元素随机分配1或-1
                A = rand(1, dim);
                A(A>0.5) = 1;
                A(A<0.5) = -1;
                A_inv = A' * (1/(A * A'));
                x_new(Fx_range_index(i), :) = x_best + abs(x(Fx_range_index(i), :) - x_best) * A_inv;
            end
        end
        %   侦察麻雀更新，公式(5):需要说明的是，在原文中并未说明侦察麻雀是哪些，只是描述了占总数的10%~20%，原文选取10%。这里就直接按照索引进行更新
        for i = 1:SD
            f_i = fitness(x_new(i, :));
            if f_i > f_g
                x_new(i, :) = x_best + rand(1, dim) .* abs(x_new(i, :) - x_best);
            elseif f_i == f_g
                k = -1 + rand() * 2;
                x_new(i, :) = x_new(i, :) + k .* (abs(x_new(i, :) - x_worst)./(f_i - f_w + varepsilon));
            end
        end
        %   如果当前位置要好于历史位置，则更新
        for i = 1:N
            if fitness(x_new(i, :)) < fitness(x(i, :))
                x(i, :) = x_new(i, :);
            end
        end
        %   对更新后的位置进行合法性检查
        x(x<x_min) = x_min;
        x(x>x_max) = x_max;
        %   求解对应位置的适应度值
        for i = 1:N
            Fx(i) = fitness(x(i, :));
        end
        %   迭代器++
        iterate = iterate + 1;
    end
    %   对适应度值进行排列并获取按照适应度从大到小的下标排列
    [~, Fx_range_index] = sort(Fx);
    %   获取当前最优解
    optimal_value = Fx(Fx_range_index(1));
    %   获取当前全局最优位置
    optimal_site = x(Fx_range_index(1), :);
end