一、简介

受到烟花在夜空中爆炸产生火花并照亮周围区域这一自然现象的启发，北大教授谭营在2010年提出了烟花算法(Fireworks Algorithm，FWA)[1]。在该算法中，烟花被看作为最优化问题的解空间中一个可行解，那么烟花爆炸产生一定数量火花的过程即为其搜索邻域的过程。
1 FWA算法框架
烟花算法的算法流程图如图1所示。

烟花算法具体包括以下几个步骤：

在可行解空间中随机产生一定数量的烟花，每个烟花代表解空间中的一个可行解。
根据优化目标函数计算每个烟花的适应度值，并据此确定烟花质量的好坏，以在不同爆炸半径下产生不同数量的火花。在烟花算法中，作者使用了两种形式的火花，分别是爆炸火花和高斯变异火花。其中爆炸火花主要负责对烟花邻近区域的搜索，适应度值好的烟花在较小的邻近区域内产生较多的火花，反之，适应度值差的烟花在较大的邻近区域内产生较少的火花。相对于爆炸火花，高斯火花的引入增强了种群的多样性。
判定是否满足终止条件。如果满足则停止搜索，否则在爆炸火花、高斯变异火花和烟花中选择一定数量的个体作为烟花进入下一代的迭代。烟花算法具有局部搜索能力和全局搜索能力自调节机制。烟花算法中每个烟花的爆炸半径和爆炸火花数是不同的，适应度值差的烟花的爆炸半径较大，使其具有更大的“探索能力”——开发性。适应度值好的烟花的爆炸半径较小，使其能够在该位置周围具有更大的“挖掘能力”——利用性。此外，高斯变异火花的引入可以进一步增加种群的多样性。

2 算子分析

2.1 爆炸算子
在可行域Ω内初始化一定数量烟花，对烟花位置的适应度值进行评估。为了差异化不同位置的烟花，一般适应度值较好的烟花能够获取更多的资源，在较小的范围内产生更多的火花，具有对于该烟花位置的强大的局部搜索能力。反之，适应度值较差的烟花只能获取相对较少的资源，在较大的范围内产生数量较少的火花，具有一定的全局搜索能力。实际上该过程中模拟不同厂家生产的不同质量的烟花，如图2所示。

（a）质量好

（b）质量差
图2 不同质量的烟花燃放效果


2.2 变异算子

2.3 选择策略
为使烟花种群中优秀的信息能够传递到下一代种群中，在产生爆炸火花和高斯变异火花后，算法会在候选者集合(包括烟花、爆炸火花和高斯变异火花)中选择一定数量的个体作为下一代的烟花。假设候选者集合为K，烟花种群大小为N。候选者集合中适应度值最小的个体会被确定性地选择到下一代作为烟花（精英策略），而对剩下的N-1个烟花的选择使用轮盘赌的方法在候选者集合中进行选择。对于候选者xi，其被选择概率的计算公式为：

二、源代码

%烟花算法进行函数优化
%fitness适应度函数，N烟花数，D变量维数，M变异火花数，Er爆炸半径，En爆炸数目
%LB,UB分别为变量上下界，T为迭代次数，a,b为爆炸数目限制因子
clear;clc
N=100;   % N烟花数
D=2;     % D变量维数
M=5;     % M变异火花数
En=6;    % En爆炸数目
Er=5;    % Er爆炸半径
a=0.3;   % a,b为爆炸数目限制因子
b=0.6;
T=500;   % T为迭代次数

%求最大值变量上下界
LB=[-5.12,-5.12];
UB=[5.12,5.12];

%随机在解空间初始化N个烟花位置
x = zeros(N,D);
for i=1:N
    x(i,:)=LB+rand(1,D).*(UB-LB);
end
%循环迭代
E_Spark=zeros(T,D,N);
Fit = zeros(1,N);
F = zeros(1,T);
for t=1:T
    %计算每个烟花适应度值
    for i=1:N
        Fit(i)=fitness(x(i,:));
    end
    [F(t),~]=min(Fit);
    Fmin=min(Fit);
    Fmax=max(Fit);
    %计算每个烟花的爆炸半径E_R和爆炸数目E_N以及产生的爆炸火花
    E_R = zeros(1,N);
    E_N = zeros(1,N);
    for i=1:N
        E_R(i)=Er*((Fit(i)-Fmin+eps)/(sum(Fit)-N*Fmin+eps));  %爆炸半径
        E_N(i)=En*((Fmax-Fit(i)+eps)/(N*Fmax-sum(Fit)+eps));  %爆炸数目
        if E_N(i)<a*En    % 爆炸数目限制
            E_N(i)=round(a*En);
        elseif E_N(i)>b*En
            E_N(i)=round(b*En);
        else
            E_N(i)=round(E_N(i));
        end
        %产生爆炸火花 E_Spark
        for j=2:(E_N(i)+1)              % 第i个烟花共产生E_N(i)个火花
            E_Spark(1,:,i)=x(i,:);      % 将第i个烟花保存为第i个火花序列中的第一个，爆炸产生的火花从序列中的第二个开始存储（即烟花为三维数组每一页的第一行）
            h=E_R(i)*(-1+2*rand(1,D));  % 位置偏移
            E_Spark(j,:,i)=x(i,:)+h;    % 第i个烟花（三维数组的i页）产生的第j（三维数组的j行）个火花
            for k=1:D   %越界检测
                if E_Spark(j,k,i)>UB(k)||E_Spark(j,k,i)<LB(k)  %第i个烟花（三维数组的i页）产生的第j个火花（三维数组的j行）的第k个变量（三维数组的k列）
                    E_Spark(j,k,i)=LB(k)+rand*(UB(k)-LB(k));   %映射规则
                end
            end
        end
    end
    %产生高斯变异火花Mut_Spark,随机选择M个烟花进行变异
    Mut=randperm(N);          % 随机产生1-N内的N个数
    for m1=1:M                % M个变异烟花
        m=Mut(m1);            % 随机选取烟花
        for n=1:E_N(m)
            e=1+sqrt(1)*randn(1,D); %高斯变异参数，方差为1，均值也为1的1*D随机矩阵
            E_Spark(n,:,m)=E_Spark(n,:,m).*e;
            for k=1:D   %越界检测
                if E_Spark(n,k,m)>UB(k)||E_Spark(n,k,m)<LB(k)  %第i个烟花（三维数组的i页）产生的第j个火花（三维数组的j行）的第k个变量（三维数组的k列）
                    E_Spark(n,k,m)=LB(k)+rand*(UB(k)-LB(k));   %映射规则
                end
            end
        end
    end
    

三、运行结果

四、备注

版本：2014a