Code covered by the BSD License  

Highlights from
Particle Swarm Optimization Toolbox

image thumbnail

Particle Swarm Optimization Toolbox

by

 

22 Apr 2005 (Updated )

With Trelea, Common, and Clerc types along with ...

[OUT,varargout]=pso_Trelea_vectorized(functname,D,varargin)
% pso_Trelea_vectorized.m
% a generic particle swarm optimizer
% to find the minimum or maximum of any 
% MISO matlab function
%
% Implements Common, Trelea type 1 and 2, and Clerc's class 1". It will
% also automatically try to track to a changing environment (with varied
% success - BKB 3/18/05)
%
% This vectorized version removes the for loop associated with particle
% number. It also *requires* that the cost function have a single input
% that represents all dimensions of search (i.e., for a function that has 2
% inputs then make a wrapper that passes a matrix of ps x 2 as a single
% variable)
%
% Usage:
%  [optOUT]=PSO(functname,D)
% or:
%  [optOUT,tr,te]=...
%        PSO(functname,D,mv,VarRange,minmax,PSOparams,plotfcn,PSOseedValue)
%
% Inputs:
%    functname - string of matlab function to optimize
%    D - # of inputs to the function (dimension of problem)
%    
% Optional Inputs:
%    mv - max particle velocity, either a scalar or a vector of length D
%           (this allows each component to have it's own max velocity), 
%           default = 4, set if not input or input as NaN
%
%    VarRange - matrix of ranges for each input variable, 
%      default -100 to 100, of form:
%       [ min1 max1 
%         min2 max2
%            ...
%         minD maxD ]
%
%    minmax = 0, funct minimized (default)
%           = 1, funct maximized
%           = 2, funct is targeted to P(12) (minimizes distance to errgoal)
%
%    PSOparams - PSO parameters
%      P(1) - Epochs between updating display, default = 100. if 0, 
%             no display
%      P(2) - Maximum number of iterations (epochs) to train, default = 2000.
%      P(3) - population size, default = 24
%
%      P(4) - acceleration const 1 (local best influence), default = 2
%      P(5) - acceleration const 2 (global best influence), default = 2
%      P(6) - Initial inertia weight, default = 0.9
%      P(7) - Final inertia weight, default = 0.4
%      P(8) - Epoch when inertial weight at final value, default = 1500
%      P(9)- minimum global error gradient, 
%                 if abs(Gbest(i+1)-Gbest(i)) < gradient over 
%                 certain length of epochs, terminate run, default = 1e-25
%      P(10)- epochs before error gradient criterion terminates run, 
%                 default = 150, if the SSE does not change over 250 epochs
%                               then exit
%      P(11)- error goal, if NaN then unconstrained min or max, default=NaN
%      P(12)- type flag (which kind of PSO to use)
%                 0 = Common PSO w/intertia (default)
%                 1,2 = Trelea types 1,2
%                 3   = Clerc's Constricted PSO, Type 1"
%      P(13)- PSOseed, default=0
%               = 0 for initial positions all random
%               = 1 for initial particles as user input
%
%    plotfcn - optional name of plotting function, default 'goplotpso',
%              make your own and put here
%
%    PSOseedValue - initial particle position, depends on P(13), must be
%                   set if P(13) is 1 or 2, not used for P(13)=0, needs to
%                   be nXm where n<=ps, and m<=D
%                   If n<ps and/or m<D then remaining values are set random
%                   on Varrange
% Outputs:
%    optOUT - optimal inputs and associated min/max output of function, of form:
%        [ bestin1
%          bestin2
%            ...
%          bestinD
%          bestOUT ]
%
% Optional Outputs:
%    tr    - Gbest at every iteration, traces flight of swarm
%    te    - epochs to train, returned as a vector 1:endepoch
%
% Example:  out=pso_Trelea_vectorized('f6',2)

% Brian Birge
% Rev 3.3
% 2/18/06

function [OUT,varargout]=pso_Trelea_vectorized(functname,D,varargin)

rand('state',sum(100*clock));
if nargin < 2
   error('Not enough arguments.');
end

% PSO PARAMETERS
if nargin == 2      % only specified functname and D
   VRmin=ones(D,1)*-100; 
   VRmax=ones(D,1)*100;    
   VR=[VRmin,VRmax];
   minmax = 0;
   P = [];
   mv = 4;
   plotfcn='goplotpso';   
elseif nargin == 3  % specified functname, D, and mv
   VRmin=ones(D,1)*-100; 
   VRmax=ones(D,1)*100;    
   VR=[VRmin,VRmax];
   minmax = 0;
   mv=varargin{1};
   if isnan(mv)
       mv=4;
   end
   P = [];
   plotfcn='goplotpso';   
elseif nargin == 4  % specified functname, D, mv, Varrange
   mv=varargin{1};
   if isnan(mv)
       mv=4;
   end
   VR=varargin{2}; 
   minmax = 0;
   P = [];
   plotfcn='goplotpso';   
elseif nargin == 5  % Functname, D, mv, Varrange, and minmax
   mv=varargin{1};
   if isnan(mv)
       mv=4;
   end    
   VR=varargin{2};
   minmax=varargin{3};
   P = [];
   plotfcn='goplotpso';
elseif nargin == 6  % Functname, D, mv, Varrange, minmax, and psoparams
   mv=varargin{1};
   if isnan(mv)
       mv=4;
   end    
   VR=varargin{2};
   minmax=varargin{3};
   P = varargin{4}; % psoparams
   plotfcn='goplotpso';   
elseif nargin == 7  % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn
   mv=varargin{1};
   if isnan(mv)
       mv=4;
   end    
   VR=varargin{2};
   minmax=varargin{3};
   P = varargin{4}; % psoparams
   plotfcn = varargin{5}; 
elseif nargin == 8  % Functname, D, mv, Varrange, minmax, and psoparams, plotfcn, PSOseedValue
   mv=varargin{1};
   if isnan(mv)
       mv=4;
   end    
   VR=varargin{2};
   minmax=varargin{3};
   P = varargin{4}; % psoparams
   plotfcn = varargin{5};  
   PSOseedValue = varargin{6};
else    
   error('Wrong # of input arguments.');
end

% sets up default pso params
Pdef = [100 2000 24 2 2 0.9 0.4 1500 1e-25 250 NaN 0 0];
Plen = length(P);
P    = [P,Pdef(Plen+1:end)];

df      = P(1);
me      = P(2);
ps      = P(3);
ac1     = P(4);
ac2     = P(5);
iw1     = P(6);
iw2     = P(7);
iwe     = P(8);
ergrd   = P(9);
ergrdep = P(10);
errgoal = P(11);
trelea  = P(12);
PSOseed = P(13);

% used with trainpso, for neural net training
if strcmp(functname,'pso_neteval')
   net = evalin('caller','net');
    Pd = evalin('caller','Pd');
    Tl = evalin('caller','Tl');
    Ai = evalin('caller','Ai');
     Q = evalin('caller','Q');
    TS = evalin('caller','TS');
end


% error checking
 if ((minmax==2) & isnan(errgoal))
     error('minmax= 2, errgoal= NaN: choose an error goal or set minmax to 0 or 1');
 end

 if ( (PSOseed==1) & ~exist('PSOseedValue') )
     error('PSOseed flag set but no PSOseedValue was input');
 end

 if exist('PSOseedValue')
     tmpsz=size(PSOseedValue);
     if D < tmpsz(2)
         error('PSOseedValue column size must be D or less');
     end
     if ps < tmpsz(1)
         error('PSOseedValue row length must be # of particles or less');
     end
 end
 
% set plotting flag
if (P(1))~=0
  plotflg=1;
else
  plotflg=0;
end

% preallocate variables for speed up
 tr = ones(1,me)*NaN;

% take care of setting max velocity and position params here
if length(mv)==1
 velmaskmin = -mv*ones(ps,D);     % min vel, psXD matrix
 velmaskmax = mv*ones(ps,D);      % max vel
elseif length(mv)==D     
 velmaskmin = repmat(forcerow(-mv),ps,1); % min vel
 velmaskmax = repmat(forcerow( mv),ps,1); % max vel
else
 error('Max vel must be either a scalar or same length as prob dimension D');
end
posmaskmin  = repmat(VR(1:D,1)',ps,1);  % min pos, psXD matrix
posmaskmax  = repmat(VR(1:D,2)',ps,1);  % max pos
posmaskmeth = 3; % 3=bounce method (see comments below inside epoch loop)

% PLOTTING
 message = sprintf('PSO: %%g/%g iterations, GBest = %%20.20g.\n',me);
 
% INITIALIZE INITIALIZE INITIALIZE INITIALIZE INITIALIZE INITIALIZE
 
% initialize population of particles and their velocities at time zero,
% format of pos= (particle#, dimension)
 % construct random population positions bounded by VR
  pos(1:ps,1:D) = normmat(rand([ps,D]),VR',1);
  
  if PSOseed == 1         % initial positions user input, see comments above
    tmpsz                      = size(PSOseedValue);
    pos(1:tmpsz(1),1:tmpsz(2)) = PSOseedValue;  
  end

 % construct initial random velocities between -mv,mv
  vel(1:ps,1:D) = normmat(rand([ps,D]),...
      [forcecol(-mv),forcecol(mv)]',1);

% initial pbest positions vals
 pbest = pos;

% VECTORIZE THIS, or at least vectorize cost funct call 
 out = feval(functname,pos);  % returns column of cost values (1 for each particle)
%---------------------------
 
 pbestval=out;   % initially, pbest is same as pos

% assign initial gbest here also (gbest and gbestval)
 if minmax==1
   % this picks gbestval when we want to maximize the function
    [gbestval,idx1] = max(pbestval);
 elseif minmax==0
   % this works for straight minimization
    [gbestval,idx1] = min(pbestval);
 elseif minmax==2
   % this works when you know target but not direction you need to go
   % good for a cost function that returns distance to target that can be either
   % negative or positive (direction info)
    [temp,idx1] = min((pbestval-ones(size(pbestval))*errgoal).^2);
    gbestval    = pbestval(idx1);
 end

 % preallocate a variable to keep track of gbest for all iters
 bestpos        = zeros(me,D+1)*NaN;
 gbest          = pbest(idx1,:);  % this is gbest position
   % used with trainpso, for neural net training
   % assign gbest to net at each iteration, these interim assignments
   % are for plotting mostly
    if strcmp(functname,'pso_neteval')
        net=setx(net,gbest);
    end
 %tr(1)          = gbestval;       % save for output
 bestpos(1,1:D) = gbest;
 
% this part used for implementing Carlisle and Dozier's APSO idea
% slightly modified, this tracks the global best as the sentry whereas
% their's chooses a different point to act as sentry
% see "Tracking Changing Extremea with Adaptive Particle Swarm Optimizer",
% part of the WAC 2002 Proceedings, June 9-13, http://wacong.com
 sentryval = gbestval;
 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
rstflg = 0; % for dynamic environment checking
% 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
 iwt(1) = iw1;
for i=1:me  % start epoch loop (iterations)

     out        = feval(functname,[pos;gbest]);
     outbestval = out(end,:);
     out        = out(1:end-1,:);

     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];
     
     %assignin('base','bestpos',bestpos(i,1:D+1));
   %------------------------------------------------------------------------      
   % 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

    % 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; % for dynamic environment checking

    if chkdyn==1
     threshld = 0.05;  % percent current best is allowed to change, .05 = 5% etc
     letiter  = 5; % # of iterations before checking environment, leave at least 3 so PSO has time to converge
     outorng  = abs( 1- (outbestval/gbestval) ) >= threshld;
     samepos  = (max( sentry == gbest ));

     if (outorng & samepos) & rem(i,letiter)==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; % reset personal bests to current positions
        pbestval  = out; 
        vel       = vel*10; % agitate particles a little (or a lot)
        
       % recalculate best vals 
        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,:);
        
        % used with trainpso, for neural net training
        % assign gbest to net at each iteration, these interim assignments
        % are for plotting mostly
        if strcmp(functname,'pso_neteval')
           net=setx(net,gbest);
        end
     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
     %[size(out),size(pbestval)]
    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,:);
            % used with trainpso, for neural net training
            % assign gbest to net at each iteration, these interim assignments
            % are for plotting mostly
             if strcmp(functname,'pso_neteval')
                net=setx(net,gbest);
             end
        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,:);
            % used with trainpso, for neural net training
            % assign gbest to net at each iteration, these interim assignments
            % are for plotting mostly
             if strcmp(functname,'pso_neteval')
                net=setx(net,gbest);
             end
        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,:);
           % used with trainpso, for neural net training
            % assign gbest to net at each iteration, these interim assignments
            % are for plotting mostly
             if strcmp(functname,'pso_neteval')
                net=setx(net,gbest);
             end
        end
     end
    end
    
    
 %   % 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
% 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 by at least ',...
              num2str(ergrd),' 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

 %    % convert back to inertial frame
 %     pos = pos - repmat(gbestoffset,ps,1);
 %     pbest = pbest - repmat(gbestoffset,ps,1);
 %     gbest = gbest + gbestoffset;
  

end  % end epoch loop

%% clear temp outputs
% evalin('base','clear temp_pso_out temp_te temp_tr;');

% output & return
 OUT=[gbest';gbestval];
 varargout{1}=[1:te];
 varargout{2}=[tr(find(~isnan(tr)))];
 
 return

Contact us