Ind. Eng. Chem. Prod. Res. Dev. 1980, 79, 430-434
430
Parameter Estimation in Linear Chemical Reaction Systems Erhard G. Christoffel Lehrstuhl fur Technische Chemie, Ruhr-Universtat Bochum, 4630 Bochum 1, West Germany
The kinetic analysis of simultaneous isomerization and cracking of n- and isohexanes on a bifunctional Pt/A1203 catalyst, a reaction system which has been shown to follow pseudo-mass-action kinetics, has been accomplished with different parameter estimation methods. It is demonstrated that the approach of Wei-Prater, a method described by Gavalas, and a library computer program which contains parameter estimation and statistical testing procedures are equally well suited for the kinetic analysis of the complex linear reaction system.
Introduction Reaction rates in heterogeneous catalysis are usually formulated in the form of the law of mass action or in the form of a hyperbolic model, where the denominator expresses the competition of the reacting species for adsorption sites. The hyperbolic model decomposes into a pseudo-mass-action system in cases where the denominator is the same function of composition and time in the rate equations of all reacting species and can thus be factorized. Chemical reactor modeling with the goals to scale up from laboratory reactors to pilot plants and industrial reactors, to select optimal operation conditions for different feedstocks and catalysts, or to contribute to the understanding of basic reaction mechanisms consists in principle of two parts: constructing of a suitable model and estimation of the parameters of the model. In order to keep the interpretation of kinetic experiments of heterogeneous catalytic reactions as simple as possible, reaction conditions in a first step are usually selected in such a way that the measured conversions are largely determined by the rates of the chemical reactions; that means the reactor model is identical with one of the above-mentioned reaction rates models. Many chemical rate processes can be satisfactorily represented or at least approximated by first-order reactions. The estimation of model parameters of the corresponding rate equations particularly for the kinetic analysis of complex reactions is nowadays mostly done by use of computer packages which contain parameter estimation and statistical testing procedures (Froment, 1975). For the kinetic analysis of linear reaction networks other methods have been described in the literature, which however have not been widely used in experimental investigations. In a previous paper one of these methods, the approach of Wei and Prater (1962), has been used for the kinetic analysis of a six-component system (interconversion and cracking of the five hexane isomers) (Christoffel, 1979). In the present investigation the results obtained from this analysis are compared with results from a kinetic analysis of the same reaction system by use of the method described by Gavalas (1973) and by use of a library computer program (Bard, 1967). Theoretical The governing differential equations for the time behavior of simultaneous interconversion and cracking of the five hexane isomers are iu = KW
where
(1)
r-eI 0,-methylpenwe
= a3-methyipentane
a,,3 -dimethylbutane % , z -dimethylbutane
is the composition vector and K is the rate constant matrix. The assumption of first-order reactions is justified in the case of constant hydrogen partial pressure and high hydrogen/hydrocarbon partial pressure ratios (Surjo and Christoffel) (1978). The Wei-Prater method elaborates via a transformation of the pure component coordinates in the composition space into the characteristic directions of the rate constant matrix. These eigenvectors must be experimentally located, which is easily achieved in the case of three-component systems and much more complicated in the case of large systems of first-order reactions. The time behavior of the transformed system is described by uncoupled first-order differential equations. The rate constants for the reactions of the hypothetical species, which are defined along the characteristic directions, are determined from a parametric representation of any curved reaction path in the composition space which contains sufficient amounts of all hypothetical species. The rate constant matrix of the original component system can then be computed from the relation K = XAX-', where X and X-' are the matrix of the eigenvectors and its inverse and A is the diagonal matrix of the eigenvalues A; of the rate constant matrix. In the method of Wei and Prater the problem of estimating a set of parameters from experimental data decomposes into an independent estimation of each eigenvector and its corresponding eigenvalue. A similar approach has been described by Gavalas (19731, but instead of determining the eigenvectors by a sequence of experiments the decomposition is achieved by computation. The principal ideas of the method are the following. The solution to eq 1 is a = Xe-AtX-'a(0)
(2) where a(0) is the initial composition at t = 0. After some simple rearrangements the following expression is obtained
? r ~ ( a ~ k-,a~l ( )l , ~ ) e - * j =~ lo
(3)
1=1
where rG) is the j t h row of the matrix X-', Iz = 1, ..., N ( N = number of measurement times), 0 = 1, ..., m ( m = number of initial states), 1 = 1, ..., n ( n = number of components in the state vector). Due to Gavalas this equation suggests the formulation of a sum of error squares J
0196-4321/80/1219-0430501.00/00 1980 American Chemical Society
N
m
11,
J = C C 3 1 L r M k , P ) - aL1,P)e-Xtk)12 k = l p = l l-:1
(4)
1
1 = 435°C p = 10,5 bar mole ratio H2ihydrocarbon = 20
WhI
where A is now a variable. Equation 4 can be rewritten as a quadratic form ,I = yTW(A)y,where W(X) is a symmetric matrix. For any given X the minimum p(A) of the quadratic form with respect to X and y is the smallest eigenvalue of W(X), which can be proved by YT 7 p ( ~ =) minl.+l +"(X)y = miny+o-W(X)-IYI IYI introducing a new variable y = yo + t6, where 6 E R" and t E R , a function
is obtained which can be differentiated and which attains a minumum at t = 0. Differentiating the function f(t) and rearranging yields GT(W(X)y0- p(A)yo) = 0 for all 6 E R" or W(A)yo- p(A)y, = 0. The function p(X) attains minimum values at points X1, ..., A, which are the eigenvalues of the matrix K. The corresponding vectors ?(A1), ..., y (A,) are the rows of the matrix X-l. The procedure using the method of Gavalas is the following. First starting from a number of initial states-the number should correspond to the number of components in the state vector-which must span the space, reaction paths in the composition space are measured. From these composition-time da.ta the matrix W(A) is computed for a given A. With any standard routine-in the present study the method of Householder and Givens was usedthe smallest eigenvalue and the corresponding eigenvector of W(X) are computed. Then A is varied over a special range and for each A the smallest eigenvalue of W(X) is again calculated. A plot of the smallest eigenvalues of W(X) against X yields a graph which attains local minima at AI, ..., A,. In the method of Gavalas the problem of estimating a set of parameters from experimental data thus decomposes into a one-dimensional search for the eigenvalues of the rate constant matrix. Standard routines included in computeir packages for parameter estimation are extensively described in the books of Bard (1974) and Himmelblau (1970). Experimental methods (a) Materials. I't/A1203 catalyst (0.35 wt 5%) was provided from Kali-Chemie Engelhard. The pure grade reactants n-hexane, 2-methylpentane, 3-methylpentane, 2,3-dimethylbutane, and 2,2-dimethylbutane were further purified by drying and distillation. Hydrogen was dried over molecular sieves. (b) Apparatus and Operating Conditions. The experiments were carriied out in an isothermally operated fixed bed reactor (i.d 30 mm, length 150 mm); 2-16 g of catalyst (particle diameter 0.4-0.5 mm) was packed between two layers of inert material. At the end of the reaction zone a small product gas flow was withdrawn and fed directly into a gas chromatograph with a squalane capillary column for analysis. The kinetic experiments were performed a t 4:35 "C, a total pressure of 10.5 bar, a hydrogen/hydrocarbon molar ratio of 20, and space velocities of 6-60 h-l. Detailed operation conditions to achieve a constant activity and selectivity level of the catalyst have been dlescribed earlier (Christoffel, 1979). Experifnental Results and Discussion The experimental results for determining the eigenvectors and eigenvalues which are required to compute the rate constant matrix using the method of Wei and Prater
f
h -
Figure 1. Evaluation of the eigenvalues of the rate constant matrix for simultaneous isomerization and cracking of the five hexane isomers with the method of Gavalas.
are given elsewhere (Christoffel, 1979). Composition-reaction time data, which are needed for the Gavalas method and the library parameter estimation program, were measured starting from the following five initial compo-
0.785 0.032 0.179
cy2=
0 0 0
1 0
cy3=
0.4 0 0.3 0 0.3
LYa
0.1 0.1 = 0.7 0
0.1
0.4192 0.4096 0.4341 0.4625 -0.249 -0.0396 -0.0271 0.5483 X - ' = 0.4626 -0.3208 -0.4482 0.3711 0.0179 -0.6305 0.7118 -0.1855 -0.0081 0.0705 -0.053 -0.8965
1 cy5
0 = 0 0 0
0.5041 0.7969 0.591 0.2472 0.4339
(5)
432
Ind. Eng. Chem. Prod. Res. Dev., Vol. 19, No. 3, 1980 reactant
alol=[01/03/07/0/011
lelCfant
-101 : 101101107101311
TzL35’C
cr
p
p i 10 5 bar
105 bar
mole
rol-0 H21 h y d m o i b o n
:2 0
0.6 0.2
0.5
20
10
30 WlF I g h l r n ~ l e ~ ~
10
20
30
WIF I g hlmo;?
Figure 2. Comparison of experimental and calculated (solid lines) reaction paths; rate constants determined with the Wei-Prater method.
Figure 3. Comparison of experimental and calculated (solid lines) reaction paths; rate constants determined with the method of Gavalas. reactant
=IO1
-
101101107101011
T = L35’t p i 105 bor mole rot10 H21hydrocarbon i 20
reaction rate constants, assuming that the catalyst is a constant entity in the total composition region of products and reactants. The rate constants obtained from the kinetic analysis with the Wei-Prater method were used as initial guesses. The results of the parameter estimation are given in Table I together with the standard deviations of the parameters, which were calculated according to Bard (1967). The correlation coefficients between the individual paramaters obtained from the covariance matrix are close to zero, indicating that the parameters are not correlated with each other. A comparison of the results obtained from the three parameter estimation methods can be accomplished on the basis of the experimental expense, of the goodness of fits, and of the contribution to the understanding of the reaction mechanism. Each of these aspects is rather complex, and only some basic problems will be discussed. The Wei-Prater method connects the experimental part and the parameter estimation part of the problem. The rate constant matrix is computed from the equation K = XAX-’. Each eigenvector of the matrix X and each eigenvalue of the diagonal matrix A is determined by fitting experimental composition-reaction time data by straight lines. The search for the straight line reaction paths consists of an experimentally performed iteration procedure, which may be assumed to be tedious, especially if several straight line reaction paths must be located. This problem, however, did not prove to be really severe, because especially on reaction hyperplanes, where the straight lines are connected with large A’s, only a few trajectories are needed to determine the required initial compositions of the straight line reaction paths. Using the Wei-Prater method the experimental difficulties do not originate from the search for the straight line reaction paths. The main limitation arises from the possibility of keeping a constant activity and selectivity level of the catalyst which is the same for all measurements. It will be nearly impossible to locate the straight lines if the selectivity is changing during the experiments. The benefit of the method is the optimal design by which the experiments are performed. Using the library parameter estimation computer program, the experimental part and parameter estimation part may be separated and realized independently of each other. As Froment (1975) discussed, however, no fitting technique will compensate for a poor experimental design. This implies that sequential methods for optimal design of experiments, which take advantage from the information
0.7
0,2
A-
t 0.1 22-OM0
0
10
30
20
W I F I 9 himole1
Figure 4. Comparison of experimental and calculated (solid lines) reaction paths; rate constants determined with a library computer program. Initial composition: a(0) = [n-Cs:0.1/2-MP:0.1/3-
MP:0.7/2.3-DMB:O/2,2-DMB:O.l]. reactant mae
flO.1
o(IO1
=i04/0/03/0/031
T = L35’C p = 105 bar rno e ratio H 2 / hydrocarbon = 20
0-
CL
02
33
02
10
25
30
LO
WIF 1 0 h l m a l e 50 l
Figure 5. Comparison of experimental and calculated (solid lines) reaction paths; rate constants determined with the Wei-Prater method.
of the previous experiments, are favorable in performing the experiments. With or without experimental designs the experimental expense will be smaller than with the Wei-Prater method. The same holds for the Gavalas method, where the experimental iteration procedure is replaced by a mathematical search for local minima of the function p ( X ) . With regard to the agreement of calculated and experimental reaction paths, all three sets of parameters satisfactorily fit the trajectories as Figures 2-7 reveal, although
Ind. Eng. Chem. Prod. Res. Dev., Vol. 19, No. 3, 1980 433
Table I. Relative Rate Constants for Simultaneous Isomerization and Cracking of the Five Hexane Isomers Estimated by Use of the Gavalas Method, the Method of Wei-Prater, and a L i b r n y Parameter Estimation Program (Bard, 1967) library comp. prog. a
Gavalas Wei-Prater n-C, -t 3-MP n-C, 2-MP C,-C, n-C, 3-MP -+ n-C, 3-MP + 2-MP 3-MP + 2,3-DMB 3-MP 2,2-DMB 3-MP + C,-C, 2-MP n-C, 2-MP -+ 3-MP 2-MP 2,3-DMB 2-MP -+ 2,2-DMB 2-MP C,-C, 2,3DMB -t 3-MP 2,3DMB 2-MP 2.3DMB 2,2DMB 2;3DMB C,-C, 2.2DMB --L 3-MP 212DMB + 2-MP 2,2DMB 2,3DMB 2,2DMB -+ C,-C, + -+
-+
-+
-+
-+
-+
--f
-+
--f
a
0.61 0.85 0.97 0.65 2.93 -0.15 0.37 0.82 1 2.2 0.54 - 0.1 1.16 -0.1 3.43 6.2 1.1 0.69 0.21 4.26 0.56
1.16 1.52 0.87 0.95 3.15 0.05 0.06 1.4 1 2.53 0.54 0.31 0.89 0.23 2.84 6.06 0.68 0.15 0.93 3.48 0.87
0.96 (0.1) 1.53 (0.09) 0.52 (0.11) 0.78 (0.08) 3.09 (0.08) O(0.13) 0.23 (0.07) 1.13 (0.08) l(0.06) 2.49 (0.06) 0.75 (0.05) O(0.2) 1.03 (0.07) O(0.13) 3.91 (0.13) 6.83 (0.14) 1.1( 0 . 1 5 ) ’ 0.55 (0.17) O(0.2) ’ 3.92 (0.08) 1.12 (0.19)
The numbers in parentheses are standard deviations.
0
20
30
50
40 W I F Ip.hImolfll
Figure 6. Comparison of experimental and calculated (solid lines) reaction paths; rate constants determined with the method of Gavalas.
single values of corresponding parameters are not identical (Table I). With the Wei-Prater method and the library computer program the thermodynamic equilibrium constants are incorporated in the reaction rate model. With the exception of the rate constants of the slow reactions by which 3-methylpentane is transformed to 2,3-dimethylbutane and 3- and 2-methylpentane to 2,2-dimethylbutane, the two sets of rate constants obtained from the Wei-Prater method and the library computer program are very similar. With the Gavalas method three of the rate constants have small negative values and the ratios of forward and back rate constants of the reversible reactions differ more or less from the thermodynamic equilibrium constants. A discussion of the physical significance of the obtained rate constants presupposes knowledge of the uniqueness of the rate constants and of the physical relevance of the reaction rate model. With the methods of Wei and Prater and Gavalas the rate coefficient matrix is calculated by K = XAX-’. The eigenvectors of the matrix X are linearly independent, which means that K is unique. Furthermore, the system of differential equations iu = K a , by which the
WIF
Ig
h mole
Figure 7. Comparison of experimental and calculated (solid lines) reaction paths; rate constants ‘determined with a library computer program. Initial composition: a ( 0 ) = [n-C6:0.4/2-MP:0/3-
MP:O.3/2,3-DMB:O/2,2DMB:0.3].
kinetics is described, will yield a unique solution for K if the set of initial composition vectors is linearly independent, which is the case in the present investigation. The uniqueness, however, is affected by experimental and numerical error. The calculation of confidence intervals for the parameters is rather complicated. With the Wei-Prater method the eigenvectors and eigenvalues needed for calculation of the rate constant matrix are determined independently from each other by fitting the experimental data by straight lines; i.e., the goodness of these fits is an indirect measure for the reliability of the rate constants. In all cases the linear regression coefficients were 0.999. The reliability of the parameter estimation in case of the Gavalas method has been extensively discussed in the paper of Gavidas (1973). Here the sharpness of the local minima of the function k ( X ) is a measure for the reliability of the eigenvalues. The computation of the smallest eigenvalue K(X) of W(X) entails considerable numerical error. Thus with the Gavalas method the uniqueness of the rate constants is not only affected by experimental error but also by numerical error. The library parameter estimation program provides approximately computed standard devialions of the parameters (Bard, 1967) which contain experimental, numerical, and model error, and some other statistical characteristics. As already stated, the entries of the covariance matrix are close to zero, indicating that the parameters are not correlated. From the point of uniqueness, the rate constants obtained from the Wei-Prater method are most reliable. The last question to be discussed is the physical significance of the reaction rate model. In the individual single rate constants of the interconversion of the hexane isomers a large number of elementary reactions is incorporated. A simplified scheme of possible elementary reactions involved in the int erconversion of n-hexane and 3-methylpentane is given i n Figure 8. With the reaction conditions applied in this investigation it has been proved that the platinum will not display any isomerization and hydrogenolysis activity (Christoffel, 1978). Furthermore, product distributions frorn ethylbenzene isomerization (Robschlager and Christoffel, 1979) and 1,2-dimethylcyclopentane ring opening (Christoffel and Robschlager, 1978) on the same 0.35 w t % Pt-A1203 catalyst, which was used in this study, were in agreement with predictions from a reaction mechanism in which nonclassical carbocations (Olah, 1974) were involved. Thus in the reaction scheme
434
-
Ind. Eng. Chem. Prod. Res. Dev., Vol. 19, No. 3, 1980
/\/\/ +H2 l l - H 2
r
II
Ioca1 zai
,o
P
l +
Pt
H3C -CH2-CH2-CH2-CH2-CH3
11
' 1 r
C-H-de
i+r
\
.-CH
l+
+H2
Pt F i g u r e 8. Proposed reaction mechanism for n-hexane-3-methylpentane isomerization.
(Figure 8) two routes are considered by which carbenium ions are formed, one route via the classical bifunctional mechanism and one route via an electrophilic attack of a hydrogen ion at a C-H bond of an alkane and subsequent H2abstraction. The rearrangements of the carbenium ions are formulated with nonclassical carbocations as intermediates in accordance with the earlier results, whereby only C-H delocalizations are taken into account. The usual approach of the kinetic analysis of reaction systems like this is the concept of the rate-determining step. In the above reaction scheme the rearrangements of carbenium ions on acidic sites are commonly identified as rate-determining steps, which have the same stoichiometric numbers as the overall reactions. Chevalier et al. (1976) pointed out that the stability of carbenium ions decreases from tertiary toward primary ones and that protonated cyclopropanes are more stable than primary and less stable than secondary carbenium ions. Thus the rate constants connected with the interconversion of the dimethylbutanes should be larger than the rate constants for the interconversion of the methylpentanes. Skeletal isomerizations of n-hexane yielding 2- and 3-methylpentanes and of 2-methylpentane yielding 2,2- and 2,3dimethylbutanes may proceed via nonclassical carbocations without forming a primary carbenium ion, whereas the interconversion of 3-methylpentane and 2,2-dimethyl-
C-C-deloca
CH3
1 , zo t , o n
J
L
:
II
1
-3
H3
i l +
F i g u r e 9. Proposed reaction mechanism for 3-methylpentane-2,2dimethylbutane interconversion v i a nonclassical carbocations.
butane even via protonated cyclopropanes involves a primary carbenium ion intermediate (Figure 9). Skeletal isomerization of 3-methylpentane yielding 2,3-dimethylbutane requires more than only on alkyl shift. As the results of Table I reveal, all three sets of parameters reflect the predictions of a reaction mechanism in which rearrangements of carbocations are the rate-determining steps, although some of the corresponding rate constants have different values. The rate constants obtained from the Wei-Prater method are most reliable because they are only affected by experimental error. The rate constants of the slow interconversions of the methylpentanes into dimethylbutanes agree well with the proposed reaction mechanism. The relative rate constant of 3-methylpentane-2,3-dimethylbutane isomerization, which should be zero, has a very small value of 0.05. 3-Methylpentane-2,2 dimethylbutane isomerization which proceeds via the formation of a primary carbenium ion is connected with a very small rate constant, too, whereas the interconversions of 2-methylpentane and the dimethylbutanes, which proceed without formation of primary carbenium ions, are connected with larger rate constants. Literature Cited Bard, Y., "Nonlinear Estimation and Programming", 360 D 13.6.003., 1967. Bard, Y., "Nonlinear Parameter Estimation", Academic Press, New York, 1974. Chevalier, F.,Guisnet, M., Maurei, R., Proc. 6th Int. Congr. Catal. London, A22 (1976). Christoffel, E. G., Ind. Eng. Chem. Prod. Res. Dev., 18, 143 (1979). Christoffel, E. G., Chem. Z . , 102, 391 (1978). Christoffel, E. G.,Robschlager, K. H.. Ind. Eng. Chem. Prod. Res. Dev., 17, 331 (1978). Froment. G. F., AICHE J., 21, 1041 (1975). Gavalas, G. R., AICHE J., 19, 214 (1973). Himmeiblau, D. M., "Process Anatysis by Statistical Methods", Wiiey, New York, 1970. Oiah, G. A,, "Carbokationen und elektrophile Reaktionen", Verlag Chemie, Weinheim, 1974. Robschlager, K. H.,Christoffel. E. G., Ind. Eng. Chem. Prod. Res. Dev., 18, 347 (1979). Surjo, I., Christoffei, E. G., Chem. Ins. Tech., 50, 882 (1978). Wei, J. Prater, C. D., Adv. Catal., 13, 206 (1962).
Received for review July 2, 1979 Accepted April 7, 1980
