实现粒子群算法的MATLAB代码 - EDA365电子论坛网

function [x, fval, exitflag, output] =pso(objfunc, nvars, options)
% PSO Particle Swarm Optimization Algorithm
%
% PSO attempts to solve continuousoptimization problem of the form
%
% min OBJFUNC(X)
% X
%
% No constraints on the values of entriesof X can be imposed, not ewen box contraints in form of
% the lower and upper bound. Any desiredconstraints should be dealth with by modification of the
% objective function: by introduction ofpenalty functions, by reformulation of the problem, or some
% other way.
%
% In the present implementation, onlyclassical 'star topology' PSO is considered, where each
% particle is interconnected with eachother particle, and there are no subswarms. Also, no
% additional, GA like, operators areapplied to the particles during the search.
%
% PSO(OBJFUNC, NVARS) Optimizes objective OBJFUNC defined as a real-valuedfunction of a single,
% real-valued vector argument. Dimensionality of the search space, that isthe number of entries
% in the argument of OBJFUNC, is defined in NVARS. PSO expects at leastthese two arguments.
% Failure to provide any of them will result in error.
%
% PSO(OBJFUNC, NVARS, OPTIONS) Enables the user to specify a number ofsettings in order to
% customize or fine-tune the performance of the algorithm. OPTIONS is astructure containing
% desired values of these settings. Detailed description of availableparameters is given in the
% sequel.
%
% X= PSO(...) Returns vector at which objective attains its minimal value. Due tothe nature of
% PSO algorithm, this is not guaranteed to be neither exact global norlocal optimum, yet PSO is
% robust optimizer, well suited for complex, multimodal problems. Welltuned, it can give VERY
% GOOD and quite commonly EXCELENT solution to the problem at hand. Inmost cases, this is all
% that is needed.
%
% [X, FVAL] = PSO(...) In addition to the optimal X, it returns the valueof the objective at X,
% FVAL = OBJFUNC(X).
%
% [X, FVAL, EXITFLAG] = PSO(...) Returns indication concerning a reasonwhy the algorithm stopped.
% In the current implementation the only supported return values are 0 and1. EXITFLAG 0 denotes
% that maximum number of iterations has been achieved, while EXITFLAG 1 isused in testing mode,
% if the value of global minimum has been found prior to achieving themaximal number of
% iteration.
%
% [X, FVAL, EXITFLAG, OUTPUT] = PSO(...) Returns OUTPUT structurecontaining valuable information
% concerning performance of the algorithm in the present run. The exactmembers of this structure
% are variable, and depend on the settings specified in the OPTIONSstructure. In general, there
% are four levels of detail that can be specified. The 'LOW' detail levelmeans that the algorithm
% keeps track of objective value for the best, the mean and the worstparticle in the swarm in
% each generation. If the 'MEDIUM' level is specified, exact position andindex of the global best
% particle are tracked for during each iteration. If the 'HIGH' level isspecified, exact position
% of each particle in the swarm is logged during the entire searchprocess. If the swarm is large
% and the number of generations is great, high level output logging canresult in considerable
% memmory consumption. Finally, the output log level can be set to 'NONE',meaning that no data is
% wanted within the OUTPUT structure. The exact way for specifying outputlogging level is given
% below.
%
% OPTIONS = PSO('options') Returns default options structure. Usefull whenone desires to change
% values of only a handfull of options.
%
% The OPTIONS structure contains the following entries
%
% options.npart The numberof particles in the swarm.
% options.niter Themaximal number of iterations.
% options.cbi Initialvalue of the individual-best acceleration factor.
% options.cbf Finalvalue of the individual-best acceleration factor.
% options.cgi Initialvalue of the global-best acceleration factor.
% options.cgf Finalvalue of the global-best acceleration factor.
% options.wi Initialvalue of the inertia factor.
% options.wf Finalvalue of the inertia factor.
% options.vmax Absolutespeed limit. If specified, the speed is clamped to the range
% [-options.vmax,options.vmax]. It is the primary speed limit, if set
% to NaN, theVMAXSCALE options is used.
% options.vmaxscale Relativespeed limit. Used only if absolute limit is unspecified, i.e.
% set to NaN. Ifused, must be a scalar quantity, and denotes the amount
% of initialpopulation span (the INITSPAN option) used as the speed
% limit.
% options.vspaninit Theinitial velocity span. Initial velocities are initialized
% uniformly in[-VSPANINIT, VSPANINIT]. This option must be specified.
% options.initoffset Offset ofthe initial population. Can be scalar or column-vector of
% dimension NVARS.
% options.initspan Span ofthe initial population. Can be scalar or column-vector of
% dimension NVARS.
% options.trustoffset If set to1 (true) and offset is vector, than the offset is
% believed to be agood solution candidate, so it is included in
% the initial swarm.
% options.initpopulation Theuser-suplied initial population. If this is set to something
% meaningfull, thenINITSPAN, INITOFFSET and TRUSTOFFSET are
% ignored. If set to NaN then the abovementioned offset is used.
% options.verbose_period Theverbose period, i.e. the number of iterations after which the
% results areprompted to the user. If set to 0, then verbosing is
% skipped.
% options.plot If set to1, evolution of the global best is ploted to the user after
% the optimizationprocess. The objective value of the best, mean
% and worseparticle troughout the optimization process are plotted
% in the singlegraph.
% options.output_level The outputlog level. Possible values are: 'none', 'low',
% 'medium', 'high'.Each log level denotes a specific amount of
% information to bereturned to the end user. If less than 4 output
% arguments arespecified log level is ignored, since it would only
% occupate (possibly)large amount of useless memory.
% options.globalmin Globalminimum, used for testing only. If specified, the algorithm
% will stop as soonas the difference between GLOBALMIN option and
% current globalbest becomes less than TOL option.
% options.tol Precisiontolerance, used for testing only. It is maximal difference
% between currentglobal best and GLOBALMIN option at which the
% algorithm stops.If GLOBALMIN options is set to NaN this option is
% ignored, and thealgorithm stops after the maximal number of iteration
% has been achieved.
%
% The OUTPUT structure is the fallowing
%
% output.itersno The actualnumber of iterations. Usualy the same as NITER options, but
% can be less if intesting mode (if GLOBALMIN option is specified).
%
% If at least LOW level output logging is specified:
%
% output.gbest_array
% output.gmean_array
% output.gworst_array Arrayscontaining the objective value of the best, mean and worst
% particle in eachiteration. In fact, gmean_array does not contain
% objective valuefor any concrete particle, but instead it contains the
% mean objectivevalue for the entire swarm in each iteration.
%
% If at least MEDIUM level output logging isspecified:
%
% output.gbes**x_array The arrayontaining indices of the best particle in each iteration.
% output.Xbest Matrix ofdimension NITERxNVARS, containing, as rows, global best
% particles in eachiteration.
%
% Only if HIGH level output logging is specified:
%
% output.X 3D matrixof dimension NPARTxNVARSxNITER containing the entire
% population in eachiteration.
%
%
% The following examples ilustrate the use of PSO function in severalcommon cases.
%
% Suppose we are attempting to optimize 5-dimensional objective f(x) =x*x'. Since, it is assumed
% that objective receives a row-vector, the objective is in fact a 5Dparaboliod. The easiest way
% to optimize would be to write
%
% obj = @(x) x*x';
% [xopt, fopt] = pso(obj, 5);
%
% The preseding code lines would yield XOPT as the found near-optimalsolution and FOPT as the
% objective value in such point. Of course, other outputs can be obtainedas described earlier.
%
% If one should choose to change the default options (say to specifydifferent number of particles
% and iterations), the code would look something like this
%
% obj = @(x) x*x';
% options = pso('options');
% options.niter = 500;
% options.npart = 60;
% [xopt, fopt] = pso(obj, 5);
%
% Other options can be specified as well on the exactly the same way. Thebest way to explore the
% option structure is to experiment.
%
% PSO is compatible with MATLAB 7 and higher. Some modifications areneeded for it to work under
% lower versions of MATLAB.
% Copyright (C) Milan Rapaic(repmilan@yahoo.com), 2007-05-10 ... 2008-10-12
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the number of input and outputarguments. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
msg = nargchk(1, 3, nargin);
if ~isempty(msg)
error('mrr:myoptim:pso:pso:narginerr', 'Inadequate number of inputarguments.');
end
msg = nargchk(0, 4, nargout);
if ~isempty(msg)
error('mrr:myoptim:pso:pso:nargouterr', 'Inadequate number of outputarguments.');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The body of the algorithm. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin==1 && ischar(objfunc)&& strcmp(objfunc, 'options')
%User desired only to access the default OPTIONS structure.
if nargout<=1
x = getDefaultOptions();
else
% The user required multiple outputs, yet only default options can bereturned.
error('mrr:myoptim:pso:pso:nargouterr', ...
'Cannot expext more than one output when only OPTIONS are required.');
end
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%User requested optimization to be conducted on a given objective. %
%The following code deals with initializations and validations of optionsstructure. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%If no options are specified, use the default ones.
if nargin<3, options=getDefaultOptions(); end
%Determination of output level, that is of amount of data to be collected inOUTPUT structure.
if nargout == 4
% User supplied four output arguments, therefore output level isdetermined from the OPTIONS
% structure. The output_level variable is used to code the desired loglevel: 0 ('none'), 1
% ('low'), 2 ('medium') and 3 ('high).
if strcmp(options.output_level, 'none')
if options.plot == 0
output_level = 0;
else
output_level = 1;
end
elseif strcmp(options.output_level, 'low')
output_level = 1;
elseif strcmp(options.output_level, 'medium')
output_level = 2;
elseif strcmp(options.output_level, 'high')
output_level = 3;
else
error('mrr:myoptim: pso:pso: optionserr: output_level', ...
'Invalid value of theOUTPUT_LEVEL options specified.');
end
else
% User has not supplied forth output argument. The only reason to loginformation during the
% run is to be able to plot to to the user after the optimizationprocess. Therefore, if
% ploting is requested low level logging is used.
if options.plot == 1
output_level = 1;
else
output_level = 0;
end
end
%Maximum velocity can be specified in absolute amount, or relative to theinitial span.
%If both values are specified, the absolute velocity limit is taken intoaccount, while the
%relative is ignored. Whatever the initial specification of the maximalvelocity, the curent
%code block will generate a column vector, vmax, containing maximal velocityalong each
%dimension of search space.
if ~all(isnan(options.vmax))
% It is not allowed to let some of the entries of the VMAX option to beNaN and others to
% have numerical or Inf values.
if any(isnan(options.vmax))
error('mrr:myoptim: pso: pso: optionserr:vmax', ...
'VMAX option cannot have someInf and some numerical (or Inf) values.');
end
% Warning of the confusing entries within the OPTIONS structure.
if ~isnan(options.vmaxscale)
warning('mrr:myoptim: pso:pso: optionserr:vmaxconflict', ...
'Both relative andabsolute velocity limit are specified. The relative limit is ignored.');
end
if length(options.vmax) == 1
vmax = options.vmax*ones(nvars, 1);
elseif length(options.vmax) == nvars
% Maximal velocity should be a column-vector or a scalar.
if size(options.vmax, 1) ~= length(options.vmax)
error('mrr:myopim: pso:pso: optionserr:vmax', ...
'VMAX option should bespecified as column-vector, or as a scalar value.');
end
vmax = options.vmax;
else
error('mrr:myoptim: pso:pso: optionserr:vmax', ...
'Inadequate dimension of VMAXoption. Should be a scalar, or a column vector with NVARS elements.');
end
else
% It is not valid to specify both VMAX and VMAXSCALE option as NaN.
if isnan(options.vmaxscale)
error('mrr:myoptim:pso:pso: optionserr:vmaxscale', ...
'Either VMAX or VMAXSCALEoptions should be different than NaN.');
end
% Contrary to the VMAX options, VMAXSCALE option must allways be ascalar. The initial span
% should take into account the different scaling among the cooedinatesof the search space.
if length(options.vmaxscale) == 1
if length(options.initspan) == 1
vmax =options.vmaxscale*options.initspan*ones(nvars, 1);
else
% If the dimension of INITSPANoption is not correct, the function will break later,
% therefore, no need to checkvalidity now.
vmax =options.vmaxscale*options.initspan;
end
else
error('mrr:myoptim:pso:pso: optionserr:vmax', ...
'Inadequate dimension ofVMAXSCALE option. Must be a scalar.');
end
end
vmax = repmat(vmax', options.npart, 1);
%Initial population.
%If the initial population is not supplied by the user, each particle of theinitial population
%is spred in [INITOFFSET-INITSPAN, INITOFFSET+INITSPAN] where both INITOFFSETand INITSPAN
%are specified within the OPTIONS structure. Both of these options are eitherscalars or
%column-vectors of appropriate size. If INITPOPULATION option is specified, bothINITOFFSET and
%INITSPAN options are ignored.
if ~isnan(options.initpopulation)
% The user supplied complete initial population within the OPTIONSstructure.
% The size of the supplied population must be consistent with populationsize and number of
% variables. If no, an error is reported.
[pno, pdim] = size(options.initpopulation);
if (pno ~= options.npart) || (pdim ~= nvars)
error('mrr:myoptim: pso: pso: optionserr:initpopulation', ...
['The format of initialpopulation is inconsistent with desired population', ...
'size or dimension of searchspace - INITPOPULATION options is invalid']);
end
X = options.initpopulation;
elseif (length(options.initoffset) == 1) &&(length(options.initspan) == 1)
% The same offset and span is specified for each dimension of the searchspace
X = (rand(options.npart, nvars)-0.5)*2*options.initspan +options.initoffset;
elseif (length(options.initoffset) ~= size(options.initoffset, 1)) ||...
(length(options.initspan) ~= size(options.initspan, 1))
error('mrr:myoptim:pso:pso: optionserr:initoffset_initspan', ...
'Both INITOFFSET and INITSPAN options must be either scalars orcolumn-vectors.');
elseif (length(options.initoffset) ~= nvars) ||(length(options.initspan) ~= nvars)
error('mrr:myoptim: pso: pso: optionserr:init', ...
'Both INITOFFSET and INITSPAN options must be scalars or column-vectorsof length NVARS.');
else
initoffset = repmat(options.initoffset', options.npart, 1);
initspan =repmat(options.initspan', options.npart, 1);
X = (rand(options.npart, nvars)-0.5)*2.*initspan + initoffset;
% TRUSTOFFSET option is used when OFFSET option is, in fact, previouslyknown good (or very
% good) solution to the problem at hand. When set to logical true (1),offset is inserted in
% the initial population. Thus, it is guaranteed that objective value atsolution is not
% greater than objective value at that, previously known, good point.
if (options.trustoffset)
X(1, : ) = options.initoffset';
end
end
%Initial velocities.
%Velocities are initialized uniformly in [-VSPANINIT, VSPANINIT].
if any(isnan(options.vspaninit))
error('mrr:myoptim: pso:pso: optionserr:vspaninit', ...
'VSPANINIT option must notcontain NaN entries.');
elseif isscalar(options.vspaninit)
V = (rand(options.npart, nvars)-0.5)*2*options.vspaninit;
else
if (length(options.vspaninit) ~= size(options.vspaninit, 1)) || ...
(length(options.vspaninit) ~= nvars)
error('mrr:myoptim: pso:pso: optionserr:vspaninit', ...
'VSPANINIT option must beeither scalar or column-vector of length NVARS');
end
V = (rand(options.npart,nvars)-0.5)*2.*repmat(options.vspaninit', options.npart, 1);
end
%Initial scores (objective values).
%Initialization of the best personal score and position, as well as global bestscore and
%position.
Y= calcobjfunc(objfunc, X);
Ybest = Y; %The best individual score for each particle - initialization.
Xbest = X; %The best individual position for each particle -
%initialization.
[GYbest, gbest] = min(Ybest); %GYbest is the best score within the entire swarm.
% gbest isthe index of particle that achived YGbest.
gbest = gbest(1); %In case when more than one particle achieved the best
% score, wechoose the one with the lowest index as the
% best one.
%These variables are used in testing mode only.
tolbreak = ~isnan(options.globalmin);
foundglobal = 0;
if tolbreak && ~isscalar(options.globalmin)
error('mrr:myoptim:pso:pso: optionserr:globalmin', ...
'globalmin option, if specified, option must be a scalar value equal tothe global minimum of the objective function');
end
%Initialization of the OUTPUT structure.
%The output structure is filled and initialized differently depending on theOUTPUT_LEVEL
%options, or equivalently depending on the output_level variable.
if output_level >= 0
% NONE log level
output.itersno = options.niter;
if output_level >= 1
% LOW log level
output.gbest_array = NaN*ones(options.niter+1, 1);
output.gmean_array = NaN*ones(options.niter+1, 1);
output.gworst_array = NaN*ones(options.niter+1, 1);
output.gbest_array(1) = GYbest;
output.gmean_array(1) = mean(Ybest);
output.gworst_array(1) = max(Ybest);
if output_level >= 2
% MEDIUM log level
output.gbes**x_array =NaN*ones(options.niter+1, 1);
output.Xbest =NaN*ones(options.niter+1, nvars);
output.gbes**x_array(1) =gbest;
output.Xbest(1, : ) = X(gbest,: );
if output_level == 3
% HIGH log level
output.X =NaN*zeros(options.npart, nvars, options.niter+1);
output.X(:,:,1) = X;
end
end
end
end
if options.verbose_period ~= 0
disp 'PSO algorithm: Initiating the optimization process.'
end
%Denotes normal algorithm termination.
exitflag = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The main loop. %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for iter = 1: options.niter
% Verbosing, if neccessary.
if options.verbose_period ~= 0
if rem(iter, options.verbose_period) == 0
disp(['iteration ', int2str(iter),'. best criteria = ', num2str(GYbest)]);
end
end
% Calculating PSO parameters
w = linrate(options.wf, options.wi, options.niter, 0, iter);
cp = linrate(options.cbf, options.cbi, options.niter, 0, iter);
cg = linrate(options.cgf, options.cgi, options.niter, 0, iter);
% For later calculations only
GXbest = repmat(Xbest(gbest, : ), options.npart, 1);
% Calculating speeds
V = w*V + cp*rand(size(V)).*(Xbest-X) + cg*rand(size(V)).*(GXbest-X);
V = min(vmax, abs(V)).*sign(V);
% Population is moving
X = X + V;
Y = calcobjfunc(objfunc, X);
% Calculating new individually best values
mask = Y<Ybest;
mask = repmat(mask, 1, nvars);
Xbest = mask.*X +(~mask).*Xbest;
Ybest = min(Y,Ybest);
% Calculating new globally best value
[GYbest, gbest] = min(Ybest);
gbest = gbest(1);
% Filling in the OUTPUT structure.
if output_level >= 0
% NONE log level
if output_level >= 1
% LOW log level
output.gbest_array(iter+1) = GYbest;
output.gmean_array(iter+1) = mean(Ybest);
output.gworst_array(iter+1) =max(Ybest);
if output_level >= 2
% MEDIUM log level
output.gbes**x_array(iter+1)= gbest;
output.Xbest(iter+1, : ) =X(gbest, : );
if output_level == 3
% HIGH log level
output.X(:,:,iter+1) =X;
end
end
end
end
% The code used in testing mode only.
if tolbreak && abs(GYbest - options.globalmin)<options.tol
output.itersno = iter;
foundglobal = 1;
break
end
end
if options.verbose_period ~= 0
disp 'Optimization process finished.'
end
%Setting up the output variables.
%X is set to be the final global best position, since that is the best positionever achieved
%by any of the particles in the swarm. FVAL is set to be the value of theobjective in X.
x= Xbest(gbest, : ); x = x(: );
fval = GYbest;
%The global moptimum has been found prior to achieving the maximal number ofiteration.
if foundglobal, exitflag = 1; end;
%Plotting the algorithm behavior at each iteration.
if options.plot
r = 0: options.niter;
figure
plot(r, output.gbest_array, 'k.', r, output.gmean_array, 'r.', r,output.gworst_array, 'b.');
str = sprintf('Best objective value : %g', fval);
title(str);
legend({'best objective', 'mean objective', 'worst objective'})
end
end
function Y = calcobjfunc(func, X)
% CALCOBJFUNC A helper function used to calculateobjective function value for a series of points.
np = size(X,1);
Y = zeros(np,1);
for i = 1:np
Y(i) = func(X(i,: ));
end
function opts = getDefaultOptions
% GETDEFAULTOPTIONS Returns a structure containing the defaultoptions.
%
% This function, in fact, defines default values of the options within theoptions structure.
opts.npart = 30; % The number of particles.
opts.niter = 100; % The number of iterations.
opts.cbi = 2.5; % Initial value of the individual-bestacceleration factor.
opts.cbf = 0.5; % Final value of the individual-bestacceleration factor.
opts.cgi = 0.5; % Initial value of the global-bestacceleration factor.
opts.cgf = 2.5; % Final value of the global-bestacceleration factor.
opts.wi = 0.9; % Initial value of the inertia factor.
opts.wf = 0.4; % Final value of the inertia factor.
opts.vmax = Inf; % Absolute speed limit. It is the primaryspeed limit.
opts.vmaxscale = NaN; % Relative speed limit. Used only ifabsolute limit is unspecified.
opts.vspaninit = 1; % The initial velocity span. Initialvelocities are initialized
% uniformly in[-VSPANINIT, VSPANINIT].
opts.initoffset = 0; % Offset of the initial population.
opts.initspan = 1; % Span of the initial population.
opts.trustoffset = 0; % If set to 1 (true) and offset is vector, than the offset is
% believed tobe a good solution candidate, so it is included in
% the initialswarm.
opts.initpopulation = NaN; % The user-suplied initial population. Ifthis is set to something
% meaningfull,then INITSPAN, INITOFFSET and TRUSTOFFSET are
% ignored.
opts.verbose_period = 10; % The verbose period, i.e. the number ofiterations after which the
% results areprompted to the user. If set to 0, then verbosing is
% skipped.
opts.plot = 0; % If set to 1, evolution of the gbestis ploted to the user after
% theoptimization process. The objective value of the best, mean
% and worseparticle troughout the optimization process are plotted
% in the singlegraph.
opts.output_level = 'low'; % The output log level. Possible values are: 'none', 'low',
% 'medium','high'.
opts.globalmin = NaN; % Global minimum, used for testing only
opts.tol = 1e-6; % Precision tolerance, used for testingonly
function x = linrate(xmax, xmin, tmax,tmin, t)
% LINRATE Linear interpolation of value X in instant T, defined by previouslyknown points
% (tmin, xmin), (tmax, xmax)
x = xmin + ((xmax-xmin)/(tmax-tmin))*(tmax-t);' F6 A/ i4 x) t6 A6 T6 S9 B