www.pudn.com > 1111111registrationbasedoncorrelation.rar > PSO.m

function [OUT,varargout]=PSO(structure) 
D=3; 
rand('state',sum(100*clock)); 
if nargin < 1 
   error('Not enough arguments.'); 
end 
 
     
% PSO PARAMETERS 
   VRmin=ones(D,1)*-20;  
   VRmax=ones(D,1)*20;     
   VR=[VRmin,VRmax]; 
   minmax = 1; 
 
P =[1 2000 20 4 2 2 0.9 0.2 1500 2 1e-5 20 1]; 
df  = P(1); 
me  = P(2); 
ps  = P(3); 
mv  = P(4); 
ac1 = P(5); 
ac2 = P(6); 
iw1 = P(7); 
iw2 = P(8); 
iwe = P(9); 
flagg=P(10); 
ergrd=P(11); 
ergrdep=P(12); 
plotflg=P(13); 
 
% PLOTTING 
 message = sprintf('PSO: %%g/%g iterations, GBest = %%g.\n',me); 
 
pos=40*rand(ps,D)-20; 
vel=8*rand(ps,D)-4; 
 
  
% initial pbest positions vals 
 pbest=pos; 
  
 for j=1:ps  % start particle loop 
    numin='0'; 
    for i=1:D 
        numin=strcat(numin,',',num2str(pos(j,i))); 
    end 
%     evstrg=strcat('feval(''',functname,'''',numin(2:end),',structure',')');   
evstrg=strcat('feval(''myMI''',numin(2:end),',structure',')'); 
    out(j)=eval(evstrg);     % evaluate desired function with particle j       
 end 
  
 pbestval=out;   % initially, pbest is same as pos 
 
% assign initial gbest here also (gbest and gbestval) 
 if minmax==1 
    [gbestval,idx1]=max(pbestval);  % this picks gbestval when we want to maximize the function 
 elseif minmax==0 
    [gbestval,idx1]=min(pbestval);  % this works for straight minimization 
 end 
 gbest=pbest(idx1,:);  % this is gbest position 
 tr(1)=gbestval;       % save for output 
 
% start PSO iterative procedures 
cnt=0; % counter used for updating display according to df in the options 
cnt2=0; % counter used for the stopping subroutine based on error convergence 
 
 
for i=1:me  % start epoch loop (iterations) 
   if flagg==0   % randimization control, one random set for each epoch 
       rannum1=rand(1);   
       rannum2=rand(2); 
   end 
    
   for j=1:ps  % start particle loop 
        
     if flagg==1   % randomization control, one random set for each particle at each epoch 
         rannum1=rand(1); 
         rannum2=rand(1); 
     end 
 
     numin='0'; 
     for dimcnt=1:D 
         numin=strcat(numin,',',num2str(pos(j,dimcnt))); 
     end 
%      evstrg=strcat('feval(''',functname,'''',numin(2:end),',structure',')');  
evstrg=strcat('feval(''myMI''',numin(2:end),',structure',')'); 
     out(j)=eval(evstrg);     % evaluate desired function with particle j   
     e(j) = out(j);              % use to minimize or maximize function to unknown values 
 
     %SSEhist(j) = sumsqr(e);    % sum squared 'error' for jth particle (averages if there is more than one output) 
      
    % update pbest to reflect whether searching for max or min of function 
     if minmax==0 
       if pbestval(j)>=e(j); 
          pbestval(j)=e(j); 
          pbest(j,:)=pos(j,:); 
       end 
     elseif minmax==1 
       if pbestval(j)<=e(j); 
           pbestval(j)=e(j); 
           pbest(j,:)=pos(j,:); 
       end 
     end 
  
           
    % assign gbest by finding minimum of all particle pbests  
     if minmax==1 
       [iterbestval,idx1]=max(pbestval);  % this picks gbestval when we want to maximize the function 
       if gbestval<=iterbestval 
          gbestval=iterbestval; 
          gbest=pbest(idx1,:); 
       end        
     elseif minmax==0  
       [iterbestval,idx1]=min(pbestval);  % this works for straight minimization and for minimizing error to target 
       if gbestval>=iterbestval 
          gbestval=iterbestval; 
          gbest=pbest(idx1,:); 
       end        
     end     
     
     tr(i+1)=gbestval; % keep track of global best val 
     te=i;           % this will return the epoch number to calling program when done 
 
    % get new velocities, positions (this is the heart of the PSO algorithm) 
     if i<=iwe 
        iwt(i)=((iw2-iw1)/(iwe-1))*(i-1)+iw1; % get inertia weight, just a linear funct w.r.t. epoch parameter iwe 
     else 
        iwt(i)=iw2; 
     end 
     
     if flagg==2              % each component of each particle gets a different random number set 
        for dimcnt=1:D 
           rannum1=rand(1); 
           rannum2=rand(1); 
           vel(j,dimcnt)=iwt(i)*vel(j,dimcnt)... 
                    +ac1*rannum1*(pbest(j,dimcnt)-pos(j,dimcnt))... 
                    +ac2*rannum2*(gbest(1,dimcnt)-pos(j,dimcnt)); 
        end 
     else                     % this velocity update is for flagg= 0 or 1 
        vel(j,:)=iwt(i)*vel(j,:)... 
                 +ac1*rannum1*(pbest(j,:)-pos(j,:))... 
                 +ac2*rannum2*(gbest(1,:)-pos(j,:)); 
     end 
 
    % update new position 
     pos(j,:)=pos(j,:)+vel(j,:);     
 
    % limit velocity/position components to maximums 
     for dimcnt=1:D 
        if vel(j,dimcnt)>mv 
           vel(j,dimcnt)=mv; 
        end 
        
        if vel(j,dimcnt)<-mv 
          vel(j,dimcnt)=-mv; 
        end     
         
        if pos(j,dimcnt)>=VR(dimcnt,2) 
           pos(j,dimcnt)=VR(dimcnt,2); 
        end 
        
        if pos(j,dimcnt)<=VR(dimcnt,1) 
           pos(j,dimcnt)=VR(dimcnt,1); 
        end        
         
         
     end 
        
   end         % end particle loop 
    
  %----------------------------------------------------------------------------------------------------------------------  
    
    
   % check for stopping criterion based on speed of convergence to desired error    
    tmp1=abs(tr(i)-gbestval); 
    if tmp1>ergrd 
       cnt2=0; 
    elseif tmp1<=ergrd 
       cnt2=cnt2+1; 
       if cnt2>=ergrdep 
         if plotflg==1 
          fprintf(message,i,gbestval);            
%           disp(' '); 
%           disp('***** global error gradient too small for too long'); 
%           disp('***** this means you''ve got the solution or it got stuck'); 
         end 
         break 
       end        
    end      
end  % end epoch loop 
      OUT=[gbest';gbestval]; 
 OUT(1:3)=round(OUT(1:3)); 
 varargout{1}=[0:te]; 
 varargout{2}=[tr];