Two-Step Neural-Network-Based Fault Isolation for Stochastic Systems

Abstract

This paper studies a fault isolation method for an optical fiber vibration source detection and early warning system. We regard the vibration sources in the system as faults and then detect and isolate the faults of the system based on a two-step neural network. Firstly, the square root B-spline expansion method is used to approximate the output probability density functions. Secondly, the nonlinear weight dynamic model is established through a dynamic neural network. Thirdly, the nonlinear filter and residual generator are constructed to estimate the weight, analyze the residual, and estimate the threshold, so as to detect, diagnose, and isolate the faults. The feasibility criterion of fault detection and isolation is given by using some linear matrix inequalities, and the stability of the estimation error system is proven according to the Lyapunov theorem. Finally, simulation experiments based on a optical fiber vibration source system are given to verify the effectiveness of this method.

1. Introduction

With the development of industrial automation, the scale and complexity of modern control systems are increasing rapidly. During operation, the sensor, actuator, and some components might fail due to aging, wearing, or other reasons [1,2,3,4]. If the faults of the system cannot be detected and eliminated in time in the process of operation, this may lead to system failure, paralysis, and even catastrophic consequences [5]. Therefore, in the past few decades, the theory of fault detection and isolation (FDI) has developed rapidly [6,7,8,9,10].

The purpose of fault detection and isolation is that, when multiple faults occur in the system, we can find and isolate these faults in time and estimate the size of each fault [6]. At present, many research achievements have been made in fault detection and fault isolation [11,12,13,14,15,16,17]. Classical fault detection and isolation methods include the T-S fuzzy model [11], the independent component analysis fault isolation method [12], the detection filter method [13], minimum variance fault estimation [14], etc.

In fact, actual industrial systems are stochastic systems [18,19,20,21]. Due to the nonlinearity of the system, even if the input obeys a Gaussian distribution, the output of the system will also obey a non-Gaussian distribution [18]. At present, there are many research results that focus on fault diagnosis and isolation for stochastic distribution systems. For example, in [8], an adaptive fault diagnosis observer was used to diagnose actuator faults and accurately estimate the fault size. In [19], fault diagnosis schemes for stochastic distributed systems were studied based on the minimum entropy principle. In [20], a fault isolation method based on entropy optimization filtering was studied for nonlinear non-Gaussian systems with unmeasurable outputs.

As we know, in order to monitor the operation of oil or gas pipelines, a commonly used method is to embed the oil or gas pipelines together with optical fibers in the same ditch. The conditions of the oil or gas pipelines can be monitored by observing the optical fiber’s output data in one end of the fiber if we input light at the other end [22,23,24]. However, In the optical fiber vibration source detection and early warning system, the sensors are quite sensitive due to inevitable non-stationary interferences from the environment, which often brings false alarms. To reduce the false alarm rate, we regard the vibration sources, such as pickaxe planing of mechanical excavation, as faults, and the objective is then transferred into detecting and isolating the faults. Different from the literature in the optical fiber vibration source detection field [22,23], we use a two-step neural network to model the optic data collected from the end of the fiber. The false alarm rate is controlled below the threshold, and the different vibration sources are estimated and isolated by using filter designing, convex optimization, Lyapunov theorems, etc.

This paper is organized as follows: In Section 2, the static modeling is carried out with the help of square root B-splines, and the nonlinear weighted dynamic model is established through the neural network to describe the FDI problem. In Section 3, based on the adaptive filter, the fault detection problem of the transformed nonlinear weighted dynamic model is studied. In Section 4, according to the adaptive fault diagnosis method, the size of different faults is estimated to achieve the purpose of fault isolation. The simulation is given in Section 5 to illustrate the feasibility of the results.

2. Static Modeling and Weighted Dynamic Modeling

For the optical fiber vibration source detection system, although the output value cannot be measured directly, the output PDFs can be estimated by the kernel density method [23]. Next, the output PDF at each sampling time is approximated by B-splines, and the weight of each B-splines is calculated.

2.1. Output PDFs’ Static Modeling

As in Figure 1, let $u\left(t\right)\in R^m$ be the input of the nonlinear stochastic optical fiber vibration detection system, $y\left(t\right)\in [a,b]$ be the output light intensity, $F_1$ and $F_2$ be the fault vectors (vibration sources) to be detected and diagnosed, and the probability distribution of the output $y\left(t\right)$ of the stochastic dynamic system in the interval $[a,\xi ]$ satisfy the condition $P(a\le y(t)\le \xi )$=${\int }_a^{\xi }\gamma (z,u\left(t\right),F_1,F_2)dz$. We can use the square root B-spline model with approximate error as shown in (1) to statically model the output PDFs [8]:

$$\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}=\sum _{i=1}^n{v}_i(u\left(t\right),F_1,F_2)b_i\left(z\right)+{\omega }_0(z,u\left(t\right),F_1,F_2)$$

(1)

where $v_i(u\left(t\right),F_1,F_2)(i=1,2,\cdots ,n)$ are the corresponding weights of B-spline expansion and $b_i\left(z\right)(i=1,2,\cdots ,n)$ are the pre-specified basis functions on interval $[a,b]$. For all $\{z,u\left(t\right),F_1,F_2\}$, the term ${\omega }_0(z,u\left(t\right),F_1,F_2)$ represents the model uncertainty or the error term on the approximation of the PDFs, which is supposed to satisfy $|{\omega }_0(z,u\left(t\right),F_1,$$F_2\left)\right|\le {\delta }_0$, where ${\delta }_0$ is assumed to be a known positive number. Denote

$$B\left(z\right)=[b_1\left(z\right),b_2\left(z\right),\cdots ,b_{n-1}\left(z\right)]$$

$$V\left(t\right)=[v_1(u\left(t\right),F_1,F_2),v_2(u\left(t\right),F_1,F_2),\cdots ,v_{n-1}(u\left(t\right),F_1,F_2)]$$

and

$${\Lambda }_1={\int }_a^b{B}^T\left(z\right)B\left(z\right)dz$$

(2)

$${\Lambda }_2={\int }_a^b{B}^T\left(z\right)b_n\left(z\right)dz$$

(3)

$${\Lambda }_3={\int }_a^b{\left(b_n\right)}^2\left(z\right)dz$$

(4)

where ${\Lambda }_1\in R^{(n-1)\times (n-1)}$, ${\Lambda }_2\in R^{(n-1)\times 1}$, and ${\Lambda }_3\in R^{1\times 1}$ are known matrices or constants. In the following, for simplicity, $V(u\left(t\right),F_1,F_2)$ is written as $V\left(t\right)$. According to the properties of the probability density function, for any $\gamma (z,u\left(t\right),F_1,F_2)$, the constraint condition ${\int }_a^{\xi }\gamma (z,u\left(t\right),F_1,F_2)dz=1$ is always true, and only $n-1$ weight vectors in the weight $\{v_i(u\left(t\right),F_1,F_2):i=1,2,\cdots ,n\}$ are independent. Let ${\Lambda }_0={\Lambda }_1{\Lambda }_3-{\Lambda }_2^T{\Lambda }_2$, then the following inequality holds [8]:

$$(1-{\omega }_1(z,u\left(t\right),F_1,F_2)){\Lambda }_3-V^T\left(t\right)V\left(t\right)\ge 0$$

(5)

where

$$\begin{array}{cc}\hfill {\omega }_1\left(z,u\left(t\right),F_1,F_2\right)& =2{\omega }_0^2\left(z,u\left(t\right),F_1,F_2\right)\times \left[\left({\int }_a^{b}B\left(z\right)dz\right)V\left(t\right)\right.\hfill \\ \hfill & \left.+\left({\int }_a^b{b}_n\left(z\right)dz\right)v_n\left(t\right)\right]+{\omega }_0^2\left(z,u\left(t\right),F_1,F_2\right)\left(b-a\right)\hfill \end{array}$$

From (5), we can obtain:

$$V^T\left(t\right){\Lambda }_{0}V\left(t\right)\le \left(1-{\omega }_1\left(z,u\left(t\right),F_1,F_2\right)\right){\Lambda }_3$$

(6)

where

$${\Lambda }_0>0,1-{\omega }_1\left(z,u\left(t\right),F_1,F_2\right)>0$$

(7)

According to (6), we can know $V\left(t\right)$ is bounded and ${∥V\left(t\right)∥}^2\le \tilde{\delta }=∥{\Lambda }_0^{-1}∥{\Lambda }_3$ holds. Therefore, (1) can be rewritten as:

$$\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}=B\left(z\right)V\left(t\right)+h_0(V\left(t\right),{\omega }_1)B_n\left(z\right)$$

(8)

where $h_0(V\left(t\right),{\omega }_1)$ is a function of $V\left(t\right)$ and ${\omega }_1(z,u\left(t\right),F_1,F_2)$, as shown in (9):

$$h_0\left(V\left(t\right),{\omega }_1\right)=\frac{1}{{\Lambda }_3}\sqrt{\left(1-{\omega }_1\right){\Lambda }_3-V^T\left(t\right){\Lambda }_{0}V\left(t\right)}-\frac{{\Lambda }_{2}V\left(t\right)}{{\Lambda }_3}$$

(9)

In order to simplify the B-spline model represented by (8), $h_0(V\left(t\right),{\omega }_1)$ can be further approximated as:

$$\begin{array}{c}\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}=B\left(z\right)V\left(t\right)+h\left(V\left(t\right)\right)b_n\left(z\right)+\omega (z,u\left(t\right),F_1,F_2)\end{array}$$

(10)

Equation (10) is a nonlinear output equation with uncertainty, and the term $h\left(V\right(t\left)\right)$ satisfies

$$h\left(V\left(t\right)\right)=\frac{1}{{\Lambda }_3}(-{\Lambda }_2\pm \sqrt{{\Lambda }_3-V^T\left(t\right){\Lambda }_{0}V\left(t\right)})$$

(11)

As $V\left(t\right)$ is bounded and $\left|{\omega }_0(z,u\left(t\right),F_1,F_2)\right|$$\le {\delta }_0$, it can be concluded that $\left|\omega (z,u\left(t\right),F_1,F_2)\right|\le \delta$ holds for all $\{z,u\left(t\right),F_1,F_2\}$. For $h\left(V\right(t\left)\right)$ in (11), it is supposed that, for any $V_1\left(t\right)$ and $V_2\left(t\right)$, there exists a known matrix $U_1$, such that:

$$∥h\left(V_1\left(t\right)\right)-h\left(V_2\left(t\right)\right)∥\le ∥U_1(V_1\left(t\right)-V_2\left(t\right))∥$$

(12)

2.2. Nonlinear Dynamic Weight Model

After the B-spline expansion of the output PDFs, the next task is to find the dynamic relationship between $u\left(t\right)$ and $\gamma$. As $V\left(t\right)$ is a nonlinear function of $u\left(t\right)$, we perform the modeling with the help of a neural network as in [25] and study the following weight dynamic system:

$$\left\{\begin{array}{c}\dot{x}\left(t\right)=Ax\left(t\right)+Gg\left(x\left(t\right)\right)+Hu\left(t\right)+F_1+F_2\hfill \\ V\left(t\right)=Dx\left(t\right)\hfill \end{array}\right.$$

(13)

where $x\left(t\right)\in R^m$ is the state vector, A, G, H, and D represent the known parameter matrices, and $g\left(x\right(t\left)\right)$ is a nonlinear function. Supposing $g\left(0\right)=0$ and for any $x_1\left(t\right)$, $x_2\left(t\right)$, the following inequality holds [8]:

$$∥g\left(x_1\left(t\right)\right)-g\left(x_2\left(t\right)\right)∥\le ∥U_2(x_1\left(t\right)-x_2\left(t\right))∥$$

(14)

where $U_2$ is a known matrix. According to (13), (10) can be rewritten as:

$$\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}=B\left(z\right)Dx\left(t\right)+h\left(Dx\left(t\right)\right)b_n\left(z\right)+\omega (z,u\left(t\right),F_1,F_2)$$

(15)

Because there exist nonlinear terms in the weighted dynamic system, the design of the nonlinear filter is the key in the process of fault detection and isolation. In this paper, the fault detection and fault diagnosis filters are designed, respectively, according to $u\left(t\right)$, as well as the output PDFs $\gamma (z,u\left(t\right),F_1,F_2)$, so as to achieve the task of fault isolation.

3. Fault Detection Filter Design

In order to detect the faults based on the changes of the output PDFs, we construct the following nonlinear filter:

$$\left\{\begin{array}{c}\dot{\widehat{x}}\left(t\right)=A\widehat{x}\left(t\right)+Gg\left(\widehat{x}\left(t\right)\right)+Hu\left(t\right)+L\epsilon \left(t\right)\hfill \\ \epsilon \left(t\right)={\int }_a^b\sigma \left(z\right)[\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}-\sqrt{\widehat{\gamma }(z,u\left(t\right))}]dz\hfill \\ \sqrt{\widehat{\gamma }(z,u\left(t\right))}=B\left(z\right)D\widehat{x}\left(t\right)+h\left(D\widehat{x}\left(t\right)\right)b_n\left(z\right)\hfill \end{array}\right.$$

(16)

where $\widehat{x}\left(t\right)$ is the estimated value of the state vector $x\left(t\right)$, $L\in R^{m\times p}$ is the gain of the detection observer to be determined, $\sigma \left(z\right)\in R^{p\times 1}$ is the pre-specified weighting vector defined on $[a,b]$, and the residual $\epsilon \left(t\right)$ represents the integral of the difference between the measured PDFs $\gamma (z,u\left(t\right),F_1,F_2)$ and the estimated PDFs $\widehat{\gamma }(z,u\left(t\right))$ [26]. Denote $e\left(t\right)=x\left(t\right)-\widehat{x}\left(t\right)$, then the first derivative of $e\left(t\right)$ with respect to time t is:

$$\begin{array}{c}\dot{e}\left(t\right)=\dot{x}\left(t\right)-\dot{\widehat{x}}\left(t\right)=(A-L{\Gamma }_1)e\left(t\right)+[Gg\left(x\left(t\right)\right)-Gg\left(\widehat{x}\left(t\right)\right)]\\ -L{\Gamma }_2[h\left(Dx\left(t\right)\right)-h\left(D\widehat{x}\left(t\right)\right)]+F_1+F_2-L\Delta \left(t\right)\end{array}$$

(17)

where

$$\left\{\begin{array}{c}{\Gamma }_1={\int }_a^b\sigma \left(z\right)B\left(z\right)dz\hfill \\ {\Gamma }_2={\int }_a^b\sigma \left(z\right)b_n\left(z\right)dz\hfill \\ \Delta \left(t\right)={\int }_a^b\sigma \left(z\right)\omega (z,u\left(t\right),F_1,F_2)dz\hfill \end{array}\right.$$

(18)

As shown in (10) and (16), it is clear that the residual $\epsilon \left(t\right)$ is a nonlinear function of $e\left(t\right)$, $x\left(t\right)$, and $\widehat{x}\left(t\right)$. According to (18), the residual $\epsilon \left(t\right)$ can be further expressed as

$$\begin{array}{cc}\hfill \epsilon \left(t\right)& ={\int }_a^b\sigma \left(z\right)B\left(z\right)De\left(t\right)dz+{\int }_a^b\sigma \left(z\right)[h\left(Dx\left(t\right)\right)-h\left(D\widehat{x}\left(t\right)\right)]b_n\left(z\right)dz\hfill \\ \hfill & +{\int }_a^b\sigma \left(z\right)\omega (z,u\left(t\right),F_1,F_2)dz\hfill \\ \hfill & ={\Gamma }_{1}e\left(t\right)+{\Gamma }_2[h\left(Dx\left(t\right)\right)-h\left(D\widehat{x}\left(t\right)\right)]+\Delta \left(t\right)\hfill \end{array}$$

(19)

Recall $\left|\omega (z,u\left(t\right),F_1,F_2)\right|<\delta$, and combine it with (18); we can obtain that

$$\left|\Delta \left(t\right)\right|<\tilde{\delta }$$

(20)

where $\tilde{\delta }=\delta ∥{\int }_a^b\sigma \left(z\right)dz∥$. In (16), $\sigma \left(z\right)$ can be any constant vector, and it is required that $(A,{\Gamma }_1)$ is observable. In the fault detection stage, according to Theorem 1 in [8], if there exist parameters ${\lambda }_i(i=1,2)$, matrices $P>0$, R, and constant $\eta >0$, then if $F_1=F_2=0$, the system (17) with gain $L=P^{-1}R$ is stable, and the error satisfies:

$$∥e\left(t\right)∥\le {\alpha }_0=max\left\{∥e\left(0\right)∥,2{\eta }^{-1}\tilde{\delta }∥R∥\right\}$$

(21)

In order to detect the faults, we select $\epsilon \left(t\right)$ as the residual signal. According to [8], the faults can be detected as follows:

$$\epsilon \left(t\right)>\alpha ={\alpha }_0(∥{\Gamma }_1∥+∥{\Gamma }_2∥∥U_1∥)+\tilde{\delta }$$

(22)

4. Fault Isolation Filter Design

Once the faults are detected according to the method in Section 2, it is necessary to carry out fault diagnosis and estimate the size of different faults, respectively, so as to achieve fault isolation. For this purpose, we construct the following adaptive filter as shown in (23):

$$\left\{\begin{array}{c}\dot{\widehat{x}}\left(t\right)=A\widehat{x}\left(t\right)+Gg\left(\widehat{x}\left(t\right)\right)+Hu\left(t\right)+\left[\begin{array}{cc}1& 1\end{array}\right]\left[\begin{array}{c}{\widehat{F}}_1\\ {\widehat{F}}_2\end{array}\right]+L\epsilon \left(t\right)\hfill \\ \left[\begin{array}{c}{\dot{\widehat{F}}}_1\\ {\dot{\widehat{F}}}_2\end{array}\right]=\left[\begin{array}{cc}-{\Lambda }_5& 0\\ 0& -{\Lambda }_7\end{array}\right]\left[\begin{array}{c}{\widehat{F}}_1\\ {\widehat{F}}_2\end{array}\right]+\left[\begin{array}{c}{\Lambda }_6\\ {\Lambda }_8\end{array}\right]\epsilon \left(t\right)\hfill \\ \epsilon \left(t\right)={\int }_a^b\sigma \left(z\right)[\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}-\sqrt{\widehat{\gamma }(z,u\left(t\right))}]dz\hfill \\ \sqrt{\widehat{\gamma }(z,u\left(t\right))}=B\left(z\right)D\widehat{x}\left(t\right)+h\left(D\widehat{x}\left(t\right)\right)b_n\left(z\right)\hfill \end{array}\right.$$

(23)

where ${\widehat{F}}_1$ and ${\widehat{F}}_2$ are the estimates of faults $F_1$ and $F_2$. In (23), ${\Lambda }_i(i=5,6,7,8)$ is the learning operator with respect to ${\widehat{F}}_1$, ${\widehat{F}}_2$, and the fault estimation errors.

Let $e\left(t\right)=x\left(t\right)-\widehat{x}\left(t\right)$, and define the fault estimation errors as ${\tilde{F}}_1=F_1-{\widehat{F}}_1$, ${\tilde{F}}_2=F_2-{\widehat{F}}_2$, then the estimation errors can be shown as in (24):

$$\left\{\begin{array}{c}\dot{e}\left(t\right)=(A-L{\Gamma }_1)e\left(t\right)+[Gg\left(x\left(t\right)\right)-Gg\left(\widehat{x}\left(t\right)\right)]\hfill \\ -L{\Gamma }_2[h\left(Dx\left(t\right)\right)-h\left(D\widehat{x}\left(t\right)\right)]+{\tilde{F}}_1+{\tilde{F}}_2-L\Delta \left(t\right)\hfill \\ {\dot{\tilde{F}}}_1=-{\Lambda }_5{\tilde{F}}_1+{\Lambda }_5{F}_1-{\Lambda }_6\epsilon \left(t\right)\hfill \\ {\dot{\tilde{F}}}_2=-{\Lambda }_7{\tilde{F}}_2+{\Lambda }_7{F}_2-{\Lambda }_8\epsilon \left(t\right)\hfill \\ \epsilon \left(t\right)={\Gamma }_{1}e\left(t\right)+{\Gamma }_2[h\left(Dx\left(t\right)\right)-h\left(D\widehat{x}\left(t\right)\right)]+\Delta \left(t\right)\hfill \end{array}\right.$$

(24)

Supposing $∥F_1∥\le M_1/2$, $∥F_2∥\le M_2/2$, the following theorem will show that, by selecting the appropriate filter gains ${\Lambda }_5$, ${\Lambda }_6$, ${\Lambda }_7$, and ${\Lambda }_8$, the fault estimation errors can be controlled in a small range.

Theorem 1. 

If there exist ${\lambda }_i>0(i=1,2)$, matrices $P>0$, R, and ${\Lambda }_i(i=5,6,7,8)$, and constants $\kappa >0$, ${\theta }_1>0$, ${\theta }_2>0$, and ${\theta }_3$ satisfying:

$$\left[\begin{array}{ccccccc}{\Pi }_0+\kappa I& P-{\Gamma }_1^T{\Lambda }_6^T& P-{\Gamma }_1^T{\Lambda }_8^T& {\Pi }_2& 0& 0& D^T{U}_1^T\\ P-{\Lambda }_6{\Gamma }_1& -2{\Lambda }_5^T& 0& 0& {\Pi }_3& 0& 0\\ P-{\Lambda }_8{\Gamma }_1& 0& -2{\Lambda }_7^T& 0& 0& {\Pi }_4& 0\\ {\Pi }_2^T& 0& 0& -I& 0& 0& 0\\ 0& {\Pi }_3^T& 0& 0& -I& 0& 0\\ 0& 0& {\Pi }_4^T& 0& 0& -I& 0\\ U_{1}D& 0& 0& 0& 0& 0& -\frac{1}{2}{\theta }_3^{2}I\end{array}\right]<0$$

(25)

where

$$\begin{array}{c}{\Pi }_2=\left[\begin{array}{ccc}{\lambda }_{1}R{\Gamma }_2& {\lambda }_{2}PG& {\theta }_1\end{array}R\right]\\ {\Pi }_3=\left[\begin{array}{cc}{\theta }_2{\Lambda }_6& {\theta }_3{\Lambda }_6{\Gamma }_2\end{array}\right]\\ {\Pi }_4=\left[\begin{array}{cc}{\theta }_2{\Lambda }_8& {\theta }_3{\Lambda }_8{\Gamma }_2\end{array}\right]\\ {\Pi }_0=(PA-R{\Gamma }_1)+{(PA-R{\Gamma }_1)}^T+\frac{1}{{\lambda }_1^2}{D}^T{U}_1^T{U}_{1}D+\frac{1}{{\lambda }_2^2}{U}_2^T{U}_2\end{array}$$

then, with gain $L=P^{-1}R$, the error system (24) is stable in the presence of $F_1$, $F_2$, and the estimation error satisfies

$${∥e\left(t\right)∥}^2\le min\{{∥e\left(0\right)∥}^2,{\kappa }^{-1}(({\theta }_1^{-2}+2{\theta }_2^{-2}){\tilde{\delta }}^2+∥{\Lambda }_5∥{M_1}^2+∥{\Lambda }_7∥{M_2}^2)\}$$

(26)

for all $t\in [0,+\infty )$.

Proof. 

Consider the following Lyapunov function:

$$\Pi \left(e\left(t\right),x\left(t\right),\widehat{x}\left(t\right),{\tilde{F}}_1,{\tilde{F}}_2,t\right)=\Phi \left(e\left(t\right),x\left(t\right),\widehat{x}\left(t\right),t\right)+{\tilde{F}}_1^T{\tilde{F}}_1+{\tilde{F}}_2^T{\tilde{F}}_2$$

(27)

where

$$\begin{array}{cc}\hfill \Phi (e\left(t\right),x\left(t\right),\widehat{x}\left(t\right),t)& =e^T\left(t\right)Pe\left(t\right)+\frac{1}{{\lambda }_2^2}{\int }_0^t[{∥U_{2}e(\tau ∥}^2-∥gx\left(\tau \right)-g(\widehat{x}\left(\tau \right)∥]d\tau \hfill \\ \hfill & +\frac{1}{{\lambda }^2}{\int }_0^t[{∥U_{1}De\left(\tau \right)∥}^2-{∥h\left(Dx\left(\tau \right)\right)-h(D\widehat{x}\left(\tau \right)∥}^2]d\tau \hfill \end{array}$$

(28)

To simplify the proof, we abbreviate $\Phi \left(e\left(t\right),x\left(t\right),\widehat{x}\left(t\right),t\right)$ to $\Phi$. Let $R=PL$; we obtain:

$$\begin{array}{c}\dot{\Phi }\le e^T\left(t\right){\Psi }_{0}e\left(t\right)-2e^T\left(t\right)PL\Delta \left(t\right)+2e^T\left(t\right)P{\tilde{F}}_1+2e^T\left(t\right)P{\tilde{F}}_2\\ \le e^T\left(t\right){\Psi }_{1}e\left(t\right)+{\theta }_1^{-2}{\Delta }^T\left(t\right)\Delta \left(t\right)+2e^T\left(t\right)P{\tilde{F}}_1+2e^T\left(t\right)P{\tilde{F}}_2\end{array}$$

(29)

where

$$\begin{array}{c}{\Psi }_1={\Psi }_0+{\theta }_1^{2}R{R}^T\hfill \\ {\Psi }_0=P\left(A-L{\Gamma }_1\right)+{\left(A-L{\Gamma }_1\right)}^{T}P+{\lambda }_2^{2}PG{G}^{T}P\hfill \\ +{\lambda }_1^{2}PL{\Gamma }_2{\Gamma }_{2}LP+\frac{1}{{\lambda }_1^2}{D}^T{U}_1^T{U}_{1}D+\frac{1}{{\lambda }_2^2}{U}_2^T{U}_2\hfill \end{array}$$

From (24) and (29), the first derivative of $\Pi$ is:

$$\begin{array}{cc}\hfill \dot{\Pi }& =e^T\left(t\right){\Psi }_{1}e\left(t\right)+{\theta }_1^{-2}{\Delta }^T\left(t\right)\Delta \left(t\right)+2e^T\left(t\right)P{\tilde{F}}_1+2e^T\left(t\right)P{\tilde{F}}_2+2{\tilde{F}}_1^T{\dot{\tilde{F}}}_1+2{\tilde{F}}_2^T{\dot{\tilde{F}}}_2\hfill \\ \hfill & \le e^T\left(t\right)[{\Psi }_1+2{\theta }_3^{-2}{D}^T{U}_1^T{U}_{1}D]e\left(t\right)+2e^T\left(t\right)P{\tilde{F}}_1+2e^T\left(t\right)P{\tilde{F}}_2-2{\tilde{F}}_1^T{\Lambda }_5^T{\tilde{F}}_1\hfill \\ \hfill & +{\theta }_2^2{\tilde{F}}_1^T{\Lambda }_6{\Lambda }_6^T{\tilde{F}}_1+{\theta }_2^2{\tilde{F}}_2^T{\Lambda }_8{\Lambda }_8^T{\tilde{F}}_2-2{\tilde{F}}_2{\Lambda }_7^T{\tilde{F}}_2+({\theta }_1^{-2}+2{\theta }_2^{-2}){\Delta }^T\left(t\right)\Delta \left(t\right)\hfill \\ \hfill & -2{\tilde{F}}_1^T{\Lambda }_6{\Gamma }_{1}e\left(t\right)-2{\tilde{F}}_2^T{\Lambda }_8{\Gamma }_{1}e\left(t\right)+{\theta }_3^2{\tilde{F}}_1^T{\Lambda }_6{\Gamma }_2{\Gamma }_2^T{\Lambda }_6^T{\tilde{F}}_1+{\theta }_3^2{\tilde{F}}_2^T{\Lambda }_8{\Gamma }_2{\Gamma }_2^T{\Lambda }_8^T{\tilde{F}}_2\hfill \\ \hfill & +2F_1^T{\Lambda }_5^T{\tilde{F}}_1+2F_2^T{\Lambda }_7^T{\tilde{F}}_2\hfill \\ \hfill & =\left[\begin{array}{ccc}{e}^T\left(t\right)& {\tilde{F}}_1^T& {\tilde{F}}_2^T\end{array}\right]\overline{\Psi }\left[\begin{array}{c}e\left(t\right)\\ {\tilde{F}}_1\\ {\tilde{F}}_2\end{array}\right]+({\theta }_1^{-2}+2{\theta }_2^{-2}){\Delta }^T\left(t\right)\Delta \left(t\right)+2F_1^T{\Lambda }_5^T{\tilde{F}}_1+2F_2^T{\Lambda }_7^T{\tilde{F}}_2\hfill \end{array}$$

where

$$\begin{array}{c}{\overline{\Psi }}_{22}=-2{\Lambda }_5^T+{\theta }_2^2{\Lambda }_6{\Lambda }_6^T+{\theta }_3^2{\Lambda }_6{\Gamma }_2{\Gamma }_2^T{\Lambda }_6^T\hfill \\ {\overline{\Psi }}_{33}=-2{\Lambda }_7^T+{\theta }_2^2{\Lambda }_8{\Lambda }_8^T+{\theta }_3^2{\Lambda }_8{\Gamma }_2{\Gamma }_2^T{\Lambda }_8^T\hfill \\ \overline{\Psi }=\left[\begin{array}{c}\begin{array}{ccc}{\Psi }_1+2{\theta }_3^{-2}{D}^T{U}_1^T{U}_{1}D& P-{\Gamma }_1^T{\Lambda }_6^T& P-{\Gamma }_1^T{\Lambda }_8^T\\ P-{\Lambda }_6{\Gamma }_1& {\overline{\Psi }}_{22}& 0\\ P-{\Lambda }_8{\Gamma }_1& 0& {\overline{\Psi }}_{33}\end{array}\end{array}\right]\hfill \end{array}$$

Using the Schur complement lemma, we can obtain that (25) is equivalent to $\overline{\Psi }\le diag\{-\kappa I,0\}$, from which we can obtain:

$$\begin{array}{c}\dot{\Pi }<-\kappa {∥e\left(t\right)∥}^2+\left({\theta }_1^{-2}+2{\theta }_2^{-2}\right){\tilde{\delta }}^2+2F_1^T{\Lambda }_5^T{\tilde{F}}_1+2F_2^T{\Lambda }_7^T\tilde{F}\\ \le -\kappa {∥e\left(t\right)∥}^2+\left({\theta }_1^{-2}+2{\theta }_2^{-2}\right){\tilde{\delta }}^2+∥{\Lambda }_5∥M_1^2+∥{\Lambda }_7∥M_2^2\end{array}$$

When $\kappa {∥e\left(t\right)∥}^2>\left({\theta }_1^{-2}+2{\theta }_2^{-2}\right){\tilde{\delta }}^2+∥{\Lambda }_5∥M_1^2+∥{\Lambda }_7∥M_2^2$, $\Pi >0$, $\dot{\Pi }<0$, so (22) holds. □

5. Simulation

In the optical fiber vibration source detection system [22] shown in Figure 1, the optical fiber is buried underground in the same ditch with oil or gas pipelines. The PDFs of the output light intensity are affected by the input light intensity $u\left(t\right)$, false alarms (environmental interferences, such as a vehicle passing), as well as various destructive vibration sources (theft, geological disasters, etc.). Our goal is to detect whether these vibration sources are false alarms or real alarms. It is supposed that the PDFs of the output light intensity can be approximately expressed by the square root B-spline basis function described by $\sqrt{\gamma (z,u\left(t\right),F_1,F_2)}={\sum }_1^3{v}_i(z,u\left(t\right),F_1,F_2)b_i\left(z\right)$, where $F_1$ represents pickaxe digging and $F_2$ represents machine excavation. Suppose

$$b_i\left(z\right)=\left\{\begin{array}{c}\left|sin2\pi z\right|,z\in [0.5\left(i-1\right),0.5i]\hfill \\ 0,z\in [0.5\left(j-1\right),0.5j]\hfill \end{array}\right.\hspace{1em}\hspace{1em}i\ne j$$

(30)

for $i=1,2,3$. Recalling (2)–(4), it can be calculated that

$${\Lambda }_1=\left[\begin{array}{cc}0.25& 0\\ 0& 0.25\end{array}\right],{\Lambda }_2=\left[\begin{array}{cc}0& 0\end{array}\right],{\Lambda }_3=0.25.$$

In this simulation, the parameters were selected as follows:

$$\begin{array}{c}A=\left[\begin{array}{cc}-0.5& 0\\ 0& -1.3\end{array}\right],G=\left[\begin{array}{cc}0& 0\\ 0& 0.1\end{array}\right],H=\left[\begin{array}{cc}0.2& 0\\ 0& -0.3\end{array}\right]\hfill \\ G\left(V\left(t\right)\right)=\left[\begin{array}{c}0\\ 0.25\sqrt{v_1^2\left(t\right)+v_2^2\left(t\right)}\end{array}\right],D=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\hfill \end{array}$$

Consequently, it can be calculated that $U_1=\left[\begin{array}{cc}1& 1\end{array}\right]$ and $U_2=\left[\begin{array}{cc}0& 0\\ 0& 0.5\end{array}\right]$. According to (16), a nonlinear detection filter can be constructed. Select $\sigma =1$; it can be verified that $(A,{\Gamma }_1)$ is observable and

$${\Gamma }_1=\left[\begin{array}{cc}\frac{1}{\pi }& \frac{1}{\pi }\end{array}\right],{\Gamma }_1=\frac{1}{\pi },\Delta \left(t\right)\le 0.15.$$

The initial values are supposed to be $x\left(0\right)=\left[\begin{array}{cc}0.5& 0.25\end{array}\right]$, $\widehat{x}\left(0\right)=\left[\begin{array}{cc}0& 0\end{array}\right]$. For ${\lambda }_1={\lambda }_2=1$, ${\theta }_1={\theta }_2=2$, ${\theta }_3=-2$, $\eta =2$, and $\kappa =0.1$. From (22), the thresholds are ${\alpha }_0=0.4719$ and $\alpha =0.5751$. By solving (25) with the LMI toolbox, we obtain:

$$\begin{array}{c}P=\left[\begin{array}{cc}3.9269& -4.5402\\ -4.5402& 6.4692\end{array}\right]R=\left[\begin{array}{c}0.1106\\ 0.1104\end{array}\right]L=\left[\begin{array}{c}0.0527\\ 0.0212\end{array}\right]\hfill \\ {\Lambda }_5=\left[\begin{array}{cc}7.1833& -8.5760\\ -4.3862& 9.5318\end{array}\right]{\Lambda }_6=\left[\begin{array}{c}0.4681\\ -0.2479\end{array}\right]\hfill \\ {\Lambda }_7={10}^6\times \left[\begin{array}{cc}30.330& -0.1912\\ -0.1912& 3.7291\end{array}\right]{\Lambda }_8=\left[\begin{array}{c}6.6378\\ -2.0262\end{array}\right]\hfill \end{array}$$

Figure 2 shows the PDFs of the output light intensity when there is no destructive vibration source (i.e., $F_1=F_2=0$), and Figure 3 shows the PDFs of the output light intensity when some destructive vibration source occurs (i.e., $F_1\ne 0$ and $F_2\ne 0$). Comparing Figure 2 and Figure 3, it can be seen that the PDFs of the output light intensity have changed significantly. In the simulation, it was assumed that the first destructive vibration source (pickaxe digging) starts from the 10th second and the second destructive vibration source (machine excavation) starts from the 30th second. Based on the fault detection filter (16) in this paper, the response of the residual $\epsilon \left(t\right)$ is shown in Figure 4. It can be seen from Figure 4 that, when these faults occur, the value of the residual $\epsilon \left(t\right)$ will change, but it always satisfies $\epsilon \left(t\right)>0.5751$. When $\epsilon \left(t\right)>0.5751$, an alarm should be sounded in the optical fiber vibration source detection system.

After the fault diagnosis filter is designed based on Theorem 1, we can obtain Figure 5, Figure 6 and Figure 7. Figure 5 is the estimation of $F_1$ when pickaxe digging occurs. It can be seen from Figure 5 that the fault diagnosis filter can quickly track the change of the fault 10 s after the first destructive vibration source $F_1$ appears. Figure 6 is the response of the fault diagnosis observer when machine excavation $F_2$ exists. It can be seen from Figure 6 that, after the second vibration source appears, the fault diagnosis filter can quickly track fault $F_2$. Figure 7 is the response of the fault diagnosis observer with multiple faults $F_1$ (pickaxe digging) and $F_2$ (machine excavation), which can be detected, respectively, based on the methods in this paper. It can be seen from Figure 7 that, when faults $F_1$ and $F_2$ exist at the same time, the diagnostic filter can clearly diagnose $F_1$ and $F_2$. The corresponding vibration source strength can thus be estimated and helps to identify which type of damage.

In the man–machine interface designed for the optical fiber vibration source detection and early warning system, we carried out a comparative experiment. We used the same set of data to compare the PDF thresholds. The experiment showed that, before adjusting the threshold value, the detection effect of the optical fiber signal under the influence of vibration is as shown in Figure 8. The platform detects that three columns have vibration alarm signals, and this section flashes to remind the user. The effect after using the method in this paper is shown in Figure 9, which proves that the method of fault isolation can reduce the false alarm rate and make the detection effect of the optical fiber vibration source more stable and accurate.

6. Conclusions

In this paper, in order to reduce the false alarm rate in oil or gas pipeline monitoring, we regarded the different vibrations in the system as faults and considered faults’ isolation to identify different types of vibrations for a non-Gaussian stochastic distribution control model because of inevitable non-stationary interferences from the environment. The fault isolation started by using square root B-spline expansion and a nonlinear weighted dynamic model. The faults in the system were estimated separately to achieve the purpose of fault isolation. The output value of the system in this paper is not measurable, but the output PDFs can be measured by optical instruments. Firstly, the PDFs of the output signal were approximated by square root B-spline expansion. Secondly, the nonlinear dynamic model between the control input and the weights of PDFs was established by a dynamic neural network. Thirdly, based on the measured output PDFs and the input of the system, a filter-based residual generator was constructed to detect and diagnose the faults. Through LMIs, the feasibility criterion for the detection and isolation system faults was given, as well as the steady error ranges. Finally, the effectiveness of this method was verified for the optical fiber vibration source detection system.

7. Future Work

In the optical fiber vibration source detection system, there will be more types of faults, such as time-varying faults, multiplicative faults, etc. This paper only focused on the detection and isolation of constant faults in optical systems. In the future work, we will continue to study the detection and isolation of multiplicative faults and time-varying faults in stochastic distributed systems.