Mathematical Problems in Engineering

Volume 2018, Article ID 6909151, 8 pages

https://doi.org/10.1155/2018/6909151

## On the Use of Interval Extensions to Estimate the Largest Lyapunov Exponent from Chaotic Data

^{1}Control and Modelling Group (GCOM), Department of Electrical Engineering, Federal University of São João del-Rei, 36307-352 São João del-Rei, MG, Brazil^{2}Department of Electronic Engineering, School of Engineering, Federal University of Minas Gerais, 31270-901 Belo Horizonte, MG, Brazil

Correspondence should be addressed to Márcio J. Lacerda; rb.ude.jsfu@adrecal

Received 28 November 2017; Revised 12 January 2018; Accepted 28 January 2018; Published 4 March 2018

Academic Editor: Gisele Mophou

Copyright © 2018 Erivelton G. Nepomuceno et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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].