Partially Coupled Stochastic Gradient Estimation for Multivariate Equation-Error Systems

Abstract

This paper researches the identification problem for the unknown parameters of the multivariate equation-error autoregressive systems. Firstly, the original identification model is decomposed into several sub-identification models according to the number of system outputs. Then, based on the characteristic that the information vector and the parameter vector are common among the sub-identification models, the coupling identification concept is used to propose a partially coupled generalized stochastic gradient algorithm. Furthermore, by expanding the scalar innovation of each subsystem model to the innovation vector, a partially coupled multi-innovation generalized stochastic gradient algorithm is proposed. Finally, the numerical simulations indicate that the proposed algorithms are effective and have good parameter estimation performances.

1. Introduction

Parameter estimation is an important part in the field of system identification, which usually identifies the unknown parameters in the models according to the input and output data of the system when the model structure of a system is known [1,2,3]. Parameter estimation has been used in many fields in recent years, including chemistry, mechanics, engineering and so on [4,5,6]. For example, in chemical engineering field, Khalik et al. applied the parameter estimation approach by using current/voltage data in achieving physically meaningful parameters of the Doyle–Fuller–Newman model for Lithium-ion batteries [7]. In mechanics field, Shamrao et al. estimated the terramechanics parameters by applying dynamic Bayesian estimation techniques on the measurements from simple single-wheel tests [8]. In engineering, Calasan et al. proposed two algorithms for the transformer parameter estimation to improve the estimation process and prevent inaccuracies and parameter mismatch with the real parameters of the transformer [9].

In many identification objects of parameter estimation, multivariate systems is a very common class [10,11,12]. The majority of industrial processes are the multivariate systems [13,14]. It is not easy to estimate the unknown parameters of multivariate systems, because such systems are quite complex; for example, there are many variables, there is coupling between variables, or there is time delay [15,16,17]. In recent years, the identification of the multivariate systems has attracted the attention of many scholars. Shafin et al. studied the angle and delay estimation for 3D massive MIMO systems under a parametric channel modeling, which was crucial for such systems to realize the predicted capacity gains [18]. Kawaria et al. designed a Levy shuffled frog leaping algorithm with high parameter estimation efficiency to estimate the parameters of multiple-input–multiple-output bilinear systems [19]. Roy et al. developed an online plant-parameter identification method for the multi-input–multi-output linear time-invariant systems [20].

The coupling identification concept is a useful identification method that has been developed in recent decades. It is suitable for parameter estimation problems for the multivariate systems with the same parameters among subsystems [21,22]. The identification algorithm based on the coupling concept has the advantage of lower computation amounts compared with the traditional identification algorithm for multivariate systems. Cui et al. combined the Kalman filtering principle and the coupling identification concept to derive a Kalman filtering-based partially coupled recursive least squares algorithm for jointly estimating the parameters and the states of the multivariable state–space system; their algorithm had high computational efficiency [23]. The multi-innovation identification theory is also an effective method that has been used to improve parameter estimation accuracy in recent years [24,25,26]. This method not only utilizes the data information of the current time, but also makes full use of the data information of the previous time [27,28,29]. Chaudhary et al. presented a multi-innovation fractional least mean square adaptive algorithm for the input-nonlinear systems by expanding the scalar innovation into a vector innovation by using the multi-innovation identification theory [30]. This method can also be applied to the parameter estimation of the multivariate systems.

We have studied the parameter identification problems of multivariate systems in the past. In [31], the colored noise of the original system was filtered into white noise by using the data filtering technique. In [32], the original identification system was decomposed into two sub-identification systems, where one contains the system model parameter vector and the other contains the noise model parameter matrix by applying the decomposition method. However, the coupling identification concept used in this paper is different from the previous two identification methods, as it decomposes the original system according to the number of system outputs to obtain several subsystems with partial parameter vectors and information vectors coupled. It has significant advantages in reducing the computational cost of the algorithm because it makes the original complex system simplified. Therefore, this paper presents new identification methods to estimate the parameters of the multivariate equation-error autoregressive system, which have high computational efficiency and high parameter estimation accuracy. The main contributions of this paper lie in the following.

• This paper decomposes the multivariate equation-error autoregressive system into several sub-identification models according to the number of the system outputs.
• A multivariate partially coupled generalized stochastic gradient (M-PC-GSG) algorithm is proposed for the multivariate equation-error system by utilizing the coupling identification concept, which can reduce the computation amounts compared with the traditional stochastic gradient algorithm.
• A multivariate partially coupled multi-innovation generalized stochastic gradient (M-PC-MI-GSG) algorithm is proposed by using the multi-innovation identification theory, which has higher parameter estimation accuracy than the M-PC-GSG algorithm.

The rest of this paper is organized as follows. Section 2 presents a multivariate equation-error autoregressive system and describes its identification difficulties. Section 3 proposes a partially coupled generalized stochastic gradient algorithm and gives its schematic diagram. Section 4 proposes a partially coupled multi-innovation generalized stochastic gradient algorithm. Section 5 presents two numerical examples to indicate that the proposed algorithms are effective. Finally, we offer some concluding remarks in Section 6.

2. System Description and Identification Model

First of all, we give some notation in this paper. ${\mathbf{I}}_m$ denotes an identity matrix of size $m\times m$; ${\mathbf{1}}_n$ stands for an n-dimensional column vector whose elements are 1, that is ${\mathbf{1}}_n=[1,1,\cdots ,1{]}^{\mathrm{T}}\in {\mathrm{ℝ}}^n$; ${\mathbf{1}}_{m\times n}$ represents a matrix of size $m\times n$ whose elements are 1; the norm of a matrix $\mathit{X}$ is defined by ${\parallel \mathit{X}\parallel }^2:=\mathrm{tr}\left[\mathit{X}{\mathit{X}}^{\mathrm{T}}\right]$, the superscript T stands for the matrix/vector transpose, the symbol ⊗ represents Kronecker product, for example, $\mathit{A}:=\left[a_{ij}\right]\in {\mathrm{ℝ}}^{m\times n}$, $\mathit{B}:=\left[b_{ij}\right]\in {\mathrm{ℝ}}^{p\times q}$, $\mathit{A}\otimes \mathit{B}=\left[a_{ij}\mathit{B}\right]\in {\mathrm{ℝ}}^{\left(mp\right)\times \left(nq\right)}$, in general, $\mathit{A}\otimes \mathit{B}\ne \mathit{B}\otimes \mathit{A}$; $\mathrm{col}\left[\mathit{X}\right]$ is defined as a vector consisting of all columns of matrix $\mathit{X}$ arranged in order, for example, $\mathit{X}:=[{\mathit{x}}_1{\mathit{x}}_2\cdots {\mathit{x}}_n]\in {\mathrm{ℝ}}^{m\times n}$, ${\mathit{x}}_i\in {\mathrm{ℝ}}^m$$(i=1,2,\cdots ,n)$, $\mathrm{col}\left[\mathit{X}\right]:=[{\mathit{x}}_1^{\mathrm{T}}$,${\mathit{x}}_2^{\mathrm{T}}$,⋯,${\mathit{x}}_n^{\mathrm{T}}{]}^{\mathrm{T}}\in {\mathrm{ℝ}}^{mn}$.

Consider the following multivariate equation-error autoregressive system,

$$\mathit{y}\left(t\right)={\mathit{\Phi }}_{\mathrm{s}}\left(t\right)\mathit{\theta }+{\mathit{C}}^{-1}\left(z\right)\mathit{v}\left(t\right),$$

(1)

where $\mathit{y}\left(t\right):={[y_1\left(t\right),y_2\left(t\right),\cdots ,y_m\left(t\right)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^m$ is the system output vector, ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)\in {\mathrm{ℝ}}^{m\times n}$ is the system information matrix consisting of the input-output data, $\mathit{\theta }\in {\mathrm{ℝ}}^n$ is the system parameter vector to be identified, $\mathit{v}\left(t\right):=[v_1\left(t\right)$,$v_2\left(t\right)$,⋯,$v_m{\left(t\right)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^m$ is a white noise process with zero mean, $\mathit{C}\left(z\right)\in {\mathrm{ℝ}}^{m\times m}$ is a polynomial matrix in the unit backward shift operator $[z^{-1}\mathit{y}\left(t\right)=\mathit{y}(t-1)]$:

$$\begin{array}{c}\hfill \mathit{C}\left(z\right):={\mathit{I}}_m+{\mathit{C}}_1{z}^{-1}+{\mathit{C}}_2{z}^{-2}+\cdots +{\mathit{C}}_{n_c}{z}^{-n_c},{\mathit{C}}_i\in {\mathrm{ℝ}}^{m\times m}.\end{array}$$

Define the noise model,

$$\mathit{w}\left(t\right):={\mathit{C}}^{-1}\left(z\right)\mathit{v}\left(t\right)\in {\mathrm{ℝ}}^m.$$

(2)

Assume that the orders m, n and $n_c$ are known, and $\mathit{y}\left(t\right)=\mathit{0}$, ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)=\mathit{0}$, $\mathit{w}\left(t\right)=\mathit{0}$ and $\mathit{v}\left(t\right)=\mathit{0}$ for $t⩽0$.

Define the parameter matrix $\mathit{\gamma }$ and the information vector $\mathit{\varphi }\left(t\right)$ as

$$\begin{array}{ccc}\hfill {\mathit{\gamma }}^{\mathrm{T}}& :=& [{\mathit{C}}_1,{\mathit{C}}_2,\cdots ,{\mathit{C}}_{n_c}]\in {\mathrm{ℝ}}^{m\times \left(m{n}_c\right)},\hfill \\ \hfill \mathit{\varphi }\left(t\right)& :=& {[-{\mathit{w}}^{\mathrm{T}}(t-1),-{\mathit{w}}^{\mathrm{T}}(t-2),\cdots ,-{\mathit{w}}^{\mathrm{T}}(t-n_c)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^{m{n}_c}.\hfill \end{array}$$

From Equation (2), we have

$$\begin{array}{ccc}\hfill \mathit{w}\left(t\right)& =& [{\mathit{I}}_m-\mathit{C}\left(z\right)]\mathit{w}\left(t\right)+\mathit{v}\left(t\right),\hfill \\ & =& (-{\mathit{C}}_1{z}^{-1}-{\mathit{C}}_2{z}^{-2}-\cdots -{\mathit{C}}_{n_c}{z}^{-n_c})\mathit{w}\left(t\right)+\mathit{v}\left(t\right)\hfill \\ & =& -{\mathit{C}}_1\mathit{w}(t-1)-{\mathit{C}}_2\mathit{w}(t-2)-\cdots -{\mathit{C}}_{n_c}\mathit{w}(t-n_c)+\mathit{v}\left(t\right).\hfill \\ & =& {\mathit{\gamma }}^{\mathrm{T}}\mathit{\varphi }\left(t\right)+\mathit{v}\left(t\right).\hfill \end{array}$$

(3)

Then, the multivariate equation-error autoregressive system in (1) can be transformed into the following identification model,

$$\begin{array}{ccc}\hfill \mathit{y}\left(t\right)& =& {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)\mathit{\theta }+\mathit{w}\left(t\right)\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill & \hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}\hspace{0.5em}=& {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)\mathit{\theta }+{\mathit{\gamma }}^{\mathrm{T}}\mathit{\varphi }\left(t\right)+\mathit{v}\left(t\right).\hfill \end{array}$$

(5)

For the identification model in (5), the objective of this paper is to identify the unknown parameters $\mathit{\theta }$ and $\mathit{\gamma }$ by researching some identification methods. Currently, the observable data are $\mathit{y}\left(t\right)$ and ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$. Certainly, the most direct identification method is to integrate the information vector $\mathit{\theta }$ and the information matrix $\mathit{\gamma }$ into a new information matrix $\mathit{\vartheta }$, and the parameter matrix ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$ and the parameter vector $\mathit{\varphi }\left(t\right)$ into a new parameter matrix $\mathit{\Phi }\left(t\right)$, then we can obtain the following identification model:

$$\begin{array}{ccc}\hfill \mathit{y}\left(t\right)& =& \mathit{\Phi }\left(t\right)\mathit{\vartheta }+\mathit{v}\left(t\right),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \mathit{\Phi }\left(t\right)& :=& [{\mathit{\Phi }}_{\mathrm{s}}\left(t\right),{\mathit{I}}_m\otimes {\mathit{\varphi }}^{\mathrm{T}}\left(t\right)]\in {\mathrm{ℝ}}^{m\times n_0},n_0:=n+m^2{n}_c,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \mathit{\vartheta }& :=& \left[\begin{array}{c}\mathit{\theta }\\ \mathrm{col}\left[\mathit{\gamma }\right]\end{array}\right]\in {\mathrm{ℝ}}^{n_0}.\hfill \end{array}$$

(8)

However, the information matrix $\mathit{\Phi }\left(t\right)$ in model (6) contains a large number of zero elements because it is calculated from the Kronecker product, which results in redundant computation amounts in the identification processes. Therefore, it is necessary to find an alternative method with less computation amounts to estimate the identification model in (5). In this paper, two efficient algorithms with good performances are proposed to solve this problem by applying the coupling identification concept and the multi-innovation identification theory.

Suppose that $\widehat{\mathit{X}}\left(t\right)$ is the estimate of $\mathit{X}$ at time t. That is to say $\widehat{c_i}\left(t\right)$, $\widehat{\mathit{\gamma }}\left(t\right)$ and $\widehat{\mathit{w}}\left(t\right)$ are estimates of $c_i$, $\mathit{\gamma }$ and $\mathit{w}$ at time t, respectively.

3. The Partially Coupled Stochastic Gradient Algorithm

First of all, according to the number of the system outputs, decomposing the identification model in (5) into m sub-identification models:

$$\left[\begin{array}{c}{y}_1\left(t\right)\\ y_2\left(t\right)\\ ⋮\\ y_m\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathit{\psi }}_1^{\mathrm{T}}\left(t\right)\\ {\mathit{\psi }}_2^{\mathrm{T}}\left(t\right)\\ ⋮\\ {\mathit{\psi }}_m^{\mathrm{T}}\left(t\right)\end{array}\right]\mathit{\theta }+\left[\begin{array}{c}{\mathit{\gamma }}_1^{\mathrm{T}}\\ {\mathit{\gamma }}_2^{\mathrm{T}}\\ ⋮\\ {\mathit{\gamma }}_m^{\mathrm{T}}\end{array}\right]\mathit{\varphi }\left(t\right)+\left[\begin{array}{c}{v}_1\left(t\right)\\ v_2\left(t\right)\\ ⋮\\ v_m\left(t\right)\end{array}\right],i=1,2,\cdots ,m,$$

(9)

where ${\mathit{\gamma }}_i^{\mathrm{T}}\in {\mathrm{ℝ}}^{1\times m{n}_c}$ are the ith row of the parameter matrix ${\mathit{\gamma }}^{\mathrm{T}}$, and ${\mathit{\psi }}_i^{\mathrm{T}}\left(t\right)\in {\mathrm{ℝ}}^{1\times n}$ are the ith row of the information matrix ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$:

$$\begin{array}{ccc}\hfill {\mathit{\gamma }}^{\mathrm{T}}& :=& {[{\mathit{\gamma }}_1,{\mathit{\gamma }}_2,\cdots ,{\mathit{\gamma }}_m]}^{\mathrm{T}}\in {\mathrm{ℝ}}^{m\times \left(m{n}_c\right)},\hfill \\ \hfill {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)& :=& {[{\mathit{\psi }}_1\left(t\right),{\mathit{\psi }}_2\left(t\right),\cdots ,{\mathit{\psi }}_m\left(t\right)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^{m\times n}.\hfill \end{array}$$

For the m sub-identification models in (9), it can be seen that the information vector $\mathit{\theta }$ and the parameter vector $\mathit{\varphi }\left(t\right)$ are common among the m subsystems. Thus, Equation (5) is a coupling identification model of partial parameter vectors and partial information vectors. Equation (9) can be represented as

$$\begin{array}{ccc}\hfill y_i\left(t\right)& =& {\mathit{\psi }}_i^{\mathrm{T}}\left(t\right)\mathit{\theta }+{\mathit{\gamma }}_i^{\mathrm{T}}\mathit{\varphi }\left(t\right)+v_i\left(t\right)\hfill \\ & =& {\mathit{\psi }}_i^{\mathrm{T}}\left(t\right)\mathit{\theta }+{\mathit{\varphi }}^{\mathrm{T}}\left(t\right){\mathit{\gamma }}_i+v_i\left(t\right),i=1,2,\cdots ,m.\hfill \end{array}$$

(10)

For the m sub-identification models in (10), define a gradient criterion function,

$$J_1(\mathit{\theta },{\mathit{\gamma }}_i):={[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right)\mathit{\theta }-{\mathit{\varphi }}^{\mathrm{T}}\left(t\right){\mathit{\gamma }}_i]}^2,i=1,2,\cdots ,m.$$

Assuming that $1/r_{\theta ,i}\left(t\right)$ and $1/r_{\gamma ,i}\left(t\right)(i=1,2,\cdots ,m)$ are the step-size, using the negative gradient search [33] and minimizing $J_1(\mathit{\theta },{\mathit{\gamma }}_i)$, we have the gradient relationships:

$$\begin{array}{ccc}\hfill \widehat{\mathit{\theta }}\left(t\right)& =& \widehat{\mathit{\theta }}(t-1)+\frac{{\mathit{\psi }}_i\left(t\right)}{r_{\theta ,i}\left(t\right)}[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right)\widehat{\mathit{\theta }}(t-1)-{\mathit{\varphi }}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)],\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill r_{\theta ,i}\left(t\right)& =& r_{\theta ,i}(t-1)+{\parallel {\mathit{\psi }}_i\left(t\right)\parallel }^2,r_{\theta ,i}\left(0\right)=1,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\gamma }}}_i\left(t\right)& =& {\widehat{\mathit{\gamma }}}_i(t-1)+\frac{\mathit{\varphi }\left(t\right)}{r_{\gamma ,i}\left(t\right)}[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right)\widehat{\mathit{\theta }}(t-1)-{\mathit{\varphi }}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)],\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill r_{\gamma ,i}\left(t\right)& =& r_{\gamma ,i}(t-1)+{\parallel \mathit{\varphi }\left(t\right)\parallel }^2,r_{\gamma ,i}\left(0\right)=1.\hfill \end{array}$$

(14)

The problem of identification in (11)–(14) is that the estimates $\widehat{\mathit{\theta }}\left(t\right)$ and $\widehat{{\mathit{\gamma }}_i}\left(t\right)$ can not be computed because $\mathit{\varphi }\left(t\right)$ contains the unmeasurable noise terms $\mathit{w}(t-i)$. The method to solve this problem is to replace the unmeasurable variables $\mathit{w}(t-i)$ with their corresponding estimates $\widehat{\mathit{w}}(t-i)$. Thus, the estimate of $\mathit{\varphi }\left(t\right)$ can be computed by $\widehat{\mathit{w}}(t-i)$:

$$\begin{array}{ccc}\hfill \widehat{\mathit{\varphi }}\left(t\right)& :=& {[-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-1),-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-2),\cdots ,-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-n_c)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^{m{n}_c}.\hfill \end{array}$$

(15)

It can be seen that the parameter vector $\mathit{\theta }$ is estimated for m times in (11) and (12), which will lead to many redundant estimates. Using the estimate ${\widehat{\mathit{\theta }}}_i$ to replace the estimate $\widehat{\mathit{\theta }}$ in (11) can reduce the redundant computation. At the same time, replacing the unknown information vector $\mathit{\varphi }\left(t\right)$ with its estimate $\widehat{\mathit{\varphi }}\left(t\right)$, we can obtain the new gradient relationships:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\theta }}}_i\left(t\right)& =& {\widehat{\mathit{\theta }}}_i(t-1)+\frac{{\mathit{\psi }}_i\left(t\right)}{r_{\theta ,i}\left(t\right)}[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_i(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)],\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill r_{\theta ,i}\left(t\right)& =& r_{\theta ,i}(t-1)+{\parallel {\mathit{\psi }}_i\left(t\right)\parallel }^2,r_{\theta ,i}\left(0\right)=1,\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\gamma }}}_i\left(t\right)& =& {\widehat{\mathit{\gamma }}}_i(t-1)+\frac{\widehat{\mathit{\varphi }}\left(t\right)}{r_{\gamma ,i}\left(t\right)}[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_i(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)],\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill r_{\gamma ,i}\left(t\right)& =& r_{\gamma ,i}(t-1)+{\parallel \widehat{\mathit{\varphi }}\left(t\right)\parallel }^2,r_{\gamma ,i}\left(0\right)=1.\hfill \end{array}$$

(19)

Additionally, according to Equations (4), we have

$$\begin{array}{ccc}\hfill \widehat{\mathit{w}}\left(t\right)& =& \mathit{y}\left(t\right)-{\mathit{\Phi }}_s\left(t\right){\widehat{\mathit{\theta }}}_m\left(t\right).\hfill \end{array}$$

(20)

It is generally believed that the parameter estimate ${\widehat{\mathit{\theta }}}_{i-1}\left(t\right)$ of the $(i-1)$th subsystem at time t is closer to the true value $\mathit{\theta }$ than the parameter estimate ${\widehat{\mathit{\theta }}}_i(t-1)$ of the ith subsystem at time $(t-1)$. Replacing ${\widehat{\mathit{\theta }}}_i(t-1)$ on the right-hand of Equation (16) with ${\widehat{\mathit{\theta }}}_{i-1}\left(t\right)$, and replacing ${\widehat{\mathit{\theta }}}_1(t-1)$ of Equation (16) when $i=1$ with ${\widehat{\mathit{\theta }}}_m(t-1)$. Combing Equations (15) and (20), we can obtain the multivariate partially coupled generalized stochastic gradient (M-PC-GSG) algorithm as

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\theta }}}_1\left(t\right)& =& {\widehat{\mathit{\theta }}}_m(t-1)+\frac{{\mathit{\psi }}_1\left(t\right)}{r_{\theta ,1}\left(t\right)}[y_1\left(t\right)-{\mathit{\psi }}_1^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_1(t-1)],\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill r_{\theta ,1}\left(t\right)& =& r_{\theta ,m}(t-1)+{\parallel {\mathit{\psi }}_1\left(t\right)\parallel }^2,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\gamma }}}_1\left(t\right)& =& {\widehat{\mathit{\gamma }}}_1(t-1)+\frac{\widehat{\mathit{\varphi }}\left(t\right)}{r_{\gamma ,1}\left(t\right)}[y_1\left(t\right)-{\mathit{\psi }}_1^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_1(t-1)],\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill r_{\gamma ,1}\left(t\right)& =& r_{\gamma ,m}(t-1)+{\parallel \widehat{\mathit{\varphi }}\left(t\right)\parallel }^2,\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\theta }}}_i\left(t\right)& =& {\widehat{\mathit{\theta }}}_{i-1}\left(t\right)+\frac{{\mathit{\psi }}_i\left(t\right)}{r_{\theta ,i}\left(t\right)}[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)],\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill r_{\theta ,i}\left(t\right)& =& r_{\theta ,i-1}\left(t\right)+{\parallel {\mathit{\psi }}_i\left(t\right)\parallel }^2,\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\gamma }}}_i\left(t\right)& =& {\widehat{\mathit{\gamma }}}_i(t-1)+\frac{\widehat{\mathit{\varphi }}\left(t\right)}{r_{\gamma ,i}\left(t\right)}[y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)],\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill r_{\gamma ,i}\left(t\right)& =& r_{\gamma ,i-1}\left(t\right)+{\parallel \widehat{\mathit{\varphi }}\left(t\right)\parallel }^2,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)& =& {[{\mathit{\psi }}_1\left(t\right),{\mathit{\psi }}_2\left(t\right),\cdots ,{\mathit{\psi }}_m\left(t\right)]}^{\mathrm{T}},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \widehat{\mathit{\varphi }}\left(t\right)& =& {[-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-1),-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-2),\cdots ,-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-n_c)]}^{\mathrm{T}},\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill \widehat{\mathit{w}}\left(t\right)& =& \mathit{y}\left(t\right)-{\mathit{\Phi }}_s\left(t\right){\widehat{\mathit{\theta }}}_m\left(t\right),\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill \widehat{\mathit{\gamma }}\left(t\right)& =& [{\widehat{\mathit{\gamma }}}_1\left(t\right),{\widehat{\mathit{\gamma }}}_2\left(t\right),\cdots ,{\widehat{\mathit{\gamma }}}_m\left(t\right)].\hfill \end{array}$$

(32)

The schematic diagram of the M-PC-GSG algorithm in (21)–(32) is shown in Figure 1, and the computation procedures are listed as follows.

• Let $t=1$, set the initial values ${\widehat{\mathit{\theta }}}_m\left(0\right)={\mathit{1}}_n/p_0$, $r_{\theta ,i}\left(0\right)=r_{\gamma ,i}\left(0\right)=1$, ${\widehat{\mathit{\gamma }}}_i\left(0\right)={\mathit{1}}_{m{n}_c}/p_0$, $i=1,2,\cdots ,m$, $\widehat{\mathit{w}}(t-j)=\mathit{0}$, $j=0$, 1, ⋯, $n_c$, $p_0={10}^6$, and set the data length K.
• Collect the observation data ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$ and $\mathit{y}\left(t\right)$, read ${\mathit{\psi }}_i\left(t\right)$ from ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$ using (29).
• Construct $\widehat{\mathit{\varphi }}\left(t\right)$ using (30).
• Compute the step-size $r_{\theta ,1}\left(t\right)$ and $r_{\gamma ,1}\left(t\right)$ using (22) and (24).
• Update the parameter estimates ${\widehat{\mathit{\theta }}}_1\left(t\right)$ and $\widehat{{\mathit{\gamma }}_1}\left(t\right)$ using (21) and (23).
• When $i=2,3,\cdots ,m$, compute $r_{\theta ,i}\left(t\right)$ and $r_{\gamma ,i}\left(t\right)$ using (26) and (28), and update the parameter estimates ${\widehat{\mathit{\theta }}}_i\left(t\right)$ and ${\widehat{\mathit{\gamma }}}_i\left(t\right)$ using (25) and (27).
• Compute $\widehat{\mathit{w}}\left(t\right)$ using (31) and construct $\widehat{\mathit{\gamma }}\left(t\right)$ using (32).
• Increase t by 1 if $t<K$, and then go to Step 2. Otherwise, obtain parameter estimates $\widehat{\mathit{\theta }}\left(t\right)$ and $\widehat{\mathit{\gamma }}\left(t\right)$ and stop computing.

4. The Partially Coupled Multi-Innovation Stochastic Gradient Algorithm

In this section, we apply the multi-innovation identification theory based on the M-PC-GSG algorithm to further improve the parameter estimation accuracy. We introduce an innovation length $p⩾1$ to expand the innovation scalars of each subsystem models to the innovation vectors. According to the M-PC-GSG algorithm in (21)–(32), define the subsystem stacked information matrices ${\mathit{\Lambda }}_1(p,t)$, ${\mathit{\Lambda }}_i(p,t)$$(i=2,3,\cdots ,m)$, $\mathit{\Xi }(p,t)$ and the subsystem stacked output vectors ${\mathit{Y}}_1(p,t)$ and ${\mathit{Y}}_i(p,t)$$(i=2,3,\cdots ,m)$ as

$$\begin{array}{ccc}\hfill {\mathit{\Lambda }}_1(p,t)& :=& [{\mathit{\psi }}_1\left(t\right),{\mathit{\psi }}_1(t-1),\cdots ,{\mathit{\psi }}_1(t-p+1)]\in {\mathrm{ℝ}}^{n\times p},\hfill \\ \hfill {\mathit{\Lambda }}_i(p,t)& :=& [{\mathit{\psi }}_i\left(t\right),{\mathit{\psi }}_i(t-1),\cdots ,{\mathit{\psi }}_i(t-p+1)]\in {\mathrm{ℝ}}^{n\times p},i=2,3,\cdots ,m,\hfill \\ \hfill \mathit{\Xi }(p,t)& :=& [\widehat{\mathit{\varphi }}\left(t\right),\widehat{\mathit{\varphi }}(t-1),\cdots ,\widehat{\mathit{\varphi }}(t-p+1)]\in {\mathrm{ℝ}}^{m{n}_c\times p},\hfill \\ \hfill {\mathit{Y}}_1(p,t)& :=& {[y_1\left(t\right),y_1(t-1),\cdots ,y_1(t-p+1)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^p,\hfill \\ \hfill {\mathit{Y}}_i(p,t)& :=& {[y_i\left(t\right),y_i(t-1),\cdots ,y_i(t-p+1)]}^{\mathrm{T}}\in {\mathrm{ℝ}}^p,i=2,3,\cdots ,m.\hfill \end{array}$$

In general, define the subsystem innovation scalars $e_1\left(t\right)$ and $e_i\left(t\right)$$(i=2,3,\cdots ,m)$ in the M-PC-GSG algorithm in (21)–(32) as

$$\begin{array}{ccc}\hfill e_1\left(t\right)& :=& y_1\left(t\right)-{\mathit{\psi }}_1^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_1(t-1)\in \mathrm{ℝ},\hfill \\ \hfill e_i\left(t\right)& :=& y_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)\in ,\mathrm{ℝ}i=2,3,\cdots ,m.\hfill \end{array}$$

According to the multi-innovation identification theory, expand the subsystem innovation scalars $e_1\left(t\right)$ and $e_i\left(t\right)$$(i=2,3,\cdots ,m)$ into the subsystem innovation vectors ${\mathit{E}}_1(p,t)$ and ${\mathit{E}}_i(p,t)$$(i=2,3,\cdots ,m)$:

$$\begin{array}{ccc}\hfill {\mathit{E}}_1(p,t)& :=& \left[\begin{array}{c}{e}_1\left(t\right)\\ e_1(t-1)\\ ⋮\\ e_1(t-p+1)\end{array}\right]\hfill \\ & =& \left[\begin{array}{c}{y}_1\left(t\right)-{\mathit{\psi }}_1^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_1(t-1)\\ y_1(t-1)-{\mathit{\psi }}_1^{\mathrm{T}}(t-1){\widehat{\mathit{\theta }}}_m(t-2)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-1){\widehat{\mathit{\gamma }}}_1(t-2)\\ ⋮\\ y_1(t-p+1)-{\mathit{\psi }}_1^{\mathrm{T}}(t-p+1){\widehat{\mathit{\theta }}}_m(t-p)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-p+1){\widehat{\mathit{\gamma }}}_1(t-p)\end{array}\right]\in {\mathrm{ℝ}}^p,\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\mathit{E}}_i(p,t)& :=& \left[\begin{array}{c}{e}_i\left(t\right)\\ e_i(t-1)\\ ⋮\\ e_i(t-p+1)\end{array}\right](i=2,3,\cdots ,m)\hfill \\ & =& \left[\begin{array}{c}{y}_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)\\ y_i(t-1)-{\mathit{\psi }}_i^{\mathrm{T}}(t-1){\widehat{\mathit{\theta }}}_{i-1}(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-1){\widehat{\mathit{\gamma }}}_i(t-2)\\ ⋮\\ y_i(t-p+1)-{\mathit{\psi }}_i^{\mathrm{T}}(t-p+1){\widehat{\mathit{\theta }}}_{i-1}(t-p+1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-p+1){\widehat{\mathit{\gamma }}}_i(t-p)\end{array}\right]\in {\mathrm{ℝ}}^p.\hfill \end{array}$$

(34)

Normally, we reach an agreement that the estimates ${\widehat{\mathit{\theta }}}_m(t-1)$ and ${\widehat{\mathit{\gamma }}}_i(t-1)$$(i=1,2,\cdots ,m)$ at time $(t-1)$ are closer to the true values ${\mathit{\theta }}_m$ and ${\mathit{\gamma }}_i$ than the estimates ${\widehat{\mathit{\theta }}}_m(t-j)$ and ${\widehat{\mathit{\gamma }}}_i(t-j)$ at time $(t-j)$$(j⩾2)$. Similarly, the estimate ${\widehat{\mathit{\theta }}}_{i-1}\left(t\right)$ at time t is closer to the true value ${\mathit{\theta }}_{i-1}$ than the estimate ${\widehat{\mathit{\theta }}}_{i-1}(t-j)$ at time $(t-j)$$(j⩾2)$. Therefore, replacing the terms ${\widehat{\mathit{\theta }}}_m(t-j)$, ${\widehat{\mathit{\gamma }}}_i(t-j)$ and ${\widehat{\mathit{\theta }}}_{i-1}(t-j)$ ($j⩾2$) with ${\widehat{\mathit{\theta }}}_m(t-1)$, ${\widehat{\mathit{\gamma }}}_i(t-1)$ and ${\widehat{\mathit{\theta }}}_{i-1}\left(t\right)$ in (33) and (34), and the subsystem innovation vectors ${\mathit{E}}_1(p,t)$ and ${\mathit{E}}_i(p,t)$$(i=2,3,\cdots ,m)$ are modified into

$$\begin{array}{ccc}\hfill {\mathit{E}}_1(p,t)& :=& \left[\begin{array}{c}{y}_1\left(t\right)-{\mathit{\psi }}_1^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_1(t-1)\\ y_1(t-1)-{\mathit{\psi }}_1^{\mathrm{T}}(t-1){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-1){\widehat{\mathit{\gamma }}}_1(t-1)\\ ⋮\\ y_1(t-p+1)-{\mathit{\psi }}_1^{\mathrm{T}}(t-p+1){\widehat{\mathit{\theta }}}_m(t-1)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-p+1){\widehat{\mathit{\gamma }}}_1(t-1)\end{array}\right]\hfill \\ & =& {\mathit{Y}}_1(p,t)-{\mathit{\Lambda }}_1^{\mathrm{T}}(p,t){\widehat{\mathit{\theta }}}_m(t-1)-{\mathit{\Xi }}^{\mathrm{T}}(p,t){\widehat{\mathit{\gamma }}}_1(t-1)\in {\mathrm{ℝ}}^p,\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill {\mathit{E}}_i(p,t)& :=& \left[\begin{array}{c}{y}_i\left(t\right)-{\mathit{\psi }}_i^{\mathrm{T}}\left(t\right){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}\left(t\right){\widehat{\mathit{\gamma }}}_i(t-1)\\ y_i(t-1)-{\mathit{\psi }}_i^{\mathrm{T}}(t-1){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-1){\widehat{\mathit{\gamma }}}_i(t-1)\\ ⋮\\ y_i(t-p+1)-{\mathit{\psi }}_i^{\mathrm{T}}(t-p+1){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\widehat{\mathit{\varphi }}}^{\mathrm{T}}(t-p+1){\widehat{\mathit{\gamma }}}_i(t-1)\end{array}\right],\hfill \\ & =& {\mathit{Y}}_i(p,t)-{\mathit{\Lambda }}_i^{\mathrm{T}}(p,t){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\mathit{\Xi }}^{\mathrm{T}}(p,t){\widehat{\mathit{\gamma }}}_i(t-1)\in {\mathrm{ℝ}}^p,i=2,3,\cdots ,m.\hfill \end{array}$$

(36)

Thus, based on the the M-PC-GSG algorithm in (21)–(32), we can obtain the following multivariate partially coupled multi-innovation generalized stochastic gradient (M-PC-MI-GSG) algorithm with an innovation length p:

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\theta }}}_1\left(t\right)& =& {\widehat{\mathit{\theta }}}_m(t-1)+\frac{{\mathit{\Lambda }}_1(p,t)}{r_{\theta ,1}\left(t\right)}{\mathit{E}}_1(p,t),\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill {\mathit{Y}}_1(p,t)& =& {[y_1\left(t\right),y_1(t-1),\cdots ,y_1(t-p+1)]}^{\mathrm{T}},\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill {\mathit{\Lambda }}_1(p,t)& =& [{\mathit{\psi }}_1\left(t\right),{\mathit{\psi }}_1(t-1),\cdots ,{\mathit{\psi }}_1(t-p+1)],\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill \mathit{\Xi }(p,t)& =& [\widehat{\mathit{\varphi }}\left(t\right),\widehat{\mathit{\varphi }}(t-1),\cdots ,\widehat{\mathit{\varphi }}(t-p+1)],\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill {\mathit{E}}_1(p,t)& =& {\mathit{Y}}_1(p,t)-{\mathit{\Lambda }}_1^{\mathrm{T}}(p,t){\widehat{\mathit{\theta }}}_m(t-1)-{\mathit{\Xi }}^{\mathrm{T}}(p,t){\widehat{\mathit{\gamma }}}_1(t-1),\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill r_{\theta ,1}\left(t\right)& =& r_{\theta ,m}(t-1)+{\parallel {\mathit{\psi }}_1\left(t\right)\parallel }^2,\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\gamma }}}_1\left(t\right)& =& {\widehat{\mathit{\gamma }}}_1(t-1)+\frac{\mathit{\Xi }(p,t)}{r_{\gamma ,1}\left(t\right)}{\mathit{E}}_1(p,t),\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill r_{\gamma ,1}\left(t\right)& =& r_{\gamma ,m}(t-1)+{\parallel \widehat{\mathit{\varphi }}\left(t\right)\parallel }^2,\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\theta }}}_i\left(t\right)& =& {\widehat{\mathit{\theta }}}_{i-1}\left(t\right)+\frac{{\mathit{\Lambda }}_i(p,t)}{r_{\theta ,i}\left(t\right)}{\mathit{E}}_i(p,t),\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill {\mathit{Y}}_i(p,t)& =& {[y_i\left(t\right),y_i(t-1),\cdots ,y_i(t-p+1)]}^{\mathrm{T}},\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill {\mathit{\Lambda }}_i(p,t)& =& [{\mathit{\psi }}_i\left(t\right),{\mathit{\psi }}_i(t-1),\cdots ,{\mathit{\psi }}_i(t-p+1)],\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill {\mathit{E}}_i(p,t)& =& {\mathit{Y}}_i(p,t)-{\mathit{\Lambda }}_i^{\mathrm{T}}(p,t){\widehat{\mathit{\theta }}}_{i-1}\left(t\right)-{\mathit{\Xi }}^{\mathrm{T}}(p,t){\widehat{\mathit{\gamma }}}_i(t-1),\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill r_{\theta ,i}\left(t\right)& =& r_{\theta ,i-1}\left(t\right)+{\parallel {\mathit{\psi }}_i\left(t\right)\parallel }^2,\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill {\widehat{\mathit{\gamma }}}_i\left(t\right)& =& {\widehat{\mathit{\gamma }}}_i(t-1)+\frac{\mathit{\Xi }(p,t)}{r_{\gamma ,i}\left(t\right)}{\mathit{E}}_i(p,t),\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill r_{\gamma ,i}\left(t\right)& =& r_{\gamma ,i-1}\left(t\right)+{\parallel \widehat{\mathit{\varphi }}\left(t\right)\parallel }^2,\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)& =& {[{\mathit{\psi }}_1\left(t\right),{\mathit{\psi }}_2\left(t\right),\cdots ,{\mathit{\psi }}_m\left(t\right)]}^{\mathrm{T}},\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill \widehat{\mathit{\varphi }}\left(t\right)& =& {[-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-1),-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-2),\cdots ,-{\widehat{\mathit{w}}}^{\mathrm{T}}(t-n_c)]}^{\mathrm{T}},\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill \widehat{\mathit{w}}\left(t\right)& =& \mathit{y}\left(t\right)-{\mathit{\Phi }}_s\left(t\right){\widehat{\mathit{\theta }}}_m\left(t\right),\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill \widehat{\mathit{\gamma }}\left(t\right)& =& [{\widehat{\mathit{\gamma }}}_1\left(t\right),{\widehat{\mathit{\gamma }}}_2\left(t\right),\cdots ,{\widehat{\mathit{\gamma }}}_m\left(t\right)].\hfill \end{array}$$

(55)

The computation procedures of the M-PC-MI-GSG algorithm in (37)–(55) are listed as follows.

• Let $t=1$, choose an innovation length p, set the initial values ${\widehat{\mathit{\theta }}}_m\left(0\right)={\mathit{1}}_n/p_0$, $r_{\theta ,i}\left(0\right)=r_{\gamma ,i}\left(0\right)=1$, ${\widehat{\mathit{\gamma }}}_i\left(0\right)={\mathit{1}}_{m{n}_c}/p_0$, $i=1,2,\cdots ,m$, $\widehat{\mathit{w}}(t-j)=\mathit{0}$, $j=0$, 1, ⋯, $n_c$, $p_0={10}^6$, and set the data length K.
• Collect the observation data ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$ and $\mathit{y}\left(t\right)$, read ${\mathit{\psi }}_i\left(t\right)$ from ${\mathit{\Phi }}_{\mathrm{s}}\left(t\right)$ using (52).
• Construct $\widehat{\mathit{\varphi }}\left(t\right)$ using (53).
• Construct ${\mathit{Y}}_1(p,t)$, ${\mathit{\Lambda }}_1(p,t)$ and $\mathit{\Xi }(p,t)$ using (38)–(40).
• Compute ${\mathit{E}}_1(p,t)$ using (41).
• Compute $r_{\theta ,1}\left(t\right)$ and $r_{\gamma ,1}\left(t\right)$ using (42) and (44).
• Update the parameter estimates ${\widehat{\mathit{\theta }}}_1\left(t\right)$ and $\widehat{{\mathit{\gamma }}_1}\left(t\right)$ using (37) and (43).
• When $i=2,3,\cdots ,m$, construct ${\mathit{Y}}_i(p,t)$ and ${\mathit{\Lambda }}_i(p,t)$ using (46) and (47).
• Compute ${\mathit{E}}_i(p,t)$ using (48).
• Compute $r_{\theta ,i}\left(t\right)$ and $r_{\gamma ,i}\left(t\right)$ using (49) and (51).
• Update the parameter estimates ${\widehat{\mathit{\theta }}}_i\left(t\right)$ and ${\widehat{\mathit{\gamma }}}_i\left(t\right)$ using (45) and (50).
• Compute $\widehat{\mathit{w}}\left(t\right)$ using (54) and construct $\widehat{\mathit{\gamma }}\left(t\right)$ using (55).
• Increase t by 1 if $t<K$, and then go to Step 2. Otherwise, obtain parameter estimates $\widehat{\mathit{\theta }}\left(t\right)$ and $\widehat{\mathit{\gamma }}\left(t\right)$ and stop computing.

It can be easily seen that the M-PC-MI-GSG algorithm with the innovation length $p=1$ is equal to the M-PC-GSG algorithm. Since the past data information of the system are utilized sufficiently, the M-PC-MI-GSG algorithm has higher parameter estimation accuracy than the M-PC-GSG algorithm.

5. The Simulation Examples

In this section, two numerical simulations are given to show the good performances of the newly proposed algorithms.

Example 1.

Consider the following multivariate equation-error autoregressive systems,

$$\begin{array}{ccc}\hfill \mathit{y}\left(t\right)& =& {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)\mathit{\theta }+{\mathit{C}}^{-1}\left(z\right)\mathit{v}\left(t\right),\hfill \\ \hfill {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)& =& \left[\begin{array}{cc}sin\left(u_1(t-2)\right)y_1(t-1)& cos\left(u_1(t-1)\right)y_1(t-2)\\ cos\left(u_1(t-2)\right)y_2(t-1)& sin\left(u_1(t-2)\right)y_2(t-2)\end{array}\right.\hfill \\ & & \left.\begin{array}{cc}sin\left(u_2(t-2)\right)y_2(t-1)& cos\left(u_2(t-1)\right)y_2(t-2)\\ cos\left(u_2(t-2)\right)y_1(t-1)& sin\left(u_2(t-1)\right)y_1(t-2)\end{array}\right],\hfill \\ \hfill \mathit{\theta }& =& {[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4]}^{\mathrm{T}}={[0.42,0.93,0.56,0.31]}^{\mathrm{T}},\hfill \\ \hfill \mathit{C}\left(z\right)& =& {\mathit{I}}_2+\left[\begin{array}{cc}\hfill c_{11}& \hfill c_{12}\\ \hfill c_{21}& \hfill c_{22}\end{array}\right]z^{-1}=\left[\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right]+\left[\begin{array}{cc}\hfill -0.25& \hfill 0.68\\ \hfill -0.33& \hfill 0.44\end{array}\right]z^{-1}.\hfill \end{array}$$

The parameter vector to be estimated is

$$\begin{array}{ccc}\hfill \mathit{\vartheta }& =& {[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,c_{11},c_{12},c_{21},c_{22}]}^{\mathrm{T}}\hfill \\ & =& {[0.42,0.93,0.56,0.31,-0.25,0.68,-0.33,0.44]}^{\mathrm{T}}.\hfill \end{array}$$

In this simulation, $\mathit{u}\left(t\right)=\left[\begin{array}{c}{u}_1\left(t\right)\\ u_2\left(t\right)\end{array}\right]\in {\mathrm{ℝ}}^2$ is the input vector, which is a random sequence with zero mean and variance one; $\mathit{y}\left(t\right)=\left[\begin{array}{c}{y}_1\left(t\right)\\ y_2\left(t\right)\end{array}\right]\in {\mathrm{ℝ}}^2$ is the output vector; $\mathit{v}\left(t\right)=\left[\begin{array}{c}{v}_1\left(t\right)\\ v_2\left(t\right)\end{array}\right]\in {\mathrm{ℝ}}^2$ is a white noise vector with zero mean; ${\sigma }_1^2$ and ${\sigma }_2^2$ are the variances of $v_1\left(t\right)$ and $v_2\left(t\right)$. Taking the noise variances ${\sigma }_1^2=0.{40}^2$ and ${\sigma }_2^2=0.{30}^2$, using the M-PC-GSG algorithm (i.e., the M-PC-MI-GSG algorithm with $p=1$) and the M-PC-MI-GSG algorithm with $p=2$, $p=4$ and $p=8$ to estimate the parameters of this example system. We obtain the parameter estimates and their errors $\delta :=\parallel \widehat{\mathit{\vartheta }}\left(t\right)-\mathit{\vartheta }\parallel /\parallel \mathit{\vartheta }\parallel$ shown in

Table 1. Parameter estimates and errors (${\sigma }_1^2={0.40}^2,{\sigma }_2^2={0.30}^2$).

Algorithms	t	100	200	500	1000	2000	3000	True Values
M-PC-GSG	${\theta }_1$	0.19022	0.28281	0.33788	0.36637	0.36787	0.36879	0.42000
	${\theta }_2$	0.77045	0.84529	0.88400	0.91050	0.90956	0.90832	0.93000
	${\theta }_3$	0.42624	0.49137	0.54550	0.57360	0.57695	0.57314	0.56000
	${\theta }_4$	0.25045	0.25999	0.29032	0.31802	0.32379	0.32917	0.31000
	$c_{11}$	−0.16279	−0.13896	−0.12959	−0.13205	−0.13608	−0.13540	−0.25000
	$c_{12}$	−0.06264	0.01024	0.07393	0.10879	0.14383	0.16385	0.68000
	$c_{21}$	−0.18783	−0.21131	−0.22760	−0.22982	−0.23502	−0.23697	−0.33000
	$c_{22}$	0.06255	0.05905	0.07305	0.09172	0.10850	0.12021	0.44000
	$\delta (\%)$	60.04604	53.54881	48.52141	45.66110	43.08714	41.59002
M-PC-MI-GSG	${\theta }_1$	0.24413	0.35242	0.38045	0.39205	0.38805	0.38809	0.42000
$(p=2)$	${\theta }_2$	0.80411	0.89433	0.91071	0.92629	0.92058	0.91797	0.93000
	${\theta }_3$	0.49437	0.54833	0.57802	0.58115	0.58224	0.57102	0.56000
	${\theta }_4$	0.30721	0.28679	0.29391	0.32301	0.32347	0.32964	0.31000
	$c_{11}$	−0.05715	−0.06016	−0.08608	−0.10747	−0.13040	−0.13492	−0.25000
	$c_{12}$	0.12099	0.22052	0.30492	0.34234	0.38260	0.40361	0.68000
	$c_{21}$	−0.23138	−0.26019	−0.26756	−0.26447	−0.26951	−0.27084	−0.33000
	$c_{22}$	0.11762	0.11798	0.15107	0.18454	0.20984	0.22701	0.44000
	$\delta (\%)$	47.51695	39.87369	33.61530	30.01753	26.59647	24.80351
M-PC-MI-GSG	${\theta }_1$	0.26824	0.39655	0.39206	0.40527	0.40117	0.40344	0.42000
$(p=4)$	${\theta }_2$	0.81764	0.91822	0.91009	0.93027	0.92223	0.92506	0.93000
	${\theta }_3$	0.53011	0.56224	0.57291	0.58006	0.58163	0.56221	0.56000
	${\theta }_4$	0.33781	0.28934	0.27962	0.32670	0.31900	0.33021	0.31000
	$c_{11}$	−0.02859	−0.07525	−0.14398	−0.17985	−0.21158	−0.20979	−0.25000
	$c_{12}$	0.40063	0.49444	0.57042	0.57794	0.60699	0.61726	0.68000
	$c_{21}$	−0.27689	−0.30473	−0.29946	−0.29170	−0.29725	−0.29909	−0.33000
	$c_{22}$	0.22503	0.21310	0.25734	0.29682	0.32439	0.34271	0.44000
	$\delta (\%)$	30.60589	22.77807	16.20218	12.94875	9.87303	8.55799
M-PC-MI-GSG	${\theta }_1$	0.22875	0.41123	0.39312	0.41067	0.40847	0.41377	0.42000
$(p=8)$	${\theta }_2$	0.82538	0.91675	0.89962	0.92660	0.91916	0.93461	0.93000
	${\theta }_3$	0.50062	0.56129	0.56967	0.58480	0.57905	0.55554	0.56000
	${\theta }_4$	0.34784	0.29455	0.27015	0.32493	0.31639	0.33324	0.31000
	$c_{11}$	−0.08032	−0.14448	−0.22262	−0.25289	−0.27162	−0.25101	−0.25000
	$c_{12}$	0.60874	0.66577	0.71758	0.67480	0.69550	0.69273	0.68000
	$c_{21}$	−0.31148	−0.33276	−0.30958	−0.29949	−0.30552	−0.31357	−0.33000
	$c_{22}$	0.32292	0.26569	0.31977	0.36227	0.39067	0.40749	0.44000
	$\delta (\%)$	20.99172	13.61264	9.45336	5.90521	4.39162	3.04341

. The parameter estimation errors versus t are shown in Figure 2. The parameter estimates ${\widehat{\theta }}_1\left(t\right)$, ${\widehat{\theta }}_2\left(t\right)$, ${\widehat{\theta }}_3\left(t\right)$, ${\widehat{\theta }}_4\left(t\right)$ and ${\widehat{c}}_{11}\left(t\right)$, ${\widehat{c}}_{12}\left(t\right)$, ${\widehat{c}}_{21}\left(t\right)$, ${\widehat{c}}_{22}\left(t\right)$ versus t of the M-PC-MI-GSG algorithm with $p=8$ are shown in Figure 3 and Figure 4.

Example 2.

Consider another multivariate equation-error autoregressive systems,

$$\begin{array}{ccc}\hfill \mathit{y}\left(t\right)& =& {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)\mathit{\theta }+{\mathit{C}}^{-1}\left(z\right)\mathit{v}\left(t\right),\hfill \\ \hfill {\mathit{\Phi }}_{\mathrm{s}}\left(t\right)& =& \left[\begin{array}{ccc}{u}_1(t-1)y_1(t-2)& u_1(t-2)u_2(t-2)& u_2(t-1)cos\left(t\right)\\ u_1(t-2)y_1(t-1)& y_1(t-2)y_2(t-2)& y_2(t-1)sin\left(t\right)\end{array}\right.\hfill \\ & & \left.\begin{array}{ccc}{cos}^2(t-1)& u_2(t-2)sin(t-1)& cos\left(u_1(t-1)\right)\\ {sin}^2(t-1)& y_2(t-2)cos(t-1)& sin\left(y_1(t-1)\right)\end{array}\right],\hfill \\ \hfill \mathit{\theta }& =& {[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6]}^{\mathrm{T}}={[-0.36,0.22,0.34,0.45,0.25,0.11]}^{\mathrm{T}},\hfill \\ \hfill \mathit{C}\left(z\right)& =& {\mathit{I}}_2+\left[\begin{array}{cc}\hfill c_{11}& \hfill c_{12}\\ \hfill c_{21}& \hfill c_{22}\end{array}\right]z^{-1}=\left[\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right]+\left[\begin{array}{cc}\hfill -0.41& \hfill -0.48\\ \hfill 0.35& \hfill -0.31\end{array}\right]z^{-1}.\hfill \end{array}$$

The parameter vector to be estimated is

$$\begin{array}{ccc}\hfill \mathit{\vartheta }& =& {[{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6,c_{11},c_{12},c_{21},c_{22}]}^{\mathrm{T}}\hfill \\ & =& {[-0.36,0.22,0.34,0.45,0.25,0.11,-0.41,-0.48,0.35,-0.31]}^{\mathrm{T}}.\hfill \end{array}$$

Here, the simulation conditions are similar to Example 1. Taking the noise variances ${\sigma }_1^2={\sigma }_2^2={0.50}^2$, using the M-PC-GSG algorithm (i.e., the M-PC-MI-GSG algorithm with $p=1$) and the M-PC-MI-GSG algorithm with $p=3$ and $p=5$ to estimate the parameters of this example system. The parameter estimates and their errors are shown in

Table 2. Parameter estimates and errors (${\sigma }_1^2={\sigma }_2^2={0.50}^2$).

Algorithms	t	100	200	500	1000	2000	3000	True Values
M-PC-GSG	${\theta }_1$	−0.38109	−0.34290	−0.35197	−0.35358	−0.35240	−0.35140	−0.36000
	${\theta }_2$	0.17154	0.19202	0.22318	0.20571	0.20894	0.20290	0.22000
	${\theta }_3$	0.52049	0.46240	0.42063	0.40418	0.38885	0.38130	0.34000
	${\theta }_4$	0.36772	0.39041	0.40431	0.43059	0.44228	0.45084	0.45000
	${\theta }_5$	0.32582	0.28799	0.27265	0.25708	0.26433	0.26468	0.25000
	${\theta }_6$	−0.05956	−0.02720	−0.02478	−0.02008	−0.00934	−0.00876	0.11000
	$c_{11}$	−0.32562	−0.33413	−0.35093	−0.36678	−0.37198	−0.37339	−0.41000
	$c_{12}$	−0.42387	−0.42801	−0.43950	−0.44549	−0.45522	−0.45784	−0.48000
	$c_{21}$	−0.07449	−0.04447	0.00357	0.05541	0.08891	0.10599	0.35000
	$c_{22}$	0.13315	0.09440	0.04896	0.00007	−0.03206	−0.04982	−0.31000
	$\delta (\%)$	62.46630	55.62215	48.66825	41.81593	37.21977	34.99402
M-PC-MI-GSG	${\theta }_1$	−0.34351	−0.30662	−0.36709	−0.35952	−0.35025	−0.35053	−0.36000
$(p=3)$	${\theta }_2$	0.13616	0.21845	0.24620	0.20140	0.21882	0.20575	0.22000
	${\theta }_3$	0.46447	0.34810	0.32286	0.33605	0.34131	0.33736	0.34000
	${\theta }_4$	0.53467	0.51841	0.49198	0.50896	0.49490	0.49143	0.45000
	${\theta }_5$	0.34052	0.24404	0.24150	0.23199	0.26443	0.26536	0.25000
	${\theta }_6$	−0.08458	−0.00717	−0.00412	0.01422	0.05155	0.05032	0.11000
	$c_{11}$	−0.40747	−0.39750	−0.42113	−0.42870	−0.41996	−0.41446	−0.41000
	$c_{12}$	−0.44953	−0.45067	−0.47122	−0.47509	−0.49160	−0.49275	−0.48000
	$c_{21}$	−0.00019	0.05273	0.13849	0.22633	0.25865	0.27774	0.35000
	$c_{22}$	−0.10619	−0.14079	−0.18170	−0.23788	−0.25384	−0.26009	−0.31000
	$\delta (\%)$	45.05598	34.25524	25.50852	16.96953	12.12756	10.74341
M-PC-MI-GSG	${\theta }_1$	−0.35745	−0.30970	−0.38425	−0.36288	−0.34050	−0.34946	−0.36000
$(p=5)$	${\theta }_2$	0.12136	0.23250	0.25401	0.20572	0.22745	0.20850	0.22000
	${\theta }_3$	0.40415	0.28647	0.30611	0.33618	0.34589	0.33781	0.34000
	${\theta }_4$	0.51978	0.48285	0.46617	0.48824	0.46462	0.46024	0.45000
	${\theta }_5$	0.38782	0.22572	0.22904	0.22473	0.27031	0.26832	0.25000
	${\theta }_6$	−0.02008	0.06621	0.05334	0.06901	0.10837	0.09507	0.11000
	$c_{11}$	−0.41690	−0.39157	−0.43170	−0.43710	−0.42419	−0.41720	−0.41000
	$c_{12}$	−0.46440	−0.45497	−0.48430	−0.48122	−0.50181	−0.49994	−0.48000
	$c_{21}$	0.09086	0.14050	0.23086	0.31994	0.33401	0.34857	0.35000
	$c_{22}$	−0.24541	−0.23841	−0.25163	−0.30483	−0.29672	−0.29278	−0.31000
	$\delta (\%)$	32.60246	22.30989	14.47746	6.91170	4.31185	3.73785

. The parameter estimation errors versus t are shown in Figure 5. For model validation, we use the estimated model obtained by the M-PC-MI-GSG algorithm with $p=5$ to predict the system outputs with 200 samples from $t=3001$ to $t=3200$. The true output $y_1\left(t\right)$ and the predicted output ${\widehat{y}}_1\left(t\right)$ as well as their errors are shown in Figure 6. The true output $y_2\left(t\right)$ and the predicted output ${\widehat{y}}_2\left(t\right)$ as well as their errors are shown in Figure 7.

From Table 1 and Table 2 and Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7, the following conclusions are obtained.

• From Table 1 and Table 2, Figure 2 and Figure 5, it can be shown that the parameter estimation errors of the M-PC-GSG and the M-PC-MI-GSG algorithms become smaller as the data length t increases, which means that the proposed algorithms are effective in parameter estimation for the multivariate autoregressive system.
• Figure 2 and Figure 5 show that the M-PC-MI-GSG algorithm has higher parameter estimation accuracy than the M-PC-GSG algorithm under the same noise variances and same data length. Introducing the innovation length p can effectively improve the parameter estimation accuracy for the M-PC-GSG algorithm, and the parameter estimates can be more stationary as the innovation length p increases.
• It can be seen from Figure 3 and Figure 4 that the M-PC-MI-GSG algorithm can obtain accurate parameter estimates.
• From Figure 6 and Figure 7, we can see that the predicted outputs of the M-PC-MI-GSG algorithm are very close to the true outputs, which indicates that the estimated model can capture the dynamics of the system.

Remark 1.

In order to show the advantages of the identification performances of the proposed algorithm in this paper, the forgetting factor stochastic gradient identification method proposed in [34] is compared with the proposed algorithm. The forgetting factor stochastic gradient identification method is applied to identify the multivariate equation-error autoregressive systems in this paper, and the multivariate forgetting factor generalized stochastic gradient (M-FF-GSG) algorithm is obtained. The M-FF-GSG algorithm is compared with the M-PC-MI-GSG algorithm with $p=5$ through simulation, and the simulation experimental conditions are the same as those in Example 2. The comparison results are shown in Figure 8. It can be seen that the proposed algorithm in this paper has faster identification speed and higher estimation accuracy.

6. Conclusions

Coupling identification concept is an emerging method in the field of system identification in recent decades, one which is usually used in the parameter estimation of the multivariate systems. Its main idea is to utilize the parameter-coupling characteristic, and to identify the parameter of each subsystem model separately and then connected, which can greatly reduce the calculation amounts in the estimation processes. This paper combines the coupling identification concept with the stochastic gradient identification method to propose a new identification algorithm for multivariate equation-error systems. The proposed algorithm has the advantage of a lower computation than the traditional stochastic gradient identification algorithm. Additionally, the multi-innovation identification theory is also a promising identification method, which makes full use of the data information collected in the past to identify the unknown parameters. Based on the partially coupled stochastic gradient algorithm, this paper then introduces the innovation length by applying the multi-innovation identification theory, and proposes the partially coupled multi-innovation stochastic gradient algorithm. The new algorithm also has the advantage of higher parameter estimation accuracy.

The proposed coupling and multi-innovation-based identification methods can be extended to study other multivariate systems with different structures and disturbance noises. Meanwhile, the idea of the algorithms can be utilized when the system identification model has the coupled terms. The future research opportunities are to apply the proposed algorithms to actual engineering production systems to improve the computational efficiency and accuracy of the system identification in production practice. Additionally, the proposed methods in the paper can also combine other mathematical tools and statistical strategies to research the performances of some parameter estimation algorithms for other linear or nonlinear systems with colored noises.