**2.1. Chaos, Thom's catastrophes and bifurcations**

Bifurcations just described in the previous section can be modeled by Thom's catastrophe theory. This theory, [1]-[4], describes sudden changes in the system behavior under slightly changing (usually) external conditions. These changes, depending on one or more parameters, can be modeled like the special N dimensional surfaces in the so-called parameter space. Each point in the parameter space represents one of the possible system configurations. Those points which are part of the so-called catastrophic fold (surface, plane, ...) are related to system parameter configuration when a system is changing its behavior (moment when bifurcation occurs). When a system control parameter is changed then the point in the parameter space moves and the moment when it crosses through the catastrophic fold changes, then behavior of the system is changed. This change can mean in reality changes in periodicity as well as switching to chaotic dynamic and/or also more drastic changes in the systemâ's physical structure and behavior. Such changes then can lead, in reality, to real catastrophes like aircrafts crashing, dam failure, collapse of a building or power network, etc. As demonstrated in [1] on the Lorenz system (weather model), born of chaos can be understood as a way through the series of bifurcations Thom's catastrophes. Mutual relation between Thom's catastrophes and bifurcations is thus clear.

**Figure 1.** Catastrophe surface - Fold.

## **2.2. Experiment setup**

4 Name of the Book

This part introduces an overview of possible use of SA on bifurcation, i.e. catastrophic events, detection (for a more detailed description, see [43] ). Catastrophic events here means Thom's catastrophes that can be used under certain conditions to model chaotic dynamics and bifurcations, that appear in the nonlinear behavior of various dynamical systems. The main aim of this work is to show that SA are capable of the bifurcations of chaotic system identification without any partial knowledge of internal structure, i.e. based only on measured data. In this part we will discuss mainly SA results. The system selected for numerical experiments here is the well-known system logistic equation derived from predator - prey system. For each algorithm and its version, simulations have been repeated 50 times. Our world, mostly consisting of nonlinear systems, is full of our i.e. human, technology that is less or more reliable. Technological systems are mostly, like their natural counterparts, nonlinear and complex and very often show chaotic as well as catastrophic behavior according to Thom's catastrophe theory, [1], [3] and [4], that describes sudden changes in the dynamical system (well developed and described for so called gradient systems) behavior under slightly changing (usually) external conditions. These changes, depending on one or more parameters, can be modeled like the special N dimensional surfaces in the so-called parameter space, see for example [1]. As an example of such systems (and catastrophic events), we can mention systems such as electrical networks (blackout,...), economic systems (black Friday, NY stock market 1929,...), weather systems (Lorenz model of weather born via series of bifurcations modeled by Thom's catastrophes, see [1]), civil construction failure (bridge collapse, etc), complex systems (self-criticality and spontaneous system reconstruction leading to the better energetic stability) and more. Different mathematical models, of which one possibility is the aforementioned Thom's catastrophe theory, model such events. Our aim in this paper is to show that it is possible to use SA to identify such events on mathematical models of such systems and/or it is possible to use SA to design technological systems in such a way that the possibility to reach regimes exhibiting sudden changes in their behavior (i.e. catastrophe events) is minimized. This study is an extension of our previous research in [9] and [10].

Identification of bifurcations has been done in this research so that Lyapunov exponent is calculated in its absolute value, see Figure 6. Zero values on this cost function landscape (multiple global extremes) then indicate for what parameter A bifurcations occur. In order to locate all these zero values SA is used. They advantage is to overcome multiple local extremes that are present in used cost function. They cannot be wiped out, because of chaotic nature of

Bifurcations just described in the previous section can be modeled by Thom's catastrophe theory. This theory, [1]-[4], describes sudden changes in the system behavior under slightly changing (usually) external conditions. These changes, depending on one or more parameters, can be modeled like the special N dimensional surfaces in the so-called parameter space. Each

**end**

x\* is an approximation of the optimal solution

**2. Case study 1: Evolutionary identification of bifurcations**

**end** *T* := *α*(*T*) ; **until** *T* < *Tf* ;

studied problem.

**2.1. Chaos, Thom's catastrophes and bifurcations**

Experiments have been designed as described in the Table 1. Each experiment was repeated 50 times and the results are reported in the Results section. All simulations have been done on a special grid computer. This grid computer consists of 16 XServers, each 2x2 GHz Intel Xeon, 1 GB RAM, 80 GB HD i.e. 64 CPUs.

### *2.2.1. Used algorithms*

For the experiments described here (SA) [16] had been used. The control parameter settings have been found empirically and are given in Table 1 (SA).

The structure of the exponents can help assess whether chaotic behavior is present in the system or not. Lyapunov exponent is a measure of the extension or contraction of the *ith* semimajor axis of the ellipsoid. For graphic reasons, Lyapunov exponents are arranged by magnitude, i.e. *λ*<sup>1</sup> ≥ *λ*<sup>2</sup> ≥ ... ≥ *λm*, where *m* is the dimension of the phase space; this is referred to as the Lyapunov exponents spectrum. For a chaotic trajectory, one Lyapunov exponent at least must be positive, although, in addition, the existence of any asymptotic periodicity must be ruled out to confirm the chaotic nature - see, e.g., [44]. In other words, the possibility that the trajectory converges to some periodic orbit with t → ∞ must be eliminated. But it is just this requirement that can pose a problem in practice if the dynamic system is investigated during a finite time interval only. Chaotic systems with more than one Lyapunov exponent are referred to as hyperchaotic, see [2]. Lyapunov exponents are closely related to Kolmogorov-Sinai entropy and kind of fractal dimensions also called Kaplan-Yorke

Simulated Annealing in Research and Applications 75



**Figure 4.** Different behavior can be observed when both trajectories will start in different part of state space. Despite its bigger difference in starting possition (*x*<sup>1</sup> = {0.4, 0.4} and *x*<sup>2</sup> = {0.8, 0.4}) trajectories

The main idea of this part is to identify bifurcation (in the general sense related to Thom's catastrophes, see [1]) by means of EAs. For this purpose a logistic equation has been selected (simplified model of predator-prey system), see (6) and [5]. Logistic equation is very simple model of predator - prey system, that exhibit rich set of complex behavior for different values

*xn*+<sup>1</sup> = *Axn* (1 − *xn*) (6)

**Figure 3.** State space trajectory for a dynamic system with 2 singular points *s*<sup>1</sup> and *s*2. On the position *s*<sup>1</sup> = {0, 0} is repellor and at the position *s*<sup>2</sup> = {−1, 0} saddle. Start points of both trajectories diverge

despite fact that this coordinates (*x*<sup>1</sup> = {-1.56, 0.92} and *x*<sup>2</sup> = {-1.57, 0.92}) are very close.



of control parameter A. For more details see [5].

merge together after certain time.

*2.2.3. The cost function*


0

0.5

1

(Lyapunov) dimension.

**Figure 2.** Catastrophe surface - Swallowtail.


#### **Table 1.** SA setting

#### *2.2.2. Lyapunov exponents*

Lyapunov exponents are another member of the family of universal features of deterministic chaos. They are numbers which basically express the divergence (or also convergence) of the state trajectories of a dynamic system. The exponents can be calculated relatively simply, both for discrete-time systems and for continuous-time systems. As will be explained later, Lyapunov exponents are closely related to the structure of the state space, which (in dynamic systems theory) is represented by an array of arrows determining the time development of the system at each point of the space. The development of the system in this space is then represented by a (usually) continuous curve [25]. Illustrative example of different behavior is depicted in Figure 3 and 4. Figure 3 shows the state space of a simple dynamic system along with two different time developments starting from two different initial conditions, which only differ by *x* = 0.01 in the *x*-axis. The behavior in the two cases is entirely different. Figure 4 shows different behavior. Hence, the behavior of a dynamical system is determined by its physical structure, which in the mathematical description is represented by the state space whose quantifiers can be Lyapunov exponents. If one is to follow colored arrows in Fig. 3, it can be noticed that they are separating with increasing time. On the other hand in Fig. 4 they after certain time occupy the same set of points in the state space, in this case called limit cycle. This observation can be described in a mathematical way by the Lyapunov exponents *λi*, see Equation (5), see [13].

$$\lambda\_i = \lim\_{t \to \infty} \frac{1}{t} \ln \frac{l\_i(t)}{l\_i(0)} \tag{5}$$

The structure of the exponents can help assess whether chaotic behavior is present in the system or not. Lyapunov exponent is a measure of the extension or contraction of the *ith* semimajor axis of the ellipsoid. For graphic reasons, Lyapunov exponents are arranged by magnitude, i.e. *λ*<sup>1</sup> ≥ *λ*<sup>2</sup> ≥ ... ≥ *λm*, where *m* is the dimension of the phase space; this is referred to as the Lyapunov exponents spectrum. For a chaotic trajectory, one Lyapunov exponent at least must be positive, although, in addition, the existence of any asymptotic periodicity must be ruled out to confirm the chaotic nature - see, e.g., [44]. In other words, the possibility that the trajectory converges to some periodic orbit with t → ∞ must be eliminated. But it is just this requirement that can pose a problem in practice if the dynamic system is investigated during a finite time interval only. Chaotic systems with more than one Lyapunov exponent are referred to as hyperchaotic, see [2]. Lyapunov exponents are closely related to Kolmogorov-Sinai entropy and kind of fractal dimensions also called Kaplan-Yorke (Lyapunov) dimension.

**Figure 3.** State space trajectory for a dynamic system with 2 singular points *s*<sup>1</sup> and *s*2. On the position *s*<sup>1</sup> = {0, 0} is repellor and at the position *s*<sup>2</sup> = {−1, 0} saddle. Start points of both trajectories diverge despite fact that this coordinates (*x*<sup>1</sup> = {-1.56, 0.92} and *x*<sup>2</sup> = {-1.57, 0.92}) are very close.

**Figure 4.** Different behavior can be observed when both trajectories will start in different part of state space. Despite its bigger difference in starting possition (*x*<sup>1</sup> = {0.4, 0.4} and *x*<sup>2</sup> = {0.8, 0.4}) trajectories merge together after certain time.

#### *2.2.3. The cost function*

6 Name of the Book

No. of particles 100 *σ* 0.5 Tmin 0.000 Tmax 10000 *α* 0.98

Lyapunov exponents are another member of the family of universal features of deterministic chaos. They are numbers which basically express the divergence (or also convergence) of the state trajectories of a dynamic system. The exponents can be calculated relatively simply, both for discrete-time systems and for continuous-time systems. As will be explained later, Lyapunov exponents are closely related to the structure of the state space, which (in dynamic systems theory) is represented by an array of arrows determining the time development of the system at each point of the space. The development of the system in this space is then represented by a (usually) continuous curve [25]. Illustrative example of different behavior is depicted in Figure 3 and 4. Figure 3 shows the state space of a simple dynamic system along with two different time developments starting from two different initial conditions, which only differ by *x* = 0.01 in the *x*-axis. The behavior in the two cases is entirely different. Figure 4 shows different behavior. Hence, the behavior of a dynamical system is determined by its physical structure, which in the mathematical description is represented by the state space whose quantifiers can be Lyapunov exponents. If one is to follow colored arrows in Fig. 3, it can be noticed that they are separating with increasing time. On the other hand in Fig. 4 they after certain time occupy the same set of points in the state space, in this case called limit cycle. This observation can be described in a mathematical way by the Lyapunov exponents

> *λ<sup>i</sup>* = lim *t*→∞

1 *t*

ln *li* (*t*)

*<sup>l</sup>* (0) (5)

**Figure 2.** Catastrophe surface - Swallowtail.

**Table 1.** SA setting

*2.2.2. Lyapunov exponents*

*λi*, see Equation (5), see [13].

The main idea of this part is to identify bifurcation (in the general sense related to Thom's catastrophes, see [1]) by means of EAs. For this purpose a logistic equation has been selected (simplified model of predator-prey system), see (6) and [5]. Logistic equation is very simple model of predator - prey system, that exhibit rich set of complex behavior for different values of control parameter A. For more details see [5].

$$\mathbf{x}\_{n+1} = A\mathbf{x}\_n \left(\mathbf{1} - \mathbf{x}\_n\right) \tag{6}$$

#### 8 Name of the Book 76 Simulated Annealing – Single and Multiple Objective Problems Simulated Annealing in Research and Applications <sup>9</sup>

The bifurcation diagram, [5] (i.e. system dependence on control parameter), of this system is depicted in Figure 7. The Lyapunov exponents of that system, related to parameter A, are depicted in 5 and 6. Simply: when the Lyapunov exponent is negative, the system has deterministic behaviour, when it is positive, the system is chaotic. When the Lyapunov exponent *λ* = 0 then bifurcation can be observed in the system behaviour. Then, to find these moments, we need to locate those parameters A for which the Lyapunov exponent *λ* = 0. It is enough to use absolute value and then we get Figure 5 and 6. Figure 6 represents the cost function landscape, where values with Lyapunov exponent = 0 are related settings of A for which bifurcation occurs. One can see that this surface is very erratic chaotic, and thus suitable candidates to find cost values with 0 are heuristics such as evolutionary algorithms. Cost function, used for depicting Figure 6, is given by (8) which is an absolute value of (7), see [5], variable *d* in (7) is the difference between two nearby trajectories. Function *f*(*x*) in Equation 7 is in this case logistic equation 6. Graphically is Equation 7 and 8 depicted in Figure 5 and 6. EAs were used to search for *fcost* = 0. Difference between Equation 5 and 7 is mainly in fact that Equation 5 calculate *λ* from function describing given system, while Equation 7 is based on recorded numerical data.

$$\begin{aligned} \lambda &= \frac{1}{n} \ln \frac{d\_n}{d\_0} \\ d\_0 &= \left| \mathbf{x}\_i - \mathbf{x}\_j \right| \\ \dots & \dots \\ d\_n &= \left| \mathbf{x}\_{i+n} - \mathbf{x}\_{j+n} \right| \end{aligned} \tag{7}$$

$$f\_{\text{cost}} = |(\lambda)| \tag{8}$$

**Figure 6.** An absolute value of dependence of Lyapunov exponent on parameter A Equation (6), see (8) Maximum 63358 Average 34728 Minimum 9946

Simulated Annealing in Research and Applications 77

Min Average Max <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>9</sup> 9.2 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>−<sup>2</sup>

Value of arameter A 3.543 3.568 3.632 3.742 3.847 Hit 23 12 32 19 14

**Figure 7.** Composed picture show identified bifurcations (red dashed lines) that are in exact correlation

with Lyapunov exponent, compare Figure 5 and 6 and Table 4, also [9], [10].

**Table 4.** The coordinates (value of parameter *A*) of bifurcations (upper row) and number of its

**Table 2.** SA Cost Function Evaluations

successful (approximative) identifications (bottom row).

**Table 3.** Cost values

**Figure 5.** Dependence of Lyapunov exponent on parameter A of Equation (6), see (7)

#### *2.2.4. Experimental results*

Results obtained in all experiments are reported in Table 2 (SA) and Table 3 that show minimal, maximal and average cost values given by experimentation. Table 4 shows localized positions of bifurcations.

**Figure 6.** An absolute value of dependence of Lyapunov exponent on parameter A Equation (6), see (8)

Maximum 63358 Average 34728 Minimum 9946

**Table 2.** SA Cost Function Evaluations


**Table 3.** Cost values

8 Name of the Book

The bifurcation diagram, [5] (i.e. system dependence on control parameter), of this system is depicted in Figure 7. The Lyapunov exponents of that system, related to parameter A, are depicted in 5 and 6. Simply: when the Lyapunov exponent is negative, the system has deterministic behaviour, when it is positive, the system is chaotic. When the Lyapunov exponent *λ* = 0 then bifurcation can be observed in the system behaviour. Then, to find these moments, we need to locate those parameters A for which the Lyapunov exponent *λ* = 0. It is enough to use absolute value and then we get Figure 5 and 6. Figure 6 represents the cost function landscape, where values with Lyapunov exponent = 0 are related settings of A for which bifurcation occurs. One can see that this surface is very erratic chaotic, and thus suitable candidates to find cost values with 0 are heuristics such as evolutionary algorithms. Cost function, used for depicting Figure 6, is given by (8) which is an absolute value of (7), see [5], variable *d* in (7) is the difference between two nearby trajectories. Function *f*(*x*) in Equation 7 is in this case logistic equation 6. Graphically is Equation 7 and 8 depicted in Figure 5 and 6. EAs were used to search for *fcost* = 0. Difference between Equation 5 and 7 is mainly in fact that Equation 5 calculate *λ* from function describing given system, while

*λ* = <sup>l</sup>

*d*<sup>0</sup> = *xi* − *xj* 

... *dn* = 

**Figure 5.** Dependence of Lyapunov exponent on parameter A of Equation (6), see (7)

Results obtained in all experiments are reported in Table 2 (SA) and Table 3 that show minimal, maximal and average cost values given by experimentation. Table 4 shows localized positions

*2.2.4. Experimental results*

of bifurcations.

*<sup>n</sup>* ln *dn d*0

*xi*<sup>+</sup>*<sup>n</sup>* − *xj*<sup>+</sup>*<sup>n</sup>*

 

*fcost* = |(*λ*)| (8)

(7)

Equation 7 is based on recorded numerical data.


**Table 4.** The coordinates (value of parameter *A*) of bifurcations (upper row) and number of its successful (approximative) identifications (bottom row).

**Figure 7.** Composed picture show identified bifurcations (red dashed lines) that are in exact correlation with Lyapunov exponent, compare Figure 5 and 6 and Table 4, also [9], [10].

10 Name of the Book 78 Simulated Annealing – Single and Multiple Objective Problems Simulated Annealing in Research and Applications <sup>11</sup>

**Figure 8.** Identified and non-identified bifurcations, see also [9], [10].

**Figure 10.** Non-identified bifurcations, zoom from Figure 9, see also [9], [10].

real systems is possible.

set for *λ* < 10−1.

**3. Case study 2: Chaos synthesis**

domain.

3. **Chaotic systems.** Besides well known logistic equation, other systems such as Ikeda, Burger, Delayed Logistic, etc. can be selected. The logistic equation system (derived from predator prey system) has been chosen because it is well known, and well investigated. Other systems (like those mentioned above) are now in the process of investigation. 4. **Simulation environment and software.** Software used for all calculations, numerical simulations and visualizations was the well known *Mathematica*. Thanks to fact that Mathematica is an integrated environment, the time scale for each simulation was in minutes. It is clear that when C++ or another fast programming language would be used, then time scale of simulation should be much shorter and thus real use of this approach on

Simulated Annealing in Research and Applications 79

5. **Localized bifurcations.** It is important to note that only main bifurcations have been localized. Bifurcations at exact positions were undiscovered, which implies space for future research, i.e. better algorithm setting, cost function improvement, etc. Figures 7-10 clearly show, that there are a lot of other bifurcations that were not localized by the SA, so the question is: what is the limit for EAs to be used on this kind of task, or what will be better settings for it? The bifurcations were determined by the same SA in repeated simulations. Localized bifurcation implies that *λ* = 0, however because SA is numerical algorithm, working with certain preciseness, successful localization has been

In the light of our previous experiments described in [11], [12] and [13] it is clear that evolutionary techniques are capable of solving a wider class of problems in the chaos research

This part introduces the notion of chaos synthesis by means of evolutionary algorithms and develops a new method for chaotic systems synthesis, [13]. This method is similar to genetic programming and grammatical evolution and is being applied along with evolutionary algorithms: differential evolution, self-organizing migrating, genetic algorithm, simulated annealing and evolutionary strategies. The aim of this investigation is to synthesize new

**Figure 9.** Identified and non-identified bifurcations, zoom from Figure 8, see also [9], [10].

SA in basic versions has been used in a simple system (6), called logistic equation, to localize bifurcations. Simulations were repeated 50 times. Based on results reported in the previous section it can be stated that:


**Figure 10.** Non-identified bifurcations, zoom from Figure 9, see also [9], [10].

10 Name of the Book

**Figure 8.** Identified and non-identified bifurcations, see also [9], [10].

section it can be stated that:

experiments.

bifurcation event, see Table 4.

**Figure 9.** Identified and non-identified bifurcations, zoom from Figure 8, see also [9], [10].

SA in basic versions has been used in a simple system (6), called logistic equation, to localize bifurcations. Simulations were repeated 50 times. Based on results reported in the previous

1. **Simulations.** All simulations has been successfully finished by localization of the

2. **Precision.** Located bifurcations were localized with high precision, however, not exactly, see Figure 7-10. Geometrically it can be interpreted so that the system is at the edge of a catastrophic model or almost before crossing the fold on the surface of Thom's catastrophe manifold. Our opinion is that such non-precise localization can be useful, when applied to a real system and bifurcation can cause real damage, which means that the evolutionary search for bifurcations on a real system can be stopped very near to the bifurcation position. Figures 7-10 clearly show, that there are a lot of other bifurcations that were not localized by the SA that was used, which imply better settings of SA or acceptable value of *λ* (see Localized bifurcations below), is probably possible. This is open question for further


In the light of our previous experiments described in [11], [12] and [13] it is clear that evolutionary techniques are capable of solving a wider class of problems in the chaos research domain.
