SECTION 16 GLOBAL OPTIMIZATION
SECTION 16.1 HYPOTHESIS TEST FOR A GLOBAL OPTMIMUM (GLOB1)
Assume that a local optimum has been found and it is desired to test the hypothesis that it is a global optimum. SUBROUTINE GLOB1 evaluates the function to be optimized n times randomly in a region specified by the user and tests the hypothesis that the previous local optimum is global. See Veall, M., "Testing for a Global Maximum in an Econometric Context," Econometrica, 58 (1990), 1459-1466. The call to GLOB1 is
CALL GLOB1(X,NP,F,XL,XH,N,P,MX,L1,L2,LP,FUNC)
where
X = the np-dimensional vector representing the point that is
to be tested for being a global optimum;
NP = the number of variables in X;
F = the function value at X;
XL = an np-dimensional vector representing the lower bounds on
the search area;
XH = an np-dimensional vector representing the upper bounds on
the search area;
N = the number of times the function will be evaluated;
P = the significance level (e.g.,0.1 or 0.05 or 0.01);
MX = 1 if function is the be maximized, =2 if minimized;
L1 = the largest function value encountered;
L2 = the second largest function value encountered;
LP = the criterion computed; to be compared with F;
FUNC = the name of the function;
IER = error code.
L1, L2, and LP have to be declared REAL*8 in the calling routine. The calling routine must also contain
EXTERNAL FUNC
GLOB1 prints results according to the values of IPT and JPT (see Section 1.1). If the user wishes to override the defaults, he should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
in the calling routine and reset the values of IPT and JPT. The SUBROUTINE FUNC should obey the same conventions as described in Section 1. GLOB1 automatically searches for candidate points over the rectangular region bounded by the vectors XL and XH. If the user wishes to search over some other region, he should write a SUBROUTINE YGEN(Y,NP) which generates np-dimensional points Y over the desired region. The user's subroutine will be loaded over the subroutine present in GQOPT.
Error codes:
IER = -80 Function cannot be evaluated at point passed
IER = -81 Function error at first search point
IER = -82 Function error at second search point
IER = -83 Function error at subsequent search points
SECTION 16.2 GLOBAL OPTIMIZATION BY SIMULATED ANNEALING
This routine implements the simulated annealing global optimization algorithm described in Corana, A., M. Marchesi, C. Martini, and S. Ridella, "Minimizing Multimodal Functions of Continuous Variables with the 'Simulated Annealing' Algorithm," ACM Transactions on Mathematical Software, 13 (1987), 262-280. The authors of this implementation describe their modifications of the original Corana, etal., paper and their own experience with it in W. L. Goffe, G. Ferrier, and J. Rogers, "Global Optimization of Statistical Functions," Souther Methodist University, May 1990. Please communicate comments to William L. Goffe, Department of Economics, Southern Methodist University, Dallas, TX 76275; e-mail to H2ZR1001@SMUVM1 ; in addition to Richard E. Quandt at Princeton. This algorithm is incorporated in GQOPT with permisssion of its authors. To call the algorithm, issue the command
CALL GLOB2(N,X,MAX,RT,EPS,NS,NT,NEPS,MAXEVL,LB,UB,C,
1 T,VM,XOPT,FOPT,NACC,NFCNEV,NOBDS,IER, 2 FSTAR,XP,NACP,FUNC)
In describing the various input and output parameters, the following conventions will be employed: DP denotes double precision (REAL*8) and SP denotes single precision (REAL*4). INT denotes integer. DP(N) denotes a double precision array with N eleme
nts. GLOB2 prints results according to the values of IPT and JPT (see Section 1.1). If the user wishes to override the defaults, he should include
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
in the calling routine and reset the values of IPT and JPT.
Other Common Blocks:
COMMON/BGLOB2/IERCTL
where IERCTL gives the number of FUNC errors permitted before termination (Default=10)
The SUBROUTINE FUNC should obey the same conventions as described in Section 1.
Input Parameters: NOTE The suggested values generally come from Corana, etal. To drastically reduce runtime, see Goffe, etal., pp. 17-18 for suggestions for choosing the appropriate RT and NT.
N = the number of variables in the function to be optimized
(INT)
X = starting values for the variables in the function to be
optimized (DP(N))
MAX = 1 if maximizing, = 2 if minimizing (INT)
RT = the temperature reduction factor (e.g. 0.85) (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 currnt temperature differs
from the current optimal function value by less than EPS,
execution terminates and IER=0 is returned. (DP)
NS = Number of cycles. After NS*N function evaluations, 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 is reduced by the
factor RT. Value suggested by Corona, etal., is MAX(100,5*N).
(INT)
NEPS = Number of final function values used to decide upon
termination. See EPS. Suggested value is 4. (INT)
MAXEVL = The
maximum number of function evaluations. Note that out of bounds
trials do not count, since the function is not actually evaluated
(see LB and UB). If it is exceeded, IER is set to -1(INT).
LB = The lower bound for the allowable solution values.(DP(N))
UB = The upper bound for the allowable solution values.(DP(N)) If
the algorithm chooses X(I).LT.LB(I) or X(I).GT.UB(I), I=1,N, the
point is rejected and the function is not evaluated. This focuses
the algorithm on the region inside UB and LB. Unless the user
wishes to concentrate the search to a particular region, UB and
LB should be set to very large positive and negative values
respectively. Also note that UB and LB are fixed in position,
wheras VM is centered on the last accepted trial set of variables
that optimizes the function.
C = Vector that controls the step length adjustment. Suggested
value is 2.0. (DP(N))
IPT = IPT and JPT print on logical units NFILE and MFILE
respectively. JPT Default values are 2. IPT and JPT can be reset
by including in the MAIN program the common block
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
IPT, JPT = 0 nothing printed = 1 function value for the
starting value and summary results before temperature reduction;
notice given upon chieving termination criteria. = 2 each new
step length (VM) = 3 each function evaluation, its acceptance or
rejection, and new optima (for many problems, this option will
produce massive output). This option is most useful for learning
about the algorithm when used with a small value of MAXEVL. See
also Section 1.1.
IERCTL = number of FUNC errors permitted before termination
(Default=10)
FUNC = the name of the SUBROUTINE that computes the function to
be maximized or minimized. The name of this subroutine must be
declared in an EXTERNAL statement in the MAIN program. For the
conventions governing the function SUBROUTINE, see GENERAL
INTRODUCTION in SECTION 1. INPUT/OUTPUT PARAMETERS
T = On input, the initial temperature. On output, the final
temperature. (DP)
VM = The step length vector. On input, it should encompass the
region of interest given the starting value of X. For point X(I),
the next trial point is selected from X(I)-VM(I) to X(I)+VM(I).
Since VM is adjusted so that about half of all points is
accepted, the input value is not very important (i.e., if the
value is off, GLOB2 adjusts to the correct value). (DP(N))
Output Parameters:
XOPT = The variable values that optimize the function. (DP(N))
FOPT = The optimal value of the function. (DP)
NACC = The number of accepted steps. (INT)
NFCNEV = The total number of function evaluations. (INT)
NOBDS = The total number of trial function evaluations out of
bounds relative to LB and UB. The function is not evaluated and
this is not counted in NFCNEV. (INT)
IER = Error return. (INT) = 0 Normal return, termination criteria
achieved.
IER = -1 Number of function evaluations (NFCNEV) is greater than
the maximum number (MAXEVL).
IER = -84 Function error on first FUNC call.
IER =-85 Persistent function errors -- may wish to reset IERCTL.
IER = -86 The starting X-vector is out of bounds (X(I) must be
between LB(I) and UB(I) for all I=1,N).
IER = -99 Should not be seen; only used internally.
Work Arrays that Have to be Dimensioned in the Calling Program:
FSTAR (DP(NEPS))
XP (DP(N))
NACP (INT(N))
Functions used:
EXPREP = Replaces the function EXP to avoid under- and
overflows. It may have to be modified for non-IBM mainframes.
(DP)
RNDY5 = Random number generator in GQOPT (DP).
Subroutines used:
PRTVCA, PRTVCB, PRT0A,...,PRT8A, PRTB0,...,PRT8N = printing
routines
FUNC = user specified
Machine Specific Features:
1. EXPREP may have to be modified for non-IBM type machines
2. Some FORMAT statements use G25.18; this may be excessive for
some machines.
General Hints:
The parameter T is crucial in using GLOB2 successfully. It influences VM, the step length over which the algorithm searches for optima. For a small starting T, it is possible for the step length to be too small; thus not enough of the function might be evaluated for finding the global maximum or minimum. The relationship between the initial temperature and the resulting step length is function dependent. To determine the starting temperature that is consistent with sampling a lot of the function, it is worth while to run a trial run first. Set RT=1.5 and T=1.0. With RT .GT. 1.0, the temperature increases and VM rises as well. Then select that T that produces a large enough VM as the initial T for an optimizing run. Basically, the way this works is as follows. The algorithm samples within VM. In a maximization problem, all uphill moves are accepted, while downhill moves use the Metropolis criterion, which depends on T. As T falls, a given downhill move is less likely to be accepted. Each element of VM is adjusted so that half of all function evaluations are accepted. Thus, as T declines, downhill moves become less likely to be accepted and the percentage of rejections rises. As a result, VM falls. One hopes that this occurs near the global optimum.