- Email: [email protected]
Prior-knowledge-based subspace identification for batch processes
Journal of Process Control 82 (2019) 22–30
Contents lists available at ScienceDirect
Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont
Prior-knowledge-based subspace identification for batch processes Jie Hou a,b,∗ , Fengwei Chen c , Penghua Li a,b , Zhiqin Zhu a,b a
College of Automation, Chongqing University of Posts and Telecommunications, Chongqing 400065, China Key Laboratory of Industrial Internet of Things & Networked Control, Ministry of Education, Chongqing 400065, China c Department of Automation, School of Electrical Engineering and Automation, Wuhan University, Wuhan 430072, China b
a r t i c l e
i n f o
Article history: Received 15 July 2018 Received in revised form 30 June 2019 Accepted 4 July 2019 Available online 22 August 2019 Keywords: Subspace identification Prior knowledge Disturbance Batch process Variance analysis Consistent analysis
a b s t r a c t In this paper, a prior-knowledge-based subspace identification method (SIM) is proposed for batch processes subject to repeatable disturbances. The proposed method is a two-step procedure for state-space model identification: in the first step, the extended observability and triangular Toeplitz matrices are estimated simultaneously from a parity space of the experimental data and, based on which, the corresponding system matrices are retrieved in the second step. More specifically, A and C are retrieved from the estimated extended observability matrix, while B and D are retrieved from the estimated triangular Toeplitz matrix. The proposed method shows several superiorities in the following aspects. Firstly, it is able to provide unbiased parameter estimation in the presence of repeatable disturbances, thanks to the proposed difference operator which eliminates the disturbance effect. Secondly, it shows better robustness to measurement noise compared with the existing SIMs using parity space, due to the inherent instrumental variable mechanism and the new technique to build the instrument, which greatly enhance the estimated model efficiency/accuracy. Lastly, by taking the auxiliary static-gain information into account in the identification procedure, the variance properties of the parameters can be improved, especially for the system matrices B and D. All the above-mentioned developments are analyzed with strict mathematical proofs, along with two illustrative examples to confirm the effectiveness and merits. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction Batch processes are an important class of manufacturing techniques in various modern industries characterized by a reaction carried out over a finite duration, and have already been applied in the pharmaceutical, biotech, specialty chemical, and microelectronics industries [1], to mention a few. For efficient operation of batch processes such as advanced control [2–4], state estimation [5,6], and quality monitoring and fault diagnosis [7–9], an accurate and reliable process model is always desired. Therefore, identification of batch processes has attracted great attentions in the control community. The existing methods to infer models necessary for optimization and control can be categorized into two groups. One is the conventional mechanism modeling approach that establishes models based on the first principles that have been known to govern the system behavior [10]. While advancements in hardware and software are continuously increasing the available computational power and efficiency, the development and implementation of
∗ Corresponding author. E-mail addresses: [email protected] (J. Hou), [email protected] (F. Chen), [email protected] (P. Li), [email protected] (Z. Zhu). https://doi.org/10.1016/j.jprocont.2019.07.002 0959-1524/© 2019 Elsevier Ltd. All rights reserved.
mechanism modeling remains still very challenging [11]. The other is the so-called data-based modeling approach. Compared to the previous first-principle modeling approach, it is more amenable to yield batch process models that are simple enough to model-based optimization and control [11], given the increased availability of process data. Several data-driven methodologies have been applied to identify batch processes, e.g., artificial neural networks methods [12], recursive least squares methods [13–15], and subspace identification methods (SIMs) [11,16–19]. SIMs have been broadly developed to identify state-space models of various continuous multivariable systems [20–22], but only a few are dedicated to batch processes. For example, state-space models were established for the purposes of batch-to-batch control and monitoring [16,17], product quality and model predictive control [18], handling multi-rate and missing data [11], characterizing batch particulate processes [19]. Note that, for the above-mentioned SIMs to generate consistent estimation, it is generally required to perform identification tests without external disturbances, which is hard to meet in realistic batch processes where disturbances are inevitable. For example, wind turbine batch systems often suffer from periodic air turbulences during continuous operation [23–27]; batch injection molding processes often suffer from disturbances arising from the periodic mold cavity pres-
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
sure change [5,6]. Such disturbances give rise to biased estimation, when applying the existing SIMs including the aforementioned ones, where the disturbances are simply treated as stochastic noise. The disturbances of batch systems are batch- and time-varying, which need to be carefully dealt with for consistent estimation. It usually has certain special “structures”, however, fastly attenuates with time within a batch and slightly changes with batch in adjacent batches. In this perspective, the disturbance can be regarded as a periodic time-varying but batch-invariant disturbance. The special “structures” of repetitiveness for disturbance have been widely adopted for further improving the performance of estimation [5,6], control [23–25], and identification [26,27]. In this paper, we propose a consistent prior knowledge-based SIM for batch processes subject to repeatable disturbances as well as white measurement noises, with the main contributions and the novelties being stated below: 1. Firstly, motivated by the fact that the disturbance is repeatable from batch to batch in adjacent batches, a difference operator is proposed to eliminate the disturbance effect at the current batch by using the data from the previous batch. Also, to preserve robustness to the measurement noise, the instrumental variable technique is applied to mitigate the noise effect. A singular value decomposition (SVD) based method, which uses the disturbance-free dynamical data along the time direction, is then proposed to estimate the extended observability matrix and the triangular Toeplitz matrix simultaneously. After that, the system matrices A and C are retrieved from the estimated extended observability matrix. Note that the proposed method is an instrumental variable method, where the instrument is built by using not only past input/output data but also future input data, and thus can enhance the estimated model efficiency/accuracy, in contrast to the exsiting SIM using parity space, like SIMPCA [28] and SIMPCA-Wc [29], where only past input/output data are used. This advantage will be demonstrated with a strict proof. 2. Secondly, the proposed difference operator to eliminate the disturbance effect may enlarge the variances of the estimated quantities. To circumvent this problem, the auxiliary static-gain information is taken into account in the computation of the system matrices B and D from the estimated triangular Toeplitz matrix. The usefulness of this technique will be verified using the variance analysis and simulations in Sections 4 and 5, where it is shown that the proposed method leads to smaller variances of the estimated matrices B and D compared to SIMs where no prior knowledge is used. As a by-product, the variance analysis shown in the paper gives a hint to improve the variance properties of the conventional SIMs [30–36] by taken the priori knowledge into account. Furthermore, a method is presented to compute the unknown static-gain information from batch process data in the presence of repeatable disturbances. The remainder of this paper is organized as follows. In Section 2 the identification problem is presented. In Section 3 the proposed method is derived in detail, followed by an analysis of estimation consistency and error variances in Sections 4 and 5, respectively. Section 6 demonstrates the performance of the proposed method, followed by some conclusions in Section 7. 2. Problem statement Consider a batch process subject to arbitrary repeatable disturbances in the following state-space form
xk (t + 1) = Axk (t) + Buk (t) + wk (t) yk (t) = Cxk (t) + Duk (t) + vk (t)
(1)
23
where xk (t) ∈ Rnx , uk (t) ∈ Rnu , and yk (t) ∈ Rny are the state vector, the input signal, and the output signal at the tth time instant and the kth batch, respectively. vk (t) ∈ Rny and wk (t) ∈ Rny are the measurement noise and disturbance at the tth time instant and the kth batch, respectively. (A, B, C, D) are the system matrices of appropriate dimensions. Regarding the system in (1), we impose the following assumptions. • All the eigenvalues of A are strictly within the unit circle. • The system is controllable and observable. • The disturbance wk (t) is unknown but repeatable. Mathematically, it is only a function of time and k-invariant [5,6]. • vk (t) is assumed to be zero-mean white noise with covariance matrix E vj (i)v (t) = Qıj,k;i,t , where ıj,k;i,t is a two-dimensional k Kronecher delta function: ıj,k;i,t = I if and only if j = k and i = t, otherwise ıj,k;i,t = 0, and Q is a positive definite matrix. • vk (t) and uk (t) are mutually independent. The identification task is to estimate the deterministic system matrices (A, B, C, D) up to a similarity transformation from N input–output data {u(t), y(t)} collected at instants t = 1, . . ., N. 3. Proposed method 3.1. Disturbance elimination Since the disturbance is unknown but repeatable (as imposed in A3), its effect on the state and output of the system at the kth batch can be eliminated by using the information of the auxiliary dynamical data from the lth batch, l = / k. Define a batch difference operator between the kth and lth batch dynamical data rk (t) = r˜k (t) = rk (t) − rl (t)
(2)
where r can be x, y, u, v or w. It should be noted that wk (t) = 0, due to the assumption of repeatable disturbances. Then, by applying the operator to the system (1), we get a new formulation of the system without disturbances
x˜ k (t + 1) = A˜xk (t) + Bu˜ k (t)
(3)
y˜ k (t) = C x˜ k (t) + Du˜ k (t) + v˜ k (t)
where x˜ k (t) = xk (t). Under these assumptions, and based on (3), we have the following results Proposition 1. Under the assumptions A3–A4, v˜ k (t) is zero mean white noise with covariance matrix E
v˜ k (i)˜vk (t)
= E (vk (i) − vl (i)) (vk (t) − vl (t))
= E vk (i)v (t) − vk (i)v (t) − vl (i)v (t) + vl (i)v (t) k l k l
= E vk (i)v (t) + vl (i)v (t) k l
(4)
= E vk (i)v (t) + E vl (i)v (t) k l = 2Qıi,t
where ıi,t is the Kronecher delta function: ıi,t = I if and only if i = t, otherwise ıi,t = 0, and I is an identity matrix of appropriate dimension. 䊏 Proposition 2. Under the assumptions A3–A5, v˜ k (t)is mutuallyindependent of u˜ k (i) and x˜ k (i), i.e., E v˜ k (t)u˜ (i) = 0, E v˜ k (t)˜xk (i) = 0, k while v˜ k (t) is correlated with y˜ k (i), with the covariance matrix given by
24
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
to eliminate the noise Post-multiplying both sides of (21) by Z˜ k,p term, we obtain
E v˜ k (t)˜yk (i)
= E (vk (t) − vl (t)) (yk (i) − yl (i))
1 1 Z˜ Z˜ = N k,f k,p N
= E (vk (t) − vl (t)) (Cxk (i) + Duk (i) + vk (i) − Cxl (i) − Dul (i) − vl (i))
= E (vk (i) − vl (i)) (vk (t) − vl (t))
x˜ k,f = Ap x˜ k,p + Lp u˜ k,p
(7)
y˜ k,f = f x˜ k,f + Hf u˜ k,f + v˜ k,f
(8)
where p and f are the past and future window sizes, respectively. The extended controllability matrix Lp , the extended observability matrices f and p , and the triangular Toeplitz matrix Hf are defined as Lp = [Ap−1 B, . . ., B]
(9)
) )]
p = [C , . . ., (CAp−1 ) )]
⎡
D
⎢ CB ⎢ Hf = ⎢ . ⎣ .. CAf −2 B
(10) (11)
0
···
0
D
···
0
⎤
.. .
⎥ ⎥ ⎥ . .. . .. ⎦
CAf −3 B
···
u˜ k,f
I
1 N
+ Z˜ k,p
v˜ k,f Z˜ k,p .
0
⎤⎡ ⎦⎣
k,f
⎡ ⎤ 1 − lim (vk,f − vl,f )(v − v )⎣ k,p l,p N→∞ N
I
f
Hf
0
I
1 Z˜ k,f Z˜ k,p = USV = U1 N
rank
lim
1 N
N→∞
⎦ =0
U2
V1 S1 0 0
S2
V2
Z˜ k,f Z˜ k,p
= nx + fnu
(25)
= nx .
(26)
f
Hf
0
I
= fnu + rank(f ) = fnu + nx .
rank
u˜ k,p = [u˜ k,p (t − p), u˜ k,p (t − p + 1), . . ., u˜ k,p (t − p + N − 1)]
(13)
Therefore, combining (24), (25) and (27), we obtain
u˜ k,f = [u˜ k,f (t), u˜ k,f (t + 1), . . ., u˜ k,f (t + N − 1)]
(14)
u˜ k,f (t) = [u˜ (t), . . ., u˜ (t + f − 1)] . k k
(16)
The Hankel matrices y˜ k,p , v˜ k,p , y˜ k,f and v˜ k,f are similarly defined to u˜ k,p and u˜ k,f . The state sequences are defined as x˜ k,p = [˜xk (t − p), x˜ k (t − p + 1), . . ., x˜ k (t − p + N − 1)] x˜ k,f = [˜xk (t), x˜ k (t + 1), . . ., x˜ k (t + N − 1)] . 3.2. Estimation of f and Hf
Z˜ k,f = y˜ k,f u˜ k,f
(19)
f
Hf
(20)
0
I
= U1 S1 V1 Z˜ k,f Z˜ k,p
f
Hf
0
I
(28)
T = U1
(29)
where T is a nonsingular matrix that can be simply chosen as an identity matrix. Since U2 U1 = 0, it follows that
U2
f
Hf
0
I
= 0.
(30)
U2 = [U21 U22 ]
x˜ k,f u˜ k,f
(31)
corresponds to the first fn columns of U , while U where U21 y 2 22 corresponds to the last fnu columns of U2 . By substituting (31) into
(30), f and Hf 1 = D(CB) · · ·
CAf −2 B
can be estimated as
[39]
(8) can then be rewritten as 1 1 Z˜ = N k,f N
1
(27)
Let
By defining
lim
N→∞ N
(17) (18)
(24)
Based on the above results, the following rank condition can be established by using the Schur complement in [38]
(15)
(23)
which will be proved in Theorem 1 in next section. Simultaneously, it follows from assumption A2 that
D
u˜ k,p (t − p) = [u˜ (t − p), . . ., u˜ (t − 1)] k k
⎦
where S1 is the largest nx + fny eigenvectors of S. U1 and V1 contain the left and right eigenvectors of U and V , respectively. Note that
⎤
0
and Hp is defined similarly to Hf by replacing f by p. The input Hankel matrices are defined as
u˜ u˜ Z˜ k,p = y˜ k,p k,p k,f
(22)
by performing an SVD onto (22), we have
rank f (12)
x˜ k,f
1 1 v˜ Z˜ = lim (v − vl,f ) ⎣ N k,f k,p N→∞ N k,f u˜
N→∞
(6)
y˜ k,p = p x˜ k,p + Hp u˜ k,p + v˜ k,p
f = [C , . . ., (CA
I
x˜ k,f
⎡
where ıi,t is the Kronecher delta function. 䊏 By iterating (3), we have the following relations in SIMs [37]
0
is uncorrelated with v ˜ k,p as shown in Proposition 2, i.e., Since Z˜ k,p
lim
f −1
Hf
(5)
= 2Qıi,t
f
1 + N
I
0
v˜ k,f .
(21)
ˆ f = (U )⊥ 21
(32)
ˆ f 1 = ( )−1 H
(33)
J. Hou et al. / Journal of Process Control 82 (2019) 22–30 ⊥
) is a basis for the orthogonal complement of the row where (U21 , and space of U21
⎡
1
⎢ ⎢ ⎢ = ⎢. ⎢ .. ⎣
2
f
=
···
2
f −1
···
3
.. .
..
0
···
.
f
0
f
.. .
.. .
0
0
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
(34)
(35)
and U , respecwhere i and i are the ith column blocks of −U21 22 tively, as
=[
1 , . . .,
f]
(36)
U22 = [1 , . . ., f ].
(37)
3.3. Estimation of system matrices ˆ can be As we know, the consistent solution of matrices Cˆ and A ˆ f as follows [40,41] extracted directly from ˆ f (1 : ny , :) Cˆ =
(38)
Substituting (45) to (46) yields † ˆ f 1 ) − (Lˆ Lˆ 2 ) = Lˆ 2 vec(H 2
ˆ vec()
× L3 (Lˆ 2 Lˆ 2 )
[1 , . . ., f ]
−U21
†
ˆ f (1 : (f − 1)ny , :) ˆ f (1 + ny : fny , :) ˆ= A
(39)
† where Lˆ 2 = (Lˆ 2 Lˆ 2 ) according to (41).
Hf 1 = L1
L1 =
0
0
f (1 : (f − 1)ny , :)
D
.
B
2 ˆ f 1 − Lˆ 2 vec()) ˆ = min vec(H
(43)
s.t.vec(K) = L3 vec()
vec(Hˆ f 1 ) − Lˆ 2 vec()2 + L3 vec() − vec(K)
(44)
where is a suitably chosen Lagrange multiplier weight. Taking partial derivatives of J with respect to and , and equating them to zero, we have [31]
L3 (Lˆ 2 Lˆ 2 )
−1 L3
× L3 (Lˆ 2 Lˆ 2 ) ˆ = (Lˆ Lˆ 2 ) vec() 2
−1 Lˆ 2 .
(47)
ˆ can be retrieved from ˆ Finally, Bˆ and D
(48)
{yk (t)}N t=1
= D + CB + · · · + CAf B u1 +
f
(49)
CAi wk (t − i) + vk (t)
i=0
= D + CB + · · · + CAf B u2 +
f
(50)
CAi wl (t − i) + vl (t).
Then N {yk (t)}N t=1 − {yl (t)}t=1
= D + · · · + CAf B (u1 − u2 ) + vk (t) − vl (t)
−1
−1 ˆ f 1 ) − vec(K) Lˆ 2 vec(H
−1 −1 ˆ f 1 ) − (Lˆ Lˆ 2 ) Lˆ . Lˆ 2 vec(H 2 3
(45)
(46)
(51)
= K(u1 − u2 ) + vk (t) − vl (t) and N
N {yk (t)}N t=1 − {yl (t)}t=1
t=f
where Lˆ 2 = Inu ⊗ Lˆ 1 and L3 = Inu ⊗ [Iny , C(Inx − A)−1 ]. The operator vec stands for the operation of forming a vector from a matrix by stacking its columns one under another. ⊗ denotes the Kronecker product. By using the Lagrange multipliers method, the problem (43) becomes an optimization problem of the form
=
−1 ˆ f 1 ) − vec(K) Lˆ 2 vec(H
The static gain information is necessary to reduce the error variances of parameter estimates. In this subsection, a method is proposed to estimate the unknown static gain based on two step tests, each having a different input amplitude. Assume that the particular static data for input signals u1 and u2 (u1 = / u2 ) are N denoted by {yk (t)}N and {y (t)} , respectively. Under assumption l t=1 t=1 A1, choose a sufficiently large future horizon f for computation
(42)
Subsequently, the system matrices B and D can be obtained by solving following constrained least squares problem
−1
i=0
K = D + CB + · · · + CA∞ B = [Iny , C(Inx − A)−1 ].
min
−1 L3
3.4. Estimation of static gain
(41)
Assuming that the auxiliary static gain is known as a constant value K, the following constraint can be imposed to compute the system matrices of B and D
ˆ =
L3 (Lˆ 2 Lˆ 2 )
†
{yl (t)}N t=1
, =
ˆ f 1 ). vec(ˆ ls ) = Lˆ 2 vec(H
(40)
Iny
−1 L3
Remark 1. For SIMs without any prior knowledge of the static gain, the system matrices B and D can be estimated by solving a standard least squares problem as follows
where the operator (:) used for partitioning matrices is the same than that in Matlab. It follows from the expression of Hf1 that
where
25
= (N − f + 1)K(u1 − u2 ) +
N
(52) (vk (t) − vl (t)) .
t=f
Under the assumptions A3–A5, the static gain can then be estimated as
N
Kˆ ≈
t=f
N {yk (t)}N t=1 − {yl (t)}t=1
(N − f + 1)(u1 − u2 )
.
(53)
Finally, the complete algorithm of the proposed method presented in the previous sections is summarized below Algorithm 1. 1. 2. 3. 4. 5. 6.
Construct the matrices y˜ k,p , u˜ k,p , y˜ k,f , and u˜ k,f . Estimate static gain parameter K using (53). Eliminate the noise using (22). Estimate f and Hf1 using (32) and (33). Estimate A and C using (38) and (39). Estimate B and D using (47).
26
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
4. Consistency analysis The estimation of the parameters of f and Hf are consistent if the rank condition in (25) is satisfied. Once the estimated f and Hf are consistent, the resulting system matrices are also consistent (see [28,29]). Note that, although the proposed method, and the existing SIMs using parity space like SIMPCA and SIMPCA-Wc, share the same algorithm framework, but they use different instruments. The instruments in the proposed method, SIMPCA and SIMPCA
Theorem 1. Under the assumptions A1–A5 and the following conditions, the rank condition of (25) holds. 1. The input signal u˜ k (t) is persistently exciting of order n + p + f; 2. p ≥ nx and f ≥ nx . Substituting Z˜ k,p and Z˜ k,f into the left side of (25) yields
rank
lim
1 N
N→∞
Z˜ k,f Z˜ k,p
⎛ ⎜ ⎝
rank ⎜
1
lim
f
Hf
0
I
=
x˜ k,f
I
+
N→∞ N
u˜ k,f
0
⎤ ⎞ ⎡ ! y˜ k,p ⎥ ⎟ ⎢ ⎥ ⎟ v˜ k,f ⎢ ⎣ u˜ k,p ⎦ ⎠ . u˜ k,f (54)
Since the input and noise are jointly stationary and uncorrelated with each other, (54) can be rewritten as
rank
lim
N→∞
1 N
⎛ ⎜ 1 rank ⎜ lim ⎝N→∞ N
Z˜ k,f Z˜ k,p
f
Hf
0
I
x˜ k,f u˜ k,f
⎤ ⎞ ⎡ y˜ k,p ⎥ ⎟ ⎢ ⎥ ⎟ ⎢ ⎣ u˜ k,p ⎦ ⎠ .
(55)
u˜ k,f
rank
lim
N→∞
⎛
1 N
Z˜ k,f Z˜ k,p
⎡
x˜ k,p
⎤⎡
f
Hf
0
I
0
l = ⎣0
I
0⎦
0
0
0
l˜xk,p
u˜ k,f
(57)
which is of rank nx + pnu . Then the rank condition in SIMPCA and SIMPCA-Wc is nx rank
lim
N→∞
1 N
Z˜ Z˜ k,p
⎤ ⎞ ⎡
p
0
min(nx + fnu , nx + pnu ).
(58)
It is obvious that the rank condition for SIMPCA and SIMPCA-Wc to be consistent may not be fulfilled, and this explains why the new instrument used in the proposed method enhances the estimation efficiency/accuracy compared to that in SIMPCA and SIMPCA-Wc.
5. Variance analysis The estimation error variance properties of the system matrix are investigated in this section. Compared with the SIM without prior knowledge, for example, the one in (48), the proposed method is able to reduce the estimation error variance of . This result is strictly proved in the following theorem. Theorem 2. Under the assumptions A1–A5, the proposed method results in smaller variances of parameter estimates than SIMs without prior knowledge. That is
E tr vec(˜ ls )vec(˜ ls
)
(59)
˜ = vec() ˆ − vec() and vec(˜ ls ) = vec(ˆ ls ) − vec() where vec() stand for, respectively, the estimation errors of the proposed method and SIMs without prior knowledge. †
Ap
Lp
0
0
0
=
⎜ ⎥⎢ ⎥ ⎟⎢ 1 ⎢ ⎢ ⎥⎢ ⎥ ⎟ rank ⎜ lim ⎝N→∞ N ⎣ u˜ k,p ⎦ ⎣ u˜ k,p ⎦ ⎠ ⎣ 0 u˜ k,f
⎤
Hp
Proof. By substituting Lˆ 2 = L2 + L˜ 2 into Lˆ 2 which is defined in (47), we have
By substituting (6), (7) and (8) into (55), we have
p
˜ ˜ E tr vec()vec( )
=
⎡
−1/2
] , Z = [˜yk,p ] , and Z (ZZ /N) , u˜ u˜ u˜ Wc are [˜yk,p k,p k,f k,p respectively (see (24) in [28], and (28) in [29]). These instruments are all uncorrelated with noise, but only the one in the proposed method, constructed by using both the past and future data, can strictly guarantee the rank condition of (25) under some conditions, see the following theorem.
Proof.
Note that replacing the instrument Z˜ k,p in the proposed method with the instruments in SIMPCA and SIMPCA-Wc, the last term of (56) becomes
I
⎤
Hp
0
I
⎥ . 0⎦
0
I
(56)
The rank of the first two factors in the right-hand side of (56) are equal to nx + fnu , and nx + fnu if p nx and f nx , respectively. The rank of the third factor in the right-hand side of (56) is equal to nx + (p + f)nu if the input signal u˜ k (t) is persistently exciting of order nx + p + f. The rank of the last factor in the right-hand side of (56) is nx + (p + f)nu if p nx and f nx . The details can be found in [39,42–45]. From these results, it is clear that the rank of (56) is equal to nx + fnu , indicating that the rank conditions in (25) hold under the conditions in Theorem 1. This completes the proof. 䊏
† Lˆ 2 =
=
L2 + L˜ 2
L2 + L˜ 2
−1
L2 + L˜ 2
L2 L2 + L2 L˜ 2 + L˜ 2 L2 + L˜ 2 L˜ 2
−1
L2 + L˜ 2
(60) .
As is√well-known in asymptotic analysis, only the terms of order O(1/ N) (in probability) contribute to the asymptotic variance. √ Hence, since L˜ 2 converges to zero, their product is o(1/ N) and can be neglected [46]. With the following relation [47]
L2 L2 + L2 L˜ 2 + L˜ 2 L2 + L˜ 2 L˜ 2
(L2 L2 )
−1
− (L2 L2 )
−1
−1
(L2 L˜ 2 + L˜ 2 L2 )(L2 L2 )
−1
(61)
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
where we have used the symbol to denote that the matrix is approximated by its first order Taylor expansion. Further, (60) can readily be elaborated to † Lˆ 2
(L2 L2 )
−1
− (L2 L2 )
× L2 + L˜ 2 (L2 L2 )
†
−1 L2
− (L2 L2 )
−1
†
= L2 + L2 (L2 ) †
†
−1
+ (L2 L2 )
−1 L˜ 2
†
†
= L2 + L2 (L2 )
†
†
− L2 L˜ 2 L2 − L2 (L2 )
−1 L˜ 2 (I
(62)
−1 L˜ 2 (I
− L2 L2 )
†
†
†
†
−1 L˜ 2 (I
†
†
+ L2 (L2 )
−1 L˜ 2 (I
−1 L3
†
† † L2 (L2 ) Pr
†
(63)
†
†
.
†
†
†
†
E tr L2 (L2 ) Pr Pr
=
†
†
= E tr L2 (L2 )
− L2 L2 ) − L˜ 2 L2 vec(Hf 1 )
†
†
From (65) and (70), we finally have
= vec(H˜ f 1 ) +
−1 (L2 ) L˜ 2 (I
† † − L2 L2 ) − L˜ 2 L2
vec(Hf 1 ).
E
tr vec(˜ ls )vec(˜ ls )
= E tr (L2 L2 )
= E tr
−1 −1 L2 L2 (L2 L2 )
† † L2 (L2 )
(64)
(65)
.
= vec(ˆ ls ) − (Lˆ 2 Lˆ 2 )
−1 L3
L3 (Lˆ 2 Lˆ 2 )
−1 L3
× L3 vec(ˆ ls ) − vec(K)
−1 −1 = vec(ˆ ls ) − (Lˆ 2 Lˆ 2 ) L3 L3 (Lˆ 2 Lˆ 2 ) L3
−1 −1 vec(ˆ ls ) − (L2 L2 ) L3 L3 (L2 L2 ) L3
†
vec() + L2 − (L2 L2 )
% I
−1 L3
−1 − (L2 L2 ) L3
† × L2
= vec() + (I − Pr)
L3 (L2 L2 )
†
†
= −E tr PrL2 (L2 )
−1
(70) .
(71)
.
The inequality of Theorem 1 holds since the term in (71) is positive definite. This completes the proof. 䊏 Remark 2. The asymptotic variances of the estimated transfer functions can be computed by using the estimation errors of the system matrices. The detailed expressions can be found in [48–52], and thus they will not be reproduced in this paper. It is clear that under the same level of estimation errors of the system matrices A and C, the smaller variances of the estimated system matrices B and D, the lower variances of the transfer function can be observed. 6. Illustration
−1 (66)
× L3 vec(˜ ls )
vec() +
−1
× L3 (vec() + vec(˜ ls )) − vec(K)
† × L3 L2
−1
(69)
− E tr(vec(˜ ls )vec(˜ ls ) )
†
†
= −E tr Pr(L2 L2 )
= −E tr PrL2 (L2 )
˜ can be derived In a similar way, by substituting (62) into (47), vec() as ˆ vec()
˜ ˜ ) E tr(vec()vec( )
The variance of vec(˜ ls ) in the SIM without prior knowledge can then be computed as follows
− E tr PrL2 (L2 )
†
which is equal to the third term on the right side of (68), and thus
= vec() + L2 where
(68)
E tr L2 (L2 ) Pr
˜ ˜ E tr vec()vec( )
†
=
The fourth term on the right side of (68) can be written as E tr PrL2 (L2 ) Pr
†
†
†
†
− E tr PrL2 (L2 )
+ E tr PrL2 (L2 ) Pr
†
(67)
†
†
− E tr
− L2 L2 ) − L˜ 2 L2
†
−1 L3
†
† † ˜ f 1) L2 vec(Hf 1 ) + L2 vec(H
˜ ˜ tr vec()vec( )
˜ f 1) × vec(Hf 1 ) + vec(H
L3 (L2 L2 )
− L2 L2 ) − L˜ 2 L2 .
L2 + L2 (L2 )
= E tr L2 (L2 )
†
−1 L3
= E tr (I − Pr)L2 (L2 ) (I − Pr)
−1 † L˜ 2 L2 L2
ˆ f 1) = Lˆ 2 vec(H
Pr = (L2 L2 )
E
ˆ f 1 = Hf 1 + H ˜ f 1 into (48) and neglecting the Substituting (62) and H ˜ f 1 and L˜ 2 immediately give product of H vec(˜ ls )
where
is a projection matrix satisfying PrPr = Pr. With the above results ˜ in the proposed method can be in mind, the variance of vec() computed as follows
−1 (L2 L˜ 2 + L˜ 2 L2 )(L2 L2 ) L2
†
†
−1
−1 L˜ 2
= L2 − L2 L˜ 2 L2 + L2 (L2 ) †
(L2 L˜ 2 + L˜ 2 L2 )(L2 L2 )
27
−1 L3
−1 L3 (L2 L2 ) L3
−1
−1 & L3
In this section, two examples are used to demonstrate the effectiveness and merit of the proposed method. One is a benchmark example from [53] and the other is an industrial injection molding process studied in [54]. Example 1. Consider the following multiple-input, multipleoutput system which is the nominal model adopted for numerical simulation study in [53]:
xk (t + 1) = Axk (t) + Buk (t) + wk (t) yk (t) = Cxk (t) + Duk (t)u + vk (t)
(72)
28
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
Fig. 1. NMAE of estimated Bode magnitude plots using different methods for example 1.
Fig. 2. Variance of estimated Bode magnitude plots using different methods for example 1.
with
⎡
0.8
A = ⎣0
C=
0 0.5 0
−0.4
0.2
⎤
⎡
0
−0.5 ⎦ , B = ⎣ 0
0.3 0 0.5
0.5
0
0
1
,D =
0
without utilizing prior information [see (52)], and SIMPCA-Wc with
⎤
u˜ u˜ prior information (replacing y˜ k,p k,p k,f
−0.6 ⎦ ,
0.5 0 0.7139 0.1174
the same figures. From these results, we can conclude that
0.3131 0.2876
1 ˆ |Gij (ω) − Gij (ω)| 100
by Z ) are shown in
(73)
where the input signal uk (t) is chosen as a zero-mean, Gaussian, random sequence with unit variance. The noise vk (t) is assumed to be zero-mean, Gaussian, white noise of variance 0.5. The disturbance wk (t) is a sinusoid signal wk (t) = sin(2 /3). The statistical properties of the proposed method are investigated by means of a Monte-Carlo (MC) simulation of 100 runs. In each run, two sets of batch dynamical data are generated to obtain the disturbance-free data, each having a different input sequence. The static data are generated by using two constant inputs: u1 = 0.3 and u2 = 0.95. The performance index for the identification methods, i.e. the normalized mean absolute error (NMAE), is computed based on the frequency response data of the 100 estimated models. The NMAE is defined as 100
NMAE =
• Compared with the conventional SIMPCA-Wc where only past input/output data are used to build the instrument, the proposed method and the SIM without utilizing prior information generally enhance the estimated model efficiency/accuracy, thanks to the instrument which is built based on both the past input/output and future input data, see Section 4 for the reason on this improvement. • While the proposed method and the SIM without utilizing prior information give consistent estimated results, the proposed method reduces the variance of the estimated models, thanks to the prior knowledge of static gain which contains adequate information on the low-frequency band of the system. A detailed explanation can be found in Section 5.
Example 2. in [54]:
Consider an injection molding batch process studied
⎧ ⎨
(74)
i=1
where ω ∈ [0, ] is the frequency variable, and Gij (ω) is the magnitude from jth input to ith output. Letting N = 500 and f = p = 10, the NMAE and variance are compared in Figs. 1 and 2. For ease of comparison, the results of the SIM
xk (t + 1) =
1.582
−0.592
xk (t) + 1 0 ⎩ yk (t) = 1.69 1.419 xk (t) + vk (t)
1 0
uk (t) + wk (t)
(75)
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
Fig. 3. NMAE of estimated Bode magnitude plots using different methods for Example 2.
29
Subsequently, the system matrices A and C are extracted from the estimated extended controllability matrix, while the system matrices B and D are extracted from the estimated triangular Toeplitz matrix. The proposed method can be regarded as an improved SIM in relative to the existing SIM using parity space, e.g. SIMPCA and SIMPCA-Wc. Due to the instrument used in the proposed method, built based on the past and future input/output to eliminate the impact of measurement noise, is better correlated with system matrix than the ones used in SIMPCA and SIMPCA-Wc, the estimation efficiency/accuracy can be greatly enhanced. This property has been analyzed with a strict proof. As the second contribution of this paper, a technique has been proposed to retrieve the prior knowledge of the system static-gain information using the auxiliary static data, which helps to improves the quality and variance properties of the estimated B and D. Two illustrative examples have shown the effectiveness and merit of the proposed techniques and identification method. Acknowledgments
Fig. 4. Variance of estimated Bode magnitude plots using different methods for Example 2.
where the noise vk (t) is zero-mean, Gaussian and white, with variance of 0.5. The input is designed as a random binary sequence in Matlab as follows u = 0.5 + 0.3 ∗ idinput(Ts, rbs ) ; The disturbance wk (t) that accounts for the mold cavity pressure will affect the injection velocity. Just as an example, it is assumed that wk (t) has the expression wk (t) =
−a1 z −1 − a−2 2 1 − b1 z −1
*
where b1 =
*
0.4
0.3 t < 5000 0.4 t ≥ 5000
(76)
* , a1 =
0.05 t < 5000 0.15 t ≥ 5000
, and
0.08 t < 5000 . 0.15 t ≥ 5000 The statistical properties of the proposed method are illustrated by means of an MC simulation of 100 runs. The simulation settings remain the same as in the previous example, except for u1 = 0.8 and u2 = 0.95. With N = 2000 and f = p = 10, the NMAE and variance of the estimated models using the same identification methods presented in Example 1 are shown in Figs. 3 and 4. Again, the proposed method results in very accurate estimates with smaller variances than the other two identification methods, especially in the presence of repeatable disturbances. a2 =
7. Conclusion For identification of batch processes in the presence of repeatable disturbances as well as output measurement noises, a prior-knowledge-based subspace identification method has been proposed to achieve consistent estimation of the system parameters. As the first contribution of this paper, a new operator has been proposed to eliminate the influence of repeatable disturbances, with the underlying rationale being to collect one more data set to recover the knowledge of the disturbances. Then, by applying the proposed difference operator to eliminate the disturbances in the data, consistent estimates of the extended controllability and triangular Toeplitz matrices can be obtained using the proposed SVD-based method from the disturbance-free batch data.
This work was supported by the National Natural Science Foundation of China under grants 61903057, 61803061, 61703347, 61703066 and 61703311; the Chongqing Natural Science Foundation under grants cstc2019jcyj-msxmX0129, cstc2016jcyjA0428 and cstc2018jcyjAX0536; the Science and Technology Research Program of Chongqing Municipal Education Commission under grant KJQN201800603. References [1] D. Bonvin, B. Srinivasan, D. Hunkeler, Control and optimization of batch processes, IEEE Control Syst. Mag. 26 (6) (2006) 34–45. [2] S. Hao, T. Liu, W. Paszke, K. Galkowski, Robust iterative learning control for batch processes with input delay subject to time-varying uncertainties, IET Control Theory Appl. 10 (15) (2016) 1904–1915. [3] S. Hao, T. Liu, F. Gao, PI based indirect-type iterative learning control for batch processes with time-varying uncertainties: a 2D FM model based approach, J. Process Control 78 (2019) 57–67. [4] D. Li, S. He, Y. Xi, T. Liu, F. Gao, Y. Wang, J. Lu, The synthesis of ILC-MPC controller with data-driven approach for constrained batch processes, IEEE Trans. Ind. Electron. (2019), http://dx.doi.org/10.1109/TIE.2019.2910034. [5] Z. Cao, J. Lu, R. Zhang, F. Gao, Iterative learning Kalman filter for repetitive processes, J. Process Control 46 (2016) 92–104. [6] Z. Cao, R. Zhang, Y. Yang, J. Lu, F. Gao, Discrete-time robust iterative learning Kalman filtering for repetitive processes, IEEE Trans. Autom. Control 61 (1) (2016) 270–275. [7] C. Zhao, F. Wang, M. Jia, Dissimilarity analysis based batch process monitoring using moving windows, AIChE J. 53 (5) (2007) 1267–1277. [8] J. Liu, T. Liu, J. Chen, Window-based stepwise sequential phase partition for nonlinear batch process monitoring, Ind. Eng. Chem. Res. 55 (34) (2016) 9229–9243. [9] J. Liu, T. Liu, J. Chen, Sequential local-based Gaussian mixture model for monitoring multiphase batch processes, Chem. Eng. Sci. 181 (2018) 101–113. [10] D. Shi, N.H. El-Farra, M. Li, P. Mhaskar, P.D. Christofides, Predictive control of particle size distribution in particulate processes, Chem. Eng. Sci. 61 (1) (2006) 268–281. [11] M.M. Rashid, P. Mhaskar, C.L.E. Swartz, Handling multi-rate and missing data in variable duration economic model predictive control of batch processes, Chem. Eng. Sci. 63 (7) (2017) 2705–2718. [12] A.Y.D. Tsen, D.S.H. Jang, S.S. Wong, B. Joseph, Predictive control of quality in batch polymerization using hybrid ANN models, AIChE J. 42 (2) (1996) 455–465. [13] Z. Cao, J. Lu, R. Zhang, F. Gao, Online average-based system modelling method for batch process, Comput. Chem. Eng. 108 (4) (2018) 128–138. [14] Z. Cao, Y. Yang, H. Yi, F. Gao, Priori knowledge-based online batch-to-batch identification in a closed loop and an application to injection molding, Ind. Eng. Chem. Res. 55 (32) (2016) 8818–8829. [15] Z. Cao, Y. Yang, J. Lu, F. Gao, Constrained two dimensional recursive least squares model identification for batch processes, J. Process Control 24 (6) (2014) 871–879. [16] A.W. Dorsey, J.H. Lee, Building inferential prediction models of batch processes using subspace identification, J. Process Control 13 (5) (2003) 397–406. [17] Y. Yao, F. Gao, Subspace identification for two-dimensional dynamic batch process statistical monitoring, Chem. Eng. Sci. 63 (13) (2008) 3411–3418. [18] B. Corbett, P. Mhaskar, Subspace identification for data-driven modeling and quality control of batch processes, AIChE J. 62 (5) (2016) 1581–1601.
30
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
[19] A. Garg, P. Mhaskar, Subspace identification based modeling and control of batch particulate processes, Ind. Eng. Chem. Res. 56 (26) (2017) 7491–7502. [20] S.J. Qin, An overview of subspace identification, Comput. Chem. Eng. 30 (10-12) (2006) 1502–1513. [21] G. van der Veen, J.W. van Wingerden, M. Bergamasco, M. Lovera, M. Verhaegen, Closed-loop subspace identification methods: an overview, IET Control Theory Appl. 7 (10) (2013) 1339–1358. [22] T. Katayama, Subspace Methods for System Identification, Springer Science & Business Media, Berlin, 2006. [23] S.T. Navalkar, J.W. van Wingerden, E. van Solingen, T. Oomen, E. Pasterkamp, G.A.M. Van Kuik, Subspace predictive repetitive control to mitigate periodic loads on large scale wind turbines, Mechatronics 24 (8) (2014) 916–925. [24] I. Houtzager, J.W. van Wingerden, M. Verhaegen, Rejection of periodic wind disturbances on a smart rotor test section using lifted repetitive control, IEEE Trans. Control Syst. Technol. 21 (2) (2013) 347–359. [25] I. Houtzager, J.W. van Wingerden, M. Verhaegen, Wind turbine load reduction by rejecting the periodic load disturbances, Wind Energy 16 (2) (2013) 235–256. [26] P.M. Gebraad, J.W. van Wingerden, P.A. Fleming, A.D. Wright, LPV identification of wind turbine rotor vibrational dynamics using periodic disturbance basis functions, IEEE Trans. Control Syst. Technol. 21 (4) (2013) 1183–1190. [27] N.E. Goodzeit, M.Q. Phan, System identification in the presence of completely unknown periodic disturbances, J. Guid. Control Dyn. 23 (2) (2000) 251–259. [28] J. Wang, S.J. Qin, A new subspace identification approach based on principal component analysis, J. Process Control 12 (8) (2002) 841–855. [29] J. Wang, S.J. Qin, Closed-loop subspace identification using the parity space, Automatica 42 (2) (2006) 315–320. [30] P. Trnka, V. Havlena, Subspace like identification incorporating prior information, Automatica 45 (4) (2009) 1086–1091. [31] A. Alenany, H. Shang, M. Soliman, I. Ziedan, Improved subspace identification with prior information using constrained least squares, IET Control Theory Appl. 5 (13) (2001) 1568–1576. [32] A. Alenany, H. Shang, Recursive subspace identification with prior information using the constrained least squares approach, Comput. Chem. Eng. 54 (11) (2013) 174–180. [33] I. Markovsky, Mercère, Subspace identification with constraints on the impulse response, Int. J. Control 90 (8) (2017) 1728–1735. [34] Y. Wang, L. Zhang, Y. Zhao, Improved closed-loop subspace identification with prior information, Int. J. Syst. Sci. 49 (9) (2018) 1821–1835. ˇ ˇ L. Ferkl, M. Sebek, [35] S. Privara, J. Cigler, Z. Vána, Subspace identification of poorly excited industrial systems, 49th IEEE Conference on Decision and Control (CDC) (2010) 4405–4410. [36] K. Li, H. Luo, C. Yang, S. Yin, Subspace-aided closed-loop system identification with application to DC motor system, IEEE Trans. Ind. Electron. (2019), http:// dx.doi.org/10.1109/TIE. 2019.2907447.
[37] P. van Overschee, B.L. De Moor, Subspace Identification for Linear Systems: Theory-Implementation-Applications, Springer Science & Business Media, 2012. [38] M. Jansson, B. Wahlberg, On consistency of subspace methods for system identification, Automatica 34 (12) (1998) 1507–1519. [39] J. Hou, F. Chen, P. Li, Z. Zhu, F. Liu, An improved consistent subspace identification method using parity space for state-space models, Int. J. Control Autom. Syst. 17 (5) (2019) 1167–1176. [40] C. Yu, M. Verhaegen, Subspace identification of distributed clusters of homogeneous systems, IEEE Trans. Autom. Control 62 (1) (2016) 463–468. [41] M. Inoue, Subspace identification with moment matching, Automatica 99 (2019) 22–32. [42] J. Hou, T. Liu, B. Wahlberg, M. Jansson, Subspace Hammerstein model identification under periodic disturbance, IFAC-PapersOnLine 51 (15) (2018) 335–340. [43] J. Hou, T. Liu, Closed-loop subspace model identification using innovation estimation and orthogonal projection, Acta Autom. Sin. 42 (11) (2016) 1657–1663. [44] J. Hou, T. Liu, Q.G. Wang, Subspace identification of Hammerstein-type nonlinear systems subject to unknown periodic disturbance, Int. J. Control (2019), http://dx.doi.org/10.1080/00207179.2019.1621385. [45] C. Yu, L. Ljung, M. Verhaegen, Identification of structured state-space models, Automatica 90 (2018) 54–61. [46] A. Chiuso, G. Picci, Asymptotic variance of subspace methods by data orthogonalization and model decoupling: a comparative analysis, Automatica 40 (10) (2004) 1705–1717. [47] M. D“ohler, L. Mevel, Efficient multi-order uncertainty computation for stochastic subspace identification, Mech. Syst. Signal Process. 38 (2) (2013) 346–366. [48] M. Jansson, Asymptotic variance analysis of subspace identification methods, IFAC Proc. Vol. 33 (15) (2000) 91–96. [49] A. Chiuso, G. Picci, The asymptotic variance of subspace estimates, J. Econom. 118 (1-2) (2004) 257–291. [50] J. Hou, T. Liu, Q.G. Wang, Recursive subspace identification subject to relatively slow time-varying load disturbance, Int. J. Control 91 (3) (2018) 622–638. [51] J. Hou, F. Chen, P. Li, Z. Zhu, Fixed point iteration-based subspace identification of Hammerstein state-space models, IET Control Theory Appl. 13 (8) (2019) 1173–1181. [52] J. Hou, T. Liu, F. Chen, Orthogonal projection based subspace identification against colored noise, Control Theory Technol. 15 (1) (2017) 69–77. [53] G. Mercère, L. Bako, S. Lecœuche, Propagator-based methods for recursive subspace model identification, Signal Process. 88 (3) (2008) 468–491. [54] T. Liu, B. Huang, S.J. Qin, Bias-eliminated subspace model identification under time-varying deterministic type load disturbance, J. Process Control 25 (2015) 41–49.
Contents lists available at ScienceDirect
Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont
Prior-knowledge-based subspace identification for batch processes Jie Hou a,b,∗ , Fengwei Chen c , Penghua Li a,b , Zhiqin Zhu a,b a
College of Automation, Chongqing University of Posts and Telecommunications, Chongqing 400065, China Key Laboratory of Industrial Internet of Things & Networked Control, Ministry of Education, Chongqing 400065, China c Department of Automation, School of Electrical Engineering and Automation, Wuhan University, Wuhan 430072, China b
a r t i c l e
i n f o
Article history: Received 15 July 2018 Received in revised form 30 June 2019 Accepted 4 July 2019 Available online 22 August 2019 Keywords: Subspace identification Prior knowledge Disturbance Batch process Variance analysis Consistent analysis
a b s t r a c t In this paper, a prior-knowledge-based subspace identification method (SIM) is proposed for batch processes subject to repeatable disturbances. The proposed method is a two-step procedure for state-space model identification: in the first step, the extended observability and triangular Toeplitz matrices are estimated simultaneously from a parity space of the experimental data and, based on which, the corresponding system matrices are retrieved in the second step. More specifically, A and C are retrieved from the estimated extended observability matrix, while B and D are retrieved from the estimated triangular Toeplitz matrix. The proposed method shows several superiorities in the following aspects. Firstly, it is able to provide unbiased parameter estimation in the presence of repeatable disturbances, thanks to the proposed difference operator which eliminates the disturbance effect. Secondly, it shows better robustness to measurement noise compared with the existing SIMs using parity space, due to the inherent instrumental variable mechanism and the new technique to build the instrument, which greatly enhance the estimated model efficiency/accuracy. Lastly, by taking the auxiliary static-gain information into account in the identification procedure, the variance properties of the parameters can be improved, especially for the system matrices B and D. All the above-mentioned developments are analyzed with strict mathematical proofs, along with two illustrative examples to confirm the effectiveness and merits. © 2019 Elsevier Ltd. All rights reserved.
1. Introduction Batch processes are an important class of manufacturing techniques in various modern industries characterized by a reaction carried out over a finite duration, and have already been applied in the pharmaceutical, biotech, specialty chemical, and microelectronics industries [1], to mention a few. For efficient operation of batch processes such as advanced control [2–4], state estimation [5,6], and quality monitoring and fault diagnosis [7–9], an accurate and reliable process model is always desired. Therefore, identification of batch processes has attracted great attentions in the control community. The existing methods to infer models necessary for optimization and control can be categorized into two groups. One is the conventional mechanism modeling approach that establishes models based on the first principles that have been known to govern the system behavior [10]. While advancements in hardware and software are continuously increasing the available computational power and efficiency, the development and implementation of
∗ Corresponding author. E-mail addresses: [email protected] (J. Hou), [email protected] (F. Chen), [email protected] (P. Li), [email protected] (Z. Zhu). https://doi.org/10.1016/j.jprocont.2019.07.002 0959-1524/© 2019 Elsevier Ltd. All rights reserved.
mechanism modeling remains still very challenging [11]. The other is the so-called data-based modeling approach. Compared to the previous first-principle modeling approach, it is more amenable to yield batch process models that are simple enough to model-based optimization and control [11], given the increased availability of process data. Several data-driven methodologies have been applied to identify batch processes, e.g., artificial neural networks methods [12], recursive least squares methods [13–15], and subspace identification methods (SIMs) [11,16–19]. SIMs have been broadly developed to identify state-space models of various continuous multivariable systems [20–22], but only a few are dedicated to batch processes. For example, state-space models were established for the purposes of batch-to-batch control and monitoring [16,17], product quality and model predictive control [18], handling multi-rate and missing data [11], characterizing batch particulate processes [19]. Note that, for the above-mentioned SIMs to generate consistent estimation, it is generally required to perform identification tests without external disturbances, which is hard to meet in realistic batch processes where disturbances are inevitable. For example, wind turbine batch systems often suffer from periodic air turbulences during continuous operation [23–27]; batch injection molding processes often suffer from disturbances arising from the periodic mold cavity pres-
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
sure change [5,6]. Such disturbances give rise to biased estimation, when applying the existing SIMs including the aforementioned ones, where the disturbances are simply treated as stochastic noise. The disturbances of batch systems are batch- and time-varying, which need to be carefully dealt with for consistent estimation. It usually has certain special “structures”, however, fastly attenuates with time within a batch and slightly changes with batch in adjacent batches. In this perspective, the disturbance can be regarded as a periodic time-varying but batch-invariant disturbance. The special “structures” of repetitiveness for disturbance have been widely adopted for further improving the performance of estimation [5,6], control [23–25], and identification [26,27]. In this paper, we propose a consistent prior knowledge-based SIM for batch processes subject to repeatable disturbances as well as white measurement noises, with the main contributions and the novelties being stated below: 1. Firstly, motivated by the fact that the disturbance is repeatable from batch to batch in adjacent batches, a difference operator is proposed to eliminate the disturbance effect at the current batch by using the data from the previous batch. Also, to preserve robustness to the measurement noise, the instrumental variable technique is applied to mitigate the noise effect. A singular value decomposition (SVD) based method, which uses the disturbance-free dynamical data along the time direction, is then proposed to estimate the extended observability matrix and the triangular Toeplitz matrix simultaneously. After that, the system matrices A and C are retrieved from the estimated extended observability matrix. Note that the proposed method is an instrumental variable method, where the instrument is built by using not only past input/output data but also future input data, and thus can enhance the estimated model efficiency/accuracy, in contrast to the exsiting SIM using parity space, like SIMPCA [28] and SIMPCA-Wc [29], where only past input/output data are used. This advantage will be demonstrated with a strict proof. 2. Secondly, the proposed difference operator to eliminate the disturbance effect may enlarge the variances of the estimated quantities. To circumvent this problem, the auxiliary static-gain information is taken into account in the computation of the system matrices B and D from the estimated triangular Toeplitz matrix. The usefulness of this technique will be verified using the variance analysis and simulations in Sections 4 and 5, where it is shown that the proposed method leads to smaller variances of the estimated matrices B and D compared to SIMs where no prior knowledge is used. As a by-product, the variance analysis shown in the paper gives a hint to improve the variance properties of the conventional SIMs [30–36] by taken the priori knowledge into account. Furthermore, a method is presented to compute the unknown static-gain information from batch process data in the presence of repeatable disturbances. The remainder of this paper is organized as follows. In Section 2 the identification problem is presented. In Section 3 the proposed method is derived in detail, followed by an analysis of estimation consistency and error variances in Sections 4 and 5, respectively. Section 6 demonstrates the performance of the proposed method, followed by some conclusions in Section 7. 2. Problem statement Consider a batch process subject to arbitrary repeatable disturbances in the following state-space form
xk (t + 1) = Axk (t) + Buk (t) + wk (t) yk (t) = Cxk (t) + Duk (t) + vk (t)
(1)
23
where xk (t) ∈ Rnx , uk (t) ∈ Rnu , and yk (t) ∈ Rny are the state vector, the input signal, and the output signal at the tth time instant and the kth batch, respectively. vk (t) ∈ Rny and wk (t) ∈ Rny are the measurement noise and disturbance at the tth time instant and the kth batch, respectively. (A, B, C, D) are the system matrices of appropriate dimensions. Regarding the system in (1), we impose the following assumptions. • All the eigenvalues of A are strictly within the unit circle. • The system is controllable and observable. • The disturbance wk (t) is unknown but repeatable. Mathematically, it is only a function of time and k-invariant [5,6]. • vk (t) is assumed to be zero-mean white noise with covariance matrix E vj (i)v (t) = Qıj,k;i,t , where ıj,k;i,t is a two-dimensional k Kronecher delta function: ıj,k;i,t = I if and only if j = k and i = t, otherwise ıj,k;i,t = 0, and Q is a positive definite matrix. • vk (t) and uk (t) are mutually independent. The identification task is to estimate the deterministic system matrices (A, B, C, D) up to a similarity transformation from N input–output data {u(t), y(t)} collected at instants t = 1, . . ., N. 3. Proposed method 3.1. Disturbance elimination Since the disturbance is unknown but repeatable (as imposed in A3), its effect on the state and output of the system at the kth batch can be eliminated by using the information of the auxiliary dynamical data from the lth batch, l = / k. Define a batch difference operator between the kth and lth batch dynamical data rk (t) = r˜k (t) = rk (t) − rl (t)
(2)
where r can be x, y, u, v or w. It should be noted that wk (t) = 0, due to the assumption of repeatable disturbances. Then, by applying the operator to the system (1), we get a new formulation of the system without disturbances
x˜ k (t + 1) = A˜xk (t) + Bu˜ k (t)
(3)
y˜ k (t) = C x˜ k (t) + Du˜ k (t) + v˜ k (t)
where x˜ k (t) = xk (t). Under these assumptions, and based on (3), we have the following results Proposition 1. Under the assumptions A3–A4, v˜ k (t) is zero mean white noise with covariance matrix E
v˜ k (i)˜vk (t)
= E (vk (i) − vl (i)) (vk (t) − vl (t))
= E vk (i)v (t) − vk (i)v (t) − vl (i)v (t) + vl (i)v (t) k l k l
= E vk (i)v (t) + vl (i)v (t) k l
(4)
= E vk (i)v (t) + E vl (i)v (t) k l = 2Qıi,t
where ıi,t is the Kronecher delta function: ıi,t = I if and only if i = t, otherwise ıi,t = 0, and I is an identity matrix of appropriate dimension. 䊏 Proposition 2. Under the assumptions A3–A5, v˜ k (t)is mutuallyindependent of u˜ k (i) and x˜ k (i), i.e., E v˜ k (t)u˜ (i) = 0, E v˜ k (t)˜xk (i) = 0, k while v˜ k (t) is correlated with y˜ k (i), with the covariance matrix given by
24
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
to eliminate the noise Post-multiplying both sides of (21) by Z˜ k,p term, we obtain
E v˜ k (t)˜yk (i)
= E (vk (t) − vl (t)) (yk (i) − yl (i))
1 1 Z˜ Z˜ = N k,f k,p N
= E (vk (t) − vl (t)) (Cxk (i) + Duk (i) + vk (i) − Cxl (i) − Dul (i) − vl (i))
= E (vk (i) − vl (i)) (vk (t) − vl (t))
x˜ k,f = Ap x˜ k,p + Lp u˜ k,p
(7)
y˜ k,f = f x˜ k,f + Hf u˜ k,f + v˜ k,f
(8)
where p and f are the past and future window sizes, respectively. The extended controllability matrix Lp , the extended observability matrices f and p , and the triangular Toeplitz matrix Hf are defined as Lp = [Ap−1 B, . . ., B]
(9)
) )]
p = [C , . . ., (CAp−1 ) )]
⎡
D
⎢ CB ⎢ Hf = ⎢ . ⎣ .. CAf −2 B
(10) (11)
0
···
0
D
···
0
⎤
.. .
⎥ ⎥ ⎥ . .. . .. ⎦
CAf −3 B
···
u˜ k,f
I
1 N
+ Z˜ k,p
v˜ k,f Z˜ k,p .
0
⎤⎡ ⎦⎣
k,f
⎡ ⎤ 1 − lim (vk,f − vl,f )(v − v )⎣ k,p l,p N→∞ N
I
f
Hf
0
I
1 Z˜ k,f Z˜ k,p = USV = U1 N
rank
lim
1 N
N→∞
⎦ =0
U2
V1 S1 0 0
S2
V2
Z˜ k,f Z˜ k,p
= nx + fnu
(25)
= nx .
(26)
f
Hf
0
I
= fnu + rank(f ) = fnu + nx .
rank
u˜ k,p = [u˜ k,p (t − p), u˜ k,p (t − p + 1), . . ., u˜ k,p (t − p + N − 1)]
(13)
Therefore, combining (24), (25) and (27), we obtain
u˜ k,f = [u˜ k,f (t), u˜ k,f (t + 1), . . ., u˜ k,f (t + N − 1)]
(14)
u˜ k,f (t) = [u˜ (t), . . ., u˜ (t + f − 1)] . k k
(16)
The Hankel matrices y˜ k,p , v˜ k,p , y˜ k,f and v˜ k,f are similarly defined to u˜ k,p and u˜ k,f . The state sequences are defined as x˜ k,p = [˜xk (t − p), x˜ k (t − p + 1), . . ., x˜ k (t − p + N − 1)] x˜ k,f = [˜xk (t), x˜ k (t + 1), . . ., x˜ k (t + N − 1)] . 3.2. Estimation of f and Hf
Z˜ k,f = y˜ k,f u˜ k,f
(19)
f
Hf
(20)
0
I
= U1 S1 V1 Z˜ k,f Z˜ k,p
f
Hf
0
I
(28)
T = U1
(29)
where T is a nonsingular matrix that can be simply chosen as an identity matrix. Since U2 U1 = 0, it follows that
U2
f
Hf
0
I
= 0.
(30)
U2 = [U21 U22 ]
x˜ k,f u˜ k,f
(31)
corresponds to the first fn columns of U , while U where U21 y 2 22 corresponds to the last fnu columns of U2 . By substituting (31) into
(30), f and Hf 1 = D(CB) · · ·
CAf −2 B
can be estimated as
[39]
(8) can then be rewritten as 1 1 Z˜ = N k,f N
1
(27)
Let
By defining
lim
N→∞ N
(17) (18)
(24)
Based on the above results, the following rank condition can be established by using the Schur complement in [38]
(15)
(23)
which will be proved in Theorem 1 in next section. Simultaneously, it follows from assumption A2 that
D
u˜ k,p (t − p) = [u˜ (t − p), . . ., u˜ (t − 1)] k k
⎦
where S1 is the largest nx + fny eigenvectors of S. U1 and V1 contain the left and right eigenvectors of U and V , respectively. Note that
⎤
0
and Hp is defined similarly to Hf by replacing f by p. The input Hankel matrices are defined as
u˜ u˜ Z˜ k,p = y˜ k,p k,p k,f
(22)
by performing an SVD onto (22), we have
rank f (12)
x˜ k,f
1 1 v˜ Z˜ = lim (v − vl,f ) ⎣ N k,f k,p N→∞ N k,f u˜
N→∞
(6)
y˜ k,p = p x˜ k,p + Hp u˜ k,p + v˜ k,p
f = [C , . . ., (CA
I
x˜ k,f
⎡
where ıi,t is the Kronecher delta function. 䊏 By iterating (3), we have the following relations in SIMs [37]
0
is uncorrelated with v ˜ k,p as shown in Proposition 2, i.e., Since Z˜ k,p
lim
f −1
Hf
(5)
= 2Qıi,t
f
1 + N
I
0
v˜ k,f .
(21)
ˆ f = (U )⊥ 21
(32)
ˆ f 1 = ( )−1 H
(33)
J. Hou et al. / Journal of Process Control 82 (2019) 22–30 ⊥
) is a basis for the orthogonal complement of the row where (U21 , and space of U21
⎡
1
⎢ ⎢ ⎢ = ⎢. ⎢ .. ⎣
2
f
=
···
2
f −1
···
3
.. .
..
0
···
.
f
0
f
.. .
.. .
0
0
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
(34)
(35)
and U , respecwhere i and i are the ith column blocks of −U21 22 tively, as
=[
1 , . . .,
f]
(36)
U22 = [1 , . . ., f ].
(37)
3.3. Estimation of system matrices ˆ can be As we know, the consistent solution of matrices Cˆ and A ˆ f as follows [40,41] extracted directly from ˆ f (1 : ny , :) Cˆ =
(38)
Substituting (45) to (46) yields † ˆ f 1 ) − (Lˆ Lˆ 2 ) = Lˆ 2 vec(H 2
ˆ vec()
× L3 (Lˆ 2 Lˆ 2 )
[1 , . . ., f ]
−U21
†
ˆ f (1 : (f − 1)ny , :) ˆ f (1 + ny : fny , :) ˆ= A
(39)
† where Lˆ 2 = (Lˆ 2 Lˆ 2 ) according to (41).
Hf 1 = L1
L1 =
0
0
f (1 : (f − 1)ny , :)
D
.
B
2 ˆ f 1 − Lˆ 2 vec()) ˆ = min vec(H
(43)
s.t.vec(K) = L3 vec()
vec(Hˆ f 1 ) − Lˆ 2 vec()2 + L3 vec() − vec(K)
(44)
where is a suitably chosen Lagrange multiplier weight. Taking partial derivatives of J with respect to and , and equating them to zero, we have [31]
L3 (Lˆ 2 Lˆ 2 )
−1 L3
× L3 (Lˆ 2 Lˆ 2 ) ˆ = (Lˆ Lˆ 2 ) vec() 2
−1 Lˆ 2 .
(47)
ˆ can be retrieved from ˆ Finally, Bˆ and D
(48)
{yk (t)}N t=1
= D + CB + · · · + CAf B u1 +
f
(49)
CAi wk (t − i) + vk (t)
i=0
= D + CB + · · · + CAf B u2 +
f
(50)
CAi wl (t − i) + vl (t).
Then N {yk (t)}N t=1 − {yl (t)}t=1
= D + · · · + CAf B (u1 − u2 ) + vk (t) − vl (t)
−1
−1 ˆ f 1 ) − vec(K) Lˆ 2 vec(H
−1 −1 ˆ f 1 ) − (Lˆ Lˆ 2 ) Lˆ . Lˆ 2 vec(H 2 3
(45)
(46)
(51)
= K(u1 − u2 ) + vk (t) − vl (t) and N
N {yk (t)}N t=1 − {yl (t)}t=1
t=f
where Lˆ 2 = Inu ⊗ Lˆ 1 and L3 = Inu ⊗ [Iny , C(Inx − A)−1 ]. The operator vec stands for the operation of forming a vector from a matrix by stacking its columns one under another. ⊗ denotes the Kronecker product. By using the Lagrange multipliers method, the problem (43) becomes an optimization problem of the form
=
−1 ˆ f 1 ) − vec(K) Lˆ 2 vec(H
The static gain information is necessary to reduce the error variances of parameter estimates. In this subsection, a method is proposed to estimate the unknown static gain based on two step tests, each having a different input amplitude. Assume that the particular static data for input signals u1 and u2 (u1 = / u2 ) are N denoted by {yk (t)}N and {y (t)} , respectively. Under assumption l t=1 t=1 A1, choose a sufficiently large future horizon f for computation
(42)
Subsequently, the system matrices B and D can be obtained by solving following constrained least squares problem
−1
i=0
K = D + CB + · · · + CA∞ B = [Iny , C(Inx − A)−1 ].
min
−1 L3
3.4. Estimation of static gain
(41)
Assuming that the auxiliary static gain is known as a constant value K, the following constraint can be imposed to compute the system matrices of B and D
ˆ =
L3 (Lˆ 2 Lˆ 2 )
†
{yl (t)}N t=1
, =
ˆ f 1 ). vec(ˆ ls ) = Lˆ 2 vec(H
(40)
Iny
−1 L3
Remark 1. For SIMs without any prior knowledge of the static gain, the system matrices B and D can be estimated by solving a standard least squares problem as follows
where the operator (:) used for partitioning matrices is the same than that in Matlab. It follows from the expression of Hf1 that
where
25
= (N − f + 1)K(u1 − u2 ) +
N
(52) (vk (t) − vl (t)) .
t=f
Under the assumptions A3–A5, the static gain can then be estimated as
N
Kˆ ≈
t=f
N {yk (t)}N t=1 − {yl (t)}t=1
(N − f + 1)(u1 − u2 )
.
(53)
Finally, the complete algorithm of the proposed method presented in the previous sections is summarized below Algorithm 1. 1. 2. 3. 4. 5. 6.
Construct the matrices y˜ k,p , u˜ k,p , y˜ k,f , and u˜ k,f . Estimate static gain parameter K using (53). Eliminate the noise using (22). Estimate f and Hf1 using (32) and (33). Estimate A and C using (38) and (39). Estimate B and D using (47).
26
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
4. Consistency analysis The estimation of the parameters of f and Hf are consistent if the rank condition in (25) is satisfied. Once the estimated f and Hf are consistent, the resulting system matrices are also consistent (see [28,29]). Note that, although the proposed method, and the existing SIMs using parity space like SIMPCA and SIMPCA-Wc, share the same algorithm framework, but they use different instruments. The instruments in the proposed method, SIMPCA and SIMPCA
Theorem 1. Under the assumptions A1–A5 and the following conditions, the rank condition of (25) holds. 1. The input signal u˜ k (t) is persistently exciting of order n + p + f; 2. p ≥ nx and f ≥ nx . Substituting Z˜ k,p and Z˜ k,f into the left side of (25) yields
rank
lim
1 N
N→∞
Z˜ k,f Z˜ k,p
⎛ ⎜ ⎝
rank ⎜
1
lim
f
Hf
0
I
=
x˜ k,f
I
+
N→∞ N
u˜ k,f
0
⎤ ⎞ ⎡ ! y˜ k,p ⎥ ⎟ ⎢ ⎥ ⎟ v˜ k,f ⎢ ⎣ u˜ k,p ⎦ ⎠ . u˜ k,f (54)
Since the input and noise are jointly stationary and uncorrelated with each other, (54) can be rewritten as
rank
lim
N→∞
1 N
⎛ ⎜ 1 rank ⎜ lim ⎝N→∞ N
Z˜ k,f Z˜ k,p
f
Hf
0
I
x˜ k,f u˜ k,f
⎤ ⎞ ⎡ y˜ k,p ⎥ ⎟ ⎢ ⎥ ⎟ ⎢ ⎣ u˜ k,p ⎦ ⎠ .
(55)
u˜ k,f
rank
lim
N→∞
⎛
1 N
Z˜ k,f Z˜ k,p
⎡
x˜ k,p
⎤⎡
f
Hf
0
I
0
l = ⎣0
I
0⎦
0
0
0
l˜xk,p
u˜ k,f
(57)
which is of rank nx + pnu . Then the rank condition in SIMPCA and SIMPCA-Wc is nx rank
lim
N→∞
1 N
Z˜ Z˜ k,p
⎤ ⎞ ⎡
p
0
min(nx + fnu , nx + pnu ).
(58)
It is obvious that the rank condition for SIMPCA and SIMPCA-Wc to be consistent may not be fulfilled, and this explains why the new instrument used in the proposed method enhances the estimation efficiency/accuracy compared to that in SIMPCA and SIMPCA-Wc.
5. Variance analysis The estimation error variance properties of the system matrix are investigated in this section. Compared with the SIM without prior knowledge, for example, the one in (48), the proposed method is able to reduce the estimation error variance of . This result is strictly proved in the following theorem. Theorem 2. Under the assumptions A1–A5, the proposed method results in smaller variances of parameter estimates than SIMs without prior knowledge. That is
E tr vec(˜ ls )vec(˜ ls
)
(59)
˜ = vec() ˆ − vec() and vec(˜ ls ) = vec(ˆ ls ) − vec() where vec() stand for, respectively, the estimation errors of the proposed method and SIMs without prior knowledge. †
Ap
Lp
0
0
0
=
⎜ ⎥⎢ ⎥ ⎟⎢ 1 ⎢ ⎢ ⎥⎢ ⎥ ⎟ rank ⎜ lim ⎝N→∞ N ⎣ u˜ k,p ⎦ ⎣ u˜ k,p ⎦ ⎠ ⎣ 0 u˜ k,f
⎤
Hp
Proof. By substituting Lˆ 2 = L2 + L˜ 2 into Lˆ 2 which is defined in (47), we have
By substituting (6), (7) and (8) into (55), we have
p
˜ ˜ E tr vec()vec( )
=
⎡
−1/2
] , Z = [˜yk,p ] , and Z (ZZ /N) , u˜ u˜ u˜ Wc are [˜yk,p k,p k,f k,p respectively (see (24) in [28], and (28) in [29]). These instruments are all uncorrelated with noise, but only the one in the proposed method, constructed by using both the past and future data, can strictly guarantee the rank condition of (25) under some conditions, see the following theorem.
Proof.
Note that replacing the instrument Z˜ k,p in the proposed method with the instruments in SIMPCA and SIMPCA-Wc, the last term of (56) becomes
I
⎤
Hp
0
I
⎥ . 0⎦
0
I
(56)
The rank of the first two factors in the right-hand side of (56) are equal to nx + fnu , and nx + fnu if p nx and f nx , respectively. The rank of the third factor in the right-hand side of (56) is equal to nx + (p + f)nu if the input signal u˜ k (t) is persistently exciting of order nx + p + f. The rank of the last factor in the right-hand side of (56) is nx + (p + f)nu if p nx and f nx . The details can be found in [39,42–45]. From these results, it is clear that the rank of (56) is equal to nx + fnu , indicating that the rank conditions in (25) hold under the conditions in Theorem 1. This completes the proof. 䊏
† Lˆ 2 =
=
L2 + L˜ 2
L2 + L˜ 2
−1
L2 + L˜ 2
L2 L2 + L2 L˜ 2 + L˜ 2 L2 + L˜ 2 L˜ 2
−1
L2 + L˜ 2
(60) .
As is√well-known in asymptotic analysis, only the terms of order O(1/ N) (in probability) contribute to the asymptotic variance. √ Hence, since L˜ 2 converges to zero, their product is o(1/ N) and can be neglected [46]. With the following relation [47]
L2 L2 + L2 L˜ 2 + L˜ 2 L2 + L˜ 2 L˜ 2
(L2 L2 )
−1
− (L2 L2 )
−1
−1
(L2 L˜ 2 + L˜ 2 L2 )(L2 L2 )
−1
(61)
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
where we have used the symbol to denote that the matrix is approximated by its first order Taylor expansion. Further, (60) can readily be elaborated to † Lˆ 2
(L2 L2 )
−1
− (L2 L2 )
× L2 + L˜ 2 (L2 L2 )
†
−1 L2
− (L2 L2 )
−1
†
= L2 + L2 (L2 ) †
†
−1
+ (L2 L2 )
−1 L˜ 2
†
†
= L2 + L2 (L2 )
†
†
− L2 L˜ 2 L2 − L2 (L2 )
−1 L˜ 2 (I
(62)
−1 L˜ 2 (I
− L2 L2 )
†
†
†
†
−1 L˜ 2 (I
†
†
+ L2 (L2 )
−1 L˜ 2 (I
−1 L3
†
† † L2 (L2 ) Pr
†
(63)
†
†
.
†
†
†
†
E tr L2 (L2 ) Pr Pr
=
†
†
= E tr L2 (L2 )
− L2 L2 ) − L˜ 2 L2 vec(Hf 1 )
†
†
From (65) and (70), we finally have
= vec(H˜ f 1 ) +
−1 (L2 ) L˜ 2 (I
† † − L2 L2 ) − L˜ 2 L2
vec(Hf 1 ).
E
tr vec(˜ ls )vec(˜ ls )
= E tr (L2 L2 )
= E tr
−1 −1 L2 L2 (L2 L2 )
† † L2 (L2 )
(64)
(65)
.
= vec(ˆ ls ) − (Lˆ 2 Lˆ 2 )
−1 L3
L3 (Lˆ 2 Lˆ 2 )
−1 L3
× L3 vec(ˆ ls ) − vec(K)
−1 −1 = vec(ˆ ls ) − (Lˆ 2 Lˆ 2 ) L3 L3 (Lˆ 2 Lˆ 2 ) L3
−1 −1 vec(ˆ ls ) − (L2 L2 ) L3 L3 (L2 L2 ) L3
†
vec() + L2 − (L2 L2 )
% I
−1 L3
−1 − (L2 L2 ) L3
† × L2
= vec() + (I − Pr)
L3 (L2 L2 )
†
†
= −E tr PrL2 (L2 )
−1
(70) .
(71)
.
The inequality of Theorem 1 holds since the term in (71) is positive definite. This completes the proof. 䊏 Remark 2. The asymptotic variances of the estimated transfer functions can be computed by using the estimation errors of the system matrices. The detailed expressions can be found in [48–52], and thus they will not be reproduced in this paper. It is clear that under the same level of estimation errors of the system matrices A and C, the smaller variances of the estimated system matrices B and D, the lower variances of the transfer function can be observed. 6. Illustration
−1 (66)
× L3 vec(˜ ls )
vec() +
−1
× L3 (vec() + vec(˜ ls )) − vec(K)
† × L3 L2
−1
(69)
− E tr(vec(˜ ls )vec(˜ ls ) )
†
†
= −E tr Pr(L2 L2 )
= −E tr PrL2 (L2 )
˜ can be derived In a similar way, by substituting (62) into (47), vec() as ˆ vec()
˜ ˜ ) E tr(vec()vec( )
The variance of vec(˜ ls ) in the SIM without prior knowledge can then be computed as follows
− E tr PrL2 (L2 )
†
which is equal to the third term on the right side of (68), and thus
= vec() + L2 where
(68)
E tr L2 (L2 ) Pr
˜ ˜ E tr vec()vec( )
†
=
The fourth term on the right side of (68) can be written as E tr PrL2 (L2 ) Pr
†
†
†
†
− E tr PrL2 (L2 )
+ E tr PrL2 (L2 ) Pr
†
(67)
†
†
− E tr
− L2 L2 ) − L˜ 2 L2
†
−1 L3
†
† † ˜ f 1) L2 vec(Hf 1 ) + L2 vec(H
˜ ˜ tr vec()vec( )
˜ f 1) × vec(Hf 1 ) + vec(H
L3 (L2 L2 )
− L2 L2 ) − L˜ 2 L2 .
L2 + L2 (L2 )
= E tr L2 (L2 )
†
−1 L3
= E tr (I − Pr)L2 (L2 ) (I − Pr)
−1 † L˜ 2 L2 L2
ˆ f 1) = Lˆ 2 vec(H
Pr = (L2 L2 )
E
ˆ f 1 = Hf 1 + H ˜ f 1 into (48) and neglecting the Substituting (62) and H ˜ f 1 and L˜ 2 immediately give product of H vec(˜ ls )
where
is a projection matrix satisfying PrPr = Pr. With the above results ˜ in the proposed method can be in mind, the variance of vec() computed as follows
−1 (L2 L˜ 2 + L˜ 2 L2 )(L2 L2 ) L2
†
†
−1
−1 L˜ 2
= L2 − L2 L˜ 2 L2 + L2 (L2 ) †
(L2 L˜ 2 + L˜ 2 L2 )(L2 L2 )
27
−1 L3
−1 L3 (L2 L2 ) L3
−1
−1 & L3
In this section, two examples are used to demonstrate the effectiveness and merit of the proposed method. One is a benchmark example from [53] and the other is an industrial injection molding process studied in [54]. Example 1. Consider the following multiple-input, multipleoutput system which is the nominal model adopted for numerical simulation study in [53]:
xk (t + 1) = Axk (t) + Buk (t) + wk (t) yk (t) = Cxk (t) + Duk (t)u + vk (t)
(72)
28
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
Fig. 1. NMAE of estimated Bode magnitude plots using different methods for example 1.
Fig. 2. Variance of estimated Bode magnitude plots using different methods for example 1.
with
⎡
0.8
A = ⎣0
C=
0 0.5 0
−0.4
0.2
⎤
⎡
0
−0.5 ⎦ , B = ⎣ 0
0.3 0 0.5
0.5
0
0
1
,D =
0
without utilizing prior information [see (52)], and SIMPCA-Wc with
⎤
u˜ u˜ prior information (replacing y˜ k,p k,p k,f
−0.6 ⎦ ,
0.5 0 0.7139 0.1174
the same figures. From these results, we can conclude that
0.3131 0.2876
1 ˆ |Gij (ω) − Gij (ω)| 100
by Z ) are shown in
(73)
where the input signal uk (t) is chosen as a zero-mean, Gaussian, random sequence with unit variance. The noise vk (t) is assumed to be zero-mean, Gaussian, white noise of variance 0.5. The disturbance wk (t) is a sinusoid signal wk (t) = sin(2 /3). The statistical properties of the proposed method are investigated by means of a Monte-Carlo (MC) simulation of 100 runs. In each run, two sets of batch dynamical data are generated to obtain the disturbance-free data, each having a different input sequence. The static data are generated by using two constant inputs: u1 = 0.3 and u2 = 0.95. The performance index for the identification methods, i.e. the normalized mean absolute error (NMAE), is computed based on the frequency response data of the 100 estimated models. The NMAE is defined as 100
NMAE =
• Compared with the conventional SIMPCA-Wc where only past input/output data are used to build the instrument, the proposed method and the SIM without utilizing prior information generally enhance the estimated model efficiency/accuracy, thanks to the instrument which is built based on both the past input/output and future input data, see Section 4 for the reason on this improvement. • While the proposed method and the SIM without utilizing prior information give consistent estimated results, the proposed method reduces the variance of the estimated models, thanks to the prior knowledge of static gain which contains adequate information on the low-frequency band of the system. A detailed explanation can be found in Section 5.
Example 2. in [54]:
Consider an injection molding batch process studied
⎧ ⎨
(74)
i=1
where ω ∈ [0, ] is the frequency variable, and Gij (ω) is the magnitude from jth input to ith output. Letting N = 500 and f = p = 10, the NMAE and variance are compared in Figs. 1 and 2. For ease of comparison, the results of the SIM
xk (t + 1) =
1.582
−0.592
xk (t) + 1 0 ⎩ yk (t) = 1.69 1.419 xk (t) + vk (t)
1 0
uk (t) + wk (t)
(75)
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
Fig. 3. NMAE of estimated Bode magnitude plots using different methods for Example 2.
29
Subsequently, the system matrices A and C are extracted from the estimated extended controllability matrix, while the system matrices B and D are extracted from the estimated triangular Toeplitz matrix. The proposed method can be regarded as an improved SIM in relative to the existing SIM using parity space, e.g. SIMPCA and SIMPCA-Wc. Due to the instrument used in the proposed method, built based on the past and future input/output to eliminate the impact of measurement noise, is better correlated with system matrix than the ones used in SIMPCA and SIMPCA-Wc, the estimation efficiency/accuracy can be greatly enhanced. This property has been analyzed with a strict proof. As the second contribution of this paper, a technique has been proposed to retrieve the prior knowledge of the system static-gain information using the auxiliary static data, which helps to improves the quality and variance properties of the estimated B and D. Two illustrative examples have shown the effectiveness and merit of the proposed techniques and identification method. Acknowledgments
Fig. 4. Variance of estimated Bode magnitude plots using different methods for Example 2.
where the noise vk (t) is zero-mean, Gaussian and white, with variance of 0.5. The input is designed as a random binary sequence in Matlab as follows u = 0.5 + 0.3 ∗ idinput(Ts, rbs ) ; The disturbance wk (t) that accounts for the mold cavity pressure will affect the injection velocity. Just as an example, it is assumed that wk (t) has the expression wk (t) =
−a1 z −1 − a−2 2 1 − b1 z −1
*
where b1 =
*
0.4
0.3 t < 5000 0.4 t ≥ 5000
(76)
* , a1 =
0.05 t < 5000 0.15 t ≥ 5000
, and
0.08 t < 5000 . 0.15 t ≥ 5000 The statistical properties of the proposed method are illustrated by means of an MC simulation of 100 runs. The simulation settings remain the same as in the previous example, except for u1 = 0.8 and u2 = 0.95. With N = 2000 and f = p = 10, the NMAE and variance of the estimated models using the same identification methods presented in Example 1 are shown in Figs. 3 and 4. Again, the proposed method results in very accurate estimates with smaller variances than the other two identification methods, especially in the presence of repeatable disturbances. a2 =
7. Conclusion For identification of batch processes in the presence of repeatable disturbances as well as output measurement noises, a prior-knowledge-based subspace identification method has been proposed to achieve consistent estimation of the system parameters. As the first contribution of this paper, a new operator has been proposed to eliminate the influence of repeatable disturbances, with the underlying rationale being to collect one more data set to recover the knowledge of the disturbances. Then, by applying the proposed difference operator to eliminate the disturbances in the data, consistent estimates of the extended controllability and triangular Toeplitz matrices can be obtained using the proposed SVD-based method from the disturbance-free batch data.
This work was supported by the National Natural Science Foundation of China under grants 61903057, 61803061, 61703347, 61703066 and 61703311; the Chongqing Natural Science Foundation under grants cstc2019jcyj-msxmX0129, cstc2016jcyjA0428 and cstc2018jcyjAX0536; the Science and Technology Research Program of Chongqing Municipal Education Commission under grant KJQN201800603. References [1] D. Bonvin, B. Srinivasan, D. Hunkeler, Control and optimization of batch processes, IEEE Control Syst. Mag. 26 (6) (2006) 34–45. [2] S. Hao, T. Liu, W. Paszke, K. Galkowski, Robust iterative learning control for batch processes with input delay subject to time-varying uncertainties, IET Control Theory Appl. 10 (15) (2016) 1904–1915. [3] S. Hao, T. Liu, F. Gao, PI based indirect-type iterative learning control for batch processes with time-varying uncertainties: a 2D FM model based approach, J. Process Control 78 (2019) 57–67. [4] D. Li, S. He, Y. Xi, T. Liu, F. Gao, Y. Wang, J. Lu, The synthesis of ILC-MPC controller with data-driven approach for constrained batch processes, IEEE Trans. Ind. Electron. (2019), http://dx.doi.org/10.1109/TIE.2019.2910034. [5] Z. Cao, J. Lu, R. Zhang, F. Gao, Iterative learning Kalman filter for repetitive processes, J. Process Control 46 (2016) 92–104. [6] Z. Cao, R. Zhang, Y. Yang, J. Lu, F. Gao, Discrete-time robust iterative learning Kalman filtering for repetitive processes, IEEE Trans. Autom. Control 61 (1) (2016) 270–275. [7] C. Zhao, F. Wang, M. Jia, Dissimilarity analysis based batch process monitoring using moving windows, AIChE J. 53 (5) (2007) 1267–1277. [8] J. Liu, T. Liu, J. Chen, Window-based stepwise sequential phase partition for nonlinear batch process monitoring, Ind. Eng. Chem. Res. 55 (34) (2016) 9229–9243. [9] J. Liu, T. Liu, J. Chen, Sequential local-based Gaussian mixture model for monitoring multiphase batch processes, Chem. Eng. Sci. 181 (2018) 101–113. [10] D. Shi, N.H. El-Farra, M. Li, P. Mhaskar, P.D. Christofides, Predictive control of particle size distribution in particulate processes, Chem. Eng. Sci. 61 (1) (2006) 268–281. [11] M.M. Rashid, P. Mhaskar, C.L.E. Swartz, Handling multi-rate and missing data in variable duration economic model predictive control of batch processes, Chem. Eng. Sci. 63 (7) (2017) 2705–2718. [12] A.Y.D. Tsen, D.S.H. Jang, S.S. Wong, B. Joseph, Predictive control of quality in batch polymerization using hybrid ANN models, AIChE J. 42 (2) (1996) 455–465. [13] Z. Cao, J. Lu, R. Zhang, F. Gao, Online average-based system modelling method for batch process, Comput. Chem. Eng. 108 (4) (2018) 128–138. [14] Z. Cao, Y. Yang, H. Yi, F. Gao, Priori knowledge-based online batch-to-batch identification in a closed loop and an application to injection molding, Ind. Eng. Chem. Res. 55 (32) (2016) 8818–8829. [15] Z. Cao, Y. Yang, J. Lu, F. Gao, Constrained two dimensional recursive least squares model identification for batch processes, J. Process Control 24 (6) (2014) 871–879. [16] A.W. Dorsey, J.H. Lee, Building inferential prediction models of batch processes using subspace identification, J. Process Control 13 (5) (2003) 397–406. [17] Y. Yao, F. Gao, Subspace identification for two-dimensional dynamic batch process statistical monitoring, Chem. Eng. Sci. 63 (13) (2008) 3411–3418. [18] B. Corbett, P. Mhaskar, Subspace identification for data-driven modeling and quality control of batch processes, AIChE J. 62 (5) (2016) 1581–1601.
30
J. Hou et al. / Journal of Process Control 82 (2019) 22–30
[19] A. Garg, P. Mhaskar, Subspace identification based modeling and control of batch particulate processes, Ind. Eng. Chem. Res. 56 (26) (2017) 7491–7502. [20] S.J. Qin, An overview of subspace identification, Comput. Chem. Eng. 30 (10-12) (2006) 1502–1513. [21] G. van der Veen, J.W. van Wingerden, M. Bergamasco, M. Lovera, M. Verhaegen, Closed-loop subspace identification methods: an overview, IET Control Theory Appl. 7 (10) (2013) 1339–1358. [22] T. Katayama, Subspace Methods for System Identification, Springer Science & Business Media, Berlin, 2006. [23] S.T. Navalkar, J.W. van Wingerden, E. van Solingen, T. Oomen, E. Pasterkamp, G.A.M. Van Kuik, Subspace predictive repetitive control to mitigate periodic loads on large scale wind turbines, Mechatronics 24 (8) (2014) 916–925. [24] I. Houtzager, J.W. van Wingerden, M. Verhaegen, Rejection of periodic wind disturbances on a smart rotor test section using lifted repetitive control, IEEE Trans. Control Syst. Technol. 21 (2) (2013) 347–359. [25] I. Houtzager, J.W. van Wingerden, M. Verhaegen, Wind turbine load reduction by rejecting the periodic load disturbances, Wind Energy 16 (2) (2013) 235–256. [26] P.M. Gebraad, J.W. van Wingerden, P.A. Fleming, A.D. Wright, LPV identification of wind turbine rotor vibrational dynamics using periodic disturbance basis functions, IEEE Trans. Control Syst. Technol. 21 (4) (2013) 1183–1190. [27] N.E. Goodzeit, M.Q. Phan, System identification in the presence of completely unknown periodic disturbances, J. Guid. Control Dyn. 23 (2) (2000) 251–259. [28] J. Wang, S.J. Qin, A new subspace identification approach based on principal component analysis, J. Process Control 12 (8) (2002) 841–855. [29] J. Wang, S.J. Qin, Closed-loop subspace identification using the parity space, Automatica 42 (2) (2006) 315–320. [30] P. Trnka, V. Havlena, Subspace like identification incorporating prior information, Automatica 45 (4) (2009) 1086–1091. [31] A. Alenany, H. Shang, M. Soliman, I. Ziedan, Improved subspace identification with prior information using constrained least squares, IET Control Theory Appl. 5 (13) (2001) 1568–1576. [32] A. Alenany, H. Shang, Recursive subspace identification with prior information using the constrained least squares approach, Comput. Chem. Eng. 54 (11) (2013) 174–180. [33] I. Markovsky, Mercère, Subspace identification with constraints on the impulse response, Int. J. Control 90 (8) (2017) 1728–1735. [34] Y. Wang, L. Zhang, Y. Zhao, Improved closed-loop subspace identification with prior information, Int. J. Syst. Sci. 49 (9) (2018) 1821–1835. ˇ ˇ L. Ferkl, M. Sebek, [35] S. Privara, J. Cigler, Z. Vána, Subspace identification of poorly excited industrial systems, 49th IEEE Conference on Decision and Control (CDC) (2010) 4405–4410. [36] K. Li, H. Luo, C. Yang, S. Yin, Subspace-aided closed-loop system identification with application to DC motor system, IEEE Trans. Ind. Electron. (2019), http:// dx.doi.org/10.1109/TIE. 2019.2907447.
[37] P. van Overschee, B.L. De Moor, Subspace Identification for Linear Systems: Theory-Implementation-Applications, Springer Science & Business Media, 2012. [38] M. Jansson, B. Wahlberg, On consistency of subspace methods for system identification, Automatica 34 (12) (1998) 1507–1519. [39] J. Hou, F. Chen, P. Li, Z. Zhu, F. Liu, An improved consistent subspace identification method using parity space for state-space models, Int. J. Control Autom. Syst. 17 (5) (2019) 1167–1176. [40] C. Yu, M. Verhaegen, Subspace identification of distributed clusters of homogeneous systems, IEEE Trans. Autom. Control 62 (1) (2016) 463–468. [41] M. Inoue, Subspace identification with moment matching, Automatica 99 (2019) 22–32. [42] J. Hou, T. Liu, B. Wahlberg, M. Jansson, Subspace Hammerstein model identification under periodic disturbance, IFAC-PapersOnLine 51 (15) (2018) 335–340. [43] J. Hou, T. Liu, Closed-loop subspace model identification using innovation estimation and orthogonal projection, Acta Autom. Sin. 42 (11) (2016) 1657–1663. [44] J. Hou, T. Liu, Q.G. Wang, Subspace identification of Hammerstein-type nonlinear systems subject to unknown periodic disturbance, Int. J. Control (2019), http://dx.doi.org/10.1080/00207179.2019.1621385. [45] C. Yu, L. Ljung, M. Verhaegen, Identification of structured state-space models, Automatica 90 (2018) 54–61. [46] A. Chiuso, G. Picci, Asymptotic variance of subspace methods by data orthogonalization and model decoupling: a comparative analysis, Automatica 40 (10) (2004) 1705–1717. [47] M. D“ohler, L. Mevel, Efficient multi-order uncertainty computation for stochastic subspace identification, Mech. Syst. Signal Process. 38 (2) (2013) 346–366. [48] M. Jansson, Asymptotic variance analysis of subspace identification methods, IFAC Proc. Vol. 33 (15) (2000) 91–96. [49] A. Chiuso, G. Picci, The asymptotic variance of subspace estimates, J. Econom. 118 (1-2) (2004) 257–291. [50] J. Hou, T. Liu, Q.G. Wang, Recursive subspace identification subject to relatively slow time-varying load disturbance, Int. J. Control 91 (3) (2018) 622–638. [51] J. Hou, F. Chen, P. Li, Z. Zhu, Fixed point iteration-based subspace identification of Hammerstein state-space models, IET Control Theory Appl. 13 (8) (2019) 1173–1181. [52] J. Hou, T. Liu, F. Chen, Orthogonal projection based subspace identification against colored noise, Control Theory Technol. 15 (1) (2017) 69–77. [53] G. Mercère, L. Bako, S. Lecœuche, Propagator-based methods for recursive subspace model identification, Signal Process. 88 (3) (2008) 468–491. [54] T. Liu, B. Huang, S.J. Qin, Bias-eliminated subspace model identification under time-varying deterministic type load disturbance, J. Process Control 25 (2015) 41–49.