On the Estimation of Rate Constants for Complex Kinetic Models Himmelblau, et a/.,recently described a method for the estimation of rate constants from a large number of data points along the concentration trajectories. By using a natural cubic spline curve-fitting technique, the applicability of their method has been extended to cases where only a small number of data points are available.
I n a recent paper Himmelblau and coworkers (1967) described a n interesting method for estimating rate constants of complex kinetic models from isothermal, batch or plug flow reactor data. Their method is applicable when a large number of data points along the concentration trajectories are available. I n the present work we extend their technique to cases where the number of data points is small. The development of the method and the results from a sample problem will be discussed in this note.
second derivatives in the interval [to, t g ] ,with to < tl which minimizes the functional
The Basic Method
where ( t - tk)+ = t - t k if t 2 t k , 0 otherwise. The q coefficients are uniquely determined by solving the 9 linear equations
The mathematical model considered in this study is given by the following differential equations in which t is the time or a time-like variable
It was shown by Schoenberg that the minimizing function is a natural cubic spline having to, . . . , t , as its knots
2
ak+2
=
k=O
The terms R j are functions of concentrations. Let T = t,, . . , , t g } be a set of q 1 distinct time points, with to being the initial time; eq 1 can be integrated to yield
+
ifo,
< . . . O), the spline approaches a straight line. For each Rj we can thus fit a natural cubic spline. Equation 3 can then be integrated analytically. Alternatively, we can fit each C i by a natural cubic spline. This approach may be needed if some of the Cik are missing. Equation 3 can then be integrated by a suitable numerical quadrature scheme, for example, a Gauss quadrature formula.
+
where (3)
If C f k are known, and if Elk can be estimated, eq 2 becomes a system of linear equations in the rate constants b,. The system can be solved by some discrete minimization techniques, for example, a least-squares method, provided that the number of independent equations is not less than the number of rate constants. Estimation of
Numerical Examples
The following example is taken from Himmelblau, et al.
R,r,
When a large number of data points are available, eq 3 can be integrated by numerical quadrature, using the data points. The scheme, however, is no longer applicable when the number of data points is small. A curve-fitting procedure followed by integrating under the curves thus obtained could be used. Note t h a t the concentration trajectories are usually transcendental functions of time which, over a relatively large interval, are difficult to approximate by simple low degree polynomials. , One class of curve-fitting methods which has attracted wide attention in recent years is the spline function approximation. For the present work, we make use of a result due to Schoenberg (see Greville, 1968). Given a set of observed punction values {yo, yl, . . . , yp} over T,we seek a function f ( t ) having continuous derivatives and piecewise continuous
Initial concentrations a t to = 0 are C1 = 1 and CZ = C3 = 0. The differential equation for C3 is omitted since it is linearly dependent. Two sets of time points were studied: (I) Five-point :
( 2 ) Six-point:
T T
=
( 0 , 0.2, 0.5, 1.0, 2.0)
=
(0, 0.1, 0.2, 0.5, 1.0, 2.0)
No data smoothing was attempted (Le., g = 0), and the spline approximation for C3 was obtained by material balance. The linear system eq 2 was solved by a least-squares technique using equal weights. All computations were carried out Ind. Eng. Chem. Fundam., Vol. 10, No. 2, 1971
321
Table
I. Rate Constants Estimated from Error-Free Data Time Points
bl
31 5 6
1 0 1 013 0 994 1 000
True H-J-R This uork This work
b2
0 0 0 0
5 497 485 503
br
b3
10 10 10 9
0 125 620 855
5 4 5 4
0 990 420 962
Table II. Simulations With Random Errors
True values Five-point estimates Average Range Standard deviation Six-point estimates .Sverage Range Standard deviation
bi
bp
b3
1. o
0.5
10.0
bd
5.0
10.447 0.431 0.968 5.386 7,889 3,779 0.246 0.875 -1.020 -16.673 -9,190 -0.547 0.055
0.119
3.311
2.003
0,460 9.641 0.978 4.358 0,282 7.397 0.887 3,483 ~ 1 1 . 7 6 7-6.210 -1.047 -0.582
I n practice one should always be cautious in interpreting the estimated rate constants since observed concentrations are seldom error-free. For the present method, as well as the original method of Himmelblau, et ai., it is extremely difficult, if not impossible, to make a rigorous statistical analysis. A practical approach is to carry out a Nonte Carlo simulation study. The results from such a study using a sample size of 9 are given in Table I1 where errors in Clr and C2k were assumed to be independent, normally distributed with means of zero, and standard deviations amount to 2.5% of the respective true concentrations. (The covariance terms are omitted from the table.) It can be seen that the estimated standard deviations for the 6-point case are reasonably small. The proposed procedure should, therefore, give good results. At the least, the estimates provide a good starting point for a much more costly nonlinear regression approach (category 3a in Himmelblau , et al.). Acknowledgment
The author is indebted to the management of the Engineering Department, Union Carbide Corporation, Chemicals and Plastics, for permission to publish this paper. Nomenclature
0,059
0.139
1.826
1.264
on a digital computer with double precision arithmetic in solving the matrix equations. Results from error-free simulations are summarized in Table 1. The rate constants obtained by Himmelblau, et al., are also included for comparison. With the added time point t = 0.1, the 6-point case yields significantly better estimates. This indicates, as noted by Himmelblau, et al., that the estimates are sensitive to errors resulting from approximating the initial integrals. The use of 5 and 6 time points are for illustration purpose. For the present problem there are two independent observations a t each time point. The total number of observations excluding the initial time are then 8 for the 5-point case and 10 for the 6-point case. Since there are 4 rate constants to be estimated, the corresponding degrees of freedom are 4 and 6, respectively.
b, = rate constant for reaction 1. C I = concentration of component i. g = a smoothing factor in spline curve-fitting. m = number of components. n = number of reactions. s t j = stoichiometric coefficient for component i in reaction j . t = time. wI = weighting factor in spline curve-fitting. literature Cited Greville, T. ?T. E., Introduction to Spline Functions, “Theory and Applications of Spline Functions,” p 1, Academic Press,
1969. Himmelblau, D. AI., Jones, C. R., Bischoff, K. B., IND.ENG. CHEM.,F U N D A M . 6 , 539 (1967).
Y. P. TANG Engineering Deparfment Cnion Carbide Corporation Chemicals and Plastics South Charleston, K-.V a . $5303 RECEIVED for review March 2, 1970 ACCEPTED November 30, 1970
CORRESPONDENCE
Determination of Start-up Conditions for Chemical Reactor Stability SIR: I wish to raise an important consideration which appears to have been overlooked by Han (1970) in the determination of start-up conditions for chemical reactor stability. Han considers the start-up of a continuous stirred-tank reactor in which a first-order irreversible reaction A + I3 occurs. There are three possible steady states and Han has chosen the low conversion steady state to be the desired state of operation, since the middle steady state is unstable and the high conversion steady state is beyond the safe temperature limit of the reactor. 322
Ind. Eng. Chern. Fundam., Vol. 10, No. 2, 1971
The middle steady state, however, can be made stable very easily by the use of proportional control with a sufficiently large proportional control constant. This mode of control as applied to a reactor is discussed and analyzed in detail first by Ark and Amundson (1958), and then by Lapidus and Luus (1967) , mho considered the stability problems associated with a series of CSTR’s. If we look a t Han’s proposed conditions of operation, it is readily seen that his chosen steady state of CAS = 0.484 with CAO = 0.50 gives a yield of only 3.2% whereas the “un-
