# A New Estimation Algorithm from Measurements with Multiple-Step Random Delays and Packet Dropouts

Authors: R. Caballero-Águila; A. Hermoso-Carazo; J. Linares-Pérez
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

The least-squares linear estimation problem using covariance
information is addressed in discrete-time linear stochastic
systems with bounded random observation delays which can lead to
bounded  packet dropouts. A recursive algorithm, including  the
computation of predictor, filter, and fixed-point smoother, is
obtained  by an innovation approach. The random delays  are
modeled by introducing some Bernoulli random variables with known
distributions in the system description. The derivation of the
proposed estimation algorithm does not require full knowledge of
the state-space model generating the signal to be estimated, but
only the delay probabilities and the covariance functions of the
processes involved in the observation equation.


---

## Body

## 1. Introduction

Originally, signal estimation problems were addressed under the assumption that the sensor data are transmitted over perfect communication channels, thus being received either instantaneously or with a known deterministic delay at the data processing center which provides the estimation. However, the use of communication networks for transmitting measured data motivates the need of considering eventual transmission delays and/or possible packet losses, due to numerous causes, such as random failures in the transmission mechanism, accidental loss of some measurements, or data inaccessibility at certain times. Often, these network uncertainties are random in nature and, hence, an appropriate model for such situations consists of describing the sensor delay or multiple packet dropout by a stochastic process, whose statistical properties are included in the system description. Therefore, estimation problems with bounded random delays in the observations and packet dropouts are challenging problems in networked control systems and have attracted much research interest.Assuming that the state-space model of the signal to be estimated is known, many results have been reported on systems with random delays and packet dropouts. For example, Ray et al. [1] proposed a recursive linear filtering algorithm which modifies the conventional one to fit situations where the arrival of sensor data at the controller terminal may be randomly delayed. In [2], the state estimation problem for a model involving randomly varying bounded sensor delays is treated by reformulating it as an estimation problem in systems with stochastic parameters. More recently, Matveev and Savkin [3] have proposed a recursive minimum variance state estimator in linear discrete-time partially observed systems perturbed by white noises, when the observations are transmitted via communication channels with random transmission times and various signal measurements may incur independent delays. Wang et al. [4] have designed a robust linear filter for linear uncertain discrete-time stochastic systems with randomly varying sensor delay. Su and Lu [5] have designed an extended Kalman filtering algorithm which provides optimal estimates of interconnected network states for systems in which some or all measurements are delayed. Hermoso-Carazo and Linares-Pérez [6] have proposed a filtering algorithm based on the unscented transformation to estimate the state of a nonlinear system from randomly delayed measurements. All of the above mentioned papers on signal estimation from randomly delayed observations assume that the random delay does not exceed one sampling time. Recently, Hounkpevi and Yaz [7] have considered the estimation problem from observations coming from multiple sensors with different one-step delay rates, indicating that the model considered could be generalized to the case of multiple-sample delay.For the packet-dropout problem, many recent papers can be mentioned. For example, the optimalH2 filtering problem in networked control systems with random delays of one sampling period, multiple packet dropouts, and uncertain observations are studied under a unified framework in [8], and the optimal H∞ filtering problem in networked control systems with multiple packet dropouts is addressed in [9]. The optimal linear minimum variance estimation problem for linear discrete-time stochastic systems with possible infinite packet dropouts and the steady-state estimators are studied in [10]. Also, for these systems with possible infinite packet dropouts the optimal full-order and reduced-order estimators in the linear minimum variance sense are obtained in [11]. In practice, the number of consecutive packet dropouts cannot be infinite but bounded by a finite number and, consequently, the above mentioned papers may lead to conservative results. A novel model to describe the case when the number of consecutive packet dropouts is limited by a known upper bound is considered in [12] and the linear filtering problem is addressed assuming that the filter has a similar recursive structure to that of the Kalman one; hence, for assuming the filter to have a fixed structure, it is only a suboptimal filter. This study is completed in [13] by considering the optimal estimation problem (including filtering, prediction, and smoothing) in the linear least-mean-square sense. In the above mentioned papers only packet dropouts are considered; recently, Sun [14] has addressed the linear filtering, prediction, and smoothing problems in discrete-time linear systems with finite random measurement delays and packet dropouts, assuming that the delay and the number of consecutive dropouts do not exceed a known upper bound.On the other hand, when the state-space model of the signal to be estimated is not available, it is necessary to use alternative information, for example, about the covariance functions of the processes involved in the observation equation. In this context, the least-squares (LS) linear and second-order polynomial estimation problems from randomly delayed observations based on covariance information have been addressed in [15, 16], respectively, under the assumption that the Bernoulli random variables modeling the delays are independent; also, linear and polynomial estimation algorithms from observations featuring correlated random delays have been proposed in [17, 18], respectively. Recently, the LS linear filtering problem of discrete-time signals using one- and two-step randomly delayed observations coming from multiple sensors with different delay rates has been studied in [19] using covariance information.In this paper, the LS linear estimation problem in systems with bounded random measurement delays and packet dropouts is addressed. The proposed estimators only depend on the delay probabilities at each sampling time, but do not need to know if a particular measurement is delayed or updated. Moreover, the estimation algorithm is derived using only covariance information. Consequently, considering the case of sensors with the same delay characteristics, the current study generalizes the results in [19] to the case of measurements with bounded multiple-step random delays and packet dropouts.The paper is organized as follows. In Section2 the observation model considered and the hypotheses on the signal and noise processes are presented. The LS linear estimation problem is formulated in Section 3, where also the innovation technique used to address such problem is described. The LS linear estimation algorithm is derived in Section 4 which includes recursive formulas to obtain the estimation-error covariance matrices; these matrices provide a global measure of the LS estimators accuracy. Finally, in Section 5, a numerical simulation example is presented to show the effectiveness of the estimation algorithm proposed in the current paper.


## 2. Observation Model

In networked systems, such as telephone networks, cable TV networks, cellular networks or the Internet, among others, the system output is measured at every sampling time and the measurement is transmitted to a data processing center producing the signal estimation. In the transmission, delays and packet dropouts are unavoidable due to eventual communication failures and, to reduce the effect of such delays and packet dropouts without overloading the network traffic, each sensor measurement is transmitted for several consecutive sampling times, but only one measured output is processed for the estimation at each sampling time.In this paper it is assumed that the largest delay is upper bounded (the bound is a known constant denoted byD); hence, the current study is a generalization of that performed in [19] for the estimation problem from one- or two-step randomly delayed observations. Note that this observation model includes bounded packet dropouts since, if a measurement is not processed after D sampling times, such measurement is lost in transmission.In this section, the observation model with bounded random measurement delays and packet dropouts as well as the assumptions about the signal and noises involved, are presented.Consider a signal vector,zk, whose measured output at the sampling time k, denoted by ỹk, is perturbed by an additive noise vector vk; that is,(2.1)ỹk=zk+vk,k≥1.The measured outputỹk is transmitted during the sampling times k,k+1,…, and k+D but, at each sampling time k, only one of the measurements ỹk-D,…,ỹk is processed; consequently, at time k>D, the processed measurement can be either delayed by d=1,…,D sampling periods with a known probability pk(d), or updated with probability pk(0). At the initial time k=1, the measured output ỹ1 is always available (p1(0)=1) and, hence, the processed measurement is equal to the real measurement y1=ỹ1. At any time k≤D, the processed measurement only can be delayed by d=1,…,k-1 sampling periods, since only ỹ1,…,ỹk are available. Also, it is assumed that the delays at different times are independent.Therefore, the following model for the processed measurements to estimate the signal is considered:(2.2)yk=∑d=0min(k-1,D)γk(d)ỹk-d,k≥1,
where, for d=0,1,…,D,{γk(d);k>d} denote sequences of mutually independent Bernoulli random variables with P[γk(d)=1]=pk(d) and ∑d=0min(k-1,D)γk(d)=1.This fact guarantees that ifγk(0)=1, γk(d)=0 for all d≠0, which means that, with probability pk(0), the kth measurement is received and processed on time. If γk(0)=0, then there exists one and only one d=1,…,D such that γk(d)=1, which means that the measurement yk is delayed by d sampling periods with probability pk(d). Note that the measured output ỹk at any time k can be received on time, delayed, or lost in transmission. Also note that some measured output can be rereceived, since the output at each instant is transmitted for D+1 consecutive times (please, see Figure 1 in Section 5 for specific results when D=3).Figure 1
Processed observations forD=3 when q1=q2=q3=0.5.The signal estimation problem is addressed based on the following assumptions.Assumption 1.
Then-dimensional signal process {zk;k≥1} has zero mean and autocovariance function Kk,s=E[zkzsT]=AkBsT,s≤k, where A and B are known matrix functions.Assumption 2.
The noise process{vk;k≥1} is a zero-mean white sequence with known autocovariance function E[vkvkT]=Rk,for all k≥1.Assumption 3.
Ford=0,1,…,D,{γk(d);k>d} are sequences of independent Bernoulli random variables with P[γk(d)=1]=pk(d) and ∑d=0min(k-1,D)γk(d)=1.Moreover,{γk(d),0≤d≤min(k-1,D)} is independent of {γj(d),0≤d≤min(j-1,D)} for all k≠j.Assumption 4.
For eachd=0,1,…,D, the processes {zk;k≥1},{vk;k≥1}, and {γk(d);k>d} are mutually independent.Remark 2.1.
The estimation problem is addressed under the assumption that the evolution model of the signal need not be available and using only information about the covariance functions of the processes involved in the observation equation. Note that, although a state-space model can be generated from covariances, when only this kind of information is available, it is preferable to address the estimation problem directly using covariances, thus obviating the need for previous identification of the state-space model.Remark 2.2.
Although Assumption1 might seem restrictive, it covers many practical situations; for example, when the system matrix Φ in the state-space model of a stationary signal is available, the signal autocovariance function can be expressed as E[zkzsT]=Φk-sE[zszsT],s≤k, and Assumption 1 is clearly satisfied, taking Ak=Φk and Bs=E[zszsT](Φ-s)T. Also, processes with finite-dimensional, possibly time-variant, state-space models have semiseparable covariance functions, E[zkzsT]=∑i=1rakibsi,s≤k (see [20]), and this structure is a particular case of that assumed, just taking Ak=(ak1,ak2,…,akr) and Bs=(bs1,bs2,…,bsr). Consequently, this structural assumption on the signal autocovariance function covers both stationary and nonstationary signals.


## 3. Problem Statement

Given the observation model (2.1)-(2.2) with random delays and packet dropouts, our purpose is to find the least-squares (LS) linear estimator, ẑk/L, of the signal zk based on the observations {y1,…,yL}. Specifically, our aim is to obtain the filter (L=k) and the fixed-point smoother (k fixed and L>k) from recursive formulas.The estimatorẑk/L is the orthogonal projection of the vector zk onto ℒ(y1,…,yL), the n-dimensional linear space spanned by {y1,…,yL}. Since the observations are generally nonorthogonal vectors, we use an innovation approach, based on an orthogonalization procedure by means of which the observation process {yk;k≥1} is transformed into an equivalent one (innovation process) of orthogonal vectors {νk;k≥1}, equivalent in the sense that each set {ν1,…,νL} spans the same linear subspace as {y1,…,yL}.Since the innovations constitute a white process, this methodology allows us to find the orthogonal projection of the vectorzk onto {ν1,…,νL} by separately projecting onto each of the previous orthogonal vectors; that is, the replacement of the observation process by the innovation one leads to the following expression of the signal estimators:(3.1)ẑk/L=∑i=1LSk,izΠi-1νi,
where Sk,iz=E[zkνiT] and Πi=E[νiνiT].Hence, to obtain the signal estimators it is necessary to find previously an explicit formula for the innovations and their covariance matrices.Innovationνk and Covariance Matrix Πk
The innovation at timek is defined as νk=yk-ŷk/k-1 where ŷk/k-1 is the one-stage linear predictor of yk. Since the LS linear estimator ŷk/k-1 is the projection of yk onto ℒ(ν1,…,νk-1), from (2.1) and (2.2) and model assumptions, it is clear that
(3.2)ŷk/k-1=∑d=0min(k-1,D)pk(d)(ẑk-d/k-1+v̂k-d/k-1),k≥2,ŷ1/0=0,
where v̂i/j denotes the orthogonal projection of the vector vi onto ℒ(ν1,…,νj).Similar that to (3.1), by denoting Sk,iv=E[vkνiT], the noise estimators are expressed as follows:(3.3)v̂k/L=∑i=1LSk,ivΠi-1νi.
From the model assumptions, the noise vk is independent of {y1,…,yk-1} and hence, for i<k, we have that Sk,iv=0; thus, the one-stage predictor of the noise is v̂k/k-1=0. Consequently, the innovation is given by
(3.4)νk=yk-pk(0)ẑk/k-1-∑d=1min(k-1,D)pk(d)(ẑk-d/k-1+v̂k-d/k-1),k≥2,ν1=y1,
and it is necessary to obtain the signal predictor ẑk/k-1 and the estimators ẑk-d/k-1 and v̂k-d/k-1 for d=1,…,min(k-1,D), that is, the filters and the smoothers of the signal and noise, respectively.Remark 3.1.
From the model assumptions, it is clear thatSk,kv=pk(0)Rk and, consequently, the filter of the noise is v̂k/k=pk(0)RkΠk-1νk.To obtain the covariance matrixΠk, from (2.1), (2.2), and (3.2) and taking into account that v̂k/k-1=0, the innovation νk is expressed as follows:(3.5)νk=∑d=0min(k-1,D)(γk(d)-pk(d))zk-d+∑d=1min(k-1,D)(γk(d)-pk(d))vk-d+γk(0)vk+∑d=0min(k-1,D)pk(d)(zk-d-ẑk-d/k-1)+∑d=1min(k-1,D)pk(d)(vk-d-v̂k-d/k-1).Using this expression and the model assumptions, we have(3.6)Πk=∑d,d′=0min(k-1,D)Cov(γk(d),γk(d′))Kk-d,k-d′+∑d,d′=0min(k-1,D)pk(d)pk(d′)Pk-d,k-d′/k-1zz+∑d=1min(k-1,D)pk(d)(1-pk(d))Rk-d+pk(0)Rk+∑d,d′=1min(k-1,D)pk(d)pk(d′)Pk-d,k-d′/k-1vv+∑d=0,d′=1min(k-1,D)pk(d)pk(d′)Pk-d,k-d′/k-1zv+∑d=1,d′=0min(k-1,D)pk(d)pk(d′)Pk-d,k-d′/k-1vz,k≥2,Π1=A1B1T+R1,
where Pl,m/Nzz=E[(zl-ẑl/N)(zm-ẑm/N)] and Pl,m/Nvv=E[(vl-v̂l/N)(vm-v̂m/N)] denote the estimation-error covariance matrices of the state and noise, respectively, and Pl,m/Nzv=E[(zl-ẑl/N)(vm-v̂m/N)] the estimation-error cross-covariance matrix between the state and noise. Hence, the filtering, and smoothing error covariance and cross-covariance matrices of the state and noise must be calculated.


## 4. Recursive Estimation Algorithm

In this section we obtain the state predictor (Section4.1), the filtering, and smoothing estimators of the signal and noise (Section 4.2), and the filtering, and smoothing error covariance and cross-covariance matrices of the state and noise (Section 4.3), which together with (3.4) and (3.6) constitute the proposed recursive estimation algorithm.


### 4.1. Predictor of the Signal

From (3.1), to determine the signal predictor, ẑk/k-1, it is necessary to calculate the coefficients
(4.1)Sk,iz=E[zkνiT]=E[zkyiT]-E[zkŷi/i-1T],i<k,
so E[zkyiT] and E[zkŷi/i-1T] must be calculated for i<k.(i) On the one hand, using (2.1)-(2.2) and the model assumptions, it is clear that
(4.2)E[zkyiT]=∑d=0min(i-1,D)pi(d)E[zkzi-dT]=∑d=0min(i-1,D)pi(d)AkBi-dT=AkGBiT,1≤i≤k,
where GBi=∑d=0min(i-1,D)pi(d)Bi-d,i≥1.(ii) On the other hand, from (3.2) and taking into account that v̂i/i-1=0, we have that
(4.3)E[zkŷi/i-1T]=∑d=0min(i-1,D)pi(d)E[zkẑi-d/i-1T]+∑d=1min(i-1,D)pi(d)E[zkv̂i-d/i-1T],2≤i≤k,
so E[zkẑi-d/i-1T] and E[zkv̂i-d/i-1T] must be calculated for i=1,…,k.(a)
Using (3.1) for ẑi-d/i-1T,
(4.4)E[zkẑi-d/i-1T]=E[zk{∑j=1i-1Si-d,jzΠj-1νj}T]=∑j=1i-1Sk,jzΠj-1Si-d,jzT,2≤i≤k.(b)
Using (3.3) for v̂i-d/i-1T, and since Si-d,jv=0, for j<i-d, we obtain
(4.5)E[zkv̂i-d/i-1T]=E[zk{∑j=i-di-1Si-d,jvΠj-1νj}T]=∑j=i-di-1Sk,jzΠj-1Si-d,jvT,2≤i≤k.Substituting into (4.1) the expectations calculated in (i) and (ii), we obtain (4.6)Sk,iz=AkGBiT-∑d=0min(i-1,D)pi(d){∑j=1i-1Sk,jzΠj-1Si-d,jzT}-∑d=1min(i-1,D)pi(d){∑j=i-di-1Sk,jzΠj-1Si-d,jvT},2≤i≤k,Sk,1z=AkGB1T.
If we now introduce a function Ji satisfying that(4.7)Ji=GBiT-∑d=0min(i-1,D)pi(d){∑j=1i-1JjΠj-1Si-d,jzT}-∑d=1min(i-1,D)pi(d){∑j=i-di-1JjΠj-1Si-d,jvT},i≥2,J1=GB1T,
we conclude that(4.8)Sk,iz=AkJi,1≤i≤k.Substituting now (4.8) into (3.1) for L=k-1, the following expression for the state predictor is obtained: (4.9)ẑk/k-1=AkOk-1,k≥1,
where the vectors Ok are defined by (4.10)Ok=∑i=1kJiΠi-1νi,k≥1,O0=0,
thus satisfying the recursive relation (4.11)Ok=Ok-1+JkΠk-1νk,k≥1,O0=0.Hence, an expression forJk must be derived.Remark 4.1.
By substituting (4.8) into (3.1), it is also clear that the l-stage predictors are given by ẑk+l/k=Ak+lOk,k≥1,l≥1.Expression forJk
First, expression (4.7) is just rewritten for i=k:
(4.12)Jk=GBkT-∑d=0min(k-1,D)pk(d){∑j=1k-1JjΠj-1Sk-d,jzT}-∑d=1min(k-1,D)pk(d){∑j=k-dk-1JjΠj-1Sk-d,jvT},k≥2.
For k=2, using that S2,1z=A2J1,S1,1z=A1J1, and S1,1v=p1(0)R1, and denoting r1=J1Π1-1J1T, it is immediately clear that
(4.13)J2=GB2T-p2(0)r1A2T-p2(1){r1A1T+p1(0)J1Π1-1R1}.
For k≥3, we examine separately the sums that appear in Jk as follows.
(i) SinceSm,jz=AmJj, for j≤m, by denoting that  rk=∑j=1kJjΠj-1JjT, we have that∑d=0min(k-1,D)pk(d){∑j=1k-1JjΠj-1Sk-d,jzT}=pk(0)rk-1AkT+pk(1)rk-1Ak-1T+∑d=2min(k-1,D)pk(d)rk-dAk-dT+∑d=2min(k-1,D)pk(d)∑i=1d-1Jk-d+iΠk-d+i-1Sk-d,k-d+izT,k≥3.
(ii) SinceSm,mv=pm(0)Rm,m≥1, we have that
(4.15)∑d=1min(k-1,D)pk(d){∑j=k-dk-1JjΠj-1Sk-d,jvT}=pk(1)Jk-1Πk-1-1pk-1(0)Rk-1+∑d=2min(k-1,D)pk(d){Jk-dΠk-d-1pk-d(0)Rk-d+∑i=1d-1Jk-d+iΠk-d+i-1Sk-d,k-d+ivT},k≥3.By substituting the above two sums in Jk, it is deduced that
(4.16)Jk=GBkT-pk(0)rk-1AkT-∑d=1min(k-1,D)pk(d){rk-dAk-dT+pk-d(0)Jk-dΠk-d-1Rk-d}-∑d=2min(k-1,D)pk(d)∑i=1d-1Jk-d+iΠk-d+i-1{Sk-d,k-d+izT+Sk-d,k-d+ivT},k≥3.
Finally, note that Sk-d,k-d+iz and Sk-d,k-d+iv are the smoothing gain matrices (which will be obtained in the next section), and the matrix rk is recursively obtained by
(4.17)rk=rk-1+JkΠk-1JkT,k≥1,r0=0.


### 4.2. Filtering and Smoothing Estimators of the Signal and Noise

Clearly, in view of (3.1) and (3.3), the filters and fixed-point smoothers of the signal and noise are obtained by the following recursive expressions:(4.18)ẑk/L=ẑk/L-1+Sk,LzΠL-1νL,L≥k,k≥1,v̂k/L=v̂k/L-1+Sk,LvΠL-1νL,L≥k,k≥1,
with initial conditions given by the one-stage predictors ẑk/k-1=AkOk-1 and v̂k/k-1=0, respectively. Hence, the gain matrices Sk,Lz=E[zkνLT] and Sk,Lv=E[vkνLT] must be calculated for L≥k.Gain Matrices
Using (3.5) for k=L and the model assumptions, it is derived that
(4.19)Sk,Lz=∑d=0min(L-1,D)pL(d)E[zk(zL-d-ẑL-d/L-1)T]+∑d=1min(L-1,D)pL(d)E[zk(vL-d-v̂L-d/L-1)T],
and, since the estimation-errors are orthogonal to the estimators, it is concluded that
(4.20)Sk,Lz=∑d=0min(L-1,D)pL(d)Pk,L-d/L-1zz+∑d=1min(L-1,D)pL(d)Pk,L-d/L-1zv,
where Pk,L-d/L-1zz and Pk,L-d/L-1zv are the covariance and cross-covariance matrices of the errors z̃k/L-1=zk-ẑk/L-1 with z̃L-d/L-1 and ṽL-d/L-1=vL-d-v̂L-d/L-1, respectively.A similar reasoning leads to the following expression:(4.21)Sk,Lv=∑d=0min(L-1,D)pL(d)Pk,L-d/L-1vz+∑d=1min(L-1,D)pL(d)Pk,L-d/L-1vv,
where Pk,L-d/L-1vv and Pk,L-d/L-1vz are the covariance and cross-covariance matrices of the errors ṽk/L-1 with ṽL-d/L-1 and z̃L-d/L-1, respectively.Remark 4.2.
Note that, from (4.8), the gain matrix of the signal filter is Sk,kz=AkJk and, from (3.1), the filter is ẑk/k=AkOk.


### 4.3. Error Covariance and Cross-Covariance Matrices

From the recursive relations (4.18), the estimation-errors admit the following expressions: (4.22)z̃k/L=z̃k/L-1-Sk,LzΠL-1νL,L≥k,k≥1,ṽk/L=ṽk/L-1-Sk,LvΠL-1νL,L≥k,k≥1.
Now, since the estimation-errors are orthogonal to the estimators, we have that E[z̃k/L-1νLT]=E[zkνLT]=Sk,Lz and E[ṽk/L-1νLT]=E[vkνLT]=Sk,Lv, thus deducing the following recursive expressions for the error covariance and cross-covariance matrices:(4.23)Pk,m/Lzz=Pk,m/L-1zz-Sk,LzΠL-1Sm,LzT,L≥k,m,k,m≥1,Pk,m/Lvv=Pk,m/L-1vv-Sk,LvΠL-1Sm,LvT,L≥k,m,k,m≥1,Pk,m/Lzv=Pk,m/L-1zv-Sk,LzΠL-1Sm,LvT,L≥k,m,k,m≥1.The initial conditions of these equations, the prediction error covariance and cross-covariance matrices, are obtained as follows:(4.24)Pk+i,k/k-1zz=Ak+i[BkT-rk-1AkT],i≥0,k≥1,Pk+i,k/k-1vv=Rkδi,0,i≥0,k≥1,Pk+i,k/k-1zv=Pk,k+i/k-1zv=0,i≥0,k≥1,
where δ denotes the Kronecker delta function.In fact, by writingPk+i,k/k-1zz=E[zk+izkT]-E[ẑk+i/k-1ẑk/k-1T] and using that ẑk+i/k-1=Ak+iOk-1 and E[Ok-1Ok-1T]=rk-1, the expression for Pk+i,k/k-1zz is immediate from Assumption 1. Analogously, using that v̂k+i/k-1=0, for i≥0, and Assumptions 2 and 4, the expressions for Pk+i,k/k-1vv, Pk+i,k/k-1zv, and Pk,k+i/k-1zv are derived.


### 4.4. Computational Procedure

At the sampling timek, once the (k-1)th iteration is finished and the new observation yk is available, the proposed estimation algorithm operates as follows.(i)
Compute the innovationνk and its covariance matrix Πk by (3.4) and (3.6), respectively.(ii)
ComputeJk by (4.16) and, from it, the filtering gain matrix Sk,kz=AkJk.(iii)
Compute the filterẑk/k=ẑk/k-1+Sk,kzΠk-1νk and the filtering error covariance matrix Pk,k/kzz=Pk,k/k-1zz-Sk,kzΠk-1Sk,kzT.(iv)
To implement the above steps at timek+1, we need to
(a)
computeOk=Ok-1+JkΠk-1νk and, from it, the signal predictor ẑk+1/k=Ak+1Ok,(b)
compute the noise filterv̂k/k=pk(0)RkΠk-1νk,(c)
ford=1,…,D-1, compute the smoothing gain matrices Sk-d,kz and Sk-d,kv from (4.20) and (4.21), respectively, and the smoothers ẑk-d/k=ẑk-d/k-1+Sk-d,kzΠk-1νk and v̂k-d/k=v̂k-d/k-1+Sk-d,kvΠk-1νk.Then, calculate the innovationνk+1. To obtain its covariance matrix, Πk+1, compute the error covariance and cross-covariance matrices of the previous estimators of the signal and noise (predictors, filters, and smoothers) using the formulas established in Section 4.3. Finally, compute rk=rk-1+JkΠk-1JkT to calculate Jk+1 and, from it, the filtering gain Sk+1,k+1z, which provides the filter ẑk+1/k+1 and the error covariance matrix Pk+1,k+1/k+1zz.


## 4.1. Predictor of the Signal

From (3.1), to determine the signal predictor, ẑk/k-1, it is necessary to calculate the coefficients
(4.1)Sk,iz=E[zkνiT]=E[zkyiT]-E[zkŷi/i-1T],i<k,
so E[zkyiT] and E[zkŷi/i-1T] must be calculated for i<k.(i) On the one hand, using (2.1)-(2.2) and the model assumptions, it is clear that
(4.2)E[zkyiT]=∑d=0min(i-1,D)pi(d)E[zkzi-dT]=∑d=0min(i-1,D)pi(d)AkBi-dT=AkGBiT,1≤i≤k,
where GBi=∑d=0min(i-1,D)pi(d)Bi-d,i≥1.(ii) On the other hand, from (3.2) and taking into account that v̂i/i-1=0, we have that
(4.3)E[zkŷi/i-1T]=∑d=0min(i-1,D)pi(d)E[zkẑi-d/i-1T]+∑d=1min(i-1,D)pi(d)E[zkv̂i-d/i-1T],2≤i≤k,
so E[zkẑi-d/i-1T] and E[zkv̂i-d/i-1T] must be calculated for i=1,…,k.(a)
Using (3.1) for ẑi-d/i-1T,
(4.4)E[zkẑi-d/i-1T]=E[zk{∑j=1i-1Si-d,jzΠj-1νj}T]=∑j=1i-1Sk,jzΠj-1Si-d,jzT,2≤i≤k.(b)
Using (3.3) for v̂i-d/i-1T, and since Si-d,jv=0, for j<i-d, we obtain
(4.5)E[zkv̂i-d/i-1T]=E[zk{∑j=i-di-1Si-d,jvΠj-1νj}T]=∑j=i-di-1Sk,jzΠj-1Si-d,jvT,2≤i≤k.Substituting into (4.1) the expectations calculated in (i) and (ii), we obtain (4.6)Sk,iz=AkGBiT-∑d=0min(i-1,D)pi(d){∑j=1i-1Sk,jzΠj-1Si-d,jzT}-∑d=1min(i-1,D)pi(d){∑j=i-di-1Sk,jzΠj-1Si-d,jvT},2≤i≤k,Sk,1z=AkGB1T.
If we now introduce a function Ji satisfying that(4.7)Ji=GBiT-∑d=0min(i-1,D)pi(d){∑j=1i-1JjΠj-1Si-d,jzT}-∑d=1min(i-1,D)pi(d){∑j=i-di-1JjΠj-1Si-d,jvT},i≥2,J1=GB1T,
we conclude that(4.8)Sk,iz=AkJi,1≤i≤k.Substituting now (4.8) into (3.1) for L=k-1, the following expression for the state predictor is obtained: (4.9)ẑk/k-1=AkOk-1,k≥1,
where the vectors Ok are defined by (4.10)Ok=∑i=1kJiΠi-1νi,k≥1,O0=0,
thus satisfying the recursive relation (4.11)Ok=Ok-1+JkΠk-1νk,k≥1,O0=0.Hence, an expression forJk must be derived.Remark 4.1.
By substituting (4.8) into (3.1), it is also clear that the l-stage predictors are given by ẑk+l/k=Ak+lOk,k≥1,l≥1.Expression forJk
First, expression (4.7) is just rewritten for i=k:
(4.12)Jk=GBkT-∑d=0min(k-1,D)pk(d){∑j=1k-1JjΠj-1Sk-d,jzT}-∑d=1min(k-1,D)pk(d){∑j=k-dk-1JjΠj-1Sk-d,jvT},k≥2.
For k=2, using that S2,1z=A2J1,S1,1z=A1J1, and S1,1v=p1(0)R1, and denoting r1=J1Π1-1J1T, it is immediately clear that
(4.13)J2=GB2T-p2(0)r1A2T-p2(1){r1A1T+p1(0)J1Π1-1R1}.
For k≥3, we examine separately the sums that appear in Jk as follows.
(i) SinceSm,jz=AmJj, for j≤m, by denoting that  rk=∑j=1kJjΠj-1JjT, we have that∑d=0min(k-1,D)pk(d){∑j=1k-1JjΠj-1Sk-d,jzT}=pk(0)rk-1AkT+pk(1)rk-1Ak-1T+∑d=2min(k-1,D)pk(d)rk-dAk-dT+∑d=2min(k-1,D)pk(d)∑i=1d-1Jk-d+iΠk-d+i-1Sk-d,k-d+izT,k≥3.
(ii) SinceSm,mv=pm(0)Rm,m≥1, we have that
(4.15)∑d=1min(k-1,D)pk(d){∑j=k-dk-1JjΠj-1Sk-d,jvT}=pk(1)Jk-1Πk-1-1pk-1(0)Rk-1+∑d=2min(k-1,D)pk(d){Jk-dΠk-d-1pk-d(0)Rk-d+∑i=1d-1Jk-d+iΠk-d+i-1Sk-d,k-d+ivT},k≥3.By substituting the above two sums in Jk, it is deduced that
(4.16)Jk=GBkT-pk(0)rk-1AkT-∑d=1min(k-1,D)pk(d){rk-dAk-dT+pk-d(0)Jk-dΠk-d-1Rk-d}-∑d=2min(k-1,D)pk(d)∑i=1d-1Jk-d+iΠk-d+i-1{Sk-d,k-d+izT+Sk-d,k-d+ivT},k≥3.
Finally, note that Sk-d,k-d+iz and Sk-d,k-d+iv are the smoothing gain matrices (which will be obtained in the next section), and the matrix rk is recursively obtained by
(4.17)rk=rk-1+JkΠk-1JkT,k≥1,r0=0.


## 4.2. Filtering and Smoothing Estimators of the Signal and Noise

Clearly, in view of (3.1) and (3.3), the filters and fixed-point smoothers of the signal and noise are obtained by the following recursive expressions:(4.18)ẑk/L=ẑk/L-1+Sk,LzΠL-1νL,L≥k,k≥1,v̂k/L=v̂k/L-1+Sk,LvΠL-1νL,L≥k,k≥1,
with initial conditions given by the one-stage predictors ẑk/k-1=AkOk-1 and v̂k/k-1=0, respectively. Hence, the gain matrices Sk,Lz=E[zkνLT] and Sk,Lv=E[vkνLT] must be calculated for L≥k.Gain Matrices
Using (3.5) for k=L and the model assumptions, it is derived that
(4.19)Sk,Lz=∑d=0min(L-1,D)pL(d)E[zk(zL-d-ẑL-d/L-1)T]+∑d=1min(L-1,D)pL(d)E[zk(vL-d-v̂L-d/L-1)T],
and, since the estimation-errors are orthogonal to the estimators, it is concluded that
(4.20)Sk,Lz=∑d=0min(L-1,D)pL(d)Pk,L-d/L-1zz+∑d=1min(L-1,D)pL(d)Pk,L-d/L-1zv,
where Pk,L-d/L-1zz and Pk,L-d/L-1zv are the covariance and cross-covariance matrices of the errors z̃k/L-1=zk-ẑk/L-1 with z̃L-d/L-1 and ṽL-d/L-1=vL-d-v̂L-d/L-1, respectively.A similar reasoning leads to the following expression:(4.21)Sk,Lv=∑d=0min(L-1,D)pL(d)Pk,L-d/L-1vz+∑d=1min(L-1,D)pL(d)Pk,L-d/L-1vv,
where Pk,L-d/L-1vv and Pk,L-d/L-1vz are the covariance and cross-covariance matrices of the errors ṽk/L-1 with ṽL-d/L-1 and z̃L-d/L-1, respectively.Remark 4.2.
Note that, from (4.8), the gain matrix of the signal filter is Sk,kz=AkJk and, from (3.1), the filter is ẑk/k=AkOk.


## 4.3. Error Covariance and Cross-Covariance Matrices

From the recursive relations (4.18), the estimation-errors admit the following expressions: (4.22)z̃k/L=z̃k/L-1-Sk,LzΠL-1νL,L≥k,k≥1,ṽk/L=ṽk/L-1-Sk,LvΠL-1νL,L≥k,k≥1.
Now, since the estimation-errors are orthogonal to the estimators, we have that E[z̃k/L-1νLT]=E[zkνLT]=Sk,Lz and E[ṽk/L-1νLT]=E[vkνLT]=Sk,Lv, thus deducing the following recursive expressions for the error covariance and cross-covariance matrices:(4.23)Pk,m/Lzz=Pk,m/L-1zz-Sk,LzΠL-1Sm,LzT,L≥k,m,k,m≥1,Pk,m/Lvv=Pk,m/L-1vv-Sk,LvΠL-1Sm,LvT,L≥k,m,k,m≥1,Pk,m/Lzv=Pk,m/L-1zv-Sk,LzΠL-1Sm,LvT,L≥k,m,k,m≥1.The initial conditions of these equations, the prediction error covariance and cross-covariance matrices, are obtained as follows:(4.24)Pk+i,k/k-1zz=Ak+i[BkT-rk-1AkT],i≥0,k≥1,Pk+i,k/k-1vv=Rkδi,0,i≥0,k≥1,Pk+i,k/k-1zv=Pk,k+i/k-1zv=0,i≥0,k≥1,
where δ denotes the Kronecker delta function.In fact, by writingPk+i,k/k-1zz=E[zk+izkT]-E[ẑk+i/k-1ẑk/k-1T] and using that ẑk+i/k-1=Ak+iOk-1 and E[Ok-1Ok-1T]=rk-1, the expression for Pk+i,k/k-1zz is immediate from Assumption 1. Analogously, using that v̂k+i/k-1=0, for i≥0, and Assumptions 2 and 4, the expressions for Pk+i,k/k-1vv, Pk+i,k/k-1zv, and Pk,k+i/k-1zv are derived.


## 4.4. Computational Procedure

At the sampling timek, once the (k-1)th iteration is finished and the new observation yk is available, the proposed estimation algorithm operates as follows.(i)
Compute the innovationνk and its covariance matrix Πk by (3.4) and (3.6), respectively.(ii)
ComputeJk by (4.16) and, from it, the filtering gain matrix Sk,kz=AkJk.(iii)
Compute the filterẑk/k=ẑk/k-1+Sk,kzΠk-1νk and the filtering error covariance matrix Pk,k/kzz=Pk,k/k-1zz-Sk,kzΠk-1Sk,kzT.(iv)
To implement the above steps at timek+1, we need to
(a)
computeOk=Ok-1+JkΠk-1νk and, from it, the signal predictor ẑk+1/k=Ak+1Ok,(b)
compute the noise filterv̂k/k=pk(0)RkΠk-1νk,(c)
ford=1,…,D-1, compute the smoothing gain matrices Sk-d,kz and Sk-d,kv from (4.20) and (4.21), respectively, and the smoothers ẑk-d/k=ẑk-d/k-1+Sk-d,kzΠk-1νk and v̂k-d/k=v̂k-d/k-1+Sk-d,kvΠk-1νk.Then, calculate the innovationνk+1. To obtain its covariance matrix, Πk+1, compute the error covariance and cross-covariance matrices of the previous estimators of the signal and noise (predictors, filters, and smoothers) using the formulas established in Section 4.3. Finally, compute rk=rk-1+JkΠk-1JkT to calculate Jk+1 and, from it, the filtering gain Sk+1,k+1z, which provides the filter ẑk+1/k+1 and the error covariance matrix Pk+1,k+1/k+1zz.


## 5. Computer Simulation Results

In this section, the application of the proposed signal estimation algorithm is illustrated by a simulation example. Consider a zero-mean scalar signal{zk;k≥1} with autocovariance function given by E[zkzsT]=1.025641×0.95k-s,s≤k, which is factorizable according to Assumption 1 of the model just taking Ak=1.025641×0.95k and Bk=0.95-k. For the simulation, the signal is supposed to be generated by an autoregressive model, zk+1=0.95zk+wk, where {wk;k≥1} is a zero-mean white Gaussian noise with Var[wk]=0.1, for all k.Using the proposed filtering and fixed-point smoothing algorithms, we have estimated the signal from bounded random measurement delays and packet dropouts, assuming that the largest delay and the maximum number of successive dropouts isD=3. For this purpose, we implemented a MATLAB program which simulates the values of the signal, zk, the real measurements, ỹk, and the available ones, yk, considering different delay probabilities, and provides the filtering and fixed-point smoothing estimates of zk, as well as the corresponding error variances.As in the theoretical study, it is assumed that, at the initial time, the available measurement is equal to the real one,y1=ỹ1; at time k=2, the real measurement may be one-step randomly delayed; at time k=3 it may be randomly delayed either by one or two steps and, at any sampling time k≥4, the available measurement, yk, may be randomly delayed by one, two or three steps. The random delays in the observations have been simulated considering three independent sequences of independent Bernoulli random variables, {ζk1;k≥2}, {ζk2;k≥3}, and {ζk3;k≥4} with constant probabilities P(ζki=1)=qi,i=1,2,3, and defining the available measurements of the signal as (5.1)y1=ỹ1,y2=(1-ζ21)ỹ2+ζ21ỹ1,y3=(1-ζ31)ỹ3+ζ31(1-ζ32)ỹ2+ζ31ζ32ỹ1,yk=(1-ζk1)ỹk+ζk1(1-ζk2)ỹk-1+ζk1ζk2(1-ζk3)ỹk-2+ζk1ζk2ζk3ỹk-3,k≥4.
Hence, for d=0,1,2,3, the sequences {γk(d);k>d} in the theoretical study are given by (5.2)γk(0)=1-ζk1withpk(0)=1-q1,k>1,γ1(0)=1withp1(0)=1,γk(1)=ζk1(1-ζk2)withpk(1)=q1(1-q2),k>2,γ2(1)=ζ21withp2(1)=q1,γk(2)=ζk1ζk2(1-ζk3)withpk(2)=q1q2(1-q3),k>3,γ3(2)=ζ31ζ32withp3(2)=q1q2,γk(3)=ζk1ζk2ζk3withpk(3)=q1q2q3,k>3.Ifζk1=0, then yk=ỹk and the kth measurement is updated; if ζk1=1, the kth measurement is one-step delayed when ζk2=0 and two-step delayed when ζk2=1 and ζk3=0; finally, if ζki=1, for i=1,2,3, the kth measurement is three-step delayed.Note thatq1=1-pk(0) represents the probability of receiving a delayed observation at each sampling time. Moreover, note that pk(0), pk(1), and pk(2) are decreasing functions of q1, q2, and q3, respectively, while pk(3) is an increasing function of qi, i=1,2,3. Actually, as q1 increases, pk(0) decreases, but the delay probabilities pk(1), pk(2) and pk(3) increase; as q2 increases, pk(1) decreases, but the two- and three-step delay probabilities pk(2) and pk(3) increase; and, as q3 increases, pk(2) decreases, but the three-step delay probability pk(3) increases.Fork=1,…,25, Figure 1 displays the processed measurements, yk, taking qi=0.5, i=1,2,3; that is, when the probability that the real measurement is used in the estimation is pk(0)=1-q1=0.5 for k≥2, the one-step delay probability is pk(1)=0.5pk(0)=0.25 for k≥3, and the two- and three-step delay probabilities are both equal pk(2)=pk(3)=0.25pk(0)=0.125 for k≥4. This figure shows that, in fact, the measured output ỹk at time k can be received on time, delayed or lost in transmission. Actually, the available measurement being processed is(i)
updated (yk=ỹk) at k=1,2,4,5,7,9,11,13,15,17,18,21,22,23,25,(ii)
delayed by one sampling period (yk=ỹk-1) at k=3,8,10,20,24,(iii)
delayed by two sampling periods (yk=ỹk-2) at k=6,12,16,19,(iv)
delayed by three sampling periods (yk=ỹk-3) at k=14,whileỹ3,ỹ6,ỹ8,ỹ12,ỹ16,ỹ20, and ỹ24 are dropped out, thus not being processed at any sampling time. Also, it is worth noting that some measurements are rereceived, since the output at each instant is transmitted for D+1=4 consecutive times (for example, ỹ2 is rereceived at k=3, ỹ4 is rereceived at k=6, and ỹ11 is rereceived at k=14, among others).Next, the performance of the estimators is analyzed by investigating the error variances forqi=0.5,i=1,2,3. Figure 2 shows the error variances for the filter ẑk/k and the smoothers ẑk/k+1 and ẑk/k+2. As occurs for nondelayed observations, the error variances of the smoothing estimators are less than the filtering one; hence, the estimation accuracy of the smoothers is superior to that of the filter and also improves as the number of iterations in the fixed-point smoothing algorithm increases. The results of a simulated signal together with the filtering and smoothing estimates, ẑk/k and ẑk/k+2, are illustrated in Figure 3. According to the previous results, this figure shows that the evolution of this signal is followed more accurately by the smoothing estimates.Figure 2
Filtering and smoothing error variances forq1=q2=q3=0.5.Figure 3
Signalzk as well as filtering and smoothing estimates, ẑk/k and ẑk/k+2, when q1=q2=q3=0.5.To analyze the performance of the proposed estimators versus the delay probabilities, the error variances have been calculated for different values ofqi, i=1,2,3, which provide different values of the probabilities pk(i), i=0,1,2,3. In all the cases examined, the estimation-error variances present insignificant variation from the 15th iteration onwards and, consequently, only the error variances at a specific iteration are shown here. Figure 4 displays the filtering and smoothing error variances at k=100 versus q1 (for q2=0.3,0.5,and0.7 and q3=0.5). From this figure it is deduced that, as q1 increases (and, consequently, the nondelay probability pk(0) decreases), the estimation-error variances become greater and, hence, worse estimations are obtained.Figure 4
Error variances versusq1, when q2=0.3,0.5,and0.7 and q3=0.5.On the other hand, for each value ofq1 (keeping q3=0.5 fixed), the error variances become smaller as q2 decreases, which means that the estimations are better. This was expectable since, as q2 decreases, the one-step delay probability increases and the two- and three-step delay probabilities decrease.Also, as expected, this improvement is more significant asq1 increases (or, equivalently, as the delay probability increases), being more appreciable for the filter than for the smoother.Next, we compare the filtering error variances atk=100 versus q2, for q3=0.3,0.5,and0.7 and q1=0.9 (similar results are obtained for other fixed values of the probability of delay, q1, but as indicated above, the comparison is clearer when such probability is large). The results are displayed in Figure 5, which shows that, as q2 increases (which means that the one-step delay probability decreases, but the two- and three-step delay probabilities increase), the estimation-error variances become greater and, consequently, the accuracy of the estimators is worse. Moreover, for each value of q2, the estimators perform worse with the increasing of q3, which is reasonable since the three-step delay probability increases with q3. Nevertheless, for small values of q2, the difference is almost imperceptible since, as q2 becomes close to zero, the one-step delay probability tends to q1 (fixed), and the two- and three-step delay probabilities both tend to zero. This consideration is confirmed in Figure 6, where the filtering and smoothing error variances at k=100 versus q3 (for q2=0.3,0.5,and0.7 and q1=0.9) are displayed.Figure 5
Error variances versusq2, when q3=0.3,0.5,and0.7 and q1=0.9.Figure 6
Error variances versusq3, when q2=0.3,0.5,and0.7 and q1=0.9.


## 6. Conclusions

In this paper, a recursive least-squares linear estimation algorithm is proposed to estimate signals from observations which can be randomly delayed or lost in transmission, a realistic and common assumption in networked control systems where, generally, transmission delays and packet losses are unavoidable due to the unreliable network characteristics. The largest delay and the maximum number of consecutive dropouts are assumed to be upper bounded by a known constantD, and the random measurement delays and packet dropouts are described by introducing D+1 sequences of Bernoulli random variables, whose parameters represent the delay probabilities. Thus, the current study generalizes the results in [19] to the case of multiple-sample delays in the observations.Using an innovation approach, the estimation algorithm is derived without requiring the knowledge of the signal state-space model, but only the covariance functions of the processes involved in the observation equation, as well as the delay probabilities. To measure the performance of the estimators, the estimation-error covariance matrices are also calculated.To illustrate the theoretical results established in this paper, a simulation example is presented, in which the proposed algorithm is applied to estimate a signal from bounded random measurement delays and packet dropouts, assuming that the largest delay and the maximum number of successive dropouts isD=3.



---

*Source: 258065-2010-05-03.xml*