Abstract
A method to estimate the (positive) largest Lyapunov exponent (LLE) from data using interval extensions is proposed. The method differs from the ones available in the literature in its simplicity since it is only based on three rather simple steps. Firstly, a polynomial NARMAX is used to identify a model from the data under investigation. Secondly, interval extensions, which can be easily extracted from the identified model, are used to calculate the lower bound error. Finally, a simple linear fit to the logarithm of lower bound error is obtained and then the LLE is retrieved from it as the third step. To illustrate the proposed method, the LLE is calculated for the following well-known benchmarks: sine map, Rössler Equations, and Mackey-Glass Equations from identified models given in the literature and also from two identified NARMAX models: a chaotic jerk circuit and the tent map. In the latter, a Gaussian noise has been added to show the robustness of the proposed method.
1. Introduction
One of the most applied techniques to confirm whether a system is chaotic or not is based on the estimation of the Lyapunov exponents [1]. The calculation of Lyapunov exponents is usually based on the average divergence or convergence of close trajectories along certain directions in state space. In chaotic systems, two trajectories separate exponentially with time despite having very similar initial conditions. Since the work of [2] several numerical methods to estimate Lyapunov exponents were proposed [3–8], just to mention a few. A comparison between estimation methods can be found in [9], for instance. The authors in [5] developed the first practical algorithm to calculate the Lyapunov exponents by estimating the growth of the corresponding set of vectors as the system evolves. This method allows the estimation of the complete spectrum of Lyapunov exponents. Amongst these exponents, the (positive) largest Lyapunov exponent (LLE) is the exponent considered to be the main reason for the separation rate. Therefore the estimation of such an exponent is used to build up the chaotic nature of the data under scrutiny. The authors in [1, 7] independently used statistical properties of the local divergence rates of nearby trajectories of the systems under investigation to estimate the LLE. These properties come in handy when the goal is to get high exactness of numerical estimates of the LLE, as shown in the work [10]. Recent and different applications using the LLE can be found in [11–14].
In [15] a straightforward method to compute the LLE using interval extensions and the original equations of motion was proposed. The exponent is estimated from the slope of the line derived from the lower bound error when considering two interval extensions of the original system [16, 17]. The authors of [15] demonstrated that the method is quick and simple to be used and could be considered as an alternative to other algorithms available in the literature. However, it requires the original equations of motion of the system and it can not be applied directly to time series. In order to overcome this limitation, a system identification methodology based on the well-known polynomial NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous inputs) [18] approach is introduced. Such an approach (system representation) has been shown to adequately represent a large number of chaotic systems [18, 19]. A similar method using the same system representation was proposed in [20] but, in their procedure, the Jacobian matrix must be found and its values estimated. For polynomial NARMAX it was argued that the Jacobian matrix is easily obtained but that may not be the case when other types of representation, such as Neural Networks [20], are considered. Using the interval extensions, it is shown that it is possible to put together the advantages pointed out in [20] without the need to calculate the Jacobian matrix to estimate the LLE. Moreover, as NARMAX is based on a polynomial structure, interval extensions are easily generated by simple arithmetic manipulations, such as commutative and associative. In addition to that, as indicated in [15], the proposed method requires similar or less points to achieve convergence when compared to other methods available in the literature.
The rest of this paper is laid out as follows. Preliminary concepts are introduced in Section 2 and it gives the necessary background material to understand the main idea of the work. Then, the proposed method is presented in Section 3. Section 4 presents the results obtained by the proposed approach when applying it to three different systems that have been identified in the literature: sine map, Rössler Equations, and Mackey-Glass Equations. Moreover, two identified polynomial NARMAX models have been considered: a chaotic jerk circuit [21] and the tent map. In the end, a Gaussian noise has been added to demonstrate the robustness of the proposed technique. Finally, the conclusions are given in Section 5.
2. Preliminary Concepts
In this section a brief review on the calculation of the LLE is given along with the recent result on the lower bound error and the interval extensions.
2.1. The LLE Calculation
To review the basic ideas laid out in [7], consider the following equation:where is the set of all other delay vectors in an -neighbourhood of the vector (data from trajectories of the system under investigation) and is the number of elements in . The initial and final point of each window over the time series is given by and , respectively. The LLE can be estimated by searching for a linear scaling in plot versus . In order to find the linear scaling, it is necessary to tune the parameter , to find the neighbourhood, and to define the number of elements in . This will not be necessary for the method proposed in Section 3; however there is the need of finding the equations of motion of the system under investigation and this will be accomplished by using a system identification methodology as reviewed in the next section.
2.2. The Polynomial NARMAX
A NARMAX model can be written as follows [22]:where , , and are, respectively, the output, input, and noise terms at the discrete time . The parameters , , and are the maximum lag considered for output, input, and noise, respectively. Noise terms are frequently included in the parameter estimation process to avoid bias in the estimates. Here is assumed to be a polynomial with nonlinearity degree . Nonlinear systems are frequently modelled using particular cases of the polynomial NARMAX, such as NAR, NARX, and NARMA [16, 23–25]. In this paper, they all are treated under the name polynomial NARMAX. Further details of this mathematical representation can be found in [18].
The great advantage of using the polynomial NARMAX representation is that the problem of finding the relevant terms and the estimates of their coefficients can be cast into a linear regression problem and therefore the well-known least squares method can be applied. In the eighties, Billings [18] devised an efficient method to take the advantage of the classical methods for solving the least squares problem to include the structure selection using the so-called ERR (Error Reduction Ratio) which attempts to measure the contribution of each term to the overall variance of the output signal. Such a method was called OLS (Orthogonal Least Squares method) and has been successfully used to obtain dynamically valid models from data generated by nonlinear systems behaving chaotically (see, e.g., [26, 27]). Additionally, it is straightforward to generate interval extensions from the polynomial NARMAX, as seen in Definition 1.
2.3. The Lower Bound Error
This subsection is an adaptation from [16, 17] on the lower bound error applied to continuous nonlinear systems.
Let , a metric space ; the relationwhere , is a recursive function or a map of a state space into itself and denotes the state at the discrete time . The sequence obtained by iterating (3) starting from an initial condition is called the orbit of [28, 29].
Definition 1 (see [16]). An interval extension of is an interval valued function of an interval variable , with the propertywhere an interval is a closed set of a real number , where
Example 2. Let us consider . Some examples of interval extensions are
Equations (6) are mathematically equivalent but they present different sequence of their basic arithmetic operations. On floating point standard neither commutative nor distributive properties hold [30, 31]. Therefore, (6) may exhibit different results when a computer-based approach is used to estimate their values. We have exploited these features in this work to estimate the LLE.
Definition 3. Let represent a pseudo-orbit, which is defined by an initial condition, an interval extension of , some specific hardware and software, numerical precision standard, and discretization scheme. A pseudo-orbit is an approximation of an orbit and can be represented by or , such thatwhere is the error and .
A pseudo-orbit characterises an interval where the genuine orbit lies. Henceforth an interval related to each estimation of a pseudo-orbit is given byFrom (7) and (8) it is clear that
Theorem 4 [16] establishes that at least one of two pseudo-orbits must have an error greater than or equal to the lower bound error.
Theorem 4. Let two pseudo-orbits and be derived from two interval extensions. Let be the lower bound error of a map ; then or .
3. The Proposed Method
The method proposed in this work can be summarized by the following steps:(1)Fit a model to the data series. Polynomial NARMAX models were the choice in this work.(2)Choose two different interval extensions of the system under investigation following Definition 1.(3)With exactly the same initial conditions simulate the two interval extensions.(4)Use the least squares method to fit a line to the slope of the logarithm curve of the absolute value of the lower bound error defined in Section 2.3. The slope of the line is the LLE.
To the interested reader, a detailed description and explanation of why this simple method works may be found in [15, 32].
4. Results
In this section, some examples are given to illustrate the usefulness of the proposed methods.
4.1. Example 1: Sine Map
A unidimensional sine map is defined as
The LLE from the original equation (10) was calculated using the method proposed in [15]. To this end, two interval extensions given by the following equations were considered:
The curve of the logarithm of the error bound is depicted in Figure 1(a) and slope of the linear regression, that is, the LLE, is shown to be 1.133. This value is in good agreement with the values found in the literature such as the one given in [33].
(a) Simulation of (11), with results for and and the same initial condition; that is, . LLE calculated is the inclination of the curve given by 1.133.
(b) Simulation of (13), with results for and and the same initial condition; that is, . LLE calculated is the inclination of the curve given by 1.114.
From the data generated by iterating the model in (10) the following polynomial NARMAX given in [34] was identified:Considering two equivalent interval extensions of the identified model (12) we haveFigure 1(b) shows the logarithm of the lower bound error for the two interval extensions given by (13). Table 1 summarizes the results. Note that the similarity of the LLE can be also used to validate the model.
4.2. Example 2: Rössler Equations
Rössler’s Equations [35] are given byThis system has been identified in [33] and the resulting polynomial NARMAX model is given by the following equation:
Equation (15) was considered as the first interval extension. The second interval extension is determined by rewriting the termin (15) asfor instance. Several other interval extensions could be determined but two of them would suffice for the purposes of finding the LLE.
As can be seen the interval extensions used to simulate the systems are mathematically equivalent, but they differ in the sequence of arithmetical operations. Figure 2 shows the logarithm of the lower bound error of the simulation using the two specified interval extensions. Table 2 summarizes the result of LLE for Rössler Equations. The literature suggests a value of 1.242 and the estimated value using the model procedure suggested in [33] is . The estimated value seems slightly different, but if the confidence interval to calculate the LLE using data is taken into consideration, one could consider instead and therefore the value of LLE calculated in [33] is in good agreement with the value in the literature. Similar reasoning can be applied to the value estimated using the method presented in this paper, that is, 1.530.
4.3. Example 3: Mackey-Glass Equations
The Mackey-Glass is an interesting system used in many papers as an example of a chaotic and infinite dimensional system, since it is a time-delay system [36, 37]. The system equation is given bywhere , , , and . The authors in [20] identified the following model to this system using a sequence of 234 data points sample at :
The interval extension used in this example is obtained just by replacing the term by which means that only the order of multiplication has been changed. Figure 3 presents the logarithm of the error for the Mackey-Glass Equation. Using the method proposed here the LLE is found to be exactly the same as the value provided by [20], as reported in Table 3. It is worth emphasizing that there is no need to calculate the Jacobian matrix of the model in the proposed method.
4.4. Example 4: A Chaotic Jerk Circuit
A chaotic jerk circuit given in [21] was also used to further illustrate the proposed approach. The circuit is composed of resistors, capacitors, and operational amplifiers as shown in Figure 4. In this work, a polynomial NAR with nonlinearity degrees and was used to identify the variable of the chaotic jerk circuit, as shown in Figure 5. The data was produced by means of the simulation of the electronic circuit using LTspice. The obtained model has terms, which have been omitted in this manuscript for the sake of simplicity.
For this case, the interval extensions can be found by changing the order in which the model terms were computed. The two interval extensions were generated by the simple order commutation of elements in the polynomial, as illustrated in Example 2. The interval extensions used to simulate the systems are mathematically equivalent, but are not computationally equivalent. Figure 5 shows the logarithm of the lower bound error of the simulation using the two specified interval extensions and Table 4 summarizes the results of the LLE for the chaotic jerk circuit. For comparison purposes [21] shows that, by using experimental data, the value of LLE can be estimated as . On the other hand, the proposed methodology estimates the LLE value as which is in a good agreement with the value for the original system.
4.5. Example 5: The Tent Map
To verify the robustness of the proposed approach under the presence of noise, a tent map (see (22)) was used as a system to be identified, within an additive zero mean Gaussian noise with standard deviation .
A polynomial NAR was identified with and , yielding the model The interval extensions can be found by changing the order in which two terms were computed. In this example, the interval extensions were generated by the simple position change of the terms and Figure 6 presents the logarithm of the error for the identified model whereas Figure 7 shows the first return map of both identified and original tent map.
Despite the presence of noise affecting the system, the LLE found by the method proposed in this paper was in good agreement with the one analytically computed from the tent map, as shown in Table 5.
5. Conclusion
This paper presents a system identification approach to calculate the LLE based on the interval extensions. The method proposed here can be seen as an extension of the method developed in [15] in which the LLE is calculated by taking into consideration the properties of the exponential growing of the lower bound error [16]. The method is designed specifically for the calculation of the LLE from chaotic data, which was not possible to undertake with method proposed in [15]. A polynomial NARMAX model was used to identify the system, but other representations could also be used in the first step. Nevertheless, the polynomial NARMAX represents a very suitable mathematical representation to produce interval extensions, which can be obtained by simple manipulation upon the terms of the identified polynomial model. Compared to existing approaches, such as the one indicated in [20], the proposed method does not require the calculation of the Jacobian matrix and, more importantly, is very simple to use. We believe that the proposed method is also a good alternative for computing Lyapunov exponents from limited data. However, the computation of Lyapunov-exponent spectrum, as performed in [38], is left for future investigations.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.