Recovering Chaotic Systems Using Genetic Algorithms

Recovering Chaotic Systems using Genetic Algorithms Emanuele Galiffi & Paul M. Secular, Imperial College London Report submitted: 24th June, 2013

NOTE: In this version figures and tables are not available because of the format. A link to the full article will be added soon.

Abstract

Recovery of the parameters and initial conditions of chaotic solutions to the Lorenz equations was attempted using genetic algorithms. For discrete sets of possible parameter values (including the correct values), a genetic algorithm was shown to converge faster on average than a brute force approach. In concordance with the work of Tahersima et al. [1], it was found to be more difficult to apply this method to a continuous range of parameter values using a real-encoded genetic algorithm. This is probably due to the extreme sensitivity of the Lorenz equations to both initial conditions and parameters.

Introduction

Chaotic systems are commonplace in Nature and are now fairly well understood. However, their complex and unpredictable behaviour continues to present a number of practical difficulties for scientists and engineers. The pioneering work of Barry Saltzman and Edward Lorenz on atmospheric convection during the early 1960s led to the discovery of a system of three nonlinear differential equations with aperiodic solutions [2]. Thus was born the field of study known as Chaos theory.

The equations now bear Lorenz's name and are cited as being one of the first known examples of a chaotic dynamical system [3]. As such, this system is widely used as both an example and as a benchmark in work on chaotic systems. Although the Lorenz equations were assumed to be capable of chaotic behaviour, it was only as recently as 1999 that Warwick Tucker finally proved mathematically that they do indeed admit chaotic solutions [4].

One of the defining characteristics of chaotic dynamical systems is their extreme sensitivity to initial conditions [3]. It is this property which makes long-term data prediction so difficult (the so-called 'butterfly effect'). Estimation of a chaotic system's parameters and/or initial conditions helps increase the range of predictability. This is essential for both short-term weather and long-term climate prediction [7]. Various sophisticated techniques have been proposed to tackle this problem [5] [6] [7].

This paper describes an approach using genetic algorithms. Similar work has previously been undertaken by Tahersima et al. on the Lorenz system [1] and by Caponetto et al. on the Chua system. Caponetto et al. use a genetic algorithm alongside a powerful technique called chaotic system synchronisation to recover approximate values for the Chua system's parameters [8] [9].

In the present work, an attempt is made to recover initial conditions as well as system parameters. Preliminary results appear to agree with those of Tahersima et al. in that precise values prove difficult to recover. Because of this, an approach proposed by Leeuwenburgh [10] was taken, in which only a discrete set of possible parameter values were considered. When this set contained the correct values, the algorithm was usually able to determine them correctly with an average performance far greater than that of a brute force approach.

The Lorenz System

The Lorenz system was conceived as a simplified mathematical model of Rayleigh-Bénard con-vection, i.e. thermal convection in a horizontal layer of fluid under the influence of gravity [2] [11]. It comprises the following set of three coupled, autonomous, nonlinear, ordinary differential equations [2]: {?(x ? = ? (y - x) @ y ? = rx - y - xz @z ? = xy - bz)?

Although the equations are dimensionless in form, the variables x, y and z and the parameters ?, r and b are proportional to physical properties of the fluid flow. x is related to the overall intensity of the convection, y is related to the difference in temperature between rising and falling convection currents, and z relates to the nonlinearity of the temperature profile [2]. ? is the Prandtl number, r is the Rayleigh number, and b is related to the fluid layer's geometry [12] [13].

It turns out that the Lorenz system can also be used to model a number of other, very different, physical phenomena. For example, the Lorenz equations have been shown to be equivalent in form to the Maxwell-Bloch equations which are applicable to the behaviour of lasers under certain special conditions. For a more detailed discussion see [14].

The motivation for applying a genetic algorithm to the Lorenz system comes from the difficulty of recovering initial conditions from chaotic time series. In order to verify this difficulty, attempts were made to retrace solutions backwards using time-reversed 4th order Runge-Kutta and Euler methods (i.e. with a negative time step) as well as an implicit backward Euler method (derived using Cardano's formula for cubic equations). Even when the time step was decreased to many orders of magnitude below that used with the forwards-in-time methods, the backwards trajectories always diverged rapidly towards infinity. In other words, with the direction of time reversed, the Lorenz system appears to act as a repeller instead of an attractor, as might intuitively be expected. This is consistent with the behaviour of the geometric Lorenz attractor model of Dumortier et al. under time reversal [15] and justifies Gouyet's assertion that: "There is then, in some sense, a loss of memory, since after a certain time has elapsed the initial conditions cannot be recovered." [16] Figure 1. Lorenz Attractor: plot of a numerical solution to the Lorenz Equations using the explicit forward Euler method over a timespan of ?t = 12.5 with a time step of 0.001. Parameter values used were: ? = 10, b = 8/3, r = 28 Figure 2. Lorenz Repeller: attempting to plot the same trajectory backwards in time using the implicit backwards Euler method results in divergence from the attractor. Plot shown for ?t = -1.0 using a time step of 0.000001 Genetic Algorithms The idea of evolutionary computation was first studied in the 50s and the 60s as an alternative optimisation method. John Holland finally invented genetic algorithms in the 1960s. His original aim was to use them for a formal study of natural selection, although their evident applicability to computational problems induced several scientists, including Holland himself to develop a theoretical framework [18] for applying them to computer systems. [17] The practical need for such methods arises because the optimisation of a cost function which depends on many parameters is often ineffective using classical numerical or brute force methods. In fact, most of the former cannot escape local minima, while the latter are normally too inefficient for real-life problems. This results in a need for a pseudo-random method which can map a broad range of solutions and choose from among these the ones which best minimise the cost function, as well as escape local minima. Such methods clearly rely on the possibility of evaluating the fitness of a solution; therefore they cannot substitute brute force methods when this is not possible, namely, when the solution can only be right or wrong. [17] Since this method is a very close analogy to the evolutionist theory of Natural Selection, the technical language used when dealing with genetic algorithms includes many words which are typical of biological fields. The most important of these terms are explained in Appendix A. [20] The way genetic algorithms (GA's) work is illustrated in the flowchart in Figure 4. A set of arrays of parameters (solutions) is generated randomly and the cost function is evaluated for each one of them. Some of the very best solutions are saved, to be added to the next generation as they are. This will prevent the fitness of the solution from worsening because of random processes. Once the cost associated with each solution is known, the fittest ones (i.e. with a lower cost function) are selected and allowed to breed together. Randomly, a number of mutations are applied to some solutions, so that new parameters are added at each generation. Each mutation consists of substituting a randomly chosen parameter into a randomly chosen solution. Finally, the best solutions saved from the last generation are added as they are by substituting them for a randomly chosen solution, and a new generation of solutions is ready to repeat the cycle. [20] This is a very general descri ption of how any GA works. In fact, it is a very broad field, and there are several different methods of implementing each of the operations described above. This fact results in GA's being a very case-specific method, and different techniques (e.g. selection methods and mutation rates) are suitable for different types and "stiffness" of problems. The principle which underpins the working of GA's is the survival of the fittest, which allows them to obtain better and better solutions with time. There are two main types of GA implementation: Discrete parameter: the solutions are chosen from a discrete set of possible genes.

Continuous parameter or real-valued: the solutions are chosen from a continuum of possible solutions. In this case the solution can find the minimum valley, but it will only approach the exact value as time approaches infinity, which causes serious problems when dealing with chaotic systems, which are extremely sensitive to small differences in the parameters. This method normally relies on mutations for introducing new values from the continuous range. [20] This project focuses on the former method: an array of possible parameters is formed (see Figure 3), and the aim of the algorithm is to find the best combination. The latter method was also attempted; findings are presented in section VII. x0 y0 z0 ? b r Figure 3. Array representation of a chromosome with six genes. Each gene holds a trial value for one of the Lorenz system's parameters. Figure 4. Flowchart describing the top-level logic of the genetic algorithm. The sub-processes shaded in yellow and green are shown in more detail in figures 5 and 6, respectively. Figure 5. Flowchart showing the process of fitness evaluation for chromosomes. Figure 6. The process of creating a new generation of chromosomes using genetic operators. Application of GA's to chaotic systems The sensitivity of chaotic systems to initial conditions makes GA's a very effective technique for finding the initial conditions and the parameters given a set of data points and a set of discrete possible parameters. The cost function will thus be the sum of the residuals [10] of the curve computed by using the parameters (genes) present in a given solution (chromosome) compared to the correct one. The exponential divergence of different points of a chaotic system [3] might make it possible to reduce the calculation of the curve to a few data points; just enough to overcome the limit required and leave a wide enough range to be able to sort the chromo¬somes by their fitness. Given the wide range of possible implemen-tations of genetic algorithms it is normally very complicated to find a mathematical way of predicting the most suitable features for a given problem such as selection and crossover methods or the most suitable number of chromosomes. The analytical treatment of this efficiency calculation is beyond the scope of this study, which only relies on statistical data to decide the best implementation. Other methods involve meta-optimisation, i.e. the optimisation of an optimisation method (for an investigation of this type, see [19]). However, it is usually possible to rely on common sense and practical evidence in order to deduce which implementations are the most suitable for a given problem. To quote from Randy L. and Sue Ellen Haupt [20, page 64]: "Selecting the various genetic algorithm parameters, such as mutation rate and type of crossover, is more of an art than a science". The next section describes the features of the GA used in this investigation. It is only one possible implementation, but it was compared to different methods (see Appendix B) and found to work best. This does not imply that there are not any better implementations for this problem.

Features of the GA used in this investigation The initial population was 100, whilst the number of chromosomes per generation was 20. This allows the GA to scan a wider range of solutions, thus avoiding most of the local minima. This turns out to be particularly useful in light of the fact that the mating method (see below) allows the worse chromosomes some probability of being chosen. A total of 16 random mutations was applied to every new generation. Elitism was performed by saving the very best chromosome from the previous generation. The effectiveness of mutations is strictly correlated to elitism. In fact, since increasing the number of mutations would introduce more randomness into the GA, this results in a larger probability of accidentally decreasing the fitness of the chromosomes. However, when elitism is applied, the best chromosomes from the previous generation are always present in the new one and do not undergo mutations; therefore this approach would increase the "variety" in the population, without risk of worsening the solution. This turns out to be very helpful when the GA needs to jump out of a local minimum. The particular parent selection method used is called the "Binary Tournament Method". It consists of choosing two chromosomes randomly, among which the one with the best fitness is selected to become a parent. This is performed twice per new chromosome, so that from each pair of parents one new chromosome is produced. This method ensures a wide variety of chromosomes in the offspring, thus allowing the GA to scan a broader set of solutions. The crossover method used is called single point crossover. The way it works is by randomly choosing a position in the chromosomes and taking the genes which are before that position in the array from one parent and the following genes from the other. The part of the chromosome taken from each parent is chosen randomly. Other possible methods, which were not used in the final program but are worth mentioning, are included in Appendix B. Results and discussion The program initially computed the solution to the Lorenz equations using a numerical method with a given set of parameters and initial conditions. The numerical integrator ODEint from the SciPy library was used. A number of possible solutions (i.e. parameters and initial conditions) as well as the original ones are randomly generated.

A genetic algorithm was implemented with the properties above to recover the correct initial conditions and parameters of the Lorenz equations from the set of possible (trial) ones. Statistical data were collected with different numbers of trial values per each parameter, namely two, three, and four. The program was run 100, 200, and 450 times for two, three, and four trial parameters respectively.

The number of possible combinations of a number n of possible values for each of k parameters to recover is n^k. In this case the numbers of possible combinations for two, three and four possible values were 64, 729 and 4096 respectively.

Using the brute force method, every possible combination is evaluated, so that the probability of finding the correct one is independent of time, as the graphs show. Thus, the average time taken to find the best solution with this method will be half of the time needed to check all the possible combinations.

The results for two, three and four trial parameters are shown in Figures 7, 8 and 9 respectively (Unfortunately figures are not available on this website). Evidently, the performance of GA's compared to brute force increases dramatically as the number of possible parameters increases. However, due to the significant level of randomness in the GA, it is not always possible to obtain the best answer at all. In fact, as shown in the graphs, the total probability of finding the best solution only tends to one as the time of run approaches infinity, while brute force methods ensure that the correct answer is found in a finite time, but high numbers of possible combinations normally make these methods of no practical use.

Nevertheless, the results show that (for GA's) there is a much higher probability of finding the answer at the very beginning, which decreases very soon with time. This fact results in the convenience of interrupting GA's and running them again if they do not find the required answer (when this is known within some level of accuracy) within a given amount of time, as was done in this study.

The graphs show how the percentage probability of failure within the time allowed increased with the number of trial parameters. This arises because the distribution tends to spread out more as the number of possible combinations increases.

Although this investigation mostly dealt with a discrete set of possible values for each parameter, a more physically realistic problem would involve finding the parameters from a continuous range. This is discussed in the following section. Figure 7. Histogram showing the % probability per unit time of the program ending in the two trial parameters case. Brute force is still more convenient because of the GA's slower start. Figure 8. Histogram showing the % probability per unit time of the program ending in the three trial parameters case.

Figure 9. Histogram showing the % probability per unit time of the program ending in the four trial parameters case. The last two bins mostly contain failures in the time given. Further remarks

The Lorenz attractor has a fractal structure, which means that in theory there may exist infinitely many solutions that lie arbitrarily close together over any finite time interval. Because of this, it may actually be impossible to recover initial conditions given some finite data set, since there could exist many vastly different initial conditions that may give almost exactly the same time series data. The results shown in table 1 illustrate this fact.

A trial set of initial conditions can be assigned a high "fitness" even if the values differ greatly from the correct ones. In other words, there is no simple monotonic relationship between the average divergence of trajectories and the distance between their initial conditions. x0 y0 z0 Average divergence 4.99990000 4.99990000 4.99990000 0.40040672 9.72245016 3.63308943 4.60335410 0.20291101 0.12401190 2.20139150 5.77109152 0.02205292

Table 1. Comparison of three different "trial" trajectories to one created from the initial conditions: x0 = 5.00000000, y0 = 5.00000000, z0 = 5.00000000.

In each case, identical pa