花授粉优化算法

1. 简介

    花授粉优化算法(Flower Pollination Algorithm, FPA)是2012年由英国学者杨新社提出的一种新型的元启发式群智能优化算法。由于该算法简单、易实现,而被广泛应用。花授粉算法与粒子群算法类似,也存在易陷入局部最优且后期收敛速度慢等缺点。为此,国内外学者对其缺点进行诸多改进。

2. 花授粉算法

    在自然界中,花朵授粉过程是非常复杂的。在设计算法时,很难通过模拟花朵授粉过程来模拟授粉过程的每个细节。此外,真实模拟花授粉过程将使算法特别复杂。基于简化原则,我们从数学上对花授粉过程的主要特征进行如下描述:
(1) 花朵异花授粉是带有花粉的传播者通过Levy飞行进行全局授粉。
(2) 花朵自花授粉是局部授粉过程。
(3) 转换概率 p p p控制着全局授粉与局部授粉之间的平衡。

2.1 可行解的表示

    花粉配子是优化问题的可行解,数学描述如下:
X = ( x 1 , 1 x 1 , 2 ⋯ x 1 , d x 2 , 1 x 2 , 2 ⋯ x 2 , d ⋮ ⋮ ⋮ ⋮ x n , 1 x n , 2 ⋯ x n , d ) \begin{gathered} X= \begin{pmatrix} x_{1,1} & x_{1,2} & \cdots & x_{1,d} \\ x_{2,1} & x_{2,2} & \cdots & x_{2,d} \\ \vdots & \vdots & \vdots & \vdots \\ x_{n,1} & x_{n,2} & \cdots & x_{n,d} \\ \end{pmatrix} \end{gathered} X=x1,1x2,1xn,1x1,2x2,2xn,2x1,dx2,dxn,d
其中, n n n 为花粉配子数目; d d d 为求解问题的维度。

function population = Initial(num_pop, dim_pop, upper, lower)
%{
      函数功能: 初始化种群
       num_pop: 可行解数目,即种群数目
       dim_pop: 问题维度
 upper & lower: 问题的取值范围,即上下界
%}
    population = zeros(num_pop, dim_pop);
    num_boundary = size(upper,2);
    if num_boundary==1
        for i=1:num_pop
            population(i,:) = lower + rand(1, dim_pop)*(upper - lower);
        end
    else
        for i=1:dim_pop
            population(:, i) = lower(i) + rand(num_pop, 1)*(upper(i) - lower(i));
        end
    end
end

2.2 异花授粉公式

    在每一代种群中,FPA以固定的概率执行异花授粉操作,其数学描述如下:
X i t + 1 = X i t + L ( X i t − g ∗ t ) X_i^{t+1} = X_i^{t} + L(X_i^{t} - g_*^t ) Xit+1=Xit+L(Xitgt)
其中, X i t X_i^{t} Xit分别表示第 t t t代第 i i i个解; g ∗ t g_*^t gt为第 t t t代最优解; L L L为步长。
     在Mantegna算法中,步长为
s = u ∣ v ∣ 1 β s= \frac{u}{|v|^{\frac{1}{\beta}}} s=vβ1u
其中 u u u v v v均服从标准正态分布,即
u ∼ N ( 0 , σ u 2 ) v ∼ N ( 0 , σ v 2 ) u \sim N(0,\sigma_u^2) \\ v \sim N(0,\sigma_v^2) uN(0,σu2)vN(0,σv2)
其中,
σ u = { Γ ( 1 + β ) sin ⁡ ( π β / 2 ) Γ ( 1 + β ) / β 2 ( β − 1 ) / 2 } 1 β σ v = 1 \sigma_u = \{\frac{\Gamma(1+\beta)\sin(\pi\beta/2)}{\Gamma(1+\beta)/\beta2^{(\beta- 1)/2}}\}^{\frac{1}{\beta}} \\ \sigma_v = 1 σu={Γ(1+β)/β2(β1)/2Γ(1+β)sin(πβ/2)}β1σv=1

function L = Levy(d)
%{
   Levy: 莱维飞行
      d: 问题维数
%}
  beta = 3/2;  % 默认参数
  sigma = (gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
  u = normrnd(0, sigma, [1 d]);
  v = normrnd(0, 1, [1 d]);
  L = 0.01*u./abs(v).^(1/beta);
end

在这里插入图片描述

2.3 自花授粉公式

    自花授粉算法的设计思想是模拟自然界中同一品种花之间近距离授粉,其数学描述如下:
X i t + 1 = X i t + ϵ ( X j t − X k t ) X_i^{t+1} = X_i^{t} + \epsilon(X_j^{t} - X_k^{t} ) Xit+1=Xit+ϵ(XjtXkt)
其中, ϵ ∈ [ 0 , 1 ] \epsilon \in [0,1] ϵ[0,1] 表示随机数; X j t X_j^{t} Xjt X k t X_k^{t} Xkt分别表示第 t t t代种群中第 j j j和第 k k k个解。

3. 仿真实验

3.1 测试函数

    所有的测试函数维数为2、5和10;种群规模为50;最大迭代次数为1000.

函数名称函数表达式解空间最小值
Griewank f 1 = 1 4000 ∑ i = 1 d x i 2 − ∏ i = 1 n cos ⁡ ( x i i ) + 1 f_1 = \frac{1}{4000}\sum_{i=1}^dx_i^2-\prod_{i=1}^n\cos(\frac{x_i}{\sqrt{i}})+1 f1=40001i=1dxi2i=1ncos(i xi)+1 [ − 600 , 600 ] [-600,600] [600,600]0
Alpine f 2 = ∑ i = 1 n ∣ x i sin ⁡ ( x i ) + 0.1 x i ∣ f_2 = \sum_{i=1}^n|x_i\sin(x_i)+0.1x_i| f2=i=1nxisin(xi)+0.1xi [ − 10 , 10 ] [-10,10] [10,10]0
Ackley f 3 = − 20 e x p ( − 1 5 1 d ∑ i = 1 d x i 2 ) − e x p ( 1 d ∑ i = 1 d cos ⁡ ( 2 π x i ) ) + 20 + e f_3 = -20exp(-\frac{1}{5}\sqrt{\frac{1}{d}\sum_{i=1}^dx_i^2})-exp(\frac{1}{d}\sum_{i=1}^d\cos(2\pi x_i))+20+e f3=20exp(51d1i=1dxi2 )exp(d1i=1dcos(2πxi))+20+e [ − 32 , 32 ] [-32,32] [32,32]0
Rastrigin f 4 = ∑ i = 1 d [ x i 2 − 10 cos ⁡ ( 2 π x i ) + 10 ] f_4 = \sum_{i=1}^d{[x_i^2-10\cos(2\pi{x_i})+10]} f4=i=1d[xi210cos(2πxi)+10] [ − 5.12 , 5.12 ] [-5.12,5.12] [5.12,5.12]0
Sphere f 5 = ∑ i = 1 n x i 2 f_5 = \sum_{i=1}^nx_i^2 f5=i=1nxi2 [ − 100 , 100 ] [-100,100] [100,100]0
function [fobj, bound] = Optimizer(select)
    %% 目标函数
    switch select
        case 1
            fobj = @ Sphere;  % 效果很好  [-100, 100]
            bound = [-100, 100];
        case 2
            fobj = @ Ackley;  % 效果比较好 [-32, 32]
            bound = [-32.768, 32.768];
        case 3
            fobj = @ Rastrigin;  % 效果好 [-5.12, 5.12]
            bound = [-5.12, 5.12];
        case 4
            fobj = @ Alpine;  % 效果好 [-10,10]
            bound = [-10, 10];
        case 5
            fobj = @ Squares; % 效果好 [-10,10]
            bound = [-10, 10];
        case 6
            fobj = @ Griewank; % 效果好 [-600, 600]
            bound = [-600, 600];
        case 7
            fobj = @ Powell;  % 效果好  [-4, 5]
            bound = [-4, 5];
        case 8
            fobj = @ Zakharov; % 效果好 [-5, 10]
            bound = [-5, 10];
        case 9
            fobj = @ Sumpow;  % 效果好 [-1, 1]
            bound = [-1, 1];
        case 10
            fobj = @ Rothyp;  % 效果好 [-65, 65]
            bound = [-65, 65];
        case 11
            fobj = @ Schaffer; % [-100, 100]
            bound = [-100, 100];
        case 12
            fobj = @ schw; % [-100, 100]
            bound = [-500, 500];
        case 13
            fobj = @ shekel; % [-100, 100]
            bound = [0, 10];
        case 14
            fobj = @ trid; % [-100, 100]
            bound = [-100, 100];
    end

    function [y] = Sphere(xx)
        d = length(xx);
        sum = 0;
        for ii = 1:d
	        xi = xx(ii);
	        sum = sum + xi^2;
        end
        y = sum;
    end

    function [y] = Ackley(xx)
        d = length(xx);
        sum1 = 0;
        sum2 = 0;
        for ii = 1:d
	        xi = xx(ii);
	        sum1 = sum1 + xi^2;
	        sum2 = sum2 + cos(2*pi*xi);
        end
        term1 = -20*exp(-0.2*sqrt(sum1/d));
        term2 = -exp(sum2/d);
        
        y = term1 + term2 + 20 + exp(1);
    end

    function [y] = Rastrigin(xx)
        d = length(xx);
        sum = 0;
        for ii = 1:d
            xi = xx(ii);
            sum = sum + (xi^2 - 10*cos(2*pi*xi));
        end
    
        y = 10*d + sum;
    end

    function [y] = Alpine(xx)
        d = length(xx);
        sum = 0; 
        for i=1:d
           sum = sum + abs(xx(i)*sin(xx(i))+0.1*xx(i)); 
        end
        y = sum;
    end

    function [y] = Squares(xx)
        d = length(xx);
        sum = 0;
        for i = 1:d
            xi = xx(i);
            sum = sum + i*xi^2;
        end
        y = sum;
    end

    function [y] = Griewank(xx)
        d = length(xx);
        sum = 0;
        prod = 1;
        
        for ii = 1:d
	        xi = xx(ii);
	        sum = sum + xi^2/4000;
	        prod = prod * cos(xi/sqrt(ii));
        end
        
        y = sum - prod + 1;
    end

    function [y] = Powell(xx)
        d = length(xx);
        sum = 0;
        
        for ii = 1:(d/4)
	        term1 = (xx(4*ii-3) + 10*xx(4*ii-2))^2;
	        term2 = 5 * (xx(4*ii-1) - xx(4*ii))^2;
	        term3 = (xx(4*ii-2) - 2*xx(4*ii-1))^4;
	        term4 = 10 * (xx(4*ii-3) - xx(4*ii))^4;
	        sum = sum + term1 + term2 + term3 + term4;
        end
        
        y = sum;
    end

    function [y] = Zakharov(xx)
        d = length(xx);
        sum1 = 0;
        sum2 = 0;

        for ii = 1:d
            xi = xx(ii);
            sum1 = sum1 + xi^2;
            sum2 = sum2 + 0.5*ii*xi;
        end
    
        y = sum1 + sum2^2 + sum2^4;
    end

    function [y] = Sumpow(xx)
        d = length(xx);
        sum = 0;
        for ii = 1:d
            xi = xx(ii);
            new = (abs(xi))^(ii+1);
            sum = sum + new;
        end
    
        y = sum;
    end

    function [y] = Rothyp(xx)
        d = length(xx);
        outer = 0;
        
        for ii = 1:d
            inner = 0;
            for jj = 1:ii
                xj = xx(jj);
                inner = inner + xj^2;
            end
            outer = outer + inner;
        end
        
        y = outer;
    end

    function [y] = Schaffer(xx)
        x1 = xx(1);
        x2 = xx(2);
    
        fact1 = (sin(x1^2-x2^2))^2 - 0.5;
        fact2 = (1 + 0.001*(x1^2+x2^2))^2;
    
        y = 0.5 + fact1/fact2;
    end
    
    function y = schw(x)
        n = 2;
        s = sum(-x.*sin(sqrt(abs(x))));
        y = 418.9829*n+s;
    end

    function y = shekel(x)  % 4变量
        m = 10;
        a = ones(10,4);
        a(1,:) = 4.0*a(1,:);
        a(2,:) = 1.0*a(2,:);
        a(3,:) = 8.0*a(3,:);
        a(4,:) = 6.0*a(4,:);
        for j = 1:2
           a(5,2*j-1) = 3.0; a(5,2*j) = 7.0; 
           a(6,2*j-1) = 2.0; a(6,2*j) = 9.0; 
           a(7,j)     = 5.0; a(7,j+2) = 3.0;
           a(8,2*j-1) = 8.0; a(8,2*j) = 1.0;
           a(9,2*j-1) = 6.0; a(9,2*j) = 2.0;
           a(10,2*j-1)= 7.0; a(10,2*j)= 3.6;
        end
        c(1) = 0.1; c(2) = 0.2; c(3) = 0.2; c(4) = 0.4; c(5) = 0.4;
        c(6) = 0.6; c(7) = 0.3; c(8) = 0.7; c(9) = 0.5; c(10)= 0.5;
        s = 0;
        for j = 1:m
           p = 0;
           for i = 1:4
              p = p+(x(i)-a(j,i))^2;
           end
           s = s+1/(p+c(j));
        end
        y = -s;
    end

    function  y = trid(x)
        n = 10;
        s1 = 0;
        s2 = 0;
        for j = 1:n
            s1 = s1+(x(j)-1)^2;    
        end
        for j = 2:n
            s2 = s2+x(j)*x(j-1);    
        end
        y = s1-s2;
    end
end

3.2 算法测试

    本小节对5种测试函数进行仿真实现,在相同环境下进行30次独立实验。其中,std为解的标准差、best为最优解、mean为解的平均值、worest为最差解。最优解和平均值可以反映算法性能的上限和搜索能力。最差解和标准差分别表示算法性能的下限和算法的鲁棒性。

函数名称dimensionstdbestmeanworest
Griewank2 0.00201943.527e-050.00177350.0073195
50.037630.029930.0971710.17891
100.0385490.0856370.146690.2594
Alpine2 2.1028e-052.0679e-077.9942e-060.00011011
50.0197920.00511450.0328010.082739
100.182550.238990.510661.0728
Ackley2 4.2601e-066.7944e-077.1477e-061.6723e-05
50.34160.131870.502731.5847
100.754453.59635.06566.5802
Rastrigin2 2.7008e-068.9284e-092.1745e-069.6593e-06
50.897891.36693.01255.0079
104.164510.770719.33830.3427
Sphere2 1.9698e-131.4906e-161.0189e-131.0776e-12
56.1936e-061.8806e-068.2342e-062.8413e-05
100.0351660.0266890.0931650.1588
    通过上表可知,花授粉算法在低维情况下,优化效果显著,可以满足实际效果。花授粉算法随着维度的增加算法的性能和鲁棒性逐渐下降。

3.3 可视化

Sphere函数
在这里插入图片描述
在这里插入图片描述
Ackley函数
在这里插入图片描述
在这里插入图片描述

4. 总结

    总体而言,花授粉优化算法是一个理论简单,易于实现的算法。在编写程序过程中,唯一让我停滞的点是 L L L 步长的计算。最后,在智能优化算法与涌现计算这本书的布谷鸟搜索章节找到答案。下面谈谈花授粉算法的优缺点:
优点:理论简单、易于实现。
缺点:固定概率会影响到一些优化问题解的精度;该算法仍存在一般元启发式优化算法的病症:收敛精度低、对维度敏感等问题。
个人观点:目前,元启发式优化算法的扩展(改进算法)仍在处在精度和维度上的改进,毫无新意可言。但是,这并不意味着该类算法已经完美无瑕。相反,在离散优化问题和约束优化问题上,元启发式算法的解决方案仍存在一定的局限性。

5. 参考文献

[1] YANG XS. Flower pollination algorithm for global optimization[C]//International Conference on Unconventional Computing and Natural Computation.Berlin, Heidelberg: Springer, 2012: 240-249.
[2] 李士勇, 李研, 林永茂. 智能优化算法与涌现计算[M]. 北京: 清华大学出版社, 2019:463-465.

6. 附加

6.1 代码

6.1.1 matlab

clc;
clear;
close all;

best_f = zeros(1, 30);
for epoch=1:30
    %% 问题定义
    n = 50;  % 种群个数
    d = 10;   % 种群维数
    Option = 1;
    [fobj, bound] = Optimizer(Option);
    Lb = bound(1)*ones(1,d);  % 种群下界
    Ub = bound(2)*ones(1,d);  % 种群上界


    %% 初始化种群
    Fitness = zeros(n, 1);
    Sol = zeros(n,d);
    for i=1:n
        Sol(i,:) = Lb + (Ub-Lb).*rand(1,d);
        Fitness(i) = fobj(Sol(i,:));
    end

    %% 计算当前最好的位置
    [fmin,I] = min(Fitness);
    best = Sol(I,:);
    S = Sol;

    %% 设置最大迭代次数
    N_iter = 1000;
    p = 0.8;
    
    best_scores = [fmin];
    
    %% FPA算法迭代
    for it = 1:N_iter
        for i= 1:n
            if rand < p
                %% 异花授粉
                L = Levy(d);
                dS = L.*(Sol(i,:) - best);
                S(i,:) = Sol(i,:) + dS;
                S(i,:) = simplebounds(S(i,:),Lb,Ub); % 越界处理
            else
                %% 自花授粉
                JK = randperm(n);
                S(i,:) = S(i,:) + rand*(Sol(JK(1),:) - Sol(JK(2),:));
                %             S(i,:) = 0.5*Sol(i,:) + rand*random('normal', 0, 1)*(Sol(JK(1),:) - Sol(JK(2),:));
                S(i,:) = simplebounds(S(i,:), Lb, Ub);  % 越界处理
            end

            Fnew=fobj(S(i,:));
            if Fnew <= Fitness(i)
                Sol(i,:) = S(i,:);
                Fitness(i) = Fnew;
            end
    
            if Fnew <= fmin
                best = S(i,:);
                fmin = Fnew;
            end
        end
        best_scores = [best_scores,fmin];
        disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(fmin)]);
        %                 disp(['Bestpop ' num2str(best)]);
        %     hold on;
        %     plot(BestScore(1:it), 'b-')
        %     pause(0.01);
    end
    %         hold off;
    % % 可视化
    % save('fpaSphere.mat','Sol');
    % figure
    % semilogy(BestScore);
    % xlabel('Iteratons', 'fontsize', 10);
    % ylabel('Fitness Value', 'fontsize', 10);
    % xlim([0 1000]);
    %     disp(['FPA算法单独运行30次' num2str(epoch)])
    best_f(epoch) = fmin;
end

disp(["标准差:",num2str(std(best_f))]);
disp(["最优值:",num2str(min(best_f))]);
disp(["平均值:",num2str(mean(best_f))]);
disp(["最差值:",num2str(max(best_f))]);

%% 越界处理
function S=simplebounds(S, lb, ub)
    Tp = S>ub;
    Tm = S<lb;
    S = (S.*(~(Tp+Tm))) + ub.*Tp + lb.*Tm;
end

6.1.2 二维可视化代码

clc;
clear;
close all;


%% 问题定义
n = 50;  % 种群个数
d = 2;   % 种群维数
Option = 1;
[fobj, bound] = Optimizer(Option);
Lb = bound(1)*ones(1,d);  % 种群下界
Ub = bound(2)*ones(1,d);  % 种群上界

%% 函数云图
x = linspace(bound(1), bound(2), 100);
y = linspace(bound(1), bound(2), 100);

[X, Y] = meshgrid(x, y);

for i=1:length(x)
    for j=1:length(y)
        Z(i,j) = fobj([x(i), x(j)]);
    end
end

box on
hold on
contour(X, Y, Z)

%% 初始化种群
Fitness = zeros(n, 1);
Sol = zeros(n,d);
for i=1:n
    Sol(i,:) = Lb + (Ub-Lb).*rand(1,d);
    Fitness(i) = fobj(Sol(i,:));
end


%% 散点图
points = scatter(Sol(:,1), Sol(:,2), 'ko', 'MarkerFaceColor','r');

%% 计算当前最好的位置
[fmin,I] = min(Fitness);
best = Sol(I,:);
S = Sol;

%% 设置最大迭代次数
N_iter = 100;
p = 0.8;
best_scores = [fmin];

%% FPA算法迭代
for it = 1:N_iter
    for i= 1:n
        if rand < p
            % 异花授粉
            L = Levy(d);
            dS = L.*(Sol(i,:) - best);
            S(i,:) = Sol(i,:) + dS;
            S(i,:) = simplebounds(S(i,:),Lb,Ub); % 越界处理
        else
            % 自花授粉
            JK = randperm(n);
            S(i,:) = S(i,:) + rand*(Sol(JK(1),:) - Sol(JK(2),:));
                        S(i,:) = 0.5*Sol(i,:) + rand*random('normal', 0, 1)*(Sol(JK(1),:) - Sol(JK(2),:));
            S(i,:) = simplebounds(S(i,:), Lb, Ub);  % 越界处理
        end

        Fnew=fobj(S(i,:));
        if Fnew <= Fitness(i)
            Sol(i,:) = S(i,:);
            Fitness(i) = Fnew;
        end

        if Fnew <= fmin
            best = S(i,:);
            fmin = Fnew;
        end
    end
    best_scores = [best_scores,fmin];
    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(fmin)]);
    disp(['Bestpop ' num2str(best)]);
    reset(points);
    points = scatter(Sol(:,1), Sol(:,2), 'ko', 'MarkerFaceColor','r');
    pause(0.1);  
    drawnow();
end

6.1.3 python

import numpy as np
from matplotlib import pyplot as plt

plt.style.use('ggplot')

plt.rcParams["font.sans-serif"] = ["Microsoft YaHei"]
plt.rcParams["axes.unicode_minus"] = False


# 目标函数
def default(x):
    y = -8.956 * x[0] + 17.642 * np.log10(x[1]) - 59.182 * x[2] - 52.967
    return y


class FPA:
    def __init__(self, num_pop, dim, ub, lb, f_obj, verbose):
        self.num_pop = num_pop  # 种群数目
        self.dim = dim  # 维数
        self.ub = ub  # 上界
        self.lb = lb  # 下界
        self.pop = np.empty((num_pop, dim))  # 种群
        self.f_obj = f_obj  # 目标函数
        self.f_score = np.empty((num_pop, 1))
        self.verbose = verbose  # 显示

        self.p_best = None  # 最好个体
        self.f_best = None  # 最好适应度值

        self.iter_f_score = []

    def initialize(self):
        """
        种群初始化
        :return: 初始化种群和种群分数
        """
        num_boundary = len(self.ub)
        if num_boundary == 1:
            for i in range(self.num_pop):
                self.pop[i, :] = self.lb + (self.ub - self.lb) * np.random.rand(1, self.dim).flatten()
                self.f_score[i] = self.f_obj(self.pop[i, :])
        else:
            for i in range(self.dim):
                self.pop[:, i] = self.lb[i] + (self.ub[i] - self.lb[i]) * np.random.rand(self.num_pop, 1).flatten()
            for i in range(self.num_pop):
                self.f_score[i] = self.f_obj(self.pop[i, :])
        return [self.pop, self.f_score]

    def get_best(self):
        """
        获取最优解
        :return: 最优索引和最优得分
        """
        idx = np.argmin(self.f_score)
        self.p_best = self.pop[idx, :]
        self.f_best = np.min(self.f_score)
        return [self.p_best, self.f_best]

    def Levy(self, beta=1.5):
        """
        莱维飞行
        :param beta: 固定参数
        :return: 随机数
        """
        sigma = (np.random.gamma(1 + beta) * np.sin(np.pi * beta / 2) / (
                np.random.gamma((1 + beta) / 2) * beta * 2 ** ((beta - 1) / 2))) ** (1 / beta)
        # print(sigma)
        u = np.random.normal(0, sigma, (1, self.dim))
        v = np.random.normal(0, 1, (1, self.dim))
        levy = u / (np.abs(v) ** (1 / beta))
        return levy

    def clip(self, p):
        """
        :param p: 花粉配子
        :return: 越界处理后的花粉配子
        """
        i = p < self.lb
        p[i] = np.array(self.lb)[i]
        j = p > self.ub
        p[j] = np.array(self.ub)[j]
        return p

    def fit(self, max_generation, p=0.8):
        """
        算法迭代
        :param max_generation: 最大迭代次数
        :param p: 切换概率
        :return: 最优值和最优个体
        """
        self.pop, self.f_score = self.initialize()  # 初始化
        self.p_best, self.f_best = self.get_best()  # 当前最优个体

        self.iter_f_score.append(self.f_best)

        # 迭代
        for i in range(max_generation):
            for j in range(self.num_pop):
                if np.random.rand() < p:
                    # 异花授粉公式
                    new_pop = self.pop[j, :] + self.Levy() * (self.pop[j, :] - self.p_best)
                else:
                    # 自花授粉公式
                    idx_set = np.random.choice(range(self.num_pop), 2)
                    new_pop = self.pop[j, :] + np.random.rand(1, self.dim) * (self.pop[idx_set[0], :] -
                                                                              self.pop[idx_set[1], :])

                # 越界处理
                new_pop = self.clip(new_pop.flatten())
                new_f_score = self.f_obj(new_pop)

                # 进化算法
                if new_f_score < self.f_score[j]:
                    self.pop[j, :] = new_pop
                    self.f_score[j] = new_f_score

                if new_f_score < self.f_best:
                    self.p_best = new_pop
                    self.f_best = new_f_score

            self.iter_f_score.append(self.f_best)

            if self.verbose:
                print("============{}/{}==============".format(i, max_generation))
                print(self.f_best)

        return [self.iter_f_score, self.f_best]


# 测试
if __name__ == "__main__":
    n_pop = 50
    d = 3
    upper = [1, 2400, 1]
    lower = [0, 500, 0]

    max_iter = 1000

    fpa = FPA(n_pop, d, upper, lower, default, True)
    iter_score, _ = fpa.fit(max_iter)

    xx = [int(i) for i in range(0, 1001, 50)]
    yy = np.array(iter_score)
    print(fpa.p_best)

    plt.figure()
    plt.title("收敛曲线")
    plt.plot(xx, yy[xx], c='k', marker='.')
    plt.xlim([0, 1000])
    plt.xlabel("迭代次数")
    plt.ylabel("适应度值")
    plt.show()

6.2 流程图

在这里插入图片描述

  • 14
    点赞
  • 50
    收藏
    觉得还不错? 一键收藏
  • 9
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 9
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值