A multi-way LPV modeling method for batch processes

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Journal of Process Control xxx (2017) xxx–xxx

Contents lists available at ScienceDirect

Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont

A multi-way LPV modeling method for batch processes Zhonggai Zhao ∗ , Youqin Wang, Fei Liu Key Laboratory of Advanced Process Control for Light Industry (Ministry of Education), Jiangnan University, Wuxi 214122, China

a r t i c l e

i n f o

Article history: Received 13 November 2016 Received in revised form 9 October 2017 Accepted 15 October 2017 Available online xxx Keywords: Linear parameter varying (LPV) Batch process Multi-way Nonlinear model

a b s t r a c t By several simple local linear models a complex nonlinear process can be well approximated, so the linear parameter varying (LPV) method can largely simplify the complexity of model development and is widely used in the continuous process. Directly employing the regular LPV method in batch processes may develop a model for a single batch run; however, a whole batch process usually covers a number of batch runs. To develop a model for the whole batch process, this paper proposes a multi-way LPV (MLPV) approach for batch process modeling. In the proposed MLPV method, a sample is constructed by the measurements of all batch runs at each sampling instant. Comparing with the regular LPV model, the MLPV method extends the input, output and scheduling variables along the time dimension to along both the time dimension and the batch dimension. Then, the regular LPV method is performed on these “extended” samples according to the “extended” scheduling variable, by which both the variation among the batch runs and the dynamics within each batch run are summarized in the model development. The application of the MLPV method in three cases illustrates its advantages. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction In industrial processes, most of processes are characterized by multiple stages or strong nonlinearity, which brings significant challenges to the model development [1,2]. A meaningful model can be developed by fundamental modeling methods; however, it is obviously impossible to completely know the first principles of a complex process. Due to little requirement on the process knowledge, the data-driven modeling methods have become popular in recent years for the nonlinear process modeling [3–5]. Usually, the nonlinear AR(MA)X model, artificial neural-networks model, Block-oriented nonlinear model and nonlinear regression model are employed as the data-driven nonlinear model structure [6–9]. However, for a nonlinear process with multiple stages, it is difficult to find a single global nonlinear model to capture the dynamics of all operating stages. Even if the global model is available, the model structure may be too complicated to identify. To satisfy the demand on both the reduction of model complexity and the improvement of model accuracy, the linear parameter varying (LPV) model is proposed, and it has gained much attention due to its capability in approximating nonlinear processes through simple linear models with varying parameters [10]. Three types of model structures are usually employed for the LPV model, e.g. the input–output LPV (IO-

∗ Corresponding author. E-mail address: [email protected] (Z. Zhao).

LPV) based model, state-space LPV (SS-LPV) based model, and basis related functions LPV (BRFs-LPV) based model [11,12]. Because of the accurate description capability of the IO-LPV based model [12], this paper will employ the IO-LPV based model as the model structure for the batch process modeling. In the LPV modeling method, the complex nonlinear system is represented by a set of local simple linear models whose parameters are determined by other exogenous variables called the scheduling variables [13]. Weighted combination of the local models generates a global model for the nonlinear process, and these weights may be calculated according to a weighting function of the scheduling variables [14]. As a result, the model complexity is largely reduced. To overcome the influence of outliers on the modeling, Lu proposed a robust LPV method by denoting the measurements with mixture t-distributions instead of a Gaussian distribution [15]. Jin synthesized the local models to construct a global LPV model by an exponential function [16]. Banerjee employed the LPV model for state estimation, and combined the state estimates calculated from local models to obtain a global state estimate [17]. To design a process controller, Zhu interpolated the local linear models at various working-points to obtain a global model [18]. Chen investigated the parameter identification of LPV model with multiple correlated scheduling variables in the case of missing data [19]. By considering the dynamics of uncertain scheduling variables, Chen estimated the scheduling variables using the recursive Bayesian method and then the model parameters using the EM algorithm [20]. Considering scheduling variables

https://doi.org/10.1016/j.jprocont.2017.10.007 0959-1524/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model JJPC-2218; No. of Pages 12

ARTICLE IN PRESS Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

2

subjected to various constraints, Shamma analyzed the stabilization and optimal control of LPV systems [21]. Deng employed local nonlinear models instead of local linear models to describe the dynamics in operating stages, and identified the parameters using the EM algorithm with particle filter [22]. In the case that the structure of all local linear models was known a priori, Xu designed the input excitation around operating points to obtain the measurements for identification [23]. Actually, all the LPV methods mentioned above are fundamentally designed for the continuous process. Different from the continuous process, a batch process is usually operated repeatedly and a number of short-period batch runs are obtained. In a continuous process the measurements are sampled along the time direction, while in a batch process the measurements are obtained not only along the time direction but also along the batch direction. In other words, measurements sampled from a number of batch runs can be summarized in a three-way (multi-way) array, and those from a continuous process or a single batch run are characterized by a two-way formulation. Therefore, in a batch process, simple transplant of the regular LPV method can only obtain a model suitable for a single batch run and cannot cover the dynamics of the whole batch process [24,25]. In batch processes, a number of batch runs should be incorporated in the model development. However, it is not straightforward to introduce the regular LPV modeling method into batch processes, and how to deal with the multi-way data is critical to the batch process modeling. During the last decade, to monitor the batch process, MacGregor committed to introducing the conventional data-extraction methods, such as principal component analysis (PCA), partial least square (PLS), to batch processes with multi-way measurements, and proposed multi-way modeling methods for the batch process [26–28], such as the multi-way PCA (MPCA), multi-way PLS (MPLS). The MPCA or MPLS method transforms the multi-way data to a two-way array, and then performs PCA or PLS on this large array, where all the measurements within each batch run are treated as a sample, so these models describe the correlation of different batch runs. The objective of this paper is to build a LPV model for the whole batch process rather than a single batch run. Due to the multi-way feature of data in batch processes, motivated by the idea involved in MPCA or MPLS, this paper will extend the regular LPV modeling method and propose a multi-way LPV (MLPV) modeling method for batch processes. In the proposed MLPV method, a sample is constructed by the measurements of all batch runs at a sampling instant, and this “extended” sample is summarized in a matrix whose row vector and column vector represent the variable and the batch run, respectively. The MLPV method aims to develop a model for these “extended” samples, by which the variation among batch runs and time dynamics within a batch run can be captured in the model development. The application of the MLPV method in three cases demonstrates its advantages. The remainder of this paper is organized as follows: Section 2 introduces the LPV modeling method. Section 3 proposes a multi-way LPV method for the batch process modeling. Section 4 verifies the advantages of the proposed method by three example cases. Section 5 draws conclusions.

2. LPV modeling

and within each operating stage the process can be locally linearly characterized. Therefore, the nonlinear process may be separately described by multiple local models [29]. The transition process between two adjacent stages can be represented by weighted combination of local models. A set of linear models with local validity can largely simplify the model structure [29]. Consider a process as yk = f (zk , )xk + wk

where yk , xk and wk represent the output, the regressor and the noise at the kth sampling instant, respectively. zk denotes the scheduling variable, and f(zk , ) represents the function parameterized by . However, an appropriate nonlinear model structure is difficultly available to describe the nonlinear process. Considering that the nonlinear function f(zk , ) may be represented by linear models in local areas around operating points, the LPV modeling approach is proposed to approximate the process by employing multiple local linear models. Owing to the capability in approximating almost any linear dynamic system, the autoregressive exogenous (ARX) model is employed to denote these local linear models [30]. The ARX model is represented by (2) and (3) yk = xkT Â m + ek (k = 1, 2, . . ., N)

(2)

xk = [yk−1 yk−2 · · ·yk−na uTk−1 uTk−2 · · ·uTk−n 1]

T

b

(3)

where yk ∈ R1 , xk ∈ Rn+1 . ek ∈ R1 denotes the Gaussian distributed noise with zero mean and variance  2 .  m represents the mth local model parameters. uk ∈ Rt , N, na and nb are the input variables, the number of samples, orders of the output and input variables, respectively, where na + t · nb = n. The order of the ARX model can be determined by performing correlation analysis on the input variable and the output variable. Also, the heuristic method may be used for the order determination of ARX models. Without loss of generality, it is assumed that M different operating points denoted as  m (m = 1, 2, . . . M) are designed in the process. In the LPV modeling method, which local linear model does a sample belong to is determined by the scheduling variable [20]. In general, if a process variable, either an input variable or an output variable, obviously indicates the characteristics of different operating stages, it can be treated as a scheduling variable. 2.2. LPV model identification The output can be estimated by the weighted average of local models according to the scheduling variable as follows: ∧

yk =

M 

∧

˛km ykm , k = 1, 2. . .N, m = 1, 2. . .M

(4)

m=1 ∧

∧

where yk and ykm = xk T m represent the estimates from the global LPV model and from the mth local model, respectively. ˛km is the normalized weight denoting the degree that the output is compliant with the mth local model. To achieve smoothing combination of the estimates from M local models, an exponential weighting function is employed as the following equation [5]:



2.1. Linear parameter varying model It is very difficult to represent a nonlinear process over its full operating range using a deterministic mathematical model. However, the process will be operated according to a certain fixed operating trajectory. The operating trajectory may be composed of several operating stages determined by scheduling variables,

(1)

wkm = exp

−(zk − m )2 2(om )2

 (5)

where zk denotes the scheduling variable at the kth sampling instant.  m denotes the mth operating point. om represents the validity width of the mth local model.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

3

On the other hand, although the batch process is mainly characterized by repetitive feature along the batch dimension, different batch runs also may have different dynamics, so the information of only a single batch run is not sufficient to cover the dynamics of the whole batch process. Meanwhile, as a data-based method, the LPV modeling approach suffers from weak generalization capability. If only the measurements of a single batch run are used in the model development, the LPV model may only describe this single batch run rather than the whole batch process with different initial conditions or different dynamics. Therefore, measurements of a number of batch runs summarized in a three-way array should be used for the model development. However, it is intractable to directly employ the LPV modeling approach which is usually used to build the model based on a two-way array.

Fig. 1. LPV modeling method for continuous processes.

 The parameters to be estimated are written as  =

1 , 2 , · · ·M , o1 , o2 , · · ·oM . According to the likelihood principle, the parameter  is estimated as:



maximum



 = argmax P yN:1 , uN:1 , zN:1 |

 N

= argmax

M



3. MLPV modeling method 3.1. Multi-way modeling based on MPCA

˛km P yk |yk−1:1 , uk−1:1 , zk:1 , m



(6)

k=1 m=1

PCA is a valid tool for data-driven modeling, and it can extract the useful information from the measurements as follows: xk = tk P T + ek

2.3. Problem description In (2), yk and uk are stated in a vector formulation, and the LPV model aims to capture the correlation among these vectors at different operating stages along the time dimension. The LPV method for the continuous process is illustrated in Fig. 1, where the process is represented by several local linear models. However, in batch processes, a whole batch process typically covers a number of batch runs. Each batch run is similar to a continuous process, and the measurements can be arranged along the time dimension. Viewed apart, a batch run can be modeled according to Fig. 1, but for the whole batch process, due to multiple batch runs, measurements at each sampling instant are characterized in the form of matrix whose row vector and column vector represent the variable and the batch run, respectively, and all the measurements at different sampling instants can be blocked in a three-way (multi-way) matrix as illustrated in Fig. 2. All data of a batch process are expressed in the form of a three-dimensional matrix X (I × J × N), where I, N and J represent the batch number, sampling instant and variable number, respectively. On the contrary, in the continuous process, all the data are summarized into a two-dimensional matrix X (J × N), so the three-dimensional batch measurements are called multi-way data [27]. In the multi-way data, I rows correspond to I batch runs, and a row represents all measurements of a batch run. Therefore, along the time dimension, the batch process data can be arranged into N time slice matrices Xk (k = 1, 2. . .N).

(7)

where xk , tk = xk P, ek and P are the measurement, score (useful information), residual and loading matrix, respectively. P is derived by analyzing the data matrix X = [x1 x2 x3 . . . xN ]T , and the detailed PCA algorithm can refer to [26]. PCA describes the correlation among the vector xk . However, in batch processes, the measurements are multi-way formed, and a “sample” is in form of a matrix rather than a vector. Actually, it is difficult to develop a model to capture the correlation among matrices rather than vectors by PCA. To extract the useful information from the multi-way matrix, the multi-way PCA (MPCA) algorithm is proposed [26–28]. In MPCA, the multi-way array is unfolded along the time dimension to form a large twodimensional matrix, then the regular PCA method is performed on this large matrix. The multi-way measurements are transformed into a two-way array as shown in Fig. 3, where each time slice sized I × J is put side by side to the right and then a two-dimensional matrix sized I × JK is obtained. A row of this large two-dimensional matrix represents all measurements within a batch run. In MPCA, all measurements within the jth batch run X j ∈ RJK×1 are treated as the jth “measurement” and total I measurements available in a whole batch process are employed for modeling [26]. Therefore, the MPCA model describes the correlation among the “extended” vector X j as X j = T j P T + Ej

(8)

Fig. 2. Data formulation in batch processes.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model JJPC-2218; No. of Pages 12

ARTICLE IN PRESS Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

4

Fig. 4. Unfolding data in MLPV.

Fig. 3. Transformation of data in MPCA.

where T j and E j denote the “extended” score and the “extended” residual for the jth batch run, respectively. The data unit dealt with by the regular PCA algorithm is the sample at each sampling instant, and the PCA model captures the correlation along the time dimension, so at each sampling instant, the model can be employed to derive the score and residual. However, the data unit in the MPCA method is all the measurements of a whole batch run, and it captures the correlation along the batch dimension, so only after a whole batch run is complete or the future measurements are predicted at each sampling instant, the score and residual can be computed. 3.2. MLPV model The model of a batch run can be developed by the LPV method as yi,k = f (zi,k , )xi,k + wi,k

(9)

At the kth sampling instant in the ith batch run, the output yi,k is represented by the regression function of xi,k determined by the scheduling variable zi,k , and the noise is wi,k . In the continuous process, the model developed by the measurements sampled from the continuous process can describe this process; however, in batch processes, if the model denoted by (9) is developed by employing the measurements from the ith batch run, it may be only suitable for a single batch run. To make this model suitable for any batch run, only the measurements of a single batch run are obviously insufficient for the model development, and all the measurements of the whole batch process in the multi-way form may be employed for modeling. This paper proposes a multi-way LPV method for the batch process model development, where the measurements from multiple batch runs characterized by a multi-way formulation are employed for the model development. Motivated by the idea hidden in MPCA, the proposed MLPV method unfolds the multi-way data into a two-

way data before modeling. The unfolding method enables the MPCA model capable of describing correlation of all batch runs; however, the MPCA model cannot capture the correlation of different variables within batch runs. Therefore, it is intractable to employ this unfolding method to develop the LPV model in batch processes. To cover both the variation among batch runs and the dynamics within each batch run into the model development, in the proposed MLPV method, a novel unfolding method is proposed as Fig. 4, where an “extended” sample refers to a two-dimensional matrix denoting the measurements of all batch runs at a sampling instant. Therefore, for a whole batch process, the objective of the MLPV modeling method is to capture the correlation between the “extended” input and “extended” output through M local linear models. The ARX structure of the local model in batch processes is written as follows: Yk = XkT Zk + Ek (k = 1, 2, . . ., N)

(10)

T T T Xk = [Yk−1 Yk−2 · · ·Yk−na Uk−1 Uk−2 · · ·Uk−n 1] b

T

(11)

where Yk = {Y1,k , Y2,k , Y3,k , . . ., YI,k } ∈ RI×1 , Xk = {X1,k , X2,k , X3,k , . . ., XI,k } ∈ R(n+1)×I , Ek = {E1,k , E2,k , E3,k , . . ., EI,k } ∈ RI×1 denote the “extended” output, “extended” regressor and “extended” noise at the kth sample instant, respectively. Zk ∈ Rn+1 represents the parameter vector of the MLPV model. Uk = {U1,k , U2,k , U3,k , . . ., UI,k } ∈ RI×t is the “extended” input. Yi,k , Xi,k and Ui,k denote the output, regressor and input of the ith batch at the kth sample instant, respectively. na and nb are orders of the output and input, respectively, and na + t · nb = n. Without loss of generality, na = nb = 1 is assumed in this paper for simplicity. Comparing Fig. 4 with Fig. 1 indicates that the regular LPV model develops a model for vectors but the proposed MLPV method for matrices. Different from the regular LPV model, in the MLPV model, to incorporate multiple batch runs in the model development, the input and regressor are extended from a vector to a matrix, and the output is extended from a constant to a vector. As a result, the MLPV model can incorporate the dynamics within the whole batch process rather than a single batch run. In (7), the MPCA model describes the correlation of whole batch runs X j , but the proposed MLPV model denoted by (10) describes the correlation of “extended” samples. Therefore, comparing with the MPCA model, the MLPV model can capture not only the correlation among batch runs but also the dynamics within batch runs.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model JJPC-2218; No. of Pages 12

ARTICLE IN PRESS Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

5

3.4. Integration of local models In the MLPV method, the global model is the weighted sum of local models. In order to achieve smoothing integration of local models, the Gaussian distribution function is employed as the weighting function [16], and the weight of the mth local model at the kth sampling instant Wkm is calculated as follows:

 Wkm = exp

−(Zk − Qm )T (Zk − Qm )

 (12)

2(om )2

where Zk = {z1,k , z2,k , z3,k , · · ·, zI,k } ∈ RI × 1 represents the scheduling variables for all batch runs at the kth sampling instant, and zi,k denotes the scheduling variable of the ith batch run at the kth sample instant. om represents the validity width of the mth local model. The normalized weight ˛km is derived as: ˛km =

Wkm

M

m=1

(13)

Wkm

Therefore, the output estimate is expressed as:

Yˆ k =

M 

˛km Yˆ km , k = 1, 2. . .N, m = 1, 2. . .M

(14)

m=1

where Yˆ k and Yˆ km = XkT m represent the estimates by the global MLPV model and by the mth local model, respectively, and  m is the parameter of the mth  local model. The parameter  to be estimated is denoted as  = 1 , 2 , . . .M , o1 , o2 , . . .oM .

3.5. MLPV model identification According to the maximum likelihood method,  is derived as: 





= argmax P YN:1 , UN:1 , ZN:1 |

Fig. 5. Operating points in MLPV model.

= argmax

I M N 



˛km P Yk |Yk−1:1 , Uk−1:1 , Zk:1 , m

 (15)

k=1 i=1 m=1

3.3. Determination of operating points in MLPV model In the MLPV modeling method, the scheduling variable is determined similar to the regular LPV modeling method. The process variable which can indicate the characteristics of different operating stages may be selected as the scheduling variable. Differently, in the MLPV method, the scheduling variable is extended and constructed by the scheduling variables within multiple batch runs. This paper employs the K-mean clustering method to determine the operating points as Fig. 5, where the number of operating stages is denoted as M. The K-mean clustering is performed twice to obtain the operating points: First, M clustering centers (“extended” operating points) are obtained by performing K-mean clustering on the “extended” scheduling variables; Then, M operating points are obtained by performing K-mean clustering on the clustering centers. After the first clustering, the measurements of scheduling variables are divided into different clusters representing different local linear regions, and the clustering centers are constructed by the operating points of all batch runs. Considering for a whole batch process, due to the same operating recipe, the operating points of different batch runs are identical, the further clustering is performed, and M operating points q1 , q2 , . . ., qM are available. For multiple batch runs, the operating vector at the ith operating stage is Qi = [qi , qi , . . ., qi ]T ∈ RI × 1 .

where YN:1 , UN:1 and ZN:1 represent all the measurements of the “extended” output variables, the “extended” input variables and the “extended” scheduling variables over all the batch runs, respectively. The EM algorithm is a valid method for the maximum-likelihood estimation in the case of data missing or variables hidden [32,33]. In the EM algorithm, the optimization is achieved by alternate iteration of Expectation step (E-step) and Maximization step (M-step) [32]. A hidden variable Hk (k = 1, 2, . . ., N) is introduced to identify the operating stage at the kth sampling instant. Denoting measurements Yk , Uk , Zk as the data set Cobs and the hidden variable Hk as the missing data data set is constructed   set  Cmis , the complete  as C = Cobs , Cmis = Yk , Uk , Zk , Hk . In E-step, the expectation





of likelihood function for the complete data set Cobs , Cmis , called the Q function, is computed based on the observed data Cobs and the current estimate of parameter old . In M-step, the parameter is reestimated by maximizing the Q function. The continuous alternate iteration between E-step and M-step results in the maximization of the likelihood function. E-step: The Q function is computed based on the observed data Cobs and the current estimate of parameter old as

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

6



Q

|old





= EH|(old ,C

obs )

= EH|(old ,C



obs )



= EH|(old ,C

obs )







log P Cobs , H|



log(

obs )



log(



|old





obs )





· P Hk |Zk:1 , 

P Yk |Yk−1:1 , Uk:1 , Zk:1 , Hk:1 , 





i=1



· P UN:1 , ZN:1 | )

maxom ,m=1,2,...,M

=

I M N   

k=1 N

+

i=1

k=1

i=1

 

2

old

log(P(Yk |Xk , m ,  ) · P Hk = m|





old

log(P(Hk = m|Zk , om ) · P Hk = m|

(17)

, Cobs )

, Cobs )





log(P UN:1 , ZN:1 |





exp P(Yk |Xk , m ,  ) = √ 2

 −1  2 2

 

T T Yk − m Xk

T Yk − m Xk



(18)

where P(Hk = m|Z the weight. According to the  k , om ) = ˛km denotes  Bayes’ rule, P Hk = m|old , Cobs is derived as:



P Hk = m|old , Cobs







P Yk |Xk , Hk = m, m ,  P (Hk = m|Zk , om )





P Yk |Xk , Hk = m, m ,  P (Hk = m|Zk , om )

(19)

Then, the Q function is written as: |old



=

I M N    

P Hk = m|old , Cobs

k=1 N

i=1 I

m=1 M

 

P Hk = m|old , Cobs

+

k=1

+

i=1

P Hk = m|old , Cobs

i=1



N I =





· log P(Hk = m|Zk , om )

(23)

The estimation of om can be achieved by solving the constrained nonlinear optimization problem. By alternate iteration of E-step and M-step, the model parameters  m , noise variance  2 and validity width om are updated iteratively. The EM algorithm can be summarized as [32,33]:

 

(20)

· log C1

m=1











Hk = m|old , Cobs XkT Xk













(Cobs ,old ) lgL C, 

(24)

(21)

By taking the gradient of Q |old with respect to  2 and setting it to zero, the local ARX model noise variance is estimated as follows:



(25)

4. Iteration: Go back to E-step until convergence. 3.6. Online application of MLPV model According to Section 3.2, after the MLPV model is available, it may be suitable for each batch run based on (9). Especially, different from the MPCA, the proposed MLPV model can be employed to estimate the output variable at each sampling instant. The MPCA model describes the correlation among the whole batch runs; however, at each sampling instant, the current batch run is not finished. Therefore, to achieve the computation of scores, the future unobservable samples from the current sampling instant to the terminal of batch run should be compensated to make the samples complete [26–28]. Therefore, at the early part of the batch run, due to plenty compensations of samples, the estimation error of scores may be very large. In the MLPV method, at each sampling instant within a batch run, according to (12), the weighting function is written as wk,m = exp

· log P(Hk = m|Zk , om )

P Hk = m|old , Cobs XkT Yk

i=1 I P i=1

1. Initialization: Initialize parameter vector old . 2. E-step: Based on the current parameter old , the approximate Q function is computed as



M-step: The parameter is estimated by taking derivative of the Q function with respect to the parameter  m and  2 . The new parameter is derived as k=1 N k=1



S.t.omin ≤ om , m = 1, 2, · · ·, M ≤ omax

· log P(Yk |Xk , m ,  2 )

m=1

I M N    

k=1

New m



Yk −  New Xk

m=1

 = argmax Q |old



1

2



i=1



· P Hk = m|old , Cobs )



Q

P Hk = m|old , Cobs

T 

3. M-step: Derive the new iterative parameter by maximizing the Q function as

as a constant C1 = P UN:1 , ZN:1 | . P(Yk |Xk ,  m ,  2 ) indicates the probability distribution of the ARX model derived as:

m=1

Yk −  New Xk

mis |

m=1

M

i=1

P Hk = m|old , Cobs

Q |old = EC

Due to the independence of the parameters  or   the input UN:1 on the scheduling variable ZN:1 , P UN:1 , ZN:1 | is treated

=

i=1

I M N    

k=1





m=1

I M N   



· P (UN:1 , ZN:1 ))

m=1 M

 k=1

+

i=1 I



P Yk |Yk−1:1 , Uk:1 , Hk:1 , 

i=1

m=1

It is intractable to obtain an analytic solution of the local model validity width om . The mathematical optimization formulation of om is expressed as



I N 

k=1

P Hk = m|old , Cobs

N I M 

i=1

(16)

· P UN:1 , ZN:1 | )

I N 

k=1

log(

=

k=1

(22)



= EH|(old ,C



N I M

k=1

According to the correlation among Yk , Uk , Zk , Hk , the Q function can be further simplified as:



2





i=1

· P Hk |Hk−1:1 , Uk:1 , Zk:1 , 

Q



New

P Yk |Yk−1:1 , UN:1 , ZN:1 , HN:1 , 

k=1

= EH|(old ,C



I N 

· P Hk |Hk−1:1 , UN:1 , ZN:1 , 





log P YN:1 , UN:1 , ZN:1 , HN:1 |

− zi,k − qm

T 

zi,k − qm

2(om )2



(26)

where wk,m is the weight of the mth local model at the kth sampling instant. Combining (26) with (9) obtains the online estimate of output variable. Therefore, along the batch dimension, the MLPV model can be used for the description of a whole batch process containing a number of batch runs, and along the time dimension, the MLPV model can also be used online. In contrary, MPCA is essentially a modeling method along the batch runs by transforming multi-way data into two-way data, and it is difficult to deal with the correlation of variables within batch runs. 4. Simulation example In this section, the proposed MLPV method is demonstrated by three examples, e.g. a numerical example, a simulated penicillin fermentation example and an experimental coupled-tank example.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

7

0.5 Input

0.1 0

5

Output

0.2

0

−0.1

0

−0.2

5

0

−0.5 100 200 300 Time(Batch 1)

400

0

1

0.5

0.5 Input

1

0

−0.5

−1

−1 100

200

300

400

0

100 200 300 Time(Batch 1)

−5

400

10

−5

200

−15

300

0

100

0 −5

−10 100

100 200 300 Time(Batch 4)

5

0

0

0

10

5

0

−0.5

0

−5

100 200 300 Time(Batch 4)

Output

0

0

100

Fig. 6. Input data of 4 batches.

200

300

−10

400

200

300

Fig. 7. Output data of 4 batches.

4.1. Numerical example Consider a second order nonlinear batch process with the following transfer function: G(s) =

K(Z)s + 1 (Z)s2 + s + 1

(27)

where both the process gain K (Z) and the time constant  (Z) are the nonlinear functions of the scheduling variable Z expressed as follows: K (Z) = 0.6 + Z 2 , Z ∈ [1, 4] 3

 (Z) = 3 + 0.5Z , Z ∈ [1, 4]

(28)

3.5

3

2.5

2

1.5

1

0

50

100

150

200

250

300

350

(29) Fig. 8. Scheduling variables.

Err =

var(y − yˆ ) var(y)

(30)

where y is the true output and yˆ is the estimated output. The relative error of self-validation is 2.43%. The cross-validation results are shown in Fig. 10, where the estimated output can well follow its true value, and the relative error of cross-validation is 4.2%.

true output estimated output

5 Output

5

0

−5

0

100

200 Time

300

−10

400

10

0

100

200 Time

300

0

100

200

300

10 5

0

−10

−20

0 −5

Output

Obviously, within the whole operating range of the batch process, both the process gain and the time constant change dramatically, so only one linear model may hardly capture the dynamics of the whole process. Therefore, multiple local linear models are required here to sufficiently describe the behavior of the batch process. The proposed MLPV modeling method is employed. The operating points are Z1 = 1, Z2 = 2.5, Z3 = 4, respectively, and the random binary excitation signals with amplitude range [0.1,1] are imposed on the input to produce the measurements. The batch duration is set to be 400, and the output noise satisfies the Gaussian distribution N(0, 0.015). Changing the amplitude of the excitation signal among the range can obtain measurements of 10 batch runs. The input variable and the output variable of 4 batch runs are shown in Figs. 6 and 7, respectively. During the transition of two operating points, the scheduling variable Z increases with a fixed step length as shown in Fig. 8. Performing the proposed MLPV modeling on the input–output data obtains the self-validation results as shown in Fig. 9. It can be seen that the model can well capture the process dynamics under the training condition. To quantitatively evaluate the proposed model, the relative error is employed as an index as follows:

0 −5

0

100

200

300

400

−10

Fig. 9. Self-validation in Example 1.

4.2. Penicillin fermentation process Penicillin is a kind of widely used antibiotics. Its production process is characterized by nonlinearity, multiple stage and multiple batch runs. In this paper, the experimental data of penicillin fermentation process are produced from the Pensim simulation platform [31]. The platform’s software kernel is based on the Birol

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

8

Table 2 Setting and changing range of manipulated variables.

True output Estimated output 20

10

0

Control variables

Set point

Range

Ventilation rate (L/h) Stirring power (W) Bottom flow rate (L/h) Substrate temperature (K) pH Reactor temperature (K)

8.6 30 0.042 296 5 298

3–10 20–50 0.035–0.045 208–296 5–6 298–300

−10

−20

45 1.2 −30

40 0

50

100

150

200

250

300

350

1

Fig. 10. Cross-validation in Example 1.

penicillin conc(g/L)

35 30 25

Table 1 Initial conditions of Pensim platform.

part1 part2 part3 center1 center2 center3

20

Measurement variables

Range

Substrate concentration (g/L) Penicillin concentration (g/L) Biomass concentration (g/L) Dissolved oxygen concentration (mmol/L) Volume of culture medium (L) CO2 concentration pH Reactor temperature (K)

5–50 0 0–0.2 1–1.2 100–200 0.5–1 4–6 298–300

Default value

15

15 0 0.1 1.16 100 0.5 5 298

10 5 0

0

200

Yp/s

1

2 X

dt

V dt

0

800

0

200

400

600

0.5

3

0

0

100

200

300

400

500

600

700

Fig. 12. Self-validation result.

(31)

where

⎧ CS Do ⎪ ⎪ ⎨  = max Kc CX + CS Kox CX + Do ⎪ ⎪ ⎩ PP = P

part1 part2 part3

0.2

1

 PP CS dV dDo ⎪ ⎪ =− CX − CX − mo CX + Kla (Do max − Do ) − ⎪ ⎪ Yx/o Yp/o V dt dt ⎪ ⎪ ⎪ ⎪ ⎪ dV ⎪ ⎪ = F − Floss ⎪ ⎪ dt ⎪ ⎪ ⎪ ⎪ ⎩ dCO2 = ˛ dCX + ˛ C + ˛ dt

0.4

True output Estimated output

⎧ dC CX dV X ⎪ = CX − ⎪ ⎪ V dt dt ⎪ ⎪ ⎪ ⎪ CP dV dCP ⎪ ⎪ = PP CX − KCP − ⎪ ⎪ V dt dt ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ dCS = −  CX − PP CX − mx CX + F − CS dV ⎨ Yx/s

600

0.6

Fig. 11. Operating points of penicillin fermentation process.

model. The mechanism model is shown as [31]:

dt

400

0.8

CS

(32)

p

Do

p

Kp + CS + CS2 /KI Kop CX + Do

where CX , CP , CS , Do and V represent biomass concentration, penicillin concentration, substrate concentration, dissolved oxygen concentration and volume of culture medium, respectively. Initial conditions of the Pensim platform measurement variables are given in Table 1. The set point and change range of manipulated variables are given in Table 2. During the process, different batch runs may have the different initial conditions. Cold water flow rate can indicate the operating stages of the fermentation process, so it is selected as the scheduling variable, however, it is difficult to define the operating stages directly according to the continuously varying cold water flow rate. This paper employs the K-mean clustering method to obtain the clustering centers and operating points of the cold water flow rate. According to the growth rate of biomass, a penicillin fermentation process

usually experiences three stages in each batch run, e.g. slowly growing phase, quickly growing phase and the stable phase, so 3 clusters are set for the K-mean clustering. The clustering results of the cold water flow rate and the corresponding penicillin concentration are illustrated in Fig. 11. The penicillin concentration and the cold water flow rate are selected as the output variable and the input variable, respectively, and the traditional LPV model is employed to develop the process model. Setting the sampling interval as 0.5 h, 800 samples are obtained in case of the default conditions. Three operating points (the flow rate of cold water), e.g. 2.22 L/min, 26.66 L/min, 44.99 L/min, are derived. The self-validation result of the model is shown in Fig. 12, where the estimated output can tightly follow the true output. The initial substrate concentration is set as 11 g/L to produce one more batch run for the cross-validation, and the cross-validation result is shown in Fig. 13. Clearly, there is a large deviation from the 80th sampling instant to the 140th sampling instant. This period is the quickly growing stage of biomass. The model based on a single batch run cannot result in an ideal cross-validation performance during this period. The reason is that the initial condition for cross-validation is dif-

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

9

Penicillin conc(g/L)

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

True output Estimated output

1 1

0.5

1

0.5

0 0.5

1.5 Penicillin conc(g/L)

1.5

1 0

0

100

200

300

400

500

600

700

0.5 Fig. 13. Cross-validation result.

1

0.5

0

0

Fig. 15. Penicillin concentration of 4 batches.

Penicillin conc(g/L)

true output estimated output

1.5 1 0.5 0

0

200 400 600 Time(Batch 1)

1.5 1 0.5 0

800

1.5 1 0.5 Fig. 14. Cold water flow rate of 4 batches.

0

200 400 600 Time(Batch 4)

0

200

2 Penicillin conc(g/L)

2

0

0

200

400

600

1.5 1 0.5

800

0

400

600

Fig. 16. Self-validation result of the MLPV model in Example 2.

ferent from that for the model development, and the model cannot perfectly cover the dynamics of the batch run for cross-validation. From Fig. 13, the different initial conditions will lead to the large deviation of estimation at the beginning, and then the estimation results tend better. To solve the problem, the MLPV modeling method is employed. By setting the initial values of substrate concentration as 10 g/L, 12 g/L, 14 g/L, 16 g/L, 18 g/L, 20 g/L, 22 g/L, 24 g/L, 26 g/L, 28 g/L, respectively, measurements of total 10 batch runs are obtained. The input signals of 4 batches randomly chosen from the 10 batches are shown in Fig. 14 and the corresponding output measurements are shown in Fig. 15. To determine the operating points, the K-mean clustering method is performed on the “extended” scheduling variables. 3 clustering centers are obtained and the operating points are 3.58 L/h, 30.53 L/h, 46.11 L/h. The self-validation results are shown in Fig. 16, and the new initial value of substrate concentration is set as 11 g/L to generate a new batch run for the cross-validation and the results is shown in Fig. 17. In Fig. 17, since multiple batch runs with different initial conditions are incorporated in the model development, the estimates in the new batch run can quickly follow the true values at the quickly growing stage of biomass. Comparing Fig. 13 with Fig. 17 indicates that the regular LPV model in batch processes can be only used for a single batch run and the proposed MLPV model can be used for the whole batch process.

Fig. 17. Cross-validation result of the MLPV model in Example 2.

4.3. Experimental coupled-tank plant The experimental coupled-tank plant consists of a pump, a water container and two tanks, and it is illustrated by a statecoupled SISO system as described in Fig. 18. In this system, the pump transports water vertically from the water container to the upper tank, and then the water will flow

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

10

Table 4 Main parameters of the two-tank plant. Variables

Description

Value

Kp Vpmax Vppeak Dout1 Dout2 L1max Dt1 KL1 L2max Dt2 KL2 Vbias Prange Dso Dmo D1o g

Pump flow constant Pump maximum continuous voltage Pump peak voltage “Out1” orifice diameter “Out2” orifice diameter Tank1 height Tank1 inside diameter Tank1 water level sensor sensitivity Tank2 height Tank2 inside diameter Tank2 water level sensor sensitivity Tank1 and tank2 pressure sensor power bias Tank1 and tank2 pressure sensor power bias Small outflow orifice diameter Medium outflow orifice diameter Large outflow orifice diameter Gravitational constant on earth

3.3 12 22 0.635 0.47625 30 4.445 6.1 30 4.445 6.1 ±12 0-6.89 0.31750 0.47625 0.55563 981

Fig. 18. Schematic of the coupled-tank plant.

Table 3 The main symbols in the couple-tank system mathematical modeling. Symbol

Description

Kp Vp L1 L2 Do1 Do2 Fi1 Fo1 Fi2

Pump flow constant Actual pump input voltage Tank1 water level Tank2 water level Tank1 outlet diameter Tank2 outlet diameter Inflow rate of tank1 Outflow rate of tank1 Inflow rate of tank2

Fig. 19. Pump input voltage of four batches.

into the lower tank and the water container under the action of gravity through the outlet port “Out1” and the outlet port “Out2”, respectively. Two variables can be measured during the experiment by pressure sensors. Different inlet and outlet diameters in Tank1 and Tank2 may be adjusted. The mathematical model of the system is shown as:

⎧ ∂ ⎪ L1 = f1 (L1 , Vp ) ⎪ ⎪ ∂t ⎪ ⎪ ⎪ ⎨ ∂ L2 = f1 (L2 , L1 )

∂t ⎪ ⎪ ⎪ ⎪ Fi1 = Kp Vp ⎪ ⎪ ⎩

(33)

Fi2 = Fo1

The symbols in (33) and Fig. 18 are listed in Table 3, and the system parameters are illustrated in Table 4. The pump power may determine the level of Tank1 and Tank2, so the pump input voltage Vp is chosen as the input variable and the scheduling variable. The water level of Tank2 is chosen as the output variable. Three operating points are simulated by setting the initial value of Vp as 5 V, 6.5 V and 8 V, respectively. Setting the threshold of Vp as ±0.1 V, ±0.2 V, ±0.3 V and ±0.4 V obtains measurements of 4 batch runs. The measurements are shown in Figs. 19 and 20, respectively.

Fig. 20. Estimation result of output variable by LPV in four batches.

The proposed MLPV modeling method is employed and the selfvalidation results are shown in Fig. 21. In the cross-validation, the threshold of Vp is set as ±0.5 V. Both the input variable and the output variable are shown in Fig. 22, and the cross-validation results

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model

ARTICLE IN PRESS

JJPC-2218; No. of Pages 12

Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

11

5. Conclusion The regular LPV modeling method is a valid tool to describe the severely nonlinear process with several local linear models, and it is confined to the model development of continuous processes with two-way measurements. To apply the LPV method to model the nonlinear batch process, a multi-way LPV modeling method is proposed by considering the multi-way nature of data in batch processes. In the proposed MLPV method, different from the regular LPV method, the inputs, outputs and scheduling variables are extended. Both the information along the batch dimension and the dynamics along the time dimension are incorporated in the model development. As a result, different from the regular LPV model in batch processes which can only describe a single batch run, the proposed MLPV model can be used for the whole batch process including a number of batch runs. By three application cases, the proposed MLPV modeling method is proven the significant advantages over the regular LPV method. Fig. 21. Estimation result of output variable by MLPV in four batches (model selfvalidation).

Acknowledgements The authors would like to acknowledge the support by the National Natural Science Foundation of China (NSFC 61573169) and the Fundamental Research Funds for the Central Universities (JUSRP51407B). References

Fig. 22. The voltage and the tank2 water level.

true output estimated output 20

15

10

5

0

0

200

400

600

800

Fig. 23. The true and estimated water level (model cross-validation).

are shown in Fig. 23. The model is built by taking into account different batch runs with different input variables, which makes the model more adaptative. Obviously, from these figures the output estimation curve is very close to its true trajectory. These results demonstrate the effectiveness of the proposed MLPV model.

[1] M. Stosch, R. Oliveira, J. Peres, S.F. Azevedo, Hybrid semi-parametric modeling in process systems engineering: past, present and future, Comput. Chem. Eng. 60 (2014) 86–101. [2] Y. Yao, F.R. Gao, A survey on multistage/multiphase statistical modeling methods for batch processes, Ann. Rev. Control 33 (2009) 172–183. [3] J.F. MacGregor, H. Yu, M.S. Garci, J.F. Cerrillo, Data-based latent variable methods for process analysis, monitoring and control, Comput. Chem. Eng. 29 (2005) 1217–1223. [4] S.J. Qin, Survey on data-driven industrial process monitoring and diagnosis, Ann. Rev. Control 36 (2012) 220–234. [5] H. Chen, B.R. Bakshi, P.K. Goel, Toward Bayesian chemometrics – a tutorial on some recent advances, Anal. Chim. Acta 602 (2007) 1–16. [6] A.M. Shaw, F.J. Doyle III, J.S. Schwaber, A dynamic neural network approach to nonlinear process modeling, Comput. Chem. Eng. 21 (4) (1997) 371–385. [7] F. Ding, T. Chen, Identification of Hammerstein nonlinear ARMAX systems, Automatica 41 (9) (2005) 1479–1489. [8] T. Proll, M.N. Karim, Model-predictive pH control using real-time NARX approach, AIChE J. 40 (2) (1994) 269–282. [9] J. Schoukens, R. Pintelon, Y. Rolain, M. Schoukens, K. Tiels, et al., Structure discrimination in block-oriented models using linear approximations: a theoretic framework, Automatica 53 (2015) 225–234. [10] J.S. Shamma, Analysis and design of gain scheduled control systems, Ph.D. Thesis), Department of Mechanical Engineering, MIT, 1988. [11] V. Laurain, M. Gilson, R. Tth, H. Garnier, Refined instrumental variable methods for identification of LPV Box–Jenkins models, Automatica 46 (6) (2010) 959–967. [12] B. Bamieh, L. Giarre, Identification of linear parameter varying models, Int. J. Robust Nonlinear Control 12 (9) (2002) 841–853. [13] R. Murray-smith, T. Johansen, Multiple Model Approaches to Nonlinear Modelling and Control, Taylor & Francis, 1997. [14] J. Huang, G. Ji, Y. Zhu, Identification of multi-model LPV models with two scheduling variables, J. Process Control 22 (7) (2012) 1198–1208. [15] Y. Lu, B. Huang, Robust multiple-model LPV approach to nonlinear process identification using mixture t distributions, J. Process Control 24 (2014) 1472–1488. [16] X. Jin, B. Huang, D.S. Shook, Multiple model LPV approach to nonlinear process identification with EM algorithm, J. Process Control 21 (2011) 182–193. [17] A. Banerjee, Y. Arkun, B. Ogunnaike, R. Pearson, Estimation of nonlinear systems using linear multiple models, AIChE J. 43 (5) (1997) 1204–1226. [18] Y. Zhu, Z. Xu, A method of LPV model identification for control, in: 17th IFAC World Congress (IFAC’08), Seoul, Korea, July 6–11, 2008, pp. 5018–5023. [19] L. Chen, S. Khatibisepehr, B. Huang, F. Liu, Y. Ding, Nonlinear process identification in the presence of multiple correlated hidden scheduling variables with missing data, AIChE J. 61 (10) (2015) 3270–3287. [20] L. Chen, A. Tulsyan, B. Huang, Multiple model approach to nonlinear system identification with an uncertain scheduling variable using EM algorithm, J. Process Control 23 (10) (2013) 1480–1496. [21] J.S. Shamma, D. Xiong, Set-valued methods for linear parameter varying systems, Automatica 35 (1999) 1081–1089.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007

G Model JJPC-2218; No. of Pages 12 12

ARTICLE IN PRESS Z. Zhao et al. / Journal of Process Control xxx (2017) xxx–xxx

[22] J. Deng, B. Huang, Identification of nonlinear parameter varying systems with missing output data, AIChE J. 58 (11) (2012) 3454–3467. [23] Z. Xu, J. Zhao, J. Qian, Y. Zhu, Nonlinear MPC using identified LPV model, Ind. Eng. Chem. Res. 6 (48) (2009) 3043–3051. [24] Z. Zhao, B. Huang, F. Liu, Parameter estimation in batch process using EM algorithm with particle filter, Comput. Chem. Eng. 57 (41) (2013) 159–172. [25] N. Lu, Y. Yuan, F. Gao, et al., Two-dimensional dynamic PCA for batch process monitoring, AIChE J. 51 (12) (2005) 3300–3304. [26] P. Nomikos, J.F. Macgregor, Multivariate SPC charts for monitoring batch processes, Technometrics 37 (1) (1995) 41–59. [27] T. Kourti, P. Nomikos, J.E. MacGregor, Analysis: monitoring and fault diagnosis of batch processes using multiblock and multi-way PLS, J. Process Control 5 (4) (1995) 277–284.

[28] P. Nomikos, J.F. Macgregor, Multi-way partial least squares in monitoring batch processes, Chemom. Intell. Lab. Syst. 30 (1) (1995) 97–108. [29] A. Boukhris, G. Mourot, J. Ragot, Non-linear dynamic system identification: a multi-model approach, Int. J. Control 72 (7-8) (1999) 591–604. [30] L. Ljung, System Identification: Theory for the User. New Jersey, 1987. [31] G. Birol, C. Undey, A. Cinar, A modular simulation package for fed-batch fermentation: penicillin production, Comput. Chem. Eng. 26 (2002) 1553–1565. [32] L. Zhou, G. Li, Z. Song, S.J. Qin, Autoregressive dynamic latent variable models for process monitoring, IEEE Trans. Control Syst. Technol. 25 (1) (2017) 366–373. [33] A.P. Dempster, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc. 39 (1) (1977) 1–38.

Please cite this article in press as: Z. Zhao, et al., A multi-way LPV modeling method for batch processes, J. Process Control (2017), https://doi.org/10.1016/j.jprocont.2017.10.007