function [x, fval, exitflag, output] = pso3(objfunc, nvars)
% e=xlsread('WellData1.1.xls','n score');
x=[1 2 3;4 5 6;1 2 6];
obj =x.^3-2*x.^2+5;
% 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
if nargin==1 && ischar(objfunc) && strcmp(objfunc, 'options')
if nargout<=1
x = getDefaultOptions();
else
error('mrr:myoptim:pso:pso:nargouterr', ...
'Cannot expext more than one output when only OPTIONS are required.');
end
else
if nargin<3, options=getDefaultOptions(); end
if nargout == 4
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
ptionserr
utput_level', ...
'Invalid value of the OUTPUT_LEVEL options specified.');
end
else
if options.plot == 1
output_level = 1;
else
output_level = 0;
end
end
if ~all(isnan(options.vmax))
if any(isnan(options.vmax))
error('mrr:myoptim:pso:pso
ptionserr:vmax', ...
'VMAX option cannot have some Inf and some numerical (or Inf) values.');
end
if ~isnan(options.vmaxscale)
warning('mrr:myoptim:pso:pso
ptionserr: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
if size(options.vmax, 1) ~= length(options.vmax)
error('mrr:myopim:pso:pso
ptionserr:vmax', ...
'VMAX option should be specified as column-vector, or as a scalar value.');
end
vmax = options.vmax;
else
error('mrr:myoptim:pso:pso
ptionserr:vmax', ...
'Inadequate dimension of VMAX option. Should be a scalar, or a column vector with NVARS elements.');
end
else
if isnan(options.vmaxscale)
error('mrr:myoptim:pso:pso
ptionserr:vmaxscale', ...
'Either VMAX or VMAXSCALE options should be different than NaN.');
end
if length(options.vmaxscale) == 1
if length(options.initspan) == 1
vmax = options.vmaxscale*options.initspan*ones(nvars, 1);
else
vmax = options.vmaxscale*options.initspan;
end
else
error('mrr:myoptim:pso:pso
ptionserr:vmax', ...
'Inadequate dimension of VMAXSCALE option. Must be a scalar.');
end
end
vmax = repmat(vmax', Ns, 1);
if ~isnan(options.initpopulation)
[pno, pdim] = size(options.initpopulation);
if (pno ~= Ns) || (pdim ~= nvars)
error('mrr:myoptim:pso:pso
ptionserr:initpopulation', ...
['The format of initial population is inconsistent with desired population', ...
'size or dimension of search space - INITPOPULATION options is invalid']);
end
X = options.initpopulation;
elseif (length(options.initoffset) == 1) && (length(options.initspan) == 1)
X = (rand(Ns, 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
ptionserr:initoffset_initspan', ...
'Both INITOFFSET and INITSPAN options must be either scalars or column-vectors.');
elseif (length(options.initoffset) ~= nvars) || (length(options.initspan) ~= nvars)
error('mrr:myoptim:pso:pso
ptionserr:init', ...
'Both INITOFFSET and INITSPAN options must be scalars or column-vectors of length NVARS.');
else
initoffset = repmat(options.initoffset', Ns, 1);
initspan = repmat(options.initspan', Ns, 1);
X = (rand(Ns, nvars)-0.5)*2.*initspan + initoffset;
if (options.trustoffset)
X(1, :) = options.initoffset';
end
end
if any(isnan(options.vspaninit))
error('mrr:myoptim:pso:pso
ptionserr:vspaninit', ...
'VSPANINIT option must not contain NaN entries.');
elseif isscalar(options.vspaninit)
V = (rand(Ns, nvars)-0.5)*2*options.vspaninit;
else
if (length(options.vspaninit) ~= size(options.vspaninit, 1)) || ...
(length(options.vspaninit) ~= nvars)
error('mrr:myoptim:pso:pso
ptionserr:vspaninit', ...
'VSPANINIT option must be either scalar or column-vector of length NVARS');
end
V = (rand(Ns, nvars)-0.5)*2.*repmat(options.vspaninit', Ns, 1);
end
Y = calcobjfunc(objfunc, X);
Ybest = Y;
Xbest = X;
[GYbest, gbest] = min(Ybest);
gbest = gbest(1);
tolbreak = ~isnan(options.globalmin);
foundglobal = 0;
if tolbreak && ~isscalar(options.globalmin)
error('mrr:myoptim:pso:pso
ptionserr:globalmin', ...
'globalmin option, if specified, option must be a scalar value equal to the global minimum of the objective function');
end
if output_level >= 0
output.itersno = K;
if output_level >= 1
output.gbest_array = NaN*ones(K+1, 1);
output.gmean_array = NaN*ones(K+1, 1);
output.gworst_array = NaN*ones(K+1, 1);
output.gbest_array(1) = GYbest;
output.gmean_array(1) = mean(Ybest);
output.gworst_array(1) = max(Ybest);
if output_level >= 2
output.gbestndx_array = NaN*ones(K+1, 1);
output.Xbest = NaN*ones(K+1, nvars);
output.gbestndx_array(1) = gbest;
output.Xbest(1, :) = X(gbest, :);
if output_level == 3
output.X = NaN*zeros(Ns, nvars, K+1);
output.X(:,:,1) = X;
end
end
end
end
if options.verbose_period ~= 0
disp 'PSO algorithm: Initiating the optimization process.'
end
exitflag = 0;
for iter = 1:K
if options.verbose_period ~= 0
if rem(iter, options.verbose_period) == 0
disp(['iteration ', int2str(iter), '. best criteria = ', num2str(GYbest)]);
end
end
w = linrate(options.wf, options.wi, K, 0, iter);
cp = linrate(options.cbf, options.cbi, K, 0, iter);
cg = linrate(options.cgf, options.cgi, K, 0, iter);
GXbest = repmat(Xbest(gbest, :), Ns, 1);
V = w*V + cp*rand(size(V)).*(Xbest-X) + cg*rand(size(V)).*(GXbest-X);
V = min(vmax, abs(V)).*sign(V);
X = X + V;
Y = calcobjfunc(objfunc, X);
mask = Y<Ybest;
mask = repmat(mask, 1, nvars);
Xbest = mask.*X +(~mask).*Xbest;
Ybest = min(Y,Ybest);
[GYbest, gbest] = min(Ybest);
gbest = gbest(1);
if output_level >= 0
if output_level >= 1
output.gbest_array(iter+1) = GYbest;
output.gmean_array(iter+1) = mean(Ybest);
output.gworst_array(iter+1) = max(Ybest);
if output_level >= 2
output.gbestndx_array(iter+1) = gbest;
output.Xbest(iter+1, :) = X(gbest, :);
if output_level == 3
output.X(:,:,iter+1) = X;
end
end
end
end
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
x = Xbest(gbest, :); x = x(:);
fval = GYbest;
if foundglobal, exitflag = 1; end;
if options.plot
r = 0:K;
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 [func]= calcobjfunc(func, X)
np = size(X,1);
Y = zeros(np,1);
for i = 1:np
Y(i) = func(X(i,
;
end
function opts = getDefaultOptions
%This function, in fact, defines default values of the options within the options structure.
opts.npart = 40; % The number of particles.
opts.niter = 20; % The number of iterations.
opts.cbi = 2.5; % Initial value of the individual-best acceleration factor.
opts.cbf = 0.5; % Final value of the individual-best acceleration factor.
opts.cgi = 0.5; % Initial value of the global-best acceleration factor.
opts.cgf = 2.5; % Final value of the global-best acceleration 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 primary speed limit.
opts.vmaxscale = NaN; % Relative speed limit. Used only if absolute limit is unspecified.
opts.vspaninit = 1; % The initial velocity span. Initial velocities 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 to be a good solution candidate, so it is included in
% the initial swarm.
opts.initpopulation = NaN; % The user-suplied initial population. If this is set to something
% meaningfull, then INITSPAN, INITOFFSET and TRUSTOFFSET are
% ignored.
opts.verbose_period = 10; % 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.
opts.plot = 1; % If set to 1, evolution of the gbest 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.
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 testing only
