(5) Carman, P. C., “Flow of Gases through Porous Media,” Academic Press, New York, 1956. (6) de Boer, J. H., “The Dynamical Character of Adsorption,” Clarendon Press. Oxford. 1953. (7) Emmett, P. H., “Catalysis,” Vol. I, Reinhold, New York, 1956. (8) Funnel, W. S., Hoover, G. I., J . Phys. Chem. 31,1099 (1927). (9) Gifford, P. H., M.S. thesis, Ohio State University, 1958.
(19) Rivarola, J. B., Smith, J. M., Znd. Eng. Chem. Fundamentals 3,308 (1965). (20) Sashihara, T. F., Ph.D. dissertation, Ohio State University, 1957.
(21) Sashihara, T. F., Syverson, Aldrich, IND. ENC. CHEM. PROCESS DESIGN DEVELOP. 5,392 (1966). (22) Satterfield, C. N., Sherwood, T. K., “Role of Diffusion in Catalysis,” Addison-Wesley, New York, 1963. (10) Gilliland, E. R., Baddour, R. F., Russell, J. L., A.Z.Ch.E. J . (23) Sutherland, K. L., Winfield, M. E., Australian J . Chem. 6, 234 4, 90 (1958). (1953). (11) Krasuk, J . H., Smith, J. M., IND. ENG.CHEM.FUNDAMENTALS (24) Ibid.. D. 244. 4, 102 (1965). (25) Tam&, K., Advan. Catalysis 15, 65 (1964). (12) Livingston, R., J . Phys. Chem. 33, 955 (1929). (26) Tamaru, K., Trans. Faraday SOC. 55, 824 (1959). (13) Low, M. J. D., Chem. Rev. 60,267 (1960). (27) Zbzd., p. 1191. (14) Macarus, D. P., Ph.D. dissertation, Ohio State University, (28) Zbid., 57, 1410 (1961). (29) Taylor, H. A,, Thon, N., J. Am. Chem. SOC.74, 4169 (1952). 1959. (15) Masumune, S., Smith, J. M., A.Z.Ch.E. J . 11,34(1965). (30) Winfield, M. E., Australian J . Chem. 6, 221 (195.3). (16) Ibid., p. 41. (17) Mennessier, A., Boucher, R., Compt. Rend. 226, 914 (1948). (18) Porter, F. P., Bardwell, D. C., Lind, S. C., Znd. Eng. Chem. RECEIVED for review December 13, 1965 18, 1086 (1926). ACCEPTEDJuly 5, 1966
OPTIMIZATION OF COMPLEX CHEMICAL PROCESSES Comparison of Variational and Steepest Descent Methods M. T. K U O , D. I . R U B I N , AND B. S. W R I G H T American Cyanamid Go., Wayne, N . J .
Two of many approaches for solving optimization problems arising from chemical process improvement are discussed: steepest descent and the variational method. The results obtained from the steepest descent technique are extended to the variational problem, and the computational algorithms are discussed. Finally, both methods are assessed. HE numerical techniques for the solution of optimization Tproblems can be classified into direct and indirect methods. I n the direct approach, such as the steepest descent techniques, the system differential equations and the constraints are satisfied and a n iteration is made on the policy function such that each new iterate improves the function to be maximized (or minimized). The indirect approach involves the development of an iterative technique for the solution of the system and the adjoining differential equations. Detailed discussions on various aspects of optimization, both theoretical and numerical, have been published (7-6, 8, 7 7). I n this paper, the approximate solution is first obtained by the steepest descent technique and the result is extended to the variational problem. We consider a semiflow batch reactor in which the following complex homogeneous reaction is involved (7). A fixed amount of reactants is continuously fed into the reactor for a certain period of time. We seek the best feed addition rate (policy function) such that the yield is maximum.
kn
B+C----tG n
E + G d I
+
Component (F H) is the product and component I is the by-product. Component E is the hydrogen ion. We could analyze for:
F
+H
(product)
I
(by-product)
A+D B+C
E Based on the above reaction sequences (7), the following rate, stoichiometric, and equilibrium relationships are derived :
ki
AfB-F
rz
A + E d D trl I
B-~--~c+E f-I, id
E + F d H t I6 1
404
Present address, FMC Co., Princeton, N. J. I & E C PROCESS D E S I G N A N D DEVELOPMENT
c 5
= 10-=
.
where Ci is the concentration of component i (z' = 1, 2, . . 9 ) , the subscripts being ordered as the alphabet. The rate, stoichiometric, and equilibrium relationships represent seven equations in the seven unknowns: C1, C2, CS,C4, CS,(Ce Cd, and Cg. From small scale batch reactor experiments over the expected range of temperature, pH, and initial concentrations, the rate and equilibrium constants were estimated, so that the sum of the squares of discrepancy between the measured and predicted concentrations was a minimum. We need not concern ourselves with parameter estimations. The sums of discrepancy were small and well within experimental error.
where eS is the Inole fraction of component Cgin the feed. (Ce
(Cs)Il = 0
+
where
d -
dt
V(C1 V(Cz
+ Cdl
[VCg]
Bo
=
Y1
=
Yz =
We will consider the semifiow batch reactor in which B is added over a prescribed time. The amount of C in the B feed is negligible because the p H of solution B is kept low. However, the reaction mixture must be maintained at a considerably higher pH, so that the reaction proceeds at a reasonable rate. Let us now consider a semiflow batch reactor, depicted in Figure 1. To component C1 an aqueous solution of Cz is added. The solution of Cz contains a small amount of byproduct, Cg. The reaction is over the range 0 6 t 6 T and the feed addition time is from t = 0 to t = 7,where 7 6 T. The batch reactor Equations 2 may be modified for the semiflow batch reactor as follows: dt
vo
V(0) =
System Equations
d - [V(C,
+ Cdo = 0
=
VklC1Cz; (Ce
VkzC2C3
+
+ c4 + Ce + Cd + Ce + Ca + 2
f '23
egQ;
=
+ Cs)o
(3)
=
Equations 5 through 6 indicate that the performance of the system depends on the addition schedule of
Jo Feed addition rate, O(t), is to be determined to maximize V(Ce Ca) at t = T.
+
Technique of Calculus of Variations
The computational scheme is the calculus of variations methods developed by Katz ( 5 ) . The essence of the technique as applied to this problem is that for a given set of differential equations in the dependent variables 2 = (XI, xZ, . . . x N )
(Cg)o = 0
VO[(Cl)O
C9) =
= 0
(I1
so'
+ (C4)Ol
(Qz
i
x~(O) =
+ 2 QJdE
In addition, we need an equation which describes the volume change of the system due to the addition of feed.
=
1, 2, . . . N
the policy function O(t), over the range 0 6 t 6 mined that maximizes xl( T ) . Define the set of adjoints, Z = (21, 12, properties that
T,is to be deterI
. .zN), with
the
where m t is the molar specific volume of the ith component and Q = Q Z 2 Q g . , Rearranging the equation gives
+
d
+
- [~(ce ca)] = 2 [ P O dt
V
[l
Qd€
V(Ce
+ C ~ ) xI
- 2 VCg - v(c6
For each t ,
e is determined so as t o maximize the Hamiltonian N
Ca)]
(5) If 0 is restricted to lie within given bounds-Le., then for each t , the value of 0 that is selected is:
(Y
6 e < @--
o = a
or
e=p or the solution, 0, of
Figure 1. Diagram of semiflow batch reactor
whichever maximizes H. The 2N differential equations defined by Equations 8 and 9 are to be treated as a boundary value problem by virtue of the VOL. 5
NO. 4 O C T O B E R 1 9 6 6 405
initial conditions in 8 and the final conditions in 9. The policy function, e(t), is determined by integrating the 2N equations using the iterative initial value procedure described below. Let
Q =
(41,49, . . .h), where 4f
= 0
i = 1, 2,
...N
z =
fZ(xl, x Z ;
(12)
Choose an initial guesstimate to Q.
B. Integrate the 2N differential equations from t = T to t = 0, where x t ( T ) = &
xl]z =
(mz
Xl(0)
= 0
XZ(0)
= 0
e)
(15)
subject to the conditions
The solution to Equation 12 is found in the following manner: A.
+ + ms)B [(l - 2 4 e - 2 x 2 -
{
Q is the solution of the nonlinear system of N equations
- at
YZ
VO
dt
= xf(T)
4
gt(@ = xi@)
-dxz =
dxi./dt, dzt/dt
O(t) =
L
t = O
monotonic increasing
0
t >I 0.8 is the solution of hH/&3 = 0. This reduces to the soiutign of a quadratic equation in e at the end of each integration step. However, 0 must satisfy the monotonicity h). constraint 0 6 e ( t ) 6 e(t If the solution of the quadratic yields two positive roots, both satisfying the monotonicity constraints, the one giving a larger Hamiltonian is selected. If there are no positive roots, 6 is unchanged for that step. F o r t = 0, e = 0.
+
Steepest Descent
and =
Equations 14 and 17 are the system of differential equations as adapted for the variational method.
The basis for the steepest descent solution is Marquardt’s algorithm (9). The range 0 t r is divided into m equal intervals defined at
J‘ Q W E
<
