Combination Method of Principal Component and Wavelet Analysis for

Principal component analysis (PCA) has been widely used in multivariate process monitoring for its ability to reduce process dimensions. PCA and other...
0 downloads 0 Views 233KB Size
Download PDF
4198

Ind. Eng. Chem. Res. 2003, 42, 4198-4207

Combination Method of Principal Component and Wavelet Analysis for Multivariate Process Monitoring and Fault Diagnosis Ningyun Lu,†,‡ Fuli Wang,‡ and Furong Gao*,† Department of Chemical Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, and School of Information Science and Engineering, Northeastern University, Shenyang, Liaoning, China

Product quality and operation safety are important aspects of industrial processes, particularly those with large numbers of correlated process variables. Principal component analysis (PCA) has been widely used in multivariate process monitoring for its ability to reduce process dimensions. PCA and other statistical techniques, however, have difficulties in differentiating faults with similar time-domain process characteristics. A wavelet-based time-frequency approach is developed in this paper to improve PCA-based methods by extending the time-domain process features into time-frequency information. Subsequently, a similarity measure is presented to compare process features for on-line process monitoring and fault diagnosis. Simulation results show that the proposed multivariate time-frequency process feature is effective in both fault detection and diagnosis, illustrating the potentials for real-world application. 1. Introduction The demands for product quality and operation safety in the process industry have spurred the recent development and application of process monitoring and fault diagnosis methods. The traditional statistical process monitoring approaches, founded on Shewhart control charts, are ill-suited for modern chemical processes. Multivariate statistical methods based on correlation models such as principal component analysis (PCA) and partial least-squares (PLS) have been widely used to handle processes with large numbers of correlated process variables.1-10 Process monitoring based on these methods conducts statistical hypothesis tests on two indices, the Hotelling T2 and Q statistics in principal component and residual subspaces, respectively. These methods are effective in fault detection, but not diagnosis. Contribution plots,11 widely used as a diagnosis tool in PCA/PLS-based process monitoring, can only identify a group of variables impacted by a process malfunction. To improve fault diagnosis, many advances have been reported. According to Yoon and MacGregor,12 approaches for fault detection and diagnosis might be classified into three categories: methods based on causal models, methods based on knowledge, and methods based on multivariate statistics. Improved fault detection and diagnosis performance have been obtained by combining methods of different categories, particularly, a statistical method, for example, PCA, with other methods. Qin and his colleges13-15 have extended PCA for sensor fault identification and reconstruction for both unidimensional and multidimensional cases. The integration of PCA and discriminant analysis was proposed by Raich and Cinar16,17 to improve PCA’s ability in diagnosis. PCA was combined with parity relation by Gertler et al.,18 to employ concepts of analytical redundancy for improving isolation properties. The PCA-based contribu* To whom correspondence should be addressed. Tel.: +8522358-7139. Fax: +852-2358-0054. E-mail: [email protected]. † Hong Kong University of Science and Technology. ‡ Northeastern University.

tion plot was enhanced by Vedam and Venkatasubramanian19 using signed digraphs. The combination of PCA with expert knowledge has also been reported for fault diagnosis improvement for PCA-based methods.20,21 All of the above methods can enhance the fault diagnosis of PCA-based approaches, provided that the fault features are differentiable in the time domain. Certain process faults might show similar time-domain features because of an underlying correlation and/or the employed closed-loop controls. In this case, it would be desirable to extend the time-domain information into the time-frequency domain to improve fault identification. It is well-known that relations between process conditions and frequency characteristics do, indeed, exist. Changes in the spectrum characteristics over an operation period can provide useful information for fault diagnosis. It is, therefore, possible to utilize the frequency-domain information to improve PCA-based process monitoring and diagnosis.22-27 Wavelet transformation is a powerful tool for transforming time-domain signals into the time-frequency domain. Combinations of wavelet transforms and PCA have begun to emerge recently. For example, process data can be calibrated first by a Harr wavelet transform (HWT) before PCA is applied.22 PCA was extended to wavelet-based multiscale PCA (MSPCA) by Bakshi,23 to capture information in both the time and frequency domains. This results in improvements in fault detection. Using MSPCA, multiscale fault identification was proposed by Misra et al.27 to analyze contribution plots at the corresponding scales. This paper proposes to employ wavelet transforms to improve PCA-based fault identification for cases where process faults with similar time-domain features can be differentiated by their frequency information. The highdimensional and correlated process data are first preprocessed by PCA to be projected onto the low-dimensional principal component subspace. The principal component scores are transformed into the timefrequency domain, and decomposed into an approximation signal and the corresponding detail signals of different scales. The energy distribution over selected

10.1021/ie0207313 CCC: $25.00 © 2003 American Chemical Society Published on Web 07/30/2003

Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003 4199

scales is defined as a feature pattern to be used for process monitoring and fault diagnosis. A similarity measure, introduced to differentiate fault features from normal process features, is also employed for feature mapping in diagnosis. The remainder of this paper is organized as follows. Section 2 briefly introduces the techniques of PCA and wavelet analysis. Section 3 first presents the definition of the process feature pattern to represent the multivariate time-frequency information, followed by the development of the proposed process monitoring and diagnosis approach. Application examples of the proposed scheme are demonstrated and discussed in section 4. Finally, conclusions are given in section 5.

ψj,k(t) ) 2-jψ(2-jt - k), j ) 1, ..., J; k ∈ Z

(5)

where j is the scale factor and k is the translation factor. The approximation coefficients a(J,k) and the detail coefficients b(j,k) can be computed by the Mallat algorithm.31,32 The approximation signal AJ(t) and the detail signals Dj(t) (j ) 1, 2, ..., J) are defined as

AJ(t) ) Dj(t) )

a(J,k) φJ,k(t) ∑ k∈z

(6)

∑ b(j,k) ψj,k(t)

(7)

k∈Z

2. Background 2.1. Principal Component Analysis (PCA). PCA28 reduces the redundant information in a data block by projecting original process measurements onto a lowerdimensional subspace defined by a few orthogonal latent variables that contain most of the variance of the original data. Let X ) (x1, x2, ..., xm) be an n × m-dimensional data set. X can be decomposed as

AJ(t) and Dj(t) are constructed from the approximation coefficients a(J,k) and the detail coefficients b(j,k), respectively, representing the information on the corresponding scales. The original signal is actually the sum of AJ(t) and the Dj(t)’s J

f(t) ) AJ(t) +

Dj(t) ∑ j)1

(8)

m

X ) TPT )

tipiT ∑ i)1

(1)

where pi is defined as the principal component loading vector and ti is the corresponding score vector. Less important components, which mostly describe noise information in the data, can be abandoned without the loss of significant information.8 By doing so, matrix X can be reconstructed as the sum of an estimation value, X ˆ , and its residual, E k

X)X ˆ + E, with X ˆ )

tipiT ∑ i)1

(2)

where k represents the number of principal components retained. The principal component score vectors, t1, ..., tk, span a lower-dimensional subspace used for further analysis. Process monitoring and fault diagnosis is conducted via statistical hypothesis tests on two indices, Hotelling T2 and Q statistics, in the principal component subspace and residual subspace, respectively. Details on these methods can be found in many references.1-10 2.2. Wavelet Analysis. Wavelet analysis,29,30 a powerful signal-processing tool, can transform the timedomain signals into the time-frequency domain. Under the framework of multiresolution analysis, the original signal space V0 can be decomposed into a hierarchical set consisting of an approximation space VJ [the space spanned by scale functions φJ,k(t), k ∈ Z] and detail spaces Wj [the spaces spanned by wavelet functions ψj,k(t), j ) 1, ..., J, k ∈ Z], where J is the coarsest scale, normally called the decomposition level. According to the above description, a signal f(t) ∈ L2(R) can be decomposed as follows J

f(t) )

∑ a(J,k)φJ,k(t) + ∑ ∑ b(j,k) ψj,k(t) k∈Z j)1k∈Z

φJ,k(t) ) 2-Jφ(2-Jt - k), k ∈ Z

(3) (4)

Wavelet analysis uses the scale factor j instead of the frequency ω. In the decomposition procedures, the “frequency” is scaled down by 2j. Assuming that the original signal has a frequency of 1, the initial space V0 is ultimately decomposed into the low-frequency subspace VJ with frequency band [0, 1/2J] and the highfrequency subspaces Wj with frequency band [1/2j, 1/2j-1]. This indicates that a time-domain signal can be decomposed by wavelet analysis into a sum of components of different frequency bands. The energy over certain frequency bands can represent certain process information, and changes in the energy distribution over scales can provide useful information for fault diagnosis. 3. Process Monitoring and Diagnosis Based on a Combination of PCA and Wavelet Analysis 3.1. Overview. Industrial processes, particularly chemical processes, are usually complicated, with many correlated process variables. It is desirable and computationally economical to reduce the process to a lowerdimensional space before applying a process monitoring and diagnosis scheme. PCA-based process monitoring and diagnosis is enhanced by wavelet transform by extending the time-domain fault features into the timefrequency domain, for better fault interpretation and differentiation. The question is then how to effectively extract quantitative feature patterns from wavelet coefficients of the principal component scores. As discussed before, the time-domain signal can be decomposed as the sum of its approximation and detail signals, each of which contains information in a corresponding frequency band. The energy of the original signal can also be decomposed as the sum of the energies of the approximation and details. The energy distribution over the selected frequency bands can be regarded as a quantitative feature. Process features can then be defined in terms of the energy distributions of the principal components scores. To compare two process features, a similarity measure is subsequently introduced for process monitoring and diagnosis.

4200 Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003

The proposed method is limited to multivariate continuous processes with “time-invariant” normal process features, that is, with time-invariant means and variances. Such processes have stable energy distributions at normal operation,33 which makes it possible to conduct on-line process monitoring by comparing the current process feature with its historical behavior. For time-varying processes, such as those encountered in batch processes, the proposed method would have to be improved. This is one focus of our future research. 3.2. Quantitative Feature Extraction by Wavelet Analysis. For a signal, f(t) ∈ L2(R), in the time domain, the energy can be defined as

Ef )

∫-∞∞f 2(t) dt

(9)

As discussed in section 2.2, f(t) can be decomposed as the sum of an approximation signal, AJ, and a set of detail signals, Dj (j ) 1, ..., J). AJ contains information in frequency band [0, 1/2J], whereas Dj contains information in frequency band [1/2j, 1/2j-1]. Normally, the finest-detail scale contains noise information. The energies of AJ and Dj take the form24

EA J ) EDj )

comparison of the two process features. For a process with only one principal component retained, the extracted process feature is the normalized energy distribution of the principal component score. The similarity measure can be defined as

Scurr,normal )

(Ξcurr)T(Ξnormal) ||Ξcurr||‚||Ξnormal||

(13)

which, in fact, takes the form of the cosine of the vectors’ angle. The energy distribution is a vector of positive elements; the similarity will fall within the range (0,1]. For a process with k principal components [t1, t2, ..., tk], the similarity measure is defined as

Scurr,normal ) g1St1,normal + ‚‚‚ + gkStk,normal Sti,normal )

(Ξti)T(Ξnormal) ||Ξti||‚||Ξnormal||

(14)

where

N/2

1

a(J,k + m)2 ∑ Nm)-N/2+1 1

N/2

∑

Nm)-N/2+1

(10)

b(j,k + m)2

where N is the window length. The energy distribution can be regarded as a feature pattern that contains both time- and frequency-domain information. With a given decomposition level J, the normalized energy distribution, E h f, is defined as a quantitative feature as follows

Ef ) [ED1, ..., EDJ, EAJ] E h f ) Ef/||Ef||

(11)

Relations exist between a process fault and its energy distribution. For example, for a slowly drifting signal, the energy in the low-frequency subspaces, particularly in the approximation subspace AJ, will be significant. High-frequency signal changes will result in changes of the energy in the finer-detail subspaces. For a signal with a cyclic disturbance, the energy in the corresponding scale will be significant. Analyzing changes in the energy distribution can be useful for fault identification. For a process having k principal components, t1, ..., tk, the process feature, Ξ, is a set of the quantitative feature patterns for each retained principal component

h t2, ..., E h tk]T Ξ ) [E h t 1, E

(12)

3.3. On-line Process Monitoring and Fault Diagnosis. As for existing statistical process monitoring methods, a reference data set, X, must be collected from the normal operation of the process for modeling. The normal feature is then computed for comparison on-line with the feature of the evolving process. A moving window is used for feature extraction and comparison in the on-line implementation. Process monitoring is conducted by comparing the current process feature, Ξcurr, with the normal process feature, Ξnormal. A similarity measure is defined for the

gi )

λi k

, i ) 1, ..., k

λj ∑ j)1 and λi represents the eigenvalues of XTX in descending order. This similarity measure also falls within the range (0,1]. The maximum similarity of 1 indicates an exact match. The closer a similarity value is to 0, the more unlikely two faults are to be similar. It should be noted that the process data in the moving window might vary with time because of noise, resulting in small variations of the process features from window to window even under normal operation. A similarity threshold, S*, is introduced to reduce the possibility of false alarms. With a similarity degree above the predefined threshold, the two process features are assumed to represent the same operating condition. The similarity threshold is set to allow 99% of normal data in the training data set above the threshold. This corresponds to a level of significance R ) 0.01 considering the probability distribution of the similarity degree for normal data.34 The threshold set in this way can reduce the chances of false alarms, as it is insensitive to normal process noises. It is, however, sensitive to the changes of measurement noise patterns, because the energy distribution covers up to the high-frequency bands, as illustrated in the second example in section 4. Fault diagnosis is, in fact, a procedure of pattern matching, by comparing the current fault feature with the existing ones in the fault database. The current fault is assumed to match an existing fault in the database if the similarity value is above the predetermined threshold. If the current fault feature matches none in the database, a new fault is deemed to have occurred. To summarize, the proposed process monitoring and diagnosis scheme consists of modeling and on-line twostep monitoring and diagnosis, as listed below: Modeling Procedure. 1. Select a proper wavelet function and determine the decomposition level, J, and the window length N (see discussion in section 4).

Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003 4201

2. Collect a period of process data under normal operation; develop the PCA model and calculate the principal component scores. 3. Form the normal process feature, Ξnormal, by applying the wavelet transform to the principal component scores. 4. Form the fault features, Ξf(i) (i ) 1, ..., Nf), where Nf is the number of known faults, following steps 2 and 3 for fault operations. On-line Process Monitoring and Diagnosis. 5. Compute the principal component scores of the process data in the evolving moving windows. Apply the wavelet transform to the scores and calculate the current process feature, Ξcurr. 6. Compare the current feature with the normal feature. If the similarity value, Scurr,normal, is above the predetermined threshold, S*, continuing the process monitoring by going back to step 5; otherwise, go to step 7. 7. Compare the current process feature with the existing ones in the fault database. If the similarity with an existing fault feature is above the threshold, the diagnosis procedure ends. If the current process feature does not match any existing one, a new fault type has occurred. The database is updated by extending it with the new fault feature, Ξcurr, and a description of this new feature. 4. Illustrative Examples The proposed process monitoring and diagnosis approach is tested with two processes, a coupled threetank process and the well-known Tennessee Eastman (TE) process. The first example demonstrates the ability of the proposed algorithm to differentiate faults with similar time-domain characteristics, whereas the second example focuses on the application of the proposed scheme to a more “complex” industrial process. The proper wavelet function, decomposition level, and moving window length have to be determined prior to the application of the present method. The Daubechies (dbN) wavelet is used in our work for its nice properties.33 The decomposition level is determined by the signal’s frequency bandwidth. The signals with abundant high-frequency information need larger numbers of decomposition levels. For the two illustrated examples, the decomposition level is set to be 5. The moving window method results in new problems: determination of the window length and the “border distortion” for finite-length signal. Because of dyadic down-sampling, the length of the wavelet coefficients is being reduced by a factor of 2j. To ensure the accuracy of signal reconstruction, it is desirable to limit the minimum length of the coefficients in the coarsest level J. With a minimum length of Lmin J, the length of the window for the original signals should be larger than 2J × Lmin J.23,27 Methods are also available to treat border distortion.29,32 Simpler schemes based on signal extension at the boundaries are preferred as they involve the computation of fewer extra coefficients at each stage of the decomposition. Symmetric extension32 is used here. This work focuses on the application of wavelet transforms, rather than on wavelet research itself. Further information on wavelet function selection, border distortion treatment, and so on can be found in the literature.29,32,33 For the first example, the energy distribution is defined as [ED1, ED2, ED3, ED4, ED5], covering all of the

Figure 1. Schematic of three-tank process.

detail subspaces. The approximate subspace is not included, because the simulated faults of the three-tank process have similar time responses. Inclusion of the approximate subspace will result in a single large component in the energy feature vector and, hence, reduce the differentiability of those process faults. For the TE process, the energy distribution includes both the detail subspaces and the approximate subspace, [ED1, ED2, ED3, ED4, ED5, EA5]. The selection of the frequency bands should consider the fault properties. In general, if the time-frequency features of two (or more) faults become difficult to identify because of their overwhelming high energy in a particular frequency band, it is recommended that that particular frequency band be temporarily omitted and the differentiation be focused on comparing information in the remaining frequency bands. 4.1. Three-Tank Process. A three-tank system, as shown in Figure 1, is used to demonstrate the proposed process monitoring and diagnosis. In this coupled system, the flow rates of the two pumps serve as the manipulated variables in controlling the levels of tanks 1 and 2, whereas the level of tank 3 is left to float to reflect the interaction between tanks 1 and 2. A decoupling control strategy has been developed for the simultaneous closed-loop control of the levels of tanks 1 and 2. Two types of typical faults are simulated: leakage faults caused by the opening of the the leakage valve, SLi (i )1, 2, 3), at the bottom of each tank and sensor failures caused by programming. Electrical white noises are purposely introduced into all measurements in the simulations. Five sensors are installed with this process: two flow sensors and three level sensors. They are sampled every second. All process measurements are first normalized with zero mean and unit variance. Three principal components are retained after analysis based on the widely accepted cross-validation algorithm.35 The normal process feature, extracted from the normal process data in the training data sets, is shown in Figure 2, where the energy distribution over five selected frequency bands is plotted against the three retained principal components. All together, nine faults are introduced, as described in Table 1, with the process features shown in Figure 3. Fault features are calculated from the data in the time span, (Tf - N/2, Tf + N/2), where Tf is the fault occurrence time for the reference data sets in the fault database and N is the moving window’s length. The similarity degrees between

4202 Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003

Figure 2. Normal process feature for the three-tank process. Table 1. Simulated Nine Faults for the Three-Tank Process fault

description of simulated fault

1 2 3 4 5 6 7 8 9

leakage of tank 1 leakage of tank 2 leakage of tank 3 leakage of tanks 1 and 2 leakage of tanks 1 and 3 leakage of tanks 2 and 3 sensor failure of tank 1 sensor failure of tank 2 sensor failure of tank 3

pairs of process faults are calculated by eq 14. The threshold of similarity is set to 0.975, corresponding to the level of significance R ) 0.01, as discussed in section 3. Table 2 shows the differentiability of the proposed monitoring and diagnosis method.

The results obtained by the PCA-based contribution plot method are also presented for comparison in Figure 4 and Table 3. The contribution plots have difficulty in differentiating the first, fourth, seventh, and ninth faults, all of which actually lead to an increase in Q1 and a decreasing level in tank 1. Similarly, the contribution plots cannot differentiate the second, third, sixth, and eighth faults, all of which lead to an increase in Q2 and a decreasing level in tank 2. The bold values in Table 3 indicate the cases that are difficult to differentiate on the basis of their contribution plots. A comparison between Figures 3 and 4, or Tables 2 and 3, indicates that the proposed method outperforms the contribution plot, with the superiority of the proposed method lying in its time-frequency analysis. The use of the energy distribution to represent process features is effective for differentiating faults with similar time-domain characteristics. 4.2. Tennessee Eastman Process. The proposed scheme was also applied to a well-known benchmark process, Tennessee Eastman (TE) industrial process, which is widely used for comparison of the various process monitoring methods. The TE process was first introduced by Downs and Vogel.36 The MATLAB simulation procedure follows that of Ashish37 with the plantwide control strategy of McAvoy and Ye.38 The TE process is schematically shown in Figure 5. The simulation for this paper is based on the operation in the base case, corresponding to a 50/50 G/H product ratio. Among the 41 measurements and 12 manipulated variables, 52 are selected (all but the agitation speed of the reactor’s stirrer). The sampling interval is set to be 3 min, the same as the recommended value in ref 34. Data collected from 96 h of normal operation are used to develop the PCA model. Six principal components are

Figure 3. Process features for the simulated nine faults (x axis is the process variable).

Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003 4203

Figure 4. Contribution plots for the simulated nine faults (x axis is the process variable). Table 2. Similarity Degrees between Fault Features Represented by Energy Distribution for the Three-Tank Process fault fault

1

2

3

4

5

6

7

8

9

1 2 3 4 5 6 7 8 9

1.0000 -

0.9022 1.0000 -

0.8803 0.8742 1.0000 -

0.9198 0.9610 0.8960 1.0000 -

0.7967 0.8722 0.9275 0.8688 1.0000 -

0.9101 0.8519 0.8158 0.9014 0.7060 1.0000 -

0.7184 0.6343 0.8351 0.6773 0.7864 0.7429 1.0000 -

0.9617 0.9479 0.9324 0.9614 0.8864 0.8581 0.7137 1.0000 -

0.8938 0.9519 0.7917 0.9425 0.7458 0.8162 0.4664 0.9303 1.0000

Table 3. Similarity Degrees between Fault Features Represented by Contribution Rate to SPE for the Three-Tank Process fault fault

1

2

3

4

5

6

7

8

9

1 2 3 4 5 6 7 8 9

1.0000 -

0.0924 1.0000 -

0.1031 0.9997 1.0000 -

0.9994 0.0606 0.0711 1.0000 -

0.9549 0.3781 0.3893 0.9444 1.0000 -

0.2828 0.9776 0.9774 0.2531 0.5434 1.0000 -

0.9969 0.0210 0.0310 0.9900 0.9292 0.2164 1.0000 -

0.0963 1.0000 0.9995 0.0646 0.3813 0.9789 0.0251 1.0000 -

0.9983 0.0386 0.0488 0.9997 0.9364 0.2326 0.9998 0.0427 1.0000

retained, capturing major correlations in the process variables. The normal process feature is extracted from the 96 h of training data. Additional data are generated with three different process upsets, IDV(1), IDV(4), and IDV(11). These three process faults are the typical cases used in the work of Chiang et al.,34 where detailed reviews and comparisons were conducted with various methods of process monitoring and fault diagnosis (PCA, DPCA, CVA, etc.) on the TE process. These three runs start with no faults, and the faults are introduced at

the 36th simulation hour. Each fault simulation operates for about 48 h, and 960 measurements are collected for monitoring. For normal operation, a “steady-state” condition, measurements show little autocorrelation, most of which can be reflected in the first few principal components. The energy distributions of the first few principal components show high energy in the low-frequency bands, whereas those of the remaining principal components show higher energy in the higher-frequency

4204 Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003

Figure 5. Tennessee Eastman process.

Figure 6. Normal process feature for the TE process.

bands, as illustrated in Figure 6. With closed-loop control, most of upsets can introduce significant autocorrelation among process measurements. The process features represented by the proposed energy distribution can change correspondingly in different manners. This is the basis of the proposed process monitoring method. For on-line application, the minimum length of the coefficients in the coarsest level is set to be Lmin5 ) 8. The moving window’s length is therefore 25 × 8 ) 256. The moving step is set to 16, so that the first 16 data in the old window are deleted and 16 new data are inserted in the new window. Principal component scores for measurements in each moving window are subsequently processed by wavelet analysis to capture the current process feature, which is then compared to the normal

Figure 7. On-line monitoring for normal data.

process feature. Figure 7 shows the similarity values of the normal process data in the training data set, which will be used to set the threshold for on-line monitoring and diagnosis. As shown in Figure 7, the similarity varies with time because of process noise. The threshold of similarity is set to eliminate any false alarms caused by the normal noise. There are 104 moving windows for the normal training data, and the threshold is set to be the second lowest similarity value, which corresponds approximately to the level of significance R ) 0.01, that is, S* ) 0.9382 is selected for the following process monitoring and fault diagnosis. If the similarity value for the evolving window is less than S*, an alarm will be given to indicate an abnormality. Figure 8 shows the on-line monitoring charts for faults IDV(1), IDV(4), and IDV(11). All three faults can

Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003 4205

Figure 8. On-line monitoring charts for faults IDV(1), IDV(4), and IDV(11).

Figure 9. On-line diagnosis of fault IDV(4).

be detected by the proposed method, shortly after their occurrence at 36 h. Comparisons with the fault detection performance of other methods are conducted as in the work of Chiang et al.,34 where the missed detection rate is used as the criterion. Table 4 gives the missed detection rates of faults IDV(1), IDV(4), and IDV(11) for PCA, DPCA, CVA, and the proposed multivariate timefrequency monitoring methods. According to Table 4, the

proposed method performs better than the others, with the lowest missed detection rates among all compared methods. Figures 9 and 10 illustrate the diagnosis ability of the proposed method for faults IDV(4) and IDV(11), respectively. Fault IDV(4) is caused by a step change in the reactor cooling water temperature, whereas fault IDV(11) is caused by a random variation in the reactor

4206 Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003

Figure 11. Features of faults IDV(4) and IDV(11). Table 4. Missed Detection Rates for Faults IDV(1), IDV(4), and IDV(11)a method

statistics

IDV(1)

IDV(4)

IDV(11)

Q

0.008 0.003

0.956 0.038

0.794 0.356

DPCA

T2 Q

0.006 0.005

0.939 0

0.801 0.193

CVA

Ts2 Tr2 Q

0.001 0 0.003

0.688 0 0.975

0.515 0.196 0.669

0

0

0.133

T2

PCA

proposed method a

similarity

Compared with the results of Chiang et

al.34

Table 5. Similarity Degrees between Features of Faults IDV(1), IDV(4), and IDV(11)

Figure 10. On-line diagnosis of fault IDV(11).

cooling water temperature. The closed-loop control system makes these two faults look similar in timedomain behavior, which results in most existing fault diagnosis methods failing to differentiate them.34 Figures 9 and 10 show that the proposed time-frequency method can differentiate these two faults without any problem, because they have different energy distributions, as shown in Figure 11. For fault IDV(4), the energy distributions weight to low-frequency bands,

Ξnormal IDV(1) IDV(4) IDV(11)

Ξnormal

IDV(1)

IDV(4)

IDV(11)

1.0000 -

0.2948 1.0000 -

0.7327 0.6674 1.0000 -

0.8748 0.4363 0.7426 1.0000

whereas the energy distributions of fault IDV(11) weight to mid-frequency bands, a difference that can readily be explained by the different disturbance types. Table 5 shows the similarity degrees for the three selected faults. The on-line diagnosis charts, as shown in Figures 9 and 10, actually monitor the similarity degrees between

Ind. Eng. Chem. Res., Vol. 42, No. 18, 2003 4207

the evolving process feature and the known fault features in the fault database. If the similarity between the evolving feature and the feature of a certain fault, Ξf(i), become larger than the threshold, the fault diagnosis procedure declares the occurrence of that fault. The on-line diagnosis charts of Figure 9 clearly show the occurrence of fault IDV(4), whereas Figure 10 clearly demonstrates the diagnosis of fault IDV(11). Fault IDV(11) is caused by random variations, a pattern similar to measurement noise. This example indicates that a proper similarity threshold, introduced to reduce false alarms, does not reduce sensitivity. 5. Conclusions A new process monitoring and diagnosis method has been proposed that combines principal component analysis to reduce the process dimension and wavelet analysis to extract the time-frequency process feature. The process feature is based on the energy distribution at the selected frequency bands, and a similarity measure is defined to compute the similarity degrees between pairs of process features. Simulations conducted on the three-tank system and the TE process have shown that the proposed method can not only detect abnormal conditions, but also differentiate faults with similar time-domain characteristics, overcoming the major difficulty associated with the existing contribution plot method. The proposed method exploits the advantages of both statistical methods and frequency methods, extracts process features from time-frequency information, and provides an effective tool for process monitoring and diagnosis. Literature Cited (1) Kresta, J. V.; MacGregor, J. F.; Marlin, T. E. Multivariate statistical monitoring of process operating performance. Can. J. Chem. Eng. 1991, 69, 35. (2) MacGregor, J. Statistical process control of multivariate processes. Control Eng. Pract. 1994, 3, 403. (3) MacGregor, J. F.; Jaeckle, C.; Kiparissides, C.; Koutoudi, M. Process monitoring and diagnosis by multiblock PLS methods. AIChE J. 1994, 40, 826. (4) Nomikos, P.; MacGregor, J. F. Monitoring batch processes using multiway principal component analysis. AIChE J. 1994, 40, 1361. (5) Nomikos, P.; MacGregor, J. F. Multivariate SPC charts for monitoring batch process. TECHNOMETRICS 1995, 37, 41. (6) Martin, E. B.; Morris, A. J.; Zhang J. Process performance monitoring using multivariate statistical process control. IEE Proc.-Control Theory Appl. 1996, 143, 132. (7) Wise, B. M.; Gallagher, N. B. The process chemometrics approach to process monitoring and fault detection. J. Process Control 1996, 6, 329. (8) Zhang, J.; Martin, E. B.; Morris, A. J. Fault detection and diagnosis using multivariate statistical process control. Chem. Eng. Res. Des. 1996, 74, 89. (9) Chen, G.; McAvoy, T. J. Predictive on-line monitoring of continuous processes. J. Process Control 1998, 8, 409. (10) Chen, J.; Liu, K. C. On-line batch process monitoring using dynamic PCA and dynamics PLS models. Chem. Eng. Sci. 2002, 57, 63. (11) Miller, P.; Swanson, R. E.; Heckler, C. E. Contribution plots: A missing link in multivariate quality control. Appl. Math. Comput. Sci. 1998, 8, 775. (12) Yoon, S.; MacGregor, J. F. Statistical and causal modelbased approaches to fault detection and isolation. AIChE J. 2000, 46, 1813. (13) Dunia, R.; Qin, S. J.; Edgar, T. F.; McAvoy, T. J. Identification of fault sensors using principal component analysis. AIChE J. 1996, 42, 2797.

(14) Dunia, R.; Qin, S. J. A unified geometric approach to process and sensor fault identification and reconstruction: The unidimensional fault case. Comput. Chem. Eng. 1998, 22, 927. (15) Dunia, R.; Qin, S. J. Subspace approach to multidimensional fault identification and reconstruction. AIChE J. 1998, 44, 1813. (16) Raich, A.; Cinar, A. Multivariate statistical methods for monitoring continuous processes: Assessment of discrimination power of disturbance models and diagnosis of multiple disturbances. Chemom. Intell. Lab. Syst. 1995, 30, 37. (17) Raich, A.; Cinar, A. Diagnosis of process disturbances by statistical distance and angle measures. Comput. Chem. Eng. 1997, 21, 661. (18) Gertler, J.; Li, W.; Huang, Y.; McAvoy, T. Isolation enhanced principal component analysis. AIChE J. 1999, 45, 323. (19) Vedam, H.; Venkatasubramanian, V. PCA-SDG based process monitoring and fault diagnosis. Control Eng. Pract. 1999, 7, 903. (20) Norvilas, A.; Negiz, A.; Decicco, J.; Cinar, A. Intelligent process monitoring by interfacing knowledge-based systems and multivariate statistical monitoring. J. Process Control 2000, 10, 341. (21) Leung, D.; Romagnoli, J. An integration mechanism for multivariate knowledge-based fault diagnosis. J. Process Control 2002, 12, 15. (22) Kosanovich, K. A.; Piovoso, M. PCA of wavelet transformed process data for monitoring. Intel. Data Anal. 1997, 1, 85. (23) Bakshi, B. R. Multiscale PCA with application to multivariate statistical process monitoring. AIChE. J. 1998, 40, 1596. (24) Zhao, Z.; Jiang, W. S.; Gu, X. S. Process monitoring based on wavelet Transform. Control Decision 1999, 14, 19. (25) Daiguji, M.; Kudo, O.; Wada, T. Wavelet-based fault detection and identification in oil refinery. Presented at the 14th IFAC World Congress, Beijing, China, Jul 4-9, 1999. (26) Zhang, J. Q.; Yan, Y. A wavelet-based approach to abrupt fault detection and diagnosis of sensors. IEEE Trans. Instrum. Meas. 2001, 50, 1389. (27) Misra, M.; Yue, H. H.; Qin, S. J.; Ling, C. Multivariate process monitoring and fault diagnosis by multi-scale PCA. Comput. Chem. Eng. 2002, 26, 1281. (28) Jackson, J. E. A User’s Guide to Principal Components; John Wiley & Sons: New York, 1991. (29) Strang, G.; Nguyen, T. Wavelets and Filter Banks; Wellesley-Cambridge Press: Wellesley, MA, 1996. (30) Mallat, S.; Hwang, W. L. Singularity detection and processing with wavelets. IEEE Trans. Inf. Theory 1992, 38, 617. (31) Mallat, S. A theory of multiresolution signal decomposition: The wavelet representation. IEEE Trans. Pattern Anal. Machine Intell. 1989, 11, 674. (32) Misiti, M.; Misiti, Y.; Oppenheim, G.; Poggi, J. M. Wavelet Toolbox for Use with MATLAB; The Mathworks, Inc.: Natick, MA, 1997. (33) Motard, R. L.; Joseph, B. Wavelet Applications in Chemical Engineering; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1994 (34) Chiang, L. H.; Russell, E. L.; Braatz, R. D. Fault Detection and Diagnosis in Industrial Systems; Springer: Hong Kong, 2001. (35) Wold, S. Cross-validatory estimation of the number of components in factors and principal component models. TECHNOMETRICS 1978, 20, 397. (36) Downs, J. H.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245. (37) Singhal, A. Tennessee Eastman plant simulation with base control system of McAvoy and Ye. Available at http://chemengr.ucsb.edu/∼ceweb/computing/TE/tesim.pdf, accessed in May 2002. (38) McAvoy, T. J.; Ye, N. Base level control for the Tennessee Eastman problem. Comput. Chem. Eng. 1994, 18, 383.

Received for review September 16, 2002 Revised manuscript received April 15, 2003 Accepted June 14, 2003 IE0207313