Global Optimization Methods

based on Tabu Search

Dissertation zur Erlangung des

naturwissenschaftlichen Doktorgrades
der Julius-Maximilians-Universität Würzburg

vorgelegt von

Svetlana Stepanenko

aus Ivanovo

Würzburg 2008
Eingereicht am: 30.10.2008________________________________________________

bei der Fakulät für Chemie und Pharmazie

1. Gutachter: Prof. Dr. Bernd Engels_______________________________________

2. Gutachter: PD Dr. Reinhold Fink________________________________________

der Dissertation

1. Prüfer: Prof. Dr. Bernd Engels_______________________________________

2. Prüfer: PD Dr. Reinhold Fink________________________________________

3. Prüfer: Prof. Dr. Christoph Sotriffer___________________________________

des Öffentlichen Promotionskolloquiums

Tag des Öffentlichen Promotionskolloquiums: 05.12.2008_________________________

Doktorurkunde ausgehändigt am: __________________________________

dedicated to the memory of my father,
Nikolay Valentinovich Kuznetsov
Die vorliegende Arbeit wurde unter Anleitung von Prof. Dr. Bernd Engels von Mai 2004 bis
Oktober 2008 am Institut für Organische Chemie der Julius-Maximilians-Universität
Würzburg angefertigt.

Teilergebnisse dieser Arbeit waren Gegenstand von Publikationen333,338 sowie von Postern
und Kurzvorträgen.
Acknowledgements
The preparation of this thesis would not have been succeeded the way it did without
the support of many people contributed directly or indirectly. First of all, I extend my
gratitude and appreciation to my supervisor, Prof. Dr. Bernd Engels, whose suggestions, help,
valuable insights, motivation, and encouragement helped me in all the time of research and
writing of this thesis, whose comments and criticism ensure the quality of it. My sincerest
thanks also go to Priv.-Doz. Dr. Reinhold Fink, who kindly accepted the task of reviewing the
present thesis, for his inspiring conversations, for the time he invested in helping me with all
sorts of theoretical chemistry problems. I thank Prof. Dr. K. Baumann for fruitful discussions.
My additional thanks go to all members of the group (past and current) for all their help,
support, interest and helpful hints, for a friendly and humorous climate.
Ongoing financial support by the Deutsche Forschungsgesellschaft in the framework
of the SPP 1178 and the SFB 630 is gratefully acknowledged.
Finally, I would like to express my dearest thanks to my parents, Nikolay Kuznetsov
and Nadezhda Kuznetsova, who supported me for better for worse, for guidance in sending
me down the right path in my education and in life, for their continuous love, encouragement,
and understanding during my life. I wish to thank my family, Oleg and Vladimir Stepanenko
for their insight, support, and patience, which allowed me to concentrate and finish my work
on this thesis. Thank you for always standing by my side in all situations and believing in me
when I didn’t believe in myself.
“Die Phantasie arbeitet in einem schöpferischen Mathematiker nicht weniger als in einem
erfinderischen Dichter.“
Jean-Baptist le Rond D'Alembert (1717 - 1783)

“Es ist nicht genug, zu wissen, man muß auch anwenden; es ist nicht genug, zu wollen, man
muß auch tun.“
Johann Wolfgang von Goethe (1749-1832)
Table of Contents

Chapter 1 Introduction ...................................................................................................... 3

1.1 Global Optimization (GO) Problem........................................................................... 4
1.2 Heuristics.................................................................................................................... 6
1.2.1 Metaheuristic Features ....................................................................................... 7
1.2.2 Population metaheuristic .................................................................................. 10
1.2.2.1 Swarm methods ............................................................................................ 10
1.2.2.1.1 Ant Colony Optimization (ACO)............................................................. 11
1.2.2.1.2 Particle Swarm Optimization (PSO) ........................................................ 14
1.2.2.2 Genetic Algorithms (GAs) ........................................................................... 18
1.2.2.3 Scatter Search (SS)....................................................................................... 22
1.2.3 Simulated Annealing (SA) ............................................................................... 25
1.2.4 Variable Neighbourhood Search (VNS) .......................................................... 30
1.2.5 Tabu Search (TS) ............................................................................................. 37
1.3 Conclusions .............................................................................................................. 46

Chapter 2 TS based methods........................................................................................... 48

2.1 Gradient Tabu Search............................................................................................... 48
2.1.1 Algorithm description ...................................................................................... 48
2.1.1.1 Modest Ascent and Local Search Strategy................................................... 50
2.1.1.2 Diversification Strategy................................................................................ 55
2.1.2 TS Memory Elements....................................................................................... 56
2.1.3 Implementation and Experiments..................................................................... 59
2.1.3.1 Parameters .................................................................................................... 60
2.1.3.2 Numerical Results ........................................................................................ 63
2.1.4 Conclusions ...................................................................................................... 70
2.2 TSPA and GOTS...................................................................................................... 70
2.2.1 Description of algorithms................................................................................. 71
2.2.2 Investigation of the efficiency of the new approaches..................................... 75
2.2.2.1 Parameters .................................................................................................... 76
2.2.2.2 Tests of the efficiency .................................................................................. 80
2.2.3 Conclusions ...................................................................................................... 91
2.3 Conclusions .............................................................................................................. 92

-1-
Chapter 3 Application and Discussion ........................................................................... 94
3.1 Introduction .............................................................................................................. 94
3.1.1 Conformational search techniques ................................................................... 95
3.2 GOTS application................................................................................................... 106
3.2.1 Variable choice and numbering rules............................................................. 106
3.2.2 Adaptation of the GOTS to the conformational search.................................. 109
3.2.3 Comparing the influence of different parameters .......................................... 115
3.3 Experimental results............................................................................................... 117
3.3.1 Conformational studies of amino acids.......................................................... 121
3.3.2 Conformational studies of acetylcholine........................................................ 126
3.3.3 Conformational studies of ACE and HIV-1 protease inhibitors. ................... 128
3.4 Conclusions ............................................................................................................ 133

Chapter 4 Summary ....................................................................................................... 134

Chapter 5 Zusammenfassung........................................................................................ 139

Appendix A Unconstrained Test Problems....................................................................... 145

Appendix B Notation ..................................................................................................... 148

References ............................................................................................................................. 152

-2-
Chapter 1 Introduction

One of the foundations of our world is the search for an optimal state. Since human
beings exist, we exert for improvement. We have a strong wish to reach the maximum success
with minimal efforts. Therefore, optimization is one of the oldest skills which even extends
into daily life 1 and many problems in science, engineering, business and economics, such as
acoustics equipment design, cancer therapy planning, chemical process modelling, data
analysis, classification and visualization, economic and financial forecasting, environmental
risk assessment and management, industrial product design, laser equipment design, model
fitting to data, optimization in numerical mathematics, optimal operation of “closed”
engineering or other systems, packing and other object arrangement problems, portfolio
management, potential energy models in computational physics and chemistry, process
control, robot design and manipulations, systems of nonlinear equations and inequalities, and
waste water treatment systems management, travelling salesman problems (TSP) 2 , the global
optimization of Artificial Neural Networks (ANNs) 3 and telecommunication networks 4 and
some applications appeared also in questions of the chemical industry5,6 , can be formulated as
computing globally optimal solutions 7,8,9,10,11 . In the case of a complex nonlinear system the
associated decision model will in general have an enormous amount of local optima whose
number is normally unknown. Typically, most of the local solutions are also unacceptable as
compared with the global one. Therefore, general local optimization strategies are not
applicable to the problems. Instead, a global search approach is required and one needs
appropriate global optimization (GO) ideas and techniques.
The objective of GO is to find the best solution of a created mathematical model which
corresponds to the global minimum (or maximum) of a suitable objective function within a
given collection of feasible constraints. The objective function is a characteristic property of
the system, such as profit, utility, damage, risk or error. Constrains may be given by existence
conditions of the physical, technical, economic or some other system.
Due to the complexity of many optimization problems, particularly of high dimensional
ones encountered in most practical settings, exact algorithms generally perform very poorly.
Actually, metaheuristics9,12,13,14,15 is conspicuously preferable in practical applications and it

-3-
was shown to obtain highly accurate solutions in many cases. The power of metaheuristics is
robustness and success for a wide range of problems. It is generally the method of choice if:
• Calculation of the objective function is very expensive or time consuming.
• The exact gradient of the objective function cannot be computed, or its numerical
approximation is very expensive or time consuming.
• The values of the objective function may contain noise.
• The search space is very large.

If obtaining any feasible solution is not sufficient and the quality of solutions is critical, it is
very important to investigate effective procedures to obtain the best possible solutions within
a given time. In this thesis recent developments of intelligent search methods like Tabu
Search and its application to problems arising in the chemical area are discussed and a new
approach applicable to chemical problem is developed.
The purpose of this chapter is to introduce basic heuristic concepts of approaches that
generate feasible solutions and to show how they can be applied to combinatorial
optimization problems. In Chapter 2 new nonlinear global optimization routines based on the
Tabu Search strategy are described. They try to determine the global minimum of a function
by the steepest descent - mildest ascent strategy. The new algorithms are explained and their
efficiencies are compared with other approaches by determining the global minima of several
well-known test functions with varying dimensionality. This includes an investigation about
the influence of user-defined parameters. The efficiency of the new approaches is also studied
by comparisons with other approaches for the test cases. In Chapter 3 one of our methods is
adapted to conformational search problems. It is tested by locating the global minimum
energy conformation of amino acids, two angiotensin converting enzyme (ACE) inhibitors 16 ,
2-acetoxy-N,N,N-trimethylethanaminium, and HIV-1 protease inhibitor 17 . The last chapter
summarized this work.

1.1 Global Optimization (GO) Problem

GO aims to find optimal configurations for a given problem. Thus, it is first necessary
to define what an optimum is 18 . In the case of an objective function F, an optimum is either a
maximum or a minimum. To formulate the problem of the global optimization it is required
that F is a function of the variable vector xi, the component-wise bounds xl and xr related to
the decision variable vector xi are finite, and the feasible set D is nonempty. If this is fulfilled,
the extreme value theorem 19 guarantees that the global optimization model is a well-posed
problem. Formally, a general optimization problem can be written as

-4-
 min F(xi ) or max F(xi )
(1)
subject to xi ∈ D
where F(xi) is the continuous objective function with D ={xi : lower bound ≤ xi ≤ upper
bound}. Depending on the problem xi are binary, integer, or continuous variables. A well-
known example from computational chemistry is the conformational search for a large and
very flexible molecule. This task comprises all typical ingredients of an optimization problem.
It possesses a large number of possible solutions with similar quality and the handling of the
problem necessitates the scan over a large space. Finally, as in many optimization problems of
econometrics, it remains uncertain if the optimal solution was really found.

-100
Local maximum

-150

-200

-250

1000
-100

-300
500 z
-200
0
-100
-350
Local minimum
y -200
-300
-400

-300
x

-400 -450

-400 -450 -400 -350 -300 -250 -200 -150 -100

Global minimum

Figure 1-1. Three-dimensional and contour plot of the Schwefel function.

Figure 1-1 shows the Schwefel function defined over a two-dimensional search space
xi = ( x1 , x 2 ) . Local and global optima are discriminated. A global one is an optimum of the

whole feasible set D while a local optimum is only an optimum of one of its subsets. The
definitions of local and global optima 20 are:
Local maximum definition: A local maximum xl∗ ∈ D of an objective function F(xi) is a

solution with F ( xl∗ ) ≥ F ( xi ) for all xi neighbouring xl∗ .

If D ⊆ ℜ , the definition can be written:

xl∗ : ∃ε > 0 : F ( xl∗ ) ≥ F ( xi ) ∀xi ∈ D xi − xl∗ < ε

Local minimum definition: A local minimum xl∗ ∈ D of an objective function F(xi) is a

solution with F ( xl∗ ) ≤ F ( xi ) for all xi neighbouring xl∗ .

If D ⊆ ℜ , the definition can be written:

-5-
xl∗ : ∃ε > 0 : F ( xl∗ ) ≤ F ( xi ) ∀xi ∈ D xi − xl∗ < ε

Global maximum definition: A global maximum xl∗ ∈ D of an objective function F(xi) is a

solution with F ( xl∗ ) ≥ F ( xi ) ∀xi ∈ D .

Global minimum definition: A global minimum xl∗ ∈ D of an objective function F(xi) is a

solution with F ( xl∗ ) ≤ F ( xi ) ∀xi ∈ D .

Depending on the context the term optimum can replace either one of the terms maximum or
minimum. Some optimization problems search for a global maximum value while others
search for a global minimum value. To find the global optimum and not just a local one, it is
necessary to search globally in the search space. But many nonlinear optimization problems
are non-convex and are likely to have many local optima. If traditional local search methods
are used to solve such a problem, the outcome depends on the starting point of the search, i.e.
as a rule only the next local optimal solution is found which can be of rather bad quality. In
order to find the global optimum, a global search strategy is needed. The most important
global optimization heuristic strategies are described below, together with short descriptions
and references. Most global optimization software implementations are based on one of these
approaches, possibly combining ideas from several strategies.

1.2 Heuristics
The word heuristic stems from Latin heuristicus or Greek heuriskein that means "to
find, to discover". In computer science, a heuristic is a technique designed to solve a problem.
The technique usually produces a good solution and helps to solve a simpler problem that
contains or is in close agreement with the solution of the more complex problem. Heuristics
increases computational performance or simplicity, at the expense of precision. There are a
great variety of definitions in the literature 21,22,23,24 . Heuristics is particularly used to find a
solution that is usually reasonably close to the best possible answer. More accurate methods,
especially when they are applied to complex problems, tend to show slow convergence that
goes along with a high computational cost. The main reason for this slow convergence is that
these methods explore the global search space by creating random movements without using
much local information about promising search directions or already visited areas. In contrast,
local search methods converge much faster as they use local information to determine the
most promising search direction. However, they are easily entrapped in local optima.
The basic concept of the heuristic search as a support in solving of problems was
introduced by Polya 25 . A classic example of a heuristics is hill climbing. Hill climbing is a
local search optimization technique. It can be applied to a large variety of optimization

-6-
problems. The algorithm starts with a random (potentially bad) solution and applies a local
search to find an improved solution. If such a solution is found, the search moves to it and the
local search starts again. The method stops when the solution is not improved upon already
obtained solutions after ordinary local search. Ideally, at that point a solution closed to the
optimal solution is found, but it is not guaranteed that hill climbing will ever come close to
the optimal solution. The main disadvantage of a hill-climbing method is its incapacity to
escape from a local optimum. A way to improve the performance of this simple heuristic
procedure considerably is metaheuristic.

1.2.1 Metaheuristic Features

In the beginning, heuristics was typically applied to the solution of complex
combinatorial optimization problems. The term “metaheuristics” for solving such tasks was
introduced by Glover 26 . It has been widely applied in the literature, both in comparative
studies as well as in research papers. Since the introduction of metaheuristics by Glover a lot
of books and monographs were published on this subject9,12,14,15,24,27,28,29,30,31,32 .
Metaheuristics includes all heuristic methods which allow finding a good quality solution and
also local optimization techniques. The picture has changed drastically. Modern powerful
computers and parallel platforms allow successful applications of metaheuristics to real-time
problems within acceptable time 33,34 . However, metaheuristics gives no guarantee of
obtaining the global solutions either.
The use of the metaheuristic search strategies makes it possible to escape from a local
optimum during the iterative procedure. Metaheuristics has been developed to solve complex
optimization problems in various fields. The efficiency of the procedures is due to the use of
information about the environment. The method is a result of adjusting metaheuristic
strategies to specific optimization problems. It is important to note that metaheuristic and
“local optimality” approaches have essential differences. Typical heuristic concepts are "rules
of thumb", educated guesses, intuitive opinions or simply common sense. They ignore
whether the solution to the problem can be proven to be correct. Such iterative heuristics
belong to local search methods such as descent or ascent methods. Therefore, the
development of methods which allow moving aside from this classical design and broadening
the application area is of great importance.
During the last few years, metaheuristics developed substantially. Now it is based on
“intelligent” search methods that can be classified with respect to the three basic ideas:
• Use of adaptive memory and responsive exploration
• Neighbourhood search

-7-
 • Number of current solutions carried from one state to the next.
These options help to classify metaheuristic methods in the form X/Y/Z, which was proposed
by Glover and Laguna 35 . The possible values of X, Y, Z elements are noted below:
A method employs adaptive memory
X
M method is “memoryless”

N method uses some systematic neighbourhood search

Y
S method relies on random sampling

I method with point-to-point strategy

Z
P population-based approach

This simple classification scheme gives us a possibility to categorize various metaheuristics

methods. Such a classification must stay ambiguous because of large quantity of method
modifications; some of them are closed to “standard concepts”, others are more or less
expanded developments of them. A metaheuristic method may be classified in two or more
X/Y/Z forms. Examples are shown Table 1.

Table 1. Metaheuristics classification.

Metaheuristic method Classification 1 Classification 2

Genetic algorithms M/S/P M/N/P
Simulated annealing M/S/I M/N/I
Tabu Search A/N/I A/N/P
Scatter search M/N/P A/N/P

Surely, metaheuristics can include other strategies to find the global optimum. A
metaheuristic method may modify the search strategy on the basis of the change of the
objective function value during the search procedure. The same modification may be carried
out during a neighbourhood search by excluding some members and introducing others. It can
be illustrated at the example of the strategic oscillation approach of the Tabu Search. The
scope of activity of the standard neighbourhood strategy that includes moves only among
feasible solutions is extended in order to involve unfeasible solutions as well. Accordingly,
the search overcomes the feasibility boundary in order to proceed into the unfeasible region.
Another classification of the metaheuristics, which differentiates between population-
based strategies and single-solution metaheuristics, is also often found in the literature. In the
latter methods, a search step requires only the information of a single preceding solution
which the next iteration starts. On the other hand, population-based strategies invoke a
collection of solutions at each stage. Such procedures are usually characterized by the class of
evolutionary methods. Well-known examples are genetic algorithms, the scatter search, and
the path relinking method35 (a useful combination of intensification and diversification

-8-
strategies) based on strategies for “combining” solutions. It is also a conspicuous subclass of
the metaheuristic methods which use multiple heuristics to generate new population members
instead of relying on a single rule.
It is also possible to classify metaheuristic methods based on the use of adaptive
memory. Patterns, whose present state depends on the sequence of previously visited
solutions, and therefore includes a covert form using “memory” term, reside in practically all
heuristics except those that use complete randomization. All previous choices are
“remembered” and are inherited by the current one. The term “memory” has represented
primarily by Tabu Search and its variations are sometimes called the “adaptive memory
programming” 36 , but a number of other metaheuristics use mechanisms that can also be
considered as memories. In recent years, other approaches have made attempts to unify
various aspects of such memory structures and strategies, however, typically in only rather
simple form. Developments, which produce hybrids of the tabu search with other approaches
at a more advanced level, have become an important way to embed an adaptive memory into
other methods, and have established an active area of research. In genetic algorithms and the
scatter search the memory is employed to store a population of the solutions, where the mode
of combination more clearly lends itself to transmit features of selected past solutions to
current solutions. However, such an implicit memory is not an intelligent memory
construction. It lacks any mechanism to save the solutions prior to the last ones as well as
methods to compare and improve the current solutions with the preceding generations.
Neural network approaches introduce another memory based distinction. Artificial
neural networks have roots in our understanding of the human brain. Initial concepts were
based on attempts to mimic the brain's way of processing information. Subsequent efforts
gave rise to various models of the biological neural network structures and learning
algorithms. Such methods accentuate an associative form of memory. Neural network
approaches implicitly involve a form of optimization, and they have been applied to several
optimization problems in recent years. In spite of a variable success of the performance,
neural networks are often regarded as being appropriate to be included within the
metaheuristic classification. Many research problems have been solved using neural network
techniques (in particular the Hopfield network 37 ) including graph partitioning 38,39 ,
knapsack 40 , and constraint satisfaction problems 41 as well as linear and nonlinear 42
programming. Looi 43 mentioned that “although there is a large collection of operations
research based approach and other methods for solving all of these problems, comparisons
between neural network methods with existing methods have been lacking". Some researchers

-9-
made up these deficiencies by successfully combining neural networks with Simulated
Annealing, Genetic Algorithms, and most recently Tabu Search.
Metaheuristics is often considered to be a set of intelligent components, but it is
important to note that the intelligence depends more on the underlying design than the
particular property (or behaviour) of the method itself. For this reason, it is not necessary for
a procedure to qualify as intelligent in a rigorous sense in order to grant its membership in the
category of metaheuristics. Figure 1-2 shows a graphical representation of various
metaheuristics whose principles are presented below.

Metaheuristic

Population
Scatter search

Genetic algorithms
Particle swarm
optimization
Ant colony
optimization
Simulated
Annealing
Memoryless

Trajectory
Variable
Tabu search neighbourhood search

Figure 1-2. Metaheuristic methods.

1.2.2 Population metaheuristic

This class of metaheuristics uses populations of solutions. It provides better results
than single-solution metaheuristics. New solutions are obtained by an implicit or explicit
combination of different solutions.

1.2.2.1 Swarm methods

There are two popular swarm inspired methods in the area of computational
intelligence: Ant Colony Optimization (ACO) and Particle Swarm Optimization (PSO). ACO

-10-
was inspired by the behaviour of ants and has many successful applications in discrete
optimization problems. The particle swarm concept originates from a simulation of simplified
social systems and can be used as an optimizer.

1.2.2.1.1 Ant Colony Optimization (ACO)

The ACO algorithm introduced by Marco Dorigo in his PhD thesis 44 and originally
developed by him, Maniezzo and Colorni in the early 1990s 45,46 , is a probabilistic technique
for solving combinatorial optimization problems that can be realized by finding good paths
through graphs. The behaviour of ants 47 in finding paths from the colony to food sources gave
rise to this idea. Even if a single ant has only restricted abilities, the behaviour of a whole ant
colony is highly structured as the consequence of coordinated interactions. The way of ant’s
communication is a chemical compound, known as pheromone. A moving ant lays various
quantities of these compounds on the ground, thus it marks the path its moves with a
pheromone trace. Usually a single ant begins to move randomly, but by detecting a
pheromone trail, the ant will follow it with higher probability intensifying it with its own
pheromone. Hence the probability that ants will follow a path correlates with the number of
ants that did so before. This is a form of autocatalytic behaviour — or allelomimesis.
Pheromones also vanish with time if they are not refreshed. If all of the food is taken away
from a particular place, ants will stop putting pheromones onto the respective track since they
cannot find any food at this location anymore. The process is thus characterized by a positive
feedback loop.
In Figure 1-3, an example is presented that this mechanism can lead to the shortest
paths (as the pheromone tends to accumulate faster on these paths). It starts with a given path
from the ant hill to food source. If this path is cut off by an obstacle, the ants have to pass the
obstacle along the right or left path. Each ant makes a choice on the basis of some heuristic
evaluations and the intensity of the pheromone trails left by previous ants. The path with a
good heuristic evaluation and a high level of pheromone is more likely to be selected, so it
gives the following ant a stronger stimulus and thus a higher probability to turn right. The first
ant forced to decide which path must to be taken has the same possibility to turn right or left
(since there is no pheromone trail on the alternative paths). If the right path is shorter than left
one, the ant that took it will be faster than the ant following the left (long) path. Due to the
shorter distance they can move more often. The next ants will find a stronger trail on the right
path and it will become preferred (in probability) to the left path. The probability with which
an ant decides to move along the path to follow is more and more prejudiced towards the
shorter one because the intensity of pheromones increases faster on the shorter path. The final

-11-
result is that all ants will quick-witted in choosing of the shorter right path. However, it is
significant that the decision is never deterministic, thus there remains always a possibility to
explore alternative routes.
Computational models have been developed, to simulate this mechanism48,49 : the
problems are visualized as (directed) graphs containing various junctions. In a preparation
step, a few ant individuals perform randomized walks through the graphs. Pheromones are
laid out in accordance with solution profit, so the probabilities to walk along given paths
increase. In the following steps, the ants move again through the graph. They cannot only
move along the already visited paths, but also can choose other routes, because randomization
is used to allow the construction of a variety of different solutions. Again pheromones are laid
out in accordance with solution profit. When all ants have constructed a complete solution,
the procedure is restarted with the updated pheromone level. This is repeated for the number
of allowed iterations or until the solution does not improve after a number of iterations. Thus,
a minimal level of individual complexity can explain a complex collective behaviour. With a
more complex prescription for the single steps it is possible to change the algorithm such that
it escapes from local optima and copes better with environmental changes. However, a lowest
limit is required to establish the desired behaviour.
ACO is now developed to a powerful, many-sided optimization tool with a lot of
publications and numerous applications in diverse areas of operations research, management
science and technology. A number of refinements have been integrated into this general
iterative scheme. The approach of Dorigo and Di Caro29 is an alternative to the Tabu Search,
Genetic and Evolutionary Algorithms, Simulated Annealing, the Iterated Local Search, the
Variable Neighbourhood Search, and some other search methods. Dorigo and Di Caro
reviewed successful implementations of the ACO method which are applied to a number of
important and difficult combinatorial optimization problems, such as:
• Quadratic assignment32,50,51
• Travelling salesman 52,53,54,55,56,57
• Scheduling 58
• Connection-oriented network routing 59,60,61
• Connection-less network routing 62,63,64
• Vehicle routing29,32,65,66
• Sequential ordering 67
• Graph colouring 68
• Shortest common supersequence 69
• Conformational analysis of flexible drug-like molecules 70 .
-12-
 1 2 3
Figure 1-3. Allelomimesis example.
1) Some ants are walking on a path from the ant hill to food source
2) An obstacle suddenly appears and the ants must get around it
3) At steady-state the ants choose the shorter path

A large amount of experimental work on ACO was carried out in the last years. In their recent
survey Dorigo and Blum used the ACO theory 71 . Meanwhile, the results of different ACO
implementations have converged if applied to standard optimization problems 72,73,74,75,76 , and
also first investigations modelling the dynamics and the finite-time dynamics of ant colony
optimization are carried out 77,78 . Gutjahr reported an article 79 that addresses the important
question how the (expected) runtime needed to obtain a solution of a given quality scales with
the size of the problem. The first results with increased complexity of the iteration steps were
given for two ACO algorithms. One of them is GBAS with lower pheromone bounds
(GBAS/lb) which also forms the backbone of the Simulation based Ant Colony Optimization
(S-ACO) algorithm 80,81 for stochastic combinatorial optimization, which has already proven
successful in applications 82,83 .
Besides ant activity ACO can include pheromone trail evaporation and daemon
procedures. These techniques are useful to avoid a too strong convergence of the algorithm
towards a sub-optimal region. Pheromone trail evaporation decreases automatically the
pheromone intensity as a function of time. Useful form of forgetting is implemented in order
to favour exploration of new areas of the search space. Individual daemon actions are used to

-13-
centralize actions, which can be performed by single ants. Examples of such actions are local
procedures or the collections of global information that can be used to analyse the usefulness
of laying additional pheromone to bias the search process from non-local perspective. The
daemon can monitor the path which was found by each ant and add extra pheromone to the
shortest path. Such pheromone updates are called offline pheromone updates.
Korošec and Šilc 84 introduced that even real vector optimizations can be treated as
graph problems. As shown by the simulation results, the approach can be extended to a
broader class of problems. The general idea of this model is that of an autocatalytic process
pressed by a "greedy force". The greedy force alone is incapable to find anything but a
suboptimal tour. The autocatalytic process alone leads to convergence to a suboptimal path
with exponential speed. In their combination greedy force steers the autocatalytic processes
towards the best available local optimum and lets it converge quickly to very good, often
optimal solutions.

1.2.2.1.2 Particle Swarm Optimization (PSO)

PSO is a stochastic, population-based evolutionary optimization algorithm. It was
developed in 1995 by Kennedy and Eberhart 85,86 . It uses the so called swarm intelligence that
is based on social-psychological principles and comprehension biological social system. By
analogy with a swarm, each swarm member (particle) takes advantage of own memory and
has a degree of freedom, or randomness, in its movement as well as knowledge gained by the
whole swarm to find the best available food source 87,88,89,90 . After the first algorithm
introduction, it has been studied by a number of different authors 91,92,93,94 who mostly
concentrated on multimodal mathematics.
A problem of a food search can be solved by optimizing a fitness function. The
definition of the communication structure (social network) is realized by assigning the
neighbours for each swarm. Each single solution is a "swarm member" called particle in the
search space. All particles have fitness values which are evaluated by the fitness function to
be optimized, and have velocities which direct their motion in the multidimensional search
space. Any single particle iteratively evaluates the fitness of the solutions it meets and keeps
track of the coordinates in the search space associated with the best solution (fitness) it has
achieved so far as well as the fitness values. Each particle remembers the information about
its best solution, its position in the search space and both are available to its neighbours. Thus,
there is a local best solution given by the best solution of all particles in the neighbourhood.
When all particles in the swarm are considered as neighbours, this best solution is the
optimum presently available.

-14-
 The swarm is typically simulated by particles in multidimensional space that have a
~
position x(p) ∈ X ⊆ ℜ n and a velocity v(p) ⊆ ℜ n . The position x(p) of a particle p indicates a
possible solution for the problem whereas its velocity v(p) determines the direction of the
subsequent search. Both the position and the velocity are real vectors. This ensured that PSO
is especially suitable for a numerical optimization. The particles move through the solution
vector hyperspace (i.e. ℜ ) and know:
n

• their own best solution

• the presently best available solution,
where "best" simply means the solution with the optimal fitness function value. In the
beginning, the positions and velocities of all individuals are randomly initialized. The swarm
members swap for information about good positions with neighbours and update first their
own velocity and then the position based on this information. In order to carry out
appropriate changes of its position and velocity, each particle p has a memory holding the
following information:
• best(p) – “Particle best” position which is the best solution the particle has seen
itself
• best(N(p)) – “Local (neighbourhood) best” position that the particle obtains by
communicating with a subset of the swarm
• best(Xpop) – The optimal (“global best”) position ever obtained by any
individual in the population that is known to all and immediately updated
when a new best position is found by any particle in the swarm.
In the original PSO two different kinds of neighbourhoods were defined. To adjust the
velocity of the particle p the PSO algorithm use either the best(N(p)) or the best(Xpop)
positions.
• In the “global best” swarm, all the particles are neighbours of each other; thus, the
position of the best overall particle in the swarm is used in the social term of the
velocity update equation. Exploitation only the “global best” position leads to rapid
convergence of the algorithm as all the particles are attracted simultaneously to the
best part of the search space but decreases the probability to find the global optimum
when the global optimum is not close to the “particle best”. It may even become
impossible for the swarm to explore other areas. It means that the swarm is trapped in
local optima.

-15-
 • In the “local best” swarm, only a certain number of particles can affect the velocity of
a given particle. With recourse of neighbourhood information the global optimum is
found more likely for the prize of slower convergence.
The ith particle velocity update equations for the cases of global and local PSO are given by

v(p)i = v( p ) i + (ci ∗ randomn () ∗ (best ( p) i − x( p) i ) +

(2)
(d i ∗ randomn () ∗ (best ( X pop ) i − x( p ) i ))

v(p)i = v( p) i + (ci ∗ randomn () ∗ (best ( p) i − x( p) i ) +

(d i ∗ randomn () ∗ (best ( N ( p)) i − x( p) i )) (3)
In both cases the ith particle positions are updated by
x(p)i = x( p) i + v( p) i (4)

The learning rate vectors c and d have further influence of the convergence speed and the
ability of the swarm to find the optimum. Moreover, the values of all dimensions of x(p) are
normally kept within the bounds of the search space. If the inertial coefficient of the velocity
is small, all particles could slow down until they approach zero velocity at the “global best”.
The fitness of the “global best” solution improves with each swarm iteration. It could also
happen that all particles, which are to be influenced by the “global best” swarm, move in
close proximity to the “global best” in the search space without exploring the rest of search
space. In this case the fitness never improves in spite of number of made PSO iterations. The
way to avoid this situation is to reinitialize the particle positions at a certain interval or after
the detection convergence. Algorithm illustrated the native form of the Particle Swarm
Optimization can easily be generalized for multi-objective optimization and for returning sets
of optimal solutions.
Both “global best” and “local best” can be seen as "social" neighbourhoods, as the
relations among particles does not depend on their positions in the search space, but on
"external" relationships that are not dependent on the problem that is being solved. Kennedy
and Mendes 95,96 investigated the alternative "social" neighbourhood topologies. The
following additional neighbourhood topologies were tested:
• Random
• “von Neumann”, a two dimensional grid with neighbours to the north, east, west, and
south
• Pyramid, a three-dimensional triangular grid
• Star, all the particles connected to a central particle

-16-
 • Heterogeneous, particles are grouped in several cliques
• Hypercube
• Ring
Illustrations of some of the topologies are presented on Figure 1-4.

a) Global b) Ring (local) c) Random d) Star

e) „von Neumann“ g) Pyramid

f) „von Neumann“
(two-dimensional)
(three-dimensional)

Figure 1-4. Neighbourhood topologies.

There are two factors offered by Watts 97,98 , which can be used to characterize the different
neighbourhoods:
• The degree of connectivity that measures the number of neighbours of a particle
• The amount of clustering that measures the number of neighbours of a particle that are
also neighbours of each other
Despite the fact that the results depend on the selected variable, the “von Neumann” and
Pyramid neighbourhoods prove to be the best while the star and “global best” are the worst.
The concept of a "dynamic" neighbourhood was also explored 99 . In the work of Li 100 , a
dynamic neighbourhood topology was used for a multimodal function optimization.
Recently, PSO has been successfully applied in many areas. Few applications of the
algorithm to structural and multidisciplinary optimization are known. Fourie and Groenwold
suggested exploitation of particle swarm in size and shape optimization 101 and an application
to topology optimization 102 . Some example areas of application of particle swarm
optimization are:
• training of artificial neural networks86,103
• training of hidden Markov models 104
-17-
 • global optimization of mathematical functions86,105
• in antenna or filter design 106
• water resource and quality management 107
• quantitative structure-activity relationship (QSAR) modelling in chemistry 108,109
The PSO algorithm was successfully applied to a continuous and integer/discrete
structural optimization problem. It is inherently a continuous algorithm although it requires
much higher computational cost than gradient-based optimization methods. The results show
that the PSO algorithm is more suited for integer/discrete and discontinuous problems where
use of a gradient-based optimizer may not be appropriate. PSO has many similarities with
evolutionary computation techniques such as Genetic Algorithms (GAs). Both algorithms
start a population of random solutions and search for optima by updating generations and
using random techniques. They don’t give guaranteed success. However, unlike GA, PSO has
no genetic operators like crossover and mutation. The advantages are that PSO is easy to
implement and possessing only a few fitting parameters.

1.2.2.2 Genetic Algorithms (GAs)

The idea of applying computer-aided simulations of biological principle of natural
evolution to study artificial systems was first introduced in the 1950s by biologists like
Barricelli. Some pioneering works, which traced back to the 1960s made by Bremermann 110
and Bledsoe 111,112 , use evolutionary approaches based on binary genomes to function
optimization and balance weights for neural networks. Later on, Bagley 113 , Cavicchio 114 , and
Frantz 115 made a further important research on these genomes that crowns with the main
presentation of the GAs of Holland in 1975 116 . Works related to the subject are from De
Jong 117 , Baker 118 , Goldberg27, Mühlenbein and Schlierkamp-Voosen 119 , Chipperfield et
al. 120 , Reeves24, Michalewicz30, and Berthiau and Siarry 121 .
GAs, evolution strategies, evolutionary programming, and genetic programming are
grouped under the term evolutionary algorithms or evolutionary computation. They have been
successfully applied to numerous problems from different domains, including optimization,
automatic programming, machine learning, economics, ecology, and population genetics,
studies of evolution and learning, and social systems. Although GA is a subclass of
evolutionary algorithms, it was restrictedly applied until multipurpose presentation of
Goldberg27 in search, optimization, design, and machine learning areas. Nowadays, GAs are
considered to be the most widely known and applicable type of metaheuristics30,122,123,124 with
many applications of genetic algorithms in science, economy, research, and development 125 .

-18-
 A GA is a search technique to find exact or approximate solutions of various
combinatorial optimization problems. GAs use techniques inspired by natural processes based
on the Darwinian principle of natural selection such as “mutation”, “selection”, and
“crossover” (also called recombination). Through the “selection” process, only the best
solutions are allowed to become “parents” and to generate “offspring”. This probabilistically
biases the algorithm towards the best elements in the population. “Crossover” is the mating
process. With a given probability, it takes two selected individuals, called “parents” and
combines their most desirable features by exchanging parts of their genomes solutions to
create one or two new individuals, called “offspring”. In the simplest form, substrings are
exchanged after a randomly selected “crossover” point. This operator tends to enable the
evolutionary process to move toward “promising” regions of the search space. “Mutation” is
performed for a few offspring: for such offspring, one variable is altered by a small
perturbation, for instance the change of one bit in the binary coding case. It is introduced to
prevent premature convergence to local optima by randomly sampling new points in the
search space. “Mutation” entails flipping bits at random, with some small probability.
GAs are implemented as computer simulation of an optimization problem with a
population of abstract representations, called chromosomes, genotypes or genomes, and of
candidate solutions, called individuals, creatures, or phenotypes. For the definition of a
typical genetic algorithm, one needs the followings:
1. genetic representation of the solution space
2. fitness function to evaluate the solution space.
Bit arrays are a standard solution representation, but also arrays of other types enabling the
crossover operation can be used. A fitness function is a particular type of objective function
defined over the genetic representation. It quantifies the optimality of a solution so that a
particular solution may be ranked against all the other ones. The fitness function is always
problem dependent. Optimal solutions or at least more promising solutions are allowed to
breed and to mix for producing a new generation that will hopefully be even better.
The standard GA starts with an initial population that is randomly generated. Every
evolutionary step is called a generation. The predefined quality criterion, the fitness or fitness
function is then evaluated for each individual. To create a new population (the next
generation), individuals are selected according to their fitness. Many selection procedures are
currently in use, one of the simplest being Holland's original fitness-proportionate selection,
where individuals are selected proportional to their relative fitness. This ensures that the
expected number of times an individual is chosen is approximately proportional to its relative
performance in the population. Thus, high-fitness individuals possess a better chance of

-19-
“reproducing”, while low-fitness ones are more likely to disappear. Parents are combined, by
means of the recombination operator, to produce offspring. Before replacing the old
population, the members of the new population receive the small random perturbations by
means of the mutation operator. The offsprings mute with a given probability, and the fitness
and the objective function value of the resulting offspring are computed. Then the next
generation is created. The cycle is iterated until some optimization criteria are reached or
when a maximum number of generations has been produced. Figure 1-5 sketches out the
structure of a simple genetic algorithm.

Random initialization of the

population of individuals

Evaluation of the
fitness function for
each individual

Is optimization Yes Final result

criterion arrived? (best individual)

No

Formation of new population

Selection

Crossover

Mutation

Figure 1-5. Simple genetic algorithm.

There are a large number of different types of genetic algorithms. The Continuous
Genetic Algorithm (CGA) 126 , which will be of some importance later, was proposed by
Chelouah and Siarry as an adaptation of GA to the global optimization problems of
continuous multiminima functions. It uses a type of real coding, which is as close as possible
to Holland’s approach using binary coding. The main contribution of the CGA is the

-20-
introduction of two concepts also widely used in Tabu Search: diversification and
intensification. In the diversification phase, CGA starts with a large population, and a high
mutation probability, to cover the whole search space homogeneously, and detect a promising
area. The intensification phase is then performed inside a promising area, after having
reduced the search domain, the population size and the mutation probability. The “selection”,
the “crossover”, and the “mutation” are performed using the decimal code. There are also
combinations between GAs and derivative based methods 127 like the quasi-Newton method to
solve difficult unconstrained optimization problems.
Genetic algorithms find application in bioinformatics, phylogenetics, computer
science, engineering, economics, chemistry, manufacturing, mathematics, physics, and other
fields. Some example areas of application of genetic algorithms are:
• scheduling applications, including job-shop scheduling 128,129,130,131
• chemistry and chemical manufacturing 132,133,134
• medicine 135,136,137,138
• data mining and data analysis 139,140,141,142,143
• geometry 144,145,146,147,148
• finance and trade 149
• optimizing distributed protocols 150,151
• building phylogenetic trees 152
• chemical kinetics (gas and solid phases)
• design of water distribution systems
• distributed computer network topologies
• game theory equilibrium resolution
• gene expression profiling analysis 153
• linguistic analysis 154
• marketing mix analysis
• mobile communications infrastructure optimization
• molecular structure optimization 155,156
• multiple sequence alignment 157
• operon prediction 158
• protein folding and protein/ligand docking 159
• RNA structure prediction 160
• timetabling problems 161
• training artificial neural networks

-21-
 • finding hardware bugs. 162,163
Genetic algorithms are a very effective way of quickly finding a reasonable solution to
a complex problem. They don’t give instantaneous and exact effects, but do an excellent job
of searching through a large and complex search space. GAs are most effective in the search
space with little knowledge about a problem to be solved, but can also produce solutions that
only work within the test environment.

1.2.2.3 Scatter Search (SS)

Scatter search 164,165,166,167 is a population-based metaheuristic that constructs new
solutions by intelligently combining and improving previous solutions, called reference set.
The update is made with regard to the results of the improvements. The methodology is very
flexible, since each of the scatter search elements can be implemented in a variety of ways
and degrees of sophistication.
SS algorithms apply “crossover” or combination mechanisms that allow the sharing of
information between solutions to create new solutions 168 similar to the GAs. The main idea is
a combination of the characteristics of two parent vectors to create several offspring, with the
assumption that good solutions may generate better ones. The main search operator of the
GA 169,170 and SS 171 approaches is handling of available information to influence future
searches. Unlike the large population in GAs, the reference set of solutions in SS is rather
small. In SS, the selection of the parents is not random as in traditional GAs, but the two or
more elements of the reference set are chosen in systematic way. Due to this selection process
that considers at least all pairs of solutions in the reference set, the reference set can be kept
small. The methods applied in the SS to improve the solutions range from simple local
searches to very specialized searches. For example, simple local searches are based on the
selection of the most improving move or of the first improving move during the selection.
Tabu Search35, a Variable Neighbourhood Search 172 or various sophisticated hybrid heuristic
searches15 are also applied as variable neighbourhoods, intermediate memory, or hashing
scanning methods of the neighbourhood. From this standpoint, SS procedures can be
considered as memetic algorithms (MAs)29,173 , which combine local search heuristics with
crossover operators.
From the implementation standpoint the SS method can be divided into the following
subroutines based on the well-known templates165,174 :
• The Diversification Generation Method is used to generate a collection of diverse trial
solutions from one or more arbitrary seed solutions as an input.

-22-
 • The Improvement Method transforms a trial solution into one or more enhanced
solutions. Neither the input nor the output solutions are required to be feasible, though
the output solutions will usually be expected to be so. If no improvement occurs for
the given trial solution, the resulting solution is considered to be the same as the one
submitted for improvement.
• The Reference Set Update Method creates and maintains a set of reference solutions
consisting of the best solutions found according to the criteria of providing efficient
accessing by other parts of the method. The goal is to ensure diversity while keeping
high-quality solutions.
• The Subset Generation Method generates subsets of the reference set as a basis for
creating combined solutions.
• The Solution Combination Method uses weighted structured combinations to
transform each subset of solutions produced by the Subset Generation Method into
one or more combined solutions
The reference set is a collection of high quality solutions and diverse solutions. Both
are required by the Solution Combination Method. Its goal is to produce weighted centres of
selected subregions and to project these centres into regions of the solution space. This space
shall be explored by auxiliary heuristic procedures. The reference set is generated from
diverse solutions that improves during the search and will provide the information for the
search process. The reference set update is based on comparisons between new and already
visited solutions. To build the large set of diverse solutions the Diversification Generation
Method is used. The size of the set of diverse solutions is typically smaller than the size of the
reference set. The initial reference set is built according to the Reference Set Update Method.
The following simple mechanism demonstrates the initialization of the reference set
and its updating procedure during the optimization search. The initialization of the reference
set starts with the selection of the best solutions from the set of diverse solutions built by
Diversification Generation Method. These best solutions are added to the reference set and
deleted from the set of diverse solutions. For each improved solution of the updated diverse
set the minimum of the Euclidean distances to the solutions of the reference set is computed
and the solution with the maximum of these distances is added to the reference set. The
process is repeated until the set of the diverse solution becomes empty. Solutions in the
reference set are ordered according to their quality. The best solution is the first one in the list.
The simple Subset Generation Method creates all pairs of reference set solutions and
puts them in a list to apply the Solution Combination Method. This method tries to combine
good characteristics of the selected solutions intelligently to get new high quality solutions

-23-
after the improvement. At each iteration the number of subsets is generated depends on the
number of new added solutions of the reference set. These trial solutions are subjected to the
Improvement Method. The Improvement Method is the range of approaches from the simplest
local searches to very specialized procedures. The reference set is updated according to the
quality and the dispersion of the improved found solutions. The process is iterated with the
new reference set until a stop condition is met. Finally, the set of disperse and high quality
solutions in the reference set is provided by the method. A general template for the SS
algorithm can be organized in two phases outlined as follows on Figure 1-6.

Initial Phase
No
Diversification Generation

Is it feasible
Improvement solution?

Reference Set Update

Yes

No Scatter Search Phase

Subset Generation

Is optimization
criteria arrived? Solution Combination

Yes
Improvement

Final result
Reference Set Update

Figure 1-6. Scatter Search algorithm.

Recent studies confirm the practical advantages of this approach for solving a various
optimization problems. SS differs from other evolutionary procedures, such as GAs, by
providing combined principles for the solutions based on generalised path constructions in the
Euclidean space. Furthermore, the search for the local optimum in SS is guided with less

-24-
space of randomization. Additional advantages of this approach are provided by
intensification and diversification mechanisms. It uses adaptive memory elements that link SS
to Tabu Search. The success of SS and related strategies such as continuous version of SS 175 ,
which works directly with vectors of real components, is applied in a variety of application
areas, such as:
• dense wavelength division multiplexing 176
• graph colouring 177
• data mining 178
• linear ordering problem29
• network design problem 179
• history-matching problem 180
• chemical and bio-process optimization 181
• 3D image registration problem 182
• neural network training 183
• multi-objective routing problem 184
• commercial implementation 185
• classical vehicle routing 186
SS is an evolutionary metaheuristic usually presented as a non-nature inspired one. However,
most of the implementation that can be found in the literature sources have nature-inspired
elements 187 in its components or are very similar to those used in standard nature-inspired
metaheuristics.

1.2.3 Simulated Annealing (SA)

Simulated annealing (SA) is a generic probabilistic technique for locating good
approximations to the global optimum in a large search space. It is based on the principles of
thermodynamics. Name and idea of the method come from the annealing process in
metallurgy, a technique involving heating, maintaining a suitable temperature and controlled
cooling of a material. The goal is to induce softness, to relieve internal stresses, and to refine
the structure so that the working properties of the material are improved. Annealing induces
the diffusion of atoms within a solid material, so that the material progresses towards its
equilibrium state. Heat increases the rate of diffusion by providing the energy needed to break
bonds so that the atoms can move from their initial positions (a local minimum of the internal
energy) and can wander randomly through states of higher energy. The slow cooling opens
the chances to arrive at some more stable (global or local minimal potential energy)
equilibrium. This general principle is applicable to both discrete and continuous global

-25-
optimization problems under mild structural requirements14,28,188 . It forms the basis of an
optimization technique for combinatorial and other problems.
The original ideas of SA come from a paper published by Metropolis et al. in 1953189
that introduced an efficient algorithm to simulate the equilibrium of a collection of atoms at a
given temperature. The algorithm simulates the cooling process by gradually lowering the
temperature of the system until it converges to a steady, “frozen” state. The adaptation of the
Metropolis-Hastings algorithm in optimization as Simulated Annealing was first proposed by
Kirkpatrick, Gelatt and Vecchi (1983) 190 , and by Cerny (1985) 191 .
By analogy with this physical process, in a combinatorial optimization context a
solution corresponds to a state of the physical system and the solution value corresponds to
the energy of the system. At each iteration the current solution is replaced by a randomly
selected trial solution. The probability to take the new or the old solution depends on the
difference between the corresponding function values and on a main control parameter in the
cooling schedule called temperature T. This means that the new solution is accepted according
to the so-called Metropolis criterion. The difference between the function value of the
resulting solution F ( posi +1 ) and the function value of the current solution F ( posi ) is defined
as:
ΔF = F ( pos i +1 ) − F ( posi ) (5)

If ∆F is negative (i.e. the function value of the resulting solution is better than the function
value of the current solution), resulting solution is accepted and becomes the new solution in
the chain. If this difference is positive (i.e. the old solution is better than the new one), the
resulting solution is only accepted on the basis of a comparison of some probability BF (Eq.
6) with a random generated number (RGN) between 0 and 1. BF is related to the magnitude of
the cost increase and a parameter T. If BF ≥ RGN, this new conformation is accepted;
otherwise, it is rejected.

⎛ ΔF ⎞
BF = exp⎜ − ⎟ (6)
⎝ kT ⎠
where ∆F is the increase in F and T is the main control parameter.
The main control parameter T is progressively lowered during the cooling process
according to the given schedule and a certain number of iterations is performed at each
temperature level. In general, a move is the most probable to be accepted if the temperature is
high and the cost increase is low. The dependency is such that the current solution changes
almost randomly when T is large. If T goes to zero only improving moves are accepted,

-26-
“downhill”, and the method stops at a local optimum, The allowance for "uphill" moves,
which avoids to become trapped in local minima, is the major advantage over other methods,
such as Tabu Search, Simulated Annealing. The illustration of the main principle of SA is
represented in Figure 1-7. The new solution with better function value is always accepted, but
the new state with worse function value is only accepted with a certain probability. The
probability of accepting a solution is high at the high temperature and decreases at the drop in
temperature.

F(x)

New solution accepted

with probability

Start solution
Local minimum

New solution
always accepted

Global minimum

a b x

Figure 1-7. Selection of new states in simulated annealing.

The implementation of the basic SA algorithm is straightforward. The following

Figure 1-8 shows its structure. It has been proved that with proper cooling schedule, SA
asymptotically converges to a global optimum. However, this requires infinite number of
iterations. Finite-time implementations do not provide such a guarantee.

Since the presentation of Kirkpatrick et al., a lot of studies that exploit SA have
appeared. The theory as well as the applications of SA has been widely studied12,15,192,193 . The
subsequent developments were carried out by Aarts, Korst, and van Laarhoven28; Aarts and
Ten Eikelder9; Henderson, Jacobson, and Johnson15 and have been focused on the following
topics:

• convergence based on more general forms of acceptance rule than the Metropolis one
• deterministic variants 194,195,196
• different forms of static and dynamic annealing schedules 197,198

-27-
 • parallel annealing 199,200,201,202
• thermostatistical persistency 203
• hybridizations with other metaheuristics 204,205,206,207 .

Initialization

Generate new solution

Assess new solution

Accept new No
solution?

Yes

Store as a best solution

Reduce T

No Terminate
search?

Yes

Result

Figure 1-8. The flowchart of the SA algorithm.

SA becomes an algorithm of choice in such application areas as:

• financial instruments and trading 208,209
• complex portfolio selection problem 210
• econometric problems 211,212

-28-
 • job shop scheduling problem (JSS) 213
• travelling salesman problem (TSP) 214,215,216
•
vehicle routing problems196
• global minimization of the Lennard-Jones energy function 217,218
• molecular conformation problem for peptide 219
• determination of silicon clusters 220
• prediction of crystal structures 221,222
• design of the composite structures for aircraft 223
• plastic design of circular steel plates 224
• water distribution in irrigation canals 225
• path planning216,226,227
• paper cutting waste optimization 228
• seismic waveform inversion 229 and etc.
The ability of easily escaping of stagnations in local minima by accepting "uphill"
moves through a probabilistic procedure especially in the earlier stages of the search is one of
the most powerful features of SA. It is a robust and general technique and its main advantages
are its flexibility and its ability to approach global optimality. It is able to handle more classes
of constraints, highly nonlinear models, and chaotic or noisy data. This makes it, as neural
nets and GAs, highly promising for portfolio optimization and asset allocation problems.
However, it appears less appropriate for modelling financial time series. Due to the fact that
the algorithm does not rely on any restrictive properties of the model, it is quite versatile. The
possibility of “tuning” SA methods is also very significant, especially for any reasonably
difficult nonlinear or stochastic system. This optimization algorithm can be tuned to enhance
the performance. The ability to tune SA for more than one problem should be considered as
an important feature of an algorithm.
A disadvantage is that the SA methods are computation-intensive. Since SA is
metaheuristic, rather delicate work is needed to account for different classes of constraints and
to fine-tune the parameters of the algorithm to adjust them to the actual problem. There is a
clear tradeoff between the quality of the solutions and the time required to compute them. The
precision of the numbers used in the implementation of SA have also a strong influence on the
result quality. The main disadvantages that have been noticed on SA are its general slow
convergence and its wandering around the optimal solution in the case where high accuracy is
needed.

-29-
 For some problems211,212, SA algorithm performs better, e.g. more reliably in finding
more optima, than other numerical techniques such as GAs. The technique appears well suited
to asset allocation problems with constraints. It is often used when the search space is
discrete. SA does not use gradient information and makes relatively few assumptions about
the problem being solved.

1.2.4 Variable Neighbourhood Search (VNS)

VNS is relatively recent metaheuristic designed for solving optimization problems. It
was proposed by Mladenović and Hansen9,15,230 and is based on the principle to explore
neighbourhoods of growing size iteratively.
VNS is based on the following three rules:
1. The global optimum solution is a local optimum in all possible neighbourhood
structures.
2. The local optimum solution with respect to a single neighbourhood structure does not
imply that it is a local optimum in any other neighbourhood structure.
3. For many problems local optima of one or more neighbourhood structures are
relatively close to each other.
The last item is an empirical observation. It suggests that a local optimum often provides
information about the global one. Therefore an intensive study of the neighbourhood of this
local optimum enables the finding of a better.
To describe the VNS algorithm a set of neighbourhood structures is denoted by Ni
with i = 1,…, imax and the set of solutions in the ith neighbourhood of the solution x is denoted
by Ni(x) (see Figure 1-9).

Nimax
…
N2
N1
x

Figure 1-9. Set of neighbourhood structures.

-30-
Neighbourhoods Nk may be determined by one or more metric (or quasimetric) functions that
define a distance between the elements of a set in the solution space. In order to solve general
optimization problem (Eq. 1) using several neighbourhoods, facts 1 to 3 can be implemented
in the following ways:
a) deterministic approach
b) stochastic approach
c) combination of deterministic and stochastic approaches.
(a) Variable Neighbourhood Descent (VND) is a popular deterministic approach
where the best neighbour of the current solution automatically becomes the new current
solution if an improvement is obtained. The search is then restarted from this point. If the
current solution is better than all its neighbours the next neighbourhood is considered. The
search stops when all neighbourhoods are visited and no improvement is possible which
means that the solution is a local optimum for all neighbourhood structures. The steps of the
resulting basic VND scheme are given below:
Initialization:
Select the set of neighbourhood structures Ni , for i =1, . . . , imax, that shall be used
in the descent;
Find an initial solution x;
WHILE (no improvement is obtained){
Set i=1;
WHILE (i<imax ){
Exploration of neighbourhood: Find the best neighbour x' of x ( x' ∈ N i' (x) );
Moving: IF (F(x') is better than F(x)){
Set x = x' and i = 1;
}ELSE Set i = i + 1;
}
}
Most local search heuristics use in their descents strategies imax ≤ 2 neighbourhoods. Since the
final solution should be a optimum with respect to all imax neighbourhoods, the opportunities
to reach a global optimum are larger than by using a single neighbourhood. A VND heuristic
can easily be constructed for a given problem by combining it with available heuristics from
the literature. The order of its applying and its ranking by non-decreasing size of the
neighbourhoods has a significant influence. In addition to the sequential order of
neighbourhood structures in VND above, one can go further and split neighbourhoods used in

-31-
a heuristic into a set of nested ones of increasing size (nested strategy). Such an approach is
applied in Ref. 231, 232, 233.
To increase the efficiency of this scheme, one should consider the complexity of the
different moves, the order of their application, the sufficiency of the considered moves to
ensure a thorough exploration of the region containing x, and desired solution precision. As a
consequence, if selecting and ranking moves involve a lot of elementary changes, the
resulting heuristic may be very slow and often take more time than the exact algorithm on
small or medium size problems. A frequent implementation involves ranking moves in the
order of the complexity of their application and a return to the easiest type each time a new
descent is found and explored. If some neighbourhoods are much easier to explore than
others, the algorithm returns to the first neighbourhood as soon as an improved local optimum
is found, instead of repeating all steps in sequence. The sufficiency of the considered moves is
a crucial question: for some problems elementary moves are not sufficient to leave a narrow
valley. Hence such heuristics can achieve poor results since they become trapped in local
minima. In the beginning, one will strive to obtain the best possible solution within the
allocated computing time by the exploitation of only the VND method. Afterwards, it is
preferred to get good solutions fairly quickly by the VND and to improve it later by the faster
stochastic search in VNS.
(b) Local optima of many combinatorial and global optimization problems tend to be
close with each other and situated in one or sometimes several parts of D. Reached local
optimum contains implicit information in its neighbourhood about better ones and even
perhaps global optimum. Hence it is reasonable to explore the vicinity of the discovered
optima.
The Reduced Variable Neighbourhood Search (RVNS) 234 is the method of choice for
these aims. The random solutions are selected from Ni(x) without being followed by the
search for the best neighbour. The neighbourhoods Ni(x) with i = 1,…, imax are nested in most
cases, i.e. each neighbour contains information about the previous one. It was observed that
imax = 2 is the best value. A solution is randomly chosen in the neighbourhood. If this solution
is better than the previous one (i.e., if F(x') < F(x) for minimization), the search is transferred
to the new position (x = x'). Otherwise, the search is proceeded in the next neighbourhood.
After visiting of all neighbourhoods the search continues again in the first one, until a
stopping condition is met. Usually, a maximum number of iterations or a maximal computing
time between two improvements is set as stopping criteria. Owing to the nested form the size
of successive neighbourhoods increases. Therefore, close neighbourhoods of x are explored

-32-
more thoroughly than farther ones, even if no further improvements were observed within the
first. The steps of RVNS are following:
Initialization:
Select the set of neighbourhood structures Ni, for i =1, . . . , imax, that will be used in
the search;
Find an initial solution x;
Choose a stopping condition;
WHILE (stopping condition is not met){
Set i = 1;
WHILE (i<imax){
Shaking: Generate a random point x' from the ith neighbourhood of x
( x ′ ∈ N i (x) );
Moving: IF (F(x’) is better than F(x)){
Move to x' (x=x'):
Continue the search with N1 (i =1):
} ELSE Set i =i + 1.
}
}
RVNS is useful for cases where the local search is expensive. It is akin to a Monte-Carlo
method, but more systematic.
(c) In the previous two ways, directions, which use variable neighbourhoods to reach
a local optimum and to find promising regions for near-optimal solutions, were examined.
Combination of the tools for both tasks leads to the General Variable Neighbourhood Search
(GVNS) scheme.
Merging the local search with systematic changes of the neighbourhoods around the
current local optimum leads to the Basic VNS method230. It is based on the idea to change the
neighbourhood structure systematically within a local search heuristics, rather than staying in
a single neighbourhood structure. In this basic scheme, a set of neighbourhood structures are
selected. They are often nested, i.e. they define the neighbourhoods around the current
solution x ∈ D of the solution space. A solution x' is randomly generated in the first
neighbourhood of the current solution. A local search is performed yielding the local optimum
x" of this neighbourhood. If x" is worse or equal to the best previous solution, the procedure is
iterated using the next neighbourhood. The search is recentered around x" and restarted using
this neighbourhood as the first one. If no better solution has been found or if every
neighbourhood structure has been explored the search begins again at the first neighbourhood

-33-
N1(x) until a stopping condition is met. The stopping criteria may be the allowed maximum
CPU time, maximum number of iterations, or maximum number of iterations between two
improvements.

N1(x)

x
x' x"

F(x")< F(x) F(x") ≥ F( x)

N1(x)

x N1(x)
x' x"
x
Ni(x")
N2(x)

Figure 1-10. Variable neighbourhood search for minimization.

The illustration of the basic VNS is presented in Figure 1-10. The steps of this
approach are following:
Initialization:
Select the set of neighbourhood structures Ni, for i =1, . . . , imax, that will be used in
the search;
Find an initial solution x;
Choose a stopping condition;
WHILE (stopping condition is met){
Set i=1;
WHILE ( i< imax) {
Shaking: Generate a random point x' from the ith neighbourhood of x
( x ′ ∈ N i (x) );
Local search: Apply some local search method with x' as initial solution;
denote with x" the so obtained local optimum;
Moving: IF (this local optimum is better than the incumbent){
Move there (x=x");
Continue the search with N1 (i=1);

-34-
 }ELSE Set i=i+1
}
}

The GVNS scheme is obtained if the VND is used for the local search of x" and if the
initial solution is improved by RVNS. This scheme has the following steps:
Initialization:
Select the set of neighbourhood structures Ni, with i = 1, . . . ,imax, that will be used in
the shaking phase, and the set of neighbourhood structures Nl with l = 1, . . . , lmax that
will be used in the local search;
Find an initial solution x and improve it by using RVNS;
Choose a stopping condition;
WHILE (stopping condition is not met){
Set i=1;
WHILE ( i< imax){
Shaking: Generate a point x' at random from the ith neighbourhood Ni(x) of x;
Local search by VND:
Set l=1;
WHILE (l =lmax ){
Exploration of neighbourhood: Find the best neighbour
x" of x' in Nl(x');
Moving: IF (F(x") better than F(x')){
Set x'= x" and l=1;
}ELSE Set l=l+1;
}
Moving: IF (this local optimum is better than the incumbent){
Move there( x=x");
Continue the search with N1 (i=1);
}ELSE Set i=i+1;
}
}
The VNS heuristic is able to find the best valleys rather quickly. To avoid trapping in
a valley, the set of the neighbourhoods around all feasible solutions x should contain the
whole feasible set D:
D ⊆ N i ( x) ∪ N 2 ( x) ∪ ... ∪ N i max ( x), ∀x ∈ D .

-35-
These neighbourhood sets may cover D without its partitioning. It is easy to implement, e.g.
by nested neighbourhoods (i.e., each one contains the previous):
N 1 (x) ⊂ N 2 (x) ⊂ ... ⊂ N imax (x), D ⊂ N imax (x), ∀x ∈ D .
It is also possible to explore D completely by investigating the small neighbourhoods around
solution along some path, but such approach is not guaranteed for more efficient search.
Nested neighbourhoods are easily obtained for many combinatorial problems. For this
purpose it is necessary to define a first neighbourhood N1(x) by a type of move and then to
iterate it i times to get neighbourhoods Ni(x) for i = 2, . . . , imax. The sizes of the subsequent
neighbourhoods increase. Since the investigation of the whole search space one goes many
times through the whole sequence of neighbourhoods. In this case the first neighbourhoods
will be explored more thoroughly than the last ones. This meets the requirements of third rule,
i.e. local optima with respect to one or more neighbourhood structures are relatively close to
each other. Moves to the feasible set D may be constrained, particularly if this set is
disconnected. Introducing some or all constraints in the objective function with Lagrangian
multipliers allows moving to infeasible solutions 235 .
Three improvements of the basic VNS for solving large instances are now considered:
• the Variable Neighbourhood Decomposition Search (VNDS) method234 extends the
basic VNS into a two-level VNS scheme based upon decomposition of the problem;
• the Skewed Variable Neighbourhood Search (SVNS) method 236 addresses the problem
of exploring valleys far from the previously found best solution;
• the Parallel Variable Neighbourhood Search (PVNS) method15,237,238 is a way for
parallelizing VNS
As the change of the neighbourhood during the search for good solutions to (Eq. 1) is a simple
but very powerful tool, several authors have added such a feature to other metaheuristics
rather than VNS:
• Tabu search 239,240,241,242
• GRASP 243,244,245
• Constraint programming 246,247 .
VNS has been successfully applied in many areas; some of them are listed below:
• Vehicle routing problem242
• Arc routing problem 248
• Travelling salesman problem 249
• Linear ordering problem 250
• Graph coloring 251
• Protein side chain placement problem 252

-36-
 • Molecular Distance Geometry Problem 253
• Fuzzy clustering problem 254,255
• Finding the three-dimensional structure of a molecule 256
• p-Median problem 257
• Simple plant location problem 258
• One- and multi-dimensional knapsack problem 259,260
• Oil pipeline design problem 261
• Maximum clique problem 262
• Degree-constrained minimum spanning tree problem 263
• Max-cut problem 264
• Cable layout problem 265
• k-Cardinality tree problem 266,267
• Generalized minimum spanning trees problem 268
• Nurse rostering problems 269
• Resource-constrained scheduling problem 270
• Generalized minimum edge biconnected network problem 271,272
• Pooling problem 273 .
VNS metaheuristic is based on a simple and clear principle, which is easy to
understand, and most important, easy to use. This approach provides optimal or near-optimal
solutions for solving the problems of several benchmarks within reasonable computing times.
Furthermore, the performance of the VNS appears to be robust. It has a few parameters or
sometimes none. Using VNS one could obtain good or better results than most other
metaheuristics on many problems. In view of the aforesaid VNS fulfils the requirements of
desirable properties of metaheuristics and may play an important role in the design of new
heuristics for new types of applications.

1.2.5 Tabu Search (TS)

Let’s dwell up on a description of the TS method that forms the basis of our new
approaches and which developed to one of the most powerful methods for solving
combinatorial optimization problems.
TS is a metaheuristic algorithm that belongs to the class of local search techniques
with the possibility to accept worse-cost local solution in order to escape from local optimum.
It was first proposed in its present form by Glover in 1986 274 . Seminal related ideas were also
developed by Hansen 275 which he called steepest ascent/mildest descent method. It is
significant that in Glovers view TS was a metaheuristic, i.e., a general strategy for guiding
-37-
and controlling interior heuristics specifically adjusted to the problems. Additional efforts of
formalization and further developments for combinatorial optimization problems were
reported later35,276,277,278 . A didactic presentation of TS and some of its applications have been
collected in a book of Glover, Taillard, Laguna, and de Werra 279 .
The basic principle of TS is to guide a local heuristic search procedure for exploration
of the solution space beyond local optimality. TS works on a single solution candidate. The
local search procedure uses an operation called move to create a neighbouring individual of
this candidate. The definition of the moves is highly problem-specific. The combination of the
current solution with the possible moves defines the neighbourhood of the current solution. Its
size simpliciter depends on the number and the kinds of defined moves. Consequently, if
every neighbour has to be evaluated at each iteration, the algorithm may become quite slow.
Therefore, usually only a subset of the neighbourhood is considered. Random selection of a
fixed number of neighbours for consideration is a commonly employed method.

It seems to be familiar.

Figure 1-11. Adaptive memory and responsive exploration.

The main feature of TS is its use of an adaptive memory, which creates a more flexible
search behaviour, and responsive exploration. Mountain climbing can be used as an
illustrative example of both features. An adaptive memory term can be explained by the
example of an alpinist who selectively remembers the path elements during the climb of the
mountain in memory. Using this knowledge to make a strategic choice along the way
represents the responsive exploration. It is illustrated in Figure 1-11. The adaptive memory
feature of TS allows the implementation of procedures that are capable of searching the
solution space more economically and effectively. Using the mountaineer analogy, one must

-38-
analyze current alternatives using previously obtained attainments during ascents in similar
regions. The responsive exploration in TS, whether in a deterministic or probabilistic
implementation, rests on the assumption that a bad strategic choice yields more information
than a good random choice. It integrates the basic principle of intelligent search i.e., using
features of previous good solutions to explore new promising regions.
Simple TS combines a local search procedure with anti-cycling memory-based rules to
prevent the search from getting trapped in local minima, so that a global optimization can be
conducted. It restricts the return to recently visited solutions by constructing a list of them
called Tabu List (TL). The last visited solution is placed on a TL, and the reverse move is
forbidden during a certain number of iterations. Certainly, after a local optimum is found, the
next neighbour is reached by an uphill move. In most cases, during the next iterations, the
search will try to go downhill again, but this reverse of an already accepted move is tabu, so
that the algorithm has to continue to go uphill along the modest ascent direction, i.e. it will
leave the local optimum. An important parameter is the TL size, which determinates the
number of moves that are forbidden, i.e. how long the visited solution will remain prohibited.
TL is usually managed with respect to FIFO (First In First Out) principle. The longer TL is,
the smaller the chances that the algorithm gets back to already visited solutions are, i.e. the
smaller the chance that it is trapped in local minima is. But with the increasing size of TL the
search gets more and more limited since feasible solutions could be missed because the moves
leading to them are tabu for a too long time.
Specific properties of TL or attributes of the tabu solutions can be made more effective
for some application domains. However, it can also put new problems and can lead to the
complication of the algorithms, because it may prohibit a new promising solution even when
there is no danger of cycling or may lead to an overall stagnation of the searching process.
Therefore, to avoid such problems aspiration criteria can be defined which revoke TL, so that
otherwise excluded solutions are included in the allowed set. The simplest and most
commonly used aspiration criterion includes solutions, if their objective values are better than
the current best-known solution. More complicated and therefore rarely used aspiration
criteria are the neglect of all tabus if cycling cannot occur. So, the update of the memory-
based TL can be modified and controlled by the following concepts:
• Tabu tenure: the length of time during which a certain move is classified as tabu. It
can be kept constant or varied dynamically throughout the search.
• Aspiration criteria: to accept improving solution the existing tabu rules can be
overridden.

-39-
In each iteration TS generates the neighbourhood of the current solution. After that the
selection of the subset is carried out by the remove of the tabu solutions from the
neighbourhood. Each solution in this new neighbourhood subset is evaluated with regards to
the objective function. Certainly, a generation process that avoids generating a trial solution
that is already recently visited is preferred. The best solution in the neighbourhood of the
current solution is selected as the new current solution and becomes the “initial” solution for
the next iteration. This will also include uphill movements which are necessary to escape
local minimum and to avoid getting trapped in an optimum. Such uphill movement go along
the modest ascent direction. In Figure 1-12 an algorithm scheme of a TS algorithm is defined.

-40-
 Initialization
Initial solution x0,
Best solution x*,
Tabu List

Neighbourhood generation

Evaluation of the
neighbourhood solutions

Selection of the best

neighbour x/
Delete x/

Is x/ the member Yes

of the Tabu List?

Update No
Tabu List

Is x/ better than No
x*?

Yes

Store as a best solution

x * = x/

x = x/

No Required criterion
satisfied?

Yes

Result solution x*

Figure 1-12. Flowchart of Tabu Search.

-41-
 Two very important components of TS are intensification and diversification
strategies. The basic purpose of the diversification search is to explore unvisited regions of
the search space that have not been investigated previously, whereas intensification search
intensifies the search within a limited space, e.g. the promising region. This main difference
between intensification and diversification is illustrated in Figure 1-13.

-100 Diversification solutions

Intensification solutions

-200

New solutions of
unvisited regions
-300

-400

-500

-500 -400 -300 -200 -100

Exploration of the neighbours of elite solutions

Figure 1-13. Intensification and diversification.

Typically, TS stops after a fixed number of iterations or a maximum number of

consecutive iterations without any improvement of the incumbent solution. During these
iterations TS keeps sufficient information about the search in memory. Such a simple TS
structure is called short-term memory TS. The short-term memory used in the search
intensification holds only the recent history of all moves or subset of them. Its basic idea is to
determine the attributes of performed moves. This information is used to encourage move
combinations and solution features, which were found good in previous stages of the search.
It is also possible to return to previous promising regions to search them more thoroughly.
The intensification search focuses on investigation of the neighbourhood of the best solutions
found during the search which are called elite solutions. In intensification stage the term
neighbour has a broader meaning as in the usual context of neighbourhood search. During the
intensification search neighbours are also produced by the combination of the components of

-42-
good solutions or by modified evaluation strategies that favour the introduction of such
components into a current solution.
Long-term memory has been suggested to achieve a better performance and to keep a
more important search information such as the recency, the quality, and the frequency9.
Specifically, long-term memory in high-level TS saves attributes of special characters like
elite and frequently visited solutions. Accordingly, the second highly important component of
the TS, the diversification strategy, may be based on some long term memory, such as
frequency-based memory. This frequency-based memory keeps the total number of iterations
that various “solution components” have been presented in the current solution or have been
involved in the selected moves. By means of this algorithmic mechanism, one makes the
attempt of forcing the search into previously unexplored areas of the search space and
generates solutions that differ in various significant ways from those seen before.
There are some major diversification strategies. The simplest form of diversification is
certainly the classical random restart diversification, in which the search is periodically
stopped (according to some predetermined criterion) and then restarted from a randomly
generated solution. More interesting but also more complicated approaches select restart
solutions by using historical information collected in a long-term memory. Examples may be
rarely used components in the current solution or the incumbent solution. Such information
can be used in different ways, for instance, to avoid highly frequent moves or, conversely, to
encourage the use of moves having very low frequencies and finally to design restart
mechanisms within TS. This principle can also be exploited in a continuous diversification
which integrates diversification considerations directly into the regular searching process.
This avoids interruptions and the loss of information accompanied with the restart scheme.
This is achieved by the modified evaluation of the possible moves which added a small term
related to the component frequencies 280 to the objective function. In this way, neighbourhood
search can be directed into yet unexplored regions where TL operation is restarted.
The TS ideas of intensification and diversification are beginning to find their way into
other metaheuristics. Since the appearance of the simple TS scheme described above, many
new developments and refinements have been proposed over the last few years. They are
briefly mentioned in the followings:
• intensification and diversification techniques: In addition to previously described methods
the following forms of long-term memories35,280 are often implemented:
– adaptive memory procedure referred to as ‘‘probabilistic diversification and
intensification”, which was first proposed by Rochat and Taillard 281 for the vehicle
routing problem. Adaptive memory involves an attribute-based focus which depends

-43-
 intimately on the elements of recency, frequency, influence, and logic. Its main idea
is that attributes of high quality solutions are used to construct other high quality
solutions, as it was done in GAs. The method therefore stores a continually updated
pool of previously generated elite solutions, which can be extracted from the pool
and used to restart the search. The part of the pool that represents better solutions
has a larger probability of being selected. Usually, different fragments of elite
solutions are taken and combined to generate a new starting solution. If the
components are taken from elite solutions, which are situated in a common area of
the search space, this is the case of intensification, otherwise it is a question of
diversification. Adaptive memories provide a generic paradigm for guiding local
searches and can be coupled with different types of metaheuristics.
– strategic oscillation or tabu tunneling35 is the method to achieve an effective
cooperation between intensification and diversification. It manages moves in
relation to a critical level called the oscillation boundary which often represents a
point where the manner would normally stop. The boundary is identified at various
stage of the construction or it represents a chosen interval of functional values. The
search is allowed to go for a specified depth beyond the boundary before turning
around. When the boundary is crossed again from the opposite direction, the search
goes beyond it for a predefined depth and after that turning around again. The
process of repeatedly approaching and crossing the critical level from different
directions creates an oscillatory behaviour. It is possible to control this behaviour by
varying the amplitude of the oscillation, by generating modified evaluations and
rules of movement to explore a region of the search space depending on the type of
the region and the search direction.
– path relinking35,282 generates new solutions by exploring connecting trajectories of
elite solutions. It starts from one of these solutions, called initiating solution, and
generates a path in the neighbourhood space leading to another solution, called the
guiding solution. This can be done by selecting moves that introduce attributes
contained in the guiding solution. This technique can be used in both diversification
and intensification strategies, depending on the path generation mechanism and the
choice of the initiating and guiding solutions.
• reactive TS 283 : (Battiti and Tecchiolli) It provides a technique which dynamically fits the
search parameters based on the search history. The continuous reactive TS is a
generalization of the reactive TS developed by the same authors. The continuous reactive

-44-
 TS method tries to locate most promising boxes. In a second step, refined solutions are
generated within these boxes.
• hybrid TS: It was developed by Al-Sultan and Al-Fawzan 284 . It represents a hybrid
method that combines TS with the local search method of Hooke and Jeeves 285 .
• probabilistic TS35: It provides a selection scheme of the moves which is based on elite
solution recovery. It also introduces randomization into the search method. This can be
attained by converting the evaluation and the references to tabu status into probabilities of
selection, strongly made over to favour higher evaluations. Such randomization can also
be implemented in the design of a candidate list strategy. It only considers a subset of
neighbourhood solutions, because a complete neighbourhood evaluation is often too
expensive.
• directed TS 286 : It uses direct-search-based strategies. These strategies are based on the
well-known Nelder-Mead method 287 and a new pattern search procedure called adaptive
pattern search. The role of direct search methods is to make the search robust especially in
the neighbourhood of a local optimum. Instead of using completely blind random searches
for the generation of the neighbourhood, appropriate direct search strategies are exploited.
Moreover, a TL conception with new anticycling rules called tabu regions (TR) and semi-
TR were introduced.
• unified TS 288 : It aims to produce simpler and more flexible TS codes. The major benefits
of the approach are its speed, its simplicity, and flexibility coupled with a dynamic
adjustment of only a few parameters. A similar approach is the granular TS 289 , which is
based on the use of drastically restricted neighbourhoods. They exclude moves with
elements that are not likely to belong to good feasible solutions. Such restricted
neighbourhoods are called granular, and may be seen as an efficient implementation of
candidate-list strategies proposed for TS algorithms.
• TS for continuous multiminima problems 290 : Cvijovic and Klinowski 291,292 generalized
TS to functions with continuous variables but this approach also allows the optimization
of integer values. Franzé and Speciale 293 introduced an algorithm which explores a grid of
points with a dynamically defined distance. Chelouah and Siarry 294 presented a new
algorithm called enhanced continuous TS. It results from an adaptation of combinatorial
TS, which aims to follow Glover’s basic approach as close as possible. To cover a wide
domain of possible solutions, this algorithm firstly performs the diversification. It locates
the most promising areas, by fitting the size of the neighbourhood structure to the
objective function and its definition domain. When the most promising areas are located,

-45-
 the algorithm continues the search by intensification within one promising area of the
solution space.
There are some examples of TS application areas:
• general manufacturing problems 295
• travelling salesman problem (TSP) 296,297
• quadratic assignment problem 298
• general combinatorial problems24,299
• optimization of artificial neural networks 300
• resonance assignment in NMR (Nuclear Magnetic Resonance) spectroscopy 301
• test design in education 302
• vehicle routing problems 303
• chemical process optimization 304
• identification of the good variable subsets within the framework of QSAR
techniques 305 .
The above examples have shown that TS is a powerful algorithmic approach. It has
been applied with great success to many difficult problems and has now become an
established optimization technique which can compete with almost all known techniques. It
can even exceed many classical procedures due to its flexibility. A great advantage of TS is
that, as in the case of all approaches based on a Local Search, it can quite easily deal with
really complicated constraints that are typical for real-life problems. Local search algorithms
move from solution to solution in the search space until a solution considered as optimal is
found or a time bound is elapsed. This ensured that it is a really expedient, but not a universal
approach. Significant problem knowledge is absolutely required to perform the most basic
steps of the development of any TS procedure, viz. the choice of the search space of an
effective neighbourhood structure and the moves between the solutions. The fine tuning of an
apparently large collection of parameters is also a problem. Consequently, this choice must be
sufficient. To be successful, all metaheuristics need to explore the search space “in depth” and
“in width”. Depth is usually not a problem for TS, which commonly achieves sufficiently
acceptable good solutions in an early stage of the search. In this respect TS is quite aggressive
and effective, but width can be a hard problem. This disadvantage can only be removed by an
effective diversification scheme.

1.3 Conclusions
The optimization aim is to find a discrete mathematical object that maximizes or
minimizes an objective function specified by the user of the metaheuristic. Such discrete

-46-
mathematical objects called states and the search space are problem-specific. User usually
provides the function to be optimized, called objective function, as a procedure that evaluates
the function on a given state.
Many computer scientists point out that no metaheuristic is better than another when
its performance is averaged over all possible discrete functions (No-free-lunch theorem 306 ).
Therefore, at best, a specific metaheuristic can be highly efficient only for classes of objective
functions. However, all these imperfections are characteristic also for other methods of the
global optimization as well.
Metaheuristics are best and generally used for non-differentiable objective function
where one has no possibility to exploit analytical tools as Hessians and derivatives, or if the
objective function has no analytically closed form (for instance, if it is only the output data of
another algorithm process). But even for all other cases metaheuristics achieve very fast and
very accurate solutions if techniques from classical optimization such as local search and
gradients are implemented in the metaheuristic approaches.
Although many criticisms exist, metaheuristics are powerful algorithmic approaches
which have been applied with great success to many difficult combinatorial optimization
problems. As previously observed, many well-known metaheuristics seem to converge
towards a combined structure. New opportunities for combining the merits of these methods
appear in this unifying process leading to even more powerful search approaches for difficult
combinatorial optimization problems.
Based on an analysis of foregoing metaheuristic algorithms which evaluated their
strengths, their weaknesses, and their appropriateness for the solving of set of tasks the TS
method was chosen to create improved approaches. It forces also the fact that in contrast to
most of all metaheuristic optimization techniques TS are extremely effective for the standards
of continuous nonlinear optimization that confirms by last method developments. Detailed
descriptions of the new approaches are introduced in the next chapters.

-47-
Chapter 2 TS based methods

2.1 Gradient Tabu Search

In this work a modification of the TS method was developed, which for clarity will be
called Gradient Tabu Search (GTS). Its main goal is to enhance the efficiency of the TS by
adapting strategies developed to find minimum energy paths of the potential energy surfaces.
It exploits the analytical first and second derivatives. The GTS uses gradients for a fast
location of the next local minimum while second derivatives are used in the mildest ascent
part. For an efficient implementation novel concepts had to be introduced for storage of
previous moves. The GTS can be classified as a multi-start method. In contrast to least square
fit procedures, which are widely used for such problems, the GTS allows to find solutions
which do not represent the next local minimum.

2.1.1 Algorithm description

Three basic search procedures are used in the GTS method: the Local, the Modest
Ascent, and the Diversification searches. In the first search procedure one of the Local Search
methods or combinations of them are used to find the nearest local minimum. After
determining a local minimum, the Modest Ascent Strategy is employed to escape from this
local minimum and to move to the next one along the shortest way. Moreover, Tabu List,
Tabu Region, and Tabu Direction restriction rules are applied to avoid return to recently
visited solutions or to get caught in a local optimum. The Diversification Search (DS) is
needed in order to jump to other areas of the solution space that may be more promising for
the seeking global minimum.
Figure 2-1 shows the flowchart of the GTS. After the initialization one starts with a
search for improved minima. This part consists of local searches for the next minimum using
gradients and mildest ascent searches to escape to the next valley. Weighted Hessians are
used for the mildest ascent search. In between an update of the solution vector and the TL is
performed. If the solution does not improve after a number of iterations or if all
neighbourhood solutions are already set tabu, the search for improved minima is aborted and
the DS is performed to obtain new starting points. A pseudo-code of the GTS algorithm is
given in Figure 2-2.

-48-
 Initialization

Local search for next

minimum

Solution Update

Yes Terminate
Result
search?

No

Yes Is this region

Diversification Search to completely
obtain new starting point investigated?

No

Modest Ascent Search to

escape local minima

Figure 2-1. Flowchart of the new Tabu Search methods.

-49-
Initialization Locate+Store
//Start from randomly chosen point //Find next local minimum
Locate+Store Steepest Descent and/or Quasi-Newton
/or Powell’s methods
Main Loop Add minimum solution to the Tabu List
FOR ( iter = 0; iter < ITER _ MAX ; iter + + ) { Testing lower and upper bounds
//Escape local minimum IF (variables out of range) {
Mildest Ascent Search Diversification search
Locate+Store } GOTO Locate+Store
}ELSE{ Store minimum in the Result List}
IF (solution does not improve
after a number of iterations){
Mildest Ascent Search Diversification search
//Select the general move direction GOTO Locate+Store }
+, if Fzi+ < Fzi− IF (we obtained the same best results){
Di = Diversification search
−, if Fzi+ > Fzi− GOTO Locate+Store }
WHILE( F ( xi ) > F ( xi −1 ) ){
//Move in the direction of the lowest Di Next X
Next X (xi+1) FOR ( i = 1; i ≤ NDIM ; i + + ){
IF (xi+1 not TABU ){ Fzi = F ( x1 , x2 ,..., x j + D j × Δx 0j , x j +1 ,..., x NDIM )
GOTO Next X
}ELSE { next non-tabu Mildest Ascent} //Compute new variable vector
} ⎛ Fz max − Fzi ⎞
ranki = rank min + (rank max − rank min )⎜⎜ ⎟⎟
⎝ Fz max − Fz min ⎠
xi +1 = xi + Δxi0 × ranki
}

Diversification search Depend on problem

Figure 2-2. Pseudo-code of the Gradient Tabu Search algorithm.

2.1.1.1 Modest Ascent and Local Search Strategy

In the GTS, the definition of the neighbourhood is crucial. For a point
r
x = ( x1 ...xi ...x N ) the neighbourhood is often defined by those points for which all xi are
increased by a given ∆xi. From this neighbourhood, the GTS either selects the move that leads
to the largest decrease in F(xi) or accepts the one that leads to the mildest increase in F(xi). A
reduction of the number of steps is obviously obtained if within one step all xi’s are
appropriately modified and if computations of the sizes of the ∆xi’s take into account the
shape of F(x1 . . . xi . . . xN). Both improvements can be achieved if gradients and/or second
derivatives are employed.

-50-
 y 1 – start solution

2 – local minimum

3 – neighborhood search solution

1
4 – next local minimum

2 3

Figure 2-3. Schematic search strategy illustration.

The rough search scheme is presented on Figure 2-3. Moves from the solution 1 to 2
and 3 to 4 are carried out by Local Search methods, which are listed above, whereas
displacement 2 to 3 is realized with the Modest Ascent strategy.
The Steepest Descent and the Quasi-Newton method are known to be very efficient for
local optimization problems. Details of the both methods are described in Ref. 22-25. Steepest
Descent is an algorithm for finding the nearest local minimum of a function which
presupposes that the gradient of the function can be computed. This method also called the
Gradient Descent method, starts at a point xo and, as many times as needed, moves from xi to
xi+1 by minimizing along the line extending from xi in the direction of − ∇F(x i ) , the local
downhill gradient. The problem with the Steepest Descent method is that it performs many
extremely small steps in going down a long, narrow valley. Quasi-Newton method uses the
Quasi-Newton Broyden Fletcher Goldfarb and Shanno (BFGS) approximation to built up the
Hessian updates based on past steps. There are many variants of Quasi-Newton methods. All
of them approximate the Hessian matrix from the function and gradient values of some or all
of the steps previously taken. In the present work the BFGS method was used, in which the
model is not stored explicitly, but is calculated by gradients and step directions stored from
the past steps. The Quasi-Newton method was chosen because it converges typically quite fast
and does not require computation of the Hessian matrix, which can be quite expensive both,
in terms of the symbolic computation and numeric evaluation. The disadvantage of the
method is that it may diverge away from the minimum

-51-
 To reinforce the advantages and minimize the disadvantages a combination of both
methods was implemented in the GTS to locate the next local minimum. The search starts
with the Steepest Descent approach. When the gradient value becomes small enough the
program switches to the Quasi-Newton method. Using gradients reduces the number of steps
necessary to locate the next local minimum. Another advantage is that many terms which
need to be computed can be reused for the calculation of the gradients and Hessians. This
accelerates the computation in comparison to approaches which solely depend on function
evaluations.

It is well known that leaving the valley of the already found local minimum or the
determination of a transition state is considerably more difficult than finding the next local
minimum 307 . In so called eigenvector-following techniques 308 eigenvalues and associated
eigenvectors of the Hessian are employed to solve this problem. This method evaluates the
mildest ascent exactly by computation and by diagonalization of the complete Hessian. The
necessary effort to compute the Hessian increases roughly with NDIM∗(NDIM+1)/2 (NDIM
being the dimensionality of the problem) and additionally a NDIM∗NDIM matrix has to be
diagonalized. Therefore, such vector-following techniques could become more demanding
than strategies which are based on less expensive function evaluations. In order to reduce the
costs our approach only uses the diagonal elements of the Hessian to determine the direction
of the mildest ascent.

Figure 2-4 illustrates the connection between the objective function, its gradient and
second derivative in a local optima. The value of a second derivative gives the convexity or
concavity of the function on the given interval. At a local minimum the first derivative is
equal to zero, while the second derivative is positive, since the function is convex at this
point. The diagonal part of the Hessian provides information along which coordinate the
function is less convex. To determine the direction the functional values of each coordinate
around the local minimum were computed at the beginning of the modest ascent strategy.

-52-
 Figure 2-4. Connection between objective function Sin(x), gradient and second derivative.

Figure 2-5 illustrates an example of the direction choice. The local minimum is at the solution
(3.14; 4.438). The partial second derivative of F with respect to y with value 0.5005 is less
than the partial second derivative of F with respect to x which equals 1.005. It means that the
function is less convex along the y-coordinate so the modest ascent should go either in the
positive or in the negative direction along this coordinate according to the calculated function
values.

Figure 2-5. Griewangk function. Modest ascent strategy illustration,

where local minimum -(x;y)=(3.14; 4.438),

∂2 f ∂2 f
≈ 1 . 0005 , ≈ 0 . 5005 , - negative move direction, - positive move direction.
∂x 2 ∂y 2
Thus, in order to escape from local minima via the modest ascent it is necessary to determine
a new solution along a direction of the smallest function value with respect to information
about partial second derivative. A schematic illustration of the idea used to determine the
-53-
direction of the modest ascent from the diagonal part of the Hessian is presented on Figure
2-6.

4
4

3
2

0
2

-2

1
-4

-4 -2 0 2 4 6

Figure 2-6. Griewangk function. Schematic illustration of moving from one local minimum to another,
where 1 – start solution, 2 – local minimum, 3- neighbourhood solution, 4 – next local minimum
and Ø - schematic movement direction.

The values for the single ∆xi’s are determined according to the following formula
xinew = xiold + Δxi0 ∗ rank i
⎛ ∂2F −∂ F 2
2
⎞
⎜ ∂x max
2
∂xi ⎟ (7)
rank i = rank min + (rank max − rank min )⎜ 2 ⎟
⎜⎜ ∂ F 2 − ∂ F 2
2
⎟⎟
⎝ ∂x max ∂x min ⎠
The value for Δxi0 is user defined and can be tuned for the problem under

consideration. The actual step size is determined from Δxi0 and ranki, which is computed
from rankmin and rankmax. Variable with smallest partial second derivative value has the
largest step size while step sizes for all others variables are calculated with respect to partial
second derivative value ranks. The second derivative values are indexed before ranking
giving the minimal element the index 0 and the maximal one the index (NDIM-1), where
NDIM is the dimension of the problem. A linear ranking procedure gives the minimal list
element the maximum ranked value rankmax and the maximal list element the minimum
ranked value rankmin. The parameters rankmin and rankmax could also be used as fine tuning
parameters but considering their definitions the most meaningful value for rankmax is 1. Using
rankmin = 0, the direction possessing the smallest second derivative would not contribute to

-54-
the overall direction of the next step. Therefore rankmin = 0.1 is used. ∂ F and ∂ F
2 2

∂x 2
max ∂xmin
2

represent the highest and lowest second derivative at the point under consideration.
If the differences between the second derivatives become very small, this ranking
scheme leads to values too small for the higher second derivatives. Therefore, in cases with
∂2F −∂ F
2
< NDIM rankmin is redefined as
∂x 2
max ∂x min
2

⎛ ∂2F −∂ F 2
2
⎞
⎜ ∂x max
2
∂x min ⎟
rank min = rank max −⎜ ⎟ (8)
⎜⎜ NDIM ⎟⎟
⎝ ⎠

This modest ascent strategy is followed until a positive second derivative indicates that the
barrier to the next valley is crossed. From this point, the next local minimum is located using
the minimization strategy explained earlier.
This approach lies in between the eigenvector-following techniques and the so-called
low mode search (LMOD) 309 . The latter is often used for conformational analysis. In this
approach, eigenvalues and eigenvectors of the Hessian are computed at the minimum. To find
a new minimum, the LMOD-procedure follows the eigenvector of one of the low-lying
frequency modes until the increase in functional value exceeds a user-defined threshold.
During this following, the Hessian is not updated. A minimization starting from this point
leads to new minima in most cases. The gradient-only algorithm 310 , which was developed to
find saddle points, was also tested. This procedure uses the gradients to follow the weakest
ascent along a valley to the next saddle point. The procedure was successful in some test
cases but turned out to be extremely dependent on the chosen parameters. Therefore, it was
not adopted.

2.1.1.2 Diversification Strategy

The last element of the GTS represents the diversification search (DS), which is used
to select new starting points in the solution space. A DS becomes necessary if the solution
does not improve after a number of iterations that means that neighbours at hand are worse
than the current solution or if all neighbourhood solutions are already set tabu. One way to
enhance the convergence is to ensure that the search switches to regions that were not already
investigated. To explore the whole search domain, and localize promising areas several
strategies can be chosen.
The standard algorithm starts from the N initial solutions (N being the number of
parameters) which are randomly generated and uniformly distributed in the solution space. To
-55-
avoid returning to already visited areas and to force the search towards solutions sufficiently
far away from previous ones all solutions are excluded which already belong to Tabu Regions
(TRs). The objective function to be minimized is evaluated for each accepted solution, and the
best of them (with the lowest F(xi) value) becomes the new start solution, even if it is worse
than the previous solution.
In a second variant 3N solutions are randomly generated. 2N solutions are
concentrated in two user-defined subspaces of the whole search domain (e.g. positive and
negative areas if possible) and additional N points are randomly generated without any
restrictions. The solutions are excluded which already belong to TRs as in the previous
standard scheme. If a new promising solution with the smallest objective function value is
accepted the exploration process starts inside this newly detected area of the promising
solution. Depending on the problem Simulated Annealing methods may also be successful
variants of the diversification search. However, for the present problems they turned out to be
less efficient than the previously mentioned approaches.

2.1.2 TS Memory Elements

The basic idea of the approach is to keep sufficient information of the search within
some memory. In the Tabu List (TL) either complete moves or single attributes are set tabu
as soon as their complements have been part of selected moves. The moves or attributes stay
tabu for a distinct time, i.e. for a number of iterations, until the probability that a solution is
revisited is small. The efficiency of the algorithm depends on the choice of the tabu status
duration. Using an effective memory conception GTS behaves as an intelligent search
technique35. Optimization search methods are divided into two categories of optimization
search methods: point-to-point methods and population-based methods. GTS belongs to the
first category.
Despite its efficiency at times, the TL seems to be rather limited approach. The
standard TL can’t avoid the zigzagging move situation near local optima during the descent
method since the visited solutions in the TL don’t show exact tabu direction. This is no longer
sufficient since moves into already visited regions become allowed. To avoid such problems
for the new approach new techniques as: Tabu Region (TR) and Tabu Direction (TD) were
developed. These techniques enables a more efficient memory use then the standard TL, i.e.
with equal memory volume more previous solutions can be set tabu and not only single
visited solutions are set tabu but also areas around them and a complete visited route.
Reverse moves and cycles are avoided by the use of a Tabu List (TL), where the
moves previously done are memorized. If a new current solution can be found, it is added to

-56-
the end of the list until it is full. Thereafter, the added solutions must replace an existing
solution in the TL. The TL is managed with the First In First Out (FIFO) principle. The
elements in TL are ranked in ascending order according to their recency. The most ancient
element in the TL has index equal to 1, the most recent element has index equal to N, where N
is the number of TL elements. If the TL is full, the newest solution replaces the tabu solution
which has a smallest recency index that satisfies the requirements of the FIFO principle.
Experiments with Multi-Ranked Tabu List286 and with modification of Multi-Ranked Tabu
List, where also frequency ranking was added, were also carried out, but they offered no
particular advantage over the traditional FIFO updating method. For approaches which move
in small steps through the function F(xi) a simple storing of all visited points is sufficient to
block already visited regions effectively. In the present approach this is no longer sufficient
since the number of steps decreases while the step sizes increase. The specific property of the
TL that is implemented in the GTS approach is that not all previously visited solutions are
kept in memory but only the local minima solutions and start solutions from which they were
found.
To overcome the resulting problems Tabu Regions (TR) and Tabu Directions (TD) are
introduced.

4
4

3
2

2
0

-2
1

-4

-4 -2 0 2 4 6

1, 2, 3 and 4 – Tabu List solutions

- Tabu Region with radius RTR

Figure 2-7. Tabu List solutions and Tabu Regions on the example of surface contour of the Griewangk function.

The conception of TD is closely connected with the Tabu Direction Vector (TDV) and New
Move Direction Vector (NMDV) terms. To block already visited regions more effectively, the
visited minima a min are stored in the TL together with the starting points of the gradient

-57-
optimization from which the respective minimum was found ( bstart ). To ensure that the

search for the transition state starting from the current minimum a min does not move back, not
only the already visited points but the whole direction is set tabu. The TDV and the NMDV
are used for this purpose. The first vector, the TDV, is direction from the already visited local
minimum solution a min to the start solution bstart from which this local minimum was reached.

The second vector, the NMDV, is the vector a min c new , which connects a min with the next trial

point c new . Figure 2-8 illustrates TDV a min bstart and NMDV a min c new on the surface of the
Rastrigin function.

1 – start solution
2 – local minimum
3 – new solution
Ø - direction vector

Figure 2-8. Surface of the Rastrigin function. Tabu Direction Vector 21 from the local minimum to the start
solution and New Move Direction Vector 23 from the local minimum to the new solution.

To decide of the move to c new both vectors (TDV and NMDV) are computed and normalized.
If for the scalar product of both normalized vectors the equation

(1 − α ) < a min bstart ⋅a min cnew with ( α > 0) (9)

holds the new direction is tabu. Depending on α, which is given as input for fine tuning, the
region which is set tabu varies. Visualization is given in Figure 2-9.

-58-
 z
Tabu Direction

1 bstart
i
c 1i

c 2i

i
amin 1 y

1
x
Figure 2-9. 3D geometrical interpretation of TD. For α=0.29 the shaded part is tabu.

For α=0 it is only forbidden to move exactly in the direction of a min bstart . With increasing α, a

cone is becoming tabu. In Figure 2-9, α was set equal to 0.29 . Here the solution c1 is tabu,

since the vector a min c1 lies in the cone while solution c 2 is allowed.

Since only a few points are visited on the way from bstart to a min it is also necessary to
define Tabu Regions (TR) around those points for the diversification strategy. Figure 2-7
provides some examples of TL solutions and TRs. For this instance, a radius RTR is defined
for each point in the TL. The radii RTR are computed according to
( R up − R low )
R TR = (10)
coeff

Rup and Rlow define the range in which a parameter is allowed to vary. coeff is a user
defined variable which accounts for the range in which F(xi) is defined and the density of
local minima. The actual numbers employed for the test cases are given in Appendix A
together with the test cases. Please note that the TR’s are only needed for the diversification
strategy.

2.1.3 Implementation and Experiments

The efficiency of the present implementation was investigated with the help of various
test functions286,291,292,311,312,313,314,315,316,317,318 . The tests comprise the Ackley (AKn), the
Branin (BR), the Goldstein-Price (GP), the Griewangk (Gn), the Rastrigin (Rn), the Hansen
(H), and the Levy (Ln) functions. The corresponding expressions are given in appendix A,
which also reports the known global minima and the search ranges for each of the functions.
A characterization of the computational efforts, which were necessary to determine the global

-59-
minima, can be seen in Table 2. As mentioned earlier, the GTS algorithm in its present
implementation contains parameters. For a solid assessment of the success of the new
algorithm their influence on the efficiency of the search has to be investigated.

Table 2. Characterization of the computational effort necessary to determine the absolute minima of the various
test functions. All values are averages over 100 trial runs. For more information see text.

Number Number
Number of
Number of of first of second
Test function NDIM function
steps derivative derivative
evaluations
calculations calculations
2 41.35 42.07 6.88 2.76
10 191.64 143.95 29.98 61.09
Rastrigin (Rn) 20 273.38 220.40 33.30 67.72
30 347.59 293.36 34.23 69.15
50 530.28 470.43 38.17 76.06
2 155.81 138.4 63.54 24.71
10 527.29 499.05 185.12 38.58
Ackley (AKn) 20 839.96 793.17 214.9 62.14
30 1009.73 949.12 227.52 80.74
50 1522.36 1335.84 228.74 114.19
2 1446.92 269.32 132.49 1236.36
Griewangk (Gn) 10 1490.95 944.21 669.85 568.78
20 1978.73 1062.07 207.7 942.18
4 5917.58 7166.64 4586.28 868.9
Levy (Ln)
7 5740.58 5716.27 4931.76 711.74
Branin (BR) 2 34 26.83 11.67 9.4
Goldstein-Price (GP) 2 56.38 49.27 16.13 10.08
Hansen (H) 2 563.86 387.52 117.46 211.17

2.1.3.1 Parameters
The parameters are summarized in the Table 3 together with brief explanations.
Parameters as the dimension of the problem (NDIM), or the upper and low bounds of the
search space (Rup, Rlow) are determined by the concrete task. Taken the example given in the
introduction (conformational search for large molecules) NDIM would be the number of
freely rotatable bonds and Rup and Rlow would be the physically meaningful range in which
they can be rotated. The parameters rankmax and rankmin are ranking parameters from the

-60-
linear ranking procedure [Eqs. 9 and 10]. They were set to 0.1 and 1.0 respectively, as
discussed earlier. The parameters Itermain, Iterworst, Iterloc, and IterMAS are loop termination
numbers. They restrict the number of allowed iterations but do not influence the convergence
of the new algorithm. The number of elements in the TL (LSIZE) should not be set too small
to avoid that the search spins around but with LSIZE= 10*NDIM no problems were found.

Table 3. Description of the parameters of the approaches.

Recommended
Parameter Purpose
values
Parameters determined by the problem
NDIM dimension
Rup upper bound for each variable (see appendix A)
Rlow lower bound for each variable (see appendix A)
determined by the
Itermain
problem
Iterworst
loop termination numbers
Iterloc
IterMAS
Parameters with negligible influence
LSIZE number of elements in tabu list 10*NDIM
rankmax maximum recency ranked value 1.0
rankmin minimum recency ranked value 0.1
Ntrial number of trial points at the diversification search 3*NDIM
Parameters with strong influence a
calculates with
coeff coefficient used to compute the radius RTR (Eq. 7)
respect to RTR=0.1
alam first step size at the line search 1.0
Δxi step size at the mildest ascent strategy 0.1
α TD coefficient, used in formula (Eq. 8) 0.4

The parameters Ntrial and coeff are only connected with the diversification search. In
our test cases, Ntrial was always set to three times of the dimensionality of the problem
(3*NDIM), but not larger than 250. The parameter RTR, which is computed with the help of
coeff, Rup, and Rlow (Eq. 7), controls the size of the TR within diversification runs. If RTR is too
large, the TR become so large that minima lying close by already visited points are
overlooked since they lie in a TR. Naturally, this happens most frequently for problems with

a
investigated in Table 4 (for GTS) and Table 8, Figure 2-15, Figure 2-16 (for GOTS and TSPA)

-61-
many close lying minima. Since some information about the overall shape of the surface can
easily be seen in a first scan, RTR can be adjusted quite easily to the problem at hand. To
ensure that no important minima are overlooked, one could make a test run with a quit small
value of RTR starting the search from the lowest minima found in previous runs. For the
Rastrigin and the Ackley functions, the influence of RTR is depicted in Table 4. It is seen that
RTR influences the convergence to some extent. It is also seen that if RTR becomes too large,
the absolute minimum is missed. Therefore, a smaller value for RTR should be chosen. For the
present cases, RTR =0.1 was always chosen.

Table 4. Influence of the user-defined parameters for thirty dimensional functions observed during 10 trial runs.

Rastrigin Ackley

Parameter Value Number of Number of

function Success function Success
evaluations evaluations

1.0 412 100% 1479 50%

RTR 0.2 478 100% 1269 100%
0.1 293 100% 949 100%

2 277 100% 984 100%

alam 1 293 100% 949 100%
0.2 302 100% 1330 100%

0.5 - 0% 1570 80%

Δxi 0.2 281 100% 1384 100%
0.1 293 100% 949 100%

0.8 671 100% 1095 100%

α 0.4 293 100% 949 100%
0.1 218 100% 1011 100%

The influence of the fine tuning parameters alam, ∆xi, and α, which directly influences
the search part, is also depicted in Table 4 using the Rastrigin and the Ackley function (each
with a dimension of 30) as examples. The parameter alam represents the first step size in the
line search performed in the Steepest Descent and Quasi-Newton part. For small values, the
local minimum search becomes more accurate, but more steps are necessary. Additionally, for
a larger step size narrower and shallower minima are missed. This is no disadvantage, since
the deeper minima is of main interest. Decreasing alam from 2 to 0.2 the absolute minima are

-62-
always found but the necessary number of function evaluations increases (8% for Rastrigin,
35% for Ackley). As expected, the correlation between the number of steps and the value of
alam is very similar for all cases. Although not optimal for all test computations discussed
alam = 1.0 (full Newton step) was used for all of them. The parameter ∆xi gives the step size
during the mildest ascent search. Its influence is not surprising, since the search for a saddle
point is a quite delicate task. One would expect that the number of steps necessary to reach
the global minimum increase if a smaller step size, ∆xi, is chosen. But this is only found for
the Rastrigin function. Additionally, while for increasing alam the global minimum was
always found, this is not the case for ∆xi. Therefore, ∆xi should not be set too large, since the
next saddle point could be missed. For the following test cases ∆xi was always set to 0.1.
The parameter α (0.0 < α ≤ 1.0) determines the size of the cone of the TD (Eq. 8,
Figure 2-9). One would expect that the number of steps decreases if a wider cone is set tabu;
however, the opposite is found. This unexpected behaviour results since large cones restrict
the flexibility of the search, i.e. the search must take detours to the global minimum since
more direct ways are set tabu by accident. This is in line with the finding that for quite large
cones, the minimum is not found anymore. For all test cases discussed later, α was set equal
to 0.4. Again, this value was not optimized by preliminary test runs.
In summary, Table 4 indicates that the new approach is quite stable with respect to the
actual values of the fine tuning parameters. In a wide range, it finds the global minimum also
for nonoptimal fine tuning parameters; only (for some functions) the number of necessary
steps increases.

2.1.3.2 Numerical Results

Table 2 summarizes the computational efforts, which were necessary to determine the
global minima of the various test functions. In most real life applications of optimization
routines, it remains unclear if the global minimum is actually found. Therefore, most
optimization routines possess stop criterions, which control when a diversification is made
and when the program stops. In the present implementation, the program performs a
diversification run if after three new local minima no better solution was found. The program
finally stops after Itermain runs were performed. For the test cases for which the global minima
are well known, Table 2 lists the number of steps that were made until the global minimum
was reached for the first time. This number was determined by an analysis of the run. For the
computation double precision variables were used. The global minima were always reached
within numerical accuracy.

-63-
Table 5. Comparisons with the literature for the Rastrigin function. The column “Evaluations” gives the function
evaluations. For the present work it also comprises the sum of function, gradient and Hessian computations as
second number.

NDIM Algorithm Evaluations Source

Pure random search 5964 315
Fast annealing evolutionary algorithm 544 314
Random tunneling algorithm 383 314
TRUST 59 314

2 Tabu Search 540 291

Differential evolution 489 324
Random search 610 324
Simulated Annealing 260 324
Gradient Tabu Search 42/61 333
Tabu Search-MS 400 11624 316
Tabu Search-REM 1981 316
10
Genetic Algorithm 4488 319
Gradient Tabu Search 143/1053 333
Genetic Algorithm >100000 321
Cooperative Coevolutionary GA 20000 321
Tabu Search-MS 200 60446 316
Tabu Search-REM 7961 316
20
Parallel Genetic Algorithm 9900 322
EASY Genetic Algorithm 6098 323
Breeder Genetic Algorithm 3608 323
Gradient Tabu Search 220/2240 333

30 Gradient Tabu Search 293/3394 333

50 Gradient Tabu Search 470/6181 333

The number of variables of the various functions was varied to study the behaviour of
the new approach with respect to the dimensionality of the problem. In each case, the present
approach found the absolute minimum. Naturally, the number of steps needed to determine
the absolute minimum depends on the starting value. To obtain statistical results, for each
case (function and dimension of the problem), 100 independent trials were carried out. Each
trial corresponds to a new starting solution, chosen randomly within the search space.

-64-
 For the Rastrigin function with two variables, the present approach needs about 42
steps to find the global minimum. Within the search, about 42 function evaluations and seven
first derivative evaluations are performed. The Hessians are computed about 3 times. Please
note that in each first and second derivative evaluation, a vector of the dimension of the
problem has to be computed. As expected, the number of steps increases with the dimension
of the problem; however, the increase is considerably less than linear. If the dimension is
increased by a factor of 25, the number of steps increases only by a factor of about 13. If the
number of variables is below 30, the wall clock times needed to determine the global
minimum are well below 1 s. For 30 and 50 variables, 1–2 s were needed.

Table 6. Comparison with the literature for the Ackley function. The column “Evaluations” gives the function
evaluations. For the present work it also comprises the sum of function, gradient and Hessian computations as
second number.

Wall clock
NDIM Algorithm Evaluations Source
time [sec]
Random search 1257 324 3.3
2 Differential evolution 311 324 0.4
Gradient Tabu Search 138/315 333 <<0.1
Differential evolution 4214 324 14.8
10 Gradient Tabu Search 499/2736 333 <<0.1
Differential evolution 9957 324 80
20 Gradient Tabu Search 793/6334 333 <<0.1
Genetic Algorithm >100000 321
Cooperative Coevolutionary GA 50000 321
Breeder Genetic Algorithm 19420 323
Easy Genetic Algorithm 13997 323
30 Tabu Search-MS 100 22842 316
Tabu Search-REM 17941 316
Differential evolution 14957 324 222.7
Gradient Tabu Search 949/10197 333 ≈1
Differential evolution 22557 324 888.9
50 Gradient Tabu Search 1336/18482 333 ≈2

For the Ackley function, the new approach needs more steps to find the global
minimum but the increase of the number of steps with respect to the dimensions is again
considerably lower than linear. It is important to note that the number of first and second
derivate computations even increases less than the number of function evaluations. The Levy
function seems to be the most complicated case. To determine the global minimum for four
variables, the present approach needs already about 5917 steps. This is considerably more

-65-
than it needed for all other test cases. For the Branin, the Goldstein-Price, and the Hansen
functions, a similar behaviour as for the Ackley and Rastrigin function is obtained. For the
Griewangk function, the number of second derivative calculations decreases considerably if
one goes from 2 to 10 dimensions. Additionally, the number of steps increases quite slowly.
This happens since during the gradient optimization, the search jumps over various local
minima. It seems that these local minima are surrounded by quite narrow valleys.
Nevertheless, the global minimum was found in each run.
Table 5 - Table 7 summarize comparisons with other optimization routines. An
objective comparison of the computational efforts needed to find the global minima is only
given by wall clock times taken on the same machine since the computational overhead for
the steps varies from method to method. Since this is nearly impracticable, the number of
function evaluation is normally taken for comparison. These are summarized in the column
Evaluations of the tables. The data obtained within the present study were averaged over 100
trial runs, but the numbers are rounded to integer. To estimate the effort for the present
approach it has to be considered that a gradient or Hessian evaluation comprises the
computation of a vector with the size of the dimension of the problem. Therefore, for the
present approach, this column also gives the sum of all computed elements (=number of
function evaluations + gradient evaluations * dimension + Hessian evaluations * dimension)
to get a better insight into the computational effort. This is upper limit for the effort since
various terms already computed for the function evaluation can be reused for the analytical
gradients and Hessians. Additionally to the numbers given in Table 5 - Table 7 the
convergence of the present approach is also compared with convergence of the differential
evolution algorithm (Figure 2-10, Figure 2-11).
The Rastrigin function was often employed to test the efficiency of optimization
routines291,292,314,315,316,319,320,321,322,323 . On the basis of the number of function evaluations, the
present implementation represents the most efficient algorithm if a Rastrigin function with
two variables is used. Even if the number of evaluations used in the present approach includes
gradient and Hessian computations only the Terminal repeller unconstrained subenergy
transformation algorithm used by Jiang et al.314 needs fewer evaluations. Figure 2-10
compares the convergences of the GTS and of the Differential evolution algorithm for the two
dimensional Rastrigin function. The Differential evolution algorithm runs were performed
with the Mathematica 324 . It clearly demonstrates the excellent convergence of the GTS, which
shows hardly any oscillation. The objective function of Differential evolution algorithm
decreases much slower. Additional quite a few oscillations can be seen.

-66-
Table 7. Comparison with the literature for the Branin, the Goldstein-Price, and the Hansen functions. The
column “Evaluations” gives the function evaluations. For the present work it also comprises the sum of function,
gradient and Hessian computations as second number.

Function Algorithm Evaluations Source

Tabu Search-CJ 492 291
Pure random search 4850
Simulated Annealing 1 2700
Multi-start 1600 317
Simulated Annealing 2 505
Tabu Search 492
Paper work 144
BR TRUST 55 318
Tabu Search-REM 192
316
Tabu Search-MS20 333
Mod. Annealing 144 317
Simulated Annealing 324
Random search 15 324
Differential evolution 10
Gradient Tabu Search 27/69 333
Simulated Annealing 1 5439
317
Pure random search 5125
Multi-start 4400 318
Tabu Search-REM 894 316
Simulated Annealing 2 563
Tabu Search 486 318
Mod. Annealing 486
GP
Enhanced Continuous Tabu Search 231 290
Directed Tabu Search 230 286
Tabu Search-MS 500 195 316
Paper work 112 317
Tabu Search-CJ 112 291
TRUST 103 318
Gradient Tabu Search 49/102 333
Random search 420 324
Simulated Annealing 412 324
H
Differential evolution 291 324
Gradient Tabu Search 217/425 333

-67-
 The test runs made for a 10 and 20 dimensional Rastrigin function allow a comparison
with previous modifications of the TS291,316 and various genetic algorithms319,321,322,323. The
so-called Tabu Search-REM method316 is much better than the TS-MS 400 approach316;
however, the number of evaluations needed by the present approach is about a factor of 2–3
less. The present method also outperforms the genetic algorithms. For even higher
dimensional Rastrigin functions, no other results could be found.

- Differential
50
evolution
Objective function

40 - GTS

30

20

10

100 200 300 400 500

Number of evaluations

Figure 2-10. Comparison of the convergence of the algorithms for two dimensional Rastrigin function. The stars
indicate the values of the objective function for the local minima which were found during the GTS.

The thirty dimensional Ackley function was used as test case for some genetic
algorithms321,323 and improved TS approaches316. To get more information, test runs were also
performed with some algorithms implemented into the Mathematica 5.1 program package.
These tests were performed with the Random search7,325,326,327 and the Differential
evolution 328,329,330, 331,332
approaches. While the former were only used for the two
dimensional case, the latter were employed for all dimensionalities.
As already discussed the Ackley function seems to be a more difficult test case for the
present approach since it needs considerably more steps to determine the global minimum of
the Ackley function than for the Rastrigin function. Indeed, while for the Rastrigin function,
the present approach outperformed nearly all methods if the number evaluations were
compared, for the Ackley function it is still better than most other approaches but the
differences are smaller. For thirty dimensions, the Differential evolution approach
implemented in the Mathematica 5.1 needs 14957 function evaluations to find the global
minimum. For the easy genetic algorithm proposed by Voigt323 13997 function evaluations
were necessary. With 10197 evaluations (number of function evaluations + gradient

-68-
 evaluations * dimension + Hessian evaluations * dimension), the present approach needs only
a slightly smaller number. Please note that the number of function evaluation (949) and also
the approximated number of steps (1010) are much smaller. The graphical plot given in
Figure 2-11 compares the number of evaluation and function evaluation of the GTS with the
number of function evaluations used by the Differential evolution algorithm. Again, the
convergence of the GTS is better even if the number of all evaluation of the GTS is compared
with the number of the function evaluations of the Differential evolution algorithm.

- Differential
20
evolution

- GTS (func. eval.)

15
Objective function

- GTS (eval.)

10

2000 4000 6000 8000 10000 12000 14000

Number of evaluations

Figure 2-11. Comparison of the convergence of the algorithms for thirty dimensional Ackley function. The stars
indicate the values of the objective function for the local minima which were found during the GTS.

While the present approach and the differential evolution ansatz taken from the
Mathematica324 program package seem to be similarly efficient on the basis of the number of
evaluations, if real wall clock times are compared, our approach seems to be much faster (1s
in comparison to 222 s). This difference may result since for the present approach, the number
of additional operations (beside the evaluations themselves) is quite small. Furthermore, it
may result since many terms which appear in a function evaluation also occur in the gradients,
i.e. they have to be computed only once. However, possible overhead within the Mathematica
program, which represents a more general code for all kinds of computations, may also
contribute. The present approach also seems to be more stable with respect to the parameters
than the differential evolution ansatz. With the differential evolution ansatz sometimes the
global minimum could not be reached at all. For some problems, this already happened with
the standard parameters included in the Mathematica. In contrast, the present algorithm finds
the global minimum also for nonoptimal fine tuning parameters.
For the Hansen (H) and the Goldstein-Price (GP) functions, which are also often used
as test cases286,291,311,312,314,315,316,317,318, a comparison between the various methods leads to

-69-
similar results as found for Rastrigin and Ackley functions. For the Branin (BR) function, the
random search and the differential evolution show faster convergence than the GTS approach.

2.1.4 Conclusions
The GTS represents a metaheuristic, which tries to determine a global minimum of a
function by the steepest descent—mildest ascent strategy. It combines algorithms developed
to locate stationary points on quantum chemical hypersurfaces with the strategy of the TS. It
uses a combination of the Steepest Descent and Quasi-Newton approach for a fast
minimization to the next local minimum. To escape local minima via the mildest ascent, the
diagonal elements of the Hessian are employed. For the steepest descent and the mildest
ascent in each step, all variables of the function F(xi) are varied simultaneously. The actual
direction is determined by a weighting procedure, which exploits the gradient vector elements
and the diagonal elements of the Hessian. To ensure an efficient blocking of already visited
regions despite the lower number of steps the GTS introduces Tabu Regions and Tabu
Directions as elements for the Tabu List.
Test computations with the Ackley, the Branin, the Goldstein- Price, the Griewangk,
the Rastrigin, the Hansen, and the Levy functions in up to 50 dimensions show that the new
approach outperforms most previous approaches on the basis of the number of function
evaluations. Comparison with the differential evolution ansatz proves the efficiency of the
present approach on the basis of wall clock timings. For these tests, the differential evolution
implementation of the Mathematica324 program package was used.
The present approach uses analytical first and second derivations. Consequently, its
area of application is restricted to differentiable functions. Possible applications could lie for
example in minimization routines for force field parameters or conformational searches with
force fields.

2.2 TSPA and GOTS

In a preceding chapter the Gradient Tabu Search (GTS) 333 was introduced. It uses
analytical gradients for a fast minimization to the next local minimum and the diagonal
elements of the analytical Hessian to escape local minima. For the minimization a
combination of the Steepest Descent and the Quasi-Newton methods is used 334,335,336,337 . To
follow the modest ascent only the diagonal elements of the Hessian are employed but to
determine the direction they are weighted by a linear ranking procedure. Moreover, novel
concepts for the Tabu Directions (TD) were introduced to ensure an efficient blocking of
already visited regions. The GTS can be classified as a multi-start method. Test computations

-70-
for well known functions showed that this new approach converges much faster than other
methods taken from the literature. However, while the exploitation of gradients and Hessians
speed up the convergence of the TS considerably especially the Hessian and in some cases
also the gradients may be difficult to compute or are even not available. This can be tackled
by the two new approaches which are given in the present chapter. The GOTS (Gradient Only
Tabu Search) 338 still uses the gradients to achieve a fast convergence to the next local minima
but instead of the Hessian the approach uses a grid of function evaluations to follow the
modest ascent. In the TSPA (Tabu Search with Powell’s Algorithm)338 also the gradients are
replaced by a grid of function evaluation.
These two new nonlinear global optimization routines are also based on the Tabu
Search strategy which tries to determine the global minimum of a function by the steepest
descent - mildest ascent strategy. The new algorithms are explained and their efficiency is
again compared with other approaches by determining the global minima of various well-
known test functions with varying dimensionality. These tests show in most cases that the
GOTS possesses a much faster convergence than global optimizers taken from the literature.
The efficiency of the TSPA method is comparable to that of genetic algorithms.

2.2.1 Description of algorithms

The search is started from a current solution F(x1, x2,…, xNDIM), where NDIM
represents the dimensionality of the problem. A set of neighbours is generated by applying
“moves” according to a GTS scheme. Combinations of the Steepest Descent and the Quasi-
Newton methods335 which are known to be very efficient for locating the next local minimum
were already successfully applied in the GTS approach and were consequently also adopted in
the GOTS. As the GTS, the GOTS starts with the Steepest Descent approach and switches to
the Quasi-Newton method if the gradient value becomes small enough. Since gradient may
sometimes not be available an approach solely based on function evaluation was also
developed. In the present implementation the Direction Set Method with Powell’s
Algorithm337 was used, whose choice of successive directions does not involve explicit
computations of the function’s gradient. The main idea of all Direction Set methods is using
the some instructions for updating the set of directions as the method proceeds. This set either
includes some very good directions that leads far along narrow valleys, or includes some
number of “non-interfering” directions with the special feature (more conventionally called
conjugate directions) that minimization along one direction is not “damaged” by subsequent
minimization along another one. In the case of a set of linearly independent, mutually
conjugate directions a boundless cycling through the set of directions can be avoided. Then,

-71-
one pass of line minimizations will put it exactly at the minimum of a quadratic form
function. For functions that are not exactly of quadratic form, it does not exactly lead to the
minimum; but repeated cycles of line minimizations in due course converge quadratically to
the minimum.
To reduce the costs the previously developed GTS only exploits the diagonal
elements of the Hessian to estimate step size and direction according to a linear ranking
procedure. In the GOTS and TSPA only function evaluations are used to estimate step size
and directions to reduce expenses concerned with necessity of Hessian calculation. For the
first step of the modest ascent all functional values Fz i+ and Fz i− are computed to determine
the general direction of the following moves. They are obtained by varying each dimension
by a user defined step size Δxi0 :

(
Fz i+ = F x1 , x 2 , K , x j + Δx 0j , K , x NDIM )
(
Fz i− = F x1 , x 2 , K , x j − Δx 0j , K , x NDIM ) (11)

The most advantageous direction for each variable is then determined according to:

+ , if Fz i+ < Fz i−
Di = (12)
−, if Fz i+ > Fz i−

The functional values Fzi needed for the linear ranking procedure and for further correction of
the direction are then computed with

(
Fz i = F x1 , x 2 , K , x j + D j × Δx 0j , K , x NDIM ) (13)

Then step size and direction are estimated according to

⎛ Fz min − Fz i ⎞
rank i = rank min + (rank max − rank min )⎜⎜ ⎟⎟
⎝ Fz max − Fz min ⎠ (14)
xi +1 = xi + Δxi × rank i
0

Please note that the Di’s are only needed for the first step of the modest ascent part since the
general direction is clear afterwards. Figure 2-12 illustrates the choice process of the most
promising direction for each variable and moving with the ranked steps that gives the
possibility to discard the unpromising solutions and decreases the number of function
evaluations needed to escape the optimality valley.

-72-
 -350

-400

-450
-200

2
4
-500

-550

-300

3
-600

-500 -450 -400 -350 -300

b
-350
-400
3
-400

-450

2
-500

-500

1 -550

-600

-600 -500 -400 -300 -200

a -600

-500 -450 -400 -350 -300

c
Figure 2-12. Modest ascent strategy of GOTS and TSPA on contour plot of the Schwefel function. a) the search
path: start solution 1 → local minimum 2 →neighbourhood solution 3→next local minimum 4; b) computation
of the functional values around the local minimum at the beginning of the modest ascent strategy; c) moving
with the ranked steps within the promising area.

In Eq. 14 rankmax and rankmin can be user defined parameters but are normally set to
0.1 and 1. Fzmax and Fzmin are the maximal and minimal Fzi values, respectively. If the
differences between maximal and minimal function values become too small this ranking
scheme leads to too small values for the larger factors ranki. Therefore, in cases with
Fzmin-Fzmax < NDIM, rankmin is redefined as
⎛ Fz − Fz max ⎞
rank min = rank max − ⎜ min ⎟ (15)
⎝ NDIM ⎠
This strategy is followed until the new calculated function value is smaller than the
previous one which indicates that the barrier to the next valley is crossed. From this point the
next local minimum is located using the minimization strategy explained above.
To avoid reverse moves and cycles within the search the TL, the TR, and the TD as it
shows Figure 2-13 are used which were already successfully employed in the GTS method333
and which related to ideas introduced by F. Glover 339 . The TR and the TD are defined as in
the GTS method (Eqs.7, 8).

-73-
 ( Rup − Rlow )
1, 2, 3 and 4: RTR =
coeff
Tabu List solutions

4
4

3
2

0 2

-2
1

-4

-4 -2 0 2 4 6

(1 − α ) < amin bstart ⋅amin cnew ≤ 1

Figure 2-13. TL, TD, TR.

The last element of the GOTS and the TSPA represents the DS for the detection of
promising areas to select new starting points in the solution space. For GOTS and TSPA
diversification is performed in the same way as in GTS.
Initialization Locate+Store
//Start from randomly chosen point //Find next local minimum
Locate+Store Steepest Descent and/or Quasi-Newton
/or Powell’s methods
Main Loop Add minimum solution to the Tabu List
FOR ( iter = 0; iter < ITER _ MAX ; iter + + ) { Testing lower and upper bounds
//Escape local minimum IF (variables out of range) {
Mildest Ascent Search Diversification search
Locate+Store } GOTO Locate+Store
}ELSE{ Store minimum in the Result List}
IF (solution does not improve
after a number of iterations){
Mildest Ascent Search Diversification search
//Select the general move direction GOTO Locate+Store }
+, if Fzi+ < Fzi− IF (we obtained the same best results){
Di = Diversification search
−, if Fzi+ > Fzi− GOTO Locate+Store }
WHILE( F ( xi ) > F ( xi −1 ) ){
//Move in the direction of the lowest Di Next X
Next X (xi+1) FOR ( i = 1; i ≤ NDIM ; i + + ){
IF (xi+1 not TABU ){ Fzi = F ( x1 , x2 ,..., x j + D j × Δx 0j , x j +1 ,..., x NDIM )
GOTO Next X
}ELSE { next non-tabu Mildest Ascent} //Compute new variable vector
} ⎛ Fz max − Fzi ⎞
ranki = rank min + (rank max − rank min )⎜⎜ ⎟⎟
⎝ Fz max − Fz min ⎠
xi +1 = xi + Δxi0 × ranki
}

Diversification search Depend on problem

Figure 2-14. Pseudo-code of GOTS and TSPA.

The flowchart of both methods is given in Figure 2-1. Sometimes it is efficient to use
DS as an initialization unit for reasons given above. After the initialization one starts with the
search for improved minima. This part consists of local searches that are applied to obtain the
nearest local minima and the modest ascent searches used to escape to the next valley. For the
modest ascent search weighted function values can be used. In between the solution vector
and the Tabu List are updated. If the solution does not improve after a given number of
iterations or if all neighbourhood solutions are already set tabu the search for improved
minima is aborted and the diversification search is performed to obtain new starting points. A
pseudo-code of the GOTS and the TSPA algorithms is presented in Figure 2-14.

2.2.2 Investigation of the efficiency of the new approaches

Various well known multi-minima functions were
285,291,292,311,312,313,314,315,316,317,318,340
used to prove the efficiency of the GOTS and the TSPA
approaches. It also allows comparing their abilities with already known methods. In addition
to test functions which were used for the investigation of the GTS method the Rosenbrock

-75-
(Rbn), the DeJoung (DJ), the Zakharov (Zn), and the Trid (Trn) functions were tested.
Additionally, a functional form which was introduced to mimic conformational searches341
was used to test the efficiency of the GOTS method for such problems. The corresponding
expressions are given in Appendix A where the global minima of the functions are also given.
One measure for the efficiency of a global optimization routine is the number of steps
(operations) being necessary to reach the global minimum of the above mentioned functions.
For the present study such numbers were obtained by running the simulations over a specific
period and determining afterwards how many steps were made until the global minimum was
passed for the first time. For the computation double precision variables were used and the
global minima were always reached within numerical accuracy. The convergence of the
optimization routines which is another measure for the efficiency of the algorithm was also
tested.

2.2.2.1 Parameters
The proposed algorithms have several parameters and for a solid assessment of the
success of the new algorithms their influence on the convergence of the optimization routine
has to be tested. The parameters are reported in Table 3. The description of parameters and
their influence on the problem are already given in the preceding chapter on the GTS
approach. The test results for the GOTS and the TSPA are given below.
The parameters rankmax and rankmin are ranking parameters from linear ranking
procedure (Eqs. 14, 15). They were always set to 0.1 and 1.0. A noticeable influence on the
convergence was only found for parameters RTR, alam, Δxi, and α. Their effects are studied in
Table 8 and illustrated in Figure 2-15 - Figure 2-18 using the Rastrigin and the Ackley
function as examples. The values for the GTS which were taken from the previous chapter333
are given for an easier comparison. To estimate the influence on the convergence the tables
show how often the function, its gradient, and the diagonal elements of its Hessian had to be
computed until the global minimum was reached. The numbers are averaged over 10 test runs.
The way of the calculation or choice of the parameter RTR is identical to the manner in the
GTS approach. Due to these reasons a value of RTR = 0.1 seems again to be a good choice for
a large variety of different problems in the cases of the GOTS and TSPA methods.
A similar behaviour is found for the parameter α (0.0 < α ≤ 1.0) which determines the
size of the cone of the tabu directions333. Similar to large tabu regions too large cones restrict
the flexibility of the search. The search must take detours to the global minimum since more
direct ways are set tabu by accident. This is in line with the finding that for quite large cones
the minimum is not found anymore (not shown in Table 3). Within our tests α =0.4 was found

-76-
to be a good choice. This value was used for all test cases for GOTS and TSPA discussed
below.
The value of the first step size in the line search performed in the Steepest Descent and
the Quasi-Newton part is defined by the parameter alam. It therefore influences only the GTS
and the GOTS. For small values the local minimum search becomes more accurate, but more
steps are necessary. Increasing alam from 0.2 to 2.0 the necessary number of function and
gradient evaluations decreases considerably while the absolute minimum of both functions are
still always found. The speed up may result since for a larger step size narrower and shallower
minima are missed. Since deeper minima are of main interest this is no disadvantage. With a
value of alam = 1.0 (full Newton step) one is on the safe side for all problems tested so far
without a tremendous loss of efficiency.

-77-
 Rastrigin function
(thirty dimension)
Influence of Δxi parameter on Influence of α parameter on
total number of evaluations total number of evaluations

15000 15000

10000 0.1 10000 0.1

0.3 0.4
5000 5000
0.5
0.5 0.8
0.8
0 0.3 0.4
0.1 0
0.1
TS PA GO TS TS PA
GTS GO TS GTS

Figure 2-15. Influence of Δxi and α parameters on the sum of evaluations for the thirty dimensional Rastrigin
function.

Ackley function
(thirty dimension)
Influence of Δxi parameter on Influence of α parameter on
total number of evaluations total number of evaluations

Success 2%
15000 15000

10000 0.1 10000 0.1

0.3 0.4
5000 5000
0.5
0.5 0.8
0.8
0 0.3 0.4
0
0.1 0.1
TS PA GO TS TS PA GO TS
GTS GTS

Figure 2-16. Influence of Δxi and α parameters on the sum of evaluations for the thirty dimensional Ackley
function.

-78-
Table 8. Influence of user-defined parameters for the thirty dimensional functions.

TSPA GOTS GTS

Function Parameter Value
Function Function Gradients Function Gradients Hessian
1.0 20095 1398 39 863 73 49
RTR 0.2 15977 1386 40 366 34 42
0.1 16128 1295 40 293 34 69
Rastrigin
2 1423 48 388 38 62

alam 1 - 1295 40 293 34 69

0.2 1362 74 387 56 62

1.0 15974 1509 233 1596(60%) 277 60

RTR 0.2 13769 1211 188 1173 219 108
0.1 13033 3230 317 949 228 81
Ackley
2 1273 206 1095 271 58

alam 1 - 1189 186 949 228 81

0.2 1789 464 1177 406 114

-79-
 The dependency of the total number of calculations (number of function evaluations +
number of gradient evaluations) on Δxi is shown in Figure 2-15 and Figure 2-16. The
parameter Δxi gives the step size during the modest ascent search. With a larger step size
fewer steps are necessary to leave a local minimum; however, the number of steps needed to
reach the global minimum could rise anyway since the search for the lowest transition state to
the next local minimum becomes less accurate. This is seen if the GOTS is used to find the
global minimum of the thirty-dimensional Rastrigin function. For Δxi= 0.1 1295 function
evaluations and 40 gradient evaluations are necessary to reach the global minimum. For Δxi=
0.3 only 647 function evaluations and 33 gradient evaluations are necessary but for Δxi= 0.5
the effort again increases to 1396 and 79. The influence of Δxi on the speed and the cost of the
optimum reaching is large; the global minimum can even be missed as seen for the GTS
method. The value Δxi= 0.1 is recommended as a standard value.
For the Ackley function similar trends are found for all discussed parameters. In
summary, Table 8, Figure 2-15 and Figure 2-16 indicate that the previous GTS and the new
approaches GOTS and TSPA seem to be quite stable with respect to the actual values of the
fine tuning parameters. Only the number of necessary steps varies, but the global minimum is
always found.

2.2.2.2 Tests of the efficiency

Table 9 allows a comparison of the three different Tabu Search algorithms the TSPA,
the GOTS, and the GTS. The values for the GTS333 are given for an easier comparison. As
test cases different functions often applied for such purposes were used and the
dimensionality of the functions was varied. As Table 8, Table 9 gives how often the function,
its gradient, or the diagonal elements of its Hessian had to be computed until the global
minimum was reached. The computation of the gradient or the diagonal elements of the
Hessian represents a calculation of a vector of the dimension NDIM. To take this effort into
account Table 9 also contains the columns “Sum” which give the sum of all necessary
evaluations. For the GOTS “Sum” is given as the number of function evaluation plus the
number of gradient evaluation times NDIM. For the GTS “Sum” additionally includes the
number of Hessian evaluations times NDIM. To perform the tests on an equal footing the
same parameter set were used for all approaches (α=0.4, Δxi=0.1, RTR=0.1, alam = 1.0). For
each case (function and dimension of the problem) 100 independent trials were carried out
and averaged to obtain statistical results. Each trial corresponds to a new starting solution,
chosen randomly in the search space.

-80-
Table 9. Comparison of the three different Tabu-Search algorithms. All values are averaged over 100 trial runs. Percentages in brackets indicate the probability to find the
global minimum.

TSPA GOTS GTS

Test function NDIM
Function Function Gradient Sum Function Gradient Hessian Sum
2 230.4 46.18 10.95 68 42.07 6.88 2.76 61
10 3765.24 304.07 28.21 586 143.95 29.98 61.09 245

Rastrigin (Rn) 20 9084.35 820.17 33.92 1498 220.40 33.30 67.72 2241
30 16128.21 1295.08 40.36 2506 293.36 34.23 69.15 3394
50 29604.20 2412.14 50.08 4916 470.43 38.17 76.06 6182
2 265.33 214.45 78.50 372 138.40 63.54 24.71 315
10 3366.52 826.98 196.91 2796 499.05 185.12 38.58 2736

Ackley (AKn) 20 9100.18 1662.23 253.16 6725 793.17 214.90 62.14 6333
30 13033.77 3230.11 317.05 12742 949.12 227.52 80.74 10197
50 20764.25 3808.75 225.42 15080 1335.84 228.74 114.19 18482
2 2364.70 1310.44 165.67 1641 269.32 132.49 1236.36 3007

Griewangk (Gn) 10 48960.60 1831.97 200.35 3835 944.21 669.85 568.78 13330
20 61902.11 540.60 13.60 812.6 1062.07 207.7 942.18 24060
4 21994.00 (5%) 11611.08 (75%) 5817.01 34879.1 7166.64 4586.28 868.90 28987
Levy (Ln)
7 37778.22 (9%) 15000.81 (54%) 10064.83 85454.6 5716.27 4931.76 711.74 45221

Branin (BR) 2 294.45 284.03 23.65 331 26.83 11.67 9.40 69

Goldstein-Price (GP) 2 3833.19 79.87 20.15 120 49.27 16.13 10.08 102

Hansen (H) 2 2872.55 614.14 154.31 923 387.52 117.46 211.17 1046

-81-
 Overall Table 9 shows that the GTS and the GOTS are considerably more efficient
than the TSPA. This could be expected since gradient-based approaches are known to be
much more efficient for minimization problems. The difference between both the GOTS and
the GTS depends on the function and its dimensionality (NDIM). Generally, the GTS seems to
be more efficient for small NDIM but for larger NDIM the GOTS becomes the most efficient
approach. The differences between the numbers of gradient evaluations indicate that both
approaches take different ways to the global minimum (see below).

35000

30000

25000
Sum of Evaluations

20000

15000

10000

5000

0
0 10 20 30 40 50
Dimension

TSPA (∆x=0.1) GOTS (∆x=0.1) GTS(∆x=0.1)

Figure 2-17. Dependency of sum of evaluations on the dimension of the Rastrigin function.

For the Rastrigin function the GTS requires a lower number of evaluations than the
GOTS for NDIM=2 and 10 but for the higher dimensionalities the GOTS becomes more
effective. In all cases the dependency of the sum over all evaluations (column “Sum”) on
NDIM is about linear but with varying prefactors. For the Rastrigin function (Figure 2-17) the
sum of all evaluations needed by the GOTS increases by a factor of about 8 if the dimension
increases from 2 to 10. Going from NDIM=10 to 50 the sum of all evaluations increases from
586 to 4916, i.e. again by a factor of about 8. The GTS approach behaves very similar. For the
TSPA algorithm the dependency is much stronger. Going from 2 to 10 dimensions the sum of
evaluations increases by a factor of 16, i.e. the increase is nearly quadratic. However, going
from 10 to 50 dimensions the sum of evaluations increases only by a factor of 8.
-82-
 25000

20000
Sum of Evaluations

15000

10000

5000

0
0 10 20 30 40 50
Dimension

TSPA (∆x=0.1) GOTS (∆x=0.1) GTS (∆x=0.1)

Figure 2-18. Dependency of sum of evaluations on the dimension of the Ackley function.

Differential Evolution GOTS

GTS TSPA

Figure 2-19. Comparison of the convergence of the algorithms for thirty dimensional Ackley function. The
symbols indicate the values of the objective function at the local minima which were found during the GOTS.

For the Ackley function similar trends are found but the slope of the dependency on
the dimensionality (Figure 2-18) is considerably steeper indicating that this function

-83-
represents a more complicated problem. This is also seen from the absolute numbers given in
Table 9. An insight into the convergence of the various approaches is also given in Figure
2-19. It sketches the computed function values of the Ackley function along a given search.
All approaches start from the same point. Of course the taken routes to the global minimum
depend on the starting points but Figure 2-19 gives a representative example. It is directly
seen that although starting from the same starting point GTS, GOTS, and TSPA take quite
different routes. The route taken from the GOTS converges most rapidly and already after
about 3000 evaluations the functional values decreases to about 1. For the GTS about 3500
evaluations are needed and for the TSPA more than 6000 evaluations are necessary. Figure
2-19 also proves that all approaches are capable to escape low lying regions quite efficiently
and shows that they converge more rapidly than the differential evolution
algorithm328,329,330,331,332 implemented in the Mathematica324 (see below).
Table 9 shows that the Griewangk function seems to be ideally suited for the GOTS. It
finds the global minimum very efficiently. Going from 10 to 20 dimensions the number of
necessary evaluations even decreases. This happens since the search jumps over various local
minima which are very narrow. This does not represent a disadvantage since only the global
minimum, which was found in each run, is be of interest. For the Griewangk function the
GTS and the TSPA need considerably more evaluations than for the Rastrigin or the Ackley
functions. This again shows that the three approaches take completely different routes
although similar strategies are used.
The Levy function seems to be a difficult case for all approaches. However, while the
GTS still finds the global minimum in all cases for NDIM=4 the TSPA finds the minimum
only in 5 % of all test runs. For NDIM=7 it finds the global minimum only in 9 % of the
cases. The GOTS finds the minimum in only 75 % (NDIM=4) and 54 % (NDIM=7). If
Δxi=1.0 instead of 0.1 is used the GOTS and the TSPA also always find the global minimum.
In this case for NDIM=4 the GOTS needs 3403 evaluations to reach the global minimum
while 10584 are needed for NDIM=7. The TSPA needs 18507 and 38930 evaluations,
respectively.
For the Branin function the outcomes are very similar to those for the Rastrigin
function. The Goldstein-Price function seems again to be more difficult for the TSPA. The
numbers indicate that the Hansen function is a difficult test function for all approaches. For
the two latter functions the influence of the step size is quite strong for the TSPA approach.
Increasing Δxi to about 1.0 the number of function evaluations drops to about 335 and 533,
respectively. GOTS tests for the DeJoung, the Rosenbrock, the Zakharov, and the Trid
functions also underline the efficiency of this approach (Table 13).

-84-
 Table 10 - Table 13 summarize comparisons with optimization routines taken from the
literature. An objective comparison of the computational effort needed to find the global
minima is only given by wall clock times taken on the same machine since the computational
overhead for the steps varies from method to method. Since this is nearly impracticable the
number of function evaluations is normally taken for comparison. They are summarized in the
column Evaluations of the tables. For GOTS and GTS the tables also give the sum over all
necessary evaluations (see Table 9). The data obtained within the present study were averaged
over 100 trial runs and the numbers are rounded to integer. A comparison of the convergences
of the GTS, the GOTS, the TSPA, and the differential evolution algorithm as implemented in
the Mathematica program324 were already given in Figure 2-19.

-85-
Table 10. Comparisons with the literature for the Rastrigin function. All values are averaged over 100 trial runs.
The column “Evaluations” gives the function evaluations. For the GTS it comprises the sum of function, gradient
and Hessian computations as second number. For the GOTS it also comprises the sum of function and gradient
computations as second number.

NDIM Algorithm Evaluations Source

Pure random search 5964 315
Fast annealing evolutionary algorithm 544 314
Random tunneling algorithm 383
318
TRUST 59
Tabu Search 540 291
2 Differential evolution 489
Random search 610 324
Simulated Annealing 260
GTS 42/61 333
GOTS 45/66
338
TSPA 230
Tabu Search-MS 400 11624
Tabu Search-REM 1981 316
Genetic Algorithm 4488
Differential Evolution 14017 331
10
Differential Evolution – Quasi Newton 13955 344
GTS 143/1053 333
GOTS 276/656
338
TSPA 3765
Genetic Algorithm >100000 321
Cooperative Coevolutionary Genetic Algorithm 20000 322
Tabu Search-MS 200 60446
316
Tabu Search-REM 7961
Parallel Genetic Algorithm 9900 322
20 EASY Genetic Algorithm 6098
323
Breeder Genetic Algorithm 3608
Differential Evolution 23870 331
Differential Evolution – Quasi Newton 15200 344
GTS 220/2240 333
GOTS 553/1387
338
TSPA 9084
GTS 293/3394 333
30 GOTS 1397/3776
338
TSPA 16128
GTS 470/6181 333
50 GOTS 1876/5079
338
TSPA 29604

-86-
Table 11. Comparisons with the literature for the Ackley function. All values are averaged over 100 trial runs.
All tests were carried out on the same computer. The column “Evaluations” gives the function evaluations. For
the GTS it comprises the sum of function, gradient and Hessian computations as second number. For the GOTS
it also comprises the sum of function and gradient computations as second number.

Wall clock
NDIM Algorithm Evaluations Source
time [sec]
Random search 1257 3243.3
324
Differential evolution 311 0.4
2 GTS 138/315 333 <<0.1
GOTS 164/300 <<0.1
338
TSPA 265 <<0.1
Differential evolution 4214 324 14.8
GTS 499/2736 333 <<0.1
10
GOTS 649/2536 <<0.1
338
TSPA 3367 <<0.1
Differential evolution 9957 324 80
GTS 793/6333 333 <<0.1
20
GOTS 968/5187 <1
338
TSPA 9100 <1
Genetic Algorithm >100000 321
Cooperative Coevolutionary GA 50000 322
Breeder Genetic Algorithm 19420
323
Easy Genetic Algorithm 13997
Tabu Search-MS 100 22842
30 316
Tabu Search-REM 17941
Differential evolution 14957 324 222.7
GTS 949/10197 333 ≈1
GOTS 1189/6768 ≈1
338
TSPA 13034 ≈1
Differential evolution 22557 324 888.9
GTS 1336/18482 333 ≈2
50
GOTS 1792/11785 ≈1
338
TSPA 20764 ≈3

-87-
Table 12. Comparisons with the literature for the Branin and the Goldstein-Price functions. All values are
averaged over 100 trial runs. The column “Evaluations” gives the function evaluations. For the GTS it comprises
the sum of function, gradient and Hessian computations as second number. For the GOTS it also comprises the
sum of function and gradient computations as second number.

Function Algorithm Evaluations Source

Tabu Search-CJ 492 291
Pure random search 4850
Simulated Annealing 1 2700
Multi-start 1600
317
Simulated Annealing 2 505
Tabu Search 492
Paper work 144
TRUST 55 318
BR Tabu Search-REM 192
316
Tabu Search-MS20 333
Mod. Annealing 144 317
Simulated Annealing 324
Random search 15 324
Differential evolution 10
GTS 32/69 333
GOTS 44/65
338
TSPA 203
Simulated Annealing 1 5439
Pure random search 5125 317
Multi-start 4400
Tabu Search-REM 894 316
Simulated Annealing 2 563
Tabu Search 486 317
Mod. Annealing 486
Enhanced Continuous Tabu Search 231 290
Differential Evolution 18123 331
Differential Evolution-Quasi Newton 14781 344
Simulated Annealing 563 320
Directed Tabu Search 230 286
GP
Tabu Search-MS 500 195 316
Paper work 112 317
TRUST 103 318
GTS 49/102 333
GOTS 52/73
338
TSPA 336
C-RTS-min 171 283
ECTS 231
DOPE 248 290
ESA 783
CTSSmultiple 119 312
CHA 259 342

-88-
 The Rastrigin function was often employed to test the efficiency of optimization
routines291,292,314,315,316,319,321,322,323. According to the number of function evaluations or also to
the number of total evaluations GTS and GOTS performs much better than the approaches
taken from the literature (Table 10). The Terminal Repeller Unconstrained Subenergy
Transformation algorithm318 includes gradient computations as GTS or GOTS but still needs
20 % more function evaluations. The number of total evaluations was not given in the
literature. The TSPA needs more evaluations but is still considerably faster than e.g. the Fast
annealing evolutionary algorithm or the original Tabu Search implementation291. The test runs
made for the 10- and 20-dimensional Rastrigin function allow a comparison with several
previous modifications of the Tabu Search291,316 and various genetic algorithms319,321,322,323.
The so called Tabu Search-REM method316 is much better than all previous Tabu Search
approaches316. Nevertheless, the GTS and the GOTS are still considerably faster. The TSPA
approach needs more function evaluations than the Tabu Search-REM. GOTS and GTS also
outperform the genetic algorithms. For even higher dimensional Rastrigin functions no other
tests could be found.
The thirty dimensional Ackley function (Table 11) was also used as test case for
genetic algorithms321,323 and improved Tabu search approaches316. Please note, that the GOTS
results were obtained with xi=0.5 while xi=0.1 was used for GTS and TSPA. All other settings
are equal to those employed in Table 9. To get more information the Random search7,325,326,327
and the Differential evolution328,329,330,331,332 approach implemented in the Mathematica 5.1
program package were also used. While the former were only used for the two dimensional
case the latter were employed for all dimensionalities. The differences between the various
approaches are similar to those found for the 20-dimensional Rastrigin function. For all
approaches the Ackley function seems to be a more difficult problem than the Rastrigin
function.

-89-
Table 13. Comparisons with the literature for the Hansen, the DeJoung, the Rosenbrock and the Zakharov
functions. All values are averaged over 100 trial runs. The column “Evaluations” gives the function evaluations.
For the GTS it comprises the sum of function, gradient and Hessian computations as second number. For the
GOTS it also comprises the sum of function and gradient computations as second number.

Function NDIM Algorithm Evaluations Source

Random search 420
Simulated Annealing 412 324
Differential evolution 291
H 2
GTS 217/425 333
GOTS 321/496
338
TSPA 534
ECTS 338 311
DOPE 131 294
DJ 3 CTSSmultiple 155 312
CHA 371 342
GOTS 23/62 338
ECTS 15720 (85%) 311
DOPE 8695 294
ESA 12403 290
DTSAPS 9037 (85%) 286
Rbn 10 Differential Evolution 67310 331
Differential Evolution-QN 16349 344
Simulated Annealing 12403 320
CGA 21563 (80%) 343
GOTS 379/1629 338
ECTS 4630 311
DOPE 5133 294
Differential Evolution 10088 331
Zn 10 Differential Evolution-QN 10013 344
Simulated Annealing 15820 320
CGA 6991 343
GOTS 61/161 338
GOTS 38/93 338
5
EM without LOCAL procedure 968 345
Trn
GOTS 125/745 338
20
EM without LOCAL procedure 43354 345
Deterministic algorithm 168 341
7
GOTS 131/271 338
Test function for Deterministic algorithm 497 341
conformational 10
search GOTS 192/412 338
Deterministic algorithm 174857 341
20
GOTS 2038/5378 338

-90-
On the basis of the number of evaluations GTS and TSPA seem to be similar in efficiency to
the differential evolution ansatz taken from the Mathematica 5.1 program package. However,
if real wall clock times are compared our approaches needs considerably less time (1 second
in comparison to 222 seconds). This difference may arise since for the present approaches the
number of additional operations (beside the evaluations themselves) is quit small, but it may
also result from the implementation and the employed computer. Another advantage of our
approaches is their stability with respect to the fine tuning parameters. They find the global
minimum also for nonoptimal fine tuning parameters. In contrast the differential evolution
ansatz sometimes misses the global minimum if the default parameters implemented in the
Mathematica 5.1 are used. The differences in the convergence were already discussed in the
context of Figure 2-19.
The investigations performed with the Branin, the Hansen, the Goldstein-Price, the
Trid, the Rosenbrock, the DeJoung, and the Zakharov functions lead to results similar to those
found for Rastrigin and Ackley function (Table 12-Table 13). For the Branin, the Hansen and
the Goldstein-Price functions was employed Δxi= 1.0 as more efficient for GOTS and TSPA.
The GOTS method was additionally applied to a functional form which was
introduced to mimic conformational searches341. The results are convincing (Table 13). The
global minimum is always found quite effectively. The larger numbers for NDIM=20 result
since the number of computations vary considerably. They range from about 900 to about
8000 within the 100 trial runs.

2.2.3 Conclusions
Two new global optimization routines; the Gradient only Tabu Search (GOTS) and the
Tabu Search with Powell’s Algorithm (TSPA) are described in this chapter. They are based
on the Tabu Search strategy which tries to determine a global minimum of a function by the
steepest descent - mildest ascent strategy. In this strategy the steepest descent is followed to
find the next local minimum. To escape the local minimum the mildest ascent is used.
Already visited regions are set tabu using a Tabu List which memorizes the already visited
regions. The GOTS uses a combination of the Steepest Descent and the Quasi-Newton
method to find the next local minimum. A step-wise approach is employed to realize the
modest ascent strategy. The TSPA uses the same implementation for the modest ascent
strategy but Powell’s algorithm for the steepest descent part, i.e. the TSPA only needs
function evaluations but no gradients. The algorithms are explained and their efficiency is
compared to other approaches by determining the global minima of various well-known test
functions. These tests show that the GOTS is more efficient than previous global optimizers.

-91-
For functions with higher dimensionality it is even more efficient than the previously
developed GTS which uses the diagonal elements of the Hessian to escape from a local
minimum via the modest ascent. The efficiency of the TSPA is comparable to Genetic
Algorithms. Test computations also indicate that the GOTS is an efficient approach for
conformational searches.

2.3 Conclusions
This chapter shows that the new approaches: GTS, GOTS, and TSPA can be
efficiently applied to the optimization of continuous multiminima functions. These techniques
can be useful for engineers or scientists, who are interested in the global solution but also
nearby local optima. Efficient results were obtained within a reasonable amount of
computational time. The proposed approaches are flexible and easy to use. Moreover, the
proposed approaches could be enhanced even further if more sophisticated descent method
and diversification strategies are applied instead.
The Figure 2-20 gives the graphical illustration of convergence of the all new
approaches described above in comparison with Differential Evolution, Linear and Quadratic
convergences. The dependency of sum of evaluations on the dimension of the Ackley
function for GOTS it is close to linear. It is obvious that the GOTS method achieves better
solutions than already known pure and hybrid global methods, as well as GTS and TSPA
methods.

GOTS
GTS
20260 TSPA
Diff. Evol.
15260 Linear
Evaluations

Quadratic

10260

5260

260
2 10 20 30 40 50
Dimension

Figure 2-20. Graphical illustration of convergence of the algorithms.

-92-
 In order to test the appropriateness of the GOTS method for molecular geometry
problems, the optimization of a function with a functional form similar to general potential
energy functions and whose global minimum is known was carried out. The number of its
local minima increases exponentially with the size of the problem, which characterizes the
principal complexity of the minimization of the molecular potential energy functions. The
utilization of a GOTS guarantees the global optimality and shows the difficulty in obtaining
the global minimum of such function. It was successfully applied to the n-chain problem (n =
7, … , 20) with decided superiority over Deterministic algorithm that uses the Branch and
Bound scheme with techniques of interval analysis to provide the lower bounds. In all the
cases the global minimum was obtained and the amount of function evaluations required
determining the global minima increases by a factor of about 11 if the dimension increases
from 10 to 20 whereas this factor for Deterministic algorithm is 352.
The examples show that the GOTS is able to provide solutions to difficult global
optimization problems with minimal computational effort. The algorithm can be suitably
modified to incorporate more rigorous mechanisms to handle very large-scale problems.
Because of problem-specific knowledge GOTS deals with variables and constraints
effectively. However, the trade-off between effectiveness and generality requires careful
consideration, and the junction with other techniques is necessary in order to built general
methodology. The recommended default parameters were given as appropriate for several test
cases. Large-scale test cases show that the initial parameter settings for GOTS may be
overcautious; however, in any case, they locate the global optima quite rapidly and with high
probability. The variable choice and parameter settings of GOTS is always problem
dependent. The optimal parameter values may require some adjustment to starting values
indicated in this chapter. As is obvious from the preceding discussion, GOTS determines the
best energy minimum of the molecular geometry problems efficiently. The next chapter
presents details for applying GOTS to chemical optimization problems on potential
hypersurfaces and provides a systematic approach to variables choice and setting the
adjustable parameters.

-93-
Chapter 3 Application and Discussion

3.1 Introduction
An efficient search for the global minimum of a highly dimensional function with
many local minima is central for the solution of many problems in computational chemistry.
Well known examples for such global optimization problems10 are the conformational search
for molecules with a high number of freely rotatable bonds 346,347 , the optimization of the
parameters of a force field or the determination of all possible reaction paths between
reactants and products 348,349,350 . Conformational analysis is the study of molecule
conformations and their influence on molecular properties. Fundamental for all
conformational analysis is the search for the most important conformers of the molecule. This
requires the location of the corresponding energy minima309,351 on the potential energy
surface (PES) 352 . In most cases the energetically lowest lying conformers are the most
important ones, since they are populated at e.g. room temperature and, hence determine most
of the properties of the molecule 353,354 .
Search for the energetically lowest lying conformer means to find the global energy
minimum of the energy hypersurface of the molecule. Hence it is a global optimization
problem in which the potential energy function is the objective function and the coordinates,
that are used to represent the conformation of the molecule, are the variables. For larger
molecules with many freely rotatable single bonds the number of minima may be so large that
it is nearly impossible to find all of them. In order to locate the global energy minimum or to
locate the lowest minima, a large number of starting conformations that are equally
distributed on the energy surface is generated. Each of them is minimized using local
optimization techniques 355 to the nearest local minima, and then all duplicated structures are
rejected. Even for small molecules the disclosure of the global energy minimum or all low
energy conformations in a multi-dimensional PES is a problem that requires considerable
computational effort. A complete search is unrealizable since the search space increases
exponentially with the number of degrees of freedom (e.g. torsion angle), which is typically
proportional to the size of the molecule. It can be referred to nondeterministic polynomial-
time hard problems 356 , where the total number of possible conformations increases
exponential with the total number degrees of freedom. The time required to solve it increases

-94-
exponentially with the size of the molecule. This is known as combinatorial explosion 357 .
Therefore, to explore the complete conformational space and to obtain all the low energy
conformations at tractable computational cost some specialized conformational search
algorithms are needed. Over the past several years, a multitude of conformational search
techniques have been developed for this purpose 358,359,360,361,362,363,364,365 . A brief review of
some of the methods is given below.

3.1.1 Conformational search techniques

Conformational search methods can be roughly divided into stochastic and
deterministic techniques. Stochastic methods are optimization algorithms which incorporate
probabilistic (random) elements to locate the global minimum either in the problem data or in
the algorithm itself, or in both. Hence there is no natural endpoint of the procedure. In the
deterministic approach the values of the objective function are assumed to be exact, and the
computation is completely determined by the values sampled so far.
• Systematic (or grid or exhaustive) search methods 366,367 posses as first step the generation
of a large number of starting geometries using combinations of selected values of variable
internal coordinates. In most cases only the torsion angles, which describe the rotation
around the single bonds are used. Each torsion angle is systematically incremented by a
fixed amount until all possible combinations of torsion angles for the chosen increment
have been generated. Then each structure is subjected to an energy minimization to yield
the corresponding low energy final conformations. The total number of search steps of
such searches must be predefined at the beginning of the search when the starting
geometries are generated. Except for simple molecules, the number of steps that are
necessary to find all conformations is unknown at the outset and may be so large that the
energy minimization stage of the search becomes prohibitive. The necessity of using
many values for each torsion angle reinforces the problem since low energy
conformations of flexible molecules commonly incorporate relatively high-energy
torsional arrays. Thus, the method is limited to very small molecules or molecule
segments.
• Build-up procedures 368 are based on the assumption that the conformations of different
fragments of the molecule are independent of the other ones. Comprehensive energy
minimizations are carried out for all fragments, and only those local minima whose
energies lie within an appropriately chosen upper bound of the lowest minimum are
retained. To obtain starting points for the energy minimization of larger fragments the
restricted set of local minima of one fragment is combined with the set from another

-95-
 fragment. If the assumption holds the global minimum of the whole chain can be built up
from its components by subsequent processes without the need to explore more than a
small part of the conformational space of the whole molecule. However, in the realization,
the build-up procedure has some difficulties. The main one is the exponential increase of
the number of minima which have to be retained at each step to obtain an acceptable value
for the cut-off energy. Additionally, two atoms can come too close if the fragments are
combined in an arbitrary manner. A final difficulty concerns the correct value for the cut-
off energy.
• Monte Carlo methods (the Metropolis algorithm) 369,370,371,372,373 are stochastic techniques
that are able to generate a random sample of low-energy conformations for molecules that
are too large and flexible to be explored systematically. In these methods the energy E0 is
calculated first for an arbitrary conformation. In each step of the algorithm, several torsion
angles (coordinates in Cartesian space) are randomly varied leading to a new
conformation with energy E. This new conformation is accepted as new starting point or
rejected depending on the rules described in the first chapter (Chapter 1 - Eq. 6). For
macrocyclic molecules, random variations of torsional angles often result in
conformations for which the ring-closure constraints are violated. Such structures cannot
be used for minimization and are rejected. The Monte Carlo step has to be repeated until a
set of low energy conformers has been generated. Several modified Monte Carlo methods
have been applied for peptides and proteins 374 .
• The Simulated annealing method219,375,376 is based on a connection between statistical
mechanics and the crystallization process. The main principles of this technique were
described in the first chapter. By analogy with this physical process, in a combinatorial
optimization context, a solution corresponds to a state of the physical system and the
solution cost to the energy of the system. At each iteration the current conformation is
replaced by a randomly selected trial conformation accepted according to the so-called
Metropolis criterion. The series of accepted steps is an “energy directed random walk”
which explores the conformational space. The simulated annealing technique exploits
both the energy and temperature dependency of the Boltzmann distribution. At a given
temperature, the Metropolis algorithm is used to imitate the conformational equilibrium
using the energy dependency of the Boltzmann distribution. Then at progressively
lowered temperature during the process according to some set cooling schedule is used to
profit by the temperature dependency of the Boltzmann distribution. The method is widely
used for the conformational search problems of molecules219,376.

-96-
• Internal coordinate tree search 377 is the more systematic version of an internal coordinate
Monte Carlo search and a rapid method for generating molecular geometries which
approximate the low-lying conformers of small to medium-sized organic molecules. It
generates starting conformations for subsequent minimizations by varying the torsion
angles of the rotatable bonds in series of possible values (e.g., 0°, 60°, 120°, 180°, . . .).
Bond lengths and angles are kept fixed. In the case of N rotatable bonds and d possible
values of each torsion angle up to dN conformations may be generated and minimized. In
worst case, every combination of all torsion angle values must be investigated. A torsion
angle ordering allows the complete conformer generation process to be graphed as a tree
structure. At the lowest level of the tree, the leaves stand for the resulting conformations
of the generation process. Each edge stands for the rotation of one torsional angle.
Intermediate nodes in the tree represent partially completed structures. For cyclic
molecules the conformations which violate the ring-closure constraint have to be
excluded. This method can be quite expensive for macrocyclic molecules but can lead to
the disclosure of previously non-described low-energy conformers. The internal
coordinate conformer generator has a number of advantages:
1. It is quite fast because it demands no energy gradient evaluations and no solution of
matrix of distance equations that are required e.g. for the distance geometry methods
(see below). That increases the speed of conformer generation which ranges from 10-
to 100-fold over molecular dynamics or distance geometry methods
2. Internal coordinates handling assures the sampling over all accessible regions of
conformational space
3. It is easily applied as a tree search, which enable analysis of acyclic and cyclic
molecules.
However, since this method is a probabilistic one it gives no guarantee that all low-energy
minima have been found.
• Genetic algorithms 378,379,380 follow standard procedures that are derived from the
principles of natural evolution. The main principles of the algorithms were already
described in the first part of the work. Here the explanation of the implementation
specificity of such methods for the conformation search is given. An initial (parent)
population is built by assigning each main-chain torsion angle to a value randomly chosen
from allowed ones. The fitness of each conformation is then calculated using the energy
of the minimized structure, i.e. each of the initial conformations is minimized using the
local minimization procedure. The optimized conformations form the first parent
generation as well as a pool of random conformations. New child conformations are

-97-
 produced from the parent population using standard genetic algorithm operations. The
probability that a parent conformation is selected for the next generation is proportional to
the Boltzmann factor of its energy so that low-energy conformations are more likely to be
selected than high-energy conformations. Each parent conformation is selected only once.
To keep up the diversity, before each crossover operation, one of the two selected parent
conformations can be replaced by a conformation, which has a small probability,
randomly chosen from the pool of random conformations. The parents selection and
crossover process is repeated until all conformations of the parent generation are selected.
The reason combining the parents and newly generated child conformations, which were
already sampled in the previous populations, with each other is the attempt to obtain
promising structures. This shall assure that the fitness of the new generation is not worse
than that of the parents. Then a “cleaning” process is preformed that consist in rejecting
conformations with a higher energy from any pair of structures which is similar to each
other. The next process iteration uses the new generation as parents. The cycle is repeated
until convergence. The best fit population is finally obtained or a maximum number of
iterations has been reached. The application of the genetic algorithm in peptides and
protein structure prediction was reviewed by Le Grand and Merz379, and Schulze-
Kremer380.
• Distance geometry methods 381,382,383,384,385 represent molecules as distance matrixes
containing lower and upper bounds of the distances between every pair of atoms. These
distance bounds consist of bonded constraints, such as bond lengths and angles, as well as
nonbonded constraints, such as van der Waals radii. It is significant that the distance
information alone is unable to differentiate between a structure and its mirror image.
Hence, chirality constraints are included to eliminate this disadvantage. The constraints
are given beforehand on the basis of the known stereochemistry of the chiral centres. The
distance geometry representation can be applied to cyclic molecules as well as to acyclic
molecules. At the first stage a distance matrix, which describes the structure of the
molecule, is created whose elements dij are the distances between atoms i and j. Then a
matrix with upper and lower bounds for each interatomic distance is determined. New
conformations are typically generated by assigning new interatomic distances which lie
between the upper and lower bounds. Then those conformations are generated that best
approximate these distances. The resulting conformations are finally energy-minimized.
Taylor and Aszodi381 applied this method to polypeptides.
• Smoothing/deformation methods rely on the assumption that the position of the global
minimum of a deformed energy hyper surface can be related to the position of the global

-98-
 minimum of the original function. The algorithm is based on the deformation of the
original potential energy hypersurface in such a way that the number of the local minima
is reduced. Only the deepest one is retained, which, in most cases, is related to the global
one. A local minimization procedure applied to the modified energy function leads from
any starting point to this single minimum instead of ending up in some higher-energy
local minimum. The hypersurface of a molecule is deformed (or smoothed) by applying
certain mathematical operators (e.g. diffusion equation) with the original form of the
hypersurface having the meaning of the initial concentration (or temperature) distribution.
However, the position of the minimum in this deformed function may have changed
during deformation and is usually different from that of the global minimum in the
original function, therefore a reversing deformation procedure is used to trace back from
the found minimum in the final deformed function to the related minimum in the original
function which is in most cases the true global minimum. An extremely important
characteristic of this method is the abilities of the deformation of the hypersurface without
investigation of large quantity of local minima or any information about their positions
beforehand. A number of smoothing procedures have been developed and applied to
peptides and proteins 386 .
• Molecular dynamics 387,388,389 is a stochastic method that simulates the physical behaviour
of molecules in a thermal bath. Molecular dynamic procedures search the conformational
superlattice of energy minima using an MD-like sampling strategy. It exploits the
additional information about the conformational space using the MD velocities. They
reflect, to varying degrees, all components of the system together that increases the
efficiency of conformational sampling in more complicated systems. Molecular dynamics
investigate the conformational surface incrementally in three-dimensional space.
Multidimensional conformational space is described in terms of a reduced representation
of conformational states. The algorithm initiates its move from a conformation state that is
a local minimum-energy configuration by stepping along the MD velocities. The atoms in
the molecule are usually constrained using a force field defined as a function of atom type,
bond type, torsion type and interatomic distance. Typically for this technique is the
application of Newton’s second law of movement to all atomic degrees of freedom under
which the new positions and velocities of the atoms are calculated. The atoms are moved
to these new positions and the cycle is repeated. Thus, the dynamic behaviour of the
molecule at the desired temperature can be reproduced by performing this process for
some time space. Conformations are selected from the trajectory and minimized at each
time step. An elevated temperature is used in order to overcome potential energy barriers

-99-
 and to reach new regions of the PES that may contain better lower-energy minima
conformers than the current region. Since the Molecular dynamics require a time step of
the order of 1 fs, it is computationally expensive and is often limited to
peptides358,360,361,362,363,364,365,357,390 .
• Corner flapping 391 produces a new conformation that is then subjected to energy
minimization. It is an intuitional flap movement of ring atoms of a cyclic molecule
relative to the average ring plane, ring atoms occupying the "corners" of the bonds. It is
preferable to define a local plane based on several ring atoms near the corner and to move
the corner atom vertically to the mirrored position on the other side of the plane (Figure
3-1). The new conformation is optimized and compared with the saved conformers relying
on usual criteria such as symmetry considerations, dihedral angle pattern, and steric
energy. Only new conformers are saved while repeated conformers are rejected. When all
corners have been flapped, this molecular structure is marked as “investigated”. After that,
the most stable of the “uninvestigated” saved conformers take part in the subsequent
flapping/minimization process. Since only local deformation of the cyclic molecule takes
place the subsequent energy minimization process is rapid. It is significant that in contrast
to other systematic algorithms this approach does not produce a tree structure, since every
branch is instantly pruned to retain the new conformer only. The search terminates when
all the saved conformers are “investigated”. Corner flapping sometimes leads to
excessively strained structure when the flap movement shifts the corner atom into the ring.
These weaknesses result since the environment of the corners being flipped are not into
account.

a b c

guache → -gauche guache → anti

Figure 3-1. Illustration of the a) Corner flapping, b) and c) Edge flipping.

• Edge flipping 392 is a modification of Corner flapping (Figure 3-1). This new local
perturbation mode flaps a pair of neighbouring corners together in opposite directions so
that the bond between the two corners changes its rotation. There are six change
possibilities, which are illustrated below:

-100-
 A

G G/

The gauche bond (G) changes its sign (G/) or becomes anti (A), while an anti bond (A)
transforms into the gauche bond with both signs (G and G/). In action, to prevent the
forming of too highly strained structures the movements of corner atoms in a ring
structure are carried out according to the rotation patterns of the two neighbouring bonds.
• Torsion flexing 393 is a term of a new stochastic (Monte Carlo) procedure that has been
developed for the conformational search of cyclic molecules. It produces a local distortion
of a ring bond in a cyclic molecule. Due to the fact that this approach does not evoke large
atomic movements, following energy minimization processes are usually carried out
rapidly even in cases where several bonds are distorted at the same time. The torsional
flexing method is able to overcome potential energy barriers so that the subsequent energy
minimization often yields another conformer of the cyclic molecule. It is significant that
such structural perturbations are in most cases not arduous that avoids rapid minimization
process. Figure 3-2 illustrates one example for torsional flexing. It is a local torsional
rotation about one ring bond which does not change the atomic position of most of the
ring atoms. To describe torsional flexing a torsion angle in a ring X−A-B-C-D−Y are
considered in which the ring bond X−A, as well as D−Y are imagined to have been
cleaved, resulting in two fragments. Both rotations are carried out by random angles. If
the B-C bond is a fusion bond between fused rings both rings and sets of atoms participate
in rotation process then a rotation around the B-C bond is performed. A (and its
substituents) and B (and its substituents) are rotated clockwise while C (and its
substituents) and D (and its substituents) are rotated anticlockwise.

Figure 3-2. Torsional flexing, where the arrows show the direction of torsional rotation. Note that the right side
of the resulting molecule structure is remained without changes after the torsional flexing procedure.

-101-
 The structure is generated by reclosure of the ring. The resulting structure undergoes an
energy minimization. During the single torsional flexing operation seven torsion angles
are changed and the related bond lengths and bond angles are also modified. The number
of torsion angles of torsional flexing operation is increased in the case of fused rings.
Since bond length and angles may change during the energy minimization the torsional
flexing is not completely a torsional coordinate method. However, the dominant structural
modifications of a molecule take place in the torsional space. This approach is an efficient
procedure for locating most of the low-energy conformers of a variety of cyclic
molecules. Moreover, it is easy to implement and can easily be combined with existing
molecular mechanics programs. Using of more than one conformational search procedure
increases the guarantee of a complete search of the conformational space available to
cyclic molecules. Nevertheless, torsional flexing is a method of choice for conformational
search for cyclic molecules.
• Low-Mode (LMOD) methods309,394,395 are based on the eigenvector-following techniques,
which were developed for locating saddle points on molecular potential energy surfaces.
The saddle point search is initiated at or near a local minimum. The eigenvalues and the
associated eigenvectors called the “normal modes” of the vibrations are determined by the
diagonalization of the Hessian matrix. At the minimum the eigenvectors with the smallest
eigenvalues point into directions of the potential energy surface with the modest ascent.
Assuming that the lowest transition state lies in these directions the basic eigenvector-
following (mode-following) idea is to select one or multiple low modes of the starting
structure and follow the corresponding Hessian eigenvectors uphill till the energy barrier
is crossed. The crossing point is reached if the next step is slightly lower in energy than
the preceding one. A subsequent energy minimization will usually lead to the bottom of
the neighbouring minimum-energy well. The LMOD procedure starts with an initial
molecular model, which is energy minimized. The minimized structure is then subjected
to normal mode analysis described above, and a user-specified number of low-frequency
modes are stored. LMOD moves in both directions of the Hessian eigenvectors
corresponding to the stored low-frequency modes until the algorithm finds a barrier
crossing defined by the following criteria:
1. the end-structure of the trajectory is lower than the energy of the starting structure,
or
2. the structure is at least lower than it was in the previous iteration step and the
molecule has also moved further away from the starting structure.

-102-
 The resulting perturbed initial structure is minimized afterwards. Naturally, it is possible
to reach the minimum that is not connected to the initial minimum of the mode-following
procedure as a result of energy minimization as well, but substantially, LMOD
concentrates the search to the local neighbourhood of a minimum on the potential energy
surface. Due to the fact that LMOD always follows the direction computed at the
minimum, no re-evaluation of the Hessian has to be performed. A set of conformers,
which are used as starting structures for structural perturbation along their low-frequency
modes, is accumulated during the search process. In spite of the fact that the new minima
found during an LMOD search allow to explore the complete potential energy surface
(with the exception of the case of the high-energy conformers search that necessitates the
higher frequency modes) LMOD is a systematic search procedure. But it is confined by
the number of low-frequency modes considered. LMOD must be switched to a stochastic
or Monte Carlo procedure when all possible systematic search directions are already
explored for a particular minimum. LMOD was developed originally for macromolecules
but it can be also successfully applied to flexible docking and to any kind of molecular
systems including complexes. It carries out the search in conformation space of cyclic and
acyclic molecules equally efficient. This approach runs neither in the torsion nor in the
Cartesian space, nor in any other user-defined search space since it generates its own
search low-dimensional search space, which is spanned by the low-frequency mode
eigenvectors of the Hessian matrix. Moreover, there is no need for special handling of
rings. The advantage of conformational searches performed with LMOD is that it
generates structures automatically. Only the user-defined threshold for the low-frequency
modes and the energy threshold for energy minimization have to be provided.
• MOLS conformational search technique 396 is based on the technique of using Mutually
Orthogonal Latin Squares (MOLS) to perform enhanced sampling of the conformational
space. The main idea of this conformational search technique consists of putting all
conformational variables of the molecule to a specific set of values and calculating the
potential energy. Obtained values of the potential energy correspond to conformations
systematically chosen to sample the entire conformational space. The torsion angles are
usually chosen as conformational variables, but other variables are also possible. Each
variable in the molecule is capable of taking up some different values in a range. Although
this range could be restricted by various factors, in the case of torsion angle space it is
taken to be 0–360°. Each cycle of the MOLS procedure consists of four steps as it is
shown in Figure 3-3. The set of all possible values of the conformational variables define
the complete conformational space of the molecule. All combinations of these values take

-103-
part in the sampling, each combination specifying one conformation of the molecule.
MOLS is used to pick up such combinations. Each torsion angle is set to correspond to
one Latin square to identify the conformations at which the potential energy calculations
are to be carried out. It means that a set of MOLS forms by arrangement of the possible
values of each torsion angle in the form of a Latin square. The energy hypersurface at
these locations in torsion angle space is sampled on the second step. This is achieved by
calculating the potential energy at each of the points. Calculation of average potential
energy of the each value of the each parameter made on the third step is to recover the
energy map of the conformational space. For this purposes one-dimensional
representations of the variation of the potential along each of the torsion angles are
produced. To estimate the efficacy of setting a particular torsion to a specific value,
irrelative of the values of the other torsion angles, the average of potential over points in
the MOLS is calculated, where a Boltzmann weighting function is used to identify the
optimum value for each variable, and arrive directly at the best structure. The final fourth
step is an examination of each one-dimensional representation. Boltzmann weighted
averages of different sets of variables are taken to identify the optimum value for each
parameter corresponding to the optimal conformation. This completes one cycle of
calculations during which one low energy structure is identified. The structure obtained in
one cycle of MOLS is either the global minimum energy structure, or one of the various
low energy structures of the molecule. To locate another low energy structure the
procedure is repeated using a different set of MOLS, however, it is not improbable that
the search repeatedly achieves the previously obtained low energy structures. If the
solution structure does not improve after a number of iterations it supposes that all the low
energy structures of the molecule are identified. The MOLS method outlined above is thus
a systematic sampling of the variable space to identify a library of low-energy three-
dimensional structures. The search is unconstrained and easily parallelized. This method
can be successfully applied to search low energy structures for a variety of small peptides
at negligible computational time.

-104-
 Search space
parameterization

Initialization of the set of

MOLS to globally sample
the space

Potential energy
calculation at each sub
square of the set of MOLS

Calculation of average
potential energy of each
value of each parameter
Yes

Selection of the optimum

Another low energy value for each parameter
conformation is identifying the set of
obtained? values corresponding to the
optimal conformation

No

Exit

Figure 3-3. Flowchart of the MOLS procedure.

• Ant Algorithm70 The original ant system, which was already described in the first part of
the work, was developed for optimization in the continuous search space of
conformational analysis. The conformational analysis problem is formulated as a
multimodal function minimization problem in n dimensions, where n is the number of
rotatable bonds in the flexible molecule. The aim is to find the global minimum of the
function e=f(θ1,θ2,…, θn), where e is the energy of the molecule and θ1,θ2,…, θn are the
torsion angles of the rotatable bonds. A roulette wheel selection algorithm is used to
generate the torsion angles. The name of this selection method arises since a roulette
wheel with slots sized according to the value of the probability function p(θi,t)=(1-

-105-
 ß)τ(θi,t)+ßη(θi) is used. p(θi,t) is associated with each torsion angle interval that gives the
probability that an ant will choose the value for the ith torsion angle in iteration t of the
algorithm. ß is a parameter between 0 and 1. τ(θi,t) corresponds to the pheromone trail and
is updated in each iteration of the algorithm. η(θi) corresponds to the visibility and keeps
constant. Artificial ants move at n intervals between 0° and 360°, and the path each ant
follows in one iteration of the algorithm determines the configuration of the molecule. The
pheromone trail is updated in each iteration of the algorithm. At each iteration all ants
provide data about conformations of the molecule, and for each ant the corresponding
strain energy is calculated using the molecular mechanics. It was tuned and tested for the
conformational analysis of flexible drug-like molecules, with the objective of obtaining an
accurate estimate of the lowest possible conformational energy.

3.2 GOTS application

On the basis of our tests the GOTS turned out to be the most efficient of the
developed algorithms. Hence, it was adapted for conformation searches. For this purpose
some modifications of the algorithm and a smart definition of the variables used for the
GOTS were necessary. These adaptations and first applications are described in the following.

3.2.1 Variable choice and numbering rules

The choice of the coordinates used to describe a molecule is critical for the efficiency
of the optimization. Since gradients and Hessians are usually calculated in Cartesian space,
Cartesian coordinates would be the most straightforward choice. But without the initial
curvature information, optimization in Cartesian coordinates is extremely inefficient,
especially for large systems 397 . Additionally, many studies have shown that Cartesian
coordinates are generally outperformed by a well-chosen set of so-called “natural internal”
coordinates 398 - consisting of bond lengths, bond angles, and torsions. Internal coordinate
systems can be classified as redundant or nonredundant. In a nonredundant internal
coordinates system, the number of coordinates is equal to the number of internal degrees of
freedom 3N-6 (where N is the number of atoms in the system); in a redundant set the number
of coordinates exceeds the number of internal degrees of freedom. So internal coordinates are
frequently used to represent conformations of molecules. These coordinates, first introduced
into geometry optimization by Pulay et al. 399 , involve the use of individual bond
displacements as stretching coordinates and linear combinations of bond angles and torsions
as deformational coordinates.

-106-
 The Z-matrix coordinates are a typical example of nonredundant internal
coordinates 400 . They describe the molecular geometry in terms of individual distances, angles,
and torsions. Each line of a Z-matrix gives the internal coordinates for one of the atoms
within the molecule. The Z-matrix format used in the Gaussian program package uses the
following syntax: Element label, atom 1, bond length, atom 2, bond angle, atom 3, and
dihedral angle. As an example the Z-matrix of n-Butanol is shown in Table 14. It provides a
description of each atom in a molecule in terms of its atomic number, bond length, bond
angle, and dihedral angle.

Table 14. A part of the internal coordinates in Z-matrix form of the n-Butanol molecule.

Element Bond Bond Dihedral

Atom 1 Atom 2 Atom 3
label length angle angle
c
… … …
c 3 cc4 2 ccc4 1 dih4
… … … … … … …
h 2 hc9 3 hcc9 4 dih9
h 2 hc10 3 hcc10 4 dih10
… … … … … … …

The illustration of the n-Butanol molecule is presented on Figure 3-4.

Figure 3-4. Atom numeration of the n-Butanol molecule.

In difference to other Z-matrix-like coordinates normally used in the

literature 401,402,403,404,405 , instead of altering of a bond length, a bond angle, and a dihedral
angle of each atom, only dihedral angles are used in the present work since the bond length
and angles of different conformers vary relatively little. Thus, they are defined to be stiff. If
Z-matrix coordinates are used for optimizations, in particular for ring systems, it is difficult to
find which an appropriate set of internal coordinates for the optimization process. Problems

-107-
occurs that various dihedral angles have to be changed together to perform reasonable
rotation. To rotate the terminal methyl group consisting of atoms C1H6H7H8 around the C1C2
axis it is necessary to modify three dihedral angles (C3C2C1H6, C3C2C1H7, C3C2C1H8)
although there is only one rotation. A crucial point is that these torsion angles must not be
varied independently of each other. To demonstrate it on the example lets take a typical
conformation of the molecule and attempt to move only one of dihedral angles at a time while
keeping the rest of the variables constant. Moving C3C2C1H6 only while keeping C3C2C1H7,
C3C2C1H8 constant distorts the internal structure of the methyl group. It is necessary to move
various dihedral angles synchronously to achieve proper rotations around a freely rotatable
single bond. The approximate separation of hard and soft movements of the system allows to
create the structure of the molecule by setting the hard coordinate to constant or to particular
dependence on the soft coordinates 406,407,408 . For this one partition the angles in hard and soft
dihedral angles is carried out. The soft ones are varied while the hard ones are changed
accordingly or keep constant. The problem was recognized and the concept of “related
dihedrals” was introduced by G. A. Chass 409 .
In this work the rules proposed by Echenique and Alonso 410 are used to account for
this problem. To perform a rotation around a given bond only the “main torsions” are
independently varied. The variations of the “dependent torsions” are not free but make sure
that proper rotations take place. In Table 15 the Z-matrix from Table 14 is reorganized
relative to the angle (4,3,2,1). This angle is a “main torsion” while the other two are
“dependent torsions”.

Table 15. Information from Z-matrix of the n-Butanol molecule.

Atom number Atom name Bond length Bond angle Dihedral angle
… … … … …
4 C (4,3) (4,3,2) (4,3,2,1)
… … … … …
9 H (9,2) (9,2,3) (9,2,3,4)
10 H (10,2) (10,2,3) (10,2,3,4)
… … … … …

To ensure a proper behavior the following scheme is used:

1. Write all dihedral angles of the Z-matrix into the list called “torsions”. Empty the lists:
“main torsions” and “dependent torsions”.
2. The first dihedral angle of the Z-matrix is chosen and written in the list “main torsions”.
3. Take the next dihedral angle from the “torsions” list.

-108-
4. Compare it with all previously found “main torsions”. If the second and third atoms of the
dihedral angle under consideration agree to those of a “main torsions” it is a “dependent
torsion”. It is related to the “main torsion”. If no agreement exists, it is a next “main
torsion”, and is written to the main torsions list.
5. If the “torsions” list is empty, the selection is finished. Otherwise, the process go back to
step 3.
It is significant that the vector direction does not count in this case (see Figure 3-5).

b2
b2 b3 d1 d3 c3 b2
b1 b c1

1 2 3
Figure 3-5. 1, 2 and 3 are related dihedral angles defined by three bond vectors connecting four atoms.

If the scheme is applied to the n- Butanol molecule (Table 14) Table 16 is obtained.

Table 16. “Main torsions” and “dependent torsions” of the n-Butanol molecule.

“Main torsions” “Dependent torsions”

(4,3,2,1) (9,2,3,4)
(10,2,3,4)
(11,3,2,1)
(12,3,2,1)
(5,4,3,2) (13,4,3,2)
(14,4,3,2)
(6,1,2,3) (7,1,2,3)
(8,1,2,3)
(15,5,4,3) -

Since only the “main torsions” are varied independently the number of degrees of freedom
decreases.

3.2.2 Adaptation of the GOTS to the conformational search

To adopt the GOTS approach for the conformation search some modification were
implemented. Its starts with the local minimization approach which is carried out as geometry
optimization of the computational chemistry environment ChemShell using the Universal
Force Field (UFF) 411 . To solve conformation search problems the variable choice is carried

-109-
out as indicated above. For the GOTS338 only function evaluations should be used to estimate
step size and directions. For the first step of the modest ascent all functional values Fz i+ and

Fz i− (Chapter 2 - Eq. 11) are computed using ChemShell single-point calculations to

determine the general direction of the following moves. They are obtained by varying each
dimension by a user defined step size Δxi0 . Step size and direction are estimated according to
Chapter 2 - Eq. 14.
As described in the second chapter a ranking procedure of the internal degrees of
freedom is used to find the optimal direction for the modest ascent. Applying this ranking
procedure to the conformational search it turned out that the computed energy differences are
too small. To obtain a reasonable differentiation between the various degrees, dynamical
calculations of the ranking coefficients rankmax and rankmin are necessary. It means that
rankmax and rankmin are not set to 0 and 1 but calculated for each case separately using the
following approach (Eq. 16-19):
Fz max ⋅ 100
coefficient = − 100 (16)
min Fz min

100 − coefficient
rank min = , if 0 < coefficient < 100
100 min (17)
rankmin = 0.1 (or 0.0 - depends on problem) , otherwise

rankmin = 0.1, value difference is large, therefore we have to redefine rankmax and after that
we use the formula (Chapter 2 - Eq. 14) to calculate ranki
rankmin ≠ 0.1, value difference is small, we directly use the formula (Chapter 2 - Eq. 14) to
calculate ranki

To calculate rankmax the same percentage analysis as for rankmin is made but for each value:
Fz i ⋅ 100
coefficient = − 100 (18)
max Fz min

100 − coefficient
rank max = , if 0 < coefficient < 100
100 max (19)
rankmax= 1.0 , otherwise

The resulting modest ascent is followed until the next calculated function value is smaller
than the previous one. This indicates that the barrier to the next valley is crossed. From this
point the next local minimum is located using local minimization.
-110-
 Other neighborhood strategies were also implemented and tested. In the first variant
(see Chapter 2 - Eqs. 11-13) not only the preferential direction determined at the minimum is
tested but at each step of this neighbourhood search all functional values Fz i+ and Fz i− are
computed to determine the direction of the next move (Chapter 2 - Eq. 11, Figure 3-6). This is
different to the standard approach in which the main promising moving area is determined at
the local optimum and kept fixed during the modest ascent strategy. The moves are carried
out with the step sizes for each variable that are estimated according to the linear ranking
procedure. It turned out that the selected moves now often lead to areas close to already
visited solutions. In consequence, the search starts to cycle and needs a considerable higher
number of function evaluations to come to the next valley. Such ideas were rejected since this
increases the computational costs without having many advantages.

3
2

-1

-1 0 1 2 3

Figure 3-6. Illustration of neighbourhood strategy with testing all direction at each step.

Ranking and the neighbourhood strategies were also considered. In one version the ranking
procedure was skipped, i.e. the move is carried out along the coordinate with the best function
value (see Figure 3-7b). However, also in this case more steps are needed to escape from a
valley. In another approach (Figure 3-7) a number of steps are tested for the selection of the
most promising move. This included steps with ranking and steps without ranking.
Furthermore, also small additional steps in positive and negative direction were performed
after the step without ranking. Comparing the functional values of all these steps, the one
from the step without ranking was often the best one. This means that the outcome of this
version is closed to the previous one. Based on test results the original neighbourhood
strategy turned out to possess the best cost-benefit ratio.

-111-
 2 2

1.5 1.5

1 1

0.5 0.5

0 0

-0.5 -0.5

-1 -1
-1 -0.5 0 0.5 1 1.5 2 -1 -0.5 0 0.5 1 1.5 2

a) step with ranking b) step without ranking

2

1.5

0.5

-0.5

-1
-1 -0.5 0 0.5 1 1.5 2

c) small steps in positive

and negative direction

Figure 3-7. Illustration of the neighbourhood search with various kinds of steps.

The TS is a mathematical optimization method, which belongs to the class of local

search techniques. To include more probability aspects, ideas from the Simulated Annealing
(SA) method have been implemented. It shall help to focus on promising areas and to locate a
good approximation to the global optimum of a given function in a large search space. In our
approach, the SA ideas are exploited to determine if a new minimum is taken or not. If the
new minimum is lower in energy it is always taken. If it is higher than the previous one the
criteria used in the SA (Chapter 1 – Eq. 6) are taken to decide if it is taken as the new starting
point or not. If the new minimum is not accepted the algorithm returns to the previous
minimum and continues the search along the next modest ascent direction.

-112-
 Initialisation

Local search for next

minimum

Simulated Annealing testing

accepted not accepted

Solution Update
Restore the previous found
minimum
Modest Ascent Search to
escape local minima

Diversification Search to
obtain new starting point

Figure 3-8. Flowchart of GOTS with SA elements.

To avoid reverse moves and cycles within the Tabu Search strategy a Tabu List (TL)
is built up which memorizes the moves previously made. The elements of the TL are ordered
according to the chronology of their appearance. If the TL is full the newest solution replaces
the oldest one employing the FIFO principle. For approaches which move in small steps
through the function F(xi) a simple storing of all visited points is sufficient to block already
visited regions effectively. To overcome the problems resulting from making large steps Tabu
Regions (TR) are used which were already successfully employed in the GTS and TSPA
methods and which are related to ideas introduced by F. Glover. Tabu Direction (TD)338
turned out to be less efficient for the conformational search as a result of the circular motion
of the the torsion angles from 0° to 360°. Examples of this are presented in Figure 3-9. The
start point is 135°. The step size is 90°. The first movements are shown in part a) at the Figure

3-9. After these two steps the TL values are 135°, 225°, 315°. Vectors 2,1 and 3,2 are TDs.
However, the next steps are already problematic if the TD approach is used. It is illustrated
with the parts b) and c). In the case of selecting 360 as the number of degrees the next
solution is 45° (see example b)). This unvisited solution is set “tabu” through the TD because

the motion vector 3,4 coincides with the TD vector 3,2 . Such approach leads to exclusion of

-113-
some unvisited solutions from the search. In the case of the unbounded number of degrees
(see example c)) all subsequent solutions (405°, 495°, 585°, 675°, and etc.) are allowed
because there are no such values in the TL. It results in cycling of the search process to the
rotation of a single torsion angle. So, the solution is “tabu” if all of the variables are in the
TR.

- start solution
- next solutions
135° 1 - problematic solutions
- Tabu Directions
0°

2 3
315°
225°

a)
135° 1 4 45° 135° 495° 405°
5 4

0° 0°

6 7
2 3 585° 675°
225° 315° 225° 315°

b) c)

Figure 3-9. Rotation of the torsion angle.

The last element of the TS represents the diversification search (DS) which is used to
select new starting points in the solution space. A DS becomes necessary if the solution does
not improve after a number of iterations or if all neighborhood solutions are already set tabu.
One way to enhance the convergence is to ensure that the search switches to regions that were
not already investigated. In the present approach the following strategy is used. There are N-3
(where N is the number of atoms in the system) “main and dependent torsions”. “Dependent
torsions” have to be changed together with the “main torsions”. To find a new starting point
the “main torsions” are changed with steps of 30°, 60° or 120°. The rest of the variables are
kept constant. To obtain reasonable geometries all “dependent torsions” are moved

-114-
accordingly. All points are excluded which already belong to Tabu Regions and the new
search starts from the point with the lowest F(xi) value.
Figure 3-8 gives the flowchart of the GOTS method with SA elements. At first a start
molecule structure in Z-matrix form is initialized. After that one starts with the search for
improved minima. This part consists of local optimizations that are applied to obtain nearest
local minima and the modest ascent searches used to escape to the next valley. The simulated
annealing criterion is used to confirm or reject a new local minimum as the new starting point.
In the case of rejection of a current minimum one returns to the previously found minimum
solution and starts the search along the second promising direction based on the already
calculated ranking coefficients. For the modest ascent search weighted function values can be
used. In between an update of the solution vector and the Tabu List is performed. If the
solution does not improve after a given number of iterations or if all neighborhood solutions
are already set tabu the search for improved minima is aborted and the diversification search
is performed to obtain new starting points.

3.2.3 Comparing the influence of different parameters

The proposed algorithm has several parameters. Some of them were already discussed
in chapter 2. For the application of the method to conformational searches some of them are
especially important for the fast convergence and efficiency of the optimization. The standard
values have also to be adapted to the conformational search. The most important parameters
are reported in Table 17.

Table 17. Description of the parameters of the approaches.

Parameter Purpose Recommended values

Δxi step size at the mildest ascent strategy 45°
number of not improved minima
BADMAX 5
after which a DS is performed
number of trial solutions for the diversification depend on number of
Ntrial
search “main torsions”
Diversification
step size used during the diversification strategy frequency 60°
step
rankmax default maximum recency ranked value 1.0
rankmin default minimum recency ranked value 0
TR ½ of the Tabu Region diameter 10°
>> ∆E normally
T control parameter (analog of temperature)
encountered

-115-
 The parameter Δxi gives the step size during the mildest ascent search. One could
assume that a small step size increases the accuracy, but this is misleading. For small step
sizes a wrong mildest ascent direction is taken. This is due to the fact that the energy
differences between the directions are so small that the ranking procedure does not work
anymore. Based on our experiences Δxi= 45° is recommended as a standard value. If even
larger step sizes are used fewer steps are necessary to leave a local minimum, but sometimes
minima are missed.
The BADMAX value was set to 5. This parameter also plays an important role during
the optimization process. Together with the Simulated Annealing approach the parameter
BADMAX controls when the search switches to the diversification search, i.e. it aborts the
search in the current region. If the solution does not improve after BADMAX local minima a
diversification search is started to select new starting points in the solution space.
The global time-varying parameter T, with which a higher lying minimum is taken as
the new starting point, influences the probability (see Eq. 5-6). T, called the temperature, is
gradually decreased during the optimization process. To determine this parameter
experimentation is required. The starting T value is set considerably larger than the largest ∆E
normally encountered and decreased by 10 percent with each new found minimum. The
dependency is such that a new local minimum solution is practically always accepted when T
is large even if it is quite high with respect to the last one. If T goes to zero the probability that
higher lying minima are accepted as new starting point decreases considerably.
During the diversification search steps equal to 60° were used. It means that a new
start solution is obtained by changing all “main torsions” by 60°, 120°, 180°, 240° or 300° in
random consecution. The “dependent torsions” are changed accordingly. The parameter Ntrial
is also connected with the diversification search. In our test cases Ntrial was set to the number
of “main torsions”.
The parameters rankmax and rankmin are ranking parameters from the linear ranking
procedure (Eqs. 14-19). They are calculated dynamically in the range from 0.0 to 1.0 that
means separately for each case.
The parameter TR controls the size of the tabu regions within the whole optimization
process. One would expect that small TR do not efficiently block already visited regions so
that the effort decreases with increasing TR. However, in most cases the effort increases for
larger TR. This counter-intuitive behaviour may result because the tabu regions become so
large that minima lying close by already visited points are overlooked or because the optimal
path to the global minimum is blocked. A value of TR = 10°, however, seems to be a good

-116-
choice for a large variety of different problems. So, if the new solution variables differ from
the variables of the already visited solution within ±10° it consider as a tabu solution.

3.3 Experimental results

The developed GOTS was applied to detailed conformational studies of Lysine (LYS),
Valine (VAL), Arginine (ARG), 2-acetoxy-N,N,N-trimethylethanaminium (Acetylcholine). In
addition two ACE inhibitors and a HIV-1 protease inhibitor were included. All these
molecules are presented in Table 18. LYS has 21 dihedral angles, ARG - 23, VAL – 16,
Acetylcholine – 23, the first ACE inhibitor - 43, the second ACE inhibitor called Fosinopril -
82 and HIV-1 protease inhibitor – 82. According to our definition of “main torsions” and
“dependent torsions” the number of variables decreases and there are only 8 “main torsions”
in the case of LYS molecule, 11 for ARG, 6 for VAL, 8 for Acetylcholine, 17 for the first
ACE inhibitor, 34 for Fosinopril and 41 for HIV-1 protease inhibitor to optimize.
As noted above, the program employs the computational chemistry environment
ChemShell under Linux for the local geometry optimization. To determine function values the
single-point calculations were carried out. The Universal Force Field (UFF)411, a molecular
mechanics force field, was used in the GOTS procedure. The parameters used to generate the
UFF include a set of hybridization dependent atomic bond radii, a set of hybridization angles,
van der Waals parameters, torsional and inversion barriers, and a set of effective nuclear
charges. The ability of UFF to correctly reproduce the experimental rotational barriers and
conformational energies of a variety of organic molecules was investigated by Casewit,
Colwell, and Rappé 412 . Good agreement with experiment was observed when UFF was
applied to the conformational equilibrium of compounds for which charge evidently does not
play an important role.

-117-
Table 18. List of the test molecules.

Lysine (LYS)

Valine (VAL)

Arginine (ARG)

Acetylcholine

Peptidomimetic HIV-1
protease inhibitor

Inhibitor of angiotensin
converting enzyme (ACE)

Fosinopril

The atomic numberings used in our test runs are shown in Figure 3-10.

-118-
LYS
ARG

VAL

Acetylcholine

Inhibitor of ACE

-119-
 Fosinopril

Peptidomimetic HIV-I protease inhibitor

Figure 3-10. Used numbering systems.

-120-
3.3.1 Conformational studies of amino acids
In the structure shown in the Figure 3-11, R represents a side chain specific to each
amino acid. The central carbon atom is a chiral carbon atom (except for glycine) to which the
two termini and the R-group are attached.

Figure 3-11. The general structure of a α-amino acid, with the amino group on the left and the carboxyl group
on the right.

Peptides, and their building blocks, the amino acids, possess many conformers each of them
differs in the location and type of bend in the backbone, and in the conformational variability
space in the side chains. The similarity of all amino acids, they differ only in the R-group, can
be used to simplify the conformational search for these molecules. This can be realized by the
division of the molecule structures into parts. One of the parts is common for all amino acids
while another one includes the rest R. The divisions for molecules LYS, ARG, and VAL are
indicated in Figure 3-12.

VAL
LYS ARG

Figure 3-12. Division used in the conformational search of amino acids.

-121-
 The conformational search is started for the common part. The best solution of this is
kept constant within the modest ascent part of the search. During the local search for the next
minimum all variables take part, i.e. local optimizations are carried out without any
restriction. Variable dihedral angles for the first and second part of the whole optimization
process and for the optimization without parting are demonstrated with the example of LYS
(see Figure 3-13, where all unmarked torsions keep constant during modest ascent search).

a b c

Figure 3-13. Variable torsions a) of the first part b) of the second part c) employed in local search.

To proof the stability and efficiency of the GOTS method for conformational search tasks,
test runs were made from four different starting structures (with the exception of VAL
because of its size). Dihedral angles were modified to obtain the various start structures of the
molecules. The results are given in Table 19. It convincingly demonstrates that partings of
amino acids increase the efficiency of the optimization search. Additionally, better optimal
structures are obtained since less number of variables participates in modest ascent strategy.

Table 19. Test results for LYS, ARG, and VAL.

Best minimum
Start First minimum Best minimum without
Molecule
structure (Hartree) (Hartree) constants
(Hartree)

1 0.0170 0.0120 0.0121

LYS 2 0.0170 0.0120 0.0123

3 0.0148 0.0120 0.0120
4 0.0151 0.0120 0.0120
1 0.0391 0.0168 0.0179
2 0.0353 0.0170 0.0192
ARG
3 0.0336 0.0168 0.0189
4 0.0355 0.0175 0.0203
VAL 1 0.0153 0.0130 0.0131

-122-
Figure 3-14 - Figure 3-16 superimpose the first, best, and worst found energy minimum
structures for the LYS, the ARG, and the VAL.

Start structure

Worst minimum
First minimum

Best minimum

Figure 3-14. First, best, and worst minimum structures of the LYS.

To demonstrate the behaviour of the search optimization process Figure 3-17 monitors the
local minima found along the search starting from the first test structure. The minimum
ordinal number is given on the abscissa while the corresponding energy value is given on the
ordinate.

-123-
 Start structure

Worst minimum

First minimum

Best minimum

Figure 3-15. First, best, and worst minimum structures of the ARG.

-124-
 Start structure

First minimum Worst minimum

Best minimum

Figure 3-16. First, best, and worst minimum structures of the VAL.

Obviously the course of the search depends strongly on the molecule under
consideration. For example, in the case of the VAL molecule, there are only two “main
torsions” in the second part of the optimization process. Hence, one frequently returns to
conformations with the similar energy values. However, it is important to note that equal
energy values can characterize different molecule structures. In such cases a comparison of all
optimized structures is necessary to decide if a new structure is found. The course found for
Valine nicely shows how our procedure can move to a complete different region (steps 27-35)
but returns to the most favourable one. The most simple case is the LYS molecule. However,
the curve for the ARG molecule structure shows that to reach the best minimum the DS and
SA were successfully applied to improve the search outcome.

-125-
 0,050
LYS
0,045

Energy (Hartree)
0,040
0,035
0,030
0,025
0,020
0,015
0,010
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41
Minimum number
0,045
ARG
0,040
Energy (Hartree)

0,035

0,030

0,025

0,020

0,015
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41
Minimum number
0,019 VAL
0,018
Energy (Hartree)

0,017
0,016
0,015
0,014
0,013
0,012
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41
Minimum number

Figure 3-17. Illustration of the optimization process for LYS, ARG, and VAL molecules.

3.3.2 Conformational studies of acetylcholine

The Acetylcholine molecule was treated analogous to the amino acids. The molecule
was optimized with and without a division of the molecule. The small part was optimized
first and again kept constant during the following modest ascent searches. The partitioning
used for the latter and the dihedral angles of both parts of Acetylcholine molecule are
presented on the Figure 3-18.

-126-
 a b c
Figure 3-18. Variable torsions a) of the first part b) of the second part c) employed in local search.

Table 20. Test results for Acetylcholine.

Best minimum
Start First minimum Best minimum without
Molecule
structure (Hartree) (Hartree) constants
(Hartree)
1 0.0478 0.0387 0.0387
Acetylcholine 2 0.0459 0.0387 0.0387
3 0.0460 0.0387 0.0387

Three different start structures were created to test the abilities of the method. In all cases, our
method converged to the same minimum. It is depicted in the
Figure 3-19 together with the first and worst found minima. The energies are collected at the
Table 20.
As for the LYS, ARG, and VAL molecules the curve on the Figure 3-20 shows the
behaviour of the search optimization with the minimum ordinal number on the absciss and
corresponding energy value on the ordinate. In the case of Acetylcholine molecule, the larger
part possesses six “main torsions” ( total number of “main torsions” is eight). The best
minimum value was found already during the first few iterations, but DS and SA strategies
were often applied to try to explore other areas of the search space.

-127-
 Start structure

Worst minimum

First minimum

Best minimum
Figure 3-19. Start structure, first, best, and worst minimum structures of Acethylcholine.

0,048 Acetylcholine
0,047
0,046
Energy (Hartree)

0,045
0,044
0,043
0,042
0,041
0,040
0,039
0,038
1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41
M inimum numbe r
Figure 3-20. Illustration of the optimization process for Acetylcholine molecule.

3.3.3 Conformational studies of ACE and HIV-1 protease inhibitors.

To test our approach for the optimizations of larger molecules conformational
searches were carried out for two angiotensin converting enzyme (ACE) inhibitors and HIV-1

-128-
protease inhibitor. ACE inhibitors are used primarily in treatment of congestive heart failure
and hypertension whereas protease inhibitors are applied to treat or prevent infection by
viruses, including HIV and Hepatitis C.

a b

Figure 3-21. Variable torsions. a)ACE inhibitor b) Fosinopril c) HIV-1 protease inhibitor.

0,2

0,18
Energy (Hartree)

0,16

0,14

0,12

ACE inhibitor
0,1
1 4 7 10 13 16 19 22 25 28 31 34 37 40
Minimum number

Figure 3-22. Illustration of the optimization process for ACE inhibitor molecule.

-129-
 For these molecules the parting approach was not employed. The Figure 3-21
illustrates the dihedral angles that are varied during the optimization process. The behaviour
of the optimization process for the ACE inhibitor is given in Figure 3-22. During this process
only the indicated dihedral angles were varied. The ordinal number of the minima is on the
abscissa and the corresponding energy value is on the ordinate. The results, collected in Table
21, show that the new approach can achieve low lying conformers. Figure 3-23 - Figure 3-25
show the first, best and worst found energy minimum structures for ACE and HIV-1 protease
inhibitors.

Start structure

Worst minimum

First minimum

Best minimum

Figure 3-23. Start structure, first, best, and worst minimum structures of ACE inhibitor.

Table 21. Test results for ACE and HIV-1 protease inhibitors.

Best minimum
First minimum
Molecule without constants
(kcal/mol)
(kcal/mol)

Peptidomimetic HIV-1 protease inhibitor 0.757 0.729

Inhibitor of ACE 0.124 0.113
Fosinopril 0.201 0.183

-130-
 Start structure

Worst minimum

First minimum

Best minimum
Figure 3-24. Start structure, first, best, and worst minimum structures of Fosinopril.

-131-
 Start structure
Worst minimum

Best minimum
First minimum

Figure 3-25. Start structure, first, best, and worst minimum structures of HIV-1 protease inhibitor.

-132-
3.4 Conclusions
In this chapter the global optimization routine Gradient only Tabu-Search (GOTS) was
adapted and applied to conformational search problems. It provides the same global minimum
in the tested molecules independently which starting structure was taken. Finally, the method
is quite general and can work with more sophisticated geometry optimization and single-point
calculation methods of ChemShell. The Z-matrix form is used for the definition of a
molecular structure and to manipulate the chemical structures. Implementation of Simulated
Annealing ideas gives a possibility to delocalize the search, to investigate more promising
areas and to locate a good approximation to the global optimum of a given function within a
large search space.
The GOTS algorithm was tested for several molecular structures. It can be also applied
to cyclic molecules under condition that some of torsions are kept constant to save a ring
structures. The test calculations indicate that the GOTS method is an efficient approach for
conformational searches.

-133-
Chapter 4 Summary

The optimization is a part of mathematics. The optimization aim is to find a discrete

mathematical object that maximizes or minimizes an objective arbitrary function specified by
the user. The nature of such mathematical objects and the search space are usually problem-
specific. The optimized system is typically only optimal in one application. Many activities
can be formulated as optimization problems. Because of the computational complexity of
these problems, exact optimization techniques of operations research like linear programming
or dynamic programming are for the most part unfeasible for large-scale problems. Therefore,
researchers, engineers, and scientists have used metaheuristic search techniques to find near
optimal, good-enough optimal or even the global solutions. Of course, the ultimate aim is to
apply metaheuristic search techniques like GAs, SA, and TS to solve real-world problems.
This work encompasses three parts. The first part provides a concise review of the
most prominent metaheuristic concepts currently available and gives essential preliminaries
together with definition of the combinatorial optimization problems. It substantiates the
choice of the investigation direction and basis idea of the developed methods.
In the second part the new nonlinear global optimization routines based on the TS
strategy are described. The new approaches are the Gradient Tabu Search (GTS), the Gradient
Only Tabu Search (GOTS), and the Tabu Search with Powell’s Algorithm (TSPA). All try to
determine the global minimum of a function by the steepest descent - mildest ascent strategy.
The steepest descent is followed to find the next local minimum. To escape a local minimum
the algorithm moves along the modest ascent till a next valley is reached. Then again the
steepest ascent strategy is used to determine the minimum of this valley and so on. GTS and
GOTS use the combination of the Steepest Descent and the Quasi-Newton methods for the
steepest descent part. To realize the modest ascent strategy the GTS evaluate the diagonal
elements of the Hessian while the GOTS scans the neighbourhood by function evaluations.
Also the TSPA uses the latter method for the modest ascent strategy. For the steepest descent
part, however, the Direction Set Method with Powell’s Algorithm is employed. To ensure an
efficient blocking of already visited regions despite the lower number of steps these
approaches introduce Tabu Regions and Tabu Directions as new elements for the Tabu List.
The new algorithms are compared with other approaches using various well-known test

-134-
functions with varying dimensionality. These test computations with up to 50 variables show
that the new approaches outperform most previous approaches on the basis of the number of
function evaluations. The surfaces of the test function are shown on the Figure 4-1.
Comparison with the differential evolution ansatz proves the efficiency of the present
approaches on the basis of wall clock timings especially for a high number of variables. The
influence of user-defined parameters on the efficiency of the new approaches is also
investigated. Moreover, the proposed approaches could be enhanced even further if more
sophisticated descent method and diversification strategies are applied instead.
Ackley function Rastrigin function Griewangk function

Levy function Hansen function Branin function

Rosenbrock function Zakharov function Goldstein-Price function

Figure 4-1. Test function.

In order to test the appropriateness of the GOTS method for molecular geometry
problems, the optimization of a function with a functional form similar to general potential
energy functions (see Eq. 20) was carried out.

-135-
 n
(−1) i
E = ∑ (1 + cos 3ω i ) + ,
i =1 10.60099896 − 4.141720682(cos ω i )
(20)
where ω i is the torsion angle ω i ( i +3) ,

n+3 is the number of atoms or beads in the given system

Illustration of the surface of this function is presented on the Figure 4-2. It is significant that
the global minimum of the function is known. The number of its local minima increases
exponentially with the size of the problem, which characterizes the principal complexity of
the minimization of the molecular potential energy functions.

4
3
2
z
1
0 00

50 50

y 100 100 x

150 150

Figure 4-2. Surface of the function with functional form similar to general PES function.

The GOTS was successfully applied to this the n-chain problem (n = 7, … , 20), i.e. the global
minimum was always found. It was much efficient than the Deterministic algorithm that uses
Branch and Bound scheme with techniques of interval analysis to provide the lower bounds.
Figure 4-3 demonstrates the efficiency of the GOTS method in comparison to the
Deterministic algorithm on the basis of number of function evaluation that introduced in the
diagram form for function of 7, 10, and 20 variables. The amount of function evaluations
required determining the global minima increases by a factor of about 11 if the dimension
increases from 10 to 20 whereas this factor for Deterministic algorithm is 352.

-136-
 500 Deterministic
450 algorithm
400 GOTS
180000
350
160000
300 140000
250 120000

200 100000
80000
150 60000
100 40000
50 20000
2038
0
0 20
7 10
Figure 4-3. Diagrams for function of 7, 10, and 20 variables.

GOTS provides excellent solutions for hard problems with minimal computational
effort. The algorithm can be suitably modified to incorporate more rigorous mechanisms to
handle very large-scale problems. Because of problem-specific knowledge GOTS deals with
variables and constraints effectively. However, the trade-off between effectiveness and
generality requires careful consideration, and the junction with other techniques is necessary
in order to built even more general methodology. In this connection the recommended default
parameters were given at this work to be destined with several test cases. The large-scale test
cases show that the initial parameter settings for the GOTS may be overly cautious; however,
in most cases, the global minimum is found independently which parameter values are used.
Nevertheless, these calculations also show that parameters settings are problem dependent.
The optimal parameter values may require some adjustment to starting values indicated in this
work. Strategies to estimate optimal parameter values are discussed.
An efficient search for the global minimum of a highly dimensional function with
many local minima is central for the solution of many problems in computational chemistry.
Well known examples for such global optimization problems are the conformational search
for molecules with a high number of freely rotatable bonds or the optimization of parameters
of a force field. Conformational analysis is the study of a molecule conformations and the
influence of them on the molecule properties. The basis of a conformational analysis is the
conformational search which has the goal to identify the energetically lowest lying
conformers of the molecule. This usually requires the location of the global energy minimum
of the PES. These structures are very important because they determine most of the properties
of the molecule. Over the past several years, a multitude of conformational search techniques
have been developed for this purpose.

-137-
 In the last part of the work the GOTS is applied for such chemical optimization
problems. The chapter provides a systematic approach how the variables are chosen and the
adjustable parameters are set. As test cases the global minimum energy conformation of some
amino acids, of two angiotensin converting enzyme (ACE) inhibitors, of 2-acetoxy-N,N,N-
trimethylethanaminium, and of a HIV-1 protease inhibitor is determined.
The major advantages of the exploitation of GOTS method are that it provides the
achievement of the global minimum without depending on the starting structure. The concept
of main and dependent torsions was used to allow easy implementation of the internal
rotations. It decreases the number of the variables of the optimization process and gives the
possibility to speed up the optimization. The mildest ascent strategy does not require the
computation of gradients, thereby leading to an additional saving in computational time.
Finally, the method is quite general and can work with more sophisticated geometry
optimization and single-point calculation methods of ChemShell. Some alterations were made
at the original version of the method. The diversification strategy was modified. In the
original GOTS the diversification strategy becomes necessary if the solution does not improve
after a number of iterations that means that three found local minima at hand are worse than
the current solution or if all neighbourhood solutions are already set tabu. Implementation of
Simulated Annealing ideas concerning found local minima gives a possibility to delocalize
the search, to investigate more promising areas, and to locate a good approximation to the
global optimum of a given function in a large search space. If the new minimum is not
accepted according to the so-called Metropolis criterion the algorithm returns to the previous
minimum and continues the search along the next modest ascent direction. This increases the
efficiency of the GOTS method.
The GOTS algorithm can be applied to cyclic molecules under condition that some of
torsions keep constant to save a wholeness of a structure. So, the test calculations indicate that
the GOTS is an efficient approach for conformational searches but also provides widespread
fields of action.

-138-
Chapter 5 Zusammenfassung

Die Optimierung ist ein Teilgebiet der Mathematik. Das Optimierungsziel ist das
globale Maximum oder Minimum einer vom Benutzer bestimmte Zielfunktion. Der Suchraum
und die Natur solcher Objekte sind problemspezifisch und die mathematische Behandlung
meist recht kompliziert. In der Regel sind zur Lösung entwickelte methematische Modelle
zudem nur begrenzt verwendbar, d.h. schon kleine Veränderungen in den realen
Fragestellungen erfordern recht unterschiedliche Modelle. Viele Fragestellungen können als
Optimierungsprobleme formuliert werden. Wegen des hohen Rechenaufwandes sind die
genauen Optimierungsmethoden der Operationsforschung, wie z. B. lineare oder dynamische
Programmierung, zum größten Teil jedoch undurchführbar. Deshalb wurden Metaheuristik-
Suchmethoden entwickelt, um möglichst gute Näherungslösungen zu finden. Die
Metaheuristik-Suchmethoden wie genetische Algorithmen, das Simulated Annealing oder die
Tabu-Suche haben sich zu wichtigen Hilfsmitteln in den Naturwissenschaften und der
Wirtschafsmathematik entwickelt.
Die Arbeit umfasst drei Kapitel. Das erste Kapitel stellt eine kurze Zusammenfassung
über die bekanntesten, zurzeit verwendeten Metaheuristischen-Konzepte dar und gibt
notwendige Einleitungen zusammen mit der Definition der kombinatorischen
Optimierungsprobleme. Das Kapitel begründet die Wahl des Tabu-Ansatzes und diskutiert die
Basisideen der entwickelten Methoden.
Im zweiten Kapitel werden die neuen entwickelten, nichtlinearen
Optimierungsroutinen beschrieben, die auf Tabu-Suchstrategien beruhen. Die neuen
Algorithmen sind Gradient Tabu Search (GTS), Gradient Only Tabu Search (GOTS) und
Tabu Search with Powell’s Algorithm (TSPA). Bei diesen Methoden wird das globale
Minimum einer Funktion durch die „steilsten Abstieg - schwächste Aufstieg“ - Strategie
bestimmt. Dem steilste Abstieg folgt man, um den nächsten lokalen Minimum zu finden. Um
von einem lokalen Minimum zum nächsten zu kommen bewegt sich der Algorithmus entlang
der geringsten Aufstiegs bis das nächste Tal erreicht ist. Dann wird wieder die steilste Abstieg
Strategie verwendet, um das Minimum dieses Tales zu finden und so weiter. GTS und GOTS
verwenden der Kombination die Steepest-Descent- und die Quasi-Newton-Methoden für den
steilste-Abstieg-Bereich. Der GTS wertet die diagonalen Elemente des Hessians zur

-139-
Realisierung der schwächste-Aufstieg-Strategie aus während die GOTS Methode die
Nachbarschaft mit Funktionswertberechnungen abscannt. Auch der TSPA Ansatz verwendet
die letztere Methode für die schwächste Aufstieg Strategie. Für den steilsten Abstieg wird aber
die Direction Set Method with Powell’s Algorithm verwendet. Um eine effiziente Sperrung
der bereits besuchten Regionen trotz der geringeren Anzahl von Schritte zu erreichen, führen
alle Konzepte Tabu-Regionen und Tabu-Richtungen als Elemente der Tabu-Liste ein. Sie
blockieren bereits besuchte Regionen sehr effizient und sind flexibel und leicht zu verwenden.
Anschließend folgt eine Diskussion über Effizienz der Algorithmen. Hierbei werden die
Ergebnisse mit anderen Ansätzen verglichen. Der Vergleich wird an Hand gut bekannten
Testfunktionen mit unterschiedlichen Dimensionalität durchgeführt. Untersucht werden
hierzu Funktionen, die von 2, 10, 30 oder 50 Variablen abhängen. Berechnungen mit häufig
verwendeten Test-Funktionen (Abbildung 5-1) zeigen, dass die neuen Methoden effizienter
sind als bisherige Ansätze. Insbesondere bei hochdimensionalen Problemen ergeben sich
erhebliche Vorteile. Ein Vergleich mit dem Differential-Evolution-Ansatz beweist ebenfalls
die Effizienz der neuen Methoden auf der Basis von „wall clock timings“. Ferner beinhaltet
das Kapitel eine Untersuchung des Einflusses benutzerdefinierter Parameter auf die Effizienz
der neuen Ansätze. Weiterentwicklung wären durch anspruchsvollere Methode für das
„steilsten Abstieg“ - Verfahren oder Diversifikationsstrategien möglich.

-140-
 Ackley function Rastrigin function Griewangk function

Levy function Hansen function Branin function

Rosenbrock function Zakharov function Goldstein-Price function

Abbildung 5-1. Test Funktionen.

Um die Eignung der GOTS-Methode für Problemen der Moleküloptimierung zu

prüfen, wurde auch die Optimierung einer Funktion durchgeführt, die den in der Chemie
auftretenden Potentialhyperflächen ähnelt (siehe Gl. 20).
n
(−1) i
E = ∑ (1 + cos 3ω i ) + ,
i =1 10.60099896 − 4.141720682(cos ω i )
(20)
wobei ω i der Verdrehungswinkel ω i ( i +3) ,

n+3 die Zahl von Atomen oder Teilchen im gegebenen System ist.

Abbildung 5-2 illustriert den Funktionsgraphen des Funktion von 2 Variablen. Die Zahl ihrer
lokalen Minima nimmt exponential mit der Größe des Problems zu, was der Komplexität der
Funktionen der molekularen potentiellen Energie entspricht.

-141-
 4
3
2
z
1
0 00

50 50

y 100 100 x

150 150

Abbildung 5-2. Oberfläche der Funktion mit funktionalen Form ähnlich allgemeinen potentielle Energie
Oberfläche Funktion.

In den Tests für dieses so genannte n-Kette-Problemen (n = 7,…, 20) findet die GOTS
Methode immer das globale Minimum, allerdings treten auch die Probleme, die bei der Suche
des globalen Minimums solcher Funktionen entstehen, klar hervor. Ein Vergleich mit dem
Deterministischen Algorithmus der die „Branch-and-Bound“-Schema mit „Intervall-
Analyse“-Methoden verwendet, um die niedrigsten Grenzen zu erreichen zeigt, dass die
GOTS Methode auch für diese Anwendung sehr effizient ist. Abbildung 5-3 vergleicht den
für beide Methoden notwendigen Aufwand auf Basis der zur Bestimmung des globalen
Minimums erforderliche Funktionswertberechnungen. Untersucht wurden Funktionen, die
von 7, 10 und 20 Variablen abhängen. In allen Fällen wurde das benannte globale Minimum
gefunden. Dabei steigt die Zahl der notwendigen Funktionswertberechnungen bei der GOTS
Methode um einen Faktor von etwa 11, wenn die Dimension von 10 auf 20 vergrößern wird.
Für den Deterministischen Algorithmus beträgt der Faktor 352.

-142-
 500 Deterministic
450 algorithm
400 GOTS
350 180000
160000
300
140000
250 120000
200 100000

150 80000
60000
100 40000
50 20000
2038
0 0
20
7 10

Abbildung 5-3. Anzahl der Funktionswertberechnungen, die zur Bestimmung des globalen Minimums der in Gl.
20 definierten Funktion notwendig waren.

Alle unsere Untersuchungen zeigen, dass der GOTS die Minima komplexer Probleme
sehr effizient findet. Allerdings wird die Effizienz von dem Parameter des Algorithmus
beeinflusst. Aufbauend auf unseren Tests werden daher Werte für diese Parameter
vorgeschlagen. Unabhängig von der Wahl der Parameter werden die globalen Optima aber in
den weitaus meisten Fällen gefunden. Die Tests zeigen, dass die optimalen Großen der
Parameter problemabhängig sind. Daher werden ebenfalls Wege diskutiert, wie man optimale
Werte abschätzen kann.
Eine effiziente Suche nach dem globalen Minimum einer hochdimensionalen Funktion
mit mehreren lokalen Minima ist auch von zentraler Bedeutung für viele Probleme der
Computerchemie. Bekannte Beispiele für solche globalen Optimierungsprobleme sind die
Konformationsanalyse bei Molekülen mit einer hohen Zahl frei drehbarer Bindungen oder die
Optimierung der Parameter eines Kraftfeldes. Die Grundlage einer Konformationsanalyse ist
die Suche nach dem energetisch tiefstliegenden Konformer des Moleküls. Mathematisch ist
die gleich bedeutend mit der Suche des globalen Energieminimums der Energiehyperfläche
des Moleküles. Diese Konformeren sind sehr wichtig, da sie die Eigenschaften des Moleküls
bestimmen. In den letzten Jahren ist eine Vielzahl von konformativen Suchmethoden
entwickelt worden.
Das letzte Kapitel der Arbeit beschreibt die Anwendung der GOTS Methode auf
dieses Problem. Diskutiert werden die Auswahl der Variablen und die Einstellung der
justierbaren Parameter. Die Effizienz der GOTS Methode wird an Hand einiger Aminosäuren,

-143-
zwei Angiotensin-Derivaten (ACE-Hemmer), des Acetylcholin und eines HIV-1-Protease-
Hemmstoff gezeigt.
Die wichtigsten Vorteile der Verwendung der GOTS Methode bestehen darin, dass das
Erreichen des globalen Minimums nicht von der Ausgangsstruktur abhängt. Das Konzept, von
Haupt- und abhängigen Torsionen wurde verwendet, um interne Rotation um
Einfachbindungen zu ermöglichen. Hierdurch wird ebenfalls die Anzahl Variablen reduziert,
was die Optimierung effizienter macht. Die „mildeste Aufstieg“-Strategie verlangt keine
Berechnung der Gradienten, was ebenfalls Zeit spart. Der neue Ansatz wurde mit dem
ChemShell-Programmpaket verknüpft. Zur Anpassung an die neue Aufgabe waren einige
Modifizierungen notwendig. Insbesondere die Diversifikationsstrategie wurde angepasst. Im
ursprünglichen GOTS wurde nach einem neuem Startpunkt gesucht (Diversifikation), wenn
nacheinander drei höher liegende lokale Minima gefunden wurden. Zur
Konformationsanalyse wurde eingebaut, dass man mit einer gewissen Wahrscheinlichkeit
zum tiefer liegenden Minimum zurückgeht, um dort die Suche nach neuen Minima
fortzuführen. Gesucht wird dann in die Richtung, die den zweitniedrigsten Anstieg besitzt.
Dies vergrößert die Effizienz der GOTS Methode, da vielversprechende Bereiche intensiver
untersucht werden.
Die beschriebene Methode kann auch für zyklische Moleküle verwenden werden.
Allerdings müssen einige der Torsionswinkel konstant gehalten werden, damit die
Ringstrukturen nicht aufgebrochen werden. Die Testberechnungen deuten darauf hin, dass der
GOTS nicht nur ein effizientes Konzept für Konformationsanalyse ist, sondern weit
verbreitete Anwendungsbereiche finden kann.

-144-
Appendix A Unconstrained Test Problems

A.1 Branin function (BR)

1 1
Definition: BR( x) = (2 x 2 + sin(4πx 2 ) − x1 ) 2 + ( x 2 − sin(2πx1 )) 2 .
20 2
Feasible search space: − 10.0 ≤ xi ≤ 10.0 , i = 1,2 .

Global minimum: x* = (0,0) ; BR( x* ) = 0 .

A.2 Goldstein and Price function (GP)

GP( x) = (1 + ( x1 + x2 + 1) 2 (19 − 14 x1 + 3 x12 − 14 x2 + 6 x1 x2 + 3x22 ))
Definition: .
(30 + (2 x1 − 3 x2 ) 2 (18 − 32 x1 + 12 x12 + 48 x2 − 36 x1 x2 + 27 x22 ))

Feasible search space: − 2.0 ≤ xi ≤ 2.0 , i = 1,2 .

Global minimum: x* = (0,−1) ; GP( x* ) = 3 x*=(0,-1); GP(x*)=3.

A.3 Hansen function (H)

5 5
Definition: H ( x) = ∑ i cos((i − 1) x1 + i )∑ j cos(( j + 1) x2 + j ) .
i =1 j =1

Feasible search space: − 10.0 ≤ xi ≤ 10.0 , i = 1,..., n .

⎧(−7.589,−7.708) ⎫
⎪(− 7.589,−1.425)⎪
⎪ ⎪
⎪(− 7.589, 4.858) ⎪
⎪ ⎪
⎪(−1.307,−7.708) ⎪
⎪ ⎪
Global minimum: x * = ⎨(−1.307,−1.425) ⎬ ; H ( x* ) = −176.542 .
⎪(−1.307,4.858) ⎪
⎪ ⎪
⎪(4.976,−7.708) ⎪
⎪(4.976,−1.425) ⎪
⎪ ⎪
⎪⎩(4.976,4.858) ⎪⎭

-145-
A.4 Ackley function (AKn)
1 1 1
∑i=1 xi2 ∑i=1 cos( 2πxi )
n n
−
Definition: AK n ( x) = 20 + e − 20e 5 n
−e n
.

Feasible search space: − 30.0 ≤ xi ≤ 30.0 , i = 1,..., n .

Global minimum: x* = (0,...,0) ; AK n ( x* ) = 0 .

A.5 Rastrigin function (Rn)

n
Definition: Rn = 10n + ∑ ( xi2 − 10 cos(2πxi )) .
i =1

Feasible search space: − 5.12 ≤ xi ≤ 5.12 , i = 1,..., n .

Global minimum: x* = (0,...,0) ; Rn ( x* ) = 0 .

A.6 Griewangk function (Gn)

n
xi2 n
x
Definition: Gn ( x) = ∑ − ∏ cos( i ) + 1 .
i =1 4000 i =1 i
Feasible search space: − 600.0 ≤ xi ≤ 600.0 , i = 1,..., n .

Global minimum: x* = (0,...,0) ; Gn ( x* ) = 0 .

A.7 Levy function (Ln)

n −1
Definition: Ln ( x) = sin 2 (3πx1 ) + ∑ ( xi − 1) 2 (1 + sin 2 (3πxi +1 )) + ( xn − 1)(1 + sin 2 (2πxn )) .
i =1

Feasible search space: − 10.0 ≤ xi ≤ 10.0 , i = 1,...,4 ,

− 5.0 ≤ xi ≤ 5.0 , i = 1,...,7 .

Global minimum: for i=4 x* = (1,1,1,−9.752) ; Ln ( x* ) = −21.502 ;

for i=7 x* = (1,...,1,−4.754) ; Ln ( x* ) = −11.504 .

A.8 Rosenbrock function (Rbn)

n/2
Definition: Rbn ( x) = ∑100( x 2i − x 22i −1 ) 2 + (1 − x 2i −1 ) 2
i =1

Feasible search space: − 10.0 ≤ xi ≤ 10.0 , i = 1,..., n

Global minimum: x * = (1,...,1) ; Rbn ( x * ) = 0

-146-
A.9 Zakharov function (Zn)
n n n
Definition: Z n ( x) = ∑ xi2 + (∑ 0.5 ⋅ i ⋅ xi ) 2 + (∑ 0.5 ⋅ i ⋅ xi ) 4
i =1 i =1 i =1

Feasible search space: − 5.0 ≤ xi ≤ 10.0 , i = 1,..., n

Global minimum: x * = (0,...,0) ; Z n ( x * ) = 0

A.10 DeJoung function (DJ)

Definition: DJ ( x) = x12 + x 22 + x32

Feasible search space: − 2.56 ≤ xi ≤ 5.12 , i = 1,2,3

Global minimum: x * = (0,...,0) ; DJ ( x * ) = 0

A.11 Trid function (Trn)

n n
Definition: Trn ( x) = ∑ ( xi − 1) 2 − ∑ xi xi −1
i =1 i=2

Feasible search space: − n 2 ≤ xi ≤ n 2 , i = 1,..., n

Global minimum: for i=5 Trn ( x * ) = −30

for i=10 xi = i ⋅ (11 − i ) ; Trn ( x * ) = −210

for i=20 Trn ( x * ) = −1520

-147-
Appendix B Notation

TSP travelling salesman problem

GO global optimization
ACE angiotensin converting enzyme
F(xi) continuous objective function
xi decision variable vector
xl lower bound
xr upper bound
D feasible set
ACO Ant colony optimization
PSO Particle swarm optimization
GBAS Graph-based Ant System
S-ACO Simulation based Ant Colony Optimization
ℜn multidimensional space
x(p) particle position in multidimensional space
ν ( p) velocity of the particle in multidimensional space
best(p) “particle best” particle position
best(N(p)) “local (neighbourhood) best” particle position
best(Xpop) “global best” particle position
QSAR quantitative structure-activity relationship
GA Genetic Algorithm
CGA Continuous Genetic Algorithm
RNA Ribonucleic Acid
SS Scatter Search
MA Memetic Algorithm
SA Simulated Annealing
T main control parameter in the cooling schedule (temperature)
∆F function difference
F ( posi +1 ) function value of the resulting solution

-148-
F ( posi ) function value of the current solution
ρ probability
JSS job shop scheduling problem
VNS Variable Neighbourhood Search
Ni set of neighbourhood structures
Ni(x) set of solutions in the ith neighbourhood of solution x
VND Variable Neighbourhood Descent
RVNS Reduced Variable Neighbourhood Search
CPU Central Processing Unit
VNDS Variable Neighbourhood Decomposition Search
SVNS Skewed Variable Neighbourhood Search
PVNS Parallel Variable Neighbourhood Search
GVNS General Variable Neighbourhood Search
GRASP Greedy Randomized Adaptive Search Procedure
TS Tabu Search
TL Tabu List
FIFO First In First Out
GTS Gradient Tabu Search
TD Tabu Direction
TR Tabu Region
TDV Tabu Direction Vector
NMDV New Move Direction Vector
RTR Radius of the Tabu Region
Rup upper bound for each variable
Rlow lower bound for each variable
coeff coefficient used to compute the radius RTR
Δxi step size at the mildest ascent strategy
α Tabu Direction coefficient
alam first step size at the line search
Ntrial number of trial points at the diversification search
LSIZE number of elements in tabu list
ranki recency ranked value
rankmax maximum recency ranked value
rankmin minimum recency ranked value

-149-
NDIM dimension
Itermain
Iterworst
loop termination numbers
Iterloc
IterMAS
i
amin visited minima
i
bstart starting points of the optimization from which the respective minimum
i
c new next trial point
BFGS Broyden Fletcher Goldfarb and Shanno
LMOD Low Mode search
DS Diversification Search
AKn n-dimensional Ackley function
BR Branin function
GP Goldstein-Price function
Gn n-dimensional Griewangk function
Rn n-dimensional Rastrigin function
H Hansen function
Ln n-dimensional Levy function
Rbn n-dimensional Rosenbrock function
DJ DeJoung function
Zn n-dimensional Zakharov function
Trn n-dimensional Trid function
TRUST Terminal Repeller Unconstrained Subenergy Tunneling
GOTS Gradient Only Tabu Search
TSPA Tabu Search with Powell’s Algorithm
Fz i+ function value of the position with increasing

Fz i− function value of the position with decreasing

Di direction for each variable
Fzi functional value needed for the linear ranking procedure
Fzmax maximal Fzi value
Fzmin minimal Fzi value
N number of parameters
ECTS Enhanced Continuous Tabu Search

-150-
 ESA Extended Simulated Annealing

CHA Correlation Height Analysis

CTSS Continuous Tabu Simplex Search

PES Potential Energy Surface

BF Boltzmann Factor

E Energy

RGN Random Generated Number

d number of possible values of each torsion angle

dN number of conformations to be generated and minimized

MOLS Mutually Orthogonal Latin Squares

ANN Artificial Neural Networks

UFF Universal Force Field

coefficient minimal function values percentage
min

coefficient maximal function values percentage

max

LYS Lysine

VAL Valine

ARG Arginine

ACE Angiotensin Converting Enzyme

TR ½ of the Tabu Region diameter

BADMAX number of local minima which did not improved the solution
ωi torsion angle

-151-
References

1
Neumaier A., Global optimization and constraint satisfaction., In: Bomze I., Emiris I., Neumaier A., Wolsey
L. (eds.), Proceedings of GICOLAG workshop (of the research project Global Optimization, Integrating
Convexity, Optimization, Logic Programming and Computational Algebraic Geometry), 2006.
2
Gendreau, M.; Laporte G.; Semet, F., A tabu search heuristic for the undirected selective travelling salesman
problem, European Journal of Operational Research, 1998, 106, pp. 539-545.
3
Sexton, R. S.; Alidaee, B.; Dorsey R. E.; Johnson, J. D., Global optimization for artificial neural networks: A
tabu search application, European Journal of Operational Research, 1998, 106, pp. 570-584.
4
Costamagna, E.; Fanni A.; Giacinto, G., A Tabu Search algorithm for the optimisation of telecommunication
networks, European Journal of Operational Research, 1998, 106, pp. 357-372.
5
Teh, Y. S.; Rangaiah, G. P., Tabu search for global optimization of continuous functions with application to
phase equilibrium calculations., Computers and Chemical Engineering, 2003, 27, pp. 1665-1679.
6
Lin, B. ; Miller, D. C., Tabu search algorithm for chemical process optimization, Comp. Chem. Engineering,
2004, 22, 2287-2306.
7
Horst R. and Pardalos P. M.(eds), Handbook of Global Optimization, Kluwer Academic Publishers, Boston,
MA, 1995.
8
Horst R., Pardalos P.M. and Thoai N.V., Introduction to Global Optimization, Second Edition, Kluwer
Academic Publishers, Boston, MA, 2000.
9
Pardalos P. M. and Resende M. G. C. (eds.), Handbook of Applied Optimization, Oxford University Press,
Oxford, 2002.
10
Pardalos P.M. and Romeijn H. E (eds.), Handbook of Global Optimization, Kluwer Academic Publishers,
Boston, MA, 2002.
11
Pardalos P. M., Romeijn H. E and Tuy H., Recent developments and trends in global optimization, J. Comput.
Appl. Math., 2000, 124, pp. 209–228.
12
Ribeiro C.C. and Hansen P. (eds.), Essays and Surveys in Metaheuristics, Kluwer Academic Publishers,
Boston, MA, 2002.
13
Voß S. Intelligent Search TH Darmstadt Habilitationsschrift, Darmstadt 1993.
14
Osman I. H. and Kelly J. P., Meta-Heuristics: Theory and Applications, Kluwer Academic Publishers,
Boston, MA, 1996.
15
Glover F. and Kochenberger G. (eds.), Handbook of Metaheuristics, Kluwer Academic Publishers, Boston,
MA, 2002.
16
LaValle, S. M.; Finn, P. W.; Kavraki, L. E.; Latombe, J.-C. Journal of Computational Chemistry, 2000, 21,
pp. 731-747.
17
Seonggu R.; Baek S.-G.; Bogman L.; Chihyo P.; Nakyen C.; Chang S.L.; Young C.S.; Hoil C.; Jong S.K.;
Heungsik Y.; Sung C.K.; Jong H.O. Bioorganic and Medicinal Chemistry Letters, 1998, 8, pp. 2423-2426.
18
Gavrilets S., Fitness Landscapes and the Origin of Species. Number MPB-41 in Monographs in Population
Biology. Princeton University Press, 2004.
19
Renze J., Weisstein E. W., Extreme Value Theorem, MathWorld, A Wolfram Web Resource.
http://mathworld.wolfram.com/ExtremeValueTheorem.html
20
Weise T., Global Optimization Algorithms: Theory and Application, 2008 (e-book)
21
Müller-Merbach H., Heuristic and their design: a survey., European Journal of Operational Research, 1981,
8, pp. 1-23.
22
Pearl J., Heuristic: Intelligent Search Techniques for Computer Problem Solving. Addison Wesley, Reading,
1984.
23
Zanakis, S. H., Evans, J. R., Vazacopoulos A. A., Heuristic methods and applications: a categorized survey.,
European Journal of Operational Research, 1989, 415, pp. 298-311.
24
Reeves C. R., Modern Heuristic Techniques for Combinatorial Problems. Blackwell, Oxford, 1993.
25
Polya G., How to solve it. Princeton University Press, Princeton, 1945.
26
Glover F., Future paths for integer programming and links to artificial intelligence, Computers and
Operations Research, 1986, 13, pp. 533–549.
27
Goldberg D. E., Genetic Algorithms in Search, Optimization & Machine Learning, Reading: Addison-
Wesley, 1989.

-152-
28
Aarts, E. H. L., Lenstra J.K. (eds.), Local Search in Combinatorial Optimization. Chichester: Wiley, 1997.
29
Corne D., Dorigo M., Glover F. (eds.), New Ideas in Optimization, London: McGraw-Hill, 1999.
30
Michalewicz Z., Genetic Algorithms+Data Structures=Evolution Programs. Third Edition, Berlin: Springer,
1996.
31
Laporte G., Osman I.H. (eds.), Metaheuristics in Combinatorial Optimization., Annals of Operations Research
63, 1996.
32
Voss S., Martello S., Osman I.H., Roucairol C. (eds.), Meta-Heuristics: Advances and Trends in Local Search
Paradigms for Optimization. Boston: Kluwer, 1999.
33
Gendreau M., Guertin F., Potvin J.-Y., Taillard E.D., Parallel Tabu Search for Real-Time Vehicle Routing
and Dispatching., Transportation Science, 1999, 33, pp. 381–390.
34
Gendreau M., Laporte G., Semet F., A Dynamic Model and Parallel Tabu Search Heuristic for Real-Time
Ambulance Relocation, Parallel Computing, 2001, 27, pp. 1641–1653.
35
Glover F., Laguna, M., Tabu Search Kluwer, Academic Publishers, MA, USA, 1997.
36
Taillard É. D., Gambardella L.M., Gendreau M., Potvin J.-Y., Adaptive memory programming: A unified
view of metaheuristics., European Journal of Operational Research, 2001, 135, pp. 1-16.
37
Hopfield J. J., Tank D. W., Neural Computation of Decisions in Optimization Problems, Biological
Cybernetics, 1985, 52, pp. 141-152.
38
Van den Bout D. E., Miller T. K., III Graph partitioning using annealed neural networks., IEEE Transactions
on Neural Networks, 1990, 1 , pp.192-203.
39
Peterson C., Anderson J., Neural Networks and NP-complete optimization problems: a performance study on
the graph bisection problem, Complex Systems, 1988, 2, pp.59-89.
40
Lee H.-M., Hsu C.–C., Neural network processing through energy minimization with learning ability to the
multiconstraint zero-one knapsack problem., Proceedings Tools for AI Conference, Fairfax, Virginia, 1990,
pp. 548-555.
41
Tagliarini G. A., Page E.W., Solving constraint satisfaction problems with neural networks, Proceedings
IEEE International Conference on Neural Networks, 1987, 3, pp.741-747.
42
Kennedy M., Chua L., Neural networks for linear and nonlinear programming., IEEE Transactions on
Circuits and Systems, 1988, 35, pp.554-562.
43
Looi C.–K., Neural network methods in combinatorial optimization, Computers and Operations Research,
1992, 19, pp. 191-208
44
Dorigo M., Optimization, Learning and Natural Algorithms, PhD thesis, Politecnico di Milano, Italy, 1992.
45
Dorigo, M., Maniezzo, V., Colorni, A., The Ant System: An autocatalytic optimization process, Technical
Report, Dept. of Electronics, Politecnico di Milano, Italy, 1991.
46
Dorigo M., Maniezzo V., Colorni A., The Ant System: Optimization by a colony of cooperating agents, IEEE
Trans. on Systems, Man, and Cybernetics, 1996, 26, pp. 1-13.
47
Colorni A., Dorigo M., Maniezzo V., Distributed Optimization by Ant Colonies, Proceedings of the First
European Conference on Artificial Life, Varela F.J., Bourgine P. (eds.), MIT Press, Cambridge, MA, 1992, pp.
134-142.
48
Denebourg J.L., Goss S., Collective patterns and decision-making., Ethology, Ecology & Evolution, 1989, 1,
pp. 295-311.
49
Goss S., Beckers R., Deneubourg J.L., Aron S., Pasteels J.M., How Trail Laying and Trail Following Can
Solve Foraging Problems for Ant Colonies., In Hughes R.N. (eds.), NATO ASI Series, Vol. G 20, Behavioural
Mechanisms of Food Selection, Springer-Verlag, Berlin, 1990.
50
Gambardella L.M., Taillard E.D., Dorigo M., Ant colonies for the QAP. Technical Report 4-97, IDSIA,
Lugano, Switzeland, 1997.
51
Maniezzo V., Colorni A., The ant system applied to the quadratic assignment problem., IEEE Trans.
Knowledge and Data Engineering, 1999, pp.769 – 778.
52
Dorigo M., Gambardella L. M., Ant colonies for the travelling salesman problem., BioSystems, 1997, 43, pp.
73-81.
53
Bullnheimer B., Hartl R. F., Strauss C., A new rank-based version of the ant system: a computational study.,
Central European Journal for Operations Research and Economics, 1999, 7, pp. 25-38.
54
Stützle T., Hoos H., The MAX-MIN ant system and local search for the travelling salesman problem., In:
Baeck T., Michalewicz Z., Yao, X. (eds.), IEEE International Conference on Evolutionary Computation and
Evolutionary Programming Conference, IEEE Press, 1997, pp.309-314.
55
Dorigo M., Gambardella L. M., Ant colony system: A cooperative learning approach to the travelling
salesman problem., IEEE Transactions on Evolutionary Computation, 1997, 1, pp. 53-66.
56
Gambardella L. M., Dorigo M., Ant-Q: A reinforcement learning to the travelling salesman problem., In:
Proceeding of the Twelfth International Conference on Machine Learning, Palo Alto, CA: Morgan Kaufmann,
1995, pp.252-260.
57
Gambardella L. M., Dorigo M., Solving symmetric and asymmetric TSPs by ant colonies., In: Proceeding of
the IEEE Conference on Evolutionary Computation, ICEC96, 1996, pp. 622-627.

-153-
58
Merkle D., Middendorf M., Schmeck H., Ant colony optimization for resource-constrained project
scheduling., In: Proceedings of the Genetic and Evolutionary Computation Conference (GECCO-2000), 2000,
pp 893–900.
59
Bonabeau E., Henaux F., Guerin S., Snyers D., Kuntz P., Therulaz G., Routing in telecommunication
networks with “Smart” ant-like agents., Workshop on Intelligent Agents for Telecommunication Applications,
Springer Verlag, 1998.
60
Schoonderwoerd R., Holland O., Bruten J., Rothkrantz L., Ant-based load balancing in telecommunications
networks., Adaptive Behaviour, 1996, 5, pp. 169-207.
61
White T., Pagurek B., Oppacher F., Connection management using adaptive mobile agents. In: Arabnia H.R.
(eds.), Proceedings of the International Conference on Parallel and Distributed Processing Techniques and
Applications, CSREA Press, 1998, pp. 802-809.
62
Di Caro G., Dorigo M., AntNet: Distributed stigmergetic control for communications networks., Journal of
artificial Intelligence Research, 1998, 9, pp.317-365
63
Saleh M., Ghani A. A. A., Adaptive Routing in Packet-Switched Networks Using Agents Updating Methods.
Malaysian Journal of Computer Science, 2003, 16, pp. 1-10.
64
Heusse M., Snyers D., Guerin S., Kuntz P., Adaptive agent-driven routing and load balancing in
communication networks, Adv. Complex Syst., 1998, 1, pp. 237-254.
65
Bullnheimer B., Hartl R.F., Strauss C., An improved Ant System algorithm for the Vehicle Routing Problem.,
Annals of Operations Research, 1999, 89, pp. 319-328.
66
Bullnheimer B., Strauss C., Tourenplanung mit dem Ant System., Technical Report 6, Institut für
Betriebwirtschaftlehre, Universität Wien, 1996.
67
Gambardella L. M., Dorigo M., HAS-SOP: An hybrid ant system for the sequential ordering problem.,
Technical Report 11-97, IDSIA, Lugano, CH, 1997.
68
Costa D., Hertz A., Ants can colour graph., Journal of Operational Research Society, 1997, 48, pp. 295-305.
69
Michel R., Middendorf M., An island model based ant system look ahead for the shortest supersequence
problem., In: Eiben A.E., Back T., Schoenauer M., Schwefel H.-P. (eds.), Proceedings of PPSN-V, Fifth
International Conference on Parallel Problem Solving from Nature, Springler–Verlag, 1998, pp.692-701.
70
Daeyaert F., De Jonge M., Koymans L., Vinkers M., An Ant Algorithm for the Conformational Analysis of
Flexible Molecules., J. Comput. Chem., 2007, 28, pp. 890–898.
71
Dorigo M., Blum C., Ant colony optimization theory: a survey., Theoretical Computer Science, 2005, 344,
pp. 243-278.
72
Gutjahr W.J., A graph-based ant system and its convergence., Future Generation Computer Systems, 2000,
16, pp. 873-888.
73
Gutjahr W.J., ACO algorithms with guaranteed convergence to the optimal solution., Information Processing
Letters, 2002, 82, pp. 145-153.
74
Gutjahr W.J., A generalized convergence result for the Graph-based Ant System., Probability in the
Engineering and Informational Sciences, 2003, 17, pp. 545-569.
75
Stützle T., Dorigo M., A short convergence proof for a class of ACO algorithms., IEEE Trans. Evol. Comput.,
2002, 6, pp. 358-365.
76
Sebastiani G., Torrisi G.L., An extended ant colony algorithm and its convergence analysis., Methodology
and Computing in Applied Probability, 2005, 7, pp. 249-263.
77
Merkle D., Middendorf M., Modelling the Dynamics of Ant Colony Optimization., Evolutionary
Computation, 2002, 10, pp. 235-262.
78
Gutjahr W.J., On the finite-time dynamics of ant colony optimization., Methodology and Computing in
Applied Probability, 2006, 8, pp. 105-133.
79
Gutjahr W.J., First steps to the runtime complexity analysis of ant colony optimization., Computers and
Operations Research, 2008, 35, pp. 2711-2727.
80
Gutjahr W.J., A converging ACO algorithm for stochastic combinatorial optimization., Proc. SAGA 2003
(Stochastic Algorithms: Foundations and Applications), Springer LNCS 2827, 2003, pp. 10-25
81
Gutjahr W.J., S-ACO: An ant-based approach to combinatorial optimization under uncertainty, Proc. ANTS
2004 (4th International Workshop on Ant Colony Optimization and Swarm Intelligence), 2004, Springer
LNCS 3172, pp. 238-249.
82
Rauner M., Brailsford S.C., Gutjahr W.J., Zeppelzauer W., Optimal screening policies for diabetic retinopathy
using a combined discrete-event simulation and ant colony optimization approach., Proc. Int. Conference on
Health Sciences Simulation, Western MultiConference '05, eds.: Anderson J.G., Katzper M., 2005, pp. 147-
152.
83
Brailsford S.C., Gutjahr W.J., Rauner M., Zeppelzauer W., Combined Discrete-event Simulation and Ant
Colony Optimisation Approach for Selecting Optimal Screening Policies for Diabetic Retinopathy.,
Computational Management Science, 2007, 4, pp. 59-83.
84
Korošec P., Šilc J., Real-parameter optimization using stigmergy., In: Proceedings of the Second
International Conference on Bioinspired Optimization Methods and their Application, 2006, pp 73–84.

-154-
85
Russel C. E., Kennedy J., A new optimizer using particle swarm theory.. In: Proceedings of the Sixth
International Symposium on Micro Machine and Human Science, 1995, pp 39–43.
86
Kennedy J., Russel C. E., Particle swarm optimization., In: Proceedings of IEEE International Conference on
Neural Networks, 1995, 4, pp. 1942–1948.
87
Venter G., Sobieszczanski-Sobieski J., Particle swarm optimization., AIAA Journal, 2003, 41, pp. 1583–1589.
88
Cai T., Pan F., Chen J., Adaptive particle swarm optimization algorithm., In: Proceedings of Fifth World
Congress on Intelligent Control and Automation, 2004, 3, pp. 2245–2247,
89
Russel C. E., Yuhui Shi, Comparison between genetic algorithms and particle swarm optimization., In:
Proceedings of the 7th International Conference on Evolutionary Programming, 1998, pp. 611–616.
90
Angeline P. J., Evolutionary optimization versus particle swarm optimization: Philosophy and performance
differences., In: Proceedings of the 7th International Conference on Evolutionary Programming VII, 1998, pp.
601–610.
91
Kennedy J., Spears W. M., Matching Algorithms to Problems: An Experimental Test of the Particle Swarm
and Some Genetic Algorithms on the Multimodal Problem Generator., IEEE International Conference on
Evolutionary Computation, Anchorage, Alaska, 1998.
92
Shi, Y., Eberhart R. C., A Modified Particle Swarm Optimizer., IEEE International Conference on
Evolutionary Computation, Anchorage, Alaska, 1998.
93
Yuhui Shi, Russel C. E., Parameter Selection in Particle Swarm Optimization, Evolutionary Programming
VII, Lecture Notes in Computer Science, 1998, pp. 591-600.
94
Clerc M., The Swarm and the Queen: Towards a Deterministic and Adaptive Particle Swarm Optimization.,
IEEE Congress on Evolutionary Computation, Washington D.C., 1999, pp. 1951-1957.
95
Kennedy J., Mendes R., Population structure and particle swarm performance., IEEE Congress on
Evolutionary Computation (CEC), 2002.
96
Kennedy J, Mendes R., Neighborhood topologies in fully-informed and best-of-neighborhood particle
swarms., IEEE International Workshop on Soft Computing in Industrial Applications, 2003, pp. 45-50.
97
Watts D. J., Small Worlds: The Dynamics of Networks Between Order and Randomness., Princeton
University Press, 1999.
98
Watts D. J., Strogatz S. H., Collective dynamics of small-world networks., Nature, 1998, 393, pp. 440-442.
99
Suganthan P. N., Particle swarm optimiser with neighbourhood operator., IEEE Congress on Evolutionary
Computation (CEC), 1999, pp. 1958-1962.
100
Li X., Adaptively choosing neighbourhood bests using species in a particle swarm optimizer for multimodal
function optimization., Lecture Notes on Computer Science, 2004, 3102, pp. 105-116.
101
Fourie P. C., Groenwold A. A., Particle Swarms in Size and Shape Optimization., Proceedings of the
International Workshop on Multidisciplinary Design Optimization, Pretoria, South Africa, 2000, pp. 97-106.
102
Fourie P. C., Groenwold A. A., Particle Swarms in Topology Optimization., Extended Abstracts of the Fourth
World Congress of Structural and Multidisciplinary Optimization, Dalian, China, 2001, pp. 52-53.
103
Meissner M., Schmuker M., Schneider G., Optimized particle swarm optimization (opso) and its application
to artificial neural network training., BMC Bioinformatics, 2006, 7, 125.
104
Rasmussen T.K., Krink T., Improved hidden markov model training for multiple sequence alignment by a
particle swarm optimization-evolutionary algorithm hybrid., BioSystems, 2003, 72, pp. 5–17.
105
Parsopoulos K. E., Vrahatis M. N., Recent approaches to global optimization problems through particle
swarm optimization., Natural Computing, 2002, 1, pp. 235–306.
106
Wu G., Hansen V., Kreysa E., Gemünd H.-P., Optimierung von fss-bandpassfiltern mit hilfe der
schwarmintelligenz (particle swarm optimization)., Advances in Radio Science, 2006, 4, pp. 65–71.
107
Baltar A. M., Fontane D. G., A generalized multiobjective particle swarm optimization solver for spreadsheet
models: application to water quality., In: Proceedings of AGU Hydrology Days 2006, Fort Collins, Colorado,
2006, pp. 1–12.
108
Cedeño W., Agrafiotis D. K., Using particle swarms for the development of qsar models based on k-nearest
neighbour and kernel regression., Journal of Computer-Aided Molecular Design, 2003, 17, pp. 255–263.
109
Shen Q., Jiang J.-H., Jiao C.-X., Huan S.-Y., li Shen G., Yu R.-Q., Optimized partition of minimum spanning
tree for piecewise modeling by particle swarm algorithm. qsar studies of antagonism of angiotensin ii
antagonists., Journal of Chemical Information and Modeling / J. Chem. Inf. Comput. Sci., 2004, 44, pp. 2027–
2031.
110
Bremermann H. J., Optimization through evolution and recombination., Self-Organizing systems, 1962, pp.
93–10.
111
Woodrow “Woody” Wilson Bledsoe, The evolutionary method in hill climbing: Convergence rates.,
Technical report, Panoramic Research Inc., Palo Alto, California, USA, 1962.
112
Woodrow “Woody” Wilson Bledsoe, An analysis of genetic populations., Technical report, Panoramic
Research Inc., Palo Alto, California, USA, 1962.
113
Bagley J. D., The Behaviour of Adaptive Systems which employ Genetic and Correlation Algorithms., PhD
thesis, The University of Michigan, College of Literature, Science, and the Arts, Computer and
Communication Sciences Department, Ann Arbor, MI, USA, 1967.

-155-
114
Cavicchio D. J. Jr., Adaptive Search using Simulated Evolution., PhD thesis, The University of Michigan,
College of Literature, Science, and the Arts, Computer and Communication Sciences Department, Ann Arbor,
Michigan, USA, 1970.
115
Frantz D. R., Nonlinearities in Genetic Adaptive Search., PhD thesis, The University of Michigan, Ann
Arbor, MI, USA, 1972.
116
Holland J., Adaptation in Natural and Artifical Systems, University of Michigan Press, Ann Arbor, 1975.
117
De Jong, K.A., An Analysis of the Behavior of a Class of Genetic Adaptative Systems., University of
Michigan, Ann Arbor, MI, Ph.D. Thesis, 1975.
118
Baker, J.E., Reducing Bias and Inefficiency in the Selection Algorithm., In: Proceedings of the Second
International Conference on Genetic Algorithms and their Applications, New Jersey, USA, 1985.
119
Mühlenbein H., Schlierkamp-Voosen D., Analysis of Selection, Mutation and Recombination in Genetic
Algorithms., Technical Report 93/94, GMD, 1993.
120
Chipperfield A.J., Fleming P.J., Pohleim H., Fonseca C.M., Genetic Algorithm Toolbox User’s Guide., ACSE
Research Report No. 512, University of Sheffield, 1994.
121
Berthiau G., Siarry P., A Genetic Algorithm for Globally Minimizing Functions of Several Continuous
Variables., In: Second International Conference on Metaheuristics, Sophia-Antipolis (France), 1997.
122
Bäck T., Fogel D.B., Michalewicz Z. (eds.), Evolutionary Computation 1: Basic Algorithms and Operators.,
Institute of Physics Publishing, 2000.
123
Bäck T., Fogel D.B., Michalewicz Z. (eds.), Evolutionary Computation 2: Advanced Algorithms and
Operators., Institute of Physics Publishing, 2000.
124
Holland J. H., Adaptation in Natural and Artificial Systems., The University of Michigan Press, Ann Arbor,
1975.
125
Quagliarella D., Periaux J., Poloni C., Winter G., Genetic Algorithms and Evolution Strategy in Engineering
and Computer Science: Recent Advances and Industrial Applications., John Wiley & Sons Ltd, 1998.
126
Chelouah R., Siarry P., A Continuous Genetic Algorithm Designed for the Global Optimization of
Multimodal Functions., Journal of Heuristics, 2000, 6, pp. 191–213.
127
Sekhon J. S., Mebane W. R., Genetic Optimization Using Derivatives., Draft, 1998.
128
Levine D., Application of a hybrid genetic algorithm to airline crew scheduling., Computers & Operations
Research, 1996, 23, pp. 547–558.
129
Cleveland G. A., Smith S. F., Using genetic algorithms to schedule flow shop releases., In: ICGA,
Proceedings of the third International Conference on Genetic Algorithms, 1989, pp. 160–169.
130
Kwan R. S. K., Kwan A. S. K., Wren A., Driver scheduling using genetic algorithms with embedded
combinatorial traits., In: Proceedings of 7th International Conference on Computer-Aided Scheduling of
Public Transport, Cambridge/Boston, MA, USA, 1997, 471, pp. 81–102.
131
Beaty S. J., Genetic algorithms for instruction sequencing and scheduling., In: Workshop on Computer
Architecture Technology and Formalism for Computer Science Research and Applications, 1992.
132
Xiao Y. L., Williams D. E., Game: Genetic algorithm for minimization of energy, an interactive program for
three-dimensional intermolecular interactions., Computers & Chemistry, 1994, 18, pp. 199–201.
133
Lopes da Cunha A. G., A multi-objective evolutionary algorithm for solving traveling salesman problems:
Application to the design of polymer extruders., In: Proceedings of the 7th International Conference on
Adaptive and Natural Computing Algorithms, 2005, pp. 189–193.
134
Deaven D. M., Ho K.M., Molecular geometry optimization with a genetic algorithm., Physical Review
Letters, 1995, 75, pp. 288–291.
135
Cagnoni A., Dobrzeniecki A., Poli R., Yanch J., Genetic algorithm-based interactive segmentation of 3D
medical images., Image and VisionComputing, 1999, 17, pp. 881–895.
136
Smigrodzki R., Goertzel B., Pennachin C., Coelho L., Prosdocimi F., Parker W. D. Jr., Genetic algorithm for
analysis of mutations in parkinson’s disease., Artificial Intelligence in Medicine, 2005, 35, pp. 227–241.
137
Yan H., Jiang Y., Zheng J., Peng C., Xiao S., Discovering critical diagnostic features for heart diseases with a
hybrid genetic algorithm., In: Valafar F., Valafar H. (eds.), Proceedings of the International Conference on
Mathematics and Engineering Techniques in Medicine and Biological Scienes, Las Vegas, Nevada, USA,
2003, pp. 406–409.
138
Vinterbo S. A., Ohno-Machado L., A genetic algorithm approach to multi-disorder diagnosis., Artificial
Intelligence in Medicine, 2000, 18, pp. 117–132.
139
Min H., Smolinski T. G., Boratyn G. M., A genetic algorithm-based data mining approach to profiling the
adopters and non-adopters of e-purchasing., In: Smari W. W. (eds.), Information Reuse and Integration, Third
International Conference, 2001, pp. 1–6.
140
Kamrani A., Rong W., Gonzalez R., A genetic algorithm methodology for data mining and intelligent
knowledge acquisition., Computers & Industrial Engineering, 2001, 40, pp. 361–377.
141
Gopalan J., Alhajj R., Barker K., Discovering accurate and interesting classification rules using genetic
algorithm., In: Crone S. F., Lessmann S., Stahlbock R. (eds.), Proceedings of the 2006 International
Conference on Data Mining, Las Vegas, Nevada, USA, 2006, pp. 389–395.

-156-
142
Szpiro G. G., A search for hidden relationships: Data mining with genetic algorithms., Computational
Economics, 1997, 10, pp. 267–277.
143
Cordón O., Herrera F., Sánchez L., Evolutionary learning processes for data analysis in electrical engineering
applications., In: Quagliarella D., Périaux J., Poloni C., Winter G. (eds.), Genetic Algorithms and Evolution
Strategy in Engineering and Computer Science, 1998, pp. 205–224.
144
Chai J., Ma S., Robust epipolar geometry estimation using genetic algorithm., Pattern Recognition Letters,
1998, 19, pp. 829–838.
145
Jedidi A., Caminada A., Finke G., 2-objective optimization of cells overlap and geometry with evolutionary
algorithms., In: Applications of Evolutionary Computing, EvoWorkshops, 2004, pp. 130–139.
146
Hsieh I., Chen K.-C., Wang C. A., A genetic algorithm for the minimum tetrahedralization of a convex
polyhedron., In: Proceedings of the 15th Canadian Conference on Computational Geometry, Halifax, Canada,
2003, pp. 115–119.
147
Kundu S., Seto K., Sugino S., Genetic algorithm based design of passive elements for vibration control., In:
Proceedings of 4th MOVIC Conference, 1998, 3, pp. 1183–1188.
148
Kundu S., Seto K., Sugino S., Genetic algorithm application to vibration control of tall flexible structures., In:
Proceedings of The First IEEE International Workshop on Electronic Design, Test and Applications,
Christchurch, New Zealand, 2002, pp. 333–337.
149
Yuret D., Maza M., A genetic algorithm system for predicting the OEX., Technical Analysis of Stocks &
Commodities, 1994, pp. 58–64.
150
El-Fakihy K., Yamaguchiz H., v. Bochmann G., A method and a genetic algorithm for deriving protocols for
distributed applications with minimum communication cost., In: Proceedings of Eleventh IASTED
International Conference on Parallel and Distributed Computing and Systems, Boston, USA, 1999.
151
Yamamoto L. A. R., Tschudin C., Genetic evolution of protocol implementations and configurations., In:
IFIP/IEEE International workshop on Self-Managed Systems and Services, Nice, France, 2005.
152
Hill T, Lundgren A, Fredriksson R, Schiöth H. B., Genetic algorithm for large-scale maximum parsimony
phylogenetic analysis of proteins., Biochimica et Biophysica Acta, 2005, 1725, pp. 19-29.
153
To C. C., Vohradsky J., A parallel genetic algorithm for single class pattern classification and its application
for gene expression profiling in Streptomyces coelicolor., BMC Genomics, 2007, 8, 49.
154
Herrera F., López E., Rodríguez M. A., A linguistic decision model for promotion mix management solved
with genetic algorithms., Fuzzy Sets and Systems, 2002, 131, pp. 47 – 61.
155
Green M.L., Miller R., Evolutionary molecular structure determination using Grid-enabled data mining.,
Cluster Computing and the Grid, 2004, pp. 328 – 335.
156
Wehrens R., Pretsch E., Buydens L. M. C., Quality Criteria of Genetic Algorithms for Structure
Optimization., J. Chem. Inf. Comput. Sci., 1998, 38, pp. 151-157.
157
Gondro C., Kinghorn B.P., A simple genetic algorithm for multiple sequence alignment., Genetics and
Molecular Research, 2007, 6, pp. 964-982.
158
Wang S., Wang Y., Du W., Sun F., Wang X., Zhou C., Liang Y., A multi-approaches-guided genetic
algorithm with application to operon prediction., Artificial Intelligence in Medicine, 2007, 41, pp. 151-159.
159
Willett P., Genetic algorithms in molecular recognition and design., Trends in Biotechnology, 1995, 13, pp.
516-521.
160
van Batenburg F. H., Gultyaev A. P., Pleij C. W., An APL-programmed genetic algorithm for the prediction
of RNA secondary structure., Journal of Theoretical Biology, 1995, 174, pp. 269-280.
161
Colorni A., Dorigo M., Maniezzo V., A genetic algorithm to solve the timetable problem., Tech. rep. 90-060
revised, Politecnico di Milano, Italy, 1992.
162
Iba H., Akiba S., Higuchi T., Sato T., BUGS: A Bug-Based Search Strategy using Genetic Algorithms.,
Journal of Japanese Society for Artificial Intelligence, 1992, 8, pp. 797-809.
163
Ibrahim W., Amer H., An Adaptive Genetic Algorithm for VLSI Test Vector Selection., Applied Simulation
and Modelling, 2006, 522, pp. 86-92.
164
Glover F., Heuristics for integer programming using surrogate constraints., Decision Sciences 1977, 8, pp.
156–166.
165
Glover F., A template for scatter search and path relinking., In: Hao J.-K., Lutton E., Ronald E., Schoenauer
M., Snyers D. (eds.), Artificial Evolution, Lecture Notes in Computer Science, 1998, 1363, pp. 13–54.
166
Laguna M., Martı R.., Scatter Search. Methodology and Implementations in C., Kluwer Academic Publishers,
2003.
167
Martı R., Laguna M., Glover F., Principles of scatter search. Special issue on scatter search and path
relinking., European Journal of Operational Research, 2006, 169, pp. 359–372.
168
Glover F., Genetic algorithms and scatter search: Unsuspected potentials., Statistics and Computing, 1994, 4
pp. 131– 140.
169
De Jong K.A., Spears W.M., A formal analysis of the role of multi-point crossover in genetic algorithms.,
Annals of Mathematics and Artificial Intelligence, 1992, 5, pp. 1–26.
170
Kita H., A comparison study of self-adaptation in evolution strategies and real-coded genetic algorithms.,
Evolutionary Computation Journal, 2001, 9, pp. 223–241.

-157-
171
Glover F., Genetic algorithms, evolutionary algorithms and scatter search: Changing tides and untapped
potentials., INFORMS Computer Science Technical Section Newsletter, 1998, 19, pp. 7–14.
172
Hansen P., Mladenovic N., Variable neighborhood search: Principles and applications., European Journal of
Operational Research, 2001, 130, pp. 449–467.
173
Moscato P.A., On evolution, search, optimization, genetic algorithms and martial arts: Towards memetic
algorithms., Technical Report Caltech Concurrent Computation Program Report 826, Caltech, Caltech,
Pasadena, CA, 1989.
174
Rego C., Alidaee B., Metaheuristic Optimization via memory and evolution: Tabu search and Scatter Search.,
Kluwer Academic Publishers Boston/Dorndrecht/London, 2005.
175
Herrera F., Lozano M., Molina D., Continuous scatter search: An analysis of the integration of some
combination methods and improvement strategies., European Journal of Operational Research, 2006, 169, pp.
450–476.
176
Perez J. A. M., Batista B. M., Laguna M., Scatter Search Based Metaheuristic for Robust Optimization of the
Deploying of “DWDM” Technology on Optical Networks with Survivability., Yugoslav Journal of
Operations Research, 2005, 1, pp. 65-77.
177
Hamiez J.-P., Hao J.-K., Scatter Search for Graph Coloring., Artificial Evolution, 2002, pp. 195 – 213.
178
García del Amo I. J., López F. G., Torres M. G., Batista B. M., Pérez J. A. M., Vega J. M. M., Martn R. R.,
Data Mining with Scatter Search., Computer Aided Systems Theory, 2005, pp. 199-204.
179
Alvarez A. M., González-Velarde J. L., De-Alba K., Scatter Search for Network Design Problem., Annals of
Operations Research, 2005, 138.
180
Sousa S.H.G., Maschio C., Schiozer D.J., Scatter Search Metaheuristic Applied to the History-Matching
Problem., Society of Petroleum Engineers 102975, 2006.
181
Egea J. A., Rodríguez-Fernández M., Banga1 J. R., Martí R., Scatter search for chemical and bio-process
optimization., Journal of Global Optimization, 2007, 37, pp. 481-503.
182
Cordón O., Damas S., Santamaría J., A Scatter search algorithm for the 3D image registration problem.,
Parallel Problem Solving from NatureVIII, 2004, pp. 471-480.
183
Laguna, M., Marti R., Neural Network Prediction in a System for Optimizing Simulations., IIE Transaction
on Operations Engineering, 2000.
184
Corberán A., Fernandez E., Laguna M., Marti R., Heuristic Solutions to the Problem of Routing School Buses
with Multiple Objectives., Journal of the Operational Research Society, 2001, 53, pp. 427-435.
185
Glover F., Laguna M., Marti R., New Ideas and Applications of Scatter Search and Path Relinking New
Optimization Techniques in Engineering, Onwubolu and Babu (eds.), Springer-Verlag 2004.
186
Rego, C., Leao P., A Scatter Search for the Vehicle Routing Problem., Research Report HCES-08-02, Hearin
Center for Enterprise Science, School of Business Administration, University of Mississippi, 2002.
187
David B. C., Julio B. S., Clara C. R., Nature-inspires components of the Scatter Search., First European
Symposium on Nature-inspired Smart Information Systems, Algarve, Portugal, 2005.
188
Van Laarhoven P. J. M., Aarts, E. H. L., Simulated Annealing: Theory and Applications. Dordrecht,
Netherlands: Kluwer, 1987.
189
Metropolis M., Rosenbluth A., Rosenbluth M., A. Teller and E. Teller, Equation of state calculation by fast
computing machines., Journal of Chemical Physics, 1953, 21, pp. 1087–1092.
190
Kirkpatrick S., Gelatt C.D. Jr., Vecchi M.P., Optimisation by simulated annealing., Science, 1983, 220, pp.
671–680.
191
Cerny V., Thermodynamical Approach to the Traveling Salesman Problem: An Efficient Simulation
Algorithm., Journal of Optimization Theory and Applications, 1985, 45, pp. 41–51.
192
Laarhoven P. J., Theoretical and Computational Aspects of Simulated Annealing., Stichting Mathematisch
Centrum, Amsterdam, 1988.
193
Laarhoven P. J., Aarts E. H., Simulated Annealing: Theory and Applications., D. Reidel Publishing
Company, Dordrecht, Holland, 1987.
194
Dueck G., Scheuer T., Threshold Accepting: A General Purpose Optimization Algorithm., Journal of
Computational Physics, 1990, 90, pp. 161–175.
195
Perttunen C.D., Ramasurbramanian N., Deterministic simulated annealing for system optimization Systems.,
Man and Cybernetics; Decision Aiding for Complex Systems, Conference Proceedings., IEEE International
Conference, 1991, pp. 603 - 608.
196
Kaku I., Xiao Y., Xia G., The Deterministic Annealing Algorithms for Vehicle Routing Problems.,
International Journal of Smart Engineering System Design, 2003, 5, pp. 327-339.
197
Ingber L., Very fast simulated re-annealing., Mathl. Comput. Modelling, 1989, 12, pp. 967-973.
198
Ingber L., Adaptive simulated annealing (ASA): Lessons learned., Control and Cybernetics, 1996, 25, pp. 33-
54.
199
Azencott, R. (ed.), Simulated Annealing: Parallelization Techniques., Chichester: Wiley, 1992.
200
Verhoeven M.G.A., Aarts E.H.L., Parallel Local Search Techniques., Journal of Heuristics, 1996, 1, pp. 43–
65.

-158-
201
Onbaşoğlu E., Özdamar L., Parallel Simulated Annealing Algorithms in Global Optimization., Journal of
Global Optimization, 2001, 19, pp. 27-50.
202
Ram D.J., Sreenivas T.H., Subramaniam K.G., Parallel Simulated Annealing Algorithms., Journal of Parallel
and Distributed Computing, 1996, 37, pp. 207-212.
203
Chardaire P., Lutton J.L., Sutter A., Thermostatistical persistency: A Powerful Improving Concept for
Simulated Annealing Algorithms., European Journal of Operational Research, 1995, 86, pp. 565– 579.
204
Hedar A.-R., Fukushima M., Hybrid Simulated Annealing and Direct Search Method for Nonlinear
Unconstrained Global Optimization., Optimization Methods and Software, 2002, 17, pp. 891-912.
205
Wah B.W., Chen Y., Hybrid constrained simulated annealing and genetic algorithms for nonlinear
constrained optimization., Evolutionary Computation, Proceedings of the 2001 Congress, 2001, 2, pp. 925 –
932.
206
Chiu Y.-C., Chang L.-C., Chang F.-J., Using a hybrid genetic algorithm-simulated annealing algorithm for
fuzzy programming of reservoir operation., Hydrological Processes, 2007, 21, pp. 3162 – 3172.
207
Chen D., Lee C.-Y., Park C. H., Hybrid genetic algorithm and simulated annealing (HGASA) in global
function optimization., Proceedings of the 17th IEEE International Conference on Tools with Artificial
Intelligence, 2005, pp. 126 – 133.
208
Ingber L., Simulated annealing: practice versus theory., Mathl. Comput. Modelling, 1993, 18, pp. 29-57.
209
Deboeck G. J. (ed.) Trading On The Edge, Wiley, 1994.
210
Crama Y., Schyns M., Simulated annealing for complex portfolio selection problems., European Journal of
Operational Research, 2003, 150, pp. 546-571.
211
Goffe W.L., Ferrier G.D., Rogers J., Global optimisation of statistical functions with simulated annealing., J.
Econometrics, 1994, 60, pp. 65-100.
212
Ingber L., Wehner M.F., Jabbour G.M., Barnhill T.M., Application of statistical mechanics methodology to
term-structure bond-pricing models., Mathl. Comput. Modelling, 1991, 15, pp. 77-98.
213
Yamada T., Nakano R., Job-Shop Scheduling by Simulated Annealing Combined with Deterministic Local
Search., In: Kluwer Academic Publishers, USA, 1996, pp. 237-248.
214
Ohlmann J. W., Thomas B. W., A Compressed-Annealing Heuristic for the Traveling Salesman Problem with
Time Windows., INFORMS Journal on Computing, 2007, 19, pp. 80-90.
215
Černý V., Thermodynamical approach to the travelling salesman problem: An efficient simulation algorithm.,
Journal of Optimization Theory and Applications, 1985, 45, pp. 41–51.
216
Buckham B. J., Lambert C., Simulated annealing applications., Seminar presentation: MECH620 Quantitative
Analysis, Reasoning and Optimization Methods in CAD/CAM and Concurrent Engineering, Department of
Mechanical Engineering, University of Victoria, 1999.
217
Wales D.J., Doye J.P.K., Global Optimization by Basin-Hopping and the Lowest Energy Structures of
Lennard-Jones Clusters Containing up to 110 Atoms., Chem. Phys. Lett., 1997, 269, pp. 408-412.
218
Coleman T., Shalloway D., Wu Z., A parallel build-up algorithm for global energy minimizations of
molecular clusters using effective energy simulated annealing., Journal of Global Optimization, 1994, 4, pp.
171-85.
219
Wilson S., Cui W., Applications of simulated annealing to peptides., Biopolymers, 1990, 29, pp. 225-235.
220
Biswas R., Hamann D.R., Simulated Annealing of silicon atom clusters in Langevin molecular dynamics.,
Phys. Rev. B, 1986, 34, pp. 895-901.
221
Bollweg W., Kroll H., Maurer H., Numerical prediction of crystal structures by simulated annealing., In:
Bomze I.M., Csendes T., Horst R., Pardalos P.M. (eds.), Developments in Global Optimization, 1997, pp. 253-
288.
222
Pannetier J., Bassas-Alsina J., Rodriguez-Carvajal J., Calgnaert V., Prediction of crystal structures from
crysal chemistry by simulated annealing., Nature, 1990, 346, pp. 343-345.
223
Romeijn H.E., Zabinski Z.B., Graesser D.L. Neogi S., New Reflection Generator for Simulated Annealing in
Mixed-Integer/Continuous Global Optimization., Journal of Optimization Theory and Applications, 1999, 101,
pp. 403-427.
224
Laplume D., Lamblin D., Simulated annealing applied to the optimal plastic design of circular steel plates.,
International Journal for Numerical Methods in Engineering, 2008, 73, pp. 1047-1060.
225
Monem M. J., Namdarian R., Application of simulated annealing (SA) techniques for optimal water
distribution in irrigation canals., Irrigation and Drainage, 2005, 54, pp. 365-373.
226
Martinez-Alfaro H., Flugrad D. R., Collision-free path planning for mobile robots and/or agvs using
simulated annealing., In: Proceedings of the IEEE International Conference on Humans, Information and
Technology, 1994, 1, pp. 270–275.
227
Malhotra A., Oliver J. H., Tu W., Synthesis of spatially and intrinsically constrained curves using simulated
annealing. ASME Transactions., Journal of Mechanical Design, 1996, 118, pp. 53–61.
228
Martínez-Alfaro H., Valenzuela-Rend´on M., Using simulated annealing for paper cutting optimization., In:
Advances in Artificial Intelligence, Lecture Notes in Computer Science, 2004, 2972, pp. 11–20.
229
Sen M. K., Stoffa P. L., Nonlinear one-dimensional seismic waveform inversion using simulated annealing.,
Geophysics, 1991, 56, pp. 1624–1638.

-159-
230
Mladenović N., Hansen P., Variable neighborhood search., Computers Oper. Res., 1997, 24, pp. 1097–1100.
231
Brimberg J., Hansen P., Mladenović N., Taillard E., Improvements and comparison of heuristics for solving
the Multisource Weber problem., Oper. Res., 2000, 48, pp. 444-460.
232
Hansen P., Mladenović N., J-Means: A new local search heuristic for minimum sum-of-squares clustering.,
Pattern Recognition, 2001, 34, pp. 405-413.
233
Belacel N., Mladenović N., Hansen P., Fuzzy J-Means: A new heuristic for Fuzzy clustering., Pattern
Recognition, 2002, 35, pp. 2193-2200.
234
Hansen P., Mladenović N., Perez-Britos D., Variable Neighborhood Decomposition Search., Journal of
Heuristics, 2001, 7, pp. 335-350.
235
Avanthay C., Hertz A., Zufferey N., Variable neighborhood search for graph coloring., European J.O pernl.
Res., 2003, 151, pp. 379–388.
236
Hansen P., Jaumard B., Mladenović N., Parreira A., Variable neighbourhood search for Weighted maximum
satisfiability problem., Les Cahiers du GERAD G- 2000-62, Montr.eal, Canada, 2000.
237
Lopez F.G., Batista B.M., Moreno Perez J.A., Moreno Vega J.M., The parallel Variable Neighborhood
Search for the p-Median problem., Research Report, University of La Laguna, Spain, 2000.
238
Crainic T., Gendreau M., Hansen P., Hoeb N., Mladenović N., Parallel Variable neighborhood search for the
p-Median, MIC’2001, 2001, pp. 595-599.
239
Rodriguez I., Moreno-Vega M., Moreno-Perez J., Heuristics for Routing-median problems., SMG Report,
Universite Libre de Bruxelles, Belgium, 1999.
240
Burke E., De Causmaecker P., Petrovićand S., Berghe G.V., Variable neighborhood search for nurse rostering
problem, MIC’2001, 2001, pp. 755-760.
241
Brimberg J., Hansen P., Mladenović N., Taillard E., Improvements and comparison of heuristics for solving
the Multisource Weber problem., Oper. Res, 2000, 48, pp. 444-460.
242
Paraskevopoulos D. C., Repoussis P. P., Tarantilis C. D., Ioannou G., Prastacos G. P., A reactive variable
neighborhood tabu search for the heterogeneous fleet vehicle routing problem with time windows Journal of
Heuristics, Journal of Heuristics, 2007.
243
Martins S.L., Resende M.G.C., Ribeiro C.C., Pardalos P., A parallel GRASP for the Steiner tree problem in
graphs using a hybrid local search strategy., J. of Global Optimization, 2000, 17, pp. 267-283.
244
Ribeiro C., Uchoa E., Werneck R., A hybrid GRASP with perturbations for the Steiner problem in graphs.,
Technical report, Computer Science Department, Catholic University of Rio de Janeiro, 2001.
245
Festa P., Pardalos P., Resende M., Ribeiro C., GRASP and VNS for Max-cut, Proceedings of MIC’2001,
2001, pp. 371-376.
246
Loudin S., Boizumault P., VNS/LDS + CP: A hybrid method for constraint optimization in anytime contexts,
MIC’2001, 2001, pp. 761-765.
247
Shaw P., Using constraint programming and local search methods to solve vehicle routing problems., In:
Principles and practice of constraint programming, 1998, pp. 417-431.
248
Polacek M., Doerner K. F., Hartl R. F., Maniezzo V., A Variable Neighborhood Search for the Capacitated
Arc Routing Problem with Intermediate Facilities, Journal of Heuristics, 2008.
249
Burke E.K., Cowling P., Keuthen R., Effective local and guided Variable neighborhood search methods for
the asymmetric traveling salesman problem., Lecture Notes in Computer Science, 1999, pp. 203-212.
250
Gonzales C.G., Brito D.P., A Variable neighborhood search for solving the linear ordering problem,
MIC’2001, 2001, pp. 181-185.
251
Avanthay C., Hertz A., Zufferey N., A variable neighborhood search for graph coloring., European Journal
of Operational Research, 2003, 151, pp. 379-388.
252
Santana R., Larranaga P., Lozano J. A., Side chain placement using estimation of distribution algorithms.,
Artificial Intelligence in Medicine, 2007, 39, pp. 49-63.
253
Liberti L., Lavor C., Maculan N., Double VNS for the Molecular Distance Geometry Problem., Proc. of Mini
Euro Conference on Variable Neighbourhood Search, Tenerife, 2005.
254
Belacel N., Cuperlovic-Culf M., Ouellette R., Boulassel M.R., The Variable Neighborhood Search
Metaheuristic for Fuzzy Clustering cDNA Microarray Gene Expression Data., Proceedings of IASTED-AIA-04
Conference. Innsbruck, Austria, 2004.
255
Belacel N., Mladenović N., Hansen P., Fuzzy J-Means: A new heuristic for Fuzzy clustering., Pattern
Recognition, 2002, 35, pp. 2193-2200.
256
Dražić M., Lavor C., Maculan N., Mladenović N., A continuous variable neighborhood search heuristic for
finding the three-dimensional structure of a molecule., European Journal of Operational Research 185, 2008,
pp. 1265–1273.
257
Mladenovic N., Brimberg J., Hansen P., Moreno-Perez J., The p-Median Problem: A Survey of Metaheuristic
Approaches., European Journal of Operational Research, 2007, 179, pp. 927-939.
258
Hansen P., Brimberg J., Urosevic D., Mladenovic N., Primal-Dual Variable Neighbourhood for the Simple
Plant Location Problem., INFORMS Journal on Computing, 2007, 19, pp. 552-564.
259
Puchinger J., Raidl G. R., Pferschy U., The Multidimensional Knapsack Problem: Structure and Algorithms.,
Technical Report No. 006149, National ICT Australia, Melbourne, Australia, 2007.

-160-
260
Fleszar K., Hindi K. S., New heuristics for one-dimensional bin-packing, Computers and Operations
Research, 2002, 29, pp. 821-839.
261
Brimberg J., Hansen P., Lih, K.-W., Mladenović N., Breton M., An Oil Pipeline Design Problem., Operations
Research, 2003, 51, pp. 228-239.
262
Hansen P., Mladenović N., Urosevic D., Variable Neighbourhood Search for the Maximum Clique Problem.,
Discrete Applied Mathematics, 2004, 145, pp. 117-125.
263
Ribeiro C., Souza C., Variable neighborhood descent for the degree-constrained minimum spanning tree
problem., Discrete Applied Mathematics, 2002, 118, pp. 43-54.
264
Festa P., Pardalos P., Resende M., Ribeiro C., GRASP and VNS for Max-cut, Proceedings of MIC’2001,
2001, pp. 371-376.
265
Costa M.-C., Monclar F.-R., Zrikem M., Variable neighborhood search for the optimization of cable layout
problem., MIC’2001, 2001, pp. 749-753.
266
Mladenović N., Urošević D., Variable neighborhood search for the k-cardinality tree., MIC’2001, 2001, pp.
749-753.
267
Brimberg J., Urosevic D., Mladenović N., Variable neighborhood search for the vertex weighted k., European
Journal of Operational Research, 2006, 171, pp. 74-84.
268
Hu B., Leitner M., Raidl G. R., Computing generalized minimum spanning trees with variable neighborhood
search., In: Hansen P., Mladenović N., Moreno Pérez J. A., Batista B. M., Moreno-Vega J. M. (eds.),
Proceedings of the 18th Mini Euro Conference on Variable Neighborhood Search, Tenerife, Spain, 2005.
269
Burke E., De Causmaecker P., Petrovic S., Berghe G.V., Variable neighbourhood search for nurse rostering
problem, MIC’2001, 2001, pp. 755-760.
270
Fleszar K., Hindi K.S., Solving the resource-constrained project scheduling problem by a variable
neighborhood search., European Journal of Operational Research, 2004, 155, pp. 402-413.
271
Leitner M., Hu B., Raidl G.R., Variable neighborhood search for the generalized minimum edge biconnected
network problem., In: Fortz B. (ed.), Proceedings of the International Network Optimization Conference,
2007, Spa, Belgium.
272
Hu B., Leitner L., Raidl G.R., The generalized minimum edge biconnected network problem: Efficient
neighborhood structures for variable neighborhood search., Technical Report TR 186-1-07-02, Institute of
Computer Graphics and Algorithms, Vienna University of Technology, 2007.
273
Audet C., Brimberg J., Hansen P., Le Digabel S., Mladenović N., Pooling Problem: Alternate Formulations
and Solution Methods., Management Science, 2004, 50, pp. 761-776.
274
Glover F., Future paths for integer programming and links to artificial intelligence., Computers and
Operations Research, 1986, 13, pp. 533–549.
275
Hansen P., The Steepest Ascent Mildest Descent Heuristic for Combinatorial Programming., Congress on
Numerical Methods in Combinatorial Optimization, Capri, Italy, 1986.
276
Glover F., Tabu search–Part I., ORSA Journal on Computing, 1989, 1, pp. 190–206.
277
Glover F., Tabu search–Part II., ORSA Journal on Computing, 1990, 2, pp. 4–32.
278
de Werra D., Hertz A., Tabu Search Techniques: A tutorial and an application to neural networks, OR
Spektrum, 1989, pp. 131-141.
279
Glover F., Taillard E., de Werra D., A user's guide to Tabu Search, Annals of Operations Research, 1993, 41,
pp. 3-28.
280
Soriano P., Gendreau M., Diversification Strategies in Tabu Search Algorithms for the Maximum Clique
Problems., Annals of Operations Research, 1996, 63, pp. 189-207.
281
Rochat Y., Taillard E. D., Probabilistic Diversification and Intensification in Local Search for Vehicle
Routing., Journal of Heuristics, 1995, 1, pp. 147–167.
282
Glover F., Laguna M., Marti R., Fundamentals of Scatter Search and Path Relinking., Control and
Cybernetics, 2000, 39, pp. 653–684.
283
(a) Battiti R, Tecchiolli G., The reactive tabu search., ORSA Journal on Computing, 1994, 6, pp.126 -140; (b)
Battiti R, Tecchiolli R., The continuous reactive tabu search: blending combinatorial optimization and
stochastic search for global optimization., Annals of Operations Research, 1996, 63, pp.153-188.
284
Al-Sultan K.S., Al-Fawzan M.A., A tabu search Hooke and Jeeves algorithm for unconstrained optimization.,
European Journal of Operational Research, 1997, 103, pp. 198-208.
285
Hooke R., Jeeves T. A., Direct search solution of numerical and statistical problems, J. Ass. Comput. Mach.,
1961, 8, pp. 212-229.
286
Hedar A., Fukushima M., Tabu Search directed by direct search methods for non-linear global optimization.,
European Journal of Operational Research, 2006, 170, pp. 329-349.
287
Nelder J. A., Mead R., A simplex method for function minimization, Comput. J., 1965, 7, pp. 308-313.
288
Cordeau J.-F., Laporte G., Mercier A., A Unified Tabu Search Heuristic for Vehicle Routing Problems with
Time Windows., Journal of the Operational Research Society, 2001, 52, pp. 928–936.
289
Toth P., Vigo D., The Granular Tabu Search and Its Application to the Vehicle-Routing Problem., INFORMS
Journal on Computing, 2003, 15, pp. 333–346.

-161-
290
Siarry P, Berthiau G., Fitting of tabu search to optimize functions of continuous variables., International
Journal for Numerical Methods in Engineering, 1997; 40, pp. 2449-2457.
291
Cvijovic D., Klinowski J., Taboo search: An approach to the multiple minima problem, Science, 1995, 667,
pp. 664-666.
292
Cvijovic D., Klinowski J., Handbook of Global Optimization, Kluwer: Boston, MA, 2002, pp. 387–406.
293
Franzè F., Speciale N., A tabu-search-based algorithm for continuous multiminima problems., Int. J. Numer.
Meth. Engng, 2001, 50, pp. 665-680.
294
Chelouah R., Siarry P., Tabu search applied to global optimization., European Journal of Operational
Research, 2000, 123, pp. 256-270.
295
Fink A., Vos S., Generic application of tabu search methods to manufacturing problems., In: Proceedings of
the IEEE International Conference on Systems, Piscataway, 1998, 3, pp. 2385–2390.
296
Cao B., Glover F., Tabu search and ejection chains-application to a node weighted version of the cardinality-
constrained TSP., Management Science, 1997, 43, pp. 908–921.
297
Schneider J. J., Kirkpatrick S. (eds.), Tabu Search Applied to TSP, Scientific Computation, 2006, pp. 441–
447.
298
McLoughlin J. F. III, Cedeno W., The enhanced evolutionary tabu search and its application to the quadratic
assignment problem., In: Proceedings of the 2005 conference on Genetic and evolutionary computation, 2005,
pp. 975–982.
299
Glover F., Laguna M., Tabu search., In; Pardalos P. M., Du D.-Z. (eds.), Handbook of Combinatorial
Optimization, Kluwer Academic Publishers, Springer Netherlands, 1998, 3, pp. 621–757.
300
Sexton R. S., Alidaee B., Dorsey R. E., Johnson J. D., Global optimization for artificial neural networks: A
tabu search application., European Journal of Operational Research, 1998, 106, pp. 570–584.
301
Szachniuk M., Popenda L., Gdaniec Z., Adamiak R.W., Błażewicz J., NMR analysis of RNA bulged
structures: Tabu search application in noe signal assignment., In: Proceedings of the 2005 IEEE Symposium on
Computational Intelligence in Bioinformatics and Computational Biology, San Diego, USA, 2005, pp. 172–
178.
302
Hwang G.-J., Yin P.-Y., Yeh S.-H., A tabu search approach to generating test sheets for multiple assessment
criteria., IEEE Transactions on Education, 2006, 49, pp. 88–97.
303
Toth P., Vigo D., The granular tabu search and its application to the vehicle-routing problem., INFORMS
Journal on Computing, 2003, 15, pp. 333–346.
304
Lin B., Miller D. C., Tabu search algorithm for chemical process optimization., Comp. Chem. Eng., 2004, 28,
pp. 2287-2306.
305
a) Baumann K., Albert H., von Korff M., A systematic evaluation of the benefits and hazards of variable
selection in latent variable regression. Part I. Search algorithm, theory and simulations., Journal of
Chemometrics, 2002, 16, pp. 339-350. b) Baumann K., von Korff M., Albert H., A systematic evaluation of
the benefits and hazards of variable selection in latent variable regression. Part II. Practical applications.,
Journal of Chemometrics, 2002, 16, pp. 351-360.
306
Whitley D., Watson J.P., Complexity theory and the no free lunch theorem., In: Burke E.K., Kendall G.
(eds.), Search Methodologies: Introductory Tutorials in Optimization and Decision Support Techniques.
Kluwer, Boston, 2005, pp. 317-339.
307
Jensen F., Introduction to Computational Chemistry, Wiley, New York, 1999.
308
a) Cerjan C. J., Miller W. H., On finding transition states., J. Chem. Phys., 1981, 75, pp. 2800-2806; b)
Simons J., Jorgensen P., Taylor H., Ozment J., Walking on potential energy surfaces., J. Phys. Chem., 1983,
87, pp. 2745-2753; c) Banerjee A., Adams N., Simons J., Shepard R., Search for stationary points on surfaces,
J. Phys. Chem., 1985, 89, pp. 52-57.
309
Kolossvary I., Guida W. C., Low mode search. An efficient, automated computational method for
conformational analysis: Application to cyclic and acyclic alkanes and cyclic peptides., J. Am. Chem. Soc.,
1996, 118, pp. 5011-5019.
310
a) Quapp W., A gradient-only algorithm for tracing a reaction path uphill to the saddle of a potential energy
surface., Chem. Phys. Lett., 1996, 253, pp. 286-292. b) Quapp W., Reduced gradient methods and their
relation to reaction paths, J. Theor. Comput. Chem., 2003, 3, pp. 385-417.
311
Chelouah R., Siarry P., Meta Heuristics Advances and Trends in Local Search Paradigms for Optimization,
Kluwer: Dordrecht, 1999, pp. 49–61.
312
Chelouah R., Siarry P., A hybrid method combining continuous tabu search and Nelder-Mead simplex
algorithms for the global optimization of multiminima functions., Eur. J. Oper. Res., 2005, 161, pp. 636-654.
313
Floudas C. A., Pardalos P. M., Adjiman C. S., Esposito W. R., Gumus Z., Harding S. T., Klepeis J. L., Meyer
C. A., Schweiger C. A. (eds.), Handbook of Test Problems for Local and Global Optimization; Kluwer:
Boston, MA, 1999.
314
Jiang H., Cai W., Shao X., A random tunneling algorithm for the structural optimization problem., Phys.
Chem. Chem. Phys., 2002, 4, pp. 4782-4788.
315
Andressen R. S., Optimization, University of Queensland Press: St. Lucia, Australia, 1972; p. 27.

-162-
316
Korff M., Neue Algorithmen zur Variablenselektion und Variablenwichtung in der Chemometrie.,
Dissertation; Bayerischen Julius-Maximilians-Universität Würzburg, 2001.
317
Stanton A. F, Bleil R. E., Kais S., A New Approach to Global Minimization., J. Comp. Chem., 1997, 18, pp.
594-599.
318
Barhen J., Protopopescu V., Reister D., TRUST: A Deterministic Algorithm for Global Optimization.,
Science, 1997, 276, pp.1094-1097.
319
Barbulescu L., Watson J., Whitley L., Dynamic representations and escaping local optima: Improving genetic
algorithms and local search., The seventeenth National Conference on Artificial Intelligence, Austin, Texas,
2000.
320
Cai W., Shao X., A fast annealing evolutionary algorithm for global optimization., J. Comput. Chem., 2002,
4, pp. 427-435.
321
Potter M., De Jong K., The Third Parallel Problem Solving from Nature III., Springer-Verlag: Berlin, 1994,
pp. 249–257.
322
Mühlebein N., Schomisch M., Born J., The parallel genetic algorithm as function optimizer, J. Parallel
Comput., 1991, 17, pp. 619-632.
323
Voigt H.-M., Soft Genetic Operators in Evolutionary Algorithms., Evolution and Biocomputation, 1995, pp.
123–141.
324
Wolfram Research, Inc. Mathematica Version 5.1; Wolfram Research, Inc.: Champaign, IL, 2004.
325
Zhigljavsky A. A., Theory of Global Random Search., Kluwer: Dordrecht, Netherlands, 1991.
326
Pintér J. D., Global Optimization in Action., Kluwer: Dordrecht, Netherlands, 1996.
327
Törn A., Zilinskas A., Global Optimization., Springer-Verlag: New York, 1989.
328
Ilonen J., Kamarainen J. K., Lampinen J., Differential Evolution Training Algorithm for Feed-Forward
Neural Networks., J. Neural Process Lett., 2003, 17, pp. 93-105.
329
Plagianakos V. P., Vrahatis M. N., Parallel evolutionary training algorithms for. ‘hardware-friendly’ neural
networks., Nat. Comput., 2002, 1, pp. 307-322.
330
Price K., Storn R., Differential Evolution., Dr. Dobb’s J., 1997, 264, pp. 18-24.
331
Storn R., Price K., Differential Evolution - a Simple and Efficient Heuristic for Global Optimization over
Continuous Spaces., J. Global Optimization, 1997, 11, pp. 341-359.
332
Storn R., Sytem Design by Constraint Adaptation and Differential Evolution., IEEE Trans. on Evolutionary
Computation, 1999, 3, pp. 22 - 34.
333
Stepanenko S., Engels B., Gradient Tabu Search, J. Comput. Chem., 2007, 28, 601-611.
334
Dennis J. E., Schnabel R. B., Numerical Methods for Unconstrained Optimization., Philadelphia: SIAM,
1996.
335
Gill P.E., Murray W., Wright M. H., Practical Optimization., San Diego: Academic Press, 1981.
336
Nocedal J., Wright S. J., Numerical Optimization., New York: Springer, 1999.
337
Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., Numerical Recipes in C++, 2nd ed.., New
York: Cambridge University Press, 2002.
338
Stepanenko S., Engels B., New Tabu Search Based Global Optimization Methods Outline of Algorithms and
Study of Efficiency., J. Comput. Chem., 2008, 29, pp. 768–780.
339
Glover F., Tabu Search for Nonlinear and Parametric Optimization (with Links to Genetic Algorithms),
Discrete Applied Mathematics, 1994, 49, pp. 231-255.
340
Laguna M., Martí R., Experimental Testing of Advanced Scatter Search Designs for Global Optimization of
Multimodal Functions., Journal of Global Optimization, 2005, 33, pp. 235-255.
341
Lavor C., Maculan N., A Function to Test Methods Applied to Global Minimization of Potential Energy of
Molecules: Theory and Practice in Optimization., In: Martínez J. M., Yuan J. Y. (guest eds.), Numerical
Algorithms, 2004, 35, pp. 287-300.
342
Chelouah R., Siarry P., Genetic and Nelder–Mead algorithms hybridized for a more accurate global
optimization of continuous multiminima functions., European Journal of Operational Research, 2003, 148,
pp. 335–348.
343
Chelouah R., Siarry P., A continuous genetic algorithm designed for the global optimization of multimodal
functions., Journal of Heuristics, 2000, 6, pp. 191-213.
344
Roger L.S., Srinivas M. S., Rangaiah G.P., Global Optimization of Benchmark and Phase Equilibrium
Problems Using Differential Evoluation., The Journal of the Institution of Engineers, 46.
345
Birbil S.I., Fang S.-C., An Electromagnetism-like Mechanism for Global Optimization., Journal of Global
Optimization, 2003, 25, pp. 263-282.
346
Schlund S., Müller R., Graßmann C., Engels B., Conformational analysis of arginine in gas phase - A strategy
for scanning the potential energy surface effectively., J. Comput. Chem., 2008, 29, pp. 407–415.
347
Schlund S., Schmuck C., Engels B., How Important Is Molecular Rigidity for the Complex Stability of
Artificial Host–Guest Systems? A Theoretical Study on Self-Assembly of Gas-Phase Arginine., Chem. Eur. J.,
2007, 13, pp. 6644-6653.

-163-
348
Helten H., Schirmeister T., Engels B., Theoretical Studies about the Influence of Different Ring Substituents
on the Nucleophilic Ring Opening of Three-Membered Heterocycles and Possible Implications for the
Mechanisms of Cysteine Protease Inhibitors., Journal of Organic Chemistry, 2005, 70, pp. 233-237.
349
Adam W., Bottke N., Engels B., Krebs O., An Experimental and Computational Study on the Reactivity and
Regioselectivity for the Nitrosoarene Ene Reaction: Comparison with Triazolinedione and Singlet Oxygen., J.
Am. Chem. Soc., 2001, 123, pp. 5542-5548.
350
Musch P., Engels B., The Importance of the Ene Reaction for the C2-C6 Cyclization of Enyne-Allenes., J. Am.
Chem. Soc., 2001, 123, pp. 5557-5562.
351
Anet F.A. L., Inflection Points and Chaotic Behavior in Searching the Conformational Space of
Cyclononane., J. Am. Chem. Soc., 1990, 112, pp. 7172-7178.
352
a) Mezey P.G., Potential Energy Hypersurfaces., Elsevier: New York, 1987; b) Boswell D.R, Coxon E.E.,
Coxon J. M., Molecular Modeling of Carbohydrates., Advances in Molecular Modeling, 1995, 3, pp. 195-224.
353
Schlund S., Schmuck C., Engels B., Knock-Out" Analogues as a Tool to Quantify Supramolecular Processes:
A Theoretical Study of Molecular Interactions in Guanidiniocarbonyl Pyrrole Carboxylate Dimers., J. Am.
Chem. Soc., 2005, 127, pp. 11115-111124.
354
Suter H. U., Engels B., Theoretical-Study of Electron-Spin-Resonance Parameters - H2CN and H2CO+.,
Journal of Chemical Physics, 1994, 100, pp. 2936-2942.
355
Schlick T., Optimization methods in computational chemistry., In: Lipkowitz K. B., Boyd D. B. (eds.),
Reviews in Computational Chemistry, VCH Publishers, New York, 1992, 3, pp. 1–71.
356
Ngo J. T., Marks J., Computational complexity of a problem in molecular structure prediction., Protein Eng.,
1992, 5, pp. 313–321.
357
Leach A. R., A survey of methods for searching the conformational space of small and medium-sized
molecules., In: Lipkowitz K. B., Boyd D. B. (eds.), Reviews in Computational Chemistry, VCH Publishers,
New York, 1991, 2, pp. 1–55.
358
Howard A.E., Kollman P.A., An analysis of current methodologies for conformational searching of complex
molecules., J. Med. Chem., 1988, 31, pp. 1669-1675.
359
Leach A. R., Molecular Modelling Principles and Applications., 2nd Edition, Prentice Hall, Pearson
Education, Harlow, England, 2001.
360
Vásquez M., Nemethy G., Scheraga, H. A., Conformational energy calculations on polypeptides and protein.,
Chem. Rev., 1994, 94, pp. 2183–2239.
361
Eisenhaber F., Persson B., Argos P., Protein structure prediction: Recognition of primary, secondary, and
tertiary structural features from amino acid sequence., Crit. Rev. Biochem. Mol. Biol., 1995, 30, pp. 1–94.
362
Bohm G., New approaches in molecular structure prediction., Biophys.Chem., 1996, 59, pp. 1–32.
363
Neumaier A., Molecular modeling of proteins and mathematical prediction of protein structure., SIAM Rev.,
1997, 39, pp. 407–460.
364
Scheraga H. A., Lee J., Pillardy J., Ye Y. J., Liwo A., Ripoll D., Surmounting the multiple-minima problem
in protein folding., J. Global Optim., 1999, 15, pp. 235–260.
365
Floudas C. A., Klepeis J. L., Pardalos P. M., Global optimization approaches in protein folding and peptide
docking., In: DIMACS Series in Discrete Mathematics and Theoretical Computer Science, American
Mathematical Society, New Jersey, 1999, 47, pp. 141–171.
366
Goodman J. M., Still W. C., An Unbounded Systematic Search of Conformational Space, J. Comput. Chem.,
1991, 12, pp. 1110-1117.
367
Bruccoleri R. E., Karplus M., Prediction of the folding of short polypeptide segments by uniform
conformational sampling., Biopolymers, 1987, 26, pp. 137–168.
368
Gibson K. D., Scheraga H. A., Revised algorithms for the build-up procedure for predicting protein
conformations by energy minimization., J. Comput. Chem., 1987, 8, pp. 826–834.
369
Chang G., Guida W.C., Still, W.C., An internal coordinate Monte Carlo method for searching conformational
space, J. Am. Chem. Soc., 1989, 111, pp. 4379-4386.
370
Li Z., Scheraga H.A., Monte Carlo-minimization approach to the multiple-minima problem in protein
folding., Proc. Natl. Acad. Sci. USA, 1987, 4, pp. 6611-6615.
371
von Freyberg B., Braun W., Efficient search for all low energy conformations of polypeptides by Monte
Carlo methods., J. Comput. Chem., 1991, 12, pp. 1065-1076.
372
Senderowitz H., Guarnieri F., Still W.C., A smart Monte Carlo techniques for free energy simulations of
multiconformational molecules. Direct calculations of the conformational populations of organic molecules.,
J. Am. Chem. Soc., 1995, 117, pp. 8211-8219.
373
Ozkan S.B., Meirovitch H., Conformational search of peptides and proteins: Monte Carlo minimization with
an adaptive bias method applied to the heptapeptide deltorphin., J. Comput. Chem., 2004, 25, pp. 565-572.
374
Ripoll D. R., Scheraga H. A., The multiple-minima problem in the conformational analysis of polypeptides.
III. An electrostatically driven Monte Carlo method: Tests on enkephalin., J. Protein Chem., 1989, 8, pp. 263–
287.
375
Wilson S. R., Cui W., Moskowitz J.W., Applications of simulated annealing to the conformational analysis of
flexible molecules., J. Comput. Chem., 1991, 12, pp. 342-349.

-164-
376
Morales L.B., Garduno-Juárez R., Romero D., Applications of simulated annealing to the multiple-minima
problem in small peptides., J. Biomol. Struct. Dyn., 1991, 8, pp. 721-735.
377
Liptom M., Still W.C., The multiple minimum problem in molecular modeling. Tree searching internal
coordinate conformational space., J. Comput. Chem., 1998, 9, pp. 343-355.
378
Yang Y., Liu H., Genetic algorithms for protein conformation sampling and optimization in a discrete
backbone dihedral angles space., J. Comput. Chem., 2006, 27, pp. 1593-1602.
379
Le Grand S.M., Merz K.M., The genetic algorithm and protein structure prediction., In: Merz K. M., Le
Grand S. M. (eds.), the Protein Folding Problem and Tertiary Structure Prediction, Birkhäuser, Boston, 1994,
pp.109-124.
380
Schulze-Kremer S., Genetic algorithm and protein folding., In: Webster D. M. (ed.), Protein Structure
Prediction – Methods and Protocols, Humana Press, New Jersey, 2000, pp. 175-222.
381
Taylor W.R., Aszódi A., Building protein folds using distance geometry: Towards a general modeling and
prediction method., In: Merz K.M., Le Grand S.M. (eds.), The Protein Folding Problem and Tertiary
Structure Prediction, Birkhäuser, Boston, 1994, pp. 165-192.
382
Grippen G.M., Havel, T.F., Distance Geometry and Molecular Conformation, Wiley, New York, 1988.
383
Preishoff C.E., Dixon J.S., Improvements to the distance geometry algorithms for conformational sampling of
cyclic structures., J. Comput. Chem., 1992, 13, pp. 565-569.
384
Weiner P.K., Profeta S., Wipff G., Havel T., Kuntz I.D., Langridge R., Kollman P.A., A distance geometry
study of ring systems: application to cyclooctane, 18-crown-6, cyclododecane and androstanedione.,
Tetrahedron, 1983, 38, pp. 1113-1121.
385
Grippen G. M., Exploring the conformational space of cycloalkanes by linearized embedding., J. Comput.
Chem., 1992, 13, pp. 351-361.
386
a)Piela L., Kostrowicki J., Scheraga H. A., The multiple-minima problem in the conformational analysis of
molecules: Deformation of potential energy hyper surface by diffusion equation method., J. Phys. Chem.,
1989, 93, pp. 3339–3346. b) Kostrowicki J., Scheraga H. A., Application of diffusion equation method for
global optimization to oligopeptides., J. Phys. Chem., 1992, 96, pp. 7442–7449.
387
Sun Y., Kollman P.A., Conformational sampling and ensemble generation by molecular dynamics
simulations: 18-Crown-6 as a test case., J. Comput. Chem., 1992, 13, pp. 33-40.
388
Byrne D., Li J., Platt E., Robson B., Weiner P., Novel Algorithms for Searching Conformational Space., J.
Comput.-aided. Mol. Des., 1994, 8, pp. 67-82.
389
Li Z., Laidig K. E., Daggett V., Conformational Search Using a Molecular Dynamics-Minimization
Procedure: Applications to Clusters of Coulombic Charges, Lennard-Jones Particles, and Waters., Journal of
Computational Chemistry, 1998, 19, pp. 60-70.
390
Gibson K. D., Scheraga H. A., Variable step molecular dynamics: An exploratory technique for peptides with
fixed geometry., J. Comput. Chem., 1990, 11, pp. 468–486.
391
a) Goto H., Osawa E., Corner flapping: a simple and fast algorithm for exhaustive generation of ring
conformations., J. Am. Chem. Soc., 1989, 111, pp. 8950-8951. b) Goto H., Further Developments in the
Algorithm for Generating Cyclic Conformers. Test with Cycloheptadecane., Tetrahedron Lett., 1992, 33, pp.
1343-1346.
392
Goto H., Osawa E., An Efficient Algorithm for Searching Low-energy Conformers of Cyclic and Acyclic
Molecules, J. Chem. Soc., 1993, Perkin Trans. 2, pp. 187-198.
393
Kolossvary I., Guida W.C., Torsional flexing: Conformational searching of cyclic molecules in biased
internal coordinate space., J. Comput. Chem., 1993, 14, pp. 691-698.
394
Kolossvary I., Guida, W.C., Low mode conformational search elucidated: application to C39H80 and flexible
docking of 9-deazaguanine inhibitors into pnp., J. Comp. Chem., 1999, 20, pp. 1671-1684.
395
Kolossvary I., Keseru G., Hessian-free low mode conformational search for large-scale protein loop
optimization: Application to c-jun N-terminal kinase JNK3., J. Comp. Chem., 2001, 22, pp. 21-30.
396
Vengadesan K., Gautham N., A new conformational search technique and its applications., Current Science,
2005, 88, pp. 1759-1770.
397
a) Baker J., Hehre W. J., Geometry optimization in cartesian coordinates: The end of the Z-matrix?, J.
Comput. Chem., 1991, 12, pp. 606-610; b) Baker J., Techniques for geometry optimization: A comparison of
cartesian and natural internal coordinates., J. Comput. Chem., 1993, 14, pp. 1085-1100.
398
Cramer C. J., Essentials of computational chemistry: Theories and models., John Wiley& Sons: Chichester,
2004, 2nd.ed.
399
Pulay P., Fogarasi G. Pang F., Boggs J. E., Systematic ab initio gradient calculation of molecular geometries,
force constants, and dipole moment derivatives., J. Am. Chem. Soc., 1979, 101, pp. 2550-2560.
400
Hehre W. J., Radom L., v. Schleyer P. R., Pople J. A., Ab initio molecular orbital theory., Wiley, New York,
1986.
401
Chass G. A., Sahai M. A., Law J. M. S., Lovas S., Farkas Ö., Perczel A., Rivail J.-L., Csizmadia I. G.,
Toward a computed peptide structure database: The role of a universal atomic numbering system of amino
acids in peptides and internal hierarchy of database., Int. J. Quantum. Chem., 2002, 90, pp. 933-968.

-165-
402
Sahai M. A., Lovas S., Chass G. A., Botond P., Csizmadia I. G., A modular numbering system of selected
oligopeptides for molecular computations: using pre-computed amino acid building blocks., J. Mol. Struct.,
2003, 666-667, pp. 169-218.
403
Frey R. F., Coffin J., Newton S. Q., Ramek M., Cheng V. K. W., Momany F. A., Schäfer L., Importance of
correlation-gradient geometry optimization for molecular conformational analyses., J. Am. Chem. Soc., 1992,
114, pp. 5369-5377.
404
Jalkanen K. J., Suhai S., N-Acetyl-L-Alanine N’-Methylamide: A Density Functional Analysis of the
Vibrational Absorption and Vibrational Circular Dichroism Spectra., Chem. Phys., 1996, 208, pp. 81-116.
405
Jalkanen K. J., Nieminen R. M., Knapp-Mohammady M., Suhai S., A vibrational analysis of various
isotopomers of L-Alanyl-L-Alanine in aqueous solution: vibrational absorption (VA), vibrational circular
dichroism (VCD), Raman and Raman optical activity (ROA) spectra., Int. J. Quantum. Chem., 2003, 92, pp.
239-259.
406
Mazur A. K., Abagyan R. A., New Methodology For Computer-Aided Modelling of Biomolecular Structure
and Dynamics. 1. Non-Cyclic Structures., J. Biomol. Struct. Dyn., 1989, 6, pp. 815-832.
407
Abagyan R. A., Mazur A. K., New Methodology for Computer-Aided Modelling of Biomolecular Structure
and Dynamics. 2. Local Deformations and Cycles., J. Biomol. Struct. Dyn., 1989, 6, pp. 833-845.
408
Abagyan R. A., Totrov M. M., Kuznetsov D. A., Icm: A New Method For Protein Modeling and Design:
Applications To Docking and Structure Prediction From The Distorted Native Conformation., J. Comp.
Chem., 1994, 15, pp. 488-506.
409
Chass A. G., Toward a computed structure database: methodology for effective molecular orbital
computations., Journal of Molecular Structure (Theochem), 2003, 666-67, pp. 61–67.
410
Echenique P., Alonso J. L., Definition of systematic approximately separable and modular internal
coordinates (SASMIC) for macromolecular simulation., J. Comput. Chem., 2006, 27, pp. 1076-1087.
411
Rappé A. K., Casewit C. J., Colwell K. S., Goddard III W. A., Skid, W. M., UFF, a Full Periodic Table Force
Field for Molecular Mechanics and Molecular Dynamics Simulations, J. Am. Chem. Soc., 1992, 114, pp.
10024-10039.
412
Casewit C. J., Colwell K. S., Rappé A. K., Application of a Universal Force Field to Organic Molecules., J.
Am. Chem. Soc., 1992, 114, pp. 10035-10046.

-166-