Comments extracted from simann.for

SIMANN COMMENTS

The simann.for comments are a lesson in themselves, including even citations. They are extracted from the body of the program below. The papers cited are worth reading. You should print this out and look at it while going over the program code. --D. Baker 7/11/2007

C ABSTRACT:
C Simulated annealing is a global optimization method that distinguishes
C between different local optima. Starting from an initial point, the
C algorithm takes a step and the function is evaluated. When minimizing a
C function, any downhill step is accepted and the process repeats from this
C new point. An uphill step may be accepted. Thus, it can escape from local
C optima. This uphill decision is made by the Metropolis criteria. As the
C optimization process proceeds, the length of the steps decline and the
C algorithm closes in on the global optimum. Since the algorithm makes very
C few assumptions regarding the function to be optimized, it is quite
C robust with respect to non-quadratic surfaces. The degree of robustness
C can be adjusted by the user. In fact, simulated annealing can be used as
C a local optimizer for difficult functions.
C
C This implementation of simulated annealing was used in "Global Optimization
C of Statistical Functions with Simulated Annealing," Goffe, Ferrier and
C Rogers, Journal of Econometrics, vol. 60, no. 1/2, Jan./Feb. 1994, pp.
C 65-100. Briefly, we found it competitive, if not superior, to multiple
C restarts of conventional optimization routines for difficult optimization
C problems.
C
C For more information on this routine, contact its author:
C Bill Goffe, bgoffe@whale.st.usm.edu
C
PROGRAM SIMANN
C This file is an example of the Corana et al. simulated annealing
C algorithm for multimodal and robust optimization as implemented
C and modified by Goffe, Ferrier and Rogers. Counting the above line
C ABSTRACT as 1, the routine itself (SA), with its supplementary
C routines, is on lines 232-990. A multimodal example from Judge et al.
C (FCN) is on lines 150-231. The rest of this file (lines 1-149) is a
C driver routine with values appropriate for the Judge example. Thus, this
C example is ready to run.
C
C To understand the algorithm, the documentation for SA on lines 236-
C 484 should be read along with the parts of the paper that describe
C simulated annealing. Then the following lines will then aid the user
C in becomming proficient with this implementation of simulated
C annealing.
C
C Learning to use SA:
C Use the sample function from Judge with the following suggestions
C to get a feel for how SA works. When you've done this, you should be
C ready to use it on most any function with a fair amount of expertise.
C 1. Run the program as is to make sure it runs okay. Take a look at
C the intermediate output and see how it optimizes as temperature
C (T) falls. Notice how the optimal point is reached and how
C falling T reduces VM.
C 2. Look through the documentation to SA so the following makes a
C bit of sense. In line with the paper, it shouldn't be that hard
C to figure out. The core of the algorithm is described on pp. 68-70
C and on pp. 94-95. Also see Corana et al. pp. 264-9.
C 3. To see how it selects points and makes decisions about uphill
C and downhill moves, set IPRINT = 3 (very detailed intermediate
C output) and MAXEVL = 100 (only 100 function evaluations to limit
C output).
C 4. To see the importance of different temperatures, try starting
C with a very low one (say T = 10E-5). You'll see (i) it never
C escapes from the local optima (in annealing terminology, it
C quenches) & (ii) the step length (VM) will be quite small. This
C is a key part of the algorithm: as temperature (T) falls, step
C length does too. In a minor point here, note how VM is quickly
C reset from its initial value. Thus, the input VM is not very
C important. This is all the more reason to examine VM once the
C algorithm is underway.
C 5. To see the effect of different parameters and their effect on
C the speed of the algorithm, try RT = .95 & RT = .1. Notice the
C vastly different speed for optimization. Also try NT = 20. Note
C that this sample function is quite easy to optimize, so it will
C tolerate big changes in these parameters. RT and NT are the
C parameters one should adjust to modify the runtime of the
C algorithm and its robustness.
C 6. Try constraining the algorithm with either LB or UB.

SUBROUTINE SA(N,X,MAX,RT,EPS,NS,NT,NEPS,MAXEVL,LB,UB,C,IPRINT,
1 ISEED1,ISEED2,T,VM,XOPT,FOPT,NACC,NFCNEV,NOBDS,IER,
2 FSTAR,XP,NACP)

C Version: 3.2
C Date: 1/22/94.
C Differences compared to Version 2.0:
C 1. If a trial is out of bounds, a point is randomly selected
C from LB(i) to UB(i). Unlike in version 2.0, this trial is
C evaluated and is counted in acceptances and rejections.
C All corresponding documentation was changed as well.
C Differences compared to Version 3.0:
C 1. If VM(i) > (UB(i) - LB(i)), VM is set to UB(i) - LB(i).
C The idea is that if T is high relative to LB & UB, most
C points will be accepted, causing VM to rise. But, in this
C situation, VM has little meaning; particularly if VM is
C larger than the acceptable region. Setting VM to this size
C still allows all parts of the allowable region to be selected.
C Differences compared to Version 3.1:
C 1. Test made to see if the initial temperature is positive.
C 2. WRITE statements prettied up.
C 3. References to paper updated.
C
C Synopsis:
C This routine implements the continuous simulated annealing global
C optimization algorithm described in Corana et al.'s article
C "Minimizing Multimodal Functions of Continuous Variables with the
C "Simulated Annealing" Algorithm" in the September 1987 (vol. 13,
C no. 3, pp. 262-280) issue of the ACM Transactions on Mathematical
C Software.
C
C A very quick (perhaps too quick) overview of SA:
C SA tries to find the global optimum of an N dimensional function.
C It moves both up and downhill and as the optimization process
C proceeds, it focuses on the most promising area.
C To start, it randomly chooses a trial point within the step length
C VM (a vector of length N) of the user selected starting point. The
C function is evaluated at this trial point and its value is compared
C to its value at the initial point.
C In a maximization problem, all uphill moves are accepted and the
C algorithm continues from that trial point. Downhill moves may be
C accepted; the decision is made by the Metropolis criteria. It uses T
C (temperature) and the size of the downhill move in a probabilistic
C manner. The smaller T and the size of the downhill move are, the more
C likely that move will be accepted. If the trial is accepted, the
C algorithm moves on from that point. If it is rejected, another point
C is chosen instead for a trial evaluation.
C Each element of VM periodically adjusted so that half of all
C function evaluations in that direction are accepted.
C A fall in T is imposed upon the system with the RT variable by
C T(i+1) = RT*T(i) where i is the ith iteration. Thus, as T declines,
C downhill moves are less likely to be accepted and the percentage of
C rejections rise. Given the scheme for the selection for VM, VM falls.
C Thus, as T declines, VM falls and SA focuses upon the most promising
C area for optimization.
C
C The importance of the parameter T:
C The parameter T is crucial in using SA successfully. It influences
C VM, the step length over which the algorithm searches for optima. For
C a small intial T, the step length may be too small; thus not enough
C of the function might be evaluated to find the global optima. The user
C should carefully examine VM in the intermediate output (set IPRINT =
C 1) to make sure that VM is appropriate. The relationship between the
C initial temperature and the resulting step length is function
C dependent.
C To determine the starting temperature that is consistent with
C optimizing a function, it is worthwhile to run a trial run first. Set
C RT = 1.5 and T = 1.0. With RT > 1.0, the temperature increases and VM
C rises as well. Then select the T that produces a large enough VM.
C
C For modifications to the algorithm and many details on its use,
C (particularly for econometric applications) see Goffe, Ferrier
C and Rogers, "Global Optimization of Statistical Functions with
C Simulated Annealing," Journal of Econometrics, vol. 60, no. 1/2,
C Jan./Feb. 1994, pp. 65-100.
C For more information, contact
C Bill Goffe
C Department of Economics and International Business
C University of Southern Mississippi
C Hattiesburg, MS 39506-5072
C (601) 266-4484 (office)
C (601) 266-4920 (fax)
C bgoffe@whale.st.usm.edu (Internet)
C
C As far as possible, the parameters here have the same name as in
C the description of the algorithm on pp. 266-8 of Corana et al.
C
C In this description, SP is single precision, DP is double precision,
C INT is integer, L is logical and (N) denotes an array of length n.
C Thus, DP(N) denotes a double precision array of length n.
C
C Input Parameters:
C Note: The suggested values generally come from Corana et al. To
C drastically reduce runtime, see Goffe et al., pp. 90-1 for
C suggestions on choosing the appropriate RT and NT.
C N - Number of variables in the function to be optimized. (INT)
C X - The starting values for the variables of the function to be
C optimized. (DP(N))
C MAX - Denotes whether the function should be maximized or
C minimized. A true value denotes maximization while a false
C value denotes minimization. Intermediate output (see IPRINT)
C takes this into account. (L)
C RT - The temperature reduction factor. The value suggested by
C Corana et al. is .85. See Goffe et al. for more advice. (DP)
C EPS - Error tolerance for termination. If the final function
C values from the last neps temperatures differ from the
C corresponding value at the current temperature by less than
C EPS and the final function value at the current temperature
C differs from the current optimal function value by less than
C EPS, execution terminates and IER = 0 is returned. (EP)
C NS - Number of cycles. After NS*N function evaluations, each
C element of VM is adjusted so that approximately half of
C all function evaluations are accepted. The suggested value
C is 20. (INT)
C NT - Number of iterations before temperature reduction. After
C NT*NS*N function evaluations, temperature (T) is changed
C by the factor RT. Value suggested by Corana et al. is
C MAX(100, 5*N). See Goffe et al. for further advice. (INT)
C NEPS - Number of final function values used to decide upon termi-
C nation. See EPS. Suggested value is 4. (INT)
C MAXEVL - The maximum number of function evaluations. If it is
C exceeded, IER = 1. (INT)
C LB - The lower bound for the allowable solution variables. (DP(N))
C UB - The upper bound for the allowable solution variables. (DP(N))
C If the algorithm chooses X(I) .LT. LB(I) or X(I) .GT. UB(I),
C I = 1, N, a point is from inside is randomly selected. This
C This focuses the algorithm on the region inside UB and LB.
C Unless the user wishes to concentrate the search to a par-
C ticular region, UB and LB should be set to very large positive
C and negative values, respectively. Note that the starting
C vector X should be inside this region. Also note that LB and
C UB are fixed in position, while VM is centered on the last
C accepted trial set of variables that optimizes the function.
C C - Vector that controls the step length adjustment. The suggested
C value for all elements is 2.0. (DP(N))
C IPRINT - controls printing inside SA. (INT)
C Values: 0 - Nothing printed.
C 1 - Function value for the starting value and
C summary results before each temperature
C reduction. This includes the optimal
C function value found so far, the total
C number of moves (broken up into uphill,
C downhill, accepted and rejected), the
C number of out of bounds trials, the
C number of new optima found at this
C temperature, the current optimal X and
C the step length VM. Note that there are
C N*NS*NT function evalutations before each
C temperature reduction. Finally, notice is
C is also given upon achieveing the termination
C criteria.
C 2 - Each new step length (VM), the current optimal
C X (XOPT) and the current trial X (X). This
C gives the user some idea about how far X
C strays from XOPT as well as how VM is adapting
C to the function.
C 3 - Each function evaluation, its acceptance or
C rejection and new optima. For many problems,
C this option will likely require a small tree
C if hard copy is used. This option is best
C used to learn about the algorithm. A small
C value for MAXEVL is thus recommended when
C using IPRINT = 3.
C Suggested value: 1
C Note: For a given value of IPRINT, the lower valued
C options (other than 0) are utilized.
C ISEED1 - The first seed for the random number generator RANMAR.
C 0 .LE. ISEED1 .LE. 31328. (INT)
C ISEED2 - The second seed for the random number generator RANMAR.
C 0 .LE. ISEED2 .LE. 30081. Different values for ISEED1
C and ISEED2 will lead to an entirely different sequence
C of trial points and decisions on downhill moves (when
C maximizing). See Goffe et al. on how this can be used
C to test the results of SA. (INT)
C
C Input/Output Parameters:
C T - On input, the initial temperature. See Goffe et al. for advice.
C On output, the final temperature. (DP)
C VM - The step length vector. On input it should encompass the
C region of interest given the starting value X. For point
C X(I), the next trial point is selected is from X(I) - VM(I)
C to X(I) + VM(I). Since VM is adjusted so that about half
C of all points are accepted, the input value is not very
C important (i.e. is the value is off, SA adjusts VM to the
C correct value). (DP(N))
C
C Output Parameters:
C XOPT - The variables that optimize the function. (DP(N))
C FOPT - The optimal value of the function. (DP)
C NACC - The number of accepted function evaluations. (INT)
C NFCNEV - The total number of function evaluations. In a minor
C point, note that the first evaluation is not used in the
C core of the algorithm; it simply initializes the
C algorithm. (INT).
C NOBDS - The total number of trial function evaluations that
C would have been out of bounds of LB and UB. Note that
C a trial point is randomly selected between LB and UB.
C (INT)
C IER - The error return number. (INT)
C Values: 0 - Normal return; termination criteria achieved.
C 1 - Number of function evaluations (NFCNEV) is
C greater than the maximum number (MAXEVL).
C 2 - The starting value (X) is not inside the
C bounds (LB and UB).
C 3 - The initial temperature is not positive.
C 99 - Should not be seen; only used internally.
C
C Work arrays that must be dimensioned in the calling routine:
C RWK1 (DP(NEPS)) (FSTAR in SA)
C RWK2 (DP(N)) (XP " " )
C IWK (INT(N)) (NACP " " )
C
C Required Functions (included):
C EXPREP - Replaces the function EXP to avoid under- and overflows.
C It may have to be modified for non IBM-type main-
C frames. (DP)
C RMARIN - Initializes the random number generator RANMAR.
C RANMAR - The actual random number generator. Note that
C RMARIN must run first (SA does this). It produces uniform
C random numbers on [0,1]. These routines are from
C Usenet's comp.lang.fortran. For a reference, see
C "Toward a Universal Random Number Generator"
C by George Marsaglia and Arif Zaman, Florida State
C University Report: FSU-SCRI-87-50 (1987).
C It was later modified by F. James and published in
C "A Review of Pseudo-random Number Generators." For
C further information, contact stuart@ads.com. These
C routines are designed to be portable on any machine
C with a 24-bit or more mantissa. I have found it produces
C identical results on a IBM 3081 and a Cray Y-MP.
C
C Required Subroutines (included):
C PRTVEC - Prints vectors.
C PRT1 ... PRT10 - Prints intermediate output.
C FCN - Function to be optimized. The form is
C SUBROUTINE FCN(N,X,F)
C INTEGER N
C DOUBLE PRECISION X(N), F
C ...
C function code with F = F(X)
C ...
C RETURN
C END
C Note: This is the same form used in the multivariable
C minimization algorithms in the IMSL edition 10 library.
C
C Machine Specific Features:
C 1. EXPREP may have to be modified if used on non-IBM type main-
C frames. Watch for under- and overflows in EXPREP.
C 2. Some FORMAT statements use G25.18; this may be excessive for
C some machines.
C 3. RMARIN and RANMAR are designed to be protable; they should not
C cause any problems.

subroutine RMARIN(IJ,KL)
C This subroutine and the next function generate random numbers. See
C the comments for SA for more information. The only changes from the
C orginal code is that (1) the test to make sure that RMARIN runs first
C was taken out since SA assures that this is done (this test didn't
C compile under IBM's VS Fortran) and (2) typing ivec as integer was
C taken out since ivec isn't used. With these exceptions, all following
C lines are original.

C This is the initialization routine for the random number generator
C RANMAR()
C NOTE: The seed variables can have values between: 0 <= IJ <= 31328
C 0 <= KL <= 30081
SUBROUTINE PRT1
C This subroutine prints intermediate output, as does PRT2 through
C PRT10. Note that if SA is minimizing the function, the sign of the
C function value and the directions (up/down) are reversed in all
C output to correspond with the actual function optimization. This
C correction is because SA was written to maximize functions and
C it minimizes by maximizing the negative a function.

SUBROUTINE PRTVEC(VECTOR,NCOLS,NAME)
C This subroutine prints the double precision vector named VECTOR.
C Elements 1 thru NCOLS will be printed. NAME is a character variable
C that describes VECTOR. Note that if NAME is given in the call to
C PRTVEC, it must be enclosed in quotes. If there are more than 10
C elements in VECTOR, 10 elements will be printed on each line.