Statistical Process Control Charts for Batch Operations Based on

Ind. Eng. Chem. Res. 2004, 43, 6731-6741

6731

PROCESS DESIGN AND CONTROL Statistical Process Control Charts for Batch Operations Based on Independent Component Analysis Hamza Albazzaz and Xue Z. Wang* Department of Chemical Engineering, The University of Leeds, Leeds LS2 9JT, United Kingdom

A significant step forward in recent years, in regard to multivariate statistical process control (MSPC) for operational condition monitoring and fault diagnosis, has been the introduction of principal component analysis (PCA) for the compression of process data. An alternative technique that has been studied more recently for data compression is independent component analysis (ICA). Published work has shown that, in some applications of statistical process monitoring, ICA-based methods have exhibited advantages over those based on other data compression techniques. However, it is inappropriate to use ICA in the same way as PCA to derive Hotelling’s T2 and SPE (squared prediction error) charts, because the independent components are separated by maximizing their non-Gaussianity, whereas the satisfying Gaussian distribution is the basis of T2 and SPE monitoring charts, as well as univariate statistical process control (SPC) charts. In this paper, we propose a new method for deriving SPC charts based on ICA, which can overcome the aforementioned limitation of non-Gaussianity of the independent components (ICs). The method generates a smaller number of variables, i.e., ICs to monitor, each with time-varying upper and lower control SPC limits, and, therefore, can be used to monitor the evolution of a batch run from one time point to another. The method is illustrated in detail by reference to a simulated semibatch polymerization reactor. To test its capability for generalization, it is also applied to a data set that has been collected from industry and proved to be able to detect all seven faults in a straightforward way. A third case study that was studied in the literature for batch statistical monitoring is used in this work, to compare the performance of the current approach with that of other methods. It proves that the new approach can detect the faults earlier than a similar PCA-based method, the PCA-based T2 approach, and the SPE approach. Comparison with a recently proposed multi-way ICA method in the literature was also made. 1. Introduction Tracking the batch-to-batch variation and detecting abnormal events at early stages of a batch run is critically important in the pharmaceutical, materials, and food industries, as well as many other industries that utilize batchwise operations. Multivariate statistical process control (MSPC) that is based on Hotelling’s T2 and squared prediction error (SPE) charts has been a key technique for this purpose, and significant progress has been made since Kresta et al.1 introduced the idea of using principal component analysis (PCA) to preprocess the process variables. PCA provides a smaller number of latent variables that can be used to replace the original observed variables in calculating the T2 and SPE charts. Each principal component (PC) is a linear combination of the original variables; the first PC captures the majority of variation in the data, the second PC captures the majority of the remaining variation and is orthogonal to the first PC, and so forth. There has been a wealth of publications in this area over the last 10 years, and the work has been reviewed by several researchers.2-7 More recently, an alternative approach to PCA has attracted attention for data compression; this method * To whom correspondence should be addressed. Tel: +44 113 343 2427. Fax: +44 113 343 2405. E-mail: che6xzw@ leeds.ac.uk.

is based on independent component analysis (ICA).8-17 In contrast to PCs, which are linear combinations of the observed variables, the independent components (ICs) are independent variables or source events that constitute (through linear or nonlinear combinations) the observed variables; ICA is the process of blind separation of these source or independent variables. Published work has demonstrated that, in some statistical process monitoring applications, techniques based on ICA have advantages over those based on PCA.8,9,11,14 However, a critical observation that we would like to make here is that the ICs do not satisfy a Gaussian distribution, because non-Gaussianity is the basis of ICA operation, whereas the T2 and SPE control limits are set based on the assumption that the variables, either original or latent, follow a Gaussian distribution. As a result, ICs cannot be used in the same way as PCs to derive T2 and SPE monitoring charts, as well as univariate SPC charts. In this paper, we propose a new SPC method for batch processes that can overcome this limitation of ICA. The paper is organized as follows. In the next section, the well-known PCA-based MSPC method will be described, but only very briefly, and the focus will be on introducing the basics, instead of reviewing the literature. The ICA algorithm will be introduced in Section 3, and published work on applying ICA to process

10.1021/ie049582+ CCC: $27.50 © 2004 American Chemical Society Published on Web 09/15/2004

6732

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004

monitoring and control will be reviewed in detail in the same section. The proposed approach for batch process monitoring will be presented in Section 4. Application of the proposed method to three case studies is described in Section 5, which is followed by final remarks in Section 6. 2. Multivariate Statistical Process Control Based on Principal Component Analysis Traditionally, univariate SPC charts (e.g., the Shewhart chart) were used in industry to separately monitor either a few process variables or key measurements on the final product that, in some way, define the quality of the product. For common cause variation (e.g., random disturbances), the values of a variable (typically, a measure of product quality) should satisfy a Gaussian distribution centered around the mean µ with a standard deviation σ. If the upper and lower control limits (UCL and LCL, respectively) are set as +σ and -σ, then the probability that the value of the variable goes beyond UCL and LCL is 5% is an indication of the potential occurrence of special cause variation, or faults. Plant supervisors and operators should then examine the process to determine the origin of the special cause. The method can be extended to simultaneously consider multiple variables for designing multivariate Shewhart charts, which are called Hotelling’s T2 charts. T2 charts are similar to the univariate Shewhart charts and a 95% control limit can also be calculated, but these charts have taken the multivariate variables into consideration. Considering the fact that, although a process may have several tens to several hundreds of variables being monitored simultaneously, they are often correlated. A significant development in MSPC in recent years has been the introduction of PCA for compression of the historical operational data prior to application of the two MSPC statistical models: Hotelling’s T2 and SPE models.1,18,19 PCA is a technique that reduces the dimensionality of a data set that consists of a large number of inter-related variables, while retaining as much of the variation present in the data set as possible. The PCA approach uses all of the original variables to obtain a smaller set of variables (which are called PCs) that can be used to approximate the original variables. The PCs are uncorrelated and are ordered so that the first few retain most of the variation present in the original data set. PCA can be applied to process the data before the multivariate Shewhart charts are used. This means that the latent variables (i.e., the first few PCs) are used rather than the original variables. The Hotelling’s T2 value can then be calculated, by A

2

T )

∑ i)1

ti2

sti2

(1)

where sti2 is the estimated variance of ti, ti is the ith PC, and A is the number of PCs selected. However, monitoring the product quality via a T2 value that is based on the first A number of PCs is not sufficient. If a totally new type of special event occurs that was not present in the reference data used to develop the in-control PCA model, a new statistic must

be used, i.e., the squared prediction error (SPE) of the residual of a new observation: k

SPE )

(ynew,i - yˆ new,i)2 ∑ i)1

(2)

where yˆ new is computed from the reference PCA model. SPE is also called the Q-statistic or distance to the model. It represents the squared perpendicular distance of a new multivariate observation from the projection space. The historical operational data of continuous processes used to develop MSPC charts are two-way arrays, X(I×J), where J is the number of variables (or PCs after PCA processing) and I is the number of observations (i.e., time). In batch processes, for each variable at each observation, because its values correspond to a trajectory that spans the entire batch run (or campaign), the data to be analyzed is a three-way array X(I×J×K), where J is the number of variables measured at K time intervals throughout a batch, and I is the number of batch runs. The MSPC was adapted to such three-way data, using a technique called multi-way PCA.18 In multi-way PCA, the data is unfolded to a two-way array (I×JK), so that PCA analysis can be conducted. There are several variations of the MSPC models for process monitoring, depending on if the PCA method is linear or nonlinear,20,21 the way the data are scaled,22 as well as whether it is used in conjunction with other techniques such as wavelets23 and inductive data mining.24 For batch processes, other methods for three-way data compression have also been investigated, using parallel factor analysis (PARAFAC) and Turker3 and compared with multi-way PCA.25-30 Comprehensive review of the work has been performed by many researchers, and readers are referred to these review papers and books for further details of the methods.2-7,31,32 3. Independent Component Analysis 3.1. The Algorithm. Statistical independence is a different concept from, and has stricter conditions than, the uncorrelatedness; the latter is the basis for PCA. Denote x1, x2, ..., xm as random variables with the joint probability density function (pdf) of p(x1, x2, ..., xm) and assume that these variables have zero-mean, then they are said to be mutually statistically independent if the following condition holds:33-35

p(x1, x2, ..., xm) ) p1(x1)‚p2(x2) ‚ ‚ ‚ pm(xm)

(3)

where pi(xi) (for i ) 1, 2, ..., m) denotes the marginal pdf of xi, i.e., the pdf of xi when it is considered alone, (e.g., p1(x1) ) ∫‚ ‚ ‚∫p(x1, x2, ‚ ‚ ‚, xm) dx2 dx3 ‚ ‚ ‚ dxm). The uncorrelatedness implies that E{xixj} ) E{xi}‚E{xj} for i * j, where E{‚ ‚ ‚} represents the mathematical expectation. If xi and xj are independent, for any functions fi and fj, they must have

∫ ∫fi(xi)fj(xj)p(xi,xj) dxi dxj ) ∫ ∫fi(xi)fj(xj)pi(xi)pj(xj) dxi dxj ) ∫fi(xi)pi(xi) dxi ∫fj(xj)pj(xj) dxj

E{fi(xi)fj(xj)} )

) E{fi(xi)} E{fj(xj)}

(4)

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004 6733

If we denote the random column vectors as x ) [x1 x2...xl]T and s ) [s1s2...sm]T and we denote M as the matrix with elements mij (i ) 1, 2, ..., l; j ) 1, 2, ..., m) (namely, the mixing matrix), the aforementioned equation can be written as

random (i.e., unpredictable and unstructured) the variable, the larger its entropy. The differential entropy (H) of a random variable y with density p(y) is defined as H(y) ) - ∫p(y) log p(y) dy. A fundamental result of information theory is that a Gaussian variable has the greatest entropy among all random variables of equal variance. This means that entropy could be used as a measure of non-Gaussianity. To obtain a measure of non-Gaussianity that is zero for a Gaussian variable and always non-negative, one often uses a slightly modified version of the definition of differential entropy: negentropy, which is defined as J(y) ) H(ygauss) - H(y), where ygauss is a Gaussian random variable of the same covariance matrix as y. Because of the aforementioned properties, negentropy is always non-negative, and it is zero if and only if y has a Gaussian distribution. The negentropy J can have the following approximation:

x ) Ms

J(y) ≈ k[E{G(y)} - E{G(ν)}]2

and, therefore, if the variables are independent, they are uncorrelated; however, the converse is not true. In other words, the blindly separated ICs are statistically independent and also uncorrelated. In contrast, the PCs obtained from the linear combination of the observed variables are uncorrelated but do not satisfy the condition of statistical independence. Assume that there are l linear mixtures x1, x2, ..., xl of m independent (source) components s1, s2, ..., sm, i.e.,

xi ) mi1s1 + mi2s2 + ... + mimsm (i ) 1, 2, ..., l, where l g m) (5)

(6)

(9)

The solution of ICA consists of the estimation of the mixing matrix M and/or the independent source vector s from the observed mixed vector x. This solution is equivalent to finding a separating matrix W that satisfies the relation

where G is practically any nonquadratic function, k an irrelevant constant, and ν a Gaussian variable of zeromean and unit variance. Therefore, to find one independent component, one can maximize the function JG:

sˆ ) Wx

JG(w) ) [E{G(wTx)} - E{G(ν)}]2

(7)

where sˆ is the estimation of s. The separating matrix W in eq 7 can be trained as the weight matrix of a two-layer feed-forward neural network in which x is the input and sˆ is the output. The constraint of the neural network is that the elements of sˆ are statistically independent, which can be reflected by non-Gaussianity. Thus, a measure of nonGaussianity can be used as the objective (contrast) function when training the neural network. Note here that non-Gaussianity is a requirement for the elements of sˆ , i.e., the source signals or independent components. There is no such requirement for the sensor data. There are two common measures of non-Gaussianity: one is kurtosis, and the other is negentropy. A brief introduction to these two measures is given below. The kurtosis of a random variable y is conventionally defined by

kurt(y) ) E{y4} - 3(E{y2})2

(8)

For simplicity, assume that the variable y has been preprocessed to have zero-mean and unit variance; therefore, the kurtosis is simplified to kurt(y) ) E{y4} - 3. This shows that kurtosis is simply a normalized version of the fourth moment E{y4}. For a Gaussian y, the fourth moment is equal to 3(E{y2})2. Thus, kurtosis is zero for a Gaussian random variable. For most (but not quite all) non-Gaussian random variables, kurtosis is nonzero. Kurtosis can be either positive or negative. Random variables that have a negative kurtosis are called sub-Gaussian (meaning that the pdf is flatter than that of Gaussian), and those with positive kurtosis are called super-Gaussian (meaning that the pdf is more undulating than that of Gaussian). The kurtosis measure is a special case of negentropy, which is described in the following discussion. Negentropy is based on the information-theory quantity of (differential) entropy. The entropy of a random variable can be interpreted as the degree of information that the observation of the variable gives. The more

(10)

where w, which is an l-dimensional (weight) vector in the weight matrix W, is constrained so that E{(wTx)2} ) 1. Several ICs can then be estimated one by one; these components will be introduced later. According to Hyva¨rinen,36,37 the function G has the following choices:

G1(u) )

1 log cosh(a1u) a1

(11)

a2u2 1 exp a2 2

(12)

G2(u) ) -

G3(u) )

(

(41)u

)

4

(13)

where 1 e a1 e 2 and a2 ≈ 1. Among these three choices, G1 is a good general-purpose function. When the ICs are highly super-Gaussian, or when robustness is very important, G2 may be better. G3 is actually based on the kurtosis and is justified on statistical grounds only for estimating sub-Gaussian ICs when there are no outliers. Hyva¨rinen and Oja34 developed a fast fixed-point algorithm (FastICA) that can be described as follows: (1) Choose an initial random weight vector w(0) and let k ) 1. (2) Let w(k) ) E{vg(w(k - 1)Tx)} - E{g′(w(k - 1)Tx)}w(k - 1), where g is the first-derivative of the function G. The expectations can be estimated using a large sample of x vectors by computing over successive points. (3) Let w(k) ) w(k)/||w(k)||. (4) If the value of |w(k)Tw(k - 1)| is not close enough to 1, let k ) k + 1 and go back to step 2. Otherwise, output the vector w(k). This methodology describes the procedure to estimate just one of the ICs. To estimate several ICs, the procedure should be performed as many times as required. To prevent different vectors from converging to the same maxima, the outputs wT1 x, wT2 x, ..., wTmx

6734


after the every run of the procedure must be decorrelated. After the p vectors have been estimated, w1, w2, ..., wp, this just simply needs two more steps for wp+1: p

let wp+1 ) wp+1 let wp+1 )

T wp+1 wj wj ∑ j)1

wp+1 T wp+1 xwp+1

(14) (15)

The independence of random variables can be reflected by non-Gaussianity, which can be measured by kurtosis or negentropy.33,34 The process of ICA is an optimizing process of maximizing the non-Gaussianity. A few algorithms have been developed based on this principle.33-36,38 In this work, we used the fast fixed-points (FastICA) algorithm developed by Hyva¨rinen and Oja.34 The details of the algorithm can be found in their report34 and, therefore, will not be repeated here. Although ICA can be considered to be a useful extension of PCA, in fact, in the first step of ICA analysis (i.e., whitening of the signals), PCA can be used. However, the objective of ICA is very different from that of PCA. PCA attempts to determine PCs that are uncorrelated and are linear combinations of the observed variables, and PC1 is often more important than PC2, PC2 is more important than PC3, and so forth. However, ICA is designed to separate the ICs that are independent and constitute the observed variables, and there is no assumption that one IC is more important than others. In practice, before the ICA algorithm is applied, the observed data are often preprocessed to remove the correlation between the observed variables, which is called whitening. Several methods for this purpose have been developed, including PCA and factor analysis. The FastICA algorithm uses PCA as the whitening method. Currently, there is no better method available to determine the optimum number of ICs automatically. In this paper, the number of ICs is determined by the number of PCs chosen, using various methods as reviewed by Valle et al.39 Various linear and nonlinear ICA algorithms have been proposed in the literature in the past few years. The main difference between them is the way that the nonlinearity of the model and the criteria to optimize the solution are represented. Some of the proposed approaches include the use of parametric sigmoidal nonlinearity and higher-order polynomials,40 linearization approach as well as multilayer perceptron neural networks,41 and self-organizing maps.42 3.2. Previous Work on Applying Independent Component Analysis to Process Monitoring. In our early work, we described the application of ICA to chemical processes for the purpose of dimension reduction of the monitored variables43 and explored its use in chemometrics for near-infrared (NIR) spectral data analysis.44 Investigations on ICA for process monitoring were recently reported by Xia and Howell,8 Xia,9 Lee at al.,11,16 Yao and co-workers,12,13 and Kano et al.14 Xia and Howell8 and Xia9 proposed an interesting approach for process monitoring; this method is called spectral ICA. The method first transforms the process measurements from the time domain into the frequency domain, using discrete Fourier transform. The main advantage of analyzing the problem in the frequency

domain is that the power spectrum is invariant to time delays or phase shifts that are caused by process dynamics. ICA then is applied to the power spectra to identify ICs that correspond to the major oscillations such as faults. The spectral ICA approach was inspired by the spectral PCA method that was proposed by Thornhill et al.45 but was found to be more powerful than spectral PCA in discriminating multiple oscillations. Xia and Howell8 and Xia9 presented a convincing case study that demonstrated that spectral ICA can completely separate two oscillations imposed on the process, whereas spectral PCA cannot. Kano et al.14 proposed the idea of process monitoring based on observing the ICs instead of the original measurements, and devised SPC monitoring charts for each IC. The difference from traditional univariate SPC charts is that the ICs were used instead of the original variables, and the way of devising the upper and lower control limits (UCL and LCL, respectively) remained the same as that for traditional SPC. In traditional univariate SPC charts (and, in fact, also for multivariate SPC charts) are based on the assumption that the variables monitored follow Gaussian distribution under normal operation. Based on this assumption, UCL and LCL can be defined as (σ, (2σ, or (3σ, where σ is the standard deviation. Since the basis for ICA separation is that the ICs satisfy non-Gaussianity, they have relaxed the use of ICA. Nevertheless, they show that ICA-based SPC gave good results. In several publications, Lee et al.11,16,46 proposed a method for statistical process monitoring based on ICA that uses kernel density estimation to define the control limits of the ICs that do not satisfy a Gaussian distribution. Three statistics are used, including Id2 and Ie2 (equivalent to T2 in PCA-based MSPC) as well as the squared prediction error (SPE). Lee et al.16 extended their method to multi-way ICA for monitoring batch processes; the method combines ICA and kernel estimation to estimate the control limits. Their methodology requires unfolding the original matrix two times. The normal batch data is unfolded in the same way as the approach of Nomikos and MacGregor18 and then mean centering is performed to remove the batch process trajectory. After the batch trajectory is eliminated, the unfold matrix is rearranged into a form previously described by Wold et al.47 (i.e., (batch × time) × variables). Yao and co-workers12,13 investigated ICA for the purpose of representing the state space using ICs, and they subsequently applied the method to process monitoring. 4. The Method for Statistical Batch Process Monitoring Based on Independent Component Analysis In batch process operations, unlike in continuous processes, the process variables do not satisfy a Gaussian distribution along the time axis. This is because, in batch operations, there are different stages, including steps such as charge and discharge of feed and intermediate products, valve closing and opening, and variations in equipment topology. As a result, the trajectory of an operational variable does not follow Gaussian distribution along the time axis. However, if multiple normal batches that produce the same products using the same feeds and operational strategies are considered, the trajectories of a variable


for all the batches should follow a Gaussian distribution, centered around a mean, which we can call “the mean trajectory”. To simplify the discussion here, we have assumed all the batches have the same batch length. Following this line of thinking, it is also fair to assume that, for a specific time point t, the values at t for all the trajectories of a variable corresponding to all batches also satisfy a Gaussian distribution. As a result, a UCL and LCL, as well as a mean µ, can be obtained for a time t. Similarly, UCL, LCL, and µ can be calculated for all other times. By plotting all UCL, LCL, and µ values together, time-varying curves for UCL, LCL, and µ can be obtained. The proposed SPC monitoring charts using ICA is based on the aforementioned idea. However, the ICs are used instead of the original variables. The idea does not require the variable or IC to satisfy a Gaussian distribution along the time axis (in fact, neither the original variables, nor the ICs satisfy a Gausian distribution along the time axis). However, the values at a specific time t for an IC for all batches satisfy a Gaussian distribution; therefore, UCL, LCL, and µ can be calculated. By plotting all the UCLs, LCLs, and µ values together for that IC, we can obtain time-varying UCL, LCL, and µ curves that cover the time span of a batch run. In the following, the method is introduced in more detail by reference to Figure 1. The original three-way data matrix (J×I×K), where J is the number of variables, I the number of batch runs, and K the sampling points, is unfolded to a two-way matrix X(J×IK). The data are then normalized for each variable, i.e., for each row of the matrix (J×IK), using mean centering by subtracting the mean. ICA is applied to the matrix (J×IK) to compress it to Y(L×IK), where L is the number of ICs (L < J), and to obtain the separating matrix W. The new matrix Y(L× IK) then is split into L matrices, each corresponding to one IC and having the dimension of (I×K), where I is the number of batch runs and K the sampling points. Each column of the matrix (I×K), yik (for k ) 1, 2, ..., K) is assumed to satisfy a Gaussian distribution and, therefore, UCL, LCL, and µ can be calculated. The UCLs and LCLs over all the sampling points will form the UCL and LCL curves. For the purpose of monitoring a new batch operation, the Xnew measurements are normalized by subtracting the mean previously obtained in developing the ICAbased SPC model, and then the values for the ICs are calculated by W×Xnew, where W is the separating matrix obtained during the ICA-based SPC model development stage. Finally, the instant values that correspond to each IC is plotted in the SPC diagram to show if the control limits have been exceeded.

Figure 1. Schematic depiction of the proposed approach. Legend is as follows: J, number of process variables; I, number of batch runs; K, number of sampling points, assuming that all batch runs have the same number of sampling events; ICA, independent component analysis; ICs, independent components; UCL, upper control limit; and LCL, lower control limit.

5. Case Studies 5.1. Case Study 1. The methodology is illustrated by reference to a case study of a semibatch polymerization reactor that was manufacturing polyol lubricant by condensation of an alkylene oxide in the presence of an alcohol, according to the following reaction:

C4H9OH + (n + 1)C3H6O f C4H9(OC3H6)nOCH2CHOHCH3 + heat The reaction is highly exothermic, and the reactor contains large quantities of volatile oxide. A schematic representation of the process is shown in Figure 2. The

Figure 2. Schematic of the polymerization reactor.

batch contents are recycled through an external cooler with sufficient capacity to maintain a constant return temperature to the reactor. The reactor is charged with the catalyzed starting material, which only occupies a very small proportion of the final batch volume. Oxide

6736


Table 1. Process Variables Monitored for the Semibatch Polymerization Reactor: Case Study 1 variable number

process variable

1

C

2 3 4 5 6 7 8 9

P T Two Q M X Ff Tc

10

Fw

description alkylene oxide concentration in reactor (kg/kg) reactor pressure (bar) reactor temperature (°C) outlet temperature of coolant (°C) rate of cooling (kJ/min) mass in reactor (kg) amount of reacted oxide feed rate of oxide (kg/min) outlet temperature of the circulating reaction mixture (°C) mass flow rate of coolant (kg/min)

is then continuously fed into the reactor until the batch is made and the reactor is full. Very careful control of the reaction temperature, pressure, and heat-removal rate is vital to ensure that catastrophic “runaway” polymerization does not occur.24,48 Monte Carlo simulation is combined with the dynamic simulator to generate the data set that corresponds to 40 batch runs with random disturbances being introduced into the data. The 10 process variables (Table 1) were sampled at 240 intervals for batch monitoring. Therefore, the data form a three-way matrix (J×I×K) ) (10 × 40 × 240). X is unfolded to a two-way matrix X(10 × 9600). Each batch run has three different stages. At the first stage, reactor pressure increases quickly, because of the increase of alkylene oxide concentration. The second stage is called reaction retaining, in which the temperature in reactor is kept stable. At the last stage, the measurement trajectories show obvious changes that are due to depletion of the alkylene oxide feed (Figure 3). 5.2. Results and Discussion. The unfolded matrix X(J×IK) ) X(10 × 9600) is normalized for each variable, i.e., for each row using mean entering, and then processed by ICA to obtain the compressed matrix Y(L×IK) ) Y(L × 9600), where L is the number of ICs (L < 10), and the separating matrix W. It was determined that four ICs must be used; i.e., the new matrix is Y(L×IK) ) Y(4 × 9600). Because the data preprocessing step of the FastICA algorithm is PCA, the number of ICs can be determined by the number of PCs chosen, using an appropriate method. Currently, there are no better methods for selecting the number of ICs. Valle

Figure 3. Plot of heat release versus degree of batch completion; Q1 (triangles) represents experimental data, whereas Q2 (circles) denotes data generated by the simulator. Legend: I, stage one; II, stage two; and III, stage three.

Figure 4. Plot of the variance values corresponding to the 10 principal components (PCs).

Figure 5. Gaussian probability test for the values of 40 batches for a time point of IC3.

et al.39 reviewed and compared various methods for determining the number of PCs and concluded that there is no definitive answer in regard to the selection of approaches. For instance, one method recommends that the selected PCs should represent no less than 97% of the variation in the data. However, in some cases, this may lead to too many PCs being selected, including those representing noise. For the current case, Figure 4 shows that the four PCs can retain 97% of the variation in data; therefore, four ICs were used. In the second step, the matrix Y(4 × 9600) (i.e., (4 × 40 × 240)) is split into four matrices, with each matrix corresponding to one IC and having a size of (40 × 240, the number of batches × the number of time points). For this matrix of (40 × 240), the UCL and LCL limits and µ value can be calculated for each column, i.e., corresponding to each point in time. Linking all the time points together, a UCL curve, a LCL curve, and a µ curve can be obtained, which cover the evolution time points of a batch run. To test the assumption that the elements of each column of the matrix (40 × 240, the number of batches × the number of time points) for an IC satisfy a Gaussian distribution, the Gaussian probability plot test was conducted. In a Gaussian probability plot, the data are plotted against a theoretical normal distribution in such a way that the points should form an approximately straight line. Departures from this straight line indicate departures from normality. Figure 5 shows the test of a Gaussian distribution plot of one time point for IC3. The data can be said to satisfy a Gaussian


Figure 6. Upper control limit (UCL), lower control limit (LCL), and mean between UCL and LCL (µ) for all independent components (ICs), plotted for three stages of the batch operation: (a) IC1, (b) IC2, (c) IC3, and (d) IC4.

distribution reasonably well, because all the data points are within the confidence limit lines of 95%. The resulting control charts for the four ICs are shown in Figure 6. This figure reflects the differences of the operations in three stages. Time point 0 to time point 60 represents the first stage, in which reactor pressure increases quickly, because of the increase of alkylene oxide concentration. The second stage, which is covered by time points 61-225, is the reaction retaining stage. In the third stage, from time point 226 to time point 240, the measurement trajectories have obvious changes, which are due to depletion of the alkylene oxide feed. In panels a-d of Figure 6, the UCL and LCL limits are calculated based on (3σ. To illustrate how to use the ICA-based SPC charts for process monitoring, we test another batch run, which we will call batch 41. During the operation, a fault was introduced. Display of the batch 41 data on the four SPC charts is given as dotted curves in panels a-d of Figure 7. Batch 41 clearly behaves within the UCL and LCL limits at all stages in the SPC charts for IC1 and IC2 but goes beyond the UCL and LCL limits in the second and third stages of IC3 and only in the first stage for IC4.

Figure 7. Batch run 41 (dotted curve) was detected as abnormal on the IC3 control chart at stages 2 and 3 and on the IC4 control chart at stage 1: (a) IC1, (b) IC2, (c) IC3, and (d) IC4.

To find the cause of abnormal operation, contribution plots were used. The contribution plot was computed following the method of Wold et al.47 by taking into account any normal batch as a reference point or average point for contribution plotting and comparing it with the abnormal batch at a local batch time point. For example, at local time point 90 and considering normal batch 1 as a reference, we see the difference between the normal batch 1 and abnormal batch 41; the difference between variables is great and has a greater contribution to abnormal conditions. Here, only the contribution plots for IC3 and IC4 are needed, because no abnormality was observed on the IC1 and IC2 charts. Three sampling pointssNo. 17 at stage 1, No. 70 at stage 2, and No. 230 at stage 3swere selected for illustration purposes. Figure 8 is the contribution plot for IC3 at sampling point 70 (at stage 2) which indicates that the major contributing original variable is variable 9, i.e., the outlet temperature of the circulating reaction mixture (the numbering of the variables is given in Table 1). Figure 9 shows the contributing plot for IC3 at sampling point 230 (at stage 3), and Figure 10 is the contribution plot for IC4 at sampling point 17 (i.e., stage 1). Figures 8-10 suggest the same conclusion: i.e., the major contributing original variable is variable 9 (i.e., the outlet temperature of the circulating reaction mixture). In fact, the abnormal operation was due to a

6738


Figure 8. Contribution plot at sampling time point 70 of the control chart for IC3; variable 9 is the outlet temperature of the circulating reaction mixture.

Figure 11. Batch run 49 (dotted curve) was detected as abnormal on the IC1 chart at time points 58-62, because the data were above the upper control limit (UCL), and at time points 63 and 64, because the data were below the lower control limit (LCL).



biased error in the sensor measuring the outlet temperature of the circulating reaction mixture. 5.3. Case Study 2. To test the capability for generalization of the method to industrial data, we also applied it to the analysis of a set of industrial data. The data set involves an industrial batch polymerization reactor of Dupont and had been studied previously by Nomikos and MacGregor49 and Eriksson et al.50 The process consists of two stages, and the operation at each stage lasts ∼2 h. Ingredients were loaded into the reactor to begin the first stage, and the reactor medium flows are adjusted to establish proper control of the pressure and the rate of temperature change. The solvent used to bring the ingredients to the reactor was vaporized and removed from the reactor vessel. The second stage began ∼1 h after the solvent was removed. During this stage, the ingredients completed their reaction to yield the final product. The vessel pressure and the rate of the temperature change also were controlled during this stage. The batch finishes by pumping the polymer product from the reaction vessel at the end of the second stage. The data set contains 18 batches, which have been identified in the literature as normal operations, and another 7 batches (which have been named batches 4955), which were identified as abnormal batches. Each batch had a duration of 100 time intervals and 10

variables. Variables 1, 2, and 3 are temperature measurements inside the reactor; variables 6 and 7 are temperature measurements in the heating-cooling medium; variables 4, 8, and 9 are pressure measurements; and variables 5 and 10 are flow rates of materials added to the reactor. Our purpose here is to demonstrate if the proposed method can quickly identify the fault batches. The threeway data about normal operations were unfolded to twoway data (10 × 1800) and normalized for each variables (i.e., for each row), using mean centering. ICA then was applied. It was determined that three ICs were needed, which captured 85% of the variation in the data. In the next step, UCL and LCL limits and the µ value then were determined for each time point. New batch data (for bathes 49-55) can be projected to the model and tested to determine if these batches exceed the control limit. The UCL, LCL, and µ values were used to test if the seven batchessbatches 49-55swere normal. We found that all seven batches were abnormal, which is consistent with the findings of Eriksson et al.50 For illustrative purposes, only an example is shown here. Figure 11 shows the UCL, LCL, the µ value for IC1, and the trajectory of batch 49 in the control chart for IC1. The figure clearly shows that the operation went above the UCL around the time points 55 and 60, and went below the LCL around time points 63 and 64. 5.4. Case Study 3. Case study 3 is a simulated process for batch-fed penicillin fermentation, which was used by Yoo et al.51 in presenting their multi-way ICA method for batch process monitoring. The simulator, PenSim v2.0, was developed by the Monitoring and Control Group of the Illinois Institute of Technology and can be downloaded from the website http:// www.chee.iit.edu/∼cinar. It simulates the concentrations of biomass, CO2, H+ ion, penicillin, carbon source, and oxygen, as well as heat generation, during the production of penicillin under various operating conditions. Table 2 gives a summary of the 10 variables that are recorded. We have used the simulator to generate 15 batch runs that correspond to normal operation, and another two batch runs that have a step decrease in the glucose feed rate (named fault 1) and a linear decrease in the glucose feed rate (named fault 2). Details for the data are given in Table 3. Each batch had a duration of 400 time points. The data contain random noise. Although we could have removed the noise components, as Yoo et al.51 did, we conducted the analysis based on

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004 6739 Table 2. Variables Recorded for the Penicillin Fermentation Simulation: Case Study 3 variable number

process variable

1 2 3 4 5 6 7 8 9 10

aeration rate (L/h) agitator power (W) glucose feed temperature (K) dissolved oxygen concentration (% saturation) culture volume (L) carbon dioxide concentration (mmol/L) pH temperature (K) generated heat (kcal) cooling water flow rate (L/h)

Table 3. Data for the Penicillin Fermentation Process: Case Study 3 case

explanation

normal batches fault 1: step fault

15 batches glucose feed rate is suddenly step-decreased by 15% at 50 h and maintained at that lower rate to the end of batch operation (400 h) glucose feed rate is linearly decreased from 0.04 L/h to 0.03 L/h from 60 h to the end of batch operation (400 h)

fault 2: ramp fault

the raw data without the noise removed, because it did not affect the performance of our method, as will be illustrated. We will use this case study to compare the performance of the proposed method in this paper to a PCAbased method previously proposed by Wold et al.47 and the multi-way ICA method that was proposed by Yoo et al.51 The procedures for developing the UCL and LCL control limits, finding the mean µ, and determining the number of ICs are the same as those for the previous two case studies; therefore, they will not be repeated. Here, we focus the discussion on presenting some representative results and the comparison to other approaches. Two ICs were chosen to develop the monitoring charts. Figure 12 shows the monitoring charts based on ICA proposed in this paper (Figure 12a) and the approach proposed by Wold et al.47 (Figure 12b). The fault, which is a step decrease by 15% in the glucose feed rate at 50 h and then maintained at that lower rate until the end of fermentation, was detected at ∼124 h by the proposed method in this paper (see Figure 12a); however, the fault was detected by the PCA approach only after ∼245 h and the deviation is not as clear as in the ICA method, as evidenced by Figure 12b. A similar result was observed for the second fault, as indicated in Figure 13a, which shows the monitoring chart based on the method of this paper using ICA, and Figure 13b, which shows the chart based on PCA. The fault of a linear decrease in glucose feed rate was initiated at 60 h and was detected by the method in this paper after ∼122 h, and that detected by the PCA method was observed after ∼165 h. In other words, the current method can detect the fault 23 h earlier than the PCA-based approach. The data used here were not exactly the same as those used by Yoo et al.:51 we simulated 15 normal data cases and Yoo et al.51 used 60. Therefore, we can only make a qualitative comparison between our approach and the method of Yoo et al. Nevertheless, the two faults used here were also discussed in the work of Yoo et al. It was determined that both methods could detect the fault

Figure 12. Monitoring charts based on (a) the method in this paper using ICA and (b) PCA by Wold et al.47 These plots indicate that the method in this paper detects the fault of a step decrease in the glucose feed rate earlier than the PCA-based approach. The UCL and LCL bounds were set at 99.7% for both methods.

Figure 13. Monitoring charts based on (a) the method in this paper using ICA and (b) PCA by Wold et al.47 These plots indicate that the method in this paper detects the fault of a linear decrease in the glucose feed rate earlier than the PCA-based approach. The UCL and LCL bounds were set at 99.7% for both methods.

almost at the same time. They compared their approach with PCA-based T2 and SPE models and concluded that their approach detected the faults earlier than PCAbased T2 and SPE models. Thus, it is reasonable to conclude that our method detected the fault earlier than PCA-based T2 and SPE models. Although both the current approach and the method of Yoo et al. detected the two faults almost at the same time, the current approach has an advantage. Yoo et al.’s approach51 actually did not detect the two faults in their I2 statistic. The two faults were only detected

6740


by the SPE statistic. This leads to a slight disadvantage, because Louwerse and Smilde,27 and ourselves,52 have found that, unlike T2 charts (or I2 charts), the SPE control limit is very difficult to set to either 99% or 99.9%. Therefore, depending on SPE to detect the fault will be disadvantageous. 6. Final Remarks We have described a new method for statistical process control (SPC) of batch operations which allows the real-time monitoring of batch evolution using a reduced number of variables, i.e., the independent components (ICs), instead of the original variables. The number of ICs is dependent on the complexity of the problem and can be determined quantitatively. The upper and lower control limits (UCL and LCL, respectively), as well as the mean (µ) of the SPC charts for each IC are calculated based on the assumption that the trajectories for a specific IC over all the batches of normal operation satisfy a Gaussian distribution, and so do the values at any given time point for all trajectories of the IC. This is considered to be a reasonable assumption as long as all the batches were created for manufacturing the same products under the same operation conditions, although the values of an IC for a batch run do not follow a Gaussian distribution along the time axis. The method is different from the multivariate SPC method of Nomikos and MacGregor18 and some other methods derived from Nomikos and MacGregor’s idea, not only because the latter are based on principal component analysis (PCA) and always uses two statistics (Hotelling’s T2 parameter and squared prediction error (SPE)), but also because the latter are more focused on batch-to-batch tracking (i.e., on classifying the completed batches as normal or abnormal). The proposed approach in this work has some common grounds with the method of Wold et al.,47 because both are used to track the batch evolution from time point to time point; however, the method of Wold et al. is based on the scores of PLS analysis of data. In the current method, the number of ICs that are monitored is dependent on the complexity of the problem, and, for each IC, the monitoring chart has UCL and LCL bounds. The method of Lee et al.,11,16,46 which is based on kernel estimation, used a fixed number of three statistics, i.e., Id2 and Ie2 and SPE. Id2 and Ie2 are equivalent to the T2 chart and also only have a UCL. Although the method of Kano et al.14 also uses SPC charts and UCL and LCL control limits for ICs, it has been developed for continuous process operations and compromises the use of ICs, because they used the traditional univariate SPC approach to calculate the UCL and LCL limits for the ICs, which essentially assumes that each IC follows a Gaussian distribution along the time axis. In addition, the work of Kano et al. gave UCL and LCL bounds that are constant over the time dimension. The method we proposed has a different way of calculating the UCL and LCL limits, and the control limits are time-varying over the batch evolution. The difference between the current method and the spectral independent component analysis (ICA) idea of Xia and Howell8 and Xia9 is more obvious, because the spectral ICA is based on analysis of the power spectra and from the work of Yao and co-workers12,13 is to represent the operation in a reduced dimensional space. The new approach was illustrated by reference to three case studies. Case study 1 was a semibatch

polymerization process that we simulated, which was described in detail to illustrate the procedure. Case study 2 was a data set collected from industry for a semibatch polymerization process that had been studied previously; it proved that the new method can detect the seven faults in a straightforward manner. Case study 3 was a simulated process that was also studied in the literature; it proved that the proposed approach can detect the two faults earlier than the method of Wold et al.47 Although both the current approach and the multi-way ICA method of Yoo et al.51 detected the two faults almost at the same time points, we argue that the current method has a slight advantage, because Yoo et al.’s approach actually did not detect the fault using the I2 statistic, but instead relied on SPE, which is known to have difficulty in deciding the most appropriate control limit. Literature Cited (1) Kresta, J. V.; MacGregor, J. F.; Marlin, T. E. Multivariate Statistical Monitoring of Process Operation Performance. Can. J. Chem. Eng. 1991, 69, 35-47. (2) Cinar, A.; Undey, C. Statistical Process and Controller Performance Monitoring: A Tutorial on Current Methods and Future Directions. In Proceedings of the American Control Conference (San Diego, CA, 1999); Vol. 4, pp 2625-2639. (3) MacGregor, J. F.; Kourti, T. Statistical Process Control of Multivariate Processes. Control Eng. Pract. 1995, 3, 403-414. (4) Wise, B. M.; Gallagher, N. B. The Process Chemometrics Approach to Process Monitoring and Fault Detection. J. Process Control 1996, 6, 329-348. (5) Russell, E.; Chiang, L. H.; Braatz, R. D. Data-Driven Methods for Fault Detection and Diagnosis in Chemical Processes; Springer: London, 2000. (6) Qin, S. J. Statistical Process Monitoring: Basics and Beyond. J. Chemom. 2003, 17, 480-502. (7) Edgar, T. F. Control and Operations: When Does Controllability Equal Profitability?. In Process Systems Engineering 2003: 8th International Symposium on Process Systems Engineering, China; Chen, B., Westerberg, A. W., Eds.; Computer-Aided Chemical Engineering 15; Elsevier: Amsterdam, 2003; pp 4861. (8) Xia, C.; Howell, J. Isolating Multiple Sources of Plant-Wide Oscillations via Independent Component Analysis. In SAFEPROCESS2003-5th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes, Washington, DC, 2003. (9) Xia, C. Control Loop Measurement Based Isolation of Faults and Disturbances in Process Plants, Ph.D. Thesis, University of Glasgow, Glasgow, U.K., 2003. (10) Lee, J. M.; Yoo, C.; Lee, I. B. Statistical Process Monitoring with Multivariate Exponentially Weighted Moving Average and Independent Component Analysis. J. Chem. Eng. Jpn. 2003, 36, 563-577. (11) Lee, J. M.; Yoo, C. K.; Lee, I. B. New Monitoring Technique with an ICA Algorithm in the Wastewater Treatment Process. Water Sci. Technol. 2003, 47, 49-56. (12) Yao, Z. X.; Qian, Y.; Li, X. X.; Jiang, Y. B. A Description of Chemical Processes Based on State Space. In Process Systems Engineering 2003: 8th International Symposium on Process Systems Engineering, China; Chen, B., Westerberg, A. W., Eds.; Computer-Aided Chemical Engineering 15; Elsevier: Amsterdam, 2003; pp 1112-1117. (13) Yao, Z. X. Statistical Modelling of the State Space and Monitoring of Chemical Process Systems, Ph.D. Thesis, South China University of Technology, Guangzhou, PRC, 2003. (14) Kano, M.; Tanaka, S.; Hasebe, S.; Hashimoto, I.; Ohno, H. Monitoring Independent Components for Fault Detection. AIChE J. 2003, 49, 969-976. (15) Li, R. F. Advanced Process Monitoring and Control Using Principal and Independent Component Analysis, Ph.D. Thesis, University of Leeds, Leeds, U.K., 2003. (16) Lee, J. M.; Yoo, C.; Lee, I. B. On-line Batch Process Monitoring Using Different Unfolding Method and Independent Component Analysis. J. Chem. Eng. Jpn. 2003, 36, 1384-1396.

Ind. Eng. Chem. Res., Vol. 43, No. 21, 2004 6741 (17) Lee, J. M.; Yoo, C. K.; Lee, I. B. Statistical Process Monitoring with Independent Component Analysis. J. Process Control 2004, 14, 467-485. (18) Nomikos, P.; MacGregor, J. F. Monitoring Batch Processes Using Multi-way Principal Component Analysis. AIChE J. 1994, 40, 1361-1375. (19) Kourti, T. a. J. F. M. Process Analysis, Monitoring and Diagnosis, Using Multivariate Projection Methods. Chemom. Intell. Lab. Syst. 1995, 28, 3-21. (20) Kramer, M. A. Nonlinear Principal Component Analysis Using Autoassociative Neural Networks. AICHE J. 1991, 37, 233. (21) Dong, D.; McAvoy, T. J. Nonlinear Principal Component AnalysissBased on Principal Curves and Neural Networks. Comput. Chem. Eng. 1996, 20, 65-78. (22) Cox, T. F. Multidimensional Scaling Used in Multivariate Statistical Process Control. J. Appl. Stat. 2001, 28, 365-378. (23) Bakshi, B. R. Multiscale PCA with Application to Multivariate Statistical Process Monitoring. AICHE J. 1998, 44, 15961610. (24) Yuan, B.; Wang, X. Z. Multilevel PCA and Inductive Learning for Knowledge Extraction from Operational Data of Batch Processes. Chem. Eng. Commun. 2001, 185, 201-221. (25) Wise, B. M.; Gallagher, N. B.; Butler, S. W.; White, D. D., Jr.; Barna, G. G. A Comparison of Principal Component Analysis, Multi-way Principal Component Analysis, Trilinear Decomposition and Parallel Factor Analysis for Fault Detection in a Semiconductor Etch Process. J. Chemom. 1999, 13, 379-396. (26) Westerhuis, J. A.; Kourti, T.; MacGregor, J. F. Comparing Alternative Approaches for Multivariate Statistical Analysis of Batch Process Data. J. Chemom. 1999, 13, 397-413. (27) Louwerse, D. J.; Smilde, A. K. Multivariate Statistical Process Control of Batch Processes Based on Three-Way Models. Chem. Eng. Sci. 2000, 55, 1225-1235. (28) Kiers, H. A. L.; Berge, J. M. F. T.; Bro, R. PARAFAC2s Part 1. A Direct Fitting Algorithm for the PARAFAC2 Model. J. Chemom. 1999, 13, 275-294. (29) Bro, R.; Andersson, C. A.; Kiers, H. A. L. PARAFAC2s Part II. Modelling Chromatographic Data with Retention Time Shifts. J. Chemom. 1999, 13, 295-309. (30) Bro, R.; Kiers, A. L. A New Efficient Method for Determining the Number of Components in PARAFAC Models. J. Chemom. 2003, 17, 274-286. (31) Chiang, L. H.; Russell, E. L.; Braatz, R. D. Fault Detection and Diagnosis in Industrial Systems; Springer: London, 2001. (32) Wang, X. Z. Data Mining and Knowledge Discovery for Process Monitoring and Control; Springer: London, 1999. (33) Hyvarinen, A. The Fixed-Point Algorithm and Maximum Likelihood Estimation for Independent Component Analysis. Neural Process. Lett. 1999, 10, 1-5. (34) Hyvarinen, A.; Oja, E. Independent Component Analysis: Algorithms and Applications. Neural Networks 2000, 13, 411430. (35) Hyvarinen, A.; Oja, E. A Fast Fixed-Point Algorithm for Independent Component Analysis. Neural Comput. 1997, 9, 14831492. (36) Hyvarinen, A. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Trans. Neural Networks 1999, 10, 626-634.

(37) Hyva¨rinen, A. Survey on Independent Component Analysis. Neural Comput. Surveys 1999, 2, 94. (38) Hyvarinen, A. Independent Component Analysis in the Presence of Gaussian Noise by Maximizing Joint Likelihood. Neurocomputing 1998, 22, 49-67. (39) Valle, S.; Li, W. H.; 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. (40) Lee, T. W. In Stochastic Dynamics and Pattern Formation in Biological and Complex Systems; American Institute of Physics (AIP): Melville, NY, 2000; Vol. 501, pp 302-316. (41) Sole, J.; Jutten, C.; Taleb, A. Parametric Approach to Blind Deconvolution of Nonlinear Channels. Neurocomputing 2002, 48, 339-355. (42) Lin, J. K.; Grier, D. G.; Cowan, J. D. Faithful Representation of Separable Distributions. Neural Comput. 1997, 9, 13051320. (43) Li, R. F.; Wang, X. Z. Dimension Reduction of Process Dynamic Trends Using Independent Component Analysis. Comput. Chem. Eng. 2002, 26, 467-473. (44) Chen, J.; Wang, X. Z. A New Approach to Near-Infrared Spectral Data Analysis Using Independent Component Analysis. J. Chem. Inf. Comput. Sci. 2001, 41, 992-1001. (45) Thornhill, N. F.; Shah, S. L.; Huang, B.; Vishnubhotla, A. Spectral Principal Component Analysis of Dynamic Process Data. Control Eng. Pract. 2002, 10, 833-846. (46) Lee, M. J.; Yoo, C.; Lee, I. B. Statistical Process Monitoring with Multivariate Exponentially Weighted Moving Average and Independent Component Analysis. J. Chem. Eng. Jpn. 2003, 36, 563-577. (47) Wold, S.; Kettaneh, N.; Friden, H.; Holmberg, A. Modelling and Diagnosis of Batch Processes and Analogous Kinetics Experiments. Chemom. Intell. Lab. Syst. 1998, 44, 331-340. (48) Yuan, B. Process Data Mining Using Neural Networks and Inductive Learning, Ph.D. Thesis, The University of Leeds, Leeds, U.K., 2002. (49) Nomikos, P.; MacGregor, J. F. Multivariate SPC Charts for Monitoring Batch Processes. Technometrics 1995, 37, 41-59. (50) Eriksson, L.; Johansson, E.; Kettaneh-Wold, N.; Wold, S. Introduction to Multi- and Megavariate Data Analysis Using Projection Methods (PCA and PLS); Umetrics AB: Windsor, U.K., 1999. (51) Yoo, C. K.; Lee, J. M.; Vanrolleghem, P. A.; Lee, I. B. Online Monitoring of Batch Processes Using Multiway Independent Component Analysis, Chemom. Intell. Lab. Syst. 2004, 71, 151163. (52) Albazzaz, H.; Wang, X. Z.; Marhoon, F. Multidimensional Visualisation for Process Historical Data Analysis: A Comparative Study with Multivariate Statistical Process Control. J. Process Control 2004, in press.

Received for review May 17, 2004 Revised manuscript received August 6, 2004 Accepted August 9, 2004 IE049582+

Statistical Process Control Charts for Batch Operations Based on

Recommend Documents