teachingpso.m

PSO算法

%W为惯性权因子。
%测试函数选择y=0.5-(sin(sqrt(x1.^2+x2.^2)).^2-0.5)./(1+0.001*(x1.^2+x2.^2)).^2;
%此函数有无数个局部极大点，只有一个（0，0）为全局最大点，最大值为1，函数最大峰周围有两圈脊，取值分别为0.990284和0.962776,优化过程
%很容易停滞在这些局部极大点，x1和x2的取值范围都为[-100，100]；
%x1,x2初始位置为[-100，100]之间的随机取值
clear all
clc
tic  
   wmax=0.9;    %设置各项参数；
   wmin=0.3;
   itmax=500;   %粒子群迭代次数为40；
   c1=1.2;
   c2=1.2;
   V1max=100;
   V2max=100;
   Vmax=sqrt((V1max).^2+(V2max).^2);
   x1max=100;
   x2max=100;
   xmax=sqrt((x1max).^2+(x2max).^2);
   
for iter=1:itmax
    W(iter)=wmax-((wmax-wmin)/itmax)*iter; %设置惯性权因子W，使它的值逐渐减小。
end
%     rand('state',0);
    rand('state',sum(100*clock))
    a=-100;  
    b=200;
    N=40;  %粒子群大小为20个粒子；
    e1=rand(N,1,1);
    r1=rand(N,1,1);
for i=1:N
     x1(i,1)= a+b*e1(i);  %对各个粒子位置进随机初始化，Xi的取值范围是-100到100；
     x2(i,1)= a+b*r1(i);
     x(i,1)=x1(i,1)+j*x2(i,1); %第i个粒子的位置矢量和初始值
end 
    m=-100;
    n=200;
    e2=rand(N,1,1);
    r2=rand(N,1,1);
for i=1:N
     V1(i,1)=m+n*e2(i);
     V2(i,1)=m+n*r2(i); %对各个粒子速度进随机初始化，取值范围是-100到100；
     V(i,1)=V1(i,1)+j*V2(i,1); %第i个粒子的速度矢量和初始值
end
for i=1:N
    F(i,1)=0.5-(sin(sqrt(x1(i,1).^2+x2(i,1).^2)).^2-0.5)./(1+0.001*(x1(i,1).^2+x2(i,1).^2)).^2; %各个粒子初始函数值 ； 
end

   [C,I]=max(F(:,1)); %找出最大的那个粒子;
   gbest1(1)=x1(I,1);%初始化全局粒子最好位置；
   gbest2(1)=x2(I,1);
   gbest(1)=gbest1(1)+j*gbest2(1);
   Fbest(1)=0.5-(sin(sqrt(gbest1(1).^2+gbest2(1).^2)).^2-0.5)./(1+0.001*(gbest1(1).^2+gbest2(1).^2)).^2; %Fbest(1)为初始化时，全局最小函数值；
  for i=1:N
       pbest1(i,1)=x1(i,1); %设置第i个粒子的个体最好位置x1,x2的初始值；
       pbest2(i,1)=x2(i,1);
       pbest(i,1)=pbest1(i,1)+j*pbest2(i,1); %第i个粒子的二维矢量和初始值
       Fpbest(1)=0.5-(sin(sqrt(pbest1(i,1).^2+pbest2(i,1).^2)).^2-0.5)./(1+0.001*(pbest1(i,1).^2+pbest2(i,1).^2)).^2;
   end
%-------------------------------------------------------------------
   V1(:,2)=W(1)*V1(:,1)+c1*rand*(pbest1(:,1)-x1(:,1))+c2*rand*(gbest1(1)-x1(:,1)); %进化一次，即第二代粒子的速度和位置；
   V2(:,2)=W(1)*V1(:,1)+c1*rand*(pbest2(:,1)-x2(:,1))+c2*rand*(gbest2(1)-x2(:,1));
   x1(:,2)=x1(:,1)+V1(:,2); 
   x2(:,2)=x2(:,1)+V2(:,2);
   V(i,2)=V1(i,2)+j*V2(i,2); %第i个粒子的速度矢量和进化到第二代时的值；
   x(i,2)=x1(i,2)+j*x2(i,2); 
%限制x的范围在-xmax和xmax之间；----------------------------------；
if x1(i,2)>x1max             
   x1(i,2)=x1max;
else if x1(i,2)<-x1max
        x1(i,2)=-x1max;
   end
end
if x2(i,2)>x2max             
   x2(i,2)=x2max;
else if x2(i,2)<-x2max
        x2(i,2)=-x2max;
   end
end

%限制速度V在-Vmax和Vmax之间；--------------------------------------；

if V1(i,2)>V1max                 
   V1(i,2)=V1max;
else if V1(i,2)<-V1max
        V1(i,2)=-V1max;
   end
end
if V2(i,2)>V2max                 
   V2(i,2)=V2max;
else if V2(i,2)<-V2max
        V2(i,2)=-V2max;
   end
end
for k=2:itmax 
for i=1:N;
    F(i,k)=0.5-(sin(sqrt(x1(i,k).^2+x2(i,k).^2)).^2-0.5)./(1+0.001*(x1(i,k).^2+x2(i,k).^2)).^2;%计算全局粒子从第二代到第(itmax-1)=39代的函数值；
end
    [C,I]=max(F(:,k));     % 当进化到第k代时，找出使函数值为最大的那个粒子，为第I个粒子；
    gbest1(k)=x1(I,k);   % 进化到第k代时的全局粒子最好位置；
    gbest2(k)=x2(I,k);
    Fbest(k)=0.5-(sin(sqrt(gbest1(k).^2+gbest2(k).^2)).^2-0.5)./(1+0.001*(gbest1(k).^2+gbest2(k).^2)).^2; %第k代的全局粒子最大函数值；
    [D,P]=max(Fbest(:)); 
if  Fbest(k)>=D                 %使第k代的全局粒子最大函数值位置为gbest1(k),gbest2(k)
    gbest1(k)=gbest1(k);
    gbest2(k)=gbest2(k);   
else
     gbest1(k)=gbest1(P);
     gbest2(k)=gbest2(P);
end  
     gbest(k)=gbest1(k)+j*gbest2(k);
     Fbest(k)=0.5-(sin(sqrt(gbest1(k).^2+gbest2(k).^2)).^2-0.5)./(1+0.001*(gbest1(k).^2+gbest2(k).^2)).^2; %进化至第j代时的全局大函数值；
 for i=1:N;                  
     [C,Q]=max(F(i,:));
    if F(i,k)>=C
       pbest1(i,k)=x1(i,k);  
       pbest2(i,k)=x2(i,k);   
   else
       pbest1(i,k)=x1(i,Q);   %第i个粒子从第1代进化到第k代时的个体最好位置为pbest1(i,k);pbest2(i,k);
       pbest2(i,k)=x2(i,Q);
    end
       pbest(i,k)=pbest1(i,k)+j*pbest2(i,k);

       Fpbest(i,k)=0.5-(sin(sqrt(pbest1(i,k).^2+pbest2(i,k).^2)).^2-0.5)./(1+0.001*(pbest1(i,k).^2+pbest2(i,k).^2)).^2;
  end
    V1(:,k+1)=W(k)*V1(:,k)+c1*rand*(pbest1(:,k)-x1(:,k))+c2*rand*(gbest1(k)-x1(:,k)); %从第k代到第k+1代的速度v1更新；
    V2(:,k+1)=W(k)*V2(:,k)+c1*rand*(pbest2(:,k)-x2(:,k))+c2*rand*(gbest2(k)-x2(:,k)); %从第k代到第k+1代的速度v2更新
    x1(:,k+1)=x1(:,k)+V1(:,k+1);  %从第k代到第k+1代的位置x1更新；
    x2(:,k+1)=x2(:,k)+V2(:,k+1); 
    V(:,k+1)=V1(:,k+1)+j*V2(:,k+1);
    x(:,k+1)=x1(:,k+1)+j*x2(:,k+1);

%限制x的范围在-xmax和xmax之间；----------------------------------；
for i=1:N
if x1(i,k+1)>x1max             
   x1(i,k+1)=x1max;
else if x1(i,k+1)<-x1max
        x1(i,k+1)=-x1max;
   end
end
if x2(i,k+1)>x2max             
   x2(i,k+1)=x2max;
else if x2(i,k+1)<-x2max
        x2(i,k+1)=-x2max;
   end
end

%限制速度V在-Vmax和Vmax之间；--------------------------------------；

if V1(i,k+1)>V1max                 
   V1(i,k+1)=V1max;
else if V1(i,k+1)<-V1max
        V1(i,k+1)=-V1max;
   end
end

if V2(i,k+1)>V2max                 
   V2(i,k+1)=V2max;
else if V2(i,k+1)<-V2max
        V2(i,k+1)=-V2max;
   end
end
end
if (1-Fbest(k))<9.0*1e-3
  Fgoal=Fbest(k)
  break
end

   Fgoal=Fbest(k);

end
Fgoal
k


k=1:1:40;
 Gbest(k)=sqrt((gbest1(k)).^2+(gbest2(k)).^2);    %计算全局最好位置的矢量gbest(k)=gbest1(k)+j*gbest2(k)的绝对值Gbest(k).
 Pbest(4,k)=sqrt((pbest1(4,k)).^2+(pbest2(4,k)).^2);  %任意取4个粒子的个体最好位置的矢量pbest(4,k)=pbest1(4,k)+j*pbest2(4,k),求它的绝对值Pbest(4,k).
 X(4,k)=sqrt((x1(4,k)).^2+(x2(4,k)).^2);  %任意取4个粒子的位置矢量x(4,k)=x1(4,k)+j*x2(4,k),求它的绝对值X(4,k).

% scrsz=get(0,'ScreenSize');
% figure('Position',[1 scrsz(4)/5 scrsz(3)/1.1 4*scrsz(4)/5])
figure(1)

plot3(gbest1(k),gbest2(k),k,'r.')          % 画出进化至每一代时的全局粒子最好位置的实部(代指测试函数中的x1)和虚部(代指测试函数中的x2)的趋势图。实部和虚部的理论值都是0
hold on 
title('粒子的空间运行轨迹')

plot3(pbest1(4,k),pbest2(4,k),k,'b.')     %画出进化至每一代时第4个粒子个体最好位置的实部(代指测试函数中的x1)和虚部(代指测试函数中的x2)的趋势图。实部和虚部的理论值都是0
plot3(x1(4,k),x2(4,k),k,'k.')              %画出进化至每一代时第4个粒子个体位置的实部(代指测试函数中的x1)和虚部(代指测试函数中的x2)的趋势图。实部和虚部的理论值都是0
hold on
legend('全局最优粒子的空间位置','第4个粒子个体最优空间位置','第4个粒子的空间位置')
xlabel('自变量x1')
ylabel('自变量x2')
zlabel('进化代数')
grid on 
axis([-100 100 -100 100 1 40])
% figure('Position',[1 scrsz(4)/3 scrsz(3)/1.2 2*scrsz(4)/3])
figure(2)

plot(k,Fbest(k),'r.')  %画出全局最好位置的矢量gbest(k)=gbest1(k)+j*gbest2(k)的绝对值Gbest(k)的进化趋势。 理论值是0
hold on  
title('进化中粒子对应的函数值')

plot(k,Fpbest(4,k),'b.') % 画出第4个粒子的个体最好位置的矢量pbest(4,k)=pbest1(4,k)+j*pbest2(4,k)的绝对值Pbest(4,k)f进化趋势。理论值是0
plot(k,F(4,k),'k.')    % 画出第4个粒子的位置矢量x(4,k)=x1(4,k)+j*x2(4,k)的绝对值X(4,k)的进化趋势。理论值是0

legend('全局最优粒子对应的函数值','第4个粒子个体最优位置对应的函数值','第4个粒子对应的函数值')
xlabel('进化代数')
ylabel('函数值')
grid on 
disp('第4个粒子的矢量位置')
x(4,k)       % 第4个粒子的矢量位置
disp('第4个粒子的个体最好矢量位置')
pbest(4,k)   % 第4个粒子的个体最好矢量位置
disp('全局粒子的最好矢量位置')
gbest(k)     % 全局粒子的最好矢量位置
disp('函数极值,理论值为1.0')
Fbest(40)
toc