本文已参加「新人创造礼」活动，一同敞开创造之路

读书使人充分，评论使人机敏，笔记使人精确…。凡有所学，皆成性情。 ———— （英国）培根

粒子群算法理论

鸟类在捕食进程中，鸟群成员能够通过个别之间的信息沟通与共享取得其他成员的发现与飞翔阅历。在食物源零星分布而且不行猜测的条件下，这种协作机制所带来的优势是决定性的，远远大于对食物的竞争所引起的下风。粒子群算法受鸟类捕食行为的启示并对这种行为进行仿照，将优化问题的查找空间类比于鸟类的飞翔空间，将每只鸟抽象为一个粒子，粒子无质量、无体积，用以表征问题的一个可行解，优化问题所要查找到的最优解则等同于鸟类寻找的食物源。粒子群算法为每个粒子制定了与鸟类运动类似的简单行为规则，使整个粒子群的运动表现出与鸟类捕食类似的特性，然后能够求解杂乱的优化问题

关键参数说明

[==求最小值==]核算函数f(x)=∑i=1nxi2(−20≤xi≤20)f(x)=\sum_{i=1}^nx_{i}^2\quad(-20 \leq x_{i} \leq 20)f(x)=∑i=1n​xi2​(−20≤xi​≤20)​ 的最小值，其间个别x的维数n=10。这是一个简单的平方和函数，只有一个极小点x=（0，0，…，0），理论最小值f（0，0，…，0）=0

解：仿真进程如下：

（1）初始化集体粒子个数为N=100，粒子维数为D=10，最大迭代次数为T=200，学习因子c1=c2=1.5，惯性权重为w=0.8，方位最大值为xmax=20x_{max}=20xmax​=20，方位最小值为xmin=20x_{min}=20xmin​=20，速度最大值为vmax=10v_{max}=10vmax​=10，速度最小值为vmin=−10v_{min}=-10vmin​=−10。

（2）初始化种群粒子方位x和速度v，粒子个别最优方位p和最优值pbest，以及粒子群大局最优方位g和最优值gbest。

（3）更新方位x和速度值v，并进行边界条件处理，判别是否替换粒子个别最优方位$p$和最优值pbestp_{best}pbest​、粒子群大局最优方位$g$和最优值gbestg_{best}gbest​​。

（4）判别是否满意终止条件：若满意，则结束查找进程，输出优化值；若不满意，则继续进行迭代优化。

优化结束后，其习惯度进化曲线如图所示，优化后的成果为x=［-0.63250.1572 -0.4814 0.1091 -0.3154 0.2236 -0.3991 0.5907 0.0221 -0.1172］10−410^{-4}10−4​，函数$f(x)$​的最小值为1.3410−81.3410^{-8}1.3410−8​。

MATLAB源程序如下：

%%%%%%%%%%%%%%%%%粒子群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %铲除一切变量
close all;              %清图
clc;                    %清屏
N=100;                  %集体粒子个数
D=10;                   %粒子维数
T=200;                  %最大迭代次数
c1=1.5;                 %学习因子1
c2=1.5;                 %学习因子2
w=0.8;                  %惯性权重
Xmax=20;                %方位最大值
Xmin=-20;               %方位最小值
Vmax=10;                %速度最大值
Vmin=-10;               %速度最小值
%%%%%%%%%%%%%%%%初始化种群个别（限制方位和速度）%%%%%%%%%%%%%%%%
x=rand(N,D) * (Xmax-Xmin)+Xmin;
v=rand(N,D) * (Vmax-Vmin)+Vmin;
%%%%%%%%%%%%%%%%%%初始化个别最优方位和最优值%%%%%%%%%%%%%%%%%%%
p=x;
pbest=ones(N,1);
for i=1:N
    pbest(i)=func1(x(i,:));
end
%%%%%%%%%%%%%%%%%%%初始化大局最优方位和最优值%%%%%%%%%%%%%%%%%%
g=ones(1,D);
gbest=inf;
for i=1:N
    if(pbest(i)<gbest)
        g=p(i,:);
        gbest=pbest(i);
    end
end
gb=ones(1,T);
%%%%%%%%%%%依照公式顺次迭代直到满意精度或许迭代次数%%%%%%%%%%%%%
for i=1:T
    for j=1:N
        %%%%%%%%%%%%%%更新个别最优方位和最优值%%%%%%%%%%%%%%%%%
        if (func1(x(j,:))<pbest(j))
            p(j,:)=x(j,:);
            pbest(j)=func1(x(j,:));
        end
        %%%%%%%%%%%%%%%%更新大局最优方位和最优值%%%%%%%%%%%%%%%
        if(pbest(j)<gbest)
            g=p(j,:);
            gbest=pbest(j);
        end
        %%%%%%%%%%%%%%%%%跟新方位和速度值%%%%%%%%%%%%%%%%%%%%%
        v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
            +c2*rand*(g-x(j,:));
        x(j,:)=x(j,:)+v(j,:);
        %%%%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D
            if (v(j,ii)>Vmax)  |  (v(j,ii)< Vmin)
                v(j,ii)=rand * (Vmax-Vmin)+Vmin;
            end
            if (x(j,ii)>Xmax)  |  (x(j,ii)< Xmin)
                x(j,ii)=rand * (Xmax-Xmin)+Xmin;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%记载历代大局最优值%%%%%%%%%%%%%%%%%%%%%
    gb(i)=gbest;
end
g;                         %最优个别         
gb(end);                   %最优值
figure
plot(gb)
xlabel('迭代次数');
ylabel('习惯度值');
title('习惯度进化曲线')
%%%%%%%%%%%%%%%%%%%习惯度函数%%%%%%%%%%%%%%%%%%%%
function result=func1(x)
summ=sum(x.^2);
result=summ;
end

[==求最小值==]函数f(x,y)=3cos(xy)+x+y2f(x,y)=3cos(xy)+x+y^2f(x,y)=3cos(xy)+x+y2​​的最小值，其间x的取值规模为［-4，4］，y的取值规模为［-4，4］。这是一个有多个部分极值的函数，其函数值图形如图所示，其MATLAB完成程序如下

%%%%%%%%%f(x，y)=3*cos(x*y)+x+y*y%%%%%%%%%%
clear all;              %铲除一切变量
close all;              %清图
clc;                    %清屏
x=-4:0.02:4;
y=-4:0.02:4;
N=size(x,2);
for i=1:N
    for j=1:N
        z(i,j)=3*cos(x(i)*y(j))+x(i)+y(j)*y(j);
    end
end
mesh(x,y,z)
xlabel('x')
ylabel('y')

解：仿真进程如下：

（1）初始化集体粒子个数为N=100，粒子维数为D=2，最大迭代次数为T=200，学习因子c1=c2=1.5，惯性权重最大值为wmax=0.8w_{max}=0.8wmax​=0.8，惯性权重最小值为wmin=0.4w_{min}=0.4wmin​=0.4，方位最大值为xmax=4x_{max}=4xmax​=4，方位最小值为xmin=−4x_{min}=-4xmin​=−4，速度最大值为vmax=1v_{max}=1vmax​=1，速度最小值为vmin=−1v_{min}=-1vmin​=−1。

（2）初始化种群粒子方位x和速度v，粒子个别最优方位p和最优值pbest，粒子群大局最优方位g和最优值gbest。

（3）核算动态惯性权重值w，更新方位x和速度值v，并进行边界条件处理，判别是否替换粒子个别最优方位$p$和最优值pbestp_{best}pbest​，以及粒子群大局最优方位$g$和最优值gbestg_{best}gbest​。

（4）判别是否满意终止条件：若满意，则结束查找进程，输出优化值；若不满意，则继续进行迭代优化

优化结束后，其习惯度进化曲线如图6.4所示。优化后的成果为：在x=-4.0000，y=-0.7578时，函数$f(x)$取得最小值-6.408

MATLAB源程序如下：

%%%%%%%%%%%%%%%%%粒子群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %铲除一切变量
close all;              %清图
clc;                    %清屏
N=100;                  %集体粒子个数
D=2;                    %粒子维数
T=200;                  %最大迭代次数
c1=1.5;                 %学习因子1
c2=1.5;                 %学习因子2
Wmax=0.8;               %惯性权重最大值
Wmin=0.4;               %惯性权重最小值
Xmax=4;                 %方位最大值
Xmin=-4;                %方位最小值
Vmax=1;                 %速度最大值
Vmin=-1;                %速度最小值
%%%%%%%%%%%%%%%%初始化种群个别（限制方位和速度）%%%%%%%%%%%%%%%%
x=rand(N,D) * (Xmax-Xmin)+Xmin;
v=rand(N,D) * (Vmax-Vmin)+Vmin;
%%%%%%%%%%%%%%%%%%初始化个别最优方位和最优值%%%%%%%%%%%%%%%%%%%
p=x;
pbest=ones(N,1);
for i=1:N
    pbest(i)=func2(x(i,:));
end
%%%%%%%%%%%%%%%%%%%初始化大局最优方位和最优值%%%%%%%%%%%%%%%%%%
g=ones(1,D);
gbest=inf;
for i=1:N
    if(pbest(i)<gbest)
        g=p(i,:);
        gbest=pbest(i);
    end
end
gb=ones(1,T);
%%%%%%%%%%%依照公式顺次迭代直到满意精度或许迭代次数%%%%%%%%%%%%%
for i=1:T
    for j=1:N
        %%%%%%%%%%%%%%更新个别最优方位和最优值%%%%%%%%%%%%%%%%%
        if (func2(x(j,:))<pbest(j))
            p(j,:)=x(j,:);
            pbest(j)=func2(x(j,:));
        end
        %%%%%%%%%%%%%%%%更新大局最优方位和最优值%%%%%%%%%%%%%%%
        if(pbest(j)<gbest)
            g=p(j,:);
            gbest=pbest(j);
        end
        %%%%%%%%%%%%%%%%核算动态惯性权重值%%%%%%%%%%%%%%%%%%%%
        w=Wmax-(Wmax-Wmin)*i/T;
        %%%%%%%%%%%%%%%%%跟新方位和速度值%%%%%%%%%%%%%%%%%%%%%
        v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
            +c2*rand*(g-x(j,:));
        x(j,:)=x(j,:)+v(j,:);
        %%%%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D
            if (v(j,ii)>Vmax)  |  (v(j,ii)< Vmin)
                v(j,ii)=rand * (Vmax-Vmin)+Vmin;
            end
            if (x(j,ii)>Xmax)  |  (x(j,ii)< Xmin)
                x(j,ii)=rand * (Xmax-Xmin)+Xmin;
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%记载历代大局最优值%%%%%%%%%%%%%%%%%%%%%
    gb(i)=gbest;
end
g;                         %最优个别         
gb(end);                   %最优值
figure
plot(gb)
xlabel('迭代次数');
ylabel('习惯度值');
title('习惯度进化曲线')
%%%%%%%%%%%%%%%%%%%%%习惯度函数%%%%%%%%%%%%%%%%%%%%%%%
function value=func2(x)
value=3*cos(x(1)*x(2))+x(1)+x(2)^2;
end

[==求最小值==]用离散粒子群算法求函数$f(x)=x+6sin(4x)+9cos(5x)$的最小值，其间x的取值规模为［0，9］。这是一个有多个部分极值的函数，其函数值图形如图所示，其MATLAB完成程序如下：

%%%%%%%%%f(x)=x+6sin(4x)+9cos(5x)%%%%%%%%%%
clear all;              %铲除一切变量
close all;              %清图
clc;                    %清屏
x=0:0.01:9;
y=x+6*sin(4*x)+9*cos(5*x);
plot(x,y)
xlabel('x')
ylabel('f(x)')
title('f(x)=x+6sin(4x)+9cos(5x)')

解：仿真进程如下：

(1）初始化集体粒子个数为N=100，粒子维数为D=20，最大迭代次数为T=200，学习因子c1=c2=1.5，惯性权重最大值为wmax=0.8w_{max}=0.8wmax​=0.8​​​，惯性权重最小值为wmin=0.4w_{min}=0.4wmin​=0.4​​​，方位最大值为xmax=9x_{max}=9xmax​=9​​​，方位最小值为xmin=0x_{min}=0xmin​=0，速度最大值为vmax=10v_{max}=10vmax​=10，速度最小值为vmin=−10v_{min}=-10vmin​=−10​​​。

（2）初始化速度v和二进制编码的种群粒子方位x，核算习惯度值，取得粒子个别最优方位p和最优值pbest，以及粒子群大局最优方位g和最优值gbest。

（3）核算动态惯性权重值w，更新方位x和速度值v，并进行边界条件处理，判别是否替换粒子个别最优方位$p$和最优值pbestp_{best}pbest​，以及粒子群大局最优方位$g$和最优值gbestg_{best}gbest​。

（4）判别是否满意终止条件：若满意，则结束查找进程，输出优化值；若不满意，则继续进行迭代优化。

其习惯度进化曲线如图所示。优化后的成果为：当x=4.37时，函数$f(x)$取得最小值-10.42。

MATLAB源程序如下：

%%%%%%%%%%%%%%%%%离散粒子群算法求函数极值%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;              %铲除一切变量
close all;              %清图
clc;                    %清屏
N=100;                  %集体粒子个数
D=20;                   %粒子维数
T=200;                  %最大迭代次数
c1=1.5;                 %学习因子1
c2=1.5;                 %学习因子2
Wmax=0.8;               %惯性权重最大值
Wmin=0.4;               %惯性权重最小值
Xs=9;                   %方位最大值
Xx=0;                   %方位最小值
Vmax=10;                %速度最大值
Vmin=-10;               %速度最小值
%%%%%%%%%%%%%%%%初始化种群个别（限制方位和速度）%%%%%%%%%%%%%%%%
x=randi([0,1],N,D);        %随机取得二进制编码的初始种群
v=rand(N,D) * (Vmax-Vmin)+Vmin;
%%%%%%%%%%%%%%%%%%初始化个别最优方位和最优值%%%%%%%%%%%%%%%%%%%
p=x;
pbest=ones(N,1);
for i=1:N
    pbest(i)= func3(x(i,:),Xs,Xx);
end
%%%%%%%%%%%%%%%%%%%初始化大局最优方位和最优值%%%%%%%%%%%%%%%%%%
g=ones(1,D);
gbest=inf;
for i=1:N
    if(pbest(i)<gbest)
        g=p(i,:);
        gbest=pbest(i);
    end
end
gb=ones(1,T);
%%%%%%%%%%%依照公式顺次迭代直到满意精度或许迭代次数%%%%%%%%%%%%%
for i=1:T
    for j=1:N
        %%%%%%%%%%%%%%更新个别最优方位和最优值%%%%%%%%%%%%%%%%%
        if (func3(x(j,:),Xs,Xx)<pbest(j))
            p(j,:)=x(j,:);
            pbest(j)=func3(x(j,:),Xs,Xx);
        end
        %%%%%%%%%%%%%%%%更新大局最优方位和最优值%%%%%%%%%%%%%%%
        if(pbest(j)<gbest)
            g=p(j,:);
            gbest=pbest(j);
        end
        %%%%%%%%%%%%%%%%核算动态惯性权重值%%%%%%%%%%%%%%%%%%%%
        w=Wmax-(Wmax-Wmin)*i/T;
        %%%%%%%%%%%%%%%%%跟新方位和速度值%%%%%%%%%%%%%%%%%%%%%
        v(j,:)=w*v(j,:)+c1*rand*(p(j,:)-x(j,:))...
            +c2*rand*(g-x(j,:));
        %%%%%%%%%%%%%%%%%%%%边界条件处理%%%%%%%%%%%%%%%%%%%%%%
        for ii=1:D
            if (v(j,ii)>Vmax)  |  (v(j,ii)< Vmin)
                v(j,ii)=rand * (Vmax-Vmin)+Vmin;
            end
        end    
        vx(j,:)=1./(1+exp(-v(j,:)));
        for jj=1:D
            if vx(j,jj)>rand
                x(j,jj)=1;
            else
                x(j,jj)=0;
            end
        end      
    end
    %%%%%%%%%%%%%%%%%%%%记载历代大局最优值%%%%%%%%%%%%%%%%%%%%%
    gb(i)=gbest;
end
g;                       %最优个别 
m=0;
for j=1:D
    m=g(j)*2^(j-1)+m;
end
f1=Xx+m*(Xs-Xx)/(2^D-1); %最优值
figure
plot(gb)
xlabel('迭代次数');
ylabel('习惯度值');
title('习惯度进化曲线')
%%%%%%%%%%%%%%%%%%%%%%%%%习惯度函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result=func3(x,Xs,Xx)
m=0;
D=length(x);
for j=1:D
    m=x(j)*2^(j-1)+m;
end
f=Xx+m*(Xs-Xx)/(2^D-1);            %译码成十进制数
fit= f+6*sin(4*f)+9*cos(5*f);
result=fit;
end

[1]包子阳 余继周 杨彬.智能优化算法及其MATLAB实例[M]电子工业出版社，2021