sce.m

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

            
            % select simplex by sampling the complex
            
            LOCATION(1) = 1; % the LOCATION of the selected point in the complex

            for l = 2:OPTIONS.nPOINTS_SIMPLEX;
                ALREADY_PARENT = 0;
                while ~ALREADY_PARENT,
                    DUMMY = 1 + floor(OPTIONS.nPOINTS_COMPLEX+0.5-sqrt((OPTIONS.nPOINTS_COMPLEX+0.5)^2 - OPTIONS.nPOINTS_COMPLEX*(OPTIONS.nPOINTS_COMPLEX+1)*rand));
                    if isempty(find(LOCATION(1:l-1) == DUMMY,1,'first')),
                        ALREADY_PARENT = 1;
                    end
                end
                LOCATION(l) = DUMMY;
            end
            
            LOCATION = sort(LOCATION);
            
            % construct the simplex
            
            SIMPLEX = COMPLEX(LOCATION,:);
            SIMPLEX_FITNESS = COMPLEX_FITNESS(LOCATION);
            
            % generate new point for simplex
            
            %% first extrapolate by a factor -1 through the face of the simplex
            %% across from the high point,i.e.,reflect the simplex from the high point
            [SFTRY,SXTRY] = AMOTRY(FUN,SIMPLEX,-1,LB,UB,varargin{:});
            
            %% check the result
            if SFTRY <= SIMPLEX_FITNESS(1),
                %% gives a result better than the best point,so try an additional
                %% extrapolation by a factor 2
                [SFTRYEXP,SXTRYEXP] = AMOTRY(FUN,SIMPLEX,-2,LB,UB,varargin{:});
                if SFTRYEXP < SFTRY,
                    SIMPLEX(end,:) = SXTRYEXP;
                    SIMPLEX_FITNESS(end) = SFTRYEXP;
                    ALGOSTEP = 'reflection and expansion';
                else
                    SIMPLEX(end,:) = SXTRY;
                    SIMPLEX_FITNESS(end) = SFTRY;
                    ALGOSTEP = 'reflection';
                end
            elseif SFTRY >= SIMPLEX_FITNESS(NDIM),
                %% the reflected point is worse than the second-highest, so look
                %% for an intermediate lower point, i.e., do a one-dimensional
                %% contraction
                [SFTRYCONTR,SXTRYCONTR] = AMOTRY(FUN,SIMPLEX,-0.5,LB,UB,varargin{:});
                if SFTRYCONTR < SIMPLEX_FITNESS(end),
                    SIMPLEX(end,:) = SXTRYCONTR;
                    SIMPLEX_FITNESS(end) = SFTRYCONTR;
                    ALGOSTEP = 'one dimensional contraction';
                else
                    %% can't seem to get rid of that high point, so better contract
                    %% around the lowest (best) point
                    SX_HELP = ones(NDIM,NDIM)*diag(SIMPLEX(1,:));
                    SIMPLEX(2:end,:) = 0.5*(SIMPLEX(2:end,:)+SX_HELP);
                    for k=2:NDIM,
                        SIMPLEX_FITNESS(k) = CALCULATE_COST(FUN,SIMPLEX(k,:),LB,UB,varargin{:});
                    end
                    ALGOSTEP = 'multiple contraction';
                end
            else
                %% if ytry better than second-highest point, use this point
                SIMPLEX(end,:) = SXTRY;
                SIMPLEX_FITNESS(end) = SFTRY;
                ALGOSTEP = 'reflection';
            end
            
            % replace the simplex into the complex
            
            COMPLEX(LOCATION,:) = SIMPLEX;
            COMPLEX_FITNESS(LOCATION) = SIMPLEX_FITNESS;
            
            % sort the complex;
            
            [COMPLEX_FITNESS,idx] = sort(COMPLEX_FITNESS);
            COMPLEX = COMPLEX(idx,:);
            
        end
        
        % replace the complex back into the population
        
        POPULATION(k2,:,i) = COMPLEX(k1,:);
        POPULATION_FITNESS(k2,i) = COMPLEX_FITNESS(k1);
        
    end
    
    % At this point, the population was divided in several complexes, each of which
    % underwent a number of iteration of the simplex (Metropolis) algorithm. Now,
    % the points in the population are sorted, the termination criteria are checked
    % and output is given on the screen if requested.
    
    % sort the population
    
    [POPULATION_FITNESS(:,i),idx] = sort(POPULATION_FITNESS(:,i));
    POPULATION(:,:,i) = POPULATION(idx,:,i);
    
    % give user feedback on screen if requested
    
    if strcmp(OPTIONS.DISPLAY,'iter'),
        if nITERATIONS == 1,
            disp(' Nr Iter  Nr Fun Eval    Current best function    Current worst function       Best function');
            disp(sprintf(' %5.0f     %5.0f             %12.6g              %12.6g           %15.6g',nITERATIONS,nFUN_EVALS,min(POPULATION_FITNESS(:,i)),max(POPULATION_FITNESS(:,i)),min(min(POPULATION_FITNESS))));
        else
            disp(sprintf(' %5.0f     %5.0f             %12.6g              %12.6g           %15.6g',nITERATIONS,nFUN_EVALS,min(POPULATION_FITNESS(:,i)),max(POPULATION_FITNESS(:,i)),min(min(POPULATION_FITNESS))));
        end
    end
    
    % end the optimization if one of the stopping criteria is met
    %% 1. difference between best and worst function evaluation in population is smaller than TOLFUN 
    %% 2. maximum difference between the coordinates of the vertices in simplex is less than TOLX
    %% 3. no convergence,but maximum number of iterations has been reached
    %% 4. no convergence,but maximum time has been reached
    
    if abs(max(POPULATION_FITNESS(:,i))-min(POPULATION_FITNESS(:,i))) < OPTIONS.TOLFUN,
        if strcmp(OPTIONS.DISPLAY,'iter'),
            disp('Change in the objective function value less than the specified tolerance (TOLFUN).')
        end
        EXITFLAG = 1;
        break;
    end
    
    if max(max(abs(diff(POPULATION(:,:,i),1,1)))) < OPTIONS.TOLX,
        if strcmp(OPTIONS.DISPLAY,'iter'),
            disp('Change in X less than the specified tolerance (TOLX).')
        end
        EXITFLAG = 2;
        break;
    end
    
    if (i >= OPTIONS.MAX_ITER*NDIM) || (nFUN_EVALS >= OPTIONS.MAX_FUN_EVALS*NDIM*(NDIM+1)),
        if strcmp(OPTIONS.DISPLAY,'iter'),
            disp('Maximum number of function evaluations or iterations reached.');
        end
        EXITFLAG = 0;
        break;
    end
    
    if toc/60 > OPTIONS.MAX_TIME,
        if strcmp(OPTIONS.DISPLAY,'iter'),
            disp('Exceeded maximum time.');
        end
        EXITFLAG = -1;
        break;
    end
    
end

% return solution

X = POPULATION(1,:,i);
FVAL = POPULATION_FITNESS(1,i);

% store number of function evaluations

OUTPUT.nFUN_EVALS = nFUN_EVALS;

% store number of iterations

OUTPUT.nITERATIONS = nITERATIONS;

% store information on the population at each iteration

OUTPUT.POPULATION = POPULATION(:,:,1:nITERATIONS);
OUTPUT.POPULATION_FITNESS = POPULATION_FITNESS(:,1:nITERATIONS);

% store the amount of time needed in OUTPUT data structure

OUTPUT.TIME = toc;
return

% ==============================================================================

% AMOTRY FUNCTION
% ---------------

function [YTRY,PTRY] = AMOTRY(FUN,P,FAC,LB,UB,varargin)
% Extrapolates by a factor FAC through the face of the simplex across from 
% the high point, tries it, and replaces the high point if the new point is 
% better.

global NDIM

% calculate coordinates of new vertex
PSUM = sum(P(1:NDIM,:))/NDIM;
PTRY = PSUM*(1-FAC)+P(end,:)*FAC;

% evaluate the function at the trial point.
YTRY = CALCULATE_COST(FUN,PTRY,LB,UB,varargin{:});

return

% ==============================================================================

% COST FUNCTION EVALUATION
% ------------------------

function [YTRY] = CALCULATE_COST(FUN,PTRY,LB,UB,varargin)

global NDIM nFUN_EVALS

% add one to number of function evaluations
nFUN_EVALS = nFUN_EVALS + 1;

for i = 1:NDIM,
    % check lower bounds
    if PTRY(i) < LB(i),
        YTRY = 1e12+(LB(i)-PTRY(i))*1e6;
        return
    end
    % check upper bounds
    if PTRY(i) > UB(i),
        YTRY = 1e12+(PTRY(i)-UB(i))*1e6;
        return
    end
end

% calculate cost associated with PTRY
YTRY = feval(FUN,PTRY,varargin{:});

return