Performance-relevant kernel independent component analysis based operating performance assessment for nonlinear and non-Gaussian industrial processes

Chemical Engineering Science 209 (2019) 115167

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

Performance-relevant kernel independent component analysis based operating performance assessment for nonlinear and non-Gaussian industrial processes Yan Liu a,c,⇑, Fuli Wang a,b, Yuqing Chang a,b, Furong Gao c, Dakuo He a,b a b c

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 Department of Chemical and Biological Engineering, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong

h i g h l i g h t s
A novel PR-KICA method is proposed to overcome the defects of conventional methods.
PR-KICA is applied to operating performance assessment of industrial processes.
Two new criteria are designed to rank the order and determine the number of ICs.
Online assessment results include both performance grades and conversions.
Cause identification result helps to improve the process operating performance.

a r t i c l e

i n f o

Article history: Received 15 August 2018 Received in revised form 5 August 2019 Accepted 21 August 2019 Available online 23 August 2019 Keywords: Comprehensive economic index Industrial processes Nonoptimal cause identification Operating performance assessment Performance-relevant kernel independent components

a b s t r a c t The operating performance assessment of industrial processes becomes increasingly important in manufacturing production. A novel operating performance assessment method based on performancerelevant kernel independent component analysis is proposed here for nonlinear and non-Gaussian processes. The proposed method accounts for the comprehensive economic index, and the objectives are simultaneously to maximize the non-Gaussianity of independent components as well as the correlations between them and the comprehensive economic index. When applying it to online assessment, it demonstrates stronger robustness and higher sensitivity than the traditional methods do, which is attributed to its capacity in highlighting the performance-relevant variation information in modeling. Furthermore, both the performance grades and the conversions can be evaluated, which enhances the interpretability of the results. For the nonoptimality, the variable contributions are used to find the possible cause. Finally, the efficiency of the proposed method is illustrated by a case of gold hydrometallurgy process. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction In order to pursue higher economic benefits, the process operating optimization is a crucial issue hence has attracted tremendous attention both in the academic and the industrial environments (Edgar et al., 2001; Srour et al., 2008; Hwang et al., 2013; del Rio-Chanona et al., 2016; Wang et al., 2017). In general, processes usually operate at some optimum level in the early stage of the production. However, along with the flow of time, the pro⇑ Corresponding author at: College of Information Science & Engineering, Northeastern University, 3 Lane 11, Wenhua Road, Heping District, Shenyang, Liaoning 110819, China. E-mail address: [email protected] (Y. Liu). https://doi.org/10.1016/j.ces.2019.115167 0009-2509/Ó 2019 Elsevier Ltd. All rights reserved.

cess operating performance may deviate from the optimum level due to process disturbances and varying model uncertainties, which often cancels the benefits of preliminary designs for process optimization and results in degraded operating performance. It is necessary therefore to implement an effective operating performance assessment strategy of industrial processes in order to ensure improvements of the operating performance and economic benefits. In the past several decades, many process monitoring methods, such as principal component analysis (PCA) (Jackson and Wiley, 1991; Dunia et al., 1996) and Fisher discriminant analysis (FDA) (Chiang et al., 2000; Chiang et al., 2004), have been applied to industrial processes (Plovoso and Kosanovich, 1994; Kourti and MacGregor, 1996; Qin, 2003). Nevertheless, they are inappropriate

2

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

for quality-relevant process analysis because the underlying models do not incorporate quality variables. For this reason, some quality-relevant monitoring method are developed (Zhou et al., 2010; Qin and Zheng, 2013; Haghani et al., 2014; Yin et al., 2015; Wang et al., 2017), and partial least squares (PLS) (Mejdell and Skogestad 1991) is one of the most popular approaches. However, since PLS focuses on the linear dependencies between process variables and only takes second-order statistics into account, it lacks the ability to deal with nonlinear and non-Gaussian data, To deal with this issue, kernel independent component analysis (KICA) (Kocsor and Csirik, 2001; Lee et al., 2007) is extended to process monitoring. KICA can extract not only nonlinear features but also non-Gaussian features through maximizing the approximation of negentropy within the high-dimensional feature space. Nevertheless, it belongs to quality-irrelevant modeling strategy for not considering the quality information when estimating independent components (ICs). Although the quality-relevant modeling strategy is rarely reported for the process with the coexistence of nonlinearity and non-Gaussianity, the pioneering work has provided abundant theoretical bases for our following work. It is well known that the main task of process monitoring is to distinguish whether the process is operating under ‘‘normal” or ‘‘fault” conditions, however, it cannot satisfy the quest by enterprises for profits any longer. The process operating performance assessment is to get a measure on how far is the gap between the current operating condition and the optimal operating level (Ye et al., 2009; Liu et al., 2014, 2016, 2015; Sedghi and Huang, 2017; Zou et al., 2017; Chu et al., 2018), and the premise of which is that the operating condition must be normal. According to process characteristics and plant personnel’s attitudes of the operating performance, the performance level in the normal operating range can be further divided into several grades including optimal, suboptimal, general, and poor, etc., and the conversion process between two grades. Thereby, managers and operators can grasp the operating level by operating performance assessment and propose suggestions on operating adjustment and performance improvement. Despite a rich body of literature in process monitoring, studies on operating performance assessment of industrial processes are still in its infancy. Generally, the process operating performance is closely associated with the comprehensive economic index (CEI), and higher CEI usually corresponds to better operating performance. Since it is often difficult to calculate the CEI online, the CEI-prediction based assessment strategies are developed (Ye et al., 2009; Liu et al., 2015) for Gaussian and non-Gaussian processes, respectively. Taking as an alternative, the assessment methods based on process variation information have been developed for linear and nonlinear industrial processes (Liu et al., 2014, 2017; Chu et al., 2018). These methods take the relationship between process variables and the CEI into account so as to establish the model to infer the CEI by indirectly evaluating the process variation information. Thus, they are also called performance-relevant variation information based assessment methods. Nevertheless, since existing techniques are oriented to the process with only nonlinearity or nonGaussianity, their effectiveness will be discounted when facing to the process with the coexistence of both nonlinearity and nonGaussianity. For nonlinear and non-Gaussian industrial processes, a novel operating performance assessment method based on performance-relevant KICA (PR-KICA) is proposed here. The proposed assessment strategy consists of three aspects: offline modeling, online assessment, and nonoptimal cause identification. The proposed method belongs to the performance-relevant variation information based assessment strategies. Hence, offline modeling is to extract performance-relevant variation information

from the training data of different performance grades so as to provide the reference standard for online assessment. PR-KICA is as an improvement of the traditional KICA and is committed to simultaneously maximize the non-Gaussianity of independent components (ICs) and the correlations between ICs and the CEI. This ensures that the ICs estimated by PR-KICA can not only characterize the nonlinearity and non-Gaussianity of the training data but also reflect the characteristics of process operating performance. In addition, to ensure the accuracy of the assessment model, retaining a small number of significant ICs rather than all is more suitable for the subject of operating performance assessment. In view of this, two new criteria, named as contribution of IC and cumulative percent contribution of ICs, are designed to rank the order and determine the number of ICs, respectively. Online assessment is to evaluate the current operating performance with the new observed data in real time. Since there are multiple performance grades and conversion processes, it is essentially a multi-class classification problem. For this, many methods can be used in theory, such as Support Vector Machines (Vapnik, 1998), Fuzzy Clustering (Yang, 1993), etc. However, because there is no assessment model of conversion processes, the online assessment strategies based on the methods mentioned above will lack generalization ability for new samples from the conversion process and are liable to cause misclassification in online assessment. Being similar to the online process monitoring, i.e., constructing statistics to achieve the online classification, our strategy is to construct an assessment index based on the similarities of the performance-relevant variation information for online assessment. Then even if there is no offline assessment model of conversion processes, one can also judge the operating performance according to the trend and value of the assessment index. For the nonoptimal operating performance, the cause identification strategy based on variable contributions is developed to find the possible cause. The idea of variable selection is used to define the variable contributions. The identified cause variables provide important reference for subsequent operation adjustment and process performance improvement. The contributions of this paper are summarized as follows: (i) the PR-KICA method is propose to deal with the operating performance assessment of nonlinear and non-Gaussian industrial processes; (ii) two new criteria are designed to rank the order and determine the number of ICs based on the correlations between ICs and the CEI; (iii) the assessment results include not only the performance grade but also the conversion process, which is more closely fit the actual situation; (iv) to identify the nonoptimal cause, a variable contribution calculation method based on the idea of variable selection is developed. The remainder of this paper is organized as below. Section 2 briefly reviews some preparatory theoretical supports. Then the proposed PR-KICA based operating performance assessment method is introduced in Section 3. In Section 4, the effectiveness of the proposed assessment approach is demonstrated by a case of gold hydrometallurgy process, and its superiority is also verified as opposed to PLS and KICA approaches. Finally, the paper ends with some conclusions and acknowledgements.

2. Preliminaries 2.1. Kernel independent component analysis The traditional KICA is derived to estimate ICs so as to identify process nonlinearity and non-Gaussianity. The basic idea is to nonlinearly map the process data into a high-dimensional feature space where the data have a more linear structure (Schölkopf et al., 1998) and then perform ICA in that high-

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

dimensional feature space to make the estimated ICs as independent as possible. Denote the process data as X ¼ ½x1 ; x2 ; :::; xN  2 RJN , where N is the number of training samples and J is the number of process variables. Given a nonlinear mapping U : RJ ! F, the original process data X is transformed into the high-dimensional feature space F as H ¼ ½Uðx1 Þ; Uðx2 Þ; :::; UðxN Þ, where the dimension of feature space F can be arbitrarily large or even infinite. Then the covariance matrix C F in the feature space is expressed as N 1X 1 C ¼ Uðxn ÞUðxn ÞT ¼ HHT : N n¼1 N F

ð1Þ

According to the KPCA method proposed by Schölkopf et al. (1998), we can find principal components using the ‘‘kernel trick”

3

addition, it should be calculated quickly and consider the order in which the components are estimated. More details about how to wisely choose the function G are introduced in literatures (Hampel et al., 1986; Hyvärinen, 1998, 1999), where they also suggested a number of commonly used functions for G. The detailed procedures of FastICA algorithm can be found in Refs. (Hyvärinen, 1999; Hyvärinen and Oja, 2000). When the algorithm converges, we can obtain the demixing matrix W ¼ ½w1 ; w2 ; :::; wD T . The ICs are then calculated from

sn ¼ W T zn ¼

pffiffiffiffi T 1 T NW KD A kðxn Þ ¼ Ukðxn Þ;

ð6Þ

pffiffiffiffi T where U ¼ NW T K1 D A . 2.2. Nonlinearity coefficient based on mutual information

instead of eigenvalue decomposition on C F directly. The Gram matrix K can be expressed as K ¼ HT H 2 RNN by ½K ij ¼    Uðxi Þ; U xj ¼ kðxi ; xj Þ, where kð; Þ is a kernel function, and then the nonlinear mapping is avoided through computing inner products in the feature space. This step leads to computational savings when the dimensionality of the feature space is large compared to the number of samples. A commonly used kernel function is Gaussian radial basis function kðx; yÞ ¼ expðjjx  yjj2 =2rÞ; and r is the Gaussian kernel width. Then the mean centered Gram matrix K and mean centered kernel vector kðxÞ ¼ HT UðxÞ ¼ ½kðx1 ; xÞ; . . . ; kðxN ; xÞT can be obtained as in Ref. (Zhang, 2008). Implement eigenvalue decomposition on the mean centered Gram matrix K

ka ¼ K a:

ð2Þ

Then we obtain D largest positive eigenvalues k1 P k2 P . . . P kD of K and their corresponding eigenvectors a1 ; a2 ; . . . ; aD . Accordingly, the eigenvalues and eigenvectors of C F can be expressed as k1 =N P k2 =N P . . . P kD =N and b1 ; b2 ; . . . ; bD , respectively. Here, bd is calculated as follows:

1 bd ¼ pffiffiffiffiffi Had ; d ¼ 1; 2; . . . ; D: kd

ð3Þ

Furthermore, they construct the eigenvector matrix B ¼ ½b1 ; b2 ; . . . ; bD  ¼ HAK1=2 with A ¼ ½a1 ; a2 ; . . . ; aD  2 RND and KD ¼ D diagðk1 ; k2 ; . . . ; kD Þ. pffiffiffiffi Let P ¼ BðKD =N Þ1=2 ¼ N HAK1 D , then the mapped data in feature space F can be whitened by the following transformation:

zn ¼ P T Uðxn Þ ¼

pffiffiffiffi 1 T NKD A kðxn Þ;

ð4Þ

where Efg denotes the expectation, and zn satisfies Efzn zTn g ¼ I. Define a demixing matrix W 2 RDD so that the estimated ICs, T

sn ¼ W zn , are as independent of each other as possible. With the requirement of Efsn sTn g ¼ I, it is clear that the demixing matrix W is an orthogonal matrix and satisfies W T W ¼ I. To calculate W, Hyvärinen (1999) and Hyvärinen and Oja (2000) introduced a simple and efficient FastICA algorithm. In FastICA algorithm, each column vector w of W is randomly initialized and updated iteratively so that the IC has the maximized non-Gaussianity. The objective function is an approximation of negentropy as follows:

JðyÞ ¼ ½EfGðyÞg - EfGðv Þg2 ;

ð5Þ

where y is assumed to be of zero mean and unit variance; v is a Gaussian variable of zero mean and unit variance; G is any nonquadratic function (Hyvärinen, 1999; Hyvärinen and Oja, 2000). To obtain a good approximation of negentropy, the function G should be chosen to ensure that each column vector w has the characteristics of consistency, asymptotic variance, and robustness. In

In order to select an appropriate data analysis method, it is necessary to analyze the data characteristics in advance, such as the dependency between variables, the distribution characteristic of data, etc. The dependency between variables is to test whether the process characteristic is linear or nonlinear. Pearson correlation coefficient is one of the most commonly used coefficient

covðX; YÞ

qX;Y ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; varðXÞ varðYÞ

ð7Þ

where X and Y are two discrete random variables concerned, and the value of qX;Y is between 1 and 1. Pearson correlation coefficient can well represent the linear dependencies between variables, but it is not sensitive to the nonlinear dependencies. That is, even if qX;Y ¼ 0, it is not able to say that two variables are independent because they may have nonlinear dependencies. Mutual information (Wang et al., 2005) is a measurement of the general dependency between two random variables, which is defined as

IðX; YÞ ¼ HðXÞ þ HðYÞ  HðX; YÞ;

ð8Þ

where HðXÞ and HðX; YÞ are the information entropy of X and the joint information entropy between X and X, respectively. Although the mutual information is sensitive to any dependencies, it cannot distinguish whether the variable dependency is linear or nonlinear. Thus, Zhang et al. (2016) proposed a new nonlinearity coefficient, combining Pearson correlation coefficient and the mutual information, i.e.,

  rX;Y ¼ v X;Y ð1  qX;Y Þ;

ð9Þ

where v X;Y ¼ IðX; YÞ=minðHðXÞ; HðYÞÞ is normalization of the mutual information. The values of nonlinearity coefficient r X;Y are located in the interval of [0, 1]. It indicates the degree of nonlinear dependencies only. If the nonlinearity coefficient r X;Y is close to 1, the variables have strong nonlinearly with each other; otherwise, if r X;Y approaches to 0, the variables have weak nonlinear dependencies. To measure the nonlinear dependencies of variable population, a population nonlinearity coefficient is constructed as follows:

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi uPJ u i;j¼1 r 2X ;X t i j : r¼ J2  J

ð10Þ

The value of r is also between 0 and 1. The closer r is to 1, the stronger the dependency is; on the contrary, the weaker the dependency is. Given a threshold Trð0 < Tr < 1Þ, and if the nonlinearity coefficient r is greater than the threshold Tr, nonlinear method should be chosen for process analysis. The threshold is always defined based on practical experience and the specialist’s

4

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

knowledge. In this study, the nonlinearity coefficient proposed by Zhang et al. (2016) is used to test the dependency between variables before process analysis. 2.3. F-straight method based on Mahalanobis distance To test the distribution characteristic of data, the F-straight method based on Mahalanobis distance (FSMD) (Zhang et al., 2016) is used in this study. The basic idea of FSMD is that, when the data approximately follow a multivariate Gaussian distribution, the distribution function of the squared Mahalanobis distances of samples, DðnÞ , is a specific F distribution. For a random matrix X ¼ ½x1 ; x2 ; :::; xN T 2 RNJ , if the samples follow a Gaussian distribution X  NJ ðl; RÞ and the population mean and population covariance matrix are unknown, the squared Mahalanobis distances of samples follow F distribution  T



Dðxn ; lÞ ¼ ðxn  lÞ S þ ðxn  lÞ 

pðN2  1Þ Fðp; N  pÞ; NðN  pÞ

ð11Þ

t  0:5 ¼ rn ðn ¼ 1; 2; :::; NÞ: N

ð12Þ

In addition, the quantile of F distribution corresponding to probability r n is given by

Fn ¼

8 <

JðN2 1Þ F ðJ; N NðNJÞ r n

: pðN

2

1Þ F ðp; N NðNpÞ r n

 JÞ; rankðSÞ ¼ J  pÞ; rankðSÞ ¼ p

:

ð13Þ

If hypothesis test H0 is valid, DðnÞ  F n . Thus, the plot of ðDðnÞ ; F n Þ should scatter on a line passing through the original point with a slope of 1. Let the linear regression equation between DðnÞ and     P F n be F n ¼ a þ bDðnÞ , where a ¼ F b D, b ¼ Nn¼1 ðDðnÞ  DÞðF n  F Þ=  2   PN P PN and F ¼ Nn¼1 F n =N. To test n¼1 ðDðnÞ  DÞ ,D ¼ n¼1 DðnÞ =N whether the linear regression equation reflects the linear dependencies between DðnÞ and F n , perform the significance test on it, and the regression standard deviation s is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  PN 2 s¼ n¼1 ðF n  ða þ bDðnÞ ÞÞ =ðN  2Þ. If s= F > 0:15, the significance test is invalid, and reject the null hypothesis. Otherwise, the significance test is valid and the linear regression equation can reflect the real dependencies between DðnÞ and F n . Furthermore, compare the intercept a with 0 and regression coefficient b with 1. It is usually considered that the regression line passes through the original point with a slope of 1 as long as the following conditions are satisfied

(



jaj < r

;

j b  1j < b 









s= F 6 0:15, jaj < r, and jb  1j < b are simultaneously fulfilled, the data follow a Gaussian distribution. 3. PR-KICA based operating performance assessment of nonlinear and non-Gaussian industrial processes Because the traditional KICA algorithm takes no account of the relationship between ICs and the process operating performance, it is difficult for the assessment models to reflect the essential characteristic of the process operating performance. In view of this, a novel process operating performance assessment method based on PR-KICA is developed for nonlinear and non-Gaussian industrial processes. The framework of the assessment method consists of three parts: offline modeling, online assessment, and nonoptimal cause identification. In this section, we will introduce these contents in detail. 3.1. Offline assessment model establishment based on PR-KICA



where l and S are the estimated values of population mean and population covariance matrix, respectively; S+ is the MoorePenrose pseudoinverse (Mardia, 1977) of the covariance matrix S; p ¼ rankðSÞ < J. It has been further proven that D is independent of the form of the pseudoinverse. Since the distribution function can be replaced by its empirical distribution functions, the multivariate normality test can be transformed into the problem of whether the empirical distribution function of statistic D is equivalent to a specific F distribution. First, the order statistics of the squared Mahalanobis distance are Dn : Dð1Þ 6 Dð2Þ 6 ::: 6 DðNÞ : Then the empirical distribution function of Dn is

F N ðDðnÞ Þ ¼

and the data are non-Gaussian. In summary, only if conditions

ð14Þ

where r and b are the thresholds. If the conditions in Eq. (14) are satisfied, the hypothesis test that the data follow a Gaussian distribution should not be rejected. Otherwise, reject the null hypothesis,

Assuming

that

there

are

C

performance

grades,

Xc ¼

½xc;1 ; xc;2 ; :::; xc;Nc  2 RJNc is the process data of performance grade c, c ¼ 1; 2; :::; C, and N c is the number of training samples. yc ¼ ½yc;1 ; yc;2 ; :::; yc;Nc  2 RNc 1 is composed of samples of the CEI corresponding to X c . The CEI may be one of important production indices, such as cost, profits, total revenue and production efficiency, etc., or a weighted integration of several important production indices. X c and yc together constitute the training data of performance grade a. In process operating performance assessment, the estimated ICs should not only have the maximum non-Gaussianity but also the greatest correlations with the CEI. Thus, the objective of PR-KICA method consists of two parts, one is to maximize the nonGaussianity of ICs, and the other is to maximize the correlations between ICs and the CEI. As in the original FastICA algorithm, maximizing the nonGaussianity can be realized by the approximation of negentropy, and the objective function is as follows:

 2 max J1 ðwc Þ ¼ EfGðwTc zc;n Þg - EfGðvÞg ; s:t: wTc wc ¼ 1

ð15Þ

where wc is a column vector of the demixing matrix W c 2 RDc Dc , Dc is the number of reserved eigenvectors in KPCA whiting, and pffiffiffiffiffiffi T zc;n ¼ Nc K1 Dc Ac kðxn;c Þ is the whitening score vector, n ¼ 1; 2; :::; N c . With the development of statistical techniques, many statistics have been proposed to measure the correlations between random variables, such as Pearson product-moment covariance and correlation, F-statistics, energy statistics (Székely and Rizzo, 2013), and so on. In this study, the covariance is used. Then the objective is to find unit vector w so that the correlation between the IC and the CEI is maximized,

max J2 ðwc Þ ¼ cov2 ðwTc zc;n ; yc;n Þ s:t: wTc wc ¼ 1;

ð16Þ

where covð; Þ is the covariance between two random variables. Based on the above analysis, the objective function of the PRKICA algorithm can be constructed by transforming J 1 ðwc Þ and J 2 ðwc Þ into an aggregated form as follows:  2 maxJðwc Þ ¼ q1 EfGðwTc zc;n Þg  EfGðvÞg þ q2 cov2 ðwTc zc;n ; yc;n Þ ð17Þ s:t: wTc wc ¼ 1; where qi ¼ bi =sf i ði ¼ 1; 2Þ is the weighting factor of the ith objective. bi 2 ½0; 1 is a preference factor and reveals the relative

5

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

importance of each objective and satisfies b1 þ b2 ¼ 1. sf i is an adaptive scaling factor, i.e., the magnitude of the gradient vector, which is used to control and adjust the convergence speed in real time. By constructing a Lagrange function, the original constrained extremum problem as in Eq. (17) is transformed to an unconstrained extremum problem, i.e.,

 2 Lðwc ; kc Þ ¼ q1 EfGðwTc zc;n Þg - EfGðvÞg þ q2 cov2 ðwTc zc;n ; yc;n Þ þ kc ðwTc wc - 1Þ; ð18Þ where kc is a Lagrange multiplier. Calculate the first partial derivatives of Lðwc ; kc Þ with respect to wc and kc , and set them equal to zero, one can obtain the following mathematical expressions:





rLwc ¼ 2q1 EfGðwTc zc;n Þg - EfGðv Þg  Efzc;n gðwTc zc;n Þg þ2q2 covðwTc zc;n ; yc;n Þ 

PNc



n¼1



ðzc;n - z c Þðyc;n yc Þ Nc - 1

þ 2kc wc ¼ 0;

ð19Þ



J F ðwc Þ ¼ 2q1 EfGðwTc zc;n Þg - EfGðv Þg  Efzc;n zTc;n g 0 ðwTc zc;n Þg  T þ2q1 Efzc;n gðwTc zc;n Þg Efzc;n gðwTc zc;n Þg PNc PNc T     2q ðzc;n - z c Þðyc;n yc Þ ðzc;n - z c Þðyc;n yc Þ n¼1 n¼1 þ 2 þ 2kc I; 2 ðN - 1Þ c

ð27Þ 0

where g is the second-order derivative of G. To simplify the inversion of J F ðwc Þ, we have to approximate the first term in Eq. (27). Since the data is sphered, a reasonable approximation as in the original FastICA algorithm (Hyvärinen, 1999) is

Efzc;n zTc;n g 0 ðwTc zc;n Þg  Efzc;n zTc;n gEfg 0 ðwTc zc;n Þg ¼ Efg 0 ðwTc zc;n ÞgI:

ð28Þ

Then J F ðwc Þ can be simplified as follows:

   J F ðwc Þ ¼ 2q1 EfGðwTc zc;n Þg - EfGðv Þg  Efg 0 ðwTc zc;n Þg þ 2kc I  T þ2q1 Efzc;n gðwTc zc;n Þg Efzc;n gðwTc zc;n Þg PNc  T PNc     2q2 ðzc;n - z c ÞðyTc;n yc Þ ðzc;n - z c Þðycn yc Þ n¼1 n¼1 þ : ðN - 1Þ2 c

ð20Þ

ð29Þ

where g is the first-order derivative of G. The Lagrange multiplier kc is calculated by left multiplying Eq. (19) with wTc , i.e.,

Since the Jacobian matrix becomes diagonal, it is easily inverted. Then the orthogonal vector wc is calculated iteratively as below:

rLkc ¼ wTc wc - 1 ¼ 0;

 kc ¼ - q1 EfGðwTc zc;n Þg - EfGðv Þg  wTc  Efzc;n gðwTc zc;n Þg PNc   ðzc;n - z c Þðyc;n yc Þ ; - q2 covðwTc zc;n ; yc;n Þ  wTc  n¼1 Nc - 1 PNc

 zc

ð21Þ

PN c



where ¼ n¼1 zc;n =Nc and yc ¼ n¼1 yc;n =N c are the means of the whitening score vectors and the CEI of performance grade c, respectively. Moreover, substituting Eq. (21) into Eq. (19), rLwc can be rewritten as follows:

 

   Efzc;n gðwTc zc;n Þg - wTc  Efzc;n gðwTc zc;n Þg  wc

PNc   ðzc;n - z c Þðyc;n yc Þ n¼1 þ q2 2covðwTc zc;n ; yc;n Þ  Nc - 1 PNc n¼1





ðzc;n - z c Þðyc;n yc Þ Nc - 1

 wc

 :

Seeing from Eq. (22), one can obtain the gradient vectors and the scaling factors of each objective,

rLwc ;1



¼ 2 EfGðwTc zc;n Þg - EfGðv Þg   Efzc;n gðwTc zc;n Þg - wTc  Efzc;n gðwTc zc;n Þg  wc ; PNc

sf 2 ¼

1 k rLwc ;1 k2 1 k rLwc ;2 k2

;

:

n¼1 ðzc;n



ð23Þ



- zc Þðyc;n  yc Þ Nc - 1 !  PNc  ðzc;n - zc Þðyc;n  yc Þ  wTc  n¼1  wc ; Nc - 1

rLwc ;2 ¼ 2covðwTc zc;n ; yc;n Þ 

sf 1 ¼

wc;d

ð24Þ

ð25Þ

ð26Þ

According to the Kuhn–Tucker conditions (Luenberger, 1969), the optima of Jðwc Þ under the constraint wTc wc ¼ 1 will be obtained at the points of Eq. (19). Here Newton’s method is applied to solve this equation. The Jacobian matrix J F ðwc Þ of the left-hand side of Eq. (19) is expressed as

ð30Þ

wc;d -

d1 X ðwTc;k wc;k Þwc;k ;

ð31Þ

k¼1

wc;d ð22Þ

 1 wc;d - J F ðwc;d Þ  rLwc ;

where wc;d represents the dth column vector of W c , d ¼ 1; 2; :::; Dc . Performing orthogonalization on wc;d so as to exclude the information contained in the solutions already found,



rLwc ¼ q1 2 EfGðwTc zc;n Þg - EfGðv Þg

wTc 

wc;d

wc;d : k wc;d k

ð32Þ

If the old and new values of wc;d tend to the same direction, i.e. their dot-product is almost equal to 1, the algorithm is convergent, and repeat the above procedure to calculate the next one. After obtaining W c ¼ ½wc;1 ; wc;2 ; :::; wc;Dc , we can calculate sn;c by

sn;c ¼ W Tc zn;c ¼ U c kðxn;c Þ; n ¼ 1; 2; :::; Nc ;

ð33Þ

pffiffiffiffiffiffi T where U c ¼ Nc W Tc K1 Dc Ac . As can be seen from Eq. (17), J 2 ðwc Þ is a quadric form, so it is a convex function and its addition does not change the properties of the solution of the optimization problem as in Eq. (15). In other words, whether or not Jðwc Þ can obtain the global optima depends on the solution of J 1 ðwc Þ. In Ref. (Hyvärinen, 1999), Hyvärinen proved that the original FastICA algorithm does converge to the right extrema for J 1 ðwc Þ and the convergence is generally local. Furthermore, in the special case of GðuÞ ¼ u4 , he proved that the convergence is globally. Therefore, we can conclude that the convergence of the solution for Jðwc Þ is local; moreover, if GðuÞ ¼ u4 , it is global convergence. Remark 1. For multi-objective optimization problems, due to the constraints of each objective, the improvement of one objective performance is often at the cost of losing the performance of other objectives (Tsionas, 2019). Therefore, there cannot be an optimal solution for all objectives, and the solution is usually a set of noninferior solutions, i.e., Pareto solution set. It is difficult to determine which solution is more desirable if there is no more information about the problem. Thus, the most important task for multiobjective optimization problems is to find as many Pareto optimal

6

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

solutions as possible, and the following two tasks are mainly accomplished: (i) find a set of solutions as close as possible to Pareto front to ensure the accuracy of solutions; (ii) find a set of solutions as different as possible, that is to make the distribution of solutions as uniform as possible. However, how to obtain the Pareto optimal solution set is closely related to many factors, such as the nonlinear characteristics of multi-objective optimization problems, the number of decision variables, the coupling degree of constraints and the number of constraints, et al.

Remark 2. Multi-objective optimization algorithms can be summed up into two categories: intelligent optimization algorithm and traditional optimization algorithm. Intelligent optimization algorithms include evolutionary algorithms (EA) and particle swarm optimization (PSO), etc. Traditional optimization algorithms include weighting method, constraint method and linear programming method, etc, and these algorithms essentially transform the multi-objective function into a singleobjective function, and then solving the multi-objective function by using the single-objective optimization method (Zadeh, 1963). In this study, the traditional weighting method is used. Since it is impossible for each objective to achieve the optimal at the same time, each objective must have its own weight (Kim and Weck, 2005). However, it is difficult to determine the weights, and how to allocate such weights has become a research hotspot. It is worth noting that the ICs that have great correlations with the CEI usually have strong ability to distinguish between different operating performances, so they are named as performancerelevant ICs (PR-ICs). On the contrary, ICs that have little correlations with the CEI are difficult to distinguish the operating performances. They are called performance-irrelevant ICs (PI-ICs). In order to establish the assessment models to distinguish the process operating performances, only the PR-ICs are used and PI-ICs are actually interferences. The existence of PI-ICs greatly affects the sensitivity of the assessment models to the change of the performance-relevant variation information and the robustness to the disturbance of the performance-irrelevant variation information. Hence, in operating performance assessment, retaining too many ICs not only increases the computational load but also reduces the accuracy of the assessment model. To exclude those PI-ICs, one has to rank the order and determine the number of ICs. Over the past decade, some criteria have been proposed for this issue from the perspective of process monitoring (Lee et al., 2006; Liu et al., 2008). Recently, Zhang et al. (2018) develop a qualityrelevance contribution index and max-dependency to prioritize ICs for regression analysis. The pioneer work gives us a lot of inspiration and theoretical bases for our work. In this study, two new criteria are separately designed to determine the order and number of ICs for assessment modeling. First, we define the contribution of IC based on the correlation between each IC and the CEI, i.e.,

     corðsc;d ; yc Þ   ; d ¼ 1; 2; :::; Dc ; gd ¼ PD    c d¼1 corðs c;d ; y c Þ 

  where corðsc;d ; yc Þ ¼ covðsc;d ; yc Þ=k sc;d kk yc k   between sc;d and yc , and sc;d ¼ Z Tc wc;d 2 RNc 1 is

ð34Þ

is the correlation

the dth IC. Sort the ICs in descending order according to the values of gd , and represent the reordered ICs and the corresponding contributions as    sc;ð1Þ ; sc;ð2Þ ; :::; sc;ðDc Þ

and gð1Þ ; gð2Þ ; :::; gðDc Þ , respectively. Based on the reordered ICs, the cumulative percent contribution of ICs is calculated as follows:

Cgc;K ¼ 100

K X

gðdÞ =

d¼1

Dc X

!

gðdÞ %:

ð35Þ

d¼1

It is a measure of the cumulative percentage of the contribution captured by the first K ICs. Set the threshold of CgK as lð0 < l < 1Þ, and the number of PR-ICs can be determined by

n o   K c ¼ arg min Cgc;K Cgc;K P l; K ¼ 1; 2; :::; Dc :

ð36Þ

K

In theory, both the assessment models of performance grades and conversion processes should be built. However, since the conversion process usually shows time-varying and dynamic characteristic, it is hard to be represented directly. Fortunately, the process characteristic of the conversion process is closely associated with those of two adjacent performance grades and can be expressed indirectly by them. Hence, only the assessment models for each performance grade rather than the conversion processes are needed. In addition, under the assumption that normal historical data already cover all possible operating performances, the training data of different performance grades are readily available. Combining the process knowledge with expert experiences, the normal historical data can be further divided into several performance grades, and then it is easy to obtain the training data of different performance grades. With the proposed PR-KICA algorithm, the assessment models of each performance grade are developed. The independent component matrix S c and the demixing matrices W c , c ¼ 1; 2; :::; C are used as the model parameters to provide reference standards for online assessment.

3.2. Online operating performance assessment With the recognition that the same operating performance contains similar performance-relevant variation information, we can evaluate the current process operating performance based on the similarities between the performance-relevant variation information of the online data and that of each performance grade. Here the performance-relevant variation information is essentially the PR-ICs estimated by PR-KICA. Then an assessment index on the basis of the similarities between PR-ICs is constructed for online assessment. From the points of engineering and practice, a single sample is not enough to characterize the process operating performance and is susceptible to disturbance. Thus a data window with width H is introduced as the basic analysis unit. It serves as smoothing filter, and the width is dependent on the actual production situation. If the process runs smoothly with few outliers, a short window is appropriate. If the process is susceptible to interference, it is better to choose a long window to reduce the influence of disturbance and improve the reliability of the assessment result. Denote the online data window as X new ¼ ½xnew;tHþ1 ; xnew;tHþ2 ; :::; xnew;t  2 RJH , where xnew;t represents the sample collected at the t sampling instant. To measure the similarities of the performance-relevant variation information, the PR-ICs of the online data must be calculated with respect to each performance grade, i.e.,

snew;c ¼ W Tc znew;c ¼

pffiffiffiffiffiffi T 1 T   Nc W c KDc Ac kðxnew;c Þ ¼ U c kðxnew;c Þ;

c ¼ 1; 2; :::; C; 

ð37Þ 

where xnew is the mean vector of X new ; kðxnew;c Þ ¼ h iT    is the kernel vector of xnew with kðxc;1 ; xnew Þ; . . . ; kðxc;Nc ; xnew Þ 

respect to performance grade c. It should be noted that kðxnew;c Þ has been mean centered.

7

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

Furthermore, the similarities between snew;c and those of each performance grade is measured by squared Euclidean distance as follows: 

2

dnew;c ¼ k snew;c  sc k ;

ð38Þ

  P c P c where sc ¼ Nn¼1 sn;c =N c ¼ Nn¼1 U c kðxn;c Þ=Nc is the mean of PR-ICs of performance grade c. Since kernel vectors kðxn;c Þ; n ¼ 1; 2; :::; Nc 

have been mean centralized in offline modeling, sc is a zero vector. Then the squared Euclidean distance can be rewritten as follows:

dnew;c ¼ k snew;c k2 :

ð39Þ

To facilitate the assessment, an assessment index cnew;c is constructed as below:

cnew;c ¼

8 < :

1=

1=dnew;c P C c¼1

dnew;c

;

ifdnew;c –0

cnew;c ¼ 1 and cnew;q ¼ 0ðq–cÞ; if dnew;c ¼ 0

ð40Þ

where cnew;c ; c ¼ 1; 2; :::; C satisfies the conditions of 0 6 cnew;c 6 1 P and Cc¼1 cnew;c ¼ 1. The closer cnew;c is to 1, the more similar the PR-ICs between the online data and the training data of performance grade c are, which indicates that the current process is operating on performance grade c. On the contrary, if cnew;c is close to 0, it means that the current process is almost impossible to be performance grade c. Being similar to the assessment rules proposed in Refs. (Liu et al., 2013, 2016), we introduce an assessment index threshold h for making a strict distinction between the performance grade and the conversion process. The values range of h is from 0 to 1 and affects the assessment result to some extent. If it is too high, the optimal operating performance may be evaluated as nonoptimal for the effect of interferences. Then the robustness of the assessment method is reduced. However, if it is too low, the nonoptimal operating performance may be evaluated as optimal, which increases the rate of error evaluation and reduces the accuracy of the assessment result. In practical application, h is usually determined using the historical data by trial and error through crossvalidation so that the numbers of online false and missing alarms are minimized. Then the process operating performance can be evaluated by following rules: Case 1: if the maximum value of the assessment indices relative to all performance grades is greater than the threshold h, n o i.e.,cnew;c ¼ max cnew;c jc ¼ 1; 2; :::; C P h, it indicates that the current operating performance is consistent with that of performance 

grade c . Case 2: if case 1 is not set up but the condition cnew;c ðt  l þ 1Þ 6 ::: 6 cnew;c ðtÞ is satisfied, it means that the current process is gradually converting from one performance grade to performance grade c. Here l is a positive integer, and cnew;c ðtÞ is the assessment index value of performance grade c at the tth sampling. Since there may be more than one performance grade meeting the above condition, the target performance grade can be determined by

 c ¼ arg maxfcnew;c ðtÞcnew;c ðt  l þ 1Þ 6 ::: 6 cnew;c ðtÞ; c ¼ 1; 2; :::; Cg;



c

ð41Þ 

where c is the argument that cnew;c ðtÞ has the maximum value under the condition of cnew;c ðt  l þ 1Þ 6 ::: 6 cnew;c ðtÞ. Case 3: if both case 1 and 2 are not satisfied, it can be attributed to the influence of noise or uncertainties, and the assessment result remains consistent with that of the previous sampling instant.

Thereby, the online assessment results include not only the performance grade but also the conversion process, which improves the practicability of the assessment strategy. 3.3. Nonoptimal cause identification In fact, except performance grade optimal, all the other performance grades and the conversion processes belong to the nonoptimal operating performance. They are unexpected in practice. However, once it appears, we need to identify the cause quickly and make a response, which is very similar to the issue of process fault diagnosis (Qin et al., 2001; Du et al., 2013;Peng et al., 2013; Mnassri et al., 2015; Li et al., 2018). Contribution plot is a popular approach in fault diagnosis, and in linear cases one can calculate the variable contributions through decomposing the monitoring statistics into the weighted sum of different process variables. Nevertheless, it is not so simple and easy in nonlinear cases. Due to the introduction of the nonlinear kernel function, the variable contributions cannot be obtained by decomposing the assessment index directly. To overcome this challenge, a novel variable contributions based nonoptimal cause identification strategy is developed in this study. We draw lessons from the idea of variable selection that is to extract the relevant variables making significant contributions to the deviation of the predefined criterion, and define the variable contributions in a sense of meaning. Here, the search for cause variable is to select the variables having relatively large contributions on the deviation of the assessment index with respect to performance grade optimal. Denote performance grade optimal as c . It is noted that the value of assessment index cnew;c depends on the value of the distance dnew;c . Hence, we can calculate the variable contributions to the distance dnew;c instead of the assessment index cnew;c . Here the distance dnew;c can be rewritten as follows:

  

T  dnew;c ¼ jjsnew;c jj2 ¼ tr U c kðxnew;c Þk xnew;c U Tc :

ð42Þ

Then the contribution of the jth variable is defined as follows:

Contrnew;c ;j

where

    @d    new;c ¼   Dxnew;c ;j   @xnew;c ;j   2 3

T       Þk xnew;c kðx new;c    6 T 7 ¼ tr 4U c  U D x 5 new;c ;j ;  c   @x new;c ;j  

h iT    kðxnew;c Þ ¼ kðxc ;1 ; xnew Þ; . . . ; kðxc ;Nc ; xnew Þ ,

    xnew;c ;j xc ;j ; xnew;c ;j and xc ;j  respectively; xc is the mean



ð43Þ



Dxnew;c ;j ¼ 

are the jth variables of xnew;c and xc ,

vector of training data of performance grade optimal. Obviously, the variable contribution is composed of 



two items, i.e., @dnew;c =@xnew;c ;j and Dxnew;c ;j . One is the gradient of  xnew;c ;j ,

and it essentially represents the with respect to variable change rate of dnew;c along the direction of the jth variable. The 

other is the deviation item, Dxnew;c ;j , between the current sampling and the training data. Their product reflects the influence of a variable on the deviation of assessment index with respect to performance grade optimal. Therefore, it is reasonable to define the variable contributions in this form. Since the greater the contribution is, the larger the deviation is, and the variables with relatively large contributions are determined as the possible cause variables responsible for the nonoptimality. To calculate the variable contribution Contrnew;c ;j analytically, 

we should first calculate the gradient of kðxc ;n ; xnew;c Þ with respect  toxnew;c ;j

as below:

8

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167 

@kðxc ;n ; xnew;c Þ 

@xnew;c ;j

¼

1

r





ðxc ;n;j  xnew;c ;j Þkðxc ;n ; xnew;c Þ:

ð44Þ

Then the gradient of the nth row and mth column element of 

T  to xnew;c ;j is kðxnew;c Þk xnew;c 





@kðxc ;n ;xnew;c Þkðxc ;m ;xnew;c Þ 

@xnew;c ;j 1







¼ ðxc ;n;j þ xc ;m;j  2xnew;c ;j Þkðxc ;n ;xnew;c Þkðxc m ; xnew;c Þ:

r

ð45Þ

Substituting Eq. (45) into Eq. (43), one can obtain the precise form of the variable contributions. Owing to the fact that the feature mapping may affect variables unevenly, the scales of variable contributions may be inconsistent even in the optimal operating performance, and it is meaningless to compare their absolute magnitudes. For this issue, a general practice is to normalize the raw contributions so that all variables give the same contributions under optimal operating performance statistically. Thus, the variable contributions are scaled as follows: 

Contrnew;c ;j ¼

Contrnew;c ;j 

Contrc ;j

; j ¼ 1; 2; :::; J;

Fig. 1. Flow diagram of the PR-KICA based operating performance assessment method.

ð46Þ

9

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

where 

Contrc ;j ¼

 

    T    PNc   kðxn;c Þk xn;c =@xn;c ;j U Tc  Dxn;c ;j =Nc is n¼1 tr U c

the mean of contributions of the jth variable of performance grade optimal. Then the variables with relevant large contributions are identified as the latent cause variables. 3.4. Implementation procedures In summary, the step-by-step procedures of the proposed PRKICA based assessment method are given as follows and the flow diagram is shown in Fig. 1. 3.4.1. Offline modeling procedure (1) Acquire the training data of different performance grades respectively and normalize the data using the mean and standard deviation of each process variable. (2) Compute the Gram matrix of the training data and perform mean centering on the Gram matrix. (3) Implement eigenvalue decomposition on the Gram matrix and calculate the whitening score matrix. (4) Solve the demixing matrix using PR-KICA and establish the assessment models for each performance grade. 3.4.2. Online application procedure Fig. 2. Schematic diagram of gold hydrometallurgy process.

(1) Construct the online data window and normalize the online data with the means and standard deviations of different performance grades, respectively. (2) Compute the kernel vector of the online data and perform mean centering on the kernel vector. (3) Calculate PR-ICs of the online data with respect to each performance grade. (4) Calculate the Euclidean distances and assessment indices between the PR-ICs of online data and those of each performance grade. (5) Carry out the online assessment based on the assessment index. (6) Identify the nonoptimal cause variable based on variable contribution. 4. Case study: gold hydrometallurgy process Gold is a kind of precious metal with great value and widely used in modern production and daily life. Nevertheless, with the decrease of mineral resources, people have to make more effort for extracting gold from the low-grade ores. Gold hydrometallurgy process is an effective technique to deal with the low-grade ores. In this process, gold ore is dissolved in the solution and precipitated in the form of a new solid phase for metal separation, enrichment and extraction (Jackson, 1986; Abbruzzese et al., 1995; Leão and Ciminelli, 2000; Kumar et al., 2001). Nonoptimal operating performance will have a serious impact on the recovery rate of gold and the comprehensive economic benefits. Thus, it is very important to implement the online operating performance assessment for gold hydrometallurgy process. In this study, it is introduced as the background to verify the effectiveness of the proposed approach. 4.1. Process description Taking a production line of gold hydrometallurgy process in Shandong, China as an example, it consists of five units: first cyanide leaching, first washing, second cyanide leaching, second

washing and cementation. Fig. 2 is the schematic diagram of this process. First, flotation pulp is sent into the first cyanide leaching unit, meanwhile, both sodium cyanide (NaCN) and oxygen are continuously added into the leaching tanks, where NaCN is the reagent and Table 1 Process variables of gold hydrometallurgy process. No.

Operating units

Process variables

Units

1

First cyanide leaching

NaCN flow 1

2

NaCN flow 2

3 4 5 6 7 8

Air flow 1 Air flow 2 Dissolved oxygen concentration CN concentration 1 CN concentration 2 Buffer level

ml/ min ml/ min m3/h m3/h mg/l mg/l mg/l m

9 10 11

First washing

Feed pressure Extrusion pressure Hydraulic pressure

MPa MPa MPa

12

Second cyanide Leaching

NaCN flow 1

13

NaCN flow 2

14 15 16 17 18 19

Air flow 1 Air flow 2 Dissolved oxygen concentration CN concentration 1 CN concentration 2 Buffer level

ml/ min ml/ min m3/h m3/h mg/l mg/l mg/l m

20 21 22

Second washing

Feed pressure Extrusion pressure Hydraulic pressure

MPa MPa MPa

23

Cementation

Vacuum degree of deaeration tower [Au(CN)2]concentration Zinc powder adding amount Hydraulic pressure

kPa

24 25 26

mg/l Kg/h MPa

10

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

(a) Nonlinearity coefficient of grade poor

(b) Nonlinearity coefficient of grade general 1

10 15

0.8

5

0.6

10

0.4

20

0.8 Variables

Variables

5

0.6 15 0.4 20

0.2

0.2

25

25 5

10 15 Variables

20

25

5

(c) Nonlinearity coefficient of grade suboptimal

10 15 Variables

20

25

(d) Nonlinearity coefficient of grade optimal 1 5

0.8

10

0.6

15

0.4

20

0.2

25

0.8 Variables

Variables

5

10 0.6 15 0.4 20 0.2 25

5

10 15 Variables

20

25

5

10 15 Variables

20

25

Fig. 3. Nonlinearity coefficients between variables of performance grades (a) poor; (b) general; (c) suboptimal; (d) optimal.

oxygen provides the oxidation-reduction potential to promote the reaction. After that, the pulp flows into the first washing unit. It is to separate the high-concentration [Au(CN)2] solution from the pulp through vertical filter presses. High solid-liquid separation efficiency is the goal pursued in this unit, and the feed pressure, extrusion pressure and hydraulic pressure are the key parameters. Separated high-concentration [Au(CN)2] solution flows into the cementation unit for recovering gold from liquid phase to solid phase. In this unit, zinc powder is the deoxidizer, and the frame filter press is used as the carrier to provide the place for reaction. The second cyanide leaching and washing units are the same to the first ones correspondingly, and the duplicate units are to extract gold as much as possible and avoid gold loss. The process variables of gold hydrometallurgy process are listed in Table 1. According to the process mechanism, the buffer levels (variable 8 and 19) do not affect the process operating performance and only used to simulate the interferences of performanceirrelevant variant information. In addition, the CEI used in this process is the comprehensive economic benefit under the consideration of many factors, including energy, material, labour cost, environmental protection, etc. Based on the expert experience

and historical data, this process is divided into four performance grades, i.e. poor, general, suboptimal, and optimal. The nonoptimal causes of performance grade poor, general, and suboptimal are the lower setting of the zinc powder adding amount (variable 25) than its optimal operating point. For each performance grade, we collect 700 samples from the historical data and used as the training data to establish the assessment model. The sampling interval is 1 min. In order to evaluate the applicability of the proposed method before offline modeling, we first test the characteristics of the training data by the method proposed in Ref. (Zhang et al., 2016). This data characteristic test method consists of two parts: tests of dependencies among variables and the distribution characteristic of data. Fig. 3 shows the nonlinear coefficients between variables of each performance grade. It can be seen that the nonlinear dependencies are similar across different performance grades. One can determine which variables have strong nonlinear dependencies and which variables have weak dependencies. For example, there are strong nonlinear dependencies for the first seven variables, while those between variable 8 and other 25 variables are weak. This result is in line with the actual situation. It is because that

Table 2 Test results of data characteristics in each performance grade. s F

ð< 0:15Þ

Performance grade



Poor General Suboptimal Optimal

0.0827 0.0651 0.0554 0.0711



jajðP F 5% = 3.4231)

jb  1jðP 0:2Þ

Distribution characteristic

rðP 0:4Þ

Correlation

17.637 15.516 13.065 15.558

0.4766 0.4125 0.3070 0.3964

Non-Gaussian Non-Gaussian Non-Gaussian Non-Gaussian

0.4397 0.4427 0.4442 0.4450

Nonlinear Nonlinear Nonlinear Nonlinear

11

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

(a) Poor

(b) Grneral

70

70 Point

Point 60

Fitting line

Fractile of F distribution

Fractile of F distribution

60 50 40 30 20 10

0

20

40

60

80

50 40 30 20 10

100

Fitting line

0

Mahalanobis distance D 2

20

40

(c) Suboptimal

100

70 Point

60

Point 60

Fitting line

Fractile of F distribution

Fractile of F distribution

80

(d) Optimal

70

50 40 30 20 10

60

Mahalanobis distance D2

0

20

40

60

80

Mahalanobis distance D2

Fitting line

50 40 30 20 10

0

20

40

60

80

Mahalanobis distance D2

Fig. 4. Data distribution test of performance grades (a) poor; (b) general; (c) suboptimal; (d) optimal.

there are three groups of significant nonlinear dependencies within the first cyanide leaching unit, i.e., NaCN flow and CN concentration, air flow and dissolved oxygen concentration, dissolved oxygen concentration and CN concentration; however, the buffer level has little dependency with others. The second cyanide leaching unit has a similar situation. In addition, the variables within washing and cementation units also show strong nonlinear dependencies, respectively. Besides the variable dependencies within each unit, some variables show nonlinear dependencies across units, such as variables 1, 2, 6, 7, 12, 13, 17, 18, 24, and 25. The population nonlinearity coefficients of different performance grades are listed in Table 2, and the threshold Tr is set as 0.4. Since all the population nonlinearity coefficients are larger than the threshold of 0.4, the process data are nonlinear. For nonlinear processes, the test of data distribution characteristic is implemented on the mapped data in the high-dimensional 





feature space. r and b are separately set as F 5% and 0.2, and the test results are illustrated in Fig. 4 and Table 2, respectively. It is obvious that the points are not scattered on a line passing through the original point with slope equal to 1, thus the process data are non-Gaussian. In summary, the proposed assessment method is applicable to the training data of gold hydrometallurgy process. The parameters used in offline modeling are set as follows: the Gaussian kernel width r ¼ 100, the threshold of the cumulative percent contribution of ICs l ¼ 0:95, the preference factors b1 ¼ b2 ¼ 0:5. The assessment models of performance grade poor, general, suboptimal, and optimal are established by the PR-KICA

method, respectively. The numbers of PR-ICs retained in each performance grade are 6, 7, 6, and 6 correspondingly.

4.2. Analysis and discussion To verify the efficiency of the proposed PR-KICA based operating performance assessment method, we consider two production scenarios here. The first scenario is that, when the process operates on performance grade optimal, the buffer level (variable 8) gradually increases within a safe range. Since the function of the buffer is to buffer the pulp and prevent the risk of groove accident, it should be considered as a fluctuation of the performance-irrelevant variation information and has no influence on the process operating performance. In the second scenario, NaCN flows (variables 1 and 2) gradually decrease and deviate from their optimal set points due to human error under the optimal operating performance. Since the NaCN flow directly influences the concentration of the reactant and the leaching rate of gold, it is actually a change of the performance-relevant variation information, resulting in the nonoptimal operating performance. Each test dataset include 2205 samples, and the scenario is introduced at the 1901th sampling. Additionally, the order of performance grades in both scenarios is poor, general, suboptimal, and optimal. The parameters used in online assessment are set as follows: the width of data windowH ¼ 20, the threshold of assessment index h ¼ 0:85, the continuous incremental number l ¼ 5.

12

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

Assessment index

(a) Poor 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400

1600 1800 2000 2200

Samplings Assessment index

(b) General 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400

1600 1800 2000 2200

Samplings Assessment index

(c) Suboptimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400

1600 1800 2000 2200

Samplings Assessment index

(d) Optimal 1

Assessment index Threshold

0.5

At the 1996 th sampling 0 200

400

600

800

1000 1200 1400

1600 1800 2000 2200

Samplings Fig. 5. PLS based assessment results for the first scenario.

For comparison, PLS and KICA are also applied to both scenarios. Figs. 5–7 are the online assessment results for the first scenario based on PLS, KICA, and PR-KICA, respectively. From Fig. 5, we can see that the PLS method makes many wrong evaluations before the 1901th sampling. After introducing the interference of buffer level, it also responds to the interference and gives wrong assessment results, i.e., the assessment result converts from optimal to nonoptimal at the 1996th sampling. On the one hand, it is because that the PLS method lacks the ability to deal with nonlinear and non-Gaussian features, thus it cannot accurately extract the major process variation information from the training data; on the other hand, as Zhou et al. (2010) have proved that there is also performance-irrelevant variation information in the PLS score space, it has no contribution for inferring the CEI and reduces the anti-interference ability of the assessment model. As can be seen from Fig. 6, KICA also gives a poor evaluation performance. With respect to the interference of buffer level, the KICA method makes an unwarranted response, and the assessment result converts from

optimal to nonoptimal at the 1983th sampling, which is not expected in practice. Tracing it to its cause, in process operating performance assessment, only those PR-ICs really play a role in distinguishing different operating performances. Since KICA cannot classify its ICs into PR-ICs and PI-ICs, both types of information are included in the KICA assessment models. Therefore, when only the performance-irrelevant variation information within the online data deviations from previous operating performance, KICA will respond and further make a wrong judgment. In comparison, PRKICA outperforms PLS and KICA as shown in Fig. 7. It indicates that, for the interference of performance-irrelevant variation information, the PR-KICA method has stronger robustness than PLS and KICA methods. Detailed comparisons of the online assessment results by the three methods are shown in Table 3, where the assessment results are compared against the actual situation. As shown in Table 3, three methods having in common is that the online assessment results have time delay compared with the actual situation, which may be caused by the data window used

Assessment index

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

13

(a) Poor 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(b) General 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(c) Suboptimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(d) Optimal 1

Assessment index Threshold

0.5

At the 1983 th sampling 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Fig. 6. KICA based assessment results for the first scenario.

in online assessment. Nevertheless, PR-KICA captures the change of operating performance earlier than KICA and PLS and is more accurate. Defining the accuracy rate of the assessment results as the percentage between the numbers of correct evaluated samplings and the total test samplings, it is up to 98.04% for PR-KICA, while they are 83.67% and 84.85% for PLS and KICA, respectively. To investigate the reason for the nonoptimal operating performance, the variable contributions are calculated separately under the frameworks of PLS, KICA and PR-KICA. In order to avoid the smearing effect on the unrelated variables, the variable contributions are calculated at the first sampling of performance grade poor and shown in Fig. 8, where the contributions of [Au(CN)2]concentration (variable 24) and zinc powder adding amount (variable 25) are all greater than those of other variables. Known from the process mechanism, zinc powder adding amount is an operation variable and affects [Au(CN)2]concentration, so it is the true cause

responsible for the nonoptimality. The cause identification result is consistent with the actual situation. In addition, variable contributions at the 1996th and 1983th samplings are separately calculated based on PLS and KICA and exhibited in Fig. 9. It shows that the buffer level (variable 8) is uniquely identified as the responsible cause variable with the largest contribution. Those result may mislead persons who are not familiar with the actual production process and also illustrates the shortage of the PLS and KICA methods from another angle. In the second scenario, the online assessment results based on the PLS method is shown in Fig. 10. Although PLS can identifies the change of the process operating performance at the 2024th sampling, it is with a big time delay compared with the actual situation. This is probably due to the inadequacy of PLS in processing nonlinear and non-Gaussian information and the degradation sensitivity for the existence of the performance-irrelevant variation

14

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

Assessment index

(a) Poor 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(b) General 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(c) Suboptimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Assessment index

(d) Optimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Fig. 7. PR-KICA based assessment results for the first scenario.

Table 3 Comparisons between the actual situation and online assessment results for the first scenario. Operating performance

Actual situation (sampling)

PLS (sampling)

KICA (sampling)

PR-KICA (sampling)

Poor Conversion General Conversion Suboptimal Conversion Optimal Conversion Accuracy rate (%)

1–460 461–579 580–1081 1082–1203 1204–1654 1655–1789 1790–2205 – –

1–474 475–638 639–1101 1102–1235 1236–1672 1673–1797 1798–1995 1996–2205 83.67

1–479 480–610 611–1092 1093–1218 1219–1680 1681–1798 1799–1982 1983–2205 84.85

1–463 464–607 608–1083 1084–1203 1204–1661 1662–1792 1793–2205 – 98.04

information within the PLS assessment model. As shown in Fig. 11, the KICA method evaluates the change of the process operating performance at the 2004th sampling, which is also later than the

actual situation. This may be due to the fact that the existence of performance-irrelevant variation information reduces the sensitivity of KICA to the change of the operating performance. In compar-

15

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

(a) PLS

(b) KICA

7

(c) PR-KICA

10

14

6

12

4 3 2

Vrairable constructions

5

Vrairable constructions

Vrairable constructions

8

6

4

10 8 6 4

2 1 0

2

0

10 20 Variable No.

0

0

10 20 Variable No.

0

0

10 20 Variable No.

Fig. 8. Cause identification results for the first scenario at the 1st sampling (a) PLS method; (b) KICA method; (c) PR-KICA method.

(a) PLS

(b) KICA

6

7

At the 1983th sampling

At the 1996th sampling 6 Vrairable constructions

Vrairable constructions

5

4

3

2

1

0

5 4 3 2 1

0

10 20 Variable No.

0

0

10 20 Variable No.

Fig. 9. Cause identification results for the first scenario (a) PLS method; (b) KICA method.

ison, the PR-KICA method captures the deviation at the 1944th sampling and performs significantly better than PLS and KICA as shown in Fig. 12. Since PR-KICA only retains PR-ICs in the assessment models, it focuses more on the change of the performancerelevant variation information. Thus, the PR-KICA method has higher sensitivity than other two methods with respect to the change of performance-relevant variation information and can capture this kind of changes much earlier. Furthermore, comparisons of the online assessment results against the actual situation are tabulated in Table 4, where the accuracy rate of the

PR-KICA method is significantly higher than those of PLS and KICA methods. The nonoptimal cause identification results based on three methods are shown in Fig. 13. Seen from that, the contributions of NaCN flows (variables 1 and 2) and CN concentration (variables 6 and 7) are greater than those of other variables. Since the NaCN flow is an operation variable and affects the CN concentration, it is the true cause variable responsible for the nonoptimality, which is consistent with the actual situation.

16

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

Assessment index

(a) Poor 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Assessment index

(b) General 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Assessment index

(c) Suboptimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Assessment index

(d) Optimal 1

Assessment index Threshold

0.5 At the 2024 th sampling 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Fig. 10. PLS based assessment results for the second scenario.

Through the above comparative analyses, we can conclude that, for the subject of process operating performance assessment, the proposed PR-KICA method has stronger antiinterference ability for the interferences of performanceirrelevant variation information and higher sensitivity for the changes in performance-relevant variation information than PLS and KICA methods. Generally, when generalizing the proposed assessment method to other processes, the first step is to test whether the process data are nonlinearity and non-Gaussian, testing the appropriateness of the proposed method. If applicable, standardize the training data of each performance grade to remove the influence of noise and

dimension, preparing for offline modeling. Then use the PR-KICA method to establish the assessment models for each performance grade, respectively. In online assessment, calculate the value of the assessment index in real time according to the sampling frequency. For the case of inconsistent sampling frequencies across multiple variables, appropriate measures need to be predefined in practice. When the process operating performance is optimal, only the online assessment result is output; if it is nonoptimal, both the online assessment result and the cause identification result are output, and the operators can determine the adjustment scheme according to the cause identification result and their own production experience.

Assessment index

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

17

(a) Poor 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(b) General 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(c) Suboptimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(d) Optimal 1

Assessment index Threshold

0.5

At the 2004 th sampling

0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Fig. 11. KICA based assessment results for the second scenario.

5. Conclusions In this study, we proposed a novel operating performance assessment method based on PR-KICA for nonlinear and non-Gaussian industrial processes. The objectives of PR-KICA are simultaneously to maximize the non-Gaussianity of ICs as well as the correlations between them and the CEI. Two new criteria are designed separately to rank the order and determine the number of ICs. It helps to remove the ICs that are weakly related to the CEI and select some dominated ICs in improving the evaluation performance of the assessment method. In online assessment, the assessment results include not only performance grades but also conversion processes, which is more in line with the actual situation. To find the nonopti-

mal cause, the idea of variable selection is used to construct the variable contributions to overcome the nonlinearity of the assessment index. Finally, the effectiveness of the PR-KICA method is verified by the process data from gold hydrometallurgy process. Due to emphasis on performance-relevant characteristics, PR-KICA has strong robustness to the interference of the performance-irrelevant variation information and high sensitivity with respect to the change of the performance-relevant variation information compared to other assessment methods. As a future research, we need to consider the impact of time-varying and dynamic on the assessment model, as well as the intelligent decision-making and control (Li, et al., 2018) for the nonoptimal operating performance.

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

Assessment index

18

(a) Poor 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Assessment index

Samplings

(b) General 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Assessment index

(c) Suboptimal 1

Assessment index Threshold

0.5 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Assessment index

(d) Optimal 1

Assessment index Threshold

0.5 At the 1944 th sampling 0 200

400

600

800

1000 1200 1400 1600 1800 2000 2200

Samplings Fig. 12. PR-KICA based assessment results for the second scenario.

Table 4 Comparisons between the actual situation and online assessment results for the second scenario. Operating performance

Actual situation (sampling)

PLS (sampling)

KICA (sampling)

PR-KICA (sampling)

Poor Conversion General Conversion Suboptimal Conversion Optimal Conversion Accuracy rate (%)

1–454 455–576 577–1076 1077–1200 1201–1642 1643–1784 1785–1900 1901–2205 –

1–473 474–639 640–1088 1089–1236 1237–1671 1672–1793 1794–2023 2024–2205 86.8

1–468 469–600 601–1091 1092–1213 1214–1675 1676–1793 1794–2003 2004–2205 90.43

1–464 465–578 579–1078 1079–1209 1210–1655 1656–1788 1789–1943 1944–2205 96.23

19

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

(a) PLS

(b) KICA

8

(c) PR-KICA

10

12

At the 2024th sampling

At the 2004th sampling

At the 1944th sampling

7

10

5 4 3 2

Vrairable constructions

6

Vrairable constructions

Vrairable constructions

8

6

4

2

0

10 20 Variable No.

0

6 4

2

1 0

8

0

10 20 Variable No.

0

0

10 20 Variable No.

Fig. 13. Cause identification results for the second scenario (a) PLS method; (b) KICA method; (c) PR-KICA method.

Declaration of Competing Interest None.

Acknowledgments This work was supported by the National Natural Science Foundation of China [grant numbers 61703078, 61533007, 61873053, 61773105, and 61973057]; Natural Foundation Funding Schemes of Liaoning Province of China [grant number 2019-MS-115]; the Fundamental Research Funds for the Central Universities [grant number N180404009]; the Hong Kong Research Grant Council (RGC) on Synchronized Process Control of Hot-Runner Temperatures [grant number 16233316]; the Funds for Creative Research Groups of China [grant number 61621004]; and the Stat Key Laboratory of Synthetical Automation for Process Industries Fundamental Research Funds [grant number 2013ZCX02-04].

References Abbruzzese, C., Fornari, P., Massidda, R., Veglio, F., Ubaldini, S., 1995. Thiosulphate leaching for gold hydrometallurgy. Hydrometallurgy 39, 265–276. Chiang, L.H., Kotanchek, M.E., Kordon, A.K., 2004. Fault diagnosis based on Fisher discriminant analysis and support vector machines. Comput. Chem. Eng. 28, 1389–1401. Chiang, L.H., Russell, E.L., Braatz, R.D., 2000. Fault diagnosis in chemical processes using Fisher discriminant analysis, discriminant partial least squares, and principal component analysis. Chemom. Intell. Lab. Syst. 50, 243–252. Chu, F., Dai, W., Shen, J., Ma, X., Wang, F.L., 2018. Online complex nonlinear industrial process operating optimality assessment using modified robust total kernel partial M-regression. Chinese J. Chem. Eng. 26, 775–785. 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. Dunia, R., Qin, S.J., Edgar, T.F., McAvoy, T.J., 1996. Identification of faulty sensors using principal component analysis. AIChE J. 42, 2797–2812. Edgar, T.F., Himmelblau, D.M., Lasdon, L.S., 2001. Optimization of Chemical Processes. McGraw-Hill, New York. Haghani, A., Jeinsch, T., Ding, S.X., 2014. Quality-related fault detection in industrial multimode dynamic processes. IEEE Trans. Ind. Electron. 61, 6446–6453. Hampel, F.R., Ronchetti, E.M., Rousseuw, P.J., Stahel, W.A., 1986. Robust Statistics. Wiley, New York.

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. Hyvärinen, A., 1998. Independent componenet analysis by general nonlinear Hebbian-like learning rules. Signal Process. 64, 301–313. Hyvärinen, A., 1999. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural Netw. 10, 626–634. Hyvärinen, A., Oja, E., 2000. Independent component analysis: algorithms and applications. Neural Networks 13, 411–430. Jackson, E., 1986. Hydrometallurgical Extraction and Reclamation. Wiley, New York. Jackson, J.E., Wiley, J., 1991. A User’s Guide to Principal Components. Wiley Online Library. Kim, Y., Weck, O., 2005. Adaptive weighted-sum method for bi-objective optimization: pareto front generation. Struct. Multidisc. Optim. 29, 149–158. Kocsor, A., Csirik, J., 2001. Fast independent component analysis in kernel feature spaces. Theor. Pract. Inform. 2234, 271–281. Kourti, T., MacGregor, J.F., 1996. Multivariate SPC methods for process and product monitoring. J. Qual. Technol. 28, 409–428. Kumar, A., Haddad, R., Sastre, A., 2001. Integrated membrane process for gold recovery from hydrometallurgical solutions. AIChE J. 47, 328–340. Leão, V.A., Ciminelli, V.S., 2000. Application of ion exchange resins in gold hydrometallurgy: a tool for cyanide recycling. Solvent Extr. Ion Exch. 18, 567– 582. Lee, J.M., Qin, S.J., Lee, I.B., 2006. Fault detection and diagnosis based on modified independent component analysis. AIChE J. 52, 3501–3514. Lee, J.M., Qin, S.J., Lee, I.B., 2007. Fault detection of non-linear processes using kernel independent component analysis. Can. J. Chem. Eng. 85, 526–536. Li, Q.K., Dimirovski, G.M., Fu, J., Wang, J., 2018. Switching strategy in traking constant references for linear time-varying-delay systems with actuator failure. Int. J Control 90. https://doi.org/10.1080/00207179.2017.1415464. Liu, X., Xie, L., Kruger, U., Littler, T., Wang, S., 2008. Statistical-based monitoring of multivariate non-Gaussian systems. AIChE J. 54, 2379–2391. Liu, Y., Wang, F.L., Chang, Y.Q., 2013. Online fuzzy assessment of operating performance and cause identification of non-optimal grades for industrial processes. Ind. Eng. Chem. Res. 52, 18022–18030. Liu, Y., Chang, Y.Q., Wang, F.L., 2014. Online process operating performance assessment and nonoptimal cause identification for industrial processes. J. Process Control 24, 1548–1555. 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. Liu, Y., Wang, F.L., Chang, Y.Q., 2017. Operating optimality assessment and cause identification for nonlinear industrial processes. Chem. Eng. Res. Des. 117, 472– 487. Liu, Y., Wang, F.L., Chang, Y.Q., Ma, R.C., 2015. Operating optimality assessment and nonoptimal cause identification for non-Gaussian multimode processes with transitions. Chem. Eng. Sci. 137, 106–118. Luenberger, D.G., 1969. Optimization by Vector Space Methods. Wiley, New York. Mardia, K.V., 1977. Mahalanobis distances and angles. Mautivariate Anal., 495–511 Mejdell, T., Skogestad, S., 1991. Estimation of distillation compositions from multiple temperature measurements using partial-least-squares regression. Ind. Eng. Chem. Res. 30, 384–390.

20

Y. Liu et al. / Chemical Engineering Science 209 (2019) 115167

Mnassri, B., Adel, E.M.E., Ouladsine, M., 2015. Reconstruction-based contribution approaches for improved fault diagnosis using principal component analysis. J. Process Control 33, 60–76. Peng, K., Zhang, K., Li, G., Zhou, D., 2013. Contribution rate plot for nonlinear quality-related fault diagnosis with application to the hot strip mill process. Control Eng. Pract. 21, 360–369. Plovoso, M., Kosanovich, K., 1994. Applications of multivariate statistical methods to process monitoring and controller design. Int. J Control 59, 743–765. Qin, S.J., 2003. Statistical process monitoring: basics and beyond. J. Chemom. 17, 480–502. Qin, S.J., Valle, S., Piovoso, M.J., 2001. On unifying multiblock analysis with application to decentralized process monitoring. J. Chemom. 15, 715–742. Qin, S.J., Zheng, Y., 2013. Quality-relevant and process-relevant fault monitoring with concurrent projection to latent structures. AIChE J. 59, 496–504. Schölkopf, B., Smola, A., Müller, K.R., 1998. Nonlinear component analysis as a kernel eigenvalue problem. Neural Comput. 10, 1299–1319. Sedghi, S., Huang, B., 2017. Real-time assessment and diagnosis of process operating performance. Engineering 3, 214–219. Tsionas, G., 2019. Multi-objective optimization using statistical models. Eur. J. Oper. Res. 276, 364–378. Srour, M., Gomes, V.G., Altarawneh, I., 2008. Optimal operating strategies for emulsion terpolymerisation. Chem. Eng. Sci. 63, 4257–4268. Székely, G.J., Rizzo, M.L., 2013. Energy statistics: a class of statistics based on distances. J Stat. Plan. Infer. 143, 1249–1272. Vapnik, V.N., 1998. Statistical Learning Theory. John Wiley and Sons INC, New York. Wang, J., Zhong, B., Zhou, J., 2017a. Quality-relevant fault monitoring based on locality preserving partial least squares statistical models. Ind. Eng. Chem. Res. 56, 7009–7020.

Wang, J.Q., Ye, Z.F., Hu, Z.Z., Wu, X., Dimirovski, G.M., Yue, H., 2017b. Model-based acceleration control of turbofan engines with a Hammerstein-Wiener pepresentation. Int. J. Turbo Jet. Eng. 34, 141–148. Wang, Q., Shen, Y., Zhang, J., 2005. A nonlinear correlation measure for multivariable data set. Phys. D. 200, 287–295. Yang, M.S., 1993. A survey of fuzzy clustering. Math. Comput. Model. 8, 1–16. Ye, L., Liu, Y., Fei, Z., 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. Yin, S., Zhu, X., Kaynak, O., 2015. Improved PLS focused on key-performanceindicator-related fault diagnosis. IEEE Trans. Ind. Electron. 62, 1651–1658. Zadeh, L., 1963. Optimality and non-scalar-valued performance criteria. IEEE T. Automat. Contr. 8, 59–60. Zhang, S.M., Wang, F.L., Zhao, L.P., Wang, S., Chang, Y., 2016. A novel strategy of the data characteristics test for selecting a process monitoring method automatically. Ind. Eng. Chem. Res. 55, 1642–1654. Zhang, X., Kano, M., Li, Y., 2018. Quality-relevant independent component regression model for virtual sensing application. Comput. Chem. Eng. 115, 141–149. Zhang, Y., 2008. Fault detection and diagnosis of nonlinear processes using improved kernel independent component analysis (KICA) and support vector machine (SVM). Ind. Eng. Chem. Res. 47, 6961–6971. Zhou, D., Li, G., Qin, S.J., 2010. Total projection to latent structures for process monitoring. AIChE J. 56, 168–178. Zou, X., Wang, F.L., Chang, Y.Q., Zhang, B., 2017. Process operating performance optimality assessment and non-optimal cause identification under uncertainties. Chem. Eng. Res. Des. 120, 348–359.