Identification of robotic systems with hysteresis using Nonlinear AutoRegressive eXogenous input models

Abstract

Identification of robotic systems with hysteresis is the main focus of this article. Nonlinear AutoRegressive eXogenous input models are proposed to describe the systems with hysteresis, with no limitation on the nonlinear characteristics. The article introduces an efficient approach to select model terms. This selection process is achieved using an orthogonal forward regression based on the leave-one-out cross-validation. A sampling rate reduction procedure is proposed to be incorporated into the term selection process. Two simulation examples corresponding to two typical hysteresis phenomena and one experimental example are finally presented to illustrate the applicability and effectiveness of the proposed approach.

Introduction

Hysteresis, a memory-dependent, multivalued relation between input and output, is often observed in many robotic systems. The system may exhibit a path-dependent pattern, where multiple outputs are associated with increasing or decreasing but the same input and form a loop under cyclic excitation. It exists in many applications, such as actuators and sensors involving smart materials (e.g. piezoelectrics ^1,2 and magnetostrictive materials ^3,4 ) which possess the property of hysteresis in the reaction, and some special robotic systems with hysteretic dynamics like aerial vehicles. ⁵ The control of these robots is difficult due to the presence of the high nonlinearity. Such nonlinearity turns to be a limitation of open-loop operations in high-precision applications, results in instabilities in closed-loop operation, and degrades the tracking performance even with the use of feedback control in tracking control application. ^6,7 It presents challenges in both analysis and controller design of robotic systems with hysteresis. A mathematical model, therefore, is required to predict and control the behavior of the robotic systems containing hysteresis.

Modeling of hysteresis has in recent years attracted increasing attention in various areas of robotics research, such as friction compensation, control of rubber tube actuator, elastic robot joints, and so on. Many researchers have studied this phenomenon, and many mathematical models have been developed to grasp the dynamic features of hysteresis phenomenon. ^{8 –11} Two of the most popular models are explained in detail in the following sections.

Preisach model

It was originally developed in the 1930s for magnetic hysteresis and is widely used to describe the hysteresis characteristics of smart materials ^{12 –14} in recent years. It has attracted considerable interest. Preisach model is a weighted superposition of simple independent delayed relays γ_αβ[u] with α, β corresponding to up and down switching values, respectively. ¹⁵ This model can be mathematically expressed as

$y(t)={\displaystyle {\iint }_{\alpha \ge \beta }\rho (\alpha ,\beta ){\gamma }_{\alpha \beta }[u](t)}d\alpha d\beta$

1

where y(t) and u(t) are the output and input at time t, respectively. ρ(α, β) is known as Preisach density (distribution) function, which weights the single relay units in the α, β - plane and defines the shape of hysteresis curve. The model is regarded as a general description of hysteresis phenomenon with two properties: the casual property and rate-independent property. That is, the output of the model has no relation with future inputs and does not depend on derivatives of the input.

Bouc–Wen model

It was originally proposed by Bouc (1967) and subsequently generalized by Wen (1976). The model has been widely used in the field of structural engineering, which provides reasonable accuracy in the deterministic and stochastic dynamic analyses. The Bouc–Wen model is formulated from mathematical analysis of the characteristic response properties. It is given by the following differential equation

$\dot{y}=\alpha \dot{x}-\beta \left|\dot{x}\right||y{|}^{\mu -1}y-\gamma \dot{x}|y{|}^{\mu }$

2

with $\alpha >0,\beta +\gamma >0,\beta -\gamma \ge 0,\mu >1$ . The shapes of hysteresis loops depend on the choice of the loop parameters. In general, α, β, γ influence the loop size, μ the smoothness. The Bouc–Wen model has many advantages. A notable advantage of this model is its capability to capture large numbers of patterns of hysteresis loops with various physical characteristics related to the hysteretic behavior, such as degradation of strength and stiffness, ¹⁶ pinching effect, ¹⁷ and asymmetry of the peak restoring force. ¹⁸ Another advantage of the model is its computational simplicity because only one auxiliary nonlinear differential equation is needed to describe the hysteretic behavior, and the model can be used to analyze the response of the robotic system under any excitation once the parameters have been identified.

Since the nonlinearity of hysteresis may show totally different properties for robotics of different areas, it is impossible to find one accurate model suitable for all types of robotic systems with hysteresis, which not only exists in robot sensors with smart materials but also robot dynamics. The hysteresis can be generally classified into two major categories: static hysteresis and dynamic hysteresis. The former is rate independent, which means there is no correlation between the behavior and the variation rate of input, like the phenomenon corresponding to the Preisach model; while the latter is rate dependent, which depends on the variation rate of input, like the features captured by the Bouc–Wen model.

The purpose of this article is to propose a generalized model with a clear structure in some straightforward way, which can be used to model robotic systems with hysteresis with no limitations on the nonlinear characteristics. The existing models are mostly continuous-time models, ^10,12,15 and the mathematical expressions are complicated, such as (1) and (2). One major shortcoming to the use of the continuous-time models is the difficulty of controller design due to the complicated form of the robotic model. Another shortcoming is lack of generalization capability. The nonlinearity and mechanism of hysteresis is different for robotics of different areas. The continuous-time models are given based on parameters with physical meaning, which analyze the mechanism for hysteresis of the particular robot. In addition, continuous-time model is not always available after theoretical analysis of the complex system. Then, discrete-time models can be considered. In practical application, discrete-time observations are obtained and the identification technique is realized by digital computer, so discrete-time models are more convenient in the identification procedure. The authors come up with the idea to use the discrete-time Nonlinear AutoRegressive with eXogenous input (NARX) model to describe the input–output relationship of robotic systems with hysteresis. The NARX model can exhibit a wide range of nonlinear behaviors with different properties such as chaos and bifurcations ¹⁹ and can be easily identified. The model structure is constructed in a linear-in-parameter form, which solves the difficulty of controller design.

The article is organized as follows: The NARX model is reviewed in the second section, and analysis for modeling robotic systems with hysteresis is made afterward. The third section introduces the orthogonal forward regression (OFR) algorithm for term selection. The fourth section is devoted to show the limitations of the OFR algorithm, and then an improved model selection procedure with new criteria is developed to solve the problems. Typical simulation examples are discussed in the fifth section, together with the detailed derivations and performance analysis. Experimental example of unmanned aerial vehicle is given in the sixth section. Finally, some concluding remarks are drawn and the limitations of potentially applying the polynomial NARX model to hysteresis identification are discussed in the seventh section. For the sake of easy implementation of digital computers, all the signal processing derived here are based on the discrete-time case.

The NARX model

The NARX representation has attracted considerable interest in modeling nonlinear systems, and many relevant analysis tools and identification algorithms have been developed in recent years. ^{20 –23} The NARX model is an extension of the linear ARX model. The AR model is used when current output is dependent only on the previous outputs, and the ARX model is used when there is exogenous input given to the AR model, as shown in Figure 1.

Figure 1.

The ARX model.

The NARX model is defined as

$y[n]=f(y[n-1],\dots ,y[n-n_y],u[n-1],\dots ,u[n-n_u])+e[n]$

3

where y[n] and u[n] are the output and input of the system, respectively; n_y and n_u are the maximum lags for system output and input, respectively; f(⋅) is a nonlinear function which needs to be identified from given observed data; e[n] is the prediction error, which is thought to be a zero mean noise sequence when f(⋅) gives the reasonable description of the nonlinear relation between the output and input. As mentioned in the Introduction, hysteresis is a multivalued and nonsmooth relation between input and output; however, when the input is expanded from u[n] to x [n], where $\mathit{x}[n]=[y[n-1],\dots ,y[n-n_y],u[n-1],\dots ,u[n-n_u]]$ , corresponding to the definition of NARX model with high-dimensional input spaces, the relationship f(⋅) will turn to be a smooth single-valued mapping, which brings a lot of convenience in the determination of the nonlinear function.

The smooth single-valued function f(⋅) is often constructed by a linear-in-parameter form, using a variety of basic functions ${\phi }_i(\cdot )$ , which can be expressed in the regression form

$y[n]={\displaystyle \sum _{i=1}^{n_m}{\theta }_i{\phi }_i(\mathit{x}[n])+e[n]}$

4

where $\mathit{x}[n]=[y[n-1],\dots ,y[n-n_y],u[n-1],\dots ,u[n-n_u]]$ , ${\theta }_i(i=1,2,\dots ,n_m)$ are unknown parameters, and n_m is the number of model terms potentially involved. The linear-in-parameter model is widely applied for system identification in many industrial areas, because simple algorithms like least squares method can be used for parameter estimation. If φ_i(⋅) is a polynomial function, the model can be given as ²⁴

$y[n]={\displaystyle \sum _{m=0}^l{\displaystyle \sum _{p=0}^m{\displaystyle \sum _{n_1=1}^{n_y}\dots }}}{\displaystyle \sum _{n_m=0}^{n_u}{\theta }_{p,m-p}[n_1,\dots ,n_m]}\times {\displaystyle \prod _{i=1}^{p}y[n-n_i]{\displaystyle \prod _{i=p+1}^{m}u[n-n_i]+e[n]}}$

5

where l is the degree of the polynomial model defined as the maximum order of model terms. Corresponding to (4), φ(⋅) is thus of the form

$\phi (\cdot )={\displaystyle \prod _{i=1}^{p}y[n-n_i]{\displaystyle \prod _{i=p+1}^{m}u[n-n_i]}}\left(p=0,1,\dots m;n_i=\{\begin{array}{cc}1,2,\dots ,n_y& i\le p\\ 0,1,\dots ,n_u& i>p\end{array};m=0,1,\dots ,l\right)$

6

Most nonlinear models (e.g. the Volterra series model) only consider the input data to the system as the input to the model. The amount of model terms is sometimes very large; otherwise, the model cannot guarantee the approximation accuracy. The NARX model, however, adds the historical output data to the input variables of the model. The model terms decrease significantly, because the historical output information contains some nonlinear characteristics of the system. So, the NARX model can give a simpler structure for the complex systems. In continuous-time models (1) and (2), integral and differential calculus is needed which adds noise to the identification procedure, as a result algorithm complexity is increased to ensure the accuracy. The NARX model avoids the problem. The discrete-time data can be used directly to approximate the continuous-time system with good performance. The NARX model decomposes the system into a linear segment and a nonlinear segment. The linear segment represents the influence of historical data which deals with problem of rate dependence, while the nonlinear segment (e.g. polynomial form) corresponds to the static nonlinear mapping. The clear structure makes it easier for control design.

OFR algorithm

From (6), it can be easily found out that the total number of terms increases rapidly with the maximum lag $n_u,n_y$ and degree l. In most application situations, however, only a small percent of the total terms are confirmed to be significant to the performance of the model, so it is necessary to find a method to select model terms. On one hand, redundant terms must be abandoned to avoid overfitting, on the other hand, the model should be as simple as possible, on the basis of containing key model characteristics of the robotic system. Model structure selection turns out to be a key task in the identification process. The main step of the method is to define a criterion for indicating the significance of each term φ(⋅). Several criteria have been proposed in the literature for NARX models. ²⁵ One of the most widely used is the error reduction ratio (ERR) based on the OFR algorithm. ²⁶ The authors of this article also use it as a reference of the final chosen term detection method. The vector Φ _j is constructed using the result of the jth model term ${\phi }_j(\cdot )$ from given data at each time. Then the jth error reduction ratio, ERR _j (also called the squared correlation coefficient), is defined as

${\text{ERR}}_j=C\left(\mathit{Y},{\Phi }_j\right)=\frac{{〈\mathit{Y},{\Phi }_j〉}^2}{〈\mathit{Y},\mathit{Y}〉〈{\Phi }_j,{\Phi }_j〉}=\frac{{\left({\mathit{Y}}^{\text{T}}{\Phi }_j\right)}^2}{\left({\mathit{Y}}^{\text{T}}\mathit{Y}\right)\left({\Phi }_j{}^{\text{T}}{\Phi }_j\right)}=\frac{{\left({\displaystyle {\sum }_{i=1}^N{y}_i{\phi }_j^i}\right)}^2}{{\displaystyle {\sum }_{i=1}^N{y}_i{}^2{\displaystyle {\sum }_{i=1}^N{\left({\phi }_j^i\right)}^2}}}$

7

The ratio provides an effective means to measure the dependency between the output and each term of the model. Then the significant terms can be selected out gradually based on the OFR algorithm. The procedure is described as follows.

Step 1. Select out the term with the largest ERR.

$l_1=\arg \underset{1\le j\le n_m}{\max }\left\{C\left(\mathit{Y},{\Phi }_j\right)\right\}$

8

Then the first significant term can be selected as $w_1={\phi }_{l_1}$ .

Step j. Let r _j represent the residual output vector of the model in the jth step. It is given by

${\mathit{r}}_m={\mathit{r}}_{m-1}-\frac{{\mathit{r}}_{m-1}{\mathit{w}}_{m-1}}{{\mathit{w}}_{m-1}^T{\mathit{w}}_{m-1}}{\mathit{w}}_{m-1}$

9

In the first step, r ₀ = Y . Select out the term with the largest ERR in the remaining terms.

$l_j=\arg \underset{i\ne l_k(1\le k\le j-1)}{\max }\left\{C\left({\mathit{r}}_j,{\Phi }_i\right)\right\}$

10

Then the jth significant term can be selected as $w_j={\phi }_{l_j}$ .

The procedure terminates at the Mth step. M is determined by the terminating condition. For simple situations, it can be given as

$1-{\displaystyle \sum _{i=1}^M{\text{ERR}}_i}<\rho$

11

where ERR _i equals $C\left({\mathit{r}}_i,{\Phi }_{l_i}\right)$ derived in each step, and ρ is the desired error tolerance.

Extension to the OFR algorithm

Limitations of the OFR identification approach

As shown in the third section, the form only considers the tolerance in the identification procedure when all the data are used for fitting. The model usually shows a bad performance when used on future data sets. ²⁷ However, in order to be generally used, a model must have good extrapolation properties. It is necessary to split the data into two subsamples: a fitting sample and a validation sample. The fitting sample is used to measure the significance of model terms, while the validation sample is used to terminate the procedure. The fitting sample must contain the key characteristics of the robotic system to ensure that the model identified using the data has the ability to perform the hysteretic behavior with expected properties.

Since the NARX model takes the output history as part of variables in the model, it results in some drawbacks. The model works only when the historical outputs of the robotic system are available. If the output of the robotic system is not measurable in the procedure, the model terms that contain the output need to be calculated by the estimated result of the previous steps, then the system error will be accumulated, leading to inaccuracy in the model.

The data used for identification procedure is obtained from the continuous-time system through periodic sampling

$u[n]=u_c(nT),y[n]=y_c(nT)$

12

where u_c (t), y_c (t) are the continuous-time input and output signals, respectively, and T is the sampling period. T has a tremendous influence on the effectiveness of the final model. The coefficients of the overall model are dependent on the choice of the sampling rate. Practical experience has shown that the mainly influence of sampling rate on modeling is the coefficients of model terms but not the model structure apart from the effects of oversampling. ²⁸ Data oversampling will bring numerical problems for nonlinear model structure selection. When the T is chosen to be too small, the fitting data will be so intensive that x [n] ( $\mathbf{x}[n]=[y[n-1],\dots ,y[n-n_y],u[n-1],\dots ,u[n-n_u]]$ ) are highly correlated, which will cause potential problems in distinguishing the significance of model terms. In particular, the output can be extrapolated by the outputs of earlier time as

$\begin{array}{l}y[n]=y[n-1]\text{\hspace{1em}or}\\ y[n]=2y[n-1]-y[n-2]\end{array}$

13

Then the model finally identified is always a combination of expressions in (13) with a random rule. As a result of that, the model only displays the linearized characteristics at every point, incapable of representing the nonlinear properties of robotic systems with hysteresis. Even if the final model doesn’t end up in the linear form for T_s is small, the ERR of term y[n − 1] is close to 1, while the terms determined by input history, like u[n − 1], u[n − 2], have extremely low ERRs, which will lead to wrong results in the following steps of the OFR procedure after y[n − 1] is selected, since the precision of orthogonalization is influenced profoundly by the noise of data, owing to the strong correlation between model terms. The sampling period, however, cannot be large, in case some important nonlinear information is missed resulting in low approximation accuracy to continuous-time system, especially when rapid change exists in the robotic system.

The PRESS statistic

A commonly used criterion in data splitting is the prediction error sum of squares (PRESS). ²⁹ The procedure of leave-one-out cross validation goes like this: remove one observation at a time, use the removed observation as validation point and the remaining N − 1 observations as fitting sample, then estimate the coefficients and evaluate the deleted response ${\widehat{y}}_n^{(-i)}[i]$ from the estimated model at x = x _i , repeat this process on all points and finally get PRESS residual, defined as

$\begin{array}{l}\text{PRESS}[n]={\displaystyle \sum _{i=1}^N{\left[y[i]-{\widehat{y}}_n^{(-i)}[i]\right]}^2}\\ ={\displaystyle \sum _{i=1}^N{\left[{\epsilon }_n^{(-i)}[i]\right]}^2}\end{array}$

14

where ${\epsilon }_n^{(-i)}[i]$ is the predicted residual evaluated at the ith point with the fitting sample of size N − 1. Obviously, the computation is complex for the model is to be fitted for N times. Researchers have done much work to simplify the procedure. When using the relationship between the PRESS residual and the ordinal residual, the PRESS reduces to

$\text{PRESS}[n]={\displaystyle \sum _{i=1}^N{\left[\frac{\epsilon [i]}{1-h_{ii}}\right]}^2}$

15

where h_ii represents the prediction variance, given by

$h_{ii}={\mathit{x}}_i\prime {(X\prime X)}^{-1}{\mathit{x}}_T$

16

An important property of h_ii can be easily derived as

$\begin{array}{c}{\displaystyle \sum _{i=1}^N{h}_{ii}}={\displaystyle \sum _{i=1}^N{\mathit{x}}_i\prime {(X\prime X)}^{\text{T}}{\mathit{x}}_i}\\ ={\displaystyle \sum _{i=1}^N\mathrm{tr}{\mathit{x}}_i\prime {(X\prime X)}^{\text{T}}{\mathit{x}}_i}\\ ={\displaystyle \sum _{i=1}^N\mathrm{tr}{\mathit{x}}_i{\mathit{x}}_i\prime {(X\prime X)}^{\text{T}}}\\ =\mathrm{tr}({\displaystyle \sum _{i=1}^N{\mathit{x}}_i{\mathit{x}}_i\prime ){(X\prime X)}^{\text{T}}}\\ =n\end{array}$

17

When N >> n, ³⁰

$\text{PRESS}[n]\approx {\displaystyle \sum _{i=1}^N\frac{{\epsilon }^2[i]}{(1-n/N)}}$

18

This approximation significantly reduces computational complexity of PRESS.

Considering the problems for model structure selection, a model identified using a finite data set may not have good performance over the fitting data, so the measure of model accuracy has to depend on an additional data set. The cross-validation way helps a lot for the model generalization when substituting the terminating condition of expression (11) with the PRESS statistic in the ORF algorithm. Moreover, the estimated function $\widehat{f}$ of (3) provides higher accuracy as the complexity of $\widehat{f}$ increases, which is mainly determined by n, the total number of the polynomial NARX model terms. This may cause overfitting to the noise in y[n]. From the simplified expression of PRESS in (18), we can see that the criterion avoids the overfitting phenomenon efficiently because the value of PRESS will increase as n increases owing to the existence of −n in the denominator.

Sampling rate reduction

Following the discussion in “The PRESS statistic” section, the sampling rate must be determined under full consideration of the dynamic characteristics of the data to guarantee that the OFR procedure can capture the main nonlinear effects of the robotic system. When the sampling period is too small, the nonlinear effects will appear to be local linearization, and the OFR procedure will be disturbed by the form of (13). So in situations of very small sampling period, the authors suggest using the improved algorithm followed by sampling rate reduction (SRR) procedure. A decimator in Figure 2, that is, a system with a low-pass filter followed by compression, is required for SRR by integer factor L.

Figure 2.

Decimator for sampling rate reduction by L.

The maximum relative deviation (MRD) and the maximum relative error (MRE) are is proposed as the criteria for choice of L in this article. For the new discrete-time data after decimation, with the new sampling period of T′ (=LT), MRD is given as

$\text{MRD}[T\prime ]=\underset{1\le k\le N/L}{\max }\left\{\frac{\left(y\prime [k+2]+y\prime [k]\right)/2-y\prime [k+1]}{y\prime [k+1]}\right\}$

19

where y′[k] = y[Lk] is the output of the new discrete-time data. Then the inferior limit of L can be determined by

$1-\text{MRD}[T\prime ]<\alpha$

20

where 0 < α < 1 is the expected confidence level. If MRD doesn’t satisfy the inequality, the output y′[n] can be directly derived from the past two points y′[n − 1] and y′[n − 2], that is, the linear relationship in (13) can already reach the requirement of modeling accuracy. The nonlinear information, thus, is not necessary to be contained in the model. However, models which capture little nonlinear characteristics cannot present the nonlinear properties of the robotic system and, certainly, have no significance for controller design in practice.

Then the procedure with SRR is given. First, evaluate the sampling rate 1/T of the original data set D ₀ by the value of MRD. If the value is smaller than 1 − α or very close to it, we may use the following procedure to handle the oversampling data. The data rate is reduced by increasing integer factor L_i (L_i = 2, 3, 4,…) using a decimator, and the value of MRD _i for each new data set D _i with the data rate 1/T_i (T_i = L_iT) is calculated. Note that the MRD _i value is getting bigger as L_i increases, which means local linear regression has less influence on model fitting, whereas some important nonlinear characteristics may be abandoned in the remaining data by the decimator when L_i increases to some extent. So the authors propose the MRE to measure the importance of the information missed due to decimation by L_i , defined as

$\text{MRE}\left[L_i\right]=\{\begin{array}{ll}\underset{1\le k\le N-L_i}{\max }\left|\frac{\left(y[k]+y\left[k+L_i\right]\right)/2-y\left[k+L_i/2\right]}{y\left[k+L_i/2\right]}\right|\hfill & \hfill L_i=2,4,6,\dots \\ \underset{1\le k\le N/L_i-1}{\max }\{\left|\frac{\left(y[k]+y\left[k+L_i\right]\right)/2-\left(y\left[k+(L_i+1)/2\right]+y\left[k+(L_i-1)/2\right]\right)/2}{\left(y\left[k+(L_i+1)/2\right]+y\left[k+(L_i-1)/2\right]\right)/2}\right|,\hfill & L_i=3,5,7,\dots \hfill \end{array}$

21

where y[k] is the output of original data set D ₀. The superior limit of decimation factor L is determined by

$1-\text{MRE}\left[L_i\right]>\alpha$

22

where 0 < α < 1 is the expected confidence level. We can choose an appropriate value of L based on the two restrictions in (20) and (22). Then the sampling rate is adjusted to 1/(LT). The model structure is finally selected by the OFR algorithm with the PRESS statistic following this SRR procedure.

Simulation studies

This section investigates the efficiency and performance of the polynomial NARX model for the identification of robotic systems with hysteresis, by applying the OFR algorithm with the PRESS statistic following SRR procedure to two typical examples. The first example is a static hysteresis, given by a simulated Preisach model, while the latter is a dynamic hysteresis, given by a simulated Bouc–Wen model.

A simulated Preisach model

Consider a Preisach model described by the expression as follows

$y(t)={\displaystyle \underset{\alpha \ge \beta }{\iint }\rho (\alpha ,\beta ){\gamma }_{\alpha \beta }[u](t)}d\alpha d\beta +\xi (t)$

23

where ${\gamma }_{\alpha \beta }[u]$ is assumed to be bounded as Figure 3(a), α, β - plane is given as Figure 3(b), ρ(α, β) is assumed to be constant, namely, given as a uniform distribution, and ξ(t) is a Gaussian white noise of zero mean and variance σ² = 0.01. The model is simulated by setting the input signal u(t) as an increasing sequence from −0.5 to 0.5 followed by a decreasing sequence from 0.5 to −0.5 with period 0.01 s, and 202 input–output data point are collected.

Figure 3.

Elements for model in (23): (a) ${\gamma }_{\alpha \beta }[u]$ and (b) α, β - plane.

The model terms φ(⋅) in (4) are chosen to be polynomial functions in (6) determined by the following element: $\mathbf{x}[n]=[y[n-1],y[n-2],y[n-3],u[n],u[n-1],u[n-2],u[n-3]]$ , which is the input to the NARX model. First, the SRR procedure is used to adjust the data rate. The values of MRD and MRE in the decimation process are summarized in Table 1.

Table 1.

Evaluation for L in the SRR procedure.

Li	MRD	MRE	Li	MRD	MRE
2	0.0129	0.0028	9	0.0425	0.0296
3	0.0155	0.0050	10	0.0480	0.0343
4	0.0192	0.0091	11	0.0534	0.0383
5	0.0233	0.0124	12	0.0589	0.0429
6	0.0276	0.0170	13	0.0631	0.0467
7	0.0324	0.0208	14	0.0711	0.0510
8	0.0372	0.0256	15	0.0745	0.0547

SRR: sampling rate reduction; MRD: maximum relative deviation; MRE: the maximum relative error.

The confidence level is usually expected to be 95% in practice, corresponding to α in expressions (20) and (22), then the decimation factor L can be set to be 11, 12, and 13 according to the results of the SRR procedure. The PRESS values (Table 2), by the OFR algorithm, with different L, over the data set of different size accordingly, are calculated for different model length n.

Table 2.

PRESS values for different model length.

n	PRESS (L = 11)	PRESS (L = 12)	PRESS (L = 13)
6	0.0114	0.0087	0.0062
7	0.0107	0.0085	0.0063
8	0.0106	0.0080	0.0060
9	0.0114	0.0080	0.0059
10	0.0117	0.0088	0.0062
11	0.0126	0.0096	0.0064
12	0.0131	0.0109	0.0081
13	0.0149	0.0137	0.0110

PRESS: prediction error sum of squares.

The PRESS statistic suggests choosing 8 model terms over training data with L = 11 in the SRR procedure, 8 or 9 model terms for L = 12, and 9 model terms for L = 13, and the model terms selected out are almost the same as shown in Table 3.

Table 3.

The term selection results based on the data with different L.

Iter.	Regressors	Parameters (L = 11)	Parameters (L = 12)	Parameters (L = 13)
1	$y^2(t-1)$	0.0631	0.0657	−0.2019
2	u(t)	0.3232	0.3627	0.3871
3	$u(t-3)$	−0.1703	−0.1924	−0.1921
4	$u^2(t)$	−0.3914	−0.0936	−0.1171
5	$u(t)u(t-3)$	−0.2066	0.1854	0.2293
6	$u^2(t-3)$	–	−0.0967	−0.1056
7	$u^3(t)$	0.5929	−0.0513	−0.0758
8	$u^2(t-1)$	0.1577	–	0.0433
9	$y(t-3)$	0.9664	1.0881	1.1669
Average relative error		0.0015	0.0014	3.3633E−4

The effect of data rate in modeling the robotic systems with hysteresis is sufficiently shown in the results in Table 3. The parameters in the model, but not the model structure, have the strong dependence on the sampling rate. The process also verifies the significance of the SRR procedure and finally shows a good performance for modeling the simulated Preisach model using a polynomial NARX model, which proves that the proposed method is effective for modeling robotic systems with static hysteresis like magnetostrictive actuators.

A simulated Bouc–Wen model

Consider a Bouc–Wen model described in (2), where $\alpha =1,\beta =1,\gamma =1,\mu =4.$ The model is simulated by setting the input signal u[n] as a sinusoidal wave and 101 input–output data point are collected. The SRR procedure is not needed, since the sampling rate is low, and the value of MRE doesn’t satisfy the inequality of (22) at L = 2. The leave-one-out cross-validation is carried out directly without the SRR procedure. The input to the model is chosen to be $\mathbf{x}[n]=[y[n-1],y[n-2],y[n-3],u[n],u[n-1],u[n-2],u[n-3]]$ , and the polynomial degree is chosen to be l = 3. The result of the PRESS statistic is shown in Figure 4.

Figure 4.

The PRESS statistic versus the model length. PRESS: prediction error sum of squares

It is clear that the model length is best to be 11. The model by choosing 11 model terms is finally identified as

$y[n]=0.3236y^3[n-1]-0.4363u[n-3]+0.0183u^2[n-3]y[n-3]-0.2970u[n-2]+0.6885y[n-3]-0.1873u[n-1]+0.1577u^2[n-2]y[n-2]+0.0251u^3[n-3]+0.7538u[n]-0.1152u^2[n-1]y[n-1]-0.0148u[n-1]y^2[n-1]$

24

The model structure identified by the OFR algorithm based on the PRESS statistic is in agreement with traditional modeling assumptions. Figure 5 displays the performance of (24) over the test data set consisting of 120 points. The continuous line is the real values from the simulated Bouc–Wen model, while the dashed line is the estimated output from the NARX model in (24). The relative error of each point is calculated to test the model validity, and the results are given in Figure 6, where the two horizontal lines indicate the desired error tolerance of 5%. The jumping phenomena can be ignored directly, because it is the result of zero values of y[n] in the denominator. Clearly, the model validity tests are well satisfied.

Figure 5.

Model performance on the test data.

Figure 6.

Model validity tests for (24).

Experimental result

This section aims at illustrating the effectiveness of the NARX model and the accuracy of above explained identification method. For this purpose, we employ an unmanned aerial vehicle for research. The mass and geometric characteristics are shown in Table 4.

Table 4.

Mass and geometric characteristics of the model.

Parameter	Value
Mass	≤8 kg
Length	1.1825 m
Wingspan	0.8475 m
Wing area	0.304688 m2
Mean aerodynamic chord	0.4269 m
Center of Gravity (CG) location	33.37 (% m.a.c.)

In general situations of steady flights with small angle of attack, the models for aerodynamic forces and moments are linear to the state variables. However, for the flights with large angles of attack and sideslip or high angular rates, the forces and moments demonstrate hysteresis effects. The test data of longitudinal large-amplitude sinusoidal motions is investigated. Preisach model cannot apply to this situation, because the hysteresis is relevant to the angular rate. Bouc–Wen model is also inappropriate, because the differential operator on aerodynamic forces and moments will lead to an accumulation of errors. In the study by Greenwell, ⁵ a reduced-frequency model is discussed, which was proposed for modeling the hysteresis nonlinearity in aerial vehicle. Comparison of the NARX model and the reduced-frequency model is shown in Figure 7 based on the test data of the normal force coefficient C_N and the angle of attack α with the frequency of 0.4 Hz. The NARX model and the proposed identification approach can be applied to the unmanned aerial vehicle with a high precision.

Figure 7.

Comparison result.

Conclusions

The polynomial NARX model has been considered for modeling robotic systems with hysteresis. The model term selection problem has been investigated for using polynomial NARX models. A critical analysis of the standard OFR algorithm has shown some limitations, particularly when the sampling rate is high. The terminating condition for the OFR algorithm has been modified. The sampling rates over a reasonable range affect the parameter estimates but not affect the model structure. A SRR procedure is proposed based on calculating the MRD and MRE. The applicability and effectiveness of the SRR procedure for term selection have been demonstrated by two simulated examples and one experimental example. Another important property of using polynomial NARX models for hysteresis identification is that few assumption or priori knowledge about the robotic system is needed. When the robotic system is appeared to have a complex noise model, the NARX models should be extended to the nonlinear autoregressive moving average with exogenous variables models.

As a final remark, it should be pointed out that the proposed modeling approach for robotic systems containing hysteresis is not viable for the situation that the previous output data are not available in the process. One solution to this problem is to use the estimated values from the previous steps. However, the solution will cause a stepwise accumulation of errors, owing to the sensitivity of the iteration process to initial condition, so a compensation measure is needed to guarantee the accuracy.