Applying GP-EKF and GP-SCKF for non-linear state estimation and fault detection in a continuous stirred-tank reactor system

Abstract

This paper proposes two novel methods for fault detection in non-linear processes. These methods apply a Gaussian process (GP) to model the underlying process, and then the extended Kalman filter (EKF) and square root cubature Kalman filter (SCKF) are used to detect faults. Accordingly, two approaches called the Gaussian process–extended Kalman filter (GP-EKF) and Gaussian process–square root cubature Kalman filter (GP-SCKF) are proposed. The most important characteristic of these proposed methods is that there is no need for an accurate model of the system. Therefore, these methods are considered non-parametric approaches of fault detection in non-linear systems. To illustrate the performance of these algorithms in fault detection, they have been used in a continuous stirred-tank reactor system (CSTR). Both proposed methods are able to detect sensor faults at an early stage.

Keywords

Extended Kalman filter fault detection Gaussian process regression non-linear systems square root cubature Kalman filter

Introduction

Progress in different sciences and industries make their applied systems bigger, more complicated and more intelligent. In these systems, the occurrence of faults is inevitable. As designing healthy systems with optimal performances is one of the most basic concerns of industrial system designers, a fault detection and diagnosis unit is an inseparable part of industrial and dynamic systems. Among fault detection and diagnosis systems, the best is the one that can detect faults within the shortest possible time and with the least loss.

Until now, several fault detection methods have been proposed for non-linear systems (Blanke and Schröder, 2003; Chen and Patton, 2012; Zhang D et al., 2013). One of the most successful approaches is the Kalman filter (KF), which has been used in dynamic systems (Doraiswami and Cheded, 2012; Safarinejadian et al., 2013; Liu J et al., 2014; Wang and Qi, 2015). One of the main disadvantages of the KF is that its application is restricted to linear systems, whereas physical systems are non-linear, inherently; therefore, the standard Kalman filter is unable to detect the faults in these systems. Consequently, the non-linear version of KF that is known as the extended Kalman filter (EKF) is used to detect faults of non-linear systems (Efimov et al., 2011; Foo et al., 2013; Zhang X et al., 2013). In this method, due to the linearization error that exists in the filter implementation, the state estimation does not show sufficient accuracy. Therefore, the present method is unable to perform well, especially in those systems that are strictly non-linear, or in those systems that need higher terms of the Taylor series expansion for the linearization of process and measurement equations. To overcome these problems, another method known as the unscented Kalman filter (UKF) has been proposed for fault detection of the non-linear systems (Cao and Ayalew, 2015; Liu H et al., 2014; Lu et al., 2014; Masoumnezhad et al., 2015). In this method, a set of sample points known as sigma points are used for the mean and the covariance estimation. Using these sigma points in the state estimation can reduce the computational complexity. However, there is more than one free parameter in the implementation of UKF and setting these parameters to reach the desired estimation is rather difficult to apply. Moreover, the performance of the filter may be interrupted in higher dimensions because, which is one of the regulative free parameters of UKF, may be negative (Julier et al., 2000). Due to the mentioned UKF drawbacks in detecting the fault of non-linear systems, the cubature Kalman filter (CKF) is proposed (Zarei and Shokri, 2014). The CKF that has the same computational load as UKF is able to generate more accurate first and second moments than UKF (Arasaratnam, 2009; Mu and Cai, 2013).

In the above-mentioned methods, an accurate model of the system is needed to detect faults. However, an accurate mathematical parametric model is never accessible, as considering all aspects of the system model is difficult; besides, some factors of modelling involved such as friction are often neglected to make the task easier. Therefore, most non-linear systems either have no parametric model or have an approximated model. In both cases, parametric methods of fault detection are neither useful nor have a desired performance. To address these advantages, researchers have proposed a non-parametric method known as the Gaussian process (GP) as a substitution for the model of the system (Rasmussen and Williams, 2006; Shi and Choi, 2011). GPs have been used in the fault detection of dynamic systems (Juricic and Kocijan, 2006). Recently, a new algorithm for the non-linear state estimation has been proposed by combining the GP and parametric methods (Deisenroth et al., 2012; Ko and Fox, 2011).

In this paper, two new methods of fault detection, called the Gaussian process–extended Kalman filter (GP-EKF) and Gaussian process–square root cubature Kalman filter (GP-SCKF), are proposed for fault detection of non-linear systems. These methods have been developed by combining the EKF and the square root of the CKF with GPs. The main advantage of the proposed approaches is that there is no need for an accurate model of the system. Furthermore, these methods have a sufficient accuracy in generating the residual signal. Moreover, by applying the square root of CKF, this method has less computational complexity and is also more stable numerically (Hao et al., 2012). Consequently, applying a highly accurate method, which has some lower computational complexity, can generate a residue that is sensitive to the faults of the system and can detect them in the minimum possible time.

The present study is structured as follows: the statement of the problem and the required necessities of presenting the proposed algorithms are described in the next section. Then, the proposed algorithms and the methods of their application in fault detection are presented. Simulation results are given and, finally, the conclusion is provided.

Problem statement

The Gaussian regression model

The main characteristic of GPs is their flexibility in modelling, as they can model the unknown functions when there is an uncertainty in the function. By using training data, a GP is able to learn the noise imposed on the function and to smooth parameters (Rasmussen and Williams, 2006). More information can be found in Safarinejadian and Kowsari (2014), Shi and Choi (2011), and references therein.

The GP is briefly introduced in the following. The non-linear regression model is given as:

y_{i} = f (x_{i}) + ε_{i}

(1)

where $ε_{i}$ is a zero-mean Gaussian noise, $ε_{i} ~ N (0, σ_{n}^{2})$ . The aim of GP regression (GPR) is to determine the Gaussian posterior distribution over function f by using the training data set $D_{f}$ . $D_{f} = 〈 X, y 〉$ , in which $X = {x_{1}, x_{2}, \dots, x_{n}}$ , includes d-dimensional $x_{i}$ s and $y = {y_{1}, y_{2}, \dots, y_{n}}$ includes $y_{i}$ observations. Function f maps $x_{i}$ s onto $y_{i}$ s by applying Equation (1). As a Gaussian random variable is completely determined with its mean and covariance, it is sufficient to select an appropriate mean and a covariance function to define an unknown function by GPR. In the GP, the prior distribution over the function is like $f = N (m_{f}, k_{f})$ , in which $m_{f} (.)$ and $k_{f} (., .)$ show the mean and covariance functions, respectively. The covariance function, which is called the kernel, is located between two x and $x'$ points and is defined as:

\begin{matrix} k_{f} (x, x') & : = E_{f} [(f (x) - m_{f} (x)) (f (x') - m_{f} (x'))] \\ = {cov}_{f} [f (x), f (x')] \end{matrix}

(2)

In a GP, a form that is selected as the covariance function is the square exponential (SE) covariance function, which is determined by using Equation (3) (Rasmussen and Williams, 2006).

k_{SE} (x_{p}, {x_{q}}^{'}) = σ_{f}^{2} e^{- \frac{1}{2} (x_{p} - {x_{q}}^{'}) W {(x_{p} - {x_{q}}^{'})}^{T}}

(3)

In this equation, $x_{p}, x_{q} \in R^{d}$ and $σ_{f}^{2}$ is the signal variance, which basically controls the predictive uncertainty in those areas where the density of the training data is low and, $W = [\frac{1}{l_{1}^{2}}, \frac{1}{l_{2}^{2}}, \dots \frac{1}{l_{d}^{2}}]$ is a diagonal matrix that includes length scale parameters $l_{i}$ , which shows how far do you need to move (along a particular axis) in input space for the function values to become uncorrelated. Length scale parameters reflect the relative smoothness of the process along the different input dimensions. The noise covariance function, which is added to function f, is determined by $k_{noise} : = δ_{pq} σ_{ε}^{2}$ . In this equation, $δ$ indicates the Kronecker delta function, and the last GP parameter is $σ_{n}^{2}$ , which is used to control the noise of the process. Finally, the prior covariance function of GP is written as $k_{f} = k_{SE} + k_{noise}$ . The reasons for selecting the mean and covariance functions have been mentioned in Deisenroth et al. (2012). By determining the prior distribution of function f, the posterior distribution over the function can also be determined. As is mentioned in Deisenroth et al. (2012), the values of the function f have a Gaussian joint distribution. Conclusively, it is possible to show observations with a Gaussian distribution as $y ~ N (0, K (X, X) + σ_{n}^{2} I)$ , in which $K (X, X)$ is the covariance matrix, which itself is calculated by the kernel function of Equation (3). In other words, $K [i, j] = [k_{SE} (x_{i}, x_{j})]$ . At first, before calculating the posterior distribution over function f, Lemma 1 and its proof is explained as a preface to computing the mentioned distribution.

Lemma 1. Consider two jointly Gaussian random variables of x and z with probability density function:

p (x, z) = N ((\begin{matrix} \bar{x} \\ \bar{z} \end{matrix}), (\begin{matrix} P_{xx} & P_{xz} \\ P_{zx} & P_{zz} \end{matrix}))

(4)

The conditional density $p (x | z)$ has a Gaussian distribution as:

p (x | z) = \frac{p (x, z)}{p (z)} = N (\hat{x}, P_{xx | z})

(5)

and its corresponding conditional mean and covariance are as follows:

\begin{matrix} \hat{x} = E [x | z] = \bar{x} + P_{xz} P_{zz}^{- 1} (z - \bar{z}) \\ P_{xx | z} = cov [x | z] = P_{xx} - P_{xz} P_{zz}^{- 1} P_{zx} \end{matrix}

(6)

Proof. As $(x, z)$ are jointly Gaussian, the marginal density $p (z)$ is also Gaussian. Based on the Bayes’ rule, the conditional distribution of $p (x | z)$ is as follows:

p (x | z) = \frac{p (x, z)}{p (z)} = \sqrt{| 2 π P_{z z} |} \exp (- \frac{1}{2} {(y - \bar{y})}^{T} P_{y y}^{- 1} (y - \bar{y})) \times 1 / \sqrt{| 2 π P_{y y} |} \exp (- \frac{1}{2} {(z - \bar{z})}^{T} P_{z z}^{- 1} (z - \bar{z}))

(7)

where $y = {[x^{T} z^{T}]}^{T}$ . To substitute non-zero mean random variables with zero mean ones, the following changes of variables are used:

\begin{matrix} ξ = x - \bar{x} \\ η = z - \bar{z} \end{matrix}

(8)

The exponent on the right-hand side of Equation (7) has the below square shape:

\begin{matrix} p & = {(\begin{matrix} ξ \\ η \end{matrix})}^{T} {(\begin{matrix} P_{xx} & P_{xz} \\ P_{zx} & P_{zz} \end{matrix})}^{- 1} (\begin{matrix} ξ \\ η \end{matrix}) - η^{T} P_{zz}^{- 1} η \\ = {(\begin{matrix} ξ \\ η \end{matrix})}^{T} (\begin{matrix} T_{xx} & T_{xz} \\ T_{zx} & T_{zz} \end{matrix}) (\begin{matrix} ξ \\ η \end{matrix}) - η^{T} P_{zz}^{- 1} η \end{matrix}

(9)

where

\begin{matrix} T_{xx}^{- 1} = P_{xx} - P_{xz} P_{zz}^{- 1} P_{zx} \\ P_{zz}^{- 1} = T_{zz} - T_{zx} T_{xx}^{- 1} T_{zx} \end{matrix}

(10)

T_{xx}^{- 1} T_{xz} = - P_{xz} P_{zz}^{- 1}

(11)

Expanding Equation (9) results in:

\begin{matrix} p & = ξ^{T} T_{xx} ξ + ξ^{T} T_{xz} η + η^{T} T_{zx} ξ + η^{T} T_{zz} η - η^{T} P_{zz}^{- 1} η \\ = {(ξ + T_{xx}^{- 1} T_{xz} η)}^{T} T_{xx} (ξ + T_{xx}^{- 1} T_{xz} η) \\ + η^{T} (T_{zz} - T_{zx} T_{xx}^{- 1} T_{xx}) η - η^{T} P_{zz}^{- 1} η \\ = {(ξ + T_{xx}^{- 1} T_{xz} η)}^{T} T_{xx} (ξ + T_{xx}^{- 1} T_{xz} η) \end{matrix}

(12)

By substituting Equations (8) and (11) with $(ξ + T_{xx}^{- 1} T_{xz} η)$ , Equation (13) is calculated.

ξ + T_{xx}^{- 1} T_{xz} η = x - \bar{x} - P_{xz} P_{zz}^{- 1} (z - \bar{z})

(13)

Therefore, the conditional density of $p (x | z)$ has the mean and the covariance given as:

\begin{matrix} \hat{x} = E [x | z] = \bar{x} + P_{xz} P_{zz}^{- 1} (z - \bar{z}) \\ P_{xx | z} = cov [x | z] = T_{xx}^{- 1} = P_{xx} - P_{xz} P_{zz}^{- 1} P_{zx} \end{matrix}

(14)

In the following, the mean and the covariance functions of the posterior distribution over function f are computed based on Lemma 1. The joint distribution of values of observations y and the predicted values of function $f (x_{*})$ are given by Equation (15) for the desired input of $x_{*}$ .

[\begin{matrix} y \\ f (x_{*}) \end{matrix}] ~ N (0, [\begin{matrix} K (X, X) + σ_{n}^{2} I & k (X, x_{*}) \\ k (x_{*}, X) & k (x_{*}, x_{*}) \end{matrix}])

(15)

G P_{μ} (x_{*}, D_{f}) = k_{*}^{T} {(K + σ_{n}^{2} I)}^{- 1} y

(16)

G P_{Σ} (x_{*}, D_{f}) = k (x_{*}, x_{*}) - k_{*}^{T} {(K + σ_{n}^{2} I)}^{- 1} k_{*}

(17)

$G P_{μ} (x_{*}, D_{f})$ and $G P_{Σ} (x_{*}, D_{f})$ , which are computed by using the training data set $D_{f}$ , show the mean and the covariance functions of GP model in a desired input of $x_{*}$ , respectively. In these equations, $K = K (X, X)$ is a kernel matrix of $n \times n$ , which is located between the training inputs and $k_{*} = k (X, x_{*})$ . Considering the matrix of GP covariance function, it is observed that the uncertainty of the model, which has the covariance $G P_{Σ}$ , is dependent on both the process noise and the correlation, which is between the input of the test and the input of the training data. Ultimately, the posterior distribution over the function is given as:

p (f (x_{*}) | x_{*}, D_{f}) ~ N (G P_{μ} (x_{*}, D_{f}), G P_{Σ} (x_{*}, D_{f}))

(18)

Learning hyper parameters

In order to calculate a suitable GP, which can approximate the unknown function f appropriately and optimally, GP parameters or those parameters that are in the kernel function must also be determined optimally. To fulfil this aim, the $θ = [W, σ_{f}, σ_{n}]$ vector is defined as the hyper parameter of GP. This parameter is learnt by maximizing the marginal likelihood of training outputs.

θ_{max} = \arg max {\log (p (y | X, θ))}

(19)

The logarithmic term of (19) can be written as:

\log (p (y | X, θ)) = - 0.5 y^{T} (K (X, X) + σ_{n}^{2} I) y - 0.5 \log | K (X, X) + σ_{n}^{2} I | - 0.5 n \log 2 π

(20)

The optimization problem is carried out by optimization technique such as the conjugate gradient. Partial derivatives, which are needed for optimization, are given as:

\partial \log (p (y | X, θ)) / \partial \log (p (y | X, θ)) \partial θ_{k} \partial θ_{k} = 0.5 tr [(K^{- 1} y) {(K^{- 1} y)}^{T} \partial K / \partial K \partial θ_{k} \partial θ_{k}]

(21)

\begin{matrix} \partial k (i, j) / \partial k (i, j) \partial σ_{f} \partial σ_{f} = 2 σ_{f} e^{- \frac{1}{2} (x_{i} - x_{j}) W {(x_{i} - x_{j})}^{T}} \\ \partial k (i, j) / \partial k (i, j) \partial σ_{n} \partial σ_{n} = 2 σ_{n} δ \\ \partial k (i, j) / \partial k (i, j) \partial W_{ii} \partial W_{ii} = - \frac{1}{2} {(x_{i} [i] - x_{j} [i])}^{2} σ_{f}^{2} e^{- \frac{1}{2} (x_{i} - x_{j}) W {(x_{i} - x_{j})}^{T}} \end{matrix}

(22)

Proposed algorithms for fault detection

In the next subsection, combining the EKF and SCKF with the GP are described. Then, the computational complexity of the proposed algorithms will be investigated. Finally, these algorithms will be used to detect faults of the system.

The combination of the Gaussian regression model with the extended Kalman filter and the square root cubature Kalman filter

The main aim of this subsection is to model $f, h$ , the process and the observation functions separately, and then to combine the Gaussian models with EKF and SCKF algorithms. Consequently, the state estimation is done with appropriate accuracy. In these methods, the covariance matrices of the process and the observation noises $Q, R$ are learnt by applying training data. The process model maps the state and the control input $(x_{k - 1}, u_{k})$ into the transient state $x_{k}$ , and the observation model maps the state $x_{k}$ into the observation $z_{k}$ . In order to learn these two models, two separate sets of training data are needed as:

\begin{matrix} D_{f} = 〈 (X, U), X' 〉 \\ D_{h} = 〈 X, Z 〉 \end{matrix}

(23)

Here, X is the matrix of states, $X' = [Δ x_{1}, Δ x_{2}, \dots, Δ x_{n}]$ is the transient state matrix and Z is a matrix obtained from observed inputs. In GP estimation, functions $f, h$ are shown by $G P^{h}$ and $G P^{f}$ . In the SCKF, the aim of predicting step is to determine $p (x_{k} | z_{1 : k - 1})$ or the posterior distribution of state $x_{k}$

p (x_{k} | z_{1 : k - 1}) = \int p (x_{k} | x_{k - 1}, u_{k - 1}) p (x_{k - 1} | z_{1 : k - 1}) d x_{k - 1}

(24)

In this equation, $p (x_{k - 1} | z_{1 : k - 1})$ is calculated from the prior step of performing the filter. It is assumed that the GPR-based function $f_{k - 1}$ has the Gaussian distribution in a desired test input as in Equation (18); therefore, the integral has a Gaussian shape, ultimately.

p (x_{k} | x_{k - 1}, u_{k - 1}) \approx N (G P_{μ} ([x_{k - 1}, u_{k - 1}], D_{f}), G P_{Σ} ([x_{k - 1}, u_{k - 1}], D_{f}))

(25)

By substituting this distribution in the EKF and SCKF:

\begin{matrix} x_{k} = {GP}_{μ}^{f} ([x_{k - 1}, u_{k - 1}], D_{f}) + ε_{k} \\ ε_{k} ~ N (0, {GP}_{Σ}^{f} ([x_{k - 1}, u_{k - 1}], D_{f})) \end{matrix}

(26)

In the measurement update step of Bayesian filters, the aim is to determine $p (x_{k} | z_{1 : k})$ distribution based on observations of the previous times and the present time that is indicated in Equation (27).

p (x_{k} | z_{1 : k}) = \frac{p (z_{k} | x_{k}) p (x_{k} | z_{1 : k - 1})}{p (z_{k} | z_{1 : k - 1})}

(27)

The distribution of function $h_{k}$ is as:

p (z_{k} | x_{k}) \approx N ({GP}_{μ}^{h} (x_{k}, D_{h}), {GP}_{Σ}^{h} (x_{k}, D_{h}))

(28)

The prior density of $p (x_{k} | z_{1 : k - 1})$ is calculated in the time prediction step. Often, the updating equation does not have the closed-form solution because of having $p (z_{k} | z_{1 : k - 1})$ , which is formulated as:

p (z_{k} | z_{1 : k - 1}) = \int p (z_{k} | x_{k}) p (x_{k} | z_{1 : k - 1}) d x_{k}

(29)

However, by applying the GPs, considering the fact that the below functions are Gaussian integrals:

\begin{matrix} z_{k} = {GP}_{μ}^{h} (x_{k}, D_{h}) + δ_{k} \\ δ_{k} ~ N (0, {GP}_{Σ}^{h} (x_{k}, D_{h})) \end{matrix}

(30)

\begin{matrix} x_{k} = {GP}_{μ}^{f} ([x_{k - 1}, u_{k - 1}], D_{f}) + ε_{k} \\ z_{k} = {GP}_{μ}^{h} (x_{k}, D_{h}) + δ_{k} \\ ε_{k} = N (0, {GP}_{Σ}^{f} ([x_{k - 1}, u_{k - 1}], D_{f})) \\ δ_{k} = N (0, {GP}_{Σ}^{h} (x_{k}, D_{h})) \end{matrix}

(31)

In order to describe GP-EKF and GP-SCKF, the steps of these algorithms are shown as Tables 1 and 2.

Table 1.

Gaussian process–extended Kalman filter (GP-EKF) algorithm.

Initialization

{\hat{x}}_{0} = E [x_{0}], P_{0} = (x_{0} - {\hat{x}}_{0}) {(x_{0} - {\hat{x}}_{0})}^{T}

Prediction:

1: Calculate the state estimate

{\hat{x}}_{k}

and process noise covariance matrix

Q_{k}

[μ_{k} Q_{k}] = GPR (({\hat{x}}_{k - 1}, u_{k - 1}), D_{f})

{\hat{x}}_{k} = {\hat{x}}_{k - 1} + μ_{k}

2: Calculate Jacobean matrix from process model

\begin{matrix} [A tmp] = \frac{\partial}{\partial x_{k - 1}} GPR (({\hat{x}}_{k - 1}, u_{k - 1}), D_{f}) \\ G_{k} = I + A \end{matrix}

3: Calculate covariance matrix

P_{k} = G_{k} P_{k - 1} G_{k}^{T} + Q_{k}

Measurement Update:

4: Calculate the estimate

{\hat{z}}_{k}

and process noise covariance matrix

R_{k}

[{\hat{z}}_{k} R_{k}] = GPR ({\hat{x}}_{k}, D_{h})

5: Calculate Jacobean matrix and the gain matrix

[H_{k} tmp] = \frac{\partial}{\partial x_{k}} GPR ({\hat{x}}_{k}, D_{h})

K_{k} = P_{k} H_{k}^{T} {(H_{k} P_{k} H_{k}^{T} + R_{k})}^{- 1}

6: Calculate the estimated state and the corresponding covariance matrix

x_{k} = x_{k} + K_{k} (z_{k} - {\hat{z}}_{k}), P_{k} = (I - K_{k} H_{k}) P_{k}

Table 2.

Gaussian process–square root cubature Kalman filter (GP-SCKF) algorithm.

Initialization

{\hat{x}}_{0} = E [x_{0}], P_{0} = S_{0} S_{0}^{T}

2: Compute the cubature points and weights

according to Table 1

Prediction

3: Evaluate the cubature points

X_{l, k - 1 | k - 1} = S_{k - 1 | k - 1} η_{l} + {\hat{x}}_{k - 1 | k - 1} l = 1, 2, \dots, m

4: Evaluate the propagated cubature points

For

i = 1, \dots 2 n_{x}

[μ Q_{k}] = GPR (X_{i, k - 1 | k - 1}, D_{f})

{X^{*}}_{i, k - 1 | k - 1} = X_{i, k - 1 | k - 1} + μ

End

5: Estimated the predicted state

{\hat{x}}_{k | k - 1} = \frac{1}{m} \sum_{l = 1}^{m} {X^{*}}_{l, k | k - 1}

6: Estimate the square-root of predicted error covariance

S_{k | k - 1} = Tria ([{χ^{*}}_{k | k - 1} S_{Q}])

where

{χ^{*}}_{k | k - 1}

\begin{matrix} {χ^{*}}_{k | k - 1} = \\ \frac{1}{\sqrt{m}} [{X^{*}}_{1, k | k - 1} - {\hat{x}}_{k | k - 1} {X^{*}}_{2, k | k - 1} - {\hat{x}}_{k | k - 1} \dots {X^{*}}_{m, k | k - 1} - {\hat{x}}_{k | k - 1}] \end{matrix}

and

S_{Q}

square-root of

Q_{k}

Update

7: Evaluate the cubature points

X_{l, k | k - 1} = S_{k | k - 1} η_{l} + {\hat{x}}_{k | k - 1} l = 1, 2, \dots, m

Where

P_{k | k - 1 |} = S_{k | k - 1} {S^{T}}_{k | k - 1}

8: Evaluate propagated cubature points

For

i = 1, \dots 2 n_{x}

[{Z^{*}}_{l, k | k - 1} R_{k}] = GPR (X_{i, k | k - 1}, D_{h})

End

9: Estimate the predicted output

{\hat{z}}_{k | k - 1} = \frac{1}{m} \sum_{l = 1}^{m} {Z^{*}}_{l, k | k - 1}

10: Estimate the square root of innovation

covariance matrix

S_{k | k - 1} = Tria ([φ_{k | k - 1} S_{R}])

Where

φ_{k | k - 1}

\begin{matrix} φ_{k | k - 1} = 1 / 1 \sqrt{m} \sqrt{m} \times \\ [{Z^{*}}_{1, k | k - 1} - {\hat{z}}_{k | k - 1} {Z^{*}}_{2, k | k - 1} - {\hat{z}}_{k | k - 1} \dots {Z^{*}}_{m, k | k - 1} - {\hat{z}}_{k | k - 1}] \end{matrix}

11: Estimate the cross-covariance matrix

P_{xz, k | k - 1} = χ_{k | k - 1} {φ^{T}}_{k | k - 1}

Where

χ_{k | k - 1}

\begin{matrix} χ_{k | k - 1} = \\ \frac{1}{\sqrt{m}} [X_{1, k | k - 1} - {\hat{x}}_{k | k - 1} X_{2, k | k - 1} - {\hat{x}}_{k | k - 1} \dots X_{m, k | k - 1} - {\hat{x}}_{k | k - 1}] \end{matrix}

12: Estimate the Kalman gain

K = (P_{XZ, K | K - 1} / {S^{T}}_{zz, k | k - 1}) / S_{zz, k | k - 1}

13: Estimate the update state and the square-root factor of the corresponding error covariance

\begin{matrix} {\hat{x}}_{k | k} = {\hat{x}}_{k | k - 1} + K (z_{k} - {\hat{z}}_{k | k - 1}) \\ S_{k | k} = Tria ([χ_{k | k - 1} - K φ_{k | k - 1} K S_{R}]) \end{matrix}

The cubature points and their corresponding weights can be calculated as presented in Table 3.

Table 3.

Compute the cubature points and their corresponding weights.

1: set the number of required cubature point,

m = 2 n_{x}

2: for

i = 1, \dots, m

a: set

η_{i} = \sqrt{\frac{m}{2}} {[1]}_{i}

where

{[1]}_{i}

is ith column of matrix

{[1]}_{n_{x} \times 2 n_{x}} = (\begin{matrix} 1 & 0 & \dots & 0 & - 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 & 0 & - 1 & \dots & 0 \\ ⋮ & ⋱ \\ 0 & 0 & \dots & 1 & 0 & 0 & \dots & - 1 \end{matrix})

w_{i} = \frac{1}{m}

End for

An investigation into the computational complexity of GP-EKF and GP-SCKF

In the time update and measurement update steps, the GP-EKF and the GP-SCKF are similar to their corresponding filters or the EKF and the SCKF. Substituting GPR with the unknown functions of the process and the observation is their only difference. Therefore, in this section only the computational complexity arisen from substituting GPR in SCKF is being investigated. At each step, the computational cost is shown by C. In the proposed method, the computational cost divided into the kernel cost $C_{kern}$ and the multiplying cost $C_{mult}$ . In each step of performing the GP-EKF algorithm, one mean, one covariance and one expansion of the Taylor series must be calculated.

C_{GP - EKF} = C_{μ} + C_{Σ} + C_{TSE}

The GP-SCKF requires the computation of one mean for each cubature point; besides, in each step of performing the GP-SCKF algorithm, one covariance must also be calculated.

C_{GP - SCKF} = 2 n_{x} C_{μ} + C_{Σ}

(32)

For more accurate analysis, the computational complexity of $C_{μ}$ and $C_{Σ}$ are broadly investigated as follows. To calculate the GPs mean function, the kernel function must be analysed for n number of training data $x_{i}$ with test input $x_{*}$ . Moreover, as one GP is trained separately for each output dimension, the GPs mean function must be repeated for each output dimension.

C_{μ} = nd (C_{\ker n} + C_{mult})

(33)

Furthermore, it must be considered that $K^{- 1} y$ is independent of $x_{*}$ and it is possible to neglect it within the cost of calculations. In order to calculate the computational complexity of the covariance function, it is assumed that $k_{*}$ was calculated. Calculating $k_{*}^{T} K^{- 1} k_{*}$ needs $n (d + 1)$ for each output dimension and $k (x_{*}, x_{*})$ is merely an analysis of one kernel function.

C_{Σ} = nd (d + 1) C_{mult} + d C_{\ker n}

(34)

Finally, the linearization of GPs needs computing the Jacobin of kernel functions along with $K^{- 1} y$ multiplication. Computing each Jacobin element requires two multiplication operations. The partial derivative of the Jacobin matrix can be calculated by applying the following method.

\partial (k (x_{*}, x)) / \partial (k (x_{*}, x)) \partial (x_{*} [i]) \partial (x_{*} [i]) = - w_{ii} (x_{*} [i] - x [i]) k (x_{*}, x)

Therefore, for every output dimension, the Jacobin complexity of the kernel function is equal to $2 nd C_{mult}$ along with $K^{- 1} y$ multiplication and its cost is $nd C_{mult}$ for each dimension. Therefore, the cost of linearization is $C_{TSE} = 3 n d^{2} C_{mult}$ .

Fault detection in non-linear systems

Assuming that there is a fault in the measurement sensor, the total system can be represented as:

\begin{matrix} x_{k} = f (x_{k - 1}, u_{k - 1}) + v_{k - 1} \\ z_{k} = h (x_{k}) + w_{k} + L_{k} f_{k} \end{matrix}

(35)

where $f_{k}$ is the fault occurring at time k, and $L_{k}$ is the known matrix of the fault. In other words, an additive fault has occurred in the sensors. One of the most famous techniques in fault detection of dynamic systems is comparing the output of the healthy mode with the output of the faulty mode of the system, and using the generated residual signal from this comparison to detect the fault. Totally, the process of the fault detection can be explained in two steps. The first step is generating the residue in a system and the second step is devoted to the analysis of this generated residual.

Residual generation

The first step in the fault detection process is generating a signal, which is sensitive to the fault and is able to show the fault at the early stages. The simplest residual signal is defined as:

r_{k} = y_{k} - {\hat{y}}_{k}

(36)

where $r_{k}$ is the generated residual signal, $y_{k}$ is the output of the real system and ${\hat{y}}_{k}$ is the estimated output.

The residual analysis

Ideally speaking, the residue is zero if there is no fault in a system; otherwise, the residue is non-zero. In this case, the residual analysis is defined as:

{\begin{matrix} r_{k} = 0 fault free \\ r_{k} \neq 0 faulty \end{matrix}

(37)

However, in many systems, the residue is non-zero even at fault-free conditions. Therefore, the decisive function that is written in Equation (37) is not an appropriate residual analysis. To analyse the residual properly, it is possible to define a statistical decisive function as:

{\begin{matrix} r_{k} ≺ \bar{r} - λ σ_{r} faulty \\ \bar{r} - λ σ_{r} \leq r \leq \bar{r} + λ σ_{r} fault free \\ r_{k} ≻ \bar{r} + λ σ_{r} faulty \end{matrix}

(38)

where

\begin{matrix} \bar{r} = E [r_{k}] = \frac{1}{m} \sum_{k = 1}^{m} r_{k} \\ σ_{r}^{2} = E [{(r_{k} - \bar{r})}^{2}] = \frac{1}{m} \sum_{k = 1}^{m} {(r_{k} - \bar{r})}^{2} \end{matrix}

In selecting $λ$ , there must be a trade-off between the maximum possibility of fault occurrence and the minimum possibility of the false alarm.

Simulation results

To show the performance of the two suggested algorithms and to compare their performances, the implementation and the simulation results of the GP-EKF and the GP-SCKF will be explained as follows. The studied scenario is on a continuous stirred-tank reactor (CSTR).

As the CSTR systems transform raw materials into valuable chemicals, they are the most important part of many chemical, biochemical, polymer and petroleum processes. They can operate at low or high temperatures, in batch, fed batch or continuous mode.

The CSTR is widely used in many chemical plants. Due to changes in process parameters, the accuracy of final product can be reduced. In order to obtain an accurate final product, the faults occurring in the CSTR during the chemical reaction need to be diagnosed. If not, the faults may lead to degradation of the performance of the system

The process and measurement equation of this system is as follows:

\begin{matrix} {\overset{\cdot}{x}}_{1} = \frac{U_{f_{0}}}{x_{4}} (C_{a 0} - x_{1}) - x_{1} K_{0} e \frac{- x_{6}}{R x_{2}} + w_{1} \\ {\overset{\cdot}{x}}_{2} = \frac{U_{f_{0}}}{x_{4}} (T_{0} - x_{2}) - \frac{β x_{1} K_{0} e \frac{- x_{6}}{R x_{2}}}{ρ C_{p}} - \frac{x_{5} (x_{2} - x_{3})}{x_{4} ρ C_{p}} + w_{2} \\ {\overset{\cdot}{x}}_{3} = \frac{U_{f_{j}}}{V_{j}} (T_{cin} - x_{3}) + \frac{x_{5} (x_{2} - x_{3})}{V_{r} ρ_{j} C_{j}} + w_{3} \\ {\overset{\cdot}{V}}_{r} = \frac{ρ_{0} U_{f_{0}}}{ρ} - U_{f} + w_{4} \\ {\overset{\cdot}{x}}_{5} = 0 \\ {\overset{\cdot}{x}}_{6} = 0 \end{matrix}

(39)

where $x (t) = {[x_{1} x_{2} x_{3} x_{4} x_{5} x_{6}]}^{T} = {[C_{a} T_{r} T_{j} V_{j} U A_{j} E]}^{T}$ . The steady-state data for this process are shown in Table 4.

y (t) = [\begin{matrix} C_{a} \\ T_{r} \\ T_{j} \\ V_{r} \end{matrix}] + v

(40)

Table 4.

Steady state parameters of continuous stirred-tank reactor system (CSTR).

Process variable	Steady-state
Product concentration ( $C_{a}$ )	0.3684 km/m³
Reactor temperature ( $T_{r}$ )	333K
Coolant flow rate ( $F_{j}$ )	9.937 × 10⁻³ m³/s
Product flow rate (F)	2.72 × 10⁻³m³/s
Feed flow rate ( $F_{0}$ )	2.65 × 10⁻³ m³/s
Feed concentration ( $C_{a 0}$ )	7.128 kmol
Feed temperature ( $T_{0}$ )	296.89K
Inlet coolant temperature ( $T_{Cin}$ )	288K
Outlet coolant temperature ( $T_{j}$ )	319.76K
Reactor volume ( $V_{r}$ )	6.739 m³
Jacket volume ( $V_{j}$ )	0.4467 m³
Heat transfer term ( $U A_{j}$ )	10⁵ W/K
Reaction rate constant ( $K_{0}$ )	1.696 × 10¹³ $1 / s$
Activation energy term ( $β$ )	⁻9 × 10⁷ J/kmol
Feed density ( $ρ_{0}$ )	936.7 kg/m³
Product density ( $ρ$ )	912.9 kg/m³
Coolant density ( $ρ_{j}$ )	1008 kg/ m³
Heat capacity of product ( $C_{P}$ )	3368 J/kgK
Heat capacity of coolant ( $C_{j}$ )	4203 J/kgK

w and R are the process and measurement noises with covariance Q and R, respectively. Q and R are consider as follows in this simulations: $Q = diag ([5.43 e - 5 0.36 0.2186 0.0182 0.10.01])$ ; $R = ([0.023 0.36 0.2186 0.0182])$

In these simulations, 500 random points have been selected for the training data set. The simulation results have been tested after 50 times of performing Monte Carlo simulation on 50 random data.

State estimation

To show the performance of GP-EKF and GP-SCKF and also to check the accuracy of these algorithms, the root mean square error (RMSE) of these methods for 50 Monte Carlo runs has been calculated and results are shown in Figures 1 and 2, and Table 5. The accuracy of the GP-SCKF method is better than the GP-EKF in state estimation as seen in Table 5.

Table 5.

Root mean square error of Gaussian process–extended Kalman filter (GP-EKF) and Gaussian process–square root cubature Kalman filter (GP-SCKF) in 50 Monte Carlo simulations.

	$x (1)$	$x (2)$	$x (3)$	$x (4)$	$x (5)$	$x (6)$
GP-EKF	0.0000	0.0124	0.0823	0.0000	0.0000	0.0000
GP-SCKF	0.0000	0.0110	0.0351	0.0000	0.0000	0.0000

A comparison between the computational complexities of the two algorithms is shown in Table 6. In the table, it can be seen that the computational complexity of the GP-EKF algorithm is approximately six times more than the computational complexity of the GP-SCKF algorithm.

Table 6.

Computational complexity in 50 Monte Carlo runs.

Algorithm	Computation time (in seconds)
$C_{GP - EKF}$	1.2277
$C_{GP - SCKF}$	0.1971

Fault detection

To show the performance of these methods in non-linear fault detection, three kinds of faults in the sensor of the system has been analysed. These faults are known as abrupt, incipient and intermittent, and all are considered additive faults. In these simulations, $L_{k} = [1 1 1 1]^{T}$ is considered and $f_{k}$ is imposed on the system in the form of three different faults as:

\begin{matrix} f_{k, abrupt} = {\begin{matrix} 520 \leq k \leq 50 \\ 0 otherwise \end{matrix} \\ f_{k, incipient} = {\begin{matrix} 5 k 20 \leq k \leq 40 \\ 0 otherwise \end{matrix} \\ f_{k, intermmitent} = {\begin{matrix} 10 20 \leq k \leq 30 \\ 10 40 \leq k \leq 50 \\ 0 otherwise \end{matrix} \end{matrix}

(41)

The additive fault detection process by applying the GP-EKF and GP-SCKF algorithm

In this simulation, $λ = 3$ is chosen. In other words, when the fault occurs, one is 98.4% sure that a fault has occurred. Figure 3 indicates the residual signal of the system that is generated by the GP-EKF algorithm, assuming that an abrupt fault is imposed upon the system. Furthermore, Figures 4 and 5 show the GP-EKF generated residual signals when there are incipient and intermitted faults in the system, respectively. As observed, when each of these mentioned faults is imposed upon the system, the residual signal exceeds the limit of the fault threshold and the amplitude of the residual signal makes detecting the fault possible in the early stages. The GP-EKF algorithm has a very good performance in fault detection; however, this algorithm has a high bulk of mathematical computations, which is considered its main shortcoming especially when the system is strictly non-linear.

Figure 6 indicates the residual signal that is generated by the GP-SCKF algorithm when an abrupt fault has occurred within the system. Figure 7 shows the system residual signal when an incipient fault is imposed upon the system and Figure 8 manifests the residual signal generated by the GP-SCKF algorithm when there is an intermittent fault in the system. As observed in Figures 6, 7 and 8, the GP-SCKF algorithm can detect the fault as well as the GP-EKF algorithm. However, the main advantage of the GP-SCKF algorithm is that it has fewer mathematical computations than the GP-EKF algorithm Moreover, unlike the GP-EKF algorithm, it is possible to use the GP-SCKF algorithm in strictly non-linear systems.

The delay times of these methods for fault detection have been shown in Table 7. According to this table, occurrence of fault can be detected in the shortest time or early stage.

Table 7.

Delay time of fault detection.

	Abrupt fault	Incipient fault	Intermittent fault
GP-EKF	4.1053e-05	4.4623e-05	4.3284e-05
GP-SCKF	3.8376e-05	4.1053e-05	4.1499e-05

GP-EKF, Gaussian process–extended Kalman filter; GP-SCKF, Gaussian process–square root cubature Kalman filter.

Non-additive fault detection process by applying GP-EKF and GP-SCKF algorithms

In the proposed fault detection methods, a residual signal will be obtained from the difference between measured output and estimated output ( $r_{k} = y_{k} - {\hat{y}}_{k}$ ). As in any change in the functions, the states and the parameters of the system can result in a change in the estimated output, and the residual signal will change compared with the fault-free mode. Therefore, the proposed methods can be used not only for additive faults, but also non-additive faults can be detected. To show the performance of the two proposed methods in the case of non-additive fault detection, a fault was considered defined by Equation (42) and its simulation results are shown in Figures 9 and 10. As it can be seen in these figures, both of the algorithms were able to detect the fault correctly.

y_{fault} = {\begin{matrix} y 0 < t < 20 \\ 1.5 y 20 \leq t \leq 50 \end{matrix}

(42)

Discussion

It is illustrated that without needing an accurate model of the system, the state estimation can be done with appropriate accuracy. The simulation results show that GP-SCKF is more accurate than the GP-EKF. Furthermore, it has lower computational complexity compared with GP-EKF. Therefore, it can detect faults faster than the GP-EKF. A brief comparison between the fault detection methods is available in Table 8.

Table 8.

Advantages and disadvantages fault detection methods.

Fault detection methods	Advantages and disadvantages
Propose methods	These methods have good performance in fault detection and state estimation although they do not need healthy and faulty data as well as the accurate model of the system
Parametric methods	In order to having good performance, these methods need the accurate model of the system
Non-parametric methods	Need both healthy and faulty data of the system however collecting these data points are difficult in practice in many systems

Conclusion

In this paper, two new methods are proposed for the purpose of non-linear fault detection systems. The algorithm of the proposed methods, their computational complexity, advantages and disadvantages have all been analysed; furthermore, these methods have applicably been tested to find out their strength in the fault detection and fault diagnosis domain. The most important advantage of these algorithms compared with the previous ones is that they do not need to use accurate equations of the system and have high accuracy in generating the residual signal. Ultimately, in order to implement the proposed method in detecting the fault of non-linear systems, a comparison has been made between these algorithms. The results obtained from this comparison have proved the superiority of the GP-SCKF algorithm, especially due to its low complexity of computations and high accuracy.

Footnotes

Appendix

Conflict of interest statement

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Arasaratnam

(2009) Cubature Kalman filtering: theory & applications. Ph.D. dissertation, Department of Electrical & Computer Engineering, McMaster University.

Blanke

Schröder

(2003) Diagnosis and Fault-tolerant Control. Berlin: Springer.

Cao

Ayalew

(2015) Estimation and predictive control of nonlinear diffusion processes with application to drying of coatings. Proceedings of the Institute Mechanical Engineers, Part J: Systems and Control Engineering 229: 235–249. doi:10.1177/0959651814559759

Chen

Patton

(2012) Robust Model-based Fault Diagnosis for Dynamic Systems. Berlin: Springer.

Deisenroth

Turner

Huber

. (2012) Robust filtering and smoothing with Gaussian processes. IEEE Transactions on Automatic Control 57: 1865–1871.

Doraiswami

Cheded

(2012) Kalman filter for parametric fault detection: an internal model principle-based approach. IET Control Theory & Applications 6: 715–725.

Efimov

Zolghadri

Simon

(2011) Fault detection improving by extended Kalman filter adjustment for oscillatory failure case in aircrafts. EUCASS.

Foo

GHB

Zhang

Vilathgamuwa

(2013) A sensor fault detection and isolation method in interior permanent-magnet synchronous motor drives based on an extended Kalman filter. IEEE Transactions on Industrial Electronics 60: 3485–3495. doi:10.1109/TIE.2013.2244537

Hao

Yang

Chen

. (2012) Square root cubature Kalman filter. Journal of Projectiles, Rockets, Missiles and Guidance 32: 169–172.

10.

Julier

Uhlmann

Durrant-Whyte

(2000) A new method for the nonlinear transformation of means and covariances in filters and estimators. IEEE Transactions on Automatic Control 45: 477–482.

11.

Juricic

Kocijan

(2006) Fault detection based on Gaussian process model. In:Proceedings of the 5th Vienna Symposium on Mathematical Modeling-MathMod, Vienna, Austria.

12.

Fox

(2011) Learning GP-Bayes Filters via Gaussian process latent variable models. Autonomous Robots 30: 3–23.

13.

Liu

. (2014) Fault diagnosis of hydraulic servo system using the unscented Kalman filter. Asian Journal of Control 16(6): 1713–1725.

14.

Liu

Wong

DSH

Chen

D-S

(2014) Bayesian filtering of the smearing effect: fault isolation in chemical process monitoring. Journal of Process Control 24: 1–21.

15.

van Eykeren

van Kampen

. (2014). Adaptive hybrid unscented Kalman filter for aircraft sensor fault detection, isolation and reconstruction. In: Throckmorton

(ed.) Proceedings of the AIAA Guidance, Navigation and Control Conference, pp. 1–18. Reston, VA: AIAA.

16.

Masoumnezhad

Jamali

Nariman-zadeh

(2015) Robust GMDH-type neural network with unscented Kalman filter for non-linear systems. Transactions of the Institute of Measurement and Control 0142331215573066. doi:10.1177/0142331215573066

17.

Cai

(2013) Likelihood-based iteration square-root cubature Kalman filter with applications to state estimation of re-entry ballistic target. Transactions of the Institute of Measurement and Control 35: 949–958. doi:10.1177/0142331212459880

18.

Rasmussen

Williams

(2006) Gaussian Processes for Machine Learning. Cambridge, MA: MIT Press.

19.

Safarinejadian

Kowsari

(2014) Fault detection in non-linear systems based on GP-EKF and GP-UKF algorithms. Systems Science and Control Engineering Open Access Journal 2: 610–620.

20.

Safarinejadian

Zarei

Ramezani

(2013) Fault diagnosis of induction motors using a recursive Kalman filtering algorithm. International Review of Electrical Engineering8.

21.

Shi

Choi

(2011) Gaussian Process Regression Analysis for Functional Data. Boca Raton, FL: CRC Press.

22.

Wang

(2015) Fault diagnosis for flight control systems using subspace method and adaptive two-stage Kalman filter. Transactions of the Institute of Measurement and Control 0142331215596805. doi:10.1177/0142331215596805

23.

Zarei

Shokri

(2014) Robust fault detection based on nonlinear unknown input observer. Measurement 48: 355–367.

24.

Zhang

Wang

Q-G

. (2013) Fuzzy-model-based fault detection for a class of nonlinear systems with networked measurements. IEEE Transactions of Measurement 62: 3148–3159. doi:10.1109/TIM.2013.2272865

25.

Zhang

Foo

Don Vilathgamuwa

. (2013) Sensor fault detection, isolation and system reconfiguration based on extended Kalman filter for induction motor drives. IET Electrical Power Applications 7: 607–617. doi:10.1049/iet-epa.2012.0308