Gaussian and non-Gaussian Double Subspace Statistical Process

Dec 28, 2014 - Motion process monitoring using optical flow–based principal ... Slow feature analysis based on online feature reordering and feature...
2 downloads 0 Views 6MB Size
Article pubs.acs.org/IECR

Gaussian and non-Gaussian Double Subspace Statistical Process Monitoring Based on Principal Component Analysis and Independent Component Analysis Jian Huang and Xuefeng Yan* 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: This study proposes a new statistical process monitoring method based on variable distribution characteristic (VDSPM) with consideration that variables submit to different distributions in chemical processes and that principal component analysis (PCA) and independent component analysis (ICA) are, respectively, suitable for processing data with Gaussian and nonGaussian distribution. In VDSPM, D-test is first employed to identify the normality of process variables. The process variables under Gaussian distribution are classified into Gaussian subspace and the others belong to non-Gaussian subspace. PCA and ICA models are respectively built for fault detection in Gaussian and non-Gaussian subspaces. Bayesian inference is used to combine the monitoring results of the two subspaces to create a final statistic. The proposed method is applied to a numerical system and to the Tennessee Eastman benchmark process. Results proved that the proposed system outperformed the PCA and ICA methods.

1. INTRODUCTION In modern industries, process monitoring aims to detect and diagnose faults immediately and accurately to ensure production safety and high product quality.1,2 With the development of distributed control systems and data storage technology, vast amounts of process data have been collected and utilized, thus providing a foundation for data-based process monitoring. Traditional multivariate statistical process monitoring methods, such as principal component analysis (PCA), are widely used in modern industries to map high-dimensional correlated process data onto a low-dimensional subspace to decorrelate such data.3−5 The low-dimensional subspace extracts the main variance information that represents the original data. However, the successful PCA application has numerous limitations: (1) the process data should follow a Gaussian distribution (the computation of threshold is beyond this), and (2) the process data should be linearly related. The variables in real industries often exhibit different characteristics (i.e., Gaussian or nonGaussian distribution, linear or nonlinear relation, and so on). Satisfying one or two conditions is generally difficult to achieve. Numerous scholars have recently proposed multiblock or multisubspace theory to group important information with similar characteristics. Thus, the monitoring model can detect faults easily and accurately when abnormal events occur. A linear subspace method was proposed by Ge et al. to solve the nonlinear issue. In this method, the original space is divided into k+1 overlapped subspaces, where k is the number of retained principal components (PCs), and rest 1 is the residual subspace (RS).3 Tong et al. proposed a four-subspace construction monitoring method that divides the original process variables into four subspaces on the basis of their relevance to the PC subspace (PCS) and RS.6 An adaptive multiblock PCA method was developed by Lee et al. This method was applied to the batch process to overcome the problem of changing process conditions by recursively updating the covariance structure.7 Wang et al. © 2014 American Chemical Society

proposed multiblock PCA based on Kullback−Leibler Divergence to group variables with similar characteristics into the same block.8 The relationships among variables are widely known to be complicated (independent, linear, and nonlinear). Moreover, variable distribution is recognized as complex (Gaussian and non-Gaussian). Mining information hiding in the process data is a crucial task in process monitoring. That is, the result of relation recognition determines the monitoring result. Previous scholars focused on different perspectives and proposed different multisubspace or multiblock monitoring methods. No official data-driven multisubspace or multiblock method has been developed. These previously proposed multiblock methods did not consider the variable distribution characteristics. Thus, each block contains both Gaussian and non-Gaussian variables. Moreover, variables with different distribution characteristics are processed by using the same method. Independent component analysis (ICA) is applied to extract non-Gaussian components and to monitor the non-Gaussian process.9,10 However, ICA cannot work well for the Gaussian process. In practice, not all variables follow a Gaussian distribution. The process data contains Gaussian variables and non-Gaussian variables. Numerous extensions of PCA and ICA have been proposed to address the non-Gaussian issue. Liu and Xie proposed non-Gaussian system monitoring that applies PCA to capture the Gaussian and non-Gaussian source signals and ICA to extract the non-Gaussian components from the retained PCs.11 The combination of ICA and PCA was proposed in a twostep strategy that extracts Gaussian and non-Gaussian information for fault detection and diagnosis.12 Zhao et al. applied kernel ICA−PCA, which captured the nonlinear characteristics and Received: Revised: Accepted: Published: 1015

June 24, 2014 December 16, 2014 December 27, 2014 December 28, 2014 DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

where Λ = diag(λ1, λ2,···, λm) is the eigenvalue matrix. λi is in descending order. The first k PCs are retained in the PCS. The others are selected moved to RS. Through the cumulative percent variance (CPV) method, the determination of k is calculated as follows:

explored both nonlinear Gaussian and non-Gaussian features for nonlinear batch processes.13 Zhang proposed a nonlinear process using kernel PCA (KPCA), kernel ICA (KICA), and support vector machine (SVM) for fault detection; this process combines the advantages of both KPCA and KICA.14 The above methods considered both Gaussian and non-Gaussian information on the basis of original data mapping but disregarded the distribution characteristics of original variables. The characteristics combining different original variables of extracted components were analyzed. These methods made the analysis more complicated. Statistical theory presents several methods for testing the normality of variables on the basis of sample data. Chi-square test is the most commonly used for any distribution. This test considers the sample frequency and theoretical frequency of the hypothesis test distribution in the same interval and then constructs a χ2 test statistic.15 A normality test method was proposed to build a statistic based on kurtosis and skewness, which gained acceptance among economists.16,17 Shapiro and Wilk proposed a W test that divides the square of an appropriate linear combination of the sample order statistics by using the usual symmetric estimate of variance.18 The Shapiro−Wilk W test is an omnibus method but fails to work well with a sample size larger than 50. D’Agostin then proposed a D-test that applies the constant ratio of Downton’s linear unbiased estimator of population standard deviation to the sample standard deviation to detect deviations from normality attributed to either skewness or kurtosis.19 This study focuses on the Gaussian and non-Gaussian behavior of variables. The normality test divides the original variables into two subspaces. We propose a novel statistical process monitoring method based on variable distribution characteristic and Bayesian inference (VDSPM). VDSPM divides the original data into two subspaces by testing the normality of every variable. A variable tested as Gaussian is classified into the Gaussian subspace. Otherwise, the variable is classified into the non-Gaussian subspace. After classification, PCA and ICA models can be built in the Gaussian and nonGaussian subspaces, respectively. For online monitoring, the PCA and ICA models are used for fault detection. Bayesian inference combines the results of both models. The VDSPM method considers the normality information, as well as addresses the concurrence issue of Gaussian and non-Gaussian variables. The remainder of this paper is organized as follows: Section 2 briefly reviews PCA, ICA, and D-test. Sections 3 and 4 provide a detailed description of the proposed VDSPM method and its monitoring procedure. Section 5 discusses the applications of the method to both cases, as well as the corresponding results. Section 6 presents the conclusions.

k i=1

(2)

i=1

Over 85% of the information represents the total information. Data matrix X can be expressed as follows: ̂ ̂T + E = X̂ + E (3) X = TP m×k n×k ̂ ̂ where P ∈ R which is the first k columns of P and T ∈ R are the loading and score matrixes, respectively; and E is the residual matrix. With consideration of an observation vector x ∈ Rm×1, two monitoring statistics, T2 and SPE, are constructed by using the following equations to monitor PCS and RS: T2 = x T P̂ Λk −1P̂Tx = t Λk −1t T

(4)

̂ ̂ T )x SPE = x T(I − PP

(5)

where Λk = diag(λ1, λ2,···, λk). The thresholds are based on the assumption that the observations follow a Gaussian distribution given by T2 ≤ T2lim =

k(n − 1) Fk , n − k , α n−k

SPE ≤ SPE lim

(6)

⎡ ⎤1/ h0 2 2 θ C h θ ( − 1) h h α 2 0 ⎥ = θ1⎢ + 1 + 2 0 02 ⎢ ⎥ θ θ 1 1 ⎣ ⎦ (7)

where θi = ∑j m= k + 1 λij (i = 1, 2, 3), h0

= 1− ((2θ1θ3)/(3θ22)), Cα is

a normal deviate with the (1 − α) percentile, and Fk,n−k,α is an F-distribution with k and n − k degrees of freedom with confidence limits α. 2.2. ICA. ICA is a statistical and computational technique used to reveal hidden factors that underlie sets of random variables, measurements, or signals.9,21 ICA is applied to extract nonGaussian components to monitor the non-Gaussian process. Given the process data, X = [x1, x2,···, xm]T ∈ Rm×n with m observations and n samples can be represented as linear combinations of d unknown independent components (ICs) S = [s1, s2,···, sd]T ∈ Rd×n(d ≤ m), which can be expressed as follows: (8)

X = AS + E

where A ∈ Rm×d is the mixing matrix and E is the residual matrix. To estimate the mixing matrix A and the IC matrix S, a demixing matrix W is calculated by the FastICA method to reconstruct matrix Ŝ

2. PRELIMINARIES This section briefly reviews PCA, ICA, and D-test for the proposed VDSPM method. 2.1. PCA. PCA is a dimensional reduction method that transposes the original data onto a lower dimensional subspace, which is a representation of the original space.20 We consider the original data matrix X = [x1, x2,···, xm] ∈ Rn×m, which is scaled to zero mean and unit variance. In this matrix, m and n are the numbers of variables and samples, respectively. The covariance matrix C can be calculated through X, and PCA loading matrix P can be obtained by using the eigenvalue decomposition of C as follows: XTX C= = PΛPT n−1

m

∑ λi /∑ λi × 100% ≥ 85%

Ŝ = WX

(9)

that exhibits an independence property. The first step in FastICA is whitening to eliminate the cross-correlations among the measurement variables.9,21 In the whitening step, eigenvalue decomposition is generally used for the covariance C.22 Given m-dimensional vector x(k), covariance is expressed as C = E(x(k)xT(k)) = UΛUT, where E is expectation. The whitening transformation can be given by

z(k) = Qx(k) −1/2

(10)

where Q = Λ U . Evidently, R = E(z(k)z (k)) is an identity matrix. Equation 10 is expressed as

(1) 1016

T

T

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

The main idea of ICA dimensional reduction is that the observations are linear combinations of the extracted ICs. We must choose the appropriate number of ICs to achieve high variation with less computation synchronously. Generally, L2 norm assumes that the rows of W with the larger sum of squares capture the larger variation of the process.9 The selected b rows of W reconstruct a low-dimensional matrix Wb, and the corresponding selected columns of B reconstruct matrix Bb. Similar to PCA, two statistics I2 and SPE at sample k are constructed as

I2 = s(̂ k)T s(̂ k)

(15)

SPE(k) = e(k)T e(k) = (x(k) − x̂(k))T (x(k) − x̂(k)) (16)

where ŝ(k) = Wbx(k), and the predicted x̂(k) = Q−1BbWbx(k). Control limits are calculated on the basis of kernel density estimation (KDE) because latent variables are not under Gaussian distribution.9 2.3. D-test. The D-test was proposed by D’Agostin for normality test to address the issue of large samples (sample number > 50).19 For n samples, the test statistic was proposed as n

Figure 1. Illustration of VDSPM method.

D=

Table 1. Y Values of Five Variables Y

1

2

3

4

5

1.729

1.960

0.746

1.773

0.373

z(k) = Qx(k) = QAs(k) = Bs(k)

Y=

(12)

Subsequently, the problem of finding a full-rank matrix A is replaced by that of finding an orthogonal matrix B = QA, which is less complicated than A and matrix Q is known. ŝ(k) is inferred by eq 11 as s(̂ k) = BTz(k) = BTQx(k)

)x(i)

n



( n )3 ∑i = 1 (x(i) − x )2

(17)

(D − 0.28209479) n 0.02998598

(18)

3. VARIABLE DISTRIBUTION CHARACTERISTIC OF THE STATISTICAL PROCESS MONITORING SCHEME This section provides a detailed description of the proposed VDSPM scheme. 3.1. Classification of Original Variables. PCA and ICA are dimensional reduction methods used to handle Gaussian and non-Gaussian issues, respectively. The computation of the PCA thresholds is beyond the Gaussian hypothesis. However, not all variables follow Gaussian distribution in real industries. The

(13)

The demixing matrix is expressed as W = BTQ

n+1 2

where n samples range in ascending order as x(1) ≤ x(2) ≤ ··· ≤ x(n), and x̅ = (1/n)∑i n= 1x(i). On the basis of the given significance level α and number of samples n, the α/2 and ((1 − α)/2) quantiles of Y, Yα/2, and Y((1−α)/2) can be acquired. When Y < Yα/2 or Y > Y((1−α)/2), the normality hypothesis is rejected; or accepted.

(11)

where B is an orthogonal matrix, which can be proven by E(z(k)z T(k)) = BE(s(k)s T(k))BT = BBT = I

(

∑i = 1 i −

(14)

in accordance with eqs 9 and 13.

Figure 2. Probability plots of original variables: (a) variable 5; (b) variable 4. 1017

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research Table 2. Detection Rates of Numerical Example PCA

ICA

DPCA

DICA

VDSPM

fault no.

T2

SPE

I2

SPE

T2

SPE

I2

SPE

BIC

1 2

0.698 0.928

0.566 0.833

0.826 0.485

0.736 0.996

0.799 0.994

0.439 0.560

0.638 0.921

0.828 0.898

0.900 0.990

Figure 3. Monitoring results of fault 1: (a) PCA; (b) ICA; (c) DPCA; (d) DICA; (e) VDSPM.

The original data matrix X = [x1, x2,···, xm]T ∈ Rn×m with m observations and n samples is considered. D-test examines the normality of every observation xi (i = 1, 2,···, m) by calculating Y in eq 18. Extremely small or large Y will result in a non-Gaussian judgment. When Y < Yα/2 or Y > Y(1−(α/2)), xi belongs to the nonGaussian block; otherwise it belongs to the Gaussian block. After

division and management of Gaussian and non-Gaussian information appear particularly important. Dividing into subspaces aims at obtaining similar information (Gaussianity or non-Gaussianity in this paper) and grouping it together. Moreover, the VDSPM method is proposed on the basis of distribution characteristics. 1018

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

Figure 4. Monitoring results of fault 2: (a) PCA; (b) ICA; (c) DPCA; (d) DICA; (e) VDSPM.

the classifying operation, the original data X can be divided into two blocks as

X = [XPCA , XICA]

The PC number is calculated by using the CPV method. The thresholds of T2 and SPE can also be calculated by using eqs 6 and 7. An ICA model is built on the non-Gaussian block XICA, as shown in the following equation:

(19)

XTICA = A ICA SICA + E ICA

where XPCA and XICA represent the Gaussian block and nonGaussian block, respectively. Scaling XPCA with zero mean and unit variance allows the following PCA model: T XPCA = TPCA PPCA + E PCA

(21)

The thresholds of I2 and SPE can be calculated by KDE. 3.2. Bayesian Inference. Each model built in every subspace has individual thresholds. Thus, every subspace can be easily monitored. However, combining the results to obtain a final decision is difficult. In this study, a Bayesian inference is employed to

(20) 1019

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

Figure 5. Contribution plots of faults 1: (a) VDSPM-ICA; (b) ICA; (c) DICA.

Figure 6. Contribution plots of faults 2: (a) VDSPM-PCA; (b) PCA; (c) DPCA.

perform the combination operation. This inference is based on the probabilistic approach. The Bayesian inference is similar to that of Ge et al.,3 who provide more details. Current monitoring sample x is divided into two blocks in accordance with the classification step. Then, the fault probability to T2 of xPCA is expressed as pT2(x |F ) pT2F PCA pT2(F|x ) = PCA p T2 x (22)

represents abnormal conditions. pT2N and pT2F are previously set as α and 1 − α, respectively, where α is the confidence level. pT2(xPCA|N) and pT2(xPCA|F) are calculated by the following equations: pT2(x

PCA |N )

PCA

and p T2 x

pT2(x PCA

= pT2(x

PCA |N )

pT2N + pT2(x

PCA |F )

p T2 F

(23)

PCA |F )

⎛ T 2⎞ x = exp⎜⎜ − PCA2 ⎟⎟ ⎝ Tlim ⎠

(24)

⎛ T 2 ⎞ = exp⎜⎜ − lim 2 ⎟⎟ ⎝ TxPCA ⎠

(25)

2

where Tlim is the threshold of Gaussian subspace. Similar to the above equations, pSPE(F|xPCA), pI2(F|xICA), and pSPE(F|xICA) can be calculated. The final statistic, which combines

where xPCA is the current monitoring sample x selected in Gaussian subspace, N represents normal conditions, and F 1020

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research Table 3. Y Values of 52 Variables variable no.

Y

variable no.

Y

variable no.

Y

variable no.

Y

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

−0.924 3.191 2.732 0.503 −0.788 0.357 2.599 0.512 −3.467 −4.444 1.258 −1.411 2.443

14 15 16 17 18 19 20 21 22 23 24 25 26

−1.075 0.484 2.055 0.231 −11.354 −9.423 −2.521 1.733 0.198 −2.372 −0.348 0.860 2.328

27 28 29 30 31 32 33 34 35 36 37 38 39

2.276 −1.460 0.655 −2.120 1.874 −2.166 1.259 3.279 −0.371 −2.273 0.485 −1.047 0.671

40 41 42 43 44 45 46 47 48 49 50 51 52

2.137 0.456 1.473 1.621 −1.283 −1.563 1.954 −2.875 −1.411 0.484 −9.445 0.909 0.302

Table 4. Variables Division in Each Subspace variables no.

Gaussian subspace

Non-Gaussian subspace

1,4,5,6,8,11,12,14,15,17,21,22,24,25,28,29,30,33,35,37,38,39,41,42,43,44,45,48,49,51,52

2,3,7,9,10,13,16,18,19,20,23,26,27,31,32,34,36,40,46,47,50

Table 5. Missing Detection Rates of TE Process PCA

ICA

DPCA

DICA

ICA-PCA

VDSPM

faultno.

T2

SPE

I2

SPE

T2

SPE

I2

SPE

I2

T2

SPE

BIC

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

0.008 0.020 0.998 0.456 0.775 0.008 0 0.025 0.994 0.666 0.794 0.029 0.060 0 0.988 0.843 0.259 0.113 0.861 0.685 0.736

0.003 0.014 0.991 0.038 0.746 0 0 0.025 0.981 0.659 0.356 0.025 0.045 0.011 0.973 0.755 0.108 0.101 0.720 0.398 0.570

0.003 0.011 0.973 0.225 0 0 0 0.031 0.989 0.254 0.303 0.003 0.046 0 0.960 0.196 0.121 0.093 0.352 0.221 0.575

0.003 0.014 0.923 0.113 0 0 0 0.013 0.969 0.245 0.446 0.003 0.045 0.001 0.848 0.251 0.096 0.104 0.626 0.258 0.510

0.007 0.014 0.971 0.879 0.732 0.009 0 0.026 0.969 0.636 0.801 0.007 0.05 0.001 0.959 0.824 0.222 0.111 0.801 0.639 0.56

0 0.011 0.756 0 0.471 0 0 0.026 0.762 0.429 0.13 0.03 0.042 0.001 0.781 0.414 0.025 0.1 0.425 0.305 0.401

0.001 0.012 0.98 0.001 0 0 0 0.021 0.974 0.101 0.205 0 0.045 0 0.887 0.135 0.027 0.092 0.091 0.107 0.556

0.002 0.017 0.974 0.005 0 0 0 0.021 0.989 0.096 0.245 0.001 0.046 0 0.919 0.249 0.044 0.094 0.097 0.186 0.626

0.005 0.018 0.996 0.032 0 0 0.050 0.028 0.994 0.212 0.712 0.024 0.016 0.042 0.412 0.185 0.206 0.100 0.346 0.412 0.622

0.003 0.015 0.993 0.020 0 0 0.035 0.022 0.993 0.156 0.526 0.005 0.011 0.013 0.268 0.146 0.107 0.085 0.238 0.247 0.549

0 0.008 0.985 0 0 0 0 0.011 0.991 0.162 0.182 0 0 0 0.196 0.125 0.195 0.098 0.208 0.306 0.526

0.002 0.011 0.877 0 0.565 0 0 0.018 0.871 0.135 0.201 0.005 0.036 0 0.648 0.119 0.041 0.081 0.144 0.166 0.448

4.1. Fault Detection. Figure 1 illustrates the VDSPM method. The detailed steps are shown below. (1). Offline Modeling. Step 1: Training data X are collected from the process under normal conditions. Step 2: The Y value of every column in X is calculated by using eqs 17 and 18. Step 3: X is divided into two blocks (Gaussian and nonGaussian) by determining whether the Y value is within the normality range. Step 4: The PCA model is built in the Gaussian subspace, meanwhile the ICA model is constructed in the non-Gaussian subspace. Step 5: The control limits of every subspace, Tlim2 and SPElim in Gaussian subspace, and Ilim2 and SPElim in non-Gaussian subspace are obtained.

the results in different subspaces, is determined by using the following weighted equation: BIC =

p T 2 (F |x

2 PCA )

p T 2 (F |x

+ pSPE(F|x

PCA )

2 PCA )

+ pSPE(F|x

PCA )

2

+ p I 2(F|x

ICA )

+ p I 2(F|x

ICA )

+ pSPE(F|x + pSPE(F|x

2 ICA )

ICA )

(26)

The BIC threshold is equal to 1 − α. In particular, when the BIC value is over 1 − α, the process can be abnormal; otherwise, the process is normal.

4. PROCESS MONITORING This section presents the VDSPM procedure for process monitoring. 1021

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research Table 6. False Alarm Rates/Detection Delays of TE Process PCA

ICA

DPCA

DICA

VDSPM

fault no.

T2

SPE

I2

SPE

T2

SPE

I2

SPE

BIC

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

0.013/21 0.019/51 0.013/− 0.019/− 0.019/48 0/30 0/3 0/69 0.038/− 0/288 0.006/18 0.019/66 0.013/147 0/12 0.006/− 0.063/936 0.013/87 0.006/279 0.006/− 0.006/261 0.006/1689

0.069/9 0.081/36 0.094/− 0.069/9 0.069/3 0.056/3 0.031/3 0.081/60 0.094/− 0.044/147 0.075/33 0.100/24 0.056/111 0.094/3 0.031/− 0.075/591 0.113/75 0.119/252 0.100/− 0.056/261 0.213/855

0.025/9 0.019/33 0.025/− 0.025/3 0.025/3 0.025/3 0.031/3 0.019/60 0.031/− 0.025/81 0.025/33 0.031/9 0.019/111 0.025/6 0.019/− 0.025/33 0.025/66 0.025/57 0.025/36 0.031/201 0.025/801

0.019/9 0.019/45 0.031/− 0.025/9 0.025/6 0.025/3 0.025/3 0.025/21 0.019/− 0.019/69 0.025/18 0.025/9 0.031/111 0.019/6 0.019/− 0.025/42 0.025/75 0.038/255 0.025/30 0.019/216 0.025/771

0/18 0.013/48 0.013/− 0/− 0/6 0/33 0.006/3 0/69 0.025/− 0/303 0.006/585 0.013/9 0/135 0.006/18 0.006/− 0.075/597 0/84 0.013/279 0/− 0/267 0.006/1566

0.088/15 0.100/39 0.044− 0.101/9 0.101/6 0.050/3 0.057/3 0.044/63 0.101/− 0.031/150 0.069/21 0.088/24 0.069/120 0.101/3 0.037/− 0.088/588 0.063/72 0.069/252 0.075/246 0.069/252 0.166/858

0.019/3 0.031/36 0.019/− 0.019/0 0.019/0 0.031/0 0.031/0 0.019/51 0.025/− 0.031/60 0.025/15 0.019/0 0.025/60 0.025/0 0.031/− 0.025/24 0.025/57 0.025/54 0.031/30 0.019/195 0.031/750

0.031/6 0.019/42 0.025/− 0.031/6 0.031/0 0.025/0 0.031/0 0.025/30 0.025/− 0.013/60 0.025/18 0.031/3 0.025/78 0.038/0 0.025/− 0.019/27 0.025/57 0.013/228 0.025/12 0.031/210 0.025/861

0.025/6 0.038/24 0.025/− 0.031/0 0.031/0 0.038/0 0.025/0 0.056/45 0.056/− 0.031/69 0.050/15 0.031/6 0.056/60 0.069/0 0.038/− 0.031/18 0.031/51 0.050/33 0.063/6 0.025/186 0.044/624

r

(2). Online Monitoring. Step 1: Current process data x are collected. x is divided into two blocks on the basis of the classification results in offline modeling step 3. Step 2: Statistics T2 and SPE, as well as I2 and SPE, are respectively calculated in Gaussian subspace and non-Gaussian subspace. Step 3: The monitoring results are combined by Bayesian inference. Subsequently, a new statistic BIC is established. The BIC control limit is 1 − α. 4.2. Fault Diagnosis. The monitoring charts cannot indicate which variable causes the fault. Once a fault is detected, the root cause of such fault must be diagnosed. In traditional PCA or ICA method, the contribution plot is widely used.9,23,24 To build the contribution plot of each subspace, we first determine which subspace is out of control. We have three possible situations: (1) a fault occurs in Gaussian subspace; (2) a fault occurs in nonGaussian subspace; and (3) both Gaussian and non-Gaussian subspaces are out of control. (1). A Fault Occurs in Gaussian Subspace. The VDSPMPCA contribution plot is generated when the Gaussian subspace detects the fault. In the scaled monitoring sample, xPCA = [x1,PCA,x2,PCA,···,xmPCA,PCA], where mPCA is the number of variables in Gaussian subspace. When the T2 statistic is over the threshold Tlim2, out-of-control PC is decided by t 2(i)/λ(i) > Tlim 2 /k

cont j =

t (i ) p xj ,PCA λ (i ) i , j

(29)

i=1

When the SPE statistic is over the threshold SPElim, the contribution of jth variable xj, PCA to SPE is calculated by the following: Q j = (xj ,PCA − xĵ ,PCA )2

(30)

where x̂j,PCA is the jth element of x̂PCA. x̂PCA can be expressed as

̂ ̂ x̂PCA = xPCAPP

T

(31)

(2). A Fault Occurs in the non-Gaussian Subspace. When the I2 statistic is over the threshold I2lim, for mICA-dimensional sample xICA = [x1,ICA, x2,ICA, ···, xmICA,ICA] with q out-of-control ICs, the total contribution of jth variable xj,ICA is calculated as q

cont j =

∑ s(i)xj ,ICA

(32)

i=1

where s(i) is the ith IC and the out-of-control IC is determined by s 2(i) > Ilim 2/b

(33)

When the SPE statistic is over the threshold SPElim, the contribution of jth variable to SPE is calculated by the following:

(27)

Q j = (xj ,ICA − xĵ ,ICA )2

where t(i) is the ith PC, and λ(i) is the corresponding eigenvalue. The number of out-of-control PCs is assumed as r. The contribution of the jth variable xj,PCA to the ith PC can be expressed as follows: cont i , j =

∑ conti ,j

(34) −1

where x̂j,ICAis the jth element of x̂ICA = Q BbWbxICA. A VDSPMICA plot is then created. (3). Both Gaussian and non-Gaussian Subspaces Are out of Control. The PCA and ICA contribution plots should be given when the Gaussian and non-Gaussian subspaces are out of control.

(28)

5. APPLICATION In this section, the proposed VDSPM method is applied to a numerical system and to the Tennessee Eastman (TE)

where pi, j is the (i, j)th element of loading matrix P. The total contribution of jth variable is computed as 1022

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

Figure 7. Monitoring results of fault 10: (a) PCA; (b) ICA; (c) VDSPM.

Figure 8. Monitoring results of fault 11: (a) PCA; (b) ICA; (c) VDSPM.

benchmark system. Monitoring results are presented. Moreover, different methods are compared. 5.1. Numerical System. A simple numerical multivariate system is used to test the performance of the VDSPM method. This system was generated by Ku et al.25 and modified by Lee at al.9

where x1 and x2 are input vectors, in which each element of x1 is uniformly distributed over the interval (−2,2). Each element of x2 is Gaussian distributed with zero mean and variance of 0.25. v is a random noise vector and each element has zero mean and variance of 0.1. The input u and output y are measurable, but t and x are not. To conduct analysis, 960 samples of five variables (y1, y2 y3 u1, u2) are obtained. The normality test of all five variables is implemented by D-test. Table 1 presents the detailed Y values of five variables. Significance level α is chosen as 0.05, whereas Yα/2 and Y1−α/2 are equal to −2.16 and 1.75, respectively. The test result shows that variables 1, 3, and 5 are Gaussian distributed, but the others are not. Figure 2 presents the probability plots of variables 5 and 4. The plots reveal that variable 4 seriously deviates from the dotted line representing the Gaussian distribution. Subsequently, variable 4 is judged to be

⎡ 0.118 −0.191 0.287 ⎤ ⎡1 2 ⎤ ⎥ ⎢ ⎢ ⎥ t(i) = ⎢ 0.847 0.264 0.943 ⎥t(i − 1) + ⎢ 3 −4 ⎥u ⎢⎣−0.333 0.514 −0.217 ⎥⎦ ⎣ −2 1 ⎦ (i − 1) y(i) = t(i) + v(i) ⎡ 0.193 0.689 ⎤ ⎡ 0.811 −0.226 ⎤ u(i) = ⎢ ⎥x(i − 1) ⎥u(i − 1) + ⎢ ⎣ 0.477 0.415 ⎦ ⎣−0.320 −0.749 ⎦ 1023

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

Figure 9. Monitoring results of fault 19: (a) PCA; (b) ICA; (c) VDSPM.

Figure 10. Monitoring results of fault 20: (a) PCA; (b) ICA; (c) VDSPM.

non-Gaussian. The probability plots and results of normality test complement each other. Thus, XPCA contains variables 1, 3, and 5. Moreover, XICA contains variables 2 and 4. The simulated faults are given by the following: Fault 1, a step change of x1 by 3 is introduced at sample 161; Fault 2, a step change of x2 by 1.3 is introduced at sample 161. Table 2 shows the detection rates of PCA, ICA, dynamic PCA(DPCA),25 DICA,26 and VDSPM methods. All confidence limits are set as 97%. Figure 3 shows the monitoring results of fault 1. We compare the monitoring performances among PCA shown in Figure 3a, ICA shown in Figure 3b, DPCA shown in Figure 3c, DICA shown in Figure 3d, and VDSPM shown in Figure 3e. VDSPM can detect more fault data than the other four methods. Fault detection rates for PCA, ICA, DPCA, and DICA are lower than that for VDSPM. Figure 4 presents the monitoring

results of fault 2. Both monitoring charts illustrate that the proposed VDSPM method outperforms the original PCA, ICA, DPCA, and DICA methods in handling a concurrence issue of Gaussian and non-Gaussian variables. Regardless of whether the fault is generated by non-Gaussian input (x1) or Gaussian input (x2), the proposed VDSPM can detect faults efficiently. Hence, VDSPM combines the advantages of PCA and ICA. For fault 1, the contribution plots of VDSPM-ICA, ICA, and DICA are presented in Figure 5 panels a, b, and c, respectively. Only non-Gaussian subspace detects the fault in the proposed VDSPM method, thereby generating the VDSPM-ICA plot. Compared with the traditional ICA plot and DICA plot, the VDSPM-ICA plot can diagnose the fault accurately while decreasing the number of fault variables. Figure 6 illustrates the 1024

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

Figure 11. Contribution plots of faults 10 and 11: (a) VDSPM-ICA for fault 10; (b) ICA for fault 10; (c) VDSPM-PCA for fault 11; (d) PCA for fault 11; (e) VDSPM-ICA for fault 11; (f) ICA for fault 11.

root cause of fault 2. The Gaussian subspace T2 statistic can detect fault 2. The PCA T2 contribution plot and DPCA T2 contribution plot are given in Figure 6 panels b and c, respectively. The VDSPM-PCA plot has a similar diagnoses as the traditional PCA and DPCA methods. The VDSPM plot can diagnose the faults with less fault variables compared to both contribution plots. 5.2. TE Benchmark System. The TE benchmark process was developed by Downs and Vogel.27 The second structure of the process was described by Lyman and Georgakis.28 A detailed description of the TE process is provided in the Supporting Information. The control system, measurements, and faults for TE process are presented in Figure 1, Table 1, and Table 2 of the Supporting Information, respectively. The first 52 variables are selected to build the process data set. Additionally, 960 samples are created for each fault by introducing a fault from sample 161.

The proposed VDSPM, PCA, and ICA methods are built for comparative analysis. First, the D-test is used for subspace division. Table 3 presents the detailed Y values of 52 variables. Significance level α is set as 0.05, whereas Y1−α and Y1−α/2 are equal to −2.16 and 1.75, respectively. On the basis of the division rules, the original variables are divided into two blocks. Table 4 shows the detailed classification results. Table 5 presents the missing detecting rates of the PCA, ICA, DPCA, DICA, ICA-PCA,12 and VDSPM. The VDSPM method is seen to perform as well as the other methods for faults 1, 2, 4, 6, 7, 8, 12, 13, 14, and 18. This comparable performance is caused by the large fault magnitude. Small faults, such as faults 3, 9, and 15, are not detected by VDSPM, PCA, and ICA. VDSPM shows good performance when handling faults 10, 11, 16, 17, 19, and 20. These faults have low missing detection rates. Table 6 illustrates the false alarm rates and detection delays 1025

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research

distribution. The D-test is used to divide original data into two blocks automatically by testing the normality of every variable. After the division step, PCA and ICA models are built in two subspaces, and Bayesian inference is employed to combine the monitoring results of the two subspaces. VDSPM is applied to monitor the TE process. The simulation results show that the VDSPM method outperforms the PCA and ICA methods. This study provides a solution to the concurrence issue of Gaussian and non-Gaussian variables. However, nonlinear and dynamic processes are not considered. Further research should focus on extending VDSPM monitoring.

of 21 faults. The false alarm rates of VDSPM are similar to those of the other four methods, thus demonstrating efficient performance when detecting normal data. VDSPM detection delays are less than those of other methods, such as those in faults 2, 13, 16, 17, 18, 19, 20, and 21. VDSPM can divide Gaussian and nonGaussian variables into two blocks for each subspace, as well as build the corresponding monitoring model. These characteristics facilitate efficient fault detection for both Gaussian and nonGaussian variables. Further analyses of faults 10, 11, 19, and 20 prove that the proposed VDSPM outperforms the traditional PCA and ICA methods. For all simulations, 97% confidence level is adopted. Case Study on Fault 10. In fault 10, a random temperature change of stream 4 results in a variation in stripper pressure. Figure 7 shows the monitoring results of three methods. Neither PCA-T2 nor PCA-SPE can detect the fault. The ICA method performs better than PCA and shows low missing detection rates of 0.254 for I2 and 0.245 for SPE. The proposed VDSPM can detect more fault samples than ICA. Case Study on Fault 11. Fault 11 is a random variation. Moreover, large oscillations in the reactor cooling water inlet temperature result in a fluctuation in reactor temperature. Figure 8 shows the performance monitoring of fault 11 using PCA, ICA, and VDSPM methods. Both PCA and ICA can detect the fault. However, the missing detection rates are higher than VDSPM. Therefore PCA and ICA regard more abnormal samples as normal than VDSPM, which may cause accidents in real industries. The VDSPM method is more reliable than PCA and ICA. Case Study on Faults 19 and 20. In the TE process, faults 19 and 20 are unknown faults. Figures 9 and 10 present the monitoring results of faults 19 and 20, respectively. Both figures show that the proposed VDSPM can detect faults efficiently and has missing detection rates lower than those of PCA and ICA. Contribution Plots for Faults 10 and 11. Once a fault is detected, the root cause of the fault must be diagnosed. For fault 10, a random temperature change in stream 4 resulted in a variation in stripper pressure. Based on the primary knowledge of the TE process, it is determined that fault 10 mainly affects tripper temperature (variable 18), stripper steam flow (variable 19), and stripper steam valve (variable 50). Figure 11a shows that only non-Gaussian subspace I2 statistic can detect the fault, which illustrates that variables 18, 19, and 50 make more contributions than other variables. Compared with the traditional ICA I2 contribution plot presented in Figure 11b, the VDSPM-ICA plot can diagnose the root cause and contains less fault variables than ICA. Both Gaussian and non-Gaussian subspaces can detect the fault for fault 11. Hence, the two VDSPM contribution plots are built. In the VDSPM-PCA plot shown in Figure 11c, variable 51 (condenser cooling water flow) is determined as the outstanding variable. Similarly, variable 9 (reactor temperature) is a significant indicator in the VDSPM-ICA plot, as presented in Figure 11e. PCA and ICA contribution plots provide the same diagnosis conclusions as VDSPM. That is, the VDSPM method can diagnose faults accurately.



ASSOCIATED CONTENT

S Supporting Information *

Detailed information on the measured variables, simulation faults, and the control scheme for the TE process. This material 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. Address: P.O. Box 293, MeiLong Road No. 130, Shanghai 200237, P. R. China. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors gratefully acknowledge the support from the following foundations: 973 project of China (2013CB733605), National Natural Science Foundation of China (21176073), and the Fundamental Research Funds for the Central Universities.



REFERENCES

(1) Qin, S. J. Statistical process monitoring: Basics and beyond. J. Chemom. 2003, 17, 480−502. (2) Kresta, J. V.; MacGregor, J. F.; Marlin, T. E. Multivariate statistical monitoring of process operating performance. Can. J. Chem. Eng. 1991, 69, 35−47. (3) Ge, Z. Q.; Zhang, M. G.; Song, Z. H. Nonlinear process monitoring based on linear subspace and Bayesian inference. J. Process Contr. 2010, 20, 676−688. (4) Ge, Z. Q.; Yang, C. J.; Song, Z. H. Improved kernel PCA-based monitoring approach for nonlinear processes. Chem. Eng. Sci. 2009, 64, 2245−2255. (5) Kano, M.; Hasebe, S.; Hashimoto, I.; Ohno, H. A new multivariate statistical process monitoring method using principal component analysis. Comput. Chem. Eng. 2001, 25, 1103−1113. (6) Tong, C. D.; Song, Y.; Yan, X. F. Distributed statistical process monitoring based on four-subspace construction and bayesian inference. Ind. Eng. Chem. Res. 2013, 52, 9897−9907. (7) Lee, D. S.; Vanrolleghem, P. A. Monitoring of a sequencing batch reactor using adaptive multiblock principal component analysis. Biotechnol. Bioeng. 2003, 82, 489−497. (8) Wang, B.; Jiang, Q.; Yan, X. Fault detection and identification using a Kullback−Leibler divergence based multi-block principal component analysis and Bayesian inference. Korean J. Chem. Eng. 2014, 1−14. (9) Lee, J. M.; Yoo, C. K.; Lee, I. B. Statistical process monitoring with independent component analysis. J. Process Control 2004, 14, 467−485. (10) Jiang, Q. C.; Yan, X. F.; Tong, C. D. Double-weighted independent component analysis for non-Gaussian chemical process monitoring. Ind. Eng. Chem. Res. 2013, 52, 14396−14405. (11) Liu, X. Q.; Xie, L.; Kruger, U.; Littler, T.; Wang, S. Q. Statisticalbased monitoring of multivariate non-Gaussian systems. AIChE J. 2008, 54, 2379−2391.

6. CONCLUSION This study proposed a new chemical process monitoring method called VDSPM to address the concurrence issue of Gaussian and non-Gaussian variables. When faced with the concurrence problem, traditional PCA may cause a high missing detection rate because the computation of PCA monitoring statistics is based on the assumption that variables follow a Gaussian 1026

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027

Article

Industrial & Engineering Chemistry Research (12) Ge, Z. Q.; Song, Z. H. Process monitoring based on independent component analysis-principal component analysis (ICA-PCA) and similarity factors. Ind. Eng. Chem. Res. 2007, 46, 2054−2063. (13) Zhao, C. H.; Gao, F. R.; Wang, F. L. Nonlinear batch process monitoring using phase-based kernel-independent component analysis−principal component analysis (KICA−PCA). Ind. Eng. Chem. Res. 2009, 48, 9163−9174. (14) Zhang, Y. W. Enhanced statistical analysis of nonlinear processes using KPCA, KICA and SVM. Chem. Eng. Sci. 2009, 64, 801−811. (15) Gumbel, E. J. On the reliability of classical chi-square test. Ann. Math. Stat. 1943, 14, 253−263. (16) D’Agostino, R. B.; Pearson, E. S. Tests for departure from normality. Empirical results for the distributions of b2 and √b1. Biometrika 1973, 60, 613−622. (17) Jarque, CarlosM.; Bera, AnilK. Efficient tests for normality, homoscedasticity and serial independence of regression residuals. Econ. Lett. 1980, 6, 255−259. (18) Shapiro, S. S.; Wilk, M. B. An analysis of variance test for normality (complete samples). Biometrika 1965, 52, 591−611. (19) D’Agostino, R. B. An omnibus test of normality for moderate and large size samples. Biometrika 1971, 58, 341−348. (20) Chiang, L. H.; Braatz, R. D.; Russell, E. L. Fault Detection and Diagnosis in Industrial Systems; Springer: New York, 2001. (21) Hyvärinen, A.; Oja, E. Independent component analysis: algorithms and applications. Neural Networks 2000, 13, 411−430. (22) Hyvarinen, A. Survey on independent component analysis. Neural Comput. Surveys 1999, 2, 94−128. (23) Westerhuis, J. A.; Gurden, S. P.; Smilde, A. K. Generalized contribution plots in multivariate statistical process monitoring. Chemom. Intell. Lab. 2000, 51, 95−114. (24) Miller, P.; Swanson, R. E.; Heckler, C. E. Contribution plots: A missing link in multivariate quality control. Appl. Math. Comput. Sci. 1998, 8, 775−792. (25) Ku, W.; Storer, R. H.; Georgakis, C. Disturbance detection and isolation by dynamic principal component analysis. Chemom. Intell. Lab. 1995, 30, 179−196. (26) Lee, J. M.; Yoo, C.; Lee, I. B. Statistical monitoring of dynamic processes based on dynamic independent component analysis. Chem. Eng. Sci. 2004, 59, 2995−3006. (27) Downs, J. J.; Vogel, E. F. A plant-wide industrial process control problem. Comput. Chem. Eng. 1993, 17, 245−255. (28) Lyman, P. R.; Georgakis, C. Plant-wide control of the Tennessee Eastman problem. Comput. Chem. Eng. 1995, 19, 321−331.

1027

DOI: 10.1021/ie5025358 Ind. Eng. Chem. Res. 2015, 54, 1015−1027