粒子群（PSO）算法学习及MATLAB实现

粒子群算法的简介及理解

基本思想

粒子群算法的思想源于对鸟群觅食行为的研究，鸟群通过集体的信息共享使群体找到最优的目的地。设想这样一个场景：鸟群在森林中随机搜索食物，它们想要找到食物量最多的位置。但是所有的鸟都不知道食物具体在哪个位置，只能感受到食物大概在哪个方向。每只鸟沿着自己判定的方向进行搜索，并在搜索的过程中记录自己曾经找到过食物且量最多的位置，同时所有的鸟都共享自己每一次发现食物的位置以及食物的量，这样鸟群就知道当前在哪个位置食物的量最多。在搜索的过程中每只鸟都会根据自己记忆中食物量最多的位置和当前鸟群记录的食物量最多的位置调整自己接下来搜索的方向。鸟群经过一段时间的搜索后就可以找到森林中哪个位置的食物量最多（全局最优解）。

MATLAB实现

c1=2;%学习因子
c2=2;%学习因子
Dimension=20;
Size=30;
Tmax=500;
Velocity_max=1200;%粒子最大速度

F_n=2;%测试函数名

Fun_Ub=600;%函数上下界
Fun_Lb=-600;
Position=zeros(Dimension,Size);%粒子位置
Velocity=zeros(Dimension,Size);%粒子速度
Vmax(1:Dimension)=Velocity_max;%粒子速度上下界
Vmin(1:Dimension)=-Velocity_max;
Xmax(1:Dimension)=Fun_Ub;%粒子位置上下界，即函数自变量的上下界
Xmin(1:Dimension)=Fun_Lb;

%初始化
%产生一个20行60列的矩阵;
%前30列表示随机产生的粒子位置、后30列表示随机产生的粒子速度;
[Position,Velocity]=Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin);
Pbest_position=Position;%粒子的历史最优位置，初始值为粒子的起始位置，存储每个粒子的历史最优位置
Gbest_position=zeros(Dimension,1);%全局最优的那个粒子所在位置，初始值认为是第1个粒子

for j=1:Size
    Pos=Position(:,j);%取第j列，即第j个粒子的位置
    fz(j)=Fitness_Function(Pos,F_n,Dimension);%计算第j个粒子的适应值
end
[Gbest_Fitness,I]=min(fz);%求出所有适应值中最小的那个适应值，并获得该粒子的位置
Gbest_position=Position(:,I);%取最小适应值的那个粒子的位置，即I列

for itrtn=1:Tmax
time(itrtn)=itrtn;

Weight=1;
r1=rand(1);
r2=rand(1);
%速度更新公式
for i=1:Size
   Velocity(:,i)=Weight*Velocity(:,i)+c1*r1*(Pbest_position(:,i)-Position(:,i))+c2*r2*(Gbest_position-Position(:,i));
end
%限制速度边界
for i=1:Size
    for row=1:Dimension
        if Velocity(row,i)>Vmax(row)
            Veloctity(row,i)=Vmax(row);
        elseif Velocity(row,i)<Vmin(row)
            Veloctity(row,i)=Vmin(row);
        else
        end
    end
end

Position=Position+Velocity;

%限制位置边界
for i=1:Size
    for row=1:Dimension
        if Position(row,i)>Xmax(row)
            Position(row,i)=Xmax(row);
        elseif Position(row,i)<Xmin(row)
            Position(row,i)=Xmin(row);
        else
        end
    end
end

  for j=1:Size
     P_position=Position(:,j)';%取一个粒子的位置
     fitness_p(j)=Fitness_Function(P_position,F_n,Dimension);
     if fitness_p(j)< fz(j) %粒子的适应值比运动之前的适应值要好，更新原来的适应值
         Pbest_position(:,j)=Position(:,j);
         fz(j)=fitness_p(j);
     end
     if fitness_p(j)<Gbest_Fitness
         Gbest_Fitness=fitness_p(j);
     end
  end
  [Gbest_Fitness_new,I]=min(fz);%更新后的所有粒子的适应值，取最小的那个，并求出其编号
   Best_fitness(itrtn)=Gbest_Fitness_new; %记录每一代的最好适应值
   Gbest_position=Pbest_position(:,I);%最好适应值对应的个体所在位置
end
plot(time,Best_fitness);
xlabel('迭代的次数');ylabel('适应度值P_g');

%% 初始化
function [Position,Velocity] = Initial_position_velocity(Dimension,Size,Xmax,Xmin,Vmax,Vmin)
  for i=1:Dimension
      Position(i,:)=Xmin(i)+(Xmax(i)-Xmin(i))*rand(1,Size); % 产生合理范围内的随机位置，rand(1,Size)用于产生一行Size个随机数
      Velocity(i,:)=Vmin(i)+(Vmax(i)-Vmin(i))*rand(1,Size);
  end
end

%% 适应值计算
function Fitness=Fitness_Function(Pos,F_n,Dimension)
 switch F_n
    case 1
        Func_Sphere=Pos(:)'*Pos(:);
        Fitness=Func_Sphere;
    case 2
        res1=Pos(:)'*Pos(:)/4000;
        res2=1;
        for row=1:Dimension
            res2=res2*cos(Pos(row)/sqrt(row));
        end
        Func_Griewank=res1-res2+1;
        Fitness=Func_Griewank;
 end
end

参考链接：
https://cloud.tencent.com/developer/article/1064332
https://zhuanlan.zhihu.com/p/346355572
https://blog.csdn.net/daaikuaichuan/article/details/81382794