Multimode Dynamic Process Monitoring Based on Mixture Canonical

Jan 8, 2015 - In this article, a novel mixture canonical variate analysis (MCVA) model is proposed to model and monitor multimode dynamic processes. F...
1 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/IECR

Multimode Dynamic Process Monitoring Based on Mixture Canonical Variate Analysis Model Qiaojun Wen, Zhiqiang Ge,* and Zhihuan Song State Key Laboratory of Industrial Control Technology, Institute of Industrial Process Control, Department of Control Science and Engineering, Zhejiang University, Hangzhou 310027, Zhejiang, People’s Republic of China ABSTRACT: For complex industrial processes with multiple operating conditions and dynamic characteristics, the traditional dynamic process monitoring techniques (e.g., canonical variate analysis, CVA) are not well-suited, because the fundamental assumption that the operating data follow a unimodal Gaussian distribution usually becomes invalid. In this article, a novel mixture canonical variate analysis (MCVA) model is proposed to model and monitor multimode dynamic processes. First, the augmented process data are assumed to be many different clusters, each of which corresponds to an operating mode and can be characterized by a Gaussian component. Then, singular value decomposition of the covariance matrices is implemented in each Gaussian cluster and the corresponding canonical variates are obtained. For process monitoring purposes, the local statistics in each cluster are calculated and the integrated monitoring indices are obtained in a probabilistic manner. The validity and effectiveness of the proposed monitoring approach are illustrated through two examples: a numerical process and the Tennessee Eastman challenge problem. The comparison of monitoring results demonstrates that the proposed approach is superior to conventional CVA and Gaussian mixture models (GMM).

1. INTRODUCTION In modern chemical engineering applications, process safety and high-quality products are highly pursued, which leads to an increasing requirement for fault detection and diagnosis techniques. Since large amount of process data are available, multivariate statistical process monitoring (MSPM) methods are widely used in both academic research and industrial applications.1−5 Process dynamics and nonlinearity are common in industrial processes. Because of the intrinsic characteristics of the plants and the closed-loop control systems, process measurements are usually temporally dependent. To monitor dynamic processes, different MSPM methods have been developed.6 Among these methods, dynamic principal component analysis (DPCA) is a popular approach.7 DPCA is a dynamic generalization of the conventional principal component analysis (PCA) method. In DPCA, first, the augmented vectors are defined, which contain time-lagged process measurements. Therefore, dynamic information is captured in the subsequent PCA model. To reduce the increase in false alarms in DPCA, univariate ARMA filters were proposed, to remove autocorrelation from the DPCA score variables.8 Besides DPCA, the state space representations are also widely used to model dynamic processes.9,10 Most of the research on the state space representations used subspace identification methods (e.g., numerical subspace state space system identification (N4SID), error in variable (EIV) identification, canonical variate analysis (CVA)).11−15 In the CVA model, the past information vectors (pt), as well as the future information vectors (ft), are defined; in addition, the dynamic processes are modeled by capturing the relationship between the two vectors. The calculation of canonical variates is based on singular value decomposition of the covariance matrices of the past information vectors and the future information vectors. © 2015 American Chemical Society

Many industrial processes are nonlinear and operate in different operation regimes, because of different product specifications, working environments, and economic considerations.16 To monitor general nonlinear processes, kernel methods have been proposed, such as kernel PCA (KPCA),17 kernel independent component analysis (KICA),18 and kernel partial least-squares (KPLS).19 For multimode processes, specifically, mixture methods, e.g., Gaussian mixture models (GMM), are powerful methods.20,21 Besides GMM, similar mixture models have been developed in process monitoring, such as mixture probabilistic principal component analysis (MPPCA)22,23 and mixture factor analysis (MFA).24 However, most of the previous research in multimode process monitoring are based on the assumption that the process measurements are independent and identically distributed (i.i.d.). In other words, process dynamics are not incorporated in the statistical models. In this paper, a new approach based on mixture canonical variate analysis (MCVA) model is proposed to monitor processes with both dynamics and multimode characteristics. The multimode process data are clustered into different operation modes via GMM modeling at first. The covariance matrices of different modes are estimated via GMM as well. Then, based on the covariance matrices estimated by GMM, a canonical variate analysis model is constructed in each operation mode and the corresponding monitoring statistics in the operation mode are computed separately. Lastly, a Bayesian inference procedure is utilized to compute the posterior probabilities of each sample and integrate the statistics in different modes to obtain the integrated monitoring Received: Revised: Accepted: Published: 1605

August 21, 2014 November 19, 2014 January 8, 2015 January 8, 2015 DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research

Ts2 = ptTJkT Jk pt

indices. In MCVA, it is unnecessary to assign a cluster label to each sample in both the calculation of canonical variates and the calculation of the global monitoring indices. Instead, the posterior probabilities of each sample are calculated. Therefore, the proposed MCVA approach avoids the potential risk of false identification induced by misclassification of the samples. The rest of this article is organized as follows. In the next section (section 2), the preliminary knowledge of CVA and GMM are briefly reviewed, followed by the proposed mixture canonical variate analysis model, as well as the corresponding process monitoring scheme (section 3). In section 4, two case studies are provided to demonstrate the efficiency of the proposed method. Finally, conclusions are made in section 5.

If the process variables are normally distributed, the control limit of T2s can be calculated via the following F-distribution:25 Ts2 ≈

Tn2 ≈

...,

(10)

G

p(x) =

∑ πg 5(x|μg , Σg ) (11)

g=1

pt = [ z̅ tT, z̅ tT− 1, ..., z̅ tT− l , u̅ tT, u̅ tT− 1 , ..., u̅ tT− l]T ft =

q(n2 − 1) F(q , n − q) n(n − q)

Both eqs 8 and 10 are based on the assumption that the process variables are Gaussian-distributed. Therefore, in the case that the process is nonlinear or multimode, conventional CVA is not suitable for process monitoring. 2.2. Gaussian Mixture Models. The distribution of random variables x ∈ ℜD from Gaussian mixture models can be written as a linear superposition of Gaussian components in the form26

(1)

z̅ tT+ h]T

(9)

Similarly, the control limit of the noise space statistics, T2n, can be calculated via the following F-distribution:1

where ηt, ut, and zt denote the process state, input, and output vectors, respectively; the term wt denotes the process noise (wt ≈ 5(0, Q)) and the term vt denotes the observation noise (vt ≈ 5(0, R )); each is independent and identically distributed. The centered stacked past and future vectors at time instance t, pt and ft, are represented as follows:

z̅ tT+ 2,

(8)

Tn2 = ptTJqT Jq pt

ηt + 1 = Aηt + Bu t + wt

[ z̅ tT+ 1,

k(n2 − 1) F (k , n − k ) n(n − k)

To monitor the process variation out of the state space (e.g., the process noise), T2n are calculated as

2. PRELIMINARY KNOWLEDGE 2.1. Canonical Variate Analysis. The state space model of a dynamic process can be represented as zt = Cηt + Du t + vt

(7)

where G is the number of Gaussian components, each of which has mean μg and covariance matrix Σg, g = 1, 2, ..., G; and the non-negative weights πg, which are called mixing probabilities, satisfy the constraint ∑Gg=1πg = 1. Let θg = {μg,Σg}, in the gth component Cg, the Gaussian density function can be expressed as

(2)

where l and h are the number of the lag and the lead, respectively, and zk̅ and uk̅ are the centered process input and output vectors. Canonical variate analysis (CVA) models the relationship of the past vectors (p) and the future vectors (f). In CVA, canonical variates s and r are calculated as the linear combinations of p and f, respectively.

p(x|x ∈ Cg ) = N (x|θg ) =

1 1 exp − (x − μg )T Σ−g 1(x − μg ) 2 (2π )D /2 |Σg |1/2

{

}

s = Jp

(12)

r = Lf

The model parameters can be written as Θ = {θ1, ..., θG, π1, ..., πG}. Given a set of N independent and identical distributed samples X = {x1, x2, ..., xN}, the log-likelihood of the model is

(3)

where J and L are the transformation matrices that maximize the correlation between s and r, respectively. According to CVA, the loading matrix J is obtained by solving the following singular value decomposition problem:1

i=1

Σ−pp1/2Σpf Σ−ff1/2 = UDVT

st = Jk pt

i=1

g=1

(13)

(5)

To estimate the model parameters, the maximum likelihood estimate (MLE) or maximum a posteriori (MAP) can be used. The MLE is achieved by

In eqs 4 and 5, Σpp = Σpf = and Σff = E(ftfTt ) are the covariance matrices of pt and ft. In eq 4, D = diag(ρ1, ..., ρd, 0, ..., 0) where ρ1, ..., ρd are the correlation of different canonical variates, e.g., ρ1 = corr(s1, r1), and ρ1 ≥ ··· ≥ ρd.1 Let k be the order of states, Jk be the first k rows of J, and Jq be the other q rows of J. In CVA, the states at time instance t, st ∈ ℜk , are calculated as E(ptpTt ),

G

∑ log ∑ πgp(x i|θg )

(4)

and J = UTΣ−pp1/2

N

N

log p(X|Θ) = log ∏ p(x i|Θ) =

Θ̂ = arg max log p(X|Θ)

E(ptfTt ),

(14)

Θ 26

The solution of eq 14 cannot be found analytically. Instead of ordinary MLE methods, the expectation−maximization (EM) algorithm can be used to solve the problem. The EM algorithm is an iterative algorithm that finds the local maxima of loglikelihood function. However, the objective function of EM algorithm is different from eq 13. According to the EM algorithm for GMM,26 X is treated as an incomplete dataset in which the missing part of GMM can be interpreted as N tags Z

(6)

and the monitoring statistics for the states, T2s , are calculated as 1606

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research

⎡ μg , p ⎤ μg = ⎢ ⎥ ⎢⎣ μg , f ⎥⎦

= {z1, z2, ..., zN}. Each missing part zi is a G-dimensional binary random variable, representing the Gaussian component in which the sample xi is generated. In cooperation with the missing part, the objective function is the complete loglikelihood Θ̂ = arg max log p(X, Z|Θ)

and ⎡ Σg , pp Σg , pf ⎤ ⎥ Σg = ⎢ ⎢⎣ Σg , fp Σg , ff ⎥⎦

(15)

Θ

With a previously known G, the parameters Θ can be learned from eq 15 via the EM algorithm, as follows: E-Step: At the lth iteration, evaluate the posterior probabilities using the current parameters Θ(l): p(l) (x i ∈ Cg ) = p(l) (Cg |x i) =

G

∑ j = 1 πj(l)p(x i|θj(l)) (16)

p(yt ∈ Cg |pt , f t) = p(Cg |pt , f t) =

M-Step: Re-estimate the parameters using the current posterior probabilities: ∑i = 1 p(l) (Cg |x i)x i N ∑i = 1 p(l) (Cg |x i)

(17)

N

N

∑i = 1 p(l) (Cg |x i)

(18)

N

∑i = 1 p(l) (Cg |x i)

πg(l + 1) =

G

pp (p|θg ) =

N

∑ j = 1 ∑i = 1 p(l) (Cj|x i)

(24)

where p(yt|θg) is calculated via the Gaussian distribution yt ≈ 5(μg , Σg ). In online monitoring, since the future vector f is unknown, the calculation of the corresponding posterior probabilities to each Gaussian component is based on the marginal probability of p,

∑i = 1 p(l) (Cg |x i)(x i − μg(l + 1) )(x i − μg(l + 1) )T

Σ(gl + 1) =

πgp(yt |θg ) G ∑ j = 1 πjp(yt |θj)

t = 1 , ..., N ; g = 1 , ..., G

N

μg(l + 1) =

(23)

where μg,p and μg,f denote the mean vectors of p and f in the gth model, Σg,pp and Σg,ff denote the covariance matrices, Σg,pf and Σg,fp denote the cross-covariance matrices. The corresponding posterior probabilities to each Gaussian component can be calculated using the Bayes formula. In offline modeling, specifically, since both p and f are known,

πg(l)p(x i|θg(l))

i = 1 , ..., N ; g = 1 , ..., G

(22)

∫ p(p, f|θg) df

(25)

(19)

Since the joint distribution is p(p, f|θg ) ≈ 5(μg , Σg ), pp(p|θg)

The above iterations will stop when the parameters converge. For some multimode processes, the number of Gaussian components can be determined with process knowledge. For general nonlinear processes, however, it is not easy to set an appropriate G for GMM. With an unknown G, the process parameters {G, Θ} can be estimated using several methods, such as variational Bayesian GMM26 and the F−J algorithm.27 Recently, Feital et al. proposed an efficient algorithm to identify the parameters of the GMM model, as well as to select the number of Gaussian components.28 Moreover, Feital’s approach is computationally simple and suitable for GMM clustering in large-scale processes.

can be calculated from the distribution p|θg ≈ 5(μg , p , Σg , pp). Therefore, the corresponding posterior probabilities to each Gaussian component are calculated as p(yt ∈ Cg |pt ) = p(Cg |pt ) = t = 1 , ..., N ; g = 1 , ..., G

T −1/2 Σ−g ,1/2 pp Σg , pf Σ g , ff = Ug Dg Vg

sg = Jg , k (p − μg , p )

∑ πg 5(y|μg , Σg ) (20)

(28)

where Jg,k denotes the first k columns of the loading matrix Jg and

where

⎡ p⎤ y=⎢ ⎥ ⎣f ⎦

(27)

The corresponding CVA states are calculated as follows:

G g=1

(26)

Since the calculation of loading matrices in CVA is only dependent on the covariance matrices of the cluster, therefore, the subsequent identification procedure based on CVA is made in the GMM probabilistic framework and it is unnecessary to assign a cluster label to each sample. Similar to eq 4, CVA model can be used to model the relationship of p and f in the gth Gaussian component (g = 1, ..., G),

3. MIXTURE CANONICAL VARIATE ANALYSIS 3.1. Mixture Canonical Variate Analysis Formulation. To model multimode dynamic processes, the CVA model can be generalized to a mixture form. Consider the joint probability distribution of past and future vectors, p and f, in a mixture of Gaussian distributions, p(y) =

πgpp (pt |θg ) G ∑ j = 1 πjpp (pt |θj)

Jg = UTg Σ−g ,1/2 pp (21)

(29)

In CVA, the order of the states k can be determined via Akaike information criteria (AIC).1 In this article, the correlation explanation ratio is proposed to determine the order of states sg in the gth Gaussian component. In eq 27, Dg = diag(ρ1,g, ..., ρdg,g, 0, ..., 0), where ρ1,g, ..., ρdg,g are the correlation

The GMM model then is trained with the collected multimode past and future matrices. Note that the mean vector μg and covariance matrix Σg in the gth Gaussian model can be partitioned as 1607

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research of different canonical variates in the gth Gaussian component and ρ1,g ≥···≥ ρdg,g. Since CVA maximizes the correlation between the canonical variates s and r, the cumulative correlation explanation (CCE) of the first k canonical variates in the gth Gaussian component can be defined as

contribution of each local Gaussian components to the overall probabilistic index. Without classifying the monitored sample into a single cluster deterministically, such a probabilistic strategy can avoid the potential risk of false detection induced G by misclassification. As 0 ≤ plocal g,s (Pt) ≤ 1 and ∑g=1p(yt∈Cg|pt) = 1, we have

k

CCE(k) =

∑i = 1 ρi , g

0 ≤ BIPs ≤ 1

d

∑ j =g 1 ρj , g

(30)

Under a prespecified confidence level 1 − α (typically, 1 − α = 99%), the process is determined within normal operation if

Then, a proper value of k can be determined via the cumulative correlation explanation. For example, k is determined to make CCE(k) >90%. 3.2. MCVA-Based Monitoring Index. With the constructed MCVA model, it is necessary to further derive the confidence boundary around the normal operating regions for process monitoring and fault detection. A two-step method is applied in the process monitoring scheme. The first step is to calculate the local monitoring statistics in the different Gaussian components separately, while the second step is to integrate the local statistics and obtain the integrated indices. In the gth model, the statistics of the state space is calculated as Ts2, g

= (pt −

μg , p )T JkT, g Jk , g (pt

− μg , p )

BIPs ≤ 1 − α

k(ng2 − 1) ng (ng − k)

(36)

Similarly, the statistic of the noise space is calculated as Tn2, g = (pt − μg , p )T JqT, g Jq , g (pt − μg , p )

(37)

where the corresponding distribution can be calculated as Tn2

≈

q(ng2 − 1) ng (ng − q)

F(q , ng − q) (38)

For the monitored sample pt, the local probability index in the noise space, relative to the gth Gaussian component, is defined as

(31)

If the distribution of Pt is non-Gaussian, the distribution of T2s,g can be estimated via Monte Carlo simulation. In this article, it is assumed that observations in each cluster follow a Gaussian distribution. In the gth Gaussian component, the estimation of the distribution of T2s,g is simpler. Similar to eq 8, the distribution of T2s,g can be approximated as Ts2, g ≈

(35)

pglocal (pt ) = p{Tg2, n(p) ≤ Tg2, n(pt )} ,n

∀p

(39)

in normal operation of the gth component, which can be calculated by integrating the probability density function for the F-distribution in eq 38. In addition, the integrated index is given by G

F(k , ng − k)

BIPn =

(32)

g=1

where ng = Nπg denotes the effective number of samples assigned to the gth Gaussian component. To calculate the integrated indices, Bayesian inference is implemented, similar to the Bayesian inference probability (BIP) in a GMM-based monitoring scheme.21 For the monitored sample pt, the local probability index in the state space relative to the gth Gaussian component is defined as pglocal (pt ) = p{Tg2, s(p) ≤ Tg2, s(pt )} ,s

∀p

(pt ) ∑ p(yt ∈ Cg |pt )pglocal ,n

(40)

The online monitoring scheme of the MCVA-based monitoring approach is illustrated in Figure 1.

4. CASE STUDIES In this section, two case studies are demonstrated to show the effectiveness of the proposed approach, one of which is a

(33)

in normal operation of the gth component, which can be calculated by integrating the probability density function of the F-distribution in eq 32. Under a given confidence level, this index serves as an indication of whether the monitored sample is normal or faulty, provided that it belongs to the corresponding Gaussian component. Considering the random characteristic that each monitored sample may come from multiple Gaussian components with the corresponding posterior probabilities, an integrated index is further defined to combine the local probability metrics across all the Gaussian clusters. The formulation of the integrated BIP index for the monitored sample pt is given by G

BIPs =

(pt ) ∑ p(yt ∈ Cg |pt )pglocal ,s g=1

(34)

where the posterior probability p(yt ∈ Cg|pt) calculated via Bayesian inference in eq 26 is used to incorporate the

Figure 1. Online monitoring scheme of MCVA-based process monitoring approach. 1608

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research

that the dimensionalities of both p and f are 8. Two Gaussian components were chosen in the mixture model. For performance comparison, two process monitoring models, which were based on conventional CVA and static GMM, respectively, were also constructed using the same training data samples, respectively. Similar to the training data, a total of 2000 samples in the two normal operation modes were used to test the false alarm rates of the models. The false alarm rates of the three models are listed in Table 1. It is showed that the false alarm rates in

numerical case while the other is the well-known benchmark simulation of Tennessee Eastman (TE) process. 4.1. Numerical Example. The numerical example is constructed with two operation modes. Specifically, the model is constructed using the following equations: x(t ) = Bx(t − 1) + ξ(t ) y(t ) = Cix(t ) + ζ (t )

(41)

where i represents the two different operation modes (i = 1, 2), t is the sample number of process data (t = 1, 2, ..., n), x(t) represents the process states (x(t ) ∈ ℜ2 ), and y(t) represents the measured process variables (y(t ) ∈ ℜ4 ). ξ and ζ are Gaussian residuals with a zero mean and a variance of 0.01. (ξ ∈ ℜ2 and ζ ∈ ℜ4 .) The state transfer matrix is given as

Table 1. False Alarm Rates in the Numerical Process MCVA

⎡ 0.8 −0.1⎤ B=⎢ ⎥ ⎣ 0.4 0.6 ⎦

CVA

BIPs

BIPn

T2s

0.010

0.010

0.023

GMM T2n

BIP

0.014

0.008

both the MCVA model and GMM are no larger than 1%. However, in the conventional CVA model, the false alarm rates are slightly higher than expected. To test the fault detection performance of the proposed monitoring approach, three different faulty datasets from the first operation mode (each contains 2000 data samples) were also generated, which are listed as follows: (1) A step change of the second process variable by 0.5 was introduced starting from sample 1001; (2) A ramp change of the third process variable from sample 1001 to 2000 was introduced by adding 0.002(i − 1000) to each data sample, where i was the sample number; (3) Process structure changed from sample 1001 to 2000 by changing matrix B to

The observation matrix in the first operation mode is ⎡ −0.731 −0.019 −0.153 1.028 ⎤T C1 = ⎢ ⎥ ⎣−2.417 0.902 0.764 0.486 ⎦

while the observation matrix in the second operation mode is ⎡ 0.731 − 0.019 0.153 1.028 ⎤T C2 = ⎢ ⎥ ⎣ 2.417 0.902 − 0.764 0.486 ⎦

For model construction, 1000 successive data samples were generated in the first operation mode and another 1000 samples were generated in the second operation mode. All of the data were collected under stationary operation conditions. Therefore, a total of 2000 samples were used to construct the model. Figure 2 illustrates the distribution of the first two

⎡ 0.8 −0.1⎤ B=⎢ ⎣−0.6 0.6 ⎥⎦

The first two testing datasets (faults 1 and 2) represent some abnormal events with the process sensors: fault 1 represents sensor bias and fault 2 represents sensor drift. The last testing dataset (fault 3) represents some process fault, which is a change in process structure. The detailed fault detection results with fault 1 are shown in Figure 3. Figure 3b shows the performance of conventional CVA. One can see that the fault detection rate of both the T2s and T2n statistics in CVA are not large. Specifically, the fault detection rate of the T2n statistics is only 18.0%. As showed in Figure 3c, the fault detection rate of GMM-based BIP index is 76.8%, which means that several faulty samples are still missed (not detected). In contrast, the proposed MCVA performs much better than these two methods. In MCVA, the fault detection rate of the BIPs and BIPn indices are 99.3% and 86.3% respectively, which means almost all of the faulty samples are identified in fault 1. The fault detection rates of the three models are listed in Table 2. In each fault scenario, the one with the best detection performance is shown in boldface font. It is shown that the fault detection rates of both state space monitoring indices BIPs and noise space monitoring indices BIPn in MCVA model are larger than that in the GMM-based model and the conventional CVA model. 4.2. Tennessee Eastman (TE) Process. As a benchmark simulation process, the Tennessee Eastman (TE) process has been widely used for algorithm testing and evaluation over the past several decades. This process consists of five major unit

Figure 2. Scatter plot of the first two variables in different operation modes.

process measurements. Since Gaussian distribution is unimodal and the scatter plot of a bivariate Gaussian distribution should look like an ellipse, it is obvious that the process measurements are not Gaussian. Because of the first-order Markov property in the process, values of l = 1 and h = 2 were chosen to form the past information vectors p and the future information vectors f, so 1609

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research

manipulated variables are measured. A more-detailed description of this process can be found in ref 29. In the present paper, 16 continuous process variables are selected for monitoring purposes, similar to that reported in the work of Ge and Song.24 In the TE benchmark process, six different operation modes can be simulated. In this case study, for the sake of simplicity, process data from three operation modes are incorporated. For model construction purposes, 800 successive data samples have been collected in each operation mode. For testing of fault detection performance, a set of 21 programmed faults was simulated in the first operation mode, which include 7 step faults, 5 random faults, 3 sticking and slow change faults, and 6 unknown process faults.1 Figure 4 illustrates the distribution of the first four process measurements. It is obvious that the process measurements are in three different operation modes. In MCVA modeling, l = 4 and h = 5 was chosen to form the past information vectors (p) and the future information vectors (f), so that the dimensionalities of both p and f are 80. Three Gaussian components are chosen in the mixture model.

Figure 3. Fault detection performance with fault 1 in the numerical process: (a) MCVA, (b) CVA, and (c) GMM.

Table 2. Fault Detection Rates in the Numerical Process MCVA fault 1 fault 2 fault 3

CVA

GMM

BIPs

BIPn

T2s

T2n

BIP

0.993 0.898 0.311

0.862 0.853 0.278

0.545 0.828 0.256

0.180 0.801 0.184

0.768 0.815 0.213

operations: a reactor, a condenser, a compressor, a separator, and a stripper. There are four reactants (identified as A, C, D, and E) and one inert component (identified as B) which are fed into the reactor; then, products (identified as G and H) are formed and, at the same time, a byproduct (identified as F) is also produced. A total of 41 process variables and 12

Figure 4. Scatter plot in different operation modes: (a) the first and second variables, and (b) the third and fourth variables. 1610

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research To test the normal operation monitoring, 500 successive data samples have been collected in each operation mode. Figures 5,

Figure 6. Fault detection performance with normal operation mode 2 in the TE process: (a) MCVA, (b) CVA, and (c) GMM. Figure 5. Fault detection performance with normal operation mode 1 in the TE process: (a) MCVA, (b) CVA, and (c) GMM.

that conventional CVA does not work well in multimode process monitoring. The fault detection performances of MCVA, CVA, and GMM were compared in all of the 21 fault scenarios. The detailed fault detection results with fault 11 and fault 19 are shown in Figures 8 and 9, respectively. Fault 11 is a fault with random variation in the reactor cooling water inlet temperature, while fault 19 is an unknown fault. Figure 8c shows that the monitoring performance of conventional GMM is not good with fault 11, because it fails to identify faulty samples after approximately the 800th sample. Similarly, the GMM-based BIP index has a large missed detection rate with fault 19. Figures 8b and 9b show the performance of CVA. CVA

6, and 7 show the monitoring results in mode 1, mode 2 and mode 3, respectively. One can see that both MCVA and GMM have small false alarm rates. However, the modeling performance of conventional CVA is not very good. Although the false alarm rates of CVA in mode 2 and mode 3 are small, the false alarm rate of CVA in mode 1 is very large. The average false alarm rates with the three operation modes are listed in Table 3. The average false alarm rates of both of the CVA-based statistics are >10%. The large false alarm rate of CVA indicates 1611

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research

Figure 7. Fault detection performance with normal operation mode 3 in the TE process: (a) MCVA, (b) CVA, and (c) GMM.

Figure 8. Fault detection performance with fault 11 in the TE process: (a) MCVA, (b) CVA, and (c) GMM.

Table 3. False Alarm Rates in the Tennessee Eastman (TE) Process MCVA

CVA

BIPs

BIPn

T2s

0.007

0.031

0.157

The fault detection rates of MCVA and GMM are listed in Table 4. Since conventional CVA is not suitable in multimode process monitoring, the monitoring results of CVA are excluded in Table 4. In each fault scenario, the ones with the best detection rates are shown in boldface font. It is shown that the fault detection performance in the MCVA model are much better than that in the GMM-based model. Specifically, the average fault detection rate of the BIPs index is 68.6%, while that of the BIPn index is 62.0%. However, the average fault detection rate of the GMM-based BIP index is only 57.0%.

GMM T2n

BIP

0.144

0

identified almost all the faulty samples; however, the false alarm rate is also very large, so that conventional CVA is not suitable for monitoring multimode processes. In contrast, the proposed MCVA performs much better than GMM and CVA. In MCVA, almost all the faulty samples in fault 11 and fault 19 are identified by both the BIPs and BIPn indices, while the false alarm rate is small. 1612

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research Table 4. Fault Detection Rates in the TE Process MCVA

GMM

fault

BIPs

BIPn

BIP

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

0.999 0.988 0.058 0.019 0.260 1 0.441 0.980 0.064 0.938 0.795 0.999 0.955 0.999 0.093 0.698 0.975 0.905 0.944 0.760 0.538

0.999 0.974 0.044 0.050 0.139 0.991 0.228 0.934 0.033 0.880 0.790 0.961 0.954 1 0.035 0.304 0.978 0.908 0.905 0.664 0.249

0.995 0.971 0.019 0.005 0.239 1 0.414 0.978 0.020 0.749 0.200 0.989 0.946 0.999 0.036 0.355 0.923 0.901 0.219 0.590 0.426

average

0.686

0.620

0.570

Gaussian components, through which the local monitoring statistics are further integrated into global BIP indices. Consequently, the multimode dynamic process can be monitored using the BIP control charts. Two case studies, including a numerical process and the TE benchmark process, are used to demonstrate the suitability and effectiveness of the proposed process monitoring method for multimode dynamic processes. Compared with conventional CVA and GMM methods, the proposed MCVA model performs well and shows its superiority in multimode dynamic process monitoring.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

Figure 9. Fault detection performance with fault 19 in the TE process: (a) MCVA, (b) CVA, and (c) GMM.

ACKNOWLEDGMENTS This work was supported in part by the National Natural Science Foundation of China (NSFC) (No. 61273167), Project National 973 (No. 2012CB720500), and the Fundamental Research Funds for the Central Universities (No. 2013QNA5016).

5. CONCLUSIONS A novel mixture canonical variate analysis model and the corresponding process monitoring approach is successfully developed in this article. Aimed at the complex industrial processes with multiple operating conditions and process dynamics, the proposed method works when the normal operating data follow a multimodal distribution, which can be characterized by a Gaussian mixture model. Canonical variates are obtained in each Gaussian component, and the local monitoring statistics in both the state space and noise space are calculated. With the constructed mixture model, a Bayesian inference procedure is utilized to compute the posterior probabilities of each monitored sample belonging to all the

■

REFERENCES

(1) Chiang, L. H.; Braatz, R. D.; Russell, E. L. Fault Detection and Diagnosis in Industrial Systems; Springer: London, 2001. (2) Ge, Z.; Song, Z. Multivariate Statistical Process Control: Process Monitoring Methods and Applications; Springer: London, 2013. (3) MacGregor, J.; Cinar, A. Monitoring, fault diagnosis, faulttolerant control and optimization: Data driven methods. Comput. Chem. Eng. 2012, 47, 111−120.

1613

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614

Article

Industrial & Engineering Chemistry Research

(26) Bishop, C. M.; Nasrabadi, N. M. Pattern Recognition and Machine Learning, Vol. 1; Springer: New York, 2006. (27) Figueiredo, M. A. T.; Jain, A. K. Unsupervised learning of finite mixture models. IEEE Trans. Pattern Anal. Machine Intell. 2002, 24 (3), 381−396. (28) Feital, T.; Kruger, U.; Dutra, J.; Pinto, J. C.; Lima, E. L. Modeling and performance monitoring of multivariate multimodal processes. AIChE J. 2013, 59 (5), 1557−1569. (29) Downs, J. J.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17 (3), 245−255.

(4) Ge, Z. Q.; Song, Z. H.; Gao, F. R. Review of Recent Research on Data-Based Process Monitoring. Ind. Eng. Chem. Res. 2013, 52 (10), 3543−3562. (5) Kruger, U.; Xie, L. Statistical Monitoring of Complex Multivariate Processes; John Wiley & Sons, Ltd.: West Sussex, U.K., 2012. (6) Ge, Z. Q.; Kruger, U.; Lamont, L.; Xie, L.; Song, Z. H. Fault detection in non-Gaussian vibration systems using dynamic statisticalbased approaches. Mech. Syst. Signal Process. 2010, 24 (8), 2972−2984. (7) Ku, W. F.; Storer, R. H.; Georgakis, C. Disturbance detection and isolation by dynamic principal component analysis. Chemom. Intell. Lab. Syst. 1995, 30 (1), 179−196. (8) Kruger, U.; Zhou, Y. Q.; Irwin, G. W. Improved principal component monitoring of large-scale processes. J. Process Control 2004, 14 (8), 879−888. (9) Ding, S. X. Data-driven design of monitoring and diagnosis systems for dynamic processes: A review of subspace technique based schemes and some recent results. J. Process Control 2014, 24 (2), 431− 449. (10) Schubert, U.; Kruger, U.; Arellano-Garcia, H.; Feital, T.; Wozny, G. Unified model-based fault diagnosis for three industrial application studies. Control Eng. Pract. 2011, 19, 479−490. (11) Simoglou, A.; Martin, E. B.; Morris, A. J. Dynamic multivariate statistical process control using partial least squares and canonical variate analysis. Comput. Chem. Eng. 1999, 23, S277−S280. (12) Juricek, B. C.; Seborg, D. E.; Larimore, W. E. Fault detection using canonical variate analysis. Ind. Eng. Chem. Res. 2004, 43 (2), 458−474. (13) Lee, C.; Choi, S. W.; Lee, I. B. Variable reconstruction and sensor fault identification using canonical variate analysis. J. Process Control 2006, 16 (7), 747−761. (14) Treasure, R. J.; Kruger, U.; Cooper, J. E. Dynamic multivariate statistical process control using subspace identification. J. Process Control 2004, 14, 279−292. (15) Xie, L.; Kruger, U.; Lieftucht, D.; Littler, T.; Chen, Q.; Wang, S. Q. Statistical monitoring of dynamic multivariate processesPart 1. Modeling autocorrelation and cross-correlation. Ind. Eng. Chem. Res. 2006, 45 (5), 1659−1676. (16) Haghani, A.; Ding, S. X.; Esch, J.; Hao, H. Y. Data-Driven Quality Monitoring and Fault Detection for Multimode Nonlinear Processes. In 2012 IEEE 51st Annual Conference on Decision and Control (CDC), Maui, HI, Dec. 10−13, 2012; pp 1239−1244. (17) Choi, S. W.; Lee, I. B. Nonlinear dynamic process monitoring based on dynamic kernel PCA. Chem. Eng. Sci. 2004, 59 (24), 5897− 5908. (18) Zhang, Y. W.; An, J. Y.; Zhang, H. L. Monitoring of time-varying processes using kernel independent component analysis. Chem. Eng. Sci. 2013, 88, 23−32. (19) Hu, Y.; Ma, H. H.; Shi, H. B. Enhanced batch process monitoring using just-in-time-learning based kernel partial least squares. Chemom. Intell. Lab. Syst. 2013, 123, 15−27. (20) Thissen, U.; Swierenga, H.; De Weijer, A. P.; Wehrens, R.; Melssen, W. J.; Buydens, L. M. C. Multivariate statistical process control using mixture modelling. J. Chemom. 2005, 19 (1), 23−31. (21) Yu, J.; Qin, S. J. Multimode process monitoring with Bayesian inference-based finite Gaussian mixture models. AIChE J. 2008, 54 (7), 1811−1829. (22) Ge, Z. Q.; Gao, F. R.; Song, Z. H. Mixture probabilistic PCR model for soft sensing of multimode processes. Chemom. Intell. Lab. Syst. 2011, 105 (1), 91−105. (23) Zhu, J. L.; Ge, Z. Q.; Song, Z. H. Robust modeling of mixture probabilistic principal component analysis and process monitoring application. AIChE J. 2014, 60 (6), 2143−2157. (24) Ge, Z. Q.; Song, Z. H. Maximum-likelihood mixture factor analysis model and its application for process monitoring. Chemom. Intell. Lab. Syst. 2010, 102 (1), 53−61. (25) Negiz, A.; Cinar, A. Statistical monitoring of multivariable dynamic processes with state-space models. AIChE J. 1997, 43 (8), 2002−2020. 1614

DOI: 10.1021/ie503324g Ind. Eng. Chem. Res. 2015, 54, 1605−1614