Time-Partitioned Piecewise Affine Output Error Model for Batch

Article Cite This: Ind. Eng. Chem. Res. 2018, 57, 1560−1568

pubs.acs.org/IECR

Time-Partitioned Piecewise Affine Output Error Model for Batch Processes Zuhua Xu, Yaobo Huang, Jun Zhao,* Chunyue Song, and Zhijiang Shao National Laboratory of Industrial Control Technology, College of Control Science and Engineering, Zhejiang University, Hangzhou 310027, People’s Republic of China ABSTRACT: In this paper, a time-partitioned piecewise affine output error (PWA-OE) model for batch processes is proposed, in which time index is used to simplify the partition of regression domain by utilizing the repetitive nature of batch processes. The identification problem involves both continuous and discrete variables, and the derivative information on discrete variables is unavailable. Thus, an identification algorithm based on separable nonlinear least-squares is developed to reduce the complexity of nonlinear minimization. Furthermore, the consistent information criterion (CIC) is used to determine model structure. Results show that the time-partitioned PWA-OE model can effectively describe batch processes.

1. INTRODUCTION Batch processes have been used increasingly to manufacture of low-volume and high-value-added products, such as specialty chemicals, pharmaceutical products, and polymers.1 Compared with continuous processes, batch processes have the following unique features: absence of steady state, strong nonlinearity, limited cycle time, and repetitive nature.2 Batch processes have received increasing attention in academics and industries because of their unique properties and superior advantages. Due to the wide operating range of batch processes, process variables inherently exhibit significant nonlinearities in process dynamics. The linear time-invariant models often fail to accurately describe the dynamics of batch processes, and nonlinear identification must be employed. Among various choices, piecewise affine (PWA) models represent an attractive model structure, because they are the simplest extension of linear models and can approximate nonlinear dynamics with arbitrary accuracy.3 PWA models are a special class of nonlinear models, obtained by partitioning the state-input domain into a finite number of polyhedral regions, and by considering linear/affine subsystems in each region.4 This structure has the ability to model a wide range of nonlinear systems, while maintaining the simplicity of linear systems. The identification of the PWA model involves the estimation of the parameters of the affine submodels and the partition of state-input space defined by hyperplanes.5,6 The main difficulty lies in the fact that the identification problem includes a classification problem; that is, each data point must be associated with the most suitable submodel. In the past few years, a number of studies on PWA system identification have been reported. Ferrari-Trecate et al.7 developed an identification method that combines clustering, pattern recognition and linear identification systems. Roll et al.8 formulated two subclasses of PWA models and solved the identification © 2017 American Chemical Society

problem via mixed-integer linear or quadratic programming. Juloski et al.9 presented a Bayesian approach for identifying piecewise affine ARX (PWARX) models. Bemporad et al.10 presented a bounded-error approach for the identification of PWARX models from input−output data. Xin et al.11 formulated and solved the PWARX system identification problem under a robust EM algorithm framework. Lauer et al.12 proposed a continuous optimization framework for hybrid system identification that benefits from the advantages of algebraic and SVR approaches. Lassoued et al.13 proposed the use of Chiu’s clustering algorithm for data classification to overcome the main drawbacks of the existing methods which are the poor initialization and the presence of outliers. Lauer14 discussed complexity issues for PWA regression and showed that the global minimization of the error is NP-hard in general. Pillonetto15 proposed a hybrid stable spline algorithm for the identification of PWA systems, in which the classification variables are considered as hyperparameters and estimated via ML optimization. Breschi et al.16 developed a numerically efficient two-stage approach for PWA regression based on the combined use of recursive multimodel least-squares techniques and linear multicategory discrimination. The above-mentioned identification methods of the PWA model can be divided into two main classes. Methods in the first class alternate between solving the estimation problem for fixed model parameters and solving the estimation problem for a fixed classification, such as the clustering-based approach and Bayesian approach. These methods are prone to lead to local minima and are sensitive to initialization. In the second class, submodel Received: Revised: Accepted: Published: 1560

September 13, 2017 December 17, 2017 December 23, 2017 December 24, 2017 DOI: 10.1021/acs.iecr.7b03792 Ind. Eng. Chem. Res. 2018, 57, 1560−1568

Article

Industrial & Engineering Chemistry Research

where ek(t) is a zero-mean independent and identically distributed (iid) sequence with respect to k. This nonstationary behavior violates the basic requirement of the identification method. Hence, we difference the output and input between consecutive batch indexes, and the following time-varying model is obtained:

parameter estimation and hyperplanes estimation are simultaneously optimized, such as the mixed-integer programming approach and bounded-error approach. In this class, a combinatorial optimization problem is solved, which can be prohibitively timeconsuming. Moreover, most of the existing methods consider the PWA system as a PWARX system, which assumes that each affine subsystem is a linear ARX model. It is well-known that the ARX model has a large bias in the colored noise case. Recently, some contributions have focused on the piecewise affine output error (PWA-OE) model. The identification of the PWA-OE models with known discrete mode was implemented in Rosenqvist et al.17 Juloski et al.18 developed a Bayesian approach for the identification of the PWA-OE models. Canty et al.19 developed a clusteringbased identification algorithm for the PWA-OE system. However, batch processes can be characterized by frequent repetition, which has not been fully utilized in the existing process modeling approaches. By utilizing the repetitive nature of batch processes, Xu et al.20 first proposed the time-parametrized LTV model for batch processes, in which model coefficient is parametrized as nonlinear functions of time index. In this work, we extend this idea to the PWA model of batch processes to simplify the partition of regression domain, and propose the timepartitioned PWA-OE model of batch processes. This model structure can greatly reduce the complexity of the PWA model. The identification problem involves both continuous and discrete variables, and the derivative information on discrete variables does not exist, which makes direct minimization of the loss function extremely difficult. Subsequently, an identification algorithm based on separable nonlinear least-squares is developed by separating the estimation of OE submodels and switching points at each iteration. The remainder of this paper is organized as follows. In Section 2, a time-partitioned PWA-OE model is proposed to describe the nonlinear behavior of batch processes. In Section 3, an identification algorithm based on separable nonlinear least-squares is developed. Case studies are illustrated in Section 4, followed by conclusions drawn in Section 5.

yk̅ (t ) = G(q , t )uk̅ (t ) + ek(t )

where yk̅ (t ) = yk (t ) − yk − 1(t )

yk̅ (t ) = f (φk (t )) + ek(t )

(4)

where ⎧ φ T (t )θ φ (t ) ∈ χ 1 k 1 ⎪ k ⎪ T ⎪ φ (t )θ2 φk (t ) ∈ χ2 f (φk (t )) = ⎨ k ⎪⋮ ⎪ ⎪ φ T (t )θ φ (t ) ∈ χ ⎩ k s k s

θi = [ai1 , ai2 , ···, ain , bi1 , bi2 , ···, bin]T φk(t ) = [− yk (t − 1), − yk (t − 2), ···, − yk (t − n), uk(t − 1), uk(t − 2), ···, uk(t − n)]T

In eq 4, the following apply: f(·) is a piecewise affine map, φk(t) is the regression vector, θi is the parameter vector of each affine submodel, s is the number of submodels, and n is the order of submodel. The regions {χi}i s= 1 forms a complete partition of regression domain χ (i.e., χ = ∪is = 1χi,χi∩χj = ⌀), assumed to be a convex polyhedron, which can be described by {χi }is= 1 = {Hiφk (t ) ≤ 0}

It is well-known that the ARX model has a large bias for the OE structure.23 Moreover, the number of convex polyhedra grows with the number of data, and the PWA model is practically applicable only when the regression domain is not large. Another difficulty of the PWARX model is how to express the hyperplane such that the regions must form a complete partition of the regression domain. Different from continuous processes, batch processes have the following distinct features: a batch process (a) has finite operation duration and (b) repeats itself until a specified amount of products have been made. For batch processes, control problems are generally formulated as tracking of given fixed time-varying reference trajectories over a fixed time period. In identification for control, it is sufficient to have a model that can approximately represent the process behavior in its reference trajectory. Due to the repetitive nature of batch processes, the nonlinear behavior of batch processes around its reference trajectory can be determined by time index t. We can use the time index to simplify the partition of regression domain χ. Hence, we consider the following time-partitioned PWA-OE model, which is given by

k = 1, ..., K (1)

where t and k represent the discrete-time and cycle/batch index, respectively; T is the time duration of each cycle; K is the number of cycles; yk(t),uk(t), and vk(t) are the output, input, and disturbance at time t in the kth cycle, respectively; G(q, t) is the timevarying transfer function from uk(t) to yk(t). In comparison with continuous processes, the main characteristic of batch processes is the frequent repetition of batch runs. As a result, error trajectory due to bias in the input trajectory will repeat itself in the next batch. The same phenomenon can be observed for the disturbance, which exhibits a strong batch-wise correlation. However, the disturbance commonly exhibits a drifting characteristic along the batch index and not just random fluctuations around a stationary mean. This nonstationary behavior can be well-represented by an integrated white-noise process along the batch index:21,22 vk(t ) = vk − 1(t ) + ek(t )

uk̅ (t ) = uk(t ) − uk − 1(t )

Here, model 3 is referred to as output error (OE) model structure. The PWA model is an ideal candidate that can used for dealing with nonlinear processes with wide operating ranges. This structure has the ability to model a wide range of nonlinear systems, while maintaining the simplicity of linear systems. The traditional formulation for the PWA system is the PWARX model, which is described as follows:

2. FORMULATION OF TIME-PARTITIONED PWA-OE MODEL FOR BATCH PROCESSES For simplicity, a single-input single-output plant is considered, which repetitively performs a given task over a finite time duration, called a batch or cycle, as described by the following linear time-varying model: yk (t ) = G(q , t )uk(t ) + vk(t )t = 1, ..., T

(3)

yk̅ (t ) = wk(t ) + ek(t ) wk(t ) = f (φk (t ))

(2) 1561

(5) DOI: 10.1021/acs.iecr.7b03792 Ind. Eng. Chem. Res. 2018, 57, 1560−1568

Article

Industrial & Engineering Chemistry Research where f(·) is a piecewise affine map defined by ⎧ φ T (t )θ 1 ⎪ k ⎪ T ⎪ φ (t )θ2 f (φk (t )) = ⎨ k ⎪⋮ ⎪ ⎪ φ T (t )θ ⎩ k s

For the time-partitioned PWA-OE model (eq 7), estimated parameters Γ and Θ can be determined by minimizing the loss function:

0 < t ≤ T1 T1 < t ≤ T2

V (Γ, Θ) =

Ts − 1 < t ≤ T

θi = [ ai1 , ai2 , ..., ain , bi1 , bi2 , ..., bin ]T φk (t ) = [− wk(t − 1), − wk(t − 2), ···, − wk(t − n), uk̅ (t − 1), uk̅ (t − 2), ···, uk(t − n)]T

Denoting Bi (q) = bi1q−1 + ... + binq−n (6)

we can rewrite model 5 as ⎧ B1(q) uk̅ (t ) + ek(t ) ⎪ ⎪ A1(q) ⎪ B (q) ⎪ ⎪ 2 uk̅ (t ) + ek(t ) yk̅ (t ) = ⎨ A 2 (q) ⎪ ⎪⋮ ⎪ B (q) ⎪ s u (t ) + e (t ) k ⎪ A (q) k̅ ⎩ s

0 < t ≤ T1 T1 < t ≤ T2

Ts − 1 < t ≤ T

(7)

Vr1+ 1(Θr + 1) =

We remark that, by utilizing the repetitive nature of batch processes, we reduce the complexity of PWA model for batch processes; i.e., the number of integer variables is equal to the number of model, which allows its complexity to scale only linearly with the number of data. Moreover, this time-partitioned PWA-OE model structure can solve the scheduling problem regarding which model is switched in the online implementation of the control strategy.

Γ = [T1 , T2 , ..., Ts − 1]T

Vr2+ 1(Γr + 1) =

(8)

Bδ(t )(q)

uk(t ) = φkT (t )θδ(t ) Aδ(t )(q) ̅

(9)

where ⎧1 ⎪ ⎪2 δ(t ) = ⎨ ⎪⋮ ⎪ ⎩s

⎤

k=1

⎢⎣ T

t=1

⎥⎦

(10)

1 K

K

⎡1

T

⎤

t=1

⎥⎦

∑ ⎢ ∑ εk(t |Θr+ 1, Γr)2 ⎥ k=1

⎢⎣ T

(11)

1 K

K

⎡1

T

⎤

∑ ⎢ ∑ εk(t |Θr+ 1, Γr+ 1)2 ⎥ k=1

⎢⎣ T t = 1

⎥⎦

(12)

Step 3: Go to Step 1, and terminate the iteration when convergence occurs. We remark here that the proposed algorithm is a special case of the so-called separable nonlinear least-squares problem. The convergence of the separable nonlinear least-squares problem has been studied in Golub and Pereyra,24 and Ruhe and Wedin.25 When the predictor is bilinear or linear−nonlinear in the parameters, it is a minimization procedure that will converge to a local minimum of the original loss function. However, the proposed algorithm does not belong to the above situation, and the convergence analysis is our future ongoing work. 3.1. Parameter Estimation. This part describes the implementation of the algorithm in detail. We use the Levenberg− Marquardt method to estimate the parameter of submodels and MDS method to estimate the switching points. Moreover, the initial parameters are given by the PWARX model. 3.1.1. Parameter Estimation of Submodels. For the nonlinear least-squares problem (eq 11), there exist a large variety of numerical methods like Gauss−Newton, Levenberg− Marquardt, and Newton−Raphson. Given that the Levenberg− Marquardt algorithm can effectively alleviate the matrix illconditioned or badly scaled problem, it is used to solve this nonlinear least-squares problem. Θhr is denoted as the estimate of

The one-step-ahead optimal predictor is yk̂ (t |Θ, Γ) =

T

Step 2: Calculate parameters Γr+1 for fixed Θr+1 by minimizing the loss function

3. IDENTIFICATION OF THE TIME-PARTITIONED PWA-OE MODEL In this section, the identification of the time-partitioned PWA-OE model for batch processes is addressed. For the time-partitioned PWA-OE model (eq 7), the identification problem can be converted into the estimation of the parameters of submodels Θ and the switching points Γ, which are defined as follows: Θ = [θ1T , θ2T , ..., θsT ]T

⎡1

K

∑ ⎢ ∑ εk(t |Θ, Γ)2 ⎥

where εk (t|Θ,Γ) = yk̅ (t) − y ̂k (t|Θ,Γ) The optimization problem (eq 10) includes both continuous variables Θ and discrete variables Γ. Moreover, the derivative information on discrete variables Γ is unavailable, which renders most derivative-based methods of little or no use. Therefore, direct minimization of the loss function can be costly and may run into numerical problems. It is highly desirable to reduce the complexity by looking for some simpler numerical schemes. However, the optimization (eq 10) can become a standard identification problem for OE submodels if switching points can be fixed. If OE submodels are fixed, then it becomes a lowdimension gradient-free optimization problem. Using decoupling between OE submodels and switching points, we develop the following PWA-OE model identification algorithm based on separable nonlinear least-squares to reduce the complexity of nonlinear minimization: Initialization: Determine the initial values for switching points and submodels. Iteration: Denote Θr and Γr as the estimates at iteration r. Step 1: Calculate parameters Θr+1 for fixed Γr by minimizing the loss function

and

Ai (q) = 1 + ai1q−1 + ... + ainq−n

1 K

0 < t ≤ T1 T1 < t ≤ T2 Ts − 1 < t ≤ T 1562

DOI: 10.1021/acs.iecr.7b03792 Ind. Eng. Chem. Res. 2018, 57, 1560−1568

Article

Industrial & Engineering Chemistry Research Step 2: Rotate the simplex around the best vertex vk0

parameter vector at iteration h. Thus, we can derive the following iteration process: Θhr + 1

rik = v0k − (vik − v0k)

⎡ ⎤ = Θhr + ⎢∑ ∑ ψk(t , Θhr )ψkT(t , Θhr ) + μI ⎥ ⎢⎣ k = 1 t = 1 ⎥⎦ K

T

⎡K T ⎤ × ⎢∑ ∑ ψk(t , Θhr )εk(t |Θhr )⎥ ⎢⎣ k = 1 t = 1 ⎥⎦

−1

If f kr = min{f(rki ): i = 1,..., n} < f k0, then go to step 3; otherwise, go to step 4. Step 3: Expand the rotated simplex eik = v0k − μ(vik − v0k)

(13)

If = min{f(eki ): i = 1, ..., n} < f kr , then accept the expanded simplex vki ← eki , i = 1, ..., n; otherwise, accept the rotated simplex vki ← rki , i = 1, ..., n. Step 4: Shrink the simplex around the best vertex vk0

∂y ̂ (t |Θ) ∂εk(t |Θ) = k ∂Θ ∂Θ

sik = v0k + λ(vik − v0k)

In this work, step size μ is chosen to ensure that the current value of the criterion function is not larger than that of the previous one. To apply the above-mentioned method in estimating parameters, we also need the gradient ψk(t, Θ) of prediction. Equation 9 can be rearranged as follows:

(14)

By differentiating the equation with respect to aji and bjim respectively, we have Aδ(t )(q)

Aδ(t )(q)

∂yk̂ (t |Θ) ∂aij ∂yk̂ (t |Θ) ∂bi

j

⎧ ̂ if δ(t ) = i ⎪ − y (t − j) k =⎨ ⎪ if δ(t ) ≠ i ⎩0 ⎧ ⎪ uk ̅ (t − j) if δ(t ) = i =⎨ ⎪ if δ(t ) ≠ i ⎩0

⎧ φ T (t )θ + e (t ) 1 k ⎪ k ⎪ ⎪ φ T (t )θ2 + ek(t ) yk̅ (t ) = ⎨ k ⎪⋮ ⎪ ⎪ φ T (t )θ + e (t ) ⎩ k s k

(15)

We remark here that, for the parameter estimation of submodels, we cannot decompose it into several individual parameter estimations for submodel i using data segment [Ti−1 + 1, Ti], i.e. ⎡ ⎤ Ti 1 1 2⎥ ⎢ ε | θ V (θi) = ( t ) ∑ ∑ k i ⎥ K k = 1 ⎢⎣ T t = T + 1 ⎦ i−1

i = 1, ..., n

Accept the shrunken simplex vki ← ski , i = 1, ···, n. Step 5: Update the number of iterations k ← k + 1, go to step 1. To generate an initial simplex of problem 12 at iteration r + 1, we use the estimate value Γr at iteration r to generate a regular simplex. Typical values for λ and μ are 1/2 and 2, respectively. A stopping criterion can consist of terminating the run when the diameters of the simplex stop changing. 3.1.3. Determination of Initial Parameters. To apply a separable nonlinear least-squares method for the PWA-OE model, a good initial estimate for switching points and submodels should be provided. In this work, we consider the following PWARX model to generate the initial estimate:

yk̂ (t |Θ) + aδ1(t )yk̂ (t − 1|Θ) + ... + aδn(t )yk̂ (t − n|Θ) = bδ1(t )uk̅ (t − 1) + ... + bδn(t )uk̅ (t − n)

i = 1, ..., n

fke

where ψkT(t , Θ) = −

i = 1, ..., n

0 < t ≤ T1 T1 < t ≤ T2 Ts − 1 < t ≤ T

(17)

where φk (t ) = [− yk̅ (t − 1), − yk̅ (t − 2), ···, − yk̅ (t − n), uk̅ (t − 1),

K

uk̅ (t − 2), ···, uk̅ (t − n)]T

For the PWARX model (eq 17), the optimization problem (eq 11) in step 1 will lead to a linear least-squares problem, and the analytic solution can be derived as follows:

(16)

because εk(t|Θ) is not only related to submodel parameter θi of the current data segment [Ti−1 + 1,Ti], but also related to submodel parameter θj of the previous data segment [Tj−1 + 1, Tj] (j < i). 3.1.2. Parameter Estimation of Switching Points. The switching points of the PWA-OE model are integer variables, and the derivation information is unavailable, which requires a derivative-free method to solve the optimization problem (eq 12). The Nelder−Mead method26,27 is one of the most popular derivative-free methods. However, the Nelder−Mead method only searches in a single direction at each iteration, which can stagnate and fail to converge to a stationary point due to the deterioration of the simplex geometry. One modification of the Nelder− Mead method is a multidirectional search (MDS) algorithm proposed by Dennis and Torczon,28 which can be globally convergent to stationary points by imposing a sufficient decrease in the acceptance conditions. In this work, we adopt the MDS method to estimate switching points. The algorithm is detailed as follows: Step 0: Choose an initial simplex V0 with vertices {v00,v01,..., v0n} and 0 < λ < 1 < μ. Step 1: Order the n + 1 vertices of simplex Vk so that f(vk0) ≤ f(vk1) ≤ ... ≤ f(vkn).

Figure 1. Step response of each submodel for simulation example. 1563


Article

Industrial & Engineering Chemistry Research Table 1. Switching Points and Loss Function for n = 1 submodels

T1

2 3 4 5 6

346 93 76 93 97

T2 314 337 217 138

T3

423 237 314

T4

580 454

T5

V

531

16.3640 15.6882 15.6134 15.5564 15.3817

3.2. Model Structure Selection. In the above identification algorithm, the number of submodels and the order of OE model need to be determined. Overparameterization can lead to unnecessary moves by the controller in response to disturbances. Underparametrization can result in the lack of control action due to model inadequacy in predicting process conditions. In this work, we employ the following information criterion to select model structure29

Table 2. Switching Points and Loss Function for n = 2 submodels

T1

2 3 4 5 6

420 80 80 80 80

T2 418 210 210 210

T3

420 420 236

T4

517 420

T5

V

517

2.5416 0.5930 0.0462 0.0461 0.0460

W = N log V + γ(ρ , N )

where N = KT is the number of identification data, and d = (2n + 1)*s − 1 is the number of estimated parameters, including switching points and parameters of submodels. In eq 18, the factor γ(ρ, N) is used to penalize model structures that are too complex in view of the parsimony principle: for Akaike information criterion (AIC), γ(ρ,N) = 2p; for Bayesian information criterion (BIC), γ(ρ,N) = p log N; for consistent information criterion (CIC), γ(ρ,N) = p log2 N. It is well-known that the CIC criterion is more sensitive to model complexity, and it is adopted to determine the model structure in this work.


T1

2 3 4 5 6

420 80 80 80 80

T2 419 210 210 210

T3

420 420 326

T4

517 420

T5

V

517

2.5412 0.5910 0.0461 0.0460 0.0458

(18)

4. CASE STUDIES In this section, the performance of the proposed time-partitioned PWA-OE identification algorithm is demonstrated through two case studies. 4.1. Simulation Example. In the first case, it is assumed that the batch process is described by the following PWA linear timevarying model with sampling period ts = 10 ms.

⎡ K T1 ⎤−1⎡ K T1 ⎤ T ⎢ θ1 = ∑ ∑ φk (t )φk (t )⎥ ⎢∑ ∑ φk (t )yk̅ (t )⎥ ⎢⎣ k = 1 t = 1 ⎥⎦ ⎢⎣ k = 1 t = 1 ⎥⎦ ⎡ K T2 ⎤−1⎡ K T2 ⎤ T ⎢ θ2 = ∑ ∑ φk (t )φk (t )⎥ ⎢∑ ∑ φk (t )yk̅ (t )⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥⎦ k = 1 t = T1+ 1 k = 1 t = T1+ 1

⎧ 0.02z + 0.28 uk(t ) ⎪ 2 ⎪ z − 1.62z + 0.66 ⎪ 0.17z + 0.13 u (t ) ⎪ ⎪ z 2 − 1.85z + 0.88 k yk (t ) = ⎨ ⎪ 0.15z + 0.12 u (t ) ⎪ z 2 − 1.82z + 0.85 k ⎪ ⎪ 0.18z + 0.15 u (t ) ⎪ 2 ⎩ z − 1.72z + 0.76 k

⋮ ⎡K ⎤−1⎡ K ⎤ T T T ⎢ θs = ∑ ∑ φk (t )φk (t )⎥ ⎢∑ ∑ φk (t )yk̅ (t )⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥⎦ k = 1 t = Ts − 1+ 1 k = 1 t = Ts − 1+ 1

Hence, the identification algorithm becomes an iterative optimization, which only needs to be performed with respect to the switching points.

+ vk(t )

0 < t ≤ 80

+ vk(t )

80 < t ≤ 210

+ vk(t )

210 < t ≤ 420

+ vk(t )

420 < t ≤ 600 (19)

Figure 2. CIC criterion for simulation example. 1564


Article


Figure 3. Actual output and prediction output for simulation example.

Figure 4. Simplified diagram of injection molding machine. Figure 6. Outputs and inputs of the GBN experiments.

shown in Figure 1, which indicates time-varying dynamics at different time periods. An identification experiment is performed by a GBN signal30 with an average switch time of 5 samples. The test amplitude of the input signal is adjusted so that the noise-to-signal ratios are 8% in power. To guarantee persistent excitation of identification data, we carry out 11 batches in this experiment. We then difference the output and input between two consecutive batches to meet the identification requirements. The PWA-OE models of the number of submodels between 2 and 6, the order of submodels between 1 and 3, are estimated. Tables 1−3show the estimated switching points and loss function for different model structures. Figure 2 shows the value of CIC criterion for different model structures. As shown in Figure 2, the number of submodels and Figure 5. Injection molding machine for experiments.

Table 4. Switching Points and Loss Function for n = 1

Here, vk(t) is the unmeasurable disturbance and is described as follows: vk(t ) = v0(t ) + ek(t )

v0(t ) = 0.5 cos(0.04πt )

In this equation, ek(t) is a zero-mean Gaussian iid sequence with standard deviation 0.3. The step response of each submodel is 1565

submodels

T1

T2

2 3 4 5 6

329 148 148 147 147

329 329 217 279

T3

577 426 305

T4

502 429

T5

V

502

0.9067 0.9011 0.9002 0.8924 0.8885


Article

Industrial & Engineering Chemistry Research Table 5. Switching Points and Loss Function for n = 2 submodels

T1

2 3 4 5 6

274 294 152 152 152

T2 416 222 222 220

T3

420 441 301

T4

504 424

T5

V

504

0.2096 0.1991 0.1572 0.1273 0.1260


T1

2 3 4 5 6

297 297 152 152 168

T2 501 220 221 197

T3

500 424 223

T4

497 441

T5

V

497

0.2027 0.1917 0.1463 0.1186 0.1208

the order of submodel are chosen as 4 and 2 by minimizing the CIC criterion. We then obtain the estimated switching points (80, 210, 420) from Table 2 and the following PWA model: ⎧ 0.0240z + 0.2768 uk(t ) ⎪ 2 ⎪ z − 1.6204z + 0.6604 ⎪ 0.1666z + 0.1331 u (t ) ⎪ 2 ⎪ z − 1.8501z + 0.8801 k yk (t ) = ⎨ ⎪ 0.1479z + 0.1220 ⎪ z 2 − 1.8202z + 0.8502 uk(t ) ⎪ ⎪ 0.1783z + 0.1533 u (t ) ⎪ 2 ⎩ z − 1.7186z + 0.7587 k

Figure 8. Estimated switching points for experimental study.

4.2. Experimental Study. Injection molding is a typical batch process, which transforms thermoplastic into products of various shapes and types. The entire process consists of three main stages: filling, packing-holding, and cooling. Figure 4 shows a simplified diagram of a reciprocating-screw injection molding machine with instrumentations. Packing-holding is an important stage in an injection molding process. Adding too much material to the cavity results in highly stressed parts and mold flash. Insufficient packing-holding may cause molding problems such as poor surfaces, sink marks, weld line, and shrinkage. Hence, holding pressure is a key controlled variable during this stage. Gao et al.31,32 investigated the evolution of cavity pressure dynamics during the filling and packing stages. The cavity pressure dynamics were found to be time-varying and nonlinear. The machine used in this work is a HaiTai reciprocating-screw injection molding machine, model HTL68/JD, which is shown in Figure 5. The machine has a maximum clamping tonnage of 680 KN, and a maximum shot weight of 93 g. The sampling rate of the injection molding machine is set as 10 ms, and the entire time of the packing-holding stage is 6000 ms. The identification experiment is carried out in two steps. In the first step, a PID controller is used to obtain the nominal input and

0 < t ≤ 80 80 < t ≤ 210 210 < t ≤ 420 420 < t ≤ 600 (20)

Under the above model structure, the identification algorithm is converged in 3 iterations. To investigate the performance of the proposed method, it is compared with that of a linear OE model, in which a two-order model is identified. The prediction output of the PWA-OE model, linear OE model, and actual output of the process are shown in Figure 3. It can be clearly seen that the identification method accurately estimates the model. Furthermore, we run the identification algorithm in different noise levels with 10% and 12% to test robustness. The algorithm works well, and the result is also satisfactory in these noise levels.

Figure 7. CIC criterion for experimental study. 1566


Article

Industrial & Engineering Chemistry Research ⎧ − 0.0263z + 0.0544 uk(t ) ⎪ 2 ⎪ z − 1.7438z + 0.7736 ⎪ − 0.0021z + 0.0164 uk(t ) ⎪ 2 ⎪ z − 1.8682z + 0.8800 ⎪ − 0.0347z + 0.0651 uk(t ) yk (t ) = ⎨ 2 ⎪ z − 1.7459z + 0.7745 ⎪ − 0.0082z + 0.0262 ⎪ 2 uk(t ) ⎪ z − 1.8454z + 0.8607 ⎪ − 0.0264z + 0.0609 ⎪ u (t ) ⎩ z 2 − 1.7373z + 0.7678 k

0 < t ≤ 152 152 < t ≤ 222 222 < t ≤ 441 441 < t ≤ 504 504 < t ≤ 600 (21)

As shown in Figure 8, two extra submodels are needed to accurately describe the dynamics of the injection molding process due to the large variation of operating conditions in the ascending phase (152, 222) and descending phase (441, 504). Figure 9 plots the step responses of each submodel, which indicates time-varying and nonlinear dynamics of injection molding process. The settling time of each submodel is 0.33, 0.52, 0.34, 0.46, 0.33 s, respectively. To investigate the performance of the proposed method, it is compared with that of the two-order linear OE model. The prediction output of the PWA-OE model and the linear OE model and actual output of the process are shown in Figure 10, which illustrates that the identified PWA-OE model can effectively capture the dynamic behavior of the actual process.

Figure 9. Step responses of PWA-OE submodels.

output trajectories. In the second step, a GBN signal with an average switch time of 7 samples is superimposed on the nominal trajectories, and the resulting signal is applied to the injection molding process, as shown in Figure 6. In the identification experiment, we carry out 11 batches, and each batch contains 600 sample data points. Then, the output and input between two consecutive batches are differenced. The initial switching points are distributed uniformly among [0, 600], and the initial submodel parameters are obtained by PWARX model identification. Then, the PWA-OE models of the number of submodels between 2 and 6, the order of submodel between 1 and 3 are estimated. Tables 4− 6 show the estimated switching points and loss function for different model structures. The CIC criterion is applied on different model structures, and the results are shown in Figure 7. Hence, the number of submodels and the order of submodel are selected by minimizing the CIC criterion: s = 5 and n = 2. Under the above model structure, the PWA-OE identification algorithm is converged in 4 iterations. We then obtain the estimated switching points (152, 222, 441, 504) and the estimated PWA model as follows:

5. CONCLUSION In this work, the time-partitioned PWA-OE model for batch processes is presented, in which the time index is used to simplify the partition of regression domain and reduce the complexity of the PWA model. The identification problem involves both continuous and discrete variables, and the derivative information on discrete variables is unavailable. Hence, an identification algorithm based on separable nonlinear least-squares is developed. The proposed PWA model is also successfully tested through a simulation and a typical batch processinjection molding process. Case studies clearly demonstrate the merits of the timepartitioned PWA model and the effectiveness of the identification algorithm.

Figure 10. Actual output and prediction output for experimental study. 1567


Article


■

(19) Canty, N.; O’Mahony, T.; Cychowski, M. T. An Output Error Algorithm for Piecewise Affine System Identification. Control Eng. Pract. 2012, 20 (4), 444−452. (20) Xu, Z. H.; Zhao, J.; Yang, Y.; Shao, Z. J.; Gao, F. R. Optimal Iterative Learning Control based on Time-Parametrized LTV Model for Batch Processes. Ind. Eng. Chem. Res. 2013, 52 (18), 6182−6192. (21) Lee, J. H.; Lee, K. S.; Kim, W. C. Model-based Iterative Learning Control with a Quadratic Criterion for Time-Varying Linear Systems. Automatica 2000, 36 (5), 641−657. (22) Chin, I.; Qin, S. J.; Lee, K. S.; Cho, M. A Two-Stage Iterative Learning Control Technique Combined with Real-Time Feedback for Independent Disturbance Rejection. Automatica 2004, 40 (11), 1913− 1922. (23) Ljung, L. System Identification: Theory for the User, 2nd ed.; Prentice-Hall, Inc: NJ, 1999. (24) Golub, G. H.; Pereyra, V. The Differentiation of Pseudo-Inverse and Nonlinear Least Squares Problems Whose Variables Separate. SIAM J. Numer. Anal. 1973, 10 (2), 413−432. (25) Ruhe, A.; Wedin, P. A. Algorithms for Separable Nonlinear Least Squares Problems. SIAM Rev. 1980, 22 (3), 318−337. (26) Andrew, R. C.; Katya, S.; Luis, N. V. Introduction to Derivative-Free Optimization; SIAM: Philadelphia, 2009. (27) Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J. 1965, 7 (4), 308−313. (28) Dennis, J.E., Jr.; Torczon, V. Direct Search Methods on Parallel Machines. SIAM J. Optimiz. 1991, 1 (4), 448−474. (29) Soderstrom, T.; Stotica, P. System Identification; Prentice-Hall, Inc: London, 1989. (30) Zhu, Y. C. Multivariable System Identification for Process Control; Elsevier Science Ltd.: Oxford, 2001. (31) Gao, F. R.; Patterson, W. I.; Kamal, M. R. Cavity Pressure Dynamics and Self-Tuning Control for Filling and Packing Phases of Thermoplastics Injection Molding. Polym. Eng. Sci. 1996, 36 (9), 1272− 1285. (32) Yang, Y.; Chen, X.; Lu, N. Y.; Gao, F. R. Injection Molding Process Control, Monitoring, and Optimization; Hanser Publications: Cincinnati, 2016.

AUTHOR INFORMATION

Corresponding Author

*Tel.: +86-13003693119. E-mail: [email protected]. ORCID

Zuhua Xu: 0000-0002-7873-9521 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the NSFC-Zhejiang Joint Fund for the Integration of Industrialization and Informatization (No. U1509209), National Key Research and Development Program (No. 2017YFB0603703), and National Science Foundation of China (No. 61773340).

■

REFERENCES

(1) Korovessi, E.; Linninger, A. A. Batch Processes; CRC Press: FL, 2005. (2) Bonvin, D.; Srinivasan, B.; Hunkeler, D. Control and Optimization of Batch Processes - Improvement of Process Operation in the Production of Specialty Chemicals. IEEE Control Syst. Mag. 2006, 26 (6), 34−45. (3) Roll, J. Local and Piecewise Affine Approaches to System Identification. Ph.D. Dissertation, LinÖ ping University: Sweden, 2003. (4) Sontag, E. D. Nonlinear Regulation:The Piecewise Linear Approach. IEEE Trans. Autom. Control 1981, 26 (2), 346−358. (5) Paoletti, S.; Juloski, A. L.; Ferrari-Trecate, G.; Vidal, R. Identification of Hybrid Systems: A Tutorial. Eur. J. Control 2007, 13 (2−3), 242−260. (6) Garulli, A.; Paoletti, S.; Vicino, A. A Survey on Switched and Piecewise Affine System Identification. Proceedings of the 16th IFAC Symposium on System Identification 2012, 45, 344. (7) Ferrari-Trecate, G.; Muselli, M.; Liberati, D.; Morari, M. A Clustering Technique for the Identification of Piecewise Affine Systems. Automatica 2003, 39 (2), 205−217. (8) Roll, J.; Bemporad, A.; Ljung, L. Identification of Piecewise Affine Systems via Mixed-Integer Programming. Automatica 2004, 40 (1), 37− 50. (9) Juloski, A. L.; Weiland, S.; Heemels, W. P. M. H. A Bayesian Approach to Identification of Hybrid Systems. IEEE Trans. Autom. Control 2005, 50 (10), 1520−1533. (10) Bempora, A.; Garulli, A.; Paoletti, S.; Vicino, A. A Bounded-Error Approach to Piecewise Affine System Identification. IEEE Trans. Autom. Control 2005, 50 (10), 1567−1580. (11) Jin, X.; Huang, B. Robust Identification of Piecewise/Switching Autoregressive Exogenous Process. AIChE J. 2010, 56 (7), 1829−1844. (12) Lauer, F.; Bloch, G.; Vidal, R. A Continuous Optimization Framework for Hybrid System Identification. Automatica 2011, 47 (3), 608−613. (13) Lassoued, Z.; Abderrahim, K. An Experimental Validation of a Novel Clustering Approach to PWARX Identification. Eng. Appl. Artif. Intel. 2014, 28, 201−209. (14) Lauer, F. On the Complexity of Piecewise Affine System Identification. Automatica 2015, 62, 148−153. (15) Pillonetto, G. A New Kernel-based Approach to Hybrid System Identification. Automatica 2016, 70, 21−31. (16) Breschi, V.; Piga, D.; Bemporad, A. Piecewise Affine Regression via Recursive Multiple Least Squares and Multicategory Discrimination. Automatica 2016, 73, 155−162. (17) Rosenqvist, F.; Karlstrom, A. Realisation and Estimation of Piecewise-Linear Output-Error Models. Automatica 2005, 41 (3), 545− 551. (18) Juloski, A.; Weiland, S. A. Bayesian Approach to the Identification of Piecewise Linear Output Error Models. Proceedings of the 14th IFAC Symposium on System Identification 2006, 39, 374. 1568


Time-Partitioned Piecewise Affine Output Error Model for Batch

Recommend Documents