Online Probabilistic Assessment of Operating Performance Based on

10912

Ind. Eng. Chem. Res. 2009, 48, 10912–10923

Online Probabilistic Assessment of Operating Performance Based on Safety and Optimality Indices for Multimode Industrial Processes Lubin Ye, Yuming Liu, Zhengshun Fei, and Jun Liang* State Key Laboratory of Industrial Control Technology, Institute of Industrial Control, Zhejiang UniVersity, Hangzhou 310027, People’s Republic of China

Operating performance of industrial process on safety and optimality may deteriorate with time due to process characteristic variation, and it is crucial to develop strategies for online operating performance assessment. Although there have been some studies and applications on process safety assessment, optimality assessment has not yet been paid sufficient attention. This paper proposes a probabilistic framework of online operating assessment for industrial processes. First, a Gaussian mixture model (GMM) is used to characterize multiple operating modes. Considering the distribution of process variables, safety and optimality indices (SI and OI) are defined and calculated by two successive nonlinear mappings. A hierarchical-level classification method is then presented to divide these indices into different performance levels, and margin analysis on each level is introduced. Finally, performance prediction and preliminary suggestions for improvement are provided. The proposed assessment strategy is then applied in two examples: Tennessee Eastman Process (TEP) and polypropylene (PP) production process, which indicate the efficiency of the proposed approach. 1. Introduction Process operating safety and optimality are two crucial issues in industry that have attracted tremendous attention from both academia and industry for the past several decades. Nevertheless, it is well-known that the industrial process operating performance of safety and optimality may deteriorate as the process characteristics drift with time, which cancels the benefits of preliminary designs for process safety and optimization and results in a degraded operating behavior. Therefore, assessment of operating performance on both safety and optimality are crucial for process maintenance. Conventional industrial process safety analysis consists of qualitative methods like process hazards analysis (PHA), hazard and operability analysis (HAZOP), and failure mode and effects analysis (FMEA)1-3 and quantitative approaches like fault tree analysis (FTA),4,5 Markov analysis, and reliability block diagram.6 Recently, some interdisciplinary methods including Bayesian analysis,7,8 complex networks,9 and fuzzy logic10 were adopted to assist the safety assessment. The detection and diagnosis of abnormal events has also been an active research area,11-13 especially on the study of data driven methods.14,15 The data driven approaches are usually termed as multivariate statistical process control (MSPC), which can avoid the extremely difficult and time-consuming work of constructing complex mechanistic models for industrial processes. On process operating optimization, many algorithms and applications have been presented.16-23 However, because of the existence of disturbance, noise, and other uncertain causes, real processes are not likely to be maintained all the time on the optimal point given by optimization approaches. Hence it is important to get a measure on how far the current operating condition is from the optimum (or how optimal the current operating state is). According to the knowledge of the authors, very little work has been published on the optimality assessment for industrial process operating. * To whom correspondence should be addressed. Tel.: (+86571)8795-1071-424. Fax: (+86-571)8795-1068. E-mail: jliang@ iipc.zju.edu.cn.

It is not enough to only give whether the process operating is safe or unsafe or optimal or bad for the work of performance assessment. And the question “how safe and how optimal is the current operation?” has to be answered. A hierarchical-level method is proposed in this work, which classifies the operating performance into several levels. Following the idea of stability margin in classical controller design, the margin, representing the relative position on each safety or optimality level, is formulated to give a more particular characterization. The practical industrial processes may run under substantially different operating conditions due to load change, feedstock variation, parameter drift, etc., which lead to a complex multimodal distribution of the operating data. Gaussian mixture model (GMM) is a widely used descriptive model for multimodal distributions in data mining discipline with the parameter estimated via the expectation maximization (EM) algorithm.24,25 It has recently been validated to be effective in multimode process monitoring.26-28 In the current work, GMM is used to characterize the multiple operating modes of industrial processes. In this paper, a probabilistic framework for online operating performance assessment for the multimode process is proposed, consisting of two procedures: (1) off-line training and (2) online assessment. In off-line training, steady-state historical data are collected and the GMM is constructed characterizing the multiple operating modes with the parameters estimated via the EM algorithm. The appropriateness of the Gaussian model is subsequently verified by a normality test. Then, safety index (SI) and optimality index (OI) are separately defined by probabilistic methods to evaluate the operating performance of a local operating mode. Global indices are then integrated weighted by the posterior probabilities belonging to each operating mode, which are calculated by Bayesian inference. The global indices can then be interpreted as two successive nonlinear mappings. Subsequently, to get a definite and comparable assessment of the current operating state, several hierarchical SI and OI levels are obtained. In the mapped lower-dimensional D space, observations of different levels are divided by dis-

10.1021/ie801870g CCC: $40.75  2009 American Chemical Society Published on Web 10/15/2009

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009

criminant hyperplanes. In online performance assessment implementation, the monitored samples are first classified into a certain SI or OI level and then the margins on the levels are calculated, based on which predictions and preliminary suggestions for practical operating improvement are given. The remainder of the paper is organized as follows. Section 2 gives a brief introduction of the application of GMM in multimode process and EM algorithm. The specific operating performance assessment strategy including definitions and calculations of safety and optimality indices, level division, and margin analysis are presented in section 3. The effect of the proposed method is demonstrated by two examples: Tennessee Eastman process (TEP) and polypropylene production process in section 4. Lastly, conclusions are outlined in section 5.

1. E-step: The posterior probabilities of the ith training sample in the sth iteration are calculated by Bayesian law: (s) (s) ω(s) k pk(xi |µk , Σk )

P(s)(Ck |xi) )

K

∑

(s) (s) where ω(s) k , µk , and Σk are the estimated parameters in the sth iteration and Ck denotes the kth Gaussian component. 2. M-step:

n

∑ ω p (x|θ ) k k

k

(1)

k)1

where x ∈ Rm denotes the m-dimensional vector of measured process variable, ωk represents the weight of the kth component with ∑kK) 1ωk ) 1, K is the number of Gaussian components, and θk ){µk,Σk} represents the local mean and covariance of the components. The local Gaussian distribution probability density is expressed by pk(x|θk) )

∑P

(s)

µ(s+1) k

)

(Ck |xi)xi

i)1 n

(4)

∑P

(s)

K

p(x|θ) )

1 1 exp - (x - µk)TΣ-1 k (x - µk) 2 (2π)m/2 det(Σk)1/2 (2)

[

]

Each component characterizes the local variation of the corresponding operating mode. The weight ωk can be interpreted as a prior probability that an arbitrary sample comes from the kth Gaussian component. To construct a GMM, the parameters ωk,µk,Σk have to be estimated. Denote X ){x1,x2,...,xn}; the training data set and expectation-maximization (EM) algorithm is used here for the estimation problem, which has two iterative steps including expectation calculation and maximization as follows:25,28

(3)

(s) (s) ω(s) k pk(xi |µk , Σk )

k)1

2. Gaussian Mixture Model Most real industrial processes are running at multiple operating conditions and have a complex multimodal distribution, as driven by local operating modes. If theses modes are regarded as classes, from which the multimodal distributed process data observations are obtained, the characterization of the process operating conditions can be seen as a classification problem with unknown labels. It was stated that an arbitrary probability density can be approximated by a mixture model,26,29 which is also a methodology that addresses the clustering problem when class labels are missing.24 So a mixture model is quite fit here to characterize the operating data distribution of industrial processes. It is crucial to choose the probability density function structure for the components of the mixture model. According to central limit theorem, the mean of N random variables asymptotically follows a Gaussian distribution. In each local operating region, the measurements of process variables are thought to be made up of the sum of multiple independent causes and thus multivariate Gaussian distribution is a reasonable assumption for the process data from a local operating mode.24 Therefore, GMM is adopted in this paper for characterizing the data driven by multiple operation modes. The GMM probability density can be formulated by a weighted sum of local Gaussians25-28

10913

(Ck |xi)

i)1

n

∑P

(s)

) Σ(s+1) k

(Ck |xi)(xi - µ(s+1) )(xi - µ(s+1) )T k k

i)1

(5)

n

∑P

(s)

(Ck |xi)

i)1 n

∑P

(s)

ω(s+1) k

)

(Ck |xi)

i)1

n

(6)

The above two steps are iterated until the likelihood converges. The GMM representation eqs 1 and 2 and EM algorithm eqs 3-6 are based on the fundamental requirement that the process is at steady state in each operating region. So it is necessary to identify the steady state before the application of GMM and EM estimation. An F test type approach30 is adopted, where two types of estimates for variance are calculated by the exponentially weighted moving average. The ratio of the two estimates (denoted by R) is taken as the testing statistic, which is compared with a threshold determined by its probability distribution. Finally, Dempster’s rule of combination is used to get the overall index.31 To get sufficient accuracy of estimates for GMM parameters, a simple criterion is formulated to provide guidance on the minimum number of observations required for training sets. The Crame´r-Rao inequation gives the lower limit of estimation covariance,32 which is formulated as ∧ 1 Cov(θ) g M-1 n

(7)

where θˆ is the estimate of parameter θ, n is the number of observations, and M is the Fisher information matrix. For each Gaussian component, only the mean estimation is considered here for simplicity, and the Fisher information matrix can be obtained as M ) Σk-1. Therefore, the estimation covariance satisfies Cov(µˆ ) g 1/n Σk, implying that Var(µi) g 1/n σi,i,k2 for the ith mean estimator, where σi,i,k2 is the ith element on the diagonal of Σk. Then a relative index for estimation variance is formulated as Var(µi) σii,k2

g

1 n

(8)

It can be easily inferred that a large number of observations can reduce the lower bound of estimation variance. Then for a

10914


practical problem, the minimum number of training observations for a prescribed estimation precision can be obtained through eq 8. In the current work, the relative value of the estimation variance is fixed at 0.005. So the minimum number of training observations for each Gaussian component is 200. A critical question is whether the obtained Gaussian components in GMM are appropriate to represent the local distributions of real process data. In the current work, this can be verified through a normality testing based on higher order statistics of kurtosis and skewness. Theoretically, kurtosis and skewness are 0 for a Gaussian distribution.33,34 The further their values are from 0, the more the probability density deviates from Gaussian distribution. In normality testing, kurtosis and skewness are calculated and compared with their thresholds determined by their standard deviation.33,34 If either statistic exceeds its confidence threshold, the distribution is considered to be non-Gaussian. It will be shown later in examples that the data from a real polypropylene production process belonging to the same operating mode do not strictly follow Gaussian distributions. However, the deviations are slight and the Gaussian approximation is considered to be appropriate here.24 3. Proposed Approach Because the practical processes seldom operate precisely on the given operating points, it is essential to consider the data distributions while assessing the operating performance.35 When the process is running at an operating mode, the local probability distribution is characterized by the kth Gaussian component of the previously obtained GMM as x 0 N(µk, Σk)

(9)

where x ∈ Rm is the vector of process variables with µk ∈ Rm and Σk ∈ Rm×m representing the mean and covariance, respectively. The m-dimensional space spanned by x is defined as the operating space. Operating performance assessment approach is developed below. 3.1. Safety Index of Individual Operating Modes. Without loss of generality, the safety here is described by the safety limits for process variables, within which the process operating is considered to be safe and acceptable. The safety limits are obtained according to the process constraints and represented by equations in the m-dimensional hyperspace gi(x) ) 0

(10)

where i ) 1,2,...,l and l is the total number of process safety limits. The expressions gi(x) are linear or nonlinear functions and are standardized so that gi(x) e 0 denotes the variable x being within the safety limit. For the kth operating mode, contours of probability density are distributed around the operating mean point and can be depicted by the annular hyperellipsoids in the high-dimensional space, as seen in Figure 1. The first contour reaching the safety limits is considered as the critical contour, denoted by Bk. The process is thought to be safe if the variables are kept inside Bk. Thus the safety index of the kth operating mode is defined by the probability SIk ) Pr{x is inside Bk}

(11)

It can be inferred from eq 2 that under the same operating condition, the probability density is determined only by the exponential term (x - µk)T(Σk)-1(x - µk), which is equivalent

Figure 1. Contours of probability density (solid lines) and safety limits (dashed lines).

to the Mahalanobis distance from the center to the points on the hyperellipsoid as dM(x, µk) ) (x - µk)T(Σk)-1(x - µk)

(12)

So the hyperellipsoids in Figure 1 are also contours of Mahalanobis distances. Let δ denote the Mahalanobis distance between the critical contour and the mean point, then the safety index in eq 11 can be rewritten as SIk ) Pr{x is inside Bk} ) Pr[dM(x, µk) e δ]

(13)

In the high-dimensional space, it is intractable to find the critical contour, which reaches the safety limits first, by geometric methods. However, it can be seen from Figure 1 that this task is equivalent to finding the minimum Mahalanobis distance on the safety limits. Therefore, this work is transformed to a nonlinear optimization problem with its objective function formulated as min J ) (x - µk)T(Σk)-1(x - µk)

(14)

subject to gi(x) ) 0. Let x* denote the optimal solution of the optimization problem, and δ can then be obtained by δ ) (x*)TΣk-1(x*)

(15)

The probability in eq 13 is calculated by the following method. First, the principal component analysis of the process variables is obtained after centering and scaling as t ) PTx

(16)

where t is the score vector and P is the loading matrix. Then the process variable vector can then be reconstructed by x ) Pt, and the expression of Mahalanobis distance becomes dM(x, µk) ) tT(Λk)-1t

(17)

where Λ, which is diagonal, represents the covariance of t. Let u ) Λ-0.5t, and the Mahalanobis distance can be rewritten by dM(x, µk) ) uTu ) u12 + u22 + ... + um2

(18)


where u1,u2,...,um are the components of vector u. It is easy to obtain that these components are independent and follow the same standard Gaussian distribution. Therefore, the statistic dM(x,µk) follows a χ2 distribution with m degrees of freedom. Then the safety index in eq 13 can be calculated through χ2 distribution law by SIk ) Pr{tTΛk-1t e δ} ) Pr{χ2 e δ}

OIk )

∫

x∈Ck

pk(x|θk) J(x) dx

j

(21)

j

The term Pr(x ∈ ∆j) is the probability of variables being in the micro-unit ∆j, which can be estimated by the frequency Pr(x ∈ ∆j) ) nj /Nk

(22)

where nj is the number of samples in ∆j and Nk is the total number of samples with maximum memberships to the kth mode. Then the sample mean of objective function in ∆j, denoted by jJ(x ∈ ∆j), is calculated by

∑ J(x )

jJ(x ∈ ∆j) ) 1/nj

k

(26)

t

where the index GI here is representing SI/OI for safety/ optimality. The samples are not classified into deterministic operating modes as such a probabilistic strategy can avoid the risk of misclassification and the subsequent biased assessment results. Substituting the P(Ck|xt) term with the expression of eq 25 in eq 26, GI(xt) can then be formulated as K

∑ω

k

exp(-dk /2)GIk |det(Σk)| -1/2

k)1 K

GI(xt) )

(27)

∑ω

k

exp(-dk /2)|det(Σk)|

-1/2

where dk denotes the Mahalanobis distance between xt and µk. LetD ) [d1,d2,...,dK]T, and eq 27 can be simplified to GI(xt) )

∑

(28)

Φ

Ψ

Φ

Φ

Rm f Rk f R1 x

1 J(xi) OIk ) Nk x ∈C

rT exp(-D/2) βT exp(-D/2)

where r and β are the vectors of multipliers for the sum calculations in the numerator and denominator, respectively. Herein the K-dimensional space spanned by vector D is named as the D space. It can be seen that the global SI and OI calculation is essentially made up of two successive nonlinear mappings. First, Mahalanobis distances between the sample and each operating mean point are obtained to form the vector D, which is actually a mapping from the m-dimensional operating space to the K-dimensional D space. Then SI and OI are calculated by eq 28, implying a mapping from the K-dimensional D space to the 1-dimensional index space. The mapping structure is shown below:

(23)

t

Substitute eqs 22 and 23 in eq 21, then OIk is estimated by

f D f SI or OI

where

(24)

Φ(x) ) [d1(x), d2(x), ..., dK(x)]T

k

3.3. Global Safety and Optimality Indices. Considering that the samples may come from multiple operating modes, the posterior probabilities of an arbitrary monitored sample xt coming from each component can then be obtained by the constructed GMM and Bayesian inference: P(Ck |xt) )

k

k)1

xt∈∆j

j

∑ GI P(C |x )

k)1

∑ Pr(x ∈ ∆ ) jJ(x ∈ ∆ ) j

K

GI(xt) )

(20)

where pk(x|θk) is the probability function of the local operating mode. The multidimensional integral in the above equation is difficult to calculate and thus a discretization approach is utilized. The integral space is segmented into micro-units represented by ∆1,∆2,...,∆j,... and eq 20 is approximated by the discretized form: OIk )

belonging to the operating regions. Then the global safety and optimality indices, that is, SI and OI, are defined as the sum of SIk and OIk metrics across all operating modes weighted by corresponding posterior probabilities. Then the global indices are calculated by

(19)

3.2. Optimality Index of Individual Operating Modes. For simplicity of presentation, it is assumed that the objective function to be optimized, such as plant cost, profits, and product quality, is formulated by J(x). Agarwal and Huang35 used expected profits as an assessment of MPC performance, and the probabilistic approach is adopted here for definition of optimality index. When the process is running under a certain operating mode with the mean and covariance estimated by EM algorithm, the optimality index of the kth operating mode (OIk) is defined as the expectation of the objective function:

10915

ωkpk(xt |µk, Σk)

(25)

K

∑ ω p (x |µ , Σ ) k k

t

k

k

k)1

where the weight ωk is taken as prior probability. It is easily derived that the sum of the posterior probabilities of a certain sample satisfies ∑kK ) 1P(Ck|xt) ) 1, and the posterior probabilities are interpreted as memberships of the samples

and Ψ(D) )

rT exp(-D/2) βT exp(-D/2)

3.4. Level Division of Safety and Optimality Indices and Margin Analysis. The safety and optimality indices calculated above are absolute values, which lack comparability. For example, a value of 0.9 for the safety index is considered very safe for one process, but might be not safe enough in another process, which emphasizes safety much more. As for the optimality index, there are various targets for different processes and thus no uniform OI reference exists. Furthermore, the calculated SI and OI using the above methods are varying due to perturbation caused by uncertain factors like disturbance

10916


and noise. To give a definite evaluation of the current operating performance with references, a performance level division approach that expresses different grades of safety or optimality is proposed here. A general procedure can be implemented to determine the hierarchical performance levels. The range of the calculated SI and OI values of historical data is separated into several intervals, corresponding to which the operating performance levels are defined. It means that the process is running on the relevant performance level when SI or OI values are within a certain interval. The performance levels may represent several safety grades including the safest, barely safe, and unsafe, or optimality grades including optimal, suboptimal, and nonoptimal. The demarcation points for these intervals are decided according to the process characteristics and plant personnel’s attitudes of these concepts. However, there are no unified criteria for the determination of performance levels since the calculated results of the SI/OI do not have the same range for different industrial processes and there exist differences in process requirements and plant personnel’s preferences. The data belonging to different performance levels can be seen as being from different classes, separated by division planes, which are determined by a classification approach utilizing the training data. While in the high-dimensional operating space, these classes are linearly inseparable due to the complicated mapping in eq 27. However, the problem becomes quite easier in the lower-dimensional D space, where the D samples belonging to different performance levels (i.e., classes) are located in respective approximate conical zones illustrated by later experiments and can be easily divided by hyperplanes. In D space, {Ω1,Ω2,...,Ωp} is denoted as the classes corresponding to the p levels for SI or OI. It is inferred that p!/(2!(p - 2)!) hyperplanes are needed to separate p classes. The minimum perceptron criterion function approach is utilized here to determine these hyperplanes, which is described as follows. For conical regions, the angles between vectors are the most obvious features and are thus taken as separating criterion. Let {D1,D2,...,DN} represent the set of N training samples in D space from two classes Ωi and Ωj, and let vector h ) [h1,h2,...,hK]T denote the normal direction of the division hyperplane. Without loss of generality, the norm of h is restricted to 1. The sine of the angle between Dn and the hyperplane (denoted as β), which is also the cosine of the angle between Dn and h (denoted by R), is calculated by sin β ) cos R )

hTDn h 2· Dn

||| |

Y)

[| | | | | | ] y1

y1

,

2

y2

y2

, ...,

2

yN

yN

T

2

A perceptron criterion function36 is then constructed representing the error ratio by JP ) Yh - b - |Yh - b|

|

|

(32)

2

where b is the allowance vector. The vector h is obtained by minimizing eq 32 subject to the constraint |h|2 ) 1. Then the division hyperplane, which is also the limit between the two levels of SI or OI is determined by h and the coordinate origin. The margin is adopted here to measure the distance from the sample point to the limit of the current level (or class). Since angles are distinct features in conical regions, the margin is defined by the sine value of the angle between the sample and the division plane as λij or γij ) |sin βij | )

|| | | | | hTij Dn

hij 2· Dn

(33)

2

where hij denotes the normal direction of division hyperplane that separates the classes Ωi and Ωj in D space, and λij and γij are margins of the data sample for safety and optimality, respectively. Further, relative values are used instead of absolute ones to normalize the margins to a range of [0, 1]. The above constructed margins give measurements on the reliability that the operation remains on the current SI or OI level. If the margin is close to 0, a small change in operating may push the process beyond the boundary resulting in rising to a higher level or falling to a lower level. In contrary, the values approaching 1 imply that the process may tolerate a large variation while staying on the current SI or OI level. In the work of online process operating assessment, the performance levels and the margins are obtained, on the basis of which preliminary suggestions are available for operating performance improvements. Plant personnel are supposed to

(29) 2

and it is assumed that the sine values satisfy sin β > 0, sin β < 0,

D n ∈ Ωi D n ∈ Ωj

(30)

Then the samples are standardized to [y1,y2,...,yN] that yn ) Let

{

Dn, for Dn ∈ Ωi -Dn, for Dn ∈ Ωj

(31) Figure 2. Schematic diagram of the proposed framework for operating performance assessment.


10917

Figure 3. Schematic diagram of TEP.

manipulate the process to adjust the D vector to the regions corresponding to higher performance levels, based on the relationship between the process variables x and the Mahalanobis distance vector D ) Φ(x) ) [d1(x),d2(x),...,dK(x)]T. However, there exists a problem of how to manipulate the process to a better condition. This problem is actually the inverse approach of the work in this paper. Although the number of operable variables is smaller than that of process variables, the inverse mapping will have infinite mathematical solutions. This problem deserves further study and will be discussed in future work. The schematic diagram of the above proposed framework for process operating performance assessment is shown in Figure 2, and the detailed procedures are summarized as follows: (1) Collect a set of historical training data under all possible operating conditions. (2) Verify whether the collected data are at steady state. Only the steady-state data are used as training data to estimate GMM parameters. (3) Construct a Gaussian mixture model and estimate the parameters via EM algorithm in eqs 3-5. (4) Check whether the obtained Gaussian components are suitable for characterizing each operating mode through normality testing based on kurtosis and skewness. (5) Obtain individual safety and optimality indices (SIk and OIk) based on the local distribution by probabilistic approaches through eqs 19 and 24. (6) Calculate global safety and optimality indices (SI and OI) of the training samples by the weighted sum in eq 26. (7) The SI and OI are divided into several hierarchical levels and corresponding classes are separated in D space with the hyperplanes determined by minimizing the criterion of eq 32. (8) In online monitoring implementation, first the Mahalanobis distances between the testing sample and mean points of each operating mode are calculated. (9) SI and OI values of the current process operating are given by eq 28. (10) The current process operating performance is classified into a certain level and the corresponding margin is calculated by eq 33. (11) Finally, performance trend predictions are obtained in margin plots and some preliminary useful suggestions for improvement are provided.

Figure 4. TEP example: steady-state identification results of training data.

4. Examples 4.1. Tennessee Eastman Challenge Problem. The Tennessee Eastman process (TEP) benchmark was developed by Downs and Vogel,37 based on an industrial process of Tennessee Eastman Chemical Company, and is widely used for the test of process control design, monitoring, optimization, etc. The schematic diagram is illustrated by Figure 3. The process consists of five units: a reactor, a product condenser, a vapor-liquid separator, a recycle compressor, and a product stripper. Four gaseous reactants A, C, D, and E are fed into the reactor, where two products G, H with a byproduct F are generated by two gas-liquid exothermic reactions. The process has 12 manipulated variables, 22 continuous measurements, and 19 composition measurements sampled less frequently. All the measurements include Gaussian-type noise. For the sake of safe operation, some constraints for process variables are provided as the safety limits. The operating cost was given as a function of the process variables. There are six operating modes at different G/H mass ratios.37 Since the open-loop process operation is unstable, several control strategies have been designed.38-40 The base control structure developed by McAvoy and Ye38 is adopted in this paper. Because of the feedback

10918


Figure 5. TEP example: scatter plot of training data. Figure 7. TEP example: scatter plot of conical regions of optimality levels in D space and division planes (dotted lines).

Figure 8. TEP example: OI levels of the testing samples.

Figure 6. TEP example: normality testing of the training data in each operating region.

control actions, the random noise affects the entire process system. So the system has stochastic components and is appropriate to be characterized by probabilistic approach. In this example, the TEP running under five operating regions in the vicinity of the base mode is studied with the G/H mass ratio at 50/50. The sampling time is 0.05 h, and the 41 measured variables are selected for assessment work. There are 750 samples under each operating condition that are collected as the training data set. The data steady-state identification results are shown in Figure 4. For most of the samples, the R statistic is below the threshold and is thus verified to be at steady state. Another set of 750 × 5 samples are used as testing data. In the training step, the data projected into 2 of the 41 axes are exhibited in Figure 5, indicating that the variables of the five operating conditions are very near to each other. However, the five Gaussian components obtained in the constructed GMM can well characterize operating conditions in each individual mode. The normality testing results are shown in Figure 6, from which it can be seen that the kurtosis and skewness of all the

41 variables in each mode are below their confidence thresholds. Therefore, the Gaussian model is fit to characterize the data distribution. The given variable constraints are taken as the safety limit for SI calculation, and the operating cost is regarded as the optimization objective to obtain OI. Whereafter, global SI and OI are acquired by eq 26. Since all the operating modes are very far away from the process safety limits and there are insignificant differences in SI values, safety level division is not necessary and will not be discussed in this case. The OI is divided into three ranked levels (level 3, 2, and 1 denoting the high, medium, and low grade, respectively) with the distributions in D space illustrated in Figure 7. The separating hyperplanes are determined by minimizing the criterion in eq 32. In the online testing, after obtaining Mahalanobis distances between the monitored samples and the five operating mean points, OI is calculated. The results of OI level classification are plotted in Figure 8. In each level, margins of the testing data in the optimality levels are then obtained and exhibited in Figure 9. The first row shows that the optimality of process operating is on the level 3, and the operating conditions of sample 1500-2250 have larger margins than the sample 3000-3750, which indicates that the first 750 samples are more reliable on level 3. In the second row, although all the 1500 samples express a condition on level 2, the first 750 samples are closer to level 1 while the latter 750 ones are closer to level 3. So the process operating state of the latter 750 samples is more advisable than the first 750 ones. Moreover, it can be seen from Figure 9 that in each operating stage, the TEP operation


10919

Figure 9. TEP example: Optimality margin of the testing data.

Figure 10. Schematic process flow of polypropylene production process. Table 1. Process Variables for Modeling of PP Production Process variable 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

description concentration of hydrogen in reactor 201 concentration of hydrogen in reactor 202 reactor 201 density reactor 202 density recycle compressor gaseous flow TEAL flow electron donor feed flow reactor 201 propylene feed flow reactor 200 propylene feed flow pump 200 flush-propylene flow reactor 202 propylene feed flow flush-propylene flow in flash evaporator pump 201 flush-propylene flow pump 301 outlet propylene flow propylene feed flow inside the battery limits pump 201 power

units -6

10 10-6 kg/m3 kg/m3 kg/h kg/h kg/h t/h kg/h kg/h t/h kg/h kg/h t/h t/h kW

has no evident trend in optimality. Finally it is suggested that the process is manipulated to the operating region of sample 1500-2250.

variable

description

units

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

pump 202 power tank 202 level pressure of hydrogen main pipe reactor 200 pressure tank 200 pressure reactor 201 pressure emergency emission pressure from reactor 201 to 202 reactor 202 pressure temperature of tank 201 jacket water condenser 201 outlet temperature reactor 200 temperature reactor 201 temperature temperature of reactor 201 jacket water reactor 202 temperature temperature of reactor 202 jacket water

kW % MPa MPa MPa MPa bar MPa °C °C °C °C °C °C °C

4.2. Polypropylene Production Process. The polypropylene (PP) production process adopts the Spheripol technology. As illustrated in Figure 10, the process has three cascade reactors

10920


Table 2. Process Variable Constraints of PP Production Process constraints variables

low limit high limit

concentration of hydrogen in reactor 201 concentration of hydrogen in reactor 202 reactor 201 density reactor 202 density compressor 301 gaseous feed flow TEAL flow DONOR flow reactor 201 propylene feed flow reactor 200 propylene feed flow pump 200 flush-propylene flow reactor 202 propylene feed flow pump 201 power pump 202 power tank 202 level pressure of hydrogen main pipe reactor 200 pressure tank 200 pressure reactor 201 pressure reactor 202 pressure reactor 200 temperature reactor 201 temperature reactor 202 temperature flow ratio of TEAL/ propylene feed in reactor 201

50(10-6) 50(10-6) none none 1000 kg/h 1.0 kg/h none 6 t/h 2000 kg/h 370 kg/h 3 t/h none none 35% 2.5 MPa 3.1 MPa 3.1 MPa 3.1 MPa 3.1 MPa none none none 0.1

4200(10-6) 4200(10-6) 700 kg/m3 700 kg/m3 3500 kg/h 4.5 kg/h 2.5 kg/h 17 t/h none 430 kg/h 15 t/h 275 kW 275 kW 45% 3.5 MPa 3.7 MPa 3.5 MPa 3.65 MPa 3.65 MPa 25 °C 73 °C 73 °C 0.3

including two liquid-phase tubular loop reactors for homopolymerization and one gas-phase fluidized bed reactor for copolymerization. The catalyst, propylene, and hydrogen are fed to the first reactor. Then a portion of the slurry together with pure propylene is injected to the second reactor. The products of homopolymerization is delivered to the third reactor and mixed with ethylene to produce copolymers. The real PP production process is affected by many uncertainties including noise, disturbance, parameter drifts, etc. Hence, the data of process variables have extremely complex distributions. These uncertain factors are regarded as being caused by stochastic components and can thus be formulated by probabilistic methods. The process is able to produce various grades of PP resins. T30S is one of the primary products and is here taken as the example to demonstrate the effectiveness of the assessment strategy. There are 31 reflective process variables selected, listed in Table 1. The process safety limits are acquired from the plant operating handbook, and the details are listed in Table 2. The most important goal of the PP production is that the melt index (MI) of the product satisfies the quality specification. Thus the square of the difference between the actual MI and the target value is adopted here as the objective function. A total of 4500 samples of the variables from four producing phases are collected to constitute the training data set. The steady-state identification results are shown in Figure 11, indicating that not all the data are coming from stationary conditions. Then only the samples verified to be at steady state are used as training data to construct the GMM. The data distributions in original operating space are illustrated in Figure 12, from which it can be seen that data in different operating modes are severely overlapped. Four Gaussian components are obtained in GMM, and the normality testing results are shown in Figure 13, which demonstrates that the kurtosis and skewness of some variables exceed their Gaussian limits and the process data in each mode do not follow Gaussian distributions strictly. However, the excesses of most variables are slight and the Gaussian assumption is considered to be appropriate. After the work of SI and OI calculation, 3 levels for safety and optimality are divided,

Figure 11. PP production process example: steady-state identification results of training data.

Figure 12. PP production process example: scatter plot of training data.

respectively. Here the levels 3, 2, and 1 represent the high, medium, and low grades of safety or optimality, similar with that of TEP case. The data belonging to different levels or classes in D space are demonstrated in Figure 14 and 15, with the division planes determined by eq 32. A testing data set of 4500 samples is used to validate the efficiency of online assessment and margin analysis. Performance level results are illustrated in Figure 16. It can be seen that there are some glitches in the chart. This is not an accidental phenomenon and can be explained that on these samples, the process operating is very near to the level limit and jumping changes may easily happen on these samples. Safety and optimality margins are then calculated, shown in Figure 17 and 18, respectively. In Figure 17a, the process is first running at safety level 1, but quite near the limits between level 1 and 2. On sample 1080, it climbs to level 2. However, the margin curve shows that the process has a trend of falling and will not remain on level 2 for long. Then on sample 1150, it touches the limit again and falls back to level 1. Then during the period of sample 1339-1350, the process is oscillating between level 1 and 2, which accords with the fact that the margin on either level is near zero. Thereafter, the process falls to level 2 and remains on that level stably since the margin curve is leaving the limit gradually. From Figure 17b, it can be shown SI is first on level 2 and meanwhile approaching level 3. Then on sample 4103, the process switches to level 3 after a short period of oscillation. Thereafter, the process keeps a far distance from the boundary, which ensures the stable stay on level 3.


10921

Figure 13. PP production process example: normality testing of the training data in each operating region.

Figure 14. PP production process example: scatter plot of conical regions of safety levels in D space and division planes (dotted lines).

Figure 15. PP production process example: scatter plot of conical regions of optimality levels in D space and division planes (dotted lines).

In Figure 18, analogous results are gained for optimality. In column a, the process is first at the optimal state (level 3 of OI). However the margin plot shows it is vibrating near the

Figure 16. PP production process example: SI/OI levels of testing data.

limit between level 3 and 1 and finally surpasses the boundary on sample 742. At the same time, the process operating performance falls directly to level 1. In column b, the process is initially at OI level 3 and it begins to degrade after sample 2280 toward level 2. Finally the operating optimality condition falls to level 2 on sample 2363. The two cases are both deterioration in optimality, but the later situation in column b is milder than that of column a since it only falls to level 2. The proposed strategy provides not only online assessment of the current operating performance, but also predictions of the performance trend, and a preliminary advice for improving process operating performance is available. It is suggested that the process is manipulated to higher levels (safety or optimality). While in D space, the D vector is supposed to be pulled away from the lower level region and meanwhile be driven toward the region with higher levels. Since there is a tight relation

10922


are available based on the online assessment, and preliminary suggestions for operating performance improvement are also available. Acknowledgment This research was supported by the National Science Foundation of China under Grant No. 60574047, the Research Fund for the Doctoral Program of Higher Education in China under Grant No. 20050335018, and the National High Technology Research and Development Program of China under Grant No. 2007AA04Z168. Literature Cited

Figure 17. PP production process example: SI levels and corresponding margins.

Figure 18. PP production process example: OI levels and corresponding margins.

between D vector and process variables linked by the Mahalanobis distance formulation, the corresponding regulation of process variables is included implicitly. 5. Conclusion A framework for online assessment of process operating performance based on safety and optimality indices is formulated for multimode operating conditions in complex industrial processes. Gaussian mixture model is constructed on the basis of historical training data with the parameters estimated by EM algorithm. Next, operating performance assessments are formulated on the basis of safety and optimality indices, which are constructed by two successive nonlinear mappings. Further, the safety and optimality indices are divided into several hierarchical levels separated by hyperplanes in D space. Finally, an angular margin is proposed to give specific measure on the relative position on those levels. The proposed approach is applied in two examples of TEP and polypropylene production process. It is concluded from the examples that the proposed framework, based on safety and optimality indices, is able to give an efficient assessment of operating conditions for an industrial process by hierarchical level description and margin analysis. Performance predictions

(1) Crowl, D. A.; Louvar, J. F. Chemical Process Safety: Fundamentals with Applications; 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2002. (2) Venkatasubramanian, V.; Zhao, J.; Viswanathan, S. Intelligent Systems for HAZOP Analysis of Complex Process Plants. Comput. Chem. Eng. 2000, 24, 2291. (3) Rouvroye, J. L.; van den Bliek, E. G. Comparing Safety Analysis Techniques. Reliab. Eng. Syst. Safety 2002, 75, 289. (4) Khan, F. I.; Husain, T.; Abbasi, S. A. Design and Evaluation of Safety Measures Using a Newly Proposed Methodology “SCAP”. J. Loss PreV. Process Ind. 2002, 15, 129. (5) Arena, U.; Romeo, E.; Mastellone, M. L. Recursive Operability Analysis of a Pilot Plant Gasifier. J. Loss PreV. Process Ind. 2008, 21, 50. (6) Lisnianski, A.; Levitin, G. Multi-State System Reliability: Assessment, Optimization, and Application; World Scientific: Singapore, 2003. (7) Meel, A.; Seider, W. D. Plant-Specific Dynamic Failure Assessment Using Bayesian Theory. Chem. Eng. Sci. 2006, 61, 7036. (8) Meel, A.; Seider, W. D. Real-Time Risk Analysis of Safety Systems. Comput. Chem. Eng. 2008, 32, 827. (9) Jiang, H.; Gao, J.; Gao, Z.; Li, G. In Safety Analysis of Process Industry System Based on Complex Networks Theory; International Conference on Mechatronics and Automation, IEEE: Harbin, China, 2007; p 480. (10) Zhou, J.; Chen, G.; Chen, Q. Real-Time Data-Based Risk Assessment for Hazard Installations Storing Flammable Gas. Process Safety Prog. 2008, 27, 205. (11) Venkatasubramanian, V.; Rengaswamy, R.; Yin, K.; Kavuri, S. N. A Review of Process Fault Detection and Diagnosis Part 1: Quantitative Model-Based Methods. Comput. Chem. Eng. 2003, 27, 293. (12) Venkatasubramanian, V.; Rengaswamy, R.; Yin, K.; Kavuri, S. N. A Review of Process Fault Detection and Diagnosis Part 2: Qualitative Models and Search Strategies. Comput. Chem. Eng. 2003, 27, 313. (13) Venkatasubramanian, V.; Rengaswamy, R.; Yin, K.; Kavuri, S. N. A Review of Process Fault Detection and Diagnosis Part 3: Process History Data-Based Methods. Comput. Chem. Eng. 2003, 27, 327. (14) Russell, E. L.; Chiang, L. H.; Braatz, R. D. Data-DriVen Methods for Fault Detection and Diagnosis in Chemical Processes; Springer-Verlag: London, 2000. (15) Kourti, T. Application of Latent Variable Methods to Process Control and Multivariate Statistical Process Control in Industry. Int. J. AdaptiVe Control Signal Process. 2005, 19, 213. (16) Biegler, L. T.; Grossmann, I. E. Retrospective on Optimization. Comput. Chem. Eng. 2004, 28, 1169. (17) Edgar, T. F.; Himmelblau, D. M.; Lasdon, L. S. Optimization of Chemical Processes; McGraw-Hill: New York, 2001. (18) McAuley, K. B.; MacGregor, J. F. Optimal Grade Transition in a Gas Phase Polyethylene Reactor. AIChE J. 1992, 38, 1564. (19) Tarafder, A.; Rangaiah, G. P.; Ray, A. K. A Study of Finding Many Desirable Solutions in Multiobjective Optimization of Chemical Processes. Comput. Chem. Eng. 2007, 31, 1257. (20) Pontes, K. V.; Maciel, R.; Embirucu, M.; Hartwich, A.; Marquardt, W. Optimal Operating Policies for Tailored Linear Polyethylene Resins Production. AIChE J. 2008, 54, 2346. (21) Li, P.; Wendt, M.; Wozny, G. A Probabilistically Constrained Model Predictive Controller. Automatica 2002, 38, 1171. (22) Li, P.; Wendt, M.; Wozny, G. Optimal Planning for Chemical Engineering Processes Under Uncertain Market Conditions. Chem. Eng. Technol. 2004, 27, 641. (23) Li, P.; Garcia, H. A.; Wozny, G. Chance Constrained Programming Approach to Process Optimization under Uncertainty. Comput. Chem. Eng. 2008, 32, 25.

Ind. Eng. Chem. Res., Vol. 48, No. 24, 2009 (24) Hand, D.; Mannila, H.; Smyth, P. Principles of Data Mining; MIT Press: Cambridge, MA, 2001. (25) Paalanen, P.; Kamarainen, J. K.; Ilonen, J.; Kalviainen, H. Feature Representation and Discrimination Based on Gaussian Mixture Model Probability DensitiessPractices and Algorithms. Pattern Recognit. 2006, 39, 1346. (26) Choi, S. W.; Park, J. H.; Lee, I. Process Monitoring Using a Gaussian Mixture Model via Principal Component Analysis and Discriminant Analysis. Comput. Chem. Eng. 2004, 28, 1377–1387. (27) Thissen, U.; Swierenga, H.; Weijer, A. P. d.; Wehrens, R.; Melssen, W. J.; Buydens, L. M. C. Multivariate Statistical Process Control Using Mixture Modelling. J. Chemom. 2005, 19, 23. (28) Yu, J.; Qin, S. J. Multimode Process Monitoring with Bayesian Inference-Based Finite Gaussian Mixture Models. AIChE J. 2008, 54, 1811. (29) Bishop, C. M. Neural Networks for Pattern Recognition; Oxford University Press: New York, 1995. (30) Cao, S.; Rhinehart, R. R. An Efficient Method for on-Line Identification of Steady State. J. Process Control 1995, 5, 363. (31) Jiang, T.; Chen, B.; He, X.; Stuart, P. Application of Steady-State Detection Method Based on Wavelet Transform. Comput. Chem. Eng. 2003, 27, 569. (32) Goodwin, G. C.; Payne, R. L. Dynamic System Identification: Experiment Design and Data Analysis; Academic Press: New York, 1977.

10923

(33) Zhao, Y.; Zhuang, X.; Ting, S.-J. Gaussian Mixture Density Modeling of non-Gaussian Source for Autoregressive Process. IEEE Trans. Signal Process. 1995, 43, 894. (34) Gao, H. Statistical Computation; Peking University Press: Beijing, 1995. (35) Agarwal, N.; Huang, B.; Tamay, E. C. Assessing Model Prediction Control (MPC) Performance. 1. Probabilistic Approach for Constraint Analysis. Ind. Eng. Chem. Res. 2007, 46, 8101. (36) Duda, R. O.; E.Hart, P.; Stork, D. G. Pattern Classification, 2nd ed.; John Wiley & Sons, Inc: New York, 2001. (37) Downs, J. J.; Vogel, E. F. A Plant-Wide Industrial Process Control Problem. Comput. Chem. Eng. 1993, 17, 245. (38) McAvoy, T. J.; Ye, N. Base Control for the Tennessee Eastman Problem. Comput. Chem. Eng. 1994, 18, 383. (39) Lyman, P. R.; Georgakis, C. Plant-Wide Control of the Tennessee Eastman Problem. Comput. Chem. Eng. 1995, 19, 321. (40) Ricker, N. L. Decentralized Control of the Tennessee Eastman Challenge Process. J. Process Control 1996, 6, 205.

ReceiVed for reView December 4, 2008 ReVised manuscript receiVed September 28, 2009 Accepted September 28, 2009 IE801870G

Online Probabilistic Assessment of Operating Performance Based on

Recommend Documents