W’ kV,
= oil flow rate, Ib.;hr. = coolant flow rate. lb.,’hr.
Literature Cited
(1) Bollinger, R. E., Lamb, D. E., IXD. ENG. CHEM.FUNDAM E N T A L S 1, 245-52 (1962). (2) Haskins. D. E., ”Synthesis of Invariance Principle Control Systems for Chemical Processes,” Ph.D. dissertation, University of Oklahorna. 1964. ( 3 ) Kulebakin. V. S., .‘The Theory of Invariance of Regulating and Control Systems.” “’Automatic and Remote Control,” Proceedings of the First Interriational Congress of the International Federation of Autoinatic Control, Moscow, CSSR, 1960. Vol. I. pp. 106-16. Butterworths. London, 1961. (4) Petrov, B. N.. “The In\-ariance Principle and the Conditions for Its Application during the Calculation of Linear and Nonlinear Systems.” Ibid.. pp. 117-25. (5) Stewart. IV. S.. “Dynamics of Heat Removal from a Jacketed, Agitated Vessel.” Ph.D. dissertation, University of Oklahoma, 1960.
GREEKSYMBOLS = a characteristic matrix [ 41 = density of coolant, lb./cu. ft. pc = density of oil, Ib.icu. ft. p, = density of wall metal, Ib./cu. ft. pLD SUBSCRIPTS = coolant variable = oil variable
C
f zu
= =
ss
wall variable steady-state variable
SUPERSCRIPT
*
RECEIVED for review June 26, 1964 ACCEPTED .January 15, 1965
= transient variable
GREEN’S FUNCTIONS AND OPTIMAL SYSTEMS Complex Interconnected Structures M . M. DEN”
AND
RUTHERFORD A R l S
Cnireisitj of .Iriinnesota. dlznnrapoiis, .Mznn.
The solution of first-order variational equations b y means of Green’s functions i s utilized to develop necessary conditions and computational methods for the optimization of plants with complex interconnections. A concrete example of a staged system and a system with both discrete and continuous elements i s completely solved.
has recently been focused on the opof large systems whose units are interconnected in a complex fashion. ‘The system of units linked in a straight chain or sequence has long been studied and is amenable to optimization in a variety of ways--notably. by calculus with the use of Lagrange multipliers, Pontryagin’s principle, and dynamic programming. These methods have been extended to chains with feedback and to more general systems ( 7 , 2, 779, 72--7R). I n some cases the adaptation of Pontryagin’s methods and of dynamic programming has been incorrect and the work of Jackson (76) and Jackson and Horn (7,l. 7.5) has been noteworthy in discovering both gross and subtle fallacies. We are not aware of any published calculations 011 complex systems that have completely foundered on these shoals. but this is in part due to the scant attention that has been paid to real computational feasibility. T h e only systems that have been studied in detail are t ~ v o -and three-stage co- and countercurrent autothermic reactors (27) and systems for which the policy is disjoint or equivalent to a straight-chain optimization and a one-dimensional search ( J , 72. 27). I n this paper we use the solution of the linear variational equations by Green‘s functions to obtain both the necessary conditions and computational techniques for systems of arbitrary topology containing both discrete and continuous elements. .A concrete numerical example is included to demonstrate the workings of the computational algorithm. UCH ATTENTION
M timization
Topological Structure
By a serial or “straight-chain“ process we mean a sequence of stages such as that shobvn in Figure 1. Its stages may be numbered consecutively in the direction of flow of information, the input to the nth stage consists solely of the entire effluent of the ( n - 1)st stage. and no component of the decision vector a t one stage appears in the transformation equation for any other stage. A continuous process of duration 0 whose finite difference approximations form straight chains will also be called a serial process. T h e recycle systems considered previously (7, g! 12) are examples of nonserial processes, and other basic structures have been discussed ( 7 , 2, 73, 77). Every complex system can be broken down into a number of subunits, each a serial process, together with matching conditions imposed by continuity relations. An optimal system will not, in general, be made u p of the set of optimal serial processes which satisfies the continuity relations. This fact has been dramatically demonstrated by .Jackson ( 7 , 76). and is clear from the form of the new boundary value problem caused by recycle ( 7 , 9 ) . \Ve can, however. extend all of the results which we have developed for straight-chain systems to complex rystcms of general topological structure ( 7 . 9. 70).
STAGE STATE DEClSlO N
Present address. Department of University of Delaware, Newark. Del. I
248
I&EC FUNDAMENTALS
Chemical
Engineering,
p,
- r k ,---pFll ---iFk qn
41
Figure 1 .
Serial process
qH
Variational Equations
T h e kth seiial subprocess of a complex system consists of a sequence of stages with transformation eqlrations pnk
= Tnk
(Pnk-1,
qnk);
nh
=
Ih,
211,
5
A v h
If we substitute Equation 8 into Equations 4 and 5 and sum over all subsystems, we obtain
(11
or a continuous process with differential transformation stages
p is the vector with Cartesian basis whose components are sufficient to describe the state of the system and q is a decision vector which is to be chosen in such a way as to maximize a given profit function. 'l'here is no loss of generality in assuming that the right-hand side of Equation 2 is independent of tk (7). It is convenient to suppose that the dimension, S.of p is the same throughout the system, and this ma?- always be realized Mithout loss of generality by defining additional state components bvith suitable transformations. The dimension of q may be different for each stage or subprocess. The input, pok,to the kth subprocess is of the form
where prt is the ith external input to the total system and P.,-~ is the output from thr j t h staged or continuous serial subprocess. Equations 3 are the matching conditions which must be satisfied. If we assume a particular set of decision vectors, we may solve Equations 1 to 3 . 'The first-order equations of variation about this solution may be obtained by Taylor's theorem. and Green's identity written for each serial subprocess (9. 7 9 ) . as follows :
all continuous subsystems
Indices k and rn in the second summation on the left have the same range, and we may interchange order of summation and then dummy indices to obtain
c k
[%l
-
E ~ o " m i . ~ i ~ k , &Lvk3 m)ll = rn
c ani b T,'
6q,'
----
a11 n
bqn'
cuo,i Cfii(lk)j 6Pf13 k
+
s
i
u,
$26 q j dt
+ (12)
all t
T h e notation on the last sum denotes the sum over all discrete stages, and on the last term, integrals for all continuous elements. Weak Maximum Principle
We will assume without loss of generality that the profit is always expressed in the form
Thus if we define the mixed boundary conditions for the Green's vectors as
(4) or, for continuous systems and the staged and continuous Hamiltonians by
Hn (The summation convention, in which summation is indicated over all values of an index repeated once in the upper and once in the lower part of an expression, is used throughout.) Subcripts and ok refer to the output and input of the kth staged or continuous subsystem. 6 p is the variation in the state about the solution and results from the variations 6q in the decision and 6 p o , in the input. T h e Green's vectors satisfy the equations
(7) T h e variation in pOkis
where
(9)
=
H(t)
u,, Tnz
(15)
T1
(161
=
i t follows from Equation 12 that the variation in the profit is
s ~2 bH
6q3 dt
(17)
all 1
We then obtain the weak version of the maximum principle of Pontryagin (20) by proofs identical to those of Denn and .4ris (7: 9 ) . The optimal decision vector has the property that the stage or continuous Hamiltonian is made stationary by components which lie a t the interior of the admissible set and a maximum by components at a constraint. Furthermore, the continuous Hamiltonian is constant over each continuous serial subprocess. The transversality and natural boundary conditions for the Green's vector are obtained in the same way. as follows: If the effluents lie a t equality constraints of the form Qkt{Px,])
0
(18)
then the boundary conditions for the Green's vectors are
6 p f , is a variation in the lth input to the complex system, the first summation in Equation 8 is over all inputs, and the second over all serial subsystems.
tvhere X is a vector multiplier VOL. 4
NO. 3
AUGUST
1965
249
*---e *---Nn
5; POI
Figure 2.
Diverging branch
If any feed element pfl’ is unspecified, then the Green’s vector must satisfy the additional condition
c
aoktni(6,k)j =
0
(20)
k
If the input vectors must satisfy constraints of the form Q’({pn}) = 0
bb,
where 7 is a vector miiltiplier. For systems with constraints on the allowable state space or on both the allowable decisions and states, Equations 4 and 5 and the weak maximum principle remain unchanged for independent decisions if appropriate modifications are made in the equations defining the Green’s vectors, as described (3, 6, 9). In particular, the “jump condition” described (3, 6, 20) may be looked upon as a particular application of Equation 19. Approximation to Problem
We may identify a particular Green’s vector ack)Iwith the ith component of the kth output by the boundary condition
and we may then write Equation 12 as
If we make the assumption, as in (9),that all decisions are chosen optimally and changes in decisions between neighboring optimal policies are small, Equation 24 becomes approximately =
fi3(m,i)r
1
8bfmr
(25)
m
If all effluents are free and all inputs to the system fixed, and if we assume a set of effluents and use the weak maximum principle to determine the optimal policy, then Equation 25 may be used as in ( 9 ) to correct the assumed outputs. Such a directed iteration is needed for the type of solution proposed by Jackson (77, 73) in order to avoid the “curse of dimensionality .” Certain obvious difficulties arise in this type of solution. For example, in the diverging branch of Figure 2 we must assume values for both pNIr and pNIII. If the result from system I1 is used to specify pXI,we must iterate on branch I11 as in ( 9 ) until the same value of p2vIis reached. Only then may the computation be carried back to the input and Equation 25 applied. 250
l&EC FUNDAMENTALS
’ (27)
(21)
the Green’s vector must satisfy the additional relations
SP,,~
For the recycle loop shown in Figure 3 Equation 14 reduces to
It is thus necessary to assume both pAvIand til,vI. T h e optimal path with respect to the assumed Green’s vector may then be determined along each part of the process, but in general Equation 27 will not be satisfied. If the optimal trajectory were independent of the values of the Green’s vector, Equations 6 and 7 would be linear, and a single correction would produce the proper value of uN,, as follows: Assume a.vf* and compute to Use Equation 27 to and compute to f c o r r * . T h e correct value of determine aNff*
mwr is then found from the linearity to be U N r = 2 tiw,* - U O I I *
- bP ~
bpNf In reality, however, the optimal trajectory will change with changes in axf and Equation 28 is only the expression for a first-order iteration procedure. T h e sensitivity of the optimal trajectory to changes in the initial values of the Green’s vector which we have noted (9) and which have been observed by others ( 7 7) would lead one to expect slow convergence. When convergence has been obtained, we have a n optimal trajectory for the assumed p w r which gives a value 6p f r . We may then apply Equation 25 and repeat the process. Because the optimal trajectory has been computed once and for all, the equations for the Green’s tensor defined by Equation 23 are truly linear and the correct Green’s tensor may be determined with a single adjustment. This over-all iteration may also be expected to exhibit some instabilities, however, because the initial value u ,, will now change as a result of the new pXI. The potential difficulties which we have indicated for branches and recycle are prototypes of those which might be expected in complex systems. In general, branches will lead to greatly increased calculation (“increase in dimensionality”) and feedback to potential instability in the iterative process. Approximation to the Solution
Let us assume that all inputs are either fixed or may be absorbed into a single stage as a decision variable. With the boundary condition (Equation 14) we may then write Equation 12 as
If the decision space is parametrized by a length variable, s, then dividing by 6s and letting 6s + 0 we obtain
dP
-
bT,,'
dq,J
r
bT'
we represent. all decisions in the entire complex system as elements of a vector q ? the transformation equations are of the form
dqj
\.\-e seek the direction in decision space of the greatest rate of increase of d P ds subject to the metric constraint (70) If the Green's vectors are defined by Equations 6 and 14 and all unspecified inputs are considered to be decisions, then in the usual manner it follow7s that the variation in profit is I t easily follo\vs from the Xveak maximum principle that the gradient direction a t each stage is
(35) Thus the weak maximum principle must be applied to the sum of all stage Hamiltonians. If a n element q' appears a t only one stage, this reduces to the previous case. Similarly. for the gradient method with metric constraint
4,' ~~
dj
(32)
the gradient direction is
and at each point of a continuous system
d4'-
(37)
-
ds
A Discrete Example
(33) T h e elements of the associated metric tensors g z j and gn23 are formed from the inverse matrices of the elements of gi, and grit,, respectively. 'The extension to constrained systems is carried out exactly as in (5) and (701, using the Green's tensors defined by Equation 23. For the case where all subprocesses are continuous, Equation 33 is a generalization of the result of Kelley ( 7 9 ) . T h e optimization problem is then solved by computing improved values of the decision by moving in the direction indicated by Equations 32 and 33: and repeating the process until the maximum is attained. For the diverging branch of Figure 2 the calculation of the state may be carried out in the fonvard direction, and there is no need for an iterative solution which leads to a n increase in dimensionality. For the: recycle process sho\vn in Figure 3 it is necessary to solve the state equations by some iterative method, and in general Equation 25 \vi11 be the most effective. Because the decision vectors are specified, Equation 25 is exact to first-order, and Equation 28 is exact and need be applied only once. Unlike "approximation to the problem," the solution of the recycle state equations is not subject to the instabilities inherznt in changing decision vectors, and the over-all stability of the gradient solution of the optimization problem appears more certain ( 7 7 ) . I t is for these reasons that Lve feel that the gradient method is more suitable for the solution of complex problems than the method of the previous section.
As a n example of the solution of a complex staged process we have considered the system shown in Figure 4. T h e external feed is pure A. T h e reaction A + B + C is carried out in three stirred-tank reactors of residence time 0: in which the temperature is the decision variable. This is followed by a two-stage countercurrent extraction system, with solvent addition to each stage as the decision. T h e solvent and reaction stream are immiscible and the solvent is assumed to extract only pure A. which is recycled to the feed to the first reactor. SOLVENT 7
T SOLVENT REMOVAL
EXTRACTION OF A (CHOOSE AMOUNT OF SOLVENT)
PURE A
\Ve will assume that the reaction q s t e m is the same as the one studied previousl) ( 9 , 70). \+ith kinetic equations a,-1
Interacting Decisions
\Ve have assumed until now that each decision appears explicitly in the transformation equation of only one stage. This assumption is unduly restrictive, and in particular will not be valid if Lve seek to consider a n unspecified input as a decision ivhen such an input is a feed to more than one stage. For convenience we will consider only staged systems. If
5
4
b,-:
-
bfl
-
a, -
0 kl(T,) an2 = 0 ; n
=
1. 2, 3
+ 0 k : (T,) a n a 0 k2
(T,) 6,
=
0 ; n
= 1, 2. 3
(38)
(39)
In addition. i t is convenient to define a fictitious state variable, u,. such that u,
= U,-I
VOL. 4
; n = 1, 2, 3 NO. 3
AUGUST
(40) 1965
251
We will define the ratio of solvent input a t stage n LO flow rate of the reactant stream by u,, and we will denote the equilibrium distribution between the concentrations of a, in the solvent and reactant streams by $ ( a n ) . We may then write the transformation equations a t stages 4 and 3 as u, =
+
u,-l
6 , = b,-,; a4
+ us
n
=
=
4, 5
+
-
c j $(aj)
a4 =
(42) -
as =
0
(43)
0
(44)
T h e boundary conditions (Equation 14) are
Equation 43 may be rewritten as
aq
T h e external feed to the system is pure A of concentration and the feed to stage 1 is bo
..
a3, -
a ~ ,
b'&t
-
--
-2331
Y j=
qD5
bj31 b't)3j
A
sat
-
-
bp5
dP0 -= 0 dp33
.
~
bP
~-'j
~
or in terms of components
0
= u, =
$ 'G 4 )
(41)
4, 5
- (us - u4) $ ( a s )
$(aq)
a:
n
0,;
013 =
(46) 013
A
+-^+a,
cy5
= p
(62)
U5
or using Equations 41 and 45 a, =
pj
+ a3
af
(48)
- a5
A
1'5
We seek to choose the temperatures and solvent flow rates to maximize
P
= bg
+
+
(63)
= 1
cy3
A
$
-A-
(a4) =
-'y
US cy3
- yuj
pas
0.3-
(49)
-cy,=o Ub
T h a t is, we seek to maximize the production of B while allowing for the costs of raw materials and extraction. The information flow or functional diagram (2) for the system defined by Equations 38 to 42. 44. 46, and 48 is sho\vn in Figure 5. ?'his diagram differs from the flolv diagram of
P3
-
P3
=
0
u3
-
;3
=
0
T h e decision variables are subject to constraints of the form
T* 5 TI,TP, TI 5 7'*
r I
2
Figure 5.
t
r
tl
3
4
5
c4:
I t can be shown from the optimal value of 01 from our point of view tions. Let LIS suppose zero. ' h e n
Process signal flow diagram
l'j
(68)
20
physical must be to shoic that the
(691
arguments that if $ ' ( a ) > 0, zero. hut it is more interesting this from the necessary condioptimal u 4 is finite and is not
Figure 4 and represents the movement of information rather than material. 'The t\co serial subprocesses consist of stages 1 to 3 and stages 4 and 5. I t is convenient to avoid renumbering and rather to denote components at the latter two stages by a caret ( A ) , I n particular, Equation 45 becomes
I,;($
- L;,
=
0
(50)
where
and from Equation 70 it follows that a n d Equation 48 becomes
b _- H5 - dVj
If we denote the components of the Green's vector corresponding to a,, b,. and u, by cyrL, p,. and v,. respectively. then the transformation Equations 6 are as follows :
cy4
=
6,
=
1 252
1
+ u j 4' (&)
#O
Thus the optimal u s must be a t the constraint
US
=
0 and
(73) or
pn-l
li. (25)
A
cy5
&>0
(74)
; n = 4, 5
+
fl5 Z'j
IBEC FUNDAMENTALS
$'(id)
(75j
T h e same physical parameters were used as in previous papers ( 9 , 70):
and Equation 70 becomes
Thus. ti4 < 0. which is impossible, and the optimal value of u p must al\vays be zero. If we add Equations 6 2 and 6 5 ive obtain ,j
+
a3 = p
(77)
k10 = 5 X 101'Jlitersmole-lmin.? E1 = 18 kcal. T* = 335' K . 8 = 2 min: a, = 1 mole liter
= 3.33 X 10'7 min.-l E P = 30 kcal. T * = 355' K . p = 0.3 k20
T h e function + ( a ) is shown in Figure 6 and has the analytical form
5 a 5 0.6 1.08 - 1 . 1 ~ 2 ; 0.6 < (I
$(a) = 2.5a - 2u2; 0
Tow. let us suppose that we assume a value a3* and compute the sct { irL*. a n * / . T h e correct solution is then of the form {in* BEYn. a,* ~ c u , , ] . From Equation 77 we have
=
+
< 0.9
+
+
6ff5
= -8ff3
(78)
and thus from Equations 5 5 and 56
where u denotes
I'~.
.4lso, from Equation 58
But it folloivs from Equation 65 that
and combining Equations 79 to 81, it follows that the correct value of a3 is 013
=
as*
Figure 6.
+ a0* -
a3*
0
(1
Solute distribution function
&* +7
Once we have chosen the temperatures and solvent ratio Equation 82 is exact. and the correct Green's vectors may be determined from Equations 53 to 67 by a single application of E;quarion 82. Thv metric tensor at each stage reduces to a single scalar component, which we will call l / G n . If we assume that k i ( T,) has Arrhenius form k l , exp ( - E,/RT,) and let d =
a
.4
344
.3
343
.2
342
klan2
+ 20 k i ~ , ) R 7 ' , ~
and
the gradient directions are
.I
0
!
341.
3401 I I 0
Figure 7.
I I 5 10 ASCENT NUMBER
I
15
Successive approximations for y = VOL. 4
NO. 3
AUGUST
0.30
1965
253
T h e cost of solvent, 7, was allowed to range over all values, and GI, Gf: GI were set equal to unity. As in (70)>the initial step size was arbitrarily set at 2, with the step size halved for each move not resulting in an increase in P and a step size of less than taken as the convergence criterion. T h e rate of convergence was sensitive to the value of Gs used. For example, at y = 0.19, a change of G 5 from 0.05 to 0.01 resulted in convergence in 32 ascents instead of more than 140 from the same starting conditions. Preliminary experimentation indicated that a value of G5 = 0.005 Lvould he satisfactory, and this value was used for most of the calculations. T h e recycle material balance was solved by a onedimensional direct search, and a complete step, including the material balance, calculation of the Green’s vector and gradient direction, and determination of the new decisions, required 29/60 second \vith a Fortran program on a CDC 1604 computer; somewhat less time was required for a re-evaluationof the state with a smaller step size. T h e successive approximations shown in Figure 7. for y = 0.30, are typical, with 17 ascents and 28 evaluations of the state required for convergence. using less than 14 seconds of computing time. T h e profit is plotted in Figure 8 as a function of 7. T h e horizontal line corresponds to the optimal nonrecycle solution found previously (9, 70)> and in the neighborhood of the intersection two locally optimal solutions were found. Figure 9, for example, shows another set of successive approximations for y = 0.30, starting at the same initial temperature policy as in Figure 7 but at a different solvent ratio. T h e resulting temperatures are those for the optimal three-stage serial process. T h e optimal temperature and solvent policies as functions of y are shown in Figures 10 and 11, For sufficiently inexpensive solvent the temperatures go to T* and for sufficiently costly solvent the optimal policy is nonrecycle. Multiple
solutions are shown by dashed lines. and the slopes of the unconstrained solutions appear to become infinite. S o unconstrained solution could he obtained beyond a value of y = 0.31381. I t is clear from Equation 49 that there is a mathematical solution to the optimization problem corresponding to c = - m , and this may he verified with the necessary conditions. T h e sharp transition Lvould appear to he an indication that beyond a critical value of y the unconstrained solution becomes unique, and the constraint has the effect of translating this solution from LN = - rn to I = 0. Alternatively, we might say that beyond a critical value of y there is no physically meaningful unconstrained solution. A Mixed Example
T h e process shown in Figures 4 and 5 may he converted to a mixed system by substituting a plug flow tubular reactor of residence time 0 = 6 minutes for stages 1 to 3. T h e kinetic equations for the continuous serial process are then
(88)
u = -k,(T)a2
b
=
kl(T)a2- ka(T)b
(89)
(90)
k = 0
and the corresponding equations for the Green’s vector are
P)
dc = 2 k l ( T ) a(a -
P
=
(91)
k2(T) P
(92) (93)
u = o
If we retain the numbering n = 4, 5 for the extraction stages and denote the effluent at t = 0 by the subscript n = 3, then Equations 41 to 57 and 60 to 79 remain valid. T h e correct
Y
0
0
Ga c
W z
30 v)
.75
.50
.25
I
.L I
.I5
.20
.25
.30
COST OF SOLVENT
Figure 8. 254
Profit vs. cost of solvent
I&EC FUNDAMENTALS
!---.----.---.----.--.--.V
342
0
.35
I
I
I
I
I
I
I
I
0
1
2
3
4
5
6
7
ASCENT
9. Successive approximations for 0.30
Figure
y
=
NUMBER
--r
3451
I
0
.IO
$15
.20
.25
2 3 4 RESIDENCE TIME (MIN.)
5
6
Figure 1 2. Successive approximations to tubular temperature profile for y = 0.25
.30
COST OF SOLVENT
Figure 1 0. Optimal three-stage temperatures vs. cost of solvent 1.11
I
.I5
I
I
.20
.30
.25 COST OF SOLVENT
Figure 1 1 . solvent
Amount
ASCENT
NUMBER
Figure 13. Successive approximation of solvent for y = 0.25
of solvent vs. cost of
value of the Green’s vector may still be obtained with a single correction after choice of z! and the temperature profile T ( t ) , but instead of Equation 82 we must use a(0)
=
a*(0)
+ a*(O)
- .*(e)
+9
._
1+-U$’(Uq)
1
[l
+
v$’(a4)]
and the gradient for solvent from Equation 32 as
-* ff3
- exp [ -
Gsds
du
f
(97)
(94)
2akl(T) dt]
If we let l j G ( t ) denote the metric at each point of the continuous process a n d
then the gradient direction a t each point of the tubular reactor follows from Equation 33 as
where dj is defined by Equation 84. T h e initial step size was arbitrarily set a t 8>and G ( t ) set equal to unity. T h e state and temperature were stored at intervals of 0.1 minute, with linear interpolation used for intrrmediate values, and Equation 96 was applied at only these points. All integrations were carried out using a RungeKutta-Gill routine. A single ascent calculation rrquired slightly less than 18 seconds of computing time and a reVOL. 4
NO. 3
AUGUST
1965
255
’
evaluation of the state with a smaller step size required approximately 14 seconds, with over 11 seconds needed in each case for the one-dimensional search for the solution of the material balance to a n accuracy of 10-j. Different values of the nominal decisions for the start of the successive approximation solution usually resulted in slightly different solutions, although the resulting profits never differed by more than O.lyo. This fact could be due in part to the effect of errors of the order of 10-5 in the material balance equations, as well as to the usual difficulties associated with the gradient method in the neighborhood of a solution Little difference was observed in the rate of convergence between values of Gb of 0.005 and 0.0005, although the smaller value tended to result in slightly higher profits. A typical convergence sequence is shown in Figures 12 and 13 for y = 0.25. T h e solid nominal curve is the optimal nonrecycle solution found in ( g ) , while the shape of the dashed nominal curve was dictated by the fact that the unconstrained solution for the temperature profile can be shown from the necessary conditions to require a n infinite slope a t t = 0. Convergence was obtained in 9 ascents from the solid nominal curve and 11 ascents from the dashed curve, with an additional 13 state evaluations required in each case for changes in step size. T h e solutions are not identical, and the apparent discontinuity in slope a t t = 0 1 is a consequence of the use of a finite number of storage locations and linear interpolation. T h e profit is shown as a function of in Figure 8, the solvent allocation in Figure 1 1, and the temperature profiles in Figure 14. Convergence was always obtained within 8 to 15 ascents from either of the nominal solutions shown in Figures 12 and 13, and as in the staged system, multiple solutions were found. T h e same profit was obtained for y = 0.325 by the curve shown in Figure 14 and solvent in Figure 11 and by the nonrecycle solution, shown as y = co . No attempt was made to determine the exact range of multiple solutions. but it was found to extend to a t least y = 0.33 and did not extend back so far as y = 0.30. T h e failure of the curve for y = 0.05to show an infinite initial slope is a consequence of the discrete storage of temperatures. 355
m
formulation of the problem, in order to facilitate solution of the boundary value problem for the Green’s vector. As we have demonstrated, the use of the weak maximum principle may simplify the problem somewhat. Multiple solutions must be expected in complex systems in which the amount of recycle or bypass is itself a decision which has monetary value. Thus, solutions must be sought from a number of starting points, increasing the required computing time. If, however, as in the problem which we have studied in this paper, the suspected additional solution can be identified, it may be used as a starting point itself and the number of calculations significantly reduced. Acknowledgment
We are grateful to the Numerical Analysis Center of the University of Minnesota for the generous grant of computer time. During the period of this work Morton M. Denn was a National Science Foundation Cooperative Fellow. Nomenclature = concentration of = concentration of
a
b d
=
= = = = = = = = = =
Ei
= = =
9 Qk
R
=
S
= = =
S t
T
=
T
=
T*!T c
=
U
= =
U
.4 B
unnormalized Euclidean gradient activation energy metric tensor associated metric tensor inverse of scalar metric Hamiltonian temperature dependent term in reaction rate constant in Arrhemius expression number of stages state vector state of kth feed to the system profit function decision vector constraint gas constant distance in decision space dimension of p timelike variable temperature differential or stage transformation bounds on T augmented state variable solvent ratio
GREEK = component of a corresponding to a = assumed value of a = component of tj corresponding to 6
a a*
B
relative value of solvent to B error in a 6p. 6q, 6P variations n = multiplier = process or holding time e x = multiplier ,iz(1,k)j = dpoi‘ d@P\,’ U = Green’s vector m(h)l = Green’s vector associated with rth component of kth effluent = relative value of A to B P = component of u corresponding to u L’ $ = solute distribution function
Y
I
I
0
I
2
3
4
5
6
RESIDENCE TIME (MIN.)
Figure 14. of solvent
Optimal temperature profilss for various costs
6a
= = =
62‘(l.k,j
=
Discussion
We have shown how complex chemical plants may be optimized with a reasonable expenditure of computing time and small storage requirements, with the greatest amount of time required for the solution of material balance equations. T h e gradient method provides a powerful technique, but care must be used in the choice of the proper metric and even in the 256
l&EC
FUNDAMENTALS
~ P Q ‘ ab/,’
Literature Cited
( I ) , Ark, R., “Discrete Dynamic Programming,” Blaisdell Publishing Co., New York, 1964. (2) Aris, R., Nemhouser, G. L., Wilde, D. J., A.Z.Ch.E. J . 10, 913
I1 964). - , (3) Bryson, A. E., Denham, W. F., Dreyfus, S. E., A.Z.A.A. J . 1, 2544 (1963). \ - -
‘4) Chien. H., personal communication. t5) Denham. W. F., Bryson, A. E., A.I.A.A. J . 2, 25 (1964). : 6 ) Denn. M . M., Ph.D. thesis. University of Minnesota, 1964. 17) Denn, M . M., .4ris, R., A.I.Ch. E . J.. in press. 18) Denn. M. M.: ilris. R., Chem. Eng. Scz., in press. 4, :9) Denn. M . M . , Aris, R., IND. Elic. CHEM.FUNDAMENTALS 7 (1965\. [IO) ‘Zbtd.: p. 213. (11) Dreyfus. S. E.. “Variational Problems with State Variable Inequality Constraints,” Rand Rept. P-2605-1 (ilugust 1963). (12) Fan, L. T., Wang. C. S., J . Appl. Math. Phys. ( Z A M P ) 15, 46 11964). (13) Fa;: L. T., Wang, C. S., “The Discrete Maximum Principle, Wiley, New York, 1964. (14) Horn, F., Jackson, R., IND. ENG. CHEM.FUNDAMENTALS 4 , 110 (1965).
(15) Horn, F.: Jackson. R.. J . Elec. Control, in press. 116) Jackson. R.. Chem. Ene. Scz. 18. 215 11963). (17j Ibzd., 19, 19 (1964). (18) Ibzd.. p. 253. (19) Kellcy. H. J., in “Optimization Techniques with Applications to .kro-saace Svstems.” G. Leitmann.. ed.., .Academic Press, New York‘, 1962. ’ (20) Pontryagin, L. S.,Boltyanskir, V. A., Gamkrelidze, R. V.. Mishchenko, E. F., “The Mathematical Theory of Optimal Processes,” IGiley, S e w York, 1962. (21) Yesberg, D., Ph.D. thesis, University of Micnesota, 1964. \
,
RECEIVED for review September 22, 1964 ACCEPTED December 28, 1964 Third in a series on Green’s Functions and Optimal Systems.
OPTIMAL REACTION SYSTEMS WITH RECYCLE G . S. G. B E V E R I D G E ’ A N D R . S. SCHECHTER
Department of Chemical Engineering, The University of Texas, Austin, Tex.
A mathematical technique has been developed for the determination of the optimal operating policy for the system involving a tubular chemical reactor, a separator, and a recycle stream. The analysis uses the concepts proposed b y Pontryagin, and the results show that one need only modify the definition of the adjoint functions in order to treat such problems. The method is illustrated b y an example to find the maximum profit for an irreversible reaction carried out in a tubular reactor with recycle,
HE use of recycle as a means of enhancing plant profitT a b i l i t y is a technique often applied in the process industries. T h e mathematical problem of determining a n optimal operating policy for a recycle system is, however, difficult and has only recently started to receive attention (7-5, 70). I t is intended not to review these contributions critically, but to demonstraLe the solution of the general recycle problem using a modification of Pontryagin’s maximum principle ( 9 ) . T h e system considered consists of a mixer, a tubular reactor, and a separator arranged as shown in Figure 1. Fresh feed is introduced into the reaction system a t a constant rate, q, with a fixed composition, X ( Z ) . [ X ( p )is used here to denote the set of N weight fractions (XI? x p ? , . . , x S ) which completely 1) component stream, specifies the chemical state of a n (A’ p . ] T h e separator following the reactor produces a product stream and a recycle stream a t the respective flow rates, q and I, with compositions A’(!) and X ( r ) . I t is assumed that the system already exists, so that the reactor operation can be varied only by manipulating the external operating variables, Y ( I ) . There are k of these variables which can be altered independently along the length of the reactor. T h e split in the separator can be characterized by the nonnegative separation factors, S, where
+
5,
= x , ( ~ ) ! ’ x ~ ( f )j; =
I, 2, . . . N
(11
so that from a mass balance over the separator we can write the composition of the j t h component in the recycle stream, x,(r) as X,(‘)
’
(I
+ q/s,)
= (r
+ q)x,(L)
S denotes a set of separation factors, SI,
52,
..
(2)
.sN.
Present address, Department of Chemical University of Edinburgh. Edinburgh, Scotland.
Engineering,
Since the feed to the reactor, X ( O ) , is the sum of the contributions from the fresh feed and the recycle streams, it follows from Equation 2 that X,(O)
= __ x,(O
4+r
+ + Y/J,) 7
x,(L); j
= 1,
2,
. . .N
(3)
(1.
T h e composition within the reactor. X ( I ) , will vary along its length, the variation being determined by the AV kinetic equations
T h e reaction rate depends upon the local compositions, X ( I ) , the level of the k control variables, Y(1):and the mass flow rate, ( r q). For a given feed rate of a fixed composition, a given set of the operating variables, recycle rate, and separation factors will then uniquely determine the product composition. W e may call such a set an operating policy, P, which we denote as
+
P [ Y ( O ,s:
(3
71
Having completed a design basis for the problem, it is now necessary to obtain the best design. T h e criterion for the determination of a ”best’‘ design varies; it may be a maximum yield, a maximum cost: a maximum profit, but in general it can be represented for the present system by a functional, I , defined in the form
Z(P) =
l-
c[A’(i)> Y ( i ) ,r . Sldl -t C [ X ( L ) :S. 71
(6)
Z is a functional of the policy, P. An illustration of the typical components in this relation is given in a n example belolc, VOL. 4
NO. 3
AUGUST
1965
257
