Page 346 - Applied Numerical Methods Using MATLAB

UNCONSTRAINED OPTIMIZATION 335
irregularities and defects in the crystal structure, analogous to being too hasty
to find the global minimum. The simulated annealing process can be imple-
mented using the Boltzmann probability distribution of an energy level E(≥0)
at temperature T described by
p(E) = α exp(−E/KT ) with the Boltzmann constant K and α = 1/KT
(7.1.21)
Note that at high temperature the probability distribution curve is almost flat over
a wide range of E, implying that the system can be in a high energy state as
equally well as in a low energy state, while at low temperature the probability
distribution curve gets higher/lower for lower/higher E, implying that the system
will most probably be in a low energy state, but still have a slim chance to be
in a high energy state so that it can escape from a local minimum energy state.
The idea of simulated annealing is summarized in the box below and cast
into the MATLAB routine “sim_anl()”. This routine has two parts that vary
with the iteration number as the temperature falls down. One is the size of step
x from the previous guess to the next guess, which is made by generating a
random vector y having uniform distribution U[−1, +1] and the same dimension
−1
as the variable x and multiplying µ (y) (in a termwise manner) by the difference
vector (u − l) between the upper bound u and the lower bound l of the domain
−1
of x.The µ -law
(1 + µ) |y| − 1
−1
g (y) = sign (y) for |y|≤ 1 (7.1.22)
µ
µ
implemented in the routine “mu_inv()” has the parameter µ that is increased
according to a rule

µ = 10 100 (k/k max ) q with q> 0: the quenching factor (7.1.23)
as the iteration number k increases, reaching µ = 10 100 at the last iteration k =
k max . Note the following:

ž The quenching factor q> 0 is made small/large for slow/fast quenching.
−1
ž The value of µ -law function becomes small for |y| < 1as µ increases
(see Fig. 7.7a).

The other is the probability of taking a step x that would result in change
f > 0 of the objective function value f(x). Similarly to Eq. (7.1.21), this is
determined by

q
k f

p(taking the step x) = exp − for f > 0 (7.1.24)
k max |f(x)|ε f