粒子群算法的函数极值求解

粒子群算法（PSO）是一种基于群体的随机优化技术，与其他基于群体的进化算法相比，他们均初始化为一组随机解，通过迭代搜寻最优解。不同的是，进化算法遵循适者生存原则，而PSO模拟社会，将每个可能产生的解表述为群中的一个微粒，每个微粒都具有自己的位置向量和速度向量，以及一个由目标函数决定的适应度，所有微粒在搜索空间中以一定的速度飞行，通过追随当前搜索到的最优值来搜索全局最优值。
——《Matlab智能算法 超级学习手册》

1. 一般粒子群算法的基本原理

设在一个S维的搜索空间中，有m个粒子组成一个群体。每个粒子拥有自己的位置坐标以及飞行速度，分别记作和，。每个粒子所对应的适应值为。记第i个粒子迄今搜素到的最优位置为：,所有粒子至今搜索到的最优位置为。

设所求目标为的极大值，则微粒i当前最好位置为：

迭代时，使用下式对粒子进行操作：

为t+1步时，第i个粒子的速度。1式第一项中，w为非负常量，称作动力常量，控制前一速度对后一速度的影响；第二项中，c1为自我学习因数，控制粒子飞向自身搜索最优解；第三项中，c2为社会影响因数，控制粒子飞向粒子社会迄今搜索到的最优解。为了防止粒子飞行飞出定义域，通常为粒子速度和粒子位置设置限制，即和。

综上，一般PSO求解极值问题的步骤为：

1，初始化粒子个数m，各粒子初始坐标，速度。
2，计算每个粒子的适应值。
3，对每个粒子将其适应值和其经历过的最好的位置的适应值进行比较，若较好，则将其本次位置作为 当前的最好位置。
4，对每个粒子将其适应值和全局经历过的最好位置的适应值进行比较，若较好，则将其本次位置作为 当前的全局最好位置。
5，根据式1和式2分别对粒子的速度和位置进行更新。
6，若满足终止条件，则输出解，否则返回第二步。
——《Matlab智能算法 超级学习手册》

PSO求解极值问题流程图

*图片感谢胡哥

1. 使用粒子群算法求解函数极值

在工程实际问题中，经常需要求解一些非线性方程的根，这些非线性方程根很难找出解析解，粒子群算法则很好的解决了这个问题。

• 使用一般粒子群算法求解一元函数最大值。

clc
clear
close all;
%% 目标函数定义
f=@(x)11.*sin(2.*x)+7*cos(5.*x);    %函数定义
lb=-3;                              %函数定义域下界
ub=3;                               %函数定义域上界
figure(1);                          %开辟图形窗口
subplot(1,2,1);                     %分割图形窗口
ezplot(f,[lb ub]);                  %绘制函数图像
xlabel('x');                        %设置坐标轴x
ylabel('y');                        %设置坐标轴y
hold on;                            %绘图保持命令
%% 算法参数初始化
m=10;                               %粒子数
vlimit=1.2;                         %速度界限
N=100;                              %迭代次数
w=0.9;                              %惯性影响因数
c1=0.5;                             %自我学习因数
c2=0.5;                             %社会影响因数

x=rand(m,1)*(ub-lb)+lb;             %粒子坐标初始化
v=-vlimit+2*vlimit*rand(m,1);       %粒子速度初始化
y=f(x);                             %粒子初始化后适应值
xbest=x;                            %粒子自身学习到的最优坐标
ybest=y;                            %粒子自身学习到的最优值
ysbest=max(y);                      %粒子群社会当前最优值
xsbest=x(y==max(y));                %粒子群社会当前最优值坐标
BestY=[];                           %存放每次迭代的社会最优适应值
BestX=xsbest;                       %存放每次迭代的社会最优适应值对应的x坐标
h=scatter(x,y,80,'*r');             %绘制散点图并返回图形句柄
%% 迭代计算
while N
   v=w.*v+c1.*(xbest-x)+c2.*(BestX-x); %粒子速度更新
   for i=1:m                           %粒子速度上下界约束
       if     v(i,1)>vlimit            %若粒子速度正向过大
              v(i,1)=vlimit;           %粒子速度为正向最大值
       elseif vlimit<-v(i,1)           %若粒子速度反向过大
              v(i,1)=-vlimit;          %粒子速度为反向最大值
       end
   end
   x=x+v;                              %粒子位置更新
   for i=1:m                           %粒子位置上下界约束
       if     x(i,1)>ub                %若粒子位置正向过大
              x(i,1)=ub;               %粒子位置为定义域上界
       elseif x(i,1)<lb                %若粒子位置反向过大
              x(i,1)=lb;               %粒子位置为定义域下界
       end
   end
   
   y=f(x);                             %重新计算各粒子适应值
   
   for i=1:m                           %遍历所有粒子
       if y(i)>ybest(i)                %若当前结果较好
          ybest(i)=y(i);               %更新第i个粒子的自身最好经验值
          xbest(i)=x(i);               %更新第i个粒子的自身最好经验坐标
       end
   end
   
   BestY=[BestY,max(ybest)];           %更新每次迭代的社会最优适应值
   for i=1:m
       if BestY(end)==max(BestY)       %本次是最优则更新,本次不是最优不更新
           if ybest(i)==max(ybest)
               BestX=xbest(i);         %更新每次迭代的社会最优适应值对应的x坐标
               break;
           end
       end
   end                          
   pause(0.01);                        %延时观察散点图趋势
   h.XData=x;                          %更新散点图x坐标
   h.YData=y;                          %更新散点图y坐标
   N=N-1;                              %循环计数
end
%% 结果输出
plot(BestX,BestY(end),'bo')           %绘制最终结果
subplot(1,2,2);                       %分割图形窗口
plot(BestY);                          %绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');                    %设置坐标轴x
ylabel('适应度');                      %设置坐标轴y
disp('最佳的位置是：'); disp(BestX)      %结果输出
disp('此时最优值是：'); disp(BestY(end)) %结果输出

正常的运行结果是这样：

当然，由于只是一般的PSO，算法容易陷入局部最优解。

• 使用一般粒子群算法求解二元函数最大值。

clc
clear
close all;
%% 目标函数定义
xlb=-5;                              %函数x定义域下界
xub=5;                               %函数x定义域上界
ylb=-5;                              %函数y定义域下界
yub=5;                               %函数y定义域上界
lx=xlb:0.05:xub;                     %x向量
ly=ylb:0.05:yub;                     %y向量
[X,Y]=meshgrid(lx,ly);               %x,y向量构建网格采样点
F=Y.*sin(10.*X)-X.*sin(Y);           %二维函数定义
%F=50-X.^2-Y.^2;
figure(1);                           %开辟图形窗口
subplot(1,2,1);                      %分割图形窗口
mesh(X,Y,F);                         %二维函数绘制
xlabel('x');                         %设置坐标轴x
ylabel('y');                         %设置坐标轴y
zlabel('z');                         %设置坐标轴z
hold on;                             %绘图保持命令
%% 算法参数初始化
m=50;                                 %粒子数
s=2;                                  %维数
vxlimit=(xub-xlb)*0.05;               %x速度界限
vylimit=(yub-ylb)*0.05;               %y速度界限
N=100;                                %迭代次数
w=0.9;                                %惯性影响因数
c1=0.5;                               %自我学习因数
c2=0.5;                               %社会影响因数

x=rand(m,1)*(xub-xlb)+xlb;            %粒子x坐标初始化
y=rand(m,1)*(yub-ylb)+ylb;            %粒子y坐标初始化
vx=-vxlimit+2*vxlimit*rand(m,1);      %粒子x向速度初始化
vy=-vylimit+2*vylimit*rand(m,1);      %粒子y向速度初始化

z=y.*sin(10.*x)-x.*sin(y);                       %粒子初始化后适应值(m个粒子，m个适应值)

xbest=x;                              %粒子自身学习到的最优x坐标
ybest=y;                              %粒子自身学习到的最优y坐标
zbest=z;                              %粒子自身学习到的最优z值

zsbest=max(z);                        %粒子群社会当前最优值
xsbest=x(z==max(z));                  %粒子群社会当前最优值x坐标
ysbest=y(z==max(z));                  %粒子群社会当前最优值y坐标

BestZ=[];                             %存放每次迭代的社会最优适应值
BestX=xsbest;                         %存放每次迭代的社会最优适应值对应的x坐标
BestY=ysbest;                         %存放每次迭代的社会最优适应值对应的y坐标
h=scatter3(x,y,z,80,'*r');            %绘制散点图并返回图形句柄
%% 迭代计算
while N
   vx=w.*vx+c1.*(xbest-x)+c2.*(BestX-x);   %粒子x速度更新
   vy=w.*vy+c1.*(ybest-y)+c2.*(BestY-y);   %粒子x速度更新  
   
   for i=1:m                               %粒子x速度上下界约束
       if     vx(i,1)>vxlimit              %若粒子x速度正向过大
              vx(i,1)=vxlimit;             %粒子x速度为正向最大值
       elseif vxlimit<-vx(i,1)             %若粒子x速度反向过大
              vx(i,1)=-vxlimit;            %粒子x速度为反向最大值
       end
       
       if     vy(i,1)>vylimit              %若粒子y速度正向过大
              vy(i,1)=vylimit;             %粒子y速度为正向最大值
       elseif vylimit<-vy(i,1)             %若粒子y速度反向过大
              vy(i,1)=-vylimit;            %粒子y速度为反向最大值
       end       
   end
   
   x=x+vx;                                 %粒子x位置更新
   y=y+vy;                                 %粒子y位置更新
   
   for i=1:m                               %粒子位置上下界约束
       if     x(i,1)>xub                   %若粒子位置正向过大
              x(i,1)=xub;                  %粒子位置为定义域上界
       elseif x(i,1)<xlb                   %若粒子位置反向过大
              x(i,1)=xlb;                  %粒子位置为定义域下界
       end
       
       if     y(i,1)>yub                   %若粒子y位置正向过大
              y(i,1)=yub;                  %粒子y位置为定义域上界
       elseif y(i,1)<ylb                   %若粒子y位置反向过大
              y(i,1)=ylb;                  %粒子y位置为定义域下界
       end       
   end
   
   z=y.*sin(10.*x)-x.*sin(y);              %重新计算各粒子适应值
   
   for i=1:m                               %遍历所有粒子
       if z(i)>zbest(i)                    %若当前结果较好
          zbest(i)=z(i);                   %更新第i个粒子的自身最好经验z值
          xbest(i)=x(i);                   %更新第i个粒子的自身最好经验x坐标
          ybest(i)=y(i);                   %更新第i个粒子的自身最好经验y坐标          
       end
   end
   
   BestZ=[BestZ,max(zbest)];               %更新每次迭代的社会最优适应值
   for i=1:m
       if BestZ(end)==max(BestZ)           %本次是最优则更新,本次不是最优不更新
           if zbest(i)==max(zbest)
               BestX=xbest(i);             %更新每次迭代的社会最优适应值对应的x坐标
               BestY=ybest(i);             %更新每次迭代的社会最优适应值对应的y坐标               
               break;
           end
       end
   end                          
%    pause(0.01);                        %延时观察散点图趋势
   h.XData=x;                          %更新散点图x坐标
   h.YData=y;                          %更新散点图y坐标
   h.ZData=z;                          %更新散点图y坐标   
   N=N-1;                              %循环计数
end
%% 结果输出
plot3(BestX,BestY,BestZ(end),'bo')     %绘制最终结果
subplot(1,2,2);                        %分割图形窗口
plot(BestZ);                           %绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');                     %设置坐标轴x
ylabel('适应度');                       %设置坐标轴y
disp('最佳的x位置是：'); disp(BestX)      %结果输出
disp('最佳的y位置是：'); disp(BestY)      %结果输出
disp('此时最优z值是：'); disp(BestZ(end)) %结果输出

运行结果如下：

但仍易于陷入局部最优解: