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

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.

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 optimal 𝐻₂ 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 𝐻_∞ 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 Section 2 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 by 𝐷); 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 𝐷 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, 𝑧_𝑘, whose measured output at the sampling time 𝑘, denoted by 𝑦_𝑘, is perturbed by an additive noise vector 𝑣_𝑘; that is,𝑦_𝑘=𝑧_𝑘+𝑣_𝑘,𝑘≥1. (2.1)

The measured output 𝑦_𝑘 is transmitted during the sampling times 𝑘,𝑘+1,…, and 𝑘+𝐷 but, at each sampling time 𝑘, only one of the measurements 𝑦_𝑘−𝐷,…,𝑦_𝑘 is processed; consequently, at time 𝑘>𝐷, the processed measurement can be either delayed by 𝑑=1,…,𝐷 sampling periods with a known probability 𝑝_𝑘^(𝑑), or updated with probability 𝑝_𝑘⁽⁰⁾. At the initial time 𝑘=1, the measured output 𝑦₁ is always available (𝑝₁⁽⁰⁾=1) and, hence, the processed measurement is equal to the real measurement 𝑦₁=𝑦₁. At any time 𝑘≤𝐷, the processed measurement only can be delayed by 𝑑=1,…,𝑘−1 sampling periods, since only 𝑦₁,…,𝑦_𝑘 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:𝑦_𝑘=^{min(𝑘−1,𝐷)}_𝑑=0𝛾_𝑘^(𝑑)𝑦_𝑘−𝑑,𝑘≥1, (2.2) where, for 𝑑=0,1,…,𝐷,{𝛾_𝑘^(𝑑);𝑘>𝑑} denote sequences of mutually independent Bernoulli random variables with 𝑃[𝛾_𝑘^(𝑑)=1]=𝑝_𝑘^(𝑑) and ∑^{min(𝑘−1,𝐷)}_𝑑=0𝛾_𝑘^(𝑑)=1.

This fact guarantees that if 𝛾_𝑘⁽⁰⁾=1, 𝛾_𝑘^(𝑑)=0 for all 𝑑≠0, which means that, with probability 𝑝_𝑘⁽⁰⁾, the 𝑘th measurement is received and processed on time. If 𝛾_𝑘⁽⁰⁾=0, then there exists one and only one 𝑑=1,…,𝐷 such that 𝛾_𝑘^(𝑑)=1, which means that the measurement 𝑦_𝑘 is delayed by 𝑑 sampling periods with probability 𝑝_𝑘^(𝑑). Note that the measured output 𝑦_𝑘 at any time 𝑘 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 𝐷+1 consecutive times (please, see Figure 1 in Section 5 for specific results when 𝐷=3).

Figure 1 

Processed observations for $D=3$ when $q^1=q^2=q^3=0.5$.

The signal estimation problem is addressed based on the following assumptions.

Assumption 1. The 𝑛-dimensional signal process {𝑧_𝑘;𝑘≥1} has zero mean and autocovariance function 𝐾_𝑘,𝑠=𝐸[𝑧_𝑘𝑧^𝑇_𝑠]=𝐴_𝑘𝐵^𝑇_𝑠,𝑠≤𝑘, where 𝐴 and 𝐵 are known matrix functions.

Assumption 2. The noise process {𝑣_𝑘;𝑘≥1} is a zero-mean white sequence with known autocovariance function 𝐸[𝑣_𝑘𝑣^𝑇_𝑘]=𝑅_𝑘,for all 𝑘≥1.

Assumption 3. For 𝑑=0,1,…,𝐷,{𝛾_𝑘^(𝑑);𝑘>𝑑} are sequences of independent Bernoulli random variables with 𝑃[𝛾_𝑘^(𝑑)=1]=𝑝_𝑘^(𝑑) and ∑^{min(𝑘−1,𝐷)}_𝑑=0𝛾_𝑘^(𝑑)=1.

Moreover, {𝛾_𝑘^(𝑑),0≤𝑑≤min(𝑘−1,𝐷)} is independent of {𝛾_𝑗^(𝑑),0≤𝑑≤min(𝑗−1,𝐷)} for all 𝑘≠𝑗.

Assumption 4. For each 𝑑=0,1,…,𝐷, the processes {𝑧_𝑘;𝑘≥1},{𝑣_𝑘;𝑘≥1}, and {𝛾_𝑘^(𝑑);𝑘>𝑑} 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 Assumption 1 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 𝐸[𝑧_𝑘𝑧^𝑇_𝑠]=Φ^𝑘−𝑠𝐸[𝑧_𝑠𝑧^𝑇_𝑠],𝑠≤𝑘, and Assumption 1 is clearly satisfied, taking 𝐴_𝑘=Φ^𝑘 and 𝐵_𝑠=𝐸[𝑧_𝑠𝑧^𝑇_𝑠](Φ^−𝑠)^𝑇. Also, processes with finite-dimensional, possibly time-variant, state-space models have semiseparable covariance functions, 𝐸[𝑧_𝑘𝑧^𝑇_𝑠]=∑^𝑟_𝑖=1𝑎^𝑖_𝑘𝑏^𝑖_𝑠,𝑠≤𝑘 (see [20]), and this structure is a particular case of that assumed, just taking 𝐴_𝑘=(𝑎¹_𝑘,𝑎²_𝑘,…,𝑎^𝑟_𝑘) and 𝐵_𝑠=(𝑏¹_𝑠,𝑏²_𝑠,…,𝑏^𝑟_𝑠). 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, 𝑧_𝑘/𝐿, of the signal 𝑧_𝑘 based on the observations {𝑦₁,…,𝑦_𝐿}. Specifically, our aim is to obtain the filter (𝐿=𝑘) and the fixed-point smoother (𝑘 fixed and 𝐿>𝑘) from recursive formulas.

The estimator 𝑧_𝑘/𝐿 is the orthogonal projection of the vector 𝑧_𝑘 onto ℒ(𝑦₁,…,𝑦_𝐿), the 𝑛-dimensional linear space spanned by {𝑦₁,…,𝑦_𝐿}. Since the observations are generally nonorthogonal vectors, we use an innovation approach, based on an orthogonalization procedure by means of which the observation process {𝑦_𝑘;𝑘≥1} is transformed into an equivalent one (innovation process) of orthogonal vectors {𝜈_𝑘;𝑘≥1}, equivalent in the sense that each set {𝜈₁,…,𝜈_𝐿} spans the same linear subspace as {𝑦₁,…,𝑦_𝐿}.

Since the innovations constitute a white process, this methodology allows us to find the orthogonal projection of the vector 𝑧_𝑘 onto {𝜈₁,…,𝜈_𝐿} 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:𝑧_𝑘/𝐿=^𝐿_𝑖=1𝑆^𝑧_𝑘,𝑖Π_𝑖⁻¹𝜈_𝑖, (3.1) where 𝑆^𝑧_𝑘,𝑖=𝐸[𝑧_𝑘𝜈^𝑇_𝑖] and Π_𝑖=𝐸[𝜈_𝑖𝜈^𝑇_𝑖].

Hence, to obtain the signal estimators it is necessary to find previously an explicit formula for the innovations and their covariance matrices.

Innovation 𝜈_𝑘 and Covariance Matrix Π_𝑘
The innovation at time 𝑘 is defined as 𝜈_𝑘=𝑦_𝑘−𝑦_{𝑘/𝑘−1} where 𝑦_{𝑘/𝑘−1} is the one-stage linear predictor of 𝑦_𝑘. Since the LS linear estimator 𝑦_{𝑘/𝑘−1} is the projection of 𝑦_𝑘 onto ℒ(𝜈₁,…,𝜈_𝑘−1), from (2.1) and (2.2) and model assumptions, it is clear that 𝑦_{𝑘/𝑘−1}=^{min(𝑘−1,𝐷)}_𝑑=0𝑝_𝑘^(𝑑)(𝑧_{𝑘−𝑑/𝑘−1}+𝑣_{𝑘−𝑑/𝑘−1}),𝑘≥2,𝑦_1/0=0, (3.2) where 𝑣_𝑖/𝑗 denotes the orthogonal projection of the vector 𝑣_𝑖 onto ℒ(𝜈₁,…,𝜈_𝑗).

Similar that to (3.1), by denoting 𝑆^𝑣_𝑘,𝑖=𝐸[𝑣_𝑘𝜈^𝑇_𝑖], the noise estimators are expressed as follows:𝑣_𝑘/𝐿=^𝐿_𝑖=1𝑆^𝑣_𝑘,𝑖Π_𝑖⁻¹𝜈_𝑖. (3.3) From the model assumptions, the noise 𝑣_𝑘 is independent of {𝑦₁,…,𝑦_𝑘−1} and hence, for 𝑖<𝑘, we have that 𝑆^𝑣_𝑘,𝑖=0; thus, the one-stage predictor of the noise is 𝑣_{𝑘/𝑘−1}=0. Consequently, the innovation is given by 𝜈_𝑘=𝑦_𝑘−𝑝_𝑘⁽⁰⁾𝑧_{𝑘/𝑘−1}−^{min(𝑘−1,𝐷)}_𝑑=1𝑝_𝑘^(𝑑)(𝑧_{𝑘−𝑑/𝑘−1}+𝑣_{𝑘−𝑑/𝑘−1}),𝑘≥2,𝜈₁=𝑦₁, (3.4) and it is necessary to obtain the signal predictor 𝑧_{𝑘/𝑘−1} and the estimators 𝑧_{𝑘−𝑑/𝑘−1} and 𝑣_{𝑘−𝑑/𝑘−1} for 𝑑=1,…,min(𝑘−1,𝐷), that is, the filters and the smoothers of the signal and noise, respectively.

Remark 3.1. From the model assumptions, it is clear that 𝑆^𝑣_𝑘,𝑘=𝑝_𝑘⁽⁰⁾𝑅_𝑘 and, consequently, the filter of the noise is 𝑣_𝑘/𝑘=𝑝_𝑘⁽⁰⁾𝑅_𝑘Π_𝑘⁻¹𝜈_𝑘.

To obtain the covariance matrix Π_𝑘, from (2.1), (2.2), and (3.2) and taking into account that 𝑣_{𝑘/𝑘−1}=0, the innovation 𝜈_𝑘 is expressed as follows:𝜈_𝑘=^{min(𝑘−1,𝐷)}_𝑑=0(𝛾_𝑘^(𝑑)−𝑝_𝑘^(𝑑))𝑧_𝑘−𝑑+^{min(𝑘−1,𝐷)}_𝑑=1(𝛾_𝑘^(𝑑)−𝑝_𝑘^(𝑑))𝑣_𝑘−𝑑+𝛾_𝑘⁽⁰⁾𝑣_𝑘+∑^{min(𝑘−1,𝐷)}_𝑑=0𝑝_𝑘^(𝑑)(𝑧_𝑘−𝑑−𝑧_{𝑘−𝑑/𝑘−1})+∑^{min(𝑘−1,𝐷)}_𝑑=1𝑝_𝑘^(𝑑)(𝑣_𝑘−𝑑−𝑣_{𝑘−𝑑/𝑘−1}). (3.5)

Using this expression and the model assumptions, we haveΠ_𝑘=^{min(𝑘−1,𝐷)}_{𝑑,𝑑′=0}Cov(𝛾_𝑘^(𝑑),𝛾^(𝑑′)_𝑘)𝐾_{𝑘−𝑑,𝑘−𝑑′}+^{min(𝑘−1,𝐷)}_{𝑑,𝑑′=0}𝑝_𝑘^(𝑑)𝑝^(𝑑′)_𝑘𝑃^𝑧𝑧_{𝑘−𝑑,𝑘−𝑑′/𝑘−1}+^{min(𝑘−1,𝐷)}_𝑑=1𝑝_𝑘^(𝑑)(1−𝑝_𝑘^(𝑑))𝑅_𝑘−𝑑+𝑝_𝑘⁽⁰⁾𝑅_𝑘+^{min(𝑘−1,𝐷)}_{𝑑,𝑑′=1}𝑝_𝑘^(𝑑)𝑝^(𝑑′)_𝑘𝑃^𝑣𝑣_{𝑘−𝑑,𝑘−𝑑′/𝑘−1}+^{min(𝑘−1,𝐷)}_{𝑑=0,𝑑′=1}𝑝_𝑘^(𝑑)𝑝^(𝑑′)_𝑘𝑃^𝑧𝑣_{𝑘−𝑑,𝑘−𝑑′/𝑘−1}+^{min(𝑘−1,𝐷)}_{𝑑=1,𝑑′=0}𝑝_𝑘^(𝑑)𝑝^(𝑑′)_𝑘𝑃^𝑣𝑧_{𝑘−𝑑,𝑘−𝑑′/𝑘−1}Π,𝑘≥2,₁=𝐴₁𝐵^𝑇₁+𝑅₁, (3.6) where 𝑃^𝑧𝑧_{𝑙,𝑚/𝑁}=𝐸[(𝑧_𝑙−𝑧_𝑙/𝑁)(𝑧_𝑚−𝑧_𝑚/𝑁)] and 𝑃^𝑣𝑣_{𝑙,𝑚/𝑁}=𝐸[(𝑣_𝑙−𝑣_𝑙/𝑁)(𝑣_𝑚−𝑣_𝑚/𝑁)] denote the estimation-error covariance matrices of the state and noise, respectively, and 𝑃^𝑧𝑣_{𝑙,𝑚/𝑁}=𝐸[(𝑧_𝑙−𝑧_𝑙/𝑁)(𝑣_𝑚−𝑣_𝑚/𝑁)] 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 (Section 4.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, 𝑧_{𝑘/𝑘−1}, it is necessary to calculate the coefficients 𝑆^𝑧_𝑘,𝑖=𝐸[𝑧_𝑘𝜈^𝑇_𝑖]=𝐸[𝑧_𝑘𝑦^𝑇_𝑖]−𝐸[𝑧_𝑘𝑦^𝑇_{𝑖/𝑖−1}],𝑖<𝑘, (4.1) so 𝐸[𝑧_𝑘𝑦^𝑇_𝑖] and 𝐸[𝑧_𝑘𝑦^𝑇_{𝑖/𝑖−1}] must be calculated for 𝑖<𝑘.

(i) On the one hand, using (2.1)-(2.2) and the model assumptions, it is clear that 𝐸[𝑧_𝑘𝑦^𝑇_𝑖]=^{min(𝑖−1,𝐷)}_𝑑=0𝑝_𝑖^(𝑑)𝐸[𝑧_𝑘𝑧^𝑇_𝑖−𝑑]=^{min(𝑖−1,𝐷)}_𝑑=0𝑝_𝑖^(𝑑)𝐴_𝑘𝐵^𝑇_𝑖−𝑑=𝐴_𝑘𝐺^𝑇_𝐵𝑖,1≤𝑖≤𝑘, (4.2) where 𝐺_𝐵𝑖=∑^{min(𝑖−1,𝐷)}_𝑑=0𝑝_𝑖^(𝑑)𝐵_𝑖−𝑑,𝑖≥1.

(ii) On the other hand, from (3.2) and taking into account that 𝑣_{𝑖/𝑖−1}=0, we have that 𝐸[𝑧_𝑘𝑦^𝑇_{𝑖/𝑖−1}]=^{min(𝑖−1,𝐷)}_𝑑=0𝑝_𝑖^(𝑑)𝐸[𝑧_𝑘𝑧^𝑇_{𝑖−𝑑/𝑖−1}]+^{min(𝑖−1,𝐷)}_𝑑=1𝑝_𝑖^(𝑑)𝐸[𝑧_𝑘𝑣^𝑇_{𝑖−𝑑/𝑖−1}],2≤𝑖≤𝑘, (4.3) so 𝐸[𝑧_𝑘𝑧^𝑇_{𝑖−𝑑/𝑖−1}] and 𝐸[𝑧_𝑘𝑣^𝑇_{𝑖−𝑑/𝑖−1}] must be calculated for 𝑖=1,…,𝑘.(a)Using (3.1) for 𝑧^𝑇_{𝑖−𝑑/𝑖−1}, 𝐸[𝑧_𝑘𝑧^𝑇_{𝑖−𝑑/𝑖−1}]=𝐸[𝑧_𝑘{^𝑖−1_𝑗=1𝑆^𝑧_{𝑖−𝑑,𝑗}Π_𝑗⁻¹𝜈_𝑗}^𝑇]=^𝑖−1_𝑗=1𝑆^𝑧_𝑘,𝑗Π_𝑗⁻¹𝑆^𝑧𝑇_{𝑖−𝑑,𝑗},2≤𝑖≤𝑘. (4.4)(b)Using (3.3) for 𝑣^𝑇_{𝑖−𝑑/𝑖−1}, and since 𝑆^𝑣_{𝑖−𝑑,𝑗}=0, for 𝑗<𝑖−𝑑, we obtain 𝐸[𝑧_𝑘𝑣^𝑇_{𝑖−𝑑/𝑖−1}]=𝐸[𝑧_𝑘{^𝑖−1_{𝑗=𝑖−𝑑}𝑆^𝑣_{𝑖−𝑑,𝑗}Π_𝑗⁻¹𝜈_𝑗}^𝑇]=^𝑖−1_{𝑗=𝑖−𝑑}𝑆^𝑧_𝑘,𝑗Π_𝑗⁻¹𝑆^𝑣𝑇_{𝑖−𝑑,𝑗},2≤𝑖≤𝑘. (4.5)

Substituting into (4.1) the expectations calculated in (i) and (ii), we obtain 𝑆^𝑧_𝑘,𝑖=𝐴_𝑘𝐺^𝑇_𝐵𝑖−^{min(𝑖−1,𝐷)}_𝑑=0𝑝_𝑖^(𝑑){^𝑖−1_𝑗=1𝑆^𝑧_𝑘,𝑗Π_𝑗⁻¹𝑆^𝑧𝑇_{𝑖−𝑑,𝑗}}−^{min(𝑖−1,𝐷)}_𝑑=1𝑝_𝑖^(𝑑){^𝑖−1_{𝑗=𝑖−𝑑}𝑆^𝑧_𝑘,𝑗Π_𝑗⁻¹𝑆^𝑣𝑇_{𝑖−𝑑,𝑗}𝑆},2≤𝑖≤𝑘,_𝑧𝑘,1=𝐴_𝑘𝐺^𝑇_𝐵1. (4.6) If we now introduce a function 𝐽_𝑖 satisfying that𝐽_𝑖=𝐺^𝑇_𝐵𝑖−^{min(𝑖−1,𝐷)}_𝑑=0𝑝_𝑖^(𝑑){^𝑖−1_𝑗=1𝐽_𝑗Π_𝑗⁻¹𝑆^𝑧𝑇_{𝑖−𝑑,𝑗}}−^{min(𝑖−1,𝐷)}_𝑑=1𝑝_𝑖^(𝑑){^𝑖−1_{𝑗=𝑖−𝑑}𝐽_𝑗Π_𝑗⁻¹𝑆^𝑣𝑇_{𝑖−𝑑,𝑗}𝐽},𝑖≥2,₁=𝐺^𝑇_𝐵1, (4.7) we conclude that𝑆^𝑧_𝑘,𝑖=𝐴_𝑘𝐽_𝑖,1≤𝑖≤𝑘. (4.8)

Substituting now (4.8) into (3.1) for 𝐿=𝑘−1, the following expression for the state predictor is obtained: 𝑧_{𝑘/𝑘−1}=𝐴_𝑘𝑂_𝑘−1,𝑘≥1, (4.9) where the vectors 𝑂_𝑘 are defined by 𝑂_𝑘=^𝑘_𝑖=1𝐽_𝑖Π_𝑖⁻¹𝜈_𝑖,𝑘≥1,𝑂₀=0, (4.10) thus satisfying the recursive relation 𝑂_𝑘=𝑂_𝑘−1+𝐽_𝑘Π_𝑘⁻¹𝜈_𝑘,𝑘≥1,𝑂₀=0. (4.11)

Hence, an expression for 𝐽_𝑘 must be derived.

Remark 4.1. By substituting (4.8) into (3.1), it is also clear that the 𝑙-stage predictors are given by 𝑧_{𝑘+𝑙/𝑘}=𝐴_𝑘+𝑙𝑂_𝑘,𝑘≥1,𝑙≥1.

Expression for 𝐽_𝑘
First, expression (4.7) is just rewritten for 𝑖=𝑘: 𝐽_𝑘=𝐺^𝑇_𝐵𝑘−^{min(𝑘−1,𝐷)}_𝑑=0𝑝_𝑘^(𝑑){^𝑘−1_𝑗=1𝐽_𝑗Π_𝑗⁻¹𝑆^𝑧𝑇_{𝑘−𝑑,𝑗}}−^{min(𝑘−1,𝐷)}_𝑑=1𝑝_𝑘^(𝑑){^𝑘−1_{𝑗=𝑘−𝑑}𝐽_𝑗Π_𝑗⁻¹𝑆^𝑣𝑇_{𝑘−𝑑,𝑗}},𝑘≥2. (4.12) For 𝑘=2, using that 𝑆^𝑧_2,1=𝐴₂𝐽₁,𝑆^𝑧_1,1=𝐴₁𝐽₁, and 𝑆^𝑣_1,1=𝑝₁⁽⁰⁾𝑅₁, and denoting 𝑟₁=𝐽₁Π₁⁻¹𝐽^𝑇₁, it is immediately clear that 𝐽₂=𝐺^𝑇_𝐵2−𝑝₂⁽⁰⁾𝑟₁𝐴^𝑇₂−𝑝₂⁽¹⁾{𝑟₁𝐴^𝑇₁+𝑝₁⁽⁰⁾𝐽₁Π₁⁻¹𝑅₁}. (4.13) For 𝑘≥3, we examine separately the sums that appear in 𝐽_𝑘 as follows.
(i) Since 𝑆^𝑧_𝑚,𝑗=𝐴_𝑚𝐽_𝑗, for 𝑗≤𝑚, by denoting that 𝑟_𝑘=∑^𝑘_𝑗=1𝐽_𝑗Π_𝑗⁻¹𝐽^𝑇_𝑗, we have that^min(𝑘∑^−1,𝐷)_𝑑=0𝑝_𝑘^(𝑑){∑^𝑘−1_𝑗=1𝐽_𝑗Π_𝑗⁻¹𝑆^𝑧𝑇_{𝑘−𝑑,𝑗}}=𝑝_𝑘⁽⁰⁾𝑟_𝑘−1𝐴^𝑇_𝑘+𝑝_𝑘⁽¹⁾𝑟_𝑘−1𝐴^𝑇_𝑘−1+∑^{min(𝑘−1,𝐷)}_𝑑=2𝑝_𝑘^(𝑑)𝑟_𝑘−𝑑𝐴^𝑇_𝑘−𝑑+∑^{min(𝑘−1,𝐷)}_𝑑=2𝑝∑_𝑘^{(𝑑)𝑑−1}_𝑖=1𝐽_{𝑘−𝑑+𝑖}Π⁻¹_{𝑘−𝑑+𝑖}𝑆^𝑧𝑇_{𝑘−𝑑,𝑘−𝑑+𝑖},𝑘≥3.
(ii) Since 𝑆^𝑣_𝑚,𝑚=𝑝_𝑚⁽⁰⁾𝑅_𝑚,𝑚≥1, we have that ^{min(𝑘−1,𝐷)}_𝑑=1𝑝_𝑘^(𝑑){^𝑘−1_{𝑗=𝑘−𝑑}𝐽_𝑗Π_𝑗⁻¹𝑆^𝑣𝑇_{𝑘−𝑑,𝑗}}=𝑝_𝑘⁽¹⁾𝐽_𝑘−1Π⁻¹_𝑘−1𝑝⁽⁰⁾_𝑘−1𝑅_𝑘−1+^{min(𝑘−1,𝐷)}_𝑑=2𝑝_𝑘^(𝑑){𝐽_𝑘−𝑑Π⁻¹_𝑘−𝑑𝑝⁽⁰⁾_𝑘−𝑑𝑅_𝑘−𝑑+^𝑑−1_𝑖=1𝐽_{𝑘−𝑑+𝑖}Π⁻¹_{𝑘−𝑑+𝑖}𝑆^𝑣𝑇_{𝑘−𝑑,𝑘−𝑑+𝑖}},𝑘≥3. (4.15)By substituting the above two sums in 𝐽_𝑘, it is deduced that 𝐽_𝑘=𝐺^𝑇_𝐵𝑘−𝑝_𝑘⁽⁰⁾𝑟_𝑘−1𝐴^𝑇_𝑘−^{min(𝑘−1,𝐷)}_𝑑=1𝑝_𝑘^(𝑑){𝑟_𝑘−𝑑𝐴^𝑇_𝑘−𝑑+𝑝⁽⁰⁾_𝑘−𝑑𝐽_𝑘−𝑑Π⁻¹_𝑘−𝑑𝑅_𝑘−𝑑}−^{min(𝑘−1,𝐷)}_𝑑=2𝑝_𝑘^{(𝑑)𝑑−1}_𝑖=1𝐽_{𝑘−𝑑+𝑖}Π⁻¹_{𝑘−𝑑+𝑖}{𝑆^𝑧𝑇_{𝑘−𝑑,𝑘−𝑑+𝑖}+𝑆^𝑣𝑇_{𝑘−𝑑,𝑘−𝑑+𝑖}},𝑘≥3. (4.16) Finally, note that 𝑆^𝑧_{𝑘−𝑑,𝑘−𝑑+𝑖} and 𝑆^𝑣_{𝑘−𝑑,𝑘−𝑑+𝑖} are the smoothing gain matrices (which will be obtained in the next section), and the matrix 𝑟_𝑘 is recursively obtained by 𝑟_𝑘=𝑟_𝑘−1+𝐽_𝑘Π_𝑘⁻¹𝐽^𝑇_𝑘,𝑘≥1,𝑟₀=0. (4.17)

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:𝑧_𝑘/𝐿=𝑧_{𝑘/𝐿−1}+𝑆^𝑧_𝑘,𝐿Π_𝐿⁻¹𝜈_𝐿𝑣,𝐿≥𝑘,𝑘≥1,_𝑘/𝐿=𝑣_{𝑘/𝐿−1}+𝑆^𝑣_𝑘,𝐿Π_𝐿⁻¹𝜈_𝐿,𝐿≥𝑘,𝑘≥1, (4.18) with initial conditions given by the one-stage predictors 𝑧_{𝑘/𝑘−1}=𝐴_𝑘𝑂_𝑘−1 and 𝑣_{𝑘/𝑘−1}=0, respectively. Hence, the gain matrices 𝑆^𝑧_𝑘,𝐿=𝐸[𝑧_𝑘𝜈^𝑇_𝐿] and 𝑆^𝑣_𝑘,𝐿=𝐸[𝑣_𝑘𝜈^𝑇_𝐿] must be calculated for 𝐿≥𝑘.

Gain Matrices
Using (3.5) for 𝑘=𝐿 and the model assumptions, it is derived that 𝑆^𝑧_𝑘,𝐿=^{min(𝐿−1,𝐷)}_𝑑=0𝑝_𝐿^(𝑑)𝐸[𝑧_𝑘(𝑧_𝐿−𝑑−𝑧_{𝐿−𝑑/𝐿−1})^𝑇]+^{min(𝐿−1,𝐷)}_𝑑=1𝑝_𝐿^(𝑑)𝐸[𝑧_𝑘(𝑣_𝐿−𝑑−𝑣_{𝐿−𝑑/𝐿−1})^𝑇], (4.19) and, since the estimation-errors are orthogonal to the estimators, it is concluded that 𝑆^𝑧_𝑘,𝐿=^{min(𝐿−1,𝐷)}_𝑑=0𝑝_𝐿^(𝑑)𝑃^𝑧𝑧_{𝑘,𝐿−𝑑/𝐿−1}+^{min(𝐿−1,𝐷)}_𝑑=1𝑝_𝐿^(𝑑)𝑃^𝑧𝑣_{𝑘,𝐿−𝑑/𝐿−1}, (4.20) where 𝑃^𝑧𝑧_{𝑘,𝐿−𝑑/𝐿−1} and 𝑃^𝑧𝑣_{𝑘,𝐿−𝑑/𝐿−1} are the covariance and cross-covariance matrices of the errors 𝑧_{𝑘/𝐿−1}=𝑧_𝑘−𝑧_{𝑘/𝐿−1} with 𝑧_{𝐿−𝑑/𝐿−1} and 𝑣_{𝐿−𝑑/𝐿−1}=𝑣_𝐿−𝑑−𝑣_{𝐿−𝑑/𝐿−1}, respectively.

A similar reasoning leads to the following expression: 𝑆^𝑣_𝑘,𝐿=^{min(𝐿−1,𝐷)}_𝑑=0𝑝_𝐿^(𝑑)𝑃^𝑣𝑧_{𝑘,𝐿−𝑑/𝐿−1}+^{min(𝐿−1,𝐷)}_𝑑=1𝑝_𝐿^(𝑑)𝑃^𝑣𝑣_{𝑘,𝐿−𝑑/𝐿−1}, (4.21) where 𝑃^𝑣𝑣_{𝑘,𝐿−𝑑/𝐿−1} and 𝑃^𝑣𝑧_{𝑘,𝐿−𝑑/𝐿−1} are the covariance and cross-covariance matrices of the errors 𝑣_{𝑘/𝐿−1} with 𝑣_{𝐿−𝑑/𝐿−1} and 𝑧_{𝐿−𝑑/𝐿−1}, respectively.

Remark 4.2. Note that, from (4.8), the gain matrix of the signal filter is 𝑆^𝑧_𝑘,𝑘=𝐴_𝑘𝐽_𝑘 and, from (3.1), the filter is 𝑧_𝑘/𝑘=𝐴_𝑘𝑂_𝑘.

4.3. Error Covariance and Cross-Covariance Matrices

From the recursive relations (4.18), the estimation-errors admit the following expressions: 𝑧_𝑘/𝐿=𝑧_{𝑘/𝐿−1}−𝑆^𝑧_𝑘,𝐿Π_𝐿⁻¹𝜈_𝐿𝑣,𝐿≥𝑘,𝑘≥1,_𝑘/𝐿=𝑣_{𝑘/𝐿−1}−𝑆^𝑣_𝑘,𝐿Π_𝐿⁻¹𝜈_𝐿,𝐿≥𝑘,𝑘≥1. (4.22) Now, since the estimation-errors are orthogonal to the estimators, we have that 𝐸[𝑧_{𝑘/𝐿−1}𝜈^𝑇_𝐿]=𝐸[𝑧_𝑘𝜈^𝑇_𝐿]=𝑆^𝑧_𝑘,𝐿 and 𝐸[𝑣_{𝑘/𝐿−1}𝜈^𝑇_𝐿]=𝐸[𝑣_𝑘𝜈^𝑇_𝐿]=𝑆^𝑣_𝑘,𝐿, thus deducing the following recursive expressions for the error covariance and cross-covariance matrices:𝑃^𝑧𝑧_{𝑘,𝑚/𝐿}=𝑃^𝑧𝑧_{𝑘,𝑚/𝐿−1}−𝑆^𝑧_𝑘,𝐿Π_𝐿⁻¹𝑆^𝑧𝑇_𝑚,𝐿𝑃,𝐿≥𝑘,𝑚,𝑘,𝑚≥1,_{𝑣𝑣𝑘,𝑚/𝐿}=𝑃^𝑣𝑣_{𝑘,𝑚/𝐿−1}−𝑆^𝑣_𝑘,𝐿Π_𝐿⁻¹𝑆^𝑣𝑇_𝑚,𝐿𝑃,𝐿≥𝑘,𝑚,𝑘,𝑚≥1,_{𝑧𝑣𝑘,𝑚/𝐿}=𝑃^𝑧𝑣_{𝑘,𝑚/𝐿−1}−𝑆^𝑧_𝑘,𝐿Π_𝐿⁻¹𝑆^𝑣𝑇_𝑚,𝐿,𝐿≥𝑘,𝑚,𝑘,𝑚≥1. (4.23)

The initial conditions of these equations, the prediction error covariance and cross-covariance matrices, are obtained as follows: 𝑃^𝑧𝑧_{𝑘+𝑖,𝑘/𝑘−1}=𝐴_𝑘+𝑖[𝐵^𝑇_𝑘−𝑟_𝑘−1𝐴^𝑇_𝑘𝑃],𝑖≥0,𝑘≥1,_{𝑣𝑣𝑘+𝑖,𝑘/𝑘−1}=𝑅_𝑘𝛿_𝑖,0𝑃,𝑖≥0,𝑘≥1,_{𝑧𝑣𝑘+𝑖,𝑘/𝑘−1}=𝑃^𝑧𝑣_{𝑘,𝑘+𝑖/𝑘−1}=0,𝑖≥0,𝑘≥1, (4.24) where 𝛿 denotes the Kronecker delta function.

In fact, by writing 𝑃^𝑧𝑧_{𝑘+𝑖,𝑘/𝑘−1}=𝐸[𝑧_𝑘+𝑖𝑧^𝑇_𝑘]−𝐸[𝑧_{𝑘+𝑖/𝑘−1}𝑧^𝑇_{𝑘/𝑘−1}] and using that 𝑧_{𝑘+𝑖/𝑘−1}=𝐴_𝑘+𝑖𝑂_𝑘−1 and 𝐸[𝑂_𝑘−1𝑂^𝑇_𝑘−1]=𝑟_𝑘−1, the expression for 𝑃^𝑧𝑧_{𝑘+𝑖,𝑘/𝑘−1} is immediate from Assumption 1. Analogously, using that 𝑣_{𝑘+𝑖/𝑘−1}=0, for 𝑖≥0, and Assumptions 2 and 4, the expressions for 𝑃^𝑣𝑣_{𝑘+𝑖,𝑘/𝑘−1}, 𝑃^𝑧𝑣_{𝑘+𝑖,𝑘/𝑘−1}, and 𝑃^𝑧𝑣_{𝑘,𝑘+𝑖/𝑘−1} are derived.

4.4. Computational Procedure

At the sampling time 𝑘, once the (𝑘−1)th iteration is finished and the new observation 𝑦_𝑘 is available, the proposed estimation algorithm operates as follows.(i) Compute the innovation 𝜈_𝑘 and its covariance matrix Π_𝑘 by (3.4) and (3.6), respectively.(ii) Compute 𝐽_𝑘 by (4.16) and, from it, the filtering gain matrix 𝑆^𝑧_𝑘,𝑘=𝐴_𝑘𝐽_𝑘.(iii) Compute the filter 𝑧_𝑘/𝑘=𝑧_{𝑘/𝑘−1}+𝑆^𝑧_𝑘,𝑘Π_𝑘⁻¹𝜈_𝑘 and the filtering error covariance matrix 𝑃^𝑧𝑧_{𝑘,𝑘/𝑘}=𝑃^𝑧𝑧_{𝑘,𝑘/𝑘−1}−𝑆^𝑧_𝑘,𝑘Π_𝑘⁻¹𝑆^𝑧𝑇_𝑘,𝑘.(iv) To implement the above steps at time 𝑘+1, we need to(a) compute 𝑂_𝑘=𝑂_𝑘−1+𝐽_𝑘Π_𝑘⁻¹𝜈_𝑘 and, from it, the signal predictor 𝑧_𝑘+1/𝑘=𝐴_𝑘+1𝑂_𝑘,(b) compute the noise filter 𝑣_𝑘/𝑘=𝑝_𝑘⁽⁰⁾𝑅_𝑘Π_𝑘⁻¹𝜈_𝑘,(c) for 𝑑=1,…,𝐷−1, compute the smoothing gain matrices 𝑆^𝑧_{𝑘−𝑑,𝑘} and 𝑆^𝑣_{𝑘−𝑑,𝑘} from (4.20) and (4.21), respectively, and the smoothers 𝑧_{𝑘−𝑑/𝑘}=𝑧_{𝑘−𝑑/𝑘−1}+𝑆^𝑧_{𝑘−𝑑,𝑘}Π_𝑘⁻¹𝜈_𝑘 and 𝑣_{𝑘−𝑑/𝑘}=𝑣_{𝑘−𝑑/𝑘−1}+𝑆^𝑣_{𝑘−𝑑,𝑘}Π_𝑘⁻¹𝜈_𝑘. Then, calculate the innovation 𝜈_𝑘+1. To obtain its covariance matrix, Π_𝑘+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 𝑟_𝑘=𝑟_𝑘−1+𝐽_𝑘Π_𝑘⁻¹𝐽^𝑇_𝑘 to calculate 𝐽_𝑘+1 and, from it, the filtering gain 𝑆^𝑧_{𝑘+1,𝑘+1}, which provides the filter 𝑧_{𝑘+1/𝑘+1} and the error covariance matrix 𝑃^𝑧𝑧_{𝑘+1,𝑘+1/𝑘+1}.

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 {𝑧_𝑘;𝑘≥1} with autocovariance function given by 𝐸[𝑧_𝑘𝑧^𝑇_𝑠]=1.025641×0.95^𝑘−𝑠,𝑠≤𝑘, which is factorizable according to Assumption 1 of the model just taking 𝐴_𝑘=1.025641×0.95^𝑘 and 𝐵_𝑘=0.95^−𝑘. For the simulation, the signal is supposed to be generated by an autoregressive model, 𝑧_𝑘+1=0.95𝑧_𝑘+𝑤_𝑘, where {𝑤_𝑘;𝑘≥1} is a zero-mean white Gaussian noise with 𝑉𝑎𝑟[𝑤_𝑘]=0.1, for all 𝑘.

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 is 𝐷=3. For this purpose, we implemented a MATLAB program which simulates the values of the signal, 𝑧_𝑘, the real measurements, 𝑦_𝑘, and the available ones, 𝑦_𝑘, considering different delay probabilities, and provides the filtering and fixed-point smoothing estimates of 𝑧_𝑘, 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, 𝑦₁=𝑦₁; at time 𝑘=2, the real measurement may be one-step randomly delayed; at time 𝑘=3 it may be randomly delayed either by one or two steps and, at any sampling time 𝑘≥4, the available measurement, 𝑦_𝑘, 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, {𝜁¹_𝑘;𝑘≥2}, {𝜁²_𝑘;𝑘≥3}, and {𝜁³_𝑘;𝑘≥4} with constant probabilities 𝑃(𝜁^𝑖_𝑘=1)=𝑞^𝑖,𝑖=1,2,3, and defining the available measurements of the signal as 𝑦₁=𝑦₁,𝑦₂=(1−𝜁¹₂)𝑦₂+𝜁¹₂𝑦₁,𝑦₃=(1−𝜁¹₃)𝑦₃+𝜁¹₃(1−𝜁²₃)𝑦₂+𝜁¹₃𝜁²₃𝑦₁,𝑦_𝑘=(1−𝜁¹_𝑘)𝑦_𝑘+𝜁¹_𝑘(1−𝜁²_𝑘)𝑦_𝑘−1+𝜁¹_𝑘𝜁²_𝑘(1−𝜁³_𝑘)𝑦_𝑘−2+𝜁¹_𝑘𝜁²_𝑘𝜁³_𝑘𝑦_𝑘−3,𝑘≥4. (5.1) Hence, for 𝑑=0,1,2,3, the sequences {𝛾_𝑘^(𝑑);𝑘>𝑑} in the theoretical study are given by 𝛾_𝑘⁽⁰⁾=1−𝜁¹_𝑘with𝑝_𝑘⁽⁰⁾=1−𝑞¹,𝑘>1,𝛾₁⁽⁰⁾=1with𝑝₁⁽⁰⁾𝛾=1,_𝑘(1)=𝜁¹_𝑘(1−𝜁²_𝑘)with𝑝_𝑘⁽¹⁾=𝑞¹(1−𝑞²),𝑘>2,𝛾₂⁽¹⁾=𝜁¹₂with𝑝₂⁽¹⁾=𝑞¹,𝛾_𝑘⁽²⁾=𝜁¹_𝑘𝜁²_𝑘(1−𝜁³_𝑘)with𝑝_𝑘⁽²⁾=𝑞¹𝑞²(1−𝑞³),𝑘>3,𝛾₃⁽²⁾=𝜁¹₃𝜁²₃with𝑝₃⁽²⁾=𝑞¹𝑞²,𝛾_𝑘⁽³⁾=𝜁¹_𝑘𝜁²_𝑘𝜁³_𝑘with𝑝_𝑘⁽³⁾=𝑞¹𝑞²𝑞³,𝑘>3. (5.2)

If 𝜁¹_𝑘=0, then 𝑦_𝑘=𝑦_𝑘 and the 𝑘th measurement is updated; if 𝜁¹_𝑘=1, the 𝑘th measurement is one-step delayed when 𝜁²_𝑘=0 and two-step delayed when 𝜁²_𝑘=1 and 𝜁³_𝑘=0; finally, if 𝜁^𝑖_𝑘=1, for 𝑖=1,2,3, the 𝑘th measurement is three-step delayed.

Note that 𝑞¹=1−𝑝_𝑘⁽⁰⁾ represents the probability of receiving a delayed observation at each sampling time. Moreover, note that 𝑝_𝑘⁽⁰⁾, 𝑝_𝑘⁽¹⁾, and 𝑝_𝑘⁽²⁾ are decreasing functions of 𝑞¹, 𝑞², and 𝑞³, respectively, while 𝑝_𝑘⁽³⁾ is an increasing function of 𝑞^𝑖, 𝑖=1,2,3. Actually, as 𝑞¹ increases, 𝑝_𝑘⁽⁰⁾ decreases, but the delay probabilities 𝑝_𝑘⁽¹⁾, 𝑝_𝑘⁽²⁾ and 𝑝_𝑘⁽³⁾ increase; as 𝑞² increases, 𝑝_𝑘⁽¹⁾ decreases, but the two- and three-step delay probabilities 𝑝_𝑘⁽²⁾ and 𝑝_𝑘⁽³⁾ increase; and, as 𝑞³ increases, 𝑝_𝑘⁽²⁾ decreases, but the three-step delay probability 𝑝_𝑘⁽³⁾ increases.

For 𝑘=1,…,25, Figure 1 displays the processed measurements, 𝑦_𝑘, taking 𝑞^𝑖=0.5, 𝑖=1,2,3; that is, when the probability that the real measurement is used in the estimation is 𝑝_𝑘⁽⁰⁾=1−𝑞¹=0.5 for 𝑘≥2, the one-step delay probability is 𝑝_𝑘⁽¹⁾=0.5𝑝_𝑘⁽⁰⁾=0.25 for 𝑘≥3, and the two- and three-step delay probabilities are both equal 𝑝_𝑘⁽²⁾=𝑝_𝑘⁽³⁾=0.25𝑝_𝑘⁽⁰⁾=0.125 for 𝑘≥4. This figure shows that, in fact, the measured output 𝑦_𝑘 at time 𝑘 can be received on time, delayed or lost in transmission. Actually, the available measurement being processed is(i)updated (𝑦_𝑘=𝑦_𝑘) at 𝑘=1,2,4,5,7,9,11,13,15,17,18,21,22,23,25, (ii)delayed by one sampling period (𝑦_𝑘=𝑦_𝑘−1) at 𝑘=3,8,10,20,24, (iii)delayed by two sampling periods (𝑦_𝑘=𝑦_𝑘−2) at 𝑘=6,12,16,19, (iv)delayed by three sampling periods (𝑦_𝑘=𝑦_𝑘−3) at 𝑘=14,

while 𝑦₃,𝑦₆,𝑦₈,𝑦₁₂,𝑦₁₆,𝑦₂₀, and 𝑦₂₄ 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 𝐷+1=4 consecutive times (for example, 𝑦₂ is rereceived at 𝑘=3, 𝑦₄ is rereceived at 𝑘=6, and 𝑦₁₁ is rereceived at 𝑘=14, among others).

Next, the performance of the estimators is analyzed by investigating the error variances for 𝑞^𝑖=0.5,𝑖=1,2,3. Figure 2 shows the error variances for the filter 𝑧_𝑘/𝑘 and the smoothers 𝑧_𝑘/𝑘+1 and 𝑧_𝑘/𝑘+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, 𝑧_𝑘/𝑘 and 𝑧_𝑘/𝑘+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 for $q^1=q^2=q^3=0.5$.

Figure 3 

Signal $z_k$ as well as filtering and smoothing estimates, ${\widehat{z}}_{k/k}$ and ${\widehat{z}}_{k/k+2}$, when $q^1=q^2=q^3=0.5$.

To analyze the performance of the proposed estimators versus the delay probabilities, the error variances have been calculated for different values of 𝑞^𝑖, 𝑖=1,2,3, which provide different values of the probabilities 𝑝_𝑘^(𝑖), 𝑖=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 𝑘=100 versus 𝑞¹ (for 𝑞²=0.3,0.5,and0.7 and 𝑞³=0.5). From this figure it is deduced that, as 𝑞¹ increases (and, consequently, the nondelay probability 𝑝_𝑘⁽⁰⁾ decreases), the estimation-error variances become greater and, hence, worse estimations are obtained.

Figure 4 

Error variances versus $q^1$, when $q^2=0.3,\mathrm{\hspace{0.17em}\hspace{0.17em}}0.5,\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{and}\mathrm{\hspace{0.17em}\hspace{0.17em}}0.7$ and $q^3=0.5$.

On the other hand, for each value of 𝑞¹ (keeping 𝑞³=0.5 fixed), the error variances become smaller as 𝑞² decreases, which means that the estimations are better. This was expectable since, as 𝑞² decreases, the one-step delay probability increases and the two- and three-step delay probabilities decrease.

Also, as expected, this improvement is more significant as 𝑞¹ 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 at 𝑘=100 versus 𝑞², for 𝑞³=0.3,0.5,and0.7 and 𝑞¹=0.9 (similar results are obtained for other fixed values of the probability of delay, 𝑞¹, but as indicated above, the comparison is clearer when such probability is large). The results are displayed in Figure 5, which shows that, as 𝑞² 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 𝑞², the estimators perform worse with the increasing of 𝑞³, which is reasonable since the three-step delay probability increases with 𝑞³. Nevertheless, for small values of 𝑞², the difference is almost imperceptible since, as 𝑞² becomes close to zero, the one-step delay probability tends to 𝑞¹ (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 𝑘=100 versus 𝑞³ (for 𝑞²=0.3,0.5,and0.7 and 𝑞¹=0.9) are displayed.

Figure 5 

Error variances versus $q^2$, when $q^3=0.3,\mathrm{\hspace{0.17em}\hspace{0.17em}}0.5,\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{and}\mathrm{\hspace{0.17em}\hspace{0.17em}}0.7$ and $q^1=0.9$.

Figure 6 

Error variances versus $q^3$, when $q^2=0.3,\mathrm{\hspace{0.17em}\hspace{0.17em}}0.5,\text{and}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}0.7$ and $q^1=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 constant 𝐷, and the random measurement delays and packet dropouts are described by introducing 𝐷+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 is 𝐷=3.

Acknowledgments

This research is supported by Ministerio de Educación y Ciencia (Grant no. MTM2008-05567) and Junta de Andalucía (Grant no. P07-FQM-02701).