Online Monitoring of Multivariate Processes Using Higher-Order

Article pubs.acs.org/IECR

Online Monitoring of Multivariate Processes Using Higher-Order Cumulants Analysis Youqing Wang,*,† Jicong Fan,† and Yuan Yao‡ †

College of Information Science and Technology, Beijing University of Chemical Technology, Beijing 100029, China Department of Chemical Engineering, National Tsing Hua University, Hsinchu 30013, Taiwan

‡

ABSTRACT: In this study, a novel approach is developed for online state monitoring based on higher-order cumulants analysis (HCA). This approach applies higher-order cumulants to monitor a multivariate process, and while conventional approaches such as independent components analysis (ICA) uses variance to monitor process. Variance is lower-order statistics and is only sensitive to amplitude. In contrast, higher-order cumulants, the typical higher-order statistics, carry important information and are sensitive to both amplitude and phase, particularly for non-Gaussian distributions. The main idea of this novel approach is to monitor the cumulants of dominant independent components and residuals of the ICA model. Therefore, higher-order statistical information of multivariate processes can be monitored online. Furthermore, a variable contribution analysis scheme is developed for HCA to diagnose faults. The proposed approach is applied to the Tennessee Eastman (TE) process to exhibit its effectiveness. The results demonstrate that HCA outperforms ICA and dynamic ICA significantly and the variable contribution analysis of HCA can diagnose faults successfully. to extract ICs.16,17 Lee et al.18 first applied ICA to a nonGaussian process and verified its superiority compared with PCA. In the ICA approach, I2 and SPE (sometimes including Ie2) are used for detecting faults and variable contribution is applied to isolate faults. However, in ICA, repetitive computation may lead to different solutions because of the random initiation. Moreover, the ICs cannot be sorted according to variance because their variances are all unit. Lee et al.19 proposed a modified ICA that could determine consistent solutions and sort ICs according to the variances of the PCA scores. Wang et al.20 developed a method to sort ICs and compared it with conventional methods. A number of extensional methods of ICA have also been developed to enhance the monitoring performance.21−32 Particularly, dynamic ICA (DICA)21,26 was proposed to deal with process dynamics by including several previous sampled data to augment the measured vector. When indices of PCA and ICA detect a fault, a variable contribution analysis is usually implemented for fault diagnosis. The main idea of the contribution analysis is that the variables giving the largest contributions to the index that detects a fault are the fault points. Alcala and Qin33 summarized and compared the existing different methods of variable contribution analysis and suggested a relative contribution, which can set a level ground to compare the contributions of different variables. In the step of dimension reduction and feature extraction, ICA uses negentropy or mutual information (higher-order statistics) to extract independent latent variables; however, in the step of online monitoring, ICA monitors the mean values

1. INTRODUCTION In order to ensure the safety of industrial processes and improve the quality of products, online state monitoring is becoming increasingly important.1 In multivariable and complex processes, multivariable statistical analysis (MSA) is widely used because MSA not only takes advantage of the large number of process measurements but also avoids considering the complex systematic models2−5 that are difficult to identify. Principal component analysis (PCA)6 is a well-known typical MSA technique. Using orthogonal projection and maximizing variance, PCA extracts relatively fewer uncorrelated latent variables to replace the measured variables. In general, PCA can deal with high-dimensional, high-correlative, and noisy data. In the PCA approach,7,8 two indices, Hotelling-T2 (T2) and squared prediction error (SPE), are applied for fault detection, where T2 and SPE are applied to monitor the selected principal components and the residuals of the PCA model, respectively. A number of extensional methods of PCA have been developed to enhance the monitoring abilities for processes with different characters.9−14 As such, PCA is a second-order method that only exploits the mean value and variance to give uncorrelated components. The mean value and variance are sufficient to describe a Gaussian distribution but not sufficient to describe a non-Gaussian distribution. Therefore, the independent component analysis (ICA)15−17 method was developed. Compared with PCA, ICA is a higher-order method that exploits higher-order statistics to determine independent components (ICs). As mentioned above, PCA extracts uncorrelated components, while ICA extracts independent components. Because the term “uncorrelated” means only partly independent,17 hence ICA can be considered an enhanced version of PCA. Because independence can be derived by minimizing the mutual information, by maximizing non-Gaussianity, and so on, higher-order statistics such as entropy and kurtosis are utilized © 2014 American Chemical Society

Received: Revised: Accepted: Published: 4328

June 11, 2013 February 22, 2014 February 23, 2014 February 24, 2014 dx.doi.org/10.1021/ie401834e | Ind. Eng. Chem. Res. 2014, 53, 4328−4338

Industrial & Engineering Chemistry Research

Article

2. PRELIMINARY: INDEPENDENT COMPONENT ANALYSIS Given that m measured variables x = [x1,x2,...,xm]T are linear combinations of d independent latent variables s = [s1,s2,...,sd]T and d ≤ m, the following relationship can be obtained: (1) x = As where A denotes the mixed matrix. The objective of ICA is to find a demixed matrix W to obtain the source variables: (2) s = Wx It is difficult to calculate W and s directly. PCA is usually utilized to whiten x in order to obtain the uncorrelated variables as follows: z = Qx (3)

and variances (lower-order statistics) of the latent variables and the modeling residuals, and hence it cannot monitor the higherorder statistics. In short, ICA is a higher-order method in terms of extracting latent variables, but it is a lower-order method in terms of monitoring. Therefore, both PCA and ICA cannot take advantage of the higher-order statistics (HOS)34 to monitor higher-order statistical information of processes. In the past several years, HOS have been widely used for solving the problems of non-Gaussianity, nonminimum phase, colored noise, and nonlinearities. Compared with first- and secondorder statistics that include mean variance, autocorrelation, and power spectrum, HOS include higher-order moments, higherorder cumulants, and their respective polyspectras. It should be noticed that all moments and their polyspectras are only sensitive to amplitude but unable to describe phase information.34,35 Here, the concept “phase” is borrowed from Fourier transform and also can be understood as the continuous changes of a time series within several sampling times. In contrast, higher-order cumulants and their polyspectras are more widely used because they are sensitive to not only amplitude but also phase. In addition, higher-order cumulants with different time lags are exactly higher-order correlations, while second-order cumulants are autocorrelation and cross-correlation. PCA and ICA monitor only the variances of selected features and model residuals; their indices only include the information of the current samples, and the information of past several samples is neglected. Hence, PCA and ICA neither can monitor higher-order statistics nor are sensitive to the changes in the phase. Furthermore, dynamic PCA and dynamic ICA consider low-order correlations rather than high-order correlation. On the other hand, HOS, particularly higher-order cumulants, have been widely applied on signal detection and classification, but these applications mainly focus on one-dimensional (1-D) signals.36−38 Therefore, in this study, we developed a novel approach for fault detection and isolation based on a higher-order cumulants analysis (HCA). In this approach, first, ICA is utilized to extract non-Gaussian features; second, the higher-order cumulants of the dominant ICs and the residuals of the ICA model are monitored; finally, a contribution analysis is developed for HCA to isolate faults. The Tennessee Eastman (TE)39−41 process is used for exhibiting the superiority of HCA compared with ICA and DICA. The contributions of this study can be summarized as follows: First, to the best knowledge of the authors, this is the first time that HOS are used for the state monitoring of multivariate processes. Second, a novel combination of ICA and higher-order cumulants analysis, abbreviated as HCA, is proposed to deal with this issue. Finally, a new variable contribution analysis method is proposed for HCA and hence can be used for fault isolation. Conventional methods (PCA, DPCA, ICA, and DICA) use low-order statistical information of a multivariate process for monitoring, but HCA uses the higher-order statistical information for monitoring and further to take advantage of that to diagnose faults. In addition, compared with DPCA and DICA that consider the previous sampled data to exploit second-order time dependence, HCA not only considers the previous sampled data but exploits higher-order autocorrelation (cumulants). Theoretically, HCA outperforms PCA, DPCA, ICA, and DICA, because HCA is considerably more sensitive to the variations of both amplitude and phase influenced by faults, while the other methods are only sensitive to amplitude variations.

where z denotes the whitened scores vector and Q denotes the whitening matrix, which is calculated by Q = Λ−1/2P

(4)

and x can be represented by x = PTΛ1/2z

(5)

P is calculated by using eigenvalue decomposition as follows: 9 = PTΛ−1P

(6) T

T

where 9 = E(xx ) = (1/n)XX denotes the covariance matrix of x; X (m × n), a matrix containing n samples of x; Λ, a diagonal matrix of eigenvalues; and P, the corresponding eigenvector matrix. It will be easier to determine independent variables from the uncorrelated variables than from the measured variables x. Assuming that (7)

z = Bs

where B (m × d) denotes an orthogonal matrix, we obtain the following: s = BTz

(8)

The independence of s can be evaluated by mutual information, non-Gaussianity, and so on.17 In terms of non-Gaussianity, according to the central-limit theorem, a larger non-Gaussianity implies more independence. Hence, s and B can be derived by maximizing the non-Gaussianity of BTz. In this study, the modified ICA algorithm proposed by Lee et al.19,22 is used for calculating s and B. The modified ICA algorithm is actually a modified form of the FastICA algorithm developed by Hyvärinen.16 After B is obtained, by combining eqs 2, 3, and 8, we obtain the following demixed matrix: W = BTQ

(9)

.By combining eqs 5 and 7, we obtain the following mixed matrix: A = PTΛ1/2B

(10)

In order to reduce the number of ICs, d dominant ICs can be selected by applying the average value criterion42 to the norms of the columns of A.20 The d columns of A with the largest norms form Ad, the selected ICs form sd, and the corresponding rows of W form Wd. Hence, the total projection of d dominant ICs to x is as follows:

x ̂ = Adsd = AdWdx 4329

(11)

dx.doi.org/10.1021/ie401834e | Ind. Eng. Chem. Res. 2014, 53, 4328−4338


Article

However, both categories of the diagnosis method only utilize lower-order statistics. Therefore, using higher-order statistics to diagnose faults should be explored. 3.2. Higher-order Monitoring Statistics. As such, higherorder (order > 2) statistics include higher-order moments, higher-order cumulants, and their respective polyspectras. Higher-order moments and their polyspectras are only sensitive to amplitude, but higher-order cumulants and their polyspectras are sensitive to both amplitude and phase. This leads to the following question: why the kurtosis or entropy is not monitored when the ICA algorithm uses kurtosis or entropy to extract ICs? The answer for this question is that kurtosis is based on a fourth-order moment but not fourth-order cumulant and the online estimation of entropy is difficult. Therefore, higher-order cumulants and their polyspectras are better choices to monitor higher-order statistical information. In this research, only the characteristics in the time domain are studied for the sake of simplicity and convenience. Hence, the polyspectras will not be discussed. In the physical sense, cumulants indicate correlation or dependency; i.e., a secondorder cumulant is a second correlation or dependency and higher-order cumulants are higher-order correlations or dependency.36 Cumulants present the cumulative effect of continuous variation if the time lags are different. In general, higher-order cumulants are widely used for solving problems where either a non-Gaussianity or a nonminimum phase exists.34 For example, higher-order cumulants are very effective in non-Gaussian signal detection and classification.36−38 Since the estimation of a cumulant with a higher order is not accurate, the third- and fourth-order cumulants are the two most popular cumulants, which for a zero-mean variable y(i) are defined as follows:34−37

and the projection residual (also called prediction error) is given by e = x − x ̂ = (I − AdWd)x = Lx

(12)

where I denotes an identify matrix and L = I − AdWd. The first index of ICA is d

I 2 = sdTsd =

∑ sp2 p=1

(13)

which is used for monitoring the magnitude variations of all the dominant ICs. Another index of ICA is m

SPE = e Te =

∑ eq2 q=1

(14)

which is used for monitoring the magnitude variations of the ICA model residuals for all the measured variables. Because the distributions of I2 and SPE are unknown, their control limits should be calculated by using kernel density estimation18,43 under a restraint (0,+∞) on the normal process data. Equation 13 denotes the quadratic sum of d components, and it can monitor the variances of components, of which the mean values are zeros in a normal process. Equation 14 denotes the quadratic sum of m modeling residuals, and it can monitor the variances of the modeling residuals. These two indices calculate the sum of deviations of all component/residuals from their expectations (zeros). Hence, one can say that theses two indices use sampled variances to monitor their expectations.

3. HIGHER-ORDER CUMULANTS ANALYSIS 3.1. Motivations. In the ICA approach, according to eqs 13 and 14, I2 and SPE are just the quadratic sum of dominant components and the squared sum of modeling residuals, respectively, and therefore they are exactly to monitor mean values and variances. In other words, ICA monitors the loworder statistical information (mean, first-order moment, and variance, second-order moment) but ignores higher-order statistics. ICA extracts the most systematic non-Gaussian information by retaining several ICs. However, the monitoring statistic I2 cannot appropriately summarize the non-Gaussian information in the retained ICs. Further, the residuals contain some non-Gaussian information (and some Gaussian information if possible). Dynamic ICA also has these above-mentioned defects. As mentioned in the Introduction, for non-Gaussian signals, higher-order statistics include very important information that cannot be described by lower-order statistics. If a fault makes a Gaussian variable be non-Gaussian, or makes a nonGaussian variable get different characteristics, the variations of higher-order statistics are considerably more evident than those of lower-order statistics. Hence, it is necessary to monitor the higher-order statistics of the retained ICs and model residuals. Moreover, the monitoring indices of ICA only reflect the information at the current time but cannot reflect the cumulative effect of the information at several previous time points. Compared with ICA, DICA considers the previously sampled data and utilizes second-order time dependence to monitor dynamic processes but does not exploit higher-order time dependence. Therefore, the cumulative effect and higherorder time dependence should be considered if possible. Finally, after one index of ICA or DICA detects a fault, a variable contribution plot or classification/discrimination based on the selected features is usually applied for fault diagnosis.

cum3, y(τ1 , τ2) = E{y(i) y(i + τ1) y(i + τ2)}

(15)

cum4, y(τ1 , τ2 , τ3) = E{y(i) y(i + τ1) y(i + τ2) y(i + τ3)} − E{y(i) y(i + τ1)} E{y(i + τ2) y(i + τ3)} − E{y(i) y(i + τ2)} E{y(i + τ1) y(i + τ3)} − E{y(i) y(i + τ3)} E{y(i + τ1) y(i + τ2)}

(16)

where τ1, τ2, and τ3 denote the time lags. In general, third- and fourth-order cumulants are sufficient to solve a non-Gaussian problem and widely used to present higher-order statistical information. In this work, only a third-order cumulant is utilized for illustration. The sampled estimation of the thirdorder cumulant is defined as follows: cum3, y(τ1 , τ2) =

1 N

N

∑ y(i) y(i + τ1) y(i + τ2) i=1

(17)

3.3. Determining Time Lags of Cumulants. It is important to determine the time lags of cumulants. However, there is no criterion to determine the time lags of cumulants for the exploration of a new sequence or process. One of the properties of cumulants is that the cumulants of independent, identically distributed random sequence data are δ functions;34 i.e., only for τ1 = τ2 = ... = τk−1 = 0, cumk,y(τ1,τ2,...,τk−1) ≠ 0. In fact, this property is a special case in which cumk,y(τ1,τ2,...,τk−1) = 0 if a subset of k random variables {yi} is independent of the rest. For a random process, a larger sample interval indicates 4330



Article

getting more independent samples. Hence, if one of τ1,τ2,...,τk−1 is sufficiently large, cumk,y(τ1,τ2,...,τk−1) = 0, which is not a wise choice. Meanwhile, this is the basis that cumulants can be used for identifying the ARMA model. If |τi − τj| ≥ p (if there exists a pair of i and j, p is of the AR order), cumk,y(τ1,τ2,...,τk−1) = 0. Therefore, the first principal for determining the time lags of cumulants is that |τi − τj| < p. A considerably useful parameter setting of a kth-order cumulant in signal processing is the 1-D slice, τ2 = τ3 = ... = τk−1.34 However, for a third-order cumulant, a 1-D slice makes no sense because only τ1 and τ2 exist in the third-order cumulant. In the work of Giannakis and Tsatsanis,36 a thirdorder cumulant is used for detecting and classifying a signal and the time lags are set to be τ1 = τ2 = 0. Hence, cum3,y(0,0) is just the third-order center moment. In this study, it is suggested that τ1 and τ2 have different values. Moreover, in order to include more information of a past time, it is proposed that both τ1 and τ2 are not zeros. In general, because |τi − τj| < p and the AR order of most processes are p = 1 or 2,12,21 in this study, the time lags of the third-order cumulant are finally determined as τ1 = 1 and τ2 = 2. Moreover, for the kth-order cumulant cumk,y(τ1,τ2,...,τk−1), τi = i, and k = p, where p is of the AR order and can be determined by identification or the self-acting method proposed in ref 8. The physical meaning of cum3,y(1,2) is a cumulative effect of three continuous samples, also known as the third-order correlation (time dependency) containing phase information. It should be noted that if a process is not dynamic or the samples are completely independent, i.e., there is no time dependence, the monitoring of cumulants is unnecessary. But, in that case, the cumulants are zero for variables having a zero mean, and hence, the monitoring of cumulants is still useful and will identify to the monitoring of mean values, because the cumulants will not be zero if a fault makes a variable’s mean value be nonzero (as the first paragraph of this section tells that cumk,y(τ1,τ2,...,τk−1) = 0 if there is no time dependence). 3.4. Monitoring Cumulants of Multivariate Process. Conventionally, cumulants are used for detecting and classifying a single signal when a number of samples have been obtained and the cumulants are estimated;36−38 which is inclined to be an off-line method. In multivariate statistical process monitoring (MSPM),2,3,7 sampled variances, shown by eqs 13 and 14, are used to monitor low-order statistical information, which is an online method. Analogously, in order to develop an online and real-time method, the sampled cumulants are used in this study to monitor higher-order statistical information. In terms of multivariate processes or multisignals, it is not advisable to monitor the cumulant of every measured variable separately because the measured variables may not be independent from each other and the high dimension is also a problem. Therefore, ICA is proposed to extract independent components (ICs) and reduce the dimension; then, the cumulants of the selected components and model residuals can be monitored. As mentioned in Preliminary: Independent Component Analysis, the number of ICs is reduced by applying the average value criterion to the norms of the columns of the mixing matrix A. This leads to the following question: Why are kurtosis or cumulants not used for reducing the number of ICs when HOS are monitored? The answer to this question is that the retaining ICs should contain the major variation or energy of a process, which cannot be ensured by a kurtosis- or cumulant-based dimension reduction. ICA is effective in

extracting non-Gaussian features and works well when the source signals contain only one or no Gaussian signal.17 If a Gaussian IC is extracted and its cumulant is monitored, the monitoring of cumulants for the multivariate process is still significant because a fault may make the Gaussian IC nonGaussian. There is some Gaussian information in the residuals of the ICA model. The monitoring for the cumulants of Gaussian information will identify with the monitoring for variances or mean values because the calculated cumulants will change if the mean values are no longer zero. After ICA extracts ICs, the sampled third-order cumulant of the pth ICs at the ith point is defined as follows: hsp(i) = sp(i) sp(i − 1) sp(i − 2) = wpx(i) wpx(i − 1) wpx(i − 2)

(18)

where wp denotes the pth row of Wd and p = 1, 2, ...,d. The physical meaning of hsp(i) is the cumulative effect of three continuous samples of the pth component. Since the ranking of time lags does not influence a cumulant34 (e.g., cum3,y(τ1,τ2) = cum3,y(τ2,τ1)), eq 18 will lead to the same result as that obtained by defining hsp(i+2) = cum3,sp,sampled(1,2) = sp(i) sp(i +1) sp(i+2) and replacing i by i + 2. In order to monitor the third-order cumulants of all of the selected ICs, the first index of HCA is defined as follows: d

HS(i) =

∑ p=1

hsp(i) − mhsp vhsp

(19)

where mhsp and vhsp denote the mean value and the standard deviation of hsp, respectively. HS is in fact the total deviation of the sampled third-order cumulants of all of the retained ICs from their third-order cumulants, and it is to monitor the variation of the cumulants of whole dominant components. This leads to the following question: Why is the sum of the absolute cumulants used and not the quadratic sum or the direct sum of the sampled cumulants? The quadratic sum used in PCA or ICA is used for monitoring the variances. The absolute sum can be used for monitoring the standard deviation in PCA or ICA. However, the quadratic sum will be more sensitive than the absolute sum, and the distribution of the quadratic sum for the Gaussian distribution is well-known. Giannakis and Tsatsanis36 utilized |cum3,y| to detect a nonGaussian signal because cum3,y = 0 for a Gaussian signal. If the quadratic sum of the cumulants is used, the variances of the sampled cumulants are inclined to be monitored. Since a thirdorder cumulant is more sensitive to the amplitude than the variance, the absolute sum is sufficient to monitor the cumulants of the multivariate process. On the other hand, if the quadratic sum is used, a fault diagnosis based on the variable contribution analysis will be very difficult to implement. Hence, in this work, the absolute sum is applied for the multivariate cumulant monitoring. As shown in eq 19, the statistical meaning of HS is the weighted sum of the absolute distance of all the online cumulants from their expectations, where the weighting factors are the reciprocals of the standard deviations of the corresponding cumulants. Analogously, the sampled third-order cumulant of the prediction error for the qth variable at the ith point is defined as follows: 4331



Article d

heq(i) = eq(i) eq(i − 1) eq(i − 2) = lqx(i) lqx(i − 1) lqx(i − 2)

p=1

(20)

× [wp , jxj(i) sp(i − 1) sp(i − 2) − cmhsp , j]/vhsp}

where lq denotes the qth row of L and q = 1, 2, ..., m. The physical meaning of heq(i) is the cumulative effect of three continuous samples of the qth model residual. In order to monitor the third-order cumulants of all of the prediction errors, the second index of HCA is defined as follows: m

(24)

and cmhsp,j denotes the contribution of the jth variable to the mean value of hsp. By using eq 18, we obtain the following: cmhsp , j = E{wp , jxj(i) sp(i − 1) sp(i − 2)}

heq(i) − mheq

∑

HE(i) =

∑ {sign[sp(i) sp(i − 1) sp(i − 2) − mhsp]

CjHS(i) =

vheq

q=1

(21)

where mheq and vheq denote the mean value and the standard deviation of heq, respectively. HE is in fact the total deviation of the sampled third-order cumulants of the ICA model residuals from their third-order cumulants. Since the distributions of HS and HE are unknown, their control limits can be calculated by using the kernel density estimation under a constraint (0,+∞) based on the history normal data. In this study, a confidence limit of 99% is assumed to be the control limit. 3.5. Diagnosing Using Variable Contribution Analysis for HCA. Being different from I2 and SPE of ICA, which only include the current information of the dominant ICs and the model residuals, HS and HE of HCA use the cumulant information and, hence, are considerably more sensitive to the variations of amplitude and/or phase. After either HS or HE has detected a fault, the contribution analysis method33 can be used for isolating the fault. The variable contribution to a monitoring index can be defined as follows:

∑ wp,jxj(i) sp(i − 1) sp(i − 2) i=3

(25)

m

∑ |[eq(i) eq(i − 1) eq(i − 2) − mheq]/vheq|

HE(i) =

q=1 m

∑ |[lqx(i) eq(i − 1) eq(i − 2) − mheq]/vheq|

=

q=1 m

∑ {sign[eq(i) eq(i − 1) eq(i − 2) − mheq]

=

q=1

× [lqx(i) eq(i − 1) eq(i − 2) − mheq]/vheq} m


=

q=1

×

m

∑ [lq , jxj(i) eq(i − 1) eq(i − 2) − cmheq , j]/vheq} j=1

m

∑ Cjindex

m

∑ ∑ {sign[eq(i) eq(i − 1) eq(i − 2) − mheq]

= (22)

j=1

n

By combining eqs 20 and 21, we can represent HE as follows:

m

index =

1 n−2

≐

j=1 q=1

× [lq , jxj(i) eq(i − 1) eq(i − 2) − cmheq , j]/vheq}

Cindex j

where denotes the contribution of the jth variable to index. In HCA, by combining eqs 2, 18, and 19, we can represent HS as follows:

m

∑ CjHE(i)

=

j=1 d

HS(i) =

(26)

∑ |[sp(i) sp(i − 1) sp(i − 2) − mhsp]/vhsp|

where

p=1

m

d

=

CjHE(i) =

∑ |[wpx(i) sp(i − 1) sp(i − 2) − mhsp]/vhsp|

q=1

p=1 d

=

× [lq , jxj(i) eq(i − 1) eq(i − 2) − cmheq , j]/vheq}

∑ {sign[sp(i) sp(i − 1) sp(i − 2) − mhsp] d

cmheq , j = E{lq , jxj(i) eq(i − 1) eq(i − 2)}

∑ {sign[sp(i) sp(i − 1) sp(i − 2) − mhsp] p=1 m

×

≐

∑ [wp,jxj(i) sp(i − 1) sp(i − 2) − cmhsp,j]/vhsp} j=1

m

=

1 n−2

n

∑ lq,jxj(i) eq(i − 1) eq(i − 2) i=3

(28)

In order to obtain a level ground for comparing the contribution of different variables, every variable’s contribution should be normalized by its control limit:

d

∑ ∑ {sign[sp(i) sp(i − 1) sp(i − 2) − mhsp] j=1 p=1

× [wp , jxj(i) sp(i − 1) sp(i − 2) − cmhsp , j]/vhsp}

CjHS = CjHS/limit_CjHS

m

=

(27)

and cmheq,j denotes the contribution of the jth variable to the mean value of heq. By using eq 20, we obtain the following:

p=1

× [wpx(i) sp(i − 1) sp(i − 2) − mhsp]/vhsp} =


∑ CjHS(i)

CjHE = CjHE/limit_CjHE

j=1

(23)

(29)

HE where limit_CHS can be calculated by using the j and limit_Cj kernel density estimation and normal historical data. In this

where 4332



Article

Figure 1. Flowchart of HCA.

contribution analysis, the variables with the largest contributions are assumed to be the fault points. Concerning the above-mentioned contribution analysis (three-order cumulant based HCA), it is worth noticing that the contribution analysis of fourth-order cumulant based HCA will be much more complicated, not to mention that of the higher-orders (≥5) cumulant based HCA. Moreover, the

estimation of cumulants with much higher orders (≥5) are not accurate. Therefore, three-order cumulant is considerably appropriate for both fault detection and fault diagnosis, and fourth-order cumulant will be preferred if fault diagnosis is not needed. The flowchart of HCA is shown in Figure 1. In both offline modeling and online monitoring, PCA prepares for ICA with 4333



Article

Figure 2. Structure chart of the TE process20

Table 1. Monitored Variables in the TE Process22

uncorrelated scores and ICA prepares for HCA with independent components. By using the historical normal data, we obtain the HCA model including the projection model and the control limits. The online measurements are substituted into the HCA model for evaluating the degree of deviation. If the degree is beyond the confidence limit, we deduce that a fault has occurred and then the variable contribution analysis is implemented to diagnose the fault.

(1) (2) (3) (4)

A feed (stream 1) D feed (stream 2) E feed (stream 3) total feed (stream 4)

(5) recycle flow (stream 8) (6) reactor feed rate (stream 6) (7) reactor pressure (8) reactor level (9) reactor temperature (10) purge rate (stream 9) (11) product separator temperature (12) product separator level

4. CASE STUDY The Tennessee Eastman44,45 is a widely used benchmark to test and evaluate the performances of process control and monitoring strategies.3,8,21,22,26 Its structure chart is shown in Figure 2.20 This process includes five major parts, namely, the reactor, condenser, compressor, separator, and the stripper, through which four feeds A/C/D/E and one inert B generate two products G/H and one byproduct F. There are 12 manipulated variables and 41 measured variables (22 continuous process measurements and 19 composition measurements) in the entire process. Since the 19 composition measurements are difficult to measure in real time and one manipulated variable, the agitation speed, is not manipulated, in this study, 33 variables are monitored in all, which are listed in Table 1. The TE process has 21 predefined faults, which are described in Table 2. A history normal data set including 480 samples is used for training the monitoring models of ICA and HCA. Because it has been proved that ICA outperforms PCA in many studies, in this study, HCA is only compared with ICA and DICA. The testing data set of every fault includes 960 samples, and every fault is introduced from sample 160. In both the monitoring models of ICA and HCA, eight dominant ICs are selected by applying the average value criterion42 to the norms of the columns of the mixed matrix A.20 In the DICA model, the time lag is 1 and two different numbers (8 and 16 respectively) of dominant ICs are selected to exhibit the performance of DICA. All of these methods’ control limits are based on the confidence limit of 99% that is

(13) product separator pressure (14) product separator underflow (stream 10) (15) stripper level (16) stripper pressure (17) stripper underflow (stream 11)

(18) stripper temperature (19) stripper steam flow (20) compressor work (21) reactor cooling water outlet temperature (22) separator cooling water outlet temperature (23) D feed flow valve (stream 2) (24) E feed flow valve (stream 3) (25) A feed flow valve (stream 1) (26) total feed flow valve (stream 4) (27) compressor recycle valve (28) purge valve (stream 9) (29) separator pot liquid flow valve (stream 10) (30) stripper liquid product flow valve (stream 11) (31) stripper steam valve (32) reactor cooling water flow (33) condenser cooling water flow

calculated by the kernel density estimation.22,43 For a normal testing data set, the false alarm rates are listed in Table 3, which shows that all of the false alarm rates are acceptable. Figure 3 shows a comparison of the fault detection rates among the first indices of four methods. The detection rates of every first index for faults 1, 5, 6, 7, 12, and 14 are all close to 100% because the magnitudes of these faults are very large. Except for these faults, HCA and DICA(16 ICs) have considerably higher detection rates than ICA and DICA(8 ICs), particularly for faults 4, 10, 11, 16, 19, 20, and 21. For faults 17, 19, and 20, DICA(16 ICs) gives slightly higher detection rates than HCA, but for faults 3, 4, 9, 11, 15, and 21, HCA outperforms DICA(16 ICs) markedly. Figure 4 shows the performance of the second index of each method. It can be 4334



Article

Table 2. Description of Faults in the TE Process22 no. 1

description

type

8

A/C feed ratio, B composition constant (stream 4) B composition, A/C ratio constant (stream 4) D feed temperature (stream 2) reactor cooling water inlet temperature condenser cooling water inlet temperature A feed loss (stream 1) C header pressure loss−reduced availability (stream 4) A, B, and C feed composition (stream 4)

9

D feed temperature (stream 2)

10

C feed temperature (stream 4)

11

reactor cooling water inlet temperature

12

condenser cooling water inlet temperature

13 14 15 16−20 21

reaction kinetics reactor cooling water valve condenser cooling water valve unknown valve for stream 4 was fixed at steady-state position

2 3 4 5 6 7

step step step step step step step random variation random variation random variation random variation random variation slow drift sticking sticking unknown constant position

Figure 4. Detection rate of the second index of each method.

HS chart of HCA, there is almost no point under the control limit. In the HE of HCA, after fault 4 occurs, there is no point under the control limit, but in the SPE charts of ICA, DICA(8 ICs), and DICA(16 ICs), a few points are regrettably under the control limit. In Figure 6, similarly to Figure 5, HCA performs better than the other three methods for fault 19. In general, these results demonstrate that HCA is more effective than other methods, DICA(16 ICs) is better than DICA(8 ICs), and ICA should be the worst. What’s more, HCA has fewer dominant ICs and a simpler model compared with DICA(16 ICs). The following facts might be reasons for the superior performance of HCA: first, HCA takes advantage of higher-order statistical information that is more sensitive to faults than low-order statistics; second, HCA uses higher-order cumulants to obtain higher-order autocorrelation, while DICA just obtains second-order time dependence by augmenting the measured vector with previous sampled data. Table 4 shows the variances and the absolute values of the cumulants of eight dominant ICs in the normal process and fault 4. It is found that fault 4 has considerably larger influences on the third-order cumulants than on the variances of the dominant ICs. In order to demonstrate the effectiveness of the proposed contribution analysis, faults 4, 6, and 12 are considered. In Figure 7 (the contribution plot of HCA for fault 4 at sample 161), variable 9 (reactor temperature) has the largest contribution to both HS and HE, which implies that the reactor cooling water inlet temperature has a fault; i.e., fault 4 occurs. In Figure 8 (the contribution plot of HCA for fault 6 at sample 161), variable 1 (A feed) and variable 25 (A feed flow valve) have the largest contribution to both HS and HE, which indicates that feed A has a fault; i.e., fault 6 occurs. In Figure 9 (the contribution plot of HCA for fault 12 at sample 164 when the indices are completely separated from their control limits), variables 11 (product separator temperature), 22 (separator cooling water outlet temperature), and 33 (condenser cooling water flow) have the largest contribution to both HS and HE, which indicate that the separator temperature or the condenser cooling water has a fault. In fact, as shown by the TE structure chart in Figure 2, the condenser can affect the separator directly and the function is irreversible. Hence, we deduce that the condenser cooling water has a fault, which further influences the separator temperature. All of the above have demonstrated

Table 3. False Alarm Rates of ICA and HCA method index false alarm rate (%)

ICA(8 ICs)

DICA(8 ICs)

DICA(16 ICs)

HCA(8 ICs)

I2

SPE

I2

SPE

I2

SPE

HS

HE

0.4

1.5

0.2

1.9

0.4

1.8

0.4

1.6

Figure 3. Detection rate of the fist index of each method.

found that, for every fault, HCA has the highest detection rate and the superiority is much more remarkable than the other methods especially for faults 3, 9, 11, 15, 16, 19, 20, and 21. To give more detailed comparisons, the monitoring charts for faults 4 and 19 are shown in Figures 5 and 6, respectively. In Figure 5, for the I2 charts of ICA, DICA(8 ICs), and DICA(16 ICs), there are approximately 45%, 90%, and 30%, respectively, of the points that are under the control limit after the reactor cooling water inlet temperature has a step change, but in the 4335



Article

Figure 5. Monitoring charts of fault 4.

Figure 6. Monitoring charts of fault 19.

Table 4. Variance and Third-Order Cumulant of Eight Dominant ICs in Normal Situation and When There Is Fault 4 s1

s2

s3

normal fault

1 1.12

1 1.24

1 1.13

normal fault

0.07 31.17

0.39 2.41

0.015 0.051

s4 Variance 1 1.39 |cum3(1,2)| 0.04 0.35

s5

s6

s7

s8

1 1.03

1 1.11

1 1.28

1 1.22

0.08 7.07

0.04 4.47

0.04 2.76

0.12 14.94

higher-order statistics (information), which are more sensitive to faults. The case study on the TE process demonstrates that HCA is considerably more effective than ICA and DICA in fault detection and the proposed contribution analysis for HCA could diagnose faults exactly. Nevertheless, more studies are required to further improve HCA. For example, the multivariate cumulants monitoring index, HS, does not consider the different importance of the selected ICs and it might be better if

that the proposed contribution analysis is effective in diagnosing faults.

5. CONCLUSION A novel approach, HCA, was developed for online state monitoring of multivariate processes, based on higher-order cumulants and ICA. Compared with ICA that uses secondorder statistics to monitor processes, HCA takes advantage of 4336



Article

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (Grant 61374099), the Beijing Nova Program (Grant 2011025), the Program for New Century Excellent Talents in University (Grant NCET-13-0652), and the National Science Council of Taiwan under Grant No. NSC 102-2221-E-007-130.

■ Figure 7. Variable contribution plot for fault 4 at sample 161.

Figure 8. Variable contribution plot for fault 6 at sample 161.

Figure 9. Variable contribution plot for fault 12 at sample 164.

a weight factor is added to every cumulant. Moreover, it is still an open problem whether the polyspectras of higher-order cumulants can be well-applied to multivariate fault detection and diagnosis.

■

REFERENCES

(1) Venkatasubramanian, V.; Rengaswamy, R.; Kavuri, S. N.; Yin, K. A review of process fault detection and diagnosis Part III: Process history based methods. Comput. Chem. Eng. 2003, 27, 327−346. (2) Qin, S. J. Statistical process monitoring: basics and beyond. J. Chemom. 2003, 17, 480−502. (3) Zhou, D. H.; Li, G.; Qin, S. J. Total projection to latent structures for process monitoring. AIChE J. 2011, 56, 168−178. (4) Yu, J.; Qin, S. J. Multiway Gaussian mixture model based multiphase batch process monitoring. Ind. Eng. Chem. Res. 2009, 48, 8585−8594. (5) Li, G.; Alcala, C. F.; Qin, S. J.; Zhou, D. H. Generalized reconstruction-based contributions for output-relevant fault diagnosis with application to the Tennessee Eastman process. IEEE Trans. Control Syst. Technol. 2011, 19, 1114−1127. (6) Wold, S.; Esbensen, K.; Geladi, P. Principal component analysis. Chemom. Intell. Lab. Syst. 1987, 2, 37−52. (7) Nomikos, P.; MacGregor, J. F. Monitoring batch processes using multiway principal component analysis. AIChE J. 1994, 40, 1361− 1375. (8) Ku, W.; Storer, R. H.; Georgakis, C. Disturbance detection and isolation by dynamic principal component analysis. Chemom. Intell. Lab. Syst. 1995, 30, 179−196. (9) Chen, J.; Liao, C.-M. Dynamic process fault monitoring based on neural network and PCA. J. Process Control 2002, 12, 277−289. (10) Lee, J. M.; Yoo, C. K.; Choi, S. W.; Vanrolleghem, P. A.; Lee, I. B. Nonlinear process monitoring using kernel principal component analysis. Chem. Eng. Sci. 2004, 59, 223−234. (11) Alcala, C. F.; Qin, S. J. Reconstruction-based contribution for process monitoring with kernel principal component analysis. Ind. Eng. Chem. Res. 2010, 49, 7849−7857. (12) Cheng, C. Y.; Hsu, C. C.; Chen, M. C. Adaptive kernel principal component analysis (KPCA) for monitoring small disturbances of nonlinear processes. Ind. Eng. Chem. Res. 2010, 49, 2254−2262. (13) Shams, M. A. B.; Budman, H. M.; Duever, T. A. Fault detection, identification and diagnosis using CUSUM based PCA. Chem. Eng. Sci. 2011, 66, 4488−4498. (14) Zhang, Y.; Ma, C. Fault diagnosis of nonlinear processes using multiscale KPCA and multiscale KPLS. Chem. Eng. Sci. 2011, 66, 64− 72. (15) Comon, P. Independent component analysis, a new concept. Signal Process. 1994, 36, 287−314. (16) Hyvärinen, A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Networks 1999, 10, 626−634. (17) Hyvärinen, A.; Oja, E. Independent component analysis: Algorithms and applications. Neural Networks 2000, 13, 411−430. (18) Lee, J. M.; Yoo, C. K.; Lee, I. B. Statistical process monitoring with independent component analysis. J. Process Control 2004, 14, 467−485. (19) Lee, J. M.; Qin, S. J.; Lee, I. B. Fault detection and diagnosis of multivariate process based on modified independent component analysis. AIChE J. 2006, 52, 3501−3514. (20) Wang, J.; Zhang, Y.; Cao, H.; Zhu, W. Dimension reduction method of independent component analysis for process monitoring based on minimum mean square error. J. Process Control 2012, 22, 477−487.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. 4337



Article

(21) Lee, J.-M.; Yoo, C.; Lee, I.-B. Statistical monitoring of dynamic processes based on dynamic independent component analysis. Chem. Eng. Sci. 2004, 59, 2995−3006. (22) Lee, J. M.; Qin, S. J.; Lee, I. B. Fault detection of non-linear processes using kernel independent component analysis. Can. J. Chem. Eng. 2007, 85, 526−536. (23) Ge, Z.; Yang, C.; Song, Z.; Wang, H. Robust Online Monitoring for Multimode Processes Based on Nonlinear External Analysis. Ind. Eng. Chem. Res. 2008, 47, 4775−4783. (24) Zhang, Y. Enhanced statistical analysis of nonlinear processes using KPCA, KICA and SVM. Chem. Eng. Sci. 2009, 64, 801−811. (25) Zhao, C.; Gao, F.; Wang, F. Nonlinear batch process monitoring using phase-based kernel-independent component analysis-principal component analysis (KICA−PCA). Ind. Eng. Chem. Res. 2009, 48, 9163−9174. (26) Stefatos, G.; Hamza, A. B. Dynamic independent component analysis approach for fault detection and diagnosis. Expert Syst. Appl. 2010, 37, 8606−8617. (27) Rashid, M. M.; Yu, J. Hidden Markov model based adaptive independent component analysis approach for complex chemical process monitoring and fault detection. Ind. Eng. Chem. Res. 2012, 51, 5506−5514. (28) Yang, Y.; Chen, Y.; Chen, X.; Liu, X. Multivariate industrial process monitoring based on the integration method of canonical variate analysis and independent component analysis. Chemom. Intell. Lab. Syst. 2012, 116, 94−101. (29) Rashid, M. M.; Yu, J. A new dissimilarity method integrating multidimensional mutual information and independent component analysis for non-Gaussian dynamic process monitoring. Chemom. Intell. Lab. Syst. 2012, 115, 44−58. (30) Fan, J.; Qin, S. J.; Wang, Y. Online monitoring of nonlinear multivariate industrial process using filtering KICA-PCA. Control Eng. Pract. 2014, 22, 205−216. (31) Fan, J.; Wang, Y. Fault detection and diagnosis of non-linear non-Gaussian dynamic processes using kernel dynamic independent component analysis. Inf. Sci. 2014, 259, 369−379. (32) Fan, J.; Wang, Y.; Qin, S. J. Combined indices for ICA and their applications to multivariate process fault diagnosis (in Chinese). Acta Autom. Sin. 2013, 39, 494−501. (33) Alcala, C. F.; Qin, S. J. Analysis and generalization of fault diagnosis methods for process monitoring. J. Process Control 2011, 21, 322−330. (34) Mendel, J. M. Tutorial on higher-order statistics (spectra) in signal processing and system theory: Theoretical results and some applications. IEEE Trans. Signal Process. 1991, 79, 278−305. (35) Giannakis, G. B.; Mendel, J. M. Cumulant-based order determination of non-Gaussian ARMA models. IEEE Trans. Acoust., Speech, Signal Process. 1990, 38, 1411−1423. (36) Giannakis, G. B.; Tsatsanis, M. K. Signal detection and classification using matched filtering and higher order statistics. IEEE Trans. Acoust., Speech, Signal Process. 1990, 38, 1284−1296. (37) Sadler, B. M.; Giannakis, G. B.; Lii, K.-S. Estimation and detection in non-Gaussian noise using higher order statistics. IEEE Trans. Signal Process. 1994, 42, 2729−2741. (38) Choudhury, M. A. A. S.; Shah, S. L.; Thornhill, N. F. Diagnosis of poor control-loop performance using higher-order statistics. Automatica 2004, 40, 1719−1728. (39) Downs, J. J.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245−255. (40) Lyman, P. R.; Georgakis, C. Plant-wide control of the Tennessee Eastman problem. Comput. Chem. Eng. 1995, 19, 321−331. (41) Chiang, L. H.; Russell, E.; Braatz, R. D. Fault detection and diagnosis in industrial systems; Springer: London, 2001. (42) Valle, S.; Li, W.; Qin, S. J. 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.

(43) Martin, E. B.; Morris, A. J. Non-parametric confidence bounds for process performance monitoring charts. J. Process Control 1996, 6, 349−358. (44) Downs, J. J.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245−255. (45) Lyman, P. R.; Georgakis, C. Plant-wide control of the Tennessee Eastman problem. Comput. Chem. Eng. 1995, 19, 321−331.

4338


Online Monitoring of Multivariate Processes Using Higher-Order

Recommend Documents