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