# Unbiased Minimum Variance Estimation for Discrete-Time Systems with Measurement Delay and Unknown Measurement Disturbance

Authors: Yu Guan; Xinmin Song
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper addresses the state estimation problem for stochastic systems with unknown measurement disturbances whose any prior information is unknown and measurement delay resulting from the inherent limited bandwidth. For such complex systems, the Kalman-like one-step predictor independent of unknown measurement disturbances is designed based on the linear unbiased minimum variance criterion and the reorganized innovation analysis approach. One simulation example shows the effectiveness of the proposed algorithms.


---

## Body

## 1. Introduction

In recent years, networked control systems have attracted much attention and much work has been done due to the wide applications in communication systems [1, 2], fault detection [3, 4], and sensor [5]. However, these networks are usually unreliable and may lead to measurement delays due to inherent limited bandwidth. According to different kinds of time delay, many results have sprung up, such as control input delay [6, 7], state-dependent delay [8], state-independent delay [9–11], output delay [12, 13], communication delay [14], distributed-delay [15–17], and time-varying delay [18, 19]. In addition, the unknown disturbances in system modeling or external environment are ubiquitous feature in the piratical systems. The measurement delays and unknown observation disturbances can influence the performance or even results in systems instability. For these reasons, it is not surprising that the study of the state estimation problem for systems with time-delays and unknown observation disturbances has been an enthusiasm for a large number of scholars.The early work on the discrete-time systems with time delay has been investigated by system augmentation [20] or partial difference Riccati equation approach [21]. In [21], the partial difference Riccati equation approach is used to settle the measurement delay problem for linear systems. In order to lessen the computational cost (compared with state augmentation approach and partial difference Riccati equation approach), [22] proposes the reorganized innovation approach, by calculating two standard Riccati difference equations of the same dimension as the original system; the authors solve the finite horizon estimation problem for measurement delayed systems. In [23], the linear minimum mean square estimation filter for systems with measurement delay is calculated in terms of two Riccati difference equations and one Lyapunov difference equation. It should be pointed that all the aforementioned estimators do not consider the disturbance in observation.In practice, the unknown disturbances in system modeling or external environment are another ubiquitous feature [24–26]. The early works on unbiased minimum variance problems with observation disturbances can be traced back to [27]. For these unknown disturbances without any prior knowledge, unbiased minimum variance filter instead of Kalman filter is more effective in tracing the true state. Reference [28] obtains a more general unbiased estimator using the approach in [27] and presents the convergence analysis of the estimator. In [29], the authors considers global optimality of unbiased estimator based on [27, 28]. In [30], the authors considers event-based state estimation of linear dynamic systems with unknown inputs. Different from the aforementioned methods about linear discrete-time systems with observation disturbance, [31, 32] solve the state estimation with partially observed inputs problem without adding unbiased constraint. However, [28–33] only consider the state estimation problem of the state equation with unknown input. It should be pointed that measurement disturbances are ubiquitous feature in practical systems and the estimation problems for discrete-time systems with unknown measurement disturbance are also important. Motivated by the preceding works about measurement delay and unknown measurement disturbance, we study linear unbiased estimation for discrete-time systems with measurement delay and unknown measurement disturbance in this paper. By employing the linear unbiased minimum variance approach and the reorganized innovation analysis approach, one Kalman-like one-step ahead predictor is derived in terms of two Riccati difference equations of the same dimension with the state model, which may reduce the computation when delay is large, compared with the classic state augmentation approach [20] and partial difference Riccati equation approach [21, 34]. Further, just like the stability of the classical Kalman filter [35] and the mean square stability analysis [36], we develop a parallel to obtain the stability properties of the proposed unbiased estimator under standard conditions.The organization of this paper is as follows. In Section2, we present the problem statement, some assumptions and remarks. In Section 3, we deduce the state estimation according to the reorganized innovation analysis approach, then we obtain the finite and infinite horizon filter based on the linear unbiased minimum variance criterion, respectively. In Section 4, a numerical example is given to illustrate the effectiveness of the proposed approach. In Section 5, we provide some concluding remarks.Notation.Throughout this paper, the superscripts “-1” and “T” represent the inverse and transpose of a matrix. Rn denotes the n-dimensional Euclidean space. Rn×m is the set of all n×m real matrices. δij=0 for i≠j and δii=1. L{y(s)}s=0k denotes the linear subspace spanned by the measurement sequence {y(0),…,y(k)}.  λi(X) denotes the ith eigenvalue of a square matrix X. Furthermore, the mathematical expectation operator is denoted by E.


## 2. Problems Statement and Preliminary

Consider the following linear system:(1)xk+1=Axk+nk(2)y0k=B0xk+C0u0k+v0k(3)y1k=B1xk-d+C1u1k+v1k,k≥dwhere x(k)∈Rn, y0(k)∈Rm0, and y1(k)∈Rm1 are the state, current, and delayed measurement, respectively. n(k) is the process noise, v0(k) and v1(k) are the measurement noises. u0(k)∈Rp0 and u1(k)∈Rp1 are the unknown measurement disturbances without any prior knowledge. For simplicity of presentation, we assume that A, B0, B1, C0, and C1 are constant matrices with suitable dimensions even though the later development and results can be easily adapted to the time-varying case.Assumption 1.
n ( k ), v0(k), and v1(k) are uncorrelated white noises of zero mean and covariance as(4)EnknTj=Qδkj,Ev0kv0Tj=R0δkj,Ev1kv1Tj=R1δkjAssumption 2.
The initial statex(0) is uncorrelated with n(k), v0(k), and v1(k) and satisfies(5)Ex0=μ0,Ex0-μ0x0-μ0T=P0Assumption 3.
r a n k C 0 = p 0 < m 0;  rankC1=p1<m1.To the estimation problems of stochastic systems, Assumptions1 and 2 are general. Assumption 3 guarantees the existence of the following predictors designed.For convenience, the measurementy(k) can be rewritten as follows:(6)yk=y0k,0≤k<dy0ky1k,k≥dProblem. For the given measurements {y(k)}k=0M, our aim is to design a minimum variance unbiased Kalman-like one-step predictor x^(k+1). Further, we will consider the infinite horizon predictor design x^(k+1).Remark 4.
As for time delay systems, we can settle the estimation problem by using the state augmentation approach or partial difference Riccati equation approach. However the augmented approach or partial difference Riccati equation approach may bring expensive computational cost when the delayd is large [22]. In the following, we will deduce the estimation problem based on the reorganized innovation approach and linear minimum variance unbiased criterion instead of the augmented approach.


## 3. Main Results

### 3.1. Finite Horizon Estimation

We note thaty1(k) is an additional measurement of the state x(k-d) which is obtained at time instant k with time delay d, so the measurement y(k) consists of time delay when k≥d. Apparently, the linear space L{y(s)}s=0k contains the same information as L{Y1(s)}s=0k-d,{Y0(s)}s=k-d+1k, where the new observations Y1(s) and Y0(s) are provided as follows:(7)Y1s=y0sy1s+d,0≤s≤k-d(8)Y0s=y0s,k-d<s≤kObviously, Y1(s) and Y0(s) satisfy(9)Y1s=Bxs+CUs+Vs(10)Y0s=B0xs+C1u0s+V0swhere(11)B=B0B1,C=C000C1,Us=u0su1s+d,Vs=v0sv1s+d,V0s=v0sIt is obvious that the new measurements Y0(s) and Y1(s) are delay-free and the associated measurement noises V0(s) and V(s) are white noises with zero mean and covariance matrices RV0(s)=R0, RV(s)=R=diag⁡{R0,R1}.Theorem 5.
Whenk≥d, based on linear minimum variance unbiased criterion, we produce a recursive state estimator decoupling with the disturbance for system (1), (7), and (8) in the Kalman-like form:(12)x^s+1,0=Ax^s,0+K0sY0s-B0x^s,0,k-d<s<kx^k-d+1,0=x^k-d+1,1where(13)K0s=G0s-Λ0TsC0TD0-1sD0s=B0P0sB0T+R0G0s=AP0sB0TΛ0Ts=C0TD0-1sC0-1C0TD0-1sG0TsThe prediction error covariance matrix P0(s+1) is computed by(14)P0s+1=K0sD0sK0Ts-G0sK0Ts-K0sG0Ts+AP0sAT+Q(15)P0k-d+1=P1k-d+1As for x^(k-d+1,1) and P1(k-d+1), they are obtained by(16)x^s+1,1=Ax^s,1+K1sY1s-Bx^s,1,0≤s≤k-d(17)x^0,1=μ0,where(18)K1s=G1s-Λ1TsCTD1-1s(19)D1s=BP1sBT+R(20)G1s=AP1sBT(21)Λ1Ts=CTD1-1sC-1CTD1-1sG1TsThe prediction error covariance matrix P1(s+1) is computed by(22)P1s+1=K1sD1sK1Ts-G1sK1Ts-K1sG1Ts+AP1sAT+Q(23)P10=P0.Proof.
When0<s≤k-d, from (1) and (16), we have the prediction error equation as follows:(24)x~s+1,1=xs+1-x^s+1,1=Axs+ns-Ax^s,1-K1sY1s-Bx^s,1=A-K1sBx~s,1+ns-K1sVs-K1sCUsAssume x^(s,1) to be unbiased; in order to guarantee that x^(s+1,1) be an unbiased estimate of x(s+1), we must have Ex~(s+1,1)=0, then we have(25)K1sC=0Substituting (25) into (24), (24) can be rewritten as(26)x~s+1,1=A-K1sBx~s,1+ns-K1sVsFrom (26), the prediction error covariance matrix P1(s+1) is computed by(27)P1s+1=Ex~s+1,1x~Ts+1,1=A-K1sBP1sA-K1sBT+Q+K1sRK1Ts=K1sD1sK1Ts-G1sK1Ts-K1sG1Ts+AP1sAT+Qwhere D1(s) and G1(s) are given by (19) and (20). In order to minimize the estimation error variance (27) under the constraint (25), we introduce an auxiliary equation(28)Js=trP1s+1+trK1sCΛs+trΛTsCTK1Tswhere Λ(s) is the Lagrange multipliers. Taking the derivatives with respect to K1(s) equal to zero yields(29)K1sD1s+ΛsCT-G1s=0Combining (25) and (29) gives the matrix equation(30)D1sCCT0K1TsΛs=G1Ts0
Equation (30) has a unique solution if and only if the coefficients is nonsingular. Due to Assumption 3 and the fact that D1(s) is nonsingular, the coefficient matrix of (30) is nonsingular. Obviously, premultiplying left- and right-hand sides of (30) by the inverse of coefficient matrix yields (19) and (21). When k-d<s≤k, the proof is similar to the case 0≤s≤k-d, so we omit it here.Remark 6.
The system dimension will get more and more higher with the increase of time delay d; thus, computational complexity will get higher when using the classical state augmentation approach. But when we adopt the reorganized innovation analysis approach, we only need to solve two Riccati equations which have the same order of the state equation, so it can greatly reduce the computational complexity compared with the classical augmentation approach when time delay d is large.


### 3.2. Infinite Horizon Estimation

In this subsection, we consider the steady unbiased estimator design. First, let us present the lemma below as the initial step for the convergence analysis.Lemma 7.
If there existK∈Rn×(m0+m1) such that(31)KC=0,λiA-KB<1for i=1,2,…,n, then the sequence {P1(s+1)} abided by (22) is bounded for any s given any initial condition 0≤{P1(0)}≤∞.Proof.
Let us consider a suboptimal unbiased filter as follows:(32)x^subs+1,1=Ax^s,1+KsubY1s-Bx^s,1,0≤s≤k-dwhere KsubC=0,λi[A-KsubB]<1. Then the following state estimation error x~sub(s+1,1)=x(s+1)-x^sub(s+1,1) is given by(33)x~subs+1,1=A-KsubBx~subs,1+ns-KsubVswith the associated covariance matrix P1sub(s+1) being given by(34)P1subs+1=A-KsubBP1subsA-KsubBT+Q+KsubRKsubTThus P1sub(s+1) is bounded for any nonnegative initial condition due to the fact that λi[A-KsubB]<1. Comparing the above suboptimal estimator to the designed optimal estimator, the optimality tell us that P1(s+1)≤P1sub(s+1) for the same initial value. This proves the boundedness of P1sub(s+1).Theorem 8.
If there existK∈Rn×(m0+m1) such that KC=0,λi[A-KB]<1, and A,Q1/2 is stabilizable, then {P1(s+1)} abided by (22) converges to a unique fixed point P1 for any initial condition, where P1 is computed by(35)P1=K1D1K1T-G1K1T-K1G1T+AP1AT+Qand(36)K1=G1-Λ1TCTD1-1(37)D1=BP1BT+R(38)G1=AP1BT(39)Λ1T=CTD1-1C-1CTD1-1G1T.When 0<s≤k-d, we have the state predictor as follows:(40)x^s+1,1=Ax^s,1+K1Y1s-Bx^s,1(41)x^0,1=μ0When k-d<s≤k, the state predictor is as follows:(42)x^s+1,0=Ax^s,0+K0sY0s-B0x^s,0(43)x^k-d+1,0=x^k-d+1,1,where(44)K0s=D0-1sG0Ts-C0Λ0TsTD0s=B0P0sB0T+R0G0s=AP0sB0TΛ0Ts=C0TD0-1sC0-1C0TD0-1sG0TsThe prediction error covariance matrix P0(s+1) is computed by(45)P0s+1=K0sD0sK0Ts-G0sK0Ts-K0sG0Ts+AP0sAT+Q(46)P0k-d+1=P1.Proof.
By simple calculation, the Riccati equation in (22) can be formulated as(47)P1s+1=A-K1sBP1sA-K1sBT+Q+K1sRK1Ts.Then according to [35, 37] and Lemma 7, the convergence analysis can be readily obtained and we omitted here.Remark 9.
Under the condition thatKC=0,λi[A-KB]<1, and A,Q1/2 is stabilizable, we have obtained the steady-state filter (40). This is important when one desires to replace the time-varying filter with the corresponding steady-state version to reduce estimator complexity. On the other hand, the filter in (42) only iterates d-1 steps; hence we can implement the estimator in finite horizon on the basis of the steady-state filter (40).


## 3.1. Finite Horizon Estimation

We note thaty1(k) is an additional measurement of the state x(k-d) which is obtained at time instant k with time delay d, so the measurement y(k) consists of time delay when k≥d. Apparently, the linear space L{y(s)}s=0k contains the same information as L{Y1(s)}s=0k-d,{Y0(s)}s=k-d+1k, where the new observations Y1(s) and Y0(s) are provided as follows:(7)Y1s=y0sy1s+d,0≤s≤k-d(8)Y0s=y0s,k-d<s≤kObviously, Y1(s) and Y0(s) satisfy(9)Y1s=Bxs+CUs+Vs(10)Y0s=B0xs+C1u0s+V0swhere(11)B=B0B1,C=C000C1,Us=u0su1s+d,Vs=v0sv1s+d,V0s=v0sIt is obvious that the new measurements Y0(s) and Y1(s) are delay-free and the associated measurement noises V0(s) and V(s) are white noises with zero mean and covariance matrices RV0(s)=R0, RV(s)=R=diag⁡{R0,R1}.Theorem 5.
Whenk≥d, based on linear minimum variance unbiased criterion, we produce a recursive state estimator decoupling with the disturbance for system (1), (7), and (8) in the Kalman-like form:(12)x^s+1,0=Ax^s,0+K0sY0s-B0x^s,0,k-d<s<kx^k-d+1,0=x^k-d+1,1where(13)K0s=G0s-Λ0TsC0TD0-1sD0s=B0P0sB0T+R0G0s=AP0sB0TΛ0Ts=C0TD0-1sC0-1C0TD0-1sG0TsThe prediction error covariance matrix P0(s+1) is computed by(14)P0s+1=K0sD0sK0Ts-G0sK0Ts-K0sG0Ts+AP0sAT+Q(15)P0k-d+1=P1k-d+1As for x^(k-d+1,1) and P1(k-d+1), they are obtained by(16)x^s+1,1=Ax^s,1+K1sY1s-Bx^s,1,0≤s≤k-d(17)x^0,1=μ0,where(18)K1s=G1s-Λ1TsCTD1-1s(19)D1s=BP1sBT+R(20)G1s=AP1sBT(21)Λ1Ts=CTD1-1sC-1CTD1-1sG1TsThe prediction error covariance matrix P1(s+1) is computed by(22)P1s+1=K1sD1sK1Ts-G1sK1Ts-K1sG1Ts+AP1sAT+Q(23)P10=P0.Proof.
When0<s≤k-d, from (1) and (16), we have the prediction error equation as follows:(24)x~s+1,1=xs+1-x^s+1,1=Axs+ns-Ax^s,1-K1sY1s-Bx^s,1=A-K1sBx~s,1+ns-K1sVs-K1sCUsAssume x^(s,1) to be unbiased; in order to guarantee that x^(s+1,1) be an unbiased estimate of x(s+1), we must have Ex~(s+1,1)=0, then we have(25)K1sC=0Substituting (25) into (24), (24) can be rewritten as(26)x~s+1,1=A-K1sBx~s,1+ns-K1sVsFrom (26), the prediction error covariance matrix P1(s+1) is computed by(27)P1s+1=Ex~s+1,1x~Ts+1,1=A-K1sBP1sA-K1sBT+Q+K1sRK1Ts=K1sD1sK1Ts-G1sK1Ts-K1sG1Ts+AP1sAT+Qwhere D1(s) and G1(s) are given by (19) and (20). In order to minimize the estimation error variance (27) under the constraint (25), we introduce an auxiliary equation(28)Js=trP1s+1+trK1sCΛs+trΛTsCTK1Tswhere Λ(s) is the Lagrange multipliers. Taking the derivatives with respect to K1(s) equal to zero yields(29)K1sD1s+ΛsCT-G1s=0Combining (25) and (29) gives the matrix equation(30)D1sCCT0K1TsΛs=G1Ts0
Equation (30) has a unique solution if and only if the coefficients is nonsingular. Due to Assumption 3 and the fact that D1(s) is nonsingular, the coefficient matrix of (30) is nonsingular. Obviously, premultiplying left- and right-hand sides of (30) by the inverse of coefficient matrix yields (19) and (21). When k-d<s≤k, the proof is similar to the case 0≤s≤k-d, so we omit it here.Remark 6.
The system dimension will get more and more higher with the increase of time delay d; thus, computational complexity will get higher when using the classical state augmentation approach. But when we adopt the reorganized innovation analysis approach, we only need to solve two Riccati equations which have the same order of the state equation, so it can greatly reduce the computational complexity compared with the classical augmentation approach when time delay d is large.


## 3.2. Infinite Horizon Estimation

In this subsection, we consider the steady unbiased estimator design. First, let us present the lemma below as the initial step for the convergence analysis.Lemma 7.
If there existK∈Rn×(m0+m1) such that(31)KC=0,λiA-KB<1for i=1,2,…,n, then the sequence {P1(s+1)} abided by (22) is bounded for any s given any initial condition 0≤{P1(0)}≤∞.Proof.
Let us consider a suboptimal unbiased filter as follows:(32)x^subs+1,1=Ax^s,1+KsubY1s-Bx^s,1,0≤s≤k-dwhere KsubC=0,λi[A-KsubB]<1. Then the following state estimation error x~sub(s+1,1)=x(s+1)-x^sub(s+1,1) is given by(33)x~subs+1,1=A-KsubBx~subs,1+ns-KsubVswith the associated covariance matrix P1sub(s+1) being given by(34)P1subs+1=A-KsubBP1subsA-KsubBT+Q+KsubRKsubTThus P1sub(s+1) is bounded for any nonnegative initial condition due to the fact that λi[A-KsubB]<1. Comparing the above suboptimal estimator to the designed optimal estimator, the optimality tell us that P1(s+1)≤P1sub(s+1) for the same initial value. This proves the boundedness of P1sub(s+1).Theorem 8.
If there existK∈Rn×(m0+m1) such that KC=0,λi[A-KB]<1, and A,Q1/2 is stabilizable, then {P1(s+1)} abided by (22) converges to a unique fixed point P1 for any initial condition, where P1 is computed by(35)P1=K1D1K1T-G1K1T-K1G1T+AP1AT+Qand(36)K1=G1-Λ1TCTD1-1(37)D1=BP1BT+R(38)G1=AP1BT(39)Λ1T=CTD1-1C-1CTD1-1G1T.When 0<s≤k-d, we have the state predictor as follows:(40)x^s+1,1=Ax^s,1+K1Y1s-Bx^s,1(41)x^0,1=μ0When k-d<s≤k, the state predictor is as follows:(42)x^s+1,0=Ax^s,0+K0sY0s-B0x^s,0(43)x^k-d+1,0=x^k-d+1,1,where(44)K0s=D0-1sG0Ts-C0Λ0TsTD0s=B0P0sB0T+R0G0s=AP0sB0TΛ0Ts=C0TD0-1sC0-1C0TD0-1sG0TsThe prediction error covariance matrix P0(s+1) is computed by(45)P0s+1=K0sD0sK0Ts-G0sK0Ts-K0sG0Ts+AP0sAT+Q(46)P0k-d+1=P1.Proof.
By simple calculation, the Riccati equation in (22) can be formulated as(47)P1s+1=A-K1sBP1sA-K1sBT+Q+K1sRK1Ts.Then according to [35, 37] and Lemma 7, the convergence analysis can be readily obtained and we omitted here.Remark 9.
Under the condition thatKC=0,λi[A-KB]<1, and A,Q1/2 is stabilizable, we have obtained the steady-state filter (40). This is important when one desires to replace the time-varying filter with the corresponding steady-state version to reduce estimator complexity. On the other hand, the filter in (42) only iterates d-1 steps; hence we can implement the estimator in finite horizon on the basis of the steady-state filter (40).


## 4. Numerical Example

In this section, we present a numerical example to manifest the proposed approach about linear optimal estimation. Consider the linear discrete-time system with measurement delay and unknown observation disturbance(48)xk+1=0.780.400.300.60xk+nk,y0k=1221xk+11u0k+v0k,y1k=2112xk-10+1001u1k+v1k,where x(0)=-22, x^(0,1)=00, Q=1001, u0(k)=0.2, u1(k)=0.250.3,R0=R1=I2×2, n(k), v0(k), and v1(k) are whites noises with zero mean and covariances Q, R1 and R2, respectively.According to Theorem5, the simulation results are given in Figures 1 and 2, respectively. From the simulation results, we observe that the estimator can track the true state x(k) well, which proves that our proposed approach in this paper is effective. According to Theorem 8, we obtain the following steady estimator:(49)x^s+1,1=Ax^s,1+K1Y1s-Bx^s,1x^0,1=μ0,where(50)P=10.75467.04477.04476.1737.Figure 1
The first state componentx1(k) and the filter x^1(k) for two channels.Figure 2
The second state componentx2(k) and the predictor x^2(k) for two channels.


## 5. Conclusion

In this paper, we have proposed a linear minimum unbiased predictor for discrete-time systems with measurement delay and unknown measurement disturbance. Firstly, we have used the reorganized innovation analysis approach to deal with the measurement delay. In this way, one has avoided the giant computation brought by the augmentation approach or partial difference Riccati equation approach. Then based on the linear unbiased minimum variance criterion, we have calculated the minimum variance unbiased predictor, which is designed by calculating two Riccati equations with the same dimension as the state model. The future study direction is to consider the linear unbiased estimation for discrete-time systems with packet dropping and unknown disturbance, where the unknown disturbance appears in both the state equation and the measurement equation.



---

*Source: 2831561-2018-09-25.xml*