Reconstruction-Based Multivariate Process Fault Isolation Using

of statistical fault isolation to the variable selection in regression analysis and solves it in a Bayesian framework. ... 3. 1. Introduction. With th...
1 downloads 8 Views 950KB Size
Subscriber access provided by Kaohsiung Medical University

Article

Reconstruction-Based Multivariate Process Fault Isolation Using Bayesian Lasso Zhengbing Yan, Yuan Yao, Tsai-Bang Huang, and Yi-Sern Wong Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b05189 • Publication Date (Web): 02 Mar 2018 Downloaded from http://pubs.acs.org on March 4, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Reconstruction-Based Multivariate Process Fault Isolation Using Bayesian Lasso Zhengbing Yan a, Yuan Yao b*, Tsai-Bang Huang c, Yi-Sern Wong c a

College of Physics and Electronic Information Engineering, Wenzhou University, Wenzhou 325035, China

b

Department of Chemical Engineering, National Tsing Hua University, Hsinchu 31003, Taiwan c

Kaohsiung Factory, Chang Chun Plastics Co., Ltd., No.14 Gongye 1st Rd., Renwu District, Kaohsiung 81469, Taiwan

Abstract: To ensure process safety and product quality, multivariate statistical process monitoring (MSPM) has been widely used in industry for decades. Among the steps of MSPM, fault isolation is an important link between the fault detection and the root-cause diagnosis, which identifies the variables closely related to the detected process abnormality. The existing methods for fault isolation often suffer from the smearing effect or rely on the impractical requirement of a sufficient amount of historical fault data. To solve these problems, a reconstruction method based on the Bayesian Lasso (short for least absolute shrinkage and selection operator) is proposed in this work, which transforms the problem of statistical fault isolation to the variable selection in regression analysis and solves it in a Bayesian framework. If the posterior distribution of a Lasso coefficient changes significantly before and after the occurrence of the fault, the corresponding process 1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 41

variable has a large probability to be faulty. Furthermore, the Bayesian framework permits the tracking of fault propagation, thereby facilitating the subsequent root-cause diagnosis. The feasibility of the proposed paper is illustrated by case studies.

Keywords: fault isolation, multivariate statistical process monitoring, Bayesian Lasso, variable selection, fault reconstruction.

*

Corresponding

author:

Tel:

886-3-5713690,

Fax:

[email protected]

2

ACS Paragon Plus Environment

886-3-5715408,

Email:

Page 3 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

1. Introduction With the development of sensor and computer technologies, large amounts of data are collected in industrial processes, providing opportunities for data-driven process monitoring. In recent decades, multivariate statistical process monitoring (MSPM) approaches have been widely studied in the literature and implemented in different branches of industry 1-4. The entire procedure of MSPM comprises of four major stages, including fault detection, fault isolation, root-cause diagnosis, and process recovery 5. This work focuses on fault isolation which aims to identify the critical variables closely correlated to the process abnormality detected in the first stage. In the literature, this task is sometimes called “fault identification” or “fault diagnosis”. In order to avoid confusion with the concept of root-cause diagnosis, the term “fault isolation” is used in this article. In previous research of MSPM, the problem of fault detection has been extensively studied. However, achievements in fault isolation are relatively limited, although fault isolation is an important link between fault detection and root-cause diagnosis. The main reason is that it is difficult to analyze the simultaneous impact of multiple variables on monitoring statistics. Traditionally, contribution plots are adopted as standard tools for identifying faulty variables, which calculate the contribution of each variable to a multivariate statistic, such as T2 or squared prediction error (SPE), and compare it to a predetermined confidence limit

6, 7

. The logic of contribution plots seems easy to

understand but not necessarily correct. As pointed out in 8, the distribution of the contributions of non-faulty variables may change before and after a process abnormality occurs. Therefore, it is improper to conclude which variables are faulty based on the control limits derived from the data collected under the normal operating conditions

3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 41

(NOC). This is the reason why the “smearing” effect, i.e. the influence of faulty variables on the contributions of non-faulty variables, is often observed in contribution plots 9 and causes misleading results. In recent years, researchers attempted to resolve the smearing problem in fault isolation by adopting a new philosophy named reconstruction reconstruction is missing variable analysis

12

10, 11

. The core of fault

. The variable(s) along each candidate fault

direction are treated as if they were missing and the monitoring statistic is re-calculated. After repeating the procedure for all the directions, the variable(s) reducing the monitoring statistic most significantly are considered to be critical to the fault. However, the effectiveness of reconstruction relies on an assumption that the complete information of candidate fault directions is known or sufficient amounts of historical fault data are available. Unfortunately, this assumption is not easy to satisfy in practical cases. In order to relieve the dependence on prior knowledge or historical data, the reconstruction-based contribution (RBC) approach was proposed

13

. Nevertheless, RBC only guarantees

correct isolation in the situation of single sensor failures with large fault magnitudes 14. In handling multivariate faults, the smearing effect can be observed when the confidence limits are implemented in the RBC plots 15. In later research, the Bayesian decision was adopted as a post-processing step of RBC to avoid smearing

14, 15

. Nevertheless, the

performance of these methods depends on the parameters selected for the conditional probability distributions. Most recently, the concept of variable selection has attracted increasing attention in MSPM community reconstruction

20-22

16-19

and been integrated with fault

. The methods belonging to this type formulate the problem of fault

isolation in the form of L0 or L1 regularization and compute the solutions using state-of-

4

ACS Paragon Plus Environment

Page 5 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

the-art solvers. The variables corresponding to non-zero regression coefficients are believed to be faulty. The L1 regularization (also known as Lasso (short for least absolute shrinkage and selection operator)) is preferred because it is more computationally efficient. Most recently, the similar idea has been extended to fault isolation of nonstationary processes 23. Despite the good performance of the variable selection-based reconstruction method in many cases, Lasso only provides a point estimate and does not explore the effect of the estimation uncertainty. Furthermore, the isolation results may be affected by the choice of the tuning parameter. In this research work, a Bayesian Lasso 24 framework is adopted to solve these problems. If the posterior distribution of a Bayesian Lasso regression coefficient changes significantly before and after the occurrence of the fault, the corresponding process variable is likely to be faulty. Furthermore, the Bayesian framework facilitates the tracking of the fault propagation paths.

2. Problem formulation Let X denote the process data matrix with the dimensions of n × p, which contains the historical process data collected under NOC. Here n is the number of observations and p is the number of variables. Without loss of generality, suppose that X is normalized by auto-scaling, such that all variables have zero mean and unit variance. Such a preprocessing step is required by most MSPM methods to eliminate the influence of units of measurement. As discussed in

11

, a family of monitoring statistics can be unified

described as follows:

φ = xT PA+ PT x .

(1)

5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 41

In the above equation, the vector x represents a process observation and P is a matrix containing the model information. For process monitoring based on principal component analysis (PCA) or partial least squares (PLS)

25

, P is the loading matrix calculated by

decomposing X. For other multivariate methods, such as independent component analysis (ICA)

26

, the monitoring statistics can also be represented in a similar way by slightly

modifying the equation. The diagonal matrix A with the dimensions of p × p has different elements for different statistics. For example, in the PCA-based monitoring scheme, a combined index 11 integrating the information of T2 and SPE can be defined by choosing A = diag(λ1 χ l2 L λl χ l2

δ 2 L δ 2) ,

(2)

where l is the number of principal components (PCs), λ1, …, λl are the eigenvalues of the sample covariance matrix of X sorted in descending order, χ l2 and δ2 represent the control limits for T2 and SPE, respectively. In online fault detection, the monitoring statistic is calculated using the formula in (1) and compared to a predetermined threshold. If the threshold is exceeded, then the process is considered as out of statistical control and fault isolation should be conducted to identify the variables significantly contributing to the fault. The previous research

22

shows that

the task of reconstruction-based isolation can be formulated as a quadratic optimization problem: min( x − f )T PA + P T (x − f ) ,

(3)

f

where f is a vector consisting of the fault magnitudes in the directions of different process variables such that (x – f) represents the normal values obtained by eliminating the influence of the process abnormality. The core idea of this method is that the underlying

6

ACS Paragon Plus Environment

Page 7 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

fault

direction

should

minimize

the

reconstructed

monitoring

statistic

(x − f )T PA+ PT (x − f ) .

3. Bayesian Lasso-based reconstruction for fault isolation If a complete candidate set of the potential fault directions is available, the solution of the optimization problem described in (3) can be easily found by exhaustive search. In practice, however, such information is usually incomplete due to the process complexity and the diversity of fault types. In such cases, we only get the trivial solution, i.e. f = x. In the context of reconstruction analysis, this means that all process variables are regarded to be abnormal, which is unreasonable in most circumstances. The practical experience indicates that it is unlikely to observe the simultaneous shift of all the process variables when a fault occurs. Instead, a particular type of fault is usually dominated by a subset of variables 16. In other words, the set of faulty variables is often sparse. Guided by the engineering practices, the cost function in (3) can be modified by adding an L1 penalty term. In doing this, the problem of reconstruction is reformulated:

min((x − f )T PA + PT (x − f ) + θ f 1 ) ,

(4)

f

where the L1 norm introduces a shrinkage effect biasing the estimate of f and eventually leads to sparsity, and θ is the tuning factor. This formulation can be further converted to a standard form of Lasso regression 27 min(( y − Zf ))T ( y − Zf )) + θ f 1 )

(5)

f

7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 41

by conducting a Cholesky decomposition on A+, where A + = LLT , y = (PL)T x and

Z = (PL)T . In

22

, the problem in (5) is solved using the algorithm of least angle

regression (LARS) 28. Nevertheless, the previous study only reports point estimates of f, while the estimation uncertainty has not been well explored. To address this issue, the Bayesian Lasso estimates are adopted in this work.

3.1. Gibbs sampler for Bayesian Lasso-based reconstruction As pointed out in

27

, the Lasso estimates of the regression coefficients are equivalent to

the posterior mode estimates with independent and identical Laplace priors. Here, the conditional Laplace prior proposed in 24 is adopted: p

π (f | σ 2 ) = ∏ j =1

θ 2 σ2

e

−θ f j / σ 2

,

(6)

together with the non-informative scale-invariant marginal prior π (σ 2 ) = 1

σ2

on the

variance of the prediction error of the regression model, where fj (j = 1,…,p) is the jth entry in the vector f, i.e. the fault magnitude associated with the jth variable. As known, the Laplace distribution can be represented as a scale mixture of Gaussians 29: ∞ a −a z 1 − z 2 /(2 s ) a 2 − a2 s /2 e =∫ e e ds, a > 0 . 0 2 2 2π s

(7)

Accordingly, the full model is expressed in a hierarchical form: y | Z, f , σ 2 ~ Gaussian( Zf , σ 2 1), f | σ 2 ,τ 12 ,...,τ 2p ~ Gaussian(0, σ 2 D), p

σ 2 ,τ 12 ,...,τ p2 ~ π (σ 2 )dσ 2 ∏ j =1

θ2 2

e

−θ 2τ 2j / 2

(8) dτ 2j ,

8

ACS Paragon Plus Environment

Page 9 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

where D = diag(τ 12 ,...,τ 2p ) and σ 2 ,τ 12 ,...,τ 2p > 0 . A Gibbs sampler can be constructed to T

fit the model based on the full conditional distributions of f, σ2, and τ 2 = τ 12 L τ p2  . In online applications, the posterior distribution of f calculated for one process observation can be used as the prior for the next, so as to facilitate the tracking of fault propagation paths. The detailed steps are as follows. Step 0 (initialization). Choose the number of the Gibbs sampling iterations K. Usually, K is a large positive integer. Specify the values for r and δ which are the shape and scale parameters of the gamma prior on θ2.

π (θ 2 ) =

δr Γ( r )

(θ 2 )r −1 e −δθ , θ 2 > 0 (r > 0, δ > 0) . 2

(9)

Let m = 1 and k = 0. Give the initial values of σ2 and τ2 which are denoted as σ 12 and τ12 . (0)

(0)

The subscript under the main symbol is the index of iterations, while the right-subscript is the index of observations. It is noted that the selection of the initial parameters is not critical to the final results but affects the convergence rate of the algorithm. Step 1. Generate a random vector f m

from the following multivariate Gaussian

( k +1)

distribution:

fm ~ Gaussian(Bm−1 ZT y m , σ m2 Bm−1 ) ,

( k +1)

(k )

(k )

(10)

k)

where, Bm = ZT Z + Dm−1 , D −m1 = diag ( τ m2 ) , y m = ( PL )T x m , and xm is the mth process (k )

(k )

(k )

(k )

observation. σ 2 from the following inverse-gamma distribution: Step 2. Generate a random value ( m,k +1)

9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

σ m2 ~ Inv-Gamma(

( k +1)

p −1 p + , Cm ) , 2 2 ( k +1)

Page 10 of 41

(11)

where C m = ( y m − Z f m )T ( y m − Z f m ) / 2 + f mT D m−1 f m / 2 . ( k +1)

( k +1)

( k +1)

( k +1) ( k ) ( k +1)

Step 3. Generate a random value θ m2 from the following Gamma distribution: ( k +1)

θ m2 ~ Gamma( p + r , ∑ j =1 p

( k +1)

τ m2 , j +δ ) ,

(k )

2

(12)

where τ m2 , j is the jth entry of τ2m . Step 4. Generate a random vector τ 2m from the following inverse Gaussian distribution: ( k +1)

1

τ m2 , j ( k +1)

θ m2 σ m2

~ Inv-Gaussian(

( k +1) ( k +1) 2 m 2 ( k +1) m, j ( k +1)

f

, θ ),

(13)

where f m , j is the jth entry of fm. Step 5. Update k = k +1. Go to Step 1 until k = K + 1. Step 6. Collect the vectors f m computed in the last k% iterations to generate a sample (i )

from the posterior distribution of the fault magnitudes of the mth process observation, i.e. fm, where i = K − k% + 1,K, K . When the sample size is large enough, the probability

density function of each fj (j = 1,…,p) can be approximated by the histogram. Step 7. Update m = m + 1, σ m2 = σ m2 −1 , and τ 2m = τ m2 −1 . Go to Step 1 until m = M + 1, where (0)

(K )

(0)

(K )

M is the total number of observations to be reconstructed.

10

ACS Paragon Plus Environment

Page 11 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

3.2. Distribution-based fault isolation The next step is to identify the faulty variables from the posterior density for each fj. A natural idea is to perform the hypothesis testing on each fj to check whether the null hypothesis, H0: fj = 0, can be rejected. If not, then the corresponding variable is considered to be critical to the detected fault. However, according to our practice and experience, the hypothesis testing-based criterion is usually too conservative and tends to accept H0 in a large number of cases. As a result, missing alarms occur frequently. In order to avoid such a problem, a distribution-based fault isolation method is proposed in this research work, the basic idea of which is as follows. The Bayesian Lasso procedure can be applied to each observation recorded during process operation. When the sample is collected under NOC, each variable behaves as expected and the corresponding Lasso estimate of fj (j = 1,…,p) follows an approximately symmetry distribution with the mode at 0. When a fault occurs, the distribution of each fj changes, while those corresponding to the faulty variables become asymmetry and spread widely, with the mode deviated from 0, and change much more significantly than those corresponding to the non-faulty variables. Based on the above findings, the proposed isolation procedure consists of two stages. In the first stage, the distribution of the Lasso estimate of fj is calculated for each online collected process observation and compared to those obtained based on the normal samples. The similarities between the distributions are quantified with the β statistic, while a small value of β implies that two distributions are dissimilar. Then in the second stage, the β values corresponding to different process variables are divided into two groups through clustering analysis. The group with smaller β values indicates the set of faulty variables. The detailed procedure is described below.

11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 41

Step 1. Given M process observations that are recorded under NOC, the Bayesian Lasso estimate of fm is calculated by using the algorithm described in Section 3.1, where m = 1,…,M is the index of observation, and fm = [fm1 … fmp]T is the vector of fault magnitudes corresponding to the mth observation in the historical dataset. Step 2. Identify the (1 - α) confidence interval of fmj (m = 1,…,M, j = 1,…,p). Denote the upper and lower confidence bounds as Umj and Lmj, respectively. Step 3. For each fault sample, calculate the Bayesian Lasso estimate of f by using the algorithm described in Section 3.1, where f = [f1 … fp]T. Step 4. Set m = 1 and j = 1. Step 5. Calculate βmj = Prob(Lmj < fj < Umj). Step 6. Set m = m + 1 and go to Step 5, until m = M. Step 7. Set m = 1, j = j + 1 and go to Step 5, until j = p. Step 8. Implement clustering analysis on the set {β1 … βp} to divide the entire set into two groups, where βj = [β1j … βMj]T (j = 1,…,p). Denote the cluster centers as c1 and c2, respectively. In the following of this paper, the K-means clustering

30

is adopted for

illustration. Step 9. Compare ||c1||2 and ||c2||2. Choose the variables corresponding to the group with a center having a larger L2 norm as the isolated faulty variables.

4. Case studies

In this section, the Tennessee Eastman (TE) process

31

, which is a well-known

benchmarking process in the field of MSPM, is used to illustrate the feasibility of the proposed fault isolation method. This process consists of five main units, including a

12

ACS Paragon Plus Environment

Page 13 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

reactor, a condenser, a separator, a stripper, and a compressor. Four reactants (A, C, D and E) and an inert tracer (B) are fed to the reactor to produce two products (G and H) and a byproduct (F). Totally 52 process variables, including 11 manipulated variables and 41 measured variables, are recorded in both normal operation and abnormal situations triggered by 21 different types of faults. The sampling interval is 3 min. In each scenario, the fault occurs at the 161st sampling time point. The detailed description of the dataset can be found in the literature

32

. The flowchart of the TE process is shown in Fig. 1. In

the following of this paper, the process variables are denoted as xj, where j is the variable index. A complete list of the variables can be found in the literature 31. For fault detection, a PCA model with 9 retained PCs was constructed using the normal operation data, while the combined index was selected as the monitoring statistic. In online monitoring, the statistic was calculated for each process observation and compared to the predetermined 95% control limit. When the statistic was beyond the limit, the process status was considered out of statistical control and the proposed fault isolation approach was implemented to identify the variables critical to the detected fault. In order to obtain the Bayesian Lasso estimates of the fault magnitudes, the Gibbs sampler was adopted as described in Section 3.1, where the parameters were set as follows. K = 2000,

k% = 1000, r = 1.5, δ = 1, σ 12 = 0.1, and τ12 = [0.1 … 0.1]T. (0)

(0)

4.1. Scenario 1: Fault 1 The first case is about Fault 1 described in the original literature 31, where the A/C feed ratio in Stream 4 was changed. As a result, the amount of A in the recycle stream

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

decreased, which further led to a decrease of the A composition in Stream 6. The feedback control system attempted to compensate this disturbance by increasing x44 which is the valve opening related to the A feed flow in Stream 1. As a sequence, the flow rate of Stream 1, x1, increased, which raised the level in the reactor and hence affected the status of multiple variables. After a period of time, the system settled down. Most variables returned to their set-points or achieved a new steady state closed to the original one, while the magnitudes of x1 and x44 changed significantly because of the control actions. A number of variable trajectories are plotted in Fig. 2. For each observation, the Bayesian Lasso estimate of the fault magnitudes was obtained by using the Gibbs sampler. For illustration, the normalized histograms, which represent the distributions of f, calculated for the 160th and 230th observations are shown in Fig. 3 and Fig. 4, respectively. Please note the different scales in the x-axes. It can be observed that when the process is under NOC, the distribution of each fj is relatively symmetric and with a small standard deviation. After the fault occurs, the distributions of the fj corresponding to the faulty variables become significantly asymmetric and change much more significantly than those corresponding to the normal variables. The β statistic was used to quantify the changes of the distributions, based on which clustering analysis was conducted to give the isolation results as shown in Fig. 5. The results well reflect the mechanism of this fault and are consistent with the results obtained by using the combination of RBC and Bayesian decision 14. An advantage of the proposed method is that the experience of selecting the distribution parameters is not required. For comparison, the conventional RBC was also calculated. The results are shown in Fig. 6, which are very difficult to interpret because of the smearing effects.

14

ACS Paragon Plus Environment

Page 14 of 41

Page 15 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

4.2. Scenario 2: Fault 2 The second illustrative case is Fault 2. The cause of this fault was a step change in the composition of the inert B in Stream 4. The A/C ratio did not change accordingly. Because of the material and energy balances, the downstream reactions and operations were affected by the disturbance in the feedstock. For example, the B composition in Stream 9 changed, which triggered the control action. As a result, the opening of the purge value, x47, was adjusted. Consequently, the purge rate x10 was also increased. According to the discussions in the literature

33

, x1, x28, x34, and x47 are the major faulty

variables. The trajectories of these four variables are shown in Fig. 7. The results of Gibbs sampling are represented by Fig. 8 and Fig. 9. Significant changes in the distributions of fj corresponding to the faulty variables can be observed from the comparison between the estimates made for the normal and fault samples. After quantifying such changes using the β statistic, clustering analysis was implemented to identify the variables most critical to the fault. The results are shown in Fig. 10, which are consistent with the process understanding. Similar results can be achieved by using the Bayesian decision technique

14

. However, as discussed before, this method requires

more experience in the parameter selection. The RBC isolation results are plotted in Fig. 11 for comparison. It is difficult to get an indication of the cause of the fault based on such a messy plot.

15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4.3. Scenario 3: Fault 4 The third example is about Fault 4, a step change in the inlet temperature of the reactor cooling water. This abnormality directly affected the reactor temperature x9. It had not been long before the feedback controller actioned. The flow rate of the reactor cooling water, x51, was adjusted to bring x9 back to its set-point. Because of the good controller performance, the adjustment was efficient and x9 returned to normal in a short time. As a result, the fault did not propagate to other variables. The only variable reaching a new steady state is x51. The trajectories of four process variables are plotted in Fig. 12. Fig. 13 shows the estimate of the posterior distribution of f at the 160th sampling time point which is under NOC, while the estimate corresponding to a fault sample is shown in Fig. 14. Obvious changes in the distribution of f9 and f51 can be observed by comparing these two figures. After calculating the β statistic and conducting clustering analysis, the Bayesian Lasso-based fault isolation results were obtained. As shown in Fig. 15, x51 is the only process variable contributing consistently to this fault, while x9 is only significant at the first time point after detecting the fault. Fig. 16 plots the RBC results which are again polluted by smearing. The comparison results of the computation time of RBC and that of the proposed method are shown in Table 1. It is noted that the proposed Bayesian-Lasso isolation method takes more time in the computation. However, the computation time is still acceptable for industrial applications. To obtain the results shown in Table 1, the programs were run ten times and the mean values of the calculation time were recorded.

16

ACS Paragon Plus Environment

Page 16 of 41

Page 17 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

5. Conclusions

Fault isolation is a critical step in MSPM, identifying the variables contributing most to the detected fault and provides information for the following root-cause diagnosis step. In this work, a Bayesian Lasso-based reconstruction method is proposed, which better considers the estimation uncertainty and has the advantages of both variable selection and Bayesian inference. Comparing to the existing methods, the proposed method well solves the problem of smearing effects and does not require much experience in the determination of the prior conditional probabilities. Such characteristics make this method a useful tool in industrial applications. In addition, most statistical fault isolation methods may suffer from the measurement noise which affects the signal-to-noise ratio contained in the data and influences the accuracy of the isolation results. The method proposed in this paper adopts a Bayesian inference framework which can integrate the information contained in multiple samples and provide more robust results. The proposed method may have a potential limitation which is inherited from Lasso. Due to the nature of the L1 penalty, the performance of Lasso degrades when the predictors are highly correlated. In the context of multivariate fault isolation, the Lasso (or Bayesian Lasso) based method may not find all the faulty variables when strong cross-correlation exists among them. In such situations, the Bayesian elastic net

34

may be a good

alternative. This issue will be addressed in future research.

Acknowledgement

This work was supported in part by Ministry of Science and Technology, ROC, under Grant no. MOST 106-2622-8-007-017. Yan was supported by the National Natural

17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Science Foundation of China (No. 61703309) and Zhejiang Provincial Natural Science Foundation of China under Grant No. LY18F030014.

References

1.

Yao, Y.; Gao, F., A survey on multistage/multiphase statistical modeling methods

for batch processes. Annual Reviews in Control 2009, 33, (2), 172-183. 2.

Qin, S. J., Survey on data-driven industrial process monitoring and diagnosis.

Annual Reviews in Control 2012, 36, (2), 220-234. 3.

Yin, S.; Ding, S. X.; Xie, X.; Luo, H., A review on basic data-driven approaches

for industrial process monitoring. IEEE Transactions on Industrial Electronics 2014, 61, (11), 6418-6428. 4.

Ge, Z., Review on data-driven modeling and monitoring for plant-wide industrial

processes. Chemometrics and Intelligent Laboratory Systems 2017, 171, 16-25. 5.

Chiang, L. H.; Kotanchek, M. E.; Kordon, A. K., Fault diagnosis based on Fisher

discriminant analysis and support vector machines. Computers & Chemical Engineering 2004, 28, (8), 1389-1401.

6.

Westerhuis, J.; Gurden, S.; Smilde, A., Generalized contribution plots in

multivariate statistical process monitoring. Chemometrics and Intelligent Laboratory Systems 2000, 51, (1), 95-114. 7.

Conlin, A.; Martin, E.; Morris, A., Confidence limits for contribution plots.

Journal of Chemometrics 2000, 14, 725-736. 8.

Liu, J.; Wong, D. S. H.; Chen, D.-S., Bayesian filtering of the smearing effect:

Fault isolation in chemical process monitoring. Journal of Process Control 2014, 24, (3), 1-21. 18

ACS Paragon Plus Environment

Page 18 of 41

Page 19 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

9.

Van den Kerkhof, P.; Vanlaer, J.; Gins, G.; Van Impe, J. F. M., Analysis of

smearing-out in contribution plot based fault isolation for statistical process control. Chemical Engineering Science 2013, 104, 285-293. 10.

Dunia, R.; Qin, S., Subspace approach to multidimensional fault identification

and reconstruction. AIChE Journal 1998, 44, (8), 1813-1831. 11.

Yue, H. H.; Qin, S. J., Reconstruction-based fault identification using a combined

index. Industrial & Engineering Chemistry Research 2001, 40, (20), 4403-4414. 12.

Kariwala, V.; Odiowei, P. E.; Cao, Y.; Chen, T., A branch and bound method for

isolation of faulty variables through missing variable analysis. Journal of Process Control 2010, 20, (10), 1198-1206. 13.

Alcala, C. F.; Qin, S. J., Reconstruction-based contribution for process monitoring.

Automatica 2009, 45, (7), 1593-1600. 14.

Zheng, Y.; Mao, S.; Liu, S.; Wong, D. S. H.; Wang, Y. W., Normalized relative

RBC-based minimum risk Bayesian decision approach for fault diagnosis of industrial process. IEEE Transactions on Industrial Electronics 2016, 63, (12), 7723-7732. 15.

Liu, J., Fault diagnosis using contribution plots without smearing effect on non-

faulty variables. Journal of Process Control 2012, 22, (9), 1609-1623. 16.

Wang, K.; Jiang, W., High-dimensional process monitoring and fault isolation via

variable selection. Journal of Quality Technology 2009, 41, (3), 247-258. 17.

Zou, C.; Jiang, W.; Tsung, F., A LASSO-based diagnostic framework for

multivariate statistical process control. Technometrics 2011, 53, (3), 297-309. 18.

Kuang, T.-H.; Yan, Z.; Yao, Y., Multivariate fault isolation via variable selection

in discriminant analysis. Journal of Process Control 2015, 35, 30-40.

19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

19.

Yan, Z.; Kuang, T.-H.; Yao, Y., Multivariate fault isolation of batch processes via

variable selection in partial least squares discriminant analysis. ISA Transactions 2017, 70, 389-399. 20.

He, B.; Yang, X.; Chen, T.; Zhang, J., Reconstruction-based multivariate

contribution analysis for fault isolation: A branch and bound approach. Journal of Process Control 2012, 22, (7), 1228-1236. 21.

He, B.; Zhang, J.; Chen, T.; Yang, X., Penalized reconstruction-based

multivariate contribution analysis for fault isolation. Industrial & Engineering Chemistry Research 2013, 52, (23), 7784-7794. 22.

Yan, Z.; Yao, Y., Variable selection method for fault isolation using least absolute

shrinkage and selection operator (LASSO). Chemometrics and Intelligent Laboratory Systems 2015, 146, 136-146. 23.

Sun, H.; Zhang, S.; Zhao, C.; Gao, F., A sparse reconstruction strategy for online

fault diagnosis in nonstationary processes with no priori fault information. Industrial & Engineering Chemistry Research 2017, 56, (24), 6993-7008. 24.

Park, T.; Casella, G., The Bayesian Lasso. Journal of the American Statistical

Association 2008, 103, (482), 681-686. 25.

Wise, B.; Gallagher, N., The process chemometrics approach to process

monitoring and fault detection. Journal of Process Control 1996, 6, (6), 329-348. 26.

Lee, J. M.; Yoo, C. K.; Lee, I. B., Statistical process monitoring with independent

component analysis. Journal of Process Control 2004, 14, (5), 467-485. 27.

Tibshirani, R., Regression shrinkage and selection via the lasso. Journal of the

Royal Statistical Society. Series B (Methodological) 1996, 58, (1), 267-288.

20

ACS Paragon Plus Environment

Page 20 of 41

Page 21 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

28.

Efron, B.; Hastie, T.; Johnstone, I.; Tibshirani, R., Least angle regression. The

Annals of Statistics 2004, 32, (2), 407-451. 29.

Andrews, D. F.; Mallows, C. L., Scale Mixtures of Normal Distributions. Journal

of the Royal Statistical Society. Series B (Methodological) 1974, 36, (1), 99-102. 30.

Jain, A. K., Data clustering: 50 years beyond K-means. Pattern Recognition

Letters 2010, 31, (8), 651-666. 31.

Downs, J.; Vogel, E., A plant-wide industrial process control problem. Computers

& Chemical Engineering 1993, 17, (3), 245-255. 32.

Russell, E. L.; Chiang, L. H.; Braatz, R. D., Fault detection in industrial processes

using canonical variate analysis and dynamic principal component analysis. Chemometrics and Intelligent Laboratory Systems 2000, 51, (1), 81-93. 33.

Chiang, L. H.; Pell, R. J., Genetic algorithms combined with discriminant analysis

for key variable identification. Journal of Process Control 2004, 14, (2), 143-155. 34.

Li, Q.; Lin, N., The Bayesian elastic net. Bayesian Analysis 2010, 5, (1), 151-170.

21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure list:

Figure 1. Flowchart of TE process Figure 2. Variable trajectories in the case of Fault 1 Figure 3. Fault magnitudes distribution at the 160th sampling time point in the case of Fault 1 Figure 4. Fault magnitudes distribution at the 230th sampling time point in the case of Fault 1 Figure 5. Fault isolation results of Fault 1: Bayesian Lasso Figure 6. Fault isolation results of Fault 1: RBC Figure 7. Variable trajectories in the case of Fault 2 Figure 8. Fault magnitudes distribution at the 160th sampling time point in the case of Fault 2 Figure 9. Fault magnitudes distribution at the 230th sampling time point in the case of Fault 2 Figure 10. Fault isolation results of Fault 2: Bayesian Lasso Figure 11. Fault isolation results of Fault 2: RBC Figure 12. Variable trajectories in the case of Fault 4 Figure 13. Fault magnitudes distribution at the 160th sampling time point in the case of Fault 4

22

ACS Paragon Plus Environment

Page 22 of 41

Page 23 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 14. Fault magnitudes distribution at the 230th sampling time point in the case of Fault 4 Figure 15. Fault isolation results of Fault 4: Bayesian Lasso Figure 16. Fault isolation results of Fault 4: RBC

23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. Flowchart of TE process

24

ACS Paragon Plus Environment

Page 24 of 41

Page 25 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 2. Variable trajectories in the case of Fault 1

25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3. Fault magnitudes distribution at the 160th sampling time point in the case of Fault 1

26

ACS Paragon Plus Environment

Page 26 of 41

Page 27 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 4. Fault magnitudes distribution at the 230th sampling time point in the case of Fault 1

27

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5. Fault isolation results of Fault 1: Bayesian Lasso

28

ACS Paragon Plus Environment

Page 28 of 41

Page 29 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 6. Fault isolation results of Fault 1: RBC

29

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7. Variable trajectories in the case of Fault 2

30

ACS Paragon Plus Environment

Page 30 of 41

Page 31 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 8. Fault magnitudes distribution at the 160th sampling time point in the case of Fault 2

31

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9. Fault magnitudes distribution at the 230th sampling time point in the case of Fault 2

32

ACS Paragon Plus Environment

Page 32 of 41

Page 33 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 10. Fault isolation results of Fault 2: Bayesian Lasso

33

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 11. Fault isolation results of Fault 2: RBC

34

ACS Paragon Plus Environment

Page 34 of 41

x 51

x9

x 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

x5

Page 35 of 41

Figure 12. Variable trajectories in the case of Fault 4

35

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 13. Fault magnitudes distribution at the 160th sampling time point in the case of Fault 4

36

ACS Paragon Plus Environment

Page 36 of 41

Page 37 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 14. Fault magnitudes distribution at the 230th sampling time point in the case of Fault 4

37

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Zoom in

Figure 15. Fault isolation results of Fault 4: Bayesian Lasso

38

ACS Paragon Plus Environment

Page 38 of 41

Page 39 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 16. Fault isolation results of Fault 4: RBC

39

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 1. Computation time (Intel i5-4590 3.3GHz, 8GB RAM, Win10, Matlab)

Methods RBC Bayesian Lasso (proposed method)

40

ACS Paragon Plus Environment

Time (s) 0.17 6.31

Page 40 of 41

Page 41 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

41

ACS Paragon Plus Environment