tvpso.m

时变参数的粒子群优化的matlab实现

function [x, fval, exitflag, output] = pso(objfunc, nvars, options)
% PSO   Particle Swarm Optimization Algorithm
%
% PSO attempts to solve continuous optimization problem of the form
%
%       min OBJFUNC(X)
%        X
%
% No constraints on the values of entries of X can be imposed, not ewen box contraints in form of
% the lower and upper bound. Any desired constraints should be dealth with by modification of the 
% objective function: by introduction of penalty functions, by reformulation of the problem, or some
% other way.
%
% In the present implementation, only classical 'star topology' PSO is considered, where each
% particle is interconnected with each other particle, and there are no subswarms. Also, no
% additional, GA like, operators are applied to the particles during the search.
%
%   PSO(OBJFUNC, NVARS) Optimizes objective OBJFUNC defined as a real-valued function of a single,
%   real-valued vector argument. Dimensionality of the search space, that is the number of entries
%   in the argument of OBJFUNC, is defined in NVARS. PSO expects at least these two arguments.
%   Failure to provide any of them will result in error.
%
%   PSO(OBJFUNC, NVARS, OPTIONS) Enables the user to specify a number of settings in order to
%   customize or fine-tune the performance of the algorithm. OPTIONS is a structure containing
%   desired values of these settings. Detailed description of available parameters is given in the
%   sequel.
%
%   X = PSO(...) Returns vector at which objective attains its minimal value. Due to the nature of
%   PSO algorithm, this is not guaranteed to be neither exact global nor local optimum, yet PSO is
%   robust optimizer, well suited for complex, multimodal problems. Well tuned, it can give VERY 
%   GOOD and quite commonly EXCELENT solution to the problem at hand. In most cases, this is all 
%   that is needed.
%
%   [X, FVAL] = PSO(...) In addition to the optimal X, it returns the value of the objective at X,
%   FVAL = OBJFUNC(X).
%
%   [X, FVAL, EXITFLAG] = PSO(...) Returns indication concerning a reason why the algorithm stopped.
%   In the current implementation the only supported return values are 0 and 1. EXITFLAG 0 denotes
%   that maximum number of iterations has been achieved, while EXITFLAG 1 is used in testing mode,
%   if the value of global minimum has been found prior to achieving the maximal number of
%   iteration.
%
%   [X, FVAL, EXITFLAG, OUTPUT] = PSO(...) Returns OUTPUT structure containing valuable information
%   concerning performance of the algorithm in the present run. The exact members of this structure
%   are variable, and depend on the settings specified in the OPTIONS structure. In general, there
%   are four levels of detail that can be specified. The 'LOW' detail level means that the algorithm
%   keeps track of objective value for the best, the mean and the worst particle in the swarm in
%   each generation. If the 'MEDIUM' level is specified, exact position and index of the global best
%   particle are tracked for during each iteration. If the 'HIGH' level is specified, exact position
%   of each particle in the swarm is logged during the entire search process. If the swarm is large
%   and the number of generations is great, high level output logging can result 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 output logging level is given
%   below.
%
%   OPTIONS = PSO('options') Returns default options structure. Usefull when one desires to change
%   values of only a handfull of options.
%
%   The OPTIONS structure contains the following entries
%
%     options.npart           The number of particles in the swarm.
%     options.niter           The maximal number of iterations.
%     options.cbi             Initial value of the individual-best acceleration factor.
%     options.cbf             Final value of the individual-best acceleration factor.
%     options.cgi             Initial value of the global-best acceleration factor.
%     options.cgf             Final value of the global-best acceleration factor.
%     options.wi              Initial value of the inertia factor.
%     options.wf              Final value of the inertia factor.
%     options.vmax            Absolute speed limit. If specified, the speed is clamped to the range
%                             [-options.vmax, options.vmax]. It is the primary speed limit, if set
%                             to NaN, the VMAXSCALE options is used.
%     options.vmaxscale       Relative speed limit. Used only if absolute limit is unspecified, i.e.
%                             set to NaN. If used, must be a scalar quantity, and denotes the amount
%                             of initial population span (the INITSPAN option) used as the speed
%                             limit.
%     options.vspaninit       The initial velocity span. Initial velocities are initialized 
%                             uniformly in [-VSPANINIT, VSPANINIT]. This option must be specified.
%     options.initoffset      Offset of the initial population. Can be scalar or column-vector of
%                             dimension NVARS.
%     options.initspan        Span of the initial population. Can be scalar or column-vector of
%                             dimension NVARS.
%     options.trustoffset     If set to 1 (true) and offset is vector, than the offset is 
%                             believed to be a good solution candidate, so it is included in 
%                             the initial swarm.
%     options.initpopulation  The user-suplied initial population. If this is set to something
%                             meaningfull, then INITSPAN, INITOFFSET and TRUSTOFFSET are
%                             ignored. If set to NaN then the above mentioned offset is used.
%     options.verbose_period  The verbose period, i.e. the number of iterations after which the 
%                             results are prompted to the user. If set to 0, then verbosing is
%                             skipped.
%     options.plot            If set to 1, evolution of the global best is ploted to the user after
%                             the optimization process. The objective value of the best, mean
%                             and worse particle troughout the optimization process are plotted
%                             in the single graph.
%     options.output_level    The output log level. Possible values are: 'none', 'low', 
%                             'medium', 'high'. Each log level denotes a specific amount of
%                             information to be returned to the end user. If less than 4 output
%                             arguments are specified log level is ignored, since it would only
%                             occupate (possibly) large amount of useless memory.
%     options.globalmin       Global minimum, used for testing only. If specified, the algorithm
%                             will stop as soon as the difference between GLOBALMIN option and
%                             current global best becomes less than TOL option.
%     options.tol             Precision tolerance, used for testing only. It is maximal difference
%                             between current global best and GLOBALMIN option at which the
%                             algorithm stops. If GLOBALMIN options is set to NaN this option is
%                             ignored, and the algorithm stops after the maximal number of iteration
%                             has been achieved.
%
%   The OUTPUT structure is the fallowing
%
%     output.itersno          The actual number of iterations. Usualy the same as NITER options, but
%                             can be less if in testing mode (if GLOBALMIN option is specified).
%
%     If at least LOW level output logging is specified:
%
%     output.gbest_array
%     output.gmean_array
%     output.gworst_array     Arrays containing the objective value of the best, mean and worst
%                             particle in each iteration. In fact, gmean_array does not contain
%                             objective value for any concrete particle, but instead it contains the
%                             mean objective value for the entire swarm in each iteration.
%
%     If at least MEDIUM level output logging is specified:
%
%     output.gbestndx_array   The array ontaining indices of the best particle in each iteration.
%     output.Xbest            Matrix of dimension NITERxNVARS, containing, as rows, global best 
%                             particles in each iteration.
%
%     Only if HIGH level output logging is specified:
%
%     output.X                3D matrix of dimension NPARTxNVARSxNITER containing the entire
%                             population in each iteration.
%
%
%   The following examples ilustrate the use of PSO function in several common 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 5D paraboliod. 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-optimal solution and FOPT as the
%   objective value in such point. Of course, other outputs can be obtained as described earlier.
%
%   If one should choose to change the default options (say to specify different 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. The best way to explore the
%   option structure is to experiment.
%
%   PSO is compatible with MATLAB 7 and higher. Some modifications are needed for it to work under
%   lower versions of MATLAB.

% Copyright (C) Milan Rapaic ([e-mail]), 2007-05-10 ... 2008-10-12


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Checking the number of input and output arguments.                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

msg = nargchk(1, 3, nargin);
if ~isempty(msg) 
    error('mrr:myoptim:pso:pso:narginerr', 'Inadequate number of input arguments.'); 
end

msg = nargchk(0, 4, nargout);
if ~isempty(msg)
    error('mrr:myoptim:pso:pso:nargouterr', 'Inadequate number of output arguments.');
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 be returned.
        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 options structure.          %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % 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 in OUTPUT structure.
    if nargout == 4 
        % User supplied four output arguments, therefore output level is determined from the OPTIONS
        % structure. The output_level variable is used to code the desired log level: 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 the OUTPUT_LEVEL options specified.');
        end
    else
        % User has not supplied forth output argument. The only reason to log information during the
        % run is to be able to plot to to the user after the optimization process. 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 the initial span.
    % If both values are specified, the absolute velocity limit is taken into account, while the 
    % relative is ignored. Whatever the initial specification of the maximal velocity, the curent
    % code block will generate a column vector, vmax, containing maximal velocity along 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 be NaN and others to
        % have numerical or Inf values.
        if any(isnan(options.vmax))
            error('mrr:myoptim:pso:pso:optionserr:vmax', ...
                'VMAX option cannot have some Inf 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 and absolute 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 be specified as column-vector, or as a scalar value.');
            end
            vmax = options.vmax;
        else
            error('mrr:myoptim:pso:pso:optionserr:vmax', ...
                'Inadequate dimension of VMAX option. 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 VMAXSCALE options should be different than NaN.');
        end
        % Contrary to the VMAX options, VMAXSCALE option must allways be a scalar. The initial span
        % should take into account the different scaling among the cooedinates of the search space.