The Lower Bound Error for polynomial NARMAX using an Arbitrary Number of Natural Interval Extensions
Abstract
Abstract. The polynomial NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous input) is a model that represents the dynamics of physical systems. This polynomial contains information from the past of the inputs and outputs of the process, that is, it is a recursive model. In digital computers this generates the propagation of the rounding error. Our procedure is based on the estimation of the maximum value of the lower bound error considering an arbitrary number of pseudoorbits produced from different natural interval extensions, and a posterior Lyapunov exponent calculation. We applied successfully our technique for two identified models of the systems: sine map and DuffingUeda oscillator.
Keywords. Polynomial NARMAX, Lower Bound Error, Natural Interval Extension, Interval Arithmetic.
The Lower Bound Error for polynomial NARMAX using an Arbitrary Number of Natural Interval Extensions
Priscila F. S. Guedes^{†}^{†}thanks: pri12_guedes@hotmail.com Márcia L. C. Peixoto^{†}^{†}thanks: marciapeixoto93@hotmail.com Alípio M. Barbosa^{†}^{†}thanks: alipiomonteiro@yahoo.com.br Samir A. M. Martins^{†}^{†}thanks: martins@ufsj.edu.br Erivelton G. Nepomuceno^{†}^{†}thanks: nepomuceno@ufsj.edu.br ^{1,2,4,5}Control and Modelling Group (GCOM), Department of Eletrical Engineering, Federal University of São João delRei, MG, Brazil ^{3}Centro Universitário Newton Paiva, Belo Horizonte, MG, Brazil
1 Introduction
System identification is one of the most consolidated and relevant fields of study in science. One of the aims of this science is to obtain mathematical models analogous to the phenomena observed in nature. By analogous systems is meant a system capable of reproducing the characteristics observed in nature. With the identification of systems it is possible to model and investigate systems in an attempt to find some pattern in the observations [3, 1]. To identify a system, is necessary to assume a model capable of representing the linear and nonlinear characteristics of the system. Models are mathematical equations that try to describe an approximation of the real system. In the literature there are several ways to identify the same system [10]. Different mathematical and computational representations are used, it can be mentioned the neural networks, fuzzy logic, NARMAX (Nonlinear AutoRegressive Moving Average model with eXogenous input) models, among others. The representation of nonlinear systems can be obtained by means of the polynomial NARMAX [4]. The nonlinear NARMAX polynomials are linear in the parameters, which allows the use of parameter estimation algorithms for linear models [5]. This mathematical representation can be seen as a wellorganized recursive function in which the parameters are cautiously chosen [5].
In general, little attention has been given to the propagation of error in the area of system identification, especially in situations that present the polynomial NARMAX. One of the first works related to this subject was [8]. The authors have presented a theorem to estimate the lower bound error for the polynomial NARMAX. In that work, two pseudoorbits from two different interval extensionswere used to estimate the lower bound error. And in [9], it was found that the basin of attraction and the invariant distribution were not preserved, the authors shows that from two natural interval extensions may result in differents trajectories, using the lower bound error. Thus, this work aims to estimate the lower bound error for pseudoorbits derived from different interval extensions. The proposed method is applied in two identified models of the systems: sine map and DuffingUeda oscillator. Afterwards, the Lyapunov exponent, which is a parameter that characterizes the attractor dynamics. It measures the rate of divergence of neighboring orbits within the attractor, quantifying the dependence or sensitivity of the system relative to the initial conditions [11], is calculated and compares this result with the present in the literature for the method presented here and for the one proposed by [6].
The rest of the paper is organized as follows. In Section 2 we recall some preliminary concepts of lower bound error and representation the nonlinear systems. Then, in Section 3, we present the developed method. Section 4 is devoted to present the results, then the final remarks are given in Section 5.
2 Preliminary concepts
2.1 The polynomial NARMAX
The NARMAX model is a representation for nonlinear systems. This model can be represented as [4]
(1) 
where , e are, respectively, the output, the input and the noise terms at the discrete time . The parameters , e are their maximum lag. And is a nonlinear function of degree .
2.2 Recursive functions
In recursive functions is possible to calculate the state , at a give time, from an earlier state
(2) 
where is a recursive function and is a function state at the discrete time n. Given an initial condition , successive applications of the function it is possible to know the sequence . This sequence can be represented by and is defined as an orbit.
Using the computer to calculated the recursive functions, numeric errors are propagated during successive calculations, then the true orbit is not calculated but a representation of the same that is called pseudoorbit.
(3) 
where is an error and . So, we may define an interval associated with each value of a pseudoorbit Thus
(4) 
2.3 Natural interval extension
The natural interval extension is achieved by changing the sequence of arithmetic operation [7], that is, the extensions are mathematically equivalents.
Furthermore, two extension which algebraically is the same function may not be equivalent in interval arithmetic.
2.4 Lower bound error
The lower bound error was proposed by [8]. It is a practical tool capable of increasing the reliability of the computational simulation of dynamic systems.
Theorem 1.
Let two pseudoorbits and derived from two interval extensions. Let be the lower bound error of a map , then or .
The proof of this theorem can be found in [8].
3 Methods
This section is an extension of the work of [8]. The authors developed the Lower bound error theorem for two pseudoorbits from two different interval extensions. But, for a same map may exist more than two natural interval extensions, so the objective is to investigate the behavior of the natural interval extensions in the computer, exploring the effect of interval dependence, due to the repeated presence of a same interval variable in an algebraic expression, then check on the lower bound error of the pseudoorbits derived from the different natural interval extensions and calculate the Lyapunov exponent of the lower bound error from pseudoorbits and compare this value. It was clear that the function has more than two extensions, that is, can be rewritten in different ways.
The proposed method can be summarized in the following steps:

Choose the natural interval extensions;

Calculate the sequence of points of each system from the chosen extensions;

Determine the lower bound error from the combination of two functions;

Determine the maximum lower bound error;

Calculate the Lyapunov exponent and compare the result with the present in literature.
3.1 Generalization of the lower bound error
The generalization of the lower bound error is presented in the following theorem.
Theorem 2.
Let an arbitrary number of pseudoorbits derived from interval extensions. is the lower bound error, subject to , , and .
Proof.
The proof is conducted by reduction ad absurdum. Conversely, let the distance between two pseudoorbit given by and let us assume that it is possible to have a lower bound error described by
for all and . If it is true, considering the two pseudoorbits, let us say, and , for which we have maximum distance between them, it implies that which is a contradiction. And that completes the proof. ∎
It is clear that for the case of two pseudoorbits, this theorem is equivalent to that one proposed by [8].
4 Results
In this section, we present the lower bound error applied for two case studies, which exhibit nonlinear dynamics. We select some natural interval extensions that are equivalent. The two chosen models are for the systems sine map and DuffingUeda and all routines are performed in Matlab R2016a in a double precision. We used a computer with a processor Dual core @ 2.7GHz and a Windows 8.1 Professional operating system.
4.1 Sine map
A unidimensional sine map is defined as
(5) 
where . A polynomial NARMAX identified for this system is given by [10]
(6) 
Let us consider four equivalent interval extensions of the model 6:
(7)  
(8)  
(9)  
(10) 
Equations (7)(10) are mathematically equivalent, but they represent a different sequence of arithmetic operations. These extensions were simulated using as an initial condition. Figure 1(a) shows the evolution of the maximum lower bound error for the sime map and still shows the Lyapunov exponent associated with these extensions using the method developed in [6]. The literature indicates a Lyapunov exponent equals to 1.15 bits/s [10]. It is clear that the result presented in Figure 1(a) is in good agreement.
4.2 DuffingUeda
Considering a damped, periodically forced nonlinear DuffingUeda oscillator [3]:
(11) 
where is the cubic stiffness parameter, is a linear damping and is the amplitude of excitation. A polynomial NARMAX for the DuffingUeda oscillator was identified by [2].
(12)  
where , and . Let us consider four interval extensions of Eq. (12):
(13)  
(14)  
(15)  
(16)  
Figure 1(b) shows the evolutions of the maximum lower bound error for the DuffingUeda oscillator with the Lyapunov exponent associated. This largest Lyapunov was calculated by method developed in [6] with value 0.1202 for the maximum of the pseudoorbits. Oncemore, the computation are in good agreement with the values found in literature, which for this model was calculated in 0.115 [2].
5 Conclusion
The errors present in numerical simulations can lead to an erroneous result that does not correspond to the real situation of the problem. These errors can be of the representation of the model, of the insertion of data, of the numerical algorithm, due to the simplifications, of truncation, of rounding and among others. Thus, some methods have been proposed in order to measure these error, but they are complex from the computational point of view and from the mathematical approximation.
We presented a method to calculate a lower bound error for freerun simulation of the polynomial NARMAX. Our method is based on the comparison of pseudoorbits of the same models, but derived from different extension intervals. It makes use of recursive functions, which increases the relevance of this observation.
The methodology was applied in two cases, which are examples of identified systems obtained from literature, by means of the polynomial NARMAX. The sine map and Duffing Ueda oscillator are well known chaotic systems and have been identified using the polynomial NARMAX.
When we compare pseudoorbits that are equivalent from the point of view of interval analysis, we proved a theorem that the maximum distance of the pseudoorbits is greater than the distance of two pseudoorbits. This maximum value represents a small difference respect to the lower bound error for two pseudoorbits, but reduces the lower bound error. To prove this statement, the Lyapunov exponent was calculated for the maximum lower bound error, and in both models (sine map and DuffingUeda) the result was very close and in good agreement with values calculated in literature.
Acknowledgments
We thank CAPES, CNPq/INERGE, FAPEMIG and UFSJ for their support.
References
 [1] L. A. Aguirre, Introdução à identificação de sistemas  Técnicas lineares e nãolineares aplicadas a sistemas reais, Editora UFMG, Edição 4, (2015).
 [2] L. A. Aguirre and S. A. Billings, Validating identified nonlinear models with chaotic dynamics, International Journal of Bifurcation and Chaos, vol. 4, 109125, (1994).
 [3] S. A. Billings, Nonlinear system identification: NARMAX methods in the time, frequency, and spatiotemporal domains, John Wiley & Sons, (2013).
 [4] S. Chen and S. Billings, Representations of nonlinear systems: the NARMAX model, International Journal of Control, vol. 49, 10131032, (1989).
 [5] M. Korenberg, S. A. Billings, Y. P. Liu and P. J. Mcllroy, Orthogonal parameter estimation algorithm for nonlinear stochastic systems, International Journal of control, vol. 48, 193210, (1988), DOI: 10.1080/00207178808906169.
 [6] E. M. A. M. Mendes and E. G. Nepomuceno, A very simple method to calculate the (positive) Largest Lyapunov Exponent using interval extensions, International Journal of Bifurcation and Chaos, vol. 26, (2016), DOI: 10.1142/S0218127416502266.
 [7] R. E. Moore, R. B. Kearfott and M. J. Cloud, Introduction to interval analysis, Society for Industrial and Applied Mathematics, (2009).
 [8] E. G. Nepomuceno and S. A. M. Martins, A lower bound error for freerun simulation of the polynomial NARMAX, Systems Science & Control Engineering, vol. 4, 5058, (2016), DOI: 10.1080/21642583.2016.1163296.
 [9] E. G. Nepomuceno and E. M. A. M. Mendes, On the analysis of pseudoorbits of continuous chaotic nonlinear systems simulated using discretization schemes in a digital computer, Chaos, Solitons & Fractals, vol. 95, 2132, (2017), DOI: 10.1016/j.chaos.2016.12.002.
 [10] E. G. Nepomuceno, R. H. C. Takahashi, G. F. V. Amaral and L. A. Aguirre, Nonlinear identification using prior knowledge of fixed points: a multiobjective approach, Journal of Bifurcation and Chaos, vol. 13, 12291246, (2003), DOI: 10.1142/S0218127403007187.
 [11] A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Determining Lyapunov exponents from a time series, Physica D: Nonlinear Phenomena, vol. 16, 285317, (1985).