# Subspace Method Aided Data-Driven Fault Detection Based on Principal Component Analysis

Authors: Lingling Ma; Xiangshun Li
Journal: Journal of Control Science and Engineering
 (2017)

---

## Abstract

The model-based fault detection technique, which needs to identify the system models, has been well established. The objective of this paper is to develop an alternative procedure instead of identifying the system models. In this paper, subspace method aided data-driven fault detection based on principal component analysis (PCA) is proposed. The basic idea is to use PCA to identify the system observability matrices from input and output data and construct residual generators. The advantage of the proposed method is that we just need to identify the parameterized matrices related to residuals rather than the system models, which reduces the computational steps of the system. The proposed approach is illustrated by a simulation study on the Tennessee Eastman process.


---

## Body

## 1. Introduction

The safety, stability, and efficiency of dynamic systems have always been a matter of great concern in the field of complex industrial processes. Fault detection plays an extremely important role in improving the safety of processes and attracts more and more attention in the field of process monitoring systems.During the past three decades, a model-based fault detection technique for linear time invariant (LTI) systems has been well established [1–4]. But in some complex systems and large-scale industries, it is difficult to obtain accurate mathematical models. With the development of the computer and information industry, Big Data has attracted more and more attention. Actually, industrial processes have large amounts of operating data which contain abundant information about the system. Subspace model identification which uses process data to identify system models comes into being under this kind of background. There are several methods of subspace algorithms, such as MOESP [5], N4SID [6], and CVA [7]. In a typical SMI, the identification procedure comprises two steps. Firstly, the extended observability matrix Γf and a block triangular Toeplitz matrix Hf are identified from input and output data. Then, A, B, C, and D are calculated from the identified observability matrix and the Toeplitz matrix. Fault detection using this method usually comprises two steps: (1) system identification and (2) model-based fault detection.Residual generation and residual evaluation are essential for model-based fault detection. Fault detection can be achieved if residuals are obtained. Parity space approach (PSA) [8] is the simplest and most widely used among model-based fault detection methods. We noticed that there is a close bond between the two methods by comparing PSA and SMI. The residual could be calculated by some parameterized matrices related to Γf and Hf. So, it is feasible that we only need to use the first step of SMI to identify parameterized matrices instead of the system models, which reduces the calculation steps. To the best of our knowledge, Ding and his coworkers subtly associated SMI with fault detection and proposed a subspace aided approach [9–11]. But most of these SMI methods are biased under the errors-invariables (EIV) situation. Wang and Qin [12] proposed a subspace identification approach based on PCA that gives consistent model estimates under the EIV situation. Inspired by their work, we propose subspace method aided data-driven fault detection based on PCA in this paper. Figure 1 shows the difference between the proposed approach and the classic model-based approach.Figure 1
Model-based approach and the proposed approach.In this paper, the main contribution of this proposed approach lies in the direct identification of the parameterized matrices by PCA. As long as these matrices related to residuals are identified, we can construct residual generators. The proposed approach is illustrated by a simulation study on the Tennessee Eastman process.The rest of the paper is organized as follows. Section2 gives the preliminaries and problem formulation. The identification of parameterized matrices is presented in Section 3. Section 4 illustrates a simulation study on the Tennessee Eastman process. Finally, the conclusions are presented in Section 5.


## 2. Preliminaries and Problem Formulation

### 2.1. Process Descriptions

The form of the state space representation of a discrete time LTI system is given by(1)xk+1=Axk+Bu∗k+pky∗k=Cxk+Du∗k,where x∈Rn, u∗(k)∈Rl, y∗(k)∈Rm, and p(k)∈Rn denote the vectors of the state variables, noise-free inputs, noise-free outputs, and process noise, respectively. The available input measurements u(k) and output measurements y(k) are described by [13](2)uk=u∗k+okyk=y∗k+vk,where o(k)∈Rl and v(k)∈Rm represent the input and output noise. The following assumptions are introduced: (i)
The system is observable and controllable.(ii)
o(k), v(k), and p(k) are assumed to be zero-mean, normally distributed white noise; that is, (3)EpkokvkpjTojTvjT=RppRpoRpvRpoTRooRovRpvTRovTRvvδkjδkj=1,k=j0,k≠j.(iii)
All the noise is assumed to be independent of the past noise-free inputu∗(k) and initial state x(1). They satisfy (4)Eu∗kpjojvjT=0,forj⩾k.(iv)
System matricesA, B, C, and D and system order n are unknown.We define the following vectors and Hankel matrices. The subscriptsp and f represent the past and the future, where p⩾f⩾n. (5)upk=uk-puk-p+1⋮uk-1∈Rlp,ufk=ukuk+1⋮uk+f-1∈Rlf,ypk=yk-pyk-p+1⋮yk-1∈Rmp,yfk=ykyk+1⋮yk+f-1∈Rmf,Up=upkupk+1⋯upk+N-1∈Rlp×N,Uf=ufkufk+1⋯ufk+N-1∈Rlf×N,Yp=ypkypk+1⋯ypk+N-1∈Rmp×N,Yf=yfkyfk+1⋯yfk+N-1∈Rmf×N,Xk=xkxk+1⋯xk+N-1∈Rn×N,Zp=YpUp∈Rlp+mp×N,Zf=YfUf∈Rlf+mf×N.By iterating (1) and (2), we can get(6)yfk=Γfxk+Hfufk-Hfvfk+Gfpfk+ofk.The vectorsvf(k)∈Rlf, of(k)∈Rmf, pf(k)∈Rnf, and uf(k)∈Rlf have similar structures to yf(k). If we use Hankel matrix instead of the vectors to describe the system, (6) can be written as(7)Yf=ΓfXk+HfUf+EfEf=-HfVf+GfPf+Of.The matricesVf and Of are defined similarly to Yf. (8)Γf=CCA⋮CAf-1∈Rmf×n,Hf=D0⋯0CBD⋯0⋮⋯⋱⋮CAf-2BCAf-3B⋯D∈Rmf×lf,Gf=00⋯0C0⋯0⋮⋯⋱⋮CAf-2CAf-3⋯0∈Rmf×nf.


### 2.2. Residual Generation and Evaluation

Residual generation and residual evaluation are essential to design a fault detection system. The main idea of this paper is to generate residual from input and output data. How to generate and evaluate residual is introduced in the following section.MovingHfUf to the left of the equation, (7) can be rewritten as (9)Yf-HfUf=ΓfXk+Ef.TakingΓf⊥∈Rmf×(mf-n) as the orthogonal complement of Γf, we can obtain (10)Γf⊥TΓf=0.SinceX(k) is unknown, we multiply (Γf⊥)T on both sides of the equation to eliminate the effect of the unknown state [14]. Then, we can get(11)Γf⊥TYf-HfUf=Γf⊥TEf.So, the residual signal(12)R=Γf⊥TYf-HfUf=Γf⊥TYf-Γf⊥THfUf.We noticed that the residual is closely related to the extended observabilityΓf and a block triangular Toeplitz matrix Hf. Therefore, if we want to obtain the residual, the key step is to identify (Γf⊥)T and -(Γf⊥)THf. In Section 3, how to identify these parameterized matrices will be introduced in detail.Once these parameterized matrices are identified, we can get residual signal. Then, the fault detection can be completed by residual evaluation. Letr(k) be the residual vector at the instant k, defining the testing statistic [15](13)J=rkTΣres-1rkΣres=RRTN-1,whereΣres is the covariance matrix of the residual matrix in the fault-free case. Note that when there is no fault, J~χ2(h), where χ2(h) is the chi-squared distribution with h degrees of freedom. We set Jth=χα2(h) as threshold and define the detection logic as (14)J=<Jth,no fault>Jth,a fault is detected.Besides the chi-squared distribution, kernel density estimation is also utilized to determine the threshold [16]. It is assumed that the measurements follow a Gaussian distribution. The threshold Jth can be calculated by P(J<Jth(α))=α, where α is confidence level. (15)Px<a=∫-∞apxdx=αpx=1Nh∑k=1NKx-xkh,where xk is the kth sample of x, N is the number of samples, h is a smoothing parameter, and K(·) is a kernel function.Given a proper confidence level, we can get the thresholdJth. In this paper, we set α=0.95. After the statistics and threshold are deduced, the fault detection can be realized according to the above detection logic.


## 2.1. Process Descriptions

The form of the state space representation of a discrete time LTI system is given by(1)xk+1=Axk+Bu∗k+pky∗k=Cxk+Du∗k,where x∈Rn, u∗(k)∈Rl, y∗(k)∈Rm, and p(k)∈Rn denote the vectors of the state variables, noise-free inputs, noise-free outputs, and process noise, respectively. The available input measurements u(k) and output measurements y(k) are described by [13](2)uk=u∗k+okyk=y∗k+vk,where o(k)∈Rl and v(k)∈Rm represent the input and output noise. The following assumptions are introduced: (i)
The system is observable and controllable.(ii)
o(k), v(k), and p(k) are assumed to be zero-mean, normally distributed white noise; that is, (3)EpkokvkpjTojTvjT=RppRpoRpvRpoTRooRovRpvTRovTRvvδkjδkj=1,k=j0,k≠j.(iii)
All the noise is assumed to be independent of the past noise-free inputu∗(k) and initial state x(1). They satisfy (4)Eu∗kpjojvjT=0,forj⩾k.(iv)
System matricesA, B, C, and D and system order n are unknown.We define the following vectors and Hankel matrices. The subscriptsp and f represent the past and the future, where p⩾f⩾n. (5)upk=uk-puk-p+1⋮uk-1∈Rlp,ufk=ukuk+1⋮uk+f-1∈Rlf,ypk=yk-pyk-p+1⋮yk-1∈Rmp,yfk=ykyk+1⋮yk+f-1∈Rmf,Up=upkupk+1⋯upk+N-1∈Rlp×N,Uf=ufkufk+1⋯ufk+N-1∈Rlf×N,Yp=ypkypk+1⋯ypk+N-1∈Rmp×N,Yf=yfkyfk+1⋯yfk+N-1∈Rmf×N,Xk=xkxk+1⋯xk+N-1∈Rn×N,Zp=YpUp∈Rlp+mp×N,Zf=YfUf∈Rlf+mf×N.By iterating (1) and (2), we can get(6)yfk=Γfxk+Hfufk-Hfvfk+Gfpfk+ofk.The vectorsvf(k)∈Rlf, of(k)∈Rmf, pf(k)∈Rnf, and uf(k)∈Rlf have similar structures to yf(k). If we use Hankel matrix instead of the vectors to describe the system, (6) can be written as(7)Yf=ΓfXk+HfUf+EfEf=-HfVf+GfPf+Of.The matricesVf and Of are defined similarly to Yf. (8)Γf=CCA⋮CAf-1∈Rmf×n,Hf=D0⋯0CBD⋯0⋮⋯⋱⋮CAf-2BCAf-3B⋯D∈Rmf×lf,Gf=00⋯0C0⋯0⋮⋯⋱⋮CAf-2CAf-3⋯0∈Rmf×nf.


## 2.2. Residual Generation and Evaluation

Residual generation and residual evaluation are essential to design a fault detection system. The main idea of this paper is to generate residual from input and output data. How to generate and evaluate residual is introduced in the following section.MovingHfUf to the left of the equation, (7) can be rewritten as (9)Yf-HfUf=ΓfXk+Ef.TakingΓf⊥∈Rmf×(mf-n) as the orthogonal complement of Γf, we can obtain (10)Γf⊥TΓf=0.SinceX(k) is unknown, we multiply (Γf⊥)T on both sides of the equation to eliminate the effect of the unknown state [14]. Then, we can get(11)Γf⊥TYf-HfUf=Γf⊥TEf.So, the residual signal(12)R=Γf⊥TYf-HfUf=Γf⊥TYf-Γf⊥THfUf.We noticed that the residual is closely related to the extended observabilityΓf and a block triangular Toeplitz matrix Hf. Therefore, if we want to obtain the residual, the key step is to identify (Γf⊥)T and -(Γf⊥)THf. In Section 3, how to identify these parameterized matrices will be introduced in detail.Once these parameterized matrices are identified, we can get residual signal. Then, the fault detection can be completed by residual evaluation. Letr(k) be the residual vector at the instant k, defining the testing statistic [15](13)J=rkTΣres-1rkΣres=RRTN-1,whereΣres is the covariance matrix of the residual matrix in the fault-free case. Note that when there is no fault, J~χ2(h), where χ2(h) is the chi-squared distribution with h degrees of freedom. We set Jth=χα2(h) as threshold and define the detection logic as (14)J=<Jth,no fault>Jth,a fault is detected.Besides the chi-squared distribution, kernel density estimation is also utilized to determine the threshold [16]. It is assumed that the measurements follow a Gaussian distribution. The threshold Jth can be calculated by P(J<Jth(α))=α, where α is confidence level. (15)Px<a=∫-∞apxdx=αpx=1Nh∑k=1NKx-xkh,where xk is the kth sample of x, N is the number of samples, h is a smoothing parameter, and K(·) is a kernel function.Given a proper confidence level, we can get the thresholdJth. In this paper, we set α=0.95. After the statistics and threshold are deduced, the fault detection can be realized according to the above detection logic.


## 3. The Identification of Parameterized Matrices

Based on the knowledge that has been introduced in the previous sections, we know that the identification of matrices related to residual is the key step in this approach. In this section, we discuss the issue of how to identify these matrices. Wang and Qin proposed a subspace identification approach based on PCA. We will briefly introduce it. When the process is corrupted by process noise and measurement simultaneously, they use instrumental variablesZp to remove the noise, because the past data Zp is uncorrelated with future noise. According to these assumptions, we can get(16)limN→∞⁡1NEfZpT=limN→∞⁡1N-HfVfZpT+GfPfZpT+OfZpT=0.If we utilizeZf to represent residual instead of Yf and Uf, (11) can be expressed as (17)Γf⊥TI∣-HfZf=Γf⊥TEf.Multiplying1/NZpT on both sides of the equation, we can get (18)1NΓf⊥TI∣-HfZfZpT=1NΓf⊥TEfZpT.From (16), we know(19)1NΓf⊥TI∣-HfZfZpT=0.Equation (19) indicates that 1/NZfZpT has zero scores. From the rank condition (20)rank⁡limN→∞⁡1NEfZpT=lf+n,we know that when N→∞, 1/NZfZpT has mf-n zero scores. Performing PCA on the process data, (21)limN→∞⁡1NZfZpT=PTT+P~T~T.When N→∞,(22)Γf⊥-HfTΓf⊥=P~M,P~=P~yP~u,where M∈R(mf-n)(mf-n) is a nonsingular matrix, P~y is the first mf rows of P~, and P~u is the last lf rows of P~. The nonsingular matrix M is unit matrix in our paper. The matrices Γf⊥ and -HfTΓf⊥ can be obtained according to (22), and then the residual can be generated. The algorithm from data to realize fault detection can be summarized as the steps shown in Algorithm 1.Algorithm 1:SMI-PCA fault detection.
(
1
) Generate Zf and Zp with fault-free data
(
2
) Perform PCA on 1/NZfZpT
(
3
) Calculate Γf⊥ and -HfTΓf⊥ according to (22)
(
4
) Construct residual R according to (12)
(
5
) Calculate Σres=RRT/(N-1) in the fault-free case
(
6
) Generate Yp and Up with on-line operating data
(
7
) Construct residual R according to (12), then calculate
the testing statisticJ by (13)
(
8
) Determine the threshold Jth
(
9
) Evaluate the residual according to the detection logic


## 4. Simulation Study on the Tennessee Eastman Process

In this section, we apply the proposed approach to Tennessee Eastman (TE) process. Three indices are used—fault detection rate (FDR), false alarm rate (FAR) [17], and fault detection time (FDT)—to demonstrate the efficiency of the proposed approach. The fault detection time is the first time instance with the testing statistic J above the thresholds Jth. (23)FDR=No. of samplesJ>Jth∣f≠0total samplesf≠0×100FAR=No. of samplesJ>Jth∣f=0total samplesf=0×100.


### 4.1. TE Process

TE process model is a realistic simulation program of a chemical process that is widely accepted as a benchmark for control monitoring studies [18]. The TE model that we used was downloaded from http://depts.washington.edu/control/LARRY/TE. The flow diagram of the process is shown in Figure 2. The process contains five major units named reactor, condenser, compressor, separator, and stripper.Figure 2
The Tennessee Eastman process.There are 53 measurements, 12 of which are manipulated variables and 41 are process variables. The 12 manipulated variables are listed in Table1. In this paper, the other 9 manipulated variables except for XMV(5), XMV(9), and XMV(12) are taken as inputs. The dimension of the TE dataset is very large. To overcome the curse of dimensionality, the 41 process variables are divided into eight blocks [19]. We choose the input feed block as outputs, which are shown in Table 2. Table 3 shows the 20 process faults of TE process.Table 1
Process manipulated variables.
Variable name
Variable number
D feed flow
XMV
(
1
)
E feed flow
XMV
(
2
)
A feed flow
XMV
(
3
)
A and C feed flow
XMV
(
4
)
Compressor recycle valve
XMV
(
5
)
Purge valve
XMV
(
6
)
Separator pot liquid flow
XMV
(
7
)
Stripper liquid product flow
XMV
(
8
)
Stripper steam valve
XMV
(
9
)
Reactor cooling water flow
XMV
(
10
)
Condenser cooling water flow
XMV
(
11
)
Agitator speed
XMV
(
12
)Table 2
Input feed block process measurements.
Variable name
Variable number
A feed (stream1)
XMEAS
(
1
)
D feed (stream1)
XMEAS
(
2
)
E feed (stream1)
XMEAS
(
3
)
A and C feed
XMEAS
(
4
)Table 3
Process faults.
Variable number
Process variable
Type
IDV(
0)
Normal operation
Step
IDV(
1)
A
/
C feed ratio, B composition constant
Step
IDV(
2)
B composition, A/C ratio constant
Step
IDV(
3)
D feed temperature
Step
IDV(
4)
Reactor cooling water inlet temperature
Step
IDV(
5)
Condenser cooling water inlet temperature
Step
IDV(
6)
A feed loss
Step
IDV(
7)
C header pressure loss-reduced availability
Step
IDV(
8)
A
,
B
,
C feed composition
Random variation
IDV(
9)
D feed temperature
Random variation
IDV(
10)
C feed temperature
Random variation
IDV(
11)
Reactor cooling water inlet temperature
Random variation
IDV(
12)
Condenser cooling water inlet temperature
Random variation
IDV(
13)
Reaction kinetics
Random variation
IDV(
14)
Reactor cooling water valve
Sticking
IDV(
15)
Condenser cooling water valve
Sticking
IDV(
16)
Unknown
Unknown
IDV(
17)
Unknown
Unknown
IDV(
18)
Unknown
Unknown
IDV(
19)
Unknown
Unknown
IDV(
20)
Unknown
UnknownIn our simulation, these data are acquired in the case of Mode 1. The operating time is set to be 36 h. The total number of samples isN=3601. The faults are introduced at 8 h, which means that the fault occurs after 800 samples. p and f are both set as 10.


### 4.2. Simulation Results

The process is simulated with 20 different faults. For the sake of simplicity, we only take some typical faults to show. The first type of fault is step change inA and C feed ratio (IDV1). Figure 3 shows the testing statistic J based on residual vectors for input feed block. The blue line is the testing statistic J and the threshold Jth is shown by a red line. The second type of fault is a random variation change in A, B, and C feed composition (IDV8). The simulated result is shown in Figure 4. Figure 5 shows the testing statistic of the third type of fault which is a sticking fault in reactor cooling water valve (IDV14). The last fault (IDV17) shown in Figure 6 is an unknown fault. The type and process variables of fault are not known.Figure 3
Detection ofIDV(1).Figure 4
Detection ofIDV(8).Figure 5
Detection ofIDV(14).Figure 6
Detection ofIDV(17).Based on the fault detection method, if the testing statistics exceed the threshold, the fault is detected. We can see from Figures3–6 that all the testing statistics J exceed the threshold Jth after 800 samples. The four faults have been detected.In addition to the simulation results, the FDR, FAR, and FDT shown in Table4 are also obtained to demonstrate the efficiency. Besides the comparison between SMI-PCA and SMI-PCA KDE, we also compare them with the classic PCA in order to better reflect the advantages of this proposed method. According to the study of Mahadevan and Shah [20], the FDR of T2 and Q for PCA are listed in the table. The maximum fault detection rate has been highlighted in boldface. We can see that most of the fault detection rate of our proposed method is higher than PCA. The fault detection rate of using kernel density estimation to calculate the threshold is slightly higher than the chi-squared distribution. Most of the faults can be detected well and they have high FDRS and low FARS. But for some faults, that is, IDV(3), (5), (9), and (15-16), the efficiency of fault detection is very poor. This may be caused by a very small change in the variables, which is also a common problem in the TE process fault detection. There is no meaning to calculate the fault detection time if the faults cannot be detected, so the fault detection time is indicated with “∗” in the table. The fault detection time of SMI-PCA and SMI-PCA KDE is the same because of their similar detection rate.Table 4
TE process: fault detection for test data (%).
Fault
PCA
SMI-PCA
SMI-PCA KDE
FDR
(
T
2
)
FDR
(
Q
)
FDR
FAR
FDT
FDR
FAR
FDT
IDV
(
1
)
99.2
99.8
100
4.25
800
100
5.63
800
IDV
(
2
)
98.0
98.6
99.46
3.88
800
99.5
5.25
800
IDV
(
3
)
0
0
11.97
3.75
∗
14.52
5.12
∗
IDV
(
4
)
4.4
96.2
100
4.62
800
100
6
800
IDV
(
5
)
22.5
25.4
14.92
4.37
∗
17.43
5.75
∗
IDV
(
6
)
98.9
100
100
4.62
800
100
6.12
800
IDV
(
7
)
91.5
100
100
4.75
800
100
6.12
800
IDV
(
8
)
96.6
97.6
97.70
3.62
811
97.81
5.12
811
IDV
(
9
)
0
0
12.4
3.62
∗
14.49
5.12
∗
IDV
(
10
)
33.4
34.1
81.56
3.62
811
82.82
5.12
811
IDV
(
11
)
20.6
64.4
98.67
3.62
811
98.78
5.12
811
IDV
(
12
)
97.1
97.5
43.35
3.62
811
46.76
5.12
811
IDV
(
13
)
94.0
95.5
98.81
3.62
811
98.85
5.12
811
IDV
(
14
)
84.2
100
99.82
3.62
801
99.86
5
801
IDV
(
15
)
0
0
5.18
3.62
∗
6.25
5.12
∗
IDV
(
16
)
16.6
24.5
4.46
3.62
∗
5.9
5.12
∗
IDV
(
17
)
74.1
89.2
91.41
3.62
811
91.88
5.12
811
IDV
(
18
)
88.7
89.9
30.63
3.62
811
32.49
5.12
811
IDV
(
19
)
0.4
12.7
99.75
3.75
801
99.75
5.12
801
IDV
(
20
)
26.4
43
89.25
3.62
811
90.08
5.12
811


## 4.1. TE Process

TE process model is a realistic simulation program of a chemical process that is widely accepted as a benchmark for control monitoring studies [18]. The TE model that we used was downloaded from http://depts.washington.edu/control/LARRY/TE. The flow diagram of the process is shown in Figure 2. The process contains five major units named reactor, condenser, compressor, separator, and stripper.Figure 2
The Tennessee Eastman process.There are 53 measurements, 12 of which are manipulated variables and 41 are process variables. The 12 manipulated variables are listed in Table1. In this paper, the other 9 manipulated variables except for XMV(5), XMV(9), and XMV(12) are taken as inputs. The dimension of the TE dataset is very large. To overcome the curse of dimensionality, the 41 process variables are divided into eight blocks [19]. We choose the input feed block as outputs, which are shown in Table 2. Table 3 shows the 20 process faults of TE process.Table 1
Process manipulated variables.
Variable name
Variable number
D feed flow
XMV
(
1
)
E feed flow
XMV
(
2
)
A feed flow
XMV
(
3
)
A and C feed flow
XMV
(
4
)
Compressor recycle valve
XMV
(
5
)
Purge valve
XMV
(
6
)
Separator pot liquid flow
XMV
(
7
)
Stripper liquid product flow
XMV
(
8
)
Stripper steam valve
XMV
(
9
)
Reactor cooling water flow
XMV
(
10
)
Condenser cooling water flow
XMV
(
11
)
Agitator speed
XMV
(
12
)Table 2
Input feed block process measurements.
Variable name
Variable number
A feed (stream1)
XMEAS
(
1
)
D feed (stream1)
XMEAS
(
2
)
E feed (stream1)
XMEAS
(
3
)
A and C feed
XMEAS
(
4
)Table 3
Process faults.
Variable number
Process variable
Type
IDV(
0)
Normal operation
Step
IDV(
1)
A
/
C feed ratio, B composition constant
Step
IDV(
2)
B composition, A/C ratio constant
Step
IDV(
3)
D feed temperature
Step
IDV(
4)
Reactor cooling water inlet temperature
Step
IDV(
5)
Condenser cooling water inlet temperature
Step
IDV(
6)
A feed loss
Step
IDV(
7)
C header pressure loss-reduced availability
Step
IDV(
8)
A
,
B
,
C feed composition
Random variation
IDV(
9)
D feed temperature
Random variation
IDV(
10)
C feed temperature
Random variation
IDV(
11)
Reactor cooling water inlet temperature
Random variation
IDV(
12)
Condenser cooling water inlet temperature
Random variation
IDV(
13)
Reaction kinetics
Random variation
IDV(
14)
Reactor cooling water valve
Sticking
IDV(
15)
Condenser cooling water valve
Sticking
IDV(
16)
Unknown
Unknown
IDV(
17)
Unknown
Unknown
IDV(
18)
Unknown
Unknown
IDV(
19)
Unknown
Unknown
IDV(
20)
Unknown
UnknownIn our simulation, these data are acquired in the case of Mode 1. The operating time is set to be 36 h. The total number of samples isN=3601. The faults are introduced at 8 h, which means that the fault occurs after 800 samples. p and f are both set as 10.


## 4.2. Simulation Results

The process is simulated with 20 different faults. For the sake of simplicity, we only take some typical faults to show. The first type of fault is step change inA and C feed ratio (IDV1). Figure 3 shows the testing statistic J based on residual vectors for input feed block. The blue line is the testing statistic J and the threshold Jth is shown by a red line. The second type of fault is a random variation change in A, B, and C feed composition (IDV8). The simulated result is shown in Figure 4. Figure 5 shows the testing statistic of the third type of fault which is a sticking fault in reactor cooling water valve (IDV14). The last fault (IDV17) shown in Figure 6 is an unknown fault. The type and process variables of fault are not known.Figure 3
Detection ofIDV(1).Figure 4
Detection ofIDV(8).Figure 5
Detection ofIDV(14).Figure 6
Detection ofIDV(17).Based on the fault detection method, if the testing statistics exceed the threshold, the fault is detected. We can see from Figures3–6 that all the testing statistics J exceed the threshold Jth after 800 samples. The four faults have been detected.In addition to the simulation results, the FDR, FAR, and FDT shown in Table4 are also obtained to demonstrate the efficiency. Besides the comparison between SMI-PCA and SMI-PCA KDE, we also compare them with the classic PCA in order to better reflect the advantages of this proposed method. According to the study of Mahadevan and Shah [20], the FDR of T2 and Q for PCA are listed in the table. The maximum fault detection rate has been highlighted in boldface. We can see that most of the fault detection rate of our proposed method is higher than PCA. The fault detection rate of using kernel density estimation to calculate the threshold is slightly higher than the chi-squared distribution. Most of the faults can be detected well and they have high FDRS and low FARS. But for some faults, that is, IDV(3), (5), (9), and (15-16), the efficiency of fault detection is very poor. This may be caused by a very small change in the variables, which is also a common problem in the TE process fault detection. There is no meaning to calculate the fault detection time if the faults cannot be detected, so the fault detection time is indicated with “∗” in the table. The fault detection time of SMI-PCA and SMI-PCA KDE is the same because of their similar detection rate.Table 4
TE process: fault detection for test data (%).
Fault
PCA
SMI-PCA
SMI-PCA KDE
FDR
(
T
2
)
FDR
(
Q
)
FDR
FAR
FDT
FDR
FAR
FDT
IDV
(
1
)
99.2
99.8
100
4.25
800
100
5.63
800
IDV
(
2
)
98.0
98.6
99.46
3.88
800
99.5
5.25
800
IDV
(
3
)
0
0
11.97
3.75
∗
14.52
5.12
∗
IDV
(
4
)
4.4
96.2
100
4.62
800
100
6
800
IDV
(
5
)
22.5
25.4
14.92
4.37
∗
17.43
5.75
∗
IDV
(
6
)
98.9
100
100
4.62
800
100
6.12
800
IDV
(
7
)
91.5
100
100
4.75
800
100
6.12
800
IDV
(
8
)
96.6
97.6
97.70
3.62
811
97.81
5.12
811
IDV
(
9
)
0
0
12.4
3.62
∗
14.49
5.12
∗
IDV
(
10
)
33.4
34.1
81.56
3.62
811
82.82
5.12
811
IDV
(
11
)
20.6
64.4
98.67
3.62
811
98.78
5.12
811
IDV
(
12
)
97.1
97.5
43.35
3.62
811
46.76
5.12
811
IDV
(
13
)
94.0
95.5
98.81
3.62
811
98.85
5.12
811
IDV
(
14
)
84.2
100
99.82
3.62
801
99.86
5
801
IDV
(
15
)
0
0
5.18
3.62
∗
6.25
5.12
∗
IDV
(
16
)
16.6
24.5
4.46
3.62
∗
5.9
5.12
∗
IDV
(
17
)
74.1
89.2
91.41
3.62
811
91.88
5.12
811
IDV
(
18
)
88.7
89.9
30.63
3.62
811
32.49
5.12
811
IDV
(
19
)
0.4
12.7
99.75
3.75
801
99.75
5.12
801
IDV
(
20
)
26.4
43
89.25
3.62
811
90.08
5.12
811


## 5. Conclusions

In this paper, subspace method aided data-driven fault detection based on PCA has been presented. This method is to identify the parameterized matrices using PCA and construct residual generators with the input and output data. A simulation study on the TE process demonstrates the availability of this method. It indicates that the method proposed in this paper has better effects than PCA on fault detection rate and is suitable for linear systems which are observable and controllable. Moreover, the problem of threshold is also discussed in this paper. The fault detection rate of using kernel density estimation is slightly superior to the chi-squared distribution.Compared with the traditional model-based fault detection method, prior knowledge of model mechanisms is not needed and only the parameterized matrices related to residuals rather than the system models need to be identified.



---

*Source: 1812989-2017-11-27.xml*