Least square algorithm based on bias compensated principle for parameter estimation of canonical state space model

Abstract

Due to the existence of system noise and unknown state variables, it is difficult to realize unbiased estimation with minimum variance for the parameter estimation of canonical state space model. This paper presents a new least squares estimator based on bias compensation principle to solve this problem, transforms canonical state space into the form suitable for the least square algorithm, introduces an augmented parameter vector and an auxiliary variable, derives parameter estimation formula based on noise compensation, realizes the unbiased estimation, and gives the specific algorithm. A simulation example is provided to verify the effectiveness of the estimator.

Introduction

State space model is an effective tool to describe systems. Due to its simplicity of equations and ease of understanding, state space model has been widely used in modeling and forecasting in economic, control, biological, and other fields. However, the state space model contains both unknown state vectors and unknown parameter vectors, and the nonlinear product relationship between parameter vectors and state vectors, which makes the identification of the model very complicated. ¹ Therefore, the parameter estimation of the state space model has always been a research hotspot in control field. At present, the main methods to estimate the parameters of state space are subspace method, ² Kalman filter method, ^3,4 gradient search method, ^5,6 least squares method, ^7,8 etc.

Literature review

Least squares (LS) estimator is a mathematical optimization technique. The minimum variance of least squares estimator explained by Gauss Markov theorem shows that the effect of ordinary LS estimator is the best one compared with any other linear unbiased estimators. Therefore, the LS method has been widely used. However, the LS method can only use part of observed data, so it depends on the distribution of noise and the statistical properties of missing data, and it is difficult to realize unbiased estimation of parameters. An LS algorithm based on bias compensated principle (BCP) was proposed in this paper ⁹ in view of the limitations of LS algorithm, which achieved the goal that online estimation can be realized when the order and variance of noise are unknown, and this method was widely used in parameter estimation of control systems. ^10–12 Kenji Ikeda et al. analyzed the asymptotic bias of recursive least squares (RLS) estimation under closed-loop conditions, and presented compensation method under the assumption that the noise is white noise. ¹³ Wang et al. proposed a parameter and state estimation algorithm based on bias compensation for observable canonical state space systems with colored noise, which can generate unbiased parameter estimation by the bias correction term. ¹⁴ For Hammerstein autoregressive moving average (ARMAX) systems with autoregressive variables, Zhang et al. proposed an RLS identification method based on the BCP. ¹⁵ Wu et al. combined bias compensation technology with LS estimation algorithm with forgetting factor to estimate parameters of output error model with moving average noise, ¹⁶ then, proposed an LS algorithm based on BCP for parameter estimation of MISO system under white noise. ¹⁷ Mu et al. proposed a globally convergent identification algorithm, modified BCP, and gave a consistent estimation of noise variance rather than introducing auxiliary vectors. ¹⁸ Jia et al. presented a unified framework based on the BCP, which introduces a nonsingular matrix and a noise-independent auxiliary vector to identify the system affected by correlated noise. Due to the rich possibility of selecting matrix and auxiliary vector, the framework has great flexibility. ¹⁹

In the numerous literature studies on system parameter identification, there are not many cases of parameter identification for canonical state space model. By using hierarchical identification principle, Zhuang et al. proposed a new parameter and state estimation algorithm for canonical state space model by replacing noise term with estimation residual in an information vector and achieved good results. ²⁰ However, even if the estimated residuals take the place of the noise item, the other items of algorithm in Ref. 20 still require system noise to be known. In actual systems, the noise is usually random noise, which is very difficult to be observed. When the noise is unknown, it is difficult to estimate parameters unbiased. In this paper, LS based on BCP is used to estimate parameters of canonical state space model, so that model transformation, bias analysis, and bias compensation are carried out to obtain unbiased estimation of system parameters.

Least squares algorithm based on bias compensated principle

Parameter estimation of canonical state space model based on least squares

By LS method, unknown data can be easily obtained, and square sum of errors between obtained data and actual data is the minimum (the least variance). In this section, LS method is applied to estimate parameters of canonical state space model, and state space model in the form of equation (1) is transformed into vector functions suitable for LS method.

Consider this state space equation (1)

$$\{\begin{array}{c}x\left(t+1\right)=Ax\left(t\right)+bu\left(t\right)\\ y\left(t\right)=cx\left(t\right)+v\left(t\right)\end{array}$$

(1)

where, $\mathit{A}=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1\\ -a_n& -a_{n-1}& -a_{n-2}& \cdots & -a_1\end{array}\right]\in {\mathrm{ℝ}}^{n\times n}$ is companion matrix $\mathit{ b}=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]\in {\mathrm{ℝ}}^n$ , $\mathit{ c}=\left[\begin{array}{ccccc}1& 0& 0& \dots & 0\end{array}\right]\in {\mathrm{ℝ}}^{1\times n}$ and $y\left(t\right)\in \mathrm{ℝ}$ are input and output of the system, respectively, $\mathit{x}\left(t\right)\in {\mathrm{ℝ}}^n$ is the unobservable state vector, and $v\left(t\right)\in \mathrm{ℝ}$ is an uncorrelated random white noise with zero mean. Assume that the order $n$ is known, $u\left(t\right)=0,y\left(t\right)=0$ and $v\left(t\right)=0$ when $t\le 0$ . $\mathit{A},\mathit{b}$ is the system parameter matrix and parameter vector to be identified, respectively.

Since $\mathit{A}$ is a companion matrix and $\mathit{c}$ has a standard structure, that is, the first element of $\mathit{c}$ is 1 and the other elements are zeros, it can be seen that the model is observable canonical state space model, that is, parameters of observable canonical model are estimated.

Transform the form of equation (1). The equation of vector $\mathit{Y}\left(t+n\right)$ can be obtained by deducing $\begin{array}{c}y\left(t\right)=cx\left(t\right)+v\left(t\right)\end{array}$ in equations (1) [20]:

$y\left(t+1\right)=cx\left(t+1\right)+v\left(t+1\right)=c\left(\mathit{Ax}\left(t\right)+\mathit{b}u\left(t\right)\right)+v\left(t+1\right)=cAx\left(t\right)+cbu\left(t\right)+v\left(t+1\right),y\left(t+2\right)=cAx\left(t+1\right)+cbu\left(t+1\right)+v\left(t+2\right)=cA\left[\mathit{Ax}\left(t\right)+\mathit{b}u\left(t\right)\right]+cbu\left(t+1\right)+v\left(t+2\right)=c{\mathit{A}}^{2}x\left(t\right)+cAbu\left(t\right)+cbu\left(t+1\right)+v\left(t+2\right),\cdots y\left(t+n-1\right)=c{\mathit{A}}^{n-1}x\left(t\right)+c{\mathit{A}}^{n-2}bu\left(t\right)+c{\mathit{A}}^{n-3}bu\left(t+1\right)+\cdots +cbu\left(t+n-2\right)+v\left(t+n-1\right),$

We have

$$\{\begin{array}{c}y\left(t\right)=cx\left(t\right)+v\left(t\right)\\ y\left(t+1\right)=cAx\left(t\right)+cbu\left(t\right)+v\left(t+1\right)\\ y\left(t+2\right)=c{\mathit{A}}^{2}x\left(t\right)+cAbu\left(t\right)+cbu\left(t+1\right)+v\left(t+2\right)\\ \cdots \\ y\left(t+n-1\right)=c{\mathit{A}}^{n-1}x\left(t\right)+c{\mathit{A}}^{n-2}bu\left(t\right)\\ +c{\mathit{A}}^{n-3}bu\left(t+1\right)+\dots \\ +cbu\left(t+n-2\right)+v\left(t+n-1\right)\\ y\left(t+n\right)=c{\mathit{A}}^{n}x\left(t\right)+c{\mathit{A}}^{n-1}bu\left(t\right)+c{\mathit{A}}^{n-2}bu\left(t+1\right)+\dots \\ +cbu\left(t+n-1\right)+v\left(t+n\right)\end{array}$$

(2)

Define

$$\mathit{Y}\left(t+n\right)=\left[\begin{array}{c}y\left(t\right)\\ y\left(t+1\right)\\ ⋮\\ y\left(t+n-1\right)\end{array}\right]\in {\mathrm{ℝ}}^n$$

(3)

$$\mathit{U}\left(t+n\right)=\left[\begin{array}{c}u\left(t\right)\\ u\left(t+1\right)\\ ⋮\\ u\left(t+n-1\right)\end{array}\right]\in {\mathrm{ℝ}}^n$$

(4)

$$\mathit{V}\left(t+n\right)=\left[\begin{array}{c}v\left(t\right)\\ v\left(t+1\right)\\ ⋮\\ v\left(t+n-1\right)\end{array}\right]\in {\mathrm{ℝ}}^n$$

(5)

$$\mathit{M}=\left[\begin{array}{ccccc}0& 0& \cdots & 0& 0\\ \mathit{cb}& 0& \cdots & 0& 0\\ \mathit{cAb}& \mathit{cb}& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ \mathit{c}{\mathit{A}}^{n-2}\mathit{b}& \mathit{c}{\mathit{A}}^{n-3}\mathit{b}& \cdots & \mathit{cb}& 0\end{array}\right]\in {\mathrm{ℝ}}^{n\times n}$$

(6)

$${\mathit{Q}}_0=\left[\begin{array}{c}\mathit{c}\\ \mathit{cA}\\ ⋮\\ \mathit{c}{\mathit{A}}^{n-1}\end{array}\right]\in {\mathrm{ℝ}}^{n\times n}$$

(7)

Due to the special structures of $\mathit{c}$ and $\mathit{A}$ , we have ${\mathit{Q}}_0={\mathit{I}}_n$ , which is an unit matrix. The first $n$ equations of equation (2) can be rewritten as (8)

$$\begin{array}{c}Y\left(t+n\right)={\mathit{Q}}_{0}x\left(t\right)+MU\left(t+n\right)+V\left(t+n\right)\\ =x\left(t\right)+MU\left(t+n\right)+V\left(t+n\right)\end{array}$$

(8)

Transform (8) to get (9)

$$\mathit{x}\left(t\right)=\mathit{Y}\left(t+n\right)-\mathit{MU}\left(t+n\right)-\mathit{V}\left(t+n\right)$$

(9)

Substituting (9) into equation (2) gives (10)

$$y\left(t+n\right)=\mathit{c}{\mathit{A}}^n\mathit{x}\left(t\right)+\mathit{c}{\mathit{A}}^{n-1}\mathit{b}u\left(t\right)+\mathit{c}{\mathit{A}}^{n-2}\mathit{b}u\left(t+1\right)+\cdots +\mathit{cb}u\left(t+n-1\right)+v\left(t+n\right)=\mathit{c}{\mathit{A}}^n\left[\mathit{Y}\left(t+n\right)-\mathit{MU}\left(t+n\right)-\mathit{V}\left(t+n\right)\right]+\mathit{c}{\mathit{A}}^{n-1}\mathit{b}u\left(t\right)+\mathit{c}{\mathit{A}}^{n-2}\mathit{b}u\left(t+1\right)+\cdots +\mathit{cb}u\left(t+n-1\right)+v\left(t+n\right)=\mathit{c}{\mathit{A}}^n\left[\mathit{Y}\left(t+n\right)-\mathit{V}\left(t+n\right)\right]-\mathit{c}{\mathit{A}}^n\mathit{MU}\left(t+n\right)+\mathit{c}{\mathit{A}}^{n-1}\mathit{b}u\left(t\right)+\mathit{c}{\mathit{A}}^{n-2}\mathit{b}u\left(t+1\right)+\cdots +\mathit{cb}u\left(t+n-1\right)+v\left(t+n\right)=\mathit{c}{\mathit{A}}^n\left[\mathit{Y}\left(t+n\right)-\mathit{V}\left(t+n\right)\right]-\mathit{c}{\mathit{A}}^n\mathit{MU}\left(t+n\right)+\left[\mathit{c}{\mathit{A}}^{n-1}\mathit{b},\mathit{c}{\mathit{A}}^{n-2}\mathit{b},\cdots ,\mathit{cb}\right]\left[\begin{array}{c}u\left(t\right)\\ u\left(t+1\right)\\ ⋮\\ u\left(t+n-1\right)\end{array}\right]+v\left(t+n\right)=\left[{\mathit{Y}}^T\left(t+n\right)-{\mathit{V}}^T\left(t+n\right)\right]{\left(\mathit{c}{\mathit{A}}^n\right)}^T-{\mathit{U}}^T\left(t+n\right){\left(\mathit{c}{\mathit{A}}^n\mathit{M}\right)}^T+\left[\mathit{c}{\mathit{A}}^{n-1}\mathit{b},\mathit{c}{\mathit{A}}^{n-2}\mathit{b},\cdots ,\mathit{cb}\right]\left[\begin{array}{c}u\left(t\right)\\ u\left(t+1\right)\\ ⋮\\ u\left(t+n-1\right)\end{array}\right]+v\left(t+n\right)=\left[{\mathit{Y}}^T\left(t+n\right)-{\mathit{V}}^T\left(t+n\right)\right]{\left(\mathit{c}{\mathit{A}}^n\right)}^T-{\mathit{U}}^T\left(t+n\right){\left(\mathit{c}{\mathit{A}}^n\mathit{M}\right)}^T+{\mathit{U}}^T\left(t+n\right)\left[\begin{array}{c}\mathit{c}{\mathit{A}}^{n-1}\mathit{b}\\ ⋮\\ \mathit{cb}\end{array}\right]+v\left(t+n\right)=\left[{\mathit{Y}}^T\left(t+n\right),{\mathit{U}}^T\left(t+n\right)\right]\left[\begin{array}{c}{\left(\mathit{c}{\mathit{A}}^n\right)}^T\\ {\left(-\mathit{c}{\mathit{A}}^n\mathit{M}+\left[\mathit{c}{\mathit{A}}^{n-1}\mathit{b},\cdots ,\mathit{cb}\right]\right)}^T\end{array}\right]-{\mathit{V}}^T\left(t+n\right){\left(\mathit{c}{\mathit{A}}^n\right)}^T+v\left(t+n\right).$$

(10)

Define

$$\mathit{\varphi }\left(t+n\right)=\left[\begin{array}{c}\mathit{Y}\left(t+n\right)\\ \mathit{U}\left(t+n\right)\end{array}\right]\in {\mathrm{ℝ}}^{2n\times 1}$$

(11)

$${\mathit{\theta }}_1={\left[\mathit{c}{\mathit{A}}^n\right]}^T$$

(12)

$${\mathit{\theta }}_2={\left[-\mathit{c}{\mathit{A}}^n\mathit{M}+\left[\mathit{c}{\mathit{A}}^{n-1}\mathit{b},\cdots ,\mathit{cb}\right]\right]}^T$$

(13)

$$\tilde{v}\left(t+n\right)=v\left(t+n\right)-{\mathit{V}}^T\left(t+n\right){\left(\mathit{c}{\mathit{A}}^n\right)}^T$$

(14)

Obtain

$${\mathit{\varphi }}^T\left(t+n\right)=\left[\begin{array}{c}{\mathit{Y}}^T\left(t+n\right),{\mathit{U}}^T\left(t+n\right)\end{array}\right]$$

(15)

$$\tilde{\mathit{V}}\left(t+n\right)=\left[\begin{array}{c}\tilde{v}\left(t\right)\\ \tilde{v}\left(t+1\right)\\ ⋮\\ \tilde{v}\left(t+n-1\right)\end{array}\right]\in {\mathrm{ℝ}}^n$$

(16)

So, equation (10) can be written as equation (17)

$$y\left(t+n\right)={\mathit{\varphi }}^T\left(t+n\right)\left[\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\end{array}\right]+\tilde{v}\left(t+n\right)={\mathit{\varphi }}^T\left(t+n\right)\mathit{\theta }+\tilde{v}\left(t+n\right)$$

(17)

Replacing $t$ in equation (17) with $t-n$ yields equation (18)

$$y\left(t\right)={\mathit{\varphi }}^T\left(t\right)\mathit{\theta }+\tilde{v}\left(t\right)$$

(18)

Define

$$\mathbf{\Phi }\left(t+n\right)=\left[\begin{array}{c}{\mathit{\varphi }}^T\left(t\right)\\ \cdots \\ {\mathit{\varphi }}^T\left(t+n-1\right)\end{array}\right]$$

(19)

So, the vector function form suitable for LS method is obtained

$$\mathit{Y}\left(t+n\right)=\mathbf{\Phi }\left(t+n\right)\mathit{\theta }+\tilde{\mathit{V}}\left(t+n\right)$$

(20)

Replacing $t$ in equation (20) with $t-n$ yields equation (21)

$$\mathit{Y}\left(t\right)=\mathbf{\Phi }\left(t\right)\mathit{\theta }+\tilde{\mathit{V}}\left(t\right)$$

(21)

According to the principle of LS, the loss function can be defined as (22)

$$J\left(\mathit{\theta }\right)={\displaystyle \sum }_{i=1}^t{\tilde{v}}^2\left(i\right)={\displaystyle \sum }_{i-1}^t{\left[y\left(i\right)-{\mathit{\varphi }}^T\left(i\right)\mathit{\theta }\right]}^2$$

(22)

LS estimation can be obtained by solving the minimization of loss function

$${\widehat{\mathit{\theta }}}_{\mathit{LS}}={\left({\mathbf{\Phi }}^{\text{T}}\left(t\right)\mathbf{\Phi }\left(t\right)\right)}^{-1}{\mathbf{\Phi }}^{\text{T}}\left(t\right)\mathit{Y}\left(t\right)$$

(23)

If the above equations are used to estimate ${\widehat{\mathit{\theta }}}_{\mathit{LS}}$ directly, and then $\mathit{A}$ and $\mathit{b}$ are deduced, parameter estimation of equation (1) can be achieved. However, it can be seen from equation (14) that noise $\tilde{v}\left(t+n\right)$ may no longer be random white noise with zero mean, so that estimation of parameters by LS method may not be the least variance, that is, unbiased estimation of parameters may not be possible. Therefore, bias compensation method will be used to eliminate the bias; the above LS estimator will be modified in order to achieve unbiased estimation.

Least squares estimator based on bias compensated principle

In this section, equation (21) is continuously to be deduced by LS gives (24)

$${\widehat{\mathit{\theta }}}_{\mathit{LS}}\left(n\right)={\left(\frac{1}{n}{\mathbf{\Phi }}^{\text{T}}\left(t\right)\mathbf{\Phi }\left(t\right)\right)}^{-1}\frac{1}{n}{\mathbf{\Phi }}^{\text{T}}\left(t\right)\mathit{Y}\left(t\right)={\left(\frac{1}{n}\left(\mathit{\varphi }\left(t-n\right)\cdots \mathit{\varphi }\left(t-1\right)\right)\left(\begin{array}{c}{\mathit{\varphi }}^T\left(t-n\right)\\ ⋮\\ {\mathit{\varphi }}^T\left(t-1\right)\end{array}\right)\right)}^{-1}\left(\frac{1}{n}\left(\mathit{\varphi }\left(t-n\right)\cdots \mathit{\varphi }\left(t-1\right)\right)\left(\begin{array}{c}y\left(t-n\right)\\ ⋮\\ y\left(t-1\right)\end{array}\right)\right)={\left(\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}\left(\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)y\left(i\right)\right)\right)={\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}^{-1}\left(n\right){\widehat{\mathit{R}}}_{\mathit{\varphi y}}\left(n\right).$$

(24)

where ${\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}^{-1}\left(n\right)={\left(\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}$ , ${\widehat{\mathit{R}}}_{\mathit{\varphi y}}\left(n\right)=\left(\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)y\left(i\right)\right)\right)$ .

Substituting $y\left(t\right)={\mathit{\varphi }}^T\left(i\right)\mathit{\theta }+\tilde{v}\left(t\right)$ into (24) gives

$${\widehat{\mathit{\theta }}}_{\mathit{LS}}\left(n\right)={\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)\left({\mathit{\varphi }}^T\left(i\right)\mathit{\theta }+\tilde{v}\left(t\right)\right)\right)\right)={\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}\left(\mathit{\theta }\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)+\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)\tilde{v}\left(t\right)\right)\right)=\mathit{\theta }{\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)+{\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)\tilde{v}\left(t\right)\right)\right)=\mathit{\theta }+{\left(\left(\frac{1}{n}\right){\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)\right)}^{-1}\left(\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)\tilde{v}\left(t\right)\right)\right)=\mathit{\theta }+{\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}^{-1}\left(n\right){\widehat{\mathit{r}}}_{\mathit{\varphi v}}\left(n\right).$$

(25)

Taking probability limit of (25) yields

$$pli{m}_{n\to \infty }{\widehat{\mathit{\theta }}}_{\mathit{LS}}\left(n\right)=\mathit{\theta }+pli{m}_{n\to \infty }{\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}^{-1}\left(n\right){\widehat{\mathit{r}}}_{\mathit{\varphi v}}\left(n\right)=\mathit{\theta }+{\mathit{R}}_{\mathit{\varphi \varphi }}^{-1}{\mathit{r}}_{\mathit{\varphi v}}$$

(26)

where

$${\mathit{R}}_{\mathit{\varphi \varphi }}=pli{m}_{n\to \infty }{\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}\left(n\right)=pli{m}_{n\to \infty }{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right){\mathit{\varphi }}^T\left(i\right)\right)$$

(27)

$${\mathit{r}}_{\mathit{\varphi v}}=pli{m}_{n\to \infty }{\widehat{\mathit{r}}}_{\mathit{\varphi v}}\left(n\right)=pli{m}_{n\to \infty }{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\varphi }\left(i\right)\tilde{v}\left(i\right)\right)$$

(28)

It can be seen from equation (26) that LS estimator is biased, and bias term is ${\mathit{R}}_{\mathit{\varphi \varphi }}^{-1}{\mathit{r}}_{\mathit{\varphi v}}\mathit{ }$ caused by noise. In order to eliminate the influence of bias term, equation (29) is obtained

$$\widehat{\mathit{\theta }}\left(n\right)={\widehat{\mathit{\theta }}}_{\mathit{LS}}\left(n\right)-{\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}^{-1}\left(n\right){\widehat{\mathit{r}}}_{\mathit{\varphi v}}\left(n\right)$$

(29)

In order to obtain an unbiased estimate of $\widehat{\mathit{\theta }}$ , a constant is needed to replace ${\widehat{\mathit{r}}}_{\mathit{\varphi v}}\left(n\right)$ to compensate for the influence of bias. Assuming that system operates in an open-loop mode, it means that inputs $u\left(i\right)$ and $v\left(i\right)$ are not statistically relevant, then equation (28) can be rewritten as equation (30)

$${\mathit{r}}_{\mathit{\varphi v}}=pli{m}_{n\to \infty }\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\begin{array}{c}\mathit{Y}\left(i\right)\\ \mathit{U}\left(i\right)\end{array}\right)\tilde{v}\left(i\right)=pli{m}_{n\to \infty }\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\begin{array}{c}\mathit{Y}\left(i\right)\tilde{v}\left(i\right)\\ \mathit{U}\left(i\right)\tilde{v}\left(i\right)\end{array}\right)=\left[\begin{array}{c}{\mathit{I}}_n\\ 0_n\end{array}\right]pli{m}_{n\to \infty }\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{Y}\left(i\right)\tilde{v}\left(i\right)\right)=\left[\begin{array}{c}{\mathit{I}}_n\\ 0_n\end{array}\right]{\mathit{r}}_{\mathit{yv}}$$

(30)

where, ${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)=\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{Y}\left(i\right)\tilde{v}\left(i\right)\right),$ $\mathit{ }{\mathit{r}}_{\mathit{yv}}= pli{m}_{n\to \infty }\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{Y}\left(i\right)\tilde{v}\left(i\right)\right).\mathit{ }$

Equation (29) can be written as equation (31)

$$\widehat{\mathit{\theta }}\left(n\right)={\widehat{\mathit{\theta }}}_{\mathit{LS}}\left(n\right)-{\widehat{\mathit{R}}}_{\mathit{\varphi \varphi }}^{-1}\left(n\right)\left[\begin{array}{c}{\mathit{I}}_n\\ 0_n\end{array}\right]{\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)$$

(31)

It can be seen from (31) that the key is to find ${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right),$ so augmented regression vector $\mathit{\psi }\left(i\right)$ and augmented parameter vector $\overline{\mathit{\theta }}$ are introduced to satisfy ${\mathit{\psi }}^T\left(i\right)\overline{\mathit{\theta }}={\widehat{\mathit{\varphi }}}^T\left(i\right)\mathit{\theta }$ , which is commonly used in BCP. It is worth noting that no matter how different the structures of $\overline{\mathit{\theta }}$ and $\mathit{\psi }\left(i\right)$ are, the calculation of ${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)$ is similar. ¹⁹ Therefore, the method described below is universal.

Introduce an auxiliary vector $\mathit{\xi }\left(i\right)\in {\mathrm{ℝ}}^m(m\ge n)$ which satisfies equation (32) and a nonsingular matrix $\mathbf{ ℱ}\in {\mathrm{ℝ}}^{\left(2n+m\right)\times \left(2n+m\right)}$

$$pli{m}_{n\to \infty }\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\left(\mathit{\xi }\left(i\right)\tilde{v}\left(i\right)\right)=0$$

(32)

Construct an augmented regression vector $\mathit{ \psi }\left(i\right)$ and an augmented parameter vector $\overline{\mathit{\theta }}$ by introducing auxiliary vector $\mathit{\xi }\left(i\right)$ and nonsingular matrix $\mathbf{ ℱ}$ , this step is the most important step, the selection of $\mathit{\psi }\left(i\right)$ and $\overline{\mathit{\theta }}$ makes calculation of ${\widehat{\mathit{r}}}_{\mathit{yv}}$ becomes possible

$$\mathit{\psi }\left(i\right)=\mathbf{ℱ}\left[\begin{array}{c}\mathit{\varphi }\left(i\right)\\ \mathit{\xi }\left(i\right)\end{array}\right]\in {\mathrm{ℝ}}^{\left(2n+m\right)}$$

(33)

$$\overline{\mathit{\theta }}={\left({\mathbf{ℱ}}^{-1}\right)}^T\left[\begin{array}{c}\mathit{\theta }\\ 0\end{array}\right]\in {\mathrm{ℝ}}^{\left(2n+m\right)}$$

(34)

The $\overline{\mathit{\theta }}$ and $\mathit{\psi }\left(\mathit{i}\right)$ satisfy equation (35)

$${\mathit{\psi }}^T\left(i\right)\overline{\mathit{\theta }}={\mathit{\varphi }}^T\left(i\right)\mathit{\theta }$$

(35)

Rewrite $y\left(t\right)={\mathit{\varphi }}^T\left(i\right)\mathit{\theta }+\tilde{v}\left(t\right)$ as $y\left(t\right)={\mathit{\psi }}^T\left(t\right)\overline{\mathit{\theta }}+\tilde{v}\left(t\right).$ It is easy to yield equation (36) from equation (34)

$${\mathbf{ℱ}}^T\overline{\mathit{\theta }}=\left[\begin{array}{c}\mathit{\theta }\\ 0\end{array}\right]$$

(36)

Let

$$\mathbf{ℱ}=\begin{array}{c}\left[{\mathbf{ℱ}}_1\mathit{ }{\mathbf{ℱ}}_2\mathit{ Z}\right]\end{array}=\begin{array}{c}\left[{\mathbf{ℱ}}_0\mathit{ Z}\right]\end{array}$$

(37)

where ${\mathbf{ℱ}}_1,{\mathbf{ℱ}}_2\in {\mathrm{ℝ}}^{\left(2n+m\right)\times n},{\mathbf{ℱ}}_0\in {\mathrm{ℝ}}^{\left(2n+m\right)\times 2n},\mathit{Z}\in {\mathrm{ℝ}}^{\left(2n+m\right)\times m}.$

Equation (38) can be obtained

$${\mathit{Z}}^T\overline{\mathit{\theta }}=0$$

(38)

Similar to equations (24) and (26), equations (39) and (40) can be obtained in the same way

$${\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}\left(n\right)={\widehat{\mathit{R}}}_{\mathit{\psi \psi }}^{-1}\left(n\right){\widehat{\mathit{r}}}_{\mathit{\psi y}}\left(n\right)$$

(39)

$$pli{m}_{n\to \infty }{\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}\left(n\right)=\overline{\mathit{\theta }}+\left(pli{m}_{n\to \infty }{\widehat{\mathit{R}}}_{\mathit{\psi \psi }}^{-1}\right){\mathbf{ℱ}}_1{\mathit{r}}_{\mathit{yv}}$$

(40)

where

$${\widehat{\mathit{R}}}_{\mathit{\psi \psi }}\left(n\right)=\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\mathit{\psi }\left(i\right){\mathit{\psi }}^T\left(i\right)$$

(41)

$${\widehat{\mathit{r}}}_{\mathit{\psi y}}\left(n\right)=\frac{1}{n}{\displaystyle \sum }_{i=t-n}^{t-1}\mathit{\psi }\left(i\right)y\left(i\right)$$

(42)

Pre-multiplying ${\mathit{Z}}^T$ on both sides of equation (40) and using (38) yield

$$\begin{array}{l}{\mathit{Z}}^{T}pli{m}_{n\to \infty }{\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}\left(n\right)={\mathit{Z}}^T\overline{\mathit{\theta }}+{\mathit{Z}}^T\left(pli{m}_{n\to \infty }{\widehat{\mathit{R}}}_{\mathit{\psi \psi }}^{-1}\left(n\right)\right){\mathbf{ℱ}}_1{\mathit{r}}_{\mathit{yv}}\\ ={\mathit{Z}}^{T}pli{m}_{n\to \infty }{\widehat{\mathit{R}}}_{\mathit{\psi \psi }}^{-1}\left(n\right){\mathbf{ℱ}}_1{\mathit{r}}_{\mathit{yv}}\end{array}$$

(43)

Define

$$\widehat{\mathit{D}}={\mathit{Z}}^T{\widehat{\mathit{R}}}_{\mathit{\psi \psi }}^{-1}\left(n\right){\mathbf{ℱ}}_1\in {\mathrm{ℝ}}^{m\times n}$$

(44)

Thus, the estimation of ${\mathit{r}}_{\mathit{yv}}\mathit{ }$ can be determined by the weighted LS as equation (45)

$${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)={\left({\widehat{\mathit{D}}}^T{\mathit{W}}_{\mathit{LS}}\widehat{\mathit{D}}\right)}^{-1}{\widehat{\mathit{D}}}^T{\mathit{W}}_{\mathit{LS}}{\mathit{Z}}^T{\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}\left(n\right)$$

(45)

where ${\mathit{W}}_{\mathit{LS}}$ is a positive definite weighting matrix. The choice of ${\mathit{W}}_{\mathit{LS}}$ does not affect the consistency properties of ${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)$ , but may have a considerable impact on the statistical accuracy of ${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)$ . One simple choice is ${\mathit{W}}_{\mathit{LS}}=\mathit{I}$ . ¹⁹ In the case of $m=n$ , if $\widehat{\mathit{D}}$ is non-singular, ${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)$ simply becomes as equation (46)

$${\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)={\widehat{\mathit{D}}}^{-1}{\mathit{Z}}^T{\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}\left(n\right)$$

(46)

Similar to equations (24) and (26), equation (47) can be obtained from (40) in the same way

$$\widehat{\overline{\mathit{\theta }}}\left(n\right)={\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}\left(n\right)-{\widehat{\mathit{R}}}_{\mathit{\psi \psi }}^{-1}\left(n\right){\mathbf{ℱ}}_1{\widehat{\mathit{r}}}_{\mathit{yv}}\left(n\right)$$

(47)

By substituting (45) into (47), consistent estimation of the augmented parameter vector $\overline{\mathit{\theta }}$ can be obtained, finally, the global consistent estimation of original parameter $\mathit{\theta }$ can be obtained by equation (36) and rewritten as equation (48)

$$\widehat{\mathit{\theta }}\left(n\right)={\mathbf{ℱ}}_0^T\widehat{\overline{\mathit{\theta }}}\left(n\right)$$

(48)

Above is the unified BCP method. Through choosing appropriately $\mathit{\xi }\left(i\right)$ and $\mathbf{ℱ}$ , various BCP-based methods can be obtained to achieve a global consistent estimation of parameters.

Parameter estimation of canonical state space model based on BCP-LS

For canonical state space model, this section uses the method introduced in 2.2 to derive a specific parameter estimation algorithm.

Choose

$$\mathbf{ℱ}={\mathit{I}}_{3n},\mathit{ \xi }\left(i\right)=\underset{¯}{\mathit{U}}\left(i\right)\in {\mathrm{ℝ}}^n$$

(49)

where $\underset{¯}{\mathit{U}}\left(i\right)=\left[\begin{array}{c}u\left(t-2n\right)\\ u\left(t-2n+1\right)\\ ⋮\\ u\left(t-n-1\right)\end{array}\right]$ ,weobtain $\overline{\mathit{\theta }}=\left[\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\\ 0\end{array}\right]\in {\mathrm{ℝ}}^{3n},\mathit{\psi }\left(i\right)=\left(\begin{array}{c}\mathit{Y}\left(i\right)\\ \mathit{U}\left(i\right)\\ \underset{¯}{\mathit{U}}\left(i\right)\end{array}\mathit{ }\right]\in {\mathrm{ℝ}}^{3n}.$

Then, ${\widehat{\mathit{r}}}_{\mathit{yv}}$ can be obtained according to equation (45), and the global unbiased estimate of the parameter $\mathit{\theta }$ can be obtained by equations (47) and (48). Furthermore, parameters of the system can be estimated by input and output data. Since the relationship between parameter vector $\mathit{\theta }$ and parameters $\mathit{A},\mathit{b},\mathit{c}$ has been established in parameter estimation of canonical state space model based on least squares, the algorithm can be directly derived from algebraic operations.

Post-multiplying ${\mathit{Q}}_0=\left[\begin{array}{c}\mathit{c}\\ \mathit{cA}\\ ⋮\\ \mathit{c}{\mathit{A}}^{n-1}\end{array}\right]={\mathit{I}}_n$ by $\mathit{b}$ gives $\left[\begin{array}{c}\mathit{cb}\\ \mathit{cAb}\\ ⋮\\ \mathit{c}{\mathit{A}}^{n-1}\mathit{b}\end{array}\right]=\mathit{b}.$ Observing the structure of vector $\mathit{b}$ to get equation (50)

$$\mathit{c}{\mathit{A}}^{k-1}\mathit{b}=b_k,\mathit{ }k=\mathrm{1,2},\dots ,n$$

(50)

The matrix $\mathit{M}$ in parameter estimation of canonical state space model based on least squares can be simplified by (50) as

$$\mathit{M}=\left[\begin{array}{ccccc}0& 0& \cdots & 0& 0\\ b_1& 0& \cdots & 0& 0\\ b_2& b_1& \ddots & ⋮& ⋮\\ ⋮& ⋮ & \ddots & 0& 0\\ b_{n-1}& b_{n-2}& \cdots & b_1& 0\end{array}\right]$$

(51)

Post-multiplying $\left[\begin{array}{c}\mathit{c}\\ \mathit{cA}\\ ⋮\\ \mathit{c}{\mathit{A}}^{n-1}\end{array}\right]={\mathit{I}}_n$ by $\mathit{A}$ gives

$$\left[\begin{array}{c}\mathit{cA}\\ \mathit{c}{\mathit{A}}^2\\ ⋮\\ \mathit{c}{\mathit{A}}^n\end{array}\right]=\mathit{A}$$

(52)

Post-multiplying $\mathit{c}{\mathit{A}}^n=\begin{array}{c}[-a_n,-a_{n-1},\cdots ,-a_1]\end{array}$ by $\mathit{M}$ gives

$$\mathit{c}{\mathit{A}}^n\mathit{M}={\left[\begin{array}{c}-\left(b_1{a}_{n-1}+b_2{a}_{n-2}+\cdots +b_{n-1}{a}_1\right)\\ -\left(b_1{a}_{n-2}+b_2{a}_{n-3}+\cdots +b_{n-2}{a}_1\right)\\ ⋮\\ -b_1{a}_1\\ 0\end{array}\right]}^T$$

(53)

Equation (54) can be obtained according to equations (52) and (12)

$${\mathit{\theta }}_1={\left[\mathit{c}{\mathit{A}}^n\right]}^T=\left[\begin{array}{c}-a_n\\ -a_{n-1}\\ ⋮\\ -a_1\end{array}\right]$$

(54)

Equation (55) can be obtained according to equations (53) and (13)

$$\begin{array}{l}{\mathit{\theta }}_2={\left[-\mathit{c}{\mathit{A}}^n\mathit{M}+\left[\mathit{c}{\mathit{A}}^{n-1}\mathit{b},\cdots ,\mathit{cb}\right]\right]}^T\\ =\left[\begin{array}{c}{b}_1{a}_{n-1}+b_2{a}_{n-2}+\cdots +b_{n-1}{a}_1+b_n\\ b_1{a}_{n-2}+b_2{a}_{n-3}+\cdots +b_{n-2}{a}_1+b_{n-1}\\ ⋮\\ b_1{a}_1+b_2\\ b_1\end{array}\right]\\ =\left[\begin{array}{ccccc}{a}_{n-1}& a_{n-2}& \cdots & a_1& 1\\ a_{n-2}& a_{n-3}& \cdots & 1& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_1& 1& \cdots & 0& 0\\ 1& 0& \cdots & 0& 0\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_{n-1}\\ b_n\end{array}\right]\end{array}$$

(55)

Since the estimated value $\widehat{\mathit{\theta }}\left(n\right)$ of $\mathit{\theta }=\left[\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\end{array}\right]$ has been obtained previously, the estimated value of $a_i$ has been obtained, that is, the estimated value $\widehat{\mathit{A}}$ of $\mathit{A}$ can be obtained. Then equation (56) can be obtained due to the expression of ${\mathit{\theta }}_2$

$$\widehat{\mathit{b}}=\left[\begin{array}{c}{\widehat{b}}_1\\ {\widehat{b}}_2\\ ⋮\\ {\widehat{b}}_{n-1}\\ {\widehat{b}}_n\end{array}\right]={\left[\begin{array}{ccccc}{\widehat{a}}_{n-1}& {\widehat{a}}_{n-2}& \cdots & {\widehat{a}}_1& 1\\ {\widehat{a}}_{n-2}& {\widehat{a}}_{n-3}& \cdots & 1& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\widehat{a}}_1& 1& \cdots & 0& 0\\ 1& 0& \cdots & 0& 0\end{array}\right]}^{-1}{\widehat{\mathit{\theta }}}_2$$

(56)

A consistent estimated value $\widehat{\mathit{b}}$ of $\mathit{b}$ can be obtained from (56). Now we have realized the global consistent unbiased estimation of parameters of canonical state space model (1).

The specific algorithm is summarized as follows:

1) Transform equation (1) into the form of $\mathit{Y}\left(t\right)=\mathbf{\Phi }\left(t\right)\mathit{\theta }+\tilde{\mathit{V}}\left(t\right)$ according to equations (2)–(21);

2) collect the input/output data of $\left\{u\left(t\right),y\left(t\right)\right\}$ ;

3) introduce auxiliary vector $\mathit{\xi }$ and matrix $\mathbf{ℱ}$ , then construct $\mathit{\psi }$ and $\overline{\mathit{\theta }}$ according to equations (33) and (34);

4) calculate ${\widehat{\mathit{R}}}_{\mathit{\psi \psi }}$ by equation (41);

5) calculate ${\widehat{\mathit{r}}}_{\mathit{\psi y}}$ by equation (42);

6) calculate ${\widehat{\overline{\mathit{\theta }}}}_{\mathit{LS}}$ by equation (39);

7) calculate $\widehat{\mathit{D}}$ by equation (44);

8) define ${\mathit{W}}_{\mathit{LS}}=\mathit{I},$ calculate ${\widehat{\mathit{r}}}_{\mathit{yv}}\mathit{\hspace{1em}}$ by equation (45);

9) calculate $\widehat{\overline{\mathit{\theta }}}\left(n\right)$ by equation (47);

10) calculate $\widehat{\mathit{\theta }}\left(n\right)$ by equation (48), then the global consistent least squares estimation of $\mathit{\theta }$ is obtained;

11) calculate $\widehat{\mathit{A}},\widehat{\mathit{b}}$ by equations (54) and (56), algorithm is end.

Simulation research

In this section, several examples of a canonical state space model are used to verify performance of presented estimator.

Example 1

For equation (57), $\left\{u\left(t\right)\right\}$ is taken as an independent random signal sequence with zero mean and unit variance, and $\left\{v\left(t\right)\right\}$ as a white noise sequence with zero mean and variances ${\sigma }_1^2={0.1}^2$ and ${\sigma }_2^2={1.0}^2$ , respectively. Time is represented as $t$ and parameter estimation error is calculated by $\epsilon =\widehat{\mathit{\theta }}-{\mathit{\theta }}_2$ , which is similar to the example in Ref. 20

$$\{\begin{array}{c}x\left(t+1\right)=\left[\begin{array}{cc}0& 1\\ -0.8& 1.6\end{array}\right]x\left(t\right)+\left[\begin{array}{c}1\\ -1\end{array}\right]u\left(t\right)\\ y\left(t\right)=\begin{array}{c}[1 0\end{array}]x\left(t\right)+v\left(t\right). \end{array}$$

(57)

First, true value of $\mathit{\theta }$ is calculated:

$$\mathit{\theta }=\left[\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\end{array}\right]=\left[\begin{array}{c}-0.8\\ 1.6\\ -2.6\\ 1\end{array}\right],$$

Let

$$\left[\begin{array}{c}{\theta }_{11}\\ {\theta }_{12}\\ {\theta }_{21}\\ {\theta }_{22}\end{array}\right]=\left[\begin{array}{c}-0.8\\ 1.6\\ -2.6\\ 1\end{array}\right].$$

Then algorithm in 2.3 is used to estimate parameters. When $t$ = 3500, the estimated value and error of $\mathit{\theta }$ are shown in Table 1:

Table 1.

Estimated value and error of parameter $\mathit{\theta }$ .

	${\mathit{\theta }}_{11}$	${\mathit{\theta }}_{12}$	${\mathit{\theta }}_{21}$	${\mathit{\theta }}_{22}$	$\mathit{\epsilon }$
True value	−0.80000	1.60000	−2.60000	1.00000
Estimated value ( ${\mathit{\sigma }}_1$ )	−0.80039	1.59986	−2.60101	1.00049	0.0011966
Estimated value ( ${\mathit{\sigma }}_2$ )	−0.80079	1.60029	−2.60211	1.00149	0.0027167

After obtaining estimated value of $\mathit{\theta }$ , the estimated values of parameters $\mathit{A}$ and $\mathit{b}$ are obtained, as shown in Table 2.

Table 2.

Estimated values and true values of parameters $\mathit{A}$ and $\mathit{b}$ .

	$\mathit{A}$	$\mathit{b}$
True value	$\left[\begin{array}{cc}0& 1\\ -0.80000& 1.60000\end{array}\right]$	$\left[\begin{array}{c}1\\ -1\end{array}\right]$
Estimated value ( ${\mathit{\sigma }}_1$ )	$\left[\begin{array}{cc}0& 1\\ -0.80039& 1.59986\end{array}\right]$	$\left[\begin{array}{c}1.00049\\ -1.00037\end{array}\right]$
Estimated value ( ${\mathit{\sigma }}_2$ )	$\left[\begin{array}{cc}0& 1\\ -0.80079& 1.60029\end{array}\right]$	$\left[\begin{array}{c}1.00149\\ -0.99944\end{array}\right]$

The curve of parameter $\mathit{\theta }$ estimation error $\epsilon$ with time $t$ is shown in Figure 1.

Figure 1.

Parameter estimation error $\epsilon$ with time $t$ .

Example 2

For equations (58), in the same way, $\left\{u\left(t\right)\right\}$ is taken as an independent random signal sequence with zero mean and unit variance, and $\left\{v\left(t\right)\right\}$ as a white noise sequence with zero mean and variances ${\sigma }_1^2={0.1}^2$ and ${\sigma }_2^2={1.0}^2$ , respectively. Time is represented as $t$ , and parameter estimation error is calculated by $\epsilon =\widehat{\mathit{\theta }}-{\mathit{\theta }}_2.$

$$\{\begin{array}{c}x\left(t+1\right)=\left[\begin{array}{cc}0& 1\\ -1.2& 0.6\end{array}\right]x\left(t\right)+\left[\begin{array}{c}1\\ -1\end{array}\right]u\left(t\right)\\ y\left(t\right)=\begin{array}{c}[1 0\end{array}]x\left(t\right)+v\left(t\right) \end{array}$$

(58)

First, true value of $\mathit{\theta }$ is calculated.

$\mathit{\theta }=\left[\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\end{array}\right]=\left[\begin{array}{c}-1.2\\ 0.6\\ -1.6\\ 1\end{array}\right],$

Let

$\left[\begin{array}{c}{\theta }_{11}\\ {\theta }_{12}\\ {\theta }_{21}\\ {\theta }_{22}\end{array}\right]=\left[\begin{array}{c}-1.2\\ 0.6\\ -1.6\\ 1\end{array}\right]$ .

Then algorithm in 2.3 is used to estimate parameters. When $t$ = 3500, the estimated value and error of $\mathit{\theta }$ are shown in Table 3:

Table 3.

Estimated value and error of parameter $\mathit{\theta }$ .

	${\mathit{\theta }}_{11}$	${\mathit{\theta }}_{12}$	${\mathit{\theta }}_{21}$	${\mathit{\theta }}_{22}$	$\mathit{\epsilon }$
True value	−1.20000	0.60000	−1.60000	1.00000
Estimated value ( ${\mathit{\sigma }}_1$ )	−1.20041	0.59787	−1.60103	1.00037	0.0024296
Estimated value ( ${\mathit{\sigma }}_2$ )	−1.20088	0.60129	−1.60215	1.00109	0.0028721

After obtaining estimated value of $\mathit{\theta }$ , the estimated values of parameters $\mathit{A}$ and $\mathit{b}$ are obtained, as shown in Table 4.

Table 4.

Estimated values and true values of parameters $\mathit{A}$ and $\mathit{b}$ .

	$\mathit{A}$	$\mathit{b}$
True value	$\left[\begin{array}{cc}0& 1\\ -1.20000& 0.60000\end{array}\right]$	$\left[\begin{array}{c}1\\ -1\end{array}\right]$
Estimated value ( ${\mathit{\sigma }}_1$ )	$\left[\begin{array}{cc}0& 1\\ -1.20041& 0.59787\end{array}\right]$	$\left[\begin{array}{c}1.00042\\ -1.00031\end{array}\right]$
Estimated value ( ${\mathit{\sigma }}_2$ )	$\left[\begin{array}{cc}0& 1\\ -1.20088& 0.60129\end{array}\right]$	$\left[\begin{array}{c}1.00138\\ -0.99896\end{array}\right]$

The curve of parameter $\mathit{\theta }$ estimation error $\epsilon$ with time $t$ is shown in Figure 2.

Figure 2.

Parameter estimation error $\epsilon$ with time $t$ .

Example 3

For equation (59), in the same way, $\left\{u\left(t\right)\right\}$ is taken as an independent random signal sequence with zero mean and unit variance, and $\left\{v\left(t\right)\right\}$ as a white noise sequence with zero mean and variances ${\sigma }_1^2={0.1}^2$ and ${\sigma }_2^2={1.0}^2$ , respectively. Time is represented as $t$ , and parameter estimation error is calculated by $\epsilon =\widehat{\mathit{\theta }}-{\mathit{\theta }}_2.$ .

$$\{\begin{array}{c}x\left(t+1\right)=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -1& 0.2& 0.8\end{array}\right]x\left(t\right)+\left[\begin{array}{c}1\\ 1\\ -1\end{array}\right]u\left(t\right)\\ y\left(t\right)=\left[\begin{array}{ccc}1& 0& 0\end{array}\right]x\left(t\right)+v\left(t\right) \end{array}$$

(59)

First, true value of $\mathit{\theta }$ is calculated.

$\mathit{\theta }=\left[\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\end{array}\right]=\left[\begin{array}{c}-1\\ 0.2\\ 0.8\\ -2\\ 0.2\\ 1\end{array}\right],$

Let

$\left[\begin{array}{c}{\theta }_{11}\\ {\theta }_{12}\\ {\theta }_{13}\\ {\theta }_{21}\\ {\theta }_{22}\\ {\theta }_{23}\end{array}\right]=\left[\begin{array}{c}-1\\ 0.2\\ 0.8\\ -2\\ 0.2\\ 1\end{array}\right]$ .

Then algorithm in 2.3 is used to estimate parameters. When $t$ = 3500, the estimated value and error of $\mathit{\theta }$ are shown in Table 5.

Table 5.

Estimated value and error of parameter $\mathit{\theta }$ .

	${\mathit{\theta }}_{11}$	${\mathit{\theta }}_{12}$	${\mathit{\theta }}_{13}$	${\mathit{\theta }}_{21}$	${\mathit{\theta }}_{22}$	${\mathit{\theta }}_{23}$	$\mathit{\epsilon }$
True value	−1.00000	0.20000	0.80000	−2.00000	0.20000	1.00000
Estimated value ( ${\mathit{\sigma }}_1$ )	−1.00139	0.20076	0.80109	−2.00101	0.20049	1.00320	0.0038397
Estimated value ( ${\mathit{\sigma }}_2$ )	−1.00288	0.19930	0.80300	−2.00411	0.19149	1.00527	0.0116133

After obtaining estimated value of $\mathit{\theta }$ , the estimated values of parameters $\mathit{A}$ and $\mathit{b}$ are obtained, as shown in Table 6.

Table 6.

Estimated values and true values of parameters $\mathit{A}$ and $\mathit{b}$ .

	$\mathit{A}$	$\mathit{b}$
True value	$\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -1& 0.2& 0.8\end{array}\right]$	$\left[\begin{array}{c}1\\ 1\\ -1\end{array}\right]$
Estimated value ( ${\mathit{\sigma }}_1$ )	$\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -1.00139& 0.20076& 0.80109\end{array}\right]$	$\left[\begin{array}{c}1.00027\\ 1.00019\\ -1.00033\end{array}\right]$
Estimated value ( ${\mathit{\sigma }}_2$ )	$\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -1.00288& 0.19930& 0.80300\end{array}\right]$	$\left[\begin{array}{c}1.00052\\ 1.00043\\ -1.00069\end{array}\right]$

The curve of parameter $\mathit{\theta }$ estimation error $\epsilon$ with time $t$ is shown in Figure 3.

Figure 3.

Parameter estimation error $\epsilon$ with time $t$ .

It can be seen from the three examples and Tables 1, 2, 3, 4, 5, and 6 that estimation error of parameters increases with the increase of noise, but it is still acceptable, so that the accuracy of estimator is relatively high. It is easy to see from Figures 1, 2, and 3 that the estimation error of parameter $\mathit{\theta }$ decreases with the increase of time, and the smaller the noise variance $\sigma$ is, the faster the convergence speed of error curve is. The error $\epsilon$ will fluctuate greatly when noise variance ${\sigma }^2={1.0}^2$ and $t<500$ , but with the increase of time $t$ , it will show a downward trend, and finally achieve satisfactory results.

Conclusion

In this paper, an LS algorithm based on BCP is presented for parameter estimation of canonical state space model and a satisfactory result is achieved. The canonical state space model is transformed into vector function form suitable for LS algorithm, and then, augmented vector and auxiliary variable are introduced to realize unbiased estimation based on noise compensation; finally, the effectiveness of the estimator is verified by an example. This method of parameter estimation is universally useful, which can be extended to parameter estimation of other systems with noise, for example, parameter estimation of closed-loop system can be realized by selecting different augmented vectors, this is also our future research work. Therefore, the method presented in our paper has widespread application and practical significance.