diffevolve.m

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

function [sol, fval, evals] = diffevolve(varargin)
%DIFFEVOLVE                          Differential Evolution Optimization
%
%   Usage:
%   sol = DIFFEVOLVE(PROBLEM)
%   sol = DIFFEVOLVE(func, popsize, lb, ub)    
%   sol = DIFFEVOLVE(func, popsize, lb, ub, 'option1', value1, ...)   
%
%   [sol, fval] = DIFFEVOLVE(...)
%   [sol, fval, evals] = DIFFEVOLVE(...)
%   
%   DIFFEVOLVE(func, popsize, lb, ub) tries to find the global optimum of 
%   the fitness-function [func], using a transversal differential evolution 
%   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] = DIFFEVOLVE(...) 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 DIFFEVOLVE uses, are 
%
%       -1 < F < 1          scaling parameter for the difference vector.
%                           By default, it is randomly selected from the
%                           range [-1, 1] for each individual in each
%                           iteration. 
%
%       crossconst  = 0.9   Probability for crossover.
%
%       convergence = 100   Maximum amount of iterations the best
%                           individual ever found should remain the best
%                           individual before the algorithm converges.
%
%   These values, as well as additional options, may be given manually in 
%   any extra input variables. Additionally, DIFFEVOLVE may be called as 
%   DIFFEVOLVE(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, DIFFEVOLVE creates two 
%   global variables, [DIFFEVOLVE_bestind] and [DIFFEVOLVE_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, simanneal, godlike.

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

%%  initialize & parse input 
       
    % initial (default) values
    skippop = false;
    options = heurset;
    [pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...
     convvalue, crossconst, Flb, Fub, n] = parseprob(options);
 
    % temporary globals
    global DIFFEVOLVE_bestind DIFFEVOLVE_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, dum5, grace, display, maxfevals, convmethod,...
        convvalue, crossconst, Flb, Fub, n] = 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, maxfevals, convmethod, ...
        convvalue, crossconst, Flb, Fub, n] = 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  
        [pop, func, popsize, lb, ub, grace, display, maxfevals, convmethod, ...
        convvalue, crossconst, Flb, Fub, n] = parseprob(problem);
    end
    
    % initialize convergence method
    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('DIFFEVOLVE needs a population size of at least 3.')
    end  
    
%%  initial values

    % iteration-zero values
    oldbestfit  = inf;
    improvement = maxiterinpop;
    iters       = 0;
    evals       = 0; 
    sizepop     = [popsize, dims];    
    fitnew      = inf(popsize, 1);
    firsttime   = true;
    converging1 = false;
    converging2 = false;
    
    % 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
    newpop  = pop;
    prevpop = pop;     
    
    % Display strings
    if display   
        fprintf(1, ['\nNote that when DIFFEVOLVE is interrupted, the best solution and function\n',...
                    'value are saved in global variables DIFFEVOLVE_bestind and DIFFEVOLVE_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';
        converge1bck = [fitbackspace, '\b\b\b\b\b\b\b\b\b\b\b'];
        converge2bck = [fitbackspace, '\b\b\b\b\b\b\b'];
        convergestr  = ', possibly converging: ';
        itersstr     = ', iterations left: ';
    end    
    
%%  Differential Evolution loop

    % loop while one of the convergence criteria is not met
    while (improvement > 0 && iters <= maxiters)
                
        % evaluate fitnesses and adjust population
        fitold  = fitnew;
        try
            fitnew  = feval(func, newpop);
        catch
            error('diffevolve:fevalerror', ...
                 ['Diffevolve cannot continue because the supplied cost function ',...