Chemical Engineering Science 60 (2005) 969 – 980 www.elsevier.com/locate/ces
Feedback control of hyperbolic distributed parameter systems Huilan Shanga,∗ , J. Fraser Forbesb , Martin Guayc a School of Engineering, Laurentian University, Sudbury, Ont., Canada P3E 2C6 b Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Al., Canada T6G 2G6 c Department of Chemical Engineering, Queen’s University, Kingston, Ont., Canada K7L 3N6
Received 1 December 2003; received in revised form 26 August 2004; accepted 7 September 2004 Available online 23 November 2004
Abstract Hyperbolic distributed parameter systems (DPS) represent a large number of industrial processes with spatially nonuniform operating variable profiles. Research has been conducted to develop high-performance control strategies for these systems by exploiting their highfidelity models. In this paper, a feedback control method that yields improved performance is proposed for DPS modelled by first-order hyperbolic partial differential equations (PDEs) using the method of characteristics. Simulation results show that this method can provide effective control for the systems modelled by a scalar PDE as well as a system of PDEs. Further, it can efficiently compensate the effect of model-plant mismatch and effectively reject the disturbances. 䉷 2004 Published by Elsevier Ltd. Keywords: Feedback control; Distributed parameter systems; Method of characteristics; Hyperbolic PDE; Process control
1. Introduction Distributed parameter systems (DPS) are a class of important processes in which process variables vary in space as well as in time. It was claimed that almost all natural and industrial processes were distributed in nature (Gay and Ray, 1995). Herein DPS refer to those systems whose obviously nonuniform spatial profiles are of interest in process and control engineering (e.g., steel-making, fluid heat exchangers, some chemical reactors, polymer processing operations, etc.) (Butkowskii, 1969; Ray, 1978). The description of DPS, obtained using conservation laws, often takes the form of partial differential equations (PDEs). The PDEs modelling DPS can be hyperbolic, parabolic or elliptic. The hyperbolic PDEs represent the dynamics of industrial processes involved in convection with negligible diffusion
∗ Corresponding author. Tel.: +1 705 675 1151X2176; fax: +1 705 675 4862. E-mail address:
[email protected],
[email protected] (H. Shang).
0009-2509/$ - see front matter 䉷 2004 Published by Elsevier Ltd. doi:10.1016/j.ces.2004.09.067
effect (e.g., fluid heat exchanger, plug-flow reactors (PFRs) and fibre spinlines). The control of DPS can be performed by approximating DPS into lumped parameter systems (LPS), in which only the variation of process variables in time is considered, then applying the abundance of the existing control theory and techniques for LPS (Shirvani et al., 1995). The control performance using such approximation techniques can be limited due to inaccuracy of the models. Since the control development based on high-fidelity models has the potential to improve control performance, research has become increasingly active on developing efficient control using accurate PDE models (Christofides, 2000; Godasi et al., 2002; Neittaanmaki and Tiba, 1994; Hernández and Arkun, 1992; Laroche et al., 1998; Alvarez-Ramirez, 2001). In some of these approaches, the infinite-dimensional DPS were controlled by representing DPS with their dominant finite-dimensional modes such as eigenfunctions or singular functions (Ray, 1981; Park and Cho, 1996; Chakravarti and Ray, 1999; Sadek and Bokhari, 1998). The recent control development using finite-dimensional approximation includes finite-dimensional nonlinear control based
970
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
on ordinary differential equations (ODEs) obtained from a combination of Galerkin’s method and approximation inertial manifold, and from Galerkin’s method with empirical eigenfunction basis functions (Christofides and Daoutidis, 1997; Baker and Christofides, 2000). These control methods were shown to be able to yield adequate performance mainly for DPS modelled by parabolic PDEs. Hyperbolic systems, however, are considered to possess infinite-dimensional modes of almost equal energy and cannot be accurately represented by a finite number of modes. Having recognized that the infinite nature of these processes has to be included in the control development, some control approaches have been proposed for hyperbolic DPS. A control approach based on a combination of the method of characteristics and sliding mode techniques was proposed for processes modelled by a firstorder quasi-linear hyperbolic PDE (Sira-Ramirez, 1989). This method was further developed to synthesize a state feedback controller for a nonlinear hyperbolic PDE model with continuous control action (Hanczyc and Palazoglu, 1995a,b). Nonlinear distributed output feedback controller was proposed for a quasi-linear hyperbolic PDE system using geometric control methods (Christofides and Daoutidis, 1996, 1998). The ideas from geometric control theory were also extended to develop velocity control of hyperbolic PDE systems (Gundepudi and Friedly, 1998). The control of the hot spot temperature in a PFR, a specific hyperbolic DPS, was addressed using a nonlinear controller (Karafyllis and Daoutidis, 2002). These existing control methods can usually produce satisfactory performance; however, the control of DPS is a broad area and further research is required to develop control strategies for more general processes which can be easily implemented in industrial practice. In this paper, a novel feedback control method is proposed for systems modelled by first-order hyperbolic PDEs. Using the method of characteristics, the control law is developed such that the closed-loop response of the process output moves towards the desired setpoint. After the work by Sira-Ramirez (1989) and Hanczyc and Palazoglu (1995b), the method of characteristics is recalled in feedback control development. Compared to these existing control methods, the proposed approach does not require a continuous spatial profile for the setpoints, which can be hard to compute in complex systems. It makes full use of the measurement information without demanding infinite-dimensional state estimation. The proposed controller can provide effective control with high performance for processes having perfect models. The control algorithm is further developed to be able to tolerate model-plant mismatch and numerical computation error. The control method in this paper can deal with more complex systems than the existing feedback control methods. Simulation results show that the proposed control method is able to deliver control performance to a higher level.
2. Method of characteristics The method of characteristics is an essential mathematical technique employed to develop the proposed control approach. It provides a geometric solution method to compute integral surfaces of a certain class of DPS by considering integral curves (or characteristics) of a specific vector field (Duchateau and Zachmann, 1989). In this section, the method of characteristics is introduced for quasi-linear and nonlinear scalar first-order PDEs. 2.1. Quasi-linear PDEs Consider the quasi-linear equation for a function v(t, x1 , . . . , xn ) on the manifold M × R of R n+1 n
*v *v + ai (x, v, u) = f (x, v, u), *t *xi i=1
(1)
where t is the time, x = [x1 , . . . , xn ] is a point in the manifold M ∈ R n , v is the distributed variable in the infinite-dimensional Hilbert space H[(0, L), R]. a1 (x, v, u), . . . , an (x, v, u) and f (x, v, u) are continuous (C 0 ) in x and v, and u ∈ R m are the parameters or the process inputs. Given v(t, x) as the solution of Eq. (1), consider the graph z = v(t, x). This graph has a normal vector *v *v *v N0 = − , − ,...,− ,1 (2) *t *x1 *x n at the point (t0 , x0 , v(t0 , x0 )). Let z0 = v(t0 , x0 ). Eq. (1) implies that the vector 0 =[1, a1 (x0 , z0 , u), . . . , an (x0 , z0 , u), f (x0 , z0 , u)] is perpendicular to the normal vector N0 , and hence, must lie in the tangent plane to the graph of Eq. (1) at the point (t0 , x0 , z0 ). In other words,
(t, x, v, u) = [1, a1 (x, v, u), . . . , an (x, v, u), f (x, v, u)] (3) defines a vector field in R n+1 , tangent to the solution graph at each point. This vector field is called the characteristic vector field of the quasi-linear equation (1). The integral curves of the characteristic vector field are called the characteristics of the quasi-linear equation. The ODEs defined by the vector field (t, x, v, u) are called the characteristic equations and the characteristic equations of the quasi-linear equation (1) have the form: t˙ = 1, x˙ = a(x, v, u), v˙ = f (x, v, u). (4) If the graph v = v(t, x, u) is a smooth surface S, which is a union of such characteristic curves, then at each point (t, x, v), the tangent plane contains the vector (x, v, u); hence, S must be an integral surface. In other words, a smooth union of characteristic curves is an integral surface of the characteristic vector field. If the given initial condition is noncharacteristic (i.e., the curve from the initial condition
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
is nowhere tangent to the vector field), a simple procedure for solving the first-order PDE problem is to flow out from each point of the initial curve along the characteristic curve through that point, thereby sweeping out an integral surface. This is the method of characteristics (McOwen, 1996). It is trivial to show that the function so constructed satisfies the original differential equation. By differentiating along the characteristics, it can be obtained: n
n
*v *v *v *v t˙ + + v˙ = f = x˙i = ai (x, v, u) . *t *x i *t *x i i=1 i=1 (5) This result shows that along the characteristics, the partial derivative terms in Eq. (1) reduce to a directional derivative of v in that direction. 2.2. Nonlinear PDEs The characteristics are more complicated for a first-order nonlinear PDE system represented by: *v *v + v, x, , u = 0. (6) *t *x Eq. (6) can be interpreted as a hypersurface in the manifold M 2n+3 = J 1 (V n+1 , R) of 1-jets equipped with the standard contact structure. Let (x, t) be the local coordinates on V n+1 and v be the coordinate in R, p = *v/*x and q = *v/*t. Denoting the corresponding local coordinates by (x, t, v, p, q) in the space of 1-jets, the differential equation (6) can be written in the form: q + (v, x, p, u) = 0.
(7)
By the method of characteristics, the characteristic equations of PDE (6) in the manifold M 2n+3 can be written as (Arnold, 1988): t˙ = 1, q˙ = −q v , x˙ = p , p˙ = −x − p v , v˙ = p p + q,
(8)
971
the characteristic vector field. Exploiting this technique in the control development has the potential to produce control algorithms with simple computation and desirable control performance. 3. Feedback control development By representing the hyperbolic PDE models using their characteristics, a control algorithm for hyperbolic DPS can be developed to drive the closed-loop output response towards its setpoint along the characteristic curves. In this paper, we confine ourselves to processes with one spatial dimension. The output to be controlled is assumed to be the value of infinite-dimensional variables at outlet locations. Without involving much complexity, the proposed control method is extendable to a wide variety of cases (e.g., the spatial average of infinite-dimensional variables) and multiinput multi-output processes. 3.1. Scalar quasi-linear systems Consider a dynamic system described by the scalar firstorder hyperbolic PDE:
*v *v + a(v, x, u) = f (v) + g(v)u, *t *x v|xb = vb , y = h(v)|xo ,
(10)
where t is the time, x is the spatial coordinate, v ∈ H[(0, L); R] is the distributed state variable, u ∈ R is the manipulated variable, a(v, x), f (v) and g(v) are continuous functions, vb is the boundary value of the state variable v at the boundary location xb , the output y to be controlled is the value of an infinite-dimensional function h(v) at outlet location xo . Since it is often possible to develop the coordinate system to have a(v, x, u) positive, we assume that a(v, x, u) > 0 for every admissible input throughout the paper. By the method of characteristics, the quasi-linear PDE system in Eq. (10) can be transformed into a nonlinear ODE system in the space of (t, x, v):
where the subscript denotes a partial derivative (e.g., = */*). Substitution of Eq. (7) can reduce the order of Eq. (8) by 1:
t˙ = 1, x˙ = a(v, x, u), v˙ = f (v) + g(v)u.
t˙ = 1, x˙ = p , p˙ = −x − p v , v˙ = p p − .
For the process to be controllable, it requires that g(v) and *a(v, x, u)/*u cannot both be zero. The variation of (t, x, v) along the characteristic curves can be described by the vector field = [1, a(v, x, u), f (v) + g(v)u]. The Lie derivative of the function h(v) along the vector field is:
(9)
Notice that, in contrast to a quasi-linear hyperbolic PDE, the characteristic equations of the nonlinear PDE contains the state variable as well as it spatial derivatives. As seen above, the method of characteristics allows one to transform hyperbolic PDE models to ODE forms along
(11)
*h *h f (v) + g(v)u. (12) *v *v Assuming that the output (value of h(v) at xo ) has the setpoint r, the control action is to be formulated in such a way L h(v) =
972
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
that h(v) reaches the setpoint r while x reaches xo along the characteristic curves. The controller is obtained by satisfying: xo L h(v) dx, (13) r − h(vb ) = xb a(v, x, u)
this subsection addresses the feedback control design for vector quasi-linear PDE systems under the condition that u does not appear in a(v, x, u) (i.e., a(v, x, u) = a(v, x)). The variation of (t, x, v) can be described by the vector field = [1, a(v, x), f(v) + g(v)u] and the corresponding characteristic ODEs are:
where h(vb ) is calculated using the given boundary condition. In Eq. (13), the control objective is achieved by adjusting L h(v) and/or the residence time of the characteristic curves through the manipulated variable u. The expression for u can be derived analytically if u does not explicitly exist in a(v, x, u), i.e., a(v, x, u) = a(v, x), or if a(v, x, u) can be written as a(v, x, u) = a (v, x)u. If a(v, x, u) = a(v, x), substituting Eq. (12) into Eq. (13) yields: xo *h/*v f (v) + *h/*v g(v)u r − h(vb ) = dx. (14) a(v, x) xb
x˙ = a(v, x), v˙ = f(v) + g(v)u.
Then the manipulated variable can be obtained as: x r − h(vb ) − xbo *h/*v f (v)/a(v, x) dx xo u= . xb *h/*v g(v)/a(v, x) dx
(15)
If a(v, x, u) = a (v, x)u, Eq. (13) can be solved and u takes the form of: xo xb *h/*v f (v)/a (v, x) dx xo u= . (16) r − h(vb ) − xb *h/*v g(v)/a (v, x) dx If the function a(v, x, u) does not take the above two forms, no explicit formula for u can be derived and numerical techniques, e.g., Newton’s method, have to be used to solve the nonlinear equation (13) for u. In the controller derived from Eq. (13), integration of the infinite-dimensional functions in spatial direction can be estimated using numerical integration schemes based on current measurement or estimates of the state variables at a finite number of spatial points. If the process model is perfect and the numerical integration error is negligible, the obtained controller is able to quickly drive the closed-loop response of the output to the setpoint, which will be illustrated in the simulation results. 3.2. Vector quasi-linear systems
*v *v + a(v, x, u) = f(v) + g(v)u, *t *x v|xb = vb , y = h(v)|xo ,
For control development, the following assumption is needed. Assumption 1. There exists an integer such that the Lie derivatives of the output function satisfy −2
Lg h(v) = Lg L1f h(v) = · · · = Lg Lf
−1 Lg Lf h(v)
where v = [v1 , v2 , . . . , vn ] is the vector of the distributed parameter variables with v1,2,...,n ∈ H[(0, L); R]. If *a(v, x, u)/*u = 0, the feedback control can be developed in the same way as in the last subsection. Therefore,
= 0,
(19)
When Assumption 1 is satisfied, the Lie derivatives of the output function h(v) along the characteristic vector field can be written as: Li h(v) = Lif h(v), i = 1, . . . , − 1,
−1
L h(v) = Lf h(v) + Lg Lf
h(v)u.
(20)
Then the process is said to have the relative degree of , using the idea from nonlinear ODE systems (Isidori, 1995). Note that the similar concept of characteristic index was introduced for PDE systems in Christofides and Daoutidis (1996). In this paper, the term relative degree is preferred over characteristic index because, along the characteristic vector field, the directional derivative of the output function h(v) is the same as in an ODE system. When the DPS modelled by systems of PDEs have the relative degree of , the controller can be developed for the system in Eq. (17) such that r − h(vb )
= Lf h(v)|x=xb −1
(17)
h(v) = 0,
∀ v, where Lf and Lg indicate Lie derivative or directional derivative.
+ Lf
Consider a dynamic system described by the system of first-order hyperbolic PDEs:
(18)
xo xb
h(v)|x=xb
1 dx1 + · · · a(v, x) xo x−2 1 ··· xb a(v, x) xb
1 dx−1 · · · dx1 a(v, x) xo x−1 Lf h(v) 1 + ··· dx · · · dx1 a(v, x) xb a(v, x) xb x−1 xo −1 Lg Lf h(v)u 1 ··· + a(v, x) xb a(v, x) xb × dx · · · dx1 . ×
(21)
Eq. (21) is to adjust the directional derivatives of the output function h(v) and drive it towards the setpoint r at xo .
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
The control action derived from the above equation has the following form:
973
For the scalar quasi-linear PDE system (10), a robust version of the proposed control approach is developed under the following assumption.
x r − h(vb ) − Lf h(v)|x=xb xbo 1/a(v, x) dx1 − u= x − ··· x−1 −1 o Lg Lf h(v)/a(v, x) dx · · · dx1 xb 1/a(v, x) · · · xb −
xo
x x−1 xo 1/a(v, x)· · · xb−2 1/a(v,x) dx−1 · · · dx1 xb 1/a(v, x)· · · xb Lf h(v)/a(v, x) dx · · · dx1 − x . xo x−1 x−1 −1 −1 o xb 1/a(v, x)· · · xb Lg Lf h(v)/a(v, x) dx · · · dx1 xb 1/a(v, x)· · · xb Lg Lf h(v)/a(v, x) dx · · · dx1
−1
Lf
h(v)|x=xb
xb
(22) This equation gives a general formula of the feedback controller for processes modelled by a system of hyperbolic quasi-linear PDEs. The controller can enforce the output towards the setpoint in the time scale of the residence time. The dynamics of the closed-loop response under the controller (22) is strongly affected by the relative degree . When = 1, the process is expected to have relatively simple closed-loop response and the output reaches the setpoint in approximately the residence time or less. As increases, the complexity of the closed-loop response increases and it often requires more time than the residence time for the output to settle on the setpoint. In Eq. (22), control calculation requires multiple integration of the functions based on state measurement or estimation, inevitably leading to numerical errors. Therefore, it is necessary to develop a robust technique to correct the effect of numerical errors. 3.3. Robust control development For complex or highly nonlinear processes, integration using state measurement or estimates at a limited number of spatial locations may lead to numerical errors and cause offset in output response. Moreover, if the models are not perfect, model-plant mismatch can also result in offset using either controller (13), (15), (16) or (22). Therefore, it is of interest to develop a modified version of the proposed controller so that it can mitigate the effects of model-plant mismatch and numerical integration errors. Research has been conducted on robust control of PDE systems. In Christofides and Daoutidis (1998), the robust technique for hyperbolic PDEs was developed to address uncertain variables or unmodelled dynamics. In this paper, a robust technique is required to simultaneously attenuate the effect of the numerical integration error, time-varying parameter uncertainty, unmodelled dynamics and disturbances of unknown sources. Due to the multiplicity of the uncertain sources, the robust problem discussed here cannot be solved using disturbance decoupling or singular perturbation framework. In this section, the control methods discussed in Sections 3.1 and 3.2 are further developed to address the effect of model-plant mismatch and numerical integration error using the modification of conventional integration technique.
Assumption 2. Assume that *h/*v g(v) is always positive (or negative) if it is not zero along its spatial profile in model (10). Assumption 2 can often be satisfied in real processes. If Assumption 2 is satisfied, from Eq. (12), h(v) at all spatial locations tends to vary in the same direction if the control action varies. For example, under the constant boundary condition, if *h/*v g(v) is positive, h(v) at all points along the characteristic curves increases when u increases. Therefore the control law in Eq. (13) can be modified such that r − h(vb ) =
xo xb
L h(v) dx + a(v, x, u)
t
(y − r) d,
(23)
0
where is a positive tuning parameter. is taken to be zero when there is no model-plant mismatch or numerical integration error. Otherwise, can be tuned to get the desired output response. When a(v, x, u) = a(v, x), from Eq. (12) and (23), the controller u is obtained:
u=
r−h(vb )−
t
xo 0 (y−r) d− xb *h/*v f (v)/a(v, x) dx xo . xb *h/*v g(v)/a(v, x) dx (24)
If a(v, x, u) = a (v, x)u, from Eq. (23), the expression of u takes the form of: xo
u=
x *h/*v f (v)/a (v, x) dx . tb xo r−h(vb )− 0 (y−r) d− xb *h/*v g(v)/a (v, x) dx (25)
If a(v, x, u) is not in these two forms, the value of u has to be solved using numerical techniques. When the manipulated variable u does not appear in the characteristic of the spatial coordinate x, i.e., a(v, x, u) = a(v, x), the effect of the output error integral in controller (24) can be investigated by examining the variation of h(v)
974
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
along the characteristic vector field : L h(v) =
*h f (v) *v
of current state variables and output variables. It does not contain derivatives of the state variables, as do some of the existing feedback control methods for quasi-linear PDE
t *h/*v g(v)(r−h(vb ))−*h/*v g(v) xxbo *h/*v f (v)/a(v, x) dx *h/*v g(v) xo + (y−r) d. − xo xb *h/*v g(v)/a(v, x) dx 0 xb *h/*v g(v)/a(v, x) dx (26)
x Note that *h/*vg(v)/ xbo *h/*vg(v)/a(v, x) dx is positive in Equation (26). If y > r, the effect of the time integral term is to let L h(v) decrease. That would gradually result in the decrease of h(v) at all locations except at xb (see Fig. 1). Therefore, the output error integral tends to tune the variation of h(v) along the characteristic vector field such that h(v) reaches setpoint at xo . −1 Similarly, under the assumption that Lg Lf h(v) is always positive (or negative) along its spatial profile, the robust control for systems of PDEs (17) can be obtained from the modification of Eq. (21): r − h(vb )
= Lf h(v)|x=xb
xo xb
−1
1 dx1 + · · · a(v, x) xo 1 ··· xb a(v, x)
+ Lf h(v)|x=xb x−2 1 dx−1 · · · dx1 a(v, x) xb xo x−1 Lf h(v) 1 + ··· dx · · · dx1 a(v, x) xb a(v, x) xb xo x−1 −1 Lg Lf h(v)u 1 + ··· dx · · · dx1 a(v, x) xb a(v, x) xb t + (y − r) d (27) 0
and it takes the form of:
systems. This makes the proposed feedback control easy to implement. 3.4. Nonlinear systems Consider the first-order nonlinear PDE control system: * * + , x, , u = 0, *t *x v|xb = vb , y = h(v)|xo , (29) where is a smooth function. In contrast to a quasi-linear system, the nonlinear PDE system requires its state variable and its first-order derivative to describe its characteristics, which complicates the control design problems. As described in Section 2, the characteristics of the nonlinear PDE can be described by the variables (t, x, v, p). In the local coordinate of (t, x, v, p), the PDE in Eq. (29) has a characteristic vector field =[1, p , pp − , −x −p v ]. The Lie derivative of the function h(v) along the vector field is L h(v) =
*h [p p (v, x, p, u) − (v, x, p, u)]. *v
(30)
Then the control law can be designed for the system in Eq. (29) to satisfy: r − h(vb ) xo *h/*v[pp (v, x, p, u) − (v, x, p, u)] dx (31) = p (v, x, p, u) xb
t x r − h(vb ) − 0 (y − r) d − Lf h(v)|x=xb xbo 1/a(v, x) dx1 − u= − ··· xo x−1 −1 Lg Lf h(v)/a(v, x) dx · · · dx1 xb 1/a(v, x) · · · xb −
xo
x−1 x xo 1/a(v, x)· · · xb−2 1/a(v, x) dx−1 · · · dx1 xb 1/a(v, x)· · · xb Lf h(v)/a(v, x) dx · · · dx1 − x . xo x−1 x−1 −1 −1 o xb 1/a(v, x)· · · xb Lg Lf h(v)/a(v, x) dx · · · dx1 xb 1/a(v, x)· · · xb Lg Lf h(v)/a(v, x) dx · · · dx1
−1
Lf
h(v)|x=xb
xb
(28) Using the output error integral in Eqs. (23), (24), (27) and (28), the controller is able to provide effective feedback control for processes with model-plant mismatch and numerical errors. This integral term can also relax accuracy requirements for state variable measurement or estimation. The controller in Eqs. (23), (24), (27) and (28) is a function
and the control action can be calculated from the above equation using a numerical technique such as Newton’s method. The control development based on the above equation is to adjust L h(v) and/or residence time along the characteristic curves so that h(v) reaches the setpoint at xo . The obtained control law expects to provide effective output tracking
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
975
h(v)
the feedback control proposed by Hanczyc and Palazoglu (1995b). In this section, the performance of the proposed control algorithm is evaluated using this example and compared with that by Hanczyc and Palazoglu. The dynamic model of the process in deviation form can be expressed as:
*T *T +u + H (T − Tj ) = 0, *t *x T |x=0 = Tb ,
r
xb
xo x
Fig. 1. Evolution of the output function along characteristics.
control for nonlinear processes with no model-plant mismatch and numerical integration error. Similar to quasi-linear systems, controller (31) for the nonlinear system can be modified to mitigate the effect of model-plant mismatch, disturbance and numerical integration errors. To achieve this, the following assumption is needed. Assumption 3. *L h(v)/*u is always positive (or negative) if it is not zero in the range of possible v of process (29). When Assumption 3 is satisfied, a controller robust to model-plant mismatch can be designed such that r − h(vb ) xo *h/*v[pp (v, x, p, u) − (v, x, p, u)] dx = p (v, x, p, u) xb t + (y − r) d, (32) 0
(33)
where T (x, t) denotes the temperature of the heat exchanger within the tube, x ∈ [0, 1], u, the manipulated variable, is the fluid velocity in the exchanger, H is a positive constant, and Tj denotes the jacket temperature. For simplicity, the steam temperature is considered to be spatially uniform. Temperature at the exit, T (x = 1) is the variable to be controlled. The characteristic vector field for Eq. (33) is = [1, u, −H (T − Tj )]. In other words, the variation of the state variable T along can be written as: L T = −H (T − Tj ).
(34)
In this example, the performance of the proposed controller for processes with perfect models is investigated. The controller is developed such that Tsp − Tb =
1 0
(L T )
1 dx, u
(35)
where Tsp is the setpoint of T (x = 1). The manipulated variable can be obtained from Eqs. (34) and (35): u=
H (Tj − T ) , Tsp − Tb
(36)
where is a positive tuning parameter. The control law from Eqs. (31) and (32) for the nonlinear PDE systems is a function of the state variable v as well as its first-order spatial derivative p. Based on the state variable measurement, the value of p can be approximated by the finite difference of v. The possible error due to the finite difference approximation can also be corrected by the output error integral term.
where T indicates the average of the current temperature profile. Note that the control law in Eq. (36) only requires simple algebraic computation. This same process was investigated by Hanczyc and Palazoglu using the combination of sliding mode control and the method of characteristics, resulting in the following control algorithm:
4. Simulation
u=
In this section, the performance of the proposed feedback controller is investigated through two examples: a heat exchanger modelled by a single PDE and a PFR by a system of first-order PDEs. 4.1. Heat exchanger A forced-flow, steam-jacketed, tubular heat exchanger with uniform jacket temperature was used to illustrate
x − H (T − T ) T − Tsp j x |x=1 Tsp x=0
,
(37)
x represents the setwhere is a tuning parameter and Tsp point spatial profile of the temperature. Note that the controlled output in Hanczyc and Palazoglu (1995a,b) was the temperature spatial profile, while the outlet temperature is controlled in the proposed approach. Recognizing that this control law may produce process response with a steady state offset, Hanczyc and Palazoglu also proposed a control law with integral action. The control law with integral action
976
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980 3.2
3
Tsp
T (x=1)
T (x=1)
3
2.4
Tsp
2.8
proposed control control I (Hanczyc et al.)
proposed control
control II (Hanczyc et al.)
control II (Hanczyc et al)
2.6
1.8 0
1
2
0
3
1
2
3
4
5
t
t Fig. 2. Output response to a setpoint change in the heat exchanger.
Fig. 3. Output response to a disturbance in Tb from 0 to 0.5 at t = 1 in the heat exchanger.
takes the form of: u=
x + (+1/i )T −Tsp i
x dt−H (T −T ) T −Tsp j
x |x=1 Tsp x=0
,
proposed control
3.1
(38)
T (x=1)
where i is the integration tuning parameter. It is noted that the controllers in Eqs. (37) and (38) by Hanczyc and Palazoglu have more involved formulation than the proposed controller in Eq. (36). They also require estimation of the setpoint spatial profile which, in many other applications, may require solving complex PDEs or ODEs. A simulation was performed to evaluate the performance of the feedback controller in Eq. (36). The parameters used for simulation were: H = 2, Tj = 5. The initial temperature profile was T (x, 0) = 2x and the boundary condition were Tb =0. For simulation purposes, the finite difference method was used to derive a finite-dimensional approximation of the original PDE equation, with a choice of 100 discretization points. To highlight the performance, a large setpoint change was applied. With the specification of the exit temperature as Tsp = 3, the output response of the process using the controller in Eq. (36) was investigated and compared with that x by Hanczyc and Palazoglu. The setpoint spatial profile Tsp required in the controllers by Hanczyc and Palazoglu can be calculated analytically from the steady state modelling ODE (i.e., the steady state of PDE (33)) with two constraint conx| x 1 ditions: Tsp x=0 =Tb =0 and Tsp |x=1 =Tsp =3. From Fig. 2, it is observed that both the proposed control and controller (38) by Hanczyc and Palazoglu yield smooth and quick convergence of the process output to the setpoint without offset, while controller (37) by Hanczyc and Palazoglu produces output response with offset. The performance of the proposed control was also investigated for the case of disturbances in boundary condition and operating variables. Fig. 3 shows the output response of the proposed control and that of the control with inte-
control II (Hanczyc et al.)
3
0
1
2
3
4
5
t Fig. 4. Output response to a disturbance in Tj from 5 to 6 at t = 1 in the heat exchanger.
gral action by Hanczyc and Palazoglu. The result indicates that the proposed control rejects the disturbance in boundary condition quickly and maintains the output at setpoint. Similar results were also obtained for disturbances in operating variable Tj , as shown in Fig. 4. This process has a relatively simple PDE model and the average temperature of current temperature profile can be estimated accurately using a finite set of temperature measurement. Therefore, controller (36) may produce accurate control without requiring the output error integral term. If the process has model-plant mismatch, the output error integral term is needed to eliminate the output offset. The performance of the controller with the output error integral term is investigated in the next example.
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
977
4.2. PFR
The system in Eq. (39) can be described by ODEs along the vector field:
In a PFR, reactants pass through the vessel with little mixing (Aris, 1969). In the ideal case, each element of the reaction mixture would have a reaction time precisely equal to the residence time of the reactor. In this section, the proposed feedback control will be applied on a PFR with uniform heating in the jacket. Consider the nonisothermal PFR where two first-order reactions in series take place:
t˙ = 1, x˙ = vl , C˙ A = −k10 e−E1 /RT r CA ,
A −→ B −→ C, where A is the reactant species, B is the desired product and C is an undesired product. The reaction rates follow an Arrhenius expression: r1 = −k10 e
−E1 /RT r
CA ,
r2 = −k20 e
−E2 /RT r
CB ,
where k10 , k20 , E1 , E2 are Arrhenius constants and the activation energies of the reactions. The reactions are assumed to be endothermic. The reactor is heated with a jacket. Under the assumptions of no radial concentration gradients in the reactor, constant volume of the liquid in the reactor, constant density and heat capacity of the reactants, and negligible diffusion and dispersion, the following model can be obtained using material and energy balances (Christofides and Daoutidis, 1996):
*C A *CA = −vl − k10 e−E1 /RT r CA , *t *x *C B *C B = −vl + k10 e−E1 /RT r CA − k20 e−E2 /RT r CB , *t *x (−Hr1 ) *Tr *T r = − vl + k10 e−E1 /RT r CA m cpm *t *x (−Hr2 ) + k20 e−E2 /RT r CB m cpm Uw + (Tj − Tr ), m cpm Vr
(39)
subject to the boundary conditions: CA (0, t) = CA0 ,
CB (0, t) = 0,
Tr (0, t) = Tr0 ,
(40)
where CA and CB are the concentrations of the species A and B in the reactor, Tr is the temperature of the reactor, Hr1 and Hr2 are enthalpies of the two reactions, m and cpm is the density and heat capacity of the fluid in the reactor, Vr is the volume of the reactor, Uw is the heat-transfer coefficient, CA0 and TA0 are the concentration and temperature of the inlet stream in the reactor, Tj is the spatially uniform temperature in the jacket and manipulated to control the concentrations, x is the spatial coordinate along the reactor and x ∈ [0, 1].
C˙ B = k10 e−E1 /RT r CA − k20 e−E2 /RT r CB , (−Hr1 ) k10 e−E1 /RT r CA T˙r = m cpm (−Hr2 ) + k20 e−E2 /RT r CB m cpm Uw + (Tj − Tr ) m cpm Vr
(41)
and the corresponding vector field in the space of (t, x, CA , CB , Tr ) is denoted as . The output function is the value of CB at x = 1: y = CB (1, t).
(42)
Then the first- and second-order Lie derivatives of CB along the vector field are: L CB = k10 e−E1 /RT r CA − k20 e−E2 /RT r CB , 2 −2E1 /RT r L L CB = − k10 e CA
2 −2E2 /RT r + k20 e CB − k10 k20 e−(E1 +E2 )/RT r CA k10 E1 CA −E1 /RT r k20 E2 CB −E2 /RT r + e − e RT 2r RT 2r (−Hr1 ) × k10 e−E1 /RT r CA m cpm (−Hr2 ) + k20 e−E2 /RT r CB m cpm Uw + (Tj − Tr ) . (43) m cpm Vr
The input appears in the second-order Lie derivative of the output function, and thus, the system can be said to have relative degree 2. Using the proposed control method in Eq. (27), the controller for this system can be developed such that: 1 x 1 1 r − CB0 = L CB + L L CB 2 dx dx vl x=0 vl 0 0 t + (y − r) d. (44) 0
Note that this equation contains double integral of the current state variable profiles. In addition, the process is highly nonlinear and the state variable profiles might be relatively complex. Therefore, the numerical integration of current state variable profiles could produce numerical error and the outt put error integral term 0 (y − r) d is necessary to eliminate any offset in the output. Substituting the expression of L CB and L L CB in Eq. (43) to the above equation,
978
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
Table 1 Model parameters for a PFR vl = 1 m/min L = 1.0 m Vr = 10.0 l E1 = 20000.0 kcal/kmol E2 = 50000.0 kcal/kmol k10 = 5.0 × 1012 min−1 k20 = 5.0 × 102 min−1 Hr1 = 0.5480 kcal/kmol
R = 1.987 kcal/(min K) m = 0.09 kg/l Uw = 0.20 kcal/(min K) Cpm = 0.231 kcal/(kg K) CA0 = 4 mol/l CB0 = 0 mol/l Tr0 = 320 K Hr2 = 0.9860 kcal/kmol
can be seen that the controller with = 0.6 produces smooth output convergence to the setpoint and eliminates offset due to model-plant mismatch. This result indicates that the proposed control approach can provide robust control for PDE systems with model-plant mismatch. The performance of the proposed control method for disturbances in operating variables was investigated. Fig. 7 illustrates the output response of the reactor when there is disturbances in initial concentration of CA . Simulation results show that the developed control method can reject the effect of the disturbance and regulate the process output at setpoint.
the formula for the manipulated variable v can be obtained: t 1x (r − CB0 − 0 (y − r) d)vl2 − vl L CB |x=0 − 0 0 1 dx dx , Tj = 1x 0 0 2 dx dx
(45)
where 2 −2E1 /RT r 2 −2E2 /RT r
1 = − k10 e CA + k20 e CB
1.4
Simulation was run to evaluate the performance of the controller in Eq. (45). The model parameters used are listed in Table 1. The initial state variable profiles used are CA (x, 0)=4 mol/l, CB (x, 0)=0 mol/l, Tr (x, 0)=320 K and the setpoint of CB at x = 1 is specified as r = 1. To investigate the performance of the controller in the presence of numerical integration error, the double integral term in Eq. (45) was evaluated using the value of distributed variables at 10 spatially discrete points. Fig. 5 compares the output setpoint tracking response using the proposed control method without the output error integral term ( = 0) and that with = 0.25. It is noted that the controller without the output error integral term ( = 0) exhibits offset due to numerical error. The controller with the output error integral term ( = 0.25) can effectively compensate for the numerical error, eliminate the offset and provide satisfactory setpoint tracking response. The performance of the controller in the presence of model-plant mismatch was studied for the process with a parameter error in heat transfer coefficient Uw . Assuming Uw = 0.2 has an error of 0.02 (the true value of Uw is assumed to be 0.18), the output response of the process using controller (45) was investigated, as illustrated in Fig. 6. It
1
0.6 λ = 0.25 λ=0
0.2
0
2
4
6
t, min Fig. 5. Setpoint tracking performance in PFR in the presence of numerical error.
1.4
1 CB at outlet, mol/l
− k10 k20 e CA k10 E1 CA −E1 /RT r k20 E2 CB −E2 /RT r + e − e RT 2r RT 2r (−Hr1 ) × k10 e−E1 /RT r CA m cpm (−Hr2 ) U w Tr −E2 /RT r , + k20 e CB − m cpm m cpm Vr k10 E1 CA −E1 /RT r k20 E2 CB −E2 /RT r
2 = e − e RT 2r RT 2r Uw . (46) × m cpm Vr
CB at outlet, mol/l
−(E1 +E2 )/RT r
0.6 λ = 0.6 λ=0
0.2
0
2
4
6
t, min Fig. 6. Setpoint tracking performance in PFR in the presence of model-plant mismatch.
CB at outlet, mol/l
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
979
as well as systems with uncertainty. To develop infinitedimensional closed-loop stability for all these cases is very important but challenging, which will be addressed in the future work. The proposed method applies to hyperbolic PDE systems with single characteristics (only one characteristic curve passes each point of the solution surface). Further research is required to extend the control method to hyperbolic PDE systems with multiple characteristics and some parabolic systems.
1.2
1
0.8 0
(a)
4
8
12
16
20
5
The authors gratefully acknowledge financial support provided by the Mechanical Wood Pulps Network Centre of Excellence and Natural Sciences and Engineering Research Council of Canada (NSERC).
4
References
6
CA at outlet, mol/l
Acknowledgements
t, min
3
0 (b)
4
8
12
16
20
t, min
Fig. 7. Output response of the PFR in the presence of disturbances in CA0 .
5. Conclusions In this paper, a feedback control method was proposed for first-order hyperbolic PDE systems. The proposed control method exploits the geometric characteristics of the PDE systems using the method of characteristics and is able to provide high-performance control for hyperbolic PDE systems. With the addition of output error integral term, the control method can tolerate numerical error and eliminate the output error caused by model-plant mismatch. Simulation results indicate that the proposed control method can produce satisfactory output response for the processes with accurate state estimation, as well as highly nonlinear processes with model-plant mismatch and inaccurate state estimation. It is shown that the proposed controller can effectively reject the effect of disturbances and maintain the processes at desirable operation. The proposed control method does not require complex computation. It can incorporate all the available measurement information to improve the control performance without the demanding requirement of infinite-dimensional observation. Therefore, it can easily be implemented in industrial practice. In this paper, focus has been on controller synthesis for hyperbolic DPS including quasi-linear PDE systems (both scalar PDEs and vector PDEs) and nonlinear PDE systems
Alvarez-Ramirez, J., 2001. Linear boundary control for a class of nonlinear pde processes. Systems and Control Letters 44, 395–403. Aris, R., 1969. Elementary Chemical Reactor Analysis, Prentice-Hall Inc., Englewood Cliffs, NJ. Arnold, V.I., 1988. Geometric Methods in the Theory of Ordinary Differential Equations, Springer, New York. Baker, J., Christofides, P.D., 2000. Finite-dimensional approximation and control of non-linear parabolic pde systems. International Journal of Control 73 (5), 439–756. Butkowskii, A.G., 1969. Distributed Parameter Systems, American Elsevier, Chakravarti, S., Ray, W.H., 1999. Boundary identification and control of distributed parameter systems using singular functions. Chemical Engineering Science 54, 1181–1204. Christofides, P.D., 2000. Nonlinear and Robust Control of PDE Systems, Birkhauser, Boston. Christofides, P.D., Daoutidis, P., 1996. Feedback control of hyperbolic pde systems. A.I.Ch.E. Journal 42 (11), 3063–3086. Christofides, P.D., Daoutidis, P., 1997. Finite-dimensional control of parabolic pde systems using approximate inertial manifolds. Journal of Mathematical Analysis and Applications 216 (2), 398–420. Christofides, P.D., Daoutidis, P., 1998. Robust control of hyperbolic pde systems. Chemical Engineering Science 53 (1), 85–105. Duchateau, P., Zachmann, D., 1989. Applied Partial Differential Equations, Harper & Row Publisher, New York. Gay, D.H., Ray, W.H., 1995. Identification and control of distributed parameter systems by means of the singular value decomposition. Chemical Engineering Science 50 (10), 1519–1539. Godasi, S., Karakas, A., Palazoglu, A., 2002. Control of nonlinear distributed parameter processes using symmetry groups and invariance conditions. Computers & Chemical Engineering 26 (7–8), 1023–1036. Gundepudi, P.K., Friedly, J.C., 1998. Velocity control of hyperbolic partial differential equation systems with single characteristic variable. Chemical Engineering Science 53 (24), 4055–4072. Hanczyc, E.M., Palazoglu, A., 1995a. Nonlinear control of a distributed parameter process: the case of multiple characteristics. Industrial Engineering & Chemistry Research 34, 4406–4412. Hanczyc, E.M., Palazoglu, A., 1995b. Sliding mode control of nonlinear distributed parameter chemical processes. Industrial Engineering & Chemistry Research 34, 557–566. Hernández, E., Arkun, Y., 1992. Design of pole placement controllers for systems described by two-dimensional models. Chemical Engineering Science 47 (2), 297–310.
980
H. Shang et al. / Chemical Engineering Science 60 (2005) 969 – 980
Isidori, A., 1995. Nonlinear Control Systems, Springer, Berlin. Karafyllis, I., Daoutidis, P., 2002. Control of hot spots in plug flow reactors. Computers & Chemical Engineering 26 (7–8), 1087–1094. Laroche, B., Martin, P., & Rouchon, P., 1998. Motion planning for a class of partial differential equations with boundary control. Proceedings of the 37th IEEE Conference on Design & Control. McOwen, R., 1996. Partial Differential Equations, Prentice-Hall Inc., Englewood Cliffs, NJ. Neittaanmaki, P., Tiba, D., 1994. Optimal control of Nonlinear Parabolic Systems: Theory, Algorithms and Applications, Marcel Dekker Inc., New York. Park, H.M., Cho, D.H., 1996. The use of the Karhunen–Loeve decomposition for the modeling of distributed parameter systems. Chemical Engineering Science 51 (1), 81–98.
Ray, W.H., 1978. Some recent applications of distributed parameter systems theory—a survey. Automatica 14, 281–287. Ray, W.H., 1981. Advanced Process Control, McGraw-Hill, New York. Sadek, I.S., Bokhari, M.A., 1998. Optimal control of a parabolic distributed parameter system via orthogonal polynomials. Optimal Control Applications & Methods 19, 205–213. Shirvani, M., Inagaki, M., Shimizu, T., 1995. Simplification study on dynamic models of distributed parameter systems. A.I.Ch.E. Journal 41 (12), 2658–2660. Sira-Ramirez, H., 1989. Distributed sliding mode control in systems described by quasilinear partial differential equations. Systems and Control Letters 13, 117–181.