IDENTIFICATION OF NOISY INPUT-OUTPUT SYSTEM USING BIAS ...
IDENTIFICATION OF NOISY INPUT-OUTPUT SYSTEM USING BIAS-COMPENSATED LEAST-SQUARES METHOD Masato Ikenoue ∗ Shunshoku Kanae ∗∗ Zi-Jiang Yang ∗∗ Kiyoshi Wada ∗∗ ∗ Department of Electrical Engineering, Ariake National College of Technology, 150, Higashihagio-machi, Omuta, Fukuoka 836-8585, Japan Email: [email protected] ∗∗ Department of Electrical and Electronic System Engineering, Kyushu University, 6-10-1, Hakozaki, Higashi-ku, Fukuoka 812-8581, Japan Email: {jin, yoh, wada}@ees.kyushu-u.ac.jp
Abstract: In this paper, a new bias-compensated least-squares (BCLS) based algorithm is proposed for identification of noisy input-output system. It is well known that BCLS method is based on compensation of asymptotic bias on the least-squares (LS) estimates by making use of noise variances estimates. The main feature of the proposed algorithm is to introduce a generalized least-squares type estimator in order to obtain the good estimates of noise variances. The results of a simulated example indicate that the proposed algorithm provides good estimates. c 2005 IFAC Copyright ° Keywords: Estimation, Identification, Least-squares method
1. INTRODUCTION Many identification methods are based on the assumption that input measurement is noise-free. However, this condition is not satisfied in most practical situations. In the presence of input noise, those methods have been shown to give erroneous results. Several methods have been proposed to estimate unknown parameters of linear discrete-time system in the presence of input and output noises. Joint Output (JO) method (S¨oderstr¨ om, 1981) and Koopmans-Levin (KL) method (Fernando and Nicholson, 1985) require a priori knowledge about the values of variances or the ratio to measurements noises. Bias-compensated least-squares (BCLS) method is proposed by Sagara et al. (Sagara and Wada,
1977) and it has been extended by Wada et al. (Wada et al., 1990) to the input-output noise case without any a priori knowledge of noise variances. BCLS method based on compensation of asymptotic bias on the least-squares (LS) estimates by making use of noise variances estimates is very efficient method for estimation of noisy inputoutput system parameters. In recent years, BCLS method has been developed to improve the estimation accuracy and several recursive algorithms have been proposed (Eguchi et al., 1992; Jia et al., 2001). On the other hand, another method named biaseliminated least-squares (BELS) method has been proposed by Zheng et al. (Zheng and Feng, 1989) in which the different estimation method of asymptotic bias is used and further developed to
E[dt ] = 0, E[et ] = 0
be the efficient method (Zheng, 1999; Zheng, 2000; Zheng, 2001; Zheng, 2002) to treat bias problem in noisy input-output system identification. In this paper, a new BCLS based algorithm is proposed for identification of linear discrete-time system in the case where input and output measurements are corrupted by white noise. Since the unknown noise variances estimates are required for compensation of asymptotic bias of LS estimates, the estimation of these noise variances plays an important role in BCLS method. If the good estimates of noise variances are obtained, the estimation accuracy of the resulting BCLS estimates can be improved. For this purpose, a generalized least-squares type estimator is introduced in order to obtain the good estimates of noise variances. It is demonstrated that the improvement in the estimate accuracy of noise variances estimates (and the resulting BCLS estimates) is achieved. This paper is organized as follows. In section 2 and section 3, the problem statement is presented and the asymptotic bias of LS estimator is described. In section 4, the BCLS algorithm is derived for estimating unknown parameters of linear discretetime system in the presence of input and output noises and it can be learned that the unknown noise variances must be estimated in order to obtain consistent estimates of parameters. In section 5, the recursive algorithm of unknown noise variances are derived by introducing a generalized least-squares type estimator and in section 6, Jia et al.’s BCLS method (Jia et al., 2001) and Zheng’s BELS method (Zheng, 2001) are briefly described. The simulation results are presented in section 7 and finally section 8 gives the conclusion.
E[di dj ] = σd2 δi,j ,
(5)
E[ei ej ] =
σe2 δi,j
E[di ej ] = 0
(6) (7)
where δi,j is Kronecker’s delta. The true input ut is zero-mean stationary random process with finite variance, and ut , dt and et are assumed to be statistically independent of each other. Substituting (4) into (1) yields A(q −1 )zt = B(q −1 )wt + vt
(8)
where vt is a composite noise defined by vt = A(q −1 )et − B(q −1 )dt .
(9)
Define some vectors as θ T = [aT , bT ] = [a1 · · · an , b1 · · · bn ] pTt
= [−z Tt ,
(10)
wTt ]
= [−zt−1 · · · − zt−n , wt−1 · · · wt−n ] (11) q Tt = [−y Tt , uTt ] = [−yt−1 · · · − yt−n , ut−1 · · · ut−n ] (12) r Tt
= [−eTt , dTt ] = [−et−1 · · · − et−n , dt−1 · · · dt−n ]
(13)
then (4), (8) and (9) can be written as pt = q t + r t
(14)
zt = pTt θ
+ vt
(15)
vt = et −
r Tt θ
.
(16)
b of θ Let the equation error ξt for an estimate θ be defined as
2. PROBLEM STATEMENT Consider the parameter estimation problem of single-input single-output linear discrete-time system described as follows: A(q −1 )yt = B(q −1 )ut
(1)
b −1 )zt − B(q b −1 )wt ξt = A(q b = zt − pT θ t
(17)
b −1 ) and B(q b −1 ) are where the polynomials A(q defined by b −1 ) = 1 + b A(q a1 q −1 + · · · + b an q −n bb1 q −1 + · · · + bbn q −n b −1 ) = B(q
where ut and yt are the true input and output, q −1 is shift operator, q −1 ut = ut−1 , and the polynomials A(q −1 ) and B(q −1 ) are defined by
(18) (19)
and A(q
−1
B(q
−1
) = 1 + a1 q
−1
)=
−1
b1 q
+ · · · + an q
−n
+ · · · + bn q
−n
(2) .
(3)
Let zt and wt be the noise-corrupted measurements of yt and ut , respectively, i.e.
T bT = [b θ aT , b b ] = [b a1 · · · b an , bb1 · · · bbn ] . (20)
The least-squares estimate is given by Ã
zt = yt + et , wt = ut + dt
(4)
where et and dt are the output and input measurement noises, respectively. The measurement noises et and dt have the following statistical properties
bLS,N = θ
N X t=1
!−1 pt pTt
N X
pt zt .
(21)
t=1
From the assumption of et and dt , the composite noise vt defined by (9) is not white. Hence the
bLS,N has a bias asympleast-squares estimate θ totically. In the next section, the asymptotic bias induced by least-squares estimator is derived.
3. ASYMPTOTIC BIAS OF LEAST-SQUARES ESTIMATOR
N X
p t vt
where P N is the product moment matrix as follows:
PN =
!−1 pt pTt
.
(23)
t=1
Taking probability limit of above equation yields bLS,N = θ + h p lim θ N →∞
N 1 X pt vt N →∞ N t=1
h = R−1 pp p lim
1 N →∞ N
(29)
T b bLS,N = θ bLS,N−1 + P N−1 pN (zN −pN θ LS,N−1 ) θ T 1+pN P N−1 pN (30) T P N−1 pN pN P N−1 P N = P N−1 − . (31) 1+pTN P N−1 pN
Practically the variances of input and output noises σe2 and σd2 in (29) are unknown, it is necessary to estimate them.
5. ESTIMATION OF NOISE VARIANCES To estimate the noise variances of input and output noises σe2 and σd2 , a filter α(q −1 ) defined by the following equation is introduced
(25)
where Rpp = E[pt pTt ]. Using the assumption of et , dt and (14), (16), it is easily shown that
p lim
bBC,N = θ bLS,N − h bN θ bLS,N + N P N D θ bBC,N −1 . =θ
(24)
where h is the asymptotic bias of the least-squares bLS,N defined as estimate θ
N X
Hence bias-compensated least-squares estimate bBC,N is given by θ
(22)
t=1
ÃN X
(28)
bLS,N and P N are The recursive algorithm of θ obtained by the ordinary recursive least-squares algorithm
Substituting (15) into (21) yields ˆ LS,N = θ + P N θ
b N = −N P N Dθ . h
α(q −1 ) =
l X
αi q −i , (l ≥ n) .
(32)
i=0
Now, let the filtered signal for equation error ξt be defined as pt vt = E[pt vt ]
t=1
= E[(q t + r t )(et −
r tT θ)]
= −E[r t r Tt ]θ = −Dθ
(26)
where D = diag{σe2 I n ; σd2 I n } and I n is n × n identity matrix. From (24), (25) and (26), the asymptotic bias h can be expressed as follows:
ξet = α(q −1 )ξt =
l X
αi ξt−i .
Minimizing the sum of squared ξet yields the estib N of θ mator ϕ
bN = ϕ
ÃN X
eTt et p p
!−1 Ã N X
t=1
bLS,N − θ h = p lim θ N →∞ = R−1 pp {−Dθ}
.
(27)
4. BIAS-COMPENSATED LEAST-SQUARES METHOD
! et zet p
(34)
t=1
et is the filtered signal for pt and zet is the where p filtered signal for zt as et = α(q −1 )pt = p
l X
αi pt−i
(35)
αi zt−i .
(36)
i=0
zet = α(q −1 )zt = From (24), it can be expected that a consistent estimate of θ can be obtained by compensating for the asymptotic bias h in the least-squares bLS,N . From (27), estimate of the asympestimate θ totic bias h becomes
(33)
i=0
l X i=0
b N can be considered as generalThe estimator ϕ ized least-squares type estimator, and the recurb N is obtained by sive algorithm of ϕ
bN = ϕ b N −1 + ϕ
e N −1 p eN (e eTN ϕ b N −1 ) P zN − p (37) T e e P N −1 p eN 1+p
Hn =
N
eN P
e e N −1 eN p eTN P e N −1 − P N −1 p =P T e e P N −1 p eN 1+p
l X
αi2 I n
i=0
(38)
+
N
l−j n−1 XX
h i αi αi+j (S n )(j) + (S Tn )(j)
j=1 i=0
·
where eN = P
ÃN X
0Tn−1 0 Sn = I n−1 0n−1
!−1 et p eTt p
.
(39)
t=1
¸ (49)
and 0n is an n × 1 zero vector. Finally,
It follows from (14) that eTt θ + vet zet = p
1 b fN p lim N→∞ N ( l ) ´ ³ X 2 2 T b a,N =σe αi +ρ a+p lim ϕ
(40)
where vet is the filtered signal for vt as l X
vet = α(q −1 )vt =
N→∞
i=0
b Ta,N H n a+σd2 p lim ϕ b Tb,N H n b +σe2 p lim ϕ N→∞ N→∞ αi vt−i .
(41)
i=0
Now, let ½ αi =
=
N 2 X
N X
t=1
t=1
be ξt =
N X
b N )2 eTt ϕ (e zt − p
b TN zet vet − ϕ
t=1
N X
et vet p
gbN =
bLS,N . ξbt = zt − pTt θ
N
p lim
N→∞
N→∞
N→∞
From (50) and (54), the estimates of input and output noise variances σe2 and σd2 can be obtained by solutions of system of equations
N 1 X et vet = E[e p pt vet ] N →∞ N t=1 · ¸ ρ 2 e = −σe − Dθ 0n
"
p lim
(45)
where l−n X
# αi αi+n
i=0
(46) e = diag{σ 2 H n ; σ 2 H n } D e d
1 gbN N
b TLS,N a+σd2 p lim b =σe2 +σe2 p lim a bLS,N b . (54)
(44)
and
i=0
(53)
T
i=0
i=0
(52)
It follows from (50) that
1 X zet vet = E[e zt vet ] N →∞ N t=1 Ã l ! X = σe2 αi2 + ρT a
αi αi+2 , · · · ,
ξbt2
where ξbt is the residual of the least-squares estibLS,N defined by mate θ
p lim
αi αi+1 ,
N X t=1
(43)
Taking probability limit of (42) yields
l−2 X
(51)
(42)
t=1
be bN . eTt ϕ ξ t = zet − p
ρT =
1 , i=0 0 , 1≤i≤l
b N in (34) becomes the leastthen the estimator ϕ bLS,N , and fbN becomes the sum squares estimate θ of squared residual defined by
b bN where ξet is the residual of the estimator ϕ defined by
" l−1 X
(50)
b TN = [b b Tb,N ]. where ϕ ϕTa,N , ϕ
Define the the sum of squared residual fbN as fbN =
(48)
(47)
T b b BC,N−1 1+b aTLS,N a bLS,N b bBC,N−1 T T b a,N H n a b BC,N−1 ϕ b b,N H n b α bN + ϕ bBC,N−1 · ¸ 1 gbN = N fbN
#"
c2 σ e c2 σ d
#
(55)
where α bN =
l X
¢ ¡ b BC,N−1 + ϕ b a,N (56) αi2 + ρT a
i=0 T T bT θ aBC,N , b bBC,N ] . BC,N = [b
(57)
6. THE PREVIOUSLY PROPOSED BIAS-COMPENSATION PRINCIPLE BASED METHODS
b ¯ LS,N formula to Rp¯p¯, the last two components of θ may have two expressions. Then the following expression can be established from (64) and (65)
In this section, Jia et al.’s BCLS method (Jia et al., 2001) and Zheng’s BELS method (Zheng, 2001) are briefly described, which are different estimation methods of noise variances. 6.1 Jia et al.’s BCLS method b defined by Consider an auxiliary estimator φ N Ã b = φ N
N X
!−1 pt−1 pTt−1
t=1
N X
pt−1 zt .
(58)
N X b )zt−1 (zt − pTt−1 φ N
(59)
then the following expression can be obtained
N→∞
(66)
N →∞
where Rpx = E[pt xTt ], r xz = E[xt zt ]. R−1 pp,1 and R−1 are composed of the first n and the last n pp,2 −1 columns of Rpp , respectively. From (54) and (66), the estimates of input and output noise variances σe2 and σd2 can be obtained by solutions of system of equations T b BC,N−1 b 1+b aTLS,N a bLS,N b bBC,N−1 T −1 b bT R b −1 a b b b R R R px pp,1 BC,N−1 px pp,2 bBC,N−1 " # gbN /N = bLS,N bT θ bxz − R r
#"
c2 σ e c2 σ d
#
(67)
px
t=1
p lim
bLS,N = r xz − RTpx p lim θ
"
t=1
Define b jN by b jN =
−1 2 T σe2 RTpx R−1 pp,1 a + σd Rpx Rpp,2 b
Zheng’s BELS method can be extended for colored output noise (Zheng, 2002).
1b jN N
b T a + σ 2 p lim φ b T b (60) = σe2 p lim φ a,N d b,N N→∞
b T = [φ bT , φ b T ]. where φ N a,N b,N From (54) and (60), the estimates of input and output noise variances σe2 and σd2 can be obtained by solutions of system of equations "
b BC,N−1 1+b aTLS,N a T b b BC,N−1 φa,N a · ¸ 1 gbN = jN N b
T b bLS,N b bBC,N−1 T b b φ bBC,N−1 b,N
#"
c2 σ e c2 σ d
#
(61)
6.2 Zheng’s BELS method Define an argumented parameter vector as ¯ T = [θ T , cT ], cT = [bn+1 , bn+2 ] = 0T . (62) θ 2 The corresponding auxiliary regression vector is given by ¯ Tt = [pTt , xTt ], xTt = [wt−n−1 , wt−n−2 ] . (63) p ¯ is given by The least-squares estimate of θ b ¯ LS,N = R−1 r pz p lim θ p¯p¯ ¯
(64)
b ¯ LS,N = θ ¯ + R−1 D ¯ ¯θ p lim θ p¯p¯
(65)
N →∞ N →∞
7. SIMULATION RESULTS
N→∞
¯ = ¯ Tt ], r pz where Rp¯p¯ = E[¯ pt p pt zt ] and D ¯ = E[¯ 2 diag{D; σd I 2 }. Applying the matrix inversion
By computer simulation, the proposed BCLS algorithm is compared with Jia et al.’s BCLS algorithm, Zheng’s BELS algorithm and LS algorithm. Consider the following second-order system: 0.169901q −1 +0.143831q −2 B(q −1 ) = . (68) A(q −1 ) 1−1.575157q −1 +0.606531q −2 The noise free input ut is white signal with variance σu2 = 1. The noise variance on input side is set as σd2 = 0.1 which yields SNR = 10 log10 (σu2 /σd2 ) = 10 [dB]. The noise variance on output side is set as σe2 = 0.3987 which yields SNR = 10 log10 (σy2 /σe2 ) = 10 [dB]. The filter α(q −1 ) is designed as l = 3, α0 = 1, α1 = 2, α2 = 2, α3 = 1. Computer simulation for comparison are carried out through M = 100 independent runs with a data length of 4000. Fig. 1 gives a plot of the RMSE which is defined by v u M b u 1 X kθ k,t −θk2 RMSE = 20 log10 t [dB] (69) M kθk2 k=1
bk,t denotes the estimate of θ at time step where θ t in the kth independent run. The mean values of the estimates of σe2 and σd2 are shown in Fig. 2 and Fig. 3, respectively. Simulation results indicate that LS method gives biased results. On the contrary, the proposed BCLS method, Jia et al.’s method and Zheng’s
LS algorithm
10
0
10
RMSE[dB]
RMSE[dB]
Zheng's BELS algorithm
20
10
-20
10
10
-40
10
10
10 0
1000
2000
3000
20
0
-20
-40
4000
0
Jia et al.'s BCLS algorithm 20
10
0
10
RMSE[dB]
RMSE[dB]
10
-20
10
10
10
-40
10
10 0
1000
2000
1000
2000
3000
4000
Proposed BCLS algorithm
3000
20
0
REFERENCES -20
-40
4000
0
1000
2000
3000
4000
Fig. 1. RMSE of parameter estimates. Zheng's BELS algorithm 0.6 0.4 0.2 0 0
500
1000
1500 2000 2500 Jia et al.'s BCLS algorithm
3000
3500
4000
0
500
1000
1500 2000 2500 Proposed BCLS algorithm
3000
3500
4000
0
500
1000
1500
3000
3500
4000
0.6 0.4 0.2 0
0.6 0.4 0.2 0 2000
2500
Fig. 2. The mean values of the estimates of σe2 . Zheng's BELS algorithm 0.2 0.1 0 -0.1 -0.2 0
500
1000
1500 2000 2500 Jia et al.'s BCLS algorithm
3000
3500
4000
0
500
1000
1500 2000 2500 Proposed BCLS algorithm
3000
3500
4000
0
500
1000
1500
3000
3500
4000
0.2 0.1 0 -0.1 -0.2
0.2 0.1 0 -0.1 -0.2 2000
new BCLS based algorithm has been proposed by introducing a generalized least-squares type estimator. Since the proposed approach can give the good estimates of the noise variances, the estimation accuracy of the resulting BCLS estimates can be improved. It is demonstrated that the proposed method can give consistent parameter estimates via simulation results.
2500
Fig. 3. The mean values of the estimates of σd2 . method can give consistent estimates. Especially, since the proposed algorithm provides the good estimates of noise variances σe2 and σd2 compared with Jia et al.’s algorithm and Zheng’s algorithm, the resulting BCLS estimates are more accurate than those obtained with Jia et al.’s algorithm and Zheng’s algorithm.
8. CONCLUSION In this paper, the method of consistent estimation of noisy input-output system has been studied. A
Eguchi, M., K. Wada and S. Sagara (1992). Identification of pulse transfer function in the presence of input and output noise. IFAC international workshop on ACQP’92 2, 463–470. Fernando, K. V. and H. Nicholson (1985). Identification of linear systems with input and output noise: the koopmans-levin method. IEE Proc. Control Theory and Applications 132, 30–36. Jia, L. J., M. Ikenoue, C. Z. Jin and K. Wada (2001). On bias compensated least squares method for noisy input-output system identification. Proc. of 40th IEEE Conf. on Decision and Control pp. 3332–3337. Sagara, S. and K. Wada (1977). On-line modified least-squares parameter estimation on linear discrete dynamic systems. Int. J. Control 25(3), 329–343. S¨oderstr¨om, T. (1981). Identification of stochastic linear systems in presence of input noise. Automatica 17(5), 713–725. Wada, K., M. Eguchi and S. Sagara (1990). Estimation of pulse transfer function via biascompensated least-squares method in the presence of input and output noise. Systems Science 16(3), 57–70. Zheng, W. X. (1999). On least-squares identification of stochastic linear systems with noisy input-output data. Int. J. Adaptive Control and Signal Processing 13, 131–143. Zheng, W. X. (2000). Unbiased identification of stochastic linear systems from noisy input and output measurements. Proc. of 39th IEEE Conf. on Decision and Control pp. 2710–2715. Zheng, W. X. (2001). Fast adaptive iir filtering with noisy input and output data. Proc. of 6th International Symposium on Signal Processing and Its Applications 1, 307–310. Zheng, W. X. (2002). A bias correction method for identification of linear synamic errorsin-variables models. IEEE Trans. Autmatic Control 47(7), 1142–1147. Zheng, W. X. and C. B. Feng (1989). Unbiased estimation of linear systems in the presence of input and output noise. Int. J. Adaptive Control and Signal Processing 3, 231–251.
Abstract: In this paper, a new bias-compensated least-squares (BCLS) based algorithm is proposed for identification of noisy input-output system. It is well known that BCLS method is based on compensation of asymptotic bias on the least-squares (LS) estimates by making use of noise variances estimates. The main feature of the proposed algorithm is to introduce a generalized least-squares type estimator in order to obtain the good estimates of noise variances. The results of a simulated example indicate that the proposed algorithm provides good estimates. c 2005 IFAC Copyright ° Keywords: Estimation, Identification, Least-squares method
1. INTRODUCTION Many identification methods are based on the assumption that input measurement is noise-free. However, this condition is not satisfied in most practical situations. In the presence of input noise, those methods have been shown to give erroneous results. Several methods have been proposed to estimate unknown parameters of linear discrete-time system in the presence of input and output noises. Joint Output (JO) method (S¨oderstr¨ om, 1981) and Koopmans-Levin (KL) method (Fernando and Nicholson, 1985) require a priori knowledge about the values of variances or the ratio to measurements noises. Bias-compensated least-squares (BCLS) method is proposed by Sagara et al. (Sagara and Wada,
1977) and it has been extended by Wada et al. (Wada et al., 1990) to the input-output noise case without any a priori knowledge of noise variances. BCLS method based on compensation of asymptotic bias on the least-squares (LS) estimates by making use of noise variances estimates is very efficient method for estimation of noisy inputoutput system parameters. In recent years, BCLS method has been developed to improve the estimation accuracy and several recursive algorithms have been proposed (Eguchi et al., 1992; Jia et al., 2001). On the other hand, another method named biaseliminated least-squares (BELS) method has been proposed by Zheng et al. (Zheng and Feng, 1989) in which the different estimation method of asymptotic bias is used and further developed to
E[dt ] = 0, E[et ] = 0
be the efficient method (Zheng, 1999; Zheng, 2000; Zheng, 2001; Zheng, 2002) to treat bias problem in noisy input-output system identification. In this paper, a new BCLS based algorithm is proposed for identification of linear discrete-time system in the case where input and output measurements are corrupted by white noise. Since the unknown noise variances estimates are required for compensation of asymptotic bias of LS estimates, the estimation of these noise variances plays an important role in BCLS method. If the good estimates of noise variances are obtained, the estimation accuracy of the resulting BCLS estimates can be improved. For this purpose, a generalized least-squares type estimator is introduced in order to obtain the good estimates of noise variances. It is demonstrated that the improvement in the estimate accuracy of noise variances estimates (and the resulting BCLS estimates) is achieved. This paper is organized as follows. In section 2 and section 3, the problem statement is presented and the asymptotic bias of LS estimator is described. In section 4, the BCLS algorithm is derived for estimating unknown parameters of linear discretetime system in the presence of input and output noises and it can be learned that the unknown noise variances must be estimated in order to obtain consistent estimates of parameters. In section 5, the recursive algorithm of unknown noise variances are derived by introducing a generalized least-squares type estimator and in section 6, Jia et al.’s BCLS method (Jia et al., 2001) and Zheng’s BELS method (Zheng, 2001) are briefly described. The simulation results are presented in section 7 and finally section 8 gives the conclusion.
E[di dj ] = σd2 δi,j ,
(5)
E[ei ej ] =
σe2 δi,j
E[di ej ] = 0
(6) (7)
where δi,j is Kronecker’s delta. The true input ut is zero-mean stationary random process with finite variance, and ut , dt and et are assumed to be statistically independent of each other. Substituting (4) into (1) yields A(q −1 )zt = B(q −1 )wt + vt
(8)
where vt is a composite noise defined by vt = A(q −1 )et − B(q −1 )dt .
(9)
Define some vectors as θ T = [aT , bT ] = [a1 · · · an , b1 · · · bn ] pTt
= [−z Tt ,
(10)
wTt ]
= [−zt−1 · · · − zt−n , wt−1 · · · wt−n ] (11) q Tt = [−y Tt , uTt ] = [−yt−1 · · · − yt−n , ut−1 · · · ut−n ] (12) r Tt
= [−eTt , dTt ] = [−et−1 · · · − et−n , dt−1 · · · dt−n ]
(13)
then (4), (8) and (9) can be written as pt = q t + r t
(14)
zt = pTt θ
+ vt
(15)
vt = et −
r Tt θ
.
(16)
b of θ Let the equation error ξt for an estimate θ be defined as
2. PROBLEM STATEMENT Consider the parameter estimation problem of single-input single-output linear discrete-time system described as follows: A(q −1 )yt = B(q −1 )ut
(1)
b −1 )zt − B(q b −1 )wt ξt = A(q b = zt − pT θ t
(17)
b −1 ) and B(q b −1 ) are where the polynomials A(q defined by b −1 ) = 1 + b A(q a1 q −1 + · · · + b an q −n bb1 q −1 + · · · + bbn q −n b −1 ) = B(q
where ut and yt are the true input and output, q −1 is shift operator, q −1 ut = ut−1 , and the polynomials A(q −1 ) and B(q −1 ) are defined by
(18) (19)
and A(q
−1
B(q
−1
) = 1 + a1 q
−1
)=
−1
b1 q
+ · · · + an q
−n
+ · · · + bn q
−n
(2) .
(3)
Let zt and wt be the noise-corrupted measurements of yt and ut , respectively, i.e.
T bT = [b θ aT , b b ] = [b a1 · · · b an , bb1 · · · bbn ] . (20)
The least-squares estimate is given by Ã
zt = yt + et , wt = ut + dt
(4)
where et and dt are the output and input measurement noises, respectively. The measurement noises et and dt have the following statistical properties
bLS,N = θ
N X t=1
!−1 pt pTt
N X
pt zt .
(21)
t=1
From the assumption of et and dt , the composite noise vt defined by (9) is not white. Hence the
bLS,N has a bias asympleast-squares estimate θ totically. In the next section, the asymptotic bias induced by least-squares estimator is derived.
3. ASYMPTOTIC BIAS OF LEAST-SQUARES ESTIMATOR
N X
p t vt
where P N is the product moment matrix as follows:
PN =
!−1 pt pTt
.
(23)
t=1
Taking probability limit of above equation yields bLS,N = θ + h p lim θ N →∞
N 1 X pt vt N →∞ N t=1
h = R−1 pp p lim
1 N →∞ N
(29)
T b bLS,N = θ bLS,N−1 + P N−1 pN (zN −pN θ LS,N−1 ) θ T 1+pN P N−1 pN (30) T P N−1 pN pN P N−1 P N = P N−1 − . (31) 1+pTN P N−1 pN
Practically the variances of input and output noises σe2 and σd2 in (29) are unknown, it is necessary to estimate them.
5. ESTIMATION OF NOISE VARIANCES To estimate the noise variances of input and output noises σe2 and σd2 , a filter α(q −1 ) defined by the following equation is introduced
(25)
where Rpp = E[pt pTt ]. Using the assumption of et , dt and (14), (16), it is easily shown that
p lim
bBC,N = θ bLS,N − h bN θ bLS,N + N P N D θ bBC,N −1 . =θ
(24)
where h is the asymptotic bias of the least-squares bLS,N defined as estimate θ
N X
Hence bias-compensated least-squares estimate bBC,N is given by θ
(22)
t=1
ÃN X
(28)
bLS,N and P N are The recursive algorithm of θ obtained by the ordinary recursive least-squares algorithm
Substituting (15) into (21) yields ˆ LS,N = θ + P N θ
b N = −N P N Dθ . h
α(q −1 ) =
l X
αi q −i , (l ≥ n) .
(32)
i=0
Now, let the filtered signal for equation error ξt be defined as pt vt = E[pt vt ]
t=1
= E[(q t + r t )(et −
r tT θ)]
= −E[r t r Tt ]θ = −Dθ
(26)
where D = diag{σe2 I n ; σd2 I n } and I n is n × n identity matrix. From (24), (25) and (26), the asymptotic bias h can be expressed as follows:
ξet = α(q −1 )ξt =
l X
αi ξt−i .
Minimizing the sum of squared ξet yields the estib N of θ mator ϕ
bN = ϕ
ÃN X
eTt et p p
!−1 Ã N X
t=1
bLS,N − θ h = p lim θ N →∞ = R−1 pp {−Dθ}
.
(27)
4. BIAS-COMPENSATED LEAST-SQUARES METHOD
! et zet p
(34)
t=1
et is the filtered signal for pt and zet is the where p filtered signal for zt as et = α(q −1 )pt = p
l X
αi pt−i
(35)
αi zt−i .
(36)
i=0
zet = α(q −1 )zt = From (24), it can be expected that a consistent estimate of θ can be obtained by compensating for the asymptotic bias h in the least-squares bLS,N . From (27), estimate of the asympestimate θ totic bias h becomes
(33)
i=0
l X i=0
b N can be considered as generalThe estimator ϕ ized least-squares type estimator, and the recurb N is obtained by sive algorithm of ϕ
bN = ϕ b N −1 + ϕ
e N −1 p eN (e eTN ϕ b N −1 ) P zN − p (37) T e e P N −1 p eN 1+p
Hn =
N
eN P
e e N −1 eN p eTN P e N −1 − P N −1 p =P T e e P N −1 p eN 1+p
l X
αi2 I n
i=0
(38)
+
N
l−j n−1 XX
h i αi αi+j (S n )(j) + (S Tn )(j)
j=1 i=0
·
where eN = P
ÃN X
0Tn−1 0 Sn = I n−1 0n−1
!−1 et p eTt p
.
(39)
t=1
¸ (49)
and 0n is an n × 1 zero vector. Finally,
It follows from (14) that eTt θ + vet zet = p
1 b fN p lim N→∞ N ( l ) ´ ³ X 2 2 T b a,N =σe αi +ρ a+p lim ϕ
(40)
where vet is the filtered signal for vt as l X
vet = α(q −1 )vt =
N→∞
i=0
b Ta,N H n a+σd2 p lim ϕ b Tb,N H n b +σe2 p lim ϕ N→∞ N→∞ αi vt−i .
(41)
i=0
Now, let ½ αi =
=
N 2 X
N X
t=1
t=1
be ξt =
N X
b N )2 eTt ϕ (e zt − p
b TN zet vet − ϕ
t=1
N X
et vet p
gbN =
bLS,N . ξbt = zt − pTt θ
N
p lim
N→∞
N→∞
N→∞
From (50) and (54), the estimates of input and output noise variances σe2 and σd2 can be obtained by solutions of system of equations
N 1 X et vet = E[e p pt vet ] N →∞ N t=1 · ¸ ρ 2 e = −σe − Dθ 0n
"
p lim
(45)
where l−n X
# αi αi+n
i=0
(46) e = diag{σ 2 H n ; σ 2 H n } D e d
1 gbN N
b TLS,N a+σd2 p lim b =σe2 +σe2 p lim a bLS,N b . (54)
(44)
and
i=0
(53)
T
i=0
i=0
(52)
It follows from (50) that
1 X zet vet = E[e zt vet ] N →∞ N t=1 Ã l ! X = σe2 αi2 + ρT a
αi αi+2 , · · · ,
ξbt2
where ξbt is the residual of the least-squares estibLS,N defined by mate θ
p lim
αi αi+1 ,
N X t=1
(43)
Taking probability limit of (42) yields
l−2 X
(51)
(42)
t=1
be bN . eTt ϕ ξ t = zet − p
ρT =
1 , i=0 0 , 1≤i≤l
b N in (34) becomes the leastthen the estimator ϕ bLS,N , and fbN becomes the sum squares estimate θ of squared residual defined by
b bN where ξet is the residual of the estimator ϕ defined by
" l−1 X
(50)
b TN = [b b Tb,N ]. where ϕ ϕTa,N , ϕ
Define the the sum of squared residual fbN as fbN =
(48)
(47)
T b b BC,N−1 1+b aTLS,N a bLS,N b bBC,N−1 T T b a,N H n a b BC,N−1 ϕ b b,N H n b α bN + ϕ bBC,N−1 · ¸ 1 gbN = N fbN
#"
c2 σ e c2 σ d
#
(55)
where α bN =
l X
¢ ¡ b BC,N−1 + ϕ b a,N (56) αi2 + ρT a
i=0 T T bT θ aBC,N , b bBC,N ] . BC,N = [b
(57)
6. THE PREVIOUSLY PROPOSED BIAS-COMPENSATION PRINCIPLE BASED METHODS
b ¯ LS,N formula to Rp¯p¯, the last two components of θ may have two expressions. Then the following expression can be established from (64) and (65)
In this section, Jia et al.’s BCLS method (Jia et al., 2001) and Zheng’s BELS method (Zheng, 2001) are briefly described, which are different estimation methods of noise variances. 6.1 Jia et al.’s BCLS method b defined by Consider an auxiliary estimator φ N Ã b = φ N
N X
!−1 pt−1 pTt−1
t=1
N X
pt−1 zt .
(58)
N X b )zt−1 (zt − pTt−1 φ N
(59)
then the following expression can be obtained
N→∞
(66)
N →∞
where Rpx = E[pt xTt ], r xz = E[xt zt ]. R−1 pp,1 and R−1 are composed of the first n and the last n pp,2 −1 columns of Rpp , respectively. From (54) and (66), the estimates of input and output noise variances σe2 and σd2 can be obtained by solutions of system of equations T b BC,N−1 b 1+b aTLS,N a bLS,N b bBC,N−1 T −1 b bT R b −1 a b b b R R R px pp,1 BC,N−1 px pp,2 bBC,N−1 " # gbN /N = bLS,N bT θ bxz − R r
#"
c2 σ e c2 σ d
#
(67)
px
t=1
p lim
bLS,N = r xz − RTpx p lim θ
"
t=1
Define b jN by b jN =
−1 2 T σe2 RTpx R−1 pp,1 a + σd Rpx Rpp,2 b
Zheng’s BELS method can be extended for colored output noise (Zheng, 2002).
1b jN N
b T a + σ 2 p lim φ b T b (60) = σe2 p lim φ a,N d b,N N→∞
b T = [φ bT , φ b T ]. where φ N a,N b,N From (54) and (60), the estimates of input and output noise variances σe2 and σd2 can be obtained by solutions of system of equations "
b BC,N−1 1+b aTLS,N a T b b BC,N−1 φa,N a · ¸ 1 gbN = jN N b
T b bLS,N b bBC,N−1 T b b φ bBC,N−1 b,N
#"
c2 σ e c2 σ d
#
(61)
6.2 Zheng’s BELS method Define an argumented parameter vector as ¯ T = [θ T , cT ], cT = [bn+1 , bn+2 ] = 0T . (62) θ 2 The corresponding auxiliary regression vector is given by ¯ Tt = [pTt , xTt ], xTt = [wt−n−1 , wt−n−2 ] . (63) p ¯ is given by The least-squares estimate of θ b ¯ LS,N = R−1 r pz p lim θ p¯p¯ ¯
(64)
b ¯ LS,N = θ ¯ + R−1 D ¯ ¯θ p lim θ p¯p¯
(65)
N →∞ N →∞
7. SIMULATION RESULTS
N→∞
¯ = ¯ Tt ], r pz where Rp¯p¯ = E[¯ pt p pt zt ] and D ¯ = E[¯ 2 diag{D; σd I 2 }. Applying the matrix inversion
By computer simulation, the proposed BCLS algorithm is compared with Jia et al.’s BCLS algorithm, Zheng’s BELS algorithm and LS algorithm. Consider the following second-order system: 0.169901q −1 +0.143831q −2 B(q −1 ) = . (68) A(q −1 ) 1−1.575157q −1 +0.606531q −2 The noise free input ut is white signal with variance σu2 = 1. The noise variance on input side is set as σd2 = 0.1 which yields SNR = 10 log10 (σu2 /σd2 ) = 10 [dB]. The noise variance on output side is set as σe2 = 0.3987 which yields SNR = 10 log10 (σy2 /σe2 ) = 10 [dB]. The filter α(q −1 ) is designed as l = 3, α0 = 1, α1 = 2, α2 = 2, α3 = 1. Computer simulation for comparison are carried out through M = 100 independent runs with a data length of 4000. Fig. 1 gives a plot of the RMSE which is defined by v u M b u 1 X kθ k,t −θk2 RMSE = 20 log10 t [dB] (69) M kθk2 k=1
bk,t denotes the estimate of θ at time step where θ t in the kth independent run. The mean values of the estimates of σe2 and σd2 are shown in Fig. 2 and Fig. 3, respectively. Simulation results indicate that LS method gives biased results. On the contrary, the proposed BCLS method, Jia et al.’s method and Zheng’s
LS algorithm
10
0
10
RMSE[dB]
RMSE[dB]
Zheng's BELS algorithm
20
10
-20
10
10
-40
10
10
10 0
1000
2000
3000
20
0
-20
-40
4000
0
Jia et al.'s BCLS algorithm 20
10
0
10
RMSE[dB]
RMSE[dB]
10
-20
10
10
10
-40
10
10 0
1000
2000
1000
2000
3000
4000
Proposed BCLS algorithm
3000
20
0
REFERENCES -20
-40
4000
0
1000
2000
3000
4000
Fig. 1. RMSE of parameter estimates. Zheng's BELS algorithm 0.6 0.4 0.2 0 0
500
1000
1500 2000 2500 Jia et al.'s BCLS algorithm
3000
3500
4000
0
500
1000
1500 2000 2500 Proposed BCLS algorithm
3000
3500
4000
0
500
1000
1500
3000
3500
4000
0.6 0.4 0.2 0
0.6 0.4 0.2 0 2000
2500
Fig. 2. The mean values of the estimates of σe2 . Zheng's BELS algorithm 0.2 0.1 0 -0.1 -0.2 0
500
1000
1500 2000 2500 Jia et al.'s BCLS algorithm
3000
3500
4000
0
500
1000
1500 2000 2500 Proposed BCLS algorithm
3000
3500
4000
0
500
1000
1500
3000
3500
4000
0.2 0.1 0 -0.1 -0.2
0.2 0.1 0 -0.1 -0.2 2000
new BCLS based algorithm has been proposed by introducing a generalized least-squares type estimator. Since the proposed approach can give the good estimates of the noise variances, the estimation accuracy of the resulting BCLS estimates can be improved. It is demonstrated that the proposed method can give consistent parameter estimates via simulation results.
2500
Fig. 3. The mean values of the estimates of σd2 . method can give consistent estimates. Especially, since the proposed algorithm provides the good estimates of noise variances σe2 and σd2 compared with Jia et al.’s algorithm and Zheng’s algorithm, the resulting BCLS estimates are more accurate than those obtained with Jia et al.’s algorithm and Zheng’s algorithm.
8. CONCLUSION In this paper, the method of consistent estimation of noisy input-output system has been studied. A
Eguchi, M., K. Wada and S. Sagara (1992). Identification of pulse transfer function in the presence of input and output noise. IFAC international workshop on ACQP’92 2, 463–470. Fernando, K. V. and H. Nicholson (1985). Identification of linear systems with input and output noise: the koopmans-levin method. IEE Proc. Control Theory and Applications 132, 30–36. Jia, L. J., M. Ikenoue, C. Z. Jin and K. Wada (2001). On bias compensated least squares method for noisy input-output system identification. Proc. of 40th IEEE Conf. on Decision and Control pp. 3332–3337. Sagara, S. and K. Wada (1977). On-line modified least-squares parameter estimation on linear discrete dynamic systems. Int. J. Control 25(3), 329–343. S¨oderstr¨om, T. (1981). Identification of stochastic linear systems in presence of input noise. Automatica 17(5), 713–725. Wada, K., M. Eguchi and S. Sagara (1990). Estimation of pulse transfer function via biascompensated least-squares method in the presence of input and output noise. Systems Science 16(3), 57–70. Zheng, W. X. (1999). On least-squares identification of stochastic linear systems with noisy input-output data. Int. J. Adaptive Control and Signal Processing 13, 131–143. Zheng, W. X. (2000). Unbiased identification of stochastic linear systems from noisy input and output measurements. Proc. of 39th IEEE Conf. on Decision and Control pp. 2710–2715. Zheng, W. X. (2001). Fast adaptive iir filtering with noisy input and output data. Proc. of 6th International Symposium on Signal Processing and Its Applications 1, 307–310. Zheng, W. X. (2002). A bias correction method for identification of linear synamic errorsin-variables models. IEEE Trans. Autmatic Control 47(7), 1142–1147. Zheng, W. X. and C. B. Feng (1989). Unbiased estimation of linear systems in the presence of input and output noise. Int. J. Adaptive Control and Signal Processing 3, 231–251.
