萤火虫算法（Firefly Algorithm）

1.萤火虫的行为

萤火虫闪烁的光芒在热带和温带地区的夏季天空中是一道令人惊叹的风景。大约有两千种萤火虫，大多数萤火虫会发出短暂而有节奏的闪光。闪光的模式对于特定物种来说往往是独一无二的。闪光是由生物发光过程产生的，这种信号系统的真正功能仍在争论中。然而，这种闪光的两个基本功能是吸引交配伙伴（交流）和吸引潜在猎物。此外，闪光还可以作为一种保护性警告机制。

有节奏的闪光、闪光的频率和时间构成了将两性联系在一起的信号系统的一部分。在同一物种中，雌性萤火虫对雄性萤火虫独特的闪光模式做出反应，而在某些物种中，如photuris（萤火虫中的一种肉食动物），雌性萤火虫可以模仿其他物种的交配闪光模式，从而引诱并吃掉可能误认为闪光是潜在合适配偶的雄性萤火虫。

我们知道，距离光源特定距离处的光强度服从平方反比定律。也就是说，光照强度 $I$ 随着距离 $r$ 的增加而减小。此外，空气吸收的光随着距离的增加而变得越来越弱。这两个综合因素使得大多数萤火虫只在有限的距离内可见，通常在夜间几百米，这通常足以让萤火虫交流。

2.萤火虫算法（FA）

（1）萤火虫算法的三个理想状态

a.所有萤火虫都是中性的，因此任意一只萤火虫会被其他萤火虫所吸引，而不管它们的性别（all fireflies are unisex so that one firefly will be attracted to other fireflies regardless of their sex）；

b.吸引力与它们的亮度成正比，因此，对于任何两个闪烁的萤火虫来说，亮度越低的萤火虫会向亮度越高的萤火虫移动。吸引力与亮度成正比，两者都随着距离的增加而减小。如果没有比特定萤火虫更亮的萤火虫，它会随机移动(Attractiveness is proportional to their brightness, thus for any two flashing fireflies, the less brighter one will move towards the brighter one. The attractiveness is proportional to the brightness and they both decrease as their distance increases. If there is no brighter one than a particular firefly, it will move randomly)；

c.萤火虫的亮度受目标函数的地形影响或决定。对于最大化问题，亮度可以简单地与目标函数的值成正比（The brightness of a firefly is affected or determined by the landscape of the objective function. For a maximization problem, the brightness can simply be proportional to the value of the objective function）。

（2）伪代码

基于这三条规则，萤火虫算法（FA）的基本步骤可以总结为以下所示的伪代码。

（3）吸引力

在萤火虫算法中，有两个重要的问题：光强度的变化和吸引力的公式。为了简单起见，我们可以假设萤火虫的吸引力是由其亮度决定的，而亮度又与编码的目标函数有关。

在最大优化问题的最简单情况下，萤火虫在特定位置 $x$ 的亮度 $I$ 可以选择为 $I_{(x)}$ ∝ $f_{(x)}$ 。然而，吸引力 $\beta$ 是相对的，它应该在旁观者的眼睛中看到，或者由其他萤火虫判断。因此，它会随着萤火虫 $i$ 和萤火虫 $j$ 之间的距离 $r_{ij}$ 而变化。此外，光强度随着距离光源的距离而降低，光也会被介质吸收，因此我们应该允许吸引力随吸收程度而变化。在最简单的形式中，光强度 $I_{(r)}$ 根据平方反比定律 $I_{(r)}=I_{(s)}/r^{2}$ 变化，其中 $I_{(s)}$ 是光源处的强度。对于具有固定光吸收系数 $\gamma$ 的给定介质，光强度 $I$ 随距离 $r$ 而变化。 $I={I_{0}}^{e^{-\gamma r}}$ ，其中 $I_{0}$ 是原始光照强度。为了避免表达式中 $I_{(s)}/r^{2}$ 在 $r=0$ 处的奇点，可以使用以下高斯形式来近似平方反比定律和吸收的组合效应：

有时，我们可能需要一个以较慢速度单调递减的函数。在这种情况下，我们可以使用以下近似值：

在较短的距离内，上述两种形式基本相同。这是因为级数展开式约为 $r=0$ ：

在 $O_{(r^{3})})$ 的顺序下彼此等价。

由于萤火虫的吸引力与相邻萤火虫看到的光强度成正比，我们现在可以通过

来定义萤火虫的吸引力，其中 $\beta _{0}$ 是在 $r=0$ 处的吸引力，由于计算 $1/(1+r^{2})$ 通常比指数函数更快，因此如果需要，上述函数可以方便地替换为 $\beta =\beta _{0}/\gamma ^{r^{2}}$ 。等式 $\beta _{(r)}=\beta _{0}e^{-\gamma r^{2}}$ 定义了特征距离 $\tau =1/\gamma ^{1/2}$ ，其吸引力从 $\beta _{0}$ 到 $\beta _{0}e^{-1}$ 发生显著变化。

在实现中，吸引力函数 $\beta _{(r)}$ 的实际形式可以是任何单调递减函数，如以下广义形式：

对于固定的 $\gamma$ ，特征长度变为 $\tau =\gamma ^{-1/m}\rightarrow 1$ ，当 $m\rightarrow \infty$ ，相反，对于优化问题中给定的长度标度 $\tau$ ，参数 $\gamma$ 可以用作典型的初始值。也就是 $\gamma =\frac{1 }{\tau ^m}$ 。

（4）距离和移动

任意两个萤火虫 $i$ 和 $j$ 在 $x_{i}$ 和 $x_{j}$ 之间的距离分别是笛卡尔距离：

其中 $x_{i,k}$ 是萤火虫空间坐标的第k个分量。在二维情况下，我们得到了：

一只萤火虫被另一只更吸引人（更亮）的萤火虫吸引的运动取决于：

其中第二项是由于吸引力，而第三项是随机化， $\alpha$ 是随机化参数。rand是一个均匀分布在[0,1]中的随机数生成器。对于我们实现中的大多数情况，我们可以取 $\beta _{0}=1$ 和 $\alpha \epsilon [0,1]$ . 此外，随机化项可以很容易地扩展到正态分布 $N(0,1)$ 或其他分布。此外，如果刻度在不同维度上有显著差异，例如 $-10^{5}$ 到 $10^{5}$ 在一个维度，又比如： $-0.001$ 到 $0.01$ ，另一方面，用 $\alpha S_{k}$ 代替 $\alpha$ 是一个好主意，其中d维中的标度参数 $S_{k}(k=1,2,...,d)$ 应根据感兴趣问题的实际标度确定。

现在，参数 $\gamma$ 表示了吸引力的变化，其值对于确定收敛速度和F A算法的行为至关重要。理论上， $\gamma \epsilon [0,\infty ]$ , 但实际上， $\gamma =O_{(1)}$ 是由待优化系统的特征长度 $\tau$ 决定的。因此，在大多数应用程序中，它通常从0变化 $0.01-100$ 。

附：Yang分享的萤火虫算法MATLAB代码

fa_mincon.m

% ======================================================== %
% Files of the Matlab programs included in the book:       %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,  %
% Second Edition, Luniver Press, (2010).   www.luniver.com %
% ======================================================== %    % -------------------------------------------------------- %
% Firefly Algorithm for constrained optimization using     %
% for the design of a spring (benchmark)                   %
% by Xin-She Yang (Cambridge University) Copyright @2009   %
% -------------------------------------------------------- %function fa_mincon
% parameters [n N_iteration alpha betamin gamma]
para=[40 500 0.5 0.2 1];
format long;help fa_mincon.m
% This demo uses the Firefly Algorithm to solve the
% [Spring Design Problem as described by Cagnina et al.,
% Informatica, vol. 32, 319-326 (2008). ]% Simple bounds/limits
disp('Solve the simple spring design problem ...');
Lb=[0.05 0.25 2.0];
Ub=[2.0 1.3 15.0];% Initial random guess
u0=Lb+(Ub-Lb).*rand(size(Lb));[u,fval,NumEval]=ffa_mincon(@cost,@constraint,u0,Lb,Ub,para);% Display results
bestsolution=u
bestojb=fval
total_number_of_function_evaluations=NumEval%%% Put your own cost/objective function here --------%%%
%% Cost or Objective functionfunction z=cost(x)
z=(2+x(3))*x(1)^2*x(2);% Constrained optimization using penalty methods
% by changing f to F=f+ \sum lam_j*g^2_j*H_j(g_j)
% where H(g)=0 if g<=0 (true), =1 if g is false%%% Put your own constraints here --------------------%%%
function [g,geq]=constraint(x)
% All nonlinear inequality constraints should be here
% If no inequality constraint at all, simple use g=[];
g(1)=1-x(2)^3*x(3)/(71785*x(1)^4);
% There was a typo in Cagnina et al.'s paper,
% the factor should 71785 insteady of 7178 !
tmpf=(4*x(2)^2-x(1)*x(2))/(12566*(x(2)*x(1)^3-x(1)^4));
g(2)=tmpf+1/(5108*x(1)^2)-1;
g(3)=1-140.45*x(1)/(x(2)^2*x(3));
g(4)=x(1)+x(2)-1.5;% all nonlinear equality constraints should be here
% If no equality constraint at all, put geq=[] as follows
geq=[];%%% End of the part to be modified -------------------%%%%%% --------------------------------------------------%%%
%%% Do not modify the following codes unless you want %%%
%%% to improve its performance etc                    %%%
% -------------------------------------------------------
% ===Start of the Firefly Algorithm Implementation ======
% Inputs: fhandle => @cost (your own cost function,
%                   can be an external file  )
%     nonhandle => @constraint, all nonlinear constraints
%                   can be an external file or a function
%         Lb = lower bounds/limits
%         Ub = upper bounds/limits
%   para == optional (to control the Firefly algorithm)
% Outputs: nbest   = the best solution found so far
%          fbest   = the best objective value
%      NumEval = number of evaluations: n*MaxGeneration
% Optional:
% The alpha can be reduced (as to reduce the randomness)
% ---------------------------------------------------------% Start FA
function [nbest,fbest,NumEval]...=ffa_mincon(fhandle,nonhandle,u0, Lb, Ub, para)
% Check input parameters (otherwise set as default values)
if nargin<6, para=[20 50 0.25 0.20 1]; end
if nargin<5, Ub=[]; end
if nargin<4, Lb=[]; end
if nargin<3,
disp('Usuage: FA_mincon(@cost, @constraint,u0,Lb,Ub,para)');
end% n=number of fireflies
% MaxGeneration=number of pseudo time steps
% ------------------------------------------------
% alpha=0.25;      % Randomness 0--1 (highly random)
% betamn=0.20;     % minimum value of beta
% gamma=1;         % Absorption coefficient
% ------------------------------------------------
n=para(1);  MaxGeneration=para(2);
alpha=para(3); betamin=para(4); gamma=para(5);% Total number of function evaluations
NumEval=n*MaxGeneration;% Check if the upper bound & lower bound are the same size
if length(Lb) ~=length(Ub),disp('Simple bounds/limits are improper!');return
end% Calcualte dimension
d=length(u0);% Initial values of an array
zn=ones(n,1)*10^100;
% ------------------------------------------------
% generating the initial locations of n fireflies
[ns,Lightn]=init_ffa(n,d,Lb,Ub,u0);% Iterations or pseudo time marching
for k=1:MaxGeneration,     %%%%% start iterations% This line of reducing alpha is optionalalpha=alpha_new(alpha,MaxGeneration);% Evaluate new solutions (for all n fireflies)
for i=1:n,zn(i)=Fun(fhandle,nonhandle,ns(i,:));Lightn(i)=zn(i);
end% Ranking fireflies by their light intensity/objectives
[Lightn,Index]=sort(zn);
ns_tmp=ns;
for i=1:n,ns(i,:)=ns_tmp(Index(i),:);
end%% Find the current best
nso=ns; Lighto=Lightn;
nbest=ns(1,:); Lightbest=Lightn(1);% For output only
fbest=Lightbest;% Move all fireflies to the better locations
[ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,nbest,...Lightbest,alpha,betamin,gamma,Lb,Ub);end   %%%%% end of iterations% -------------------------------------------------------
% ----- All the subfunctions are listed here ------------
% The initial locations of n fireflies
function [ns,Lightn]=init_ffa(n,d,Lb,Ub,u0)% if there are bounds/limits,
if length(Lb)>0,for i=1:n,ns(i,:)=Lb+(Ub-Lb).*rand(1,d);end
else% generate solutions around the random guessfor i=1:n,ns(i,:)=u0+randn(1,d);end
end% initial value before function evaluations
Lightn=ones(n,1)*10^100;% Move all fireflies toward brighter ones
function [ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,...nbest,Lightbest,alpha,betamin,gamma,Lb,Ub)
% Scaling of the system
scale=abs(Ub-Lb);% Updating fireflies
for i=1:n,
% The attractiveness parameter beta=exp(-gamma*r)for j=1:n,r=sqrt(sum((ns(i,:)-ns(j,:)).^2));% Update moves
if Lightn(i)>Lighto(j), % Brighter and more attractivebeta0=1; beta=(beta0-betamin)*exp(-gamma*r.^2)+betamin;tmpf=alpha.*(rand(1,d)-0.5).*scale;ns(i,:)=ns(i,:).*(1-beta)+nso(j,:).*beta+tmpf;endend % end for jend % end for i% Check if the updated solutions/locations are within limits
[ns]=findlimits(n,ns,Lb,Ub);% This function is optional, as it is not in the original FA
% The idea to reduce randomness is to increase the convergence,
% however, if you reduce randomness too quickly, then premature
% convergence can occur. So use with care.
function alpha=alpha_new(alpha,NGen)
% alpha_n=alpha_0(1-delta)^NGen=10^(-4);
% alpha_0=0.9
delta=1-(10^(-4)/0.9)^(1/NGen);
alpha=(1-delta)*alpha;% Make sure the fireflies are within the bounds/limits
function [ns]=findlimits(n,ns,Lb,Ub)
for i=1:n,% Apply the lower boundns_tmp=ns(i,:);I=ns_tmp<Lb;ns_tmp(I)=Lb(I);% Apply the upper boundsJ=ns_tmp>Ub;ns_tmp(J)=Ub(J);% Update this new movens(i,:)=ns_tmp;
end% -----------------------------------------
% d-dimensional objective function
function z=Fun(fhandle,nonhandle,u)
% Objective
z=fhandle(u);% Apply nonlinear constraints by the penalty method
% Z=f+sum_k=1^N lam_k g_k^2 *H(g_k) where lam_k >> 1
z=z+getnonlinear(nonhandle,u);function Z=getnonlinear(nonhandle,u)
Z=0;
% Penalty constant >> 1
lam=10^15; lameq=10^15;
% Get nonlinear constraints
[g,geq]=nonhandle(u);% Apply inequality constraints as a penalty function
for k=1:length(g),Z=Z+ lam*g(k)^2*getH(g(k));
end
% Apply equality constraints (when geq=[], length->0)
for k=1:length(geq),Z=Z+lameq*geq(k)^2*geteqH(geq(k));
end% Test if inequalities hold
% H(g) which is something like an index function
function H=getH(g)
if g<=0,H=0;
elseH=1;
end% Test if equalities hold
function H=geteqH(g)
if g==0,H=0;
elseH=1;
end
%% ==== End of Firefly Algorithm implementation ======

fa_ndim.m

% ======================================================== %
% Files of the Matlab programs included in the book:       %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,  %
% Second Edition, Luniver Press, (2010).   www.luniver.com %
% ======================================================== %    % -------------------------------------------------------- %
% Firefly Algorithm for constrained optimization using     %
% for the design of a spring (benchmark)                   %
% by Xin-She Yang (Cambridge University) Copyright @2009   %
% -------------------------------------------------------- %function fa_ndim
% parameters [n N_iteration alpha betamin gamma]
para=[20 500 0.5 0.2 1];help fa_ndim.m% Simple bounds/limits for d-dimensional problems
d=15;
Lb=zeros(1,d);
Ub=2*ones(1,d);% Initial random guess
u0=Lb+(Ub-Lb).*rand(1,d);[u,fval,NumEval]=ffa_mincon(@cost,u0,Lb,Ub,para);% Display results
bestsolution=u
bestojb=fval
total_number_of_function_evaluations=NumEval%%% Put your own cost/objective function here --------%%%
%% Cost or Objective functionfunction z=cost(x)
% Exact solutions should be (1,1,...,1)
z=sum((x-1).^2);%%% End of the part to be modified -------------------%%%%%% --------------------------------------------------%%%
%%% Do not modify the following codes unless you want %%%
%%% to improve its performance etc                    %%%
% -------------------------------------------------------
% ===Start of the Firefly Algorithm Implementation ======
%         Lb = lower bounds/limits
%         Ub = upper bounds/limits
%   para == optional (to control the Firefly algorithm)
% Outputs: nbest   = the best solution found so far
%          fbest   = the best objective value
%      NumEval = number of evaluations: n*MaxGeneration
% Optional:
% The alpha can be reduced (as to reduce the randomness)
% ---------------------------------------------------------% Start FA
function [nbest,fbest,NumEval]...=ffa_mincon(fhandle,u0, Lb, Ub, para)
% Check input parameters (otherwise set as default values)
if nargin<5, para=[20 500 0.25 0.20 1]; end
if nargin<4, Ub=[]; end
if nargin<3, Lb=[]; end
if nargin<2,
disp('Usuage: FA_mincon(@cost,u0,Lb,Ub,para)');
end% n=number of fireflies
% MaxGeneration=number of pseudo time steps
% ------------------------------------------------
% alpha=0.25;      % Randomness 0--1 (highly random)
% betamn=0.20;     % minimum value of beta
% gamma=1;         % Absorption coefficient
% ------------------------------------------------
n=para(1);  MaxGeneration=para(2);
alpha=para(3); betamin=para(4); gamma=para(5);% Total number of function evaluations
NumEval=n*MaxGeneration;% Check if the upper bound & lower bound are the same size
if length(Lb) ~=length(Ub),disp('Simple bounds/limits are improper!');return
end% Calcualte dimension
d=length(u0);% Initial values of an array
zn=ones(n,1)*10^100;
% ------------------------------------------------
% generating the initial locations of n fireflies
[ns,Lightn]=init_ffa(n,d,Lb,Ub,u0);% Iterations or pseudo time marching
for k=1:MaxGeneration,     %%%%% start iterations% This line of reducing alpha is optionalalpha=alpha_new(alpha,MaxGeneration);% Evaluate new solutions (for all n fireflies)
for i=1:n,zn(i)=fhandle(ns(i,:));Lightn(i)=zn(i);
end% Ranking fireflies by their light intensity/objectives
[Lightn,Index]=sort(zn);
ns_tmp=ns;
for i=1:n,ns(i,:)=ns_tmp(Index(i),:);
end%% Find the current best
nso=ns; Lighto=Lightn;
nbest=ns(1,:); Lightbest=Lightn(1);% For output only
fbest=Lightbest;% Move all fireflies to the better locations
[ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,nbest,...Lightbest,alpha,betamin,gamma,Lb,Ub);end   %%%%% end of iterations% -------------------------------------------------------
% ----- All the subfunctions are listed here ------------
% The initial locations of n fireflies
function [ns,Lightn]=init_ffa(n,d,Lb,Ub,u0)% if there are bounds/limits,
if length(Lb)>0,for i=1:n,ns(i,:)=Lb+(Ub-Lb).*rand(1,d);end
else% generate solutions around the random guessfor i=1:n,ns(i,:)=u0+randn(1,d);end
end% initial value before function evaluations
Lightn=ones(n,1)*10^100;% Move all fireflies toward brighter ones
function [ns]=ffa_move(n,d,ns,Lightn,nso,Lighto,...nbest,Lightbest,alpha,betamin,gamma,Lb,Ub)
% Scaling of the system
scale=abs(Ub-Lb);% Updating fireflies
for i=1:n,
% The attractiveness parameter beta=exp(-gamma*r)for j=1:n,r=sqrt(sum((ns(i,:)-ns(j,:)).^2));% Update moves
if Lightn(i)>Lighto(j), % Brighter and more attractivebeta0=1; beta=(beta0-betamin)*exp(-gamma*r.^2)+betamin;tmpf=alpha.*(rand(1,d)-0.5).*scale;ns(i,:)=ns(i,:).*(1-beta)+nso(j,:).*beta+tmpf;endend % end for jend % end for i% Check if the updated solutions/locations are within limits
[ns]=findlimits(n,ns,Lb,Ub);% This function is optional, as it is not in the original FA
% The idea to reduce randomness is to increase the convergence,
% however, if you reduce randomness too quickly, then premature
% convergence can occur. So use with care.
function alpha=alpha_new(alpha,NGen)
% alpha_n=alpha_0(1-delta)^NGen=10^(-4);
% alpha_0=0.9
delta=1-(10^(-4)/0.9)^(1/NGen);
alpha=(1-delta)*alpha;% Make sure the fireflies are within the bounds/limits
function [ns]=findlimits(n,ns,Lb,Ub)
for i=1:n,% Apply the lower boundns_tmp=ns(i,:);I=ns_tmp<Lb;ns_tmp(I)=Lb(I);% Apply the upper boundsJ=ns_tmp>Ub;ns_tmp(J)=Ub(J);% Update this new movens(i,:)=ns_tmp;
end%% ==== End of Firefly Algorithm implementation ======

firefly_simple.m

% ======================================================== %
% Files of the Matlab programs included in the book:       %
% Xin-She Yang, Nature-Inspired Metaheuristic Algorithms,  %
% Second Edition, Luniver Press, (2010).   www.luniver.com %
% ======================================================== %    % =========================================================%
% Firefly Algorithm by X S Yang (Cambridge University)     %
% Usage: firefly_simple([number_of_fireflies,MaxGeneration])
%  eg:   firefly_simple([12,50]);                          %
% ======================================================== %
% This is a demo for 2D functions; for higher dimenions,   %
% you should use fa_ndim.m or fa_mincon.m                  %
% Parameters choice:
% Gamma should be linked with scales. Otherwise, the FA    %
% the efficiency will be significantly reduced because     %
% the beta term may be too small.                          %
% Similarly, alpha should also be linked with scales,      %
% the steps should not too large or too small, often       %
% steps are about 1/10 to 1/100 of the domain size.        %
% In addition, alpha should be reduced gradually           %
% using alpha=alpha_0 delta^t during eteration t.          %
% Typically, delta=0.9 to 0.99 will be a good choice.      %
% ======================================================== %function [best]=firefly_simple(instr)
% n=number of fireflies
% MaxGeneration=number of pseudo time steps
if nargin<1,   instr=[12 50];     end
n=instr(1);  MaxGeneration=instr(2);
% Show info
help firefly_simple.m
rand('state',0);  % Reset the random generator
% ------ Four peak functions ---------------------
str1='exp(-(x-4)^2-(y-4)^2)+exp(-(x+4)^2-(y-4)^2)';
str2='+2*exp(-x^2-(y+4)^2)+2*exp(-x^2-y^2)';
funstr=strcat(str1,str2);
% Converting to an inline function
f=vectorize(inline(funstr));
% range=[xmin xmax ymin ymax];
range=[-5 5 -5 5];% ------------------------------------------------
alpha=0.2;      % Randomness 0--1 (highly random)
gamma=1.0;      % Absorption coefficient
delta=0.97;      % Randomness reduction (similar to % an annealing schedule)
% ------------------------------------------------
% Grid values are used for display only
Ngrid=100;
dx=(range(2)-range(1))/Ngrid;
dy=(range(4)-range(3))/Ngrid;
[x,y]=meshgrid(range(1):dx:range(2),...range(3):dy:range(4));
z=f(x,y);
% Display the shape of the objective function
figure(1);    surfc(x,y,z);% ------------------------------------------------
% generating the initial locations of n fireflies
[xn,yn,Lightn]=init_ffa(n,range);
% Display the paths of fireflies in a figure with
% contours of the function to be optimizedfigure(2);
% Iterations or pseudo time marching
for i=1:MaxGeneration,     %%%%% start iterations
% Show the contours of the functioncontour(x,y,z,15); hold on;
% Evaluate new solutions
zn=f(xn,yn);% Ranking the fireflies by their light intensity
[Lightn,Index]=sort(zn);
xn=xn(Index); yn=yn(Index);
xo=xn;   yo=yn;    Lighto=Lightn;
% Trace the paths of all roaming  fireflies
plot(xn,yn,'.','markersize',10,'markerfacecolor','g');
% Move all fireflies to the better locations
[xn,yn]=ffa_move(xn,yn,Lightn,xo,yo,Lighto,alpha,gamma,range);
drawnow;
% Use "hold on" to show the paths of fireflieshold off;% Reduce randomness as iterations proceed
alpha=newalpha(alpha,delta);end   %%%%% end of iterations
best(:,1)=xo'; best(:,2)=yo'; best(:,3)=Lighto';% ----- All subfunctions are listed here ---------
% The initial locations of n fireflies
function [xn,yn,Lightn]=init_ffa(n,range)
xrange=range(2)-range(1);
yrange=range(4)-range(3);
xn=rand(1,n)*xrange+range(1);
yn=rand(1,n)*yrange+range(3);
Lightn=zeros(size(yn));% Move all fireflies toward brighter ones
function [xn,yn]=ffa_move(xn,yn,Lightn,xo,yo,...Lighto,alpha,gamma,range)
ni=size(yn,2); nj=size(yo,2);
for i=1:ni,
% The attractiveness parameter beta=exp(-gamma*r)for j=1:nj,
r=sqrt((xn(i)-xo(j))^2+(yn(i)-yo(j))^2);
if Lightn(i)<Lighto(j), % Brighter and more attractive
beta0=1;     beta=beta0*exp(-gamma*r.^2);
xn(i)=xn(i).*(1-beta)+xo(j).*beta+alpha.*(rand-0.5);
yn(i)=yn(i).*(1-beta)+yo(j).*beta+alpha.*(rand-0.5);
endend % end for j
end % end for i
[xn,yn]=findrange(xn,yn,range);% Reduce the randomness during iterations
function alpha=newalpha(alpha,delta)
alpha=alpha*delta;% Make sure the fireflies are within the range
function [xn,yn]=findrange(xn,yn,range)
for i=1:length(yn),if xn(i)<=range(1), xn(i)=range(1); endif xn(i)>=range(2), xn(i)=range(2); endif yn(i)<=range(3), yn(i)=range(3); endif yn(i)>=range(4), yn(i)=range(4); end
end
%  ============== end =====================================

firefly_simple.m是Yang提供的一个简单的二维空间示例。该函数的目的是寻找二元函数

在 $x\epsilon [-5,5],y\epsilon [-5,5]$ 范围内的最大值。

在计算范围内有四个极值:

参考文献：

1.Xin-She Yang Firefly Algorithms for Multimodal Optimization

2.Firefly Algorithm - File Exchange - MATLAB Central