Optimal Fusion Estimation with Multi-Step Random Delays and Losses in Transmission

Abstract

This paper is concerned with the optimal fusion estimation problem in networked stochastic systems with bounded random delays and packet dropouts, which unavoidably occur during the data transmission in the network. The measured outputs from each sensor are perturbed by random parameter matrices and white additive noises, which are cross-correlated between the different sensors. Least-squares fusion linear estimators including filter, predictor and fixed-point smoother, as well as the corresponding estimation error covariance matrices are designed via the innovation analysis approach. The proposed recursive algorithms depend on the delay probabilities at each sampling time, but do not to need to know if a particular measurement is delayed or not. Moreover, the knowledge of the signal evolution model is not required, as the algorithms need only the first and second order moments of the processes involved. Some of the practical situations covered by the proposed system model with random parameter matrices are analyzed and the influence of the delays in the estimation accuracy are examined in a numerical example.

1. Introduction

Over the last few decades, research on the estimation problem for networked stochastic systems has gained considerable attention, due to the undeniable advantages of networked systems, whose applicability is encouraged, among other causes, by the development and advances in communication technology and the growing use of wireless networks. As it is well known, the Kalman filter provides a recursive algorithm for the optimal least-squares estimator in stochastic linear systems, assuming that the system model is exactly known and all of the measurements are instantly updated. The development of sensor networks motivates the necessity of designing new estimation algorithms that integrate the information of all the sensors to achieve a satisfactory performance; thus, using different fusion techniques, the measurements from multiple sensors are combined to obtain more accurate estimators than those obtained when a single sensor is used. In this framework, important extensions of the Kalman filter have been proposed for conventional sensor networks in which the measured outputs of the sensors always contain the actual signal contaminated by additive noises, and the transmissions are carried through perfect connections (see, e.g., [1,2,3] and references therein).

However, in a network environment, usually the standard observation models are not suitable due to the existence of network-induced uncertainties that can occur in both of the sensor measured outputs, and during the data transmission through the network. Accordingly, the consideration of appropriate observation models is vitally important to address the estimation problem in networked systems. Random failures in the transmission of measured data, together with the inaccuracy of the measurement devices, cause often the degradation of the estimator performance in networked systems. In light of these concerns, the estimation problem with one or even several network-induced uncertainties is recently attracting considerable attention, and the design of new fusion estimation algorithms has become an active research topic (see, e.g., [4,5,6,7,8,9,10,11,12,13] and references therein). In addition, some recent advances on the estimation, filtering and fusion for networked systems with network-induced phenomena can be reviewed in [14,15], where a detailed overview of this field is presented.

One of the most common network-induced uncertainties in the measured outputs of the different sensors is the presence of multiplicative noise, due to different reasons, such as interferences or intermittent sensor failures. Specifically, in situations involving random observation losses (see [5]), sensor gain degradation (see [6]), missing or fading measurements (see [13,16], respectively), the sensor observation equations include multiplicative noises. A unified framework to model these random phenomena is provided by the use of random measurement matrices in the sensor observation model. For this reason, the estimation problem in networked systems with random measurement matrices has become a fertile research subject, since this class of systems allow for covering different networked-induced random uncertainties as those mentioned above (see, e.g., [17,18,19,20,21,22,23], and references therein).

In relation to the network-induced uncertainties during the data transmission, it must be indicated that sudden changes in the environment and the unreliability of the communication network, together with the limited bandwidths of communication channels, cause unavoidable random failures during the transmission process. Generally, random communication delays and/or transmission packet dropouts are two essential issues that must be taken into account to model the measurements, which, being available after transmission, will be used for the estimation. Several estimation algorithms have been proposed in multisensor systems considering either transmission delays or packet losses (see, e.g., [24,25]) and also taking into account random delays and packet dropouts simultaneously (see, e.g., [4,23,26]). By using the state augmentation method, systems with random delays and packet dropouts can be transformed into systems with random parameter matrices (see, e.g., [8,9,10,20]). Hence, systems with random parameter measurement matrices also provide an appropriate unified context for modelling these random phenomena in the transmission.

Nevertheless, it must be indicated that the state augmentation method leads to a rise of the computational burden, due to the increase of the state dimension. Actually, in models with more than one or two-step random delays, the computational cost can be excessive and alternative ways to model and address the estimation problem in this class of systems need to be investigated. Recently, a great variety of models have been used to describe the phenomena of multi-step random delays and packet losses during the data transmission in networked systems, and fusion estimation algorithms have been proposed based on different approaches—for example, the recursive matrix equation method in [6], the measurement reorganisation approach in [27], the innovation analysis approach in [28] and the state augmentation approach in [29,30,31]. It should be noted that, in the presence of multi-step random delays and packet losses during the data transmission, many difficulties can arise in the design of the optimal estimators when the state augmentation approach is not used.

In view of the above considerations, this paper is concerned with the optimal fusion estimation problem, in the least-squares linear sense, for sensor networks featuring stochastic uncertainties in the sensor measurements, together with multi-step random delays and packet dropouts during the data transmission. The derivation of the estimation algorithms will be carried out without using the evolution model generating the signal process. The uncertainties in the measured outputs of the different sensors are described by random measurement matrices. The multi-step random delays in the transmissions are modeled by using a collection of Bernoulli sequences with known distributions and different characteristics at each sensor; the exact value of these Bernoulli variables is not required, and only the information about the probability distribution is needed. To the best of the authors’ knowledge, the optimal estimation problem (including prediction, filtering and fixed-point smoothing) has not been investigated for systems involving random measurement matrices and transmission multi-step random delays simultaneously, and, therefore, it constitutes an interesting research challenge. The main contributions of this research can be highlighted as follows: (a) even though our approach, based on covariance information, does not require the signal evolution model, the proposed algorithms are also applicable in situations based on the state-space model (see Remark 1); (b) random measurement matrices are considered in the measured outputs, thus providing a unified framework to address different network-induced phenomena (see Remark 2); (c) besides the stochastic uncertainties in the sensor measurements, simultaneous multi-step random delays and losses with different rates are considered in the data transmission; (d) unlike most papers about multi-step random delays, in which only the filtering problem is considered, we propose recursive algorithms for the prediction, filtering an fixed-point smoothing estimators under the innovation approach, which are computationally very simple and suitable for online applications; and (e) optimal estimators are obtained without using the state augmentation approach, thus reducing the computational cost in comparison with the augmentation method.

The rest of the paper is organized as follows. In Section 2, we present the sensor network and the assumptions under which the optimal linear estimation problem will be addressed. In Section 3, the observation model is rewritten in a compact form and the innovation approach to the least-squares linear estimation problem is formulated. In Section 4, recursive algorithms for the prediction, filtering and fixed-point smoothing estimators are derived. A simulation example is given in Section 5 to show the performance of the proposed estimators. Finally, some conclusions are drawn in Section 6.

Notations. 

The notations used throughout the paper are standard. ${\mathbb{R}}^n$ and ${\mathbb{R}}^{m\times n}$ denote the n-dimensional Euclidean space and the set of all $m\times n$ real matrices, respectively. For a matrix A, $A^T$ and $A^{-1}$ denote its transpose and inverse, respectively. The shorthand $Diag(A_1,\dots ,A_m)$ stands for a block-diagonal matrix whose diagonal matrices are $A_1,\dots ,A_m$. ${\mathbf{1}}_n={(1,\dots ,1)}^T$ denotes the all-ones $n\times 1$-vector and $I_n$ represent the $n\times n$-identity matrix. If the dimensions of matrices are not explicitly stated, they are assumed to be compatible with algebraic operations. The Kronecker and Hadamard product of matrices will be denoted by ⊗ and ∘, respectively. ${\delta }_{k,s}$ denotes the Kronecker delta function. For any $a,b\in \mathbb{R}$, $a\wedge b$ is used to mean the minimum of a and b.

2. Observation Model and Preliminaries

In this paper, the optimal fusion estimation problem of multidimensional discrete-time random signals from measurements obtained by a sensor network is addressed. At each sensor, the measured outputs are perturbed by random parameter matrices and white additive noises that are cross-correlated at the same sampling time between the different sensors. The estimation is performed in a processing centre connected to all sensors, where the complete set of sensor data is combined, but due to eventual communication failures, congestion or other causes, random delays and packet dropouts are unavoidable during the transmission process. To reduce the effect of such delays and packet dropouts without overloading the network traffic, each sensor measurement transmits a fixed number of consecutive sampling times and, when several packets arrive at the same time, the receiver discards the oldest ones, so that only one measured output is processed for the estimation at each sampling time.

2.1. Signal Process

The optimal estimators will be obtained using the least-squares (LS) criterion and without requiring the evolution model generating the signal process. Actually, the proposed estimation algorithms, based on covariance information, only need the mean vectors and covariance functions of the processes involved, and the only requirement will be that the signal covariance function must be factorizable according to the following assumption:

(A1)

The $n_x$-dimensional signal process $\{x_k;k\ge 1\}$ has zero mean and its autocovariance function is expressed in a separable form, $E\left[x_k{x}_s^T\right]=A_k{B}_s^T,s\le k,$ where $A_k,B_s\in {\mathbb{R}}^{n_x\times n}$ are known matrices.

Remark 1 (on assumption (A1)).

The estimation problems based on the state-space model require new estimation algorithms when the signal evolution model is modified; therefore, the algorithms designed for stationary signals driven by $x_{k+1}=\mathsf{\Phi }{x}_k+{\xi }_k$ cannot be applied for non-stationary signals generated by $x_{k+1}={\mathsf{\Phi }}_k{x}_k+{\xi }_k$, and these, in turn, cannot be used in uncertain systems where $x_{k+1}=({\mathsf{\Phi }}_k+{ϵ}_k{\widehat{\mathsf{\Phi }}}_k)x_k+{\xi }_k$. A great advantage of assumption (A1) is that it covers situations in which the signal evolution model is known, for both stationary and non-stationary signals (see, e.g., [23]). In addition, in uncertain systems with state-dependent multiplicative noise, as those considered in [6,32], the signal covariance function is factorizable, as it is shown in Section 5. Hence, assumption (A1) on the signal autocovariance function provides a unified context to deal with different situations based on the state-space model, avoiding the derivation of specific algorithms for each of them.

2.2. Multisensor Observation Model

Assuming that there are m different sensors, the measured outputs before transmission, $z_k^{\left(i\right)}\in {\mathbb{R}}^{n_z}$, are described by the following observation model:

$$z_k^{\left(i\right)}=H_k^{\left(i\right)}{x}_k+v_k^{\left(i\right)},k\ge 1;i=1,\dots ,m,$$

(1)

where the measurement matrices, $H_k^{\left(i\right)}$, and the noise vectors, $v_k^{\left(i\right)}$, satisfy the following assumptions:

(A2)

$\left\{H_k^{\left(i\right)};k\ge 1\right\}$, $i=1,\dots ,m$, are independent sequences of independent random parameter matrices, whose entries have known means and second-order moments; we will denote ${\overline{H}}_k^{\left(i\right)}\equiv E\left[H_k^{\left(i\right)}\right],k\ge 1$.

(A3)

$\left\{v_k^{\left(i\right)};k\ge 1\right\},i=1,\dots ,m$, are white noise sequences with zero mean and known second-order moments, satisfying $E\left[v_k^{\left(i\right)}{v}_s^{\left(j\right)T}\right]=R_k^{\left(ij\right)}{\delta }_{k,s},i,j=1,\dots ,m.$

Remark 2 (on assumption (A2)).

Usually, in network environments, the measurements are subject to different network-induced random phenomena and new estimation algorithms must be designed to incorporate the effects of these random uncertainties. For example, in systems with stochastic sensor gain degradation or missing measurements as those considered in [6,7], respectively, or in networked systems involving stochastic multiplicative noises in the state and measurement equations (see, e.g., [31,32]), new estimation algorithms are proposed since the conditions necessary to implement the conventional ones are not met. The aforementioned systems are particular cases of systems with random measurement matrices, and, hence, assumption (A2) allows for designing a unique estimation algorithm, which is suitable to address all of these situations involving random uncertainties. In addition, based on an augmentation approach, random measurement matrices can be used to model the measured outputs of sensor networks with random delays and packet dropouts (see, e.g., [8,9,10,20]). Therefore, assumption (A2) provides a unified framework to deal with a great variety of network-induced random phenomena as those mentioned above.

2.3. Measurement Model with Transmission Random Delays and Packet Losses

Assuming that the maximum time delay is D, the measured output of the i-th sensor at time r, $z_r^{\left(i\right)}$, is transmitted during the sampling times $r,r+1,\cdots ,r+D$, but, at each sampling time $k>D$, only one of the measurements $z_{k-D}^{\left(i\right)},\dots ,z_k^{\left(i\right)}$ is processed. Consequently, at any time $k>D$, the measurement processed can either arrive on time or be delayed by $d=1,\dots ,D$ sampling periods, while at any time $k\le D$, the measurement processed can be delayed only by $d=1,\dots ,k-1$ sampling periods, since only $z_1^{\left(i\right)},\dots ,z_k^{\left(i\right)}$ are available. Assuming, moreover, that the transmissions are perturbed by additive noises, the measurements received at the processing centre, impaired by random delays and packet losses, can be described by the following model:

$${\displaystyle y_k^{\left(i\right)}=\sum _{d=0}^{(k-1)\wedge D}{\gamma }_{d,k}^{\left(i\right)}{z}_{k-d}^{\left(i\right)}+w_k^{\left(i\right)},k\ge 1,}$$

(2)

where the following assumptions on the random variables modelling the delays, ${\gamma }_{d,k}^{\left(i\right)}$, and the transmision noise, $w_k^{\left(i\right)}$, are required:

(A4)

For each $d=0,1,\dots ,D$, $\left\{{\gamma }_{d,k}^{\left(i\right)};k>d\right\}$, $i=1,\dots ,m,$ are independent sequences of independent Bernoulli random variables with $P[{\gamma }_{d,k}^{\left(i\right)}=1]={\overline{\gamma }}_{d,k}^{\left(i\right)}$ and ${\displaystyle \sum _{d=0}^{(k-1)\wedge D}{\gamma }_{d,k}^{\left(i\right)}\le 1},k\ge 1.$

(A5)

$\left\{w_k^{\left(i\right)};k\ge 1\right\},i=1,\dots ,m$, are white noise sequences with zero mean and known second-order moments, satisfying $E\left[w_k^{\left(i\right)}{w}_s^{\left(j\right)T}\right]=Q_k^{\left(ij\right)}{\delta }_{k,s},i,j=1,\dots ,m.$

Remark 3 (on assumption (A4)).

For $i=1,\dots ,m,$ when ${\gamma }_{0,k}^{\left(i\right)}=1$, the transmission of the i-th sensor is perfect and neither delay nor loss occurs at time k; that is, with probability ${\overline{\gamma }}_{0,k}^{\left(i\right)}$, the k-th measurement of the i-th sensor is received and processed on time. Since ${\displaystyle \sum _{d=0}^{(k-1)\wedge D}{\gamma }_{d,k}^{\left(i\right)}\le 1}$, if ${\gamma }_{0,k}^{\left(i\right)}=0$, there then exists at most one $d=1,\dots ,D$, such that ${\gamma }_{d,k}^{\left(i\right)}=1$. If there exists d such that ${\gamma }_{d,k}^{\left(i\right)}=1$ (which occurs with probability ${\overline{\gamma }}_{d,k}^{\left(i\right)}$), then the measurement is delayed by d sampling periods. Otherwise, ${\gamma }_{d,k}^{\left(i\right)}=0$ for all d, meaning that the measurement gets lost during the transmission at time k with probability ${\displaystyle 1-\sum _{d=0}^{(k-1)\wedge D}{\overline{\gamma }}_{d,k}^{\left(i\right)}}$.

Finally, the following independence hypothesis is assumed:

(A6)

For $i=1,\dots ,m$ and $d=0,1,\dots ,D$, the processes $\{x_k;k\ge 1\}$, $\{H_k^{\left(i\right)};k\ge 1\}$, $\{v_k^{\left(i\right)};k\ge 1\}$, $\{w_k^{\left(i\right)};k\ge 1\}$ and $\left\{{\gamma }_{d,k}^{\left(i\right)};k>d\right\}$ are mutually independent.

3. Problem Statement

Given the observation Equations (1) and (2) with random measurement matrices and transmission random delays and packet dropouts, our purpose is to find the LS linear estimator, ${\widehat{x}}_{k/L}$, of the signal $x_k$ based on the observations from the different sensors $\{y_1^{\left(i\right)},\dots ,y_L^{\left(i\right)},i=1,\dots ,m\}$. Specifically, our aim is to obtain recursive algorithms for the predictor ($L<k$), filter $(L=k)$ and fixed-point smoother (k fixed and $L>k$).

3.1. Stacked Observation Model

Since the measurements coming from the different sensors are all gathered and jointly processed at each sampling time k, we will consider the vector constituted by the meaurements from all sensors, $y_k={\left(y_k^{\left(1\right)T},\dots ,y_k^{\left(m\right)T}\right)}^T$. More specifically, the observation Equations (1) and (2) of all sensors are combined, yielding the following stacked observation model:

$$\begin{array}{c}{z}_k=H_k{x}_k+v_k,k\ge 1,\hfill \\ {\displaystyle y_k=\sum _{d=0}^{(k-1)\wedge D}{\mathsf{\Gamma }}_{d,k}{z}_{k-d}+w_k,k\ge 1,}\hfill \end{array}$$

(3)

where $z_k={\left(z_k^{\left(1\right)T},\dots ,z_k^{\left(m\right)T}\right)}^T$, $H_k={\left(H_k^{\left(1\right)T},\dots ,H_k^{\left(m\right)T}\right)}^T$, $v_k={\left(v_k^{\left(1\right)T},\dots ,v_k^{\left(m\right)T}\right)}^T$, $w_k={\left(w_k^{\left(1\right)T},\dots ,w_k^{\left(m\right)T}\right)}^T$ and ${\mathsf{\Gamma }}_{d,k}=Diag\left({\gamma }_{d,k}^{\left(1\right)},\dots ,{\gamma }_{d,k}^{\left(m\right)}\right)\otimes I_{n_z}$.

Hence, the problem is to obtain the LS linear estimator of the signal, $x_k$, based on the measurements $\left\{y_1,\dots ,y_L\right\}$, given in the observation Equation (3). Next, we present the statistical properties of the processes involved in Equation (3), from which the LS estimation algorithms of the signal will be derived; these properties are easily inferred from the assumptions (A1)–(A6).

(P1)

$\left\{H_k;k\ge 1\right\}$ is a sequence of independent random parameter matrices with known means, ${\overline{H}}_k\equiv E\left[H_k\right]={\left({\overline{H}}_k^{\left(1\right)T},\dots ,{\overline{H}}_k^{\left(m\right)T}\right)}^T$, and

$$\begin{array}{cc}\hfill E\left[H_k{x}_k{x}_s^T{H}_s^T\right]=E\left[H_k{A}_k{B}_s^T{H}_s^T\right]={\left(E\left[H_k^{\left(i\right)}{A}_k{B}_s^T{H}_s^{\left(j\right)T}\right]\right)}_{i,j=1\dots ,m}\hfill & \hfill ,s\le k,\hfill \end{array}$$

where $E\left[H_k^{\left(i\right)}{A}_k{B}_s^T{H}_s^{\left(j\right)T}\right]={\overline{H}}_k^{\left(i\right)}{A}_k{B}_s^T{\overline{H}}_s^{\left(j\right)T}$, for $j\ne i$ or $s\ne k$, and the entries of $E\left[H_k^{\left(i\right)}{A}_k{B}_k^T{H}_k^{\left(i\right)T}\right]$ are computed as follows:

$${\left(E\left[H_k^{\left(i\right)}{A}_k{B}_k^T{H}_k^{\left(i\right)T}\right]\right)}_{pq}=\sum _{a=1}^{n_x}\sum _{b=1}^{n_x}E\left[h_{pa}^{\left(i\right)}\left(k\right)h_{qb}^{\left(i\right)}\left(k\right)\right]{\left(A_k{B}_k^T\right)}_{ab},p,q=1,\dots ,n_z,$$

where $h_{pq}^{\left(i\right)}\left(k\right)$ denotes the $(p,q)-$entry of the matrix $H_k^{\left(i\right)}$.

(P2)

The noises $\left\{v_k;k\ge 1\right\}$ and $\left\{w_k;k\ge 1\right\}$ are zero-mean sequences with known second-order moments given by the matrices $R_k\equiv {\left(R_k^{\left(ij\right)}\right)}_{i,j=1,\dots ,m}$ and $Q_k\equiv {\left(Q_k^{\left(ij\right)}\right)}_{i,j=1,\dots ,m}$.

(P3)

$\{{\mathsf{\Gamma }}_{d,k};k>d\}$, $d=0,1,\dots ,D$, are sequences of independent random matrices with known means, ${\overline{\mathsf{\Gamma }}}_{d,k}\equiv E\left[{\mathsf{\Gamma }}_{d,k}\right]=Diag\left({\overline{\gamma }}_{d,k}^{\left(1\right)},\dots ,{\overline{\gamma }}_{d,k}^{\left(m\right)}\right)\otimes I_{n_z}$, and if we denote ${\gamma }_{d,k}={\left({\gamma }_{d,k}^{\left(1\right)},\dots ,{\gamma }_{d,k}^{\left(m\right)}\right)}^T\otimes {\mathbf{1}}_{n_z}$ and ${\overline{\gamma }}_{d,k}=E\left[{\gamma }_{d,k}\right]$, the covariance matrices ${\mathsf{\Sigma }}_{d,d^{\prime }}^{{\gamma }_k}\equiv E\left[({\gamma }_{d,k}-{\overline{\gamma }}_{d,k}){({\gamma }_{d^{\prime },k}-{\overline{\gamma }}_{d^{\prime },k})}^T\right]$, for $d,d^{\prime }=0,1,\dots ,D$, are also known matrices. Specifically,

$${\mathsf{\Sigma }}_{d,d^{\prime }}^{{\gamma }_k}=Diag\left(Cov\left({\gamma }_{d,k}^{\left(1\right)},{\gamma }_{d^{\prime },k}^{\left(1\right)}\right),\dots ,Cov\left({\gamma }_{d,k}^{\left(m\right)},{\gamma }_{d^{\prime },k}^{\left(m\right)}\right)\right)\otimes \left({\mathbf{1}}_{n_z}{\mathbf{1}}_{n_z}^T\right),$$

(4)

with $Cov\left({\gamma }_{d,k}^{\left(i\right)},{\gamma }_{d^{\prime },k}^{\left(i\right)}\right)=\left\{\begin{array}{cc}{\overline{\gamma }}_{d,k}^{\left(i\right)}(1-{\overline{\gamma }}_{d,k}^{\left(i\right)}),\hfill & d^{\prime }=d,\hfill \\ -{\overline{\gamma }}_{d,k}^{\left(i\right)}{\overline{\gamma }}_{d^{\prime },k}^{\left(i\right)},\hfill & d^{\prime }\ne d.\hfill \end{array}\right.$ Moreover, for any deterministic matrix S, the Hadamard product properties guarantee that

$$E\left[\left({\mathsf{\Gamma }}_{d,k}-{\overline{\mathsf{\Gamma }}}_{d,k}\right)S\left({\mathsf{\Gamma }}_{d^{\prime },k}-{\overline{\mathsf{\Gamma }}}_{d^{\prime },k}\right)\right]={\mathsf{\Sigma }}_{d,d^{\prime }}^{{\gamma }_k}\circ S.$$

(P4)

For $d=0,1,\dots ,D$, the signal, $\{x_k;k\ge 1\}$, and the processes $\{H_k;k\ge 1\}$, $\{v_k;k\ge 1\}$, $\{w_k;k\ge 1\}$ and $\{{\mathsf{\Gamma }}_{d,k};k>d\}$ are mutually independent.

Remark 4 (on the observation covariance matrices).

From the previous properties, it is clear that the observation process $\{z_k;k\ge 1\}$ is a zero-mean sequence whose covariance function, ${\mathsf{\Sigma }}_{k,s}^z\equiv E\left[z_k{z}_s^T\right]$, is obtained by the following expression:

$${\mathsf{\Sigma }}_{k,s}^z=E\left[H_k{A}_k{B}_s^T{H}_s^T\right]+R_k{\delta }_{k,s},s\le k,$$

(5)

where $E\left[H_k{A}_k{B}_s^T{H}_s^T\right]$ and $R_k$ are calculated as it is indicated in properties (P1) and (P2), respectively.

3.2. Innovation Approach to the LS Linear Estimation Problem

The proposed covariance-based recursive algorithms for the LS linear prediction, filtering and fixed-point smoothing estimators will be derived by an innovation approach. This approach consists of transforming the observation process $\{y_k;k\ge 1\}$ into an equivalent one of orthogonal vectors called an innovation process, which will be denoted $\{{\mu }_k;k\ge 1\}$ and defined by ${\mu }_k=y_k-{\widehat{y}}_{k/k-1}$, where ${\widehat{y}}_{k/k-1}$ is the orthogonal projection of $y_k$ onto the linear space spanned by $\left\{{\mu }_1,\dots ,{\mu }_{k-1}\right\}$. Since both processes span the same linear subspace, the LS linear estimator of any random vector ${\alpha }_k$ based on the observations $\left\{y_1,\dots ,y_N\right\}$, denoted by ${\widehat{\alpha }}_{k/N}$, is equal to that based on the innovations $\left\{{\mu }_1,\dots ,{\mu }_N\right\}$, and, denoting ${\mathsf{\Pi }}_h=E\left[{\mu }_h{\mu }_h^T\right]$, the following general expression for the LS linear estimators of ${\alpha }_k$ is obtained

$${\widehat{\alpha }}_{k/N}=\sum _{h=1}^{N}E\left[{\alpha }_k{\mu }_h^T\right]{\mathsf{\Pi }}_h^{-1}{\mu }_h.$$

(6)

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

Innovation ${\mu }_L$ and Covariance Matrix ${\mathsf{\Pi }}_L$. Applying orthogonal projections in Equation (3), it is clear that the innovation ${\mu }_L$ is given by

$${\displaystyle {\mu }_L=y_L-\sum _{d=0}^{(L-1)\wedge D}{\overline{\mathsf{\Gamma }}}_{d,L}{\widehat{z}}_{L-d/L-1},L\ge 2;{\mu }_1=y_1,}$$

(7)

so it is necessary to obtain the one-stage predictor ${\widehat{z}}_{L/L-1}$ and the estimators ${\widehat{z}}_{L-d/L-1}$, for $d=1,\dots ,(L-1)\wedge D$, of the observation process.

In order to obtain the covariance matrix ${\mathsf{\Pi }}_L$, we use Equation (3) to express the innovations as:

$${\displaystyle {\mu }_L=\sum _{d=0}^{(L-1)\wedge D}\left(({\mathsf{\Gamma }}_{d,L}-{\overline{\mathsf{\Gamma }}}_{d,L})z_{L-d}+{\overline{\mathsf{\Gamma }}}_{d,L}(z_{L-d}-{\widehat{z}}_{L-d/L-1})\right)+w_L,L\ge 2}$$

(8)

and, taking into account that

$$E\left[({\mathsf{\Gamma }}_{d,L}-{\overline{\mathsf{\Gamma }}}_{d,L})z_{L-d}{(z_{L-d^{\prime }}-{\widehat{z}}_{L-d^{\prime }/L-1})}^T{\overline{\mathsf{\Gamma }}}_{d^{\prime },L}\right]=0,\forall d,d^{\prime },$$

we have

$$\begin{array}{c}{\displaystyle {\mathsf{\Pi }}_L=\sum _{d,d^{\prime }=0}^{(L-1)\wedge D}\left({\mathsf{\Sigma }}_{d,d^{\prime }}^{{\gamma }_L}\circ {\mathsf{\Sigma }}_{L-d,L-d^{\prime }}^z+{\overline{\mathsf{\Gamma }}}_{d,L}{P}_{L-d,L-d^{\prime }/L-1}^z{\overline{\mathsf{\Gamma }}}_{d^{\prime },L}\right)+Q_L,L\ge 2;}\hfill \\ {\mathsf{\Pi }}_1=\left({\mathsf{\Sigma }}_{0,0}^{{\gamma }_1}+{\overline{\gamma }}_{0,1}{\overline{\gamma }}_{0,1}^T\right)\circ {\mathsf{\Sigma }}_{1,1}^z+Q_1,\hfill \end{array}$$

(9)

where the matrices ${\mathsf{\Sigma }}_{d,d^{\prime }}^{{\gamma }_L}$ and ${\mathsf{\Sigma }}_{L-d,L-d^{\prime }}^z$ are given in the Equations (4) and (5), respectively, and $P_{L-d,L-d^{\prime }/L-1}^z\equiv E\left[(z_{L-d}-{\widehat{z}}_{L-d/L-1}){(z_{L-d^{\prime }}-{\widehat{z}}_{L-d^{\prime }/L-1})}^T\right]$.

4. Least-Squares Linear Signal Estimators

In this section, we derive recursive algorithms for the LS linear estimators, ${\widehat{x}}_{k/L},k\ge 1,$ of the signal $x_k$ based on the observations $\left\{y_1,\dots ,y_L\right\}$ given in Equation (3); namely, a prediction and filtering algorithm ($L\le k$) and a fixed-point smoothing algorithm (k fixed and $L>k$) are designed.

4.1. Signal Predictor and Filter ${\widehat{x}}_{k/L},L\le k$

From the general expression given in Equation (6), to obtain the LS linear estimators ${\displaystyle {\widehat{x}}_{k/L}=\sum _{h=1}^{L}E\left[x_k{\mu }_h^T\right]{\mathsf{\Pi }}_h^{-1}{\mu }_h}$, $L\le k$, it is necessary to calculate the coefficients

$${\mathcal{X}}_{k,h}\equiv E\left[x_k{\mu }_h^T\right]=E\left[x_k{y}_h^T\right]-E\left[x_k{\widehat{y}}_{h/h-1}^T\right],h\le k.$$

• On the one hand, using Equation (3) together with the independence hypotheses and assumption (A1) on the signal covariance factorization, it is clear that

  $${\displaystyle E\left[x_k{y}_h^T\right]=\sum _{d=0}^{(h-1)\wedge D}E\left[x_k{(H_{h-d}{x}_{h-d}+v_{h-d})}^T\right]{\overline{\mathsf{\Gamma }}}_{d,h}=A_k\sum _{d=0}^{(h-1)\wedge D}{B}_{h-d}^T{\overline{H}}_{h-d}^T{\overline{\mathsf{\Gamma }}}_{d,h},h\le k.}$$
• On the other hand, since ${\displaystyle {\widehat{y}}_{h/h-1}=\sum _{d=0}^{(h-1)\wedge D}{\overline{\mathsf{\Gamma }}}_{d,h}{\widehat{z}}_{h-d/h-1},h\ge 2,}$ and taking into account that, from Equation (6), ${\displaystyle {\widehat{z}}_{h-d/h-1}=\sum _{j=1}^{h-1}{\mathcal{Z}}_{h-d,j}{\mathsf{\Pi }}_j^{-1}{\mu }_j}$ with ${\mathcal{Z}}_{h-d,j}=E\left[z_{h-d}{\mu }_j^T\right]$, the following identity holds

  $${\displaystyle E\left[x_k{\widehat{y}}_{h/h-1}^T\right]=\sum _{d=0}^{(h-1)\wedge D}\left({\displaystyle \sum _{j=1}^{h-1}{\mathcal{X}}_{k,j}{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{h-d,j}^T}\right){\overline{\mathsf{\Gamma }}}_{d,h}.}$$

Therefore, it is easy to check that ${\mathcal{X}}_{k,h}=A_k{\mathcal{E}}_h,1\le h\le k,$ where ${\mathcal{E}}_h$ is a function satisfying that

$$\begin{array}{c}{\displaystyle {\mathcal{E}}_h=\sum _{d=0}^{(h-1)\wedge D}{B}_{h-d}^T{\overline{H}}_{h-d}^T{\overline{\mathsf{\Gamma }}}_{d,h}-\sum _{d=0}^{(h-1)\wedge D}\left({\displaystyle \sum _{j=1}^{h-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{h-d,j}^T}\right){\overline{\mathsf{\Gamma }}}_{d,h},h\ge 2,}\hfill \\ {\mathcal{E}}_1=B_1^T{\overline{H}}_1^T{\overline{\mathsf{\Gamma }}}_{0,1}.\hfill \end{array}$$

(10)

Hence, it is clear that the signal prediction and filtering estimators can be expressed as

$${\widehat{x}}_{k/L}=A_k{e}_L,L\le k,k\ge 1,$$

(11)

where the vectors $e_L$ are defined by ${\displaystyle e_L=\sum _{h=1}^L{\mathcal{E}}_h{\mathsf{\Pi }}_h^{-1}{\mu }_h}$, for $L\ge 1$, with $e_0=0$, thus obeying the recursive relation

$$e_L=e_{L-1}+{\mathcal{E}}_L{\mathsf{\Pi }}_L^{-1}{\mu }_L,L\ge 1;e_0=0.$$

(12)

Matrices ${\mathcal{E}}_L$. Taking into account the above relation, an expression for ${\mathcal{E}}_L,L\ge 1,$ must be derived. For this purpose, Equation (10) is rewritten for $h=L$ as

$${\displaystyle {\mathcal{E}}_L=\sum _{d=0}^{(L-1)\wedge D}{B}_{L-d}^T{\overline{H}}_{L-d}^T{\overline{\mathsf{\Gamma }}}_{d,L}-\sum _{d=0}^{(L-1)\wedge D}\left({\displaystyle \sum _{j=1}^{L-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{L-d,j}^T}\right){\overline{\mathsf{\Gamma }}}_{d,L},L\ge 2,}$$

and we examine the cases $d=0$ and $d\ge 1$ separately:

−

For $d=0$, using Equation (3), it holds that ${\mathcal{Z}}_{L,j}={\overline{H}}_L{\mathcal{X}}_{L,j}={\overline{H}}_L{A}_L{\mathcal{E}}_j$, for $j<L$, and, by denoting ${\displaystyle {\mathsf{\Sigma }}_L^e=\sum _{h=1}^L{\mathcal{E}}_h{\mathsf{\Pi }}_h^{-1}{\mathcal{E}}_h}$, $L\ge 1$, we obtain that

$$\left({\displaystyle \sum _{j=1}^{L-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{L,j}^T}\right){\overline{\mathsf{\Gamma }}}_{0,L}=\left({\displaystyle \sum _{j=1}^{L-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{E}}_j^T}\right)A_L^T{\overline{H}}_L^T{\overline{\mathsf{\Gamma }}}_{0,L}={\mathsf{\Sigma }}_{L-1}^e{A}_L^T{\overline{H}}_L^T{\overline{\mathsf{\Gamma }}}_{0,L}.$$

−

For $d\ge 1$, since ${\mathcal{Z}}_{L-d,j}={\overline{H}}_{L-d}{A}_{L-d}{\mathcal{E}}_j$, for $j<L-d$, we can see that

$$\left({\displaystyle \sum _{j=1}^{L-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{L-d,j}^T}\right){\overline{\mathsf{\Gamma }}}_{d,L}={\mathsf{\Sigma }}_{L-d-1}^e{A}_{L-d}^T{\overline{H}}_{L-d}^T{\overline{\mathsf{\Gamma }}}_{d,L}+\left({\displaystyle \sum _{j=L-d}^{L-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{L-d,j}^T}\right){\overline{\mathsf{\Gamma }}}_{d,L}.$$

By substituting the above sums in ${\mathcal{E}}_L$, it is deduced that

$$\begin{array}{c}{\displaystyle {\mathcal{E}}_L=\sum _{d=0}^{(L-1)\wedge D}{\left(B_{L-d}-A_{L-d}{\mathsf{\Sigma }}_{L-d-1}^e\right)}^T{\overline{H}}_{L-d}^T{\overline{\mathsf{\Gamma }}}_{d,L}-\sum _{d=1}^{(L-1)\wedge D}\left({\displaystyle \sum _{j=L-d}^{L-1}{\mathcal{E}}_j{\mathsf{\Pi }}_j^{-1}{\mathcal{Z}}_{L-d,j}^T}\right){\overline{\mathsf{\Gamma }}}_{d,L},L\ge 2;}\hfill \\ {\mathcal{E}}_1=B_1^T{\overline{H}}_1^T{\overline{\mathsf{\Gamma }}}_{0,1},\hfill \end{array}$$

(13)

where the matrices ${\mathcal{Z}}_{L-d,j},j\ge L-d$, will be obtained in the next subsection, as they correspond to the observation smoothing estimators, and the matrices ${\mathsf{\Sigma }}_L^e$ are recursively obtained by

$${\mathsf{\Sigma }}_L^e={\mathsf{\Sigma }}_{L-1}^e+{\mathcal{E}}_L{\mathsf{\Pi }}_L^{-1}{\mathcal{E}}_L^T,L\ge 1;{\mathsf{\Sigma }}_0^e=0.$$

(14)

Finally, from assumption (A1) and since the estimation errors are orthogonal to the estimators, we have that the error covariance matrices, $P_{k/L}^x\equiv E\left[(x_k-{\widehat{x}}_{k/L}){(x_k-{\widehat{x}}_{k/L})}^T\right]$, are given by

$$P_{k/L}^x=A_k{\left(B_k-A_k{\mathsf{\Sigma }}_L^e\right)}^T,L\le k,k\ge 1.$$

(15)

4.2. Estimators of the Observations ${\widehat{z}}_{k/L},k\ge 1$

As it has been already indicated, the Equation (7) require obtaining the observation estimators (predictor, filter and smoother). From the general expression for the estimators given in Equation (6), we have that ${\displaystyle {\widehat{z}}_{k/L}=\sum _{j=1}^L{\mathcal{Z}}_{k,j}{\mathsf{\Pi }}_j^{-1}{\mu }_j}$, with ${\mathcal{Z}}_{k,j}=E\left[z_k{\mu }_j^T\right]$. Next, recursive expressions will be derived separately for $L<k$ (predictors) and $L\ge k$ (filter and smoothers).

Observation Prediction Estimators. Since ${\mathcal{Z}}_{k,j}={\overline{H}}_k{A}_k{\mathcal{E}}_j$, for $j<k$, we have that the prediction estimators of the observations are given by

$${\widehat{z}}_{k/L}={\overline{H}}_k{A}_L{e}_L,L<k,k\ge 1.$$

(16)

Observation Filtering and Fixed-Point Smoothing Estimators. Clearly, the filter and fixed-point smoothers of the observations are obtained by the following recursive expression:

$${\widehat{z}}_{k/L}={\widehat{z}}_{k/L-1}+{\mathcal{Z}}_{k,L}{\mathsf{\Pi }}_L^{-1}{\mu }_L,L\ge k,k\ge 1,$$

(17)

with initial condition given by the one-stage predictor ${\widehat{z}}_{k/k-1}={\overline{H}}_k{A}_k{e}_{k-1}$.

Hence, the matrices ${\mathcal{Z}}_{k,L}$ must be calculated for $L\ge k$. Since the innovation is a white process, $E\left[{\widehat{z}}_{k/L-1}{\mu }_L^T\right]=0$ and hence ${\mathcal{Z}}_{k,L}=E\left[z_k{\mu }_L^T\right]=E\left[(z_k-{\widehat{z}}_{k/L-1}){\mu }_L^T\right]$. Now using Equation (8) for ${\mu }_L$ and, taking into account that $E\left[(z_k-{\widehat{z}}_{k/L-1})z_{L-d}^T({\mathsf{\Gamma }}_{d,L}-{\overline{\mathsf{\Gamma }}}_{d,L})\right]=0,\forall d,$ we have

$${\displaystyle {\mathcal{Z}}_{k,L}=\sum _{d=0}^{(L-1)\wedge D}{P}_{k,L-d/L-1}^z{\overline{\mathsf{\Gamma }}}_{d,L},L\ge k,}$$

(18)

where $P_{k,L-d/L-1}^z\equiv E\left[(z_k-{\widehat{z}}_{k/L-1}){(z_{L-d}-{\widehat{z}}_{L-d/L-1})}^T\right]$.

Consequently, the error covariance matrices $P_{k,h/m}^z=E\left[(z_k-{\widehat{z}}_{k/m}){(z_h-{\widehat{z}}_{h/m})}^T\right]$ must be derived, for which the following two cases are analyzed separately:

*

For $m\ge k\wedge h$, using Equation (17) and taking into account that ${\mathcal{Z}}_{k,m}=E\left[(z_k-{\widehat{z}}_{k/m-1}){\mu }_m^T\right]$, it is easy to see that

$$P_{k,h/m}^z=P_{k,h/m-1}^z-{\mathcal{Z}}_{k,m}{\mathsf{\Pi }}_m^{-1}{\mathcal{Z}}_{h,m}^T,m\ge k\wedge h.$$

(19)

*

For $m<h\le k$, using Equation (16), assumption (A1) and the orthogonality between the estimation errors and the estimators, we obtain

$$P_{k,h/m}^z=A_k{B}_h^T-{\overline{H}}_k{A}_k{\mathsf{\Sigma }}_m^e{A}_h^T{\overline{H}}_h^T,m<h\le k.$$

(20)

4.3. Signal Fixed-Point Smoother ${\widehat{x}}_{k/L},L>k$

Starting with the filter, ${\widehat{x}}_{k/k}$, and the filtering error covariance matrix, $P_{k/k}^x$, it is clear that the signal fixed-point smoother ${\widehat{x}}_{k/L}$, $L>k$, and the corresponding error covariance matrix, $P_{k/L}^x\equiv E\left[(x_k-{\widehat{x}}_{k/L}){(x_k-{\widehat{x}}_{k/L})}^T\right]$, are obtained by

$$\begin{array}{c}{\widehat{x}}_{k/L}={\widehat{x}}_{k/L-1}+{\mathcal{X}}_{k,L}{\mathsf{\Pi }}_L^{-1}{\mu }_L,L>k,k\ge 1,\hfill \\ P_{k/L}^x=P_{k/L-1}^x-{\mathcal{X}}_{k,L}{\mathsf{\Pi }}_L^{-1}{\mathcal{X}}_{k,L}^T,L>k,k\ge 1.\hfill \end{array}$$

(21)

An analogous reasoning to that of Equation (18) leads to the following expression for the matrices ${\mathcal{X}}_{k,L}$:

$${\displaystyle {\mathcal{X}}_{k,L}=\sum _{d=0}^{(L-1)\wedge D}{P}_{k,L-d/L-1}^{xz}{\overline{\mathsf{\Gamma }}}_{d,L},L>k,}$$

(22)

where $P_{k,L-d/L-1}^{xz}\equiv E\left[(x_k-{\widehat{x}}_{k/L-1}){(z_{L-d}-{\widehat{z}}_{L-d/L-1})}^T\right]$.

The derivation of the error cross-covariance matrices $P_{k,h/m}^{xz}=E\left[(x_k-{\widehat{x}}_{k/m}){(z_h-{\widehat{z}}_{h/m})}^T\right]$ is similar to that of the matrices $P_{k,h/m}^z$, and they are given by

$$\begin{array}{cc}{P}_{k,h/m}^{xz}=P_{k,h/m-1}^{xz}-{\mathcal{X}}_{k,m}{\mathsf{\Pi }}_m^{-1}{\mathcal{Z}}_{h,m}^T,\hfill & m\ge k\wedge h,\hfill \\ P_{k,h/m}^{xz}=A_k{(B_h-A_h{\mathsf{\Sigma }}_m^e)}^T{H}_h^T,\hfill & m<h\le k,\hfill \\ P_{k,h/m}^{xz}=(B_k-A_k{\mathsf{\Sigma }}_m^e)A_h^T{H}_h^T,\hfill & m<k\le h,\hfill \end{array}$$

(23)

where ${\mathcal{X}}_{k,m}=A_k{\mathcal{E}}_m$, for $m\le k$, and ${\mathcal{Z}}_{h,m}={\overline{H}}_h{A}_m{\mathcal{E}}_m$, for $m<h$, and otherwise these matrices are given by Equations (22) and (18), respectively.

4.4. Recursive Algorithms: Computational Procedure

The computational procedure of the proposed prediction, filtering and fixed-point smoothing algorithms can be summarized as follows:

(1)

Covariance Matrices. The covariance matrices ${\mathsf{\Sigma }}_{d,d^{\prime }}^{{\gamma }_k}$ and ${\mathsf{\Sigma }}_k^z$ are obtained by Equations (4) and (5), respectively; these matrices only depend on the system model information, so they can be calculated offline before the observed data packets are available.

(2)

LS Linear Prediction and Filtering Recursive Algorithm. At the sampling time k, once the $(k-1)$-th iteration is finished and ${\mathcal{E}}_{k-1}$, ${\mathsf{\Pi }}_{k-1}$, ${\mathsf{\Sigma }}_{k-1}^e$${\mu }_{k-1}$ and $e_{k-1}$ are known, the proposed prediction and filtering algorithm operates as follows:

(2a)

Compute ${\mathcal{Z}}_{k,k-1}={\overline{H}}_k{A}_k{\mathcal{E}}_{k-1}$ and ${\mathcal{Z}}_{k-d,k-1}$, for $d=1,\dots ,(k-1)\wedge D$, by Equation (18). From these matrices, we obtain the observation estimators ${\widehat{z}}_{k-d/k-1}$, for $d=0,1,\dots ,(k-1)\wedge D$, by Equation (19) and (20), and the observation error covariance matrices $P_{k-d,k-d^{\prime }/k-1}^z$, for $d,d^{\prime }=0,1,\dots ,(k-1)\wedge D$, by Equation (19) and (20).

(2b)

Compute ${\mathcal{E}}_k$ by Equation (13) and use $P_{k-d,k-d^{\prime }/k-1}^z$ to obtain the innovation covariance matrix ${\mathsf{\Pi }}_k$ by Equation (9). Then, ${\mathsf{\Sigma }}_k^e$ is obtained by Equation (14) and, from it, the prediction and filtering error covariance matrices, $P_{k/k-s}^x$ and $P_{k/k}^x$, respectively, are obtained by Equation (15).

(2c)

When the new measurement $y_k$ is available, the innovation ${\mu }_k$ is computed by Equation (7) using ${\widehat{z}}_{k-d/k-1}$, for $d=0,1,\dots ,(k-1)\wedge D$, and, from the innovation, $e_k$ is obtained by Equation (12). Then, the predictors, ${\widehat{x}}_{k/k-s}$ and the filter, ${\widehat{x}}_{k/k}$ are computed by Equation (11).

(3)

LS linear fixed-point smoothing recursive algorithm. Once the filter, ${\widehat{x}}_{k/k}$, and the filtering error covariance matrix, $P_{k/k}^x$ are available, the proposed smoothing estimators and the corresponding error covariance matrix are obtained as follows:

For $L=k+1,k+2,\dots ,$ compute the error cross-covariance matrices $P_{k,L-d/L-1}^{xz}$, for $d=0,1,\dots ,(k-1)\wedge D$, using Equation (23) and, from these matrices, ${\mathcal{X}}_{k,L}$ is derived by Equation (22); then, the smoothers ${\widehat{x}}_{k/L}$ and their error covariance matrices $P_{k/L}^x$ are obtained from Equation (21).

5. Computer Simulation Results

In this section, a numerical example is presented with the following purposes: (a) to show that, although the current covariance-based estimation algorithms do not require the evolution model generating the signal process, they are also applicable to the conventional formulation using the state-space model, even in the presence of state-dependent multiplicative noise; (b) to illustrate some kinds of uncertainties which can be covered by the current model with random measurement matrices; and (c) to analyze how the estimation accuracy of the proposed algorithms is influenced by the sensor uncertainties and the random delays in the transmissions.

Signal Evolution Model with State-Dependent Multiplicative Noise. Consider a scalar signal $\left\{x_k;k\ge 0\right\}$ whose evolution is given by the following model with multiplicative and additive noises:

$$x_k=(0.9+0.01{ϵ}_{k-1})x_{k-1}+{\xi }_{k-1},k\ge 1,$$

where $x_0$ is a standard Gaussian variable and $\{{ϵ}_k;k\ge 0\}$, $\{{\xi }_k;k\ge 0\}$ are zero-mean Gaussian white processes with unit variance. Assuming that $x_0$, $\{{ϵ}_k;k\ge 0\}$ and $\{{\xi }_k;k\ge 0\}$ are mutually independent, the signal covariance function is given by $E\left[x_k{x}_s\right]={0.9}^{k-s}{D}_s,s\le k,$ where $D_s=E\left[x_s^2\right]$ is recursively obtained by $D_s=0.8101D_{s-1}+1,$ for $s\ge 1$, with $D_0=1;$ hence, the signal process satisfies assumption (A1) taking, for example, $A_k={0.9}^k$ and $B_s={0.9}^{-s}{D}_s$.

Sensor Measured Outputs. As in [22], let us consider scalar measurements provided by four sensors with different types of uncertainty: continuous and discrete gain degradation in sensors 1 and 2, respectively, missing measurements in sensor 3, and both missing measurements and multiplicative noise in sensor 4. These uncertainties can be described in a unified way by the current model with random measurement matrices; specifically, the measured outputs are described according to Equation (1):

$$z_k^{\left(i\right)}=H_k^{\left(i\right)}{x}_k+v_k^{\left(i\right)},k\ge 1,i=1,2,3,4,$$

with the following characteristics:

• For $i=1,2,3$, $H_k^{\left(i\right)}=C^{\left(i\right)}{\theta }_k^{\left(i\right)}$ and $H_k^{\left(4\right)}=\left(C^{\left(4\right)}+C^{\left(4^{\prime }\right)}{\rho }_k\right){\theta }_k^{\left(4\right)}$, where $C^{\left(1\right)}=C^{\left(3\right)}=0.8$, $C^{\left(2\right)}=C^{\left(4\right)}=0.75$, $C^{\left(4^{\prime }\right)}=0.95$, and $\{{\rho }_k;k\ge 1\}$ is a zero-mean Gaussian white process with unit variance. The sequences $\{{\rho }_k;k\ge 1\}$ and $\{{\theta }_k^{\left(i\right)};k\ge 1\},i=1,2,3,4,$ are mutually independent, and $\{{\theta }_k^{\left(i\right)};k\ge 1\},i=1,2,3,4$, are white processes with the following time-invariant probability distributions:

  −

  ${\theta }_k^{\left(1\right)}$ is uniformly distributed over the interval $[0.1,0.9]$;

  −

  $P[{\theta }_k^{\left(2\right)}=0]=0.3,P[{\theta }_k^{\left(2\right)}=0.5]=0.3P[{\theta }_k^{\left(2\right)}=1]=0.4;$

  −

  For $i=3,4,$ ${\theta }_k^{\left(i\right)}$ are Bernoulli random variables with the same time-invariant probabilities in both sensors $P[{\theta }_k^{\left(i\right)}=1]=p$.

• The additive noises are defined by $v_k^{\left(i\right)}=c_i{\eta }_k^v,i=1,2,3,4$, where $c_1=0.5$, $c_2=c_3=0.75$, $c_4=1$, and $\{{\eta }_k^v;k\ge 1\}$ is a zero-mean Gaussian white process with unit variance. Clearly, the additive noises $\{v_k^{\left(i\right)};k\ge 1\}$, $i=1,2,3,4$, are correlated at any time, with $R_k^{\left(ij\right)}=c_i{c}_j,i,j=1,2,3,4.$

Observations with Bounded Random Delays and Packet Dropouts. Next, according to the theoretical observation model, let us suppose that bounded random measurement delays and packet dropouts, with different delay probabilities, exist in the data transmission. Specifically, assuming that the largest delay is $D=3$, let us consider the observation Equation (2):

$${\displaystyle y_k^{\left(i\right)}=\sum _{d=0}^{(k-1)\wedge 3}{\gamma }_{d,k}^{\left(i\right)}{z}_{k-d}^{\left(i\right)}+w_k^{\left(i\right)},k\ge 1,}$$

where, for $i=1,2,3,4$ and $d=0,1,2,3$, $\{{\gamma }_{d,k}^{\left(i\right)};k>d\}$, are sequences of independent Bernoulli variables with the same time-invariant delay probabilities for the four sensors ${\overline{\gamma }}_{d,k}^{\left(i\right)}={\overline{\gamma }}_d$, where ${\displaystyle \sum _{d=0}^3{\overline{\gamma }}_d\le 1}$. Hence, the packet dropout probability is ${\displaystyle 1-\sum _{d=0}^3{\overline{\gamma }}_d}$. The transmission noise is defined by $w_k^{\left(i\right)}=c_i{\eta }_k^w$, $i=1,2,3,4$, where $\{{\eta }_k^w;k\ge 1\}$ is a zero-mean Gaussian white process with unit variance.

Finally, in order to apply the proposed algorithms, and according to (A5), we will assume that all of the processes involved in the observation equations are mutually independent.

To illustrate the feasibility and effectiveness of the proposed algorithms, they were implemented in MATLAB (R2011b 7.13.0.564, The Mathworks, Natick, MA, USA) and one hundred iterations of the prediction, filtering and fixed-point smoothing algorithms have been performed. In order to analyze the effect of the network-induced uncertainties on the estimation accuracy, different values of the probabilities p of the Bernoulli random variables that model the uncertainties of the third and fourth sensors, and several values of the delay probabilities ${\overline{\gamma }}_d$, $d=0,1,2,3$, have been considered.

Performance of the Prediction, Filtering and Fixed-Point Smoothing Estimators. Considering the values $p=0.5$, ${\overline{\gamma }}_0=0.6$ (packet arrival probability) and ${\displaystyle {\overline{\gamma }}_d=\frac{1-{\overline{\gamma }}_0}{4},d=1,2,3}$ (delay probabilities), Figure 1 displays a simulated trajectory along with the prediction, filtering and smoothing estimations, showing a satisfactory and efficient tracking performance of the proposed estimators. Figure 2 shows the error variances of the predictors ${\widehat{z}}_{k/k-2}$ and ${\widehat{z}}_{k/k-1}$, the filter ${\widehat{z}}_{k/k}$ and the smoothers ${\widehat{z}}_{k/k+1}$ and ${\widehat{z}}_{k/k+2}$. Analogously to what happens for non-delayed observations, the performance of the estimators becomes better as more observations are used; that is, the error variances of the smoothing estimators are less than the filtering ones which, in turn, are less than those of the predictors. Hence, the estimation accuracy of the smoothers is superior to that of the filter and predictors and improves as the number of iterations in the fixed-point smoothing algorithm increases. The performance of the proposed filter has also been evaluated in comparison with the standard Kalman filter; for this purpose, the filtering mean square error (MSE) at each sampling time was calculated by considering one thousand independent simulations and one hundred iterations of each filter. The results of this comparison are displayed in Figure 3, which shows that the proposed filter performs better than the Kalman filter, a fact that was expected since the latter does not take into account either the uncertainties in the measured outputs or the delays and losses during transmission.

Influence of the Missing Measurements. To analyze the sensitivity of the estimation performance on the effect of the missing measurements phenomenon in the third and fourth sensors, the error variances are calculated for different values of the probability p. Specifically, considering again ${\overline{\gamma }}_0=0.6$ and ${\overline{\gamma }}_d=0.1,d=1,2,3,$ Figure 4 displays the prediction, filtering and fixed-point smoothing error variances for the values p = 0.5 to p = 0.9. This Figure shows that, as p increases, the estimation error variances become smaller and, hence, as it was expected, the performance of the estimators improves as the probability of missing measurements $1-p$ decreases.

Influence of the Transmission Random Delays and Packet Dropouts. Considering a fixed value of the probability p, namely, $p=0.5$, in order to show how the estimation accuracy is influenced by the transmission random delays and packet dropouts, the prediction, filtering and smoothing error variances are displayed in Figure 5 when ${\overline{\gamma }}_0$ is varied from 0.1 to 0.9, considering again ${\displaystyle {\overline{\gamma }}_d=\frac{1-{\overline{\gamma }}_0}{4},d=1,2,3}$. Since the behaviour of the error variances is analogous from a certain iteration on, only the results at the iteration $k=50$ are shown in Figure 5. From this figure, we see that the performance of the estimators (predictor, filter and smoother) is indeed influenced by the transmission delay and packet dropout probabilities and, as it was expected, it is confirmed that the error variances become smaller, and hence the performance of the estimators improves, as the packet arrival probability ${\overline{\gamma }}_0$ increases. Moreover, for the filter and the smoothers, this improvement is more significant than for the predictors. In addition, as it was deduced from Figure 2, it is observed that, for all the values of ${\overline{\gamma }}_0$, the performance of the estimators is better as more observations are used, and this improvement is more significant as ${\overline{\gamma }}_0$ increases.

Finally, it must be remarked that analogous conclusions are deduced for other values of the probabilities p and ${\overline{\gamma }}_d,d=0,1,2,3$, and also when such probabilities are different at the different sensors.

6. Conclusions

This paper makes valuable contributions to the optimal fusion estimation problem in networked stochastic systems with random parameter matrices, when multi-step delays or even packet dropouts occur randomly during the data transmission. By an innovation approach, recursive prediction, filtering and fixed-point smoothing algorithms have been designed, which are easily implementable and do not require the signal evolution model, but only the mean and covariance functions of the system processes.

Unlike other estimation algorithms proposed in the literature, where the estimators are restricted to obey a particular structure, in this paper, recursive optimal estimation algorithms are designed without requiring a particular structure on the estimators, but just using the LS optimality criterion. Another advantage is that the current approach does not resort to the augmentation technique and, consequently, the dimension of the designed estimators is the same as that of the original signal, thus reducing the computational burden in the processing centre.