1

Distributed Dimensionality Reduction Fusion
Estimation with Communication Delays in
Cyber-Physical Systems
areas such as intelligent transportation systems [14], multirobot systems [15], and smart grid systems [16]. As one of
important issues in CPSs, real-time state estimation based
on sensor measurements has recently attracted considerable
research interests because state estimate can provide a CPS
with the real-time monitoring and control capability [17],
[18]. For example, estimating the real voltage from sensor
information must be completed before taking certain actions
to regulate the voltage into some desired range in a power
grid [16]. It is noted that the accuracy of state estimation has
an important impact on computing control commands for safe
and efficient operation of a CPS [17]–[19]. Therefore, it is of
theoretical significance and practical relevance to investigate
the problem of information fusion estimation for the CPSs
[20], [21].

C1(t)
C2(t)
C3(t)

Index Terms—Distributed Fusion Estimation, Kalman Filtering, Bandwidth Constraints, Communication Delays, Stability
Analysis, Cyber-Physical Systems.

I. I NTRODUCTION
Information fusion has attracted considerable research interest during the past decades, and has found applications in
a variety of areas, including internet of things [1] and cyberphysical systems (CPSs) [2]. Particularly, multi-sensor fusion
estimation utilizes useful information contained in multiple
sets of data for the purpose of estimating a quantity or parameter in a process [3]. It is widely used in practical applications
because it can potentially improve estimation accuracy and
enhance reliability and robustness against faults [3]–[5].
Many fusion estimation approaches have been presented in
the literature (see [6]–[12], and the references therein). At
the same time, advances in embedded computing, communication, and related hardware technologies have recently
brought the paradigm of CPSs to a new research frontier [13].
Moreover, CPSs have found applications in a broad range of
B. Chen was with the Department of Mathematics, City University of
Hong Kong, Hong Kong, 999077. He is now with the School of Electrical
and Electronic Engineering, Nanyang Technological University, 639798
Singapore (email:bchen@aliyun.com).
D. W. C. Ho is with the Department of Mathematics, City University of
Hong Kong, Hong Kong, 999077. (email:madaniel@cityu.edu.hk).
G. Hu is with the School of Electrical and Electronic Engineering, Nanyang Technological University, 639798 Singapore (email:
gqhu@ntu.edu.sg).
L. Yu is with the College of Information Engineering, Zhejiang University
of Technology, HangZhou 310023, China (email:lyu@zjut.edu.cn).

Cn(t)

ensor!no"e

pGmG
j

Abstract—This paper studies the distributed dimensionality
reduction fusion estimation problem with communication delays
for a class of cyber-physical systems (CPSs). The raw measurements are preprocessed in each sink node to obtain the local
optimal estimate (LOE) of a CPS, and the compressed LOE
under dimensionality reduction encounters with communication
delays during the transmission. Under this case, a mathematical
model with compensation strategy is proposed to characterize
the dimensionality reduction and communication delays. This
model also has the property to reduce the information loss
caused by the dimensionality reduction and delays. Based on
this model, a recursive distributed Kalman fusion estimator
(DKFE) is derived by optimal weighted fusion criterion in the
linear minimum variance sense. A stability condition for the
DKFE, which can be easily verified by the exiting software,
is derived. In addition, this condition can guarantee that
estimation error covariance matrix of the DKFE converges
to the unique steady-state matrix for any initial values, and
thus the steady-state DKFE (SDKFE) is given. Notice that
the computational complexity of the SDKFE is much lower
than that of the DKFE. Moreover, a probability selection
criterion for determining the dimensionality reduction strategy
is also presented to guarantee the stability of the DKFE. Two
illustrative examples are given to show the advantage and
effectiveness of the proposed methods.

Communication!!Networks

arXiv:1802.03122v1 [] 9 Feb 2018

Bo Chen, Member, IEEE, Daniel W. C. Ho, Fellow, IEEE, Guoqiang Hu, Member, IEEE, Li Yu, Member, IEEE

ink!no"e Ci(t) #artia$!com%onents!o&!t'e!%'(sica$!%rocess!)*t+

Fig. 1. Information fusion estimation for a class of spatially distributed
physical systems over communication networks: i) x(t) is the state of the
physical process, where x(t)=col{c1(t),c2(t),c3(t),...,cn(t)}; ii) sensor node
only measures the target information; iii) sink node is a gateway, which
is responsible for receiving measurements, computing the local optimal
estimate (LOE) and sending the LOE to an information fusion center via
communication networks.

There mainly exist two kinds of fusion architectures:
centralized fusion structure and distributed fusion structure.
However, the distributed fusion structure is generally more
robust and fault-tolerant as compared with the centralized
fusion structure [4]–[8]. This motivates us to consider the distributed fusion estimation problem in this paper for a class of
CPS architecture (see Fig.1), where system state is spatially
distributed in the physical space. When the local estimates
are transmitted to the fusion center (FC) via communication
channels, bandwidth constrains and communication delays
are unavoidable in communication networks [22]. Moreover,
the above two factors can degrade the fusion estimation
performance because of the information loss caused by bandwidth and delay constrains [23]–[25]. Thus, how to design
distributed fusion methods in the presence of bandwidth and
delay constraints is essential for real-time state estimate of
CPSs.

2

A. Related Work
Regarding the problem of bandwidth constraints in multisensor systems, as pointed out in [26], there are mainly
two approaches to reduce the communication traffic: the
quantization method (see [27]–[29], and references therein)
and the dimensionality reduction method (see [30]–[32],
and the references therein). Particularly, the dimensionality
reduction method to the original multi-sensor observations
was designed in [33] based on the principal components
analysis, while the dimensionality reduction strategy with the
quantization error was developed in [34] to deal with stable
multi-sensor fusion systems. Under the distributed fusion
structure, when the physical state x(t) as shown in Fig.1 is
multi-dimensional (or even high-dimensional) in a CPS, it is
unrealistic to completely send the local estimate of the state
x(t) to the FC via a bandwidth-constrained communication
channel. In this sense, bandwidth constraint in the CPSs is
the primary consideration when designing a distributed state
fusion estimator. Notice that, to reduce the communication
traffic, the idea of the dimensionality reduction method is that
a multi-dimensional signal is directly converted into a lowdimensional signal, while the idea of the quantization method
is that the number of coding bits for each component of a
multidimensional signal is reduced before being transmitted.
Meanwhile, the quantization usually results in nonlinear
dynamics, and it is difficult to find a data compression
operator analytically, particularly, for the multidimensional
signals. Therefore, the dimensionality reduction method can
provide an attractive alternative to solve the distributed fusion
estimation problem with bandwidth constraints in the CPSs.
Though the dimensionality reduction fusion estimation
algorithms have been proposed in [30]–[34] to reduce the
communication traffic, the communication delays, which
occurs during the transmission, were not taken into account.
With the communication delays, the dimensionality reduction fusion estimation must solve two challenging issues:
one is how to compensate the information loss caused by
the communication delays and bandwidth constraints under
a unified mathematical model; The other one is how to
fuse the asynchronous local compressed estimates because
of communication delays. Notice that the centralized and
distributed fusion estimation algorithms have been proposed
in [23], [35]–[40] based on different communication delay
models, however, the main results in [23], [35]–[40] cannot
be extended to the case of the dimensionality reduction
estimation with communication delays. The reason is that the
data compression and information compensation in dimensionality reduction may change the property of the original
measurements (e.g., the statistical correlation in [30], [32]
has been changed under the Kalman fusion structure). Under
this case, we have studied the information fusion estimation
problem in [20], [24] for the CPSs with bandwidth constraints and communication delays. It should be pointed out
that the steady-state fusion estimator with simple calculation
cannot be obtained based on the proposed communication
model in [20], while the covariance intersection (CI) fusion
strategy in [24] was suboptimal because fusion estimator was
determined by minimizing an upper bound of estimation error
covariance.
B. Contributions
Motivated by the aforementioned analysis, we study the
distributed stochastic dimensionality reduction fusion esti-

mation problem with communication delays for the CPSs.
Notice that the information loss is inevitable because of
the dimensionality reduction and communication delays, and
such a fusion estimation with incomplete information will
degrade the estimation performance. Since the delays are
caused by communication channels, the key issue is how
to design an efficient dimensionality reduction strategy to
guarantee the stability of the distributed fusion estimator. Although our previous works in [20], [24], [30] have studied the
related stochastic dimensionality fusion estimation problems,
there are still fundamental problems that cannot be solved up
to now. In detail,
• When only considering stochastic dimensionality reduction strategy, the stable probability selection criterion
in [30] was derived from the inequality relaxation of
the matrix trace. However, the inequality relaxation will
lead to certain conservatism, thus how to find a new
derivation idea to reduce the conservatism is very important for the application of the proposed dimensionality
reduction strategy. Notice that the stability conditions
in [20] were directly derived from the similar derivation
in [30], and thus the corresponding conservatism cannot
also be avoided in [20].
• When considering stochastic dimensionality reduction
strategy under communication delays, the distributed
CI fusion estimator in [24] was suboptimal because
the corresponding optimization objective was an upper
bound of the estimation error covariance matrix. Particularly, the CI fusion results in [24] needed to solve nonconvex nonlinear optimization problems online at each
time, which may lead to a large number of calculation.
Though the distributed fusion estimator in [20] was optimal based on the optimal weighted fusion criterion, the
model of communication delays cannot be applicable to
the case of time-varying delays. More importantly, the
computational complexity of the fusion estimator in [20]
was also high. Obviously, the common disadvantage of
the results in [20] and [24] is the high computation cost,
and the optimal weighed fusion criterion can provide
the optimal and analytic solutions. Therefore, based on
the optimal weighted fusion criterion, how to design
steady-state dimensionality reduction fusion estimators
with simple calculation is of great significance in the
presence of communication delays.
We shall solve the above two problems, and the main
contributions of this paper can be summarized as follows:
• An optimal distributed Kalman fusion estimator (DKFE)
is derived in the linear minimum variance sense when
there are bandwidth and communication delay constraints in CPSs, and each weighting fusion matrix is
calculated by the analytic form.
• A delay-dependent and probability-dependent stability
condition is derived such that the fusion estimation
error covariance matrix of the DKFE converges to a
unique steady-state matrix for any initial values. Under
this condition, the steady-state DKFE, which has much
lower computational complexity as compared with the
DKFE, is given. Moreover, when each communication
delay is known, the probability selection criterion for determining dimensionality reduction strategy is presented
to guarantee the stability of the DKFE.
• Compared with the fusion estimation method in [20], the

3

•

model of communication delays in this paper does not
require that each sink node knows the communication
delay in advance, and the steady-state DKFE with
simple calculation is derived (see Remark 1). Since the
covariance intersection fusion criterion in [24] is suboptimal, the estimation performance of the designed DKFE
must be better than that of the fusion estimator in [24]
when each communication delay is constant. Moreover,
the computation cost of the steady-state DKFE must be
lower than that of the CI fusion estimator in [24] (see
Remark 2).
When there is no communication delay for the scenario
described in Fig.1, it is shown that the stability condition
in this paper has less conservatism than the result in
[30]. This is because a new derivation idea without any
inequality relaxation is proposed to design the stochastic dimensionality reduction strategy. Moreover, when
considering communication delays, the corresponding
stability analysis is also based on this new derivation
idea. Notice that it is difficult to obtain the stability
condition by using the derivation idea in [30] when
the communication delay is modeled in this paper (see
Remarks 7-8).

The rest of this paper is organized as follows. Section II
presents the problem formulation. The finite-horizon DKFE
is designed in Section III. In Section IV, the stability condition and the steady-state DKFE are derived, and the probability selection criteria are given to determine satisfactory compression operators. Two illustrative examples are presented
in Section V to show the advantage and effectiveness of the
proposed approaches, and then the conclusions are drawn in
Section VI.
Notations: The notations used throughout the paper are
fairly standard. The superscript ′ T′ represents the transpose,
and E{·} is the mathematical expectation. Im represents the
identity matrix of size m × m, while diag{·} stands for
a block diagonal matrix. Prob{A} means the occurrence
probability of the event A, while Tr(B) denotes the trace
of the matrix B. ||A||2 represent the 2-norm of the matrix
A. x⊥y denotes that x and y are orthogonal vectors, and
col{a1 , · · · , aL } represents the column vector that is composed of the elements a1 , · · · , aL . The symbol lcm(a, b) is
the least common multiple of a and b, while rank(A) denotes
the rank of the matrix A. The function f◦~ (t) is defined
∆
by f◦~ (t) = f (f (· · · (f (t)) · · · )), and X > (<)0 denotes
| {z }
~ times

from each sensor are sent to sink nodes, the ith sink node’s
measurement yi (t) ∈ Rqi is modeled by:
yi (t) = Ci x(t) + vi (t)(i = 1, 2, · · · , L),

(2)

where Ci is the measurement matrix with appropriate dimension, and vi (t) is the measurement noise. Moreover, w(t)
and vi (t) are uncorrelated zero-mean Gaussian white noises
satisfying
E{[wT (t) viT (t)]T [wT (t1 ) vjT (t1 )]}
, (3)
= δt,t1 diag{Qw , δi,j Qvi }
where δt,t1 is defined by:

1
δt,t1 =
0

if t = t1
.
if t =
6 t1

(4)

Then, based on the measurements {yi (1), · · · , yi (t)}, the
local optimal estimate (LOE) x̂i (t) is given by the Kalman
filter:
x̂i (t) = GKi (t)Ax̂i (t − 1) + Ki (t)yi (t),
where

(5)

∆

GKi (t) = In − Ki (t)Ci .

(6)

∆

Define x̃i (t) = x(t) − x̂i (t). Then, the optimal gain matrix Ki (t) and the local estimation error covariance ma∆
trix Pii (t) = E{x̃i (t)x̃T
i (t)} are calculated by

∗
 Ki (t) = Pii (t)CiT [Ci Pii∗ (t)CiT + Qvi ]−1
Pii (t) = GKi (t)Pii∗ (t)
.
(7)
 ∗
Pii (t) = APii (t − 1)AT + Qw
Moreover, it follows from (1), (5) and (7) that the lo∆
cal estimation error cross-covariance matrix Pij (t) =
T
E{x̃i (t)x̃j (t)}(i 6= j) is calculated by:
Pij (t) = GKi (t)[Qw + APij (t − 1)AT ]GT
Kj (t).

(8)

Physical Process: x(t)
y1 ( t )

yi (t )

Sink node 1

...

LOE x1 (t )

Sink node i

y L (t )

...

LOE xi (t )

 s1 (t )
Selected ASC x

Sink node L
LOE xL (t )

 si (t )
Selected ASC x

 sL (t )
Selected ASC x

Communication Channels with Communication Delays

 s1 (t d1)
x

...

 si (t di )
x

...

 sL (t dL)
x

a positive-definite (negative-definite) matrix.

II. P ROBLEM F ORMULATION

Fusion Center

Fig. 2. Distributed dimensionality reduction fusion estimation with communication delays in CPSs

A. Dimensionality Reduction and Communication Delays
Consider the physical process in Fig.1 described by the
following discrete state-space model:
x(t + 1) = Ax(t) + w(t),

(1)

where x(t) ∈ Rn (n > 1) is the state of the process,
w(t) is the system noise, and A is a constant matrix with
appropriate dimension. As pointed out in [18], the model
(1) is widely adopted for describing state dynamics of CPSs
including power systems, smart grid infrastructures, and
building automation systems, etc. When the measurements

Under the distributed fusion structure, each LOE x̂i (t)
must be sent to the FC to design an optimal fusion estimator.
However, it is unrealistic to send the complete information
included in x̂i (t)(∈ Rn ) to the FC over communication
networks because almost all communication network can
only carry a finite amount of information per unit time. This
problem is especially prominent in the fusion estimation for
the large-scale CPSs integrated by wireless sensor networks.
To reduce communication traffic, only ri (1 ≤ ri < n) components of the ith LOE x̂i (t) are allowed to be transmitted
to the FC at each time, and other components are discarded.

4

Compared with the original LOE x̂i (t), the dimension of the
transmitted signal is reduced. In this sense, the above method
can be viewed as one of the dimensionality reduction strategies. According to this dimensionality reduction strategy, the
allowed sending components (ASC) of
,x̂i (t) has ∆i possible
riQ
−1
ri
Q
ℓri . Then, at a
(n − ℓni )
cases, where ∆i =
ℓri =1

ℓn
i =0

particular time, only one vector signal, which is taken from
one group of the above ∆i cases, is selected and transmitted
to the FC, and this selected signal is denoted by x̂si (t) ∈ Rri .
When x̂si (t) is sent to the FC by the sink node, the FC
will receive the data packet containing x̂si (t) at time t + di
because of communication delay. Let x̄si (t) denote the local
estimation information received by the FC at time t. Then,
x̄si (t) in the FC is given by:
x̄si (t) = x̂si (t − di ).

(9)

Up to now, the problem of dimensionality reduction and
communication delays has been presented, and the process
diagram is shown in Fig.2.
It is noted that the signal x̂si (t) only takes one element
from the following finite set:
Si (t) = {x̂~sii (t)|~i = 1, 2, · · · , ∆i },

(10)

where x̂~sii (t) ∈ Rri represents one group of ASCs. To
characterize the determining process of x̂si (t), we introduce
the following indicator functions:

1 if x̂si (t) = x̂~sii (t)
i
σ~i (t) =
,
(11)
0 if x̂si (t) 6= x̂~sii (t)
where σ~i i (t)(~i = 1, 2, · · · , ∆i ) are required to satisfy
X∆ i
σ~i i (t)σ~i 0 (t) = 0(~i 6= ~0i ) ,
σ~i i (t) = 1 (12)
~i =1

i

such that x̂si (t) only takes one ASC from the set (10) at
time t, i.e.,
X∆ i
x̂si (t) =
(13)
σ~i i (t)x̂~sii (t).
~i =1

Then, it is derived from (9) and (13) that
X∆ i
x̄si (t) =
σ~i i (t − di )x̂~sii (t − di ).
~i =1

(14)

At time t, if the fusion estimate of x(t) is directly designed
based on x̄si (t), the fusion estimation performance must
be poor because of the communication delays and the untransmitted component of x̂i (t). In this case, the compensating state estimate (CSE) of x(t), denoted by x̂ci (t), can be
modeled as follows:
x̂ci (t) = Adi Hi (t − di )x̂i (t − di )
, (15)
+Adi [In − Hi (t − di )]Ax̂ci (t − di − 1)
where Hi (t − di ) is determined by
X∆ i
Hi (t) =
σ~i i (t)H~i i = diag{γ1i (t), · · · , γni (t)}.(16)
~i =1

H~i i

Here,
represents a diagonal matrix that contains ri diagonal elements “1” and n − ri diagonal elements “0”. Then
it follows from (11) and (12) that
Xn
γℓi (t) = ri (i = 1, · · · , L), (17)
γℓi (t) ∈ {0, 1},
ℓ=1

where γℓi (t) = 1 means that the ℓth component of x̂i (t) is
selected and sent to the FC, while γℓi (t) = 0 means that the
ℓth component of x̂i (t) is discarded. Particularly, at time t,

the compensation strategy in the CSE model (15) is reflected
by the following aspects:
• The un-transmitted components of x̂i (t) are compensated by the one-step prediction based on x̂ci (t− di − 1).
• The delayed information x̄si (t) is compensated by the
di -step prediction based on x̂ri (t−di ), where x̂ri (t−di ) =
Hi (t − di )x̂i (t − di ) + [I − Hi (t − di )]Ax̂ci (t − di − 1).
Remark 1. In [20], at the ith sink node, the di -step
prediction
based on the local estimate x̂i (t) was given by

x̂di i (t) = Adi x̂i (t) . Due to the bandwidth constraints, only

ri components of x̂di i (t) were allowed to be sent. Then, the
CSE of x(t), denoted as x̂cdi (t), was given by (i.e., the model
(18) in [20]):
x̂cdi (t) = Hic (t − di )Adi x̂i (t − di )
+[I − Hic (t − di )]Ax̂f (t − 1),

(18)

x̄di i (t) = x̂si i (t − di (t)),

(19)

x̄di i (t) = x̂si i (t − dui )

(20)

where the definition of Hic (t − di ) is the same as that of
Hi (t−di ), and x̂f (t−1) denotes the fusion estimate designed
by [20]. For the CSE model (18), the di -step prediction
x̂di i (t) must be completed at the sink node, which implies
that each sink node must know the communication delay
from the sink node to the FC in advance. Under this case,
when the communication delay is unknown for the sink node
or time-varying, the model (18) will be invalid. Different
from the modeling method in [20], the CSE model (15) does
not require that each sink node knows the communication
delay in advance, and thus the model (15) can be more
easily implemented in a practical system. Particularly, when
considering the time-varying communication delay di (t), the
local estimation information received by the FC, denoted as
x̄di i (t), is given by:
where x̂si i (t) denotes the selected ASC at the sink node.
Meanwhile, it is reasonable to consider that the time-varying
delay di (t) is bounded in practical applications, and satisfies
di (t) ≤ dui . Then, by resorting to the buffers at the FC, each
time-varying delay can be prolonged to its upper bound dui
at each time, i.e., the model (19) is reduced to:

Since the structure of (20) is the same as that of (9), the case
of time-varying delays can still be modeled by (15). Notice
that the CSE model (18) in [20] will not be applicable to
this case, because the time-varying communication delays
are only known to the FC, and each sink node impossibly
know the time-varying delays a priori. On the other hand,
the stability condition in [20] could only guarantee the MSE
of the fusion estimator converged to a steady-state value.
It should be pointed out that the computational complexity
of the fusion estimator in [20] is a slightly high, yet the
corresponding steady-state fusion estimator cannot be derived
from the stability condition in [20]. In contrast, the steadystate DKFE with simple calculation can be designed based
on the stability condition in Theorem 3.
Remark 2. For the case of time-varying delays, the
estimation error cross-covariance matrices cannot be obtained
under the dimensionality reduction strategy in this paper.
Fortunately, the covariance intersection (CI) fusion criterion
does not need the cross-covariance matrices. Therefore, the
distributed CI fusion estimation algorithm was developed in
[24] to deal with the time-varying delays. Notice that the

5

CI fusion criterion is not optimal because the optimization
objective is an upper bound of estimation error covariance
matrix, and each weighting matrix is obtained by solving
non-convex nonlinear optimization problems at each time.
Different from the fusion criterion in [24], the optimal
weighted fusion criterion with analytic solutions is used to
design the DKFE in this paper. Thus, when considering the
constant communication delays, the estimation performance
of the DKFE is better than that of the fusion estimator in
[24]. On the other hand, as pointed out in Remark 1, the
designed fusion estimation algorithms in this paper can be
also applicable to the case of time-varying communication
delays. However, it is difficult to show whose estimation
performance is optimal between the DKFE in this paper
and the fusion estimator in [24] when dealing with timevarying delays. This is because the conservatism in this
paper is introduced from the delay model (i.e., prolonging the
time-varying delay to its upper bound at each time), while
the conservatism in [24] is introduced from the CI fusion
criterion (i.e., minimizing an upper bound of the fusion
estimation error covariance). However, from the perspective
of computational complexity, the steady-state DKFE in this
paper is better than the CI fusion estimator in [24] whenever
considering the constant delays or time-varying delays.
B. Problem of Interest
It is concluded from (13) that the selected ASC x̂si (t)(∈
Rri ) at the sink node is determined by the binary variables σ~i i (t)(~i = 1, 2, · · · , ∆i ). On the other hand, it is
known from (15) that the design of optimal σ~i i (t)(~i =
1, 2, · · · , ∆i ) must be completed at the FC, because the
communication delay (from the sink node to the FC) and
each CSE x̂ci (t) are only obtained by the FC, but these
information are unknown to each sink node. Therefore, an
optimal x̂si (t) may be difficult to be designed at the sink
node. Based on the above consideration, let each binary
variable σ~i i (t) be generated in a random way at the sink
i
node, and let random variables {σ1i (t), σ2i (t), · · · , σ∆
(t)}
i
obey the categorical distribution satisfying

∆

i
}. This means that when
where ςi = col{π1i , π2i , · · · , π∆
i
i
each selection probability π~i is given by (21–22), γℓi in
(23) will be determined by (24). Notice that the selection
probabilities ςi (i = 1, 2, · · · , L) are to be designed in this
paper for guaranteeing the stability of the DKFE.
∆
Let x̃ci (t) = x(t) − x̂ci (t) denote the estimation error of
each CSE. Then, it follows from (1) and (15) that

x̃ci (t) = Adi Hi (t − di )x̃i (t − di )
+Adi [In − Hi (t − di )]Ax̃ci (t − di − 1)
(25)
+Adi [In − Hi (t − di )]w(t − di − 1) + Fw (di , t),
where Fw (di , t) is determined by the following function:
Xg
∆
Aθ−1 w(t − θ).
(26)
Fw (g, t) =
θ=1

When E{x̂ci (−dni )} = E{x(−dni )}(dni = 0, 1, · · · , di ), it
is concluded from (3), (25) and the fact E{x(t)} = E{x̂i (t)}
that each CSE x̂ci (t) is unbiased, i.e.,
E{x̂ci (t)} = E{x(t)}(i = 1, 2, · · · , L).

According to the CSEs x̂ci (t)(i = 1, 2, · · · , L) in the FC, the
DKFE for the addressed CPSs is given by:
x̂(t) =


i
if i = j, t = t1

 δ~i ,~i0i E{σ~i (t)}
.(21)
i
E{σ
(t)}E{σ
(t
)}
if
i = j, t 6= t1
~i
=
~0i 1

 E{σ i (t)}E{σ j 0 (t1 )} if i =
6 j
~i
~
j

Under this case, a group of ASC x̂~sii (t)(~i ∈ {1, 2, · · · , ∆i })
in the set (10) is randomly selected as the x̂si (t) at time t.
Moreover, the occurrence probabilities of the cases σ~i i (t) =
1 and σ~i i (t) = 0 are given by Prob{σ~i i (t) = 1} = π~i i and
Prob{σ~i i (t) = 0} = 1 − π~i i , where the selection probability
π~i i ≥ 0 satisfies:
X∆ i
(22)
π~i i = 1 (i ∈ {1, 2, · · · , L}).
~i =1

Then, it is concluded from (16) and (21) that the binary variables γℓi (t)(ℓ = 1, · · · , n) in (17) are independent Bernoulli
∆
distributed white noise sequences with Prob{γℓi (t) = 1} =
∆
γℓi and Prob{γℓi (t) = 0} = 1 − γℓi , which yields
∆

Hi = E{Hi (t − di )} = diag{γ1i , γ2i , · · · , γni }.

(23)

γℓi = Uℓi ςi (ℓ = 1, · · · , n; i = 1, · · · , L),

(24)

From (16) and (23), there must exist a constant matrix Uℓi ∈
R1×∆i such that

XL

i=1

Ωi (t)x̂ci (t),

(28)

PL
where i=1 Ωi (t) = In , and combining (27) yields that the
DKFE x̂(t) is unbiased if E{x̂ci (dni )} = E{x(dni )}(dni =
0, 1, · · · , di ).
Consequently, the problems to be solved in this paper are
described as follows:
1) When the selection probabilities π~i i (~i
=
1, · · · , ∆i ; i = 1, · · · , L) satisfying (22) are given in
advance, the aim is to design optimal weighting matrices
Ω1 (t), · · · , ΩL (t) such that the MSE of the DKFE x̂(t) is
minimal at each time step, i.e.,
{Ω1 (t), · · · , ΩL (t)}
= arg P min
L
i=1

E{σ~i i (t)σ~j 0 (t)}
j

(27)

Ωi (t)=I

E{[x(t) − x̂(t)]T [x(t) − x̂(t)]}.(29)

2) Find stability conditions, which are dependent on the
communication delay di in (9) and the selection probability
π~i i in (22), such that the estimation error covariance matrix
of the DKFE converges to a unique positive matrix, i.e.,
lim E{[x(t) − x̂(t)][x(t) − x̂(t)]T } = P,

t→∞

(30)

and P is independent of the initial values.
Remark 3. When the ith sink node knows the selection
probability ςi in advance, the binary variables σ~i i (t)(~i =
1, · · · , ∆i ) obeying the categorical distribution will be randomly generated at each time step, and then the selected ASC
x̂si (t) can be determined by (13) at the sink node. Under
this case, one of the important issues in this paper is how
to design the satisfactory probability selection criteria, which
will be solved in Section IV. On the other hand, when the
result (30) holds, the limit of each weighting matrix Ωi (t)
must exist, and will be independent of the initial values.
This is because the estimation error covariance matrix of the
DKFE is dependent on each time-varying matrix Ωi (t). In
such a case, the steady-state DKFE with simple calculation
will be given in this paper.

6

III. F INITE -H ORIZON DKFE

FOR THE

CPS S

In this section, the recursive DKFE will be derived by
using the optimal fusion criterion weighted by matrices in the
∆
linear minimum variance sense. Define x̃(t) = x(t)−x̂(t) and
nL×n
Ia = col{In , · · · , In } ∈ R
. Then, from the results in
[7], [8], the optimal weighting matrices Ω1 (t), · · · , ΩL (t) in
(29) and the corresponding fusion estimation error covariance
∆
matrix P (t) = E{x̃(t)x̃T (t)} can be calculated by:

ΦF
xi (t1 , g, t2 ) are given by:
Q

t1 −t2 −2
Φw
ΦKi (t1 − ϕi )
xi (t1 , t2 ) = Co (t1 , t2 )
ϕi =0
×GKi (t2 + 1)Qw
 Q

t1 −t2 −1

ΦKi (t1 − ϕi ) Pij (t2 )

ϕi =0
Φxxji (t1 , t2 ) =
if t1 ≥ t2

 [Φxi (t , t )]T if t < t
1
2
xj 2 1
Xg
θ−1 T
w
F
Φxi (t1 , t2 − θ)[A ]
Φxi (t1 , g, t2 ) =
θ=1

[Ω1 (t), Ω2 (t), · · · , ΩL (t)] = (IaT Ξ−1 (t)Ia )−1 IaT Ξ−1 (t)(31)
P (t) = (IaT Ξ−1 (t)Ia )−1 ,

(32)

where the weighting matrices Ωi (t)(iP= 1, 2, · · · , L) deterL
mined by (31) satisfy the constraint i=1 Ωi (t) = In , and

Ξ(t) = (Ξij (t))nL×nL , Ξij (t) = E{x̃ci (t)(x̃cj (t))T }. (33)

It is concluded from (31) and (33) that if the computation
procedure of Ξ(t) is given, then the optimal weighting
matrices Ωi (t)(i = 1, · · · , L) in (31) can be thus obtained.
In what follows, six lemmas will be given before deriving
the recursive form of Ξij (t). For notational convenience, the
following indicator function is introduced:
Co (t1 , t2 ) =



1 if t1 > t2
.
0 if t1 ≤ t2

(34)

Meanwhile,
if τ1 > Pτ2 , it will be specified that
Qτ2
τ2
τ =τ1 F (τ ) = Im and
τ =τ1 G(τ ) = 0, where F (τ ) ∈
Rm×m and G(τ ) ∈ Rn×n represent different matrix functions with respect to the variable τ .
Lemma 1 [30] For stochastic matrices U , B, G, where

Φw
F (g, t1 , t2 )

t2 −1

= Co (t1 , t2 )Co (t2 , t1 − g − 1)A

(36)

(37)
(38)

Qw ,(39)

where δi,j and Co (t1 , t2 ) are determined by (4) and (34),
respectively. Pij (t2 ) in (37) is calculated by (7) or (8).
Proof: See A.1 in Appendix.
Lemma 3 Define

∆
∆
0

fi (t) = t − di − 1, fio
(t) = t



∆
χ (t ,t )
 χ (t , t ) =
min{χi (t1 , t2 )|fioi 1 2 (t2 ) − t1 ≤ 0}
i 1 2
.(40)
∆
c
T

Θw

xci (t1 , t2 ) = E{x̃i (t1 )w (t2 )}


∆
 F
Θxci (t1 , g, t2 ) = E{x̃ci (t1 )FT
w (g, t2 )}
v

j
F
Then, Θw
xci (t1 , t2 ), Θxc (t1 , t2 ), Θxci (t1 , g, t2 ) are given by:
i

Θw
t2 )
xci (t1 , n

= Co (t1 − di , t2 )
Pχi (t1 ,t2 )−1
{δf ~+1 (t1 ),t2 H~Adi H̄Adi Qw
×
~=0
io
(41)
~
w
, t2 )
+Hdi Φw
xi (fio (t1 ) − di , t2 )} + ΦF (di , t1o
Pχi (t1 ,t2 )−1
~
+ ~=1
HAdi Φw
F (di , fio (t1 ), t2 )
Xg
θ−1 T
Θw
) , (42)
ΘF
xci (t1 , t2 − θ)(A
xci (t1 , g, t2 ) =
θ=1

where



Hdi = Adi Hi , H̄Adi = Adi − Hdi
,
HAdi = Adi [In − Hi ]A

(43)

where “⊗” is defined as [G1 ⊗ G2 ]ij = G1ij G2ij , and the
product “⊙” for the matrices U and B is defined by

where Hi is determined by (23), while Φw
xi (t1 , t2 ) and
and
(39).
Φw
(d
,
t
,
t
)
are
calculated
by
(36)
i
1
2
F
Proof: See A.2 in Appendix.
The statistical correlation between x̃i (t) and w(t) is presented in Lemma 2, while Lemma 3 gives the statistical correlation between x̃ci (t) and w(t). Additionally, to obtain the
Ξij (t), it is still required to know the statistical correlations
between x̃i (t1 ) and x̃ci (t2 ), and between x̃ci (t1 ) and x̃ci (t2 ),
which will be derived in Lemmas 4–6.
Lemma 4 Define
(
∆
Γij (t) = E{x̃i (t)[x̃cj (t)]T }
.
(44)
∆
Γij (t1 , t2 ) = E{x̃i (t1 )[x̃cj (t2 )]T }(t1 > t2 )


u 1 b1 · · · u 1 bn


U ⊙ B =  ... . . . ...  .
u n b1 · · · u n bn

Then, Γij (t) is calculated by the following recursive form:
Q

dj
Γij (t) =
(t
−
ϕ
)
Γij (t − dj − 1)HT
Φ
i
K
i
ϕj =0
Adj
xj
F
,(45)
+Φxi (t, t − dj )HT
(t,
d
,
t)
+
Φ
j
x
dj
i
T
w
+Φxi (t, t − dj − 1)H̄Adj

∆

∆

U = diag{u1 , · · · ,u
 1 , · · · , bn }
 n }, B = diag{b
g11 · · · g1n
.
∆ 

G =  ... . . . ... 
gn1 · · · gnn

If each random variable gij in G is independent of any
random variables of uk and bk (k = 1, 2, · · · , n), then
E{U GB} = E{U ⊙ B} ⊗ E{G},



Lemma 2 Define

∆
∆
T

ΦKi (t) = GKi (t)A, Φw

xi (t1 , t2 ) = E{x̃i (t1 )w (t2 )}


∆
 xj
Φxi (t1 , t2 ) = E{x̃i (t1 )x̃T
j (t2 )}
,(35)
∆
F

Φxi (t1 , g, t2 ) = E{x̃i (t1 )FT

w (g, t2 )}


 w
∆
ΦF (g, t1 , t2 ) = E{Fw (g, t1 )wT (t2 )}

where GKi (t) and Fw (g, t2 ) are determined by (6) and (26),
vj
xj
respectively. Then, Φw
xi (t1 , t2 ), Φxi (t1 , t2 ), Φxi (t1 , t2 ) and

x

where Hdj , H̄Adj , HAdj are given by (43), while Φxji (t, t −
w
dj ), ΦF
xi (t, dj , t), Φxi (t, t − dj − 1) are computed by (36),
(37) and (38). In this case, Γij (t1 , t2 ) is calculated by:
Y

t1 −t2 −1
Γij (t1 , t2 ) =
ΦKi (t1 − ϕi ) Γij (t2 ). (46)
ϕi =0

Proof: See A.3 in Appendix.
Lemma 5 Define
∆

Ψij (t) = E{x̃i (t − di )[x̃cj (t − dj − 1)]T }.

(47)

7

For i = j, Ψii (t) is calculated by
Ψii (t) = ΦKi (t − di )Γij (t − di − 1).

(48)

∆

For i 6= j, let ηij = min{ηij |ηij (dj + 1) − di ≥ 0}. Then,
Ψij (t) is calculated by
Pηij −1
x
κ
(t) − dj )HT
Ψij (t) = κ=1
{[Φxji (t − di , fjo
dj
κ+1
T
w
+Φxi (t − di , fjo (t))H̄Adj
T
κ
(t − di , dj , fjo
(t))](Hκ−1
+ΦF
Adj ) }
i
xQ
 , (49)
ηij (dj +1)−1−di
+
ΦKi (t − di − ϕi )
ϕi =0
η

η −1

ij
)T
×Γij (fjoij (t))(HAd
j

where Hdj , H̄Adj , HAdj are given by (43); Γij (t) is computed
xj
κ+1
κ
by (45), while Φw
xi (t − di , fjo (t)), Φxi (t − di , fjo (t) − dj )
κ
F
and Φxi (t − di , dj , fjo (t)) are calculated by (36), (37), (38).
Proof: See A.4 in Appendix.
Remark 4. It should be pointed out that Ψij (t) in Lemma
5 is not a special case of Γij (t1 , t2 ) in Lemma 4 because
there may exist the case di > dj + 1. Notice that when di ≤
d
according to the definitions of
j + 1, one has ηij = 1.
P
QThen,
τ2
τ2
G(τ
)
=
0
and
F
τ =τ1
=τ1 (τ ) = I for τ1 > τ2 ,Ψij (t) is
τ Q
dj −di
calculated by Ψij (t) =
ϕi =0 ΦKi (t − di − ϕi ) Γij (t −
dj − 1).
Lemma 6 Define

∆
∆

τij = τji = lcm(di + 1, dj + 1), τdi = τij /(di + 1)



∆
 Υ (t) =
E{x̃ci (t − di − 1)[x̃cj (t − dj − 1)]T }
ij
,(50)
∆
w
T

Υxij (t) = E{x̄w

fi (t)[x̄fj (t)] }


 c
∆
c
T
Υij (t) = E{x̄w
fi (t)[x̃j (t − τij )] }

where x̄w
f (t) is defined as follows:

τ −1
∆
1

(t) − di ), · · · , x̃i (fiodi (t) − di )}
x̃fi (t) = col{x̃i (fio



τ
∆

2
(t)), · · · , w(fiodi (t))}
wfi (t) = col{w(fio
(51)
.
τ −1
∆
1

(t)), · · · , Fw (di , fiodi (t))}
Ffi (t) = col{Fw (di , fio



∆
 w
x̄f (t) = col{x̃fi (t), wfi (t), Ffi (t)}

Then, Υxij (t) can be calculated by (36–38) (see Lemma 2),
while Υcij (t) can be calculated by (41–42) (see Lemma 3)
and (46). In this case, Υij (t) is calculated by:

τd −1
τd −1

Υ (t) = HAdii Ξij (t − τij )[HAdjj ]T + Υ̂ij (t)

 ij

 Υ̂ij (t) = (1 − δ1,τ )(1 − δ1,τ )Σi Υx (t)ΣT
di
dj
ij
j
, (52)
τdj −1 T
c

(t)[H
]
+(1
−
δ
)Σ
Υ

1,τdi
i ij
Adj


τd −1

+(1 − δ1,τdj )HAdii [Υcji (t)]T ΣT
j
where δ1,τdi is defined in (4), and

τd −2
∆

 Σ1i = [Hdi , HAdi Hdi , · · · , HAdii Hdi ]


τd −2
∆

Σ2i = [H̄Adi , HAdi H̄Adi , · · · , HAdii H̄Adi ]
.
τ
−2
∆
d

Σ3i = [In , HAdi , · · · , HAdii ]




∆
Σi = [Σ1i Σ2i Σ3i ]

form of Ξij (t), one of the key issues is how to obtain the
relationship between E{x̃ci (t − di − 1)[x̃cj (t − dj − 1)]T } and
Ξij (t), which has been solved by Lemma 6.
According to the results in Lemmas 1–6, the recursive
form of Ξij (t) in (33) will be given by Theorem 1.
Theorem 1 Define
 ∆

 Λi = E{Hi (t) ⊙ Hi (t)}
∆
.
(54)
Vi = E{Hi (t) ⊙ [In − Hi (t)]}


∆
Wi = E{[In − Hi (t)] ⊙ [In − Hi (t)]}
∆

Then, the local estimation error covariance matrix Ξii (t) =
E{x̃ci (t)[x̃ci (t)]T } for each CSE x̂ci (t) is given by:

Ξii (t) = Adi [Wi ⊗ (AΞii (t − di − 1)AT )](Adi )T
+Adi [Λi ⊗ Pii (t − di ) + Wi ⊗ Qw ](Adi )T
+Adi [Vi ⊗ (ΦKi (t − di )Γii (t − di − 1)AT
+GKi (t − di )Qw )](Adi )T
,(55)
T
+Adi [ViT ⊗ (AΓT
ii (t − di − 1)ΦKi (t − di )
(t − di ))](Adi )T
+Qw GT
Pdi Ki θ−1
T
+ θ=1 A Qw [Aθ−1 ]

where Pii (t − di ) and Γii (t − di − 1) are calculated by (7)
and (45) (see Lemma 4). On the other hand, the estimation
∆
error cross-covariance matrix Ξij (t) = E{x̃ci (t)[x̃cj (t)]T } is
given
by:

τd
τ di
Ξij (t) = HAd
Ξij (t − τij )(HAdjj )T + Ξ̂ij (t)


i


 Ξ̂ij (t) = HAd Υ̂ij (t)HT + Pmin{di ,dj } Aθ Qw (Aθ )T

i

Adj
θ=1

xj

T

(t
−
d
,
t
−
d
{Φ
+H
x
i
j )Hdj
d

i
i


T
F

+Ψij (t)HAdj + Φxi (t − di , dj , t)




T
w

+Φ

xi (t − di , t − dj − 1)H̄Adj }


T

+HAdi {Θw

xci (t − di − 1, t − dj − 1)HAdj
T
T
F
(56)
+Ψji (t)Hdj + Θxci (t − di − 1, dj , t)}


T
di T

H̄
)
+
C
(d
,
d
)Q
(A
Q
{δ
+
H̄

o j
i
w
Adi di ,dj w Adj


T T

+(Φw

xj (t − dj , t − di − 1)) Hdj



T T

+(Θw

xci (t − dj − 1, t − di − 1)) HAdj }


T
T
F

+(Φxj (t − dj , di , t)) Hdj




F
T T

+(Θ

xcj (t − dj − 1, di , t)) HAdj


T
dj
+Co (di , dj )A Qw H̄Adj

where δdi ,dj is determined by (4), and Co (di , dj ) is determined by (34); Hdi , H̄Adi , HAdi are defined by (43), while
τij , τdi and τdj are defined in (50). Υ̂ij (t) is calculated by
(52) (see Lemma 6), and Ψij (t) is calculated by (49) (see
F
Lemma 5); Θw
xci (t1 , t2 ), Θxci (t1 , t2 ) are calculated by (41–
x
F
42) (see Lemma 3), while Φxji (t1 , t2 ), Φw
xi (t1 , t2 ), Φxi (t1 , t2 )
are calculated by (36–38) (see Lemma 2). Moreover, the
relationship between the CSE x̂ci (t) and the DKFE x̂(t) is
Tr{P (t)} ≤ Tr{Ξii (t)}(i ∈ {1, 2, · · · , L}).

(53)

Proof: See A.5 in Appendix.
Remark 5. Notice that the structure of Υxij (t) consists
xj
F
F
of Φw
xi (t1 , t2 ), Φxi (t1 , t2 ), Φxi (t1 , g, t2 ) and Φxi (t1 , g, t2 ),
x
and thus Υij (t) can be calculated by Lemma 2. Meanwhile,
Υcij (t) can be calculated by Lemma 3 and (46), because
its structure consists of Γij (t1 , t2 )(t1 > t2 ), Θw
xci (t1 , t2 )
and ΘF
(t
,
g,
t
).
On
the
other
hand,
in
most
cases,
the
c
2
xi 1
delay di is not equal to dj . Thus, to design the recursive

(57)

Proof: See A.6 in Appendix.
From Theorem 1, Ξ(t) can be calculated by (55–56),
then the optimal weighting matrices Ω1 (t), · · · , ΩL (t) are
obtained by (31). Moreover, the computation procedures for
the DKFE x̂(t) can be summarized by Algorithm 1.
Remark 6. According to (55–56), each covariance matrix
Ξij (t) is independent of the measurement yi (t) and the LOE
x̂i (t). Thus, Ξij (t) can be calculated at the FC when the
initial values are given. In this case, only if each selected
ASC x̂si (t) ∈ Rri is sent to the FC, Algorithm 1 will be
implemented in practical applications. On the other hand,
when the communication delays and the number of local

8

Algorithm 1 For the given selection probabilities π~i i (~i =
1, · · · , ∆i ; i = 1, · · · , L) satisfying (22)
1: At each sink node:
2: for i := 1 to L do
3:
Calculate GKi (t) and Ki (t) by (6) and (7);
4:
Calculate the LOE x̂i (t) by (5);
5:
Generate the binary variables σ~i i (t)(~i = 1, · · · , ∆i )
satisfying the categorical distribution, then determine
the selected ASC x̂si (t) by (13);
6: end for
7: At the FC:
8: for i := 1 to L do
9:
Calculate GKi (t) by (6);
10:
Calculate the CSE x̂ci (t) by (15);
11:
for j := i to L do
12:
Calculate Pij (t) by (6–8);
13:
Calculate Γij (t),Ψij (t),Υij (t) by (45), (49), (52);
14:
Calculate Ξij (t) by (55–56);
15:
end for
16: end for
17: Calculate Ω1 (t), Ω2 (t), · · · , ΩL (t) by (31);
18: Calculate the DKFE x̂(t) by (28).

estimates increase slightly, the computational complexity
of Algorithm 1 will be high. In this case, the steady-sate
DKFE with time-invariant weighting matrices can reduce
the amount of computation. Therefore, to obtain the steadysate DKFE, we should find the stability conditions satisfying
the following two points: i) The covariance matrix of the
recursive DKFE converges to a positive-definite matrix; ii)
The limit of the covariance matrix P (t) is independent of the
initial values. Following this idea, the steady-state DKFE will
be derived in the next section.
IV. S TABILITY A NALYSIS

FOR THE

DKFE

The stability condition and the steady-state DKFE will be
given in this section.
A. Stability Condition of Each Local CSE
The estimation performance of each CSE x̂ci (t) will be
discussed in this subsection. First, it is considered that the
ith subsystem satisfies
p
(A, Qw ) is stabilizable and (A, Ci ) is detectable.(58)

When the condition (58) holds, it is well known that the
estimation error covariance matrix Pii (t) in (7) will converge
from any initial conditions Pii (0) > 0 to the unique positive
semi-definite solution Pii . This means that
lim Pii (t) = Pii , lim ΦKi (t) = ΦKi , lim Ki (t) = Ki (59)

t→∞

t→∞

t→∞

where the limits Pii , ΦKi and Ki are independent of the
initial values. Moreover, ΦKi is a stable matrix. Thus, there
must exist an integer NPi > 0 such that, for t > NPi , the
estimation error system (94) reduces to:
x̃i (t) = ΦKi x̃i (t − 1) + GKi w(t − 1) − Ki vi (t). (60)
Then, it follows from (60) that
x̃i (t + 1) = ΦdKii+1 x̃i (t − di ) + ζio (t),

(61)

where
Pd +1
i −1
ζio (t) = αii =1 Φα
Ki GKi w(t − αi + 1)
Pdi
i −1
− αi =0 Φα
Ki GKi Ki vi (t − αi + 1).

(62)

Meanwhile, it is derived from (25) and (60) that

x̃ci (t + 1) = Ai1 (t)x̃ci (t − di ) + Ai2 (t)x̃i (t − di ) + ζic (t)(63)
where

Ai1 (t) = Adi [In − Hi (t − di + 1)]A




 Ai2 (t) = Adi Hi (t − di + 1)ΦKi
.(64)
ζic (t) = Adi Hi (t − di + 1)GKi w(t − di )

di

H
(t
−
d
+
1)w(t
−
d
)
+
F
(d
,
t
+
1)
+A

i
i
i
w
i


−Adi Hi (t − di + 1)vi (t − di + 1)
∆

∆

Define ξ̃i (t) = col{x̃ci (t), x̃i (t)}, ζi (t) = col{ζic (t), ζio (t)}.
Then, combining (61) and (63) yields that
ξ˜i (t + 1) = Ai (t)ξ˜i (t − di ) + ζi (t),

(65)

where
Ai (t) =




Ai1 (t) Ai2 (t)
.
0 ΦdKii+1

(66)

It is concluded from the definition of ξ̃i (t) that if the
∆
covariance matrix Ξ̃ξi (t) = E{ξ˜i (t)[ξ˜i (t)]T } converges to
the unique matrix Ξ̃ξi , there must exist the unique limit of
the estimation error covariance matrix Ξii (t) for the ith CSE.
Lemma 7 Define
∆

f (B) = E{AT
(67)
i (t)BAi (t)},


B11 B12
. Then, f (B) is calculated by:
where B =
B21 B22


B̃11 B̃12
,
(68)
f (B) =
B̃21 B̃22
where

B̃11




B̃12







 B̃
21





B̃22








= AT {Wi ⊗ [(Adi )T B11 Adi ]}A
= AT {ViT ⊗ [(Adi )T B11 Adi ]}ΦKi
di +1
+AT [In − Hi ](Adi )T B12 ΦK
i
T
di
di T
= ΦKi {Vi ⊗ [(A ) B11 A ]}A
(69)
+[ΦdKii+1 ]T B21 Adi [In − Hi ]A
di
di T
= ΦT
Ki {Λi ⊗ [(A ) B11 A ]}ΦKi
di +1 T
di
+[ΦKi ] B21 A Hi ΦKi + [ΦdKii+1 ]T B22 ΦdKii+1
di +1
di T
+ΦT
Ki Hi (A ) B12 ΦKi

where the probability selection matrix Hi is given by (23),
and Wi , Λi , Vi are given by (54). Moreover, for any matrices
B1 ,B2 , there will be:
f (B1 + B2 ) = f (B1 ) + f (B2 ).

(70)

Proof: (68) can be obtained from (16), (21), Lemma 1
and the definition of f (B), while (70) is derived from (68).
This completes the proof.
Based on Lemma 7, the delay-dependent stability condition of the CSE x̂ci (t) will be given in Theorem 2.
Theorem 2 For the communication delay di and selection
probabilities π~i i (~i = 1, · · · , ∆i ) in (68), if there exist Di >
0, Xi , Yi , Zi and Si > 0 such that


Xi Yi
≥ 0,
(71)
YiT Zi


Mi (1, 1)
−Yi − di Zi Ai
< 0, (72)
Mi =
−YiT − di AT
i Zi f (Di ) + di f (Zi ) − Si

9

where Mi (1, 1) = −Di + Xi + YiT + Yi + di Zi + Si and
Ai = E{Ai (t)}, while f (Di ) and f (Zi ) are calculated by
(68) in Lemma 7, then the covariance matrix Ξii (t) (55)
converges to the unique matrix, i.e.,
lim Ξii (t) = Ξii ,

t→∞

(74)

where fˆ(Di ) is calculated by f (Di ) (i.e., (68) in Lemma 7)
for di = 0, then the corresponding covariance matrix Ξii (t)
will converge to the unique steady-state value. Meanwhile,
it is concluded from the result in [30] that if
λmax (AT (In − Hi )A) < 1,

(75)

where Hi is determined by (23), then the lim Tr{Ξii (t)}
t→∞
will exist for di = 0. It should be pointed out that, under the
condition (75), one cannot prove the following results: a) the
limit of the covariance matrix Ξii (t) exists, and b) the limit
of Ξii (t) is independent of the initial conditions. Since the
results (a) and (b) are necessary before deriving the steadysate DKFE, the steady-state fusion estimator cannot be given
under the condition (75). Moreover, when only considering
the convergence of the sequence {Tr{Ξii (t)}}, the condition
(74) has less conservatism than the condition (75). This is
because the relaxation technique of matrix trace inequality is
introduced to derive (75), but the condition (74) is derived
by the stability theory without any relaxation. This result has
been demonstrated by Example 1.
B. Steady-State DKFE Design
According to (32) and (33), the stability of the DKFE is
also dependent on each estimation error cross-covariance matrix Ξij (t) (see (56)). Thus, the convergence of the sequence
{Ξij (t)} will first be discussed in this subsection, and then
combining Theorem 2 leads to the delay-dependent stability
condition of the DKFE and the steady-state DKFE (SDKFE).
The main result will be presented in Theorem 3.
Theorem 3 Consider the CPSs (1–2) under the condition (58), if the selection probabilities π~i i (~i =
1, 2, · · · , ∆i , i = 1, 2, · · · , L) and the communication
delays di (i = 1, 2, · · · , L) satisfy (71), (72), and
ρ(Adi (In − Hi )A) < 1(i = 1, 2, · · · , L),

lim P (t) = P

t→∞

(73)

and the limit Ξii is independent of the initial values.
Proof: See A.7 in Appendix.
Remark 7. When di 6= 0, it is calculated from (54) and
(55) that Tr{Adi [Vi ⊗ (ΦKi (t − di )Γii (t − di − 1)AT +
6 0 and Tr{Adi [ViT ⊗ (AΓT
GKi (t − di )Qw )](Adi )T } =
ii (t −
T
di T
di − 1)ΦKi (t − di ) + Qw GT
)
}
=
6
0,
which
(t
−
d
))](A
i
Ki
are different from the results (101) and (102) in [30]. This
implies that the stability condition for each CSE is difficult to
be obtained by the derivation method based on the property
of the operator Tr{•} in [30]. In contrast, by adopting a novel
derivation idea in this paper, the stability conditions (71) and
(72) are linear matrix inequalities (LMIs), and thus they can
be verified by resorting to Matlab LMI Toolbox [42].
Remark 8. When di = 0, the Lyapunov function candidate
can be chosen as Vξi (t) = E{ξiT (t)Di ξi (t)}. Then it is
concluded from the similar derivation of Theorem 2 that if
there exists Di > 0 such that
fˆ(Di ) − Di < 0,

where ρ(•) denotes the spectral radius, and Hi is given by
(23), then the fusion estimation error covariance matrix P (t)
(32) will converge to the unique matrix P , i.e.,

(76)

(77)

with P independent of the initial values. Moreover, the
steady-state weighing matrices Ωi (i = 1, 2, · · · , L) are
calculated by
[Ω1 , Ω2 , · · · , ΩL ] = P −1 IaT Ξ−1 ,

(78)

where lim Ξ(t) = Ξ, and the limit Ξ is independent of the
t→∞
initial values. In this case, the SDKFE x̂s (t) at the FC side
is given by:
XL
Ωi x̂ci (t),
(79)
x̂s (t) =
i=1

x̂ci (t)

where
is calculated by (15).
Proof: See A.8 in Appendix.
According to Theorem 3, the computation procedures for
the SDKFE x̂s (t) can be summarized by Algorithm 2.
Algorithm 2 For the given selection probabilities π~i i (~i =
1, · · · , ∆i ; i = 1, · · · , L) satisfying (22)
1: Determine the weighting matrices Ωi (i = 1, · · · , L) by
(78);
2: At each sink node:
3: for i := 1 to L do
4:
Calculate the LOE x̂i (t) by (5);
5:
Generate the binary variables σ~i i (t)(~i = 1, · · · , ∆i )
satisfying the categorical distribution, then determine
the selected ASC x̂si (t) by (13);
6: end for
7: At the FC:
8: for i := 1 to L do
9:
Calculate the CSE x̂ci (t) by (15);
10: end for
11: Calculate the SDKFE x̂s (t) by (79).
Remark 9. It has been proved in Theorem 3 that when
the conditions (71), (72) and (76) hold, the estimation error
covariance matrix P (t) can converge to the unique steadystate values for any initial conditions. Thus, the steadystate weighting matrices (78) can be obtained off-line by
implementing the Step 7–Step 16 of Algorithm 1. It is noted
that the computational complexity of the SDKFE obtained by
Algorithm 2 is much lower than that of the DKFE obtained
by Algorithm 1.
Notice that the stability condition in Theorems 2–3 are
dependent on the communication delays and the selection
probabilities of dimensionality reduction. Since each communication delay is determined by the property of communication channel, it is difficult to adjust the parameter di to
satisfy the stability condition. In this case, from the result
(57) in Theorem 1 and Theorems 2–3, how to determine the
selection probabilities (22) such that the MSE of the DKFE
is bounded or convergent will be presented in Theorem 4.
Theorem 4 For the CPSs (1–2), when each communication
delay di is known in prior, two probability criteria to determine the dimensionality reduction strategy are presented as
follows:

10

(C.1): To guarantee that the MSE of the DKFE is bounded,
one needs to determine one group of the selection probability
∆
i
} in (24) by (71–72).
ςi = col{π1i , · · · , π∆
i
(C.2): To guarantee the existence of the SDKFE, one needs
to determine the L selection probabilities ςi (i = 1, · · · , L)
in (24) by (71–72) and (76).

TABLE I
C OMPARISON OF THE RESULT IN T HEOREM 2 AND THE RESULT IN [30]

V. N UMERICAL E XAMPLES
In this section, two illustrative examples are presented
to show the advantage and effectiveness of the proposed
dimensionality reduction fusion estimation methods.
Example 1: Consider a CPS (1) with the following system
parameters [18]:




1.25 0
20 0
A=
, Qw =
,
(80)
1 1.1
0 20
where the parameters of the first measurement equation in
(2) are taken as:
C1 = [1 0], Qv1 = 2.5.

γ11

Inequality (83)

Inequality (84) derived by [30]

0

False

False

0.1

False

False

0.2

False

False

0.3

False

False

0.4

True

False

0.5

True

False

0.6

True

False

0.7

True

False

0.8

True

False

0.9

False

False

1.0

False

False

2000

(81)
1500

Therefore, it is concluded from Theorem 2 that the limit
of Tr{Ξ11 (t)} exists, however, the condition (75) derived
by [30] does not hold for this example. Moreover, when
choosing γ11 from 0 to 1, the effectiveness of the conditions
(74) and (75) is shown in Table I, Fig.3 and Fig.4. It can be
seen from Fig.3 and Fig.4 that the sequence of Tr{Ξ11 (t)}
is convergent when γ11 ∈ {0.4, 0.5, 0.6, 0.7, 0.8}. This result
can be directly obtained by the judgement condition (74),
however, the judgment condition (75) derived by [30] is invalid for the above cases. Therefore, the judgement condition
(74) in this paper has much less conservatism than the result
in [30], and thus is applicable to more fusion systems under
the dimensionality reduction.

(t)}

!11=0.7

!11=0.8

Tr{

11

40

60

80

100
120
t/step

140

160

180

200

Fig. 3. The estimation performances of the CSE (i.e., Tr{Ξ11 (t)}) with
different selection probabilities γ11 .

Example 2: Consider a CPS described by power grid
with a 4-bus model for the distribution test feeders [43]. To
monitor the work status of the grid, two sink nodes collect
their sensor measurements, and the local estimates computed
by the sink nodes are transmitted to the FC (e.g., monitoring
center or control center). According to the continuous-time
smart grid system in [16], and setting the sampling time
T0 = 10−4 s, the discretized system matrix in (1) is given

!10

40
Tr{! 11(t)}

4
2

!11=0

0

50

100

150

200

4
2

!10

"1

!11=0.1

0

50

t/step

4
2

!10

21

!11=0.2

0

50

100

150

200

4
2

!10

1
0

!10

50

!11=0.#

t/step

100

150

200

150

200

t/step

10

100

200

!11=0."

0

10

50

150

10

t/step

2

100
t/step

Tr{! 11(t)}

λmax (AT (I2 − H1 )A) = 1.5887 > 1.

20

Tr{! 11(t)}

while it is calculated from (80) and (83) that

!11=0.6

500

Tr{! 11(t)}

where 0 ≤ γ11 ≤ 1.
To demonstrate the advantage of the designed stability
condition, it is assumed that there is no communication delay
for this example, and the selection probability γ11 is taken as
γ11 = 0.5. Then, by using LMI Toolbox in Matlab to solve
the inequality matrix (74), one has


1.3284 0.1730 −0.0731 −0.0395
 0.1730 0.3727 0.0009 −0.0782 

D1 = 
 −0.0731 0.0009 1.0315 −0.5204  ,
−0.0395 −0.0782 −0.5204 1.9647

!11=0.5

1000

Tr{! 11(t)}

According to the dimensionality reduction strategy, it is
considered in this example that only one component of x̂1 (t)
is allowed to be transmitted to the FC at each time step. In
this case, it is calculated from (54) that

 Λ1 = diag{γ
11 , 1 −γ11 }, W1 = diag{1 − γ11 , γ11 }

0
γ11
(83)
, H1 = diag{γ11 , 1 − γ11 }
 V1 =
1 − γ11 0

!11=0.4

Tr{! 11(t)}

Then, √it is √calculated from (80) and (81) that
rank{[ Qw A Qw ]} = 2 and rank{col{C1 ,C1 A}} = 2,
which means that (58) holds. In this case, one has by (59)
that


0.4760 −0.8573
GKi =
.
(82)
0.0314 0.0376

150

200

!10

18

!11=1

5
0

50

100
t/step

Fig. 4. The estimation performances of the CSE (i.e., Tr{Ξ11 (t)}) with
different selection probabilities γ11 .

11

by:

3


1.0156 0.0139 0.0457 0.0971
 −0.0353 0.9997 −0.0008 −0.0017 

A=
 −0.0526 −0.0448 0.9625 −0.0797  ,
−0.008 −0.0505 −0.0903 0.9011


2"5

(84)
2
1"5

where λmax (A) = 1.0441 > 1 means that this 4-bus smart
grid system is unstable, and the covariance of the process
noise is taken as:


0.04 0.1 0.06 0.08
 0.1 0.25 0.15 0.2 

Qw = 
(85)
 0.06 0.15 0.09 0.12  .
0.08 0.2 0.12 0.16

Then, the measurement matrices in (2) are given by




1010
1001
0 1 0 0
1 0 1 0



C1 = 
 0 1 1 0  , C2 =  0 0 0 1  ,
1010
1010

(86)

which means that the measurement information on the
fourth component of “x(t)” cannot be obtained by the
first sink node, while the measurement information on the
second component of “x(t)” cannot be obtained by the
second sink node. The covariances of vi (t)(i = 1, 2)
in (2) are taken as Qv1 = diag{0.9, 0.6, 0.9, 0.4} and
Qv2 = diag{0.3, 0.4, 0.5, 0.2}, respectively. Then it is calcu2
3
lated from (84–86) that√rank(col{C
√ i , Ci2A,
√ Ci A ,3C
√i A }) =
4 (i = 1, 2) and rank([ Qw , A Qw , A Qw , A Qw ]) =
4, which means that the condition (58) holds. Thus, the limits
in (59) exist, and one has



0.7915 −0.0778 −0.1333 0.0886



 −0.3288 0.6887 −0.4510 0.0035 




ΦK 1 = 


 −0.1573 −0.2237 0.6459 −0.0659 



 −0.2790 −0.2490 −0.4294 0.9001  . (87)
0.7324 0.025 −0.0683 −0.0919



 −0.6848 1.0244 −0.2209 −0.4579 



 ΦK 2 = 


 −0.4547 −0.0287 0.5831 −0.1738 



−0.5093 −0.0291 −0.3937 0.6284
20

#r$% 11(t)&

1

0

10

20

0
20

DKFE for x 2(t)

30

40

50

60

70

t/step

50

60

70

10

DKFE for x 3(t)

0
60

70

-10
20

30

40

50

t/step

The trajectories of the DKFE x̂(t) and the state “x(t)”.

For this example, according to the dimensionality reduction strategy, it is considered that only two components of x̂i (t) are allowed to be transmitted to the FC for
satisfying the finite bandwidth, and thus r1 = r2 = 2

70

0

!0

100

(88)

Hi (t) = diag{σ1i (t) + σ2i (t) + σ3i (t), σ1i (t) + σ4i (t)
,(89)
+σ5i (t), σ2i (t) + σ4i (t) + σ6i (t), σ3i (t) + σ5i (t) + σ6i (t)}
where σ~i i (t)(~i = 1, 2, 3, 4, 5, 6) are determined by (11–12),
and each stochastic process {σ~i i (t)} obeys the categorical
distribution. To determine the signal x̂si (t) (see (13)), the
selection probabilities in (22) are taken as follows:
 1
 π1 = 0.3, π21 = 0.2, π31 = 0.1, π41 = 0.1
π 1 = 0.1, π61 = 0.2, π12 = 0.2, π22 = 0.1 .
 52
π3 = 0.2, π42 = 0.1, π52 = 0.3, π62 = 0.1

(90)

Thus, the selection probability matrices H1 and H2 (see (23))
are given by:
H1 = diag{0.6, 0.5, 0.5, 0.4}
.
H2 = diag{0.5, 0.6, 0.3, 0.6}

(91)

When each selected ASC x̂si (t) is transmitted to the FC, the
communication delays are taken as d1 = 1 and d2 = 2. In
this case, it is calculated from (84) and (91) that


x 3(t)

DKFE for x 3(t)

t/step

Fig. 5.

40

20

-20
-30
50

30

t/step
x 3(t)

-10

-10
20

60

Then it follows from (16) and (88) that

10
x 2(t)

50
t/step

 i
 H1 = diag{1, 1, 0, 0}, H2i = diag{1, 0, 1, 0}
H i = diag{1, 0, 0, 1}, H4i = diag{0, 1, 1, 0} .
 3i
H5 = diag{0, 1, 0, 1}, H6i = diag{0, 0, 1, 1}



-5

40

and ∆1 = ∆2 = 6. In this case, the diagonal matrices
H~i i (t)(i = 1, 2; ~i = 1, 2, 3, 4, 5, 6) in (16) are given by

0

5

30

Fig. 6. Comparison of estimation performance for the CSEs, DKFE and
ODKFE.

5
DKFE for x 1(t)

#r$'o(t)&

#r$'(t)&

0"5

x 1(t)

15

#r$% 22(t)&

ρ(A(I4 − H1 )A) = 0.5759 < 1
,
ρ(A2 (I4 − H2 )A) = 0.6631 < 1

(92)

which means that the condition (76) holds in Theorem 3.
Meanwhile, by using LMI Toolbox in Matlab, the variables
Di , Xi , Yi , Zi and Si (i = 1, 2) are obtained by solving the
matrix inequalities (71) and (72), i.e., the conditions (71–
72) hold for the two local CSEs with different selection
probabilities and communication delays. Under this case,
it is concluded from Theorem 3 that the fusion estimation
covariance matrix P (t) for this example converges to a
unique matrix, and the SDKFE exists. Then, implementing
Algorithm 1 obtains the steady-state weighting matrices as

12

follows:



0.6254 0.0921 0.3294 0.044



 0.0585 0.7874 0.2587 0.2107 




Ω1 = 


 0.1654 0.0271 0.6857 0.2670 



 0.0065 0.0765 0.2257 0.6729
 . (93)
0.3746 −0.0921 −0.3294 −0.044



 −0.0585 0.2126 −0.2587 −0.2107 




Ω2 = 


 −0.1654 −0.0271 0.3143 −0.2670 



−0.0065 −0.0765 −0.2257 0.3271

Thus, the SDKFE x̂s (t) for this example is obtained by
substituting (93) into (79).
By using Algorithm 1, the trajectories of the DKFE “x̂(t)”
and the state “x(t)” are plotted in Fig.5, which shows that
the designed DKFE is able to estimate the original state
“x(t)” well. Meanwhile, let Po (t) denote the original DKFE
(ODKFE) under the dimensionality reduction when there are
no communication delays between the sink nodes and the FC.
Then, the estimation performances (assessed by the trace of
the estimation error covariance matrix) of the local CSEs,
DKFE and ODKFE are shown in Fig.6. It is seen from this
figure that the estimation performance of the DKFE is better
than that of each CSE at each time-step, which is in line with
the result (57). However, the estimation performance of the
DKFE is worse than that of the ODKFE, which implies that
the communication delays can affect the fusion estimation
performance. Moreover, it is known from the this figure that
the MSEs of the DKFE and CSEs all converge to the steadysate values, which accords with the results (73) and (77).
To demonstrate the effectiveness of the SDKFE for this
example, the matrix 2-norms of P (t) and Ωi (t)(i = 1, 2)
are shown in Fig.7 under different initial values. It is seen
from this figure that ||P (t)||2 and ||Ωi (t)||2 (i ∈ {1, 2}) can
converge to the unique steady-state values under different
initial values, which is in line with the result in Theorem
∆
3. On the other hand, let Eri (t) = x̂(i, t) − x̂s (i, t), where
x̂(i, t) represents the ith component of the DKFE x̂(t), and
the meaning of x̂s (i, t) is the same as that of x̂(i, t). Then,
implementing Algorithms 1–2, the trajectories of Eri (t)(i =
1, 2, 3, 4) are depicted in Fig.8, where the measurement
sequences {yi (t)(i = 1, 2)} are the same when computing
the DKFE x̂(t) and the SDKFE x̂s (t). It is shown from this
figure that the errors between the DKFE and SDKFE will
converge to zero as t increases, which is in line with the
property of the SDKFE. It should be pointed out that the
SDKFE is much easier to implement as compared with the
DKFE in practical applications.
VI. CONCLUSIONS
As CPSs are being widely integrated in various critical
infrastructures and running on wired or wireless communication networks, however, bandwidth constraints and communication delays are usually unavoidable. Notice that state
estimation plays an essential role in the monitoring and
supervision of CPSs, and its importance has made the robustness and estimation performance a major concern. Therefore,
to guarantee the satisfactory estimation performance in CPSs,
the distributed dimensionality reduction fusion estimation
problem with communication delays has been studied in this
paper. Based on the stochastic dimensionality reduction strategy, a mathematical model was proposed to establish the relationship between the dimensionality reduction and communication delays, and then the recursive DKFE was obtained by

2.5
||P1(t)||2
||P2(t)||2

2

|| 1
(t)||2
1
|| 2
(t)||2
1
|| 1
(t)||2
2

1.5

|| 2
(t)||2
2

1

0.5
5

10

15

20

25

30
35
t/step

40

45

50

55

60

Fig. 7. The trajectories of ||P i (t)||2 , ||Ωi1 (t)||2 and ||Ωi2 (t)||2 (i = 1, 2),
where P i (t), Ωi1 (t) and Ωi2 (t) represents the covariance matrix P (t) and
the weighting matrices Ω1 (t), Ω2 (t) under different initial values for i 6= j.

0.4
$%1(t)

0.2

$%2(t)

$%3(t)

$%4(t)

0
#0.2
#0.4
#0.6
#0.!
#1
10

Fig. 8.

20

30

40

50

60
t/step

0

!0

"0

100

The trajectories of Eri (t) (i = 1, 2, 3, 4).

resorting to the optimal weighted fusion criterion. A delaydependent and probability-dependent condition, which can be
easily judged by using Matlab LMI Toolbox, was derived for
the DKFE such that the fusion estimation error covariance
matrix P (t) converges to a unique steady-state matrix. This
result is very important to derive the SDKFE, and the computational complexity of the SDKFE is much lower than that
of the DKFE. Meanwhile, when the communication delay
di is known in advance, the selection probability criterion
to determine the dimensionality reduction strategy has also
been presented. Moreover, it has been shown that the stability
condition in this paper has less conservatism than the exiting
ones. Finally, two examples were given to demonstrate the
advantage and effectiveness of the proposed methods.
Along this line of work, the design of distributed dimensionality reduction fusion estimator with communication
delays for nonlinear CPSs is one of our future works.
APPENDIX
A.1: The proof of Lemma 2
Proof: It follows from (1) and (5) that
x̃i (t) = ΦKi (t)x̃i (t − 1)
+GKi (t)w(t − 1) − Ki (t)vi (t)
For t1 ≥ t2 , it is derived from (94) that

(94)

13


Q
t1 −t2 −1
(t
−
ϕ
)
x̃i (t2 )
Φ
x̃i (t1 ) =
1
i
K
i
ϕi =0

Pt1 −t2 nQαi −2
+ αi =1
ϕi =0 ΦKi (t1 − ϕi )
×GKi (t1 − αi + 1)w(t1 − αi )}

Pt1 −t2 −1 nQαi −1
− αi =0
(t
−
ϕ
)
Φ
1
i
K
i
ϕi =0
×Ki (t1 − αi )vi (t1 − αi )}

Meanwhile, it follows from (25) that
(95)

On the other hand, it is concluded from (3) and the geometric
meaning of x̃i
(t) that
x̃i (t1 )⊥w(t2 )(t2 ≥ t1 )



x̃i (t1 )⊥vi (t2 )(t2 > t1 )
(96)
x̃

i (t1 )⊥vj (t2 )(i 6= j, ∀t1 , t2 )


w(t1 )⊥vi (t2 )(∀i, t1 , t2 )

One has by (34) that Co (t1 , t2 − 1) = 1 for t1 ≥ t2 , and thus
the results (36)-(37) are obtained by (95) and (96). Moreover,
the result (38) is directly obtained by the definitions of
F
Φw
xi (t1 , t2 ) and Φxi (t1 , g, t2 ) in (35), while the result (39)
is obtained from (3) and (26).
A.2: The proof of Lemma 3
Proof: Let us define
(
∆
HAdi (t) = Adi [In − Hi (t − di )]A
∆

Hdi (t) = Adi Hi (t − di ), H̄Adi (t) = Adi − Hdi (t)

(97)

For t1 ≥ t2 , it follows from (25) that
Pχi (t1 ,t2 )−1
~
x̃ci (t1 ) = ~=0
(t1 ) − di )
{Hdi (fi~ (t1 ))x̃i (fio
Q~−1
µ
~+1
~
(t1 ))w(fio
(t1 ))}
+ µ=0 HAdi (fi◦ (t1 ))H̄Adi (fio
Pχi (t1 ,t2 )−1
~−1
~
+ ~=1
{HAdi (fio (t1 ))Fw (di , fio (t1 ))} (98)
Qχi (t1 ,t2 )−1
χ (t ,t )
µ
(t1 ))x̃ci (fioi 1 2 (t1 ))
HAdi (fio
+ µ=0
0
(t1 ))
+(1 − δt1 ,t2 )Fw (di , fio
where fi (t) and χi (t1 , t2 ) are defined in (40). Notice that
χ (t1 ,t2 )

fi◦i

χ (t1 ,t2 )−1

(t1 ) < fi◦i

0
(t1 ) < · · · < fi◦
(t1 )

(99)

Γij (t) = E{x̃i (t)[x̃cj (t − di − 1)]T }HT
Adj
T
T
+E{x̃i (t)x̃T
j (t − dj )}Hdj + E{x̃i (t)Fw (dj , t)} (104)
T
T
+E{x̃i (t)w (t − dj − 1)}H̄Adj
Therefore, the result (45) is derived from (103–104) and
Lemma 2. Meanwhile, according to (95), (46) can be derived
from the similar derivation of (45).
A.4: The proof of Lemma 5
Proof: (48) can be derived from (103). On the other
hand, it follows from (25) that
x̃cj (t − dj −n
1) = x̃cj (fj (t))

Pηij −1 Qκ−1
υ
(f
(t))
Hd (f κ (t))
= κ=1
H
Ad
j
jo
υ=1
nQ j jo

P
ηij −1
κ−1
υ
κ
×x̃j (fjo
(t) − dj )} + κ=1
υ=1 HAdj (fjo (t))
(105)
κ
(t))w(fjκ+1
×H̄Adj (fjo
(t))} 
◦
Pηij −1 Qκ−1
υ
κ
+ κ=1
υ=1 HAdj (fjo (t)) Fw (dj , fjo (t))
Q

ηij −1
ηij
κ
c
+
κ=1 HAdj (fjo (t)) x̃j (fjo (t))

where Hdj (t), HAdj (t) and H̄Adj (t) are defined in (97).
Meanwhile, it follows from the similar derivation of (103)
that
η

E{x̃i (t
− di )[x̃cj (fjoij (t))]T }
Q

(106)
ηij (dj +1)−1−di
η
=
ΦKi (t − di − ϕi ) Γij (fjoij (t))
ϕi =0

Therefore, (49) is obtained from (100), (105), (106) and
Lemma 2.
A.5: The proof of Lemma 6
Proof: To establish the relationship between Υij (t) and
Ξij (t), the least common multiple of di + 1 and dj + 1 is
introduced, and thus one has
τd

τd

fi◦ i (t) = fi◦ j (t) = t − τij
(107)
Then taking the statistical property of γℓi (t) into account
yields:
where fi (t) is defined in (40). On the other hand, for i 6= j,
( ~−1
!)
it
is concluded from the statistical property of Hi (t) that
Y
µ
~
E
(t1 ))
= H~Adi H̄Adi (100)
(t1 ))H̄Adi (fi◦
HAdi (fi◦
Υij (t) = E{x̃ci (fi (t))[x̃cj (fj (t))]T }
µ=0
(108)
= E{x̄ci (fi (t))[x̄cj (fj (t))]T }
where HAdi and H̄Adi are given by (43). Moreover, it is
where
concluded from (3), (25) and (96) that
Pτdi −1 κ−1
 c
κ
x̄ci (fi (t)) = κ=1
(t) − di )
H i Hdi x̃i (fio
 x̃i (t1 )⊥w(t2 )(t2 ≥ t1 − d1 )
Pτdi −1Adκ−1
κ+1
c
+
w(f
H̄
H
x̃ (t1 )⊥vi (t2 )(t2 > t1 − d1 )
(101)
io (t))
Adi Adi
(109)
Pτκ=1
 ci
di −1
κ
x̃i (t1 )⊥vj (t2 )(i 6= j, ∀t1 , t2 )
+ κ=1
Hκ−1
Adi Fw (di , fio (t))
τd −1
+HAdii x̃ci (t − τij )
Thus, the result (41) is derived from (98–101). On the
other hand, (42) is directly obtained from the definitions of Notice that when τ = 1, x̄c (f (t)) = x̃c (t − τ ). Then,
ij
di
i i
i
Θw
xci (t1 , t2 ) and Fw (g, t).
(109) can be written as:
A.3: The proof of Lemma 4
Proof: It follows from (95) that
Q

Pdj +1
dj
x̃i (t) =
αi =1
ϕi =0 ΦKi (t − ϕi ) x̃i (t − dj − 1) +
o

nQ
αi −2
(102)
(t
−
α
+
1)w(t
−
α
)
(t
−
ϕ
)
G
Φ
i
i
i
K
K
i
i
ϕi =0

Pdj Qαi −1
− αi =0
ϕi =0 ΦKi (t − ϕi ) Ki (t − αi )vi (t − αi )

Then one has by (101) and (102) that

E{x̃i (t)[x̃cj (t− dj − 1)]T }

Qdj
(103)
=
ϕj =0 ΦKi (t − ϕi ) Γij (t − dj − 1)

τ

−1

di
c
x̄ci (fi (t)) = (1 − δ1,τdi )Σi x̄w
fi (t) + HAdi x̃i (t − τij()110)

where x̄w
fi (t) and Σi are given by (51) and (53), respectively.
Therefore, (52) is derived from (108) and (110).
A.6: The proof of Theorem 1
Proof: It is concluded from (3), (96) and (101) that
x̃i (t − di )⊥Fw (di , t), x̃ci (t − di )⊥Fw (di , t)
(111)
where Fw (di , t) is defined by (26). Notice that

T
E{x̃ci (t − di − 1)x̃T
i (t − di )} = Ψii (t)
T
E{[In − Hi (t)] ⊙ Hi (t)} = Vi

(112)

14

Then, (55) is derived from Lemma 1, (25), (96), (101), (45),
(48) and (111–112). On the other hand, (25) is rewritten as:
x̃ci (t) = Hdi (t)x̃i (t − di ) + HAdi (t)x̃ci (t − di − 1)
(113)
+H̄Adi (t)w(t − di − 1) + Fw (di , t)
where Hdi (t), HAdi (t) and H̄Adi (t) are defined in (97).
Meanwhile, for i 6= j, one has by (16) and (21) that
E{Hi (t)Hj (t)} = E{Hi (t)}E{Hj (t)}}

(114)

Moreover, it is derived from (3) and (26) that

di T
E{wi (t − di − 1)FT
w (dj , t)}Co (dj , di )Qw (A )
Pmin{di ,dj } θ
T(115)
T
E{Fw (di , t)Fw (dj , t)} = θ=1
A Qw (Aθ )
where Co (dj , di ) is determined by (34). Therefore, (56)
is derived from (113–115) and the results in Lemmas 2–
6. Furthermore, at a particular time, the optimal fusion
estimation error covariance matrix of (28) can be calculated
by (32), while each local estimation error covariance matrix
of (15) is calculated by (55), then (57) is obtained from the
results in [7], [8].
A.7: The proof of Theorem 2
Proof: Consider the following stochastic system:
ξi (t + 1) = Ai (t)ξi (t − di ),

(116)

where Qζi = E{ζi (t)ζiT (t)}, and f (Ξ̃ξi (t−di )) is calculated
by (68) in Lemma 7. In what follows, we will prove that,
under the conditions (71–72), the sequence {Ξ̃ξi (t)} obtained
from (122) is convergent, and the limit is independent of the
initial values.
∆
Define ∆Ξ̃ξi (t) = Ξ̃ξi (t) − Ξ̃ξi (t − 1). Then, it is derived
from (122) that
∆Ξ̃ξi (t + 1) = f (Ξ̃ξi (t − di )) − f (Ξ̃ξi (t − di − 1)).(123)
Combining (70) in Lemma 7 yields that
∆Ξ̃ξi (t + 1) = f (∆Ξ̃ξi (t − di )).

(124)

Notice that the recursive form of (124) is the same as that
of (120), thus one has lim ∆Ξ̃ξi (t) = 0, which leads to
t→∞

lim Ξ̃ξi (t) = Ξ̃ξi .

t→∞

(125)

For the recursive equation (122), let Ξ̃1ξi and Ξ̃2ξi (t) denote
any matrices with different initial conditions, and define
∆
Ξ̂ξi (t) = Ξ̃1ξi (t) − Ξ̃2ξi (t). In this case, it is derived from
(122) and (124) that Ξ̂ξi (t + 1) = f (Ξ̂ξi (t − di )), whose
recursive form is similar to (120). Then, it is concluded that
lim Ξ̂ξi (t) = 0, which implies
t→∞

∆

where Ai (t) is defined by (66). Define ηi (t) = ξi (t) − ξi (t −
(126)
lim Ξ̃1 (t) = lim Ξ̃2ξi (t),
t→∞
t→∞ ξi
1). Then, one chooses a Lyapunov function candidate for the
system (116) as follows:
i.e., the limit Ξ̃ξi in (125) is unique. Moreover, for t > NPi ,
Pt−1
T
T
it follows from the definition of Ξ̃ξi (t) that
Vξi (t) = E{ξi (t)Di ξi (t)} + κ=t−di E{ξi (κ)Si ξi (κ)}
Pt
P−1
(117)
.


+ βi =−di κ=t+βi +1 E{ηiT (κ)Zi ηi (κ)}
Ξii (t) ΓT
(t)
ii
,
(127)
Ξ̃ξi (t) =
Γii (t) Pii
Notice that
(
P
where Pii is given by (59), while Γii (t),Ξii (t) are calculated
ξi (t − di ) = ξi (t) − tκ=t−di +1 ηi (κ)
Pt
.(118)
by
(45) and (55). In this case, it can be concluded from (125–
ξi (t + 1) = Ai (t)ξi (t) − Ai (t) κ=t−di +1 ηi (κ)
127) that lim Ξii (t) = Ξii , and the limit Ξii is independent
t→∞
From the similar derivation of Theorem 1 in [41], it can be of the initial
values.
derived from (116–118) that
A.8: The proof of Theorem 3
∆Vξi (t) = Vξi (t + 1) − Vξi (t)
T
T
≤ ξi (t){−Di + Xi + Yi + Yi + di Zi
Proof: When the CPSs (1–2) satisfy the condition (58),
+Si }ξi (t) + ξiT (t − di )E{AT
it is concluded from the result in [8] that
i (t)Di Ai (t) (119)
+di E{AT
i (t)Zi Ai (t)} − Si }ξi (t − di )
lim Pij (t) = Pij ,
(128)
+ξiT (t){−Yi − di Zi E{Ai (t)}}ξi (t − di )
t→∞
T
T
T
+ξi (t − di ){−Yi − di E{Ai (t)}Zi }ξi (t),
where Pij (t) is calculated by (8), and the limit Pij is
where Xi , Yi and Zi are required to satisfy (71). Under independent of the initial values. Then, it follows from
this case, according to Lyapunov stability theory (see [42]: Lemmas 2–3, (59) and (128) that
p.131), when the condition (72) holds, the system (116) is

x
x
mean-square stable. This means that the state covariance
lim Φxji (t, t − ε) = Φxji (ε)


t→∞

∆

F

matrix Ξξi (t) = E{ξi (t)ξiT (t)} of (116) converges to zero

lim ΦF
xi (t, g, t − ε) = Φxi (ε)
 t→∞

under the conditions (71) and (72) (i.e., lim Ξξi (t) = 0),
w
lim Φw
t→∞
,
(129)
xi (t, t − ε) = Φxi (ε)
t→∞
where


w
w

lim Θ c (t, t − ε) = Θxci (ε)


t→∞ xi

Ξξi (t + 1) = E{Ai (t)ξi (t − di )ξiT (t − di )AT
i (t)} (120)

F
 lim ΘF
c (t, g, t − ε) = Θxc (ε)
= f (Ξξi (t − di )).
i
t→∞ xi
On the other hand, one has by (16),(21),(96),(101) that
where these limits are independent of the initial values.
∆
x

F
Define Γ̂ij (t) = Φxji (t, t − dj )HT
dj + Φxi (t, dj , t) +
ξ˜i (t − di )⊥ζi (t)
(121)
,
T
w
Φxi (t, t − dj − 1)H̄Adj . Then, one has by (59) and (129)
E{Ai (t)ξ˜i (t − di )ξ˜iT (t − di )AT
i (t)} = f (Ξ̃ξi (t − di ))
that



where Ξ̃ξi (t) = E{ξ˜i (t)ξ˜iT (t)}. Then, it follows from (65)
 lim Qdj ΦK (t − ϕj ) = Φdj +1
and (121) that
i
ϕ
=0
Ki
j
t→∞
,
(130)
 lim Γ̂ij (t) = Γ̂ij
(122)
Ξ̃ξi (t + 1) = f (Ξ̃ξi (t − di )) + Qζi ,
t→∞

15

where Γ̂ij is independent of the initial values. Thus, there
must exist an integer NΓi (> NPi ) such that, for t > NΓi ,
(45) reduces to:
d +1

Γij (t) = ΦKji Γij (t − dj − 1)HT
Adj + Γ̂ij ,

(131)

where HAdj is defined by (43). In this case, it follows from
(131) that
d +1
~−~0
~−~0
Γij (t) = [ΦKji ]~−~0 Γij (fjo
(t))[HT
Adj ]
(132)
P~−~0 −1 dj +1 κ
κ
+ κ=0
[ΦKi ] Γ̂ij [HT
Adj ] ,

where fj (t) is defined in (40), and the variable ~0 is
determined by:
~∗
~0 = min{~∗ |fjo
(t) − NΓi ≥ 0},

(133)

which implies that t → ∞ ⇔ ~ → ∞. Meanwhile, it is
known from (59) and (76) that ΦKi and HAdj are stable
matrices, which means that
lim ΦtKi = 0, lim HtAdj = 0.

(134)

t→∞

t→∞

Then, it is concluded from (133) and (134) that

~0
~−~0
 lim [ΦdKj +1 ]~−~0 Γij (fjo
(t))[HT
=0
Ad ]
j

i

t→∞

~−~0 −1
 lim [ΦdKji+1 ]~−~0 −1 Γ̂ij [HT
=0
Adj ]

. (135)

~→∞

Therefore, one has by (132) and (135) that
lim Γij (t) = Γij .

(136)

t→∞
d +1

2
Let Γ1ij (t) = ΦKji Γ1ij (t − dj − 1)HT
Adj + Γ̂ij and Γij (t) =
d +1

1
ΦKji Γ2ij (t − dj − 1)HT
Adj + Γ̂ij , where Γij (t − dj − 1) and
Γ2ij (t − dj − 1) denote the different initial values. Then, it is
known from (136) that lim Γ1ij (t) = Γ1ij and lim Γ2ij (t) =
t→∞

t→∞

∆

Γ2ij . Meanwhile, defining ∆Γij (t) = Γ1ij (t) − Γ2ij (t) yields
d +1
~0
~−~0
(t))[HT
that ∆Γij (t) = [ΦKji ]~−~0 ∆Γij (fjo
, and
Adj ]
thus it follows from (135) that lim ∆Γij (t) = 0, i.e.,
t→∞

Γ1ij

= Γ2ij .

(137)

This implies that the limit Γij in (136) is independent of the
initial values.
According to the computation formulas of Ψij (t) and
Υ̂ij (t), it is obtained from (129) and (136) that
lim Ψij (t) = Ψij , lim Υ̂ij (t) = Υ̂ij ,

t→∞

t→∞

(138)

where Ψij and Υ̂ij are independent of the initial values.
Then, combining (129) and (138) yields that
lim Ξ̂ij (t) = Ξ̂ij ,

t→∞

(139)

where Ξ̂ij is independent of the initial values. Moreover,
from (139), there must exist an integer NΞi (> NΓi ) such
that Ξij (t) (56) reduces to:
τd

τd

Ξij (t) = HAdii Ξij (t − τij )[HAdjj ]T + Ξ̂ij (t > NΞi ).(140)
When the condition (76) holds, HAdi and HAdj in (140) are
stable matrices. In this case, the form of (140) is the same
as that of (131), and thus it is obtained from the similar
derivation of (136–137) that
lim Ξij (t) = Ξij ,

t→∞

and Ξij is independent of the initial values.

(141)

Therefore, when the conditions (71), (72) and (76) hold,
the result (77) can be obtained from (73) and (141). Moreover, the steady-state weighted matrices (78) can be derived
from (31), the definition of Ξ(t) (see (33)) and (77). Notice
that the results (77) and (78) have shown that the designed
DKFE is independent of the initial values, the SDKFE (79)
can be thus obtained.
R EFERENCES
[1] P.Y. Chen, S.M. Cheng, K.C. Chen, Information fusion to defend
intertional attack in internet of things, IEEE Internet of Things Journal,
vol. 1, no. 4, 2014, pp. 337-348.
[2] A. Miloslavov, M. Veeraraghavan, Sensor data fusion algorithms for
vehicular cyber-physical systems, IEEE Transactions on Paralle and
Distributed Systems, vol. 23, no. 9, 2012, pp. 1762-1774.
[3] X.R. Li, Y. Zhu, J. Wang, C. Han, Optimal linear estimation fusionpart I: unified fusion rules, IEEE Transactions on Information Theory,
vol. 49, no.9, 2003, pp. 2192-2208.
[4] J.A. Roecker, C.D. McGillen, Comparison of two-sensor tracking
methods based on state vector fusion and measurement fusion, IEEE
Transactions on Aerospace and Electronic Systems, vol. 24, no. 4,
1988, pp. 447-449.
[5] Y. Bar-Shalom, X.R. Li, T. Kirubarajan, Estimation with applications
to tracking and navigation, John Wilely and Sons, Inc., 2001.
[6] N.A. Carlson, Federated square root filter for decentralized parallel
processes, IEEE Transactions on Aerospace and Electronic Systems,
vol. 26, no. 3, 1990, pp. 517-525.
[7] S. Sun, Z. Deng, Muti-sensor optimal information fusion Kalman filter,
Automatica, vol. 40, 2004, pp. 1017-1023.
[8] Z. Deng, Y. Gao, L. Mao, Y. Li, G. Hao, New approach to information
fusion steady-state Kalman filtering, Automatica, vol. 41, 2005, pp.
1695-1707.
[9] Y. Wang, X.R. Li, Distributed estimation fusion with unavaliable crosscorrelation, IEEE Transactoins on Aerospace and Electronic Systems,
vol. 48, no.1, 2012, pp. 259-278.
[10] E. Song, J. Xu, Y. Zhu, Optimal distributed Kalman filtering fusion
with singular covariances of filtering errors and measurement noises,
IEEE Transactions on Automatic Control, vol. 59, 2014, pp. 12711282.
[11] W. Zhang, B. Chen, M.Z.Q. Chen, Hierarchical fusion estimation
for clustered asynchronous sensor networks, IEEE Transactions on
Automatic Control, vol. 61, no. 10, 2016, pp. 3064-3069.
[12] B. Chen, G. Hu, D.W.C. Ho, W. Zhang, L. Yu, Distributed robust
fusion estimation with application to state monitoring systems, IEEE
Transactions on Systems, Man, and Cybernetics: Systems, Article in
press, 2016, doi:10.1109/TSMC.2016.2558103.
[13] K.H. Johansson, G.J. Pappas, P. Tabuada, C.J. Tomlin, Special issue
on control of cyber-physical systems, IEEE Transactions on Automatic
Control, vol. 59, no. 12, 2014, pp. 3120-3121.
[14] A. Gokhale, M. P.Mcdonald, A cyber physical systems perspective on
the real-time and reliable dissemination of information in intelligent
transportation systems, Network Protocls and Algorithms, vol. 2, no.
3, 2010, pp. 116-136.
[15] J. Fink, A. Ribeiro, V. Kumar, Robust control for mobility and wireless
communication in cyber-physical systems with application to robot
teams, Proceedings of the IEEE, vol. 100, no. 1, 2012, pp. 164-178.
[16] H. Li, L. Lai, V. Poor, Multicast routing for decentralized control of
cyber physical systems with an application in smart grid, IEEE Journal
on Selected Areas in Communication, vol. 30, 2012, pp. 1097-1107.
[17] Y.F. Huang, S. Werner, J. Huang, N. Kashyap, V. Gupta, State
estimation in electric power girds: meeting new challenges presented
by requirements of the future grid, IEEE Signal Processing Magazine,
vol. 29, no. 5, 2012, pp. 33-43.
[18] X. Cao, P. Cheng, J. Chen, S.S. Ge, Y. Cheng, Y. Sun, Cognitive radio
based state estimation in cyber-physical systems, IEEE Journal on
Selected Areas in Communicaiton, vol. 32, no. 3, 2014, pp. 489-502.
[19] S. Deshmukh, B. Natarajan, A. Pahwa, State estimation over a lossy
network in spatially distributed cyber-physical systems, IEEE Transactions on Signal Processing, vol. 62, no. 15, 2014, pp. 3911-3923.
[20] B. Chen, G. Hu, W. Zhang, L. Yu, Information fusion estimation
for spatially distributed cyber-physical systems with communication
delay and bandwidth constraints, Proceedings of the American Control
Conference, Chicago, USA, 2015, pp. 5152-5157.
[21] B. Chen, D.W.C. Ho, G. Hu, L. Yu, Secure fusion estimation for
bandwidth constrained cyber-physical systems under replay attacks,
IEEE Transactions on Cybernetics, Article in press, 2017, Doi:
10.1109/TCYB.2017.2716115.
[22] Z. Wang, Y. Niu, Distributed estimation and filtering for sensor
networks, International Journal of System Science: Speical Issue, vol.
42, 2011, pp. 1421-1425.

16

[23] B. Chen, W. Zhang, G. Hu, L. Yu, Networked fusion kalman filtering
with multiple uncertainties, IEEE Transactions on Aerospace and
Electronic Systems, vol. 51, no. 3, 2015, pp. 2332-2349.
[24] B. Chen, G. Hu, D.W.C. Ho, L. Yu, Distributed covariance intersection
fusion estimation for cyber-physical systems with communication
constraints, IEEE Transactions on Automatic Control, vol. 61, no. 12,
2016, pp. 4020-4026.
[25] S. Sun, H. Lin, J. Ma, X. Li, Multi-sensor distributed fusion estimation
with applications in networked systems: a review paper, Information
Fusion, vol. 38, 2012, pp. 122-134.
[26] J.J. Xiao, A. Ribeiro, Z.Q. Luo, G.B. Giannakis, Distributed
compression-estimation using wireless sensor networks, IEEE Signal
Processing Magazine, vol. 23, no. 4, 2006, pp. 27-41.
[27] J. Fang, H. Li, Hyperplane-based vector quantization for distributed
estimation in wireless sensor networks, IEEE Transactions on Information Theory, vol. 55, no. 12, 2009, pp. 5682-5699.
[28] X. Shen, P. K.Varshney, Y. Zhu, Robust distributed maximum likelihood estimation with dependent quantized data, Automatica, vol. 50,
no. 1, 2014, pp. 169-174.
[29] D. Li, S. Kar, F.E. Alsaadi, S. Cui, Distributed Kalman filtering with
quantized sensing state, IEEE Transactions on Signal Processing, vol.
63, no. 19, 2015, pp. 5180-5193.
[30] B. Chen, W.A. Zhang, L. Yu, Distributed finite-horizon fusion Kalman
filtering for bandwidth and energy constrained wireless sensor networks, IEEE Transactions on Signal Processing, vol. 62, no. 4, 2014,
pp.797-812.
[31] H. Ma, Y.H. Yang, Y. Chen, K.J. Ray Liu, Q. Wang, Distributed state
estimation with dimension reduction preprocessing, IEEE Transactions
on Signal Processing, vol. 62, no. 12, 2014, pp. 3098-3110.
[32] B. Chen, W.A. Zhang, L. Yu, G. Hu, H. Song, Distributed fusion estimation with communication bandwidth constraints, IEEE Transactions
on Automatic Control, vol. 60, no. 5, 2015, pp. 1398-1403.
[33] I. D. Schizas, G. B. Giannkis, Z.Q. Luo, Distributed estimation using
reduced-dimensionlaity sensor observations, IEEE Transactions on
Signal Processing, vol. 55, no. 8, 2007, pp. 4284-4299.
[34] B. Chen, G. Hu, W. Zhang, L. Yu, Distributed mixed H2 /H∞ fusion
estimation with limited communication capacity, IEEE Transactions
on Automatic Control, vol. 61, no. 3, 2016, pp. 805-810.
[35] A. Chiuso, L. Schenato, Information fusion strategies and performance
bounds in packet-drop networks, Automatica, vol. 47, no. 7, 2011, pp.
1304-1316.
[36] J. Ma, S. Sun, Centralized fusion estimators for multisensor systems
with random sensor delays, multiple packet dropouts and uncertain
observations, IEEE Sensors Journal, vol. 13, 2013, pp. 1228-1235.
[37] Y. Xia, J. Shang, J. Chen, G.P. Liu, Networked data fusion with packet
losses and variable delays, IEEE Transactions on Systems, Man and
Cybernetics–Part B: Cybernetics, vol. 39, no. 5, 2009, pp. 1107-1120.
[38] B. Chen, W.A. Zhang, L. Yu, Distributed fusion estimation with missing measurements, random transmission delays and packet dropouts,
IEEE Transactions on Automatic Control, vol. 59, no. 7, 2014, pp.
1961-1967.
[39] R. Caballero-Aguila, A. Hermoso-Carazo, J. Linares-Perez, Fusion
estimation using measured outputs with random parameter matrices
subject to random delays and packet dropouts, Siganl Processing, vol.
127, 2016, pp. 12-23.
[40] Z. Xing, Y. Xia, Distributed federated kalman filter fusion over multisensor unreliable networked systems, IEEE Transactions on Circuits
and Systems I: Regular Papers, vol. 63, no. 10, 2016, pp. 1714-1725.
[41] Y. S. Lee, W. H. Kwon, Delay-dependent robust stabilization of
uncertain discrete-time state-delayed systems, Proceedings of the 15th
IFAC World Congress, Barcelona, Spain, 2002, pp. 610-615.
[42] S.P. Boyd, L.E. Ghaoui, E. Feron, V. Balakrishan, Linear Matrix
Inequalities in System and Control Theory, SIAM, Philadelphia, PA,
1994.
[43] IEEE PES Distribution System Analysis Subcommittee, Distributed
Test Feeders: 4-Bus Test Feeder Cases, [Online], Available:
http://ewh.ieee.org/soc/pes/dsacom/testfeeders/index.html

