COMPARISONS AMONG STOCHASTIC
OPTIMIZATION ALGORITIIMS
By
DEBAO CHEN
Bachelor of Science
Shanghai Institute ofEducation
Shanghai, China
1984
Submitted to the Faculty of the
Graduate College ofthe
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
December 1997
COMPARISONS AMONG STOCHASTIC
OPTIMIZATION ALGORITHMS
Thesis Approved:
11
d
ACKNOWLEDGMENTS
First and foremost I am deeply grateful to Dr. John P. Chandler, my supervisor,
for introducing me to such a beautiful area. for his encouragement, patience, and
assistance throughout my ,studies in the Computer Science Department. and for his
guidance in this research and the writing ofthis thesis.
My appreciation extend to my other committee members, Dr. K. M. 'George and
Dr. H. Lu, as welL Their patience and guidance in reviewing this thesis were invaluable..
I am forever ,indebted to my wife Yijun Huang and our daq.ghter Xuejing Chen,
who gave me wholehearted support in order to make this work possible.
11l
'(
TABLE OF CONTENTS
Chapter Page
1. INTRODUCTION , 1
2. SOME STOCHASTIC ALGORITHl\.1S 6
2.1 Random Pattern Search Algorithm ' 6
2.2 Torus Algorithm 7
3. SIMULATED ANNEALING ALGORITHMS 9
3.1 General Simulated Annealing Algorithms ' 9
3.2 Boltzmann Annealing (BA) ' 11
3.3 Fast Annealin~ (FA) 11
1
4. VERY FAST SIMULATED REANNEALING (VFSR) AND THE ASA CODE 13
4.1 Very Fast Annealing " 13
4.2 Reannealing 15
4.3 ASA Code " " 18
5. TEST PROBLEMS 20
6. TEST RESULTS " 25
6.1 Results ofRandom Pattern Search Algorithm 25
6.2 Results ofTorus Algorithm 26
6.3 Results ofSA :.29
16.4 Results ofASA · 29
6.5 Comparisons 32
7. CONCLUSIONS 37
REFERENCES 39
APPENDIX: PROGRAM LIST FOR TORUS.F .43
IV
LIST OF TABLES
TABLE PAGE
TABLE 6-1 RESULTS OF RANDOM SEARCH ALGORITHM 25
TABLE 6-2 RESULTS OF TORUS PROGRAM .27
TABLE 6-3 RESULTS ON OSBORNE FUNCTION 1 28
TABLE 6-4 RESULTS ON OSBORNE FUNCTION 2 .28
TABLE 6..5 RESULTS OF CORANA'S ALGORITHM : .29
TABLE 6-6 RESULTS OF ASA PROGRAM .31
TABLE 6-7 RESULTS ON OSBORNE FUNCTION 1 .31
TABLE 6-8 RESULTS ON OSBORNE FUNCTION 2 .32
TABLE 6-9 RESULTS ON ROSENBROCK FUNCTION .32
TABLE 6-10 RESULTS ON MODIFIED ROSENBROCK FUNCTION 1.. 33
TABLE 6- I I RESULTS ON BOHACHEVSKY FUNCTION 33
TABLE 6-12 RESULT ON POWELL FUNCTION 34
TABLE 6-13 RESULTS ON WOOD FUNCTION 34
TABLE6-14 RESULTS ON BEALE FUNCTION .34
TABLE 6-15 RESULTS ON ENGVALL FUNCTION 35
TABLE 6-16 RESULTS ON OSBORNE FUNCTION 1 35
TABLE 6-17 RESULTS ON OSBORNE FUNCTION 2 36
v
1. Introduction
Global optimization is concerned with the determination of global optima, either
maxima or minima, of a function. Such problems occur frequently in numerous
disciplines which model real world systems. Mathematically, the global optimization
problem can be defined as [15,41]:
Given a set D c 9ln
, and given a real function f: D -) ~, find'
minf(x) .
xeD
Here/is called the objective function, or cost function. Since maximizing/is equivalent
to minimizing to -f, this definition sufficiently includes- the search for global maxima as
well as global minima. Throughout this thesis, the function / will be referred to as, our
cost function, un]ess we mention .otherwise. , I ~.
Besides its importance, global optimization is also an extremely difficult problem.
Various methods have been proposed to solv~ global optimization problems. However,
there are no efficient algorithms which solve aU general global optimization problems. In
general, optimization algorithms can be classified as either stochastic or deterministic'. In
this thesis, we mainly consider the stochastic methods, which evaluate the cost function!
at randomly sampled points from the feasible region D.
There are various stochastic methods [43]. We mainly discuss the simulated
annealing method [3,5,27,30] and its variants [16,42]. Two other stoch~ic optimization
algorithms are also discussed. One is the random pattern search algorithm [28,29],
1
which is one of the early stochastic algorithms. The other is the toros algorithm [39],
which is a recent stochastic algorithm.
Originally, the simulated annealing algorithm was based on the analogy between
the simulation of the annealing of solids and the problem of solving large combinatorial
optimization problems [27]. In condensed matter physics, annealing denotes a physical
process in which a solid in a heat bath ]s heated up by increasing the temperature of the
heat bath to a maximum value at which aU particles ofllie solid to some extent randomly
arrange themselves, followed by cooling thiough slowly lowering the temperature of the
heat bath. In this way, all particles arrange themselves in the low energy ground state of
a corresponding lattice, provided the maximum temperature is sufficiently high and the
cooling is carried out sufficiently slowly. At each temperature -value T, the solid is
allowed to reach thennal equilibrium, characterized by a probability of being in a state
with energy given by the Boltzmann distribution [30,32].
As the temperature decreases, the Boltzmann distribution concentrates on the
states with lowest energy and finally, when the temperature approaches zero, only the
minimum energy states have a non-zero probability of occurrence. However, it is 'well
known [26] that if the cooling is too rapid, i.e., if the solid is not allowed to reach thermal
equilibrium for each temperature,. metastable amorphous structures can be re~ched rather
than the low energy crystalline lattice structure.
The above ideas can be applied to combinatorial optimization. problems in the
following way [27]: First, a new point is randomly sampled. If it generates a new
2
minimum value, the new point is always accepted. If not, it is accepted provided that a
random number between 0 and I is less than a probability defined by a mathematical
function, usually the Boltzmann equation [32]. Early in the iterative process, the
mathematical function generates values near unity, and most points are accepted. By
adjusting a parameter in the probability function, usually referred to as the temperature,
the probability function generates smaller values across successive iterations (cooling),
and eventually, only pomts that produce better solutions are accepted. The procedure
generates approximately optimal solutions for combinatorial problems, such as the
traveling salesman problem, for which exact solutions are presently mathematically
intractable [27].
Bohachevsky et al. [3], Corana et al. [5], and others extended simulated annealing
ideas from combinatorial problems to. the .optimization of functions defined in a
continuous domain. The mathematical model of the simulated annealing algorithm for
continuous optimization based on the ergodic theory of Markov chains can be found in
[8].
8zu and Hartley [42] introduced a fast simulated annealing. Fast anne.aling (FA)
uses the Cauchy distribution, and is often superior to that of Boltzmann annealing. The
fatter tail of the Cauchy distribution allows it to test states farther from the current local
minima during the search process. In addition, fast annealing has an annealing schedule
exponentially faster than the method ofBoltzmann annealing.
3
Boltzmann annealing and fast annealing have distributions which sample infinite
ranges, and there is no provision for considering differences in each parameter-dimension.
Ingber [16] proposed a new probability distribution to accommodate these
desired features. This algorithm is called very fast annealing, and is another variant of
simulated annealing and is exponentially faster than fast annealing. The details are
discussed in Chapter 4.
,
Various other stochastic algorithms have appeared in the literature [43]. Lee [29]
studied a stochastic algorithm, called the random pattern search algorithm, which was
first described by Lawrence and Steiglitz [28]. Rabinowitz [39] presented a stochastic
algorithm called the torus algorithm, for fmding the global optimum of a function of n
variables. The performance of this algorithm was compared to that of the NeIder-Mead
simplex algorithm [34] and the simulated annealing algorithm on a variety of nonlinear
functions.
Many comparisons among different algorithms have been discussed by authors,
for example, see [20,21,24,29,35,39,40]. There is no known optimization algoritlnn
which is better than all other algorithms. For example, if one knows that a cost function
to be optimized is unimodal and smooth everywhere, then a simple Newton iteration [41]
might well be faster than most of the other optimization methods, detez:ninistic or
stochastic. Many stochastic algorithms perform poorly on ill-conditioned smooth
functions, and are most useful on discontinuous, nonsmooth, or multi-modal problems
that are not too ill-conditioned.
4
The main purpose of this thesis is to compare the simulated annealing algorithm,
very fast simulated annealing algorithm, random pattern search .algorithm, and torus
algorithm on a variety of multi-modal, discon$luous, and/or iU-conditioned functions. In
particular, we will test the Osborne functions [36], which are moderately ill-conditioned,
smooth practical problems easily solved by deterministic methods [36]. We have not
seen a satisfactory solution of optimizing Osborne functions given by any of the early
stochastic algorithms [3,5]. Our test results show that these problems can be solved by
the very fast annealing algorithm, the random pattern search algorithm, and the torus
algorithm.
In Chapter 2, the random pattern search algorithm and the torus algorithm are
explained. In Chapter 3, we discuss the general simulated annealing algorithms as well as
Boltzmann annealing and fast annealing. In Chapter 4, we discuss very fast annealing in
detail. In Chapter 5, we list the functions which we use to compare the algorithms. The
comparisons, results, and conclusions are smnmarized in Chapter 6 and Chapter 7 of this
thesis.
5
2. Some Stochastic Algorithms
In the past decade, simulated annealing algorithms have been studied in detail. In
Chapter 3 and Chapter 4, we discuss the simulated annealing method and its variants.
Besides that, there are various other stochastic programming methods proposed by
various authors. There are collected more than two thousands references in [43], though
it is not exhaustive, in this subject. In this chapter we study two stochastic algorithms,
the random pattern search algorithm [28,29] and the torus algorithm [39].
2.1 Random Pattern Search Algorithm
One of the early stochastic optimization algorithms is the random pattern search
algorithm, which was fIrst described by Lawrence and Steiglitz [28]. It was successfully
applied to the optimization ofa variety of chemical engineering problems [12].
Lee [29J studied the random pattern search, and modified the random search
procedure. He tested this procedure on various functions and compared its perfonnance
with the pattern search method [14] and DFP [6,9] method on several aspects, such as the
final value of the objective fimction, the number of iterations of the algorithm, an~ the
number of times the function is evaluated. Recently, Chandler [4] has tried
unsuccessfully to reproduce Lee's results using Lee's program.
The random search algorithm came from a detenninistic method, called pattern
search, which was devised by Hooke and Jeeves [14].
The pattern search algorithm consists of two major phases, an exploratory move
and a pattern move. The exploratory phase, moving from the base point, is designed to
6
explore the local behavior of the objective function. The pattern phase steps along the
approximate negative gradient direction, which is determined from the results of the
exploratory phase.
During the exploratory phase, if there are instances of successive successful
searches, the step size will be extended. Should the next search with this expanded step
size be successful, it is retained; otherwise the step size before extension win be applied
again. If an unsuccessful search is encountered, the step size is decreased; if it becomes
less than some small preset tolerance, then convergence is assumed to be ac~eved.
The random pattern search basically has the same searching procedures as pattern
search, except that it searches in pseudorandom directions uniformly distributed over the
surface of a sphere or hypersphere.
2.2 Torus Algorithm
Rabinowitz proposed a stochastic algorithm, called the torus algorithm, for
fmding the global optimum of a function ofn variables [39].
Three computer functions (Controlling function, Multidimensional function,. and
Single-dimension function) constitute the core of the algorithm. The very detailed
pseudo-codes of these three functions were given (39]. The algorithm is based on using
an adaptive n-dimensional torus to surround and isolate the global minimum. In the
controlling function, an n-dimensional torus moves in n-space and monotonically shrinks
in size over trials, isolating the region containing the global minimum. This shrinkage is
gradual, mimicking slow cooling in a multidimensional function and single-dimension
7
function which are repeatedly called from the controlling function; new points are
randomly sampled around the currently best-fitting point. The pennitted range of these
sampled points shrinks logarithmically over iterations (mimicking rapid cooling).
Detailed descriptions of the relevant user-specified parameters were given. However, it
seems that it would be difficult to implement the algorithm independently [33]. The
original program was written in Common Lisp [44]. In fact, there are some differences
between the pseudo-codes in [39] and the Lisp program [44]. We translated the Lisp
program into a FORTRAN program (see Appendix).
8
3. Simulated Annealing Algorithms
3.1 General Simulated Annealing Algorithms
In general, simulated annealing consists ofthree functional relationships [16]:
1. g(x): The probability density function of the state space x = {xi: ;=1,2, ... , D}.
2. hex): The probability density function for accepting a new value given the just
previous value.
3. T(k): An annealing temperature (T) schedule in annealing-time step k.
General simulated armealing optimization methods choose new points at various
distances from their current point x. Each new point Xnew is generated probabilistically
according to a given distribution g. These algorithms calculate the function value E =
.f{x), and then probabilistically decide to accept or reject it. If accepted, the new point
becomes the current point. The new point may be accepted even if it is worse and has a
larger function value than the current point. The criterion for acceptance is determined by
the acceptance function h, the temperature parameter T, and the difference in the function
values of the two points. Initially, the temperature T is large. As the algorithm
progresses, T is reduced, thus lowering the probability that the acceptance function will
accept a new point if its fimction value is greater than that ofthe current point.
Let Ek = I(x). The acceptance probability is based on the chances of obtaining a
new state with "energy" Ek +\ relative to a previous state with "energy" Ek .
heM) = exp(- Ek+1IT)
exp(-Ek+I/T) + exp(-Ek/T)
9
1 =-----
1+exp(M,/T)
:.:::: exp(-MjT), (1)
where !J.E represents the "energy" difference between the present and previous values of
the energies (considered here as cost functions) appropriate to the physical problem, i.e.,
M =Ek+1 - Ek • 1bis essentially is the Boltzmann distribution contributing to the
statistical mechanical partition function of the system [32]. However, one may choose
another function as the acceptance function [22].
Suppose the function g is given. Let the state-generating probability at the
cooling temperature T(k) at the annealing-time k and within a neighborhood be ;:.:: gk; then
the probability ofnot generating a state in the neighborhood is obviously ~ (1- gk). Our
purpose here is to choose suitable T(k) such that it will suffice to give a global minimum
ofthe cost function. In order to statistically guarantee to obtain the global minimum, we
require that any point in x-space can be sampled infinitely often in annealing-time (IOT).
It suffices to prove that the products of probabilities of not generating a state x lOT for all
annealing-times successive to kO yield zero,
co IT(1- gk) =0 .
k=ku
This is equivalent to
10
(2)
(3)
If the probability density function g is given, then the problem reduces to finding T(k) to
satisfY Equation (3).
3..2 Boltzmann Annealing (BA)
The Boltzmann algorithm chooses a Gaussian probability density function [32],
g(x) = (27CTr.DJ2 exp[-L~.i/(27)], (4)
where ~ = Xnew - x is the deviation ofXnew from the current accepted point x. It has been
proven [10] that it suffices to obtain a global minimum of/if T is selected to be not faster
than
T(k) = To.
Ink
(5)
One can prove that this cooling schedule satisfies Equation (3) in the D-dimensional
neighborhood for an arbitrary size 1Lh"1 and k. In fact we have
~ ~ ~ Lgk
~ :Lexp(- In k) =LYk =00.
~ ~ ~
3.3 Fast Annealing (FA)
(6)
There are sound physical principles underlying the choices of Equations (4) and
(5) [27,32]. It was noted that this method of finding the global minimum in x-space is not
limited to physics examples requiring "temperatures" and "energies". Rather this
methodology can be readily extended to any problem for which a reasonable probability
density can be formulated [42]. It was also noted this methodology' can be readily
extended to use any reasonable generation function g.
11
8zu and Hartley [42] introduced a fast annealing method which uses the following
Cauchy distribution as the generation function:
(7)
The Cauchy distribution has some definite advantages over the Boltzmann form [42]. It
has a "fatter" tail than the Gaussian form ofthe Boltzmann distribution, permitting easier
access far from a local minimum in the search for the desired global minimum.
On the other hand, we can set the annealing schedule as
T(k) = To.
k
(8)
Then one can prove that this cooling schedule satisfies Equation (3). In fact we have
(9)
The method of FA is thus statistically seen to have an annealing schedule exponentially
faster than the method ofBA. This method has been tested on a variety ofproblems [42].
12
4. Very Fast Simulated Reannealing (VFSR) and the ASA Code
In a variety of physical problems, we have a D-dimensional parameter-space.
Different parameters have different finite ranges, fixed by physical considerations, and
different annealing-time-dependent sensitivities. BA and FA have g distributions wlllch
sample infinite ranges, and there is no provision for considering differences in each
parameter-dimension. For example, different sensitivities might require different
annealing schedules.
One might choose a D-product of one-dimensional Cauchy distributions because
the one-dimensional Cauchy distribution has a few quick algorithms. This would also
pennit different Td s to account for different sensitivities
But then we would require an annealing schedule:
T (k) = 1;0
I 1/
kiD
which, although faster than BA, is still quite slow.
(10)
(11)
Motivated by the above considerations, Ingber introduced a very fast simulated
annealing method [16].
4.1 Very Fast Annealing
13
As we mentioned previously, different parameters may have different annealing-time-
dependent sensitivities. We consider a parameter a~ in dimension i generated at
annealing-time k with the range
a~ e[A; ,B;J
The pararneter a~+1 can be calculated from the random variable yi
y; E[-l,l].
Defme the generating function
DID . .
gT =TI I il == TIg~(yl),
;=1 2(y + 1;) In(l + 1/1;} ;=1
(12)
where the subscript i on 1; specifies the parameter index, and the k-dependence in 1; (k)
for the annealing schedule has been dropped for brevity. Its cumulative probability
distribution is
D
Y Y D
GT(y) = r· Jdy'!· ··dy,DgT(Y') == TI G~(yi)" (13)
-I -1 i=1
Gi i =J:. + sgn(y;) In(1 +IiI/I;) (14)
T(Y) 2 2 In(l+ l/r) .
/ is generated from a d from the uniform distribution
ul E U[O, 1],
/ =sng (ui -J:.rZ;[(l+ 1/T,)12ul
-II _1].
2
By a straightforward calculation, one can set the annealing schedule for 1; as
14
(15)
(16)
T; (k) = 1;0 exp (-cikl/D
). (17)
A global minimum statistically can be obtained; that is, the above annealing schedule
satisfies Equation (3):
(18)
It seems sensible to choose control over Ci such that
Ci =m; exp(-n; ID),
(19)
(20)
(21)
where m; and n; can be considered "fr,ee" parameters to help tune ASA for specific
problems. Here ASA refers to Ingber's adaptive simulated annealing code, which we will
explain briefly in Section 4.3.
It has proven fruitful to use the same type of annealing schedule for the
acceptance function h as is used for the generating function g, i.e., Equations (17) and
(19), but with the number of acceptance points, instead of the number of generated points,
used to detennine the k for the acceptance temperature.
In one implementation of this algorithm, new parameters a~+1 are generated from
old parameters a~ by generating the yi until a set ofD is obtained satisfying the range
constraints.
4.2 Reannealing
15
Whenever doing a multi-dimensional search in the course of solving a real-world
nonlinear physical problem, inevitably one must deal with different changing sensitivities
of the ai in the search. At any given annealing-time, it seems sensible to attempt to
"stretch out" the range over which the relatively insensitive parameters are being
searched, relative to the ranges of the more sensitive parameters.
It has proven fruitful to accomplish this by periodically rescaling the annealingtime
k, essentially reannealing, every hundred or so acceptance-events (or at some userdefmed
modules of the number of accepted or generated states), in terms of the
sensitivities s; calculated at the most current minimum value of the cost function, f,
In terms of the largest Si =smax' a default rescaling is perfonned for each k, of
parameter dimension, whereby a new index k: is calculated from each kj •
k: =(In(7;0 / I;k' ) / Ci ) D •
I'iO is set to unity to begin the search, which is ample to span each parameter dimension.
Recall that we use the Boltzmann acceptance criterion as the acceptance criterion.
That is, if
exp(- tJ.f/I'coSl
) > v ,
16
the new point is accepted as the new saved point for the next iteration. Otherwise, the
last saved point is retained. Here Tcost is the ''temperature'' used in this test, and v is from
the unifonn distribution
V E U[O,I].
The annealing schedule for the cost temperature (or, acceptance temperature) is
developed similarly to the parameter temperature. However, the Boltzmann acceptance
criterion uses an exponential distribution which is not as fat-tailed as the distribution used
for the parameters. The index for reannealing the cost function, kcost, is determined by the
nmnber of accepted points, instead of the number of generated points as used for the
parameters.
There is still an unanswered question: How to choose the initial acceptance
temperature? Ingber said: "The initial acceptance temperature is set equal to an initial
trial value off' [16,21]. Here and in the following, we understand that only the quantities
of the temperatures and the function values are compared.
The initial trial value of f is "typically very large relative to the current best
minimum, which may tend to distort the scale of the region currently being sampled".
"Therefore, when this rescaling is perfonned, the initial acceptance temperature is reset to
the maximwn of the most current minimum and the best current minimum ofJ, and the
annealing-time index associated with this temperature is reset to give a new temperature
equal to the minimum of the current cost-function and the absolute values of the current
best and last minima" [16,21].
17
We discussed the above problem with Chandler [4] and he made the following
comments: "However, the cost function f may have units, such as Mlsec., that are
inappropriate for temperature. Further, the scaling (units) of f is arbitrary. Last, all
values of/could be negative. Therefore, these remarks ofIngber's make no sense to me,
although his prescription no doubt will work in most cases."
Also generated are the "standard deviations" of the theoretical forms, calculated
as [?f /(&i)2 r~, for each parameter ai • This gives an estimate of the "noise" that
accompanies fits to stochastic data or functions. At the end of the nm, the off-diagonal
elements of the "covariance matrix" are calculated for all parameters. This inverse
curvature of the theoretical cost function can provide a quantitative assessment of the
relative sensitivity of parameters to statistical errors in fits to stochastic systems.
4.3 ASACode
The adaptive simulated annealing (ASA) code was first developed by Lester
Ingber in 1987 as Very Fast Simulated Reannealing (VFSR) to deal with the necessity of
performing adaptive global optimization on multivariate nonlinear stochastic systems
"[16]. Since 1993, many features have been added, leading to the current ASA code [22].
"Adaptive" in Adaptive Simulated Annealing refers to adaptive options available to a
user to tune the ASA algorithm to optimize the code for application to specific systems.
While the default options may suffice for many applications, this is not intended to imply
.
that the code will automatically adaptively seek the best tuning options. There are many
user options in the ASA code. Among them, there are only a few options which are very
18
influential [22]. However, it seems that it is not easy to grasp all the user options. In our
testing we only use the most influential option "Temperature_Ratio_Scale". For aU other
options we simply use the default values. Of course, we may not obtain the best results
for some problems.
19
5. Test Problems
In this chapter we list the functions which we use to compare SA algorithm,
VFSR algorithm, the random pattern search algorithm, and the torus algorithm. These
functions exhibit moderate ilI-conditioning, nonsmoothness, and multi-modality in
various forms. For the detailed description ofthese functions, the readers are referred to
[3,12,28,35,36,38]. We specify the range for each variable and the starting point for
each function to be minimized. We also state the actual minimwn or approximate
minimum ofeach function.
1. Rosenbrock function:
f(x) =lOO(x2 - xi)2 + (1- XI)2
variable range: -2000 S xp x2 S 2000 [39].
starting point: x = (-1.2, I).
The actual minimum oftrus function is 0 at (1,1).
2. Modified Rosenbrock's function 1 ("flat-groW1d bent knife-edge function"):
f(x) =100lx2 - x~ 1+ (1- X I )2
variable range: -2000 S Xl ,X2 S 2000.
starting point: x = (-1.2, 1).
The actual minimum ofthis function is aat (1,1). This function is not smooth.
3. Modified Rosenbrock's function 2 ("hollow-ground bent knife-edge
function"):
20
variable range: -2000 ~ XI~X2 ~ 2000.
starting point: x =(-1.2, 1).
The actual minimum ofthis function is 0 at (1,1). This function is not smooth.
4. Bohachevsky function:
I(x) = x~ +2x; - 0.3cOS(37lX}) - 0.4cos(47lX2 ) + 0.3+ 0.4
variable range: -2000 ~ Xl 'X2
~ 2000.
starting point: x = (-1, I).
The actual minimum of this function is 0 at (0, 0).
5. Powell function
starting point: X =(3, -1, 0, 1).
The actual minimum ofthis function is 0 at (0,0,0,0). The Hessian matrix of this
function is singular at the minimum.
6. Wood function
I(x) =lOO(x2 - X1
2
)2 +(1- XI )2 +90(x4 - X;)2 +(1- X3)2
+lO.l[(x2_1)2 +(x4 -1)2]+19.8(x2-1)(x4 -1).
starting point: x = (-3, -1, -3, -1).
The actual minimum of this function is 0 at (I, 1, 1, 1).
7. Beale Function
21
variable range: -2000 ~ x],x2
~ 2000.
starting point: x = (1, 0.8).
The actual minimum ofthis function is 0 at (3, 0.5).
8. Engvall function
!(x)=x: +x; +2x:xi -4x] +3.
variable range: -2000 ~ xi ,x2
~ 2000.
starting point: x = (0.5, 2.0).
The actual minimum of this function is 0 at (l, 0).
The following two functions are the Osborne functions [36). Osborne [36]
studied a general method for IIlininrizing a sum of squares which has the property that a
linear least squares problem is solved at each stage and which includes the GaussNewton,
Levenberg, Marquardt, and Morrison methods as particular special cases.
The problem ofminimizing a sum of squares arises naturally from the problem of
determining parameters Xi, i = 1,2, ...,p in the model equation
yet) =F(t, x)
from observations
Yi =y(tJ + &;0 (i =1, 2,...,n ),
where the &i (the experimental errors) are independent, normally distributed random
variables with mean zero and standard deviation cr. In the case n > p the appropriate
maximum likelihood analysis indicates that x should be estimated by minimizing .f(x) =
22
and
n
I(x) = Ilg(x)ll
z
= Lg;(X)2
i=l
This problem will be referred to as the model problem, and it is stressed that we have
offered a statistical justification for minimizing a sum of squares. Osborne's two test
problems are classic practical nonlinear least squares problems.
9. Osborne function 1
In this example, the data values {(t,.,YI)' (l::;; i::;; 33)}, which are given in [36],
are fitted by the model
Osborne's original method is a detenninistic one, which does not have to specify the
range of each variable. Based on the results in each stage of Osborne's problem, we may
choose:
variable range: 0::;; XI ::;; 3, a::;; Xz ~ 3, - 3 ~ x3 ::;; 0, 0 .::;; x4 ::;; 3, 0 ::;; Xs ::;; 3.
starting point: x = (0.5, 1.5, -1, 0.01, 0.02).
The approximate minimum is 0.546E-4 at (0.3753, 1.9358, -1.4647, 0.01287,
0.02212). However, if we choose such a range and starting point, it is very difficult to
obtain the global minimum by using a stochastic method. We cannot figure out what is
the reason for this phenomenon.
23
In order to obtain a global minimum of Osborne function 1 by a stochastic
method, we choose the following ranges and starting point to avoid the possible local
lnImmum.
variable range: 0 ~ Xl ~ 3, - 0.95 ~ x2
~ 1.95, - 3.45 ~ x3
~ -1.45, 0 ~ x4
~ 3,
starting point: (0.5, 1.5, -2, 0.01, 0.02).
10. Osborne function 2
In this example, the model has the fonn
F(t ,x) =Xl exp(-xst) + x2 exp[-x6 (t - X 9 )2]
+ X 3 exp[-x7(t - X IO )2] + X4 exp[-x8 (t - XlI )2]
The data values {(ti 'Yi)' 1~ i ~ 65} are also given in [36]. Osborne function 2 is easier
than Osborne function 1 to m:inimize by stochastic methods. Based on the data in [36],
we choose:
starting point: (1.3,0.65,0.65,0.7,0.6,3,5, 7, 2, 4.5,5.5).
The approximate global minimum is 0.0402 at (1.3100,0.4315,0.6336,0.5993,0.7539,
0.9056, 1.3651,4.8248,2.3988,4.5689,5.6754).
24
, i
6. Test Results
In this chapter we run the random pattern search program, the torus program, the
Corana's SA program, and the ASA program on the test functions we list in the previous
chapter, and. compare the results of these four programs.
6.1 Results ofRandom Pattern Search Algorithm
The random pattern search program in the M.S. report of Daniel Lee does not
solve the modified Rosenbrock 1 and 2 problems, and it is clear that the algorithm used
by Lee should not be able to handle nonsmooth problems. Perhaps Lee used some other
version of his program to solve these two problems.
Lee's algorithm takes random steps only parallel to the coordinate axis. This is
not consistent with the random pattern search algorithms of Lawrence and Steiglitz [28]
and ofBeltrami and Indusi [2], on which Lee's algorithm is supposed to be based.
Chandler [4] has programmed random pattern search [2,28]. Without trying many
more search directions than prescribed, this algorithm should not be able to solve the two
modified Rosenbrock problems, and it does not. The results of this algorithm on the
other t,est problems are shown below.
Table 6-1 Results of Random Search Algorithm
Function NF j{x)
Rosenbrock 443 4.6398068E-9
Bohachevsky 237 2.287497
Powell 8093 1.5708663E-8
Wood 13587 2.5537496E-7
Beale 365 1.1568832E-15
Engvall 299 4.3298698E-14
Osborne1 23151 5.4730076E-5
Osbome2 30403 4.0137737E-2
25
Explanations:
j{x): The minimum value we obtained.
NF: Number of function evaluation.
Random pattern search cannot be recommended in general for constrained or
nonsmooth problems, which are the kinds of problems for which it was designed.
6.2 Results ofTorus Algorithm
The toms algorithm can nm on parallel processors. As the author pointed out, this
approach is more of a Monte Carlo approach, and results are given by conducting
function evaluations in a group to simulate parallel performance, rather than by actually
running on a parallel computational system. For easily comparing with other algorithms,
we only consider one processor.
We set all parameters to the default values except the parameters "scalar2" and
"exit". "Exit" is a stopping criterion with the default value 10.6
. For some problems we
need to set "exit" to be smaller to get more accurate results. "Scalar2" is the weighting
factor for the number of multiple-variable iterations per trial. The user controls the
number of function cans of j{x) made on each pass in multidimensional functions by
.adjusting the "scalar2" parameters. The "scalar2" parameter has the greatest impact on
the outcome, and it should be increased for difficult problems.
There are four stopping criteria in the torus algorithm:
1. The number of trial blocks ,exceeds a pre-defined number "counter2". The
default value of "counter2" is 40.
26
2. The number of successive failures exceeds a pre-defined number "Flagz". The
default value of "F}agz" is 36.
3. The nwnber of consecutive successes exceeds a pre-defined number "flag-count".
The default value of "flag-count" is 24.
4. The difference "last-minimum - best-minimum" is between zero and "exit"
(success). Here, "last-minimum" refers to the smallest values returned by all prior calls
of the single-dimension function, while "best-minimum" refers to the value returned by
the current call ofthe single-dimension function.
The third criterion was not stated in the paper [39], but was actually coded in the
author's Lisp program [44].
We ran our FORTRAN torus program on our test functions and list the results in
the following table:
Table 6-2 Results ofTorus Program
Function scalar2 exit .f{x) NF' Stop Type
Rosenbrock 4 1.0E-6 1.62542E-8 11120 4
Rosenbrockl 5 LOE-6 1.6606£-7 15220 4
Rosenbrock2 6.5 1.0E-6 8.71213E-3 27680 1
Bohachevsky 1 1.0E-6 3.67665E-7 3880 4
Powell 1 1.0E-6 4.51138E-7 9800 4
Wood 8 1.0E-6 3.19705E-7 73400 4
Beale 6 1.0E-6 7.30086E-8 11220 4
Engvall 8 1.0E-6 8.22563E-7 10720 4
Osbomel 14 1.0E-8 5.46544E-5 110450 ' 4
Osbome2 6 1.0E-6 0.0401409 140030 4
I
Explanations:
Rosenbrockl and 2 mean modified Rosenbrock function 1 and 2.
Stop type: Stopping type ofthe problem, as we explained before.
27
From the results in Table 6.2, we see that the torus problem converges for all test
functions except the modified Rosenbrock function 2. The results are satisfactory. In
particular, it converges for both Osborne functions.
Roughly speaking, for certain problem, the larger the value of "scalar2", the larger
the number of function evaluations, and the more accurate the results. However, if
"scalar2" is set too large for a fixed problem, we may not get a more accurate answer, or
may even not get the correct answer. In the following tables we list the results of the
Osborne functions 1 and 2 for different values of "scalar2".
_ Table 6-3 Results on Osborne Function 1
scalar2 j(x) NF stop type
13 5.67855E-5 73000 4
14 5.46544E-5 110450 4
15 5.46579E-5 110650 4
16 5.58468E-5 85400 4
20 5.48989E-5 152000 4
Table 6-4 Results on Osborne Function 2
scalar2 j(x) NF stop type
1 0.0413549 50490 4
2 0.107364 68310 3
3 0.0401507 63910 4
4 0.0401473 94820 4
5 0.0401482 154000 3
6 0.0401409 140030 4
7 0.0401395 231990 4
8 0.0401401 176550 I[ 4
9 0.0401428 187330 4
10 0.041452 183700 4
15 0.0401839 256190 4
20 0.0401398 389290 4
28
From the above tables we see that we get our best results for Osborne function I
and Osborne function 2 when "scalar2" equals 14 or 7, respectively.
6.3 Results ofSA
In [35], Ohm compared Bohachevsky's simulated annealing algorithm and
Corana's simulated annealing algorithm on several test functions. The author concluded
that the results of Corana's program are more accurate than the results of Bohachevsky's
program.
In this section we only test Corana's program on our test functions. The results
are listed in following table. In the table, To and ~ are starting temperature and starting
step vector, respectively.
As the author of [35] suggest,ed, we choose best parameters To and Vo for each
:function. The first four functions were also tested in [35]. Since we use different ranges
ofthe variables, our results are slightly different than the results in [3?l
Table 6-5 Results of Corana's Algorithm
Function ~ Vo NF .f{x)
Rosenbrock 1000 0.01 220000 1.2776709E-8
Rosenhrock1 1000 0.01 216000 8.434226E-7
Rosenbrock2 1000 0.7 368000 5.6360820E-2
Bohachevsky 1000 0.7 180000 1.0778106E-8
Powell 1000 0.01 440000 5.2603830E-7
Wood 1000 0.01 500000 2.228336E-5
Beale 1000 0.01 124000 4.8416447E-9
Engvall 1000 0.01 152000 4.7172453E-8
The results are satisfactory except for the modified Rosenbrock function 2.
6.4 Results ofASA
29
The ASA code and its related documents are updated ·continually by the author,
Lester Ingber. We use the most recent version, Version 15.10, which was released on
June 20, 1997, to test our functions.
To use the ASA code, one has to set up the ASA interface. The program should
be divided into two basic modules. (1) The user calling procedure, containing the cost
function to be minimized, is contained in user.c, user.h and user_cst.h. (2) The ASA
optimization procedure is contained in asa.c and asa.h. The file asa user.h contains
definitions and macros common to both asa.h and user.h. We simply defined our cost
function in user cst.h.
There are many user options in the ASA code, which allow the user to minimize
very different functions. However, we cannot grasp all of the options. One of the very
influential options is Temperature_Ratio_Scale, which determines the scale of parameter
annealing. The default value of Temperature_Ratio_Scale is 10-5
. One may set a larger
value than the default to slow down the annealing, or set a smaller va1ue than the default
to speed up the annealing.
In our test, we set different values of Temperature_Ratio_Scale for different
functions to be minimized. For all other options, we simply use the default values.
For convenience, the author used a Makefile in the ASA code. In Makefile we
set the following options:
DASA TEST = FALSE
DOPTIONS FILE =TRUE
DOPTIONS FILE DATA =TRUE
30
The fitst option tells the program to run our cost function, not the author's test
function. The other two options tell the program to read the parameter values from the
We list our ASA test results for some functions in the following table.
Table 6-6 Results ofASA Program
I Function TRS j{x) NF
Rosenbrock . 0.2 1.136695E-9 59613
Rosenbrock. 0.1 1.598384E-7 34217
Rosenbrockl 0.9 2.144271E-2 132875
Rosenbrock2 0.9 8.271849 125141
Beale 0.1 3.034363E-17 32202
Beale 0.08 6.620662E-18 26832
Beal,e 0.07 1.978987E-18 25760
EngvaU 0.001 4.88498E-15 3704
EngvaU 0.00001 5.52931E-11 1265
Osborne} 0.0001 5.4658E-5 86731
Osborne2 1.0E-1O 0.04013813 312260
In the above table TRS means Temperature_Ratio_Scale.
For some functions the results are satisfactory. In particular., the ASA program
also converges for both Osborne functions.
We also ran the program on the Osborne functions using different values of
Temperature_Ratio_Scale, and list the results in the following tables:
Table 6-7 Results on Osborne Function 1
TRS j{x) NF
1.0E-2 5.474716E-5 268714
1.0E-3 5.464924E-5 147654
1.0E-4 5.4658E-5 86731
1.0E-5 5.465131E-5 59560
1.0E-6 5.468168E-5 33289
1.0E-7 8.742937E-5 12967
31
Table 6-8 Results on Osborne Function 2
TRS fix) , NF
1.0E-5 1.04042136 !l 360069
1.0E-6 0.04017192 247986
I 1.0E-7 0.04020761 203693 i
1.0E-8 0.0402596 148501
1.0E-9 0.04017329 166660
1.0£-10 0.04013813 312260
1.0E-11 0.04017536 68672
1.0E-12 0.0405114 17414
1.0E-13 0.05284764 7848
Roughly speaking, the smaller the value of Temperature_Ratio_Scale, the smaller
the number of function evaluations. However, if we set the value of
Temperature_Ratio_Scale too small or too large for any particular problem, the results
may not be correct.
From the above tables we see that the best TRS for Osborne function 1 is 10'3,
which is larger than the default value, while the best IRS for Osborne function 2 is 10-1°,
which is much smaller than the default value.
6.5 Comparisons
In this section, we compare the results of the four programs.
For all the functions we tested, the simulated annealing program is much slower
than the other three programs. Hence, in the foHowing we mainly compare the results of
the random pattern search program, the torus program, and the ASA program.
Since we do not know very well how to tune the ASA program, we cannot get
satisfactory results for some functions, and we simply omit them in our comparisons.
Table 6-9 Results on Rosenbrock Function
Program j{x) NF
Random Pattern Search 4.6398068E-9 443
32
Torus 1.62542E-8 11120
CoranaSA 1.2776709E-8 220000
ASA 1.136695E-9 59613
For the Rosenbrock function, the fom programs have the similar accuracy. The
random pattern search program is faster than the other three programs, while the Corana
SA program is much slower than others. It seems that the torus program is faster than the
ASA program on the Rosenbrock function. On one hand, the author of [39] tuned his
program mainly based on the Rosenbrock function. On the other hand, we did not tune
all the option parameters in the ASA code to get the best result for the Rosenbrock
function. Actually, we do not know how to tune ASA optimally. These two reasons may
explain why the torus program is superior to ASA on the Rosenbrock function.
Table 6-10 Results on Modified Rosenbrock Function 1
Program j{x) NF
Torus 1.6606E-7 15220
Corana SA 8.434226E-7 216000
For the modified Rosenbrock function I, the torus program and the Corana SA
program have the similar accuracy. The torus program is much faster than the Corana SA
program.
All four programs do not solve the modified Rosenbrock function 2 very well.
Table 6-11 Results on Bohachevsky Function
Program j{x) NF
Random Pattern Search 2.287497 237
Torus 3.67665E-7 3880
CoranaSA 1.0778106E-8 180000
33
The random pattern search program does not solve the Bobachevsky function,
which has several local minimum. It may stops at local minimum. Both the torus
program and Corana SA program so,lve the Bohachevsky function with the similar
accuracy. But the tours program is much faster than the Corana program.
Table 6-12 Result on Powell Function
Program j{x) NF
Random Pattern Search 1.5708663E-8 8093
Torus 4.51138E-7 9800
CoranaSA 5.2603830E-7 440000
For the Powell function, the random pattern search program, the torus program,
and the Corana SA program have the similar accuracy. The first two programs are much
faster than the Corana SA program.
Table 6-13 Results on Wood Function
Program j{x) NF
Random Pattern Search 2.5537496E-7 13587
Torus 3.1 9705E-7 73400
CoranaSA 2.228336E-5 500000
For the Wood function, the random pattern search program and the torus program
are more accurate and much faster than the Corana SA program. The random pattern
search program is even faster than the torus program.
Table 6-14 Results on Beale Function
Program j(x) NF
Random Pattern Search 1.1568832E-15 365
Torus 7.30086E-8 11220
SA 4.8416447£-9 124000
ASA 1.978987E-18 2576u
34
For the Beale function, the results of the random pattern search program and the
ASA program are much more accurate than the results of the other programs. The
random pattern search program is much faster than the other three programs. The torus
program was tuned. to produce median final function values in the range of the other
algorithm, for example, NeIder-Mead simplex method and Corana SA algorithm [39].
The ASA program can produce final function values with arbitrary accuracy. If we adjust
the related parameter in ASA, then ASA can produce a result similar to that of the torus
program, but ASA is faster than the torus program on the Beale function.
Table 6-15 Results on Engvall Function
Program j{x) NF
Random Pattern Search 4.3298698E-14 299
Torus 8.22563E-7 10720
CoranaSA 4.7172453E-8 152000
ASA 5.52931E-ll 1265
For the Engvall function, the random pattern search program is more accurate and
faster than the other three programs. The ASA is much more accurate and much faster
than the torus program.. This may be an example where ASA is more robust than the
torus program.
Table 6-16 Results on Osborne Function 1
Program fix) NF
Random Pattern Search 5.4730076E-5 23151
Torus 5.46544E-5 110450
ASA 5.4658E-5 86731
We are particularly interested in solving both Osborne functions. Except for the
Corana SA program, all other three programs solve the Osborne function 1. The random
35
pattern search program is not aocurate as the other two programs, but is faster than them.
ASA and the torus program have similar accuracy, but ASA is faster than the torus
program.
Table 6-17 Results on Osborne Function 2
Program .f{x) NF
Random Pattern Search 0.040137737 30403
Torus 0.0401409 140030
ASA 0.04013813 312260
Except for the Corana SA program, all other three programs also solve the
Osborne function 2 with similar accuracy. The random pattern search program is faster
than the other two programs.
Finally, we mention again that, for the modified Rosenbrock function 2, all four
programs do not find satisfactory results. Actually, the modified Rosenbrock function 2
is a very difficult function to minimize by a stochastic algorithm [35]. We do not know if
we can solve this problem using ASA, even if we tune it accordingly. ~ortunately,
Rosenbrock 2 is not similar to functions that arise often in practical problems.
Minimizing the Lp-norm of the residuals in a fitting problem, for p=O.5, would yield a
problem similar to the modified Rosenbrock function 2. Values ofp less than 1.0 seem
never to have been used in the statistical literature, as far as we know.
36
:~ . t
I~
''I
7. Conclusions
To solve the global optimization problems, many stochastic optimization
algorithms have been proposed in the past decades. In this thesis, we compare four such
algorithms: the random pattern search algorithm, the torus algorithm, the Corana
simulated annealing (SA) algorithm, and the ASA program. Brief explanations of these
algorithms are given. Ten functions are chosen to test these algorithms. In particular, we
test these algorithms on both Osborne functions.
AIl four programs fail to solve th.e modified Rosenbrock function 2, which is a
very difficult function to minimize by any algorithm. In the remaining comments, we
only consider the other nine functions.
The random pattern search program solved seven functions, but not the modified
Rosenbrock function I and Bohachevsky function. The Corana simulated annealing
program solved the seven functions, but solved neither Osborne function. Both the torus
program and the ASA program solve nine functions.
For all functions it solved, the random pattern search program is faster or much
faster than the other programs.
The torus program, the Corana simulated annealing program, and ASA program
use the annealing principle and can jump out of a local minimum of a function.
Theoretically, simulated annealing algorithm and the very fast annealing algorithm can
find a global minimum of a function, and very fast simulated annealing is much faster
than simulated annealing. The mathematical foundations ofthe torus algorithm. were not
given.
37
In our testing, the Corana simulated annealing program is much slower than the
torus program and the ASA program, as one predicted. Since we do not know very wen
how to tune the ASA program, we failed to solve some of our test functions. For some
functions, like the Rosenbrock function, the torus program is faster than the ASA
program. For some other functions, like the Engvall function, the ASA program is much
faster than the torus program. However, we believe that if one tuned the ASA program
accordingly, the ASA program may be more rebust and faster than the torus program.
In particular, the random pattern search program, the torus program, and the ASA
program, solved both Osborne functions. The random pattern search program is much
faster than the torus program and the ASA program on both Osborne functions. The torus
program is faster than the ASA program on Osborne :function 2, while the ASA program
is faster than the torus program on Osborne function 1.
Suggestions for further study:
Design and/or find a stochastic optimization algorithm to solve the modified
Rosenbrock function 2.
Give a mathematical foundation for the torus algorithm.
Find the reason we should specify the variable ranges when we use the torus
program and the ASA program to test the Osborne function 1.
38
•
REFERENCES
[1] E. Aarts and J. Korst, Simulated Annealing and Boltzmann Machines, John Wiley &
Sons, New York, 1989.
[2] E. J. Beltrami and J. P. Indusi, An adaptive random search algorithmfor constrained
minimization, IEEE Transactions on Computers 21 (1972) 1004-1008.
[3] 1. O. Bohachevsky, M. E. Johnson, and M. L Stein, Generalized annealing for
function optimization, Technometrics 28 {1986) 209-217.
[4] 1. P. Chandler, private communication, 1997.
[5] A. Corana, M. Marchesi, C. Martini, and S. Ridella, Minimizing multimodalfunctions
of continuous variables with the "simulated annealing" algorithms, ACM Transactions
on Mathematical Software 13 (1987) 262-280.
[6] W. C. Davidon, Variable metric method for minimization, ABC Research and
Development Report, ANL-5990, Argonne National Laboratory, Lemont, Illinois, 1959.
[7] M. H. A. Davis, Markov Models and Optimization, Chapman & Hall, London, 1993.
[8] A. Dekkers and E. Arts, Global optimization and simulated annealing, Mathematical
Programming 50 (1991) 367-393.
[9] R. Fletcher and M. 1. D. Powell, A rapid descent method for minimization, Computer
Journal 6 (1963) 163-168.
[10] W. 1. Goffe, G. D. Ferrier, and Jolhn Rogers, Global optimization of statistical
functions with simulated annealing, J. Econometrics 60 (1994) 65-99.
[11] S. Geman and D. Geman, Stochastic relaxation, Gibbs distributions and the
Bayesian restoration of images, IEEE Transactions on Pattern Analysis and Machine
ill Intelligence 6 (1984) 721-741.
[12] M. W. Heuckroth, J. 1. Gaddy, and L. D. Gains, An examination ofthe adaptive
random search technique, AICHE Journal 22 (1976) 744-750.
[13J W. Hock and K. Schittkowski, Test Examples for Nonlinear Programming Codes,
Springer-Verlag, New York, 1980
[14] R. Hooke and T. A. Jeeves, Direct search solution of numerical and statistical
problems, Journal of Association for Computing Machinery 8 (1961) 212-229.
39
[15] R. Horst and P. M. Pardalos, eds., Handbook of Global Optimization, Kluwer
Academic Publications, Boston, 1995.
[16] L. Ingber, Very fast simulated re-annealing, Mathl. Comput. Modelling 12 (1989)
967-973.
[17] L. Ingber, Statistical mechanical aids to calculating term structure models, Phys.
Rev. A 42 (1990 ) 7057-7064.
[18] L. Ingber, Statistical mechanics of neocortical interactions: A scaling paradigm
applied to electroencephalography, Phys. Rev. A 44 (1991) 4017-4060.
[19] L. Ingber, Generic Mesoscopic neural networks based on statistical mechanics of
neocortical interactions, Phys. Rev. A 45 (1992) 2183-2186.
[20] L. Ingber, Simulated annealing: Practice versus theory, Mathl. Comput. Modelling
28 (1993) 29-57.
[21] L. Ingber, Adaptive simulated annealing (ASA): Lessons learned, Control and
Cybernetics 25 (1996) 33-54.
[22] L. Ingber, Adaptive Simulated Annealing (ASA), [ftp.alumni.caltech.edu:
Ipub/ingber/ASA-shar, ASA-shar.Z, ASA.tar.Z, ASA.tar.gz, ASA.zip], Mclean, VA,
Lester Ingber Research.
[23] L. Ingber, H. Fujio, and M. F. Wehner, Mathematical comparison of combat
computer models to exercise data, Mathl. Comput, Modelling 15 (199i) 65-90.
[24] L. Ingber, and B. Rosen, Genetic algorithms and very fast simulated reannealing: A
comparison, Mathl. Comput. Modelling 16 (1992) 87-100.
[25] L. Ingber and D. D. Sworder, Statistical mechanics of combat with human factors,
Math!. Comput. Modelling 15 (1991) 99-127.
[26] L. Ingber, M. F. Wehner, G. M. Jabbour, and T. M. Barnhill, Application of
statistical mechanics methodology to term-structure bond-pricing models, Math!.
Comput. Modelling 15 (1991) 77-98.
[27] S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, Optimization by simulated
annealing, Science 220 (1983) 671-680.
[28] J. P. Lawrence and K. Steiglitz, Randomized pattern search, IEEE Transactions on
Computers 21 (1972) 382-385.
40
[29] D. C. Lee, Computer Experimentation and Comparison of Random Search and
Other Non-linear Optimization Methods, M. S. Report, Computer Science Department,
Oklahoma State University, 1984.
[30] P. J. M. van Laarhoven and E. H. L. Aarts, Simulated Annealing: Theory and
Applications, Kluwer Academic Publishers,. Boston, 1988.
{31] S. F. Masri and G. A. Bekey, A global optimization algorithm using adaptive
random search, Applied Mathematics and Computation 7 (19'80) 353-375.
[32] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller,
Equation of state calculations by fast computing machines, 1. Chern. Phys. 21 (1953)
1087-1092.
[33] M. Mintoff, Algorithm 744: A stochastic algorithm for global optimization with
constraints, Computing Reviews, 37 (1996) 267-268.
[34] J. A. NeIder and R. Mead, A simplex methodfor function minimization, Comput. J.
7 (1965) 308-313.
[35] D. Y. Ohm, Generalized Simulated Annealing for Function Optimization over
Continuous Variables, M. S. Thesis, Computer Science Department, Oklahoma State
University, 1994.
[36] M. R. Osborne, Some aspects of non-linear least squares calculations, Numerical
Methods for Non-Linear Optimization, F. A. Lootsma ed., Academic:: Press, New York,
1971, 171-189.
[37] R. H. J. M. Otten and L. P. P. P. van Ginneken, The Annealing Algorithm, KIuwer
Academic Publishers, Boston, 1989.
[38] J. D. Pinter,. Global Optimization in Action, Kluwer Academic Publications, Boston,
1996.
[39] F. M. Rabinowitz, Algorithm 744: A stochastic algorithm for global optimization
with constraints, ACM Transactions on Mathematical Software 21 (1995) 194-213.
[40] B. Rosen, Function optimization based on advanced simulated annealing, IEEE
Workshop on Physics and Computation, IEEE Press, Dallas, 1992,289-293.
[41] L. E. Seales, Introduction To Non-Linear Optimization, Springer-Verlag, New
York, 1985.
[42] H. Szu and R. Hartley, Fast simulated annealing, Phys. Lett. A 122 (1987) 157-162.
41
[43] URL, http://pc_mally.eco.rug.nl:80Ibiblio/stoprog.html, Stochastic programming
bibliography
[44] URL, http://www.acm.org!contents/journals/toms
42
APPENDIX: Program List for TORUS.F
C
C
C
C
C
C
C
C
C
C
C
*
*
***
*
C
ALGORITHM 744: A STOCHASTIC ALGORITHM FOR GLOBAL OPTIMIZATION
WITH CONSTRAINTS, BY MICHAEL RABINOWITZ. PUBLISHED IN ACM
TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL 21, NO.2, JUNE
1995, PAGES 194-213.
PROGRAM MAIN
THIS MAIN PROGRAN IS A DRIVER OF THE TORUS ALGORITHM.
IMPLICIT REAL*S(A-H,O-Z)
INTEGER N,II,COUNT1,COUNT2,FUNCNU,MAXFUN,SDITER,MDITER,
SDCONT,MDCONT,LOGSNG,LOGMUL,SEED,LP
DOUBLE PRECISION BUMP, EPS, HIT, TRIAL, TORUS, SCALI,
SCAL2,DATAN, START (20) ,UPPER (20) ,LOWER(20),
CUTOFF (20) ,PI,Y1,Y2,SI,TI,CR,DI2,DI4,CUTALT(20},
MINTRS(20) ,CUT(20),RANGE(20) ,NRANGE(20),BSTSET(20),
LSTSET(20) ,BMPSET(20) ,TM1SET(20) ,TMPSET(20) ,TMP2(20) ,
VAR2(20),VARSET(20) ,DELTAS(20),TEMP(20)
COMMON /AO/LP
COMMON /A1/COUNT1,COUNT2,BUMP,EPS,HIT,TRIAL,TORUS,SCAL1,SCAL2
COMMON /A2/FONCNU, MAXFUN, PI
COMMON /A3 / SDITER, MDITER, SDCONT, MDCONT , LOGSNG, LOGMUL
COMMON /COR2/SI,TI,CR,DI2(2) ,DI4(4)
COMMON /OSB/Y1 (33) , Y2 {65)
C
C CALL TRSET TO SET PARAMETER VALUES IN TABLE 1.
C
CALL TRSET
C
C CALL OSBSET TO SET THE PARAMETERS OF THE TWO OSBORNE FUNCTIONS.
C
CALL OSBSET(Y1,Y2)
C
C MAXIMUM NUMBER OF FUNCTIONS TO BE MINIMIZED.
C
MAXFUN = 13
C
C SET THE UNIT NUMBER FOR THE PRINTING.
C
LP = 6
C
PI=4.0DO*DATAN{1.0DO)
C
DO 10 11=1,13
C
C SET THE VALUE OF THE SEED.
C
43
SEED=567~
~O CONTINUE
C
C
**
FUNCNU = II
CALL INIT (N, START, LOWER, UPPER, CUTOFF)
CALJ., CONTRO (SEED, START, N, LOWER, UPPER, CUTOFF, CUTALT, MINTRS ,
CUT, RANGE,NRANGE, BSTSET, LSTSET,BMPSET,TM1SET,TMPSET,
TMP2,VAR2,VARSET,DELTAS,TEMP)
C
STOP
END
SUBROUTINE OSBSET(Y~,Y2)
CC
THIS SUBROlITlNE IS USED TO SET THE PARAMETERS OF THE TWO
C OSBORNE FUNCTIONS.
C
IMPLICIT REAL*8(A-H,O-Z)
INTEGER J
DOUBLE PRECISION Y~(33),Y2(6S),YY1(33) ,YY2(65)
CC
DATA OF OSBORNE FUNCTION 1.
C
DATA YY~/0.844DO,O.908DO,O.932DO,O.936DO/O.925DO,O.908DO,
* O.88~DO,O.B50DO,O.818DO,O.784DO,O.751DO,O.71BDO,O.685D0,
* O.6SBDO,O.628DO,O.603DO,O.S80DO,O.558DO , 0.538DO,O.522D0,
* 0.S06DO,O.490DO,O.478DO,O.467DO,O.457DO,O.448DO,O.438D0,
* 0.43~DO,O.424DO,O.420DO,0.414DO,O.4~IDO,O.406DO/
CC
DATA OF OSBORNE FUNCTION 2.
C
DATA YY2/~.366DO,l.191DO,1.~12DO,I.013DO,O.991DO/0.BB5DO,
* O.831DO,O.847DO,O.786DO,O.725DO,O.746DO,O.679DO,O.608D0,
* O.6S5DO,O.6~6DO,O.606DO,O.602DO,O.626DO,O.6S~DO,O.724D0,
* O.649DO,O.649DO,O.694DO,0.644DO,O.624DO,O.66~DO,0.612D0,
* O.S5BDO,O.533DO,O.495DO,O.SOODO,O.423DO,O.395DO,O.37SD0,
* 0.372DO,O.391DO,0.396DO,O.40SDO,O.42BDO,O.429DO,O.523D0,
* O.562DO,O.607DO,O.6S3DO,O.672DO,O.70BDO,O.633DO,~.668Da,
* 0.645DO,O.632DO,O.S91DO,O.559DO,O.597DO,0.625DO,O.739D0,
* O.710DO,O.729DO,O.720DO,O.636DO,O.581DO,0.428DO,O.292D0,
* 0.162DO,O.098DO,O.OS4DO/
C
DO 10 J=l,33
Yl (J)=YYl (J)
10 CONTINUE
C
DO 20 J=I, 65
Y2 (J) =YY2 (J)
20 CONTINUE
RETURN
END
SUBROUTINE TRSET
CC
THIS SUBROUTINE SET THE DEFAULT VALUES OF THE PARAMETERS IN
C TABJ.,E ~ (PAGE 197) .
CC
THIS PROGRAM CAN SIMULATE THE PARALLEL PROCESSORS.
C HOWEVER, WE RUN OUR PROGRAM ONLY FOR ONE PROCESSOR. HENCE, WE
C SET COUNTl=1 AND SCAL2=4.
CC
THESE DEFAULT VALUES CAN BE ADJUSTED FOR EACH FUNCTION TO BE
C MINIMIZED IN THE SUBROUTINE INIT.
CCFOR MOST OF THE PROBLEMS, WE HAVE TO ADJUST SCAL2. FOR A FEW
44
C PROBLEMS WE NEED TO ADJUST EPS.
C
IMPLICIT REAL*S(A-H,O-Z)
INTEGER COONI'l, COUNT2
DOUBLE PRECISION BUMP, EPS, HIT, TRIAL, TORUS, SCALI, SCAL2
C
C
C
COMMON!A1!COUNT1,COUNT2,BUMP,EPS,HIT,TRIAL,TORUS,SCAL1,SCAL2
BUMP = 5.0D-l
COUNTl = 1
COUNT2 = 40
EPS = 1.OD-06
SCAl..l = 1. ODO
SC1I.L2 = 4.0DO
HIT = 1.5DO
TRIAL 1.SDO
TORUS 4.0D+3
RETURN
END
SUBROUTINE INIT (N,X,LOWER,UPPER,CUTOFF)
CC
IN THIS SUBROUTINE WE SET THE DIMENSIONI START POINT, LOWER
C BOUND, UPPER BOUND, AND CUTOFF FOR EACH FUNCTION TO BE MINIMIZED.
C WE ALSO CHANGE THE DEFAUTE VALUES IN SUBROUTINE TRSET IF NESESSARY.
C HOWEVER, IF THE VALUE OF A PARAMETER IN TRSET IS CHANGED FOR ONE
C FUNCTION, WE HAVE TO RESET THIS VALUE FOR THE NEXT FUNCTION.
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER N, J, MAXFUN, FUNCNU,LP,COUNTI,COUNT2
DOUBLE PRECISION X(20),LOWER(20) ,UPPER(20) , CUTOFF (20) ,PI,
* SI,TI,CR,DI2,DI4,BUMP,EPS,HIT,TRIAL,TORUS,SCAL1,SCAL2
C
COMMON !AO!LP
COMMON !A1/COUNT1,COUNT2,BUMP,EPS,HIT,TRIAL,TORUS,SCALl,SCAL2
COMMON !A2!FUNCNU, MAXFUN, PI
COMMON !COR2/ SI,TI,CR,DI2(2),DI4(4)
C
IF(FUNCNU .LT.1 .OR. FUNCNU .GT. MAXFUN) STOP
GO TO (10,40,70,100,130,160,190,220,250,280,310,340,370) ,FUNCNU
CC
FUNCNU = 1
C ROSENBROCK FUNCTION
C
10 N=2
DO 20 J=l, N
LOWER(J)=-2.0D3
UPPER{J)=2.0D3
CUTOFF(J)=1.0D-7
20 CONTINUE
X(1)=-1.2DO
X (2) =1. 000
WRITE (LP,30)
30 FORMAT(!' ROSENBROCK TEST FUNCTION')
GO TO 400
CC
FUNCNU = 2
C MODIFIED ROSENBROCK FUNCTION WITH OBLIQUE CREASE
C
40 N = 2
DO SO J=l, N
LOWER(J)=-2.0D3
UPPER(J)=2.0D3
CUTOFF(J)=1.0D-7
45
50 CONTINUE
X(l) =-1.2DO
X(2)=1.0DO
CC
SET SCAL,2 IN TRSET.
C
SCAL2 = 5.0DO
WRITE (LP, 60)
60 FORMAT(/' MODIFIED ROSENBROCK TEST FUNCTION WITH
* OBLIQUE CREASE')
GO TO 400
CC
FUNCNU = 3
C MODIFIED ROSENBROCK WITH CUSP
C
70 N=2
DO 80 J=.l, N
LOWER(J)=-2.0D3
UPPER(J)=2.0D3
CUTOFF (J) =1. OD-7
80 CONTINUE
X(l) =-1.2DO
X(2)=1.0DO
CC
SET SCAL2 IN TRSET.
C
SCAL2 = 6.5DO
WRITE (LP, 90)
90 FORMAT(/' MODIFIED ROSENBROCK TEST FUNCTION WITH CUSP')
GO TO 400
C
C FUNCNU = 4
C BOHACHEVSKY FUNCTION
C
100 N=2
DO 110 J=l, N
LOWER(J)=-2.0D3
UPPER(J)=2.0D3
CUTOFF(J)=1.0D-7
110 CONTINUE
X(1) =1. ODD
X(2)=1.0DO
CC
SET SCAL2 IN TRSET
C
SCAL2 = 1.000
WRITE (LP,120)
120 FORMAT(/' BOHACHEVSKY TEST FUNCTION')
GO TO 400
C
C FUNCNU = 5
C OSBORNE 1 FUNCTION
C THE MINIMUM IS APPROXIMATELY F(0.3754,1.9358,-1.4647,O.01287,
C 0.02212)=0.546D-4
130 N=S
C
LOWER (1) =0. ODD
UPPER(l)=3.0DO
LOWER(2)=-O.95DO
UPPER (2) =1. 95DO
LOWER (3) =-3.4500
UPPER (3) =-1.45DO
LOWER(4)=O.ODO
UPPER(4)=3.0DO
LOWER(S)=O.ODO
46
UPPER(5)=3.0DO
C
DO 140 J=l,N
CUTOFF(J)=1.OD-7
1.4 0 CONTINUE
C
X(I)=0.5DO
X(2)=1.5DO
X(3)=-2.0DO
X(4)=1.0D-2
X(5)=2.0D-2
CC
SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 14.0DO
EPS = 1.0D-08
WRITE (LP, 150)
150 FORMAT{j fIX, I TEST FUNCTION OF OSBORNE 1')
GO TO 400
C
C FUNCNU=6
C OSBORNE 2 FUNCTION
CC
THE MINIMUM IS APPROXIMATELY F(I.3100,O.4315,O.6336,
C 0.5993,0.7539,1.3652,4.8248,2.3988,4.5689,5.6754)=0.0402
C
160 N=ll
C
LOWER(I)=O.ODO
UPPER(I)=3.0DO
LOWER (2) =0. ODD
UPPER(2)=3.0DO
LOWER(3)=0.ODO
UPPER (3 ) =3 .. ODO
LOWER(4)=0.ODO
UPPER(4)=3.0DO
LOWER(5)=O.ODO
UPPER(5)=3.0DO
LOWER (6) =0. ODO
UPPER(6)=3.0DO
LOWER(7)=O.ODO
UPPER(7)=5.0DO
LOWER(8)=4.0DO
UPPER(8)=7.0DO
LOWER(9)=0.ODO
UPPER(9)=3.0DO
LOWER{lO)=2.0DO
UPPER(20)=5.0DO
LOWER (IJ.) =3 . DDO
UPPER(1l)=6.0DO
C
DO 270 J=I, N
CUTOFF(J)=1.0D-7
170 CONTINUE
C
X(I)=1.3DO
X(2)=6.5D-2
X(3)=6.SD-2
X {4) =7. OD-I
X(S)=6.0D-I
X(6)=3.DDO
X(7)=5.0DO
X(8)=7.0DO
X(9)=2.0DO
X(10)=4.5DO
47
ii
I,
X(1l)=5.5DO
c
C SET SCAL2 AND EPS IN TRSET
C
SCAL,2 = 6. 000
EPS = 1.0D-06
WRITE (L,P, 180)
180 FORMAT (j I TEST FUNCTION OF OSBORNE 2 I )
GO TO 400
C
C FUNCNU = 7
C CORANA FUNCTION WITH DIMENSION N=2
C
190 N=2
SI = 2.00-1
TI = 5.00-2
CR = 1.50-1
DI2(1) = 1.0D+0
DI2(2) = 1.00+3
DO 200 J=l, N
LOWER (J) =-1. 004
UPPER(J)=1.0D4
CUTOFF(J) =1.00-4
200 CONTINUE
X(l)=1.1D+3
X(2)=1.1D+3
WRITE (LP, 210)
CC
SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 4.0DO
EPS = 1. 00- 06
210 FORMAT (//lX, I CORANA FUNCTION, N=2')
GO TO 400
C
C FUNCNU = 8
C CORANA FUNCTION WITH DIMENSION N=4
C
220 N=4
SI = 2.0D-1
TI = 5.00-2
CR = 1.50-1
DI4 (1) 1. OD+O
DI4(2) = 1.0D+3
DI4 (3) = 1. OD+2
014(4) = l.OD+1
DO 230 J=l, N
LOWER(J)=-1.0D4
UPPER(J)=l.OD4
CUTOFF(J) =1.00-4
230 CONTINUE
X(1)=-1.0D+03
X(2)=1.0D+03
X(3)=-1.0D+3
X(4) =1. OD+03
WRITE (LP,240)
CC
SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 6.5DO
EPS = 1.0D-06
240 FORMAT (//lX, ' CORANA FUNCTION, N=4')
GO TO 400
CC
FUNCNU = 9
48
C POWELL'S SINGULAR TEST FUNCTION
C
C THE MINIMUM IS F(O.O,O.O,O.O,O.O)=O.O
C
250 N=4
X(I)=3.0DO
X(2) =-1. ODO
X(3)=O.ODO
X(4)=1.0DO
C
C SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 1.0000
EPS = 1.0D-06
WRITE (LF, 260)
260 FORMAT (' SINGULAR· TEST FUNCTION OF POWELL I )
DO 270 J=I, N
LOWER(J)=-2.0D3
UPPER(J)=2.0D3
CUTOFF(J)=1.0D-7
270 CONTINUE
GO TO 400
CC
FUNCNU=lO
C WOOD'S TEST FUNCTION
CC
THE MINIMUM IS F(l.0,1.0,1.0,1.0)=0.0
C
280 N=4
X(l)=-3.0DO
X(2) =-1.0DO
X(3)=-3.0DO
X(4) =-1. ODO
CC
SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 8.00DO
EPS = 1.0D-06
C
WRITE (LP ,290)
290 FORMAT (I TEST FUNCTION OF WOOD I )
DO 300 J=1, N
LOWER(J)=-2.0D3
UPPER (J) =2 .. OD3
CUTOFF(J)=1.0D-7
300 CONTINUE
GO TO 400
CC
FUNCNU=ll
C HELICAL VALLEY TEST FUNCTION OF FLETCHER AND POWELL
C THE MINIMUM IS F(l.0,0.O,O.ol=0.0 .
c
310 N=3
X(1) =-1.. ODO
X(2)=0.ODO
X(3)=0.ODO
CC
SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 1. OODO
EPS = 1.0D-06
C
WRITE (LP, 320)
320 FORMAT (' HELICAL VALLEY TEST FUNCTION OF FLETCHER AND POWELL I)
DO 330 J=1, N
49
LOWER (J) =-2. OD3
UPPER(J)=2.0D3
CUTOFF (J) =1. OD-7
330 CONTINUE
GO TO 400
C
C FUNCNU=12
C BEALE'S TEST FUNCTION
C
C THE MINIMUM IS F(3.0,O.5)=0.0
C
340 N=2
X(l)=O.lDO
X(2)=0.lDO
CC
SET SCAL2 AND EPS IN.TRSET
C
SCAL2 = 6.00DO
EPS = 1. OD- 06
C
WRITE(LP,350)
350 FORMAT (I BEALE TEST FUNCTION')
DO 360 J=l, N
LOWER(J)=-2.0D3
UPPER(J)=2.0D3
CUTOFF(J)=1.0D-7
360 CONTINUE
GO TO 400
CC
FUNCNU =13
C ENGVALL TEST FUNCTION
C
370 N = 2
X(l) = S.OD-l
X(2) = 2.0DO
C
C SET SCAL2 AND EPS IN TRSET
C
SCAL2 = 8.00DO
EPS = 1. OD- 06
C
WRITE (LP,380)
380 FORMAT(' ENGVALL TEST FUNCTION')
DO 390 J=l, N
LOWER(J)=-2.0D3
UPPER(J)=2.0D3
CUTOFF(J)=1.0D-7
390 CONTINUE
GO TO 400
CC
PRINT FINAL RESULT
C
400 WRITE (LP, 410)
410 FORMAT (j 21X,' THE INITIAL VALUES ARE ',20X)
WRITE (LP, 420)
420 FORMAT (21X, '====================',20X)
WRITE{LP,430) (X(J) ,J=l,N}
430 FORMAT(' X =',lPG14.6,4G14.6/(4X,5G14.6})
C
RETURN
END
DOUBLE PRECISION FUNCTION TESFNC(N,X)
C
C FUNCTION TESFNC IS USED TO SET THE FUNCTIONS TO BE MINIMIZED
C
50
IMPLICIT REAL*S(A-H,O-Z)
C
INTEGER N, FUNCNU, MAXFUN, KI (10) , INDEX, I, J
DOUBLE PRECISION X(N} ,PI,TEMP(4} ,Y1,Y2,T(65} ,DABS,DATAN,IDINT,
* DCOS,DEXP,DSQRT,R,S,FTX,SI,TI,DI2,DI4,ZI(10) ,CR
C
C
COMMON /A2/FUNCNU, MAXFON, PI
COMMON /COR2/ SI,TI,CR,DI2(2) ,DI4(4)
COMMON /OSB/Y1(33),Y2(65}
IF{FUNCNU .LT.1 .OR. FUNCNU .GT. MAXFUN} STOP
GO TO (10,20,30,40,50,80,120,190,260,270,280,290, 300), FONCNU
C
C FUNCNU = 1
C ROSENBROCK'S TEST FUNCTION
C
10 TESFNC=1.0D2*{X(2)-X(l}**2)**2+(1.ODO-X(1»**2
RETURN
C
C FUNCNU = 2
C MODIFIED ROSENBROCK WITH OBLIQUE CREASE
C
20 TESFNC=1.0D2*DABS(X(2}-X(l}**2)+(1.0DO-X(1»**2
RETURN
C
C FUNCNU = 3
C MODIFIED ROSENBROCK WITH CUSP
C
30 TESFNC=1.0D2*DSQRT(DABS(X(2)-X(1)**2})+(1.ODO-X(1»**2
RETURN
C
C FUNCNU = 4
C BOHACHEVSKY FUNCTION
C
40 TESFNC=X(1)**2+2*X(2}**2-3.0D-l*DCOS(3.0DO*PI*X(l}}
* -4.0D-l*DCOS(4.0DO*PI*X(2»+3.0D-1+4.0D-l
RETURN
CC
FUNCNU=5
C OSBORNE FUNCTION 1
C
50 TESFNC = O.ODO
DO 60 J=l, 33
T(J)=lO.ODO*(J-l)
Ii 0 CONTINUE
DO 70 J=1,33
R=DEXP«-1)*X(4)*T(J}}
S=DEXP({-1)*X(5}*T(J)}
FTX=X(I)+X{2)*R+X(3)*S-Yl(J)
TESFNC = TESFNC+FT'X* *2
70 CONTINUE
RETURN
CC
FUNCNU=6
C OSBORNE FUNCTION 2
C
80 TESFNC=O.ODO
DO 90 J=l, 65
T(J)=O.lDO*(J-l)
90 CONTINUE
DO 110 J=l, 65
TEMP (I) -X(5)*T(J)
TEMP (2) -X(6)*(T{J}-X(9»**2
TEMP(3) -X(7)*(T(J)-X(lO»**2
TEMP (4) -X(B)*{T(J)-X(11)}**2
51
100
110
IF (TEMP (1) . LT. -69)
IF (TEMP (2) . LT. -69)
IF (TEMP (3) . LT. -69)
IF (TEMP (4) .LT. -69)
TEMP(l)=DEXP(TEMP(l»
TEMP(2)=DEXP(TEMP{2»
TEMP (3) =DEXP (TEMP (3) )
TEMP(4)=DEXP(TEMP(4»
FTX=O.ODO
DO 100 1=1, 4
FTX=FTX+X(I) *TEMP (I)
CONTINUE
FTX=FTX-Y2(J}
TESFNC = TESFNC+FTX**2
CONTINUE
RETURN
TEMP (1)
TEMP (2)
TEMP (3)
TEMP (4)
-69
-69
-69
-69
C
C FUNCNU=7
C CORANA FUNCTION, N=2
C
120 TESFNC = O.ODO
INDEX = 2
DO 130 1=1, 2
IF(X(I) .GT. O.ODO) THEN
KI(I)=IDINT(X(I)/SI + O.5D+0}
ELSE IF(X(I) .LT. O.ODO) THEN
KI(I)=IDINT(X(I)/SI - 0.5D+0)
ELSE
KI(I)=O
END IF
IF(KI(I) .EQ. 0) INDEX= INDEX-1
130 CONTINUE
c
IF (INDEX .EQ. 0) THEN
DO 140 I = 1,2
TESFNC = TESFNC+DI2(I)*X(I)**2
140 CONTINUE
GO TO 180
END IF
C
INDEX = 2
DO 150 1=1, 2
IF(DABS(KI(I)*SI - XII»~ .LT. TI) INDEX=INDEX - 1
150 CONTINUE
IF (INDEX .EQ. oj THEN
DO 160 I = 1,2
IF(KI (I) .LT. 0) THEN
ZI(I) = KI(I)*SI + TI
ELSE IF(KI(I) .GT. 0) THEN
ZI(I) KI(I)*SI -TI
ELSE
ZI{I) = O.ODO
END IF
TESFNC TESFNC+CR*DI2{I)*ZI(I)**2
160 CONTINUE
ELSE
DO 1.70 I = 1,2
TESFNC = TESFNC+DI2{I)*X(I)**2
170 CONTINUE
END' IF
1.80 RETURN
CC
FUNCNU=B
C CORANA FUNCTION, N=4
C
52
190 TESFNC = 0.000
INDEX = 4
DO 200 1=1, 4
IF(X(I) .GT. O.ODO) THEN
KI(I)=IDINT(X(I)/SI + 0.5D+0)
ELSE IF{X(I) .LT. 0.000) THEN
KI(I)=IDINT(X(I)/SI - 0.5D+0)
ELSE
KI(I)=O
END IF
IF(KI(I) .EQ. 0) INDEX= INDEX-l
200 CONTINUE
c
IF (INDEX .EQ. 0) THEN
DO 210 I = 1,4
TESFNC = TESFNC+DI4(I)*X(I)**2
210 CONTINUE
GO TO 250
END IF
C
INDEX = 4
DO 220 1=1, 4
IF(DABS(KI(I)*SI - XCI»~ .LT. TI) INDEX=INDEX - 1
220 CONTINUE
IF (INDEX .EQ. 0) THEN
DO 230 I = 1,4
IF(KI(I) .LT. 0) THEN
ZI(I) = KI(I)*SI + '1'1
ELSE 1F(KI(1) .GT. 0) THEN
ZI (I) KI (I) *SI -'1'1
ELSE
ZI(I) = O.ODO
END IF
TESFNC TESFNC+CR*DI4(I)*ZI(I)**2
230 CONTINUE
ELSE
DO 240 I = 1,4
TESFNC = TESFNC+DI4(I)*X(I)**2
240 CONTINUE
END IF
250 RETURN
c
C FUNCNU=9
C POWELL'S SINGULAR TEST FUNCTION
C
260 TESFNC=(X(1)+lO.ODO*X(2»**2+5.0DO*(X(3)-X(4»**2+
* (X(2)-2.0DO*X(3»**4+10.DDO*(X(1)-X(4»**4
RETURN
CC
FUNCNU=lO
C WOOD'S TEST FUNCTION
C
270 TESFNC=100.0DO*(X{2)-X(1)**2)**2+(1.0DO-X(1»**2+
* 90.0DO*(X(4)-X(3)**2}**2+(1.ODO-X(3)}**2+
* 10.lDO* «X(2) -1.0DO)**2+(X(4) -1.000) **2)+
* 19.8DO*(X(2)-1.ODO}*(X(4)-1.0DO)
RETURN
CC
FUNCNU=ll
C HELICAL, VALLEY TEST FUNCTION OF FLETCHER AND POWELL
C
280 R=DSQRT(X(1)**2+X(2)**2)
c
1F(X(l) .EQ.O.ODO) THEN
S=0.25DO
53
ELSE
S=DATAN(X(2)/X(1»/(2.0DO*PI)
ENDIF
C
IF(X(l) .LT.O.ODO) S=S+O.SDO
TESFNC=100. ODO* ( (X (3) -10. ODO*S) **2+ (R-l. ODO) **2) +X {3) **2
RETURN
C
C FUNCNU=12
C BEALE'S TEST FUNCTION
C
290 TESFNC =(1.5DO-X(I)*(1.ODO-X(2»))**2+
* (2.2SDO-X(1)*(1.0DO-X(2)**2»)**2+
* (2.62SDO-X(1)*(I.DDO-X(2)**3})**2
RETURN
C
C FUNCNU = 13
C ENGVALL FUNCTION
C
300 TESFNC = (X(I)**2+X(2)**2)**2-4.0DO*X(1)+3.0DO
RETURN
C
END
DOUBLE PRECISION FUNCTION RAND(SEED)
C
C THIS FUNCTION IS USED TO GENERATE A RANDOM NUMBER BETWEEN -1 AND 1.
C
IMPLICIT REAL*S(A-H,O-Z)
C
BUMP ---- *BUMP*. USED TO DISPLACE A VARIABLE FROM THE BEST-THIS
IS THE CONTROLLING FUNCTION (PAGE 199) .
THE ORIGINAL PROGRAM WAS WRITTEN IN COMMON LISP.
WE TRANSLATE IT INTO A FORTRAN PROGRAM.
VARIABLES:
THE FOLLOWING VARIABLES ARE GIVEN IN TABLE 1 (PAGE 197) :
AND FUNCTIONS:
SET THE PARAMETERS OF THE TWO OSBORNE FUNCTIONS.
SET THE PARAMETERS.
SET THE DIMENSION, LOWER, UPPER, CUTOFF, START
POINT FOR EACH FUNCTION TO BE MINIMIZED.
DEFINE THE FUNCTIONS TO BE MINIMIZED.
GENERATE A RANDOM NUMBER BETWEEN -1 AND 1.
THE DRIVER OF THE MULTIDlMENSION FUNCTION.
MULTIDIMENSION FUNCTION.
SINGER-DIMENSION FUNCTION.
THE INNER LOOP OF THE SINGER-DIMENSION FUNCTION.
INTEGER SEED
SEED = 2045*SEED + 1
SEED = SEED - (SEED/I048576)*1048576
RAND = 2*(SEED+l)/104B577.0 -1.0
RETURN
END
SUBROUTINE CONTRO (SEED, START', N, LOWER, UPPER, CUTOFF, CUTALT,
* MINTRS, CUT, RANGE,NRANGE, BSTSET, LSTSET, BMPSET,TMISET, TMPSET,
* TMP2,VAR2,VARSET, DELTAS, TEMP)
ALGORITHM 744: A STOCHASTIC ALGORITHM FOR GLOBAL OPTIMIZATION
WITH CONSTRAINTS, BY MICHAEL RABINOWITZ. PUBLISHED IN ACM
TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOL 21, NO.2, JUNE
1995, PAGES 194-213.
SUBROUTINES
OSBSET( )
TRSET ()
INIT()
TESFNC ()
RANDO
MULTDM()
MULTI ()
SINGER{)
SNG10
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
54
THE FOLLOWING VARIABLES ARE GIVEN IN TABLE 3 (PAGE 201) :
FITTING VALUE ON THE LAST ITERATION
PROPORTION OF THE RANGE.
*COUNTERI*. THE NUMBER OF PARALLEL PROCESSORS
SIMULAT.
*COUNTER2*. MAXIMUM NUMBER OF TRIAL BLOCKS.
*EXIT*. EXIT CRITERION.
*SCALAR1*. WEIGHTING FACTOR FOR THE NUMBER OF
SINGLE-VARIABLE ITERATIONS PER TRIAL.
*SCALAR2*. WEIGHTING FACTOR FOR THE NUMBER OF
MULTIPLE-VARIABLE ITERATIONS PER TRIAL.
*SHRINK-HIT*. RATE THE TORUS COLLAPSES AFTER A HIT.
*SHRINK-TR.IAL*. RATE THE TORUS COLLAPSES AFTER A
TRIAL BLACK.
*TORUS*. USED TO PARTIALLY DETERMINE THE SIZE OF
THE HOLE IN THE TORUS.
COUNT1
COUNT2
EPS
SCALl
SCAL2
CUTALT
RANGE
MINTRS
TORUS
HIT
TRIAL
START
N
SDI'TER
MOlTER
LSTMIN
BSTSET
INTRAL
INTRCT
COUNT
NEXTVA
FLGDIR
CUT
NRANGE
DOWNUP
BUMPV
BMPSE'T
BSTMIN
BSTSET
FLAGZ
FLGZCT
LAST-MINIMUM.
BEST-SET.
WITHIN-TRIAL.
WITHIN-TRIAL-COUNT.
COUN'I'.
NEXT-VARIABLE.
FLAG-DIRECTION.
CUT.
NEW-RANGE.
DOWN-UP.
BUMP.
BUMP-SET.
LAST-SCORE.
LAST-SET.
FLAGZ.
FLGZCT. (WAS NOT DEFINED IN THE PAPER, BUT WAS
DEFINED IN LISP PROGRAM) .
IMPLICIT REAL*8lA-H,O-Z}
INTEGER N,I,COUNTI,COUNT2,SDITER,MDITER,LP,MOD,
* SDCONT,MDCONT,LOGSNG,LOGMUL,SEED,INTRAL,INTRCT,
* FLAGZ, FLGZCT, FLGDIR, COUNT,NEXTVA, TOTAL,STPVAR
DOUBLE PRECISION BUMP, EPS, HIT, TRIAL, TORUS, SCALI, SCAL2,
* START(N),UPPER(N),LOWER(N),
* CUTOFF(N) , CUTALT(N) ,MINTRS(N) , BSTMIN,
.. CUT(N),RANGE(N) ,NRANGE(N) ,BUMPV,
RANGE. FOR EACH VARIABLE, THE UPPER-LOWER BOUND
SUPPLIED BY THE USER.
CUTOFF-ALT. FOR EACH VARIABLE, THE LARGER OF THE
RANGE/*TORUS* AND THE MINIMUM MEANINGFUL VALUE
SUPPLIED BY. THE USER.
MINlMUM-OUTSIDE-OF-THE-TORUS. IF THE VARIABLE IS AN
INTEGER, THEN TWICE THE MINIMUM SPECIFIED BY THE
USER. IF THE VARIABLE IS A SINGLE FLOAT, THEN 32
TIMES THE MINIMUM SPECIFIED. IF THE VARIABLE IS
A DOUBLE FLOAT, THEN 64 TIMES THE MINIMUM
SPECIFIED. HOWEVER, IN THIS FORTRAN PROGRAM, WE
ONLY CONSIDER THE DOUBLE FLOAT.
START-VALUE. FOR BACH VARIABLE, T.HE STARTING VALUE
SUPPLIED BY THE USER.
NUMBER-OF-VARIABLES. NUMBER OF VARIABLES.
SINGER-DIMBNSION-ITERATIONS. ROOND(10 *SCALAR1*).
MULTIPLE-DIMENSION-ITERATIONS. ROUND(10 *SCALAR2*
TIMES NUMBER-OF-VARIABLES SQUARED) .
PARAMETERS IN CONTROLLING FUNCTION (PAGE 199-200) :
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
55
C
C
BSTSET{N),LSTSET{N),BMPSET(N),TMP2(N) ,
LSTMIN, DFMIN, DATAN, TM1SET (N) , TMPSET (N) ,
VAR2(N) ,VARSET(N) ,DELTAS(N) ,TEMP{N) ,
DLOG, DMAX1, DTEMP
LOGICAL DOWNUP
****
COMMON /AO/LP
COMMON /~/COUNTl,COUNT2,BUMP,EPS,HIT,TRIAL,TORUS.SCAL1,SCAL2
COMMON / A3 / SDITER, MOlTER, SDCONT ,. MDCONT, LOGSNG, LOGMUL
DO 10 1=1, N
MINTRS{I) = 64*CtITOFF{I)
RANGE (I) = UPPER(I)-LOWER(I}
CUTALT(I) = DMAX1(RANGE(IJ/TORUS,CUTOFF(I»
10 CONTINUE
CC
ROUND TO NEAREST INTEGER.
C
SDITER = 10.0DO*SCALI +O.SDO
MDITER = 10.0DO*SCAL2*N**2 + O.5DO
CC
NUMBER OF FUNCTION EVALUATIONS IN EACH CALL OF
C SINGER-DIMENSION FUNCTION.
C
SDCONT = N*COUNT1*SDITER
c
C NUMBER OF FUNCTION EVALUATIONS IN EACH CALL OF
C MULTIDlMENSION FUNCTION.
C
MDCONT = COUNTl*MDITER
CC
SET LOG-CONSTANTS:
C SDITER = LN(SINGLE-DIMENSION-ITERATIONS} (TABLE 5, PAGE 204)
C MDITER = LN(MULTIPLE-DIMENSION-ITERATIONS) (TABLE 4, PAGE 203)
C
DTEMP=SDITER
LOGSNG=DLOG(DTEMP)
DTEMP=MDITER
LOGMUL=DLOG(DTEMP)
CC
PASS START-VALUE TO THE MULTIDIMENSION FUNCTION, AND GET
C THE FIRST LAST-MINIMUM AND THE FIRST BEST-SET.
C
CALL MULTDM{N,START,BSTMIN,BSTSET, LOWER,UPPER,
* RANGE, CUTALT, SEED, TMPSET,VARSET, DELTAS, TEMP)
LSTMIN = BSTMIN
CC
SET THE STARTING VALUES OF EACH VARIABLE IN THE LOOP.
C
INTRAL = 1
INTRCT = 0
FLGZCT = 0
COUNT = 0
CC
IN LISP, THE INDEX OF THE FIRST ELEMENT OF AN ARRAY IS O.
C IN FORTRAN, THE INDEX OF THE FIRST ELEMENT OF AN ARRAY IS 1.
C SO, THE STARTING VALUE OF NEXTVA IS 1, NOT O.
C
NEXTVA = 1
FLGD1R = 0
C
DO 20 1=1, N
BMPSET(I} = BSTSET(I)
LSTSET(I} = BMPSET(I)
CUT (I) = CUTALT (I)
56
: I
NRANGE(I) = RANGE (I)
20 CONTINUE
c
DOWNUP = . FALSE.
C
CALL SINGER (N, BMPSET, BSTMIN, LSTSET, LOWER, UPPER, NRANGE,
* CUT,SEED,NEXTVA,FLGDIR,TMP2,VAR2)
C
IF(BSTMIN .LT. LSTMIN} THEN
FLAGZ 0
ELSE
FLAGZ = ~
END IF
C
TOTAL = MDCONT + SDCONT
CC
LOOP
C
30 CONTINUE
CC
STEP ~. SET THE WITH-TRIAL PARAMETER (STARTING VALUE=1 ) .
C
IF(FLAGZ .EQ. 0 .AND. INTRAL .GT. 1) THEN
INTRAL = INTRAL
ELSE IF (IN'I"RAL . EQ. 3 THEN
INTRAL ~
ELSE
INTRAL = INTRAL + 1
END IF
CC
STEP 2. SET THE WITHIN-TRIAL-COUNT PARAMETER (STARTING VALUE 0).
C
IF (INTRAL .GT. ~) THEN
INTRCT INTRCT + ~
ELSE
INTRCT = 0
END IF
CC
STEP 2+ IN THE PSEUDOCODE [PAGE 199] I THIS STEP WAS MISSING.
C SET THE FLGZCT PARAMETER (STARTING VALUE = 0).
C
IF(FLAGZ .EQ. 0) THEN
FLGZCT FLGZCT + ~
ELSE
FLGZCT = 0
END IF
'I II;'
CC
STEP 3. SET THE COUNT PARAMETER (STARTING VALUE
C
oj .
r
Ii
IF(INTRAL .EQ. ~) COUNT = COUNT + 1
CC
STEP 4. SET THE NEXT-VARIABLE PARAMETER (STARTING VALUE = 1) .
C
NEXTVA = MOD (COUNT, N) + 1
CC
STEP 5. SET THE FLAG-DIRECTION PARAMETER (STARTING VALUE 0) .
C
IF{FLGDIR .EQ. OJ THEN
FLGDIR 1
ELSE
FLGDIR = 0
END IF
C
C STEP 6 AND STEP 7.
C SET THE CUT PARAMETER (STARTING VALUE
57
CUTOFF) .
C
C
SET THE NEW-RANGE PARAMETER (STARTING VALUE RANGE)
DO 40 1=1, N
IF(INTRCT .GT. 1 .AND. FLAGZ .EQ. 0) THEN
CUT(I) = DMAX1(CUT(I)!HIT, CUTOFF(I»
NRANGE(I) = DMAX1(NRANGE(I)!HIT, MINTRS(I»
ELSE IF(INTRAL .EQ. 1) THEN
CUT(I) = DMAX1(CUT(I)/TRIAL, CUTOFF(I»
NRANGE(I) = DMAXl{NRANGE (I) !TRIAL, MINTRS(I»
END IF
40 CONTINUE
C
C STEP 8. SET THE DOWN-UP PARAMETER (STARTING VAL:uE = 0).
C
IF{FLAGZ .EQ. 0 .AND. INTRCT .GT. 1) THEN
IF (BSTSET (NEXTVA) . LT. LSTSET (NEXTVA) J THEN
DOWNUP . FALSE.
ELSE
DOWNUP = .TRUE .
END IF
ELSE IF(INTRAL .EQ. 2) THEN
DOWNUP = .TRUE .
ELSE IF (DOWNUP) THEN
DOWNUP . FALSE.
ELSE
DOWNUP = . TRUE .
END IF
CC
STEP 9. SET THE BUMP PARAMETER (STARTING VALUE 0) .
C
IF (DOWNUP) THEN
BUMPV NRANGE (NEXTVA) * BUMP
ELSE
BUMPV - NRANGE (NEXTVA) * BUMP
END IF
c
C IN LISP PROGRAM, STEP 10 IS AFTER STEP 11 AND STEP 12.
C IT SEEMS IT WORKS BETTER.
CC
STEP 11 AND STEP 12.
C SET THE LAST-MINIMUM PARAMETER (STARTING VALUE = FIRST
C LAST-MIN ).
C SET THE BEST-SET PARAMETER (STARTING VALUE= FIRST
C BEST-SET) .
C
IF{FLAGZ .EQ. 0) THEN
LSTMIN = BSTMIN
DO SO 1=1, N
BSTSET (I) =LSTSET(I}
SO CONTINUE
END IF
CC
STEP 10. SET THE BUMP-SET PARAMETER (STARTING VALUE =
C FIRST BEST-SET)
C
DO 60 1=1., N
BMPSET(I) = BSTSET(I)
60 CONTINUE
IF (INTRAL . NE. 1) THEN
IF(BSTSET(NEXTVA)+BUMPV .LT. UPPER (NEXTVA) .AND.
* BSTSET (NEXTVA) +BUMPV . GT. LOWER (NEXTVA) ) THEN
BMPSET(NEXTVA) = BSTSET(NEXTVA)+BUMPV
ELSE IF(BSTSET(NEXTVA)-BUMPV .LT. UPPER (NEXTVA) .AND.
* BSTSET (NEXTVA) - BUMPV . GT. LOWER (NEXTVA) ) THEN
BMPSET (NEXTVA) = BSTSET (NEXTVA) -BUMPV
58
.'
*
*
*
END IF
END IF
TOTAL ~ TOTAL + SDCONT
C
C STEP l3. SET THE LAST-SCORE AND LAST-SET PARAMETERS (STARTING
C VALUES RETURNED BY A CALL, TO THE SINGLE -DIMENSION FUNCTION) .
C
IF(INTRAL .EQ. l) THEN
CALL SINGER (N,BMPSET, BSTMIN,LSTSET, LOWER,UPPER,NRANGE,
CUT,SEED,NEXTVA,FLGDIR,TMP2,VAR2)
ELSE
CALL MULTDM(N,BMPSET,BSTMIN,TMISET, LOWER, UPPER,
NRANGE,CUT,SEED,TMPSET,VARSET,DELTAS,TEMP)
CALL SINGER(N,TMlSET,BSTMIN,LSTSET,LOWER,UPPER,NRANGE,
CUT, SEED, NEXTVA,FLGDIR,TMP2,VAR2)
TOTAL = TOTAL + MDCONT
END IF
C
C
C
C
C
STEP 14. SET THE FLAGZ PARAMETER (STARTING VALUE = 0 IF NEW
LAST-MIN OBTAINED BY THE FIRST CALL TO THE
SINGLE-DIMENSION FUNCTION OTHERWISE, STARTING VALUE
IF(BSTMIN .LT. LSTMIN) THEN
FLAGZ 0
ELSE
FLAGZ ~ FLAGZ + 1
END IF
1) .
CC
STEP 15.
C
C l.
C 2.
C 3.
C
C 4.
C
EXIT TEST: RETURN TO THE BEGINING OF THE LOOP UNLESS ONE
OF FOLLOWING CRITERIA IS MET:
COUNT=COUNT2 (TOO MANY COMPLETE TRIALS) .
FLAGZ=36 (TOO MANY SUCCESSIVE FAILURES) .
FLGZCT=24 (TOO MANY CONSECUTIVE SUCCESSES) .
THIS THIRD CRITERION WAS NOT STATED IN THE PAPER.
O<LSTMIN - BSTMIN<EPS (SUCCESS).
C
STPVAR ~ 0
IF(COUNT .E.Q. COUNT2) STPVAR 1
IF(FLAGZ .EQ. 36) STPVAR 2
IF(FLGZCT .EQ. 24) STPVAR=3
DFMIN = LSTMIN-BSTMIN
IF(DFMIN .GT. 0 .AND. DFMIN .LT. EPS)
IF (STPVAR .G'T. 0) GO TO 70
GO TO 30
STPVAR 4
C
c
70 WRITE(LP,80)
80 FORMAT(j// 21X, 'THE FINAL RESULT IS ',20X)
WRITE (LP, 90)
90 FORMAT (21X, I~==~==========~=~=~=',20X)
WRITE(LP,100) BSTMIN,DFM.IN
lOO FORMAT (' FUNCTION VALUE =', IG14 .6, 2X, 'DFMIN= ' , IG14 . 6)
WRITE (LP, 110) TOTAL, STPVAR,COUNT,FLAGZ,FLGZCT
110 FORMAT (' NF=', 16, 3X, 'STPVAR=' ,12, 3X, I COUNT=' ,12, 3X,
* ,FLAGZ= I, 12, 3X, 'FLGZCT=' , 12)
WRITE(LP,120) (LSTSET(I) ,I=I,N)
120 FORMAT(' X =',IPGI4.6,4GI4.6/(4X,5G14.6»
RETURN
END
SUBROUTINE MULTDM(N,BMPSET,BSTVAL,BSTSET,LOWER,UPPER,
* NRANGE,CUTALT,SEED,TMPSET,VARSET,DELTAS,TEMP)
CCTHIS SUBROUTINE IS A DRIVER OF THE FOLLOWING MULTIDIMENSION
C FUNCTION. COUNTI IS THE NUMBER OF PARALLEL PROCESSORS SIMULATED.
59
C
IMPLICIT REAL*S (A-H,O-Z,)
INTEGER N,I,J,SEED,COONT1,COUNT2
DOUBLE PRECISION BUMP, EPS, HIT, TRIAL, TORUS, SCALI, SCAL2,
* BMPSET (N) , BSTVAL, BSTSET (N) ,NRANGE (N) , UPPER (N) ,
* LOWER (N) ,CUTALT (N) , TMPBST, TMPSET (N) ,
* VARSET(N),DELTAS(N),TEMP(N)
COMMON /Al/COUNTl., COUNT2, BUMP, EPS, HIT, TRIAL, TORUS, SCALI, SCAL2
C
*
10
20
DO 30 1=1, COUNT1
CALL MULTl(N,BMPSET,TMPBST,TMPSET,LOWER,UPPER,NRANGE,
CUTALT,SEED,VARSET,DELTAS,TEMP)
IF(I .EQ.l) THEN
BSTVAL = TMPBST
DO 10 J=l, N
BSTSET(J) = TMPSET(J)
CONTINUE
ELSE IF(TMPBST .LT. BSTVAL) THEN
BSTVAL = TMPBST
DO 20 J=l, N
BSTSET{J) = TMPSET(J)
CONTINUE
END IF
30 CONTINUE
RETURN
END
C
SUBROUTINE MULTI (N,BMPSET,BSTVAL,BSTSET,LOWER,UPPER,
* NRANGE,CUTALT,SEED,VARSET,DELTAS,TEMP)
IMPLICIT REAL*S(A-H,O-Z)
C
CC
THIS IS A MONTE CARLO ALGORITHM FOR CHANGING THE VALUES OF
C ALL THE VARIABLES. ANNEAL,ING PRINCIPLES ARE USED IN
C CALCULATING THE MAXIMUM RANGE OF THE VARIABLES ON EACH
C ITERATION. FAST COOLING IS IMPLEMENTED AS THE RANGE IS
C LOGARITHMICALLY REDUCED. RANDOMNESS IS INTRODUCED INTO THE
C FIT BY MULTIPLYING THE MAXIMUM RANGE OF THE FIT VARIABLE BY
C A RANDOM NUMBER BETWEEN -1. AND 1 TO COMPUTE THE DEL-TA VALUE
C FOR EACH ITERATION. THE BEST FITTING SCORE AND RELATED
C VARIABLE SET ARE RETURNED BY THE FUNCTION.
C
INTEGER I, N, ITER, FUNCNU, MAXFUN, SDITER,MDITER,
* SDCONT,MDCONT,LOGSNG,LOGMUL,SEED
DOUBLE PRECISION VALUE,TESFNC,RAND,BSTVAL,BMPSET(N),
* BSTSET (N) ,NRANGE (N) , UPPER (N) , LOWER (N) ,
* CUTALT (N) , VARSET (N) , ITERDT ,
* DELTAS(N), TEMP(N),DTEMP
LOGICAL VARFLG,VALFLG
COMMON /A3/SDITER,MDITER,SDCONT,MDCONT,LOGSNG,LOGMUL
C
ITER=O
VARFLG = .FALSE.
VALFLG = . FALSE.
C
DO 10 1=1, N
BSTSET(I) BMPSET(I)
VARSET(I) = BSTSET(I)
10 CONTINUE
CC
LOOP. IN THE FIRST ITERATION (ITER=O), WE SET THE STARTING VALUES
C OF EACH VARIABLES. SOME VARIABLES, SUCH THAT ITERATE-DELTA,
C DELTA, ARE NOT NECESSARILY TO BE INITIALIZED.
C
20 CONTINUE
60
C
C STEP 1: SET THE ITERATE PARAMETER TO COUNT THE NUMBER OF ITERATION
C
IF(.NOT. VARFLG) ITER=ITER+1
IF(ITER .EQ. 1) GO 'TO 50
C
C STEP 2: SET THE ITERATE-DELTA PARAMETER TO WEIGHT THE VALUES
C OF THE VARIABLES.
C
DTEMP=ITER
ITERDT=1.0DO-DLOG(DTEMP}/LOGMUL
C
C STEP 3: SET THE DELTA PARAMETER TO CALCULATE A SET OF BOUNDED
C STOCHASTIC VALUES BY WHICH THE VALUES IN VARIABLE-SET
C CAN BE MODIFIED.
C THERE IS A·SLIGHT DIFFERENCE BETWEEN PSEUDOCODE AND
C LISP PROGRAM.
C
DO 30 1=1, N
DELTAS (I) = I TERDT*NRANGE (I) *RAND (SEED)
IF (DABS (DELTAS (I)} .LT. CUTALT(I» THEN
DELTAS (I) = 4 * CUTAL~(I) * RAND (SEED)
IF(DABS(DELTAS(Ill .LT. CUTALT(I)} THEN
IF (DELTAS (I) .LT. 0) THEN
DELTAS (I) -CUTALT (I)
ELSE
DELTAS (I) = CUTALT (I)
END IF
END IF
END IF
CC
STEP 4: SET THE VARIABLE-SET PARAMETER TO CALCULATE A SET
C OF BOUNDED STOCHASTIC VALUES BY WHICH THE VALUES
C IN VARIABLE-SET CAN BE MODIFIED.
C
TEMP (I) = BSTSET(I) + DELTAS (I)
IF (TEMP (I) .GT;LOWER(I) .AND. TEMP (I) .LT.UPPER(1») THEN
VARSET(I) = TEMP (I)
ELSE
VARSET(I) = BSTSET(I)
END IF
30 CONTINUE
CC
STEP 5: SET THE VARIABLE-FLAG PARAMETER TO FLAG IF THE VALUES
C IN VARIABLE-SET AND VALUE-SET ARE IDENTICAL FOR ALL
C VARIABLES.
C
DO 40 1=1, N
IF (VARSET(I) .NE. BSTSET(I» THEN
VARFLG = . FALSE.
GO TO 50
END IF
40 CONTINUE
CC
IF THE VALUES IN VARIABLE-SET AND VALUE-SET ARE IDENTICAL FOR
C ALL VARIABLES, THE FOLLOWING STEPS ARE NOT NECESSARY.
C
VARFLG = . TRUE.
GO TO 20
CC
STEP 6: SET THE VALUE PARAMETER TO COMPUTE A VALUE FOR EACH
C FUNCTION CALL.
C
50 VALUE = TESFNC(N,VARSET)
C
61
C STEP 7-STEP 9:
C STEP 7: SET THE VALUE-FLAG PARAMETER TO FLAG IF A NEW MINIMUM
C OBTAINED.
C STEP S: SET THE BEST-VALUE PARAMETER TO STORE THE MINIMUM
C VALUE OBTAINED ACROSS ITERATIONS.
C STEP 9: SET VALUE-SET PARAMETER TO STORE THE VALUES OF 'THE
C VARIABLE SET THAT GENERATES THE BEST VALUE.
C
IF (ITER . EQ. 1) BSTVAL = VALUE
IF (VALUE . LT. BSTVAL) THEN
VAL,FLG = . TRUE .
BSTVAL = VALUE
DO 60 1=1, N
BSTSET(I) VARSET(I)
60 CONTINUE
END IF
CC
STEP 10: EXIT TEST: RETURN TO THE BEGINING OF THE LOOP UNLESS
C MULTIPLE-DIMENSION-ITERATIONS = ITERATE.
C
IF{ITER .LT. MDITER) GO TO 20
C
RETURN
END
C
SUBROUTINE SINGER(N,START,BSTVAL,BSTSET,LOWER,UPPER,
* NRANGE, CUT, SEED, NEXTVA, FLGDIR, TMP2, VAR2)
CC
THIS IS A MONTE CARLO ALGORITHM FOR CHANGING THE VALUES OF
C ONE VARIABLES. ANNEALING PRINCIPLES ARE USED IN
C CALCULATING THE MAXIMUM RANGE OF THE VARIABLES ON EACH
C ITERATION. FAST COOLING IS IMPLEMENTED AS THE RANGE IS
C LOGARITHMICALLY REDUCED. RANDOMNESS IS INTRODUCED INTO THE
C FIT BY MULTIPLYING THE MAXIMUM RANGE OF THE FIT VARIABLE BY
C A RANDOM NUMBER BETWEEN -1 AND 1 TO COMPtITE THE DELTA VALUE
C FOR EACH ITERATION. THE BEST FITTING SCORE AND RELATED
C VARIABLE SET ARE RETURNED BY THE FUNCTION.
CC
THIS IS THE OUTER LOOP OF THE SINGER-DIMENSION FUNCTION.
C
IMPLICIT REAL*S(A-H,O-Z)
C
C
C
C
INTEGER I,J,N,COUNT1,COUNT2,CNT1,CNT2,NEXTVA,FLGDIR,MOD,
* SDITER,MDITER, SDCONT,MDCONT, LOGSNG,LOGMUL, SEED
DOUBLE PRECISION BUMP, EPS, HIT, TRIAL, TORUS, SCAL1, SCAL2,
* BSTVAL, START (N) , BSTSET (N) , NRANGE (N) , VAR2 (N) ,
* UPPER(N), LOWER(N) ,CUT(N) ,TMP2(N) ,TMPBST
LOGICAL VARFLG,VALFLG
COMMON /A1/COUNT1,COUNT2,BUMP,EPS,HIT,TRIAL,TORUS,SCAL1,SCAL2
COMMON /A3/SDITER,MDITER,SDCONT,MDCONT,LOGSNG,LOGMUL
DO 50 J=l, COUNT1
DO 20 CNTl=l, N
IF( FLGDIR .EQ. 0) THEN
CNT2 MOD ( (NEXTVA+CNT1), N) + 1
ELSE
CNT2 = MOD ( (NEXTVA-CNT1+N), N) + 1
END IF
CALL SNG1{CNT2,N,START,TMPBST,TMP2,LOWER,UPPER,
* NRANGE,CUT,SEED,VAR2)
DO 10 1=1, N
START(I)=TMP2{I)
10 CONTINUE
62
20 CONTINUE
IF(J .EQ.1) THEN
BSTVAL = TMPBST
DO 30 1=1, N
BSTSET(I) = TMP2(I)
30 CONTINUE
ELSE IF (TMPBST . LT. BSTVAL) THEN
BSTVAL = TMPBST
DO 40 1=1, N
BSTSET{I) = TMP2{I)
40 CONTINUE
END IF
SO CONTINUE
RETURN
END
SUBROUTINE SNG1(CNT2,N,START,BSTVAL,BSTSET,LOWER,UPPER,
* NRANGE,CUT,SEED,VAR2)
C
C THIS IS THE INNER LOOP OF THE SINGLE-DIMENSION FUNCTION. IT IS
C IDENTICAL TO THE MULTIDIMENSION FUNCTION EXCEPT THAT LOG-CONSTANT
C IS SET TO LN(SINGLE-DIMENSION-ITERATIONS), AND ONLY THE VARIABLES
C DESIGNATED BY CNT2 IS CHANGED ON EACH PASS.
C
IMPLICIT REAL* 8 (A-H,O-Z)
C
INTEGER I,N,CNT2,ITER,FUNCNU,MAXFUN,SDITER,
* MDITER,SDCONT,MDCONT,LOGSNG,LOGMUL,SEED
DOUBLE PRECISION VALUE, TESFNC, RAND/ BSTVAL, START (N) ,
* DABS, BSTSE'T IN) ,NRANGE (N) , UPPER (N) , LOWER (N) ,
* DLOG,CUT (N) ,VAR2(N) ,ITERDT/DELTA,TEMP,DTEMP
LOGICAL VARFLG,VALFLG
c
COMMON /A3jSDITER,MDITER/SDCONT,MDCONT,LOGSNG,LOGMUL
C
ITER = 0
VARFLG = . FALSE.
VALFLG = . FALSE.
C
DO 10 1=1, N
BSTSETII) = START (I)
VAR2(I} = BSTSET(I)
10 CONTINUE
C
20 CONTINUE
1F(.NOT. VARFLG) ITER = ITER + 1
IF(ITER .EQ. 1) GO TO 30
c
C
C
DTEMP=ITER
ITERDT=1.0DO-DLOG(DTEMP)/LOGSNG
DELTA = ITERDT * NRANGE(CNT2} * RAND(SEED)
IF((DABS(DELTA» .LT. CUT(CNT2)} THEN
DELTA = 1.6 * CUT (CNT2) * RAND (SEED)
IF(DABS(DELTA» .LT. CUT(CNT2}) THEN
IF(DELTA .LT. 0) THEN
DELTA -CUT (CNT2)
ELSE
DELTA = CUT (CNT2 )
END IF
END IF
END IF
TEMP = BSTSET(CNT2) + DELTA
IF (TEMP.GT.LOWER(CNT2) .AND. TEMP.LT.UPPER(CNT2}) THEN
VAR2(CNT2) = TEMP
63
C
C
C
C
ELSE
VAR2(CNT2) = BSTSET(CNT2)
END IF
IF (VAR2 (CNT2) .NE. BSTSET(CNT2)) THEN
VARFLG = . FALSE.
GO TO 30
END IF
VARFLG = .TRUE .
GO TO 20
30 VALUE = TESFNC(N,VAR2)
IF(ITER .EQ. 1.) BSTVAL = VALUE
IF (VALUE . LT. BSTVAL) THEN
VALFLG = . TRUE .
BSTVAL = VALUE
BSTSET(CNT2) = VAR2(CNT2)
END IF
IF(ITER .LT. SDITER) GO TO 20
RETURN
END
64
VITA
Debao Chen
Candidate for the Degree of
Master of Science
Thesis: COMPARISONS AMONG STOCHASTIC OPTIMIZATION ALGORITHMS
Major Field: Computer Science
Biographical:
Education: Graduated from Shanghai Second Institute of Education, China, and
Shanghai Institute of Education, China, in July 1980 and July 1984,
respectively.. Graduated from Fudan University, China, in January 1986;
received a diploma in graduate program in Mathematics for Assistant
Professors. Graduated from The University of Texas at Austin in August
1994 and received a Philosophy of Doctor degree in Mathem?ltics. Completed
the requirements for the Master of Science degree with a major in Computer
Science at Oklahoma State University in December 1997.
Professional Experienc,e: Employed as an assistant professor by Shanghai Second
Institute of Education, China, August 1980 to August 1988. Employed as a
graduate teaching assistant and a graduate research assistant by the
Mathematics Department, The University of Texas at Austin, August 1988 to
August 1994. Employed as an assistant professor by the Mathematics
Department, Oklahoma State University, August 1994 to August 1995.
Professional Membership: America Mathematical Society.