1938
Ind. Eng. Chem. Res. 2003, 42, 1938-1947
Identification and Control of Processes with Periodic Operations or Disturbances Yangdong Pan School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907-1283
Jay H. Lee* School of Chemical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100
Periodically time-varying (PTV) systems are commonly used to model natural and forced periodic phenomena encountered in many engineering applications. We propose an approach in which the subspace identification method is used to identify the “lifted” form of a PTV system. Then, the identified lifted system is reformulated as a PTV system in an augmented form for the design of a PTV Kalman filter and a PTV linear quadratic regulator that can eliminate periodic error without offset. We touch upon some important practical issues of identifying a lifted system and also of designing an optimal controller with the identified model. We illustrate the application of the proposed approach through two examples involving a polyamine reactor and a distillation column. 1. Introduction Periodically time-varying (PTV) systems are often used to model natural or forced periodic phenomena occurring in various engineering applications. In the process industries, many processes experience periodic disturbances due to natural cycle times of upstream processes or other cyclical environmental influences such as diurnal temperature fluctuations. In wastewater treatment plants, for example, the feed flow rate and composition can exhibit strong diurnal variations.1 In some applications, forced periodic operations can be used either to improve selectivity and yields2 or to make a continuous operation feasible.3,4 A potential drawback of adopting such a mode of operation is that operation becomes more complex and difficult to control. Often, processes with periodic characteristics show strong nonstationary or cyclostationary behavior. Also, because of the wide envelope of operating space covered by a periodic operation, the process dynamics can be strongly nonlinear. These factors make the control of periodic processes more complicated than that of steady processes. As an alternative to a full-blown nonstationary, nonlinear model, which can lead to difficult identification and controller design problems, a linear timevarying model valid along some nominal periodic trajectory offers a good compromise, given the existence of rich theory and design techniques for the control of linear PTV systems. Many such techniques are based on the idea of “lifting”, which allows the PTV system to be represented as an equivalent system in time-invariant form.5-7 To use the available tools, however, one needs an accurate PTV system description of the underlying process. If a fundamental nonlinear model is available, one can linearize the model along a nominal trajectory to obtain a PTV system approximation around it. In most practical situations, however, a fundamental * To whom all correspondence should be addressed. Phone: (404)385-2148. Fax: (404)894-2866. E-mail: Jay.Lee@ che.gatech.edu.
model is not readily available, and system identification has to be performed. Hence, it is fair to say that the practicability of the available control methods hinges on the ability to identify an accurate PTV system model in an efficient manner. In contrast to the large volume of research on analysis and control of PTV systems, the literature on PTV system identification is minimal. In Bittani et al.,8 the prediction-error method is used to obtain a periodic ARMA model. Verhaegen and Yu8 discuss a subspace identification method for general time-varying MIMO systems. This technique, which is based on the MOESP algorithm developed originally for linear time-invariant (LTI) systems, can be applied to PTV systems with only minor modifications. However, for a PTV system, the complexity is significantly increased by the calculation of multiple Hankel matrices and subsequent QR and singular value decompositions. Also, if a combined deterministic and stochastic system description is desired, the procedure would be further complicated. Given the equivalence of a PTV system and its lifted form, one can take the approach of identifying the system in the lifted form first. Then, the identified lifted system needs to be realized as a PTV system for the controller design. One apparent advantage of this approach is that, because a lifted system is in the form of an LTI system, any MIMO identification method for LTI systems can be used. For a similar problem, Li et al.10 proposed this approach for the identification of multirate systems. In their method, subspace identification was used to obtain a model for the lifted form of a multirate system. Furthermore, they developed a procedure to extract a fast-rate LTI system model from the identified lifted system, which could be useful for inferential control. The problem of identifying multirate systems is considerably simpler than the problem of identifying general PTV systems, because the underlying system in the former problem is assumed to be time-invariant, with multirate sampling being solely responsible for the system’s periodic characteristic. The more challenging aspect of the latter approach is the realization of the lifted system as a PTV system.
10.1021/ie020313y CCC: $25.00 © 2003 American Chemical Society Published on Web 04/01/2003
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1939
For the purpose of designing periodic controllers, several researchers have developed methods to realize a transfer matrix as a PTV state-space model.11-13 In particular, Lin and King11 proposed an algorithm to obtain a periodic realization from a lifted state-space model. Their algorithm can be adopted here. However, the realization algorithm is computationally complex, involving multiple rank decompositions and generalized inverse calculations. In addition, the realized periodic system is usually nonminimal, with a state dimension larger than that of the lifted system. Given the complexity involved in directly identifying a PTV system model or realizing an identified lifted system model as a PTV system, a worthy question to explore is whether the lifted form can be directly used for on-line control. We show in this paper that such an approach is indeed possible if the lifted system model is rewritten in recursive form. What results from this simple reformulation of the model is a standard PTV system structure, to which the wealth of existing optimal estimation and control methods can readily be applied. We also show that the augmented PTV system structure derived from the lifted model is particularly useful for designing a controller with periodwise integral action, which enables the rejection of periodic errors without offset. Additionally, we touch upon a few important practical issues of identifying a lifted system and designing an optimal controller based on the identified model. Finally, we apply the proposed procedure to a polyamine reactor and a distillation column in which the feed flow rates experience significant periodic swings.
2.2. Lifted-System Description. We first note that the PTV system of eqs 1 and 2 is equivalent to the following lifted system, which describes the system dynamics over the period index alone
xk+1(0) ) Φxk(0) + Γ1Uk + Γ2Ek
where the lifted vectors Uk, Yk, and Ek are defined as
[ ] [ ]
uk(0) yk(0) l l Uk } , Yk } , l l uk(N-1) yk(N-1)
[ [
(5)
(A(N-1) ‚‚‚ A(1) B(0)) (A(N-1) ‚‚‚ A(2) B(1)) Γ1 } ‚‚‚ (A(N-1) B(N-2)) B(N-1)
(6)
[
C(0) C(1) A(0) l Π } l C(N-2) A(N-3) ‚‚‚ A(0) C(N-1) A(N-2) ‚‚‚ A(0)
yk(t) ) C(t) xk(t) + D(t) uk(t) + ek(t), t ) 0, ..., N - 1 (2)
[
D1 }
D(0) (C(1) B(0)) (C(2) A(1) B(0)) l l (C(N-1) A(N-2) ‚‚‚ A(1) B(0))
[
] ]
Φ } A(N-1) A(N-2) ‚‚‚ A(0)
(A(N-1) ‚‚‚ A(1) K(0)) (A(N-1) ‚‚‚ A(2) K(1)) Γ2 } ‚‚‚ (A(N-1) K(N-2)) K(N-1)
(1)
In the above equations, x is the state, u is the input, and y is the output. The system is given in innovation form, where ek(t), which represents the prediction error, is a white noise sequence of covariance Re(t). k is the period index, and t ∈{0, ..., N - 1} is the time index. Hence, we assume that there are N time sample points within each period. The dimensions of u and y can be allowed to vary with t for more generality. This way, the above equations can describe multirate sampled data systems, in which case only B and C are timevarying. In this paper, for the convenience of exposition, we assume that the dimensions remain constant with t. The extension of the methodology to the case of varying dimension will be self-apparent. We assume that the system is stabilizable and detectable. For precise definitions of structural properties of a PTV system, such as stabilizability and detectability, readers are referred to a review paper by Bittanti.14 The objective is to identify, from collected inputoutput data, a model equivalent to the above system in the input-output sense for the purpose of on-line prediction and control.
(4)
and the system matrices are given as
2.1. Problem Statement. We assume that underlying the problem is the following discrete-time PTV linear stochastic system
xk+1(0) ) xk(N)
[ ]
ek(0) l Ek } l ek(N-1)
2. Preliminaries
xk(t+1) ) A(t) xk(t) + B(t) uk(t) + K(t) ek(t)
(3)
Yk ) Πxk(0) + D1Uk + D2Ek
0 D(1) C(2) B(1) l l ‚‚‚
‚‚‚ 0 D(2) l l ‚‚‚
(7)
]
‚‚‚ ‚‚‚ 0 D(3) l
(8)
‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚
0 0 l 0 l D(N-1)
]
(9)
D2 } I (C(1) K(0)) (C(2) A(1) K(0)) l l (C(N-1) A(N-2) ‚‚‚ A(1) K(0))
0 I C(2) K(1) l l ‚‚‚
‚‚‚ 0 I l l ‚‚‚
‚‚‚ ‚‚‚ 0 I l
‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚
0 0 l 0 l I
]
(10)
If the original PTV system is stabilizable and detectable, so is the lifted system. Note that eq 3 can be
1940
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003
rewritten in the standard innovation form as
xk+1(0) ) Φxk(0) + Γ1Uk + Γ ˜ 2Ek Yk ) Πxk(0) + D1Uk + Ek
(11)
noise Ek as an integrated white noise (of dimension N). In this case, the standard subspace identification method that assumes a stationary system can give poor results. A practical remedy is to difference the data with respect to the period index k. If we define ∆Uk ) Uk - Uk-1, we can express eq 11 equivalently as
where Ek ) D2Ek and Γ ˜ 2 ) Γ2D2-1.
∆xk+1(0) ) Φ∆xk(0) + Γ1∆Uk + Γ ˜ 2∆Ek ∆Yk ) Π∆xk(0) + D1∆Uk + ∆Ek
3. Identification of PTV Systems In this section, we first discuss a few issues regarding the identification of lifted systems. The PTV realization algorithm in ref 11 is briefly discussed, and some disadvantages are pointed out. We then show how the identified lifted system model can be rewritten to allow for a recursive update of the state at each sample time. 3.1. Identification of the Lifted System. Note that the lifted system in eq 11 is in the standard innovation form of an LTI stochastic system. Therefore, the subspace identification method can be applied to identify it. However, the size of the lifted input and output vectors can be huge, which complicates the identification. A few other issues also require attention. 3.1.1. Assuring Causality. Note that, because we are dealing with causal systems, D1 in eq 11 has to be restricted to a block lower triangular form. This restriction is easy to implement in the least-squares step of subspace identification. For example, in N4SID,15 states are first constructed from the input-output data, and then the system matrices are obtained through linear least squares. In estimating the system matrices Π and D1, one can perform the least-squares calculation for one block row (corresponding to one sample time) at a time. In this way, the elements of the lifted input vector that correspond to zero blocks can be dropped from the regressor vector to impose the block lower triangular structure on D1. 3.1.2. Test Signal Design. The lifted form also offers some useful insights into the design of the test signal. The test input signals are to be designed as perturbations to the respective nominal trajectories. By viewing the lifted system as a MIMO discrete-time LTI system, one can use available guidelines and techniques for designing test signals for MIMO system identification. For example, we can easily see that a PRBS of period N (or any periodic signal of period N) will result in a Uk that repeats over k. Within the MIMO LTI system analogy, this would correspond to using a step signal of fixed magnitude and randomly chosen sign for each input channel. Hence, we see that, with such a signal, we can expect to obtain a lifted model that is accurate only at low frequencies. To identify a lifted system over all frequencies of interest, a PRBS of much larger period (many multiples of N) is needed in general. We can also see that identifying a PTV system model with a large period is tantamount to identifying an LTI MIMO system of very large dimension, which is inherently a difficult problem requiring a large set of test data. 3.1.3. Differenced Form. Disturbances for periodic systems are often assumed to be cyclostationary, but in many cases, their mean dynamic behavior over a period can change from one period to another. For example, a step change in the feed to a periodically operated process would manifest itself as a sustained change in the converged periodic behavior. Mathematically, such a case would be better described by a lifted system in the same form as eq 11 but with the external
(12)
Notice that, if Ek is an integrated white noise sequence, ∆Ek is a white noise. Hence, the assumption of a stationary system is more reasonable for the differenced form of the model, and the standard subspace identification methods can be used effectively for this form. 3.2. PTV Realization of a Lifted Form. With the lifted system model in its current form, one can perform the measurement update only at the end of each period because measurements of an entire cycle are needed. When the period is large compared to the sample interval, this will limit the ability to reject disturbances that occur in the middle of a period. Hence, it is desired to convert the lifted system into a PTV system that can be used for prediction and control at every sample time. The natural inclination would be to attempt to extract the PTV model matrix sequences [A(i), B(i), C(i), D(i), K(i), Re(i), i ) 0, ..., N - 1] from the lifted system model matrices (Φ,Γ1, Γ ˜ 2, Π, D1, RE). This would be a formidable task in general, as the inverse step is nonunique and the relationship is nonlinear. The realization algorithm by Lin and King11 can be used to convert the lifted system into an equivalent PTV system. Even though their algorithm was developed in the context of deterministic systems, it extends to the case of combined deterministic-stochastic systems, which we have, in a straightforward manner. The main drawback is that the procedure is computationally complex, involving multiple rank decompositions and calculations of generalized inverses. Furthermore, a minimal realization is generally not possible, and uncontrollable/unobservable modes will result, thus complicating the subsequent controller design. Finally, Re(i) extracted from RE will not be exact, because RE obtained from the identification generally will not be block diagonal. 3.3. Converting the Lifted System into a PTV System. As an alternative, we look for a way to use the lifted model as is. To this end, we define new state z as follows
[ ] [ ]
Γ0+ (t) xk+1(0) zk(t) } - 10+ U0+(t) Yk D1 (t) k where
(13)
[ ]
uk(t) ) l U0+ k uk(N-1)
Γ (t) ‚‚‚ Γ1(N-1) ] Γ0+ 1 (t) ) [ 1
(14)
D (t) ‚‚‚ D1(N-1) ] D0+ 1 (t) ) [ 1 and Γ1(i) and D1(i) denote the ith block columns of Γ1 and D1 [corresponding to uk(i)], respectively. From this
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1941
definition, the following recursive relationship can be established
zk(t+1) ) zk(t) +
ζk(t+1) ) ζk(t) +
[ ]
ζk+1(0) )
[ ]
[ ]
Γ1(t) u (t), t ) 0, ..., N - 1 (15) D1(t) k
[ ]
xk+1 Also, because zk(N) ) Y according to the definition, k the lifted system equation of eq 11 is translated into
[ ]
[ ]
Γ ˜ Φ 0 zk+1(0) ) zk(N) + 2 Ek+1 Π 0 I
[ ]
Φ 0 K ζ (N) + Ek+1 Π I k I
(19)
where
ζk(t) ) (16)
[
] [ ]
∆xk+1(0) Γ0+(t) - 10+ ∆U0+ k (t) Yk(t) D1 (t)
(20)
and
The output equation for each sample time can be written as
yk(t) ) [0 H(t) ]zk(t) + D(t)uk(t)
Γ1(t) ∆uk(t), t ) 0, ..., N - 1 D1(t) (18)
(17)
where H(t) is simply a matrix that extracts yk(t) from Yk. Equations 15-17 are equivalent to the original (or lifted) system in the input-output sense and comprise a standard PTV state-space system that can be used for controller design. 4. Control For a periodic system, the standard control objective is to track periodic reference trajectories and/or to reject disturbances having periodic effects on the controlled variables. The effectiveness of conventional feedback controllers such as the proportional-integral-derivative (PID) type or dynamic matrix control (DMC) type on such problems can be quite limited in general because the system might not settle down to a constant condition. For servo-mechanical systems that are designed to track periodic reference trajectories, the technique called repetitive control (RC) has seen much use.16 The idea is based on the internal model principle, which tell us to embed into the controller an element that has poles on all of the harmonic and subharmonic frequencies of the period. However, the assumption of an unconstrained linear time-invariant system inherent in the development of RC is rather limiting for chemical process applications. Aside from RC, extensions of optimal control methods such as LQR and LQG to periodic systems have been developed.17 However, most optimal controllers, when formulated as in the literature, would lead to periodic offsets. Lee and co-workers18 addressed this issue and proposed an approach based on model predictive control (MPC). This technique, which they referred to as R-MPC (R for repetitive), embeds into the controller the integral action over period index and thereby ensures convergence to a periodic set-point trajectory and offset-free rejection of periodic disturbances. Here, we briefly describe the R-MPC approach, but we develop it as an LQG controller rather than an MPC controller. 4.1. Embedding Periodwise Integrators in the Model. Integral action allows for offset-free control. Because we are interested in periodic reference trajectories and disturbances, the integral action should be created with respect to run index k. To design a controller with the desired integral action, it is convenient to embed integrators directly into the model. The PTV system in eqs 15-17 does not have such integrators, but it can be easily re-expressed as
yk(t) ) [0 H(t) ]ζk(t) + D(t) ∆uk(t) )H h ζk(t) + D(t) ∆uk(t)
(21)
Notice the integrators in the transition matrix in eq 19. These integrators are detectable and also stabilizable if the original system is stabilizable and there are at least as many linearly independent inputs as outputs, i.e., the steady-state gain matrix of the lifted system in eq 3 has full row rank. It is worth mentioning that the stochastic part of the model might or might not come from the system identification. If the disturbances experienced during the identification experiment are indeed similar in nature to those experienced during real operation and enough data were collected to allow for a good estimation of the stochastic part, the identified model can be used directly. In this case, K is simply equal to Γ ˜ 2. If this is not the case, the system matrices for the stochastic part can be chosen to reflect those disturbances one is most interested in rejecting. For example, with K ) 0, periodic disturbances, the values of which at various sample times change from one period to the next in a completely random manner, are assumed to be added directly to the deterministic output trajectories. This choice enables the rejection of any persistent periodic disturbance without offset, if the embedded integrators are stabilizable. For better transient results, the covariance of E also can be chosen to express the expected correlation among the intraperiod samples. For example, if the deviation from the previous run’s trajectory is assumed to be similar to an integrated white noise (over the index t), the covariance could be chosen in the form
[
I I RE ) R l I
I 2I l 2I
‚‚‚ ‚‚‚ l ‚‚‚
l 2I l N‚I
]
(22)
where I is an identity matrix of size ny, which is the output dimension. 4.2. PTV Kalman Filter for State Estimation. The time-varying Kalman filter equation for the periodic system in eqs 18-21 is of the form
[ ]
Γ1(t) u (t) + Kk(t){yk(t) D1(t) k [H h ζˆ k(t) + D(t) uk(t)]}, t ) 0, ..., N - 1
ζˆ k(t+1) ) ζˆ k(t) +
[ ]
Φ 0 ζˆ k+1(0) ) ζˆ (N) Π I k
(23)
1942
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 ∞ N-1
where
Kk(t) ) Pk(t) H h T(t) [H h (t) Pk(t) H h T(t)]-1, t ) 0, ..., N - 1 (24) and the error covariance matrix sequence is calculated from
∑ ∑ {||ek(t)||Q + ||∆uk(t)||R ) || ζ˜ k(t)||Q (t) + ∆u (t)k)1 t)0 min
ζ˜
k
|∆uk(t)|R} (33)
where the notation ||x||Q represents xTQx and Ek ) Yk - Rk. Also
h (t) Pk(t), t ) 0, ..., N - 1 Pk(t+1) ) Pk(t) - Kk(t) H Pk+1(0) )
[ ] [ ] [ ] [ ] T
Φ 0 Φ 0 P (N) Π I k Π I
+
(25) T
K K RE I I
Assuming that the system is detectable, the covariance should converge to a periodic solution.19 Hence, for simplicity, one can choose to implement the periodic steady-state solution
K∞(t) ) P∞(t) H h T(t)[H h (t) P∞(t) H h T(t)]-1, t ) 0, ..., N - 1 (26)
Qζ˜ (t) ) [0 H(t) ]TQ[0 H(t) ] Note that, because
[ ]
Γ ζ˜ k(N) ) ζ˜ k(0) + D1 ∆Uk 1
[ ]
[ ] [ ] [ ] Γ ˜ Γ ˜ + 2 RE 2 I I
T
Φ 0 K ζ + Ek+1 Π I k I
∞
min
∑
∆Uk k)1
the solution P is equal to P∞(0).19 Once P∞(0) is obtained in this manner, the rest of the values P∞(1), ..., P∞(N) can be easily calculated by propagating eq 25 for one period. An alternative is to run the time-varying Riccati equation until it converges to the periodic solution. 4.3. Periodic Linear Quadratic Regulator. Let us redefine the state to express the output as a deviation from a periodic reference trajectory
[ ]
(29)
ζ˜ k+1(0) )
[ ]
[ ]
Γ1(t) ∆uk(t), t ) 0, ..., N - 1 D1(t) (30)
[ ]
[]
Φ 0 K 0 ζ˜ (N) + Ek+1 ∆Rk+1 (31) Π I k I I
The following quadratic control objective is adopted
Ek
Q
+ ∆Uk
R
)
ζ˜ k(0) ∆Uk
Σ
(37)
and
Σ)
[
ˆ TQD1 C ˆ TQC C D1TQC D1TQD1 + R
]
(38)
The solution can be obtained by solving an algebraic Riccati equation that has been modified to account for the cross-weighting term.20 From this solution, S∞, the solution of the Riccati equation, is obtained. Step 2. Because the cost of all runs from the (k + 1)th on under the LQR can be expressed as ζ˜ k+1(0)S∞ζ˜ k+1(0), the infinite-horizon problem is equivalent to
where ∆Rk+1 ) Rk+1 - Rk. Note that
ek(t) } yk(t) - rk(t) ) [0 H(t) ]ζ˜ k(t), t ) 0, ..., N - 1 (32)
{|| || | | ||[ ]||
where
where Rk represents the lifted vector containing the values of the reference trajectory at the sample times over period k. Then, the PTV system of eqs 18-21 can be rewritten in terms of the new state as
ζ˜ k(t+1) ) ζ˜ k(t) +
[]
which is expressed in terms of the k index only To solve the above infinite-horizon control problem, the following procedure can be used: Step 1. Solve the following infinite-horizon problem for the lifted system of eqs 36
(28)
0 ζ˜ k(t) ) ζk(t) - R k
[] []
Ek ) C ˆ ζ˜ k(0) + D1∆Uk
[ ] [ ] Yk ) [0 I ]ζk
[ ][ ] [ ]
(27)
It turns out that, if one solves the algebraic Riccati equation corresponding to the Kalman filter for the lifted system
ζk+1 )
[ ]
Φ 0 Φ 0 Γ1 K ζ˜ (0) + ∆Uk + Ek+1 Π I k Π I D1 I 0 ∆Rk+1 I K 0 )A ˆ ζ˜ k(0) + B ˆ ∆Uk + Ek+1 ∆Rk+1 (36) I I Ek ) [0 I ]ζ˜ k(0) + D1∆Uk
ζ˜ k+1(0) )
P∞(t+1) ) P∞(t) - K∞(t) H h (t) P∞(t), t ) 0, ..., N - 1 T
(35)
eq 31 can be written as
where the error covariance matrix sequence is the solution to the following PTV Riccati equation
Φ I Φ I P (N) P∞(0) ) Π 0 ∞ Π 0
(34)
N-1
min
∑ {||ζ˜ k(t)||Q (t) + ||∆uk(i)||R} + ζ˜ k+1(0)S∞ζ˜ k+1(0)
∆uk(t) t)0
Also, note that
ζ˜
(39)
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1943
include only the primary variables
ζ˜ k+1T(0)S∞ζ˜ k+1(0) ) ζ˜ kT(N)A ˆ TS∞A ˆ ζ˜ k(N) ) ζ˜ kT(N) S(N) ζ˜ k(N)
(40)
Hence, the above infinite-horizon problem reduces to a finite-horizon problem of a single period. We can run the Riccati difference equation for the above finite horizon problem and the PTV system in eq 30, starting with the terminal weight of S(N) to obtain S(N-1), ..., S(1) and also the optimal periodic state feedback law for ∆u(N-1), ..., ∆u(0) sequentially. Although the control problem is formulated as an unconstrained tracking problem here, this need not necessarily be so. Constraints can be entered explicitly into the formulation, as in MPC,21 which leads to an optimization problem that must be solved on-line at each time step. An economic objective and quality constraints can also be considered, as in ref 22. 4.4. Inferential Control. In many chemical processes, property variables cannot be measured at all or can be measured only with large delays and/or at unacceptably low sampling frequencies. In such cases, secondary variables, such as pressure and temperature, that are more easily measured can be used to infer the property variables. This is referred to as the inferential control approach. Here, we show that inferential control also can be implemented within the proposed PTV system identification and control framework. First, the following lifted model is identified for both the primary and secondary variables
˜ 2Ek ∆xk+1(0) ) Φxk(0) + Γ1∆Uk + Γ ∆Yk ) Πxk(0) + D1∆Uk + Ek where
[ ]
ypk(0) ysk(0) Yk ) l ypk(N-1) ysk(N-1)
(41)
(42)
ypk are the primary variables, and ysk are the secondary variables. The Kalman filter can be designed by proceeding as before. The measurement equations can be defined for those samples for which measurements become available at each sample time. For example, if only secondary measurements are to be used for real-time estimation, the measurement equation becomes
ysk(t) ) [0 Hs(t) ]ζk(t) + D(t) uk(t), t ) 0, ..., N - 1 (43) where Hs(t) is a matrix that extracts ysk(t) from Yk. If the primary measurements with delays are also used, the state can be further augmented with the samples of previous cycles so that model estimates for the delayed measurements can be extracted from the current state to form the innovation term. In the LQR design, the system in such a situation will typically be nonstabilizable because an independent integrator is added for each ypk(t) and ysk(t). Let us assume that there is no control requirement for ysk(t). Then, for the LQR design, we can redefine the state to
[
where
] [
]
∆xk+1(0) Γ0+(t) - 1p0+ ∆Uk0+(t) p p Yk - Rk D1 (t)
ζ˜ pk(t) }
[ ] [ ]
ysk(0) l , Rpk ) Ypk ) l ysk(N-1)
rpk(0) l l rpk(N-1)
(44)
(45)
Then, the controller can be constructed using the following lifted model form p ζ˜ k+1 (0) )
[ ]
[ ][ ]
Φ 0 p Φ 0 Γ1 ζ˜ k(0) + p ∆Uk p I Π Πp I D 1
)A ˆ ζ˜ pk(0) + B ˆ ∆Uk Epk ) [0 I ]ζk(0) + Dp1∆Uk
(46)
)C ˆ ζ˜ k(0) + Dp1∆Uk where Πp and Dp1 represent the rows of Π and D1, respectively, corresponding to the primary variables ysk. Using the above model, the controller design can proceed as before. It is noteworthy that, for inferential control, it is essential to use a stochastic model derived from data to reflect the correlation between the primary and secondary variables. The choice of K ) 0, which is effective in the case that all outputs are measured efficiently, does not enable inferential control, as this would make the assumed disturbance for each output channel be independent. In the case that plant test data are not collected long enough to support accurate identification of a stochastic model, one can use regular operation data for the model and then combine the resulting stochastic model with the deterministic model identified from the plant test data.23 5. Simulation Result To demonstrate the effectiveness of the proposed approach, two example processes were chosen. The first example is a polyamine reactor, and the second is a moderate-purity distillation column. 5.1. Example 1: Polyamine Reactor. The polyamine process (illustrated in Figure 1 produces an organic intermediate through a polycondensation reaction. The polyamine reactor contains a single well-mixed liquid phase in which a set of competing reactions occur. Although oligomerization reactions result in products of various chain lengths, only the product with specific end groups and a chain length of 2 is desired. The original model was developed by Bayer AG24 and contains 16 components and 53 irreversible second-order reactions. The model used in this paper is an extension of the original model by Till,25 who included some additional features such as the energy balance to express the effect of cooling on the reactor temperature. Recently, a dynamic optimization was performed with the original polyamine reactor model by Cruse et al.26 It was found that periodic switching between some high and low levels of the feed leads to better selectivity than
1944
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003
Figure 3. Control result for the case where the controller is initialized with a constant-cooling profile. No disturbance was assumed. (Sample time ) 55 s.)
Figure 1. Schematic representation of the polyamine reactor.
Figure 2. Optimal flow rate profile and heat of reaction for the polyamine reactor.
steady-state operation. The selectivity was improved by 8% while the same average production rate was maintained as in the steady-state operation. The optimal feed profile was obtained using the assumption of a constant reactor temperature. Figure 2 shows the optimal trajectories of the input feed streams of excess reactant fin1 and monomer fin2 and of the product stream fout, as well as the heat produced by the chemical reactions in the isothermal operation scenario. The optimal period was found to be 1100 s. The optimal temperature was found to be 110 °C. Till et al.27 performed a sensitivity analysis to show that a deviation of just a few degrees in the reactor temperature could cancel the selectivity improvement afforded by the periodic operation. Motivated by this finding, they proceeded to design a model predictive controller (repetitive MPC) to control the reactor temperature at the fixed set point of 110 °C. The controller manipulated the coolant flow rate to the jacket to control the temperature. They showed that, starting with a constant cooling rate that led to a temperature fluctua-
tion of more than (5 °C, the repetitive MPC controller was able to shape the cooling profile gradually to control the temperature very tightly. The performance was clearly superior to that of a PID controller or an MPC controller of the conventional type. Even though it delivered a good performance, their controller was based on linearization of a fundamental reactor model along the open-loop periodic trajectory. In practice, such a fundamental model might not exist or might be too complex for direct use in control. Hence, we examined the possibility of using an identified model. The proposed PTV system identification and control approach is applied here, and its performance is checked under different scenarios. First, we consider the scenario of starting with a constant cooling profile and then manipulating it to control the temperature at the constant set-point value of 110 °C. Note that, with the constant cooling profile, the temperature fluctuates by more than (5 °C because of the periodic forcing of the feed flow rate. The identification data are generated by perturbing the constant cooling profile with a random signal of magnitude 6 g/s. From the data, the nominal trajectory values (corresponding to the constant cooling rate) are subtracted to express the data in terms of deviations from the nominal trajectory. The deviation data are then used to identify a lifted model by using the modified N4SID method as we described earlier. It is found that a five-state lifted model captures most of the data trend. After converting the lifted model into the augmented PTV model form, a Kalman filter and an LQR are designed and combined into a controller. In the design, the noise input matrix K is set to 0; hence, a simple periodwise integrating disturbance model is assumed in the design. The designed controller is an unconstrained controller; any input outside the allowed range is scaled to the maximum or minimum value. The period was chosen to be 1100 s, and the sample time was set to be 55 s (hence, N ) 20). Figure 3 compares the closed-loop control result to the open-loop result with the constant cooling rate for 15 runs. It is shown that the cooling rate converges to a fixed periodic trajectory after just a few runs. Note that some periodic tracking error persists because of the input constraints. However, the tracking error is substantially smaller than in the open-loop case. The performance achieved
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1945
Figure 4. Control result for the case a step disturbance of +10 °C in the inlet temperature.
Figure 6. Control result for the case a stochastic disturbance in the inlet temperature.
Figure 5. Control result for the case a step disturbance of -10 °C in the inlet temperature.
with the identified model is comparable to the performance obtained by Till et al.27 with a fundamental model. The control performance under a disturbance was also tested. First, a step disturbance was assumed to occur in the inlet temperature of the first feed stream. Figures 4 and 5 show the control performance under disturbances of +10 and -10 °C, respectively. Again, compared to the open-loop result, the controller does a good job of reducing the error, even though perfect control is not achieved because of the input constraints. Finally, we tested for the case of a stochastic disturbance occurring in the temperature of the first inlet stream. The control performance is illustrated in Figure 6. Notice that, here, the input trajectory, instead of converging to a periodic trajectory, continues to change in response to the stochastic disturbance. The open-loop result is not shown here, because the process ran away under constant-cooling-rate operation. 5.2. Example 2: Distillation Column. We next tested a binary distillation column of moderate purity with 20 trays. The usual assumptions of equimolar overflow and negligible vapor-phase dynamics were made in simulating the column. For details of the model equations, readers are referred to Appendix A of Morari and Zafiriou.28 The molar reflux flow rate L and boilup vapor flow rate V were treated as the manipulated
Figure 7. Profiles for the binary distillation column with periodic feed flow rate and constant L, V, and xf. (Sample time ) 2 min.)
variables (denoted by u1 and u2, respectively), and the compositions of the light component in the bottom and top products (denoted by y1 and y2, respectively) were regarded as the controlled variables. The composition of the light component in the feed was assumed to be the main disturbance, which was unmeasured. It was shown in ref 29 that steady-state operation of a distillation column, in general, might not be optimal and that, by periodically forcing the feed, product, reflux, and vapor flow rates, substantial energy savings can be achieved. The feed flow rate and composition might also fluctuate in a periodic manner as a result of some natural cycling of an upstream unit. To simulate a periodic operation, the inlet flow rate F was assumed to switch between a high and a low value, while other input variables, such as the inlet concentration xf, reflux flow rate L, and boilup flow rate V, were kept constant. The period was assumed to be 20 min, and the sample time was set to 2 min (hence, N ) 10). Figure 7 shows the trajectories of the input and output variables under these conditions, which exhibit distinct periodic behavior. The identification data were generated by simulating the nonlinear column model after introducing random perturbations to L and V about their nominal values,
1946
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003
Figure 8. Control result for the distillation column when the product compositions are directly measured and the initial input profiles are constant. No disturbance was assumed.
as well as to xf around its nominal value. The perturbations for L and V were drawn from independent uniform distributions in the range (-0.5, 0.5). The perturbations for xf were drawn from a uniform distribution in the range (-0.03, 0.03). In total, data for 200 periods were generated for identification, and additional data for 10 periods were generated for testing. All of the data used for identification were expressed in terms of deviations from the respective nominal trajectory values. Then, the lifted model was obtained by using the modified N4SID method. It was found that a three-state model was sufficient to capture most of the data trends. The lifted model was then converted into the augmented PTV model, on the basis of which a PTV Kalman filter and an LQR were designed. Here, we used the identified matrix Γ ˜ 2 for K. Figure 8 shows the on-line control performance of the controller for 30 runs when the inputs are initialized as constant trajectories. The controller continuously change the input profiles, and after several runs, zero tracking errors are achieved. Meanwhile, a step disturbance in the feed composition is tested, and zero tracking errors are also achieved. We also tested the case of a stochastic disturbance. A stochastic disturbance in xf, with the same statistics as that used for generating the identification data, was assumed. The on-line control performance is shown in Figure 9. Note that, in response to the stochastic disturbance, the control input trajectory is shown to be adjusted constantly around some normal trajectory, rather than converging to a fixed trajectory. As seen from Table 1, the error is substantially reduced compared to the open-loop case. Finally, we tested the efficacy of inferential control in this example. Here, two midtray temperatures (of trays 6 and 14) were assumed to be measured, whereas the controlled variables, the top and bottom compositions, were not measured. A lifted model was identified on the basis of the data from all four output variables. Whereas the LQR is designed for the two controlled variables only, the PTV Kalman filter was designed for the full output model to infer the product compositions from the midtray temperature measurements. Figure 10 shows the control result for the same disturbance as used in Figure 9. Because of some estimation error, the
Figure 9. Control result for the distillation column when the product compositions are directly measured and a stochastic disturbance in the feed composition is assumed.
Figure 10. Control result for the distillation column when the product compositions are inferred from the midtray temperatures and a stochastic disturbance in the feed composition is assumed. Table 1. Sums of Squared Errors for the Case of a Stochastic Disturbance SSE
open-loop
product composition measured directly
inferential control
y1 y2
0.7424 0.3097
0.0914 0.0610
0.1892 0.1188
tracking errors are a bit larger than those for the case when the product compositions are measured. Nevertheless, satisfactory control performance is achieved. The sums of squared errors of the product compositions for different controller designs are shown in Table 1. Conclusion Periodic disturbances and periodic operation of a process present us with an interesting control challenge because of the strong nonlinear and nonstationary characteristics. Theories for optimal control of periodic systems are already well in place, but the bottleneck lies in obtaining an accurate PTV system model in a reliable, efficient manner. The proposed route of identifying a lifted system model and rewriting it as a PTV system model represents a simple and generally applicable alternative to the approach of identifying a PTV
Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1947
system model directly. The form of the resulting model also enables us to conveniently design a controller that has integral action over the period index and therefore gives good converged behavior. Its application was illustrated in simulated studies involving a polyamine reactor and a binary distillation column. Acknowledgment The authors gratefully acknowledge the financial support from the National Science Foundation (CTS#0096326). Thanks also go to Jochen Till, Andreas Curse, and Wolfgang Marquardt for providing the simulation code of the polyamine reactor. Literature Cited (1) Butler, D.; Friedler, E.; Gatt, K. Characterizing the quantity and quality of domestic wastewater inflows. Water Sci. Technol. 1995, 31 (7), 13-24. (2) Bailey, J. E. Periodic operation of chemical reactors: A review. Chem. Eng. Commun. 1973, 1, 111-124. (3) Ruthven, D. M.; Farooq, S.; Knaebel, K. S. Pressure Swing Adsorption; VCH Publishers: New York, 1994. (4) Wu, D. J.; Xie, Y.; Ma, Z.; Wang, N.-H. L. Design of SMB chromatography for amino acid separations. Ind. Eng. Chem. Res. 1998, 37, 4023-4035. (5) Bittanti, S.; Coaneri, P.; Nicolao, G. D. The difference periodic Riccati equation for the periodic prediction problem. IEEE Trans. Autom. Control 1988, AC-33, 706-712. (6) Grasselli, O. M.; Longhi, S. Robust tracking and regulation of linear periodic discrete-time systems. Int. J. Control 1991, 54, 613-633. (7) Mayer, R. A.; Burrus, C. S. A unified analysis of multirate and periodically time-varying digital filters. IEEE Trans. Circuits Syst. 1975, 22, 162-168. (8) Bittanti, S.; Bolzern, P.; Nicolao, G.; Piroddi, L. Representation, prediction, and identification of cyclostationary processess A state-space approach. In Cyclostationarity in Communications and Signal Processing; Bardner, W. A., Ed.; IEEE Press: Piscataway, NJ, 1994; pp 267-294. (9) Verhaegen, M.; Yu, X. A class of subspace model identification algorithms to identify periodically and arbitrarily time-varying systems. Automatica 1995, 31 (2), 201-216. (10) Li, D.; Shah, S. L.; Chen, T. Identification of fast-rate models from multirate data. Int. J. Control 2001, 74, 680-689. (11) Lin, C.-A.; King, C.-W. Minimal periodic realizations of transfer matrices. IEEE Trans. Autom. Control 1993, 38, 462266. (12) Sanchez, E. V.; Hernandez, V.; Bru, R. Minimal realization for discrete-time periodic systems. Linear Algebra Appl. 1992, 162-164, 685-708. (13) Colaneri, P.; Longhi, S. The realization problem for linear periodic systems. Automatica 1995, 31 (5), 775-779. (14) Bittanti, S. Deterministic and stochastic linear periodic systems. In Time Series and Linear Systems; Bittanti, S., Ed.; Springer-Verlag: Berlin, 1986; pp 141-182. (15) Van Overschee, P.; De Moor, B. Subspace Identification for Linear Systems, Theory-Implementation-Applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (16) Hara, S.; Yamamoto, Y.; Omata, T.; Nakano, N. Repetitive control system. IEEE Trans. Autom. Control 1988, 33, 659-668. (17) Bittanti, S.; Guardabassi, G. Optimal cyclostationary control: The LQG approach. In IEEE Conference on Decision and Control; IEEE Press: Piscataway, NJ, 1981; pp 166-167.
(18) Lee, J. H.; Natarajan, S.; Lee, K. S. A model-based predictive control approach to repetitive control of continuous processes with periodic operations. J. Process Control 2001, 11, 195-207. (19) Amit, N. Optimal control of multirate digital control systems. Ph.D. Thesis, Stanford University, Stanford, CA, 1980. (20) Arnold, W. F., III; Laub, A. J. Generalized eigenproblem algorithms and software for algebraic Riccati equations. Proc. IEEE, 1984, 72, 1746-1754. (21) Lee, J. H.; Cooley, B. L. Recent advances in model predictive control and other related areas. AIChE Symp. Ser. 1997, 93, 201-216. (22) Natarajan, S.; Lee, J. H. Repetitive model predictive control applied to a simulate moving bed chromatography system. Comput. Chem. Eng. 2000, 24, 1127-1133. (23) Amrithalingam, R.; Su, W. S.; Lee, J. H. Two-step procedure for databased modeling for inferential control applications. AIChE J. 2000, 46 (10), 1974-1988. (24) Russel, J. A. Model of a Stirred Tank Reactor for Polyamine Synthesis; Publication ZT-TE7; Bayer-AG: Leverkusen, Germany, 1998. (25) Till, J. Repetitive model based predictive control of continuous large scale process with periodic operations. Master Thesis, Rheinisch-Westfa¨lischen Hochschule Aachen (RWTH Aachen), Aachen, Germany, 2001. (26) Cruse, A.; Marquardt, M.; Allgor, R. J.; Kussi, J. Integrated conceptual design of stirred tank reactors by periodic dynamic optimization. Comput. Chem. Eng. 2000, 24. (27) Till, J. E.; Lee, J. H.; Cruse A.; Marquardt, W. Repetitive model predictive control of a continuous polyamine process with periodic operations. Presented at the AIChE Annual Conference, Reno, NV, Nov 2001. (28) Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall, Englewood Cliffs, NJ, 1989. (29) Bausa, J.; Tsatsaronis, G. Reducing the energy demand of continuous distillation processes by optimal controlled forced periodic operation. Comput. Chem. Eng. 2001, 25, 359-370. (30) Bolzen, P.; Colaneri, P.; Scattolini, R. Zeros of discretetime linear periodic systems. IEEE Trans. Autom. Control 1986, AC-31, 1057-1058. (31) Colaneri, P. Zero-error regulation for discrete-time linear periodic systems. Syst. Control Lett. 1990, 15, 377-380. (32) Dahleh, M. A.; Voulgaris, P. G.; Valavani, L. S. Optimal and robust controllers for periodic and multirate systems. IEEE Trans. Autom. Control 1992, AC-37, 90-99. (33) Grasselli, O. M.; Longhi, S. Zeros and poles of linear periodic discrete-time systems. Circuit Syst. Singal Process. 1988, 7, 361-380. (34) Khargonekar, P. P.; Poolla, K.; Tannenbaum, A. Robust control of linear time-invariant plant using periodic compensators. IEEE Trans. Autom. Control 1985, 30, 1088-1096. (35) Kuijper, M. A periodically time-varying minimal partial realization algorithm. Automatica 1999, 35, 1543-1548. (36) Sterman, L. E.; Ydstie, B. E. The steady-state process with periodic perturbations. Chem. Eng. Sci 1990, 45 (3), 721-736. (37) Verhaegen, M.; Devilde, P. Subspace model identification. Part I: The output-error state space model identification algorithm. Int. J. Control 1992, 56, 1197-1210. (38) Verhaegen, M. Identification of the deterministic part of MIMO state space models given in innovations form from inputoutput data. Automatica 1994, 30 (1), 61-74.
Received for review April 25, 2002 Revised manuscript received January 13, 2003 Accepted February 3, 2003 IE020313Y
