Online identification for batch processes in closed loop incorporating priori controller knowledge

Computers and Chemical Engineering 90 (2016) 222–233

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Online identification for batch processes in closed loop incorporating priori controller knowledge Zhixing Cao a , Ridong Zhang a , Jingyi Lu a , Furong Gao a,b,∗ a b

Department of Chemical & Biomolecular Engineering, Hong Kong University of Science & Technology, Clear Water Bay, Kowloon, Hong Kong Fok Ying Tung Research Institute, Hong Kong University of Science & Technology, Guangzhou, China

a r t i c l e

a b s t r a c t

i n f o

Article history: Received 14 December 2015 Received in revised form 4 April 2016 Accepted 18 April 2016 Available online 20 April 2016 Keywords: Batch processes Closed-loop system identification Priori knowledge Process modeling Trust region method Two-time dimensional

It is of great importance to develop an online modeling method for chemical processes operated in closed loop for better understanding, monitoring the process or other purposes without endangering the system. This paper intends to devise an online system identification method, particularly for the batch process, by fully exploiting its intrinsic repetitiveness. It properly uses the information from the time direction and the batch direction, thus leading to a gradual performance enhancement. In addition, the identification method formulates the priori controller knowledge such as closed-loop stability as optimization constraints to refine the parameter estimates. A trust region method is employed to overcome the significant computation burden of directly handling these constraints such as solving Lyapunov inequalities. An adaptive filter is introduced to further smooth the parameter estimates. Finally, the effectiveness of the approach is verified by three numerical examples including a two-tank system. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Dynamic mathematical models play a pivotal role in chemical engineering, especially in process control, as they provide engineers a precise evolutional picture of the process of interests, which renders the future performance improvement possible, such as redesigning controllers or set points. Generally speaking, dynamic mathematical models are mechanism-based, data-based or a mixture of the both. Mechanism-based models link the process input and output by invoking fundamental physical or chemical principles, e.g., mass conservation, with a bunch of ordinary differential equations (ODEs) or partial differential equations (PDEs). These equations are generally difficult to solve, or it lacks efficient methods to solve them; that makes some advanced control or monitoring strategy intractable, like adaptive control (Åström and Wittenmark, 2013; Ioannou and Sun, 2012; Slotine and Li, 1987; Craig et al., 1987; Mosca, 1995), model predictive control (Mayne et al., 2000; Qin and Badgwell, 2003; Camacho and Alba, 2013; Morari and Lee, 1999; Chen and Allgöwer, 1998), etc. Unlike mechanism-based models, data-based models take a different approach to circumvent these problems, by only focusing on the relationship between process input and output instead of emphasizing the real mechanism between them. In the discipline of process control, the methods

∗ Corresponding author. E-mail address: [email protected] (F. Gao). http://dx.doi.org/10.1016/j.compchemeng.2016.04.025 0098-1354/© 2016 Elsevier Ltd. All rights reserved.

to build dynamic data-based model are called system identification. There are two categories of system identification methods according to the operation condition of the plant under investigation – open loop and closed loop. In the past decades, lots of interesting results have contributed to the open-loop system identification, no matter in theory (Ljung, 1987; Söderström and Stoica, 1988) or in application (Zhu, 2001; Huang and Shah, 2012; Ikonen and Najim, 2001). However, in reality, there are many situations where openloop identification is difficult to implement or even not a feasible option. The simplest example to illustrate this is a plant exhibits integral behaviour, or has open-loop pole(s) in the right half complex plane (in the continuous-time sense). In this particular case, the plant’s output always has an unbounded trend for almost any input signal; that compromises the identifiability of the open-loop identification and necessitates the closed-loop identification without jeopardizing the plant’s normal operation (Landau and Karimi, 1997). The closed-loop identification also provides room for controller tuning or controller redesign to further improve control performance. Motivated by these reasons, closed-loop identification has aroused lots of attention in control community. From the aspect of frequency analysis, Gustavsson, Ljung and Söderström deeply discussed the fundamental problem of closed-loop identification – identifiability (Gustavsson et al., 1977; Ljung, 1987; Söderström and Stoica, 1988). Zang, Hjalmarsson, Van den hof, and their colleagues studied the connection between identification in closed loop and controller redesign in the frequency domain

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

(Zang et al., 1995; Hjalmarsson et al., 1996; Van Den Hof and Schrama, 1995). Forssell and Ljung categorized and analyzed prediction error method in a unified framework (Forssell and Ljung, 1999). Huang and Shah developed a two-stage method and implemented it on a pilot-scale process (Huang and Shah, 1997). On the other hand, from the perspective of time domain, Landau and Karimi developed various recursive identification methods by duplicating the system and comparing the output error (Landau and Karimi, 1997, 1999). Huang, Qin, and their coworkers developed closed-loop subspace identification method based on whitening the error (Huang et al., 2005; Qin and Ljung, 2003). In this paper, our primary contribution is to devise a closedloop identification algorithm, particularly for batch process. It is known that batch process, different from continuous process, has its unique dynamic properties such as time variation (Yang and Gao, 2000). Thus, directly applying the aforementioned methods to batch process may not yield satisfactory identification results. For instance, suppose that a plant y(t) = G(t, q−1 )r(t) is operated perfectly under a certain controller with r(t) as its reference signal. If using prediction error method to identify the plant, the prediction error is comprised of two components: one is from the mismatch of parameter estimates, the other is from the dynamics variation. This can be seen from (t) = y(t) − yˆ (t|t − 1) = [G(t, q−1 ) − ˆ G(t − 1, q−1 )]r(t) + [G(t − 1, q−1 ) − G(t|t − 1, q−1 )]r(t). A plausible approach to resolve the problem is to exploit the repetitiveness of batch process, which borrows the idea from iterative learning control (ILC), using the information from two directions – time direction and batch direction. To be specific, the prediction error will be k (t) = yk (t) − yˆ k (t, k|k − 1) = [Gk (t, q−1 ) − ˆ k|k−1 (t, q−1 )]r(t); the first term Gk−1 (t, q−1 )]r(t) + [Gk−1 (t, q−1 ) − G −1 will disappear since Gk (t, q ) = Gk−1 (t, q−1 ).1 Similar literatures are: Ma and Braatz studied offline identification methods for batch process (Ma and Braatz, 2003); Tayebi (2004) and Chi et al. (2008) used the ILC idea but without considering transient identification performance; our previous paper (Cao et al., 2014) has only discussed the open-loop case. The second contribution of this paper is that we use the priori closed-loop knowledge to refine the parameter estimates. As pointed out in our previous work (Cao et al., 2014), the identification results directly generated by two-time dimensional identification algorithm entails severe variation on parameter estimates; to some extent, it contradicts the fact that most chemical processes have slow dynamic variations. The technique – soft constraint – has been employed to tackle this problem in our previous work (Cao et al., 2014); whereas in this paper, we intend to impose “hard” constraints instead. These constraints are formulated from the priori closed-loop knowledge. A straightforward example is that the closed-loop poles are all within the unit disk, which follows from the closed-loop stability. Unfortunately, these constraints usually appear as Lyapunov inequalities, a non-convex form or computationally unfriendly. Thus, a trust region method is taken advantage of to circumvent this challenge. To ensure the recursive feasibility, the “size” of the trust region is recursively estimated by the techniques developed in robust control. Interestingly, the “size” estimation process and the identification process can run in parallel. In addition, an adaptive low-pass filter is introduced to further neutralize the variation; in the mean time, the filter does not compromise the established recursive feasibility. The paper is organized as follows: Section 2 derives the unconstrained closed-loop identification methods based on minimum prediction error methods (MPE); Section 3 refines the parameter estimates by imposing constraints formulated from the priori

1

k stands for batch index, which can also be called iteration index or period index.

223

closed-loop knowledge; Section 4 analyzes the projection property and recursive feasibility; Section 5 verifies the proposed method by presenting three numerical examples; Section 6 draws a conclusion. Notations: t and k stand for time and batch index respectively. √ q−1 represents a unit time backward shift operator. M =  T M is the Euclidean norm. The hat (ˆ•) means estimation or something associated with estimates. ⊗ is the Kronecker product. vec is the vectorization operator. T signifies transposition. max (•), (•) are the maximum eigenvalue and spectral radius of (•) respectively. 2. Unconstrained closed-loop identification 2.1. Closed-loop output prediction Within the paper, we assume that a batch process can be delineated by the following discrete-time autogressive exogenous (ARX) model, A(t, q−1 )yk (t) = B(t, q−1 )uk (t) + ek (t),

(1)

where yk (t) and uk (t) are, respectively, the plant’s output and input at time t of the kth iteration. ek (t) is a two dimensional white noise, with E[ek (t)em (n)] =  2 ık,m ıt,n . ık,m , ıt,n are Kronecker deltas, and ık,m = 1 if and only if k = m. A(t, q−1 ) and B(t, q−1 ) are scalar polynomials with time varying coefficients as shown in the following equations, A(t, q−1 ) = 1 + a1 (t)q−1 + a2 (t)q−2 + · · · + ana (t)q−na , −1

B(t, q

−1

) = b1 (t)q

−2

+ b2 (t)q

−nb

+ · · · + bnb (t)q

.

(2a) (2b)

Here a1 (t), a2 (t), . . ., ana (t) and b1 (t), b2 (t), . . ., bnb (t) are the time varying coefficients to be identified. na, nb are the polynomial orders. In the paper, it is assumed that these orders are exactly known. It is necessary to put two remarks on the plant model in (1). First, in reality, most batch processes exhibit certain nonlinearity. However, in many cases, the plants are required to track a given set point trajectory; the objective is usually achieved by a feedback controller. This fact provides process control engineers sufficient rationales to approximate the real dynamics in the proximity of the trajectory with a linear time varying (LTV) model, i.e., (1). Second, from (2b), it is assumed that the input delay order (nd) is equal to 1. This assumption does not impose any restriction. Because it is trivial to extend to nd = d case by treating first d − 1 coefficients in (2b) as zero. The ARX model has been assumed, because it is simple and enough to capture batch processes’ dynamics in practice, for example, the work of Yang and Gao on injection molding (Yang and Gao, 2000). Furthermore, assume that the plant in (1) is operated in closed loop to track a given reference by an R-S-T controller (a two-degreeof-freedom controller). This type of controller has been reported widely applied in controller analysis, i.e., Landau and Karimi (1997, 1999). The block diagram of the closed loop is shown in Fig. 1. The associated control law is given as uk (t) = −

R(q−1 ) T (q−1 ) yk (t) + rk (t). S(q−1 ) S(q−1 )

(3)

Here rk (t) is the reference signal; it can be decomposed into two components like rk (t) = rr (t) + rnr,k (t). rr (t) is the repetitive component that the system is required to follow; rnr,k (t) is the non-repetitive part with relatively small magnitude to provide the system with sufficient excitation. nr, ns and nt correspond to the orders of scalar polynomials R(q−1 ), S(q−1 ) and T(q−1 ). Without loss of generality, the polynomial S(q−1 ) is given in the canonical form, i.e., S(q−1 ) = 1 + S* (q−1 ), where S* (q−1 ) = s1 q−1 + · · · + sns q−ns .

224

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

According to (5), we have 1 + H(t) = A(t)S + B(t)R. Then, it follows from the last line in (7) that

k [t|t − 1; (t)] = S −1 [1 + H(t)]yk (t) − S −1 B(t)Tr k (t) = S −1 [A(t)S + B(t)R]yk (t) − S −1 B(t)Tr k (t) = A(t)yk (t) + B(t)S −1 [Ryk (t) − Tr k (t)], (from 1 + H(t) = A(t)S + B(t)R) or it can be reformulated in linear regression form

k [t|t − 1; (t)] = yk (t) −

T (t)(t), k

(8)

with Fig. 1. Block diagram of the closed loop.

k (t)

f

f

T

= [−yk (t − 1), . . ., −yk (t − na), −uk (t − 1), . . ., −uk (t − nb)] .

uf signifies the filtered input S−1 [Ryk (t) − Trk (t)].

If S(q−1 ) = s0 + S* (q−1 ), letting R (q−1 ) = R(q−1 )/s0 , T (q−1 ) = T(q−1 )/s0  and S  (q−1 ) = 1 + S  (q−1 ), then the yielded R -S -T controller can exactly generate the same control input as the R-S-T controller. For simple notations, denote A(t, q−1 ) = 1 + A* (t, q−1 ). In addition, in the remainder of the paper, the argument q−1 will be dropped to further simplify the notations. Then, the closed-loop system can be derived from (1) and (3) as

Remark 2.1. It is noted from Fig. 1 that −uf is equal to the plant input, i.e., uk (t). Hence, the equivalent form of (8) is k [t|t − 1; (t)] = A(t)yk (t) − B(t)uk (t). Comparing with (1), the prediction error estimator (8) basically attempts to whiten the prediction error. Additionally, uf does not contain any parameter to be identified so that it can be exactly computed.

yk (t) = −H(t)yk (t) + B(t)Tr k (t) + S ∗ ek (t) + ek (t),

2.2. Minimum prediction error method

(4)

with H(t) = A∗ (t) + S ∗ + A∗ (t)S ∗ + B(t)R = 1 + A∗ (t) + S ∗ [1 + A∗ (t)] + B(t)R − 1 = A(t) + S ∗ (t)A(t) + B(t)R − 1 = A(t)S + B(t)R − 1.

(5)

The MPE method has played an important role in the field of system identification (Ljung, 1987; Söderström and Stoica, 1988). Its basic mechanism is to adjust the parameter estimates by minimizing a proper cost function of the prediction error. In classic system identification literatures (Ljung, 1987; Söderström and Stoica, 1988; Mosca, 1995), the cost function is selected as 1  t−i 2  k [i|i − 1; (i)]. 2 t

From (4), it clearly shows that the closed-loop system preserves both the linearity and time variation of the original system (1). Moreover, it also suggests that the innovation ek (t) can be estimated from the plant’s previous input and output, c.f., {yk (1), yk (2), . . ., yk (t)} and {uk (0), uk (1), . . ., uk (t − 1)} via ek (t) = S −1 [1 + H(t)]yk (t) − S −1 B(t)Tr k (t).

(6)

Thus, it leads to a reasonable one-step output prediction yˆ k (t|t − 1; ) based on the input and output up to time t within the current batch (kth batch) and the true parameters (t) = [a1 (t), . . ., ana (t), b1 (t), . . ., bnb (t)]T . ∗

Jk (t) =

Here  ∈ (0, 1] is the forgetting factor. This index is the summation of all the prediction error up to now within the kth batch (normally, it is the current batch). In this case, when the optimization problem solved with respect to Jk (t) above, the time-varying dynamics incurs a bias on the parameter estimates as argued before. Thus, solving the optimization problem is to find the best tradeoff between the mismatch on (i) and the noise impact. This motivates us to revise the cost function based on the repetitiveness of the process,

⎡ ⎤ k k   1 ⎣ Jk (t) = j ⎦ 2i [t|t − 1; (t)], 2

(4)

i=1

yˆ k [t|t − 1; (t)] = −H(t)yk (t) + B(t)Tr k (t) + S ek (t) =yk (t) (6)

− ek (t) =yk (t) − S −1 [1 + H(t)]yk (t) + S −1 B(t)Tr k (t). (7)

(10)

j=i+1

k

with the convention j=k+1 j = 1. Fig. 2 illustrates how the weight i enters Jk (t). It should be noted that the sequence {i } starts at 2,

It should be remarked that yˆ k [t|t − 1; (t)] is the minimum meansquare error (MMSE) prediction of yk (t) provided that  represents the real parameters of the system (1). This can be easily seen from E[yk (t) − yˆ k (t)2 ] = E{ˆyk [t|t − 1; (t)] − yˆ k (t) + ek (t)2 } = E{ˆyk [t|t − 1; (t)] − yˆ k (t)2 } + E[ek (t)2 ], where yˆ k (t) is the arbitrary output prediction based on the previous I/O data (Mosca, 1995). It is clear that the minimum of E[yk (t) − yˆ k (t)2 ] is only attained at yˆ k (t) = yˆ k [t|t − 1; (t)]. The prediction error will be denoted as

k [t|t − 1; (t)] = yk (t) − yˆ k [t|t − 1; (t)].

(9)

i=1

Fig. 2. Illustration of weights adding.

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

i.e. 2 ; although it seems weird, it is for easy memorizing. i stands for the “discount” selected in the ith batch for all the previous data; new data enters with the default weight 1. The underlying idea of the cost function is that the summation direction alters from the time direction to the batch direction. The motivation behind this is that at a given time t, all the outputs yi (t) (or equivalently i [t|t − 1; (t)]) can be considered as sampled from an oracle (data generator) parameterized with (t). If the summation includes terms of a different time, i.e. t , then minimizing Jk (t) means to find the best ˆ fitting parameter (t) to describe the data generated by the two oracles – (t), (t ); there has to be a tradeoff between fitting the data of (t) and that of (t ), thus compromising the identification algorithm’s consistency. In simple words, with the index in (10), an LTV identification problem is turned into an “LTI” problem. In (10), the forgetting factor j is given in a more general form; j is allowed to select different values between 0 and 1. The motivation to introduce j is to handle the slow dynamic drifts between batches and expedite the convergence rate. Since the algorithm is proposed for online identification, there is no requirement on minimum number of batches; but the online function can be terminated if the result satisfies the identification requirement. As seen in (8), (10) can be rewritten as 1 2 k

Jk (t) =

i=1



k 

T T (t)(t)] i

{[

j

− 2yi (t)

T (t)(t) i

+ yi2 (t)}



d2 Jk (t) d 2 (t)

d2 Jk−1 (t) d2 Jk (t) = +  k d 2 (t) d 2 (t)

 ˆ k−1 (t)

−1

dJ k (t) d(t)

 ˆ k−1 (t)

,

(11)

= k

T (t), k

k (t)

dJ k−1 (t) d(t)

(12a)

−

k (t)k [t|t

ˆ k−1 (t)

Introduce the notation Pk (t)  [d2 Jk (t)/d 2 (t)] from (12a) that Pk−1 (t)

=

−1 k (t)]

−1 k Pk−1 (t) +

In the previous section, the identification algorithm is carried out based on solving an unconstrained optimization problem. In this section, constraints will be included into the optimization problem to further enhance the identification performance. The constraints confine the estimates to the region, where the closedloop characteristic polynomial associated with the estimates is stable or back off certain distance to the unit circle. Prior to that, a simple example will be presented to give readers a straightforward impression of the necessity to introduce these constraints. 3.1. Example 1 Consider the following ARX system and control law,

System 1:

yk (t) − 1.1yk (t − 1) = uk (t − 1) + ek (t)

t ∈ [1, 150]

yk (t) − 1.2yk (t − 1) = uk (t − 1) + ek (t)

t ∈ (150, 400],

k (t)

T (t). k

− 1; ˆ k−1 (t)].

Controller 1: uk (t) = −0.3yk (t) + rk (t).

(12b) −1

. Then, it follows

d(t)

ˆ k−1 (t)

=−

k (t)k [t|t

PTE(k) =

L 

ˆ k (t) − (t)2 ,

t=1

(13)

− 1; ˆ k−1 (t)].

(14)

By injecting the definition of Pk (t) and (14) into (11), the updating equation on the parameter estimates ˆ k (t) is ˆ k (t) = ˆ k−1 (t) + Pk (t)

k (t)k [t|t

− 1; ˆ k−1 (t)].

(15)

Notice that it involves much computation in (13) due to the existence of matrix inversion. Taking advantage of matrix inversion lemma (Zhu, 2001), (13) can be transformed into

(17b)

Here the repeatable part of rk (t) is 0.5 before time t = 150 and 1 afterwards. The unrepeatable part is the pseudorandom binary sequence (PRBS) with the magnitude equal to 0.1. It is clear that System 1 is not identifiable due to its instability; the system will diverge within a finite time. Aided by Controller 1, the process identification becomes implementable. Moreover, it is easy to see the closed-loop poles are within the circle of radius 0.95 on the z-plane. Figs. 3 and 4 show the identification results corresponding to the algorithm derived in Section 2. The identification of the first batch is carried out based on solving the optimization problem with respect to the index defined in (9). Qualitatively speaking, it is shown in Fig. 3 that the identification based on (9) cannot respond fast enough to track the dynamic change at time t = 150, whereas the estimates, in the 50th batch, is able to follow the time-varying dynamics. On the other hand, the estimates corresponding to the 50th batch are “jumping” around the true value; the first batch’s is quite smooth. For the sake of quantitative description, define the following two indices

Furthermore, as mentioned before, ˆ k−1 (t) is the minimizer with respect to Jk−1 (t). Thus, (12b) is simplified to

dJ (t)  k

T (t)Pk−1 (t)} k

(17a)

where ˆ k−1 (t) is the parameter estimates based on I/O data up to the (k − 1)th batch, but also a minimizer to the cost function of the batch, i.e., Jk−1 (t). Therefore, for the sake of online implementation, it is necessary to have the following recursive relationships of the cost function Jk (t) and its derivatives,

dJ (t)  k

T (t)Pk−1 (t) k

3. Constrained identification with priori knowledge

j=i+1

ˆ k (t) = ˆ k−1 (t) −

+

To this end, we have completed the unconstrained algorithm derivation; (8), (15) and (16) constitute the identification algorithm.

 T (t)(t) i

k (t)[k

(16)



The right-hand side of the equation above is a quadratic form of (t) so that it is not difficult to calculate its first and second-order derivatives. The next step is to derive a recursive minimizing algorithm suitable for online implement with relatively small computation. Using Newton–Raphson method to construct the next-step estimates (Söderström and Stoica, 1988), it yields that

d(t)

Pk (t) = 1/k {Pk−1 (t) − Pk−1 (t)

225

Fig. 3. The identification results of the parameter a1 (t).

226

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

assumed to be obtained as priori knowledge. Thereafter, for any k and i, at every time t the identification problem can be reformulated as minˆ

Jk (t)

s.t.

ˆ clo (q−1 , k, t)] < . ri [P

k (t)

(18)

Note that there are lots of difficulties to handle the constraint (18), since ri is not a continuous function not to mention the convexity. Fortunately, the constraint in nature is a stability condition; that renders us to use the powerful tool – Lyapunov equation theorem. Before invoking the theorem, the characteristic equation ˆ clo (q−1 , k, t) should be rewritten into matrix form like P

⎡

Fig. 4. The spectral radius of the closed-loop characteristic polynomial associated with the parameter estimates.

SI(k) =

L 

−hˆ 1,k (t) −hˆ 2,k (t) · · · −hˆ nh−1,k (t) −hˆ nh,k (t)

ˆ clo (q P

−1

ˆ k (t)  , k, t)∼H

ˆ k (t) − ˆ k (t − 1)2 .

⎢ ⎢ ⎢ ⎢ ⎣

t=2

The first index is called parameter tracking error (PTE), which quantifies an identification algorithm’s process dynamics tracking performance. The second one is named smoothing index (SI); it describes the smoothness of the identification results. The comparison results are summarized in the following table. Surprisingly, Table 1 shows that the PTE of the fiftieth batch almost triples that of the first batch, which counters the direct intuitive that the 50th batch is better. The index SI explains the reason attributing to the conclusion; it is the severe “jumping” that accounts for the observation. In fact, this result is quite understandable, because the algorithm yielded index (9) intrinsically has the capability to smooth its estimates (Mosca, 1995), i.e., limt→∞ ˆ k (t) − ˆ k (t − 1) = 0. Unfortunately, the algorithm based on (10) does not have such an ability. Therefore, it is necessary to suppress the “jumping”; and Fig. 4 inspires us with a possible solution. Fig. 3 interprets the result from the frequency domain; it suggests imposing a closed-loop constraint on the estimates. To be specified, supposing that the closed-loop poles are within the circle of 0.95 on the z-plane, it is sensible to require that the corresponding poles of the estimates should be confined to the circle. These circles are easy to obtain, because they are a measure often used by process control engineers to describe the robustness of the closedloop system. The most popular choice is the unit circle, which is the boundary between stable and unstable. So they can be imposed as constraints into the optimization problem. 3.2. Estimates refinement Before proceeding to the main results, the following notations will be defined. Seen from (4) and (5), the closed-loop characteristic polynomial is A(t)S + B(t)R. Let Pclo (q−1 , t) denote the polynomial ˆ clo (q−1 , k, t) correequation, i.e., A(t)S + B(t)R = 0. Analogously, P sponds to the polynomial equation generated by the parameter estimates ˆ k (t). Suppose that all the closed-loop poles are falling in the disk with radius  on the z-plane. To be precise, it means that ri [P(q−1 , t)] ≤  for any time t and any i, where ri (•) is the modulus of the ith root of the polynomial equation (•). This information is Table 1 Identification comparison based on PTE and SI. Index

PTE SI

1

0

···

0

0

0

1

···

0

0

. . . 0

. . . 0

..

. . . 1

. . . 0

. ···

⎤

⎥ ⎥ ⎥, ⎥ ⎦

where hˆ 1,k (t), . . ., hˆ nh,k (t) are the coefficients of the polynomial H(t) defined in (5) with the parameters of A* (t), B(t) substituted by the estimate ˆ k (t). nh is the order of the polynomial, defined as ˆ k (t) is the controller canonical realizanh = max(na + ns, nb + nr). H ˆ clo (q−1 , k, t). tion of P Then, it is easy to see that ˆ k (t)] < . ˆ clo (q−1 , k, t)] <  ⇔ [H ri [P ˆ k (t) can be decomposed by eigenvalue decomposition, i.e., Since H ˆ k (t), it is obviˆ k (t) = U −1 U, containing all the eigenvalues of H H ˆ k (t)/ will satisfy ous that the matrix H ˆ k (t)/] < 1. [H Applying Lyapunov equation theorem, it follows that minˆ

Jk (t)

s.t.

ˆ k (t)/2 − M = −Q ˆ k (t)MH H

k (t),M,Q

T

(19)

M > 0, Q > 0. Here M > 0, Q > 0 mean that M, Q are positive definite matrices. The optimization problem (19) is yet nonconvex; in other words, it cannot be solved efficiently yet (in polynomial time), despite that the cost function and the constraints are all quadratic. This is because T ˆ k (t)MH ˆ k (t) in the constraint is comprised of the matrix product H

two optimization variables. To make it implementable online, we are interested in finding a sensible suboptimal solution to the optimization problem (19). The plausible way to construct such a suboptimal solution is based on the trust region method, a method widely used in nonlinear programming (Bertsekas, 1999; Hiriart-Urruty and Lemaréchal, 2013). The basic idea is, instead of searching the whole optimization variable space, to seek the solution in the proximity of the previous feasible solution. Mathematically, the optimization will be restated as minˆ

Jk (t)

s.t.

ˆ k (t) − H ˆ k−1 (t)M (t) ≤ k (t). H k

k (t)

(20)

ˆ k−1 (t), Mk (t)] is a feasible pair to the problem (19). X2 = Here [H M

Batch 1st

50th

23.0670 0.7214

60.4782 109.9085

Tr(XMXT ); it is defined in a Frobenius norm fashion. Recalled from ˆ k (t), the constraint in (20) can be further reduced the definition of H to ˆ k (t) − H ˆ k−1 (t)M (t) ≤ k (t) ⇔ ˆ H,k (t) − ˆ H,k−1 (t)M (t) ≤ k (t), H k−1 k

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

ˆ k (t) and H ˆ k−1 (t) are the same except the first line. ˆ H,k (t) is since H the vector that collects the elements {hˆ 1,k (t), . . ., hˆ nh,k (t)}. k (t) is the radius of the corresponding trust region; how to appropriately design it will be detailed later. One of the efficient way to solve (20) is the projection method. With ˆ H,k (t), Mk (t) and k (t), it automatically admits a domain (trust region) Dk (t) = {x ∈ Rnh : x − ˆ H,k−1 (t)Mk (t) ≤ k (t)} Thanks to the convexity of the norm and finite dimension of x, Dk (t) is a compact domain, thus making the projection possible. Henceforth, the mechanism can be summarized as follows. First, solving an unconstrained version of the problem (20) yields the solution denoted as ˆ ku (t), whose corresponding closed-loop coefficient vector is ˆ u (t). As a matter of fact, ˆ u (t) is generated recursively by H,k

k

the algorithm addressed in Section 2. The second step is to project u (t) on to the boundary of D (t) with Proj[ ˆ u (t)] as the proˆ H,k k H,k jection. It is tangible that Dk (t) is a norm ball; this fact makes the projection step quite easy, because the projection can be a convex u (t) and  ˆ H,k−1 (t), i.e., combination of ˆ H,k u Proj[ˆ H,k (t)] = ˆ H,k−1 (t) + ˛k (t)[ˆ H,k (t) − ˆ H,k−1 (t)]

(21)

u (t) ∈ D (t) ˆ H,k k

H,k

k (t) (t) − ˆ H,k−1 (t)M

k (t)

,

u (t) ∈ ˆ H,k / Dk (t).

(22)

u (t)]. In general, the The last step is to recover ˆ k (t) from Proj[ˆ H,k recovery step is not trivial, since it needs to solve a set of overdetermined equations. However, it is easy to verify that the map ϕ : ˆ k (t) → ˆ H,k (t) is a linear mapping, since the control blocks R and S are fixed. The linearity of ϕ renders us great convenience, because it preserves the topology between two spaces, shown in Fig. 5. This point can be seen as:

ϕ(ˆ k (t)) =

To this end, the missing piece of the whole puzzle is how to properly design the trust region radius k (t). We will take advantage of the techniques developed by the robust control community to resolve the problem. From (19), given a positive definite matrix Mk (t), it is known that only T

ˆ k (t)Mk (t)H ˆ k (t)/2 − Mk (t) < 0 H ˆ k (t) denote H ˆ k (t) − H ˆ k−1 (t). Then, the inequality is required. Let ıH above becomes T

ˆ k−1 (t) + ıH ˆ k (t)]Mk (t)[H ˆ k−1 (t) + ıH ˆ k (t)] − 2 Mk (t) < 0 [H

(23)

The given matrix Mk (t) is solved from the following Lyapunov equation ˆ k−1 (t)Mk (t)H ˆ k−1 (t) − 2 Mk (t) = −Q. H

⎧ ⎪ ⎨ 1,

⎪ ⎩ ˆ u

a linear programming; second, the step size (˛k (t)) in ours is tuned adaptively, whereas in FWA, it is chosen from a sum-divergent  sequence, i.e., 1/k. Additionally, ours also resembles to the adaptive solution filter proposed by Gao and Engell (2005), Marchetti et al. (2009) and Singhal et al. (2015). Their method has been proposed to solve a robust optimization problem rather than system identification problem.

T

with ˛k (t) ∈ (0, 1) defined as

˛k (t) =

227

ϕ([1 − ˛k (t)]ˆ k−1 (t) + ˛k (t)ˆ ku (t))

= [1 − ˛k (t)]ϕ(ˆ k−1 (t)) + ˛k (t)ϕ(ˆ ku (t)) = ˆ H,k (t) The second equality arises due to the linearity of ϕ; the second equality follows from the definition of ϕ. Remark 3.1. The approach developed above to solve the optimization problem is quite similar to the famous Frank–Wolfe algorithm (FWA). But there are two points that differentiate them: first, the new solution (ˆ ku (t)) in our algorithm is obtained from an unconstrained optimization problem, while in FWA, it is generated from

(24)

In this particular case, the positive definite matrix Q is fixed as the ˆ k−1 (t)/ can be guaridentity matrix to simplify calculation. Since H anteed to be a stable matrix by recursive arguments, it ensures the existence of Mk (t) to (24). Actually, Mk (t) can be solved via −1 ˆ k−1 (t) ⊗ H ˆ k−1 (t) − 2 I] vec(I). vec[Mk (t)] = −[H On the other hand, introducing a new parameter  ∈ (0, +∞), one can have that

 (



ˆ k−1 (t) −1 H √



ˆ k (t) − ıH

Mk (t) −





√

ˆ k (t) ıH

Mk (t))(



ˆ k−1 (t) −1 H

Mk (t)

T

Mk (t)) ≥ 0,

whose equivalent form is T

T

T

T

ˆ k−1 (t)Mk (t)H ˆ k−1 (t) + ıH ˆ k (t)Mk (t)ıH ˆ k (t) −1 H ˆ k (t)Mk (t)H ˆ k−1 (t) + H ˆ k−1 (t)Mk (t)ıH ˆ k (t) ≥ ıH Injecting the above equation into (23), (23) will hold if T

ˆ k (t)]Mk (t)[H ˆ k−1 (t) + ıH ˆ k (t)] − 2 Mk (t) ˆ k−1 (t) + ıH [H T

T

ˆ k−1 (t)Mk (t)H ˆ k−1 (t) + (1 + )ıH ˆ k (t)Mk (t)ıH ˆ k (t) ≤ (1 + −1 )H − 2 Mk (t) = −1 2 Mk (t) − (1 + −1 )I + (1 + )ˆ H,k (t) − ˆ H,k−1 (t)2M

k (t)

< 0.

The third equality follows from (24). A more restrictive version will be

−1 2 max [Mk (t)] − (1 + −1 ) + (1 + ) 2k (t) < 0,

(25)

There will be a solution k (t) to (25) if and only if

 > −1 + 2 max [Mk (t)],

Fig. 5. Map ϕ preserves the topology between two spaces.

since 2 Mk (t) > I from (24). Now k (t) is readily solvable. But in order to reduce the conservatism, it needs to enlarge the feasible domain, namely, to maximize k (t). Equivalently, let

228

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

Fig. 6. Workflow of the proposed algorithm.

 =−1 + 2 max [Mk (t)] + ˇ, ˇ > 0. Thereafter, 2k (t) will be a function of ˇ, signified as f(ˇ). Then,

max 2k (t) = maxf (ˇ) ˇ

ˇ

= max ˇ

 (−1 +

2 

ˇ 2 max [Mk (t)] + ˇ)( max [Mk (t)] + ˇ)

 .

To this end, it can be ensured that the imposed constraints will never be violated. The constraints formulated from the priori closed-loop knowledge can to some extent mitigate the undesired estimate variations. To take one step forward, the following adaptive filter is introduced to further suppress the variations.



f f 1 ˆ k (t) + (1 − 1 )ˆ k (t − 1), ˆ k (t) − ˆ k (t − 1) > 1 + 2 3k

f ˆ k (t) =

f

2 ˆ k (t) + (1 − 2 )ˆ k (t − 1), Otherwise

It can be directly concluded that

(27)

max 2k (t) =



22 max [Mk (t)] − 1 + 2

1 (2 max [Mk (t)] − 1)max [Mk (t)]

.

(26)

f

Here ˆ k (t) is the filtered parameter estimate. The parameter 1 ∈ (0, 1) in the top line of (27) should be selected relatively small so that the first filter will be granted a strong smooth effect

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

to remove the excessively severe “jumping”. Following this idea,

1 > 0 should be chosen a little larger than the maximum real dynamics drifting magnitude, which is required as a priori process knowledge. 2 > 0, 3 > 0 are two tuning knobs for transient performance enhancement. Since in the first several batches the identification moving steps are larger that in the following batches due to their larger distance to the true, 2 3k permits such large moving steps in the beginning. But this term should be decaying, because the identification performance is expected to improve gradually from batch to batch. Hence, 3 should be less than 1. Analogously, 2 should be chosen relatively large to capture the real dynamics change while neutralizing the variations with small magnitudes.

229

Fig. 7. The identification results of the parameter a1 (t) by the proposed method.

3.3. Summary The proposed identification method essentially employs a recursive unconstrained closed-loop identification algorithm as its “innovation” generator to steer the parameter estimates to fit the real dynamics. The estimates generated by the unconstrained algorithm is refined by projection and two adaptively low-pass filter in the sequel. The approach’s workflow is shown in Fig. 6. Interestingly, Fig. 6 also suggests the possibility to run in parallel; it further reduces the computation load. 4. Property analysis In this section, some theoretical analysis will be presented to the proposed identification algorithm, concerning about two crucial issues – the benefits of constraints imposed and recursive feasibility. Theorem 1. If the linear mapping ϕ : ˆ k (t) → ˆ H,k (t) is both oneto-one and onto, letting ˆ k−1 (t) − ˆ k (t − 1) = dk (t) be bounded, then there exists an upper bound on the variation ˆ k (t) − ˆ k (t − 1). Furthermore, if dk (t) ≤ k (t), then

ˆ k (t) − ˆ k (t − 1) ≤ ˆ ku (t) − ˆ k (t − 1). Proof.

(28)

According to the projection step, we have

ˆ k (t) − ˆ k (t − 1) = ˛k (t)[ˆ ku (t) − ˆ k−1 (t)] + [ˆ k−1 (t) − ˆ k (t − 1)].

Fig. 8. The spectral radius of the closed-loop characteristic polynomial under the proposed method.

ˆ k+1 (t)] <  as well. Note that all the filters in (27) are convex [H f (t − 1). Due to the linearity of ϕ, combinations of ˆ k+1 (t) and ˆ k+1

f it automatically generates ˆ H,k (t) = 1 ˆ H,k (t) + (1 − 1 )ˆ H,k (t − 1)

f or similarly for 2 corresponding to ˆ k (t). Then, it is directly from f Kharitonov theorem that ˆ (t) is a feasible solution. That com-

pletes the proof. 䊐

H,k

5. Illustration

From the triangle inequality, it follows that ˆ k (t) − ˆ k (t − 1) ≤ ˛k (t)ˆ ku (t) − ˆ k−1 (t) + dk (t). Since ϕ is bijective, there exists a positive constant ¯ such that for arbitrary 1 , 2 ∈ Rna+nb , 1 − 2  ≤ ϕ( ¯ 1 ) − ϕ(2 ). Then, the upper bound is

5.1. Revisit Example 1 Let us revisit Example 1 and apply the proposed method to it. In this case, k is selected as constant and equal to 0.9.  is chosen to be 0.95. 1 , 2 , 1 , 2 , 3 are respectively 0.2, 0.8, 1, 1.04, 0.98. Comparing Figs. 7 and 8 with Figs. 3 and 4, the parameter estimate tracks the true parameter trajectory closely; the corresponding

ˆ k (t) − ˆ k (t − 1) ≤ ˛ ¯ k (t) k (t) + dk (t) ≤ 

¯ k (t) + dk (t). Since our projection is linear and bijective, the projection on ˆ H,k (t) preserves on ˆ k (t). Because dk (t) ≤ k (t), ˆ k (t − 1) ∈ Dk (t). Then, (28) follows from Lemma 1 in the paper of Nedic´ et al. (2010). 䊐

Theorem 2 (Recursive feasibility). If the parameter estimates in the first batch are a set of feasible solutions, then the following estimates will never violate the imposed constraint. Proof. We will use inductive arguments to show it. It is trivial f to show the first batch. Assume that ˆ k (t) is a feasible solution. There must exist an Mk (t) to the Lyapunov equation (24). It follows that the projected estimate ˆ k+1 (t) satisfies the constraint

Fig. 9. MSPE comparison among the three methods on Example 1.

230

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

Fig. 11. MSPE comparison among the three methods on Example 2.

Fig. 10. PTE comparison among the three methods on Example 1.

closed-loop spectral radius is uniformly contained within the disk with the radius 0.95. Moreover, the abrupt “jumping” is significantly reduced by the proposed method. Fig. 9 shows that the performance of both the proposed method and the unconstrained method based on index (10) improves gradually from batch to batch, while the approach based on index (9) only maintains at a level of 0.31. From the aspect of mean-squared-prediction-error (MSPE), in the first several batches, the unconstrained method performs slightly better than the proposed method; this is understandable because the step size is restricted in our method. Finally, the proposed method outperforms the unconstrained one. Moreover, the proposed method is able to reduce PTE from batch to batch as well; but the constrained cannot according to Fig. 10.

Fig. 12. PTE comparison among the three methods on Example 2.

5.3. Example 3 5.2. Example 2 In this part, a more complicated example will be presented, adapted from Example 1 in Laudau and Karimi’s work (Landau and Karimi, 1997). The system is set up as follows:

⎧ 1 − 1.5z −1 + 0.7z −2 , ⎪ ⎨



t ∈ [0, 150)

0.32 A(z −1 ) = 1 + −1.5 + × (t − 150) z −1 + 0.7z −2 , t ∈ [150, 300) 150 ⎪

⎩

1 − 1.18z

 B(z −1 ) =

−1

+ 0.7z

−2

,

z −1 + 0.5z −2 ,

t ∈ [0, 150)

0.9z −1

t ∈ [150, 600].

+ 0.5z −2 ,

t ∈ [300, 600].

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 . k , , 1 , 2 , 1 , 2 , 3 are respectively chosen as 0.98, 0.97, 0.2, 0.7, 0.8, 1.25, 0.8. The reference signal is selected as a step change from 0.5 to −0.5 at time 150. The unrepeatable part is designed as PRBS with magnitude of 0.25. The white noise variance is 0.09. Fig. 11 shows a similar result as in Example 1. Fig. 12 shows that the proposed method is the best and the gap between the proposed method and the unconstrained one is quite significant. In Fig. 13, “conventional” stands for the unconstrained method and “1st” means the method based on index (9). Fig. 13 suggests the gradual performance improvement of the proposed one and the constrained one. Fig. 14 indicates that the estimates of our method can always satisfy the imposed constraint.

As a final example, we will test our algorithm on a two-tank process (Ikonen and Najim, 2001). The process is shown in Fig. 15; its dynamics is

 dY 1 (t) 1 = (Q (t) − k1 Y1 (t)) A1 dt   dY 2 (t) 1 = (k1 Y1 (t) − k2 Y2 (t)), A2 dt where the parameters above is listed in Table 2. The system is controlled by the proportional feedback with a gain −0.1. The sampling time is equal to 1. The variance of measurement noise is 0.04. The reference signal is a step change from 10 to 7 at time 300; its nonrepetitive component is PRBS with magnitude 0.5. k , , 1 , 2 , 1 ,

2 , 3 are chosen as 0.98, 1, 0.4, 0.8, 0.05, 0.078, 0.8 respectively. According to Figs. 16 and 17, the proposed method achieves a competing MSPE as the unconstrained one, but the index SI improves a lot, which matches our understanding of the process more. Fig. 18 demonstrates again that our method can always fulfill the imposed constraint. 5.4. Discussion It should be remarked that Examples 2 and 3 are two typical benchmarks for testing system identification method (Ikonen and Table 2 Parameter list of the two-tank system. Parameter

Physical meaning

Value

Y1 (t) Y2 (t) k1 k2 A1 A2

Tank 1 liquid level Tank 2 liquid level Tank 1 restriction coefficient Tank 2 restriction coefficient Tank 1 cross-section area Tank 2 cross-section area

n/a n/a 1 0.9 1 1

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

231

Fig. 13. Identification results comparison among the three methods on Example 2.

Najim, 2001; Landau and Karimi, 1997); so the simulation results demonstrate the potential for applying the proposed method to a real batch process. It is also noted that, in Example 1, it needs almost 50 batches to reach a steady state; it might be significant. However, the number can be reduced by choosing a relatively small i . Although in the previous derivations and analysis, it implicitly involves an underlying assumption that the dynamics within a batch should be highly repeatable, our method does not limit to that. The approach has the capability to track slow-drifting inter-batch dynamics, also by properly selecting i . Regarding the

fast-changing inter-batch dynamics, the approach may not be able to produce satisfactory results; but it is understandable, since fastchanging inter-batch dynamics undermines the repetitiveness of batch processes. The simple implementation of the projection step is the key to our method. Since the step relies upon a linear scaling, it suggests that our method may be extendable to some general model classes, i.e., ARMAX, OE, but not to BJ. Because in the BJ, the poles are decided by the multiplication of two polynomials, which is not linear in terms of estimated parameters.

232

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233

Fig. 14. The spectral radius evolution comparison on Example 2.

Fig. 18. The spectral radius evolution comparison on Example 3.

6. Conclusions An online closed-loop identification algorithm has been developed for batch processes in this paper. The algorithm incorporates the priori closed-loop process knowledge into the identification procedure to refine the parameter estimates. A trust region method is adopted to recursively solve the constrained optimization problem; the trust region size is obtained by the techniques developed in robust control. Three numerical examples have demonstrated the effectiveness of the proposed method. Acknowledgements

Fig. 15. The setup of the two-tank system.

This project is supported in part by Hong Kong Research Grant Council Project No. 612512, Guangzhou Scientific and Technology Project No. 1219007, the National Natural Science Foundation of China Project No. 61227005, and Guangdong Academician Workstation Project No. 2012B090500010, Guangdong Scientific and Technological Project No. 2013B011301012, and Guangdong Scientific and Technological Project No. 2013B051000004. References

Fig. 16. PTE comparison on Example 3.

Fig. 17. SI comparison on Example 3.

Åström, K.J., Wittenmark, B., 2013. Adaptive Control. Courier Corporation. Bertsekas, D.P., 1999. Nonlinear Programming. Athena Scientific. Camacho, E.F., Alba, C.B., 2013. Model Predictive Control. Springer Science & Business Media. Cao, Z., Yang, Y., Lu, J., Gao, F., 2014. Constrained two dimensional recursive least squares model identification for batch processes. J. Process Control 24, 871–879. Chen, H., Allgöwer, F., 1998. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica 34, 1205–1217. Chi, R., Hou, Z., Xu, J., 2008. Adaptive ILC for a class of discrete-time systems with iteration-varying trajectory and random initial condition. Automatica 44, 2207–2213. Craig, J.J., Hsu, P., Sastry, S.S., 1987. Adaptive control of mechanical manipulators. Int. J. Robot. Res. 6, 16–28. Forssell, U., Ljung, L., 1999. Closed-loop identification revisited. Automatica 35, 1215–1241. Gao, W., Engell, S., 2005. Iterative set-point optimization of batch chromatography. Comput. Chem. Eng. 29, 1401–1409. Gustavsson, I., Ljung, L., Söderström, T., 1977. Identification of processes in closed loop – identifiability and accuracy aspects. Automatica 13, 59–75. Hiriart-Urruty, J.-B., Lemaréchal, C., 2013. Convex Analysis and Minimization Algorithms I: Fundamentals, vol. 305. Springer Science & Business Media. Hjalmarsson, H., Gevers, M., De Bruyne, F., 1996. For model-based control design, closed-loop identification gives better performance. Automatica 32, 1659–1673. Huang, B., Shah, S.L., 1997. Closed-loop identification: a two step approach. J. Process Control 7, 425–438. Huang, B., Shah, S.L., 2012. Performance Assessment of Control Loops: Theory and Applications. Springer Science & Business Media. Huang, B., Ding, S.X., Qin, S.J., 2005. Closed-loop subspace identification: an orthogonal projection approach. J. Process Control 15, 53–66. Ikonen, E., Najim, K., 2001. Advanced Process Identification and Control, vol. 9. CRC Press. Ioannou, P.A., Sun, J., 2012. Robust Adaptive Control. Courier Corporation. Landau, I.D., Karimi, A., 1997. Recursive algorithms for identification in closed loop: a unified approach and evaluation. Automatica 8, 1499–1523. Landau, I.D., Karimi, A., 1999. A recursive algorithm for ARMAX model identification in closed loop. IEEE Trans. Autom. Control 44, 840–843.

Z. Cao et al. / Computers and Chemical Engineering 90 (2016) 222–233 Ljung, L., 1987. System Identification: Theory for the User. Prentice Hall, NJ. Ma, D.L., Braatz, R.D., 2003. Robust identification and control of batch processes. Comput. Chem. Eng. 27, 1175–1184. Marchetti, A., Chachuat, B., Bonvin, D., 2009. Modifier-adaptation methodology for real-time optimization. Ind. Eng. Chem. Res. 48, 6022–6033. Mayne, D.Q., Rawlings, J.B., Rao, C.V., Scokaert, P.O., 2000. Constrained model predictive control: stability and optimality. Automatica 36, 789–814. Morari, M., Lee, J.H., 1999. Model predictive control: past, present and future. Comput. Chem. Eng. 23, 667–682. Mosca, E., 1995. Optimal, Predictive, and Adaptive Control. ´ A., Ozdaglar, A., Parrilo, P., 2010. Constrained consensus and optimization in Nedic, multi-agent networks. IEEE Trans. Autom. Control 55, 922–938. Qin, S.J., Badgwell, T.A., 2003. A survey of industrial model predictive control technology. Control Eng. Pract. 11, 733–764. Qin, S.J., Ljung, L., 2003. Closed-loop Subspace Identification with Innovation Estimation.

233

Söderström, T., Stoica, P., 1988. System Identification. Prentice-Hall Inc. Singhal, M., Faulwasser, T., Bonvin, D., 2015. On handling cost gradient uncertainty in real-time optimization. In: Proceedings of 9th International Symposium on Advanced Control of Chemical Processes (ADCHEM). Slotine, J.-J.E., Li, W., 1987. On the adaptive control of robot manipulators. Int. J. Robot. Res. 6, 49–59. Tayebi, A., 2004. Adaptive iterative learning control for robot manipulators. Automatica 40, 1195–1203. Van Den Hof, P.M., Schrama, R.J., 1995. Identification and control – closed-loop issues. Automatica 31, 1751–1770. Yang, Y., Gao, F., 2000. Adaptive control of the filling velocity of thermoplastics injection molding. Control Eng. Pract. 8, 1285–1296. Zang, Z., Bitmead, R.R., Gevers, M., 1995. Iterative weighted least-squares identification and weighted LQG control design. Automatica 31, 1577–1594. Zhu, Y., 2001. Multivariable System Identification for Process Control. Elsevier.