1function [X,FVAL,EXITFLAG,OUTPUT] = simpsa(FUN,X0,LB,UB,OPTIONS,varargin) 2 3% Finds a minimum of a function of several variables using an algorithm 4% that is based on the combination of the non-linear smplex and the simulated 5% annealing algorithm (the SIMPSA algorithm, Cardoso et al., 1996). 6% In this paper, the algorithm is shown to be adequate for the global optimi- 7% zation of an example set of unconstrained and constrained NLP functions. 8% 9% SIMPSA attempts to solve problems of the form: 10% min F(X) subject to: LB <= X <= UB 11% X 12% 13% Algorithm partly is based on paper of Cardoso et al, 1996. 14% 15% X=SIMPSA(FUN,X0) start at X0 and finds a minimum X to the function FUN. 16% FUN accepts input X and returns a scalar function value F evaluated at X. 17% X0 may be a scalar, vector, or matrix. 18% 19% X=SIMPSA(FUN,X0,LB,UB) defines a set of lower and upper bounds on the 20% design variables, X, so that a solution is found in the range 21% LB <= X <= UB. Use empty matrices for LB and UB if no bounds exist. 22% Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is 23% unbounded above. 24% 25% X=SIMPSA(FUN,X0,LB,UB,OPTIONS) minimizes with the default optimization 26% parameters replaced by values in the structure OPTIONS, an argument 27% created with the SIMPSASET function. See SIMPSASET for details. 28% Used options are TEMP_START, TEMP_END, COOL_RATE, INITIAL_ACCEPTANCE_RATIO, 29% MIN_COOLING_FACTOR, MAX_ITER_TEMP_FIRST, MAX_ITER_TEMP_LAST, MAX_ITER_TEMP, 30% MAX_ITER_TOTAL, MAX_TIME, MAX_FUN_EVALS, TOLX, TOLFUN, DISPLAY and OUTPUT_FCN. 31% Use OPTIONS = [] as a place holder if no options are set. 32% 33% X=SIMPSA(FUN,X0,LB,UB,OPTIONS,varargin) is used to supply a variable 34% number of input arguments to the objective function FUN. 35% 36% [X,FVAL]=SIMPSA(FUN,X0,...) returns the value of the objective 37% function FUN at the solution X. 38% 39% [X,FVAL,EXITFLAG]=SIMPSA(FUN,X0,...) returns an EXITFLAG that describes the 40% exit condition of SIMPSA. Possible values of EXITFLAG and the corresponding 41% exit conditions are: 42% 43% 1 Change in the objective function value less than the specified tolerance. 44% 2 Change in X less than the specified tolerance. 45% 0 Maximum number of function evaluations or iterations reached. 46% -1 Maximum time exceeded. 47% 48% [X,FVAL,EXITFLAG,OUTPUT]=SIMPSA(FUN,X0,...) returns a structure OUTPUT with 49% the number of iterations taken in OUTPUT.nITERATIONS, the number of function 50% evaluations in OUTPUT.nFUN_EVALS, the temperature profile in OUTPUT.TEMPERATURE, 51% the simplexes that were evaluated in OUTPUT.SIMPLEX and the best one in 52% OUTPUT.SIMPLEX_BEST, the costs associated with each simplex in OUTPUT.COSTS and 53% from the best simplex at that iteration in OUTPUT.COST_BEST, the amount of time 54% needed in OUTPUT.TIME and the options used in OUTPUT.OPTIONS. 55% 56% See also SIMPSASET, SIMPSAGET 57 58 59% Copyright (C) 2005 Henning Schmidt, FCC, henning@fcc.chalmers.se 60% Copyright (C) 2006 Brecht Donckels, BIOMATH, brecht.donckels@ugent.be 61% Copyright (C) 2013-2017 Dynare Team. 62% 63% This file is part of Dynare. 64% 65% Dynare is free software: you can redistribute it and/or modify 66% it under the terms of the GNU General Public License as published by 67% the Free Software Foundation, either version 3 of the License, or 68% (at your option) any later version. 69% 70% Dynare is distributed in the hope that it will be useful, 71% but WITHOUT ANY WARRANTY; without even the implied warranty of 72% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 73% GNU General Public License for more details. 74% 75% You should have received a copy of the GNU General Public License 76% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 77 78% handle variable input arguments 79 80if nargin < 5 81 OPTIONS = []; 82 if nargin < 4 83 UB = 1e5; 84 if nargin < 3 85 LB = -1e5; 86 end 87 end 88end 89 90% check input arguments 91 92if ~ischar(FUN) 93 error('''FUN'' incorrectly specified in ''SIMPSA'''); 94end 95if ~isfloat(X0) 96 error('''X0'' incorrectly specified in ''SIMPSA'''); 97end 98if ~isfloat(LB) 99 error('''LB'' incorrectly specified in ''SIMPSA'''); 100end 101if ~isfloat(UB) 102 error('''UB'' incorrectly specified in ''SIMPSA'''); 103end 104if length(X0) ~= length(LB) 105 error('''LB'' and ''X0'' have incompatible dimensions in ''SIMPSA'''); 106end 107if length(X0) ~= length(UB) 108 error('''UB'' and ''X0'' have incompatible dimensions in ''SIMPSA'''); 109end 110 111% declaration of global variables 112 113global NDIM nFUN_EVALS TEMP YBEST PBEST 114 115% set EXITFLAG to default value 116 117EXITFLAG = -2; 118 119% determine number of variables to be optimized 120 121NDIM = length(X0); 122 123% set default options 124DEFAULT_OPTIONS = simpsaset('TEMP_START',[],... % starting temperature (if none provided, an optimal one will be estimated) 125 'TEMP_END',.1,... % end temperature 126 'COOL_RATE',10,... % small values (<1) means slow convergence,large values (>1) means fast convergence 127 'INITIAL_ACCEPTANCE_RATIO',0.95,... % when initial temperature is estimated, this will be the initial acceptance ratio in the first round 128 'MIN_COOLING_FACTOR',0.9,... % minimum cooling factor (<1) 129 'MAX_ITER_TEMP_FIRST',50,... % number of iterations in the preliminary temperature loop 130 'MAX_ITER_TEMP_LAST',2000,... % number of iterations in the last temperature loop (pure simplex) 131 'MAX_ITER_TEMP',10,... % number of iterations in the remaining temperature loops 132 'MAX_ITER_TOTAL',2500,... % maximum number of iterations tout court 133 'MAX_TIME',2500,... % maximum duration of optimization 134 'MAX_FUN_EVALS',20000,... % maximum number of function evaluations 135 'TOLX',1e-6,... % maximum difference between best and worst function evaluation in simplex 136 'TOLFUN',1e-6,... % maximum difference between the coordinates of the vertices 137 'DISPLAY','iter',... % 'iter' or 'none' indicating whether user wants feedback 138 'OUTPUT_FCN',[]); % string with output function name 139 140% update default options with supplied options 141 142OPTIONS = simpsaset(DEFAULT_OPTIONS,OPTIONS); 143 144% store options in OUTPUT 145if nargout>3 146 OUTPUT.OPTIONS = OPTIONS; 147end 148 149% initialize simplex 150% ------------------ 151 152% create empty simplex matrix p (location of vertex i in row i) 153P = zeros(NDIM+1,NDIM); 154% create empty cost vector (cost of vertex i in row i) 155Y = zeros(NDIM+1,1); 156% set best vertex of initial simplex equal to initial parameter guess 157PBEST = X0(:)'; 158% calculate cost with best vertex of initial simplex 159YBEST = CALCULATE_COST(FUN,PBEST,LB,UB,varargin{:}); 160 161% initialize temperature loop 162% --------------------------- 163 164% set temperature loop number to one 165TEMP_LOOP_NUMBER = 1; 166 167% if no TEMP_START is supplied, the initial temperature is estimated in the first 168% loop as described by Cardoso et al., 1996 (recommended) 169 170% therefore, the temperature is set to YBEST*1e5 in the first loop 171if isempty(OPTIONS.TEMP_START) 172 TEMP = abs(YBEST)*1e5; 173else 174 TEMP = OPTIONS.TEMP_START; 175end 176 177% initialize OUTPUT structure 178% --------------------------- 179if nargout>3 180 OUTPUT.TEMPERATURE = zeros(OPTIONS.MAX_ITER_TOTAL,1); 181 OUTPUT.SIMPLEX = zeros(NDIM+1,NDIM,OPTIONS.MAX_ITER_TOTAL); 182 OUTPUT.SIMPLEX_BEST = zeros(OPTIONS.MAX_ITER_TOTAL,NDIM); 183 OUTPUT.COSTS = zeros(OPTIONS.MAX_ITER_TOTAL,NDIM+1); 184 OUTPUT.COST_BEST = zeros(OPTIONS.MAX_ITER_TOTAL,1); 185end 186% initialize iteration data 187% ------------------------- 188 189% start timer 190tic 191% set number of function evaluations to one 192nFUN_EVALS = 1; 193% set number of iterations to zero 194nITERATIONS = 0; 195 196% temperature loop: run SIMPSA till stopping criterion is met 197% ----------------------------------------------------------- 198 199while 1 200 201 % detect if termination criterium was met 202 % --------------------------------------- 203 204 % if a termination criterium was met, the value of EXITFLAG should have changed 205 % from its default value of -2 to -1, 0, 1 or 2 206 207 if EXITFLAG ~= -2 208 break 209 end 210 211 % set MAXITERTEMP: maximum number of iterations at current temperature 212 % -------------------------------------------------------------------- 213 214 if TEMP_LOOP_NUMBER == 1 215 MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_FIRST*NDIM; 216 % The initial temperature is estimated (is requested) as described in 217 % Cardoso et al. (1996). Therefore, we need to store the number of 218 % successful and unsuccessful moves, as well as the increase in cost 219 % for the unsuccessful moves. 220 if isempty(OPTIONS.TEMP_START) 221 [SUCCESSFUL_MOVES,UNSUCCESSFUL_MOVES,UNSUCCESSFUL_COSTS] = deal(0); 222 end 223 elseif TEMP < OPTIONS.TEMP_END 224 TEMP = 0; 225 MAXITERTEMP = OPTIONS.MAX_ITER_TEMP_LAST*NDIM; 226 else 227 MAXITERTEMP = OPTIONS.MAX_ITER_TEMP*NDIM; 228 end 229 230 % construct initial simplex 231 % ------------------------- 232 233 % 1st vertex of initial simplex 234 P(1,:) = PBEST; 235 Y(1) = CALCULATE_COST(FUN,P(1,:),LB,UB,varargin{:}); 236 237 % if output function given then run output function to plot intermediate result 238 if ~isempty(OPTIONS.OUTPUT_FCN) 239 feval(OPTIONS.OUTPUT_FCN,transpose(P(1,:)),Y(1)); 240 end 241 242 % remaining vertices of simplex 243 for k = 1:NDIM 244 % copy first vertex in new vertex 245 P(k+1,:) = P(1,:); 246 % alter new vertex 247 P(k+1,k) = LB(k)+rand*(UB(k)-LB(k)); 248 % calculate value of objective function at new vertex 249 Y(k+1) = CALCULATE_COST(FUN,P(k+1,:),LB,UB,varargin{:}); 250 end 251 252 % store information on what step the algorithm just did 253 ALGOSTEP = 'initial simplex'; 254 255 % add NDIM+1 to number of function evaluations 256 nFUN_EVALS = nFUN_EVALS + NDIM; 257 258 % note: 259 % dimensions of matrix P: (NDIM+1) x NDIM 260 % dimensions of vector Y: (NDIM+1) x 1 261 262 % give user feedback if requested 263 if strcmp(OPTIONS.DISPLAY,'iter') 264 if nITERATIONS == 0 265 disp(' Nr Iter Nr Fun Eval Min function Best function TEMP Algorithm Step'); 266 else 267 disp(sprintf('%5.0f %5.0f %12.6g %15.6g %12.6g %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,'best point')); 268 end 269 end 270 271 % run full metropolis cycle at current temperature 272 % ------------------------------------------------ 273 274 % initialize vector COSTS, needed to calculate new temperature using cooling 275 % schedule as described by Cardoso et al. (1996) 276 COSTS = zeros((NDIM+1)*MAXITERTEMP,1); 277 278 % initialize ITERTEMP to zero 279 280 ITERTEMP = 0; 281 282 % start 283 284 for ITERTEMP = 1:MAXITERTEMP 285 286 % add one to number of iterations 287 nITERATIONS = nITERATIONS + 1; 288 289 % Press and Teukolsky (1991) add a positive logarithmic distributed variable, 290 % proportional to the control temperature T to the function value associated with 291 % every vertex of the simplex. Likewise,they subtract a similar random variable 292 % from the function value at every new replacement point. 293 % Thus, if the replacement point corresponds to a lower cost, this method always 294 % accepts a true down hill step. If, on the other hand, the replacement point 295 % corresponds to a higher cost, an uphill move may be accepted, depending on the 296 % relative COSTS of the perturbed values. 297 % (taken from Cardoso et al.,1996) 298 299 % add random fluctuations to function values of current vertices 300 YFLUCT = Y+TEMP*abs(log(rand(NDIM+1,1))); 301 302 % reorder YFLUCT, Y and P so that the first row corresponds to the lowest YFLUCT value 303 help = sortrows([YFLUCT,Y,P],1); 304 YFLUCT = help(:,1); 305 Y = help(:,2); 306 P = help(:,3:end); 307 308 if nargout>3 309 % store temperature at current iteration 310 OUTPUT.TEMPERATURE(nITERATIONS) = TEMP; 311 312 % store information about simplex at the current iteration 313 OUTPUT.SIMPLEX(:,:,nITERATIONS) = P; 314 OUTPUT.SIMPLEX_BEST(nITERATIONS,:) = PBEST; 315 316 % store cost function value of best vertex in current iteration 317 OUTPUT.COSTS(nITERATIONS,:) = Y; 318 OUTPUT.COST_BEST(nITERATIONS) = YBEST; 319 end 320 321 if strcmp(OPTIONS.DISPLAY,'iter') 322 disp(sprintf('%5.0f %5.0f %12.6g %15.6g %12.6g %s',nITERATIONS,nFUN_EVALS,Y(1),YBEST,TEMP,ALGOSTEP)); 323 end 324 325 % if output function given then run output function to plot intermediate result 326 if ~isempty(OPTIONS.OUTPUT_FCN) 327 feval(OPTIONS.OUTPUT_FCN,transpose(P(1,:)),Y(1)); 328 end 329 330 % end the optimization if one of the stopping criteria is met 331 %% 1. difference between best and worst function evaluation in simplex is smaller than TOLFUN 332 %% 2. maximum difference between the coordinates of the vertices in simplex is less than TOLX 333 %% 3. no convergence,but maximum number of iterations has been reached 334 %% 4. no convergence,but maximum time has been reached 335 336 if (abs(max(Y)-min(Y)) < OPTIONS.TOLFUN) && (TEMP_LOOP_NUMBER ~= 1) 337 if strcmp(OPTIONS.DISPLAY,'iter') 338 disp('Change in the objective function value less than the specified tolerance (TOLFUN).') 339 end 340 EXITFLAG = 1; 341 break 342 end 343 344 if (max(max(abs(P(2:NDIM+1,:)-P(1:NDIM,:)))) < OPTIONS.TOLX) && (TEMP_LOOP_NUMBER ~= 1) 345 if strcmp(OPTIONS.DISPLAY,'iter') 346 disp('Change in X less than the specified tolerance (TOLX).') 347 end 348 EXITFLAG = 2; 349 break 350 end 351 352 if (nITERATIONS >= OPTIONS.MAX_ITER_TOTAL*NDIM) || (nFUN_EVALS >= OPTIONS.MAX_FUN_EVALS*NDIM*(NDIM+1)) 353 if strcmp(OPTIONS.DISPLAY,'iter') 354 disp('Maximum number of function evaluations or iterations reached.'); 355 end 356 EXITFLAG = 0; 357 break 358 end 359 360 if toc/60 > OPTIONS.MAX_TIME 361 if strcmp(OPTIONS.DISPLAY,'iter') 362 disp('Exceeded maximum time.'); 363 end 364 EXITFLAG = -1; 365 break 366 end 367 368 % begin a new iteration 369 370 %% first extrapolate by a factor -1 through the face of the simplex 371 %% across from the high point,i.e.,reflect the simplex from the high point 372 [YFTRY,YTRY,PTRY] = AMOTRY(FUN,P,-1,LB,UB,varargin{:}); 373 374 %% check the result 375 if YFTRY <= YFLUCT(1) 376 %% gives a result better than the best point,so try an additional 377 %% extrapolation by a factor 2 378 [YFTRYEXP,YTRYEXP,PTRYEXP] = AMOTRY(FUN,P,-2,LB,UB,varargin{:}); 379 if YFTRYEXP < YFTRY 380 P(end,:) = PTRYEXP; 381 Y(end) = YTRYEXP; 382 ALGOSTEP = 'reflection and expansion'; 383 else 384 P(end,:) = PTRY; 385 Y(end) = YTRY; 386 ALGOSTEP = 'reflection'; 387 end 388 elseif YFTRY >= YFLUCT(NDIM) 389 %% the reflected point is worse than the second-highest, so look 390 %% for an intermediate lower point, i.e., do a one-dimensional 391 %% contraction 392 [YFTRYCONTR,YTRYCONTR,PTRYCONTR] = AMOTRY(FUN,P,-0.5,LB,UB,varargin{:}); 393 if YFTRYCONTR < YFLUCT(end) 394 P(end,:) = PTRYCONTR; 395 Y(end) = YTRYCONTR; 396 ALGOSTEP = 'one dimensional contraction'; 397 else 398 %% can't seem to get rid of that high point, so better contract 399 %% around the lowest (best) point 400 X = ones(NDIM,NDIM)*diag(P(1,:)); 401 P(2:end,:) = 0.5*(P(2:end,:)+X); 402 for k=2:NDIM 403 Y(k) = CALCULATE_COST(FUN,P(k,:),LB,UB,varargin{:}); 404 end 405 ALGOSTEP = 'multiple contraction'; 406 end 407 else 408 %% if YTRY better than second-highest point, use this point 409 P(end,:) = PTRY; 410 Y(end) = YTRY; 411 ALGOSTEP = 'reflection'; 412 end 413 414 % the initial temperature is estimated in the first loop from 415 % the number of successfull and unsuccesfull moves, and the average 416 % increase in cost associated with the unsuccessful moves 417 418 if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START) 419 if Y(1) > Y(end) 420 SUCCESSFUL_MOVES = SUCCESSFUL_MOVES+1; 421 elseif Y(1) <= Y(end) 422 UNSUCCESSFUL_MOVES = UNSUCCESSFUL_MOVES+1; 423 UNSUCCESSFUL_COSTS = UNSUCCESSFUL_COSTS+(Y(end)-Y(1)); 424 end 425 end 426 427 end 428 429 % stop if previous for loop was broken due to some stop criterion 430 if ITERTEMP < MAXITERTEMP 431 break 432 end 433 434 % store cost function values in COSTS vector 435 COSTS((ITERTEMP-1)*NDIM+1:ITERTEMP*NDIM+1) = Y; 436 437 % calculated initial temperature or recalculate temperature 438 % using cooling schedule as proposed by Cardoso et al. (1996) 439 % ----------------------------------------------------------- 440 441 if TEMP_LOOP_NUMBER == 1 && isempty(OPTIONS.TEMP_START) 442 TEMP = -(UNSUCCESSFUL_COSTS/(SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES))/log(((SUCCESSFUL_MOVES+UNSUCCESSFUL_MOVES)*OPTIONS.INITIAL_ACCEPTANCE_RATIO-SUCCESSFUL_MOVES)/UNSUCCESSFUL_MOVES); 443 elseif TEMP_LOOP_NUMBER ~= 0 444 STDEV_Y = std(COSTS); 445 COOLING_FACTOR = 1/(1+TEMP*log(1+OPTIONS.COOL_RATE)/(3*STDEV_Y)); 446 TEMP = TEMP*min(OPTIONS.MIN_COOLING_FACTOR,COOLING_FACTOR); 447 end 448 449 % add one to temperature loop number 450 TEMP_LOOP_NUMBER = TEMP_LOOP_NUMBER+1; 451 452end 453 454% return solution 455X = transpose(PBEST); 456FVAL = YBEST; 457 458if nargout>3 459 % store number of function evaluations 460 OUTPUT.nFUN_EVALS = nFUN_EVALS; 461 462 % store number of iterations 463 OUTPUT.nITERATIONS = nITERATIONS; 464 465 % trim OUTPUT data structure 466 OUTPUT.TEMPERATURE(nITERATIONS+1:end) = []; 467 OUTPUT.SIMPLEX(:,:,nITERATIONS+1:end) = []; 468 OUTPUT.SIMPLEX_BEST(nITERATIONS+1:end,:) = []; 469 OUTPUT.COSTS(nITERATIONS+1:end,:) = []; 470 OUTPUT.COST_BEST(nITERATIONS+1:end) = []; 471 472 % store the amount of time needed in OUTPUT data structure 473 OUTPUT.TIME = toc; 474end 475 476return 477 478% ============================================================================== 479 480% AMOTRY FUNCTION 481% --------------- 482 483function [YFTRY,YTRY,PTRY] = AMOTRY(FUN,P,fac,LB,UB,varargin) 484% Extrapolates by a factor fac through the face of the simplex across from 485% the high point, tries it, and replaces the high point if the new point is 486% better. 487 488global NDIM TEMP 489 490% calculate coordinates of new vertex 491psum = sum(P(1:NDIM,:))/NDIM; 492PTRY = psum*(1-fac)+P(end,:)*fac; 493 494% evaluate the function at the trial point. 495YTRY = CALCULATE_COST(FUN,PTRY,LB,UB,varargin{:}); 496% substract random fluctuations to function values of current vertices 497YFTRY = YTRY-TEMP*abs(log(rand(1))); 498 499return 500 501% ============================================================================== 502 503% COST FUNCTION EVALUATION 504% ------------------------ 505 506function [YTRY] = CALCULATE_COST(FUN,PTRY,LB,UB,varargin) 507 508global YBEST PBEST NDIM nFUN_EVALS 509 510for i = 1:NDIM 511 % check lower bounds 512 if PTRY(i) < LB(i) 513 YTRY = 1e12+(LB(i)-PTRY(i))*1e6; 514 return 515 end 516 % check upper bounds 517 if PTRY(i) > UB(i) 518 YTRY = 1e12+(PTRY(i)-UB(i))*1e6; 519 return 520 end 521end 522 523% calculate cost associated with PTRY 524YTRY = feval(FUN,PTRY(:),varargin{:}); 525 526% add one to number of function evaluations 527nFUN_EVALS = nFUN_EVALS + 1; 528 529% save the best point ever 530if YTRY < YBEST 531 YBEST = YTRY; 532 PBEST = PTRY; 533end 534 535return 536