Novel Direct and Accurate Identification of Kalman Filter for General Systems Described by a Box-Jenkins Model

Abstract

A novel robust Kalman filter (KF)-based controller is proposed for a multivariable system to accurately track a specified trajectory under unknown stochastic disturbance and measurement noise. The output is a sum of uncorrelated signal, disturbance and measurement noise. The system model is observable but not controllable while the signal one is controllable and observable. An emulator-based two-stage identification is employed to obtain a robust model needed to design the robust controller. The system and KF are identified and the signal and output error estimated. From the identified models, minimal realizations of the signal and KF, the disturbance model and whitening filter are obtained using balanced model reduction techniques. It is shown that the signal model is a transfer matrix relating the system output and the KF residual, and the residual is the whitened output error. The disturbance model is identified by inverse filtering. A feedback-feedforward controller is designed and implemented using an internal model of the reference driven by the error between the reference and the signal estimate, the feedforward of reference and output error. The successful evaluation of the proposed scheme on a simulated autonomously-guided drone gives ample encouragement to test it later, on a real one.

1. Introduction

In conventional Kalman filter applications, the system involved is typically linearized and then identified. Based on the identified system, the Kalman filter is then identified. In this chapter, we propose a novel approach in that (a) the system is represented by a more general model, termed multi-input multi-output Box-Jenkins (MIMO BJ model) which subsumes all previous classical models, such as ARMA models and their derivatives, and (b) the associated Kalman filter identification is carried out directly, i.e. it does not necessitate the prior identification of the system involved. The various tools involved in our proposed approach are all explained below.

1.1 Box-Jenkins model and its applications

Identification of a class of system described by MIMO BJ model, and the associated Kalman filter directly from the input-output data is proposed [1, 2]. There is no need to specify the covariance of the disturbance and the measurement noise, thereby avoiding the use the Riccati equation to solve for the Kalman gain. The output is the desired waveform, termed signal, corrupted by a stochastic disturbance and zero-mean white measurement noise. The state-space BJ model is an augmented system formed of the signal and disturbance model. The signal model and the disturbance models are driven respectively by a user-defined accessible input, and an inaccessible zero-mean white noise process. The signal model is generally a cascade, parallel and feedback combinations of subsystems such as controllers, actuators, plants, and sensors [3]. Unlike the ARMA model, the Box-Jenkins model is observable but not controllable while the signal model is both controllable and observable. In other words, the transfer matrix of the system is non-minimal whereas that of the signal is minimal. This issue will need to be addressed in the identification and implementation of the Kalman filter.

1.2 Kalman filter and its key properties

The structure of the Kalman filter is determined using the internal model principle which establishes the necessary and sufficient condition for the tracking of the output of a dynamical system [3, 4]. In accordance with this principle, the Kalman filter consists of (a) a copy of the system model driven by the residuals, and (b) a gain term, termed the Kalman gain, to stabilize the filter. The Kalman gain is determined such that the residual of the Kalman filter is a zero-mean white noise process with minimum variance. The Kalman filter enjoys the following key properties:

Tracking a signal: The estimate of the Kalman filter tracks a given signal if and only if the model that generates the signals including those of the noise and disturbances is embodied in the Kalman filter. In other words, the Kalman filter tracks the input, thanks to its internal model-based structure [3, 4].

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 version embodied in the Kalman filter, and its variance is minimum [4].

Optimality: 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 [5].

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 [5].

Model-mismatch: If there is a model mismatch, then the residual will not be a zero-mean white noise process and an additive term termed fault-indicative term will occur. The fault-indicative term is a filtered version of the deviation in the linear regression model of the system or that of the signal [6, 7, 8].

1.3 Identification using residual model

The equation error in the regression model of the system is a colored noise process, and hence a direct identification of the system model from the input-output data by minimizing the equation error will not ensure that the estimates are consistent, unbiased, and efficient. The fundamental requirement of identification is that the leftover signal from identification, namely the residual is a zero-mean white noise process that contains no information. To meet this requirement, both the input and output of the system are filtered. Among the class of all linear whitening filters, the Kalman filter is the best. The system model is indirectly identified by minimizing the residual generated by the Kalman filter instead of the equation error. The Subspace Method (SM) uses the structure of the state-space model of the Kalman filter, whereas the prediction error method (PEM), which is the gold standard in system identification, is developed from the residual model [1, 9, 10, 11].

1.4 Emulator-based two-stage identification

The static and dynamic behavior of a physical system change as a result of variations in the parameters of some of its subsystems such as sensors, actuators, plant, disturbance models, and controllers. As the parameters of these subsystems are not generally accessible to generate data, instead, emulators, which are hardware or software devices, are connected in cascade to the output, input or both, of the subsystems. An emulator is a transfer function block which mimics the variations in the associated subsystems including the disturbance model. An emulator takes the form of a static gain or an all-pass filter to induce gain or phase variations in the subsystem it is connected to. Emulator parameters are perturbed to mimic various normal and abnormal, or faulty, operating scenarios resulting from variations in these subsystems. The emulator-generated data is employed in (a) the identification of robust systems and signal models and their associated Kalman filters using the two-stage identification scheme [2, 3, 6, 7, 8].

A two-stage identification is used in various applications including the non-parametric identification of impulse response, estimation of Markov parameters in the SM, in model predictive control, identification of a signal model and in system identification. The use of the two-stage identification is inspired by the seminal paper by [12] for an accurate estimation of the parameters of an impulse response from measurements in an additive white noise. It is shown via simulation that the variance of the parameter estimation error approaches the Cramer-Rao lower bound [13]. Further, it is shown analytically that using a high-order model (with an order several times larger than the true order) improves significantly the accuracy of the parameter estimates. The two-stage scheme has not received much attention in system identification although it has been mentioned as an alternative scheme to the PEM [1, 14], and has been successfully employed in identification in [15, 16, 17].

It should be emphasized that the prediction error method (PEM), viewed as a gold standard for system identification, is not geared for the estimation of the signal buried in the output, i.e. it is developed for the ARMA model and not for the Box-Jenkins one. A two-stage identification of the Box-Jenkins model is proposed as the system model is observable but not controllable while the signal model is both controllable and observable:

• In the first stage, the robust system model and the associated Kalman filter are identified using the emulator-generated data using PEM, and the signal and the output error are both estimated. Further, the whitening filter that relates the output error and the residual is obtained.

• In the second stage, minimal realizations of the signal model and the associated Kalman filter are obtained using model reduction method [18].

The high- order for the first stage and the reduced-order for the second are both selected using the Akaike information criterion (AIC) and are cross-checked by verifying the whiteness of the associated residual. The two-stage identification has also been successfully employed in the identification model.

The question arises as to how to obtain the system model and the signal model from the identified high-order model Kalman filters. A key property of the Kalman filter is established here, namely that the transfer matrix of the signal and the system is the matrix fraction description model derived from the Kalman filter residual model of the system. This property is exploited to derive the signal and the system transfer matrices. The state-space models of the signal and the system models are derived from the identified state-space models of the Kalman filters. Thus, the proposed scheme identifies (a) the Kalman filter for the system, (b) the Kalman filter for the signal first. Then, the system model and the signal model are separately obtained.

The proposed scheme is further extended to identify the signal model to complement the PEM. In the first stage, a very high-order model is identified using PEM. In the second stage, the signal model is identified using a balanced model reduction of the high-order identified model obtained in the first stage. The PEM and state-space method (SM) are both tailored to identify the signal model and estimate the signal by employing the proposed version of the two-stage identification scheme. The results of the comparison of the performance of these methods in identifying the system and signal models are presented.

1.5 Highlights of the contributions

• The Auto-Regressive (AR), The Moving Average (MA), and the Auto-Regressive and Moving Average (ARMA) models are all special cases of the proposed Box-Jenkins model. As this model is more general and hence has wider applications, including robust controller design; estimation of latent variables; monitoring of the status of the system, fault diagnosis, development of condition-based maintenance programs and design of fault-tolerant systems; filtering of signals, speech enhancement, noise and echo cancelation in communication; 2-D image filtering and tracking of moving objects.

• The state-space models of the system and the signal models are derived from the identified Kalman filters, by invoking the (causal) invertibility of the output error and the residual [5].

• An efficient scheme to monitor the status of the system may be implemented from the proposed scheme. First the status of the system is monitored by analyzing the residual of the Kalman filter of the system model. If there is a variation, then the residual of the Kalman filter of the signal model is analyzed to ascertain whether a fault has occurred.

• In practice, disturbances are inevitable, and can negatively affect the system performance. When the system is in an abnormal state, it is not in general easy to determine whether the abnormal operation is the result of variations in the disturbance or the occurrence of a fault. The proposed scheme provides a simple solution by analyzing the residuals of both Kalman filters as the residual of the Kalman filter for the system captures the variations in both the system and disturbance models, while that of the Kalman filter for the signal, and captures only the variations in the signal model. This is crucial for reducing false alarm, and all its concomitant risks and costs, resulting from variations in the disturbance and not in the signal model [19].

• The PEM (SM) may be tailored to identify the signal model and estimate the signal itself by using the proposed two-stage identification scheme.

1.6 Applications

Applications include monitoring the status of the system and the signal models, distinguishing between the variations in the disturbance model and those in the signal model to help diagnose a fault in the system and ensure a low false alarm probability, estimating the latent variable, namely the signal, developing a framework for applications including robust controller design; fault diagnosis; speech and biological signal processing; tracking of moving objects, design of soft sensors to replace maintenance-prone hardware sensors, evaluate and monitor product quality, meeting the ever-increasing need for fault-tolerant systems for mission-critical systems found in aerospace, the nuclear power systems, and autonomous vehicles.

2. Problem formulation

The output $\mathit{y}\left(k\right)\in {\mathbb{R}}^q$ is an additive sum of the signal $\mathit{s}\left(k\right)\in {\mathbb{R}}^q$ , disturbance, $\mathit{d}\left(k\right)\in {\mathbb{R}}^q$ and the measurement noise $\mathit{v}\left(k\right)\in {\mathbb{R}}^q$ where $\mathbb{R}$ is real scalar field.

$$\mathit{y}\left(k\right)=\mathit{s}\left(k\right)+\mathit{d}\left(k\right)+\mathit{v}\left(k\right)$$ E1

Where the signal and the disturbance models are:

$$\mathit{s}\left(z\right)={\mathit{G}}_s\left(z\right)\mathit{u}\left(z\right)$$ E2

$$\mathit{d}\left(z\right)={\mathit{G}}_w\left(z\right)\mathit{w}\left(z\right)$$ E3

Where $\mathit{u}\left(k\right)\in {\mathbb{R}}^q$ is the input; $\mathit{w}\left(k\right)\in {\mathbb{R}}^p$ is zero-mean white noise process that generates the disturbance $\mathit{d}\left(k\right)\in {\mathbb{R}}^p$ , and is uncorrelated with the measurement noise $\mathit{v}\left(k\right)$ ; ${\mathit{G}}_s\left(z\right)=D_s^{-1}\left(z\right){\mathit{N}}_s\left(z\right)$ and ${\mathit{G}}_w\left(z\right)=D_w^{-1}\left(z\right){\mathit{N}}_w\left(z\right)$ are $\mathit{qxp}$ transfer matrix of order $n_s$ and $n_w$ respectively; $\mathit{\vartheta }\left(k\right)\mathit{=}\mathit{d}\left(k\right)+\mathit{v}\left(k\right)$ is the output error.

The signal model ${\mathit{G}}_s\left(z\right)$ is formed of cascade and parallel combinations of the subsystems such as actuators, plant and the sensors. Let the state space model of the signal and the disturbance models be respectively $\left({\mathit{A}}_s,{\mathit{B}}_s,{\mathit{C}}_s\right)$ and $\left({\mathit{A}}_w,{\mathit{B}}_w,{\mathit{C}}_w\right)$ .

Figure 1 shows the input-output model relating the input, the signal model, the signal, the disturbance model, the disturbance, the measurement noise and the output.

Figure 1.

System: signal, the disturbance and the measurement noise.

Linear regression model:

Using, (1)–(3), the expression for the linear regression becomes:

$$\begin{array}{l}{D}_{\mathit{sw}}\left(z\right)\mathit{y}\left(z\right)=D_w\left(z\right){\mathit{N}}_{\mathit{sw}}\left(z\right)\mathit{u}\left(z\right)+\mathit{\upsilon }\left(z\right)\\ \mathit{\upsilon }\left(z\right)=D_s\left(z\right){\mathit{N}}_w\left(z\right)\mathit{w}\left(z\right)+D_{\mathit{sw}}\left(z\right)\mathit{v}\left(z\right)\end{array}$$ E4

Where $\mathit{\upsilon }\left(z\right)=D\left(z\right)\mathit{\vartheta }\left(z\right)$ is the equation error; $D_{\mathit{sw}}\left(z\right)=D_s\left(z\right)D_w\left(z\right)$ and ${\mathit{N}}_{\mathit{sw}}\left(z\right)=D_w\left(z\right){\mathit{N}}_s\left(z\right)$ are respectively the denominator and numerator polynomials. The model is termed Box-Jenkins model.

Note that the model that generates the equation error $\mathit{\upsilon }\left(z\right)$ is a Moving Average (MA) model, whereas the one that generates the output error $\mathit{\vartheta }\left(k\right)$ is an Auto-Regressive Moving-Average (ARMA) model.

Augmented state-space model: The augmented state-space representation of the multi-input and multi-output (MIMO) system $\left(\mathit{A},\mathit{B},\mathit{C},\mathit{D}\right)$ formed of the signal model $\left({\mathit{A}}_s,{\mathit{B}}_s,{\mathit{C}}_s,{\mathit{D}}_s\right)$ and $\left({\mathit{A}}_w,{\mathit{B}}_w,{\mathit{C}}_w,{\mathit{D}}_w\right)$ representing a p-input, q-output system, is given by:

$$\begin{array}{l}\mathit{x}\left(k+1\right)=\mathit{Ax}\left(k\right)+\mathit{Bu}\left(k\right)+{\mathit{E}}_w\mathit{w}\left(k\right)\\ \mathit{s}\left(k\right)={\mathit{C}}_s\mathit{x}\left(k\right)+{\mathit{D}}_s\mathit{u}\left(k\right)\\ \mathit{y}\left(k\right)=\mathit{Cx}\left(k\right)+\mathit{Du}\left(k\right)+\mathit{v}\left(k\right)\end{array}$$ E5

Where $\mathit{A}=\left[\begin{array}{ll}{\mathit{A}}_s& \mathbf{0}\\ \mathbf{0}& {\mathit{A}}_w\end{array}\right];\mathit{B}=\left[\begin{array}{l}{\mathit{B}}_s\\ \mathbf{0}\end{array}\right];{\mathit{E}}_w=\left[\begin{array}{l}\mathbf{0}\\ {\mathit{B}}_w\end{array}\right];\mathit{C}=\left[{\mathit{C}}_s{\mathit{C}}_w\right];$ $\mathit{A}\in {\mathbb{R}}^{\mathit{nxn}}$ is an augmented state transition matrix formed of ${\mathit{A}}_s\in {\Re }^{n_s{\mathit{xn}}_s}$ and ${\mathit{A}}_w\in {\Re }^{n_w{\mathit{xn}}_w}$ ; $\mathit{B}\mathit{=}\left[\begin{array}{l}{\mathit{B}}_1 {\mathit{B}}_2 . {\mathit{B}}_p\end{array}\right]\in {\mathbb{R}}^{\mathit{nxp}}$ ; $\mathit{C}\mathit{=}{\left[{\begin{array}{l}{\mathit{C}}_1 {\mathit{C}}_2 . \mathit{C}\end{array}}_q\right]}^T\in {\mathbb{R}}^{\mathit{qxn}}$ ; ${\mathit{E}}_w\in {\mathbb{R}}^{\mathit{nxp}}$ is a disturbance entry matrix;

$\mathit{x}\left(k\right)={\left[\begin{array}{lllll}{x}_1\left(k\right)& x_2\left(k\right)& x_3\left(k\right)& \dots & x_n\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^n$ ; $\mathit{s}\left(k\right)={\left[\begin{array}{lllll}{s}_1\left(k\right)& s_2\left(k\right)& s_3\left(k\right)& \dots & s_q\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^q$ ;

$\mathit{u}\left(k\right)={\left[\begin{array}{l}{u}_1\left(k\right) u_2\left(k\right) u_3\left(k\right) \dots u_p\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^p$ ; $\mathit{y}\left(k\right)={\left[\begin{array}{l}{y}_1\left(k\right) y_2\left(k\right) y_3\left(k\right) \dots y_q\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^q$ are respectively the state, the input and output; $n=n_s+n_w$ is the order; $\mathit{w}\left(k\right)\in {\mathbb{R}}^n$ and $\mathit{v}\left(k\right)\in {\mathbb{R}}^q$ are respectively the disturbances and measurement noise; $D_{\mathit{sw}}\left(z\right)=\left|\left(z\mathit{I}-{\mathit{A}}_s\right)\right|\left|\left(z\mathit{I}-{\mathit{A}}_w\right)\right|$ where $\left|\left(.\right)\right|$ is the determinant of $\left(.\right)$ . Using $D_{\mathit{sw}}\left(z\right)=D_s\left(z\right)D_w\left(z\right)$ and ${\mathit{N}}_{\mathit{sw}}\left(z\right)=D_w\left(z\right){\mathit{N}}_s\left(z\right)$ we get:

$$\begin{array}{l}\mathit{G}\left(z\right)=\mathit{C}{\left(z\mathit{I}-\mathit{A}\right)}^{-1}\mathit{B}\mathit{+}\mathit{D}\mathit{=}{D}_{\mathit{sw}}^{-1}\left(z\right){\mathit{N}}_{\mathit{sw}}\left(z\right)\\ \mathit{G}\left(z\right)=D_s^{-1}\left(z\right){\mathit{N}}_s\left(z\right)={\mathit{G}}_s\left(z\right)\end{array}$$ E6

The augmented transfer matrix is not a minimal realization of the system output model as there is (stable) pole-zero cancelation since the polynomial $D_w\left(z\right)$ , which is common to both the numerator ${\mathit{N}}_{\mathit{sw}}\left(z\right)$ and the denominator $D_{\mathit{sw}}\left(z\right)$ . In other words, ${\mathit{N}}_{\mathit{sw}}\left(z\right)$ and $D_{\mathit{sw}}\left(z\right)$ are not coprime. The signal model $\left({\mathit{A}}_s,{\mathit{B}}_s,{\mathit{C}}_s,{\mathit{D}}_s\right)$ associated with signal model ${\mathit{G}}_s\left(z\right)$ is controllable and observable while $\left(\mathit{A},\mathit{B},\mathit{C},\mathit{D}\right)$ , associated with $\mathit{G}\left(z\right)$ , is merely an observable. $\left({\mathit{A}}_s,{\mathit{B}}_s,{\mathit{C}}_s,{\mathit{D}}_s\right)$ ( ${\mathit{G}}_s\left(z\right)$ ) is a minimal realization of $\left(\mathit{A},\mathit{B},\mathit{C},\mathit{D}\right)$ ( $\mathit{G}\left(z\right)$ ).

Assumptions: It is assumed that (a) the disturbance $\mathit{w}\left(k\right)$ and the measurement noise $\mathit{v}\left(k\right)$ are independent zero-mean Gaussian white noise processes with unknown but finite covariance, $\mathit{Q}=E\left[\mathit{w}\left(k\right){\mathit{w}}^T\left(k\right)\right]$ and $\mathit{R}=E\left[\mathit{v}\left(k\right){\mathit{v}}^T\left(k\right)\right]$ , respectively, and are inaccessible, (b) $\left(\mathit{A},\mathit{C}\right)$ is observable, (c) the signal and disturbance models are both minimal, $\left({\mathit{A}}_s,{\mathit{B}}_s,{\mathit{C}}_s,{\mathit{D}}_s\right)$ and $\left({\mathit{A}}_w,{\mathit{B}}_w,{\mathit{C}}_w,{\mathit{D}}_w\right)$ are both controllable and observable, (d) The initial conditions $\mathit{x}\left(0\right)$ , $\mathit{w}\left(k\right)$ and $\mathit{v}\left(k\right)$ are mutually uncorrelated. However, the signal $\mathit{s}\left(z\right)$ and the disturbance $\mathit{w}\left(z\right)$ may have spectral overlap, (e) the output error is bounded.

2.1 Kalman filter

Predictor form: A robust Kalman filter of the identified system $\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)$ relating the system input $\mathit{u}\left(z\right)$ and system output $\mathit{y}\left(k\right)$ to the estimated output $\widehat{\mathit{y}}\left(k\right)$ is:

$$\begin{array}{l}\widehat{\mathit{x}}\left(k+1\right)=\left({\mathit{A}}^0-{\mathit{K}}^0{\mathit{C}}^0\right)\widehat{\mathit{x}}\left(k\right)+\left({\mathit{B}}^0-{\mathit{K}}^0{\mathit{D}}^0\right)\mathit{u}\left(k\right)+{\mathit{K}}^0\mathit{y}\left(k\right)\\ \widehat{\mathit{y}}\left(k\right)={\mathit{C}}^0\widehat{\mathit{x}}\left(k\right)+{\mathit{D}}^0\mathit{u}\left(k\right)\\ \mathit{e}\left(k\right)=\mathit{y}\left(k\right)-\widehat{\mathit{y}}\left(k\right)\end{array}$$ E7

Where $\widehat{\mathit{x}}\left(k\right)={\left[\begin{array}{lllll}{\widehat{x}}_1\left(k\right)& {\widehat{x}}_2\left(k\right)& {\widehat{x}}_3\left(k\right)& \dots & {\widehat{x}}_n\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^n$ and $\widehat{\mathit{y}}\left(k\right)={\left[\begin{array}{lllll}{\widehat{y}}_1\left(k\right)& {\widehat{y}}_2\left(k\right)& {\widehat{y}}_3\left(k\right)& \dots & {\widehat{y}}_q\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^q$ are respectively the best estimate of the state $\mathit{x}\left(k\right)$ , and of the output $\mathit{y}\left(k\right)$ ; $\mathit{e}\left(k\right)={\left[\begin{array}{lllll}{e}_1\left(k\right)& e_2\left(k\right)& e_3\left(k\right)& \dots & e_q\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^q$ is the residual or the innovation sequence; the Kalman gain ${\mathit{K}}^0\in {\mathbb{R}}^{\mathit{nxq}}$ ensures the asymptotic stability of the Kalman filter, i.e. ( ${\mathit{A}}^0-{\mathit{K}}^0{\mathit{C}}^0$ ) is strictly Hurwitz having all its eigenvalues strictly inside the unit circle.

Innovation form: There is duality between the predictor, and the innovation forms of the Kalman filter [5]. The output $\mathit{y}\left(k\right)$ and the residual $\mathit{e}\left(k\right)$ are (causally) invertible. In other words, $\mathit{e}\left(k\right)$ can be generated from the output $\mathit{y}\left(k\right)$ and $\mathit{r}\left(k\right)$ using the (causal) predictor form, and $\mathit{y}\left(k\right)$ can be generated from $\mathit{e}\left(k\right)$ and $\mathit{r}\left(k\right)$ using the innovation form. The Kalman filter given by (7) is termed the predictor form and can be expressed in an alternative form, termed the innovation form, given by:

$$\begin{array}{l}\widehat{\mathit{x}}\left(k+1\right)={\mathit{A}}^0\widehat{\mathit{x}}\left(k\right)+{\mathit{B}}^0\mathit{u}\left(k\right)+{\mathit{K}}^0\mathit{e}\left(k\right)\\ \widehat{\mathit{y}}\left(k\right)={\mathit{C}}^0\widehat{\mathit{x}}\left(k\right)+{\mathit{D}}^0\mathit{u}\left(k\right)\end{array}$$ E8

Figure 2 shows the system and the Kalman filter which embodies the system model $\left(\mathit{A},\mathit{B},\mathit{C}\right)$ . The inputs to the Kalman filter are the input $\mathit{r}\left(k\right)$ and the output $\mathit{y}\left(k\right)$ which is corrupted by the noise $\mathit{v}\left(k\right)$ and affected by the disturbance $\mathit{w}\left(k\right)$ .

Figure 2.

The system and the Kalman filter model.

2.2 Residual model

The frequency-domain expression relating the input $\mathit{u}\left(z\right)\in {\mathbb{R}}^p$ and the output $\mathit{y}\left(z\right)\in {\mathbb{R}}^q$ to the residual $\mathit{e}\left(z\right)\in {\mathbb{R}}^q$ is given by the following model termed the residual model:

$$\mathit{e}\left(z\right)=F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)\mathit{y}\left(z\right)-F^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)\mathit{u}\left(z\right)$$ E9

where $\overline{\mathit{D}}\left(z\right)$ and $\overline{\mathit{N}}\left(z\right)$ are matrix polynomials, $F\left(z\right)$ is the scalar characteristic polynomial termed Kalman polynomial, $F\left(z\right)=\left|z\mathit{I}-{\mathit{A}}^0\mathit{+}{\mathit{K}}^0{\mathit{C}}^0\right|$ ; $D\left(z\right)=\left|\left(z\mathit{I}-{\mathit{A}}^0\right)\right|$ ; $\overline{\mathit{D}}\left(z\right)=F\left(z\right)\left(\mathit{I}-{\mathit{C}}^0{\left(z\mathit{I}-{\mathit{A}}^0\mathit{+}{\mathit{K}}^0{\mathit{C}}^0\right)}^{-1}{\mathit{K}}^0\right)$ is $\mathit{qxq}$ matrix; $\overline{\mathit{N}}\left(z\right)=F\left(z\right)\left({\mathit{C}}^0{\left(z\mathit{I}-{\mathit{A}}^0\mathit{+}{\mathit{K}}^0{\mathit{C}}^0\right)}^{-1}\left({\mathit{B}}^0-{\mathit{K}}^0{\mathit{D}}^0\right)+{\mathit{D}}^0\right)$ is $\mathit{qxp}$ matrix; $\mathit{I}\in {\mathbb{R}}^{\mathit{qxq}}$ is an identity matrix;

$$\overline{\mathit{D}}\left(z\right)=\left[\begin{array}{l}{\overline{\mathit{D}}}_1\left(z\right)\\ {\overline{\mathit{D}}}_2\left(z\right)\\ .\\ {\overline{\mathit{D}}}_q\left(z\right)\end{array}\right]=\left[\begin{array}{llll}{\overline{D}}_{11}\left(z\right)& {\overline{D}}_{12}\left(z\right)& .& {\overline{D}}_{1q}\left(z\right)\\ {\overline{D}}_{21}\left(z\right)& {\overline{D}}_{22}\left(z\right)& .& {\overline{D}}_{2q}\left(z\right)\\ .& ;& .& .\\ {\overline{D}}_{q1}\left(z\right)& {\overline{D}}_{q2}\left(z\right)& .& {\overline{D}}_{\mathit{qq}}\left(z\right)\end{array}\right];\overline{\mathit{N}}\left(z\right)=\left[\begin{array}{l}{\overline{\mathit{N}}}_1\left(z\right)\\ {\overline{\mathit{N}}}_2\left(z\right)\\ .\\ {\overline{\mathit{N}}}_q\left(z\right)\end{array}\right]=\left[\begin{array}{llll}{\overline{N}}_{11}\left(z\right)& {\overline{N}}_{12}\left(z\right)& .& {\overline{N}}_{1p}\left(z\right)\\ {\overline{N}}_{21}\left(z\right)& {\overline{N}}_{22}\left(z\right)& .& {\overline{N}}_{2p}\left(z\right)\\ .& ;& .& .\\ {\overline{N}}_{q1}\left(z\right)& {\overline{N}}_{q2}\left(z\right)& .& {\overline{N}}_{\mathit{qp}}\left(z\right)\end{array}\right]$$ E10

${\overline{D}}_{\mathit{ij}}\left(z\right)=\sum _{\mathrm{\ell }=0}^n{\overline{a}}_{\mathit{ij}\mathrm{\ell }}{z}^{-\mathrm{\ell }}$ ; ${\overline{N}}_{\mathit{ij}}\left(z\right)=\sum _{\mathrm{\ell }=1}^n{\overline{b}}_{\mathit{ij}\mathrm{\ell }}{z}^{-\mathrm{\ell }}$ ; ${\overline{a}}_{\mathit{ij}\mathrm{\ell }}$ and ${\overline{b}}_{\mathit{ij}\mathrm{\ell }}$ are the coefficients of the polynomials ${\overline{D}}_{\mathit{ij}}\left(z\right)$ and ${\overline{N}}_{\mathit{ij}}\left(z\right)$ , respectively. The rational polynomials $F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)$ and $F^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)$ associated with the system output $\mathit{y}\left(z\right)$ and the input $\mathit{u}\left(z\right)$ are termed as an output IIR filter, and an input IIR filter, respectively. The estimate of the Kalman filter $\widehat{\mathit{y}}\left(k\right)$ is:

$$\widehat{\mathit{y}}\left(z\right)=\left(\mathit{I}-F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)\right)\mathit{y}\left(z\right)+F^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)\mathit{u}\left(z\right)$$ E11

The residual model of the Kalman filter forms the backbone of the proposed identification scheme.

2.3 The key properties of the Kalman filter

The map relating the signal and its model, and the output IIR filter and an input IIR filter of residual model is developed next.

The following lemmas are developed by invoking the key property namely that 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 [4], that is, the identified model embodied in the Kalman filter is identical to that of the actual model:

2.3.1 Derivation of the signal and the signal model

The following Lemma 1 shows that (a) the estimate of the signal model is the matrix fraction description relating the transfer matrices relating the residual of the Kalman filter to the input, and the output of the system; (b) the estimate of the signal is its output generated by the system input; and the Kalman filter whitens the output error. Lemma 1

(a) The left-matrix description of the MIMO signal model derived from the state-space model $\left({\mathit{A}}_s^0,{\mathit{B}}_s^0,{\mathit{C}}_s^0,{\mathit{D}}_s^0\right)$ , namely, ${\mathit{G}}_s\left(z\right)={\mathit{C}}_s^0{\left(z\mathit{I}-{\mathit{A}}_s^0\right)}^{-1}{\mathit{B}}_s^0\mathit{=}{D}_s^{-1}\left(z\right){\mathit{N}}_s\left(\mathrm{z}\right)$ and the left-matrix description of the Kalman filter derived from the residual model, $\overline{\mathit{G}}\left(z\right)={\overline{\mathit{D}}}^{-1}\left(z\right)\overline{\mathit{N}}\left(\mathrm{z}\right)$ are identical. The signal model ${\mathit{G}}_s\left(z\right)$ and the signal $\mathit{s}\left(z\right)$ are:

$${{\displaystyle \widehat{\mathit{G}}}}_s\left(z\right)={\displaystyle \widehat{\mathit{G}}}\left(z\right)={\displaystyle \overline{\mathit{G}}}\left(z\right){\displaystyle \widehat{\mathit{s}}}\left(z\right)={{\displaystyle \widehat{\mathit{G}}}}_s\left(z\right)\mathit{u}\left(z\right)={\displaystyle \overline{\mathit{G}}}\left(z\right)\mathit{u}\left(z\right)$$ E12

Proof:

(a) Consider the residual model (9). Substituting for $\mathit{y}\left(z\right)$ yields:

$$F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)\left(D_s^{-1}\left(z\right){\mathit{N}}_s\left(z\right)\mathit{u}\left(z\right)-{\overline{\mathit{D}}}^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)\mathit{u}\left(z\right)+\mathit{\vartheta }\left(z\right)\right)=\mathit{e}\left(z\right)$$ E13

Since the residual is a zero-mean, white noise process and is uncorrelated with $\mathit{u}\left(z\right)$ and $\mathit{\upsilon }\left(z\right)$ , correlating both sides with the input $\mathit{u}\left(z^{-1}\right)$ yields:

$$\left(F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right){\widehat{D}}_s^{-1}\left(z\right){\widehat{\mathit{N}}}_s\left(z\right)-F^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)\right)E\left[\mathit{u}\left(z\right)\mathit{u}\left(z^{-1}\right)\right]=0$$ E14

Assuming that the input correlation is not identically equal to zero, i.e. $E\left[\mathit{u}\left(z\right)\mathit{u}\left(z^{-1}\right)\right]\ne \mathbf{0}$ yields:

$$\left(F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right){\widehat{D}}_s^{-1}\left(z\right){\widehat{\mathit{N}}}_s\left(z\right)-F^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)\right)=0$$ E15

Simplifying we get:

$${\widehat{D}}_s^{-1}\left(z\right){\widehat{\mathit{N}}}_s\left(z\right)={\overline{\mathit{D}}}^{-1}\left(z\right)\overline{\mathit{N}}\left(z\right)$$ E16

Hence ${\widehat{\mathit{G}}}_s\left(z\right)=\overline{\mathit{G}}\left(z\right)$ holds. Since $D_{\mathit{sw}}\left(z\right)=D_s\left(z\right)D_w\left(z\right)$ and ${\mathit{N}}_{\mathit{sw}}\left(z\right)=D_w\left(z\right){\mathit{N}}_s\left(z\right)$ are not coprime as $D_w\left(z\right)$ is a common factor, then ${\widehat{\mathit{G}}}_s\left(z\right)=\widehat{\mathit{G}}\left(z\right)=\overline{\mathit{G}}\left(z\right)$ .

(b) Substituting (16) in (13) we get

$$\mathit{e}\left(z\right)=F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)\mathit{\vartheta }\left(z\right)=F^{-1}\left(z\right)D^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)\mathit{\upsilon }\left(z\right)$$ E17

where $F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)=\left(\mathit{I}-\mathit{C}{\left(z\mathit{I}-\mathit{A}\mathit{+}\mathit{KC}\right)}^{-1}\mathit{K}\right)$ _.

Proof: follows from $\mathit{\vartheta }\left(k\right)\mathit{=}\mathit{d}\left(k\right)+\mathit{v}\left(k\right)$ , (9) and

2.3.2 Derivation of the output error and its model

The following Lemma 2 shows that the output error is the difference between the output and the estimate of the signal; the estimate of the output error model is the matrix fraction description of the transfer matrices relating the residual of the Kalman filter to the input; the estimate of the output error is obtained as its output when its input is the residual. It is assumed that the transfer matrix of the Kalman filter $F^{-1}\left(z\right)\overline{\mathit{D}}\left(z\right)=\left(\mathit{I}-{\mathit{C}}^0{\left(z\mathit{I}-{\mathit{A}}^0\mathit{+}{\mathit{K}}^0{\mathit{C}}^0\right)}^{-1}{\mathit{K}}^0\right)$ relating the output $\mathit{y}\left(z\right)$ and the residual $\mathit{e}\left(z\right)$ is minimum-phase, that is both the numerator and the denominator polynomials are asymptotically stable. Lemma 2.

If the matrix $\overline{\mathit{D}}\left(z\right)$ is invertible, the output error $\mathit{\vartheta }\left(z\right)=\mathit{d}\left(z\right)+\mathit{v}\left(z\right)$ is given by:

$$\begin{array}{l}\widehat{\mathit{\vartheta }}\left(k\right)=\mathit{y}\left(k\right)-\widehat{\mathit{s}}\left(k\right)\\ \widehat{\mathit{\vartheta }}\left(z\right)=F\left(z\right){\overline{\mathit{D}}}^{-1}\left(z\right)\mathit{e}\left(z\right)\end{array}$$ E18

where $F\left(z\right){\overline{\mathit{D}}}^{-1}\left(z\right)={\left(\mathit{I}-{\mathit{C}}^0{\left(z\mathit{I}-{\mathit{A}}^0\mathit{+}{\mathit{K}}^0{\mathit{C}}^0\right)}^{-1}{\mathit{K}}^0\right)}^{-1}$ is termed the disturbance model estimate which generates the output error when excited by the residual.

Proof: Using (17) we get (18).

As the input $\mathit{w}\left(k\right)$ driving the disturbance model is not accessible, then by substituting the actual input $\mathit{w}\left(k\right)$ by the residual $\mathit{e}\left(k\right)$ , although both are zero-mean white noise processes, only the denominator polynomial of the disturbance model can be identified. Hence the term “disturbance model estimate”.

Minimum realization of the output error model is obtained using balanced model reduction method by treating $\mathit{e}\left(z\right)$ as the input and $\widehat{\mathit{\vartheta }}\left(k\right)$ as the output of a model [7].

2.3.3 Minimal realization of the signal model

There are two approaches to identifying the signal model and the signal. One approach is by deriving them from the residual model of the Kalman filter as shown in Lemma 1 given by (12) and the other approach is to invoke the duality between the predictor form (7) and the innovation form (8) of the Kalman filter. The latter approach may be more convenient.

In view of (12), the system model $\mathit{G}\left(z\right)$ and the signal model ${\mathit{G}}_s\left(z\right)$ is derived from the identified Kalman filter (7) by simply replacing the transition matrix of the Kalman filter ${\mathit{A}}^0-{\mathit{K}}^0{\mathit{C}}^0$ by the system transition matrix ${\mathit{A}}^0$ . Lemma 3.

$$\begin{array}{l}\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)=\left({\mathit{A}}^0\mathit{+}{\mathit{K}}^0{\mathit{C}}^0,{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^{=}\right)\\ \left({\mathit{A}}_s,{\mathit{B}}_s,{\mathit{C}}_s,{\mathit{D}}_s\right)=\mathit{\text{minreal}}\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)\end{array}$$ E19

Where $\mathit{\text{minreal}}\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)$ is the minimal realization of $\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)$ .

Proof

There is duality between the predictor form (7) and the innovation form (8) of the Kalman filter [5]. The output $\mathit{y}\left(k\right)$ and the residual $\mathit{e}\left(k\right)$ are (causally) invertible. In other words, $\mathit{e}\left(k\right)$ can be generated from the output $\mathit{y}\left(k\right)$ and $\mathit{r}\left(k\right)$ using the (causally) predictor form, and $\mathit{y}\left(k\right)$ can be generated from $\mathit{e}\left(k\right)$ and $\mathit{r}\left(k\right)$ using the innovation form [3]. Moreover, $\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)$ and $\left({\mathit{A}}_s^0,{\mathit{B}}_s^0,{\mathit{C}}_s^0,{\mathit{D}}_s^0\right)$ are associated with the system transfer matrices $\mathit{G}\left(z\right)=D^{-1}\left(z\right)\mathit{N}\left(z\right)$ and ${\mathit{G}}_s\left(z\right)=D_s^{-1}\left(z\right){\mathit{N}}_s\left(z\right)$ respectively, as shown in (4), implying that $\left({\mathit{A}}_s^0,{\mathit{B}}_s^0,{\mathit{C}}_s^0,{\mathit{D}}_s^0\right)$ is a minimum realization of $\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)$ .

The minimum realization of the system $\left({\mathit{A}}_s^0,{\mathit{B}}_s^0,{\mathit{C}}_s^0,{\mathit{D}}_s^0\right)$ is obtained from the balanced model reduction method by treating $\mathit{u}\left(k\right)$ as the input and $\widehat{\mathit{s}}\left(k\right)$ as the output of a model [7].

3. Emulator-based two-stage identification

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 that system. To overcome this problem, the system model is identified by performing a number of emulator parameter-perturbed experiments proposed in [7, 8, 9]. Each experiment consists of perturbing one or more emulator parameters. A robust model is identified as the best fit to the input–output data from the set of emulated perturbations. The robust 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 (that is, the conventional approach assumes that the model of the system remains unperturbed at every operating point). In [7, 8, 9], it is theoretically shown that the identification errors resulting from the variations in the emulator parameters are significantly lower compared to those of the conventional ones based on performing a single experiment (that is, without using emulators). The emulator-based identification scheme is inspired from the model-free artificial neural network approach which captures the static and dynamic behaviors by presenting the neural network with data covering likely operating scenarios. The PEM identifies the robust model of the plant, and the Kalman filter associated with the plant is then derived from the identified model without any a-priori knowledge of the statistics, such as covariance of the disturbance and measurement noise affecting the input-output data.

An accurate emulator-based model identification scheme is proposed and employed here. An emulator, which is modeled as a product of first-order all-pass filters and which induces phase and gain changes, is connected in cascade to the input, output or both, of the signal model to emulate a set of likely operating regimes around the nominal operating point. The identified model is obtained as the best fit over all emulated operating regions, thereby ensuring both accuracy and robustness of the identified model.

3.1 Two-stage identification

• In the first stage, a robust model of the system $\left({\mathit{A}}^0{\mathit{B}}^0{\mathit{C}}^0{\mathit{D}}^0\right)$ and its associated Kalman filter $\left({\mathit{A}}^0-{\mathit{K}}^0{\mathit{C}}^0,\left[\left({\mathit{B}}^0-{\mathit{K}}^0{\mathit{D}}^0\right){\mathit{K}}^0\right],{\mathit{C}}^0{\mathit{D}}^0\right)$ are identified using PEM from the set of the emulator-generated input-output data. Then the estimate ${\mathit{s}}^0\left(k\right)$ of the signal $\mathit{s}\left(k\right)$ and the estimate $\widehat{\mathit{\vartheta }}\left(k\right)$ of the output error $\mathit{\vartheta }\left(k\right)$ are derived.

• In the second stage, using the key properties established in Lemmas 1–3, the robust signal model $\left({\mathit{A}}_s^0,{\mathit{B}}_s^0,{\mathit{C}}_s^0,{\mathit{D}}_s^0\right)$ and its associated Kalman filter $\left({\mathit{A}}_s^0-{\mathit{K}}_s^0{\mathit{C}}_s^0,\left[{\mathit{B}}_s^0-{\mathit{K}}_s^0{\mathit{D}}_s^0{\mathit{K}}_s^0\right],{\mathit{C}}_s^0,{\mathit{D}}_s^0\right)$ are obtained using balanced model reduction method and the PEM.

Akaike Information Criterion: To select an appropriate order for the identified system model in the first stage, and for the signal model in the second stage, the widely popular Akaike Information Criterion (AIC) is used, which weights both the parameter estimation error and the complexity of the model so as to arrive at an optimal order [1].

3.2 Signal model and the Kalman filter

Similar to the Kalman filter for the system (7), the Kalman filter for the signal is:

$$\begin{array}{l}{{\displaystyle \widehat{\mathit{x}}}}_s\left(k+1\right)=\left({\mathit{A}}_s^0-{\mathit{K}}_s^0{\mathit{C}}_s^0\right){{\displaystyle \widehat{\mathit{x}}}}_s\left(k\right)+\left({\mathit{B}}_s^0-{\mathit{K}}_s^0{\mathit{D}}_s^0\right)\mathit{u}\left(k\right)+{\mathit{K}}_s^0{\mathit{s}}^0\left(k\right)\\ {\displaystyle \widehat{\mathit{s}}}\left(k\right)={\mathit{C}}_s^0{{\displaystyle \widehat{\mathit{x}}}}_s\left(k\right)+{\mathit{D}}_s^0\mathit{u}\left(k\right)\\ {\mathit{e}}_s\left(k\right)={\mathit{s}}^0\left(k\right)-\widehat{\mathit{s}}\left(k\right)\end{array}$$ E20

Where ${\widehat{\mathit{x}}}_s\left(k\right)\in {\mathbb{R}}^{n_s}$ ; $\widehat{\mathit{s}}\left(k\right)\in {\mathbb{R}}^{n_s}$ ; the residual ${\mathit{e}}_s\left(k\right)={\left[\begin{array}{lllll}{e}_{s1}\left(k\right)& e_{s2}\left(k\right)& e_{s3}\left(k\right)& \dots & e_{\mathit{sq}}\left(k\right)\end{array}\right]}^T\in {\mathbb{R}}^q$ is the residual; and ${\mathit{K}}_s\in {\mathbb{R}}^{n_s\mathit{xq}}$ is the Kalman gain.

Status monitoring: The residuals $\mathit{e}\left(k\right)$ and ${\mathit{e}}_s\left(k\right)$ of the Kalman filters (7) and (20) are employed to monitor the status of the overall system and to detect and isolate faults in the signal and disturbance models and the sensors. The proposed scheme provides a sound framework for developing fault-tolerant systems and condition-based maintenance systems as well.

4. Evaluation on the illustrative example

The proposed two-stage identification scheme and the key properties of the Kalman filter established in the lemmas in Sections 2.3.1–2.3.3 are verified using the illustrative example given in Section 3.1. The results of this illustration are shown below in Figure 3a and b.

Figure 3.

(a) Signal and its estimate; output error and (b) autocorrelations of equation error and the residual and its estimate.

Subfigures A and B, of Figure 3a compare the true step response of the signal and its Kalman filter estimate; subfigures C and D show the output error $\mathit{\vartheta }\left(k\right)$ and its estimate.

Remarks: These subfigures confirm the accuracy of the estimates of the signal and the output error (18) established in Lemmas 1 and 2. Subfigures A and B, of Figure 3b show the autocorrelation of the equation error whereas subfigures C and D show the autocorrelations of the residual.

Moreover, these subfigures clearly confirm that the equation error is a colored noise that is whitened by the KF, thus confirming (17) of Lemma 1 and making the KF residual a zero-man white noise process.

Table 1 compares the true and estimated poles of the signal and disturbances models. The estimated poles are obtained from the model reduction techniques employed in the second stage of the two-stage identification scheme.

	True poles	Identified poles
Signal ${\widehat{\mathit{G}}}_s\left(z\right)$	$0.7500\pm \mathrm{j}0.3708$; $0.8500\pm j0.2784$	$0.7510\pm \mathrm{j}0.3715$; $0.8483\pm \mathrm{j}0.2769$
Disturbance ${\widehat{\mathit{G}}}_w\left(z\right)$	$0.1980\pm \mathrm{j}0.8737$; $0.5663\pm \mathrm{j}0.4114$	$0.2031\pm \mathrm{j}0.8752$; $0.5822\pm \mathrm{j}0.3746$

Table 1.

Poles of the signal and disturbance models.

Remarks: The estimated poles are close to the true ones, especially those of the signal.

5. Evaluation of the proposed scheme

The management of leakage faults in fluid systems is becoming increasingly important in recent years from the point of view of economy, potential hazards, pollution, and conservation of scarce resources. Leakage in pipes and storage tanks occurs due to faulty joints, aging, excessive loads, holes caused by corrosion and accidents and the like. The process control system is a MIMO system that exhibits turbulence and is modeled as a combination of a signal, (which includes an ideal noise-free height, flow rate, and control input), a disturbance that includes effects of turbulence, and a measurement noise. The augmented model of the signal, and the disturbance, whose output is a sum of the signal, the disturbance, and the measurement noise described by Box-Jenkins model. The transfer matrices of the signal and the disturbance may be totally different from those of the ARMA model, where the signal and the disturbance have identical denominator polynomials.

Physical systems are subject to model uncertainties and are affected to unknown stochastic disturbances such as turbulence and measurement noise. The proposed scheme covers a wider class of systems compared to the laminar flow model proposed. The laminar flow is the flow of a fluid when each particle of the fluid follows a smooth path which results in the velocity of the fluid being constant. The turbulent flow is an irregular flow that exhibits tiny whirlpool regions.

It is assumed that the disturbance is a Gaussian stochastic process and the measurement noise is a zero-mean Gaussian white noise process. The measurement output is, in general, an additive combination of the signal, disturbance and measurement noise. The output error, which is a sum of the disturbance and measurement noise, is assumed to be bounded. The signal and the disturbance are both modeled as outputs of linear time-invariant systems driven by some known input, and a Gaussian zero-mean white noise process, respectively. It is assumed that the signal, disturbance and measurement noise are mutually uncorrelated with each other.

5.1 Physical two-tank fluid system

A benchmark model is a cascade connection of a dc motor and a pump relating the reference input $r\left(t\right)$ and the flow rate $f\left(t\right)$ , the outflow $q_0\left(t\right)$ and leakage $q_{\mathrm{\ell }}\left(t\right)$ is a fourth-order system. The linearized signal model of the nonlinear SIMO system is:

$$\left[\begin{array}{l}\dot{h}\\ {\dot{h}}_2\\ \dot{u}\\ \dot{f}\end{array}\right]=\left[\begin{array}{llll}-a_1-\alpha & a_1& 0& b_1\\ a_2& -a_2-\beta & 0& 0\\ -1& 0& 0& 0\\ -b_m{k}_p& 0& b_m{k}_I& -a_m\end{array}\right]\left[\begin{array}{l}h\\ h_1\\ u\\ f\end{array}\right]+\left[\begin{array}{l}0\\ 0\\ 1\\ b_m{k}_p\end{array}\right]r\left(t\right)$$ E21

Where $h$ , $h_2$ , $u$ and $f$ are respectively the height of tank 1, the height of tank 2, the control input and the flow rate; $a_1$ , $a_2$ , $\alpha$ and $\beta$ are parameters associated with the linearization process, $\alpha$ is the leakage flow rate, $q_{\mathrm{\ell }}=\mathit{\alpha h}$ , and $\beta$ is the output flow rate, $q_o=\beta h_2$ . The output is given by:

$$\begin{array}{l}s\left(k\right)={\left[h\left(k\right) f\left(k\right) u\left(k\right)\right]}^T\\ y\left(k\right)=s\left(k\right)+d\left(k\right)+v\left(k\right)\end{array}$$ E22

5.2 Simulation of faults in the system

The process control system is interfaced to National Instruments LABVIEW as shown below in Figure 4a. The controller is implemented in LABVIEW.

Figure 4.

(a) Two-tank fluid system controlled by LABVIEW interfaced to a PC and (b) block diagram of process control system.

The two-tank system formed of subsystems and whose faults are to be isolated is shown in Figure 4b. There are four subsystems whose faults are to be isolated. Subsystem 1 is the flow rate sensor ${\gamma }_{s1}$ , subsystem 2 the height sensor ${\gamma }_{s2}$ , subsystem 3 the actuator $G_1=G_1^0{\gamma }_a$ where $G_1^0$ is the fault-free transfer function, and subsystem 4 the leakage fault sensor gain ${\gamma }_{\mathrm{\ell }}$ . The fault-free cases correspond to ${\gamma }_{\mathit{si}}=1:i=0,1,2$ , ${\gamma }_a=1$ and ${\gamma }_{\mathrm{\ell }}=1$ . The various subsystems and sensor blocks are all shown in Figure 4b. The first two blocks $G_0$ and $G_1=G_1^0{\gamma }_a$ , represent the controller and the actuator sub-systems, respectively. The leakage is modeled by the gain ${\gamma }_{\mathrm{\ell }}$ which is used to quantify the amount of flow lost from the first tank. Thus, the net outflow from tank 1 is quantified by the gain ( $1-{\gamma }_{\mathrm{\ell }}$ ). Since the two blocks $G_2^0$ and ( $1-{\gamma }_{\mathrm{\ell }}$ ) cannot be dissociated from each other, they are fused into a single block labeled $G_2=G_2^0\left(1-{\gamma }_{\mathrm{\ell }}\right)$ . The physical two-tank system is controlled using LABVIEW which acquires the flow rate, and the height sensor outputs. The controller is implemented in LABVIEW and the controller output drives the actuator, namely the DC motor and pump combination. A fault in the sensor is introduced by including the emulator block, ${\gamma }_{\mathit{si}}:i=0,1,2$ in the control input, flow rate, the height sensors, respectively in LABVIEW software. Similarly, an actuator fault is introduced by including an emulator ${\gamma }_a$ between the controller output and the input to the DC motor. The leakage fault is simulated by opening the drainage valve of the first tank. The amount by which the valve is opened is modeled by the emulator ${\gamma }_{\mathrm{\ell }}$ .

The height, flow rate, and control input profiles under various types of faults, are shown in Figure 5. Subfigures A, B and C show profiles for the leakage; subfigures D, E and F show the profiles for the actuator fault; subfigures G, H, and I show the profiles for sensor fault. The fault was simulated by varying the appropriate emulator parameters ${\gamma }_{\mathrm{\ell }}$ , ${\gamma }_a$ and ${\gamma }_{s2}$ , by 0.25, 0.5 and 0.75 times their nominal values representing ‘small’, ‘medium’ and ‘large’ fault sizes respectively.

Figure 5.

(a) Height, flow rate, control: nonlinear and (b) height, flow rate, control: linearized.

The height, the flow rate and the control input profiles under various types of faults are shown in Figure 5a for the nonlinear model, namely the dead-band effect of the actuator on the flow measurements. The measurement outputs are corrupted by the disturbance and measurement noise. Figure 5b show the outputs of the two-stage identification of the linearized signal model. Subfigures A, B and C on the top show height profiles, and subfigures D, E and F in the middle show the flow rate profiles, and G, H, and I at the bottom show the control input profiles under leakage, actuator and sensor faults, respectively. The faults are induced by varying the appropriate emulator parameters to 0.25, 0.5 and 0.75 times the nominal values to represent ‘small’, ‘medium’ and ‘large’ faults. However, by 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. Also, 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 and saturation of the flow.

Remarks: The two-stage identification is employed to estimate the height, flow rate and the control input; their estimates are shown in Figure 5b. Comparing subfigures D, E, F confirms the superior performance of the identified estimates, thanks to the use of emulators.

Figure 6a shows the residuals and their test statistics, and Figure 6b shows the autocorrelations of the residuals when the system is subject to leakage, actuator, and sensor faults of various degrees such a small, medium and large fault sizes. Subfigures A, B, and C; D, E, and F; and G, H, and I of Figure 6a show the residuals and their statistics when there is a leakage, actuator and sensor faults, respectively. Subfigures A, B, and C; D, E, and F; and G, H, and I of Figure 6b show the corresponding auto-correlations for different fault types.

Figure 6.

(a) The residuals and test statistics and (b) autocorrelations of the residuals.

Remarks: The Bayes decision strategy was employed to assert the fault type, i.e., to decide whether it is either a leakage or an actuator or sensor fault, respectively, using the fault isolation scheme proposed in [8]. The variance of the residual, which is the maximum value of the autocorrelation function evaluated at the origin (i.e. at zero delay), indicates the fault size.

The proposed Kalman-filter-based scheme can detect and isolate small and nascent faults and estimate the fault size. Thanks to the emulator-generated data, it can also provide an accurate prognosis of the status of the system.

6. Conclusions

Emulator-based identification of a wider class of multiple-input and multiple-output system governed by Box-Jenkins model and the associated Kalman filter directly from the input-output data without a-priori knowledge of the disturbance and measurement noise statistics, and the establishment of key properties of estimation of the signal, the output error and their models are developed. The applications include monitoring the status of the system including faults, distinguishing between the variations in the disturbance model and those in the signal model to help diagnose a fault in the system and ensure a low false alarm probability, developing a framework for controlling autonomous vehicles, and meeting the ever-increasing need for fault-tolerant systems. The proposed emulator-based two-stage identification and estimation of the signal and its model were evaluated physical laboratory-scale process control system so as to estimate the signal corrupted by disturbance such as turbulence. Thanks to emulator-based identification, the estimates of the signal were accurate, the detection and isolation of leakage faults were promising and, as such, provide sufficient encouragement and impetus to try the proposed scheme on real-life processes in our future work.

Acknowledgments

The first author acknowledges the support of the Department of Electrical and Computer Engineering, The University of New Brunswick, and both authors acknowledge the help and suggestions of Dr. Haris Khalid.