萤火虫算法

1. 萤火虫优化算法背景

受萤火虫发光强度的启发，2008年，英国剑桥大学学者Xin-She Yang提出萤火虫算法(Firefly Algorithm, FA)。自然界中，萤火虫可以发出短促、有节奏的闪光。通常这种闪光仅在一定范围内可见。萤火虫通过闪光可以吸引异性和猎取食物。为了使算法更加简单，该算法只考虑了萤火虫强度的变化和吸引力这两个因素。

2. 萤火虫优化算法理想化数学模型

依照萤火虫发光的特性，给出以下理想化规则:
(1) 萤火虫不分雌雄，每个萤火虫都会被比它发光更亮的萤火虫吸引；
(2) 吸引力与发光强度成正比；
(3) 萤火虫的亮度由目标函数值决定。

3. 萤火虫优化算法的更新过程

3.1 绝对亮度的定义

为了表示萤火虫(i)的亮度随距离(r)的变化，定义如下绝对亮度：
萤火虫(i)绝对亮度为距离(r=0)时的亮度，记为(I_i).
注意：为了降低算法的复杂度，假定萤火虫(i)的绝对亮度(I_i)与(x_i)的目标函数值相等。

3.2 相对亮度的定义

为了表示萤火虫(i)对萤火虫(j)的吸引大小，定义如下相对亮度：
萤火虫(i)在萤火虫(j)位置的光强度，记为(I_{ij})

[I_{ij}(r_{ij}) = I_ie^{-gamma{r_{ij}^2}} ]

其中，(gamma)为光吸收系数，(r_{ij})为萤火虫(i)到萤火虫(j)的距离.

3.3 吸引力的定义

假设萤火虫(i)对萤火虫(j)的吸引力和萤火虫(i)对萤火虫(j)的相对亮度成比例，所以萤火虫(i)对萤火虫(j)的吸引力可表示为：

[eta_{ij}(r_{ij}) = eta_0e^{-gamma{r_{ij}^2}} ]

其中，(eta_0)为最大吸引力，当距离(r=0)时，吸引力最大。通常(eta_0=1),(gammain{[0.01,100]}).

3.4 萤火虫位置更新公式

萤火虫(i)吸引着萤火虫(j)，因此萤火虫(j)的位置更新公式：

[x_j(t+1) = x_j(t) + eta_{ij}(r_{ij})(x_i(t) - x_j(t)) + alphavarepsilon_j ]

其中，(t)为算法的迭代次数；(x_i)、(x_j)分别为萤火虫(i)和萤火虫(j)所处的空间位置；(alphain{[0,1]}), (varepsilon)是高斯分布得到的随机向量。

% -------------------------------------------------------------------------
% 名  称：萤火虫算法
% 代  码：编码雪人
% 时  间：2021-6-5
% -------------------------------------------------------------------------

%% 清空运行环境
clc
clear
close all

%% 定义目标函数
fobj = @ sphere;

%% 参数定义
n = 30;  % 种群数目
d = 2;   % 种群维数

beta0 = 1.0;             % 最大吸引力
gamma = 0.01;            % 光强吸收系数
alpha = 0.2;             % 步长因子
theta = 0.97;            % alpha衰减因子

Max_iter = 1000;        % 最大迭代步数
ub =  100*ones(1, d);  % 变量的上界
lb = -100*ones(1, d);  % 变量的下界

%% 初始化种群和目标函数值作为各自最大发光亮度
I = zeros(n, 1);
X = zeros(n, d);
for p=1:n
    X(p, :) = lb + (ub - lb).*rand(1, d);  % 初始化萤火虫的位置
    I(p) = fobj(X(p, :));                  % 最大亮度I_0
end
[Y, idx] = sort(I);
X = X(idx, :);
X1 = X;  % 复制一份X
Y1 = Y;  % 复制一份Y

%% 算法寻优
for it=1:Max_iter
    for i=1:n
        for j=1:n 
            if Y(i)>Y1(j)
                % 计算萤火虫i与萤火虫j之间的距离
                r = sqrt(sum((X(i, :) - X1(j, :)).^2, 2));
                
                % 计算萤火虫i与萤火虫j的吸引力
                beta = beta0*exp(-gamma*r^2);
                
                % 搜索精度
                steps = alpha.*(rand(1, d) - 0.5);
                
                % 更新萤火虫的位置
                X(i, :) = X(i, :) + beta*(X1(j, :) - X(i, :)) + steps;
                
                % 越界处理
                X(i, :) = simplebound(X(i, :), ub, lb);
            end 
        end
    end
    % 计算更新后的萤火虫亮度
    I = fobj(X);
        
    % 对萤火虫亮度进行排序
    [Y, idx] = sort(I);
    X = X(idx, :);
    X1 = X;       % 复制一份X
    Y1 = Y;       % 复制一份Y
    
    % 输出
    bestScore = min(Y(1));
    disp(['迭代次数：' num2str(it) '最优值：' num2str(bestScore)])
end
%% 越界函数处理
function x = simplebound(x, ub, lb)
    d1 = length(x);
    for i=1:d1
        if x(i)>ub(i)
            x(i) = ub(i);
        elseif x(i)<lb(i)
            x(i) = lb(i);
        end
    end
end

%%  目标函数

function  z  = sphere(x)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
    z = sum(x.^2,2);
end