A lower bound error for free-run simulation of the polynomial NARMAX

ABSTRACT A lower bound error for free-run simulation of the polynomial NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous input) is introduced. The ultimate goal of the polynomial NARMAX is to predict an arbitrary number of steps ahead. Free-run simulation is also used to validate the model. Although free-run simulation of the polynomial NARMAX is essential, little attention has been given to the error propagation to round off in digital computers. Our procedure is based on the comparison of two pseudo-orbits produced from two mathematical equivalent models, but different from the point of view of floating point representation. We apply successfully our technique for three identified models of the systems: sine map, Chua's circuit and Duffing–Ueda oscillator. This technique may be used to reject a simulation, if a required precision is greater than the lower bound error, increasing the numerical reliability in free-run simulation of the polynomial NARMAX.


Introduction
The polynomial NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous input) is widely applied to represent all sorts of systems (Billings, 2013;Chen & Billings, 1989). This mathematical representation may be seen as a well-structured type of recursive function, in which terms and parameters are carefully selected (Korenberg, Billings, Liu, & Mcilroy, 1988). Recursive functions are also used to describe and solve many types of systems and problems (Feigenbaum, 1978) and play an essential role in solving the majority of nonlinear dynamics systems (Hammel, Yorke, & Grebogi, 1987;Lozi, 2013;Nepomuceno, 2014).
Numerous researchers are confident in subscribing a chaotic behaviour to various systems through numerical solutions (Ott, 2002). These solutions make use of popular software and computers easily accessible to most researchers. However, Lozi (2013, p. 63-64) says that 'there are many published works whose reliability numerical results is not carefully evaluated'. For example, Lozi (2013, p. 64) says that 'In the simple case of a dynamic discrete system (of Henon map) there are doubts as to the nature of the 'computational results: long unstable pseudo-orbits or strange attractors? ' Nepomuceno (2014) shows, using double precision in a popular software, that a simple sequence of iterations may generate wrong results in steady state, that CONTACT E.G. Nepomuceno nepomuceno@ufsj.edu.br Supplemental data for this article can be accessed here. http://dx.doi.org/10. 1080/21642583.2016.1163296 is, converged to the wrong answer. In that work, a careful way is presented to calculate error propagation in the computation of recursive function. Investigation of propagation error is not a recent issue (Galias, 2013;Faranda, Mestre, & Turchetti, 2012;Goldberg, 1991;Hammel et al., 1987). In fact, there are many works based on deterministic or stochastic tools that provide some confidence in simulation of recursive functions. Many of those techniques are not practical for nonlinear recursive functions with many terms, such as the polynomial NARMAX, or for non-expert users, who like to measure or at least estimate the error.
To deal with the numerical questions in nonlinear dynamic systems, many researchers have used the shadowing property (Chow & Palmer, 1992, 1991Hammel et al., 1987;Sauer, Grebogi, & Yorke, 1997). This property ensures that a pseudo-orbit (a sequence of points calculated from a recursive function map) is a homeomorphism that remains close to the real orbit (Faranda et al., 2012). What fewer researchers report is the fact that this property is valid only for uniformly hyperbolic systems (Lozi, 2013). We present two works that have been developed to address this situation. First, the work of Hammel et al. (1987) used computational effort to prove a theorem which states that a pseudo-orbit of the logistic map (May, 1976) remains at a distance of 10 −7 from the true orbit up to 10 million iterations using double precision. The computational effort is still a strong tool to support results of analysis of nonlinear dynamic systems. Liao & Wang (2014) used 1200 CPUs and an integration algorithm based on an expansion order of 3500 Taylor series, and 4180 digits of precision to perform simulation of the Lorenz equations. The second example is a variant of the shadowing property, in which the domain is analysed by intervals, resulting in property parameter-shifted shadowing (Lozi, 2013).
Despite advances in numerical computation for nonhyperbolic cases, the two cases presented still have limitations. In the first case (Hammel et al., 1987), an issue if the computer test constitutes a sufficient condition can be raised, and therefore the theorem is proved for all cases, or whether it is a necessary condition, the result may not be valid for all cases. This issue of necessary or sufficient condition on computer simulation has recently come up in Nepomuceno (2014). In the second example, Lozi (2013) states that the parameter-shifted shadowing applied to the map tent (Monteiro, 2002) is valid for almost all values of the parameter a, but not for all values. Hammel et al. (1987) present a technique that may be used to estimate an upper bound error of a simulation. In that work, the authors show that for a specific parameter and initial condition for the logistic map, the computed orbit is within a distance of 10 −8 for 10 7 iterates. This is an important computational result, as the logistic map is a non-hyperbolic system. However, this result is presented for few specific initial conditions and parameters, as well in a specific hardware, the Cray X-MP computer. It is easy to show that the logistic map is a particular case of a polynomial NARMAX. It is indeed interesting to find the upper bound error for computational simulations, but it is not always possible. We can now point out two main facts to justify the relevance of this paper. First, we could not find any paper that deals with the error propagation due to finite precision for the polynomial NARMAX. Second, we present an approach to evaluate a lower bound error for recursive functions, as alternative tool to increase the reliability of computational simulation. The technique is based on the fact that the interval extensions (Moore & Moore, 1979) are mathematically equivalent but may exhibit different computer simulation results. The result of using multiple extensions brings a new concept of 'multiple pseudo-orbits'. The technique presents a tool that may be used to reject a simulation. It means that when a simulation presents a lower bound error greater than the expected precision, we have to stop this simulation and undertake other procedures, such as an increase in the precision of software. We applied our technique in three identified models of the systems: sine map, Chua's circuit and Duffing-Ueda oscillator, focusing on systems which present nonlinear behaviour. It was possible to show a moment of the simulation in which there is no more numerical reliability regarding a required precision.

The polynomial NARMAX
A NARMAX model can be written as Chen & Billings (1989): where y n , u n and e n are, respectively, the output, input and noise terms at the discrete time n ∈ N. The parameters k y , k u and k e are the maximum lag considered for output, input and noise. Terms of e n are frequently included in the parameter estimation process, to avoid bias of an estimator. In this work, F [·] is assumed to be a polynomial with non-linearity degree ∈ Z + . Nonlinear systems are frequently modelled using particular cases of the polynomial NARMAX, such as NAR, NARX and NARMA. In this paper, we name them all as polynomial NARMAX.
Although structure selection and parameter detection are decisive steps in system identification, they are not the main aim of this paper. Billings (2013) presents a extensive list of these techniques. In general, model structure is selected using orthogonal techniques as the error reduction ratio and parameters are estimated using prediction error minimization-based techniques. Freerun simulation is often used as a criterion to validate the polynomial NARMAX (Barbosa, Aguirre, Martinez, & Braga, 2011;Billings, 2013). In many cases, n step-ahead free predictions are performed and compared to a specific validation data set. To the best of our knowledge, little attention has been given to the error propagation in such simulations.

Recursive functions
Let n ∈ N, a metric space M ⊂ R, the relation where f : M → M is a recursive function or a map of a state space into itself and x n denotes the state at the discrete time n. The sequence {x n } obtained by iterating Equation (2) starting from an initial condition x 0 is called the orbit of x 0 (Gilmore & Lefranc, 2012). It is clear that if we take x ∈ R n , F of Equation (1) can be seen as a specific case of f in Equation (2).

Interval extension
Let f be a function of real variable x. Moore & Moore (1979) present the following definition.

Definition 2.1: An interval extension of f is an interval valued function F of an interval variable X, with the property
where by an interval we mean a closed set of real numbers Let us give an example.
We present a few examples of interval extension of f : = rX n − (rX n )X n .
In Example 2.2, Equations (4)-(7) are mathematically equivalent, but they have different sequences of its basic arithmetic operations (Rudolph-Lilith & Muller, 2014;Yabuki & Tsuchiya, 2013). As it is well known that on floating point standard, we do not have commutative and distributive operations, they may exhibit different pseudo-orbits (Goldberg, 1991; Institute of Electrical and Electronics Engineers (IEEE), 2008; Overton, 2001). This difference is also explained when we deal with an interval representation of a number. In such cases, there are interval extensions that are equivalent. This leads us to the following definition. Example 2.2: Let the following extension intervals: In this example, only G(X) and H(X) are equivalent interval extensions.

Orbits and pseudo-orbits
Associated with a map we may define an orbit as follows.

Definition 2.3: An orbit is a sequence of values of a map, represented by
The calculation of this orbit is generally carried out by a finite-precision computer. There is no unique pseudoorbit, as there are different hardware, software and numerical precision standard, such as IEEE 754-2008, which may yield different output for each extension interval.
Definition 2.4: Let i ∈ N represents a pseudo-orbit, which is defined by an initial condition, an interval extension of f, some specific hardware, software and numerical precision standard. A pseudo-orbit is an approximation of an orbit and we represent as where δ i,n ∈ R is the error and δ i,n ≥ 0.
A pseudo-orbit define an interval where the true orbit relies. Hence, we may define an interval associated with each value of a pseudo-orbit From Equations (8) and (9), it is clear that

Lower bound error
From what has been presented on the previous section, we may establish the following lemma.
Proof: Let us assume that there is an intersection equals to empty set for some a,b,n. It means that at least one point n belongs to only a pseudo-orbit a or only a pseudoorbit b, such that I a,n ∩ I b,n = ∅. Hence, x n ∈ I a,n and x n / ∈ I b,n or x n ∈ I b,n and x n / ∈ I a,n . But it implies that either I a,n or I b,n is not an pseudo-orbit according to Equation (10), which is a contradiction. Now we present the following theorem.
Proof: Let I a,n = [x a,n − δ a,n ,x a,n + δ a,n ] and Equation (11) may be rewritten aŝ And that completes the proof. Now, we define the lower bound error in the following theorem.
Theorem 3.3: Let two pseudo-orbits {x a,n } and {x b,n } derived from two interval extensions. Let δ α,n = |x a,n − x b,n |/2 be the lower bound error of a map f (x), then δ a,n ≥ δ α,n or δ b,n ≥ δ α,n .
Proof: Suppose that δ a,n < δ α,n and δ b,n < δ α,n . It implies that which is a contradiction to Theorem 3.2. And that completes the proof.
Theorem 3.2 establishes that at least one of the two pseudo-orbits must have an error greater or equal to the lower bound error. This has a practical meaning. If this lower bound error is greater than the required precision, we cannot go on with the simulation without a further analysis. To guarantee the confidence in the simulation, we must prove that one of the pseudo-orbits has the required accuracy. We believe that the lower bound error is a simple and practical tool to increase the reliability of computational simulation for dynamical systems.
Besides that, if we take G and H as equivalent interval extension there is no reason to state that δ a,n δ b,n neither δ b,n δ a,n . So we may expect that there are simulations which both pseudo-orbits present an error greater than lower bound error, and that was what happened in three studied cases in this paper at a specific time n, as we can see in Section 4.

Stop simulation criterion
One of direct consequences of our approach is to yield information that allows the user to stop the simulation when a required numerical precision is no longer satisfied. We present here an example of such a criterion based on the relative error of a sequence of values. This criterion may be easily adapted for different user needs. Let ε α,n be the relative precision at iteration n defined by where n ∈ N. A user must define the minimum precision, ε, required by simulation. It implies that the simulation must be interrupted at any moment when ε α,n > ε. In this work, we adopt ε = 0.001.

Routines
All routines are performed in Scilab 5.5.2 in a double precision and the results, up to discrete time n = 10, are verified by means symbolic computation performed in wxMaxima 13.04.2 with 1000 digits of accuracy. We used a computer with a processor Intel i5-4440 @ 3.1GHz and a Windows 8 operating system. Scilab and wxMaxima are free and may be downloaded from the internet. In the supplemental material, we present all routines used to produce the results of this paper.

Results
In this section, we present the lower bound error applied for three case studies, which exhibit nonlinear dynamics, such as chaos. The idea is to produce two pseudoorbits from two different interval extensions of models obtained from the literature. We select two interval extensions that are equivalent, according to Definition 2.2. The three chosen models are for the systems sine map (Nepomuceno, Takahashi, Amaral, & Aguirre, 2003), Chua's circuit (Aguirre, 1997) and Duffing-Ueda (Aguirre & Billings, 1994).

Sine map
A unidimensional sine map is defined as where α = 1.2π. A polynomial NARMAX identified for this system is given by (Nepomuceno et al., 2003) y n+1 = 2.6868y n − 0.2462y 3 n . Let us consider two equivalent interval extensions of the model (14): in which the underline was used only to stress the difference between G(X n ) and H(X n ). The difference between the two natural extensions was also underlined in the examples of Chua's circuit and Duffing-Ueda oscillator. As already mentioned in Definition 2.3 and Example 2.4, Equations (15) and (16) are mathematically equivalent, but they represent a different sequence of arithmetical operations, which may produce different outcome. These extensions were simulated using X 0 = 0.1 as an initial condition. By iterating Equations (15) and (16), we produce two pseudo-orbits {ŷ G,n } and {ŷ H,n } of the model (14), respectively. Table 1 shows the first 10 values of a simulation performed in Scilab and the lower bound error, whereas Figure 1(a) shows the results for n up to 60. The relative error is presented in Figure 1(b), in which values are in a logarithm scale. The precision required is no longer satisfied when n = 42. The lower bound error is checked up to the 10th iteration by means the method suggested in Section 3.3. We found δ G,10 ≥ δ α,10 and δ H,10 ≥ δ α,10 , which agrees with Theorem 3.3 and additionally points out for the possibility of cases in which the lower bound error occurs for both pseudo-orbits.

Chua's circuit
Chua's circuit (Chua, Wu, Huang, & Zhong, 1993) is an electrical circuit able to reproduce different regimes, such as periodic and chaotic oscillations. The circuit is composed of a nonlinear Chua's diode, which plays the nonlinear role. Equations (17) and (18) present the white-box model of Chua's circuit  (15) and (16), with results for G(X n ) (−×−) and H(X n ) (−o−) and the same initial condition, that is, X 0 = 0.1. n stands for the number of iterations. (b) Evolution of relative error ε α,n of Equation (15) compared to required precision ε = 0.001. The values are plotted using the log 10 . When n ≥ 44 the ε α,n > ε = log 10 (0.001) = −3, and the simulation no longer satisfies the adopted precision.
These interval extensions were simulated using the initial condition X p = 1 for p = 0 . . . 3. By iterating Equations (20) and (21), we produce two pseudo-orbits {ŷ G,n } and {ŷ H,n } of the model (14), respectively. Table 2 presents 10 first values and the lower bound error of a simulation performed using Scilab. The results up to n = 10 were verified using 1000 digits precision and are presented as a supplemental material. Figure 2(a) shows a specific window for the performed simulation,  (20) and (21), with results for G(X n ) (−×−) and H(X n ) (−o−) and the same initial condition, that is, X p = 1, for p = 0 . . . 3. n stands for the number of iterations. The interval of each iteration is 12 × 10 −6 s. We present in this figure the interval of simulation between 7.2 and 12 ms. (b) Evolution of relative error ε α,n of Equation (15) compared to required precision ε = 0.001. The values are plotted using the log 10 . When n ≥ 464 the ε α,n > ε, and the simulation no longer satisfies this requirement. This means that after around 5.6 ms the simulation does not satisfy the required. Table 2. Numerical computation in Scilab of the first 10 values of G(X n ) and H(X n ) as presented in Equations (20) and (21).

G(X n )
H(X n ) δ α,n Note: The third column is the lower bound error of these two pseudo-orbits according to Theorem 3.3.
using the criterion presented in Section 3.2, when n = 464 the relative error is greater than the minimum precision (see Figure 2(b)). In other words, using double precision and a relative error of ε = 0.001 we may simulate the polynomial NARMAX equation (20) up to around 5.6 ms in order to guarantee the precision adopted. In this case, we also found that δ G,10 ≥ δ α,10 and δ H,10 ≥ δ α,10 .

Duffing-Ueda oscillator
Let us consider a damped, periodically forced nonlinear Duffing-Ueda oscillator (Billings, 2013): A polynomial NARMAX for the Duffing-Ueda oscillator was identified by (Aguirre and Billings, 1994) y n+1 = 2.1579y n − 1.3203y n−1 + 0.16239y n−2 + 0.0003416u n + 0.0019463u n−1 − 0.0048196y 3 n + 0.003523y 2 n y n−1 − 0.0012162y n y n−1 y n−2 + 0.0002248y 3 n−2 , where u = A cos(kT s ) , n ∈ N and T s = π/60. Let us consider two equivalent extensions of the model (23), which were simulated using X p = 2 for p = 0 . . . 2 as initial conditions given by Equations (24) and (25). By iterating Equations (24) and (25), we produce two pseudo-orbits {ŷ G,n } and {ŷ H,n } of the model (14), respectively. Table 3 presents the first 10 values of a simulation of these interval extensions, performed in Scilab and verified using 1000 digits precision in wxMaxima. The codes are provided as a supplemental material. Figure 3(a) presents a specific window from n = 5000 to n = 5500. Figure 3(b) presents the evolution of relative error, which reaches the adopted precision when n = 3345. Thus, after 175.1 s the simulation of Duffing-Ueda by means of Equation (24) presents a relative error superior than the minimum precision adopted of 0.001. Again, δ G,10 ≥ δ α,10 and δ H,10 ≥  (24) and (25), with results for G(X n ) (−×−) and H(X n ) (−o−) and the same initial condition, that is, X p = 2, for p = 0 . . . 2. n stands for the number of iterations. The interval of each iteration is π/60 s. We present in this figure the interval of simulation between approximately 262.8 and 288.0 s. (b) Evolution of relative error ε α,n of Equation (15) compared to required precision ε = 0.001. The values are plotted using the log 10 . When n ≥ 3345, the ε α,n > ε, and the simulation no longer satisfies this requirement. This means that after around 175.1 s the simulation does not satisfy the required precision. δ α,10 are noticed for this example. − 0.0012162X n X n−1 X n−2 + 0.0002248X 3 n−2 ; H(X n ) = 0.0003416U n + 0.0019463U n−1 + 2.1579X n −1.3203X n−1 + 0.16239X n−2 − 0.0048196X 3 n + 0.003523X 2 n X n−1 − 0.0012162X n X n−1 X n−2 + 0.0002248X 3 n−2 ;

Final remarks
We presented a method to calculate a lower bound error for free-run simulation of the polynomial NARMAX.
Although there are other methods to estimate the error propagation, based on the Taylor expansion or shadowing properties, these methods are complex from the point of view of computation and mathematical approach. In general, these methods try to calculate a value of errors instead of a lower bound. In our case, a lower bound is useful as it consists in a way to interrupt simulation when a minimum precision is no longer satisfied. Our method is based on the comparison of two pseudo-orbits of the same models, but derived from two extension intervals. This brings a new idea to be further explored, the very fact that a polynomial NARMAX has multiple pseudo-orbits. It is clear that it also happens with other mathematical representation, whatever it makes use of recursive functions, which increases the relevance of this observation. When we compare two pseudo-orbits that are equivalent from the point of view of interval analysis, we proved a theorem that the sum of absolute error of each pseudo-orbit is at least equal to the distance of these two pseudo-orbits. This theorem works for each iteration n. We also presented a theorem that states that at least one of the two pseudo-orbits presents an error greater than a value, called as lower bound error, which can be easily calculated as a half of the distance of the two pseudo-orbits in each iteration n. This strategy may also be explored for more than two pseudo-orbits, which belongs to the set of future works of our research group.
We applied the methodology in three cases, which are examples of identified systems obtained from literature, by means of the polynomial NARMAX. The sine map, Chua's circuit and Duffing-Ueda oscillator are well known chaotic systems and have been identified using the polynomial NARMAX. These models have been validated using free-run simulation, but little attention has been given to the error propagation of such models due to round off and truncation operations, which are intrinsic at floating point software that uses IEEE 754-2008 standard. Although, one can think that we can simulate these models for an arbitrary n steps ahead, our results show that it is not the case. For the sine map, we see that when n = 42 iterations the relative error reaches the minimum precision adopted, that is, ε = 0.001. Also, for Chua's circuit, after 464 iterations, which means a time of t ≈ 5.6 ms, the relative error is superior than ε. Finally, it is not different from the Duffing-Ueda system, in which n = 3345, or approximately t = 175.1 s, takes the simulation to a relative error greater than ε. It is interesting to notice that for all these three cases, as the number of iterations increases the relative error reaches the same magnitude of the pseudo-orbit, which means that we have lost the significance of all digits. Moreover, the relative error increases exponentially, as can be seen in Figures 1(b), 2(b) and 3(b). This fact also deserves further investigation.
We have found, for the three studied cases, that both pseudo-orbits present an error greater than the lower bound error at n = 10. We have used wxMaxima to check it. This is a stronger result than that established in Theorem 3.3. We believe that it is related to the fact that we have chosen two interval extensions that are equivalent according to Definition 2.3. We aim at investigating if it would be possible to generalize the results of Theorem 3.3. Finally, we also intend to analyse if the concepts present in this paper may be used to develop a strategy to calculate an upper bound for the simulation error of polynomial NARMAX.
It is important to emphasize that in a great number of published works, the double precision is used to perform the simulation, without taking into account the limits of the precision of these recursive simulations. Therefore to guarantee the reliability of the use of these models, it is recommended that a required precision should be fixed, which may imply in a maximum number of iterations for a specific software using double precision or any other finite precision system. Another important aspect, which explains the fact that in the Supplemental Material we present all of software codes, is that without the access of the software and algorithm it may be impossible to reproduce the results.