- Email: [email protected]
Modeling and monitoring for handling nonlinear dynamic processes
Information Sciences 235 (2013) 97–105
Contents lists available at SciVerse ScienceDirect
Information Sciences journal homepage: www.elsevier.com/locate/ins
Modeling and monitoring for handling nonlinear dynamic processes Yingwei Zhang a,⇑, Jiayu An a, Zhiming Li a, Hong Wang b a b
State Laboratory of Synthesis Automation of Process Industry, Northeastern University, Shenyang, Liaoning 110004, PR China Univ. Manchester, Control Syst. Ctr., Sch. Elect. & Elect. Engn., Manchester M60 1QD, Lancs, England, UK
a r t i c l e
i n f o
Article history: Available online 21 April 2012 Keywords: Subspace separation Kernel method Common subspace Process monitoring
a b s t r a c t In this paper, a new online monitoring approach is proposed for handling the dynamical multimode problem in the industrial processes. The contributions are as follows: (1) extracting method of the common characteristics from different modes is proposed; (2) nonlinear dynamic monitoring method is proposed; and (3) a new model analysis method is proposed. There are both similarity and dissimilarity in the underlying correlations of different modes. After two different subspaces are separated, models of the common and specific subspaces are built respectively. Then the common subspace and specific subspace are analyzed, where the monitoring process is carried out in each subspace. When the mode switches, the specific monitoring model is changed. The corresponding confidence regions are constructed according to their models respectively. The effectiveness of the proposed method has been demonstrated via simulated examples. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction In order to ensure safety of the equipment operation and quality of product, the monitoring of the process performance has become a key issue. Multivariate statistical process control (MSPC) has been intensively studied in the last few decades. In particular, principal component analysis (PCA) and partial least squares (PLS) have been widely applied in industrial processes [1,3,8,9,11,16,20,21,26], where some improved methods such as kernel principal component analysis (KPCA) has shown great success in process monitoring and fault diagnosis [4,6,7,17–19]. Recently, monitoring the batch processes is needed for various reasons such as safety, waste-stream reduction, consistency and quality improvement. In this context, multi-way principal component analysis (MPCA) has been developed to deal with such batch processes [2,9,14]. However, in many industrial processes, the same production line is often used to produce different products. Therefore, there are often different production modes in the same production line, which makes dynamical multimode batch processes more complicated [8,10]. However, the MSPM methods are not available for the dynamical multimode processes. These methods may cause false alarms even when the process is operating under another nominal steady-state mode. Recently, recursive or adaptive PCA and PLS methods have been proposed [5,15,22–24]. Although these methods can be applied to treat the online process changes, they still lack the ability of coping with processes subjected to multiple operating modes [12,13]. Alternatively, model library based methods have been introduced [25,26], where predefined models match their corresponding operating modes. However, the effect of this method is not satisfactory since the nonlinearity can mot be considered. In this work, an online monitoring method is proposed for handling the problem of dynamical multimode in batch processes. The nonlinear similarity and dissimilarity of different modes are analyzed. When the mode switches, the specific
⇑ Corresponding author. Tel.: +86 24 83684946; fax: +86 24 83681006. E-mail address: [email protected] (Y. Zhang). 0020-0255/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.ins.2012.04.023
98
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
monitoring model is changed accordingly. By analyzing both common part and specific part, the different operating modes are identified and the faults of dynamical multimode process are diagnosed. The rest of this paper is organized as follows. The statistical analysis approach and monitoring method are proposed in Section 2. The simulation results are given to show the effectiveness of the proposed method in Section 3. Finally, conclusions are summarized in Section 4. 2. Statistical analysis and monitoring method In the measurement data, suppose that there are M industrial production patterns in the same production line. Therefore m m m T the multiple datasets X m ¼ xm 2 ðN m JÞ denote different modes, where J denotes the number of vari1 ; x2 ; . . . ; xl ; . . . ; xNm ables and m = 1, 2, . . . , M. The proposed algorithm maps the data from the original space to the feature space as Xm ? U(Xm). m m T In this context, UðX m Þ ¼ U xm are centered nonlinear mapping of the input variables and are supposed 1 ; U x2 ; . . . ; U xn P m to satisfy Nj¼1 Uðxm Þ ¼ 0. j The cross-mode relationship should be considered when the common underlying correlations are extracted. In each measurement dataset, it is always possible to find out a subset of vectors, which are representative enough to the other samples and can therefore substitute all samples by their linear combinations. Actually, in KPCA method they are equivalent to the KPCA loadings. Here they are called basis vectors. The major underlying correlations in the original measurement space are also represented by them and the associated distribution variances. Since any sub-basis in each dataset space, pm,j(j = 1, 2, . . . , J; m = 1, 2, . . . , M) (where m denotes different modes) must lie in h i m m m the span of the input observations, there exists combination coefficients am j ¼ a1;j ; a2;j ; . . . ; an;j , such that
pm;j ¼
Nm X m am ¼ UðX m ÞT am n;j U xi j
ð1Þ
i¼1
where xm i is a sample of Xm(Nm J). Therefore, the sub-basis vector pm,j is the function of the original observations in each dataset. Here, the similarity of variable correlations over sets can be obtained through the introduction of a global and common basis vector pg. That is, these real sub-basis vectors should be able to be comprehensively described and even substituted by the global basis. To figure out the common bases, the global and common vector pg is defined and then can be extracted. After solving the optimization problem, the global and common vector pg can be obtained. The optimization solutions and constraints are used for the extracting process. First, the common basis, pg, is defined as an orthogonal group configuration, which is a J dimensional basis. To make the correlation between pg and all original measurement data sets as close as possible, the maximum of the polynomial: 2 2 2 e1 pTg p1 þ e2 pTg p2 þ þ em pTg pm (where ei is a constant scalar) needs to be obtained. Then the object function is given as follows: 2
max R ¼ max
2 em pTg pm
M X
! ð2Þ
m¼1
By substituting Eq. (1) is substituted into Eq. (2), it can be obtained that 2
max R ¼ max
em pTg UðX
M X
m T
m
Þ a
2
! ð3Þ
m¼1
To obtain the common basis, the certain constraints are given as follows:
( s:t:
pTg pg ¼ 1 ðam ÞT am ¼ 1
where m = 1, 2, . . . , M, and the combination coefficient vector am is set to unit length. Indeed, it can be seen that em pTg UðX m ÞT am actually build the covariance model between the sub-basis vector emU(Xm)Tam and global basis vector pg. Therefore, the objective function involves the covariance information which is much better than the pure correlation analysis. The initial objective function is defined as follows, which can be expressed as an extreme value problem.
Fðpg ; a; kÞ ¼
M X m¼1
em pTg UðX m ÞT am
2
M X kg pTg pg 1 km ððam ÞT am 1Þ m¼1
where kg and km are constant scalars. To obtain the common basis pg, the following conditions should be satisfied:
ð4Þ
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
@Fðpg ; am ; kg ; kM Þ ¼0 @pg
99
ð5Þ
@Fðpg ; am ; kg ; kM Þ ¼0 @ am m ¼ 1; 2; . . . ; M
ð6Þ ð7Þ
Eqs. (5) and (6) represent that the partial derivative of pg and am are calculated respectively. From Eqs. (5) to (6), we have M X pffiffiffiffiffiffi em UðX
m T
Þ
pffiffiffiffiffiffi
em UðX m Þ pg ¼ kg pg
ð8Þ
m¼1
pffiffiffiffiffiffi pffiffiffiffiffiffi Þ is the covariance model between em UðX m ÞT and em UðX m Þ. Moreover, the two matrixes PM pffiffiffiffiffiffi ffiffiffiffiffi ffi p T are expanded into vectors form, the polynomial m¼1 em UðX m Þ em UðX m Þ can be rewritten as follows:
where
pffiffiffiffiffiffi em UðX
m T
Þ
pffiffiffiffiffiffi em UðX
m
Nm M X X pffiffiffiffiffiffi m pffiffiffiffiffiffi m T C¼ em U xi em U xi m¼1
! ð9Þ
i¼1
where U xm is the ith sample from the dataset of the mode m, then Eq. (4) can be rewritten as: i Nm M X X pffiffiffiffiffiffi m pffiffiffiffiffiffi m T em U xi em U xi m¼1
!!
pg ¼ kg pg
ð10Þ
i¼1
From Eq. (10), it can be obtained that
pffiffiffiffiffi 1 pffiffiffiffiffi 1 hpffiffiffiffiffi T pffiffiffiffiffi T T iT pffiffiffiffiffiffi pffiffiffiffiffiffi e1 U x1 ; e1 U x2 ; . . . ; em UðxCNm Þ e1 U x11 ; e1 U x12 ; . . . ; em U xCNm pg ¼ kg pg
ð11Þ
From Eq. (11), we have H X pffiffiffiffiffiffi pffiffiffiffiffiffi em Uðxi Þ em Uðxi ÞT pg ¼ kg pg
ð12Þ
i¼1
where H is the sample number of new dataset and can be formulated that
kg pg ¼
H X pffiffiffiffiffiffi
pffiffiffiffiffiffi
em Uðxi Þ em Uðxi ÞT pg ¼
i¼1
where have
pffiffiffiffiffiffi em Uðxi Þ is one sample from the new dataset in mode m. From Eq. (12), it
H X pffiffiffiffiffiffi
em Uðxi Þ; pg
pffiffiffiffiffiffi
em Uðxi Þ
ð13Þ
i¼1
pffiffiffiffiffiffi
pffiffiffiffiffiffi em Uðxi Þ; pg denotes the inter product between em Uðxi Þ and pg. Multiplying both sides of Eq. (10) by U(xk), we
kg
pffiffiffiffiffiffi
em Uðxk Þ; pg ¼
H X
pffiffiffiffiffiffi
pffiffiffiffiffiffi pffiffiffiffiffiffi em Uðxi Þ; pg em Uðxk Þ; em Uðxi Þ
ð14Þ
i¼1
Because the global basic vector pg must lie in the span of the input observations, there exists linear combination coefficients bi(i = 1, . . . , H), such that
pg ¼
H X pffiffiffiffiffiffi bi em Uðxi Þ
ð15Þ
i¼1
According to Eqs. (14) and (15), it can be shown that
kg
H H H X X pffiffiffiffiffiffi
X pffiffiffiffiffiffi
pffiffiffiffiffiffi
pffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffi bi em Uðxk Þ; em Uðxi Þ ¼ bi em Uðxj Þ; em Uðxi Þ em Uðxk Þ; em Uðxi Þ i¼1
j¼1
ð16Þ
i¼1
kx x k2 The Gaussian kernel is used in the method, which is defined as K i;j ¼ exp i c j , where i and j stand for the ith row and the jth column of x, respectively, and c is the width constant of the Gaussian function. The kernel matrix K is defined as follows
½Kij ¼ K ij ¼ ðUðxi Þ; Uðxj ÞÞ ¼ kðxi ; xj Þ
ð17Þ
where i = 1, 2, . . . , H, j = 1, 2, . . . , H. Before the kernel matrix K is used, it should be centered as follows:
e ¼ K 1H K K1H þ 1H K1H K
ð18Þ
100
2
1 . where 1H ¼ H1 4 .. 1
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
.. .
3 1 .. 5 HH . Eq. (16) can be rewritten as: . 2R 1
eb ¼ K e 2b kg K
ð19Þ T
where b = [b1, . . . , bH] . To find the solution of Eq. (19), we only need to calculate the following equation
eb kg b ¼ K
ð20Þ
From Eq. (20), the eigenvalues of the kernel matrix kg,1 P kg,2 P P kg,H and the eigenvectors b1,b2, . . . bH can be obtained. After the dominant eigenvalues bk(k = 1, 2, . . . , p) are retained, the global basis vector pg is computed according to Eq. (15) as follows:
pg;k ¼
H X pffiffiffiffiffiffi bi;k em Uðxi Þ
ð21Þ
i¼1
Once the p global basis vectors are obtained, the basis of the final common global subspace is given by Pg(J p) (where Pg ðJ pÞ ¼ pg;1 pg;2 pg;p ). Thus, a common subspace can be obtained within each mode. The control limits for the SPE are based on Box’s equation and are obtained by fitting a weighted v2-distribution to the reference distribution generated from normal operating data. By calculating the Hotelling-T2 of the common features, the faults in common information can be detected. If there are no faults in the common features, the calculation of the specific Hotelling-T2 can reveal the faults in the specific features. Further, the calculation of SPE can help to detect the process faults. 3. Simulation results The proposed method is used to the electro-fused magnesia furnace. The electro-fused magnesia furnace is one of the main equipment used to produce electro-fused magnesia, which belongs to a kind of hot electric arc furnace. The industrial process is shown in Fig. 1. In the melting process, furnace leakage is the serious fault, which is selected to be monitored in this paper. The furnace takes the light-burned magnesia as the raw material. It makes use of the heat generated by both the burden resistance when the current through the burden and the arc between the electrodes and the burden to melt the raw material, and then obtain the fused magnesia crystals with higher purity. In this example, there are two modes: mode A and mode B. The materials are powdery magnesium in mode A. The materials are massive magnesium in mode B. There are 300 process observations and three variables in each mode. First, a normal dataset is used to test the feasibility of the method. The test result is shown in Fig. 2a. The common T2 statistics are enclosed by the confidence region. Then the mode can be further identified by monitoring specific subspace as shown in Fig. 2. Specific T2 statistics of modes A and B are shown in Fig. 2b and c, respectively. Specific SPE statistics of modes A and B are shown in Fig. 2d and e, respectively. Obviously, due to the use of different specific parts, the out-of-control indications in mode B are much obvious. Combining their monitoring results in both subspaces, the affiliation of the two modes is definitely fixed and the operation status is also checked in both subspaces. Finally, the residual information is supervised and shown in Fig. 2e and f, which also indicates that the current batch is operating normally.
A
Operating Board A
A
1
2
4
3 5 7 8
6
Fig. 1. Diagram of electro-fused magnesium furnace: 1 – transformer, 2 – short circuit network, 3 – electrode holder, 4 – electrode, 5 – furnace shell, 6 – trolley, 7 – electric arc, and 8 – burden.
101
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
(a) common T2
6 4 2 0
0
10
20
30
40
50
60
70
Samples
specific T2
(b)
0.995 0.99 0.985 0.98
0
10
20
30
40
50
60
70
50
60
70
50
60
70
50
60
70
Samples
specific T2
(c)
1 0.98 0.96 0.94
0
10
20
30
40
Samples
(d) SPE
0.7
0.6
0.5
0
10
20
30
40
Samples
SPE
(e)
0.6
0.55
0.5
0
10
20
30
40
Samples Fig. 2. Monitoring results for (a) common T2 statistics, (b) specific T2 statistics of mode A, (c) specific T2 statistics of mode B, (d) SPE statistics of mode A, and (e) SPE statistics of mode B.
At the same time, from above analysis, the electro-fused magnesia furnace is working in which mode can be judged. From Fig. 2, the T2 of the common and the specific part are enclosed by the confidence region in mode A, and SPE of mode A also shows there are no faults. While T2 of the specific go beyond the confidence region in mode B, and SPE shows the test data is normal. From the test results of mode A, the test data is a normal data. From the test results of mode B, the specific part of mode B is abnormal. Therefore it can be known that the test data is a normal data belongs to mode A. Then a new abnormal batch run is applied in the algorithm, which belongs to mode A. Process faults occurred from the 60th observation. As it is shown in Fig. 3a, in the common part there are some samples which are not enclosed by the confidence region, and some part of the faults is not very obvious. Further checking via monitoring specific subspace shows that there are faults which start from 60 approximately in Fig. 3b. Also, the residual information is supervised which reveals the
102
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
common T2
(a) 150 100 50
0
0
20
40
60
80
100
120
Samples
specific T2
(b) 150 100 50 0
SPE
(c)
0
20
40
0
20
40
60
80
100
120
60
80
100
120
Samples
1
0.9
0.8
Samples Fig. 3. Monitoring results in mode A (a) common T2 statistics, (b) specific T2 statistics, and (c) SPE statistics.
common T2
(a) 150 100 50 0
0
20
40
60
80
100
120
80
100
120
80
100
120
Samples
specific T2
(b) 150 100 50 0
0
20
40
60
Samples
SPE
(c)
1
0.9
0.8
0
20
40
60
Samples Fig. 4. Monitoring results in mode B (a) common T2 statistics, (b) specific T2 statistics, and (c) SPE statistics.
103
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
MPCA T2
(a)
80 60 40 20 0
0
5
10
15
20
25
15
20
25
Samples
MPCA SPE
(b)
60 40 20 0
0
5
10
Samples Fig. 5. The T2 and SPE monitoring results by MPCA method.
Common T2
(a) 60 40 20 0
0
5
10
15
20
25
15
20
25
15
20
25
Samples
Specific T2
(b) 20 15 10 5 0 0
5
10
Samples
SPE
(c)
8 6 4 2 0
5
10
Samples Fig. 6. Monitoring results by the kernel statistical analysis and monitoring for (a) common T2 statistics, (b) specific T2 statistics, and (c) SPE statistics.
faults in the batch in Fig. 3c. Generally, clear and trustable alarms are revealed especially by the T2 monitoring system in specific systematic subspace, which means that the abnormal behavior mainly disturbs the specific part of information. Applied similar methods can detect a new abnormal batch which belongs to mode B. As it is shown in Fig. 4, the results reveal good fault detection performance. Generally, clear and stable alarms are revealed especially by the T2 monitoring system in common part, which means that the abnormal behavior mainly disturbs the common part of information. This
104
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
method can be achieved by switching between different modes, where the results shows that the faults can be effectively detected using the proposed method. After identifying the current operating mode, the proposed method and MPCA method are applied to testing faults. There is a fault data in mode B, which contains 24 samples and four variables. The faults happened from the 10th samples to the end. As it is shown in Fig. 5, the T2 of MPCA method has one fault that cannot be detected, and the SPE of MPCA method did not detect all fault either. Therefore, the MPCA method did not detect fault accurately. As it is shown in Fig. 6, the common part T2 of the proposed method has some fault pots that were not beyond the control limits, but the specific part T2 has detect all the faults. It can therefore be deduced that the faults mainly occurred at the specific part, showing that the SPE of the proposed method also detect all the faults. Therefore the proposed method has a strong detection capability and accuracy. 4. Conclusions In this work, a new method is provided for analysis of multimode batches. Since the cross-mode correlations are considered, the dynamical multimodes can be separated correctly using the proposed method. The similarity and dissimilarity of different modes are first analyzed. By performing dataset decomposition and subspace separation, the common information can be obtained from the specific information in each mode. The common information is the similar variable correlations over modes, and the specific information is the difference in each mode. By performing dataset decomposition and subspace separation, the underlying variations of different modes can be analyzed comprehensively. The strengths of the proposed strategy lie in not only the effective monitoring but also the appealing analysis results and comprehension for multi-mode problem. Fault detection can be accurate and obvious. The proposed method has been applied to the monitoring of an electro-fused magnesia furnace, where its effectiveness has been clearly seen. The monitoring for between-mode transition and the irregular transition dynamics over batches is a challenging issue for the future work. It can first reveal the problem of irregular transition dynamics over batches. Instead of identifying the definite boundary between mode and transition, all process patterns are regarded as possible transitions to accommodate their irregular dynamics. By using between-mode modeling, the underlying mode information is decomposed meaningfully for between-mode transition analysis. Acknowledgements The work is supported by China’s National 973 Program (2009CB320602 and 2009CB320604) and the NSF (60974057 and 61020106003). References [1] A. AlGhazzawi, B. Lennox, Monitoring a complex refining process using multivariate statistics, Control Engineering Practice 16 (2008) 294–307. [2] J.D. Carroll, Generalization of canonical correlation analysis to three or more sets of variables, in: Proceeding of the 76th Convention of the American Psychological Association, vol. 3, 1968, pp. 227–228. [3] L.H. Chiang, F.L. Russell, R.D. Braatz, Fault Detection and Diagnosis in Industrial Systems, Springer, London, 2001. [4] L.H. Chen, S.Y. Chang, Adaptive learning algorithm for principal component analysis, IEEE Transactions on Neural Networks 6 (5) (1995) 1255–1263. [5] B.S. Dayal, J.F. MacGregor, Recursive exponentially weighted PLS and its applications to adaptive control and prediction, Journal of Process Control 7 (3) (1997) 169–179. [6] T. Dune, L.B. Leopold, On relative convergence properties of principal component analysis algorithms, IEEE Transactions on Neural Networks 9 (2) (1998) 319–329. [7] Y. Gu, Y. Liu, Y.A. Zhang, Selective KPCA algorithm based on high-order statistics for anomaly detection in hyperspectral imagery, IEEE Geoscience and Remote Sensing Letters 5 (1) (2008) 43–47. [8] Z. Guo, H. Yue, H. Wang, A modified PCA based on the minimum error entropy, in: Proceedings of American Control Conference, Boston, 2004, pp. 3800–3801. [9] N. He, S.Q. Wang, L. Xie, An improved adaptive multi-way principal component analysis for monitoring streptomycin fermentation process, Chinese Journal of Chemical Engineering 12 (2004) 96–101. [10] D.H. Hwang, C. Han, Real-time monitoring for a process with multiple operating modes, Control Engineering Practice 7 (1999) 891–902. [11] U. Kruger, G. Dimitriadis, Diagnosis of process faults in chemical systems using a local partial least squares approach, AIChE Journal 54 (2008) 2581– 2596. [12] W. Li, H.H. Yue, S. Valle-Cervantes, S.J. Qin, Recursive PCA for adaptive process monitoring, Journal of Process Control 10 (2000) 471–486. [13] N. Lu, Y. Yi, F. Wang, F. Gao, A stage-based monitoring method for batch processes with limited reference data, in: 7th International Symposium on Dynamics and Control of Process Systems, (Dycops-7), Boston, USA, 2004. [14] P. Nomikos, Detection and diagnosis of abnormal batch operations based on multi-way principal component analysis, ISA Transactions 35 (1996) 259– 266. [15] S.J. Qin, Recursive PLS algorithms for adaptive data monitoring, Computers and Chemical Engineering 22 (1998) 503–514. [16] S.J. Qin, Statistical process monitoring: basics and beyond, Journal of Chemometrics 17 (2003) 480–502. [17] J.C. Sun, X.H. Li, Y. Yang, J.G. Luo, Y.H. Bai, Scaling the kernel function based on the separating boundary in input space: a data dependent way for improving performance of kernel methods, Information Sciences 184 (2012) 140–154. [18] Q. Wu, R. Law, Complex system fault diagnosis based on a fuzzy robust wavelet support vector classifier and an adaptive Gaussian particle swarm optimization, Information Sciences 180 (2010) 4514–4528. [19] L. Zhang, Q.X. Cao, A novel ant-based clustering algorithm using the kerne method, Information Sciences 181 (2011) 4658–4672. [20] Y.W. Zhang, S.J. Qin, Improved nonlinear fault detection technique and statistical analysis, AIChE Journal 54 (2008) 3207–3220. [21] Y.W. Zhang, H. Zhou, S.J. Qin, T.Y. Chai, Decentralized fault diagnosis of large-scale processes using multiblock kernel partial least squares, IEEE Transactions on Industrial Informatics 6 (2010) 3–12. [22] Y.W. Zhang, S.J. Qin, Nonlinear fault detection technique and statistical analysis, AIChE Journal 54 (12) (2008) 3207–3220. [23] Y.W. Zhang, Y. Zhang, Fault detection of non-Gaussian processes based on modified independent component analysis, Chemical Engineering Science 65 (16) (2010) 4630–4639.
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
105
[24] Y.W. Zhang, Y.D. Teng, Y. Zhang, Complex process quality prediction using modified kernel partial least squares, Chemical Engineering Science 65 (5) (2010) 2153–2158. [25] S.J. Zhao, Y.M. Xu, J. Zhang, A multiple PCA model based technique for the monitoring of processes with multiple operating modes, in: A. BarbosaPovoa, H. Matos (Eds.), European Symposium on Computer-Aided Chemical Engineering, vol. 18, 2004, pp. 865–870. [26] S.J. Zhao, J. Zhang, Y.M. Xu, Monitoring of processes with multiple operation modes through multiple principle component analysis models, Industrial and Engineering Chemistry Research 43 (2004) 7025–7035.
Contents lists available at SciVerse ScienceDirect
Information Sciences journal homepage: www.elsevier.com/locate/ins
Modeling and monitoring for handling nonlinear dynamic processes Yingwei Zhang a,⇑, Jiayu An a, Zhiming Li a, Hong Wang b a b
State Laboratory of Synthesis Automation of Process Industry, Northeastern University, Shenyang, Liaoning 110004, PR China Univ. Manchester, Control Syst. Ctr., Sch. Elect. & Elect. Engn., Manchester M60 1QD, Lancs, England, UK
a r t i c l e
i n f o
Article history: Available online 21 April 2012 Keywords: Subspace separation Kernel method Common subspace Process monitoring
a b s t r a c t In this paper, a new online monitoring approach is proposed for handling the dynamical multimode problem in the industrial processes. The contributions are as follows: (1) extracting method of the common characteristics from different modes is proposed; (2) nonlinear dynamic monitoring method is proposed; and (3) a new model analysis method is proposed. There are both similarity and dissimilarity in the underlying correlations of different modes. After two different subspaces are separated, models of the common and specific subspaces are built respectively. Then the common subspace and specific subspace are analyzed, where the monitoring process is carried out in each subspace. When the mode switches, the specific monitoring model is changed. The corresponding confidence regions are constructed according to their models respectively. The effectiveness of the proposed method has been demonstrated via simulated examples. Ó 2012 Elsevier Inc. All rights reserved.
1. Introduction In order to ensure safety of the equipment operation and quality of product, the monitoring of the process performance has become a key issue. Multivariate statistical process control (MSPC) has been intensively studied in the last few decades. In particular, principal component analysis (PCA) and partial least squares (PLS) have been widely applied in industrial processes [1,3,8,9,11,16,20,21,26], where some improved methods such as kernel principal component analysis (KPCA) has shown great success in process monitoring and fault diagnosis [4,6,7,17–19]. Recently, monitoring the batch processes is needed for various reasons such as safety, waste-stream reduction, consistency and quality improvement. In this context, multi-way principal component analysis (MPCA) has been developed to deal with such batch processes [2,9,14]. However, in many industrial processes, the same production line is often used to produce different products. Therefore, there are often different production modes in the same production line, which makes dynamical multimode batch processes more complicated [8,10]. However, the MSPM methods are not available for the dynamical multimode processes. These methods may cause false alarms even when the process is operating under another nominal steady-state mode. Recently, recursive or adaptive PCA and PLS methods have been proposed [5,15,22–24]. Although these methods can be applied to treat the online process changes, they still lack the ability of coping with processes subjected to multiple operating modes [12,13]. Alternatively, model library based methods have been introduced [25,26], where predefined models match their corresponding operating modes. However, the effect of this method is not satisfactory since the nonlinearity can mot be considered. In this work, an online monitoring method is proposed for handling the problem of dynamical multimode in batch processes. The nonlinear similarity and dissimilarity of different modes are analyzed. When the mode switches, the specific
⇑ Corresponding author. Tel.: +86 24 83684946; fax: +86 24 83681006. E-mail address: [email protected] (Y. Zhang). 0020-0255/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.ins.2012.04.023
98
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
monitoring model is changed accordingly. By analyzing both common part and specific part, the different operating modes are identified and the faults of dynamical multimode process are diagnosed. The rest of this paper is organized as follows. The statistical analysis approach and monitoring method are proposed in Section 2. The simulation results are given to show the effectiveness of the proposed method in Section 3. Finally, conclusions are summarized in Section 4. 2. Statistical analysis and monitoring method In the measurement data, suppose that there are M industrial production patterns in the same production line. Therefore m m m T the multiple datasets X m ¼ xm 2 ðN m JÞ denote different modes, where J denotes the number of vari1 ; x2 ; . . . ; xl ; . . . ; xNm ables and m = 1, 2, . . . , M. The proposed algorithm maps the data from the original space to the feature space as Xm ? U(Xm). m m T In this context, UðX m Þ ¼ U xm are centered nonlinear mapping of the input variables and are supposed 1 ; U x2 ; . . . ; U xn P m to satisfy Nj¼1 Uðxm Þ ¼ 0. j The cross-mode relationship should be considered when the common underlying correlations are extracted. In each measurement dataset, it is always possible to find out a subset of vectors, which are representative enough to the other samples and can therefore substitute all samples by their linear combinations. Actually, in KPCA method they are equivalent to the KPCA loadings. Here they are called basis vectors. The major underlying correlations in the original measurement space are also represented by them and the associated distribution variances. Since any sub-basis in each dataset space, pm,j(j = 1, 2, . . . , J; m = 1, 2, . . . , M) (where m denotes different modes) must lie in h i m m m the span of the input observations, there exists combination coefficients am j ¼ a1;j ; a2;j ; . . . ; an;j , such that
pm;j ¼
Nm X m am ¼ UðX m ÞT am n;j U xi j
ð1Þ
i¼1
where xm i is a sample of Xm(Nm J). Therefore, the sub-basis vector pm,j is the function of the original observations in each dataset. Here, the similarity of variable correlations over sets can be obtained through the introduction of a global and common basis vector pg. That is, these real sub-basis vectors should be able to be comprehensively described and even substituted by the global basis. To figure out the common bases, the global and common vector pg is defined and then can be extracted. After solving the optimization problem, the global and common vector pg can be obtained. The optimization solutions and constraints are used for the extracting process. First, the common basis, pg, is defined as an orthogonal group configuration, which is a J dimensional basis. To make the correlation between pg and all original measurement data sets as close as possible, the maximum of the polynomial: 2 2 2 e1 pTg p1 þ e2 pTg p2 þ þ em pTg pm (where ei is a constant scalar) needs to be obtained. Then the object function is given as follows: 2
max R ¼ max
2 em pTg pm
M X
! ð2Þ
m¼1
By substituting Eq. (1) is substituted into Eq. (2), it can be obtained that 2
max R ¼ max
em pTg UðX
M X
m T
m
Þ a
2
! ð3Þ
m¼1
To obtain the common basis, the certain constraints are given as follows:
( s:t:
pTg pg ¼ 1 ðam ÞT am ¼ 1
where m = 1, 2, . . . , M, and the combination coefficient vector am is set to unit length. Indeed, it can be seen that em pTg UðX m ÞT am actually build the covariance model between the sub-basis vector emU(Xm)Tam and global basis vector pg. Therefore, the objective function involves the covariance information which is much better than the pure correlation analysis. The initial objective function is defined as follows, which can be expressed as an extreme value problem.
Fðpg ; a; kÞ ¼
M X m¼1
em pTg UðX m ÞT am
2
M X kg pTg pg 1 km ððam ÞT am 1Þ m¼1
where kg and km are constant scalars. To obtain the common basis pg, the following conditions should be satisfied:
ð4Þ
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
@Fðpg ; am ; kg ; kM Þ ¼0 @pg
99
ð5Þ
@Fðpg ; am ; kg ; kM Þ ¼0 @ am m ¼ 1; 2; . . . ; M
ð6Þ ð7Þ
Eqs. (5) and (6) represent that the partial derivative of pg and am are calculated respectively. From Eqs. (5) to (6), we have M X pffiffiffiffiffiffi em UðX
m T
Þ
pffiffiffiffiffiffi
em UðX m Þ pg ¼ kg pg
ð8Þ
m¼1
pffiffiffiffiffiffi pffiffiffiffiffiffi Þ is the covariance model between em UðX m ÞT and em UðX m Þ. Moreover, the two matrixes PM pffiffiffiffiffiffi ffiffiffiffiffi ffi p T are expanded into vectors form, the polynomial m¼1 em UðX m Þ em UðX m Þ can be rewritten as follows:
where
pffiffiffiffiffiffi em UðX
m T
Þ
pffiffiffiffiffiffi em UðX
m
Nm M X X pffiffiffiffiffiffi m pffiffiffiffiffiffi m T C¼ em U xi em U xi m¼1
! ð9Þ
i¼1
where U xm is the ith sample from the dataset of the mode m, then Eq. (4) can be rewritten as: i Nm M X X pffiffiffiffiffiffi m pffiffiffiffiffiffi m T em U xi em U xi m¼1
!!
pg ¼ kg pg
ð10Þ
i¼1
From Eq. (10), it can be obtained that
pffiffiffiffiffi 1 pffiffiffiffiffi 1 hpffiffiffiffiffi T pffiffiffiffiffi T T iT pffiffiffiffiffiffi pffiffiffiffiffiffi e1 U x1 ; e1 U x2 ; . . . ; em UðxCNm Þ e1 U x11 ; e1 U x12 ; . . . ; em U xCNm pg ¼ kg pg
ð11Þ
From Eq. (11), we have H X pffiffiffiffiffiffi pffiffiffiffiffiffi em Uðxi Þ em Uðxi ÞT pg ¼ kg pg
ð12Þ
i¼1
where H is the sample number of new dataset and can be formulated that
kg pg ¼
H X pffiffiffiffiffiffi
pffiffiffiffiffiffi
em Uðxi Þ em Uðxi ÞT pg ¼
i¼1
where have
pffiffiffiffiffiffi em Uðxi Þ is one sample from the new dataset in mode m. From Eq. (12), it
H X pffiffiffiffiffiffi
em Uðxi Þ; pg
pffiffiffiffiffiffi
em Uðxi Þ
ð13Þ
i¼1
pffiffiffiffiffiffi
pffiffiffiffiffiffi em Uðxi Þ; pg denotes the inter product between em Uðxi Þ and pg. Multiplying both sides of Eq. (10) by U(xk), we
kg
pffiffiffiffiffiffi
em Uðxk Þ; pg ¼
H X
pffiffiffiffiffiffi
pffiffiffiffiffiffi pffiffiffiffiffiffi em Uðxi Þ; pg em Uðxk Þ; em Uðxi Þ
ð14Þ
i¼1
Because the global basic vector pg must lie in the span of the input observations, there exists linear combination coefficients bi(i = 1, . . . , H), such that
pg ¼
H X pffiffiffiffiffiffi bi em Uðxi Þ
ð15Þ
i¼1
According to Eqs. (14) and (15), it can be shown that
kg
H H H X X pffiffiffiffiffiffi
X pffiffiffiffiffiffi
pffiffiffiffiffiffi
pffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffi bi em Uðxk Þ; em Uðxi Þ ¼ bi em Uðxj Þ; em Uðxi Þ em Uðxk Þ; em Uðxi Þ i¼1
j¼1
ð16Þ
i¼1
kx x k2 The Gaussian kernel is used in the method, which is defined as K i;j ¼ exp i c j , where i and j stand for the ith row and the jth column of x, respectively, and c is the width constant of the Gaussian function. The kernel matrix K is defined as follows
½Kij ¼ K ij ¼ ðUðxi Þ; Uðxj ÞÞ ¼ kðxi ; xj Þ
ð17Þ
where i = 1, 2, . . . , H, j = 1, 2, . . . , H. Before the kernel matrix K is used, it should be centered as follows:
e ¼ K 1H K K1H þ 1H K1H K
ð18Þ
100
2
1 . where 1H ¼ H1 4 .. 1
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
.. .
3 1 .. 5 HH . Eq. (16) can be rewritten as: . 2R 1
eb ¼ K e 2b kg K
ð19Þ T
where b = [b1, . . . , bH] . To find the solution of Eq. (19), we only need to calculate the following equation
eb kg b ¼ K
ð20Þ
From Eq. (20), the eigenvalues of the kernel matrix kg,1 P kg,2 P P kg,H and the eigenvectors b1,b2, . . . bH can be obtained. After the dominant eigenvalues bk(k = 1, 2, . . . , p) are retained, the global basis vector pg is computed according to Eq. (15) as follows:
pg;k ¼
H X pffiffiffiffiffiffi bi;k em Uðxi Þ
ð21Þ
i¼1
Once the p global basis vectors are obtained, the basis of the final common global subspace is given by Pg(J p) (where Pg ðJ pÞ ¼ pg;1 pg;2 pg;p ). Thus, a common subspace can be obtained within each mode. The control limits for the SPE are based on Box’s equation and are obtained by fitting a weighted v2-distribution to the reference distribution generated from normal operating data. By calculating the Hotelling-T2 of the common features, the faults in common information can be detected. If there are no faults in the common features, the calculation of the specific Hotelling-T2 can reveal the faults in the specific features. Further, the calculation of SPE can help to detect the process faults. 3. Simulation results The proposed method is used to the electro-fused magnesia furnace. The electro-fused magnesia furnace is one of the main equipment used to produce electro-fused magnesia, which belongs to a kind of hot electric arc furnace. The industrial process is shown in Fig. 1. In the melting process, furnace leakage is the serious fault, which is selected to be monitored in this paper. The furnace takes the light-burned magnesia as the raw material. It makes use of the heat generated by both the burden resistance when the current through the burden and the arc between the electrodes and the burden to melt the raw material, and then obtain the fused magnesia crystals with higher purity. In this example, there are two modes: mode A and mode B. The materials are powdery magnesium in mode A. The materials are massive magnesium in mode B. There are 300 process observations and three variables in each mode. First, a normal dataset is used to test the feasibility of the method. The test result is shown in Fig. 2a. The common T2 statistics are enclosed by the confidence region. Then the mode can be further identified by monitoring specific subspace as shown in Fig. 2. Specific T2 statistics of modes A and B are shown in Fig. 2b and c, respectively. Specific SPE statistics of modes A and B are shown in Fig. 2d and e, respectively. Obviously, due to the use of different specific parts, the out-of-control indications in mode B are much obvious. Combining their monitoring results in both subspaces, the affiliation of the two modes is definitely fixed and the operation status is also checked in both subspaces. Finally, the residual information is supervised and shown in Fig. 2e and f, which also indicates that the current batch is operating normally.
A
Operating Board A
A
1
2
4
3 5 7 8
6
Fig. 1. Diagram of electro-fused magnesium furnace: 1 – transformer, 2 – short circuit network, 3 – electrode holder, 4 – electrode, 5 – furnace shell, 6 – trolley, 7 – electric arc, and 8 – burden.
101
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
(a) common T2
6 4 2 0
0
10
20
30
40
50
60
70
Samples
specific T2
(b)
0.995 0.99 0.985 0.98
0
10
20
30
40
50
60
70
50
60
70
50
60
70
50
60
70
Samples
specific T2
(c)
1 0.98 0.96 0.94
0
10
20
30
40
Samples
(d) SPE
0.7
0.6
0.5
0
10
20
30
40
Samples
SPE
(e)
0.6
0.55
0.5
0
10
20
30
40
Samples Fig. 2. Monitoring results for (a) common T2 statistics, (b) specific T2 statistics of mode A, (c) specific T2 statistics of mode B, (d) SPE statistics of mode A, and (e) SPE statistics of mode B.
At the same time, from above analysis, the electro-fused magnesia furnace is working in which mode can be judged. From Fig. 2, the T2 of the common and the specific part are enclosed by the confidence region in mode A, and SPE of mode A also shows there are no faults. While T2 of the specific go beyond the confidence region in mode B, and SPE shows the test data is normal. From the test results of mode A, the test data is a normal data. From the test results of mode B, the specific part of mode B is abnormal. Therefore it can be known that the test data is a normal data belongs to mode A. Then a new abnormal batch run is applied in the algorithm, which belongs to mode A. Process faults occurred from the 60th observation. As it is shown in Fig. 3a, in the common part there are some samples which are not enclosed by the confidence region, and some part of the faults is not very obvious. Further checking via monitoring specific subspace shows that there are faults which start from 60 approximately in Fig. 3b. Also, the residual information is supervised which reveals the
102
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
common T2
(a) 150 100 50
0
0
20
40
60
80
100
120
Samples
specific T2
(b) 150 100 50 0
SPE
(c)
0
20
40
0
20
40
60
80
100
120
60
80
100
120
Samples
1
0.9
0.8
Samples Fig. 3. Monitoring results in mode A (a) common T2 statistics, (b) specific T2 statistics, and (c) SPE statistics.
common T2
(a) 150 100 50 0
0
20
40
60
80
100
120
80
100
120
80
100
120
Samples
specific T2
(b) 150 100 50 0
0
20
40
60
Samples
SPE
(c)
1
0.9
0.8
0
20
40
60
Samples Fig. 4. Monitoring results in mode B (a) common T2 statistics, (b) specific T2 statistics, and (c) SPE statistics.
103
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
MPCA T2
(a)
80 60 40 20 0
0
5
10
15
20
25
15
20
25
Samples
MPCA SPE
(b)
60 40 20 0
0
5
10
Samples Fig. 5. The T2 and SPE monitoring results by MPCA method.
Common T2
(a) 60 40 20 0
0
5
10
15
20
25
15
20
25
15
20
25
Samples
Specific T2
(b) 20 15 10 5 0 0
5
10
Samples
SPE
(c)
8 6 4 2 0
5
10
Samples Fig. 6. Monitoring results by the kernel statistical analysis and monitoring for (a) common T2 statistics, (b) specific T2 statistics, and (c) SPE statistics.
faults in the batch in Fig. 3c. Generally, clear and trustable alarms are revealed especially by the T2 monitoring system in specific systematic subspace, which means that the abnormal behavior mainly disturbs the specific part of information. Applied similar methods can detect a new abnormal batch which belongs to mode B. As it is shown in Fig. 4, the results reveal good fault detection performance. Generally, clear and stable alarms are revealed especially by the T2 monitoring system in common part, which means that the abnormal behavior mainly disturbs the common part of information. This
104
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
method can be achieved by switching between different modes, where the results shows that the faults can be effectively detected using the proposed method. After identifying the current operating mode, the proposed method and MPCA method are applied to testing faults. There is a fault data in mode B, which contains 24 samples and four variables. The faults happened from the 10th samples to the end. As it is shown in Fig. 5, the T2 of MPCA method has one fault that cannot be detected, and the SPE of MPCA method did not detect all fault either. Therefore, the MPCA method did not detect fault accurately. As it is shown in Fig. 6, the common part T2 of the proposed method has some fault pots that were not beyond the control limits, but the specific part T2 has detect all the faults. It can therefore be deduced that the faults mainly occurred at the specific part, showing that the SPE of the proposed method also detect all the faults. Therefore the proposed method has a strong detection capability and accuracy. 4. Conclusions In this work, a new method is provided for analysis of multimode batches. Since the cross-mode correlations are considered, the dynamical multimodes can be separated correctly using the proposed method. The similarity and dissimilarity of different modes are first analyzed. By performing dataset decomposition and subspace separation, the common information can be obtained from the specific information in each mode. The common information is the similar variable correlations over modes, and the specific information is the difference in each mode. By performing dataset decomposition and subspace separation, the underlying variations of different modes can be analyzed comprehensively. The strengths of the proposed strategy lie in not only the effective monitoring but also the appealing analysis results and comprehension for multi-mode problem. Fault detection can be accurate and obvious. The proposed method has been applied to the monitoring of an electro-fused magnesia furnace, where its effectiveness has been clearly seen. The monitoring for between-mode transition and the irregular transition dynamics over batches is a challenging issue for the future work. It can first reveal the problem of irregular transition dynamics over batches. Instead of identifying the definite boundary between mode and transition, all process patterns are regarded as possible transitions to accommodate their irregular dynamics. By using between-mode modeling, the underlying mode information is decomposed meaningfully for between-mode transition analysis. Acknowledgements The work is supported by China’s National 973 Program (2009CB320602 and 2009CB320604) and the NSF (60974057 and 61020106003). References [1] A. AlGhazzawi, B. Lennox, Monitoring a complex refining process using multivariate statistics, Control Engineering Practice 16 (2008) 294–307. [2] J.D. Carroll, Generalization of canonical correlation analysis to three or more sets of variables, in: Proceeding of the 76th Convention of the American Psychological Association, vol. 3, 1968, pp. 227–228. [3] L.H. Chiang, F.L. Russell, R.D. Braatz, Fault Detection and Diagnosis in Industrial Systems, Springer, London, 2001. [4] L.H. Chen, S.Y. Chang, Adaptive learning algorithm for principal component analysis, IEEE Transactions on Neural Networks 6 (5) (1995) 1255–1263. [5] B.S. Dayal, J.F. MacGregor, Recursive exponentially weighted PLS and its applications to adaptive control and prediction, Journal of Process Control 7 (3) (1997) 169–179. [6] T. Dune, L.B. Leopold, On relative convergence properties of principal component analysis algorithms, IEEE Transactions on Neural Networks 9 (2) (1998) 319–329. [7] Y. Gu, Y. Liu, Y.A. Zhang, Selective KPCA algorithm based on high-order statistics for anomaly detection in hyperspectral imagery, IEEE Geoscience and Remote Sensing Letters 5 (1) (2008) 43–47. [8] Z. Guo, H. Yue, H. Wang, A modified PCA based on the minimum error entropy, in: Proceedings of American Control Conference, Boston, 2004, pp. 3800–3801. [9] N. He, S.Q. Wang, L. Xie, An improved adaptive multi-way principal component analysis for monitoring streptomycin fermentation process, Chinese Journal of Chemical Engineering 12 (2004) 96–101. [10] D.H. Hwang, C. Han, Real-time monitoring for a process with multiple operating modes, Control Engineering Practice 7 (1999) 891–902. [11] U. Kruger, G. Dimitriadis, Diagnosis of process faults in chemical systems using a local partial least squares approach, AIChE Journal 54 (2008) 2581– 2596. [12] W. Li, H.H. Yue, S. Valle-Cervantes, S.J. Qin, Recursive PCA for adaptive process monitoring, Journal of Process Control 10 (2000) 471–486. [13] N. Lu, Y. Yi, F. Wang, F. Gao, A stage-based monitoring method for batch processes with limited reference data, in: 7th International Symposium on Dynamics and Control of Process Systems, (Dycops-7), Boston, USA, 2004. [14] P. Nomikos, Detection and diagnosis of abnormal batch operations based on multi-way principal component analysis, ISA Transactions 35 (1996) 259– 266. [15] S.J. Qin, Recursive PLS algorithms for adaptive data monitoring, Computers and Chemical Engineering 22 (1998) 503–514. [16] S.J. Qin, Statistical process monitoring: basics and beyond, Journal of Chemometrics 17 (2003) 480–502. [17] J.C. Sun, X.H. Li, Y. Yang, J.G. Luo, Y.H. Bai, Scaling the kernel function based on the separating boundary in input space: a data dependent way for improving performance of kernel methods, Information Sciences 184 (2012) 140–154. [18] Q. Wu, R. Law, Complex system fault diagnosis based on a fuzzy robust wavelet support vector classifier and an adaptive Gaussian particle swarm optimization, Information Sciences 180 (2010) 4514–4528. [19] L. Zhang, Q.X. Cao, A novel ant-based clustering algorithm using the kerne method, Information Sciences 181 (2011) 4658–4672. [20] Y.W. Zhang, S.J. Qin, Improved nonlinear fault detection technique and statistical analysis, AIChE Journal 54 (2008) 3207–3220. [21] Y.W. Zhang, H. Zhou, S.J. Qin, T.Y. Chai, Decentralized fault diagnosis of large-scale processes using multiblock kernel partial least squares, IEEE Transactions on Industrial Informatics 6 (2010) 3–12. [22] Y.W. Zhang, S.J. Qin, Nonlinear fault detection technique and statistical analysis, AIChE Journal 54 (12) (2008) 3207–3220. [23] Y.W. Zhang, Y. Zhang, Fault detection of non-Gaussian processes based on modified independent component analysis, Chemical Engineering Science 65 (16) (2010) 4630–4639.
Y. Zhang et al. / Information Sciences 235 (2013) 97–105
105
[24] Y.W. Zhang, Y.D. Teng, Y. Zhang, Complex process quality prediction using modified kernel partial least squares, Chemical Engineering Science 65 (5) (2010) 2153–2158. [25] S.J. Zhao, Y.M. Xu, J. Zhang, A multiple PCA model based technique for the monitoring of processes with multiple operating modes, in: A. BarbosaPovoa, H. Matos (Eds.), European Symposium on Computer-Aided Chemical Engineering, vol. 18, 2004, pp. 865–870. [26] S.J. Zhao, J. Zhang, Y.M. Xu, Monitoring of processes with multiple operation modes through multiple principle component analysis models, Industrial and Engineering Chemistry Research 43 (2004) 7025–7035.