萤火虫算法

1.1 萤火虫的生物现象

在热带的夏夜中，萤火虫会聚集在一起产生短暂而有节奏的光，不同种类的萤火虫的闪光模式往往是不同的。这种闪光的基本功能有三个：①是吸引交配的异性伙伴 ②则是吸引潜在的猎物 ③保护性的警告机制，告诉捕猎者其毒性或苦味

在光源远处r距离的光强服从平方反比定律: 光强I随着距离r的增加而减小，I∝1/r^2
并且，空气也会吸收光，随着距离增大，光会变得越来越弱。

这两个综合的因素使大多数萤火虫视觉到一个有限的距离，通常在晚上几百米，但这也足以让萤火虫交流。

闪烁光的表述可以与要优化的目标函数相关联，这使得制定新的优化算法成为可能。

2.1 设计萤火虫算法

我们可以由萤火虫的一些闪烁特性，从而开发萤火虫启发的算法。为了简单地描述萤火虫算法(FA)，现在使用以下三个理想的规则：
①所有的萤火虫雌雄同体的，所以一只萤火虫有可能会被其他任意的萤火虫所吸引（只要这只萤火虫光强更大），无论它们的性别。
②吸引力与它们的亮度成正比，因此对于任何两只闪烁的萤火虫，越不亮的萤火虫就会向越亮的萤火虫移动。吸引力与亮度成正比，它们都随着距离的增加而减小。如果没有比某只萤火虫更亮的萤火虫，它会随机移动；
③萤火虫的亮度受到目标函数的影响。也就是说，萤火虫本身的亮度是由目标函数决定的，而萤火虫对别的萤火虫的吸引是由别的萤火虫看到他的光强决定的，这个光强会随着距离增大而减少。

2.2 伪代码 Pseudo code of the firefly algorithm (FA)

基于以上三条规则
萤火虫算法(FA)的基本步骤可以总结为如图所示的伪代码。

3.1 光的强度和吸引力

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

萤火虫在特定位置x的亮度I可以选择为I(x)∝f(x)。然而，吸引力是相对的，它应该由其他萤火虫判断，并且随着萤火虫i和萤火虫j之间的距离的变化而变化。距离越大，光强被空气吸收的越多，吸引力也因此而减少。

①平方反比定律和吸收的联合效应可以近似为以下高斯形式：

其中 i₀是初始的光强，γ为固定光吸收系数，r为两个萤火虫之间的距离。

②由于萤火虫的吸引力与相邻萤火虫看到的光强度成正比，我们现在可以定义萤火虫的吸引力β为

其中，β₀是r=0的吸引力，β₀一般取为1。任意两个萤火虫i和j在xi和xj之间的距离分别为笛卡尔距离。

③并且为了增大跳出局部最优的能力，还加入了随机性因素

其中α为随机化参数，α ∈ [0, 1]， 是由高斯分布或均匀分布中抽取的随机数向量，在最简单的二维空间下，使用rand − 1/2进行表示，即此向量的每个分量是[-0.5,0.5]之间的随机数。

萤火虫i被另一个更有吸引力（更亮）的萤火虫j所吸引，i向j移动的过程综合上述三条，表示为以下公式：
其中，坐标xj-坐标xi意味着i的移动方向为j。

值得指出的是，这是一种偏向于向着更明亮的萤火虫的随机行走。如果β₀=0，它将变成一个简单的随机游走。此外，随机化项可以很容易地扩展到其他分布，如L‘evy flight。

参数γ现在表征了吸引力的变化，它的值在决定收敛速度和FA算法的行为方面至关重要。理论上，γ∈[0，∞)，但在实践中，对于大多数应用程序，它通常从0.1到10不等。

需要注意的是：在算法的实际代码中，萤火虫之间的吸引并不是只有最亮那个才能对别的萤火虫产生吸引，而是对于萤火虫i而言，其他所有的萤火虫都会比较自己与萤火虫i的光强，如果i更亮，那么别的萤火虫再去计算他们之间的吸引力，然后根据吸引力向着萤火虫i进行移动，也就是说，只有一个萤火虫的光强大于别的任何萤火虫，就会对别的萤火虫产生吸引，微观上来看是最亮那个吸引力最强而已。

4.1 参数的理解

当γ→0和γ→∞时分别有两种重要的极限情况。

①对于γ→0来说：吸引力是恒定的β=β0这相当于说在一个理想化的天空中，光的强度不会降低。因此，在该领域的任何地方都可以看到一只闪烁的萤火虫。因此，可以很容易地达到一个单一的（通常是全局的）最佳状态。如果我们去掉FA伪代码中j的内环（inner loop），并将xj替换为当前的全局最佳g∗，那么萤火虫算法就成为粒子群优化算法(PSO)的特殊情况。因此，这种特殊情况的效率与PSO相同。

② 极限情况γ→∞ 这意味着在其他萤火虫的视线中，吸引力几乎为零。这相当于萤火虫在一个非常厚的有雾的区域随机漫游的情况。看不到其他的萤火虫，每只萤火虫都以一种完全随机的方式漫游。因此，这与完全随机搜索方法相对应。

由于萤火虫算法通常处于这两个极端之间，因此可以调整参数γ和α，使其性能优于随机搜索和PSO。实际上，FA可以同时有效地找到全局最优和局部最优。

5.1 改进措施

基础的萤火虫算法是非常有效的，但我们可以预料到，随着迭代次数增加，生成的解逐渐接近最优解，但是解决方案仍在变化并且可能由于随机项的突然增大反而偏离最优解。这种情况下，我们可以通过随着迭代次数增加，随机性减少这种方法来提高解决方案的质量。

有两种方案可供选择：
①

其中t_max为最大迭代次数。α₀为初始随机化参数，α_∞为最终值。
很明显，算法迭代的开始阶段的α是接近于α₀的，之后便逐渐减少，直到趋近于α_∞

②

这是类似于模拟退火算法的方法，θ∈(0,1]是随机性还原常数，也即是一个0~1的随机数。
很明显，算法迭代的开始阶段的α是接近于α₀的，之后便逐渐减少，直到趋近于0。

6.1 实际应用及代码

为了证明全局最优和局部最优可以同时找到，我们现在使用下面的四峰函数，在实现过程中，参数的值分别为α=0.2、γ=1和β₀=1，delta=0.97并且x、y的取值范围都是[-5,5]

由Yang xin she教授给出的Matlab代码：

% ======================================================== % 
% 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 %%不知道这个nargin是干嘛的 也不知道是这个if是干嘛的，应该有无都无所谓 确定的是instr=[12 50]中12是萤火虫个数，50是迭代次数
    instr=[12 50];
end
n=instr(1);
MaxGeneration=instr(2);
% Show info
help firefly_simple.m %%没啥用
rand('state',0);  % Reset the random generator 产生一个随机数，"state"的意思有下标地创建随机数，如果换一个下标比如从0换成1则重新随机，如果不换下标则随机数无论怎么rand都还是这个不会变化
% ------ 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));%%inline把串转化为函数  vectorize将函数的标量运算转换为函数运算：如矩阵相乘*转换为.*点乘 ^2转换为.^2 前者是矩阵乘矩阵 后者是对矩阵中每个元素进行平方
% 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));%%生成一个三维图像 并对x、y轴-5~5进行分格 间隔为0.1 再加上零点所以一共101个点
z=f(x,y);%z为计算此四峰函数在不同分隔点上对应的函数值
% Display the shape of the objective function 作出此四峰函数的图像
figure(1);
surfc(x,y,z);

% ------------------------------------------------
% generating the initial locations of n fireflies 对12个萤火虫的初始位置进行初始化
[xn,yn,Lightn]=init_ffa(n,range);%%使用init_ffa函数进行初始化并返回值，xn为萤火虫的x轴初始坐标、yn为萤火虫的y轴初始坐标，lightn为
% Display the paths of fireflies in a figure with
% contours of the function to be optimized
 figure(2);
% Iterations or pseudo time marching
for i=1:MaxGeneration,     %%%%% start iterations 开始迭代
% Show the contours of the function
 contour(x,y,z,15); hold on;%创建等高线图 15条等高线，函数x轴为x矩阵，y轴为y矩阵，函数值大小为z矩阵，颜色由蓝变浅，越浅说明函数值越大
% Evaluate new solutions 计算xn，yn坐标下的函数值zn
zn=f(xn,yn);

% Ranking the fireflies by their light intensity 比较不同萤火虫的亮度
[Lightn,Index]=sort(zn);%lightn中存储的是zn进行升序排序之后的矩阵，Index存储的是进行升序排序之后lightn中每个元素对应于排序前zn中每个元素的位置如A={19,23,12} 排序后lightn={12,19,23},Index={3,1,2}
xn=xn(Index); yn=yn(Index);%排序之后xn=xn(Index); yn=yn(Index);意味着xn中存储的是对应函数值升序排序的x坐标，比如第一个xn(1)存的就是函数值zn最小的那个萤火虫的x坐标，y同理
xo=xn;   yo=yn;    Lighto=Lightn;%使用新的xo、yo、lightno存储，作用是作为参数传入ffa_move来进行萤火虫每次迭代的吸引与随机移动工作，移动之后的新坐标再传回xn、yn
% 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);%通过ffa_move函数进行吸引与随机移动
drawnow;
% Use "hold on" to show the paths of fireflies
    hold off;
    
% Reduce randomness as iterations proceed 随着迭代次数增加减少随机性alpha
alpha=newalpha(alpha,delta);
    
end   %%%%% end of iterations
%迭代完成之后，获取最终的答案
%best中存的就是萤火虫（x,y）处对应的函数值第一列存的是x，第二列存的y，第三列存的函数值
best(:,1)=xo';
best(:,2)=yo';
best(:,3)=Lighto';

end %%end for firfly_simple

% ----- 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));
end

% 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);
end
    end % end for j
end % end for i
[xn,yn]=findrange(xn,yn,range);
end

% Reduce the randomness during iterations
function alpha=newalpha(alpha,delta)
alpha=alpha*delta;
end 

% 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); end
   if xn(i)>=range(2), xn(i)=range(2); end
   if yn(i)<=range(3), yn(i)=range(3); end
   if yn(i)>=range(4), yn(i)=range(4); end
end
end
%  ============== end =====================================

执行结果：

7.1 带约束的优化问题——Spring Design Problem

一个带约束实际优化问题 弹簧设计问题：拉开弹簧和压缩弹簧的设计是一个著名的基准优化问题。其主要目的是根据偏转、应力、喘振频率和几何形状的限制来最小化重量。它涉及三个设计变量：线圈直径x₁，线圈直径x₂和线圈的数量/长度x₃。这个问题可以概括为
受以下约束条件的约束

MATLAB代码：from Dr.Yang

% ======================================================== % 
% 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 随机猜一个解（x1，x2，x3）
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 function
 function 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=[];
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)
%nargin为关键字：代表本函数输入的参数个数
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);%ns为初始化后的各个萤火虫坐标

% Iterations or pseudo time marching
for k=1:MaxGeneration,     %%%%% start iterations

% This line of reducing alpha is optional
 alpha=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);% lightn是排序后的函数值
ns_tmp=ns;
for i=1:n,
 ns(i,:)=ns_tmp(Index(i),:);%对ns也按照lightn排序后的顺序重新对应索引
end

%% Find the current best
nso=ns; Lighto=Lightn;
nbest=ns(1,:); Lightbest=Lightn(1);%nbest为最佳的萤火虫的坐标 Lightbest为最佳的萤火虫的适应度函数值

% 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);7

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)%n萤火虫个数，d维度，u0随机猜的一个解，ns初始解空间，lightn是排序后的函数值
  % 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 guess 无上下限就在随机解附近去生成所以随机解
   for 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 attractive
   beta0=1; beta=(beta0-betamin)*exp(-gamma*r.^2)+betamin;%这里beta为什么是需要计算的呢，而且与firefly_simple不同？
   tmpf=alpha.*(rand(1,d)-0.5).*scale;
   ns(i,:)=ns(i,:).*(1-beta)+nso(j,:).*beta+tmpf;%移动
      end
   end % end for j

end % 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 bound
  ns_tmp=ns(i,:);
  I=ns_tmp<Lb;
  ns_tmp(I)=Lb(I);

  % Apply the upper bounds
  J=ns_tmp>Ub;
  ns_tmp(J)=Ub(J);
  % Update this new move
  ns(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)
%%这里有约束的情况其实就是把约束带上，直接求整体的函数值，因为约束前面的系数极大10^15所以如果要整体最小的话，那么必然要让约束趋近于0才可以，那么当整体最小的时候，约束也必然被满足，这就是惩罚法
Z=0;
% Penalty constant >> 1
lam=10^15; lameq=10^15; 
% Get nonlinear constraints
[g,geq]=nonhandle(u)%g是非等式约束求出的约束值 geq是什么等式约束求出的约束值

% Apply inequality constraints as a penalty function
%%加上不等式约束条件
for k=1:length(g),
    Z=Z+ lam*g(k)^2*getH(g(k))%lam应该是对应的惩罚法中的ui（u），lameq对应vi（v）
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;
else
    H=1;
end

% Test if equalities hold
function H=geteqH(g)
if g==0,
    H=0;
else
    H=1;
end
%% ==== End of Firefly Algorithm implementation ======