- Email: [email protected]
Direct integral method via splne-approximation for estimating rate constants
139
Applied Catalysis, 2 (1982)
139-164 EbevierScientificPublishiigCompany,AmsterdamF’rintedinTheNetherlands
ICNFoRBSTIMN'INGPATB~N'7TS DIRBCIIhTECPALWEXHCDVIASPLINB-APPKXIWXr A.
yEI(MAKovA*, S. vMDA**and P. vALKC**
U.S.S.R. *;fnstitute of Catalysis,630090Novosibirsk, *"B$tG
LorandUniwrsity, Laboratoryfor &emical Cybernetics, tieurnkrt. 6-8,
H-1088Budapest,Hungary. (Received 26 Way 1981,accepted21 September1981)
The problem of estimating kinetic parameters from integral reactor data is considered. The direct integral method is extended to rate models, where the parameters appear nonlinearly (e.g., Langmuir-Hinshelwood and Hougen-Watson kinetics). The main advantage is that no repeated numerical solution of the differential equations is required. Results show that in svite of its simplicity the direct integral method yields reasonable good estimates when applied to some kinetic models of practical importance.
In reaction engineering applications it is frequently necessary to estimate kinetic parameters from batch or integral data.
reactor
The problem can be posed as follows. The system is described
by the ordinary vectorial differential equation
(1) g(O)
= Lro,
where x is n-dimensional state vector and k is the p-dimensional parameter vector. Measurements on I are made m times at
lfi
tlJtZ,...,tm:
= c(t$
+ Q’
(2)
OlSS-9834/82/000+0000/$02.76 8 1982RbevierScientificPubli~1hingCompany
140
where
is
E.
The probkm
the vector of,experimental errors at
is to fit the model (1)
t = ti.
to the observations
Yi' i = 1,2,...,m in some optimal way. The most frequently used criterion is that of the least squares of deviations. Since in integral reactors we do not measure the values of the derivatives, Eqs. (I) cannot be used directly. This difficulty,may be overcome in one of the following ways Cl]: (i) fitting the solution of the rate equations (1) (indirect method 1;
(ii)
method); or
differentiation of data (direct differential
(iii) fitting the equivalent integral equations
to data (direct integral method). The classical and well-documented approach is the indirect integration of the rate equations for use in a method, i.e.,* nonlinear least squares algorithm (see. e.g., [I-31). of this method requires repeated numerUsually, application ical integration of rate equations, a procedure that frequently exhibits considerable numerical difficulties (see, e.g., [4j). If we
wish to use a gradient-type method for the minimization
of the objective function, we must compute also.the derivatives of r with respect to the parameters. Therefore, each Step in the iterative procedure requires complex and troublesome calculations. Because of the difficulties associated with the indirect approach, the direct differential method is also frequently used in applications. The direct differential method is based on the approximation of the derivatives
5 by numerical differentiation of the experimental data, a highly inaccurate procedure. Differentiation usually involves some smoothing procedure (see, e.e.,C51) and the parameter estimates are sensitive to the particular choice of the algorithm. The main advantage of the direct differential approach is that the estimation is performed using Eqs. (1) directly, hence the computation can be performed much faster. The idea of the direct integral approach has been-intro-
duced by Himmelblau et al., ~61
for cases in which rate expressions are linear in the parameters. The method has
been improved applying cubic spline interpolation ~'71. The purpose of this work is to qive
a straightforward'
extension of the method proposed by Hiasaelblau et al., to the
141
problem of estimating rate constants in general kinetic models (e.g., Langmuire-Hinshelwood-type
kinetics) from batch or integral .
reactor data.
THE DIRECT INTEGRAL METHOD The differential equations (1) may be transformed into the equivalent integral equations
The basic idea of the direct integral method is to approximate the integrands in Eq. (3) by functions interpolating the points f(yi,k),
i=1,2, ....m. Then the integrals on the right hand
side may be evaluated and relations (3) algebraic equations in the parameters Let
2%
can be regarded as
k.
-denote the
n-vector of natural cubic splines interpolating the values f(yi,k), i=7,2,... ,m, over the interval
to 2 t 2 tm CBI.
Define the vectors
g(k)=
t2
I
f
~kC~)d~
(4)
t0
tm f ! gk(T)dT tO
According to the least squares principle we shall minimize the objective function
J = (2
-
?~(k))*(x --
Q(k)) --
(51
where superscript T stands for transposed. As is well known Cll, one of the most efficient (in terms of total number of function evaluation) methods to minimize
142 quadratic objective functions is the Gauss-Marquardt procedure 191.
If we wish to use this method, we must compute also the
Spline-interpolation gives a great advanQ(k). -ag/a&. Changing the tage when computing the derivatives
derivatives of
order of differentiation and integration, and taking into account that cubic splines over fixed nodes form a linear space, we obtain the relation t.
= fz gkaf’akj 0
CT)
a~,
(6)
af/ak;
denotes the n-vector of natural cubic Skv spline functions interpolating the values af(Yi,&)/a&jj, i=l,Z,...,m over the interval to 2 t 2 t,,,. Let XC k) where
denote the Jacobian of
tl aflak, f -Sk tO
X(k) --
=
then by virtue of Eq. (6) _
Q(k), --
(T)dt
. . .
tl
/
aflak _sk
P(T)dT
tO
t2 aflak (T)dT l Sk +
t2 .
.
.
bO
/
,kaf’akP
(T1d-c
(7)
tO
%I f s afm -k
1 (T)dT.
.
.
tin /
afm
Sk
p(T)dT
tO
tO
Thus the Jacobian tion technique as
X(k) can be computed by the same interpola-used to evaluate the functions Q(k).
The proposed algorithm of the estimation procedure is as fol lows: 1.
Select a first quess f(Yi'k)
and
af(yi,k)/akj
k,
compute the values i=l,Z,...,m;
j=l,Z,.
. .,p.
143 2.
Determine the interpolating natural cubic splines 2f/ak
and
j J
. .,p.
j=1,2,.
82
3.
Compute the integrals and form according to Eq. (4)
4.
and
8(k) --
X(k) --
and Eq. (71, respectively.
Generate the new estimate
_k^
by the Gauss-Marquardt
method
;; - &
= C&k)X(k)+X~21-1 ---
&&
,Cu
where the selection of the parameters
-
Cp(_k )I,
AandB2
(8)
is
discussed in the literature Cll. 5.
Return to step 2 with until
Ii - 51
g
in place of
6 and continue
is less than some specified treshold.
Notice that the algorithm given by Himmelblau et al. C21 for linear estimation problems is a special case of the above procedure. Another specialization can be developed for cases in a& the right hard sidesof EC+ (I) are homogeneous functions of first order in the parameters. Then setting
X=0
Eq. (8)
is reduced to the fixed-point iteration process
;; = c&k)x(k)+ ---
&)Y.
-
(9)
Further details on this special case are given in ClOI.
EXAMPLES To test the effectiveness of the direct integral approach, scdnenumerical experiments have been performed using rate equations of practical importance. Example
I
The kinetic model of the selective hydrogenation of o-cresol suggested by Zwicky et al., c11,12lhas the following form:
144 dc -
=
dt’
kl cl
-VI
de2 dt=
k
C~~Q~C~+Q2c 3
mk
klCl - k2QlC2 cI+QI~2+Q2~3
“k
c~+Q~c~+Q~c~
de
(10)
k2QIC2
-zG=
where k
i
-
i=l,
k;U’)W,pL
K @(T,p)
= +
k;
=
‘i
= QiM
KH + T p)-1
p(1
(5
-
&)
($
-
+q)
($
-
R
23630
C
-
R
2.
We used the
(11)
1
A)]
exP
The 8 parameters to estimate are i=l,
1
EQi exp
KH = 0.0151 exp
H = 0.3466
2, ;
kTM,
QiM,
Eke,
and E
Qi’
"data"
generated by numerkal integration of the rate equations (10) with parameter values given by Zwicky
Cl21
and shown in
Table
1. For
each temperature, 5 runs of 10 observations were simulated with 1% relative error drawn from normal distribution. The analysis of "data" were carried out in two steps. First, the rate constants were estimated in each run by the direct integral method. The summary of results is given in Table
2.
secondly, activation energies and preexponential
145 TABLE
1
"True" values of consaxntsevaluated from relations (11) p = 40 bar,
mk = 5 %.
T, R Constant 366.2
393.2
423.2
mk * kl
0.0098
0.032
0.09
mk * k2
0.014
0.027
0.046
Ql
0.057
0.087
0.13
Q2
0.011
0.029
0.072
factors were evaluated by fitting a line to the Arrhenius-plot. Results are shown in
Table
3.
It can be seen that
the empirical statistical properties of the estimates are reasonably good.
Exampk?e
2
Consider the reaction system
Areaction
mechanism of this form has been proposed to
describe the selective hydrogenation of cotton seed oils ~13,147,
The mass-action type kinetic model is
% at
_ - -
(kl
+ $1
clcR2
2
calstmtsinmodel
S-
0.08701
0.02848
4.89
1.18
o.am33
0.01087
4
0.02680
0.00251
4
-1.00
0.05421
0.032cO
O.ooO23
O-31
relative mean value bias (%)
yk2,min -1 0.01414
deviaticm
O.COlO8
0.00309
O.CXXA6
o.azJO15
0.07259
0.12943
-0.01
1.79
0.04631
0.09004
mean value
0.74
0.0
standard relative deviaticn bias (%I
(5 mns)
-0.44
a.82
0.00182
-0.67
-0.05
0.00143
o.OcO42
OXCO25
stadard relative deviat1cxl bias (%I
T = 423.2 K
simulated data with 1% error by
T = 393.2 K (5 runs)
O.oooo4
mean value
(5 runs)
(lo), obtained fmn
0.00977
~*$,min-l
cu-mant
T = 366.2 K
the&reztintegralGauss~ardtmethod
Isothmalestimatesof
TABLE
147 TABLE
3
Estha+sdvalues
of&-rheniusparanetersfrunrateccxxtants&tained
bythedirEctintsgralmethcdandArrhEniusp10ts
constant
"true I' value
estimated value
relative bias (%)
klM
0.00994
0.01007
-1.31
Ekl
57530
58352
-1.43
* k2M
0.00838
0.00844
-0.72
Ek2
34650
34978
-0.95
*lM
0.0868
0.08701
-0.24
19660
-4.52
l
18810
EQ1
*2M
0.0289
0.02848
42320
EQ2
dC2 -x
= klCICH
dc3 dt
= k2clcH2
1.45
42958
2
- (kg
-1.51
+ k4cH2 )c 2+k
+ k3c2
-
(k_3
-3c3
+ k5cH
)c3
(12)
2 dc4 aF
Ci(U)
= k4c2cH
+ k5c3cH
2 = cio
;
2 i=l,2,3,4.
Taking into account mass transfer, the concentration of in the solution is described by
H2
in m&l
0.0281
0.3050
B
k5
0.0154
-1.65
0.0214
0.5014
k4
0.8912
-0.27
0.0543
0.3943
k3
0.98
1.43
-0.35
0.0218
1.2544
k2
0.63
relative bias (%)
0.0191
st. dev.
l%error(5nms)
of prmeters
0.9937
Value
msan
values
!l
anstant
Vn=thaJ
Estimated
TABLE 4
0.8908
0.3060
0.4960
0.0163
0.0452
0.0396
0.1091
0.0035
1.2544
0.3940
0.0302
Sit. dev.
2% ermr
obtained
1.0156
man v&he
(12),(14)
1.02
-2.3
0.39
1.51
-0.44
-1.56
relative bias (%)
(5 runs)
0.8895
0.3336
0.4680
0.3066
1.2310
1.0184
V&LIE
mean
0.0099
0.0587
0.0674
0.1859
0.0638
0.0429
1.17
-11.19
6.40
22.84
1.52
-1.84
relative bias (%)
(5 runs)
st. dev.
3% error
frun simulated data by the direct integral
0.90
0.30
0.50
0.40
1.25
1.00
"true" values
!z m
149 dc H2
-
=
,6(c
Y B2
dt
-
c~2)-C(kl+k2)Clik~C2i~
SC3
(13)
%2
Applying the quasi steady-state approximation to Eq. (131, we obtain * 'X2 =
8
‘HZ
(14)
k4c2 + k5c3
B + (kl+k2)5
The equilibrium constant
k_3/k3 = 1.5 is supposed to be known from the literature. The problem is to find 6 unknown
constants
and
kIJk2,k3,k4,k5
from the concentrations C
2.0,
2.5,
3.0,
4.0,
with initial' conditions
i3 with known i=1,2,3,4, t=O.
25,
5.0,
CI(0) = 1.0,
= 1.0
observid. "Data"
0.5,
6.0,
~1 0.75,
7.0,
8.0,
Ci(0) = 0,
1.0, 10.0 i = 2,3,4.
To "exact" values l%, 2% and 3% relative error was added simulating 5 runs for each case. The summary of results is given in
TabLe 4.
The estimates at error level less or
equal 2% are only slightly biased with a moderate standard deviation. At 3% error some of the estimates become strongly biased. It seems reasonable to require experimental error less than 2% to obtain reliable estimates of parameters in kinetic models of similar ccunplexity.
Example
3
Here the sensitivity of the method to the initial guess is studied. Data obtained in Four different initial
EotampLe
guesses
2
with 2% error were used.
were generated using uniformly
distributed random numbers in the interval O-10 (see TabLe
5.)
In two cases the iteration converged to the global minimum, in one case a local minimum was found, and in one case the program was unable to find a finite minimum (estimates runned out of reasonable range). It is seen that the direct integral method is faced with the usual convergence problem of nonlinear para6 an interesting example meter estimation. However, in Tab& is shown where convergence was obtained in spite of negative in the course of the search. This k5
values of the parameter
could be hardly the case if numerical solution of differential equations were involved.
160
151 TABLE 6 A typical iteration process (model (12),(14), 2% error)
cbjective
kl
k4
k5
6
function
0
9,2703
0,1749
4,1427
8,158O
4,0036
9,8939
1
7,7265
0,3318
3,7286
6,708O
3,0232
6,8528
771,72
2
7,0097
0,3988
3,5928
5,922O
2,2514
4,6367
438,06
4
4,7771
2,0927
3,237O
2,6293
-0,5876
0,4575
1,5807
6
2,4587
3,4384
3,2914
1,7724
-0,6177
0,5298
0,154O
8
1,5862
2,9448
3,0875
1,5882
-0,5334
0,5822
0,09626
12
0,8972
1,3999
1,3327
0,8488
-0,0751:
0,8041
0,04682
14
0,9534
1,2432
0,5989
0,5744
0,2075
0,9032
0,005383
18
0,9763
1,2135
0,4l38
0,5077
0,2827
0,9241
0,003212
20
0,9772
1,2123
0,4064
0,5051
0,2857
0,9248
0,003203
25
0,9775
1,2121
0,4044
0,5044
0,2866
0,9251
o,ax!O2
30
0,9775
1,2121
3,4044
0,5044
0,2866
0,9251
o,ax!O2
DISCUSSION In spite of its simplicity, has been
somewhat
are as follows: observed,
overlooked
the concentration
the method
statistically
the direct
biased
requires estimates
integral
in the literature. of every
dense
data,
method
Usual
component
claims must
it gives. rise
and unassessable
errors.
to
be
152
The principal condition is the observability of all the concentrations. Notice, however, that typical nonlinear kinetic models (e.g., Langmuir-Hinshelwood, Hougen-Watson, Temkin kinetics)
are obtained by eliminating all the unmeasurable
variables. To obtain satisfactory approximations of integrals, the observations must cover the region of interest. It was, however, shown that integration through cubic splines results in rather good estimates even from a limited number of data points
C71.
From a statistical point of view, the basic shortcoming of the direct integral approach is that the estimates are biased due to evaluating the integrals from data instead of "true" values. Notice, however, that the structure of the rate equations (3)
is usually much simpler than that of the integrated
rate equations. Therefore, in a number of cases we may obtain estimates of the parameters with considerably smaller mean square error than the estimates given by the indirect method. AS usual in nonlinear estimation theory, no exact statistical statements can be made on the estimates. In addition, the right hand side of Eq. (3) is approximated from data and hence we obtain only estimated values of the residual error. Exact residuals may be evaluated by solving Eqs.
(1)
with the
estimated parameters _. ^k Let 8r denote the variance of the residuals, then the covariance matrix of the estimates can be roughly approximated by
In principle, the direct integral method can be applied to any
kind of nonlinearity. It is, however, necessary to
note that reliable parameter estimates have been obtained only for moderate nonlinearities. Estimating activation energies from non-isothermal reactor data or reaction orders is beyond the comprehension of the method. It should be noted that the estimates of rate constants frcanmultiresponse data usually may be improved by minimizing the determinant of the matrix of sums of squares and cross products Cl51. The determinant criterion is superior in most cases to the sum of squares criterion, because it allows the
163 responses to have unequal variances, and takes into account correlation among the responses. The extension of the direct integral method for this objective function is straightforward. In the examples of the present paper one could, however, expect no significant advantage from using the determinant criterion as measurement errors are uncorrelated ~161. Finally, we note that the relative performance of the different methods obviously deserves to be more deeply investigated. This problem is,however, far from being simple because the performance of a particular method is much depending upon the structure of the kinetic model, the design of experiments and the statistics of measurement errors. Therefore, comparison of the methods is not considered in full in the present paper and requires further experience with the proposed integral method.
CONCLUSION The resultsof this work show that the general direct integral method is a useful and easy device for parameter estimation. The method has been effectively applied to a number of simulated kinetic problems. Since the procedure do not require numerical solution of differential equations, the computation can be performed much faster than in case of the indirect method. The direct integral method is particularly useful when the computation capacity is limited. A perspective field of application is discriminating several kinetic models, where the great advantage of computational simplicity may outweigh the disadvantage of somewhat moderate accuracy. NOTATIONS ( in Example 1: ) T I1'M
tenperature, K nean tenperature, K
(TM= 393.2)
t
tine,
min
P c. 2
mole fraction
ki
rate constant,(min * % cat)-1
Qi
quotientofequilibriumconstants
pressure,bar
164 Henry's ccnstant,m3 (bar * ml)
H
m
catalystccncentraticm, w %
n
-1 gas ccnstant,J(sol - K)
R
(
(R = 8.314)
in Example 2: )
t
tine,min
c. 2
ccncentration, kn-ol- I&-~
ki
ki
* cH2 8
-1
(i=1,2,4,5)
rate constants,d
(i=-3,3)
-1 rate ccnstants,min
(km31 . min)-l
equilibriirn H2 concentratimin solution,knol .
m-3
-1 transfercoefficient, min
REFERENCES 1
Y. Bard, Nonlinear Parameter Estimation.
2
Academic Press N.Y. 1974. W.E. Ball, L.C.D. Groenweghe Ind. Eng. Chem. Fundamentals, 5 (1966) 181.
3 4
I.H. Seinfeld, L. Edsberg,
Ind. Eng. Chem. (1970) 32.
Numerical Methods for Mass Action Kinetics.
Int*Numerical Methods for Differential Systems Ed. by L. Lapidus and W.E. Schiesser
Academic Press, N.Y. 1976. AIChe Journal 12 (1966), 1110.
5
N. De Nevers,
6
D.M. Himmelblau,
C.R. Jones, K.B. Bischoff,
Ind. Eng. Chem. Fundamentals, (1967) 539. 7 8
Y.P. Tang,
Ind. Eng. Chem. Fundamentals
10 (1971) 321. Spline Functions, Interpolation and Numerical Quadrature. Mathematical Methods for Digital T.N.E. Greville,
Computers, Volume 2. John Wiley & Sons, N.Y. 1967. 9
D.W. Marguardt,
10
P. Valk6, A. Yermakova, S. Vajda,
11
J.J. zwickyi
12
14
J.J. Zwicky, Th. Biihlmann, M. Egli, G. Gut, Chimia Kenji Hashimoto, Masaaki Teremoto, Shinhi Nagota, J. Chem. Eng. Japan 4 (1971) 150. React. Kinet. Catd.Lett. A. Yermakova, A.S. Umbetov,
15 16
14 (1980) 187. G.E.P. Box, N.R. Draper, Biometrika 52 (1965) 355. J. Erjavec, Ind. Eng. Chem. Fundamentals 9 (1970) 187.
13
J. SIAM 11 (1963) 431.
Hungarian Journal of Ind. Chem. (1980) 437. G. Gut, Chem. Eng. Sci.(33) 1978
1363. 31 (1977) 22
Applied Catalysis, 2 (1982)
139-164 EbevierScientificPublishiigCompany,AmsterdamF’rintedinTheNetherlands
ICNFoRBSTIMN'INGPATB~N'7TS DIRBCIIhTECPALWEXHCDVIASPLINB-APPKXIWXr A.
yEI(MAKovA*, S. vMDA**and P. vALKC**
U.S.S.R. *;fnstitute of Catalysis,630090Novosibirsk, *"B$tG
LorandUniwrsity, Laboratoryfor &emical Cybernetics, tieurnkrt. 6-8,
H-1088Budapest,Hungary. (Received 26 Way 1981,accepted21 September1981)
The problem of estimating kinetic parameters from integral reactor data is considered. The direct integral method is extended to rate models, where the parameters appear nonlinearly (e.g., Langmuir-Hinshelwood and Hougen-Watson kinetics). The main advantage is that no repeated numerical solution of the differential equations is required. Results show that in svite of its simplicity the direct integral method yields reasonable good estimates when applied to some kinetic models of practical importance.
In reaction engineering applications it is frequently necessary to estimate kinetic parameters from batch or integral data.
reactor
The problem can be posed as follows. The system is described
by the ordinary vectorial differential equation
(1) g(O)
= Lro,
where x is n-dimensional state vector and k is the p-dimensional parameter vector. Measurements on I are made m times at
lfi
tlJtZ,...,tm:
= c(t$
+ Q’
(2)
OlSS-9834/82/000+0000/$02.76 8 1982RbevierScientificPubli~1hingCompany
140
where
is
E.
The probkm
the vector of,experimental errors at
is to fit the model (1)
t = ti.
to the observations
Yi' i = 1,2,...,m in some optimal way. The most frequently used criterion is that of the least squares of deviations. Since in integral reactors we do not measure the values of the derivatives, Eqs. (I) cannot be used directly. This difficulty,may be overcome in one of the following ways Cl]: (i) fitting the solution of the rate equations (1) (indirect method 1;
(ii)
method); or
differentiation of data (direct differential
(iii) fitting the equivalent integral equations
to data (direct integral method). The classical and well-documented approach is the indirect integration of the rate equations for use in a method, i.e.,* nonlinear least squares algorithm (see. e.g., [I-31). of this method requires repeated numerUsually, application ical integration of rate equations, a procedure that frequently exhibits considerable numerical difficulties (see, e.g., [4j). If we
wish to use a gradient-type method for the minimization
of the objective function, we must compute also.the derivatives of r with respect to the parameters. Therefore, each Step in the iterative procedure requires complex and troublesome calculations. Because of the difficulties associated with the indirect approach, the direct differential method is also frequently used in applications. The direct differential method is based on the approximation of the derivatives
5 by numerical differentiation of the experimental data, a highly inaccurate procedure. Differentiation usually involves some smoothing procedure (see, e.e.,C51) and the parameter estimates are sensitive to the particular choice of the algorithm. The main advantage of the direct differential approach is that the estimation is performed using Eqs. (1) directly, hence the computation can be performed much faster. The idea of the direct integral approach has been-intro-
duced by Himmelblau et al., ~61
for cases in which rate expressions are linear in the parameters. The method has
been improved applying cubic spline interpolation ~'71. The purpose of this work is to qive
a straightforward'
extension of the method proposed by Hiasaelblau et al., to the
141
problem of estimating rate constants in general kinetic models (e.g., Langmuire-Hinshelwood-type
kinetics) from batch or integral .
reactor data.
THE DIRECT INTEGRAL METHOD The differential equations (1) may be transformed into the equivalent integral equations
The basic idea of the direct integral method is to approximate the integrands in Eq. (3) by functions interpolating the points f(yi,k),
i=1,2, ....m. Then the integrals on the right hand
side may be evaluated and relations (3) algebraic equations in the parameters Let
2%
can be regarded as
k.
-denote the
n-vector of natural cubic splines interpolating the values f(yi,k), i=7,2,... ,m, over the interval
to 2 t 2 tm CBI.
Define the vectors
g(k)=
t2
I
f
~kC~)d~
(4)
t0
tm f ! gk(T)dT tO
According to the least squares principle we shall minimize the objective function
J = (2
-
?~(k))*(x --
Q(k)) --
(51
where superscript T stands for transposed. As is well known Cll, one of the most efficient (in terms of total number of function evaluation) methods to minimize
142 quadratic objective functions is the Gauss-Marquardt procedure 191.
If we wish to use this method, we must compute also the
Spline-interpolation gives a great advanQ(k). -ag/a&. Changing the tage when computing the derivatives
derivatives of
order of differentiation and integration, and taking into account that cubic splines over fixed nodes form a linear space, we obtain the relation t.
= fz gkaf’akj 0
CT)
a~,
(6)
af/ak;
denotes the n-vector of natural cubic Skv spline functions interpolating the values af(Yi,&)/a&jj, i=l,Z,...,m over the interval to 2 t 2 t,,,. Let XC k) where
denote the Jacobian of
tl aflak, f -Sk tO
X(k) --
=
then by virtue of Eq. (6) _
Q(k), --
(T)dt
. . .
tl
/
aflak _sk
P(T)dT
tO
t2 aflak (T)dT l Sk +
t2 .
.
.
bO
/
,kaf’akP
(T1d-c
(7)
tO
%I f s afm -k
1 (T)dT.
.
.
tin /
afm
Sk
p(T)dT
tO
tO
Thus the Jacobian tion technique as
X(k) can be computed by the same interpola-used to evaluate the functions Q(k).
The proposed algorithm of the estimation procedure is as fol lows: 1.
Select a first quess f(Yi'k)
and
af(yi,k)/akj
k,
compute the values i=l,Z,...,m;
j=l,Z,.
. .,p.
143 2.
Determine the interpolating natural cubic splines 2f/ak
and
j J
. .,p.
j=1,2,.
82
3.
Compute the integrals and form according to Eq. (4)
4.
and
8(k) --
X(k) --
and Eq. (71, respectively.
Generate the new estimate
_k^
by the Gauss-Marquardt
method
;; - &
= C&k)X(k)+X~21-1 ---
&&
,Cu
where the selection of the parameters
-
Cp(_k )I,
AandB2
(8)
is
discussed in the literature Cll. 5.
Return to step 2 with until
Ii - 51
g
in place of
6 and continue
is less than some specified treshold.
Notice that the algorithm given by Himmelblau et al. C21 for linear estimation problems is a special case of the above procedure. Another specialization can be developed for cases in a& the right hard sidesof EC+ (I) are homogeneous functions of first order in the parameters. Then setting
X=0
Eq. (8)
is reduced to the fixed-point iteration process
;; = c&k)x(k)+ ---
&)Y.
-
(9)
Further details on this special case are given in ClOI.
EXAMPLES To test the effectiveness of the direct integral approach, scdnenumerical experiments have been performed using rate equations of practical importance. Example
I
The kinetic model of the selective hydrogenation of o-cresol suggested by Zwicky et al., c11,12lhas the following form:
144 dc -
=
dt’
kl cl
-VI
de2 dt=
k
C~~Q~C~+Q2c 3
mk
klCl - k2QlC2 cI+QI~2+Q2~3
“k
c~+Q~c~+Q~c~
de
(10)
k2QIC2
-zG=
where k
i
-
i=l,
k;U’)W,pL
K @(T,p)
= +
k;
=
‘i
= QiM
KH + T p)-1
p(1
(5
-
&)
($
-
+q)
($
-
R
23630
C
-
R
2.
We used the
(11)
1
A)]
exP
The 8 parameters to estimate are i=l,
1
EQi exp
KH = 0.0151 exp
H = 0.3466
2, ;
kTM,
QiM,
Eke,
and E
Qi’
"data"
generated by numerkal integration of the rate equations (10) with parameter values given by Zwicky
Cl21
and shown in
Table
1. For
each temperature, 5 runs of 10 observations were simulated with 1% relative error drawn from normal distribution. The analysis of "data" were carried out in two steps. First, the rate constants were estimated in each run by the direct integral method. The summary of results is given in Table
2.
secondly, activation energies and preexponential
145 TABLE
1
"True" values of consaxntsevaluated from relations (11) p = 40 bar,
mk = 5 %.
T, R Constant 366.2
393.2
423.2
mk * kl
0.0098
0.032
0.09
mk * k2
0.014
0.027
0.046
Ql
0.057
0.087
0.13
Q2
0.011
0.029
0.072
factors were evaluated by fitting a line to the Arrhenius-plot. Results are shown in
Table
3.
It can be seen that
the empirical statistical properties of the estimates are reasonably good.
Exampk?e
2
Consider the reaction system
Areaction
mechanism of this form has been proposed to
describe the selective hydrogenation of cotton seed oils ~13,147,
The mass-action type kinetic model is
% at
_ - -
(kl
+ $1
clcR2
2
calstmtsinmodel
S-
0.08701
0.02848
4.89
1.18
o.am33
0.01087
4
0.02680
0.00251
4
-1.00
0.05421
0.032cO
O.ooO23
O-31
relative mean value bias (%)
yk2,min -1 0.01414
deviaticm
O.COlO8
0.00309
O.CXXA6
o.azJO15
0.07259
0.12943
-0.01
1.79
0.04631
0.09004
mean value
0.74
0.0
standard relative deviaticn bias (%I
(5 mns)
-0.44
a.82
0.00182
-0.67
-0.05
0.00143
o.OcO42
OXCO25
stadard relative deviat1cxl bias (%I
T = 423.2 K
simulated data with 1% error by
T = 393.2 K (5 runs)
O.oooo4
mean value
(5 runs)
(lo), obtained fmn
0.00977
~*$,min-l
cu-mant
T = 366.2 K
the&reztintegralGauss~ardtmethod
Isothmalestimatesof
TABLE
147 TABLE
3
Estha+sdvalues
of&-rheniusparanetersfrunrateccxxtants&tained
bythedirEctintsgralmethcdandArrhEniusp10ts
constant
"true I' value
estimated value
relative bias (%)
klM
0.00994
0.01007
-1.31
Ekl
57530
58352
-1.43
* k2M
0.00838
0.00844
-0.72
Ek2
34650
34978
-0.95
*lM
0.0868
0.08701
-0.24
19660
-4.52
l
18810
EQ1
*2M
0.0289
0.02848
42320
EQ2
dC2 -x
= klCICH
dc3 dt
= k2clcH2
1.45
42958
2
- (kg
-1.51
+ k4cH2 )c 2+k
+ k3c2
-
(k_3
-3c3
+ k5cH
)c3
(12)
2 dc4 aF
Ci(U)
= k4c2cH
+ k5c3cH
2 = cio
;
2 i=l,2,3,4.
Taking into account mass transfer, the concentration of in the solution is described by
H2
in m&l
0.0281
0.3050
B
k5
0.0154
-1.65
0.0214
0.5014
k4
0.8912
-0.27
0.0543
0.3943
k3
0.98
1.43
-0.35
0.0218
1.2544
k2
0.63
relative bias (%)
0.0191
st. dev.
l%error(5nms)
of prmeters
0.9937
Value
msan
values
!l
anstant
Vn=thaJ
Estimated
TABLE 4
0.8908
0.3060
0.4960
0.0163
0.0452
0.0396
0.1091
0.0035
1.2544
0.3940
0.0302
Sit. dev.
2% ermr
obtained
1.0156
man v&he
(12),(14)
1.02
-2.3
0.39
1.51
-0.44
-1.56
relative bias (%)
(5 runs)
0.8895
0.3336
0.4680
0.3066
1.2310
1.0184
V&LIE
mean
0.0099
0.0587
0.0674
0.1859
0.0638
0.0429
1.17
-11.19
6.40
22.84
1.52
-1.84
relative bias (%)
(5 runs)
st. dev.
3% error
frun simulated data by the direct integral
0.90
0.30
0.50
0.40
1.25
1.00
"true" values
!z m
149 dc H2
-
=
,6(c
Y B2
dt
-
c~2)-C(kl+k2)Clik~C2i~
SC3
(13)
%2
Applying the quasi steady-state approximation to Eq. (131, we obtain * 'X2 =
8
‘HZ
(14)
k4c2 + k5c3
B + (kl+k2)5
The equilibrium constant
k_3/k3 = 1.5 is supposed to be known from the literature. The problem is to find 6 unknown
constants
and
kIJk2,k3,k4,k5
from the concentrations C
2.0,
2.5,
3.0,
4.0,
with initial' conditions
i3 with known i=1,2,3,4, t=O.
25,
5.0,
CI(0) = 1.0,
= 1.0
observid. "Data"
0.5,
6.0,
~1 0.75,
7.0,
8.0,
Ci(0) = 0,
1.0, 10.0 i = 2,3,4.
To "exact" values l%, 2% and 3% relative error was added simulating 5 runs for each case. The summary of results is given in
TabLe 4.
The estimates at error level less or
equal 2% are only slightly biased with a moderate standard deviation. At 3% error some of the estimates become strongly biased. It seems reasonable to require experimental error less than 2% to obtain reliable estimates of parameters in kinetic models of similar ccunplexity.
Example
3
Here the sensitivity of the method to the initial guess is studied. Data obtained in Four different initial
EotampLe
guesses
2
with 2% error were used.
were generated using uniformly
distributed random numbers in the interval O-10 (see TabLe
5.)
In two cases the iteration converged to the global minimum, in one case a local minimum was found, and in one case the program was unable to find a finite minimum (estimates runned out of reasonable range). It is seen that the direct integral method is faced with the usual convergence problem of nonlinear para6 an interesting example meter estimation. However, in Tab& is shown where convergence was obtained in spite of negative in the course of the search. This k5
values of the parameter
could be hardly the case if numerical solution of differential equations were involved.
160
151 TABLE 6 A typical iteration process (model (12),(14), 2% error)
cbjective
kl
k4
k5
6
function
0
9,2703
0,1749
4,1427
8,158O
4,0036
9,8939
1
7,7265
0,3318
3,7286
6,708O
3,0232
6,8528
771,72
2
7,0097
0,3988
3,5928
5,922O
2,2514
4,6367
438,06
4
4,7771
2,0927
3,237O
2,6293
-0,5876
0,4575
1,5807
6
2,4587
3,4384
3,2914
1,7724
-0,6177
0,5298
0,154O
8
1,5862
2,9448
3,0875
1,5882
-0,5334
0,5822
0,09626
12
0,8972
1,3999
1,3327
0,8488
-0,0751:
0,8041
0,04682
14
0,9534
1,2432
0,5989
0,5744
0,2075
0,9032
0,005383
18
0,9763
1,2135
0,4l38
0,5077
0,2827
0,9241
0,003212
20
0,9772
1,2123
0,4064
0,5051
0,2857
0,9248
0,003203
25
0,9775
1,2121
0,4044
0,5044
0,2866
0,9251
o,ax!O2
30
0,9775
1,2121
3,4044
0,5044
0,2866
0,9251
o,ax!O2
DISCUSSION In spite of its simplicity, has been
somewhat
are as follows: observed,
overlooked
the concentration
the method
statistically
the direct
biased
requires estimates
integral
in the literature. of every
dense
data,
method
Usual
component
claims must
it gives. rise
and unassessable
errors.
to
be
152
The principal condition is the observability of all the concentrations. Notice, however, that typical nonlinear kinetic models (e.g., Langmuir-Hinshelwood, Hougen-Watson, Temkin kinetics)
are obtained by eliminating all the unmeasurable
variables. To obtain satisfactory approximations of integrals, the observations must cover the region of interest. It was, however, shown that integration through cubic splines results in rather good estimates even from a limited number of data points
C71.
From a statistical point of view, the basic shortcoming of the direct integral approach is that the estimates are biased due to evaluating the integrals from data instead of "true" values. Notice, however, that the structure of the rate equations (3)
is usually much simpler than that of the integrated
rate equations. Therefore, in a number of cases we may obtain estimates of the parameters with considerably smaller mean square error than the estimates given by the indirect method. AS usual in nonlinear estimation theory, no exact statistical statements can be made on the estimates. In addition, the right hand side of Eq. (3) is approximated from data and hence we obtain only estimated values of the residual error. Exact residuals may be evaluated by solving Eqs.
(1)
with the
estimated parameters _. ^k Let 8r denote the variance of the residuals, then the covariance matrix of the estimates can be roughly approximated by
In principle, the direct integral method can be applied to any
kind of nonlinearity. It is, however, necessary to
note that reliable parameter estimates have been obtained only for moderate nonlinearities. Estimating activation energies from non-isothermal reactor data or reaction orders is beyond the comprehension of the method. It should be noted that the estimates of rate constants frcanmultiresponse data usually may be improved by minimizing the determinant of the matrix of sums of squares and cross products Cl51. The determinant criterion is superior in most cases to the sum of squares criterion, because it allows the
163 responses to have unequal variances, and takes into account correlation among the responses. The extension of the direct integral method for this objective function is straightforward. In the examples of the present paper one could, however, expect no significant advantage from using the determinant criterion as measurement errors are uncorrelated ~161. Finally, we note that the relative performance of the different methods obviously deserves to be more deeply investigated. This problem is,however, far from being simple because the performance of a particular method is much depending upon the structure of the kinetic model, the design of experiments and the statistics of measurement errors. Therefore, comparison of the methods is not considered in full in the present paper and requires further experience with the proposed integral method.
CONCLUSION The resultsof this work show that the general direct integral method is a useful and easy device for parameter estimation. The method has been effectively applied to a number of simulated kinetic problems. Since the procedure do not require numerical solution of differential equations, the computation can be performed much faster than in case of the indirect method. The direct integral method is particularly useful when the computation capacity is limited. A perspective field of application is discriminating several kinetic models, where the great advantage of computational simplicity may outweigh the disadvantage of somewhat moderate accuracy. NOTATIONS ( in Example 1: ) T I1'M
tenperature, K nean tenperature, K
(TM= 393.2)
t
tine,
min
P c. 2
mole fraction
ki
rate constant,(min * % cat)-1
Qi
quotientofequilibriumconstants
pressure,bar
164 Henry's ccnstant,m3 (bar * ml)
H
m
catalystccncentraticm, w %
n
-1 gas ccnstant,J(sol - K)
R
(
(R = 8.314)
in Example 2: )
t
tine,min
c. 2
ccncentration, kn-ol- I&-~
ki
ki
* cH2 8
-1
(i=1,2,4,5)
rate constants,d
(i=-3,3)
-1 rate ccnstants,min
(km31 . min)-l
equilibriirn H2 concentratimin solution,knol .
m-3
-1 transfercoefficient, min
REFERENCES 1
Y. Bard, Nonlinear Parameter Estimation.
2
Academic Press N.Y. 1974. W.E. Ball, L.C.D. Groenweghe Ind. Eng. Chem. Fundamentals, 5 (1966) 181.
3 4
I.H. Seinfeld, L. Edsberg,
Ind. Eng. Chem. (1970) 32.
Numerical Methods for Mass Action Kinetics.
Int*Numerical Methods for Differential Systems Ed. by L. Lapidus and W.E. Schiesser
Academic Press, N.Y. 1976. AIChe Journal 12 (1966), 1110.
5
N. De Nevers,
6
D.M. Himmelblau,
C.R. Jones, K.B. Bischoff,
Ind. Eng. Chem. Fundamentals, (1967) 539. 7 8
Y.P. Tang,
Ind. Eng. Chem. Fundamentals
10 (1971) 321. Spline Functions, Interpolation and Numerical Quadrature. Mathematical Methods for Digital T.N.E. Greville,
Computers, Volume 2. John Wiley & Sons, N.Y. 1967. 9
D.W. Marguardt,
10
P. Valk6, A. Yermakova, S. Vajda,
11
J.J. zwickyi
12
14
J.J. Zwicky, Th. Biihlmann, M. Egli, G. Gut, Chimia Kenji Hashimoto, Masaaki Teremoto, Shinhi Nagota, J. Chem. Eng. Japan 4 (1971) 150. React. Kinet. Catd.Lett. A. Yermakova, A.S. Umbetov,
15 16
14 (1980) 187. G.E.P. Box, N.R. Draper, Biometrika 52 (1965) 355. J. Erjavec, Ind. Eng. Chem. Fundamentals 9 (1970) 187.
13
J. SIAM 11 (1963) 431.
Hungarian Journal of Ind. Chem. (1980) 437. G. Gut, Chem. Eng. Sci.(33) 1978
1363. 31 (1977) 22