4292
Ind. Eng. Chem. Res. 2001, 40, 4292-4301
Analysis and Design of a Linear Input/Output Data-Based Predictive Control In-Hyoup Song, Kee-Youn Yoo, and Hyun-Ku Rhee* School of Chemical Engineering and Institute of Chemical Processes, Seoul National University, Kwanak-ku, Seoul 151-742, South Korea
In this work, a subspace identification algorithm is reformulated from a control point of view. The proposed algorithm is referred to as an input/output data-based predictive control, in which an explicit model of the system to be controlled is not calculated at any point in the algorithm. First, the state estimation obtained by the subspace identification algorithm is analyzed in comparison with the receding-horizon-based estimation. For the subspace algorithm, it is wellknown that a Kalman filter state is calculated by simple linear algebra under specific conditions. In general, however, it is shown that the present state estimation scheme gives a biased state estimate and that it has a structure similar to the best linear unbiased estimation (BLUE) filter obtained by solving the least-squares problem analytically. With such an interpretation of the state estimation, we augment the integrated white noise model to add integral action to a linear input/output data-based predictive controller and use each the BLUE filter and the Kalman filter as a stochastic observer for the unmeasured disturbance. The proposed linear input/ouput data-based predictive controller is applied to the property control of a continuous styrene polymerization reactor to demonstrate its improved performance. 1. Introduction A number of modern control methods, such as model predictive control, use an artificial model to design a control system.1 A model predictive controller (MPC), highlighted as an advanced controller in recent years, can demonstrate outstanding control performance when the controller is based on a good model and a reasonably accurate state estimate. When such a model is not available, the design of the model predictive control requires the first step of system identification from input/output data. Once a specific model is available, the state of the system must be estimated using the second step of a stochastic observer design such as a Kalman filter. Finally, the third step consists of designing an on-line optimization for the model predictive controller. Most researchers concentrate either on finding a controller given the plant model or on identifying a model from the data. There is no doubt that model-based control and system identification are closely related, simply because one strongly depends on the other. The system identification step is, in fact, nothing more than a vehicle for controller design. Thus, the system identification step could be reformulated by devising a method that allows for the calculation of a controller directly from the input/output data. The idea of computing a linear quadratic Gaussian (LQG) controller directly from input/output data was first introduced to a technique called iterative feedback tuning (IFT).2,3 In the context of subspace identification, a number of researchers have used the subspace identification method to obtain the state space model for the design of a model predictive controller.4,5 These studies followed the conventional three steps to designing a * To whom correspondence should be addressed. Tel.:+822-880-7405. Fax:+82-2-888-7295. E-mail:
[email protected].
controller. Amirthalingam et al.6 also proposed a twostep procedure for data-based modeling for inferential predictive control system design. These authors noticed that the numerical algorithms for the subspace state space system identification (N4SID) method gave biased results in cases when the input data were correlated with the residual and proposed a method for the selection of an appropriate data set from a huge historical database. Favoreel et al.7,8 suggested a model-free subspace-based LQG controller that allows for the construction of an optimal controller directly from measurements of the input/output data, and Favoreel and De Moor also developed a subspace predictive control.9 In subspace-based controller design, the state estimation in subspace identification is obtained by a receding(moving-) horizon data-based type of method. The state estimate in the subspace identification is known to be equal to a Kalman filter (KF) state under specific conditions.10 However, the specific conditions, such as infinite number of available data points and adequately chosen initial state, are not usually satisfied in real plants. Therefore, we focus our attention on the receding-horizon type of estimator to analyze the state estimate in the subspace identification. MacGregor and co-workers11,12 introduced statistical monitoring using on-line process measurements for batch and semibatch processes. Some researchers extended MacGregor’s approach to recursive data-based prediction in a batch product control problem. Robertson et al.13 developed a moving-horizon-based approach for least-squares estimation and compared the performance of the proposed estimator with that of other stochastic observers such as KF, extended Kalman filter (EKF), and optimizationbased estimators. Also, Russel and co-workers14 proposed another data-based predictor by applying the orthogonal projection method. In these studies, the estimator is constructed to calculate the control input,
10.1021/ie0010159 CCC: $20.00 © 2001 American Chemical Society Published on Web 06/13/2001
( )
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4293
rather than calculating the control input directly from the input/output data. Ruscio and Foss8 defined a leastsquares estimate of the current state in a formanalogous to that of the best linear unbiased estimation (BLUE) filter, but they did not perform a comparison with the state estimation in the subspace identification. Here, we use a subspace technique to design an input/ output data-based predictive controller and explore the kind of state estimate that is obtained. It is shown that the state estimation scheme in subspace identification results in a biased state estimate, even if the input data are not correlated with the residual, and that the state estimation scheme has a structure similar to that of the BLUE filter, which is known to be the most efficient within the class of linear filters. With this estimation scheme, we add integral action to a linear input/output data-based predictive control. There are a number of methods for introducing an integrator into a control law. For subspace-based predictive control, the integral property was added to the control law using an objective function containing the velocity form of the control input.8,15,16 In this study, we use a different approach, i.e., we introduce the integrated white-noise disturbance and use each the BLUE filter and the Kalman filter as a stochastic observer for the unmeasured disturbance. The control performance is corroborated by treating the property control of a continuous styrene solution polymerization reactor.
2. An Input/Output Data-Based Prediction Model: An Overview We are interested in the design of model predictive controller for linear time-invariant systems that can be described by the following state space equations
xk+1 ) Axk + Buk + ωk yk ) Cxk + Duk + νk
(1)
where xk ∈ Rn is the state variable, uk ∈ Rm is the input vector, and yk ∈ Rl is the output vector. Favoreel et al.7,8 developed an input/output databased prediction model, the formulation of which is summarized in this section. For more detail, it is recommended that the reader consult their articles.7,8 According to these authors, the subspace identification problem can be interpreted as follows: Given the past inputs and outputs Wp ) (YpT UpT)T and the future inputs Uf, find an optimal prediction of the future outputs Y ˆ f. Here, a matrix input/output equation is used, i.e.
Yf ) ΓNpXf + HNpUf + HNpsZf + Vf ≈ LwWp + LuUf (2) The system-related matrices included in eq 2 are defined as follows
C CA 2 ΓNp ) CA , HNp ) · · · Np-1 CA
(
)
D CB D · , ·· ·· · · · · Np-2 CA B CANp-3 · · · D
(
0 ·· C s ) HN · ··· ·· p · · · · Np-2 · · · C 0 CA
)
and
(
Xf ) (xM xM+1 · · · xM+Np-1 ), u0 u1 · · · u1 u2 · · · Up ) · ·· · · ·u M-1 uM · · · uM uM+1 uM+1 uM+2 Uf ) · · ·u u
(
M+Np-1
M+Np
uj-1 uj
)
,
uM+j-2 · · · uM+j-1 · · · uM+j ·· · · · · uM+Np+j-2
)
(3)
where M and Np denote the backward receding horizon and the prediction horizon, respectively. In a similar way, Yp and Yf are defined for the outputs yk, and Zf and Vf are defined for the process disturbances ωk and the measurement noises νk, respectively. The leastsquares prediction Y ˆ f of Yf can be found from the following least-squares problem:
||
( )||
W minLw,Lu Yf - (Lw Lu ) U p f
2
(4)
F
The solution to this problem is the orthogonal projection of the row space of Yf on the row space spanned by Wp and Uf, i.e.
( )
W Y ˆ f ) Yf/ U p ) Yf/UfWp + Yf/WpUf ) LwWp + LuUf (5) f where the input signal is persistently exciting of order 2Np and is not related to the residual of the process. The first term on the right-hand side of eq 5 is interpreted as follows:
Yf/UfWp ) ΓNpXf/UfWp + HNpUf/UfWp + HNpsZf/UfWp + Vf/UfWp
(6)
) ΓNpX ˆf The last two projections on the right-hand side are equal to zero under the assumption that the deterministic input uk is persistently exciting and uncorrelated with both the process disturbance and the measurement noise. In the simulation study, we use a pseudo-random four-level signal. The numerical implementation of the projection in eq 5 can be done in an efficient way by making the following QR decomposition:
4294
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001
)(
()(
Q1T Wp R11 Uf ) R21 R22 Q2T Yf R31 R32 R33 Q T 3
)
(7)
By posing
(
R L ) (R31 R32 ) R11 R 21 22
)
3. Analysis of State Estimate Obtained from Subspace Identification
†
(8)
where the superscript indicates the pseudo-inverse, it is straightforward to show that eq 5 can be written as
( ) ( )
W W Yf / U p ) L U p f f
(9)
Using Matlab notation, we then have
Lw ) L(:, 1:M(m + l)) Lu ) L(:, M(m + l) + 1:end)
(10)
These two matrices in the prediction model form the backbone of the prediction into the future. These matrices are obtained in an open-loop manner. If they can be obtained in a closed loop, then one can obtain adaptive Lw and Lu matrices and, consequently, implement adaptive I/O data-based predictive controller. From eqs 5 and 9, we notice that the optimal predicted output sequence yˆ f can be expressed as
yˆ f ) Lwwp + Luuf
one another by the specific choice of the weighting matrices W1 and W2, but such algorithms have the common basis that the Kalman filter state can be obtained.
In subspace identification theory, it is well-known that the stochastic system can be transformed into the forward innovation form if the state covariance matrix is chosen to be equal to P, the solution of the algebraic Riccatti equation. In that case, when the number of available data points is infinite, the product between the extended observability matrix ΓNp and the nonsteady-state Kalman filter estimate X ˆ f is equivalent to the oblique projection of the row space of Yf along the row space of Uf on the row space of Wp.10 Here, we are interested in the state estimate obtained by the subspace identification algorithm in the case when the state covariance matrix is not equal to P, the number of available data points is finite, and the initial state is chosen arbitrarily. First, the state estimation obtained by the subspace identification of a purely stochastic system with no external input (uk ≡ 0) will be analyzed. Such a system can be analyzed rather easily by virtue of its simple structure. Consider the following stochastic system
xk+1 ) Axk + wk
where xk ∈ Rn is the state variable and yk ∈ Rl is the output vector. wk and vk represent zero-mean, white noise sequences with covariances Rw and Rv, respectively. Without the external input, eq 2 is reduced to the form
where
( )
y wp ) up p
with M last known values of the inputs up ∈ R and the outputs yp ∈ R, which are defined as
( ) ( )
yt-M+1 · yp ) ·· y
ut-M+1 · and up ) ·· u
t-1
(12)
t-1
yt
ut
The extended observability matrix ΓNp and the state estimate X ˆ f can be found from the weighted singular value decomposition of Lw, i.e.
(
W1LwW2 ) (U1 U2 )
S1 S2
)(
)
V1T ≈ U1S1V1T V2T
(13)
from which we obtain -1
(15)
yk ) Cxk + vk
(11)
1/2
ΓNp ) W1 U1S1
X ˆ f ) S11/2 V1TW2-1Wp
(14)
Note that the above singular value decomposition is weighted by two weighting matrices W1 and W2 and that an appropriate choice for W1 and W2 leads to a balanced controller.10 In this study, we focus on the N4SID method to analyze the property of the state estimate. Various subspace algorithms, such as N4SID, multivariable output-error state space (MOESP), and canonical variate analysis (CVA), can be distinguished from
Yf ) ΓNpXf + HNpsZf + Vf
(16)
By using different system-related matrices, the present state vector can be expressed as a linear combination of the past input/output data vectors. Consider the following matrix input/output equation
h MX f + H h MU p + H h MsZp + Vp Yp ) Γ
( ) (
(17)
where
CA-M CA-M+1 , Γ hM ) · · · -1 CA CA-1B CA-2B · · · CA-1B · · · H hM) ·· ·
(
CA-1
H h Ms )
)
CA-MB CA-M+1B , · · · -1 CA B CA-2 · · · CA-M CA-1 · · · CA-M+1 ·· · · ·· CA-1
)
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4295
The matrix input/output equation corresponding to the purely stochastic system is given by eliminating the second term on right-hand side of eq 17, i.e.
Yp ) Γ h MX f + H h MsZp + Vp
(18)
Van Overschee and De Moor10 showed that the orthogonal projection of the row space of Yf onto the row space of Yp is equal to the product of the extended observability matrix and the forward Kalman filter state sequence for a purely stochastic system when both the initial state estimate and the initial covariance of the state estimate are equal to zero. Now we are interested in a more general system that one might encounter in real chemical processes. Let us calculate the state estimate for a purely stochastic process using the projection method. The first step is to calculate the orthogonal projection Θ of the row space of Yf onto the row space of Yp, i.e. -1
Θ ≡ Yf/Yp ) CMLM Yp
(19)
where CM and LM are defined as
CM ) E{YfYp } T
(20)
LM ) E{YpYpT}
Substituting eqs 16 and 18 into eq 20, we obtain CM and LM in terms of the state covariance and noise covariance matrices as
CM ) E{(ΓNpXf + HNp Zf + Vf)(Γ h MXf + H h M Zp + Vp) } s
s
) ΓNpE{XfXfT}Γ h MT ) ΓNpPsΓ h MT
T
(21)
LM ) E{(Γ h MX f + H h MsZp + Vp)(Γ h MXf + H h MsZp + Vp)T} h MT + H h Ms diag(Rw, ..., Rw)H h MsT + )Γ h ME{XfXfT}Γ diag(Rv, ..., Rv) (22) h MT + ΞM )Γ h MPsΓ h Ms diag(Rw, ..., Rw)H h MsT + diagwhere ΞM represents H (Rv, .., Rv) and Ps is the block matrix consisting of the state covariance matrix. It is assumed that the process noise, the measurement noise, and the state are not correlated with one another. Substituting these equations into eq 19, we obtain
Θ ) CMLM-1Yp ) ΓNpPsΓ h MT(Γ h MPsΓ h MT + ΞM)-1Yp (23) ˆ f, According to Van Overschee and De Moor,10 Θ ) ΓNpX and hence
h MT(Γ h MP sΓ h MT + ΞM)-1Yp ) HsYp X ˆ f ) PsΓ
The linear receding-horizon filter can be expressed as a linear combination of past input/output data
X ˆ f ) LyYp + LpUp
(26)
If X ˆ f is to be unbiased, then the following equation must hold
E{X ˆ f} ) E{LyYp + LpUp} ) E{Ly(Γ h MXf + LpH h MU p + h MsZp + LyVp) + LpUp} (27) LyH h ME{Xf} + (LyH h M + Lp)E{Up} ) Ly Γ ) E{Xf} For every set of Up values, eq 27 is satisfied if and only if
h M ) IM LyΓ h M + Lp ) 0 L yH
(28)
For purely stochastic systems without external input, hM one can disregard the second condition, and thus LyΓ ) IM is the only condition that must be satisfied to obtain an unbiased estimate. h M can be expressed in the form The product of Hs and Γ
Hs Γ h M ) (Γ h MT ΞM-1Γ h M + Ps-1)-1Γ h MT ΞM-1Γ h M (29) The right-hand side of eq 29 is not equal to the identity matrix, which indicates that the state estimation scheme in the subspace identification gives a biased estimate. Therefore, a biased estimate is obtained even if the input data are not correlated with the residual. It is to be noted that the system under consideration is purely stochastic and that no assumptions are made other than that the measurement noise and process noise are not correlated with each other. The right-hand side of eq 29 approaches the identity matrix if the norm of the state covariance matrix is very large or the norm of the noise covariance matrix is very small. Thus, when the system is corrupted by white noise with a small amplitude or when the amplitude of state is much larger than that of the noise, a better state estimate can be obtained. Here, we have not prescribed any initial conditions such as the initial state, the initial state covariance, and so on. In the case of a combined system consisting of a stochastic part and a deterministic part having external input, it can easily be shown that the state estimation scheme in the subspace identification gives a biased estimate if certain conditions are not satisfied. By applying the oblique projection of Yp along the row space of Up on the row space of Yp, we obtain
(24)
Because the terms in parentheses in eq 24 are diagonal, one can proceed as follows:
h MT(Γ h MP sΓ h MT + ΞM)† ) {Γ h M + ΞM(Γ h MT)†Ps†}† Hs ) PsΓ ) (Γ h MT ΞM†Γ h M + Ps†)†Γ h MT ΞM† )
(25)
h M + Ps-1)-1Γ h MT ΞM-1 (Γ h MT ΞM-1Γ
Note that the past output data are correlated with the past noise. If we assume that Γ h MLy is equal to the identity matrix, then eq 30 yields
4296
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001
Yp ) Yp + H h MsZp/UpYp + Vp/UpYp
(31)
This contradictory result indicates that Γ h MLy is not equal to the identity matrix. Consequently, the state estimation scheme in subspace identification always gives a biased estimate. It appears to be natural because, in the combined system, a state can be considered as the sum of a stochastic part and a deterministic part. Even if the deterministic part can be estimated by its true value in the subspace estimation scheme, there still exists a bias with respect to the stochastic part. Finally, we compare the structure of the gain matrix Hs to that of BLUE filter. The BLUE filter uses the finite impulse response (FIR) structure matrix that transforms the system-related matrix to estimate the future state using the past input/output data with the minimum error variance constraint.17 The BLUE filter also gives the state estimate by processing the input/ output data linearly. The gain matrix HB is designed in such a way that the minimum error variance constraint and the unbiasedness constraint are satisfied, i.e.
ˆ fB)(Xf - X ˆ fB)T] HB ) argminHB E[(Xf - X
(32)
H BΓ h M ) IM
(33)
The structure of the BLUE filter can be expressed as
h MUp) X ˆ fB ) HB(Yp - H
(34)
process. We assume that dk is generated through the stochastic difference equation
xk+1w ) Awxkw + Bwω jk dk ) Cwxkw
where ω j k is the discrete-time white noise with covariance Rwj . It is also possible that the measurements of yk are corrupted by measurement noise νjk that is white noise with covariance of Rvj. These integral white noise terms are augmented to the prediction model as follows:
( ) (
)( ) ( ) ( )
xk+1 A BdCw xk 0 B ) uk + B ω jk w + A 0 xk+1w x 0 w w k yˆ k ) (C CdCw )
( )
xk + Duk + νjk xkw
˜ NpΓNpw ) Γ h Np ) (ΓNp H (35)
ΞM ) H h Ms diag(Rw, ..., Rw)H h MsT + diag(Rv, ..., Rv) Comparing the gain matrix HB in the BLUE filter with that of the subspace identification scheme Hs of eq 25, we find that their structures are analogous to each other. If the state covariance term can be neglected in Hs, then we obtain exactly the same gain filter by the subspace identification algorithm as the gain matrix of the BLUE filter. This implies that the state estimate in the subspace identification is not the simple pseudoinverse of Γ h M but close to that of BLUE filter obtained by solving the least-squares objective function analytically. Therefore, one can conclude that the state estimation scheme in the subspace identification yields a biased estimate regardless of the characteristics of the input data, but in cases when the ratio of the norm of the noise error covariance to the norm of the state covariance is very small, the gain matrix of the subspace identification algorithm becomes approximately equal to that of the BLUE filter. 4. Addition of Integrated White Noise Disturbance Model to the Prediction Model To the basic structure of system we add the capability of estimating the effect of disturbances on future outputs. For the linear disturbance model, we will use the “type 1” model that has been used in model predictive control design.18 It is common to express the unmeasured disturbance dk in terms of a stochastic
(37)
This type of disturbance description is able to capture a constant bias between the plant and the prediction model. Therefore, if the combined effects of the external disturbance and the prediction model mismatch reach a constant value, the type 1 assumption, along with a controller with integral action, yields an offset-free performance. When stochastic observers estimate the states of the above augmented system, it is well-known that the integral action occurs and a controller compensates the steady-state offset.18 The extended observability matrix of the augmented system Γ h Np is derived by using the extended observability matrix ΓNp as follows
where
HB ) (Γ h MT ΞM)-1
(36)
where
( )
0(i-1)l×nw H ˜ Np ) (h1 h2 · · · hNp) and hi ) Cd ΓNp-iBd
(38)
(39)
Let us consider the disturbance model in eq 36. We will use the notation xk|lw to denote the optimal estimate for xkw based on measurements up to time l. The confidence in this estimate is expressed through the conditional covariance matrix Σk|l as
Σk|l ) E{[(xkw - xk|lw)(xkw - xk|lw)T]}
(40)
The problem of recursive disturbance estimation can be viewed as that of computing the new estimate xk|kw, the covariance matrix Σk|k based on the previous estimate xk-1|k-1w, the covariance matrix Σk-1|k-1, and the new measurement yˆ k. The following Kalman filter provides the optimal disturbance estimate
xk|k-1w ) Awxk-1|k-1w xk|kw ) xk|k-1w + Kkdk
(41)
where
Kk ) Σk|k-1(CdCw)T(CdCwΣk|k-1(CdCw)T + Rv)-1 Σk|k-1 ) AwΣk-1|k-1AwT + BwRwj BwT Σk|k ) (I - KkCdCw)Σk|k-1
(42)
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4297
Combining eq 41 with the optimal prediction eq 11, we obtain the following augmented prediction model for yfa:
yˆ fa ) Lwwp + H ˜ NpΓNpw xk|kw + Luuf
∆u jf ) D hu jf - u jt
(
(48)
) ()
by introducing the auxiliary quantities
Im 0 · · · -Im Im D h ) ·· ·· 0 · · 0 · · · -Im
(43)
Also we obtain the state estimate xk|kw using the BLUE filter. The BLUE filter estimate is represented in an iterative form by
0 · · · 0 Im
ut 0 and u jt ) · · · 0
xˆ k|kw ) ΩM-1ηˆ k|k
Then the optimal control criterion can be written as
Ωj+1 ) [I + Aw-T(Ωj +
J ) (rf - yˆ fa)TQ ˆ (rf - yˆ fa) + u j fTR ˆu jf + j t)TR h ∆(D hu jf - u j t) (49) (D hu jf - u
CwT Rvj-1Cw)Aw-1BwRwj BwT]-1 Aw-T(Ωj + CwT Rvj-1Cw)Aw-1 (44)
where
ηˆ k-M+j+1|k ) [I + Aw-T(ΩM +
Q ˆ ) diag(Q · · · Q) ∈ RNl×Nl
CwT Rvj-1Cw)Aw-1BwRwj BwT]-1 Aw-T(ηˆ k-M+j|k +
R ˆ ) diag(R · · · R) ∈ RNum×Num
Cw Rvj-1dk-M+j-1) T
where Ω0 ) 0, ηˆ k-M ) 0, and j ) 0, ..., M - 1.
In this section, we study the design of a model predictive controller from another point of view. Instead of solving for the performance index J in the classical way, we start from the input/output equation. Suppose that the corresponding estimated future output sequence yˆ fa, the future reference trajectory, and the future input sequence uf are defined as
() () () rt+1 rt+2 · · ·r
t+Np
ut+1 ut+2 , and uf ) · · ·u
()
(47)
( )
Im 0 · · Λ) · 0 · · · 0
0 Im · · · 0 · · · 0
C hu j f e cj
( ) ( INu×m -INu×m D h -D h
1Nu X umax -1Nu X umin and cj ) 1 X ∆umax + u jt Nu -1Nu X ∆umin - u jt
(52)
)
(53)
The symbol X denotes the Kronecker delta product and 1Nu represents the column vector with all the elements equal to 1. We thus obtain a quadratic programming (QP) problem that can be efficiently solved using the standard numerical optimization software. 6. An Implementation Algorithm
t+Nu
where
can be reformulated as
C h ) (46)
(51)
∆umin e ∆ut+k e ∆umax
(45)
If we take into account the facts that the control horizon Nu is typically shorter than the prediction horizon Np and that the condition in eq 46 is required, it is wellknown that this also has the effect of producing less aggressive controllers. uf can be written as
ut+1 ut+2 jf ) Λ · uf ) Λu · ·u
umin e ut+k e umax
where
t+Np
∆ut+k-1 ) 0, k > Nu
R ˆ ∆ ) diag(R∆ · · · R∆) ∈ R
For the unconstrained case, the minimizing control sequence is thus obtained explicitly by means of the ordinary least-squares theory. For the constrained case, the conditions
5. An Input/Output Data-Based Predictive Control
yˆ t+1a yˆ t+2a a , rf ) yˆ f ) · · · yˆ t+Npa
(50)
Num×Num
· · · 0 · · · 0 ·· · · ·· · · · Im · · · · · · Im
The sequence of control increments can also be expressed in vector form
The input/output data-based predictive control algorithm can now be implemented as follows: (1) Construct the data block Hankel matrices Yf, Uf, and Wp from the input/output data. (2) Calculate the R component of the QR decomposition (eq 7). (3) Calculate L of eq 8, and derive the predictive controller parameters given by eq 10. (4) Calculate the weighted SVD of eq 13, and derive ˜ Np given by the extended observability matrix ΓNp and H eqs 14 and 39, respectively. (5) Construct the controller inputs
wp ) (yt-M+1T · · · ytT ut-M+1T · · · utT )T and the optimal disturbance estimate xk|kw by using the measurement correction through the Kalman filter or the BLUE filter. (6) Implement the first input ut+1 of the predictive control sequence uf by minimizing the
4298
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001
Figure 1. Input and output data sequences used for the identification of the controller parameters Lu and Lw.
performance index of eq 49 subject to the constraints of eq 52, and go to step 5. Here, we have presented an algorithm for the calculation of predictive controllers based on the input/output data only without any plant model. Although the derivation is based on expressions from the subspace identification theory, the implementation bypasses the identification step. Given the input/output data of the system, one can directly derive the prediction model from the single-step QR factorization and the SVD decomposition. 7. Simulation Study: Application to a Continuous Styrene Solution Polymerization Reactor To perform the simulation study, we consider the mathematical model of a continuous stirred-tank reactor in which styrene solution polymerization occurs. The reaction kinetics is assumed to follow the free-radical polymerization mechanism. The method of moments is adopted to calculate the number-average molecular weight (Mn) and the weight-average molecular weight (Mw). The dynamic behavior of the continuous reactor is described by the first-principles model in Na and Rhee,19 and the physical properties and kinetic paremeters are also taken from their work. The empirical gel effect correlation by Hamer et al.20 is included in the reactor model. We consider a multi-input multioutput styrene polymerization reactor system for which the monomer conversion y1 and the weight-average molecular weight y2 are controlled by manipulating the feed flow rate u1 and the jacket inlet temperature u2. Figure 1 presents the input/output data sequence used for the identification of the controller parameters Lu and Lw. The sampling period is 1 min. The input signals are drawn from a uniform distribution with four
levels. The pseudo-random four-level input signals excite nonlinear modes that pseudo-random binary signals (PRBS) cannot.19 Because of these features, we use four-level input signals rather than the PRBS. The output data for the pseudo-random four-level input signals are obtained by the mathematical model. All of the simulation studies were conducted with the same tuning parameters under the same conditions. Reference conditions for the identification and controlcan be found in Na and Rhee.19 The upper and lower bounds of the jacket inlet temperature were 90 and 50 °C, respectively, and those of the feed flow rate are 25 and 5 mL/min, respectively. The sampling time of the discrete-time control was 1 min. The prediction horizon and control horizon were chosen to be 15 and 7 min, respectively, by trial and error. The output and input velocity weighting factors were chosen as Q ) diag(1, 1) and R∆ ) diag(4, 4), respectively. The disturbance model matrices were introduced as design variables, and all of them except for Cd were chosen as identity matrices. The matrix Cd was chosen as diag(1.5, 1.5) by trial and error. The number of disturbance states was taken as 2, which is equal to the number of outputs. When the system reached the steady state, the controller designed a priori by the input/output data started to perform the tracking control. All of the controllers had the same tuning parameters and the same control parameters Lw and Lu. Figure 2 shows the servo response of the input/output data-based predictive controller without integral action for a multilevel setpoint change. We observe that steadystate offsets exist for both conversion and Mw. This poor control performance is natural because the controller is obtained on the basis of the averaged linear model identified under various operating conditions. To reduce the offset, we added disturbances of integrating characteristics to the states. Figures 3 and 4
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4299
Figure 2. Tracking control performance of the input/output data-based predictive controller without integral action.
Figure 3. Tracking control performance of the input/output data-based predictive controller with integral action using the Kalman filter.
present the servo responses of the two controllers using the Kalman filter and the BLUE filter, respectively, to estimate the integrated white-noise disturbances. The initial value of the unmeasured disturbance was unknown. Thus, we arbitrarily chose the mean value and covariance as the zero vector and 0.01I2, respectively. In the case of the BLUE filter, the information about
the initial disturbance state is unnecessary, and the other system matrices of the augmented system are chosen to be the same for both filters. The steady-state offset of Mw is eliminated, whereas that of conversion remains but is reduced. This phenomenon is due to the high nonlinearity of the styrene polymerization reactor system. When the nonlinear input/output predictive
4300
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001
Figure 4. Tracking control performance of the input/output data-based predictive controller with integral action using the BLUE filter.
controller is applied to the same styrene polymerization reactor, the steady-state offset is satisfactorily eliminated.21 The method of adding integrated white-noise disturbances to the state with the Kalman filter is wellestablished.18 We also observe that the BLUE filter shows an equivalent effect. 8. Conclusions The state estimation scheme in subspace identification is analyzed to show that its estimate is biased but that the gain matrix has a structure similar to that of the BLUE filter obtained by solving the least-squares objective function analytically. It turns out that the state estimate obtained by the subspace identification algorithm indeed approaches the estimate obtained by BLUE filter when the norm of the state covariance matrix becomes large or the norm of the noise covariance matrix becomes small. On the basis of such an interpretation of the state estimation, a linear input/output data-based predictive control with integral action was designed by introducing an integrated white noise model and using each the BLUE filter and the Kalman filter as a stochastic observer for the unmeasured disturbance. The advantage of such a controller is that there is no need to calculate the state estimate explicitly from an input/ output data-based prediction model and thus the controller can be calculated directly from the input/output data in a single step. However, it is necessary to correct the future output sequence yˆ f by using the integrated white-noise disturbance because the controller parameters, Lw and Lu, contain inaccurate information about the state. The linear input/output data-based predictive controller proposed in this study is corroborated by applying it to the property control of a styrene solution polymerization reactor. As a result, it is clearly demon-
strated that the control performance is significantly improved with each the BLUE filter and the Kalman filter. Acknowledgment Financial aid from the Brain Korea 21 Program supported by the Ministry of Education of South Korea is gratefully acknowledged. Literature Cited (1) Lee, J. H.; Morari, M.; Garcia, C. E. State-Space Interpretation of Model Predictive Control. Automatica 1994, 30, 707. (2) Hjalmarsson, H.; Gevers, M.; Gunnarsson, S.; Lequin, O. Iterative Feedback Tuning: Theory and Applications. IEEE Control Syst. Mag. 1998, 18, 26. (3) Hjalmarsson, H. Efficient Tuning of Linear Multivariable Controllers Using Iterative Feedback Tuning. Int. J. Adapt. Control 1999, 13 (7), 553. (4) Amirthalingam, R.; Lee, J. H. Subspace Idedtification Based Inferential Control of a Continuous Pulp Digester. Comput. Chem. Eng. 1997, 21, 1143. (5) Ruscio, D. D.; Foss, B. On state space model based predictive control. Presented at the 5th IFAC Symposium on Dynamics and Control of Process Systems, Corfu, Greece, June 8-10, 1998. (6) Amirthalingam, R.; Lee, J. H.; Sung; S. W. A Two-Step Procedure for Data-Based Modeling for Inferential Predictive Control System Design. AIChE J., manuscript submitted. (7) Favoreel, W.; De Moor, B.; Gevers, M.; Van Overschee, P. Model-free Subspace-based LQG-design; Technical Report ESATSISTA-TR-1998-34; Department of Electrical Engineering, Katholieke Universiteit: Leuven, Belgium, 1998. (8) Favoreel, W.; De Moor, B.; Van Overschee, P.; Gevers, M. Model-free Subspace-based LQG-design. Proceedings of the American Control Conference; American Institue of Chemical Engineers: New York, 1999; 3372. (9) Favoreel, W.; De Moor, B. SPC: Subspace Predictive Control; Technical Report 98-84; Katholieke Universiteit: Leuven, Belgium, 1998. (10) Van Overschee, P.; De Moor, B. Subspace Identification for Linear Systems: Theory, Implementation, Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996.
Ind. Eng. Chem. Res., Vol. 40, No. 20, 2001 4301 (11) MacGregor, J. F.; Kaeckle, C.; Kiparissides, C.; Kouto, M. Process Monitoring and Diagnosis by Multiblock PLS Methodology. AIChE J. 1994, 40, 826. (12) MacGregor, J. F.; Kourti, T.; Nomikos, P. Analysis, Monitoring and Fault Diagnosis of Industrial Processes Using Multivariable Statistical Projection Methods. Proceedings of IFAC World Congress; Elsevier: New York, 1996; Vol. M7a-05(1). (13) Robertson, D. G.; Lee, J. H.; Rawlings, J. B. A Moving Horizon-Based Approach for Least-Squares Estimation. AIChE J. 1996, 42 (8), 2209. (14) Russel, S. A.; Kesavan, P.; Lee, J. H.; Ogunnaike, B. A. Recursive Data-based Prediction and Control of Batch Product Quality. AIChE J. 1998, 44 (10), 2442. (15) Ruscio, D. D. Model predictive control and identification: A linear state space model approach. (A) Proceedings of the 36th Conference on Decision and Control; Wiley & Sons: New York, 1997; p 3202. (16) Ruscio, D. D. Model based predictive control: An extended state space approach. (B) Proceedings of the 36th Conference on Decision and Control; Wiley & Sons: New York, 1997; p 3210. (17) Kwon, W. H.; Kim, P. S.; Park, P. A Receding Horizon
Kalman FIR filter for discrete time-invariant systems. IEEE Trans. Autom. Control 1999, 44 (9), 1787. (18) Lee, J. H.; Ricker, N. L. Extended Kalman Filter Based Nonlinear Model Predictive Control. Ind. Eng. Chem. Res. 1994, 33, 1530. (19) Na, S. S.; Rhee, H.-K. Polynomial ARMA Model Identification for a Continuous Styrene Polymerization Reactor Using Online Measurements of Polymer Properties. J. Appl. Polym. Sci. 2000, 76, 1889. (20) Hamer, J. W.; Akarmov, T. A.; Ray, W. H. The Dynamic Behavior of Continuous Polymerization. II. Nonisothermal Solution Homopolymerization and Copolymerization in a CSTR. Chem. Eng. Sci. 1981, 36, 1897. (21) Yoo, K.-Y.; Rhee H.-K. An Input/Output Data-based Predictive Controller Design for Nonlinear Systems. AIChE J., manuscript submitted.
Received for review November 30, 2000 Revised manuscript received April 16, 2001 Accepted April 18, 2001 IE0010159