simpsa.m~

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

function [X,FVAL,EXITFLAG,OUTPUT] = SIMPSA(FUN,X0,LB,UB,OPTIONS,varargin)
%SIMPSA finds a minimum of a function of several variables using an algorithm 
% that is based on the combination of the non-linear smplex and the simulated 
% annealing algorithm (the SIMPSA algorithm, Cardoso et al., 1996). 
% In this paper, the algorithm is shown to be adequate for the global optimi-
% zation of an example set of unconstrained and constrained NLP functions.
%
%   SIMPSA attempts to solve problems of the form:
%       min F(X) subject to: LB <= X <= UB
%        X
%
% Algorithm partly based on section 10.4 and 10.9 in "Numerical Recipes in C",
% ISBN 0-521-43108-5, and the paper of Cardoso et al, 1996.
%                                                                             
%   X=SIMPSA(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=SIMPSA(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=SIMPSA(FUN,X0,LB,UB,OPTIONS) minimizes with the default optimization
%   parameters replaced by values in the structure OPTIONS, an argument 
%   created with the SIMPSASET function. See SIMPSASET for details. 
%   Used options are TEMP_START, TEMP_END, COOL_RATE, INITIAL_ACCEPTANCE_RATIO, 
%   MAX_ITER_TEMP_FIRST, MAX_ITER_TEMP_LAST, MAX_ITER_TEMP, MAX_ITER_TOTAL, MAX_TIME,
%   MAX_FUN_EVALS, TOLX, TOLFUN, DISPLAY and OUTPUT_FCN.
%   Use OPTIONS = [] as a place holder if no options are set.
%   
%   X=SIMPSA(FUN,X0,LB,UB,OPTIONS,varargin) is used to supply a variable 
%   number of input arguments to the objective function FUN.
%   
%   [X,FVAL]=SIMPSA(FUN,X0,...) returns the value of the objective 
%   function FUN at the solution X.
%   
%   [X,FVAL,EXITFLAG]=SIMPSA(FUN,X0,...) returns an EXITFLAG that describes the 
%   exit condition of SIMPSA. 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]=SIMPSA(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 temperature profile in OUTPUT.TEMPERATURE,
%   the simplexes that were evaluated in OUTPUT.SIMPLEX and the best one in 
%   OUTPUT.SIMPLEX_BEST, the costs associated with each simplex in OUTPUT.COSTS and 
%   from the best simplex at that iteration in OUTPUT.COST_BEST, the amount of time 
%   needed in OUTPUT.TIME and the options used in OUTPUT.OPTIONS.
% 
%   See also SIMPSASET, SIMPSAGET



% Copyright (C) 2006 Brecht Donckels, BIOMATH, [e-mail]
% 
% inspired by:
% Systems Biology Toolbox for MATLAB
% Copyright (C) 2005 Henning Schmidt, FCC, [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 ''SIMPSA''');
end
if ~isfloat(X0),
    error('''X0'' incorrectly specified in ''SIMPSA''');
end
if ~isfloat(LB),
    error('''LB'' incorrectly specified in ''SIMPSA''');
end
if ~isfloat(UB),
    error('''UB'' incorrectly specified in ''SIMPSA''');
end
if length(X0) ~= length(LB),
    error('''LB'' and ''X0'' have incompatible dimensions in ''SIMPSA''');
end
if length(X0) ~= length(UB),
    error('''UB'' and ''X0'' have incompatible dimensions in ''SIMPSA''');
end

% declaration of global variables

global NDIM nFUN_EVALS TEMP YBEST PBEST

% 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 = SIMPSASET('TEMP_START',[],...  % starting temperature (if none provided, an optimal one will be estimated)
             'TEMP_END',1,...                    % end temperature
             'COOL_RATE',10,...                  % small values (<1) means slow convergence,large values (>1) means fast convergence
             'INITIAL_ACCEPTANCE_RATIO',0.95,... % when initial temperature is estimated, this will be the initial acceptance ratio in the first round
             'MAX_ITER_TEMP_FIRST',50,...        % number of iterations in the preliminary temperature loop
             'MAX_ITER_TEMP_LAST',50,...         % number of iterations in the last temperature loop (pure simplex)
             'MAX_ITER_TEMP',10,...              % number of iterations in the remaining temperature loops
             'MAX_ITER_TOTAL',2500,...           % maximum number of iterations tout court
             'MAX_TIME',2500,...                 % maximum duration of optimization
             'MAX_FUN_EVALS',2500,...            % maximum number of function evaluations
             'TOLX',1e-6,...                     % 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 = SIMPSASET(DEFAULT_OPTIONS,OPTIONS);

% store options in OUTPUT

OUTPUT.OPTIONS = OPTIONS;

% initialize simplex
% ------------------

% create empty simplex matrix p (location of vertex i in row i)
P = zeros(NDIM+1,NDIM);
% create empty cost vector (cost of vertex i in row i)
Y = zeros(NDIM+1,1);
% set best vertex of initial simplex equal to initial parameter guess
PBEST = X0(:)';
% calculate cost with best vertex of initial simplex
YBEST = CALCULATE_COST(FUN,PBEST,LB,UB,varargin{:});

% initialize temperature loop
% ---------------------------

% set temperature loop number to one
TEMP_LOOP_NUMBER = 1;

% if no TEMP_START is supplied, the initial temperature is estimated in the first
% loop as described by Cardoso et al., 1996 (recommended)

% therefore, the temperature is set to YBEST*1e5 in the first loop
if isempty(OPTIONS.TEMP_START),
    TEMP = abs(YBEST)*1e5;
else
    TEMP = OPTIONS.TEMP_START;
end

% initialize OUTPUT structure
% ---------------------------

OUTPUT.TEMPERATURE = zeros(OPTIONS.MAX_ITER_TOTAL,1);
OUTPUT.SIMPLEX = zeros(NDIM+1,NDIM,OPTIONS.MAX_ITER_TOTAL);
OUTPUT.SIMPLEX_BEST = zeros(OPTIONS.MAX_ITER_TOTAL,NDIM);
OUTPUT.COSTS = zeros(OPTIONS.MAX_ITER_TOTAL,NDIM+1);
OUTPUT.COST_BEST = zeros(OPTIONS.MAX_ITER_TOTAL,1);

% initialize iteration data
% -------------------------

% start timer
tic
% set number of function evaluations to one
nFUN_EVALS = 1;
% set number of iterations to zero
nITERATIONS = 0;

% temperature loop: run SIMPSA till stopping criterion is met
% -----------------------------------------------------------

while 1,
    
    % detect if termination criterium was met
    % ---------------------------------------
    
    % 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
    
    % set MAXITERTEMP: maximum number of iterations at current temperature
    % --------------------------------------------------------------------
    
    if TEMP_LOOP_NUMBER == 1,
        MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_FIRST*NDIM;
        % The initial temperature is estimated (is requested) as described in 
        % Cardoso et al. (1996). Therefore, we need to store the number of 
        % successful and unsuccessful moves, as well as the increase in cost 
        % for the unsuccessful moves.
        if isempty(OPTIONS.TEMP_START),
            [SUCCESSFUL_MOVES,UNSUCCESSFUL_MOVES,UNSUCCESSFUL_COSTS] = deal(0);
        end
    elseif TEMP < OPTIONS.TEMP_END,
        TEMP = 0;
        MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_LAST*NDIM;
    else
        MAXITERTEMP = OPTIONS.MAX_ITER_TEMP*NDIM;
    end
    
    % construct initial simplex
    % -------------------------
    
    % 1st vertex of initial simplex
    P(1,:) = PBEST;
    Y(1) = CALCULATE_COST(FUN,P(1,:),LB,UB,varargin{:});
    
    % 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
    
    % remaining vertices of simplex
    for k = 1:NDIM,
        % copy first vertex in new vertex
        P(k+1,:) = P(1,:);
        % alter new vertex
        P(k+1,k) = LB(k)+rand*(UB(k)-LB(k));
        % calculate value of objective function at new vertex
        Y(k+1) = CALCULATE_COST(FUN,P(k+1,:),LB,UB,varargin{:});
    end
    
    % store information on what step the algorithm just did
    ALGOSTEP = 'initial simplex';
    
    % add NDIM+1 to number of function evaluations
    nFUN_EVALS = nFUN_EVALS + NDIM;
    
    % note:
    %  dimensions of matrix P: (NDIM+1) x NDIM
    %  dimensions of vector Y: (NDIM+1) x 1
    
    % give user feedback if requested
    if strcmp(OPTIONS.DISPLAY,'iter'),
        if nITERATIONS == 0,