烟花算法FWA的理论知识以及python代码实现

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：烟花算法FAW对Ackely函数的研究
获取代码方式2：
通过Performer_Cherry博客主页开通优快云会员，凭支付凭证，私信博主，可获得此代码。

获取代码方式3：
通过订阅Performer_Cherry博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：开通优快云会员，仅只能免费获得1份代码（有效期为开通日起，三天内有效）；
订阅紫极神光博客付费专栏，可免费获得2份代码（有效期为订阅日起，三天内有效）；

二、烟花算法介绍

FWA 是一种新型群体智能(Swarm Intelligence, SI)优化算法，对于待求解的优化问题 $minf(x)\in R$， $x\in \Omega$ 采用 FWA 算法求解该优化问题，具体步骤 如下：

① 种群初始化。在特定解空间中随机产生一 些烟花，每一个烟花个体 $x_i$代表解空间中的一个解， 即 $x_i\in \Omega$。

② 计算适应度值。对初始种群中的每一个烟 花个体 $xi$根据适应度函数 $f(x)$计算适应度值 $f(x_i)$，并 根据式(5)和式(6)计算每个烟花爆炸产生烟花的个 数 $Si$和爆炸半径 $Ai$。
$$S_i=c\times \frac{y_{max}-f(x_i)+\epsilon }{{\sum }_{i=1}^n(y_{max}-f(x_i))+\epsilon }(5)$$

$$A_i=d\times \frac{f(x_i)-y_{min}+\epsilon }{{\sum }_{i=1}^n(f(x_i)-y_{min})+\epsilon }(6)$$

式中，$y_{max}=max(f(x_i))(i=1,2,...,N)$是当前种群中所有 烟花个体适应度值最差烟花个体的适应度值； $y_{min}=min(f(x_i))(i=1,2,...,N)$为当前种群中所有烟花 个体适应度最优烟花个体的适应度值；$c$ 和 $d$ 为常 数，分别用来限制火花产生的总数以及表示最大的 爆炸半径，$\epsilon$ 为一个用来避免分母为零的常数。

③ 生成火花。随机选取 $z$ 个维度组成集合 $Z$， 其中，$z=rand(1,d\times rand(0,R_i))$​，且 $rand(0,R_i)$为爆炸 半径 $A_i$内产生的随机数。在集合 $Z$ 中，对于每个维 度 $k$，用式(8)和式(9)对烟花进行爆炸变异操作，通 过式(10)中的高斯变异映射规则对超出边界的火花 进行映射处理并保存在火花种群中。
$$h=A_i\times rand(-1,1)(7)$$

$$c{x}_{i}j=x_{i}j+h(8)$$

$$c{x}_{ik}=x_{ik}\times r(9)$$

式中，$A_i$为第 $i$ 个烟花的爆炸半径；$h$ 为位置偏移量； $x_{ik}$为种群中第 i 个烟花的第 k 维；$e{x}_{ik}$为第 i 个烟花 经过爆炸后产生的火花；$c{x}_{ik}$ 为 $x_{ik}$ 经过高斯变异后 的高斯变异火花，并且 $r$ 服从 $rN(-1,1)$的高斯分布。

④ 选择下一代群体。应用选择策略选择得到 下一代烟花群体，即从烟花、爆炸火花及高斯变异 火花种群中选择$N$个烟花个体形成下一代候选烟花 种群。对于候选烟花种群 $K$，选择策略如下：选择 适应度值最小$min(f(x_i))$的个体$x_k$直接为下一代烟花 种群个体，其余的$N-1$个烟花个体采取轮盘赌方式， 对于候选个体 $x_i$其被选择的概率式(10)。
$$p(x_i)=\frac{R(x_i)}{{\sum }_{j\in K}R(x_j)}(10)$$
式中，R(xi)为烟花个体 xi与其他个体的距离之和， 通过式(11)计算得出。
$$R(x_i)=\sum _{j=1}^k||x_i-x_j||(11)$$
⑤ 判断终止条件。若满足终止条件则停止迭 代，否则继续执行步骤②。

三、烟花算法的python实现

以测试函数Ackley函数的寻优为例
$$f(x)=-a\times exp(-b\sqrt{\frac{1}{d}\sum _{i=1}^d{x_i}^2})-exp(\frac{1}{d}\sum _{i=1}^d\cos c{x}_i)+a+exp(1)$$

Ackley测试函数的代码实现

def function(X):
    y1 = 0
    for i in X:
        y1 += i*i
    y2 = 0
    for i in X:
        y2 += np.cos(2*np.pi*i)
    y = -20 * np.exp(-0.2 * np.sqrt(1 / 30 * y1))-np.exp(1/30*y2)+20+np.exp(1)
    return y

初始化参数

#初始化参数
upperInit=32#最大值
lowerInit=-32#最小值
upperBound=32;lowerBound=-32
dim = 30
seednum = 5 #个数
sonnum = 50
maxEva = 1000
gaussianNum	= 5
total_fiteval_times = 0
Coef_Explosion_Amplitude = 40
Max_Sparks_Num = 40
Min_Sparks_Num = 2
Coef_Spark_Num = 50

相关数据的初始占位

SeedsMatrix = np.zeros([seednum, dim]) #种群矩阵 
fitness_best_array = np.zeros([maxEva,1+dim]) #记录适应度最好的个体
SeedsFitness = np.zeros([seednum,1]) #记录种群适应度
fitness_ave=np.zeros(maxEva) #记录平均适应

种群的初始化

for i in range(seednum):
    SeedsMatrix[i,:]=np.random.rand(1,dim)*(upperInit-lowerInit)+lowerInit

计算初始种群的适应

for i in range(seednum):
    SeedsFitness[i]=test_fwa_F.function(SeedsMatrix[i,:])

初始种群的相关数据

bestfit = np.min(SeedsFitness)#保存最优适应度(最小值)
[bestindex,n]=np.where(SeedsFitness==bestfit)
fitness_best_array[0,0]=bestfit # 把初始的最好适应度放在(0,0)位置
fitness_best_array[0,1:]=SeedsMatrix[bestindex,:] 

FWA算法的迭代寻优

for i in range(maxEva):
    #计算数量Si，计算每个烟花的火花数
    sonsnum_array = test_fwa_S.sonsnum_cal(SeedsFitness,Max_Sparks_Num,Min_Sparks_Num,Coef_Spark_Num,seednum)
    #计算火星振幅Ai
    scope_array = test_fwa_S.scope_cal(SeedsFitness,Coef_Explosion_Amplitude,seednum);
    #根据每个烟花的火花数和爆炸幅度产生火花
    [SonsMatrix,SonsFitness] = test_fwa_S.sons_generate(sonsnum_array, scope_array, SeedsMatrix,seednum,dim,upperBound,lowerBound)
    #产生高斯火花
    [SeedsMatrixGauss,SeedsFitGaussMutation] = test_fwa_S.seedGaussMutation(SeedsMatrix, gaussianNum,dim,seednum,upperBound,lowerBound)
    #合并种群及适应度
    AllMatrix=np.vstack((SeedsMatrix,SonsMatrix,SeedsMatrixGauss))#vstack在列上合并（要求列数一样）
    AllFitness = np.vstack((SeedsFitness,SonsFitness,SeedsFitGaussMutation));
    #保存最优
    bestfit=np.min(AllFitness)
    [bestindex,_]=np.where(AllFitness==bestfit)
    fitness_best_array[i,0] = bestfit    # 放最好适应度
    # 放最好个体
    fitness_best_array[i,1:]=AllMatrix[bestindex[0],:]
    #选择 
    [SeedsMatrix,SeedsFitnessCurrentIte] = test_fwa_S.selectNextIterationOnEntropy(AllMatrix,AllFitness,dim,seednum)
    fitness_ave[i]=np.mean(AllFitness)
    SeedsFitness = SeedsFitnessCurrentIte
    print("第",i,"次训练","最好的适应度为：",bestfit,"平均适应度为：",fitness_ave[i])

烟花个体进行爆炸（计算数量Si，计算每个烟花的火花数）

#计算数量Si
def sonsnum_cal(fitness_array,Max_Sparks_Num,Min_Sparks_Num,Coef_Spark_Num,seednum):
    fitness_max = max(fitness_array)
    fitness_sub_max = abs(fitness_max - fitness_array)
    fitness_sub_max_sum = sum(fitness_sub_max)
    sonsnum_array = np.zeros([seednum,1])
    for i in range(seednum):
        sonnum_temp = (fitness_sub_max[i] + eps) / (fitness_sub_max_sum + eps)
        sonnum_temp = np.rint( sonnum_temp * Coef_Spark_Num)#取整
        if sonnum_temp > Max_Sparks_Num:
            sonnum_temp = Max_Sparks_Num
        elif sonnum_temp < Min_Sparks_Num:
            sonnum_temp = Min_Sparks_Num
        sonsnum_array[i] = sonnum_temp
    return sonsnum_array

计算火星振幅Ai

#计算火星振幅Ai
def scope_cal(fitness_array, Coef_Explosion_Amplitude,seednum):
    fitness_best = min(fitness_array) # 最好适应度
    fitness_sub_best = abs(fitness_best - fitness_array) # 适应度与最好适应度的差距
    fitness_sub_best_sum = sum(fitness_sub_best) # 适应度差距的求和
    scope_array = np.zeros([seednum,1]) # [seednum,1]
    for i in range(seednum):
        scope_array[i] = Coef_Explosion_Amplitude * (fitness_sub_best[i] + eps) / (fitness_sub_best_sum + eps)
    return scope_array

根据每个烟花的火花数和爆炸幅度产生火花

def sons_generate(sonsnum_array, scope_array, seeds_matrix,seednum,dim,upperBound,lowerBound):
    fiteval_time = int(sum(sonsnum_array)) # 所有烟花产生的火花总数
    sons_matrix = np.zeros([fiteval_time, dim]) # [fiteval_time, dim]，火花矩阵
    sons_fitness_array = np.zeros([fiteval_time,1]) # 适应度核矩阵占位
    sons_index = 0
    #新火花
    for i in range(seednum):
    	 # i 来选择每次用来产生火花的烟花
        for j in range(int(sonsnum_array[i])):
            # j用来确定每次烟花产生火花的个数
            seed_position = np.copy(seeds_matrix[i,:])
            allowed_dim_array = np.ones([1,dim])    #用于记录某一维度是否。。
            dimens_select = np.ceil(np.random.rand()*dim) #从D维中选z维,z=dimens_select
            offset = (np.random.rand()*2-1) * scope_array[i] #计算Ai*rand(-1,1)
            for k in range(int(dimens_select)):
                rand_dimen = int(np.ceil(np.random.rand()*dim))-1;
                while allowed_dim_array[0,int(rand_dimen)]==0:
                    rand_dimen = int(np.ceil(np.random.rand()*dim))-1;

                allowed_dim_array[0,rand_dimen]=0;
                seed_position[rand_dimen] = seed_position[rand_dimen] + offset; # 选择一个维度，进行随机
                if seed_position[rand_dimen] > upperBound:
                    seed_position[rand_dimen]=upperBound
                elif seed_position[rand_dimen] < lowerBound:
                    seed_position[rand_dimen]=lowerBound
            sons_matrix[sons_index,:] = seed_position
            sons_index = sons_index + 1
    #新火花适应度
    for i in range(fiteval_time):
        sons_fitness_array[i] = test_fwa_F.function(sons_matrix[i,:])
    return sons_matrix, sons_fitness_array

高斯火花，轮盘对赌选择

#高斯火花，爆炸变异操作
def seedGaussMutation(seedMatrix, gaussianNum,dim,seednum,upperBound,lowerBound):
    gaussian_number = gaussianNum
    # 随机一个大小一样的高斯矩阵
    seed_gaussian_matrix = np.zeros([gaussian_number, dim])
    # 适应度占位
    fit_seed_gaussian_array = np.zeros([gaussian_number,1])

    # i选择个体，j选择维度，
    for i in range(gaussian_number):
        #在火花中选几个来产生特别的
        seed_select = int(np.ceil(np.random.rand() * seednum))-1#选中的火花标号
        rand_dimens = int(np.ceil(np.random.rand() * dim))#选维度个数
        seed_gaussian_matrix[i,:]= np.array(seedMatrix[seed_select,:])
        coef_gaussian = np.random.normal(1,1,1)#randGauss(1,1)
        allow_dimension_array = np.ones([1,dim])

        for j in range(rand_dimens):
            dim_mutation = int(np.ceil(np.random.rand() * dim))-1#选维度
            while allow_dimension_array[0,dim_mutation]==0:
                dim_mutation = int(np.ceil(np.random.rand() * dim))-1

            allow_dimension_array[0,dim_mutation] = 0
            seed_gaussian_matrix[i, dim_mutation] = seed_gaussian_matrix[i, dim_mutation] * coef_gaussian
        fit_seed_gaussian_array[i] = test_fwa_F.function(seed_gaussian_matrix[i,:])
    return seed_gaussian_matrix,fit_seed_gaussian_array

轮盘赌选择（简单的，不考虑距离）

def selectNextIterationOnEntropy(seedMatrix,seedFitness,dim,seednum):
    seedNextMatrix = np.zeros([seednum,dim])#输出矩阵
    seedNextMatrixFitness = np.zeros([seednum,1])
    #最佳直接保留
    bestfit = np.min(seedFitness)
    [bestindex,n]=np.where(seedFitness==bestfit)
    seedNextMatrix[-1,:]=seedMatrix[bestindex[0],:]
    #最小值留下概率大
    fitness=(max(seedFitness)-seedFitness)/(max(seedFitness)-min(seedFitness))
    #轮盘赌选择
    sumFits = sum(fitness)
    rndPoint = np.random.uniform(0, sumFits, seednum-1)
    rndPoint = np.sort(rndPoint)
    
    accumulator = 0.0
    k=0
    for ind, val in enumerate(fitness):#ind下标，val值
        accumulator += val
        if accumulator >= rndPoint[k]:
            seedNextMatrix[k,:]=seedMatrix[ind,:]
            k=k+1
            if k>=seednum-1:
                break
    for i in range(seednum):
        seedNextMatrixFitness[i]=test_fwa_F.function(seedNextMatrix[i,:])
    return seedNextMatrix,seedNextMatrixFitness

图像展示训练结果

plt.plot(fitness_best_array[:,0])
plt.show()

训练结果的展示题

代码的付费链接

本文完整代码的付费下载链接