1 Introduction

Each technical system can be usually split into three parts: actuators, process and sensors. Each of these parts is affected by the so-called unknown inputs, which can be perceived as process and measurement noise as well as external disturbances acting on the system. The system may also be affected by faults. A fault can generally be defined as an unpermitted deviation of at least one characteristic property or parameter of the system from the normal condition, e.g., a sensor malfunction. All the unexpected variations that tend to degrade the overall performance of a system can also be interpreted as faults. Contrarily to the term failure, which suggests a complete breakdown of the system, the term fault is used to denote a malfunction rather than a catastrophe. Indeed, failure can be defined as a permanent interruption of the system ability to perform a required function under specified operating conditions. In the light of the above discussion, it is clear that the design of fault diagnosis schemes that prevent turning faults into failures is of paramount importance. One way to settle a challenging problem of fault diagnosis is to use model-based redundancy [3, 8, 27]. In this case, the system model quality determines the effectiveness of the Fault Detection and Isolation (FDI) [3, 8, 13, 18, 28], and consequently, the Fault-Tolerant Control (FTC) [2, 19, 20, 27].

An obvious way the obtain the model is to employ the physical relation governing the investigated system. As a result, the so called analytical model is obtained. Unfortunately, the complexity of modern industrial systems usually makes the model derivation difficult or even impossible (in the light of the quality of the achieved model). In the case of non-linear dynamic system, the Artificial Neural Networks (ANNs) constitute an elegant remedy to the above-mentioned problem [4]. Unfortunately, the ANNs have disadvantages, e.g., they are usually not available in the state-space form [11, 22, 30] frequently used for fault diagnosis. Moreover, only rare approaches ensure the stability [21] and there is a limited number of solutions that can settle the robustness problems regarding neural model uncertainty [15, 26].

This issue is very important for the model-based FDI systems which are usually based on the residual generation and constant threshold application. Neglecting the model uncertainty and measurements noise [12, 15] in the FDI system, may result in the undetected faults or false alarms. To solve such a challenging problem, a methodology of dynamic non-linear system identification on the basis of the state-space Group Method of Data Handling (GMDH) neural network [14] was proposed. Such a neural model is gradually constructed by the connection of the partial models (neurons) with the application of the appropriated selection methods [15], what result in the significant reduction of the neural model inaccuracy. The application of the Unscented Kalman Filter (UKF) [25] during the training of the neural model allows to obtain the neurons parameters estimate and the corresponding description of the neural model uncertainty. Such knowledge is necessary to calculate the neural model output adaptive thresholds which allow to perform the robust sensor fault detection [15].

Unfortunately, the above method can be only applied for sensors but not for the actuator fault detection, which means that the neural network works as a virtual sensor parallel to the one present in the system. Then, the measurements provided by the system sensor and those by the virtual sensor are compared to perform the diagnostics decisions. In order to solve such a disadvantage in this paper it is shown how to use the above-mentioned results for actuator fault diagnosis [1, 24]. In particular, the obtained GMDH neural model is used to obtain the system input estimates and the corresponding input adaptive thresholds. The achieved thresholds are then compared with the original system input signals to formulate the diagnostics decisions. The input estimation scheme is based on a chain of Robust Unknown Input Filters (RUIF) [9, 29, 31], which guaranties that the input estimates are obtained with a prescribed attenuation level while ensuring the convergence of the observers. Thus, the complete solution that enables to perform robust fault detection and isolation of the actuator fault is delivered.

The paper is organised as follows. Section 2 portrays a general scheme of the state-space GMDH neural networks along with the associated sensor fault diagnosis scheme. It also outlines the approach that can be used for estimating the parameters of such a neural network. Section 3 presents the main contribution of the paper resulting in the robust actuator fault detection and isolation scheme. Section 4 illustrates the application of the proposed approach in robust fault detection of the tunnel furnace. Finally, the last section is devoted to conclusions.

2 Robust Sensor Fault Detection with the GMDH Neural Network

The effectiveness of the model-based fault detection system mainly depends on the quality of the model of the diagnosed system, which is obtained during system identification. Several methods for improving the neural model quality can be found in the literature. However, it should be underlined that irrespective of the used identification method the neural models will never ideally mimic the identified system. Thus, the robustness of the fault detection system against model uncertainty is one of the most desirable features. The robust fault detection system requires the knowledge about the uncertainty of the model. Based on the mathematical description of the model uncertainty it is possible to calculate the output adaptive thresholds which, allow performing robust sensor fault detection according to the scheme presented in Fig. 1. The output adaptive thresholds should contain real system responses in fault-free mode. Note that in the remaining part of this section it is assumed that the system along with all actuators are fault-free. Under such an assumption it s possible to use the presented scheme for sensor fault detection and isolation. An occurrence of the sensor faults is signaled when system outputs \(\varvec{y}_k\) cross the output adaptive threshold:

$$\begin{aligned} \hat{y}_{i,k}^{m} \le y_{i,k} \le \hat{y}_{i,k}^{M}, \end{aligned}$$
(1)

where \(\hat{y}_{i,k}^{m}\) and \(\hat{y}_{i,k}^{M}\) denote the minimum and maximum value of the adaptive threshold for the \(i\)-th system output.

Fig. 1
figure 1

Scheme of the robust sensor fault detection

Full size image

The model of the diagnosed system can be obtained with the application of the GMDH approach. Such a method allows to identify non-linear dynamic system along with the description of its uncertainty. Moreover, synthesis process of such a neural network allows obtaining a model with relatively small uncertainty, which increases the sensitivity of the fault detection system. The subsequent steps of procedure of the GMDH neural model synthesis procedure [5, 14, 16, 17, 23] are presented in Fig. 2.

Fig. 2
figure 2

Synthesis process of the GMDH neural network

Full size image

During the GMDH neural model synthesis new layers of neurons are added to the network until the quality of the neural model is evaluated by the suitable criteria [15, 16] is increasing (cf. Fig. 3).

Fig. 3
figure 3

Synthesis process of the GMDH neural network

Full size image

It is worth to emphasis that during the neural model training the parameters of each neuron in the GMDH network are estimated separately. Moreover, the neurons parameters are estimated in such a way to ensure the slightest possible uncertainty. It is possible by the application of the appropriate learning algorithm and assuming adequate structure of the neuron. For this reason in the paper the following form of the state-space neuron is proposed:

$$\begin{aligned} \hat{\varvec{x}}_{k+1}^{}&=\varvec{A}\hat{\varvec{x}}_{k}^{}+\varvec{B}\varvec{u}_{k}^{}, \end{aligned}$$
(2)
$$\begin{aligned} \hat{\varvec{y}}_{k+1}^{}&=\varvec{g}_{}(\varvec{C}\hat{\varvec{x}}_{k+1}^{}), \end{aligned}$$
(3)

where \(\varvec{u}_{k}^{}\in \mathbb {R}^{n_u}\) and \(\varvec{y}_{k}^{}\in \mathbb {R}^{n_y}\) represent the inputs and outputs of the dynamic neuron created on the combination of systems inputs. \(\varvec{g}_{}(\cdot )=[g_1(\cdot ),...,g_{n_y}(\cdot )]^{T}\) where \(g_i(\cdot )\) denotes a non-linear activation functions. \(\varvec{A}\in \mathbb {R}^{n_x\times n_x}\), \(\varvec{B}\in \mathbb {R}^{n_x\times n_u}\), \(\varvec{C}\in \mathbb {R}^{n_y\times n_x}\) and \(\varvec{x}_{k}^{}\in \mathbb {R}^{n_x}\). As the matrix \(\varvec{A}\) has an upper-triangular form it means that the neuron is asymptotically stable iff all diagonal elements of matrix \(\varvec{A}\) fulfill the condition:

$$\begin{aligned} \mid a_{i,i} \mid <1, \quad i=1,...,n_z. \end{aligned}$$
(4)

Such a neuron model clearly determines the class of systems for which the proposed neural network can be used. Thus, an assumption underlying further deliberations is that the system can be modeled with a specific network structure composed of the neurons described by (2)–(3).

As it was already mentioned, the parameters of each neuron in the GMDH neural network are estimated separately. This property allows to apply the UKF in the process of the synthesis of the GMDH neural network. The UKF employs the unscented transform [6], which approximates the mean \(\hat{\varvec{y}}_{k}^{}\in \mathbb {R}^{n_y}\) and covariance \(\varvec{P}^{yy}_k\in \mathbb {R}^{n_y\times n_y}\) of so-called transformed sigma points after the non-linear transformation \(\varvec{y}_{k}^{}=\mathcal {H}_{}(\varvec{x}_{k})\), where the mean and covariance of sigma points are given as \(\hat{\varvec{x}}_{k}^{}\in \mathbb {R}^n\) and \(\varvec{P}^{xx}_k\in \mathbb {R}^{n\times n}\) (Fig. 4). The UKF [7] can be perceived a derivative-free alternative to the Extended Kalman Filter (EKF) in the framework of the state-estimation. One of the main advantage of application of the UKF to the constrained parameter estimation is that the asymptotically stable neurons are obtained. It should be underlined that the state-space GMDH neural model has a cascade structure and is asymptotically stable, when each of neurons in the network is asymptotically stable [10]. Moreover, the application of the UKF with the procedure of truncation of the probability density function [25] allows obtaining the uncertainty of the GMDH model in the form of a covariance matrix \(\varvec{P}^{xxt}\). Such knowledge allows to calculate the system output adaptive thresholds which should contain the real system responses in the fault-free mode:

$$\begin{aligned} \mathcal {F}_{i}(\varvec{c}_i\hat{\varvec{x}}_{k}^{}-t^{\alpha /2}_{n_t-n_p}\hat{\sigma }_i\sqrt{\varvec{c}_i\varvec{P}^{xxt}\varvec{c}_i^T}) \le y_{i,k} \le \mathcal {F}_{i}(\varvec{c}_i\hat{\varvec{x}}_{k}^{}+t^{\alpha /2}_{n_t-n_p}\hat{\sigma }_i\sqrt{\varvec{c}_i\varvec{P}^{xxt}\varvec{c}_i^T}). \end{aligned}$$
(5)

where \(\varvec{c}_i\) stands for the \(i\)-th row (\(i=1,...,n_y\)) of the matrix \(\varvec{C}{}\) of the output neuron, \(\hat{\sigma }_i\) represents the standard deviation of the \(i\)-th fault-free residual and \(t^{\alpha /2}_{n_t-n_p}\) denotes the t Student distribution quantile. The sensors faults are signaled when system outputs \(\varvec{y}_k\) crosses the output adaptive threshold (5) what is presented in Fig. 5. Finally, it is also worth to underline, that the faulty sensors can be replaced by the response of a neural network, which can be perceived as a virtual sensor.

Fig. 4
figure 4

Scheme of the UKF algorithm

Full size image
Fig. 5
figure 5

Robust sensors fault detection with the application of the output adaptive thresholds obtained via the UKF

Full size image

3 Robust Actuators Fault Detection and Isolation with the GMDH Neural Network and RUIF

The approach presented in Sect. 2 enables performing robust sensor fault detection. Unfortunately, it does not allow to detect and isolate the faulty actuator. To solve such a challenging problem, a method depicted in Fig. 6 is to be suitably developed in the subsequent part of this section. To achieve this goal it is necessary to develop a methodology for calculating the input adaptive threshold:

$$\begin{aligned} \hat{u}^{m}_{i,k} \le u_{i,k} \le \hat{u}^{M}_{i,k}, \end{aligned}$$
(6)

where \(\hat{u}^{m}_{i,k}\) and \(\hat{u}^{M}_{i,k}\) represent the minimum and maximum value of the adaptive threshold for the \(i\)-th system input of the diagnosed system.

Fig. 6
figure 6

Application of the GMDH neural model and RUIF to robust fault detection and isolation of actuators

Full size image

The state-space description of the GMDH network neuron allows to develop a new RUIF-based approach which enables to estimate the input signals of the GMDH neural model (cf. Fig. 7), and in the consequence calculating the input adaptive thresholds for the robust fault diagnosis of the actuators. Let us consider a non-linear discrete-time neuron (2–3):

Fig. 7
figure 7

Estimation of the system inputs via GMDH model and RUIF

Full size image
$$\begin{aligned} \varvec{x}_{k+1}^{}&=\varvec{A}\varvec{x}_{k}^{}+\varvec{B}\varvec{u}_{k}^{}+ \varvec{W}_{1}\varvec{w}_{k}, \end{aligned}$$
(7)
$$\begin{aligned} \varvec{y}_{k+1}^{}&=\varvec{g}_{}(\varvec{C}\varvec{x}_{k+1}^{}) + \varvec{W}_{2}\varvec{w}_{k+1}, \end{aligned}$$
(8)

where \(\varvec{W}_{1}\) and \(\varvec{W}_{2}\) are known disturbance distribution matrices, \(\varvec{w}_{k} \in l_2\) is an exogenous disturbance vector, \(l_2=\left\{ \varvec{w}_{} \in \mathbb {R}^n, \Vert \varvec{w}_{}\Vert _{l_2}< + \infty , \right\} \) where \(\Vert \varvec{w}_{}\Vert _{l_2}=\left( \sum _{k=0}^{\infty }\Vert \varvec{w}_{k}\Vert ^2\right) ^{\frac{1}{2}}\). Thus, the neuron can be perceived as a system with unknown inputs. Note also that \(\varvec{w}_{k}\) may represent various sources of uncertainty, including modelling uncertainty.

Subsequently, the system output can be written as follows:

$$\begin{aligned}&\varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{}-\varvec{W}_{2}\varvec{w}_{k+1})=\varvec{C}\varvec{x}_{k+1}^{}= \varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{})+\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{},\end{aligned}$$
(9)
$$\begin{aligned}&\qquad \quad \varvec{H}[\varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{})+\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{}]=\varvec{H}\!\varvec{C}\varvec{x}_{k+1}^{}, \end{aligned}$$
(10)

where \(\varvec{v}_{k}^{} \in l_2\), and \(\tilde{\varvec{W}}_{2}\) stands for the distribution matrix of \(\varvec{v}_{k}^{}\) and must be determined by the designer. Substituting (7) into (10):

$$\begin{aligned} \varvec{H}[\varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{})+\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{}]=\varvec{H}\!\varvec{C}\!\varvec{A}\varvec{x}_{k}^{}+\varvec{H}\!\varvec{C}\!\varvec{B}\varvec{u}_{k}^{}+\varvec{H}\!\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k}, \end{aligned}$$
(11)

and satisfying \(\varvec{H}\!\varvec{C}\!\varvec{B}=\varvec{I}\), i.e. \(\varvec{H}=(\varvec{C}\varvec{B})^+\), which implies that \(\mathrm {rank}(\varvec{C}\varvec{B})=\mathrm {rank}(\varvec{B})=n_u\), the system input receives the following form:

$$\begin{aligned} \varvec{u}_{k}^{} =\varvec{H}\varvec{g}_{}^{-1} (\varvec{y}_{k+1}^{})+\varvec{H}\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{}-\varvec{H}\!\varvec{C}\!\varvec{A}\varvec{x}_{k}^{}-\varvec{H}\!\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k}. \end{aligned}$$
(12)

Based on (12), the input estimate can be defined as:

$$\begin{aligned} \hat{\varvec{u}}_{k}^{} =\varvec{H}\varvec{g}_{}^{-1} (\varvec{y}_{k+1}^{})-\varvec{H}\!\varvec{C}\!\varvec{A}\hat{\varvec{x}}_{k}^{}. \end{aligned}$$
(13)

The input estimation error can be defined as follows:

$$\begin{aligned} \varepsilon _{u,k}^{}=\varvec{u}_{k}^{}-\hat{\varvec{u}}_{k}^{}=-\varvec{H}\!\varvec{C}\!\varvec{A}\varvec{e}_{k}^{}+\varvec{H}\tilde{\varvec{w}}_{2}\varvec{v}_{k+1}^{}-\varvec{H}\!\varvec{C}\varvec{w}_{1}\varvec{w}_{k}. \end{aligned}$$
(14)

Substituting (12) into (7) gives:

$$\begin{aligned} \varvec{x}_{k+1}^{}&=\varvec{A}\varvec{x}_{k}^{}+\varvec{B}\!\varvec{H}\varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{}) +\varvec{B}\!\varvec{H}\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{}\nonumber \\&\quad -\varvec{B}\!\varvec{H}\!\varvec{C}\!\varvec{A}\varvec{x}_{k}^{}-\varvec{B}\!\varvec{H}\!\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k}, \end{aligned}$$
(15)

and denoting \(\bar{\varvec{A}}=\varvec{A}-\varvec{B}\!\varvec{H}\!\varvec{C}\!\varvec{A}\) and \(\bar{\varvec{B}}=\varvec{B}\!\varvec{H}\), (7) yield:

$$\begin{aligned} \varvec{x}_{k+1}^{}=\bar{\varvec{A}}\varvec{x}_{k}^{}+\bar{\varvec{B}}\varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{}) +\bar{\varvec{B}}\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{}-\bar{\varvec{B}}\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k}. \end{aligned}$$
(16)

Consequently, the robust unknown input observer structure is:

$$\begin{aligned} \hat{\varvec{x}}_{k+1}^{}=\bar{\varvec{A}}\hat{\varvec{x}}_{k}^{}+\bar{\varvec{B}}\varvec{g}_{}^{-1}(\varvec{y}_{k+1}^{}) +\varvec{K}(\varvec{g}_{}^{-1}(\varvec{y}_{k}^{})-\varvec{C}\hat{\varvec{x}}_{k}^{}), \end{aligned}$$
(17)

while the state estimation error is given by:

$$\begin{aligned} \varvec{e}_{k+1}^{}&= \varvec{x}_{k+1}^{}-\hat{\varvec{x}}_{k+1}^{}=\bar{\varvec{A}}\varvec{e}_{k}^{}+ \bar{\varvec{B}}\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{} \nonumber \\&-\bar{\varvec{B}}\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k} -\varvec{K}(\varvec{g}_{}^{-1}(\varvec{y}_{k}^{})-\varvec{C}\hat{\varvec{x}}_{k}^{})\nonumber \\&= \bar{\varvec{A}}\varvec{e}_{k}^{}+\bar{\varvec{B}}\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{}-\bar{\varvec{B}}\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k} -[\varvec{K}(\varvec{C}\varvec{x}_{k}^{}-\tilde{\varvec{W}}_{2}\varvec{v}_{k}^{}-\varvec{C}\hat{\varvec{x}}_{k}^{})]\nonumber \\&= \varvec{A}_1\varvec{e}_{k}^{} +\bar{\varvec{B}}\tilde{\varvec{W}}_{2}\varvec{v}_{k+1}^{} -\bar{\varvec{B}}\varvec{C}\!\varvec{W}_{1}\varvec{w}_{k} +\varvec{K}\tilde{\varvec{W}}_{2}\varvec{v}_{k}^{}\nonumber \\&= \varvec{A}_1\varvec{e}_{k}^{} +[\varvec{0} \bar{\varvec{B}}\tilde{\varvec{W}}_{2}]\bar{\varvec{v}}_{k+1}^{} + [-\bar{\varvec{B}}\varvec{C}\!\varvec{W}_{1} \varvec{0}]\bar{\varvec{v}}_{k}^{} + [\varvec{0} -\varvec{K}\tilde{\varvec{W}}_{2}]\bar{\varvec{v}}_{k}^{}\nonumber \\&= \varvec{A}_1\varvec{e}_{k}^{}+\bar{\varvec{V}}_{1}\bar{\varvec{v}}_{k+1}^{}+\bar{\varvec{V}}_{2}\bar{\varvec{v}}_{k}^{}, \end{aligned}$$
(18)

where: \(\bar{\varvec{v}}_{k}^{}=[\varvec{w}_{k}^T \varvec{v}_{k}^{T}]^{T}\), \(\varvec{A}_1=\bar{\varvec{A}}-\varvec{K}\varvec{C}\), \(\bar{\varvec{V}}_{1}=[\varvec{0} \bar{\varvec{B}}\tilde{\varvec{W}}_{2}]\), \(\bar{\varvec{V}}_{2}=[-\bar{\varvec{B}}\varvec{C}\!\varvec{W}_{1} -\varvec{K}\tilde{\varvec{W}}_{2}]\).

As a consequence, the input estimation error (14) can be redefined as follows:

$$\begin{aligned} \varepsilon _{u,k}^{}&= -\varvec{H}\!\varvec{C}\!\varvec{A}\varvec{e}_{k}^{}+ [\varvec{0} \varvec{H}\tilde{\varvec{W}}_{2}]\bar{\varvec{v}}_{k+1}^{} + [-\varvec{H}\!\varvec{C}\!\varvec{W}_{1} \varvec{0}] \bar{\varvec{v}}_{k}^{}\nonumber \\&= -\varvec{H}\!\varvec{C}\!\varvec{A}\varvec{e}_{k}^{}+ \tilde{\varvec{V}}_{1}\bar{\varvec{v}}_{k+1}^{} + \tilde{\varvec{V}}_{2} \bar{\varvec{v}}_{k}^{}. \end{aligned}$$
(19)

where \(\tilde{\varvec{V}}_{1}=[\varvec{0} \varvec{H}\tilde{\varvec{W}}_{2}]\) and \(\tilde{\varvec{V}}_{2}=[-\varvec{H}\!\varvec{C}\!\varvec{W}_{1} \varvec{0}]\).

Theorem 1

For a prescribed disturbance attenuation level \(\upsilon >0\), for the input estimation error (19), the \(\mathcal {H}_{\infty }\) observer design problem for the system (7–8) and the observer (17) is solvable if there exist \(\alpha >0\), \(\beta >0\), \(\varvec{P}\succ \varvec{0}\), \(\varvec{Q}\succ \varvec{0}\), \(\varvec{N}\) such that the following LMI is satisfied:

$$\begin{aligned} \begin{array}{r} \left[ \begin{array}{cc} -\varvec{P}+\varvec{A}^T\varvec{C}^T\varvec{H}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{2} \\ -\tilde{\varvec{V}}_{2}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} \tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{2}^T-\mu ^2\varvec{I} \\ -\tilde{\varvec{V}}_{1}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} \tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{2} \\ \varvec{P}\varvec{A}_1 &{} \varvec{P}\bar{\varvec{V}}_{2} \\ \end{array} \right. \left. \begin{array}{cc} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{1} &{} \varvec{A}_1^T\varvec{P}\\ \tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{1} &{} \bar{\varvec{V}}_{2}^T\varvec{P}\\ \tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{1}-\mu ^2\varvec{I} &{} \bar{\varvec{V}}_{1}^T\varvec{P}\\ \varvec{P}\bar{\varvec{V}}_{1} &{} -\varvec{P}\\ \end{array} \right] <\varvec{0}. \end{array} \end{aligned}$$
(20)

Proof

The problem of \(\mathcal {H}_{\infty }\) observer design [31] is to determine the gain matrix \(\varvec{K}\) such that

$$\begin{aligned} \lim _{k \rightarrow \infty } \varvec{e}_{k}^{}=\varvec{0} \quad \mathrm {for} \quad \bar{\varvec{v}}_{k}^{}=\varvec{0}, \end{aligned}$$
(21)

and

$$\begin{aligned} \Vert \varepsilon _{u,k}^{}\Vert _{l_2} \le \upsilon \Vert \bar{\varvec{v}}_{}^{}\Vert _{l_2} \quad \mathrm {for} \quad \bar{\varvec{v}}_{k}^{} \ne \varvec{0} \quad \mathrm {and} \quad \varvec{e}_{0}^{}=\varvec{0}. \end{aligned}$$
(22)

In order to settle the above problem it is sufficient to find a Lyapunov function \(V_k\) such that:

$$\begin{aligned} \Delta V_k + \varepsilon _{u,k}^{T}\varepsilon _{u,k}^{} -\mu ^2\bar{\varvec{v}}_{k+1}^{T}\bar{\varvec{v}}_{k+1}^{} -\mu ^2\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{} <0, \end{aligned}$$
(23)

where:

$$\begin{aligned} \Delta V_k&= V_{k+1}-V_{k}=\varvec{e}_{k+1}^{T}\varvec{P}\varvec{e}_{k+1}^{}-\varvec{e}_{k}^{T}\varvec{P}\varvec{e}_{k}^{}\nonumber \\&= [\varvec{e}_{k}^{T}\varvec{A}_1^T+\bar{\varvec{v}}_{k+1}^{T}\bar{\varvec{V}}_{1}^T+\bar{\varvec{v}}_{k}^{T}\bar{\varvec{V}}_{2}^T] \varvec{P}[\varvec{A}_1\varvec{e}_{k}^{}+\bar{\varvec{V}}_{1}\bar{\varvec{v}}_{k+1}^{}+\bar{\varvec{V}}_{2}\bar{\varvec{v}}_{k}^{}]\nonumber \\&-\varvec{e}_{k}^{T}\varvec{P}\varvec{e}_{k}^{}, \end{aligned}$$
(24)

Thus, if \(\bar{\varvec{v}}_{k}^{}=\varvec{0}\), (\(k=0,\dots ,\infty \)) then (23) boils down to

$$\begin{aligned} \Delta V_k+\varepsilon _{u,k}^{T}\varepsilon _{u,k}^{}<0,\, k=0,\dots \infty , \end{aligned}$$
(25)

and hence \(\Delta V_k<0\), which leads to \(\lim _{k \rightarrow \infty } \varvec{e}_{k}^{}=\varvec{0}\) for \(\bar{\varvec{v}}_{k}^{} = \varvec{0}\). If \(\bar{\varvec{v}}_{k}^{}\ne \varvec{0}\), (\(k=0,\dots ,\infty \)) then inequality (23) yields:

$$\begin{aligned} J=\sum _{k=0}^{\infty }\left( \Delta V_k+\varepsilon _{u,k}^{T}\varepsilon _{u,k}^{}-\mu ^2\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{}-\mu ^2\bar{\varvec{v}}_{k+1}^{}\bar{\varvec{v}}_{k+1}^{}\right) <0, \end{aligned}$$
(26)

which can be written as:

$$\begin{aligned} J=-V_{0}+\sum _{k=0}^{\infty }\varepsilon _{u,k}^{T}\varepsilon _{u,k}^{}-\mu ^2\sum _{k=0}^{\infty }\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{}-\mu ^2\sum _{k=0}^{\infty }\bar{\varvec{v}}_{k+1}^{T}\bar{\varvec{v}}_{k+1}^{}<0. \end{aligned}$$
(27)

Bearing in mind that:

$$\begin{aligned} \mu ^2\sum _{k=0}^{\infty }\bar{\varvec{v}}_{k+1}^{T}\bar{\varvec{v}}_{k+1}^{}=\mu ^2\sum _{k=0}^{\infty }\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{}-\mu ^2\bar{\varvec{v}}_{0}^{T}\bar{\varvec{v}}_{0}^{} \end{aligned}$$
(28)

inequality (27) can be written as:

$$\begin{aligned} J=-V_{0}+\sum _{k=0}^{\infty }\varepsilon _{u,k}^{T}\varepsilon _{u,k}^{}-2\mu ^2\sum _{k=0}^{\infty }\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{}+\mu ^2\bar{\varvec{v}}_{0}^{T}\bar{\varvec{v}}_{0}^{}<0. \end{aligned}$$
(29)

Knowing that \(V_0=0\) for \(\varvec{e}_{0}^{}=0\), (29) leads to \(\Vert \varepsilon _{k,u}^{}\Vert _{l_2} \le \upsilon \Vert \bar{\varvec{v}}_{}^{}\Vert _{l_2}\) with \(v=\sqrt{2}\mu \).

Thus, for \(\varvec{z}_{k}^{}=[\varvec{e}_{k}^{T}, \bar{\varvec{v}}_{k}^{T}, \bar{\varvec{v}}_{k+1}^{T}]^T\) the inequality (23) becomes:

$$\begin{aligned} \varvec{z}_{k}^{T}\varvec{X}_{}^{}\varvec{z}_{k}^{}<0, \end{aligned}$$
(30)

where the matrix \(\varvec{X}_{}^{}\prec 0\) has the following form:

$$\begin{aligned} \begin{array}{r} \left[ \begin{array}{c@{\quad }c} \varvec{A}_1^T \varvec{P}\varvec{A}_1 -\varvec{P}+ \varvec{A}^T\varvec{C}^T\varvec{H}^T\varvec{H}\!\varvec{C}\!\varvec{A} &{} \varvec{A}_1^T \varvec{P}\bar{\varvec{V}}_{2} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{2} \\ \bar{\varvec{V}}_{2}^T\varvec{P}\varvec{A}_1 -\tilde{\varvec{V}}_{2}^T\varvec{H}\!\varvec{C}\!\varvec{A} &{} \bar{\varvec{V}}_{2}^T\varvec{P}\bar{\varvec{V}}_{2}+\tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{2}-\mu ^2\varvec{I} \\ \bar{\varvec{V}}_{1}^T\varvec{P}\varvec{A}_1 -\tilde{\varvec{V}}_{1}^T\varvec{H}\!\varvec{C}\!\varvec{A} &{} \bar{\varvec{V}}_{1}^T\varvec{P}\bar{\varvec{V}}_{2}+\tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{2} \\ \end{array} \right. \\ \\ \left. \begin{array}{cc} \varvec{A}_1^T\varvec{P}\bar{\varvec{V}}_{1} -\varvec{A}^T\varvec{C}^T \varvec{H}^T\tilde{\varvec{V}}_{1} \\ \bar{\varvec{V}}_{2}^T\varvec{P}\bar{\varvec{V}}_{1} +\tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{1} \\ \bar{\varvec{V}}_{1}^T\varvec{P}\bar{\varvec{V}}_{1}+\tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{1}-\mu ^2\varvec{I}\\ \end{array} \right] <\varvec{0}, \end{array} \end{aligned}$$
(31)

Moreover, by applying the Schur complements, (31) is equivalent to

$$\begin{aligned} \begin{array}{r} \left[ \begin{array}{c@{\quad }c} -\varvec{P}+\varvec{A}^T\varvec{C}^T\varvec{H}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{2} \\ -\tilde{\varvec{V}}_{2}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} \tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{2}^T-\mu ^2\varvec{I} \\ -\tilde{\varvec{V}}_{1}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} \tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{2} \\ \varvec{A}_1 &{} \varvec{V}_{2} \\ \end{array} \right. \left. \begin{array}{c@{\quad }c} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{1} &{} \varvec{A}_1^T \\ \tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{1} &{} \bar{\varvec{V}}_{2}^T \\ \tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{1}-\mu ^2\varvec{I} &{} \bar{\varvec{V}}_{1}^T \\ \varvec{V}_{1} &{} -\varvec{P}^{-1} \\ \end{array} \right] <\varvec{0}, \end{array} \end{aligned}$$
(32)

Multiplying (32) from both sites by \(\mathrm {diag}( \varvec{I}, \varvec{I}, \varvec{I}, \varvec{P})\), and then substituting \(\varvec{A}_1=\bar{\varvec{A}}-\varvec{K}\varvec{C}\), \(\varvec{P}\varvec{A}_1=\varvec{P}\bar{\varvec{A}}-\varvec{P}\varvec{K}\varvec{C}=\varvec{P}\bar{\varvec{A}}-\varvec{N}\varvec{C}\), \(\varvec{A}_1^T\varvec{P}=\bar{\varvec{A}}^{T}\varvec{P}-\varvec{C}^T\varvec{N}^T\) and \(\varvec{N}=\varvec{P}\varvec{K}\), \(\varvec{P}\bar{\varvec{V}}_{2}=\varvec{P}[\varvec{0} -\varvec{K}\tilde{\varvec{W}}_{2}]=[\varvec{0} -\varvec{P}\varvec{K}\tilde{\varvec{W}}_{2}]=[\varvec{0} -\varvec{N}\tilde{\varvec{W}}_{2}]\), (32) receives the form:

$$\begin{aligned} \begin{array}{r} \left[ \begin{array}{c@{\quad }c} -\varvec{P}+\varvec{A}^T\varvec{C}^T\varvec{H}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{2} \\ -\tilde{\varvec{V}}_{2}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} \tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{2}^T-\mu ^2\varvec{I} \\ -\tilde{\varvec{V}}_{1}^T\varvec{H}\!\varvec{C}\!\varvec{A}&{} \tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{2} \\ \varvec{P}\varvec{A}_1 &{} \varvec{P}\bar{\varvec{V}}_{2} \\ \end{array} \right. \left. \begin{array}{c@{\quad }c} -\varvec{A}^T\varvec{C}^T\varvec{H}^T\tilde{\varvec{V}}_{1} &{} \varvec{A}_1^T\varvec{P}\\ \tilde{\varvec{V}}_{2}^T\tilde{\varvec{V}}_{1} &{} \bar{\varvec{V}}_{2}^T\varvec{P}\\ \tilde{\varvec{V}}_{1}^T\tilde{\varvec{V}}_{1}-\mu ^2\varvec{I} &{} \bar{\varvec{V}}_{1}^T\varvec{P}\\ \varvec{P}\bar{\varvec{V}}_{1} &{} -\varvec{P}\\ \end{array} \right] <\varvec{0}, \end{array} \end{aligned}$$
(33)

which completes the proof. \(\square \)

As the result of solving of Linear Matrix Inequality (LMI) (33), for a given disturbance attenuation level \(\mu \) the observer gain matrix \(\varvec{K}\) can be obtained:

$$\begin{aligned} \varvec{K}=\varvec{P}^{-1}\varvec{N}\end{aligned}$$
(34)

The above-presented methodology allows calculate estimates of GMDH neural network inputs. Furthermore, on the basis of (23):

$$\begin{aligned} \varepsilon _{u,k}^{T}\varepsilon _{u,k}^{} \le \mu ^2\bar{\varvec{v}}_{k+1}^{T}\bar{\varvec{v}}_{k+1}^{} +\mu ^2\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{}, \end{aligned}$$
(35)

and assuming that \(\bar{\varvec{v}}_{k}^{T}\bar{\varvec{v}}_{k}^{}=\Vert \bar{\varvec{v}}_{k}^{} \Vert _2^2<\delta \), where \(\delta >0\) is a given bound, then

$$\begin{aligned} \varepsilon _{u,k}^{T}\varepsilon _{u,k}^{} \le 2\mu ^2\delta , \end{aligned}$$
(36)

and the adaptive threshold for the inputs of the GMDH neural model can be defined sa follows:

$$\begin{aligned} \hat{u}_{i,k}- \mu \sqrt{2\delta } \le u_{i,k} \le \hat{u}_{i,k}+\mu \sqrt{2\delta }. \end{aligned}$$
(37)

During the actuator fault diagnosis, an occurrence of the fault of the \(i\)-th actuator is signaled when input \(u_{i,k}\) crosses the input threshold (37) (cf. Fig. 8).

Fig. 8
figure 8

Robust actuators fault detection and isolation with the application of the input adaptive thresholds obtained via the RUIF

Full size image

Note that the proposed scheme depicted in Fig. 6 can be used for detecting and isolating actuator under an assumption that all sensors are fault-free. To relax such an assumption it is necessary to employ the scheme presented in Fig. 9 for which the sensitivity to actuator and sensor fault is described in Table 1. The notation being used is as follows:

  • \(y_k^i\) denotes the output vector \(y_k\) without \(i\)th element,

  • \(f_{a,i}\) denotes \(i\)th actuator fault,

  • \(f_{s,i}\) stands for \(i\)th sensor fault,

  • \(\textit{RUIF}_i\) denotes the \(i\)th neural model supported with RUIF.

After providing appropriate nomenclature it is possible to analyze the fault sensitivity matrix expressed by Table 1. It is evident that it makes it possible to determine if there is either an actuator or sensor fault. If the sensor a symptom of sensor fault is obtained then the approach presented in Sect. 2 can be used for further analysis else the actuator fault isolation scheme presented in Fig. 9 should be employed.

Fig. 9
figure 9

Scheme of the actuator fault isolation

Full size image
Table 1 Fault sensitivity matrix
Full size table

4 Illustrative example

The objective of this section is to design a dynamic GMDH model of the tunnel furnace (cf. Figs. 10, 11) and apply it to the robust fault detection of the actuator and sensor fault with the input and output adaptive thresholds developed in Sects. 2 and 3.

Fig. 10
figure 10

Laboratory model of a tunnel furnace

Full size image
Fig. 11
figure 11

Interior of a tunnel furnace

Full size image

In the laboratory conditions the tunnel furnace enables to mimic the real industrial tunnel furnaces, which can be applied in the food industry. It is equipped in three electric heaters representing actuators and four temperature sensors. The temperature of the furnace can be controlled by a group regulation of voltage in the heaters with the application of the controller PACSystems RX3i manufactured by GE Fanuc Intelligent Platforms and semiconductor relays RP6 produced by LUMEL providing an impulse control with a variable impulse frequency \(fmax = 1\) Hz. The maximum power outputs of the heaters were measured to be approximately 686, 693 and \(756 \pm 20\) W, respectively. The temperature of the furnace is measured via IC695ALG600 module with Pt100 Resistive Thermal Devices (RTDs) with an accuracy of \(\pm 0.7^{\circ }\) C. The visualisation of the behaviour of the tunnel furnace is made by Quickpanel CE device from GE Fanuc Intelligent Platforms. It should be underlined that the considered system is a distributed parameter one (i.e., a system whose state space is infinite dimensional), thus any resulting model from input-output data will be at best an approximation.

The modeled furnace is a three-input and four-output system \((t_1,t_2,t_3,t_4)\!=\!f(u_1,u_2,u_3)\), where \(t_1,\ldots ,t_4\) represents the temperatures from sensors and \(u_1,\ldots ,u_3\) denotes input voltages of the electric heaters. For the modeling of the tunnel furnace purpose the state-space GMDH neural model approach [14] was applied. The parameters of the state-space dynamic neurons are estimated with the application of the UJF training algorithm [25]. The selection of the best performing neurons is realized with the application of the Soft Selection Method [14] based on the Sum Squared Error evaluation criterion. Figure 12 depicts temperature \(t_1\) of the furnace and the adaptive thresholds obtained with (5) for the validation data set (no fault case).

Fig. 12
figure 12

Temperature \(t_1\) and output adaptive threshold

Full size image

At the next stage of the experiment the estimates of the GMDH neural model inputs and corresponding input adaptive thresholds are estimated with the application of the approach based on the application of the RUIF developed in Sect. 3. Figure 13 depicts the measurements of the input voltage \(u_1\) of the electric heater and the corresponding input adaptive threshold (37) obtained with the application of the RUIF for the fault free case.

Fig. 13
figure 13

Voltage \(u_1\) and input adaptive threshold

Full size image

After the synthesis of the GMDH model, it is used for the robust fault detection of the tunnel furnace. For this reason two faults were simulated. The first fault in the sensor was simulated by its partly removing from the tunnel furnace (for \(k\,=\,400\)). Figure 14 presents the measurements of temperature \(t_1\) from the faulty sensor and the output adaptive threshold obtained with the application of the state space GMDH neural model. As it can be seen the fault is detected for \(k\,=\,400\) when the measurements of temperature \(t_1\) crosses the output adaptive threshold calculated according Eq. (5).

Fig. 14
figure 14

Detection of faulty temperature sensor via output adaptive threshold

Full size image

The second fault was simulated in the actuator by the decreasing of the input voltage by \(20\,\%\). Figure 15 presents the measurements of the input voltage \(u_1\) and the corresponding input adaptive threshold obtained with the application of the GMDH neural model and RUIF. As it can be seen the faulty actuator is detected for \(k = 300\) when the value of voltage \(u_1\) crosses the input adaptive threshold (37).

Fig. 15
figure 15

Detection of faulty electric heater via input adaptive threshold

Full size image

5 Conclusion

The main objective of this paper was to develop a novel robust FDI method on the basis of the state-space GMDH neural model. The application of the UKF enables to obtain the asymptotically stable GMDH neural model during the network synthesis. Moreover, the application of such an algorithm enables to calculate the output adaptive threshold, which can be applied for the robust sensor fault detection. Furthermore, in the paper a novel methodology of the GMDH model input estimation with RUIF approach was proposed. Such an approach allows to calculate the input adaptive thresholds and enables to perform the robust fault detection and isolation of the actuators. The main contribution of this paper is to show how to use the above-mentioned results for actuator fault diagnosis. In particular, the obtained neural model is used to obtain the input estimates. The achieved estimates are then compared with the original input signals to formulate the diagnostics decisions. The input estimation scheme is based on a chain of robust observers, which guaranties that the input estimates are obtained with a prescribed disturbance attenuation level while ensuring the convergence of the observers. The final part of the paper shows a comprehensive case study regarding the laboratory tunnel furnace, which exhibit the performance of the proposed approach.