pso.m

matlab环境下

 
function pso      
global x;
x(1,:)=[5.1,4.9,4.7,4.6,5.0,5.4,4.6,5.0,4.4,4.9,5.4,4.8,4.8,4.3,5.8,5.7,5.4,5.1,5.7,5.1,5.4,5.1,4.6,5.1,4.8,5.0,5.0,5.2,5.2,4.7,4.8,5.4,5.2,5.5,4.9,5.0,5.5,4.9,4.4,5.1,5.0,4.5,4.4,5.0,5.1,4.8,5.1,4.6,5.3,5.0];
x(2,:)=[3.5,3.0,3.2,3.1,3.6,3.9,3.4,3.4,2.9,3.1,3.7,3.4,3.0,3.0,4.0,4.4,3.9,3.5,3.8,3.8,3.4,3.7,3.6,3.3,3.4,3.0,3.4,3.5,3.4,3.2,3.1,3.4,4.1,4.2,3.1,3.2,3.5,3.1,3.0,3.4,3.5,2.3,3.2,3.5,3.8,3.0,3.8,3.2,3.7,3.3];
x(3,:)=[1.4,1.4,1.3,1.5,1.4,1.7,1.4,1.5,1.4,1.5,1.5,1.6,1.4,1.1,1.2,1.5,1.3,1.4,1.7,1.5,1.7,1.5,1.0,1.7,1.9,1.6,1.6,1.5,1.4,1.6,1.6,1.5,1.5,1.4,1.5,1.2,1.3,1.5,1.3,1.5,1.3,1.3,1.3,1.6,1.9,1.4,1.6,1.4,1.5,1.4];
x(4,:)=[0.2,0.2,0.2,0.2,0.2,0.4,0.3,0.2,0.2,0.1,0.2,0.2,0.1,0.1,0.2,0.4,0.4,0.3,0.3,0.3,0.2,0.4,0.2,0.5,0.2,0.2,0.4,0.2,0.2,0.2,0.2,0.4,0.1,0.2,0.1,0.2,0.2,0.1,0.2,0.2,0.3,0.3,0.2,0.6,0.4,0.3,0.2,0.2,0.2,0.2];
success=0;
PopSize=20;
MaxIt=70;
iter=1;
fevals=0;%function evaluation'counter 
maxw=1.2;%maximum interia weight's value
minw=0.1;%minimum interia weight's value
c1=0.5;
c2=0.5;
inertdec=(maxw-minw)/MaxIt;% weight decrease by same value, search change function of weight
w=maxw;
f='gold';
dim=4;
ErrGoal=1e-4;
%initializing swarm and velocity
vel=rand(dim,PopSize);

for i=1:4
   popul(i,:)=rand(1,PopSize)*(max(x(i,:))-min(x(i,:)))+min(x(i,:));
end
vel=rand(dim,PopSize);
%evaluation initial population
for i=1:PopSize
    x1=popul(1,i); %popul is 2D array,2D particle population
     x2=popul(2,i);  
     x3=popul(3,i);
     x4=popul(4,i);
     fpopul(i)=feval(f,x1,x2,x3,x4);   %fpopul is 1D array,function values of 2D particle population
     fevals=fevals+1;
     ppopul(:,i,iter)=popul(:,i);
end
bestpos=popul; % positions of  initial particle
fbestpos=fpopul;% func values of initial particle
%ppopul(:,:,iter)=ppopul;
% finding best particle in initial 
[fbestpart(iter),g]=min(fpopul);% min value of population and index.
lastbpf=fbestpart(iter);% current values of optimal particle

%swarm evolution loop 
while(iter<MaxIt)
    iter=iter+1;
    % update the value of inertia weight w
    if(iter<=MaxIt)
        w=maxw-(iter-1)*inertdec;
    end
    
    %velocity update
    for i=1:PopSize
        A(:,i)=bestpos(:,g);% optical position of each particle is g.
    end    
        R1=rand(dim,PopSize);
        R2=rand(dim,PopSize);
        vel=w*vel+c1*R1.*(bestpos-popul)+c2*R2.*(A-popul);%local optical is bestpos,global optical is A.
        %swarm update
        popul=popul+vel; 
        for i=1:4
            for j=1:PopSize
                 ppopul(i,j,iter)=popul(i,j);
             end
         end
            %evaluate the new swarm
        for i=1:PopSize
            x1=popul(1,i);
            x2=popul(2,i);
            x3=popul(3,i);
            x4=popul(4,i);
            fpopul(iter,i)=feval(f,x1,x2,x3,x4);
            fevals=fevals+1;
        end
        
        %update the best position for each particle
        for i=1:PopSize 
          if(fpopul(iter,i)<fbestpos(i))
           fbestpos(i)=fpopul(iter,i);
            bestpos(:,i)=popul(:,i);
        end
        end
              
        % updateing index g
        [fbestpart(iter),g]=min(fbestpos);
        %checking stoping criterion
       if (abs(fbestpart(iter)-fbestpart(iter-1))<=ErrGoal)&(abs(fpopul(iter,1)-fpopul(iter-1,1))<=ErrGoal)
           break;
       end
    end
    
    %swarm evoluation loop
    % output arguments 
    xmin=popul(:,g)
    fxmin=fbestpos(g)
    bestpos
   % plot3(bestpos(1,:),bestpos(2,:),bestpos(3,:),'b.')%bestpos
   subplot(2,2,1)
    a=iter/3;
   
    for i=1:a
       i=i*3;
       % text(xx1(i,1),yy1(i,1),zz1(i,1),['\leftarrow(',num2str(i),')'],'FontSize',6.0)
       plot3(ppopul(1,1,i),ppopul(2,1,i),ppopul(3,1,i),'linestyle','--','marker','.','markersize',3) 
       xlabel('萼片长度/cm','fontsize',8),ylabel('花瓣高度/cm','fontsize',8),zlabel('花瓣宽度/cm','fontsize',8)
       title('粒子1优化过程轨迹','FontSize',8)
       axis([2.5  7 2 5 0 4])
      % text(ppopul(1,1,i)+0.1,ppopul(2,1,i),ppopul(3,1,i),[,num2str(i)],'FontSize',7,'linestyle','-')
       hold on
       %'\leftarrow'
    end
    
    iter
    fevals
   % ppopul(:,:,10)
   % ppopul(:,:,20)
      

  function sum=gold(x1,x2,x3,x4)
   global x;
  sum=0;
     x(1,:)=[5.1,4.9,4.7,4.6,5.0,5.4,4.6,5.0,4.4,4.9,5.4,4.8,4.8,4.3,5.8,5.7,5.4,5.1,5.7,5.1,5.4,5.1,4.6,5.1,4.8,5.0,5.0,5.2,5.2,4.7,4.8,5.4,5.2,5.5,4.9,5.0,5.5,4.9,4.4,5.1,5.0,4.5,4.4,5.0,5.1,4.8,5.1,4.6,5.3,5.0];
     x(2,:)=[3.5,3.0,3.2,3.1,3.6,3.9,3.4,3.4,2.9,3.1,3.7,3.4,3.0,3.0,4.0,4.4,3.9,3.5,3.8,3.8,3.4,3.7,3.6,3.3,3.4,3.0,3.4,3.5,3.4,3.2,3.1,3.4,4.1,4.2,3.1,3.2,3.5,3.1,3.0,3.4,3.5,2.3,3.2,3.5,3.8,3.0,3.8,3.2,3.7,3.3];
     x(3,:)=[1.4,1.4,1.3,1.5,1.4,1.7,1.4,1.5,1.4,1.5,1.5,1.6,1.4,1.1,1.2,1.5,1.3,1.4,1.7,1.5,1.7,1.5,1.0,1.7,1.9,1.6,1.6,1.5,1.4,1.6,1.6,1.5,1.5,1.4,1.5,1.2,1.3,1.5,1.3,1.5,1.3,1.3,1.3,1.6,1.9,1.4,1.6,1.4,1.5,1.4];
     x(4,:)=[0.2,0.2,0.2,0.2,0.2,0.4,0.3,0.2,0.2,0.1,0.2,0.2,0.1,0.1,0.2,0.4,0.4,0.3,0.3,0.3,0.2,0.4,0.2,0.5,0.2,0.2,0.4,0.2,0.2,0.2,0.2,0.4,0.1,0.2,0.1,0.2,0.2,0.1,0.2,0.2,0.3,0.3,0.2,0.6,0.4,0.3,0.2,0.2,0.2,0.2];
     for i=1:50
     %sum=sum+0.148*abs(x11(i)-x1)+0.097*abs(x12(i)-x2)+0.346*abs(x13(i)-x3)+0.408*abs(x14(i)-x4);    		
      sum=sum+0.148*(x(1,i)-x1)^2+0.097*(x(2,i)-x2)^2+0.346*(x(3,i)-x3)^2+0.408*(x(4,i)-x4)^2;
     end