Nonlinear process fault detection and identification using kernel PCA and kernel density estimation

ABSTRACT Kernel principal component analysis (KPCA) is an effective and efficient technique for monitoring nonlinear processes. However, associating it with upper control limits (UCLs) based on the Gaussian distribution can deteriorate its performance. In this paper, the kernel density estimation (KDE) technique was used to estimate UCLs for KPCA-based nonlinear process monitoring. The monitoring performance of the resulting KPCA–KDE approach was then compared with KPCA, whose UCLs were based on the Gaussian distribution. Tests on the Tennessee Eastman process show that KPCA–KDE is more robust and provide better overall performance than KPCA with Gaussian assumption-based UCLs in both sensitivity and detection time. An efficient KPCA-KDE-based fault identification approach using complex step differentiation is also proposed.


Introduction
There has been an increasing interest in multivariate statistical process monitoring methods in both academic research and industrial applications in the last 25 years (Chiang, Russell, & Braatz, 2001;Ge, Song, & Gao, 2013;Russell, Chiang, & Braatz, 2000;Yin, Ding, Xie, & Luo, 2014;Yin, Li, Gao, & Kaynak, 2015). Principal component analysis (PCA) (Jolliffe, 2002;Wold, Esbensen, & Geladi, 1987) is probably the most popular among these techniques. PCA is capable of compressing high-dimensional data with little loss of information by projecting the data onto a lower-dimensional subspace defined by a new set of derived variables (principal components (PCs)) (Wise & Gallagher, 1996). This also addresses the problem of dependency between the original process variables. However, PCA is a linear technique; therefore, it does not consider or reveal nonlinearities inherent in many real industrial processes (Lee, Yoo, Choi, Vanrolleghem, & Lee, 2004). Hence, its performance is degraded when it is applied to processes that exhibit significant nonlinear variable correlations.
To address the nonlinearity problem, Kramer (1992) proposed a nonlinear PCA based on auto-associative neural network (NN). Dong & McAvoy (1996) suggested a nonlinear PCA that combined principal curve and NN. Their approach involved: (i) using principal curve method to obtain associated scores and the correlated data, CONTACT Yi Cao y.cao@cranfield.ac.uk (ii) using an NN model to map the original data into scores, and (iii) mapping the scores into the original variables. Nonlinear PCA methods have also been proposed by (Cheng & Chiu, 2005;Hiden, Willis, Tham, & Montague, 1999;Jia, Martin, & Morris, 2000;Kruger, Antory, Hahn, Irwin, & McCullough, 2005). However, most of these methods are based on NNs and require the solution of a nonlinear optimization problem. Scholkopf, Smola, & Muller (1998) proposed Kernel PCA (KPCA) as a nonlinear generalization of the PCA. A number of studies adopting the technique for nonlinear process monitoring have also been reported in the literature (Cho, Lee, Choi, Lee, & Lee, 2005;Choi, Lee, Lee, Park, & Lee, 2005;. KPCA is performed in two steps: (i) mapping the input data into a high-dimensional feature space, and (ii) performing standard PCA in the feature space. Usually, high-dimensional mapping can seriously increase computational time. This difficulty is addressed in kernel methods by defining inner products of the mapped data points in the feature space, then expressing the algorithm in a way that needs only the values of the inner products. Computation of the inner products in the feature space is then done implicitly in the input space by choosing a kernel that corresponds to an inner product in the feature space. Unlike NN-based methods, KPCA does not involve solving a nonlinear optimization problem; it only solves an eigenvalue problem.
In addition, KPCA does not require specifying the number of PCs to extract before building the model .
Similar to the PCA, the Hotelling's T 2 statistic and the Q statistic (also known as squared prediction error, SPE) are two indices commonly used in KPCA-based process monitoring. The T 2 is used to monitor variations in the model space while the Q statistic is used to monitor variations in the residual space. In a linear method such as PCA, computation of the upper control limits (UCLs) of T 2 and Q statistics is based on the assumption that random variables included in the data are Gaussian. The actual distribution of T 2 and Q statistics can be analytically derived based on this assumption. Hence, the UCLs can also be derived analytically. However, many real industrial processes are nonlinear. Even though the sources of randomness of these processes could be assumed as Gaussian, variables included in measured data are non-Gaussian due to inherent nonlinearities. Hence, adopting UCLs for fault detection based on the multivariate Gaussian assumption in such processes is inappropriate and may lead to misleading results .
An alternative solution to the non-Gaussian problem is to derive the UCLs from the underlying probability density functions (PDFs) estimated directly from the T 2 and the Q statistics via a non-parametric technique such as kernel density estimation (KDE). This approach has been suggested in various linear techniques, such as PCA (Chen, Wynne, Goulding, & Sandoz, 2000;Liang, 2005), independent component analysis (Xiong, Liang, & Qian, 2007), and canonical variate analysis (Odiowei & Cao, 2010). It is even more important to adopt this kind of approach to derive meaningful UCLs for a nonlinear technique such as the KPCA. This is because the Gaussian-assumption-based UCLs for latent variables obtained through a nonlinear technique will not be valid at all. Unfortunately, this issue has not attracted enough attention in the literature. In this work, the KDE approach is adopted to derive UCLs for PCA and KPCA. Then, their fault detection performances in the Tennessee Eastman (TE) process are compared with their Gaussianassumption-based equivalents. The results show that the KDE-based approaches perform better than their counterparts with Gaussian-assumption-based UCLs. The contributions of this work can be summarized as follows: • To combine the KPCA with KDE for the first time to show that it is not appropriate to use the Gaussianassumption-based UCLs with a nonlinear approach such as the KPCA. • To compare the robustness of KPCA-KDE and KPCA associated with Gaussian-assumption-based UCLs.
• Propose an efficient KPCA-KDE-based fault identification approach using complex step differentiation.
The paper is organized as follows. The KPCA algorithm, KPCA-based fault detection and identification, and the KDE technique are discussed in Section 2. Application of the monitoring approaches to the TE benchmark process is presented in Section 3. Finally, conclusions drawn from the study are given in Section 4.

Kernel PCA algorithm
Given m training samples x k ∈ n , k = 1, . . . , m, the data can be projected onto a high-dimensional feature space using a nonlinear mapping, φ : The covariance matrix in the feature space is then computed as where φ(x j ), for j = 1, . . . m is assumed to have zero mean and unit variance. To diagonalize the covariance matrix, we solve the eigenvalue problem in the feature space as where λ is an eigenvalue of C F , satisfying λ ≥ 0, and a ∈ F is the corresponding eigenvector (a = 0). The eigenvector can be expressed as a linear combination of the mapped data points as follows: Using φ(x k ) to multiply both sides of Equation (2) gives Substituting Equations (1) and (3) in Equation (4) we have Instead of performing eigenvalue decomposition directly on C F in Equation (1) and finding eigenvalues and PCs, we apply the kernel trick by defining an m × m (kernel) matrix as follows: for all i, j = 1, . . . , m. Introducing the kernel function of the form k(x, y) = φ(x), φ(y) in Equation (5) enables the computation of the inner products φ(x i ), φ(x j ) in the feature space as a function of the input data. This precludes the need to carry out the nonlinear mappings and the explicit computation of inner products in the feature space Scholkopf et al., 1998). Applying the kernel matrix, we re-write Equation (5) as Notice that k = 1, . . . , m, and therefore, Equation (7) can be represented as Equation (8) is equivalent to the eigenvalue problem Furthermore, the kernel matrix can be mean centred as follows: where U is an m × m matrix in which each element is equal to 1/m. Eigen decomposition of K ctr is equivalent to performing PCA in F . This, essentially, amounts to resolving the eigenvalue problem in Equation (9), which yields eigenvectors α 1 , α 2 , . . . , α m and the corresponding eigenvalues Since the kernel matrix, K ctr is symmetric, the derived PCs are orthonormal, that is, where δ i,j represents the Dirac delta function. The score vectors of the nonlinear mapping of meancentred training observations x j , j = 1, . . . , m, can then be extracted by projecting φ(x j ) onto the PC space spanned by the eigenvectors α k , k = 1, . . . , m, Applying the kernel trick, this can be expressed as

Fault detection metrics
The Hotelling's T 2 of the jth samples in the feature space used for KPCA fault detection is computed as where z i,j , i = 1, . . . , q represents the PC scores of the jth samples, q is the number of PCs retained and −1 represents the inverse of the matrix of eigenvalues corresponding to the retained PCs. The control limit of T 2 can be estimated from its distribution. If all scores are of Gaussian distributions, then the control limit corresponding to a significance level, α, T 2 α can be derived from the F-distribution analytically as F q,m−q,α is the value of the F-distribution corresponding to a significance level, α with degrees of freedom q and m−q for the numerator and denominator, respectively. Furthermore, a simplified computation of the Qstatistic has been proposed . For the jth samples, If all scores are of normal distributions, the control limit of the Q-statistic at the 100(1 − α)% confidence level can be derived as follows (Jackson, 1991): where The alternative method of computing control limits directly from the PDFs of the T 2 and Q statistics is explained in the next section.

Kernel density estimation
KDE is a procedure for fitting a data set with a suitable smooth PDF from a set of random samples. It is used widely for estimating PDFs, especially for univariate random data (Bowman & Azzalini, 1977). The KDE is applicable for the T 2 and Q statistics since both are univariate although the process characterized by these statistics is multivariate.
Given a random variable y, its PDF g(y) can be estimated from its m samples, y j , j = 1, . . . , m, as follows: where K is a kernel function while h is the bandwidth or smoothing parameter. The importance of selecting the bandwidth and methods of obtaining an optimum value are documented in Chen et al. (2000) and Liang (2005). Integrating the density function over a continuous range gives the probability. Thus, assuming the PDF g(y), the probability of y to be less than c at a specified significance level, α is given by Consequently, the control limits of the monitoring statistics (T 2 and Q) can be calculated from their respective PDF estimates:

On-line monitoring
For a mean-centred test observation, x tt , the corresponding kernel vector, k tt is calculated with the training samples, x j , j = 1, . . . , m as follows: The test kernel vector is then centred as shown below: The corresponding test score vector, z tt is calculated using This can be re-written as In vector form, where A = [α 1 , · · · , α m ].

Outline of KPCA-KDE fault detection procedure
Tables 1 and 2 show the outline of KPCA-KDE-based fault detection procedure.
To provide a more intuitive picture, a flowchart of the procedure is presented in Figure 1.

Fault variable identification
After a fault has been detected, it is important that the variables most strongly associated with the fault are identified in order to facilitate the location of root causes. TR1. Obtain data under normal operating conditions (NOC) and scale the data using the mean and standard deviation of the columns of the data set which represent the different variables TR2. Decide on the type of kernel function to use and determine the kernel parameter TR3. Construct the kernel matrix of the NOC data and centre it using Equation (10) TR4. Obtain eigenvalues and their corresponding eigenvectors and rearrange them in a descending order TR5. Orthonormalize the eigenvectors using Equation (11) TR6. Obtain nonlinear components using Equation (13) TR7. Compute monitoring indices (T 2 and Q) based on the kernelized NOC data using Equations (14) and (16) TR8. Determine control limits of T 2 and Q using Equations (20) and (21) Table 2. On-line monitoring.
TT1. Acquire test sample x tt and normalize using the mean and standard deviation values used in step 1 of the off-line stage TT2. Compute the kernel vector of the test sample using Equation (22) TT3. Centre the kernel vector according to Equation (23) TT4. Obtain the PC of the test sample from Equation (25) TT5. Compare the T 2 and Q of the test sample with their respective control limits obtained in the model development stage TT6. If both T 2 and Q are less than their monitoring statistics, the process is in-control. If either T 2 or Q exceeds its control limit, the process is out-of-control and therefore fault identification is carried out to identify the source of the fault Contribution plots which show the contributions of variables to the high statistical index values in a fault region is a common method that is used to identify faults. However, nonlinear PCA-based fault identification is not as straightforward as that of linear PCA due to the nonlinear relationship between the transformed and the original process variables.
In this article, fault variables were identified using a sensitivity analysis principle (Petzold, Barbara, Tech, Li, Cao, & Serban, 2006). The method is based on calculating the rate of change in system output variables resulting from changes in the problem causing parameters (Deng, Tian, & Chen, 2013). Given a test data vector x i ∈ n with n variables, the contribution of the i th variable to a monitoring index is defined by where a i = ∂T 2 /∂x i and b i = ∂Q/∂x i . In this work, the partial derivatives were obtained by differentiating the functions defining T 2 and Q at a given reference fault instant using complex step differentiation (Martins, Sturdza, & Alonso, 2003). This is an efficient generalized approach for obtaining variable contributions in fault identification studies that use multivariate statistical methods.

Tennessee Eastman process
The TE process is a simulation of a real industrial process manifesting both nonlinear and dynamic properties (Downs & Vogel, 1993). It is widely used as a benchmark process for evaluating and comparing process monitoring and control approaches (Chiang et al., 2001). The process consists of five key units: separator, compressor, reactor, stripper and condenser, and eight components coded A to H. The control structure of the process is presented in Figure 2. There are 960 samples and 53 variables, which include 22 continuous variables, 19 composition measurements sampled by 3 composition analysers, and 12 manipulated variables in the TE process. Sampling is done at 3-minute intervals while each fault is introduced at sample number 160. Information on disturbances and baseline operating conditions of the process are documented in (Downs & Vogel, 1993;McAvoy & Ye, 1994).

Application procedure
Five hundred samples obtained under NOC were used as the training data set and all 960 samples obtained under each of the faulty operating conditions were used as test data. All 22 continuous variables and 11 manipulated variables were used in this study. The agitation speed of the reactor's stirrer (the 12th manipulated variable) was not included because it is constant. A total of 20 faults in the process were studied. Descriptions of the variables and faults studied are presented in Tables 3 and 4.
Several methods have been proposed for determining the number of retained PCs. Some of these methods are scree tests, the average eigenvalue approach, crossvalidation, parallel analysis, Akaike information criterion, and the cumulative percent eigenvalue. However, none of these methods have been proved analytically to be the best in all situations (Chiang et al., 2001). In this paper, the number of PCs that explained over 90% of the total variance were retained. Based on this approach, 16 and 17 PCs were selected for PCA and KPCA, respectively.
Another, important parameter for kernel-based methods in model development for process monitoring is the choice of kernel and its width. The radial basis kernel which is a common choice for process monitoring studies Stefatos & Hamza, 2007) was used in this paper. The value of the kernel parameter c was   determined using the relation c = Wnσ 2 , where W is a constant, which is dependent on the process being monitored, n and σ 2 are the dimension and variance of the input space, respectively Mika et al., 1999). The value of W was set at 40 with validation from the training data.
The T 2 and Q statistics were used jointly for fault detection due to their complementary nature. This means that a fault detection was acknowledged when either of the monitoring statistics detected a fault. This is because detectable process variation may not always occur in both the model space and the residual space at the same time.

Fault detection rule
Since measurements obtained from chemical processes are usually noisy, monitoring indices may exceed their thresholds randomly. This amounts to announcing the presence of a fault when no disturbance has actually occurred, that is, a false alarm. In other words, a monitoring index may exceed its threshold once but if no fault is present, the monitoring index may not stay above its threshold in subsequent measurements. Conversely, a fault has likely occurred if the monitoring index stays above its threshold in several consecutive measurements. A fault detection rule is used to address the problem of spurious alarms Tien, Lin, & Jun, 2004;van Sprang, Ramaker, Westerhuis, Gurden, & Smilde, 2002). A detection rule also provides a uniform basis for comparing different monitoring methods. In this paper, successful fault detection was counted when a monitoring index exceeds its control limit in at least two consecutive observations. All algorithms recorded a false alarm rate (FAR) of zero when tested with the training data based on this criterion. Computation of the metrics for evaluating the monitoring performance of the different techniques was therefore based on this criterion.

Computation of monitoring performance metrics
Monitoring performance was based on three metrics: fault detection rates (FDRs), FARs, and detection delay. FDR is the percentage of fault samples identified correctly. It was computed as FDR = n fc n tf × 100, where n fc denotes the number of fault samples identified correctly and n tf is the total number of fault samples. FAR was calculated as the percentage of normal samples identified as faults (or abnormal) during the normal operation of the plant.
where n nf represents the number of normal samples identified as faults and n tn is the total number of normal samples. Detection delay was computed as the time that elapsed before a fault introduced was detected.

Results and discussion
KPCA-based fault detection is demonstrated using Faults 11 and 12 of the TE process. Fault 11 is a random variation in the reactor cooling water inlet temperature, while Fault 12 is a random variation in the condenser cooling   water inlet temperature. The monitoring charts for the two faults are shown in Figures 3 and 4, respectively. The solid curves represent the monitoring indices, while the dash-dot and dash lines represent the control limits at 99% confidence level based on Gaussian distribution and KDE, respectively. It can be seen that in both cases, especially in the T 2 control charts, the KDE-based control limits are below the Gaussian distribution-based control limits. That is, the monitoring indices exceed the KDE-based control limits to a greater extent compared to the Gaussian distribution-based control limits. This implies that using the KDE-based control limits with the KPCA technique gives higher monitoring performance compared to using the Gaussian distribution-based control limits. Table 5 shows the detection rates for PCA, PCA-KDE, KPCA, and KPCA-KDE for all 20 faults studied. The results show that the KDE versions have overall higher FDRs   3 3  3 3  3  2346  1656  1725  1725  4  3  3  3  3  5  3  3  3  3  6  3  3  3  3  7  3  3  3  3  8  4 8  4 8  4 8  4 8  9  2346  1665  1725  1725  10  186  186  180  180  11  15  15  15  15  12  42  30  60  42  13  102  102  111  105  1 4  6  6  6  6  15  ND  1656  1725  1725  16  81  78  81  81  17  45  45  45  45  18  24  24  24  24  19  213  24  213  36  20  105  102  105  105 Note: ND, not detected.
compared to the corresponding Gaussian distributionbased versions. Furthermore, in Table 6, it can be seen that the detection delays of the KDE-based versions are either equal to or lower than the non-KDE-based techniques. This implies that the approaches based on KDEderived UCLs detected faults earlier than their Gaussian distribution-based counterparts. Thus, associating KDEbased control limits with the KPCA technique for fault detection provides better monitoring compared to using control limits based on the Gaussian assumption. KPCA-KDE-based fault identification is demonstrated using Fault 11 as an example. The occurrence of Fault 11 induces change in the reactor cooling water flow rate, which causes the reactor temperature to fluctuate. Both the T 2 -and SPE-based contribution plots at sample 300 shown in Figures 5 and 6, respectively, identified the two fault variables correctly. Variable 9 is the reactor temperature while variable 32 corresponds to the reactor cooling water flow rate. Although it is possible for the control loops to compensate the change in the reactor temperature after a longer time has elapsed, the fluctuations in both variables affected early after the introduction of the fault were correctly identified by the contribution plots.

Test of robustness
To test the robustness of the KPCA-KDE technique, fault detection was performed by varying two parameters: bandwidth and the number of PCs retained. Figure 7 shows the monitoring charts for KPCA and KPCA-KDE with W = 40 for Fault 14. This fault represents sticking of the reactor cooling water valve, which is quit easily detected by most statistical process monitoring approaches. At W = 40, both KPCA and KPCA-KDE  recoded zero false alarms (Figure 7). However, at W = 10, KPCA recorded a false alarm rate of 8.13% while the FAR for KPCA-KDE was still zero (Figure 8). Also, Table 7 shows that the KPCA recorded a similar high FAR when 25 PCs were retained. Conversely, the KPCA-KDE approach still recorded zero false alarms. Thus, apart from generally providing higher FDRs and earlier detections, the KPCA-KDE is more robust than the KPCA technique with control limits based on the Gaussian assumption. A more sensitive technique is better for process operators since less faults will be missed. Secondly, when faults are detected early, operators will have more time to find the root cause of the fault so that remedial actions can be taken before a serious upset occurs. Thirdly, although methods are available for obtaining optimum design parameters for developing process monitoring models, there is no guarantee that the optimum values are used all the time. The reason for this may range from lack of experience of personnel to lack of or limited understanding of the process itself. Therefore, the more robust a technique is, the better it is for process operations.

Conclusion
This paper investigated nonlinear process fault detection and identification using the KPCA-KDE technique. In this approach, the thresholds used for constructing control charts were derived directly from the PDFs of the monitoring indices instead of using thresholds based on the Gaussian distribution. The technique was applied to the benchmark Tennessee Eastman process and its fault detection performance was compared with the KPCA technique based on the Gaussian assumption.
The overall results show that KPCA-KDE detected faults more and earlier than the KPCA with control limits based on the Gaussian distribution. The study also shows that the UCLs based on KDE are more robust than those based on the Gaussian assumption because the former follow the actual distribution of the monitoring statistics more closely. In general, the work corroborates the claim that using KDE-based control limits give better monitoring results in nonlinear processes than using control limits based on the Gaussian assumption. A generalizable approach for computing variable contributions in fault identification studies that centre on multivariate statistical methods was also demonstrated.