Danckwerts, P. V., Sharma, M. M., Chem. Engr. 44, No. 202, CE244CE280 (October 1966). Eickmeyer, A. G.,‘ Chem. Eng. d6, 113-16 (Aug. 25, 1958). Ellis, J. E . , Trans. Inst. Chem. Engrs. 38, 216-24 (1960). Jeffrevs, G. V., Bull, A. F., Trans. Inst. Chem. Enqrs. 42, T118-25
-~ - (moles COS absorbed per liter) I.(
[C032-10 = solution viscosity, cp.
@OH, @Am
= parameter for absorption with reaction,
Dconkos[OH-I/k2 and D C O ~ [ R R ’ K H ] / ~ Lrespectively ~,
(1964).
Killeffer, D. H., I d . Eng. Chem. 29, 1293 (1937). hlcNeil, K. M., Ph.D. dissertation, University of Cambridge, 196.5.
McNeil, K. M., Danckwerts, P. V., Trans, Inst. Chem. Engrs. 42. T294-7. ( 1 9fi4). __,
SUBSCRIPT
_ _ I -
= initial value, before COZ absorption
0
~
Pvlullowney, J. F., Petrol. Refiner. 36 (12), 149-52 (1957). Nysing, R. A. T. O., Kramers, H., Chem. Eng. Sci8,Sl-9 (1958). Pinsent. E. R., Pearson, L.. Roughton, F. J. W.. Trans. Faraday SOC.62. 1512-20 (1956).’ - ’ Kiou, P., Cartier, P., Compt. Rend. 184, 325-6 (1927); 196,
literature Cited
1727-8 (1928).
Astarita, G., lIarrucri, G , Gioia, F., Chem. Eng. Sci. 19, 95-103 (1964).
Benson, 13. E., Field, J. H., Hayes, W. P., Chem. Eng. Progr.
Sharma, M. &I., Trans. Faraday SOC.61, 681-8 (1965). Sharma, hl. bl., Danckwerts, P. V., Chem. Eng. Sci. 18, 729-35 (1963); Trans. Faraday SOC.69, 386-95 (1963). RECEIVED for review July 22, 1968. ACCEPTED January 10, 1969.
62 110). 433-8 (1956).
Ilanckwe;ts, P. Y., Gillham, A. J., Trans. Inst. Chem. Engrs. 44 (2), T42-T49 (1967).
Danrkwerts, P. Y., IIcNeil, K. >I., Chem. Eng. Sei. 22, 925-30
/.
(1067) ,-. .I.
Danckwerts, P. V., McNeil, K. X,, Trans. Inst. Chem. Engrs. 46 ( I ) , T32-T49 (1067).
Work carried out during 1965-66, while A. L. Shrier was a rstdoctoral visitor in the Department of Chemical Engineering, niversity of Cambridge, on leave from Esso Research and Engineering Coo
ESTIMATING RATE CONSTANTS An, Improvement on the Time-Elimination Procedure NORMAN R. DRAPER AND
HIROMITSU KANEMASU
Department of Statistics, University of Wisconsin, itfadison, was. 53706
RElJl
MEZAKI
Department of Engineering and Applied Science, Yale University, New Haven, Conn. 06520 In certain complex chemical systems, the exact starting time of the reaction i s not known. One method sug. gested for treating such situations is to eliminate reaction time from the differential equatiohs which describe the system, and estimate the parameters in the models which result. The deficiencies o f this method are first illustrated by two constructed examples. An alternative method of tackling the problem is then presented, and i s applied to the same two examples to demonstrate its superior properties.
N RECENT years great progress has been made in the
I development of systematic procedures for determining rate constants in complex reaction systems. The development and application of these procedures have been documented and extensively discussed (hmes, 1960, 1962; Ark, 1964; Ball and Groenweghe, 1966; Himmelblau et al., 1967; Lapidus, 1961; Wei, 1962, 1963, 1965). The modeling process typically involves the following sequence of steps: -4. Postulation of plausible reaction models based upon experimental findings and setting up of rate equations which represent the principal features of the models. B. Estimation of the rate parameters involved in the rate equations, utilizing experimental observations. C. Examination of the adequacy of the fitted models via statistical methods such as examination of residuals and analysis of variance, and physical and chemical interpretation of the parameter estimates.
To provide (for step A ) a mathematical representation of complex chemical systems, it is a common practice to postu-
late differential rate equations which describe the concentration changes of all existing major components with respect to reaction time. The analytical or numerical solution of these equations is usually possible but is often not easy, particularly with a complex system which involves many reaction components and, consequently, many unknown rate parameters. The use of a “steady-state approximation” may considerably simplify the situation. However, in many problems, such an approximation obviously fails and leads to an erroneous interpretation of rate data (Bowen et d.,1963). The elimination of an independent variable, frequently reaction time, has been suggested as an alternative procedure to obtain a simplified rate expression (Benson, 1960). The elimination can usually be performed readily by dividing all but one of a series of differential equations by the remaining equation. This suggestion appears, a t first sight, attractive, particularly when the starting time of the reaction is not known. Such a situation is often encountered, for example, when a recycling catalytic reactor is used for a kinetic experiVOL.
8
NO.
3 AUGUST
1969
423
ment. A major disadvantage of this “time-elimination” procedure is, however, that the parameter estimates which are obtained by using the resulting equations are usually much more highly correlated than the parameter estimates which result from the original equations. As a result, the precise determination of individual rate constants may be difficult, or even impossible. The principal purpose of this paper is to offer, and to illustrate, an alternative estimation procedure applicable when the starting time of a reaction is not known, and to substantiate, via the approximate correlation matrices of the estimated parameters, the advantages of this alternative suggestion. Estimating Parameters of a Nonlinear System
All the standard methods of estimating parameters (in models that are nonlinear in the parameters) have, as their basis for justification, the method of maximum likelihood, applied to a normal error distribution function. Draper and Smith (1966) give additional detail and we summarize, here, only the points that we need in our analysis. In what follows we distinguish three cases of interest. Single Response Case. Suppose that a general nonlinear reaction rate equation is given
r
=
f (21,XZ, . . , Z m ;kl, kz, . ., kp) I
(2.1)
where XI,xp, . . . , xm denote the independent variables and kl, k2, . . . , k, denote rate parameters. In vector notation, Equation 2.1 becomes
r = f(x; k )
+
where xu denotes the value of the independent variable settings which resulted in yu. An experimentally observed value, yu, differs from the true response f(xu;k) because of an experimental error, E,. It is usually assumed that these errors are independently and normally distributed. The best estimates of unknown parameters k can be obtained by minimizing the sum of squares function N
C
(yu -
(2.4)
u=l
with respect to the elements of k, where r, is obtained from 2.2 by substituting xu for x. There are several methods for carrying out the actual minimization, one of the most generally useful being due to Marquardt (Draper and Smith, 1966). We have used Gaushaus, the University of Wisconsin’s programmed version of the Marquardt procedure (Meeter, 1966), actually to obtain our estimates in the examples we have used. Most good minimization procedures also provide an approximate correlation matrix for the estimated parameters. Draper and Smith give details of the derivation of this matrix (1966). This correlation matrix, C, is a p X p matrix whose (i,j ) element can be denoted-by c , ~ . The latter is the correlation between the estimate k , and the estimate k , (so that, of course, c , ~= 1 when i 73 ) . If c,j is close to +1, the parameter estimates k , and k j are highly positively correlated, 424
I&EC
FUNDAMENTALS
If we can assume that N multiple observations yl,, yzu, and (u = 1, 2, . . ., N ) on A , B , and C, respectively, have completely independent errors attached to them, we estimate the parameters of the models by minimizing the weighted sum of squares function y3,
S(k) =
1
”
C (yiu - &I2+ u=l
1
7C JVB u-1
-
(~2u
+
(2.2)
where x = (51,XZ, . . . , zm)’and k = (kl,kz, . . . , k,)’. If the mathematical model of Equation 2.1 or 2.2 is adequate to describe a system of interest and if N experimental observations on r-namely, yl,y2, . . . , yN-are available, we can write u = 1, 2, . . . , N yu = f(xu;k) e,, (2.3)
S(k) =
and changes in both ki and kj in the same direction can be made with little effect on the sum of squares function. This usually implies a great deal of ambiguity in the “best” values of the parameters. If c i j is close to - 1 (high negative correlation), similarly ambiguity exists when the parameters are changed one in the positive and one in the negative direction. When parameters are highly correlated, it is usually difficult to obtain precise point estimates, because the sum of squares contours tend to be long ellipsoidal, or bananalike, shapes in such cases. Thus estimation situations which result in highly correlated parameter estimates are to be avoided if at all possible. It is sometimes feasible to avoid (or a t least reduce) such difficulties by a good choice of experimental design or by reparameteriaation methods (Blakemore and Hoerl, 1963; Box, 1960, 1964; Draper and Smith, 1966). Multiple Response Case. Situations frequently arise where we wish to consider several responses simultaneously. In such a case, Equation 2.2 is replaced by several models. For example, when there are three responses A , B , and C, we would write
1
N
where W A ~W, B ~and , W’c2 are weighting functions for components A , B, and C, respectively, which allow us to make appropriate adjustments when the errors in different responses have different variances. Equation 2.6 replaces 2.4 for situations of the type mentioned (Hunter, 1967; Johnson et al., 1968).
In our applications to examples in subsequent sections we use formula 2.6 with all the W’s equal to unity, since this is the appropriate formulation for the data we have generated. The appropriate correlation matrix for the estimated parameters can be obtained in a straightforward manner which is an extension of the single response procedure. In this extension we regard the N observations on the three responses as 3 N consecutive observations on a single response whose response function changes its form a t the ( N + 1)st and ( 2 N + 1)st observations. This change is easily made on Gaushaus and produces the appropriate C matrix for this situation. If the sets of observations (yl,, yz, ~ 3 and ~ )(ylv,yzU, ~ 3 are independent for u # u but ylu, y~,, y3, are not independent of one another, we have a much more complicated situation. The best parameter estimates are found by minimizing a certain determinant. The Bayesian development of this criterion is given by Box and Draper (1965) and it can also be derived by maximum likelihood methods. [See also Hunter (1967); for an example, see Mezaki and Butt (1968).] We, however, do not need this additional complication in our subsequent discussion. We present two examples, each of which is solved via the time elimination procedure and then resolved using the methods of this section.
~
)
Analysis of Kinetic Data, Starting Time of Reaction Unknown
Example 1. Let A , B , and C represent the values of the concentrations a t time t for the reaction system k?
ki
,4+B--tC
(3.1)
The concentration changes can be accounted for by the following ordinary differential equations.
-dA/dt = klA dB/dt
=
(3.2)
klA - k2B
(3.3 1
kzB
(3.4)
dC/dt
=
We assume that experiments were conducted in a batch stirred tank reactor in which a certain degree of mixing is needed to initiate the reaction. I t is impossible to determine exactly when this degree of mixing has been attained. What does this imply? It means that we must measure our time scale from some arbitrary starting point and that the exact elapsed time of reaction is not known. Table I shows a set of 18 simulated observations from the chemical reactions system 3.1. The time is assumed measured from an arbitrary origin. (How these data were actually generated is discussed below.) The problem is to estimate parameters k l and k2 using the data of Table I. We first discuss the time-elimination solution and then show how to apply, and the superiority of, the “extra parameter” method which we recommend. TIME-ELIMINATION SOLUTION.The time-elimination procedure is applied as follows. We eliminate reaction time by dividing Equations 3.3 and 3.4 by Equation 3.2. The division yields the following two differential equations
d B - klA - k2B -dA
Now, strictly speaking, it is incorrect to apply the least squares technique as stated above, via minimization of Equation 2.6, to models 3.7 and 3.8 since, for the A’s on the right-hand side, we have to substitute observations subject to error. However we minimize 2.6 even though it is only an approximation to the correct least squares treatment of such situations (which is extremely complicated). This is one, but not the major, drawback of the time-elimination system. The application of the standard least squares technique, employing the data in Table I, gives rise to the estimates = 0.2564
A
=
{exp C-kz(t+
Table 1.
(Akzikl- 1)
l)]
(3.8)
Simulated Experimental Data for Example 1 Time, -Win.
A
1 2 3 4 5 6 7 8 9
6 9 12 15 18 21 24 27 30 33 36 39 42 45 48 51
0.657 0.532 0.433 0.349 0.281 0.229 0.187 0.151 0.114 0.095 0.080 0.059 0.054 0.041 0.032 0.030 0.023 0.021
11 12 13 14 15 16 17 18
54
57
+
(3.10)
T)]
r)I- exp C--kl(t+
r)1)
(3.11) (3.12)
Via minimization of 2.6, using the data in Table I, we obtain estimates h1
= 7.0181
X lo-’
h2 =
1.9995 X
i = 3.0142
(3.13)
Run No.
10
- ( A-
exp [ - k l ( t
C=1-A-B (3.7)
:
(3.9)
“EXTRA PARAMETER’’ METHOD.Suppose, in system 3.1 we introduce an additional parameter T , which represents the unknown time lag between the actual starting time and the chosen arbitrary origin. Solving 3.2 to 3.4 analytically with this condition in mind, we obtain the solutions:
which can be solved analytically. The results are
C= - (k2 k x :
=
of k l and kz, respectively. Examination of the correlation matrix, C, shows, however, that- the correlation between kl and f , is 1-i.e., that kl and are perfectly correlated! This implies that estimation of kl and kz separately is impossible. That this is so is clear upon examination of Equations 3.7 and 3.8, both of which can be rewritten in terms of a single parameter, D = kz/kl. In other words, use of the time-elimination procedure in this example severely curtails the information provided by the data on the two parameters, making it impossible to estimate them separately. We can overcome these difficulties by using an extra parameter to designate the unknown starting time lag, and then using standard least squares procedures.
(3.5)
kiA
I& 0.0731
A, B, and C are mole fractions.
B 0.327 0.424
0.592 0.584 0.567 0.549 0.525 0.513 0.490 0.460 0.440 0.427
C
0.025 0.042 0.076 0.107 0.133 0.167 0.211 0.244 0.276 0.314 0.350 0.384 0.422 0.447 0.478 0.505 0.532 0.561
The minimum value of the sum of squares function was 1.0166 X loF5.
We now investigate how good these estimates are and how they compare with those in 3.9. The data in Table I were in fact generated by using Equations 3.10 to 3.12 with kl = 0.07, k2 = 0.02, and r = 3 and adding independent random normal deviates, all with standard deviation 0.003. Clearly 3.13 are very close to the true values, as we should hope; values 3.9 obtained through the time-elimination procedure are extremely poor estimates. In fact, the estimates obtained in this situation tend to have a strong dependence on the starting values employed, since a great many pairs of values ( k l , k z ) achieve practically the same sum of squares value. In the time elimination procedure, not only are the estimates poor but we do not have an estimate for r which the extra parameter procedure allows us to obtain. The correlation matrix for the present parameter estimates can be evaluated by the usual methods and we find correlation VOL.
8
NO.
3 AUGUST
1969
425
the equations
values as follows: Between til,
rl,
0.256
kl, i k.2,
(3.16)
Correlation -0.843 - 0.396
i
(3.14)
Clearly, the situation is much better than in the time eliminat,ion procedure, t’he only fairly high correlation being between k1 and i . A further demonstration of the superiority of the extra parameter procedure is given by an examination of the approximate 95y0 confidence regions for the parameters (Draper and Smith, 1966) which arise in the two cases. These are shown in Figures 1 and 2. Clearly the extra parameter procedure is much superior. A second example, while less striking, also shows the superiority of the extra parameter method. Example 2. We now consider ai1 example involving the system.
dB _ - k l A - knB + k3C
(3.17)
dC _ - kzB - k3C
(3.18)
at
at
Again we assume that the starting time of reaction is not known and that it is desired to estimate the k’s. We use the simulated data shown in Table 11. TIMEELIMINATION PROCEDURE. Again the independent variable reaction time is eliminated by dividing Equations 3.17 and 3.18 by Equation 3.16. The resulting equations are: dB k l A - kzB + k3C _ -
(3.19)
dC - kzB - k3C _ -
(3.20)
dA
-klA
-klA
dA (3.15)
The coilcentration changes of A , B , and C are explained by
Equations 3.19 and 3.20 have the analytical solution
B=
ki - k3 { A (kz+ka)lki ki - (kz k3)
+
POINT ESTIMATE
- AI +
C = 1- A - B
FROM TIME ELIMINATING PROCEDURE
which leads, if we adopt a minimization procedure similar to that used in the first example, to estimates 9 5 % CONFIDENCE REGION
i&=
FROM TIME ELIMINATING
0.6122
/& =
f , = 0.1322
0.4123
(3.22)
FROM EXTRA PARAMETER
The correlations between the parameter estimates are given in Table 111. EXTRA PARAMETER METHOD. Introducing, as before, a parameter to take account of the initial time lag, we can write the analytical solutions of Equations 3.16 to 3.18 as:
(.07, .02)
.M
0.10
A = exp [-kl(t
0.15 R A T E CONSTANT, hI
Figure 1. Approximate 95y0 confidence regions from the two estimation procedures
B
=
4 exp C-kl(t -k 711 4G3
exp C- (kz
+
(3.23)
T)]
+ k3) + T)I 4-GZ (t
C=l-A-B
(3.24) (3.25)
2.01EXTRA PARAMETER PROCEDURE 9 5 % CONFIDENCE REGIOI
FROM EXTRA PARAMETER PROCEDURE
FROM T I M E ELIMINATING PROCEDURE I.--
7.00
6.95
7.10
7.05 RATE CONSTANT, k, a
IO‘
Figure 2. Approximate 95% confidence regions from the two estimation procedures, situation near true values 426
l&EC
FUNDAMENTALS
Table II. Simulated Experimental Data for Example 2 Run No. Time, Min.
1 2 3 4 5 6 7 8 9 10 11 12 13 14
4 7 10 13 16 19 22 25 28 31 34 37 40 43
A
B
C
0.772 0.602 0.498 0.404 0.342 0.254
0.215 0.349 0.424 0.458 0.476 0.508 0.488 0.490 0.506 0.475 0.485 0.462 0.437 0.445
0.143 0.192 0.223 0.281 0.329 0.359 0.399 0.423 0.450 0.492 0.509
A , B , and C are mole fractions.
Acknowledgment Table 111.
Correlations of Parameter Estimates
A grant of computer time by the Wisconsin Alumni Research Foundation through the University Research Committee is gratefully acknowledged. The authors are grateful t o a referee for his suggestions for a certain improvement in this paper.
Correlation in Time-
Extra
elimination
parameter method
Between A
!I,
method
I
0.337 -0.292 0.954 -0.812 -0.309 -0.237
-0.944 -0.916 0.922
k?
$1, i 8
k?, k3 $1, . i $2, . i ka, .i
...
... ...
Nomenclature
where
- Jc2)/(Jc2 4-Jc3 - k i ) G2 = k3/ (Jcz -I- k 3 ) G3 -4 - G2 GI =
+=
2.970
(Y3u,
C U )
1
C, respectively
.
models (3.26)
= vector of parameters = estimate of parameter k
f
N
= number of observations
P
= number of parameters
r
= response
= sum of squares of deviations for a set of
(k)
parameter values k = reaction time
1
W A ~W, B ~Wc2 , = weighting functions for A , B , and C, respectively X
6.9860 X
= 2.1552
Bu
k
k2 = 4.1282 X
k3
(YZU,
of A , respectively = observed and calculated mole fraction of B , respectively = observed and calculated mole fraction of
= observed and calculated mole fraction
= element of correlation matrix kl, ICq, k ~, , ., k, = rate constants or “parameters” of the
(In fact the 14 data points in Table I1 were obtained by using Equations 3.23 to 3.26, setting kl = 0.07, kz = 0.04, k( = 0.02, T = 3, and adding, to the calculated response at each time, a random normal deviate with standard deviation 0.01.) Employing the data in Table I1 we obtained the following least squares estimates. =
Au)
Cij
(Jci
kl
(Ylu,,
= vector of independent variables
. . ., zm
~ 1 ~ x 2 ,
= independent variables
?Ju
= observed value of uth experiment
GREEK LETTERS
X (3.27)
The minimum value of the sum of squares function was 3.2621 X
= experimental error of uth experiment
EU
= time lag in reaction time
7
literature Cited
We again see excellent agreement of estimates 3.27 with the true values, as we might expect. Estimates 3.22 are, again, extremely poor and unreliable. I n the extra parameter method, we also have a n estimate of the time lag, T , which the time-elimination method cannot provide. -1s a final comparison we examine the two approximate correlation matrices. The various terms are shown side by side in Table 111. hlthough in one case (&,k3) we see a slightly increased correlation, the over-all picture clearly favors the extra parameter method and the substantial increases in correlations for the time-elimination method are clearly noticeable. The correlation arrangement in the final column of Table I11 hints to us that Equations 3.23 to 3.26 are overparameterized. Examination of 3.21 confirms this, since B (and hence C) can be represented by a model which involves only two parameters, 01 = k2/kl and 02 = k3/kl. Thus, again, the time elimination procedure makes it impossible to estimate all the original parameters.
Ames, W. F., Znd. Eng. Chem. 62, 517 (1960). Ames, W. F., IND. ENG.CHEM.F U N D A M E S T A L S 1, 214 (1962). Ark, R., IND.EKG.CHEM.FUSDAMENTALS 3, 28 (1964). Ball, W. E., Groenweghe, L. C. D., IND.ENG.CHEM.FUNDA-
Conclusions
WiC’J., Prater, C. D., Aduan. Catalysis 13, 204 (1962). Wei, J., Prater, C. D., A.1.Ch.E. J . 9, 77 (1963). Wei, J., IND.ENG.CHEM.FUNDAMENTALS 4, 161 (1965).
MENTALS
6, 181 (1966).
Benson, S. W., “Foundations of Chemical Kinetics,” McGrawHill, New York, 1960. Blakemore, J. W., Hoerl, A. E., Chem. Eng. Progr. Symp. Ser. 69, 14 (1963).
Bowen, J. R., Acrivos, A., Oppenheim, A. K., Cheni. Eng. Sci. 18, 177 (1963).
Box. G. E. P.. Ann. N . Y . Acad. Sci. 86. 792 (1960). Box; G. E. P., Dept. of Statistics, Univerdty of Wisconsin, Madison, Wis., Tech. Report 26, (1964). Box, G. E. P., Draper, N. R., Biometrika 62, 355 (1965). Draper, N. R., Smith, H., “Applied Regression Analysis,” Chap. 10, Wiley, New York, 1966. ENG. Himmelblau. D. hl.. Jones. C. R.. Bischoff. K. B.. IND. CHEM.FUNDAMENTALS 6, 539 (1967). Hunter, W. G., IND.ENG.CHEM.FUNDAMEFTALS 6, 461 (1967). Johnson, R. A., Standal, N. A,, Mezaki, R., ISD.EKG.CHEM. FUNDAMEXTALS
7, 181 (1968).
Lapidus, Leon, Chem. En9. Progr. Symp. Ser. 67, No. 36, (1961). Meeter, D. A,, “Nonlinear Least Square (Gaushaus),” “University of Wisconsin Computing Center Users’ Manual,” Vol. 4, Section 3.22, 1966. Mezaki, R . , Butt, J. B., IND.EXG.CHEM.FUNDAMENTALS 7, 120 I I U R R ) \_.._,.
Through two multiresponse examples, we have illustrated the application of methods of treating kinetic parameter estimation when the time at which the reaction starts is unknown. Our suggested extra parameter method appears to have superior characteristics and is recommended.
RECEIVED for review September 18, 1968 ACCEPTED February 24, 1969 Financial support received from the National Science Foundation under Grant GK-105.5, and from the Office of Naval Research under Contract 1202(17) Project 042 222.
VOL.
8
NO.
3
AUGUST
1969
427
