- Email: [email protected]
LQR control of an infinite dimensional time-varying CSTR-PFR system*
Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011
LQR control of an infinite dimensional time-varying CSTR-PFR system ⋆ Amir Alizadeh Moghadam, Ilyasse Aksikas, Stevan Dubljevic, J. Fraser Forbes Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB Canada T6G 2V4 (e-mails: alizadeh, aksikas, Stevan.Dubljevic, [email protected]) Abstract: This contribution addresses the development of a Linear Quadratic Regulator (LQR) for a set of time-varying hyperbolic PDEs coupled with a set of time-varying ODEs through the boundary. The approach is based on an infinite-dimensional Hilbert state-space description of the system and the well-known operator Riccati equation (ORE). In order to solve the optimal control problem, the ORE is converted to a set of equivalent differential and algebraic matrix Riccati equations. The feedback operator can then be found by solving the resulting matrix Riccati equations. The performance of the designed control policy is assessed by applying it to a system of interconnecting continuous stirred tank reactor (CSTR) and a plug flow reactor (PFR) through a numerical simulation. Keywords: Distributed parameter system (DPS), Lumped parameter system(LPS), LQR control, Infinite dimensional system, Boundary control system 1. INTRODUCTION Many chemical processes have a distributed nature. A common approach for modelling such distributed processes involves the use of PDEs. Lumping assumptions are often used to convert the PDEs to sets of ODEs, which allows the use of a rich literature. However, this approximation results in some mismatches in the dynamical properties of the original distributed parameter and the lumped parameter models, which affects the performance of the designed model-based controller. A more rigorous way to deal with distributed parameter processes is to respect the infinite-dimensional characteristic of the system. This approach has been used by some researchers to design high-performance controllers for PDEs systems (see Curtain and Zwart (1995), Bensoussan et al. (2007), Orlov and Utkin (1987), Christofides (2001), and Krstic and smyshlyaev (2008)). Occasionally, distributed chemical processes are coupled with lumped parameter processes. Such systems are modelled by a combination of PDEs and ODEs (DPS-LPS) in which the interaction can appear either in the differential equations or in the boundary conditions. For instance, in pressure swing adsorption, the mass balance on components are modelled by a set of PDEs and the adsorption rates are modelled by a set of ODEs. Another example includes a jacket-equipped fixed-bed reactor, where the reactor (DPS) is interacting with the well-mixed jacket (LPS) via its boundary (see Wang (1966), Tzafestas (1970a), and Oh (1995) for more examples). Composite PDE-ODE models can also be used to describe transportation delay process accompanied by other chemical processes such as chemical reactions or mixing where differential-difference ⋆ This work is supported by Natural Science and Engineering Research Council (NSERC).
978-3-902661-93-7/11/$20.00 © 2011 IFAC
equations fail to model the system (see Hiratsuka and Ichikawa (1969)). In such systems, the transportation delay process can be modelled by hyperbolic PDEs and the other portion of the system can be described by ODEs. Despite the importance and inherent complexities in the structure of composite lumped and distributed parameter systems, research in the area of feedback control for these systems is scarce. In Wang (1966) a sufficient condition for stability and asymptotic stability for a scalar parabolic PDE-ODE system was derived by using the maximum principle for parabolic partial differential equations. However, this work assumed a specific type of the elliptical operator to use the weak maximum principle; moreover, the coupling terms were assumed to be uniformly bounded. In Tzafestas (1970b) classical Calculus of Variations was used to solve the optimal control problem for non-linear, mixed lumped and distributed parameter systems and the associated canonical equations were derived. In addition, the optimal final-value control problem for these systems were solved in Tzafestas (1970a) and necessary optimality conditions were found by applying Green’s identity together with functional analysis techniques. In this work, the optimal control input is calculated by solving a set of non-linear canonical equations which involves a large number of coupled PDEs and ODEs. Dynamic programming was used in Thowsen and Perkins (1973) and Thowsen and Perkins (1975) to solve discrete-time optimal feedback control problem for a class of linear composite systems. A classical method in the optimal feedback controller synthesis is the well known Linear Quadratic Regulator (LQR). The main objective of this control policy is to regulate a linear system by minimizing a quadratic performance index. An important advantage of LQR control is using a state feedback law, in which the state feedback gain is independent of the system’s states. Thus, the gain is
14446
10.3182/20110828-6-IT-1002.02195
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
calculated just once which reduces the amount of on-line calculations, significantly. In solving a LQR problem for an infinite-dimensional (distributed) system, two common methods are available in the literature. The first approach is based on frequency domain description and is known as spectral factorization. In this method the control law is obtained via solving an operator Diophantine equation (Callier and Winkin (1990)). This technique is applied in Aksikas et al. (2007) to control the temperature and the concentration in a plug flow reactor. The second methodology involves solving an ORE for a given state-space model (Curtain and Zwart (1995)). This method was used in Aksikas et al. (2008) for a particular class of hyperbolic PDEs. The methodology was then extended to a more general class of hyperbolic system by using an infinitedimensional Hilbert state-space setting with distributed input and output(Aksikas et al. (2009)). When a state-space model is available, solving the optimal control problem with the ORE method requires less computational effort in comparison to the spectral factorization approach which is more convenient for transfer function models. In our previous work (Moghadam et al. (2010)), the LQR control problem for a class of linear hyperbolic distributed parameter system interacting with a set of linear lumped parameter system through a Dirichlet boundary condition was treated. In such system, control variable affects the boundary of the distributed parameter system directly or indirectly, and through the lumped parameter system. In the present work, the methodology is extended to a timevarying system which is common in chemical operations, due to the phenomena such as catalyst deactivation and heat-exchanger fouling. In addition, dynamical properties of the system including stability and stabilizability are explored. In order to solve the optimal control problem, first, the system is described as an infinite-dimensional statespace by using boundary control transformation method. The infinite-time horizon LQR control problem for the system is then formulated and the ORE is computed. For solving the ORE, it is converted to equivalent matrix Riccati equation and an iterative computational algorithm is proposed for solving the resulting matrix Riccati equation. Finally, an example is given in order to illustrate the theory. 2. FORMULATION OF THE PROBLEM The general mathematical formulation for the systems considered here is as follows: ∂xd ∂xd (t, z) = V (t, z) + M (t, z)xd (t, z) + Bd (t, z)u(t) ∂t ∂z (1) dxl (t) = A(t)xl (t) + B(t)u(t) (2) dt y(t, .) = C(t, .)[xd (t, .), xl (t)]T (3) with the following boundary and initial conditions: xd (t, 0) = xl (t) (4) xd (0, z) = xd,0 (z) (5) xl (0) = xl,0 (6) where xd (t, .) ∈ L2 (0, 1)n and xl (t) ∈ Rn denote the state variables for the distributed and the lumped parameter systems, respectively; y(t) = [yd (t, .), yl (t)]T ∈ L2 (0, 1)p ⊕ Rp ; yd (t, .) ∈ Y := L2 (0, 1)p is the output variable
for the distributed system; yl (t) ∈ Rp is the output variable for the lumped system; z ∈ [0, 1] is the spatial coordinate; t ∈ [0, ∞] is the time; u(t) ∈ Rm is the input variable; V = −υI ∈ Rn×n with υ > 0 is a symmetric matrix; M (t, z) and Bd (t, z) are continuous matrices whose entries are functions in L∞ ((0, 1) × (0, ∞)); C(t, .) = diag(Cd (t, .), Cl (t)); Cd (t, .) is a continuous matrix whose entries are functions in L∞ ((0, 1) × (0, ∞)); A(t), B(t) and Cl (t) are real matrices with entries in L∞ (0, ∞); xd,0 (z) is a real continuous space-varying vector; and xl,0 is a constant vector. System (1) to (6) can be stated as an infinite-dimensional state-space system in the Hilbert space H = L2 (0, 1)n (Curtain and Zwart (1995)): x˙ d (t) = A (t)xd (t) + Bd (t)u(t) (7) x˙ l (t) = A(t)xl (t) + B(t)u(t) (8) y(t) = C (t)[xd (t), xl (t)]T (9) Bb xd (t) = xl (t) (10) Here A (t) is a linear operator defined as: dh + M (t, z)h (11) A (t)h(z) = V dz with the following domain: dh(z) D(A (t)) = {h(z) ∈ H : h(z) and dz (12) dh(z) are a.c. cont., ∈H} dz where a.c. means absolutely continuous. Bb ∈ L(H , Rn ) is a linear boundary operator defined as: Bb h(z) = h(0) (13) D(Bb ) = {h(z) ∈ H : h(z) is a.c.} (14) Bd (t) ∈ L(Rm , H ) is given by Bd (t) = Bd (t, .)I and C (t) ∈ L(H ⊕ Rn , Y ⊕ Rp ) is given by C (t) = C(t, .)I, where I is the identity operator. The infinite-dimensional state-space system (7) to (10) with an inhomogeneous boundary condition can be transformed to a new system with a homogenous boundary condition using boundary control transformation (see Curtain and Zwart (1995) and Fattorini (1968)). We assume that there is a function B(z) such that for all xl (t), Bxl (t) ∈ D(A (t)) and: Bb Bxl (t) = xl (t) (15) n By assuming that xl (t) ∈ L2 (0, ∞) is sufficiently smooth and using the state transformation ω(t) = xd (t) − Bxl (t) (Curtain and Zwart (1995)), we have: ω(t) ˙ = x˙ d (t) − Bx˙ l (t) Then: ω(t) ˙ = F (t)ω(t) + A (t)B(z)xl (t) + Bd (t)u(t) − B(z)x˙ l (t) ω(0) = ω0 (16) where ω0 = xd,0 − B(z)xl,0 ∈ D(F (t)) and: F (t)h(z) = A (t)h(z) The domain of F (t) is defined as: D(F (t)) = D(A (t)) ∩ ker(Bb ) = {h(z) ∈ H : dh(z) dh(z) are a.c., ∈ H , (17) h(z) and dz dz and h(0) = 0} By combining (8) and (16) we obtain the new infinitedimensional state-space representation of the DPS-LPS on H ⊕ Rn as:
14447
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
ω(t) ˙ F (t) A(t) ω(t) B¯d (t) = + u(t) x˙ l (t) 0 A(t) xl (t) B(t)
y(t) = C(t)[ω(t), xl (t)]T ω(0) = ω0 , xl (0) = xl,0
(18)
where A(t) = A (t)B(z) − B(z)A(t); B¯d (t) = Bd (t) − I B(z) m B(z)B(t) ∈ L(R , H ) and C(t) = C (t) ∈ 0 I L(H ⊕ Rn , Y ⊕ Rp ) are linear operators. Remark 1. System (18) has a homogeneous boundary condition as Bb ω(t) = 0. We define the state variables for system (18) as x(t) = [ω(t), xl (t)]T ∈ H ⊕ Rn . Now, the system can be written as: x(t) ˙ = A(t)x(t) + B(t)u(t) y(t) = C(t)x(t) (19) x0 = [ω0 , xl,0 ]T ∈ H ⊕ Rn F (t) A(t) B¯d (t) where A(t) = and B(t) = are linear 0 A(t) B(t) operators. Theorem 2. if {F (t)}t∈[0,T ] is an exponentially stable family of infinitesimal generators of C0 −semigroup on H and ϕ(t, s) = eA(t−s) is exponentially stable on Rn , {A(t)}t∈[0,T ] will be an exponentially stable family of infinitesimal generators of C0 −semigroup on H ⊕ Rn . Proof. The theory can be proved by using the mild solution of system (19) according to Theorem 5.3, p. 147 in Pazy (1983). 2 3. OPTIMAL CONTROL DESIGN In this section we are interested in LQR control synthesis for the DPS-LPS system according to the infinitedimensional state-space representation (19). The design is based on the minimization of an infinite-time horizon quadratic objective function that requires the solution of an ORE (Curtain and Zwart (1995),Bensoussan et al. (2007)). The solution of the ORE can be achieved by converting it to equivalent matrix Riccati equation. The optimal feedback operator can then be found by solving the resulting matrix Riccati equation. Let us consider the following infinite-time horizon quadratic objective function: Z ∞ (hCx(t), PCx(t)i + hu(t), Ru(t)i)dt (20) J(x0 , u) = 0 P P where P = P I ∈ L(Y ⊕ Rp ); P = 11 12 ∈ R2p×2p is P21 P22 a positive semi-definite symmetric matrix; and R ∈ Rm×m is a positive symmetric matrix. The minimization of the above objective function subject to system (19) results in solving the following ORE (Bensoussan et al. (2007)): ∗ ˙ [Q(t) + A(t) Q(t) + Q(t)A(t) + C(t)∗ PC(t)− (21) Q(t)B(t)R−1 B(t)∗ Q(t)]x = 0 According to Theorem 5.2, p. 507 in Bensoussan et al. (2007), if A(t) and B(t) meet conditions (7.8) and (7.9), p. 415 in Bensoussan et al. (2007) and pair (A(t), B(t)) is C(t)−stabilizable, the above ORE has a non-negative bounded solution Q(t) and this solution is minimal among
all non-negative bounded solutions. In addition, for any initial condition x0 ∈ H ⊕ Rn the unique optimal control variable uopt is obtained on t ≥ 0 as: uopt (t) = K(t)x(t) (22) where K(t) = −R−1 B(t)∗ Q(t) (23) and the minimum cost function is given by J(x0 , uopt ) = hx0 , Q(0)x0 i. Moreover, A(t) + B(t)K(t) is an exponentially stable family of infinitesimal generators of C0 −semigroup on H ⊕ Rn . Remark 3. Following from Theorem 2, A(t) is a family of infinitesimal generators of C0 −semigroup on H ⊕ Rn , therefore, it satisfies condition (7.8), p. 415 in Bensoussan et al. (2007). On the other hand since the entries of operators B(t) and C(t) are continuous, condition (7.9) is also satisfied. Theorem 4. If the finite-dimensional system (A(t), B(t)) is exponentially stabilizable on Rn , so will be system (19) on H ⊕ Rn . Proof. The proof is based on the closed-loop form of system (19) under feedback control law which stabilizes the finite-dimensional system. 2 Remark 5. Regarding definition 2.1, p. 480 in Bensoussan et al. (2007), it can be concluded that (A(t), B(t)) will be C(t)−stabilizable for any bounded operator C(t). It can be concluded from Theorem 4 and Remarks 3 and 5 that if the finite-dimensional system is exponentially stabilizable, ORE (21) has a non-negative and self-adjoint solution Q(t) ∈ L(H ⊕ Rn ). In order to solve the ORE we are going to convert it to equivalent matrix Riccati equation. Let the solution to the ORE (21) be: Φ0 (t, .)I 0 Q(t) := (24) 0 Ψ0 (t)I where Φ0 (t, .) and Ψ0 (t) ∈ Rn×n are positive self-adjoint matrices. By substituting for A(t), B(t), and C(t) into (21), we obtain the following set of equations: ˙ Φ(t) + F ∗ (t)Φ(t)+Φ(t)F (t) + Cd (t)∗ P11 Cd (t)− (25) Φ(t)B¯d (t)R−1 B¯d (t)∗ Φ(t) = 0 Φ(t)A(t) + Cd (t)∗ P11 Cd (t)B + Cd (t)∗ P12 Cl (t)− Φ(t)B¯d (t)(t)R−1 B(t)∗ Ψ(t) = 0
(26)
˙ Ψ(t) + A∗ Ψ(t) + Ψ(t)A(t) + B∗ Cd (t)∗ P11 Cd (t)B+ Cl∗ P21 Cd (t)B + B∗ Cd (t)∗ P12 Cl + Cl∗ P22 Cl − (27) Ψ(t)B(t)R−1 B(t)∗ Ψ(t) = 0 where Φ(t) = Φ0 (t, .)I; Ψ(t) = Ψ0 (t)I; Cd (t) = Cd (t, .)I; and Cl (t) = Cl (t)I; P11 = P11 I; P12 = P12 I; P21 = P21 I; and P22 = P22 I. Equation (25) is an ORE. We assume that the matrix V = −υI, υ > 0 is diagonal with identical entries. In this condition, (25) can be converted to the following differential matrix Riccati equation (see Aksikas and Forbes (2010) for the proof): ∂Φ0 ∂Φ0 = −V + M ∗ Φ0 + Φ0 M + Cd∗ P11 Cd − ∂t ∂z ¯d R−1 B ¯d∗ Φ0 (28) Φ0 B Φ0 (t, 1) = 0
14448
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
¯d = Bd − BB. where B The following algebraic matrix Riccati equation can be obtained from (26): Φ0 M B − Φ0 BA + Cd∗ P11 Cd B + Cd∗ P12 Cl − (29) ¯d R−1 B ∗ Ψ0 = 0 Φ0 B The equivalent matrix Riccati equation for (27) is: ˙ 0 + A∗ Ψ0 + Ψ0 A + B∗ C ∗ P11 Cd B + C ∗ P21 Cd B+ Ψ d
l
B∗ Cd∗ P12 Cl + Cl∗ P22 Cl − Ψ0 BR−1 B ∗ Ψ0 = 0 (30) Equations (28), (29), and (30) form a set of partial differential, ordinary differential and algebraic equations (PDAEs) in which Φ0 (t, z) and Ψ0 (t) are the differential variables and P12 is the algebraic variable. This set of equations can be solved by using the following procedure: • Choose weighting matrices P11 , P22 , and R; • Find Ψ0 (t), Φ0 (t, z) and P12 via solving PDAEs system (28) to (30); • Check whether matrix P is positive semi-definite; • In the case that P is not positive semi-definite choose new weighting matrices and re-solve the PDAEs system; ∗ . Note that since matrix P is symmetric, we have P12 = P21 Also of note is that, based on our experience, P can be ensured to be positive semi-definite by always selecting P22 ≫ P11 . Finally, the state feedback operator can be calculated from (23) which will be: (31) K(t) = −R−1 B¯d (t)∗ Φ B(t)∗ Ψ
Fig. 1. CSTR-PFR system and indirectly, and through the CSTR. With the assumptions of negligible diffusion in the PFR, perfect level control in the CSTR, no transportation lags in the connecting lines, constant fluid velocity in the PFR with respect to spatial coordinate, and constant physical properties, and incompressible fluid, the mathematical model of the system will be: Fin in dCA = (CA − CA ) − k1 e−E1 /RT CA − dt Vs 2 k3 e−E3 /RT CA (34) Fin dCB =− CB + k1 e−E1 /RT CA − dt Vs k2 e−E2 /RT CB (35) 1 dT = [k1 e−E1 /RT CA (−∆H1 )+ dt ρcp 2 k2 e−E2 /RT CB (−∆H2 ) + k3 e−E3 /RT CA (−∆H3 )]+ Q Fin (Tin − T ) + (36) Vs ρcp Vs
4. CASE STUDY In this section we consider a CSTR-PFR configuration shown in Fig. 1 as an interacting lumped and hyperbolic distributed parameter system. This reactor configuration is recommended for some types of chemical reactions (see e.g., Fogler (2005) and Schweier and Floudas (1999)) and may be used to carry out Van de Vusse reaction to achieve the maximum conversion to the desired product (Schweier and Floudas (1999)). Here, we assume reactions and kinetics: A −→ B, −r1 = k1 e−E1 /RT CA B −→ C,
−r2 = k2 e−E2 /RT CB −E3 /RT
2 CA
(32)
2A −→ D, −r3 = k3 e where k1 , k2 and k3 are pre-exponential constants; E1 , E2 , and E3 are the activation energy; and R is the universal gas constant. The exothermic reactions take place in both CSTR and PFR and component B is the desired component. The reaction kinetics are considered to be time-varying according to the following exponential decay model (see Lie and Himmelblau (2000)): ki = k0,i + k1,i e−αi t (33) where subscripts i = 1, 2, 3 denote number of reactions and k0,i , k1,i , and αi are the decay model parameters. The objective is to control the concentration of the components and the temperature within both reactors by using inlet flow rate (Fin ) and cooling rate from the CSTR (Q) as manipulated variables. It should be noticed that when liquid level is perfectly controlled in the CSTR, variations of Fin affects the controlled variables in the PFR directly
p ∂CA ∂C p p = −υ A − k1 e−E1 /RT CA − ∂t ∂z p2 (37) k3 e−E3 /RT CA p p ∂CB ∂CB p = −υ + k1 e−E1 /RT CA − ∂t ∂z p k2 e−E2 /RT CB (38) ∂Tp ∂Tp k1 −E1 /RT p e CA (−∆H1 )+ = −υ + ∂t ∂z ρcp k2 −E2 /RT p e CB (−∆H2 )+ (39) ρcp k3 −E3 /RT p 2 e CA (−∆H3 ) ρcp p CA (t, 0) = CA (40) p (41) CB (t, 0) = CB Tp (t, 0) = T (42) where CA and CB are the concentration of the components A and B in the CSTR, respectively; T is the temperature p p in the CSTR; CA and CB are the concentration of the components A and B in the PFR, respectively; Tp is the temperature in the PFR; z ∈ [0, L] is the spatial in coordinate; t ∈ [0, ∞] is the time; Fin , CA , and Tin are the volumetric flow-rate, concentration, and temperature of the feed to the CSTR; Vs and Vp are the volumes of the CSTR and the PFR, respectively; υ is the fluid velocity in the PFR which is given by υ = FVinpL ; ∆H1 , ∆H2 , and ∆H3 are the heat of reaction for reactions 1, 2, and 3, respectively; and ρ and cp are the average fluid density
14449
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
Table 1. Model parameters Concentration (kmol/m3)
Value 150.150 × 106 sec−1 225.225 × 106 sec−1 1.759 × 106 sec−1 100.1 × 106 sec−1 25.025 × 106 sec−1 0 sec−1 1.389 × 10−2 sec−1 174.845 × 10−6 m3 /sec −1.36 kJ/sec 5.1 kmol/m3 403.15 K −4200 kJ/kmol −11000 kJ/kmol −41850 kJ/kmol 9758.3 K 9758.3 K 8560.0 K 0.01 m3 0.005 m3 934.2 kg/m3 3.01 kJ/kgK
3
419
2.5
418 417
Tp
416
1.5
415
1 0.5 0
414 0.05
0.1
0.15 z (m)
Fig. 2. Steady-state profiles in the PFR
0.2
0.25
413 0.3
Temperature (K)
Concentration (kmol/m3)
We use the proposed optimal policy to control the concentration of the components and also the temperature in both CSTR and PFR. Cd (t, z) and Cl (t) are selected to be the identity matrices to yield a state LQR problem. In order to solve the LQR control problem for this system, we follow the procedure discussed in section 3. The set of differential and algebraic equations (28) to (30) are solved numerically in gPROMS. By choosing P11 = diag(80, 80), P22 = diag(350, 350), and R = diag(0.01, 0.0001) the matrix P will be positive semi-definite. The feedback operator can then be obtained using (31).
CpB
A
C (open−loop)
8
B
C (closed−loop) B
6 4 2 100
200 300 Time (sec)
400
500
520
4.1 Controller Design
CpA
A
C (closed−loop)
Fig. 3. Concentration responses in the CSTR
and specific heat. The model parameters used in this case study are given in table 1. In order to find the equilibrium condition for the system, the modelling equations (34) to (42) have R been solved at steady-state in gPROMS (gpr (19972010)). Simulation yields the steady-state values for CA , CB , and T ; 2.89 kmol/m3 , 0.80 kmol/m3 , and 413.10 K, p p respectively. The steady-state profiles for CA , CB and Tp are shown in Fig. 2. The model equations can be linearized around the steady-state condition to describe the system as the general formulation (1) to (6).
2
C (open−loop)
10
0 0
Temperature (K)
Parameter k0,1 k0,2 k0,3 k1,1 k1,2 k1,3 α1 , α2 , α 3 Fin,ss Qss in CA Tin ∆H1 ∆H2 ∆H3 E1 /R E2 /R E3 /R Vs Vp ρ cp
12
500
T (open−loop) T (closed−loop) Tp (open−loop)
480
Tp (closed−loop)
460 440 420 400 0
100
200 300 Time (sec)
400
500
Fig. 4. Temperature responses in the CSTR and at the outlet of the PFR
4.2 Simulation Results In order to assess the performance of the control policy, the designed feedback operator is used with the original non-linear system (34) to (42). The set of coupled nonlinear PDEs and ODEs are solved using orthogonal collocation on finite element method in gPROMS. We use CA (0) = 12 kmol/m3 , CB (0) = 12 kmol/m3 , T (0) = 400 p p K, CA (0, z) = 0.7 kmol/m3 , CB (0, z) = 0.7 kmol/m3 , and Tp (0, z) = 400 K as the initial conditions. To have a measure of how good the designed controller is, the responses of open-loop and closed-loop systems for the concentrations and the temperature in the CSTR and at the outlet of the PFR are compared in Fig. 3 to Fig. 5. As it can be seen, the controlled system is able to reject the effect of the initial condition about 2 times faster (with respect to the residence time) than the open-loop system, and it converges to the desired steady-state condition, accurately. Moreover, the open-loop inverse-responses appearing at the outlet of the PFR have disappeared. The variations of the control inputs are also shown in Fig. 6. The control efforts are not particularly aggressive, and are physically realizable. 5. CONCLUSION In this work LQR control problem for a class of timevarying, composite lumped and distributed parameter system is formulated and solved. In this mixed system, the lumped parameter system interacts with the hyperbolic distributed parameter system through its boundary. The control variables can affect the boundary of the distributed parameter system directly or indirectly, and through the
14450
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
Concentration (kmol/m3)
3 p C A p C A p C B p C B
2.5 2 1.5
(open−loop) (closed−loop) (open−loop) (closed−loop)
1 0.5 0 0
100
200 300 Time (sec)
400
500
Fig. 5. Concentration responses at the outlet of the PFR −4
x 10
−1
4
Fin
3.5
Q
−1.5
−2
3
2.5
Q (kJ/sec)
Fin (m3/sec)
4.5
−2.5
2 1.5 0
100
200 300 Time (sec)
400
−3 500
Fig. 6. Control input profiles lumped parameter system. The LQR control problem is formulated based on an infinite-dimensional state-space representation of the system, which is obtained via a state transformation from the original system using boundary control transformation method. The solution of the LQR control problem is achieved by solving the matrix Riccati equation that results from the ORE of the infinitedimensional state-space representation. The designed optimal control policy was implemented on an interacting CSTR-PFR system. The performance of the controller was assessed by implementing the controller on the original non-linear system and resulted in a high performance closed-loop system. REFERENCES (1997-2010). Process systems enterprise, gPROMS, www.psenterprise.com/gproms. Acquistapace, P. and Terreni, B. (1984). Some existence and regularity results for abstract non-autonomous parabolic equations. Journal of Mathematical Analysis and Applications, 99(1), 9–64. Aksikas, I. and Forbes, J.F. (2010). Linear quadratic regulator of time-varying hyperbolic distributed parameter systems. IMA Journal of Mathematical Control and Information, 27(3), 387–401. Aksikas, I., Fuxman, A., Forbes, J.F., and Winkin, J. (2009). LQ control design of a class of hyperbolic PDE systems: Application to fixed-bed reactor. Automatica, 45, 1542–1548. Aksikas, I., Winkin, J., and Dochain, D. (2007). Optimal LQ-feedback regulation of a nonisothermal plug flow reactor model by spectral factorization. IEEE Transaction on Automatic Control, 52(7), 1197–1193.
Aksikas, I., Winkin, J., and Dochain, D. (2008). Optimal LQ-feedback for a class of first-order hyperbolic distributed parameter systems. ESAIM:COCV. Bensoussan, A., Prato, G.D., and Mitter, K. (2007). Representation and Control of Infinite Dimensional Systems. Systems and Control: Foundations & Applications. Birkh¨ auser Boston. Callier, F. and Winkin, J. (1990). Spectral factorization and LQ-optimal regulation for multivariable distributed systems. International Journal of Control, 52(1), 55–57. Christofides, P.D. (2001). Nonlinear and robust control of PDE systems: Methods and applications to transportreaction processes. Boston:Brikhauser. Curtain, R.F. and Zwart, H.J. (1995). An Introduction to Infinite Dimensional Linear Systems. Springel-Verlag. Fattorini, H.O. (1968). Boundary control systems. SIAM Journal on Control, 6, 349–385. Fogler, H.S. (2005). Elements of Chemical Reaction Engineering. Prentice-Hall. Hiratsuka, S. and Ichikawa, A. (1969). Optimal control of systems with transportation lags. IEEE Transaction on Automatic Control, 14(3), 237. Krstic, M. and smyshlyaev, A. (2008). Boundary control of PDEs: A course on backstepping designs. Society for industrial and Applied Mathematics, Philadelphia, PA. Lie, B. and Himmelblau, M. (2000). Catalyst deactivation: control relevance of model assumptions. Ind. Eng. Chem. Res, 39, 1242–1248. Moghadam, A.A., Aksikas, I., Dubljevic, S., and Forbes, J.F. (2010). LQ control of coupled hyperbolic PDEs and ODEs: Application to a CSTR-PFR system. Proceedings of the 9th International Symposium on Dynamics and Control of Process Systems, DYCOPS 2010, 713–718. Oh, M. (1995). Modelling and simulation of combined lumped and distributed processes. Ph.D. thesis, Imperial College, London. Orlov, Y. and Utkin, V. (1987). Sliding mode control in infinite-dimensional systems. Automatica, 23(753-757). Pazy, A. (1983). Semigroup of linear operators and application to partial differential equations. Springel-Verlag, New York. Schweier, C.A. and Floudas, C.A. (1999). Optimization framework for the synthesis of chemical reactor networks. Ind. Eng. Chem. Res, 38, 744–766. Thowsen, A. and Perkins, W.R. (1975). Sampled-data linear-quadratic regulator for systems with generalized transportation lags. International Journal of Control, 21, 49–65. Thowsen, A. and Perkins, W.R. (1973). Optimal discretetime feedback control of mixed distributed and lumped parameter systems. International Journal of Control, 18(3), 657–665. Tzafestas, S. (1970a). Final-value control of nonlinear composite distributed and lumped parameter systems. Franklin Inst., 290(5). Tzafestas, S. (1970b). Optimal distributed-parameter control using classical variational theory. International Journal of Control, 12(4), 593–608. Wang, P.K.C. (1966). On the stability of equilibrium of a mixed distributed and lumped parameter control system. International Journal of Control, 3(2), 139–147.
14451
LQR control of an infinite dimensional time-varying CSTR-PFR system ⋆ Amir Alizadeh Moghadam, Ilyasse Aksikas, Stevan Dubljevic, J. Fraser Forbes Department of Chemical and Materials Engineering, University of Alberta, Edmonton, AB Canada T6G 2V4 (e-mails: alizadeh, aksikas, Stevan.Dubljevic, [email protected]) Abstract: This contribution addresses the development of a Linear Quadratic Regulator (LQR) for a set of time-varying hyperbolic PDEs coupled with a set of time-varying ODEs through the boundary. The approach is based on an infinite-dimensional Hilbert state-space description of the system and the well-known operator Riccati equation (ORE). In order to solve the optimal control problem, the ORE is converted to a set of equivalent differential and algebraic matrix Riccati equations. The feedback operator can then be found by solving the resulting matrix Riccati equations. The performance of the designed control policy is assessed by applying it to a system of interconnecting continuous stirred tank reactor (CSTR) and a plug flow reactor (PFR) through a numerical simulation. Keywords: Distributed parameter system (DPS), Lumped parameter system(LPS), LQR control, Infinite dimensional system, Boundary control system 1. INTRODUCTION Many chemical processes have a distributed nature. A common approach for modelling such distributed processes involves the use of PDEs. Lumping assumptions are often used to convert the PDEs to sets of ODEs, which allows the use of a rich literature. However, this approximation results in some mismatches in the dynamical properties of the original distributed parameter and the lumped parameter models, which affects the performance of the designed model-based controller. A more rigorous way to deal with distributed parameter processes is to respect the infinite-dimensional characteristic of the system. This approach has been used by some researchers to design high-performance controllers for PDEs systems (see Curtain and Zwart (1995), Bensoussan et al. (2007), Orlov and Utkin (1987), Christofides (2001), and Krstic and smyshlyaev (2008)). Occasionally, distributed chemical processes are coupled with lumped parameter processes. Such systems are modelled by a combination of PDEs and ODEs (DPS-LPS) in which the interaction can appear either in the differential equations or in the boundary conditions. For instance, in pressure swing adsorption, the mass balance on components are modelled by a set of PDEs and the adsorption rates are modelled by a set of ODEs. Another example includes a jacket-equipped fixed-bed reactor, where the reactor (DPS) is interacting with the well-mixed jacket (LPS) via its boundary (see Wang (1966), Tzafestas (1970a), and Oh (1995) for more examples). Composite PDE-ODE models can also be used to describe transportation delay process accompanied by other chemical processes such as chemical reactions or mixing where differential-difference ⋆ This work is supported by Natural Science and Engineering Research Council (NSERC).
978-3-902661-93-7/11/$20.00 © 2011 IFAC
equations fail to model the system (see Hiratsuka and Ichikawa (1969)). In such systems, the transportation delay process can be modelled by hyperbolic PDEs and the other portion of the system can be described by ODEs. Despite the importance and inherent complexities in the structure of composite lumped and distributed parameter systems, research in the area of feedback control for these systems is scarce. In Wang (1966) a sufficient condition for stability and asymptotic stability for a scalar parabolic PDE-ODE system was derived by using the maximum principle for parabolic partial differential equations. However, this work assumed a specific type of the elliptical operator to use the weak maximum principle; moreover, the coupling terms were assumed to be uniformly bounded. In Tzafestas (1970b) classical Calculus of Variations was used to solve the optimal control problem for non-linear, mixed lumped and distributed parameter systems and the associated canonical equations were derived. In addition, the optimal final-value control problem for these systems were solved in Tzafestas (1970a) and necessary optimality conditions were found by applying Green’s identity together with functional analysis techniques. In this work, the optimal control input is calculated by solving a set of non-linear canonical equations which involves a large number of coupled PDEs and ODEs. Dynamic programming was used in Thowsen and Perkins (1973) and Thowsen and Perkins (1975) to solve discrete-time optimal feedback control problem for a class of linear composite systems. A classical method in the optimal feedback controller synthesis is the well known Linear Quadratic Regulator (LQR). The main objective of this control policy is to regulate a linear system by minimizing a quadratic performance index. An important advantage of LQR control is using a state feedback law, in which the state feedback gain is independent of the system’s states. Thus, the gain is
14446
10.3182/20110828-6-IT-1002.02195
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
calculated just once which reduces the amount of on-line calculations, significantly. In solving a LQR problem for an infinite-dimensional (distributed) system, two common methods are available in the literature. The first approach is based on frequency domain description and is known as spectral factorization. In this method the control law is obtained via solving an operator Diophantine equation (Callier and Winkin (1990)). This technique is applied in Aksikas et al. (2007) to control the temperature and the concentration in a plug flow reactor. The second methodology involves solving an ORE for a given state-space model (Curtain and Zwart (1995)). This method was used in Aksikas et al. (2008) for a particular class of hyperbolic PDEs. The methodology was then extended to a more general class of hyperbolic system by using an infinitedimensional Hilbert state-space setting with distributed input and output(Aksikas et al. (2009)). When a state-space model is available, solving the optimal control problem with the ORE method requires less computational effort in comparison to the spectral factorization approach which is more convenient for transfer function models. In our previous work (Moghadam et al. (2010)), the LQR control problem for a class of linear hyperbolic distributed parameter system interacting with a set of linear lumped parameter system through a Dirichlet boundary condition was treated. In such system, control variable affects the boundary of the distributed parameter system directly or indirectly, and through the lumped parameter system. In the present work, the methodology is extended to a timevarying system which is common in chemical operations, due to the phenomena such as catalyst deactivation and heat-exchanger fouling. In addition, dynamical properties of the system including stability and stabilizability are explored. In order to solve the optimal control problem, first, the system is described as an infinite-dimensional statespace by using boundary control transformation method. The infinite-time horizon LQR control problem for the system is then formulated and the ORE is computed. For solving the ORE, it is converted to equivalent matrix Riccati equation and an iterative computational algorithm is proposed for solving the resulting matrix Riccati equation. Finally, an example is given in order to illustrate the theory. 2. FORMULATION OF THE PROBLEM The general mathematical formulation for the systems considered here is as follows: ∂xd ∂xd (t, z) = V (t, z) + M (t, z)xd (t, z) + Bd (t, z)u(t) ∂t ∂z (1) dxl (t) = A(t)xl (t) + B(t)u(t) (2) dt y(t, .) = C(t, .)[xd (t, .), xl (t)]T (3) with the following boundary and initial conditions: xd (t, 0) = xl (t) (4) xd (0, z) = xd,0 (z) (5) xl (0) = xl,0 (6) where xd (t, .) ∈ L2 (0, 1)n and xl (t) ∈ Rn denote the state variables for the distributed and the lumped parameter systems, respectively; y(t) = [yd (t, .), yl (t)]T ∈ L2 (0, 1)p ⊕ Rp ; yd (t, .) ∈ Y := L2 (0, 1)p is the output variable
for the distributed system; yl (t) ∈ Rp is the output variable for the lumped system; z ∈ [0, 1] is the spatial coordinate; t ∈ [0, ∞] is the time; u(t) ∈ Rm is the input variable; V = −υI ∈ Rn×n with υ > 0 is a symmetric matrix; M (t, z) and Bd (t, z) are continuous matrices whose entries are functions in L∞ ((0, 1) × (0, ∞)); C(t, .) = diag(Cd (t, .), Cl (t)); Cd (t, .) is a continuous matrix whose entries are functions in L∞ ((0, 1) × (0, ∞)); A(t), B(t) and Cl (t) are real matrices with entries in L∞ (0, ∞); xd,0 (z) is a real continuous space-varying vector; and xl,0 is a constant vector. System (1) to (6) can be stated as an infinite-dimensional state-space system in the Hilbert space H = L2 (0, 1)n (Curtain and Zwart (1995)): x˙ d (t) = A (t)xd (t) + Bd (t)u(t) (7) x˙ l (t) = A(t)xl (t) + B(t)u(t) (8) y(t) = C (t)[xd (t), xl (t)]T (9) Bb xd (t) = xl (t) (10) Here A (t) is a linear operator defined as: dh + M (t, z)h (11) A (t)h(z) = V dz with the following domain: dh(z) D(A (t)) = {h(z) ∈ H : h(z) and dz (12) dh(z) are a.c. cont., ∈H} dz where a.c. means absolutely continuous. Bb ∈ L(H , Rn ) is a linear boundary operator defined as: Bb h(z) = h(0) (13) D(Bb ) = {h(z) ∈ H : h(z) is a.c.} (14) Bd (t) ∈ L(Rm , H ) is given by Bd (t) = Bd (t, .)I and C (t) ∈ L(H ⊕ Rn , Y ⊕ Rp ) is given by C (t) = C(t, .)I, where I is the identity operator. The infinite-dimensional state-space system (7) to (10) with an inhomogeneous boundary condition can be transformed to a new system with a homogenous boundary condition using boundary control transformation (see Curtain and Zwart (1995) and Fattorini (1968)). We assume that there is a function B(z) such that for all xl (t), Bxl (t) ∈ D(A (t)) and: Bb Bxl (t) = xl (t) (15) n By assuming that xl (t) ∈ L2 (0, ∞) is sufficiently smooth and using the state transformation ω(t) = xd (t) − Bxl (t) (Curtain and Zwart (1995)), we have: ω(t) ˙ = x˙ d (t) − Bx˙ l (t) Then: ω(t) ˙ = F (t)ω(t) + A (t)B(z)xl (t) + Bd (t)u(t) − B(z)x˙ l (t) ω(0) = ω0 (16) where ω0 = xd,0 − B(z)xl,0 ∈ D(F (t)) and: F (t)h(z) = A (t)h(z) The domain of F (t) is defined as: D(F (t)) = D(A (t)) ∩ ker(Bb ) = {h(z) ∈ H : dh(z) dh(z) are a.c., ∈ H , (17) h(z) and dz dz and h(0) = 0} By combining (8) and (16) we obtain the new infinitedimensional state-space representation of the DPS-LPS on H ⊕ Rn as:
14447
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
ω(t) ˙ F (t) A(t) ω(t) B¯d (t) = + u(t) x˙ l (t) 0 A(t) xl (t) B(t)
y(t) = C(t)[ω(t), xl (t)]T ω(0) = ω0 , xl (0) = xl,0
(18)
where A(t) = A (t)B(z) − B(z)A(t); B¯d (t) = Bd (t) − I B(z) m B(z)B(t) ∈ L(R , H ) and C(t) = C (t) ∈ 0 I L(H ⊕ Rn , Y ⊕ Rp ) are linear operators. Remark 1. System (18) has a homogeneous boundary condition as Bb ω(t) = 0. We define the state variables for system (18) as x(t) = [ω(t), xl (t)]T ∈ H ⊕ Rn . Now, the system can be written as: x(t) ˙ = A(t)x(t) + B(t)u(t) y(t) = C(t)x(t) (19) x0 = [ω0 , xl,0 ]T ∈ H ⊕ Rn F (t) A(t) B¯d (t) where A(t) = and B(t) = are linear 0 A(t) B(t) operators. Theorem 2. if {F (t)}t∈[0,T ] is an exponentially stable family of infinitesimal generators of C0 −semigroup on H and ϕ(t, s) = eA(t−s) is exponentially stable on Rn , {A(t)}t∈[0,T ] will be an exponentially stable family of infinitesimal generators of C0 −semigroup on H ⊕ Rn . Proof. The theory can be proved by using the mild solution of system (19) according to Theorem 5.3, p. 147 in Pazy (1983). 2 3. OPTIMAL CONTROL DESIGN In this section we are interested in LQR control synthesis for the DPS-LPS system according to the infinitedimensional state-space representation (19). The design is based on the minimization of an infinite-time horizon quadratic objective function that requires the solution of an ORE (Curtain and Zwart (1995),Bensoussan et al. (2007)). The solution of the ORE can be achieved by converting it to equivalent matrix Riccati equation. The optimal feedback operator can then be found by solving the resulting matrix Riccati equation. Let us consider the following infinite-time horizon quadratic objective function: Z ∞ (hCx(t), PCx(t)i + hu(t), Ru(t)i)dt (20) J(x0 , u) = 0 P P where P = P I ∈ L(Y ⊕ Rp ); P = 11 12 ∈ R2p×2p is P21 P22 a positive semi-definite symmetric matrix; and R ∈ Rm×m is a positive symmetric matrix. The minimization of the above objective function subject to system (19) results in solving the following ORE (Bensoussan et al. (2007)): ∗ ˙ [Q(t) + A(t) Q(t) + Q(t)A(t) + C(t)∗ PC(t)− (21) Q(t)B(t)R−1 B(t)∗ Q(t)]x = 0 According to Theorem 5.2, p. 507 in Bensoussan et al. (2007), if A(t) and B(t) meet conditions (7.8) and (7.9), p. 415 in Bensoussan et al. (2007) and pair (A(t), B(t)) is C(t)−stabilizable, the above ORE has a non-negative bounded solution Q(t) and this solution is minimal among
all non-negative bounded solutions. In addition, for any initial condition x0 ∈ H ⊕ Rn the unique optimal control variable uopt is obtained on t ≥ 0 as: uopt (t) = K(t)x(t) (22) where K(t) = −R−1 B(t)∗ Q(t) (23) and the minimum cost function is given by J(x0 , uopt ) = hx0 , Q(0)x0 i. Moreover, A(t) + B(t)K(t) is an exponentially stable family of infinitesimal generators of C0 −semigroup on H ⊕ Rn . Remark 3. Following from Theorem 2, A(t) is a family of infinitesimal generators of C0 −semigroup on H ⊕ Rn , therefore, it satisfies condition (7.8), p. 415 in Bensoussan et al. (2007). On the other hand since the entries of operators B(t) and C(t) are continuous, condition (7.9) is also satisfied. Theorem 4. If the finite-dimensional system (A(t), B(t)) is exponentially stabilizable on Rn , so will be system (19) on H ⊕ Rn . Proof. The proof is based on the closed-loop form of system (19) under feedback control law which stabilizes the finite-dimensional system. 2 Remark 5. Regarding definition 2.1, p. 480 in Bensoussan et al. (2007), it can be concluded that (A(t), B(t)) will be C(t)−stabilizable for any bounded operator C(t). It can be concluded from Theorem 4 and Remarks 3 and 5 that if the finite-dimensional system is exponentially stabilizable, ORE (21) has a non-negative and self-adjoint solution Q(t) ∈ L(H ⊕ Rn ). In order to solve the ORE we are going to convert it to equivalent matrix Riccati equation. Let the solution to the ORE (21) be: Φ0 (t, .)I 0 Q(t) := (24) 0 Ψ0 (t)I where Φ0 (t, .) and Ψ0 (t) ∈ Rn×n are positive self-adjoint matrices. By substituting for A(t), B(t), and C(t) into (21), we obtain the following set of equations: ˙ Φ(t) + F ∗ (t)Φ(t)+Φ(t)F (t) + Cd (t)∗ P11 Cd (t)− (25) Φ(t)B¯d (t)R−1 B¯d (t)∗ Φ(t) = 0 Φ(t)A(t) + Cd (t)∗ P11 Cd (t)B + Cd (t)∗ P12 Cl (t)− Φ(t)B¯d (t)(t)R−1 B(t)∗ Ψ(t) = 0
(26)
˙ Ψ(t) + A∗ Ψ(t) + Ψ(t)A(t) + B∗ Cd (t)∗ P11 Cd (t)B+ Cl∗ P21 Cd (t)B + B∗ Cd (t)∗ P12 Cl + Cl∗ P22 Cl − (27) Ψ(t)B(t)R−1 B(t)∗ Ψ(t) = 0 where Φ(t) = Φ0 (t, .)I; Ψ(t) = Ψ0 (t)I; Cd (t) = Cd (t, .)I; and Cl (t) = Cl (t)I; P11 = P11 I; P12 = P12 I; P21 = P21 I; and P22 = P22 I. Equation (25) is an ORE. We assume that the matrix V = −υI, υ > 0 is diagonal with identical entries. In this condition, (25) can be converted to the following differential matrix Riccati equation (see Aksikas and Forbes (2010) for the proof): ∂Φ0 ∂Φ0 = −V + M ∗ Φ0 + Φ0 M + Cd∗ P11 Cd − ∂t ∂z ¯d R−1 B ¯d∗ Φ0 (28) Φ0 B Φ0 (t, 1) = 0
14448
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
¯d = Bd − BB. where B The following algebraic matrix Riccati equation can be obtained from (26): Φ0 M B − Φ0 BA + Cd∗ P11 Cd B + Cd∗ P12 Cl − (29) ¯d R−1 B ∗ Ψ0 = 0 Φ0 B The equivalent matrix Riccati equation for (27) is: ˙ 0 + A∗ Ψ0 + Ψ0 A + B∗ C ∗ P11 Cd B + C ∗ P21 Cd B+ Ψ d
l
B∗ Cd∗ P12 Cl + Cl∗ P22 Cl − Ψ0 BR−1 B ∗ Ψ0 = 0 (30) Equations (28), (29), and (30) form a set of partial differential, ordinary differential and algebraic equations (PDAEs) in which Φ0 (t, z) and Ψ0 (t) are the differential variables and P12 is the algebraic variable. This set of equations can be solved by using the following procedure: • Choose weighting matrices P11 , P22 , and R; • Find Ψ0 (t), Φ0 (t, z) and P12 via solving PDAEs system (28) to (30); • Check whether matrix P is positive semi-definite; • In the case that P is not positive semi-definite choose new weighting matrices and re-solve the PDAEs system; ∗ . Note that since matrix P is symmetric, we have P12 = P21 Also of note is that, based on our experience, P can be ensured to be positive semi-definite by always selecting P22 ≫ P11 . Finally, the state feedback operator can be calculated from (23) which will be: (31) K(t) = −R−1 B¯d (t)∗ Φ B(t)∗ Ψ
Fig. 1. CSTR-PFR system and indirectly, and through the CSTR. With the assumptions of negligible diffusion in the PFR, perfect level control in the CSTR, no transportation lags in the connecting lines, constant fluid velocity in the PFR with respect to spatial coordinate, and constant physical properties, and incompressible fluid, the mathematical model of the system will be: Fin in dCA = (CA − CA ) − k1 e−E1 /RT CA − dt Vs 2 k3 e−E3 /RT CA (34) Fin dCB =− CB + k1 e−E1 /RT CA − dt Vs k2 e−E2 /RT CB (35) 1 dT = [k1 e−E1 /RT CA (−∆H1 )+ dt ρcp 2 k2 e−E2 /RT CB (−∆H2 ) + k3 e−E3 /RT CA (−∆H3 )]+ Q Fin (Tin − T ) + (36) Vs ρcp Vs
4. CASE STUDY In this section we consider a CSTR-PFR configuration shown in Fig. 1 as an interacting lumped and hyperbolic distributed parameter system. This reactor configuration is recommended for some types of chemical reactions (see e.g., Fogler (2005) and Schweier and Floudas (1999)) and may be used to carry out Van de Vusse reaction to achieve the maximum conversion to the desired product (Schweier and Floudas (1999)). Here, we assume reactions and kinetics: A −→ B, −r1 = k1 e−E1 /RT CA B −→ C,
−r2 = k2 e−E2 /RT CB −E3 /RT
2 CA
(32)
2A −→ D, −r3 = k3 e where k1 , k2 and k3 are pre-exponential constants; E1 , E2 , and E3 are the activation energy; and R is the universal gas constant. The exothermic reactions take place in both CSTR and PFR and component B is the desired component. The reaction kinetics are considered to be time-varying according to the following exponential decay model (see Lie and Himmelblau (2000)): ki = k0,i + k1,i e−αi t (33) where subscripts i = 1, 2, 3 denote number of reactions and k0,i , k1,i , and αi are the decay model parameters. The objective is to control the concentration of the components and the temperature within both reactors by using inlet flow rate (Fin ) and cooling rate from the CSTR (Q) as manipulated variables. It should be noticed that when liquid level is perfectly controlled in the CSTR, variations of Fin affects the controlled variables in the PFR directly
p ∂CA ∂C p p = −υ A − k1 e−E1 /RT CA − ∂t ∂z p2 (37) k3 e−E3 /RT CA p p ∂CB ∂CB p = −υ + k1 e−E1 /RT CA − ∂t ∂z p k2 e−E2 /RT CB (38) ∂Tp ∂Tp k1 −E1 /RT p e CA (−∆H1 )+ = −υ + ∂t ∂z ρcp k2 −E2 /RT p e CB (−∆H2 )+ (39) ρcp k3 −E3 /RT p 2 e CA (−∆H3 ) ρcp p CA (t, 0) = CA (40) p (41) CB (t, 0) = CB Tp (t, 0) = T (42) where CA and CB are the concentration of the components A and B in the CSTR, respectively; T is the temperature p p in the CSTR; CA and CB are the concentration of the components A and B in the PFR, respectively; Tp is the temperature in the PFR; z ∈ [0, L] is the spatial in coordinate; t ∈ [0, ∞] is the time; Fin , CA , and Tin are the volumetric flow-rate, concentration, and temperature of the feed to the CSTR; Vs and Vp are the volumes of the CSTR and the PFR, respectively; υ is the fluid velocity in the PFR which is given by υ = FVinpL ; ∆H1 , ∆H2 , and ∆H3 are the heat of reaction for reactions 1, 2, and 3, respectively; and ρ and cp are the average fluid density
14449
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
Table 1. Model parameters Concentration (kmol/m3)
Value 150.150 × 106 sec−1 225.225 × 106 sec−1 1.759 × 106 sec−1 100.1 × 106 sec−1 25.025 × 106 sec−1 0 sec−1 1.389 × 10−2 sec−1 174.845 × 10−6 m3 /sec −1.36 kJ/sec 5.1 kmol/m3 403.15 K −4200 kJ/kmol −11000 kJ/kmol −41850 kJ/kmol 9758.3 K 9758.3 K 8560.0 K 0.01 m3 0.005 m3 934.2 kg/m3 3.01 kJ/kgK
3
419
2.5
418 417
Tp
416
1.5
415
1 0.5 0
414 0.05
0.1
0.15 z (m)
Fig. 2. Steady-state profiles in the PFR
0.2
0.25
413 0.3
Temperature (K)
Concentration (kmol/m3)
We use the proposed optimal policy to control the concentration of the components and also the temperature in both CSTR and PFR. Cd (t, z) and Cl (t) are selected to be the identity matrices to yield a state LQR problem. In order to solve the LQR control problem for this system, we follow the procedure discussed in section 3. The set of differential and algebraic equations (28) to (30) are solved numerically in gPROMS. By choosing P11 = diag(80, 80), P22 = diag(350, 350), and R = diag(0.01, 0.0001) the matrix P will be positive semi-definite. The feedback operator can then be obtained using (31).
CpB
A
C (open−loop)
8
B
C (closed−loop) B
6 4 2 100
200 300 Time (sec)
400
500
520
4.1 Controller Design
CpA
A
C (closed−loop)
Fig. 3. Concentration responses in the CSTR
and specific heat. The model parameters used in this case study are given in table 1. In order to find the equilibrium condition for the system, the modelling equations (34) to (42) have R been solved at steady-state in gPROMS (gpr (19972010)). Simulation yields the steady-state values for CA , CB , and T ; 2.89 kmol/m3 , 0.80 kmol/m3 , and 413.10 K, p p respectively. The steady-state profiles for CA , CB and Tp are shown in Fig. 2. The model equations can be linearized around the steady-state condition to describe the system as the general formulation (1) to (6).
2
C (open−loop)
10
0 0
Temperature (K)
Parameter k0,1 k0,2 k0,3 k1,1 k1,2 k1,3 α1 , α2 , α 3 Fin,ss Qss in CA Tin ∆H1 ∆H2 ∆H3 E1 /R E2 /R E3 /R Vs Vp ρ cp
12
500
T (open−loop) T (closed−loop) Tp (open−loop)
480
Tp (closed−loop)
460 440 420 400 0
100
200 300 Time (sec)
400
500
Fig. 4. Temperature responses in the CSTR and at the outlet of the PFR
4.2 Simulation Results In order to assess the performance of the control policy, the designed feedback operator is used with the original non-linear system (34) to (42). The set of coupled nonlinear PDEs and ODEs are solved using orthogonal collocation on finite element method in gPROMS. We use CA (0) = 12 kmol/m3 , CB (0) = 12 kmol/m3 , T (0) = 400 p p K, CA (0, z) = 0.7 kmol/m3 , CB (0, z) = 0.7 kmol/m3 , and Tp (0, z) = 400 K as the initial conditions. To have a measure of how good the designed controller is, the responses of open-loop and closed-loop systems for the concentrations and the temperature in the CSTR and at the outlet of the PFR are compared in Fig. 3 to Fig. 5. As it can be seen, the controlled system is able to reject the effect of the initial condition about 2 times faster (with respect to the residence time) than the open-loop system, and it converges to the desired steady-state condition, accurately. Moreover, the open-loop inverse-responses appearing at the outlet of the PFR have disappeared. The variations of the control inputs are also shown in Fig. 6. The control efforts are not particularly aggressive, and are physically realizable. 5. CONCLUSION In this work LQR control problem for a class of timevarying, composite lumped and distributed parameter system is formulated and solved. In this mixed system, the lumped parameter system interacts with the hyperbolic distributed parameter system through its boundary. The control variables can affect the boundary of the distributed parameter system directly or indirectly, and through the
14450
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
Concentration (kmol/m3)
3 p C A p C A p C B p C B
2.5 2 1.5
(open−loop) (closed−loop) (open−loop) (closed−loop)
1 0.5 0 0
100
200 300 Time (sec)
400
500
Fig. 5. Concentration responses at the outlet of the PFR −4
x 10
−1
4
Fin
3.5
Q
−1.5
−2
3
2.5
Q (kJ/sec)
Fin (m3/sec)
4.5
−2.5
2 1.5 0
100
200 300 Time (sec)
400
−3 500
Fig. 6. Control input profiles lumped parameter system. The LQR control problem is formulated based on an infinite-dimensional state-space representation of the system, which is obtained via a state transformation from the original system using boundary control transformation method. The solution of the LQR control problem is achieved by solving the matrix Riccati equation that results from the ORE of the infinitedimensional state-space representation. The designed optimal control policy was implemented on an interacting CSTR-PFR system. The performance of the controller was assessed by implementing the controller on the original non-linear system and resulted in a high performance closed-loop system. REFERENCES (1997-2010). Process systems enterprise, gPROMS, www.psenterprise.com/gproms. Acquistapace, P. and Terreni, B. (1984). Some existence and regularity results for abstract non-autonomous parabolic equations. Journal of Mathematical Analysis and Applications, 99(1), 9–64. Aksikas, I. and Forbes, J.F. (2010). Linear quadratic regulator of time-varying hyperbolic distributed parameter systems. IMA Journal of Mathematical Control and Information, 27(3), 387–401. Aksikas, I., Fuxman, A., Forbes, J.F., and Winkin, J. (2009). LQ control design of a class of hyperbolic PDE systems: Application to fixed-bed reactor. Automatica, 45, 1542–1548. Aksikas, I., Winkin, J., and Dochain, D. (2007). Optimal LQ-feedback regulation of a nonisothermal plug flow reactor model by spectral factorization. IEEE Transaction on Automatic Control, 52(7), 1197–1193.
Aksikas, I., Winkin, J., and Dochain, D. (2008). Optimal LQ-feedback for a class of first-order hyperbolic distributed parameter systems. ESAIM:COCV. Bensoussan, A., Prato, G.D., and Mitter, K. (2007). Representation and Control of Infinite Dimensional Systems. Systems and Control: Foundations & Applications. Birkh¨ auser Boston. Callier, F. and Winkin, J. (1990). Spectral factorization and LQ-optimal regulation for multivariable distributed systems. International Journal of Control, 52(1), 55–57. Christofides, P.D. (2001). Nonlinear and robust control of PDE systems: Methods and applications to transportreaction processes. Boston:Brikhauser. Curtain, R.F. and Zwart, H.J. (1995). An Introduction to Infinite Dimensional Linear Systems. Springel-Verlag. Fattorini, H.O. (1968). Boundary control systems. SIAM Journal on Control, 6, 349–385. Fogler, H.S. (2005). Elements of Chemical Reaction Engineering. Prentice-Hall. Hiratsuka, S. and Ichikawa, A. (1969). Optimal control of systems with transportation lags. IEEE Transaction on Automatic Control, 14(3), 237. Krstic, M. and smyshlyaev, A. (2008). Boundary control of PDEs: A course on backstepping designs. Society for industrial and Applied Mathematics, Philadelphia, PA. Lie, B. and Himmelblau, M. (2000). Catalyst deactivation: control relevance of model assumptions. Ind. Eng. Chem. Res, 39, 1242–1248. Moghadam, A.A., Aksikas, I., Dubljevic, S., and Forbes, J.F. (2010). LQ control of coupled hyperbolic PDEs and ODEs: Application to a CSTR-PFR system. Proceedings of the 9th International Symposium on Dynamics and Control of Process Systems, DYCOPS 2010, 713–718. Oh, M. (1995). Modelling and simulation of combined lumped and distributed processes. Ph.D. thesis, Imperial College, London. Orlov, Y. and Utkin, V. (1987). Sliding mode control in infinite-dimensional systems. Automatica, 23(753-757). Pazy, A. (1983). Semigroup of linear operators and application to partial differential equations. Springel-Verlag, New York. Schweier, C.A. and Floudas, C.A. (1999). Optimization framework for the synthesis of chemical reactor networks. Ind. Eng. Chem. Res, 38, 744–766. Thowsen, A. and Perkins, W.R. (1975). Sampled-data linear-quadratic regulator for systems with generalized transportation lags. International Journal of Control, 21, 49–65. Thowsen, A. and Perkins, W.R. (1973). Optimal discretetime feedback control of mixed distributed and lumped parameter systems. International Journal of Control, 18(3), 657–665. Tzafestas, S. (1970a). Final-value control of nonlinear composite distributed and lumped parameter systems. Franklin Inst., 290(5). Tzafestas, S. (1970b). Optimal distributed-parameter control using classical variational theory. International Journal of Control, 12(4), 593–608. Wang, P.K.C. (1966). On the stability of equilibrium of a mixed distributed and lumped parameter control system. International Journal of Control, 3(2), 139–147.
14451