simpsa.m~

非线性单纯性和模拟退货算法的matlab程序： 本程序实现了非线性单纯性和模拟退货算法的联合优化求解

            disp(' Nr Iter  Nr Fun Eval    Min function       Best function        TEMP           Algorithm Step');
        else
            disp(sprintf('%5.0f      %5.0f       %12.6g     %15.6g      %12.6g       %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,'best point'));
        end
    end

    % run full metropolis cycle at current temperature
    % ------------------------------------------------
    
    % initialize vector COSTS, needed to calculate new temperature using cooling
    % schedule as described by Cardoso et al. (1996)
    COSTS = zeros((NDIM+1)*MAXITERTEMP,1);
    
    % initialize ITERTEMP to zero
    
    ITERTEMP = 0;
    
    % start

    for ITERTEMP = 1:MAXITERTEMP,
        
        % add one to number of iterations
        nITERATIONS = nITERATIONS + 1;
        
        % Press and Teukolsky (1991) add a positive logarithmic distributed variable,
        % proportional to the control temperature T to the function value associated with 
        % every vertex of the simplex. Likewise,they subtract a similar random variable 
        % from the function value at every new replacement point.
        % Thus, if the replacement point corresponds to a lower cost, this method always
        % accepts a true down hill step. If, on the other hand, the replacement point 
        % corresponds to a higher cost, an uphill move may be accepted, depending on the
        % relative COSTS of the perturbed values.
        % (taken from Cardoso et al.,1996)
        
        % add random fluctuations to function values of current vertices
        YFLUCT = Y+TEMP*abs(log(rand(NDIM+1,1)));
        
        % reorder YFLUCT, Y and P so that the first row corresponds to the lowest YFLUCT value
        help = sortrows([YFLUCT,Y,P],1);
        YFLUCT = help(:,1);
        Y = help(:,2);
        P = help(:,3:end);
        
        % store temperature at current iteration
        OUTPUT.TEMPERATURE(nITERATIONS) = TEMP;

        % store information about simplex at the current iteration
        OUTPUT.SIMPLEX(:,:,nITERATIONS) = P;
        OUTPUT.SIMPLEX_BEST(nITERATIONS,:) = PBEST;
        
        % store cost function value of best vertex in current iteration
        OUTPUT.COSTS(nITERATIONS,:) = Y;
        OUTPUT.COST_BEST(nITERATIONS) = YBEST;
        
        if strcmp(OPTIONS.DISPLAY,'iter'),
            disp(sprintf('%5.0f      %5.0f       %12.6g     %15.6g      %12.6g       %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,ALGOSTEP));
        end
        
        % if output function given then run output function to plot intermediate result
        if ~isempty(OPTIONS.OUTPUT_FCN),
            feval(OPTIONS.OUTPUT_FCN,P(1,:),Y(1));
        end
        
        % end the optimization if one of the stopping criteria is met
        %% 1. difference between best and worst function evaluation in simplex 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(Y)-min(Y)) < OPTIONS.TOLFUN) && (TEMP_LOOP_NUMBER ~= 1),
            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(P(2:NDIM+1,:)-P(1:NDIM,:)))) < OPTIONS.TOLX) && (TEMP_LOOP_NUMBER ~= 1),
            if strcmp(OPTIONS.DISPLAY,'iter'),
                disp('Change in X less than the specified tolerance (TOLX).')
            end
            EXITFLAG = 2;
            break;
        end
        
        if (nITERATIONS >= OPTIONS.MAX_ITER_TOTAL*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
        
        % begin a new iteration
        
        %% 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
        [YFTRY,YTRY,PTRY] = AMOTRY(FUN,P,-1,LB,UB,varargin{:});
        
        %% check the result
        if YFTRY <= YFLUCT(1),
            %% gives a result better than the best point,so try an additional
            %% extrapolation by a factor 2
            [YFTRYEXP,YTRYEXP,PTRYEXP] = AMOTRY(FUN,P,-2,LB,UB,varargin{:});
            if YFTRYEXP < YFTRY,
                P(end,:) = PTRYEXP;
                Y(end) = YTRYEXP;
                ALGOSTEP = 'reflection and expansion';
            else
                P(end,:) = PTRY;
                Y(end) = YTRY;
                ALGOSTEP = 'reflection';
            end
        elseif YFTRY >= YFLUCT(NDIM),
            %% the reflected point is worse than the second-highest, so look
            %% for an intermediate lower point, i.e., do a one-dimensional
            %% contraction
            [YFTRYCONTR,YTRYCONTR,PTRYCONTR] = AMOTRY(FUN,P,-0.5,LB,UB,varargin{:});
            if YFTRYCONTR < YFLUCT(end),
                P(end,:) = PTRYCONTR;
                Y(end) = YTRYCONTR;
                ALGOSTEP = 'one dimensional contraction';
            else
                %% can't seem to get rid of that high point, so better contract
                %% around the lowest (best) point
                X = ones(NDIM,NDIM)*diag(P(1,:));
                P(2:end,:) = 0.5*(P(2:end,:)+X);
                for k=2:NDIM,
                    Y(k) = CALCULATE_COST(FUN,P(k,:),LB,UB,varargin{:});
                end
                ALGOSTEP = 'multiple contraction';
            end
        else
            %% if YTRY better than second-highest point, use this point
            P(end,:) = PTRY;
            Y(end) = YTRY;
            ALGOSTEP = 'reflection';
        end
        
        % the initial temperature is estimated in the first loop from 
        % the number of successfull and unsuccesfull moves, and the average 
        % increase in cost associated with the unsuccessful moves
        
        if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START),
            if Y(1) > Y(end),
                SUCCESSFUL_MOVES = SUCCESSFUL_MOVES+1;
            elseif Y(1) <= Y(end),
                UNSUCCESSFUL_MOVES = UNSUCCESSFUL_MOVES+1;
                UNSUCCESSFUL_COSTS = UNSUCCESSFUL_COSTS+(Y(end)-Y(1));
            end
        end

    end

    % stop if previous for loop was broken due to some stop criterion
    if ITERTEMP < MAXITERTEMP,
        break;
    end
    
    % store cost function values in COSTS vector
    COSTS((ITERTEMP-1)*NDIM+1:ITERTEMP*NDIM+1) = Y;
    
    % calculated initial temperature or recalculate temperature 
    % using cooling schedule as proposed by Cardoso et al. (1996)
    % -----------------------------------------------------------
    
    if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START),
        TEMP = -(UNSUCCESSFUL_COSTS/(SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES))/log(((SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES)*OPTIONS.INITIAL_ACCEPTANCE_RATIO-SUCCESSFUL_MOVES)/UNSUCCESSFUL_MOVES);
    elseif TEMP_LOOP_NUMBER ~= 0,
        STDEV_Y = std(COSTS);
        COOLING_FACTOR = 1/(1+TEMP*log(1+OPTIONS.COOL_RATE)/(3*STDEV_Y));
        OLD_TEMP = TEMP;
        TEMP = TEMP*COOLING_FACTOR;
        if OLD_TEMP-TEMP < MIN_COOLING_FACTOR;
            TEMP = 0;
        end
    end
    
    % add one to temperature loop number
    TEMP_LOOP_NUMBER = TEMP_LOOP_NUMBER+1;
    
end

% return solution
X = PBEST;
FVAL = YBEST;

% store number of function evaluations
OUTPUT.nFUN_EVALS = nFUN_EVALS;

% store number of iterations
OUTPUT.nITERATIONS = nITERATIONS;

% trim OUTPUT data structure
OUTPUT.TEMPERATURE(nITERATIONS+1:end) = [];
OUTPUT.SIMPLEX(:,:,nITERATIONS+1:end) = [];
OUTPUT.SIMPLEX_BEST(nITERATIONS+1:end,:) = [];
OUTPUT.COSTS(nITERATIONS+1:end,:) = [];
OUTPUT.COST_BEST(nITERATIONS+1:end) = [];

% store the amount of time needed in OUTPUT data structure
OUTPUT.TIME = toc;

return

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

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

function [YFTRY,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 TEMP

% 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{:});
% substract random fluctuations to function values of current vertices
YFTRY = YTRY-TEMP*abs(log(rand(1)));

return

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

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

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

global YBEST PBEST NDIM nFUN_EVALS

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{:});

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

% save the best point ever
if YTRY < YBEST,
    YBEST = YTRY;
    PBEST = PTRY;
end

return