2532
Ind. Eng. Chem. Res. 2005, 44, 2532-2548
Feedforward and Feedback Tracking Control of Nonlinear Diffusion-Convection-Reaction Systems Using Summability Methods† T. Meurer* and M. Zeitz Institut fu¨ r Systemdynamik und Regelungstechnik, Universita¨ t Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany
Motion planning and design of feedforward and feedback tracking control are studied for a tubular reactor modeled by nonlinear parabolic diffusion-convection-reaction equations. The approach is based on formal power series parametrizations of the system states and inputs, whereby the domain of convergence of the formal series solutions can be greatly extended using appropriate summation methods. As a result, feedforward controls can be determined for a wide range of trajectories and system parameters. Furthermore, on the basis of a reinterpretation of the formal power series parametrization, state feedback with asymptotic tracking control and an observer can be designed for the tubular reactor. On the other hand, following the 2-degrees-of-freedom approach, the feedforward control can be supplemented by a standard output feedback to obtain robust tracking control for tubular reactors with parameter sets reflecting the unpacked tubular reactor up to the real fixed-bed reactor. 1. Introduction The control of infinite-dimensional systems described by delay equations or partial differential equations (PDEs) is an area of research comprising both theoretical impacts and practical applications, stemming, e.g., from mechanics or chemical engineering. In particular, modeling of chemical processes leads to parabolic diffusion-convection-reaction equations (DCREs) due to the appearing diffusive and convective phenomena combined with strong nonlinearities in the reactive regime and spatial variations. Well-known examples concern tubular and fixed-bed reactors with multiple steady states and complex dynamical behavior.1 In addition to the stabilization problem, the increasing demands on product quality and process efficiency require advanced model-based tracking control concepts to force the output variables to follow some predefined trajectories. This, in particular, illustrates the importance of appropriate motion planning, i.e., the design of attainable and realizable output trajectories. Typical examples concern start-up, the realization of finite-time transitions between stationary operating points, or shutdown of a reactor. Although of practical interest, traditional approaches seem to provide only indirect and often complicated solutions to the tracking control problem. Model-based control design methods for DCREs are typically based on either early- or late-lumping approaches. In the early-lumping approach, the process model is approximated first, using, for instance, finitedifference or finite-element techniques, modal approaches,2,3 or Galerkin’s method and its modifica† The authors dedicate this paper to Prof. Harmon Ray in honor of his 65th birthday, his longstanding cooperation with the Stuttgart Institute of Prof. E. D. Gilles, as well as his great commitment in the course of a 25-year student exchange program with Universita¨t Stuttgart. * To whom correspondence should be addressed. Tel.: +49(0)711-685-6291. Fax: +49(0)711-685-6371. E-mail: meurer@ isr.uni-stuttgart.de.
tions,4,5 in order to proceed with the controller synthesis based on the approximation of the infinite-dimensional system by a finite-dimensional set of ordinary differential equations (ODEs) capturing the essential dynamics. Nevertheless, while useful for system analysis, these approximations usually lead to high-dimensional and often complex model-based feedback control structures, narrowing their applicability. Furthermore, an immense computational effort is often required for the determination of the finite-dimensional process models, such as, e.g., in the case of Galerkin’s method in conjunction with so-called (approximate) inertial manifolds.5 In the late-lumping approach, the distributed nature of the system is kept as long as possible and controller synthesis is performed using the infinite-dimensional process model.6 The main drawback of this theoretically appealing approach stems from the necessary approximation of the usually obtained infinite-dimensional control laws as required for implementation. Besides the purely numerical approaches such as model predictive control with its computational drawbacks, the main focus of controller design concerns the stabilization problem, whereby only a few analytical (late-lumping) approaches exist in the design of tracking control for linear distributed parameter systems (DPSs).7 To overcome this fact, the industrially widely used 2-degrees-of-freedom (2DOF) approach8 seems appealing, which amends stabilizing control (with an observer) by a feedforward loop to obtain a desired tracking behavior, as illustrated in Figure 1. This, on the one hand, introduces additional degrees of freedom, which can be exploited, e.g., for optimization, and, on the other hand, is well-known to greatly improve the control performance. Hence, the question arises of how to derive accurate model-based feedforward control for a given DPS. In this presentation, an implementable late-lumping feedforward control design methodology based on formal power series parametrizations9-14 in conjunction with recent results on summation methods for nonlinear
10.1021/ie0495729 CCC: $30.25 © 2005 American Chemical Society Published on Web 10/16/2004
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2533
[
]
is used, and the nonlinearities are summarized in
(
)
x2 1 + x2/σ φ(x) ) x2 BDa(1 - x1) exp + βx2w 1 + x2/σ Da × Le(1 - x1) exp
Figure 1. Block diagram of a 2-degrees-of-freedom control scheme with DPS Σ∞, feedforward control Σff, feedback control Σfb, observer Σobs, and trajectory generator Σd for output tracking y(t) f yd(t).
DCREs15,16 is proposed in order to extend the 2DOF approach to nonlinear DPS. This is exemplarily studied for the tracking control design for a tubular reactor model. The proposed approach directly accounts for the infinite-dimensional structure of the system and yields analytically a feedforward control to achieve the tracking performance. In addition, the formal power series approach can be applied within the early-lumping framework because a reinterpretation of the formal power series solution allows the derivation of a finitedimensional design model for the determination of asymptotic tracking control with an observer. 1.1. Tracking Control Problem for a Tubular Reactor. As a generic example, the pseudohomogeneous axial dispersion model of a tubular reactor in the case of an irreversible first-order liquid-phase reaction A f B,1,4 as depicted in Figure 2, is considered to introduce the novel approach using formal power series and summability methods. The PDEs for the tubular reactor in dimensionless form1,4 read as
∂x - Lx ) φ(x), z ∈ (0, 1), t > 0 ∂t
(1)
for the states x(z,t) :) [x1(z,t), x2(z,t)]T with boundary (BCs) and initial conditions (ICs)
E0x(0,t) ) -u(t), t > 0
(2)
E1x(1,t) ) 0
(3)
x(z,0) ) x0(z), z ∈ [0, 1]
(4)
The dimensionless states stem from a normalization by x1(z,t) ) [c0 - c(z,t)]/c0 and x2(z,t) ) [ϑ(z,t) - ϑ0]EA/ Rϑ02 with some constant concentration c0 and temperature ϑ0, e.g., a constant stationary profile, saturation concentration, or temperature maximum. The parameters EA and R represent the activation energy and general gas constant, respectively. The boundary inputs ui(t), i ) 1 and 2, correspond to similarly normalized inflow concentration and temperature. In the PDEs, the linear operator
L)
[
Le ∂2 ∂ - Le Pe1 ∂z2 ∂z
0
0
1 ∂2 ∂ - -β Pe2 ∂z2 ∂z
]
(5)
)
The operators in Danckwert’s BCs (2) and (3) take the form
E0 )
Figure 2. Tubular reactor of length L with first-order liquidphase reaction A f B.
(
(6)
[
1 ∂ -1 Pe1 ∂z 0
] [ ]
∂ 0 ∂z , E ) 1 1 ∂ ∂ -1 0 Pe2 ∂z ∂z 0
(7)
Here, Pe1 and Pe2 denote the Peclet numbers for mass and heat transfer, respectively, Le is the Lewis number, Da is the Damko¨hler number, B is the adiabatic temperature rise, and β and γ represent characteristic quantities for heat exchange and reaction.1 The value x2w ) (ϑw - ϑ0)EA/Rϑ02 denotes the normalized wall temperature, where x2w ) 0 if ϑ0 ) ϑw with constant ϑw. Furthermore, temporal and spatial coordinates are normalized according to t ) t′vz/LepL and z ) z′/L with gas velocity vz, bed porosity p, and reactor length L. The control problem concerns the design of the boundary inputs u(t) :) [u1(t), u2(t)]T in order to track the output
y(t) ) x(1,t)
(8)
along desired trajectories t f yd(t) :) [y1,d(t), y2,d(t)]T for the robust realization of, e.g., start-up, set-point changes, transient operation, or shutdown of the tubular reactor. For the solution of the control problem, formal power series in conjunction with suitable and powerful summability methods are utilized to design feedforward control as well as asymptotic feedback tracking control with an observer. The paper is organized as follows: the notion of formal power series parametrizability and its application to feedforward control design and motion planning are established in section 2 and illustrated for the considered tubular reactor model. In section 3, appropriate summability methods are introduced and utilized for convergence acceleration and extension of the domain of convergence for the derived feedforward control. These results serve as a basis for the derivation of a differentially flat finite-dimensional control and observer design model in section 4. Simulation results for the feedback-controlled tubular reactor are presented in section 5. Finally, some remarks conclude the paper. 2. Formal Power Series Parametrization The notion of formal power series parametrizability (FPSP) basically states that any system state and input can be parametrized via formal power series in terms of a parametrizing variable. It can be shown that quite general DCRE models follow this definition.17 This is similar to the notion of differential flatness, which was introduced for nonlinear finite-dimensional systems18 and extended to DPS.9,10,12,19,20 Note that no general definition of differential flatness exists for DPS, such that for parabolic DCREs the following definition of FPSP17 is introduced:
2534
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
Definition 1. The boundary controlled set of nonlinear DCREs
N
[
]
∂x(z,t) ∂x(z,t) ∂2x(z,t) )0 , x(z,t), , ∂t ∂z ∂z2
(9)
z ∈ (0, 1), t > 0 ∂x G0 x(0,t), (0,t), u0(t) ) 0, t > 0 ∂z ∂x G1 x(1,t), (1,t), u1(t) ) 0, t > 0 ∂z x(z,0) ) x0(z), z ∈ [0, 1]
[ [
(10)
] ]
(11) (12) (13)
defined on (z, t) ∈ [0, 1] × R+ with K-dimensional state vector x(z,t) ) [x1(z,t), ..., xK(z,t)]T is called formal power series parametrizable (FPSP), if there exists a formal power series ∞
xˆ (z,t) )
∑ xˆ n(t) sn(z)
n)0
with sn(z) an appropriate polynomial of degree n in z, which formally satisfies (9)-(13) and whose coefficients xˆ n(t) can be expressed in terms of a parametrizing function
[
y(t) ) H x(0,t),
∂x ∂x (0,t), x(1,t), (1,t), u0,1(t) ∂z ∂z
]
with dim y ) dim u0 + dim u1 and its time derivatives, i.e.,
xˆ n(t) ) Ψn[y(t), y3 (t), ..., y(rx)(t)] where rx ∈ N and possibly rx f ∞. The system inputs u0,1(t) are called FPSP if (9)-(13) are FPSP and (11) and (12) can be (explicitly) solved for u0,1(t). In the sequel, it is shown that the tubular reactor model (1)-(8) satisfies the conditions of definition 1. 2.1. Parametrization of the Tubular Reactor Model. The formal power series ansatz ∞
xˆ (z,t) )
xˆ n(t) (1 - z)n ∑ n)0
(14)
with xˆ n(t) :) [xˆ 1,n(t), xˆ 2,n(t)]T is formally [series products and interchange of summation and differentiation are valid operations if the series is uniformly convergent because the set of uniformly convergent power series forms a differential algebra with respect to t and z;21 hence, the following operations on the series (14) are only applied formally and validated in sections 2.2 and 3] substituted into the PDEs (1), which yields
∞
∑
n)0
[
]
(n + 2)(n + 1) xˆ 1,n+2 - (n + 1)xˆ 1,n+1 ∞ Le dt Pe1 n (1 z) ) φ ˆ ( xˆ n(1 - z)n) dxˆ 2,n (n + 2)(n + 1) n)0 xˆ 2,n+2 - (n + 1)xˆ 2,n+1 + βxˆ 2,n dt Pe2 1 dxˆ 1,n
-
∑
where φˆ (x) ) φ1(x)/Le. Under the assumption that the nonlinear source function φˆ (x) can be expanded into a formal power series in (1 - z) ∞
φˆ [
∑ n)0
∞
xˆ n(1 - z)n] )
φˆ n(xˆ n+1)(1 - z)n ∑ n)0
(15)
with
xˆ n+1 :) [xˆ 0T, xˆ 1T, ..., xˆ n+1T]T
[
terms of equal order in (1 - z) can be sorted
∞
∑
n)0
(16)
]
(n + 2)(n + 1) xˆ 1,n+2 - (n + 1)xˆ 1,n+1 - φˆ 1,n(xˆ n+1) Le dt Pe1 (1 - z)n ) 0 dxˆ 2,n (n + 2)(n + 1) xˆ 2,n+2 - (n + 1)xˆ 2,n+1 + βxˆ 2,n - φˆ 2,n(xˆ n+1) dt Pe2 1 dxˆ 1,n
-
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2535
Because this equation has to hold for any z ∈ (0, 1), it follows directly for any n ∈ N0
1 dxˆ 1,n (n + 2)(n + 1) xˆ 1,n+2 Le dt Pe1 ˆ 1,n(xˆ n+1) ) 0 (17) (n + 1)xˆ 1,n+1 - φ
By solving for the coefficients xˆ 1,n+2 and xˆ 2,n+2, (17) and (18) allow an interpretation as differential recursions
×
(n + 2)(n + 1) 1 dxˆ 1,n - (n + 1)xˆ 1,n+1 - φ ˆ 1,n(xˆ n+1) (19) Le dt
[
xˆ 2,n+2 )
]
Pe2
× (n + 2)(n + 1) dxˆ 2,n - (n + 1)xˆ 2,n+1 + βxˆ 2,n - φ ˆ 2,n(xˆ n+1) (20) dt
[
]
Their unique solution requires the determination of two starting conditions xˆ 0(t) and xˆ 1(t). Therefore, consider first the homogeneous BCs (3) with (7) at the outlet of the tubular reactor which evaluate to T
xˆ 1(t) ) [xˆ 1,1(t), xˆ 2,1(t)] ) 0
(21)
using the formal ansatz (14). Because as a result of the unknown inputs u(t) no further information is provided by BC (2), the output equation (8) is imposed as BC (14)
x(1,t) ) y(t) ) [xˆ 1,0(t), xˆ 2,0(t)]T ) xˆ 0(t)
(22)
Under the assumption of sufficiently smooth functions y(t) ∈ C∞, the differential recursions (19) and (20) with starting conditions (21) and (22) are uniquely solvable by sequential evaluation. Note that the complexity of the recursions prevents a closed-form solution; nevertheless, their structure allows the extraction of fundamental information. As proven in lemma 1 of appendix A, the solution of the recursion satisfies
(
xˆ i,2n ) Ψi,2n xˆ ′2n-2,
xi,2n+1 ) Ψ ˆ i,2n+1(y,y3 ,...,y(n))
(28)
Hence, evaluation of the formal power series ansatz with (27) and (28) yields
xˆ (z,t) )
ˆ 2,n(xˆ n+1) ) 0 (18) (n + 1)xˆ 2,n+1 + βxˆ 2,n - φ
Pe1
(27)
∞
dxˆ 2,n (n + 2)(n + 1) xˆ 2,n+2 dt Pe2
xˆ 1,n+2 )
xˆ i,2n ) Ψ ˆ i,2n(y,y3 ,...,y(n-1),y(n) i )
)
dxˆ i,2n-2 dt
xˆ i,2n+1 ) Ψi,2n+1(xˆ ′2n)
(23) (24)
for i ) 1 and 2 and any n ∈ N where
xˆ ′2j :) [xˆ ′1,2jT, xˆ ′2,2jT]T
(25)
xˆ ′i,2j :) [xˆ i,0, xˆ i,2, ..., xˆ i,2j]T, i ) 1 and 2
(26)
This illustrates that the coefficients with odd indexes depend nonlinearly but algebraically on the coefficients with even and smaller indexes. Furthermore, because xˆ 0 ) y and xˆ 1 ) 0, it follows from (23) and (24) that for i ) 1 and 2
Ψ ˆ 2n(y,y3 ,...,y(n)) (1 - z)n + ∑ n)0 ∞
∑ Ψˆ 2n+1(y,y3 ,...,y(n)) (1 - z)n n)0 ) Ψ1(1-z,y,y3 ,...,y(n),...)
(29)
For the parametrization of the inputs u(t), formally evaluate BC (2) with (7) and use the power series ansatz (14) for i ) 1 and 2
[
∞
ui(t) )
∑ xˆ i,n(t) + n)0
n+1 xˆ i,n+1(t) Pei
]
(30)
Substitution of the parametrized series coefficients (27) and (28) into (30) provides the desired expression of the input functions u(t) in terms of the parametrizing function y(t) and its time derivatives up to infinite order, which reads after rearrangement of the series coefficients ∞
ui(t) ) yi(t) +
[
] ∑[ ]
∑ 1+ n)0
2n + 1 Pei
∞
1+
n)1
2n
Pei
Ψ ˆ i,2n+1(y,...,y(n)) +
Ψ ˆ i,2n(y,...,y(n-1),y(n) i ) (31)
Because dim y ) dim u, the conditions of definition 1 for FPSP are satisfied for the tubular reactor model (1)-(8) under assumption (15), with the parametrizing function y(t) ) x(1,t) corresponding to the output equation (8). As can be clearly seen from the parametrizations (27) and (28), series convergence greatly depends on an appropriate choice of trajectories for y(t). Hence, if convergence can be ensured, (31) immediately provides the open-loop control ud(t) :) [u1,d(t), u2,d(t)]T, which is required to track the prescribed output trajectories y(t) f yd(t) in the nominal case. This fact is studied in the following. Before doing so, assumption (15) has to be justified for the tubular reactor model. Note that this imposed restriction is rather weak because many types of nonlinearities can at least within a subset of their range of definition be approximated by a formal power series of the above form.14 For the vector φ(x) of nonlinearities (6), the following expansion of the exponential term holds locally, with constants ai, i ∈ {0, 1, ..., I} adjusted for a given range of operation for x2(z,t). Assumption 1. There exist constants ai ∈ R, i ∈ {0, 1, ..., I}, s.t.
exp
(
x2
1 + x2/σ
)
I
≈
aix2i ∑ i)0
at least locally within some range of values of x2.
(32)
2536
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
Consequently, evaluation of (6) with (32) provides
φ(x) ≈
[
]
I
[ ]
0 Da × Le (1 - x1) aix2i + βx BDa 2w i)0
∑
(33)
which can be transferred in the required form (15) by successive application of the Cauchy product formula. As outlined above, uniform convergence of the formal power series solution (29) and hence the validation of the so far formal operations such as the interchange of differentiation and summation greatly depends on appropriate motion planning. Remark 1. The proposed approach for feedforward control design using FPSP can be easily extended to several scenarios involving input and output delays: (i) In the case of input delays, i.e., when the input is given as u(t-∆tu), the input parametrization (31) can be evaluated at t + ∆tu, which directly provides the compensation of the delay ∆tu. (ii) In the case of a product recycle loop,22,23 a ratio r of the product stream at the outlet of the reactor is fed back time-delayed into the reactor inlet, such that BC (2) reads E0x(0,t) ) -(1 - r)ur(t) - ry(t-∆ty). Because the reactor output y(t) denotes the parametrizing function, simply choose ur(t) ) [u(t) - ry(t-∆ty)]/(1 - r) with u(t) from (31) to compensate the influence of the delayed output due to the recycle loop. Remark 2. For computation of the series coefficients xˆ n, n g 0, by sequential evaluation of the differential recursions (19)-(22), efficient symbolic algorithms were developed and implemented in the computer algebra system MATHEMATICA.17 As a result, analytical expressions for the series coefficients are obtained, which can be used for evaluation of the formal power series. 2.2. Feedforward Control and Motion Planning. To obtain a feedforward boundary control for the solution of the tracking control problem, the derived formal power series solution (29) and hence (31) requires the specification of smooth desired trajectories yd(t) :) [y1,d(t), y2,d(t)]T ∈ C∞ for the output or parametrizing function y(t), respectively. In addition, uniform convergence of the formal solution has to be ensured by a suitable choice of yd(t). Note that, because z ∈ [0, 1], it is sufficient to verify at least unit radii of convergence. In the following, functions of a certain Gevrey order10,12,19,24,25 are used for motion planning. Definition 2. A function y(t) is called Gevrey of order R if the bounds
sup |y(n)(t)| e M(n!)R/Rn t∈R+
(34)
hold ∀ n ∈ N0 and suitable constants M, R ∈ R+. Convergence results are available for some simplified versions of the tubular reactor model (1)-(8): (i) Several scalar PDEs with temporally and spatially varying as well as state-dependent coefficients are studied by Lynch and Rudolph.12 (ii) Convergence conditions are proven for a scalar nonlinear PDE with cubic nonlinearity by Meurer and Zeitz.26 (iii) Coupled sets of PDEs modeling a tubular reactor with quadratic nonlinearity are considered by Lynch and Rudolph11 or Meurer et al.13 For a comprehensive study, the reader is also referred to the work of Rudolph et al.27-29 and results for a nonlinear Stefan problem.30
These results share the fact that convergence greatly depends on the constants M and R of the assigned output trajectories y(t) satisfying (34) and the system parameters. Furthermore, the Gevrey order is limited to R ∈ (1, 2). On the other hand, the chosen formal power series parametrization directly provides a solution to the considered control problem of finite-time transitions between an arbitrary number IS of stationary profiles xiS(z), i ) 0, 1, ..., IS. For the tubular reactor model (1)(8), stationary profiles xS(z) are governed by the following boundary value problem:
LxiS(z) + φ(xiS(z)) ) 0, z ∈ (0, 1)
(35)
E1xiS(1) ) 0
(36)
xiS(1) ) yiS
(37)
Note that the BCs (2) at the reactor inlet are replaced by (37) as the imposed BCs (22) for the determination of the stationary profiles xiS(z). The respective stationary inputs uiS can hence be directly computed by evaluating the inhomogeneous BCs (2) for t f ∞. Variations in yiS result in different stationary profiles xiS(z). As a representative scenario for tracking control for the tubular reactor, it is assumed to start up the reactor from the initial profile x0(z) ) x0S(z), to reach an operation profile x1S(z) along a predefined trajectory in finite time, and to switch to new stationary profiles x2S(z) and x3S(z) along suitable paths. Because of (35)(37), the stationary profiles xiS(z) are parametrized by yiS, i ) 0, 1, 2, and 3, such that appropriate motion planning is required to design a desired trajectory yd(t) for the parametrizing function connecting yiS, i ) 0, 1, and 2, along predefined paths. To meet these objectives (34), following work by Fliess et al.,19 the C∞ function Φγ,T: R+ f R defined as
Φγ,T(t) )
is used, where
{
0 1
∫0tφγ,T(τ) dτ ∫0Tφγ,T(τ) dτ
{ ([
if t e 0 if t g T if t ∈ (0, T)
if t * (0, T) -1 φγ,T(t) ) exp t t γ if t ∈ (0, T) 1TT 0
(
)]
)
The function Φγ,T(t) is nonanalytical at t ∈{0, T}; i.e., (n) for n > 0 the derivatives satisfy Φγ,T (t) ) 0 for t e 0 or (n) t g T but Φγ,T(t) * 0 for t ∈(0, T). This property allows one to start exactly at Φγ,T(0) ) 0 and reach exactly the final value Φγ,T(t g T) ) 1 within the finite-time interval t ∈ [0, T]. Further properties of Φγ,T(t) are derived by Lynch and Rudolph,12 where, in particular, its Gevrey order is determined as R ) 1 + 1/γ. Note that there is no solution to the control problem using an analytical function because an analytical function being constant
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2537
Figure 3. Desired trajectories yj,d(t) for yj(t), j ) 1 and 2, to realize transitions between four operating points parametrized by yiS ) xiS(1), i ) 0, 1, 2, and 3, with parameters summarized in appendix B.
on an interval of its domain is necessarily constant. Hence, no transition is possible in this case. When these results are combined for Φγ,T(t), the desired trajectory yd(t) ) [y1,d(t), y2,d(t)]T for the parametrizing output y(t) is assigned as 0 1 0 yj,d(t) ) yj,S + (yj,S - yj,S )Φγ 1,T 1-T 0(t-T0j ) + j
j
j
2 1 - yj,S )Φγ 2,T 3-T 2(t-T2j ) (38) (yj,S j
j
j
for j ) 1 and 2. This approach can be easily generalized to an arbitrary number of transitions j
IS
yj,d(t) )
0 yj,S
+
i+1 i (yj,S - yj,S )Φγ ∑ i)0
i+1,T 2i+1-T 2i j j
j
(t-T2i j ) (39)
between stationary profiles. Note that the trajectories for y1,d(t) and y2,d(t) can be chosen independently, such that depending on the transition times Tij, i ) 0, 1, ..., 2 ISj, j ) 1 and 2, a maximum number of 1 + 2∑j)1 ISj operating profiles can be realized by (39). In addition, this approach for motion planning can be suitably modified to realize finite-time transitions between arbitrary, not necessarily stationary, profiles following the ideas of Laroche et al.10 In the sequel, the trajectories yd(t), as illustrated in Figure 3, are used for feedforward control design by evaluating (31) with yd(t) substituted for y(t). Considering tubular reactor models, the available convergence conditions11-13,26 are too restrictive for most physical situations;26 e.g., convergence for a nonisothermal tubular reactor is verified for Peclet numbers Pei < 2, i ) 1 and 2.11 Furthermore, the proof of convergency does not account for the respective speed of convergence, which might be very low for convectiondominated systems.16 On the other hand, powerful summation methods exist, which accelerate convergence of a formal series and under some circumstances also allow one to sum divergent series to a meaningful limit.32 This somehow surprising fact is explained in the next section and exploited for feedforward control design. 3. Summability Methods The solution x(z,t) of (1)-(4) is an element of the space U of real-valued functions in z and t. However,
the derived formal solution xˆ (z,t) is given as a formal power series (14) with coefficients in U. For the deduction of the solution from the formal power series, a map must be found that for a given formal power series with coefficients in U unambiguously assigns a function in U. These maps are called summation methods S, and the set of formal power series to which an element of U can thus be assigned is called the set of Usummable series. The most commonly used summation method is the limit of the partial sums as the number of coefficients in the partial sum tends to infinity, i.e., ∞
Sp(
N
lim ∑ xˆ n(t) (1 - z)n ∑ xˆ n(t) (1 - z)n) ) Nf∞ n)0 n)0
(40)
However, because power series must converge in order to assign to them an element of U, the space of series that can be summed via the limit of partial sums is relatively small. To enlarge the space of summable series, other methods have been introduced e.g., by Borel31 or Hardy.32 Summation methods and the corresponding sets of summable series must have certain properties in order to be useful for the solution of the control problem (1)-(8): (i) The summation method must be regular; i.e., it must assign to each convergent series the limit of its partial sums. (ii) The set of summable series must have the structure of a differential algebra in view of differentiation with respect to z. (iii) The formal derivative with respect to t of the formal power series solution must be uniformly summable in t. One summation method that ensures the first two properties is the so-called k-summation,21,33-35 which will be defined for scalar formal power series. Definition 3. The k-sum Sk of the (scalar) formal ∞ power series xˆ (z,t) ) ∑n)0 xˆ n(t) (1 - z)n is defined as ∞
Sk(xˆ (z,t)) ) lim ξf∞
∑n
sn(z,t) ∞
ξn Γ(1+nk) ξn
(41)
∑n Γ(1+nk)
n where sn(z,t) ) ∑j)0 xˆ j(t) (1 - z)j denotes the nth partial sum if the limit on the right-hand side of (41) exists. Whether property (iii) is fulfilled depends on the series under consideration and must be individually verified. From a theoretical point of view, the proof of k-summability or uniform k-summability is most interesting but, on the other hand, is so far only available for formal power series solutions to the Cauchy problem (infinite spatial extension) for linear PDEs.35,36 Hence, from an application point of view and because of the facts that (i) the limiting processes in (41) cannot be evaluated numerically and (ii) only a finite number of series coefficients can be computed as a solution to the differential recursions (19)-(22), an approximation of (41) must be defined:16
2538
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
Figure 4. Comparison of summation results using partial summation (40) and (N,ξ)-approximate k summation (42) for feedforward controls uj,d(t), j ) 1 and 2 in (31), with desired trajectory (38) for two sets of Pe numbers, Le ) 1, and parameters from Table 1 of appendix B.
Definition 4. The (N,ξ)-approximate k-sum SN,ξ,k of a k-summable (scalar) formal power series xˆ (z,t) ) ∑nxˆ n(t) (1 - z)n is defined as
ξn
N
SN,ξ,k(xˆ (z,t)) )
∑n sn(z,t) Γ(1+nk) N
ξ
n
chosen using an optimization approach. Because focus is set on tracking control, the optimization criteria is designed to minimize the time integral over the absolute tracking error on a given sufficiently large time horizon T, i.e.,
min J(y(t),yd(t),ud(t))
(42)
∑n Γ(1+nk) n where sn(z,t) ) ∑j)0 xˆ j(t) (1 - z)j denotes the nth partial sum. Numerical results obtained by applying SN,ξ,k elementwise on the formal power series parametrization of the inputs (31) compared to standard partial summation Sp are depicted in Figure 4 for Pe1 ) Pe2 ) 6 and Pe1 ) Pe2 ) 10. Obviously, in both cases, partial summation fails to sum the coefficients to a meaningful curvature. In view of these results, the (N,ξ)-approximate k-summation approach can be also interpreted as a filtering process, which allows one to resolve sufficiently accurate summation results for a given finite set of temporally and spatially varying series coefficients. The (N,ξ)-approximate k-sum (42) is a function of ξ and N, such that the quality of the approximation strongly depends on a suitable choice of both parameters. In general, the quality of the approximation increases with N, while for each given N, there exists an optimal value for the parameter ξ, possibly ξ f ∞.16 As a result, an optimization approach for the determination of the summation parameters seems promising, which is illustrated in the following section when applying SN,ξ,k elementwise on the determined parametrization of the tubular reactor states (29) and the corresponding inputs (31). 3.1. Choice of Summation Parameters. The summation parameters ξ :) [ξ1, ξ2] and k :) [k1, k2] are
ξ,k
(43)
where
J(y(t),yd(t),ud(t)) )
∫0T{y(τ;ud(τ)) - yd(t)}TQ{y(τ;ud(τ)) - yd(t)} dτ
with the positive-definite, symmetric weighting matrix Q ∈ R2×2. Here, y(τ;ud(τ)) denotes the solution of the feedforward-controlled tubular reactor model (1)-(4) with boundary inputs ud(‚) evaluated at time τ. Note that the optimization for the determination of the summation parameters ξ and k can be also performed when combining the feedforward control ud(t) with some feedback control, e.g., standard proportional-integralderivative (PID) control for stabilization and/or compensation of model errors and exogenous disturbances. Furthermore, the chosen optimization approach also allows one to impose additional constraints, e.g., on the input values. For the explicit determination of the summation parameters ξ and k, (43) is evaluated by numerical simulations by applying the feedforward controls ud to a high-order semidiscretization of the tubular reactor model (1)-(8). 3.2. Feedforward Control of the Tubular Reactor. To illustrate the wide applicability of the proposed formal power series approach in conjunction with the introduced (N,ξ)-approximate k-summation, feedforward control (31) of tubular reactor model (1)-(8) is considered. Therefore, the desired trajectory (38) is substituted
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2539
Figure 5. Simulation results for feedforward tracking control (31) of the tubular reactor model (1)-(8) with Peclet and Lewis numbers varied. All other parameters can be found in appendix B. From left to right: boundary inputs (2); tracking behavior for x1(1,t) and x2(1,t).
into (31) to obtain the desired feedforward controls ui,d(t), i ) 1 and 2.
[ ] ∑[ ] ∞
ui,d(t) ) yi,d(t) +
∑
1+
n)0 ∞
1+
n)1
2n + 1 Pei
2n
Pei
Ψ ˆ i,2n+1(yd,...,y(n) d ) +
(n) Ψ ˆ i,2n(yd,...,y(n-1) ,yi,d ) (44) d
For the simulations, a high-order semidiscretized version of the DPS (“method of lines” approach) is used as a plant model. The model parameters are summarized in appendix B. Remark 3. For the application of the formal power series approach, the nonlinear source (6) has to be approximated following assumption 1. Therefore, similar to work by Lynch and Rudolph,11 a first-order Taylor linearization of φ˜ (x2) ) exp[x2/(1 + x2/σ)] around some
constant reference value x02 can be performed to determine the coefficients a0 and a1 of (32), i.e.,
[
a0 ) φ˜ (x02) 1 -
(
x02
1+
)]
x02 σ
2
, a1 )
φ˜ (x02) - a0 x02
Alternatively, the linearization can be evaluated along the desired trajectory y2,d(t), which leads to much higher accuracy but time-varying coefficients a0(t) and a1(t) and more complex differential recursions. Hence, for any of the following simulations, a0(t) and a1(t) were assumed to be slowly varying, such that dnai(t)/dtn ≈ 0 for n > 0, i ) 1 and 2. This has proven to be a good tradeoff between accuracy and computational complexity. Simulation results for this scenario are depicted in Figure 5 for various sets of model parameters Pe1, Pe2, and Le. The remaining model and summation param-
2540
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
eters are summarized in Table 1 of appendix B. Following Jensen and Ray,1 the chosen parameter sets correspond to an empty tube for Pe1 ) Pe2 and Le ) 1 (Figure 5a-c) and merge into a packed fixed-bed reactor (Figure 5d) for Pe1 ) 200, Pe2 ) 100, and Le ) 500. Note that, throughout the literature on control of tubular reactors, usually low and unrealistic Pe numbers are employed (among others, e.g., by Christofides5 or Antoniades and Christofides,22 as was already summarized in work by Panfilov and Sheintuch37). However, obviously the proposed approach for feedforward control design is applicable for wide parameter ranges. Initial results on the application of divergent power series solutions to boundary control design for linear DPS are available using least-term summation.10 Leastterm summation at first determines the absolutely smallest coefficient of the series and then calculates the partial sum up to this smallest term. Nevertheless, a comparison16 of (N,ξ)-approximate k-summation with least-term summation illustrates that the latter is often not applicable to nonlinear problems. Besides the nice tracking performance for the state x1(1,t), some deviations appear between the desired and obtained behavior for x2(1,t), which are due to the approximation of the nonlinear source term (6) reflecting the exothermic reaction. Nevertheless, under the chosen severe conditions, a sufficiently accurate tracking performance is obtained using the proposed feedforward strategy, whereby the compensation of model errors and exogenous disturbances or the stabilization of operating regimes requires a feedback strategy as depicted in Figure 1. 4. Feedback Tracking Control with Observer On the basis of a reinterpretation of the formal power series results (23)-(25) used for the feedforward control design in conjunction with appropriate summability methods, the feedback control design approach proposed by Meurer and Zeitz15 for scalar and weakly nonlinear DCREs is extended to the considered general nonlinear case. This design is based on the derivation of finitedimensional generalized nonlinear controller normal form approximating the infinite-dimensional tubular reactor. 4.1. Nonlinear Controller Normal Form. Because the coefficients xˆ 2n+1, n ∈ N0, with odd indexes can be expressed algebraically by the coefficients xˆ 2n, n ∈ N0, with even indexes following lemma 1 of appendix A, it suffices to consider only the equations (19) and (20) with even indexes, i.e.,
where
ϑ1,2n(xˆ ′2n) :) Le(2n + 1)Ψ1,2n+1(xˆ ′2n) + Leφˆ 1,2n(xˆ 2n+1(xˆ ′2n)) ϑ2,2n(xˆ ′2n) :) (2n + 1)Ψ2,2n+1(xˆ ′2n) + φ ˆ 2,2n(xˆ 2n+1(xˆ ′2n)) with xˆ 2n ′ ) [xˆ 0T, xˆ 2T, ..., xˆ 2nT]T and Ψ(‚) from (24). Recall the formal power series parametrization (31) of the inputs ∞
ui(t) )
[
∑ xˆ i,n(t) + n)0
In work by Meurer and Zeitz,14 it is shown that an appropriate flat finite-dimensional design model can be obtained under the formal assumption of an at least unit radius of convergence of (14). This formally allows a truncation of the set of equations (47) and (48) at some integer n ) N - 1, N > 1. The main idea is based on the interpretation of (47) and (48) for n ∈ {0, 1, ..., N - 1} as sets of ODEs for the coefficients of the formal power series with even indexes. As can be seen from (47) and (48), any ODE is characterized by a forward coupling due to the terms xˆ i,2n+2, i ) 1 and 2, such that by truncation at n ) N - 1 the two variables xˆ i,2N, i ) 1 and 2, remain in the ODEs for xˆ i,2N-2, i ) 1 and 2. Recall that, in the previous sections, the differential recursions for the formal power series parametrization are expressed in terms of xˆ i,2n+2, i ) 1 and 2, by imposing BCs (8) or (22), respectively, and formally neglecting the BCs (2) involving the inputs ui, i ) 1 and 2. However, under the formal assumption of convergence, a finite-dimensional design model can be derived by replacing xˆ i,2N, i ) 1 and 2, in (47) and (48) for n ) N - 1 and introducing the inputs using an approximation of (49), i.e.,
On the other hand, the convergence requirement can be replaced by the less restrictive assumption of ksummability.21 Hence, instead of truncating (49) at n ) N - 1, its (2N-1,ξ)-approximate k-sum (42) is introduced for i ) 1 and 2, i.e.,
1 dxˆ 1,2n (2n + 2)(2n + 1) ) xˆ 1,2n+2 + Le dt Pe1 ˆ 1,2n(xˆ 2n+1) (45) (2n + 1)xˆ 1,2n+1 + φ dxˆ 2,2n (2n + 2)(2n + 1) ) xˆ 2,2n+2 - βxˆ 2,2n + dt Pe2
2N-1
ui ≈
dxˆ 1,2n (2n + 2)(2n + 1)Le ) xˆ 1,2n+2 + ϑ1,2n(xˆ ′2n) (47) dt Pe1 dxˆ 2,2n (2n + 2)(2n + 1) ) xˆ 2,2n+2 - βxˆ 2,2n + ϑ2,2n(xˆ ′2n) dt Pe2Le (48)
ξni
n
∑ (∑uˆ i,j)Γ(1+nk )
n)0 j)0
2N-1
i
ξni
(50)
∑n Γ(1+nk ) i
ˆ 2,2n(xˆ 2n+1) (46) (2n + 1)xˆ 2,2n+1 + φ Using (24), it follows that
]
n+1 xˆ i,n+1(t) , i ) 1 and 2 (49) Pei
When (24) and (50) are utilized, it follows after some intermediate calculations that N-1 2N xˆ i,2N ) - Ai,2nxˆ i,2n + bi,N(ξi,ki) ui Pei n)0
∑
N-1
∑ Ai,2n+1Ψi,2n+1(xˆ ′2n),
n)0
i ) 1 and 2 (51)
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2541
Figure 6. Feedback control and observer design model (52)-(54) for the tubular reactor (1)-(8) in generalized nonlinear controller normal form (here ηN ) 2N - 1).
where
Ai,n )
n Pei
( )(
σi,n-1 + 1 +
n
Pei
2N-2
1+
∑ σi,j
j)n
)
2N-1
bi,N(ξi,ki) )
∑ σi,n
n)0
with σi,n ) Γ(1+(2N-1)ki)/Γ(1+nki)/ξi n ∈ N0, σi,-1 ) 0. This allows the introduction of inputs ui, 2N-1-n,
i ) 1 and 2, into the ODEs (47) and (48), which results in a finite-dimensional approximation of the DCRE by the ODEs (52) and (53) and output equations (54) as given in Figure 6. Simple comparison with (25) allows one to rewrite (52)-(54) in the following form
dxˆ ′2N-2 ) Axˆ ′2N-2 + T(xˆ ′2N-2) + Bu dt
(55)
y ) Cxˆ ′2N-2
(56)
2542
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
Figure 7. Eigenvalue distribution in the complex s-plane using S2N-1,ξ,k with ξ ) ξ1 ) ξ2 ∈ {13, 13.5, ..., 55}, k ) k1 ) k2 ) 1, and N ) 10, for the linear tubular reactor example (57). Here O denotes the exact eigenvalues from modal analysis; triangle right and triangle left correspond to the eigenvalues for ξ ) 13 and 55, respectively.
Figure 8. Eigenvalue distribution in the complex s-plane using S2N-1,ξ,k with ξ ) ξ1 ) ξ2 ∈ {800, 820, ..., 4000}, k ) k1 ) k2 ) 2, and N ) 10, for the linear tubular reactor example (57). Here O denotes the exact eigenvalues from modal analysis; triangle right and triangle left correspond to the eigenvalues for ξ ) 800 and 4000, respectively.
with matrices
A)
[
]
[
A1 0 b 0 , B) 1 b 0 A2 0 2
]
and the vector of nonlinearities T(xˆ ′2N-2) ) [T1(xˆ ′2N-2)T, T2(xˆ 2N-2 ′ )T]T. This input-affine nonlinear multiple-input multipleoutput (MIMO) system of finite dimension allows an interpretation as a generalized nonlinear controller normal form because of the obtained band structures
in both matrices A1 and A2 and the double (with respect to the two subsystems) triangular structure of T1(xˆ 2N-2 ′ ) and T2(xˆ 2N-2 ′ ). To illustrate that the ODEs (52) and (53) represent indeed a sufficient approximation of the system dynamics if the summation parameters ξ and k are chosen appropriately for given N, consider Figures 7-9. There, the eigenvalue spectrum of a simplified linear version of the tubular reactor model (1) given by
∂x - L*x ) φ* (x) ∂t
(57)
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2543
Figure 9. Eigenvalue distribution in the complex s-plane using S2N-1,ξ,k with ξ ) ξ1 ) ξ2 ∈ {110 000, 115 000, ..., 300 000}, k ) k1 ) k2 ) 3, and N ) 10, for the linear tubular reactor example (57). Here O denotes the spectrum of the operator; triangle right and triangle left correspond to the eigenvalues for ξ ) 110 000 and 300 000, respectively.
with
L* )
[
1 ∂2 Pe1 ∂z2
0
0
1 ∂2 +2 Pe1 ∂z2
]
[
-x - x2 , φ*(x) ) 2x 1 1
]
and BCs (2) and (3) as above is depicted and compared to the corresponding spectrum of the system matrix A of the design model (52) and (53) obtained through the previous derivations for different summation parameters ξ and k with N ) 10 and Pe1 ) Pe2 ) 1. In any case, the dominating and grouped eigenvalues of the operator L*x + φ*(x), which is unstable for the present parameters, can be approximated highly accurately for suitable choices of ξ and k in the k-summation. Furthermore, for k increasing from a value of 1 (Figure 7) to a value of 3 (Figure 9), also some of the higher modes can be reconstructed in the design model. Note that each subfigure within Figures 7-9 represents a zoom denoted by Z1-Z4 into the marked parts of the complex s-plane. Because inverse systems are, in general, ill-posed,38 degenerate eigenvalues appear in the spectrum of the design model, which tend into the open right-half complex s-plane along with increasing ξ for given k. Nevertheless, k-summation parameter combinations of ξ and k exist, which provide a suitable approximation of the governing DPS. This fact is exploited in the following section for the design of flatness-based tracking control with an observer. 4.2. Design of Tracking Control with Observer. Because of the structure of ODEs (52) and (53), it can be easily verified that each output yi, i ) 1 and 2, of (52)-(54) has a relative degree ri ) N, i ) 1 and 2. Hence, the system consisting of (52) and (53) is differentially flat,18 such that any state xˆ i,2n, i ) 1 and 2, n ∈ {0, 1, ..., N - 1}, and any input ui, i ) 1 and 2, can be expressed in terms of y(t) ) x(1,t) and its time derivatives. In the sequel, the flatness property is exploited
for the design of asymptotically stabilizing feedback tracking control for the tubular reactor (1)-(8). Because flat systems are exactly feedback linearizable,18 asymptotic tracking control can be designed by means of the linear control theory; for details, the reader is referred to work by Fliess et al.18 or Rothfuss et al.39 Following this approach, a static feedback law is obtained
u ) Ω[y,...,y(N-1),y(N) d (t)-Λ(e)]
(58)
where Λ(e) ) [Λ1(e1), Λ2(e2)]T, ei ) [ei, e˘ i, ..., e(N-1) ]T, i i ) 1 and 2, can be any type of feedback control asymptotically stabilizing the tracking error ei ) yi yi,d, i ) 1 and 2. In the sequel, PIDN control is considered N
Λi(ei) ) pi,0
∫0tei(τ) dτ + ∑pi,je(j-1) i
(59)
j)1
where the parameters pi,j, i ) 1 and 2, j ) 0, 1, ..., N, are assumed to be coefficients of a Hurwitz polynomial and can be determined by eigenvalue assignment λi,j, i ) 1 and 2, j ) 0, 1, ..., N, and comparison of coefficients N
λN+1 + pi,NλN i i + ... + pi,1λ + pi,0 )
(λ - λi,j) ∏ j)0
(60)
Because full state information is necessary for the implementation of the flatness-based feedback control (58) and (59), an observer is needed to estimate the nonmeasured states. For the sake of simplicity, it is assumed that the flat output y is measured such that observability analysis is equivalent to the proof of differential flatness.39 Note that more general cases can be treated similarly if observability is preserved. Following work by Rothfuss et al.,39 a nonlinear tracking
2544
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
Figure 10. Block diagram of the flatness-based boundary tracking control scheme with spatial profile estimation for DPS (1)-(8) with yd ) [ydT, ..., y(N),T ]T. d
observer can be designed based on the finite-dimensional model (52)-(54) or equivalently (55) and (56), i.e.,
dxˆ˜ ′2N-2 ) Axˆ˜ ′2N-2 + T(xˆ˜ ′2N-2) + Bu ˜ + dt L(t) [y - Cxˆ˜ ′2N-2] (61) 0 x˜ˆ ′2N-2(0) ) xˆ˜ ′ 2N-2
(62)
N-1
(λ˜ - λ˜ i,j) ∏ j)0
(63)
4.3. Spatial Profile Estimation. In addition to the realization of state feedback, the designed observer can be used for spatial profile estimation throughout the transition process.13,14 An estimate x˜ (z,t) of the spatial profile x(z,t) is obtained by an evaluation of the first 2N terms of formal power series ansatz (14), with the coefficients xˆ n(t) being replaced by their estimated counterparts N-1
x˜ (z,t) )
∑ [x˜ˆ ′2n + Ψˆ 2n+1(x˜ˆ ′2n)(1 - z)](1 - z)n
n)0
2N-1
)
∑ x˜ˆ n(1 - z)n
ξn x˜ˆ j(1 - z)j Γ(1+nk) n)0 j)0
2N-1 n
∑∑
x˜ (z,t) )
2N-1
ξn
(65)
∑ Γ(1+nk)
n)0
0 with suitable ICs xˆ ′ 2N-2 . The time-variable gain matrix L(t) ∈ R+ × R2N×2N can be determined based on a linearization of the observer error dynamics along the (t), which are known because of desired states xˆ 2N-2,d ′ the flatness property of (52) and (53).39 This allows the application of the MIMO time-variable Ackermann formula39 for the design of L(t), such that appropriate eigenvalues λ˜ i,j, i ) 1 and 2, j ) 0, 1, ..., N - 1, can be assigned for the observer error dynamics
λ˜ N ˜ i,Nλ˜ N ˜ i,1λ˜ + p˜ i,0 ) i + p i + ... + p
significantly improved by applying the (2N-1,ξ)-approximate k-sum (42) on (64) such that
(64)
n)0
Note that the results of the profile estimation can be
Finally, the block diagram depicted in Figure 10 is obtained for the realization of the proposed flatnessbased feedback control with an observer for the boundary tracking control of nonlinear DCREs. Note that all nonmeasured states are replaced by their estimated counterpart as indicated by the tilde on each affected quantity. The control loop is extended by the profile estimation (65), which provides additional insight in the process behavior (e.g., for monitoring purposes) without further measurement or design effort. 5. Simulation Results for Tubular Reactor Control On the basis of a simulation scenario, the applicability of the proposed tracking control strategy using formal power series is illustrated for the tubular reactor model. Therefore, at first, the derived flatness-based feedback control scheme of Figure 10 with state feedback, asymptotic tracking control, observer, and spatial profile estimation is applied to a method-of-lines, high-order discretized model of the tubular reactor. Moreover, the 2DOF approach combining the feedforward part with output feedback control, as schematically depicted in Figure 1, is considered alternatively. For numerical simulation, the controlled tubular reactor model (1)-(8) is simulated “in the loop”; i.e., a sampling time ∆tsample ) 0.005 (dimensionless) is explicitly taken into account. Hence, the controls are assumed to be constant for each simulation time period tsimi+1 - tsimi ) ∆tsample, and measurements for the evaluation of the observer are only available at discrete instances of time tsimi+1. This approach directly reflects a possible implementation and allows one to determine
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2545
Figure 11. Simulation results for feedback tracking control with observer (58)-(63) of the tubular reactor model (1)-(8) with Pe1 ) Pe2 ) 6, Le ) 1, and model error, i.e., Daplant ) 1.1Dadesign. Top, from left to right: boundary inputs (2), tracking behavior for x1(1,t) and x2(1,t). Bottom: comparison of estimated profiles x˜ i(z,tj), i ) 1 and 2, and simulated profiles xi(z,tj), i ) 1 and 2, at different instances of time tj, j ) 1, 2, 3, and 4. The remaining parameters are summarized in Table 1.
effects due to the discrete and time-delayed application of the controls. 5.1. Tracking Control with Observer. To study the robustness of the proposed flatness-based feedback tracking control scheme (58)-(62) when applied to a high-order method-of-line discretization of the tubular reactor, a scenario involving a model error in the Da number is illustrated in Figure 11. The Da number for the simulation model is assumed to differ by +10% from the Da number used for the derivation of the feedback control. This yields a more intense, exothermic reaction, because of the influence of Da on both PDEs (1). Nevertheless, the desired trajectory is tracked highly accurate by the flatness-based control scheme. For comparison purposes, simulation results for the feedforward-controlled tubular reactor model are depicted additionally. The design model (52)-(54) used for controller synthesis is of dimension 10, i.e., N ) 5. The eigenvalues for the feedback control (60) and observer (63) are assigned as follows:
λi,j ) -4 - j/2, j ) 0, 1, ..., N λ˜ i,j ) -16 - j, j ) 0, 1, ..., N - 1 with i ) 1 and 2. The results shown for profile estimation (Figure 11, bottom) clearly prove the applicability of the control concept for process monitoring. In particular, for low and medium Pe numbers (e10), the proposed flatness-based control scheme with state feedback and an observer provides robust trajectory tracking. Furthermore, when the feedforward tracking control results of Figure 5 are recalled, very accurate trajectory tracking can be obtained also for severe conditions (high Pe and Le numbers) using the derived feedforward control. Hence, as an alternative to state feedback, application of the 2DOF approach seems to be appealing, by supplementing the already determined
feedforward control for the realization of the desired tracking behavior by stabilizing output feedback control. 5.2. Feedforward Control with Output Feedback. Although accurate feedforward control results can be obtained using the proposed (N,ξ)-approximate ksum, stabilizing feedback is needed to account for model errors and exogenous disturbances. For its design, it suffices to stabilize the error e(z,t) :) x(z,t) - xd(z,t) between the actual and desired profiles because the feedforward control ud(t) from (44) drives the system at least in a neighborhood of xd(z,t). Note that xd(z,t) is known because of the formal power series parametrization (29), with y(t) replaced by the desired trajectory yd(t). As a result, the error DPS for the tubular reactor (1)-(8) reads for z ∈ (0, 1) and t > 0
∂e - Le ) φ(e + xd) - φ(xd) ∂t
(66)
E0e(0,t) ) -[u(t) - ud(t)], t > 0
(67)
E1e(1,t) ) 0
(68)
e(z,0) ) e0(z), z ∈ [0, 1]
(69)
with
whereby the right-hand side of (66) also reflects the difference between the exact (6) and approximated (33) nonlinear source term φ(x). Choosing the boundary control in Figure 1 as
u(t) ) ud(t) + ∆u(t)
(70)
the output feedback control ∆u(t) has to be designed in order to asymptotically stabilize the error e(z,t). Various methods exist using, e.g., full state feedback,6 domainaveraging,40 or backstepping.41 Although theoretically appealing, these approaches (if at all applicable) require
2546
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
Figure 12. Simulation results for feedforward and PI output feedback control of the tubular reactor model (1)-(8) with Pe1 ) 200, Pe2 ) 100, Le ) 500, and model error, i.e., Daplant ) 1.1Dadesign. Top: tracking behavior for x1(1,t) and x2(1,t). Bottom: comparison of feedforward ui,d(t) and feedback ratio ∆ui(t), i ) 1 and 2, in the applied controls (70). The remaining parameters are summarized in Table 1.
full or at least partial knowledge of the profile e(z,t) for all z ∈ [0, 1] and t ∈ R+, making the control law infinitedimensional and usually not realizable. To avoid these problems and to obtain an implementable approach, proportional-integral (PI) output feedback
∆ui(t) ) KPi ei(t) + KIi
∫0tei(τ) dτ
(71)
is applied with ei(t) ) yi(t) - yi,d(t) for i ) 1 and 2. For its design, either tuning rules such as the ZieglerNichols criteria or Lyapunov’s direct method can be applied. In the following, the gains KPi and KIi , i ) 1 and 2, are determined from numerical experiments using Ziegler-Nichols. Simulation results for a parameter set corresponding to a fixed-bed tubular reactor with Pe1 ) 200, Pe1 ) 100, and Le ) 500 when applying feedforward and PI output feedback control are depicted in Figure 12. To illustrate the robustness of this 2DOF concept, similar to section 5.1, a model error in the Da number is assumed with a deviation of +10% between simulation and control design model. Hence, the initial profile x0(z) also differs from the nominal profile used for feedforward control design. Nevertheless, simple PI output feedback is sufficient to supplement the feedforward control and stabilize the system along the desired output trajectory under these severe conditions. In Figure 12 (top), the outputs (8) are shown compared to their desired values (38). The applied controls (70) are split up into the feedforward part ui,d(t) and the PI output feedback part ∆ui(t), i ) 1 and 2 (Figure 12, bottom). Obviously, the feedback ratio is of magnitudes lower than the feedforward part, which is typical for the 2DOF approach, where the feedforward control ensures the desired tracking behavior and the stabilizing feedback control accounts for model errors and exogenous disturbances. Note that these results can be improved using more sophisticated output feedback control concepts but il-
lustrate the wide range of applicability of the proposed formal power series approach. 6. Conclusions This contribution introduces a novel method to tracking control design for nonlinear parabolic diffusionconvection-reaction systems and is exemplarily studied for a generic tubular reactor model. The approach is based on a formal power series parametrization of the system states in conjunction with a suitable k-summation technique and provides a completely model-based solution for motion planning and feedforward control design. As illustrated, sophisticated summation techniques greatly enhance the convergence behavior of the power series solution, whereby the proof of uniform summability is still an open question. Nevertheless, this systematic approach extends the applicability of the 2-degrees-of-freedom concept to nonlinear DCREs, where so far only a few analytical methods are available for the design of nonlinear dynamic feedforward controls. Following this concept, it suffices to supplement highly accurate feedforward control by simple feedback control to robustify the desired tracking behavior. Simulation studies for a fixed-bed tubular reactor in a highly convective regime confirm this result. Besides feedforward control design, a reinterpretation of the formal power series approach allows one to derive an approximation of the infinite-dimensional DCRE by a finite-dimensional system of ODEs, which serves as a design model for asymptotic state feedback tracking control with observer and spatial profile estimation. For moderate convection, simulation studies illustrate the performance and robustness of the state feedback control approach. As is typical for feedback control of (nonlinear) DPS based on a finite-dimensional approximation, theoretical verificiation of the stability of the combined system, i.e., DPS with flatness-based
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2547
feedback control and an observer, is a major problem, which still lacks its solution. In summary, FPSP in conjunction with the appropriate summation methods yields a mutual basis for motion planning, feedforward, and feedback tracking control design for nonlinear diffusion-convection-reaction systems in a wide range of operation.
The authors gratefully acknowledge financial support provided by the “Deutsche Forschungsgemeinschaft” (DFG).
Lemma 1. The coefficients of the differential recursion for the tubular reactor
× (n + 2)(n + 1) 1 dxˆ 1,n - (n + 1)xˆ 1,n+1 - φ ˆ 1,n(xˆ n+1) Le dt
[
]
Pe2
× (n + 2)(n + 1) dxˆ 2,n - (n + 1)xˆ 2,n+1 + βxˆ 2,n - φ ˆ 2,n(xˆ n+1) dt
[
]
xˆ i,0 ) yi(t), i ) 1, 2
ξ1
k1
ξ2
k2
1 1 1 1 1 500 1 500
60 25 60 25 25 35000 5 35000
1.5 2 1.5 2 2 4.9 1 4.9
1650 3900 1650 3900 3900 22 5 22
3.05 2.95 3.05 3 2.95 2.11 1 2.11
5
xˆ i,2j ) Ψi,2j
(
)
dxˆ i,2j-2 xˆ ′2j-2, dt
xˆ i,2j+1 ) Ψi,2j+1(xˆ ′2j)
(72) (73)
×
(2m + 2)(2m + 1)
[
1 dxˆ 1,2m ˆ 1,2m(xˆ 2m+1) - (2m + 1)xˆ 1,2m+1 - φ Le dt Pe1
]
×
(2m + 2)(2m + 1)
[
1 dxˆ 1,2m - (2m + 1)Ψ1,2m+1(xˆ ′2m) Le dt
φ ˆ 1,2m([xˆ 0T,0T,xˆ 2T,Ψ1,3T(xˆ ′2),...,Ψ1,2m+1T(xˆ ′2m)]T)]
)
dxˆ i,2m dt
) Ψ1,2m+2 xˆ ′2m,
Pe1
× (2m + 3)(2m + 2) 1 dxˆ 1,2m+1 - (2m + 2)xˆ 1,2m+2 - φˆ 1,2m+1(xˆ 2m+2) Le dt
[
[
1
Pe1 (2m + 3)(2m + 2) m
∑
Le l)0
∂Ψ1,2m+1 dxˆ 2l ∂xˆ 2l
dt
]
×
- (2m + 2)xˆ 1,2m+2 -
φˆ 1,2m+1([xˆ 0T,0T,xˆ 2T,Ψ1,3T(xˆ ′2),...,xˆ 2m+2T]T)]
Because the differential recursion is linear in the time derivative, (72) can be explicitly solved for any dxˆ i,2l/ ′ ,xˆ i,2l+2), i ) 1 and 2. Hence, the latter dt ) Θi,2l(xˆ 2l equation reduces to
Appendix B: Model and Simulation Parameters
for i ) 1 and 2 and j ∈ N where xˆ j ) [xˆ 0T, xˆ 1T, ..., xˆ jT]T and xˆ 2j ′ ) [xˆ 0T, xˆ 2T, ..., xˆ 2jT]T. Proof. The proof follows by induction. Assume that (72) and (73) hold for j ∈ 1, 2, ..., m. This implies for n ) 2m in the differential recursion
Pe1
Furthermore, it follows for n ) 2m + 1
Similar calculations for xˆ 2,2m+2 and xˆ 2,2m+3 complete the proof.
with n ∈ N0 satisfy
(
Le
6 10 1 6 10 100 6 100
xˆ 1,2m+3 ) Ψ1,2m+3(xˆ ′2m+2)
xˆ i,1 ) 0
)
Pe2
6 10 1 6 10 200 6 200
)
Pe1
xˆ 1,2m+2 )
Pe1
4
xˆ 1,2m+3 )
Appendix A: Lemma on Recursion
xˆ 2,n+2 )
figure
11 12
Acknowledgment
xˆ 1,n+2 )
Table 1. Summary of Varied Model Parameters and Summation Parameters in the Consecutive Figures
The dimensionless parameters of the tubular reactor model (1)-(54) are chosen as such: B ) 12.0, Da ) 0.11, β ) 1.5, σ ) 20.0, x2w ) 0. Peclet numbers Pe1 and Pe2 and Lewis number Le are varied and summarized in the respective sections. The desired trajectory (38) is parametrized by y01,S ) 0.1, y11,S ) 0.7, y21,S ) 0.5, y02,S ) 0, y12,S ) 2, y22,S ) 1, γ11 ) γ21 ) γ12 ) γ22 ) 1, T01 ) 0, T11 ) 6.5, T21 ) 14, T31 ) 18, T02 ) 0.5, T12 ) 4.5, T22 ) 10, and T32 ) 14. Table 1 summarizes the Peclet Pe1 and Pe2, Lewis numbers Le, and summation paramaters ξ1, k1, ξ2, and k2 used in the respective figures. Literature Cited (1) Jensen, K. F.; Ray, W. H. The bifurcation behavior of tubular reactors. Chem. Eng. Sci. 1982, 37 (2), 199-222. (2) Georgakis, C.; Aris, R.; Amundson, N. R. Studies in the control of tubular reactorssI/II/III. Chem. Eng. Sci. 1977, 32, 1359-1387. (3) Balas, M. J. Modal control of certain flexible dynamic systems. SIAM J. Control Optim. 1978, 16 (3), 450-462. (4) Ray, H. W. Advanced Process Control; McGraw-Hill: New York, 1981. (5) Christofides, P. D. Nonlinear and Robust Control of PDE Systems; Birkha¨user: Boston, 2001. (6) Curtain, R. F.; Zwart, H. An Introduction to Infinite Dimensional Linear Systems Theory; Springer-Verlag: New York, 1995.
2548
Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005
(7) Byrnes, C. I.; Lauko´, I. G.; Gillam, D. S.; Shubov, V. I. Output Regulation for Linear Distributed Parameter Systems. IEEE Trans. Autom. Control 2000, 45 (12), 2236-2252. (8) Horowitz, I. M. Synthesis of Feedback Systems; Academic Press: New York, 1963. (9) Laroche, B.; Martin, P.; Rouchon, P. Motion planning for a class of partial differential equations with boundary control. Proceedings of the 37th Conference on Dec. Control, Tampa, FL, pp 736-741. (10) Laroche, B.; Martin, P.; Rouchon, P. Motion planning for the heat equation. Int. J. Robust Nonlinear Control 2000, 10, 629643. (11) Lynch, A. F.; Rudolph, J. Flatness-based boundary control of coupled nonlinear PDEs modelling a tubular reactor. Proceedings of NOLTA 2000, Dresden, Germany, 2000; pp 641-644. (12) Lynch, A. F.; Rudolph, J. Flatness-based boundary control of a class of quasilinear parabolic distributed parameter systems. Int. J. Control 2002, 75 (15), 1219-1230. (13) Meurer, T.; Becker, J.; Zeitz, M. Flatness-based feedback tracking control of a distributed parameter tubular reactor model. Proceedings (CD-ROM) of the European Control Conference (ECC), Cambridge, U.K., 2003. (14) Meurer, T.; Zeitz, M. A novel design approach to flatnessbased feedback boundary control for nonlinear reaction-diffusion systems with distributed parameters. In New Trends in Nonlinear Dynamics and Control, LNCIS; Kang, W., Borges C., Xiao, M., Eds.; Springer-Verlag: New York, 2003; pp 221-236. (15) Meurer, T.; Zeitz, M. Flatness-based feedback control of diffusion-convection-reactions systems via k-summable power series. Preprints (CD-ROM) of NOLCOS 2004, Stuttgart, Germany, 2004; pp 191-196. (16) Wagner, M. O.; Meurer, T.; Zeitz, M. K-summable power series as a design tool for feedforward control of diffusionconvection-reaction systems. Preprints (CD-ROM) of NOLCOS 2004, Stuttgart, Germany, 2004; pp 149-154. (17) Meurer, T. Feedforward and feedback tracking control of diffusion-convection-reaction equations. Ph.D. Thesis, planned 2004. (18) Fliess, M.; Le´vine, J.; Martin, P.; Rouchon, P. Flatness and defect of nonlinear systems: introductory theory and examples. Int. J. Control 1995, 61, 1327-1361. (19) Fliess, M.; Mounier, H.; Rouchon, P.; Rudolph, J. Syste`mes line´aires sur les ope´rateurs de Mikusin´ski et commande d’une poutre flexible. ESAIM-PROC 1997, 2, 183-193. (20) Fliess, M.; Mounier, H.; Rouchon, P.; Rudolph, J. A distributed parameter approach to the control of a tubular reactor: a multi-variable case. Proceedings of the 37th Conference on Dec. Control, Tampa, FL, pp 736-741. (21) Balser, W. Formal power series and linear systems of meromorphic ordinary differential equations; Springer-Verlag: New York, 2000. (22) Antoniades, C.; Christofies, P. D. Studies on Nonlinear Dynamics and Control of a Tubular Reactor with Recycle. Nonlinear Anal. 2001, 47, 5933-5944. (23) Berezowski, M. Effect of the delay time on the generation of chaos in continuous systems. One-dimensional model. Twodimensional model - tubular chemical reactor with recycle. Chaos, Solitons Fractals 2001, 12, 83-89. (24) Gevrey, M. La nature analytique des solutions des e´quation aux de´rive´es partielles. Ann. Sci. Ecole Normale Superieure 1918, 25, 129-190.
(25) Hua, C.; Rodino, L. Generel theory of PDE and Gevrey classes. In General Theory of PDEs and Microlocal Analysis; MinYou, Q., Rodino, L., Eds.; Addison-Wesley: Reading, MA, 1996; pp 6-81. (26) Meurer, T.; Zeitz, M. Two-degree-of-freedom tracking control design for boundary controlled distributed parameter systems using formal power series. IEEE 43rd Conf. Dec. Control 2004, in press. (27) Rudolph, J. Flatness Based Control of Distributed Parameter Systems; Shaker-Verlag: 2003. (28) Rudolph, J.; Winkler, J.; Woittennek, F. Flatness Based Control of Distributed Parameter Systems: Examples and Computer Exercises from Various Technological Domains; ShakerVerlag: 2003. (29) Rudolph, J. Beitra¨ ge zur flachheitsbasierten Folgeregelung linearer und nichtlinearer Systeme endlicher und unendlicher Dimension; Shaker-Verlag: 2003. (30) Dunbar, W. B.; Petit, N.; Rouchon, P.; Martin, P. Motion Planning for a Nonlinear Stefan Problem. ESAIM-COCV 2003, 9, 275-296. (31) Borel, E. Lecons sur les se´ ries divergentes; GauthierVillars: Paris, 1928. (32) Hardy, G. H. Divergent series; Oxford University Press: Oxford, U.K., 1949. (33) Ramis, J. P. Le se´ries k-sommable et leurs applications. In Lecture Notes in Physics 126: Complex Analysis, Micrological Calculus and Relativistic Quantum Theory; Springer-Verlag: New York, 1980; pp 178-199. (34) Malgrange, B. Sommation des se´ries divergentes. Expo. Math. 1995, 13, 163-222. (35) Balser, W.; Kostov, V. Recent progress in the theory of formal solutions for ODE and PDE. Appl. Math. Comput. 2003, 141, 113-123. (36) Balser, W. Divergent solutions of the heat equation: on an article of Lutz, Miyake and Scha¨fke. Pac. J. Math. 1999, 188 (1), 53-63. (37) Panfilov, V.; Sheintuch, M. Control strategies for front stabilization in a tubular reactor model. AIChE J. 2001, 47 (1), 187-196. (38) Isakov, V. Inverse Problems for Partial Differential Equations; Springer-Verlag: New York, 1998. (39) Rothfuss, R.; Rudolph, J.; Zeitz, M. Flatness-based control of a nonlinear chemical reactor model. Automatica 1996, 32, 14331439. (40) Boskovic´, D. M.; Krstic´, M.; Liu, W. Boundary control of an unstable heat equation via measurement of domain-averaged temperature. IEEE Trans. Autom. Control 2001, 46 (12), 20222028. (41) Boskovic´, D. M.; Krstic´, M. Backstepping control of chemical tubular reactors. Comput. Chem. Eng. 2002, 26, 1077-1085.
Received for review May 19, 2004 Revised manuscript received July 29, 2004 Accepted August 6, 2004 IE0495729