Article pubs.acs.org/IECR
Dynamic Multimode Process Modeling and Monitoring Using Adaptive Gaussian Mixture Models Xiang Xie† and Hongbo Shi*,† †
Key Laboratory of Advanced Control and Optimization for Chemical Processes (East China University of Science and Technology), Ministry of Education, Shanghai 200237, China ABSTRACT: For multimode processes, it is inevitable to encounter disturbances, such as equipment aging, catalyst deactivation, sensor drifting, reaction kinetics drifting, or adding new operating modes. The existing monitoring algorithms are established either for coping with multimode feature under time-invariant circumstance or for handling the time-varying problem of processes with single operating mode. The purpose of this article is to develop an effective modeling and monitoring approach for complex processes with both multimode and time-varying properties. We propose a novel adaptive monitoring scheme based on Gaussian Mixture Model (GMM). The new method is able to model different operating modes as well as trace process variations. The effectiveness and efficiency of the new method are validated by a numerical example and the Tennessee Eastman (TE) simulation platform in different scenarios.
1. INTRODUCTION Multimode is one of the most significant features of modern industrial processes due to the various demands of the market. To ensure the safety of production, it is meaningful to develop effective monitoring methods for multimode processes. Once abnormal events occur during the production, faulty symptoms should be reflected in the monitoring chart and faulty alarms should be triggered. Recently, this area has been intensively studied, and many statistical monitoring algorithms have been reported. For processes with multiple operating modes, the most intuitive idea is to build separate models for different modes, and this concept is adopted by most works published on this field. Since PCA (Principal Component Analysis) and PLS (Partial Least Squares) are the most mature statistical techniques in single-mode process monitoring, they are extended to multimode processes in a multiple-model way, combined with other pattern classification algorithms, and the traditional T2 and SPE statistics are used as monitoring indices.1−8 In contrast with these PCA/PLS based methods, GMM is extended to multimode process monitoring recently. GMM is one of the artificial intelligence tools for pattern recognition and is widely used in background modeling and subtraction. Choi et al. integrated GMM with PCA and DA (Discriminant Analysis) to monitor processes with nonlinearity and multimode properties.9 Along the same line, Yoo et al. used GMM combined with multiway PCA to monitor a biological batch process.10 However, these multimode methods mentioned above utilized only one local model (i.e., Gaussian component) for each online sample and neglected information from other models, which may lead to biased monitoring results. To overcome this deficiency, Yu et al. proposed a Bayesian Inference-based GMM method. In their approach, a probability monitoring index is constructed by incorporating all the Gaussians’ posterior probabilities with their Mahalanobis distances.11,12 Later, Chen extended the application of GMM to monitoring a batch semiconductor etch process.13 Yu © 2012 American Chemical Society
developed the nonlinear version of GMM for a multimode wastewater treatment process.14 Nevertheless, these multimode monitoring models are time-invariant, where the Gaussians stay static after the offline training procedure. If normal process changes occur, such as catalyst deactivation, equipment aging, sensor and reaction kinetics drifting, etc., the performance of fixed monitoring model will degrade and false alarms will be triggered in an unexpected way (Type II error). Therefore, it is meaningful to develop adaptive monitoring models to accommodate these time-varying behaviors of processes. In essence, the dynamic characteristics of an industrial process include changes in the mean and the covariance of collected data. To make the monitoring model up-to-date, Wold developed an exponentially weighted moving average (EWMA) scheme for both PCA and PLS, in which the model is dynamically updated by giving more weight to recent samples than to past samples.15 Wise et al. applied an exponentially weighted moving covariance (EWMC) method to monitor a microelectronics manufacturing process.16 Li et al. proposed two adaptive PCA algorithms for both sample-wise and blockwise recursions,17 and these algorithms have successful applications in adaptive process monitoring.18−21 However, the data used in recursive PCA is ever-growing, leading to a reduction in the speed of adaptation as the data size increases. In addition, the older samples are not representative of the current process status. To overcome the aforementioned deficiencies, Wang et al. proposed a fast moving window PCA algorithm.22 As the window slides along the data, a new process model is updated by including the latest sample and excluding the oldest one. Their method has also been extended into nonlinear filed by Liu et al., who developed a moving window kernel PCA approach.23 More recently, Jeng took Received: Revised: Accepted: Published: 5497
November 23, 2011 February 24, 2012 March 29, 2012 March 29, 2012 dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
advantage of both recursive PCA and moving window PCA and proposed a combined online updating and monitoring scheme.24 However, most of the existing adaptive monitoring methods are developed explicitly for single mode processes, and they cannot be used directly under multimode circumstances. The main contribution of this article is to build an effective statistical monitoring model for processes with both multimode and time-varying properties. The proposed approach, the socalled “moving window GMM”, is an extension of Yu et al. multimode approach under time-varying situations. The novel method is able to describe separate multimode features of processes as well as trace process changes and keep the monitoring model up-to-date. In our approach, the number of initial Gaussian components and their parameters are calculated automatically through an offline training procedure. When used online, the posterior probability vector for each sample plays the role of model indicator and decides dynamically which Gaussian component should get updated and which stay unchanged. The mean vector and covariance matrix of the chosen Gaussian are then updated online through a moving window. The rest of the paper is organized as follows. Section 2 introduces the original GMM based monitoring scheme; section 3 elaborates upon the novel adaptive GMM modeling and monitoring approach; then the effectiveness and efficiency of the proposed method are verified by a numerical example and the well-known TE (Tennessee Eastman) simulation platform in section 4 and 5, respectively; finally, section 6 gives the conclusion of the whole paper.
k=1
(2)
n
α(ks + 1)
∑ j = 1 β(s)(CK |xj)
=
(3)
n n
μ(ks + 1)
∑ j = 1 β(s)(CK |xj)xj
=
n
∑ j = 1 β(s)(CK |xj) n
Σ(ks + 1)
(4)
(
=
T
)
∑ j = 1 β(s)(CK |xj)(xj − μ(ks + 1)) xj − μ(ks + 1) n
∑ j = 1 β(s)(CK |xj)
(5)
where β (Ck|xj) denotes the posterior probability of the j-th training sample within the k-th Gaussian component Ck at the sth iteration. g(·) is the Gaussian function. α(s+1) , μ(s+1) , and k k Σ(s+1) are prior probability and the mean and covariance of the k k-th Gaussian component at the (s+1)-th iteration, respectively. By repeating E-step and M-step iteratively, the posterior probabilities and corresponding distribution parameters are computed with the specified accuracy.25 When GMM is utilized in monitoring procedure, different from the traditional T2 and SPE statistics used in single model procedure, a probability monitoring index BIP is proposed by Yu and Qin as follows11 (s)
K
BIP =
∑ β(Ck|xt )PL(k)(xt ) k=1
(6)
where β(Ck|xt) donates the posterior probability of the online sample xt to each Gaussian component Ck, and it has the same computational procedure with the E-step in eq 2. P(k) L (xt) is the local Mahalonobis distance-based probability index, which is PL(k)(xt ) = Pr{D(x , Ck)|x ∈ Ck ≤ D((xt , Ck)|x ∈ Ck)} (7)
where D(xt,Ck) donates the Mahalanobis distance of xt to Ck. As 0 ≤ P(k) L (xt) ≤ 1, we have K
0 ≤ BIP ≤
∑ β(Ck|xt ) = 1 k=1
(8)
The process faulty alarms are triggered when BIP > 0.95 or 0.99. It is noticed that the faulty symptom with BIP is not as evident as distance based indices (e.g., T2 or SPE). Moreover, the subsequent diagnostic procedure is difficult to be carried out based on probability index. In this article, we directly utilize the Mahalanobis distance as local index, which is
K
∑ αkg(x|μk , Σk)
K
∑i = 1 α(i s)g (xj|μ(ks), Σ(ks))
M-step:
2. STATIC GMM BASED PROCESS MONITORING For obtained data from a multimode process, the assumption of global multivariate Gaussian distribution is invalid due to their mean and covariance shift. However, from another point of view, the local Gaussian distribution is still appropriate to characterize each subset of observing data from the same operating modes. Therefore, the GMM is suited to describe the data generated from different operating modes. Specifically, each Gaussian component represents an individual operating mode, and the number of Gaussians represents the number of all possible modes through the entire monitoring process. The GMM algorithm adopted here is originally proposed by Figueiredo and A. K. Jain25 and is extended in the statistical process monitoring area by Yu and Qin.10 Let x∈Rm be an m-dimensional sample from a multimode process. Its probability density function can be formulated as GMM below p(x|μ , Σ) =
α(ks)g (xj|μ(ks), Σ(ks))
β(s)(Ck|xj) =
(1)
where K is the number of Gaussian components included in GMM; αk denotes the mixing weight of the k-th component and adds up to one. {μ1,Σ1,···,μK,ΣK} represents the mean vector and covariance matrix of each local Gaussian model. To estimate the maximum likelihood distribution parameters of GMM, the EM (Expectation-Maximization) algorithm is adopted. Given the training data {x1,···,xn} and an initial estimate {α(0),μ(0),Σ(0)}, the E-step (expectation-step) and M-step (maximization-step) are performed iteratively as follows E-step:
(k) Dlocal (xt ) = (xt − μk)T Σ−k 1(xt − μk)
(9)
Then the integrated global monitoring statistics, donated as Dglobal, are presented in the following form K
Dglobal =
(k) (x t ) ∑ β(Ck|xt )Dlocal k=1
(10)
The upper control limit for Dglobal, denoted as DL, can be calculated from the F-distribution, which is1 5498
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research DL =
p(n2 − 1) Fp , n − p; γ n(n − p)
Article
The probability density function and its corresponding posterior probabilities for each online sample are also updated with the latest means and covariance matrix as follows
(11)
where n and p represent the samples and variables’ number of historical data respectively. Fp,n−p;γ is an F distribution with p and n-p degrees of freedom under given significant level γ, usually 95% or 99%. If Dglobal > DL for several consecutive samples, the process is considered to be abnormal, and the faulty alarms are triggered.
gω+ 1(xω+ 1|Ck) =
3. ADAPTIVE GMM BASED PROCESS MONITORING 3.1. The Moving Window GMM. When slow and natural process changes occur in a multimode process, it is reasonable to update its corresponding monitoring component rather than the whole GMMs. Therefore, the main problem is changed into finding out the most suitable Gaussian component for online samples. Luckily, with the help of posterior probabilities obtained with eq 2, the Gaussian model with the maximum value of β(Ck|xt) can be used as the target. Then the target Gaussian is updated through a moving window, since the old data cannot represent the current status of process any more. The latest sample is augmented to the monitored data matrix, while the oldest sample is discarded, in order to keep a fixed capacity for monitored samples, i.e. the window size N. Let the ω-th monitored data matrix with window size N be Xω = [xω−N+1 xω−N+2 ··· xω]T, and the next monitored data matrix would be Xω+1 = [xω−N+2 xω−N+3 ··· xω+1]T. These two data matrices share a transient data matrix of X̃ = [xω−N+2 xω−N+3 ··· xω]T. A two-step adaptation is used to recursively update the covariance matrix. First, with the elimination of the oldest sample from Xω, the mean vector and covariance matrix associated with X̃ can be represented in terms of N 1 xω− N + 1 μ̃ = μω − (12) N−1 N−1 N − 1 ⎛⎜ 1 Σω − Σ̃ = (X ω− N + 1 − μω− N + 1)T ⎝ N−2 N−1 ⎞ (X ω− N + 1 − μω− N + 1) + Δμ̃ TωΔμ̃ω⎟ ⎠
βω+ 1(Ck|xω+ 1) =
N−1 1 μ̃ω + xω+ 1 N N
Σω+ 1 =
N−2 ̃ 1 (xω+ 1 − μω+ 1)T Σω + N−1 N−1
α(ks + 1)
(13)
(15)
1 (xω− N + 1 − μω− N + 1)T N−1
(xω− N + 1 − μω− N + 1) − ΔμTωΔμω +
K
∑i = 1 αigω+ 1(xω+ 1|Ci)
(18)
=
{
n
max 0, ∑ j = 1 β(s)(CK |xj) − K
{
n
V 2
}
∑k = 1 max 0, ∑ j = 1 β(s)(CK |xj) −
V 2
}
(19)
For each Gaussian component, the mean vector μ and the symmetric covariance matrix Σ(s+1) have m and 1/2(m2 + m) scalar parameters, respectively. Therefore, V = 1/2m2 + 3/ 2m denotes the total number of free parameters specifying each Gaussian component. An important feature of eq 20 is that it performs component annihilation. When one of the components becomes “too weak”, meaning that it is not supported by the data, it is simply annihilated. By gradually eliminating those components whose weight falls to zero in each iteration step, the number of the Gaussians is thus adaptively determined. The most significant advantage is that, even in the absence of process knowledge in prior, the algorithm is still able to automatically learn the number of Gaussian components and their parameters from a historical data set. In multimode process monitoring, once the best Gaussians’ number K is initially learned through a training procedure, the number of Gaussian components for each online sample xt could be either K or K+1, since it is possible to add a new operating mode during the practical production. 3.3. The Methodology of Adaptive Monitoring. With the presented moving window GMM algorithms, an adaptive process monitoring scheme can be implemented in practical industry processes. A historical data set of normal states under each operating mode is utilized to train the GMM model. The offline modeling steps are summarized below: (1) Collect a set of historical training data under all possible operating conditions. (2) Use the F-J algorithm to learn the best Gaussian components K, mixture model parameters {μ1,Σ1,···,μK,ΣK} and prior probabilities {α1,···,αK}. (3) Specify a confidence level (1 − α)% and compute the control boundary DL with eq 11. For the online updating and monitoring procedure, an important issue is to distinguish normal process disturbances
where Δμ̃ ω+1 = μω+1 − μ̃. Finally, by substituting eq 13 into eq 15, the recursive calculation of covariance matrix in a moving window scheme can be given as Σω+ 1 = Σω −
αkgω+ 1(xω+ 1|Ck)
(s+1)
(14)
(xω+ 1 − μω+ 1) + Δμω+ 1T Δμω+ 1
(17)
where both gω+1(xω+1|Ck) and βω+1(Ck|xω+1) for each online sample are k × 1 vectors. Through these iterations with realtime samples, the GMM adapts itself to match the latest industrial conditions, and the parameters of GMM carry the latest statistical information of each operating mode. 3.2. Determine the Number of Gaussian Components. The GMM algorithm developed by Figueiredo and Jain adopts the MML (Minimum Message Length) criterion with a variant of EM, where eq 3 in the M-step is modified as follows25
where Δμ̃ ω = μω − μ̃ . With μ̃ and Σ̃, the recursive relations of mean vector and covariance matrix associated with xω+1 are obtained, which is μω+ 1 =
1 (2π)m /2 |Σω+ 1|1/2 ⎡ 1 1 exp⎢ − (xω+ 1 − μω+ 1)T Σ−ω+ 1 ⎣ 2 ⎤ (xω+ 1 − μω+ 1)⎥ ⎦
1 N−1
(xω+ 1 − μω+ 1)T (xω+ 1 − μω+ 1) + ΔμTω+ 1Δμω+ 1 (16) 5499
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
from real faults. To avoid blind updating, some updating rules are demanded to turn on the updating procedure. There are several methods to extract “if-then” rules from chemical processes. The most commonly used method is to establish rule-library with expert knowledge. Besides, Lee et al. adopted a signed graph to generate rules for a refinery process.18 Chen et al. utilized mixture PCA networks to draw expert rules for a simulated continuous stirred tank reactor.26 In short, the general updating rules can be written in the following form: Updating Rule i: IF {change of known disturbance factor Fi is detected} THEN the process is under normal disturbance caused by Fi The updating rule library should be established before the online monitoring is implemented. If the “IF” statement in “Rule i” is not satisfied, then go to “Rule i+1”. If the current situation accords with one of the rules, the updating procedure of Section 3.1 is implemented. To avoid overfiting and turn off the updating timely, the following rule is utilized: Terminating Rule: IF {maxβ(Ck|xt) = 1} THEN the updating for GMM is terminated In addition, for multimode processes, it is possible to add a new operating mode if the current operating conditions are unable to meet the market demands. If a new mode starts with xt, a new Gaussian CK+1 should be added to the model library. The following initiations are utilized αK + 1 = 0;
μK + 1 = xt ;
ΣK + 1 = Σ 0 ;
xk + 1 = Axk + Buk ⎡ 0.67 0.67 0 ⎢ −0.67 − 0.67 0.67 =⎢ ⎢−0.67 0.67 0.67 ⎢ ⎣ 0.67 −0.67 0 ⎡ 0.66 1.55 ⎤ ⎢ ⎥ 1.97 2.38 ⎥ ⎢ + u ⎢ 4.32 − 0.73⎥ k ⎢ ⎥ ⎣−2.64 3.26 ⎦
⎤ 0 ⎥ 0 ⎥x k − 0.67 ⎥ ⎥ −0.67 ⎦
(21)
yk = Cxk + ek + μk ⎡−0.67 ⎢ 2.51 =⎢ ⎢−0.39 ⎢ ⎣−1.33
− 2.66 0.61 1.74 ⎤ ⎥ − 1.28 0.25 − 1.40 ⎥ xk + ek + μk 3.27 − 1.50 0.76 ⎥ ⎥ 1.43 1.26 − 0.93 ⎦
(22)
where xk ∈ R and yk ∈ R are the state variables and process outputs, respectively. The process noise ek ∈ R4 is assumed as Gaussian distribution with ek ∼ (0,0.12). The vector μk ∈ R4 represents the operating center, and it is defined as μk = 1 for initial conditions. The vector uk ∈ R2 denotes process inputs, and they are set to three different modes as follows, to represent three operating conditions (Figure 1) 4
4
Mode 1: uk ∼ (0, 0.052 );
Mode 2: uk ∼ (0.5, 0.032 );
Mode 3: uk ∼ ( −0.5, 0.042)
βK + 1 = 1
(23)
(20)
where Σ0 is the covariance of historical data, and K donates the number of existing Gaussian components. The online modeling and monitoring steps are summarized as follows: (1) For a new online sample xt, calculate its Mahalanobis distance-based integrated monitoring index Dglobal with eq 10. (2) If Dglobal > DL, the process data undergo the predefined updating rules. Each rule is applied to detect the normal changing according to the priority. If the process accords with one of the defined conditional statements, it is identified as normal change. Besides, if xt is the starting point of a new mode, the process is also judged as normal change. Otherwise, xt is classified as an outlier. (3) When detecting normal change, the target Cbest = arg max
Figure 1. Monitored variables under three prespecified modes.
Ck
β(Ck|xt) is updated with the procedure described in Section 3.1. For a new mode event, a new Gaussian CK+1 is initiated with eq 20. (4) If the terminating rule is satisfied, the updating is turned off. (5) If several consecutive samples are judged as outliers, the faulty alarms are triggered.
Four outputs and two inputs, altogether six process variables are measured for monitoring. Under each mode, 600 samples are generated, and altogether 1,800 samples are collected and used as historical data. The offline training results with the GMM algorithm is displayed in Figure 2. As observed, three ellipses with different mean and covariance represent the 95% confidence boundary of three operating modes, respectively. Case 1: Process Slow Drift. To simulate the dynamic feature of an industrial process, a slow drift is introduced to the input variable uk under Mode 1. The inputs uk gradually drift away from its normal operating condition (donated as u*k) from the number serial i = 301 to i = 600 by the relation of uk =
4. A NUMERICAL EXAMPLE The proposed adaptive monitoring approach is applied to a numerical example to validate its effectiveness. Consider the following 4 × 2 linear dynamic process24 5500
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
statistics of adaptive GMM stay steadily under the control boundary until the true fault occurs. As expected, the adaptive procedure keeps the monitoring model up-to-date and lower down type II errors significantly. Figure 4 shows the posterior probability plot of test samples under each operating mode. Figure 5 displays the scatter plot of
Figure 2. Training results of GMM algorithms with historical data.
u*k + 0.003(i − 300). The drifting behaviors are reflective of normal process changes such as catalyst deactivation or equipment fouling in real chemical industrial processes. Then, a step process fault is introduced such that the first system output y1 gets an increase of “3” from time t = 601 to t = 900. With the 900 online samples, the monitoring result of the proposed adaptive GMM monitoring scheme is compared with that of a static GMM method in Figure 3. Figure 3(a) shows an Figure 4. Posterior probabilities of test data with static GMM (a) under Mode 1; (b) under Mode 2; and (c) under Mode 3.
Figure 5. Some test samples and their locations to static GMMs under Case 1. Figure 3. Monitoring results of (a) the static GMM method and (b) the adaptive GMM method under Case 1.
some typical test samples, where black stars “*” represent test samples and their serial numbers. These two figures can make a more intuitive explanation of the monitoring result of static GMM. Specifically, the first test sample initially locates within the normal area of Mode 1, corresponding to posterior probability “1” in for Mode 1 and “0” for Modes 2 and 3 in Figure 4. As the process drifts, test samples gradually move closer to the 95% control boundary, as shown by Sample 384− 410, whose posterior probabilities oscillate up and down during these periods. When the drift goes further, test samples move
increasing number of violations of the 95% confidence boundary (red line) arise after 350th samples, and then Dglobal statistics fall down to the normal region for the 430−480th sample . After that, the statistics rise even higher, and the fault alarms are continuously triggered. The step fault since the 600th sample has no reflection in the monitoring chart because the real fault symptoms are buried in variations caused by process drifting. In contrast, Figure 3(b) illustrates that Dglobal 5501
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
can provide us more accurate information about process operating conditions and therefore make it easier for subsequent fault isolation and diagnosis procedures. Case 2: New Mode Occurs. Another experiment is set to validate the new-event-tolerance ability of the proposed method. From time t = 1 to t = 300, the process runs steadily under Mode 2; while from t = 301 to t = 900, a new operating mode is introduced such that the system matrix C changes into 2C. In real application this situation happens due to the various changes of product specification, feed flow rate, compositions, and set points, etc. After the process running under the new mode for 300 normal samples, a drift fault is introduced into the system from t = 601 to t = 900. The operating center drifts in the form of μk+1 = 1.003 μk + ek, where ek ∼ (0,0.12). The monitoring results of the static GMM and adaptive GMM methods are shown in Figure 8(a) and (b), respectively.
into the normal areas of Mode 2. Samples 427−521 clearly show that the posterior probabilities turn to “1” for Mode 2 and “0” for Mode 1 and 3. This explains why there is a significant drop at around 450th samples in Figure 3(a). After the fault occurs, test samples go farther from their normal areas and belong to none of the existing three modes, Samples 700 and 800 as typical. For the adaptive GMM, with the guide of maximum posterior probabilities, the corresponding model updates itself with process changes. Figure 6 illustrates that the 95%
Figure 6. Test data and adaptive GMMs---from the initial model evolve to the medium and then to the final one.
confidence ellipse changed its envelope with test samples until the model parameters become stable. In this procedure, as shown in Figure 7, besides a few confusions with Mode 2 around the 450th sample, for most of the test data, their posterior probabilities stay around “1” for Mode 1 and around “0” for Modes 2 and 3. This clear specification for online data
Figure 8. Monitoring results of (a) the static GMM method and (b) the adaptive GMM method under Case 2.
For the former approach, the monitoring statistics behave normally for the first 300 samples, and after that, they rise up immediately and stay high above the control boundary until the end of the experiment. This monitoring result is usually considered as a step fault at Sample 300, and alarms will consequently be triggered. The “false alarms” conceal symptoms of the real fault. On the contrary, with the proposed adaptive approach, the monitoring chart is displayed in Figure 10. An obvious pulse at the 301st sample is a sign of process changes, and monitoring statistics perform normally afterward until Sample 700. The 100-sample delay is attributing to the characteristic of drifting error, as the center-drift fault drives process leaving its normal conditions gradually and continuously. The two cases of this numerical experiment verify that the adaptive GMM outperforms static GMM in building monitoring models for multimode processes with dynamic behaviors.
Figure 7. Posterior probabilities of test data with adaptive GMM (a) under Mode 1; (b) under Mode 2; and (c) under Mode 3. 5502
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
5. TE SIMULATION The TE simulation program is widely applied to evaluate the efficiency of process monitoring techniques. The original version of the simulation program is provided in FORTRAN scripts, where the plant worked under open-loop and would shut down after running two hours due to the high pressure in reactor.27 In our experiment, the decentralized control strategy proposed by Ricker is implemented to generate simulating process data.28 The Matlab simulation codes can be downloaded from Ricker’s Web site: http://depts.washington.edu/ control/LARRY/TE/download.html. There are six prespecified operating modes in the platform to meet the demands of different production grades, among which three are chosen in this experiment, as shown in Table 1, where
operating mode. The impact of this disturbance is widespread within the process, and most of the monitored variables are affected. Random variation is commonplace in real industrial environment, and the control mechanism is designed to be tolerant to the impact of random variations. Figure 9(a) gives
Table 1. Three Operating Modes Used in Experiment mode
G/H mass ratio
production rate
1 2 3
50/50 10/90 90/10
7038 kg/h G and 7038 kg/h H 1111 kg/h G and 10,000 kg/h H 10,000 kg/h G and 1111 kg/h H
G and H are the final products of the TE plant. The process is set to running under each mode for 25 h. The sampling time used is 3 min; hence 500 samples are generated under each operating mode. Altogether, the whole 1500 samples from all the three modes constitute the historical data, covering a time range of 75 h. Moreover, there are 51 process variables made up of three parts, i.e. 22 continuous measurements, 17 component measurements, and 12 manipulating variables, among which the 22 continuous measurements are chosen as monitored variables in our simulation.27 The TE platform also provides 20 process disturbances under each operating mode, and detailed descriptions can be found in ref 27. Our experiments also verified that 19 of the disturbances could be suppressed successfully by Ricker’s control mechanism, and the plant is able to produce G and H normally, with the exception of IDV(6). In this simulation, two scenarios are set to verify the effectiveness of the proposed monitoring method, as displayed in Table 2.
Figure 9. Monitoring results of (a) the static GMM method and (b) the adaptive GMM method under Scenario 1.
the monitoring result of the static GMM method under the 95% confidence level. The monitoring statistics jump through the control limit and then fluctuate with the random variation since the 15th hour, which trigger continuously false alarms. The real fault is also reflected in this monitoring result by a sudden increase at the 30th hour. However, we expect to eliminate the unnecessary false alarms and reveal the real fault. Figure 9(b) shows the monitoring result of adaptive GMM, which illustrates clear and significant fault symptoms since the 30th hour. Before the valve of steam 4 get sticking, the monitoring model is robust to normal random variations, and as we expected, the Dglobal statistics stay steadily under their control limit with few false alarms. This scenario shows that the adaptive GMM method can model and monitor processes with normal random variations effectively under control mechanism. The following simulation further approves that the adaptive scheme is also efficient to processes with drifting disturbances. The second case involves a slow drift of D feed temperature in stream 2. This disturbance gradually spreads to the entire system after the 15th hour and then impacts all the process variables. Although it is regulated by controllers, the variations of monitored variables are still reflected on the monitoring chart, as shown in Figure 10(a), which is the monitoring result of the static GMM approach. Due to the nature properties of slow drift, Dglobal statistics go beyond their control limits about 1 h later after the disturbance occurred, after that they gradually go higher beyond the control boundary and give us a typical sign of fault. In fact the control system is designed to be tolerant of this slow drift and guarantee steady production already, so this kind of alarms do not make sense. Since the
Table 2. Two Scenarios Set in This Simulation mode Mode 1
Mode 3
state description steady state random variation of condenser cooling water inlet temperature the valve of stream 4 fixed in a stable position steady state slow drift of D feed temperature in stream 2 the valve of A feed get sticking
starting time initial 15th hour 30th hour initial 15th hour 30th hour
The experimental time for each scenario is 45 h with 900 test samples generated, where the first 300 samples are collected under steady state of the specified operating mode. The process disturbances begin from the 15th hour; the real faults are introduced at the 30th hour and last until the end of test time. Both of the process disturbances are adjusted successfully by the control mechanism, and the monitoring model should be robust to them while detecting faults. In the first scenario, we consider a random variation of condenser cooling water inlet temperature under the first 5503
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This research is supported by the National Natural Science Foundation of China (No. 61074079) and Shanghai Leading Academic Discipline Project (No. B504).
■
(1) Qin, S. J. Statistical process monitoring: basics and beyond. J. Chemom. 2003, 17, 480−502. (2) Chen, J. H.; Liu, J. L. Mixture principal component analysis models for process monitoring. Ind. Eng. Chem. Res. 1999, 38, 1478− 1488. (3) Zhao, S. J.; Zhang, J. Monitoring of processes with multiple operation modes through multiple principle component analysis models. Ind. Eng. Chem. Res. 2004, 43, 7025−7035. (4) Zhao, S. J.; Zhang, J.; Xu, Y. M. Performance monitoring of processes with multiple operating modes through multiple PLS models. J. Process Control. 2006, 16, 763−772. (5) Zhao, C. H.; Wang, F. L.; Lu, N. Y.; Jia, M. X. Stage-based softtransition multiple PCA modeling and on-line monitoring strategy for batch processes. J. Process Control. 2007, 17, 728−741. (6) Liu, J. L. Fault Detection and Classification for a Process with Multiple Production Grades. Ind. Eng. Chem. Res. 2008, 47, 8250− 8262. (7) Ng, Y. S.; Srinivasan, R. An adjoined multi-model approach for monitoring batch and transient operations. Comput. Chem. Eng. 2009, 33, 887−902. (8) Liu, J.; Chen, D. Nonstationary fault detection and diagnosis for multimode processes. AIChE J. 2010, 56, 207−219. (9) Choi, S. W.; Park, J. H.; Lee, I. B. Process monitoring using a Gaussian mixture model via principal component analysis and discriminant analysis. Comput. Chem. Eng. 2004, 28, 1377−1387. (10) Yoo, C. K.; Villez, K.; Lee, I. B.; Rosén, C.; Vanrolleghem, P. A. Multi-model statistical process monitoring and diagnosis of a sequencing batch reactor. Biotechnol. Bioeng. 2007, 96, 687−701. (11) Yu, J.; Qin, S. J. Multimode process monitoring with Bayesian inference-based finite Gaussian mixture models. AIChE J. 2008, 54, 1811−1829. (12) Yu, J.; Qin, S. J. Multiway Gaussian Mixture Model Based Multiphase Batch Process Monitoring. Ind. Eng. Chem. Res. 2009, 48, 8585−8594. (13) Chen, T.; Zhang, J. On-line multivariate statistical monitoring of batch processes using Gaussian mixture model. Comput. Chem. Eng. 2010, 34, 500−507. (14) Jie, Y. A nonlinear kernel Gaussian mixture model based inferential monitoring approach for fault detection and diagnosis of chemical processes. Chem. Eng. Sci. 2012, 68, 506−519. (15) Wold, S. Exponentially weighted moving principal components analysis and projections to latent structures. Chemom. Intell. Lab. Syst. 1994, 23, 149−161. (16) Wise, B. M.; Gallagher, N. B. The process chemometrics approach to process monitoring and fault detection. J. Process Control. 1996, 6, 329−348. (17) Li, W. H.; Yue, H. H.; Valle-Cervantes, S.; Qin, S. J. Recursive PCA for adaptive process monitoring. J. Process Control. 2000, 10, 471−486. (18) Lee, Y. H.; Jin, H. D.; Han, C. H. On-line process state classification for adaptive monitoring. Ind. Eng. Chem. Res. 2006, 45, 3095−3107. (19) Jin, H. D.; Lee, Y. H.; Lee, G.; Han, C. H. Robust recursive principal component analysis modeling for adaptive monitoring. Ind. Eng. Chem. Res. 2006, 45, 696−703. (20) Choi, S. W.; Martin, E. B.; Morris, A. J.; Lee, I. B. Adaptive multivariate statistical process control for monitoring time-varying processes. Ind. Eng. Chem. Res. 2006, 45, 3108−3118.
Figure 10. Monitoring results of (a) the static GMM method and (b) the adaptive GMM method under Scenario 2.
process drifts from its steady state, the monitoring model should update itself and capture the changes timely. The monitoring result of adaptive GMM method is displayed in Figure 10(b). As observed, all the monitoring statistics stay below their control boundaries until the valve of A feed get sticking. The monitoring model adapts with this slow drift and functions well from the beginning to the end of the simulation. The different test scenarios in this TE experiments demonstrate the outstanding potential of the proposed approach in building adaptive models for dynamic multimode industrial processes.
6. CONCLUSIONS A novel adaptive multimode process monitoring method based on moving window GMM has been proposed in this article. The original GMM algorithm was used to train initial models with historical data. Through the indication of maximum posterior probabilities of online samples, the best-suit Gaussian model updates itself with latest process changes, while other Gaussians stay unchanged unless the operating mode shifts to another one. With established “updating rules”, our adaptive monitoring approach minimizes the human intervening. A numerical example and the TE simulation platform were used to validate the effeciency of the proposed method. The comparison of the monitoring results demonstrated that the adaptive GMM is superior to static GMM in the circumstances of process drifing, new mode appearance, and random variation of variables. The new approach provides a promising tool to modeling and monitoring dynamic industrial processes under multiple operating conditions.
■
REFERENCES
AUTHOR INFORMATION
Corresponding Author
*Phone: +86-021-64252189. Fax: +86-021-65252189. E-mail:
[email protected]. 5504
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505
Industrial & Engineering Chemistry Research
Article
(21) Liu, H. C.; Jiang, W.; Tangirala, A.; Shah, S. An adaptive regression adjusted monitoring and fault isolation scheme. J. Chemom. 2006, 20, 280−293. (22) Wang, X.; Kruger, U.; Irwin, G. W. Process monitoring approach using fast moving window PCA. Ind. Eng. Chem. Res. 2005, 44, 5691−5702. (23) Liu, X. Q.; Kruger, U.; Littler, T.; Xie, L.; Wang, S. Q. Moving window kernel PCA for adaptive monitoring of nonlinear processes. Chemom. Intell. Lab. Syst. 2009, 96, 132−143. (24) Jeng, J. Adaptive process monitoring using efficient recursive PCA and moving window PCA algorithms. J. Taiwan Inst. Chem. Eng. 2010, 41, 475−481. (25) Figueiredo, M.; Jain, A. K. Unsupervised learning of finite mixture models. IEEE Trans. Pattern Anal. 2002, 24, 381−396. (26) Chen, J. H.; Liu, J. L. Using mixture principal component analysis networks to extract fuzzy rules from data. Ind. Eng. Chem. Res. 2000, 39, 2355−2367. (27) Downs, J. J.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245−255. (28) Ricker, N. L. Decentralized control of the Tennessee Eastman challenge process. J. Process Control 1996, 6, 205−222.
5505
dx.doi.org/10.1021/ie202720y | Ind. Eng. Chem. Res. 2012, 51, 5497−5505