Model Predictive Control for Quasilinear Hyperbolic Distributed

2140

Ind. Eng. Chem. Res. 2004, 43, 2140-2149

Model Predictive Control for Quasilinear Hyperbolic Distributed Parameter Systems Huilan Shang,*,† J. Fraser Forbes,‡ and Martin Guay§ School of Engineering, Laurentian University, Sudbury, Ontario, Canada P3E 2C6, Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6, Department of Chemical Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6

Distributed parameter systems (DPSs) constitute an important class of systems that include many industrial processes modeled by partial differential equations (PDEs). Model predictive control (MPC) techniques for DPSs have been relatively scarce because of the relative difficulty of solving PDEs. Conventional MPC approaches for DPSs have usually been developed on the basis of approximate lumped models. The resulting controllers can exhibit poor performance characteristics or require substantial on-line computation to ensure adequate performance. This paper presents a novel MPC scheme for output control of distributed parameter systems that is based on the method of characteristics. It is shown via simulation that the proposed approach can yield a high-performance controller with a comparatively small computational load. 1. Introduction The basic premise of model predictive control (MPC) is to use process models to forecast the effect of control actions on future process states and outputs. Current control actions are obtained by solving, at each sampling instant, a finite horizon open-loop optimal control problem, using the current state estimate or measurement as the initial condition. Because MPC can provide adequate performance, enforce stringent safety requirements, and ensure efficient material and energy recovery, it has become widely accepted in the field of process control. Successful applications have been found in the chemical industry since the original pioneering works by Culter and Ramaker1 and Richalet et al.2 An account of current trends in linear MPC is provided in a number of survey papers.3-7 The use of high-fidelity models in MPC is motivated by the possibility of performance enhancement resulting from improved forecasting, and this promise of performance has resulted in an active research effort in nonlinear model predictive control (NLMPC).8-12 To date, most MPC developments and applications have focused on the processes that are amenable to lumped parameter modeling. Many industrial processes are distributed parameter systems (DPSs), in which states vary in both time and space (e.g., fixed-bed-reactor, polymer-extrusion, and sheet-coating processes). Mathematical descriptions of such systems, usually obtained by applying conservation laws, often take the form of partial differential equations (PDEs). The most common approach to DPSs is to approximate the system by a set of linear lumped models and then design linear model-based control strategies using the approximation.13,14 When the processes have large spatial variations, the use of highfidelity PDE models in model-based control for DPSs can improve the control and constraint-handling per* To whom correspondence should be addressed. Fax: 1-705675-4862. E-mail: [email protected]. † Laurentian University. ‡ University of Alberta. § Queen’s University.

formance. Unfortunately, MPC techniques based on PDE models are relatively scarce, in part, because of the mathematical complexity arising from the partial differential equation models. Some researchers have addressed the design of MPC for distributed parameter systems;15-20 however, most of these strategies use an approximate ordinary differential equation (ODE) model or a discretization of the underlying PDEs into a system of ODEs, so that the existing MPC techniques can then be applied to these systems. In such an approach, highaccuracy control usually requires that the DPS be approximated with a large number of ODEs, thus dramatically increasing the complexity of the calculations. The computational cost associated with these lumping approaches can be prohibitive for some processes and might not provide an improvement in performance that warrants the additional computational cost. To reduce the computational cost for PDE systems, research has been conducted into the development of computationally efficient optimization algorithms, and approaches using advanced-order reduction techniques have been proposed.21,22 The PDEs modeling DPSs can be hyperbolic, parabolic, or elliptic. Hyperbolic PDEs represent the dynamics of a large number of industrial processes involved in convection with negligible diffusion effects (e.g., fluid heat exchangers, plug-flow reactors, and fiber spinlines). In contrast to parabolic PDEs, the hyperbolic spatial differential operator contains eigenmodes of the same, or nearly the same, amounts energy. As a result, the treatment of hyperbolic PDEs requires an analysis of the infinite-dimensional nature of the systems. Research has been conducted into the development of control approaches exploring the infinite-dimensional nature of hyperbolic systems.23-25 However, the use of MPC methods to control such systems has not been fully explored, and research is required to develop MPC schemes that exhibit higher levels of performance with better handling of constraints. In this paper, an MPC scheme that is capable of effectively driving a process output to its set point is presented for processes modeled by quasilinear hyperbolic PDEs. For these hyperbolic PDE systems, the

10.1021/ie030653z CCC: $27.50 © 2004 American Chemical Society Published on Web 03/31/2004

Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004 2141

method of characteristics provides a geometric way of viewing the solution structure and can help in providing insight into the future evolution of the process output.26,27 Therefore, the method of characteristics can be an effective tool in developing model predictive control algorithms for a class of DPSs. Further, the proposed approach, which is based on the method of characteristics, is shown to be able to provide improved control performance with a comparatively small on-line computational load. 2. Characteristics of Quasilinear Hyperbolic PDE Systems Any given hyperbolic PDE system can be transformed into an equivalent system that involves differentiation in only one direction, namely, along the characteristics. Then, a partial differential equation can be reduced to an ordinary differential equation along characteristic curves, which leads to a powerful method for solving hyperbolic PDEs.28-30 2.1. Single Characteristic. All first-order scalar PDEs have single characteristics, i.e., only one characteristic curve passes through each coordinate point. This characteristic property of scalar PDEs is also shared by some systems of PDEs. Here, for the convenience of illustration, scalar PDEs are used to demonstrate the method of characteristics for systems with single characteristics. Consider the quasilinear equation for a function υ(t,x1,...,xn) on Rn+1

∂υ

n

+

∂t

∑ i)1

ai(x,υ,u)

∂υ

) f(x,υ,u)

(1)

∂xi

where t ∈ R is the time; x ) [x1 ‚‚‚ xn] ∈ Rn is a point in the manifold; a1(x,υ,u), ..., an(x,υ,u), and f(x,υ,u) are continuous (C0) functions in x and υ; and u represents the parameters. Given υ(t,x) as the solution of eq 1, we consider the graph z ) υ(t,x). This graph has the normal vector

[

N0 ) -

∂υ ∂t

-

∂υ ∂x1

‚‚‚

-

∂υ ∂xn

1

]

(2)

‚‚‚

an(x,υ,u)

f(x,υ,u)] (3)

defines a vector field in Rn+1 that is tangent to the solution graph at each point. This vector field is called the characteristic vector field of eq 1. The integral curves of the characteristic vector field are called the characteristics of the quasilinear equation. The ordinary differential equation defined by the vector field ξ(t,x,υ,u) is called the characteristic equation, which, for eq 1, has the form

t˙ ) 1 x3 ) a(x,υ,u) υ˘ ) f(x,υ,u)

∂υ1 ∂υ1 + a1 ) f1(υ1,...,υq,u) ∂t ∂x ∂υ2 ∂υ2 + a2 ) f2(υ1,...,υq,u) ∂t ∂x l

(5)

∂υq ∂υq + aq ) fq(υ1,...,υq,u) ∂t ∂x where υ1, υ2, ..., υq, are distributed state variables; u is the parameter, a1, a2, ..., aq, are constants or functions of x; and f1, f2, ..., fq, are continuous functions. If a1 * a2 ... * aq, the system has q distinct characteristics. The characteristics for eq 4 are q curves determined by

characteristic C:

at the point (t0, x0, υ(t0,x0)). Let z0 ) υ(t0,x0). Then, 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 plane tangent to the graph of z ) υ(t,x) at the point (t0, x0, z0). We say that

ξ(t,x,υ,u) ) [1 a1(x,υ,u)

If the graph υ ) υ(t,x,u) is a smooth surface S that is a union of such characteristic curves, then the tangent plane contains the vector ξ(x,υ,u) at each point (t, x, υ); 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 defined by the initial conditions is nowhere tangent to the vector field), then 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. 2.2. Multiple Characteristics. The complexity of the method increases dramatically for a system of first-order PDEs.29 In such systems, one usually finds more than one characteristic direction at each point of the solution surface, thus leading to a set of the characteristic ODEs that are usually coupled. However, even with the increased complexity, it remains possible to use the method of characteristics to develop a relatively efficient solution technique for a given system of first-order PDEs, especially when the dimension of the PDE system is low. A general quasilinear system of first-order equations with q dependent variables υ1, υ2, ..., υq, and two independent variables t and x can be transformed into a “characteristic normal form”31

dx ) ai , dt

i ) 1, 2, ..., q

(6)

The variations of state variables υ1, ..., υq, can be described by the ordinary differential equations

dυ1 ) f1(υ1,...,υq,u) dt

along characteristic C1

dυ2 ) f2(υ1,...,υq,u) dt

along characteristic C2 l

(7)

dυq ) fq(υ1,...,υq,u) along Cq characteristic dt Thus, using the method of characteristics, the system of PDEs (eqs 5) is transformed into ODEs along the characteristic curves. In general, the characteristic ODEs along different characteristic curves are coupled. 3. Characteristic-Based MPC Design

(4)

In this section, an MPC scheme for output regulation is designed for some hyperbolic systems using the

2142 Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004

method of characteristics. Control design for systems with both scalar characteristics and multiple characteristics is addressed. 3.1. Scalar Case. For processes modeled by a scalar hyperbolic PDE, Cauchy characteristics exist, and the PDE model can be described by its characteristic ODEs, which describes the evolution of the state along the characteristic curves. The method of characteristics provides a convenient and efficient way to predict the future output from the current state-variable profile and can be used to develop a characteristic-based MPC (CBMPC) scheme. A general first-order quasilinear system of one spatial dimension can be written as

∂υ ∂υ + a(x,υ,u) ) f(x,υ,u) ∂t ∂x υ(x)0) ) υb

(8)

y ) υ(xout,t) where a(x,υ,u) and f(x,υ,u) are smooth functions, the boundary value υb can be constant or time-varying, and the output is the value of the distributed state variable υ at spatial point xout (xout is often located at the exit of a process). The PDE model in eq 7 possesses a characteristic vector field ξ ) [1 a(x,υ,u) f(x,υ,u)] that defines the characteristic ODEs, similarly to eq 5. Then, a prediction of the future output can be obtained by integrating this characteristic ODE. The formulation of CBMPC might differ depending on the complexity of the PDE model. If the PDE model is linear, a(x,υ,u) is a constant, a, and f(x,υ,u) is linear in υ and u, i.e., f(x,υ,u) ) bυ + cu. Then, passing through an initial point (t0, x0, υ0), the characteristics of the PDE model can be obtained analytically as

where ∆t1 ) ts, ∆t2 ) 2ts, ..., ∆tm ) mts, and u0, u1, ..., umc-1, are the control moves in the next mc sample instants. The control beyond mc steps ahead is set to umc-1 ) umc ) ... ) um-1. By defining the vectors

yˆ ) [y(t0+∆t1)

u ) [u0

(9)

y(t0+∆ti) ) eb∆tiυ0(xm-i+1) + j-1)

c - eb(∆ti-∆tj)) uj-1 i ) 1, ..., mc b

[

c × b b∆t1 e -1 eb∆t2 - eb(∆t2-∆t1) l

c

b(∆ti-∆tj)

-e

(eb(∆ti-∆tmc-1) - 1) umc-1, b

‚‚‚ 0

eb(∆t2-∆t1) - 1

‚‚‚ 0

]

(12)

eq 10 can be expressed in a compact form as

yˆ ) y0 + Su

(13)

This equation gives an exact prediction of future output for the linear PDE systems in the absence of disturbances. Given the specified sampling time and the current state-variable profile at the spatial nodes, y0 and S can be determined in eq 13. The future control sequence u is calculated to make yˆ approach, as closely as possible, the output set point yr via the solution of the optimization problem u

(14)

With m > mc, the control moves can be computed using the least-squares solution

u ) (STS)-1ST(yr - y0)

(15)

Then, the first element of u is implemented. In this control law, S can be calculated off-line and y0 needs to be updated according to the available measurements using eq 11. This is the linear CBMPC form. For a general quasilinear PDE system, the CBMPC design is more complicated in comparison with that for linear systems. The characteristic ODEs are nonlinear, and numerical integration is usually required to obtain the sampled future output from the value of the current state variable at the spatial nodes. Assuming that the value of the current state variable at the spatial node x0 is υ0(x0), the output at a future time instant can be obtained by simultaneously integrating the characteristic ODEs of the PDE model in eqs 8 with initial conditions t(0) ) t0, x(0) ) x0, and υ(0) ) υ0(x0)

t ) t0 + ∆t

mc-1

(e ∑ j)1

0

eb∆tmc - eb(∆tmc-∆t1) eb(∆tmc-∆t1) - eb(∆tmc-∆t2) ‚‚‚ eb(∆tmc-∆tmc-1) - 1 ... eb∆tm - eb(∆tm-∆t1) eb(∆tm-∆t1) - eb(∆tm-∆t2) ‚‚‚ eb(∆tm-∆tmc-1) - 1

y(t0+∆ti) ) eb∆tiυ0(xm-i+1) + b(∆ti-∆tj-1)

umc-1]T

S)

)

These equations indicate that the state variable υ evolves along each characteristic curve and that no coupling between characteristic curves is involved. Therefore, the variation of υ along a particular characteristic curve is determined by the initial condition on the same characteristic curve. Discretizing the space into m points with a spacing ∆x and choosing a prediction sampling time ts ) ∆x/a, the output y ) υ(xout) at the future sample instants can be formulated from eq 9 as

i

‚‚‚

u1

and the matrix

c c υ ) υ0 + u eb∆t - u b b

i

y(t0+∆tm)]T

y0 ) [eb∆t1v0(xm) eb∆t2υ0(xm-1) ‚‚‚ eb∆tmυ0(x1)]T (11)

u

x ) x0 + a∆t

(eb(∆t -∆t ∑ j)1

‚‚‚

min||yˆ - yr||2 ) min||y0 + Su - yr||2

t ) t0 + ∆t

(

y(t0+∆t2)

c

) uj-1 + b

i ) mc + 1, ..., m (10)

x)

∫tta(υ,u,x) dτ ) φx(υ0,u,x0,∆t)

υ)

∫ttf(υ,u, x) dτ ) φυ(υ0,u,x0,∆t)

0

(16)

0

where ∆t ) t - t0. Integration can proceed until x

Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004 2143

reaches xout, at which point the corresponding statevariable value is the output at time t. The output at a different future time can be obtained by varying the initial point x0. The prediction horizon time can be taken to be equal to the residence time, given that the current value of the state variable and control action can only affect the output within one residence time. This leads to a prediction horizon p equal to the number of initial discretization points m. The calculation of the CBMPC can be performed by solving an optimization problem in eq 14 or

min J ) min {(y - yˆ ) Q(y - yˆ ) + ∆u R∆u} (17) r

∆u

T

r

T

∆u

subject to eqs 16 and y ) υ(xout,t), where yr is the desired future output trajectory and Q and R are two symmetric positive-definite weighting matrices. The calculation of the control actions using eq 17 requires solution of an optimization problem subject to nonlinear integral equation constraints. In the above CBMPC design procedure, numerical integration of the characteristic ODEs is decoupled for each spatial point and prediction time instant. This decoupling property allows one to predict the future output at relatively high accuracy with a comparatively small computational load. For general quasilinear systems, the prediction horizon p and control horizon mc can significantly affect the on-line computational demand of the proposed CBMPC. It is reasonable to choose a prediction horizon time equal to the residence time (i.e., p ) m). The control horizon can be chosen to be a small number to ensure small dimensions for S and yˆ 0 and also to favor stable and smooth operation.32 3.2. Multivariate Systems. For many problems encountered in practice, multiple dependent variables are required to model the physical systems. In such cases, the integral manifold of the PDEs cannot be obtained by the solution of unique characteristics (e.g., a counterflow heat exchanger, a plug-flow reactor with a counterflow heating jacket, a counterflow gas absorber and one-dimensional flow of an ideal gas, etc). In this subsection, the characteristic-based MPC technique is developed for these more complex systems having multiple characteristics. For the purpose of illustration and brevity, we restrict our discussion to a 2 × 2 quasilinear system. The general case can be handled analogously. Consider a quasilinear 2 × 2 hyperbolic system modeled by

∂υ1 ∂υ1 + a1 ) f1(υ1,υ2,u) ∂t ∂x ∂υ2 ∂υ2 + a2 ) f2(υ1,υ2,u) ∂t ∂x

dx ) a1 dt dx ) a2 dt

characteristic C1: characteristic C2:

(19)

The variation of state variables υ1 and υ2 can be described by the ordinary differential equations

dυ1 ) f1(υ1,υ2,u) dt dυ2 ) f2(υ1,υ2,u) dt

along characteristic C1 (20) along characteristic C2

Then, the evolution of the state and output variables can be predicted from the characteristic ODEs. In most applications, the PDEs in eq 18 are nonhomogeneous, and f1(υ1,υ2,u) and f2(υ1,υ2,u) are nonlinear functions. For these systems, the characteristic ODEs are coupled with respect to the two characteristic curves, and the future state variables at one spatial point have to be obtained by simultaneously integrating both characteristic ODEs along two nonparallel characteristic curves. Predictions of future output values can be obtained by discretizing the initial state at a finite number of spatial points and then projecting the characteristic curves from each of these points, computing the values of the state variables at the intersection points. Figure 1 illustrates the calculation of the state variables at point P from the values at points Q and R. The segment QR is the domain of dependence of point P, given that the values of state variables at point P are completely defined by the state-variable values on the segment QR. The spatial coordinates of point P are obtained from eq 19 as x(P) 1 x(P) 1 dx + t(Q) ) ∫x(R) dx + t(R) ∫x(Q) a1 a2

(21)

and the time coordinate of P is

t(P) ) t(Q) +

x(P) 1 dx ∫x(Q) a1

(22)

In the case where a1 and a2 are constants, the spatial coordinate and time coordinate of P can be written as

x(P) )

a1x(R) - a2x(Q) + a1a2[t(R) - t(Q)] a1 - a2

(23) a1t(Q) - 2a2t(Q) + a2t(R) + x(R) - x(Q) t(P) ) a1 - a2 Then, the values of the state variables υ1 and υ2 at point P can be obtained from eq 20 using Euler’s method as follows

(18)

y ) h(υ1(xout),υ2(xout)) where y is the process output and h is an output function. If a1 ) a2, the system has one single characteristic. Otherwise, the system has two characteristics. The characteristics for eq 18 are two curves determined by

x(P) dx ∫x(Q) a1

υ1(P) ) υ1(Q) + f1(Q) υ2(P) ) υ2(R) + f2(R)

x(P) dx ∫x(R) a1

(24)

A more accurate estimate can be obtained by replacing f1(Q) by 1/2[f1(Q) + f1(P)] and f2(R) by 1/2[f2(R) + f2(P)] or using other more accurate numerical integration schemes.

2144 Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004

Figure 2. State-variable prediction along characteristics.

Figure 1. Characteristic curves of a system of PDEs with multiple characteristics.

By varying point P and repeating the procedure, the values of the state variables at different grid points and different future times can be calculated. If the state variable υ1 at xout is the output to be controlled, the prediction of the output can be carried out as shown in Figure 2. The outputs at the future sample times can be obtained from the values of the state variables at the intersection points. In the above output prediction procedure, the domain of dependence determining the state-variable values at one point is approximated by two points. Because the discretization affects the accuracy of the prediction, a careful discretization is needed to obtain the desired accuracy. On the other hand, the approximation of the state-variable value within segment QR by those at points Q and R reflects the true solution of hyperbolic PDE systems more closely than other numerical solution methods. This method can accommodate the use of larger spatial grids and time steps with a minimal loss of accuracy. Using the output prediction method described above, the value of the output for a prediction horizon time can be obtained for some specified control actions. The control action can be calculated using proper nonlinear MPC schemes, as discussed for scalar systems in the previous subsection. 4. Closed-Loop Stability of the CBMPC A definition of stability for infinite-dimensional systems was proposed in the literature to address the infinite nature of DPSs.24,33 This paper has been focused on output model predictive control, and therefore, only input-output stability is addressed here. For the proposed CBMPC, closed-loop stability for linear and quasilinear systems is proven differently. The stability of linear CBMPC can be analyzed unconditionally because of its off-line form, whereas the stability of general quasilinar CBMPC has to be analyzed by adding a constraint. 4.1. Linear Systems. As discussed in the preceding section, the unconstrained CBMPC for linear PDE systems has the form of an off-line control law. Hence, it is possible to analyze the closed-loop stability by examining the denominator polynomials of the closed-

loop transfer functions. This method was used to establish stability results for unconstrained MPC based on impulse response models.34 The output prediction for linear systems can be explicitly expressed in terms of the current statevariable profiles and the future output, as in eq 10. Because the current state-variable value at every spatial point can be determined by the boundary-condition value υb and the past control actions, the future output can also be formulated in terms of both past and future control actions as well as the boundary-condition value υb. Using a sampling time ts, define

c hi ) (eibts - e(i-1)bts) , b

i ) 1, ..., m

(25)

and u ) [u0 u1 ‚‚‚ umc-1]T for current and future control actions, u-1 ) [u-1 u-2 ‚‚‚ u-(m-1)]T for the control actions in the past (m - 1) sample instants, vb ) [embtsυb embtsυb ‚ ‚‚ embtsυb]T, and matrices

[

h1 0 ... 0 0 h2 h1 ... 0 0 l S1 ) hm-1 hm-2 ... hm-mc+1 hm hm-1 ... hm-mc+2

m-m hi ∑i)1 m-m +1 ∑i)1 hi

[

h2 h3 ... h3 h4 ... S2 ) ... hm 0 ... 0 0 ...

hm-1 hm hm 0 0 0

0 0

c

]

c

]

(26)

(27)

Then, the output for the next m time instants, ti ) t0 + its, i ) 1, ..., m, can be formulated as

yˆ ) vb + S1u + S2u-1

(28)

where yˆ is defined as in eq 11. Comparing eqs 28 and 13 leads to the relation

S ) S1 y0 ) vb + S2u-1

(29)

In CBMPC for linear PDE systems, a sequence of control actions is obtained explicitly, as in eq 15. The input to be implemented is obtained by selecting the first element of the sequence. Hence, the control law developed for linear PDE systems can be written as

Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004 2145

u0 ) bTu

(30)

) bT(S1TS1)-1S1T(yr - vb - S2u-1)

where bT ) [1 0 ‚‚‚ 0]T. According to the formulation of the control law in eq 30, the following theorem establishes the stability of the proposed CBMPC for linear PDE systems. Theorem. The linear CBMPC in eq 30 yields a stable closed-loop process for control horizon mc chosen to be sufficiently small. Proof. By applying Z transforms to eq 30, the control law in Z transform has the general form

Dc(z) u0(z) ) Ncr(z) yr(z) - Ncb(z)υb(z)

(31)

where Dc(z), Ncr(z), and Ncb(z) indicate polynomial functions of z and the Dc(z) polynomial can be written explicitly as

Dc(z) ) 1 + bT(S1TS1)-1S1TS2[z-1

z-2

‚‚‚

z

-(m-1) T

]

(32)

Stability of the scheme is determined by the roots of Dc(z). By defining the backward shift operator q ) z-1, the equivalent stability condition requires that the characteristic polynomial

C(q) ) 1 + bT(S1TS1)-1S1TS2[q

q2

‚‚‚ q

m-1 T

] ) 0 (33)

have roots outside the unit circle. We establish this result by demonstrating stability for mc ) 1. From eqs 26 and 27, S1TS1 and S1TS2 can be formulated. Substituting the obtained expressions into eq 33 yields m

j

∑ ∑ j)1 i)1 (

m-1 j

hi)2 +

∑ ∑ j)1 i)1

m-k j

hihj+1q + ‚‚‚ + (

hihj+k)qk + ∑ ∑ j)1 i)1

... + h1hmqm-1 ) 0 (34) Starting with the last coefficient (h1hm), one can verify that all of the terms in each coefficient are included in the next. From the definition of hi in eq 25, it is clear that all hi, i ) 1, ..., m, have the same sign. Consequently, the monotonic condition for a polynomial to have roots outside the unit circle is satisfied.35 Therefore, all of the roots of the characteristic polynomial lie outside the unit circle. This proves the stability of the CBMPC for linear PDE systems. 9 The above proof establishes the closed-loop stability for unconstrained linear CBMPC in the case of mc ) 1. As the control horizon increases, the complexity of characteristic polynomials increases significantly and stability analysis using the above approach becomes much more involved or even impossible. Because a control horizon of 1 is often adequate for achieving satisfactory performance in many applications, the above establishment of stability adds a desirable feature to the linear CBMPC formulation. 4.2. Quasilinear Systems. In contrast to the linear case, CBMPC for quasilinear PDE systems does not yield a control law in closed form. Hence, it is not possible to analyze its stability on the basis of the closedloop transfer functions. The stability issues are similar to those encountered in the analysis of finite-dimen-

Figure 3. Output prediction along characteristic curves for scalar PDEs.

sional nonlinear MPC systems. Some modifications such as a terminal constraint or terminal cost are usually added to ensure closed-loop stability.10 In this subsection, the stability of quasilinear systems is considered by adding a terminal constraint and examining the objective functions. Consider the first-order quasilinear hyperbolic system described by eq 8, with constant boundary condition at x ) 0 [i.e., υ(x)0) ) constant] and x ∈ [0, 1]. In the design of the CBMPC, the infinite-dimensional state variable υ is discretized at spatial points x1, x2, ..., xm. The output for the next p sample instants is obtained by integrating characteristic ODEs with initial values υ(xm), ..., υ(x2), υ(x1) and boundary value υb (see Figure 3). The design parameters of the CBMPC are chosen such that p g m + mc - 1, where p is the prediction horizon, mc is the control horizon, and m is the number of the discretization points. The control actions in the first mc sample instants are free variables and the remaining p - mc variables are specified as umc-1 ) umc ) ... ) up-1, or ∆umc ) ∆umc+1 ) ... ) ∆up-1 ) 0. The selection p g m + mc - 1 ensures that the effects of chosen control actions are included in the prediction horizon. On the basis of the current state-variable profile υ, the CBMPC can be designed to minimize a performance objective function. For simplicity of discussion, the objective function here takes the form

J0(υ) ) minJ(υ,u) ∆u

p

∑[yr - y(ti,υ,u)]2 ∆u i)1

) min

(35)

with no penalty on the control increments. Note, however, that the discussion below applies to the case with small cost parameters for control increments as well. Defining

l(ti,υ,u) ) [yr - y(ti,υ,u)]2 eq 35 can be written as

(36)

2146 Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004 p

∑l(ti,υ,u) ∆u i)1

J0(υ) ) min

(37)

One way of ensuring stability for finite-horizon MPC is to add a “terminal constraint” that forces the states or output to a particular value or set at the end of the prediction horizon. In this subsection, the following terminal constraint is used to ensure the stability of the CBMPC

l(tp,υ,u) e min{l(t1,υ,u1), l(t2,υ,u2), ..., l(tp-1,υ,up-1)} (38) This constraint requires the terminal cost to be the smallest among the costs within the prediction horizon. With this constraint, if eq 38 is satisfied by the control sequence (u0, u1, ..., umc-1), then it is also satisfied by the control sequence (u1, u2, ..., umc-1, umc-1) at the next sample instant. Theorem. The CBMPC that satisfies eq 35 subject to eqs 8 and 38 is stabilizing if the process in eq 8 is controllable at y ) yr. Proof. Suppose that the control problem of minimizing the cost function in eq 35 subject to eq 8 and constraint 38 is solved from the current state-variable value υ(x1), υ(x2), ..., υ(xm), yielding the optimal control sequence 0 (υ)} u0(υ) ) {u00(υ), u01(υ), ..., um c-1

(39)

where the superscript 0 indicates the optimal value. Under this optimal control sequence, the optimal objective function can be obtained as

J0(υ) ) J(υ,u0(υ))

(40)

Using the mathematical description in eq 16 of the preceding section, the output trajectory under the optimal control sequence can be written as 0 y0(ti,υ) ) φυ(υxm-i+1,u00,...,ui-1 ), 0 ), y0(ti,υ) ) φυ(υxm-i+1,u00,...,um c-1

i ) 1, ..., mc - 1 i ) mc, ..., m - 1

0 0 y0(ti,υ) ) φυ(υb,um-i ,...,um ), c-1

i ) m, ..., m + mc - 2 0 ), y0(ti,υ) ) φυ(υb,um c-1

i ) m + mc - 1, ..., p (41)

Note that the predicted output at t g tm+mc-1 is a 0 function of only um for constant υb. c-1 Implementing the first control move u00 of the optimal control sequence from t0 through t1, we obtain the output y(t1) ) y0(t1,υ) and a new state-variable profile υ′ at spatial points x′i, i ) 1, ..., m, at t1. Then the problem of minimizing the cost function subject to the constraints at t1 is solved from the new state-variable values at the discrete points υ′(x′1), υ(x′2), ..., υ′(x′m). Define a control sequence u ˜ at t1 based on the optimal control sequence obtained at t0 as

u ˜ ) {u1 , u2 , ..., umc-2 , umc-1 , umc-1 } 0

0

0

0

0

(42)

Then, the output trajectory under u ˜ at t1 can be obtained

y(υ′,u ˜ ) ) {y0(t2,υ), ...., y0(tp,υ), y0(tp+1,υ)}

(43)

From eq 41, it is easy to see that y0(tp+1,υ) ) y0(tp,υ). Then u ˜ satisfies process model eq 8 and constraint 38 at t1. Therefore, the optimal value of the objective function at t1 is not greater than the objective-function value obtained using u ˜ . The following inequality holds

˜) J0(υ′) e J(υ′,u

(44)

Because p

J(υ′,u ˜) )

[yr - y0(ti,υ)]2 + [yr - y0(tp,υ)]2 ∑ i)2

) J0(υ) - l0(t1,υ) + l0(tp,υ)

(45)

and l0(tp,υ) e l0(t1,υ)) (because of the terminal constraint at t0), the optimal objective-function values at two consecutive sample instants satisfy

J0(υ′) e J0(υ)

(46)

Therefore, the positive cost function J0 is nonincreasing, which implies that the cost function J0 converges to zero or a positive value. If J0 converges to a steady-state positive value, then J0(υ′) ) J0(υ) and l0(t1,υ) ) l0(t2,υ) ) ... ) l0(tp,υ) * 0. Therefore, the output converges to a steady-state value y0s * yr under the optimal control sequence u0s . Because the process is controllable at y ) yr, there exists a control sequence ud that drives the output y to the set point yr at the steady state and satisfies the constraint. Therefore, J(υ,ud) < J0(υ,us), which contradicts the fact that u0s is the optimal control action. Then, eq 47 becomes

J0(υ′) < J0(υ) when y0(υ) * yr J0(υ′) ) J0(υ) ) 0 when y0(υ) ) yr

(47)

Thus, J0 is a Lyapunov function and the output converges to the set point yr. This proves that the developed CBMPC is stabilizing. 9 As in other nonlinear MPC schemes, the CBMPC for general quasilinear PDE systems does not guarantee stability without a limiting constraint. In this section, a terminal constraint is added to ensure the stability of the proposed CBMPC. The limiting constraint might affect the performance and computational efficiency of the CBMPC. As a special case of CBMPC, linear CBMPC exhibits the guaranteed stability so that the terminal constraint is not necessary. 5. Simulation Results In this section, the proposed CBMPC algorithm for hyperbolic PDE systems is evaluated via simulation, using two plug-flow reactor (PFR) systems: one with uniform heating and one with counterflow heating. 5.1. Plug-Flow Reactor with Uniform Heating. Consider a nonisothermal plug-flow reactor in which two first-order reactions take place in series

AfBfC where A is the reactant species, B is the desired product, and C is an undesired product. The reactions are endothermic. The reactor is heated with a steam jacket.

Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004 2147 Table 1. Model Parameters for a PFR parameter

value

units

parameter

value

units

υl L Vr E1 E2 k10 k20 Hr1 cpmj Fmj

1 1.0 10.0 20 000.0 50 000.0 5.0 × 1012 5.0 × 102 0.5480 0.8 0.10

m/min m L kcal/kmol kcal/kmol min-1 min-1 kcal/kmol kcal/(kg K) kg/L

R Fm Uw cpm CA0 CB0 Tr0 Hr2 Vj Tj0

1.987 0.09 0.20 0.231 4 0 320 0.9860 8 375

kcal/(min K) kg/L kcal/(min K) kcal/(kg K) mol/L mol/L K kcal/kmol L K

Under the assumptions of no radial concentration gradients in the reactor, constant volume of the liquid within the reactor, constant density and heat capacity of the reactants, and negligible diffusion and dispersion, the following model for the system can be obtained using material and energy balances24 Figure 4. Simulation results of the PFR using CBMPC for a setpoint change from 0.8 to 1 mol/L.

∂CA ∂CA ) -υl - k10e-E1/RTrCA ∂t ∂x ∂CB ∂CB ) -υl + k10e-E1/RTrCA - k20e-E2/RTrCB ∂t ∂x (48) ∂Tr (-∆Hr1) ∂Tr -E1/RTr k e CA + ) -υl + ∂t ∂x Fmcpm 10 (-∆Hr2) Uw k20e-E2/RTrCB + (T - Tr) Fmcpm FmcpmVr j subject to the boundary conditions

CA(0,t) ) CA0,

CB(0,t) ) 0,

Tr(0,t) ) Tr0 (49)

where CA and CB are the concentrations of the species A and B, respectively, in the reactor; υl is the velocity of the reactants through the reactor; Tr is the temperature of the reactor; ∆Hr1 and ∆Hr2 are the enthalpies of the two reactions; Fm and cpm is the density and heat capacity, respectively, of the fluid in the reactor; Vr is the volume of the reactor; Uw is the heat-transfer coefficient; k10, k20, E1, and E2 are the Arrhenius constants and activation energies, respectively, of the reactions; CA0 and Tr0 are the concentration and temperature, respectively, of the inlet stream in the reactor; Tj is the spatially uniform temperature in the jacket, which is manipulated to control the concentrations; and x is the spatial coordinate along the reactor with x ∈ [0, 1]. The outlet concentration of product B is the variable to be controlled. The values of the constants in eq 48 used for the simulations are listed in Table 1. For simulation purposes, the process is represented by discretization of the PDE model along the spatial direction into 500 points. The control horizon mc is set to 1. The prediction horizon is taken to be equal to the residence time. The control action is calculated to minimize the objective function m

J)

[CrB - C ˆ B(ti,x)1)]2 ∑ i)1

(50)

where CrB is the set point of the outlet concentration CB, and C ˆ B(ti,x)1) is the predicted output in the next m time instants. At the initial state, the system is at a steady state with the jacket temperature Tj ) 359.5 K. Figure 4

Table 2. Computations of CBMPC vs Conventional MPC in PFR MPC using characteristics

lumped MPC

4.73 × 104 0.018 5 points

4.7 × 107 0.022 50 points

flops (Matlab) ISE state estimation

illustrates the simulation of the response of CB to an output set-point change from CrB ) 0.8 mol/L to CrB ) 1 mol/L. The simulation results indicate that CBMPC provides good tracking performance, as the distributed state variable displays a smooth and rapid response, reaching the desired set point within slightly more than 1 min (or a single residence time). The conventional MPC approach requires the discretization of the PDE model, which can result in a large number of ODEs in the approximation. To investigate whether CBMPC is a viable alternative from a computational point of view, simulations were performed using each approach. The integral of the square of the error, ISE, given by

ISE )

∫0t

ISE

[CrB - CB(τ,x)1)]2 dτ

(51)

was used to evaluate the performance of the controller. In each case, the control action was computed by minimizing the objective function in eq 50. The computation of the control action based on the finite-difference method was conducted by discretizing the PDE model into 50 spatial nodes. The ISE values and computational demands for 2.2 min of reactor simulation are listed in Table 2. Note that for equivalent control accuracy, as measured by the ISE, the CBMPC approach required 3 orders of magnitude less computation and a considerably coarser discretization. Of course, the approach taken in lumped MPC was simplistic, and more sophisticated numerical approaches might provide improved performance; however, these caveats could also be applied to CBMPC. 5.2. Counterflow PFR. In many applications, the temperature of the heating or cooling medium within a PFR in the jacket of a PFR is not spatially uniform. In such situations, the distributed state variables include the temperature within the PFR jacket. When the jacket and reactor velocities are not identical, such PFR systems have to be described along two different char-

2148 Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004

acteristics. In this subsection, the application of CBMPC for a nonisothermal PFR with a counterflow heating medium in the jacket is considered. The reactor described in the preceding subsection is used except that, in this case, the jacket fluid has a spatially varying temperature and the fluid velocity of the jacket is manipulated to control the outlet concentration of product B. The process model is given by

∂CA ∂CA ) -υl - k10e-E1/RTrCA ∂t ∂x ∂CB ∂CB ) -υl + k10e-E1/RTrCA - k20e-E2/RTrCB ∂t ∂x (52) ∂Tr (-∆Hr1) ∂Tr -E1/RTr ) -υl + k e CA + ∂t ∂x Fmcpm 10 (-∆Hr2) Uw k20e-E2/RTrCB + (T - Tr) Fmcpm FmcpmVr j ∂Tj Uwj ∂Tj )u + (T - Tj) ∂t ∂x FmjcpmjVj r subject to the boundary conditions

CA(0,t) ) CA0,

CB(0,t) ) 0,

Tr(0,t) ) Tr0, Tj(1,t) ) Tj0 (53)

where Fmj and cpmj are the density and heat capacity, respectively, of the fluid in the jacket; Vj is the volume of fluid in the jacket; Uwj is the heat-transfer coefficient from the jacket to the reactor; and u is the velocity of fluid in the jacket. In eq 52, the first three equations can be described by a system of ODEs along the characteristic vector field

ξ1 )

∂ ∂ + υl ∂t ∂x

(54)

Because the heating fluid in the jacket flows in a different direction and at a different rate than the reactants in the reactor, the PDE representing the jacket temperature has to be described by a different vector field

ξ2 )

∂ ∂ -u ∂t ∂x

(55)

Development of the CBMPC for this example follows the approach described in subsection 3.2. The controller was designed to minimize the objective function m

J)

[CrB - C ˆ B(ti,x)1)]2 + λ(∆u)2 ∑ i)1

(56)

The model parameters listed in Table 1 were used for the simulation. The simulation was performed by discretizing the PDE model using 60 spatial nodes and solving the resulting system of ODEs. The number of discretization points for the control calculation was taken to be m ) 10. The parameters of the nonlinear quadratic DMC algorithm used were p ) 10, mc ) 1, and λ ) 0.0004 as tuning parameters. The initial condition was assumed to be at steady state for a flow rate of 0.5 m/min.

Figure 5. Performance comparison of CBMPC (s) vs finitedifference-based MPC (- - -, - - -) in the counterflow PFR. Table 3. Computations of CBMPC vs Conventional MPC in Counterflow PFR control method

flops (Matlab)

CBMPC (m ) 10) CBMPC (m ) 20) finite-difference-based MPC (m ) 10) finite-difference-based MPC (m ) 20)

3.5826 × 104 1.4003 × 105 6.31 × 106 1.40 × 107

Figure 5 shows the output response of the CBMPC for a set-point change in CB(x)1) from 0.83 to 1 mol/L (solid line). The proposed CBMPC yields a smooth output set-point tracking response. The process output converges to the set point quickly without large overshoot and oscillations. This performance of the CBMPC was also compared with that of an MPC based on a finite-difference method. For the proposed CBMPC, satisfactory performance can be obtained using discretization m ) 10 (almost the same response can be obtained for m ) 20). On the other hand, finitedifference-based MPC using 10 discretization points (m ) 10) does not provide adequate performance, and finer discretization is required in finite-difference-based MPC in order to obtain satisfactory performance. Using 20 discretization points (m ) 20), the finite-differencebased MPC yields a closed-loop performance comparable to that of CBMPC with m ) 10. A slightly larger offset was observed with the finite-difference-based MPC. The computational efficiency of the CBMPC and finitedifference-based MPC are compared in Table 3, and several points are worthy of note. First, finer discretization is required to achieve controller performance for the finite-difference-based MPC similar to that obtained with the proposed CBMPC approach. Second, the online computation required in the CBMPC approach is considerably less than that required in the finitedifference-based MPC, typically by several orders of magnitude. Although there might be a number of approaches that can be used to improve the performance of any discretization-based MPC techniques, they could be used with equal effect for CBMPC. Thus, the proposed CBMPC approach has a number of advantages for industrial applications. 6. Conclusions It is well recognized that MPC can provide effective control that outperforms the standard feedback control approaches. Despite the widespread success of MPC, there is a paucity of results for distributed parameter systems. The CBMPC formulation developed here makes good use of the geometric properties of PDE models and

Ind. Eng. Chem. Res., Vol. 43, No. 9, 2004 2149

provides a control approach that exploits the dynamics of hyperbolic DPSs. The CBMPC approach overcomes some of the barriers encountered in applying MPC to DPSs and can produce high-accuracy control with efficient on-line computation. The CBMPC approach uses the method of characteristics to predict the future output. The high prediction accuracy of this method is the basis for the high accuracy of the proposed CBMPC. For linear PDE systems, the proposed control design produces off-line control laws, and closed-loop stability is verified when the control horizon is tuned to 1. Simulation results indicate that the CBMPC generates satisfactory performance to set-point changes. It greatly reduces the online computation demand in comparison to that of finitedifference-based MPC. The CBMPC discussed in this paper has been limited to unconstrained MPC with no model-plant mismatch. For hyperbolic DPSs with constraints, the developed CBMPC can be implemented by incorporating the existing constraint-handling techniques used in the MPC of lumped parameter systems. Acknowledgment The authors acknowledge financial support provided by the Mechanical Wood Pulps National Centre of Excellence and Natural Sciences and Engineering Research Council (NSERC) of Canada. Literature Cited (1) Culter, C. R.; Ramaker, B. L. Dynamic Matrix Control: A Computer Control Algorithm. Presented at the 86th National AIChE Meeting, Houston, TX, Apr 1-5, 1979. (2) Richalet, J.; Rault, A.; Testud, L.; Papon, J. Algorithmic Control of Industrial Process. in Proceedings of the Fourth IFAC Symposium on Identification and System Estimation; Elsevier Science: New York, 1976; p 1119. (3) Badgewell, T. A. Robust Stability Conditions for SISO Model Predictive Control Algorithms. Automatica 1997, 33 (7), 1357. (4) Qin, J. S.; Badgwell A Survey of Industrial Model Predictive Technology. Control Eng. Pract. 2003, 11 (7), 733. (5) Rani, K. Y.; Unbehauen, H. Study of Predictive Controller Tuning Methodology. Automatica 1997, 33 (12), 2243. (6) Rao, C. V.; Rawlings, J. B. Linear Programming and Model Predictive Control. J. Process Control 2000, 10, 283. (7) Rawlings, J. B. Tutorial Overview of Model Predictive Control. IEEE Control Syst. Mag. 2000, 20 (3), 38. (8) Chen, H.; Allgo¨wer, F. A Quasi-Infinite Horizon Nonlinear Model Predictive Control Scheme with Guaranteed Stability. Automatica 1998, 34 (10), 1205. (9) Henson, M. A. Nonlinear Model Predictive Control: Current Status and Future Directions. Comput. Chem. Eng. 1998, 23, 187. (10) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. Constrained Model Predictive Control: Stability and Optimality. Automatica 2000, 36, 789. (11) Scokaert, P. O. M.; Mayne, D. Q.; Rawlings, J. B. Suboptimal Model Predictive Control (Feasibility Implies Stability). IEEE Trans. Autom. Control 1999, 44 (3), 648. (12) Sistu, P. B.; Bequette, B. W. Nonlinear Model-predictive Control: Closed-loop Stability Analysis. AIChE J. 1996, 42 (12), 3388.

(13) Eaton, J. W.; Rawlings, J. B. Feedback Control of Chemical Processes Using Online Optimization Techniques. Comput. Chem. Eng. 1990, 14 (4-5), 469. (14) Patwardhan, A. A.; Wright, G. T.; Edgar, T. F. Nonlinear Model-predictive Control of Distributed-parameter Systems. Chem. Eng. Sci. 1992, 47 (4), 721. (15) Toure´, Y.; Biston, J.; Gills, G. Modelling of a Distributed Parameter Process with a Variable Boundary: Application to its Control. Chem. Eng. Sci. 1994, 49 (1), 61. (16) Irizarray-Rivera, R.; Seideer, W. D. Model-Predictive Control of the Czochralski Crystallization Process Part I. Conduction-Dominated Melt. J. Cryst. Growth 1997, 178 (4), 593. (17) Dufour, P.; Couenne, F.; Toure´, Y. Model Predictive Control of a Catalyptic Reverse Flow Reactor. IEEE Trans. Control Syst. Technol. 2003, 11 (5), 705. (18) Dufour, P.; Toure´, Y.; Blanc, D. On Nonlinear Distributed Parameter Model Predictive Control Strategy: Online Calculation Time Reduction and Application to an Experimental Drying Process. Comput. Chem. Eng. 2003, 27 (11), 1533. (19) Campbell, J. C.; Rawlings, J. B. Predictive Control of Sheet- and Film-forming Processes. AIChE J. 1998, 44 (8), 1713. (20) Vanwerp, J. G.; Braatz, R. D. Fast Model Predictive Control of Sheet and Film Processes. IEEE Trans. Control Syst. Technol. 2000, 8 (3), 408. (21) Bendersky, E.; Christofides, P. D.; Optimization of Transport-Reaction Processes Using Nonlinear Model Reduction. Chem. Eng. Sci. 2000, 55, 4349. (22) Armaou, A.; Christofides P. D.; Dynamic Optimization of Dissipative PDE Systems Using Nonlinear Order Reduction. Chem. Eng. Sci. 2002, 57, 5083. (23) Hanczyc, E. M.; Palazoglu, A. Sliding Mode Control of Nonlinear Distributed Parameter Chemical Processes. Ind. Eng. Chem. Res. 1995, 34, 557. (24) Christofides, P. D.; Daoutidis, P. Feedback Control of Hyperbolic PDE Systems. AIChE J. 1996, 42 (11), 3063. (25) Christofides, P. D.; Daoutidis, P. Robust Control of Hyperbolic PDE Systems. Chem. Eng. Sci. 1998, 53, 85. (26) Bryant, R. L.; Chern, S. S.; Gardner, R. B.; Goldschmidt, H. L.; Griffiths, P. A. Exterior Differential Systems; SpringerVerlag: New York, 1991. (27) Edelen, D. G. B. Applied Exterior Calculus; John Wiley & Sons Inc.: New York, 1985. (28) Arnold, V. I. Geometric Methods in the Theory of Ordinary Differential Equations; Springer-Verlag: New York, 1988. (29) Duchateau, P.; Zachmann, D. Applied Partial Differential Equations; Harper & Row Publisher: New York, 1989. (30) McOwen, R. Partial Differential Equations; Prentice Hall Inc.: Upper Saddle River, NJ, 1996. (31) Simpson, D. A Numerical Method of Characteristics for Solving Hyperbolic Partial Differential Equations. PhD Dissertation, Iowa State University, Ames, IA, 1967. (32) Brosilow C.; Joseph, B. Techniques of Model-Based Control; Prentice Hall: Upper Saddle River, NJ, 2002. (33) Curtain, R. F.; Zwart, H. J. An Introduction to InfiniteDimensional Linear Systems Theory; Springer-Verlag: New York, 1995. (34) Garcia, C. E.; Morari, M. Internal Model Control 1. A Unifying Review and Some New Results. Ind. Eng. Chem. Process Des. Dev. 1982, 21 (2), 308. (35) Jury, E. I. Theory and Application of the Z-Transform Method; Huntington: New York, 1964.

Received for review August 8, 2003 Revised manuscript received December 4, 2003 Accepted February 25, 2004 IE030653Z