function [X,FVAL,EXITFLAG,OUTPUT] = PSO(FUN,X0,LB,UB,OPTIONS,varargin)
%PSO finds a minimum of a function of several variables using the particle swarm
% optimization (PSO) algorithm originally introduced in 1995 by Kennedy and
% Eberhart. This algorithm was extended by Shi and Eberhart in 1998 through the
% introduction of inertia factors to dampen the velocities of the particles.
% In 2002, Clerc and Kennedy introduced a constriction factor in PSO, which was
% later on shown to be superior to the inertia factors. Therefore, the algorithm
% using a constriction factor was implemented here.
%
% PSO attempts to solve problems of the form:
% min F(X) subject to: LB <= X <= UB
% X
%
% X=PSO(FUN,X0) start at X0 and finds a minimum X to the function FUN.
% FUN accepts input X and returns a scalar function value F evaluated at X.
% X0 may be a scalar, vector, or matrix.
%
% X=PSO(FUN,X0,LB,UB) defines a set of lower and upper bounds on the
% design variables, X, so that a solution is found in the range
% LB <= X <= UB. Use empty matrices for LB and UB if no bounds exist.
% Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is
% unbounded above.
%
% X=PSO(FUN,X0,LB,UB,OPTIONS) minimizes with the default optimization
% parameters replaced by values in the structure OPTIONS, an argument
% created with the PSOSET function. See PSOSET for details.
% Used options are SWARM_SIZE, COGNITIVE_ACC, SOCIAL_ACC, MAX_ITER,
% MAX_TIME, MAX_FUN_EVALS, TOLX, TOLFUN, DISPLAY and OUTPUT_FCN.
% Use OPTIONS = [] as a place holder if no options are set.
%
% X=PSO(FUN,X0,LB,UB,OPTIONS,varargin) is used to supply a variable
% number of input arguments to the objective function FUN.
%
% [X,FVAL]=PSO(FUN,X0,...) returns the value of the objective
% function FUN at the solution X.
%
% [X,FVAL,EXITFLAG]=PSO(FUN,X0,...) returns an EXITFLAG that describes the
% exit condition of PSO. Possible values of EXITFLAG and the corresponding
% exit conditions are:
%
% 1 Change in the objective function value less than the specified tolerance.
% 2 Change in X less than the specified tolerance.
% 0 Maximum number of function evaluations or iterations reached.
% -1 Maximum time exceeded.
%
% [X,FVAL,EXITFLAG,OUTPUT]=PSO(FUN,X0,...) returns a structure OUTPUT with
% the number of iterations taken in OUTPUT.nITERATIONS, the number of function
% evaluations in OUTPUT.nFUN_EVALS, the coordinates of the different particles in
% the swarm in OUTPUT.SWARM, the corresponding fitness values in OUTPUT.FITNESS,
% the particle's best position and its corresponding fitness in OUTPUT.PBEST and
% OUTPUT.PBEST_FITNESS, the best position ever achieved by the swarm in
% OUTPUT.GBEST and its corresponding fitness in OUTPUT.GBEST_FITNESS, the amount
% of time needed in OUTPUT.TIME and the options used in OUTPUT.OPTIONS.
%
% See also PSOSET, PSOGET
% Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
% USA.
% handle variable input arguments
if nargin < 5,
OPTIONS = [];
if nargin < 4,
UB = 1e5;
if nargin < 3,
LB = -1e5;
end
end
end
% check input arguments
if ~ischar(FUN),
error('''FUN'' incorrectly specified in ''PSO''');
end
if ~isfloat(X0),
error('''X0'' incorrectly specified in ''PSO''');
end
if ~isfloat(LB),
error('''LB'' incorrectly specified in ''PSO''');
end
if ~isfloat(UB),
error('''UB'' incorrectly specified in ''PSO''');
end
if length(X0) ~= length(LB),
error('''LB'' and ''X0'' have incompatible dimensions in ''PSO''');
end
if length(X0) ~= length(UB),
error('''UB'' and ''X0'' have incompatible dimensions in ''PSO''');
end
% declaration of global variables
global NDIM nFUN_EVALS
% set EXITFLAG to default value
EXITFLAG = -2;
% determine number of variables to be optimized
NDIM = length(X0);
% seed the random number generator
rand('state',sum(100*clock));
% set default options
DEFAULT_OPTIONS = PSOSET('SWARM_SIZE',25,... % number of particles in swarm for each variable to be optimized
'COGNITIVE_ACC',2.8,... % cognitive acceleration coefficient (value as recommended in Schutte and Groenwold, 2006)
'SOCIAL_ACC',1.3,... % social acceleration coefficient (value as recommended in Schutte and Groenwold, 2006)
'MAX_ITER',2500,... % maximum number of iterations
'MAX_TIME',2500,... % maximum duration of optimization
'MAX_FUN_EVALS',2500,... % maximum number of function evaluations
'TOLX',1e-6,... % maximum difference between best and worst function evaluation in simplex
'TOLFUN',1e-3,... % maximum difference between the coordinates of the vertices
'DISPLAY','none',... % 'iter' or 'none' indicating whether user wants feedback
'OUTPUT_FCN',[]); % string with output function name
% update default options with supplied options
OPTIONS = PSOSET(DEFAULT_OPTIONS,OPTIONS);
% store OPTIONS in OUTPUT
OUTPUT.OPTIONS = OPTIONS;
% initialize swarm (each row of swarm corresponds to one particle)
SWARM = zeros(OPTIONS.SWARM_SIZE,NDIM,OPTIONS.MAX_ITER);
for i = 1:OPTIONS.SWARM_SIZE,
if i == 1,
SWARM(1,:,1) = X0(:)';
else
SWARM(i,:,1) = LB(:)' + rand(1,NDIM).*(UB(:)'-LB(:)');
end
end
% initialize VELOCITIES
VELOCITIES = zeros(OPTIONS.SWARM_SIZE,NDIM,OPTIONS.MAX_ITER);
% initialize FITNESS, PBEST_FITNESS, GBEST_FITNESS, INDEX_PBEST, index_gbest_particle and INDEX_GBEST_ITERATION
FITNESS = nan(OPTIONS.SWARM_SIZE,OPTIONS.MAX_ITER);
PBEST = nan(OPTIONS.SWARM_SIZE,NDIM,OPTIONS.MAX_ITER);
GBEST = nan(OPTIONS.MAX_ITER,NDIM);
PBEST_FITNESS = nan(OPTIONS.SWARM_SIZE,OPTIONS.MAX_ITER);
GBEST_FITNESS = nan(OPTIONS.SWARM_SIZE,1);
INDEX_PBEST = nan(OPTIONS.SWARM_SIZE,OPTIONS.MAX_ITER);
INDEX_GBEST_PARTICLE = nan(OPTIONS.MAX_ITER,1);
INDEX_GBEST_ITERATION = nan(OPTIONS.MAX_ITER,1);
% calculate constriction factor from acceleration coefficients
if OPTIONS.COGNITIVE_ACC+OPTIONS.SOCIAL_ACC <= 4,
% display message
disp('Sum of Cognitive Acceleration Coefficient and Social Acceleration Coefficient is less then or equal to 4.')
disp('Their values were adjusted automatically to satisfy this condition.');
disp(' ')
% the values are adjusted so that the sum is equal to 4.1, keeping the ratio COGNITIVE_ACC/SOCIAL_ACC constant
OPTIONS.COGNITIVE_ACC = OPTIONS.COGNITIVE_ACC*4.1/(OPTIONS.COGNITIVE_ACC+OPTIONS.SOCIAL_ACC);
OPTIONS.SOCIAL_ACC = OPTIONS.SOCIAL_ACC*4.1/(OPTIONS.COGNITIVE_ACC+OPTIONS.SOCIAL_ACC);
% calculate constriction factor
k = 1; % k can take values between 0 and 1, but is usually set to one (Montes de Oca et al., 2006)
OPTIONS.ConstrictionFactor = 2*k/(abs(2-(OPTIONS.COGNITIVE_ACC+OPTIONS.SOCIAL_ACC)-sqrt((OPTIONS.COGNITIVE_ACC+OPTIONS.SOCIAL_ACC)^2-4*(OPTIONS.COGNITIVE_ACC+OPTIONS.SOCIAL_ACC))));
else
% calculate constriction factor
k = 1; % k can take values between 0 and 1, but is usually set to one (Montes de Oca et al., 2006)
OPTIONS.ConstrictionFactor = 2*k/(abs(2-(OPTIONS.COGNITIVE_A
评论0