% 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 pbes
