- Email: [email protected]
Adaptive Local Outlier Probability for Dynamic Process Monitoring
Chinese Journal of Chemical Engineering 22 (2014) 820–827
Contents lists available at ScienceDirect
Chinese Journal of Chemical Engineering journal homepage: www.elsevier.com/locate/CJCHE
Process Monitor
Adaptive Local Outlier Probability for Dynamic Process Monitoring☆ Yuxin Ma, Hongbo Shi ⁎, Mengling Wang Key Laboratory of Advanced Control and Optimization for Chemical Processes of Ministry of Education, East China University of Science and Technology, Shanghai 200237, China
a r t i c l e
i n f o
Article history: Received 24 December 2013 Received in revised form 28 January 2014 Accepted 7 February 2014 Available online 19 June 2014 Keywords: Time-varying Complex data distribution Local outlier probability Multi-mode Fault detection
a b s t r a c t Complex industrial processes often have multiple operating modes and present time-varying behavior. The data in one mode may follow specific Gaussian or non-Gaussian distributions. In this paper, a numerically efficient moving window local outlier probability algorithm is proposed. Its key feature is the capability to handle complex data distributions and incursive operating condition changes including slow dynamic variations and instant mode shifts. First, a two-step adaption approach is introduced and some designed updating rules are applied to keep the monitoring model up-to-date. Then, a semi-supervised monitoring strategy is developed with an updating switch rule to deal with mode changes. Based on local probability models, the algorithm has a superior ability in detecting faulty conditions and fast adapting to slow variations and new operating modes. Finally, the utility of the proposed method is demonstrated with a numerical example and a non-isothermal continuous stirred tank reactor. © 2014 Chemical Industry and Engineering Society of China, and Chemical Industry Press. All rights reserved.
1. Introduction In industrial processes, operating conditions are usually affected by some slow variations denoted as time-varying characteristics, caused by some dynamic behavior such as seasonal fluctuation, catalyst deactivation, equipment aging, sensor or process drifting, preventive maintenance and cleaning [1]. Generally, effects of the time-varying behavior on the mean and covariance of variables cannot be neglected, so there may be many false alarms if conventional multivariate statistical process monitoring (MSPM) methods are applied directly [2]. In order to maintain process efficiency for a long period of time, numerous adaptive methods have been developed. Recursive MSPM methods and methods based on the moving window strategy are two alternative widely used approaches [3,4]. Multimodality is another important feature of industrial processes due to changes of market demands, alternations of feedstock or variations of manufacturing strategy. The difference between the characteristics of nearby operating conditions is always significant, so intensive studies have been carried out with either multiple local models or a single global model [5,6]. While it is more practical
☆ Supported by the National Natural Science Foundation of China (61374140), Shanghai Postdoctoral Sustentation Fund (12R21412600), the Fundamental Research Funds for the Central Universities (WH1214039), and Shanghai Pujiang Program (12PJ1402200). ⁎ Corresponding author. E-mail address: [email protected] (H. Shi).
to accommodate the time-varying behavior and multimode features together. The developed methods can be divided into two categories. One is the adaptive clustering methods. Teppola et al. [7] applied adaptive fuzzy C-means algorithms on the score values of principle component analysis (PCA) to monitor a wastewater treatment plant. Liu [8] used an adaptive Takagi-Sugeno fuzzy model on PCA subspace to model a large scale nonlinear system containing many operating regions. Since PCA is used as a preprocessing tool, monitoring results of these two methods more or less depend on and be restricted by the capability of PCA. Petković et al. [9] designed an on-line adaptive clustering method utilizing a generalized information potential. Although previously unseen functioning modes can be included by introducing an adaptive expert system, the method suffers from a non negligible detection delay. The other category is adaptive statistical methods. Improved recursive algorithms based on recursive PCA or the signed digraph were proposed with some if-then rules to distinguish process condition changes from disturbances [10–12]. Ge and Song [13] introduced the just-in-timelearning strategy to the modeling procedure of local least squares support vector regression and the residuals between the real output and the predicted one was analyzed by a two-step information extraction strategy. Xie and Shi [14] and Yu [15] developed two different dynamic fashions of Gaussian mixture model (GMM) separately based on the moving window strategy and a particle filter resampling method. The problem of complex data distributions in time-varying and multimode processes has scarcely been addressed. Although the moving window strategy has been proven to be effective, it still encounters some limitations when incorporated with statistical methods such as
http://dx.doi.org/10.1016/j.cjche.2014.05.015 1004-9541/© 2014 Chemical Industry and Engineering Society of China, and Chemical Industry Press. All rights reserved.
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
PCA, partial least squares (PLS) or GMM. Since the variables of an industrial process may satisfy specified Gaussian or non-Gaussian distributions, and high order statistics are usually helpful to reveal more information from the data [16–18], adaptive monitoring algorithms should be developed, which can explore both Gaussianity and nonGaussianity of process data. Local outlier probability (LoOP) is an unsupervised data mining technique proposed for outlier detection [19]. It combines the idea of local, density-based outlier scoring with a probabilistic, statistically-oriented approach, and assigns the probability of being an outlier to all data records. Since a normalization procedure is included, LoOP is independent of any specific data distribution. Therefore, a combination of LoOP and moving window strategy should be potential to tackle these problems. The main contribution of this paper is to propose a numerically efficient moving window LoOP algorithm for monitoring industrial processes with complex data distributions, time-varying property and multiple operating modes. Some designed rules are introduced and incorporated with a two-step adaption approach to ensure that the monitoring model can be updated at a high speed. To cope with the multimode features, a semi-supervised monitoring strategy is employed, and an update termination rule is developed to prevent the monitoring model contaminated by faults or disturbances. Since the method is based on local probabilistic models, the accuracy of model is higher and it will be much easier to detect faulty conditions.
Finally, by applying the Gaussian error function, the local outlier probability indicating the probability that a sample is an outlier can be calculated as: n .pffiffiffi o LoOP 1 x j ¼ max 0; erf PLOF 1 x j 2 nPLOF 1
For low computation burden and practical applications, it is fast and reasonable to only update the information of those samples whose neighbors have changed due to the insertion and discard of samples. The key problems addressed in this section are how to find the affected samples and how to update their information. 2.1. Offline initialization To make an initialization and calculate the LoOP value for each sample xj (j = 1, 2, …, L) with dimension D in the initial window W1, its k nearest neighbors are found as follows, with its neighborhood set in W1 can be recognized as knn1(xj). vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 D u uX p ≠ j and x j ; xp ∈ W 1 d x j ; xp ¼ t xjn − xpn
ð1Þ
n¼1
Assuming that samples in knn1(xj) are centered around xj, then we can define probabilistic set distance as: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ffi X ; x =k d x pdist 1 x j ¼ λ j p xp ∈ knn1 ðx j Þ
ð2Þ
where λ is a weighted factor usually taken as 2. For estimating the density around xj, the probabilistic local outlier factor (PLOF) is defined as follows with function E(.) used to compute the expectation of PLOF in the current window. PLOF 1 x j ¼ pdist 1 x j = Ex
p ∈ knn1 ðx j Þ
h i pdist 1 xp −1
ð3Þ
To achieve normalization, the aggregate value nPLOF1 which can be considered as a standard deviation of PLOF values is obtained: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nPLOF 1 ¼ λ E ðPLOF 1 Þ2
ð4Þ
ð5Þ
where erf(.) is the Gaussian error function applied to obtain a probabilistic value. 2.2. Online updating and process monitoring By applying the moving window strategy, a two-step adaption procedure is introduced to update the monitoring model. Some more details of the adaption procedure for a window size L are as follows. Step 1: discard The effect of eliminating the oldest sample xi from the previous window Wi on the mean and variance can be evaluated as follows. e ¼ ðLμ i −xi Þ=ðL −1Þ μ
ð6Þ
e ¼ μi − μ e Δμ
ð7Þ
e ðmÞ2 ¼ σ
2. Adaptive Process Monitoring Based on Moving Window Loop
821
1 2 e ðmÞ 2 −½xi ðmÞ−μ i ðmÞ2 ðL−1Þ ½σi ðmÞ − Δ μ L−2 ðm ¼ 1; 2…; DÞ
ð
e ¼ diagσ e ð1Þ; σ e ð2Þ; ⋯; σ e ðDÞ Σ
ð8Þ
ð9Þ
where diag(.) is the function used to calculate the diagonal matrix. Eq. (6) describes the updating of the variable mean while Eqs. (7)–(9) describe the updating of the variable variance. After moving all the information about xi from the current monitor0 ing model, a set Si− 1(i N 1) is constructed to store the samples, in which xi is one of their k nearest neighbors. n o 0 0 S i−1 ¼ S i−1 ∪ x j ; if i ≠ j; x j ∈ W i and xi ∈ knni x j
ð10Þ
where knni(xj) represents the neighborhood set of sample xj in window Wi. Obviously, if xj ∈ Si0− 1, due to the deletion of xi, the neighborhood set knni(xj) will change. Step 2: insertion When a new sample xi + L is judged normal and added into the data matrix, the updated mean vector and variance in Wi + 1 are computed as follows. e þ xiþL =L μ iþ1 ¼ ðL−1Þμ
ð11Þ
e Δμ iþ1 ¼ μ iþ1 − μ
ð12Þ
2
σiþ1 ðmÞ ¼
1 e ðmÞ 2 þ ðL − 1Þ Δμ iþ1 ðmÞ 2 ðL − 2Þ σ L−1 2 ðm ¼ 1; 2…; DÞ þ xiþL ðmÞ − μ iþ1 ðmÞ
Σiþ1 ¼ diag σiþ1 ð1Þ; σiþ1 ð2Þ; ⋯; σiþ1 ðDÞ
ð13Þ ð14Þ
Eq. (11) describes the updating of the mean vector while Eqs. (12)–(14) describe the updating of the variance. However, only for those with new sample xi + L among their k nearest
822
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
neighbors, their neighborhood set knni(xj) will be updated to knni+1(xj). Therefore, the set S0i− 1 can be augmented to: n o 0 0 0 S i−1 ¼ S i−1 ∪ x j ; if x j ∉ S i−1 and xiþL ∈ knniþ1 x j
ð15Þ
The constructed set S0i − 1 should consist of the samples where xi or xi+L must be one of their k nearest neighbors. Then the probabilistic set distance is calculated by: 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X > 0 <λ > d x j ; xp =k if xp ∈ S i−1 xp ∈ knniþ1 ðx j Þ pdist iþ1 x j ¼ ð16Þ > > :pdist x otherwise i j where p ≠ j and xj ∈ Wi + 1. Since the PLOF of a sample will change not only with the change of its neighbors but also with the change of its probabilistic set distance, two sets e S i−1 and Si − 1 are constructed as: n o e S i−1 ∪ x j ; S i−1 ¼ e n o S i−1 ∪ x j ; S i−1 ¼ e
0 if xp ∈ knniþ1 x j and xp ∈ S i−1 0 if x j ∉ e S i−1 and x j ∈ S i−1
ð17Þ ð18Þ
where only the samples in Si − 1 need to update their PLOF values. By incorporating the updated probabilistic set distance, the PLOF and local outlier probability for the new data window Wi + 1 are computed as: pdistiþ1 x j h i −1 if x j ∈ S i−1 PLOF iþ1 x j ¼ Ex ∈ knn ðx Þ pdistiþ1 xp p iþ1 j > > > : PLOF x otherwise i j 8 > > > <
ð19Þ
19 8 8 0 < = > PLOF iþ1 x j > < max 0; erf @pffiffiffi > A if x j ∈ S i−1 : ; 2 nPLOF LoOP iþ1 x j ¼ ð20Þ iþ1 > > > : LoOP x otherwise i j
where p ≠ j ; x j ∈ W iþ1 and nPLOF iþ1 ¼ λ
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i E PLOF iþ1 2 .
Since LoOP assumes no specific data distribution, it will not be proper to apply 0.95 or 0.99 as the confidence limit for judging an outlier as that employed in GMM or other probabilistic methods. For those algorithms without assuming Gaussian distribution, kernel density estimation (KDE) is an effective method widely used in estimating control limits [20]. However, it will be time consuming to run the KDE algorithm every time after model updating. Obviously, if the local outlier probability of a sample is zero, it must be a normal one and will have negligible effect on estimating the control limit in the current window. However, the monitoring model must be updated every time when a normal sample is inserted, because new normal samples are appropriate representations of current states of the monitored industrial process and they can always bring useful information to guarantee the accuracy of the monitoring model. Therefore, samples with LoOP value 0 will be used only to update the monitoring model but the control limit will not be reevaluated. To cope with the multi-mode problem, model updating is operated through a semi-supervised switch strategy. If a mode change is previously known to occur, an alternative approach will be enabled to fit the transient stage and the new mode by blindly accepting every new sample as a normal one. However, there should be a termination rule to make the monitoring scheme switch back to its former state for fault detection as soon as possible to prevent the model from adapting to faulty conditions. During this transition, if the local outlier probability of a new sample is 0, there must be enough data to construct an accurate
local model in the new cluster. Therefore, when the LoOP value of a new sample is 0, it is reasonable to say that the model is ready for monitoring the new operating mode and the period of blind updating can be terminated. 2.3. Methodology The flow diagram of the proposed moving window local outlier probability (MWLoOP) monitoring scheme is shown in Fig. 1, with the detailed approach as follows. The offline modeling steps are summarized below: (1) Collect L samples from the current operating condition to construct the initial window W1. (2) Based on the standardized samples, an offline model is built according to Eqs. (1)–(5). (3) Specify a confidence level (1 − α)% and apply KDE to estimate a control limit for the LoOP values in the initial window W1. By introducing the switch rule, the online semi-supervised monitoring steps are summarized below. (1) For a new sample xi+L, standardization is first done by using the mean μi and variance Σi of the window Wi, where i N 1. (2) Calculate the distance from xi + L to samples xi+ 1, xi + 2, …, and xi+L− 1 in the window Wi and compute LoOPi(xi+ L). (3) If a mode change is previously known to occur, the flag value that is initially 0 should be set to 1 by an operator and every new sample is accepted as a normal one. Until the LoOP value of a new sample is 0, the flag will be automatically set back to 0. Go to Step (5). (4) If LoOPi(xi+L) N limit − LoOPi, where limit − LoOPi is the control limit of the window Wi, xi+L is judged as an outlier, and then the circulation will go for the next new sample. Otherwise, it continues to the next step. (5) The two-step adaption strategy is adopted to update the current model according to Eqs. (6)–(20) by discarding the oldest sample xi and inserting the newest sample xi+ L. (6) If LoOPi(xi+ L) N 0, the control limit is recalculated by KDE. (7) If several consecutive samples are judged as outliers, an alarm is triggered.
Fig. 1. Flow diagram of process monitoring scheme based on MWLoOP.
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
3. Case Study 3.1. Numerical example To demonstrate the superiority of the proposed method in dealing with complex data distribution and time-varying behavior, a numerical example is employed, which is similar to that used by Lee et al. [21]. Consider three source variables as follows: h i T x1 ðkÞ ¼ 2 diag ð cosð0:08kÞÞ sinð0:006kÞ
ð21Þ
x2 ðkÞ ¼ sign½ sinð0:03kÞ þ 9 cosð0:01kÞ
ð22Þ
x3 ðkÞ ∼ U ð−1; 1Þ
ð23Þ
x4 ðkÞ ∼ Nð2; 0:1Þ
ð24Þ
where k is a sampling index and k = 1, 2, …, 2000. Totally 2000 samples are generated with the following system: 2
0:86 6 −0:55 6 y ¼ Ax þ e ¼ 6 6 0:17 4 −0:33 0:89
3 2 3 e1 0:79 0:67 0:81 2 3 x 1 6 7 0:65 0:46 0:51 7 76 x2 7 6 e2 7 6 7 6 7 0:32 −0:28 0:13 7 74 x 5 þ 6 e3 7 4 e4 5 0:12 0:27 0:16 5 3 x4 e5 −0:97 −0:74 0:82
ð25Þ
where e = [e1, e2, e3, e4, e5]T are zero-mean white noises with variance 0.02 and y = [y1, y2, y3, y4, y5]T are the monitored variables. The first 1–1000 samples are normal ones, while at the 1001st sample, a slow drift 0.001(k − 1000) is added to A(1, 2) and A(2, 2) to simulate the time-varying behavior that should be adapted by the monitoring methods. Then a step bias of x2 with magnitude 3 is introduced at the 1501st sample.
823
In the moving window strategy, choosing a proper window size is a difficult task for compromising computational efficiency and model accuracy. A smaller window size means a lower computation load while a larger window size means higher model accuracy. As for the number of nearest neighbors, a large value will diminish the difference between normal samples and outliers while a small value will lead to inaccurate expression of local density. To verify the effect of the two parameters on Type I error, two tests on the 1–1500 samples are conducted and the results are shown in Fig. 2. Fig. 2(a) demonstrates the variant tendency of k versus Type I error with a window size of 750 while Fig. 2(b) shows the variant tendency of window size versus Type I error with the number of neighbors of 30. By compromising the trade-off between computation speed and model accuracy, the window size is chosen to be 750 and the number of nearest neighbors is 30 through trial and error method. The moving window PCA (MWPCA) applied in this paper is the algorithm proposed in [22], and the number of principle components is determined by the cumulative percentage variance (≥ 85%). For all methods applied in this paper, the confidence of control limit is set to be 99%. A fault occurs due to a step change in the non-Gaussian source variable x2. As shown in Table 1 and Fig. 3, moving window PCA fails to raise an alarm as the fault occurs since it cannot figure out the changes in variables with complex distributions, while conventional LoOP can show an obvious difference between normal and faulty conditions but it makes too many consecutive false alarms from the 1300th sample because it cannot handle the time-varying behaviors. In contrast, the proposed moving window LoOP algorithm has an acceptable Type I error and the best Type II error. The fault is rapidly detected without missing alarms.
3.2. Non-isothermal continuous stirred tank reactor The proposed method is compared to moving window PCA and conventional LoOP by simulating a non-isothermal continuous stirred tank reactor (CSTR). The process is shown in Fig. 4. It is a first order reaction,
Type I error/%
15
10
5
0
5
10
15
20
25
30
35
40
k
(a)
Type I error/%
80
60
40
20
0 100
200
300
400
500
600
700
window size
(b) Fig. 2. The effect of k and window size on Type I error.
800
900
1000
824
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
Table 1 Monitoring results of the numerical example Statistics of different methods/% SPE (MWPCA)
LoOP (MWLoOP)
14.13 0.00
0.00 43.40
0.20 74.20
2.47 0.00
20
4
15
3
SPE
T2
T2 (MWPCA)
10
0
2
1
5
0 0
500
1000
1500
2000
1000
1500
2000
(a) Moving window PCA-T2
(b) Moving window PCA-SPE 1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2 0 0
500
sample
1.0
0
0
sample
LoOP
LoOP
Type I error Type II error
LoOP (LoOP)
reactant A premixing with a solvent to product B. It should be noted that only the PI control loop for temperature T is active in this simulation. More details about the simulation condition should be referred to [23]. The nine monitored process variables are: T = outlet temperature, C = outlet concentration, FC = cooling water flow rate, T0 = inlet temperature, TC = cooling water temperature, CAA = concentration of pure A, CAS = concentration of solvent, FS = solvent flow, and FA = flow rate of constant A. Consider a very slow drift in reaction kinetic constant k0 to represent the time-varying feature of catalyst deactivation. With simulation time 5000 with 2500 samples are generated. The slow drift is introduced
500
1000
1500
2000
0
500
1000
1500
sample
sample
(c) LoOP
(d) Moving window LoOP
Fig. 3. Monitoring results of the numerical example.
Fig. 4. Diagram of CSTR process.
2000
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827 Table 2 Monitoring results for the three faults of CSTR Statistics of different methods/% LoOP (LoOP) T2 (MWPCA) SPE (MWPCA) LoOP (MWLoOP) Fault 1 Type I error 16.20 Type II error 0.20 Fault 2 Type I error 16.20 Type II error 0.40 Fault 3 Type I error 16.20 Type II error 0.20
0.85 2.60 0.85 57.40 0.85 8.60
4.90 0.60 4.75 9.60 4.75 16.60
4.70 0.20 4.65 4.40 4.65 6.20
from t = 2000 min as k0 = k0initial(1 − (t − 2000) × 10−4), while three kinds of fault are introduced from t = 4000 min: Fault 1: a step bias of cooling water temperature sensor with a magnitude of 1.5 K; Fault 2: a random noise of cooling water temperature that obeys uniform distribution U(−4, 4); Fault 3: a drift in the sensor of CAA and its slope is dCAA/ dt =0.001 kmol·m−3·min−1. According to the empirical guidance described in Section 3.1 and through trial and error method, the window size is chosen to be 700 and other parameters are the same with those in the numerical sample in Section 3.1. Monitoring results for the three faults are shown in Table 2. Type I errors of conventional LoOP for three faults are as high as 16.20% because it lacks the capability of adapting to time-varying processes. Compared with moving window PCA, Type I error of MWLoOP is acceptable but it performs better in Type II error. The reason lies in that with limited data samples, the Gaussian distribution assumptions of PCA cannot be fully satisfied, while the proposed method, which is free of distribution, can achieve a more accurate model. As a result, it is more sensitive to faulty conditions.
825
Fault 2 is a random non-Gaussian noise added in Tc. Fig. 5 shows that the monitoring results of conventional LoOP start to crash after t = 3200 min because it cannot deal with the time-varying behavior in the CSTR process. The SPE statistic of moving window PCA shows an acceptable result in Type II error, but T2 does not work well, because PCA is designed for extracting the Gaussian information into its feature space, and the remaining information including non-Gaussian features and disturbances will be separated into SPE. For the proposed method, a more accurate model is built without assuming any specific distribution of data. It seems that a few normal samples are judged as outliers, but they are discontinuous, so no alarm is triggered. Table 2 shows that Type II error of moving window LoOP is the most satisfactory. Therefore, a conclusion can be drawn that the proposed method is the most effective one compared to moving window PCA and conventional LoOP. Next, the ability of the proposed monitoring scheme to deal with mode changes is tested through Scenario 1 described below. The whole operation period consists of three stages. In the first stage t = 1–2000 min, the process is operated under mode 1, with the outlet temperature T setting at 368.25 K. In the second stage t = 2000– 4000 min, the reaction kinetic constant k0 changes as k0 = k0initial(1 − (t − 2000) × 10−4). Then the set-point of T changes to 370 K from t = 4000 min and after a while the process reaches steady state and run under mode 2. In the last stage t = 4000–5000 min, a step bias of concentration of pure A with a magnitude of 2 kmol·m−3 is introduced at t = 4800 min to simulate a faulty condition. Fig. 6 shows the monitoring results of moving window PCA and moving window LoOP for scenario 1. The vertical lines represent the stable time of the process judged by the dissimilarity index [10] and the stable time of MWLoOP judged by the proposed update termination rule. The statistic of proposed method becomes stable much faster than T2 and SPE, which means that the risk of adapting to faulty conditions will be
100
70 60
80
50
SPE
T2
40 30
60 40
20
0 2000
2200
1500
2000
0 1000
2500
1500
2500
(b) Moving window PCA-SPE
1.0
1.0
0.8
0.8
LoOP
LoOP
(a) Moving window PCA-T2
0.6
0.6
0.4
0.4
0.2
0.2
1500
2000
sample
sample
0 1000
2400
20
10 0 1000
5
2000
2500
0 1000
1500
2000
sample
sample
(c) LoOP
(d) Moving window LoOP
Fig. 5. Monitoring results for Fault 2 of CSTR.
2500
826
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
60
T2
40 20 0
0
500
1000
1500
2000
2500
sample
(a) Moving window PCA-T2
t=4084min
SPE
40 20
0
0
500
1000
1500
2000
2500
sample
(b) Moving window PCA-SPE
t=4054min
LoOP
1.0
0.5
0
0
500
1000
1500
2000
2500
sample
(c) Moving window LoOP Fig. 6. Monitoring results of Scenario 1.
much lower. In terms of the 100 faulty samples, Type II error of the proposed method is 5%, while that of T2 is 32% and SPE cannot detect this fault. Thus, the local model built by the proposed method is more accurate than the global model built by moving window PCA. From this point of view, the proposed method is much more practical, since it needs fewer samples to build a local model than a global one. 4. Conclusions The proposed moving window LoOP methodology offers many peculiarities, among which its fault detection capability, adaptive online implementation, and utility for multimode processes without prior knowledge requirement are addressed in this paper. To cope with the time-varying behavior, a two-step adaption approach is introduced to update the monitoring model while some updating rules are designed to reduce the computation load. For handling multi-mode characteristics, a semi-supervised switch strategy is incorporated and an update termination rule is designed to prevent the monitoring model fouled by faulty conditions. Due to the superiority of local probabilistic models, the proposed method can achieve a more accurate model and the monitoring efficiency can be easily maintained. Through a numerical example and a non-isothermal CSTR process, the flexibility and effectiveness of the proposed method are validated compared to moving window PCA and conventional LoOP. Future work will be focused on how to isolate faulty variables and how to tackle strong nonlinear property. References [1] W.H. Li, H.H. Yue, S.V. Cervantes, S.J. Qin, Recursive PCA for adaptive process monitoring, J. Process Control 10 (2000) 471–486.
[2] C. Rosen, J.A. Lennox, Multivariate and multi-scale monitoring of wastewater treatment operation, Water Res. 35 (2001) 3402–3410. [3] Y.W. Zhang, S. Li, Y.D. Teng, Dynamic process monitoring using recursive kernel principle component analysis, Chem. Eng. Sci. 72 (2012) 78–86. [4] X.Q. Liu, U. Kruger, T. Littler, L. Xie, S.Q. Wang, Moving window kernel PCA for adaptive monitoring of nonlinear processes, Chemom. Intell. Lab. Syst. 96 (2009) 132–143. [5] F.L. Wang, S. Tan, J. Peng, Y.Q. Chang, Process monitoring based on mode identification for multi-mode process with transitions, Chemom. Intell. Lab. Syst. 110 (2012) 144–155. [6] H.H. Ma, Y. Hu, H.B. Shi, Fault detection and identification based on neighborhood standardized local outlier factor method, Ind. Eng. Chem. Res. 52 (2013) 2389–2402. [7] P. Teppola, S.P. Mujunen, P. Minkkinen, Adaptive fuzzy C-means clustering in process monitoring, Chemom. Intell. Lab. Syst. 45 (1999) 23–38. [8] J.L. Liu, Modeling a large-scale nonlinear system using adaptive Takagi-Sugeno fuzzy model on PCA subspace, Ind. Eng. Chem. Res. 46 (2007) 788–800. [9] M. Petković, M.R. Rapaić, Z.D. Jeličić, A. Pisano, On-line adaptive clustering for process monitoring and fault detection, Expert Syst. Appl. 39 (2012) 10226–10235. [10] H.D. Jin, Y.H. Lee, G. Lee, C.H. Han, Robust recursive principal component analysis modeling for adaptive monitoring, Ind. Eng. Chem. Res. 45 (2006) 696–703. [11] Y.H. Lee, H.D. Jin, C.H. Han, On-line process state classification for adaptive monitoring, Ind. Eng. Chem. Res. 45 (2006) 3095–3107. [12] S.W. Choi, E.B. Martin, A.J. Morris, I.B. Lee, Adaptive multivariate statistical process control for monitoring time-varying processes, Ind. Eng. Chem. Res. 45 (2006) 3108–3118. [13] Z.Q. Ge, Z.H. Song, Online monitoring of nonlinear multiple mode processes based on adaptive local model approach, Control. Eng. Pract. 16 (2008) 1427–1437. [14] X. Xie, H.B. Shi, Dynamic multimode process modeling and monitoring using adaptive Gaussian mixture models, Ind. Eng. Chem. Res. 51 (2012) 5497–5505. [15] J. Yu, A particle filter driven dynamic Gaussian mixture model approach for complex process monitoring and fault diagnosis, J. Process Control 22 (2012) 778–788. [16] M. Kermit, O. Tomic, Independent component analysis applied on gas sensor array measurement data, IEEE Sens. J. 3 (2003) 218–228. [17] Z.Q. Ge, Z.H. Song, Process monitoring based on independent component analysisprincipal component analysis (ICA–PCA) and similarity factors, Ind. Eng. Chem. Res. 46 (2007) 2054–2063. [18] C.H. Zhao, F.R. Gao, F.L. Wang, Nonlinear batch process monitoring using phasebased kernel independent component analysis–principal component analysis (KICA–PCA), Ind. Eng. Chem. Res. 48 (2009) 9163–9174.
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827 [19] H.P. Kriegel, P. Kröger, E. Schubert, A. Zimek, LoOP: Local Outlier Probabilities, Proceedings of the 18th ACM conference on Information and knowledge management, 2009, pp. 1649–1652. [20] J.S. Lee, B.Y. Kang, S.H. Kang, Integrating independent component analysis and local outlier factor for plant-wide process monitoring, J. Process Control 21 (2011) 1011–1021.
827
[21] J.M. Lee, S.J. Qin, I.B. Lee, Fault detection and diagnosis based on modified independent component analysis, AIChE J. 52 (2006) 3501–3514. [22] X. Wang, U. Kruger, G.W. Irwin, Process monitoring approach using fast moving window PCA, Ind. Eng. Chem. Res. 44 (2005) 5691–5702. [23] S. Yoon, J.F. MacGregor, Fault diagnosis with multivariate statistical models part I: using steady state fault signatures, J. Process Control 11 (2001) 387–400.
Contents lists available at ScienceDirect
Chinese Journal of Chemical Engineering journal homepage: www.elsevier.com/locate/CJCHE
Process Monitor
Adaptive Local Outlier Probability for Dynamic Process Monitoring☆ Yuxin Ma, Hongbo Shi ⁎, Mengling Wang Key Laboratory of Advanced Control and Optimization for Chemical Processes of Ministry of Education, East China University of Science and Technology, Shanghai 200237, China
a r t i c l e
i n f o
Article history: Received 24 December 2013 Received in revised form 28 January 2014 Accepted 7 February 2014 Available online 19 June 2014 Keywords: Time-varying Complex data distribution Local outlier probability Multi-mode Fault detection
a b s t r a c t Complex industrial processes often have multiple operating modes and present time-varying behavior. The data in one mode may follow specific Gaussian or non-Gaussian distributions. In this paper, a numerically efficient moving window local outlier probability algorithm is proposed. Its key feature is the capability to handle complex data distributions and incursive operating condition changes including slow dynamic variations and instant mode shifts. First, a two-step adaption approach is introduced and some designed updating rules are applied to keep the monitoring model up-to-date. Then, a semi-supervised monitoring strategy is developed with an updating switch rule to deal with mode changes. Based on local probability models, the algorithm has a superior ability in detecting faulty conditions and fast adapting to slow variations and new operating modes. Finally, the utility of the proposed method is demonstrated with a numerical example and a non-isothermal continuous stirred tank reactor. © 2014 Chemical Industry and Engineering Society of China, and Chemical Industry Press. All rights reserved.
1. Introduction In industrial processes, operating conditions are usually affected by some slow variations denoted as time-varying characteristics, caused by some dynamic behavior such as seasonal fluctuation, catalyst deactivation, equipment aging, sensor or process drifting, preventive maintenance and cleaning [1]. Generally, effects of the time-varying behavior on the mean and covariance of variables cannot be neglected, so there may be many false alarms if conventional multivariate statistical process monitoring (MSPM) methods are applied directly [2]. In order to maintain process efficiency for a long period of time, numerous adaptive methods have been developed. Recursive MSPM methods and methods based on the moving window strategy are two alternative widely used approaches [3,4]. Multimodality is another important feature of industrial processes due to changes of market demands, alternations of feedstock or variations of manufacturing strategy. The difference between the characteristics of nearby operating conditions is always significant, so intensive studies have been carried out with either multiple local models or a single global model [5,6]. While it is more practical
☆ Supported by the National Natural Science Foundation of China (61374140), Shanghai Postdoctoral Sustentation Fund (12R21412600), the Fundamental Research Funds for the Central Universities (WH1214039), and Shanghai Pujiang Program (12PJ1402200). ⁎ Corresponding author. E-mail address: [email protected] (H. Shi).
to accommodate the time-varying behavior and multimode features together. The developed methods can be divided into two categories. One is the adaptive clustering methods. Teppola et al. [7] applied adaptive fuzzy C-means algorithms on the score values of principle component analysis (PCA) to monitor a wastewater treatment plant. Liu [8] used an adaptive Takagi-Sugeno fuzzy model on PCA subspace to model a large scale nonlinear system containing many operating regions. Since PCA is used as a preprocessing tool, monitoring results of these two methods more or less depend on and be restricted by the capability of PCA. Petković et al. [9] designed an on-line adaptive clustering method utilizing a generalized information potential. Although previously unseen functioning modes can be included by introducing an adaptive expert system, the method suffers from a non negligible detection delay. The other category is adaptive statistical methods. Improved recursive algorithms based on recursive PCA or the signed digraph were proposed with some if-then rules to distinguish process condition changes from disturbances [10–12]. Ge and Song [13] introduced the just-in-timelearning strategy to the modeling procedure of local least squares support vector regression and the residuals between the real output and the predicted one was analyzed by a two-step information extraction strategy. Xie and Shi [14] and Yu [15] developed two different dynamic fashions of Gaussian mixture model (GMM) separately based on the moving window strategy and a particle filter resampling method. The problem of complex data distributions in time-varying and multimode processes has scarcely been addressed. Although the moving window strategy has been proven to be effective, it still encounters some limitations when incorporated with statistical methods such as
http://dx.doi.org/10.1016/j.cjche.2014.05.015 1004-9541/© 2014 Chemical Industry and Engineering Society of China, and Chemical Industry Press. All rights reserved.
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
PCA, partial least squares (PLS) or GMM. Since the variables of an industrial process may satisfy specified Gaussian or non-Gaussian distributions, and high order statistics are usually helpful to reveal more information from the data [16–18], adaptive monitoring algorithms should be developed, which can explore both Gaussianity and nonGaussianity of process data. Local outlier probability (LoOP) is an unsupervised data mining technique proposed for outlier detection [19]. It combines the idea of local, density-based outlier scoring with a probabilistic, statistically-oriented approach, and assigns the probability of being an outlier to all data records. Since a normalization procedure is included, LoOP is independent of any specific data distribution. Therefore, a combination of LoOP and moving window strategy should be potential to tackle these problems. The main contribution of this paper is to propose a numerically efficient moving window LoOP algorithm for monitoring industrial processes with complex data distributions, time-varying property and multiple operating modes. Some designed rules are introduced and incorporated with a two-step adaption approach to ensure that the monitoring model can be updated at a high speed. To cope with the multimode features, a semi-supervised monitoring strategy is employed, and an update termination rule is developed to prevent the monitoring model contaminated by faults or disturbances. Since the method is based on local probabilistic models, the accuracy of model is higher and it will be much easier to detect faulty conditions.
Finally, by applying the Gaussian error function, the local outlier probability indicating the probability that a sample is an outlier can be calculated as: n .pffiffiffi o LoOP 1 x j ¼ max 0; erf PLOF 1 x j 2 nPLOF 1
For low computation burden and practical applications, it is fast and reasonable to only update the information of those samples whose neighbors have changed due to the insertion and discard of samples. The key problems addressed in this section are how to find the affected samples and how to update their information. 2.1. Offline initialization To make an initialization and calculate the LoOP value for each sample xj (j = 1, 2, …, L) with dimension D in the initial window W1, its k nearest neighbors are found as follows, with its neighborhood set in W1 can be recognized as knn1(xj). vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 D u uX p ≠ j and x j ; xp ∈ W 1 d x j ; xp ¼ t xjn − xpn
ð1Þ
n¼1
Assuming that samples in knn1(xj) are centered around xj, then we can define probabilistic set distance as: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ffi X ; x =k d x pdist 1 x j ¼ λ j p xp ∈ knn1 ðx j Þ
ð2Þ
where λ is a weighted factor usually taken as 2. For estimating the density around xj, the probabilistic local outlier factor (PLOF) is defined as follows with function E(.) used to compute the expectation of PLOF in the current window. PLOF 1 x j ¼ pdist 1 x j = Ex
p ∈ knn1 ðx j Þ
h i pdist 1 xp −1
ð3Þ
To achieve normalization, the aggregate value nPLOF1 which can be considered as a standard deviation of PLOF values is obtained: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nPLOF 1 ¼ λ E ðPLOF 1 Þ2
ð4Þ
ð5Þ
where erf(.) is the Gaussian error function applied to obtain a probabilistic value. 2.2. Online updating and process monitoring By applying the moving window strategy, a two-step adaption procedure is introduced to update the monitoring model. Some more details of the adaption procedure for a window size L are as follows. Step 1: discard The effect of eliminating the oldest sample xi from the previous window Wi on the mean and variance can be evaluated as follows. e ¼ ðLμ i −xi Þ=ðL −1Þ μ
ð6Þ
e ¼ μi − μ e Δμ
ð7Þ
e ðmÞ2 ¼ σ
2. Adaptive Process Monitoring Based on Moving Window Loop
821
1 2 e ðmÞ 2 −½xi ðmÞ−μ i ðmÞ2 ðL−1Þ ½σi ðmÞ − Δ μ L−2 ðm ¼ 1; 2…; DÞ
ð
e ¼ diagσ e ð1Þ; σ e ð2Þ; ⋯; σ e ðDÞ Σ
ð8Þ
ð9Þ
where diag(.) is the function used to calculate the diagonal matrix. Eq. (6) describes the updating of the variable mean while Eqs. (7)–(9) describe the updating of the variable variance. After moving all the information about xi from the current monitor0 ing model, a set Si− 1(i N 1) is constructed to store the samples, in which xi is one of their k nearest neighbors. n o 0 0 S i−1 ¼ S i−1 ∪ x j ; if i ≠ j; x j ∈ W i and xi ∈ knni x j
ð10Þ
where knni(xj) represents the neighborhood set of sample xj in window Wi. Obviously, if xj ∈ Si0− 1, due to the deletion of xi, the neighborhood set knni(xj) will change. Step 2: insertion When a new sample xi + L is judged normal and added into the data matrix, the updated mean vector and variance in Wi + 1 are computed as follows. e þ xiþL =L μ iþ1 ¼ ðL−1Þμ
ð11Þ
e Δμ iþ1 ¼ μ iþ1 − μ
ð12Þ
2
σiþ1 ðmÞ ¼
1 e ðmÞ 2 þ ðL − 1Þ Δμ iþ1 ðmÞ 2 ðL − 2Þ σ L−1 2 ðm ¼ 1; 2…; DÞ þ xiþL ðmÞ − μ iþ1 ðmÞ
Σiþ1 ¼ diag σiþ1 ð1Þ; σiþ1 ð2Þ; ⋯; σiþ1 ðDÞ
ð13Þ ð14Þ
Eq. (11) describes the updating of the mean vector while Eqs. (12)–(14) describe the updating of the variance. However, only for those with new sample xi + L among their k nearest
822
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
neighbors, their neighborhood set knni(xj) will be updated to knni+1(xj). Therefore, the set S0i− 1 can be augmented to: n o 0 0 0 S i−1 ¼ S i−1 ∪ x j ; if x j ∉ S i−1 and xiþL ∈ knniþ1 x j
ð15Þ
The constructed set S0i − 1 should consist of the samples where xi or xi+L must be one of their k nearest neighbors. Then the probabilistic set distance is calculated by: 8 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X > 0 <λ > d x j ; xp =k if xp ∈ S i−1 xp ∈ knniþ1 ðx j Þ pdist iþ1 x j ¼ ð16Þ > > :pdist x otherwise i j where p ≠ j and xj ∈ Wi + 1. Since the PLOF of a sample will change not only with the change of its neighbors but also with the change of its probabilistic set distance, two sets e S i−1 and Si − 1 are constructed as: n o e S i−1 ∪ x j ; S i−1 ¼ e n o S i−1 ∪ x j ; S i−1 ¼ e
0 if xp ∈ knniþ1 x j and xp ∈ S i−1 0 if x j ∉ e S i−1 and x j ∈ S i−1
ð17Þ ð18Þ
where only the samples in Si − 1 need to update their PLOF values. By incorporating the updated probabilistic set distance, the PLOF and local outlier probability for the new data window Wi + 1 are computed as: pdistiþ1 x j h i −1 if x j ∈ S i−1 PLOF iþ1 x j ¼ Ex ∈ knn ðx Þ pdistiþ1 xp p iþ1 j > > > : PLOF x otherwise i j 8 > > > <
ð19Þ
19 8 8 0 < = > PLOF iþ1 x j > < max 0; erf @pffiffiffi > A if x j ∈ S i−1 : ; 2 nPLOF LoOP iþ1 x j ¼ ð20Þ iþ1 > > > : LoOP x otherwise i j
where p ≠ j ; x j ∈ W iþ1 and nPLOF iþ1 ¼ λ
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h i E PLOF iþ1 2 .
Since LoOP assumes no specific data distribution, it will not be proper to apply 0.95 or 0.99 as the confidence limit for judging an outlier as that employed in GMM or other probabilistic methods. For those algorithms without assuming Gaussian distribution, kernel density estimation (KDE) is an effective method widely used in estimating control limits [20]. However, it will be time consuming to run the KDE algorithm every time after model updating. Obviously, if the local outlier probability of a sample is zero, it must be a normal one and will have negligible effect on estimating the control limit in the current window. However, the monitoring model must be updated every time when a normal sample is inserted, because new normal samples are appropriate representations of current states of the monitored industrial process and they can always bring useful information to guarantee the accuracy of the monitoring model. Therefore, samples with LoOP value 0 will be used only to update the monitoring model but the control limit will not be reevaluated. To cope with the multi-mode problem, model updating is operated through a semi-supervised switch strategy. If a mode change is previously known to occur, an alternative approach will be enabled to fit the transient stage and the new mode by blindly accepting every new sample as a normal one. However, there should be a termination rule to make the monitoring scheme switch back to its former state for fault detection as soon as possible to prevent the model from adapting to faulty conditions. During this transition, if the local outlier probability of a new sample is 0, there must be enough data to construct an accurate
local model in the new cluster. Therefore, when the LoOP value of a new sample is 0, it is reasonable to say that the model is ready for monitoring the new operating mode and the period of blind updating can be terminated. 2.3. Methodology The flow diagram of the proposed moving window local outlier probability (MWLoOP) monitoring scheme is shown in Fig. 1, with the detailed approach as follows. The offline modeling steps are summarized below: (1) Collect L samples from the current operating condition to construct the initial window W1. (2) Based on the standardized samples, an offline model is built according to Eqs. (1)–(5). (3) Specify a confidence level (1 − α)% and apply KDE to estimate a control limit for the LoOP values in the initial window W1. By introducing the switch rule, the online semi-supervised monitoring steps are summarized below. (1) For a new sample xi+L, standardization is first done by using the mean μi and variance Σi of the window Wi, where i N 1. (2) Calculate the distance from xi + L to samples xi+ 1, xi + 2, …, and xi+L− 1 in the window Wi and compute LoOPi(xi+ L). (3) If a mode change is previously known to occur, the flag value that is initially 0 should be set to 1 by an operator and every new sample is accepted as a normal one. Until the LoOP value of a new sample is 0, the flag will be automatically set back to 0. Go to Step (5). (4) If LoOPi(xi+L) N limit − LoOPi, where limit − LoOPi is the control limit of the window Wi, xi+L is judged as an outlier, and then the circulation will go for the next new sample. Otherwise, it continues to the next step. (5) The two-step adaption strategy is adopted to update the current model according to Eqs. (6)–(20) by discarding the oldest sample xi and inserting the newest sample xi+ L. (6) If LoOPi(xi+ L) N 0, the control limit is recalculated by KDE. (7) If several consecutive samples are judged as outliers, an alarm is triggered.
Fig. 1. Flow diagram of process monitoring scheme based on MWLoOP.
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
3. Case Study 3.1. Numerical example To demonstrate the superiority of the proposed method in dealing with complex data distribution and time-varying behavior, a numerical example is employed, which is similar to that used by Lee et al. [21]. Consider three source variables as follows: h i T x1 ðkÞ ¼ 2 diag ð cosð0:08kÞÞ sinð0:006kÞ
ð21Þ
x2 ðkÞ ¼ sign½ sinð0:03kÞ þ 9 cosð0:01kÞ
ð22Þ
x3 ðkÞ ∼ U ð−1; 1Þ
ð23Þ
x4 ðkÞ ∼ Nð2; 0:1Þ
ð24Þ
where k is a sampling index and k = 1, 2, …, 2000. Totally 2000 samples are generated with the following system: 2
0:86 6 −0:55 6 y ¼ Ax þ e ¼ 6 6 0:17 4 −0:33 0:89
3 2 3 e1 0:79 0:67 0:81 2 3 x 1 6 7 0:65 0:46 0:51 7 76 x2 7 6 e2 7 6 7 6 7 0:32 −0:28 0:13 7 74 x 5 þ 6 e3 7 4 e4 5 0:12 0:27 0:16 5 3 x4 e5 −0:97 −0:74 0:82
ð25Þ
where e = [e1, e2, e3, e4, e5]T are zero-mean white noises with variance 0.02 and y = [y1, y2, y3, y4, y5]T are the monitored variables. The first 1–1000 samples are normal ones, while at the 1001st sample, a slow drift 0.001(k − 1000) is added to A(1, 2) and A(2, 2) to simulate the time-varying behavior that should be adapted by the monitoring methods. Then a step bias of x2 with magnitude 3 is introduced at the 1501st sample.
823
In the moving window strategy, choosing a proper window size is a difficult task for compromising computational efficiency and model accuracy. A smaller window size means a lower computation load while a larger window size means higher model accuracy. As for the number of nearest neighbors, a large value will diminish the difference between normal samples and outliers while a small value will lead to inaccurate expression of local density. To verify the effect of the two parameters on Type I error, two tests on the 1–1500 samples are conducted and the results are shown in Fig. 2. Fig. 2(a) demonstrates the variant tendency of k versus Type I error with a window size of 750 while Fig. 2(b) shows the variant tendency of window size versus Type I error with the number of neighbors of 30. By compromising the trade-off between computation speed and model accuracy, the window size is chosen to be 750 and the number of nearest neighbors is 30 through trial and error method. The moving window PCA (MWPCA) applied in this paper is the algorithm proposed in [22], and the number of principle components is determined by the cumulative percentage variance (≥ 85%). For all methods applied in this paper, the confidence of control limit is set to be 99%. A fault occurs due to a step change in the non-Gaussian source variable x2. As shown in Table 1 and Fig. 3, moving window PCA fails to raise an alarm as the fault occurs since it cannot figure out the changes in variables with complex distributions, while conventional LoOP can show an obvious difference between normal and faulty conditions but it makes too many consecutive false alarms from the 1300th sample because it cannot handle the time-varying behaviors. In contrast, the proposed moving window LoOP algorithm has an acceptable Type I error and the best Type II error. The fault is rapidly detected without missing alarms.
3.2. Non-isothermal continuous stirred tank reactor The proposed method is compared to moving window PCA and conventional LoOP by simulating a non-isothermal continuous stirred tank reactor (CSTR). The process is shown in Fig. 4. It is a first order reaction,
Type I error/%
15
10
5
0
5
10
15
20
25
30
35
40
k
(a)
Type I error/%
80
60
40
20
0 100
200
300
400
500
600
700
window size
(b) Fig. 2. The effect of k and window size on Type I error.
800
900
1000
824
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
Table 1 Monitoring results of the numerical example Statistics of different methods/% SPE (MWPCA)
LoOP (MWLoOP)
14.13 0.00
0.00 43.40
0.20 74.20
2.47 0.00
20
4
15
3
SPE
T2
T2 (MWPCA)
10
0
2
1
5
0 0
500
1000
1500
2000
1000
1500
2000
(a) Moving window PCA-T2
(b) Moving window PCA-SPE 1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2 0 0
500
sample
1.0
0
0
sample
LoOP
LoOP
Type I error Type II error
LoOP (LoOP)
reactant A premixing with a solvent to product B. It should be noted that only the PI control loop for temperature T is active in this simulation. More details about the simulation condition should be referred to [23]. The nine monitored process variables are: T = outlet temperature, C = outlet concentration, FC = cooling water flow rate, T0 = inlet temperature, TC = cooling water temperature, CAA = concentration of pure A, CAS = concentration of solvent, FS = solvent flow, and FA = flow rate of constant A. Consider a very slow drift in reaction kinetic constant k0 to represent the time-varying feature of catalyst deactivation. With simulation time 5000 with 2500 samples are generated. The slow drift is introduced
500
1000
1500
2000
0
500
1000
1500
sample
sample
(c) LoOP
(d) Moving window LoOP
Fig. 3. Monitoring results of the numerical example.
Fig. 4. Diagram of CSTR process.
2000
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827 Table 2 Monitoring results for the three faults of CSTR Statistics of different methods/% LoOP (LoOP) T2 (MWPCA) SPE (MWPCA) LoOP (MWLoOP) Fault 1 Type I error 16.20 Type II error 0.20 Fault 2 Type I error 16.20 Type II error 0.40 Fault 3 Type I error 16.20 Type II error 0.20
0.85 2.60 0.85 57.40 0.85 8.60
4.90 0.60 4.75 9.60 4.75 16.60
4.70 0.20 4.65 4.40 4.65 6.20
from t = 2000 min as k0 = k0initial(1 − (t − 2000) × 10−4), while three kinds of fault are introduced from t = 4000 min: Fault 1: a step bias of cooling water temperature sensor with a magnitude of 1.5 K; Fault 2: a random noise of cooling water temperature that obeys uniform distribution U(−4, 4); Fault 3: a drift in the sensor of CAA and its slope is dCAA/ dt =0.001 kmol·m−3·min−1. According to the empirical guidance described in Section 3.1 and through trial and error method, the window size is chosen to be 700 and other parameters are the same with those in the numerical sample in Section 3.1. Monitoring results for the three faults are shown in Table 2. Type I errors of conventional LoOP for three faults are as high as 16.20% because it lacks the capability of adapting to time-varying processes. Compared with moving window PCA, Type I error of MWLoOP is acceptable but it performs better in Type II error. The reason lies in that with limited data samples, the Gaussian distribution assumptions of PCA cannot be fully satisfied, while the proposed method, which is free of distribution, can achieve a more accurate model. As a result, it is more sensitive to faulty conditions.
825
Fault 2 is a random non-Gaussian noise added in Tc. Fig. 5 shows that the monitoring results of conventional LoOP start to crash after t = 3200 min because it cannot deal with the time-varying behavior in the CSTR process. The SPE statistic of moving window PCA shows an acceptable result in Type II error, but T2 does not work well, because PCA is designed for extracting the Gaussian information into its feature space, and the remaining information including non-Gaussian features and disturbances will be separated into SPE. For the proposed method, a more accurate model is built without assuming any specific distribution of data. It seems that a few normal samples are judged as outliers, but they are discontinuous, so no alarm is triggered. Table 2 shows that Type II error of moving window LoOP is the most satisfactory. Therefore, a conclusion can be drawn that the proposed method is the most effective one compared to moving window PCA and conventional LoOP. Next, the ability of the proposed monitoring scheme to deal with mode changes is tested through Scenario 1 described below. The whole operation period consists of three stages. In the first stage t = 1–2000 min, the process is operated under mode 1, with the outlet temperature T setting at 368.25 K. In the second stage t = 2000– 4000 min, the reaction kinetic constant k0 changes as k0 = k0initial(1 − (t − 2000) × 10−4). Then the set-point of T changes to 370 K from t = 4000 min and after a while the process reaches steady state and run under mode 2. In the last stage t = 4000–5000 min, a step bias of concentration of pure A with a magnitude of 2 kmol·m−3 is introduced at t = 4800 min to simulate a faulty condition. Fig. 6 shows the monitoring results of moving window PCA and moving window LoOP for scenario 1. The vertical lines represent the stable time of the process judged by the dissimilarity index [10] and the stable time of MWLoOP judged by the proposed update termination rule. The statistic of proposed method becomes stable much faster than T2 and SPE, which means that the risk of adapting to faulty conditions will be
100
70 60
80
50
SPE
T2
40 30
60 40
20
0 2000
2200
1500
2000
0 1000
2500
1500
2500
(b) Moving window PCA-SPE
1.0
1.0
0.8
0.8
LoOP
LoOP
(a) Moving window PCA-T2
0.6
0.6
0.4
0.4
0.2
0.2
1500
2000
sample
sample
0 1000
2400
20
10 0 1000
5
2000
2500
0 1000
1500
2000
sample
sample
(c) LoOP
(d) Moving window LoOP
Fig. 5. Monitoring results for Fault 2 of CSTR.
2500
826
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827
60
T2
40 20 0
0
500
1000
1500
2000
2500
sample
(a) Moving window PCA-T2
t=4084min
SPE
40 20
0
0
500
1000
1500
2000
2500
sample
(b) Moving window PCA-SPE
t=4054min
LoOP
1.0
0.5
0
0
500
1000
1500
2000
2500
sample
(c) Moving window LoOP Fig. 6. Monitoring results of Scenario 1.
much lower. In terms of the 100 faulty samples, Type II error of the proposed method is 5%, while that of T2 is 32% and SPE cannot detect this fault. Thus, the local model built by the proposed method is more accurate than the global model built by moving window PCA. From this point of view, the proposed method is much more practical, since it needs fewer samples to build a local model than a global one. 4. Conclusions The proposed moving window LoOP methodology offers many peculiarities, among which its fault detection capability, adaptive online implementation, and utility for multimode processes without prior knowledge requirement are addressed in this paper. To cope with the time-varying behavior, a two-step adaption approach is introduced to update the monitoring model while some updating rules are designed to reduce the computation load. For handling multi-mode characteristics, a semi-supervised switch strategy is incorporated and an update termination rule is designed to prevent the monitoring model fouled by faulty conditions. Due to the superiority of local probabilistic models, the proposed method can achieve a more accurate model and the monitoring efficiency can be easily maintained. Through a numerical example and a non-isothermal CSTR process, the flexibility and effectiveness of the proposed method are validated compared to moving window PCA and conventional LoOP. Future work will be focused on how to isolate faulty variables and how to tackle strong nonlinear property. References [1] W.H. Li, H.H. Yue, S.V. Cervantes, S.J. Qin, Recursive PCA for adaptive process monitoring, J. Process Control 10 (2000) 471–486.
[2] C. Rosen, J.A. Lennox, Multivariate and multi-scale monitoring of wastewater treatment operation, Water Res. 35 (2001) 3402–3410. [3] Y.W. Zhang, S. Li, Y.D. Teng, Dynamic process monitoring using recursive kernel principle component analysis, Chem. Eng. Sci. 72 (2012) 78–86. [4] X.Q. Liu, U. Kruger, T. Littler, L. Xie, S.Q. Wang, Moving window kernel PCA for adaptive monitoring of nonlinear processes, Chemom. Intell. Lab. Syst. 96 (2009) 132–143. [5] F.L. Wang, S. Tan, J. Peng, Y.Q. Chang, Process monitoring based on mode identification for multi-mode process with transitions, Chemom. Intell. Lab. Syst. 110 (2012) 144–155. [6] H.H. Ma, Y. Hu, H.B. Shi, Fault detection and identification based on neighborhood standardized local outlier factor method, Ind. Eng. Chem. Res. 52 (2013) 2389–2402. [7] P. Teppola, S.P. Mujunen, P. Minkkinen, Adaptive fuzzy C-means clustering in process monitoring, Chemom. Intell. Lab. Syst. 45 (1999) 23–38. [8] J.L. Liu, Modeling a large-scale nonlinear system using adaptive Takagi-Sugeno fuzzy model on PCA subspace, Ind. Eng. Chem. Res. 46 (2007) 788–800. [9] M. Petković, M.R. Rapaić, Z.D. Jeličić, A. Pisano, On-line adaptive clustering for process monitoring and fault detection, Expert Syst. Appl. 39 (2012) 10226–10235. [10] H.D. Jin, Y.H. Lee, G. Lee, C.H. Han, Robust recursive principal component analysis modeling for adaptive monitoring, Ind. Eng. Chem. Res. 45 (2006) 696–703. [11] Y.H. Lee, H.D. Jin, C.H. Han, On-line process state classification for adaptive monitoring, Ind. Eng. Chem. Res. 45 (2006) 3095–3107. [12] S.W. Choi, E.B. Martin, A.J. Morris, I.B. Lee, Adaptive multivariate statistical process control for monitoring time-varying processes, Ind. Eng. Chem. Res. 45 (2006) 3108–3118. [13] Z.Q. Ge, Z.H. Song, Online monitoring of nonlinear multiple mode processes based on adaptive local model approach, Control. Eng. Pract. 16 (2008) 1427–1437. [14] X. Xie, H.B. Shi, Dynamic multimode process modeling and monitoring using adaptive Gaussian mixture models, Ind. Eng. Chem. Res. 51 (2012) 5497–5505. [15] J. Yu, A particle filter driven dynamic Gaussian mixture model approach for complex process monitoring and fault diagnosis, J. Process Control 22 (2012) 778–788. [16] M. Kermit, O. Tomic, Independent component analysis applied on gas sensor array measurement data, IEEE Sens. J. 3 (2003) 218–228. [17] Z.Q. Ge, Z.H. Song, Process monitoring based on independent component analysisprincipal component analysis (ICA–PCA) and similarity factors, Ind. Eng. Chem. Res. 46 (2007) 2054–2063. [18] C.H. Zhao, F.R. Gao, F.L. Wang, Nonlinear batch process monitoring using phasebased kernel independent component analysis–principal component analysis (KICA–PCA), Ind. Eng. Chem. Res. 48 (2009) 9163–9174.
Y. Ma et al. / Chinese Journal of Chemical Engineering 22 (2014) 820–827 [19] H.P. Kriegel, P. Kröger, E. Schubert, A. Zimek, LoOP: Local Outlier Probabilities, Proceedings of the 18th ACM conference on Information and knowledge management, 2009, pp. 1649–1652. [20] J.S. Lee, B.Y. Kang, S.H. Kang, Integrating independent component analysis and local outlier factor for plant-wide process monitoring, J. Process Control 21 (2011) 1011–1021.
827
[21] J.M. Lee, S.J. Qin, I.B. Lee, Fault detection and diagnosis based on modified independent component analysis, AIChE J. 52 (2006) 3501–3514. [22] X. Wang, U. Kruger, G.W. Irwin, Process monitoring approach using fast moving window PCA, Ind. Eng. Chem. Res. 44 (2005) 5691–5702. [23] S. Yoon, J.F. MacGregor, Fault diagnosis with multivariate statistical models part I: using steady state fault signatures, J. Process Control 11 (2001) 387–400.