pso_trelea_vectorized.m

C++实现的PSO算法,很有用写关于PSO算法的论文少不了他

 sentry    = gbest;
 
if (trelea == 3)
% calculate Clerc's constriction coefficient chi to use in his form
 kappa   = 1; % standard val = 1, change for more or less constriction    
 if ( (ac1+ac2) <=4 )
     chi = kappa;
 else
     psi     = ac1 + ac2;
     chi_den = abs(2-psi-sqrt(psi^2 - 4*psi));
     chi_num = 2*kappa;
     chi     = chi_num/chi_den;
 end
end

% INITIALIZE END INITIALIZE END INITIALIZE END INITIALIZE END

% 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)
     t0  = clock;
     out = feval(functname,[pos;gbest]);
     outbestval = out(end,:);
     out = out(1:end-1,:);

    % check for an error space that changes wrt time/iter
    % threshold value that determines dynamic environment 
    % sees if the value of gbest changes more than some threshold value
    % for the same location
    chkdyn=1;
    rstflg=0;    
    
    if chkdyn==1
     threshld = .20;  % +/- 20% of current best
     outorng  = abs( 1- (outbestval/gbestval) ) >= threshld;
     samepos  = (max( sentry == gbest ));

     if (outorng & samepos) & rem(i,5)==0
         rstflg=1;
       % disp('New Environment: reset pbest, gbest, and vel');
       %% reset pbest and pbestval if warranted
%        outpbestval = feval( functname,[pbest] );
%        Poutorng    = abs( 1-(outpbestval./pbestval) ) > threshld;
%        pbestval    = pbestval.*~Poutorng + outpbestval.*Poutorng;
%        pbest       = pbest.*repmat(~Poutorng,1,D) + pos.*repmat(Poutorng,1,D);   

        pbest     = pos;
        pbestval  = out;
        vel       = vel*1.2;
        
        if minmax == 1
           [gbestval,idx1] = max(pbestval);
        elseif minmax==0
           [gbestval,idx1] = min(pbestval);
        elseif minmax==2 % this section needs work
           [temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2);
           gbestval    = pbestval(idx1);
        end
        
        gbest  = pbest(idx1,:);
     end  % end if outorng
     
     sentryval = gbestval;
     sentry    = gbest;
     
    end % end if chkdyn
    
    % find particles where we have new pbest, depending on minmax choice 
    % then find gbest and gbestval
    if rstflg==0
     if minmax==0
        [tempi]=find(pbestval>=out);         % new min pbestvals
        pbestval(tempi,1)=out(tempi);        % update pbestvals
        pbest(tempi,:)=pos(tempi,:);         % update pbest positions
       
        [iterbestval,idx1]=min(pbestval);
        if gbestval>=iterbestval
            gbestval=iterbestval;
            gbest=pbest(idx1,:);
        end
     elseif minmax==1
        [tempi,dum]=find(pbestval<=out);     % new max pbestvals
        pbestval(tempi,1)=out(tempi,1);      % update pbestvals
        pbest(tempi,:)=pos(tempi,:);         % update pbest positions
 
        [iterbestval,idx1]=max(pbestval);
        if gbestval<=iterbestval
            gbestval=iterbestval;
            gbest=pbest(idx1,:);
        end
     elseif minmax==2  % this won't work as it is, fix it later
        egones=errgoal*ones(ps,1);           % vector of errgoals
        sqrerr2=((pbestval-egones).^2);
        sqrerr1=((out-egones).^2);
        [tempi,dum]=find(sqerr1 <= sqrerr2); % find particles closest to targ
        pbestval(tempi,1)=out(tempi,1);      % update pbestvals
        pbest(tempi,:)=pos(tempi,:);         % update pbest positions

        sqrerr=((pbestval-egones).^2);    % need to do this to reflect new pbests
        [temp,idx1]=min(sqrerr);
        iterbestval=pbestval(idx1);
        if (iterbestval-errgoal)^2 <= (gbestval-errgoal)^2
           gbestval=iterbestval;
           gbest=pbest(idx1,:);
        end
     end
    end
    tr(i+1) = gbestval; % keep track of global best val
    te      = i;        % returns epoch number to calling program when done
    bestpos(i,1:D+1) = [gbest,gbestval];
    
 %   % build a simple predictor 10th order, for gbest trajectory
 %   if i>500
 %    for dimcnt=1:D
 %      pred_coef  = polyfit(i-250:i,(bestpos(i-250:i,dimcnt))',20);
 %     % pred_coef  = polyfit(200:i,(bestpos(200:i,dimcnt))',20);       
 %      gbest_pred(i,dimcnt) = polyval(pred_coef,i+1);
 %    end
 %    else 
%       gbest_pred(i,:) = zeros(size(gbest));
%    end
  
   gbest_pred(i,:)=gbest;    
   assignin('base','gbest_pred',gbest_pred);

 %   % convert to non-inertial frame
 %    gbestoffset = gbest - gbest_pred(i,:);
 %    gbest = gbest - gbestoffset;
 %    pos   = pos + repmat(gbestoffset,ps,1);
 %    pbest = pbest + repmat(gbestoffset,ps,1);     
 
     %PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO

      % get new velocities, positions (this is the heart of the PSO algorithm)     
      % each epoch get new set of random numbers
       rannum1 = rand([ps,D]); % for Trelea and Clerc types
       rannum2 = rand([ps,D]);       
       if     trelea == 2    
        % from Trelea's paper, parameter set 2
         vel = 0.729.*vel...                              % prev vel
               +1.494.*rannum1.*(pbest-pos)...            % independent
               +1.494.*rannum2.*(repmat(gbest,ps,1)-pos); % social  
       elseif trelea == 1
        % from Trelea's paper, parameter set 1                     
         vel = 0.600.*vel...                              % prev vel
               +1.700.*rannum1.*(pbest-pos)...            % independent
               +1.700.*rannum2.*(repmat(gbest,ps,1)-pos); % social 
       elseif trelea ==3
        % Clerc's Type 1" PSO
         vel = chi*(vel...                                % prev vel
               +ac1.*rannum1.*(pbest-pos)...              % independent
               +ac2.*rannum2.*(repmat(gbest,ps,1)-pos)) ; % social          
       else
        % common PSO algo with inertia wt 
        % get inertia weight, just a linear funct w.r.t. epoch parameter iwe
         if i<=iwe
            iwt(i) = ((iw2-iw1)/(iwe-1))*(i-1)+iw1; 
         else
            iwt(i) = iw2;
         end
        % random number including acceleration constants
         ac11 = rannum1.*ac1;    % for common PSO w/inertia
         ac22 = rannum2.*ac2;
         
         vel = iwt(i).*vel...                             % prev vel
               +ac11.*(pbest-pos)...                      % independent
               +ac22.*(repmat(gbest,ps,1)-pos);           % social                  
       end
       
       % limit velocities here using masking
        vel = ( (vel <= velmaskmin).*velmaskmin ) + ( (vel > velmaskmin).*vel );
        vel = ( (vel >= velmaskmax).*velmaskmax ) + ( (vel < velmaskmax).*vel );     
        
       % update new position (PSO algo)    
        pos = pos + vel;
    
       % position masking, limits positions to desired search space
       % method: 0) no position limiting, 1) saturation at limit,
       %         2) wraparound at limit , 3) bounce off limit
        minposmask_throwaway = pos <= posmaskmin;  % these are psXD matrices
        minposmask_keep      = pos >  posmaskmin;     
        maxposmask_throwaway = pos >= posmaskmax;
        maxposmask_keep      = pos <  posmaskmax;
     
        if     posmaskmeth == 1
         % this is the saturation method
          pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos );
          pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos );      
        elseif posmaskmeth == 2
         % this is the wraparound method
          pos = ( minposmask_throwaway.*posmaskmax ) + ( minposmask_keep.*pos );
          pos = ( maxposmask_throwaway.*posmaskmin ) + ( maxposmask_keep.*pos );                
        elseif posmaskmeth == 3
         % this is the bounce method, particles bounce off the boundaries with -vel      
          pos = ( minposmask_throwaway.*posmaskmin ) + ( minposmask_keep.*pos );
          pos = ( maxposmask_throwaway.*posmaskmax ) + ( maxposmask_keep.*pos );

          vel = (vel.*minposmask_keep) + (-vel.*minposmask_throwaway);
          vel = (vel.*maxposmask_keep) + (-vel.*maxposmask_throwaway);
        else
         % no change, this is the original Eberhart, Kennedy method, 
         % it lets the particles grow beyond bounds if psoparams (P)
         % especially Vmax, aren't set correctly, see the literature
        end

     %PSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSOPSO

 %    % convert back to inertial frame
 %     pos = pos - repmat(gbestoffset,ps,1);
 %     pbest = pbest - repmat(gbestoffset,ps,1);
 %     gbest = gbest + gbestoffset;
  %------------------------------------------------------------------------      
   % 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(['--> Solution likely, GBest hasn''t changed for ',...
                  num2str(cnt2),' epochs.']);  
          eval(plotfcn);
         end       
         break
       end
    end
   % this stops if using constrained optimization and goal is reached
    if ~isnan(errgoal)
     if ((gbestval<=errgoal) & (minmax==0)) | ((gbestval>=errgoal) & (minmax==1))  

         if plotflg==1
             fprintf(message,i,gbestval);
           disp(' ');            
             disp(['--> Error Goal reached, successful termination!']);
             
             eval(plotfcn);
         end
         break
     end
    % this is stopping criterion for constrained from both sides    
     if minmax==2
       if ((tr(i)<errgoal) & (gbestval>=errgoal)) | ((tr(i)>errgoal) ...
               & (gbestval<=errgoal))        
         if plotflg==1
             fprintf(message,i,gbestval);
           disp(' ');            
             disp(['--> Error Goal reached, successful termination!']);            
             
             eval(plotfcn);
         end
         break              
       end
     end % end if minmax==2
    end  % end ~isnan if
   
% this section does the plots during iterations   
if plotflg==1      
  if (rem(i,df) == 0 ) | (i==me) | (i==1) 
     fprintf(message,i,gbestval);
     cnt=cnt+1; % count how many times we display (useful for movies)
      
       eval(plotfcn); % defined at top of script
     
  end  % end update display every df if statement    
 end % end plotflg if statement
 
 psotime(i)=etime(clock,t0);
 assignin('base','psotime',psotime);
end  % end epoch loop

% clear temp outputs
 evalin('base','clear temp_pso_out temp_te temp_tr;');
% hold off
% outputs
 OUT=[gbest';gbestval];
 assignin('base','bestpos',[bestpos]);
 varargout{1}=[0:te];
 varargout{2}=[tr(find(~isnan(tr)))];