Operating optimality assessment and cause identification for nonlinear industrial processes

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

Contents lists available at ScienceDirect

Chemical Engineering Research and Design journal homepage: www.elsevier.com/locate/cherd

Operating optimality assessment and cause identification for nonlinear industrial processes Yan Liu a,∗ , Fuli Wang a,b , Yuqing Chang a,b a

College of Information Science & Engineering, Northeastern University, Shenyang, Liaoning 110819, China State Key Laboratory of Synthetical Automation for Process Industries (Northeastern University), Shenyang, Liaoning 110819, China

b

a r t i c l e

i n f o

a b s t r a c t

Article history:

To pursue optimal comprehensive economic benefit, the assessment on the operational per-

Received 20 May 2016

formance of the process becomes more and more critical for industrial processes. In this

Received in revised form 30 October

study, a novel assessment method regarding the operating optimality as well as the cause

2016

identification strategy are proposed for nonlinear industrial processes based on nonlinear

Accepted 1 November 2016

optimality related variation information (NORVI). The proposed method is committed to

Available online 9 November 2016

extracting the NORVI from each performance grade and using it to establish the assessment model. Since the process variation information unrelated to operating optimality is aban-

Keywords:

doned in evaluation, the accuracy of the assessment results is improved. Additionally, owing

Operating optimality assessment

to only the process data used in developing the assessment model, the time-consuming data

Nonlinear optimality related

alignment work is avoided, which improves the application efficiency of the proposed algo-

variation information

rithm. Furthermore, based on the similarities between the NORVI of the test data and those

Cause identification

of the modeling data of each performance grade, the actual operational performance of the

Nonlinear industrial process

process can be evaluated in real time. For the nonoptimal performance grade, the idea of variable selection is used to develop the cause identification strategy. Finally, the efficiency of the proposed method is illustrated with a case of gold hydrometallurgical process. © 2016 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1.

Introduction

In order to pursue the higher comprehensive economic benefits, the optimization of process operation becomes a key issue in modern industrial productions and has attracted great attention from both academia and practice recently (Del Rio-Chanona et al., 2016; Edgar et al., 2001; Hwang et al., 2013; Rossi et al., 2015; Srour et al., 2008). However, due to the presence of process disturbances and uncertainties, the operational performance of the process might deviate from its optimal setting. Thus, the development of the methods for the assessment of the optimality of a chosen operating policy becomes imperative. Process monitoring pays attention to the operating condition of normal or fault (Hu et al., 2013; Mori and Yu, 2014; Qin, 2012; Zhao et al., 2010). Comparatively speaking, the operating optimality assessment of the process is to evaluate the optimality of current operational performance under normal operating condition or get a measure on how far the operational performance is from the optimum. If the operational performance of the process can be evaluated accurately in real time and the causes responsible for the nonoptimality are identified timely, the appropriate adjustment strategy can be implemented on the process, and then some unnecessary losses are avoided. Recently, some research results on the operating optimality assessment of industrial process have emerged out gradually (Liu et al., 2014a, 2016; Ye et al., 2009). Liu et al. (2014a) developed the assessment method based on Total Projection to Latent Structure (T-PLS) algorithm. Since T-PLS (Zhou et al., 2010) can effectively extract the variation information closely related to the comprehensive economic benefits from the process data, the process variation information of each performance grade was extracted by T-PLS and used for operating optimality assessment. For simplicity, it is named as T-PLS-based assessment method. However, due to the production cycle of the products, there usually exists a time interval between the process measurements and the outputs, which leads to a deviation or misalignment between them. Additionally, the inherent prerequisite of T-PLS-based assessment is that the

∗ Corresponding author at: College of Information Science & Engineering, Northeastern University, 3 Lane 11, Wenhua Road, Heping District, Shenyang, Liaoning 110819, China. Fax: +86 24 23890912. E-mail address: xiaoyan [email protected] (Y. Liu). http://dx.doi.org/10.1016/j.cherd.2016.11.001 0263-8762/© 2016 Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

473

Fig. 1 – The relations among the original measurements, NORVI and NOUVI. process measurements and outputs have been aligned before modeling. Therefore, when the production cycle is unfixed with the influence of external environment, it is difficult to obtain accurate data alignment results, and then the accuracy of the assessment model may be affected to some extent. Moreover, the data alignment work is time-consuming and reduces the application efficiency of the algorithm. In view of this, the optimality related variation information (ORVI)-based assessment method was developed (Liu et al., 2016) to avoid the data alignment. However, for some complicated cases with nonlinear data, both the T-PLS and ORVI-based assessment methods perform poorly due to their potential linear assumption. To handle the problem posed by nonlinear characteristics, the Kernel Total Projection to Latent Structure (KT-PLS)based assessment method (Liu et al., 2014b), a nonlinear version of T-PLS-based assessment, was developed and applied to the nonlinear industrial processes. Nevertheless, the data alignment work faced by T-PLS-based assessment is also encountered in KT-PLS-based assessment. In order to deal with the issue of the operating optimality assessment for nonlinear processes, the nonlinear optimality related variation information (NORVI)-based assessment method is proposed in this study, and it is also named as NORVI-based assessment method for brevity. The basic idea of the proposed method is to extract the NORVI unique to each performance grade only by the process measurements firstly, and then establish the assessment models and further evaluate the operational performance of the process based on the NORVI. During this procedure, the tedious data alignment work can be avoided, and the application efficiency of the algorithm can be improved. In addition, because the process variation information related to operational performance is highlighted in operating optimality assessment, the proposed method has higher accuracy than the existing methods. For the nonoptimal performance grade, the cause identification strategy based on the idea of variable selection is developed for finding the responsible causes. The remainder of this paper is organized as follows. Section 2 presents the extraction of NORVI and the establishment of the assessment model. Then the new operating optimality assessment strategy is introduced in Section 3. In Section 4, the cause identification approach is proposed for finding the causes responsible for the nonoptimal performance grade. In case study, the effectiveness of the proposed approach is demonstrated by a gold hydrometallurgical process. Finally, the paper ends with some conclusions and acknowledgments.

2.

Extraction of NORVI and establishment of the assessment model

In reality, although the process variation information may be different from one performance grade to another, there also exists some similar or common process variation information among them for all of them originating from the same production process. The specific process variation information contained in each performance grade changes with the performance grade and constitutes the important feature in distinguishing the operational performance of the process, thus they are considered as the NORVI. On the contrary, the common process variation information across all the performance grades does not affect the operational performance of the process and not take effect in identifying the performance grade, so they are referred to as the nonlinear optimality unrelated variation information (NOUVI). In operating optimality assessment, the NOUVI can be considered as the interference information. If the whole process variations are used and do not make a distinction between NOUVI and NORVI, wrong judgment could be done on the actual operational performance of the process when the NOUVI changes. Thus, the NOUVI should be abandoned and only the NORVI is used in evaluation to ensure the accuracy and reliability of the assessment result. That is, extracting the NORVI accurately from each performance grade is the foundation of the offline modeling and operating optimality assessment.

2.1.

Extraction of NORVI

From the viewpoint of space decomposition, the original measurement space can be divided into two sub-spaces, as shown in Fig. 1. One is the sub-space spanned by the process variation information related to the operating optimality, and the other is with the process variation information unrelated to the operating optimality. Unlike the KT-PLS algorithm which has the comprehensive economic index as the aid to extract the NORVI directly, it is difficult to do that for the proposed method without any auxiliary information. However, this issue can be turned into extracting the NOUVI shared by all the performance grades. Hereafter, subtract the NOUVI from each performance grade, and the residual part naturally becomes the NORVI. The NOUVI can be analyzed from the following two aspects. The first aspect is to extract the common variable correlation shared by all the performance grades and constitute a common variable correlation sub-space; the second aspect, on the basis of the common variable correlation sub-space, is to further analyze the common variation information amplitude. In other words, if both the variable correlation and amplitude of the process variation information are almost the same for all of the performance grades, it can be considered as the NOUVI shared by all the performance grades. To extract the common variable correlation is essentially to reveal the common bases among multiple grades, which are also considered as the representative of those sub-basis vectors closely related across all the performance grades. Assume that there are C performance grades in a production process, and the modeling dataset of performance grade c are denoted as a matrix Xc = [xc,1 , xc,2 , ..., xc,Nc ]T ∈ RNc ×J , c = 1, 2, ..., C, where Nc and J are the numbers of samples and process variables, respectively. To deal with the nonlinear variable correlations, the computation trick proposed by Lee et al. (2004) and Schölkopf et al. (1998) is borrowed. It is that a nonlinear mapping function ˚(·) is introduced to project the vectors from the input space to a high-dimensional feature space F, i.e., xc,n ∈ RJ×1 → ˚(xc,n ) ∈ Rh×1 , n = 1, 2, ..., Nc , where h is the dimensionality of the feature

474

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

space F and can be arbitrarily large or possibly infinite. According to Cover’s theorem, the nonlinear data structure in the input space is more likely to be linear after high-dimensional nonlinear mapping (Haykin, 1999). This process is equivalent to linearizing T the nonlinear variable correlations. Then matrix Xc is mapped to ˚(Xc ) = [˚(xc,1 ), ˚(xc,2 ), ..., ˚(xc,Nc )] ∈ RNc ×h correspondingly. In each performance grade, a sub-basis vector is actually a linear combination of the observations in feature space F. Thus, if denote p¯ c,j ∈ Rh×1 as the j th sub-basis vector of ˚(Xc ), j = 1, 2, ..., h, there exists a combination coefficient vector ˛¯ c,j = [˛¯ c,j,1 , ..., ˛¯ c,j,Nc ]T , such that,

p¯ c,j =

Nc 

T

˛¯ c,j,n ˚(xc,n ) = ˚(Xc ) ˛¯ c,j .

(1)

n=1

To analyze the common variable correlation across multiple performance grades, the correlations of sub-basis vectors should be measured in terms of the closeness with each other over datasets. However, be different from the correlation measurement between two vectors, it is quite complicated that all set-to-set correlations between the sub-basis vectors are simultaneously and directly evaluated. In view of this, the correlation assessment of sub-basis vectors over sets are achieved by introducing a common basis vector pg . pg can be regarded as the pseudo (C + 1) st sub-basis vector as close to all C sub-bases as possible and even can substitute the real sub-basis vector. To obtain the common basis vector, a two-step basis vector extraction strategy is designed with the detailed mathematical derivation procedure shown in part one of Appendix A, named as ‘Extraction of Common Variable Correlation’. In Zhao et al. (2011), the linear version of the common variable correlation extraction strategy is proposed for multiple datasets. Here, it is extended to the nonlinear process for the first time. In the first step, the common basis vectors are preparatorily computed from the original measurements and then they are further refined by enhancing their correlations in the second step. In these two steps, different optimization solutions and constraints are used as shown in Eqs. (A1) and (A16), and both of them come down to the simple analytic solutions of unconstrained optimization problems as shown in Eqs. (A11) and (A18). Then the common basis matrix Pg is obtained, which is also consider as the common variable correlation sub-space spanned by the common basis vectors. The common variable correlation sub-space provides a platform for further extraction of the common variation information amplitude. The analysis of the common variation information amplitude is further carried out in the common variable correlation subspace Pg . The detailed calculation process is shown in part two of Appendix A with the heading of ‘Extraction of Common Variation Information Amplitude’. It is actually to pick out several common basis vectors from Pg , along each of which the amplitudes of the process variation information of different performance grades are all approximately equal. The equality of the amplitudes is measured by a predefined ratio , and then the common variable correlation sub-space Pg is further divided into two sub-spaces, ˘



˘

i.e., Pg = [Pg , P g ]. Pg is driven by a part of the common basis vectors, in which, the variable correlation and variation information 

amplitude are deemed to be consistent over performance grades; P g is spanned by the rest part of the common basis vectors, where the variable correlation is also the same but the variation information amplitude is unique to each definite performance ˘

grade. Therefore, if the original measurements ˚(Xc ) is directly projected to the sub-space Pg , the NOUVI is obtained and denoted ˘

˘

˘

T

˘

˙ c ) = (Xc ) − (Xc ), is the NORVI special to performance grade c. as ˚(Xc ) = ˚(Xc )Pg Pg . Then the residual part, i.e., (X

2.2.

Establishment of the assessment model based on NORVI

˙ c ) and formulated In order to remove the influence of process noise, Principal Component Analysis (PCA) is conducted on ˚(X below:

T ˙ c ) = T˙ c P˙ c + E˙ c , ˚(X

(2)

where T˙ c , P˙ c and E˙ c are the loading matrix, score matrix and residual matrix specific to performance grade c, respectively. Then the assessment models with parameters P˙ 1 , P˙ 2 , ..., P˙ C for C performance grades are established. The detailed PCA decomposition ˙ c ) is shown in Appendix B. procedure on ˚(X

3.

NORVI-based operating optimality assessment

In operating optimality assessment, owing to the fact that a single sample is difficult to characterize the operational performance of the process and is susceptible to noise, a data window constructed by a number of consecutive samples, i.e., Xon,k = [xon,k−H+1 , xon,k−H+2 , ..., xon,k ]T , is used as the basic analysis unit. The data window is actually a matrix. xon,k means the sample collected at the k sampling time, and the number of the samples included in Xon,k , i.e., H, is named as the window width. The data window plays the role of smoothing filter. Furthermore, to evaluate the operational performance of the process, an assessment index based on the similarity between the NORVI of the test data those of and the modeling data of each performance grade is define. The procedure of the operating optimality assessment based on NORVI is summarized below.

475

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

˘

˘

T

˙ x¯ on,k ) = ˚(x¯ on,k )T (I − Pg Pg ). Then First of all, calculate the mean of Xon,k , i.e., x¯ on,k and extract the NORVI from ˚(x¯ on,k ) by ˚( ˙ x¯ on,k ) with respect to performance grade c as: calculate the score vector of ˚(

˘

T

˘

T



˘

˘

T

T T

(Xc ) Q c = (kcon,k − kon,k W c )Q c , c = 1, 2, ..., C,

˙ x¯ on,k )P˙ c = ˚( ˙ x¯ on,k )˚(X ˙ c ) Q c = (x¯ on,k )(I − Pg Pg ) I − Pg Pg tcon,k = ˚(

(3)

where kon,k = [k(x1 , x¯ on,k ), ..., k(xN , x¯ on,k )] and kcon,k = [k(xc,1 , x¯ on,k ), ..., k(xc,Nc , x¯ on,k )] are the kernel vectors of ˚(x¯ on,k ) with respect T

T

to ˚(X) and ˚(Xc ) respectively, k(xn , x¯ on,k ) = ˚(xn )˚(x¯ on,k ) , n = 1, 2, ..., N,k(xc,nc , x¯ on,k ) = ˚(xc,nc )˚(x¯ on,k ) , nc = 1, 2, ..., Nc ; Q c is made −1

˘ ˘

T

−1

⁄2 DD  ⁄2 B˜ T K˜ T Kc . ˜ B ˜ ˙ c )˚(X ˙ c )T ; W c = K up of the first Mc eigenvectors of ˚(X Moreover, the similarity between the NORVI of test data and that of the modeling data of performance grade c is calculated 2 by dcon,k =tcon,k − t¯ c  , where t¯ c is the mean vector of T˙ c . Because ˚(Xc ) has been mean-centered, t¯ c = 0, and dcon,k =tcon,k 2 . It is worth noting that kcon,k and kon,k should be mean-centered before calculating the score vector t con,k , i.e., c k¯ on,k = kcon,k − 1Nc kKc,c − kcon,k 1Nc + 1Nc kKc,c 1Nc ,

(4)

k¯ on,k = kon,k − 1Nc kKTc − kon,k 1N + 1Nc kKTc 1N,

(5)

where 1Nc k is the Nc dimension row vector with all the values of element as 1/Nc . c To facilitate the assessment process, the assessment index on,k is defined below:

c on,k =

⎧ C

 ⎪ ⎪ ⎨ (1/dc )/ 1/ dq , on,k on,k

if dcon,k = / 0

⎪ ⎪ ⎩ c

if dcon,k = 0

(6)

q=1

on,k

= 1 and

q on,k

= 0(q = 1, 2, ..., C, q = / c),

c where on,k satisfies the conditions that

C 

c c on,k = 1 and 0 ≤ on,k ≤ 1.

c=1

In industrial process, it usually includes multiple performance grades and performance grade conversions (Liu et al., 2014a, 2016). To make a strict distinction between them, a assessment threshold,˛(0 < ˛ < 1), is introduced and the assessment rule q c = max {on,k } ≥ ˛, it indicates that the NORVI similar to those in Liu et al. (2014a) and Liu et al. (2016) is utilized in this study. If on,k 1≤q≤C

of test data is in line with that of performance grade c, and the process is running on performance grade c. If the above situation q q q is not met but the condition c = arg max {on,k |on,k−l+1 < ... < on,k } is satisfied, it can be determined that the current process is 1≤q≤C

gradually converting from the previous performance grade to performance grade c. l is a positive integer and used to judge the continuous increasing trend of the assessment indices. When both of the above two cases are not satisfied, it can be attributed to the influence of noise or uncertainties, and the assessment result is unchanged.

4.

Cause identification strategy

Be similar to the fault diagnosis issue under the abnormal operating conditions (Du et al., 2013; Peng et al., 2013; Qin et al., 2001; Zhao and Gao, 2015), when the process is operating on the nonoptimal performance grade, it is necessary to identify the responsible cause for further production adjustment and improvement of the operational performance of the process. In this section, the variable contribution-based cause identification method is proposed for nonlinear processes. It essentially uses the idea of variable selection problem which is to extract a few relevant variables making significant contributions to the predefined criterion. In cause identification, the search for cause variable is to select the variables having the relatively large influence on the deviation of assessment index with respect to optimal performance grade. Since the introduction of kernel function, the variable contribution cannot be decomposed directly from the assessment index, and the gradient of kernel function is used for measuring the contribution of each variable in NORVI-based assessment. Then the variables with relatively large contributions are considered as the possible cause variables for the nonoptimal operational performance. ∗ is calculated from the similarity of the test Use ‘*’ represents optimal performance grade. Because the assessment index on,k ∗ data with respect to optimal performance grade, the variable contributions to on,k is equivalent to those to d∗on,k . Here, d∗on,k can be rewritten as follows: d∗on,k



=t∗on,k 2 =t∗on,k t ∗T =tr t∗T t∗ on,k on,k on,k =tr =tr



∗ (k¯ on,k − k¯ on,k W ∗ )Q ∗

T Q T∗ (W T∗ k¯ on,k k¯ on,k W ∗



T 

∗ (k¯ on,k − k¯ on,k W ∗ )Q ∗

T ∗ − W T∗ k¯ on,k k¯ on,k



∗T − k¯ on,k k¯ on,k W ∗

∗T ∗ + k¯ on,k k¯ on,k )Q ∗



.

(7)

476

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

Then the contribution of variable j to d∗on,k can be defined below:

Contrj

=|

∂d∗on,k ∂x¯ on,k,j



Q T∗

= |tr

· x¯ on,k,j |





W T∗

¯ ∂ kon, kT k¯ on,k





W ∗ − W T∗

∂x¯ on,k,j

∗ ¯ ∂ kon, kT k¯ on,k





−

∂x¯ on,k,j

¯ k ∗ T k¯ on,k ∂ kon,





W∗ +

∂x¯ on,k,j

∗ ¯ k ∗ T k¯ on,k ∂ kon,

∂x¯ on,k,j



 Q∗

(8) · x¯ on,k,j |,

where x¯ on,k,j = x¯ on,k,j − x¯ ∗ , x¯ on,k,j and x¯ ∗ are the j variable of x¯ on,k and x¯ ∗ , and x¯ ∗ is the mean vector of X∗ ; tr{·} is the trace of a matrix. Suppose that the radial basis function, i.e., k(x, y) = exp(−x − y2 /2), is used for making kernel matrix in this paper, the parameter  is specified a priori by the user. In order to obtain the precise expression of Contrj , some preparations need to be done: (i) The gradient of a kernel function with respect to x¯ on,k,j : ∂k(xn , x¯ on,k ) 1 = (xn,j − x¯ on,k,j )k(xn , x¯ on,k ). ∂x¯ on,k,j 

(9)

(ii) The gradient of the product of two kernel functions to x¯ on,k,j : ∂k(xn , x¯ on,k )k(xm , x¯ on,k ) 1 = (xn,j + xm,j − 2x¯ on,k,j )k(xn , x¯ on,k )k(xm , x¯ on,k ). ∂x¯ on,k,j  T

(10) ∗T

∗

T

∗

(iii) The gradient of the p th row and q th column element of k¯ on,k k¯ on,k , k¯ on,k k¯ on,k and k¯ on,k k¯ on,k with respect to x¯ on,k,j , respectively. T

a. The p th row of k¯ on,k can be expressed as: T 1 k¯ on,k |p =k(xp , x¯ on,k ) − ep − k(xm , x¯ on,k ) + E∗ , N N

(11)

m=1

where ep =

N∗ 

k(x∗,n , xp )/N∗ and E∗ =

n=1

N∗ N  

k(x∗,n , xm )/NN∗ are constants.

n=1 m=1

T Thus, the p th row and q th column element of k¯ on,k k¯ on,k is,

1 k(xm , x¯ on,k )[k(xq , x¯ on,k ) + k(xp , x¯ on,k )] N N

(k¯ Ton,k k¯ on,k )

pq

= k(xp , x¯ on,k )k(xq , x¯ on,k ) + (E∗ − ep )k(xq , x¯ on,k ) + (E∗ − eq )k(xp , x¯ on,k ) −

m=1

 1  1 k(xm , x¯ on,k ) + 2 k(xm , x¯ on,k )k(xm , x¯ on,k ) + E (ep + eq − 2E∗ ) N N N

N

+

N

(12)

m=1 m =1

m=1

with E = ep eq − E∗ ep − E∗ eq + E2∗ as a constant. T Then the gradient of (k¯ k¯ on,k ) with respect to x¯ on,k,j is calculated as follows: on,k

¯ ∂(kon, kT k¯ on,k )pq ∂x¯ on,k,j

=

pq

(E∗ − eq ) (E∗ − ep ) 1 (x + xq,j − 2x¯ on,k,j )k(xp , x¯ on,k )k(xq , x¯ on,k ) + (xp,j − x¯ on,k,j )k(xp , x¯ on,k ) + (xq,j − x¯ on,k,j )k(xq , x¯ on,k )  p,j  

1  1  (xm,j + xp,j − 2x¯ on,k,j )k(xm , x¯ on,k )k(xp , x¯ on,k ) − (xm,j + xq,j − 2x¯ on,k,j )k(xm , x¯ on,k )k(xq , x¯ on,k ) N N N

−

N

m=1

m=1

(ep + eq − 2E∗ )  1  (xm,j − x¯ on,k,j )k(xm , x¯ on,k ) + (xm,j + xm ,j − 2x¯ on,k,j )k(xm , x¯ on,k )k(xm , x¯ on,k ). N N2 N

N

+

N

(13)

m=1 m =1

m=1

∗T b. The p th row of k¯ on,k can be expressed as: ∗T 1  k(x∗,n , x¯ on,k ) + E∗,∗ k¯ on,k |p =k(x∗,p , x¯ on,k ) − e∗p − N∗ N∗

n=1

(14)

477

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

with e∗p =

N∗ 

k(x∗,n , x∗,p )/N∗ and E∗,∗ =

n=1

N∗ N∗   n=1

k(x∗,n , x∗,n )/N∗2 .

n =1

∗T ∗ Furthermore, the p th row and q th column element of k¯ on,k k¯ on,k is,

1  k(x∗ ,n , x¯ on,k )[k(x∗ ,q , x¯ on,k ) + k(x∗ ,p , x¯ on,k )] N∗ N∗

¯∗ (k¯ ∗T on,k kon,k )

pq

= k(x∗,p , x¯ on,k )k(x∗,q , x¯ on,k )+(E∗,∗ − e∗p )k(x∗,q , x¯ on,k )+(E∗,∗ − e∗q )k(x∗,p , x¯ on,k )−

n=1

 1 ∗ 1  (ep + e∗q − 2E∗,∗ ) k(x∗,n , x¯ on,k ) + 2 k(x∗,n , x¯ on,k )k(x∗,n , x¯ on,k ) + E , N∗ N∗ N∗

N∗

+

N∗

(15)

n=1 n =1

n=1

where E = e∗p e∗q − E∗,∗ e∗p − E∗,∗ e∗q + E2∗,∗ is a constant. ∗T

∗

Then the gradient of (k¯ on,k k¯ on,k )pq with respect to x¯ on,k,j is expressed as follows: ∗ ¯ ∂(kon, k ∗ T k¯ on,k )pq

∂x¯ on,k,j +

−

+

(E∗,∗ − e∗p ) 

(E∗,∗ − e∗q ) 1 (x∗,p,j + x∗,q,j − 2x¯ on,k,j )k(x∗,p , x¯ on,k )k(x∗,q , x¯ on,k ) + (x∗,p,j − x¯ on,k,j )k(x∗,p , x¯ on,k )   1  (x∗,n,j + x∗,p,j − 2x¯ on,k,j )k(x∗,n , x¯ on,k )k(x∗,p , x¯ on,k ) N∗ N∗

(x∗,q,j − x¯ on,k,j )k(x∗,q , x¯ on,k ) −

n=1

N∗ 1 

N∗

=

(e∗p + e∗q − 2E∗,∗ ) 

(16)

N∗

(x∗,n,j + x∗,q,j − 2x¯ on,k,j )k(x∗,n , x¯ on,k )k(x∗,q , x¯ on,k ) +

n=1 N∗ N∗

N∗

(x∗,n,j − x¯ on,k,j )k(x∗,n , x¯ on,k )

n=1

1  (x∗,n,j + x∗,n ,j − 2x¯ on,k,j )k(x∗,n , x¯ on,k )k(x∗,n , x¯ on,k ). N∗2 n=1 n =1

T

∗

c. Based on Eqs. (11) and (14), the p th row and q th column element of k¯ on,k k¯ on,k has the following expression: 1  k(x∗,n , x¯ on,k )k(xp , x¯ on,k ) N∗ N∗

(k¯ Ton,k k¯ ∗on,k )

pq

= k(xp , x¯ on,k )k(x∗,q , x¯ on,k ) + (E∗ − ep )k(x∗,q , x¯ on,k ) + (E∗,∗ − e∗q )k(xp , x¯ on,k ) −

n=1

−

  1 1 1 1  1 k(xm , x¯ on,k )k(x∗,q , x¯ on,k )+ (ep −E∗ ) k(x∗,n , x¯ on,k )+ (e∗q − E∗,∗ ) k(xm , x¯ on,k ) + k(xm , x¯ on,k )k(x∗,n , x¯ on,k )+E , N N∗ N N∗ N N

N∗

N

m=1

n=1

m=1

N

N∗

m=1 n=1

(17) 

where E = ep e∗q − E∗,∗ ep − E∗ e∗q + E∗,∗ E∗ is a constant. T

∗

Then the gradient of (k¯ on,k k¯ on,k )pq with respect to x¯ on,k,j is represented below: ∗ ¯ ∂(kon, kT k¯ on,k )pq

∂x¯ on,k,j +

−

+

(E∗,∗ − e∗q )  N 1 

N (e∗q

m=1

(E∗ − ep ) 1 (x∗,q,j − x¯ on,k,j )k(x∗,q , x¯ on,k ) = (xp,j + x∗,q,j − 2x¯ on,k,j )k(xp , x¯ on,k )k(x∗,q , x¯ on,k ) +  

n=1

(ep − E∗ )  (xm,j + x∗,q,j − 2x¯ on,k,j )k(xm , x¯ on,k )k(x∗,q , x¯ on,k ) + (x∗,n,j − x¯ on,k,j )k(x∗,n , x¯ on,k ) N∗ N∗

(18)

n=1

N − E∗,∗ ) 

N

1  (x∗,n,j + xp,j − 2x¯ on,k,j )k(x∗,n , x¯ on,k )k(xp , x¯ on,k ) N∗ N∗

(xp,j − x¯ on,k,j )k(xp , x¯ on,k ) −

m=1

(xm,j − x¯ on,k,j )k(xm , x¯ on,k ) +

N∗ N 1 

N∗ N

(x∗,n,j + xm,j − 2x¯ on,k,j )k(x∗,n , x¯ on,k )k(xm , x¯ on,k ).

n=1 m=1

Substitute Eqs. (13), (16) and (18) into Eq. (8), the variable contribution to d∗k is obtained. In cause identification, the variables with relatively large contributions are considered as the responsible cause variables. The flowchart of the proposed operating optimality assessment and cause identification method based on NORVI is shown in Fig. 2.

5.

Case study: gold hydrometallurgical process

In view of the high value of gold, to ensure the satisfactory operational performance is very important for the gold hydrometallurgical process. In this study, the process data collected from an actual gold hydrometallurgical process are used to verify the feasibility and effectiveness of the proposed assessment approach.

478

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

Fig. 2 – Schematic diagram of the proposed approach.

Fig. 3 – Schematic diagram of gold hydrometallurgical process.

5.1.

Process description

This gold hydrometallurgical process consists of three subprocesses: cyanide leaching, washing, and cementation. The schematic diagram of this process is shown in Fig. 3. First of all, the flotation pulp flows into the cyanide leaching subprocess, meanwhile append the sodium cyanide (NaCN) and oxygen to the leaching tank. NaCN is used as the reagents, and dissolved oxygen provides oxidation-reduction potential to promote the reaction. Afterwards, the pulp is feed into the washing subprocess to separate the high-concentration [Au(CN)2 ]− solution from the pulp. For the purpose of solid–liquid separating, the vertical filter presses is used in this subprocess. Then separated high-concentration [Au(CN)2 ]− solution is feed into the cementation subprocess to recover gold. Zinc powder is used as the deoxidizer, and the frame filter press is used as the reaction carrier. The process variables used for operating optimality assessment are listed in Table 1. Through in-depth research on the technology and mechanism of gold hydrometallurgical process, the process variables 1–7 listed in Table 1 are all the important variables reflecting the operational performance of the process. Additionally, there are strong nonlinear coupling relationships within three groups of process variables: the first group includes the pulp concentration, oxygen concentration and the air

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

479

Table 1 – Process variables used for operating optimality assessment. No.

Process variables

Units

1 2 3 4 5 6 7 8 9 10

Pulp concentration in cyanide leaching NaCN flow of leaching tank 1 CN− concentration in leaching tank 1 Air flow Oxygen concentration in the pulp Zinc-powder addition amount [Au(CN)2 ]− concentration Leaching tank 1 level Leaching tank 2 level Leaching tank 3 level

% mL/min mg/L m3 /h mg/L Kg/h mg/L m m m

Fig. 4 – ORVI-based assessment results for the first case. flow, the second one contains NaCN flow and CN− concentration, and the third group has zinc-powder addition amount and [Au(CN)2 ]− concentration. All of these process variables enhance the nonlinear characteristics of the whole simulation data. However, the process variables 8–10 could not affect the optimality of the operational performance of the process and only used for simulating the interferences of the NOUVI, which will be further explained in the text below. Based on the process knowledge and expert experience, the process can be divided into three performance grades, i.e., general, suboptimal and optimal. Both general and suboptimal are the nonoptimal performance grades, and the reasons can be attributed to the lower setting of the NaCN flow (variable 2) than that of the optimal operating points. The modeling data of the three performance grades are collected from the actual process and constitute the modeling datasets X1 ∈ R2000×10 , X2 ∈ R2000×10 and X3 ∈ R2000×10 , where the subscripts 1, 2, 3 correspond to the general, suboptimal and optimal performance grade respectively for simplicity. To establish the assessment models for each performance grade, the relaxation factor ı and the parameter  are set as 0.1 and 100, respectively.

5.2.

Operating optimality assessment

To verify the efficiency of the proposed assessment method, two cases are considered. The first case is that leaching tank 1 level (variable 8) gradually increases within a safe range when the process is running on the optimal performance grade. In production, the main function of the leaching tank is to buffer the pulp, which can prevent the risk of groove accident. Known from the process knowledge and expertise, normal fluctuations of the leaching tank level in the safe range cannot affect the operational performance of the process, so it is considered as the change of NOUVI. In the second case, when the performance grade is optimal, the zinc-powder addition amount (variable 6) gradually decreases and deviates from the optimal set point due to human error. It is actually the change of NORVI and leads to nonoptimal performance grade. The test dataset for each case includes 1800 samples with the operating trend as general → conversion → suboptimal → conversion → optimal, and both of the two cases are introduced at the 1501th sampling, respectively. The parameters used for assessment are set as follows: H = 35, ˛ = 0.85 and l = 5. For comparison, the ORVI-based assessment method is also applied to the above two cases in this study. Figs. 4 and 5 are the assessment results for the first case based on ORVI and NORVI assessment methods, respectively. Fig. 4(a)–(c) reveals the assessment index with regard to three performance grades correspondingly, i.e., general, suboptimal, and optimal, and Fig. 5(a)–(c) have the same meanings. From Figs. 4 and 5, it can be seen that the erroneous evaluations occur in the ORVI-based assessment method, but the assessment results based on the NORVI seem right. As can be seen in Table 2, the accuracy rate of the assessment result is up to 95.72% for NORVI-based assessment, while that is only 89.28% for ORVI-based assessment. Here, the accuracy rate is calculated as the percentage between the numbers of correct evaluated samplings and the total test samplings. In addition,

480

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

Fig. 5 – NORVI-based assessment results for the first case. Table 2 – Accuracy rates of ORVI and NORVI-based assessment methods for the first case. Method

Numbers of correct evaluated samplings

Accuracy rates(%)

ORVI NORVI

1607 1723

89.28 95.72

Table 3 – Comparison between the actual situation and the assessment results obtained from NORVI-based assessment method. Performance grade

Actual situation (samplings)

NORVI-based assessment (samplings)

General Conversion Suboptimal Conversion Optimal

1–500 501–650 651–1150 1151–1300 1301–1800

1–520 521–660 661–1173 1174–1324 1325–1800

Fig. 6 – ORVI-based cause identification results for the first case at the 1st sampling. the comparison between the actual situation and the assessment results obtained from the NORVI-based assessment method is given in Table 3, which shows that the assessment results based on NORVI are in line with the actual situation. Through the above comparison, we can say that the accuracy of the proposed NORVI-based assessment is higher than that of ORVI-based assessment with respect to the change of NOUVI, and it can provide more reliable assessment results for the production process. To investigate the reasons for the nonoptimal performance grade, the variable contributions to the assessment index are calculated based on ORVI and NORVI-based assessment methods, and the results are shown in Figs. 6 and 7, respectively. In order to avoid the smearing effect on unrelated variables, the variable contributions at the first sampling of general performance grade are used to identify the responsible variables. From Figs. 6 and 7, it can be seen that the contributions of NaCN flow (variable

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

481

Fig. 7 – NORVI-based cause identification results for the first case at the 1st sampling.

Fig. 8 – ORVI-based assessment results for the second case.

Fig. 9 – NORVI-based assessment results for the second case. 2) and CN− concentration (variable 3) are greater than those of others and they are considered as the responsible variables for the nonoptimal performance grade. Since CN− concentration is influenced by NaCN flow, the cause identification result is consistent with the actual situation, which illustrates the effectiveness of the cause identification strategy proposed in this study. In the second case, the assessment results based on ORVI and NORVI are shown in Figs. 8 and 9, respectively. The turning point from optimal to nonoptimal are separately identified at the 1640 and 1585th samplings with ORVI and NORVI-based assessment

482

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

Table 4 – Accuracy rates of ORVI and NORVI-based assessment methods for the second case. Method

Numbers of correct evaluated samplings

Accuracy rates(%)

ORVI NORVI

1594 1670

88.56 92.78

Fig. 10 – ORVI-based cause identification results for the second case at the 1640th sampling.

Fig. 11 – NORVI-based cause identification results for the second case at the 1585th sampling. methods. Although both of the above two assessment methods are able to identify the conversion of the operational performance of the process from optimal to nonoptimal, the NORVI-based assessment method can recognize it much earlier than ORVI-based assessment method for the change of NORVI. Trace it to its cause, main reason is the proposed method in this study uses the processing technique of nonlinear variable correlations, which prompts it to extract the process variation information related to the operational performance more precise for nonlinear processes. Therefore, the accuracy rates of the assessment results based on NORVI is still higher than that of the ORVI-based assessment method for the change of NORVI as shown in Table 4. The cause identification results based on ORVI and NORVI-based assessment methods for the second case are shown in Figs. 10 and 11, respectively, where the variable contributions at the 1640 and 1585th samplings are calculated correspondingly. It can be observed from Figs. 10 and 11 that zinc-powder addition amount (variable 6) and [Au(CN)2 ]− concentration (variable 7) are with greater contributions than other variables. Because [Au(CN)2 ]− concentration is a process variable affected by zinc-powder addition amount, the cause identification result is consistent with the actual situation. By comparing with the ORVI-based assessment method for the above two cases, we can conclude that the proposed NORVIbased assessment method in this paper has higher accuracy and stronger sensitivity for its stronger ability in dealing with the nonlinear variable correlations, and the improved accuracy rate and timeliness of which illustrates our viewpoints.

6.

Conclusion

In this paper, we propose a novel operating optimality assessment and cause identification method based on NORVI for nonlinear industrial processes. The proposed approach firstly extends the linear two-step basis vector extraction algorithm to the nonlinear version to obtain the common variable correlation sub-space. Then the common basis vectors, along which the ampli-

483

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

tudes of process variation information of different performance grades are almost the same, are chosen to extract the NOUVI. Furthermore, the assessment models are developed based on the NORVI of each performance grade. In operating optimality assessment, the operational performance of the process can be evaluated accurately based on the predetermined assessment rules. Furthermore, the variable contribution-based cause identification strategy is developed for the nonlinear process to find the responsible nonoptimal causes. From the simulation results, it can be concluded that the proposed NORVI-based assessment method has higher accuracy and timeliness than ORVI-based assessment for nonlinear processes, which has been confirmed by the simulation data of the gold hydrometallurgical process.

Acknowledgments This work was supported by the National Natural Science Foundation of China (grant numbers 61533007 and 61374146), and State Key Laboratory of Synthetical Automation for Process Industries Fundamental Research Funds (grant number 2013ZCX02-04).

Appendix A. Part one: extraction of common variable correlation In the first step of basis vector extraction, the common basis vector can be found by the following constrained optimization problem:

max R2 = max

 s.t.

C  

T ¯ ¯c pgT˚(X c) ˛

2

c=1

(A1)

¯ p¯ g = 1 pgT ˛¯ Tc ˛¯ c = 1 T

Because the combination coefficient vector ˛¯ c is set to unit length, the sub-basis vector ˚(Xc ) ˛¯ c could have length more than T unit. Thence, the objective function of Eq. (A1) actually represents the covariance between the sub-basis vector ˚(Xc ) ˛¯ c and the common basis vector p¯ g rather than the pure correlation analysis. However, if analyze the correlations directly in this step, i.e., T

change the constraint ˛¯ Tc ˛¯ c = 1 to ˛¯ Tc ˚(Xc )˚(Xc ) ˛¯ c = 1, it would bring some problems difficult to be dealt with, and the detailed explanations are listed in the following second step of basis vector extraction. Then construct Lagrange function, and the objective function can be expressed as an unconstrained extremum problem:

L(p¯ g , ˛¯ c , g , c ) =

C  

T p¯ Tg ˚(Xc ) ˛¯ c

2

− g (p¯ Tg p¯ g − 1) −

c=1

C 

c (˛¯ Tc ˛¯ c − 1),

(A2)

c=1

where g and c are constant scalars. Calculate the derivatives of L(p¯ g , ˛¯ c , g , c ) with respect to p¯ g , ˛¯ c , g and c , set all of them equal to zero, and then the following mathematical expressions are derived: C  

T

T



|p¯ Tg ˚(Xc ) ˛¯ c |˚(Xc ) ˛¯ c − 2 g p¯ g = 0,

∂L/∂p¯ g = 2

(A3)

c=1 T

∂L/∂˛¯ c = 2|p¯ Tg ˚(Xc ) ˛¯ c |˚(Xc )p¯ g − 2 c ˛¯ c = 0,

(A4)

p¯ Tg p¯ g − 1 = 0,

(A5)

˛¯ Tc ˛¯ c − 1 = 0.

(A6)

Left multiply Eq. (A3) with p¯ Tg and Eq. (A4) with ˛¯ Tc respectively, and obtain the following relations in combination with Eqs. (A5) and (A6), i.e., C  

T p¯ Tg ˚(Xc ) ˛¯ c

2

= g ,

(A7)

c=1



T p¯ Tg ˚(Xc ) ˛¯ c

2

= c .

(A8)

Thereafter, Eqs. (A3) and (A4) can be modified as follows: C   c=1

T

c ˚(Xc ) ˛¯ c = g p¯ g ,

(A9)

484

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487



˚(Xc )p¯ g /

c = ˛¯ c .

(A10)

Input Eq. (A10) into Eq. (A9), and then perform the constrained optimization problem (A1) is equivalent to resolving the eigenvalue problem of Eq. (A11), C  



T ˚(Xc ) ˚(Xc ) p¯ g = g p¯ g .

(A11)

c=1

Denote X = [x1 , x2 , ...xN ]T ∈ RN×J as the dataset constituted by the modeling data from all the performance grades and p¯ g =

 N

T

ˇi ˚(xi ) = ˚(X) ˇ as the linear combination of all the observations, where N =

C 

Nc , ˇ ∈ RN×1 . Then left multiply Eq. (A11)

c=1

i=1

with ˚(X), and Eq. (A11) is modified below:



C 

 T

T

T

ˇ = g ˚(X)˚(X) ˇ.

˚(X)˚(Xc ) ˚(Xc )˚(X)

(A12)

c=1 T

T

To facilitate the calculation, define kernel matrix Kc = ˚(X)˚(Xc ) ∈ RN×Nc and K = ˚(X)˚(X) ∈ RN×N , where Kcin = ˚(xi ), ˚(xc,n ) is the i th row and nth column element of Kc , Kij = ˚(xi ), ˚(xj ) is the ith row and jth column element of K, and x, y

denotes the dot product between vectors x and y. Then Eq. (A12) is rewrote as:

K

−1

C 

Kc KTc

ˇ = g ˇ.

(A13)

c=1

To solve the above eigenvalue problem, one can avoid performing the nonlinear mapping ˚(·) and computing the dot product in the feature space F by introducing a kernel function of form k(x, y) = ˚(x), ˚(y) . According to Mercer’s theorem of functional analysis, if the kernel function is a continuous kernel of a positive integral operator, there exists a mapping into a space where a kernel function acts as a dot product. Hence, the requirement on the kernel function is that it satisfies Mercer’s theorem (Courant and Hilbert, 1953). Then R¯ eigenvectors ˇ1 , ˇ2 , ..., ˇR¯ are obtained from Eq. (A13) and used as the combination coefficients to T ¯ ¯ constitute the common basis matrix P¯ g = [p¯ g,1 , p¯ g,1 , ..., p¯ g,R¯ ] = ˚(X) B ∈ Rh×R , where B = ˇ1 , ˇ2 , ..., ˇR¯ ∈ RN×R . In addition, input Eq. (A10) into Eq. (A8), c can be calculated as follows: T

T

T

c = p¯ Tg ˚(Xc ) ˚(Xc )p¯ g = ˇT ˚(X)˚(Xc ) ˚(Xc )˚(X) ˇ = ˇT Kc KTc ˇ.

(A14)

Based on Eq. (A10), the sub-basis vector is calculated below: T p¯ c,j = ˚(Xc ) ˛¯ c,j =



T

1/ c ˚(Xc ) KTc ˇ.

(A15)





T All the sub-basis vectors of performance grade c form the sub-basis matrix P¯ c = p¯ c,1 , p¯ c,2 , ..., p¯ c,R¯ = ˚(Xc ) KTc B(c )

−1

⁄2 with  = c

¯ ¯ diag{ c,1 , c,2 , ..., c,R¯ } ∈ RR×R . From the viewpoint of sample selection, this is also equivalent to picking R¯ representatives out of the original Nc samples while keeping the variable dimension fixed. T ¯ In the second step of basis vector extraction, P¯ c ∈ RR×h is used to replace ˚(Xc ). Construct the following objective function:

max R2 = max

 s.t.

C 

(pTg P¯ c ˛c )

2

c=1

pTg pg = 1 T ˛Tc P¯ c P¯ c ˛c =

(A16)

1

Be different from the objective function of Eq. (A1), Eq. (A16) is constructed with the constraint of P¯ c ˛c  = 1. Therefore, it is inherently excludes the effect of module length of each sub-basis vector and directly maximizes the sum of the correlations between the common basis vector pg and the sub-basis vector of each performance grade. Using a Lagrange operator, the objective function of Eq. (A16) can be substituted by an unconstrained extremum problem: C  

T P¯ c (P¯ c P¯ c )

−1 T P¯ c



pg = g pg .

(A17)

c=1

Here, we can say that if the second step of basis vector extraction is directly implemented on the original measurements ˚(Xc ) T

T

In practice, because the rather than P¯ c , it hereby depends heavily on whether or not each Nc × Nc matrix ˚(Xc )˚(X  c ) is invertible.  T observations in ˚(Xc ) are high-dimensional and highly correlated, it is general that rank ˚(Xc )˚(Xc ) < Nc . Thus, the invertibility

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

485

cannot be achieved so that directly applying eigenvalue decomposition on the original measurements may lead to a rank-deficient T problem. However, sub-basis matrix P¯ c obtained from the first step of basis vector extraction can guarantee the invertibility of T matrix P¯ c P¯ c . Therefore, the first step of basis vector extraction prepares a good initial platform for achieving the final common basis matrix. −1 T ⁄2 into Eq. (A17), and then the common basis vector extraction comes down to a simple analytical Input P¯ c = ˚(Xc ) KTc B(c ) solution.



C 

 T

˚(Xc ) F c ˚(Xc ) pg = g pg ,

(A18)

c=1 −1

T

where F c = KTc B(BT Kc Kc,c KTc B) BT Kc , Kc,c = ˚(Xc )˚(Xc ) . Be similar to the first step of basis vector extraction, the common basis vector can be expressed as the linear combination of T T ¯ i.e., the observations in P¯ = [P¯ 1 , P¯ 2 , ..., P¯ C ] = [p¯ 1 , p¯ 2 , ..., p¯ R ]T ∈ RR×h ,R = CR,

pg =

R 

−1

˜ ¯ = ˚(X)T K˜ B

dr p¯ r = Pd

⁄2 d,

(A19)

r=1 ˜

˜

˜ = CN. ˜ = diag{KT1 , ...KT } ∈ RN×N ,B˜ = diag{B, ..., B} ∈ RN×R ,  = diag{1 , ..., C } ∈ RR×R , N where d = [d1 , ..., dR ]T ∈ RR×1 , K C

   C

⁄2 )T ˚(X) , one can get that,

−1

˜ Insert Eq. (A19) into Eq. (A18), and left multiply both sides of Eq. (A18) with (K˜ B

−1

−1 −1 T ⁄2 )T KK˜ B

⁄2 ] (K˜ B

⁄2 ) ( ˜ ˜

−1

˜ [(K˜ B

C 

−1

˜ Kc Fc KcT K˜ B

⁄2 )d = g d.

(A20)

c=1

−1

T

−1

−1

⁄2 ) KK˜ B

⁄2 ] is invertible, and one can find the solutions of Eq. (A20) by solv˜ ˜ ˜ B˜ and  are full column rank, [(K˜ B

Because K, ing the eigenvalue problem. Orderly, a total of A common basis vectors can be derived by Eq. (A20) in accord with the descending g , T ˜ −1⁄2 and constitute the common basis matrix Pg = [pg,1 , pg,2 , ..., pg,A ] = ˚(X) K˜ B

D ∈ Rh×A with D = [d1 , d2 , ..., dA ] ∈ RR×A . Thereafter, ˚(Xc ) can be divided into two parts: ˚(Xc )

ˆ c ) + ˚(X ˜ c) = ˚(X

(A21)

= ˚(Xc )Pg PTg + (Xc )(I − Pg PTg ),

ˆ c ) = ˚(Xc )Pg PTg includes the process variation information driven by the common variable correlation across all perforwhere ˚(X ˜ c ) = ˚(Xc )(I − Pg PTg ) contains unique variable correlation specific to performance grade c. mance grades, while ˚(X Before applying the two-step basis vector extraction, mean centering in the high dimensional feature space should be ¯ i.e., performed. This can be done by substituting the kernel matrix Kc with K¯ c and K with K, K¯ c = Kc − 1NKc − Kc 1Nc + 1NKc 1Nc ,

(A22)

K¯ = K − 1NK − K1N + 1NK1N,

(A23)

where 1N and 1Nc represent the N × N and Nc × Nc matrices whose elements are 1/N and 1/Nc , respectively.

Part two: extraction of common variation information amplitude ˜ c ) are dissimilar for different performance grades, it is actually a part of the NORVI of perSince the variable correlation in ˚(X formance grade c; additionally, it necessary to measure the amplitudes of the variation information of ˚(Xc ) along each common basis vector, and separate the common basis vectors with approximately equal amplitudes from Pg to constitute a sub-space containing the common variable correlation and variation information amplitude. The variation information of ˚(Xc ) along pg,a can be calculated as follows: −1

T ˜ tc,a = ˚(Xc )pg,a = ˚(Xc )˚(X) K˜ B

−1 ⁄2 da = KT K˜ B

⁄2 da . ˜

c

(A24)

The common basis vectors with approximately equal amplitudes can be found through the following steps. (i) Calculate the amplitude of variation information. Without loss of generality, the amplitude of t c,a is denoted as f (tc,a ), a = 1, 2, ..., A, where f (tc,a ) can be the mean or median of the elements in tc,a , or some other operators characterizing the amount of information carried by t c,a . (ii) Define an index to measure the equality of the amplitude. Chose performance grade C as a reference, and calculate the ratio c,a = f (tc,a )/f (tC,a ), c = 1, 2, ..., C − 1. Then given a relaxation factor ı(0 < ı < 1), if 1,a , 2,a , ..., C−1,a are all within [1 − ı, 1 + ı],

486

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

it means that f (t1,a ), ..., f (tC,a ) are approximately equal; otherwise, the equality does not hold. The value of ı can be set by trial and error through cross validation in order that the error rate and missing alarms are minimized in evaluation. (iii) Decompose the common variable correlation sub-space Pg . Pick out the common basis vectors satisfying the condi˘

˘

˘

˘

tion of 1 − ı ≤ 1,a , 2,a , ..., C−1,a ≤ 1 + ı from Pg and denote them as pg,1 , pg,2 , ..., pg,˘A . They construct the sub-space pg = ˘

˘

˘

[pg,1 , pg,2 , ..., pg,˘A ] containing the common variable correlation and common variation information amplitude for all the 





performance grades. The rest common basis vectors are rewritten as p g,1 , p g,2 , ..., p g,A and drive the sub-space  Pg = [ pg,1 ,  pg,2 , ...,  pg,A ],  A = A −˘A, in which, the common variable correlation but unequal variation information amplitude present in different performance grades. Through the above procedure, ˚(Xc ) is formulated below: T

= (Xc )˘Pg˘Pg + (Xc ) Pg  PTg + (Xc )(I − Pg PTg )

(Xc )

˙ c ), = ˘˚(Xc )+˚(X

(A25)

T

˙ c ) = ˚(Xc )( Pg  PTg + I − Pg PTg ) = ˚(Xc )(I − where˘˚(Xc ) = ˚(Xc )˘Pg˘Pg is the NOUVI shared by all of the performance grades and ˚(X T ˘Pg˘Pg ) is

the NORVI unique to performance grade c.

Appendix B. ˙ c ) is given below. The detailed PCA decomposition procedure on ˚(X ˙ c ), i.e., Firstly, calculate the covariance matrix of ˚(X ˙ c )T ˚(X ˙ c )/Nc . c = ˚(X

(B1)

Then, one has to solve the eigenvalue problem in the feature space F as follows: ˙ ˙ c p˙ = c p,

(B2)

where ˙ c is the eigenvalue larger than zero and p˙ is the corresponding eigenvector. ˙ c )T q, Eq. (B2) can be further written as follows: Since exist coefficient q = [q1 , q2 , ...qNc ]T ∈ RNc ×1 such that p˙ = (X ˙ c )T q = c ˚(X ˙ c )T q. ˙ c (X

(B3)

˙ c ), and Eq. (B3) is modified as: Left multiply both sides of Eq. (B3) with (X 2 ˙ c Nc K˙ c,c q = K˙ c,c q,

(B4)

˙ c ) and can be calculated below: where K˙ c,c is the kernel matrix of (X K˙ c,c

˙ c )˚(X ˙ c) = ˚(X

T T

T

= (Xc )(I −˘Pg˘Pg )(I −˘Pg˘Pg )T (Xc ) = (Xc )(I =

T

T T −˘Pg˘Pg )(Xc )

(B5)

−1 T T T −1 2 ˜ B ˜ ˜ (X)(Xc )T Kc,c − (Xc )(X)T K ˘D˘D  2 B˜ K

⁄

⁄

= Kc,c − KTc˘WKc −1

−1

⁄2˘D˘DT ⁄2 B˜ T K˜ T .˘Dconsists of several vectors of D corresponding to˘Pg . ˜ with ˘W = K˜ B

Performing eigenvalue decomposition on Eq. (B2) is equivalent to resolving the following eigenvalue problem: ˙ c Nc q = K˙ c,c q

(B6)

From Eq. (B6), Mc eigenvectors q1 , q2 , ..., qMc corresponding to the first Mc maximum and nonzero eigenvalues are obtained and compose the eigenvector matrix Q c = [q1 , q2 , ..., qMc ] ∈ RNc ×Mc . Then the loading matrix P˙ c , the score matrix T˙ c and the residual matrix E˙ c are formulated as follows: T

˙ c )T Q c = (I −˘Pg˘PTg ) (Xc )T Q c , P˙ c = (X T˙ c

(B7)

˙ c )P˙ c = (X T

T T

T

= (Xc )(I −˘Pg˘Pg )(I −˘Pg˘Pg ) (Xc ) Q c =

(Kc,c − KTc WKc )Q c ,

(B8)

chemical engineering research and design 1 1 7 ( 2 0 1 7 ) 472–487

E˙ c

T ˙ c ) − T˙ c P˙ c = (X

= (Xc ) − (Kc,c − 2KTc WKc + KTc WKWKc )Q c Q Tc (Xc ) + (Kc,c Q c Q Tc − 2KTc WKc Q c Q Tc − I + KTc WKWKc Q c Q Tc )KTc W(X).

487

(B9)

References Courant, R., Hilbert, D., 1953. Methods of Mathematical Physics, Vol. 1. Interscience, New York. Del Rio-Chanona, E.A., Zhang, D., Vassiliadis, V.S., 2016. Model-based real-time optimisation of a fed-batch cyanobacterial hydrogen production process using economic model predictive control strategy. Chem. Eng. Sci. 142, 289–298. Du, M., Scott, J., Mhaskar, P., 2013. Actuator and sensor fault isolation of nonlinear process systems. Chem. Eng. Sci. 104, 294–303. Edgar, T.F., Himmelblau, D.M., Lasdon, L.S., 2001. Optimization of Chemical Processes. McGraw-Hill, New York. Haykin, S., 1999. Neural Networks. Prentice-Hall, Englewood Cliffs, NJ. Hu, Y., Ma, H.H., Shi, H.B., 2013. Enhanced batch process monitoring using just-in-time-learning based kernel partial least squares. Chemom. Intell. Lab. Syst. 123, 15–27. Hwang, J.H., Lee, J.C., Roh, M.I., Lee, K.Y., 2013. Determination of the optimal operating conditions of the dual mixed refrigerant cycle for the LNG FPSO topside liquefaction process. Comput. Chem. Eng. 49, 25–36. Lee, J.M., Yoo, C.K., Choi, S.W., Vanrolleghem, P.A., Lee, I.B., 2004. Nonlinear process monitoring using kernel principal component analysis. Chem. Eng. Sci. 59, 223–234. Liu, Y., Chang, Y.Q., Wang, F.L., 2014a. Online process operating performance assessment and nonoptimal cause identification for industrial processes. J. Process Control 24, 1548–1555. Liu, Y., Chang, Y.Q., Wang, F.L., Ma, R.C., Zhang, H.L., 2014b. Complex process operating optimality assessment and nonoptimal cause identification using modified total kernel PLS. In: Proceedings of the 26th Chinese Control and Decision Conference, Changsha, China, May 31–June 2, pp. 1221–1227. Liu, Y., Wang, F.L., Chang, Y.Q., 2016. Operating optimality assessment based on optimality related variations and nonoptimal cause identification for industrial processes. J. Process Control 39, 11–20. Mori, J., Yu, J., 2014. Quality relevant nonlinear batch process performance monitoring using a kernel based multiway non-Gaussian latent sub-space projection approach. J. Process Control 24, 57–71. Peng, K.X., Zhang, K., Li, G., Zhou, D.H., 2013. Contribution rate plot for nonlinear quality-related fault diagnosis with application to the hot strip mill process. Control Eng. Pract. 21, 360–369. Qin, S.J., 2012. Survey on data-driven industrial process monitoring and diagnosis. Annu. Rev. Control 36, 220–234. Qin, S.J., Valle, S., Piovoso, M.J., 2001. On unifying multiblock analysis with application to decentralized process monitoring. J. Chemom. 15, 715–742. Rossi, F., Copelli, S., Colombo, A., Pirola, C., Manenti, F., 2015. Online model-based optimization and control for the combined optimal operation and runaway prediction and prevention in (fed-)batch systems. Chem. Eng. Sci. 138, 760–771. Schölkopf, B., Smola, A., Müller, K.R., 1998. Nonlinear component analysis as a kernel eigenvalue problem. Neural Comput. 10, 1299–1319. Srour, M., Gomes, V.G., Altarawneh, I., 2008. Optimal operating strategies for emulsion terpolymerisation. Chem. Eng. Sci. 63, 4257–4268. Ye, L.B., Liu, Y.M., Fei, Z.S., Liang, J., 2009. Online probabilistic assessment of operating performance based on safety and optimality indices for multimode industrial processes. Ind. Eng. Chem. Res. 48, 10912–10923. Zhao, C.H., Gao, F.R., 2015. Online fault prognosis with relative deviation analysis and vector autoregressive modeling. Chem. Eng. Sci. 138, 531–543. Zhao, C.H., Yao, Y., Gao, F.R., Wang, F.L., 2010. Statistical analysis and online monitoring for multimode processes with between-mode transitions. Chem. Eng. Sci. 65, 5961–5975. Zhao, C.H., Gao, F.R., Niu, D.P., Wang, F.L., 2011. A two-step basis vector extraction strategy for multiset variable correlation analysis. Chemom. Intell. Lab. Syst. 107, 147–154. Zhou, D.H., Li, G., Qin, S.J., 2010. Total projection to latent structures for process monitoring. AIChE J. 56, 168–178.