Fault Detection and Isolation

Abstract

Fault diagnosis of a class of linear multiple‐input and multiple‐output (MIMO) systems is developed here. An emulator‐based scheme is proposed to detect and isolate faults in a system formed by interconnected subsystems. Emulators, which are hardware or software devices, are connected to the input and measurement outputs in cascade with the subsystems whose faults are to be diagnosed. The role of an emulator is to induce variations in cascade combination of the nominal fault‐free subsystem so as to mimic the actual perturbations that may occur in the subsystem during the offline identification phase. The emulator‐generated data are employed in the reliable identification of the nominal system, the associated Kalman filter, and a map that relates the emulator parameters to the feature vector. In the operational stage, the Kalman filter residual is used to detect a fault in the system; the emulator parameter that has varied is estimated, and using the emulator‐feature vector map, the faulty subsystem is isolated. The main contributions of this work are accurate and reliable identification of the system, the fault diagnosis of multivariable systems using feature vector-emulator map fault diagnosis of multivariable systems, and the establishment of the key properties of the Kalman filter for fault detection. The proposed scheme was successfully evaluated on a number of simulated as well as physical systems.

1. Introduction

Fault detection and isolation (FDI) of physical systems—especially mission critical systems including nuclear reactors, aircraft, automotive systems, spacecraft, autonomous vehicles, and fast rail transportation—is becoming increasingly important in recent times thanks mainly to advances in sensors, computing, and communication technologies. It still poses a challenge in view of the stringent and conflicting requirements, high probability of correct detection and isolation, low false alarm probability, and timely decision on the fault status.

The identification of the system model is crucial to the performance of the fault diagnosis scheme. The more accurate the identified model, the higher is the probability of correct diagnosis and the lower is the false alarm probability. The reliability and accuracy of the identification hinges on ensuring that the identified model is captured completely and what is leftover is the information‐less zero‐mean white noise process. As the Kalman filter is a zero‐mean white noise process if and only if there is no mismatch between the identified model and the model of the system, the identification scheme should minimize the residual of the Kalman filter—instead the equation error, which in general, is a colored noise [1]. The widely popular, consistent, and efficient scheme that meets the above state requirement is the prediction error method (PEM) [2]. The PEM identifies the system by minimizing the residual of the Kalman filter.

A physical system is subject to perturbation resulting from the variations of the parameters and effects nonlinearities resulting in the deviation in the neighborhood of the nominal operating point. A model identified at a nominal operating point will not capture the static and the dynamic behavior of the perturbed system. To overcome this, an emulator, which is a hardware or a software device, is connected to either an accessible input or an accessible output in cascade with a subsystem to mimic its operating scenarios [3–5]. The powerful concept of emulators, which is employed to mimic the likely operating scenarios for single‐input and single‐output (SISO) system, is extended to multiple‐input and multiple‐output (MIMO) and multiple‐input and single‐output (MISO) system. The system is identified and the feature vector‐emulator map is estimated from the emulator‐generated data covering all likely operating scenarios including the normal and the faulty ones similar in spirit to that employed in training the neural network [6]. The identified nominal model, an optimal nominal model, is robust to model perturbation in the neighborhoods of the nominal operating point. It may be worth noting that the conventional scheme uses only the input‐output data from the system in the nominal operating scheme.

There are essentially three approaches to the failure detection and isolation problem: the non‐parametric approach, the parametric approach, and the combined approach. The non‐parametric approach is based on analyzing a residual. The residual is defined as a signal, which is ideally non‐zero in a statistical sense when there is a failure present, and zero otherwise. The residual may be generated using Kalman filters, observers, unknown‐input observers, other forms of detection filters, and parity equations [7–12]. In view of the following key properties of the Kalman filter listed below, the Kalman filter is deemed the most preferable for both fault detection and fault isolation [1]:

1. Model matching: The residual is a zero‐mean white noise process if and only if there is no mismatch between the actual model of the system and its identified model embodied in the Kalman filter, that is, and its variance is minimum.

2. Optimal estimation: The estimate is optimal in the sense that it is the best estimate that can be obtained by any estimator in the class of all estimators that are constrained by the same assumptions.

3. Robustness: Thanks to the feedback (closed‐loop) configuration of the Kalman filter with residual feedback, the Kalman filter provides the highest robustness against the effect of disturbance and model variations.

4. Model mismatch: If there is a model mismatch, the residual will not be a zero‐mean white noise process and an additive term termed fault‐indicative term. The fault‐indicative term is affine in the deviation in the linear regression or the transfer function model.

The feature vector‐emulator map relating the deviation of the feature vector and variations of the emulator parameter is used for fault isolation if a fault is detected. The influence vector, which is the partial derivative of the feature vector with respect to an emulator parameter, plays a crucial role in pinpointing the faulty subsystem and tracks its parameter variation.

The main contributions here are the development of emulator‐based system identification, and estimation of the feature vector‐emulator map and its application to performance monitoring and fault diagnosis of multivariable system. The key properties of the Kalman filter, including model matching, whitening of the equation error, and residual expression for the model‐mismatch case, are established for MIMO, MISO, and SISO systems.

The chapter is organized as follows. In Section 2, the mathematical model of the multiple‐input and multiple‐output system in state‐space, frequency‐domain, and a linear regression form is developed. The multiple‐input and single‐output and the single‐input, single‐output models are derived. Modeling of faults is also given. In Section 3, the concept of emulators, the generation of emulator‐perturbed data, and its role in the identification of the system, the estimation of the feature vector‐emulator map for fault isolation is developed. In Section 4, the identification of the system and the associated Kalman filter using prediction error method is suggested. The feature vector‐emulator map is estimated using the expression of the Kaman filter residual in the model‐mismatch case. In Section 5, the model of the Kalman filter, residual model, and the key properties of this filter are given. The key properties of the residual are established including whitening of the equation error, and expressions for the residual for the model‐mismatch case. In Section 6, Bayesian approach to fault diagnosis is explained. Finally, in Sections 7 and 8, the successful evaluation of the proposed scheme on both a simulated and physical system is given, respectively.

2. Mathematical model of the system

The MIMO state‐space model of the system denoted $(\mathit{A},\mathit{B},\mathit{C})$ is given by

$\begin{array}{l}\mathit{x}(k+1)=\mathit{A}\mathit{x}(k)+\mathit{B}\mathit{r}(k)+{\mathit{E}}_w\mathit{w}(k)\\ \mathit{y}(k)=\mathit{C}\mathit{x}(k)+\mathit{v}(k)\end{array}$E1

where $\mathit{x}(k)={\left[\begin{array}{ccccc}{x}_1(k)& x_2(k)& x_3(k)& \mathrm{\dots }& x_n(k)\end{array}\right]}^T$, $\mathit{y}(k)={\left[\begin{array}{ccccc}{y}_1(k)& y_2(k)& y_3(k)& \mathrm{\dots }& y_q(k)\end{array}\right]}^T$, $\mathit{r}(k)={\left[\begin{array}{ccccc}{r}_1(k)& r_2(k)& r_3(k)& \mathrm{\dots }& r_p(k)\end{array}\right]}^T$, $\mathit{w}(k)$ and $\mathit{v}(k)$_, are respectively, $nx1$ state vector, $qx1$ output, $px1$ input to the system, $px1$ disturbance and $qx1$ measurement noise; $\mathit{A}$, $\mathit{B}$, $\mathit{C}$, ${\mathit{E}}_w$ are $nxn$ state transition, $nxp$ input, and $qxn$ output and $nxp$ input disturbance matrices; $\mathit{A}$ and $\mathit{C}$ are block diagonal matrices; $\mathit{A}=\left[\begin{array}{cccc}{\mathit{A}}_1& \mathbf{0}& .& \mathbf{0}\\ \mathbf{0}& {\mathit{A}}_2& .& \mathbf{0}\\ .& .& .& .\\ \mathbf{0}& \mathbf{0}& .& {\mathit{A}}_q\end{array}\right];\mathit{B}=\left[\begin{array}{c}{\mathit{B}}_1\\ {\mathit{B}}_2\\ .\\ {\mathit{B}}_q\end{array}\right];{\mathit{E}}_w=\left[\begin{array}{c}{\mathit{E}}_{w1}\\ {\mathit{E}}_{w2}\\ .\\ {\mathit{E}}_{wq}\end{array}\right];\mathit{C}=\left[\begin{array}{cccc}{\mathit{C}}_1& \mathbf{0}& .& \mathbf{0}\\ \mathbf{0}& {\mathit{C}}_2& .& \mathbf{0}\\ .& .& .& .\\ \mathbf{0}& \mathbf{0}& .& {\mathit{C}}_q\end{array}\right]$, ${\mathit{A}}_j$, ${\mathit{B}}_j$, ${\mathit{E}}_{wj}$_,and ${\mathit{C}}_j$ are, respectively, $n_{j}x{n}_j$, $n_{j}xp$, $n_{j}xp$, and $1x{n}_j$ matrices. The output of the system is corrupted by disturbance $\mathit{w}(k)$ and measurement noise $\mathit{v}(k)$; $\mathit{G}(z)=\mathit{C}{\left(z\mathit{I}-\mathit{A}\right)}^{-1}\mathit{B}=D^{-1}(z)\mathit{N}(z)$; $\mathit{I}$ is an identity matrix; $D(z)=\left|\left(z\mathit{I}-\mathit{A}\right)\right|=1+{\displaystyle \sum _{\ell =1}^n{a}_{\ell }}{z}^{-\ell }$; ${\mathit{B}}_j=\left[\begin{array}{cccc}{B}_{j1}& B_{j2}& .& B_{jp}\end{array}\right]$; ${\mathit{E}}_w=\left[\begin{array}{cccc}{E}_{wj1}& E_w{}_{j2}& .& E_{wjp}\end{array}\right]$.

We assume that the system is controllable and observable, that is, $\left(\mathit{A},\mathit{C}\right)$ is observable, $\left(\mathit{A},\mathit{B}\right)$ is controllable, implying that all the states may be estimated from the input and the output data, and the input affects all the states. The disturbance $\mathit{w}(k)$ and the measurement noise $\mathit{v}(k)$ are assumed zero‐mean white noise processes. The covariance of $w(k)$ and $v(k)$ are

$E[\mathit{w}{\mathit{w}}^T]=\mathit{Q}\text{and}E[\mathit{v}{\mathit{v}}^T]=\mathit{R}$E2

where $\mathit{Q}$ and $\mathit{R}$ are positive definite and positive semi‐definite matrices, $\mathit{Q}>0$ and $\mathit{R}\ge 0$. The covariances $\mathit{Q}$ and $\mathit{R}$ are not known a priori.

The MIMO model in the frequency domain is

$\mathit{y}\left(z\right)=\mathit{G}\left(z\right)\mathit{r}\left(z\right)+\mathit{\vartheta }\left(z\right)$E3

where $\mathit{G}\left(z\right)$ is $qxp$ matrix transfer function, and $\mathit{N}(z)$ is the $qxp$ numerator matrix; $\mathit{\vartheta }(z)$ is the $qx1$ is the effect of disturbance $\mathit{w}(k)$ and the measurement noise $\mathit{v}(k)$ on the output $\mathit{y}\left(z\right)$.

$\mathit{\vartheta }(z)=\mathit{C}{\left(z\mathit{I}-\mathit{A}\right)}^{-1}{\mathit{E}}_w\mathit{w}(z)+\mathit{v}(z)$E4

2.1. Single‐input single‐output pairing

A single‐input single‐output (SISO) model derived from the state‐space model relating the input $r_i\left(z\right)$, and its associated output, termed $y_{ji}(z)$, which is the same as the output $y_j(z)$ when the input is $r_i\left(z\right)$ and the rest of the inputs $r_j\left(z\right)=0$ for $j\ne i$, is

$y_{ji}\left(z\right)=G_{ji}\left(z\right)r_i\left(z\right)+{\vartheta }_{ji}(z)$E5

where $G_{ji}\left(z\right)={\mathit{C}}_j{\left(z\mathit{I}-{\mathit{A}}_j\right)}^{-1}{\mathit{B}}_{ji}=D_j^{-1}\left(z\right){\mathit{N}}_{ji}\left(z\right)$; and ${\vartheta }_{ji}(z)={\mathit{C}}_j{\left(z\mathit{I}-{\mathit{A}}_j\right)}^{-1}{\mathit{E}}_{wji}{w}_i(z)$. The transfer function $G_{ji}\left(z\right)$ may in general be a cascade combination of subsystems $\left\{G_{ji\ell }\left(z\right)\right\}$:

$G_{ji}\left(z\right)={\displaystyle \prod _{\ell }{G}_{ji\ell }\left(z\right)}$E6

The subsystems $G_{ji\ell }\left(z\right)$ may, for example, be a transfer function of a controller, an actuator, a plant, or a sensor associated with a position control system, process control system, magnetic levitation system, or other systems [4].

Expressing the frequency‐domain model (5) in a linear regression form yields

$y_{ji}(k)={\mathit{\psi }}_{ji}^T\left(k\right){\mathit{\theta }}_{ji}+{\upsilon }_{ji}\left(k\right)$E7

where ${\upsilon }_{ij}\left(z\right)=D_j(z){\vartheta }_{ij}\left(z\right)$; ${\mathit{\psi }}_{ji}^T(k)$ is $1x2n_j$ regression vector formed of the regression vectors, formed ${\mathit{\psi }}_{yji}^T(k)$ associated with $y_{ji}(k)$, and ${\mathit{\psi }}_{ri}^T(k)$ associated with input $r_i\left(k\right)$:

${\mathit{\psi }}_{ij}^T(k)=\left[\begin{array}{cc}{\mathit{\psi }}_{yji}^T(k)& {\mathit{\psi }}_{ri}^T(k)\end{array}\right]$E8

${\mathit{\psi }}_{yji}^T(k)=\left[\begin{array}{cccc}-y_{ji}(k-1)& -y_{ji}(k-2)& .& -y_{ji}(k-n_j)\end{array}\right]$; ${\mathit{\psi }}_{ri}^T(k)=\left[\begin{array}{cccc}{r}_i(k-1)& r_i(k-2)& .& r_i(k-n_j)\end{array}\right]$; ${\mathit{\theta }}_{ji}$ is $2n_j\mathrm{\text{x}}1$ feature vector formed of the $n_j$ coefficients of the denominator polynomial $D_j\left(z\right)$ and the numerator polynomial ${\mathit{N}}_{ij}\left(z\right)$:

${\mathit{\theta }}_{ji}={\left[\begin{array}{cc}{\mathit{\theta }}_{yj}& {\mathit{\theta }}_{rji}\end{array}\right]}^T$E9

Remarks: In the operational stage, we may not have access to the output $y_j\left(k\right)$, termed $y_{ji}\left(k\right)$, generated by the input $r_i\left(k\right)$ alone when rest of the inputs are set to zero. It is estimated during the identification phase of the multi‐input and single‐output model relating the accessible output $y_j(k)$ generated by all the inputs $\mathit{r}\left(k\right)$.

2.2. Multi‐input and single‐output pairing

Using Eq. (5), the output $y_j\left(z\right)$ is the output due to all the inputs $\mathit{r}\left(k\right)$ of MISO system, which is

$y_j\left(z\right)={\displaystyle \sum _{i=1}^p{y}_{ji}(z)}={\mathit{G}}_j\left(z\right)\mathit{r}(z)+{\vartheta }_j\left(z\right)$E10

where ${\mathit{G}}_j\left(z\right)=D_j^{-1}\left(z\right){\mathit{N}}_j\left(z\right)=\left[\begin{array}{cccc}{G}_{j1}\left(z\right)& G_{j2}\left(z\right)& .& G_{jp}\left(z\right)\end{array}\right]$; ${\upsilon }_j\left(k\right)={\displaystyle \sum _{i=1}^p{\upsilon }_{ij}\left(k\right)}$.

Expressing the frequency‐domain model (10) in a linear regression form yields

$y_j(k)={\mathit{\psi }}_j^T\left(k\right){\mathit{\theta }}_j+{\upsilon }_j\left(k\right);j=1,2,3,\mathrm{\dots },q$E11

where ${\mathit{\psi }}_j^T(k)$ is $1x(n_j+n_{j}p)$ regression vector formed of the regression vectors ${\mathit{\psi }}_{yj}^T(k)$ associated with $y_j(k)$, and ${\mathit{\psi }}_r^T(k)$ associated with $\mathit{r}\left(k\right)$:

${\mathit{\psi }}_j^T(k)=\left[\begin{array}{cc}{\mathit{\psi }}_{yj}^T(k)& {\mathit{\psi }}_r^T(k)\end{array}\right]$E12

${\mathit{\psi }}_{yj}^T(k)=\left[\begin{array}{cccc}-y_j(k-1)& -y_j(k-2)& .& -y_j(k-n_j)\end{array}\right]$; ${\mathit{\psi }}_r^T(k)=\left[\begin{array}{cccc}{\mathit{\psi }}_{r1}^T(k)& {\mathit{\psi }}_{r2}^T(k)& .& {\mathit{\psi }}_{rp}^T(k)\end{array}\right]$; ${\mathit{\theta }}_j$ is $(n_j+n_{j}p)\mathrm{x}1$ feature vector formed of the $n$ coefficients of the denominator polynomial $D_j\left(z\right)$ and the $n_{j}p$ coefficients of the numerator polynomial ${\mathit{N}}_j\left(z\right)$;

${\mathit{\theta }}_j={\left[\begin{array}{cc}{\mathit{\theta }}_{yj}& {\mathit{\theta }}_{rj}\end{array}\right]}^T$E13

where ${\mathit{\theta }}_{yj}={\left[\begin{array}{cccc}{a}_{j1}& a_{j2}& .& a_{j{n}_j}\end{array}\right]}^T$; ${\mathit{\theta }}_{rj}={\left[\begin{array}{cccc}{\mathit{\theta }}_{rj1}^T& {\mathit{\theta }}_{rj2}^T& .& {\mathit{\theta }}_{rjp}^T\end{array}\right]}^T$.

2.3. Multi‐input and multiple‐output system

Extending the results of the time‐domain expression to the MIMO (3), we get

$\mathit{y}(k)={\mathit{\psi }}^T\left(k\right)\mathit{\theta }+\mathit{\upsilon }\left(k\right)$E14

where ${\mathit{\psi }}^T\left(k\right)$ is $qx(n+npq)$ regression matrix formed of the regression vectors $\left\{{\mathit{\psi }}_{ij}^T(k)\right\}$, and $\mathit{\theta }$ is $(n+npq)\mathrm{\text{x}}1$ feature vector formed of ${\mathit{\theta }}_j$, $j=1,2,\mathrm{\dots },q$ is given as follows:

${\mathit{\psi }}^T\left(k\right)=\left[\begin{array}{cccccc}{\mathit{\psi }}_{y1}^T(k)& {\mathit{\psi }}_r^T(k)& \mathbf{0}& \mathbf{0}& .& \mathbf{0}\\ {\mathit{\psi }}_{y2}^T(k)& \mathbf{0}& {\mathit{\psi }}_r^T(k)& \mathbf{0}& .& \mathbf{0}\\ {\mathit{\psi }}_{y3}^T(k)& \mathbf{0}& \mathbf{0}& {\mathit{\psi }}_r^T(k)& .& \mathbf{0}\\ .& .& .& .& .& .\\ {\mathit{\psi }}_{yq}^T(k)& \mathbf{0}& \mathbf{0}& \mathbf{0}& .& {\mathit{\psi }}_r^T(k)\end{array}\right];\mathit{\theta }=\left[\begin{array}{c}{\mathit{\theta }}_y\\ {\mathit{\theta }}_{r1}\\ {\mathit{\theta }}_{r2}\\ .\\ {\mathit{\theta }}_{rp}\end{array}\right]$E15

The regression model (14) is the time‐domain version of the frequency‐domain model (3). Expressing the time‐domain model (14) in the frequency domain, we get

$\mathit{y}(z)={\mathit{\psi }}^T\left(z\right)\mathit{\theta }+\mathit{\upsilon }\left(z\right)$E16

2.4. Interconnected system

The system is an interconnection of subsystems such as the plant, the actuator, the sensors, and the controllers shown in Figure 1. Subfigure A at the top shows that jth output of the system $y_j={\displaystyle \sum _{i=1}^p{y}_{ji}}$ is given by Eq. (10) where $y_{ij}\left(z\right)$ given in Eq. (5) is the output generated by the input $r_i$ acting alone.

Figure 1.

Pairing of the inputs and an output and the subsystem in the path $ji$.

Subfigure B at the bottom shows that the transfer function $G_{ji}\left(z\right)$ in the path from the input $r_i$ to the output $y_{ij}$ is formed of subsystems $\{G_{ijl}(z)\}$. The subsystem $G_{ijl}\left(z\right)$ is driven by the input $u_{jil}\left(z\right)$ and its output is corrupted by the disturbance $w_{jil}\left(z\right)$. The input and the output of $G_{ji}\left(z\right)$ are $r_i$ and $y_{ji}$_, respectively, $v_{ji}$ is the measurement noise, ${\vartheta }_{ji}$ given in Eq. (5) is the combined effect of the disturbances $\left\{w_{jik}\right\}$ and $\left\{v_{ji}\right\}$ on the output $y_{ji}(z)$.

2.5. Modeling of faults

There are two types of fault models, namely the additive and the multiplicative (or parametric) types. In the additive type, a fault is modeled as an additive exogenous input to the system, whereas in the multiplicative type, a fault is modeled as a change in the parameters, which completely characterize the fault behavior of the subsystems. Although the multiplicative and additive perturbation models are equivalent, the multiplicative‐type perturbation model is preferable. The multiplicative perturbation model of the cascade combination of subsystems can actually model the particular perturbation in any one of the subsystems under consideration.

3. Emulators

The emulator‐based identification scheme is motivated by the model‐free artificial neural network approach to capture the static and the dynamic behavior by presenting neural network data covering likely operating scenarios. An identified model at each operating point characterizes the behavior of the system in the neighborhood of that point. In practice, however, the system model may be perturbed because of variations in the parameters of the system. To overcome this problem, the system model is identified by performing a number of emulator parameter‐perturbed experiments proposed in [4–5]. Each experiment consists of perturbing one or more emulator parameters. A linear model, termed optimal model, is identified as a best fit to the input‐output data from the set of emulated perturbations. The optimal model thus obtained characterizes the behavior of the system over wider operating regions (in the neighborhood of the operating point), whereas the conventional model characterizes the behavior merely at the nominal operating point (i.e., the conventional approach assumes that the model of the system remains unperturbed at every operating point). The optimal model is more robust, that is, the identification errors resulting from the variations in the emulator parameters are significantly lower compared to those of the conventional one based on performing a single experiment (i.e., without using emulators).

During the system identification phase, a number of experiments are performed by (a) not perturbing the emulator parameters and (b) perturbing the emulator parameters one at a time, simultaneously perturbing two at a time, three at time, and so on till perturbing all of them. The input‐output data collected from all experiments are termed emulator‐generated data.

• Nominal system model and the Kalman filter: The emulator‐generated data are used to identify the nominal optimal model of the system and the optimal Kalman filter model using the prediction error method.

• Estimation of the influence vectors: Using the least‐squares method, the influence vectors are identified recursively using the input‐output data obtained from the emulator‐perturbed parameter experiments. First, the influence vector for influence vector for the single parameter perturbation is identified, and then using the estimated influence vector, the influence vector for the two simultaneous emulator perturbations is estimated. Generalizing, the influence vector for $m$ simultaneous perturbation is identified, and then using all previous $m$ estimates of the influence vectors, the $\left(m+1\right)\text{th}$ influence vector is identified.

The emulators are transfer functions, which are connected in cascade with the subsystems to generate likely operating scenarios including normal and faulty one for reliable and accurate identification of the system, its associated Kalman filter, and the feature vector‐emulator map.

Emulators are connected to the system during the identification phase and its parameter is varied to generate likely operating scenarios. During the operational phase, the static emulators are disconnected, as it were, by setting them to unit values. The dynamic emulator, however, is not disconnected. Its gain is set to unity and its phase made a non‐zero negligibly small value so that (a) both of these parameters have a negligible effect on the dynamic behavior of the system during the operational phase and (b) the order of the system during the identification and the operational phases remains identical to ensure mathematical tractability without causing performance degradation. The role of the emulator‐generated data includes the following:

3.1. Emulator‐generated data for MISO system

The MISO system is given by Eq. (11) relating all the inputs $\mathit{r}(k)$ and the output $y_j(k)$ identified by connecting an emulator $E_j(z)$ in cascade with $\mathit{r}(z)$. The emulator is a first‐order all‐pass filter given by

$E_j\left(z\right)={\gamma }_{j2}\left(\frac{{\gamma }_{j1}+z^{-1}}{1+{\gamma }_{j1}{z}^{-1}}\right)$E17

where $\left|{\gamma }_{j1}\right|<1$ to ensure stability. The emulators ${\gamma }_{j1}$ and ${\gamma }_{j2}$ are varied one at a time, and both simultaneously. During the identification, an emulator $E_j(z)$, which is a first‐order all‐pass filter (17), is connected to the input $r_j(k)$ in cascade with nominal model ${\mathit{G}}_{j0}(z)$. A number of experiments are performed by varying the emulator parameters ${\gamma }_{j1}$, ${\gamma }_{j2}$ one at a time and both simultaneously to acquire emulator‐generated data: it is assumed for simplicity that the same input is applied to all the experiments. Using Eq. (10), the MISO model relating $r_j(k)$ and $y_j(k)$ becomes

$y_j^{el}(z)={\mathit{G}}_{j0}(z)E_j(z)\mathit{r}(z)+{\vartheta }_j^e(z),el=1,2,\mathrm{\dots },n_{\exp },l=1,2,3$E18

where $y_j^{e1}(z)$, $y_j^{e2}(z)$, and $y_j^{e3}(z)$ denote, respectively, the output generated by varying ${\gamma }_{j1}$, ${\gamma }_{j2}$ and both ${\gamma }_{j1}$, ${\gamma }_{j2}$.

3.2. Emulator‐generated data for SISO system

The feature vector‐emulator map of the SISO system (5) is estimated for the isolation of faults in the subsystems $\left\{G_{ij\ell }\left(z\right)\right\}$. The emulators $E_{ji\ell }\left(z\right)$ are connected to an accessible input or output $\left\{u_{ji\ell }\right\}$ in cascade with the subsystems $\left\{G_{ji\ell }(z)\right\}$ to mimic their variations. In other words, the known emulator parameter variations mimic those of the unknown parameters of the associated subsystems. The accessible inputs include the tracking error, the control input, actuator input, and sensor output.

The emulator $E_{ji\ell }\left(z\right)$ may be a dynamic system, a constant gain $\left({\gamma }_{ji\ell }\right)$, a gain, and a pure delay of $d$ time instants $\left({\gamma }_{ji\ell }{z}^{-d}\right)$, a first‐order all‐pass filter $\left(\frac{{\gamma }_{ji\ell }+z^{-1}}{1+{\gamma }_{ji\ell }{z}^{-1}}\right)$ or a Blaschke product of all first‐order‐pass filters $\left({\displaystyle \prod _{\ell }\frac{{\gamma }_{ji\ell }+z^{-1}}{1+{\gamma }_{ji\ell }{z}^{-1}}}\right)$ [3]. The emulator $E_{ji}\left(z\right)$ is chosen to be a product of a static gain and a first‐order all‐pass filter to mimic the behavior of the subsystem $G_{ji}\left(z\right)={\displaystyle \prod _{\ell =1}^l{G}_{ji\ell }\left(z\right)}$ of the SISO system given by Eqs. (5) and (6)

${\gamma }_{ji2}\left(\frac{{\gamma }_{ji1}+z^{-1}}{1+{\gamma }_{ji1}{z}^{-1}}\right)$E19

In order to ensure stability of the dynamic emulator, parameter ${\gamma }_{ji1}$ is constrained by $\left|{\gamma }_{ji\ell }\right|<1$.

Connecting the emulator $E_{ji}\left(z\right)$ given in Eq. (19) to the nominal SISO model $G_{ji0}\left(z\right)$ using Eqs. (5) and (6), we get

$y_{ji}^{el}\left(z\right)=E_{ji}\left(z\right)G_{ji0}\left(z\right)r_i\left(z\right)+{\vartheta }_{ji}^e(z),e=1,2,\mathrm{\dots },n_{\exp },l=1,2,3$E20

where $y_{ji}^{e1}\left(z\right)$, $y_{ji}^{e2}\left(z\right)$_, and $y_{ji}^{e3}\left(z\right)$ denote, respectively, the output generated by varying ${\gamma }_{ji1}$, ${\gamma }_{ji2}$ and both ${\gamma }_{ji1}$ and ${\gamma }_{ji2}$.

Figure 2 shows an example of a closed‐loop position control system formed of a controller, an actuator, a plant, and a sensor in the path connecting the tracking error $e_{ri}\left(k\right)=r_i\left(k\right)-y_{ji}\left(k\right)$ and the output $y_{ji}$. Only $e_{ri}\left(k\right)$, $u_{ij1}\left(k\right)$, and $u_{ij3}\left(k\right)$ are the measurement outputs. The emulators $E_{ji1}\left(z\right)=\frac{{\gamma }_{ji1}+z^{-1}}{1+{\gamma }_{ji1}{z}^{-1}}$, and $E_{ji2}={\gamma }_{ji2}$ are connected to $u_{ji1}$, and $E_{ji3}={\gamma }_{ji3}$ is connected to $u_{ji3}$ to mimic the perturbations in the dynamic plant $G_{ji1}\left(z\right)$, the static actuator $G_{ji2}\left(z\right)=k_A$ and the static sensor $G_{ji3}\left(z\right)=k_s$, respectively, where $E_{ji1}(z)$ is dynamic, and $E_{ji2}$ and $E_{ji3}$ are static emulators.

Figure 2.

Position control system: emulators and subsystems.

The nominal static emulator is set to unit value ${\gamma }_{ijk}^0=1$. The variation $\text{Δ}{\gamma }_{jik}$ of an emulator ${\gamma }_{ijk}$ may be expressed in terms of its nominal value ${\gamma }_{jik}^0$ as $\text{Δ}{\gamma }_{jik}={\gamma }_{jik}-{\gamma }_{jik}^0$.

3.3. Feature vector‐emulator map

The feature vector‐emulator map for the SISO and the MISO systems is developed subsequently.

3.3.1. SISO system

Consider the emulator‐perturbed SISO system (20) relating the inputs $r_i\left(k\right)$ and $y_{ji}\left(k\right)$ and the associated linear regression model (7). The feature vector ${\mathit{\theta }}_{ji}$ is a nonlinear function of the emulator parameter ${\mathit{\gamma }}^{ji}=\left[\begin{array}{cc}{\gamma }_{ji1}& {\gamma }_{ji2}\end{array}\right]$. Assuming that the feature vector ${\mathit{\theta }}_{ji}$ is a continuous function of ${\mathit{\gamma }}^{ji}$, then using Weierstrass approximation theorem, the feature vector‐emulator map becomes

$\text{Δ}{\mathit{\theta }}_{ji}={\mathbf{\Omega }}_{ji1}\text{Δ}{\gamma }_{ji1}+{\mathbf{\Omega }}_{ji2}\text{Δ}{\gamma }_{ji2}+{\mathbf{\Omega }}_{ji12}\text{Δ}{\gamma }_{ji1}\text{Δ}{\gamma }_{ji2}$E21

where $\text{Δ}{\mathit{\theta }}_{ji}={\mathit{\theta }}_{ji}-{\mathit{\theta }}_{ji}^0$;$\text{Δ}{\gamma }_{j\ell }={\gamma }_{j\ell }-{\gamma }_{j\ell }^0$ is the parameter variation;${\mathit{\theta }}_{ji}^0$ is the nominal feature vector; ${\mathbf{\Omega }}_{ji1}$ is a $2n_{j}x1$ vector of partial derivatives of the feature vector ${\mathit{\theta }}_{ji}$ with respect to ${\gamma }_{ji1}$ evaluated at the unperturbed nominal emulator value ${\gamma }_{ji1}^0$. Similarly, ${\mathbf{\Omega }}_{ji2}$ is a $2n_{j}x1$ vector of partial derivatives of the feature vector ${\mathit{\theta }}_{ji}$ with respect to ${\gamma }_{ji2}$ evaluated at the unperturbed nominal emulator value ${\gamma }_{ji2}^0$, ${\mathbf{\Omega }}_{ji12}$ is the second partial derivatives with respect to ${\gamma }_{ji1}$ and ${\gamma }_{ji2}$ evaluated at the unperturbed nominal emulator value ${\gamma }_{ji1}^0$ and ${\gamma }_{ji2}^0$. The partial derivative terms ${\mathbf{\Omega }}_{ji1}$, ${\mathbf{\Omega }}_{ji2}$ ${\mathbf{\Omega }}_{ji12}$, which are the Jacobean of the feature vector ${\mathit{\theta }}_{ji}$ with respect to the emulator parameters $\left\{{\gamma }_{jik}\right\}$, are termed influence vectors. The influence vectors play a crucial role in isolating a fault occurring in any subsystem. The influence vectors ${\mathbf{\Omega }}_{ji1}$, ${\mathbf{\Omega }}_{ji2}$_, and ${\mathbf{\Omega }}_{ji12}$ track the degree of variations in the parameters of the subsystem perturbations.

Substituting for ${\mathit{\theta }}_{ji}$ in (7), the variation $\text{Δ}{y}_{ji}(k)=y_{ji}(k)-y_{ji}^0(k)$ between the actual output $y_{ji}(k)$ and the nominal fault‐free output $y_{ji}^0(k)$ becomes

$\text{Δ}{y}_{ji}(k)={\mathit{\psi }}_{ji}^T(k)\left({\mathbf{\Omega }}_{ji1}\text{Δ}{\gamma }_{ji1}+{\mathbf{\Omega }}_{ji2}\text{Δ}{\gamma }_{ji2}+{\mathbf{\Omega }}_{ji12}\text{Δ}{\gamma }_{ji1}\text{Δ}{\gamma }_{ji2}\right)+{\upsilon }_{ji}(k)$E22

Let ${\mathbf{\Omega }}_{ji}$ be an influence matrix associated with the emulators located at the path $ij$

${\mathbf{\Omega }}_{ji}=\left[\begin{array}{ccccc}{\mathbf{\Omega }}_{jik}& {\mathbf{\Omega }}_{jik\ell }& {\mathbf{\Omega }}_{jik\ell mn}& .& {\mathbf{\Omega }}_{ji\mathrm{12\dots }q}\end{array}\right]$E23

A number of emulator parameter‐perturbed experiments are performed by perturbing the parameters of the emulators (20). For each experiment, $N$ input‐output data $({\mathit{y}}_j^e(k),\mathit{r}(k))$ are obtained, $k=1,2,\mathrm{\dots },N$. The input $\mathit{r}(k)$ for each experiment is chosen to be persistently exciting. The regression models associated with the experiments and Eq. (22) are given as follows:

$\begin{array}{l}\text{Δ}{y}_{ji}^{e1}(k)={\psi }_{ji}^T(k)\text{Δ}{\gamma }_{ji1}{\mathbf{\Omega }}_{ji1}+{\mathit{\upsilon }}_{ji}^{e1}(k)\\ \Lambda y_{ji}^{e2}(k)-{\psi }_{ji}^T(k)\text{Δ}{\gamma }_{ji1}{\mathbf{\Omega }}_{ji1}={\psi }_{ji}^T(k){\mathbf{\Omega }}_{ji2}\text{Δ}{\gamma }_{ji2}+{\upsilon }_{ji}^{e2}(k)\\ \text{Δ}{y}_{ji}^{e3}(k)-{\psi }_{ji}^T(k)\left(\text{Δ}{\gamma }_{ji1}{\mathbf{\Omega }}_{ji1}+{\mathbf{\Omega }}_{ji2}\text{Δ}{\gamma }_{ji2}\right)={\psi }_{ji}^T(k)\left({\mathbf{\Omega }}_{ji12}\text{Δ}{\gamma }_{ji12}\right)+{\upsilon }_{ji}^{e3}(k)\end{array}$E24

3.3.2. MISO system

Consider the emulator‐perturbed MISO system (18) relating the inputs $\mathit{r}\left(k\right)$ and $y_j\left(k\right)$, and the associated linear regression model (11). Similar to Eqs. (21) and (24), we get

$\text{Δ}{\mathit{\theta }}_j={\mathbf{\Omega }}_{j1}\text{Δ}{\gamma }_{j1}+{\mathbf{\Omega }}_{j2}\text{Δ}{\gamma }_{j2}+{\mathbf{\Omega }}_{j12}\text{Δ}{\gamma }_{j1}\text{Δ}{\gamma }_{j2}$E25

$\begin{array}{l}\text{Δ}{y}_j^{e1}(k)={\psi }_j^T(k)\text{Δ}{\gamma }_{j1}{\mathbf{\Omega }}_{j1}+{\mathit{\upsilon }}_j^{e1}(k)\\ \Lambda y_j^{e2}(k)-{\psi }_j^T(k)\text{Δ}{\gamma }_{j1}{\mathbf{\Omega }}_{j1}={\psi }_j^T(k){\mathbf{\Omega }}_{j2}\text{Δ}{\gamma }_{j2}+{\upsilon }_j^{e2}(k)\\ \text{Δ}{y}_j^{e3}(k)-{\psi }_j^T(k)\left(\text{Δ}{\gamma }_{j1}{\mathbf{\Omega }}_{j1}+{\mathbf{\Omega }}_{j2}\text{Δ}{\gamma }_{j2}\right)={\psi }_j^T(k)\left({\mathbf{\Omega }}_{j12}\text{Δ}{\gamma }_{j12}\right)+{\upsilon }_j^{e3}(k)\end{array}$E26

4. Identification

The prediction error method can be derived from the residual model of the Kalman filter, which is presented in the next section. It is used to identify both the nominal system and the Kalman filter associated with the system without the need for a priori knowledge of the covariances of the noise and the disturbance statistics. Prediction error method is consistent, efficient, and a gold standard for system identification, and can identify open‐loop and closed‐loop systems. The variance the parameter estimates asymptotically approaches the Cramer‐Rao lower bound.

Optimal models: The optimal system and the associated Kalman filter are identified using the prediction error method using computationally efficient scheme. First, the MISO system is identified and then the SISO system is derived from the estimate of feature vector associated with the MISO system. The emulator‐generated data generated using Eq. (18) are used to identify MISO system (10) and the nominal feature vector ${\mathit{\theta }}_j^0$ for Eq. (11), which is the best least‐squared fit to set all perturbed feature vector ${\mathit{\theta }}_j$, and the Kalman gain ${\mathit{K}}_{j0}$ are estimated. Let the optimal state‐space model of the MISO system be $\left({\mathit{A}}_{j0},{\mathit{B}}_{j0},{\mathit{C}}_{j0}\right)$ and associated Kalman filter be $\left({\mathit{A}}_{j0}-{\mathit{K}}_{j0}{\mathit{C}}_{j0},\left[{\mathit{K}}_{j0}{\mathit{B}}_{j0}\right],{\mathit{C}}_{j0}\right)$. Let the optimal transfer matrix of the MISO system and the optimal estimate of the output be ${\mathit{G}}_j^{opt}\left(z\right)$ and ${\widehat{y}}_j^{opt}(k)$_, respectively. Using Eq. (10), we get

${\widehat{y}}_j^{opt}\left(z\right)={\mathit{G}}_j^{opt}\left(z\right)\mathit{r}\left(z\right)+{\vartheta }_j\left(z\right)$E27

Then, the best estimate of the feature vector ${\mathit{\theta }}_{ji}$ of the SISO system (7), denoted ${\mathit{\theta }}_{ji}^0$_, and the Kalman gain are estimated from ${\mathit{\theta }}_j^0$.

4.1. Estimation of the influence vectors

SISO system: Knowing the emulator parameter perturbations $\text{Δ}{\gamma }_{ji1}$, $\text{Δ}{\gamma }_{ji2}$, $\text{Δ}{\gamma }_{ji12}$ and the resulting emulator‐generated data, the influence vectors ${\widehat{\mathbf{\Omega }}}_{ji1}$, ${\widehat{\mathbf{\Omega }}}_{ji2}$, and ${\widehat{\mathbf{\Omega }}}_{ji12}$ are estimated recursively using the least‐squares method using Eq. (24)

$\begin{array}{l}{\widehat{\mathbf{\Omega }}}_{ji1}=\text{arg}\underset{{\mathbf{\Omega }}_{ji1}}{\min }\left\{{\Vert \text{Δ}{y}_{ji}^{ei}(k)-{\psi }_{ji}^T(k){\mathbf{\Omega }}_{ji1}\text{Δ}{\gamma }_{ji1}\Vert }^2\right\}\\ {\widehat{\mathbf{\Omega }}}_{ji2}=\text{arg}\underset{{\mathbf{\Omega }}_{ji2}}{\min }\left\{{\Vert \mathrm{\Lambda }{y}_{ji}^{e2}(k)-{\psi }_{ji}^T(k)\text{Δ}{\gamma }_{ji1}{\widehat{\mathbf{\Omega }}}_{ji1}-{\psi }_{ji}^T(k){\mathbf{\Omega }}_{ji2}\text{Δ}{\gamma }_{ji2}\Vert }^2\right\}\\ {\widehat{\mathbf{\Omega }}}_{ji12}=\text{arg}\underset{{\mathbf{\Omega }}_{jiklm}}{\min }\left\{{\Vert \text{Δ}{y}_j^{e3}(k)-{\psi }_{ji}^T(k)\left(\text{Δ}{\gamma }_{ji1}{\widehat{\mathbf{\Omega }}}_{ji1}+{\widehat{\mathbf{\Omega }}}_{ji2}\text{Δ}{\gamma }_{ji2}\right)-{\psi }_{ji}^T(k){\mathbf{\Omega }}_{ji12}\text{Δ}{\gamma }_{ji12}\Vert }^2\right\}\end{array}$E28

where ${\Vert x\left(k\right)\Vert }^2={\displaystyle \sum _{k=1}^N{x}^2\left(k\right)}$.

MISO system: Similar to Eq. (28), the influence vectors ${\widehat{\mathbf{\Omega }}}_{jk}$, ${\widehat{\mathbf{\Omega }}}_{j2}$, and ${\widehat{\mathbf{\Omega }}}_{j12}$ are estimated.

5. Model of the Kalman filter

The Kalman filter forms the backbone of the MISO and the SISO systems fault detection and for fault isolation, respectively. The Kalman filter is a closed‐loop system, which is (a) an exact copy of the identified nominal of the system driven by the residual, which is the error between the output and its estimate, and (b) is stabilized by the Kalman gain.

MISO system: Using the state‐space model $\left({\mathit{A}}_{j0},{\mathit{B}}_{j0},{\mathit{C}}_{j0}\right)$ derived from the identified nominal feature vector ${\mathit{\theta }}_j^0$. The Kalman filter $\left({\mathit{A}}_{j0}-{\mathit{K}}_{j0}{\mathit{C}}_{j0},\left[{\mathit{K}}_{j0}{\mathit{B}}_{j0}\right],{\mathit{C}}_{j0}\right)$ associated with the MISO system (10) is

$\begin{array}{l}{\widehat{\mathit{x}}}_j(k+1)=\left({\mathit{A}}_{j0}-{\mathit{K}}_{j0}{\mathit{C}}_{j0}\right){\widehat{\mathit{x}}}_j(k)+{\mathit{K}}_{j0}{y}_j\left(k\right)+{\mathit{B}}_{j0}\text{}\mathit{r}(k)\\ {\widehat{y}}_j(k)={\mathit{C}}_{j0}{\widehat{\mathit{x}}}_j(k)\\ e_j(k)=y_j(k)-{\widehat{y}}_j(k)\end{array}$E29

where ${\widehat{\mathit{x}}}_j(k)$ and ${\widehat{y}}_j(k)$ are, respectively, the minimum variance estimates of the state and the output.

Figure 3 shows the nominal fault‐free system and the Kalman filter. The structure of the Kalman filter is based on the internal model principle, which embodies the nominal system model $\left({\mathit{A}}_{j0},{\mathit{B}}_{j0},{\mathit{C}}_{j0}\right)$. The inputs to the Kalman filter are the input $\mathit{r}(k)$ and the output $y_j(k)$ which is corrupted by the disturbance $w_j(k)$ and the measurement noise $v_j(k)$.

Figure 3.

The system and its associated Kalman filter.

5.1. Expressions of the residual

The expression for the residuals for the MISO system $e_j(z)$ and the SISO system $e_{ji}(z)$ is derived from the Kalman filter (29).

MISO model: The frequency‐domain expression, relating the $n_{u}x1$ input $\mathit{r}(z)$ and output $y_j(z)$ to the residual $e_j(z)$ is given by the following model, termed residual model:

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}{y}_j(z)-\frac{{\overline{\mathit{N}}}_{j0}(z)}{F_{j0}(z)}\mathit{r}(z)$E30

where $F_{j0}(z)=\left|z\mathit{I}-{\mathit{A}}_{j0}+{\mathit{K}}_{j0}{\mathit{C}}_{j0}\right|$ is the characteristic polynomial termed Kalman polynomial; ${\overline{D}}_{j0}\left(z\right)=F_{j0}(z)\left(\mathit{I}-{\mathit{C}}_{j0}{\left(\mathit{z}\mathit{I}-{\mathit{A}}_{j0}+{\mathit{K}}_{j0}{\mathit{C}}_{j0}\right)}^{-1}{\mathit{K}}_{j0}\right)$${\overline{\mathit{N}}}_{j0}\left(z\right)=\left[\begin{array}{cccc}{\overline{N}}_{j10}\left(z\right)& {\overline{N}}_{j10}\left(z\right)& .& {\overline{N}}_{jp0}\left(z\right)\end{array}\right]=F_{j0}(z)\left({\mathit{C}}_{j0}{\left(z\mathit{I}-{\mathit{A}}_{j0}+{\mathit{K}}_{j0}{\mathit{C}}_0\right)}^{-1}{\mathit{B}}_{j0}\right)$

SISO system: The residual $e_{ji}(z)$ is derived from the residual model (30) from the map relating $e_{ji}(z)$ to $y_j(z)$ and $r_i(z)$:

$e_{ji}\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}{y}_j(z)-\frac{{\overline{N}}_{ji0}(z)}{F_{j0}(z)}{r}_i(z)$E31

where ${\overline{N}}_{0ji}(z)$ is the $i^{\text{th}}$ element of ${\overline{\mathit{N}}}_{j0}\left(z\right)$.

5.1.1. Key properties of the Kalman filter residual

The Kalman filter forms the backbone of the proposed scheme in view of its key properties proved in [1]. These properties exploited in developing the system identification using the residual model, and in unified approach to fault detection and isolation fault, where a fault is defined as an incipient fault resulting in the model mismatch.

5.2. Propositions

We establish important results, in the form of lemmas that are crucial to the development of the proposed fault diagnosis scheme. In Lemma 1, it is shown that (a) the system transfer function can be estimated from the residual model and (b) Kalman filter whitens the output error ${\vartheta }_j\left(z\right)$ given in Eq. (10). Lemma 2 shows that the residual will not be a zero‐mean white noise process if there is a model mismatch, and there will be an additive fault indicating term, which is a function of the deviation between the actual feature vector ${\mathit{\theta }}_j$ of the system model $\left({\mathit{A}}_j,{\mathit{B}}_j,{\mathit{C}}_j\right)$ and the nominal fault‐free feature vector ${\mathit{\theta }}_j^0$ of nominal fault‐free model $\left({\mathit{A}}_{j0},{\mathit{B}}_{j0},{\mathit{C}}_{j0}\right)$.

Case 1: The system and the nominal models are identical

Lemma 1:

${\mathit{G}}_{j0}(z)=D_{j0}^{-1}\left(z\right){\mathit{N}}_{j0}(z)={\overline{D}}_{j0}^{-1}\left(z\right){\overline{\mathit{N}}}_{j0}(z)$E32

where ${\mathit{G}}_{j0}(z)$ is the transfer function of the nominal fault‐free model $\left({\mathit{A}}_{j0},{\mathit{B}}_{j0},{\mathit{C}}_{j0}\right)$.

Proof: Substituting for $\mathit{y}\left(z\right)$ from Eq. (10), the residual model (30) becomes

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}\left(\left(D_j^{-1}(z){\mathit{N}}_j\left(z\right)-{\overline{D}}_{j0}^{-1}\left(z\right){\overline{\mathit{N}}}_{j0}(z)\right)\mathit{r}(z)\right)+\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}{\vartheta }_j\left(z\right)$E33

Correlating both sides with input $\mathit{r}(k)$, and invoking the orthogonality properties, the residual, namely $\mathit{r}(k)$_, is uncorrelated with both $e_j(k)$ and the output error ${\vartheta }_j(k)$ [4], we get

$\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}\left(\left(D_j^{-1}(z){\mathit{N}}_j\left(z\right)-{\overline{D}}_{j0}^{-1}\left(z\right){\overline{\mathit{N}}}_{j0}(z)\right)\mathit{r}(z)\right)=0$E34

Hence, Eq. (32) holds.

Corollary 1: The filter $\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}$ whitens the output error ${\vartheta }_j\left(z\right)$ if there is no model mismatch:

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}{\vartheta }_j\left(z\right)$E35

Proof: Consider the expression for the model‐matching case (33). Using Eq. (32), we establish Eq. (35).

Case 2: System and the nominal model mismatch

Lemma 2: If there is model mismatch, then

$e_j(z)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}\text{Δ}{G}_j\left(z\right)+{\vartheta }_{jf}$E36

$e_j(z)={\mathit{\psi }}_{jf}^T(z)\text{Δ}{\mathit{\theta }}_j+{\upsilon }_{jf}(z)$E37

where $\text{Δ}{G}_j\left(z\right)=D_j^{-1}\left(z\right){\mathit{N}}_j(z)-D_{j0}^{-1}\left(z\right){\mathit{N}}_{j0}(z)$, $\text{Δ}{\mathit{\theta }}_j={\mathit{\theta }}_j-{\mathit{\theta }}_j^0$; ${\mathit{\psi }}_{jf}^T(z)=\frac{{\overline{D}}_{j0}\left(z\right)}{D_j(z)F_{j0}(z)}{\mathit{\psi }}_j^T(z)$, ${\vartheta }_{jf}=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}{\vartheta }_j(z)$, and $e_{jf}(z)=\frac{{\overline{D}}_{j0}\left(z\right)}{D_j(z)F_{j0}(z)}{\upsilon }_j\left(z\right)$ are the filtered regression matrix ${\mathit{\psi }}_j^T(z)$ and filtered output error ${\vartheta }_j(z)$, filtered equation error ${\upsilon }_j\left(z\right)$_, respectively.

Proof:

Case 1: Consider expression (33). Using Eq. (32), we get

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}\left(\left(D_j^{-1}(z){\mathit{N}}_j\left(z\right)-D_{j0}^{-1}\left(z\right){\mathit{N}}_{j0}(z)\right)\mathit{r}(z)\right)+{\vartheta }_{jf}\left(z\right)$E38

Substituting $\text{Δ}{G}_j\left(z\right)=D_j^{-1}\left(z\right){\mathit{N}}_j(z)-D_{j0}^{-1}\left(z\right){\mathit{N}}_{j0}(z)$, we get Eq. (36).

Case 2: Expressing the residual model (30) in an alternative form:

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{F_{j0}(z)}\left(y_j(z)-{\overline{D}}_{j0}^{-1}\left(z\right){\overline{\mathit{N}}}_{j0}(z)\mathit{r}(z)\right)$E39

Using Eq. (32) and re‐arranging, we get

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{D_{j0}\left(z\right)F_{j0}(z)}\left(D_{j0}\left(z\right)y_j(z)-{\overline{\mathit{N}}}_{j0}(z)\mathit{r}(z)\right)$E40

Adding and subtracting $y_j(z)$ inside the bracket on the right‐hand side yields

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{D_{j0}\left(z\right)F_{j0}(z)}\left(y_j(z)-\left(1-D_{j0}\left(z\right)\right)y_j(z)-{\overline{\mathit{N}}}_{j0}(z)\mathit{r}(z)\right)$E41

Using the expression for the regression model (11) and substituting for the actual and the nominal fault‐free cases, we get

$e_j\left(z\right)=\frac{{\overline{D}}_{j0}\left(z\right)}{D_{j0}\left(z\right)F_{j0}(z)}{\mathit{\psi }}_j^T\left(k\right)\text{Δ}{\mathit{\theta }}_j+{\upsilon }_{jf}\left(k\right)$E42

Remarks: If there is a model mismatch because of variations in the subsystem parameters, the residual is no longer zero‐mean white noise process. The residual has an additive term, which is affine in the deviation in the system transfer function $\text{Δ}{G}_j\left(z\right)$ or equivalently affine in the feature vector ${\mathit{\psi }}_{jf}^T(z)\text{Δ}{\mathit{\theta }}_j$. The additive terms are termed fault indicators. This shows that the Kalman filter provides a unifying approach to handle both fault detection and fault isolation.

In view of the key properties, the Kalman filter is employed for identification and the fault diagnosis. In system identification, the criterion for determining whether the identified model has captured completely the dynamic behavior of the system is that the residual (error between the output and its estimate obtained using the identified model) is a zero‐mean white noise process. Consider the problem of identification of the system. Since the equation error $\mathit{\upsilon }\left(k\right)$ is a colored noise process, the parameter estimates will be biased and inefficient. To overcome this, the input and the output are whitened using the Kalman filter as shown in Eq. (35) of Corollary 1. The Kalman filter model (29) may be interpreted as an inverse system generating the innovation sequence $\mathit{e}(k)$_, or alternatively as a whitening‐filter realization of a state‐space model that is driven by both the disturbance and measurement noise.

Lemma 3

$e_{ji}\left(z\right)={\mathit{\psi }}_{jif}^T\left(z\right)\text{Δ}{\mathit{\theta }}_{ji}+{\upsilon }_{jif}\left(z\right)$E43

where $F_{ji0}(z)=\left|z\mathit{I}-{\mathit{A}}_{ji0}+{\mathit{K}}_{ji0}{\mathit{C}}_{j0}\right|$, ${\mathit{\psi }}_{jif}^T(z)=\frac{{\overline{D}}_{j0}\left(z\right)}{D_j(z)F_{ji0}(z)}{\mathit{\psi }}_{ji}^T(z)$;${\upsilon }_{jif}(z)=\frac{{\overline{D}}_{0j}\left(z\right)}{D_j(z)F_{ji0}(z)}{\upsilon }_{ji}\left(z\right)$

Proof: The proof follows from Eqs. (31) and (37).

6. Bayesian approach fault diagnosis

The objective of fault detection is to assert whether the given residual belongs to a set of fault‐free data or faulty residual data, while fault isolation is determined to which class of emulator‐perturbed residual the given data belong. The problem of fault detection and fault isolation is formulated by a pattern classification problem. Fault detection is a binary pattern classification, while the fault isolation is a multi‐class pattern classification. The Bayesian decision strategy is employed to assert appropriate class label. The Bayesian decision strategy is based on the a posteriori conditional probability of deciding a hypothesis given the data, a priori probability of the hypothesis, and a performance measure. The decision strategy is determined from the minimization of the performance measure with respect to all hypotheses.

The $Nx1$ residual $\mathit{e}(k)$ is located in a different region of the $N$‐dimensional plane depending upon the fault type. In the ideal case regions, there will not be overlaps between regions associated with different fault types. However, due to noise, disturbances, and other measurement artifacts there will be overlap between the various regions. Hence, Bayesian strategy is employed to asset an appropriate class label to ensure a high‐probability correct decision, and a low probability of false alarms.

6.1. Fault detection

Fault detection is posed as a binary hypothesis‐testing problem. The criterion to choose between the two hypotheses, namely the presence or an absence of a fault, is based on minimizing the Bayes risk, which quantifies the costs associated with correct and incorrect decisions. The $Nx1$ Kalman filter residual data $\mathit{e}(k)$ generated by Eq. (29) is employed. The minimization of the Bayes risk yields the likelihood ratio test. The decision between the two hypotheses is based on comparing the likelihood ratio, which is the ratio of the conditional probabilities under the two hypotheses, to a threshold value. The resulting binary composite hypothesis‐testing problem compares the test statistics of residual $e(k)$ with a threshold value $\eta$:

$t_s(e)\{\begin{array}{cc}\le \eta & \text{no}\text{fault}\\ >\eta & \text{fault}\end{array}$E44

The test statistics depends upon the input $r(k)$ that generates the residual $e(k)$ [4]:

$t_s(e)=\{\begin{array}{cc}\left|\frac{1}{N}{\displaystyle \sum _{i=k-N+1}^{k}e(i)}\right|& r(k)=\text{constant}\\ P_{ee}(f_0)& r(k)\text{is}\text{a}\text{sinusoid}\\ \frac{1}{N}{\displaystyle \sum _{i=k-N+1}^k{e}^2(i)}& r(k)\text{is}\text{an}\text{arbitrary}\text{signal}\end{array}$E45

6.1.1. Computationally efficient scheme

A computationally efficient scheme is employed here for the detection:

• The status of each of the MISO systems ${\mathit{G}}_j(z)$ relating all the inputs $\mathit{r}(z)$ and all the outputs $y_j(z)$ is evaluated for all $j=1,2,\mathrm{\dots },q$ using the binary hypothesis scheme (44). Using the test statistics of the residuals $e_j(k)$ given by Eq. (30) yields

  $t_s(e_j)\{\begin{array}{cc}\le {\eta }_j& \text{no}\text{fault}\\ >{\eta }_j& \text{fault}\end{array},j=1,2,3,\mathrm{\dots },q$E46

• If a fault is asserted in ${\mathit{G}}_j\left(z\right)$, then the status of each of the $p$ subsystems $G_{ji}(z)$ of the SISO system is asserted using the test statistics of the residuals $e_{ji}(k)$ (31):

$t_s(e_{ji})\{\begin{array}{cc}\le {\eta }_{ji}& \text{no}\text{fault}\\ >{\eta }_{ji}& \text{fault}\end{array},i=1,2,3,\mathrm{\dots },p$E47

Fault accommodation: If a fault is asserted, then the Kalman gain is adapted online, the system re‐identified, and the Kalman filter redesigned accordingly, thus the fault is accommodated and, in the extreme case, the system is shut down for safety reasons.

7. Evaluation on simulated system

The proposed emulator‐based system identification of the system, the associated Kalman filter, feature vector‐emulator map, and finally the fault diagnosis are illustrated using an example of a position control system formed of an actuator, a sensor, and a plant.

7.1. System model

A two‐input and two‐output system fault‐free system $\left({\mathit{A}}_0,{\mathit{B}}_0,{\mathit{C}}_0\right)$ given by Eq. (1) is considered where

${\mathit{A}}_0=\left[\begin{array}{cccc}0& -0.7& 0& 0\\ 1& 1.5& 0& 0\\ 0& 0& 0& -0.82\\ 0& 0& 1& 1.8\end{array}\right];{\mathit{B}}_0=\left[\begin{array}{cc}0.5& 1\\ 1& 0\\ 1& -0.3\\ 0& 1\end{array}\right];{\mathit{C}}_0=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]$E48

The nominal transfer matrix of the MIMO system (3) is

${\mathit{G}}_0\left(z\right)=\left[\begin{array}{cc}{G}_{11}\left(z\right)& G_{12}\left(z\right)\\ G_{21}\left(z\right)& G_{22}\left(z\right)\end{array}\right]=\left[\begin{array}{cc}\frac{1+z^{-1}}{1-1.5z^{-1}+0.7z^{-2}}& \frac{1}{1-1.5z^{-1}+0.7z^{-2}}\\ \frac{1}{1-1.8z^{-1}+0.82z^{-2}}& \frac{1-0.2z^{-1}}{1-1.8z^{-1}+0.82z^{-2}}\end{array}\right]$E49

The nominal MISO transfer matrix, ${\mathit{G}}_{j0}(z)=D_0^{-1}(z){\mathit{N}}_{j0}(z)$, $j=1,2$, of the system is

$\begin{array}{l}{D}_0(z)=\text{ 1 }-\text{3}\text{.3}{z}^{-1}\text{ +4}\text{.22}{z}^{-2}-\text{ 2}\text{.49}{z}^{-3}\text{ + 0}\text{.574}{z}^{-4}\\ {\mathit{N}}_0\left(z\right)=\left[\begin{array}{cc}{z}^{-1}-\text{1}\text{.3}{z}^{-2}& -\text{0}\text{.08}{z}^{-1}\text{ + 0}\text{.41}{z}^{-2}\\ z^{-1}-\text{1}\text{.8}{z}^{-2}& \text{1}\text{.15}{z}^{-1}-\text{ 0}{\text{.2z}}^{-1}\end{array}\right]\end{array}$E50

Figure 4a shows the emulator‐generated MISO output 1, $y_1^{e1}$ MISO output 2, $y_2^{e1}$, SISO output 11, $y_{11}^{e1}$, SISO output 12, $y_{12}^{e1}$, SISO output 21, $y_{21}^{e1}$ and SISO output 22, $y_{22}^{e1}$ given in Eqs. (18) and (20) resulting from the variations of the emulator parameters ${\gamma }_{j1}$ and ${\gamma }_{ji1}$_, respectively. Subfigures A and B show plots of the perturbed step responses $y_1^{e1}\left(k\right)$ and $y_2^{e1}\left(k\right)$ with respect to time, while subfigures C–F show plots of the perturbed outputs $y_{11}^{e1}\left(k\right)$, $y_{12}^{e1}\left(k\right)$, $y_{21}^{e1}\left(k\right)$_, and $y_{22}^{e1}(k)$ with respect to time. The outputs are in centimeters (cm) and the time is in seconds (s). The plots are generated when the emulator parameter ${\gamma }_{j1}$ is varied. The variations $\text{Δ}{\gamma }_{j1}$ are $\left\{\begin{array}{cccc}0.1& 0.5& 0.9& 1\end{array}\right\}$.

Figure 4.

(a) Emulator generated data and (b) performance of the identified model.

The mean‐squared error (or residual), namely the error between the output of the optimal model, denoted by ${\widehat{y}}_j^{\text{opt}}\left(k\right)$ and given by Eq. (27), and the perturbed outputs $y_j^{e1}\left(k\right)$ resulting from the variations of the emulator parameter ${\gamma }_{j1}$. The mean‐squared error, denoted ${\text{mse}}_j\left({\gamma }_{j1}\right)$, is computed as follows:

${\text{mse}}_j\left({\gamma }_{j1}\right)=\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left({\widehat{y}}_j^{\text{opt}}(k)-y^{e1}(k)\right)}^2}$E51

The conventional scheme identifies only the unperturbed nominal model. Let the identified model of the MISO system (10) be ${\widehat{\mathit{G}}}_{j0}\left(z\right)$, the estimated output be ${\widehat{y}}_{j0}^c(z)$. The mean‐squared error, denoted ${\text{mse}}_j^c\left({\gamma }_{j1}\right)$_, becomes

${\text{mse}}_j^c\left({\gamma }_{j1}\right)=\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left({\widehat{y}}_{j0}(k)-y^{e1}(k)\right)}^2}$E52

The mean‐squared errors ${\text{mse}}_j\left({\gamma }_{j1}\right)$ and ${\text{mse}}_j^c\left({\gamma }_{j1}\right)$ are plotted as functions of the emulator parameter perturbations $\text{Δ}{\gamma }_{j1}$. The mean‐squared profiles of both the proposed emulator‐based and the conventional identification schemes are shown in subfigures A and B of Figure 4b.

The identified state‐space model and the Kalman gain are

$\begin{array}{ll}{\widehat{A}}_0\hfill & \text{ =}\left[\begin{array}{l}\text{ 0}\text{.9843 }‐\text{0}\text{.1588 0}\text{.0213 0}\text{.0224}\hfill \\ \text{ 0}\text{.1572 0}\text{.9266 0}\text{.2317 }‐\text{0}\text{.2502}\hfill \\ ‐\text{0}\text{.0631}‐\text{0}\text{.3144 0}\text{.9090 }‐\text{0}\text{.3230}\hfill \\ \text{ 0}\text{.0171 0}\text{.0153 0}\text{.1659 0}\text{.7225}\hfill \end{array}\right],{\widehat{B}}_0=\left[\begin{array}{l}\text{ 0 0}\hfill \\ \text{ 0}\text{.1 0}\text{.1}\hfill \\ ‐\text{0}\text{.3 }‐\text{0}\text{.3}\hfill \\ \text{ 2}\text{.6 2}\text{.6}\hfill \end{array}\right]{10}^3,\hfill \\ {\widehat{C}}_0\hfill & =\left[\begin{array}{l}\text{ 21 }‐\text{564}\text{.2 }‐\text{245}\text{.6 49}\text{.5}\hfill \\ \text{ 1679}\text{.2 }‐\text{336}\text{.3 }‐\text{211}\text{.4 44}\text{.0}\hfill \end{array}\right]\hfill \end{array}$E65

The ranges of the mean‐squared errors ${\text{mse}}_j\left({\gamma }_{j1}\right)$ and ${\text{mse}}_j^c\left({\gamma }_{j1}\right)$ are given below:

$\begin{array}{l}\text{ 1}\text{.8390}\le {\text{mse}}_1\left({\gamma }_{j1}\right)\le \text{ 2}\text{.225 }\\ \text{ 0}\text{.0137}\le {\text{mse}}_1^c\left({\gamma }_{j1}\right)\le \text{7}\text{.1815}\end{array}$E53

$\begin{array}{l}\text{18}\text{.2224}\le {\text{mse}}_2\left({\gamma }_{j1}\right)\le \text{ 21}\text{.7167 }\\ \text{ 0}\text{.0007}\le {\text{mse}}_2^c\left({\gamma }_{j1}\right)\le 69.8841\end{array}$E54

Remarks: The emulator‐generated data cover the operating scenarios, including both the normal and abnormal ones, exhibiting variations of the rise time, the settling times, and the overshoots.

The identified optimal model $\left({\widehat{\mathit{A}}}_0,{\widehat{\mathit{B}}}_0,{\widehat{\mathit{C}}}_0\right)$ is different from the nominal system model $\left({\mathit{A}}_0,{\mathit{B}}_0,{\mathit{C}}_0\right)$. Even the block diagonal strictures of ${\mathit{A}}_0$ and ${\mathit{B}}_0$ are not preserved.

It can be deduced from Figure 4b on the right, Eqs. (53) and (54), that compared to the conventional scheme, the proposed emulator‐based identification is significantly more robust to variations in the operating points, which are simulated by emulator parameter perturbations.

The poles of the MISO transfer functions ${\mathit{G}}_2\left(z\right)$ of $y_2(k)$ and ${\mathit{G}}_1\left(z\right)$ of $y_1(k)$ were, respectively, $0.8500\pm \text{ j}0.3122$ and $0.7500\pm \text{ j}0.3708$. The same emulator was used for inducing phase shift to the MISO models. ${\mathit{G}}_2\left(z\right)$ with poles close to the unit circle was affected more than ${\mathit{G}}_1\left(z\right)$ with poles well inside. In view of the difference in the perturbations induced in the two models, the mean‐squared errors ${\text{mse}}_2$ and ${\text{mse}}_2^c$ are higher than ${\text{mse}}_{\text{1}}$ and ${\text{mse}}_1^c$.

7.2. Fault diagnosis

Detection of a fault: Various types of faults include (a) actuator, (b) sensor, and (c) plant, we introduced by varying the columns of ${\mathit{B}}_0$, the rows of ${\mathit{C}}_0$_, and the diagonal matrices of ${\mathit{A}}_0$. A fault is detected using appropriate test statistics depending upon the reference input waveform from Eq. (45). Since the reference input $\mathit{r}(k)$ is a constant waveform, the test statistics for the MISO and the SISO system using Eqs. (46) and (47) are

$t_s\left(e_j\right)=\left|\frac{1}{N}{\displaystyle \sum _{i=k-N+1}^k{e}_j(i)}\right|;t_s\left(e_{ji}\right)=\left|\frac{1}{N}{\displaystyle \sum _{i=k-N+1}^k{e}_{ji}(i)}\right|$E55

A visual picture of the faulty and the normal subsystems may be deduced from the autocorrelations of the residuals associated with the fault‐free, sensor fault, actuator fault, and the plant faults shown in Figure 5. Subfigures A, and B, subfigures C and D, subfigures E, and F, and subfigures G and H show respectively autocorrelations of the residual for the ideal no fault, the sensor fault, the actuator fault, and the plant fault.

Figure 5.

Autocorrelations of the residuals: ideal, sensor fault, actuator fault, and plant faults.

Remarks: The maximum value of the autocorrelation of the residual (i.e., its variance) provides an indication of the presence or an absence of the fault. In the case of the sensor fault introduced by perturbing ${\mathit{C}}_{20}$, it affects only the residual $e_2\left(k\right)$. The variance of the autocorrelation $e_2\left(k\right)$ is large while that of $e_1\left(k\right)$ indicating a fault in ${\mathit{C}}_2$. However, a fault in either the actuator or the plant, depending upon which elements of ${\mathit{B}}_0$ or ${\mathit{A}}_0$ are perturbed, may affect both residuals, and hence would be difficult to isolate.

7.2.1. Fault isolation

If a fault is asserted, and the path where the fault is located, then it is isolated using Bayesian multiple hypotheses testing scheme. The size of the fault is also estimated. The objective of fault isolation is to determine which of the emulator parameter has varied using the residual data generated or parameters using the expression for the Kalman filter residual for the model‐mismatch case given in Eq. (43). The residual $e_{ji}\left(k\right)$ is affine in the unknown emulator parameter variations $\left\{\text{Δ}{\gamma }_{ijk}\right\}$. The emulator parameter variation that is most likely to fit the perturbed residual with additive term ${\mathit{\psi }}_{jif}^T\left(z\right)\text{Δ}{\mathit{\theta }}_{ji}$ is determined sequentially by first hypothesizing single faults. If the estimates thus obtained do not fit the residual, then two simultaneous faults are hypothesized. If again the estimates do not fit the residual model, then hypothesize triple faults, and so on until the estimates fit the residual model. The maximum likelihood method, which is efficient and unbiased, is employed herein to estimate the variation $\text{Δ}\mathit{\gamma }$. The maximum likelihood estimates of the emulator parameters are obtained by minimizing the log likelihood function [13].

Let ${\mathit{H}}^{(1)}$, ${\mathit{H}}^{(2)}$_, and ${\mathit{H}}^{(3)}$ denote a hypothesis that emulator parameter ${\gamma }_{ji1}$, ${\gamma }_{ji2}$, and ${\gamma }_{ji12}$ has varied. The Kalman filter residual for ${\mathit{H}}^{(1)}$ becomes

${\mathit{H}}^{\left(1\right)}:e_{ji}^{(1)}(k)={\mathit{\psi }}_{jif}^T(k)\text{Δ}{\mathit{\theta }}_{ji}^{(1)}+{\upsilon }_{jif}(k)$E56

The least‐squares estimate $\text{Δ}{\widehat{\gamma }}_{ji1}$ from

$\text{Δ}{\widehat{\gamma }}_{ji1}=\arg \underset{\left\{\text{Δ}{\gamma }_{ji1}\right\}}{\min }\left\{{\Vert e_{ji}(k)-{\mathit{\psi }}_{jif}^T(k)\text{Δ}{\mathit{\theta }}_{ji} \Vert }^2\right\}$E57

If the estimate does not meet the criteria, then hypothesize that ${\gamma }_{ji2}$ has varied. The criteria for fitting a hypothesis are given later. The Kalman filter residual for ${\mathit{H}}^{(2)}$ becomes

${\mathit{H}}^{\left(2\right)}:e_{ji}^{(2)}(k)={\mathit{\psi }}_{jif}^T(k)\text{Δ}{\mathit{\theta }}_{ji}^{(2)}+{\upsilon }_{jif}(k)$E58

The least‐squares estimate $\text{Δ}{\widehat{\gamma }}_{ji2}$ from

$\text{Δ}{\widehat{\gamma }}_{ji2}=\arg \underset{\left\{\text{Δ}{\gamma }_{ji2}\right\}}{\min }\left\{{\Vert {\mathit{\psi }}_{jif}^T(k)\text{Δ}{\mathit{\theta }}_{ji}^{(2)}+{\upsilon }_{jif}(k)\Vert }^2\right\}$E59

If it does not meet the criteria, then hypothesize that ${\gamma }_{ji12}$ has varied. The Kalman filter residual for ${\mathit{H}}^{(3)}$ becomes

${\mathit{H}}^{(3)}:e_{ji}^{(3)}(k)={\mathit{\psi }}_{ijf}^T(k)\text{Δ}{\mathit{\theta }}_{ji}^{(3)}+{\upsilon }_{jif}(k)$E60

The least‐squares estimate $\text{Δ}{\widehat{\gamma }}_{ji12}$ from

$\text{Δ}{\widehat{\gamma }}_{ji12}=\arg \underset{\left\{\text{Δ}{\gamma }_{ji12}\right\}}{\min }\left\{{\Vert e_{ji}^{(3)}(k)-{\mathit{\psi }}_{ijf}^T(k)\text{Δ}{\mathit{\theta }}_{ji}^{(3)}\Vert }^2\right\}$E61

where $\text{Δ}{\mathit{\theta }}_{ji}^{(1)}$, $\text{Δ}{\mathit{\theta }}_{ji}^2$_, and $\text{Δ}{\mathit{\theta }}_{ji}^3$ are deviations in the feature vector when ${\gamma }_{ji1}$, ${\gamma }_{ji2}$, and ${\gamma }_{ji12}$ are assumed to have varied.

7.2.2. Criteria for asserting the hypothesis

The most likely hypotheses is determined by verifying which of the emulator parameter or parameters have varied by comparing the deviation with some threshold value

$Assert{\mathit{H}}^{\left(1\right)}if\text{Δ}{\widehat{\gamma }}_{ji1}\ge {\eta }_1$E62

$\text{Assert}{\mathit{H}}^{\left(2\right)}if\text{Δ}{\widehat{\gamma }}_{ji2}\ge {\eta }_2$E63

$\text{Assert}{\mathit{H}}^{\left(3\right)}if\text{Δ}{\widehat{\gamma }}_{ji12}\ge {\eta }_3$E64

where ${\eta }_1$, ${\eta }_2$_, and ${\eta }_3$ are threshold values. The subsystem associated with the subsystem is asserted to be faulty if the criterion is met.

8. Evaluation on physical process control system

A laboratory‐scale two‐tank physical system is formed of a controller, a DC motor, a pump, two tanks connected by a pipe, a flow rate sensor, and a liquid level sensor. The system is interfaced to a PC with the National Instruments LABVIEW for data acquisition and implementing the controller and the soft sensor [14]. The actuator, namely the pump driven by the DC motor, sends the fluid to the first tank to maintain a specified fluid level in the second tank. An evaluation of the proposed scheme for fault diagnosis was performed on a benchmark laboratory‐scale process control system using the National Instruments LABVIEW as shown below in Figure 6. The sampling period is $T_s=0.05$.

Figure 6.

Process control system: controller, actuator, and tank.

Emulator‐generated height and flow rate profiles under various types of faults are shown in under the caption Height/Flow rate Profiles for PI controller with Consumer in Fig. 7. Figures 7a–c show the height and flow rate profiles when subjected to (a) leakage fault, (b) actuator fault, and (c) sensor faults, respectively. The height profile is shown on the top and the flow rate profile is shown at the bottom of Figure 7. The faults are induced by varying the appropriate emulator parameters to 0.25, 0.5, and 0.75 time the nominal values, in order to represent “small,” “medium,” and “large” faults. However, by virtue of its control design objective, the closed‐loop PI controller will hide any fault that may occur in the system and hence will make it difficult to detect it. In addition, the physical system exhibits a highly nonlinear behavior. The flow rate saturates at 4.5 ml/s. The dead‐band effect in the actuator exhibits itself as a delay in the output response: when a step reference input is applied, the height output responds after some delay, as a minimum force is required to drive the actuator. These nonlinearities affect the steady‐state value of the height: even though there is an integral action in the closed‐loop control system, the steady‐state error is non‐zero for a constant reference input.

Figure 7.

Emulator‐generated data: height and flow.

The system is modeled as a single‐input, multi‐output system where $r$ is the reference input, and the outputs are the control input $u$, the flow rate $f$ and the height $h$. Faults were induced in the height sensor, the flow sensor, the actuator, and also as a leakage. The proposed fault diagnosis successfully detected and isolated all the faults compared to SISO scheme [14], where all the faults were detected and isolated using the reference input and the height output.

9. Conclusions

Fault detection and isolation of a class of linear multiple‐input and multiple‐output system based on the Kalman residual and the emulators were presented. The key properties of the Kalman filter, namely the residual, is a zero‐mean white noise process if and only if there is no model mismatch, drive the prediction error identification of the nominal system model, and the Kalman filter. In view of the closed‐loop configuration, the noise and the disturbance are attenuated at the estimated output. The Kalman filter is the best minimum variance estimator in the class of all linear estimators.

To handle fault isolation, the powerful and effective concept of emulators was introduced. Similar in spirit to the training of the artificial neural network, a number of emulator parameter‐perturbed experiments were performed to capture the perturbation model of the subsystems to help with fault isolation. The influence vectors of the emulator parameters, which are indirectly the associated subsystems, were estimated. The influence vectors captured the emulator perturbation model and hence that of the subsystem.

The residual of the Kalman filter was shown to have an additive fault indicating term when there is a model mismatch due to emulator perturbations. The model‐mismatch term is affine in the emulator parameter variations. Using the expression for the fault indicating term, the fault was isolated using the influence vectors and its size was estimated. The residual, being affine in the emulator parameter variation, easily lends itself to the widely used and successful composite Bayes hypothesis‐testing scheme for fault isolation.

The future work generated from this work includes its extension to a class of nonlinear multiple‐input and multiple‐output systems, and the development of a computationally efficient identification of the Kalman filter directly from the input data even for unstable systems. Although a gold standard for system identification, the prediction error method involves a nonlinear optimization problem and hence can suffer from the existence of local minima. Unlike the least‐squares approach, it does not offer a closed‐form solution to the parameter estimation problem. Instead, it relies on a recursive solution that may be time‐consuming (slow convergence rate), computationally complex, and which may also suffer from initialization problems.

Acknowledgments

The first author acknowledges support of the Department of Electrical and Computer Engineering, the University of New Brunswick, and the National Science and Engineering Research Council (NSERC) of Canada. The author is grateful to Professor C.P. Diduch of the University of New Brunswick, Mr. Jiong Tang of MDS and Mr. H.M. Khalid of KFUPM for their help and suggestions. The second author acknowledges the support of KFUPM, Saudi Arabia, and the help of Mr. H. M. Khalid of KFUPM.