function [xOpt,fval,exitflag,output,population,scores] = ...
pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB,nonlcon,options)
% Find the minimum of a function using Particle Swarm Optimization.
%
%
%
% Syntax:
% psodemo
% pso
% x = pso(fitnessfcn,nvars)
% x = pso(fitnessfcn,nvars,Aineq,bineq)
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq)
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB)
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB,nonlcon)
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB,nonlcon,options)
% x = pso(problem)
% [x, fval] = pso(...)
% [x, fval,exitflag] = pso(...)
% [x, fval,exitflag,output] = pso(...)
% [x, fval,exitflag,output,population] = pso(...)
% [x, fval,exitflag,output,population,scores] = pso(...)
%
% Description:
% psodemo
% Runs a demonstration of the PSO algorithm using test function specified
% by user.
%
% pso
% Runs a default demonstration using Rosenbrock's banana function.
%
% x = pso(fitnessfcn,nvars)
% Runs the particle swarm algorithm with no constraints and default
% options. fitnessfcn is a function handle, nvars is the number of
% parameters to be optimized, i.e. the dimensionality of the problem. x is
% a 1xnvars vector representing the coordinates of the global optimum
% found by the pso algorithm.
%
% x = pso(fitnessfcn,nvars,Aineq,bineq)
% Linear constraints, such that Aineq*x <= bineq. Aineq is a matrix of size
% nconstraints x nvars, while b is a column vector of length nvars.
%
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq)
% Linear equality constraints, such that Aeq*x = beq. 'Soft' or 'penalize'
% boundaries only. If no inequality constraints exist, set Aineq and bineq
% to [].
%
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB)
% Lower and upper bounds defined as LB and UB respectively. Both LB and UB,
% if defined, should be 1 x nvars vectors. Use empty arrays [] for Aineq,
% bineq, Aeq, or beq if no linear constraints exist.
%
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB,nonlcon)
% Non-linear constraints. Nonlinear inequality constraints in the form c(x)
% <= 0 have now been implemented using 'soft' boundaries, or
% experimentally, using 'penalize' constraint handling method. See the
% Optimization Toolbox documentation for the proper syntax for defining
% nonlinear constraints. Use empty arrays [] for Aineq, bineq, Aeq, beq,
% LB, or UB if they are not needed.
%
% x = pso(fitnessfcn,nvars,Aineq,bineq,Aeq,beq,LB,UB,nonlcon,options)
% Default algorithm parameters replaced with those defined in the options
% structure:
% Use >> options = psooptimset('Param1,'value1','Param2,'value2',...) to
% generate the options structure. Type >> psooptimset with no input or
% output arguments to display a list of available options and their
% default values.
%
% NOTE: the swarm will perform better if the PopInitRange option is defined
% so that it closely bounds the expected domain of the feasible region.
% This is automatically done if the problem has both lower and upper bound
% constraints, but not for linear and nonlinear constraints.
%
% NOTE 2: If options.HybridFcn is to be defined, and if it is necessary to
% pass a non-default options structure to the hybrid function, then the
% options structure may be included in a cell array along with the pointer
% to the hybrid function. Example:
% >> % Let's say that we want to use fmincon to refine the result from PSO:
% >> hybridoptions = optimset(@fmincon) ;
% >> options.HybridFcn = {@fmincon, hybridoptions} ;
%
% NOTE 3:
% Perez and Behdinan (2007) demonstrated that the particle swarm is only
% stable if the following conditions are satisfied:
% Given that C0 = particle inertia
% C1 = options.SocialAttraction
% C2 = options.CognitiveAttraction
% 1) 0 < (C1 + C2) < 4
% 2) (C1 + C2)/2 - 1 < C0 < 1
% If conditions 1 and 2 are satisfied, the system will be guaranteed to
% converge to a stable equilibrium point. However, whether or not this
% point is actually the global minimum cannot be guaranteed, and its
% acceptability as a solution should be verified by the user.
%
% x = pso(problem)
% Parameters imported from problem structure.
%
% [x,fval] = pso(...)
% Returns the fitness value of the best solution.
%
% [x,fval,exitflag] = pso(...)
% Returns information on outcome of pso run. This should match exitflag
% values for GA where applicable, for code reuseability between the two
% toolboxes.
%
% [x,fval,exitflag,output] = pso(...)
% The structure output contains more detailed information about the PSO
% run. It should match output fields of GA, where applicable.
%
% [x,fval,exitflag,output,population] = pso(...)
% A matrix population of size PopulationSize x nvars, with the locations of
% particles across the rows. This may be used to save the final positions
% of all the particles in a swarm.
%
% [x,fval,exitflag,output,population,scores] = pso(...)
% Final scores of the particles in population.
%
% Bibliography
% J Kennedy, RC Eberhart, YH Shi. Swarm Intelligence. Academic Press, 2001.
%
% SM Mikki, AA Kishk. Particle Swarm Optimization: A Physics-Based
% Approach. Morgan & Claypool, 2008.
%
% RE Perez and K Behdinan. "Particle swarm approach for structural
% design optimization." Computers and Structures, Vol. 85:1579-88, 2007.
%
% See also:
% PSODEMO, PSOOPTIMSET, PSOBINARY.
if ~nargin % Rosenbrock's banana function by default, as a demonstration
fitnessfcn = @(x)100*(x(2)-x(1)^2)^2+(1-x(1))^2 ;
nvars = 2 ;
LB = [-1.5,-2] ;
UB = [2,2] ;
options.PopInitRange = [[-2;4],[-1;2]] ;
options.PlotFcns = {@psoplotbestf,@psoplotswarmsurf} ;
options.Generations = 200 ;
options.DemoMode = 'on' ;
options.KnownMin = [1 1] ;
options.OutputFcns = {} ;
options.ConstrBoundary = 'penalize' ;
options.UseParallel = 'never' ;
elseif isstruct(fitnessfcn)
nvars = fitnessfcn.nvars ;
Aineq = fitnessfcn.Aineq ;
bineq = fitnessfcn.bineq ;
Aeq = fitnessfcn.Aeq ;
beq = fitnessfcn.beq ;
LB = fitnessfcn.LB ;
UB = fitnessfcn.UB ;
nonlcon = fitnessfcn.nonlcon ;
if ischar(nonlcon) && ~isempty(nonlcon)
nonlcon = str2func(nonlcon) ;
end
options = fitnessfcn.options ;
fitnessfcn = fitnessfcn.fitnessfcn ;
elseif nargin < 2
msg = 'PSO requires at least two input arguments' ;
error('%s, or a problem structure. Type >> help pso for details',...
msg)
end % if ~nargin
if ~exist('options','var') % Set default options
options = struct ;
end % if ~exist
options = psooptimset(options) ;
options.Verbosity = 1 ; % For options.Display == 'final' (default)
if strcmpi(options.Display,'off')
options.Verbosity = 0 ;
elseif strncmpi(options.Display,'iter',4)
options.Verbosity = 2 ;
elseif strncmpi(options.Display,'diag',4)
options.Verbosity = 3 ;
end
if ~exist('Aineq','var'), Aineq = [] ; end
if ~exist('bineq','var'), bineq = [] ; end
if ~exist('Aeq','var'), Aeq = [] ; end
if ~exist('beq','var'), beq = [] ; end
if ~exist('LB','var'), LB = [] ; end
if ~exist('UB','var'), UB = [] ; end
if ~exist('nonlcon','var'), nonlcon = [] ; end
% Check for swarm stability
% -------------------------------------------------------------------------
if options.SocialAttraction + options.CognitiveAttraction >= 4 && ...
options.Verbosity > 2
msg = 'Warning: Swarm is unstable and may not converge ' ;
msg = [msg 'since the sum of the cognitive and social attraction'] ;
msg = [msg ' parameters is greater than or equal to 4.'] ;
msg = [msg ' Suggest adjusting options.CognitiveAttraction and/or'] ;
sprintf('%s options.SocialAttraction.',msg)
end
% -------------------------------------------------------------------------
% Check for constraints and bit string population type
% ----------------------------------------------