Code covered by the BSD License  

Highlights from
OpenOpt

image thumbnail

OpenOpt

by

 

25 Nov 2006 (Updated )

nonSmoothSolve (similar to fsolve), non-smooth & noisy local + some global solvers; works in Octave

r=hPSO(prob)
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

Contact us