2013-10-02 15:29:56 +02:00
|
|
|
function [X,FVAL,EXITFLAG,OUTPUT] = simpsa(FUN,X0,LB,UB,OPTIONS,varargin)
|
|
|
|
|
2017-05-16 15:10:20 +02:00
|
|
|
% 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).
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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 is based on paper of Cardoso et al, 1996.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
|
|
|
% X=SIMPSA(FUN,X0) start at X0 and finds a minimum X to the function FUN.
|
2013-10-02 15:29:56 +02:00
|
|
|
% FUN accepts input X and returns a scalar function value F evaluated at X.
|
|
|
|
% X0 may be a scalar, vector, or matrix.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
|
|
|
% 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
|
2013-10-02 15:29:56 +02:00
|
|
|
% unbounded above.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
2013-10-02 15:29:56 +02:00
|
|
|
% X=SIMPSA(FUN,X0,LB,UB,OPTIONS) minimizes with the default optimization
|
2017-05-16 15:10:20 +02:00
|
|
|
% parameters replaced by values in the structure OPTIONS, an argument
|
|
|
|
% created with the SIMPSASET function. See SIMPSASET for details.
|
2013-10-02 15:29:56 +02:00
|
|
|
% Used options are TEMP_START, TEMP_END, COOL_RATE, INITIAL_ACCEPTANCE_RATIO,
|
|
|
|
% MIN_COOLING_FACTOR, 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.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
|
|
|
% X=SIMPSA(FUN,X0,LB,UB,OPTIONS,varargin) is used to supply a variable
|
2013-10-02 15:29:56 +02:00
|
|
|
% number of input arguments to the objective function FUN.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
|
|
|
% [X,FVAL]=SIMPSA(FUN,X0,...) returns the value of the objective
|
2013-10-02 15:29:56 +02:00
|
|
|
% function FUN at the solution X.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
|
|
|
% [X,FVAL,EXITFLAG]=SIMPSA(FUN,X0,...) returns an EXITFLAG that describes the
|
|
|
|
% exit condition of SIMPSA. Possible values of EXITFLAG and the corresponding
|
2013-10-02 15:29:56 +02:00
|
|
|
% exit conditions are:
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
|
|
|
% [X,FVAL,EXITFLAG,OUTPUT]=SIMPSA(FUN,X0,...) returns a structure OUTPUT with
|
2013-10-02 15:29:56 +02:00
|
|
|
% the number of iterations taken in OUTPUT.nITERATIONS, the number of function
|
|
|
|
% evaluations in OUTPUT.nFUN_EVALS, the temperature profile in OUTPUT.TEMPERATURE,
|
2017-05-16 15:10:20 +02:00
|
|
|
% 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
|
2013-10-02 15:29:56 +02:00
|
|
|
% needed in OUTPUT.TIME and the options used in OUTPUT.OPTIONS.
|
2017-05-16 15:10:20 +02:00
|
|
|
%
|
2013-10-02 15:29:56 +02:00
|
|
|
% See also SIMPSASET, SIMPSAGET
|
|
|
|
|
|
|
|
|
|
|
|
% Copyright (C) 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se
|
|
|
|
% Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be
|
2017-05-18 18:36:38 +02:00
|
|
|
% Copyright (C) 2013-2017 Dynare Team.
|
2013-10-02 15:29:56 +02:00
|
|
|
%
|
|
|
|
% This file is part of Dynare.
|
|
|
|
%
|
|
|
|
% Dynare 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 3 of the License, or
|
|
|
|
% (at your option) any later version.
|
|
|
|
%
|
|
|
|
% Dynare 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 Dynare. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
|
|
|
|
% handle variable input arguments
|
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if nargin < 5
|
2013-10-02 15:29:56 +02:00
|
|
|
OPTIONS = [];
|
2017-05-16 12:42:01 +02:00
|
|
|
if nargin < 4
|
2013-10-02 15:29:56 +02:00
|
|
|
UB = 1e5;
|
2017-05-16 12:42:01 +02:00
|
|
|
if nargin < 3
|
2013-10-02 15:29:56 +02:00
|
|
|
LB = -1e5;
|
|
|
|
end
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
% check input arguments
|
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if ~ischar(FUN)
|
2013-10-02 15:29:56 +02:00
|
|
|
error('''FUN'' incorrectly specified in ''SIMPSA''');
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
if ~isfloat(X0)
|
2013-10-02 15:29:56 +02:00
|
|
|
error('''X0'' incorrectly specified in ''SIMPSA''');
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
if ~isfloat(LB)
|
2013-10-02 15:29:56 +02:00
|
|
|
error('''LB'' incorrectly specified in ''SIMPSA''');
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
if ~isfloat(UB)
|
2013-10-02 15:29:56 +02:00
|
|
|
error('''UB'' incorrectly specified in ''SIMPSA''');
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
if length(X0) ~= length(LB)
|
2013-10-02 15:29:56 +02:00
|
|
|
error('''LB'' and ''X0'' have incompatible dimensions in ''SIMPSA''');
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
if length(X0) ~= length(UB)
|
2013-10-02 15:29:56 +02:00
|
|
|
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);
|
|
|
|
|
|
|
|
% set default options
|
|
|
|
DEFAULT_OPTIONS = simpsaset('TEMP_START',[],... % starting temperature (if none provided, an optimal one will be estimated)
|
2017-05-16 15:10:20 +02:00
|
|
|
'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
|
|
|
|
'MIN_COOLING_FACTOR',0.9,... % minimum cooling factor (<1)
|
|
|
|
'MAX_ITER_TEMP_FIRST',50,... % number of iterations in the preliminary temperature loop
|
|
|
|
'MAX_ITER_TEMP_LAST',2000,... % 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',20000,... % maximum number of function evaluations
|
|
|
|
'TOLX',1e-6,... % maximum difference between best and worst function evaluation in simplex
|
|
|
|
'TOLFUN',1e-6,... % maximum difference between the coordinates of the vertices
|
|
|
|
'DISPLAY','iter',... % 'iter' or 'none' indicating whether user wants feedback
|
|
|
|
'OUTPUT_FCN',[]); % string with output function name
|
2013-10-02 15:29:56 +02:00
|
|
|
|
|
|
|
% update default options with supplied options
|
|
|
|
|
|
|
|
OPTIONS = simpsaset(DEFAULT_OPTIONS,OPTIONS);
|
|
|
|
|
|
|
|
% store options in OUTPUT
|
2017-06-23 10:50:06 +02:00
|
|
|
if nargout>3
|
|
|
|
OUTPUT.OPTIONS = OPTIONS;
|
|
|
|
end
|
2013-10-02 15:29:56 +02:00
|
|
|
|
|
|
|
% 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
|
2017-05-16 12:42:01 +02:00
|
|
|
if isempty(OPTIONS.TEMP_START)
|
2013-10-02 15:29:56 +02:00
|
|
|
TEMP = abs(YBEST)*1e5;
|
|
|
|
else
|
|
|
|
TEMP = OPTIONS.TEMP_START;
|
|
|
|
end
|
|
|
|
|
|
|
|
% initialize OUTPUT structure
|
|
|
|
% ---------------------------
|
2017-06-23 10:50:06 +02:00
|
|
|
if nargout>3
|
|
|
|
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);
|
|
|
|
end
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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
|
|
|
|
% -----------------------------------------------------------
|
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
while 1
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% detect if termination criterium was met
|
|
|
|
% ---------------------------------------
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if EXITFLAG ~= -2
|
2013-10-02 15:29:56 +02:00
|
|
|
break
|
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% set MAXITERTEMP: maximum number of iterations at current temperature
|
|
|
|
% --------------------------------------------------------------------
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if TEMP_LOOP_NUMBER == 1
|
2013-10-02 15:29:56 +02:00
|
|
|
MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_FIRST*NDIM;
|
2017-05-16 15:10:20 +02:00
|
|
|
% 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
|
2013-10-02 15:29:56 +02:00
|
|
|
% for the unsuccessful moves.
|
2017-05-16 12:42:01 +02:00
|
|
|
if isempty(OPTIONS.TEMP_START)
|
2013-10-02 15:29:56 +02:00
|
|
|
[SUCCESSFUL_MOVES,UNSUCCESSFUL_MOVES,UNSUCCESSFUL_COSTS] = deal(0);
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
elseif TEMP < OPTIONS.TEMP_END
|
2013-10-02 15:29:56 +02:00
|
|
|
TEMP = 0;
|
|
|
|
MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_LAST*NDIM;
|
|
|
|
else
|
|
|
|
MAXITERTEMP = OPTIONS.MAX_ITER_TEMP*NDIM;
|
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% construct initial simplex
|
|
|
|
% -------------------------
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% 1st vertex of initial simplex
|
|
|
|
P(1,:) = PBEST;
|
|
|
|
Y(1) = CALCULATE_COST(FUN,P(1,:),LB,UB,varargin{:});
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% if output function given then run output function to plot intermediate result
|
2017-05-16 12:42:01 +02:00
|
|
|
if ~isempty(OPTIONS.OUTPUT_FCN)
|
2013-10-02 17:08:04 +02:00
|
|
|
feval(OPTIONS.OUTPUT_FCN,transpose(P(1,:)),Y(1));
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% remaining vertices of simplex
|
2017-05-16 12:42:01 +02:00
|
|
|
for k = 1:NDIM
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% store information on what step the algorithm just did
|
|
|
|
ALGOSTEP = 'initial simplex';
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% add NDIM+1 to number of function evaluations
|
|
|
|
nFUN_EVALS = nFUN_EVALS + NDIM;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% note:
|
|
|
|
% dimensions of matrix P: (NDIM+1) x NDIM
|
|
|
|
% dimensions of vector Y: (NDIM+1) x 1
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% give user feedback if requested
|
2017-05-16 12:42:01 +02:00
|
|
|
if strcmp(OPTIONS.DISPLAY,'iter')
|
|
|
|
if nITERATIONS == 0
|
2013-10-02 15:29:56 +02:00
|
|
|
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
|
|
|
|
% ------------------------------------------------
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% initialize vector COSTS, needed to calculate new temperature using cooling
|
|
|
|
% schedule as described by Cardoso et al. (1996)
|
|
|
|
COSTS = zeros((NDIM+1)*MAXITERTEMP,1);
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% initialize ITERTEMP to zero
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
ITERTEMP = 0;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% start
|
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
for ITERTEMP = 1:MAXITERTEMP
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% add one to number of iterations
|
|
|
|
nITERATIONS = nITERATIONS + 1;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% Press and Teukolsky (1991) add a positive logarithmic distributed variable,
|
2017-05-16 15:10:20 +02:00
|
|
|
% proportional to the control temperature T to the function value associated with
|
|
|
|
% every vertex of the simplex. Likewise,they subtract a similar random variable
|
2013-10-02 15:29:56 +02:00
|
|
|
% from the function value at every new replacement point.
|
|
|
|
% Thus, if the replacement point corresponds to a lower cost, this method always
|
2017-05-16 15:10:20 +02:00
|
|
|
% accepts a true down hill step. If, on the other hand, the replacement point
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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)
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% add random fluctuations to function values of current vertices
|
|
|
|
YFLUCT = Y+TEMP*abs(log(rand(NDIM+1,1)));
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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);
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-06-23 10:50:06 +02:00
|
|
|
if nargout>3
|
|
|
|
% store temperature at current iteration
|
|
|
|
OUTPUT.TEMPERATURE(nITERATIONS) = TEMP;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-06-23 10:50:06 +02:00
|
|
|
% store information about simplex at the current iteration
|
|
|
|
OUTPUT.SIMPLEX(:,:,nITERATIONS) = P;
|
|
|
|
OUTPUT.SIMPLEX_BEST(nITERATIONS,:) = PBEST;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-06-23 10:50:06 +02:00
|
|
|
% store cost function value of best vertex in current iteration
|
|
|
|
OUTPUT.COSTS(nITERATIONS,:) = Y;
|
|
|
|
OUTPUT.COST_BEST(nITERATIONS) = YBEST;
|
|
|
|
end
|
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if strcmp(OPTIONS.DISPLAY,'iter')
|
2013-10-02 15:29:56 +02:00
|
|
|
disp(sprintf('%5.0f %5.0f %12.6g %15.6g %12.6g %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,ALGOSTEP));
|
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% if output function given then run output function to plot intermediate result
|
2017-05-16 12:42:01 +02:00
|
|
|
if ~isempty(OPTIONS.OUTPUT_FCN)
|
2013-10-02 17:08:04 +02:00
|
|
|
feval(OPTIONS.OUTPUT_FCN,transpose(P(1,:)),Y(1));
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% end the optimization if one of the stopping criteria is met
|
2017-05-16 15:10:20 +02:00
|
|
|
%% 1. difference between best and worst function evaluation in simplex is smaller than TOLFUN
|
2013-10-02 15:29:56 +02:00
|
|
|
%% 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
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if (abs(max(Y)-min(Y)) < OPTIONS.TOLFUN) && (TEMP_LOOP_NUMBER ~= 1)
|
|
|
|
if strcmp(OPTIONS.DISPLAY,'iter')
|
2013-10-02 15:29:56 +02:00
|
|
|
disp('Change in the objective function value less than the specified tolerance (TOLFUN).')
|
|
|
|
end
|
|
|
|
EXITFLAG = 1;
|
2017-05-16 12:42:01 +02:00
|
|
|
break
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if (max(max(abs(P(2:NDIM+1,:)-P(1:NDIM,:)))) < OPTIONS.TOLX) && (TEMP_LOOP_NUMBER ~= 1)
|
|
|
|
if strcmp(OPTIONS.DISPLAY,'iter')
|
2013-10-02 15:29:56 +02:00
|
|
|
disp('Change in X less than the specified tolerance (TOLX).')
|
|
|
|
end
|
|
|
|
EXITFLAG = 2;
|
2017-05-16 12:42:01 +02:00
|
|
|
break
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if (nITERATIONS >= OPTIONS.MAX_ITER_TOTAL*NDIM) || (nFUN_EVALS >= OPTIONS.MAX_FUN_EVALS*NDIM*(NDIM+1))
|
|
|
|
if strcmp(OPTIONS.DISPLAY,'iter')
|
2013-10-02 15:29:56 +02:00
|
|
|
disp('Maximum number of function evaluations or iterations reached.');
|
|
|
|
end
|
|
|
|
EXITFLAG = 0;
|
2017-05-16 12:42:01 +02:00
|
|
|
break
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if toc/60 > OPTIONS.MAX_TIME
|
|
|
|
if strcmp(OPTIONS.DISPLAY,'iter')
|
2013-10-02 15:29:56 +02:00
|
|
|
disp('Exceeded maximum time.');
|
|
|
|
end
|
|
|
|
EXITFLAG = -1;
|
2017-05-16 12:42:01 +02:00
|
|
|
break
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% begin a new iteration
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
%% 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{:});
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
%% check the result
|
2017-05-16 12:42:01 +02:00
|
|
|
if YFTRY <= YFLUCT(1)
|
2013-10-02 15:29:56 +02:00
|
|
|
%% 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{:});
|
2017-05-16 12:42:01 +02:00
|
|
|
if YFTRYEXP < YFTRY
|
2013-10-02 15:29:56 +02:00
|
|
|
P(end,:) = PTRYEXP;
|
|
|
|
Y(end) = YTRYEXP;
|
|
|
|
ALGOSTEP = 'reflection and expansion';
|
|
|
|
else
|
|
|
|
P(end,:) = PTRY;
|
|
|
|
Y(end) = YTRY;
|
|
|
|
ALGOSTEP = 'reflection';
|
|
|
|
end
|
2017-05-16 12:42:01 +02:00
|
|
|
elseif YFTRY >= YFLUCT(NDIM)
|
2013-10-02 15:29:56 +02:00
|
|
|
%% 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{:});
|
2017-05-16 12:42:01 +02:00
|
|
|
if YFTRYCONTR < YFLUCT(end)
|
2013-10-02 15:29:56 +02:00
|
|
|
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);
|
2017-05-16 12:42:01 +02:00
|
|
|
for k=2:NDIM
|
2013-10-02 15:29:56 +02:00
|
|
|
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
|
2017-05-16 15:10:20 +02:00
|
|
|
|
|
|
|
% the initial temperature is estimated in the first loop from
|
|
|
|
% the number of successfull and unsuccesfull moves, and the average
|
2013-10-02 15:29:56 +02:00
|
|
|
% increase in cost associated with the unsuccessful moves
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START)
|
|
|
|
if Y(1) > Y(end)
|
2013-10-02 15:29:56 +02:00
|
|
|
SUCCESSFUL_MOVES = SUCCESSFUL_MOVES+1;
|
2017-05-16 12:42:01 +02:00
|
|
|
elseif Y(1) <= Y(end)
|
2013-10-02 15:29:56 +02:00
|
|
|
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
|
2017-05-16 12:42:01 +02:00
|
|
|
if ITERTEMP < MAXITERTEMP
|
|
|
|
break
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% store cost function values in COSTS vector
|
|
|
|
COSTS((ITERTEMP-1)*NDIM+1:ITERTEMP*NDIM+1) = Y;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
|
|
|
% calculated initial temperature or recalculate temperature
|
2013-10-02 15:29:56 +02:00
|
|
|
% using cooling schedule as proposed by Cardoso et al. (1996)
|
|
|
|
% -----------------------------------------------------------
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START)
|
2013-10-02 15:29:56 +02:00
|
|
|
TEMP = -(UNSUCCESSFUL_COSTS/(SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES))/log(((SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES)*OPTIONS.INITIAL_ACCEPTANCE_RATIO-SUCCESSFUL_MOVES)/UNSUCCESSFUL_MOVES);
|
2017-05-16 12:42:01 +02:00
|
|
|
elseif TEMP_LOOP_NUMBER ~= 0
|
2013-10-02 15:29:56 +02:00
|
|
|
STDEV_Y = std(COSTS);
|
|
|
|
COOLING_FACTOR = 1/(1+TEMP*log(1+OPTIONS.COOL_RATE)/(3*STDEV_Y));
|
|
|
|
TEMP = TEMP*min(OPTIONS.MIN_COOLING_FACTOR,COOLING_FACTOR);
|
|
|
|
end
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
% add one to temperature loop number
|
|
|
|
TEMP_LOOP_NUMBER = TEMP_LOOP_NUMBER+1;
|
2017-05-16 15:10:20 +02:00
|
|
|
|
2013-10-02 15:29:56 +02:00
|
|
|
end
|
|
|
|
|
|
|
|
% return solution
|
2013-10-02 17:08:04 +02:00
|
|
|
X = transpose(PBEST);
|
2013-10-02 15:29:56 +02:00
|
|
|
FVAL = YBEST;
|
|
|
|
|
2017-06-23 10:50:06 +02:00
|
|
|
if nargout>3
|
|
|
|
% 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;
|
|
|
|
end
|
2013-10-02 15:29:56 +02:00
|
|
|
|
|
|
|
return
|
|
|
|
|
|
|
|
% ==============================================================================
|
|
|
|
|
|
|
|
% AMOTRY FUNCTION
|
|
|
|
% ---------------
|
|
|
|
|
|
|
|
function [YFTRY,YTRY,PTRY] = AMOTRY(FUN,P,fac,LB,UB,varargin)
|
2017-05-16 15:10:20 +02:00
|
|
|
% 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
|
2013-10-02 15:29:56 +02:00
|
|
|
% 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
|
|
|
|
|
2017-05-16 12:42:01 +02:00
|
|
|
for i = 1:NDIM
|
2013-10-02 15:29:56 +02:00
|
|
|
% check lower bounds
|
2017-05-16 12:42:01 +02:00
|
|
|
if PTRY(i) < LB(i)
|
2013-10-02 15:29:56 +02:00
|
|
|
YTRY = 1e12+(LB(i)-PTRY(i))*1e6;
|
|
|
|
return
|
|
|
|
end
|
|
|
|
% check upper bounds
|
2017-05-16 12:42:01 +02:00
|
|
|
if PTRY(i) > UB(i)
|
2013-10-02 15:29:56 +02:00
|
|
|
YTRY = 1e12+(PTRY(i)-UB(i))*1e6;
|
|
|
|
return
|
|
|
|
end
|
|
|
|
end
|
|
|
|
|
|
|
|
% calculate cost associated with PTRY
|
2013-10-02 17:08:04 +02:00
|
|
|
YTRY = feval(FUN,PTRY(:),varargin{:});
|
2013-10-02 15:29:56 +02:00
|
|
|
|
|
|
|
% add one to number of function evaluations
|
|
|
|
nFUN_EVALS = nFUN_EVALS + 1;
|
|
|
|
|
|
|
|
% save the best point ever
|
2017-05-16 12:42:01 +02:00
|
|
|
if YTRY < YBEST
|
2013-10-02 15:29:56 +02:00
|
|
|
YBEST = YTRY;
|
|
|
|
PBEST = PTRY;
|
|
|
|
end
|
|
|
|
|
|
|
|
return
|