godlike.m

一系列好用的用户友好的启发式优化算法

function [sol, fval, evals] = godlike(varargin)
%GODLIKE                 Optimization using the GODLIKE algorithm 
%
%   Usage:
%   sol = GODLIKE(PROBLEM)
%   sol = GODLIKE(func, popsize, lb, ub)    
%   sol = GODLIKE(func, popsize, lb, ub, 'option1', value1, ...)   
%
%   [sol, fval] = GODLIKE(...)
%   [sol, fval, evals] = GODLIKE(...)
%
%   GODLIKE(func, popsize, lb, ub) tries to find the global optimum of 
%   the fitness-function [func], using the GODLIKE-algorithm. The 
%   population size set by [popsize], and the boundaries for each 
%   dimension are set by the vectors [lb] and [ub], respectively. 
%
%   [sol, fval, evals] = GODLIKE(...) returns the trial vector found to 
%   yield the global minimum in [sol], and the corresponding function value 
%   by [fval]. The total amount of function evaluations that the algorithm
%   performed is returned in [evals].
%
%   The function [func] must accept a vector argument of length [N], equal to 
%   the number of elements in the vectors [lb] or [ub]. Also, the function 
%   must be vectorized so that inserting a matrix of [popsize]x[N] will return 
%   a vector of size [popsize]x[1] containing the corresponding function values
%   for the [N] trial vectors.
%
%   GODLIKE stands for "Global Optimum Determination by Linking and
%   Interchanging Kindred Evaluators". It basically links 4 different
%   optimization algorithms and interchanges their findings iteratively, in 
%   an attempt to negate the various problems all the other algorithms 
%   encounter.
%
%   The default control parameters GODLIKE uses, are:
%
%       AlgIters   = 250        The number of iterations spent in each
%                               heuristic optimizer.
%
%       Swapiters  = 10         The number of iterations before the
%                               populations are swapped between the 
%                               algorithms.
%
%       PutbackIters = 25       The number if iterations before the best
%                               individual ever found is inserted back into
%                               one of the populations. 
%
%   These values, as well as additional options, may be given manually in 
%   any extra input variables. Additionally, GODLIKE may be called as 
%   GODLIKE(PROBLEM), where [PROBLEM] is a complete problem-structure 
%   constructed by HEURSET.
%
%   Some optimizations require repeated evaluation of a function that takes 
%   in the order of hours per evaluation. In those cases, it might be 
%   convenient to interrupt the optimization after a few days and use the 
%   results thus far obtained. Unfortunately, usually after you interrupt a 
%   function, its results are lost. For that reason, GODLIKE creates two 
%   global variables, [GODLIKE_bestind] and [GODLIKE_bestfval], which store 
%   the best individual and function value thus far found. When the
%   optimization is interrupted, the intermediate results from the
%   optmization are still available. Once an optimization has succesfully
%   converged, these variables are deleted from the workspace again. 
%
%   See also heurset, diffevolve, genetic, swarm, simanneal.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: [e-mail]
% Last edited: 02/Feb/2009


%%  initialize & parse input
    
    % initial values
    problem = heurset;
    [func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...
     convvalue, swapiters, algiters, putbackiters] = parseprob(problem);
 
    % define temporary solution globals
    global GODLIKE_bestind GODLIKE_bestfval

    % common inputs
    if (nargin >= 4)
        func     =  varargin{1};
        popsize  =  varargin{2}; 
        lb       =  varargin{3};
        ub       =  varargin{4};        
    end
    
    % with additional options
    if (nargin >= 5)
        if ~isstruct(varargin{5})
            options = heurset(varargin{5:end});
        else
            options = varargin{5};                
        end
        [dum1, dum2, dum3, dum4, grace, display, maxfevals, convmethod,...
        convvalue, swapiters, algiters, putbackiters] = parseprob(options);    
    end
    
    % if given a full problem structure
    if (nargin == 1)
        problem  = varargin{1};        
        % errortrap
        if ~isstruct(problem)
            error('PROBLEM should be a structure. Type `help heurset` for details.')
        end   
        [func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...
        convvalue, swapiters, algiters, putbackiters] = parseprob(problem);
    end
    
    % initialize convergence method, using default values if omitted     
    if strcmpi(convmethod, 'exhaustive')
        convergence  = 1;
        maxiterinpop = convvalue;
        maxiters     = inf; 
    elseif strcmpi(convmethod, 'MaxIters')
        convergence  = 2;
        maxiterinpop = inf;
        maxiters     = convvalue; 
    elseif strcmpi(convmethod, 'achieveFval')
        convergence  = 3;         
        %errortrap
        if isempty(convvalue)
            error('Please define function value to achieve.')
        end        
        maxiterinpop = inf;
        maxiters     = inf;        
    else
        convergence  = 1;
        maxiterinpop = convvalue;
    end  
    
    % problem dimensions
    dims = size(lb, 2);
    
    % errortraps
    if ( (size(lb, 2) ~= 1 || size(ub, 2) ~= 1) &&...
         (size(lb, 2) ~= dims || size(ub, 2) ~= dims) )
        error('Upper- and lower boundaries must be the same size.')
    end
    if (popsize < 3)
        error('GODLIKE needs a population size of at least 3.')
    end   

%%  initial values

    % iteration-zero values
    firsttime     = true;
    iters         = 0;
    evals         = 0;
    globalbestfit = inf;
    improvement   = maxiterinpop;    
    converging1   = false;
    converging2   = false;

    % set popsize to be divisible by 4
    while (rem(popsize, 4) ~= 0)
        popsize = popsize + 1;
    end

    % initialize populations
    pop = repmat(ub, popsize, 1) - repmat(ub-lb, popsize, 1).* rand(popsize, dims);
    
    % initial fitnesses
    fits = inf(popsize, 1);    
    
    % update options structure
    problem.costfun = func;
    problem.lb = lb;
    problem.ub = ub;
    problem.popsize = popsize/4;
    problem.grace = 0;  
    problem.display = false;
    problem.conv.method = 'MaxIters';
    problem.conv.value  = algiters; 
    
    % initialize different evaluators     
    probGA = problem;     
    probDE = problem; 
    probPS = problem; 
    probSA = problem;    
    
    % Display strings
    if display    
        fprintf(1, ['\nNote that when GODLIKE is interrupted, the best solution and function\n',...
                    'value are saved in global variables GODLIKE_bestind and GODLIKE_bestfval.\n\n']);
        fprintf(1, 'Optimizing with ');
        switch convergence
            case 1
            	fprintf(1, 'exhaustive convergence criterion...\n');
            case 2
                fprintf(1, 'maximum iterations convergence criterion...\n');
            case 3
                fprintf(1, 'goal-attain convergence criterion...\n');
        end
        fprintf(1, 'Current best solution has cost of   ------------- ');
        overfevals1  = '\n\nCould not find better solution within given maximum';
        overfevals2  = '\nof function evaluations. Increase MaxFevals option.\n\n';
        fitbackspace = '\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b';
        converge2bck = [fitbackspace, '\b\b\b\b\b\b\b'];
        converge1bck = [converge2bck, '\b\b\b\b'];        
        convergestr  = ', possibly converging: ';
        itersstr     = ', iterations left: ';
    end 

%%  optimization loop

    % loop while one of the convergence criteria is not met
    while (improvement > 0 && iters <= maxiters)

        % if [iters] exceeds [swapiters], swap the populations randomly
        if (mod(iters, swapiters) == 0)
            % randomize the row-ordering            
            [ignore, inds] = sort(rand(1, popsize));
            pop = pop(inds, :);   
            
            % assign new populations to the different evaluators
            probGA.GA.pop = pop(1:4:end, :);
            probDE.DE.pop = pop(2:4:end, :);
            probPS.PS.pop = pop(3:4:end, :);
            probSA.PS.pop = pop(4:4:end, :);
            % probSA.SA.pop = pop(4:4:end, :);  % Simulated Annealing is horribly slow!!
        end
        
        % use all heuristic evaluators
        [pop(1:4:end, :), fits(1:4:end), evals1] = genetic([], probGA);    
        [pop(2:4:end, :), fits(2:4:end), evals2] = diffevolve([], probDE); 
        [pop(3:4:end, :), fits(3:4:end), evals3] = swarm([], probPS);      
        [pop(4:4:end, :), fits(4:4:end), evals4] = swarm([], probSA);  
        
        % update function evaluations
        evals = evals + evals1 + evals2 + evals3 + evals4;
        
        % find best and worst fitnesses
        [bestfit , bind] = min(fits);        
        [worstfit, wind] = max(fits);
        
        % improving solutions
        if (bestfit < globalbestfit)
            
            % new best individual
            improvement   = maxiterinpop;
            globalbestfit = bestfit;            
            globalbestind = pop(bind, :);
            
            % store also in globals
            GODLIKE_bestind  = globalbestind;
            GODLIKE_bestfval = globalbestfit;
            
            % display improving solution
            if display                
                if converging1
                   fprintf(1, converge1bck);
                   pause(0.05)
                end
                if converging2
                   fprintf(1, converge2bck);
                   pause(0.05)
                end                   
                if (globalbestfit < 0)
                    fprintf('\b')
                end
                fprintf(1, fitbackspace);
                fprintf(1, '%1.8e', globalbestfit);
                converging1 = false;
                firsttime   = true;
                pause(0.05)
            end           
        end

        % Check for convergence
        if (evals > maxfevals)
            
            % maximum allowed function evaluations has been superceded
            fprintf(overfevals1);
            fprintf(overfevals2);
            break
            
        end        
        switch convergence
            
            % exhaustive 
            case 1
                if (bestfit >= globalbestfit) && (globalbestfit ~= inf)
                   if (improvement <= 100) && display
                        if ~converging1
                           fprintf(1, convergestr);
                           fprintf(1, '%3.0f', improvement-1);                   
                           converging1 = true;
                           pause(0.05)
                        else
                            fprintf(1, '\b\b\b%3.0f', improvement-1);
                            pause(0.05)
                        end    
                    end           
                    improvement = improvement - 1;
                end
                
            % max iterations 
            case 2                   
                if display && (maxiters-iters < 100)            
                    if firsttime                
                        fprintf(1, itersstr);
                        fprintf(1, '%3.0f', maxiters-iters);                
                        pause(0.05)
                    else
                        fprintf(1, '\b\b\b%3.0f', maxiters-iters);
                        pause(0.05)
                    end
                    firsttime   = false;
                    converging2 = true;
                end  
                iters = iters + 1;
                
            % goal-attain
            case 3
                if (globalbestfit <= convvalue)
                    break;
                end
        end
        
        % replace worst individual with best individual
        if (mod(iters, putbackiters) == 0)
            pop(wind, :) = globalbestind;  
        end  
    end
    
%%  (pre-) end values
      
    % if a solution has been found
    if isfinite(globalbestfit)
        sol  = globalbestind;
        fval = globalbestfit;
        if display
            fprintf(1, '\nGODLIKE-algorithm has Converged.\n');
        end 
        
    % no feasible value might be found
    else
        fprintf(1,'\n');
        warning('godlike:no_solution',...
              'GODLIKE was unable to find any solution that gave finite values.')
        fval = globalbestfit;
        sol  = NaN;        
    end

%%  Grace function evaluations

    if (grace > 0)
        
        % display progress
        if display
            fprintf(1, 'Performing direct-search...');
            pause(0.05)
        end
        
        % perform direct search
        options = optimset('TolX', eps, 'MaxFunEvals', grace, 'TolFun', ...
                           eps, 'MaxIter', 1e4, 'Display', 'off');
        [soltry, fvaltry] = fminsearch(func, sol, options);
        
        % enforce boundaries
        if ~any(soltry >= ub | soltry <= lb)
           sol  = soltry;
           fval = fvaltry;
        end    
        evals = evals + grace;        
    end
    
%%  finalize
    
    % display progress
    if display
        fprintf(1, 'All done.\n\n');
        pause(0.05)
    end  
    
    % clear temp globals
    clear global GODLIKE_bestfval GODLIKE_bestind
    
end

% parser function to easily parse the input arguments
function [func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...
        convvalue, swapiters, algiters, putbackiters] = parseprob(problem)

        func         = problem.costfun;
        popsize      = problem.popsize; 
        lb           = problem.lb;
        ub           = problem.ub; 
        grace        = problem.grace;
        display      = problem.display;        
        maxfevals    = problem.maxfevals;
        convmethod   = problem.conv.method;
        convvalue    = problem.conv.value;
        swapiters    = problem.GODLIKE.swapiters;  
        algiters     = problem.GODLIKE.algiters;   
        putbackiters = problem.GODLIKE.putbackiters; 
        
end