function r=hPSO(prob)
%Syntax: [x,fval,gfx,output]=hPSO(fitnessfun,nvars,options,varargin)
%uses LbUbWrapper()
%___________________________________________________________________
%
% A hybrid Particle Swarm Optimization algorithm for finding the minimum of
% the function 'fitnessfun' in the real space.
%
% x is the scalar/vector of the functon minimum
% fval is the function minimum
% gfx contains the best particle for each flight (columns 2:end) and the
% corresponding solutions (1st column)
% output structure contains the following information
% reason : is the reason for stopping
% flights: the nuber of flights before stopping
% time : the total time before stopping
% fitnessfun is the function to be minimized
% nvars is the number of variables of the problem
% options are specified in the file "PSOoptions.m"
%
%
% Reference:
% Kennedy J., Eberhart R.C. (1995): Particle swarm optimization. In: Proc.
% IEEE Conf. on Neural Networks, IV, Piscataway, NJ, pp. 19421948
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 17 Nov, 2004.
r = struct('alg', 'Kennedy J., Eberhart R.C. (1995): Particle swarm optimization. In: Proc.',...
'authors', 'Alexandros Leontitsis, me00743@cc.uoi.gr, University of Ioannina, Greece',...
'lastChanger', 'Dmitrey Kroshko, openopt@ukr.net');
options=hPSOoptions;
c1 = options.c1;
c2 = options.c2;
w = options.w;
maxv = options.maxv;
popul = options.bees;
MaxNonEnhance = 4;
ExtractRoutineParamsFromProb;
if size(maxv,1)==1; maxv=maxv*ones(1, prob.n);
elseif size(maxv,1)~=prob.n
error('The rows of options.maxv are not equal to nvars.');
end
if ~exist('maxInnerIter', 'var')
prob.warn('hPSO requires prob.hPSO.maxInnerIter field; setting to default 400')
maxInnerIter = 400;
end
prob.assert(~any(isinf(prob.lb)) || ~any(isinf(prob.ub)), 'hPSO solver can''t handle inf in lb & ub');
prob2 = prob;
prob2.iterPrint=0;
if ~isfield(prob, 'hPSO') || ~isfield(prob.hPSO, 'merit')
% prob.warn('missing prob.hPSO.merit; default (UkrOpt ralg) will be used')
merit = @ralg;
end
prob2.iterfcn = @innerIterFCN;
prob2.MaxIter = maxInnerIter;
% if ~(prob.classF>0 && prob.classC>0)
% prob.err('this solver can''t run this task')
% end
prob.assert(all(~isinf(prob.lb)) && all(~isinf(prob.ub)), 'hPSO requires finite lb-ub bounds')
% flights = options.flights;
% StallFliLimit = options.StallFliLimit;
% Goal = options.Goal;
% Initial population (random start)
ru=rand(popul,prob.n);
pop=ones(popul,1)*prob.lb'+ru.*(ones(popul,1)*(prob.ub-prob.lb)');
r.ff = inf;
if isreal(prob.x0); pop(1,:) = prob.x0'; end
% Hill climb of each solution (bee)
fxi = zeros(popul,1);
for i=1:popul
prob2.x0 = pop(i,:)';
r2 = feval(merit,prob2);
% disp('DONE')
pop(i,:) = r2.xf';
fxi(i) = r2.ff;
end
% Local minima
p=pop;
fxip=fxi;
[r.ff ind] = min(fxi);
r.xf = pop(ind,:)';
r = prob.iterfcn(prob, r);
if r.istop; return; end
% Initialize the velocities
v=zeros(popul,prob.n);
% Isolate the best solution
[Y,I]=min(fxi);
gfx(1,:)=[Y pop(I,:)];
P=ones(popul,1)*pop(I,:);
output.reason = 'Optimization terminated: maximum number of flights reached.';
counter = 0;
% For each flight
while 1
% Estimate the velocities
r1=rand(popul,prob.n);
r2=rand(popul,prob.n);
v=v*w+c1*r1.*(p-pop)+c2*r2.*(P-pop);
v=max(v,-ones(popul,1)*maxv);
v=min(v,ones(popul,1)*maxv);
% Add the velocities to the population
pop=pop+v;
% Drag the particles into the search space
pop=min(pop,ones(popul,1)*prob.ub');
pop=max(pop,ones(popul,1)*prob.lb');
% Hill climb search for the new population
pnew=p;
fxipnew=fxip;
for j=1:popul
prob2.x0 = pop(j,:)';
r2 = feval(merit,prob2);
pop(j,:) = r2.xf';
fxi(j) = r2.ff;
prob2.x0 = pop(j,:)';
r2 = feval(merit,prob2);
pnew(j,:) = r2.xf';
fxipnew(j) = r2.ff;
end
% if isfield(r2, 'df'); r2 = rmfield(r2, 'df'); end
% Min(fxi,fxip)
s=find(fxi<fxip, 1);
p(s,:)=pop(s,:);
fxip(s)=fxi(s);
% Min(fxipnew,fxip);
s=find(fxipnew<fxip);
p(s,:)=pnew(s,:);
fxip(s)=fxipnew(s);
% Isolate the best solution
Y=min(fxip);
I = find(Y==fxip,1);
gfx(i,:)=[Y p(I,:)];
P=ones(popul,1)*p(I,:);
r.x = p(I,:);
r.x = r.x(:);
r.f = Y;
if r.f<r.ff
r.xf = r.x;
r.ff = r.f;
counter = 0;
else
counter = counter+1;
if counter>MaxNonEnhance
r.istop = 5;
end
end
r = prob.iterfcn(prob, r);
if r.istop; return; end
end
end
function options=hPSOoptions(varargin)
%Syntax: options=hPSOoptions(varargin)
%_____________________________________
%
% Options definition for PSO.
%
% options is the options struct:
% 'c1': the strength parameter for the local attractors
% 'c2': the strength parapeter for the global attractor
% 'w': the velocity decline parameter
% 'maxv': the maximum possible velocity for each dimension
% 'space': the 2-column matrix with the boundaries of each
% particle
% 'bees': the number of the population particles
% 'flights': the number of flights
% 'HybridIter': the maximum number of iterations allowed for the
% @fminsearch
% 'Show': displays the progress (or part of it) on the screen
% 'StallFliLimit': the number of flights after the last improvement before
% the optimization stops
% 'Goal': a value to be reached
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110- Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 17 Nov, 2004.
options.c1 = 2;
options.c2 = 2;
options.w = 0;
options.maxv = inf;
options.space = [0 1];
options.bees = 20;
options.flights = 50;
options.HybridIter = 0;
options.Show = 'Final';
options.StallFliLimit = 50;
options.Goal = -inf;
if (nargin == 0) && (nargout == 0)
fprintf(' c1: [ real positive scalar | {2} ]\n');
fprintf(' c2: [ real positive scalar | {2} ]\n');
fprintf(' w: [ scalar in [0 1) | {0} ]\n');
fprintf(' maxv: [ real positive vector | {ones(nvars,1)*inf} ]\n');
fprintf(' space: [ real matrix | {[0 1]} ]\n');
fprintf(' bees: [ integer positive scalar | {20} ]\n');
fprintf(' flights: [ integer positive scalar | {50} ]\n');
fprintf(' HybridIter: [ integer non-negative scalar | {0} ]\n');
fprintf(' Show: [ integer non-negative scalar | {"Final"} ]\n');
fprintf(' StallFliLimit: [ real positive scalar | {50} ]\n');
fprintf(' Goal: [ real scalar | {-inf} ]\n');
end
if nargin > 1
if rem(nargin,2)==0
for i=1:nargin/2
switch varargin{i*2-1}
case {'c1','c2','bees','flights','StallFliLimit'}
% varargin{i*2} must be a scalar
if sum(size(varargin{i*2}))>2
error([varargin{i*2-1} ' must be a scalar.']);
end
% varargin{i*2} must be positive
if varargin{i*2}<=0
error([varargin{i*2-1} ' must be positive.']);
end
case {'maxv'}
% All the elements of varargin{i*2} must be positive
varargin{i*2}=varargin{i*2}(:);
if sum(any(varargin{i*2}))>0
error(['All the elements of ' varargin{i*2-1} ' must be positive.'])
end
case {'HybridIter','Show'}
if isa(varargin{i*2},'numeric')==1
% varargin{i*2} must be a scalar
if sum(size(varargin{i*2}))>2
error([varargin{i*2-1} ' must be a scalar.']);
end
% varargin{i*2} must be positive
if varargin{i*2}<=0
error([varargin{i*2-1} ' must be positive.']);
end
% varargin{i*2} must be an integer
if varargin{i*2}~=round(varargin{i*2})
error([varargin{i*2-1} ' must be an integer.']);
end
elseif strcmp(varargin{i*2-1},'Show')==1 && strcmp(varargin{i*2},'Final')==0
error('The only non-numeric value of "Show" is "Final".');
end
case 'space'
% varargin{i*2} must be a 2 dimensinal matrix
if ndims(varargin{i*2})>3
error(['"' varargin{i*2-1} '" must be a 2 dimensinal matrix.']);
end
% varargin{i*2-1} must contain exactly 2 columns
if size(varargin{i*2},2)~=2
error(['"' varargin{i*2-1} '" must contain exactly 2 columns.']);
end
case 'w'
% varargin{i*2} must be in [0 1)
if varargin{i*2}<0 || varargin{i*2}>=1
error(['"' varargin{i*2-1} '" must be in [0 1).']);
end
case 'Goal'
% varargin{i*2} must be a scalar
if sum(size(varargin{i*2}))>2
error([varargin{i*2-1} ' must be a scalar.']);
end
otherwise
error(['The field "' varargin{i*2-1} '" does not exist in the struct PSOoptions.']);
end
options = setfield(options,varargin{i*2-1},varargin{i*2}); %#ok<SFLD>
fprintf(' %s updated... OK\n',varargin{i*2-1});
end
else
error('The number of input arguments must be even.');
end
end
end
