Fault Isolation and Estimation in Networks of Linear Process Systems

Abstract

Fault detection and isolation is a ubiquitous task in current complex systems even in the linear networked case when the complexity is mainly caused by the complex network structure. A simple yet practically important special case of networked linear process systems is considered in this paper with only a single conserved extensive quantity but with a network structure containing loops. These loops make fault detection and isolation challenging to perform because the effect of fault is propagated back to where it first occurred. As a dynamic model of network elements, a two input single output (2ISO) LTI state-space model is proposed for fault detection and isolation where the fault enters as an additive linear term into the equations. No simultaneously occurring faults are considered. A steady state analysis and superposition principle are used to analyse the effect of faults in a subsystem that propagates to the sensors’ measurements at different positions. This analysis is the basis of our fault detection and isolation procedure that provides the position of the faulty element in a given loop of the network. A disturbance observer is also proposed to estimate the magnitude of the fault inspired by a proportional-integral (PI) observer. The proposed fault isolation and fault estimation methods have been verified and validated by using two simulation case studies in the MATLAB/Simulink environment.

1. Introduction

Nowadays, we live in a world that is surrounded by networks, e.g., computer networks, transportation networks, social networks, electrical networks, etc. Although they could comprise simple elements, their general large number and interconnections make them an important subclass of complex systems. As a dynamic system, a network presents many theoretical challenges during control or diagnosis method design. There have been many excellent surveys about this in systems and control literature (see, e.g., [1,2,3,4,5,6,7]).

Fault detection and isolation (FDI) is a subfield of control engineering which mainly deals with process monitoring to ensure the safety of production processes. The basic FDI notions are used similarly as presented in [8]. The fault is considered an unexpected variation of some process or environmental variable that could yield unacceptable changes in the process behaviour. During system model construction, additive faults are assumed and the faults are defined as fictive inputs. By fault detection, we mean a decision on the presence of any fault or the absence of all faults. Meanwhile, fault isolation is related to the localization of a fault. During weak fault isolation, the faults are assumed to happen one at a time. By fault estimation, we mean the method that can determine the steady state value of the additive faults based on the available measurements.

Although it is crucial to ensure safe operation in networks of process systems, the number of works related to FDI in networked systems is still limited [9].

There exist many different techniques that can be applied to solve an FDI problem. The important main categories comprise model based and data-driven methods [10,11,12]. Commonly, model-based methods use dynamic models derived from first principles that rely heavily on the knowledge of basic physical and/or chemical processes. Meanwhile, the current popular data-driven methods rely on data from system operations that are used to train machine learning (ML) or artificial intelligence (AI) data structures to detect and identify faults. However, even though data-driven methods could be easier to implement compared to model-based approaches, their reliability strictly depends on the availability and quality of the data [13].

The FDI problem in industrial networks, such as process networks, presents special challenges. For the sake of energy and/or material efficiency, such networks usually contain loops. Therefore, fault diagnosis is harder to perform as the effect of the fault is propagated back to the subsystem where it first occurred [14,15,16,17,18]. Recent research has approached the fault isolation problem in this loop using an improved deep neural network [19].

Another difficulty in the complex networked industrial systems related to FDI is the sparsity of sensors. Because of financial reasons, it is rare to install a multitude of sensors along the connections of the interconnected subsystems in the network. Commonly, they are placed at the end of a connection or at some of the subsystem’s outputs. This motivates the development of an FDI method that groups network elements for fault identification [20,21].

As mentioned before, the interconnections that define the network topology make a networked dynamic system complex even if one considers simple elements (subsystems) in the network. Therefore, a practically important subclass of process systems is considered in this work where the elements are linear process systems connected by linear (possibly dynamic) connections.

Despite being a special case, it is a very important subclass because it can describe the dynamic behaviour of processes that serve our daily basic needs, e.g., domestic heating/cooling systems. Even some cellular processes belong to this subclass [22]. The fault that occurred in one of the subsystems is considered to be generated by an external source that can be treated as an extra input in a subsystem. Physically, this kind of fault represents some leakage phenomenon common inside networked linear process systems [23].

A model-based approach is proposed in this research to handle the fault diagnosis problem in a network of linear process systems, which may contain cycles/loops. The proposed approach does not need high computational costs. The burden of high computational power is common for fault diagnosis in a network [24].

The research questions that we are going to tackle are:

• In a process network with loops, how can a fault be identified in a subsystem or group of subsystems regardless of the fault effect propagation through the network loops and branches?
• Which sensor measurements are necessary to isolate a fault in a subsystem of a process network (or group of subsystems), and how can measurement noise be handled during fault diagnosis?
• If fault isolation can be conducted, how can the magnitude of the fault be estimated?

In concordance with the addressed research questions, the main contributions of our proposed approach can be summarized as follows:

• A fault diagnosis-oriented modelling approach is developed for such process networks in which the transport mechanisms can be described by stable and positive linear dynamic systems. The model includes the sources of the fault and indicates the sensor measurements. Based on this model, an analysis of how the fault input affects the steady states of the subsystems in the network is presented.
• A model-based fault isolation approach is proposed for networks of linear process systems. The main benefit of this approach compared to previously reported ones is that, under reasonable assumptions, it can isolate the faults in a subsystem or group of subsystems regardless of the loops and ramifications in the network. The proposed algorithm also indicates which sensor measurements are necessary to perform the fault isolation tasks. Moreover, it is also applicable in the presence of measurement noise.
• Finally, to determine the magnitude of the localized fault, the design of a disturbance observer-based asymptotic fault estimator is proposed.

2. Networked Linear Process Systems

Linear process systems are systems where there are only linear source terms, such as chemical reactions and phase changes (evaporation, condensation etc.), present besides the usual transport terms that include convection, diffusion and transfer through phase boundaries [22]. Although this is a special case, linear process systems are of great practical importance, e.g., heat exchanger networks and domestic heating/cooling systems.

Most linear process systems have a complex structure. They comprise linear subsystems connected by a linear network that can provide static or dynamic connections between the subsystems. Such composite linear process systems are called networked linear process systems. This section describes the dynamic modelling of networked linear process systems where we consider that there are $j=1\dots N$ subsystems.

2.1. Linear Process Subsystems

For fault isolation and estimation, one usually applies simplified dynamic models where the effect of the considered fault(s) is also described. These dynamic models are constructed from first engineering principles (see [22]) where the state equations originate from dynamic balance equations for the conserved extensive quantities (such as energy, overall mass, component masses) in a balanced volume. However, they will be transformed to their intensive variable forms using algebraic constitutive equations. Therefore, the state variables in linear process subsystems are considered to be intensive variables, e.g., temperature, concentration, etc.

Intensive variables usually spatially vary within a balanced volume so that the conservation balance equations related to them are in the form of partial differential equations (PDEs). To obtain simple dynamic models suitable for fault isolation and estimation, the lumped form of process models is used to approximate PDEs to bring them into a set of ordinary differential equations (ODEs) [22].

It is important to note that the dynamic linear interconnections in networked process systems can be realized by PDEs to model linear process subsystems where the lumping produces a finite-dimensional model representing the distributed delay phenomenon in process systems [25].

2.1.1. Modelling Assumptions

The following modelling assumptions are used in this paper:

A1

Only the transport of a single conserved extensive quantity (such as component mass, or energy) is considered in the process systems. Thus, we have either energy-transport or mass-transport systems. Heat exchanger networks and domestic heating/cooling systems belong to the linear energy-transport class. More details can be found in [26,27].

A2

Only linear convection and transfer is considered without any linear source.

A3

Constant overall mass and constant physicochemical parameters (such as density, specific heat, heat transfer coefficient, and convective flow rate) are assumed.

A4

One inlet and one outlet flow are considered where the inputs of the systems are the intensive variable (temperature or concentration) at the inlet and that of the balance volume with which transfer is assumed. Meanwhile, the output is the intensive variable (temperature or concentration) at the outlet.

2.1.2. Model Equations

Based on the above general assumptions and considering spatially homogeneous lumping (i.e., the same parameters in every lump), the dynamic model equations of a linear process subsystem becomes a two input–single output (2ISO) LTI state-space model which has the form:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S^{\left(j\right)}:=& \left\{\begin{array}{c}{\dot{\mathit{x}}}^{\left(j\right)}=A^{\left(j\right)}{\mathit{x}}^{\left(j\right)}+B^{\left(j\right)}{\mathit{u}}_X^{\left(j\right)}\hfill \\ y^{\left(j\right)}=C{\mathit{x}}^{\left(j\right)}\hfill \end{array}\right.\hfill \\ \hfill A^{\left(j\right)}=& \left[\begin{array}{cccccc}-(v^{\left(j\right)}+k_E^{\left(j\right)})& 0& 0& \dots & \dots & 0\\ v^{\left(j\right)}& -(v^{\left(j\right)}+k_E^{\left(j\right)})& 0& \dots & \dots & 0\\ 0& v^{\left(j\right)}& -(v^{\left(j\right)}+k_E^{\left(j\right)})& \dots & \dots & 0\\ 0& \dots & \dots & \dots & \dots & 0\\ 0& \dots & \dots & \dots & v^{\left(j\right)}& -(v^{\left(j\right)}+k_E^{\left(j\right)})\end{array}\right]\hfill \\ \hfill B^{\left(j\right)}=& {\left[\begin{array}{cccc}{v}^{\left(j\right)}& 0& \dots & 0\\ k_E^{\left(j\right)}& k_E^{\left(j\right)}& \dots & k_E^{\left(j\right)}\end{array}\right]}^T,C=\left[\begin{array}{ccccc}0& 0& 0& \dots & 1\end{array}\right]\hfill \end{array}\end{array}$$

(1)

where $S^{\left(j\right)}$ is the general model of the jth subsystem with ${\mathit{x}}^{\left(j\right)}$ being the state variable, $y^{\left(j\right)}$ is the output variable, and ${\mathit{u}}_X^{\left(j\right)}={\left[u^{\left(j\right)}{u}_E^{\left(j\right)}\right]}^T$ is the input variable which consists of $u^{\left(j\right)}$ as the intensive variable of the jth subsystem at the inlet and $u_E^{\left(j\right)}$ as the intensive variable of the external balance volume (environment). Meanwhile, $v^{\left(j\right)}>0$ is the mass flow rate and $k_E^{\left(j\right)}>0$ is the transfer coefficient.

By looking at Equation (1), we can derive some important points:

• As long as $k_E^{\left(j\right)}$ is positive, the model will always be stable because the eigenvalues are already stable shown by the negative sign of the diagonal entries of the $A^{\left(j\right)}$ matrix (see Gershgorin circle theorem for details [28]).
• The $A^{\left(j\right)}\in {\mathbb{R}}^{n\times n}$ matrix is a Metzler matrix by which all of its off-diagonal elements are non-negative. Thus, the matrix represents the time delayed differential equations and positive linear dynamical systems [29]. This is understandable because between the input and output states of the intensive variables in the subsystem, a propagation inside the transport element will occur contributing to the increase in delayed processes. This also implies the stability of the model because this matrix is sign stable and Hurwitz [30].
• It is jointly controllable and observable because it is derived from conservation balance equations (it is a compartmental system) [31,32].

2.2. Interconnections and Topology of the Network

In the following, we consider process networks consisting of linear process subsystems described in Section 2.1 that are connected by static interconnections. In the case of dynamic interconnections, a special linear subsystem is used to represent the dynamics (i.e., the distributed delay) of the connection.

2.2.1. Physically Meaningful Connections

Considering the modelling assumptions and the model equations of the applied linear process subsystem models, we can make physically meaningful connections between two subsystems (the jth and ℓth, for example) by connecting part of the output flow of one (the jth) to the inlet flow of the other (the ℓth). This implies that the characterizing intensive variable $y^{\left(j\right)}$ will determine the intensive variable $u^{(\ell )}$ and the inlet of the ℓth subsystem.

2.2.2. Equations Describing the Connections

In a realistic network, it is common to have a branching phenomenon in the interconnection between subsystems. Here, the Kirchoff law applies to the flows of extensive variables entering and exiting this interconnection. Considering the ℓth junction, the sum of either the overall mass, component mass, or energy flows entering the junction ($v^{(\ell )}$) is equal to the sum of flows exiting ($v_{OUT}^{(\ell )}$), as shown in the following equation for the case of overall mass flows:

$$\sum _{k\in I{N}^{(\ell )}}{v}^{\left(k\right)}=\sum _{l\in OU{T}^{(\ell )}}{v}^{\left(l\right)}=V^{(\ell )}$$

(2)

where k runs over the branches by which the flow enters the interconnection. Meanwhile, l runs over those where the flows exits the interconnection. $I{N}^{(\ell )}$ is the input set of the interconnection junction point, $OU{T}^{(\ell )}$ is the output set of the interconnection junction point, and v is measured in $\left[\frac{k𝘨}{s}\right]$.

Similar conservation equations apply for the energy flows $\dot{Q}$ in energy-transport systems or the component mass $m_X$ in mass-transport systems. However, to obtain the relations between the intensive variables (temperature T or concentration $c_X$, respectively) among the subsystems in this networked process systems, we substitute the algebraic equations representing the relationship of the intensive–extensive variables into the mass conservation results. These algebraic relationships are in the following general form:

$$Q=M{c}_{P}T,m_X=M{c}_X$$

(3)

where M is the overall mass and $c_P$ is the specific heat in the balance volume.

Then, we can substitute relations in (3) into the Kirchhoff law in Equation (2) while taking into account that the value of the intensive variable for all the outflows is the same. This way, one obtains the linear algebraic equations for the intensive variables T or $c_X$ at the ℓth junction in the following general form:

$$\begin{array}{c}\hfill \begin{array}{ccc}& & \hfill \sum _{k\in I{N}^{(\ell )}}\frac{v^{\left(k\right)}}{V^{(\ell )}}{T}^{\left(k\right)}=T^{\left(l\right)}\mathrm{or}\sum _{k\in I{N}^{(\ell )}}\frac{v^{\left(k\right)}}{V^{(\ell )}}{c}_X^{\left(k\right)}=c_X^{\left(l\right)}\\ & & \hfill \sum _{k\in I{N}^{(\ell )}}\frac{v^{\left(k\right)}}{V^{(\ell )}}{y}^{\left(k\right)}=u^{\left(l\right)},\forall l\in OU{T}^{(\ell )}\end{array}\end{array}$$

(4)

2.2.3. Network Topology

An underlying graph ($\mathcal{G}$) can be associated with the process network. The edges of the graph represent the dynamic subsystems of the network. Meanwhile, the vertices represent the interconnections among the subsystems which have no dynamics.

Example 1. 

Consider a network composed of N subsystems. Figure 1 gives an example for this network of transport elements where $N=6$ and $v^{\left(j\right)}$ is the mass flow rate of the jth subsystem.

Based on Equation (2), we can derive some relations between $v^{\left(j\right)}$ in Figure 1 as follows: 

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 1.& v^{\left(A\right)}=v^{\left(B\right)}+v^{\left(C\right)}\hfill \\ \hfill 2.& v^{\left(B\right)}+v^{\left(C\right)}=v^{\left(D\right)}+v^{\left(E\right)}=V^{\left(1\right)}\hfill \\ \hfill 3.& v^{\left(D\right)}+v^{\left(E\right)}=v^{\left(F\right)}\hfill \\ \hfill 4.& v^{\left(F\right)}=v^{\left(A\right)}\hfill \end{array}\end{array}$$

(5)

If we want to obtain some relations in terms of temperature (either $u^{\left(j\right)}$ or $y^{\left(j\right)}$) among the subsystems from the same Figure 1, we can substitute Equation (3) into each relation in Equation (5). Thus, we obtain: 

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 1.& {\dot{Q}}^{\left(A\right)}={\dot{Q}}^{\left(B\right)}+{\dot{Q}}^{\left(C\right)},sothat{y}^{\left(A\right)}=u^{\left(B\right)}=u^{\left(C\right)}\hfill \\ \hfill 2.& {\dot{Q}}^{\left(B\right)}+{\dot{Q}}^{\left(C\right)}={\dot{Q}}^{\left(D\right)}+{\dot{Q}}^{\left(E\right)},sothat\frac{v^{\left(B\right)}}{V^{\left(1\right)}}{y}^{\left(B\right)}+\frac{v^{\left(C\right)}}{V^{\left(1\right)}}{y}^{\left(C\right)}=u^{\left(D\right)}=u^{\left(E\right)}\hfill \\ \hfill 3.& {\dot{Q}}^{\left(D\right)}+{\dot{Q}}^{\left(E\right)}={\dot{Q}}^{\left(F\right)},sothat\frac{v^{\left(D\right)}}{v^{\left(F\right)}}{y}^{\left(D\right)}+\frac{v^{\left(E\right)}}{v^{\left(F\right)}}{y}^{\left(E\right)}=u^{\left(F\right)}\hfill \\ \hfill 4.& {\dot{Q}}^{\left(F\right)}={\dot{Q}}^{\left(A\right)},sothat{y}^{\left(F\right)}=u^{\left(A\right)}\hfill \end{array}\end{array}$$

(6)

Note that in the considered process network model, all the edges are part of at least one loop in the graph.

3. Faults in the Network

3.1. Fault Modelling

In this work, the considered fault is a constant input signal which additively modifies the external intensive variable input signal $u_E^{\left(j\right)}$. The faulty external variable input signal $u_{Ef}^{\left(j\right)}$ has the form:

$$u_{Ef}^{\left(j\right)}=u_E^{\left(j\right)}+f^{\left(j\right)}$$

(7)

where $f^{\left(j\right)}$ is the fault signal in the jth subsystem.

Such fault modelling can describe several fault events: the unforeseen appearance of an unknown external source, or change in the heat transfer coefficient $k_E^{\left(j\right)}$.

If there is a change $f_k^{\left(j\right)}$ in the heat transfer coefficient, i.e., $k_{Ef}^{\left(j\right)}=k_E^{\left(j\right)}+f_k^{\left(j\right)}$, then the second input of the model (1) can be rephrased as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill (k_E^{\left(j\right)}+f_k^{\left(j\right)})u_E^{\left(j\right)}& =k_E^{\left(j\right)}(u_E^{\left(j\right)}+f^{\left(j\right)}),\hfill \\ \hfill \mathrm{where}{f}^{\left(j\right)}& =\frac{f_k^{\left(j\right)}{u}_E^{\left(j\right)}}{k_E^{\left(j\right)}}.\hfill \end{array}\end{array}$$

(8)

During the fault diagnosis algorithm design, we assumed that:

A5

The probability of multiple fault events happening at the same time in the network is negligible (weak fault isolation), i.e., there are no simultaneous faults during the fault isolation and estimation processes.

3.2. The Fault’s Effect in Different Measurement Positions in the Network

First, to investigate the fault’s effect on the subsystems of the network, the state space model in Equation (1) is converted into input–output realization as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y^{\left(j\right)}\left(s\right)& =S_1^{\left(j\right)}\left(s\right)u^{\left(j\right)}\left(s\right)+S_2^{\left(j\right)}\left(s\right)u_E^{\left(j\right)}\left(s\right)\hfill \\ \hfill S_1^{\left(j\right)}\left(s\right)& ={\left(\frac{v^{\left(j\right)}}{s+v^{\left(j\right)}+k_E^{\left(j\right)}}\right)}^n\hfill \\ \hfill S_2^{\left(j\right)}\left(s\right)& ={\displaystyle \sum _{h=1}^n}\left[\frac{{\left(v^{\left(j\right)}\right)}^{(h-1)}}{{(s+v^{\left(j\right)}+k_E^{\left(j\right)})}^h}\right]k_E^{\left(j\right)}\hfill \end{array}\end{array}$$

(9)

where $S_1^{\left(j\right)}\left(s\right)$ is the transfer function in Laplace domain from $u^{\left(j\right)}\left(s\right)$ to $y^{\left(j\right)}\left(s\right)$ and $S_2^{\left(j\right)}\left(s\right)$ is the transfer function in Laplace domain from $u_E^{\left(j\right)}\left(s\right)$ to $y^{\left(j\right)}\left(s\right)$. Zero initial states are assumed.

Figure 2 shows the proposed realization for the fault effect analysis. For the sake of convenience, the notation $S_i^{\left(j\right)}$ will be used instead of $S_i^{\left(j\right)}\left(s\right)$ from here on.

Now, consider a loop in the process network as shown in Figure 3. There, $S^{\left(j\right)}$ represents the block diagram of the jth subsystem which contains $S_1^{\left(j\right)}$ and $S_2^{\left(j\right)}$ as shown in Figure 2. Meanwhile, $i^{\left(j\right)}$, where $j=1\dots m$, represents the inflows of the subsystems that are not part of the loop (possible joining connections), and $o^{\left(j\right)}$, where $j=1\dots m$, represents the outflows from the loop in splitting connections.

It is considered that the sensors are placed at the outputs $y^{\left(l\right)}$ and $y^{\left(m\right)}$ where $m\ge 2$ and $0<l<m$.

The fault $f^{\left(k\right)}$ is represented by a constant input in the kth subsystem, where $0<k<l$, which enters the subsystem from the same channel as the external source $u_E^{\left(k\right)}$.

For fault effect analysis, the final value theorem (FVT) of the Laplace transform is applied:

$$\underset{t\to \infty }{lim}f\left(t\right)=\underset{s\to 0}{lim}sF\left(s\right)$$

(10)

where $F\left(s\right)$ is the Laplace transform of $f\left(t\right)$.

To calculate the steady state value of the lth subsystem’s output in the fault-free case $f^{\left(k\right)}=0$ (see Figure 2), the FVT is applied to the previous transfer functions along with the superposition principle. By assuming that the inputs are step functions with zero initial conditions, it yields:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{ss}^{\left(l\right)}& ={\displaystyle \sum _{j=1}^l}\left[\frac{{\parallel S_2^{\left(j\right)}\parallel }_0{\displaystyle \prod _{h=j+1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)u_E^{\left(j\right)}\hfill \\ & +{\displaystyle {\displaystyle \sum _{j=l+1}^m}}\left[\frac{{\parallel S_2^{\left(j\right)}\parallel }_0{\displaystyle \prod _{h=j+1}^m}\parallel S_1^{\left(h\right)}{\parallel }_0{\displaystyle \prod _{h=1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)u_E^{\left(j\right)}\hfill \\ & +{\displaystyle \sum _{j=1}^l}\left[\frac{{\displaystyle \prod _{h=j}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)i_{ss}^{\left(j\right)}+{\displaystyle \sum _{j=l+1}^m}\left[\frac{{\displaystyle \prod _{h=j}^m}\parallel S_1^{\left(h\right)}{\parallel }_0{\displaystyle \prod _{h=1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)i_{ss}^{\left(j\right)}\hfill \end{array}\end{array}$$

(11)

where $V^{\left(j\right)}$ is the sum of the mass flow rate passing through the jth joining/splitting connection before the jth subsystem input, $y_{ss}^{\left(l\right)}$ is the steady state value of $y^{\left(l\right)}$ when there is no fault, $i_{ss}^{\left(j\right)}$ is the steady state value of $i^{\left(j\right)}$, and ${\parallel ·\parallel }_0$ is the steady state gain of the related transfer function. $u_E^{\left(j\right)}$ is assumed to be constant.

Since the addressed subsystem class is positive (see Section 2.1), the steady state gains are also always positive.

Note that the terms of $\left(1-\prod \parallel S_1^{\left(j\right)}{\parallel }_0\right)$ appear because of the loop. Meanwhile, the terms of $\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)$ come from the mass/energy conservation balance (see Section 2.2).

When a step-like fault arises in the kth subsystem ($f^{\left(k\right)}\ne 0$), we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{fss}^{\left(l\right)}& =\left[\frac{{\parallel S_2^{\left(k\right)}\parallel }_0{\displaystyle \prod _{j=k+1}^l}{\parallel S_1^{\left(j\right)}\parallel }_0}{1-{\displaystyle \prod _{j=1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}\right]\left(\frac{v^{\left(k\right)}}{V^{\left(k\right)}}\right)f^{\left(k\right)}\hfill \\ & +{\displaystyle \sum _{j=1}^l}\left[\frac{{\parallel S_2^{\left(j\right)}\parallel }_0{\displaystyle \prod _{h=j+1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)u_E^{\left(j\right)}\hfill \\ & +{\displaystyle \sum _{j=l+1}^m}\left[\frac{{\parallel S_2^{\left(j\right)}\parallel }_0{\displaystyle \prod _{h=j+1}^m}\parallel S_1^{\left(h\right)}{\parallel }_0{\displaystyle \prod _{h=1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)u_E^{\left(j\right)}\hfill \\ & +{\displaystyle \sum _{j=1}^l}\left[\frac{{\displaystyle \prod _{h=j}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)i_{ss}^{\left(j\right)}+{\displaystyle \sum _{j=l+1}^m}\left[\frac{{\displaystyle \prod _{h=j}^m}\parallel S_1^{\left(h\right)}{\parallel }_0{\displaystyle \prod _{h=1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)i_{ss}^{\left(j\right)}\hfill \end{array}\end{array}$$

(12)

where $y_{fss}^{\left(j\right)}$ is the steady state value of $y^{\left(j\right)}$ in the presence of a fault.

Now, by subtracting Equation (11) from Equation (12), we obtain the deviation of the faulty output related to the fault-free case:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}& =\left[\frac{{\parallel S_2^{\left(k\right)}\parallel }_0{\displaystyle \prod _{j=k+1}^l}{\parallel S_1^{\left(j\right)}\parallel }_0}{1-{\displaystyle \prod _{j=1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}\right]\left(\frac{v^{\left(k\right)}}{V^{\left(k\right)}}\right)f^{\left(k\right)}\hfill \end{array}\end{array}$$

(13)

The difference between the faulty and fault-free outputs of the mth subsystem can be computed similarly:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}& =\left[\frac{{\parallel S_2^{\left(k\right)}\parallel }_0{\displaystyle \prod _{j=k+1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}{1-{\displaystyle \prod _{j=1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}\right]\left(\frac{v^{\left(k\right)}}{V^{\left(k\right)}}\right)f^{\left(k\right)}\hfill \end{array}\end{array}$$

(14)

Equations (13) and (14) show that the fault influences all subsystems in the loop. However, the fault effect on the outputs of the subsystems is calculable.

Example 2. 

Figure 4 presents two connected subsystems which form a loop. The input–output models of the subsystems are $S^{\left(P\right)}$ and $S^{\left(Q\right)}$, respectively, where each contains $S_1^{\left(P\right)}$, $S_2^{\left(P\right)}$, $S_1^{\left(Q\right)}$, and $S_2^{\left(Q\right)}$ as shown in Figure 2.

No joining and splitting connections are assumed ($i^{\left(j\right)}=o^{\left(j\right)}=0$), so $v^{\left(P\right)}=v^{\left(Q\right)}=V$ which leads to $\frac{v^{\left(j\right)}}{V^{\left(j\right)}}=1$ where $j=P,Q$. Moreover, $y^{\left(P\right)}=u^{\left(Q\right)}$ and $y^{\left(Q\right)}=u^{\left(P\right)}$.

From Equations (13) and (14), the step-like fault $f^{\left(P\right)}\ne 0$ generates deviations in the subsystems as follows: 

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{fss}^{\left(P\right)}-y_{ss}^{\left(P\right)}& =\left(\frac{\parallel S_2^{\left(P\right)}{\parallel }_0}{1-\parallel S_1^{\left(P\right)}{\parallel }_0{\parallel S_1^{\left(Q\right)}\parallel }_0}\right)f^{\left(P\right)}\hfill \\ \hfill y_{fss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}& =\left(\frac{\parallel S_2^{\left(P\right)}{\parallel }_0{\parallel S_1^{\left(Q\right)}\parallel }_0}{1-\parallel S_1^{\left(P\right)}{\parallel }_0{\parallel S_1^{\left(Q\right)}\parallel }_0}\right)f^{\left(P\right)}\hfill \end{array}\end{array}$$

(15)

4. Fault Diagnosis

4.1. Problem Formulation

Given a process network as defined in Section 2.2, consider that a fault event arises in one of the subsystem in the network. As a consequence of the branches and loops in the underlying graph ($\mathcal{G}$) of the process network, the effect of this fault could induce deviations in the outputs of all the subsystems. Moreover, due to the loops, the fault could propagate back to the input of the faulty subsystem as well. The loops and network branches make the fault effect propagation barely traceable.

The presence of a fault in the network can simply be detected by comparing the measured outputs of the subsystems in the network with the predicted outputs based on a reliable model of the corresponding subsystem. However, due to fault effect propagation, the localization of the fault source in the network is a difficult task.

Recall that the subsystems represent the edges in the underlying graph $\mathcal{G}$ of the process network. We consider a directed path of subsystems $\left\{S^{\left(1\right)}\dots S^{\left(l\right)}\right\}$ where $l>0$ such that this path may be part of at least one simple loop of $\mathcal{G}$.

Let a fault event happen in a subsystem in the network. Formulate the following fault diagnosis problems:

• Consider a path of l subsystems that can be part of a loop consisting of $m>l$ subsystems (see Figure 3). We must determine whether the fault occurred in the addressed path. Furthermore, we must determine which measurements are necessary to perform the isolation problem (sensors placement) in this path.
• If the fault has been isolated in one subsystem, an estimation algorithm must be designed that outputs the magnitude of the fault.

Note that the path can be part of more than one loop in the graph of the process network, e.g., in Figure 1, the path $\left\{S^{\left(A\right)},S^{\left(F\right)}\right\}$ is in more than one simple loop. For an algorithm that finds all the simple loops in a graph, see, e.g., [33]. In this case, such a loop should be chosen that is more representative for the fault isolation process, e.g., from the perspective of the sensor placement.

The formulated diagnosis problem can specify whether the fault happened in a group of subsystems (or one subsystem for $l=1$), or in some other part of the network. However, the diagnosis process can be repeated for different groups (or subsystems). Thus, by exhaustive search, the fault in the network can be localized.

The fault isolation algorithm is designed by applying the fault effect analysis as presented in Section 3.2. For fault estimator design, we take the PI observer approach.

4.2. Fault Isolation

To derive a fault isolation algorithm, we consider that two sensors are placed in a loop at two different locations (see Figure 3). They measure the output of the lth and mth subsystems.

The isolation logic is based on Equations (13) and (14). By subtracting these equations, we obtain:

$$(y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)})-(y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)})=\left[\left(\frac{{\parallel S_2^{\left(k\right)}\parallel }_0{\displaystyle \prod _{j=k+1}^l}{\parallel S_1^{\left(j\right)}\parallel }_0}{1-{\displaystyle \prod _{j=1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}\right)\left(1-{\displaystyle \prod _{j=l+1}^m}\parallel S_1^{\left(j\right)}{)\parallel }_0\right)\right]\left(\frac{v^{\left(k\right)}}{V^{\left(k\right)}}\right)f^{\left(k\right)}$$

(16)

We make the following assumptions on the steady state gains of the subsystems:

A6

Either every ${\parallel S_1^{\left(j\right)}\parallel }_0\in (0,1)$ or every ${\parallel S_1^{\left(j\right)}\parallel }_0>1$;

A7

${\parallel S_2^{\left(j\right)}\parallel }_0>0\forall j$.

In the view of Equation (16), if these assumptions hold, then

$$\left(\frac{{\parallel S_2^{\left(k\right)}\parallel }_0{\displaystyle \prod _{j=k+1}^l}{\parallel S_1^{\left(j\right)}\parallel }_0}{1-{\displaystyle \prod _{j=1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}\right)\left(1-{\displaystyle \prod _{j=l+1}^m}\parallel S_1^{\left(j\right)}{)\parallel }_0\right)>0$$

(17)

Note that A6 and A7 hold true for the subsystem’s model introduced in Section 2.1; if we apply the FVT to Equation (9), we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}& {\parallel S_1^{\left(j\right)}\parallel }_0=\underset{s\to 0}{lim}{S}_1^{\left(j\right)}\left(s\right)={\left(\frac{v^{\left(j\right)}}{v^{\left(j\right)}+k_E^{\left(j\right)}}\right)}^n\in (0,1)\hfill \\ & {\parallel S_2^{\left(j\right)}\parallel }_0=\underset{s\to 0}{lim}{S}_2^{\left(j\right)}\left(s\right)={\displaystyle \sum _{h=1}^n}\left[\frac{{\left(v^{\left(j\right)}\right)}^{(h-1)}}{{(v^{\left(j\right)}+k_E^{\left(j\right)})}^h}\right]k_E^{\left(j\right)}>0\hfill \end{array}\end{array}$$

(18)

Furthermore, by assumptions A6 and A7, Equation (16) leads to:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill (y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)})>(y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)})& \mathrm{f}\mathrm{o}\mathrm{r}{f}^{\left(k\right)}>0\mathrm{a}\mathrm{n}\mathrm{d}\hfill \\ \hfill (y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)})<(y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)})& \mathrm{f}\mathrm{o}\mathrm{r}{f}^{\left(k\right)}<0\hfill \\ \hfill \therefore \parallel y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel >0& \mathrm{f}\mathrm{o}\mathrm{r}{f}^{\left(k\right)}\ne 0\hfill \end{array}\end{array}$$

(19)

where $\parallel ·\parallel$ is the absolute value of the related function.

Thus, Equation (19) shows that $\parallel y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel >0$ when a fault occurs in a subsystem between $S^{\left(1\right)}$ and $S^{\left(l\right)}$.

With the same assumptions, we can perform the same derivation to obtain $\parallel y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel <0$ when a fault occurs in a subsystem between $S^{(l+1)}$ and $S^{\left(m\right)}$.

For implementation of the fault isolation algorithm, the fault-free steady state value $y_{ss}^{\left(l\right)}$ and $y_{ss}^{\left(m\right)}$ have to be known prior, or they have to be computed. In view of the relation in (11), to compute $y_{ss}^{\left(l\right)}$ and $y_{ss}^{\left(m\right)}$, the steady state value of the inputs $i_{ss}^{\left(j\right)}$ and $u_E^{\left(j\right)}$ have to be measured.

To conclude, the fault isolation can be performed according to Algorithm 1 as follows:

Algorithm 1 Fault isolation algorithm.

Measure $y_{fss}^{\left(l\right)}$ and $y_{fss}^{\left(m\right)}$ in steady state. Compute $y_{ss}^{\left(l\right)}$ and $y_{ss}^{\left(m\right)}$. Isolate the fault: - If $\parallel y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel =\parallel y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel =0$, then no fault event occurred. - If $\parallel y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel >0$, then the fault occurred before l and after m. - If $\parallel y_{fss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel <0$, then the fault occurred before m and after l.

4.3. Fault Isolation in the Presence of Measurement Noise

In a realistic environment, it has to be considered that the measurements on the subsystems are affected by signal noise. In model (1), the noise inputs that influence the model’s inputs and outputs are introduced as:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{\mathit{x}}}^{\left(j\right)}=A^{\left(j\right)}{\mathit{x}}^{\left(j\right)}+B^{\left(j\right)}\left({\mathit{u}}_X^{\left(j\right)}+{\mathit{w}}_{uX}^{\left(j\right)}\right)\hfill \\ y^{\left(j\right)}=C{\mathit{x}}^{\left(j\right)}+w_y^{\left(j\right)}\hfill \end{array}\right.\end{array}$$

(20)

where ${\mathit{w}}_{uX}^{\left(j\right)}={\left[w_u^{\left(j\right)}{w}_E^{\left(j\right)}\right]}^T$.

The following assumptions are considered for signal noise:

A8

${\parallel w_u^{\left(j\right)}\left(t\right)\parallel }_{\infty }\le w_{uM}^{\left(j\right)},\parallel w_E^{\left(j\right)}{\left(t\right)\parallel }_{\infty }\le w_{EM}^{\left(j\right)},{\parallel w_y^{\left(j\right)}\left(t\right)\parallel }_{\infty }\le w_{yM}^{\left(j\right)}$

where ${\parallel ·\parallel }_{\infty }$ is the infinity norm of the related function.

A9

$w_u^{\left(j\right)}\left(t\right),w_E^{\left(j\right)}\left(t\right),w_y^{\left(j\right)}\left(t\right)$ are locally integrable.

A10

The signal noise is unbiased with known variances:

$w_E^{\left(j\right)}\left(t\right)\sim (0,R_u^{\left(j\right)}),w_y^{\left(j\right)}\left(t\right)\sim (0,Q_y^{\left(j\right)})$.

A11

$w_u^{\left(j\right)}\left(t\right),w_E^{\left(j\right)}\left(t\right),w_y^{\left(j\right)}\left(t\right)$ are mutually uncorrelated with both each other and system states.

The input–output model in (9) in the presence of newly considered signal noise takes the form:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_n^{\left(j\right)}\left(s\right)& =S_1^{\left(j\right)}\left(s\right)\left[u^{\left(j\right)}\left(s\right)+w_u^{\left(j\right)}\left(s\right)\right]+S_2^{\left(j\right)}\left(s\right)\left[u_E^{\left(j\right)}\left(s\right)+w_E^{\left(j\right)}\left(s\right)\right]+w_y^{\left(j\right)}\left(s\right)\hfill \\ \hfill y_n^{\left(j\right)}\left(s\right)& =y^{\left(j\right)}\left(s\right)+\left[S_1^{\left(j\right)}\left(s\right)w_u^{\left(j\right)}\left(s\right)+S_2^{\left(j\right)}\left(s\right)w_E^{\left(j\right)}\left(s\right)+w_y^{\left(j\right)}\left(s\right)\right]\hfill \\ \hfill y_n^{\left(j\right)}\left(s\right)& =y^{\left(j\right)}\left(s\right)+w^{\left(j\right)}\left(s\right).\hfill \end{array}\end{array}$$

(21)

where $y_n^{\left(j\right)}$ is the subsystem’s output in the presence of noise and $w^{\left(j\right)}$ is the noise’s effect where

$${\parallel w^{\left(j\right)}\left(t\right)\parallel }_{\infty }\le w_M^{\left(j\right)}.$$

(22)

As the noise from each subsystem propagates inside the loop, the measurement at both points l and m in the network are affected. Thus, to compensate for noise in the fault isolation process, a suitable threshold needs to be specified.

By applying the FVT and superposition principle, the noise’s effect on the measurements at sensors l and m can be computed, obtaining:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill y_{fnss}^{\left(l\right)}& =y_{fss}^{\left(l\right)}+{\displaystyle \sum _{j=1}^l}\left[\frac{{\displaystyle \prod _{h=j+1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)w^{\left(j\right)}+{\displaystyle \sum _{j=l+1}^m}\left[\frac{{\displaystyle \prod _{h=j+1}^m}\parallel S_1^{\left(h\right)}{\parallel }_0{\displaystyle \prod _{h=1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)w^{\left(j\right)}\hfill \\ \hfill y_{fnss}^{\left(m\right)}& =y_{fss}^{\left(m\right)}+{\displaystyle \sum _{j=1}^{m-1}}\left[\frac{{\displaystyle \prod _{h=j+1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)w^{\left(j\right)}+\left[\frac{{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(m\right)}}{V^{\left(m\right)}}\right)w^{\left(m\right)}\hfill \end{array}\end{array}$$

(23)

where $y_{fnss}^{\left(j\right)}$ is the steady state value of $y^{\left(j\right)}$ in the presence of a fault and noise, $w^{\left(l\right)}$ is the noise’s effect on the measurement at sensor l, and $w^{\left(m\right)}$ is the noise’s effect on the measurement at sensor m.

With the addition of the noise’s effect, Equation (16) becomes:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill (y_{fnss}^{\left(l\right)}-y_{ss}^{\left(l\right)})-(y_{fnss}^{\left(m\right)}-y_{ss}^{\left(m\right)})& =\left[\left(\frac{{\parallel S_2^{\left(k\right)}\parallel }_0{\displaystyle \prod _{j=k+1}^l}{\parallel S_1^{\left(j\right)}\parallel }_0}{1-{\displaystyle \prod _{j=1}^m}{\parallel S_1^{\left(j\right)}\parallel }_0}\right)\left(1-{\displaystyle \prod _{j=l+1}^m}\parallel S_1^{\left(j\right)}{)\parallel }_0\right)\right]\left(\frac{v^{\left(k\right)}}{V^{\left(k\right)}}\right)f^{\left(k\right)}\hfill \\ & +{\displaystyle \sum _{j=1}^l}\underset{S_{wl}^{\left(j\right)}}{\underbrace{\left[\frac{{\displaystyle \prod _{h=j+1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0\left(1-{\displaystyle \prod _{h=l+1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0\right)}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)}}{w}^{\left(j\right)}\hfill \\ & +{\displaystyle \sum _{j=l+1}^m}\underset{S_{wm}^{\left(j\right)}}{\underbrace{\left[\frac{{\displaystyle \prod _{h=j+1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0\left({\displaystyle \prod _{h=1}^l}{\parallel S_1^{\left(h\right)}\parallel }_0-1\right)}{1-{\displaystyle \prod _{h=1}^m}{\parallel S_1^{\left(h\right)}\parallel }_0}\right]\left(\frac{v^{\left(j\right)}}{V^{\left(j\right)}}\right)}}{w}^{\left(j\right)}\hfill \end{array}\end{array}$$

(24)

Now, for measurements at points l and m, the threshold value $t{h}^{\left(lm\right)}$ is defined as:

$$t{h}^{\left(lm\right)}={\displaystyle \sum _{j=1}^l}\parallel S_{wl}^{\left(j\right)}{\parallel }_{\infty }{w}_M^{\left(j\right)}+{\displaystyle \sum _{j=l+1}^m}{\parallel S_{wm}^{\left(j\right)}\parallel }_{\infty }{w}_M^{\left(j\right)}$$

(25)

Then, by considering the measurement noise, Algorithm 1 is updated into Algorithm 2 as follows:

Algorithm 2 Fault isolation algorithm in the presence of measurement noise.

Measure $y_{fnss}^{\left(l\right)}$ and $y_{fnss}^{\left(m\right)}$ in steady state. Compute $y_{ss}^{\left(l\right)}$, $y_{ss}^{\left(m\right)}$, and $t{h}^{\left(lm\right)}$. Isolate the fault: - If $\parallel y_{fnss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fnss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel \in (-t{h}^{\left(lm\right)},t{h}^{\left(lm\right)})$, then no fault event occurred or the fault is negligible in comparison to the threshold. - If $\parallel y_{fnss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fnss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel >t{h}^{\left(lm\right)}$, then the fault occurred before l and after m. - If $\parallel y_{fnss}^{\left(l\right)}-y_{ss}^{\left(l\right)}\parallel -\parallel y_{fnss}^{\left(m\right)}-y_{ss}^{\left(m\right)}\parallel <-t{h}^{\left(lm\right)}$, then the fault occurred before m and after l.

Example 3. 

This example is an extension of Example 2.

Consider that in the loop shown in Figure 4, the subsystem’s sensors are affected by measurement noises $w^{\left(P\right)}$ and $w^{\left(Q\right)}$, where ${\parallel w^{\left(P\right)}\left(t\right)\parallel }_{\infty }<w_M^{\left(P\right)}$ and ${\parallel w^{\left(Q\right)}\left(t\right)\parallel }_{\infty }<w_M^{\left(Q\right)}$.

Then, using Equations (24) and (25), the threshold value can be computed as follows: 

$$t{h}^{\left(PQ\right)}={∥\frac{1-{\parallel S_1^{\left(Q\right)}\parallel }_0}{1-\parallel S_1^{\left(P\right)}{\parallel }_0{\parallel S_1^{\left(Q\right)}\parallel }_0}∥}_{\infty }{w}_M^{\left(P\right)}+{∥\frac{{\parallel S_1^{\left(P\right)}\parallel }_0-1}{1-\parallel S_1^{\left(P\right)}{\parallel }_0{\parallel S_1^{\left(Q\right)}\parallel }_0}∥}_{\infty }{w}_M^{\left(Q\right)}$$

(26)

4.4. Fault Estimation

After the fault has been isolated, under certain sensor placement assumptions, the fault estimation can be performed. A linear disturbance–observer approach is proposed to determine the magnitude of the fault given by Equation (7).

During the estimation process, consider $k=l=1$ and the fault has been isolated in subsystem $k=1$.

The state space model given by Equation (1) with fault (7) can be rewritten as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\mathit{x}}}^{\left(k\right)}& =A^{\left(k\right)}{\mathit{x}}^{\left(k\right)}+B^{\left(k\right)}{\mathit{u}}_X^{\left(k\right)}+E^{\left(k\right)}{f}^{\left(k\right)}\hfill \\ \hfill y^{\left(k\right)}& =C{\mathit{x}}^{\left(k\right)}\hfill \\ \hfill E^{\left(k\right)}& ={\left[\begin{array}{ccccccc}{k}_E^{\left(k\right)}& k_E^{\left(k\right)}& \dots & k_E^{\left(k\right)}& \dots & k_E^{\left(k\right)}& k_E^{\left(k\right)}\end{array}\right]}^T\hfill \end{array}\end{array}$$

(27)

where $E^{\left(k\right)}$ is a fault distribution column matrix and $f^{\left(k\right)}$ is assumed to be constant.

This new model in Equation (27) can be further transformed into an extended state space model by defining an extended state vector containing the fault as ${\mathit{z}}^{\left(\mathit{k}\right)}={\left[{\mathit{x}}^{\left(\mathit{k}\right)}{f}^{\left(k\right)}\right]}^T$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\mathit{z}}}^{\left(k\right)}& =\left[\begin{array}{cc}{A}^{\left(k\right)}& E^{\left(k\right)}\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathit{x}}^{\left(\mathit{k}\right)}\\ f^{\left(k\right)}\end{array}\right]+\left[\begin{array}{c}{B}^{\left(k\right)}\\ 0\end{array}\right]{\mathit{u}}_X^{\left(k\right)}=A_z^{\left(k\right)}{\mathit{z}}^{\left(k\right)}+B_z^{\left(k\right)}{\mathit{u}}_X^{\left(k\right)}\hfill \\ \hfill y^{\left(k\right)}& =\left[\begin{array}{cc}C& 0\end{array}\right]\left[\begin{array}{c}{\mathit{x}}^{\left(\mathit{k}\right)}\\ f^{\left(k\right)}\end{array}\right]=C_z{\mathit{z}}^{\left(k\right)}\hfill \end{array}\end{array}$$

(28)

As stated by [34], if the steady state value of the fault is not zero, a proportional observer cannot correctly estimate the states of the plant because there will always be a steady state error between the actual and estimated states. However, using a proportional-integral (PI) observer, the steady state error can be reduced.

The state space model of a PI observer is:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\widehat{\mathit{x}}}}^{\left(\mathit{k}\right)}& =A^{\left(k\right)}{\widehat{\mathit{x}}}^{\left(\mathit{k}\right)}+B^{\left(k\right)}{\mathit{u}}_X^{\left(k\right)}+L_P^{\left(k\right)}(y^{\left(k\right)}-{\widehat{y}}^{\left(k\right)})+E^{\left(k\right)}{\widehat{f}}^{\left(k\right)}\hfill \\ \hfill {\dot{\widehat{f}}}^{\left(k\right)}& =L_I^{\left(k\right)}(y^{\left(k\right)}-{\widehat{y}}^{\left(k\right)})\hfill \\ \hfill {\widehat{y}}^{\left(k\right)}& =C{\widehat{x}}^{\left(k\right)}\hfill \end{array}\end{array}$$

(29)

where ${\widehat{\mathit{x}}}^{\left(\mathit{k}\right)}$ is the estimated state vector, ${\widehat{f}}^{\left(k\right)}$ is the estimated fault magnitude, $L_P^{\left(k\right)}$ is the observer’s proportional gain, and $L_I^{\left(k\right)}$ is the observer’s integral gain.

Hence, a PI observer can be designed for the extended state space model in Equation (28) to not only estimate the states but also the fault. By applying Equation (29) into Equation (28) with the assumption that $(A_z^{\left(k\right)},C_z)$ is an observable pair, we obtain:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\widehat{\mathit{z}}}}^{\left(\mathit{k}\right)}& =\left(\left[\begin{array}{cc}{A}^{\left(k\right)}& E^{\left(k\right)}\\ 0& 0\end{array}\right]-\left[\begin{array}{c}{L}_P^{\left(k\right)}\\ L_I^{\left(k\right)}\end{array}\right]\left[\begin{array}{cc}C& 0\end{array}\right]\right){\widehat{\mathit{z}}}^{\left(\mathit{k}\right)}+\left[\begin{array}{c}{L}_P^{\left(k\right)}\\ L_I^{\left(k\right)}\end{array}\right]y^{\left(k\right)}+\left[\begin{array}{c}{B}^{\left(k\right)}\\ 0\end{array}\right]{\mathit{u}}_X^{\left(k\right)}\hfill \\ & =(A_z^{\left(k\right)}-L_z^{\left(k\right)}{C}_z){\widehat{\mathit{z}}}^{\left(\mathit{k}\right)}+L_z^{\left(k\right)}{y}^{\left(k\right)}+B_z^{\left(k\right)}{\mathit{u}}_X^{\left(k\right)}\hfill \\ \hfill {\widehat{y}}^{\left(k\right)}& =C_z{\widehat{\mathit{z}}}^{\left(\mathit{k}\right)}\hfill \end{array}\end{array}$$

(30)

Thus, if $L_P^{\left(k\right)}$ and $L_I^{\left(k\right)}$ are chosen such that $(A_z^{\left(k\right)}-L_z^{\left(k\right)}{C}_z)$ is Hurwitz, then ${lim}_{t\to \infty }({\mathit{z}}^{\left(k\right)}-{\widehat{\mathit{z}}}^{\left(k\right)})=0$.

Meanwhile, in the presence of noise, the widely known linear quadratic estimator (LQE) procedure can be used to compute the observer’s gain $L_z^{\left(k\right)}$. For details, see Appendix A.

5. Case Studies

To verify and validate the proposed fault isolation and fault estimation methods, two simulation case studies were examined in the MATLAB/Simulink environment. The first case study is based on Example 2. Meanwhile, for the second case study, the investigated network comprises six subsystems as presented in Example 1.

5.1. Case Study 1

In this case study, two identical subsystems $S^{\left(P\right)}$ and $S^{\left(Q\right)}$ are connected as shown in Figure 4. The model of the subsystems is given by Equation (20) where each subsystem is considered to have five state variables ($n=5$). The applied parameters and external inputs are shown in

Table 1. Case study 1—parameters and external inputs.

j	P	Q
$v^{\left(j\right)}$	10	10
$k_E^{\left(j\right)}$	3	3
$u_E^{\left(j\right)}$	400	400
$R_u^{\left(j\right)}$	20	20
$Q_y^{\left(j\right)}$	20	20

.

It is considered that the measurements $y^{\left(P\right)}$ and $y^{\left(Q\right)}$ are influenced by noise as presented in (20). To compensate for the noise that influences the system, a suitable threshold $t{h}^{\left(PQ\right)}$ is computed by applying Equation (25). From the proposed fault isolation logic (see Algorithm 2), as the loop is composed of two subsystems $S^{\left(P\right)}$ and $S^{\left(P\right)}$, the simulation results should show that $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel$ is above $t{h}^{\left(PQ\right)}$ when there is a fault in subsystem $S^{\left(P\right)}$, or $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel$ is below $-t{h}^{\left(PQ\right)}$ when there is a fault in subsystem $S^{\left(Q\right)}$.

In the first scenario, no fault is injected into the subsystems. The simulation results in Figure 5 show $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel \in (-t{h}^{\left(PQ\right)},t{h}^{\left(PQ\right)})$, i.e., no fault is injected.

In the second scenario, a constant fault signal $f^{\left(P\right)}=100$ is injected into subsystem $S^{\left(P\right)}$ at $t=20$ s. Figure 6 shows that the value of $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel$ is above the threshold after the fault event occurred. This indicates that a fault is occurring in subsystem $S^{\left(P\right)}$.

To verify the proposed fault estimation approach, an observer is designed specifically for subsystem $S^{\left(P\right)}$ based on Equation (30). By using the lqe MATLAB function (see Appendix A), the observer’s gain is computed such that the measurement noise variances are taken into consideration. Figure 7 shows the fault and its estimated value. It is seen that the fault’s magnitude is correctly estimated. Meanwhile, Figure 8 shows that the observer successfully estimates the states of subsystem $S^{\left(P\right)}$.

In the third scenario, a constant fault signal $f^{\left(Q\right)}=-50$ is injected into subsystem $S^{\left(Q\right)}$ at $t=20$ s. Figure 9 shows that the value of $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel$ is below $-t{h}^{\left(PQ\right)}$ after the fault event occurred. This indicates that a fault is occurring in subsystem $S^{\left(Q\right)}$, which is correct.

Moreover, an observer is designed specifically for subsystem $S^{\left(Q\right)}$ based on Equation (30). The observer’s gain is also computed using the same lqe MATLAB function such that the noise variances are taken into consideration. Figure 10 shows the fault and its estimated value. Here, the fault’s magnitude is also correctly estimated. Meanwhile, Figure 11 shows that the observer successfully estimates the states of subsystem $S^{\left(Q\right)}$.

5.2. Robustness Analysis (for Case Study 1)

To analyse the robustness of our proposed approach, modified parameter values are considered in the subsystem model (1) only during the simulations. Here, we have two parameters to modify: the mass flow rate and the transfer coefficient. The mass flow rate $v^{\left(j\right)}$ can be easily measured. However, the transfer coefficient $k_E^{\left(j\right)}$ can only be estimated as its value changes depending on the physical conditions and circumstances. Thus, the robustness analysis is performed by checking the performance of fault isolation and estimation in the presence of transfer coefficient parameter uncertainty.

By using case study 1 for the sake of convenience but without loss of generality, the simulation is carried out by increasing the value of $k_E^{\left(j\right)}$ in Table 1 by 25% ($k_E^{\left(j\right)}=3.75$ where $j=P,Q$). After this, a constant fault signal $f^{\left(Q\right)}=-50$ is injected into subsystem $S^{\left(Q\right)}$ at $t=20$ s. Figure 12 shows that the value of $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel$ is below $-t{h}^{\left(PQ\right)}$ after the fault event occurred. This indicates that a fault is occurring in subsystem $S^{\left(Q\right)}$, which is correct. The real and estimated fault values are shown in Figure 13. There, the fault’s magnitude is estimated with an error of less than 10%. This error is acceptable provided that the maximum tolerable value of the uncertainty in $k_E^{\left(j\right)}$ is 25%.

5.3. Incipient Fault Analysis (for Case Study 1)

We also analysed the performance of our proposed approach in the case of an incipient fault (slowly developing fault). To do this, by using case study 1 again, a linearly increasing fault signal is assumed with an initial value of 0 and a slope value of 0.1 per unit of time. This fault is injected into subsystem $S^{\left(P\right)}$ at $t=20$ s. Figure 14 shows that the value of $\parallel y_{fnss}^{\left(P\right)}-y_{ss}^{\left(P\right)}\parallel -\parallel y_{fnss}^{\left(Q\right)}-y_{ss}^{\left(Q\right)}\parallel$ is above the threshold $t{h}^{\left(PQ\right)}$ some time after the fault event occurred. This indicates that a fault is occurring in subsystem $S^{\left(P\right)}$, which is correct. Meanwhile, Figure 15 shows that the observer successfully estimates the magnitude of the incipient fault.

5.4. Case Study 2

In this case study, six subsystems $S^{\left(A\right)}$, $S^{\left(B\right)}$, $S^{\left(C\right)}$, $S^{\left(D\right)}$, $S^{\left(E\right)}$, and $S^{\left(F\right)}$ are connected as shown in Figure 1. The model of the subsystems is given by Equation (20) where each subsystem is considered to have five state variables ($n=5$). The applied parameters and external inputs are shown in

Table 2. Case study 2—parameters and external inputs.

j	A	B	C	D	E	F
$v^{\left(j\right)}$	20	4	16	8	12	20
$k_E^{\left(j\right)}$	3	3	3	3	3	3
$R_u^{\left(j\right)}$	20	20	20	20	20	20
$Q_y^{\left(j\right)}$	20	20	20	20	20	20

.

In this network, two sensors are placed to measure both $y^{\left(A\right)}$ and $y^{\left(F\right)}$ that are affected by noise as presented in (20). To compensate for the noise that influences the system, a suitable threshold $t{h}^{\left(AF\right)}$ is computed as presented in Equation (25). As both subsystem $S^{\left(A\right)}$ and $S^{\left(F\right)}$ are inside the loops of either $\left\{S^{\left(A\right)},S^{\left(B\right)},S^{\left(D\right)},S^{\left(F\right)}\right\}$ or $\left\{S^{\left(A\right)},S^{\left(B\right)},S^{\left(E\right)},S^{\left(F\right)}\right\}$, the simulation’s results should show that: $\parallel y_{fnss}^{\left(A\right)}-y_{ss}^{\left(A\right)}\parallel -\parallel y_{fnss}^{\left(F\right)}-y_{ss}^{\left(F\right)}\parallel$ is above the threshold $t{h}^{\left(AF\right)}$ when there is a fault in subsystem $S^{\left(A\right)}$, or $\parallel y_{fnss}^{\left(A\right)}-y_{ss}^{\left(A\right)}\parallel -\parallel y_{fnss}^{\left(F\right)}-y_{ss}^{\left(F\right)}\parallel$ is below $-t{h}^{\left(AF\right)}$ when there is a fault in either subsystem $S^{\left(B\right)}$, $S^{\left(C\right)}$, $S^{\left(D\right)}$, $S^{\left(E\right)}$, or $S^{\left(F\right)}$ (see Algorithm 2).

In the first scenario, no fault is injected into the subsystems. The simulation results in Figure 16 show that $\parallel y_{fnss}^{\left(A\right)}-y_{ss}^{\left(A\right)}\parallel -\parallel y_{fnss}^{\left(F\right)}-y_{ss}^{\left(F\right)}\parallel \in (-t{h}^{\left(AF\right)},t{h}^{\left(AF\right)})$, i.e., no fault is injected.

In the second scenario, a constant fault signal $f^{\left(A\right)}=200$ is injected into subsystem $S^{\left(A\right)}$ at $t=20$ s. Figure 17 shows that the value of $\parallel y_{fnss}^{\left(A\right)}-y_{ss}^{\left(A\right)}\parallel -\parallel y_{fnss}^{\left(F\right)}-y_{ss}^{\left(F\right)}\parallel$ is above the threshold after the fault event occurred. This indicates that a fault is occurring in subsystem $S^{\left(A\right)}$, which is correct.

To verify the proposed fault estimation approach, an observer is designed for subsystem $S^{\left(A\right)}$ based on Equation (30). By using the lqe MATLAB function, the observer’s gain is computed such that the measurement noise variances are taken into consideration. Figure 18 shows the fault and its estimated value. It is seen that the fault’s magnitude is correctly estimated. Meanwhile, Figure 19 shows that the observer successfully estimates the states of subsystem $S^{\left(A\right)}$.

In the third scenario, a constant fault signal $f^{\left(F\right)}=-150$ is injected into subsystem $S^{\left(F\right)}$ at $t=20$ s. Figure 20 shows that the value of $\parallel y_{fnss}^{\left(A\right)}-y_{ss}^{\left(A\right)}\parallel -\parallel y_{fnss}^{\left(F\right)}-y_{ss}^{\left(F\right)}\parallel$ is below $-t{h}^{\left(AF\right)}$ after the fault event occurred. This indicates that a fault is occurring either in subsystem $S^{\left(B\right)}$, $S^{\left(C\right)}$, $S^{\left(D\right)}$, $S^{\left(E\right)}$, or $S^{\left(F\right)}$, which is correct.

Moreover, an observer is designed specifically for subsystem $S^{\left(F\right)}$ based on the same Equation (30). The observer’s gain is computed using the same lqe MATLAB function such that the noise variances are taken into consideration. Figure 21 shows the fault and its estimated value. Here, the fault’s magnitude is also correctly estimated. Meanwhile, Figure 22 shows that the observer successfully estimates the states of subsystem $S^{\left(F\right)}$.

6. Conclusions

The case of networked linear process systems is considered in this paper with only a single conserved extensive quantity but with a network structure containing loops. It is assumed that the probability of multiple faults happening at the same time in the network is negligible.

For fault detection and isolation purposes, the network elements are described by a simple 2ISO LTI state space model where the fault enters as an additive linear term into the second input of the equations. Using the models of the network elements, a general model of the network is constructed which includes static splitting and joining connections. This results in an LTI state space model for the overall system that is suitable for fault detection and isolation.

By analysing the effect of the fault in a subsystem that propagates to the sensors’ measurements at different positions, a fault isolation algorithm (Algorithm 1) was proposed. It uses two sensors, installed at the output of the two subsystems placed along a path, to locate the fault (i.e., decide if the fault has occurred in a subsystem on the connecting path or outside it). A steady state analysis and superposition principle were used to build the algorithm. An improved version of the algorithm (Algorithm 2) was also proposed to perform the localization in the presence of measurement noise.

Having completed the fault isolation process, a PI-based disturbance observer was then proposed to estimate the magnitude of the fault.

Two simulation case studies were used in the MATLAB/Simulink environment to verify and validate the proposed fault isolation and fault estimation methods. The first case study comprised a loop with two subsystems only. Meanwhile, the investigated network of the second case study comprised six subsystems. Separate subsections on analysing the parametric robustness and effect of an incipient fault were provided using the first case study for the sake of convenience but without generality loss.

In both case studies, the simulation results show good performance both in fault detection and isolation, and in fault magnitude estimation. In addition, our proposed method also shows good robustness against the transfer coefficient $k_E$ which is difficult to estimate in practice. By using a ramp-like fault signal as an incipient fault, the simulation results show that our proposed fault estimation method can correctly estimate its magnitude.

Further work includes the extension of our method to cases of multiple faults, and verifying and validating the method on real data of, e.g., household heating/cooling systems.