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

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. 1 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, 2 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 9 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][11][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. 13 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. 14 For Hammerstein autoregressive moving average (ARMAX) systems with autoregressive variables, Zhang et al. proposed an RLS identification method based on the BCP. 15 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, 16 then, proposed an LS algorithm based on BCP for parameter estimation of MISO system under white noise. 17 Mu et al. proposed a globally convergent identification algorithm, modified BCP, and gave a consistent estimation of noise variance rather than introducing auxiliary vectors. 18 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. 19 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. 20 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)  where, 2 R n×n is and yðtÞ 2 R are input and output of the system, respectively, xðtÞ 2 R n is the unobservable state vector, and vðtÞ 2 R is an uncorrelated random white noise with zero mean. Assume that the order n is known, uðtÞ ¼ 0, yðtÞ ¼ 0 and vðtÞ ¼ 0 when t ≤ 0. A,b is the system parameter matrix and parameter vector to be identified, respectively. Since A is a companion matrix and c has a standard structure, that is, the first element of 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 Y ðt þ nÞ can be obtained by deducing yðtÞ ¼ cxðtÞ þ vðtÞ in equations (1) [20]: We have Due to the special structures of c and A, we have Q 0 ¼ I n , which is an unit matrix. The first nequations of equation (2) can be rewritten as (8) Transform (8) to get (9) Substituting (9) into equation (2) gives (10) yðt þ nÞ � V T ðtþnÞðcA n Þ T þvðtþnÞ: Obtain So, equation (10) can be written as equation (17) yðt Replacing t in equation (17) with t � n yields equation (18) yðtÞ So, the vector function form suitable for LS method is obtained Replacing t in equation (20) with t � n yields equation (21) Y ðtÞ ¼ ΦðtÞθ þṼ ðtÞ According to the principle of LS, the loss function can be defined as (22) LS estimation can be obtained by solving the minimization of loss function If the above equations are used to estimate b θ LS directly, and then A and b are deduced, parameter estimation of equation (1) can be achieved. However, it can be seen from equation (14) that noiseṽðt þ nÞ 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) Substituting ff ðnÞb r fv ðnÞ: Taking probability limit of (25) yields It can be seen from equation (26) that LS estimator is biased, and bias term is R �1 ff r fv caused by noise. In order to eliminate the influence of bias term, equation (29) In order to obtain an unbiased estimate of b θ, a constant is needed to replace b r fv ðnÞ to compensate for the influence of bias. Assuming that system operates in an open-loop mode, it means that inputs uðiÞ and vðiÞ are not statistically relevant, then equation (28) can be rewritten as equation (30) where, b r yv ðnÞ ¼ It can be seen from (31) that the key is to find b r yv ðnÞ, so augmented regression vector ψðiÞ and augmented parameter vector θ are introduced to satisfy ψ T ðiÞθ ¼ b f T ðiÞθ, which is commonly used in BCP. It is worth noting that no matter how different the structures of θ and ψðiÞ are, the calculation of b r yv ðnÞ is similar. 19 Therefore, the method described below is universal.
Introduce an auxiliary vector ξðiÞ 2 R m ðm ≥ nÞ which satisfies equation (32) and a nonsingular matrix Construct an augmented regression vector ψðiÞ and an augmented parameter vector θ by introducing auxiliary vector ξðiÞ and nonsingular matrix F , this step is the most important step, the selection of ψðiÞ and θ makes calculation of b r yv becomes possible The θ and ψðiÞ satisfy equation (35) Rewrite yðtÞ ¼ f T ðiÞθ þṽðtÞ as yðtÞ ¼ ψ T ðtÞθ þṽðtÞ: It is easy to yield equation (36) from equation (34)

Similar to equations (24) and (26), equations (39) and (40) can be obtained in the same way
Pre-multiplying Z T on both sides of equation (40) and using (38) yield Thus, the estimation of r yv can be determined by the weighted LS as equation (45) b r yv ðnÞ ¼ where W LS is a positive definite weighting matrix. The choice of W LS does not affect the consistency properties of b r yv ðnÞ, but may have a considerable impact on the statistical accuracy of b r yv ðnÞ.
One simple choice is W LS ¼ I . 19 In the case of m ¼ n, if b D is non-singular, b r yv ðnÞ simply becomes as equation (46) Similar to equations (24) and (26), equation (47) can be obtained from (40) in the same way By substituting (45) into (47), consistent estimation of the augmented parameter vector θ can be obtained, finally, the global consistent estimation of original parameter θ can be obtained by equation (36) and rewritten as equation Above is the unified BCP method. Through choosing appropriately ξðiÞ and F , 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 Then, b r yv can be obtained according to equation (45), and the global unbiased estimate of the parameter θ 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 θ and parameters A, b, 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.
The matrix M in parameter estimation of canonical state space model based on least squares can be simplified by (50) as

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), fuðtÞg is taken as an independent random signal sequence with zero mean and unit variance, and fvðtÞg as a white noise sequence with zero mean and variances σ 2 1 ¼ 0:1 2 and σ 2 2 ¼ 1:0 2 , respectively. Time is represented as t and parameter estimation error is calculated by ε ¼ b θ � θ 2 , which is similar to the example in Ref.
First, true value of θ is calculated: Then algorithm in 2.3 is used to estimate parameters. When t = 3500, the estimated value and error of θ are shown in Table 1: After obtaining estimated value of θ, the estimated values of parameters A and b are obtained, as shown in Table 2.
The curve of parameter θ estimation error ε with time t is shown in Figure 1.

Example 2
For equations (58), in the same way, fuðtÞg is taken as an independent random signal sequence with zero mean and unit variance, and fvðtÞg as a white noise sequence with zero mean and variances σ 2 1 ¼ 0:1 2 and σ 2 2 ¼ 1:0 2 , respectively. Time is represented as t, and parameter estimation error is First, true value of θ is calculated. Then algorithm in 2.3 is used to estimate parameters. When t = 3500, the estimated value and error of θ are shown in Table 3: After obtaining estimated value of θ, the estimated values of parameters A and b are obtained, as shown in Table 4.
The curve of parameter θ estimation error ε with time t is shown in Figure 2.

Example 3
For equation (59), in the same way, fuðtÞg is taken as an independent random signal sequence with zero mean and unit variance, and fvðtÞg as a white noise sequence with zero mean and variances σ 2 1 ¼ 0:1 2 and σ 2 2 ¼ 1:0 2 , respectively. Time is represented as t, and parameter estimation error is calculated by ε ¼ b θ � θ 2 :.
First, true value of θ is calculated. .
Then algorithm in 2.3 is used to estimate parameters. When t = 3500, the estimated value and error of θ are shown in Table 5.
After obtaining estimated value of θ, the estimated values of parameters A and b are obtained, as shown in Table 6.
The curve of parameter θ estimation error ε with time t is shown in Figure 3.
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 θ decreases with the increase of time, and the smaller the noise variance σ is, the faster the convergence speed of error curve is. The error ε will fluctuate greatly when noise   variance σ 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.  Table 6. Estimated values and true values of parameters A and b.