Efficient Recursive Principal Component Analysis Algorithms for

Nov 16, 2009 - These false alarms threaten the credibility of the monitoring system. Therefore, recursive PCA (RPCA) algorithms are recommended. The m...
1 downloads 0 Views 243KB Size
Download PDF
252

Ind. Eng. Chem. Res. 2010, 49, 252–259

Efficient Recursive Principal Component Analysis Algorithms for Process Monitoring Lamiaa M. Elshenawy,*,† Shen Yin,‡ Amol S. Naik,‡ and Steven X. Ding‡ Department of Industrial Electronics and Control Engineering, Faculty of Electronic Engineering, UniVersity of Minufiya, 32952 Menouf, Egypt, and Institute of Automatic Control and Complex Systems, UniVersity of Duisburg-Essen, 47057 Duisburg, Germany

Principal component analysis (PCA) has been successfully applied in large scale process monitoring. However, classical PCA has some drawbacks: one of these aspects is the inability to deal with parameter-varying processes, where it interprets the natural changes in the process as faults, resulting in numerous false alarms. These false alarms threaten the credibility of the monitoring system. Therefore, recursive PCA (RPCA) algorithms are recommended. The most important challenge faced by these algorithms is the high computation costs, due to repeated eigenvalue decomposition (EVD) or singular value decomposition (SVD). Motivated by this issue, we present two RPCA algorithms that will greatly reduce the computation cost. The first algorithm is based on first-order perturbation analysis (FOP), which is a rank-one update of the eigenvalues and their corresponding eigenvectors of a sample covariance matrix. The second one is based on the data projection method (DPM), which is a simple and reliable approach for adaptive subspace tracking. The effectiveness of the presented RPCA algorithms is evaluated with an application of monitoring a nonisothermal continuous stirred tank reactor (CSTR) system. The results show the efficiency of these approaches compared to the classical PCA. 1. Introduction Modern industrial processes are large scale interconnected systems. Thus, efficiency of any data-driven monitoring scheme depends upon its ability to compress a huge amount of process data and extract the meaningful information within. One of the most common multivariate statistical process control (MSPC) methods used for this purpose is principal component analysis (PCA). The basic strategy of PCA is to discard the noise, while preserving the most important information of the original data set. Although PCA was introduced by Pearson1 and developed by Hotelling2 in the early 20th century, the vast majority of applications were developed only in the last few decades with an increasing use of computers.3 Many successful applications of PCA have been reported in the literature. These applications involve data compression, feature extraction, image processing, pattern recognition, signal processing, factor analysis, time series analysis, chemometrices,4-6 and fault detection and diagnosis.7,8 PCA-based process monitoring, especially for chemical processes, has received much attention in the last two decades. Despite its tremendous success, classical PCA has a few major shortcomings, e.g., monitoring of parametervarying processes. Therefore, it is necessary to develop a method that adapts the underlying PCA model in a continuous and automatic manner in compliance with present process conditions. In fact, there are a few approaches that have already been developed. Wold9 proposed a so-called exponentially weighted moVing aVerage (EWMA) approach for updating of both PCA and partial least-squares (PLS) models. Li et al.10 proposed an RPCA-based adaptive process * To whom correspondence should be addressed: E-mail: [email protected]. † University of Minufiya. ‡ University of Duisburg-Essen.

monitoring scheme which includes a (i) recursive update of the data correlation matrix, a (ii) recursive calculation of the PCA model structure instead of performing SVD or EVD at every time step, and based on it, a (iii) recursive determination of the confidence limits for process monitoring purposes. Recursive approaches have also been considered in adaptive statistical process control11-14 and parity spacebased fault detection.15,16 Another promising technique combining moving window (MW) and recursive updating of a PCA model was proposed by Wang et al.17 The fast MWPCA enables online application of a generic moving window-based recursive PCA with a larger window size. To overcome the difficulty of incipient fault detection, an N-Step-Ahead predictor is also proposed. Recently, Liu et al.18 have also presented the so-called moVing window kernel PCA for adaptive monitoring of nonlinear processes. The adaptation mechanism is based on the incremental method proposed by Hall et al.19-21 The data covariance matrix is normalized after projecting it to the feature space and an approximate eigendecomposition can be computed with the help of up- and downdating procedure. In this paper, two new RPCA algorithms are proposed, motivated by the development in the field of signal processing, especially for subspace tracking applications. The major advantages of these algorithms are low computation cost and simplicity of online implementation. The first algorithm is based on first-order perturbation theory (FOP), which is a rank-one update of the eigenpairs of the data covariance matrix.22,23 The second one is based on data projection method (DPM), which serves as a simple and reliable approach for adaptive subspaces tracking.24 The paper is organized as follows: In Section 2, the problem formulation and preliminaries about the classical PCA are provided. The two algorithms to recursively update the PCA model are illustrated in Section 3. Section 4 gives

10.1021/ie900720w CCC: $40.75  2010 American Chemical Society Published on Web 11/16/2009

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

the integrated RPCA-based process monitoring scheme and a comparison of the computation costs of different algorithms. In Section 5, the performance of the proposed monitoring scheme is investigated with the help of a nonisothermal continuous stirred tank reactor (CSTR) benchmark, and a conclusion and summary are given in Section 6. 2. Problem Formulation and Preliminaries The classic PCA-based process monitoring involves three main steps: (i) collecting process data, (ii) compressing it with the help of decomposition algorithm, and (iii) choosing testing statistics to monitor process behavior. The data compression is usually done with eigenvalue decomposition of the covariance matrix of the sampled measurements. Let X ∈ Rn×m denotes the data matrix which records n samples of m variables. The matrix X is then scaled to zero mean for covariance-based PCA and further to unit variance for correlation-based PCA. The PCA model can be expressed as: X ) TPT ) TˆPˆT + T˜P˜T ) Xˆ + X˜

(1)

where Xˆ and X˜ stand for predicted and residual parts of X. The matrices T and P are called score and loading matrices, respectively. The loading vectors, i.e., columns of matrix P, are actually eigenvectors of the covariance matrix of X. The covariance matrix and its eigenvalue decomposition are expressed as S)

[

]

Λlk O 1 XTX ) [VlkVm-lk] [VlkVm-lk]T Λ n-1 O m-lk

(2)

The number of most significant eigenvalues, i.e., lk determines the acceptable degree of compressibility without losing much of information. It can be determined in various ways, and the details can be found elsewhere in the literature. Alternatively, matrix X can also be directly decomposed using singular value decomposition (SVD). To monitor the process status, a measurement sample is projected with the help of the PCA model to either the predicted or the residual subspace. Two monitoring statistics are developed based on these projections, namely T2 and Q-statistic (squared prediction error (SPE)). To complete the design, we need only to choose a suitable threshold to enable detecting the smallest faults with a minimal rate of false alarms. A process behavior can be termed “safe” if the process statistic remains below these thresholds. In reality, it is always recommended to use both T2 and SPE index, simultaneously.25 During the process operation, a temporary change such as shift in operating point, short-term disturbance, or variation in parameters could occur. These changes are not necessarily critical, yet a process monitoring system could indicate it as a malicious behavior. A classic PCA-based monitoring must then involve “switching” a new model to monitor the process during and after these changes. Instead, a possibility is to use an adaptive model that does not require hard switching or tuning to keep monitoring the process during these changes and to avoid false alarms. The recursive PCA-based process monitoring attempts to deal with such problems. To begin with, it is necessary to compute the data covariance matrix Sk and to update its eigenstructure recursively. The computation cost associated to perform it is very high and in many cases prohibitive. Therefore, based on eq 2, it is desired to have an approximate eigendecomposition of Sk:

Sk )

1 XTX k-1 k k

253

(3)

and ) [Vk-1 + ∆Vk][Λk-1 + ∆Λk][Vk-1 + ∆Vk]T

(4)

where (∆Vk, ∆Λk) can be calculated from the old eigenpairs (Vk-1, Λk-1) and the new measurement data xk. In the next section, we introduce two algorithms to effectively compute the eigendecomposition of covariance matrix, which in turn will solve the problem of designing recursive process monitoring based on PCA. 3. Recursive PCA Algorithm 3.1. RPCA Based on FOP (Algorithm 1). FOP-based recursive computation of eigendecomposition implicitly updates the covariance matrix. In simple words, the covariance matrix is interpreted as a product of eigenvalues and eigenvectors, and as new sampled data becomes available, the eigenstructure is updated as a whole, instead of the covariance matrix itself.22 Recently, Willink23 proposed a recursive approach for SVD based on FOP. In fact, the roots of both these approaches can be traced back to the following lemma.26 Lemma 1: Given the following recursive update equation of the covariance matrix Sk ) Sk-1 + ε[xkxTk - Sk-1] where Sk ∈ Rm×m, ε is a small positive number (εf0), xk ∈ Rm is the sample data vector at the kth time instant, and then the eigenvalues and eigenvectors of the updated matrix Sk can be described by: λk,i ) (1 - ε)λk-1,i + ε|fi | 2

(5)

m

Vk,i ) Vk-1,i +

∑b V

(6)

ji k-1,j

j)1

with

bji )

{

T xk fi ) Vk-1,i

(7)

0

j)i

fjfi -bij ) (λk-1,i - λk-1,j)

j*i

(8)

where i,j ) 1, · · · , m. A stabilization mechanism is necessary to prevent an undesirable behavior of the algorithm, which may occur when there exists a very small difference between two successive eigenvalues. This problem can be solved, without significantly altering the complexity of the algorithm by using the following modification: bji )

f j fi , j*i max(γλk-1,1, λk-1,i - λk-1,j)

(9)

where γ is a small positive number. Note that it can be further assumed that λk-1,i . λk-1,j for i < j to obtain a low complexity algorithm. The resulting algorithm is still robust even if this assumption is not holding.23 As a result, eq 8 is approximated for i < j as

254

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

bji )

fjfi λk-1,i

Table 2. Recursive Eigendecomposition Algorithm Based on DPM

(10)

The complete FOP-based eigenValue decomposition (EVD) algorithm is summarized in Table 1. FOP-based EVD updates entire eigenpairs with a computation complexity of (O(m2)). However, since PCA requires only a few of these eigenpairs for the process monitoring task, a significant further reduction in computation cost is also available. 3.2. RPCA Based on DPM (Algorithm 2). DPM is a fast and numerically robust approach that depends on the orthogonal iteration procedure (OIP).27 It has been recently proposed for the application in subspace tracking.24 The main advantages of this algorithm are (i) simple implementation, (ii) capability to update both signal and noise subspace with a simple formulation, and (iii) the low computation complexity of O(ma) with a < m denoting the dimension of the tracked subspace. First, we introduce the following lemma, which gives the eigenvectors corresponding to the l dominant eigenvalues of a symmetric, nonnegative definite matrix. Lemma 2: Let H ∈ Rm×m be a symmetric, nonnegative definite matrix, with eigenvalues λ1 g · · · g λl > λl+1 g · · · g λm g 0 and the corresponding eigenvectors V1, · · · , Vl, Vl+1 · · · , Vm. Consider the sequence of matrices Vk ∈ Rm×l, defined by the iteration: Vk ) orthnorm{HVk-1}

(11)

where orthnorm stands for the orthonormalization procedure using QR-decomposition or Gram-Schmidt orthogonalization method, then limkf∞Vk ) [V1, · · · , Vl]

Vk ) orthnorm{(Im + µH)Vk-1}

(13)

Vk ) orthnorm{(Im - µH)Vk-1}

(14)

where µ is a small positive scalar parameter and Im ∈ Rm×m is an identity matrix. The underlying idea of these two variants is based on the identical eigenvectors of H and I ( µH, while the eigenvalues of I ( µH are equal to 1 ( µλ. Remarks: (i) For the recursive update of the eigenpairs of a covariance matrix, H is replaced by Hk, which can be calculated as: Hk ) RHk-1 + (1 -

(15)

Table 1. Recursive Eigendecomposition Algorithm Based on FOP 1. Initial conditions: Vk-1,i and λk-1,i, for i ) 1, · · · , m, are calculated from the last sample time k - 1 2. Collect new data vector xk 3. Compute xk ) (ε)1/2xk T 4. Compute fi ) Vk-1,i xk, 5. With initial condition pV0 ) xk and qV1 ) 0, calculate V piV ) pi-1 - fiVk-1,i V qi+1 ) qiV + (fi)/(λk-1,i)Vk-1,i 6. Update the ith eigenvalue λk,i ) (1 - ε) λk.i-1 + fi2 7. Update the ith eigenvector Vk,i ) Vk-1,i + (fi)/(λk-1,i)piV - fiqiV Vk,i ) (Vk,i)/(|Vk,i|)

where 0 e R < 1 is the forgetting factor. Because the mean projector trajectory E[VkVTk ] and the steady-state projection error power limkf∞ E[|VkVTk - VVT|F2], where | · |F denotes the Frobenius norm, are independent of the value R in Hk,24 the simplest selection for Hk is the instantaneous estimate of the covariance matrix, i.e., Hk ) xkxTk . (ii) In order to guarantee the non-negativeness of matrix Im ( µHk, the positive scalar µ should be appropriately chosen such that 0 < µ < 1 ak ) rk - |rk| e, e ) [1 0 . . . 0]T ∈ Rl Zk ) Tk - (2/|ak|2){Tkak}akT Vk ) normalization of the columns of Zk 8. If subspace rank l ) 1 Vk ) (Tk)/(|Tk|)

1. 2. 3. 4. 5. 6. 7.

4. Process Monitoring Based on RPCA 4.1. Data Preprocessing. To effectively extract the information in the data relevant to process monitoring, it is necessary to pretreat the data. One of the pretreatment procedures is the so-called autoscaling, which normalizes the process variables in a way that each variable has equal weight before applying the process monitoring method. Autoscaling consists of two steps, i.e., centering the variables to zero mean and then dividing these variables by their standard deviations.3 For recursive implementation of PCA, the mean and variance should be updated according to the following equations: mk ) Rmk-1 + (1 - R)x0k

(17)

2 σ2k,i ) Rσk-1,i + (1 - R)(x0k,i - mk,i)2

(18)

where x0k is the measurement data vector before normalization, 2 are the update mean vector mk ) [mk,1, · · · , mk, m]T ∈ Rm and σk,i and variances of measurement data, and i ) 1, · · · , m. According to these updated relations, the normalized data is calculated as follows: xk ) (x0k - mk)∑-1 k

(19)

where Σ ) diag(σk,1, · · · , σk,m). 4.2. Determination of Number of PCs. A complete implementation of RPCA also requires recursive calculation of the number of PCs. There are many approaches to calculate the number of PCs, but most of them use monotonically increasing or decreasing indices. The decision to choose the index is very

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

255

subjective and depends on application requirements. Several comparative studies have been conducted on these approaches.28 Unfortunately, not all of the approaches are suitable for the recursive determination of the number of PCs. The CPV approach is implemented in our application study.10,28,29 4.3. Process Monitoring Based on Algorithm 1. Algorithm 1 updates eigenvectors and eigenvalues directly with the equations given in Table 1. The number of PCs, lk, can be determined recursively by using CPV index. Another important issue for RPCA-based process monitoring is to update the process monitoring indices and their control limits. The T2 index can be defined according to the adaptive algorithm:

where pm is the eigenvector corresponding to the smallest eigenvalue. Under the assumption that the process is corrupted by white Gaussian noise, the residual signal will follow χ2 distribution and the upper limit of the scaled residual signal can be calculated according to the following relation:

ˆT T2k ) xTk PˆlkΛ-1 lk Plkxk

4.5. Summary of Steps. We present a summary of main steps involved in integrated RPCA-based process monitoring. Algorithm for RPCA-based Process Monitoring: • The initial conditions are calculated for sample time k - 1 2 2 , λk-1,i, Vk-1,i, lk-1, Tlim and For Algorithm 1: mk-1, σk-1,i k-1 Qlimk-1 For Algorithm 2: Hk-1, Vk-1, and Jth • Collect new data vector xk and normalize it using eq 19. • Calculating monitoring indices. For Algorithm 1: T2k , Qk according to eqs 20 and 22. For Algorithm 2: rk using eq 26. • If the monitoring indices exceed their control limits for the last samples, then there is a fault. The updating algorithm is terminated, and the alarm is released. Otherwise, the update procedure can be continued with the following steps: 2 by eqs 17 and 18. Update mk, σk,i For Algorithm 1: Update λk,i, Vk,i, T2limk-1 and Qlimk-1 using Table 1, (21) and (23). Calculate lk by e.g. CPV or AE For Algorithm 2: Update Hk, Vk, using eq 15, Table 2. • Increase sample time k by one and go to the second step. 4.6. Comparison of Adaptive PCA Algorithms. In this subsection, we compare the computation complexity of the two proposed RPCA algorithms with other known adaptive algorithms, e.g., standard SVD,27 Lanczos approach,10 inverse iteration approach,27 fast MWPCA (MWPCAfa)17 and MWPCA using incremental method (MWPCAIM).18 In the case of samplewise updating, the overall computation complexity of standard SVD27 approach is 22m3 + 4m2 + 12m + 1, whereas, Lanczos tridiagonalization requires about (4lk - 1)m2 + (6.5lk - 1.5l2k )m - 2l2k flops for Hk10,18 and the overall computation cost is equal to (4lk + 3)m2 + (6.5lk - 1.5l2k + 12)m - 2l2k + 1. The inverse iteration approach needs (2)/(3)m3lk + 2lkm2 + (5lk + 12)m + lk + 1 flops to update lk eigenpairs. Moreover, fast MWPCA17 contains 6m3 + 20m2 + 11m + 17 flops and the computation cost of MWPCA using incremental method18 is equal to (2lk + 8)m2 + 18l3k + 33l2k - 4m + 33lk + 9. The computation complexity of the first RPCA algorithm, i.e., Algorithm 1, is 14m2 + 19m + 1 and that of Algorithm 2 is

(20)

where Pˆlk and Λlk are recursively calculated according to the updated PCA model structure. The upper limit for the T2 with a confidence level δ can be updated using the relation: 2 ) χδ2 (lk) Tlim k

(21)

The Q statistic is given as Qk ) xTk (I - PˆlkPˆTlk)xk

(22)

The control limit for the Q statistic will be updated as:10,30

Qlimk

[

cR√2θ2h20 θ2h0(h0 - 1) ) θ1 + +1 θ1 θ2 1

]

1/h0

(23)

where h0 ) 1 -

2θ1θ3 3θ22

(24)

m

∑

θj )

λji,

j ) 1, 2, 3

∑λ,

j ) 1, 2, 3

i)lk+1

or lk

θj )

trace(Sjk)

-

j i

(25)

i)1

4.4. Process Monitoring Based on Algorithm 2. As mentioned in the previous section, DPM approach can track both the signal and noise subspaces with a simple change of sign. Accordingly, there are two ways to apply Algorithm 2 for the purpose of process monitoring. The first one is to calculate the eigenvectors corresponding to the dominant eigenvalues and using the relations of PCA-based monitoring indices to calculate both the Hotelling’s T2 and Q statistic, given by eqs 20-25. The second one uses the eigenvector corresponding to the smallest eigenvalue to generate a residual signal. Here, we take the second way, which means to update only the eigenvector corresponding to the smallest eigenvalue; moreover, the computation complexity will decrease from O(mlk) to O(m). After this eigenvector is updated according to Table 2, a residual signal is generated as: rk )

xTk pmpmTxk

(26)

Jth ) χδ2 (F)

(27)

where F is the degrees of freedom (in our approach, F ) 1), respectively. Finally, we have the simple decision logic:

{

rk e Jth no fault rk > Jth there is a fault

Table 3. Comparison of RPCA Algorithms approach

flops

complexity

standard SVD Lanczos approach

22m + 4m + 12m + 1 (4lk + 3)m2 + (6.5lk - 1.5lk2 + 12)m 2l2k + 1 (2/3)m3lk + (2lk + 4)m2 + (5lk + 12)m + lk + 1 6m3 + 20m2 + 11m + 17 (2lk + 8)m2 + 18l3k + 33l2k - 4m + 33lk + 9 14m2 + 19m + 1 (13lk + 14)m + 7lk + 8 or 22m + 7

O(m3) O(m2lk)

inverse iteration MWPCAfa MWPCAIM FOP DPM

3

2

O(m3) O(m3) O(m2lk) O(m2) O(mlk) or O(m)

256

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010 Table 6. Simulation Cases

Figure 1. Nonisothermal CSTR process and nine measured variables: (1) coolant temperature, (2) reactant mixture temperature, (3) reactant A concentration, (4) solvent concentration, (5) solvent flow rate, (6) solute flow rate, (7) coolant flow rate, (8) outlet concentration, and (9) outlet temperature. Table 4. Model and PI Controller Parameters notation

parameters and constants

value

V F Fc Cp Cpc ∆Hr k0 E/R Kc(T) τ1(T) Kc(C) τ1(C)

volume of reaction mixture in the tank density of reaction mixture density of coolant specific heat capacity of the mixture specific heat capacity of the coolant heat of reaction exponential kinetic constant activation energy/ideal gas constant temperature controller gain integral time concentration controller gain integral time

1.0 m3 106 g/m3 106 g/m3 1 cal/(g K) 1 cal/(g K) -1.3 × 107 cal/kmol 1010 min-1 8330K -1.5 5.0 0.4825 2.0

(13lk + 14)m + 7lk + 8 for lk > 1. If the eigenvector corresponding to the smallest eigenvalue is used for the process monitoring purpose as mentioned above, the computation cost of Algorithm 2 is further reduced to 22m + 7. We summarize the findings in Table 3 given below. Here m and lk are the number of system variables and the dimension of the signal subspace, respectively. Table 3 shows that the two recursive algorithms proposed in this paper considerably reduce the computation cost. Note that both the MWPCA using incremental method and Algorithm 2 are not necessary to update the complete eigendecomposition, i.e., only the lk , m most significant eigenpairs instead of the entire eigendecomposition can be adapted. Moreover, Algorithm 2 provides the possibility to update the nonsignificant eigenpairs, which can be directly used for process monitoring purposes. This further reduces the computation complexity. Despite the higher computation cost compared to proposed algorithms, the MWPCA using incremental method can also be applied by including a sufficient number of data in the time-window that avoids the difficulty in selecting a forgetting factor. On the other hand, for RPCA algorithms such as Algorithm 1 and Algorithm 2, the covariance matrix is based on the ever-growing data set, which would render the adaptive scheme less sensitive in detecting incipient faults.17 To overcome this difficulty, the forgetting factor can be applied to provide less weight on

cases

description

magnitude

slow parameter-varying

βr(t) ) βr(t) - r1(t - 900) + r2(t - 1000)

set-point change slow parameter-varying with sensor fault

C(t) ) C(t) + r3 βr(t) ) βr(t) - r4(t - 900) + r5(t - 1000) Ti(t) ) Ti(t) + f1(t - 1500) f1(t - 1800)

r1 ) 10-3 r2 ) 10-3 r3 ) 0.05 r4 ) 10-3 r5 ) 10-3 f1 ≈ U(2, 3)

older samples. The forgetting factor can be chosen based on a priori knowledge of the process. Algorithm 1 and Algorithm 2 both reduce the computational cost compared to the first five algorithms listed in Table 3. Algorithm 2 has potential to further reduce the cost only when the number of significant eigenvalues is extremely small compared to the total number of eigenvalues, i.e., lk , m. In the special case of lk ) m, both Algorithm 1 and Algorithm 2 have the same amount of computational complexity. It should also be noted that since PCA-based monitoring usually employs both T2 and Q -statistic, Algorithm 2 may also be performed to update complete eigendecomposition. Even then, the total cost would still be of the order O(m2). 5. Benchmark Study We now apply the proposed RPCA-based monitoring schemes to the nonisothermal continuous stirred tank reactor (CSTR).13,31 The schematic diagram of CSTR and its variables are shown in Figure 1. The reaction is of the first order (A f B) where reactant A premixed with a solvent, converts into product B with a rate r ) βrk0e-E/(RT)C

(28)

The dynamic behavior of this system is described as follows: V

VFCp

dC ) F(C - Ci) - Vr dt

(29)

dT UA ) FCpF(Ti - T) (T - Tc) + dt 1 + UA/2FcFcCpc (-∆Hr)Vr

(30)

where UA is the heat transfer coefficient which is related to the coolant flow by the empirical relation UA ) βUAaFbc . The inlet reactant concentration Ci obtained from the two feed streams, the reactant and the solvent, is calculated as Ci )

(FaCa + FsCs) (Fa + Fs)

(31)

Table 5. Process Inputs, Outputs, Disturbances, and Measurement Noise variable

mean

βr βUA Tc Ti Fs Ca Cs Fa Fc C T

1.0 1.0 365 K 370 K 0.9 m3/min 19.1 kmol/m3 0.3 kmol/m3 0.1 m3/min 15 m3/min

set-point

φ

d

k

σ2V

0.25 0.2

1.9 × 10-3 9.75 × 10-4 4.75 × 10-2 4.75 × 10-2 1.9 × 10-2 4.75 × 10-2 1.875 × 10-3

0.9 0.95 0.2 0.2 0.1 2 0.2 0.8 kmol/m3 368.25 K

σe2

2.5 × 10-3 2.5 × 10-3 4.0 × 10-6 1.0 × 10-2 2.5 × 10-5 4.0 × 10-6 1.0 × 10-2 2.5 × 10-5 4.0 × 10-4

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

257

The outlet temperature (T) and concentration (C) are regulated with a proportional-integral (PI) controller by manipulating both the inlet coolant flow (Fc) and the inlet reactant flow (Fa), respectively. The details of the controller and model parameters are given in Table 4. It is assumed that the stochastic disturbances stem from poisoning of the reaction and from fouling of the cooling coils as described in Choi et al.13 These stochastic disturbances are introduced by the parameters βr and βUA, respectively. All process inputs and disturbances are generated using a first-order autoregressive models or sinusoidal signals according to the following relations: x(t) ) φx(t - 1) + V(t)

(32)

x(t) ) d sin(kt) + V(t)

(33)

or

Figure 4. Set-point change: (a) Time series plot of model parameters and mean values for the outlet concentration C. (b) Time series plot for F c, F a.

where the process noise V(t) ∼ (0, σ2V). In addition, all measured variables are contaminated with white Gaussian noise e(t) ∼ (0, σ2e ); see Table 5. We have selected the sampling time equal to 1 s for the entire simulation run. The performance of the process monitoring algorithms is tested in three simulation cases: (i) slow parameter variation, (ii) an abrupt change in one of the process set-points, (iii) parameter variation during a fault, see Table 6. The tests are completed in following steps: (i) the process variables are measured at each sampling time, (ii) 300 samples are used to build the initial PCA model, (iii) the process is simulated again and 2000 measurements are collected. The number of significant PCs is calculated by using CPV, such that the variance explained is approximately 99% of the total variance.

Figure 5. Set-point change: (a) Monitoring indices of PCA model. (b) Monitoring indices of Algorithm 1 and (c) residual term of Algorithm 2. Figure 2. Slow parameter variation: plot of model parameters and mean values of Fc, Fa.

Figure 3. Slow parameter variation: (a) monitoring indices of PCA model, (b) monitoring indices of Algorithm 1, and (c) residual term of Algorithm 2.

The first simulation case is a slow parameter variation. The parameter βr changes with a slope of (-10-3) s-1 from time, k ) 901 to k ) 1000. As a result, the reaction rate of the chemical process decreases according to eq 28. Both the reactant flow rate Fa and coolant flow rate Fc decrease as well due to the closed loop and maintain the desired outlet concentration and the temperature as seen in Figure 2. The number of PCs for Algorithm 1 is three. The process monitoring indices of the classical PCA exceed their thresholds once the parameter variation begins, i.e., the process monitoring system indicates a fault, see Figure 3a, whereas RPCA-based monitoring shows better performance in dealing with the change as shown in Figures 3b,c. In the second simulation case, an abrupt change in set-point takes place. The reference for the outlet concentration C changes from 0.8 to 0.85 kmol/m3 at time k ) 901. As a result of this change, Fa and Fc increase; see Figures 4a,b. The classical PCA releases false alarms just after the set-point change as illustrated in Figure 5a. The results obtained by the two recursive algorithm are plotted in Figures 5b,c. It is evident that the recursive algorithms more effectively deal with abrupt changes than the classical PCA. In the third case, a slow parameter variation is simulated with a fault in one of the sensors. The parameter variation is the

258

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010

monitoring results. The reconstruction of a measurement vector excluding these outliers is an often recommended solution suggested in the literature.10,13 In future work, the RPCA algorithms proposed here could be integrated with the reconstruction techniques. A further modification can be applied by intelligently terminating the adaptive computation once a fault is successfully detected. Acknowledgment

Figure 6. Slow parameter variation with sensor fault: (a) Time series plot for Fc, Fa. (b) Time series plot for Ti.

The authors thank the reviewers for their valuable comments and constructive suggestions which have significantly improved the manuscript. The first author thanks the Egyptian government for the financial support during her stay. Literature Cited

Figure 7. Slow parameter variation during sensor fault: (a) Monitoring indices of PCA model. (b) Monitoring indices of Algorithm 1 and (c) residual term of Algorithm 2.

same as that in the first case. The sensor measuring the inlet temperature Ti produces biased readings from time k ) 1501 to k ) 1800. The bias values follow uniform distribution within the interval, [2, 3]. Figure 6a,b shows the changes in the process parameters Fc, Fa, and Ti, respectively. From Figure 7a, it can be seen that the monitoring indices of classical PCA exceed their control limits during the fault. But once the process parameters reach steady states, the two recursive algorithms detect the sensor fault correctly as shown in Figures 7b,c. 6. Conclusion Because most of the real industrial processes exhibit slow parameter changes which are a part of normal process operation, it is necessary to develop an adaptive monitoring scheme to deal with it. In this paper, we have proposed two such algorithms based on first-order perturbation theory and data projections integrated into a PCA-based process monitoring scheme. These algorithms significantly reduce the online computation cost compared to other known approaches. The complete monitoring system is applied to the CSTR benchmark, and results are demonstrated with frequently occurring process scenarios. In the current work, data outliers are not taken into consideration which are occasionally misleading and could distort the

(1) Pearson, K. On Lines and Planes of Closest Fit to Systems of Points in Space. Philos. Mag. 1901, 2, 559–572. (2) Hotelling, H. Analysis of a complex of statistical variables into principal components. J. Educ. Psych. 1933, 24, 417–441. (3) Chiang, L. H.; Russel, E.; Braatz, R. Fault Detection and Diagnosis in Industrial Systems; Springer: London, 2001. (4) Liu, X.; Chen, T.; Thornton, S. Eigenspace Updating for NonStationary Process and Its Application to Face Recognition. Pattern Recognit. 2003, 36, 1945–1959. (5) Malinowski, E. R. Factor Analysis in Chemistry; Wiley: New York, 1991. (6) Jackson, J. E. A User’s Guide to Principal Components; John Wiley & Sons, Inc: New York, 1991. (7) Raich, A. C.; C¸inar, A. Multivariate statistical methods for monitoring continuous processes: assessment of discrimination power of disturbance models and diagnosis of multiple disturbances. Chemom. Intell. Lab. Syst. 1995, 30, 37–48. (8) Dong, D.; McAvoy, T. J. Batch tracking via nonlinear principal component analysis. AIChE J. 1996, 42, 2199–2208. (9) Wold, S. Exponentially weighted moving principal component analysis and projection to latent structures. Chemom. Intell. Lab. Syst. 1994, 23, 149–161. (10) Li, W.; Yue, H.; Valle-Cervantes, S.; Qin, S. Recursive PCA for adaptive process monitoring. J. Process Control 2000, 10, 471–486. (11) Wang, X.; Kruger, U.; Lennox, B. Recursive partial least squares algorithms for monitoring complex industrial processes. Control Eng. Practice 2003, 11, 613–632. (12) Jin., H.; Lee, Y.-H.; Lee, G.; Han, C. Robust Recursive Principal Component Analysis Modeling for Adaptive Monitoring. Ind. Eng. Chem. Res. 2006, 45, 696–703. (13) 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. (14) He, X. B.; Yang, Y. P. Variable MWPCA for Adaptive Process Monitoring. Ind. Eng. Chem. Res. 2008, 47, 419–427. (15) Naik, A. S.; Yin, S.; Ding, S. X. Recursive identification algorithm for parity space based fault detection systems. Preprints of the 7th IFAC SafeProcess, Jun 30-Jul 3, 2009, Barcelona, Spain. (16) Yin, S.; Naik, A. S.; Ding, S. X. Adaptive process monitoring based on parity space methods. Preprints of the 7th IFAC SafeProcess, Jun 30Jul 3, 2009, Barcelona, Spain. (17) Wang, X.; Kruger, U.; Irwin, G. Process monitoring approach using fast moving window PCA. Ind. Eng. Chem. Res. 2005, 44, 5691–5702. (18) Liu, X.; Kruger, U.; Littler, T.; Xie, L.; Wang, S. Moving window kernel PCA for adaptive monitoring of nonlinear processes. Chemom. Intell. Lab. Syst. 2009, 96, 132–143. (19) Hall, P.; Marshall, D.; Martin, R. Incrementally Computing Eigenspace Models. Proc. Br. Machine Vision Conf. 1998, 286–295. (20) Hall, P.; Marshall, D.; Martin, R. Merging and Splitting Eigenspace Models. IEEE Trans. Pattern Anal. Machine Intell. 2000, 22, 1042–1049. (21) Hall, P.; Marshall, D.; Martin, R. Adding and subtracting eigenspaces with eigenvalue decomposition and singular value decomposition. Image Vision Comput. 2002, 20, 1009–1016. (22) Champagne, B. Adaptive Eigendecomposition of Data Covariance Matrices Based on First-Order Perturbations. IEEE Trans. Signal Process. 1994, 42, 2758–2770. (23) Willink, T. Efficient Adaptive SVD Algorithm for MIMO Applications. IEEE Trans. Signal Process. 2008, 56, 615–622.

Ind. Eng. Chem. Res., Vol. 49, No. 1, 2010 (24) Doukopoulos, X. G.; Moustakides, G. V. Fast and Stable Subspace Tracking. IEEE Trans. Signal Process. 2008, 56, 1452–1465. (25) Qin, S. J. Statistical process monitoring: basics and beyond. J. Chemom. 2003, 17, 480–502. (26) Stewart, G.; Sun, J. Matrix perturbation theory; Academic Press, San Diego, 1990. (27) Golub, G.; Loan., C. V. Matrix computations, 3rd ed.; Johns Hopkins Press Ltd.: London, 1996. (28) Himes, D.; Storer, R.; Georgakis, C. Determination of the number of principal components for disturbance detection and isolation. Proc. Am. Control Conf. 1994, 1279–1283. (29) Valle, S.; Li, W.; Qin, S. Selection of the Number of Principal Components: The Variance of the Reconstruction Error Criterion with a Comparison to other Methods. Ind. Eng. Chem. Res. 1999, 38, 4389–4401.

259

(30) Jackson, J. E.; Mudholkar, G. Control procedures for residuals associated with principal component analysis. J. Process Control 1979, 21, 341–349. (31) 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.

ReceiVed for reView May 4, 2009 ReVised manuscript receiVed October 14, 2009 Accepted October 21, 2009 IE900720W