# Recursive Identification for Dynamic Linear Systems from Noisy Input-Output Measurements

Authors: Dan Fan; Kueiming Lo
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Errors-in-variables (EIV) model is a kind of model with not only noisy output but also noisy input measurements, which can be used for system modeling in many engineering applications. However, the identification for EIV model is much complicated due to the input noises. This paper focuses on the adaptive identification problem of real-time EIV models. Some derivation errors in an accuracy research of the popular Frisch scheme used for EIV identification have been pointed out in a recent study. To solve the same modeling problem, a new algorithm is proposed in this paper. A Moving Average (MA) process is used as a substitute for the joint impact of the mutually independent input and output noises, and then system parameters and the noise properties are estimated in the view of the time domain and frequency domain separately. A recursive form of the first step calculation is constructed to improve the calculation efficiency and online computation ability. Another advantage of the proposed algorithm is its applicableness to different input processes situations. Numerical simulations are given to demonstrate the efficiency and robustness of the new algorithm.


---

## Body

## 1. Introduction

In the field of engineering, modeling is an essential issue. In most cases, the systems are modeled by stochastic models in which the input signals are assumed to be measured exactly and all the disturbed noises are added to the output signals; that is, only the output measurements are noisy. These models are called “errors-in-equation models.” However, there are always signals beyond our control that also affect the input of the systems; some of them cannot be included in the output noises. Therefore, it is also necessary to consider the modeling problem for those systems with noisy input-output measurements, especially when we concern the actual physical laws of the process rather than the prediction of the future behaviour [1]. This kind of model whose input and output measurements are both containing noise is called “errors-in-variables (EIV) model [2].”The identification of EIV models has received a lot of attention during the past decades. By far, EIV models have been used in numerous applications, such as the modeling problems in econometrics, computer vision, biomedicine, chemical and image reconstruction, spectrum estimation, speech analysis, noise cancelation, and digital communications [3–8].In EIV models, the noise in input measurements cannot be equivalent to the output error, which makes the identification of EIV models much more difficult. The identifiability of EIV dynamic models was analyzed in [9, 10]. It is pointed out that EIV dynamic models cannot be uniquely identified from the second-order properties [9]. Thus specific prior knowledge is needed to achieve the identifiability. Once the identifiability is established, estimation algorithms can be developed [10]. Owing to the noisy input measurements in the EIV models, the standard least squares method for errors-in-equation models cannot yield consistent estimates anymore. To overcome this problem, a bias-compensated least squares (BCLSs) principle was proposed in [4]. On the basis of BCLS principle, various algorithms have been developed, such as the Frisch scheme-based algorithms [7], the KL algorithm [8], ECLS [9], BELS [10], and others in [11–15].Although there are such a number of approaches for identifying different EIV models, the convergence of the algorithms has always been a difficulty. Only a few literatures have tried to solve this problem [12, 15, 16]. In [16], the accuracy of Frisch scheme for EIV identification was analyzed, in which the estimates of the system parameters as well as the noise variance were both proved asymptotically Gaussian distributed by linearizing three primary equations in this scheme. This conclusion can be perceived as the theoretical support of the Frisch scheme-based algorithms. Based on this work, continued extensions and real applications have sprung up recently [17–20], which reaffirms the value of this particular analysis result. However, the analysis in [16] needs a condition that the estimates of the parameters are close to their true values, which is not clear how to be guaranteed. A counterexample that could not converge was present in [21], and in addition, some derivation errors of [16] were found and discussed at the same time. Furthermore, another method was provided to identify the EIV model in [21]. But comparing to the model concerned in [16], due to the difficulty of the identification problem, the one considered in [21] was a simpler one with a stronger condition that the input and output noise processes had the same variance, which has been hampering its application in some degree.The purpose of this paper is to consider how to avoid the restrict condition in [21]; in other words, we are trying to solve the same modeling problem as in [16], that is, to propose an identification algorithm for the modeling of dynamic EIV systems with independent input and output noises to estimate both the unknown system parameters and the noise signals. In order to achieve this purpose, we used a two-step method: in Step 1, the original model is rewritten into another form to get the system parameters in the time domain; in Step 2, the noise variances are calculated in the frequency domain. Moreover, the recursive form of the proposed method will be presented to improve its operational efficiency and enhance its online applicability.The structure of the paper is as follows. In Section2, the concerned model is described in detail. The new identification algorithm is presented in Section 3. Some simulations are given in Section 4 to illustrate the identification accuracy, the convergence rate, and the antinoise performance. Finally, conclusion remarks are given in Section 5.


## 2. Problem Formulation

A basic dynamic EIV system is shown in Figure1.Figure 1
A basic errors-in-variables system.Unlike the normal errors-in-equation model, as mentioned before, the EIV model has noise in both input measurements and output measurements. The immeasurable true input and output processesu0(t) and y0(t) are linked by a dynamic system, which can be a linear or a nonlinear system in different applications. So far, most of the related studies are focused on the linear systems, which is also the focus of research in this paper. The ARX(na,nb) model is considered here as
(1)A(z)y0(t)=B(z)u0(t),
where
(2)A(z)=1+a1z+⋯+anazna,B(z)=b1z+⋯+bnbznb,
are the polynomials in the backward shift operator z. The {a1,a2,…,ana,b1,b2,…,bnb}are the unknown system parameters to be identified, while the measured variables u(t) and y(t) are disturbed by the unknown noises u~(t) and y~(t). Thus the input and output measurements are
(3)u(t)=u0(t)+u~(t),y(t)=y0(t)+y~(t).After introducing the notations(4)θ=(a1,a2,…,ana,b1,b2,…,bnb)T,φ(t)=(-y(t-1),…,-y(t-na),u(t-1),…,u(t-nb))T,φ~(t)=(-y~(t-1),…,-y~(t-na),u~(t-1),…,u~(t-nb))T,
the EIV system can be described as the following model:
(5)y(t)=φ(t)Tθ+A(z)y~(t)-B(z)u~(t).
To ensure the identifiability, we list some assumptions first.(A1)
The EIV system is asymptotically stable, which means that there is no zero ofA(z) inside the unit circle.(A2)
The noisesu~(t) and y~(t) are mutually independent and also independent of the true input and output signals u0(t) and y0(t).(A3)
u~(t) and y~(t) are white noises with zero mean and independent variances λu and λy.The problem we need to solve is to estimate the system parameter vectorθ with the help of the measured regressor vector φ(t). Furthermore, considering that a noise process can be described by the mean and variance, to identify the zero-mean input and output noises is simplified to identify the variances. Therefore, except for the system parameters, we also want to estimate the output and input noise variances λy and λu. In the following section, we will give an algorithm with two independent steps to fulfill the two aspects of the estimate requirements.


## 3. Identification Algorithms

As mentioned before, the identification for EIV system is much more difficult because the input and output noises are unknown. For the EIV system described in Section2, to overcome the influence of the input noise, we will use another MA process {w(t)} as a substitute for the joint impact of the mutually independent input and output noises, as two mutually independent sequences of independent random variables can be represented as an MA process which has the same spectra with the two jointly sequences [22]. Then the system can be modified as an ARMAX model, and what need to do is changed to estimate the system parameters of the new model and to determine the variances of the input/output noises in terms of {w(t)}. Thus a two-step recursive estimation algorithm can be constructed to identify the system parameters θ and the noise variances λy and λu, respectively.Step  1. for the time t, we used the obtained estimation of w(t-1) to estimate the parameters and get the current estimates θ(t) and w(t).Step  2. These results are utilized to calculate the estimates of the noise variance λy(t) and λu(t).In the following, we will give the algorithm followed by proof.Step 1   (estimation for the unkown system parameterθ).
For convenience, denote the last two terms of (5) by v(t), that is,
(6)v(t)=A(z)y~(t)-B(z)u~(t),
where u~(t) and y~(t) are mutually independent with
(7)Ey~(t)=Eu~(t)=0,Ey~2(t)=λy,Ey~2(t)=λu.
Introduce an MA(nc) process
(8)w(t)=e(t)+c1e(t-1)+⋯+cnce(t-nc),
where {e(t)} is white noise with
(9)Ee(t)=0,Ee2(t)=λe,nc=max{na,nb}.
It can be shown that we can find a pair of {ci,0≤i≤nc} and λe such that {w(t)} and {v(t)} have the same spectra [22], which means that {v(t)} can be represented by {w(t)} in (8) as
(10)y(t)=φ(t)Tθ+w(t).
The {ci,0≤i≤nc} and λe are intermediate variables.
For the new model (10), denote a new parameter vector θ- and a new regressor vector φ-(t) by
(11)θ-=(θT,c1,c2,…,cnc)T,(12)φ-(t)=(φ(t)T,e(t-1),…,e(t-nc))T,
and then the EIV system (5) can be rewritten as
(13)y(t)=φ-(t)Tθ-+e(t).
In this step, we will give a recursive algorithm to identify (13).
The covariance matrix of the regressor vectorφ-(t) and output variables y(t) is denoted by
(14)Rφ-=Eφ-(i)φ-(i)T,rφ-y=Eφ-(i)y(i).
For convenience, introduce
(15)R^φ-(t)=∑i=1tφ-(i)φ-(i)T,(16)r^φ-y(t)=∑i=1tφ-(i)y(i).
Assume that the input {u(t)} is a stationary process; in the calculation, we can use the algebra means R^φ-(t)/t and r^φ-y(t)/t instead of the mathematical expectations Rφ- and rφ-y in (14), as by ergodicity, we have
(17)R^φ-(t)t→t→∞Rφ-,r^φ-y(t)t→t→∞rφ-y.Lemma 1 (Matrix Inversion Formula [23]).
For the matricesA∈Rn×n, C∈Rn×1, and D∈Rn×1, the inverse matrix of B=A+CDT is
(18)(A+CDT)-1=A-1-a-1A-1CDTA-1,
where a=1+DTA-1C.Theorem 2.
For system (13), under the assumptions (A1)–(A3), the parameter vector θ- can be estimated recursively as follows with a large P(0) and arbitrary θ-(t):(19)θ-(t)=θ-(t-1)+a-1(t)P(t-1)φ-^(t)×[y(t)-φ-^(t)Tθ-(t-1)],P(t)=P(t-1)-a-1(t)P(t-1)φ-^(t)φ-^(t)TP(t-1),ε(t)=y(t)-φ-^(t)Tθ-(t),φ-^(t)=(φ(t)T,ε(t-1),…,ε(t-nc))T,
where a(t)=1+φ-^(t)TP(t-1)φ-^(t).Proof.
Like the least squares methods, we use the covariance matrix for help. By multiplying withφ-(t) to the system model (13), we have
(20)Eφ-(t)y(t)-Eφ-(t)φ-(t)Tθ-=Eφ-(t)e(t),
with assumptions (A2) and (A3) which can be rewritten as
(21)Rφ-θ-=rφ-y.
Replacing Rφ- and rφ-y with R^φ-(t)/t and r^φ-y(t)/t in (15) and (16), respectively, as mentioned before, one has
(22)R^φ-(t)θ-=r^φ-y(t).
We note that (15) and (16) imply
(23)R^φ-(t)=R^φ-(t-1)+φ-(t)φ-(t)T,r^φ-y(t)=r^φ-y(t-1)+φ-(t)y(t).
Then on condition that R^φ-(t) is reversible, θ- is estimated as
(24)θ-(t)=R^φ--1(t)r^φ-y(t)=R^φ--1(t)(r^φ-y(t-1)+φ-(t)y(t))=R^φ--1(t)[(R^φ-(t)-φ-(t)φ-(t)T)θ-(t-1)+φ-(t)y(t)]=θ-(t-1)+R^φ--1(t)φ-(t)(y(t)-φ-(t)Tθ-(t-1)).
Using (23), we can calculate the vector θ-(t). But we note that there is an inverse operation of R^φ-(t) at each recursive step, which is a very time-consuming process. To avoid the inversing, introduce
(25)P(t)=R^φ--1(t),
and apply the Matrix Inversion Formula in Lemma 1 to (24); taking A=R^φ-(t),C=DT=φ-(t), we have
(26)P(t)=P(t-1)-P(t-1)φ-(t)φ-(t)TP(t-1)1+φ-(t)TP(t-1)φ-(t).
Moreover, by (26) it is clear that
(27)R^φ--1(t)φ-(t)=P(t-1)φ-(t)1+φ-(t)TP(t-1)φ-(t).
Taking (27) into (23), we have
(28)θ-(t)=θ-(t-1)+P(t-1)φ-(t)1+φ-(t)TP(t-1)φ-(t)·(y(t)-φ-(t)Tθ-(t-1)).
Noting that{e(t-1)} in (12) is unknown, φ-(t) cannot be constructed directly. This problem can be solved in the similar way as for RPLR (recursive pseudolinear regression) algorithm [22], that is, to form a substitute for φ-(t) as
(29)φ-^(t)=(φ(t)T,ε(t-1),…,ε(t-nc))T,
where ε(t-i)=y(i)-φ-^(i)Tθ-(i), i≥1. The proof of the substitution’s correctness is omitted (see details in [22]).
Then usingφ-^(t) defined by (29) to replace the φ-(t) in (26) and (28), we get Theorem 2 easily.Theorem2 gives a recursive algorithm to get the estimate of θ-: At time t-1, we store only the finite-dimensional information {θ-(t-1),P(t-1),φ-^(t-1)}. At time t, it is updated using (23), (27), and (29), which is done with a given fixed amount of operations, making it a high operational efficiency and suitable for online applications. Since θ-(t) is obtained, obviously θ(t) can be easily got by (11).Next we go to estimate the noise properties, which will be helpful in real applications such as the cascade system modeling.Step 2 (estimation for the noise variancesλy(t) and λu(t)).
To estimate the noise variancesλy, λu, we need to find the relationship between {v(t)} in (6) and {w(t)} in (8).
We know that the spectrumΦs(ω) of a signal {s(t)} is the Fourier transform of its covariance function Rs(τ) as
(30)Φs(ω)=∑τ=-∞∞Rs(τ)e-iτω,
where
(31)Rs(τ)=E-s(t)s(t-τ)=limN→∞1N∑t=1Ns(t)s(t-τ).
Since{w(t)} and {v(t)} have the same spectra, that is,
(32)Φw(ω)≡Φv(ω),
this means that for all τ=0,1,2,…,nc,
(33)Rv(τ)=Rw(τ).
Thus we can find the relationship between λy, λu, and λe by using the covariance functions Rv(τ) and Rw(τ).
At stept, an estimate of Rs(τ) can be used as
(34)Rsτ(t)=1t∑k=1ts(k)s(k-τ).
We switched to the notation Rvτ, Rwτ rather than Rv(τ), Rw(τ) to account for certain differences due to recursive step t.
Introduce(35)θaτ(t)=(a-τ(t),a-τ+1(t),…,a-τ+na(t))T,θbτ(t)=(b-τ+1(t),b-τ+2(t),…,b-τ+nb(t))T,θcτ(t)=(c-τ(t),c-τ+1(t),…,c-τ+nc(t))T,φy~(t)=(-y~(t),-y~(t-1),…,-y~(t-na))T,φu~(t)=(u~(t-1),…,u~(t-nb))T,φe(t)=(e(t),e(t-1),…,e(t-nc))T,
where
(36)a-k(t)={1,k=0,ak(t),0<k≤na,0,k>naork<0,b-k(t)={bk(t),0<k≤nb,0,k>nbork≤0,c-k(t)={1,k=0,ck(t),0<k≤nc,0,k>ncork<0,
are the transformation of the parameters in θ-(t).
Then according to (A2), (A3), (6), (8), and (34), we have
(37)Rwτ(t)=1t∑k=1tw(k)w(k-τ)=1t∑k=1tθc0(t)Tφe(k)θc0(t)Tφe(k-τ)=1t∑k=1tθc0(t)Tφe(k)φe(k)Tθcτ(t)=θc0(t)Tθcτ(t)λe.
Similarly we have
(38)Rvτ(t)=θa0(t)Tθaτ(t)λy+θb0(t)Tθbτ(t)λu.
Then it is clear that for τ=0,1,2,…,nc,
(39)Rvτ(t)=Rwτ(t)
is a set of linear equations about the unknown λy and λu, in which {θaτ(t),θbτ(t),θcτ(t)} is known at time t, and λe can be estimated by:
(40)λe(t)=Rε0(t)=1t∑k=1tε2(k)
by taking ε(t) in Theorem 2 as the estimate of e(t).
Then the estimatesλy(t) and λu(t) for the output/input noise variances λy and λu can be calculated by (37)–(40).


## 4. Simulation Examples

This section addresses some numerical evaluation of the identification algorithm presented in this paper. Matlab 7.7 is used to do the simulations. To demonstrate its validness to various EIV systems, we have chosen different signal processes as the true input variables{u0(t)}in each case: in Case A, a zero-mean Gaussian process is used; in Case B, a sawtooth signal is applied; in Case C, it is an ARMA process. The noise processes {u~(t)} and {y~(t)} in these cases are mutually uncorrelated white noise signals with zero-mean. The robustness of the algorithm is also tested, which is shown in Case C.Case A.
First we examine how well the algorithm works for systems with Gaussian input. Consider an EIV dynamic system withna=nb=2 and
(41)θ=(a1,a2,b1,b2)T=(-0.2,-0.15,0.3,-0.27)T.
It is easy to get the system as follows, which is denoted by System 1:
(42)System1:y0(t)-0.2y0(t-1)-0.15y0(t-2)=0.3u0(t-1)-0.27u0(t-2).
Let the input signal{u0(t)} be a zero-mean Gaussian process whose variance equals 1. Let the noise signals {u~(t)} and {y~(t)} be mutually uncorrelated white noise signals with λy=0.2, λu=0.5, which means a strong noise environment for the system.
The system is simulated forN=8000 steps. Calculation results are listed in Table 1, where the calculation error is defined by the standard deviation. Figures 2 and 3 show the system parameter and the noise variances estimates separately. Solid lines indicate the true values and dashed lines denote the corresponding estimates. Noting that the vertical coordinate scopes are very small in both figures, it can be seen easily that the estimates are converging fast to the true parameters.Table 1
Estimation results of System 1.
Parameter
True value
Estimation
a
1
-
0.2
0
-
0.1975
±
0.0025
a
2
-
0.15
-
0.148
3
±
6.9910
e
-
4
b
1
0.30
0.284
2
±
2.3776
e
-
4
b
2
-
0.27
-
0.25
74
±
7.5975
e
-
4
λ
y
0.20
0.192
2
±
4.5330
e
-
4
λ
u
0.50
0.5
125
±
0.0024Figure 2
Estimation result ofθ in System 1.Figure 3
Estimation result ofλyand λu in System 1.Case B.
Consider another system, System 2, with sawtooth input:(43)System2:y0(t)=z-2z2(1-0.9z)(1-0.8z)u0(t),
in which na=nb=nc=2, and
(44)θ=(a1,a2,b1,b2)T=(-1.7,0.72,1,-2)T.
This is the counterexample which was presented in [21] to show the unconvergency of the Frisch-based method analyzed in [16]. We use our proposed algorithm to identify this system under the same conditions; that is, the input sawtooth signal’s amplitude equals 1 and its frequency is 10 Hz; the noises’ variances are
(45)λy=λu=0.5.
The simulation results of the system parameters and the noise variances are all displayed in Figure4. The true values and estimates are also denoted by solid lines, and dotted lines respectively. We can see that the algorithm has an even better performance for this kind of EIV system. All the estimates converge to their corresponding true values consummately.Figure 4
Estimates of System 2 by the new recursive method.Case C.
The proposed algorithm in this paper is valid not only for EIV models with iid random input process but also for those whose input signals are more general such as the ARMA process. In order to examine it, we have considered the following system:(46)System3:y0(t)=2(1-0.5z)z1+0.2z-0.48z2u0(t),
where
(47)θ=(a1,a2,b1,b2)T=(0.2,-0.48,2,-1)T.
In this system, the input process{u0(t)} is an ARMA process:
(48)(1-5z)u0(t)=(1-0.3z)ξ(t)
with ξ(t) being a zero-mean iid Gaussian process whose variance is
(49)λξ=1.
The noise signals {u~(t)} and {y~(t)} are mutually uncorrelated white noise signals with
(50)λy=0.2,λu=0.5.
The estimation results are shown in Figures5 and 6. Similarly, the true values are marked by solid lines and the estimates of parameters are marked by dotted lines. We can see that for ARMA input process, the estimation can still converge to the true parameters quickly.Figure 5
Estimation result ofθ in System 3 with ARMA input process.Figure 6
Estimation ofλyand λu in System 3 with ARMA input process.In order to verify the robustness against noise, we further did several experiments with different signal-to-noise ratios (SNRs) in System 3. The noise processes used in the experiments are shown by the first three columns in Table2. Column 1 is the signal variables used to generate the input signals; columns 2 and 3 are the variables of input noise and output noise.Table 2
Comparison of the estimations for different SNRs.
λ
ξ
λ
u
λ
y
Average error
1
0.05
0.05
0.0025
1
0.2
0.2
0.0039
1
0.2
0.5
0.0040
1
0.5
0.2
0.0049
1
0.5
0.5
0.0052
1
1
1
0.0086
1
1
2
0.0064The average estimation error of System 3 (including the system parameter estimation and the noises estimation) is listed in column 4. It is clear that the performance of the proposed algorithm is keeping good when the noises increase (even when the noises are equal or larger than the input signals).


## 5.  Conclusions

This paper discussed the identification problem of dynamic errors-in-variables (EIV) systems. EIV model is very useful and has a wide engineering application range such as the modeling of cascade system, and camera calibration. In comparison with the usual errors-in-equation models, EIV model has a more troublesome noise problem with the input measurements being disturbed. Since several severe errors in the previous analysis of the attractive Frisch scheme identification approaches have been presented in the recent studies, we developed an adaptive algorithm to solve the same modeling problem. For the dynamic EIV model with mutually independent input and output noises, this two-step algorithm can not only estimate the system parameter vector as well as the noise variances with greater accuracy but also reduce the computational complexity significantly due to its recursive form. It has been shown by several simulation results that the presented algorithm demonstrates, as shown in the numerical simulations, a great accuracy, fast convergence speed, and good antinoise performance. Theoretical analysis of the proposed algorithm and the identification of some more complicated model such as EIV nonlinear models will be considered in future work.



---

*Source: 318786-2013-05-20.xml*