Multiple Hypotheses Testing-Based Operating Optimality Assessment

Article pubs.acs.org/IECR

Multiple Hypotheses Testing-Based Operating Optimality Assessment and Nonoptimal Cause Identification for Multiphase Uneven-Length Batch Processes Yan Liu,† Fuli Wang,*,†,‡ Yuqing Chang,†,‡ Ruicheng Ma,§ and Shumei Zhang† †

College of Information Science & Engineering, Northeastern University, 3 Lane 11, Wenhua Road, Heping District, Shenyang, Liaoning 110819, China ‡ Stat Key Laboratory of Synthetical Automation for Process Industries, Northeastern University, Shenyang, Liaoning 110819, China § School of Mathematics, Liaoning University, Shenyang 110036, China ABSTRACT: As a novel research issue, the operating performance assessment for industrial processes has received growing interest in recent years. In this study, a multiple hypotheses testing-based operating optimality assessment and nonoptimal cause identification method is proposed to handle the online assessment of the batch processes with multiphase and uneven-length characteristics. A new Gaussian mixture model (GMM)-based offline phase division procedure is proposed. Compared to the existing methods, it can not only ensure that each phase contains the consecutive sampling instants but can also retain the characteristics of phase-uneven-length for retaining more-accurate phase information. The extra postprocessing is not needed in phase division, and the assessment models of each phase are set up along with it. In online application, both the local and global online assessments are performed for a new batch based on the multiple hypotheses testing technology by controlling the false discovery rate. For the nonoptimal batch, the cause variables can be determined based on the variable contributions to the assessment index. Finally, the effectiveness of the proposed method is verified through a fed-batch penicillin fermentation process. batch processes,6−9 the products can be divided into qualified and unqualified ones, according to the monitoring results. In contrast, the purpose of operating optimality assessment for batch processes is to further divide the qualified products into several grades, such as optimal and nonoptimal, depending on the quality, raw material costs, and energy consumption, determine the causes for the nonoptimal goods, and provide a reasonable reference for further production adjustment. As a novel research issue, the operating performance assessment for industrial processes has received growing interest in recent years.10−14 Liu et al.10−12 proposed the operating optimality assessment methods based on feature extraction and matching for continuous processes. Also, the evaluation strategies for the continuous multimode processes are developed under probabilistic frameworks.13,14 However, because of the inherent natures, such as multiphase, finite operation durations, time-varying, and uneven-length duration, etc., the existing operating optimality assessment methods cannot be directly applied to batch processes, which makes the

1. INTRODUCTION In modern industries, batch and semi-batch processes play an important role in the production of low-volume and high-valueadded products, including the specialty chemical, semiconductor, food, and biology industries, to meet the rapidly changing market.1−5 In order to maximize economic benefits, the batch process product is usually considered synthetically from the aspects of the quality, raw material costs, energy consumption, etc. Then, a desired quality standard can be determined based on the properties and application of the product, which can not only ensure the product quality but also improve the comprehensive economic benefits. Moreover, in industrial processes, the product quality is not necessarily viewed as being the higher the better, but it is hoped that the product quality can achieve or be as close as possible to the specified quality standard. Hence, for the products that are all consistent with the specified quality standard, the corresponding comprehensive economic benefits may be still unequal for the differences in raw material costs and energy consumption. That is to say, the blind pursuit of high product quality without considering the impact factors cannot guarantee the comprehensive economic benefits, and it is needed to produce low-cost and high-quality products, from the perspective of long-term development of the enterprises. In the quality monitoring of © 2016 American Chemical Society

Received: Revised: Accepted: Published: 6133

December 15, 2015 March 23, 2016 April 13, 2016 May 17, 2016 DOI: 10.1021/acs.iecr.5b04775 Ind. Eng. Chem. Res. 2016, 55, 6133−6144

Article

Industrial & Engineering Chemistry Research

uneven-length in achieving more accurate phase information. In addition, the proposed method does not need post-processing and additional process monitoring, which makes it simpler. In online application, phase identification for the new batch is the basis of online operating optimality assessment, since the result of phase identification can indicate the right model to be used to evaluate the current operating performance. Based on the durations of each phase in each training batch obtained from the offline phase division procedure, the maximum operating ranges of each phase, i.e., the earliest start points and the latest ending points, can be determined and used for achieving the preliminary phase identification. According to the maximum operating ranges of each phase, the fuzzy region between two adjacent phases is often inevitable, where the phase ownership cannot be identified for certain. However, this can be further determined according to the online assessment process. In online operating optimality assessment, both the local and global assessments are performed on the new batch. The local assessment works to master the operating performance of a period of time within the batch, and the results are used to guide and adjust the subsequent production of the current batch. In order to avoid external interference and obtain reliable assessment result, a data window including a series of samples is used to represent a period of time and serves as the basic analysis unit in local assessment. All the samples within the data window should be seen as a whole and each sample is evaluated simultaneously to reduce the false positive rate to the full extent, which is quite consistent with the idea of multiple hypotheses testing. Multiple hypotheses testing technology involves tackling the multiplicity problems caused by a large number of hypotheses under the assumption of independent distribution by controlling the error rate at a low level, and it has been widely applied in process monitoring and fault isolation.20−23 In this study, the multiple hypotheses testing method is used in process operating optimality assessment for the first time. Each sample of the data window can be considered as a hypothesis, and the false discovery rate (FDR) is controlled by using the Benjamini−Hochberg (B−H)24 procedure, to ensure the reliability of the assessment results. Furthermore, the global assessment is carried out based on the local assessment results. It is committed to making a comprehensive evaluation of the entire batch at the end, and the results are helpful in improving and adjusting the productions of future batches. When the nonoptimal operating performance rises in batch processes, the cause variables should be further identified. Similar to the contribution plots-based fault diagnosis,5,25,26 the variable contributions are calculated and used to determine the responsible cause variables for the nonoptimal operating performance. The rest of this paper is organized as follows. Section 2 gives some preliminaries of GMM and multiple hypotheses testing technology. Then, in section 3, the detailed description of the proposed approach is introduced, which includes the offline phase division and modeling, online phase identification and assessment, and nonoptimal cause identification. In section 4, the utility of the proposed batch assessment approach is demonstrated with a simulated fed-batch penicillin fermentation example. Finally, section 5 concludes this paper.

operating optimality assessment for batch processes quite challenging. Considering the fact that different phases contain different underlying process characteristics and exhibit significantly different behaviors, the assessment models should be established for each phase respectively, and then the offline phase division of training batches becomes a key issue before the establishment of assessment models. In recent years, some phase division methods for multiphase batch processes have been developed. For simplicity, the uneven length for the total batch durations is referred to as batch-uneven-length and that occurs in phases is called as phase-uneven-length in the following content. Lu et al.15 proposed a stage-based sub-PCA method to address the multiphase batch processes, where a kmeans clustering technique, as well as some post-processing, is adopted to partition the PCA loading matrices into several clusters for offline phase division. As an improvement, Zhao et al.16 developed the stage-based soft-transition multiple PCA offline phase division method. The entire process can be properly divided into different substages and the transition regions, which overcomes the disadvantage of Lu’s hardpartition strategy. However, both methods were developed based on the assumption of equal duration and ill-suited for the uneven-length batch processes. To this end, Yu et al.17 developed the Gaussian mixture model (GMM) to batch processes and proposed the multiway Gaussian mixture modelbased phase division method, where the hybrid unfolded multiway data matrix was partitioned into different clusters by GMM, and then the post-processing procedures were implemented to adjust clusters to satisfy the following two conditions: (i) each of them only contains consecutive sampling instants (ii) all the training batches at the same sampling time belong to the same cluster. Nevertheless, condition (ii) is ill-suited for the process with phase-uneven-length, and the phase division results are inaccurate if the above approach is directly applied to this type of batch process. Lu et al.4 extended their original stagebased sub-PCA method15 to the batch processes with phaseuneven-length and batch-uneven-length. The phases were derived in series rather than in parallel, with respect to the process sequence based on the clustering algorithm, and then the durations of each phase can be determined precisely by checking whether the SPE statistic is beyond a predetermined threshold. Also, many other works have contributed their efforts in solving the phase division of batch processes.18,19 These pioneer works have provided abundant theoretical bases for our work. Under the framework of operating optimality assessment, the offline phase division issue for multiphase batch processes with both phase-uneven-length and batch-uneven-length is solved through a new thought in this study. The GMM is first estimated based on the variable-wise unfolded data matrix. The phases then are examined one by one, by checking the posterior probabilities of the training samples with the target Gaussian components. After all phases have been determined, the Gaussian component parameters, i.e., mean vectors and covariance matrices, are updated by the data of each phase and used as the parameters of the assessment models of the corresponding phase. The proposed phase division method can not only ensure each phase containing a series of consecutive sampling instants but also retain the characteristics of phase6134

DOI: 10.1021/acs.iecr.5b04775 Ind. Eng. Chem. Res. 2016, 55, 6133−6144

Article


referenced as “family-wise testing”, and testing each hypothesis of the family-wise testing simultaneously. It provides a framework for defining and controlling appropriate error rates in order to protect against wrong conclusions, from which one can conclude which tests are significant and which are not. Assume that there are W hypotheses that need to be tested at the same time, and then they are recorded as Hr (r = 1, 2, ..., W). In addition, let θr (r = 1, 2, ..., W) and Θ represent a series of parameters and their space, respectively. If Hr is true, Hr should be denoted as Hr,0: θr ∈ Θ, which is also called a null hypothesis; otherwise, Hr should be recorded as Hr,1: θr ∉ Θ with the means of alternative hypothesis. Then, all the W hypotheses constitute a family-wise testing {Hr, r = 1, 2, ..., W} with Hr = (Hr,0, Hr,1). Until now, many error rates have been proposed,32−34 and the most popular FDR-controlling multiple hypotheses testing is considered in the pioneering article by Benjamini and Hochberg.24 Sometimes, it is even referenced as the B−H procedure. The FDR is the expected proportion of falsely rejected hypotheses among all of the rejected hypotheses. It has been theoretically proven that, when the number of hypotheses is large, controlling FDR is much more powerful than controlling the traditional type I error rate. With the predefined span size W, the control level q and a suitable test statistic, the B−H procedure is implemented step by step as follows: Step 1: Sort the W p-values of the test statistics in increasing order, such that p(r) is the rth smallest p-value. Let H(r) be the null hypothesis corresponding to p(r). Step 2: Initially set r̂ = 0. Then, determine r̂, the maximum value of r (for r = 1, 2, ..., W), such that

2. PRELIMINARIES 2.1. Gaussian Mixture Model (GMM). GMM offers a natural representation in characterizing the inherent heterogeneity of a random variable that possibly comes from a finite number of latent classes.27,28 For the training data matrix X ∈ RN×J with N samples and J variables, each sample x ∈ RJ×1 can be assumed to follow a GMM with M components as follows: M

G(x|Φ) =

∑ ωmg(x|μm , Σm)

(1)

m=1

with the multivariate Gaussian density function for the mth Gaussian component (Cm) being defined as g (x|θm) =

⎤ ⎡ 1 1 exp⎢ − (x − μm )T Σ−m1(x − μm )⎥ 1/2 ⎦ ⎣ 2 (2π ) |Σm| J /2

(2)

where Φ = {ω1, ω2, ..., ωM, θ1, θ2, ..., θM} is the set of all the estimated parameters of the GMM; ω1, ω2, ..., ωM are the prior probabilities of M Gaussian components and satisfy the condition of ∑M m=1 ωm = 1; θm = (μm, Σm) is the model parameter of the mth Gaussian component Cm, μm, and Σm are the corresponding mean vector and covariance matrix, respectively. Here, the above GMM are estimated with Figueiredo−Jain (F-J) algorithm29 for its ability in automatically optimizing the number of Gaussian components and estimating their statistical distribution parameters, which involves the following two step iterations: E-step and M-step. E-step: P(s)(Cm|x n) =

ωm(s)g (x n|μm(s) , Σ(ms)) M

∑m = 1 ωm(s)g (x n|μm(s) , Σ(ms))

⎧ ⎛r ⎞ ⎫ r ̂ = max⎨r |p(r) ≤ ⎜ ⎟q⎬ ⎝W ⎠ ⎭ ⎩

(3)

M-step: μm(s + 1) =

N ∑n = 1 P(s)(Cm|x n)x n N ∑n = 1 P(s)(Cm|x n)

Step 3: If r̂ is nonzero, the hypotheses H(1), H(2), ..., H(r̂) are rejected; otherwise, all of the null hypotheses should be accepted. In the above procedure, the p-value corresponding to (marginal) test statistic is supported on the unit interval [0, 1] and represents the probability of the observation or a more extreme case under the situation that the original hypothesis is correct. In probability theory, p-value is reflected by the value of the cumulative distribution function, rather than that of the probability density function. The precise definition of the pvalue is described in ref 35. However, it is worth noting that the methods used for computing the p-value are different, depending on the problem, and should be designed by the users.

(4)

N

Σ(ms+ 1) =

∑n = 1 P(s)(Cm|x n)(x n − μm(s + 1) )(x n − μm(s + 1) )T N

∑n = 1 P(s)(Cm|x n) (5)

ωm(s + 1)

=

N

{

max 0, ∑n = 1 P(s)(Cm|x n) − M ∑m = 1 max

{0,

N ∑n = 1 P(s)(Cm|x n)

V 2

}

−

V 2

}

(7)

(6)

where s is the serial number of F-J iterations, V represents the total number of scalar parameters specifying each Gaussian component (V = J2/2 + 3J/2), and xn comes from the training data X = [x1, x2, ..., xN]T (where n = 1, 2, ..., N). The parameters (s+1) (s+1) μ(s+1) are the mean vector, covariance matrix, m , Σm , and ωm and prior probability of the mth Gaussian component Cm at the (s+1)th iteration, respectively. The F-J algorithm is implemented in an iterative way until the parameter values converge to the optimal solution.30 2.2. Multiple Hypotheses Testing by Controlling the False Discovery Rate. The more questions you ask, the more wrong answers you are expected to receive, even if every single source of your information is quite trustworthy. From the mathematical point of view, this leads to multiple hypotheses testing problems.31 Multiple hypotheses testing is committed to seeing multiple single hypothesis testing as a whole, which is

3. METHOD DEVELOPMENT Generally, the batch processes consist of a sequence of discrete phases or stages. In our opinion, the phases normally are used to describe the parts of a batch process with different process characteristics, variable relationships, or data distributions, while the stages are more often related to the physical operation units. Therefore, the phases and the stages corresponding to the same batch process may be inconsistent. For example, there are two physical operation stages in the fedbatch penicillin fermentation process, i.e., preculture stage and fed-batch stage, but there can be more than two partitioned phases, based on the variable relationships.1,2 Another example is the injection molding process, which is usually composed four stages: injection, packing−holding, plastication, and 6135


Article


Figure 1. Illustration of the variable-wise unfolding for uneven-length training batch.

cooling. However, there are up to 10 phases obtained by the phase division method proposed in ref 1. In this study, whether the focus of discussion is the offline phase division or the online phase identification, we are committed to identify the phases of the batch process according to their inherent data distribution characteristics, which may cause a different division result with the real physical operation stages but can help to learn more about the intrinsic process characteristics and enhance process understanding. 3.1. GMM-Based Phase Division and Assessment Modeling for Uneven-Length Batch Processes. In offline modeling, only the batch data under optimal operating performance are collected and used as the training batch, for which the online new batch can be compared with the optimal assessment model to clearly identify the online operating performance. If the online new batch has a good match with the assessment model, the online operating performance is likely to be optimal; otherwise, nonoptimal operating performance may occur. Assume that there are I optimal and uneven-length batches collected from the historical data set, and the ith batch is denoted as Xi = [xi,1,xi,2, ...,xi,Ki]T ∈ RKi × J (where i = 1, 2, ..., I, and Ki is the number of samples in Xi and also represents the duration of the i th batch, and xi,k indicates the kth sample in the ith batch). In order to establish the GMM, the I training batches are unfolded into a two-dimensional matrix X = [XT(k=1), XT(k=2), ..., XT(k=kmin), ..., XT(k=kmax)]T ∈ RK×J via variable-wise unfolding, as illustrated in Figure 1, where K = ∑Ii=1 Ki, Kmin is the minimum value of the durations (Kmin = min1≤i≤I {Ki}), and Kmax is the maximum value of the durations (Kmax = max1≤i≤I {Ki}); X(k=K̃ ) is the matrix constructed with the K̃ th samples of each batch and satisfies the relation X(k=K̃ ) ∈ RI×J when K̃ = 1, 2, ..., Kmin or X(k=K̃ ) ∈ RI×̃ J, in the case of K̃ = Kmin + 1, ..., Kmax and I ̃ < I. The normalization method proposed in refs 3 and 4 then is implemented on X for phase division, i.e.,

x k̃ , i , j =

x k , i , j − xj̅ sj

(8)

where xk,i,j is the jth variable of the sample at moment k of the ith batch, and xj̅ and sj are, respectively, the mean and standard deviation of the jth variable in X. For simplicity, the normalized form is still denoted as X. After the GMM is estimated by the Bayesian inference based F-J algorithm, M Gaussian components are obtained and denoted as C1, C2, ..., CM. According to the difference in data distribution characteristic, the batch processes can be divided into several phases. In the same phase, the data distribution characteristics are approximately consistent, while the data distribution characteristics are significantly different across different phases. In addition, the uneven-length feature is reflected not only in the batch but also in the phase, which is illustrated by Figure 2a. In view of this, the GMM-based phase division and assessment modeling method is proposed for multiphase uneven-length batch processes in this study.

Figure 2. Schematic of phase division for multiphase uneven-length batch processes. 6136


Article

Industrial & Engineering Chemistry Research Considering the fact that the obvious changes of the data distribution only occur after the process runs into the next phase and has nothing to do with the running time, the posterior probabilities, with respect to each Gaussian component, can be used to determine the phase ownership of the training data, to realize the offline phase division of uneven-length batch processes. Based on Bayesian inference, the posterior probabilities of the samples from phase 1, with respect to the target Gaussian component, should be much larger than those of other Gaussian components. Hence, the target Gaussian component of phase 1 can be determined according to the posterior probabilities of the first h samples from each batch, with respect to all of the Gaussian components. Using h samples rather than one can not only avoid the interference of process noise but also make the results more reliable and more consistent with the statistical law. The means of posterior probabilities with each Gaussian component then are calculated, and the Gaussian component with the maximum mean is considered as the target Gaussian component for phase 1. Thereafter, only the posterior probabilities, with respect to the target Gaussian component, are calculated from the (h + 1)th sample of each batch and are used to judge the ending point based on the predefined phase division rule, until phase 1 of all batches are divided. Through the above procedure, one can find the duration of each phase by checking whether the data can be described by the target Gaussian component with the corresponding posterior probabilities. The detailed procedure of the offline phase division and assessment modeling is described below. First of all, calculate the posterior probabilities of the first h samples of each batch with all of the respective Gaussian components as follows:

P(1), ̅ h′ =

M

(9)

where the value of h should be smaller than the shortest duration of the phases and can be determined according to the process knowledge and expert experience. The means of posterior probabilities with each Gaussian component then are given by Pm̅ ,0

h

∑ ∑ P(Cm|x i ,h′)

for m = 1, 2 , ..., M

i = 1 h ′= 1

(10)

Find the Gaussian component that satisfies the following condition: m* = arg

max

m = 1,2 , ..., M

{Pm̅ ,0}

(11)

where m* represents the number of the Gaussian component that corresponds to the maximum value of P̅m,0, so Cm* is the target Gaussian component for phase 1. Denote Cm* and P̅m,0 as C(1) and P̅(1),0, respectively. Furthermore, from the (h + 1)th sample, the posterior probabilities of the samples of each batch corresponding to C(1) and their means are respectively calculated as P(C(1)|x i , h ′) =

ω(1)g (x i , h ′|μ(1) , Σ(1)) M

∑m = 1 ωmg (x i , h |μm , Σm) ′ for i = 1, 2, ..., I ; h′ = h + 1, h + 2 , ...

(14)

where i* is the number of the batch where phase 1 is over. Note that, when P̅(1),h′ < δ, the posterior probabilities of some samples in batch i* have been less than δ. In order to determine the ending point of the phase accurately, trace back the latest sample with the posterior probability larger than δ for the batches satisfying eq 14 and denote the corresponding moment as the accurate ending point of phase 1 of batch i*. Then exclude the batches where phase 1 has been over and recalculate the means of posterior probabilities with Gaussian component C(1) of the rest of the samples until the data belonging to phase 1 of all batches have been identified. The mean vector μ(1) and covariance matrix Σ(1) of Gaussian component C(1) are updated by the data of phase 1 and used as the assessment model parameters for operating optimality assessment. Thereafter, with the information on phase 1, the corresponding data in each training batch can be removed and the remaining data are aligned in the head, as shown in Figure 2b. The same procedure can readily be applied to determine the duration of phase 2, and so on. Hence, the phase division is conducted iteratively until all phases are separated. The proposed phase division method can not only ensure that each phase contains a series of consecutive samplings but also retains the characteristics of phase-uneven-length in achieving more-accurate phase information. The step-by-step procedure of the offline phase division and assessment modeling is given as follows: (1) Collect batch data under optimal operating performance for phase division and assessment modeling. (2) Unfold the three-dimensional batch data into a twodimensional matrix X via variable-wise unfolding and make the normalization on X. (3) Use the F-J algorithm to estimate the GMM and obtain several Gaussian components. (4) Align all batches in the head and then use the first h samples of each batch to determine the target Gaussian component of the current phase. (5) At the h′th sample of the current phase (h′ = h + 1, h + 2, ...), calculate the mean of the posterior probabilities corresponding to the target Gaussian component at the same moment and compare it with δ. If it is larger than δ, repeat this step for the next sample; otherwise, go to step (6). (6) Compare the posterior probabilities of each batch, with respect to the target Gaussian component with δ, in turn, and exclude the batches that are over for the current

∑m = 1 ωmg (x i , h ′|μm , Σm)

I

(13)

i=1

i* = arg{i|P(C(1)|x i , h + 1) < δ , i = 1, 2 , ..., I }

where i = 1, 2 , ..., I ; m = 1, 2 , ..., M ; and h′ = 1, 2 , ..., h

1 = Ih

I

∑ P(C(1)|xi ,h ′)

where μ(1), ∑(1), and ω(1) are the mean vector, covariance matrix, and prior probability of the Gaussian component C(1), respectively. In order to determine the duration of phase 1 in each batch, the posterior probability threshold δ(0 < δ < 1) is introduced and can be determined by process knowledge and experience. If P̅(1),h′ > δ, it can be temporarily considered that the h′th sample of each batch still belongs to phase 1, and continues to calculate the posterior probabilities and their mean for the next sample; otherwise, it can be inferred that phase 1 of some batches is over, which leads to a decrease in P̅(1),h′. These batches can be determined as follows:

ωmg (x i , h ′|μm , Σm)

P(Cm|x i , h ′) =

1 I

(12) 6137


Article


Figure 3. Schematic diagram of online phase identification.

3.3. Online Operating Optimality Assessment Based on Multiple Hypotheses Testing. In this section, the online assessment strategy for process operating performance will be introduced from two aspects: local assessment and global assessment. In local assessment, the online data window with local samples is used as the basic analysis object, which is focused on mastering the current operating performance covered by the online data window, and the local assessment result can be used to guide and adjust the subsequent production of the current batch. By contrast, the global assessment is committed to making a comprehensive operating performance evaluation for the entire batch based on the local assessment results, and the results are helpful in improving and adjusting the productions of the future batches. 3.3.1. Local Assessment. Considering that a single sample is not enough to characterize the process operating performance and is susceptible to outliers, an online data window with width W, denoted as Xnew, is constructed as the basic analysis unit for local assessment and is shown below:

phase; then go back to step (5) until the data belonging to this phase of all batches have been identified. (7) Remove the data of the previous phase from each batch, and go back to step (4) to divide the next phase until all phases are separated. (8) Update the mean vectors and covariance matrices of the Gaussian components corresponding to each phase by the data after phase division, and use them as the assessment model parameters for online operating optimal assessment. 3.2. Online Phase Identification. In the procedure of offline phase division, the durations of each phase of each batch are obtained, and then the earliest start points and the latest ending points of each phase can be determined, which are useful for the online phase identification of the new batch Xnew. Take a batch process with four phases as an example, and the detailed phase identification process is described below. First, denote the earliest start point and the latest ending point of phase m as tm,in and tm,out (where m = 1, 2, 3, 4), respectively. A time axis then can be constructed with those points, as shown in Figure 3, where the blank and shaded parts indicate the certain intervals of phases and the fuzzy regions between two adjacent phases, respectively. Here, k represents the current moment of Xnew and the phase of the new batch can be determined as follows:

T ⎧[x if k < W ⎪ new,1, x new,2 , ..., x new, k] X new = ⎨ T ⎪[x ⎩ new, k − W + 1, x new, k − W + 2 , ..., x new, k] if k ≥ W

(15)

where xnew,k is the sample of the new batch at moment k. From section 3.2, we know that the preliminary phase identification result may be a certain phase or the fuzzy region between two adjacent phases. First, we assume that the current process belongs to phase m and take it as an example to demonstrate the online local assessment strategy based on the multiple hypotheses testing. For the case of the fuzzy region between two phases, the assessment result can be simply calculated by averaging the results obtained from those two adjacent phases. The operating optimality assessment for Xnew in phase m can be formulated into the following multiple testing of W hypotheses at moment k:

(1) If k ∈ [1, t2,in), it indicates that Xnew belongs to phase 1. (2) If k ∈ [tm,in, tm−1,out] (m = 2, 3), one can determined that Xnew is running in the fuzzy region between phases m − 1 and m, i.e., maybe in phases m − 1 or m. (3) If k ∈ (tm−1,out, tm+1,in) with m = 2, 3, Xnew must belong to phase m. (4) If k ∈ (t3,out, t4,out], Xnew is within phase 4. As with any statistical modeling method, the data used for modeling are assumed to cover all normal process disturbances. Thus, any case in which the earliest start point of phase m of a new batch is earlier than that of modeling batches, tm,in, or later than the latest ending point of phase m, tm,out, will be considered as nonoptimal operating performance by the assessment procedure. Note that the above statement is just a preliminary phase identification procedure. When the new batch is identified in the fuzzy region between two phases, according to k, the phase may be determined based on the online assessment result, which will be described in detail in the next section. In addition, the durations of some phases may be very short and the fuzzy regions could contain more than two phases. However, special treatments are not required for this case, because it only increases the number of the optional phases in the fuzzy region and does not affect the phase identification procedure. Without any loss of generality, in the following sections, we assume that there are, at most, two phases in the fuzzy region and the case of more than two phases can be treated in a similar way.

Hkk− W + 1,0: x new, k − W + 1 ∈ C(m) vs Hkk− W + 1,1: x new, k − W + 1 ∉ C(m) Hkk− W + 2,0: x new, k − W + 2 ∈ C(m) vs Hkk− W + 2,1: x new, k − W + 2 ∉ C(m) ⋮ Hkk,0:

x new, k ∈ C(m) vs Hkk,1: x new, k ∉ C(m)

If too many of the above hypotheses are rejected, the current process operating performance can be declared as nonoptimal at moment k. The B−H procedure24 is used to keep the FDR of the above multiple hypotheses testing as close as possible to the target control level q, as follows: m Step 1: Compute the Mahalanobis distances Dk−W+1 , m m Dk−W+2 , ..., Dk , respectively corresponding to xnew,k−W+1, xnew,k−W+2, ..., xnew,k as follows: 6138


Article

Industrial & Engineering Chemistry Research Dlm = (x new, l − μ(m))T Σ(−m1)(x new, l − μ(m))

duration of Knew, and then the global assessment result can be given by integrating the weighted local assessment results, i.e.,

l = k − W + 1 ,..., k (16)

Step 2: Calculate the p-values of k − W + 1, ..., k: plm =

Dml ,

with respect to

Hkl,0,

K new

l=

Rg =

m l

f (x ) d x

where Rk and wk are the assessment index value and weighted coefficient at moment k, respectively. Generally, the weighted coefficient of phase m (w̃ m) represents the degree of influence and the importance of the production operation of phase m to the final product quality, and it can be determined based on the process prior knowledge. The weighted coefficient at moment k then can be defined as w̃ wk = m Nm̃

(17)

where f(x) is the probability density function of D(m),n, which is the Mahalanobis distance of x(m),n, and D(m),n = (x(m),n − μ(m))T ∑−1 (m)(x(m),n − μ(m)), n = 1, 2, ..., N(m); x(m),n and N(m) are the nth sample and the sample number of phase m. Because x(m),n follows a Gaussian distribution, the distribution of D(m),n is described as follows: D(m), n ∼

J((N(m))2 − 1) N(m)(N(m) − J )

FJ , N(m)− J

(20)

k=1

+∞

∫D

∑ wkR k

with Ñ as the number of data windows belonging to phase m. Furthermore, the weighted coefficient corresponding to the local assessment result Rk = (Rm−1 + Rmk )/2 in the fuzzy region k can be given by the mean of the weighted coefficients of phase m − 1 and m, i.e.,

(18)

where FJ,N(m)−J is an F distribution with j and (N(m) − J) degrees of freedom. Step 3: Apply the B−H procedure. (1) Rank pmk−W+1, pmk−W+2, ..., pmk in ascending order and denote them as pm(k−W+1) ≤ pm(k−W+2) ≤ , ..., ≤ pm(k). Simultaneously, let Hk(l),0 be the null hypothesis corresponding to pm(l). (2) Initially, set lmax = k − W, and find

(wm̃ − 1/Nm̃ − 1) + (wm̃ /Nm̃ ) 2 However, if there is not enough process knowledge to determine the weighted coefficient of each phase, a commonly used method is to set all of them equal. Here, assume that all the weighted coefficients are equal, i.e., wk =

⎫ ⎧ ⎛l + W − k ⎞ ⎟, k − W + 1 ≤ l ≤ k ⎬ lmax = arg max⎨l|p(ml) ≤ q·⎜ ⎝ ⎠ ⎩ ⎭ W

(19)

w1 = w2 = ... = wK new =

(3) Construct the assessment index Rmk = 1 − (lmax + W − k)/W. From the formula of the assessment index, it is easily to see that the value of Rmk is between 0 and 1. The smaller lmax is, the larger Rmk is, and few null hypotheses are rejected. In contrast, smaller Rmk values correspond to larger lmax values, and more null hypotheses do not pass the test. In order to make a distinction between optimal and nonoptimal operating performance, an assessment threshold β is introduced. If Rmk ≥ β, it means that the process operating performance is optimal at moment k; otherwise, the local assessment is nonoptimal. The value of β can be determined according to the expert experience and a value between 0.5 and 1 is usually considered appropriate. For the case of the fuzzy region between two adjacent phases, without any loss of generality, assume that the process is running in the fuzzy region between phase m − 1 and m, and m the assessment indices Rm−1 k and Rk are calculated with respect to phase m − 1 and m, respectively. If Rm−1 ≥ β, the current k process is operating on the optimal performance and belongs to phase m − 1; similarly, when Rmk ≥ β, the current process operating performance is optimal and it belongs to phase m. If both of the above two cases are not satisfied, the deterministic phase is unknown, but it can be determined that the process operating performance is nonoptimal and Rk = (Rm−1 + Rmk )/2 k is used as the final assessment index value. 3.3.2. Global Assessment. When the new batch production is over, global assessment for the entire batch should be done, which is very useful for guiding the production and operation of the future batches. From the online local assessment for the new batch, a total of Knew local assessment results are obtained for the batch with the

1 K new

and then eq 20 is further simplified as follows: Rg =

1 K new

K new

∑ Rk k=1

(21)

If Rg ≥ β, the operating performance of the entire batch is optimal; otherwise, the nonoptimal operating performance is identified, and some adjustments must be made in the next batch production. 3.4. Contribution Analysis-Based Nonoptimal Cause Identification. When the process operating performance is nonoptimal, it is important to determine the possible cause variables for the managers and operators to take appropriate operating adjustment strategy on production improvement. In this section, a variable contributions-based nonoptimal cause identification method is proposed, which is similar to the contribution plots-based fault diagnosis.5,25,26 The variable contributions to the Mahalanobis distance, as shown in eq 16, are calculated first. Moreover, consider the fact that the variables’ contributions under optimal operating performance are different from each other; it is meaningful to normalize the raw contributions to ensure that all variables give the same contributions under optimal operating performance statistically. The variables with relatively large contributions then are considered as the possible cause variables for nonoptimal operating performance. A new dataset, denoted as Xnew ′ , is constituted by the rejected samples of X new and used in the nonoptimal cause identification. Assume that the nonoptimal operating performance occurs in phase m, and then the Mahalanobis distance of Xnew ′ to the center of phase m, as in eq 16, can be rewritten as follows: 6139


Article

Industrial & Engineering Chemistry Research T −1 Dkm = (xnew ̅ − μ(m)) Σ(m)(xnew ̅ − μ(m)) 2

= ||(Σ(−m1))1/2 (xnew ̅ − μ(m))|| J

=

2 ∑ [e j(Σ(−m1))1/2 (xnew ̅ − μ(m))]

(22)

j=1

where xn̅ ew is the mean vector of X′new, which is defined as 1 J T xnew ̅ = [xnew , ..., xnew ]

and xjnew (where j = 1, 2, ..., J) is the jth variable; ej is the unit vector of the jth variable and is defined as e j = [0, ..., 0, 1, 0, ..., 0]

Thus, the raw contribution of the jth variable to the Mahalanobis distance is defined below: j j −1 1/2 2 CI raw, (xnew m = [e (Σ(m)) ̅ − μ(m))]

Figure 4. Flow diagram of the penicillin fermentation process.2 (23)

and Control Group of the Illinois Institute of Technology in 2002. This simulator can simulate the concentrations of biomass, CO2, hydrogen ion, penicillin, carbon source, oxygen, and heat generation under various operating conditions. In this study, the pH and temperature are implemented under closedloop control, while glucose addition is performed under openloop control. The process variables used for operating optimality assessment are listed in Table 1, as well as their normal operating

The ultimate variable contributions for the online data then are divided by its normal mean: CImj =

j CI raw, m

Cm̅ j

(24)

C̅ jm

where is the mean of the raw variable contributions of the jth variable of the training samples x(m),n (n = 1, 2, ..., N(m)) and is defined as N

Cm̅ j

=

∑n =(m1) [e j(Σ(−m1))1/2 (x (m), n − μ(m))]2

Table 1. Process Variables Used for Operating Optimality Assessment

N(m)

When the process runs in the fuzzy region between phase m − 1 and m, the contribution of the jth variable is equal to (CIjm−1 + CIjm)/2. Note that the larger the variable contribution, the greater the possibly that the corresponding variable has a certain type of nonoptimality in the operating process.

4. APPLICATION STUDY In this section, the proposed operating optimality assessment method is applied to the evaluation of a well-known benchmark process: the fed-batch penicillin fermentation process.36,37 4.1. Fed-Batch Penicillin Fermentation. The production of fed-batch penicillin fermentation is one of the most important chemical industrial cultivation processes and has received widespread attention for its academic and industrial importance. It is a typical multiphase batch process. From the operational point of view, the penicillin fermentation process can be divided into two physical stages, i.e., the preculture stage (about 45 h) and the fed-batch stage (about 355 h). Most of the necessary cell mass is obtained during the initial preculture stage. When the initially added substrate has been consumed by the microorganisms, the substrate feed begins. The penicillin is generated gradually from the exponential growth procedure until the stationary procedure. During this stage, the biomass growth rate must be kept constant for the purpose of optimizing penicillin production. Thus, the substrate (glucose) is supplied continuously into the fermentor, instead of being added all at once at the beginning. A flow diagram of the penicillin fermentation process is given in Figure 4. A modular simulator (PenSim v2.0) for fed-batch fermentation has been developed by Ali Ç inar et al. from the Monitoring

No.

variable

normal operating ranges

1 2 3 4 5 6 7 8 9 10 11

aeration rate agitator power substrate feed rate substrate feed temperature substrate concentration dissolved oxygen concentration biomass concentration culture volume carbon dioxide concentration pH generated heat

3−10 L/h 20−50 W 0.035−0.045 L/h 296−298 K 5−50 g/L 1.16% 0−0.2 g/L 100−200 L 0.5−1 mol/L 4−6 0 kcal/h

ranges. Under the initial operating conditions given in Table 2, 30 batches corresponding to high penicillin concentrations are produced by running the simulator repeatedly and adjusting the aeration rate, agitator power, substrate feed rate, and pH within Table 2. Initial Operating Conditions of the Penicillin Fermentation Process

6140

variable

initial conditions

substrate concentration dissolved oxygen concentration biomass concentration penicillin concentration culture volume carbon dioxide concentration pH fermentor temperature generated heat

15 g/L 1.16% 0.1 g/L 0 g/L 100 L 0.5 mol/L 5 298 K 0 kcal/h DOI: 10.1021/acs.iecr.5b04775 Ind. Eng. Chem. Res. 2016, 55, 6133−6144

Article


The process variables and the penicillin concentration of all 31 batches are shown in Figure 5, where the trajectories of each variable labeled by the red dotted line in each subfigure correspond to the nonoptimal test batch. It can be seen that, in the nonoptimal test batch, the low substrate feed rate leads to the low substrate concentration (variable 5) initially, and then biomass concentration (variable 7), culture volume (variable 8), and generated heat (variable 11) are all deduced subsequently. Known from the process mechanism, biomass growth has a high degree of dependence on both the carbon source (glucose) and oxygen as substrates. When the substrate concentration decreases, the biomass concentration will also decrease, which would, in turn, affect the utilization of substrate and the production of penicillin. Meanwhile, the declining substrate concentration and biomass concentration cause a total volume change of the culture broth and reduce the generated heat in the fed-batch process operation. 4.2. Results and Discussion. 4.2.1. Offline Phase Division and Modeling. First of all, 29 optimal training batches are unfolded into a two-dimensional matrix via variable-wise unfolding. The GMM then is established based on this matrix, and five Gaussian components are obtained by the F-J algorithm. The parameters used in offline phase division are determined according to the process knowledge and expert experience and are listed as follows: h = 5 and δ = 0.9. The offline phase division result is shown in Figure 6, from which

small ranges given in Table 3 to model the optimal operating performance. Then, 29 batches served as the training data Table 3. Optimal Ranges for Simulating Optimal Operating Performance variable

optimal ranges

aeration rate agitator power substrate feed rate pH

8.5−10 L/h 29−31 W 0.042−0.045 L/h 4.9−5.1

under optimal operating performance and one is used for online testing. In the present simulation experiment, the durations of all batches are unequal and are 390−420 h, with a sampling interval of 1 h. In addition, a test batch with low penicillin concentration is generated to simulate the nonoptimal operating performance. The substrate feed rate is adjusted to 0.035 L/h, which is in the normal operating range but beyond the optimal range. A decrease in substrate feed rate (variable 3) directly reduced the substrate concentration (variable 5) and finally resulted in a reduction in penicillin production, since glucose is the main carbon source to be fed during the fed-batch fermentation.36 Thus, the penicillin concentration at the end of nonoptimal batch is much lower than those of optimal ones.

Figure 5. Trajectories of 11 process variables and penicillin concentration of all 31 batches. 6141


Article


Figure 6. Offline phase division results of 29 optimal training batches.

one can intuitively see that the process has been devised into five phases and the durations of the training batches are unequal in each phase. Furthermore, the mean vectors and covariance matrices are updated by the data of each phase and used as the assessment model parameters. Note that the phase division result based on the proposed method is not necessarily consistent with physical stages of the penicillin fermentation process and places more emphasis on describing the characterization of local data distribution. 4.2.2. Online Phase Identification and Assessment. On the basis of the durations of each phase of each training batch, the earliest start points and the latest ending points of each phase can be determined and are tabulated in Table 4. The preliminary phase identification for the online new batch then can be implemented accordingly.

Figure 7. Optimal test batch (a) online phase identification results and (b) online assessment results.

Table 4. Earliest Start Points and Latest Ending Points of Each Phase phase No.

earliest start points

latest ending points

1

1

32

2

32

46

3

45

120

4

106

219

5

199

420

preliminary phase identification results Phase Phase Phase Phase Phase Phase Phase Phase Phase

1: k ∈ [1,32) 1 or 2: k = 32 2: k ∈ (32,45) 2 or 3: k ∈ [45,46] 3: k ∈ (46,106) 3 or 4: k ∈ [106,120] 4: k ∈ (120,199) 4 or 5: k ∈ [199,219] 5: k ∈ (219,420]

Figure 8. Nonoptimal test batch (a) online phase identification results and (b) online assessment results.

In online application, the parameters to be used are set as follows: W = 10, q = 0.05, and β = 0.8. In order to verify the effectiveness of the proposed method, it is applied to the optimal and nonoptimal test batches of the penicillin fermentation process, respectively. For the optimal test batch with a duration of 416 h, the online phase identification and assessment results are shown in Figure 7. Since the preliminary phase identification results in the fuzzy region can be corrected by the online assessment results, one can see from Figure 7a that there is no fuzzy region between any two adjacent phases. In addition, all local assessment index values are not less than 0.8 and the global assessment index value can be up to 0.9902, as shown in Figure 7b. That is to say, both the local and global assessment results show that the process operating performance is optimal and the assessment results are consistent with the actual situation. For the nonoptimal test batch with reduced substrate feed rate at the fed-batch stage, the online phase identification and assessment results are shown in Figure 8. Figure 8a shows that

the fuzzy region at the 32nd hour can be deterministic identified for the essentially optimal operating performance at this moment. Because of the nonoptimal operating performance at the fed-batch stage, there still have the fuzzy regions from phase 3 to phase 4 and phase 4 to phase 5, and the assessment indices are calculated by averaging those two adjacent phases. The online assessment results are presented in Figure 8b, where the process operating performance changes from optimal to nonoptimal at the 47th hour and the nonoptimal operating performance continues to the end of the batch production. The global assessment index value is only 0.1203. To determine the cause variables for the nonoptimal operating performance, the variable contributions to the assessment index at the 47th hour are calculated and displayed in Figure 9. By investigating the variable contributions, substrate feed rate (variable 3) and substrate concentration (variable 5) are found to be the variables responsible for 6142


Article


Figure 9. Variable constructions of the nonoptimal test batch at the 47th hour.

Figure 11. Nonoptimal test batch with 5% noise (a) online phase identification results and (b) online assessment results.

nonoptimal behavior, which is accurate and in accord with the process mechanism. In order to test the failure limitation of the proposed online assessment strategy, some uncertainties should be added in the test data. Since the uncertainty is usually random, 5% random noise is added in the optimal and nonoptimal test batches to simulate the uncertainties in the production process. The online phase identification and assessment results of the optimal and nonoptimal test batches with uncertainties are shown in Figures 10 and 11, respectively. From Figures 10 and

the process operating performance will be considered to be nonoptimal (or optimal), to avoid the interference of uncertainties to a certain extent. In addition, choosing the appropriate window width and assessment threshold value can also reduce the error evaluations, and the systemic selection scheme for those parameters is worth further study.

5. CONCLUSIONS In this study, a novel operating optimality assessment and nonoptimal cause identification method is proposed for batch processes based on multiple hypotheses testing. The multiphase and uneven-length characteristics of the batch processes are addressed in offline phase division, which makes the method more practical and more consistent with the actual situation; besides, the online assessment results are further used to correct the preliminary online phase identification results. In online operating optimal assessment, the multiple hypotheses testing technology is used for its simultaneous consideration of the operating performance, with respect to family-wise samples rather than a single sample and commit to obtain more-reliable assessment results. In addition, both the local and global evaluations have been done, from which one can master the operating performance of the current batch and control the productions of the future batches. When the nonoptimal operating performance occurs, the variable contributions to the assessment index are calculated and the variables with large contributions are considered as the responsible causes for nonoptimal operating performance. Finally, the proposed approach is applied to the fed-batch penicillin fermentation process. Both the optimal and nonoptimal batches are tested, and the online assessment and nonoptimal cause identification results are satisfactory and are in accord with the actual situations encountered with the proposed method.

Figure 10. Optimal test batch with 5% noise (a) online phase identification results and (b) online assessment results.

■

11, we can see that both the results of online phase identification and online assessment are consistent with the actual situation under the uncertainties for the two test batches. That is to say, although the production process is affected by the uncertainties, the online assessment results are still accurate and reliable. Note that, in actual production, only continuous multiple assessment indices are less (or greater) than the threshold, so

AUTHOR INFORMATION

Corresponding Author

*Tel.: +86 024 83687434. Fax: +86 024 23890912. E-mail: [email protected]. Notes

The authors declare no competing financial interest. 6143


Article


■

(19) Yoo, C. K.; Villez, K.; Lee, I. B.; Rosén, C.; Vanrolleghem, P. A. Multi-model statistical process monitoring and diagnosis of a sequencing batch reactor. Biotechnol. Bioeng. 2007, 96, 687. (20) Park, J.; Jun, C. H. A new multivariate EWMA control chart via multiple testing. J. Process Control 2015, 26, 51. (21) Li, Y. T.; Tsung, F. G. False discovery rate-adjusted charting schemes for multistage process monitoring and fault identification. Technometrics 2009, 51, 186. (22) Lee, S. H.; Jun, C. H. A process monitoring scheme controlling false discovery rate. Commun. Stat-Simul. Comput. 2012, 41, 1912. (23) Lai, T. L. Sequential multiple hypothesis testing and efficient fault detection−isolation in stochastic systems. IEEE Trans. Inf. Theory 2000, 46, 595. (24) Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. B 1995, 57, 289. (25) Qin, S. J.; Valle, S.; Piovoso, M. J. On unifying multiblock analysis with application to decentralized process monitoring. J. Chemom. 2001, 15, 715. (26) Yu, J. A new fault diagnosis method of multimode processes using Bayesian inference based Gaussian mixture contribution decomposition. Eng. Appl. Artif. Intel. 2013, 26, 456. (27) Young, D. An overview of mixture models. Stat. Survey 2008, 0, 1−24. (28) Yu, J.; Qin, S. J. Multimode process monitoring with Bayesian inference-based finite Gaussian mixture models. AIChE J. 2008, 54, 1811. (29) Figueiredo, M. A.; Jain, A. K. Unsupervised learning of finite mixture models. IEEE Trans. Pattern. Anal. Mach. Intell. 2002, 24, 381. (30) Titsias, M. K.; Likas, A. C. Shared kernel models for class conditional density estimation. IEEE Trans. Neural Networks 2001, 12, 987. (31) Dickhaus, T. Simultaneous Statistical Inference; Springer−Verlag: Berlin, Heidelberg, Germany, 2014; Chapters 10−12 (DOI: 10.1007/ 978-3-642-45182-9). (32) Hochberg, Y.; Tamhane, A. C. Multiple Comparison Procedures; Wiley: New York, 2009. (33) Storey, J. D. A direct approach to false discovery rates. J. R. Stat. Soc. B 2002, 64, 479. (34) Storey, J. D. The positive false discovery rate: A Bayesian interpretation and the q-value. Ann. Stat. 2003, 31, 2013. (35) Lehmann, E. Testing Statistical Hypotheses; Wiley: New York, 1986. (36) Birol, G.; Ü ndey, C.; Ç inar, A. A modular simulation package for fed-batch fermentation: Penicillin production. Comput. Chem. Eng. 2002, 26, 1553. (37) Birol, G.; Ü ndey, C.; Parulekar, S. J.; Ç inar, A. A morphologically structured model for penicillin production. Biotechnol. Bioeng. 2002, 77, 538.

ACKNOWLEDGMENTS We appreciate the financial support from National Natural Science Foundation of China (Nos. 61533007, 61374146, 61174130, and 61304055), State Key Laboratory of Synthetical Automation for Process Industries Fundamental Research Funds (No. 2013ZCX02-04), and Research Program of Shanghai Municipal Science and Technology Commission (No.13dz1201700).

■

REFERENCES

(1) Zhao, C. H.; Wang, F. L.; Gao, F. R.; Lu, N. Y.; Jia, M. X. Adaptive monitoring method for batch processes based on phase dissimilarity updating with limited modeling data. Ind. Eng. Chem. Res. 2007, 46, 4943. (2) Zhao, C. H.; Wang, F. L.; Mao, Z. Z.; Lu, N. Y.; Jia, M. X. Adaptive monitoring based on independent component analysis for multiphase batch processes with limited modeling data. Ind. Eng. Chem. Res. 2008, 47, 3104. (3) Wold, S.; Kettaneh, N.; Fridén, H.; Holmberg, A. Modelling and diagnostics of batch processes and analogous kinetic experiments. Chemom. Intell. Lab. Syst. 1998, 44, 331. (4) Lu, N. Y.; Gao, F. R.; Yang, Y. H.; Wang, F. PCA-based modeling and on-line monitoring strategy for uneven-length batch processes. Ind. Eng. Chem. Res. 2004, 43, 3343. (5) Qin, S. J. Survey on data-driven industrial process monitoring and diagnosis. Annu. Rev. Control 2012, 36, 220. (6) Kourti, T.; Nomikos, P.; MacGregor, J. F. Analysis, monitoring and fault diagnosis of batch processes using multiblock and multiway PLS. J. Process Control 1995, 5, 277. (7) Chen, J. H.; Liu, K. C. On-line batch process monitoring using dynamic PCA and dynamic PLS models. Chem. Eng. Sci. 2002, 57, 63. (8) Hu, Y.; Ma, H. H.; Shi, H. B. Enhanced batch process monitoring using just-in-time-learning based kernel partial least squares. Chemom. Intell. Lab. Syst. 2013, 123, 15. (9) Mori, J.; Yu, J. Quality relevant nonlinear batch process performance monitoring using a kernel based multiway non-Gaussian latent subspace projection approach. J. Process Control 2014, 24, 57. (10) Liu, Y.; Chang, Y. Q.; Wang, F. L. Online process operating performance assessment and nonoptimal cause identification for industrial processes. J. Process Control 2014, 24, 1548. (11) Liu, Y.; Chang, Y. Q.; Wang, F. L.; Ma, R. C.; Complex process operating optimality assessment and nonoptimal cause identification using modified total kernel PLS. In Proceedings of the 26th Chinese Control and Decision Conference, Changsha, China, May 31−June 2, 2014; p 1221. (12) Liu, Y.; Wang, F. L.; Chang, Y. Q. Online fuzzy assessment of operating performance and cause identification of non-optimal grades for industrial processes. Ind. Eng. Chem. Res. 2013, 52, 18022. (13) Ye, L. B.; Liu, Y. M.; Fei, Z. S.; Liang, J. Online probabilistic assessment of operating performance based on safety and optimality indices for multimode industrial processes. Ind. Eng. Chem. Res. 2009, 48, 10912. (14) Liu, Y.; Wang, F. L.; Chang, Y. Q.; Ma, R. C. Operating optimality assessment and nonoptimal cause identification for nonGaussian multimode processes with transitions. Chem. Eng. Sci. 2015, 137, 106. (15) Lu, N. Y.; Gao, F. R.; Wang, F. L. Sub-PCA Modeling and online monitoring strategy for batch processes. AIChE J. 2004, 50, 255. (16) Zhao, C. H.; Wang, F. L.; Lu, N. Y.; Jia, M. X. Stage-based softtransition multiple PCA modeling and on-line monitoring strategy for batch processes. J. Process Control 2007, 17, 728. (17) Yu, J.; Qin, S. J. Multiway Gaussian mixture model based multiphase batch process monitoring. Ind. Eng. Chem. Res. 2009, 48, 8585. (18) Zhao, C. H.; Wang, F. L.; Mao, Z. Z.; Lu, N. Y.; Jia, M. X. Improved knowledge extraction and phase-based quality prediction for batch processes. Ind. Eng. Chem. Res. 2008, 47, 825. 6144


Multiple Hypotheses Testing-Based Operating Optimality Assessment

Recommend Documents