SHORT DESCRIPTION OF THE GLOBAL OPTIMIZATION SUBROUTINE
-------------------------------------------------------
Global optimization is a part of nonlinear optimization, it deals with
problems with (possibly) several local minima. The presented method is
stochastic (i.e. not deterministic). The framework procedure, the GLOBAL
routine gives a computational evidence, that the best local minimum found is
with high probability the global minimum. This routine calls a local search
routine, and a routine for generating random numbers.
The subroutines are written in FORTRAN. The user has to provide a main program
doing the input - output, and a subroutine capable of computing the objective
function value for any point in the parameter space. There were different
versions for IBM-compatible mainframes and personal computers. The subroutine
GLOBAL can be used for the approximate solution of the global minimization
problem, as follows:
Let F(X) be a real function of NPARM parameters and we are looking for
parameter values X(I) from the given intervals [MIN(I), MAX(I)] for each
I = 1, 2, ..., NPARM. The problem is to determine such a point X*, that the
function value F(X) is greater than or equal to F(X*) for every X in the
NPARM-dimensional interval specified by MIN(I)'s and MAX(I)'s.
The given version allows 15 parameters to be optimized. For modifying the
program to accept larger problems change the numbers 15 everywhere in the
source lines to the new value N, and change 120 in a declaration to N(N+1)/2.
I. THE ALGORITHM
----------------
The algorithm consists of the following steps:
0. step: initialization of the parameters of the algorithm. X0 is the set of
the local minima found, F0 contains corresponding values of the
objective function. The local minima are ordered increasingly
according to the function values. X1 is the set of points starting
from which the local search procedure led to a local minimum. These
points are called seed points, and as such they are used as seed
points in the clustering phase. At the start the mentioned sets are
empty. The subroutine checks the parameter bounds, and gives error
message, and stops if they are not correct, or if they do not meet
the requirements.
1. step: generate sample points with uniform distribution, and add them to the
sample. The number of generated points is specified by NSAMPL.
2. step: generate the actual sample, as the best 100 * NSEL/NSAMPL percentage
of the sample points generated so far (according to the objective
function values). This actual sample contains in general NSEL more
points than the previous one.
3. step: form clusters in the actual sample by Single Linkage method, by
growing the clusters first around the seed points (elements of the
sets X0 and X1). A new point will join a cluster, if there is a point
in the cluster with which the distance is less than a critical
distance calculated automatically by the program for the given
problem, and if the point in the cluster has a smaller objective
function value than that of the considered one. The critical distance
depends on the number of points in the whole sample, and on the
dimension of the problem NPARM. If all points of the actual sample
are successfully ordered to some of the existing clusters, then comes
the 5th step.
4. step: start local search from the actual sample points not yet clustered in
ascending order by the values of the respective function values. If
the result of the local search is close to an element of the sets X0
and X1, then the starting point will be added to the set X1. If the
result of the local search cannot be clustered to any of the existing
clusters then the point is regarded as a new local minimizer, and is
added to X0. Choose now this result point of the local search to be a
seed point, and try to find (not clustered) points in the current
sample that can be clustered to this one. If a new local minimum was
found in step 4, go to the step 1.
5. step: determine the element of the set X0 with the smallest function value.
This is the candidate of the program for the global minimizer. Stop.
The presented program is a modification of the algorithm by Boender et al.
(see [1]). The following changes were made.
1. The local search procedure is an algorithm of Quasi-Newton type which uses
the so called DFP (Davidon-Fletcher-Powell) update formula. The comparison
of this local search method with others can be found in [2]. For smooth
objective functions this seems to be a good choice, for problems with
discontinuous objective function or derivatives the robust random search
method UNIRANDI (more details in [5]) can be recommended.
2. The subroutine GLOBAL will automatically scale the parameters to be
optimized (as in [3]) in order to keep all the starting points of the
parameters between -1 and 1 for the local minimizer procedure. The user
does not have to care about the scaling, because the results are
transformed back to the original interval before each output, and before
giving the control back to the calling program.
3. We have also made use of later results [6] on this method. For example, the
condition used at clustering has been changed to another one.
4. Instead of the Euclidean distance the greatest difference in absolute
values is used in our GLOBAL routine. This change leads to a simpler and
quicker version of the algorithm.
II. HOW TO CALL THE SUBROUTINE GLOBAL
-------------------------------------
CALL GLOBAL (MIN, MAX, NPARM, M, NSAMPL, NSEL, IPR, NSIG, X0, NC, F0)
MIN and MAX are real vectors (input) of NPARM elements. They specify an
NPARM-dimensional interval that contains the starting points for the local
searches, i.e. the Ith coordinate of a random sample point, X(I) is always in
[MIN(I), MAX(I)]. The values of this vector will not be changed during the run
of the GLOBAL subroutine.
NPARM is an integer constant (input) containing the number of the parameters
to be optimized. The value will not be changed during the run of GLOBAL.
M is an integer constant (input) containing the number of the residual
functions. The n+1 dimensional function f(i,x) is a residual function, if the
objective function can be written in the form F(x) = f(1,x)**2 + f(2,x)**2 +
... + f(M,x)**2. The parameter M was used only earlier, when a Levenberg-
Marquardt method was applied as a local search procedure. With the present
routines this parameter is not used any more, yet it can be applied to pass
values from the main program to the routine FUNCT calculating the objective
function.
NSAMPL is an integer constant (input) containing the number of the sampling
points to be generated in the parameter space with uniform distribution in one
loop. The value will not be changed during the run of the GLOBAL subroutine.
NSEL is an integer constant (input). In one sampling the number of points
selected for starting points of local search routine is NSEL. These points are
those with the smallest objective function value. The value will not be
changed by the GLOBAL routine.
IPR is an integer constant (input), the FORTRAN logical file number for the
output file. The value will not be changed during the GLOBAL run. The user
should take care on the corresponding OPEN statement in the main program. In
the IBM-PC version of the GLOBAL routine each output of the routine will be
repeated for the default unit (the screen).
NSIG is an integer constant (input), a parameter for the stopping criterion of
the local search algorithm. If the value of the objective function is the same
in the last two steps in the first NSIG significant digits, then the local
search will be cancelled. The value of NSIG will not be changed during the run
of GLOBAL.
X0 is a two-dimensional array of 15 * 20 real elements (output). After the run
of GLOBAL this matrix contains the parameters of the local minimizer points
found: the J-th local minimizer vector is stored in X0(I,J) I=1,2,...,NPARM.
These minimizers are ordered in ascending order according to their objective
function values. X0 will be set by GLOBAL.
NC is an integer constant (output) containing the number of the local minima
found. This constant will be set by the GLOBAL routine.
F0 is a real array of 20 elements (output) that contains the objective
function values of the respective local minimizer points. Thus F(X0(.,J)) =
F0(J). F0 will be set by GLOBAL.
III. THE FOLLOWING SUBROUTINES ARE CALLED BY THE SUBROUTINE GLOBAL:
-------------------------------------------------------------------
URDMN, FUN, LOCAL, UPDATE, FUNCT, TIMER
The user has to provide only the routine FUNCT (see later) out of the
mentioned. URDMN generates the uniformly distributed pseudorandom numbers. If
necessary, it can be replaced with any other routine (in case the latter
fulfills the calling requirements).
The starting value for the random number generation is given by the TIMER
routine. This starting value is based on the computer clock, hence it is
expected to be different for all runs. The routine TIMER is written in the
ASSEMBLER language of the IBM-PC (MASM). It can be changed to another routine
that asks the user what should be the starting number - if it is preferred.
FUN transforms the scaled variables back before calling the FUNCT routine. It
is provided together with the GLOBAL routine. Thus, the user can write the
FUNCT routine according to the usual parameter space.
The subroutines LOCAL and UPDATE contain the Quasi-Newton local search method.
If a random search algorithm is preferred, the routine UNIRANDI should be
linked, it has the same calling form as LOCAL.
IV. THE SUBROUTINE FUNCT FOR COMPUTING THE OBJECTIVE FUNCTION
-------------------------------------------------------------
The calling sequence is:
CALL FUNCT (X, F, NPARM, M)
where the parameters should be understood as follows:
X is a real vector of NPARM elements (input) containing the values of the
parameters, i.e. the co-ordinates of the point in the parameter space. The
variables X(I) are now not scaled. (They are transformed back to the interval
given by the user.) The vector X must not be changed in the routine FUNCT.
F is a real constant (output). This variable has to contain the computed
objective function value after returning from the routine FUNCT.
M is an integer constant (input). In earlier versions of the program it
contained the number of residual functions. This variable is not used any more
in the recent version of the GLOBAL routine, thus it can be applied to
transfer values from the main program to the FUNCT routine. It can be changed
freely.
If there are more necessary information besides the values of X, NPARM and M
for the routine FUNCT, then they can be passed through COMMON variables.
V. LIMITATIONS
--------------
parameter limitation otherwise
------------------------------------------------------------------------------
MIN(I) MIN(I) = MAX(I) error message, interrupt
MAX(I) MIN(I) = MAX(I) error message, interrupt
NPARM 1 <= NPARM <= 15 error message, interrupt
M M <= 100 error message, interrupt
NSAMPL 20 <= NSAMPL <= 10000 the value of NSAMPL is changed to
20 or 10000 whichever is closer
to the previous value
NSEL 1 <= NSEL <= 20 the value of NSEL is changed to 1
or 20 whichever is closer to the
previous value
------------------------------------------------------------------------------
The memory required for the GLOBAL subroutine and for the routines called by
it is around 80 Kbyte in the present version. The use of a numerical
coprocessor reduces the computation time substantially (it depends on the
objective function very much). The recent version is written in Microsoft
FORTRAN Version 4.01, and it is compatible with the FORTRAN 77 standard.
VI. INPUT - OUTPUT
------------------
There is no separate input of the GLOBAL subroutine, and its output goes to
the given logical device (specified by IPR). This output gives a document of
the run, and provides a list of the local minima found.
VII. THE SUGGESTED VALUES FOR THE PARAMETERS OF THE ALGORITHM
-------------------------------------------------------------
The filling out of the arrays MIN and MAX does not mean generally any problem.
The program accepts the input parameters even if MAX(I) is less than MIN(I).
However, if MIN(I) = MAX(I) the program halts with an error message. According
to this description of the algorithm, the sampling points are generated in the
NPARM-dimensional interval given by vectors MIN and MAX. However, the local
search may leave this interval. If the given limits are ment as strict limits,
the objective function should be modified to incorporate some form of penalty
for points outside the interval.
The usual value for IPR is 6. In the main program the file with the logical
file number IPR has to be opened.
NSAMPL and NSEL can affect the reliability of the GLOBAL subroutine. If NSAMPL
= NSEL, then the subroutine corresponds to the so-called Multiply Starts
method (usually it is not sufficiently efficient). It is worth to choose
NSAMPL to be at least as large as the number of function evaluations used by a
single local search. The smaller the NSEL/NSAMPL ratio, the smaller can be the
region of attraction of the global minimum. Thus, decreasing this ratio will
cause the GLOBAL routine to be more reliable.
It is not worth to give small value for NSIG. Take into account that the local
search procedure is capable to determine the location of the local minimum
only to that extent what is allowed by the objective function. As a minimum, 6
is suggested for a value of NSIG. The clustering phase of the GLOBAL routine
is more effective if the local minima are well approximated by the local
search procedure.
VIII. SAMPLE PROGRAM TO SHOW THE USAGE OF THE GLOBAL SUBROUTINE ON PC-S
-----------------------------------------------------------------------
REAL X0(15,20), F0(20), MIN(15), MAX(15)
M = 1
NPARM = 1
NSAMPL = 100
NSEL = 2
IPR = 6
OPEN(6, FILE='OUTPUT')
NSIG = 6
MIN(1) = -100.0
MAX(1) = 100.0
CALL GLOBAL(MIN, MAX, NPARM, M, NSAMPL, NSEL, IPR, NSIG, X0, NC, F0)
END
SUBROUTINE FUNCT(X, VALUE, NPARM, M)
DIMENSION X(1)
VALUE = 1.0 - COS(X(1)) + (X(1)/100.0)**2
RETURN
END
The given test example is a one-dimensional problem, it has several local
minima. The global minimum is 0, and the objective function reaches this value
at the point x(1) = 0.0. The test program found 5 local minima in 48 sec. -
among them the global one - with 523 function evaluations. Note, that if the
NSEL/NSAMPL ratio is too small, then the local minima with relatively large
objective function values cannot be found by the program. This feature is
useful if we are interested mainly only in the global minimum.
In order to run this program, one should, of course, link it with the
subroutine GLOBAL and with those called by it.
If you encounter any problems using GLOBAL, please inform
Dr. Tibor Csendes,
Institute of Informatics, Jozsef Attila University
H-6701 Szeged, Pf. 652, Hungary
Phone: +36 62 310 011 (ext. 3839)
Fax: +36 62 312 292
E-mail: csendes@inf.u-szeged.hu
URL: http://www.inf.u-szeged.hu/~csendes/
You can find additional important information in the comments of the source codes.
REFERENCES
----------
1. Boender, C.G.E., A.H.G. Rinnooy Kan, G.T. Timmer, L. Stougie: A stochastic
method for global optimization, Mathematical Programming 22(1982) 125-140.
2. Gill, P.E., W. Murray, M.H. Wright: Practical Optimization (Academic Press,
London, 1981).
3. Csendes, T., B. Daroczy, Z. Hantos: Nonlinear parameter estimation by
global optimization: comparison of local search methods in respiratory
system modelling, in: System Modelling and Optimization (Springer-Verlag,
Berlin, 1986) 188-192.
4. Csendes, T.: Nonlinear parameter estimation by global optimization -
Efficiency and reliability, Acta Cybernetica 8(1988) 361-370.
5. Jarvi, T.: A random search optimizer with an application to a maxmin
problem, Publications of the Inst. of Appl. Math., Univ. of Turku, No. 3,
1973.
6. Timmer, G.T.: Global optimization: a stochastic approach, (Ph.D. Thesis,
Erasmus University, Rotterdam, 1984).
Last update: 8.6.1995