Ind. Eng. Chem. Res. 2007, 46, 8033-8043
8033
Batch Process Monitoring in Score Space of Two-Dimensional Dynamic Principal Component Analysis (PCA) Yuan Yao and Furong Gao* Department of Chemical Engineering, Hong Kong UniVersity of Science and Technology, Clear Water Bay, Kowloon, Hong Kong SAR, People’s Republic of China
Two-dimensional dynamic principal component analysis (2-D-DPCA) is a recent developed method for twodimensional (2-D) dynamic batch process monitoring. However, it only utilizes residual information in fault detection and information in score space is wasted, which may compromise the monitoring efficiency. In this paper, 2-D multivariate score autoregressive (AR) filters are designed to remove the 2-D dynamics retained in score space and make the filtered scores obey certain statistical assumptions, so that the T2 statistic can be calculated reasonably for process monitoring. Simulation shows that using the filters enhances the monitoring efficiency while reducing the chances of false alarms and missed alarms. 1. Introduction Batch processes widely exist in chemical and other manufacturing industries, such as polymerization, semiconductor processing, injection molding, and most bioprocesses, because of its flexibility with the production of high-value-added products. In response to the challenges of operation safety and product quality requirements, many research efforts have been devoted to batch process monitoring.1-4 Process monitoring can be based on process knowledge, first-principle, or historical data. Among them, the data-based multivariate statistical process monitoring methods, such as principal component analysis (PCA) and partial least-squares (PLS), are widely used in cases where process mechanisms are complicated and process knowledge is incomplete. Many researches have extended PCA/PLS methods from continuous monitoring processes to batch processes monitoring.5-9 Traditionally, multivariate statistical monitoring methods for batch process are developed with the assumption of batch independency. However, this assumption may not hold for many practical batch processes. For example, the slow property changing of feed stocks, the drift of process characteristics, process controllers designed in the manner of run-to-run adjustment, and the effects of slow response variables (observed in many cases) each results in dynamics, not only within a batch, but also from batch to batch. In response to that observation, several modeling and monitoring methods have been developed to take batch-to-batch dynamics into consideration. For example, the lifted state space model10 and the method that incorporates information from prior batches into multiway PCA (MPCA) model building11 are concerned with batchwise dynamics. For a batch process, the dynamics can be viewed with two different time axes: a within-batch time axis and a cross-batch time axis. Based on such an idea, a two-dimensional dynamic principal component analysis (2-D-DPCA) model has been developed, with within-batch time direction as one dimension, and batch direction as the second dimension.12 With a proper region of support (ROS), this method can provide efficient faults detection on-line. Even small changes in correlation or process drifts can be detected effectively. Most recently, an automatic ROS selection algorithm for the method has been developed by the authors.13 However, the 2-D-DPCA method that was proposed previously can only monitor the process in the residual space, and information in the score space was not explored. This may affect the efficiency of the fault detection.
As a complement to the 2-D-DPCA that we developed, we extend the monitoring to score space in this paper. The article is organized as follows. In the next section, process monitoring in the residual space and score space of the PCA model is reviewed. Section 3 gives a review of the 2-D-DPCA method previously proposed. In section 4, the problems and difficulties in score space monitoring of 2-D-DPCA model are described. In section 5, solutions to these problems then are provided in detail. In section 6, the complete procedures of batch process monitoring based on the 2-D-DPCA model, including both score space monitoring and residual space monitoring, are developed. In section 7, simulations are given to compare the monitoring results of the T2 statistics of 2-D-DPCA model with and without 2-D multivariate score autoregressive (AR) filtering. The comparison between the monitoring based on SPE and T2 control plots are also conducted in this section. Finally, conclusions are drawn.
* To whom correspondence should be addressed. Tel.: (852)-23587139. Fax: (852)-2358-0054. E-mail address:
[email protected].
T2 ) tTΛ-1t
2. Process Monitoring of PCA-Based Methods 2.1. Fundamentals of PCA. The PCA analysis divides the original data into two subspaces: score space, which extracts most systematical variance information, and residual space, which is expected to contain only noise. In each subspace, a statistic can be computed for process monitoring purposes. In score space, such a statistic is called T2, which is a measure of the distances between new operation points and the historical normal operation region. In residual space, this statistic is the squared prediction error (SPE), which is a measurement of model fit. Abnormal process variations (i.e., faults) can be detected by comparing their values to the corresponding control limits that have been determined from the historical data. These two statistics are concerned with two different types of information and are complementary to each other. The simultaneous use of both of them in process monitoring can avoid potential missed alarms. Mathematically, PCA is shown as A
X ) TP ) T
m
tjpj ∑ j)1
T
+
tjpjT ) Xˆ + E ∑ j)A+1
(1)
where X(m × n) is the data matrix; n is the number of samples; m is the number of variables; tj(n × 1) is the score vector; pj(n × 1) is the loading vector (which projects data into score space); T and P are the score matrix and loading matrix, respectively; and E is the residual matrix. SPE is calculated from E, and T2 can be calculated as
10.1021/ie070579a CCC: $37.00 © 2007 American Chemical Society Published on Web 11/01/2007
(2)
8034
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
where Λ is a diagonal matrix that contains the first n-largest eigenvalues of the covariance matrix of t. The control limits of SPE and T2 are calculated statistically, based on some assumptions. SPE requires a normal distribution of residuals, i.e., residuals contain only noises without any systematic variance or dynamic information. The calculation of the control limits of T2 is based on the assumption of a normal distribution of each score, which indicates an independence of time. If such assumptions are not obeyed, the calculations of control limits become somewhat unreasonable, and the efficiency and correctness of fault detection is compromised. The details of PCA can be found in the work by Jackson.14 2.2. Dynamic Process Monitoring Based on Extensions of PCA. The traditional PCA cannot handle dynamic process data well, because, in such situations, the score space of PCA cannot extract all of the systematic variance information and some portions of the dynamics are left behind. This breaks the statistical basis for control limits calculation of SPE. At the same time, dynamic information also is contained in scores, so the control limits calculated for T2 are also compromised. To enhance the ability of dynamic process monitoring, researchers developed time-series models on scores.15-17 Because the scores extract only a portion of the dynamics, the monitoring efficiency and accuracy are still doubtful. Ku et al., in their article, declared that “building a dynamic model only for the score values of the principal components obtained from an initial static PCA step will not solve the problem because the linear constraints that span the noise space will remain static”.18 Dynamic PCA18 can extract all the dynamic information with scores, and only leaves noise behind. This method improves the detection results, based on SPE. However, the dynamics that are retained in scores break the statistical assumption of T2 control limits calculation. Autoregressive moving average (ARMA) filters are designed to solve the problem.19 The dynamics in scores are filtered, so the retained portions are independent of time. The T2 statistic, and its control limits that are calculated based on the retained values of scores, then are more reasonable. The number of false alarms can be reduced. A potential problem with this method is that the univariate ARMA filters can only remove the autocorrelations in the scores, not the cross-correlations. The cross-correlation function plots in the work of Ku et al.18 clearly show that there still exist cross-correlations among scores at different steps. To counter this problem, subspace identification-based PCA (SI-PCA) has been proposed.20 In that work, a state space model instead of ARMA models is used to filter the dynamics in scores. After doing so, most autocorrelations and cross-correlations are removed from the scores. 3. Review of the 2-D-DPCA Method Dynamics are important inherent characteristics of batch processes. In some cases, such dynamics exist not only within a particular batch, but also from batch to batch. Several different reasons may cause batchwise dynamics, such as slow property changing of feed stocks, drift of the process characteristics, process controllers designed for run-to-run adjustment, the effects of slow response variables, and so on. All these reasons are common in batch processes. To monitor the 2-D dynamic batch processes better, a parsimonious 2-D model structure is needed. When 2-D dynamics exist, the current measurements are dependent not only on lagged measurements in the past time direction in the
same batch, but also on lagged measurements in some past batches. All these lagged measurements that can reflect 2-D dynamics locate in a region. We call this region the support region or the region of support (ROS). 2-D-DPCA expands the data matrix to include all the lagged measurements in ROS, together with current measurements, and then performs PCA on such a matrix. Thus, both 2-D dynamics and cross-correlation information among variables can be extracted simultaneously. Only noise remains in the residual space. Therefore, SPE can be calculated to perform process monitoring. To build a good 2-D-DPCA model, a key step is to choose the ROS reasonably. Based on the previous discussion, a proper ROS should include the necessary past data to provide good prediction of the variables’ current values. In a sense, ROS determination is a problem similar to the task of selecting proper independent variables to predict the current sample in the development of a prediction model. Enlightened by such idea, a data-based ROS selection method is designed. First, the region including candidate independent variables for current measurements prediction is selected based on the values of autocorrelation and cross-correlation functions. This candidate region is also called the initial ROS. The target proper ROS then can be obtained by eliminating unnecessary independent variables in the initial ROS. Iterative backward elimination procedures are designed to accomplish this task. For more details about the 2-D-DPCA method with autoselected ROS, please refer to our previous work.12,13 4. 2-D Autocorrelations in the Score Space of the 2-D-DPCA Model As noted by Kruger et al., “If a dynamic AR model structure is incorporated into the conventional PCA analysis, the retained score variables are auto-correlated irrespective of whether the process variables are auto-correlated or not.”19 Similar situations can be observed in 2-D dynamic cases. Because the 2-D-DPCA model utilizes a 2-D AR structure, the retained scores are 2-D autocorrelated and cross-correlated in a complicated way. The score space of 2-D-DPCA extracts both cross-correlations and 2-D dynamic information, so the score values at different sampling intervals are not independent. There are significant autocorrelations in both the time direction in a particular batch and the batch direction at each sampling time interval, and even in the cross-directions. In addition, the cross-correlations also exist among scores of different steps. In some situations, the dynamics structure in score space may be even more complicated than that in the original data space. Such features obviously break the statistical basis of the control limits calculation for T2. Therefore, using a T2 control plot directly for monitoring is not proper. Some pretreatments should be made before T2 and corresponding control limit calculations. The complexity of the dynamics structure in score space can be shown with the simulation example in section 7. 5. 2-D Multivariate AR Score Filters for 2-D-DPCA 5.1. Removing of Dynamics with an AR Filter. To make a good control plot for process monitoring, reasonable control limits should be set, which means the calculation of control limits should obey the corresponding statistical basis. Therefore, the 2-D dynamics in scores including autocorrelation and crosscorrelations should be removed as much as possible before calculation of the T2 parameter. In this paper, 2-D multivariate AR score filters are designed to eliminate the 2-D dynamics.
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8035
Figure 1. Schematic diagrams of the initial ROS and proper ROS.
The common form of a one-dimensional (1-D) univariate AR model is similar to
xt ) a1xt-1 + ... + apxt-p + zt
(3)
where zt is a purely random process with E(zt) ) 0, Var(zt) ) σz2, p is the order of the model, and the parameters a1, ..., ap are called AR coefficients. For a 1-D time-series process, the 1-D AR filter can be designed in such a way that
xˆ t ) a1xt-1 + ... + apxt-p
(4)
The filtered signals then can be calculated as
xft ) xt - xˆ t ) zt
(5)
Obviously, the dynamics are removed by the AR filter, and the filtered signals are time-independent. 5.2. Requirements of Filter Design and Related Works. The principle and purpose of using a 2-D multivariate AR scores filter is similar to those of the 1-D univariate AR filter. However, in our cases, there are several special requirements and difficulties in designing such filters. (1) Two dimensions are required. The dynamics in both the time and batch directions should be filtered simultaneously. (2) Multiple variables are required. Univariate filters may be sufficient to examine the autocorrelations, but they cannot effectively filter out cross-correlations at different steps. (3) The ROS must have an irregular shape. An AR model can be regarded as a regression model, which relates variable past values with current values. Thus, an important step of building a good AR model is to choose the variable values at past steps that can provide a good prediction of the variables’ current values. In the design of a 1-D AR model, this is called order selection. In building the 2-D AR model, a corresponding concept is the selection of the ROS. Similar to the ROS of the 2-D-DPCA model, the ROS of the 2-D AR filter is a region that contains the necessary lagged variables to provide good prediction of the variables’ current values. If a 2-D AR model is to be designed, order pairs should be selected if the ROS is assumed to be located in a quarter plane or a half plane with a regular shape (e.g., rectangular shape). However, the data of the 2-D dynamic batch process may not have a regular ROS, because of the complicated and various process characteristics. Moreover, as illustrated in section 4, the 2-D dynamics retained in scores can be more complex than
Figure 2. Procedural flowchart for the entire 2-D-DPCA-based batch process modeling and monitoring.
those in original process data. Thus, the shapes of the ROS of the 2-D multivariate AR score filters are even harder to know. Many works have discussed how to select the order pairs for a 2-D AR or 2-D ARMA model.21-23 However, these methods all assume that the shape of the ROS is known and are designed for univariate cases. To the authors’ knowledge, no order selection method for the 2-D multivariate AR model is observed in the literature. In 1-D cases, the multivariate time-series modeling methods have been studied,24 but no further extension to 2-D cases is found. 5.3. 2-D multivariate AR Score Filters with Auto-selected ROS. Usually, a 2-D AR filter is described as p1-1 p2-1
yˆ (i,k) )
∑ ∑ a(m, n)y(i - m,k - n)
(6)
m)0 n)0
where y(i,k) is the value of variable y with index pair (i,k), a(m,n) is the corresponding coefficient of y(i - m,k - n), yˆ (i,k) is the estimation of y(i,k), p1 and p2 are the orders of both directions (which determine the ROS), and (m,n) * (0,0). Obviously, this type of form cannot meet all three of the requirements described previously. Although it can describe the autocorrelations in both the time and batch directions, it still has several shortcomings. First, the estimation of a current
8036
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
Figure 3. Cross-correlation function plots of scores in the time direction.
variable is only based on the past value of the same variable, so it is a univariate form and no cross-correlations can be included. Second, the order pair (p1,p2) can only determine a ROS with a rectangular shape in the past quarter-plane but not an irregular shape in the past half-plane. To overcome these problems, a 2-D multivariate AR score filter for the jth score is formulated as J
ˆt j(i,k) )
qd1(0)
∑ (∑
d1)1 d3)1 pd1
ad1(0,d3)td1(i,k - d3) + qd1(d2)
∑ ∑
d2)1 d3)-fd1(d2)
ad1(d2, d3)td1(i - d2,k - d3)) (7)
where tj(i,k) means the value of score j at the kth sampling interval in batch i, ad1(d2,d3) is the AR coefficient that corresponds to td1(i - d2,k - d3), pd1 is the maximum number of lagged batches of score d1 included in the AR filter, qd1(d2) is the maximum number of past sampling intervals of score d1 in the (i - d2)th batch included, fd1(d2) is the maximum number of future sampling intervals of score d1 in the (i - d2)th batch included, and J is the maximum number of scores. The filtered scores then can be calculated as
tjf(i,k) ) tj(i,k) - ˆt j(i,k) which are time-independent.
(8)
From eq 7, we can observe that this AR model relates the current value of the jth score with the past value of all score variables in some regions (i.e., ROS) in the past half-plane. Because the ROS include past values of score variables in both the time and batch directions, if it is chosen reasonably, 2-D dynamic information can be extracted well. At the same time, because of the inclusion of past values of all score variables, not only autocorrelations but also cross-correlations at different steps are taken into consideration. And it can be further noticed that the ROS corresponding to different score variables can be different, and each ROS is not required to have a regular shape. Thus, the structure of this AR model ensures the satisfaction of the three requirements described in the last section. The problem is transformed to “how to select the ROS reasonably based on historical data”. After the ROS is determined properly, the AR coefficients can be calculated by regression. The techniques used in the task of ROS selection for 2-D multivariate AR score filters is similar to those used in ROS selection for the 2-D-DPCA model. The major difference is the application. The former one determines the ROS for a 2-DDPCA model. Based on that technique, a 2-D-DPCA model can be built for a batch process. Two-dimensional process dynamics are extracted by the score space and noise is left in the residual space. The latter one selects the ROS for 2-D multivariate AR score filters in the score space of the 2-D-DPCA model. The 2-D dynamics in the scores then can be filtered by the score filters that are designed based on this ROS.
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8037
Figure 4. Cross-correlation function plots of scores in the batch direction.
There are two components of work in the ROS autoselection for 2-D multivariate AR score filters. One component is called the initial ROS selection, and the other component is called the proper ROS determination. The purpose of the initial ROS selection is to provide a candidate region of past variables that is large enough to contain the entire proper ROS. At the same time, it should not be too large, so that the computation burden can be controlled. Without process knowledge, several variable selection methods can be used for the initial ROS selection, such as the simple regression (SR) method, the random correlations (RCOR) method, and the coefficients (COEF) method.25 We do the initial ROS selection based on the correlations between the current sample and the lagged variables, which can be indicated by the values of autocorrelation and cross-correlation functions, because the lagged variables that significantly correlate with the current sample may provide good prediction. When there is a significant correlation, the value of the function is ∼1; otherwise, it approaches zero. Thus, the past region that contains significant correlations to the current score values is selected as the initial ROS. This method is similar to the SR method, in regard to variable selection. The proper ROS then can be determined from this initial ROS by backward elimination. In each run of the iterative backward eliminations, a regression model is built to relate the current sample with lagged samples in the candidate region, which are
regarded as the independent variables, and one independent variable that corresponds to the smallest regression coefficient is eliminated from the candidate region. Because the coefficient of each independent variable denotes its contribution to the dependent variable, if all the variables are normalized, the independent variable that is eliminated in each run is the one that contributes the least to the prediction. This iteration continues until only one variable is retained. At the same time, an index is used as a criterion to evaluate the regression model in each run. This index could be the Akaike Information Criterion (AIC)26 or any of the other indices that are used in variable selection. Based on the AIC criterion, the smaller the index, the better the prediction. The best choice of the ROS is determined based on a comparison of such index values calculated during the iteration. Details of the method are given below. (1) Suppose there are a total of J score vectors retained in a 2-D-DPCA model. Select an initial ROS as the candidate region in the past half-plane for further selection. In the example shown in Figure 1, the initial ROS is selected as tj(i,k - 1), tj(i,k 2), ..., tj(i,k - nj(0)), tj(i - 1,k + rj(1)), tj(i - 1,k + rj(1) 1), ..., tj(i - 1,k), tj(i - 1,k - 1), ..., tj(i - 1,k - nj(1)), ..., tj(i - mj,k + rj(mj)), ..., tj(i - mj,k), tj(i - mj,k - 1), ..., tj(i - mj,k - nj(mj)), where tj(i,k) represents the value of score j at the kth
8038
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
Figure 5. Training data of 2-D multivariate AR score filter for the first score: (s9s) training data and (‚ ‚ ‚f‚ ‚ ‚) prediction.
Figure 8. Checking data of 2-D multivariate AR score filter for the second score: (s9s) checking data and (‚ ‚ ‚f‚ ‚ ‚) prediction.
[]
(2) Augment all lagged scores in the initial ROS into a data matrix as
Tm+1,n+1 l Tm+1,K-r X) l Ti,k l TI,K-r where
(9)
Ti,k ) [t1(i,vk), ..., tj(i,k), ..., tJ(i,k)] Figure 6. Training data of 2-D multivariate AR score filter for the second score: (s9s) training data and (‚ ‚ ‚f‚ ‚ ‚) prediction.
tj(i,k) ) [tj(i,k - 1), ..., tj(i,k - nj(0)), tj(i - 1,k + rj(1)), tj(i - 1,k + rj(1) - 1), ..., tj(i - 1,k), tj(i - 1,k - 1), ..., tj(i - 1,k - nj(1)), ..., tj(i - mj,k + rj(mj)), tj(i - mj, k + rj(mj) - 1), ..., tj(i - mj,k), tj(i - mj,k - 1), ..., tj(i - mj,k - nj(mj))] n ) max(n1(0), n1(1), ..., n1(m), ..., nj(0), nj(1), ..., nj(m), ..., nJ(0), nJ(1), ..., nJ(m)) r ) max(r1(1), r1(2), ..., r1(m), ..., rj(1), rj(2), ..., rj(m), ..., rj(1), rj(2), rj(m)) m ) max(m1,...,mj,...,mJ) (3) Normalize X into zero means and unit variances. (4) Choose a score th (where h ) 1, 2, ..., J) and J is the number of process variables and augment data into a matrix, which is given as
Figure 7. Checking data of 2-D multivariate AR score filter for the first score: (s9s) checking data and (‚ ‚ ‚f‚ ‚ ‚) prediction.
sampling interval in batch i, mj is the maximum number of lagged batches of score j includes in the initial ROS, nj(V) is the maximum number of past sampling intervals of score j in the (i - V)th batch included therein, rj(V) is the maximum number of future sampling intervals of score j in the (i - V)th batch in this initial ROS, and j ) 1, 2, ..., J.
[ ]
th(m + 1,n + 1) l th(m + 1,K - r) Yh ) l th(i,k) l th(I,K - r)
(5) Normalize Yh to a zero mean and unit variance.
(10)
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8039
Figure 9. Cross-correlation function plots of filtered scores in the time direction.
(6) Take X as the matrix of independent variables, Yh as the matrix of dependent variable, and build a regression model between X and Yh. (7) Calculate and store the evaluation index value, to assess this regression model. The index can be AIC or other indices that are used for model evaluation. (8) Eliminate a column in X that corresponds to the smallest coefficient in the regression model, because the smallest coefficient indicates the smallest contribution to the prediction. (9) Return to step 6 until all columns in X have been eliminated. (10) Compare all of the evaluation index values calculated and determine the one that indicates the most proper model. If AIC indices are used, to overcome the over-fitting problem, the one closer to the minimum value with fewer independent variables are retained. (11) Check X that corresponds to the best model determined in step 10. The remaining columns of X give the most proper ROS of th. The selected proper ROS is also shown in Figure 1. Similar to initial ROS, the proper ROS may have irregular shape in the past half-plane, which includes tj(i,k - 1), tj(i,k - 2), ..., tj(i,k - qj(0)), tj(i - 1,k + fj(1)), tj(i - 1,k + fj(1) - 1), ..., tj(i - 1,k), tj(i - 1,k - 1), tj(i - 1,k - qj(1)), ..., tj(i - pj,k + fj(pj)), ..., tj(i - pj,k), tj(i - pj,k - 1), tj(i - pj,k - qj(pj)), where pj is the maximum number of lagged batches of score j included in the proper ROS of score h, qj(V) is the maximum number of past sampling intervals of score j in the (i - V)th batch in this
ROS, fj(V) is the maximum number of future sampling intervals of score j in the (i - V)th batch included therein, and j ) 1, 2, ..., J. At the same time, the regression model that related this X to Yh is the 2-D multivariate AR filter for score th. (12) Return to step 4. Choose another th and find its proper ROS and corresponding 2-D multivariate AR filter. (13) If all 2-D multivariate AR score filters are designed, the procedures are ended. (14) Regress the lagged variables that have been selected in the ROS with the current variables to determine the coefficients of the 2-D multivariate AR score filters. 5.4. Application of 2-D Multivariate AR Score Filters for Better Monitoring. Because these filters catch most of the 2-D autocorrelations and cross-correlations, they are used to remove the dynamics from the scores, as shown in eq 8. After doing so, the filtered scores mainly consist of normal distributed statistical information. This means that we can use the filtered scores to calculate the T2 statistic and the corresponding control limits can be obtained based on the common statistical assumption. Monitoring based on this new T2 value is more reasonable than monitoring based on the T2 value of unfiltered scores. 6. Batch Process Monitoring Based on 2-D-DPCA Figure 2 shows the entire procedure of 2-D dynamic batch process monitoring based on 2-D-DPCA, from 2-D-DPCA model building (including the ROS selection) to online fault
8040
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
Figure 10. Cross-correlation function plots of filtered scores in the batch direction.
diagnosis in both residual and score spaces. First, the ROS selection for the 2-D-DPCA model is performed based on the normal history data. The 2-D-DPCA model then is built with such ROS and history data, and the original data space is divided into two subspaces. In residual space, the SPE control limits are calculated. In score space, the 2-D multivariate AR score filters are designed to remove the dynamics that have been retained. This design step includes the ROS selection for the filters, and coefficients calculation, which have been introduced in section 5.3. The T2 control limits then are calculated based on the filtered score values. When new process data are obtained, the SPE can be calculated with the 2-D-DPCA model. By comparing the values of SPE to its control limits, if there are any faults that are related to correlation structure changes, they can be detected. At the same time, in score space, the new score values are filtered with the score filters and used in calculation of the T2 values. These new T2 values are compared with the control limits that had been calculated previously. If the faults are with the variable magnitude changes, the T2 values will be beyond the control limits, so that the faults can be detected. After a fault is detected, fault diagnosis can be performed based on contribution plots.27 7. Simulation Results In this part, simulations are performed to show the benefit of the T2 control plot with 2-D multivariate AR score filters by
comparing the monitoring efficiency values. A batch process is described in eq 11, where a batch process is described with the formula
x1(i,k) ) 0.8x1(i - 1,k) + 0.5x1(i,k - 1) 0.33x1(i - 1,k - 1) + w1 (11a) x2(i,k) ) 0.44x2(i - 1,k) + 0.67x2(i,k - 1) 0.11x2(i - 1,k - 1) + w2 (11b) x3(i,k) ) 0.65x1(i,k) + 0.35x2(i,k) + w3
(11c)
x4(i,k) ) - 1.26x1(i,k) + 0.33x2(i,k) + w4
(11d)
where x1 and x2 are two independent signals with dynamics in the time direction; x3and x4 are functions of x1, x2, and its own value at one step before; wj (where j ) 1, 2, 3, 4) are Gaussian noises with a variance of 0.01. One hundred data batches are generated, and there are 200 samples in each batch. A 2-D-DPCA model is used to model the process. The ROS of such model is selected as
[x(i,k), x(i,k - 1), x(i - 1,k), x(i - 1,k - 1)] where
x(i,k) ) [x1(i,k), x2(i,k), x3(i,k), x4(i,k)]
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8041
This is consistent with the process. Four of the 16 score variables are chosen with a cross-validation algorithm to be retained in the 2-D-DPCA model. After performing 2-D-DPCA, only noise information is retained in the residuals, and all dynamics are extracted by the score space. 7.1. Comparison of Dynamics in Unfiltered and Filtered Scores. Now let us examine the cross-correlation function plots of the retained scores. Figure 3 shows the autocorrelations and cross-correlations along the time direction in the same batches. In addition, the autocorrelations and cross-correlations along thebatch direction at the same time intervals are plotted in Figure 4. It can be observed that the cross-correlation functions among scores at the same step are all equal to zero in both figures. This is due to the property of the PCA algorithm. However, there are significant autocorrelations in both the time and batch directions, as shown by the plots on the diagonal of each figure. The cross-correlations also exist among scores of different steps, especially between the first score vector and other scores. Although the structure of the 2-D dynamics of this process is quite simple, the 2-D autocorrelations and cross-correlations that have been retained in score space are very complex. Such features obviously break the statistical basis of the control limits calculation for T2. Thus, using the T2 control plot directly for monitoring is not proper. This can affect the fault detection efficiency, based on the T2 control plot. For comparison, our proposed 2-D multivariate AR score filters are designed for the process. The results of the ROS selection for the filters indicate that, although the original process data have a simple quarter-plane support region, the ROS of the retained scores could be much more complex. Taking the first score t1 as an example, its determined ROS for the 2-D AR filter is
[t1(i,k - 1), t1(i - 1,k + 1), t1(i - 1,k), t1(i - 1,k - 1), t1(i - 2,k), t3(i,k - 1), t3(i - 1,k + 1), t4(i,k - 1), t4(i - 1,k)] which is an irregular region in the past half-plane. This shows the complicated structure of the 2-D dynamics in the retained scores. After the ROS of each score is selected, the coefficients of the 2-D multivariate AR score filters are calculated. The filters then are applied in the scores filtering. Figures 5 and 6 show the training errors of the 2-D multivariate score AR filters of the first two scores, and the checking errors are shown in Figures 7 and 8. Figure 9 shows the autocorrelations and cross-correlations along the time direction in the same batches, and Figure 10 shows the autocorrelations and cross-correlations along the batch direction at the same time intervals. It is observed that, by comparing to Figure 3 and Figure 4, the autocorrelations and cross-correlations in score space are largely removed after the filtering. Therefore, the filtered scores are more suitable to be used in T2 control plot development. 7.2. Comparison of Monitoring Efficiencies in Score Space. To compare the fault detection efficiencies of the T2 control plots based on unfiltered and filtered scores, two different types of faults are generated to be detected. This example is the same as that in the previous work.12 However, this time, it is used to check the monitoring ability of the methods in score space. The first fault is caused by changing the variable correlation structure. From batch 61, process variable x2 is formulated as
x2(i, k) ) 0.67x2(i - 1,k) + 0.8x2(i,k - 1) 0.47x2(i - 1,k - 1) + w2 (12) As shown in previous work,12 although the autocorrelation is changed quite significantly, the changes in the magnitude of the variables’ values are not obvious until five batches later. This makes the fault detection somewhat difficult. The other fault is a type of small process drift on variable x2 from batch 61 by adding a signal that increases slowly with time and batch. The T2 control plot that is based on unfiltered scores then is used in online monitoring. Because one step of lagged measurements in the time direction is needed in modeling, a total of 199 samples are monitored in a batch. In Figure 11, we can observe that such a T2 control plot cannot clearly detect the fault of changing the correlation structure in batch 61 at all. In Figure 12, there are 10 samples (∼5.03%) beyond the 99% control limit and 49 samples (∼24.62%) beyond the 95% control limit in batch 61. Therefore, for this drift fault, the control plot can detect it in batch 61, but one also notices that the trend of gradual increase is not shown clearly. The online monitoring results of the T2 control plot of the filtered scores are displayed in Figures 13 and 14. Because there are a total of 5 steps of lagged variables in the time direction used in modeling, the monitoring results of 195 samples in each batch are plotted. No false alarms are observed in the monitoring of the normal batch. For the detection of a change in the variable correlation structure, 18 of the 195 samples (∼9.23%) are outside the 99% control limit and 39 of the 195 samples (∼20%) are outside the 95% control limit in batch 61. For the detection of process drift, the numbers of samples outside the 99% and 95% control limits are 67 (∼34.36%) and 79 (∼40.51%), respectively. This is a clearer result than that based on unfiltered scores. In addition, the increasing trend is also more obvious than previously observed. Thus, we can observe that the filtered scores give clear and efficient detection results. The result in Figure 13 is especially much better than the original one (Figure 11). 7.3. Comparison between SPE and T2, Based on Filtered Scores. As stated in section 2.1, the SPE and T2 statistics are concerned with two different types of information. When there is a fault in the magnitude of a variable’s value that does not affect the correlation structures, the T2 statistic provides morereliable fault detection results. In this section, an example is used to illustrate this observation. The process is still that same as that described in eq 11. From the start of batch 61, the value of variable x2 drifts away from its normal value significantly, which causes a change in the magnitude of x2. Such drift stops at the end of batch 61, which means, from batch 62, the process correlation structure has been recovered to normal, but the magnitudes of the values of variable x2 are still abnormal. The SPE control plot detects this drift in batch 61, because the drift of variable x2 changes the correlation structure of the process variables. However, when the variable correlation structure becomes stable, the SPE control plot cannot reflect the change in the magnitude of variable x2’s value and causes a missed alarm. All these scenarios are plotted in Figure 15. This may give the operator an incorrect signal that the system has recovered from the fault automatically and the operation has become normal already, which is false. Now, let us examine the monitoring result of the T2 control plot based on the filtered scores shown in Figure 16. It detects the drift fault efficiently in batch 61, because the drift affects the dynamic features of the process. Even in the following batches, the T2 values are significantly beyond the control limits.
8042
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007
Figure 11. Monitoring of the changing correlation structure fault with the T2 control plot, based on unfiltered scores: (- - -) 95% control limit and (s) 99% control limit.
Figure 14. Monitoring of the process drift fault with the T2 control plot, based on filtered scores: (- - -) 95% control limit and (s) 99% control limit.
Figure 12. Monitoring of the process drift fault with the T2 control plot, based on unfiltered scores: (- - -) 95% control limit and (s) 99% control limit.
Figure 15. Monitoring of the change in magnitude of the process variable’s value with the SPE control plot: (‚ ‚ ‚) 95% control limit and (s) 99% control limit.
Figure 13. Monitoring of the changing correlation structure fault with the T2 control plot, based on filtered scores: (- - -) 95% control limit and (s) 99% control limit.
Figure 16. Monitoring of the change in magnitude of the process variable’s value with the T2 control plot, based on filtered scores: (‚ ‚ ‚) 95% control limit and (s) 99% control limit.
This example shows that the T2 control plot, based on filtered scores, is a good complement of the SPE statistic, and it is more effective in the detection of the faults in the magnitudes of the variables’ values than the SPE control plot.
processes with two-dimensional (2-D) dynamics. However, the dynamic information that is retained in the score space violate the statistical basis of T2 control plot monitoring. Thus, only information in the residual space is used in online fault detection. In this work, 2-D multivariate autoregressive (AR) score filters are designed to remove the autocorrelations and cross-correlations in score space. The region of support (ROS) auto-selection method is used to determine the ROS of these filters. The
8. Conclusion In previous work, two-dimensional dynamic principal component analysis (2-D-DPCA) was proposed to monitor the batch
Ind. Eng. Chem. Res., Vol. 46, No. 24, 2007 8043
coefficients of the filters then can be obtained using the regression method. After filtering, most of the dynamic information is removed, and the monitoring efficiency, based on the information in the score space, is improved significantly. Literature Cited (1) Wold, S.; Kettaneh, N.; Friden, H.; Holmberg, A. Modeling and diagnosis of batch processes and analogous kinetic experiments. Chemom. Intell. Lab. Syst. 1998, 44, 331-340. (2) Sprange, E. N. M.; Ramaker, H. J.; Westerhuis, J.; Smilde, A. Critical evaluation of approaches for on-line batch process monitoring. Chem. Eng. Sci. 2002, 57, 3979-3991. (3) U ¨ ndey, C.; Cinar, A. Statistical monitoring of multistage, multiphase batch processes. IEEE Control Syst. Mag. 2002, 22, 40-52. (4) Kourti, T. Multivariate dynamic data modeling for analysis and statistical process control of batch processes, start-ups and grade transitions. J. Chemom. 2003, 17, 93-109. (5) Nomikos, P.; MacGregor, J. F. Monitoring of batch processes using multiway principal component analysis. AIChE J. 1994, 40, 1361-1375. (6) Nomikos, P.; MacGregor, J. F. Multiway partial least squares in monitoring batch processes. Chemom. Intell. Lab. Syst. 1995, 30, 97-108. (7) Wold, S.; Kettaneh, N.; Tjessem, K. Hierarchical multiblock PLS and PC models for easier model interpretation and as an alternative to variable selection. J. Chemom. 1996, 10, 463-482. (8) Ra¨nnar, S.; MacGregor, J. F.; Wold, S. Adaptive batch monitoring suing hierarchical PCA. Chemom. Intell. Lab. Syst. 1998, 41, 73-81. (9) Lu, N.; Gao, F.; Wang, F. Sub PCA modeling and monitoring strategy for batch processes. AIChE J. 2004, 50, 255-259. (10) Lee, J. H.; Dorsey, A. W. Monitoring of batch processes through state-space models. AIChE J. 2004, 50, 1198-1210. (11) Flores-Cerrillo, J.; MacGregor, J. F. Multivariate monitoring of batch processes using batch-to-batch information. AIChE J. 2004, 50, 12191228. (12) Lu, N.; Yao, Y.; Gao, F. Two-dimensional dynamic PCA for batch process monitoring. AIChE J. 2005, 51, 3300-3304. (13) Yao, Y.; Lu, N.; Gao, F. Two-Dimensional Dynamic PCA with Auto-Selected Support Region. In DYCOPS 2007: Proceedings of the 8th International IFAC Symposium on Dynamics and Control of Process Systems, Cancu´n, Mexico, June 6-8, 2007; International Federation of Automated Control: Preprints, Vol. 1, pp 69-74. (14) Jackson, J. E. A User’s Guide to Principal Components; Wiley: New York, 2001. (15) Wikstro¨m, C.; Albano, C.; Eriksson, L.; Friden, H.; Johansson, E.; Nordahl, A.; Rannar, S.; Sandberg, M.; Kettaneh-Wold, N.; Wold, S.
Multivariate process and quality monitoring applied to an electrolysis process. Part I. Process supervision with multivariate control charts. Chemom. Intell. Lab. Syst. 1998, 42, 221-231. (16) Wikstro¨m, C.; Albano, C.; Eriksson, L.; Friden, H.; Johansson, E.; Nordahl, A.; Rannar, S.; Sandberg, M.; Kettaneh-Wold, N.; Wold, S. Multivariate process and quality monitoring applied to an electrolysis process. Part II. Multivariate time-series analysis of lagged latent variables. Chemom. Intell. Lab. Syst. 1998, 42, 233-240. (17) Samanta, B. Multivariate control charts for grade control using principal-component analysis and time-series modeling. Trans. Inst. Min. Metall., Sect. A 2002, 111, 149-157. (18) Ku, W.; Storer, R. H.; Georgakis, C. Disturbance detection and isolation by dynamic principal component analysis. Chemom. Intell. Lab. Syst. 1995, 30, 179-196. (19) Kruger, U.; Zhou, Y.; Irwin, G. W. Improved principal component monitoring of large-scale processes. J. Process Control 2004, 14, 879888. (20) Xie, L.; Zhang, J.; Wang, S. Investigation of dynamic multivariate chemical process monitoring. Chin. J. Chem. Eng. 2006, 14, 559-568. (21) Miyazaki, N.; Hamada, N. An identification method for a twodimensional nonminimum phase ARMA system using a phase-equivalent MA system. Electron. Commun. Jpn., Part 3 1997, 80, 1454-1463. (22) Aksasse, B.; Radouane, L. Two-dimensional autoregressive (2-DAR) model order estimation. IEEE Trans. Signal Process. 1999, 47, 20722077. (23) Ojeda, S.; Vallejos, R.; Lucini, M. Performance of robust RA estimator for bidimensional autoregressive models. J. Stat. Comput. Simul. 2002, 72, 47-62. (24) Pen˜a, D.; Tiao, G. C.; Tsay, R. S. A Course in Time Series Analysis; Wiley: New York, 2000. (25) Gauchi, J.; Chagnon, P. Comparison of selection methods of explanatory variables in PLS regression with application to manufacturing process data. Chemom. Intell. Lab. Syst. 2001, 58, 179-193. (26) Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control 1974, 19, 716-723. (27) Westerhuis, J. A.; Gurden, S. P.; Smilde, A. K. Generalized contribution plots in multivariate statistical process monitoring. Chemom. Intell. Lab. Syst. 2000, 51, 95-114.
ReceiVed for reView April 25, 2007 ReVised manuscript receiVed August 30, 2007 Accepted August 31, 2007 IE070579A
