sce.m

SCE（shuffled complex evolution ）是一种相对较新的连续性问题的元启发搜索算法。非常适合于求解具有多个局部最小的全局优化问题。SCE算法的主要特征是通过竞争进化和定期洗牌来

function [X,FVAL,EXITFLAG,OUTPUT] = SCE(FUN,X0,LB,UB,OPTIONS,varargin)
%SCE finds a minimum of a function of several variables using the shuffled
% complex evolution (SCE) algorithm originally introduced in 1992 by Duan et al. 
%
%   SCE attempts to solve problems of the form:
%       min F(X) subject to: LB <= X <= UB
%        X
%
%   X=SCE(FUN,X0) start at X0 and finds a minimum X to the function FUN. 
%   FUN accepts input X and returns a scalar function value F evaluated at X.
%   X0 may be a scalar, vector, or matrix.
%   
%   X=SCE(FUN,X0,LB,UB) defines a set of lower and upper bounds on the 
%   design variables, X, so that a solution is found in the range 
%   LB <= X <= UB. Use empty matrices for LB and UB if no bounds exist. 
%   Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is 
%   unbounded above.
%   
%   X=SCE(FUN,X0,LB,UB,OPTIONS) minimizes with the default optimization
%   parameters replaced by values in the structure OPTIONS, an argument 
%   created with the SCESET function. See SCESET for details. 
%   Used options are nCOMPLEXES, nITER_INNER_LOOP, MAX_ITER,
%   MAX_TIME, MAX_FUN_EVALS, TOLX, TOLFUN, DISPLAY and OUTPUT_FCN.
%   Use OPTIONS = [] as a place holder if no options are set.
%   
%   X=SCE(FUN,X0,LB,UB,OPTIONS,varargin) is used to supply a variable 
%   number of input arguments to the objective function FUN.
%   
%   [X,FVAL]=SCE(FUN,X0,...) returns the value of the objective 
%   function FUN at the solution X.
%   
%   [X,FVAL,EXITFLAG]=SCE(FUN,X0,...) returns an EXITFLAG that describes the 
%   exit condition of SCE. Possible values of EXITFLAG and the corresponding 
%   exit conditions are:
%   
%     1  Change in the objective function value less than the specified tolerance.
%     2  Change in X less than the specified tolerance.
%     0  Maximum number of function evaluations or iterations reached.
%    -1  Maximum time exceeded.
%   
%   [X,FVAL,EXITFLAG,OUTPUT]=SCE(FUN,X0,...) returns a structure OUTPUT with 
%   the number of iterations taken in OUTPUT.nITERATIONS, the number of function
%   evaluations in OUTPUT.nFUN_EVALS, the different points in the population 
%   at every iteration and their fitness in OUTPUT.POPULATION and 
%   OUTPUT.POPULATION_FITNESS respectively, the amount of time needed in 
%   OUTPUT.TIME and the options used in OUTPUT.OPTIONS.
% 
%   See also SCESET, SCEGET



% Copyright (C) 2006 Brecht Donckels, BIOMATH, [e-mail]
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. 
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,
% USA.



% handle variable input arguments

if nargin < 5,
    OPTIONS = [];
    if nargin < 4,
        UB = 1e5;
        if nargin < 3,
            LB = -1e5;
        end
    end
end

% check input arguments

if ~ischar(FUN),
    error('''FUN'' incorrectly specified in ''SCE''');
end
if ~isfloat(X0),
    error('''X0'' incorrectly specified in ''SCE''');
end
if ~isfloat(LB),
    error('''LB'' incorrectly specified in ''SCE''');
end
if ~isfloat(UB),
    error('''UB'' incorrectly specified in ''SCE''');
end
if length(X0) ~= length(LB),
    error('''LB'' and ''X0'' have incompatible dimensions in ''SCE''');
end
if length(X0) ~= length(UB),
    error('''UB'' and ''X0'' have incompatible dimensions in ''SCE''');
end

% declaration of global variables

global NDIM nFUN_EVALS

% set EXITFLAG to default value

EXITFLAG = -2;

% determine number of variables to be optimized

NDIM = length(X0);

% seed the random number generator

rand('state',sum(100*clock));

% set default options

DEFAULT_OPTIONS = SCESET('nCOMPLEXES',5,...        % number of complexes
                         'nITER_INNER_LOOP',30,... % number of iteration in inner loop (CCE algorithm)
                         'MAX_ITER',2500,...       % maximum number of iterations
                         'MAX_TIME',2500,...       % maximum duration of optimization
                         'MAX_FUN_EVALS',2500,...  % maximum number of function evaluations
                         'TOLX',1e-3,...           % maximum difference between best and worst function evaluation in simplex
                         'TOLFUN',1e-3,...         % maximum difference between the coordinates of the vertices
                         'DISPLAY','none',...      % 'iter' or 'none' indicating whether user wants feedback
                         'OUTPUT_FCN',[]);         % string with output function name

% update default options with supplied options

OPTIONS = SCESET(DEFAULT_OPTIONS,OPTIONS);

% store OPTIONS in OUTPUT

OUTPUT.OPTIONS = OPTIONS;

% define number of points in each complex if not provided

OPTIONS.nPOINTS_COMPLEX = 2*NDIM+1;

% define number of points in each simplex if not provided

OPTIONS.nPOINTS_SIMPLEX = NDIM+1;

% define total number of points

nPOINTS = OPTIONS.nCOMPLEXES*OPTIONS.nPOINTS_COMPLEX;

% initialize population

POPULATION = nan(nPOINTS,NDIM,OPTIONS.MAX_ITER);

for i = 1:nPOINTS,
    if i == 1,
        POPULATION(i,:,1) = X0(:)';
    else
        POPULATION(i,:,1) = LB(:)'+rand(1,NDIM).*(UB(:)'-LB(:)');
    end
end

% initialize counters

nITERATIONS = 0;
nFUN_EVALS = 0;

% initialize timer

tic

% calculate cost for each point in initial population

POPULATION_FITNESS = nan(nPOINTS,OPTIONS.MAX_ITER);

for i = 1:nPOINTS;
    POPULATION_FITNESS(i,1) = CALCULATE_COST(FUN,POPULATION(i,:,1),LB,UB,varargin{:});
end

% sort the population in order of increasing function values

[POPULATION_FITNESS(:,1),idx] = sort(POPULATION_FITNESS(:,1));
POPULATION(:,:,1) = POPULATION(idx,:,1);

% for each iteration...

for i = 2:OPTIONS.MAX_ITER,
    
    % if a termination criterium was met, the value of EXITFLAG should have changed
    % from its default value of -2 to -1, 0, 1 or 2
    
    if EXITFLAG ~= -2,
        break
    end
    
    % add one to number of iterations counter
    
    nITERATIONS = nITERATIONS + 1;
    
    % The population matrix POPULATION will now be rearranged into so-called complexes.     
    % For each complex ...
    
    for j = 1:OPTIONS.nCOMPLEXES,
        
        % construct j-th complex from POPULATION
        
        k1 = 1:OPTIONS.nPOINTS_COMPLEX;
        k2 = (k1-1)*OPTIONS.nCOMPLEXES+j;
        
        COMPLEX(k1,:) = POPULATION(k2,:,i-1);
        COMPLEX_FITNESS(k1,1) = POPULATION_FITNESS(k2,i-1);
        
        % Each complex evolves a number of steps according to the competitive 
        % complex evolution (CCE) algorithm as described in Duan et al. (1992). 
        % Therefore, a number of 'parents' are selected from each complex which 
        % form a simplex. The selection of the parents is done so that the better
        % points in the complex have a higher probability to be selected as a 
        % parent. The paper of Duan et al. (1992) describes how a trapezoidal 
        % probability distribution can be used for this purpose. The implementation
        % found on the internet (implemented by Duan himself) however used a 
        % somewhat different probability distribution which is used here as well.
        
        for k = 1:OPTIONS.nITER_INNER_LOOP,