Reconstruction-Based Contribution for Process Monitoring with Kernel

Ind. Eng. Chem. Res. 2010, 49, 7849–7857

7849

Reconstruction-Based Contribution for Process Monitoring with Kernel Principal Component Analysis Carlos F. Alcala and S. Joe Qin* Mork Family Department of Chemical Engineering and Materials Science, UniVersity of Southern California, Los Angeles, California

This paper presents a new method for fault diagnosis based on kernel principal component analysis (KPCA). The proposed method uses reconstruction-based contributions (RBC) to diagnose simple and complex faults in nonlinear principal component models based on KPCA. Similar to linear PCA, a combined index, based on the weighted combination of the Hotelling’s T2 and SPE indices, is proposed. Control limits for these fault detection indices are proposed using second-order moment approximation. The proposed fault detection and diagnosis scheme is tested with a simulated CSTR process where simple and complex faults are introduced. The simulation results show that the proposed fault detection and diagnosis methods are effective for KPCA. Introduction Principal component analysis (PCA) is a multivariate statistical tool widely used in industry for process monitoring.1-4 It decomposes the measurement space into a principal component space that contains the common-cause variability and a residual space that contains the process noise. The task of process monitoring is to detect and diagnose faults when they happen in a process. In process monitoring with PCA, fault detection is performed with fault detection indices. Popular indices are the Hotelling’s T2 statistic for monitoring the principal component subspace; the SPE index, which monitors the residual subspace; and a combination of both indices that monitors the whole measurement space. Once a fault is detected, it is necessary to diagnose its cause. A popular method for fault diagnosis is the contribution plots. They are based on the idea that variables with large contributions to a fault detection index are likely the cause of the fault. However, as was demonstrated by Alcala and Qin,5 even for simple sensor faults, contribution plots fail to guarantee correct diagnosis. As an alternative to contribution plots, Alcala and Qin5 propose the use of reconstruction-based contributions (RBC) that overcome the aforementioned shortcoming. The RBC method is related to fault identification by reconstruction,6,7 but it does not require the knowledge of fault directions. The PCA-based process monitoring scheme mentioned above is based on the assumption that the process behaves linearly. However, when a process is nonlinear, the monitoring of a process using a linear PCA model might not perform properly. To address this issue, several nonlinear PCA methods have been proposed.8-11 One of these methods is kernel principal component analysis (KPCA), developed by Scho¨lkopf et al., which maps measurements from their original space to a higher dimensional feature space where PCA is performed.11 KPCA has been applied successfully for process monitoring.12-15 The fault detection with KPCA is done in a similar way as with PCA using a T2, an SPE, and combined indices. However, fault diagnosis cannot be performed with contribution plots since there seems to be no way to calculate them with KPCA. Choi et al.13 propose a method for fault diagnosis employing the reconstruction-method of Dunia et al.,7 which looks at a fault identification index obtained when a fault detection index has * To whom correspondence should be addressed. E-mail: sqin@ usc.edu.

been reconstructed along the direction of a variable. While this method applies to sensor faults or a process fault with a known direction, it cannot be applied to diagnosing faults with unknown directions. In this paper, reconstruction-based contribution for kernel PCA (RBC-KPCA) is proposed for fault diagnosis. We present a matrix formulation for KPCA and the kernel trick that is simpler than the standard kernel PCA description, and is analogous to the linear PCA formulation familiar to the process monitoring community. The T2 and SPE fault detection indices are expressed in terms of a kernel vector mapped by the kernel function, rather than the intangible vector in the feature space. Furthermore, the combined index approach of Yue and Qin16 is proposed for fault detection in the feature space, and new control limits for the fault detection indices are proposed. The proposed RBC diagnosis method is used just like the standard contribution plots for KPCA-based fault diagnosis, although it is difficult to extend the standard contribution plots for KPCA. It is shown that simple and complex faults can be diagnosed correctly using RBC-KPCA with a simulation of a CSTR process. Simple faults are defined in Yoon and MacGregor17 as ones that do not propagate to other variables, while complex faults do affect many variables. The organization for the rest of this paper is as follows. The first section presents the KPCA algorithm in a simple matrix formulation. The second section presents fault detection indices expressed as quadratic forms of the kernel vector mapped by the kernel function. Control limits for these indices are derived on the basis of second-order moments in the feature space. The third section derives reconstruction-based contribution for the KPCA model for fault diagnosis. A CSTR simulation example is given in the fourth section, followed by conclusions in the last section. Kernel PCA KPCA Background. To make the kernel PCA model of process data with n variables, m measurements, taken under normal operating conditions, form the training set X ) [x1 x2 · · · xm]T. Each of the measurements xi is an n-dimensional vector and can be mapped into an h-dimensional space via a mapping function φi ) Φ(xi). The h-dimensional space is called the feature space and the n-dimensional space is called the measurement space. An important property of the feature space is that

10.1021/ie9018947  2010 American Chemical Society Published on Web 03/30/2010

7850

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

even though the original data is nonlinear in a bidimensional space, its map to a tridimensional space is linear. For process monitoring with linear PCA, the model is obtained from the training set X. In kernel PCA, the covariance matrix of a set of m mapped vectors φi is eigen-decomposed to obtain the principal loadings of the model. Assume that the vectors in the feature space are scaled to zero mean and form the training data as X ) [φ1 φ2 · · · φm]T. Let the sample covariance matrix of the data set in the feature space be S. We have, m

(m - 1)S ) X TX )

∑φφ

T i i

(3)

i)1

KPCA in the feature space is equivalent to solving the following eigenvector equation, m

X TXν ) Figure 1. Normal and faulty data in a bidimensional space.

∑ φ φ ν ) λν T i i

(4)

i)1

Note that φi is not explicitly defined nor is φiφjT. The so-called kernel trick premultiplies eq 4 by X : XX TXν ) λXν Defining

K ) XX ) T

[

φT1 φ1 · · · φT1 φm ·

l φmT φ1

··

l

···

φmT φm

]

[

)

]

· · · k(x1, xm) · ·· l l (5) k(xm, x1) · · · k(xm, xm) k(x1, x1)

and denoting r ) Xν

(6)

the dot product of two vectors φi and φj can be calculated as a function of the corresponding vectors xi and xj, that is,

Kr ) λr

(7)

φTi φj ) k(xi, xj)

Equation 7 shows that R and λ are an eigenvector and eigenvalue of K, respectively. In order to solve ν from eq 6, we premultiply it by X T and use eq 4,

Figure 2. Normal and faulty data in a tridimensional space.

we have

(1)

The function k( · , · ) is called the kernel function, and there exist several types of these functions. The Gaussian kernel function, or radial basis function, is expressed as

[

(xi - xj)T(xi - xj) k(xi, xj) ) exp c

]

X Tr ) X TXν ) λν which shows that ν is given by

(2)

and will be used in this paper. To illustrate how a nonlinear map to an expanded dimensional space can change a nonlinear distribution to a linear distribution, we give the following illustrative example. Suppose we have a nonlinear process with two variables, x1 and x2, and there are two data sets; one set has normal measurements and the other one faulty measurements. Figure 1 shows the plots of these data sets; the normal measurements are marked with + symbols and the faulty ones with dots. In this case it is impossible to apply linear PCA to separate the normal from the faulty data. However, as shown in Figure 2, if we add a third dimension to the plot, calculated as x12 + x22, it is really easy to separate the normal and faulty measurements with linear PCA. Therefore,

ν ) λ-1X Tr

(8)

Therefore, to calculate the PCA model (λi and νi), we first perform eigen-decomposition of eq 7 to obtain λi and Ri. Then use eq 8 to find νi. To guarantee that νiTνi ) 1, eqs 6 and 4 are used to derive rTi ri ) νTi X TXνi ) νTi λiνi ) λi Therefore, Ri needs to have a norm of √λi. Let R°i be the unit norm eigenvector corresponding to λi, ri ) √λiroi The matrix with the l leading eigenvectors are the KPCA principal loadings in the feature space, denoted as Pf ) [ν1 ν2


· · · νl]. From eq 8, Pf is related to the loadings in the measurement space as

[

1 T 1 X r1 · · · X T rl Pf ) λ1 λl [X Tro1λ1-1/2 T -1/2

)

···

]

(9)

X Trol λl-1/2]

) X PΛ

where P ) [R°1 · · · R°l ] and Λ ) diag{λ1 · · · λl} are the l principal eigenvectors and eigenvalues of K, respectively, corresponding to the largest eigenvalues in descending order. Another useful relationship that can be obtained from eq 6 is XPf ) [r1 · · · rl] ) [ro1λ11/2 · · · rol λl1/2] ) PΛ

Therefore, it is straightforward to transform vectors in the feature space to zero mean before KPCA. Without loss of generality, we assume the vectors in the feature space are zero mean for the rest of the paper. Fault Detection with KPCA Fault Detection Indices. Since the key idea of KPCA is to map the measurement space into the feature space, so that data in the feature space are linearly distributed, it is natural to perform fault detection by defining statistics in the feature space. The Hotelling’s T2 index is calculated as T2(x) ) tTΛ-1t, where Λ ) diag(λ1, · · · , λl) is the covariance of the scores t in the feature space. From eq 11, the T2 is calculated using kernel functions as T2 ) k(x)TPΛ-2PTk(x) ) k(x)TDk(x)

(10)

1/2

For a given measurement x and its mapped vector φ ) Φ(x), the scores are calculated as t ) PfTφ, which, from eq 9, can be expressed as t ) Λ-1/2PTXφ ) Λ-1/2PTk(x)

(11)

7851

(17)

with D ) PΛ-2PT. The SPE index is defined as the norm of the residual vector ˜ fφ, where C ˜ f is the in the feature space and is calculated as φTC projection matrix of the residual space, which is orthogonal to the principal component space. Let ˜t be the residual components and P˜f the corresponding loading matrix, t˜ ) P˜Tf φ ) [νl+1 νl+2 · · · ]Tφ

where k(x) ) Xφ ) [φ1 φ2 · · · φm]Tφ )

[φT1 φ

φT2 φ

···

φmT φ]T

(12)

) [k(x1, x) k(x2, x) · · · k(xm, x)]

SPE ) t˜Tt˜ ) φTP˜fP˜Tf φ

T

Scaling. The calculation of the covariance matrix in eq 3 holds if the vector φ in the feature space has zero mean. If this is not the case, the vectors in the feature space have to be scaled to zero mean using the sample mean of the training data. The j is scaled vector φ 1 φ¯ ) φ m

The SPE index is calculated as the squared norm of the residual components

Since we do not know the dimension of the feature space, it is not possible to know the number of residual components there. Thus, we cannot calculate explicitly the loading matrix P˜f. However, we can calculate the product P˜fP˜fT as the projection orthogonal to the principal component space, which is ˜ f ) P˜fP˜Tf ) I - PfPTf C

m

∑φ

i

) φ - [φ1 · · · φm]1m and leads to

i

where 1m is an m-dimensional vector whose elements are 1/m. j i and φ j j is The kernel function of two scaled vectors φ

SPE ) φT(I - PfPTf )φ ) φTφ - φTPfPTf φ From eqs 1 and 9 we have

jk(xi, xj) ) φ¯ Ti φ¯ j ) k(xi, xj) - k(xi)T1m - k(xj)T1m · · ·

(13)

+ 1mT K1m Similarly, the scaling of the kernel vector k(x) is k¯(x) ) [φ¯ 1 φ¯ 2 · · · φ¯ m] φ ) F[k(x) - K1m]

(14)

F)I-E

(15)

T¯

where

In this equation, I is the identity matrix, and E is an m × m matrix with elements 1/m. Finally, the scaled kernel matrix K, j , is calculated as K ¯ ) [φ¯ 1 φ¯ 2 · · · φ¯ m]T[φ¯ 1 φ¯ 2 · · · φ¯ m] K ) FKF

(16)

SPE ) k(x, x) - φTXTPΛ-1PTXφ ) k(x, x) - k(x)PΛ-1PTk(x) ) k(x, x) - k(x)TCk(x)

(18)

where C ) PΛ-1PT. Yue and Qin16 proposed the use of a combined index for monitoring the principal and residual space simultaneously. Such an index is a combination of the T2 and SPE indices weighted by their control limits. The same concept is used here to define a fault detection index to monitor the principal and residual spaces in the feature space. A combined index for fault detection in the feature space has been proposed by Choi et al.;13 however, their definition is different to the one proposed here as they use an energy approximation concept to calculate the index. The combined index proposed for the feature space is defined as φ(x) )

SPE(x) T2(x) + 2 2 δ τ

(19)

7852


Table 1. Control Limits for the Fault Detection Indices. Limit ) gIndexχr(hIndex) Index

limit

2

2

T

τ

SPE

δ2

φ

ζ2 2

parameters T2

g 2 ) 1/(m - 1) hT ) l m gSPE ) Σi)l+1 λi2/[(m - 1)Σim)l+1λi] m m hSPE ) (Σi)l+1 λi2) λi)2/(Σi)l+1 m m gφ ) (l/τ4 + Σi)l+1 λi2/δ4/[(m - 1)(l/τ2 + Σi)l+1 λi/δ2)] m 2 4 m 2 2 4 φ 2 h ) (l/τ + Σi)l+1λi/δ ) /[l/τ + Σi)l+1λi /δ ] 2

2

where δ and τ are the control limits for the SPE and T indices, respectively. The calculation of these control limits is provided in the next section. The combined index can be calculated with kernel functions as

[

]

PΛ-2PT PΛ-1PT k(x, x) + k(x)T k(x) φ(x) ) 2 2 δ τ δ2 (20) k(x, x) T + k(x) Ωk(x) ) δ2 where Ω)

PΛ-2PT PΛ-1PT D C ) 2 - 2 2 τ δ2 τ δ

(21)

The combined index has a control limit ζ2 and its calculation is provided in the next section. Like linear PCA, one needs to select the number of principal components in the KPCA model as well. The selection methods in Qin and Dunia18 for linear PCA can be applied to KPCA in a straightforward manner. When the combined index is used for fault detection, the selection of the number of principal components is not as critical as using SPE or T2 alone, since the entire feature space is monitored. Control Limits. In the feature space, the fault detection indices have the quadratic form J ) xTAx

(22)

with A g 0. If it is considered that x is zero mean and normally distributed, the results of Box19 can be applied to calculate control limits for the fault detection indices. These limits can be calculated as gIndexχR2 (hIndex), with a confidence level of (1 - R) × 100%, and the parameters gIndex and hIndex calculated as gIndex ) hIndex )

the direction of a variable, but along an arbitrary direction. Therefore, it is possible to obtain the contribution of simple faults, as well as complex faults provided that their directions are available in a fault library. The objective of RBC is to find the magnitude fi of a vector with direction ξi such that the fault detection index of the reconstructed measurement zi ) x - ξifi

is minimized; that is, we want to find fi ) arg min Index(x ξifi). The same concept can be applied to KPCA and find fi such that fi ) arg min Index(k(x - ξifi))

(23)

[tr{SA}]2 tr{SA}2

(24)

where tr{A} is the trace of matrix A. We use Index to represent either of the SPE, T2, or φ indices. The parameters gIndex and hIndex for the detection indices are shown in Table 1, and the details of their calculation are given in Appendix A. In Table 1, the parameter λi is the ith eigenvalue of the kernel matrix K. Reconstruction-Based Contribution for Fault Diagnosis Reconstruction-based contribution (RBC), proposed by Alcala and Qin,5 defines the reconstruction of a fault detection index along the direction of a variable as the variables’ contribution for fault diagnosis. The variable with the largest amount of reconstruction is likely a major contributor to the fault. This method is useful for it cannot only obtain contributions along

(26)

The RBC value for the direction ξi is fi2. If we want to do reconstruction along the direction of a variable, the direction can be written as ξi )[0 0 · · · 1 · · · 0], where 1 is placed at the ith position. When complex faults are reconstructed, the direction of the fault can be any unitary vector ξi. To find the RBC along a direction ξi for a fault detection index, we have to perform a nonlinear search of the value of fi that minimizes Index(k(x - ξif)). This is done by obtaining the first derivative of the fault detection index with respect to fi, then equating it to zero, and solving for fi. However, the resulting expression is not an explicit solution for fi, and it has to be iterated until fi has converged. With the use of the kernel function in eq 2 and eq 17, the derived expression of fi for the T2 index is fi )

ξTi BTFDk¯(zi)

(27)

kT(zi)FDk¯(zi)

In this equation, F is the scaling matrix defined in eq 15, k(zi) j (zi) is the vector with the unscaled values of k(zi, xi), and k contains the scaled values and is calculated with eq 14. The matrix B is calculated as

B)

[

k(zi, x1)(x - x1)T k(zi, x2)(x - x2)T l k(zi, xm)(x - xm)T

]

(28)

For SPE, via eq 18, the derived expression for fi is

2

tr{SA} tr{SA}

(25)

fi )

ξTi BT[1m + FCk¯(zi)] kT(zi)[1m + FCk¯(zi)]

(29)

For the index φ, using eq 20, fi is calculated as ξTi BT fi ) T

[ [

k (zi)

1m δ2 1m δ2

] ]

- FΩk¯(zi)

(30)

- FΩk¯(zi)

A detailed derivation of the RBC values for the fault detection indices is given in Appendix B. In the next section a simulation is performed to test the fault detection and diagnosis methods discussed so far. CSTR Simulation In this section, a nonisothermal continuous stirred tank reactor (CSTR) is simulated. The process model and simulation


7853

conditions are similar to the ones provided by Yoon and MacGregor,17 and they are given in Appendix C; a diagram of the proces is shown in Figure 3. The simulation is performed according to the following model: dCA F F ) CA0 - CA - k0e-E/(RT)CA dt V V VFcp

(31)

aFCb+1 dT (T - TC,in) + ) FcpF(T0 - T) dt FC + aFCb /(2FCcpc) (-∆Hrxn)Vk0e-E/(RT)CA(32)

The process is monitored measuring the cooling water temperature TC, the inlet temperature T0, the inlet concentrations CAA and CAS, the solvent flow FS, the cooling water flow FC, the outlet concentration CA, the temperature T, and the reactant flow FA. These nine variables form the measurement vector x ) [TC T0 CAA CAS FS FC CA T FA]T

Figure 4. Fault detection indices and RBC values for Fault 1.

(33)

The variables are sampled every minute and 100 samples, taken under normal conditions, are used as the training set. Five different faults are introduced into the system. The first four faults are sensor faults similar to the ones shown in Choi et al.13 The first simulated fault, Fault 1, is a bias in the sensor of the output temperature T; the bias magnitude is 1 (K). Since T is a controlled variable, the effect of the fault will be removed by the PI controller, and its effect will propagate to other variables. This fault is considered a complex fault since it affects several variables. The second and third faults, Faults 2 and 3, are biases in the sensors of the inlet temperature To and inlet reactant concentration CAA, respectively; the bias magnitude for T0 is 1.5 (K) and for CAA is 1.0 (kmol/ m3). These faults are considered simple faults since they only affect one variable. The fourth fault, Fault 4, is a drift in the sensor of CAA and its magnitude is dCAA/dt ) 0.2 (kmol/(m3 min)); this is also a simple fault. The last fault, Fault 5, is a slow drift in the reaction kinetics. The fault has the form of an exponential degradation of the reaction rate caused by catalyst poisoning; in this case, the reaction rate coefficient will change with time as k0(t + 1) ) 0.996 × k0(t). This process fault is a complex fault that affects several variables such as the output temperature T, concentration CA, and the cooling water flow FC. For each of the five scenarios mentioned above, 100 measurements are simulated, and the fault is introduced at the 51st measurement. The KPCA model is built using four principal components, and the parameter c in the kernel equation is set to 0.65. Figures 4-8 show the fault detection and diagnosis results for the five faults. The first column in each set of plots shows



Figure 3. Diagram of the CSTR process.17

the SPE, T2, and φ indices for the simulated measurements. The second column shows the results of fault diagnosis given by RBC with KPCA for the three fault detection indices. In Figure 4, it is shown that the SPE and combined indices detect the fault as soon as it happens, not doing the same as the T2 index. The first nine bars of the RBC plots, for the three fault detection indices, show the RBC values of all variables. Since the bias fault happens in the temperature sensor, variable 8, this variable is supposed to have the largest RBC; however, due to the feedback control, the effect of the fault is transferred to the

7854


they become very large compared to their control limits; this is not the case of T2 since it does not change much. Looking only at the first nine bars of the diagnosis results, they would point to the cooling water flow, variable 6, as the origin of the fault. However, if we have information of this fault from past experience, we can calculate RBC values along its direction; the average of these values is shown in the last bar of the plot. It can be seen that the RBC values are the largest ones along this direction; thus, we can diagnose this fault as a drift in the kinetics. Conclusions



cooling water flow, variable 6, as indicated by the RBCs of SPE, T2, and φ. In this case, more information is needed to determine the cause of the fault. Once the cause of the fault has been determined and stored in a library of known faults, its direction can be used to calculate the RBC value along this fault direction. The tenth bar in the right column of Figure shows the RBC value for the direction of Fault 1. We can see that this fault has the largest RBC value with the three fault detection indices. The detection and diagnosis results for Fault 2 are shown in Figure 5. Since this is a bias fault in T0, variable 2, it is expected that this variable has the largest RBC value, which actually happens for the SPE and φ indices. Again, the T2 index does not detect nor diagnose this fault correctly. Figure 6 shows the results for the fault in the sensor CAA, variable 3. In this case all detection indices show CAA as the cause of the problem; however, the SPE and combined indices make a better detection of the fault. Figure 7 shows the results for the drift in the sensor of CAA, variable 3. In this case, only the SPE and combined indices detect the fault right after it is introduced in the system; however, even though the fault magnitude becomes very large, the T2 index only detects the first few faulty measurements. On the other hand, all three detection indices diagnose the fault correctly. In Figure 8 we can see the detection and diagnosis results for the slow drift on the reaction kinetics, Fault 5. In this case the SPE and combined indices start growing up slowly until

It is shown in this paper that a reconstruction-based contribution can be used with kernel PCA for the diagnosis of faults in nonlinear processes. The proposed RBC method works like the standard contribution plots in linear PCA, which do not require the fault direction to be known beforehand. Further, contributions along sensor directions as well as known process fault directions can be obtained with RBC. The combined index is proposed to monitor the whole feature space, and new control limits are proposed for the three fault detection indices in the feature space. The detection power and diagnosis power of the proposed methods have been tested with a simulation of a CSTR process with sensor and process faults. It was shown that, while the SPE and T2 indices sometimes disagreed on their results, the combined index was effective in correctly diagnosing simple and complex faults. Several open issues deserve further study. One is the high dimensionality of the feature space when the number of samples is very large. Another issue is the selection of the radius of the kernel function. This parameter, while being determined by trial and error at the moment, affects the nonlinearity of the mapping from the data space to the feature space. Acknowledgment We appreciate the financial support from the Roberto Rocca Education Program and the Texas-Wisconsin-California Control Consortium. Appendix A. Calculation of Control Limits for Fault Detection Indices. Control Limit for the T2 Index. In the third section, the T2 index was calculated with the scores in the feature space. Using mapped vectors, the index can be calculated as T2 ) φTPfΛ-1PfTφ, where Pf is the loading matrix in the feature space. Since the covariance matrix (m - 1)S can be decomposed as

[ ]

Λ 0 [Pf P˜f]T ˜ 0 Λ ˜ P˜Tf ) PfΛPTf + P˜fΛ

(m - 1)S ) [Pf P˜f]

we have ˜ P˜Tf )PfΛ-1PTf (m - 1)SA ) (PfΛPTf + P˜fΛ ) PfΛPTf PfΛ-1Pf ) PfPTf and (m - 1)tr{SA} ) tr{PfPTf } ) tr{PTf Pf} ) tr{I} )l

(A.1)


In the same way, with eq A.1, we can calculate the trace of {SA}2 as hSPE )

T T (m - 1) tr{SA} ) tr{PfPf PfPf } ) tr{PfPTf } ) tr{PTf Pf} ) tr{I} )l 2

2

2

2

(tr{SA})2 )l tr{SA}2

(A.3)

2

1 χ2 (l) m-1 R

(A.4)

(m - 1)SA ) )

+

PfΛ-1PTf τ2

P˜fP˜Tf

+

δ2

The matrix Φf is equivalent to the matrix A in eq 22; thus, SA can be calculated as (m - 1)SA ) (m - 1)SΦf

Control Limit for the SPE Index. Using mapped vectors, the SPE index can be calculated in the feature space as SPE ) φTP˜fP˜fTφ ) φT(I - PfPfT)φ. Therefore, the matrix A of the quadratic form in eq 22 is A ) P˜fP˜fT ) I - PfPfT, and the term SA can be calculated as (PfΛPTf ˜ P˜Tf P˜fΛ

∑λ

where

and the control limit 2

(A.6) 2

φ ) φTΦfO

Φf )

τ2 ) gT χ2R(hT ) )

i)l+1 m

Control Limit for the Combined Index. From the definitions of the T2 and SPE indices in the feature space, the combined index can be calculated with mapped vectors as

2

(A.2)

2

(tr{SA}) ) tr{SA}2

2

i

i

1 tr{SA}2 ) tr{SA} m-1

hT )

∑ λ)

(

2

i)l+1

Now, we can calculate the parameters gT and hT as gT )

7855

m

(

˜ P˜Tf ) ) (PfΛPTf + P˜fΛ )

PfPTf 2 τ

˜ P˜Tf )P˜fP˜Tf P˜fΛ

+

˜ P˜Tf P˜fΛ 2

PfΛ-1PTf

+

τ2

P˜fP˜Tf δ2

)

δ

and m

and l (m - 1)tr{SA} ) 2 + τ

h

∑λ

˜ P˜Tf P˜f} ) tr{Λ ˜} ) (m - 1)tr{SA} ) tr{Λ

i

∑λ

i

i)l+1 2

δ

i)l+1

Further, The use of Gaussian kernel functions makes h, the dimension of the feature space and size of the covariance matrix S, to be infinite. Thus, the number of residual eigenvalues λ˜ is infinite, too. However, since S is calculated with m mapped vectors φ, and m is smaller than h, then the rank of S, and number of nonzero eigenvalues λi, can be up to m. Therefore,

(m - 1)2{SA}2 )

PfPTf τ4

+

˜ 2P˜Tf P˜fΛ δ4

and we have m

m

(m - 1)tr{SA} )

∑λ

l (m - 1) tr{SA} ) 4 + τ

i

2

l+1

Noting that

2

∑λ

2 i

i)l+1 4

δ

The parameters of the control limits are (m - 1) tr{SA} ) 2

2

˜ 2P˜Tf P˜fΛ gφ )

m l/τ4 + Σi)l+1 λ2i /δ4 tr{SA}2 ) m tr{SA} (m - 1)(l/τ2 + Σi)l+1 λi /δ2)

hφ )

m (l/τ2 + Σi)l+1 λi /δ2)2 (tr{SA})2 ) m tr{SA}2 l/τ4 + Σi)l+1 λ2i /δ4

we have m

˜ 2} ) (m - 1)2tr{SA}2 ) tr{Λ

∑λ

2 i

l+1

and the limit is

Therefore, m

gSPE

tr{SA}2 ) ) tr{SA}

∑ λ /(m - 1) 2 i

i)l+1 m

2

∑ λ /(m - 1) i

i)l+1

ζ2 ) gφχR2(hφ)

m

)

∑λ

2 i

i)l+1

m

(m - 1)

∑λ

i

i)l+1

(A.5)

B. Calculation of the RBC of Fault Detection Indices. For the reconstructed measurement given in eq 25, we want to find the value of fi that minimizes Index(k(zi)). We can perform this calculation for each of the fault detection indices. For the sake of simplicity, we define a general form for the detection indices

7856


[ [

Table 2. Parameters for Fault Detection Indices in Kernel Form Index

a

M

2

0 1 1/ δ2

D -C D/τ2 - C/δ2

T SPE φ

calculated with kernel functions. From eq 17, eq 18, and eq 20 we can form a general expression for Index as Index ) akj(zi, zi) + k¯(zi)TMk¯(zi)

]

T ∂k(zi) 2 k(zi, x2)(zi - x2) ) ξi ∂fi c l k(zi, xm)(zi - xm)T . k(zi, x1)(x - x1)T - k(zi, x1)fiξTi

.

(B.2) B)

Since the kernel function is scaled as jk(zi, zi) ) 1 - 2kT(zi)1m + 1mT K1m

[

k(zi, x1)(x - x1)T k(zi, x2)(x - x2)T l k(zi, xm)(x - xm)T

]

We can now calculate the derivative of Index as ∂Index 4 ) - [a1m - FMk¯(zi)]T[Bξi - k(zi)fi] ∂fi c

we have that ∂kj(zi, zi) ∂k(zi) ) -21mT ∂fi ∂fi

(B.3)

j (zi) ) F[k(zi) - K1m], where Also, the scaled kernel vector is k F is given by eq 15. Then, we have that ∂k¯(zi) ∂k(zi) )F ∂fi ∂fi

(B.4)

From eqs B.2-B.4, the derivative of Index turns out to be ∂k(zi) ∂Index ) -2[a1m - FMk¯(zi)]T ∂fi ∂fi To calculate the derivative of k(zi) with respect to fi, we know that ∂k(zi) ∂zi ∂k(zi) ) ∂fi ∂zi ∂fi Since k(zi) ) [k(zi, x1) k(zi, x2) · · · k(zi, xm)]T and k(zi, xj) ) exp(-(zi - xj)T(zi - xj)/c), we have (zi - xj)T ∂ k(zi, xj) ) -2k(zi, xj) ∂zi c and ∂zi ) -ξi ∂fi Therefore, the vector with the derivative of k(zi) respect to fi is

(B.5)

2 ) [Bξi - k(zi)fi] c where B is calculated as

∂kj(zi, zi) ∂k¯(zi) ∂Index )a + 2k¯T(zi)M ∂fi ∂fi ∂fi

]

T T 2 k(zi, x2)(x - x2) - k(zi, x2)fiξi ) ξi c l k(zi, xm)(x - xm)T - k(zi, xm)fiξTi

(B.1)

where the parameters a and M for each detection index are shown in Table 2. Now, the derivative of Index with respect to fi is

k(zi, x1)(zi - x1)T

After setting the derivative equal to zero and solving for fi we obtain fi )

ξTi BT[a1m - FMk¯(zi)] kT(zi)[a1m - FMk¯(zi)]

(B.6)

C. CSTR Simulation Parameters and Initial Conditions. The CSTR simulation parameters are V ) 1 (m3), F ) 106 (g/ m3), FC ) 106 (g/ m3), E/R ) 8333.1 (K), cp ) 1 (cal/(g K)), cpc ) 1(cal/(g K)), b ) 0.5, k0 ) 1010 (m3 /kmole × min), a ) 1.678 × 106 (cal/(min K)), and ∆Hrxn ) -1.3 × 107 (cal/kmole). The parameters of the temperature PI controller are KC ) -1.5 and TI ) 5.0. The initial conditions for the simulation are T0 ) 370.0 (K), TC ) 365.0 (K), FC ) 15 (m3/min), T ) 368.25 (K), Fs ) 0.9 (m3/min), FA ) 0.1 (m3/min), CA ) 0.8 (kmol/ m3), CAS ) 0.1 (kmol/ m3), CAA ) 19.1 (kmol/ m3). Gaussian noise with small variance was added to the measurements and disturbances. For more details of the simulation check Yoon and MacGregor.17 Literature Cited (1) Nomikos, P.; MacGregor, J. Multivariate SPC charts for monitoring batch processes. Technometrics 1995, 37, 41–59. (2) Wise, B.; Gallagher, N. The process chemometrics approach to process monitoring and fault detection. J. Process Control 1996, 6, 329– 348. (3) Qin, S. J. Statistical process monitoring: Basics and beyond. J. Chemom. 2003, 17, 480–502. (4) Chiang, L.; Russell, E.; Braatz, R. Fault diagnosis and Fisher discriminant analysis, discriminant partial least squares, and principal component analysis. Chemom. Intell. Lab. Syst. 2000, 50, 243–252. (5) Alcala, C. F.; Qin, S. J. Reconstruction-based contribution for process monitoring. Automatica 2009, 45, 1593–1600. (6) Dunia, R.; Qin, S. J. A unified geometric approach to process and sensor fault identification: The unidimensional fault case. Comput. Chem. Eng. 1998, 22, 927–943. (7) Dunia, R.; Qin, S. J.; Edgar, T. F.; McAvoy, T. J. Identification of faulty sensors using principal component analysis. AIChE J. 1996, 42, 2797–2812.

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010 (8) Kramer, M. Nonlinear principal component analysis using autoassociative neural networks. AIChE J. 1991, 37, 233–243. (9) Dong, D.; McAvoy, T. J. Nonlinear principal component analysis based on principal curves and neural networks. Proc. Am. Control Conf. 1994, 1284–1288. (10) Hiden, H. G.; Willis, M. J.; Tham, M. T.; Turner, P.; Montague, G. A. Nonlinear Principal Components Analysis Using Genetic Programming. Genetic Algorithms in Engineering Systems: Innovations and Applications, Glasgow, UK, 1997. (11) Scho¨lkopf, B.; Smola, A.; Mu¨ller, K. Nonlinear component analysis as a kernel eigenvalue problem. Neural Comput. 1998, 10, 1299–1319. (12) Lee, J.; Yoo, C.; Lee, I. Fault detection of batch processes using multiway kernel principal component analysis. Comput. Chem. Eng. 2004, 28, 1837–1847. (13) Choi, S. W.; Lee, C.; Lee, J.; Park, J. H.; Lee, I. Fault detection and identification of nonlinear processes based on kernel PCA. Chemom. Intell. Lab. Syst. 2005, 75, 55–67. (14) Zhang, Y.; Qin, S. J. Fault detection of nonlinear processes using multiway kernel independent analysis. Ind. Eng. Chem. Res. 2007, 46, 7780– 7787.

7857

(15) Zhang, Y.; Qin, S. J. Improved nonlinear fault detection technique and statistical analysis. AIChE J. 2008, 54, 3207–3220. (16) Yue, H.; Qin, S. J. Reconstruction-based fault identification using a combined index. Ind. Eng. Chem. Res. 2001, 40, 4403–4414. (17) Yoon, S.; MacGregor, J. Fault diagnosis with multivariate statistical models, Part I: Using steady-state fault signatures. J. Process Control 2001, 11, 387–400. (18) Qin, S. J.; Dunia, R. Determining the number of principal components for best reconstruction. J. Process Control 2000, 10, 245– 250. (19) Box, G. Some theorems on quadratic forms applied in the study of analysis of variance problems, I. Effect of inequality of variance in the one-way classification. Ann. Math. Stat. 1954, 25, 290–302.

ReceiVed for reView December 1, 2009 ReVised manuscript receiVed February 22, 2010 Accepted February 24, 2010 IE9018947

Reconstruction-Based Contribution for Process Monitoring with Kernel

Recommend Documents