螢火蟲算法

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\)的吸引力可表示為：

\[\beta_{ij}(r_{ij}) = \beta_0e^{-\gamma{r_{ij}^2}} \]

其中，\(\beta_0\)為最大吸引力，當距離\(r=0\)時，吸引力最大。通常\(\beta_0=1\),\(\gamma\in{[0.01,100]}\).

3.4 螢火蟲位置更新公式

螢火蟲\(i\)吸引着螢火蟲\(j\)，因此螢火蟲\(j\)的位置更新公式：

\[x_j(t+1) = x_j(t) + \beta_{ij}(r_{ij})(x_i(t) - x_j(t)) + \alpha\varepsilon_j \]

其中，\(t\)為算法的迭代次數；\(x_i\)、\(x_j\)分別為螢火蟲\(i\)和螢火蟲\(j\)所處的空間位置；\(\alpha\in{[0,1]}\), \(\varepsilon\)是高斯分布得到的隨機向量。
第一版本：

% =========================================================================
%     Coder: Lee WenTsao
%      Time: 2022-05-14
%     Email: [e-mail]
% Reference: Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,      
%            Luniver Press, First Edition.
% =========================================================================

%% 清理運行環境
clc
clear 
close all

%% 問題定義
option = 2;                          % 選擇優化函數
dimension = 2;                       % 維數
[fobj, bound] = Optimizer(option);
lb = bound(1);                       % 下界
ub = bound(2);                       % 上界

%% 繪制雲圖
figure(1)
x = linspace(lb, ub, 101);
y = linspace(lb, ub, 101);
z = zeros(101);
for i=1:length(x)
    for j=1:length(y)
        z(i,j) = fobj([x(i), y(j)]);
    end
end

contour(x,y,z);                      % 函數雲圖顯示
hold on;

%% 螢火蟲參數設置
num_pop = 30;      % 初始種群數目
Max_gen = 1000;    % 最大迭代次數

beta_0 = 1.0;      % 最大吸引力
gamma = 0.1;       % 光強吸收系數

alpha = 1.0;       % 步長因子
theta = 0.97;      % alpha衰減因子

%% 初始化種群
Sol = zeros(num_pop, dimension);
I = zeros(num_pop, 1);
for i=1:num_pop
    Sol(i,:) = lb + (ub - lb)*rand(1,dimension);  % 初始化種群
    I(i) = fobj(Sol(i,:));                        % 適應度 
end        

[fmin, id] = min(I);
best_scores = [fmin];

%% 動態表示
points = scatter(Sol(:,1),Sol(:,2),"ro","filled");
xlim([lb,ub]);
ylim([lb,ub]);
drawnow;


%% 仿真
for iter=1:Max_gen
    alpha = alpha*theta;
    scale = abs((ub -lb)*ones(1, dimension));   % 優化問題的尺度
    for i=1:num_pop
        for j=1:num_pop
            if I(i)>I(j)
                % 計算螢火蟲i和螢火蟲j之間的距離
                r = sqrt(sum((Sol(i,:) - Sol(j,:)).^2)); 

                % 計算螢火蟲之間的吸引力
                beta = beta_0*exp(-gamma*r.^2);

                % 搜索精度
                steps = alpha.*(rand(1,dimension) - 0.5).*scale;

                % 更新螢火蟲的位置
                S = Sol(i,:) + beta*(Sol(j,:) - Sol(i,:)) + steps; 
                
                % 螢火蟲越界處理
                Tp = S>ub;
                Tm = S<lb;
                S = S.*(~(Tp+Tm)) + ub.*Tp + lb.*Tm;

                % 適應度值
                new_fun = fobj(S);

                % 進化機制
                if new_fun<I(i)
                    Sol(i,:) = S;
                    I(i) = new_fun;
                end

                % 更新最優
                if I(i)<fmin
                    fmin = I(i);
                    id = i;
                end

            end
        end
    end
    best_scores = [best_scores,fmin];
    
    reset(points);
    points = scatter(Sol(:,1), Sol(:,2), 'ro','filled');
    pause(0.1)

    %% 輸出
    if ~mod(iter,50)
        disp(['迭代次數：' num2str(iter) '||   最優值：' num2str(fmin)]);
    end
end

%% 收斂曲線可視化
figure(2)
xx = 1:50:1001;
xxx = 0:50:1000;
plot(xxx,log10(best_scores(xx)),'r','LineWidth',1.5);
hold on;
sz = 40;
scatter(xxx(2:end-1),log10(best_scores(xx(2:end-1))),sz,'ro','filled');
xlabel('$$Iteration$$','Interpreter','latex');
ylabel('$$\log_{10} (fitness)$$','Interpreter','latex');
title('Convergence curve')

xticks(0:50:1000);

楊老師原始版本：

% ---------------------------------------------------------------------- %
% The Firefly Algorithm (FA) for unconstrained function optimization     %
% by Xin-She Yang (Cambridge University) @2008-2009                      %
% Programming dates: 2008-2009, then revised and updated in Oct 2010     %
% ---------------------------------------------------------------------- %

% References -- citation details: -------------------------------------- %
% (1) Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,            %
%     Luniver Press, First Edition, (2008).                              %
% (2) Xin-She Yang, Firefly Algorithm, Stochastic Test Functions and     %
%     Design Optimisation, Int. Journal of Bio-Inspired Computation,     %
%     vol. 2, no. 2, 78-84 (2010).                                       %
% ---------------------------------------------------------------------- %

% -------- Start the Firefly Algorithm (FA) main loop ------------------ % 
function fa_ndim_new 
n=20;                   % Population size (number of fireflies)
alpha=1.0;              % Randomness strength 0--1 (highly random)
beta0=1.0;              % Attractiveness constant
gamma=0.01;             % Absorption coefficient
theta=0.97;             % Randomness reduction factor theta=10^(-5/tMax) 
d=10;                   % Number of dimensions
tMax=500;               % Maximum number of iterations
Lb=-10*ones(1,d);       % Lower bounds/limits
Ub=10*ones(1,d);        % Upper bounds/limits
% Generating the initial locations of n fireflies
for i=1:n,
   ns(i,:)=Lb+(Ub-Lb).*rand(1,d);         % Randomization
   Lightn(i)=cost(ns(i,:));               % Evaluate objectives
end

%%%%%%%%%%%%%%%%% Start the iterations (main loop) %%%%%%%%%%%%%%%%%%%%%%%
for k=1:tMax,        
 alpha=alpha*theta;     % Reduce alpha by a factor theta
 scale=abs(Ub-Lb);      % Scale of the optimization problem
% Two loops over all the n fireflies
for i=1:n,
    for j=1:n,
      % Evaluate the objective values of current solutions
      Lightn(i)=cost(ns(i,:));           % Call the objective
      % Update moves
      if Lightn(i)>=Lightn(j),           % Brighter/more attractive
         r=sqrt(sum((ns(i,:)-ns(j,:)).^2));
         beta=beta0*exp(-gamma*r.^2);    % Attractiveness
         steps=alpha.*(rand(1,d)-0.5).*scale;
      % The FA equation for updating position vectors
         ns(i,:)=ns(i,:)+beta*(ns(j,:)-ns(i,:))+steps;
      end
   end % end for j
end % end for i

% Check if the new solutions/locations are within limits/bounds
ns=findlimits(n,ns,Lb,Ub);
%% Rank fireflies by their light intensity/objectives
[Lightn,Index]=sort(Lightn);
nsol_tmp=ns;
for i=1:n,
 ns(i,:)=nsol_tmp(Index(i),:);
end
%% Find the current best solution and display outputs
fbest=Lightn(1), nbest=ns(1,:)
end % End of the main FA loop (up to tMax) 

% Make sure that new fireflies are within the bounds/limits
function [ns]=findlimits(n,ns,Lb,Ub)
for i=1:n,
  nsol_tmp=ns(i,:);
  % Apply the lower bound
  I=nsol_tmp<Lb;  nsol_tmp(I)=Lb(I);
  % Apply the upper bounds
  J=nsol_tmp>Ub;  nsol_tmp(J)=Ub(J);
  % Update this new move
  ns(i,:)=nsol_tmp;
end

%% Define the objective function or cost function
function z=cost(x)
% The modified sphere function: z=sum_{i=1}^D (x_i-1)^2
z=sum((x-1).^2); % The global minimum fmin=0 at (1,1,...,1)
% -----------------------------------------------------------------------%