Monitoring Batch Processes Using Sparse Parallel Factor

Oct 12, 2017 - Traditional multivariate statistical analysis methods used for batch process monitoring fail to build high-performance monitoring model...
2 downloads 10 Views 2MB Size
Subscriber access provided by LAURENTIAN UNIV

Article

Monitoring batch processes using sparse parallel factor decomposition Lijia Luo, Shiyi Bao, Jianfeng Mao, and Di Tang Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b02618 • Publication Date (Web): 12 Oct 2017 Downloaded from http://pubs.acs.org on October 13, 2017

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 free 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 accessible to all readers and 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 31

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

155x66mm (300 x 300 DPI)

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 31

Monitoring batch processes using sparse parallel factor decomposition Lijia Luo∗, Shiyi Bao, Jianfeng Mao, and Di Tang Institute of Process Equipment and Control Engineering, Zhejiang University of Technology, Hangzhou 310014, China ABSTRACT: Traditional multivariate statistical analysis methods used for batch process monitoring fail to build high-performance monitoring models partly because of the inability to extract sparse and interpretable latent variables from the batch data. To overcome this deficiency, we propose sparse parallel factor (SPARAFAC) decomposition for the sparse modeling of three-way data arrays of batch processes. The SPARAFAC extracts sparse factor matrices from the three-way data array. Sparse factor matrices have the potential advantage of being easier to interpret, because they eliminate

redundant

data

information

and

reveal

significant

variable

correlations.

A

SPARAFAC-based batch process monitoring method is proposed. Due to the sparsity and interpretability, the SPARAFAC monitoring model is well suitable for fault detection, analysis and diagnosis. Two fault detection indices are defined based on the SPARAFAC model, and they are capable of detecting the occurrence of process faults. Multi-level contribution plots are developed for the identification of faulty variable groups and faulty variables responsible for process faults. The effectiveness and advantages of the proposed methods are illustrated by an industrial-scale fed-batch fermentation process. 1. Introduction Batch processes are widely used in chemical industries to produce higher value added chemicals,



Corresponding Author. Tel.: +86 (0571) 88320349. E-mail address: [email protected] 1

ACS Paragon Plus Environment

Page 3 of 31

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

such as special polymers, pharmaceuticals, and biochemicals. These processes have the common trait of a repeated operation procedure starting from nearly the same initial conditions. For the batch operation, it is crucial to guarantee the consistency and reproducibility of batches. Most of the existing industrial approaches to ensure the consistency of batches are based on the precise sequencing and automation of all the stages in the batch operation.1 Batch processes, however, often suffer from the batch-to-batch disturbances, such as variations in raw material properties, changes in the startup initialization, and disturbances encountered during an individual batch execution.2 Consequently, significant variations in the final product quality are often observed between batches.2 Monitoring batch processes is therefore very important to ensure the safe operation and the production of consistent high quality products. Batch process monitoring methods have focused on the use of fundamental mathematical models (e.g. theoretical models of batch processes), detailed knowledge based models (e.g. expert systems or artificial intelligence methods), or data-driven models.3-5 The data-driven monitoring methods use only the information contained in the batch data, rather than utilizing detailed engineering knowledge about the batch process as in model-based or knowledge-based methods. This enables the data-driven methods to be easily implemented in industrial practice and well suited for monitoring complicated batch processes. Data-driven batch process monitoring methods have been widely investigated over the past decades.4-10 These methods often use multivariate statistical analysis (MSA) techniques, such as partial least squares (PLS) and principal component analysis (PCA), to extract important information from batch data sets and to build data-driven monitoring models. The historical data of batches are usually recorded in a three-way data array (with three ways of batches, variables and operating time). 2

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

The MSA techniques available for such a three-way data array are classified into two types: multiway methods and multilinear methods. The multiway methods are equivalent to performing the two-way MSA techniques on a large data matrix formed by unfolding the three-way data array in one of three ways, such as multiway PCA,3 multiway PLS,11 multiway independent component analysis,12 and multiway locality preserving projections.13 The multilinear methods factorize the three-way data array directly using tensor decomposition techniques, such as parallel factor (PARAFAC) decomposition,14 Tucker decomposition,15 generalized Tucker2 decomposition,16 and multilinear PLS17. The common feature of multiway and multilinear methods is that they map the high dimensional batch data into a low dimensional space defined by a few latent variables, such that the important information of the batch data is summarized by some latent variables. Latent variables are a set of combinations of process variables, and they account for correlations among process variables. The latent variables can be used to define monitoring indices for tracking the operation status of new batches, detecting the occurrence of process faults, and diagnosing the causes of faults. Latent variables extracted from the batch data therefore have significant effects on the fault detection and diagnosis performance. Traditional MSA methods, such as PCA, PLS and PARAFAC, can only produce dense latent variables, that is, all process variables are involved in each dense latent variable. The dense latent variables cause two drawbacks for fault detection, analysis and diagnosis. First, the dense latent variables suffer from the strong interference among process variables. As a consequence of the variable interference, the influences of faulty process variables on the latent variables can be partly counteracted by other fault-free process variables. Thus the dense latent variables may be insensitive 3

ACS Paragon Plus Environment

Page 4 of 31

Page 5 of 31

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

to the occurrence of process faults. Second, the dense latent variables usually cannot reveal significant variable correlations corresponding to control loops, unit operations or biochemical reactions in the batch process. This complicates and confounds the physical interpretation of latent variables, and hinders the application of latent variables to detailed fault analysis and diagnosis. To overcome these drawbacks, sparse MSA methods have recently been proposed to extract sparse latent variables from the historical data of continuous processes, e.g., sparse PCA,18,19 sparse global-local preserving projections (SGLPP),20 and the variable-correlation-based sparse modeling method.21 Unlike the dense latent variable, a sparse latent variable consists of only a few process variables. These sparse MSA methods can be easily extended to the aforementioned multiway-type methods for handling the three-way data arrays of batch processes. However, the multiway-type sparse MSA methods also suffer from the problems caused by unfolding the three-way data array to a large two-dimensional data matrix, such as the loss of three-way data characteristics. The multilinear methods have advantages over the multiway methods in modeling the three-way data array, because they can explicitly describe the three-way interactions in the data array. Therefore, developing sparse multilinear methods for analyzing the three-way data array is necessary to improve the batch process monitoring performance. In this paper, a novel sparse parallel factor (SPARAFAC) decomposition is proposed and used for batch process monitoring. The ordinary PARAFAC decomposition is extended to the SPARAFAC decomposition by adding L1-norm penalties on factor vectors. The L1-norm penalized optimization problem of SPARAFAC is solved by the alternating shrunken least squares (ASLS) algorithm. The SPARAFAC is able to produce sparse factor vectors with only a few nonzero elements. A batch 4

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 31

process monitoring method is proposed on the basis of the SPARAFAC decomposition. A sparse monitoring model is built by performing the SPARAFAC decomposition on the three-way batch data array. The T2 and SPE statistics are defined based on the SPARAFAC monitoring model for tracking the running status of new batches and detecting the occurrence of process faults. The multi-level contribution plots, which consist of group-wise and group-variable-wise contribution plots to the T2 statistic and the variable-wise contribution plot to the SPE statistic, are developed for identifying faulty variable groups and faulty variables responsible for process faults. The effectiveness and advantages of the proposed methods are illustrated by a case study in an industrial-scale fed-batch fermentation process The rest of this paper is organized as follows. The next section introduces the optimization problem and solving algorithm of the SPARAFAC decomposition. The SPARAFAC-based batch process monitoring method, consisting of the sparse monitoring model, fault detection indices and multi-level contribution plots, is presented in Section 3. This is followed in Section 4 with the application of the proposed methods to an industrial-scale fed-batch fermentation process. The conclusions are given in Section 5. 2. Sparse parallel factor decomposition 2.1. Problem and Solution Parallel factor (PARAFAC) decomposition factorizes a tensor into a sum of rank-one components. For a third-order tensor  ∈ ℛ ×× , the PARAFAC decomposition is  ≈ ∑  ∘  ∘ 

where R is the number of components, ∈ ℛ are scalars, ∈ ℛ  , 5

ACS Paragon Plus Environment

(1)  ∈ ℛ  and  ∈ ℛ are

Page 7 of 31

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

factor vectors with the unit norm (i.e. ‖ ‖ = ‖ ‖ = ‖ ‖ = 1), and the symbol “◦” denotes the vector outer product. The PARAFAC model in Eq. (1) can be written in the matricized form as22 () ≈ (⨀) () ≈ (⨀) () ≈ (⨀)

(2)

where  = diag(  , ⋯ ,

) is a diagonal matrix,  = &  , ⋯ , ' ∈ ℛ × is a factor matrix (and likewise for B and C), (() denotes the mode-i matricization (i.e., transforming a tensor into a

matrix along the i-th mode) of the tensor ,22 and the symbol “⊙” represents the Khatri–Rao product. The ordinary PARAFAC decomposition produces dense solutions, that is, factor matrices with all elements being nonzero. To achieve sparse solutions (i.e., most of elements of the factor matrices being exactly zero), a sparse PARAFAC decomposition is developed by adding L1-norm penalties to all factor vectors. For a third-order tensor  ∈ ℛ ×× , the sparse PARAFAC decomposition seeks the solution to the following problem 

arg min-. , . ,. ,. / − ∑  ∘  ∘  /1



subject to ‖ ‖ = ‖ ‖ = ‖ ‖ = 1 and ‖ ‖  ≤ ;. for @ = 1, ⋯ , A

(3)

where ‖∙‖1 represents the squared tensor Frobenius norm (the sum of squared elements of the tensor), ‖∙‖ or ‖∙‖ denotes the L2 or L1 norm of a vector, and ;. are penalization

parameters for penalizing the L1 norm of vectors to produce sparse solutions. 2.2. Algorithm The alternating shrunken least squares (ASLS) algorithm proposed by Rasmussen and Bro23 is used to solve the optimization problem in Eq. (3). Consider optimizing Eq. (3) with respect to 6

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 31

C( = ( ( (i = 1, …, R) for any given  = & , ⋯ ,  ',  = & , ⋯ ,  ' and C = (@ =

1, … , A and @ ≠ F). The optimization problem in Eq. (3) is simplified as



C ∘  ∘  − C( ∘ ( ∘ ( / arg min CG / − ∑ , H( 1 subject to ‖ ( ‖ = 1, ‖ ( ‖  ≤ ; j> ⁄(]` − 1) denotes the covariance

matrix of the matrix j> which is the row centered T. The T2 and SPE statistics monitor data

variations in the model subspace defined by j and in the residual subspace defined by k, respectively. Confidence limits of these two fault detection indices are needed to determine the occurrence of a process fault. The T2 statistic follows a weighted F-distribution with R and IK-R degrees of freedom 3

Mx  x ON

r  ~ (O ) y(A, ]` − A)

(16)

where R is the number of components in the SPARAFAC model, and I and K are numbers of batches 10

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 31

and samples in the training data set. Suppose the residuals are independent and normally distributed, then the SPE statistic follows a weighed χ2-distribution with h degrees of freedom3 }

 Qtu ~ z{(|) , z = ~ , ℎ =

~x }

(17)

where m and v are the mean and variance of the SPE statistics of all training samples. Confidence limits of the T2 and SPE statistics at significance level α (e.g., α = 0.01) can be easily obtained through above F- and χ2-distributions. If either the T2 or SPE statistic of a new sample exceeds the corresponding confidence limits, the occurrence of a process fault is detected. 3.3. Fault diagnosis using multi-level contribution plots Once a process fault was detected, the cause of the fault needs to be diagnosed in order to eliminate the fault in a timely manner. To utilize the sparsity of the SPARAFAC model (more specifically, the significant variable correlations revealed by the sparse factor matrix B) to improve the fault diagnosis capability, the two-level contribution plots to the T2 statistic are developed for fault diagnosis. This two-level contribution plots consist of the first level group-wise contribution plot and the second level group-variable-wise contribution plot. The group-wise contribution plot shows the contribution of each variable group (variables corresponding to nonzero elements in each column of B constitute a variable group) to the T2 statistic. The group-variable-wise contribution plot shows the contribution of each variable in a group to the T2 statistic. As for the SPE statistic, it is difficult to define the group-wise and group-variable-wise contribution plots, and therefore the variable-wise contribution plot is used for fault diagnosis. Two-level contribution plots to the T2 statistic and the variable-wise contribution plot to the SPE statistic constitute multi-level contribution plots. 11

ACS Paragon Plus Environment

Page 13 of 31

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 T2 statistic in Eq. (14) can be rewritten as r  = (p − pi ) sO (p − pi ) )

= (l − li (  

)O

€

M N

‚ƒ  Y „(x)  Y  MN‚ƒ   „(x)

Y „(x)  Y   O   „(x)

= (l − li )  †

O

‡

O

 (l − li )

…

O

( )O  (l − li )

= (l − li ) ∑ (,‰M( Qˆ(‰ ‰ N (l − li ) = ∑ ((l − li ) ( Qˆ(( ( (l − li ) + ∑ (,‰,(H‰(l − li ) ( Qˆ(‰ ‰ (l − li )

(18)

= ∑ ( r(( + ∑ (,‰,(H‰ r(‰

Y >() is the column-centered  Y () ,the matrix sŠ is where li is the mean of all training samples,  Y >()  Y >() NO, Qˆ(‰ denotes the (i,j)-th element of sŠ, and bi (i = 1, …, defined as sŠ = (]` − 1)M 

R) is the i-th column of the sparse factor matrix B. Recall that only a few variables are associated with the nonzero elements in each column of B due to the sparsity of B. These variables can be viewed as a group because of the close correlations between them. Thus, the r(( in Eq. (18)

represents the independent contribution of the i-th variable group to the T2 statistic, while the r(‰

denotes the joint contribution of i-th and j-th groups to the T2 statistic. Note that the joint contributions cause the interference between variable groups. This leads to the smearing effect and may reduce the fault diagnosis accuracy. To eliminate the smearing effect, only the independent contribution r(( is used to define the group-wise contribution

‹(, = (l − li ) ( Qˆ(( ( (l − li ) 

= Qˆ(( M∑‰ S̅‰ (‰ N =

 Qˆ(( Ž∑‰MS̅‰ (‰ N

+

∑‰,,‰HMS̅‰ (‰ ( S̅ N

(19)

where S̅‰ = S‰ − Si‰ is the centered measurement of the j-th variable in x, and (‰ is the j-th

element of the sparse vector ( . A variable group is considered faulty if it makes the largest group-wise contribution to the T2 statistic. This faulty group most likely includes the faulty variables responsible for the fault. Then, the group-variable-wise contribution is used to further determine the 12

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 14 of 31

root faulty variable in the faulty group, which is defined as ‹(‰, = Qˆ(( MS̅‰ (‰ N



(20)

where ‹(‰, denotes the independent contribution of the j-th variable in the i-th group to the T2

statistic. The joint contributions of variables, i.e., S̅‰ (‰ ( S̅ (‘ ≠ ’) in Eq. (19), are also ignored in order to eliminate the smearing effect. The SPE statistic in Eq. (15) can be rewritten as Qtu = (l − p) (l − p) = l (“ − ( )O  )l YO”  l = l l − l s   Y‰ Ns Y YO” M∑  = ∑‰ S‰ − M∑‰ S‰  ( ( S( N Y‰ s Y‰ S‰ − S‰  Y‰ s Y( S( N– YO”  YO” M∑(,(H‰  = ∑‰•S‰ − S‰ 

(21)

= ∑‰•Qtu‰‰ + ∑(,(H‰ Qtu‰( –

Y‰ is the j-th Y =  , S‰ (j = 1, …, J) is the measurement of the j-th variable in x, and  where s

Y‰ s Y‰ S‰ in Eq. (21) represents the independent YO”  row of the matrix B. The Qtu‰‰ = S‰ − S‰ 

Y‰ s Y ( S( (F ≠ ‘) denotes YO”  contribution of the j-th variable to the SPE statistic, while Qtu‰( = −S‰  the joint contribution of j-th and i-th variables to the SPE statistic. The joint contribution Qtu‰( is neglected to eliminate the smearing effect, and thus the variable-wise contribution to the SPE statistic is defined as Y‰ s Y‰ S‰ YO”  ‹‰,L—˜ = S‰ − S‰ 

(22)

The variable that makes the largest variable-wise contribution to the SPE statistic is identified as the faulty variable responsible for the fault. 3.4. Modeling and monitoring procedures Offline modeling and online monitoring procedures of the SPARAFAC-based monitoring method are summarized as follows (see Fig. 1). 13

ACS Paragon Plus Environment

Page 15 of 31

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) Offline modeling procedure Step 1: Normalize the training data set X. Step 2: Perform the SPARAFAC decomposition on the training data set. Step 3: Build the SPARAFAC monitoring model via Eq. (12). Step 4: Compute the T2 and SPE statistics of all training samples via Eq. (14) and Eq. (15). Step 5: Determine confidence limits for the T2 and SPE statistics. (2) Online monitoring procedure for a new batch Step 1: Normalize the new sample xk at time k using the mean and variance of training data. Step 2: Compute score and residual vectors of xk using the SPARAFAC monitoring model. Step 3: Compute the T2 and SPE statistics of xk via Eq. (14) and Eq. (15). Step 4: If either the computed T2 or SPE statistic exceeds the confidence limit, report a fault and go to Step 5; otherwise skip to Step 6. Step 5: Identify faulty variables and/or faulty variable groups using multi-level contribution plots. Step 6: Return to Step 1 and monitor the next sample xk+1.

14

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

Fig. 1. Modeling and monitoring procedures of the SPARAFAC-based monitoring method. 4. Case study: industrial-scale fed-batch fermentation process 4.1. Process description The performance of the proposed methods is illustrated by an industrial fed-batch fermentation process.24 The penicillin fermentation is carried out in a 100,000 L bioreactor. Fig. 2 shows the schematic of the bioreactor, with all the process inputs and outputs. The bioreactor pH is kept at 6.5 by adjusting acid/base flow rates. The cooling/heating water flow rates are adjusted to maintain the bioreactor temperature at 298 K. The aeration rate, sugar feed rate, water for injection/dilution, air head pressure, and soybean oil feed rate are adjusted by the sequential batch control strategy which standardizes the input flow rates based on a predefined optimum profile.24 Goldrick et al.24 have developed a simulator, referred to as IndPenSim, for this industrial-scale fed-batch fermentation 15

ACS Paragon Plus Environment

Page 16 of 31

Page 17 of 31

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

process. This simulator provides a platform for use in the development of monitoring strategies for industrial fermentation processes. In this case study, 40 normal batches were simulated by running the IndPenSim with randomized initial conditions and a fixed batch length of 250 hours. In each batch, 16 process variables were measured with a sample rate of 0.2 hours, as shown in Table 1. These normal batches constitute a training data set, X(40×16×1250), for process modeling. Seven faulty batches were also simulated by IndPenSim, with the fault information shown in Table 2. These faulty batches are used for testing the performance of the proposed monitoring method.

Fig. 2. Schematic of the industrial bioreactor.

16

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.Process variables in the fed-batch fermentation process No. Process variables Unit No. Process variables 1 Aeration rate L/h 9 Air head pressure 2 Agitator RPM RPM 10 Dissolved oxygen concentration 3 Sugar feed rate L/h 11 Penicillin concentration 4 Acid flow rate L/h 12 Vessel Volume 5 Base flow rate L/h 13 Bioreactor pH 6 Cooling water flow rate L/h 14 Bioreactor temperature 7 Heating water flow rate L/h 15 CO2 percent in outgas 8 Water for injection/dilution L/h 16 Dissolved CO2

Page 18 of 31

Unit bar mg/L g/L L pH K % mg/L

Table 2. Faulty batches in the fed-batch fermentation process Occurrence End No. Faulty variable Fault type time (h) time(h) 1 Aeration rate Ramp changes from +0% to +5% 20 100 2 Air head pressure Ramp changes from -0% to -5% 80 180 3 Sugar feed rate Ramp changes from +0% to +3% 90 190 4 Base flow rate Step changes by +2% 140 200 120 250 5 Cooling water flow rate Step changes by +7% 6 Bioreactor temperature Ramp changes from +0 to +0.4 K 60 250 7 Bioreactor pH Ramp changes from +0 to +0.1 40 250 4.2. Sparse monitoring model The SPARAFAC decomposition is implemented on the training data set X(40×16×1250) to build the monitoring model. The number R of components is selected to be 11. Penalization parameters are chosen to be ;. = √1250⁄4 for r = 1, …, R. The extracted sparse factor matrices A, B and C have the degrees of sparsity (i.e. the percentage of zero elements in the matrix) of 81.59%, 83.52% and 84.76%, respectively. To demonstrate the capability of SPARAFAC in revealing significant variable correlations, Table 3 shows process variables corresponding to nonzero elements in columns of the variable-wise sparse factor matrix B. It can be seen that the sparse factor matrix B reveals actual control loops or meaningful physical/chemical correlations between process variables. For example, the 4th column of B represents the control loop consisting of the acid flow rate (v4), base flow rate (v5) and bioreactor pH (v13), which is used for the 17

ACS Paragon Plus Environment

Page 19 of 31

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

control of the bioreactor pH (Fig. 2). Similarly, the 6th and 7th columns of B represent the control loop that consists of the cooling water flow rate (v6), heating water flow rate (v7) and the bioreactor temperature (v14). The first column of B reveals the physical correlation between the aeration rate (v1) and the air head pressure (v9), that is, as the aeration rate (v1) increases, the air head pressure (v9) increases as well. The third column of B indicates the chemical correlation between the sugar feed rate (v3) and the bioreactor temperature (v14), since the heat is generated with the consumption of sugar in biochemical reactions. The 5th and 10-11th columns of B also uncover meaningful physical/chemical correlations between variables in groups: {v5, v11}, {v11, v15, v16} and {v10,

v11}. These significant variable correlations are very helpful for subsequent fault analysis and diagnosis. Table 3. Nonzero variables in columns of the sparse factor matrix B Column Nonzero variables Column Nonzero variables 1 7 {v1, v9} (v7, v14) 2 8 v2, v8 v1, v8 3 9 {v3, v14} v2, v3, v5, v6, v9, v11, v12 4 10 (v4, v5, v13) {v11, v15, v16} 5 11 {v5, v11} {v10, v11} 6 / {v6, v7} / (·) and {·} denote the control loop and physical/chemical correlations between variables. 4.3. Fault detection results The monitoring performance of the SPARAFAC-based method was illustrated by 7 faulty batches in Table 2, with the comparison of the PARAFAC-based and MPCA-based methods. The PARAFAC monitoring model was built in a similar way as the SPARAFAC model, whereas the PARAFAC decomposition extracted dense factor matrices from the training data set. The MPCA monitoring model was built by unfolding the training data set X(40×16×1250) along the variable way to a matrix X(50000×16) firstly, and then modeling the unfolded matrix X(50000×16) via PCA. A total of 11 18

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

components were retained in the PARAFAC and MPCA models. Three monitoring models were used with the T2 and SPE statistics for fault detection. Confidence limits for T2 and SPE statistics were set to be 99%. The fault detection performance of T2 and SPE statistics were quantified by two metrics: fault detection rate (FDR), defined as the ratio between the number of detected faulty samples and the total number of samples collected during the faulty operating period; fault detection delay (FDD), defined as the time difference between the actual fault occurrence time and the detected fault occurrence time. Table 4 shows fault detection results of three monitoring methods for 7 faulty batches. Both the SPARAFAC-based and PARAFAC-based methods performed better than the MPCA-based method by virtue of higher FDRs and shorter FDDs, especially for Faults 4-6. This is because the PARAFAC decomposition has advantages over MPCA in modeling the three-way data array. For example, PARAFAC can explicitly describe the three-way interactions in the batch data set, while MPCA destroys the three-way data structure due to data unfolding. For all the 7 faults, the SPARAFAC-based method has comparable performance to the PARAFAC-based method, since the largest FDR of the SPARAFAC-based method for each fault was close to that of the PARAFAC-based method. The SPARAFAC T2 statistic, however, has higher FDRs than the PARAFAC T2 statistic for most faults, with the exception of Fault 6. Unlike PARAFAC, the SPARAFAC extracts sparse factor matrices from the batch data. These sparse factor matrices enable SPARAFAC not only to reveal significant variable correlations but also to reduce the interference among process variables. This improves the sensitivity of the SPARAFAC T2 statistic to the faults caused by variables corresponding to larger nonzero elements in the sparse factor matrix B, and thus 19

ACS Paragon Plus Environment

Page 20 of 31

Page 21 of 31

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

increases the FDRs for Faults 1-5. However, for the faults caused by variables corresponding to small nonzero elements in the sparse factor matrix B, most of the fault information is captured by the residual subspace (defined by k) rather than the model subspace (defined by j); consequently, the SPARAFAC T2 statistic may have lower sensitivity than the SPE statistic to these faults. For example, the SPARAFAC SPE statistic outperforms the T2 statistic in detecting Fault 6, because the faulty variable, bioreactor temperature (v14), corresponds to two very small nonzero elements (≈ 0.0556) in the sparse matrix B. It is easy to infer that increasing the sparsity of the factor matrix B may result in some inactive variables corresponding to no nonzero elements in the matrix B. Faults caused by inactive variables cannot be detected by the SPARAFAC T2 statistic, because no information on the inactive variables is passed to the T2 statistic. Such faults, however, may be well detected by the SPARAFAC SPE statistic. Table 4. Fault detection results of three methods for 7 fault batches SPARAFAC PARAFAC MPCA 2 2 2 Fault T SPE T SPE T SPE No. FDR FDD FDR FDD FDR FDD FDR FDD FDR FDD FDR FDD (%) (h) (%) (h) (%) (h) (%) (h) (%) (h) (%) (h) 1 0.75 73.8 91.52 5.0 77.31 17.4 91.02 5.0 90.77 5.0 90.02 6.4 2 91.62 6.2 0.0 / 87.43 8.4 92.42 6.0 89.22 8.4 0.0 / 3 0.4 80.6 81.84 12.6 87.62 11.2 83.43 12.2 22.95 42.8 87.62 8.8 4 1.0 10.6 80.73 0.2 2.66 2.0 75.75 0.4 18.60 1.0 82.39 0.2 5 0.2 19.35 0.2 33.64 0.2 16.74 1.0 28.57 0.2 18.89 0.2 34.87 6 40.06 4.8 96.11 1.0 99.16 1.0 33.86 6.2 74.76 1.6 84.96 1.6 7 99.05 1.0 83.25 21.6 99.05 1.0 69.36 7.6 98.48 2.2 94.29 6.4 The highest FDR for each fault is highlighted in bold. Fig. 3 and Fig. 4 illustrate the differences between monitoring charts of three methods for Faults 3 and 6. Fault 3 was caused by ramp changes in sugar feed rate (v3) in the time period between 90 h and 190 h, and therefore it leads to ramp increases in T2 or SPE statistic, as shown in all monitoring charts in Fig. 3. The SPARAFAC T2 statistic detected the occurrence of Fault 3 at 98.8 h (with a FDD 20

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

of 8.8 h), which was 2.4 hours and 3.4 hours earlier than the fault occurrence time (at 101.2 h and 102.2 h) detected by the PARAFAC SPE statistic and the MPCA T2 statistic. Since a shorter FDD implies the better fault detection performance, the SPARAFAC method is better than the PARAFAC and MPCA methods in detecting Fault 3. In Fig. 4a and Fig. 4b, both the SPARAFAC SPE statistic and the PARAFAC T2 statistic detected Fault 6 at 61.0 h (with a FDD of 1.0 hour), and they clearly indicated that Fault 6 lasted until the end of the batch. However, the MPCA T2 and SPE statistics detected the Fault 6 after a longer FDD of 1.6 hours (Fig. 4c). Besides, many samples in the period between 162 h and 250 h were not detected by either the MPCA T2 or SPE statistic, resulting in a lower FDD for Fault 6 (Table 4). The above two examples illustrate how the SPARAFAC method outperforms the PARAFAC and MPCA methods in fault detection.

21

ACS Paragon Plus Environment

Page 22 of 31

Page 23 of 31

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

Fig. 3. Monitoring charts of (a) SPARAFAC, (b) PARAFAC and (c) MPCA for Fault 3.

22

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

23

ACS Paragon Plus Environment

Page 24 of 31

Page 25 of 31

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

Fig. 4. Monitoring charts of (a) SPARAFAC, (b) PARAFAC and (c) MPCA for Fault 6. 4.4. Fault diagnosis results The faulty variables responsible for faults are identified by the multi-level contribution plots. Fault 3 and Fault 6 are used as two examples to illustrate how the multi-level contribution plots produce reliable fault diagnosis results. Fault 3 was first detected by the SPARAFAC T2 statistic at 98.8 h (Fig. 3a). Therefore, the two-level contribution plots to the T2 statistic at 98.8 h are used for fault diagnosis, as shown in Fig. 5. The group-wise contribution plot in Fig. 5a indicates that the group 3 makes a larger contribution than other groups. The group 3 is thus identified as the faulty group. The group-variable-wise contribution plot in the group 3 is shown in Fig. 5b. Clearly, the contribution of sugar feed rate (v3) is much larger than that of bioreactor temperature (v14). Therefore, the sugar feed rate is identified as the faulty variable responsible for Fault 3. This diagnosis result aligns well with the actual cause of Fault 3, listed in Table 2 as the ramp changes in sugar feed rate. Next, consider Fault 6, which was first detected by the SPARAFAC SPE statistic at 61.0 h (Fig. 4a). Fig. 6 shows the variable-wise contribution plot to the SPE statistic at 61.0 h. The bioreactor temperature (v14) 24

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

makes the largest contribution, and hence it is identified as the faulty variable responsible for Fault 6. Again, this diagnosis result is consistent with the actual cause of Fault 6, listed in Table 2 as the ramp changes in bioreactor temperature (v14).

Fig. 5. Two-level contribution plots to the T2 statistic at 98.8 h for Fault 3. (a) Group-wise contribution, (b) group-variable-wise contribution.

Fig. 6. Variable-wise contribution plot to the SPE statistic at 61.0 h for Fault 6. Table 5 summarizes the diagnosis results of multi-level contribution plots for all the 7 faults. By comparing the results in Table 5 with the fault information in Table 2, it can be seen that actual fault variables responsible for all faults are correctly identified by multi-level contribution plots. Moreover, 25

ACS Paragon Plus Environment

Page 26 of 31

Page 27 of 31

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 group-wise contribution plot to T2 statistic reveals the faulty variable groups containing the faulty variables. The significant variable correlations (Table 3) in variable groups enable us to associate faults with actual control loops or physical/chemical links (each link consists of the variables with close physical/chemical correlations) in the process, which provides more useful information for fault diagnosis, analysis and elimination, such as inferring the root causes of faults, analyzing potential influences of faults to the process, and taking special fault-removing measures. For example, the group 4 was identified as the faulty group for Faults 4 and 7 (Table 5). Because the group 4 represents the control loop of pH (Table 3), we can infer that the disturbances to the control loop of pH may be the root causes of Faults 4 and 7. The bioreactor pH may change due to the occurrence of Faults 4 and 7, which slows down the growth rate of biomass and the production rate of penicillin. A possible measure to eliminate Faults 4 and 7 is adjusting the acid/base flow rates. The results in Table 5 indicate that the multi-level contribution plots are able to produce reliable fault diagnosis results by identifying faulty variable groups and faulty variables simultaneously. Table 5. Diagnosis results of multi-level contribution plots for 7 faults Fault No. DFOT (h) Statistic Faulty groups Faulty variable 1 25.0 T2 g1 v1 (Aeration rate) 2 2 86.2 T g9 v9 (Air head pressure) 3 98.8 T2 g3 v3 (Sugar feed rate) 2 g4 v13 (Bioreactor pH),v5 (Base flow rate) 4 140.2 T 2 5 120.2 T g6 v6 (Cooling water flow rate) 6 61.0 SPE / v14 (Bioreactor temperature) 2 g4 v13 (Bioreactor pH), v5 (Base flow rate) 7 41.0 T DFOT: detected fault occurrence time. 5. Conclusions In this paper, a novel sparse parallel factor (SPARAFAC) decomposition is developed and used for batch process monitoring. The SPARAFAC adds additional L1-norm constraints on the optimization 26

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

problem of the ordinary PARAFAC to produce sparse factor vectors. Only a few variables are active in each component of the SPARAFAC model due to the sparsity of factor vectors. This enables SPARAFAC to eliminate redundant data information and to reveal significant variable correlations. Therefore, SPARAFAC has the potential advantage of being easier to interpret as compared to the ordinary PARAFAC, where all variables appear in every component. A batch process monitoring method is proposed on the basis of the SPARAFAC decomposition. The SPARAFAC model is well suitable for fault analysis and diagnosis due to the model sparsity and the enhanced model interpretability. Based on the SPARAFAC model, the T2 and SPE statistics are defined for fault detection. Moreover, to use the significant variable correlations revealed by SPARAFAC to improve the fault diagnosis capability, multi-level contribution plots, which consist of the group-wise and group-variable-wise contribution plots to the T2 statistic and the variable-wise contribution plot to the SPE statistic, are developed for fault diagnosis. The effectiveness and advantages of the proposed methods are illustrated by a case study in an industrial-scale fed-batch fermentation process. The results indicate that the variable-wise sparse factor matrix in the SPARAFAC model reveals actual control loops or meaningful physical/chemical correlations among process variables. The SPARAFAC-based monitoring method outperforms the PARAFAC-based and MPCA-based methods in terms of better model interpretability, higher fault detection rate and shorter fault detection delay. Reliable fault diagnosis results are obtained by multi-level contribution plots. In particular, the group-variable-wise and group-wise contribution plots to the T2 statistic not only identify the faulty variables responsible for faults but also the faulty variable groups associated with the control loops or physical/chemical links (each link consists of the 27

ACS Paragon Plus Environment

Page 28 of 31

Page 29 of 31

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

variables with close physical/chemical correlations) in the process. Such fault diagnosis results are help to locate the fault accurately. Appendix A. Sparse parallel factor decomposition algorithm Algorithm 1: Sparse parallel factor (SPARAFAC) decomposition 1. Input the tensor  ∈ ℛ ×× and penalization parameters ;. for r = 1, …, R 2. Initialize  ∈ ℛ × ,  ∈ ℛ × ,  ∈ ℛ × and  ∈ ℛ× 3. Repeat until converge: I. Iteratively update ( and ( (i = 1, …, R) until convergence: O

  C(,KL ← M() − ( ( I,( (a) NJ,( MJ,( J,( N

C(,KL , R