1

Particle Filter-Based Fault Diagnosis of Nonlinear
Systems Using a Dual Particle Filter Scheme*
Najmeh Daroogheh1 , Nader Meskin2 and Khashayar Khorasani1

arXiv:1606.08806v1 [] 28 Jun 2016

Abstract
In this paper, a dual estimation methodology is developed for both time-varying parameters and
states of a nonlinear stochastic system based on the Particle Filtering (PF) scheme. Our developed
methodology is based on a concurrent implementation of state and parameter estimation filters as
opposed to using a single filter for simultaneously estimating the augmented states and parameters. The
convergence and stability of our proposed dual estimation strategy are shown formally to be guaranteed
under certain conditions. The ability of our developed dual estimation method is testified to handle
simultaneously and efficiently the states and time-varying parameters of a nonlinear system in a context
of health monitoring which employs a unified approach to fault detection, isolation and identification
is a single algorithm. The performance capabilities of our proposed fault diagnosis methodology is
demonstrated and evaluated by its application to a gas turbine engine through accomplishing state
and parameter estimation under simultaneous and concurrent component fault scenarios. Extensive
simulation results are provided to substantiate and justify the superiority of our proposed fault diagnosis
methodology when compared with another well-known alternative diagnostic technique that is available
in the literature.

I. INTRODUCTION
Systems state estimation is a fundamental problem in control, signal processing, and fault
diagnosis fields [1]. Investigations on both linear and nonlinear state estimation and filtering in
stochastic environments have been an active area of research during the past several decades.
Linear state estimation methods use a simpler representation of an actual nonlinear system and
can provide an acceptable performance only locally around an operating point and in the steady
state operational condition of the system. However, as nonlinearities in the system dynamics
become dominant, the performance of linear approaches deteriorates and linear algorithms will
not necessarily converge to an accurate solution. Although an optimal state estimation solution for
linear filtering methods exists, nonlinear filtering methods suffer from generating sub-optimal or
near-optimal solutions. Consequently, investigation of nonlinear estimation and filtering problems
remain a challenging research area.
*This publication was made possible by NPRP grant No. 4-195-2-065 from the Qatar National Research Fund (a member of
Qatar Foundation). The statements made herein are solely the responsibility of the authors.
1
N. Daroogheh and K. Khorasani are with the Department of Electrical and Computer Engineering, Concordia University,
Montreal, Quebec, Canada n_daroog@encs.concordia.ca and kash@ece.concordia.ca
2
N. Meskin is with the Department of Electrical Engineering, Qatar University, Doha, Qatar

nader.meskin@qu.edu.qa
March 13, 2018

DRAFT

2

Numerous studies have been conducted in the literature to solve and analyze standard nonlinear
filtering problems [2]–[8]. These methods can be broadly categorized into [7]: (a) linearization methods (extended Kalman filters (EKF)) [2], (b) approximation methods using finitedimensional nonlinear filters [3], (c) particle filter (PF) methods as one of the most popular
Bayesian recursive methods for state estimation [4], (d) classical partial differential equation
(PDE) methods for approximating a solution to the Zakai equation [5], (e) Wiener chaos expansion methods [6], (f) moment methods [7], and (g) high dimensional nonlinear Kalman filter
methods known as the Cubature Kalman filters [8].
One of the most important recent applications of nonlinear filtering methods is in the area of
fault diagnosis of dynamical systems that can include fault detection, isolation, and identification
(FDII) modules. Diagnosis methods that are based on linearization techniques suffer from poor
detection and high rates of false alarms. Therefore, Monte Carlo filtering approach based on
particle filters was first proposed in [9] to address the fault detection and isolation problem in
nonlinear systems. In this work, the negative log-likelihood, which is calculated for a predefined
time window, is considered as a measure for the fault detection. The fault isolation was achieved
by using the augmentation of the fault parameters vector to the system states to perform the
estimation task. However, the augmented state space model tends to increase the dimensionality
of the model and as a result increases the number of required particles for achieving a sufficiently
accurate result. For decreasing the computational burden of this method, the augmented model
is used only after the fault detection stage and for only the fault isolation stage. An external
covariance adjustment loop was added to this augmented model in [10] to enable the estimation
algorithm to track changes in the system parameters in case of fault occurrences.
The combination of a particle filtering algorithm and the log-likelihood ratio (LLR) test in the
multiple model environment, has led to the development of sensor/actuator FDI scheme in [11]
for a general class of nonlinear non-Gaussian dynamical systems but with the assumption of full
state measurements. The fault detection problem recently is addressed for a mobile robot based
on the combination of the negative LLR test and particle filtering approach in [12]. However,
both of the methods in [11], and [12] suffer from the high computational burden for on-line
implementation of the algorithms. Hence, the idea of parallelized particle filters for on-line fault
diagnosis is introduced in [12] to improve the performance of the algorithm.
A PF-based robust navigation approach was proposed in [13] to address multiple and simultaneous faults occurrences in both actuators and sensors in an underwater robot where an
anomaly is modeled by a switching-mode hidden Markov system. The component and actuator
fault detection and isolation of a point mass satellite was tackled in [14] by introducing several
particle filters that run in parallel and each rejects a different subset of the faults.
Generally, the main issues with applying standard particle filters to the problem of fault
diagnosis problem can be stated as follows [15]: (i) False diagnosis decisions due to low
probabilities of transitions to fault states when there are fewer samples of states, and (ii)
The exponential growth of the required samples for accurately approximating the a posteriori
distributions as dimensionality of the estimation problem increases. The risk-sensitive PF is
March 13, 2018

DRAFT

3

introduced to address the first problem and the variable resolution PF is developed to overcome
the second problem in [16]. Moreover, the Gaussian PF (GPF) is also introduced in [17] as
an efficient algorithm for performing fault diagnosis of hybrid systems faster than traditional
methods that are based on PFs. Finally, the sample impoverishment problem in particle filters
due to fault occurrence in a hybrid system is addressed in [18]. The developed algorithm enables
the PF method to be implemented by fewer number of particles even under faulty conditions.
In this work, our objective is to address the above problems by specifically utilizing particle
filters. This is accomplished through development of a PF-based dual state/parameter estimation
scheme for the system component fault diagnosis. The developed dual estimation scheme relies
on an on-line estimation of the system states as well as the system health parameters. These
parameters are not necessarily fixed and their time variations are governed and influenced by
the fault vector that directly affects them.
The on-line estimation of the system time-varying parameters by using particle filters is a
challenging and active area of research. There are two main classes of PF-based parameter
estimation algorithms (for on-line as well as off-line implementations) [19] known as Bayesian
and maximum likelihood (ML) approaches. In the Bayesian approach, a prior distribution is
considered for the unknown parameters and the posterior distribution of the parameters is approximated given the observations [20], [21], whereas in the ML approach the estimated parameter is
the maximizing argument of the likelihood function given the observations [22]–[25]. In the ML
framework for parameter estimation, the maximization of any cost function can be performed
based on gradient-based search methods [22]. On the other hand, expectation maximization
(EM) methods are only applicable for maximization of the likelihood functions [25]. However,
EM methods are not suitable for on-line applications due to their high computational cost for
implementation. The recursive maximum likelihood method (RML) is recognized as a promising
method for on-line parameter estimation based on a stochastic gradient algorithm [23]. In order
to avoid the direct computation of the likelihood function gradient, an alternative method is
proposed in [26] that is known as the gradient-free ML parameter estimation. Despite the above,
the on-line ML methods suffer from the practical point of view of slow convergence rates and
requiring large number of particles to achieve accurate estimates [27].
In the Bayesian framework, on-line implementation of particle filter-based parameter estimation algorithms are computationally intensive [28]. A general method that is capable of
simultaneously estimating the static (i.e. , constant or fixed) parameters and time-varying states
of a system is developed in [29]. The work is based on the sequential Monte Carlo (SMC)
method in which an artificial dynamic evolution model is considered for the unknown model
parameters. In order to overcome the degeneracy concerns arising from the particle filtering,
kernel smoothing technique as a method for smoothing the approximation of the parameters
conditional density has been utilized in [20]. The estimation algorithm is further improved by
re-interpretation of the artificial evolution algorithm according to the shrinkage scaling concept.
However, the proposed method in [20] is only applicable for estimating fixed parameters of the
system and it uses the augmented state/parameter vector for the estimation task.
March 13, 2018

DRAFT

4

Our main goal in this work is to extend the Bayesian framework that is proposed in [20] for
both state and parameter estimation by invoking our proposed modified artificial evolution law
that is inspired from the prediction error (PE) concept. This modification enables the scheme to
track and estimate the time-varying parameters of the system. The PE is used for both off-line
and on-line system identification by incorporating a quadratically convergent scheme [30], [31].
The concept of the kernel mean shrinkage (KS) is also used in our proposed scheme to avoid
over-dispersion in the variance of the estimated parameters.
Component fault diagnosis of dynamical systems can be achieved through dual state/parameter
estimation methods [32], [33] by first detecting the faulty competent and then by determining its
severity and isolating its location. In this work, our proposed dual state and parameter estimation
strategy is applied for performing component fault diagnosis of a gas turbine engine [34]. This
is accomplished by utilizing a multiple-model approach as developed in [34], [35]. Although
the multiple-model based approach can detect and isolate faults based on a set of predefined
faulty models, the utilization of parameter estimation schemes embedded in our proposed dual
state/parameter estimation methodology enables the diagnostic scheme to not only detect and
isolate the faults but also accurately identify the type and severity of simultaneous fault scenarios
in the gas turbine system. The simultaneous detection, isolation, and identification of faults are
the main feature that distinguishes the current work with the works conducted in [34], [35].
Moreover, the main methodology that is developed in these references is based on Kalman
filters as opposed to the particle filtering methodology that is developed on the current work.
Based on the above discussion, the main contributions of this paper can be summarized as that
of utilizing nonlinear Bayesian and Sequential Monte Carlo (SMC) methods to develop, design,
analyze, and implement a unified framework for both states and parameters estimation as well as
fault diagnosis of nonlinear systems. Our methodology is based on solving the Bayesian recursive
relations through SMC methods. An on-line parameter estimation scheme is developed based on
a modified artificial evolution concept by using the particle filters (PF) approach. Specifically,
by using the prediction error to correct the time-varying changes in the system parameters, a
novel methodology is proposed for parameter estimation of nonlinear systems based on the PF.
Specifically, in the implementation of our proposed scheme, a dual structure for both state and
parameter estimation is developed within the PF framework. In other words, the hidden states
and variations of the system parameters are estimated through operating two concurrent filters.
Convergence and stability properties of our proposed dual estimation strategy are shown to be
guaranteed formally under certain conditions.
The remainder of this paper is organized as follows. In Section II, the statement of the
nonlinear filtering problem is presented. Our proposed dual state/parameter estimation scheme
is developed in Section III, in which state and parameter estimation methods are first developed
concurrently and subsequently integrated together for simultaneously estimating the system
states and parameters. The stability and convergence properties of the proposed schemes under
certain conditions are also provided in Section III. Our proposed fault diagnosis framework and
formulation are also provided in Section III. In Section IV, extensive simulation results and
March 13, 2018

DRAFT

5

case studies are provided to demonstrate and justify the merits of our proposed method for fault
diagnosis of a gas turbine engine under simultaneous and concurrent component faults. Finally,
the paper is concluded in Section V.
II. PROBLEM STATEMENT
The problem under consideration is to obtain an optimal estimate of states as well as timevarying parameters of a nonlinear system whose dynamics is governed by a discrete-time stochastic model,
xt+1 = ft (xt , θt , ωt ),

(1)

yt = ht (xt , θt ) + νt ,

(2)

where xt ∈ Rnx is the system state, t ∈ N, ft : Rnx × Rnθ × Rnω −→ Rnx is a known nonlinear
function, θt ∈ Rnθ is an unknown and possibly time-varying parameter vector governed by
an unknown dynamics. The function ht : Rnx × Rnθ −→ Rny is a known nonlinear function
representing the map between the states, parameters and the system measurements, and ωt and
νt are uncorrelated stochastic process and measurement noise sequences with covariance matrices
Lt and Vt , respectively. The following assumption is made regarding the dynamical system (1)
and (2).
Assumption A1. The vector {xt , θt } ranges over a compact set denoted by DN , for which
the functions ft (xt , θt , ωt ) and ht (xt , θt ) are continuously differentiable with respect to the state
xt as well as the parameter θt .
The main objective of the dual state and parameter estimation problem is to approximate the
following conditional expectations:
Z
E(φ1 (xt )|y1:t , θt−1 ) =

φ1 (xt )p(xt |y1:t , θt−1 )dxt ,

(3a)

φ2 (θt )p(θt |y1:t , xt )dθt ,

(3b)

Z
E(φ2 (θt )|y1:t , xt ) =

where y1:t = (y1 , y2 , ..., yt ) denotes the available observations up to time t, φ1 : Rnx → R and
φ2 : Rnθ → R are functions of states and parameters, respectively, that are to be estimated. The
conditional probability functions p(xt |y1:t , θt−1 )dxt and p(θt |y1:t , xt )dθt are to be approximated
by the designed particle filters (PFs) through determining the filtering distributions according to
p̂N (xt |y1:t , θt−1 )dxt =

N
X

wx(i)t δx(i) (dxt ),
t

i=1

p̂N (θt |y1:t , xt )dθt =

N
X

(4)
(j)

wθt δθ(j) (dθt ),
t

j=1

where the subscript N in p̂N (.) implies that the state/parameter conditional probability distribu(i)
(i)
tions are obtained from N particles. Each state particle xt has a weight wxt and each parameter

March 13, 2018

DRAFT

6
(j)

(j)

particle θt has a weight wθt , where δ(.) denotes the Dirac-delta function mass that is positioned
at xt or θt .
Based on the approximations used in equation (4), our goal is to address the convergence
properties of the subsequently designed estimators to their true optimal estimates and also to
develop and demonstrate under what conditions this convergence remains valid.
III. P ROPOSED D UAL S TATE /PARAMETER E STIMATION AND FAULT D IAGNOSIS
F RAMEWORK
In this section, the main theoretical framework for our proposed dual state/parameter filtering
as well as the fault diagnosis methodology of the nonlinear system (1) and (2) are introduced
and developed.
A. Dynamic Model in Presence of Time-Varying Parameters
Our first task is to represent the model (1) and (2) into another framework for our subsequent
theoretical developments. Let (Ω, F, P ) denote the probability space on which the three real
vector-valued stochastic processes X = {Xt , t = 1, 2, ...}, Θ = {Θt , t = 1, 2, ...}, and Y =
{Yt , t = 1, 2, ...} are defined. The nx -dimensional process X describes the evolution of the hidden
states, the nθ -dimensional process Θ describes the evolution of the hidden system parameters
that are conditionally independent of the states, and the ny -dimensional process Y denotes the
observation process of the system.
The processes X and Θ are Markov processes with the associated initial state and parameter
X0 and Θ0 , respectively. They are drawn from the initial distributions πx0 (dx0 ) and πθ0 (dθ0 ),
respectively. The dynamic evolution of states and parameters are modeled by the Markov transition kernels Kx (dxt |xt−1 , θt−1 ) and Kθ (dθt |θt−1 , xt ), that also admit densities with respect to
the Lebesgue measure 1 , such that
Z
P (Xt ∈ A1 |Xt−1 = xt−1 , Θt−1 = θt−1 ) =

Kx (xt |xt−1 , θt−1 )dxt ,

(5)

A1

Z
P (Θt ∈ A2 |Θt−1 = θt−1 , Xt = xt ) =

Kθ (θt |θt−1 , xt )dθt ,

(6)

A2

for all A1 ∈ B(Rnx ) and A2 ∈ B(Rnθ ), where B(Rnx ) and B(Rnθ ) denote the Borel σ-algebra
on Rnx and Rnθ , respectively. The transition kernel Kx (xt |xt−1 , θt−1 ) is a probability distribution
function (pdf) that follows the pdf of the stochastic process in process (1). The probability density
function for approximating the parameter kernel transition Kθ (θt |θt−1 , xt ) is to be provided in
the subsequent subsections.
Given the states and parameters, the observations Yt are conditionally independent and have
the marginal distribution with a density with respect to the Lebesgue measure as given by,
R
P (Yt ∈ B|Xt = xt , Θt = θt ) = B ρ(yt |xt , θt )dyt ,
(7)
1

The transition kernel K(dxt |xt−1 ) admits density with respect to the Lebesgue measure if one can write P (Xt ∈
dxt |Xt−1 = xt−1 ) = K(dxt |xt−1 ) = K(xt |xt−1 )dxt .

March 13, 2018

DRAFT

7

xˆ
t 1|t 1

Time update for
state estimation
filter

Observations

ˆ

t 1|t 1

xˆ
t|t 1

Measurement
update for state
estimation filter

xˆ
t|t

yt
First step measurement
update for parameter
estimation filter

~
t |t



Second step
measurement update
for parameter
estimation filter

ˆ

t|t

Fig. 1. The schematic of the dual particle filter.

where ρ(yt |xt , θt ) is a probability density function that follows the probability density function
of the stochastic process in equation (2).
In the dual state/parameter estimation framework, at first the state xt is estimated (which is
denoted by x̂t|t ). The estimated value at time t is then used to estimate the parameter θt at
time t (which is denoted by θ̂t|t ). In the Bayesian framework for parameter estimation, the prior
evolution of parameters are not specified, therefore it is necessary to consider a given evolution
for the parameters in order to design an estimation filter. In our proposed dual structure for
the state estimation filter, first the parameters are assumed to be constant at time t − 1 at their
estimated value θ̂t−1|t−1 , and then for the parameter estimation filter they are evolved to the
next time instant by applying an update law that is based on the prediction error (PE) method.
The details regarding our proposed methodology are presented in the subsequent subsections in
which the filtering of states and parameters are fully described and developed.
B. The Dual State/Parameter Estimation Framework
In our proposed dual state/parameter estimation framework, two filters are runing concurrently.
At every time step, the first PF-based state filter estimates the states by using the current available
estimate of the parameters, θ̂t−1|t−1 , whereas the second PF-based parameter filter estimates the
unknown parameters by using the current estimated states, x̂t|t . The developed schematic is
shown in Figure 1.
In our dual estimation framework, the well-known maximum a posteriori (MAP) solution
corresponding to the marginal estimation methods based on the decoupled approach is used
for solving the dual estimation problem [36]. In this method, the joint state/parameter marginal
density p(xt , θt |y1:t ) is expressed as
p(xt , θt |y1:t ) = p(xt |θt , y1:t )p(θt |y1:t ),

(8)

where p(xt |θt , y1:t ) and p(θt |y1:t ) denote the state and parameter marginal densities, respectively.
Assuming that the variations of parameters are slow when compared to the system state time
March 13, 2018

DRAFT

8

variations, one can use the approximation θt ≈ θt−1 , so that the joint marginal density is
approximated as
p(xt , θt |y1:t ) ≈ p(xt |θt−1 , y1:t )p(θt−1 |y1:t ).

(9)

Our ultimate goal is to maximize the two marginal distribution terms in expression (9) separately
according to the decoupled approach in [36] as follows
x̂t|t = argmaxxt p(xt |θt−1 , y1:t ), θ̂t|t = argmaxθt−1 p(θt−1 |y1:t ).

(10)

In the above decoupled methodology, one attribute is optimized at a time by keeping the other
attribute fixed and then alternating them. Associated with optimization of both marginal distributions, different cost functions can be chosen [36]. For developing a dual extended Kalman filter,
corresponding to specific cost functions of the parameter marginal density, various estimation
methods have been proposed in the literature [36], [37]. For example, the maximum likelihood
(ML) and prediction error approaches are selected for marginal estimations. The main motivation
for choosing these two approaches is due to the fact that one considers to maximize only the
marginal density p(θt−1 |y1:t ) as opposed to the joint density p(xt , θt−1 |y1:t ). However, in order
to maximize the parameter marginal density, it is also necessary to generate state estimates that
are produced by maximizing the state marginal density p(xt |θt−1 , y1:t ).
It should be noted that in marginal estimation methods no explicit cost function is considered
for maximization of the state marginal distribution, since the state estimation is only an implicit
step in marginal approaches and the joint state/parameter cost is used that may have variety of
forms in different filtering algorithms [36]. In our proposed dual particle filtering framework,
p(xt |θt−1 , y1:t ) is approximated by the state filtering distribution p̂N (xt |θt−1 , y1:t ) from equation
(4). Next, the prediction error cost function is chosen for maximization of the parameter marginal
density, where this cost function is implemented through a PE approach in order to attain a less
computational cost [30].
In the subsequent subsections, specific details regarding the concurrent state and parameter
estimation filters design and development are provided.
C. The State Estimation Problem
For designing the state and parameter estimation filters, our main objectives are to approximate the integrals in equations (5) and (6) by invoking the particle filter (PF) scheme as well
as to approximate the estimate of the conditional state and parameter distributions. Consider
R
πxt|t−1 (dxt ) = Rnx πxt−1|t−1 (dxt−1 )Kx (dxt |xt−1 , θt−1 ) to denote the a priori state estimation
distribution before the observation at time t becomes available, and πθt−1|t−1 (dθt−1 ) to denote
the marginal distribution of the parameter at time t − 1. The a posteriori state distribution after
the observation at the instant t becomes available is obtained according to the following rule,
πxt|t (dxt ) ∝ ρ(yt |xt , θt−1 )πxt|t−1 (dxt )πθt−1|t−1 (dθt−1 ).

March 13, 2018

(11)

DRAFT

9

In the above it is assumed that θ̂t−1|t−1 is known for this filter. Therefore, the last distribution
in the right hand side of equation (11) is set to one.
The particle filter (PF) procedure for implementation of the state estimation and for determining πxt|t (dxt ) consists of two main steps, namely (a) the prediction step (time update step),
and (b) the measurement update step. Consider one states in the N particles at time t. The
prediction step utilizes the knowledge of the previous distribution of the states as well as the
(i)
previous parameter estimate, these are denoted by {x̂t−1|t−1 , i = 1, ..., N } (corresponding to N
estimated state particles that follow the distribution πxt−1|t−1 (dxt−1 )) and θ̂t−1|t−1 , respectively, as
well as the process model given by equation (1). In other words, the prediction step is explicitly
governed by the following equations for i = 1, ..., N , namely
(i)

(i)

(i)

(i)

(i)

x̂t|t−1 = ft (x̂t−1|t−1 , θ̂t−1|t−1 , ωt ),

(12a)

ŷt|t−1 = ht (x̂t|t−1 , θ̂t−1|t−1 ),
Σx̂t|t−1 =

(i)
(x̂t|t−1

(12b)

N
N
1 X (i) T
1 X (i)
(i)
−
−
)(x̂
) ,
x̂
x̂
N i=1 t|t−1 t|t−1 N i=1 t|t−1

(i)

(12c)

(i)

where ωt denotes the process noise related to each particle x̂t|t−1 and is drawn from the noise
(i)
distribution with the probability distribution function pωt (.), and x̂t|t−1 denotes the independent
(i)
samples generated from equation (12a) for i = 1, ..., N particles. Moreover, ŷt|t−1 denotes the
(i)
independent samples of the predicted outputs that are evaluated at x̂t|t−1 samples, and Σx̂t|t−1
denotes the a priori state estimation covariance matrix.
For the first step, the one-step ahead prediction distribution known as the a priori state
estimation distribution is now given by,
N
1 X
N
π̃xt|t−1 (dxt ) ,
δ (i) (dxt ),
(13)
N i=1 x̂t|t−1
For the second step, the information on the present observation yt is used. This results in
approximating πxt|t (dxt ), where θ̂t−1|t−1 is considered to be given from a parameter estimation
(i)
filter and obtained from the distribution πθNt−1|t−1 (dθt−1 ). Consequently, the particle weights wxt
(i)

are updated by the likelihood function (the importance function) according to wxt ∼ pνt (yt −
(i)
(i)
ŷt|t−1 ) = ρ(yt |x̂t|t−1 , θ̂t−1|t−1 ), where pνt (.) denotes the probability distribution function of the
(i)
additive noise of the output and is evaluated at yt − ŷt|t−1 .
In this work, since our ultimate goal is in developing a fault diagnosis algorithm that is practically stable, the structure of regularized particle filters (RPF) is chosen which is quite common in
many practical applications [38]. The characteristics of the RPFs are related to the fact that they
are capable of transforming the discrete-time approximation of the a posteriori state estimation
distribution πxNt|t (dxt ) into a continuous-time one. Consequently, the resampling step is modified
in such a manner that the new resampled particles are obtained from an absolutely continuous(i)
(i)
time distribution with N different locations x̂t|t from that of x̂t|t−1 [39]. Therefore, where the
(i)

(k)

(k)

probability for taking the k-th particle is P (x̂t|t = x̂t|t−1 ) = w̃xt ,
March 13, 2018

(k)

ρ(yt |x̂t|t−1 ,θ̂t−1|t−1 )
,
PN
(k)
k=1 ρ(yt |x̂t|t−1 ,θ̂t−1|t−1 )
DRAFT

10
(k)

and w̃xt for k = 1, ..., N denotes the normalized particle weights. In other words, the particle
selection in the resampling step is performed for particles that have higher probabilities of
(k)
ρνt (yt − ŷt|t−1 ).
For the resampling step, two main choices can be considered that are known as (i) Bayesian
bootstrap, and (ii) Sampling importance resampling (SIR) [4]. Although both approaches are
applicable for this filter, the bootstrap method is chosen in this paper. Therefore, the a posteriori
state estimation distribution is approximated by π̃xNt|t (dxt ) before one performs the resampling
by using the RPF structure [39], and by πxNt|t (dxt ) after one performs the resampling that is
provided below,
Nreg N
X
X
|A−1 | 1
regl
(i)
N
π̃xt|t (dxt ) ≈
w̃x(i)t nt x K( A−1
− x̂t|t−1 )),
t (xt
b
b
l=1 i=1
(i)

w̃x(i)t

ρ(yt |x̂t|t−1 , θ̂t−1|t−1 )

, PN

(i)

i=1 ρ(yt |x̂t|t−1 , θ̂t−1|t−1 )

πxNt|t (dxt ) =

(14)

,

N
N
1 X
1 X (i)
δx̂(i) (dxt ) → x̂t|t =
x̂ ,
N i=1 t|t
N i=1 t|t

reg

where xt l , l = 1, .., Nreg denotes the regularized state vector that is evaluated at Nreg points
that are obtained from the absolutely continuous-time distribution of the particles as given by
(N )

(1)

Xt|t−1 = [x̂t|t−1 , ..., x̂t|t−1 ],
reg1

xt

= min(Xt|t−1 ) − std(Xt|t−1 ),

regN
xt reg

regNreg

dxreg = (xt
reg
xt l

=

(15)

= max(Xt|t−1 ) + std(Xt|t−1 ),
reg
xt l−1

reg1

− xt

)/(Nreg − 1),

+ dxreg , l = 2, ..., Nreg ,
(i)

where std denotes the first standard deviation of the particles from their mean. Hence, {x̂t|t }N
i=1
is obtained from the continuous-time distribution through the regularization kernel K that is
considered to be a symmetric density function on Rnx [39]. The matrix At in equation (14) is
(i)
chosen to yield a unit covariance value in the new x̂t|t population and At AT
t = Σx̂t|t−1 . The
constant b denotes the optimal bandwidth of the kernel, and x̂t|t denotes the a posteriori state
estimation at time t.
We are now in a position to introduce our overall particle filter (PF) scheme for implementing
the state estimation filter. Our goal for proposing this algorithm is to ensure that an approximation
P
N
to E(φ(xt )|y1:t , θt−1 ) by φ(xt ) = xt takes x̂t|t ∼ πxNt|t (dxt ) = N1 N
i=1 δx̂(i) (dxt ), where πxt|t (dxt )
t|t

(i)

(i)

N
denotes the a posteriori distribution of {x̂t|t }N
i=1 (after the resampling from {x̂t|t−1 }i=1 ), that is
P
(i)
given by x̂t|t = N1 N
i=1 x̂t|t . The estimated output from the state estimation filter is also given
by ŷt = ht (x̂t|t , θ̂t−1|t−1 ).

March 13, 2018

DRAFT

11

The State Estimation Particle Filter Scheme
(i)

1) Initialize the PF scheme with N particles, {x0 }N
i=1 ∼ πx0 (dx0 ) and the parameters θ0 (the
mean of the parameter initial distribution πθ0 (dθ0 )).
(i)
2) Draw ωt ∼ pωt (.), where pωt (.) denotes a given distribution for the process noise in the
(i)
filter, and then predict the state particles x̂t|t−1 according to equation (12a).
(i)
(i)
(i)
3) Compute ŷt|t−1 from equation (12b) to obtain the importance weights {wxt }N
i=1 as wxt =
(i)

4)

5)
6)
7)

(i)

w

(i)

ρ(yt |x̂t|t−1 , θ̂t−1|t−1 ), i = 1, ..., N, and normalize them to w̃xt = PN xt (i) .
i=1 wxt
Resampling: Draw N new particles with the replacement for each i = 1, ..., N , according to
(i)
(k)
(k)
(i)
P (x̂t|t = x̂t|t−1 ) = w̃xt , k = 1, ..., N , from the regularized kernel K where x̂t|t ∼ π̃xNt|t (dxt )
as given by equation (14).
Calculate x̂t|t from the conditional distribution that is given by equation (14),
P
PN (i)
(i)
1
πxNt|t (dxt ) = N1 N
i=1 δx̂(i) (dxt ) with equally weighted x̂t|t as x̂t|t = N
i=1 x̂t|t .
t|t
Update the parameters from the parameter estimation filter (to be specified in the next
subsection).
Set t := t + 1 and go to Step 2.

Following the implementation of the above state estimation filter, the parameter estimation
filter that is utilized for adjusting the parameters is now described in detail in the next subsection.

D. The Parameter Estimation Problem
One of the main contributions of this paper is to develop a novel PF-based parameter estimation
filter within our proposed dual state/parameter estimation framework by utilizing the prediction
error (PE) concept. For this methodology it is assumed that the a priori distribution of the timevarying developed is not known. Moreover, the estimated states that are generated by the state
estimation filter provided in the previous subsection will be used. Therefore, it is imperative that
one considers a dynamical model associated with the parameters evolution in order to estimate
the density function πθt|t (dθt ).
The most common dynamical model that is considered for the parameter propagation (in
case of the system with constant parameters) is the conventional artificial evolution law. In this
representation small random disturbances are added to the state particles (parameters) between
each consecutive time step [29]. However, in our work, the conventional update law for the
parameters is modified to include the output prediction error as a multiplicative term to allow
one to deal with time variations in the parameters that can affect the system output.
In order to derive the parameter update law, an algorithm based on the prediction error (PE)
¯ t−1 )
method is proposed by minimizing the expectation of a quadratic performance index J(θ
with respect to θt−1 . This is due to the fact that our parameter estimation algorithm for obtaining
the distribution of the a posteriori parameter estimate is based on the kernel smoothing that
uses the shrinkage of the particle locations. This method attempts to force the particles towards
their mean from the previous time step, i.e. the estimated value of θt−1 , and is denoted by
March 13, 2018

DRAFT

12

θ̂t−1|t−1 (before adding noise to the particles). This is also used in the state estimation filter for
approximating x̂t|t . Therefore, our goal is to investigate the convergence properties of θ̂t−1|t−1
whose boundedness ensures the boundedness of θ̂t|t . Towards this end, the performance index
R
¯ t−1 )|y1:t−1 , xt ) = J(θ
¯ t−1 )p(θt−1 |y1:t−1 , xt )dθt−1 , where the integral is
is now selected as E(J(θ
¯ t−1 )|y1:t−1 , xt ) ≈ 1 PN J(
¯ θ̂(j)
approximated in the PF by E(J(θ
j=1
t−1|t−1 ).
N
(j)
¯ θ̂
The term J(
t−1|t−1 ) now represents a quadratic function of the output prediction error related
(j)
to each particle j, j = 1, ..., N . The prediction error is now defined according to (t, θ̂t−1|t−1 ) ,
(j)
(j)
(j)
t = yt − ht (x̂t|t , θ̂t−1|t−1 ), where θ̂t−1|t−1 denotes the particle related to the estimated value
?
of the parameter whose true value is denoted by θt−1
(this is clearly assumed to be unknown).
P
(j)
(j)
τ
=t
1
¯ θ̂
Therefore, we define J(
τ =t−κ E(Q((τ, θ̂t−1|t−1 ))), in which the expectation is
t−1|t−1 ) = κ
taken over the observation sequence of κ samples. Let us now select the quadratic criterion
(j)
Q((t, θ̂t−1|t−1 )) as
1
(j)
(j)
(j)
Q((t, θ̂t−1|t−1 )) = (t, θ̂t−1|t−1 )T (t, θ̂t−1|t−1 ).
(16)
2
The following modified artificial evolution law is now proposed for the parameter update in
the particle filters for generating j = 1, ..., N parameter particles that correspondingly determine
(j)
the distribution from which the a priori parameter estimate θ̂t|t−1 is considered to be the same
(j)
as θ̂t−1|t−1 , and the a posteriori parameter estimate is obtained in two steps that are denoted by
(j)
(j)
θ̃t|t and θ̂t|t , respectively.
(j)
In the first step (the second step is described in the next page) for θ̃t|t one gets
(j)

(j)

(j)

(j)

(j)

mt = θ̂t−1|t−1 + γt Rt ψt (t, θ̂t−1|t−1 ),
(j)
θ̃t|t

where ψt =

∂ ŷt
∂ θ̂t−1|t−1

=
=

(j)
Amt

+ (I − A)m̄t−1 +

∂ht (x̂t|t ,θ̂t−1|t−1 )
∂ θ̂t−1|t−1

(j)
ζt ,

m̄t−1

(17a)
N
1 X (j)
=
,
θ̂
N j=1 t−1|t−1

(17b)
(j)

(j)

, which when evaluated at θ̂t−1|t−1 is denoted by ψt , γt
(j)

denotes the step size design parameter, ζt ∼ N (0, (I − A2 )Vθ̂t−1|t−1 ) denotes the zero-mean
normal increment particles to the parameter update law at each time step with the covariance
matrix (I −A2 )Vθ̂t−1|t−1 through the use of the kernel smoothing concept, A denotes the shrinkage
matrix, and Vθ̂t−1|t−1 denotes the covariance of the parameter estimates in the previous time
step t − 1. The kernel shrinkage algorithm attempts to force the parameter particles distribution
towards its mean in the previous time instant that was denoted by m̄t−1 , by applying the shrinkage
(j)
coefficient matrix A to the obtained mt .
(j)
(j)
The processes θ̂t−1|t−1 and ζt are conditionally independent given observations up to time
q
Pl=n (l) (j)
(j)
(j) (j)T
(j)
(j)
t. Moreover, Rt = trace(Et Et ) where Et = t (θ̂t−1|t−1 ) − n1y l=1 y t (θ̂t−1|t−1 ) and
(l)

(j)

(j)

(j)

t (θ̂t−1|t−1 ) denotes the l-th element of the vector t (θ̂t−1|t−1 ). The term Rt denotes a timevarying coefficient that determines the updating direction and is a positive scalar to ensure that
(j)
the criterion (16) can be minimized by changing mt in the steepest descent direction. Therefore,

March 13, 2018

DRAFT

13
(j)

the first step estimate of the a posteriori parameter estimation particle is denoted by θ̃t|t . The
convergence of the update law (17a)-(17b) will be shown in the Subsection III-F.
The parameter update law according to (17a)-(17b) contains a term in addition to the in(j)
dependent normal increment ζt . The estimated parameter from this update law is invoked in
the PF-based parameter estimation filter to represent the distribution from which the parameter
particle population for the next time step is chosen. Therefore, the above proposed prediction
error based modified artificial evolution law enables the PF-based estimation algorithm to handle
(j)
and cope with the time-varying parameter scenarios. The time-varying term γt Rt acts as an
adaptive step size in equations (17a)-(17b), and therefore our algorithm can also be considered
as an adaptive step size scheme.
(j)
In order to ensure that the obtained θ̃t|t from the modified artificial evolution law given
by equations (17a)-(17b) remains in DN (refer to Assumption A1), the following projection
(j)
algorithm is utilized that forces θ̃t|t to remain inside DN according to the following procedure
[30],
1)
2)
3)
4)
5)
6)

Choose a factor 0 ≤ µ ≤ 1,
(j)
(j) (j)
(j)
Compute θ̆t|t := γt Rt ψt (t, θ̂t−1|t−1 ),
(j)
(j)
(j)
Construct mt := θ̂t−1|t−1 + θ̆t|t ,
(j)
If mt ∈ DN go to Step 6, else go to Step 5,
(j)
(j)
Set θ̆t|t = µθ̆t|t , and go to Step 3,
Stop.

It should be noted that the main reason for considering the above mapping is related the fact
that the actual dynamics of the parameters are not known, therefore such mapping ensures that
the assumed dynamics for the parameters based on modified artificial evolution model does not
cause instability of the entire system.
Consequently, the a priori distribution of the parameter θt is assumed to have the same
distribution as in the previous time step. On the other hand, as the present observation yt becomes
available in the measurement update step, the a posteriori distribution of the parameter is obtained
through two steps that denoted by π̃θ̃N (dθt ) and π̃θNt|t (dθt ), respectively. In what follows, more
t|t
details related to these distributions are presented.
Consider equations (17a)-(17b). The a posteriori distribution of the parameters calculated from
(j)
the distribution of the parameter particles θ̃t|t is given by,
P
π̃θ̃N (dθt ) , N1 N
(18)
j=1 δθ̃(j) (dθt ),
t|t

t|t

and the measurement equation is expressed as,
(j)

(j)

ȳt|t = ht (x̂t|t , θ̃t|t ),

(19)

(j)

where ȳt|t denotes the evaluated output that is obtained by the parameter estimation filter that is
different from the one that is obtained by the state estimation filter, as provided in the Subsection
III-C.
March 13, 2018

DRAFT

14

Now, in the second step for estimating the a posteriori parameter estimate distribution, consider
(j)
the present observation yt , so that the particle weights wθt are updated by the likelihood function
(j)
(j)
(j)
according to wθt ∼ pνt (yt − ȳt|t ) = ρ(yt |x̂t|t , θ̃t|t ). This can now be expressed by using the
(j)
PN
ρ(yt |x̂t|t ,θ̃t|t )
(j)
(j)
(j)
(dθ
),
where
w̃
,
normalized weights w̃θt as π̃θNt|t (dθt ) =
w̃
δ
(j)
P
t
(j) .
N
θt
j=1 θt θ̃
ρ(yt |x̂t|t ,θ̃t|t )
πθNt|t (dθt ) is

t|t

j=1

Following the resampling/selection step, an equally weighted particle distribution
P
obtained as πθNt|t (dθt ) = N1 N
j=1 δθ̂(j) (dθt ) for approximating πθt|t (dθt ), and the resampled
t|t

(j)

(selected) particles that are denoted by θ̂t|t follow the distribution π̃θNt|t (dθt ). Therefore, the
a posteriori parameter estimation distribution is approximated by a weighted sum of the Diracdelta masses as π̃θNt|t (dθt ) before one performs the resampling and with an equally weighted
particle distribution approximation as πθNt|t (dθt ) according to
N
X
(j)
N
π̃θt|t (dθt ) ≈
w̃θt δθ̃(j) (dθt ),
j=1

t|t

(j)

(j)
w̃θt

ρ(yt |x̂t|t , θ̃t|t )

, PN

j=1

πθNt|t (dθt )

(j)

ρ(yt |x̂t|t , θ̃t|t )

,

(20)

N
N
1 X
1 X (j)
=
δ (j) (dθt ) → θ̂t|t =
θ̂ ,
N j=1 θ̂t|t
N j=1 t|t
(j)

(j)

where w̃θt denotes the normalized parameter particle weight, {θ̂t|t }N
j=1 is obtained from the
(j)
resampling/selection step of the scheme by duplicating the particles θ̃t|t having large weights
and discarding the ones with small values to emphasize the zones with higher a posteriori
(k)
(k)
(j)
probabilities according to P (θ̂t|t = θ̃t|t ) = w̃θt , k = 1, ..., N . In our proposed filter the residual
resampling method is used to ensure that the variance reduction among the resampled particles
is guaranteed [40].
Therefore, an approximation to E(φ(θt )|y1:t , xt ) by φ(θt ) = θt takes on the form θ̂t|t ∼
P
N
N
πθt|t (dθt ) = N1 N
j=1 δθ̂(j) (dθt ), where πθt|t (dθt ) denotes the a posteriori distribution of the
t|t

(j)

parameter estimate (after performing the resampling from θ̃t|t ). The resulting estimated output
of this filter is obtained by ŷt = ht (x̂t|t , θ̂t|t ). The explicit details for implementation of the
parameter estimation filter are now provided below.
The Parameter Estimation Filter
The particle filter for implementation of the parameter estimation is described as follows:
1) Initialize the N particles for the parameters as {θ0j }N
j=1 ∼ πθ0 (dθ0 ), and use the initial values
of the states as x0 that represents the mean of the states initial distribution πx0 (dx0 ).
(j)
2) Draw ζt ∼ N (0, (I − A2 )Vθ̂t−1|t−1 ).
(j)

3) Predict θ̃t|t , j = 1, ..., N from equations (17a)-(17b) with the projection algorithm.

March 13, 2018

DRAFT

15
(j)

(j)

(j)

4) Compute the importance weights {wθt }N
j=1 , wθt = ρ(yt |x̂t|t , θ̃t|t ), j = 1, ..., N , and normal(j)

ize them to w̃θt =

(j)

wθ
t
PN
(j) .
j=1 wθ
t

(j)

(k)

5) Resampling: Draw N new particles with replacement for each j = 1, ..., N , P (θ̂t|t = θ̃t|t ) =
P
(k)
(j)
(j)
w̃θt , k = 1, ..., N , where θ̃t|t ∼ π̃θNt|t (dθt ) = N
j=1 w̃θt δθ̃(j) (dθt ).
t|t
P
6) Construct θ̂t|t from the conditional distribution πθNt|t (dθt ) = N1 N
j=1 δθ̂(j) (dθt ) with equally
t|t
P
(j)
(j)
weighted θ̂t|t as θ̂t|t = N1 N
θ̂
.
j=1 t|t
7) Set t = t + 1 and go to Step 2 of the state estimation filter as provided in the Subsection
III-C.
(j)

As stated earlier, the kernel from which the parameter particles i.e. θ̃t|t for the next time step
(j)
is chosen is a Gaussian kernel and its mean is obtained from mt and its variance is obtained
based on the kernel smoothing consideration that is provided in the next Subsection III-E. In the
subsections below, the required conditions for boundedness of the parameter transition kernel
Kθ (.) are also investigated and developed.
E. Kernel Smoothing (KS) of the Parameters
In this subsection, the kernel smoothing (KS) approach [20] is utilized to ensure that the
variance of the normal distribution which is obtained according to the modified artificial evolution
law for the parameter estimates remains bounded.
(j)
Consider the modified artificial evolution law (17a)-(17b) in which ζt is a normal zero-mean
uncorrelated random increment to the parameter that is estimated at time t − 1. If A = I, i.e.
when there is no kernel shrinkage, as t → ∞, the variance of the added evolution increases and
(j)
can therefore yield θ̃t|t in (17b) completely unreliable. This phenomenon is known as the loss
of information that can also occur between two consequent sampling times [20]. On the other
hand, since θt is time-varying, generally there will not exist an optimal value for the variance
(j)
of the evolution noise ζt that remains suitable for all times.
Consequently, the idea of the kernel shrinkage has been proposed in [20] and later updated
in [21]. In the kernel shrinkage approach [29], for the next time step one takes the mean of
the estimated parameter distribution in the particle filter according to the following normal
distribution
(j)

(j)

Kθ (dθt |θt−1 , xt ) ∼ N (Amt + (I − A)m̄t−1 , (I − A2 )
Vθ̂t−1|t−1 ),

(21)

(j)

where mt for j = 1, ..., N, is obtained from (17a). By utilizing this kernel shrinkage rule, the
resulting normal distribution retains the mean m̄t−1 and has the appropriate variance for avoiding
over-dispersion relative to the a posteriori sample. The kernel shrinkage forces the parameter
(j)
samples towards their mean before the noise ζt is added. In our proposed approach the changes
due to the parameter variations are considered in the mean of the parameter estimate distribution
March 13, 2018

DRAFT

16

through the modified artificial evolution rule. Consequently, the mean of the distribution, i.e.
m̄t−1 , itself is time-varying and the kernel shrinkage ensures a smooth transition in the estimated
parameters even when they are subjected to changes. To eliminate the information loss effect,
by taking the variance from both sides of equation (17b) results in Vθ̂t|t−1 = A2 Vθ̂t−1|t−1 + (I −
A2 )Vθ̂t−1|t−1 = Vθ̂t−1|t−1 . This ensures that the variance of the added random evolution would not
cause over-dispersion in the parameter estimation algorithm for all time.
The following proposition specifies an upper bound on the shrinkage factor. This upper bound
is calculated in the worst case, that is when the parameter is considered to be constant but the
modified evolution law (17a)-(17b) is used in the parameter estimation filter for estimating it.
Utilization of this upper bound in the kernel shrinkage algorithm ensures the boundedness of
the variance of the estimated parameters distribution that is obtained according to the PE-based
artificial evolution update law and the kernel smoothing augmented with the shrinkage factor.
Proposition 1: Upper bound on the kernel shrinkage factor: Given the parameter update law
(17a)-(17b), the estimated parameters conditional normal distribution based on the kernel smoothing as given by equation (21), results in an upper bound for A that is obtained as A ≤ I(1 −
r
2
σmin (Pmax
ΨVy ΨT V −1 )
θ̂

2
σmax (Pmax
ΨVy ΨT V −1 )

(j)

), where Ψ denotes the ψt

in equation (17a) but considered as a constant

θ̂

parameter between the time steps t and t − 1. Moreover, σmin and σmax denote the minimum and
the maximum eigenvalues of a matrix, respectively, W denotes the upper-bound on the variance
of the added noise Wt , Vy denotes the upper-bound on the variance of the measurement noise
Rt , Vθ̂ denotes the variance of the parameters when they are constantqthat can be assumed to be
the same as the initial covariance of the parameters, and Pmax = γ0 trace(Emax Emax T ), where
γ0 denotes the initial value of the step size, and Emax Emax T is a design parameter that denotes
the maximum acceptable variance among the prediction error vector elements.
Proof: The proof is provided in the appendix.
(j)
The convergence of the estimated parameter particles θ̂t−1|t−1 , j = 1, ..., N to the local mini¯ θ̂(j)
mum of E(J(
t−1|t−1 )|y1:t−1 , x̂t ) is now investigated in the following subsection. The developed
convergence proof does not ensure the convergence of the PE-based parameter estimation method
to the true parameter value, but only to a set of zeros of the gradient of the chosen performance
index. The conditions under which the convergence of the estimated parameters to their optimal
values can be guaranteed as N → ∞ is stated in Remark 1 (refer to the next page).
F. Convergence of the PE-based Parameter Update Law
In this subsection, it will be shown that the update law (17a)-(17b) can guarantee the conver(j)
gence of the parameter estimate particles θ̂t−1|t−1 , j = 1, ..., N (after the resampling step), to a
¯ θ̂(j)
local minimum of E(J(
t−1|t−1 )|y1:t−1 , x̂t ), that is located in a compact set of {xt , θt }, denoted
by DN as per Assumption A1.
In order to investigate the convergence of our proposed PE-based modified artificial evolution
law for updating the parameter particles distribution and to achieve a local minimization of
¯ t−1 )|y1:t−1 , xt ), consider equation (17a), where γt denotes a time-varying step size such
E(J(θ
March 13, 2018

DRAFT

17

that limt→∞ γt = µ0 > 0, where µ0 is a small positive constant. The introduction of the step
size γt is necessary to transform the discrete-time model (17a)-(17b) into a continuous-time
representation as shown subsequently.
First, we need to state the following two assumptions A2-A3 according to [30], to guarantee
the convergence of our proposed algorithm as presented in our main result below in Theorem
1. Specifically, we have:
Assumption A2. The function Q((t, θ̂t−1|t−1 )) is sufficiently smooth and twice continuously
differentiable w.r.t. , and |Q ((t, θ̂t−1|t−1 ))| ≤ C for θ̂t−1|t−1 ∈ DN , where Q ((t, θ̂t−1|t−1 ))
denotes the second derivative of Q((t, θ̂t−1|t−1 )) w.r.t. .
Assumption A3. The observation sequence yt (generated from equation (2)), is such that
¯ θ̂t−1|t−1 ) and Ē[ d Q((t, θ̂t−1|t−1 ))] = −g(θ̂t−1|t−1 ) exist for all
Ē(Q((t, θ̂t−1|t−1 ))) = J(
dθ̂t−1|t−1
P
θ̂t−1|t−1 ∈ DN , where Ē(Q((t, θ̂t−1|t−1 ))) = κ1 tτ =t−κ EQ((τ, θ̂t−1|t−1 )).
It must be noted that the kernel shrinkage method, as stated earlier, attempts to retain the mean
of the parameter estimation particles at time t near the estimated parameter in the previous time
step t−1, i.e. θ̂t−1|t−1 . Therefore, in the following theorem the convergence properties of θ̂t−1|t−1
is addressed. The main result of this section is stated below.
Theorem 1. Consider the parameter estimation algorithm as specified by the equations (17a)(19). Also consider the a posteriori parameter estimate as governed by equation (20). Let
(j)
Assumptions A1 to A3 hold. It now follows that the particles θ̂t−1|t−1 , j = 1, ..., N , and consequently the distribution of the estimated parameter particles approximated by the particle filter
(j)
(j)
¯ θ̂(j)
πθNt−1|t−1 (dθt−1 ), w.p.1 converge either to the set DC = {θ̂t−1|t−1 |θ̂t−1|t−1 ∈ DN , (j)d J(
t−1|t−1 )
dθ̂t−1|t−1

= 0, j = 1, ..., N } or to the boundary of DN as t → ∞.
Proof: The proof is provided in the appendix.
The main reason that our proposed dual state/parameter estimation method for its implementation does not necessarily need more particles than the one that we needed for only the state
estimation scheme, is illustrated by the result that one can extract and obtain from Theorem
1. According to this theorem, PE-based modified artificial evolution law enables each single
particle to tend to DC . Therefore, even increasing the number of particles would not affect the
convergence properties of the filter but it can certainly result in a more accurate state/parameter
estimates.
It was indicated earlier that the above theorem can only guarantee boundedness of the estimated
parameter distribution from the particle filters and not its convergence to the optimal distribution.
However, from practical considerations it is not readily feasible to obtain an exact dynamical
behavior that fully specifies variations of the system health parameters. This is due to the fact
that these parameters are affected by fault occurrences and/or damage during the operation
of the system. Consequently, the estimated parameters optimal values cannot be guaranteed
unless Assumptions A1, A2, and A3 are satisfied. However, based on the results of Theorem
1 and Proposition 1 one can ensure that the probability density function and its related kernel
Kθ (dθt |θt−1 , xt−1 ) (in the particle filter) do remain bounded. Moreover, the convergence of the
March 13, 2018

DRAFT

18

dual state and parameter estimation algorithm can be investigated according to the result of
Theorem 3.1 in [41] as explained in the following remark.
Remark 1. Using the extended setting introduced in [41], and assuming that ρ(yt |xt , θ̂t ) < ∞,
and Kx (xt |xt−1 , θ̂t−1 ) < ∞, the boundedness of the parameter estimation transition kernel
Kθ (θ̂t |θ̂t−1 , x̂t ) is ensured from the Proposition 1 and Theorem 1. Therefore, the convergence of
our proposed dual state/parameter estimation filter to their optimal distributions, for {xt , θt } ∈
DN as N → ∞ can be investigated along the lines of Theorem 3.1 of [41]. This is left as a
topic of our future research.
G. Equivalent Flop (EF) Complexity Analysis of Dual State/Parameter Estimation Algorithm
In this subsection, the computational complexity of our proposed dual state and parameter
estimation algorithm is studied based on the equivalent flop (EF) complexity metric. This measure
is introduced in [42], in which the number of flops that result in the same computational
time corresponding to a given operation is evaluated. The EF metric is mostly evaluated for
those operations that depend on matrix and vector manipulations. The details regarding the EF
complexity analysis for (a) our proposed dual estimation scheme, (b) the conventional Bayesian
method for state and parameter estimation [29] (where Regularized particle filter structure is
utilized to implement the filter), and (c) the recursive maximum likelihood (RML) method
according to simultaneous perturbation stochastic approximation (SPSA) algorithm [26] are
presented in the Appendix C.
Let nx , nθ . ny denote the dimension of the state vector, the parameter vector, and the
output vector, respectively, and c1 , c2 , c3 denote the complexity of random number generation,
resampling, and regularization, respectively. According to the summarized results in Tables VIII,
IX, X, and XI (in Appendix C), the EF complexity of the selected algorithms are given in Table
I, where only the dominant parts of CD (nx , nθ , c1 , c2 , c3 , N ) (this represents the EF complexity of
our proposed dual estimation scheme), CB (nx , nθ , c1 , c2 , c3 , N ) (this represents the EF complexity
of the conventional Bayesian method [29]), and CM (nx , nθ , c1 , c2 , c3 , N ) (this represents the
recursive maximum likelihood (RML) method according to SPSA algorithm [26]) are provided.
This selection is justified by the fact that N  1, therefore the dominant terms are related to
N.
To achieve the same complexity level in the three schemes, the number of the required particles
in our proposed dual estimation scheme is determined based on the number of the particles
that are utilized in the other two methods that are denoted by NB and NM . Namely, we have
N = NB (1 − [2n2θ + 5nθ + 2nθ ny + 6ny + 2nθ − 6nx nθ − c3 nθ ]/[3n2x + 5n2θ + 6nθ + 2nθ ny + 7ny +
3nx + c1 (nx + nθ ) + c2 (nx + nθ ) + c3 nx ]), and N = NM (1 − [n2x + 5n2θ + 2nθ + nx + 2nθ ny + 7ny +
c2 nθ ]/[3n2x + 5n2θ + 6nθ + 2nθ ny + 7ny + 3nx + c1 (nx + nθ ) + c2 (nx + nθ ) + c3 nx ]). Considering the
first equation we have 2n2θ +5nθ +2nθ ny +6ny +2nθ −6nx nθ −c3 nθ > 3n2x +5n2θ +6nθ +2nθ ny +
7ny +3nx +c1 (nx +nθ )+c2 (nx +nθ )+c3 nx ), since the complexity corresponding to regularization
(c3 ) is assumed to be the dominant complexity term in the nominator of the coefficient of NB .
March 13, 2018

DRAFT

19

TABLE I
T HE A PPROXIMATED T OTAL E QUIVALENT F LOP (EF) C OMPLEXITY OF THE S ELECTED Three Filters.
Prediction Method
Our Dual Estimation Scheme
Bayesian PF-based Estimation Method [29]
RML PF-based Estimation Method [26]

Total Equivalent Complexity
CB (nx , nθ , c1 , c2 , c3 , N ) ≈ N (3n2x + 5n2θ + 6nθ + 2nθ ny + 7ny + 3nx + c1 (nx + nθ ) + c2 (nx + nθ ) + c3 nx )
CB (nx , nθ , c1 , c2 , c3 , N ) ≈ N (3n2x + 3n2θ + 6nx nθ + (1 + c1 + c2 + c3 )nx + (1 + c1 + c2 + c3 )nθ + ny )
CM (nx , nθ , c1 , c2 , c3 , N ) ≈ N (2n2x + 4nθ + 2nx + c1 (2nx + nθ ) + c2 nx + c3 nx )

Therefore, the coefficient of NB is greater than that of N . Therefore, this indicates that in order
to achieve the same complexity by using the dual estimation algorithm and the one by using
the conventional Bayesian method using regularized particle filter structure, the number of the
particles required in the Bayesian method should be selected less than the number of particles in
the dual estimation method. On the other hand, similar analysis for the RML method according
to the second equation above shows that since n2x + 5n2θ + 2nθ + nx + 2nθ ny + 7ny + c2 nθ <
3n2x + 5n2θ + 6nθ + 2nθ ny + 7ny + 3nx + c1 (nx + nθ ) + c2 (nx + nθ ) + c3 nx ), unlike the Bayesian
method, in the RML method one needs more particles in order to achieve the same computational
complexity as that of the dual estimation scheme.
Finally, the EF complexity results are utilized to measure the time complexity of each algorithm
considering the fact that the EF complexity is proportional to the time complexity of the
algorithm. Our proposed dual state and parameter estimation scheme is an effective methodology
for solving the fault diagnosis of nonlinear systems problem. Without loss of any generality one
practically initiate operating the system from the healthy mode of operation. During the healthy
operation of the system our proposed dual state and parameter estimation strategy can provide
an accurate and reliable information on the health parameters of the system. This information
can then be readily used to perform the task of fault detection, isolation and identification of the
system. This is due to the fact that following the presence or injection of a fault in the components
of the system, the deviation and changes in the health parameters do provide clear signals for
the presence of a fault. In the next subsection, the application of our approach developed in
the previous subsections will be demonstrated and utilized for addressing the fault diagnosis
problem of nonlinear systems.
H. The Fault Diagnosis Formulation
Determination and diagnosis of drifts in unmeasurable parameters of a system require an online parameter estimation scheme. In parametric modeling of a system anomaly or drift, generally
it is assumed that the parameters are either constant or dependent on only the system states [43].
Hence, drifts in the parameters must be estimated through estimation techniques.
In [44], various parameter estimation techniques such as least squares, instrumental variables
and estimation via discrete-time models have been surveyed. The main drawbacks and limitation
with such methods arise due to the complexity and nonlinearity of the systems that we are
considering in this paper that render the parameter estimation here a nonlinear optimization

March 13, 2018

DRAFT

20

problem that must be solved in real-time. In [45] a nonlinear least squares (NLS) optimization
scheme is developed for only the fault identification of a hybrid system.
Parameter estimation techniques that are used for fault diagnosis of system components
generate residuals by comparing the estimated parameters that are obtained by either the ordinary
least squares (OLS) or the recursive least squares (RLS) algorithms with parameters that are
estimated under the initial fault free operation of the system [46].
The fault diagnosis problem under consideration deals with simultaneously obtaining an
optimal estimate of the system states as well as the time-varying parameters of a nonlinear
system whose dynamics is governed by the discrete-time stochastic model,
xt+1 = ft (xt , θtT λ(xt ), ωt ),

(22)

yt = ht (xt , θtT λ(xt )) + νt ,

(23)

where ft : Rnx × Rnθ × Rnω −→ Rnx is a known nonlinear function, θt ∈ Rnθ is the unknown
and time-varying parameter vector that for a healthy system is set to 1, λ : Rnx −→ Rnθ is a
differentiable function that determines the relationship between the system states and the health
parameters. The function ht : Rnx × Rnθ −→ Rny is a known nonlinear function, ωt and νt are
uncorrelated noise sequences with covariance matrices Lt and Vt , respectively. According to the
formulation used in equations (22) and (23), the parameter θt is a multiplicative fault vector
whose value is considered to be set equal to 1 under the healthy mode of the system operation.
The model (22) and (23) is now used to investigate the problem of fault diagnosis (FD),
which in this work is defined as the problem of fault detection, isolation, and identification
(FDII) when the system health parameters are considered to be affected by an unknown and
potentially time-varying multiplicative fault vector θt .
The system health parameters are known functions of the system states, λ(xt ), and the
multiplicative fault vector θt is to be estimated. In other words, the a posteriori estimated
parameter θ̂t|t will be used to generate residual signals for accomplishing the fault diagnosis
goal and objective. It is worth noting that based on our proposed formulation in equations (22)
and (23) in order to capture variations of the system health parameters, these parameters are
considered to be functions of the system states while their time variations due to degradations
and anomalies are captured by introducing the fault vector. Therefore, changes due to variations
in the system states are not considered as faults and determination of thresholds for developing
the fault diagnosis decision making scheme will effectively be based on the healthy system in
which the fault vector is represented and set to one.
The required residuals to determine the fault detection criterion are obtained as the difference
between the estimated parameters under the healthy operational mode that is denoted by θ̂0 , and
the estimated parameters under the faulty operational mode of the system that is denoted by θ̂t|t
as follows
rt = θ̂0 − θ̂t|t .
(24)
It should be pointed out that the true value of the parameter is denoted by θt? , which is assumed
to be unknown.
March 13, 2018

DRAFT

21

For the implementation of our proposed fault diagnosis strategy that is constructed based
on the previously developed state/parameter estimation framework, the parameter estimates will
be considered as the main indicators for detecting, isolating, and identifying the faults in the
system components. The residuals are generated from the parameter estimates under the healthy
and faulty operational modes of the system according to equation (24). The estimation of the
parameters under the healthy operational mode is determined according to,
θ̂0 = argmax(−log(p̂(θ0 |y1:T )),

(25)

where p̂(θ0 |y1:T ) denotes the probability density (conditioned on the observations up to time
T associated with the healthy data), that is obtained from the collected estimates and fitted
to a normal distribution. The time window T is chosen according to the convergence time of
the parameter estimation algorithm. The thresholds to indicate the confidence intervals for each
parameter are obtained through Monte Carlo analysis that is performed under different singlefault and multi-fault scenarios. The estimated parameters θ̂t|t are generated by following the
procedure that was developed and proposed in previous subsections.
IV. FAULT D IAGNOSIS OF A G AS T URBINE E NGINE
In this section, the utility of our proposed dual estimation framework when applied to the
problem of fault diagnosis of a nonlinear model of a gas turbine engine is demonstrated and
investigated. The performance of our proposed state/parameter estimation scheme will be evaluated and investigated when the gas turbine is subjected to deficiencies in its health parameters
due to injection of simultaneous faults.
A. Model Overview
The mathematical model of a gas turbine used in this paper is a single spool jet engine as
developed in [34]. The four engine states are the combustion chamber pressure and temperature, PCC and TCC , respectively, the spool speed S, and the nozzle outlet pressure PNLT . The
continuous-time state space model of the gas turbine is given as follows,
ṪCC =

1
[(cp TC θmC ṁC + ηCC Hu ṁf − cp TCC θmT ṁT )−
cv ṁcc

cv TCC (θmC ṁC + ṁf − θmT ṁT )],
Ṡ =

ηmech θmT ṁT cp (TCC − TT ) − θmC ṁC cp (TC − Td )
,
π 2
JS( 30
)

PCC 1
[(cp TC θmC ṁC + ηCC Hu ṁf − cp TCC θmT ṁT )−
TCC cv ṁcc
γRTCC
cv TCC (θmC ṁC + ṁf − θmT ṁT )] +
(θmC ṁC + ṁf − θmT ṁT ),
VCC
β
TM
(θmT ṁT +
θm ṁC − ṁNozzle ).
ṖNLT =
VM
β+1 C

(26)

ṖCC =

For the physical significance of the model parameters and details refer to [34]. The five gas
turbine measured outputs are considered to be the compressor temperature (y1 ), the combustion

March 13, 2018

DRAFT

22

chamber pressure (y2 ), the spool speed (y3 ), the nozzle outlet pressure (y4 ), and the turbine
temperature (y5 ), namely
y1 = TC = Tdiffuser [1 +

1
PCC γ−1
) γ − 1]],
[(
θηC ηC Pdiffuzer

y2 = PCC , y3 = S, y4 = PNLT ,
y5 = TCC [1 − θηT ηT (1 − (

(27)

PNLT γ−1
) γ ].
PCC

In order to discretize the above model for implementation of our proposed dual state/parameter
estimation particle filters, a simple Euler Backward method is applied with the sampling period
of Ts = 10 msec.
The system health parameters are represented by the compressor and the turbine efficiency,
ηC and ηT , respectively, and the compressor and turbine mass flow capacities, ṁC and ṁT ,
respectively. A fault vector is incorporated in the above model to manifest the effects of system
health parameters that are denoted by θ = [θηC , θmC , θηT , θmT ]T . By introducing a new parameter
as θ́ηC = θη1 , the measurement equations (27) can be represented as smooth functions with
C
respect to the fault parameters. Each parameter variation can be a manifestation of changes in
the fault vector that is considered as a multiplicative fault type. All the simulations that are
conducted in this section corresponds to the cruise flight condition mode.
In order to demonstrate the effectiveness and capabilities of our proposed algorithms, we have
also conducted simulation results corresponding to two other techniques in the literature, namely
(a) the conventional Bayesian method [29], and (b) the well-known recursive maximum likelihood
(RML) parameter estimation method based on PF [23], [24], [26]. It should be noted that the
number of particles in each algorithm is chosen based on the execution time of the algorithm
such that approximately the same execution time is achieved for each algorithm. Moreover,
a third method known as the gradient free PF-based RML method [26] could not yield an
acceptable performance in this application given the large number of the tuning parameters that
are required and associated with each parameter in this method. Therefore, the RML based on
the direct gradient method is utilized for the purpose of our performance comparison analysis.
(j)
(j)
In our schemes, the adaptive step size (Pt = γt Rt ) is defined as the product of the constant
(j)
γt (γt = 0.9) with Rt , which is evaluated on-line from the trace of the prediction error covariance
matrix that is estimated from the maximum likelihood method. On the other hand, the step size
in the RML method was chosen as γt = 0.05 = const. that is obtained by trial and error. The
residuals corresponding to the parameter estimates are also obtained. Based on the percentage of
the maximum absolute error criterion, a convergence time of 2 seconds is obtained for estimating
both the states and parameters corresponding to 25 Monte Carlo runs of simultaneous faults with
severities ranging from 1% to 10% loss of effectiveness of the healthy condition.
To choose the number of particles for implementation of the state and parameter estimation
filters, a quantitative study is now conducted. Specifically, based on the mean absolute error
(MAE%) that was obtained at the estimation process steady state and by taking into account the
March 13, 2018

DRAFT

23

algorithm’s computational time, the number of particles is chosen as N = 50 for both the state
and parameter estimation filters in this application. On the other hand, by taking into account the
average execution time of 18sec for one iteration of the dual estimation algorithm, an equivalent
execution time can be achieved for the Bayesian algorithm with NB = 45 and for the RML
algorithm with NM = 150 particles.
We have confirmed through simulations that acceptable performance and convergence times
are obtained. The shrinkage matrix is selected as 0.93I. The initial distributions (i.e., the mean
and covariance matrices) of the states and parameters are selected to correspond to the cruise
flight operational condition as provided in [34]. In what follows, the two main case scenarios
for performing fault diagnosis of the gas turbine engine are presented.
Scenario I: Concurrent Faults in the Compressor and Turbine Health Parameters.
In this scenario, the input fuel flow rate to the engine is changed by decreasing it by 2%
from its nominal value one second after reaching the steady state condition. Next, the effects
of concurrent faults in both the compressor and the turbine are studied by injecting sequential
fault patterns affecting the system components. Specifically, at time t = 4 sec the compressor
efficiency is reduced by 5% (this represents the level of the fault severity), followed by at
t = 9 sec the same fault type affecting the compressor mass flow capacity, and at t = 14 sec the
same fault type affecting the turbine efficiency, and finally at t = 19 sec the same fault type is
applied to the turbine mass flow capacity.
The results corresponding to changes in the fault parameters are depicted in Figure 2. The
dotted lines depict the confidence bounds for residuals that are determined based on 50 Monte
Carlo simulation runs under various concurrent and simultaneous single and/or multiple fault
scenarios by using the PE-based method. By analyzing the residuals, the detection time of a
fault in each component and its severity can be determined and identified. It follows from this
figure that the residuals corresponding to the dual estimation method using PE-based method as
well as the RML method practically do not exceed their confidence bounds subject to changes
in the engine input (which is applied at t = 1 sec). On the other hand, the Bayesian method
shows false alarms in the residuals corresponding to the turbine heath parameters. Moreover,
this method is also not capable of tracking the changes in the fault vector in the selected time
window in all the residual signals after fault occurrence.
In order to obtain a quantitative measure on the accuracy and reliability of our proposed estimation algorithm the results related to the 5% fault severity using the mean absolute error (MAE%)
metric of estimates corresponding to the last 2 sec of simulations (following the algorithm
convergence) after each change are provided and summarized in Table II. The state/parameter
estimation MAE% for (a) our proposed dual estimation algorithm according to the PE-based
method with N = 50, (b) the conventional Bayesian method based kernel smoothing (KS-based)
with NB = 45, and (c) the RML method with NM = 150 are presented. In this table, the i-th
fault for i = 1, ..., 4 denotes the last 2 sec of simulations after the i-th fault occurrence, and the
first column refers to the healthy system before the fault occurrence.
The results shown in Table II demonstrate that for the PE-based method the maximum MAE%
March 13, 2018

DRAFT

24

Residual for the Compressor Efficiency Fault

Residual for the Compressor Mass Flow Fault
1.02

1.01
1.01
1
1
0.99

0.99

0.98

0.97

0.96

0.96

0.95

0.95

0.94
0.93

Residual upper bound
Residual lower bound
PE-based method
RML-based method
Bayesian KS-based method

0.98

0.97

0.94
0

5

10

15

20

0.93

0

5

Residual for the Turbine Efficiency Fault

10

15

20

Residual for the Turbine Mass Flow Fault

1.02
1.04

1.01
1

1.02

0.99
1
0.98
0.98

0.97
0.96

0.96

0.95
0.94

0.94
0.93

0

5

10

15

20

0.92

0

5

10

15

20

Fig. 2. Residuals corresponding to the concurrent fault scenarios in the turbine and the compressor parameters.

for the states is between 0.03%−1.06% of their nominal values. In case of the health parameters,
for ηC and ṁC the maximum MAE% is around 0.91% and for ηT and ṁT it is around 0.98% of
their nominal values. On the other hand, according to the results presented in Table II (b), the
maximum MAE% for the states corresponding to the RML method is between 0.1% − 1.16%
of their nominal values. In case of the health parameters, the maximum MAE% ranges between
0.8% − 2.8% of their nominal values, where both mass flow rates are estimated with higher
MAE%. The results corresponding to Table II (c) indicate that the maximum MAE% for the
state estimation results in the Bayesian KS-based method ranges between 0.2% − 18.3%, for the
compressor health parameters between 0.53% − 19.3%, and for the turbine health parameters
between 0.15% − 3.0%.
The MAE% for the estimated measurements (outputs) of the engine are also provided in Table
III for the PE-based, the RML, and the Bayesian KS-based methods. From the results presented
in Table III (a) one can conclude that the maximum MAE% for the temperature (of the turbine
and the compressor) corresponding to our proposed PE-based method is less than 0.3%, and for
the spool speed it is less than 0.16%, and for the compressor pressure it is less than 1.4%, and
for the turbine pressure it is less than 2.5%. On the other hand, the results presented in Table
III (b) for the RML method show that the maximum MAE% for the compressor and turbine
temperatures is less than 0.4% and 0.6%, respectively. For the spool speed the MAE% is less

March 13, 2018

DRAFT

25

than 0.2%, and for the compressor pressure it is less than 1.5% and for the turbine pressure it is
less than 2.5%. For the Bayesian KS-based method, instead of the compressor temperature and
spool speed outputs, the maximum MAE% exceeds 13% of the nominal values. Consequently,
the results presented in the above two tables confirm that the Bayesian KS-based method does
not have acceptable estimation accuracy as compared to the other two alternative methods. On
the other hand, the PE-based method outperforms the RML method significantly. The accuracy
in the measurement estimation is an important consideration and factor given that from practical
perspectives the system states and parameters are unknown. Therefore, it is generally necessary
to judge the estimation accuracy based on the output estimation error performance and behavior.
In order to demonstrate and illustrate the accuracy of our proposed fault detection algorithm
based on the dual state/parameter estimation scheme, at the end of this section a quantitative
study is conducted by performing a confusion matrix analysis [47] in presence of various fault
scenarios having different fault severities and in presence of the same level of process and
measurement noise that are stated in [34] for the PE-based method, the RML method, and the
Bayesian KS-based combined state and parameter estimation algorithm.
Scenario II: Simultaneous Faults in the Compressor and the Turbine Health Parameters.
In the second scenario, a simultaneous fault in all the 4 health parameters of the engine is
applied at t = 9 sec. The compressor and the turbine efficiencies faults follow the pattern of
a drift fault that starts at t = 9 sec and causes a 5% loss of effectiveness in the compressor
efficiency by the end of the simulation time (i.e. at t = 19 sec), and a 3% loss of effectiveness
in the turbine efficiency by the end of the simulation time. Simultaneously, the mass flow rate
capacities of both the compressor and the turbine are affected by a 5% loss of effectiveness fault.
The residuals corresponding to the three previous estimation methods are provided in Figure 3.
The simulations show that in case of changes in the engine input (applied at t = 1 sec) the RML
method residuals has high false alarm rates as compared to dual estimation method according
to PE-based algorithm, similar to the first scenario for the concurrent faults. More quantitative
analysis on the performance of the RML method that is compared to the PE-based method is
provided in the subsequent subsection. The presented results admit that the Bayesian KS-based
method is not able to track the variations in the fault vectors in the case of simultaneous fault
scenario.
The results in Table IV (a) show that for our proposed PE-based method, the maximum
MAE% for both state and parameter estimates are between 0.1% − 0.5% of their nominal values.
However, in the worst case the post fault estimated MAE% of the ṁT is 0.47% of its nominal
value. Moreover, in the results shown for the RML method in Table IV (b) it follows clearly
that the state estimation MAE% can be achieved within 0.1% − 0.8% of the nominal values,
whereas the parameter estimation MAE% is achieved within 0.7% − 3% of the nominal values
with higher error rates after the fault occurrence, specially in the compressor and turbine mass
flow rate capacities. However, for the Bayesian KS-based method in Table IV (c) the maximum
MAE% is achieved within 0.19% − 8.4% of the nominal values for the estimated states and
within 0.25% − 7.2% of the nominal values for the estimated parameters.
March 13, 2018

DRAFT

26

TABLE II
S TATE /PARAMETER MAE% IN CASE OF CONCURRENT FAULT SCENARIOS FOR ( A ) T HE D UAL ESTIMATION ALGORITHM
ACCORDING TO PE- BASED METHOD WITH N = 50, ( B ) T HE RML METHOD WITH N = 150, AND ( C ) T HE BAYESIAN
KS- BASED METHOD WITH N = 45.
State

No Fault

PCC
N
TCC
PNLT
ηC
ṁC
ηT
ṁT

0.3529
0.1473
0.2683
0.8575
0.2702
0.6621
0.2865
0.4744

(a)
1nd Fault 2rd Fault
0.2097
0.0761
0.1674
0.5325
0.1785
0.4229
0.1648
0.4557

0.3614
0.1087
0.1678
0.3978
0.2749
0.3682
0.1743
0.4889

(b)
1nd Fault 2rd Fault

State

No Fault

PCC
N
TCC
PNLT
ηC
ṁC
ηT
ṁT

0.5352
0.0995
0.2064
0.7181
0.9268
1.6338
0.9252
1.3719

0.6342
0.0912
0.2443
0.8112
1.9195
1.8037
0.7876
1.7858

State

No Fault

1nd Fault

PCC
N
TCC
PNLT
ηC
ṁC
ηT
ṁT

1.8961
0.2127
0.4789
0.6841
0.7248
0.5306
0.1445
1.4943

2.6032
0.5032
1.0029
1.3838
3.3660
1.3399
0.7394
1.6979

0.3921
0.1018
0.2574
0.7771
2.1698
1.0761
0.8714
1.6162
(c)
2rd Fault
6.1590
0.4490
1.6025
3.6558
3.0788
4.3086
0.9198
2.6633

3th Fault

4th Fault

0.4336
0.1624
0.3838
1.0614
0.3879
0.9132
0.4885
0.9757

0.2374
0.0296
0.1155
0.3213
0.2107
0.2236
0.1873
0.5037

3th Fault

4th Fault

0.8934
0.2060
0.5174
1.1603
1.4508
2.6062
1.6517
1.9666

0.5882
0.1383
0.3374
0.5666
1.4651
2.2717
1.1411
2.7653

3th Fault

4th Fault

18.2816
4.7275
8.6930
14.9288
6.3141
19.3026
1.0082
3.0198

7.9636
2.5564
6.5029
8.9250
4.9584
11.8476
1.8379
1.7125

The MAE% measurement (output) estimate error given in Table V (a) for the PE-based method
shows that after simultaneous fault occurrences the error increases when compared to their
values before the fault occurrences. This is caused due to accumulation of parameter estimation
errors while all the four parameters are affected by a fault. On the other hand, the results
corresponding to the output estimation as given in Tables V(a)-(c) show that with the exception
of the turbine pressure, our PE-based method outperforms the RML method for estimating the
other four measurement outputs. However, the maximum MAE% for the outputs from Bayesian
KS-based method performs high level of errors after fault occurrence in all measurement outputs
as compared to the other two estimation methods.

March 13, 2018

DRAFT

27

TABLE III
O UTPUT ESTIMATION MAE% IN CASE OF CONCURRENT FAULT SCENARIOS FOR ( A ) T HE D UAL ESTIMATION ALGORITHM
ACCORDING TO PE- BASED METHOD WITH N = 50, ( B ) T HE RML METHOD WITH N = 150, AND ( C ) T HE BAYESIAN
KS- BASED METHOD WITH N = 45.
Output

No Fault

TC
PC
N
TT
PT

0.2893
1.3548
0.1473
0.2034
2.2231

(a)
1nd Fault 2rd Fault
0.2319
1.2332
0.0761
0.1857
2.1696

0.2749
1.2507
0.1087
0.1911
2.1839

(b)
1nd Fault 2rd Fault

Output

No Fault

TC
PC
N
TT
PT

0.2902
1.4012
0.0995
0.2181
2.3446

0.3240
1.3755
0.0912
0.2461
2.3994

Output

No Fault

1nd Fault

TC
PC
N
TT
PT

0.9416
1.9784
0.2127
0.3428
2.2715

0.7569
2.6211
0.5032
0.7125
2.7023

0.2956
1.3030
0.1018
0.2122
2.1356
(c)
2rd Fault
0.6033
6.5484
0.4490
2.6660
3.9467

3th Fault

4th Fault

0.2805
1.4070
0.1624
0.2804
2.4577

0.2357
1.1813
0.0296
0.1322
2.0783

3th Fault

4th Fault

0.3985
1.4779
0.2060
0.5786
2.5220

0.2991
1.2902
0.1383
0.5206
2.2474

3th Fault

4th Fault

3.4434
18.0905
4.7275
13.3328
14.9107

0.3789
7.8644
2.5564
8.6219
8.7964

To summarize, our proposed PE-based fault diagnosis algorithm is capable of detecting,
isolating and estimating the component faults of a gas turbine engine with an average accuracy of
0.3% for the compressor and 0.5% for the turbine faults. In contrast the RML algorithm is capable
of achieving the performance of an average 3% for the compressor and 1.6% for the turbine
faults. The Bayesian KS-based method does not have acceptable accuracy for simultaneous fault
diagnosis application.
Fault Diagnosis Confusion Matrix Analysis
Finally, in this subsection a quantitative study is performed by utilizing the confusion matrix
analysis [47] to evaluate the increase in the false alarms and/or misclassification rates of the faults
in our considered application when the fault diagnosis algorithm is implemented by our proposed
PE-based method with N = 50 particles, the RML method with N = 150, and the Bayesian
KS-based method with N = 45 particles. The thresholds corresponding to each algorithm are
determined from 25 Monte Carlo runs on simultaneous fault scenarios that are not necessarily the
same for the three algorithms. The confusion matrix data is obtained by performing simulations
for another 35 Monte Carlo simultaneous fault scenarios having different fault severities and in
March 13, 2018

DRAFT

28

Residual for the Compressor Efficiency Fault

Residual for the Compressor Mass Flow Fault
1.02

1.01
1.01
1

1

0.99

0.99

0.98

0.98

0.97

0.96

0.95

0.95

0.94
0.93

Residual upper bound
Residual lower bound
PE-based method
RML-based method
Bayesian KS-based method

0.97

0.96

0.94
0

2

4

6

8

10

12

14

16

0.93

18

0

2

4

Residual for the Turbine Efficiency Fault

6

8

10

12

14

16

18

14

16

18

Residual for the Turbine Mass Flow Fault

1.02
1.04

1.01
1

1.02

0.99
1
0.98
0.98

0.97
0.96

0.96

0.95
0.94

0.94
0.93

0

2

4

6

8

10

12

14

16

0.92

18

0

2

4

6

8

10

12

Fig. 3. Residuals corresponding to the simultaneous fault scenarios.

TABLE IV
S TATE /PARAMETER MAE% IN CASE OF SIMULTANEOUS FAULT SCENARIOS FOR ( A ) T HE PE- BASED METHOD WITH
N = 50, ( B ) THE RML METHOD WITH N = 150, AND ( C ) THE BAYESIAN KS- BASED METHOD WITH N = 45.

March 13, 2018

State

(a)
Before Fault

After Fault

State

(b)
Before Fault

After Fault

PCC
N
TCC
PNLT
ηC
ṁC
ηT
ṁT

0.2217
0.0535
0.2086
0.3970
0.1735
0.2811
0.1016
0.4589

0.2372
0.1061
0.1928
0.4291
0.1821
0.3293
0.1485
0.4744

PCC
N
TCC
PNLT
ηC
ṁC
ηT
ṁT

0.4865
0.1025
0.2247
0.7540
0.9295
1.7291
0.7306
1.4290

0.4383
0.1053
0.1888
0.5264
1.7956
2.9729
1.0040
1.5923

State

(c)
Before Fault

After Fault

PCC
N
TCC
PNLT
ηC
ṁC
ηT
ṁT

0.5186
0.1906
0.2834
0.6506
0.5992
0.5775
0.2463
0.2861

8.3590
2.7197
4.5219
5.0157
4.8358
7.2215
4.8958
4.3080

DRAFT

29

TABLE V
O UTPUT ESTIMATION MAE% IN CASE OF SIMULTANEOUS FAULT SCENARIOS FOR ( A ) T HE PE- BASED METHOD WITH
N = 50, ( B ) T HE RML METHOD WITH N = 150, AND ( C ) T HE BAYESIAN KS- BASED METHOD WITH N = 45.
(a)

(b)

Output

Before Fault

After Fault

Output

Before Fault

After Fault

TC
PC
N
TT
PT

0.2207
1.2926
0.0535
0.1565
2.1210

0.2852
1.3729
0.1027
0.1650
2.2348

TC
PC
N
TT
PT

0.2482
1.4051
0.1025
0.2058
2.1413

0.2832
1.4369
0.1053
0.1774
2.1580

(c)

Output

Before Fault

After Fault

TC
PC
N
TT
PT

0.3250
1.3621
0.1906
0.4126
2.2252

3.5756
8.1037
2.7197
4.9311
5.0404

presence of the same process and measurement noise covariances corresponding to 50% of the
nominal values of the process and measurement noise covariances (according to [34]). In these
scenarios, at each time more than one of the system health parameters are affected by component
faults.
The results are shown in Tables VI(c)-VI(a) corresponding to PE-based method with N = 50
particles, the RML method with N = 150, the RML method with N = 150, and the Bayesian
KS-based method with N = 45 particles, respectively. In these tables the rows depict the
actual number of fault categories that are applied and the columns represent the number of
estimated fault categories. The diagonal elements represent the true positive rate (T P ) for each
fault occurrence. The accuracy (AC), precision (P ), and the false positive rate (F P ) of the
three algorithms are also evaluated from the confusion matrix results according to the following
formulae [47],
P5
P4
c
c
jj
jj
j=1
j=1 c5j
AC = P5 P5
, P j = P5
, F P = P5
,
i=1
j=1 cij
i=1 cij
j=1 c5j
where cij , i, j = 1, ..., 5 denote the elements of the confusion matrix. In Table VII, the confusion
matrix results according to the above metrics for the Tables VI(c)-VI(a) are provided. The results
demonstrate that the accuracy of the fault diagnosis for the dual PE-based estimation algorithm
outperforms RML method with 7.43% and the false positive rate of 5.71% less than RML method.
The precision of the algorithm for all the system four health parameters is more than the ones
from RML method. However, the Bayesian KS-based method indicates poor accuracy and high
false alarm rate for the fault diagnosis of the system. Consequently, the PE-based method with
N = 50 outperforms the other two methods significantly in terms of higher accuracy, lower
March 13, 2018

DRAFT

30

TABLE VI
C ONFUSION MATRIX FOR ( A ) T HE PE- BASED METHOD WITH N = 50, ( B ) T HE RML METHOD WITH N = 150, AND ( C )
T HE BAYESIAN KS- BASED METHOD WITH N = 45.
(a)

(b)

Fault

η̇C

ṁC

η̇T

ṁT

No Fault

Fault

η̇C

ṁC

η̇T

ṁT

No Fault

ηC
ṁC
ηT
ṁT
No Fault

31
0
1
1
0

0
30
1
1
0

2
2
28
3
1

2
3
4
29
1

0
0
1
1
33

ηC
ṁC
ηT
ṁT
No Fault

28
1
2
1
0

2
27
3
3
2

3
1
26
4
1

2
4
3
26
1

0
2
1
1
31

(c)

Fault

η̇C

ṁC

η̇T

ṁT

No Fault

ηC
ṁC
ηT
ṁT
No Fault

10
9
6
5
10

5
13
6
7
9

6
8
9
8
7

4
6
7
11
4

10
9
7
4
5

false positive rate, and higher precision for all the four health parameters of the gas turbine
engine.
TABLE VII
C ONFUSION MATRIX A NALYSIS RESULTS .

Noise Level

AC%

FP%

Pη C %

PṁC %

PηT %

PṁT %

PE-based (50 Part.)
RML ( 150 Part.)
Bayesian KS-based (45 Part.)

86.29
78.86
25.95

5.71
11.43
85.71

93.94
87.50
25.00

93.75
72.97
32.50

77.78
74.29
23.68

74.36
72.22
34.38

V. C ONCLUSION
In this paper, a novel dual estimation filtering scheme is proposed and developed based on
particle filters (PF) to estimate a nonlinear stochastic system states and time variations in its
parameters. The dual structure is based on the extension of the Bayesian parameter estimation
framework. A dual structure is proposed for achieving simultaneous state and parameter estimation objectives. Performance results of the application of our method to a gas turbine engine under
healthy and faulty scenarios are provided to demonstrate and illustrate the superior capability and
performance of our scheme for a challenging fault diagnostic application as compared to the wellknown recursive maximum likelihood (RML) method based on particle filters and conventional
Bayesian method for combined state and parameter estimation based on particle filters while the
computational complexity of all the algorithms remains the same. On the other hand, the false
alarm rates of our proposed dual algorithm is significantly lower than the RML and conventional
March 13, 2018

DRAFT

31

Bayesian methods. These two main characteristics justify and substantiate the observation that
our proposed algorithm is more suitable for the purpose of fault diagnosis of critical nonlinear
systems that require lower fault detection times and false alarm rates. Moreover, the estimation
results accuracy in terms of the fault identification are also provided. The obtained results are
demonstrated and validated by performing a confusion matrix analysis. The verification and
validation of our results to real gas turbine engine is a topic of our future research.
R EFERENCES
[1] J. Korbicz, Z. Fathi, and W. Fred Ramirez, “State estimation schemes for fault detection and diagnosis in dynamic systems,”
International Journal of Systems Science, vol. 24, no. 5, pp. 985–1000, 1993.
[2] H. Kushner and H. Joseph, Probability methods for approximations in stochastic control and for elliptic equations, vol. 129.
Academic Press New York, 1977.
[3] V. Beneš, “Exact finite-dimensional filters for certain diffusions with nonlinear drift,” Stochastics: An International Journal
of Probability and Stochastic Processes, vol. 5, pp. 65–92, 1981.
[4] A. Doucet, N. De Freitas, and N. Gordon, Sequential Monte Carlo Methods in Practice. Springer Verlag, 2001.
[5] A. Bensoussan, R. Glowinski, and A. Raşcanu, “Approximation of some stochastic differential equations by the splitting
up method,” Applied Mathematics & Optimization, vol. 25, pp. 81–106, 1992.
[6] S. Lototsky and B. Rozovskii, “Recursive multiple wiener integral expansion for nonlinear filtering of diffusion processes,”
Stochastic Processes and Functional Analysis, Lecture Notes in Pure and Appl. Math, vol. 186, pp. 199–208, 1997.
[7] D. Crisan, J. Gaines, and T. Lyons, “Convergence of a branching particle method to the solution of the Zakai equation,”
SIAM Journal on Applied Mathematics, vol. 58, pp. 1568–1590, 1998.
[8] I. Arasaratnam and S. Haykin, “Cubature Kalman filters,” IEEE Transactions on Automatic Control, vol. 54, no. 6, pp. 1254–
1269, 2009.
[9] V. Kadirkamanathan, P. Li, M. Jaward, and S. Fabri, “A sequential Monte Carlo filtering approach to fault detection and
isolation in nonlinear systems,” in Proceedings of the 39th IEEE Conference on Decision and Control, vol. 5, pp. 4341–
4346, IEEE, 2000.
[10] M. E. Orchard, A Particle Filtering-based Framework for On-line Fault Diagnosis and Failure Prognosis. PhD thesis,
Georgia Institute of Technology, 2006.
[11] P. Li and V. Kadirkamanathan, “Particle filtering based likelihood ratio approach to fault diagnosis in nonlinear stochastic
systems,” IEEE Transactions on Systems, Man and Cybernetics. Part C, Applications and reviews, vol. 31, no. 3, pp. 337–
343, 2001.
[12] M. Zajac, “Online fault detection of a mobile robot with a parallelized particle filter,” Neurocomputing, vol. 126, pp. 151–
165, 2014.
[13] B. Zhao, R. Skjetne, M. Blanke, and F. Dukan, “Particle filter for fault diagnosis and robust navigation of underwater
robot,” IEEE Transactions on Control Systems Technology, vol. 22, no. 6, pp. 2399–2407, 2014.
[14] Q. Zhang, F. Campillo, F. Cerou, and F. Legland, “Nonlinear system fault detection and isolation based on bootstrap particle
filters,” in 44th IEEE Conference on Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC’05,
pp. 3821–3826, IEEE, 2005.
[15] Z.-H. Duan, Z.-x. Cai, and Y. Jin-xia, “Fault diagnosis and fault tolerant control for wheeled mobile robots under unknown
environments: A survey,” in Proceedings of the 2005 IEEE International Conference on Robotics and Automation, ICRA
2005., pp. 3428–3433, IEEE, 2005.
[16] V. Verma, Tractable particle filters for robot fault diagnosis. PhD thesis, Stanford University, 2004.
[17] F. Hutter, R. Dearden, et al., “The Gaussian particle filter for diagnosis of non-linear systems,” in Proceedings of the 5th
IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes, 2003.
[18] S. Tafazoli and X. Sun, “Hybrid system state tracking and fault detection using particle filters,” IEEE Transactions on
Control Systems Technology, vol. 14, no. 6, pp. 1078–1087, 2006.
[19] A. Doucet and V. B. Tadić, “Parameter estimation in general state-space models using particle methods,” Annals of the
institute of Statistical Mathematics, vol. 55, no. 2, pp. 409–422, 2003.
[20] M. West, “Mixture models, Monte Carlo, Bayesian updating, and dynamic models,” Computing Science and Statistics,
vol. 24, pp. 325–325, 1993.

March 13, 2018

DRAFT

32

[21] T. Flury and N. Shephard, “Bayesian inference based only on simulated likelihood: particle filter analysis of dynamic
economic models,” Econometric Theory, vol. 27, no. 05, pp. 933–956, 2011.
[22] N. Kantas, A. Doucet, S. S. Singh, and J. M. Maciejowski, “An overview of sequential Monte Carlo methods for
parameter estimation in general state-space models,” in 15th IFAC Symposium on System Identification (SYSID), Saint-Malo,
France.(invited paper), vol. 102, p. 117, 2009.
[23] G. Poyiadjis, A. Doucet, and S. S. Singh, “Particle methods for optimal filter derivative: Application to parameter
estimation,” in Proceedings of IEEE International Conference on Acoustics, Speech, and Signal Processing, (ICASSP’05),
vol. 5, pp. v–925, IEEE, 2005.
[24] G. Poyiadjis, A. Doucet, and S. S. Singh, “Particle approximations of the score and observed information matrix in state
space models with application to parameter estimation,” Biometrika, vol. 98, no. 1, pp. 65–80, 2011.
[25] T. B. Schön, A. Wills, and B. Ninness, “System identification of nonlinear state-space models,” Automatica, vol. 47, no. 1,
pp. 39–49, 2011.
[26] G. Poyiadjis, S. S. Singh, and A. Doucet, “Gradient-free maximum likelihood parameter estimation with particle filters,”
in American Control Conference, 2006, pp. 6–pp, 2006.
[27] J. Westerborn and J. Olsson, “Efficient particle-based online smoothing in general hidden Markov models,” in IEEE
International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 8003–8007, 2014.
[28] A. Doucet, M. Pitt, G. Deligiannidis, and R. Kohn, “Efficient implementation of Markov chain Monte Carlo when using
an unbiased likelihood estimator,” arXiv preprint arXiv:1210.1871, 2012.
[29] J. Liu and M. West, Combined parameter and state estimation in simulation-based filtering. Institute of Statistics and
Decision Sciences, Duke University, 1999.
[30] L. Ljung and T. Söderström, Theory and Pactice of Recursive Identification. 1983.
[31] I. Collings, V. Krishnamurthy, and J. Moore, “On-line identification of hidden Markov models via recursive prediction
error techniques,” IEEE Transactions on Signal Processing, vol. 42, no. 12, pp. 3535–3539, 1994.
[32] R. Isermann, Fault-diagnosis systems. Springer, 2006.
[33] E. Sobhani-Tehrani and K. Khorasani, Fault diagnosis of nonlinear systems using a hybrid approach, vol. 383. Springer,
2009.
[34] N. Meskin, E. Naderi, and K. Khorasani, “A multiple model-based approach for fault diagnosis of jet engines,” IEEE
Transactions on Control Systems Technology, vol. 21, no. 1, pp. 254–262, 2013.
[35] E. Naderi, N. Meskin, and K. Khorasani, “Nonlinear fault diagnosis of jet engines by using a multiple model-based
approach,” Journal of Engineering for Gas Turbines and Power, vol. 134, no. 1, p. 011602, 2012.
[36] S. S. Haykin, Kalman filtering and neural networks. Wiley Online Library, 2001.
[37] F. L. Lewis and F. Lewis, Optimal estimation: with an introduction to stochastic control theory. Wiley New York et al.,
1986.
[38] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian
Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188, 2002.
[39] C. Musso, N. Oudjane, and F. LeGland, “Improving regularised particle filters,” Sequential Monte Carlo methods in
practice, pp. 247–271, 2001.
[40] A. Budhiraja, L. Chen, and C. Lee, “A survey of numerical methods for nonlinear filtering problems,” Physica D: Nonlinear
Phenomena, vol. 230, no. 1, pp. 27–36, 2007.
[41] X.-L. Hu, T. Schon, and L. Ljung, “A general convergence result for particle filtering,” IEEE Transactions on Signal
Processing, vol. 59, no. 7, pp. 3424–3429, 2011.
[42] R. Karlsson, T. Schön, and F. Gustafsson, “Complexity analysis of the marginalized particle filter,” IEEE Transactions on
Singnal Processing, vol. 53, no. 11, pp. 4408–4411, 2005.
[43] V. Venkatasubramanian, R. Rengaswamy, K. Yin, and S. N. Kavuri, “A review of process fault detection and diagnosis:
Part I: Quantitative model-based methods,” Computers & Chemical Engineering, vol. 27, no. 3, pp. 293–311, 2003.
[44] R. Isermann, “Process fault detection based on modeling and estimation methods a survey,” Automatica, vol. 20, no. 4,
pp. 387–404, 1984.
[45] D. Wang, M. Yu, C. B. Low, and S. Arogeti, Model-based Health Monitoring of Hybrid Systems. Springer, 2013.
[46] S. Simani, Model-based fault diagnosis in dynamic systems using identification techniques. PhD thesis, dell’Universit di
Ferrara, 2003.
[47] R. Kohavi and F. Provost, “Confusion matrix,” Machine Learning, vol. 30, no. 2-3, pp. 271–274, 1998.

March 13, 2018

DRAFT

33

A PPENDIX
A: Proof of Proposition 1:Let us consider the modified artificial evolution law by assuming
(j)
that A = I in equation (17b). Let us substitute mt from equation (17a) where the superscript
(j) is omitted in order to define the modified artificial evolution law in a more general form that
is also applicable to each single particle as
θ̃t|t = θ̂t−1|t−1 + Pt ψt (t, θ̂t−1|t−1 ) + ζt ,

(28)

where Pt = γt Rt . Now, let V (.|y1:t ) denote the variance of the stochastic process assuming that
the observations up to time t are available, and C(., .|y1:t ) denotes the covariance of the two
stochastic processes by assuming that the observations up to time t are available. By taking into
account the relationship between the variance of both sides of equation (28) when the covariance
matrix is assumed to be non-singular and when the prediction error at time t is uncorrelated
with the parameter estimate at time t, given that θ̂t−1|t−1 is independent of Pt Ψt , therefore we
2
T
get V (θ̃t|t |y1:t ) = V (θ̂t−1|t−1 |y1:t ) + Pp
t ΨVy Ψ + Wt + 2C(θ̂t−1|t−1 , ζt |y1:t ) + 2C(Pt Ψt , ζt |y1:t ).
Furthermore, since Pt = γt Rt = γt trace(Et EtT ) is a scalar, one can write V (Pt Ψt |y1:t ) =
(E[Pt |y1:t ])2 V (Ψt |y1:t ) + (E[Ψt |y1:t ])2 V (Pt |y1:t ) + V (Pt |y1:t )V (Ψt |y1:t ) = Pt2 ΨVt ΨT .
In order to ensure that there is no information loss (particularly in the case that θt is constant),
one must have, V (θ̃t|t |y1:t ) = V (θ̂t−1|t−1 |y1:t ) = Vθ̂t−1|t−1 , which implies that, C(θ̂t−1 , ζt |Yt ) +
C(Pt Ψt , ζt |Yt ) = − 21 Wt − 12 Pt2 ΨVt ΨT . Therefore, negative correlations are needed to remove
the effects of unwanted information loss. In case of approximate joint normality of the stochastic
process (θ̂t−1|t−1 , ζt |Yt ) and (Pt Ψt , ζt |Yt ), the conditional normal evolution is obtained as
p(θ̂t|t |θ̂t−1|t−1 ) ∼ N (θ̂t|t |At θ̃t|t + (I − At )θ̂t−1|t−1 , (I − A2t )Vθ̂t−1|t−1 ),

(29)

where the mean of this Gaussian distribution at each time step is found from equation (17a),
when θ̂t−1|t−1 is substituted by its modified version according to the shrinkage kernel. The
+ Pt2 ΨVy ΨT Vθ−1
)]. Assuming
shrinkage matrix At , is obtained as At = I − [ 12 (Wt Vθ−1
t−1|t−1
t−1|t−1
that in the kernel shrinkage method, the variance of the evolution noise is interpreted as Wt =
(I − A2t )Vθ̂t−1|t−1 , therefore by replacing Wt in the equation that was obtained for At results in
1
At = I − [(I − A2t ) + Pt2 ΨVt ΨT Vθ−1
].
(30)
t−1|t−1
2
Let us assume that our main goal is to obtain and determine the shrinkage matrix At as
A = aI, therefore the matrix equation (30) can be written as
(a2 − 2a + 1)I − Pt2 ΨVt ΨT Vθ−1
= 0.
t−1|t−1

(31)

We are interested in obtaining an upper bound for the shrinkage matrix that can be used for
all time. Assuming that the last term in the right hand side of equation
p (30) has an upper bound
2
T −1
2
T −1
T ) with γ
given by |Pt ΨVt Ψ Vθt−1|t−1 | ≤ Pmax ΨVy Ψ Vθ̂ , where Pmax = γ0 trace(Emax Emax
0
denoting the initial step size, therefore Emax is considered as the maximum acceptable variance
of the prediction error, Vy is an upper bound of the measurement noise covariance, and Vθ̂ is
the minimum bound of the parameter estimation covariance that is considered to be similar to
March 13, 2018

DRAFT

34

the initial covariance of the parameters (before adding the evolution
a
rnoise in time). Therefore,
2
ΨVy ΨT V −1 )
σmin (Pmax
bound on aI and consequently A can be obtained as A = aI ≤ I − σ (P 2 ΨV ΨT Vθ̂−1 ) I, where
max

max

y

θ̂

the normalization of ther
eigenvalue is performed to ensure that the associated fraction remains
2
σmin (Pmax
ΨVy ΨT V −1 )
less than 1. Let a = 1 − σ (P 2 ΨV ΨT Vθ̂−1 ) , therefore, the shrinkage matrix becomes A = aI.
max

max

y

θ̂

The smoothing matrix corresponding to the normal distribution variance is now obtained from
the shrinkage factor as (1 − a2 )I. This guarantees that the distribution (29) has a finite variance
as t → ∞ for both constant and time-varying parameter estimation cases. This completes the
proof of the proposition.
B: Proof of Theorem 1: The existence of the projection algorithm in the parameter estimation
(j)
scheme ensures that θ̃t|t remains inside DN . According to equation (20), the a posteriori estimate
of the parameter at time t is obtained from the resampled particles of the parameter estimate
P
(j)
(j)
(j)
(j)
θ̂t|t , as θ̂t|t = N1 N
j=1 θ̂t|t , where θ̂t|t is selected from the N particles of θ̃t|t for which ρνt (yt −
(j)
h(x̂t|t , θ̃t|t )) yields higher probabilities. In order to avoid the discontinuity that is caused by
(j)
resampling, in this procedure only the particles that are maintained from time t − 1, i.e. θ̂t−1|t−1
(j)
and are propagated to time t as θ̂t|t , are considered. However, the rest of the particles that are
to be discarded in the resampling process will be replaced by the kept particles. Therefore, the
results can be generalized to all the particles.
(j)
(j)
Consider equation (17b) for generating θ̃t|t and let us substitute mt from the PE-based update
(j)
rule of (17a)-(17b) to obtain the following expression for the resampled particles θ̂t|t , namely
(j)

N
1 X (i)
θ̂
N i=1 t−1|t−1
p
(j)
(j) (j)
(j)
+ Aγt Rt ψt (t, θ̂t−1|t−1 ) + I − A2 ζt ,
(j)

θ̂t|t =Aθ̂t−1|t−1 + (I − A)

(32)

√
(j)
where I − A2 ζt denotes the evolution noise by taking into account the kernel smoothing concept. By applying the sum operator to both sides of equation (32) to construct θ̂t|t
P
PN (j)
PN 1 PN (j)
PN 1 PN (j)
(j)
1
1
1
yields, N1 N
θ̂
=
A
θ̂
+
θ̂
−A
θ̂t−1|t−1 +
j=1 t|t
j=1 t−1|t−1 N
i=1 N
i=1 t−1|t−1
N
√
PN
P
PN (j)N i=1 N1 Pi=1
(j) (j)
(j)
(j)
(j) (j)
(j)
N
N
1
1
1
2
A N j=1 γt Rt ψt (t, θ̂t−1|t−1 )+ I − A N j=1 ζt = N j=1 θ̂t−1|t−1 +A N j=1 γt Rt ψt (t, θ̂t−1|t−1 ),
P
(j) (j)
(j)
which results in θ̂t|t = θ̂t−1|t−1 + A N1 N
j=1 γt Rt ψt (t, θ̂t−1|t−1 ). Assumptions A1 and A2
ensure that the regularity conditions are satisfied according to [30]. Consequently, a differential
equation associated with (17a)-(17b) can be obtained by considering that ∆τ is a sufficiently
P
small number and t, ť are specified such that ťk=t γk = ∆τ . Through a change of time-scales
(j)
as t → τ and ť → τ + ∆τ , for a sufficiently small ∆τ , and by assuming that θ̂t−1|t−1 = θ̌, Rt =
(j)
Ř(j) , A = aI is a constant matrix, the difference equation for θ̂t|t is now expressed as
(j)

θ̂ť ≈ θ̌(j) + a∆τ Ř(j) g(θ̌(j) ),
Pt=ť

(j)

(j)

(33)
(j)

1
where g(θ̌(j) ) = ∆τ
k=t ψt (k, θ̂t−1|t−1 ). In the above derivation it is assumed that the θ̂t|t
(j)
(j)
(j)
particle is kept after resampling (that is θ̂t−1|t−1 → θ̃t|t → θ̂t|t ). Consequently, considering
Assumption A3, the differential equation associated with the evolution of each single particle is
obtained as,

March 13, 2018

DRAFT

35
(j)

dθD
d ¯ (j) T
(j)
(j)
(j)
= aRD (τ )g(θ̂D (τ )) = −aRD (τ )[ (j) J(
θ̂D )] ,
dτ
dθ̂D

(34)

where the subscript D is used to differentiate the solution of the differential equation (34) from
the solution of the difference equation (33). Now, the required convergence analysis is reduced
to investigating the properties of the deterministic continuous-time system (34).P
(j)
N
1
¯ (j)
¯ θ̂(j)
Consider the positive definite function L(θ̂t−1|t−1 ) = E(J(
j=1 J(θ̂t−1|t−1 )
t−1|t−1 )) = N
(j)
that represents the expectation of a positive definite function through N data points for θ̂t−1|t−1 .
Our goal is to evaluate the derivative of this function along the trajectories of the system (34).
(j)
According to Assumption A2, the second derivative of Q((t, θ̂t−1|t−1 )) is bounded, therefore
(j)
the summation and derivative operations commute. According to Assumption A3 for θ̂t−1|t−1 ∈
(j)
d
d
´¯ θ̂(j) (τ )) =
¯ θ̂(j)
DN , J(
Q((t, θ̂
)), exists and is
J(
)| (j)
(j)
(j) = Ē(
D

t−1|t−1

dθ̂t−1|t−1

θ̂t−1|t−1 =θ̌

dθ̂t−1|t−1

t−1|t−1

(j)
¯ θ̂(j) (τ )), and given that a > 0,
approximated by −g(θ̌(j) ). Therefore, let us define V (θ̂D ) = E(J(
D
(j)
and RD (τ ) is a positive scalar for j = 1, ..., N (which represents the trace of a positive definite
matrix at time τ ), one gets
N
d
d ¯ (j)
1 X ´¯ (j) d (j)
(j)
V (θ̂D ) = E( J(
θ̂D (τ )) =
J(θ̂D ) θ̂D (τ )
dτ
dτ
N j=1
dτ

(35)

N

−a X
(j)
(j)
(j)
[g(θ̂D (τ ))]RD (τ )[g(θ̂D )(τ )]T ≤ 0,
=
N j=1
(j)

where the equality is obtained only for θ̂D (τ ) ∈ DC . Therefore, as t → ∞ either θ̂t−1|t−1 and
consequently, πθNt−1|t−1 w.p.1 tends to DC or to the boundary of DN , where w.p.1 is with respect
to the random variables related to the parameter estimate particles. It should be noted that for
particles that have been replaced in the resampling this equality is valid since they are replaced
by particles that have satisfied (35). This completes the proof of the theorem.

C: Complexity Analysis Results: In the following Tables VIII-X the equivalent flop (EF)
complexity associated with the three methods namely, (a) our proposed dual estimation scheme,
(b) the conventional Bayesian method for state and parameter estimation [29] (where Regularized
particle filter structure is utilized to implement the filter), and (c) the recursive maximum likelihood (RML) method according to simultaneous perturbation stochastic approximation (SPSA)
algorithm [26] are presented, respectively.

March 13, 2018

DRAFT

36

TABLE VIII
T HE E QUIVALENT C OMPLEXITY FOR THE STATE ESTIMATION STEP IN THE D UAL S TRUCTURE .

(i)

x̂t|t−1

Instruction
[U1 , T1 ] = schur(Σx̂t−1|t )
R1 = randn(nx , N )
√
(i)
ωt = (U1 T1 )R1
(i)
(i)
(i)
T
= ft (x̂t−1|t−1 , θ̂t−1|t−1 λ(x̂t|t−1 ), ωt )

(i)

Σx̂t|t−1

(i)

Mult.
−
−
n3x + N n2x
−

(i)

T
ŷt|t−1 = ht (x̂t|t−1 , θ̂t−1|t−1
λ(xt|t−1 ))
PN
(i)
(i)
1
= N −1 i=1 (x̄t|t−1 − x̂t|t−1 )(x̄t|t−1 − x̂t|t−1 )T
(i)

Func. Eval.
−
−
−
N nx

Other
10n3x
N nx c1
n2x
−
−

−

−

N ny

N n2x

2N nx

−

−

−

−

−

N nx c2 + N nx c3

nx
n3x + 2N n2x + nx

N nx
n3x + (N − 1)n2x + 2N nx

−
N (nx + ny )

−
10n3x + n2x
N nx (c1 + c2 + c3 )

(i)

Reguralization and resampling to find weights wxt and x̂t|t
PN (i)
x̂t|t = N1
i=1 x̂t|t
Total

Add
−
−
(nx − 1)n2x + N (nx − 1)nx
−

TABLE IX
T HE E QUIVALENT C OMPLEXITY FOR THE PARAMETER ESTIMATION STEP IN THE D UAL S TRUCTURE .
Instruction
(j)T

(j)

ȳt|t−1 = ht (x̂t|t , θt−1|t−1 λ(x̂t|t ))
Σθ = (I − A2 )Σθ̂t−1|t−1
(j)
(j)
t = yt − ȳt|t−1
(j)
dh
ψt = dθ |x̂ ,θ(j)

q t|t t−1|t−1
(j) (j)T
= γt ( trace(t t )
[U2 , L2 ] = schur(Σθ )
R2 = randn(nθ , N ))
√
(j)
ζt = (U2 L2 )R2
(j)
(j)
(j) (j) (j)
mt = θ̂t−1|t−1 + Pt ψt t
PN
(j)
(j)
(j)
= Amt + (I − A) N1
j=1 θ̂t−1|t−1 + ζt
(j)
Pt

(j)

θ̃t|t

(j)T

ȳt|t = ht (x̂t|t , θ̃t|t λ(x̂t|t ))
(j)

(j)

Resampling to find weights, wθt , and θ̂t|t
PN
(j)
θ̂t|t = N1
j=1 θ̂t|t
P
(j)
(j)
T
θ̂
−
θ̂
Σθ̂t|t = N 1−1 N
(
t|t )(θ̂t|t − θ̂t|t )
i=1 t|t
Total

March 13, 2018

Mult.

Add

Func. Eval.

Other

−
n3θ

−
(nθ − 1)n2θ + n2θ

N ny
−

−
−

−

N ny

−

−

−

−

ny nθ

−

N + N ny
−
−
n3θ + N n2θ
N (ny nθ + nθ )

N (ny − 1) + N ny
−
−
(nθ − 1)n2θ + N (nθ − 1)nθ
N (ny − 1)nθ + N nθ

−
−
−
n2θ
−

−
10n3θ
N nθ c1
−
−

N n2θ + nθ

N n2θ + 2N nθ + N nθ + n2θ

−

−

−

−

N ny

−

−

−

−

N nθ c2
−

nθ

N nθ

−

N n2θ

2N nθ

−

−

2n3θ + 3N n2θ
+(N + 2)nθ + N nθ ny
+N (ny + 1)

2n3θ + 2N n2θ
+5N nθ + 3N ny − N
+N nθ ny

n2θ + 2N ny
+ny nθ

10n3θ + N nθ c1 + N nθ c2

DRAFT

37

TABLE X
T HE E QUIVALENT C OMPLEXITY FOR THE BAYESIAN AUGMENTED STATE AND PARAMETER ESTIMATION SCHEME [29].
Instruction
[U1 , T1 ] = schur(Σx,θ )
R1 = randn(nx + nθ , N )
√
(i)
ωt = (U1 T1 )R1
(i)

Mult.
−
−
(nx + nθ )3 + N (nx + nθ )2

Func. Eval.
−
−
−

Other
10(nx + nθ )3
N (nx + nθ )c1
(nx + nθ )2

−
n3θ

Add
−
−
(nx + nθ − 1)(nx + nθ )2
+N (nx + nθ − 1)(nx + nθ )
−
(nθ − 1)n2θ + n2θ

−
−

−
−

−

−

N (nx + nθ )

−

−

−

N ny

−

N (nx + nθ )2

2N (nx + nθ )

−

−

−

−

−

N (nx + nθ )(c3 + c2 )

nx + nθ
2n3θ + n3x + n2x (2N + 3nθ )
+n2θ (3nx + 2N ) + 4N nθ nx
+nx + nθ

N (nx + nθ )
n3x + 2n3θ + n2x (3nθ − 1 + N )
+n2θ (N − 1 + 3nx )
+nx nθ (2N − 2)

−
N (nx + nθ )
+N ny

−
10n3x + 10n3θ
n2x (30nθ + 1)
+n2θ (30nx + 1)
+nx (N c1 + 2nθ + N c3 )
+nθ (N c1 + N c3 )

(i)

ωxt = ωt (1 : nx )
(i)
ωθt = (I − A2 )ωt (nx + 1 : nx + nθ )
(i)
(i)
state/parameter augmentation: [x̂t|t−1 ; θ̂t|t−1 ] =
(i)T

(i)

(i)

(i)

(i)

[ft (x̂t|t−1 , θ̂t|t−1 λ(x̂t|t−1 ), ωxt ); θ̂t|t−1 ]
Σx,θ =

(i)
(i)
(i)T
(i)
ŷt|t−1 = ht (x̂t|t−1 , θ̂t|t−1 λ(x̂t|t−1 ))
PN
(i)
(i)
1
i=1 ([x̄t|t−1 ; θ̄t|t−1 ] − [x̂t|t−1 ; θ̂t|t−1 ])
N −1
(i)
(i)
×([x̄t|t−1 ; θ̄t|t−1 ] − [x̂t|t−1 ; θ̂t|t−1 ])T

Reguralization to find weights and resampling
(i)
(i)
to find, x̂t|t and θ̂t|t
PN (i) 1 PN (i)
1
[x̂t|t ; θ̂t|t ] = [ N i=1 x̂t|t ; N i=1 θ̂t|t ]
Total

TABLE XI
T HE E QUIVALENT C OMPLEXITY FOR THE R ECURSIVE M AXIMUM L IKELIHOOD (RML) PARAMETER SCHEME BASED ON
PARTICLE F ILTERS U SING SPSA [26].

(i)

x̃t,+

Instruction
[U1 , T1 ] = schur(Σx )
R1 = randn(nx , N )
√
(i)
ωt = (U1 T1 )R1
Generate random perturbation vector ∆t
for i = 1, ..., N sample:
(i)
(i)
(i)
= ft (x̂t−1|t−1 , (θ̂t−1|t−1 + ct ∆t )T λ(x̂t−1|t−1 ), ωt )

(i)

(i)

(i)

(i)

x̃t,− = ft (x̂t−1|t−1 , (θ̂t−1|t−1 − ct ∆t )T λ(x̂t−1|t−1 ), ωt )
Evaluate:
(i)
(i)
(i)
(i)
aθ (yt , x̃t,+ , x̂t−1|t−1 ), and aθ (yt , x̃t,− , x̂t−1|t−1 )
Evaluate:
ˆ
ˆ
dt,µ (θ̂t−1|t−1 ) = Jt (θ̂t−1|t−1 +ct ∆t )−Jt (θ̂t−1|t−1 −ct ∆t )
∇J
2ct ∆t,µ

where:
dt,µ (θ̂t−1|t−1 ± ct ∆t ) = log{ 1 PN aθ (yt , x̃(i) , x̂(i)
∇J
t,±
i=1
t−1|t−1 )}
N
Parameter Update:
dt (θ̂t−1|t−1 ),
θ̂t|t = θ̂t−1|t−1 + γt ∇J
dt (θ̂t−1|t−1 ) = [∇J
dt,1 (θbt−1|t−1 ), ..., ∇J
dt,n (θ̂t−1|t−1 )]
∇J
θ
for i = 1, ..., N sample:
(i)
(i)
(i)
(i)
T
x̃t = ft (x̂t−1|t−1 , θ̂t|t
λ(x̂t−1|t−1 ), ωt )
Regularization to evaluate weights and resampling
Total

March 13, 2018

Mult.
−
−
n3x + N n2x
−
−
−

Add
−
−
(nx − 1)n2x + N (nx − 1)nx
−
−
−

Func. Eval.
−
−
−
−
−
−

Other
10n3x
N nx c1
n2x
nθ c1
−
N nx

−
−
−
−

−
−
−
−

−
−
−
−

N nx
N (nx + nθ )c1
2N nθ
−

nθ + 1
−
2nθ
−
n2θ
−
−
−
−
n3x + N n2x + n2 θ
+3nθ + 1

2nθ − 1
−
2N nθ
−
nθ + (nθ − 1)nθ
−
−
−
−
n3x + (N − 1)n2x − N nx
+n2θ + (2N + 2)nθ − 1

−
−
2nθ
−
−
−
−
N nx
−
N nx + 2nθ

−
−
−
−
−
−
−
−
N nx (c2 + c3 )
10n3x + n2x
+nx (2N c1 + 2N + N c2 )
+nθ (c1 + N c1 + 2N )

DRAFT

