Recursive Fault Detection and Identification for Time-Varying

Article pubs.acs.org/IECR

Recursive Fault Detection and Identification for Time-Varying Processes Liangliang Shang,†,‡ Jianchang Liu,†,§ and Yingwei Zhang*,†,§ †

College of Information Science and Engineering and §State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University, Shenyang, Liaoning 110819, PR China ‡ School of Electrical Engineering, Nantong University, Nantong, Jiangsu 226019, PR China

ABSTRACT: Canonical variate analysis (CVA) has been extensively applied in monitoring of different industrial processes. However, conventional CVA is unable to handle the characteristics of time-varying processes. It tends to interpret the natural changes of the process as faults, which would cause high false alarm rates. To solve this problem, a recursive canonical variate analysis based on the first order perturbation theory (RCVA-FOP) is proposed to detect faults in time-varying processes. Without recalling past training data, the covariance of past observation vectors is updated by the exponential weighted moving average (EWMA) method. Moreover, the first order perturbation theory is introduced to realize the recursive singular value decomposition (SVD) of the Hankel matrix, which can reduce computational time significantly compared with the conventional SVD. To identify the real reason for a fault, an EWMA contribution plot based on CVA is also proposed to enhance the fault identification rate. The proposed method is verified with simulations of the continuous stirred tank reactor. Simulation results indicate that the RCVA-FOP method not only can effectively adapt to the natural changes of time-varying processes but also can detect and identify three types of faults, which include sensor precision degradation, heat exchanger fouling fault, and sensor bias. identification (N4SID)5 are the three conventional subspace identification methods. Juricek et al.6 certified that the accuracy of CVA was higher than that of N4SID. Larimore7 clarified that only the CVA procedure has been developed on the basis of optimal statistical inference principles, and as a result it achieves optimal statistical accuracy while the others can be considerably less accurate. Many successful examples for fault detection and identification with the CVA approach have already been reported. Negiz and Cinar8 applied the canonical variate state space (CVSS) model in milk pasteurization process monitoring.

1. INTRODUCTION As the competition of the global market becomes increasingly intense, modern industry has to be more efficient and raise the quality of end products. As an essential technique, process monitoring has already become more and more important in recent years. The main process monitoring approaches consist of three categories: the knowledge-based approach, the datadriven approach, and the model-based approach.1 The disadvantage of the knowledge-based approach is that it needs a large amount of knowledge, which is not easy to obtain in practice. Although the data-driven approach has been widely applied, subspace identification methods (SIMs) have greatly attracted researchers’ attention for monitoring and modeling in the last two decades.2−4 Canonical variate analysis (CVA), multivariable output-error state space (MOESP), and numerical algorithms for subspace state space system © XXXX American Chemical Society

Received: July 12, 2016 Revised: September 16, 2016 Accepted: October 21, 2016

A

DOI: 10.1021/acs.iecr.6b02653 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Russell and Braatz et al.9 implemented CVA in the Tennessee Eastman (TE) process simulator for fault detection. Juricek et al.2 applied CVA to monitor process faults in a continuous stirred tank reactor. Stubbs et al.10 applied CVA to the Tennessee Eastman process simulator for fault detection and diagnosis. CVA shows good monitoring performance in stationary operating conditions and linear processes. However, most industrial processes have the characteristics of timevarying, such as input feed to a chemical reactor, slow drift in process parameters, and so on.11−13 Thus, it is necessary to develop an online fault detection and identification approach based on recursive canonical variate analysis for time-varying processes. One of the most important problems for recursive modeling identification is how to reduce the complexity of computation. According to different ways of getting the state variables or the observation matrices, the recursive subspace identification methods can be classified into two categories. One category is to avoid the singular value decomposition, and the other is to reduce the cost of SVD calculation or just directly use it. For the first category, Gustafusson14 originally developed a recursive subspace identification method, in which projection approximation subspace tracking (PAST) was adopted to estimate the extended observability matrices recursively. Oku and Kimura15 proposed recursive 4SID algorithms based on gradient type subspace tracking. Mercère and Lecoeuche et al.3 proposed two recursive subspace identification algorithms, in which the observability subspace can be estimated via a specific unconstrained optimization method. Li and Chu et al.16 proposed a recursive closed-loop subspace identification approach, which is based on a predictor-based subspace identification optimized version (PBSIDopt) and projection approximation subspace tracking (PAST) algorithm. Houtzager and Verhaegen et al.4 proposed a predictor-based recursive subspace identification method in either open-loop or in closed-loop. For the second category, Lee et al.17 proposed a recursive model updating method based on CVA for monitoring the Tennessee Eastman Process. Ding et al.18 presented a fault detection and isolation approach based on the data-driven parity space and observer. Naik et al.12 proposed two recursive algorithms for fault detection, which are based on the perturbation theory and the orthogonal iterations. Pan and Sheng et al.19 proposed an online virtual metrology approach based on recursive CVA with application to an industrial sputtering process. In addition, a few monitoring approaches based on the multivariate statistical theory have also been proposed. Choi and Lee et al.20 described an adaptive monitoring approach based on multivariate statistical process monitoring (MSPC) for a operating condition changes process. Elshenawy et al.21,22 presented two recursive principal component analysis (RPCA) algorithms, in which the data projection method (DPM) and the first order perturbation analysis (FOP) were adopted respectively; and they also proposed two recursive fault isolation methods for monitoring time-varying processes, in which two types of natural changes represent the time-varying characteristics of a process. After one or more monitoring statistics exceed their corresponding control limits, a fault will be detected. Subsequently, the next essential procedure needed to be done is to identify the real reason for this fault. Although many works about fault detection have been reported, only a few papers on fault identification are available. The contribution plot is the

most popular method for the Q and T2 statistics, which is based on an assumption that faulty variables cause high contributions to monitoring statistics.22 Jiang and Braatz et al.23 proposed two contribution charts based on the canonical state vectors and the corresponding residual vectors for offline fault identification. The other existing methods can be divided into the following four categories: (1) reconstructed-based indices,24 which identify faults based on whether a reconstruction can bring the process return to normal; (2) discrimination by angles,25 which extracts fault signatures compared with the known fault directions; (3) pattern matching methods-based similarity factors,26,27 which perform diagnosis by calculating similarity factors or dissimilarity factors between current snapshot data and historical fault data; (4) isolation enhanced techniques,28 in which sensor and actuator faults can be diagnosed by directional and structured residuals from model-based methods. In general, the fault identification methods are more suitable for time-invariant processes. Therefore, an efficient online faults identification method for time-varying processes is required. The main contributions of this paper are 2-fold: First, a recursive CVA based on the first order perturbation theory (RCVA-FOP) is proposed for monitoring time-varying processes. Second, an exponential weighted moving average contribution plot based on CVA is also proposed for the faults online identification. The proposed recursive CVA utilized the first order perturbation theory to realize recursive SVD on the scaled Hankel matrix. Compared with the conventional CVA, the computation cost of recursive SVD is significantly less. Three types of faults from the continuous stirred tank reactor, including sensor precision degradation, heat exchanger fouling fault, and sensor bias, are simulated to verify the proposed approach. The remainder of the paper is outlined as follows. Section 2 reviews the conventional canonical variate analysis. Section 3 presents RCVA-FOP for monitoring time-varying processes. Section 4 provides the statistics and the control limits for process monitoring, and the fault identification method based on the exponential weighted moving average contribution plot. Section 5 applies the proposed method to a continuous stirred tank reactor for simulation, illustrates the monitoring performance of the RCVA-FOP approach for time-varying processes, and provides a comparative analysis of results with those of CVA, and also plots the EWMA contribution for fault identification. Conclusions are presented in Section 6.

2. CANONICAL VARIATE ANALYSIS The multivariate stochastic state space model could be formulated as follows: xk + 1 = Axk + wk yk = Cxk + vk

(1)

where xk ∈ n and yk ∈ m are state vectors and observation

vectors respectively at the sampling time k; A ∈ n × n and C ∈ m × n are the system and output matrix, respectively; and wk ∈ n is the stochastic disturbance; vk ∈ m is the measurement noise, which is assumed to be zero-mean Gaussian white noise. The past and future observation vectors yp,r and yf,r can be expanded by using the past and future measurements, respectively, B


Article

Industrial & Engineering Chemistry Research

yp , r

yf , r

⎡ yr − 1 ⎤ ⎢ ⎥ ⎢ yr − 2 ⎥ mp =⎢ ⎥∈ ⋮ ⎢ ⎥ ⎢⎣ yr − p ⎥⎦

(2)

⎡ yr ⎤ ⎢ ⎥ ⎢ yr + 1 ⎥ mp =⎢ ⎥∈ ⋮ ⎢ ⎥ ⎢⎣ yr + p − 1⎥⎦

(3)

−1/2 where J = VTΣpp is the transformation matrix, which transforms the past measurements into canonical variate space. It consists of the main space and the residual space, and the dominant singular values can be used to determine the dimension of the main space. The first n dominant singular values determine the dimension of the main space xr ∈ n , and the remaining singular values determine the dimension of the residual space dr ∈ mp − n, where n < mp. The number of singular values n should be known a prior. Therefore, the canonical variate can be written as follows

sr = Jyp , r = [[Jx yp , r ]T drT]T = [xrT drT]T

where r denotes a generic index and p represents the window length of the past observation vectors. Following the terminologies and equations in ref 29, the past and future Hankel matrices Yp and Yf are defined as Yp = [ yp , p + 1 yp , p + 2 ... yp , p + N ] ∈ mp × N

(4)

Yf = [ yf , p + 1 yf , p + 2 ... yf , p + N ] ∈ 

(5)

mp × N

The state vectors are not only part of the estimated canonical variates, but also are a linear combination based on the past observation vectors. Hence, the state vectors can be defined as a linear combination based on the past observation vectors. xr = Jx yp , r = VTx Σ−pp1/2yp , r

Σpp =

− 1)

(6)

Σff = Yf YTf /(N − 1)

(7)

Σfp = Yf YTp /(N − 1)

(8)

er = (I − VxVTx )Σ−pp1/2yp , r

a TΣfpb (a TΣff a)1/2 (bTΣppb)1/2

(9)

1/2 By define u1 = Σ1/2 f f a and u2 = Σpp b, the optimization to determine u1 and u2 can be represented as

max u1, u 2

s. t .

∑pp , r = (1 − β)yp , r ypT, r + β ∑pp , r − 1

u1T(Σ−ff1/2ΣfpΣ−pp1/2)u 2 u1Tu1 = 1; u 2Tu 2 = 1

(10)

(11)

The canonical variate sr can be obtained by performing the SVD on eq 11 as follows: ⎡bT ⎤ ⎢ 1 ⎥ ⎢bT ⎥ 2 sr = ⎢ ⎥yp , r = VTΣ−pp1/2yp , r = Jyp , r ⎢⋮⎥ ⎢ ⎥ T ⎢⎣bmp ⎥⎦

(16)

where β is the forgetting factor. r indicates the number of every iteration in the recursive case. It starts from the sampling time instant 2p. In other words, we start performing online recursive identification after we collect measurements at sampling time 2p. 3.1. Update state vectors based on the first order perturbation theory. In this subsection, we briefly clarify to update the state vectors based on the perturbation theory recursively, which is mainly adopted in signal subspace tracking applications. We consider a multidimensional signal zr ∈ q and a total of M samples of it:

The optimization problem can be obtained via performing SVD of the scaled Hankel matrix H. H = Σ−ff1/2ΣfpΣ−pp1/2 = UΣVT

(15)

3. RECURSIVE CANONICAL VARIATE ANALYSIS BASED ON THE FIRST ORDER PERTURBATION THEORY Most industrial processes have the characteristics of timevarying, reflected in the change of the variance, mean, and correlation, which always cause change of the order of state vectors. When new measurements become available at each sampling time, both the mean and covariance should be updated as the degree of variation in the model structure. In this paper, the exponential weighted moving average11,29 is adopted to realize a recursive canonical variate analysis approach, in which there is no need to recall the past training data. The covariance and cross-covariance matrices can be estimated respectively as follows:

CVA tries to obtain the best linear combinations aTyf,r, a ∈ mp, and bTyp,r, b ∈ mp, so that the correlation ρfp(a, b) between the linear combinations is maximized. ρfp (a , b) =

(14)

where Vx consists of the first n columns of V. The prediction errors are represented in the scaled past observation space as follows.

Given a data set with M observations of m measured variables, the last element of yp,p+1 in eq 2 is y1 and the last element of yf,p+N in eq 3 would be yl. The dimension of the column in these Hankel matrices is N = M − 2p + 1. Subsequently, the cross-covariance and covariance matrices can be calculated by using YpYTp /(N

(13)

⎡ zr ,1 zr + 1,1 ⎢ ⎢ zr ,2 zr + 1,2 Zr = ⎢ ⋮ ⎢ ⋮ ⎢⎣ zr , q zr + 1, q

(12) C

⋯ zr + M − 1,1 ⎤ ⎥ ⋮ zr + M − 1,2 ⎥ q×M ⎥∈ ⋱ ⋮ ⎥ ⋯ zr + M − 1, q ⎥⎦

(17)


Article

Industrial & Engineering Chemistry Research M

Similarly, we define Zτ ∈ q × M as instrument variables, which Zτ can be written as ⎡ zr − τ ,1 zr − τ + 1,1 ⎢ ⎢ zr − τ ,2 zr − τ + 1,2 Zτ = ⎢ ⋮ ⎢ ⋮ ⎢⎣ zr − τ , q zr − τ + 1, q

⋯ zr − τ + M ,1 ⎤ ⎥ ⋮ zr − τ + M ,2 ⎥ q×M ⎥∈ ⋱ ⋮ ⎥ ⋯ zr − τ + M , q ⎥⎦

Ur̅ , i = Ur − 1, i +

q×q

j=1

Vr̅ , i = Vr − 1, i +

j=1

q×q

where Ur ∈  and Vr ∈  are matrices including left singular vectors and right singular vectors, and Σ = diag(σ1, σ1, ..., σq). Equation 21 can be reformulated in the following recursive manner Φz , r = Φz , r − 1 + E = Ur − 1(Σr − 1 +

dji ≈

σr − 1, i fijT σr − 1, i

= (1 − β)

= (1 − β)

(23)

∑ i−1

qiu =

∑ j=1

zτ̅ , jVr − 1, j

j=i+1

zr̅ , j σr − 1, j

Vr − 1, j (29)

zτ̅ , jUr − 1, j = piu− 1 − zτ̅ , iUr − 1, i zr̅ , j

σr − 1, j

Ur − 1, j = qiu− 1 +

zr̅ , i σr − 1, i

(30)

Ur − 1, i (31)

1 − β zr ∈ q and q0u = [0 ··· 0] ∈ q . Sim-

where p0u = ilarly, M

piv =

zτ̅ , jVr − 1, j = piv− 1 − zτ̅ , iVr − 1, i

∑ j=i+1 i−1

qiv =

∑ j=1

zr̅ , j σr − 1, j

Vr − 1, j = qiv− 1 +

zr̅ , i σr − 1, i

(32)

Vr − 1, i (33)

where p0v = 1 − β zr ∈ q and q0v = [0 ··· 0] ∈ q . Then, the left singular vectors can be updated as the following equations: piu = piu− 1 − zr̅ , iUr − 1, i

(34)

qiu = qiu− 1 + σr−−11, izτ̅ , iUr − 1, i

(35)

Ur̅ , i = Ur − 1, i + σr−−11, izτ̅ , ipiu − zr̅ , iqiu

(36)

The right singular vectors can be calculated as the following equations:

(24)

VrT− 1, jzτzrTUr − 1, i σr − 1, i

M

∑

j=i+1

UrT− 1, jzrzτTVr − 1, i σr − 1, i

(28)

M

piu =

where F ∈  can be directly computed F = UTr−1EVr−1. q×q E ∈  is the first order perturbation matrix. The result in ref 30 presented the origin of updating the singular values and singular vectors based on perturbation theory. The left and right singular vectors at sample time r can be formulated in terms of small perturbations of those at r − 1, U̅ r = Ur−1(I + B), and V̅ r = Vr−1(I + D), where I is the identity matrix. The overbar indicates the singular vectors are unnormalized; thus, the ith unnormalized left and right singular vectors can be denoted U̅ r,i and V̅ r,i, respectively. In our case study, offline identification offers the initial value of singular values and vectors. The state vectors can be recursively updated when the current measurements are available. We denote zr by zr = (yf,ryTf,r)yf,r, where yf,r is column vectors assembled with future measurements. The instrument variable vector zτ can be similarly established with the past measurements zτ = (yp,ryTp,r) yp,r. According to the results in ref 31, assuming the singular values are simple, their (i, j) elements of B and D, bji = UTr−1,iUr,j and dij = VTr−1,iVr,j, can be approximated as j < i, bji = −bTij , dji = −dTij ; i = j, bii = 0, dii = 0; for i < j, bji ≈

Ur − 1, j

The series of transition vectors can be defined as follow,

q×q

f ji

zr̅ , jUr − 1, j

j=i+1

σr − 1, j

σr − 1, i

i−1

M

∑

zτ̅ , j

zr̅ T, i

− zτ̅ T, i ∑

(22)

F)VrT− 1

σr − 1, i

i−1

(21)

Vr VrT = I

(27)

zτ̅ T, i

− zr̅ T, i ∑

where β is a forgetting factor with the range from 0 to 1. It is worth noting that Φz is a positive semidefinite matrix. The SVD of Φz,r then gives

UrUrT = I ;

djiVr − 1, j

Let us define z̅ r ,i and z̅ τ ,i as the ith elements of zr̅ = 1 − β UrT− 1zr ∈ q and zτ̅ = 1 − β VrT− 1zτ ∈ q , respectively. The recursive updating of singular vectors can be efficiently implemented as follows,

(18)

(20)

Φz , r = UrΣr VrT

∑ j = 1, j ≠ i

In a dynamic environment, it can be recursively calculated as follows Φz , r = β Φz , r − 1 + (1 −

(26)

M

Vr̅ , i = Vr − 1, i +

(19)

β)ZrZTτ

bjiUr − 1, j

j = 1, j ≠ i

We denote the product of Zr and Zτ for M samples by Φz ∈ q × q, and it can be computed as

Φz = ZrZTτ /M

∑

Ur̅ , i = Ur − 1, i +

(25)

The left and right singular vectors can then be updated using D

piv = piv− 1 − zτ̅ , iVr − 1, i

(37)

qiv = qiv− 1 + σr−−11, izr̅ , iVr − 1, i

(38)

Vr̅ , i = Vr − 1, i + σr−−11, izr̅ , ipiv − zτ̅ , iqiv

(39)


Article


cumulative percent variance (CPV), which is the percent variance captured by the first n singular values.

At last, the singular values can be recursively updated as Σ̅ r , i = β Σr − 1, i + (1 − β)UrT− 1, izrzτVrT− 1, i

(40)

n

The dimension of the significant space can be obtained in eq 40. The singular and singular vectors in eqs 36, 39, and 40 are only scaled, which can be normalized as follows Ur , i =

Ur̅ , i Ur̅ , i

, Σr , i = |Σ̅ r , i| , Vr , i =

CPV (n) =

(41)

Given updated the scaled singular vectors, the new canonical variates can be obtained. The order of the recursive model should be known in advance. The equations for computing the state vector and residual vector are given as xd , r = Jd yp , r = VTx Σ−pp1/2yp , r

(42)

ed , r = Je yp , r = (I − VxVTx )Σ−pp1/2yp , r

(43)

mp

∑ j = 1 σj

× 100% (44)

4. FAULT DETECTION AND IDENTIFICATION BASED ON RCVA-FOP 4.1. Statistics and control limits of process monitoring. A few monitoring statistics based on CVA have been already proposed. Lee et al.35 proposed the state and total process noise monitoring statistics that follow an F-distribution, and the corresponding control limits can be obtained by using Hotelling’s T2 statistics. Odiowei and Cao33 proposed T2 based on the state vectors and Q metrics based on the residual variables from the state-space model updated by CVA. More appropriate control limits for nonlinear dynamic process monitoring were directly estimated by the underlying probability density function (PDF) of the monitoring metrics. Stubbs et al.10 proposed a simpler state space model based on CVA for the specific process monitoring. The first k CVA states T2 statistics used by Negiz and Cinar8 and T2 and Q statistics based on the state residuals and output residuals matrix proposed by Simoglou et al.38 were employed for the benchmark Tennessee Eastman process monitoring. Two monitoring statistics from the state and noise spaces are respectively considered to monitor the abnormalities of the process. We define the T2 monitoring index based on the state spaces proposed by33

Vr̅ , i Vr̅ , i

∑ j = 1 σj

where xd,r and ed,r are the state vector and residual vector at the time point r, respectively; yp,r is the new past observation vector at the time point r in the recursive modeling process. The complexity and computation cost are only calculated for complete update of the decomposition of Φz ∈ q × q, if l denotes the rank of the subspace estimated, since usually l ≪ q. Many alternative approaches requiring fewer operations were proposed to avoid the excessively high computational cost. The complexity of several common algorithms is summarized in Table 1. Table 1. Comparison of Computation Cost

T 2 = xdT, rxd , r

Method

Flops

Complexity

SVD Inverse iteration Rank-1 update PAST RSVD

22q3 + 4q2 2/3q3l + 2q2l + 5ql + l 2q3 + 9q2 + 4q + 1 7/2q2 + 6ql + 4l2 14q2 + 7q

O(q3) O(q3l) O(q3) O(q2) O(q2)

(45)

The corresponding upper control limit for a significance level is derived as follows: 2 (α ) = TUCL

n(N 2 − 1) Fn , N − n(α) N (N − n)

(46)

where Fn,N−n(α) is the critical value of the F-distribution with n and N − n degrees of freedom for a significance level α. By comparing T2 with T2UCL(α) in real-time, an abnormal condition can be determined when T2 > T2UCL(α). The total process noise monitoring index can be represented as follows:

Compared to the SVD, the computation cost of recursive SVD based on the first order perturbation theory is significantly less. 3.2. Recursive determinate of the dimension of main space. It is indispensable to determine adaptively the number of significant singular values because it may vary over time. There are numerous proposed methods for determining the order of the CVA model. The percent of cumulative singular values to the sum of all singular values is employed to select the order of the system.32 The number of dominant singular values which represent the dimension of the canonical variate space were assigned as the order of the system.33,34 The number of state order was selected, in which over 95% of the total variance of the scaled Hankel matrix was explained.17,35 The Akaike Information Criterion (AIC), that is an informational theoretical approach, is one of the popular methods.36 Because CVA is a near maximum-likelihood (ML) procedure, likelihood ratio tests, that is the ML-based statistical selection method, can be adopted to select the state order.2 The optimal dimension n is chosen as the one which gives the best one-step-ahead, crossvalidated predictions.37 However, all the approaches are not necessarily suitable for recursive CVA. For example, because previous data cannot represent the current process, the crossvalidation approach is not suitable. In this paper, we adopt the

Q = edT, red , r

(47)

where ed,r is the total process noise, which has a normal distribution. Under the assumption of normal distribution, the threshold of the Q monitoring statistics can be estimated by ref 33 as follows 1/ h0 ⎛ h c 2θ θ2h0(h0 − 1) ⎞ 0 α 2 ⎟⎟ Q UCL(α) = θ1⎜⎜ +1+ θ1 θ12 ⎝ ⎠

(48)

where θ i = hr − n, hr is the rank of scaled Hankel matrix H, h0 = 1 − (2θ1θ3/3θ22), and cα is the normal deviate corresponding to 1 − α percentile. An abnormal condition is determined in real time when Q > QUCL(α). As a matter of fact, a fault is then detected if either the T2 > T2UCL(α) or Q > QUCL(α) condition is satisfied. 4.2. Fault identification based on the EWMA contribution. After detecting a fault, the most important E


Article

Industrial & Engineering Chemistry Research work is to find the real reason for the fault. In this section, T2 and Q monitoring statistics are adopted for computing contributions. The value of two statistics can be severely affected by deviated variables. Although state vectors do not directly identify the variables, the loadings (Jd, Fd) indicate how much of the contribution projected to the state space and the residual space. Considering the contribution based on the state and residual vectors from the recursive CVA model, a recursive weighted variable contribution chart for online fault identification was proposed. Motivated by the adaptive contribution of T2 in ref 22, the contribution based on the state vector for a new sample time vector yp,r can be defined as follows,

and each variable number are the horizontal and vertical coordinates, respectively. The contribution is shown as different degrees of gray. Much more information about 2D contribution maps can be found in refs 23 and 39. The proposed approach includes two stages: (1) off-line training and monitoring model and (2) online monitoring and fault identification. The detailed steps are given as follows: Off-Line Training and Monitoring. 1. Set the value of p which satisfies (n/m) < p, and obtain a data set Y with M observations. The order of system n should be given previously. 2. Assemble matrices Yp and Yf using eqs 4 and 5. 3. Covariance and cross-covariance matrices can be computed by using eqs 6−8. 4. Perform SVD on the Hankel matrix H in eq 11. 5. Compute the state vector and residual vector using eqs 14 and 15. 6. Compute T2, Q statistics and the corresponding control limits using eqs 45−48. If either the T2 > T2UCL(α) or Q > QUCL(α) condition is satisfied, a fault happens. Online Monitoring and Fault Identification. 1. Obtain the measurements at the sample time 2p, and normalize. 2. Calculate T2, Q monitoring statistics and the corresponding control limits. 3. Judge whether the slow drift of the process is in normal or not. If the monitoring statistics do violate the control limits, the procedure can be continually performed as follows. (a) assemble matrices Yp and Yf using eqs 4 and 5. (b) Update recursively covariance and cross-covariance using eq 16. (c) Update recursively singular vectors and singular values of the Hankel matrix using eqs 34−41. (d) Update the state vectors and residual vectors using eqs 42 and 43. (e) Go to step 1. Otherwise, if the monitoring statistics overstep their control limits, go to the next step. 4. A fault is detected if either T2 > T2UCL(α) or Q > QUCL(α) condition is satisfied. 5. After a fault is detected, the proposed EWMA contribution approach starts to identify the true faulty variable. (a) Calculate the combined contribution of all the variables using eq 53. (b) Calculate the EWMA contribution of all the variables using eq 54. (c) Check which process variable has the largest contribution.

mp

C(dr) = xdT, rxd , r = xdT, r Jd yp , r = xdT, r Jd

∑ ξiξiTyp,r i=1

mp

mp

=∑ xdT, r Jd ξiξiTyp , r = i=1

∑ xdT, r Jid ypi , r i=1

(49)

where yip,r is the ith element of the data vector yp,r, Jid is the ith row of the Jd, and ξi is the ith column of the identity matrix. Due to serial correlation contained in CVA, the contribution for variable ym in the state space is the sum of all its p lagged observations ym(r−j), (j = 1, 2, ..., p). p

C yd (r ) = m

∑ xdT,r Jd ,m ym(r− j) j

(50)

j=1

In a similar way, the process variable contribution based on the residual space can be defined as follows p

C ye (r ) = m

∑ edT,r Je ,m ym(r− j) j

(51)

j=1

A process variable that has higher contribution indicates a more severe fault with the specific variable. Because of the inversion 23 of Σ−1/2 contributions are excessively sensitive. Combining pp , the two contributions is conducive to reduce the sensitivity for fault identification. The combined contribution used in our case study is as follows, Ccom(r ) = (C yd (r ) + C ye (r ))/2 m

m

(52)

Considering the effect of serial correlation from the data and random noise, variables at the lagged sample times have different contributions. To improve the correct fault identification rate, an exponential weighted moving average (EWMA) approach is introduced to combine with the combined variable contributions. The contribution of all variables based on the EWMA contribution approach can be calculated as follows CEWMA(r ) = (1 − γ )Ccom(r ) + γCEWMA(r − 1)

(53)

where γ is the forgetting factor and CEWMA(r − 1) is the mean of contribution of all the past r − 1 variables. The initial mean of contribution is defined as CEWMA(r − 1) =

1 L

5. APPLICATION TO CONTINUOUS STIRRED TANK REACTOR SYSTEM The effectiveness of the proposed RCVA-FOP monitoring method for time-varying processes is verified with a simulated continuous stirred tank reactor (CSTR) system.40−42 The diagram of the simulated CSTR system is shown in Figure 1. Under the assumption of a first order irreversible reaction, the reactant A and inlet flow of solvent produce a single component B, which is an outlet stream. The cooling flow of

r−1

∑ k=r−L+2

Ccom(k)

(54)

where L is the width of the initial moving window and r should be larger than the width of the initial moving windows. The combined contribution and EWMA contributions are plotted as a two-dimensional gray map, in which sampling time F


Article


Table 2. Measured Variables and Noise Standard Deviation Values

Figure 1. Diagram of the CSTR process.

Q f (Tf − T ) Ah

+

UAc (Tc − T ) ρCpAh

(55)

Fault Type

1

Qc sensor precision degradation Heat exchanger fouling Tcf sensor bias fault

Qf − Qo dh = dt A

(58)

Fault description

Fault time

f k ∼ N(0,5)

700

UAC = UAC − 120 J·min−2·K f k ≈ U(2,4)

701 699

varying characteristic were considered in each case study. The number of retained CVs was determined by using the CPV method, which can approximately explain 99% of the total variance. The charts of T2 and Q monitoring statistics are demonstrated in the subsequent subsections. 5.1. Sensor precision degradation fault detection and identification. A slow variation in Caf with a slope of 0.6 × 10−3 s−1 from k = 301 to k = 700 was considered in this case study. A sensor precision degradation fault was introduced into the coolant flow rate Qc at the sample time k = 700. The monitoring statistics of CVA and RCVA-FOP for the sensor precision degradation fault are shown in Figure 2. The conventional CVA approach gave false alarms during parameter variation. At that time, the RCVA-FOP approach can handle the natural change and remarkably decreased the false alarm rate. After this sensor fault occurred, the proposed approach could immediately detect fault without any delay. The contribution plot for the combined statistics is presented as shown in Figure 3(a). And the largest variable contribution charts are shown in Figure 3(b). Observed from that, the faulty variables are 2, 5, and 10. Variable 10 has the largest contribution, which means the coolant flow rate Qc is the root cause of this sensor precision degradation fault. Variable 5 also has less contribution, because the coolant temperature was affected by the variation of the coolant flow rate Qc in the cooling jacket. The correct fault identification rate is 71.21%, which is the percentage of variable 10 among all variables after sample 700. For comparison, the EWMA contribution plot for the combined statistics is presented in Figure 4(a). And the largest EWMA variable contribution is plotted in Figure 4(b). Though the correct faulty variable is 10, variables 2 and 5 possess a large contribution at the initial phase of the fault. This indicates that the smearing effect still exists in the EWMA contribution plot

(56)

(57)

0.0024 mol·L−1 0.0024 mol·L−1 0.45 K 0.45 K 0.45 K 0.45 K 0.004 m 0.71 L·min−1 0.71 L·min−1 0.32 L·min−1

Case

2 3

Q c(Tcf − Tc) dTc UAc (T − Tc) = + dt Vc ρc CpcVc

Standard deviation values

Ca Caf T Tf Tc Tcf h Qo Qf Qc

Table 3. Sensor Fault Simulation

k e−E / RT Ca( −ΔH ) dT = 0 dt ρCp +

Variable

1 2 3 4 5 6 7 8 9 10

K·min−1. The abnormal deviation from normal conditions was simulated by three types of faults, which are sensor precision degradation fault, heat exchanger fouling fault, and sensor bias fault. Three sets of fault types data were simulated. The detailed fault types are listed in Table 3. One of the faults and the time-

the jacket is utilized to transfer the heat generated from the reaction. The cascade control strategy is adopted to regulate the temperature and liquid level of the reactor. Based on the basic principle of the energy, mass, and component balances, the CSTR system’s dynamic model can be formulated as Q f (Caf − Ca) dCa = −k 0e−E / RT Ca + Ah dt

No.

where A is the cross-sectional area of the reactor in the CSTR system, Cp is the heat capacity of the contents, Ca is the concentration of species A in the reactor, Caf is the concentration of species A in the feed stream, Cpc is the heat capacity of coolant, h is the liquid level of the reactor, E is the activation energy, k0 is preexponentional factor, Qc is the coolant flow rate, Qf is the feed flow rate into the reactor, Qo is the outlet flow rate out of the reactor, R is the gas constant, T is the temperature of the reactor, Tcf is the temperature of the coolant feed, Tc is the coolant’s temperature in the cooling jacket, Tf is the temperature of the feed stream, U is the heattransfer coefficient, ΔH is the reaction heat, ρ is the contents density, Ac is the total heat transfer area, and ρc is the density of the coolant. Data sets were generated with the sample period of 10 s under normal operating mode in the CSTR simulation system. Table 2 lists ten measured variables and the corresponding Gaussian noise standard deviation values used in the simulation procedure. The CVA model was established by using the first 200 samples of the data set. The data were scaled to zero mean and unit variance. The time-varying characteristics were described by two types of natural changes from k = 301 to k = 700:22,40 (1) a slow variation in Caf with a slope of 0.6 × 10−3 s−1 and (2) a slow variation of catalyst deactivation E/R with 3 G


Article


Figure 3. Online fault identification for a precision degradation fault: (a) contribution plot; (b) variable with the largest contribution.

Figure 2. Monitoring charts for a precision degradation fault: (a) conventional CVA; (b) RCVA-FOP approach (dashed line: 99% control limit).

method. Nevertheless, the correct fault identification rate was improved to 95.48%, which is higher than that of the conventional contribution plot method. 5.2. Heat exchanger fouling fault detection and identification. A slow variation of catalyst deactivation E/R with 3 K·min−1 was introduced as a normal change of the process from k = 301 and continued to k = 700. And a process fault was introduced at k = 701 until the end of the simulation, which is a slow variation of heat exchanger fouling UAC with −120 J·min−2·K. The charts of monitoring statistics from the conventional CVA approach are shown in Figure 5(a). After the slow change of the process, the monitoring statistics largely exceeded the corresponding control limits because the model was not updated in a timely manner. That is why it caused a higher false alarm rate. Compared with the results of CVA, the monitoring statistics of the proposed approach were more reliable during the parameter variation process, as shown in Figure 5(b). After the fault occurred, the monitoring statistics detected it with a

Figure 4. Online fault identification for a precision degradation fault: (a) EWMA contribution plot; (b) variable with the largest EWMA contribution.

short time delay. Because of the minor impact of the heat exchanger fouling fault at the beginning, it is very difficult to identify the fault immediately. The contribution plot for the combined statistics is determined as shown in Figure 6(a). And the largest variable contribution is shown in Figure 6(b). Seen from the figure, variables 5 and 10 have a relatively large contribution. This means this fault is caused by one of the variables 5 or 10. Only the contribution of variable 5 increases gradually, which indicates the root cause of the fault is related to the temperature of the coolant in the cooling jacket Tc. The correct fault identification rate is 68.04%. For comparison, the EWMA contribution plot for the combined statistics is presented in Figure 7(a). And the largest EWMA variable H


Article


Figure 6. Online fault identification for a drift fault: (a) contribution plot; (b) variable with the largest contribution.

Figure 5. Monitoring charts for a drift fault: (a) the conventional CVA; (b) RCVA-FOP approach (dashed line: 99% control limit).

contribution is plotted in Figure 7(b). The correct fault identification rate was improved to 73.85%, which is a little higher than that of the conventional contribution plot method. 5.3. Bias sensor fault detection and identification. A slow variation in Caf with a slope of 0.6 × 10−3 s−1 was regarded as a normal change of the process. The slow change of the process began at k = 301 and continued to k = 698. A sensor bias fault of the temperature of the coolant feed Tcf was introduced at the sample time 699. The monitoring statistics generated from the conventional CVA and RCVA-FOP approaches are shown in Figure 8. It seems that the RCVA-FOP can effectively deal with the timevarying characteristic, and reduces the false alarm rate remarkably, as shown in Figure 8(b). In contrast, during the process of gradual change, the monitoring charts of statistics from the conventional CVA show a number of violations in Figure 8(a). The process variable contribution plot for the combined statistics is shown in Figure 9(a). And the largest variable contribution is shown in Figure 9(b). Variables 2, 5, and 6 have

Figure 7. Online fault identification for a drift fault: (a) EWMA contribution plot; (b) variable with the largest EWMA contribution.

a relatively large contribution, because the temperature of the coolant feed Tcf is the manipulated variable; its fault effect can propagate into other related variables, such as the coolant flow rate Qc and the coolant temperature in the cooling jacket Tc. In contrast, the EWMA contribution plot for the combined statistics is presented in Figure 10(a). And the largest EWMA variable contribution is plotted in Figure 10(b). Though the correct faulty variable is 6, variables 2 and 5 also possess a large contribution, caused by the smearing effect still of fault. Nevertheless, the correct fault identification rate was improved to 84.85%, which is higher than the 63.96% of the conventional contribution plot method. I


Article


Figure 9. Online fault identification for a bias fault: (a) contribution plot; (b) variable with the largest contribution.

Figure 8. Monitoring charts for a sensor bias fault: (a) the conventional CVA; (b) RCVA-FOP approach (dashed line: 99% control limit).

The false alarm rate (FAR) was adopted to measure the capability of detecting faults. The false alarm rate can be obtained as follows: FAR =

Nv × 100% Nn

Figure 10. Online fault identification for a bias fault: (a) EWMA contribution plot; (b) variable with the largest EWMA contribution.

Table 4. False Alarm Rates of the Conventional CVA and RCVA-FOP Approach (59) CVA

where Nv is the number of fault samples and Nn is the total number of normal samples. False alarm rates based on the conventional CVA and RCVA-FOP approaches are shown in Table 4, respectively. Both T2 and Q of the RCVA-FOP approach are smaller than that of the CVA. This indicates the proposed approach can significantly reduce false alarm rates compared with the results of CVA. For comparing the correct fault identification rate, a moving windows (MW) contribution approach is also adopted in simulation. The correct fault identification rate (FIR) is determined by calculating the number of faulty variables, with

2

case

T

1 2 3

49.5714 38.2857 47.5714

RCVA-FOP 2

Q

T

Q

13.7143 43 41.7534

0.1431 0.0127 0.1423

0.5714 0.5742 1.4286

the largest contribution compared to the number of total faulty data points according to the relation as follow: FIR = J

Fv × 100 Fn

(60) DOI: 10.1021/acs.iecr.6b02653 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


toward an oil sand primary separation. Comput. Chem. Eng. 2014, 71, 281−297. (2) Juricek, B. C.; Seborg, D. E.; Larimore, W. E. Fault detection using canonical variate analysis. Ind. Eng. Chem. Res. 2004, 43, 458− 474. (3) Mercère, G.; Bako, L.; Lecœuche, S. Propagator-based methods for recursive subspace model identification. Signal Processing 2008, 88, 468−491. (4) Houtzager, I.; van Wingerden, J.-W.; Verhaegen, M. Recursive predictor-based subspace identification with application to the realtime closed-loop tracking of flutter. IEEE T. Control Syst. T. 2012, 20, 934−949. (5) Qin, S. J. An overview of subspace identification. Comput. Chem. Eng. 2006, 30, 1502−1513. (6) Juricek, B. C.; Seborg, D. E.; Larimore, W. E. Process control applications of subspace and regression-based identification and monitoring methods 2005, 2341. (7) Larimore, W. E. Canonical variate analysis in control and signal processing. Statistical methods in control and signal processing 1997, 83− 120. (8) Negiz, A.; Cinar, A. Statistical monitoring of multivariable dynamic processes with state-space models. AIChE J. 1997, 43, 2002− 2020. (9) Russell, E. L.; Chiang, L. H.; Braatz, R. D. Fault detection in industrial processes using canonical variate analysis and dynamic principal component analysis. Chemom. Intell. Lab. Syst. 2000, 51, 81− 93. (10) Stubbs, S.; Zhang, J.; Morris, J. Fault detection in dynamic processes using a simplified monitoring-specific CVA state space modelling approach. Comput. Chem. Eng. 2012, 41, 77−87. (11) Lee, C.; Lee, I.-B. Adaptive monitoring statistics with state space model updating based on canonical variate analysis. Korean J. Chem. Eng. 2008, 25, 203−208. (12) Naik, A. S.; Yin, S.; Ding, S. X.; Zhang, P. Recursive identification algorithms to design fault detection systems. J. Process Control 2010, 20, 957−965. (13) Ge, Z.; Song, Z.; Gao, F. Review of recent research on databased process monitoring. Ind. Eng. Chem. Res. 2013, 52, 3543−3562. (14) Gustafsson, T. Instrumental variable subspace tracking using projection approximation. IEEE T. Signal Proces. 1998, 46, 669−681. (15) Oku, H.; Kimura, H. Recursive 4SID algorithms using gradient type subspace tracking. Automatica 2002, 38, 1035−1043. (16) Li, Y.; Xie, L.; Su, H.; Chu, J. Recursive closed-loop subspace identification based on predictor Markov parameters 2010, 3655. (17) Lee, C.; Wook Choi, S.; Lee, I.-B. Adaptive monitoring statistics based on state space updating using canonical variate analysis. Comput.-Aided Chem. Eng. 2006, 21, 1545−1550. (18) Ding, S.; Zhang, P.; Naik, A.; Ding, E.; Huang, B. Subspace method aided data-driven design of fault detection and isolation systems. J. Process Control 2009, 19, 1496−1510. (19) Pan, T.-H.; Sheng, B.-Q.; Wong, D. S.-H.; Jang, S.-S. A virtual metrology model based on recursive canonical variate analysis with applications to sputtering process. J. Process Control 2011, 21, 830− 839. (20) Choi, S. W.; Martin, E. B.; Morris, A. J.; Lee, I.-B. Adaptive multivariate statistical process control for monitoring time-varying processes. Ind. Eng. Chem. Res. 2006, 45, 3108−3118. (21) Elshenawy, L. M.; Yin, S.; Naik, A. S.; Ding, S. X. Efficient recursive principal component analysis algorithms for process monitoring. Ind. Eng. Chem. Res. 2010, 49, 252−259. (22) Elshenawy, L. M.; Awad, H. A. Recursive fault detection and isolation approaches of time-varying processes. Ind. Eng. Chem. Res. 2012, 51, 9812−9824. (23) Jiang, B.; Huang, D.; Zhu, X.; Yang, F.; Braatz, R. D. Canonical variate analysis-based contributions for fault identification. J. Process Control 2015, 26, 17−25. (24) Yue, H. H.; Qin, S. J. Reconstruction-based fault identification using a combined index. Ind. Eng. Chem. Res. 2001, 40, 4403−4414.

where Fv is the number of faulty variables with the largest contribution and Fn is the total number of faulty samples. The detailed fault identification rates for different cases are shown in Table 5. Fault identification rates based on the EWMA Table 5. Comparison of Fault Identification Rates Case

Combined contribution

MW contribution

EWMA contribution

1 2 3

71.21% 68.04% 63.96%

93.43% 71.65% 74.11%

95.48% 73.85% 84.85%

contribution in cases 1 and 3 are much higher than that of the combined contribution approach. Because of the effect of the drift fault, fault identification rates based on the EWMA contribution are slightly higher than that of the combined contribution approach. In summary, the proposed recursive monitoring approach can not only handle the time-varying characteristics but also detect and identify three types of faults. The root cause of faults can be effectively identified by the proposed EWMA contribution approach. To a certain extent, the correct fault identification is improved.

6. CONCLUSIONS For monitoring of time-varying processes, a recursive canonical variate analysis based on the first order perturbation theory approach is proposed. Without recalling normal operating data after operating condition change, covariance and crosscovariance matrices were updated by using the exponential weighted moving average. To identify the root cause of faults, a EWMA contribution chart is also proposed to do online fault identification. The proposed approach is verified with the simulation in the continuous stirred tank reactor. Compared with the simulation results of CVA, the results of the proposed RCVA-FOP not only can adapt the time-varying characteristics of industrial processes, but also can identify three types of faults, which are sensor precision degradation fault, heat exchanger fouling fault, and sensor bias fault. The correct fault identification rate based on the EWMA contribution approach is improved to a certain extent. An online fault identification approach which can eliminate the effect of fault smearing should be further studied in the future.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (NSFC) under Grant 61374137, 61325015 and the State Key Laboratory of Integrated Automation of Process Industry Technology and Research Center of National Metallurgical Automation Fundamental Research Funds (2013ZCX02-03).

■

REFERENCES

(1) Sammaknejad, N.; Huang, B.; Fatehi, A.; Miao, Y.; Xu, F.; Espejo, A. Adaptive monitoring of the process operation based on symbolic episode representation and hidden Markov models with application K


Article

Industrial & Engineering Chemistry Research (25) Yoon, S.; MacGregor, J. F. Fault diagnosis with multivariate statistical models part I: using steady state fault signatures. J. Process Control 2001, 11, 387−400. (26) Singhal, A.; Seborg, D. E. Pattern matching in multivariate time series databases using a moving-window approach. Ind. Eng. Chem. Res. 2002, 41, 3822−3838. (27) Singhal, A.; Seborg, D. E. Evaluation of a pattern matching method for the Tennessee Eastman challenge process. J. Process Control 2006, 16, 601−613. (28) Gertler, J.; Cao, J. PCA-based fault diagnosis in the presence of control and dynamics. AIChE J. 2004, 50, 388−402. (29) Shang, L.; Liu, J.; Turksoy, K.; Shao, Q. M.; Cinar, A. Stable recursive canonical variate state space modeling for time-varying processes. Control Eng. Pract. 2015, 36, 113−119. (30) Stewart, G. W.; Sun, J.-g. Matrix Perturbation Theory; 1990. (31) Willink, T. J. Efficient adaptive SVD algorithm for MIMO applications. IEEE T. Signal Proces. 2008, 56, 615−622. (32) Chen, Z.; Zhang, K.; Hao, H.; Ding, S. X.; Krueger, M.; He, Z. A Canonical Variate Analysis based Process Monitoring Scheme and Benchmark Study 2014, 47, 10634. (33) Odiowei, P.-E.; Cao, Y. Nonlinear dynamic process monitoring using canonical variate analysis and kernel density estimations. IEEE T. Ind. Inform. 2010, 6, 36−45. (34) Yang, Y.; Chen, Y.; Chen, X.; Liu, X. Multivariate industrial process monitoring based on the integration method of canonical variate analysis and independent component analysis. Chemom. Intell. Lab. Syst. 2012, 116, 94−101. (35) Lee, C.; Choi, S. W.; Lee, I.-B. Variable reconstruction and sensor fault identification using canonical variate analysis. J. Process Control 2006, 16, 747−761. (36) Lu, J.; Liu, C.-L. Statistical modeling of dynamic multivariate process using canonical variate analysis 2006, 218. (37) Pilgram, B.; Judd, K.; Mees, A. Modelling the dynamics of nonlinear time series using canonical variate analysis. Phys. D 2002, 170, 103−117. (38) Simoglou, A.; Martin, E.; Morris, A. Statistical performance monitoring of dynamic multivariate processes using state space modelling. Comput. Chem. Eng. 2002, 26, 909−920. (39) Zhu, X.; Braatz, R. Two-Dimensional Contribution Map for Fault Identification [Focus on Education]. Control Systems, IEEE 2014, 34, 72−77. (40) Johannesmeyer, M. C.; Singhal, A.; Seborg, D. E. Pattern matching in historical data. AIChE J. 2002, 48, 2022−2038. (41) Singhal, A.; Seborg, D. E. Effect of data compression on pattern matching in historical data. Ind. Eng. Chem. Res. 2005, 44, 3203−3212. (42) Deng, X.; Tian, X. Multimode Process Fault Detection Using Local Neighborhood Similarity Analysis. Chin. J. Chem. Eng. 2014, 22, 1260−1267.

L


Recursive Fault Detection and Identification for Time-Varying

Recommend Documents