粒子群算法

粒子群算法(particle swarmoptimization，PSO)由Kennedy和Eberhart在1995年提出，该算法模拟鸟集群飞行觅食的行为，鸟之间通过集体的协作使群体达到最优目的，是一种基于Swarm Intelligence的优化方法。同遗传算法类似，也是一种基于群体叠代的，但并没有遗传算法用的交叉以及变异，而是粒子在解空间追随最优的粒子进行搜索。PSO的优势在于简单容易实现同时又有深刻的智能背景，既适合科学研究，又特别适合工程应用，并且没有许多参数需要调整。

 基本PSO算法

    D维空间中，有N个粒子；

    粒子i位置：xi=(xi1,xi2,…xiD)，将xi代入适应函数f(xi)求适应值；

    粒子i速度：vi=(vi1,vi2,…viD)

    粒子i个体经历过的最好位置：pbesti=(pi1,pi2,…piD)

    种群所经历过的最好位置：gbest=(g1,g2,…gD)

  

   

   

   粒子群算法关键参数:

    (1)粒子群规模N : 一般20-50

    (2)惯性权值w:一般为［０.８,１.２］，可以固定也可以动态变化

　(3) 加速常数(学习因子)c1和c2: 一般设置c1=c2

    (4)粒子最大速度Vmax: 控制范围[-Vmax,Vmax]

    (5)停止准则: 最大迭代次数、计算精度或最优解的最大停滞步数

    (6)边界条件处理：比如通过设置最大位置Xmax和Vmax.

  

  利用粒子群算法求函数f(x,y)=3cos(xy)+x+y^2的最小值，其中x的取值范围为[-4,4]，y的取值范围为[-4,4]。

%%%求函数f(x,y)=3cos(xy)+x+y^2的最小值，其中x的取值范围为[-4,4]，y的取值范围为[-4,4]。
%粒子个数N=100
%粒子维数D=2
%迭代次数T=200
%学习因子c1=c2=1.5
%惯性权重最大值为Wmax=0.8 惯性权重最小值为Wmin=0.4
%位置最大值为Xmax=4,位置最小值为Xmin=-4
%速度最大值为Vmax=1,速度最小值为Vmin=-1


%%%%%%粒子群算法求函数极值%%%%%%%%%%%%%%%%%%
%%%%%%初始化%%%%%%%%%%%%%%%%%%
clear all; %清除所有变量
close all; %清图
clc ;      %清屏
N=100;    %粒子个数
D=2;     %粒子维数
T=200;    %最大迭代次数
c1=1.5;   %学习因子c1
c2=1.5;   %学习因子c2
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)=func4(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(func4(x(j,:)<pbest(j)))
           p(j,:)=x(j,:);
           pbest(j)=func4(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;
    disp(i);
end

 %%%%%%输出最优值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 disp('最优个体')
 disp(g);         %最优个体
 disp('最优值')
 disp(gb(end));    %最优值
 
 %%%%%%画图%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure
 plot(gb);
 xlabel('迭代次数');
 ylabel('适应度值');
 title('适应度进化曲线');
 
 
 
 

%%%%%%%求解的函数%%%%%%%%%%%%%%%%%%
function value=func4(x)
value=3*cos(x(1)*x(2))+x(1)+x(2)^2;
end