Batch Process Monitoring Based on Multisubspace Multiway Principal

Article pubs.acs.org/IECR

Batch Process Monitoring Based on Multisubspace Multiway Principal Component Analysis and Time-Series Bayesian Inference Lv Zhaomin, Jiang Qingchao, and Yan Xuefeng* Key Laboratory of Advanced Control and Optimization for Chemical Processes of Ministry of Education, East China University of Science and Technology, Shanghai 200237, P. R. China S Supporting Information *

ABSTRACT: Multiway principal component analysis (MPCA), which is a dimensionality reduction method for process variables, has been widely used to monitor batch and fed-batch processes. However, three main factors affect the performance of MPCA monitoring: The future status of the online batch has to be predicted, the discarded principal components with small variance might contain useful information, and self-correlation and industrial noise exist in process data. Thus, a new batch process monitoring method based on multisubspace multiway principal component analysis and time-series Bayesian inference through a moving window is developed. The feasibility and effectiveness of the proposed batch process monitoring method is demonstrated using a numerical process and the fed-batch penicillin fermentation process, and its performance is compared with that of the MPCA. The results show that the proposed method is more accurate in detecting different types of batch process faults.

1. INTRODUCTION Batch and fed-batch processes serve an important function in the production of high-quality, low-volume products including specialty chemicals, food, biological industry products, and semiconductor devices. Ensuring both the consistency of the quality of products and the safe operation of the batch process is very important. The characteristics of batch processes, such as instability, finite duration, and batch-to-batch variations, differ from those of a continuous process. Establishing a monitoring model for batch processes is more difficult than building one for continuous processes. Meanwhile, biological batch processes, such as penicillin fermentation, often encounter slow recovery from faulty conditions to normal operation.1,2 A general feature of biological batch processes is that slight changes in operating conditions during critical periods can degrade the quality and yield of the final product. Some changes are initially unnoticeable and can gradually be exacerbated until they become a serious process fault. Therefore, the early detection of process faults is crucially important in batch process monitoring to enable the operator to take some control action to avoid the waste of products or serious damage to expensive process equipment.3 Multiway principal component analysis (MPCA), which is called batchwise MPCA in this article, and multiway partial least-squares (MPLS), developed by Nomikos and MacGregor, are the most popular methods used in batch process monitoring.4−8 However, when the batchwise MPCA and MPLS methods are used in the online monitoring of batch process, one main disadvantage is noted: The future status until the end of each batch must be estimated. Three methods for overcoming this disadvantage were summarized by Nomikos and MacGregor:6 (1) fill the missing values with zeros, (2) fill the missing values with the current value, and (3) use the capability of PCA to handle missing data.6 However, these estimates do not truly reflect the actual potential relationship in © 2014 American Chemical Society

the process and affect online monitoring performance. A historical batch database was used as the reference trajectory by Cho and Kim.9 To implement this database online, the distance between the known measurement data of the new batch and the corresponding parts in the historical batches has to be calculated first, after which the most similar batch is selected to fill the missing values. However, the necessary computation is very complicated and also cannot truly reflect the actual potential relationship in the process. A local MPCA model was established on marked time points by Louwerse and Smilde.10 A local MPLS prediction model was established according to the evolution of the process by Ü ndey et al.11 Both of these methods only reduce the calculation of status predictions, but cannot resolve the problem of future status prediction. A recursive hierarchical adaptive PCA algorithm was proposed by Rännar et a.12 This algorithm avoids the problem of status prediction by including the underlying time-slice PCA model and the upper PCA model based on principal component scores. However, the model structure is complex and requires a large amount of calculation. Moreover, the current model precision is seriously affected by the preceding model error in the iterative updating of the process. The algorithm exhibits numerous cumulative errors. MPCA based on variable-wise unfolding does not require future status prediction.13,14 The combination of batchwise and variable-wise methods has recently been used to build a two-dimensional modeling data unit that successfully overcomes the default prediction of future status.15−17 Aside from the problem of future status prediction, principal component selection requires further investigation. Numerous Received: Revised: Accepted: Published: 6457

October 23, 2013 March 24, 2014 March 24, 2014 March 25, 2014 dx.doi.org/10.1021/ie403576c | Ind. Eng. Chem. Res. 2014, 53, 6457−6466

Industrial & Engineering Chemistry Research

Article

2. PRELIMINARIES This section reviews batchwise MPCA6 and MPCA15 based on a combination of batchwise and variable-wise with time-varying score covariance structures for process monitoring, which is called variable-wise MPCA in this article. These methods are compared with the proposed method. 2.1. Fault Detection Using Batchwise MPCA. The data set of a batch process is a three-way array, X(I × J × K), where I represents the number of batches, J denotes the number of measurement variables, and K is the number of sampling instants within each batch. In batchwise MPCA, X is unfolded such that every slice I × J is placed along the K axis from the first time interval to the last time interval side by side. A twodimensional matrix XB(I × JK) is then constructed. All data are illustrated in Figure S1 of the Supporting Information. The variation among batches can be analyzed in this unfolding manner. Before modeling, each column in XB(I × JK) needs to be scaled to zero mean and unit variance, because less important variables with large magnitudes can overshadow important variables with small magnitudes. Then, matrix X̅ B(I × JK) is obtained. Considering the matrix X̅ B(I × JK) with I samples and JK variables, PCA is performed as follows

methods can be used to determine which principal components should be extracted, such as cumulative percent variance (CPV), cross validation, variance of reconstruction error (VRE), and fault signal-to-noise ratio (SNR).18−20 These methods select only the first several principal components with large variances, whereas the last several components with small variances are rejected. Jolliffe demonstrated that the principal components with small variances can be as important as those with large variances.21 Therefore, any abandoned components might deteriorate the process monitoring performance. Sensitive PCA (SPCA) was proposed by Jiang et al., who used the rate of change of T2 to select the principal components.22 Both components with large variances and those with small variances were selected. However, sensitive components were selected only at the moment when a fault occurred. Many components are rejected in SPCA, thus leading to information loss. A k + 1 subspace construction method was proposed by Ge et al., in which a linear subspace method is used for nonlinear process monitoring [i.e., block-specific principal component analysis (BSPCA)].23 They used the PCA decomposition method for linear subspace construction and divided the original process variables into k + 1 overlapping subsets, where k is the number of retained principal components. A four-subspace construction method was proposed by Tong et al., who used PCA to derive four subspaces from the original process variables.24 A squared Mahalanobis distance was calculated directly in each subspace for fault detection, thereby preserving the information of the original data space. However, the correlation information of some variables in different subspaces was not extracted in both the variable subspace construction methods. In this article, a new method based on multisubspace multiway PCA and time-series Bayesian inference (MMPCATSBI) through a moving window is developed for batch process monitoring. To avoid the problem of future status prediction, a combination of batchwise and variable-wise methods is employed to build two-dimensional modeling data. PCA is applied to all variables to extract all of the correlation information. To avoid information loss, all loadings are reserved in the proposed method. A denotative matrix is defined according to the loading matrix. According to the denotative matrix, all loadings are divided into several subspaces (two subspaces are employed in this article) with each subspace serving as a low-dimensional representation of the original data space. The separation rule is based on the clustering result of the denotative matrix, and the loadings in each subspace have similar characteristics. The proposed method can enhance the relevant information in the same subspace by grouping the loadings with similar characteristics. In this way, an abnormal event can be easily detected. TSBI through a moving window is employed as a moving weighting strategy to handle the selfcorrelation that exists between different intervals because the current behavior is affected by the behavior of past time in a dynamic process. Moreover, TSBI can overcome the noise in industrial processes. The remainder of this article is organized as follows: Section 2 presents two traditional methods for batch process monitoring. Section 3 presents detailed descriptions of the proposed method, MMPCA-TSBI, and the monitoring procedure. Section 4 discusses a numerical process and the fed-batch penicillin fermentation process to demonstrate the effectiveness of the proposed method in batch process monitoring. Finally, we present the conclusions of our work.

X̅ B = TPR T + E

(1)

where the loading matrix is PR(JK × R), the score matrix is T(I × R), and the residual matrix is E(I × JK). R is the number of principal components. PR(JK × R) can be obtained from the eigenvalue decomposition of the covariance matrix C

C=

X̅ B TX̅ B = PΛPT I−1

(2)

where Λ = diag(λ1, λ2, ..., λJK) is the eigenvalue matrix. The eigenvalues are arranged in descending order. The principal components are ordered in the same manner as eigenvalues, such that the first one describes the largest amount of variation in the data, the second one describes the second largest amount of variation, and so on. P is then separated into two parts: PR( JK × R) and PE(JK × ( JK − R)), which are called PCS and RS, respectively. T2 and SPE are two statistics for online monitoring that correspond to PCS and RS. These statistics can be calculated as follows T 2 = x̅BPR ΛR −1PR Tx̅B T

SPE = ee T ,

(3)

e = eB(1, (k − 1)J : kJ ),

eB = x̅B − x̅BPR PR T

(4)

where ΛR = diag(λ1, λ2, ..., λR), xB̅ (1 × JK) is the monitored sample, which has been unfolded on batchwise and each column of which is normalized, and e indicates the residuals at time k with the J variables. Proper control limits are defined to determine whether the process is operated under normal conditions as follows T2 ≤

R(I 2 − 1) FR ,(I − R), α I (I − R )

SPEα ≤ aχy, α 2 , 6458

a=

v 2m2 , y= 2m v

(5)

(6)

dx.doi.org/10.1021/ie403576c | Ind. Eng. Chem. Res. 2014, 53, 6457−6466


Article

where FR,(I−R),α is an F distribution with R and I − R degrees of freedom under confidence limit α and m and v are the estimated mean and variance, respectively, of the squared prediction error from the modeling data. 2.2. Fault Detection Using Variable-Wise MPCA with Time-Varying Score Covariance. Two different aspects can be identified between batchwise and variable-wise MPCA. One is the method of data set unfolding, and the other is the method of T2 calculation. (1) A combination of batchwise and variable-wise methods is used to build a two-dimensional modeling data unit. The detailed steps can be found in the Supporting Information. (2) PCA is performed on XV(IK × J). The score matrix TV(IK × R) and loading matrix PV(J × R) are obtained. TV(IK × R) is separated into K time internal score matrices Tk(I × R). Ck(R × R) is obtained by calculating the variance−covariance matrix of Tk(I × R) at each time k. x(1 × J) is the online batch data at time k. x(1 ̅ × J) is obtained by normalizing x(1 × J) using the data set XB(I,(k − 1)J:kJ). Then, the T2 statistic can be calculated as −1 T T T 2 = xP ̅ V Ck PV x̅

where J is the number of variables and pnj is the nth element in the jth loading vector. The number value represents the significance of the characteristic of the loading vector. A larger number value represents a more significant characteristic. The loading vectors with similar characteristics should be grouped in the same subspace, and Euclidean distance is used as a tool to measure similarity. Based on PD, the k-means clustering algorithm25 is employed to cluster similar denotative vectors, the corresponding loading vectors of which construct the subspace. J J-dimensional denotative vectors q1, q2, ..., qJ, exist in PD and belong to the h cluster. Meanwhile h = 2 in this article, and h cluster centers are represented by h J-dimensional vectors c1, c2, ..., ch. The sth cluster set is represented by Ds and s = 1, 2, ..., h. The subspaces are determined by the following procedure. (1) Set t = 1. Select h cluster centers c1, c2, ..., ch, randomly from PD. (2) Calculate the square of the Euclidean distance between the jth sample and the gth cluster center, d[qj,cg(t)] = ∥qj − cg(t)∥ for j = 1, 2, ..., J, and g = 1, 2, ..., h. If d[qj,cs(t)] = min{d[qj, cg(t)]}, then qj ∈ Ds(t), s = 1, 2, ..., h. (3) Recalculate the center cs(t + 1) of cluster Ds(t + 1), s = 1, 2, ..., h, as

(7)

3. MMPCA-TSBI In this study, a new method based on MMPCA-TSBI through a moving window is proposed for batch process monitoring. A combination of batchwise and variable-wise methods is used to build a two-dimensional modeling data unit to overcome the limitations on future status estimation. 3.1. Multisubspace Division Algorithm. The loading matrix P is obtained through the eigenvalue decomposition of the normalized covariance matrix. The principal component is a linear combination of variables, and different variables have different effects on the principal component. If a variable has a significant influence on a principal component, it is called a characteristic of the principal component. To present the characteristics of every principal component, a denotative matrix PD for the loading matrix is defined as pĵ = max(|pj|)

(8)

pnj̃ = |pnj | /pĵ

(9)

⎧ 0.05 ⎪ ⎪ 0.15 ⎪ ⎪ 0.25 ⎪ ⎪ 0.35 ⎪ ⎪ 0.45 ⎪ PD(n , j) = ⎨ ⎪ 0.55 ⎪ ⎪ 0.65 ⎪ ⎪ 0.75 ⎪ ⎪ 0.85 ⎪ ⎪ 0.95 ⎩

cs(t + 1) =

∑

q

q ∈ Ds (t )

(11)

where Ns is the total number of samples in the cluster Ds(t). (4) If cs(t + 1) = cs(t), s = 1, 2, ..., h, then the algorithm has converged, and the procedure is terminated. Otherwise, t = t + 1. Go back to step 2. (5) According to Ds and s = 1, 2, ..., h, a separation rule is defined as if q j ∈ Ds , then pj ∈ Ps ,

j = 1, 2, ..., J , s = 1, 2, ..., h (12)

where pj is the jth loading vector and Ps is the loading matrix in the sth subspace. h subspaces are then obtained. Notably, the number of subspaces should also be considered in the multisubspace division algorithm. In practice, it is difficult to set a unified approach for the selection of this parameter. It is suggested that this parameter can be determined by a priori knowledge and cross-validation on the training data set. In this article, we separate the original loading matrix into two subspaces. Figure 1 illustrates the division of the multisubspace algorithm. 3.2. TSBI through a Moving Window. After the division of the two principal subspaces, T2 is calculated in each subspace as follows:

0 ≤ pnj̃ < 0.1 0.1 ≤ pnj̃ < 0.2 0.2 ≤ pnj̃ < 0.3

−1 T T Ts 2 = xP (13) ̅ sΛs Ps x̅ where Ts2 is the T2 statistic in the sth principal component subspace, x(1 ̅ × J) is obtained by normalizing x(1 × J) using the data set XB(I,(k − 1)J:kJ), x(1 × J) is the online at time k, Λs = diag(λ1, λ2, ..., λRs), and Rs is the number of principal components in the sth subspace. The variances of the scores all vary with time in a batch process, 15 so that the T2 statistic is unstable along time. Thus, its control limit is calculated by the Gaussian kernel function-based kernel density estimation from the historical data at every moment in each subspace. Because the current status is affected by the past-time behavior, TSBI through a moving window is proposed for batch process monitoring. TBSI acts as a moving weighting strategy.

0.3 ≤ pnj̃ < 0.4 0.4 ≤ pnj̃ < 0.5 0.5 ≤ pnj̃ < 0.6 0.6 ≤ pnj̃ < 0.7 0.7 ≤ pnj̃ < 0.8 0.8 ≤ pnj̃ < 0.9 0.9 ≤ pnj̃ ≤ 1.0

1 Ns

(10) 6459



Article

where L is the window length. The control limit for PTs is 1 − α, which denotes that some abnormal event occurs when the value of PTs in the sth subspace exceeds 1 − α. Otherwise, the process is considered normal. The proposed MMPCA-TSBI approach for online batch process monitoring is illustrated in Figure 2, and its detailed step-by-step procedure is as follows:

Figure 1. Illustration of division of multisubspaces.

The influence of noise on monitoring performance can also be mitigated. The time-window size is a parameter of the moving weighting strategy. A large time-window size results in an extremely high influence of past-time behavior on current behavior, which indicates that the monitoring statistic will exhibit poor sensitivity to changes in the process. This condition will degrade monitoring performance when an abnormal event occurs. By contrast, if the time-window size is small, large noise in the industrial process will increase the false alarm rate. Therefore, an appropriate time-window size should be selected to simultaneously obtain proper past information and mitigate the influence of noise. In this method, 7−12 time intervals are proposed for the time-window size according to different processes. The subsequent Bayesian inference is similar to that described in ref 23, which can be consulted for more details. The fault probability of online data in the time window of each subspace can be calculated as b(F | Xs ,l) =

Figure 2. Schematic diagram of MMPCA-TSBI for batch process monitoring.

(A) Off-line Modeling (1) Unfold X(I × J × K) to XB(I × JK). (2) Obtain X̅ B(I × JK) by normalizing the data XB(I × JK) in each column. (3) Rearrange the scaled X̅ B(I × JK) to X̅ V(IK × J). (4) Apply the division of multisubspace algorithm to X̅ V(IK × J). (5) Calculate the control limit of Ti2 in each subspace. (6) Specify the time-window size of the moving weighting strategy. (B) Online Monitoring (1) Normalize the new batch sampling data at time k, x(1 × J), using the same moment data XB(I,(k − 1)J:kJ). (2) Project the data into the multisubspace. (3) Calculate Ts2 in the sth subspace. (4) Calculate PTs by TSBI in the current time window in the sth subspace. (5) Monitor whether the value of PTs in the sth subspace exceeds the confidence level (CL).

b(Xs ,l | F)b(F) b(Xs ,l)

(14)

b(Xs ,l) = b(Xs ,l | N)b(N) + b(Xs ,l | F)b(F)

(15)

where Xs,l represents the lth data point in the time window in the sth subspace and N and F represent normal and abnormal conditions, respectively. b(N) and b(F) can be denoted as α and 1 − α, respectively, where α is the confidence level. The conditional probabilities b(Xs,l|N) and b(Xs,l|F) are defined as ⎛ T 2 ⎞ s,l ⎟, b(Xs ,l | N) = exp⎜⎜ − 2⎟ ⎝ Ts , l ,lim ⎠

2⎞ ⎛ T s , l ,lim ⎟ b(Xs ,l | F) = exp⎜⎜ − 2 ⎟ ⎝ Ts , l ⎠

(16)

where Ts,l,lim2 is the control limit of the lth data point in the time window in the sth subspace. The b(F|Xs,l) is employed as a weighting coefficient. The final probabilistic statistic in the sth subspace can then be determined as L

PTs =

∑ l=1

4. APPLICATION EXAMPLES In this section, the proposed batch process monitoring method is tested using two examples: a simple numerical process and a fed-batch penicillin fermentation batch process. The false alarm rate, fault detection rate, and missed alarm rate are three

b(F|Xs ,l)2 L

∑l = 1 b(F|Xs ,l)

(17) 6460



Article

⎡ 0.193 0.689 ⎤ ⎡ 0.811 −0.226 ⎤ u (i ) = ⎢ ⎥w ⎥u(i − 1) + ⎢ ⎣ 0.477 0.415 ⎦ ⎣−0.320 − 0.749 ⎦

commonly used performance indexes for comparing monitoring methods. Before the fault occurs, the batch is in a period of normal operation. If one statistic crosses the control limit, it is called a false alarm. After the fault occurs, the batch is in a period of fault operation. If one statistic crosses the control limit, it is called a fault detection, or if the statistic does not cross the control limit, it is called a missed alarm. The false alarm rate, fault detection rate, and missed alarm rate are calculated as follows

(i − 1)

The input w is a random vector, each element of which is uniformly distributed over the interval (−2, 2). The output y is equal to z plus a random noise vector v. The Gaussian signal v has a zero mean and a variance of 0.1. Both input u and output y are measured, but z and w are not. Under the normal operating conditions, X(100 × 5 × 800) can be obtained from 100 cycles of the continuous process, where 100 represents the number of batches, 5 denotes the number of measurement variables (y1, y2, y3, u1, u2), and 800 is the number of sampling instants within a batch. X(100 × 5 × 800) typical batch process data serve as benchmark batches. To evaluate the performance of the proposed method, two simulated faults for monitoring batches are created. The corresponding process parameters generated for each case are listed in Table 1.

false alarm rate = [(∑ false alarm) /(total number of samples during normal operation)] (18)

× 100%

fault detection rate = [(∑ false detection) /(total number of samples during fault operation)] (19)

× 100% missed alarm rate = [(∑ missed alarm)

Table 1. Settings of Operating Conditions in the Numerical Process

/(total number of samples during fault operation)] (20)

× 100%

The fault detection rate and missed alarm rate are two complementary statistics, related by missed alarm rate = 100 − fault detection rate

case no.

type

size

duration

1 2

mean step change in w1 mean step change in w2

+2.5 +1.5

200−800 200−800

In MMPCA-TSBI, the two-dimensional modeling data of the numerical process are obtained through a combination of batchwise and variable-wise methods. The division of the multisubspace algorithm is applied, and the gray schematic diagram of denotative matrix PD of the numerical process is shown in Figure 3. This matrix presents all of the characteristics

(21)

so only one of them is used to compare the monitoring performance. To synthesize the false alarm rate and fault missed alarm rate, the error rate is calculated as error rate = [(∑ false alarm +

(25)

∑ missed alarm)

/(total number of samples during overall operation)] × 100%

(22)

In this article, the false alarm rate and missed alarm rate are used for the numerical process. The false alarm rate, missed alarm rate, and error rate are used for the fed-batch penicillin fermentation process. 4.1. Case Study on a Numerical Process. In this section, the proposed batch process monitoring method, MMPCATSBI, and the variable-wise MPCA method15 are applied to a simple numerical multivariate process that is initially a typical continuous process. To compare the monitoring effects of the two methods, the numerical process is modified to generate a certain number of cycles for batch process data. Zhao et al. generated numerical batch data using this approach.26 The simple numerical multivariate process suggested by Ku et al.27 and modified by Lee et al.28 is given by ⎡ 0.018 −0.191 0.287 ⎤ ⎢ ⎥ z(i) = ⎢ 0.847 0.264 0.943 ⎥z(i − 1) ⎢⎣−0.333 0.514 − 0.217 ⎥⎦ ⎡1 2 ⎤ ⎢ ⎥ + ⎢ 3 −4 ⎥u(i − 1) ⎣− 2 1 ⎦

(23)

y(i) = z(i) + v(i)

(24)

Figure 3. Gray schematic diagram of denotative matrix in the numerical process.

of every loading vector. The x coordinate is the sequence number of the loading vectors, whereas the y coordinate is the sequence number of variables in the numerical process. White is the most significant characteristic of the loading vector, and black is the least significant. Table 2 shows that the original loading space is separated into P1 and P2. The corresponding first, third, and fifth loading vectors combine to form a

where u is the correlated input 6461



Article

Table 2. Result of Division of the Multisubspace Algorithm in the Numerical Process subspace

loading vectors

P1 P2

1, 3, 5 2, 4

subspace. The second and fourth loading vectors combine to form another subspace. The proposed method was compared with the variable-wise MPCA approach. Both control limits of MMPCA-TSBI and variable-wise MPCA statistics are based on a significance level of 99%. The first three principal components are retained in variable-wise MPCA, and the time-window size of the proposed method is 10 time intervals. The monitoring results in terms of false alarm rate and missed alarm rate are reported in Table 3. Table 3. Comparison of the Results of the Numerical Process Using Variable-Wise MPCA and MMPCA-TSBI monitoring method

monitoring statistic

false alarm rate (%)

missed alarm rate (%)

Figure 4. Fault detection results of variable-wise MPCA for case 1 in the numerical process.

Case No. 1 variable-wise MPCA MMPCA-TSBI

T2 SPE PT1 PT2

1.0 0.5 0 10.6

94.3 92.5 12.5 77.4

0 0.5 0 5.0

77.7 75.0 4.2 15.8

Case No. 2 variable-wise MPCA MMPCA-TSBI

T2 SPE PT1 PT2

For the two test cases, we can readily observe that the false alarm rate of the proposed method is lower than that of the variable-wise MPCA method. Meanwhile, the proposed method yields a significantly lower missed alarm rate than variable-wise MPCA. From Figure 3, it can be clearly seen that the first two loading vectors capture all of the variables. The first three loading vectors are used in variable-wise MPCA, so it can capture all of the variables. For the first test case, the time-series plots of the monitoring statistics of the regular variable-wise MPCA and MMPCA-TSBI are shown in Figures 4 and 5, respectively. The monitoring performance of the proposed method is much better than that of variable-wise MPCA. 4.2. Case Study on Fed-Batch Penicillin Fermentation Batch Process. 4.2.1. Fed-Batch Penicillin Fermentation Process Description. In this section, the proposed monitoring method MMPCA-TSBI is applied to the monitoring problems of a well-known benchmark process, the fed-batch penicillin fermentation industrial process,29 to investigate its effectiveness. A flow diagram of the penicillin fermentation process is given in Figure S3 of the Supporting Information. A modular simulator (PenSim v2.0 was used in this work) was developed by the monitoring and control group of the Illinois Institute of Technology (http://mypages.iit.edu/ ∼cinar) for fed-batch fermentation in 2002. The substrate feed rate during the fed-batch stage is under open-loop operation, whereas the pH and temperature in the bioreactor are controlled using two cascade loops that manipulate the acid/base equilibrium and the hot/cold water flow ratio. Furthermore, air is continuously fed into the bioreactor to maintain a desirable oxygen level for cell growth.

Figure 5. Fault detection results of MMPCA-TSBI for case 1 in the numerical process.

In this study, the observation data on 16 measurement variables given in Table S1 of the Supporting Information were collected for batch process monitoring. Initial conditions and set points are listed in Table S2 of the Supporting Information. The sampling time used was 0.5 h, and the overall duration of each batch was 400 h, including both the batch and fed-batch stages. The normal measurement variables in the fermentation process are shown in Figure S4 of the Supporting Information. The study included 100 normal batches and 12 fault batches, which are listed in Table 4. 4.2.2. Fed-Batch Penicillin Fermentation Process Monitoring Results and Comparison. The two-dimensional modeling data of the fed-batch penicillin fermentation process were obtained through a combination of batchwise and variable-wise methods. The division of the multisubspace algorithm is applied, and a gray schematic diagram of denotative matrix PD 6462



Article

fed-batch penicillin fermentation process. The control limits of batchwise MPCA, variable-wise MPCA, and MMPCA-TSBI statistics are based on a significance level of 99%. The first eight principal components are retained in batchwise MPCA and variable-wise MPCA. The time-window size of the proposed method is 10 time intervals. For batchwise MPCA, the current deviation was used to fill the missing future observations. The monitoring results in terms of false alarm rate and missed alarm rate are listed in Tables 6 and 7, respectively. In all cases, the

Table 4. Test Cases of the Fed-Batch Penicillin Fermentation Process case no. 1 2 3 4 5 6 7 8 9 10 11 12

variable name aeration rate aeration rate aeration rate aeration rate agitator power agitator power agitator power agitator power substrate feed rate substrate feed rate substrate feed rate substrate feed rate

variable no.

type

size (%)

1 1 1 1 2 2 2 2 3

step step drift drift step step drift drift step

−2 −1.8 −0.5 −0.3 −1.8 −1.7 −0.7 −0.6 −3

fault duration (h) 100−400 100−400 150−400 150−400 100−400 100−400 150−400 150−400 200−400

Table 6. Comparison Results of the False Alarm Rates Using Batchwise MPCA, Variable-Wise MPCA, and MMPCA-TSBI batchwise MPCA

3

step

−2.8

200−400

3

drift

−0.015

200−400

3

drift

−0.01

200−400

of the fed-batch penicillin fermentation process is shown in Figure 6, which presents all of the characteristics of every

2

case no.

T

1 2 3 4 5 6 7 8 9 10 11 12

0 0 0 0 0 0 0 0 0 0 0 0

variable-wise MPCA 2

SPE

T

5.5 0 12.7 14.7 1.0 1.0 1.0 11.0 7.7 4.5 1.0 7.8

0 0 0.3 4.0 0 1.0 0 0 0 0 0 0

MMPCA-TSBI

SPE

PT1

PT2

0.5 0 3.0 1.0 0 1.0 0 2.3 3.0 4.5 0.5 0

0 0 12.3 0 0 0 0.3 0 0 5.3 0 0

0 0 0 0 0 0 0 0 0.7 0 0 0

Table 7. Comparison Results of the Missed Alarm Rates Using Batchwise MPCA, Variable-Wise MPCA, and MMPCA-TSBI batchwise MPCA

Figure 6. Gray schematic diagram of denotative matrix in the fedbatch penicillin fermentation process.

loading vector. The x coordinate is the sequence number of loading vectors, whereas the y coordinate is the sequence number of variables in the fed-batch penicillin fermentation process. White is the most significant characteristic of the loading vector, and black is the least significant. Table 5 shows that the original loading space is separated into P1 and P2. The corresponding 1st, 2nd, and 10th−16th loading vectors combine to form one subspace, whereas the 3rd−9th loading vectors combine to form the other subspace. The proposed MMPCA-TSBI approach was compared to the batchwise MPCA6 and variable-wise MPCA15 methods in the

loading vector no.

P1 P2

1, 2, 10, 12, 13, 14, 15, 16 3, 4, 5, 6, 7, 8, 9

2

case no.

T

SPE

T

1 2 3 4 5 6 7 8 9 10 11 12

100 100 100 100 100 100 100 100 49.9 100 14.5 18.5

14.0 19.8 10.0 14.8 13.1 19.1 19.2 12.8 70.3 68.6 3.3 5.0

29.0 42.9 17.6 27.8 33.8 37.4 21.4 19.0 66.6 79.8 8.5 15.7

MMPCA-TSBI

SPE

PT1

PT2

85.0 85.9 21.8 40.5 95.0 87.5 71.3 71.1 78.3 89.0 15.7 25.4

100 100 99.4 99.8 100 99.3 99.8 99.8 0.3 13.7 5.8 3.7

5.7 11.6 11.4 17.7 8.5 15.1 14.6 15.8 100 97.3 30.7 54.9

missed alarm rate of the proposed method is better than that of variable-wise MPCA. In most cases (i.e., cases 1, 2, 5−7, 9, 10, and 12), the missed alarm rate of the proposed method is slightly better than that of batchwise MPCA. In some drift faults (i.e., cases 3, 4, 8, and 11), the missed alarm rate of batchwise MPCA is slightly better than that of the proposed method, but the false alarm rate of batchwise MPCA is slightly higher than that of the proposed method. To synthesize the false alarm rate and missed alarm rate, the monitoring results in terms of error rate are listed in Table 8. In all cases except case 11, the error rate of the proposed method is lowest. Meanwhile, the first time instants to detect the fault are listed in Table 9. In this table, a dash indicates that the fault could not be detected. In all cases except cases 3 and 11, the proposed method is the

Table 5. Results of the Division of the Multisubspace Algorithm in the Fed-Batch Penicillin Fermentation Process subspace

2

variable-wise MPCA

6463



Article

Table 8. Comparison Results of the Error Rates Using Batchwise MPCA, Variable-Wise MPCA, and MMPCA-TSBI batchwise MPCA

variable-wise MPCA

MMPCA-TSBI

case no.

T2

SPE

T2

SPE

PT1

PT2

1 2 3 4 5 6 7 8 9 10 11 12

75.0 75.0 62.5 62.5 75.0 75.0 62.5 62.5 24.6 50.0 7.2 9.2

11.9 14.9 11.0 14.7 10.1 13.9 12.3 12.1 39.0 36.5 2.1 6.4

21.7 32.2 11.2 18.8 25.3 28.3 13.3 11.9 33.3 39.9 4.2 7.9

63.9 64.4 14.7 25.7 71.3 65.9 44.5 45.3 40.7 46.8 7.4 12.8

75.0 75.0 66.8 62.4 75.0 74.5 62.5 62.4 0.1 9.5 2.8 1.9

4.2 8.7 7.1 11.1 6.7 11.3 9.1 9.9 50.4 48.6 19.8 27.4

Table 10. Comparison Average Results of False Alarm Rates, Missed Alarm Rates, and Error Rates Using Batchwise MPCA, variable-Wise MPCA, and MMPCA-TSBI

variable no. 1

2

3

Table 9. Comparison Results of the First Times to Detect the Fault Using Batchwise MPCA, Variable-Wise MPCA, and MMPCA-TSBI batchwise MPCA

variable-wise MPCA

false alarm rate missed alarm rate error rate false alarm rate missed alarm rate error rate false alarm rate missed alarm rate error rate

batchwise MPCA

variable-wise MPCA

T2

SPE

T2

SPE

PT1

PT2

0

8.2

1.1

1.1

3.1

0

100

14.6

29.3

58.3

99.6

11.6

13.1 4.7

21.0 0.2

42.2 0.8

69.8 0.1

7.8 0

15.8

27.9

81.2

99.7

13.5

68.8 0

12.1 5.3

19.7 0

56.8 2.0

68.6 1.3

9.3 0.2

45.7

36.8

42.6

51.8

5.9

73.0

22.8

21

21.3

26.9

3.6

36.6

68.8 0 100

MMPCATSBI

MMPCA-TSBI

case no.

T2

SPE

T2

SPE

PT1

PT2

1 2 3 4 5 6 7 8 9 10 11 12

− − − − − − − − 258 − 229 237

100 100 173 175 100 100 190 191 203 203 209 210

100 100 175 180 100 100 192 194 216 208 217 228

107 107 167 184 109 104 220 197 211 203 229 213

− − 399 − − 152 − − 201 201 210 208

100 100 174 175 100 100 187 191 − 232 280 310

first to detect the fault. Because the multisubspace division is according to variables, the average monitoring results of the same variable with fault are listed in Table 10. It can be clearly seen that the proposed method provides an obvious improvement in error rate if the fault occurs in variable 3 (substrate feed rate), and is provides a slight improvement in error rate if the fault occurs in variable 1 (aeration rate) or variable 2 (agitator power). For case 9, in the first 200 h, the process operates normally. Then, a step decrease error in the substrate feed rate with a size of 3% arises until the end of the batch. The time-series plots of the monitoring statistics of the three methods are shown in Figures 7−9. In this case, both batchwise MPCA and variablewise MPCA perform poorly, because both of these methods cause information loss by discarding the last several principal components. In contrast, PT1 of MMPCA-TSBI can detect the fault because the principal components in the PCS of the variable-wise MPCA method are the 1st−8th, whereas those in the 1st subspace of the proposed method are the 1st, 2nd, and 10th−16th, as shown in Table 5. The principal components in the 1st subspace can detect most faults. Thus, the monitoring performance of the proposed method is better than the performances of both batchwise MPCA and variable-wise MPCA. The use of time-series Bayesian inference through a

Figure 7. Fault detection results of batchwise MPCA for case 9 in the fed-batch penicillin fermentation process.

moving window in the proposed method can reduce the impact of noise, and it also reflects the influence of past-time behavior on the current behavior in the process.

5. CONCLUSIONS This article presents a new method for batch process monitoring called MMPCA-TSBI. First, a combination of batchwise and variable-wise methods is used to build a twodimensional modeling data unit. The division of the multisubspace algorithm is then applied. Finally, Bayesian inference along a time series through a moving time window is applied to detect abnormal events in the batch process. The proposed method does not require future status prediction and presents a division of the multisubspace algorithm that can retain all of the information in the original principal component space and group the principal components with similar characteristics. The proposed approach overcomes noise and employs 6464



Article

process. This information is available free of charge via the Internet at http://pubs.acs.org/.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: 0086-21-64251036. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge the support of the following foundations: 973 Project of China (2013CB733605); National Natural Science Foundation of China (21176073); Fundamental Research Funds for the Central Universities; and Open Research Project of the State Key Laboratory of Industrial Control Technology, Zhejiang University, Hangzhou, Zhejiang, China (No. ICT1429).

■

Figure 8. Fault detection results of variable-wise MPCA for case 9 in the fed-batch penicillin fermentation process.

Figure 9. Fault detection results of MMPCA-TSBI for case 9 in the fed-batch penicillin fermentation process.

information on the self-correlation in industrial processes by using TSBI through a moving window. The proposed method, MMPCA-TSBI, was applied to both a numerical process and the fed-batch penicillin fermentation process. The monitoring results indicate that the proposed method can perform well. Therefore, the proposed method should prove useful for batch process monitoring. Future work could focus on developing a strategy for determining the number of principal component subspaces and on making the proposed method suitable for nonlinear and non-Gaussian processes.

■

REFERENCES

(1) Bakshi, B. R.; Locher, G.; Stephanopoulos, G.; Stephanopoulos, G. Analysis of operating data for evaluation, diagnosis and control of batch operations. J. Process Control 1994, 4, 179. (2) Lennox, B.; Montague, G. A.; Hiden, H. G.; Kornfeld, G.; Goulding, P. R. Process monitoring of an industrial fed-batch fermentation. Biotechnol. Bioeng. 2001, 74, 125. (3) Rashid, M. M.; Yu, J. Nonlinear and Non-Gaussian Dynamic Batch Process Monitoring Using a New Multiway Kernel Independent Component Analysis and Multidimensional Mutual Information Based Dissimilarity Approach. Ind. Eng. Chem. Res. 2012, 51 (33), 10910. (4) Kresta, J. V.; Marlin, T. E.; MacGregor, J. F. Development of inferential process models using PLS. Comput. Chem. Eng. 1994, 18, 597. (5) Wold, S.; Geladi, P.; Esbensen, K. Multi-way principal components and PLS analysis. J. Chemom. 1987, 1, 41. (6) Nomikos, P.; MacGregor, J. F. Monitoring of batch processes using multi-way principal component analysis. AIChE J. 1994, 40, 1361. (7) Nomikos, P.; MacGregor, J. F. Multivariate SPC charts for monitoring batch processes. Technometrics 1995, 37 (3), 41. (8) Nomikos, P.; MacGregor, J. F. Multiway partial least squares in monitoring batch processes. Chemom. Intell. Lab. Syst. 1995, 30, 97. (9) Cho, H. W.; Kim, K. J. A method for predicting future observations in the monitoring of a batch process. J. Qual. Technol. 2003, 35 (1), 59. (10) Louwerse, D. J.; Smilde, A. K. Multivariate statistical process control of batch processes based on three-way models. Chem. Eng. Sci. 2002, 55 (7), 1225. (11) Ü ndey, C.; Tatara, E.; Cinar, A. Intelligent realtime performance monitoring and quality prediction for batch/fed-batch cultivations. J. Biotechnol. 2004, 108 (1), 61. (12) Rännar, S.; MacGregor, J. F.; Wold, S. Adaptive batch monitoring using hierarchical PCA. Chemom. Intell. Lab. Syst. 1998, 41 (1), 73. (13) Wold, S.; Kettaneh, N.; Friden, H.; Holmberg, A. Modelling and diagnostics of batch processes and analogous kinetic experiments. Chemom. Intell. Lab. Syst. 1998, 44 (1−2), 331. (14) Ü ndey, C.; Ertunc, S.; Cinar, A. Online batch/fed-batch process performance monitoring, quality prediction, and variable-contribution analysis for diagnosis. Ind. Eng. Chem. Res. 2003, 42 (20), 4645. (15) Lee, J. M.; Yoo, C. K.; Lee, I. B. Enhanced process monitoring of fed-batch penicillin cultivation using time-varying and multivariate statistical analysis. J. Biotechnol. 2004, 110 (2), 119. (16) Albazzaz, H.; Wang, X. Z. Statistical process control charts for batch operations based on independent component analysis. Ind. Eng. Chem. Res. 2004, 43 (21), 6731.

ASSOCIATED CONTENT

S Supporting Information *

Illustrations of a batchwise data set of a batch process and batchwise unfolding transformed to variable-wise unfolding and detailed information on the fed-batch penicillin fermentation 6465



Article

(17) Lee, J. M.; Yoo, C. K.; Lee, I. B. On-line batch process monitoring using different unfolding method and independent component analysis. J. Chem. Eng. Jpn. 2003, 36 (11), 1384. (18) Valle, S.; Li, W.; Qin, S. J. Selection of the Number of Principal Components: The Variance of the Reconstruction Error Criterion with a Comparison to Other Methods. Ind. Eng. Chem. Res. 1999, 38, 4389. (19) Zwick, W. R.; Velicer, W. F. Comparison of five rules for determining the number of components to retain. Psych. Bull. 1986, 99, 432. (20) Li, Y.; Tang, X. C. Improved performance of fault detection based on selection of the optimal number of principal components. Acta Autom. Sin. 2009, 35, 1550. (21) Jolliffe, I. T. A note on the use of principal components in regression. Appl. Stat. 1982, 300. (22) Jiang, Q.; Yan, X.; Zhao, W. Fault Detection and Diagnosis in Chemical Processes Using Sensitive Principal Component Analysis. Ind. Eng. Chem. Res. 2013, 52 (4), 1635. (23) Ge, Z.; Zhang, M.; Song, Z. Nonlinear process monitoring based on linear subspace and Bayesian inference. J. Process Control 2010, 20, 676. (24) Tong, C.; Song, Y.; Yan, X. Distributed Statistical Process Monitoring Based on Four-Subspace Construction and Bayesian Inference. Ind. Eng. Chem. Res. 2013, 52 (29), 9897. (25) Hartigan, J. A.; Wong, M. A. Algorithm AS 136: A K-Means Clustering Algorithm. J. R. Stat. Soc., Ser. C (Appl. Stat.) 1979, 28 (1), 100. (26) Zhao, C.; Wang, F.; Jia, M. Dissimilarity analysis based batch process monitoring using moving windows. AIChE J. 2007, 53 (5), 1267. (27) Ku, W.; Storer, R. H.; Georgakis, C. Disturbance detection and isolation by dynamic principal component analysis. Chemom. Intell. Lab. Syst. 1995, 30, 179. (28) Lee, J. M.; Yoo, C. K.; Lee, I. B. Statistical process monitoring with independent component analysis. J. Process Control 2004, 14, 467. (29) Birol, G.; Ü ndey, C.; Cinar, A. A modular simulation package for fed-batch fermentation: Penicillin production. Comput. Chem. Eng. 2002, 26, 1553.

6466


Batch Process Monitoring Based on Multisubspace Multiway Principal

Recommend Documents