Comprehensive economic index prediction based operating optimality assessment and nonoptimal cause identification for multimode processes

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

Contents lists available at ScienceDirect

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

Comprehensive economic index prediction based operating optimality assessment and nonoptimal cause identification for multimode processes Yan Liu a , Fuli Wang a,b,∗ , Yuqing Chang a,b , Ruicheng Ma c a

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

a r t i c l e

i n f o

a b s t r a c t

Article history:

For many multimode processes, the process operating performance may deteriorate with

Received 31 August 2014

time from optimal state due to process disturbances, noise, and other uncertainties, and it

Received in revised form 15 January

is important to develop an effective operating optimality assessment method; however, it

2015

has not yet been paid sufficient attention and few researches have been reported in this area

Accepted 4 March 2015

so far. In this study, a novel comprehensive economic index (CEI) prediction based operating

Available online 14 March 2015

optimality assessment and nonoptimal cause identification method is proposed for contin-

Keywords:

transitional mode on account of their different process characteristics. In stable mode, the

Operating optimality assessment

CEI is predicted by some common methods and then the optimality index is constructed

uous multimode processes. The assessment strategies are formulated for both stable and

Nonoptimal cause identification

based on the predicted CEI. In transitional mode, the CEI of a transition is predicted by the

Continuous multimode processes

weighted average of the CEIs of the similar historical transitions, and then the optimality

Tennessee Eastman process

index is calculated for online assessment of the transitional mode. When the operating performance is nonoptimal, the responsible cause variables can be identified by the proposed nonoptimal cause identification method. Finally, the feasibility and efficiency of the proposed strategies are demonstrated through the Tennessee Eastman (TE) process. © 2015 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

1.

Introduction

In industrial processes, the process operating optimality is a crucial problem that has attracted tremendous attention from both academia and industrial areas for the past decades (Camponogara and Scherer, 2011; Shen et al., 2014). However, because of the disturbance, noise, and other uncertain causes, industrial processes are not likely to be maintained all the time on the optimal point given by optimization

Additionally, due to load change, feedstock variation and parameter drift, etc., many industrial processes, including both continuous and batch processes, have the multimode characteristic (Zhang et al., 2013; Yu and Qin, 2008; Tan et al., 2012; Wang et al., 2012). For example, the Tennessee Eastman (TE) process is a continuous multimode process and simulates a chemical production process with four kinds of reactants A, C, D, E and an inert gas B to produce two kinds of liquid products G and H. To get different product, the operating conditions should be adjusted accordingly, which causes different stable modes, and the adjustment process from one stable mode to another can be considered

approach, which discounts the benefits of preliminary designs for process optimization and results in the degraded operating behaviors. Therefore, it is crucial to achieve the online operating optimality

as the transitional process. As another example, i.e. injection-molding

assessment for the industrial processes.

process, which is a batch process and a polymer processing technique.

∗ Corresponding author at: College of Information Science & Engineering, Northeastern University, 3 Lane 11, Wenhua Road, Heping District, Shenyang, Liaoning 110819, China. Tel.: +86 024 83687434; fax: +86 024 23890912. E-mail address: [email protected] (F. Wang). http://dx.doi.org/10.1016/j.cherd.2015.03.008 0263-8762/© 2015 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.

78

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

The entire production will go through four stages: injection, packingholding, plastication, and cooling. The process characteristics within each stage are relatively stable and can be considered as stable modes,

but it is time-consuming work and only suitable for the processes with a small number of modes. Wang et al. (2012) proposed the

while the process between two stages is considered as transitional mode. In this study, only continuous multimode process is considered.

method which can save online computing time and enhance realtime performance of online mode identification.

The stable mode is the main process which occupies the most production time to yield high productivity, and the good operating performance of a stable mode is the guarantee for the improvement of the comprehensive economic index (CEI). Owing to the unstable production conditions, however, the transitional process usually produces the substandard products, which may affect the CEI. If nonoptimal transitional mode can be evaluated timely and accurately, managers and operators can take appropriate action to adjust and improve the transition process, which can shorten the transition duration, reduce scrap rates and economic losses. In short, both stable and transitional mode are the indispensable parts of a multimode process, and the online operating optimality assessment for multimode processes is meaningful for further production operating adjustment and process performance improvement. Nowadays, the operating performance assessment for multimode processes mainly focus on safety. Operating safety assessment is to ensure production safety, human health and avoid environmental pollution, etc., and some methods have been applied to this aspect recently (Ye et al., 2009; Lin et al., 2013). In particular, multimode process monitoring (Wang et al., 2012; Yu and Qin, 2009; Haghani et al., 2014; Zhang and Li, 2013), as a kind of effective means to ensure the process operating safety, is also considered to belong to the scope of safety assessment. By way of comparison, the purpose of operating safety assessment or process monitoring is to distinguish between safe and unsafe, or normal and fault, whereas operating optimality assessment is to answer “whether the operating performance is optimal or nonoptimal under normal operating conditions”. Therefore, process operating optimality assessment is above the level of operating safety and makes a deeper understanding and mastering with the process operating performance for the operators and managers. Through operating optimality assessment, they can propose reference suggestions for further operating adjustment and performance improvement. However, to the best of the authors’ knowledge, few works have been reported on operating optimality assessment for multimode industrial processes. Recently, Ye et al. proposed a probabilistic framework of operating safety and optimality assessment for continuous multimode industrial processes (Ye et al., 2009). They used a Gaussian mixture model (GMM) to characterize multiple stable modes and constructed the safety and optimality indices (SI and OI) in each mode respectively. Then, a hierarchical-level classification method was presented to divide these indices into different performance levels, and margin analysis on each level was introduced. Nevertheless, they only considered the operating performance assessment for stable modes and neglected that for transitional modes which objectively exist in multimode processes. In addition, when the operating performance is nonoptimal, only experience-based qualitative analysis is provided rather than a general quantitative analysis. That is to say, an automatic nonoptimal cause identification method has not been explored. It is true that the operating optimality assessment for stable and

mode transformation probability-based online mode identification

In this paper, only the existing methods of offline and online mode identification are used and the main efforts are focused on the online operating optimality assessment and nonoptimal cause identification for both stable and transitional modes. It is widely accepted that the process operating optimality is closely associated with the CEI which may be one of the important production indices, such as costs, profits, total revenue and production efficiency, etc. or the weighted integration of several production indices. If the CEI approaches to or reaches the history optimal level, we believe that the process operating performance is optimal. Thus, it is suitable to evaluate the process operating performance based on CEI. However, CEI is difficultly achieved online and usually determined by offline statistical analysis, which may introduce a large delay and discount the timeliness of online implementation, so the easy-measured variables can be used to predict the CEI for online assessment. To predict the quality in the stable mode, various data-driven prediction methods have been developed and applied to the industrial processes (Liu et al., 2011; Wang, 2011; Ding et al., 2011; Wang et al., 2010). Then, the online assessment method for stable mode is proposed based on the predicted CEI. Because of the dynamics of transitional modes, the CEI of a transition is derived from the accumulation effects of the whole transition. That is to say, in transitional mode assessment, the whole transition information should be used in CEI prediction, and multiway partial least squares or other prediction algorithms for batch processes are seemed to be a choice to develop the prediction model. Nevertheless, the characteristics of the transitional mode limit the application of those algorithms in CEI prediction. On the one hand, because the transitional modes do not play the primary role in the continuous industrial processes compared with the stable modes, it is difficult to collect adequate transition data to develop the model, which affects the accuracy and reliability of the prediction model. On the other hand, transitions usually have uneven-length durations within the same kind of transitional mode, and the durations of nonoptimal transitions may be much longer than those of optimal ones, which causes the prediction algorithms for batch processes difficult to be applied to transitional process. Basing on the fact that the transitions with similar operating performance usually have the similar transition trajectories, the CEI of the new transition can be predicted by the weighted average CEI of the similar transitions. Therefore, the CEI prediction method based on the similar transition trajectories is given for transitional mode, and the process operating performance of the transitional mode is evaluated on account of the predicted CEI. The proposed CEI prediction method for transitional mode can easily achieve the purpose of CEI prediction for online transition and do not need to synchronize the transitions with different lengths no matter for offline or online. Also, when the process operating performance is nonoptimal, the contribution-based nonoptimal cause identification methods for both

transitional modes as well as the nonoptimal cause identification is the core part of operating performance assessment for multimode process,

stable and transitional mode are proposed for identifying the responsible cause variables. The idea of the proposed nonoptimal cause identification is similar to that of the contribution plots based fault

but considering the integrity of the method, some issues have to be

diagnosis method (Peng et al., 2013; Liu et al., 2013). Then, the managers

addressed and are listed as follows:

and operators can take appropriate operating adjustment strategies by combining their experience and the results of nonoptimal cause

(i) Offline mode identification for modeling data The common treatment for this issue is to assume that the offline modeling data for stable and transitional modes can be separated with the help of the expert knowledge and process information before modeling. Encouragingly, some automatic mode identification algorithms have been developed in recently years (Tan et al., 2012; Wang et al., 2012). (ii) Online mode identification for online application The online mode identification is to determine the mode type of the current process and then select the right model for online implementation. To try each mode in turn is an alternative way,

identification for production improvement. The schematic diagram of proposed method is shown as Fig. 1. The rest of this paper is organized as follows. First, some preparatory theoretical supports are framed by revisiting to the existing offline and online mode identification algorithms. Subsequently, the prediction models for each of the stable modes are established and the optimality index based on the predicted CEI is calculated for online operating optimality assessment. In Section 4, the online assessment strategy for transitional mode is developed, where the CEI is predicted based on the similar transition trajectories. For nonoptimal operating performance, the nonoptimal cause identification

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

2.2.

Fig. 1 – Schematic diagram of proposed method.

approach is proposed in Section 5. In case study, a typical multimode process, Tennessee Eastman (TE) process, is studied to demonstrate the feasibility and efficiency of the proposed method. Finally, the paper ends with some conclusions and acknowledgements.

2.

Preliminaries

2.1.

Offline mode identification algorithm

Offline mode identification is to classify the modeling data into different data sets corresponding to different modes. It is the basis of operating optimality assessment for multimode processes. So far, several offline mode identification algorithms have been proposed (Tan et al., 2012; Wang et al., 2012; Lu et al., 2004). Recently, Tan et al. (2012) achieve the offline mode identification based on the underlying changing of the process characteristics. They considered that the local process characteristic will be largely similar within the same mode whereas those of different modes are significantly distinct. Starting from stable mode, they divided the offline data into a series of data windows with the same length firstly. Then, the loading matrices which indicate the underlying process characteristics of each data window were extracted by PCA, and the similarities between the loading matrices of each data window and that of the reference window were calculated sequentially using the similarity index   defined by Tan et al. (2012), (Pk , Pref ) = 1 −

J  k ref  p − pj  /2J, where Pk and Pref denote j=1  j

the loading matrices of data window k and a stable mode ref reference window, pkj and pj are the jth column of Pk and Pref , respectively, and J is the number of process variables. If the similarity is larger than a predefined threshold, they are considered to belong to the same stable mode; otherwise, it is believed that the process characteristics of these two windows are different and the process changes from stable mode to transitional mode. To highlight the details of mode type changeovers, a smaller data window is introduced and the process characteristics of each small data window are extracted again for comparing the similarities. This method can precisely describe the local process characteristics by adjusting the window length. Furthermore, not only can it deal with the large amounts of modeling data quickly but also it can characterize the changes in the transition area detailed. Thus, it will be employed in this study.

79

Online mode identification algorithm

Before online operating optimality assessment, the right mode should be selected for the current process. The general way is to try all the modes and choose the mode with the monitoring statistics (i.e. I2 , T2 and SPE) below the corresponding control limits as the target one, but it would lead to great computational burden and decrease the effectiveness of online assessment. Recently, mode transformation probability (MTP)based online mode identification is proposed to guide the online mode selection (Wang et al., 2012). MTP is a probability parameter describing the frequency of transformation between two stable modes. It is believed that if the historical data are enough, the mode transformation behaviors must be hidden in these data, and one can obtain and use the MTP to guide the online mode identification. In addition, to distinguish the transitional mode from the process fault, the process monitoring results are combined with MTP together for online mode identification. Assume A and B are two stable modes, MTP from mode A to mode B is defined as P(B|A) = nAB /nA-all , where nA-all is the frequency of all modes transforming from A, and nAB is the frequency of transforming only from A to B. If the mode type at moment k − 1 is stable mode A, and the monitored statistics at moment k do not exceed the corresponding control limits, it means that the current process is normal. Otherwise, the attempt begins with the transitional mode starting from stable mode A with the maximum MTP. Then, the sample of moment k is monitored again using the model of the corresponding transitional modes sequentially. If all models do not fit the process characteristics, a fault or an unmodeled new mode may appear. Based on the above procedures, the online mode identification is achieved and one can choose the corresponding mode for further online operating optimality assessment. This method can save online computing time and enhance real-time performance. Here, we use it for online mode identification.

3. Online prediction and assessment for stable mode 3.1. Comprehensive economic index prediction for stable mode Based on the offline mode identification procedure, the normal operating data for different stable modes are classified into several data sets. Assume that there are C stable modes, and the modeling data corresponding to stable mode c (c = 1, 2, T . . ., C) is denoted as (Xc , yc ), where yc = [yc (1), yc (2), . . ., yc (Nc )] ∈ RNc ×1 represents the CEI in stable mode, which may be one of the important production indices, such as costs, profits, total revenue and production efficiency, etc. or the weighted integration of several production indices. Xc = T [xc (1), xc (2), . . ., xc (Nc )] ∈ RNc ×J is constituted by the process variables which are closely associated with the CEI; Nc is the number of samples, and xc (i) and yc (i), i = 1, 2, . . ., Nc , represent the ith sample in Xc and yc , respectively. It is assumed that Xc and yc have been centered to zero mean and scaled to unit variance before modeling. Consider the fact that the CEI is difficult to be achieved online and a lot of process measured variables which are closely related to CEI are easily obtained, we will use the measurable process information to predict the CEI and then achieve the purpose of online operating optimality

80

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

assessment for stable mode. So far, many data-driven methods have been developed for quality prediction in the stable production processes, such as partial least squares (PLS) (Witten et al., 2010), principal component regression (PCR) (Huang and Yang, 2012), artificial neural networks (ANN) (Yu et al., 2012), and support vector machine (SVM) (To et al., 2014), etc. Without loss of generality, the prediction relationship between the process variable xc (i) and the CEI, yc (i), can be denoted as yˆ c (i) = fc (xc (i)),

i = 1, 2, . . ., Nc ,

(1)

where fc (·) represents the relationship especially for stable mode c, and yˆ c (i) is the predicted value of yc (i). For the new sample xnew ∈ RJ×1 , yˆ cnew can be calculated as yˆ cnew = fc (xnew ).

(2)

Given the ability of PLS in terms of handling large numbers of highly correlated variables, measurement errors, and missing data and reducing the dimensionality of data space, it has been widely used in many industrial processes. Taking PLS as an example, the CEI prediction model for stable mode c is given as follows: yˆ c = Xc bc ,

(3) −1

where bc = XTc Uc (T Tc Xc XTc Uc ) T Tc yc = [bc,1 , bc,2 , . . ., bc,J ]T ∈ RJ×1 is the regression coefficient vector between Xc and yc ; Tc and Uc are the score matrixes with appropriate dimension, respectively. The detailed PLS algorithm can be found in several references (Höskuldsson, 1988; Wold et al., 2001). Then, yˆ cnew = xTnew bc . It is worth noting that the focuses of this study is to achieve the purpose of online operating optimality assessment based on the idea of comprehensive economic index prediction. Although an appropriate prediction method can enhance the effectiveness of the assessment method, how to select one of the most suitable prediction methods is not the main scope of this study. To evaluate the predicted abilities of the model, several performance indices are introduced. The root-mean square error (RMSE) and the maximum relative error (MRE) are defined as follows:

 RMSE =

Ntest i=1

Step 1: At moment k, construct the online data window T Xnew = [xnew (k − H + 1), . . ., xnew (k)] ∈ RH×J . Step 2: Normalize Xnew with the mean and standard deviation of stable mode c obtained during offline modeling. Then, the normalized data window and its mean vector are T denoted as Xcnew = [xcnew (k − H + 1), . . ., xcnew (k)] and x¯ cnew =

k

xc (h)/H, respectively. h=k−H+1 new Step 3: Calculate yˆ cnew as follows:

yˆ cnew = fc (x¯ cnew ).

(y(i) − yˆ (i))

2

Ntest

,

(4)

(5)

where Ntest is the number of test samples, y(i) and yˆ (i) are the real and predicted values, respectively. The RMSE is termed the root mean square error in calibration (RMSEC) for the calibration data and the root mean square error in prediction (RMSEP) for the test data. Similarly, the MREC and MREP are the maximum relative errors for calibration and test data, respectively.

Online assessment for stable mode

Assume that based on the online mode identification method introduced in Section 2, the new sample has been assigned to one of the stable modes, and then the corresponding prediction model is used for predicting the CEI.

(6)

Step 4: Without loss of generality, we assume that the CEI used here is optimal when it is small, and then the optimality index for stable mode c, OIsc , is defined as follows:

OIsc =

   y(i) − yˆ (i)   , MRE = max  y(i)  i

3.2.

Considering that a single sample is not enough to characterize the process operating performance and is susceptible to outliers, a data window with width H is introduced as the basic analysis unit for stable mode CEI prediction and online assessment, and its width depends on the actual situation of the production process. If the process runs smoothly and contains few singular values and outliers, H should be set relatively short. In contrast, if the process is easy to be disturbed by external environment, H should be selected relatively long to reduce the influence of disturbance and improve the reliability of assessment result. The wider the window width H, the less sensitive to the process changes, which leads to the longer time delay for the assessment from one operating performance to another. On the contrary, the narrower the window, the shorter the time delay. However, the assessment result may become more sensitive to the process interference, and the reliability of which is reduced. Firstly, the CEI corresponding to online data window is predicted using the prediction model as in Eq. (1). Then, the optimality index for stable mode c is calculated based on the predicted CEI. Generally, the process operating performance can be evaluated as optimal or nonoptimal. To strictly distinguish between optimal and nonoptimal, an assessment threshold is introduced. Combined with the threshold and the optimality index, the process operating performance in stable mode c can be evaluated by the following procedures:

⎧ 1, ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩

1− 0,

if yˆ cnew < ycmin yˆ cnew − ycmin

ycmax − ycmin

,

if ycmin ≤ yˆ cnew < ycmax ,

(7)

if yˆ cnew ≥ ycmax

where ycmin and ycmax are the minimum and maximum values of yc (i), i = 1, 2, . . ., Nc . Seen from Eq. (7), we know that the value OIsc is between 0 and 1. When OIsc is close to 1, it means that the process operating performance is optimal as yˆ cnew is close to ycmin . On the contrary, if OIsc is close to 0, it indicates that the process operating performance is likely to be nonoptimal. For strict distinction between optimal and nonoptimal, an assessment threshold  is introduced, and its value can be determined based on the historical data and expert experience. According to the historical data, the maximum and minimum values of the CEI can be counted, and then a dividing value of the CEI between optimality and nonoptimality can be given by the expert. Furthermore,  can be calculated through replacing the predicted CEI with the dividing value in Eq. (7).

81

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

Fig. 3 – Data representation of the uneven-length transitions.

Fig. 2 – Schematic diagram of online prediction and assessment for stable mode. Step 5: Assessment result is obtained by the following assessment rules: Case 1: If OIsc ≥ , it means that the CEI obtained under the current operating condition is good, and the process operating performance is optimal. Case 2: If OIsc < , we can say that the process is operating on the nonoptimal performance. Usually, the values of  between 0.5 and 1 are appropriate. If it is too high, the process operating under optimal performance may be evaluated as nonoptimal for the effects of interference. That is to say, the robustness of the assessment method can be reduced. However, if it is too low, the nonoptimal operating performance can be evaluated as optimal one, which increases the rate of error evaluation and reduces the accuracy of the assessment result. The online operating optimality assessment procedures for stable mode are shown in Fig. 2. Actually, the managers always hope the process operates on optimal performance, and when the operating performance is nonoptimal, some valuable references can be provided for operating adjustment. Therefore, it is necessary to propose an effective nonoptimal cause identification method, and this problem will be discussed in Section 5.

4. Online prediction and assessment for transitional mode In multimode processes, transitional mode serves as a bridge linking two stable modes. Different from the stable mode with constant characteristics (Tan et al., 2012), the transitional mode usually has dynamic characteristics. In addition, due to process disturbances, operating adjustments and physical constraints, etc., the transition durations, for the same kind of transitional mode, do not have the same length of time. Considering the fact that most of the products during the transitional mode are wastes, enterprises usually hope the transition duration as short as possible, which would decrease the operation costs to a great degree. Generally, the durations of normal transitions should be within a certain range (Tan et al., 2012), and one with too long duration is likely to be running on a nonoptimal operating performance. Anyway, we are

able to find a CEI for operating optimality assessment in transitional mode, such as the process operation costs, transition duration, etc. In this section, the CEI prediction method based on the similar transition trajectories is proposed. Owing to the fact that the transitions with similar operating performance usually have the similar transition trajectories, the CEI of the new transition can be predicted by the weighted average CEI of those similar transitions. Firstly, the similar transitions ˜ new = [x˜ new (1), . . ., x˜ new (k)]T ∈ Rk×J , of the current transition, X ˜ new contains are selected from historical transition data set. X incomplete transition information only from the transition beginning to the moment k. Several indices, such as distance (Hu et al., 2013), integrated factor with both distance and angle (Cheng and Chiu, 2005), and correlation (Fujiwara et al., 2009), etc. have been used to evaluate the similarity between data sets. However, this is not the main scope of the present paper. The distance with the consideration of variable weighs to CEI is utilized to measure the similarity in this study, not only for its simplicity but also because of its effectiveness for most cases. If other similarity indices are used instead of the distance, the essence of the operating optimality assessment for transitional mode based on CEI prediction will not change. In addition, based on the assumption that the future information ˜ new from moment k + 1 to the end is of the current transition X consistent with those of the similar transitions, only the existing information, i.e. samples from the beginning to moment k, is used for similarity measure. Then, the CEI of the current ˜ new is predicted by those of the similar transitions. transition X The proposed CEI prediction method for transitional mode can easily achieve the purpose of CEI prediction for online transition and does not need to synchronize the transitions with different lengths no matter for offline or online, which can avoid distorting or losing the original process information or introducing interference information. Assume that there are M transitional modes, and Im ˜m ˜m )}i=1 represents the historical data set of transitional {(X i ,y i mode m (m = 1, 2, . . ., M), where Im is the number of transiIm Km ×J m T ˜m ˜m i ˜m ˜m ˜m )}i=1 ; X and y˜ m tions in {(X i ,y i = [x i (1), . . ., x i (Ki )] ∈ R i i stand for the process data and the CEI of the ith transition, m respectively; Kim is the sample number, and Kmin = min(Kim ) i

m and Kmax = max(Kim ) stand for the minimum and maximum i

transition durations, respectively. Fig. 3 shows the representation of the uneven-length transitions. To remove the impact of the dimension, it is necessary to ˜m normalize the historical transitions X i , i = 1, 2, . . ., Im , m = 1, 2, . . ., M, before similarity measure. The normalization process ˜m for transition X i is listed as follow: m (h) ← x˜ i,j

m (h) − x m ˜ i,j,min x˜ i,j m m − x˜ i,j,min x˜ i,j,max

,

j = 1, 2, . . ., J,

h = 1, 2, . . ., Kim , (8)

82

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

where x˜ new,j (h) is the jth variable of x˜ new (h). To facilitate the presentation, the normalized transition is still denoted as ˜ new . X ˜ new between normalized X Step 3: Calculate the similarity dm i m ˜ and Xi as follows: m , (1) If k ≤ Kmin

dm i

= exp

k

−

= exp

 [x˜ new (h) −

T ˜ new (h) x˜ m i (h)] W[x

−



h=1 k

J

−

h=1

x˜ m i (h)]/ˇ

[wj (x˜ new,j (h) −

2 m x˜ i,j (h)) ]/ˇ

,

i = 1, 2, . . ., Im ,

j=1

(10)

where W = diag(w1 , w2 , . . ., wJ ) is the weights coefficient matrix of process variables with respect to the CEI, wj is the weights of the jth variable to the CEI and can be determined according to the process knowledge and expert experience, or some data mining methods. In this study, they are determined based on the process knowledge. ˇ is a user-defined adjustable parameter and can be set by trial and error and expert experience; m m , only the transitions with the durations (2) If Kmin < k ≤ Kmax no shorter than k are used for similarity measure, and the similarity is calculated the same as in Eq. (10). Step 4: Select the historical transitions with the similarities larger than similarity threshold ı (0 < ı < 1) and constitute I I ˜m m ˜m m the similar transition set {(X - i , y-˜ i )}i=1 , where {(X - i , y-˜ i )}i=1 ⊂ m m Im ˜ i , y˜ )} , I ≤ Im . Then, the CEI corresponding to X ˜ new can {(X i i=1 be predicted as follows:

Fig. 4 – Schematic diagram of online prediction and assessment for transitional mode. m (h) is the jth variable of x m m (h)] ˜ i,j,min ˜m where x˜ i,j = min[˜xi,j i (h), x h

m m (h)] are the minimum and maximum valand x˜ i,j,max = max[˜xi,j h

m

˜ i , respectively. The data pre-treatment ues of variable j in X procedure shown as in Eq. (8) can not only eliminate the influence of dimension on the transitional trajectory similarity, but also avoid synchronizing the uneven-length transitions, which may introduce the interference information or loss and distort the original transition information. Before online assessment, we assume that the transitional ˜ new has been determined mode type of the current transition X based on the online mode identification method introduced in Section 2, and then the corresponding historical transitions are used for similarity measure and CEI prediction. The framework of online CEI prediction and operating optimality assessment for transitional mode is exhibited in Fig. 4, and the details are given by the following procedures: ˜ new = Step 1: Construct the current transition X T m [x˜ new (1), x˜ new (2), . . ., x˜ new (k)] . If k > Kmax , it means that ˜ new has exceeded the upper limit of the the duration of X normal range. Thus, it is evaluated as nonoptimal directly, and the assessment process is finished. Otherwise, go to step 2. ˜ new with the maximum and minimum Step 2: Normalize X values of each process variable of the historical transition to ˜ new and be compared. For example, if the similarity between X ˜m ˜ new X is to be measured, the normalization of variable j in X i can be implemented as follow:

x˜ new,j (h) ←

m x˜ new,j (h) − x˜ i,j,min m m − x˜ i,j,min x˜ i,j,max

,

h = 1, 2, . . ., k,

(9)

m

y˜ˆ new =

I

dm y˜ m i=1 i - i

I

dm i=1 i

(11)

.

Step 5: The optimal index for transitional mode m, OItm , is defined as

OItm =

⎧ ⎪ 1, ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩

1−

m

if y˜ˆ new < y˜ m min m

y˜ˆ new − y˜ m min

˜m y˜ m max − y min

,

m , if y˜ m ≤ y˜ˆ new < y˜ m max min

(12)

m

if y˜ˆ new ≥ y˜ m max

0,

˜m where y˜ m = min(y˜ m ) and y˜ m ) are the minimum max = max( y min i i i

i

and maximum values of y˜ m , i = 1, 2, . . ., Im . i Step 6: Assessment result is obtained according to the following assessment rules: Case 1: If OItm ≥ , it means that the current transition trajectory is consistent with those of the optimal transitions and the process operating performance is evaluated as optimal; Case 2: If OItm < , we can conclude that current transition is under the nonoptimal operating performance. In the above process, the parameter ˇ and ı should be further clarified. If ˇ is too large, all the similarities, calculated by Eq. (10), between the transitional trajectories from same operating performance and different ones will become large and close to 1; if ˇ is too small, both of them will become small and approach to 0. That is to say, no matter the value of ˇ is too large or too small, it will increase the difficulty in selecting the similar transitions and distinguishing the different ones with respect to the online transition. Therefore, an appropriate

83

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

parameter ˇ should be set by trial and error and expert experience. ı is used to determine how many similar transitions in the historical data can be utilized in CEI prediction with respect to the online transition, and it can be set based on the of the historical data. For the transitions from similarity dm i among the same operating performance, the similarities dm i them usually fall into a range and the lower limit can be set as the similarity threshold ı, which can not only make sure more useful historical transition information is used in online prediction and assessment but also exclude the influence of dissimilar transitions. If ı is determined by the expert experience, and then the larger the assessment threshold ı is, the more similar transitional trajectories but less transition batches will be selected, which may loss the useful transition information and make the predicted CEI only decided by a handful of historical transitions with lack of general. On the contrary, the smaller the assessment threshold ı is, the more transition batches but less similar with the online transition will be selected, which introduces irrelevant information and reduces the accuracy of predicted CEI.

5.

Nonoptimal cause identification

ref

ref

c where x¯ c,j and x¯ new,j are the jth variables of x¯ c and x¯ cnew , respectively. The larger Cscj is, the jth variable is more likely to be the cause variable leading to nonoptimal operating performance.

5.2.

Variable contributions for transitional mode

Unlike stable mode where the relationship between the process variables and the CEI can be described precisely, such as in Eq. (1), it is difficult to depict this kind of relationships in the transitional mode. Thus, the variable contributions in transitional mode cannot be expressed as the partial derivative m of OItm or (y˜ˆ new − y˜ m ) to the virtual scale factor v. However, min based on the fact that the transition trajectories are significantly different for optimal and nonoptimal transitions, the deviation between them can be attributed to some of the process variables which deviate from their optimal operating trajectories. Therefore, the variable contributions can be calculated on the basis of the deviation between the transition trajectories of optimal and nonoptimal transitions. In transitional mode m, pick out the optimal transitions Im ˜m ˜m )}i=1 and construct the optimal reference transifrom {(X i ,y i m,ref

When the process operating performance is nonoptimal, it is meaningful to find the corresponding causes. Like the contribution plots based fault diagnosis method (Peng et al., 2013; Liu et al., 2013), the nonoptimal cause identification methods based on variable contributions are proposed for both stable and transitional mode in this section. The variable contributions to the optimality index OIs (OIt) are defined, and then the variables with relatively large contributions are considered as the possible variables leading to nonoptimal operating performance.

5.1.

˜ tions {X i

m,ref

, y˜ i

Im

m,ref

}i=1 . Ki

is the number of samples of the

m,ref ith optimal transition and Kmax

m,ref

= max(Ki i

). Then, the devi-

˜ new and that ation ˜ym between the transition trajectories of X of the optimal reference transition is defined as follows: m,ref

Case 1: if k ≤ Kmax ,

˜ym

=

k

T

[x˜ new (h) − x¯ m,ref (h)] W[x˜ new (h) − x¯ m,ref (h)]

h=1

=

Variable contributions for stable mode

J k

(15) 2 m,ref [wj (x˜ new,j (h) − x¯ j (h)) ],

h=1 j=1

For stable mode c, the process data corresponding to the minimum CEI, i.e. optimal operating performance data, are collected from Xc and constitute the optimal reference ref ref ref ref T = [xc (1), xc (2), . . ., xc (Nc )] ref number in Xc .

ref data Xc

ref Nc ×J

ref , where Nc

∈R is sample From Eq. (7), we can conclude that the variable contributions to OIsc is equivalent to those to yc , and yc = yˆ cnew − ref ycmin = fc (x¯ cnew ) − fc (x¯ c ) is the deviation of the CEI between the current process and the optimal operating performance. Thus, the variable contributions can be calculated only based on yc . Denote a virtual scale factor v as v = [v1 , v2 , . . ., vJ ]T , where vj = 1, j = 1, 2, . . ., J, and x   represents the vector of [x1 v1 , x2 v2 , . . ., xJ vJ ]T , where xj vj is the variation of variable xj . For the online data Xcnew , the contribution of variable j to yc is defined as ∂fc (x¯ cnew  ) ∂yc ∂fc (x¯ c  ) = = − , ∂j ∂j ∂j ref

where x¯ c

=

Ncref i=1

ref

ref

xc (i)/Nc

and x¯ cnew =

ref

are the means vectors of Xc and Xnew , respectively. Specific to the PLS, the variable contribution has the following form: T

Cscj =

ref

T

∂(x¯ cnew  ) bc ∂(x¯ c  ) bc ref c − = (x¯ new,j − x¯ c,j )bc,j , ∂j ∂j

(14)

is the mean vector of the m,ref

(h) is the jth variable

m,ref

Case 2: if k > Kmax , m,ref Kmax −1

˜ym

=

T

[x˜ new (h) − x¯ m,ref (h)] W[x˜ new (h) − x¯ m,ref (h)]

h=1

+

k

m,ref

T

m,ref

[x˜ new (h) − x¯ m,ref (Kmax )] W[x˜ new (h) − x¯ m,ref (Kmax )]

m,ref h=Kmax m,ref Kmax −1

J



h=1 k

+

xc (h)/H h=k−H+1 new

m,ref ˜h (h)/N x˜ n=1 n

˜ h is the number of transitions with the durations of x¯ m,ref (h), N m,ref ˜ m,ref (h) is the hth sample of X , n= longer than h, and x˜ n n ˜ h; 1, 2, . . ., N

(13)

k

N˜ h

reference transitions at moment h, x¯ j

=

ref

Cscj

where x¯ m,ref (h) =

(16) [wj (x˜ new,j (h) −

2 m,ref x¯ j (h)) ]

j=1 J



m,ref

[wj (x˜ new,j (h) − x¯ j

m,ref

2

(Kmax )) ],

m,ref j=1 h=Kmax

which is based on the fact that if the current transition m,ref is optimal, it should arrive at the set point x¯ m,ref (Kmax ) m,ref at moment Kmax , and nonoptimal operating performance is just because some variables do not reach the set point m,ref m,ref within an appropriate time. x¯ j (Kmax ) is the jth variable m,ref

of x¯ m,ref (Kmax ).

84

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

The variable contributions to ˜ym for case 1 and case 2 are defined as in Eqs. (17) and (18), respectively: Ctjm

=

Variable no.

∂˜ym ∂vj



∂

Table 1 – Process variables for operating optimality assessment.

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



k

T

[x˜ new (h)   − x¯ m,ref (h)  ] W[x˜ new (h)   − x¯ m,ref (h)  ]

h=1

=

∂

k J

 h=1

=

∂vj wj (x˜ new,j (h)vj −

2 m,ref x¯ j (h)vj )





j=1

∂vj k

m,ref

=2

[wj (x˜ new,j (h) − x¯ j

2

 

(h))] vj 

h=1

k

=2

m,ref

[wj (x˜ new,j (h) − x¯ j

vj =1

2

Description

Units

A feed (stream 1) D feed (stream 2) E feed (stream 3) A and C feed (stream 4) Recycle flow (stream 8) Reactor feed rate (stream 6) Reactor temperature Purge rate (stream 9) Product separator temperature Product separator pressure Product separator underflow (stream 10) Stripper pressure Stripper temperature Reactor cooling water outlet temperature Separator cooling water outlet temperature

kscmh kg/h kg/h kscmh kscmh kscmh ◦ C kscmh ◦ C kPa m3 /h kPa ◦ C ◦ C ◦ C

(h))] ,

h=1

(17)

Ctjm

=

schematic diagram of TE process. The process consists of five major units: a reactor, a condenser, a vapor-liquid separator, a recycle compressor, and a product stripper, and it involves

∂˜ym ∂j

⎡

∂⎣

⎤

m,ref

Kmax −1

[x˜ new (h)  ␯ − x¯ m,ref (h)  ] W[x˜ new (h)   − x¯ m,ref (h)  ]⎦ T

h=1

=

∂j

⎡

⎤

T ⎢ ˜ ⎥ m,ref m,ref [xnew (h)   − x¯ m,ref (Kmax )  ] W[x˜ new (h)   − x¯ m,ref (Kmax )  ]⎦ k

∂⎣

m,ref

h=Kmax

+

∂

⎧ m,ref J ⎨Kmax

−1 ⎩

=

∂

h=1

m,ref

[wj (x˜ new,j (h)j − x¯ j

2

⎫ ⎬

(h)j ) ]

⎭

j=1

⎧ ⎪ J k ⎨

⎪ ⎩

∂j

∂j m,ref

[wj (x˜ new,j (h)j − x¯ j

2

m,ref

⎫ ⎪ ⎬

(Kmax )j ) ]

⎪ ⎭

m,ref h=Kmax j=1

+

(18)

∂j m,ref

Kmax −1

=2

m,ref

[wj (x˜ new,j (h) − x¯ j

2

 

(h))] j 

h=1

=2

⎧ m,ref ⎪ ⎨Kmax

−1 ⎪ ⎩

h=1

m,ref

[wj (x˜ new,j (h) − x¯ j

2

vj =1

(h))] +

+2

k

m,ref

[wj (x˜ new,j (h) − x¯ j

m,ref

h=Kmax k

m,ref

[wj (x˜ new,j (h) − x¯ j

m,ref h=Kmax

6.

Case study: Tennessee Eastman process

6.1.

Process description and data preparation

To demonstrate the efficiency of the proposed method, TE process is introduced as the research background in this paper. TE process is the simulation of a complex continuous multimode industrial process (Downs and Vogel, 1993). As a benchmark, it has been widely used for the test of process control (Lawrence Ricker, 1996), process monitoring (Tan et al., 2012), and operating optimization (Shen et al., 2014), etc. Fig. 5 is the well-known

2

m,ref

 

(Kmax ))] j 

m,ref

(Kmax ))]

⎫ ⎪ ⎬ 2 ⎪ ⎭

vj =1

.

eight components: A–H. The gaseous reactants A, C, D, E, and the inert B are fed to the reactor while the liquid products G and H are formed. TE process contains two blocks of variables: 12 manipulated variables and 41 measured variables. For detailed process description, one can refer to Lee and Chiang et al. (Downs and Vogel, 1993). In this paper, the process operation costs is used as the CEI, and 15 variables which are closely related to the CEI are chosen from the 41 process measured variables for operating optimality assessment and tabulated in Table 1. All the process variables are measured with an interval of 0.01 h and 100 samplings are collected in one hour.

85

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

Fig. 5 – The schematic diagram of TE process (Ge et al., 2013). To produce the reasonable simulation data, the differences between the mode and the operating performance in a mode should be clarified beforehand, which can be summarized into the following two aspects. On the one hand, the change of mode is usually required by the production or production-driven for producing the products with different specifications, models, component contents, etc., which can be implemented by adjusting the main operational variables, whereas the change of operating performance from optimal to poor is usually caused by the disturbance, noise, and other uncertainties, which would lead to poor CEI, and it is undesirable for the productions. On the other hand, the set points of the operational variables of different modes are significant different; however, the process variables which are responsible for the change in the operating performance usually fluctuate within a small range. In view of the above analysis, four stable modes (A, B, C, D) and three transitional modes (A to B, A to D, and B to C) are considered. The imposed modifications for generating different stable modes are shown in Table 2 and arise from controlled changes in reactor pressure and reactor level (Tan et al., 2012). In addition, because the reactor temperature affects the full extent of the reaction and further involves in the process operation costs, we simulate the optimal and nonoptimal operating performances in stable

mode by adjusting the set point of reactor temperature, and the ranges of the operation costs calculated from TE process for each stable mode are also revealed in Table 2. To establish the prediction models, a total of 2000 samples, including both optimal and nonoptimal operating performance, have been collected for each stable mode B A B and denoted as (XA calibration , ycalibration ), (Xcalibration , ycalibration ), C D C D (Xcalibration , ycalibration ), and (Xcalibration , ycalibration ), respectively. The characteristics of three process variables, i.e. recycle flow (variable 5), reactor feed rate (variable 6), and product separator temperature (variable 9) are exhibited in Fig. 6, which is very clear that the process data coming from four different stable modes occupy different distribution areas. Known from the mechanism of TE process, the lower the reactor temperature is, the higher the operation costs is under the normal operating conditions for all the stable modes mentioned above, which results in nonoptimal operating performance. Fig. 7 shows the operation costs respond to the changes of the reactor temperature in each stable mode. In this study, PLS is employed for establishing the prediction models of each stable mode, and the numbers of latent components are automatically determined by cross-validation (Zhang, 2000). Then, 5 latent components are reserved for each PLS prediction model, respectively.

Table 2 – Four stable modes of TE process. Stable mode

Reactor pressure (kPa)

Reactor level (%)

Reactor temperature (◦ C) Optimal

A B C D

2800 2705 2705 2500

65 65 75 65

121.6 121.7 121.5 121.9

Operation costs ($/sampling)

Nonoptimal 111.9 111.9 111.9 113.9

38.29–105.09 26.89–79.16 42.59–122.45 49.57–57.64

Operation costs (Downs and Vogel, 1993) = (purge costs) × (purge rate) + (product stream costs) × (product rate) + (compressor costs) × (compressor work) + (steam costs) × (steam rate).

86

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

Table 3 – Transitional modes of TE process. Mode type

Operating performance

Adjustment strategy

Optimal

AB

Nonoptimal

Optimal

AD

Nonoptimal

Optimal

BC

Nonoptimal

Operation costs ($/tra)a

Duration (sampling)

P: 2800–2705 L: 65 P: 2800–2750–2705 L: 65

216–240

17,688–19,653

346–399

32,904–36,560

P: 2800–2650–2500 L: 65 P: 2800–2500 L: 65

730–810

62,823–69,803

913–1010

96,198–106,890

P: 2705 L: 65–75 P: 2705 L: 65–70–75

530–585

27,696–31,076

650–717

49,732–55,257

Operation costs = [(purge costs) × (purge rate) + (product stream costs) × (product rate) + (compressor costs) × (compressor work) + (steam costs) × (steam rate)] × duration. a

$/tra is a unit symbol and means the cumulative operation costs per complete transition.

Mode Mode Mode Mode

85

100 Stable mode A 50

A B C D

Operation cost($)

0

500

1000 Samplings

1500

2000

0

500

1000 Samplings

1500

2000

0

500

1000 Samplings

1500

2000

0

500

1000 Samplings

1500

2000

120 Stable 100 80 mode B 60 40

80 Variable 9

Reactor temperature ( C)

75

Stable 100 mode C 50

70 65 36 26

34

24 22

32 20 30

Variable 6

Variable 5

18

Fig. 6 – Three-dimensional data characteristic of four different stable modes.

Fig. 7 – Operation costs change with the reactor temperature in stable mode A–D.

Generally, different adjustment strategies between two stable modes will generate different transitional behaviors and influence the operating performance of the transitional mode. Thus, the optimal and nonoptimal operating performances for each transitional mode are produced according to different

Variable 1~4

Variable 9~12

Variable 13~15

adjustment strategies which are listed in Table 3. For simplicity, the transitional mode from A to B is named as transitional mode AB and so on. As shown in Table 3, the difference of optimal and nonoptimal operating performances for transitional mode AB is attributed to the different ways for changing the

0.25

5400

800

8.5

0.2

5200

700

8

0 Variable 5~8

120 Stable 100 mode D 80 60 40 20

200

5000

0

200

600

0

200

7.5

0

200

0 0 3200

200

25

40

130

1

20

35

120

0.5

15

0

200

30

0

200

110

100

2800

18

80

2700

17.5

60

0

200

2600

0

200

17

80

110

47

60

100

46

40

90

0

200

0

200

45

0

200

3000 0

0 Samplings

200

2800

0

200

200

Fig. 8 – Optimal (blue solid) and nonoptimal (red dotted line) transition trajectories in transitional mode AB. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

87

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

MREC

0.0199 0.0207 0.0258 0.0302

0.0663 0.1015 0.1266 0.1469

0.0012

0.0103

0.0026

0.0191

0.0031

0.0265

reactor pressure (P). For an optimal transition AB, the reactor pressure is directly adjusted from 2800 kPa to 2705 kPa, which makes the operation costs of a transition fluctuate between 17,688 $ and 19,653 $. For a nonoptimal transition AB, however, the reactor pressure is first regulated to 2750 kPa, and then adjusted to 2705 kPa after the process is stable, which leads to the transition costs fluctuate between 32,904 $ and 36,560 $ per transition. Pick out one of the optimal and nonoptimal transitions from the historical data of transitional mode AB, respectively, and then the comparison of transition trajectories under optimal and nonoptimal operating performances is presented in Fig. 8, from which it can be seen that the variable trajectories are significantly different. Similarly, the differences between optimal and nonoptimal operating performances for each of the other transitional modes are also caused by the different adjustment strategies on reactor pressure (P) or reactor level (L), and the transition trajectories are different. For each transitional mode, 50 transitions are severally generated to construct the historical data set. There are both optimal and nonoptimal transitional trajectories in these 50 transitions. The difference between optimal and nonoptimal transitions is for adjustment strategy, but the differences among the transitions belonging to the same performance grade are just for random noise. The historical data sets of transitional mode AB, AD, and BC ˜ AB ˜ AD ˜ AB ˜ AD are denoted as (X ), (X ), and calibration , y calibration , y calibration calibration

Predicted CEI Measured values

100 80 60 40 200

400

600

800 1000 1200 1400 Samplings (a)

1600

1800

200

400

600

800 1000 1200 1400 Samplings (b)

1600

1800

0.1 Relative error

RMSEC

0.05 0 -0.05 -0.1

Fig. 9 – Predicted CEI and relative errors for stable mode A.

3.6 Operation cost($)

(XA , yA ) calibration calibration (XBcalibration , yBcalibration ) (XCcalibration , yCcalibration ) , yD ) (XD calibration calibration AB ˜ AB ˜ (X , y ) calibration calibration ˜ AD ˜ AD (X ) calibration , y calibration ˜ BC ˜ BC (X ) calibration , y calibration

x 10

4

Predicted CEI Measured values

3.58 3.56 3.54 3.52

2000

2050

2100

2150 2200 Samplings (a)

2250

2300

2350

2000

2050

2100

2150 2200 Samplings (b)

2250

2300

2350

0.01 Relative error

Calibration data

Operation cost($)

Table 4 – RMSE and MRE for the calibration data.

0.005 0 -0.005 -0.01

BC

For online test, additional 6560 samplings are produced according to the adjustment strategies listed in Tables 2 and 3 and constitute the test data set Xtest ∈ R6560×15 and ytest ∈ R6560×1 , which includes three stable modes and two transitional modes experienced optimal or (and) nonoptimal operating performance. A detailed description of the change in mode and operating performance in the test data is given in Table 5. Before online assessment, the test data have been assigned to a certain mode by the online mode identification algorithm proposed by Wang et al. (2012), and the online mode identification results are tabulated in Table 5, through which it can be found that the online mode identification results are almost the same with the actual situation. To ensure the accuracy of the predicted CEI is the key to online operating optimality assessment. Figs. 9–13 report the comparisons of the operation costs and the predicted

80 Operation cost($)

6.2. Online assessment and nonoptimal cause identification

Fig. 10 – Predicted CEI and relative errors for transitional mode AB.

Predicted CEI Measured values 60 40 20

2400

2600

2800

3000

3200 3400 Samplings (a)

3600

3800

4000

4200

2600

2800

3000

3200 3400 Samplings (b)

3600

3800

4000

4200

0.2 Relative error

˜ calibration , y˜ BC ), respectively. The performance indices for (X calibration the calibration data are summarized in Table 4. For each of the prediction models, the values of RMSEC and MREC are all very small, which indicates that the predictive abilities for the calibration data are satisfactory.

0.1 0 -0.1 -0.2 2400

Fig. 11 – Predicted CEI and relative errors for stable mode B.

88

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

Table 5 – Test data and online mode identification results. Mode type

Performance grade

Stable mode A Transitional mode AB Stable mode B Transitional mode BC Stable mode C

Operation cost($)

3.1

x 10

Mode identification results (sampling)

Optimal Nonoptimal

1–520 521–1984

Nonoptimal

1985–2350

1989–2353

Nonoptimal Optimal

2351–3545 3546–4217

2354–4220

Optimal

4218–4762

4221–4766

Optimal Nonoptimal

4763–5271 5272–6560

4767–6560

1–1988

4

Predicted CEI Measured values

3.05 3 2.95 2.9

4300

4400

4500 Samplings (a)

4600

4700

4300

4400

4500 Samplings (b)

4600

4700

0.04 Relative error

Actual situation (sampling)

Table 6 – RMSE and MRE for the test data.

0.02

0

-0.02

Fig. 12 – Predicted CEI and relative errors for transitional mode BC. A ˜ AB ˜ AB ), values based on PLS for test data sets (XA test , ytest ), ( Xtest , y test C C ˜ BC ˜ BC (XBtest , yBtest ), (X ), and (X , y ), respectively, together test , y test test test with their corresponding prediction relative errors. The relative error is described as REk = (yk − yˆ k )/yk . As can be seen, most of the samplings have been well estimated for

Test data

RMSEP

MREP

A (XA test , ytest ) AB AB ˜ (Xtest , y˜ test ) (XBtest , yBtest ) ˜ BC ˜ BC (X test , y test ) (XCtest , yCtest )

0.0226

0.0906

0.0015 0.0271

0.0121 0.1234

0.0057 0.0301

0.0306 0.1575

stable modes; although the prediction relative errors at the beginning of transitional mode AB and BC are relatively large, which is because the similar process characteristics of optimal and nonoptimal transitions in the early stage, the predicted CEI approaches to real values gradually with the increasing of the sampling number. Table 6 summarizes the performance indices for the test data, which shows that both RMSE and MRE are small enough and within the acceptable range for the real productions. Set the parameters H = 10,  = 0.8, ˇ = 2500, ı = 0.7, and W = diag(0.0340, 0.0117, 0.0528, 0.0510, 0.0031, 0.0606, 0.2062, 0.0996, 0.1395, 0.0560, 0.0213, 0.0086, 0.1315, 0.1186, 0.0055). Then, the online assessment and nonoptimal cause identification results for stable mode A and transitional mode AB are displayed in Figs. 14 and 15, respectively. Fig. 14(a) shows the online assessment results based on OIsA for stable mode A. At moment 523, OIsA = 0.7902 < 0.8, which indicates that the process operating performance changes from optimal to nonoptimal at this moment, and the assessment result is accordance with the actual situation. That is to say, the process operating performance can be evaluated timely and accurately using our proposed method. Furthermore, the variable contributions at moment 523 are calculated using the method proposed in Section 5.1 and shown in Fig. 14(b). It can be seen from Fig. 14(b) that variable 7 as well as variables 8, 9, 12, 13 has relatively large contributions than the others. Known from the mechanism of TE process, there are relatively strong coupling relationships between the reactor temperature and the other process variables with large contributions. The reaction rate is a function of temperature through an Arrhenius expression. Thus, if the reactor temperature deviates from the

Nonoptimal

80 60 40 4800

5000

5200

5400

5600 5800 Samplings (a)

6000

6200

0.2

Relative error

1 0.8 0.6 0.4

0.1 0 -0.1 5000

5200

5400

5600 5800 Samplings (b)

6000

6200

A

6400

Fig. 13 – Predicted CEI and relative errors for stable mode C.

OIsA

0.2

Threshold

0 200

400

600

15 10

800 1000 1200 1400 1600 1800 Samplings (a)

At moment 523

5 0 -5 -10

-0.2 4800

At moment 523: OIs =0.7902

6400 Variable contributions

Operation cost($)

Predicted CEI Measured values

100

Optimality index OIs A

Optimal

120

1

2

3

4

5

6

7 8 9 Variables (b)

10 11 12 13 14 15

Fig. 14 – Online assessment result for stable mode A and variable contributions.

89

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

1

0.6

OItAB

0.4

Threshold

0.2 0 2000

2050

2100

2150 2200 Samplings

2300

2350

(a)

200

Acknowledgments

At moment 1989

100 50 0

1

2

3

4

5

6

7 8 9 Variables (b)

10 11 12 13 14 15

Fig. 15 – Online assessment for transitional mode AB and variable contributions.

optimal set point, the reaction rate will be affected, and it will lead to the reactant flowing to the next subprocess before the ratio of the component in the reactor meeting the production requirements. To accommodate the change, the subsequent production subprocesses, i.e. vapor-liquid separator and product stripper, will make some automatic adjustments which are reflected in separator temperature (variable 9), stripper pressure (variable 12), and stripper temperature (variable 13). In addition, the change in the ratio of the component will finally affect the purge rate (variable 8). All of these variables affect the operation costs directly or indirectly. Therefore, the nonoptimal cause identification result is consistent with the actual situation. The online assessment and nonoptimal cause identification result for transitional mode AB is shown in Fig. 15. It is clear that the assessment results are always correct even though the predictions are inaccurate at the beginning of the transition. Furthermore, variables 7, 9, 13, and 14 with larger contributions are considered as the responsible cause variables for this nonoptimal operating performance, for which one can refer to Fig. 8 for further confirmation. Similarly, the online assessment and nonoptimal cause identification results for stable modes B, C, and transitional mode BC all agree with the actual facts. For the simulation results of them, please see Appendix A.

7.

Conclusions

This paper proposes a novel operating optimality assessment and nonoptimal cause identification method for continuous multimode processes based on CEI prediction. In stable mode assessment, the CEI can be predicted by some common methods for stable production processes and then the optimality index is calculated to evaluate the process operating performance in real time; in transitional mode assessment, the CEI of a new transition is predicted by weighted average of CEI using the similar transitions selected from the historical data. Also, the corresponding nonoptimal cause identification method is developed for both stable and transitional mode. Then, the proposed approach is applied to TE process. It is concluded from the example that the proposed online assessment framework is able to give an efficient operating

We are grateful to the anonymous reviewers for their constructive comments and suggestions which greatly improved the quality and presentation of this paper. Also, we appreciate the financial support from Project 863 of China (No. 2011AA060204), Project 973 of China (2009CB320601), National Natural Science Foundation of China (Nos. 61374146, 61174130 and 61304055), State Key Laboratory of Synthetical Automation for Process Industries Fundamental Research Funds (No. 2013ZCX02-04), the Fundamental Research Funds for the Central Universities (No. N120404020), and the Research Program of Shanghai Municipal Science and Technology Commission (No. 13dz1201700).

Appendix A. The online assessment and nonoptimal cause identification results for stable modes B, C, and transitional mode BC are shown in Figs. A1–A3, respectively.

Nonoptimal Optimality indexOIsB

150

Optimal

1 0.8 At moment 3553: OIs =0.7998

0.6

B

0.4

OIsB

0.2

Threshold

0 2400

Variable contributions

Contributions

2250

optimality assessment for both stable and transitional mode. When the process operating performance is nonoptimal, the responsible process variables can be identified by the proposed nonoptimal cause identification method, which can provide the reference for operators and managers to propose reasonable strategies for further production operating adjustment and process performance improvement.

2600

2800

3000

6 4

3200 3400 Samplings (a)

3600

3800

4000

4200

At moment 2354

2 0 -2

1

2

3

4

5

6

7 8 9 Variables (b)

10 11 12 13 14 15

Fig. A1 – Online assessment for stable mode B and variable contributions.

1 0.8

OItBC

OItAB

0.8

OItBC

0.6

Threshold 0.4 0.2 0 4300

4400

4500 Samplings

4600

4700

Fig. A2 – Online assessment for transitional mode BC.

90

chemical engineering research and design 9 7 ( 2 0 1 5 ) 77–90

Variable contributions

Optimality index OIs C

Optimal

Nonoptimal

1 0.8 0.6 0.4

At moment 5288: OIs =0.7961

OIsC

C

0.2

Threshold

0 4800

5000

5200

5400

10

5600 5800 Samplings (a)

6000

6200

6400

At moment 5288 5

0

-5

1

2

3

4

5

6

7 8 9 Vari ables (b)

10 11 12 13 14 15

Fig. A3 – Online assessment for stable mode C and variable contributions.

References Camponogara, E., Scherer, H.F., 2011. Distributed optimization for model predictive control of linear dynamic networks with control-input and output constraints. IEEE Trans. Autom. Sci. Eng. 8, 233–242. Cheng, C., Chiu, M., 2005. Nonlinear process monitoring using JITL-PCA. Chemometr. Intell. Lab. Syst. 76, 1–13. Ding, J., Chai, T., Wang, H., 2011. Offline modeling for product quality prediction of mineral processing using modeling error PDF shaping and entropy minimization. IEEE Trans. Neural Netw. 22, 408–419. Downs, J.J., Vogel, E.F., 1993. A plant-wide industrial process control problem. Comput. Chem. Eng. 17, 245–255. Fujiwara, K., Kano, M., Hasebe, S., Takinami, A., 2009. Soft-sensor development using correlation-based just-in-time modeling. AIChE J. 55, 1754–1765. Ge, Z., Song, Z., Wang, P., 2013. Probabilistic combination of local independent component regression model for multimode quality prediction in chemical processes. Chem. Eng. Res. Des. 92, 509–521. Haghani, A., Jeinsch, T., Ding, S.X., 2014. Quality-related fault detection in industrial multimode dynamic processes. IEEE Trans. Ind. Electron. 61, 6446–6453. Höskuldsson, A., 1988. PLS regression methods. J. Chemom. 2, 211–228. Hu, Y., Ma, H., Shi, H., 2013. Enhanced batch process monitoring using just-in-time-learning based kernel partial least squares. Chemometr. Intell. Lab. Syst. 123, 15–27. Huang, S., Yang, J., 2012. Improved principal component regression for face recognition under illumination variations. IEEE Signal Process. Lett. 19, 179–182. Lawrence Ricker, N., 1996. Decentralized control of the Tennessee Eastman challenge process. J. Process Control 6, 205–221.

Lin, Y., Chen, M., Zhou, D., 2013. Online probabilistic operational safety assessment of multimode engineering systems using Bayesian methods. Reliab. Eng. Syst. Safety 119, 150–157. Liu, Y., Wang, F., Chang, Y., Li, C., 2011. A SNCCDBAGG-based NN ensemble approach for quality prediction in injection molding process. IEEE Trans. Autom. Sci. Eng. 8, 424–427. Liu, Q., Qin, S.J., Chai, T., 2013. Decentralized fault diagnosis of continuous annealing processes based on multilevel PCA. IEEE Trans. Autom. Sci. Eng. 10, 687–698. Lu, N., Gao, F., Wang, F., 2004. Sub-PCA modeling and online monitoring strategy for batch processes. AIChE J. 50, 255–259. 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. Shen, Y., Hao, L., Ding, S.X., 2014. Real-time implementation of fault-tolerant control systems with performance optimization. IEEE Trans. Ind. Electron. 61, 2402–2411. Tan, S., Wang, F., Peng, J., Chang, Y., Wang, S., 2012. Multimode process monitoring based on mode identification. Ind. Eng. Chem. Res. 51, 374–388. To, A., Paul, G., Liu, D., 2014. Surface-type classification using RGB-D. IEEE Trans. Autom. Sci. Eng. 11, 359–366. Wang, D., 2011. Robust data-driven modeling approach for real-time final product quality prediction in batch process operation. IEEE Trans. Ind. Informat. 7, 371–377. Wang, D., Liu, J., Srinivasan, R., 2010. Data-driven soft sensor approach for quality prediction in a refining process. IEEE Trans. Ind. Informat. 6, 11–17. Wang, F., Tan, S., Peng, J., Chang, Y., 2012. Process monitoring based on mode identification for multimode process with transitions. Chemometr. Intell. Lab. Syst. 110, 144–155. Witten, J.M., Park, S., Myers, K.J., 2010. Partial least squares: a method to estimate efficient channels for the ideal observers. IEEE Trans. Med. Imag. 29, 1050–1058. Wold, S., Sjöström, M., Eriksson, L., 2001. PLS-regression: a basic tool of chemometrics. Chemometr. Intell. Lab. Syst. 58, 109–130. 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. Yu, J., Qin, S.J., 2008. Multimode process monitoring with Bayesian inference-based finite Gaussian mixture models. AIChE J. 54, 1811–1829. Yu, J., Qin, S.J., 2009. Multiway Gaussian mixture model based multiphase batch process monitoring. Ind. Eng. Chem. Res. 48, 8585–8594. Yu, Y., Choi, T., Hui, C., 2012. An intelligent quick prediction algorithm with applications in industrial control and loading problems. IEEE Trans. Autom. Sci. Eng. 9, 276–287. Zhang, J., 2000. Multivariable Statistical Process Control. Chemical Industry Press, Beijing, China. Zhang, Y., Li, S., 2013. Modeling and monitoring between-mode transition of multimodes processes. IEEE Trans. Ind. Informat. 9, 2248–2255. Zhang, Y., An, J., Ma, C., 2013. Fault detection of non-Gaussian processes based on model migration. IEEE Trans. Control Syst. Technol. 21, 1517–1526.