- Email: [email protected]
Optimal cold-shot cooling of an adiabatic tubular reactor
Pergamon Press.
Chemical Engineering Science, 1969, Vol. 24, pp. 279-289.
Printed in Great Britain.
Optimal cold-shot cooling of an adiabatic tubular reactor L. PADMANABHAN
and S. G. BANKOFF
Department of Chemical Engineering, Northwestern University, Evanston, Illinois, U.S.A. (First received 2 August 1968; in revisedform 2OAugust 1968) Abstract - Optimal design of multi-bed adiabatic reactors with cold-shot cooling is considered for the case of several independent chemical reactions. The necessary equations for optimality are derived by applying Pontryagin’s methods to the variational problem. A computational scheme based on Newton’s method in the control-parameter space is devised and applied to the consecutive reaction system A + B 4 C in a tubular reactor with a single cold-shot injection. Upper and lower bounds in the performance index of such reactors are computed by analyzing limiting configurations. A sensitivity analvsis leads to a scheme for devising an optimal feedforward controller that compensates for small perturbations in initial conditions. INTRODUCTION
THIS paper focuses its attention on a class of unconstrained optimal control problems associated with systems exhibiting discontinuous dynamics. By this we mean that the direction field associated with the differential equations of motion suffers a finite number of discontinuities and the jumps at such points are given by continuous functions of the state variables. Such discontinuities divide the composite system into a chain of staged subsystems, each subsystem having its own characteristic equation of motion. Frequently, such d&continuities are deliberately introduced into the system in order to provide feasible approximations to some complex, perhaps unattainable, optimal strategy. Analysis of such problems is of immense value to the chemical engineer concerned with optimal design of complex plants. Confronted with the prohibitive costs of implementing some optimal control configurations on the one hand, and the unacceptably poor performance of extremely simple non-optimal or sub-optimal modes of operation on the other, the designer is compelled to seek a compromise which may lead to discrete approximations to an optimal control function. A case in point is the task of designing a tubular reactor so as to maximise the outlet concentration of a valuable product. The kinetics associated with this problem may involve several exothermic reactions.
The optimal policy is to maintain a temperature profile that favors the desired reactions leading to the valuable product and suppresses the dissipative reactions. One can easily compute this optimal temperature program by techniques which are now welldeveloped[ 11. Theoretically attractive as this may be, it is beset with the difficulty of implementing such a temperature control along the reactor. Equally difficult to achieve in practice is the idealized continuous cross-feed reactor studied in detail by Horn and others[2]. However, the yields from these reactors represent upper bounds for the performance index of a variety of other reactor configurations. On the other extreme among the optimal reactors is the adiabatic tubular reactor with the feed temperature maintained at an optimal value. This system is the simplest in structure and the easiest to accomplish, but may have a poor performance, mainly due to the fact that the temperature rise near the exit in an adiabatic reactor tends to favor the side reactions, which were successfully suppressed near the entrance. The logical answer to this is the compromise afforded by the multi-bed adiabatic reactor chain with intermittent cooling[3, 41. This cooling can be accomplished indirectly through a heat exchanger, or directly via coldshot injection, wherein a fraction of relatively cold feed stock by-passes the first few beds. The method of cold-shot cooling dispenses with
279
L. PADMANABHAN
and S. G. BANKOFF
the intermediate heat exchanger equipment and therefore lowers the capital and operating costs. This makes it attractive despite the unavoidable dilution of the enriched reaction mixture. The possibility of achieving an optimal design of the cold-shot injection and the by-pass system, so as to diminish the dilution effects, is a real incentive for the study of such reactors. In this paper we attempt to answer the problems associated with this. For the case of a single reaction, the optimal design of multi-bed adiabatic reactors has been examined in depth by Aris[3] who uses dynamic programming to evolve elegant graphical procedures that provide rapid solutions to the optimization problem. Such techniques are unique to the case of a single reaction, where temperature and conversion are the only state variables, and therefore can be depicted on a two-dimensional state plane. The effectiveness of such geometrical construction derives, in part, from the fact that, in order to maximize the yield from a single-reaction system, one simply has to maximize the reaction rate everywhere. This enormous simplification is not available when there are several competing or consecutive reactions. Although the dynamic programming concept is still applicable, its embedding approach to such problems becomes computationally inefficient and expensive. The reason for this is rather dramatically termed by Bellman the “curse of dimensionality” [5]. For the same reason the use of graphical methods is precluded. One is therefore compelled to seek analytical methods which will be computationally effective with more complex systems. One approach is to employ Newton’s first-variational scheme on the variational equations derived from the application of Pontryagin’s formulation. This will be elaborated in detail below. STATEMENT
OF
THE
PROBLEM
Let us consider a plug-flow adiabatic reactor of specified length tf, with a cold-shot injection at t, < tp The treatment for the case when tf is free can be performed easily by the techniques employed here. The raw material has composition
Cj,,, j = 1, 2, m (with m species) and temperature T,,. The total feed to the reactor is F,, (in volumetric rate units). The reactions are accompanied by evolution of heat, and the consequent temperature rise is to be curbed by ideally locating the cold-shot injection point, tl. The injection point divides the reactor into two regions: the hot pre-injection zone (0 d t < tJ and the relatively cold post-injection zone (tl < t s tf), each governed by different state equations. The discontinuity at tl is determined by appropriate material and energy balances. [The axial distance is specified in terms of the residence time t, based on the total feed F,]. Since the injection point is unspecified, the overall problem breaks down into two free-time control problems for the two sections of the reactor. In an earlier paper[6] we have established effective methods for handling free-time problems using the concept of “reduced-time domain analysis”, which provides’a method of reducing the free-time problem into a fixed-time problem. We will use this concept in what follows. Returning to the reactor design under study, the optimization is to be introduced via parametric control. We use this term to emphasize the fact that the decision vector consists of the following three parameters: the inlet temperature (or the preheat temperature), T(O), the bypass ratio, ( 1 - y) , and the location of the cold-shot injection point, tl. The preheating of the feed from To to T(0) enhances production rates in the pre-injection section. We have a genuine optimization problem in the choice of the bypass ratio, as it should be large enough to quench the hot reaction mixture, but small enough to avoid excessive dilution. We now formulate the problem for the case of IZ independent reactions. It is convenient to express the local concentrations Cj in dimensionless form by referring them to the feed and writing xj = CJ 3: Cjo, SO that T xi, = 1. If we pick II of these xj’s as a basis, then we have a state vector xeEn. In addition the local temperature T is a variable. In order to distinguish between the variables at the pre-injection and the postinjection zones; we will refer to the former by
280
Optimal cold-shot cooling of an adiabatic tubular reactor
attaching by a (+) From following
dr dt=
a (-) sign to them and to the latter sign. a local material balance we get the state equations for the two zones:
ys=
1 y&(t),
“:Xx-,
T(0) +J(x-(t) --x0)1 T(O),y)
^;f-
(10)
T-)
g(x-,
(t < t,>
(1) +(I-_y)To
F=g(x+(t),yT(O)
5 =g(x+,
T+)
(2)
(t > t1)
+J(x+(t)
’
-
T-)
(t ==ct1)
(3)
dT+ = Jg(X+, T+) dt
(t > t,)
(4)
dt
= Jg(x-,
where J E En has components Ji = -AHJC, and the inner product is implied. Here (-AHi) is the heat generation accompanying unit increase in Xi, and C, is the volumetric heat capacity of the reaction mixture. For simplicity we shall take J to be a constant vector. It then follows from Eqs. (l)-(4) that T-(t)
= T(0) +J(x-(t)
T+(t)
= T+(t,)
+J(x+(t)
-x0)
= rx-(t1)
L-(x+,
T(O), y) L f+ (11)
We will restrict y to the open interval (0, 1). The cases y= 0, 1 refer to the simple adiabatic reactor to be considered later. If the objective function to be maximized is @(xf, T(O), y), wherex,=x($) and@ E C2(En, El, El), then the necessary equations for optimality are derived in the following manner. Define the Lagrangian L by L=Q(x,,
T(O),y)+l”-
h-(t)(F-g}dt
+I,”
“+(t){f+-F}
dt
(12)
(5)
-x+(tl))
where x+ (tl) and T+ (tI) are determined jump balances at t = tl: X’(h)
-xc,))
(t1 < t c tf).
where g E C2 (x; T): En+’ --f En. Similarly we get from an energy balance: dT-
(0 c t < t1)
where the adjoint variables AkeEn form the Lagrange multipliers for the differential constraints (10) and (11). Following the standard procedure, one takes variations of all the terms in (12) in order to get the first variation of L.
(6)
through
+ (1 -Yh
(7)
T+(t,) = yT-(t,)+(l-y)T,.
(8) +
Gx-dt
Here x0 is the vector of feed concentration and To is the temperature of the raw material. The latter are simply mixing equations for the injection of the bypass stream (I- y)Fo into the through feed yFo. Using (7) and (8) in (6), we get:
- h-(t,)
ifix-
Gx+dt - A+(t f )8x+(t f )
A-f-
T+(t)
= yT(0)
+ (l-y)T,+J(x+(t)
-x,,).
(9)
Now Eqs. (5) and (9) enable one to eliminate the temperature variables from (1) and (2).
-~Sy+~g~-6T(O)
+
+
l,;
(A+g,+(T(O) - To)6y+yg,+A+6T(0)}dt. (13)
281
L. PADMANABHAN
and
It follows from (7) that A-r+&) = rh-0,)
ordinate 8 and time-scaling factors p’ defined by: + (x-(G-
(14) /
-My
where Ax*&) is the total variation in x*(tJ made up of 6x*( tl) which is the local variation in x* at fixed t1 and of (dr’ldt) I& , the contribution due to a variation in t,: h’(M
A 6x+(1,) +f(X+(tl),
Ax-(t,) L 8x-(1,)++x-(r,),
T(O),y)
(19)
A-01) = rA’(t1)
(20)
(25)
tr= (p--p+)e,+p+e,.
(26)
and
(27)
so that (26) implies that p-+/3+
= 2.
(28)
Also from 0 < c1 < tr, it follows that 0 < p’ c 2. We now recast the various equations in the reduced-time domain. First, we write the state Eqs. (10) and (11) as integral equations:
x-(p-e)
=x,+1
y op-efw(t~, I
x+(q+(e))=yx-(p-el)+
wv9
7) d5
(29)
(1-Y)-G T(O), y) de.
(30)
(21) Noting that aq+ (0)/ a0 = p+, we obtain from (29) and (30)
Q 2 A+(?,) (x-(tl) -x0) +@,-I y 2 :- A-f- dt J ’ A+gT+dt=O I t,+
t1 = p-e,
+ j;:;‘f(x+(&
= 0
(24)
so that
-f(x-(G,
T(O),y))
0 > tl) 1
8, = t,/2 = ef/2
(18)
A’(b) = d$
+(T(O)-Z-,,)
71-(d) = P-R (I < t1) [ 71+(e) = (p--p+) e,+B+e,
If we take 8, to be fixed, then the variability of t1 is reflected totally in that of 0’ through (24) and (25). e1 and 0, can be arbitrary chosen, and a natural and convenient choice is
St,. (16)
(rl -c t =sf,)
T(O),y)
t=
(15)
T(O), r> at,
In the absence of any inequality constraints, maximization of @ implies that SL be zero, which in turn requires that all the independent variations in (13)‘vanish individually. Using (14), (15) and (16) in (13) and assembling the various vanishing terms, we arrive at the necessary conditions:
s A A+(&) {f(X’(G,
S. G. BANKOFF
ax-=ae
(22)
p-
_
YfT
R 4 yfi’m,, + J;‘-A-gT-df + y* J,;+ A+gT+dt = 0.
(e < 0,)
(31)
(0 > 0,).
(32)
(23) (21), (22), and (23) are the determining tions for tl, y and T(0) respectively. COMPUTATIONAL
equa-
SCHEME
Along the same lines, the adjoint equations (17) and (18) can be written as: A-(P-0)
The computational method is easily formulated through the introduction of a reduced-time co282
= rA+(P-0,)
Optimalcold-shotcooling of an adiabatic tubular reactor sEc= (b-p,
T(O),Y) d5
(34)
(0 < 01)
(35)
(6 ’ h).
(36)
+r,
f(o) - T(o))T
(42)
and JI = (&, &. I/J~)~with the elements given by:
whence we have
ax-= ae
+I= S+AA+(e1)(f+--f-)It=tl
-5
A-f,-
!!$ = - p+
A+fz+
+A+(t,)(Yf,+--f,-)l,=t,~-(e,) $2
=
Q+W’,z, %+AA+(&)(x-(h) +A+(h)
We will also need the reduced-time version of the determining equations (21), (22) and (23): S = A+(p-0,) {f(x+(P-e,), -f(x-(P-e,), Q = h+(p-8,)
--d
h-w {A-f,- hx-(6’) +AA-(e).L-}
de
A+yAx+(e) +(T(0)-To)P+ly ox+ i
T(O), Y)
T(O), Y)] = 0
(x-(p-el)
(43)
(37) +AA+(@g,+
-x0> -$/“-
I
de
(44) e,-
$3 = R + Y@&o,,~,&+
xh-f-dB+Q,+(T(O)--T.)p+y;+
P-
I
0
x A+gT+de = 0
X AA-(@ g,-+A-d$Lr-(8))
(38)
R = /3- J;‘- A-gT- de+ y2/I+ JO; A+g,+ de + YJ%YO~ = 0.
de
AA+(e) gT++Atd$&+(8)) (39)
The proposed computational scheme is essentially Newton’s method applied to S = Q = R = 0 in the (p-, y. T(0)) parameter space. At any point CL= (p-, Q, ?(0))T, we have the following first-order Taylor expansion for 6= (5, $9 2). fi=
R+n,s/_&
(40)
where 0, = gradJl and 6, = fi - CL. But from (37), (38) and (39), we have 6 = 0, so that a first-order estimate of the location of the optimum fi can be obtained by solving (40). To this end, we take variations in S, Q, R, substitute these in (40), after making use of (14) and finally get the simplified linear system of variational equations,
(45) where
b*(e)
.
= a-c(+y(e))
--xyq(e))
and AA’(O) = fi*(+j*(e)) -A’(r)‘(@): (Note that Ax+(e) involves a comparison of the state variable at the same reduced-time coordinate and not at the same real-time. This is the main feature in the reduced-time domain analysis, as the quantities Ax’( 0)) AA’( 0) are more readily obtained than F(t) -x’(t), A*(t) -A’(t), for the latter will demand use of interpolation schemes. This is because the time-scaling factors p’ map the variable intervals [0, r,], [tI, rr] into ones of fixed length, viz., [0, e,] , rh
ASp=J,
de
(41)
where A is a 3 X 3 matrix to be defined below,
e,i).
The matrix A has elements listed below: All=0
283
(46)
L. PADMANABHAN AI2 = -A+(b)
t.fr+(x-(t)
and S. G. BANKOFF
- xo) + V(O) - To)g,+I It=t,
(47) A13 =Y--h+(h) {WTf - gT-lIt=t* A,, = $
(48)
“- A-f- de+ ( T(0) - To) I0 6
X
h+g,+
de
(49)
(35), (36) backward in the reduced-time domain, using the end-condition (19) and the jump condition (20). 3. Compute P, Qck’ and Rck’. If these are smaller in absolute value than some preassigned tolerances, the computation is terminated. Otherwise proceed to step 4. 4. If k = 0, set
I e,+
A 22 =__2pY3
e’-A-f-dt9-_rB+(T(0)-Z-o)2 I0 eI X A+& df3- G’,, I
P
I
A+ [gT++Y(T(O)
- To)gbl
(50)
yuc)
(56)
@-@mo)
e,+ (51)
AB1 = - J,“’ A-g,- de + y2 Je; A+g,+ de AX = - (@‘no, + vhom
) -?‘p+
je.B
=
-
Y%o~~o~
-
p-
Joe+-
P (k+l)
s
e’
8,
Once these quantities are yields the correction in control: 8/L= EA-‘JI
Ag& de. evaluated
,Jk’+
E(k) {A(k)}-‘q,‘k’.
(57)
If the matrix Atk’ is ill-conditioned, then use (56). 5. Increase k by unity and go to step 2. THE CASE
(53)
WHEN
y = 0,l
This corresponds to a simple adiabatic reactor without any cold-shot cooling. The only freedom we have is in the choice of the preheat temperature, T(0). For this case, the state equation is simply
A-g, de - y”/l+ X
=
(52)
A+(2&+
+ y(T(0) - To)gzG} de A33
p(k)+
where P, i = 1,2,3 are small arbitrary increments in /3-, y, T(O). If k a 1, set x = tik), A = Ack), Ax’(e) = x”(~) (nfi”c’(0)) -x’(k-i)(rl+(k-l)(e)), and use a similar result for AA’(e). Compute the vector $‘k’ and the matrix A(“). Choose a suitable l E (0, 1) and obtain the correction from:
A,, = P- ‘I- A-g,- de Ti+ Y2I 0 0,
X
(k+l) =
(54) (41)
$ = g(x, T(0) +J(x-x0)) (55)
where E, 0 < E < 1, is introduced to cut down undesirable overshoots in correction. It is now clear that (55) can be employed to develop a suitable iteration process that converges to the optimum. The algorithm is outlined below: 1. Set k = 0. Choose p-(O), y(O)and PO’(O). 2. Using x-(O) =x0 as the initial condition, integrate the state equation (3 1) forward over the pre-injection section. Then, using the jump condition (7), integrate (3 1) until the exit 8 = 0, is reached. Now, integrate the adjoint equations
(58)
with the initial condition x(O) = x,,. The adjoint equation is easily obtained as: (59) with A(tf) = CD+. The determining T(0) can be shown to be:
equation
U : aRo, + ,f Ag,,, dt = 0..
for
(60)
The computational algorithm for finding T(0) is obtained, as before, by taking variations in U and setting 0 = 0, whereupon we have:
284
Optimalcold-shotcoolingof an adiabatictubularreactor
AN
ILLUSTRATIVE
EXAMPLES
The theory developed so far will now be applied to the study of optimization for the consecutive reaction system A + B --, C conducted in a tubular reactor. The objective here is to maximize
(62) kz(T)xZ
(63)
where k,(T) = kN exp (--E,IT), j = 1,2. Using (5) and (9), we eliminate the temperature dependence from (62) and (63) to obtain c&x-, T(O) +1(x--x,)1, f(x)={,(x+,yT(0)+(l--y)To
(t < tl) (64) (t ’
+Jb+--x,)1,
h).
For numerical work, the various constants appearing in Eqs. (62), (63) and (64) are chosen as follows: k,, = 0.5 x 108,
k,, = 0.5 x 1016
E, = 6900 (OK)-‘,
Ez = 13800 (OK)-’
J, =-32
J2 = 10 (“K)
x 10= 0.95,
(OK),
X 20 =
0.05.
A good choice of the nominal control policy is essential for convergence of the iterations. The nominal control used in the example was: p-(O) = 1.75, ~‘0’= 0.9, T(O)(O)= 350°K.
Starting from these values, the control vector progressively converged in 40 iterations to the optimum: fl- = O-8575, p = O-9151, f(O) = 34OeO”K. The objective function @ at this point had the value $ ( tf) = O-5176, as compared with xz(tf) = 0.2852 corresponding to the nominal policy. A fourth-order Runge-Kutta method was used to integrate the differential equations, partitioning the reactor length tf = 10 into 40 equal intervals. The various integrals appearing in Eqs. (43) to (54) were evaluated using Simpson’s rule, with a partition of the reactor length into 80 equal intervals. The computation, performed on a CDC6400 computer, took 0.57 set of central processor time per iteration. Also for the same reaction system, the optimal design of the simple adiabatic reactor without cold-shot cooling was considered. The computational procedure using Eq. (61) was employed to determine the optimal preheat temperature. Convergence was observed in 20 iterations to the optimal inlet temperature T(0) = 338”K, and the optimal yield was x,(tf) = O-5134, which, as expected, is less than the yield from the reactor with cold-shot cooling. The computation time was 0.43 set per iteration. The yield of the reactor with cold-shot cooling can be improved further by increasing the number of the injection points, which give more degrees of freedom. As the number of injection points increases, the yield of the reactor increases monotonically and ultimately attains an upper bound in performance index corresponding to the idealized continuous cross-feed reactor. However, the computational problems associated with this limiting configuration are not easily tractable, particularly in multi-reaction systems. This is due to the possibility of singular control. A more easily computed upper bound is the yield from the optimal temperature-programmed reactor. For the consecutive reaction system under study, the optimal yield from the programmed reactor was found to be x, (t,) = 0.5452. By comparing this figure with the yield from the simple adiabatic reactor, one can get a rough estimate of the improvement in yield to be expected from more
285
L. PADMANABHAN
and S. G. BANKOFF
cold-shot injection points. This improvement will be significant for those reaction systems that are strongly exothermic. LINEARIZED
(tl < t s tr) x+(t,)
FEEDFORWARD
CONTROLLER
In the course of the computation process for the reactor with cold-shot cooling, it was found that the yield was extremely sensitive to variations in the preheat temperature, bypass ratio, and feed composition. This makes it imperative to incorporate into the reactor system a feedforward controller that continually responds to perturbations in the initial conditions by providing compensating changes in the control parameters, thus restoring the system to an optimal configuration. Our approach to this problem is to approximate the surfaces Q = 0, R = 0 by their tangent planes at the point (x$, f * (0)) Y* ) . The asterisk symbol refers to the fact that the optimum (Z?(0)) 9* ) , corresponding to a particular choice of x$, is predetermined. In the neighbourhood of this point, the tangent planes are valid approximations to the optimal response surfaces. The determination of such tangent planes involve the calculation of the sensitivity of Q and R ‘with respect to x0, T(0) and y. For this, we go to the integral equations (29) and (30) and obtain by differentiation with respect to the initial condition, the sensitivity matrix X*(t) k (ax’(t)/(ax,):
X-(t) =
I+;,‘,(.&-X-(5) +.&Id5
-
dt
+ (1 --y)Z.
(69)
Analogous results for the sensitivity of the adjoint variable A’ (t) are obtained from (33) and (34). Defining A*(t) = aA’(t)/ax,,, we have:
+
a-(t,)
(0 s t < ?I)
A-f,- 1,
= YA’(&)
(71)
-v =
X+(t) (f&Y++_&,)
+
(tl < t d tf)
A+f,+,
A+(?,) = @,+,X+(G).
+_f,1
(72) (73)
To find the dependence of x*(i) on Y and T(O), we again start from the integral equations (29) and (30). Definingp’(t) = ax*(t)/aY and q*(t) = ax’(t)/aT(O), we have: dp-(0 dt
-=
1 r
fx-p-(t)
_!f-
(0 G t < tJ
Y
(74)
p-(O) = 0
W (0
-
dt
(75)
=.L+P+(t) +“fY+
(t1
c t c ?f)
P’(h) = w-(h) -+x-(4) -x0
(76) (77)
dq-(r) - dt
1 = &L-4-0)
+_&o,)
(0 G t c t1) (78)
q-(O) = 0
1
= 5 CL-X-(z)
(70)
(65)
from which we get the differential equation dx-(2)
= yx-(t1)
(68)
(79)
(0 s t < t*)
and the initial condition:
(%)
dq+(t) - dt
=.L+q+(t)
q+(t,) = N_(G). X-(O) = I where I is the n X n identity matrix. Similarlyforthepost-injectionsectionwehave
+.Go,
(fl < t s ff)
(80)
(81)
(67) It can be easily shown that the following equations hold for p*(t) = aA*(t)/dY and u*(f) = aA*(t)/aT(O): 286
Optimal cold-shot cooling of an adiabatic tubular reactor
dp+W dt =-~{p-fz-+A-(f;.,p--~fz-)} (0 c t < t1) P-(b)
= yP+(h)
+A+(&)
(82)
.
(83)
,.
where dQ/dx, and dZ?/dx, are total derivatives and take into account the dependence of (x*(t), A*(t) ) on X~ The various coefficients appearing in (90) and (91) are determined as follows:
dp+(t) -=-{P+f,++A+(fLp++fk,)~ dr
(84) (85)
2
da-(t) dt
(0
et*)
=
da+(r) _=_o+fz+dt
--I) + A+(t,) (x-(tl) -x,,)
= A+(&) (X-(t,) 0
-=-;{u-fz-+A-(f,zq
“-(*-f
s t < t1)
ya+ta A+(f Lq+ (t1 < t
c tf)
(86)
+
@)Y.x,fx+(ff)
-$
(87)
+A-(fz-X-+fG)}
_
I 0
df+ (T(O) -To)
dt
X
+f$no,)
(92)
038)
Now, suppose that for a given (xg*, Tz), the optimal control parameters ?*(O), q* and i: are determined and the trajectory i* (t; x,f, T,* ) is computed. Suppose that an adiabatic reactor is designed with cold-shot injection at t = I,*. If now a small perturbation, 6x,, = x0-x$, in the feed composition occurs, this will reduce the reactor performance, particularly if the state variables are very sensitive to x0. (We ignore the sensitivity calculations for To, for the new feed can always be preheated to the design level, T$. However, if the cost of this preheating operation becomes high, then the treatment should be modified.) Now, Eqs. (86)-(89) give precisely the information we want regarding the sensitivity problem. Based on this knowledge, the corrections in optimal control, 6f(O) = F(O) -f*(O), and ST = 9 - T* are determined by requiring that the vector (6x,, 6T(O), 6~) lie on the intersection of the tangent hyperpianes to the surfaces & = 0, Z? = 0 This requirement in quantitative terms becomes:
)l
+g;,zodt Qy =
287
P+(&)(x-@I)
--x0) +A+(h)p-(1,)
+%,r+@P,,,~p+(G)
-$
I
:-{p-f
(93)
+$
- + A-f,-p-}
I
Gi A-f-
dt
0
dt
p+&++A+
x
de+
(
UP++&
p-g,
)I + A-
dt
(94)
L. PADMANABHAN
+
@no, +Y(%o,,v + %oLqP+(ff) 1 dgT+
If
II
+Y
(
yp+g,+ + 2h+ ~Pft-gG*,
rr+
)I
a (95)
Qno, = (a,,, no, +%,xf4+(~fH+~+(~l> x W(h) -$
J
+
+ h+(4)q-(h)
The inverse used in (99) may not exist in some cases. For example when the surfaces Q = 0, d = 0 are cotangential, or nearly cotangential, at (x$, ?* (0)) +*), the system (90), (91) has no solution. For such cases, the linear approximations used in (90) and (91) are not valid, and one has to employ a second-order expansion of_ Q and R, approximating the surfaces Q = R = 0 by quadric surfaces.
;-{v-f-+h-(fX-q-+f~oj)}dt CONCLUSION
h+g,+ + (T(0) - To) cr+g,+
+A
G(O)=
--x0)
and S. G. BANKOFF
@cl-+
d (~4++&to, >I1
Y (@TtOLT(O)
(96)
+%70,,x,4+(tr))
+ de+
We close our discussion by noting that the theory and the algorithm developed here can be easily generalized for multi-bed adiabatic reactor systems. It is clear that in such cases the engineer should first weigh the possible improvement in yield that comes from the introduction of additional cold-shot injections, against the significant computational effort required. It is here that the computation of the yield for the temperatureprogrammed reactor and the simple adiabatic reactor becomes very helpful.
a+g,+ + A+ xq+
Acknowledgment-The
+
g~mo, df.
(97)
authors thank the Northwestern Computation Center for a grant of computer time and the National Science Foundation for financial support provided under Grant NSF GK-1126.
All these quantities given by Eqs. (92)-(97) are evaluated on the predetermined optimal path to the control L?*(t; x$, To) corresponding (T*(O),+*). Using these in (90) and (91) we finally arrive at the desired feedforward optimal controller: (,ro))
= (&;o))+K*Sxo
NOTATION
Cj concentration of speciesj CjO inlet concentration of speciesj f rate functions defined in (10) and (11) g rate functions defined in (1) and (2) k iteration number P sensitivity vector defined in (74)-(77) 4 sensitivity vector defined in (78)-(81) t residence time cold shot injection point t1 total residence time for the reactor tf x dimensionless concentration vector dimensionless initial concentration vector x0 A a (3 x 3) matrix employed in (4 1) Fo volumetric feed rate I identity matrix J adiabatic temperature rise factor used in (3) and (4) K* gain matrix for the linearized controller
(98)
where K* is the proportional-gain matrix associated with the controller and is given explicitly as:
Note that dQ */duo and & */dx, are n-dimensional row vectors. 288
Optimal cold-shot cooling of an adiabatic tubular reactor
Lagrangian defined in (12) quantity defined in (22) quantity defined in (23) quantity defined in (2 1) local temperature raw-feed temperature preheat temperature quantity defined in (60) sensitivity matrix generated by (66)-(69). time-scaling factors defined by (24) through-feed fraction step-size factor dummy variable functions defined in (24)
aL R S T
TO T(O)
u
p: Y
z
7’
n
reduced-time coordinate injection point in reduced-time domain total residence time in reduced-time domain adjoint vector control vector incremental vector sensitivity vector obtained from (82)-(85) sensitivity vector obtained from (86)-(89) objective function vector defined by (43)-(45) sensitivity matrix associated with the adjoint vector vector with components S, Q and R
REFERENCES 111 PADMANABHAN L., and BANKOFF S. G.,J. Optimiz Theory Applic. Accepted for publication. [2] DYSON D. C. and HORN F. J. M.,J. Optimiz. Theory Applic. 1967 140. [3] ARIS R., Optimal Design of Chemical Reactors. Academic Press 196 1. [4] DYSON D. C. et al., Can. J. Chem. Engng 1967 45 310. [5] BELLMAN R., Dynamic Programming. Princeton University Press 1957. [6] PADMANABHAN L. and BANKOFF S. G., Znt.J. Control. In press. R&arm& Un modble optimal de reactems adiabatiques a &ages multiples, avec refroidissement par jet froid, est considbt dans le cas de plusieurs reactions chimiques independantes. Les equations necessaires pour I’optimalitC sont obtenues en appliquant les methodes de Pontryagin au probleme de variation. Un programme de calcul base sur la methode de Newton dans la zone controle-parambtre est &labor& et applique au systeme de reactions consecutives A -+ B + C dans un reacteur tubulaire avec une injection simple de Boid. Les limites superieures et inferieures de l’index de performance de teis reacteurs sont calculees par l’analyse des configurations limitatives. Une analyse de la sensibilite conduit a un programme d’elaboration dun controleur optimal d’avance qui compense les faibles perturbations rencontrees dans les conditions initiales. ZusammenfassnngDie optimale Auslegung eines adiabatischen Mehrbettreaktors mit Kalteinspritzung wird fiir den Fall mehrerer, voneinander unabhtigiger chemischer Reaktionen betrachtet. Die fiir die Bestimmung des Optimums erforderlichen Gleichungen werden durch Anwendung der Methoden von Pontryagin auf das Variationsproblem abgeleitet. Ein Berechnungsschema, das sich auf die Methode von Newton im Raume de Kontrollparameter stiitzt, wird ausgearbeitet und auf das Folgereaktionssystem A + B + C in einem Rohrenreaktor mit einer einzigen Kalteinspritzung angewendet. Die unteren und oberen Grenzen der Leistungszahl solcher Reaktoren werden durch eine Analyse der Grenzanordnungen berechnet. Durch Empfindlichkeitsanalyse wird ein Schema zur Entwicklung eines optimalen, direkt wirkenden Reglers erhalten, der geringfiigige Storungen in den Anfangsbedingungen ausgleicht.
289 C.ES.Vd.24?4a2-F
Chemical Engineering Science, 1969, Vol. 24, pp. 279-289.
Printed in Great Britain.
Optimal cold-shot cooling of an adiabatic tubular reactor L. PADMANABHAN
and S. G. BANKOFF
Department of Chemical Engineering, Northwestern University, Evanston, Illinois, U.S.A. (First received 2 August 1968; in revisedform 2OAugust 1968) Abstract - Optimal design of multi-bed adiabatic reactors with cold-shot cooling is considered for the case of several independent chemical reactions. The necessary equations for optimality are derived by applying Pontryagin’s methods to the variational problem. A computational scheme based on Newton’s method in the control-parameter space is devised and applied to the consecutive reaction system A + B 4 C in a tubular reactor with a single cold-shot injection. Upper and lower bounds in the performance index of such reactors are computed by analyzing limiting configurations. A sensitivity analvsis leads to a scheme for devising an optimal feedforward controller that compensates for small perturbations in initial conditions. INTRODUCTION
THIS paper focuses its attention on a class of unconstrained optimal control problems associated with systems exhibiting discontinuous dynamics. By this we mean that the direction field associated with the differential equations of motion suffers a finite number of discontinuities and the jumps at such points are given by continuous functions of the state variables. Such discontinuities divide the composite system into a chain of staged subsystems, each subsystem having its own characteristic equation of motion. Frequently, such d&continuities are deliberately introduced into the system in order to provide feasible approximations to some complex, perhaps unattainable, optimal strategy. Analysis of such problems is of immense value to the chemical engineer concerned with optimal design of complex plants. Confronted with the prohibitive costs of implementing some optimal control configurations on the one hand, and the unacceptably poor performance of extremely simple non-optimal or sub-optimal modes of operation on the other, the designer is compelled to seek a compromise which may lead to discrete approximations to an optimal control function. A case in point is the task of designing a tubular reactor so as to maximise the outlet concentration of a valuable product. The kinetics associated with this problem may involve several exothermic reactions.
The optimal policy is to maintain a temperature profile that favors the desired reactions leading to the valuable product and suppresses the dissipative reactions. One can easily compute this optimal temperature program by techniques which are now welldeveloped[ 11. Theoretically attractive as this may be, it is beset with the difficulty of implementing such a temperature control along the reactor. Equally difficult to achieve in practice is the idealized continuous cross-feed reactor studied in detail by Horn and others[2]. However, the yields from these reactors represent upper bounds for the performance index of a variety of other reactor configurations. On the other extreme among the optimal reactors is the adiabatic tubular reactor with the feed temperature maintained at an optimal value. This system is the simplest in structure and the easiest to accomplish, but may have a poor performance, mainly due to the fact that the temperature rise near the exit in an adiabatic reactor tends to favor the side reactions, which were successfully suppressed near the entrance. The logical answer to this is the compromise afforded by the multi-bed adiabatic reactor chain with intermittent cooling[3, 41. This cooling can be accomplished indirectly through a heat exchanger, or directly via coldshot injection, wherein a fraction of relatively cold feed stock by-passes the first few beds. The method of cold-shot cooling dispenses with
279
L. PADMANABHAN
and S. G. BANKOFF
the intermediate heat exchanger equipment and therefore lowers the capital and operating costs. This makes it attractive despite the unavoidable dilution of the enriched reaction mixture. The possibility of achieving an optimal design of the cold-shot injection and the by-pass system, so as to diminish the dilution effects, is a real incentive for the study of such reactors. In this paper we attempt to answer the problems associated with this. For the case of a single reaction, the optimal design of multi-bed adiabatic reactors has been examined in depth by Aris[3] who uses dynamic programming to evolve elegant graphical procedures that provide rapid solutions to the optimization problem. Such techniques are unique to the case of a single reaction, where temperature and conversion are the only state variables, and therefore can be depicted on a two-dimensional state plane. The effectiveness of such geometrical construction derives, in part, from the fact that, in order to maximize the yield from a single-reaction system, one simply has to maximize the reaction rate everywhere. This enormous simplification is not available when there are several competing or consecutive reactions. Although the dynamic programming concept is still applicable, its embedding approach to such problems becomes computationally inefficient and expensive. The reason for this is rather dramatically termed by Bellman the “curse of dimensionality” [5]. For the same reason the use of graphical methods is precluded. One is therefore compelled to seek analytical methods which will be computationally effective with more complex systems. One approach is to employ Newton’s first-variational scheme on the variational equations derived from the application of Pontryagin’s formulation. This will be elaborated in detail below. STATEMENT
OF
THE
PROBLEM
Let us consider a plug-flow adiabatic reactor of specified length tf, with a cold-shot injection at t, < tp The treatment for the case when tf is free can be performed easily by the techniques employed here. The raw material has composition
Cj,,, j = 1, 2, m (with m species) and temperature T,,. The total feed to the reactor is F,, (in volumetric rate units). The reactions are accompanied by evolution of heat, and the consequent temperature rise is to be curbed by ideally locating the cold-shot injection point, tl. The injection point divides the reactor into two regions: the hot pre-injection zone (0 d t < tJ and the relatively cold post-injection zone (tl < t s tf), each governed by different state equations. The discontinuity at tl is determined by appropriate material and energy balances. [The axial distance is specified in terms of the residence time t, based on the total feed F,]. Since the injection point is unspecified, the overall problem breaks down into two free-time control problems for the two sections of the reactor. In an earlier paper[6] we have established effective methods for handling free-time problems using the concept of “reduced-time domain analysis”, which provides’a method of reducing the free-time problem into a fixed-time problem. We will use this concept in what follows. Returning to the reactor design under study, the optimization is to be introduced via parametric control. We use this term to emphasize the fact that the decision vector consists of the following three parameters: the inlet temperature (or the preheat temperature), T(O), the bypass ratio, ( 1 - y) , and the location of the cold-shot injection point, tl. The preheating of the feed from To to T(0) enhances production rates in the pre-injection section. We have a genuine optimization problem in the choice of the bypass ratio, as it should be large enough to quench the hot reaction mixture, but small enough to avoid excessive dilution. We now formulate the problem for the case of IZ independent reactions. It is convenient to express the local concentrations Cj in dimensionless form by referring them to the feed and writing xj = CJ 3: Cjo, SO that T xi, = 1. If we pick II of these xj’s as a basis, then we have a state vector xeEn. In addition the local temperature T is a variable. In order to distinguish between the variables at the pre-injection and the postinjection zones; we will refer to the former by
280
Optimal cold-shot cooling of an adiabatic tubular reactor
attaching by a (+) From following
dr dt=
a (-) sign to them and to the latter sign. a local material balance we get the state equations for the two zones:
ys=
1 y&(t),
“:Xx-,
T(0) +J(x-(t) --x0)1 T(O),y)
^;f-
(10)
T-)
g(x-,
(t < t,>
(1) +(I-_y)To
F=g(x+(t),yT(O)
5 =g(x+,
T+)
(2)
(t > t1)
+J(x+(t)
’
-
T-)
(t ==ct1)
(3)
dT+ = Jg(X+, T+) dt
(t > t,)
(4)
dt
= Jg(x-,
where J E En has components Ji = -AHJC, and the inner product is implied. Here (-AHi) is the heat generation accompanying unit increase in Xi, and C, is the volumetric heat capacity of the reaction mixture. For simplicity we shall take J to be a constant vector. It then follows from Eqs. (l)-(4) that T-(t)
= T(0) +J(x-(t)
T+(t)
= T+(t,)
+J(x+(t)
-x0)
= rx-(t1)
L-(x+,
T(O), y) L f+ (11)
We will restrict y to the open interval (0, 1). The cases y= 0, 1 refer to the simple adiabatic reactor to be considered later. If the objective function to be maximized is @(xf, T(O), y), wherex,=x($) and@ E C2(En, El, El), then the necessary equations for optimality are derived in the following manner. Define the Lagrangian L by L=Q(x,,
T(O),y)+l”-
h-(t)(F-g}dt
+I,”
“+(t){f+-F}
dt
(12)
(5)
-x+(tl))
where x+ (tl) and T+ (tI) are determined jump balances at t = tl: X’(h)
-xc,))
(t1 < t c tf).
where g E C2 (x; T): En+’ --f En. Similarly we get from an energy balance: dT-
(0 c t < t1)
where the adjoint variables AkeEn form the Lagrange multipliers for the differential constraints (10) and (11). Following the standard procedure, one takes variations of all the terms in (12) in order to get the first variation of L.
(6)
through
+ (1 -Yh
(7)
T+(t,) = yT-(t,)+(l-y)T,.
(8) +
Gx-dt
Here x0 is the vector of feed concentration and To is the temperature of the raw material. The latter are simply mixing equations for the injection of the bypass stream (I- y)Fo into the through feed yFo. Using (7) and (8) in (6), we get:
- h-(t,)
ifix-
Gx+dt - A+(t f )8x+(t f )
A-f-
T+(t)
= yT(0)
+ (l-y)T,+J(x+(t)
-x,,).
(9)
Now Eqs. (5) and (9) enable one to eliminate the temperature variables from (1) and (2).
-~Sy+~g~-6T(O)
+
+
l,;
(A+g,+(T(O) - To)6y+yg,+A+6T(0)}dt. (13)
281
L. PADMANABHAN
and
It follows from (7) that A-r+&) = rh-0,)
ordinate 8 and time-scaling factors p’ defined by: + (x-(G-
(14) /
-My
where Ax*&) is the total variation in x*(tJ made up of 6x*( tl) which is the local variation in x* at fixed t1 and of (dr’ldt) I& , the contribution due to a variation in t,: h’(M
A 6x+(1,) +f(X+(tl),
Ax-(t,) L 8x-(1,)++x-(r,),
T(O),y)
(19)
A-01) = rA’(t1)
(20)
(25)
tr= (p--p+)e,+p+e,.
(26)
and
(27)
so that (26) implies that p-+/3+
= 2.
(28)
Also from 0 < c1 < tr, it follows that 0 < p’ c 2. We now recast the various equations in the reduced-time domain. First, we write the state Eqs. (10) and (11) as integral equations:
x-(p-e)
=x,+1
y op-efw(t~, I
x+(q+(e))=yx-(p-el)+
wv9
7) d5
(29)
(1-Y)-G T(O), y) de.
(30)
(21) Noting that aq+ (0)/ a0 = p+, we obtain from (29) and (30)
Q 2 A+(?,) (x-(tl) -x0) +@,-I y 2 :- A-f- dt J ’ A+gT+dt=O I t,+
t1 = p-e,
+ j;:;‘f(x+(&
= 0
(24)
so that
-f(x-(G,
T(O),y))
0 > tl) 1
8, = t,/2 = ef/2
(18)
A’(b) = d$
+(T(O)-Z-,,)
71-(d) = P-R (I < t1) [ 71+(e) = (p--p+) e,+B+e,
If we take 8, to be fixed, then the variability of t1 is reflected totally in that of 0’ through (24) and (25). e1 and 0, can be arbitrary chosen, and a natural and convenient choice is
St,. (16)
(rl -c t =sf,)
T(O),y)
t=
(15)
T(O), r> at,
In the absence of any inequality constraints, maximization of @ implies that SL be zero, which in turn requires that all the independent variations in (13)‘vanish individually. Using (14), (15) and (16) in (13) and assembling the various vanishing terms, we arrive at the necessary conditions:
s A A+(&) {f(X’(G,
S. G. BANKOFF
ax-=ae
(22)
p-
_
YfT
R 4 yfi’m,, + J;‘-A-gT-df + y* J,;+ A+gT+dt = 0.
(e < 0,)
(31)
(0 > 0,).
(32)
(23) (21), (22), and (23) are the determining tions for tl, y and T(0) respectively. COMPUTATIONAL
equa-
SCHEME
Along the same lines, the adjoint equations (17) and (18) can be written as: A-(P-0)
The computational method is easily formulated through the introduction of a reduced-time co282
= rA+(P-0,)
Optimalcold-shotcooling of an adiabatic tubular reactor sEc= (b-p,
T(O),Y) d5
(34)
(0 < 01)
(35)
(6 ’ h).
(36)
+r,
f(o) - T(o))T
(42)
and JI = (&, &. I/J~)~with the elements given by:
whence we have
ax-= ae
+I= S+AA+(e1)(f+--f-)It=tl
-5
A-f,-
!!$ = - p+
A+fz+
+A+(t,)(Yf,+--f,-)l,=t,~-(e,) $2
=
Q+W’,z, %+AA+(&)(x-(h) +A+(h)
We will also need the reduced-time version of the determining equations (21), (22) and (23): S = A+(p-0,) {f(x+(P-e,), -f(x-(P-e,), Q = h+(p-8,)
--d
h-w {A-f,- hx-(6’) +AA-(e).L-}
de
A+yAx+(e) +(T(0)-To)P+ly ox+ i
T(O), Y)
T(O), Y)] = 0
(x-(p-el)
(43)
(37) +AA+(@g,+
-x0> -$/“-
I
de
(44) e,-
$3 = R + Y@&o,,~,&+
xh-f-dB+Q,+(T(O)--T.)p+y;+
P-
I
0
x A+gT+de = 0
X AA-(@ g,-+A-d$Lr-(8))
(38)
R = /3- J;‘- A-gT- de+ y2/I+ JO; A+g,+ de + YJ%YO~ = 0.
de
AA+(e) gT++Atd$&+(8)) (39)
The proposed computational scheme is essentially Newton’s method applied to S = Q = R = 0 in the (p-, y. T(0)) parameter space. At any point CL= (p-, Q, ?(0))T, we have the following first-order Taylor expansion for 6= (5, $9 2). fi=
R+n,s/_&
(40)
where 0, = gradJl and 6, = fi - CL. But from (37), (38) and (39), we have 6 = 0, so that a first-order estimate of the location of the optimum fi can be obtained by solving (40). To this end, we take variations in S, Q, R, substitute these in (40), after making use of (14) and finally get the simplified linear system of variational equations,
(45) where
b*(e)
.
= a-c(+y(e))
--xyq(e))
and AA’(O) = fi*(+j*(e)) -A’(r)‘(@): (Note that Ax+(e) involves a comparison of the state variable at the same reduced-time coordinate and not at the same real-time. This is the main feature in the reduced-time domain analysis, as the quantities Ax’( 0)) AA’( 0) are more readily obtained than F(t) -x’(t), A*(t) -A’(t), for the latter will demand use of interpolation schemes. This is because the time-scaling factors p’ map the variable intervals [0, r,], [tI, rr] into ones of fixed length, viz., [0, e,] , rh
ASp=J,
de
(41)
where A is a 3 X 3 matrix to be defined below,
e,i).
The matrix A has elements listed below: All=0
283
(46)
L. PADMANABHAN AI2 = -A+(b)
t.fr+(x-(t)
and S. G. BANKOFF
- xo) + V(O) - To)g,+I It=t,
(47) A13 =Y--h+(h) {WTf - gT-lIt=t* A,, = $
(48)
“- A-f- de+ ( T(0) - To) I0 6
X
h+g,+
de
(49)
(35), (36) backward in the reduced-time domain, using the end-condition (19) and the jump condition (20). 3. Compute P, Qck’ and Rck’. If these are smaller in absolute value than some preassigned tolerances, the computation is terminated. Otherwise proceed to step 4. 4. If k = 0, set
I e,+
A 22 =__2pY3
e’-A-f-dt9-_rB+(T(0)-Z-o)2 I0 eI X A+& df3- G’,, I
P
I
A+ [gT++Y(T(O)
- To)gbl
(50)
yuc)
(56)
@-@mo)
e,+ (51)
AB1 = - J,“’ A-g,- de + y2 Je; A+g,+ de AX = - (@‘no, + vhom
) -?‘p+
je.B
=
-
Y%o~~o~
-
p-
Joe+-
P (k+l)
s
e’
8,
Once these quantities are yields the correction in control: 8/L= EA-‘JI
Ag& de. evaluated
,Jk’+
E(k) {A(k)}-‘q,‘k’.
(57)
If the matrix Atk’ is ill-conditioned, then use (56). 5. Increase k by unity and go to step 2. THE CASE
(53)
WHEN
y = 0,l
This corresponds to a simple adiabatic reactor without any cold-shot cooling. The only freedom we have is in the choice of the preheat temperature, T(0). For this case, the state equation is simply
A-g, de - y”/l+ X
=
(52)
A+(2&+
+ y(T(0) - To)gzG} de A33
p(k)+
where P, i = 1,2,3 are small arbitrary increments in /3-, y, T(O). If k a 1, set x = tik), A = Ack), Ax’(e) = x”(~) (nfi”c’(0)) -x’(k-i)(rl+(k-l)(e)), and use a similar result for AA’(e). Compute the vector $‘k’ and the matrix A(“). Choose a suitable l E (0, 1) and obtain the correction from:
A,, = P- ‘I- A-g,- de Ti+ Y2I 0 0,
X
(k+l) =
(54) (41)
$ = g(x, T(0) +J(x-x0)) (55)
where E, 0 < E < 1, is introduced to cut down undesirable overshoots in correction. It is now clear that (55) can be employed to develop a suitable iteration process that converges to the optimum. The algorithm is outlined below: 1. Set k = 0. Choose p-(O), y(O)and PO’(O). 2. Using x-(O) =x0 as the initial condition, integrate the state equation (3 1) forward over the pre-injection section. Then, using the jump condition (7), integrate (3 1) until the exit 8 = 0, is reached. Now, integrate the adjoint equations
(58)
with the initial condition x(O) = x,,. The adjoint equation is easily obtained as: (59) with A(tf) = CD+. The determining T(0) can be shown to be:
equation
U : aRo, + ,f Ag,,, dt = 0..
for
(60)
The computational algorithm for finding T(0) is obtained, as before, by taking variations in U and setting 0 = 0, whereupon we have:
284
Optimalcold-shotcoolingof an adiabatictubularreactor
AN
ILLUSTRATIVE
EXAMPLES
The theory developed so far will now be applied to the study of optimization for the consecutive reaction system A + B --, C conducted in a tubular reactor. The objective here is to maximize
(62) kz(T)xZ
(63)
where k,(T) = kN exp (--E,IT), j = 1,2. Using (5) and (9), we eliminate the temperature dependence from (62) and (63) to obtain c&x-, T(O) +1(x--x,)1, f(x)={,(x+,yT(0)+(l--y)To
(t < tl) (64) (t ’
+Jb+--x,)1,
h).
For numerical work, the various constants appearing in Eqs. (62), (63) and (64) are chosen as follows: k,, = 0.5 x 108,
k,, = 0.5 x 1016
E, = 6900 (OK)-‘,
Ez = 13800 (OK)-’
J, =-32
J2 = 10 (“K)
x 10= 0.95,
(OK),
X 20 =
0.05.
A good choice of the nominal control policy is essential for convergence of the iterations. The nominal control used in the example was: p-(O) = 1.75, ~‘0’= 0.9, T(O)(O)= 350°K.
Starting from these values, the control vector progressively converged in 40 iterations to the optimum: fl- = O-8575, p = O-9151, f(O) = 34OeO”K. The objective function @ at this point had the value $ ( tf) = O-5176, as compared with xz(tf) = 0.2852 corresponding to the nominal policy. A fourth-order Runge-Kutta method was used to integrate the differential equations, partitioning the reactor length tf = 10 into 40 equal intervals. The various integrals appearing in Eqs. (43) to (54) were evaluated using Simpson’s rule, with a partition of the reactor length into 80 equal intervals. The computation, performed on a CDC6400 computer, took 0.57 set of central processor time per iteration. Also for the same reaction system, the optimal design of the simple adiabatic reactor without cold-shot cooling was considered. The computational procedure using Eq. (61) was employed to determine the optimal preheat temperature. Convergence was observed in 20 iterations to the optimal inlet temperature T(0) = 338”K, and the optimal yield was x,(tf) = O-5134, which, as expected, is less than the yield from the reactor with cold-shot cooling. The computation time was 0.43 set per iteration. The yield of the reactor with cold-shot cooling can be improved further by increasing the number of the injection points, which give more degrees of freedom. As the number of injection points increases, the yield of the reactor increases monotonically and ultimately attains an upper bound in performance index corresponding to the idealized continuous cross-feed reactor. However, the computational problems associated with this limiting configuration are not easily tractable, particularly in multi-reaction systems. This is due to the possibility of singular control. A more easily computed upper bound is the yield from the optimal temperature-programmed reactor. For the consecutive reaction system under study, the optimal yield from the programmed reactor was found to be x, (t,) = 0.5452. By comparing this figure with the yield from the simple adiabatic reactor, one can get a rough estimate of the improvement in yield to be expected from more
285
L. PADMANABHAN
and S. G. BANKOFF
cold-shot injection points. This improvement will be significant for those reaction systems that are strongly exothermic. LINEARIZED
(tl < t s tr) x+(t,)
FEEDFORWARD
CONTROLLER
In the course of the computation process for the reactor with cold-shot cooling, it was found that the yield was extremely sensitive to variations in the preheat temperature, bypass ratio, and feed composition. This makes it imperative to incorporate into the reactor system a feedforward controller that continually responds to perturbations in the initial conditions by providing compensating changes in the control parameters, thus restoring the system to an optimal configuration. Our approach to this problem is to approximate the surfaces Q = 0, R = 0 by their tangent planes at the point (x$, f * (0)) Y* ) . The asterisk symbol refers to the fact that the optimum (Z?(0)) 9* ) , corresponding to a particular choice of x$, is predetermined. In the neighbourhood of this point, the tangent planes are valid approximations to the optimal response surfaces. The determination of such tangent planes involve the calculation of the sensitivity of Q and R ‘with respect to x0, T(0) and y. For this, we go to the integral equations (29) and (30) and obtain by differentiation with respect to the initial condition, the sensitivity matrix X*(t) k (ax’(t)/(ax,):
X-(t) =
I+;,‘,(.&-X-(5) +.&Id5
-
dt
+ (1 --y)Z.
(69)
Analogous results for the sensitivity of the adjoint variable A’ (t) are obtained from (33) and (34). Defining A*(t) = aA’(t)/ax,,, we have:
+
a-(t,)
(0 s t < ?I)
A-f,- 1,
= YA’(&)
(71)
-v =
X+(t) (f&Y++_&,)
+
(tl < t d tf)
A+f,+,
A+(?,) = @,+,X+(G).
+_f,1
(72) (73)
To find the dependence of x*(i) on Y and T(O), we again start from the integral equations (29) and (30). Definingp’(t) = ax*(t)/aY and q*(t) = ax’(t)/aT(O), we have: dp-(0 dt
-=
1 r
fx-p-(t)
_!f-
(0 G t < tJ
Y
(74)
p-(O) = 0
W (0
-
dt
(75)
=.L+P+(t) +“fY+
(t1
c t c ?f)
P’(h) = w-(h) -+x-(4) -x0
(76) (77)
dq-(r) - dt
1 = &L-4-0)
+_&o,)
(0 G t c t1) (78)
q-(O) = 0
1
= 5 CL-X-(z)
(70)
(65)
from which we get the differential equation dx-(2)
= yx-(t1)
(68)
(79)
(0 s t < t*)
and the initial condition:
(%)
dq+(t) - dt
=.L+q+(t)
q+(t,) = N_(G). X-(O) = I where I is the n X n identity matrix. Similarlyforthepost-injectionsectionwehave
+.Go,
(fl < t s ff)
(80)
(81)
(67) It can be easily shown that the following equations hold for p*(t) = aA*(t)/dY and u*(f) = aA*(t)/aT(O): 286
Optimal cold-shot cooling of an adiabatic tubular reactor
dp+W dt =-~{p-fz-+A-(f;.,p--~fz-)} (0 c t < t1) P-(b)
= yP+(h)
+A+(&)
(82)
.
(83)
,.
where dQ/dx, and dZ?/dx, are total derivatives and take into account the dependence of (x*(t), A*(t) ) on X~ The various coefficients appearing in (90) and (91) are determined as follows:
dp+(t) -=-{P+f,++A+(fLp++fk,)~ dr
(84) (85)
2
da-(t) dt
(0
et*)
=
da+(r) _=_o+fz+dt
--I) + A+(t,) (x-(tl) -x,,)
= A+(&) (X-(t,) 0
-=-;{u-fz-+A-(f,zq
“-(*-f
s t < t1)
ya+ta A+(f Lq+ (t1 < t
c tf)
(86)
+
@)Y.x,fx+(ff)
-$
(87)
+A-(fz-X-+fG)}
_
I 0
df+ (T(O) -To)
dt
X
+f$no,)
(92)
038)
Now, suppose that for a given (xg*, Tz), the optimal control parameters ?*(O), q* and i: are determined and the trajectory i* (t; x,f, T,* ) is computed. Suppose that an adiabatic reactor is designed with cold-shot injection at t = I,*. If now a small perturbation, 6x,, = x0-x$, in the feed composition occurs, this will reduce the reactor performance, particularly if the state variables are very sensitive to x0. (We ignore the sensitivity calculations for To, for the new feed can always be preheated to the design level, T$. However, if the cost of this preheating operation becomes high, then the treatment should be modified.) Now, Eqs. (86)-(89) give precisely the information we want regarding the sensitivity problem. Based on this knowledge, the corrections in optimal control, 6f(O) = F(O) -f*(O), and ST = 9 - T* are determined by requiring that the vector (6x,, 6T(O), 6~) lie on the intersection of the tangent hyperpianes to the surfaces & = 0, Z? = 0 This requirement in quantitative terms becomes:
)l
+g;,zodt Qy =
287
P+(&)(x-@I)
--x0) +A+(h)p-(1,)
+%,r+@P,,,~p+(G)
-$
I
:-{p-f
(93)
+$
- + A-f,-p-}
I
Gi A-f-
dt
0
dt
p+&++A+
x
de+
(
UP++&
p-g,
)I + A-
dt
(94)
L. PADMANABHAN
+
@no, +Y(%o,,v + %oLqP+(ff) 1 dgT+
If
II
+Y
(
yp+g,+ + 2h+ ~Pft-gG*,
rr+
)I
a (95)
Qno, = (a,,, no, +%,xf4+(~fH+~+(~l> x W(h) -$
J
+
+ h+(4)q-(h)
The inverse used in (99) may not exist in some cases. For example when the surfaces Q = 0, d = 0 are cotangential, or nearly cotangential, at (x$, ?* (0)) +*), the system (90), (91) has no solution. For such cases, the linear approximations used in (90) and (91) are not valid, and one has to employ a second-order expansion of_ Q and R, approximating the surfaces Q = R = 0 by quadric surfaces.
;-{v-f-+h-(fX-q-+f~oj)}dt CONCLUSION
h+g,+ + (T(0) - To) cr+g,+
+A
G(O)=
--x0)
and S. G. BANKOFF
@cl-+
d (~4++&to, >I1
Y (@TtOLT(O)
(96)
+%70,,x,4+(tr))
+ de+
We close our discussion by noting that the theory and the algorithm developed here can be easily generalized for multi-bed adiabatic reactor systems. It is clear that in such cases the engineer should first weigh the possible improvement in yield that comes from the introduction of additional cold-shot injections, against the significant computational effort required. It is here that the computation of the yield for the temperatureprogrammed reactor and the simple adiabatic reactor becomes very helpful.
a+g,+ + A+ xq+
Acknowledgment-The
+
g~mo, df.
(97)
authors thank the Northwestern Computation Center for a grant of computer time and the National Science Foundation for financial support provided under Grant NSF GK-1126.
All these quantities given by Eqs. (92)-(97) are evaluated on the predetermined optimal path to the control L?*(t; x$, To) corresponding (T*(O),+*). Using these in (90) and (91) we finally arrive at the desired feedforward optimal controller: (,ro))
= (&;o))+K*Sxo
NOTATION
Cj concentration of speciesj CjO inlet concentration of speciesj f rate functions defined in (10) and (11) g rate functions defined in (1) and (2) k iteration number P sensitivity vector defined in (74)-(77) 4 sensitivity vector defined in (78)-(81) t residence time cold shot injection point t1 total residence time for the reactor tf x dimensionless concentration vector dimensionless initial concentration vector x0 A a (3 x 3) matrix employed in (4 1) Fo volumetric feed rate I identity matrix J adiabatic temperature rise factor used in (3) and (4) K* gain matrix for the linearized controller
(98)
where K* is the proportional-gain matrix associated with the controller and is given explicitly as:
Note that dQ */duo and & */dx, are n-dimensional row vectors. 288
Optimal cold-shot cooling of an adiabatic tubular reactor
Lagrangian defined in (12) quantity defined in (22) quantity defined in (23) quantity defined in (2 1) local temperature raw-feed temperature preheat temperature quantity defined in (60) sensitivity matrix generated by (66)-(69). time-scaling factors defined by (24) through-feed fraction step-size factor dummy variable functions defined in (24)
aL R S T
TO T(O)
u
p: Y
z
7’
n
reduced-time coordinate injection point in reduced-time domain total residence time in reduced-time domain adjoint vector control vector incremental vector sensitivity vector obtained from (82)-(85) sensitivity vector obtained from (86)-(89) objective function vector defined by (43)-(45) sensitivity matrix associated with the adjoint vector vector with components S, Q and R
REFERENCES 111 PADMANABHAN L., and BANKOFF S. G.,J. Optimiz Theory Applic. Accepted for publication. [2] DYSON D. C. and HORN F. J. M.,J. Optimiz. Theory Applic. 1967 140. [3] ARIS R., Optimal Design of Chemical Reactors. Academic Press 196 1. [4] DYSON D. C. et al., Can. J. Chem. Engng 1967 45 310. [5] BELLMAN R., Dynamic Programming. Princeton University Press 1957. [6] PADMANABHAN L. and BANKOFF S. G., Znt.J. Control. In press. R&arm& Un modble optimal de reactems adiabatiques a &ages multiples, avec refroidissement par jet froid, est considbt dans le cas de plusieurs reactions chimiques independantes. Les equations necessaires pour I’optimalitC sont obtenues en appliquant les methodes de Pontryagin au probleme de variation. Un programme de calcul base sur la methode de Newton dans la zone controle-parambtre est &labor& et applique au systeme de reactions consecutives A -+ B + C dans un reacteur tubulaire avec une injection simple de Boid. Les limites superieures et inferieures de l’index de performance de teis reacteurs sont calculees par l’analyse des configurations limitatives. Une analyse de la sensibilite conduit a un programme d’elaboration dun controleur optimal d’avance qui compense les faibles perturbations rencontrees dans les conditions initiales. ZusammenfassnngDie optimale Auslegung eines adiabatischen Mehrbettreaktors mit Kalteinspritzung wird fiir den Fall mehrerer, voneinander unabhtigiger chemischer Reaktionen betrachtet. Die fiir die Bestimmung des Optimums erforderlichen Gleichungen werden durch Anwendung der Methoden von Pontryagin auf das Variationsproblem abgeleitet. Ein Berechnungsschema, das sich auf die Methode von Newton im Raume de Kontrollparameter stiitzt, wird ausgearbeitet und auf das Folgereaktionssystem A + B + C in einem Rohrenreaktor mit einer einzigen Kalteinspritzung angewendet. Die unteren und oberen Grenzen der Leistungszahl solcher Reaktoren werden durch eine Analyse der Grenzanordnungen berechnet. Durch Empfindlichkeitsanalyse wird ein Schema zur Entwicklung eines optimalen, direkt wirkenden Reglers erhalten, der geringfiigige Storungen in den Anfangsbedingungen ausgleicht.
289 C.ES.Vd.24?4a2-F