PSO算法优化MCKD

% 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')%tf = strcmp（s1，s2）比较s1和s2，如果二者相同则返回1（真），否则返回0（假）。 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;%me为迭代次数 % 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 =（粒子号，维数）的格式构造由VR界定的随机群体位置 pos(1:ps,1:D) = normmat(rand([ps,D]),VR',1); if PSOseed == 1 % initial positions user input, see comments above tmpsz