simanneal.m

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

function [sol, fval, evals] = simanneal(varargin)
%SIMANNEAL                       Simulated Annealing Optimization
%
%   Usage:
%   sol = SIMANNEAL(PROBLEM)
%   sol = SIMANNEAL(func, popsize, lb, ub)    
%   sol = SIMANNEAL(func, popsize, lb, ub, 'option1', value1, ...)   
%
%   [sol, fval] = SIMANNEAL(...)
%   [sol, fval, evals] = SIMANNEAL(...)
%   
%   SIMANNEAL(func, popsize, lb, ub) tries to find the global optimum of 
%   the fitness-function [func], using a population-based simulated annealing 
%   strategy. The population size set by [popsize], and the boundaries for
%   each dimension are set by the vectors [lb] and [ub], respectively. 
%
%   [sol, fval, evals] = SIMANNEAL(...) 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.
%
%   The default control parameters SIMANNEAL uses, are 
%
%     T0       = 1      Initial temperature
%
%     minT     = 1e-8   Final temperature
%
%     k        = 1      Bolzmann constant
%
%     maxiters = 150    Maximum number of iterations before the cooling 
%                       schedule is applied
%
%     cool     = (minT/T0)^(1/maxiters);   cooling schedule                                           
%
%   These values, as well as additional options, may be given manually in 
%   any extra input variables. Additionally, SIMANNEAL may be called as 
%   SIMANNEAL(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, SIMANNEAL creates two 
%   global variables, [SIMANNEAL_bestind] and [SIMANNEAL_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, genetic, swarm, diffevolve, godlike.

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

%%  initialize & parse input
    
    % default values    
    skippop = false; 
    options = heurset;
    [pop, func, popsize, lb, ub, grace, display, ...
     convmethod, convvalue, T0, minT, k, cool] = parseprob(options);
 
    % define temporary globals
    global SIMANNEAL_bestind SIMANNEAL_bestfval
   
    % common inputs
    if (nargin >= 4)
        func     =  varargin{1};
        popsize  =  varargin{2}; 
        lb       =  varargin{3};
        ub       =  varargin{4};        
    end
    
    % additional inputs
    if (nargin >= 5) 
            if ~isstruct(varargin{5})
                options = heurset(varargin{5:end});
            else
                options = varargin{5};                
            end
            [dum1, dum2, dum3, dum4, dum5, grace, display, ...
            convmethod, convvalue, T0, minT, k, cool] = parseprob(options);
    end
    
    % if called from GODLIKE
    if (nargin == 2)  
        problem = varargin{2};  
        % errortrap
        if ~isstruct(problem)
            error('PROBLEM should be a structure. Type `help heurset` for details.')
        end 
        [pop, func, popsize, lb, ub, grace, display, convmethod, ...
         convvalue, T0, minT, k, cool] = parseprob(problem);
        
        % make sure the options are correct        
        convmethod = 'maxiters'; 
        grace   = 0;
        display = false;
        skippop = true;
    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    
        [dum1, func, popsize, lb, ub, grace, display, convmethod, ...
         convvalue, T0, minT, k, cool] = parseprob(problem);
    end    
    
    % initialize convergence method, setting the default values if omitted     
    if strcmpi(convmethod, 'MaxIters')
        convergence  = 1;
        maxiters     = convvalue; 
    else
        convergence  = 1;
        maxiters     = 100;
    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
    
%%  initial values 

    % iteration-zero values
    iters       = 0;
    evals       = 0;
    done        = false;
    skip        = false;       
    firsttime   = true;  
    tempdisp    = false;
    T           = T0;
    verybest    = inf;
    verybestval = inf;
    prevbestval = inf;
    nrg         = inf(popsize, 1);
    prevnrg     = nrg;    
    
    % boundaries
    if (numel(lb) == 1)
        mins = repmat(lb, popsize, dims);
        maxs = repmat(ub, popsize, dims); 
    end
    if (numel(lb) == numel(ub) && size(lb, 2) == dims)
        mins = repmat(lb, popsize, 1);
        maxs = repmat(ub, popsize, 1);                        
    end
    range = (maxs - mins);

    % initialize population   
    if (~skippop)
        pop = repmat(ub, popsize, 1) - repmat(ub-lb, popsize, 1).* rand(popsize, dims);
    end   
    
    % Display strings
    if display
        fprintf(1, ['\nNote that when SIMANNEAL is interrupted, the best solution and function\n',...
                    'value are saved in global variables SIMANNEAL_bestind and SIMANNEAL_bestfval.\n\n']);
        fprintf(1, 'Optimizing with ');
        switch convergence            
            case 1
                fprintf(1, 'maximum iterations convergence criterion...\n');            
        end
        fprintf(1, 'Current best solution has cost of   --------------- ');
        tempstr = ', Temperature is now: ';
        tempbck = '\b\b\b\b\b\b\b\b\b\b\b\b\b';
        fitbck  = [tempbck, '\b\b\b\b'];              
        convbck = [fitbck, fitbck, tempbck, '\b\b\b'];
    end
    
%%  Simulated Annealing loop     

    % loop while not done
    while ~done;
        
        % increase iterations
        iters = iters + 1;
    
        % convergence and decrease T criteria
        if (iters >= maxiters);
            if (T < minT)
                done = true;
                skip = true;
            else
                
                % apply cooling schedule
                T = T * cool(T, T0, minT, maxiters);
                
                % display progress
                if display    
                    if ( verybestval < prevbestval)
                        prevbestval = verybestval;
                        if tempdisp
                            fprintf(1, convbck);
                        else
                            fprintf(1, fitbck); 
                        end
                        if (verybestval < 0)
                            fprintf('\b')
                        end
                        fprintf(1, '%1.8e', verybestval);
                        firsttime = true;
                        pause(0.05)
                    end  
                    
                    % convergence
                    switch convergence
                        case 1                   
                            if display 
                                if firsttime                
                                    fprintf(1, tempstr);
                                    fprintf(1, '%1.6e', T);                
                                    pause(0.05)
                                else
                                    fprintf(1, tempbck);
                                    fprintf(1, '%1.6e', T);
                                    pause(0.05)
                                end
                                firsttime  = false; 
                                tempdisp   = true;
                            end 
                    end
                end  
            end
            
            % reset iterations
            iters = 0;
            
        end
        
        % generate new population
        rands         = randn(popsize, dims);
        [dummy, inds] = sort(rands, 2);
        locations     = (inds == 1);
        newpop        = pop + locations.*rands.*range;        
        
        % enforce boundaries
        outsiders1 = (newpop < mins);
        outsiders2 = (newpop > maxs);
        if any(outsiders1(:)) || any(outsiders2(:))
            reinit     = rand(popsize, dims).*range + mins;              
            newpop(outsiders1) = reinit(outsiders1);
            newpop(outsiders2) = reinit(outsiders2);
        end
        
        % find their energies
        if (~skip)
            prevpop = newpop;
            try
                nrg = feval(func, newpop);
            catch
                error('simanneal:fevalerror', ...
                     ['Simanneal cannot continue because the supplied cost function ',...
                      'gave the following error:\n %s'], lasterr);
            end   
            if (numel(nrg) ~= popsize)
                 error('simanneal:function_not_vectorized_or_incorrect_dimensions', ...
                 ['The user-supplied cost function does not appear to get enough arguments,\n',...
                  'or is not vectorized properly. Make sure the dimensions of the limits [Lb]\n',...
                  'and [Ub] are the same as the required input vector for the function, or that\n',...
                  'the function is vectorized properly.'])            
            end
            nrgval  = nrg;
        end
        
        % increase evals counter
        evals = evals + numel(nrg);

        % reject or accept them, according to the probabilistic rule
        nrgdiff  = (prevnrg - nrg);
        if iters > 2
            nrgdiff = nrgdiff / max(abs(nrgdiff(:)));
        end
        ind          = nrgdiff > 0;
        probind      = ~ind & rand(popsize, 1) < exp( nrgdiff /(k*T) );        
        pop(ind)     = newpop(ind);
        pop(probind) = newpop(probind);
        prevnrg(ind) = nrg(ind);
        prevnrg(probind) = nrg(probind);
        
        % best individual
        [bestval, inds] = min(nrg(:));
        bestind = newpop(inds, :);  
        if (bestval < verybestval)
            verybest    = bestind;
            verybestval = bestval; 
            
            % assign also the globals
            SIMANNEAL_bestind  = bestind;
            SIMANNEAL_bestfval = bestval;
            
        end
        
        % make sure the best individual is always in the population        
        inds      = round(rand)*(popsize-1) + 1; 
        pop(inds, :) = verybest;        
        nrg(inds) = verybestval;
        
    end

%%  (pre-) end values

    % if solution has been found
    if isfinite(verybestval)
        
        % when called normally
        if (~skippop)
            fval = verybestval;
            sol  = verybest;
            
        % when called from GODLIKE
        else            
            fval = nrgval;
            sol  = prevpop;
        end
        
        % display final message
        if display
            fprintf(1, '\nSimulated Annealing algorithm has converged.\n');
            pause(0.05)
        end
        
    % all trials might be infeasible 
    else
        fprintf(1,'\n');
        warning('simanneal:no_solution',...
              'SIMANNEAL was unable to find any solution that gave finite values.')
        fval = verybestval;
        sol  = NaN;
    end
    
%%  allow direct-search, with [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 the temp globals
    clear global SIMANNEAL_bestind SIMANNEAL_bestfval
    
end

% parser function to easily parse the input arguments
function [pop, func, popsize, lb, ub, grace, display, ...
        convmethod, convvalue, T0, minT, k, cool] = parseprob(problem)
   
        func       = problem.costfun;
        popsize    = problem.popsize; 
        lb         = problem.lb;
        ub         = problem.ub; 
        grace      = problem.grace;
        display    = problem.display;        
        convmethod = problem.conv.method;
        convvalue  = problem.conv.value; 
        pop        = problem.SA.pop;
        T0         = problem.SA.T0;
        minT       = problem.SA.minT;    
        k          = problem.SA.k;   
        cool       = problem.SA.cool; 
        
end