A Hierarchical Approach for Joint Parameter and State Estimation of a Bilinear System with Autoregressive Noise

This paper is concerned with the joint state and parameter estimation methods for a bilinear system in the state space form, which is disturbed by additive noise. In order to overcome the difficulty that the model contains the product term of the system input and states, we make use of the hierarchical identification principle to present new methods for estimating the system parameters and states interactively. The unknown states are first estimated via a bilinear state estimator on the basis of the Kalman filtering algorithm. Then, a state estimator-based recursive generalized least squares (RGLS) algorithm is formulated according to the least squares principle. To improve the parameter estimation accuracy, we introduce the data filtering technique to derive a data filtering-based two-stage RGLS algorithm. The simulation example indicates the efficiency of the proposed algorithms.

In the literature of bilinear system identification, Verdult and Verhaegen applied the subspace identification techniques for multi-input multi-output bilinear systems by using the separable least squares principle [32].Larkowski et al. handled the parameter estimation problem of the diagonal bilinear errors-in-variables system and employed the bias-compensated least squares method to estimate the parameters [33].Hizir et al. transformed a bilinear system into its equivalent linear model and used the observer Kalman filter identification algorithm to identify the bilinear system [34].Vicario et al. extended the interaction matrices to bilinear systems to present the relationship between the system states and the measurement data and developed the intersection subspace algorithm for the bilinear system identification [35].
The recursive least squares (RLS) method is regarded as the common parameter estimation approach among many different parameter estimation techniques [36][37][38].Gan et al. proposed an efficient variable projection method to solve nonlinear least squares problems based on the matrix decomposition [39].Compared with the transfer function representation [40,41], the state space model can reflect the motion of the inner mechanism and cover state estimation [42][43][44].Many estimation methods such as sequential Monte Carlo methods [45] and subspace identification methods [46] have been extensively applied in many areas.Stroud et al. proposed a Bayesian method for the joint state and parameter estimation of state space models with additive Gaussian noise using the ensemble Kalman filter [47].Urteaga et al. presented a sequential Monte Carlo method for filtering and prediction of time-varying signals by fusing the information from candidate models instead of using model selection [48].Martino et al. presented the parallel particle filters based on the Bayesian model averaging principle for sequential tracking and online model selection [49].Schön et al. derived an expectation maximization algorithm under the maximum likelihood framework and gave the parameter estimates of nonlinear state space models [50].Li and Liu derived the input-output representation of the bilinear system through eliminating the state variables and proposed the filtering-based least squares iterative algorithm [51].However, the state variables in the system were eliminated, so the state estimation of the bilinear system was not studied.
The bilinear system involves the products of control inputs and state variables in the state equation.Thus, linear system identification methods cannot be utilized to compute the parameter estimates due to the unmeasurable states.Different from the work in [51], this paper focuses on the joint parameter and state estimation problem for a bilinear state space system with colored noise.The bilinear model can be regarded as a time-varying state space model [52].Thus, this paper presents a bilinear state estimator based on the Kalman filtering algorithm for computing the unknown states.Then, the state estimator-based recursive generalized least squares (RGLS) identification algorithm is developed for estimating the system parameters according to the least squares principle.For the purpose of improving the parameter estimation accuracy, the data filtering technique is introduced by transforming the bilinear model into several filtered submodels with smaller dimensions.The main contributions of this paper are listed as follows.

•
We present a bilinear state estimator on the basis of the Kalman filtering algorithm.

•
We derive a state estimator-based RGLS algorithm for joint state and parameter estimation of bilinear state space models on the basis of the hierarchical identification principle.

•
We derive a filtering-based two-stage RGLS (F-TS-RGLS) algorithm using the state estimates for improving the parameter estimation accuracy by introducing the data filtering technique.
The remainder of this paper is organized as follows.Section 2 describes the bilinear system and gives its identification model.Section 3 derives a state estimator-based RGLS identification algorithm.Section 4 presents a state estimator-based F-TS-RGLS algorithm by using the data filtering technique.Section 5 provides an illustrative example to demonstrate the effectiveness of the proposed algorithms.Finally, some conclusions are given in Section 6.

System Description
Some symbols are introduced for convenience."A =: X" or "X: = A" stands for "A is defined as X"; θt denotes the estimate of θ at time t; the symbol I (I n ) represents an identity matrix of an appropriate size (n × n); z denotes a unit forward shift operator like zx t = x t+1 and z −1 x t = x t−1 ; the superscript T symbolizes the matrix/vector transpose.
Recently, a state observer-based multi-innovation stochastic gradient algorithm and a state observer-based recursive least squares identification algorithm were presented for a bilinear system with white noise [53]: Zhang et al. derived a state filtering-based hierarchical identification algorithm and a state filtering-based forgetting factor recursive least squares algorithm for a bilinear system with white noise [54]: On the basis of the data filtering technique, a bilinear state observer-based multi-innovation extended stochastic gradient algorithm and a data filtering-based multi-innovation extended stochastic gradient algorithm were developed for a bilinear system with moving average noise [55]: where Different from the work in [53][54][55], this paper considers a single-input single-output bilinear system with autoregressive noise: where x t : = [x 1,t , x 2,t , • • • , x n,t ] T ∈ R n is the state vector, u t ∈ R and y t ∈ R are the system input and output data, w t ∈ R is the measurement noise, and A ∈ R n×n , B ∈ R n×n , f ∈ R n , and h ∈ R 1×n are the system parameters: Without loss of generality, assume that u t = 0, x t = 0, y t = 0, and w t = 0 for t 0.
Assumption 1.This paper assumes that the stochastic noise v t ∈ R is a Gaussian noise with zero mean and variance σ 2 .The colored noise w t can be a moving average process, an autoregressive process, or an autoregressive moving average process.In this paper, w t is approximated by an autoregressive process: Assumption 2. The system stability is the basis of system identification.This paper assumes that the bilinear system in ( 1) and ( 2) is stable, that is to say the bounded input leads to the bounded output, and the system is observable and controllable.
Assumption 3. System identification contains the determination of the system dimension and parameter estimation.This paper assumes that the orders n and n c of the system are known.The parameters a i , b ij , f i , and c i are to be identified from experimental data u t and y t .
Referring to the method in [53] and from ( 1) and ( 2), we have: Define the system parameter vectors θ, θ s , and θ n and the information vectors ϕ t , ϕ s,t , and ϕ n,t as: Inserting (4) into (2) gives: Then, Equation (3) can be written as: Equation ( 7) is the identification model of the bilinear system in ( 1) and ( 2).The objective of this paper is to derive a bilinear state estimator to obtain the state estimates and explore an efficient parameter identification method to generate highly-accurate parameter estimates and to improve the computational efficiency by using available experimental data.

The RGLS Algorithm Using the Bilinear State Estimates
In this section, a state estimator-based RGLS algorithm is employed to estimate the unknown states and parameters of the considered bilinear state space system.
Define a cost function: Using the least squares principle and minimizing J(θ) give the recursive least squares algorithm: The first identification difficulty occurs because only the observation data u t and y t are available, but ϕ t , ϕ x,t and ϕ xu,t contain the unknown x t .
Remark 1.To solve this problem, one method is to eliminate the state vector x t and to obtain a new identification model that involves only the system input and output variables [51].However, this cannot work because the parameter matrix B here is not the special matrix with many zero entries like in [51].Thus, it is necessary to study new state estimation methods to obtain the state estimate.
The second identification difficulty occurs because ϕ t also contains the noise variable {w t−i , i = 1, 2, . . ., n c }, so the parameter estimates θt cannot be computed by ( 8)- (10).
Remark 2. The solution here is to introduce the hierarchical identification method for obtaining the parameter estimates and state estimates, respectively.Based on the auxiliary model identification idea, we replace the unknown disturbance w t−i and the unknown state {x t−i , i = 1, 2, . . ., n} with their estimates ŵt−i and xt−i in the identification algorithms and define the estimated information vectors φt , φs,t , φx,t , φxu,t , φn,t of ϕ t , ϕ s,t , ϕ x,t , ϕ xu,t , ϕ n,t as: and the estimated parameter vectors θt , θs,t , ât , bt , ft , and θn,t of θ, θ s , a, b, f , and θ n as: According to the structure of ( 5), the estimate ŵt can be calculated through: According to the method in [55], we construct the following bilinear state estimator: Remark 3. Another problem occurs because the parameters A, B, f are unknown.Therefore, the state estimation in ( 17)-( 19) cannot be realized.The solution is to use ât , bt , and ft to construct the parameter estimates Ât , Bt , and ft of A, B, and f as: Then, replacing A, B, f , and θ n in ( 17)- (19) with their estimates Ât , Bt , ft , and θn,t obtains the bilinear state estimator: Replacing the unknown ϕ t in ( 8)- (10) with its estimate φt , and combining the bilinear state estimator in (20)- (23) give: Equations ( 11)-( 26) form the recursive generalized least squares (RGLS) algorithm based on the state estimates.
Remark 4. In order to improve the RGLS parameter estimation accuracy, the next section will construct a linear filter to filter the experimental data for estimating the system parameters.

The F-TS-RGLS Algorithm Using the Bilinear State Estimator
For the purpose of improving the parameter estimation accuracy of the RGLS algorithm, this section presents an F-TS-RGLS algorithm for the bilinear system by constructing a linear filter C(z) to filter the observation data of the system and to decompose the identification model in (7) into two sub-identification models.Multiplying both sides of ( 1) and ( 2) by C(z) obtains: Define the filtered input u f,t , the filtered output y f,t , and the filtered state x f,t as: x f,t : = C(z) Equations ( 27) and ( 28) can be transformed as: where x f,t : = [x 1f,t , x 2f,t , • • • , x (n−1)f,t , x nf,t ] T ∈ R n .From ( 32) and (33), we have: where the filtered information vector: Equations ( 33) and ( 6) can be expressed as the following two filtered sub-identification models: On the basis of the two sub-identification models in (36) and (37), define two cost functions: Minimizing J 1 (θ s ) and J 2 (θ n ) gives the following recursive relations: However, Equations ( 39)-( 44) cannot generate the parameter estimates θs,t and θn,t .Because the noise term w t and ϕ n,t are unmeasurable, the system state x t in the system information vector ϕ s,t is unmeasurable.In addition, C(z) is unknown, then y f,t , u f,t , x f,t , and ϕ f,t cannot be obtained.
Remark 5.The difficulty is overcome by utilizing the auxiliary model identification idea to take the place of the unknown variables with their estimates.Use the estimate ŵt−i of w t−i to construct the estimate of ϕ n,t : From ( 5), replacing ϕ s,t and θ s with their corresponding estimates φs,t and θs,t−1 obtains the estimate of w t : ŵt = y t − φT s,t θs,t−1 , where: By virtue of the Kalman filter, a state estimator in the bilinear form can be established for generating the state estimates.Then, utilizing the parameter estimates: to construct the estimates of C(z) gives: Since the linear filter is unknown, we use its estimate to filter the system input u t , the system output y t , and the system state x t to give the filtered estimates of u f,t , y f,t , and x f,t : Referring to the derivation of ( 21)-( 23), we can obtain the following state estimator for obtaining the estimate of the filtered state x f,t : By replacing x 1f,t−i , x f,t−i and u f,t−i in (35) with their estimates x1f,t−i , xf,t−i and ûf,t−i , define the estimated filtered information vector: Similarly, replacing the unmeasurable information vector ϕ f,t and ϕ n,t in ( 39)-( 44) with φf,t and φn,t obtains: L n,t = P n,t−1 φn,t [1 + φT n,t P n,t−1 φn,t ] −1 , (65) The procedures of calculating θs,t and θn,t by the filtering-based recursive generalized least squares algorithm in ( 45)-( 66) are listed as follows.
Compute the gain vector L s,t using (62) and the covariance matrix P s,t using (63).Update the parameter estimates θs,t using (61).7.
Increase t by one, and turn to Step 2.
Remark 6.By using the hierarchical identification principle, the unknown parameters and states can be estimated interactively.

Remark 7.
By introducing a linear filter, the original identification model is decomposed into two filtered sub-identification models.Then, the system parameters and the noise parameters are estimated interactively.
Remark 8.For identification methods, the large dimensions of the parameter vector θ cause a heavy computational cost.Thus, this paper uses the hierarchical identification principle for improving the computational efficiency of the proposed methods based on the decomposition.

Remark 9.
The number of the multiplication and addition operations is used to measure the computational cost of an algorithm.The computational cost of the RGLS and F-TS-RGLS algorithms is shown in Table 1.
Table 1.The computation efficiency of the recursive least squares (RLS) and filtering-based two-stage RGLS (F-TS-RGLS) algorithms.

Algorithms Number of Multiplications Number of Additions Total Flop
The difference of the computational cost of the RGLS and the F-TS-RGLS algorithms is: Remark 10.Since n 1 and n c 1, the difference N 1 − N 2 is positive.That is to say, the computational cost of the F-TS-RGLS algorithm is less than that of the RGLS algorithm.

Numerical Example
This section provides an example to illustrate the effectiveness of the proposed algorithms.Here, we consider a bilinear state space model: The parameter vector to be identified is: In the simulation, the input signal {u t } was chosen as a pseudo-random binary sequence with zero mean and unit variance and {v t } as a white noise sequence with zero mean and variances σ 2 .Set the data length L = 3000, and employ the proposed algorithms to estimate the parameters and states of this example system.The simulation input and output data are shown in Figure 1, which indicate that the bounded input signal leads to the bounded output.The parameter estimates and their estimation errors of the RGLS algorithm and the F-TS-RGLS algorithm are shown in Tables 2 and 3.The RGLS estimation error δ versus t is shown in Figure 2 with σ 2 = 0.10 2 and σ 2 = 0.15 2 .The F-TS-RGLS estimation error versus t is shown in Figure 3 with different variances.The predicted outputs and the true outputs are shown in Figures 4 and 5.The states x i,t and their estimates xi,t versus t are shown in Figure 6.From the simulation results in Tables 2 and 3 and Figures 1-6, we can draw the following conclusions.

•
The estimation errors of the RGLS algorithm and the F-TS-RGLS algorithm became smaller with the data length increasing.This means that the proposed algorithms are effective.

•
Under the same noise levels, the state estimator-based F-TS-RGLS algorithm could generate more accurate parameter estimates than the state estimator-based RGLS algorithm.

•
The F-TS-RGLS algorithm could decrease the dimensions of the covariance matrices and improve the computational efficiency.

•
The state estimates obtained from the bilinear state estimator could track their true values as t increased.

Conclusions
This paper presents a recursive parameter and state estimation algorithm on the basis of the hierarchical identification principle for bilinear systems.In this approach, the system parameters and states could be estimated interactively.The bilinear state estimator was derived based on the Kalman filter.Then, the state estimator-based RGLS algorithm and the state estimator-based F-TS-RGLS algorithm were proposed based on the data filtering technique and the decomposition-coordination principle.The simulation results showed that the state estimator-based F-TS-RGLS algorithm had higher parameter estimation accuracy compared with the state estimator-based RGLS algorithm.Moreover, the state estimator-based F-TS-RGLS algorithm could greatly improve the computational efficiency.The methods proposed in this paper could be combined the particle Monte Carlo methods [56,57] and some statistical methods to study the model selection and parameter estimation [58][59][60] for different systems [61][62][63][64][65][66] and could be applied to other fields [67][68][69][70][71][72] such as communication networks [73][74][75].