Priori Knowledge-Based Online Batch-to-Batch Identification in a

Article pubs.acs.org/IECR

Priori Knowledge-Based Online Batch-to-Batch Identification in a Closed Loop and an Application to Injection Molding Zhixing Cao,†,‡ Yi Yang,¶ Hui Yi,† and Furong Gao*,†,§ †

Department of Chemical & Biomolecular Engineering, Hong Kong University of Science & Technology, Clear Water Bay, Kowloon, Hong Kong ‡ Harvard John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, United States ¶ Kunda Mold (Shenzhen) Company Ltd., Shenzhen 51800, China § Fok Ying Tung Research Institute, Hong Kong University of Science & Technology, Guangzhou, China ABSTRACT: Online identification of a plant based on the input−output data obtained under closed-loop operation conditions is a fundamental problem in many industrial applications, because it paves the way for process monitoring, controller calibration, controller redesign, etc. This paper proposes a recursive batch-to-batch identification method for modeling a batch process regulated by a within-batch controller by exploiting its intrinsic repeatable pattern. To overcome severe variations of parameter estimates, three kinds of priori plant knowledge are included, namely, time constant, static gain, and stability. The constraints formulated from priori plant knowledge are efficiently handled by a technique−sequential projection. Additionally, the approach’s robustness against interbatch dynamics drift is analyzed mathematically. Finally, the performance of the proposed method is demonstrated on a numerical benchmark and a real industrial application−injection molding.

1. INTRODUCTION Online closed-loop identification is an important and fundamental problem in process control, because it depicts an evolutionary “map” to process control engineers for better understanding the plant, also enabling further controller upgrades such as redesigning controller parameters or setpoints,1 or even implementation of advanced control strategies such as model predictive control (MPC)2−9 and adaptive control.10−14 More than that, closed-loop identification extends the feasibility of the classical identification methods already established for open-loop systems. It is well-known that, in reality, there are lots of scenarios that are impossible for openloop identification. For example, an unstable process, which always has the unbounded tendency for almost any input signal, may compromise the open-loop identifiability condition; in contrast, with a stabilizing controller, closed-loop identification can achieve the modeling objective without endangering a plant’s normal operations. On the other hand, closed-loop identification is much more challenging than open-loop identification owing to the existence of the correlation between process noise and process input. Because of its importance and challenges, many papers have contributed to closed-loop identification in the past years. Roughly speaking, all the closed-loop identification methods can be categorized into three types according to the required data, namely, direct approach, indirect approach, and joint input−output approach.15 Direct approach is a natural extension of open-loop identification methods, as it formulates the regressor directly © 2016 American Chemical Society

based on the input−output data regardless of feedback mechanism. Contrarily, the indirect approach requires the knowledge of controller structure, and feeds reference signal and process outputs to the algorithm. Unlike both of them, the joint input−output approach augments process input and output as a new output and attempts to identify the augmented system. The proposed method belongs to the indirect approach. The contributions of this paper are 3-fold. The primary contribution of this article is to develop an online closed-loop identification algorithm particular for the batch process. In the classic closed-loop identification literature, based on frequency analysis, Gustavsson, Ljung, and Söderström systematically studied the closed-loop identifiability problem;15−17 Mišković and his colleagues discussed the reference excitation problem for multivariable systems;18 Hjalmarsson and Gevers focused on the experiment design of closed-loop identification;19,20 Forssell and Ljung surveyed the application of prediction error method in closed-loop identification;21 it has been reported that a two-stage method has been applied to a pilot-scale process.22 On the basis of the time-domain analytical techniques, the effects of time delay and sampling time on the identification results have been studied by Shardt and Received: Revised: Accepted: Published: 8818

May 18, 2016 July 5, 2016 July 23, 2016 July 24, 2016 DOI: 10.1021/acs.iecr.6b01900 Ind. Eng. Chem. Res. 2016, 55, 8818−8829

Article

Industrial & Engineering Chemistry Research Huang;23 Gilson and his co-workers developed a robust instrumental variable(IV) method to the initialization step,24 and further extended the method to identify a continuous-time model based on sampled data;25 Huang, Ding, and Qin used a projection subspace identification approach to eliminate the bias effect caused by the presence of feedback.26 However, identifying a batch-type plant is not a trivial task thanks to its intrinsic nonlinearity and time variation;27 generally, directly applying these methods to a batch process may thereby not yield satisfactory results, for instance, a slow response to capture intrabatch dynamics. Fortunately, the majority of batch processes are under repeated operations over a finite duration. This justifies changing the data processing and parameter estimates updating direction from time direction to batch/ iteration direction, for example, changing from t − 1 → t to k − 1 → k (t is the time index within a batch/iteration; k is the batch/iteration index); it is quite similar to iterative learning control (ILC), where control input is updated in the iteration direction, for example, uk(t) = f [uk−1(t)].28−30 To be specific, consider a simple scalar static system yk(t) = a(t)uk(t). Collecting the data in time direction gives {yk(1), uk(1), ..., yk(t), uk(t)}, in which the data contain distribution mismatch. On the contrary, organizing the data in iterations {y1(t), u1(t), ..., yk(t), uk(t)} successfully avoids the distribution mismatch. In short, processing data along the iteration direction can be approximately interpreted as handling independent identical distribution (i.i.d.) data. As to this point, similar research can be found in ref 31 in which an offline method was proposed; in refs 32 and 33 the idea was incorporated with adaptive control but without studying transient performance; in refs 34 and 35 only the open-loop cases have been discussed. Second, three different types of priori plant knowledge have been translated into constraints to refine parameter estimates. It has been argued in papers 34−36 that directly changing the data processing and estimates updating direction may give rise to undesired severe fluctuations on the parameter estimates, which somehow violates the common sense that dynamics of chemical process are slowly varying. Thus, it necessitates the estimates refine procedure. The refinements are attained by imposing a set of constraints translated from priori plant knowledge, which are time constant, static gain, and stability. The first two can be obtained from a first-order-plus-time-delay (FOPTD) model, which is quite popular in practice.37 It is not trivial to formulate the constraints on them, because the FOPTD model is in continuous time, while our identification model is in discrete time. The stability knowledge is translated into dynamic constraints, which is updated in run-to-run fashion. All these constraints are efficiently handled by the sequential projection technique. Regarding this point, relevant literature can be found: in ref 35 only the open-loop identification problem has been investigated; in ref 36 only the closed-loop stability priori knowledge has been exploited. Third, the proposed algorithm is tested on a real industrial application−injection molding. According to our knowledge, it is the first time that a linear time varying (LTV) model has been identified on injection molding. In refs 38−40 it was only reported that a linear time invariant (LTI) model has been established; in ref 41 a set-point based linear model was identified on injection molding. The remainder of this paper is organized as follows: section 2 presents the problem setup in closed-loop fashion; section 3 gives the detailed identification method; section 4 formulates the priori plant knowledge into constraints and develops the

sequential projection for these constraints; section 5 analyzes the algorithm’s robustness; section 6 studies the algorithm on a numerical example and a real industrial application−injection molding; section 7 draws a conclusion. Notations: ∥θ∥M = θ T Mθ is the Euclidean norm. The hat (•̂) means estimation or something associated with estimates. T signifies transposition. f ◦ g(x) means the function composition defined as f(g(x)). Polynomials are represented in uppercase italic symbols, i.e., A; matrices are expressed in uppercase bold symbols, i.e., A; vectors are in lowercase italic symbols, i.e., x or θ. For two positive definite matrices A, B, A ≻̲ B(A ≺̲ B) means A − B is a positive (negative) semidefinite matrix.

2. CLOSED-LOOP SYSTEM In this paper, it is assumed that a batch-type plant can be depicted by an autogressive exogenous (ARX) model, which is A(t , q−1)yk (t ) = B(t , q−1)uk(t ) + ek(t )

(1)

Here yk(t) and uk(t) are the process output and input, respectively. The subscript k is the iteration/batch index. t stands for the time index within an iteration. Denote L as the duration, and q−1 is the one-step backward shift operator in time. The disturbance and measurement noise can be represented in ek(t). Following the standard identification setup, ek(t) is assumed to be a two-dimensional white noise. That is [ei(j)em(n)] = σ 2

if and only if i = m and j = n; otherwise, it is equal to zero. Two dimensional white noise is a standard assumption in twodimensional filtering; it represents the worst scenario that the noise has effects on the entire frequency domain. A(t, q−1) and B(t, q−1) are time-varying polynomials in the form of A(t , q−1) = 1 + a1,0(t )q−1 + a 2,0(t )q−2 + ... + ana ,0(t )q−na

(2a)

B(t , q−1) = b1,0(t )q−1 + b2,0(t )q−2 + ... + bnb ,0(t )q−nb (2b)

where na and nb are the orders of the three polynomials 2a−2b. The coefficients a1,0, ..., ana,0 and b1,0, ..., bnb,0 : [0, L] →  are time-varying real functions. Hereby we have assumed that the delay order is equal to 1; in fact, this assumption does not lose any generality, because it can be addressed by letting the first d − 1 terms in eq 2b be equal to zeros. In addition, the time variations arise from the repeated operations in the proximity of a certain operation profile, which may roughly make a nonlinear system reduce to an LTV system. To be specific, consider a system y = G(t, u) regulated to follow an operating line ur, which is a function of time t. Under a relatively good regulation, the closed-loop dynamics can be approximated by a linearized model y − yr = G(t , u) − G(t , ur ) ≈

∂G |u (u − ur ) ∂u r

∂G

where ∂u |ur is a function of time t. This justifies the LTV ARX model assumption with a relatively good controller. 8819

DOI: 10.1021/acs.iecr.6b01900 Ind. Eng. Chem. Res. 2016, 55, 8818−8829

Article

Industrial & Engineering Chemistry Research

Figure 1. Predictor structure.

As for the controller, it assumes that the plant (1) is regulated under an R−S−T controller, a two-degree-of-freedom controller, whose explicit form is given as uk(t ) = −

R(q−1) S(q−1)

yk (t ) +

T (q−1) S(q−1)

rk(t )

yk (t ) = −A*(t )yk (t ) + B(t )uk(t ) + ek(t )

where A(t) = 1 + A*(t). A*(t) contains the terms with at least one unit backward shift. It can also be written in linear regression form:

(3)

yk (t ) = ϕkT(t )θ0 + ek(t )

It is noted that the PID controller can be easily accommodated in this type. Actually, the controller is quite general, and has become a standard assumption in closed-loop identification.42 The rk(t) in eq 3 is the reference signal; it comprises two parts: one is the repeatable part (the operation profile rr(t)), the other is the nonrepeatable part to provide sufficient excitation, rnr,k(t). To be specific, it is rk(t ) = rr(t ) + rnr, k(t )

ϕkT(t ) = [−yk (t − 1) ...

(5a)

...

− yk (t − na)

uk(t − 1)

uk(t − nb)]

θ0T(t ) = [a1,0(t )

(5b)

...

ana,0(t )

b1,0(t )

...

bnb,0(t )] (5c)

Given the parameter estimates, the output prediction ŷk(t) is generated by

(4)

Furthermore, it is assumed that S(q−1) is in the canonical form whose constant term is equal to 1, that is, S(q−1) = 1 + q−1S*(q−1). It does not lose any generality; it is convenient to transform noncanonical into canonical form by simple scaling. Thereafter, the control law can be summarized as

yk̂ (t ) = φkT(t )θk̂ (t ) φkT(t ) = [−yk̂ (t − 1) ...

uk(t ) = f (yk (t ))

T

...

− yk̂ (t − na)

uk̂ (t − 1)

uk̂ (t − nb)]

θk̂ (t ) = [a1,̂ k (t )

a function of yk(t).

(6a)

(6b)

...

ana ̂ , k (t )

b1,̂ k (t )

...

bnb̂ , k (t )] (6c)

where ûk(t) is the control input based on the prediction ŷk(t), namely ûk(t) = f(ŷk(t)) = −Rŷk(t)/S + Trk(t)/S. Note that, basically, the predictor in eq 6 proposed by Landau et al. replicates the plant dynamics but based on the predictions and parameter estimates only, shown in Figure 1. Define the output error

3. IDENTIFICATION METHODS With the basics defined, we are ready to present the identification method. Basically, it includes two components: predictor and estimator. The predictor is used to generate output prediction, that is, ŷk(t), while the estimator is used to yield parameter estimates by minimizing the output prediction error, that is, ŷk(t) − yk(t). 3.1. Predictor. To keep the notations neat, the shift operator q−1 will be dropped out in the remainder of the paper without any ambiguity. From eqs 1 and 3, the closed-loop dynamics of the plant is

ϵk (t ) ≔ yk (t ) − yk̂ (t )

Then, one rewrites eq 5a as yk (t ) = φkT(t )θ0(t ) + [ϕk (t ) − φk (t )]T θ0 8820


Article


identification algorithm, such parameter estimate fluctuations are not often seen. Even in mathematics, it can be proven that the traditional RLS has an inherent smoothing effect; that is ∥θ̂k(t) − θ̂k(t − 1)∥ → 0 even ∥θ̂k(t + i) − θ̂k(t)∥ → 0 for an arbitrary positive integer i as t → ∞.43 The possible reason responsible for these undesired fluctuations is that there do not exist strong relations between the data in different iterations k, for instance, φk(t) and φk−1(t), due to the i.i.d. sampling. But in terms of time t, the regressors at different times are governed by some difference equation. For example, due to the existence of plant dynamics, there must exist a transition matrix ((t ) and an input vector ηk(t) such that

Since −yk(t) + ŷk(t) = −ϵk(t) and uk(t) − ûk(t) = −Rϵk(t)/S, it can be further simplified to yk (t ) = φkT(t )θ0(t ) − [A*(t ) + B(t )R /S]ϵk (t )

Combining with eq 6, one can get ϵk (t ) = φkT(t )[θ0(t ) − θk̂ (t )] + ek(t ) − [A*(t ) + B(t )R /S]ϵk (t )

It is equivalent to H(t )ϵk (t ) = S[−φkT(t )θk̃ (t ) + ek(t )]

(7)

where H(t) ≔ A(t)S + B(t)R is the closed-loop characteristic polynomial, which decides the closed-loop stability. The minus sign arises from the definition θ̃k(t) ≔ θ̂k(t) − θ0. To justify the predictor, letting θ̂k(t) = θ0, which suggests that for a bounded noise ek(t), the output prediction error ϵk(t) will converge to ek(t) as t → ∞. 3.2. Estimator. It is known that identification problems are usually translated into an empirical optimization problem such as minimizing the output prediction error. Therefore, how to formulate the associated cost function is of great importance. In the conventional identification methods, the cost function is selected as the summation of the output prediction error along the time t. This choice is not wise in the case of a time-varying plant, since all the data {yk(1), φk(1), ..., yk(t), φk(t)} are drawn from different distributions caused by a time variation of the plant. That may give rise to the failure of converging to the true model, that is, θ0. This problem is easily circumvented by switching the summation direction from time t to iteration k. That means instead of solving a single optimization problem, we solve a set of optimization problems differentiated by time t; and at each time t, the data {y1(t), φ1(t), ..., yk(t), φk(t)} can be approximately regarded as i.i.d. drawn from the same distribution. Specifically, at each time t, we attempt to minimize the cost function

φk (t ) = ((t )φk (t − 1) + ηk (t )

which in nature is a low-pass filter with smoothing effect. With this understanding, an intuitive solution is to impose some sensible constraints to restrict the estimate fluctuations. In this section, we are going to formulate these constraints based on priori plant knowledge. Nevertheless, practically, it may lack computationally efficient methods to solve the constrained optimization problem, therefore making online identification infeasible. Since new measurements are continuously coming, it is more important to find an efficient suboptimal solution to accommodate the new measurements than seek an optimal solution at each time. Projecting an unconstrained optimizer onto a closed feasible domain is a plausible approach for suboptimal solution seeking. Mathematically speaking, given a feasible domain Θ constructed from priori plant knowledge, the projection step is to find ProjΘ [θk̂ (t )] ∈ arg min ∥θ − θk̂ (t )∥ θ ∈Θ

This step is easy to implement if the domain Θ is geometrically simple, for example, ball, polygon, cube, etc. Section 4.1 formulates three types of priori plant knowledge into various constraints. Section 4.2 proposes an efficient approach to solve the constrained optimization problem based on sequential projection. 4.1. Priori Knowledge Formulation. In this subsection, we focus on constraints formulation of priori plant knowledge. Three types of priori plant knowledge will be discussed, namely, time constant, static gain, and stability. The first two types are the key elements of a FOPTD transfer function model; that is,

k

Jk (t ) =

∑ [yi (t ) − φiT(t )θ ]̂ 2 i=1

Invoking Newton’s methods,17,34 it is easy to arrive at θk̂ (t ) = θk̂ − 1(t ) + Kk(t )[yk (t ) − φkT(t )θk̂ − 1(t )] K k (t ) =

(8a)

Pk − 1(t )φk (t ) 1 + φkT(t )Pk − 1(t )φk (t )

Pk(t ) = Pk − 1(t ) −

G (s ) =

(8b)

Pk − 1(t )φk (t )φkT(t )Pk − 1(t ) 1 + φkT(t )Pk − 1(t )φk (t )

K e−τds τs + 1

Here K, τd, and τ are the static gain, time delay, and time constant, respectively. The FOPTD model plays a crucial role in process control, since it provides important guidance for PID or IMC (internal model controller) tuning.37 Note that our plant has been already operated in closed loop; this implies that an FOPTD model is probably ready for use. Hence, it is not risky or restrictive to assume its existence. It should be remarked that the FOPTD model is a continuous-time model, whereas the ARX model is in discrete time. It is not trivial to translate the knowledge of time constant and static gain into constraints imposed on the estimates of the ARX model. On the other hand, process stability widely exists in industry, thanks to various physicochemical conservation laws, which make the assumption on it not lose generality either. Unlike time constant or static gain, the stability constraint can be

(8c)

It is similar to eq 15 and 16 in ref 34, but here we use a different regressor φk(t) which only relies on the past predictions ŷk(t), ûk(t), rather than past outputs yk(t), uk(t).

4. ESTIMATE REFINEMENT As argued in refs 34−36 a direct employment of eq 8 cannot yield a “smooth” parameter estimate that may result in an unsatisfactory real dynamics tracking but still with good output prediction capabilities. In practice, such model fluctuations may give rise to frequent controller action for some model-based control policies, thus perhaps degrading control performance. However, in the traditional recursive least-squares (RLS) 8821


Article

Industrial & Engineering Chemistry Research na ⎧ nb ⎪ ∑ bî , k (t ) − 2K[1 + ∑ aî , k (t )] ≤ 0 ⎪ ⎪ i=1 i=1 *2: ⎨ nb na ⎪ ̂ − + b t K ( ) 0.5 [1 ⎪∑ i , k ∑ aî ,k(t )] ≥ 0 ⎪ ⎩ i=1 i=1

established in discrete time with the help of the Lyapunov equation. 4.1.1. Time Constant. It is well-known that the FOPTD model is usually a reduced-order model approximating the original high-order dynamics. There are many methods contributing to this topic, such as Taylor series approximation and Skogestad’s half rule.37,44 Assume that the original plant has the dynamics represented by

4.1.3. Stability. In contrast to the aforementioned two kinds of constraints, the stability constraint will be given as a norm ball whose size is updated in real time. This kind of constraints has been studied in Cao et al.35 To be concise, only the backbone will be presented; for more details, please refer to ref 35. For every parameter estimate θ̂k(t), there exists a corresponding state matrix; that is,

F (s )

G (s ) =

na ∏i = 1 (τis

+ 1)

where F(s) is a polynomial of s. The sequence {τi} is a nonincreasing sequence, that is, τ1 ≥ τ2 ≥ ··· ≥ τna. Note that due to exact knowledge of process order, the order na of ARX model has to be consistent with the number of terms in the product of G(s). Decomposing G(s) into several first-order terms gives na

G (s ) =

∑ i=1

θk̂ (t ) ∼ Â k (t ) ⎡−a1,̂ k (t ) ⎢ ⎢ 1 =⎢ 0 ⎢ ⎢ ⋮ ⎢⎣ 0

fi s + 1/τi

for some f i. Discretizing G(s), one can have na

G(z −1) =

∑ i=1

fi 1−e

−T / τi

=

Suppose that θ̂k−1(t) satisfies the stability constraint; to be precise, ρmax [Â k−1(t)] < η < 1, where η is the priori modulus upper bound of the largest eigenvalue. According to the Lyapunov theorem, for identify matrix, there exists a positive definite matrix Mk(t) such that

F ′(z) na ∏i = 1 (1

− e−T / τi)

for some F′(z). Therefore, the following relation can be established.

T Â k − 1(t )Mk(t )Â k − 1(t ) − η2 Mk(t ) = −I

i=1

[Â k − 1(t ) + δ Â k (t )]Mk(t )[Â k − 1(t ) + δ Â k (t )]T (t )

where T is the sampling time. It remains to find appropriate upper and lower bounds of the summation ∑ina= 1 e−T/τi. Notice that the function g(x) = exp(−T/x) is a monotonically increasing function and τ ≥ τ1. Therefore, a reasonable lower bound for a1(t) is − (na)e−T/τ. On the other hand, according to Skogestad’s half rule,44 the lower bound of the largest time constant τ2 is

− η 2 M k (t ) ≺ 0

ϵ−1η2 Mk(t ) − (1 + ϵ−1)I T + (1 + ϵ)δ Â k (t )Mk(t )δ Â k (t ) ≺ 0

Note that there exists a relation that δÂ k(t)Mk(t)δÂ k = diag{∥δθ̂A,k(t)∥2Mk(t), 0(na−1)×(na−1)}, where δθ̂A,k(t) is the first na elements of δθ̂ k (t). Based on that, maximizing 2 ∥δθ̂A,k(t)∥M , one can have k(t)

Hence, a heuristic upper bound for a1(t) is − (na)e−3T/(2τ). By now, we have obtained an interval for â1,k(t), which is *1: a1,̂ k (t ) ∈ [−(na)e−T / τ , −(na)e−3T /(2τ)]

Δk (t ) ≔ max ∥δθÂ , k(t )∥2Mk(t ) =

4.1.2. Static gain. Recall the definition of static gain, that is

In this way, the stability constraints have been translated into a weighted 2-norm ball whose center is the previous parameter estimate θ̂k−1(t). The size of the ball is controlled by eq 11. To summarize, the constraints can be written in the following compact form:

where y ̅ is the steady-state output under constant input u̅. According to eq 1, we have

i=1

1 2η2λmax [Mk(t )] − 1 + 2η (η2λmax [Mk(t )] − 1)λmax [Mk(t )]

(11)

y ̅ = Ku ̅

na

(10)

with δÂ k(t) ≔ Â k(t) − Â k−1(t). It is noticed that it is a standard setup in robust control theory. By robust parametrization, the sufficient condition for eq 10 is

2 τ1 + τ2/2 = τ < τ1 + τ1/2 ⇔ τ1 < τ 3

∑ ai)y ̅

(9)

If we want to ensure θ̂k(t) satisfying stability constraint, or equivalently, Â k(t) is stable, it has to guarantee that

na

a1(t ) = −∑ e−T / τi

(1 +

... −ana ̂ − 2, k (t ) −ana ̂ − 1, k (t ) −ana ̂ , k (t )⎤ ⎥ ⎥ 0 ... 0 0 ⎥ 1 ... 0 0 ⎥ 0 ⋱ ⋮ ⋮ ⎥ ⎥⎦ 0 ... 1 0

nb

= (∑ bi)u ̅

*3: θk̂ (t ) ∈ {θ|∥θA , k − θÂ , k − 1(t )∥2Mk(t ) ≤ Δk (t )}

i=1

na It indicates that K = (∑i nb = 1bi)/(1 + ∑i = 1ai). Given 100% variation on K, the heuristic bound can be given as

Remark 1: Note that the constraints *1, *2 are static constraints so that they can be computed offline; the constraint 8822


Article


*3 is a dynamic constraint so that it needs updating; that makes it a standard trust region method,45 where Δk(t) is the size of trust region. In the updating procedure, it is not computationally cheap to solve Mk(t) in the Lyapunov equation; however, the computation of Mk(t) and Δk(t) can be implemented in parallel with the identification procedure (also can be calculated in the resettling time of batch process). Remark 2: It is worthy remark that FOPTD model is only used for estimate refinement not the identification objective to achieve. Additionally, we do not impose any constraint on the model dimension; it does not have to be the first-order model. Admittedly, it is relatively easy to solve the optimization problem with the time constant and static gain constraints by the means of a direct optimization method; but the stability constraint is not computational friendly for direct optimization in real time, because it is nonconvex as argued in our previous paper.35 This can be seen from eq 9, which contains a product of two optimization variables Â k−1(t) and Mk(t). Thus, it motivates the usage of projection method detailed in the following section. 4.2. Sequential Projection. In the above subsection, the priori plant knowledge on time constant, static gain, and stability has been turned into constraints, namely *1, *2, *3. Further note that *1 delineates an interval on real axis denoted as Θ *1; *2 describes a space supported by two hyperplane

constraints have been formulated with stronger confidence so that they have to be satisfied more strictly. In our case, the stability constraint is the most confident one, because there exists a step response; that makes it the most outside one. Now it remains to figure out each projection step in eq 12. Then, for *1, we have for any x ∈ na + nb ⎧[− (na)e−T / τ , x T]T , x < − (na)e−T / τ a1̅ a1 ⎪ ⎪ ProjΘ (x) = ⎨ x , − (na)e−T / τ ≤ xa1 < − (na)e−3T /(2τ) *1 ⎪ ⎪[− (na)e−3T /(2τ), x T]T , x ≤ − (na)e−3T /(2τ) a1̅ a1 ⎩ (13)

where xa1 stands for the first element of x and xa¯1 means the last na + nb − 1 elements of x. Note *3 imposes a constraint on the first na elements of θ̂k(t) ; *2 can be used to impose on the last nb elements of θ̂k(t). ⎧[x T , x T ]T , x < 0.5K (1 + x ) b+ a+ ⎪ a b★ ⎪ ProjΘ (x) = ⎨ x , 0.5K (1 + xa +) ≤ xb + < 2K (1 + xa +) *2 ⎪ ⎪[x T , x T ]T , x ≥ 2K (1 + x ) ⎩ a b* b+ a+ (14)

where xa+, xb+ mean the summations of all the first na and all the last nb elements of x respectively. xa stands for the vector that contains the first na elements of x. xb★ and xb* are defined as

denoted as Θ *2 ; *3 depicts a norm ball in space na + nb, denoted as Θ *3. Hence, the projection step turns to be ProjΘ

*1 ∩Θ *2 ∩Θ *3

[θk̂ (t )] = arg

min

θ ∈Θ *1∩Θ *2 ∩Θ *3

∥θ − θk̂ (t )∥

Although every space Θ *1, Θ *2 or Θ *3 is geometrically simple enough, their intersection can be quite complex, perhaps making the projection step computationally intractable. Fortunately, recall that the primary goal of our problem is to use proper constraints to mitigate the estimates variations efficiently rather than to attain the optimality of the corresponding projection step. Bearing that in mind, an intuitive efficient approximation of the aforementioned projection operator is

and xb = *

ProjΘz (x) = [za + αk(t )(xa − za), xbT]T

(15a)

⎧1, ∥xa − za∥Mk(t ) ≤ Δk (t ) ⎪ ⎪ αk(t ) = ⎨ Δk (t ) , otherwise ⎪ ⎪ ∥xa − za∥M (t ) ⎩ k

(15b)

*3

≔ ProjΘ ◦ ProjΘ ◦ ProjΘ [θk̂ (t )] ≈ ProjΘ

2K (1 + xa +) xb xb +

xb is the vector comprising the last nb elements of x. As for the projection of *3, it is a projection on to a norm ball. So we particularly introduce the norm ball projection operator as

R θk̂ (t ) ≔ Proj[θk̂ (t )] *3

0.5K (1 + xa +) xb xb +

xb★ =

*2

*1 ∩Θ *2 ∩Θ *3

[θk̂ (t )]

*1

(12)

Thereafter, since every step the estimate projects onto a geometrically simple space, it is efficient to compute this approximation. This approach is named sequential projection. Remark 3: The first thing should be remarked is that the sequential projection does not guarantee the constraint satisfaction of *1 or *2 . However, this does not matter much, owing to two reasons. First, the constraints *1 and *2 are formulated semiheuristically; strongly sticking to them may rule out the true parameter set. Second, constraint satisfaction is not the primary goal either; the main concern here is to find an efficient method for online implementation. Remark 4: Different orders of Θ *1, Θ *2 and Θ *3 may yield different projection results. Thus, the order may affect the final estimation performance. A plausible rule of thumb is that of placing the concrete constraint outside and putting the heuristic one inside. It is based on the philosophy that the outside

Here z is the center of the norm ball. Additionally, za, at present, means the first na elements of θ̂Rk (t) or Proj[θk̂ − 1(t )]. Later, we will introduce another version of z. Up to now, the constraints based priori plant knowledge have been formulated; their corresponding efficient solving approaches have all been presented.

5. ROBUST ANALYSIS It is well-known that a batch process may exhibit slow interbatch dynamic drift, which arises from facilities aging or accumulation of residual reactants.27 Therefore, it is important to investigate the robustness of the proposed identification algorithm against interbatch dynamics drift. Consider that A(t , q−1)yk (t ) = B(t , q−1)uk(t ) + dk(t ) 8823

(16)


Article


Figure 2. Identification results comparison among Landau’s method, conventional method, and the proposed method. K

where dk(t) represents the effect of model-plant mismatch and disturbances. This setup is quite general, and widely used in identification literature, for instance Landau et al.42 The robust analysis studies, given dk(t) and the reference rk(t) are norm bounded, are used to determine the boundedness of ŷk(t) and posterior prediction error yk(t) − φTk θ̂Rk (t). Obviously, dk(t) has absorbed the noise term ek(t). Theorem 1: Consider the plant represented in eq 16 and the parameter estimates generated by algorithm in eq 8 together with estimate refinement procedures eqs 13−15. Assume that A1. the order of the plant na, nb are exactly known, A2. the controller (3) stabilizes the plant (16), A3. the reference rk(t) and unmodeled plant dynamics dk(t) are both bounded in the sense that there exist Cr1, Cr2, Cr3, Cr4 ∈ (0, + ∞) such that

≤ Ce1K + Ce2

i=1

for any K and uniformly for all t. Proof: Since the proof is quite technical, please refer to Appendix A for details.

6. NUMERICAL EXPERIMENT 6.1. Adaptive Filter. The imposition of the constraints *1, *2, *3 will reduce the magnitude of fluctuation on the estimates; in practice, in order to further improve the smoothness, it is suggested to employ the following adaptive filter:35 ⎧ ̂R ̂f ⎪ γ1θk (t ) + (1 − γ1)θk (t − 1), R f ⎪ f k θk̂ (t ) = ⎨ ∥θk̂ (t ) − θk̂ (t − 1)∥ ≥ κ1 + κ2κ3 ⎪ R ⎪ γ θ ̂ (t ) + (1 − γ )θ ̂ f (t − 1), otherwise ⎩2 k 2 k

K

∑ ri2(t ) ≤ Cr1K + Cr 2 i=1

and

(17)

K

∑

R

∑ [yi (t ) − φiT(t )θk̂ (t )]2

di2(t )

It is a low-pass filter. γ1 and γ2 are positive real numbers between 0 and 1. κ1, κ2, and κ3 are all positive real numbers; κ3 is between 0 and 1. It is suggested that the γ2 should have a stronger smoothing effect to neutralize the small-magnitude fluctuations; that means γ2 should be smaller than γ1. κ1 should be closer to the maximum dynamics change within an iteration. κ2 and κ3 are two tuning knobs for the filter; they have a decaying effect on the filter switch threshold. The goal of introducing them is to allow the algorithm to aggressively pursue minimizing the prediction error at first and gradually switch to pursue smoothing estimates later. According to experiences, γ1 is usually selected between 0.7 and 0.9, while γ2 is chosen between 0.2 and 0.5. κ2 is usually selected as 1.2κ1; κ3 is between 0.7 and 0.9. The filter is organized in a convex combination form to avoid massive violations of the constraints *3; the trust region calculation part can handle the violation by rejecting the filtered version but keep using θ̂fk−1(t). Then,

≤ Cd1K + Cd 2

i=1

for any K and uniformly for all t, A4. all the priori constraints *1, *2, *3 cover the real parameter set θ0, A5. there exists a positive real number μ > 1 such that

μ S − H (t ) 2 is strictly positive real and stable uniformly for all t. Then, there exists Cp1, Cp2, Ce1, Ce2 ∈ (0, + ∞) such that K

∑ yi2̂ (t ) ≤ Cp1K + Cp2 i=1

and 8824


Article

Industrial & Engineering Chemistry Research another version of z in eq 15a will be θ̂fk(t). The whole proposed algorithm has been summarized in Appendix B. 6.2. Example 1. In this part, we consider the adapted example from Example 1 in Laudau and Karimi’s work.42 The system is set up as follows: ⎧1 ⎪ ⎪ ⎪1 −1 A(z ) = ⎨ ⎪ ⎪ ⎪1 ⎩

− 1.5z −1 + 0.7z −2 ,

(Filtered FCSC). Figure 3 suggests that incorporating priori plant knowledge is important in batch process identification. The maximum spectral radius for the kth iteration is defined as supt∈[0,L]ρmax [Â k(t)]. Figure 4 shows the evolution of the

t ∈ [0, 150)

⎡ ⎤ 0.32 + ⎢ − 1.5 + × (t − 150)⎥z −1 + 0.7z −2 , ⎣ ⎦ 150 t ∈ [150, 300) − 1.18z−1 + 0.7z −2 ,

t ∈ [300, 600]

−1 ⎧ t ∈ [0, 150) + 0.5z −2 , ⎪z B(z −1) = ⎨ ⎪ −1 −2 ⎩ 0.9z + 0.5z , t ∈ [150, 600]

Figure 4. Evolution of maximum spectral radius within an iteration in the proposed method.

The R-S-T controller is chosen as R(z −1) = 0.8659 − 1.2763z −1 + 0.5204z −2 S(z −1) = 1 − 0.6283z −1 − 0.3717z −2

maximum spectral radius. It suggests that in the proposed method, the stability constraint has never been violated, while it is difficult for the conventional method to fulfill the stability constraint, because the priori stability level η = 0.9 has always been satisfied by the proposed method, that is, supt∈[0,L]ρmax [Â k(t)] ≤ 0.9 for any k. To see that a smooth identified model can truly improve control performance, Figure 5 shows a tracking performance comparison of minimum variance control, one of the advanced control strategies.10 The minimum variance controller uses the identified models of the proposed method and the conventional method in the 600th iteration to regulate the plant A(z−1), B(z−1). Figure 6 shows the PTE logarithmic plot of the conventional and the proposed methods under the noisy scenario (the noise level is 0.6). In this case, the PTE of the conventional method starts to increase after 100 iterations, while the proposed method is still able to reduce PTE. 6.3. Example 2. Example 2 is from a real industrial application−injection molding. Injection molding, as a typical batch process, plays a pivotal role in polymer processing. It transforms plastic granules into various polymer products with high versatility and productivity. Figure 7 shows a simplified schematics of an injection molding machine, which comprises four parts−injection unit, clamping unit, hydraulic unit, and control unit. A typical operation cycle of injection molding consists of three phasesfilling, packing, and holding and cooling. Within this paper, it attempts to identify a model between valve opening and injection velocity (filling phase). The data are collected from the injection molding machine Haitian MA-600 located in Kunda Mold Company, Shenzhen, China. The material used is PA-757 ABS (acrylonitrile butadiene styrene). The barrel temperature is regulated in a closed loop with four zones, whose set-points are 225 °C, 225 °C, 220 °C, and 205 °C (from nozzle). The machine is working in air-shot mode. Air-shot mode means that the mold is detached from the barrel; it is an important mode for studying the injection molding dynamics.46 Figure 8 shows the step response test in open loop. The step change is introduced in t = 100 ms with a magnitude of 10% valve opening. The step response is prefiltered with 12-order averaging filter. According to Figure 8, an approximate FOPTD model is obtained by the graphic method,37 which is

The system A(z−1), B(z−1) is regulated by the controller R(z−1), S(z−1) repeatedly over a finite period t ∈ [0,600]. λk, η, γ1, γ2, κ1, κ2, and κ3 are respectively chosen as 0.98, 0.9, 0.7, 0.3, 0.8, 1.25, and 0.8. The reference signal is selected as a step change from 0.5 to −0.5 at time 150. The unrepeatable part rnr,k(t) in eq4 is designed as the PRBS (pseudorandom binary sequence) with magnitude of 0.3. The white noise variance is 0.2. The priori time constant and static gain are 0.8 and 5, respectively. We use Landau’s method to initialize the algorithm (identify the first iteration/batch).42 Figure 2 shows the identification results of three different methodsLandau’s method (Figure 2a,d), conventional method in eq 8 (Figure 2b,e), and the proposed method (Figure 2c,f). In Figure 2, the blue solid line represents the estimates; the red dash line stands for the true parameter. It clearly shows that not changing data processing or estimate updating direction yields a slow tracking response. It also indicates that the proposed method can significantly reduce the fluctuations on the estimates. Figure 3 is

Figure 3. PTE comparison of 100 Monte Carlo simulations.

the logarithm plot of the PTE comparison. PTE is short for parameter tracking error.35 Its definition is PTE = ∑tL= 1∥θ̂k(t) − θ0(t)∥2. The blue line is the result of directly applying the conventional method (eq 8); the orange line represents the result with the stability constraints *3 (SC) applied; the magenta line is the result with both stability constraints *3 and first-order constraints *1, *3 (FCSC) applied; the red line is the result of the proposed method with filter and *1, *2, *3 applied 8825


Article


Figure 5. Minimum variance control results based on the identification of the proposed method and the conventional method.

Figure 8. Step response of injection velocity (IV) vs valve opening (VO).

Figure 6. PTE comparison of 100 Monte Carlo simulations for noisy data.

G (s ) =

1.16e−30s 70s + 1

with millisecond as time unit. Afterward, 80 closed-loop PRBS tests are conducted. The duration of a cycle (iteration) is 2 s; the sampling period is 10 ms. The injection velocity is controlled by PID in velocity form; that is, R = S = T = 1 − q−1. The control block diagram is shown in Figure 9. The loop of valve opening to injection velocity is closed, and valve opening is the manipulated variable. The reference is 60 mm/s on and before time t = 100 and 40 mm/s afterward. The magnitude of PRBS on reference is 6 mm/s. Since the step response test is introduced in t = 100 and in the beginning of each cycle, the injection molding system is not very steady. Also recalling that injection molding has strong time variations,39 it is reasonable

Figure 9. Diagram of closed-loop control block of Example 2.

to impose *1, *2 over the last 2/3 duration. η, γ1, γ2, κ1, κ2, κ3 are selected as 0.9, 0.7, 0.5, 0.8, 1.25, 0.8, respectively. According to the literature,39 the order is chosen as na = 1, nb = 2, nd = 3. Figure 10 shows the comparison of identification results. It suggests that the prediction ability of Landau’s method does not improve from iteration to iteration, but the

Figure 7. Simplified schematics of injection molding machine39 8826


Article

Industrial & Engineering Chemistry Research From the error dynamics (eq 7), one can have

−ϵ*k = H φkTθk̃ − ηk* * where η*k ≔ Sηk/H and ϵ*k ≔ ϵk + μφTk θ̃k/2 and H = S/H − * μ/2. Then, dividing H and multiplying ϵ*k on both side of the * above equation, we have −H −1(ϵ*k )2 = φkTθk̃ ϵ*k − H −1ηk*ϵ*k * *

(19)

Since H is positive real, so is H −1. According to Lemma D.1 * * in ref 42 or section (b) in ref 47 (pp. 560), there must exist constant γ > 0 and l > 0 such that Figure 10. Identification results comparison among the three methods on Example 2.

H −1(ϵ*k )2 ≥ −γ + l(ϵ*k )2 *

(20)

From eq 19 and eq 20, it follows that conventional eq 8 and the proposed methods can. The proposed method is slightly worse than the conventional; this is understandable, because the proposed is suboptimal with respect to the empirical cost (prediction error) with the presence of constraints. It also implies that the proposed method yields much smoother estimates than the conventional eq 8 quantified by the index SI. SI stands for smoothness index, which is defined as SI = ∑tL= 2∥θ̂k(t) − θ̂k(t − 1)∥2. Thus, the estimates given by the proposed method are much more reliable.

φkTθk̃ ϵ*k ≤ γ − l(ϵ*k )2 + H −1ηk*ϵ*k * (H −1ηk*)2 * ≤γ+ − (l − β)(ϵ*k )2 4β

(21)

The second inequality comes from Young’s inequality and β < l. Adding ∑ik= 1μφTi θ̃i/2 on both sides of eq 18, one can get that k

∑ θĩ i=1

7. CONCLUSIONS In this paper, an online closed-loop identification algorithm particular for batch process has been devised. The algorithm employs three kinds of priori plant knowledge, time constant, static gain, stability, to refine the yielded parameter estimates. The robustness of the approach has been mathematically studied. The approach can not only improve prediction ability from iteration to iteration, but also output smooth parameter estimates in the mean time. The effectiveness has been shown on a real industrial application−injection molding.

k

⎡μ 1 1⎤ 1 φi ϵ*k ≥ − V0 + ⎢ − ⎥ ∑ (φi Tθi)̃ 2 ≥ − V0 ⎦ ⎣ 2 2 2 i=1 2

T

(22)

because μ > 1. Combining eq 21 and eq 22, we have k

(l − β) ∑ (ϵ*i )2 ≤ γk + i=1

1 4β

k

∑ (H*−1ηi*)2 + i=1

1 V0 2

Since ηk is bounded and H is stable, ηk* is bounded as well. k Further from H −1 is stable; ∑i = 1 (H −1ηi*)2 is bounded too. * * Thus, there exists C1 and C2 such that ∑ik= 1 (ϵi*)2 ≤ C1k + C2. Since φTk θ̃k = H −1ϵk* + H −1ηk* and with the definition of ϵk*, * * the relation between ϵk and ϵ*k turns to be ϵk = (1 + μH −1/2)ϵ*k * − μH −1ηk*/2. Thus, ϵk is upper bounded by 2[(1 + * μH −1/2)ϵk*]2 + 2[ μH −1ηk*/2]2. The boundedness of ϵk can * * be concluded from the boundedness of ϵ*k and η*k and because H −1 is stable. Then, since yk is bounded, the boundedness of ŷk * follows. The boundedness of φk directly follows. From φTk θ̃k = −H −1ϵk* + H −1ηk* and the definition of ϵk*, it further concludes * * that θ̃k is bounded. Additionally, since θ0 ∈ Θ *1 ∩ Θ *2 ∩ Θ *3, according to the property of projection operator,48 we have

■

APPENDIX A Proof: The proof is based on the idea in ref 42. Within this section, to keep the notations neat, the argument t will be dropped without any ambiguity. For instance, θ̂k(t) will be reduced to θ̂k. From eq 8, it is known that Kk = Pkφk. Then, together with the definition of θ̃k, eq 8a becomes θk̃ = θk̃ − 1 + Pkφk ϵk

̃ Define the functional Vk = θ̃Tk P−1 k θk. Combining with eq 8, after some algebraic manipulations, the relation between Vk and Vk−1 can be established as T

Vk = Vk − 1 + 2θk̃ φk ϵk − φkTPk − 1φk ϵk2 + (φkTθk̃ )2

R ∥θk̂ − θ0∥ ≤ ∥Proj* [Proj* (θk̂ )] − θ0∥ 2

It follows that 1 [Vk − Vk − 1 − (φkTθk̃ )2 ] 2 Summing the equation from 1 to k, by telescoping, we have T

θk̃ φk ϵk ≥

k

∑ θĩ i=1

1 1 φi ϵi ≥ − V0 − 2 2

T

1

≤ ∥θk̂ − θ0∥ = ∥θk̃ ∥

k

∑ (φiTθi)̃ 2 i=1

1

≤ ∥Proj* (θk̂ ) − θ0∥

Hence θ̂Rk is bounded. The boundedness of yk − φTk θ̂Rk follows directly.

(18) 8827


Article


■

APPENDIX B

■

AUTHOR INFORMATION

(8) Zhang, R.; Xue, A.; Wang, S.; Ren, Z. An improved model predictive control approach based on extended non-minimal state space formulation. J. Process Control 2011, 21, 1183−1192. (9) Zhang, R.; Xue, A.; Wang, S. Modeling and nonlinear predictive functional control of liquid level in a coke fractionation tower. Chem. Eng. Sci. 2011, 66, 6002−6013. (10) Åström, K. J.; Wittenmark, B. Adaptive Control; Courier Corporation: 2013. (11) Ioannou, P. A.; Sun, J. Robust Adaptive Control; Courier Corporation: 2012. (12) Slotine, J.-J. E.; Li, W. On the adaptive control of robot manipulators. Int. J. Robotics Res. 1987, 6, 49−59. (13) Craig, J. J.; Hsu, P.; Sastry, S. S. Adaptive control of mechanical manipulators. International Journal of Robotics Research 1987, 6, 16−28. (14) Mosca, E. Optimal, predictive, and adaptive control. IEEE Control Syst. 1995, 15, 92. (15) Ljung, L. System Identification: Theory for the User; Prentice Hall: NJ, 1987. (16) Gustavsson, I.; Ljung, L.; Söderström, T. Identification of processes in closed loop−identifiability and accuracy aspects. Automatica 1977, 13, 59−75. (17) Söderström, T.; Stoica, P. System identification; Prentice-Hall, Inc.: 1988. (18) Mišković, L.; Karimi, A.; Bonvin, D.; Gevers, M. Closed-loop identification of multivariable systems: With or without excitation of all references? Automatica 2008, 44, 2048−2056. (19) Hjalmarsson, H. From experiment design to closed-loop control. Automatica 2005, 41, 393−438. (20) Gevers, M. Identification for Control: From the Early Achievements to the Revival of Experiment Design. European Journal of Control 2005, 11, 335−352. (21) Forssell, U.; Ljung, L. Closed-loop identification revisited. Automatica 1999, 35, 1215−1241. (22) Huang, B.; Shah, S. L. Closed-loop identification: a two step approach. J. Process Control 1997, 7, 425−438. (23) Shardt, Y. A.; Huang, B. Closed-loop identification with routine operating data: Effect of time delay and sampling time. J. Process Control 2011, 21, 997−1010. (24) Gilson, M.; Garnier, H.; Young, P. C.; Van den Hof, P. M. Optimal instrumental variable method for closed-loop identification. Control Theory & Applications, IET 2011, 5, 1147−1154. (25) Gilson, M.; Garnier, H.; Young, P. C.; Van den Hof, P. Identification of Continuous-Time Models from Sampled Data; Springer: 2008; pp 133−160. (26) Huang, B.; Ding, S. X.; Qin, S. J. Closed-loop subspace identification: an orthogonal projection approach. J. Process Control 2005, 15, 53−66. (27) Bonvin, D. Optimal operation of batch reactors−a personal view. J. Process Control 1998, 8, 355−368. (28) Bristow, D.; Tharayil, M.; Alleyne, A. G. A survey of iterative learning control. Control Systems, IEEE 2006, 26, 96−114. (29) Ahn, H.-S.; Chen, Y.; Moore, K. L. Iterative learning control: brief survey and categorization. IEEE Transactions on Systems Man and Cybernetics part C Applications and Reviews 2007, 37, 1099. (30) Zhang, R.; Xue, A.; Wang, J.; Wang, S.; Ren, Z. Neural network based iterative learning predictive control design for mechatronic systems with isolated nonlinearity. J. Process Control 2009, 19, 68−74. (31) Ma, D. L.; Braatz, R. D. Robust identification and control of batch processes. Comput. Chem. Eng. 2003, 27, 1175−1184. (32) Tayebi, A. Adaptive iterative learning control for robot manipulators. Automatica 2004, 40, 1195−1203. (33) Chi, R.; Hou, Z.; Xu, J. Adaptive ILC for a class of discrete-time systems with iteration-varying trajectory and random initial condition. Automatica 2008, 44, 2207−2213. (34) Cao, Z.; Yang, Y.; Lu, J.; Gao, F. Constrained two dimensional recursive least squares model identification for batch processes. J. Process Control 2014, 24, 871−879.

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This project is supported in part by the National Natural Science Foundation of China Project Nos. 61227005 and 61503181, Guangdong Scientific and Technological Project No. 2013B090800004, Guangdong Natural Science Foundation Project No. 2015A030310266, Guangdong Innovative and Entrepreneurial Research Team Program No. 2013G076, Guangzhou Science and Technology Bureau Project No. 201508030040.

■

REFERENCES

(1) Jeng, J.-C.; Lin, Y.-Y. Closed-loop identification of dynamic models for multivariable systems with applications to monitoring and redesign of controllers. Ind. Eng. Chem. Res. 2011, 50, 1460−1472. (2) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. Constrained model predictive control: Stability and optimality. Automatica 2000, 36, 789−814. (3) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Control engineering practice 2003, 11, 733−764. (4) Camacho, E. F.; Alba, C. B. Model Predictive Control; Springer Science & Business Media: 2013. (5) Morari, M.; Lee, J. H. Model predictive control: past, present and future. Comput. Chem. Eng. 1999, 23, 667−682. (6) Chen, H.; Allgöwer, F. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica 1998, 34, 1205−1217. (7) Zhang, R.; Xue, A.; Wang, S.; Zhang, J. An improved state-space model structure and a corresponding predictive functional control design with improved control performance. Int. J. Control 2012, 85, 1146−1161. 8828


Article

Industrial & Engineering Chemistry Research (35) Cao, Z.; Zhang, R.; Lu, J.; Gao, F. Online identification for batch processes in closed loop incorporating priori controller knowledge. Comput. Chem. Eng. 2016, 90, 222−233. (36) Cao, Z.; Zhang, R.; Lu, J.; Gao, F. Two-time dimensional recursive system identification incorporating priori pole and zero knowledge. J. Process Control 2016, 39, 100−110. (37) Seborg, D.; Edgar, T. F.; Mellichamp, D. Process Dynamics & Control; John Wiley & Sons: 2006. (38) Gao, F.; Yang, Y.; Shao, C. Robust iterative learning control with applications to injection molding process. Chem. Eng. Sci. 2001, 56, 7025−7034. (39) Yang, Y.; Gao, F. Adaptive control of the filling velocity of thermoplastics injection molding. Control Engineering Practice 2000, 8, 1285−1296. (40) Yang, Y.; Gao, F. Cycle-to-cycle and within-cycle adaptive control of nozzle pressure during packing-holding for thermoplastic injection molding. Polym. Eng. Sci. 1999, 39, 2042−2063. (41) Li, M.; Yang, Y.; Gao, F.; Wang, F. Fuzzy multi-model based adaptive predictive control and its application to thermoplastic injection molding. Can. J. Chem. Eng. 2001, 79, 263−272. (42) Landau, I. D.; Karimi, A. Recursive Algorithms for Identification in Closed Loop: A Unified Approach and Evaluation. Automatica 1997, 833, 1499−1523. (43) Mosca, E. Optimal, Predictive, and Adaptive Control; PrenticeHall, Inc.: 1995. (44) Skogestad, S. Simple analytic rules for model reduction and PID controller tuning. J. Process Control 2003, 13, 291−309. (45) Bertsekas, D. P. Nonlinear Programming; Athena Scientific: 1999. (46) Gao, F. The control of cavity pressure throughout the injection molding cycle. Ph.D. Thesis, McGill University, 1993. (47) Caines, P. Linear Stochatic Systems; John Wiley & Sons: 1988. (48) Nedić, A.; Ozdaglar, A.; Parrilo, P. Constrained consensus and optimization in multi-agent networks. IEEE Trans. Autom. Control 2010, 55, 922−938.

8829


Priori Knowledge-Based Online Batch-to-Batch Identification in a

Recommend Documents