2064
Ind. Eng. Chem. Res. 2007, 46, 2064-2076
Measurement-Based Optimization and Predictive Control for an Exothermic Tubular Reactor System Wei Wu* and Chao-Wei Chen Department of Chemical Engineering, National Yunlin UniVersity of Science and Technology, Douliou, Yunlin 64002, Taiwan, Republic of China
In this article, the analytic optimization algorithm connected to the measurement-based predictive control framework is implemented on an exothermic tubular reactor system. The two stages of design procedure include (i) the desired input/output references, which are determined by the steady-state optimization approach, and (ii) the output regulation design of nonlinear distributed parameter systems, which is addressed using nondistributed predictive controls. Under the assumption of steady-state and plug-flow characteristics, the bang-bang type of extremal control is applied. Because of the fact that inlet perturbations could trigger the hot spots or thermal runaway, we propose two dependent manipulated variables to dominate the heat exchange function for cooling devices. With respect to specified state/input constraints, the output feedback architecture, which uses a few sensors for measurement of the reactor temperature at the prescribed axial position, is successfully demonstrated. It is a specific measurement-based feedback technique used to combat the effects of heat or unknown disturbances. All tests show that the no-offset output tracking is achieved and the undesired peak temperature is removed while physical constraints and unknown disturbances are being considered simultaneously. 1. Introduction The use of partial differential equation (PDE) models for system analysis and the design of estimation and control algorithms has become increasingly important recently.1-3 However, those control strategies for a class of nonlinear distributed parameter systems were too complex to implement in practice, because of the original distributed control framework. Although tubular reactors have been largely used in chemical process industry, they cannot be operated at excessive temperatures, because of exothermic reactions and inlet perturbations, such that a undesired peak (hot-spot) temperature or thermal runaway4 would induce undesired side reactions and thermal decomposition of the products. Therefore, the problem of suppressing the magnitude of the hot-spot temperature of tubular reactors using a coolant device has been recently addressed.5,6 Because the position and magnitude of the hotspot temperature cannot be directly established via explicit nonlinear functions, the state feedback control framework or steady-state optimization approach is often applied. Recently, model predictive control (MPC) algorithms have been widely studied and applied in many chemical processes.7 Unfortunately, reports of MPC based on partial differential equation (PDE) models are relatively scarce. Shang et al.8 proposed a modified MPC scheme that is capable of effectively driving a distributed process output to its setpoint, Dufour et al.9 provided a general MPC framework for a PDE model of the autoclave curing process, and Dufour and Toure10 developed a multiple-input/multiple-output (MIMO) MPC for PDE systems. Notably, input/output constraints are handled in the optimization task using a nonlinear programming method, and the objective of approximating a PDE model is to reduce the calculation time for the online MPC procedure. In addition to input/output constraints and the characteristic of distributed parameter systems, the exothermic tubular reactor system should consider the other crucial and practical problems on hot spots, coolant devices, and process optimization. Moreover, low-order * To whom correspondence should be addressed. Tel.: 886-55342601. Fax: 886-5-5312071. E-mail:
[email protected].
predictive control of a class of PDE systems has been developed.11,12 Note that the decomposition and transformation is used to dominate the model reduction such that closed-loop stability in the face of state/input constraint should be dependent on the state estimation design and the number of output measurements. Regarding optimal control and design, the terminal-cost optimization problem was usually implemented for batch process optimization for the past several years.13-16 Notably, the path and terminal constraints are considered simultaneously, and the explicit synthesis of optimal control input was of the bangbang type. Moreover, the optimal control design was applied for tubular reactor systems. Smet et al.17 described an optimal temperature control problem of a steady-state plug-flow reactor (PFR), which has a tradeoff between process performance and heat loss. Birk et al.18 proposed a computation of optimal feed rates that can maximize the overall profit of the process. It is noted that those feedforward-like control schemes could be seriously degraded, because of model errors or disturbances. Recently, the boundary feedback control for nonlinear PDE models was presented to improve closed-loop robustness.19,20 For the some control problems on the elimination of hot spots, boundary feedback design, and physical constraints, in this article, the analytical optimization algorithm connected to the measurement-based predictive control framework is implemented by two stages of design procedure. The optimal temperature and concentration profiles of exothermic tubular reactors were determined under the assumption of steady-state and plug-flow characteristics. For the first stage of steady-state optimization approach, Pontryagin’s minimum principle21 was applied to the terminal-cost optimization problem, in which the extended terminal cost criterion may enable a tradeoff between process performance and the overall heat effect. The bangbang type of extremal control can facilitate the optimization algorithm for tubular systems. In the second stage, the on-line computation and feedback-based implementation scheme was studied according to the optimum operating policy and approximate difference-based models. Moreover, nondistributed model predictive control (NMPC) strategies were developed,
10.1021/ie0611296 CCC: $37.00 © 2007 American Chemical Society Published on Web 03/06/2007
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2065
and two dependent manipulated variables were used to dominate the heat exchange function for cooling devices. Particularly, the boundary feedback control, which uses a few sensors for measurement of the reactor temperature at the prescribed axial position, could accurately dominate the heat exchange of cooling devices such that the undesired peak temperature disappeared throughout the entire reactor. Compared to recent works for tubular reactor temperature control, our scheme is different from the design of multiple cooling zones22 or the cooling coil apparatus.23 Through theoretical analysis and closed-loop simulations, both NMPC techniques, with respect to state/input constraints, show that the no-offset output regulation of the illustrative phthalic anhydride reactor system can be achieved by tuning a single parameter. 2. Exothermic Tubular Reactor System 2.1. Process Description and Modeling. Several phthalic anhydride reactors have been presented by Froment24 and Chen and Sun.25 This reactor is packed with a V2O5 catalyst, the feed stream is o-xylene and air, and the outside of the tubular reactor is covered with a co-current cooling jacket, as depicted in Figure 1. Assume that the axial diffusion of reactants, the wall capacity, and the heat loss of the reactor are negligible. The reaction kinetics are denoted as pseudo-first-order, because of the large excess of oxygen. Under the mass balance and energy balance, four coupled nonlinear PDEs can be expressed:
∂CA(z,t) ∂CA(z,t) ) -V + µ(1 - )RA(CA,T) ∂t ∂z
(1a)
∂CB(z,t) ∂CB(z,t) ) -V + µ(1 - )RB(CA,CB,T) (1b) ∂t ∂z
δT(z,t) µ(1 - ) δT(z,t) ) -V + R ˜ (CA,CB,T) + Le δt δz Ffcpf 2UL (T (z,t) - T(z,t)) (1c) r0Ffcpf c ∂Tc(z,t) ∂Tc(z,t) 2UL (T(z,t) - Tc(z,t)) ) Vc + ∂t ∂z r′Fccpc
L)4m ) 0.35 k01 ) 2.418 × 109 s-1 k02 ) 2.706 × 109 s-1 k03 ) 1.013 × 109 s-1 E1 ) 1.129 × 105 J/kmol E2 ) 1.313 × 105 J/kmol E3 ) 1.196 × 105 J/kmol -∆H1 ) 1.285 × 106 J/kmol -∆H2 ) 3.276 × 106 J/kmol -∆H3 ) 4.561 × 106 J/kmol CAf ) 0.181 kmol/m3 CBf ) 0 kmol/m3
Ff ) 0.582 kg/m3 Fc ) 1851.456 kg/m3 Fs ) 2000 kg/m3 cpf ) 1045 J kg-1 K-1 cpc ) 483.559 J kg-1 K-1 cps ) 836.0 J kg-1 K-1 r0 ) 0.0125 m r1 ) 0.0225 m U ) 96.02 J kg-1 K-1 s-1 µ)1 V ) 2.06 m/s Tf ) 628 K
where CA and CB are the concentrations of o-xylene and phthalic anhydride, respectively, Le ) [(1 - )Fscps + Ffcpf]/(Ffcpf), and r′ ) r0[(r1/r0)2 - 1]. The reaction rates are given by
( ) ( ) ( )
R1 ) k01 exp -
E1 C RgT A
R2 ) k02 exp -
E2 C R gT B
R3 ) k03 exp -
E3 C RgT A
and
RA ) -R1 - R3
V
Table 1. Process Parameters and Operating Conditions
(1d)
with the boundary conditions
CA(0,t) ) CAf
(2a)
CB(0,t) ) CBf
(2b)
T(0,t) ) Tf
(2c)
Tc(0,t) ) Tcf
(2d)
and
and the initial conditions
R B ) R1 - R2 3
R ˜}
(-∆Hi)Ri ∑ i)1
where CA,ss, CB,ss, Tss, and Tc,ss represent steady-state variables. Most of the symbols have been specified in the Notation section, and the process parameters and operating conditions are given in Table 1. In regard to the process design and control optimization, the operating policies of this exothermic tubular system are stated as follows: (1) The coolant flow rate (Vc) and coolant temperature (Tcf) in the feed are non-negative and bounded. To increase the effectiveness of the heat transfer, both Vc and Tcf are treated as dependent manipulated variables. (2) The intitial steady-state variables CA,ss, CB,ss, Tss, and Tc,ss are treated as the optimum operating conditions. (3) Physical constraints involve (i) bounded inputs and (ii) the reactor temperature profile: bounded inputs are shown by the equations
Vc,min e Vc e Vc,max
(4a)
Tcf,min e Tcf e Tcf,max
(4b)
whereas the reactor temperature profile obeys
CA(z,0) ) CA,ss(z)
(3a)
CB(z,0) ) CB,ss(z)
(3b)
T(z,0) ) Tss(z)
(3c)
Tc(z,0) ) Tc,ss(z)
(3d)
and
Tr(z,t) e Tr(L,t)
∀z ∈ [0,L]
(5)
(4) The feed compositions (CAf) and the feed temperature of the cooling (Tcf) step could be affected by unknown disturbances. (5) The exit product concentration (CB,out) is not measurable; however, the exit reactor temperature (Tout), as the controlled variable, is actually measurable.
2066
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007
Figure 1. Illustration of the phthalic anhydride reactor system covered with a co-current cooling jacket.
Remark 1: In regard to the operating policies of the illustrated example, Vc and Tcf are considered to be variables and are treated as controls, the optimal operating points (initials) are not yet determined, the specified reactor temperature profile is restricted to follow eq 5, and so on. All specifications will be investigated in the article. 2.2. Steady-State System. Inspired by some articles for the same illustrative example,6,26 these results show that the aforementioned exothermic reactions may cause the hot-spot temperature (Tr,max) in the reactor, because of the heat effect. If the reactor has no coolant device, the wall temperature of the reactor (Tw) is used to dominate the heat-exchange effect. As for observing the characteristics of tubular systems, the steadystate model is easily investigated by
dCA,ss µ ) (1 - )RA(CA,ss,Tss) dz V dCB,ss dz
() µ ) ( )(1 - )R (C V B
A,ss,CB,ss,Tss)
(6a) (6b)
dTss µ(1 - ) 2UL R ˜ (CA,ss,CB,ss,Tss) + (T - Tss) (6c) ) dz VFfcpf Vr0Ffcpf w subject to
0ezeL
(7a)
CA,ss g 0
(7b)
CB,ss g 0
(7c)
Tss > 0
(7d)
and
Discussion: Inspired by the temperature policy for a unconvered steady-state PFR,6,27,28 assume that the reactor has perfect coolant function, which means that the overall wall temperature is uniform, because of the fast dynamics of the coolant device. Figure 2 shows that the higher wall temperature (Tw > 645 K) will cause the undesired hot spot or thermal runaway phenomena along the trail of plug flow in the reactor. However, the low wall temperature can reduce the appearance of the hot-spot temperature, but it also reduces the conversion of reactants at the outlet of the reactor. 3. Optimal Temperature Control of the Steady-State System The higher or lower wall temperature could be inappropriate for the steady-state operation of the tubular system. Generally, the stable reactor temperature associated with the high conversion rate of reactions at the outlet of reactor must be satisfied simultaneously under the feasible optimization strategy. 3.1. Optimal Design for Constrained Steady-State System. Referring to Figure 2, the physical constraints are bounded by Tr,max e 680 K and Tcf e 633 K. The off-line computing algorithm is added to obtain the distribution of wall temperature Tw(z), whereas CB is the amount of phthalic anhydride being maximized at the outlet of the reactor and Tr,max is the peak temperature being restricted under the upper limit of the reactor temperature (T < Tr,max). The usual independent variable (the time t) for optimal control problem is replaced by the spatial variable z. Moreover, the terminal-cost optimization of a steadystate tubular system, subject to physical constraints, is shown by
∫0L g(x(z)) dz
min J ) h(x(L)) + K1 uss
and initial conditions
) -x2|z)L + K1 CA,ss(0) ) CAf
(8a)
CB,ss(0) ) CBf
(8b)
and
Tss(0) ) Tf
(8c)
∫0L [uss - x3(z)]2 dz
(9)
subject to
dx ) f(x) + buss dz
(10a)
x(0) ) x0
(10b)
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2067
b)
2UL Vr0Ffcpf
Remark 2: Generally, the optimal control problem is to find an admissible distribution u* so that the trajectory of steadystate model (eq 10) can follow an admissible trajectory under the specified constraint (eq 11). Through the aforementioned minimization algorithm, the objective function (eq 9) involves both objectives of the high conversion rate of reaction and the bounded integral cost due to the heat effect. In addition, the overall temperature profile is influenced by the selection of parameter K1. 3.2. Extremal Control Profile Determination. According to Pontryagin’s minimum principle, the preceding minimization problem is equivalent to minimizing the following Hamiltonian H:
H(x,u*,λ) ) λTf(x) + (λTb)u* } φ + ψu*
(12)
where the adjoint vector λ is the solution of
δH δx
(13a)
∂h | ∂x z)L
(13b)
λ˙ T ) λ(L) )
Obviously, the boundary conditions in eqs 8 and 11 are specified at z ) 0 or z ) L. 3.3. Extremal Control Approach. According to the objective in eq 6, the bang-bang type of extremal control u*, with respect to the minimization of H, is characterized as
u*(z) ) Figure 2. Open-loop steady-state profiles of a tubular system with a uniform wall temperature: (a) concentration of A response, (b) concentration of B response, and (c) reactor temperature response.
with x g 0 and uss e 633 K. The specified physical constraint is set by
x3 e x3|z)L
umax using(z)
if ψ < 0 if ψ ) 0, for z ∈ [z*,L]
}
∂f(x) dψ ) -λT b dz ∂x
( )
2 d 2ψ T∂f(x) ∂f(x) T∂ f(x) ) λ λ b(f(x) + bu*) b ∂x ∂x dz2 ∂x2
where
x ) [CA,ss,CB,ss,Tss]T
[
uss ) Tw
µ (1 - )RA(x) V µ f(x) ) V (1 - )RB(x) µ(1 - ) 2UL R ˜ (x) x VFfcpf Vr0Ffcpf 3
() ()
(15)
If the control u does not appear in eq 15, then the second spatial derivative is written as
Tr,max ) x3|z)L
K1 g 0
(14)
where z* represents the position away from the inlet of the reactor. Because ψ remains identically zero over a finite distance interval [z*,L], the singular control using is obtained by repeatedly differentiating the function ψ until u appears explicitly. In the first step, the spatial derivative of ψ ) λTb is shown as
(11)
i.e.,
and
{
(16)
Moreover, the singular control law is obtained by setting d2ψ/dz2 ) 0:
]
λT using )
( )
∂f(x) ∂f(x) ∂2f(x) b - λT 2 bf(x) ∂x ∂x ∂x 2
f(x) 2 λ b ∂x2 T∂
(17)
Remark 3: Equations 10 and 13 show two point-boundaryvalue problems (TPBVPs). However, the extremal control that is connected to singular control (eqs 14 and 17) performs the analytical optimal control approach. Two parameters (K1, z*)
2068
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007
Figure 3. Optimal state and input profiles for the steady-state tubular system by adjusting K1 ) 1 × 10-5 (bold solid line), K1 ) 5 × 10-5 (solid line), and K1 ) 1 × 10-4 (dashed line): (a) concentration of A distribution, (b) concentration of B distribution, (c) reactor temperature distribution, and (d) wall temperature distribution.
Figure 4. Open-loop dynamic response of a tubular system with a co-current cooling jacket: (a) the steady-state coolant temperature profile is determined by the minimization algorithm, (b) concentration of B response, (c) reactor temperature response, and (d) coolant temperature response.
can be tuned according to the specified optimization algorithm for the high conversion rate of reactions, as well as the safe reactor temperature that is being achieved simultaneously. The
explicit form of using for the steady-state model was obtained from the symbolic software of Maple, which is a generalpurpose computer algebra system.
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2069
Furthermore, the optimal temperature distribution of the wall is established using the following steps: Step 1: For the first interval [0,z*], the wall temperature Tw ) umax ) 633 K is fixed at the upper limit in order to trigger the large conversion rate of reactions. Step 2: For the second interval [z*,L], the singular control Tw ) using(z) is applied to satisfy the dual objectives, i.e., the largest amount of CB at the outlet and no peak temperature in this interval. Step 3: Both tuning parameters (K1, z*) are detuned to search the wall temperature profiles appropriately while the reactor temperature distribution (eq 11) is being satisfied. Remark 4: The extremal control with two parameters (K1, z*) aims to find out the appropriately optimal trajectory of this steady-state system. Under the extremal control with fixed tuning parameters, Figure 3 shows the steady-state optimal profiles for process variables corresponding to computed value for wall temperature. When the uniform wall temperature (K1 ) 1 × 10-5) shown in Figure 3d is considered, the bold solid line of Figure 3c shows that the peak temperature is observed in the reactor. Note that the optimal control is failed and saturated, because of an inappropriate parameter K1. If a K1 value of 5 × 10-5 is tuned, the solid line in Figure 3c shows that the peak temperature is suppressed by constrained optimal design (eqs 9 and 11). In addition, we find that K1 g 5 × 10-5 induces the point z* f 0. When K1 ) 1 × 10-4 and z* ) 0 are given, we find that the smooth distribution of the wall temperature (the dashed line) shown in Figure 3c is acceptable in a real operating environment. Afterward, both parameters (K1, z*) are fixed in our control system, with respect to physical constraints. Moreover, the optimal state distribution as initial conditions are shown as
C*A,ss(z) ) CA(z,0)
(18a)
C*B,ss(z) ) CB(z,0)
(18b)
and
T*ss(z) ) T(z,0)
(18c)
with CAf ) CA,ss(0), CBf ) CB,ss(0), and Tf ) Tss(0). 4. Predictive Control Strategies Referring to the aforementioned optimization strategy for an uncovered tubular reactor, the desired wall temperature distribution is obtained via the off-line computation. According to the illustrative example, the exothermic tubular system is covered with a co-current cooling jacket. Two manipulated inputs (Vc, Tcf) have strong interactive effect on the dynamics of the coolant device. Thus, the steady-state coolant temperature is expected to follow the profile of optimal wall temperature by solving the following minimization algorithm:
(19)
dTc,ss 2UL ) (Tc,ss - Tss) dz V* cr′Fccpc
(20a)
Tc,ss(0) ) T*cf
(20b)
c
T* c,ss(0) ) T* cf
(21a)
T*c,ss(z) ) Tc(z,0)
(21b)
and
Remark 5: Under the analytic optimal algorithm for steadystate systems, the initial conditions of the tubular system with a co-current cooling jacket (eq 1) have been determined. Figure 4a shows that the steady-state coolant temperature profile is determined by the minimization algorithm between the coolant device and the reactor’s surface, and Figure 4b-d demonstrates the dynamical responses of the illustrated tubular system. Notably, the base values of Tcf and Vc are confirmed, because of smooth, no-peak dynamic responses. Assuming that unmeasurable disturbances may appear in the feed of the tubular reactor, the feedback-based scheme is naturally required for the closed-loop performance and stability. From the process control viewpoint, the MPC scheme is an efficiently optimization-based feedback control strategy. Moreover, the aforementioned optimal operating policy should be integrated into the NMPC framework for the output regulation of the distributed reactor system while state/input constraints and unknown inlet disturbances are being considered simultaneously. In our study, the one-step-ahead predictive horizon is implemented to minimize the following objective function with no penalty on the control increments:
Ji ) [ym,i(k + 1) - ym,i(k) + yi(k) - yisp(k)]2 (for i ) 1, 2, ..., N) (22) Let zi represents the spatial position along the reactor length, and zi ) (L/N)i, for 0 e i e N. N is the number of spatial points, ym,i represents the model output at the location of zi, yi represents the process output at location zi, and yisp represents the reference output at the location of zi. Regarding the previous steady-state optimization design, the desired reference output is set by yisp ) T/ss(zi). Following that piecewise formulation, the originally distributed tubular system must be reduced to a lumped difference model:
xˆ i(k + 1) ) F(xˆ i(k),u(k))
(for i ) 1, 2, ..., N) (23)
with initial conditions
xˆ i(0) ) xˆ ss(zi)
(for i ) 1, 2, ..., N)
where the process states are
∫L (Tw - Tc,ss)2 dz V*, T* 0
min
Moreover, the initial conditions of the distributed coolant system (eq 1d) are obtained and written as
cf
xˆ i(k) } [CA(zi,k),CB(zi,k),T(zi,k),Tc(zi,k)]T
subject to and
where V* c and T* cf, as initials, are obtained using iterative computation through the minimization algorithm (eq 19).
xˆ ss(k) } [CA,ss(zi),CB,ss(zi),Tss(zi),Tc,ss(zi)]T Consider that the space discretization is based on the finite difference approximation; e.g., the five-point central difference formula29 is shown for
2070
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007
∂f(xˆ i) -f(xˆ i+2) + 8f(xˆ i+1) - 8f(xˆ i-1) + f(xˆ i-2) ≈ ∂z 12∆z )
∂f(xˆ i) ∂z
(24)
such that the difference-based nonlinear model described by eq 1 is shown at location zi:
C˙ A(zi,t) ) -
C˙ B(zi,t) ) -
T˙ (zi,t) ) -
(1 - ) V ∂CA(zi,t) +µ RA(CA(zi,t),T(zi,t)) ∂z (25a)
(1 - ) V ∂CB(z,t) +µ RB(CA(zi,t),CB(zi,t),T(zi,t)) ∂z (25b)
(1 - ) V ∂T(zi,t) +µ R ˜ (CA(zi,t),CB(zi,t),T Le ∂z FfcpfLe 2UL (zi,t)) + (T (z ,t) - T(zi,t)) (25c) r0FfcpfLe c i
T˙ c(zi,t) )
Vc∂Tc(z,t) 2UL + (T(zi,t) - Tc(zi,t)) V ∂z r′FccpcV
(25d)
Let Vc ) u1 and Tc(z0,t) ) Tcf ) u2. Moreover, the time discretization of eq 25 using an implicit Euler method with a sample time of ∆t will yield the discrete-time model shown in eq 23. Remark 6: The finite-difference method is a very simple technique, such that a large number of node points are used to improve the approximate accuracy and convergence properties. In the case of the predictive control approach by eq 22, a feedback controller is described by a set of difference equations,
ξi(k + 1) ) F(ξi(k),u(k))
(for i ) 1, 2, ..., N)
(26a)
u(k) ) Ψ(ξ1(k), ..., ξN(k); y1(k), ..., yN(k))
(26b)
N (∑j)1 Jj).
to minimize the lumped objective function The control law (eq 26) is an implicit function of the measured output at location zi, and the model state ξi corresponds to the positions of zi. In the present framework, the closed-loop system is formed by augmenting the open-loop system with an output feedback controller (eq 26). Remark 7: If the lumped difference model (eq 23) is an openloop stable system, i.e., a region of attraction is bounded by a positive constant Bd,
Figure 5. Proposed nondistributed model predictive control (NMPC) scheme: (a) the first control scheme, without sensing state information, and (b) the second control scheme, with sensing state information at the prescribed location.
because of the solvability of the minimization of the lumped objective functions. In addition, the numerous outputs (distributed output) feedback design may obviously increase the cost of the overall devices, and the no-offset tracking performance for multiple outputs are hardly achieved at the finite time. 4.1. Nondistributed Model Predictive Control (NMPC). According to the previous steady-state minimization approach for cooling jacket design, the base values of the two manipulated inputs (initials) are set by
N
∑ i)1
||xˆ i(k) - xˆ ss(zi)||eBd
(as k f ∞)
kf∞
(for i ) 1, 2, ..., N) (28) N ∑j)1 Jj,
the no-
(for i ) 1, 2, ..., N)
(29)
Using the previous minimization algorithm for offset tracking performance is obtained:
lim||yi(k + 1) - yisp(k)|| ) 0
kf∞
(30a)
u2(0) ) T*cf
(30b)
and
Furthermore, the asymptotic output regulation is achieved, i.e.,
lim||ym,i(k + 1) - ym,i(k)|| ) 0
u1(0) ) V*c (27)
In fact, the control law described by eq 26 is difficult to obtain,
and the discretized reference output yisp (for i ) 1, 2, ..., N) are determined by the previous steady-state Pontryagin’s minimization approach. Figure 4b-d has demonstrated the smooth, nopeak dynamic responses for the illustrated tubular system. For the reduction of the previous controller (eq 26), the one-stepahead prediction horizon is proposed for the specified single output at the outlet, with respect to amplitude and velocity
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2071
Figure 6. Disturbances rejection using the first NMPC scheme: (a) response of the exit concentration of B, (b) response of the exit reactor temperature, (c) the corresponding manipulated coolant flow rate, (d) the corresponding manipulated inlet coolant temperature, (e) three-dimensional plot of the reactor temperature response during a +20% change in the reactant concentration in the feed, and (f) three-dimensional plot of the reactor temperature response during a -20% change in the reactant concentration in the feed.
constraints on both manipulated inputs. The NMPC algorithm is reduced by
(ii) amplitude constraints on the manipulated input,
uj,min e uj(k) e uj,max min JN ) [∆ym,N(k + 1) + yN(k) - yN (k)] + sp
(for j ) 1, 2)
(33)
2
and (iii) velocity constraints on the manipulated input,
u(k)
2
kmpc
∆uj2(k) ∑ j)1
(31)
subject to (i) a nominal discrete-time model as the dynamic model constraint,
ξi(k + 1) ) F(ξi(k),u(k))
(for i ) 1, 2, ..., N)
ym,N(k) ) CNξN(k)
(32a) (32b)
∆uj,min e ∆uj(k) e ∆uj,max
(for j ) 1, 2)
(34)
where ∆ym,N(k + 1) ) ym,N(k + 1) - ym,N(k), ∆uj(k) ) uj(k) uj(k - 1), CN ) [0 0 1 0], kmpc is the input weight factor, and ∆uj represents the movement of the manipulated input increment. Remark 7: The proposed predictive control algorithm is a nondistributed output feedback controller. It manipulates a distributed reactor system using the steady-state optimization
2072
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007
Figure 7. Disturbances rejection using the first NMPC scheme: (a) response of the exit concentration of B, (b) response of the exit reactor temperature, (c) the corresponding manipulated coolant flow rate, (d) the corresponding manipulated inlet coolant temperature, (e) three-dimensional plot of the reactor temperature response during +1% changes in the coolant temperature in the feed, and (f) three-dimensional plot of the reactor temperature response during a -1% change in the coolant temperature in the feed.
approach and an open-loop observer. Moreover, the proposed model-based control scheme is depicted in Figure 5a. Note that the lumped difference model as the output prediction, with respect to the input constraints, is treated as the feedback-based implementation scheme. The step disturbance rejection is achieved due to the error between the process output and the model prediction, i.e., yN(k) - ym,N(k), added in the prescribed objective function (eq 31). When the inlet disturbances are added, the previous steadystate approach cannot seize the regular dynamic behavior, e.g., the location of the peak temperature. In addition, the state constraint by eq 11 for the exothermic reactor system is strictly required. Obviously, the distributed state response will degrade the former NMPC scheme, because of the appearance of a moving peak temperature and the inaccurate state information
by the tubular systems. Thus, a measurement-based NMPC algorithm is proposed and written as M
min JM ) u(k)
∆yˆ m,l2(k + 1) + ∑ l)1 2
[yN(k) - yNsp(k)]2 + kmpc
∆uj2(k) ∑ j)1
(35)
subject to (i) the lumped difference process model,
xˆ i(k + 1) ) F(xˆ i(k),u(k)) + η(d(k)) yˆ m,l(k) ) Clxˆ l(k)
(for i ) 1, 2, ..., N) (36a)
(for l ) 1, 2, ..., M)
(36b)
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2073
(ii) the state constraint,
φ(xˆ i-1(k)) e φ(xˆ i(k))
(for i ) 2, 3, ..., N)
(37)
and (iii) the same conditions for amplitude and velocity constraints on the manipulated input. In eqs 35-37, ∆yˆ m,l(k + 1) ) yˆ m,l(k + 1) - yˆ m,l(k), yˆ m,l represents the measurable output using sensor measurements at the prescribed position zl, η represents the effect of unknown disturbances d, and M is the number of spatial points that settle on measurement devices. Remark 8: With the aid of a few sensors for measurement, the measurement-based predictive control algorithm controller design can induce the stable and no-offset output regulation at the outlet of the tubular reactor. Moreover, the proposed control scheme is depicted in Figure 5b. In addition, the number of sensors M should be less than the number of spatial points N for lumped difference equations. If the number of sensors is
sufficient, the prescribed state constraint (eq 37) for suppressing an undesired peak temperature could be satisfied. Remark 9: Assume that the nominal solution u(k) ) Ψ(ξ1(k), ..., ξN(k); y1(k), ..., yN(k)) to the minimization of N ∑j)1 Jj is obtained from the system without disturbances, e.g., xˆ i(k + 1) ) F(xˆ i(k),u(k)) (for i ) 1, 2, ..., N). Apparently, the feasibility of the aforementioned optimization algorithm in the face of physical constraints would be dependent on the controller parameter kmpc. It implies that the feedback control u(ym,N, yN, kmpc) can be determined by solving the aforementioned optimization problem (eq 31). Moreover, if the optimization problem (eq 31) is feasible, eq 35 also has a feasible solution, because of the aid of measurements, i.e., min JM e min JN. 4.2. Demonstration. The presented simulation tests are based on the following conditions: (a) The number of spatial points N for lumped difference equations is 100.
Figure 8. Disturbances rejection using the second NMPC scheme: (a) response of the exit concentration of B, (b) response of the exit reactor temperature, (c) the corresponding manipulated coolant flow rate, (d) the corresponding manipulated inlet coolant temperature, (e) three-dimensional plot of the reactor temperature response during a +20% change in the reactant concentration in the feed, and (f) three-dimensional plot of the reactor temperature response during a -20% change in the reactant concentration in the feed.
2074
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007
Figure 9. Disturbances rejection using the second NMPC scheme: (a) response of the exit concentration of B, (b) response of the exit reactor temperature, (c) the corresponding manipulated coolant flow rate, (d) the corresponding manipulated inlet coolant temperature, (e) three-dimensional plot of the reactor temperature response during a +1% change in the coolant temperature in the feed, and (f) three-dimensional plot of the reactor temperature response during a -1% change in the coolant temperature in the feed.
(b) The number of spatial points for sensor measurements is 20. (c) The bounds for the input constraints include (i) an amplitude constraint,
0.01 m/s e Vc e 0.05 m/s
(38a)
605 K e Tcf e 625 K
(38b)
and (ii) a velocity constraint,
and
u2(0) ) T*cf ) 614.3 K and (ii) the desired setpoint as the reference output,
T*ss(zN) ) yNsp ) 630.5 K
(41)
(e) The optimal state distribution (Figure 4) for the differencebased model (eq 25) is given as
0 m/s e ∆Vc e 10-4 m/s
(39a)
0 K e ∆Tcf e 0.1 K
(39b)
(d) The optimal input and output (steady-state optimization approach) includes (i) the initials of input,
u1(0) ) V*c ) 0.02 m/s
(40b)
(40a)
and
CA(zi,0) ) C*A,ss(zi)
(42a)
CB(zi,0) ) C*B,ss(zi)
(42b)
T(zi,0) ) T*ss(zi)
(42c)
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2075
Tc(zi,0) ) T*c,ss(zi)
(42d)
(f) The undesired operating temperature is >640 K. (g) The sample time ∆t is 1 s. (h) The tuning parameter kmpc is fixed (kmpc ) 0.01). According to the first NMPC scheme shown in Figure 5a, it will be validated by the illustrated reactor system while the state/ input constraints and inlet disturbances are being considered simultaneously. For inlet disturbances, changes of (20% in the reactant concentration in the feed, and (1% in the coolant temperature in the feed are considered in the process. Figure 6b shows that the NMPC with only output information at the outlet temperature Tout can ensure the stable and no-offset output regulation in the presence of (20%CAf. Notably, Figures 6c and 6d show both actuators Vc and Tcf being restricted under physical constraints, and the three-dimensional figures shown in Figures 6e and 6f show the complete state distribution. However, the perturbation by +20%CAf will induce the peak temperature close to the inlet bound. Under the same NMPC scheme for (1%Tcf, Figure 7b shows that the stable and nooffset output regulation can be achieved. The corresponding two controller actions are depicted in Figures 7c and 7d, and the three-dimensional temperature responses are shown in Figures 7e and 7f. Unfortunately, the perturbation by +1%Tcf will obviously induce the undesired peak temperature (>640 K) close to the inlet bound. In our opinion, this (internal) model-based control scheme can be easily realized through the fast online computation; however, the nondistributed and output feedback design cannot effectively dominate the entire temperature effect inside the reactor. According to the second NMPC scheme shown in Figure 5b, it will be validated by the illustrated reactor system with the same state/input constraints and inlet disturbances. Figure 8b shows that the NMPC with 20 sensors for measurement of the reactor temperature at the fixed interval can ensure the stable and no-offset output regulation in the presence of (20%CAf. The corresponding actuators Vc and Tcf shown in Figures 8c and 8d are being restricted under physical constraints, and the complete state distributions are shown in Figures 8e and 8f. Compared to Figure 7b, the second NMPC scheme can provide better tracking performance than the first scheme, but the moreoscillating responses are detected using the second scheme. Compared to Figure 7e, the second NMPC scheme can provide a smooth temperature distribution, because of the multisensor design. Under the same NMPC scheme for (1%Tcf, Figure 9b shows that the stable and no-offset output regulation can be achieved. The corresponding two controller actions are depicted in Figures 9c and 9d, and the three-dimensional temperature responses are shown in Figures 9e and 9f. Compared to Figure 8b, it verifies again that the second NMPC scheme can provide better tracking performance than the first scheme, and the second NMPC scheme can almost reduce the peak temperature, as shown in Figure 9e. Consequently, the proposed two NMPC schemes have been successfully validated for controlling a distributed system in the face of inlet disturbances or state/input constraints. Particularly, the second NMPC scheme almost captures the characteristic of distributed parameter systems, because of the number of measurement devices. Based on the accurate steady-state optimization approach, the aforementioned two input control schemes can ensure the stable/robust output regulation around prescribed operating conditions.
5. Conclusion For the illustrative exothermic tubular reactor system, the problems of hot spots, partial differential equation (PDE) models, boundary control, and physical constraints are challenging to traditional feedback control design. Fortunately, Pontryagin’s minimum principle, in association with the extremal control, facilitates the steady-state optimization approach, so that the modified model predictive control (MPC) strategy, using several output measurements, is added to optimize the processing temperature profile and ensure the largest amount of products in a fixed space time. The two control actuators are described as boundary control of distributed parameter systems to dominate the heat exchange ability of the cooling device. In regard to state/input constraints handling and the optimization resolution, an exterior penalty method is used to account for the controller parameter, with respect to closed-loop feasibility during on-line iterative procedures. All nondistributed model predictive control (NMPC) strategies are reduced to a one-stepahead predictive horizon. The simplest tuning procedure in the output feedback framework is quite easy to implement in practice. In regard to the control of hot spots, the proposed measurement-based NMPC is implemented to reduce the undesired peak temperature. Although the position and numbers of sensors aren’t explored by theoretical analyses, simulation tests indicate that the measurement sensors at the hot-spot zone can effectively remove the hot-spot temperature. Acknowledgment This work is supported by the National Science Council of the Republic of China, under Grant No. NSC-94-2214-E-224007. Notation CA ) concentration of A (kmol/m3) CAf ) feed concentration of A (kmol/m3) CB ) concentration of B (kmol/m3) CBf ) feed concentration of B (kmol/m3) CB,out ) exit product concentration (kmol/m3) cpc ) heat capacity of coolant (J kg-1 K-1) cpf ) heat capacity of fluid (J kg-1 K-1) cps ) heat capacity of catalyst (J kg-1 K-1) Ei ) activation energy of reaction i, for i ) 1, 2, 3 (J/kmol) -∆Hi ) heat of reaction i, for i ) 1, 2, 3 (J/kmol) L ) reactor length (m) r0 ) reactor radius (m) r1 ) coolant tube radius (m) Rg ) gas constant Ri ) rate of reaction i, for i ) 1, 2, 3 t ) time (s) T ) reactor temperature (K) Tf ) feed temperature (K) Tc ) coolant temperature (K) Tcf ) inlet coolant temperature (K) Tw ) wall temperature (K) Tout ) exit reactor temperature (K) Tr,max ) hot-spot temperature (K) ∆t ) sample time U ) overall heat transfer coefficients (J kg-1 K-1 s-1) V ) feed flow rate (m/s) Vc ) coolant flow rate (m/s) z ) axial distance (m) ) void volume fraction of reactor µ ) catalyst activity
2076
Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007
Fc ) coolant density (kg/m3) Ff ) fluid density (kg/m3) Fs ) catalyst density (kg/m3) H ) Hamiltonian function λ ) adjoint vector Superscripts * ) optimal sp ) set point Subscripts ss ) steady state Literature Cited (1) Wu, W. Finite difference output feedback control of a class of distributed parameter processes. Ind. Eng. Chem. Res. 2000, 39, 4250. (2) Christofides, P. D. Nonlinear and Robust Control of Partial Differential Equation Systems: Methods and Applications to TransportReaction Processes; Birkauser, Boston, 2001. (3) Hudon, N.; Perrier, M.; Guay, M.; Dochain, D. Adaptive extremum seeking control of a non-isothermal tubular reactor with unknown kinetics. Comput. Chem. Eng. 2005, 29, 839. (4) Wu, H.; Morbidelli, M.; Varma, A. Pseudo-adiabatic operation and runaway in tubular reactors. AIChE J. 1998, 44, 1157. (5) Karafyllis, I.; Daoutidis, P. Control of hot spots in plug flow reactors. Comput. Chem. Eng. 2002, 26, 1087. (6) Wu, W.; Huang, M. Y. Nonlinear inferential control for an exothermic packed-bed reactor: geometric approaches. Chem. Eng. Sci. 2003, 58, 2023. (7) Rossiter J. A. Model-Based PredictiVe Control: A Practical Approach; CRC Press: Boca Raton, FL, 2003. (8) Shang, H.; Forbes, J. F.; Guay, M. Model predictive control for quasilinear hyperbolic distributed parameter systems. Ind. Eng. Chem. Res. 2004, 43, 2140. (9) Dufour, P.; Michaud, D. J.; Toure, Y.; Dhurjati, P. S. A partial differential equation model predictive control strategy: application to autoclave composite processing. Comput. Chem. Eng. 2004, 28, 545. (10) Dufour, P.; Toure, Y. Multivariable model predictive control of a catalytic reverse flow reactor. Comput. Chem. Eng. 2004, 28, 2259. (11) Dubljevic, S.; El-Farra, N. H.; Mhaskar, P.; Christofides, P. D. Predictive control of parabolic PDEs with state and control constraints. Int. J. Robust Nonlinear Control 2006, 16, 749. (12) Dubljevic, S.; Christofides, P. D. Predictive control of parabolic PDEs with boundary control actuation. Chem. Eng. Sci. 2006, 61, 6239.
(13) Rahman, S.; Palanki, S. On-line optimization of batch processes in the presence of measurable disturbances. Ind. Eng. Chem. Res. 1996, 35, 3567. (14) Visser, E.; Srinivasan, B.; Palanki, S.; Bonvin, D. A feedbackbased implementation scheme for batch process optimization. J. Process Control 2000, 10, 399. (15) Smets, I. Y.; Johan, E. C.; November, E. J.; Bastin, G. P.; Van Impe, J. F. Optimal adaptive control of (bio)chemical reactors: past, present and future. J. Process Control 2004, 14, 795. (16) Welz, C.; Srinivasan, B.; Marchetti, A.; Bonvin, D. Evaluation of input parameterization for batch process optimization. AIChE J. 2006, 52, 3155. (17) Smets, H. Y.; Dochain, D.; Van Impe, J. F. Optimal temperature control of a steady-state exothermic plug-flow reactor. AIChE J. 2002, 48, 279. (18) Birk, J.; Liepelt, M.; Schittkowski, K.; Vogel, F. Computation of optimal feed rates and operation intervals for tubular reactors. J. Process Control 1999, 9, 325. (19) Boskovic, D. M.; Krstic, M. Backstepping control of chemical tubular reactors. Comput. Chem. Eng. 2002, 26, 1077. (20) Vecchio, D. D.; Petit, N. Boundary control for an industrial underactuated tubular chemical reactor. J. Process Control 2005, 15, 771. (21) Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Wiley: New York, 1975. (22) Christofides, P. D.; Daoutidis, P. Feedback control of hyperbolic PDE systems. AIChE J. 1996, 42, 3063. (23) Wu, W.; Liou, C. T. Output Regulation of Two-time-scale Hyperbolic PDE Systems. J. Process Control 2001, 11, 637. (24) Froment, G. F. Fixed-bed Catalytic Reaction. Ind. Eng. Chem. Res. 1967, 59, 18. (25) Chen, C. Y.; Sun, C. C. Adaptive Inferential Control of PackedBed Reactors. Chem. Eng. Sci. 1991, 46, 1041. (26) Hua, X.; Jutan, A. Nonlinear inferential cascade control of exothermic fixed-bed reactors. AIChE J. 2000, 36, 980. (27) Smets, I. Y.; Van Impe, J. F. Optimal Control of (Bio)Chemical Reactors: Generic Properties of Time and Space Dependent Optimization. Math. Comput. Simul. 2002, 60, 475. (28) Wu, W.; Liou, C. T. Output regulation of nonisothermal plug-flow reactors. Comput. Chem. Eng. 2001, 25, 433. (29) Constantinides, A.; Mostoufi, N. Numerical Methods for Chemical Engineers with MATLAB Applications; Prentice Hall: Englewood Cliffs, NJ, 1999.
ReceiVed for reView August 28, 2006 ReVised manuscript receiVed November 30, 2006 Accepted February 5, 2007 IE0611296