1

Fault Estimation Filter Design with Guaranteed
Stability Using Markov Parameters

arXiv:1708.09034v1 [] 29 Aug 2017

Yiming Wan, Tamás Keviczky, and Michel Verhaegen

Abstract—For additive actuator and sensor faults, we propose
a systematic method to design a state-space fault estimation
filter directly from Markov parameters identified from faultfree data. We address this problem by parameterizing a systeminversion-based fault estimation filter with the identified Markov
parameters. Even without building an explicit state-space plant
model, our novel approach still allows the filter gain design
for stabilization and suboptimal H2 performance. This design
freedom cannot be achieved by other existing data-driven fault
estimation filter designs so far. Another benefit of our proposed
design is the convenience of determining the state order: a higher
state order of the filter leads to better estimation performance,
at the cost of heavier computational burden. In contrast, order
determination is cumbersome when using an identified statespace plant model for the filter design, because of the complicated
propagation of the model mismatch into the fault estimation
errors. Simulations using an unstable aircraft system illustrate
the effectiveness of the proposed new method.
Index Terms—Data-driven method, fault estimation, system
inversion, Markov parameters.

I. I NTRODUCTION
BSERVER-based fault diagnosis techniques have been
well established during the past two decades [1]. However, an explicit and accurate system model is often unknown
in practice. In such situations, a conventional approach follows
two steps: first identifying the state-space plant model from
system input/output (I/O) data, and then designing observers
for fault diagnosis [2]–[5]. Different from the conventional
two-step approach, the data-driven approach to fault diagnosis
observer design has been investigated recently for additive
sensor or actuator faults, without explicitly identifying a statespace plant model [6]. As an alternative to multivariate statistical process monitoring such as principle component analysis
(PCA) and partial least squares [7], [8], the data-driven fault
diagnosis observer design offers a more powerful tool for
highly dynamic systems, and allows developing systematic
methods to address the same fault diagnosis performance
criteria as the existing model-based approaches [6].
Despite the recent progress in direct data-driven observer
or filter design [9], research about dealing with completely
unknown disturbances or faults in data-driven filter design
has just started. Existing data-driven fault diagnosis observer

O

This work has received funding from the European Unions Seventh Framework Programme (FP7-RECONFIGURE/2007-2013) under grant agreement
no. 314544.
Yiming Wan is with Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139. Email: ywan@mit.edu
Tamás Keviczky and Michel Verhaegen are with Delft Center for Systems
and Control, Delft University of Technology, 2628CD, The Netherlands.
Emails: t.keviczky, m.verhaegen@tudelft.nl

design methods construct observers with either the parity vector [10] or the predictor Markov parameters (MPs) [11], [12]
that can be identified from data. Data-driven fault estimation
is much more involved than fault detection and isolation,
because it is inherently related to inverting the underlying
system whose model is unavailable. A non-recursive fault
estimator was proposed in [13] and [14] by minimizing the
squared reconstructed prediction error in the residual subspace
of a latent variable model. Recently, a receding horizon leastsquares fault estimator has been proposed in [15] by using
identified predictor MPs, which enables robust design that
compensates for identification errors of the MPs. Even fewer
studies have been reported on data-driven design of recursive
fault estimation observers or filters. Ding et al. first constructed
a diagnostic observer realized with the identified parity vector,
and then estimated faults as augmented state variables, see
Chapter 10 of [6]. This augmented observer scheme, however,
imposed certain limitations on how fault signals vary with
time, thus introduced bias in fault estimates. In contrast,
without any assumptions on the dynamics of fault signals,
Dong et al. first constructed a non-recursive fault estimation
filter (FEF) in the form of finite impulse response (FIR) from
the identified MPs, then used its state-space realization as a
recursive FEF [16].
As opposed to the model-based design, it is nontrivial to
design a stable FEF directly from data without identifying
an explicit state-space model. It is well known in modelbased design that the existence of a stable inversion-based FEF
is ensured when the fault subsystem has no unstable zeros
[17]–[20]. This property, however, cannot be guaranteed in
current data-driven FEFs. For example, even under the above
condition, 1) the parity vector based fault estimation observer
in Chapter 10 of [6] needs the augmented fault state with
assumed dynamics, which unnecessarily introduces estimation
bias; and 2) the MP-based FEF in [16] might still be unstable.
This paper focuses on data-driven design of FEF with stability guarantee, for additive actuator and sensor faults whose
fault subsystem has arbitrary relative degrees. This problem
is challenging, because it requires inverting the underlying
plant dynamics without building an explicit state-space plant
model. In order to pave the way for the data-driven design, an
FEF is first constructed given the plant model in the predictor
representation, which is structured into a residual generator
and the inverse of the residual dynamics. Such a structured
system-inversion-based FEF (SI-FEF) allows us to establish
the link between the MPs of the SI-FEF and the predictor
MPs. By exploiting this link, our data-driven design method
first computes the MPs of the SI-FEF with the predictor

2

MPs identified from data, and then constructs a state-space
realization of the SI-FEF from its MPs. Even without building
an explicit state-space plant model, our data-driven design still
allows the design freedom of the filter gain for stabilization
and suboptimal H2 performance, which is missing in other
existing data-driven designs.
Another benefit of our proposed design is related to the convenience of dertermining the state order: a higher state order of
the filter leads to better estimation performance, at the cost of
heavier computational burden. However, order determination is
cumbersome when using an identified state-space plant model
for the filter design, because of the complicated propagation
of the model mismatch into the fault estimation errors.
This paper is organized as follows. In Section II, we describe
the system and formulate the data-driven FEF design problem.
In Section III, a SI-FEF is constructed given the plant model
in the system predictor representation. Then the link between
the MPs of the SI-FEF and the predictor MPs is established in
Section IV. Our proposed data-driven design is developed in
Section V. The advantages of this new method are illustrated
via a numerical simulation example of an unstable aircraft
system in Section VI. Finally, we give some concluding
remarks in Section VII.
II. P RELIMINARIES AND PROBLEM FORMULATION
A. Notations

x(k + 1) = Φx(k) + B̃u(k) + Ẽf (k) + Ky(k),
y(k) = Cx(k) + Du(k) + Gf (k) + e(k),

(5a)
(5b)

where x(k) ∈ Rn and e(k) ∈ Rny are the predictor states
and the innovation signal, respectively. K is the steady-state
Kalman gain, Φ = A−KC, B̃ = B −KD, and Ẽ = E −KG.
No assumption is made about how the fault signals f (k) evolve
with time.
Define the MPs of the predictor representation (5) as


D
i=0
0
i=0
y
,
,
H
=
Hiu =
i
CΦi−1 K i > 0
CΦi−1 B̃ i > 0

G
i=0
Hif =
.
i−1
CΦ Ẽ i > 0
(6)
For the additive fault in the jth actuator or sensor, we may
construct the predictor MPs {Hif } from {Hiu } and {Hiy } as
below according to (5) and (6):
j th actuator fault:

For the state-space model (A, B, C, D), define Markov
parameters as H0 = D and Hi = CAi−1 B for i > 0. {Hi }
represents the sequence of Markov parameters. Let Os and Ts
denote the extended observability matrix with s block elements
and the lower triangular block-Toeplitz matrix with s block
columns and rows, respectively, i.e.,

 H
0
... 0
0
 C 
. . .. 
CA
 H
. . 
H0
1
Os (A, C) =  ..  , Ts ({Hi }) = 
,
 .
.
.
.s−1
..
..
.. 0
CA
Hs−1 Hs−2 ··· H0

(1)
D

0

... 0

CB

D

. . ..
. .
.
..
. 0



or Ts (A, B, C, D) = 

..
.

CA

s−2

..
.

B CA

s−3





y T (k − l + 1)

···

y T (k)

j

th

i ≥ 0;

(7a)

sensor fault:

E = 0, G = I [j] , Hif =



I [j]
−(Hiy )[j]

i=0
, (7b)
1<i≤p

where X [j] denotes the jth column of the matrix X.
The relative degree of the fault subsystem (Φ, Ẽ, C, G)
can be determined from its MPs {Hif }, i.e., the smallest
nonnegative integer τ such that Hτf is nonzero. Note that τ = 0
for sensor faults and τ > 0 for actuator faults. We adopt the
following assumption:
Assumption 1. The τ th MP of the fault subsystem
(Φ, Ẽ, C, G) has full column rank, where τ is the relative
degree of the fault subsystem.

(3)

Assumption 1 assumes sufficient number of measured outputs (ny ≥ nf for Hτf ∈ Rny ×nf ) and no collinearity
among the fault directions to ensure the uniqueness of fault
reconstruction. This assumption is common in fault estimation
or input reconstruction literature, e.g., [11], [18], [19], [21].

B ··· D

T

E = B [j] , G = D[j] , Hif = (Hiu )[j]

(2)

Define
yk,l =

and latent faults, respectively. The process and measurement
noises w1 (k) and w2 (k) are zero-mean white Gaussian.
A, B, C, D, E, G are time-invariant matrices unavailable to the
data-driven design. Assume that the system description (4)
admits a Kalman filter for its fault-free subsystem, then this
system (4) can be equivalently represented by the following
Kalman predictor representation [5], [15]:

by stacking data vectors {y(i)} in a sliding window [k − l + 1, k]. diag(X1 , X2 , · · · , Xn ) denotes a blockdiagonal matrix. E represents the mathematical expectation.

C. Markov parameter identification
B. System description
Consider a linear discrete-time system governed by
ξ(k + 1) = Aξ(k) + Bu(k) + Ef (k) + w1 (k)
y(k) = Cξ(k) + Du(k) + Gf (k) + w2 (k)

(4)

where ξ(k) ∈ Rn , u(k) ∈ Rnu , y(k) ∈ Rny , and f (k) ∈ Rnf
represent the states, system inputs, output measurements,

When the accurate knowledge of the state-space description
(4) or (5) is unavailable, we may identify the predictor MPs
from data. It is well known that the predictor representation (5)
can be approximated by the following vector ARX (VARX)
model with arbitrary accuracy as the VARX order becomes
sufficiently high [22], [23]:
A(q −1 )y(k) = B(q −1 )u(k) + F(q −1 )f (k) + v(k)

(8)

3

where q −1 is the backward shift operator, A(q −1 ) = I −
p
p
p
P
P
P
Miy q −i , B(q −1 ) =
Miu q −i , F(q −1 ) =
Mif q −i ,

i=1

i=0

i=0

v(k) ∈ Rny represent the noise signal. Therefore, the coefficients of the high-order VARX approximation can be the
estimates of the predictor MPs, i.e.,

Mis 0 ≤ i ≤ p
s
Hi ≈
, for s represents u, y, f. (9)
0
i>p
Note that H0y = 0 is already known in (6). For more detailed
derivations, we refer to Section 2.2 of [23].
With the fault-free identification data, we can identify the
VARX coefficients {Miu } and {Miy } as the estimates of the
predictor MPs {Hiu } and {Hiy }, and then construct {Hif } for
the additive faults according to (7). The residual signal v(k) =
A(q −1 )y(k) − B(q −1 )u(k) generated from the identification
data approximates the innovation e(k) of the predictor (5), and
can be used to estimate the innovation covariance as

Σe = cov A(q −1 )y(k) − B(q −1 )u(k) .
In practice, data collected under faulty conditions may be
seldom available, or recorded without a reliable description of
the fault type [6]. Therefore, no faulty historical data is used
in our data-driven design.
Remark 1. In theory, an infinite-order VARX model is needed
to fully represent a Kalman predictor (5). Therefore, the identification of the predictor MPs actually requires identifying an
infinite-order VARX model. For this purpose, we adopt a finite
high-order VARX approximation. It should be noted that the
consistent estimation of the infinite-order VARX coefficients,
i.e., the predictor MPs, does not follow the conventional rules
in the case of identifying a finite-order VARX model [22], [24].
In this paper, an empirical approach is used to select the order
p for the high-order VARX model: first, we may identify a loworder VARX or VARMA model in order to roughly estimate the
predictor poles and the noise variance; then, we determine the
order p according to the location of the predictor poles, so
that the remaining predictor MPs can be well approximated
by zero compared to the noise level. Note that the VARX order
selection also involves a trade-off, i.e., selecting a higher order
leads to smaller bias but larger variance in the identified
predictor MPs.
D. Data-driven design of fault estimation filter
Given the predictor MPs {Hiu , Hiy , Hif } identified offline
from data as in Section II-C, the basic idea of a systeminversion-based fault estimator follows two steps:
(i) Residual generation using the online I/O data, i.e., r(k) =
A(q −1 )y(k) − B(q −1 )u(k). Then the residual dynamics
is r(k) = F(q −1 )f (k) + e(k) according to (8).
(ii) τ -delay fault estimation by processing the residual signal
with the τ -delay left inverse of F(q −1 ), i.e., fˆ(k − τ ) =
F inv (q −1 )r(k), with F inv (q −1 )F(q −1 ) = q −τ Inf .
Note that the estimation delay τ in the above step (ii) is due
to the fact that the residual signal r(k) contains only the
fault information up to the time instant k − τ according to
the definition of the relative degree τ . Finding a stable left

inverse system F inv (q −1 ) is a long-studied problem in the
literature [17]–[19], [21], [25]–[27]. The capability of placing
poles of the left inverse system is critical to the stability and
performance of the system-inversion-based fault estimation.
To achieve this capability, an explicit state-space plant model
is needed in most system inversion literature, e.g., [17], [21],
[25]–[27], Chapter 3 of [18], and the references therein. Such
a pole tuning or placement capability becomes non-trivial
if only the knowledge of an input-output plant model is
available, which is the case for the data-driven design problem
in this paper. This prevents the applicability of the data-driven
methods in many situations.
In this subsection, we will briefly review the existing
approaches for constructing the left inverse system F inv (q −1 )
from the predictor MPs {Hiu , Hiy , Hif }, and point out their
limitations that motivate our research.
One category of a fault estimator is in the form of an FIR
filter. Note that the FIR filter is actually an approximated left
inverse of F(q −1 ), and its construction avoids stability and
pole placement mentioned above. Assume the order of this FIR
filter to be L, then the residuals r(k−L+1), · · · , r(k−1), r(k)
are involved to produce the fault estimate at time instant k. By
stacking the residual signal over the time window [k−L+1, k],
we obtain the stacked residual vector
rk,L = Ψf fk,L+p + ek,L ,

(10)

according to the VARX model (8), where rk,L , fk,L+p , and
ek,L are defined as in (3), and
 Mf

Mf
··· M f M f
0
··· 0
p



Ψf = 


f
Mp+1

..
.

f
Mp+L−1

p−1

Mpf

···

..
.

..

f
Mp+L−2

1

0

M2f

M1f

.
. ..

···

f
ML

..
.

f
ML−1

M0f

..
.

f
ML−2

. 
. .. 
.

..
. 0

..

··· M0f

Since the coefficients Mif of the high-order VARX model approximate the predictor MPs {Hif }, the first p block-columns
of Ψf approximate the Hankel matrix of the fault subsystem
(Φ, Ẽ, C, G). Because of this link, the first p block-columns
of Ψf can be ill-conditioned or rank deficient for large L and
p, see Theorem 6.1 of [28]. This implies that the stacked fault
vector fk,L+p may not be uniquely or reliably reconstructed
from the residual vector rk,L . Similarly to this above reason,
the dynamic PCA based fault reconstruction in [13] might be
non-unique, and it estimates a linear combination of the true
faults, as can be seen from the equation (24) of [13]. Such
a fault estimate is obviously biased. Instead of reconstructing
the entire fault vector fk,L+p , recent research in [15], [16]
produced the τ -delay fault estimate f (k − τ ) by applying an
FIR filter on the residual signal. It was shown in [15] that
unbiased fault estimates can be achieved asymptotically as the
FIR filter order L goes to infinity, if the fault subsystem has
no unstable zeros.
The second category of a fault estimator is in the state-space
form for the benefit of efficient online recursive computation.
The conventional approach designs a state-space FEF based on
a state-space plant model identified from data. This approach
might lead to large fault estimation errors, because there is
no effective method in literature to suppress the complicated

4

nonlinear propagation of the state-space system identification
errors into the fault estimates. Another approach in the recent
data-driven design is to construct the FEF as the state-space
realization of the aforementioned FIR filter [16]. However,
such an obtained state-space FEF is not guaranteed to be
stable, and its poles cannot be tuned in the design.
The aim of this paper is to construct a state-space SIFEF with tunable stable poles by using the predictor MPs
{Hiu , Hiy , Hif } identified from data. As summarized in Figure 1, our proposed approach constructs the SI-FEF from a
residual generator, an open-loop left inverse of the residual
dynamics, and the feedback from the residual reconstruction
errors. This structure allows (i) establishing the link connecting
{Hiu , Hiy , Hif } and the MPs of the SI-FEF as in Figure 2, and
(ii) designing the feedback gain of the residual reconstruction
errors for stability and performance.
Note that the identification errors of the predictor MPs
affect the fault estimation performance. This issue has been
investigated recently in [15] for the robust data-driven design
of a receding horizon fault estimator. How to address the same
issue for a data-driven state-space FEF can be investigated only
after the stability is ensured. In this paper, we focus on the
stability guarantee, and leave the robustness issue for future
research. Because of this reason, the identification errors of
the predictor MPs are not explicitly considered in this paper,
and we will use the notations {Hiu }, {Hiy }, {Hif } for both
the true predictor MPs and their estimates.
Closed-loop left inverse
Residual
generator

Open-loop
left inverse

Feedback from residual
reconstruction error

Predictor
Markov
parameters

Markov parameters of
residual generator
open-loop left inverse
feedback from residual
reconstruction errors

x1 (k + 1) = Φx1 (k) + B̃u(k) + Ky(k)
y1 (k) = Cx1 (k) + Du(k), x1 (0) = x̂(0),

(11a)
(11b)

and
x2 (k + 1) = Φx2 (k) + Ẽf (k)

(12a)

r(k) = Cx2 (k) + Gf (k) + e(k),

(12b)

such that x(k) = x1 (k) + x2 (k) and y(k) = y1 (k) + r(k).
Starting from an initial guess of the predictor state x̂(0), the
subsystem (11) predicts the output without accounting for the
fault. As shown in Figure 1, (11) is used to generate a residual
signal r(k) = y(k) − y1 (k) from the I/O data. Then, the
subsystem (12) is the residual dynamics decoupled from the
I/O data. This will be used in the following to design a closedloop left inverse system as depicted in Figure 1.
Since the fault subsystem (Φ, Ẽ, C, G) has the relative
degree τ (see Assumption 1), the residual signal at the time
instant k +τ is needed to produce a fault estimate fˆ(k), which
introduces an estimation delay when τ > 0. Considering this
estimation delay, we construct the following equation r(k +τ )
by successively substituting (11a) and (12a) into (11b) and
(12b), respectively:
r(k + τ ) = y(k + τ ) − y1 (k + τ )
= −CΦτ x1 (k) − Bτu+1 uk+τ,τ +1
+ Bτy+1 yk+τ,τ +1

(13a)

= CΦτ x2 (k) + Hτf f (k) + e(k + τ )

(13b)

where uk+τ,τ +1 and yk+τ,τ +1 are defined in (3), Bτu+1 and
Bτy+1 are respectively defined as


Bτu+1 = Hτu Hτu−1 · · · H0u ,


Bτy+1 = −Hτy −Hτy−1 · · · −H1y I .
In (13b), we use the fact Hif = 0 for i < τ according to the
definition of relative degree τ .
From (13b), f (k) can be estimated as below by using x2 (k)
and a left inverse matrix Π of Hτf :

Fig. 1. Our proposed fault estimation filter scheme

Fault-free
I/O data

Firstly, we decompose the predictor (5) into two subsystems:

Fault estimation Filter
state-space
realization
filter gain design

Fig. 2. Basic idea of our proposed data-driven design

III. S YSTEM - INVERSION - BASED FAULT ESTIMATION
FILTER USING THE PREDICTOR REPRESENTATION

As the foundation for our data-driven design, we construct
an SI-FEF in this section by exploiting the accurate knowledge
of the predictor representation (5). Note that all system matrices in (5) are unknown to our data-driven design, but used
here for establishing the link between the predictor MPs and
the MPs of the SI-FEF in Section IV.

fˆ(k) = Π [r(k + τ ) − CΦτ x2 (k)] , ΠHτf = I.

(15)

The left inverse matrix Π is a design parameter, whose
existence is ensured by Assumption 1. Since the state x2 (k)
is actually unknown, we construct the following left inverse
of the residual dynamics (12) and (13) in the state-space form
which jointly estimates the state and the fault:
x̂2 (k + 1) = Φx̂2 (k) + Ẽ fˆ(k) + Kr r̃(k + τ )
fˆ(k) = Π [r(k + τ ) − CΦτ x̂2 (k)]
r̃(k + τ ) = r(k + τ ) − r̂(k + τ )
τ

= r(k + τ ) − CΦ x̂2 (k) −

(16a)
(16b)
(16c)

Hτf fˆ(k).

(16d)

By replacing the state x2 and the fault f with their estimates
x̂2 and fˆ, r̂(k + τ ) = CΦτ x̂2 (k) + Hτf fˆ(k) in (16c) and
(16d) follows (13b) to reconstruct the residual signal from
the state and fault estimates. Then r̃(k + τ ) = r(k + τ ) −

5

r̂(k + τ ) is the residual reconstruction error. (16a) is a copy
of the residual dynamics (12a) with a feedback term Kr r̃(k +
τ ) from the residual reconstruction error r̃(k + τ ). Similarly,
(16b) constructs the fault estimate fˆ(k) by following (15). By
substituting (16b) into (16a) and (16d), respectively, the left
inverse (16) can be equivalently rewritten as
x̂2 (k + 1) = Φ1 x̂2 (k) + B1 r(k + τ ) + Kr r̃(k + τ ) (17a)
fˆ(k) = C1 x̂2 (k) + D1 r(k + τ )
(17b)
r̃(k + τ ) = −C2 x̂2 (k) + D2 r(k + τ )

(17c)

with
Φ1 = Φ − ẼΠCΦτ , B1 = ẼΠ, C1 = −ΠCΦτ ,
D1 = Π, C2 = (I −

Hτf Π)CΦτ ,

D2 = I −

Hτf Π.

(18)
(19)

With Kr = 0, (Φ1 , B1 , C1 , D1 ) in the above left inverse
system is referred to as an open-loop left inverse. With the
feedback gain Kr , the residual reconstruction error r̃(k + τ )
is used as a feedback signal to stabilize the above left inverse.
Such a structured form of the closed-loop inverse (17), i.e.,
the combination of the open-loop left inverse and the feedback
from the residual reconstruction errors r̃(k + τ ), enables our
data-driven design in Sections IV and V.
By cascading the residual generator (11) and the left inverse
(17), we obtain the SI-FEF as below:
x̂(k + 1) = Φ1 x̂(k) + Bf uk+τ,τ +1 + Kf yk+τ,τ +1
+ Kr r̃(k + τ )
fˆ(k) = C1 x̂(k) + Df,1 uk+τ,τ +1 + Gf,1 yk+τ,τ +1 ,
r̃(k + τ ) = −C2 x̂(k) − Df,2 uk+τ,τ +1 + Gf,2 yk+τ,τ +1 .
(20)
Note that x̂(k) = x1 (k)+ x̂2 (k) is an estimate of the predictor
state x(k), because x̂2 (k) is the estimate of x2 (k) and the
predictor state is decomposed as x(k) = x1 (k) + x2 (k). In
the above SI-FEF, Φ1 , B1 , C1 , D1 , C2 and D2 are defined in
(18) and (19), respectively, and




B̃τ = B̃ 0nx ×τ nu , Kτ = K 0nx ×τ ny ,
(21)
Bf
Df,1
Df,2

=
=
=

B̃τ − B1 Bτu+1 ,
−D1 Bτu+1 ,
D2 Bτu+1 ,

Kf
Gf,1
Gf,2

B1 Bτy+1 ,

= Kτ +
= D1 Bτy+1 ,
= D2 Bτy+1 .

Next, the error dynamics of the SI-FEF (20) is analyzed by
defining the state estimation error x̃(k) = x(k) − x̂(k) and the
fault estimation error f˜(k) = f (k) − fˆ(k):
x̃(k + 1) = (Φ1 − Kr C2 ) x̃(k) − (B1 + Kr D2 ) e(k + τ )
f˜(k) = C1 x̃(k) − D1 e(k + τ ).
(22)
Therefore, if the pair (Φ1 , C2 ) is observable or detectable,
there exists a stabilizing gain Kr in (22), such that starting
from any arbitrary initial estimate x̂(0), unbiasedness of the
estimates x̂(k) and fˆ(k) can be achieved
 asymptotically, i.e.,
˜
lim E (x̃(k)) = 0 and lim E f (k) = 0.
k→∞

k→∞

Theorem 1. (Φ1 , C2 ) is observable if the fault subsystem
(Φ, Ẽ, CΦτ , Hτf ) has no invariant zeros; (Φ1 , C2 ) is detectable if all invariant zeros of (Φ, Ẽ, CΦτ , Hτf ) are stable.

The proof is given in the Appendix. Theorem 1 shows
how the observability or detectability of the pair (Φ1 , C2 )
is determined by the invariant zeros of the underlying fault
subsystem. Thus it provides a sufficient condition for the
existence of a stabilizing filter gain for the SI-FEF (20).
The SI-FEF (20) produces both the state estimate x̂(k)
and the fault estimate fˆ(k). However, it is different from the
simultaneous state and input estimation filter proposed by [29],
[30] in two aspects: (i) the condition in Theorem 1 that ensures
stabilization and asymptotic unbiasedness was not provided;
(ii) only the special cases of τ = 0 and τ = 1 were discussed
in [29], [30].
How to design Π in (15) and Kr in (22) for our data-driven
design problem will be discussed in Section V-A.
IV. M ARKOV PARAMETERS OF SYSTEM - INVERSION - BASED
FAULT ESTIMATION FILTER

As illustrated in Figure 2, after the MPs of the SI-FEF
(20) are computed, the state-space realization of the SI-FEF
can be constructed. In this section, we establish the link for
computing MPs of the SI-FEF (20) from the predictor MPs
{Hiu , Hiy , Hif }.
As the first step towards the above goal, we rewrite the
residual generator (11), the left inverse system (17), and the
SI-FEF (20) into extended forms over a time window. With
k0 = k − L + 1, we define

T
T
,
(23)
z̄k,L = zT
k0 +τ,τ +1 · · · zk+τ,τ +1

T
T
by stacking zk+τ,τ +1 = uTk+τ,τ +1 yk+τ,τ
over the
+1
time window [k0 , k]. According to (11a), (12a), and (13), the
stacked residual vector rk+τ,L over the time window [k0 , k]
can be written in the extended form
rk+τ,L = OL (Φ, −CΦτ ) x1 (k0 ) + TLz z̄k,L
τ

= OL (Φ, CΦ ) x2 (k0 ) +

TLf fk,L

with B̃τ and Kτ defined in (21),



TLz = TL Φ, B̃τ Kτ , −CΦτ , −Bτu+1


TLf = TL Φ, Ẽ, CΦτ , Hτf .

(24a)

+ ek+τ,L (24b)

Bτy+1



, (25)
(26)

Since the residual generator (11) has the initial state
x1 (k0 ) = x̂(k0 ), the closed-loop left inverse (17) then has the
initial state x̂2 (k0 ) = 0 according to x̂(k) = x1 (k) + x̂2 (k)
in (20). Hence, the closed-loop left inverse (17) can be
transformed into the following extended form over the time
window [k0 , k] to produce the stacked fault estimates f̂k,L :
f̂k,L = KL rk+τ,L = (GL + ML JL ) rk+τ,L ,

(27)

with
KL = TL (Φ1 − Kr C2 , B1 + Kr D2 , C1 , D1 ) ,

(28a)

GL = TL (Φ1 , B1 , C1 , D1 ) ,

(28b)

JL = I − TLf GL = TL (Φ1 , B1 , −C2 , D2 ) ,

(28c)

ML = TL (Φ1 − Kr C2 , Kr , C1 , 0) .

(28d)

6

The proof of KL = GL + ML JL in (27) is given in the
Appendix. Note that KL , GL , JL and ML are all blockToeplitz matrices, and can be explained as below:
(i) GL corresponds to the open-loop left inverse, i.e., (17)
with Kr = 0;
(ii) JL produces the residual reconstruction errors r̃(k + τ )
in (17) with Kr = 0;
(iii) ML corresponds to the feedback dynamics from the
residual reconstruction errors r̃(k + τ ) in the closed-loop
inverse (17).
By substituting the residual generator (24) into the extended
closed-loop inverse (27), the following extended form of the
SI-FEF (20) is obtained:
f̂k,L = OL (Φ1 − Kr C2 , C1 ) x1 (k0 ) + (RL + ML QL ) z̄k,L
(29a)
= OL (Φ1 − Kr C2 , −C1 ) x2 (k0 ) + fk,L + KL ek+τ,L
(29b)
with
RL = GL TLz = TL Φ1 , [Bf
QL = JL TLz = TL Φ1 , [Bf


Gf,1 ] ,
(30a)

Kf ], −C2 , [−Df,2 Gf,2 ] .
(30b)

Kf ], C1 , [Df,1

Similarly to GL and JL in (27), RL and QL correspond to two
subsystems of the SI-FEF (20) with Kr = 0, which produce
fˆ(k) and r̃(k + τ ) in the open loop, respectively. ML is the
same feedback dynamics as in (28d).
The extended form (29a) can be regarded as a batch
estimator which provides the estimate f̂k,L from the I/O data
z̄k,L and the initial state x1 (k0 ) = x̂(k0 ). Moreover, it can
be seen from (29b) that f̂k,L is a biased estimate of fk,L due
to the presence of unknown initial state x2 (k0 ). However, it
follows from the definition of OL (Φ1 − Kr C2 , −C1 ) in (1)
that


L−1
E fˆ(k) − f (k) = −C1 (Φ1 − Kr C2 )
x2 (k0 ),
where fˆ(k) and f (k) are the last nf entries of f̂k,L and fk,L ,
respectively. The above equation shows that fˆ(k), extracted
from f̂k,L in (29a), gives asymptotically unbiased fault estimation as L goes to infinity, if Φ1 − Kr C2 is stabilized given
the condition in Theorem 1.
In the above derivations, the block-Toeplitz matrices TLz ,
TLf , GL , JL , and QL are expressed with state-space matrices.
For the data-driven design, the next step is to construct their
corresponding MPs defined as


TLz = TL ({Hiz }) , TLf = TL {Hif } , GL = TL ({Gi }) ,
JL = TL ({Ji }) , RL = TL ({Ri }) , QL = TL ({Qi }) ,
(31)
from the predictor MPs {Hiu , Hiy , Hif }. To achieve this goal,
we first need to take a closer look at TLz , TLf and GL which
are needed in computing RL and QL . According to (25) and

(30)

Residual generator

Predictor
(32)

Open-loop left inverse

(33)

Fault estimation filter

(35)

(31)

(34)

(34)

Generator of residual
reconstruction error

Fig. 3. Link between predictor MPs and MPs of SI-FEF

(26), the MPs {Hiz } and {Hif } can be computed from the
predictor MPs {Hiu , Hiy , Hif } as below:



y
z
u


 H0 = −Bτu+1 Bτ +1 u


= −Hτ · · ·  −H0 −Hτy · · · −H1y I

τ +i−1
Hiz = −CΦ
B̃τ Kτ

 u


0
Hτy+i 0ny ×τ ny ,
H
=
−

n
y ×τ nu
τ
+i


for 1 ≤ i ≤ L − 1,
(32)

f
f
H0 = Hτ ,
(33)
Hif = CΦτ +i−1 Ẽ = Hτf+i , for 1 ≤ i ≤ L − 1.
As pointed out in the explanations below (27) and (28), GL
is a left inverse matrix with block-Toeplitz structure for TLf .
Such a left inverse matrix is non-unique, but can be computed
from the MPs {Hif }. With ΠH0f = ΠHτf = I according to
(15) and (33), one possible solution of GL is given below:

 G0 = Π,
i
P
(34)
Gi−j Hjf G0 , for 1 ≤ i ≤ L − 1.
 Gi = −
j=1

which ensures GL TLf = I. Then, according to (30), the MPs
of RL can be computed as the convolution of {Gi } in (34)
and {Hiz } in (32):
Ri =

i
X

Gi−j Hjz , for 0 ≤ i ≤ L − 1.

(35)

j=0

Similarly, the MPs {Ji , Qi } of JL in (28c) and QL in (30)
can be computed as
(
J0 = I − H0f G0 ,
Pi
(36)
f
Ji = − j=0 Hi−j
Gj , for 1 ≤ i ≤ L − 1,
Qi =

i
X

Ji−j Hjz , for 0 ≤ i ≤ L − 1.

(37)

j=0

Equations (32)-(37) reveal the link from the predictor MPs
to the MPs of the SI-FEF (20), as summarized in Figure 3.
V. FAULT ESTIMATION FILTER DESIGN USING M ARKOV
PARAMETERS

By exploiting the link between the predictor MPs and the
SI-FEF MPs, as analyzed in Section IV, the proposed MP
based data-driven design is given as below.
Algorithm A. Data-driven design of fault estimation filter

7

(i) Identify the predictor MPs {Hiu } and {Hiy } using VARX
modelling with the historical or experimental fault-free
I/O data.
(ii) Compute MPs of SI-FEF (20).
Construct the MPs {Hif }, {Hiz }, and {Hif } according
to (7), (32), and (33), respectively. Select
one left inverse
−1
(Hτf )T . Then
matrix Π of Hτf , e.g., Π = (Hτf )T Hτf
compute {Gi }, {Ji }, {Ri }, and {Qi } by following (34)(37).
(iii) State-space realization of the SI-FEF (20) from the MPs
{Ri } and {Qi }.
According to (30) and (31), the MPs {Ri } and {Qi } cor-
respond to systems Φ1 , [Bf Kf ], C1 ,[Df,1 Gf,1 ]
and Φ1 , [Bf Kf ], −C2 , [−Df,2 Gf,2 ] , respectively.
Then it is straightforward to obtain




D̂f,1 Ĝf,1 = R0 , −D̂f,2 Ĝf,2 = Q0 .
Formulate two block-Hankel matrices HR and HQ as
 W1 W2 ··· Wm 
W2

HW =  ..
.

W3

···

Wm+1

..
.

..

..
.

.

 , W represents R or Q,

Wl Wl+1 ··· Wl+m−1

(38)
then compute their singular value decomposition (SVD),
i.e.,


T 

 ΣW
VW
0
⊥

HW = UW UW
⊥ T .
0
Σ⊥
VW
W
In this above equation, the nonsingular and diagonal
matrices ΣR and ΣQ consist of the n̂ largest singular
values of HR and HQ , respectively, where n̂ is the
selected order of the fault estimation filter (20). The order
n̂ can be chosen by examining the gap among the singular
values of HR and HQ , respectively, as in subspace
identification methods [31]. Let the rank-reduced blockHankel matrices ĤR and ĤQ be
T
ĤW = UW ΣW VW
, W represents R or Q.

(39)

For ĤR defined in (39), the estimated controllability and
observability matrices can be constructed as [31]
1

1

CˆR = ΣR2 VRT , ÔR = UR ΣR2 .

(40)

Then the state-space realization of ĤR are computed as
below:
[B̂f

K̂f ] = the first nu + ny columns of CˆR ,
Ĉ1 = the first nf rows of ÔR ,

−1
T
T
Φ̂1 = CˆR,2 CˆR,1
CˆR,1 CˆR,1
,

where CˆR,1 and CˆR,2 are the matrices consisting of the
first and, respectively, the last nu (m − 1) columns of
CˆR . According to (30), the state-space realizations of
the block-Hankel matrices ĤR and ĤQ have the same
controllability matrix, i.e., CˆR obtained in (40). Then the
observability matrix in the state-space realization of ĤQ
can be computed below by using ĤQ = ÔQ CˆR :

−1
T
T
ÔQ = ĤQ CˆR
CˆR CˆR
.
(41)

Finally, we have
− Ĉ2 = the first ny rows of ÔQ .

(42)

(iv) Design the filter gain Kr by following Algorithm B in
Section V-A; and construct the SI-FEF (20) with the
identified system matrices in Step (iii).
Remark 2. The VARX model order p in Step (i) is selected
according to Remark 1. In Step (ii), the length L of the SIFEF MPs needs to be sufficiently large to ensure satisfactory
fault estimation performance. This is due to the asymptotic
unbiasedness of the batch fault estimation (29) as L goes to
infinity, which is explained in Section IV. In Step (iii), we select
the size of the block-Hankel matrix in (38) to be l + m = L,
with l and m defined in (38). By doing so, all MPs {Ri , Qi }
(i = 1, 2, · · · , L) obtained in Step (ii) are used to construct
HR and HQ in (38).
A. Suboptimal design of filter gain
The joint design of both Π in (15) and the filter gain Kr
is extremely difficult, because all system matrices in the SIFEF (20) depend on Π. Alternatively, our proposed data-driven
design selects Π in Step (ii) of Algorithm A before designing
the steady-state filter gain Kr in Step (iv) of Algorithm A.
Therefore, designing the filter gain Kr given a predefined Π is
suboptimal compared to the joint design. Recall that in some
existing model-based unknown input estimation methods, Π
and Kr were jointly designed to achieve globally unbiased
minimum-variance estimation, e.g., in [29], without discussing
stability of the obtained filter therein.
Based on the fault estimation error dynamics (22), the H2
fault estimation problem can be formulated as
1

min kĈ1 (zI − Φ̂1 + Kr Ĉ2 )−1 (B̂1 + Kr D̂2 )Σe2 k22
Kr

(43)

to find the steady-state filter gain Kr . It is well known that
the solution Kr to the problem (43) does not depend on Ĉ1 ,
and is actually the steady-state Kalman filter gain, see Section
6.5 of [32] and Section 7.3 of [33]. In this above problem
formulation, Φ̂1 , Ĉ1 , and Ĉ2 are obtained in Algorithm A as
the estimates of Φ1 , C1 , and C2 , respectively, while estimating
B̂1 and D̂2 will be explained later in Step (i) of Algorithm B.
With these above estimated matrices, the solution to the
problem (43) is discussed as below. Note that in Step (i) of
Algorithm B, we have
D̂2 = J0 = I − Hτf Π

(44)

according to (33), (34), and (36), and we have ΠD̂2 = 0
since ΠHτf = I. Then it can be seen that D̂2 is row-rank
deficient, hence the solution to (43) is non-unique. To tackle
this problem, we follow [34] to restrict the filter gain Kr to
be in the form
Kr = K̄r α,
where α ∈ Rs×ny ensures rank(D̂2 ) = rank(αD̂2 ) = s. Then
the H2 optimization problem (43) becomes
1

min kĈ1 (zI − Φ̂1 + K̄r C̄2 )−1 (B̂1 + K̄r D̄2 )Σe2 k22
K̄r

(45)

8

with C̄2 = αĈ2 and D̄2 = αD̂2 . With a proper selection of α,
the sufficient and necessary condition given below in Theorem
2 guarantees that the solution to (45), i.e., [32], [33]


(46)
K̄r = Φ̂1 P C̄2T + B̂1 Σe D̄2T Ξ−1
e

then we select α = U2T so that αD̂2 = α(I −Hτf Π) = U2T
is full row rank according to (44).
(iii) With C̄2 = αĈ2 and D̄2 = αD̂2 , compute K̄r in (46) by
solving the ARE (47). Then the filter gain is Kr = K̄r α.

stabilizes the SI-FEF (20), where P is the stabilizing solution
to the algebraic Riccati equation (ARE)

B. Comparisons and discussions

P = Φ̂1 P Φ̂T1 + B̂1 Σe B̂1T
(47a)



T
− Φ̂1 P C̄2T + B̂1 Σe D̄2T Ξ−1
Φ̂1 P C̄2T + B̂1 Σe D̄2T ,
e
Ξe = C̄2 P C̄2T + D̄2 Σe D̄2T .

(47b)

Lemma 1. The selected
  α in Step (ii) of Algorithm B ensures
α
that (i) the matrix
is nonsingular; and (ii) ΠĈ2 = 0.
Π
The proof is given in the Appendix.
Despite the identification errors in Ĉ2 compared to the true
C2 defined in (19), the condition (ii) of Lemma 1 still holds,
just like the fact that ΠC2 = 0 holds for the true C2 . This
condition will be used in the proof of Theorem 2 below.
Theorem 2. With Assumption 1 and the selection of α in Step
(ii) of Algorithm B, the ARE (47) has a unique stabilizing
solution P if and only if


Φ̂1 − λI
= n, for |λ| ≥ 1,
(48a)
rank
Ĉ2


Φ̂ − ejω I B̂1
rank 1
= n + ny , for ω ∈ [0, 2π]. (48b)
Ĉ2
D̂2
The proof is given in the Appendix.
Theorem 2 shows that the existence of a unique stabilizing
solution to the H2 estimation problem (45) depends on the
system matrices of (Φ̂1 , B̂1 , Ĉ2 , D̂2 ), despite errors in these
identified matrices. Such a suboptimal design and its stability
guarantee are not provided in recently proposed data-driven
fault estimation observers in [16] and Chapter 10 of [6].
The design method for the filter gain Kr is now summarized
in Algorithm B.
Algorithm B. Suboptimal design of filter gain
(i) Identify B1 and D2 using the MPs {Ji } identified in the
Step (ii) of Algorithm A.
From (28c) and (31), we can see that {Ji } are the
MPs of the system (Φ1 , B1 , −C2 , D2 ). It is then easy
to obtain D̂2 = J0 . Formulate the block-Hankel matrix
HJ with the MPs {Ji } by using the definition (38). With
the selected filter order n̂, we compute the rank-reduced
matrix ĤJ by following procedures similar to Step (iii) of
Algorithm A. Since the observability matrix of the statespace realization of ĤJ is the same as that of ĤQ , i.e.,
ÔQ in (41), we can compute the controllability matrix
CˆJ of ĤJ as below by using ĤJ = ÔQ CˆJ :
T
T
CˆJ = (ÔQ
ÔQ )−1 ÔQ
ĤJ .

Finally, we obtain B̂1 as the first ny columns of CˆJ .
(ii) Let the SVD of Hτf be
 

 SH
Hτf = U1 U2
V T,
0

The data-driven FIR fault estimator have been reviewed in
Section II-D. Next, we focus on further comparisons with other
existing state-space FEF designs from data.
The data-driven method of [16] considered only the openloop left inverse GL (28b) corresponding to (17) with Kr = 0.
Hence, it has no stability guarantees. In contrast, our datadriven design is based on the closed-loop left inverse (17)
and its extended form (29) which ensure the stability and
asymptotical unbiasedness under the condition specified in
Theorem 1 and explained in Section IV.
For the above reason, the data-driven filter in [16] cannot
be applied to sensor faults of an unstable open-loop plant.
It is worth noting that this difficulty cannot be solved by
simply applying the same method to the stabilized closedloop system. The reason is that the sensor faults affect not
only the output equations but also the closed-loop dynamics,
hence (7) is no longer valid for the MPs {Hif } of the closedloop fault subsystem. Section 2 of [35] provides more detailed
analysis about the effect of sensor fault propagation in closedloop dynamics. To circumvent this difficulty, Section V-B of
[16] proposed to use a special control law such that the sensor
faults did not affect the closed-loop dynamics, which is not
always possible in practice.
With the parity vector identified from data, Chapter 10 of
[6] constructed a diagnostic observer and estimated faults as
augmented state variables. This augmented observer scheme,
however, imposed certain limitations on how fault signals vary
with time, thus introduced bias in fault estimates. In contrast,
our proposed data-driven design needs no assumptions about
the time-varying fault signals. Moreover, for systems with
multiple outputs, the parity vector based approach in Chapter
10 of [6] becomes much complicated, while our proposed
method remains the same.
The suboptimality of our data-driven design is mainly
due to three reasons: firstly, the identification errors of the
predictor MPs in Step (i) of Algorithm A is neglected, as
pointed out in the last paragraph of Section II; secondly, the
unmodelled dynamics of the state-space realization in Step
(iii) of Algorithm A is not explicitly considered; and thirdly,
the matrix Π is not jointly designed with the filter gain, as
explained in the first paragraph of Section V-A.
Further comments are in order for the second reason above
in our suboptimal design. The unmodelled dynamics in our
proposed approach is the result of approximating the batch
fault estimator (29a) with a state-space filter. The higher state
order of the designed filter leads to better approximation, thus
giving better fault estimation performance. Because of this, the
obtained state-space filter cannot achieve better performance
than the batch estimator (29a) given a fixed horizon length L.
Therefore, considering also heavier online computational load
due to a higher state order, the order determination of our

9

designed filter is a simple trade-off between the fault estimation performance and the computational load. In contrast, the
order determination is cumbersome for the conventional twostep design and the parity vector based design in Chapter 10
of [6]. In these two approaches, we need to select the order of
a state-space plant model or a parity vector that represents a
residual subspace. The cumbersome issue is that the model
mismatch is introduced in the very first step of these two
approaches, and propagates through all the sequential steps.
Due to this complicated error propagation, there are no clear
guidelines for selecting the state order for the fault estimation
performance.
The two sources of uncertainties described in the first two
reasons of the suboptimality of our proposed approach are
common in most existing data-driven design methods. How to
explicitly quantify and deal with their effects remains an open
problem. One solution to the identification errors of MPs has
been investigated in [15] by using a nonrecursive receding
horizon estimator. However, it is still a challenge to address
both sources of uncertainties mentioned above in the datadriven design of a recursive FEF filter.
VI. S IMULATION STUDIES
Consider the linearized continuous-time vertical takeoff and
landing aircraft model used in [16]. The model has four states,
namely horizontal velocity, vertical velocity, pitch rate, and
pitch angle. The two inputs are collective pitch control and
longitudinal cyclic pitch control, both of which are fed through
the second-order actuator
21.3501s + 162.3867
.
s2 + 17.9994s + 162.3867
The sampling interval is 0.5 second. The process and measurement noises are zero mean white, with covariances of 10−4 I4
and 0.0016I2 respectively. Since the open-loop plant is unstable, an empirical stabilizing output feedback controller is
used, i.e.,


0 0 −0.5
0
u(k) = −
y(k) + η(k),
(49)
0 0 −0.1 −0.1
where η(k) is the reference signal. All the parameters of the
plant and the controller are the same as those in [16]. The
plant model is unknown to our data-driven design problem.
In the identification experiment, the reference signal η(k)
is zero-mean white noise with the covariance of diag (1, 1),
which ensures persistent excitation. We collect N = 100000
data samples, and select the VARX model order to be p = 12
by following Remark 1.
The simulated fault scenarios include: 1) faults in the two
actuators; 2) faults in the first two sensors. In both scenarios,
the reference signal η(k) in the control law (49) is set to be

T
2 2 , and the fault signals are
( 
T
0 0
,
0 ≤ k ≤ 500,
f (k) =

T
1 sin (0.01πk)
, k > 500.
We will compare the following methods for data-driven FEF
design:

Alg0: the SI-FEF (20) using the accurate predictor model
(5);
• Alg1: the SI-FEF (20) using the state-space model of the
predictor (5) identified from data;
• Alg2: the data-driven method proposed in [16];
• Alg3: our proposed new method in Section V.
The first algorithm is model-based, and the other three methods are all data-driven. For Alg2 and Alg3, the length of the
time window to construct the data-driven FEF is L = 90, and
the number of block rows and columns of the block-Hankel
matrix HW in (38) is l = m = 45, according to the guidelines
in Remark 2.
First, all state orders in the three data-driven designs are
set to be 8, i.e., the true state order of the underlying system.
By doing so, we focus on the stability of the obtained filters.
The estimated fault signals are shown in Figure 4(a). With
the accurate plant model, we can see that the fault subsystem
has stable invariant zeros in the actuator fault scenario, and
no invariant zeros in the sensor fault scenario. Therefore,
the model-based approach Alg0 gives stable FEFs in both
faulty scenarios according to Theorems 1 and 2. However,
Alg2 results in an unstable FEF in the sensor fault scenario
due to the reason explained in Section V-B, thus it is not
plotted for the sensor faults. In contrast to Alg2, Alg1 and
Alg3 are based on the closed-loop left inverse (17), and their
stability is guaranteed under the conditions in Theorems 1
and 2. Because of the model identification errors multiplying
with the online I/O data, the three data-driven designs, i.e.,
Alg1, Alg2, and Alg3, all give larger estimation errors than
Alg0. In comparison, Alg2 and Alg3 suffer less from model
identification errors than Alg1.
Next, we examine the state order selection in the datadriven designs. For Alg1, the state order is determined in
the phase of state-space plant model identification, which
results in unavoidable model mismatch. Due to complicated
error propagation from model mismatch, the fault estimation
performance may drastically change with different state orders,
as illustrated in Figures 4(b) and 5 when different state orders
from 8 to 20 are used. Moreover, Figure 5 shows even the
selection of the true plant state order 8 may not necessarily
give smaller fault estimation errors. In contrast, the state-space
filters in Alg2 and Alg3 are both approximations to batch fault
estimators. Then the selection of a higher state order in Alg2
and Alg3 leads to better approximation, thus gives better fault
estimation performance in Figures 4(b) and 5. This comparison
shows that it is much easier to determine the state order of
our proposed design: choose a high state order as long as the
online computational load is allowed. More detailed reasons
are explained in the second-to-last paragraph of Section V-B.
Besides the state order selection discussed above, the VARX
order p is also critical, since it is related to the accuracy
of the identified predictor MPs. For different VARX order
selections p = 10, 12, 14, we implement 100 Monte Carlo
runs for each p, with the state orders of Alg1, Alg2, and
Alg3 set to 8, 14, and 14, respectively. The results are shown
as boxplots in Figure 6. With the VARX order increasing,
the identified predictor MPs have smaller biases but larger
variances. This results in slightly worse performance for Alg2
•

10

Alg1
Alg2
Alg3

1.2

root mean square error

1

0.8

0.6

10

5

1.5
1

0.2

f1

1
f1

True
fault
Alg0
Alg1
Alg2
Alg3

Sensor faults
15

0.4

Sensor faults

Actuator faults

Actuator faults
1.4

root mean square error

and Alg3 when a larger VARX order is used. In comparison,
our proposed Alg3 gives the smallest averaged root mean
square error of fault estimates, Alg2 has a larger averaged
root mean square error but the smallest variance, and Alg1
gives the worst performance. In the above Monte Carlo runs
of the actuator fault scenario, our proposed Alg3 does not
consistently perform better than Alg2 for all selections of
the VARX order p. However, it should be noted that in all
the Monte Carlo runs of the sensor fault scenario, Alg3 can
stabilize the designed filter whereas Alg2 fails to do so.

0.5

0.5

0

0

0
500

600

700

500

600

12 14 16
state order

18

0

20

8

10

12 14 16
state order

18

20

Fig. 5. Root mean square error of fault estimates when selecting different
state orders. (The results of Alg2 are not plotted for sensor faults because it
gives unstable filters.)

1

0.5

10

700

1.5
1

8

f2

f2

0

0

−1

Alg1

−0.5

Alg2

Alg3

−2

−1
500

600
Time (k)

700

500

600
Time (k)

700

1

1

10

1

10

10

(a) State order: 8
Actuator faults

Sensor faults

1

0

0.8

0

−2
f1

0.6
f1

True
fault
Alg0
Alg1
Alg2
Alg3

0.4

0

10

10

−6

0.2
0

−8
500

600

700

500

600

700

1

−1

−1

10

−1

10

10

4

0.5
0

f2

f2

0

10

−4

2

−0.5
−1
−1.5

p=10 p=12 p=14

0
500

600
Time (k)

700

500

600
Time (k)

700

(b) State order: 18

p=10 p=12 p=14

p=10 p=12 p=14

Fig. 6. Boxplots of root mean square error of fault estimates in Monte Carlo
simulations of the actuator fault scenario: 100 Monte Carlo runs for each
different selection of the VARX order p; the state orders of Alg1, Alg2, and
Alg3 are set to 8, 14, and 14, respectively.

Fig. 4. Fault estimates given by different methods. (The result of Alg2 is not
plotted for sensor faults because it gives an unstable filter.)

VII. C ONCLUSIONS
A novel direct data-driven design method has been proposed
for FEFs by parameterizing the system-inversion-based fault
estimation filter with predictor Markov parameters. The proposed approach does not need to identify a state-space plant
model, but still allows the filter gain design for stabilization
and suboptimal H2 performance. This has not been achieved
by other existing data-driven fault estimation methods so far.
Moreover, the fault estimation performance can be improved
by simply increasing the state order of the designed filter,
at the cost of higher online computational load. A numerical
simulation example illustrates the advantages of our method

applied to actuator and sensor faults of an unstable aircraft
system. Future work will focus on the robustification of our
data-driven design against identification errors of Markov parameters, and non-Gaussian distributions in real applications.

A PPENDIX
P ROOF OF T HEOREM 1
In order to prove (Φ1 , C2 ) is detectable, we need to show
that (Φ1 , C2 ) has no unstable unobservable modes, i.e.,


Φ1 − λI
rank
= n for |λ| ≥ 1.
(50)
C2

11

By following (18) and (19), it can be derived that


 

I
0
Φ − λI Ẽ
Φ1 − λI Ẽ
. (51)
=
CΦτ
Hτf −ΠCΦτ I
C2
Hτf
With Assumption 1, if (Φ, Ẽ, CΦτ , Hτf ) has no unstable
invariant zeros, it follows that




Φ − λI Ẽ
Φ1 − λI Ẽ
= n + nf
= rank
rank
CΦτ
Hτf
C2
Hτf
for |λ| ≥ 1, which implies (50). The proof for observability
of (Φ1 , C2 ) is similar, thus is omitted.

Since KL , GL , ML , and JL are block-Toeplitz matrices
defined in (28), we can prove KL = GL + ML JL by proving
Ki = Gi +

Mi−j Jj

P ROOF OF T HEOREM 2
It is well known that the ARE (47) has a unique stabilizing
1
solution if and only if (Φ̂1 , C̄2 ) is detectable and (Fs , Qs2 ) is
controllable on the unit circle, see Appendix E of [36], where
Fs = Φ̂1 − B̂1 Σe D̄2T D̄2 Σe D̄2T

P ROOF OF KL = GL + ML JL

i
X

we have col(ÔQ ) = col(ĤQ ) ⊆ col(HQ ) for the rank-reduced
block-Hankel matrix ĤQ and the corresponding observability
matrix ÔQ , where col(X) represents the column space of a
matrix X. Then ÔQ satisfies Π̄ÔQ = 0. Consequently, Ĉ2
obtained in (42) ensures the condition (ii).

(52)

j=0

1
2

1
2

Qs = B̂1 Σe −

B̂1 Σe D̄2T

−1

C̄2 ,

−1
D̄2 Σe D̄2T

(54a)
1
2

D̄2 Σe .

(54b)

Using Lemma 1, the detectability of (Φ̂1 , C̄2 ) is given by





I 0 
Φ̂1 − λI
Φ̂1 − λI


rank
= rank 0 α
C̄2
Ĉ2
0 Π
(55)


Φ̂ − λI
= n, for |λ| ≥ 1,
= rank 1
Ĉ2

where Ki , Gi , Mi and Ji are the MPs that construct KL , GL ,
ML , and JL , respectively, as described in (31). According to
(28), these MPs are

D1
i=0
which proves (48a).
Ki =
,
1
C1 (Φ1 − Kr C2 )i−1 (B1 + Kr D2 ) i > 0
Let λ1 be an uncontrollable mode of (Fs , Qs2 ). This is


D1
i=0
D2
i=0
Gi =
, Ji =
, equivalent to the existence of a nonzero row vector ν such
i−1
C1 Φi−1
B
i
>
0
−C
Φ
B
1
2 1
1 i>0
that
1
h
i

1
ν Fs − λ1 I Qs2 = 0.
0
i=0
Mi =
.
C1 (Φ1 − Kr C2 )i−1 Kr i > 0
(53) From (54), the above equation can be rewritten as
#

"
By using (53), we can prove (52) as follows:
0
Φ̂1 − λ1 I B̂1 I
1
ν1
=0
Ki = C1 (Φ1 − Kr C2 )i−1 B1 + Mi J0
Ĉ2
D̂2 0 Σe2
h
= C1 (Φ1 − Kr C2 )i−2 (Φ1 B1 − Kr C2 B1 ) + Mi J0
−1 i
with ν1 = ν I −B̂1 Σe D̄2T D̄2 Σe D̄2T
α . Therefore, the
i−2
= C1 (Φ1 − Kr C2 ) Φ1 B1 + Mi−1 J1 + Mi J0
1
2
controllability of (Fs , Qs ) on the unite circle is equivalent to
2
X
i−3 2
(48b).
= C1 (Φ1 − Kr C2 ) Φ1 B1 +
Mi−j Jj
j=0

R EFERENCES

···
= C1 (Φ1 − Kr C2 )Φi−2
1 B1 +

i−2
X

Mi−j Jj

j=0

= Gi +

i−1
X

Mi−j Jj + M0 Ji .

j=0

P ROOF OF L EMMA 1
α
Proof of (i): Suppose that [ Π
] is singular, then there exist
ny −nf
nonzero vectors β1 ∈ R
and β2 ∈ Rnf such that β1T α +
β2T Π = 0. Under this condition, we have (β1T α+β2T Π)Hτf = 0.
The above two equations imply β1 = 0 and β2 = 0 because
αHτf = 0 (see Step (ii) of Algorithm B) and ΠHτf = I. The
contradiction with the nonzero assumption about β1 and β2
proves (i).
Proof of (ii): Using ΠHτf = I, (36) and (37), it is straightforward but tedious to verify that ΠJi = 0 and ΠQi = 0 for
i ≥ 0. Then we have Π̄HQ = 0, where HQ is defined in (38)
and Π̄ = diag(Π, Π, · · · , Π) with l diagonal blocks. Note that

[1] S. X. Ding, Model-Based Fault Diagnosis Techniques: Design Scheme,
Algorithms, and Tools, 2nd ed. London: Springer-Verlag, 2013.
[2] S. Simani, S. Fantuzzi, and R. Patton, Model-Based Fault Diagnosis in
Dynamic Systems Using Identification Techniques. London: SpringerVerlag, 2003.
[3] S. C. Patwardhan and S. L. Shah, “From data to diagnosis and control
using generalized orthonormal basis filters. Part I: development of state
observers,” Journal of Process Control, vol. 15, pp. 819–835, 2005.
[4] S. Manuja, S. Narasimhan, and S. C. Patwardhan, “Unknown input
modeling and robust fault diagnosis using black box observers,” Journal
of Process Control, vol. 19, pp. 25–37, 2009.
[5] R. Mustata, M. Verhaegen, H. Ohlsson, and F. Gustafsson, “Receding
horizon estimation of arbitrarily changing unknown inputs,” in Proceedings of the 19th IFAC World Congress, Cape Town, Sowth Africa, 2014,
pp. 5939–5944.
[6] S. X. Ding, Data-Driven Design of Fault Diagnosis and Fault-Tolerant
Control Systems. London: Springer-Verlag, 2014.
[7] S. J. Qin, “Survey on data-driven industrial process monitoring and
diagnosis,” Annual Reviews in Control, vol. 36, pp. 220–234, 2012.
[8] D. Gorinevsky, “Fault isolation in data driven multivariate process
monitoring,” IEEE Transactions on Control Systems Technology, vol. 23,
pp. 1840–1852, 2015.
[9] C. Novara, F. Ruiz, and M. Milanese, “Direct filtering: a new approach
to optimal filter design for nonlinear systems,” IEEE Transactions on
Automatic Control, vol. 58, pp. 86–99, 2013.

12

[10] S. X. Ding, P. Zhang, A. Naik, E. Ding, and B. Huang, “Subspace
method aided data-driven design of fault detection and isolation systems,” Journal of Process Control, vol. 19, pp. 1496–1510, 2009.
[11] J. Dong, M. Verhaegen, and F. Gustafsson, “Robust fault detection with
statistical uncertainty in identified parameters,” IEEE Transactions on
Signal Processing, vol. 60, pp. 5064–5076, 2012.
[12] ——, “Robust fault isolation with statistical uncertainty in identified
parameters,” IEEE Transactions on Signal Processing, vol. 60, pp. 5556–
5561, 2012.
[13] S. J. Qin and W. Li, “Detection and identification of faulty sensors in
dynamic processes,” AIChE Journal, vol. 47, pp. 1581–1593, 2001.
[14] C. Lee, S. W. Choi, and I. Lee, “Variable reconstruction and sensor
fault identification using canonical variate analysis,” Journal of Process
Control, vol. 16, pp. 747–761, 2006.
[15] Y. Wan, T. Keviczky, M. Verhaegen, and F. Gustaffson, “Data-driven
robust receding horizon fault estimation,” Automatica, vol. 71, pp. 210–
221, 2016.
[16] J. Dong and M. Verhaegen, “Identification of fault estimation filter
from I/O data for systems with stable inversion,” IEEE Transactions
on Automatic Control, vol. 57, pp. 1347–1361, 2012.
[17] P. J. Antsaklis, “Stable proper nth-order inverses,” IEEE Transactions
on Automatic Control, vol. 23, pp. 1104–1106, 1978.
[18] S. Gillijns, “Kalman filtering techniques for system inversion and data
assimilation,” Ph.D. dissertation, Katholieke University Leuven, 2007.
[19] S. Kirtikar, H. Palanthandalam-Madapusi, E. Zattoni, and D. S. Bernstein, “l-delay input and initial-state reconstruction for discrete-time
linear systems,” Circuits, Systems, and Signal Processing, vol. 30, pp.
233–262, 2011.
[20] F. Shi and R. J. Patton, “Simultaneous state and fault estimation for
descriptor systems using an augmented PD observer,” in Proceedings of
the 19th IFAC World Congress, Cape Town, Sowth Africa, 2014, pp.
8006–8011.
[21] G. Marro and E. Zattoni, “Unknown-state, unknown-input reconstruction in discrete-time nonminimum-phase systems: geometric methods,”
Automatica, vol. 46, pp. 815–822, 2010.
[22] A. Chiuso, “The role of vector autoregressive modeling in predictor
based subspace identification,” Automatica, vol. 43, pp. 1034–1048,
2007.
[23] G. van der Veen, J. van Wingerden, M. Bergamasco, M. Lovera,
and M. Verhaegen, “Closed-loop subspace identification methods: an
overview,” IET Control Theory and Applications, vol. 7, pp. 1339–1358,
2013.
[24] G. M. Kuersteiner, “Automatic inference for infinite order vector autoregressions,” Econometric Theory, vol. 21, pp. 85–115, 2005.
[25] P. J. Moylan, “Stable inversion of linear systems,” IEEE Transactions
on Automatic Control, vol. 22, pp. 74–78, 1977.
[26] M. Hou and R. J. Patton, “Input observability and input reconstruction,”
Automatica, vol. 34, pp. 789–794, 1998.
[27] Y. Xiong and M. Saif, “Unknown disturbance input estimation based on
a state functional observer design,” Automatica, vol. 39, pp. 1389–1398,
2003.
[28] T. Katayama, Subspace Methods for System Identification. London:
Springer-Verlag, 2005.
[29] S. Gillijns and B. D. Moor, “Unbiased minimum-variance input and
state estimation for linear discrete-time systems,” Automatica, vol. 43,
pp. 111–116, 2007.
[30] ——, “Unbiased minimum-variance input and state estimation for linear
discrete-time systems with direct feedthrough,” Automatica, vol. 43, pp.
934–937, 2007.
[31] A. Chiuso, “On the relation between CCA and predictor-based subspace
identification,” IEEE Transactions on Automatic Control, vol. 52, pp.
1795–1812, 2007.
[32] H. Kwakernaak and R. Sivan, Linear Optimal Control Systems. New
York: John Wiley, 1972.
[33] J. B. Burl, Linear Optimal Control: H2 and H∞ Methods. California:
Addison-Wesley, 1998.
[34] M. Darouach and M. Zasadzinski, “Unbiased minimum variance estimation for systems with unknown exogenous inputs,” Automatica, vol. 33,
pp. 717–719, 1997.
[35] Y. Wan and H. Ye, “Data-driven diagnosis of sensor precision degradation in the presence of control,” Journal of Process Control, vol. 22, pp.
26–40, 2012.
[36] T. Kailath, A. Sayed, and B. Hassibi, Linear Estimation. Englewood
Cliffs, NJ: Prentice-Hall, 2000.

