SysIdentPy: A Python package for System Identification using NARMAX models

The field of System Identification (SI) aims to build mathematical models for static and dynamic behavior from experimental data (Ljung, 1987). In particular, nonlinear system identification has become a central issue in the SI community, and from the 1950s onwards many methods have been proposed. In this respect, NARMAX (Nonlinear AutoRegressive Moving Average with eXogenous input) models are among the most well-documented and used model representation of dynamical systems (Billings, 2013).


Summary
The field of System Identification (SI) aims to build mathematical models for static and dynamic behavior from experimental data (Ljung, 1987). In particular, nonlinear system identification has become a central issue in the SI community, and from the 1950s onwards many methods have been proposed. In this respect, NARMAX (Nonlinear AutoRegressive Moving Average with eXogenous input) models are among the most well-documented and used model representation of dynamical systems (Billings, 2013).
The NARMAX model was proposed by (Billings & Leontaritis, 1981;Chen & Billings, 1989;Leontaritis & Billings, 1985) and can be described as where n y ∈ N * , n x ∈ N, n e ∈ N , are the maximum lags for the system output and input respectively; x k ∈ R nx is the system input and y k ∈ R ny is the system output at discrete time k ∈ N n ; e k ∈ R ne represents uncertainties and possible noise at discrete time k. In this case, F is some nonlinear function of the input and output regressors and d is a time delay typically set to d = 1.
Although there are many possible approximations of F(·) (e.g., Neural Networks, Fuzzy, Wavelet, Radial Basis Function), the power-form Polynomial NARMAX model is the most commonly used (Billings, 2013;Khandelwal, Schoukens, & Tóth, 2020): where p is the number of regressors, Θ i are the model parameters, and a i , m, b i , j and d i , l ∈ N are the exponents of the output, input and noise terms, respectively.
The following example is a polynomial NARMAX model where the nonlinearity degree is equal to 2, identified from experimental data of a DC motor/generator with no prior knowledge of the model form, taken from (Lacerda Junior, Almeida, & Martins, 2017): The Θ values are the coefficients of each term of the polynomial equation. Polynomial basis functions are one of the most used representations of NARMAX models due to several interesting atrributes, such as (Aguirre, 2004;Billings, 2013): • All polynomial functions are smooth in R.
• The Weierstrass approximation theorem (Weierstrass, 1885) states that any continuous real-valued function defined on a closed and bounded space [a, b] can be uniformly approximated using a polynomial on that interval. • They can describe several nonlinear dynamical systems (Billings, 2013), including industrial processes, control systems, structural systems, economic and financial systems, biology, medicine, and social systems (Aguirre, 2004;Billings, 2013;Boynton, Balikhin, Wei, & Lang, 2018;Chiras, Evans, & Rees, 2001;Fung, Wong, Ho, & Mignolet, 2003;Guo, Guo, Billings, & Wei, 2016;Kukreja, Galiana, & Kearney, 2003;Lacerda Junior, Martins, Nepomuceno, & Lacerda, 2019;Martins & Aguirre, 2016). • Several algorithms have been developed for both structure selection and parameter estimation of polynomial NARMAX models. • Polynomial NARMAX models can be used both for prediction and inference. The structure of polynomial NARMAX models are easy to interpret and can be related to the underlying system, which is not a trivial task when using, for example, neural or wavelet functions.
Estimating the parameters of NARMAX models is a simple task if the model structure is known a priori. However, usually there is no information on what terms one should include in the final model, and selecting the correct terms has to be part of the system identification procedure. Thus, the identification of NARMAX models is twofold: selecting the most significant regressors given a dictionary of candidate terms, which relies on model structure selection algorithms, and estimating their parameters.

SysIdentPy
SysIdentPy is an open-source Python package for system identification using polynomial NARMAX models. The package can handle SISO (Single-Input Single-Output) and MISO (Multiple-Inputs Single-Output) NARMAX model identification and its variants such as NARX, NAR, ARMAX, ARX, and AR models. It provides various tools for both model structure selection and parameter estimation including classical algorithms, e.g., forward regression orthogonal least squares and extended least squares orthogonal forward regression; parameter estimation using ordinary least squares, recursive algorithms and adaptative filters; the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), Khinchin's law of iterated logarithm criterion (LILC), and Final Prediction Error (FPE) methods for model order selection (Haber & Keviczky, 1999); regression metrics; and residual analysis. The reader is referred to the package documentation for further details.
SysIdentPy is designed to be easily expanded and user friendly. Moreover, the package aims to provide useful tools for researchers and students not only in the SI field, but also in correlated areas such as Machine Learning, Statistical Learning and Data Science. Recently, an R package was published (Ayala, Gritti, & Santos Coelho, 2020) with tools to model dynamic systems using NARMAX models. However, to the best of our knowledge, SysIdentPy is the first open-source package for system identification using NARMAX models in Python. Moreover, SysIdentPy includes recursive and gradient methods for parameter estimation, e.g., recursive least squares, affine least mean squares, sign-sign least mean squares and many others that are not available in the above-mentioned R package. Also, the user can choose between four methods for model order selection, which is not possible with the mentioned R package.

Example
The following is an example of how to use SysIdentPy to build a NARMAX model from data. For simplicity, the example uses simulated data with 1000 samples, generated using the method get_miso_data: import numpy as np import pandas as pd import matplotlib.pyplot as plt from sysidentpy.polynomial_basis import PolynomialNarmax from sysidentpy.metrics import root_relative_squared_error from sysidentpy.utils.generate_data import get_miso_data x_train, x_valid, y_train, y_valid = get_miso_data(n=1000, colored_noise=False, sigma=0.001, train_percentage=90) Assuming that there is no information regarding what system generated the data, a dictionary of candidate terms must be created by defining the nonlinearity degree of the polynomial function and the maximum lag of the input and output terms. These parameters are, respectively, non_degree, ylag, xlag. The Akaike Information Criterion is chosen for model order selection and the least squares method is used for parameter estimation: The user can also run a SISO example by replacing get_miso_data with get_siso_data and the xlag values with an integer or a list of integers. If one wants to estimate the parameters using, for example, the recursive least squares algorithm, just set estimator to 'recursive_least_squares'. Replacing the AIC method with BIC, for example, can be done analogously by replacing 'aic' with 'bic'.
The fit method is used to obtain the model and predict to validate the model using new data. The metric to evaluate is the relative root squared error. To get the root mean square error metric, for example, import it using from sysidentpy.metrics import root_mean _square_error and replace the root relative squared error method with it. The table below and Figure 1 are the ouput of the aforementioned example. Table 1 details the regressors chosen to compose the final model, its respective parameters and the error reduction ratio (ERR), which measure the contribution of each regressor to explain the system output. ERR values can be interpreted as a feature importance metric. Figure 1 depicts the simulation of model prediction and the validation data as well as the autocorrelation of the model residues and the cross-correlation between the input and the residues.
For more information and examples of how to build NARMAX models and its variants using different methods for parameters estimation, model order selection and many more, see the package documentation.

Future work
Future releases will include new methods for model structure selection of polynomial NAR-MAX models, new basis functions, multiobjective model structure selection and parameter Lacerda et al., (2020 estimation algorithms, new adaptative filters, frequency domain analysis, and algorithms for using NARMAX models for classification problems.