SUBSCRIPTS A N D SUPERSCRIPTS
exp = experimental L = liquid pred = predicted V = vapor O = pure component property
literature Cited
Beattie, J. A,, Marple, S. J., J . Am. Chem. Sac. 72,4143 (1950). Benedict, M., Webb, G. B., Rubin, L. C., Chem. Eng. Progr. 47, 413 (1951). Benedict, M., Webb, G. B., Rubin, L. C., J . Chem. Phys. 10, 747 (1942). Beyer, H. H., Griskey, R. G., A.Z.Ch.E. J . IO, 764 (1964). Canjar, L. N., Smith, R. F., Volianitis, E., Galluzzo, J. F., CabarCOS,M.,Znd. Eng. Chem. 47,1028 (1955). Esso Research and Engineering Files, unpublished data, 1958. Griskey, R. G., Canjar, L. N., A.1.Ch.E. J . 9,182 (1963). Joffe, J., Chem. Eng. Progr. 45, 160 (1949). Lin, M., Naphtali, L. M., A.Z.Ch.E. J. 9,580 (1963).
Meissner, H.P., Seferian, R., Chem. Eng. Progr. 47, 579 (1951). Opfell, J. B., Sage, B. H., Pitzer, K. S., Znd. Eng. Chem. 48, 2069 (1956). Pitzer, K. S., Lippmann, D. Z., API Project 50, University of California, Berkeley, 1954. Rossini, F. D., Pitzer, K. S., Arnett, R. L., Brown, R. M., Pimental, G. C., API Project 44, Carnegie Press, Pittsburgh, 1953. Sage, B. H., Lacey, W. N.,Znd. Eng. Chem. 40,1299 (1948). Scott, R. B., Myers, C. H., Rands, R. D., Brickwedde, F. G., Bekkedahl, N., J . Res. Natl. Bur. Std. 35,39 (1945). Smith, C. H., Thodos, G., A.1.Ch.E. J . 6 , 569 (1960). Su, G. S., Sc. D. thesis, Massachusetts Institute of Technology, Cambridge, Mass., 1937. Su. G. S.. Viswanath. D. S.. A.1.Ch.E. J . 11.205- 11965’1. Stark, T.’M., J . Chem’. Eng.’Daia 11,570 (1966). Stark, T. M., Joffe, J., J . Chem. Eng. Data 9 , 327 (1964). Union Oil Co. of California, Research Department, unpublished data, 1944. Vohra, S. P., Kang, T. L., Kobe, K. A., McKetta, J. J., J. Chem. Eng. Data 7,150 (1962). Waals, J. D. van der, dissertation, Leiden, 1873. Zudkevitch, D., Kaufmann, T. G., A.Z.Ch.E. J . 12, 577 (1966). - 7 -
\ -
- I .
- .
RECEIVED for review March 9, 196 ACCEPTED July 25, 196
ESTIMATION OF RATE CONSTANTS FROM MULTIRESPONSE KINETIC DATA R E l J l M E Z A K I A N D J O H N B. B U T T Department of Engineering and Applied Science, Yale University, New Haven, Conn. Alternative procedures for the estimation of several parameters from multivariate observations are evaluated for a complex system. Primary attention is given to methods based on generalization of nonlinear least squares and Bayesian estimation as proposed by Box and Draper. The latter method provides an effective means for estimation even in very complicated cases, and demonstrates a parametric sensitivity which allows a precision of estimation not possible with other techniques.
HE
establishment of systematic procedures for the deter-
T mination of constants involved in many-parameter, nonlinear models, and for choosing between alternative models, has assumed considerable importance in the current literature. Various rate schemes for chemical reactions (both homogeneous and heterogeneous) have provided the favorite means for investigation of these procedures. Development and application of both linear and nonlinear estimation procedures have been documented (Box, 1960; Hougen and Watson, 1947; Kittrell et al., 1965; Lapidus and Peterson, 1966) and extensively discussed. Nonlinear estimation of several parameters from data ordinarily available for these reaction systems is sometimes difficult. These difficulties often stem from a lack of sufficient information regarding the reaction system; however, in many reactions of practical significance a number of reactants, intermediates, and products may be involved and a large amount of information is potentially available if only it is measured. It is, thus, of interest to investigate estimation procedures applicable when more than a single response variable is measured and to compare the results of specific multiresponse techniques with those obtained by alternative methods. 120
I&EC FUNDAMEN.TALS
Estimation Procedures and Criteria
T h e principal procedures for determination of parameters in multivariate systems which may be nonlinear have been summarized as follows (Hunter, 1967): Generalization of Nonlinear Least Squares. The best fit criterion is given by :
This criterion is the proper one only under certain restrictive conditions, as discussed later (Hunter, 1967). A corresponding weighted least squares criterion can be applied if the variance of the data is known. This ordinarily is not the case, so comparison with such a weighted criterion is not included in the present work. Bayesian Estimation. When the values of variance and covariance are unknown and a number of variables have been measured, best estimation criterion is given by :
min k
\ II
u-1
T h e model postulates a series of four surface reactions:
or Minimization of the determinant of Equation 2 represents the most general criterion for fit of multiresponse data (Box and Draper, 1965; Hunter, 1967). T h e confidence region of constants so estimated is given by (Hunter, 1966) :
CzH50H
+ S$ C + HzO
( K I ~Kir) ,
C Z H ~ O C Frt ~ HC~Z H ~ O H fC
( K z ~KZr) ,
+ A'
(K3)
CzH5OCzHs f 2S-t 2C
f HzO
C + CzH4
(K4f)
and the appropriate rate equations are given by: in which the numerator of the exponential argument is the abscissa value below which is 100 (1 a) % of the area under the curve of the chi-squared distribution with m degrees of freedom. Box and Draper (1965) have shown that the confidence region of parameter estimates from multiresponse data employing the determinant minimization criterion will be smaller than that associated with alternative estimation procedures.
-
dA = RZ
- Ri
dt
dO dt
=
R3
Application to Complex Catalytic Reaction Sequence
Solomon (1966) and Solomon et al. (1967) have recently reported an extensive study of the kinetics of ethyl alcohol and diethyl ether dehydration over y-alumina, encompassing a 150' C. temperature range and pressures below 1 atm. Two principal reaction:) occur simultaneously (ether and ethylene formation) and the selectivity between them changes over the temperature range 200' to 400' C., ether formation predominating a t low trmperatures. For experimental conversion measurements a batch recycling reactor was used from which gaseous samples were removed for periodic analysis. Raw data so obtained were composition of the reacting mixture us. time of reaction; these data could then be transformed into partial pressure-time information. T h e stoichiometry of these reactions with alcohol feed is given in terms of two independent reactions:
+ HzO CzHsOCzHa + HzO
CzHQH @ CzH4
We are interested in applying the estimation procedures previously outlined to the evaluation of rate constants associated with this complex model, in comparing the results so obtained by means of the fit to experimental data, and in discussing some operational problems involved in applying these procedures to multiresponse systems. T h e computational methods employed require numerical solution of the rate equations for the set of constants, evaluation of the criteria given by Equations 1 and 2 from computed and experimental values, changing the constants as indicated by some minimization procedure, and then repeating the numerical solution.
2CzH60H[~2
Parameter Estimation Results
Over-all stoichiometry, however, does not describe the reactions on the catalytic surface, and Solomon employed a kinetic model for dehydration derived from mechanistic proposals that the reaction proceeds via the formation of an ethoxy complex on the surface after alcohol chemisorption on a surface site of the form O=Al-OH:
0
II
AI-H
H
H
I I f HO-4-CH I I H
H
0
Ii
Al-OCzHs
+ HzO
The experimental data reported (Solomon, 1966) were obtained under conditions as specified in Table I, and constants were estimated for all runs. In Table I1 are given the rate constants estimated for 329' C., a representative run. The data employed for the estimation procedures were on alcohol, ether, olefm, and water compositions in the reaction mixture as a function of time for both alcohol and ether feed experiments. The three sets of constants given are those determined by Solomon, those determined by nonlinear least squares, Equation 1, and those determined by the determinant criterion, Equation 2. VOL 7
NO.
1
FEBRUARY 1968
121
Table 1.
Catalyst.
Alcohol feed
Ether feed
Experimental Conditions
5.002 grams of ?-alumina, A1-1404 mfg. by Harshaw Chemical Co. BET area of 190 sq. meters per gram O c. Mm. Hg 220 290 initial pressure 254 155, 230, 310, 316, 319, 321 291 308 326 306 304 329 182 374 261 241 291 183 306 207 191 329
0.8
WATER
n J
f 0.6
--
4
EQN ( 2 )
-----
ALCOHOL
4
0
20
10
30
40
60
50
70
80
90
TIME, MINUTES W
I
The fit to experimental data by the model using the rate constants estimated by the two methods is given for the example run with alcohol feed a t 329' C. in Figure 1. There is little to choose between the two on the basis of visual comparison, other than perhaps a little better fit to the ether maxima by Equation 2. The rate constant estimates differ, though, as shown in Table 11, and it is shown in Table I11 that minimization of the determinant does not minimize the sum of squares. The visual comparison is qualitatively informative concerning the relative precision of the constants estimated by the various methods but is not sufficient to demonstrate the superiority of Equation 2. More quantitative evidence on this point is provided by considering in detail the sensitivity of the fitting criterion to the parameters. Minimization procedures when applied to multidimensional problems cannot generally be depended upon to determine global minima unless the response surface has particular characteristics. Particularly troublesome is the fact that some parameters are likely to be much more important in determining fit than others, and unless the criterion is fairly sensitive to all parameters, spurious values may be obtained for some constants. The present system provides a fine example of this type of parametric insensitivity, since the conversion curves are largely dictated by the values of K1,S0, K2,So, and &So, while the reverse rate constants may be varied within a relatively wide range with little visual effect on the fit to the data. In Figure 2 is given the minimization history of Equations 1 and 2 for the estimation at 329' C., using the constants reported by Solomon (1966) as starting values. T h e minimization routine employed (Sobol, 1963) was not a gradient technique, but determined the behavior of the function with respect to each parameter and adjusted individual increments separately. This routine was tested for the six-dimensional case using generated data and it was demonstrated that the method retained its effectiveness if the response surface was sensitive to some parameters and insensitive to others. The validity of the minima shown, at least as local minima, was
Criterion Solomon Equation 1 Equation 2
30
40
50
lbEC FUNDAMENTALS
60
70
a0
90
TIbiE, MINUTES
Figure 1.
Fit to conversion data at 329'C. Upper. lower.
Alcohol and water Ether and olefln
established by pattern search around the indicated final values. I t is clear from Figure 2 that the determinant of Equation 2 is considerably more sensitive over-all to the rate constant estimates than is the generalized least squares criterion; from start to finish lVIJl varied by a factor of 16.2, while the sum of squares decreased only about twofold. I n the most pronounced case (alcohol conversion at 374' C.) the determinant decreased by a factor of 10l2while the sum of squares decreased only by a factor of 5. The initial sharp decrease of the determinant shown in Figure 2 (up to about step 150) is the result of changes in the values of the more sensitive From this point on parameters, K1,So, K2,So, and K& these constants were changed relatively little by the search routine and subsequent minimization of the determinant was accomplished through changes in KlrSo,K.&,, and Kd/S,2. By contrast, the minimization of least squares is not nearly SO selective with respect to important parameters-that is, over the 500 minimization steps shown in Figure 2 for the case of least squares all six parameters are being varied by the minimization routine throughout the course of the run. Major reductions in sum of squares such as those shown a t step 400 are still due to changes in the sensitive forward rate
Table II. Estimated Values of Rate Constants Fitted to Conversion Data at 329' C. Rate Constant, (G. Catalyst - Min.)-' K1/So Kd'o K2lSo K4'o 1.96 X 8.03 x 10-2 2 62 j< 3.80 X 10-' 1.99 X 1.99 X 10-1 1.08 X 6.57 X lo-' 2.34 X low2 1.20 x 10-1 1.51 x 10-2 3.80 X 10-'
Rate Constant, ( C . Catalyst Criterion Solomon Equation 1 Equation 2
122
OLEFIN
- Min. - Mm. Hg)-l K$6 0.88 X lo1 1.63 X 10' 1.13 X 10'
KalSo2 2.20 x 10-3 4.89 x 10-4 1.58 X lo-'
T h e improvement in precision of estimation is particularly notable for the "unimportant" constants. I t is shown below that an equilibrium constant involving KI,S, and K2,So is in much better agreement with the literature when Equation 2, rather than Equation 1, estimates of these constants are used.
Table 111. Values of Sum Squares and Deteriminant at Minimum Value of Fit According to Eq, 1 Eq. 2 Equation 1 3.96 X 7.37 x 10-12 Equation 2 5.23 X lo-* 1.72 X lo-'*
25.
Precision of Estimates
EON(5) MINIMIZATION
2 0 **
- 7 . 6:"
0
-
5 w 4 3 2 - 1
L 3 ' 0 100 200 300 400 500 600 900 800 900 Id00 NO OF bINIMIZATION STEPS
Figure 2. Minimization history for Equations 1 and 2 for the fit of Figure 1
constants-all this occurring a t a point where, by the Bayesian method, these parameters had essentially been fixed a t their final (best estimate) values for 250 minimization steps. This was the general pattern in comparison of the estimation techniques under all conditions employed in the study. T h e determinant of Equation 2 decreased by a large factor in relatively few initial minimization steps, and such decrease was due to changes in values of the forward rate constants to which the fit of conversion data is particularly sensitive. O n the other hand, minimization of least squares involved an ongoing variation of all parameters by the minimization technique with no initial large change associated with the forward rate constants. In addition, if the parameters corresponding to some point along the least squares minimization path are chosen as a starting point for Equation 2 minimization, the same type of rapid initial decrease in the value of the determinant is noted; this change is due to adjustment of forward rate constant values almost exclusively, and the new estimates of constants give better fit (in the sense of the equilibrium constant calculation described below) than the least squares values a t the starting point, or after the same number of minimization iterations from that point. This behavior does not appear to be a property of the particular type of minimization technique, since experiences using gradient minimization methods have been the same, but rather of the nature of the response surfaces generated by the two criteria. O n the basis of these comparisons one argument in favor of the Bayesian method, then, is that convergence to final estimates is much more rapid than with nonlinear least squares. There are two additional points of rather general importance associated with the sensitivity and convergence properties discussed : Multiparameter estimation problems may often be simplified. One might have expedited the Bayesian estimation illustrated in Figure 2 by treating the problem, in initial stages a t least, as consisting of two, three-parameter minimizations, the second with respect to the less important constants. A procedure similar to this, using nonlinear least squares for single response data, has been suggested (Lapidus and Peterson, 1965, 1966).
T h e results of estimation by Equations 1 and 2 for the entire body of data are given as Arrhenius plots in Figures 3, 4, and 5 and corresponding activation energies are reported in Table IV. In the present work it was somewhat simpler to treat separately data a t each temperature level and then employ a regression analysis on the assembled results, though in theory there is no reason why parameters a t all levels could not be determined simultaneously. The approximate 95% confidence region of these estimates was established by means of Equation 3. Since six parameters are involved in the scheme, it is not possible to define simply the region of specified confidence for the over-all fit with respect to a given individual estimate; however, these individual confidence regions include values obtained from the regression line in almost all cases. Theoretical discussions of these alternative methods (Box and Draper, 1965; Hunter, 1967) conclude that Equation 1 can be used for obtaining the best estimates of parameters only in the restricted circumstance in which one can assume the normal distribution of responses, the equality of all the variances, and the absence of correlation between the responses. If any of these conditions is violated, or the behavior of the variance or covariances of responses is not clearly known, the use of Equation 2 provides the best estimate of Darameters from multiresponse data.
O v)
IO2:\
L Y
L
374.C
329 306 291
25424
220
I/T Figure 3. and K&,
Arrhenius plot for Equation 2 estimates of
VOL 7
NO. 1
FEBRUARY 1968
KljS.
123
Temperature Dependence of the Rate Constants Estimated Activation Energies Constant - AE, Cal./G.-Mole
Table IV.
Eq. 1
Eq. 2
23,300 12,600 18,100 3,100 27,400 25,900
19,800 8,100 18,500 3,000 25,600 25,250
It may be estimated from other data (Pines and Haag, 1960) that the temperature dependence of So may be approximated by an “activation energy” of about -3000, so true values for the reactions should be greater by this amount than those shown. The figure for Kd,So* is, of course, 6000 greater.
That such theoretical conclusions are correct is illustrated most directly in the present study by the results plotted on the upper half of Figure 5 . The line represents values of the equilibrium constant for the reaction: 2CzHkOH
CzHsOC2Hs
+ HzO
as reported in the literature (Cope, 1964). This equilibrium constant may also be computed from values of the individual rate constants by the expression :
W t Figure 4. Arrhenius plot for Equation 2 estimates of and K4,So2
KJ,
0
X*lO.O
0
0
LITERATURE (CoPE,IP64
-0
1.0 1.0
/
&
Conclusions
- KrSo
Q
‘O0l
I 1 -
L 374 *C 32f 306291
254 241 220
I/ T
Figure 5. Arrhenius plot for Equation 2 estimates of KI& and K2$, 124
The value and the temperature dependence of K,,, as computed from Equation 8, will be sensitive to the estimated constants; hence comparison with literature data can be used to provide an independent check on the precision of estimation. In Figure 5, top, the open circles represent K,, computed from individual rate constant estimates according to the Bayesian method; closed circles represent K,, computed from least square estimates. It is apparent that Equation 2 values are in much better agreement with the literature, and use of the determinant criterion with the multiresponse data has provided a significant improvement in parameter estimates.
I&EC FUNDAMENTALS
The results of this work show that the multiresponse technique of Box and Draper provides a useful and relatively easy method for parameter estimation when multiresponse data are available. There is very little difference between Equations l and 2 as far as computational effort is required, since the major time requirement and the major source of difficulty are in computing values of 9 from the reaction model and in minimization interaction. The additional effort required to compute the off-diagonal elements of the determinant and its value is justified by the improvements in estimates afforded by the method. The parametric sensitivity of the determinant criterion allows precise estimation of constants which may have only a small influence in certain regions of the data fitted, a feature which can assume considerable importance in the analysis of chemical reaction rate data where rate constants for reverse reactions may be reliably determined from experiments which reflect predominantly the effects of the forward rate. I t was not the purpose of the present study to test the validity of the kinetic model employed with respect to other mechanistic proposals, since these matters have been treated elsewhere (Pines and Haag, 1960; Solomon et al., 1967).
Acknowledgment
The authors are indebted to the Yale Computer Center for the support of computation costs, and W. G. Hunter for his assistance and advice.
t Y
=
k 8,
- 0)
xm2(1
A , E, 0, Lt‘ = fractional pressure of alcohol, ether, olefin, and water respectively, with respect to initial partial pressures c, s = concentration of complex and unoccupied sites, respectively - AE = activation energy j , I , u, u = summing indices k = vector of parameters to be estimated K1/, Klv = forward and reverse rate constants, alcoholsite reactions Kz,, Kzl = forward and reverse rate constants, alcoholcompkx reaction K3 = forward rate constant, complex decomposition reaction K4/ = forward rate constant, ether-site reaction Ke, = equilibrium constant, Equation 8 m = degrees of freedom, (r - 1) n = number of experimental points P, = initial pressure of reactant R1 . . . Rd = rate terms defined in Equations 4 to 7 r = number of responses So = surface site concentration
computed value corresponding to y determinant defined by Equation 2 fraction of surface coverage by complex with respect to So = fraction of vacant sites with respect to So = abscissa value below which is 100 (1 - a) % ’ of area under the chi-squared distribution curve with m degrees of freedom
= = =
I VU
ea
Nomenclature
time
= measured datum
literature Cited
Box, G. E. P., Ann. N . Y. Acad. Sci.86, 792 (1960). Box, G. E. P., Draper, N. R., Biometrika 52, 355 (1965). Cope, C. S., A.I.Ch.E. J . 10, 277 (1964). Hougen, 0. A., Watson, K. M., “Chemical Process Principles,” Wiley, New York, 1947. Hunter, W. G., IND.ENC.CHEM.FUNDAMENTALS 6, 461 (1967). Hunter, W. G., personal communication, 1966. Kittrell, J. R., Hunter, W. G., Watson, C. C., A.Z.Ch.E. J . 11, 1051 (1965). Lapidus, Leon, Peterson, T. I., A.I.Ch.E. J . 11, 891 (1965); Chem. Ene. Sci. 21. 655 (1966). Pine: H.. HaaP. W. 0..J.’Am. Chem. SOC.82. 2471 (1960). Soboi, M.,IBM Share Program, No. 155i; NU‘ Scoop Linear Surface Minimization Routine, 1963. Solomon, H. J., D. Eng. - dissertation, Yale University, New Haven, Conn., 1966. Solomon, H. J., Bliss, Harding, Butt, J. B., IND. ENC. CHEM. 6, 325 (1967). FUNDAMENTALS RECEIVED for review March 13, 1967 ACCEPTEDNovember 6, 1967
KINETICS OF COUPLED FIRST-ORDER REACTIONS WITH TIME-DEPENDENT RATE COEFFICIIENTS IN TERNARY SYSTEMS LOTHAR R I E K E R T AND J A M E S WE1 Research Department, Central Research Division,“Mobil Rcsearch and Development Cor)., Princeton, N . J . 08540
A method for calculating the kinetic behavior of coupled first-arder reactions in a ternary system from any given set of time-dependent rate coefficients i s based on Lie algebra. The method reduces the integration of a set of simultaneous differential equations with variable coefficients to successive integrations of ordinary differential equations. The application is illustrated by an example.
HE
kinetics of coupled first-order reactions can be handled
T by linear algebra (as described by Wei and Prater, 1962) if the rate coefficients are: constant. I n many cases of practical importance, however, the rate coefficients are functions of time ; two important examples are reactions proceeding under nonisothermal conditions and catalyst deactivation. T h e problem of optimization of a temperature profile in systems of coupled reactions has been treated by Bilous and Amundson (1956) and by Mah and Aris (1964). Froment and Bischoff (1962) have considered the effect of catalyst deactivation on the selectivity of irreversible para.llel reactions.
A method is described for calculating the behavior (reaction paths and time dependence) of any system of coupled firstorder reactions that have time-dependent rate coefficients and involve three components. A mathematically equivalent problem was encountered by Sargent (1963) in connection with the dynamics of distillation columns; he obtained numerical solutions with the help of power series. T h e present method is based on the theory of Lie algebras (Jacobson, 1962), which allows reduction of the integration of a set of simultaneous linear differential equations with variable coefficients to successive integrations of ordinary differential equations. VOL 7
NO.
1
FEBRUARY 1968
125
