Novel recursive least squares identification for a class of nonlinear multiple-input single-output systems using the filtering technique

Abstract

In this article, we consider the identification problem of a class of nonlinear multiple-input single-output output-error autoregressive systems. First, a recursive generalized least squares algorithm using the auxiliary model identification idea is developed. Then, using the filtering technique, the identification model is decomposed into a filtered sub-identification model and a noise sub-identification model. For solving the difficulties that the filter is unknown and the information vectors contain the unknown variables, the interactive estimation theory and the the idea of replacing the unknown variables with their corresponding estimates are employed: the recursive least squares method is again used for identifying the system and noise model parameters, and the parameter estimates of the noise model are used to construct the estimated filter. Finally, a nonlinear example is given to verify the effectiveness of the algorithms, and the simulation results show that the recursive least squares algorithm using the filtering technique can produce more accurate parameter estimates under larger noise variances.

Introduction

Nowadays, a wide range of recursive parameter estimation methods have been developed for system identification such as recursive least squares methods,^1,2 gradient methods,^3,4 Kalman filtering methods,^5–7 and Newton methods.^8–10 Compared with the gradient methods, recursive least squares methods take advantage of fast convergence rates, but when the system is contaminated by the colored noises, the performance of these methods, such as estimation accuracy, will decline if the colored noises are not well handled. One way to improve the estimation accuracy under the noise environment is to include the noise information into the information vector. Based on that idea, the recursive extended least squares algorithms and the recursive generalized least squares algorithms have been developed for identifying the parameters of linear systems and nonlinear systems. For example, Wang ¹¹ presented a hierarchical parameter estimation algorithm for a class of multiple-input multiple-output (MIMO) Hammerstein systems based on the reframed models, and Wang and Zhang ¹² improved the least squares identification algorithm for multivariable Hammerstein systems. Hu et al. ¹³ proposed a recursive extended least squares algorithm for Wiener nonlinear systems with moving average noises.

The filtering technique is effective in noise reduction and has been applied in various areas, ¹⁴ including fault detection, ¹⁵ signal processing,^16–18 acoustic echo cancellation, ¹⁹ and filter design.^20–22 Recently, a multi-innovation stochastic gradient identification algorithm has been presented for Hammerstein-controlled autoregressive autoregressive systems based on the filtering technique; ²³ an auxiliary model-based least squares algorithm has been presented for a dual-rate state-space system with time-delay using the data filtering. ²⁴

Multiple-input single-output (MISO) systems are a special class of multivariable systems,^25,26 which have various applications in the industrial processes, communication systems,^27,28 and so on. In the literature, in order to improve the convergence rate, Liu developed a stochastic gradient algorithm for MISO systems using the multi-innovation theory. ²⁹ The least square–based algorithms for identifying MISO systems can be found in Zhang. ³⁰ In this work, we extend the previous work from considering the identification problems for a class of single-input single-output (SISO) linear-in-parameter systems^31,32 to MISO linear-in-parameter systems with the autoregressive noises and propose a recursive least squares algorithm using the data filtering technique. For the purpose of comparison, an auxiliary model-based recursive generalized least squares algorithm is developed. The simulation results show that the filtering-based least squares algorithm can provide more highly accurate parameter estimation under the colored noises.

The rest of the article is organized as follows. Section “Problem formulation” introduces a class of linear-in-parameter MISO autoregressive systems and derives the identification model. Section “Recursive generalized least squares algorithm” proposes a recursive generalized least squares algorithm using the auxiliary model idea. Section “Filtering-based recursive least squares algorithm” presents a recursive least squares algorithm using the data filtering technique. Section “Example” provides an illustrative example to show the effectiveness of the algorithms. Finally, concluding remarks are given in section “Conclusion.”

Problem formulation

Let us define some symbols. The symbol ${\mathit{I}}_n$ denotes an identity matrix of order n; $1_n$ denotes an n-dimensional column vector whose elements are 1; the superscript T denotes the matrix/vector transpose.

Consider the following MISO stochastic system shown in Figure 1

$$y(t)=\frac{{\mathit{\vartheta }}_1^T{\mathit{\eta }}_1(u_1(t))}{A_1(z)}+\frac{{\mathit{\vartheta }}_2^T{\mathit{\eta }}_2(u_2(t))}{A_2(z)}+\cdots +\frac{{\mathit{\vartheta }}_r^T{\mathit{\eta }}_{r}x(u_r(t))}{A_r(z)}+\frac{1}{C(z)}v(t)$$

(1)

where $y(t)\in \mathbb{R}$ is the system output, $u_j(t)\in \mathbb{R},j=1,2,\dots ,r$ are the system inputs, $v(t)\in \mathbb{R}$ is the stochastic white noise with zero mean. $A_j(z)$ and $C(z)$ are polynomials, of known orders $(n_j,n_c)$ , in the unit backward shift operator $z^{-1}$ , and defined by

$$\begin{array}{c}{A}_j(z):=1+a_{j1}{z}^{-1}+a_{j2}{z}^{-2}+\cdots +a_{j{n}_j}{z}^{-n_j}\in \mathbb{R}\hfill \\ C(z):=1+c_1{z}^{-1}+c_2{z}^{-2}+\cdots +c_{n_c}{z}^{-n_c}\in \mathbb{R}\hfill \end{array}$$

where $a_{\mathrm{ji}}$ , $c_i$ , and ${\mathit{\vartheta }}_j^T$ are the unknown parameters to be estimated. Note that system (1) represents a class of linear-in-parameter systems. When ${\eta }_j(u_j(t))$ is a linear function of $u_j(t)$ and $C(z)=D(z)=1$ , system (1) reduces to a linear MISO output-error system. ²⁹ In this work, we assume that ${\eta }_j(u_j(t))$ is a nonlinear function of $u_j(t)$ , for example

$$\begin{array}{c}{\mathit{\eta }}_j(u_j(t)):=[u_j^2(t),\sin (u_j(t-1)),\dots ,\hfill \\ \sqrt{2}\cos (u_j(t-n_j)){]}^T\in {\mathbb{R}}^{n_j}\hfill \end{array}$$

and system (1) denotes a nonlinear system.

Figure 1.

Linear-in-parameter multiple-input single-output output-error autoregressive systems.

In the next paragraph, we develop new identification algorithms for estimating the parameter vector ${\mathit{\vartheta }}_j$ and the parameters of $A_j(z)$ and $C(z)$ for system (1) by utilizing the input-output measured data $\{u_j(t),y(t):t=1,2,\dots \}$ . Without loss of generality, assume that $u_j(t)=0$ , $y(t)=0$ , and $v(t)=0$ for $t\le 0$ .

Define the intermediate variables

$$\begin{array}{c}\multicolumn{1}{c}{x_j(t):=\frac{{\mathit{\vartheta }}_j^T{\mathit{\eta }}_j(u_j(t))}{A_j(z)}}\\ \multicolumn{1}{c}{=-\sum _{i=1}^{n_a}{a}_{\mathrm{ji}}{x}_j(t-i)+{\mathit{\vartheta }}_j^T{\mathit{\eta }}_j(u_j(t))}\end{array}$$

(2)

$$\begin{array}{c}\multicolumn{1}{c}{w(t):=\frac{1}{C(z)}v(t)}\\ \multicolumn{1}{c}{=-\sum _{i=1}^{n_c}{c}_{i}w(t-i)+v(t)}\end{array}$$

(3)

Define the parameter vectors

$$\mathit{\theta }:=\left[\begin{array}{c}{\mathit{\theta }}_{\mathrm{s}}\hfill \\ \multicolumn{1}{c}{{\mathit{\theta }}_{\mathrm{n}}}\end{array}\right]\in {\mathbb{R}}^{n_0+n_c}$$

$${\mathit{\theta }}_{\mathrm{s}}:=[{\mathit{\theta }}_1,{\mathit{\theta }}_2,\dots ,{\mathit{\theta }}_r{]}^T\in {\mathbb{R}}^{n_0}$$

$${\mathit{\theta }}_{\mathrm{n}}:=[c_1,c_2,\dots ,c_{n_c}{]}^T\in {\mathbb{R}}^{n_c}$$

$${\mathit{\theta }}_j:=[a_{j1},a_{j2},\dots ,a_{j{n}_j},{\mathit{\vartheta }}_j^T{]}^T\in {\mathbb{R}}^{2n_j}$$

$$n_0:=2(n_1+n_2,\dots ,+n_r)$$

and the information vectors

$$\mathit{\phi }(t):=\left[\begin{array}{c}{\mathit{\phi }}_{\mathrm{s}}(t)\hfill \\ \multicolumn{1}{c}{{\mathit{\phi }}_{\mathrm{n}}(t)}\end{array}\right]\in {\mathbb{R}}^{n_0+n_c}$$

$${\mathit{\phi }}_{\mathrm{s}}(t):=[{\mathit{\phi }}_1^T(t),{\mathit{\phi }}_2^T(t),\dots ,{\mathit{\phi }}_r^T(t){]}^T\in {\mathbb{R}}^{n_0}$$

$${\mathit{\phi }}_{\mathrm{n}}(t):=[-w(t-1),\dots ,-w(t-n_c){]}^T\in {\mathbb{R}}^{n_c}$$

$${\mathit{\phi }}_j(t):=[-x_j(t-1),\dots ,-x_j(t-n_j),{\mathit{\eta }}_j^T(u_j(t)){]}^T\in {\mathbb{R}}^{2n_j}$$

Then, equations (2) and (3) can be expressed as

$$x_j(t)={\mathit{\phi }}_j^T(t){\mathit{\theta }}_j$$

(4)

$$w(t)={\mathit{\phi }}_{\mathrm{n}}^T(t){\mathit{\theta }}_{\mathrm{n}}+v(t)$$

(5)

and system (1) can be rewritten as

$$\begin{array}{c}\multicolumn{1}{c}{y(t)=\sum _{j=1}^r{x}_j(t)+w(t)}\\ \multicolumn{1}{c}{={\mathit{\phi }}_{\mathrm{s}}^T(t){\mathit{\theta }}_{\mathrm{s}}+{\mathit{\phi }}_{\mathrm{n}}^T(t){\mathit{\theta }}_{\mathrm{n}}+v(t)}\\ \multicolumn{1}{c}{={\mathit{\phi }}^T(t)\mathit{\theta }+v(t)}\end{array}$$

(6)

Equation (6) is the identification model of system (1), and parameter vector $\theta$ consists of all the parameters of the system. In order to further develop the filter-based algorithms, in next section, we first introduce the recursive generalized least squares algorithms using the auxiliary model identification idea.

Recursive generalized least squares algorithm

Define a quadratic criterion function

$$J(\mathit{\theta }):=\sum _{i=1}^t[y(i)-{\mathit{\phi }}^T(i)\mathit{\theta }{]}^2$$

Minimizing $J(\mathit{\theta })$ and letting the partial derivative of $J(\mathit{\theta })$ with regard to $\mathit{\theta }$ be zero, we obtain the least squares estimate of $\mathit{\theta }$

$$\hat{\mathit{\theta }}(t)={[\sum _{i=1}^t{\mathit{\phi }}^T(i)\mathit{\phi }(i)]}^{-1}[\sum _{i=1}^t{\mathit{\phi }}^T(i)y(i)]$$

(7)

In order to recursively compute $\hat{\mathit{\theta }}(t)$ , we define the covariance matrix $\mathit{P}(t)$ and the vector $\mathit{\xi }(t)$ as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}^{-1}(t):=\sum _{i=1}^t{\mathit{\phi }}^T(i)\phi (i)}\\ \multicolumn{1}{c}{={\mathit{P}}^{-1}(t-1)+{\mathit{\phi }}^T(t)\mathit{\phi }(t)}\end{array}$$

(8)

$$\mathit{\xi }(t):=\sum _{i=1}^t\mathit{\phi }(i)y(i)=\xi (t-1)+\mathit{\phi }(t)y(t)$$

(9)

$$\mathit{P}(0)=p_0{\mathit{I}}_n>0,\xi (0)=0$$

Then, equation (7) can be expressed as

$$\hat{\mathit{\theta }}(t)=\mathit{P}(t)\xi (t)$$

(10)

Applying the matrix inverse lemma

$$(A+\mathrm{BC}{)}^{-1}=A^{-1}+A^{-1}B(I+{\mathrm{CA}}^{-1}B{)}^{-1}{\mathrm{CA}}^{-1}$$

to equation (8) gives

$$\mathit{P}(t)=\mathit{P}(t-1)-\frac{\mathit{P}(t-1)\mathit{\phi }(t){\mathit{\phi }}^T(t)\mathit{P}(t-1)}{1+{\mathit{\phi }}^T(t)\mathit{P}(t-1)\mathit{\phi }(t)}$$

(11)

The recursive algorithm (9)–(11) cannot be implemented because the intermediate variables $x_j(t-i)$ and $w(t-i)$ in $\mathit{\phi }(t)$ are unknown. To deal with the difficulties, we adopt the auxiliary model identification idea. ³¹ Let ${\hat{x}}_j(t)$ and $\hat{w}(t)$ be the estimates of $x_j(t)$ and $w(t)$ at time t, substituting ${\hat{x}}_j(t-i)$ into ${\mathit{\phi }}_j(t)$ and $\hat{w}(t-i)$ into ${\mathit{\phi }}_n(t)$ , we can construct the estimated vectors ${\hat{\mathit{\phi }}}_j(t)$ and ${\hat{\mathit{\phi }}}_{\mathrm{n}}(t)$ , and form $\hat{\mathit{\phi }}(t)$ . Replacing ${\mathit{\phi }}_j(t)$ in equation (4) with its estimate ${\hat{\mathit{\phi }}}_j(t)$ , we can compute ${\hat{x}}_j(t)$ through

$${\hat{x}}_j(t)={\hat{\mathit{\phi }}}_j^T(t){\hat{\mathit{\theta }}}_j(t)$$

From equation (6), we have

$$w(t)=y(t)-\sum _{j=1}^r{x}_j(t)$$

Replacing $x_j(t)$ on the right-hand side of the above equation with their estimates ${\hat{x}}_j(t)$ , we can compute $\hat{w}(t)$ through

$$\hat{w}(t)=y(t)-\sum _{j=1}^r{\hat{x}}_j(t)$$

Replacing $\mathit{\phi }(t)$ in equations (9)–(11) with its estimate $\hat{\mathit{\phi }}(t)$ , we can summarize the following auxiliary model-based recursive generalized least squares algorithm for identifying $\mathit{\theta }$ (the MISO-AM-RGLS algorithm for short)

$$\hat{\mathit{\theta }}(t)=\mathit{P}(t)\mathit{\xi }(t)$$

(12)

$$\mathit{\xi }(t)=\mathit{\xi }(t-1)+\hat{\mathit{\phi }}(t)y(t)$$

(13)

$$\mathit{P}(t)=\mathit{P}(t-1)-\frac{\mathit{P}(t-1)\hat{\mathit{\phi }}(t){\hat{\mathit{\phi }}}^T(t)\mathit{P}(t-1)}{1+{\hat{\mathit{\phi }}}^T(t)\mathit{P}(t-1)\hat{\mathit{\phi }}(t)}$$

(14)

$$\begin{array}{c}\multicolumn{1}{c}{\hat{\mathit{\phi }}(t)=[{\hat{\mathit{\phi }}}_1^T(t),{\hat{\mathit{\phi }}}_2^T(t),\dots ,{\hat{\mathit{\phi }}}_r^T(t),-\hat{w}(t-1),}\\ \multicolumn{1}{c}{-\hat{w}(t-2),\dots ,-\hat{w}(t-n_c){]}^T}\end{array}$$

(15)

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{\phi }}_j(t)=[-{\hat{x}}_j(t-1),-{\hat{x}}_j(t-2),\dots ,}\\ \multicolumn{1}{c}{-{\hat{x}}_j(t-n_j),{\mathit{\eta }}_j^T(u_j(t)){]}^T}\end{array}$$

(16)

$${\hat{x}}_j(t)={\hat{\mathit{\phi }}}_j^T(t){\hat{\mathit{\theta }}}_j(t)$$

(17)

$$\hat{w}(t)=y(t)-\sum _{j=1}^r{\hat{x}}_j(t)$$

(18)

$$\hat{\mathit{\theta }}(t)=[{\hat{\mathit{\theta }}}_1^T(t),{\hat{c}}_1(t),{\hat{c}}_2(t),\dots ,{\hat{c}}_{n_c}(t){]}^T$$

(19)

$${\hat{\mathit{\theta }}}_j(t)=[{\hat{a}}_1(t),{\hat{a}}_2(t),\dots ,{\hat{a}}_{n_a}(t),{\hat{\mathit{\vartheta }}}^T(t){]}^T$$

(20)

To initialize the algorithm (12)–(20), we set $\hat{\mathit{\theta }}(i)=1_{n_0+n_c}/p_0$ , ${\hat{x}}_j(i)=1/p_0$ , and $\hat{w}(i)=1/p_0$ for $\mathrm{ileqslant}0$ , and $\mathit{P}(t)=p_0{\mathit{I}}_{n_0}+n_c$ , $p_0={10}^6$ .

Filtering-based recursive least squares algorithm

Define the filtered variables

$$y_{\mathrm{f}}(t):=C(z)y(t)$$

$$x_{j,\mathrm{f}}(t):=C(z)x_j(t)$$

and the filtered information vectors

$${\mathit{\phi }}_{\mathrm{f}}(t):=[{\mathit{\phi }}_{1,\mathrm{f}}^T(t),{\mathit{\phi }}_{2,\mathrm{f}}^T(t),\dots ,{\mathit{\phi }}_{r,\mathrm{f}}^T(t){]}^T\in {\mathbb{R}}^{n_0}$$

$${\mathit{\phi }}_{j,\mathrm{f}}(t):=[-x_{j,\mathrm{f}}(t-1),\dots ,-x_{j,\mathrm{f}}(t-n_a),{\mathit{\eta }}_{j,\mathrm{f}}^T(t){]}^T\in {\mathbb{R}}^{2n_j}$$

$${\mathit{\eta }}_{j,\mathrm{f}}(t):=C(z)\mathit{\eta }(u_j(t))\in {\mathbb{R}}^{n_j}$$

System (1) can be rewritten as

$$y(t)=\sum _{j=1}^r{x}_j(t)+\frac{1}{C(z)}v(t)$$

Multiplying both sides of the above equation by $C(z)$ gives

$$y_{\mathrm{f}}(t)=\sum _{j=1}^r{x}_{j,\mathrm{f}}(t)+v(t)={\mathit{\phi }}_{\mathrm{f}}^T(t){\mathit{\theta }}_{\mathrm{s}}+v(t)$$

(21)

Equations (4) and (21) are the noise identification model and the filtered identification model of system (1), respectively. Minimizing two quadratic cost functions

$$J_{\mathrm{f}}({\mathit{\theta }}_{\mathrm{s}}):=\sum _{i=1}^t[y_{\mathrm{f}}(i)-{\mathit{\phi }}_{\mathrm{f}}^T(i){\mathit{\theta }}_{\mathrm{s}}{]}^2$$

$$J_{\mathrm{n}}({\mathit{\theta }}_{\mathrm{n}}):=\sum _{i=1}^t[w(i)-{\mathit{\phi }}_{\mathrm{n}}^T(i){\theta }_{\mathrm{n}}{]}^2$$

leads to the following recursive parameter updates

$${\hat{\mathit{\theta }}}_{\mathrm{s}}(t)={\mathit{P}}_{\mathrm{f}}(t){\mathit{\xi }}_{\mathrm{f}}(t)$$

(22)

$${\mathit{\xi }}_{\mathrm{f}}(t)={\mathit{\xi }}_{\mathrm{f}}(t-1)+{\mathit{\phi }}_{\mathrm{f}}(t)y_{\mathrm{f}}(t)$$

(23)

$${\mathit{P}}_{\mathrm{f}}(t)={\mathit{P}}_{\mathrm{f}}(t-1)-\frac{{\mathit{P}}_{\mathrm{f}}(t-1){\phi }_{\mathrm{f}}(t){\mathit{\phi }}_{\mathrm{f}}^T(t){\mathit{P}}_{\mathrm{f}}(t-1)}{1+{\mathit{\phi }}_{\mathrm{f}}^T(t){\mathit{P}}_{\mathrm{f}}(t-1){\mathit{\phi }}_{\mathrm{f}}(t)}$$

(24)

$${\hat{\mathit{\theta }}}_{\mathrm{n}}(t)={\mathit{P}}_{\mathrm{n}}(t){\mathit{\xi }}_{\mathrm{n}}(t)$$

(25)

$${\mathit{\xi }}_{\mathrm{n}}(t)={\mathit{\xi }}_{\mathrm{n}}(t-1)+{\mathit{\phi }}_{\mathrm{n}}(t)w(t)$$

(26)

$${\mathit{P}}_{\mathrm{n}}(t)={\mathit{P}}_{\mathrm{n}}(t-1)-\frac{{\mathit{P}}_{\mathrm{n}}(t-1){\mathit{\phi }}_{\mathrm{n}}(t){\mathit{\phi }}_{\mathrm{n}}^T(t){\mathit{P}}_{\mathrm{n}}(t-1)}{1+{\mathit{\phi }}_{\mathrm{n}}^T(t){\mathit{P}}_{\mathrm{n}}(t-1){\mathit{\phi }}_{\mathrm{n}}(t)}$$

(27)

Equations (22)–(27) cannot be implemented because the filtered vector ${\mathit{\phi }}_{\mathrm{f}}(t)$ is unknown due to the unknown filter $F(z)$ and the noise items $w(t-i)$ in ${\mathit{\phi }}_{\mathrm{n}}(t)$ is unmeasurable. Here, we employ the interactive estimation theory and the idea of replacing the unknown items with their corresponding estimates. Let ${\hat{\mathit{\theta }}}_{\mathrm{s}}(t)$ , ${\hat{\mathit{\theta }}}_{\mathrm{n}}(t)$ , and $\hat{\mathit{\theta }}(t)$ be the estimates of the parameters ${\mathit{\theta }}_{\mathrm{s}}$ , ${\mathit{\theta }}_{\mathrm{n}}$ , and $\mathit{\theta }$ , and $\hat{x}(t-i)$ and $\hat{w}(t-i)$ be the estimates of $x(t-i)$ and $w(t-i)$ . Define the estimated information vectors

$${\hat{\mathit{\phi }}}_{\mathrm{s}}(t):=[{\hat{\mathit{\phi }}}_1^T(t),{\hat{\mathit{\phi }}}_2^T(t),\dots ,{\hat{\mathit{\phi }}}_r^T(t){]}^T\in {\mathbb{R}}^{n_0}$$

$${\hat{\mathit{\phi }}}_{j,\mathrm{f}}(t):=[-{\hat{x}}_{j,\mathrm{f}}(t-1),\dots ,-{\hat{x}}_{j,\mathrm{f}}(t-n_j),{\hat{\mathit{\eta }}}_{j,\mathrm{f}}^T(t){]}^T$$

$${\hat{\mathit{\phi }}}_{\mathrm{n}}(t):=[-\hat{w}(t-1),-\hat{w}(t-2),\dots ,-\hat{w}(t-n_c){]}^T$$

Replacing ${\mathit{\phi }}_j(t)$ in equation (4) and ${\mathit{\phi }}_{\mathrm{n}}(t)$ in equation (5) with ${\hat{\mathit{\phi }}}_j(t)$ and ${\hat{\mathit{\phi }}}_{\mathrm{n}}(t)$ yields

$${\hat{x}}_j(t)={\hat{\mathit{\phi }}}_j^T(t){\hat{\mathit{\theta }}}_j(t)$$

(28)

$$\hat{w}(t)=y(t)-{\hat{\mathit{\phi }}}_{\mathrm{s}}^T(t){\hat{\theta }}_{\mathrm{s}}(t)$$

(29)

Using the parameter estimate of the noise model

$${\hat{\mathit{\theta }}}_{\mathrm{n}}:=[{\hat{c}}_1(t),{\hat{c}}_2(t),\dots ,{\hat{c}}_{n_c}(t){]}^T\in {\mathbb{R}}^{n_c}$$

to construct the estimate of $C(z)$ as

$$\hat{C}(z)=1+{\hat{c}}_1(t)z^{-1}+{\hat{c}}_2(t)z^{-2}+\cdots +{\hat{c}}_{n_c}(t)z^{-n_c}$$

Then, we compute the estimates of the filtered variables $y_{\mathrm{f}}(t)$ and $x_{j,\mathrm{f}}(t)$ through

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{y}}_{\mathrm{f}}(t):=\hat{C}(z)y(t)=\hat{y}(t)+{\hat{c}}_1(t)\hat{y}(t-1)+{\hat{c}}_2(t)\hat{y}(t-2)}\\ \multicolumn{1}{c}{+\cdots +{\hat{c}}_{n_c}(t)\hat{y}(t-n_c)}\end{array}$$

(30)

$$\begin{array}{c}\multicolumn{1}{c}{{\hat{x}}_{j,\mathrm{f}}(t):=\hat{C}(z)x_j(t)={\hat{x}}_j(t)+{\hat{c}}_1(t){\hat{x}}_j(t-1)+{\hat{c}}_2(t){\hat{x}}_j(t-2)}\\ \multicolumn{1}{c}{+\cdots +{\hat{c}}_{n_c}(t){\hat{x}}_j(t-n_c)}\end{array}$$

(31)

and form the estimates of the filtered information vectors

$${\hat{\mathit{\phi }}}_{\mathrm{f}}(t):=[{\hat{\mathit{\phi }}}_{1,\mathrm{f}}^T(t),{\hat{\mathit{\phi }}}_{2,\mathrm{f}}^T(t),\dots ,{\hat{\mathit{\phi }}}_{n_j,\mathrm{f}}^T(t){]}^T$$

$${\hat{\mathit{\phi }}}_{j,\mathrm{f}}(t):=[-{\hat{x}}_{j,\mathrm{f}}(t-1),\dots ,-{\hat{x}}_{j,\mathrm{f}}(t-n_j),{\hat{\mathit{\eta }}}_{j,\mathrm{f}}^T(t){]}^T$$

$${\hat{\mathit{\eta }}}_{j,\mathrm{f}}(t):=\hat{C}(z){\mathit{\eta }}_j(u_j(t))$$

Replacing ${\mathit{\phi }}_{\mathrm{f}}(t)$ in equations (23) and (24), ${\mathit{\phi }}_{\mathrm{n}}(t)$ in equations (26) and (27), $y_{\mathrm{f}}(t)$ in equation (23), and $w(t)$ in equation (26) with their corresponding estimates ${\hat{\mathit{\phi }}}_{\mathrm{f}}(t)$ , ${\hat{\mathit{\phi }}}_{\mathrm{n}}(t)$ , ${\hat{y}}_{\mathrm{f}}(t)$ , and $\hat{w}(t)$ leads to the following filtering-based least squares algorithms for identifying ${\hat{\mathit{\theta }}}_{\mathrm{s}}(t)$ and ${\hat{\mathit{\theta }}}_{\mathrm{n}}(t)$ (the MISO-F-RLS algorithm for short)

$${\hat{\mathit{\theta }}}_{\mathrm{s}}(t)={\mathit{P}}_{\mathrm{f}}(t){\mathit{\xi }}_{\mathrm{f}}(t)$$

(32)

$${\mathit{\xi }}_{\mathrm{f}}(t)={\mathit{\xi }}_{\mathrm{f}}(t-1)+{\hat{\mathit{\phi }}}_{\mathrm{f}}(t){\hat{y}}_{\mathrm{f}}(t)$$

(33)

$${\mathit{P}}_{\mathrm{f}}(t)={\mathit{P}}_{\mathrm{f}}(t-1)-\frac{{\mathit{P}}_{\mathrm{f}}(t-1){\hat{\mathit{\phi }}}_{\mathrm{f}}(t){\hat{\mathit{\phi }}}_{\mathrm{f}}^T(t){\mathit{P}}_{\mathrm{f}}(t-1)}{1+{\hat{\mathit{\phi }}}_{\mathrm{f}}^T(t){\mathit{P}}_{\mathrm{f}}(t-1){\hat{\mathit{\phi }}}_{\mathrm{f}}(t)}$$

(34)

$${\hat{\mathit{\theta }}}_{\mathrm{n}}(t)={\mathit{P}}_{\mathrm{n}}(t){\xi }_{\mathrm{n}}(t)$$

(35)

$${\mathit{\xi }}_{\mathrm{n}}(t)={\mathit{\xi }}_{\mathrm{n}}(t-1)+{\hat{\mathit{\phi }}}_{\mathrm{n}}(t)\hat{w}(t)$$

(36)

$${\mathit{P}}_{\mathrm{n}}(t)={\mathit{P}}_{\mathrm{n}}(t-1)-\frac{{\mathit{P}}_{\mathrm{n}}(t-1){\hat{\mathit{\phi }}}_{\mathrm{n}}(t){\hat{\mathit{\phi }}}_{\mathrm{n}}^T(t)P_{\mathrm{n}}(t-1)}{1+{\hat{\mathit{\phi }}}_{\mathrm{n}}^T(t){\mathit{P}}_{\mathrm{n}}(t-1){\hat{\mathit{\phi }}}_{\mathrm{n}}(t)]}$$

(37)

$${\hat{\mathit{\phi }}}_{\mathrm{f}}(t)=[{\hat{\mathit{\phi }}}_{1,\mathrm{f}}^T(t),{\hat{\mathit{\phi }}}_{2,\mathrm{f}}^T(t),\dots ,{\hat{\mathit{\phi }}}_{n_j,\mathrm{f}}^T(t){]}^T$$

(38)

$${\hat{\mathit{\phi }}}_{j,\mathrm{f}}(t)=[-{\hat{x}}_{j,\mathrm{f}}(t-1),\dots ,-{\hat{x}}_{j,\mathrm{f}}(t-n_j),{\hat{\mathit{\eta }}}_{j,\mathrm{f}}^T(t){]}^T$$

(39)

$${\hat{\mathit{\phi }}}_j(t)=[-{\hat{x}}_j(t-1),\dots ,-{\hat{x}}_j(t-n_j),{\hat{\mathit{\eta }}}_j^T(t){]}^T$$

(40)

$${\hat{\mathit{\phi }}}_{\mathrm{n}}(t)=[-\hat{w}(t-1),\dots ,-\hat{w}(t-n_c){]}^T$$

(41)

$${\hat{\mathit{\eta }}}_{j,\mathrm{f}}(t)=\hat{C}(z){\mathit{\eta }}_j(u_j(t))$$

(42)

$${\hat{y}}_{\mathrm{f}}(t)=\hat{y}(t)+{\hat{c}}_1(t)\hat{y}(t-1)+{\hat{c}}_2(t)\hat{y}(t-2)+\cdots +{\hat{c}}_{n_c}(t)\hat{y}(t-n_c)$$

(43)

$${\hat{x}}_{j,\mathrm{f}}(t)={\hat{x}}_j(t)+{\hat{c}}_1(t){\hat{x}}_j(t-1)+{\hat{c}}_2(t){\hat{x}}_j(t-2)+\cdots +{\hat{c}}_{n_c}(t){\hat{x}}_j(t-n_c)$$

(44)

$${\hat{x}}_j(t)={\hat{\mathit{\phi }}}_j^T(t){\hat{\mathit{\theta }}}_{\mathrm{s}}(t)$$

(45)

$$\hat{w}(t)=y(t)-{\hat{\mathit{\phi }}}_{\mathrm{s}}^T(t){\hat{\mathit{\theta }}}_{\mathrm{s}}(t)$$

(46)

The steps of computing ${\hat{\mathit{\theta }}}_{\mathrm{s}}(t)$ and ${\hat{\mathit{\theta }}}_{\mathrm{n}}(t)$ involved in the algorithm are summarized as follows:

1. Let $t=1$ , set ${\hat{\mathit{\theta }}}_{\mathrm{s}}(i)=1_{n_0}/p_0$ and ${\hat{\mathit{\theta }}}_{\mathrm{n}}(i)=1_{n_c}/p_0$ ; ${\hat{x}}_{j,\mathrm{f}}(i)=1/p_0$ , ${\hat{y}}_{\mathrm{f}}(i)=1/p_0$ , ${\hat{x}}_j(i)=1/p_0$ , and $\hat{w}(i)=1/p_0$ for $\mathrm{ileqslant}0$ , ${\mathit{P}}_{\mathrm{f}}(t)=p_0{\mathit{I}}_{n_0}$ and ${\mathit{P}}_{\mathrm{n}}(t)=p_0{\mathit{I}}_{n_c}$ , $p_0={10}^6$ .

2. Collect the input-output data $u_j(t)$ and $y(t)$ .

3. Compute the intermediate variable ${\mathit{\xi }}_{\mathrm{n}}(t)$ using equation (36) and the covariance matrix ${\mathit{P}}_{\mathrm{n}}(t)$ using equation (37), and the filtered vector ${\hat{\mathit{\eta }}}_{j,\mathrm{f}}(t)$ using equation (42).

4. Update the parameter estimate ${\hat{\mathit{\theta }}}_{\mathrm{n}}(t)$ using equation (35).

5. Compute filtered variables ${\hat{y}}_{\mathrm{f}}(t)$ using equation (43) and ${\hat{x}}_{j,\mathrm{f}}(t)$ using equation (44). Construct information vectors ${\hat{\mathit{\phi }}}_{j,\mathrm{f}}(t)$ using equation (39) and form ${\hat{\mathit{\phi }}}_{\mathrm{f}}(t)$ using equation (38).

6. Compute the intermediate variable ${\mathit{\xi }}_{\mathrm{f}}(t)$ using equation (41) and the covariance matrix ${\mathit{P}}_{\mathrm{f}}(t)$ using equation (34).

7. Update the parameter estimate ${\hat{\mathit{\theta }}}_{\mathrm{s}}(t)$ using equation (32).

8. Compute ${\hat{x}}_j(t)$ using equation (45) and $\hat{w}(t)$ using equation (46).

9. Increase t by 1 and go to step 2.

Example

An example is given to demonstrate the effectiveness of the proposed algorithms. Consider the following nonlinear MISO autoregressive system

$$y(t)=\frac{{\mathit{\vartheta }}_1^T{\mathit{\eta }}_1(u_1(t))}{A_1(z)}+\frac{{\mathit{\vartheta }}_2^T{\mathit{\eta }}_2(u_2(t))}{A_2(z)}+\frac{1}{C(z)}v(t)$$

$$A_1(z)=1+a_{11}{z}^{-1}=1-0.37z^{-1}$$

$$A_2(z)=1+a_{21}{z}^{-1}=1-0.36z^{-1}$$

$$C(z)=1+c_1{z}^{-1}=1+0.43z^{-1}$$

$${\mathit{\eta }}_1(u_1(t))=[u_1(t-1){]}^T,{\mathit{\eta }}_2(u_2(t))=[0.5u_2^2(t-1){]}^T$$

$${\mathit{\vartheta }}_1=[b_{11}]=[1.92{]}^T,{\mathit{\vartheta }}_2=[b_{21}]=[1.71{]}^T$$

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{\theta }=[a_{11},b_{11},a_{21},b_{21},c_1{]}^T}\\ \multicolumn{1}{c}{=[-0.37,1.92,-0.36,1.71,0.43{]}^T}\end{array}$$

Here, the inputs $\{u_1(t)\}$ and $\{u_2(t)\}$ are taken as uncorrelated persistent excitation signal sequences with zero means and unit variances ${\sigma }^2=1.{00}^2$ , and $\{v(t)\}$ as a white noise sequence with zero mean.

Using $L=3000$ data and applying the MISO-AM-RGLS algorithm in equations (12)–(20) and the MISO-F-RLS algorithm in equations (32)–(46) to estimate the parameters of this system, we have the parameter estimates of each algorithm and their errors shown in Table 1. The parameter estimation errors $\mathit{\delta }:=\parallel \hat{\mathit{\theta }}(t)-\theta \parallel /\parallel \mathit{\theta }\parallel$ of each algorithm are illustrated in Figure 2. We also investigate the performance of two algorithms under a relatively high noise level with noise variance ${\sigma }^2=2.{00}^2$ , and the corresponding simulation results are shown in Table 2 and Figure 3.

Table 1.

Parameter estimates and errors of two algorithms ( ${\sigma }^2=1.{00}^2$ ).

Algorithm	t	$a_1$	$b_1$	$a_2$	$b_2$	$c_1$	$\delta (\%)$
MISO-AM-RGLS	100	−0.40225	1.88121	−0.39152	1.51705	0.45924	7.67723
	200	−0.43206	1.86657	−0.42636	1.50374	0.44105	8.72618
	500	−0.38271	1.88901	−0.45826	1.50328	0.48582	8.95445
	1000	−0.37510	1.88141	−0.42254	1.54181	0.43925	6.91840
	2000	−0.38450	1.90960	−0.39575	1.57727	0.45323	5.28890
	3000	−0.37168	1.91495	−0.38149	1.63437	0.43828	2.98169
MISO-F-RLS	100	−0.43547	1.81912	−0.40881	1.49686	0.36082	9.74465
	200	−0.45862	1.82965	−0.43603	1.54089	0.38822	8.59276
	500	−0.40191	1.85752	−0.46840	1.58530	0.45769	6.83482
	1000	−0.38905	1.86563	−0.43613	1.62300	0.41933	4.87697
	2000	−0.39383	1.89119	−0.40581	1.65900	0.43787	2.95343
	3000	−0.38217	1.90506	−0.39399	1.71096	0.42758	1.47377
True values		−0.37000	1.92000	−0.36000	1.71000	0.43000

Figure 2.

Estimation errors $\delta$ versus t ( ${\sigma }^2=1.{00}^2$ ).

Table 2.

Parameter estimates and errors of two algorithms ( ${\sigma }^2=2.{00}^2$ ).

Algorithm	t	$a_1$	$b_1$	$a_2$	$b_2$	$c_1$	$\delta (\%)$
MISO-AM-RGLS	100	−0.40998	1.83268	−0.40413	1.37836	0.45969	13.14586
	200	−0.48425	1.80825	−0.47401	1.34306	0.43978	15.66428
	500	−0.39221	1.85506	−0.54720	1.31548	0.48580	16.76624
	1000	−0.37927	1.84093	−0.47480	1.39148	0.43886	13.09249
	2000	−0.39923	1.89834	−0.42556	1.45855	0.45334	9.91290
	3000	−0.37302	1.90946	−0.39826	1.56997	0.43839	5.48722
MISO-F-RLS	100	−0.46335	1.80848	−0.43694	1.32093	0.43259	15.89646
	200	−0.51580	1.78277	−0.48750	1.40339	0.42127	14.59524
	500	−0.41998	1.81901	−0.52449	1.49352	0.47407	11.19835
	1000	−0.39289	1.82251	−0.46778	1.54505	0.43058	8.31758
	2000	−0.40320	1.87261	−0.41948	1.58831	0.44654	5.57759
	3000	−0.37874	1.89413	−0.39523	1.68482	0.43475	1.93488
True values		−0.37000	1.92000	−0.36000	1.71000	0.43000

Figure 3.

Estimation errors $\delta$ versus t ( ${\sigma }^2=2.{00}^2$ ).

From Tables 1 and 2 and Figures 2 and 3, we can draw the following conclusion:

• The parameter estimation errors $\delta$ are becoming smaller as t increases.

• Both algorithms can provide high-accuracy parameter estimation, and the filtered algorithm achieves better estimation accuracy over the non-filtered algorithm under different noise variances.

Conclusion

In this article, we have presented two recursive least squares algorithms, an auxiliary model-based generalized least squares algorithm, and a filtered recursive least squares algorithm for a class of nonlinear MISO output-error autoregressive systems. The illustrative example shows that the filtering-based least squares algorithm can produce more highly accurate estimates over the recursive generalized least squares algorithm. The proposed methods can be extended to nonlinear systems with colored noise^32–35 and applied to other fields.^36–38