% 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/06function [OUT,varargout]=pso_Trelea_vectorized(functname,D,varargin)rand('state',sum(100*clock));if nargin < 2   error('Not enough arguments.');end% PSO PARAMETERSif 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 paramsPdef = [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 trainingif 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 flagif (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 hereif length(mv)==1 velmaskmin = -mv*ones(ps,D);     % min vel, psXD matrix velmaskmax = mv*ones(ps,D);      % max velelseif length(mv)==D      velmaskmin = repmat(forcerow(-mv),ps,1); % min vel velmaskmax = repmat(forcerow( mv),ps,1); % max velelse error('Max vel must be either a scalar or same length as prob dimension D');endposmaskmin  = repmat(VR(1:D,1)',ps,1);  % min pos, psXD matrixposmaskmax  = repmat(VR(1:D,2)',ps,1);  % max posposmaskmeth = 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