• No results found

Estimation of kinetical parameters in chemical initial value problems

N/A
N/A
Protected

Academic year: 2021

Share "Estimation of kinetical parameters in chemical initial value problems"

Copied!
20
0
0

Bezig met laden.... (Bekijk nu de volledige tekst)

Hele tekst

(1)

Estimation of kinetical parameters in chemical initial value

problems

Citation for published version (APA):

Molenaar, J. (1990). Estimation of kinetical parameters in chemical initial value problems. (IWDE report; Vol. 9008). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1990

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

(2)

"

t(iJ

Technische . ·~~ Unlversiteit Eindhoven Den Dolech 2

lnstituut

Wiskundige Dienstverlening

Eindhoven

REPORT IWDE 90-08

ESTIMATION OF KINETICAL PARAMETERS IN CHEMICAL INITIAL VALUE PROBLEMS

J.

Molenaar

(3)

IWDE report 90-08

Estimation of Kinetical Parameters

in Chemical Initial Value Problems

J. Molenaar

november 1990

The author wishes to acknowledge useful discussions with R. v.d. Rout, AKZO, Arnhem, and S.J.L. van Eijndhoven, TUE.

(4)

§1.

Introduction

The estimation ofkinetical parameters in chemical reactions is often a hard task. In [1] a first try is reported to tackle this problem. In this report we point out which problems commonly arise in this field and investigate possible improvements of frequently used procedures. The parameters are generally estimated by fitting experimental data with a theoretical model. In the experiment some reactants are brought together under well-defined conditions, i.e. the initial concentrations are known. After that, the concentrations of the reactants are measured at certain time points. In general it is not possible to measure the concentrations of all the reactants involved. In the mathematical model the kinetics of the reactions under consideration are described by a set of ordinary differential equations (ODE) with the kinetical constants as unknown parameters. For a given set of parameters the set of ODE's can be solved numerically. The differences between the measured and the calculated concentrations yield an indication of the reliability of the used parameters. The real parameter values are to be found from minimization of these differences. In this general scheme the following complications are often met:

- The determination of a starting point for the optimization procedure may be difficult. The complexer the model is, the closer the starting values should be to the real values in order to avoid the procedure to converge to an irrelevant local minimum.

- The adequacy of the model is not always clear in advance, so that sometimes many trials with various models are needed.

- The lack of data concerning some concentrations may introduce ambiguities, which pro-hibit the optimization procedure to converge. The less experimental data is available, the less succesful an estimation procedure will of course be.

- The measurements may be inaccurate.

In §2 we state the problem in quite a general form. Solutions have to be obtained numerically and an essential requirement from practice is to keep the calculation times reasonable. This implies that the product of the number of iteration steps in the optimization procedure and the time per step should be not too big. The reduction of computer time is an important issue in the following. The differences between measured and calculated concentration curves can be measured in several ways. In §3 we introduce and compare some alternative met-des. A general discussion on optimization techniques is given in §4. In §5 we present an extremely efficient method, which is applicable if all involved concentrations are measured. In this case the problem can be formulated as a set of algebraic equations. In §6 we cope with the common case of missing data by modifying the method of §5. In this modification numerical integration is involved, but only with respect to the evolution equations of the missing concentrations.

The system is often stiff, i.e. the reactions proceed only vehemently during a short time inter-val, the boundary layer, after which the system smoothly approaches equilibrium. In §7 we deal with a method in which the time consuming numerical integration through the boundary layer is avoided. In this section also the choice of the used norm is discussed. Furthermore, we present and discuss an alternative method to obtain the derivatives of the object function with respect to the parameters. The method is based on the idea of linearization. Section 8 contains numerical applications of the ideas in §§5-7 to a specific system, described in the Appendix. The conclusions are given in §9.

(5)

§2. The estimation problem

Here, we present the estimation problem of kinetical constants and introduce notations. The concentrations involved are denoted by

y(t)

=

(y

1

(t), ... , YN(t)l

and the kinetical parameters

by

.X

= ( A1 , ... ,AM

Chemical reactions are typically governed by an autonomous initial

value problem (IVP):

(2.1)

{

y

=

A(y).X

y(O)

=YO·

with A an N X M matrix. In the general case the r.h.s. of the ODE in (2.1) is given by a vector field, which is linear neither in

y

nor in

.X.

However, because a great deal of the practical problems is linear in

.X,

we prefer to use the present notation throughout this report. This restriction is not essential for the methods to be discussed. The solutions of (2.1.) are denoted by

y(

t, .X).

In order to describe a realistic situation, the vector field

A(y ).X

should be such, that the solutions of (2.1) exist for all t

>

0. In general A contains products of components of y. An example is given in the Appendix.

In practice some or all components of y are measured, let us say

Yi,

i

=

1, ... ,](,with]( :$

N,

at discrete time points

r;,j

= 1, .... ,L. We denote the measured values by

yf(r;).

The problem now is to determine

.X

such that the corresponding solution of the IVP defined above fits the data as well as possible. This might be formulated as:

minimize the function K L

c2.2)

F1(.X)

=

:L :L

Wij

(Rdr;

i=l i=l

with the residuals R1 given by

(2.3)

(Rl)ii

=

yf(r;)- Yi(r;,.X),

in which

Yi(t, .X)

is the i-th component of the solution of IVP (2.1) for given value of

.X.

In (2.2) we introduce a weighted least squares norm. The weights Wij ~ 0 can be used to correct for variations in the measurement accuracies of the different components and, if required, to take into account the relative, rather than the absolute errors. Throughout this report we assume all measurements to be equally accurate. For the weights we take

(2.4)

w· · -'3 - max .

ym(r·)

i 3 3

The minimization of

F1(.X)

is in practice notoriously hard. The evaluation of

F1

requires the numerical integration of the IVP, which is very time consuming. In the following we try to reduce the integration effort as much as possible.

(6)

§3. Norms

In (2.2) we introduced a certain metric in the space of concentration curves. Here, we discuss more in detail possible choices for the object function. Let us define , analogous to (2.3), the residual vector R by

(3.1)

R(t) = ym(t)- y(t,

A)

In practice this vector is used only at discrete times, but for notational ease we prefer here to deal with the continuous time case with 0

s

t

s

T. So, we assume the measurements to be continuously known at this interval. The transition from discrete to continuous time and vice versa is obvious for the formulae in this section. For the moment being we take the weights introduced in (2.2) equal to unity, thus WiJ = 1. This restriction is not at all essential. For fixed

t

we introduce the Holdernorm:

(3.2)

( )

1/p

IIR(t)llp

=

~I Ri(t)

lp ' p ~ 1

The norms IIR(t)ll are real, continuous functions oft on the interval 0

s

t

sT.

They belong to the space

Lp([O,

T]) with norm

(3.3)

IIRII,

= [ /

(IIR(t)ll,)"

dt]

11

The object function in (2.2) is nothing else but a discrete time version of this norm. It is by no means clear that the Lp norm is the most suitable one to measure the differences between the calculated solutions of IVP (2.1) and the measured concentrations. The Lp norm does not exploit all information contained in the data. To be specific, no use is made of the fact that the data are ordered in time, i.e. the first and higher dervatives of the curves under consideration are not taken into account. Let us assume that the time dervatives R(i)(t)

exist, j

=

1, ... , k . A possible norm is then

(3.4)

IIRII, ..

=

t.

(j

(IIRW(t)ll,)'

dt)

11

These norms give rise to the Sobolev spaces

WP·lc([o,

1]). These spaces are obtained by adding the norm (3.4) to Ck([O, T]), the space of k times continuously differentiable functions on [0,

T]

with continuous k-th derviative.

We note that in practice it has no sense to take for k a big value, because the data con-tain measurement errors, so that higher derivatives are usually unreliable. An appropriate choice is k

=

1. R(1) then follows from

(cf.

(2.3))

(7)

in which y(t,>.) is given by the r.h.s. of (2.1), and y""(t) is calculated numerically from the data.

If one wants to make use of information contained in the highest derivative R(le) only, an alternative norm comes into consideration,

viz.

In the next sections we shall apply (a version of) this norm with p

=

2 and k

=

1. Because R(O) = 0 by definition, this norm then takes the reduced form

T

(3.7) IIRII2

=I

(11R(l)(t)112)2 dt

0

We see, that norm (3.7) only measures the slope of R(t) and not the values of the residuals themselves. In the literature one commonly uses, just as in (2.2) with (2.3), the Lp norm with

p

=

2. It is quite straightforward to establish that the L2 norm is majorized by a modified version of norm (3.7). To show this, we start from the identy

t

(3.8) R(t)

=I

R<1>(r)dr,

0

from which it follows that t

(3.9) IIR(t)112

~I

IIR(1)(r)ll2dr. 0

The Schwarz inequality yields that t

(3.10) CIIR(t)ll2)2

~

t

I

(IIR(1)(r)ll2) 2

dr

0

For the L2 norm in (3.3) thus holds

If we interchange the order of integration in (3.11) we find

T

(8)

I

I •

with the weighting function

w(t)

given by

(3.13)

We note that the r.h.s. of (3.12) is just the norm in (3.7) apart from the weighting function

w(t).

Because

w(t)

>

0 on [0, T] the r.h.s. of (3.12) is also an appropriate norm. The appearance of

w(t)

in (3.12) has a clear interpretation. R<1

>(t)

measures the error in the

slopes of calculated concentration curves at timet. When integrating (2.1) a slope error at time t leads to erroneous function values in the integration interval

[t, T].

So, a slope error near t

=

0 is much more serious than a slope error near

t-

T.

In §§5 and 6 we shall use norm (3.7) in a discretized form. In §7a the weighting function

w(t)

is added as a refinement.

(9)

§4. Optimization

Following [7] we may roughly discern three types of optimization procedures:

1. Ad-hoc methods, e.g. the simplex method. In these methods only function values are compared. They miss a theoretical basis, but may be robust in particular cases. 2. No-derivative methods, e.g. the Powell method. In these methods the Jacobian is not

evaluated. Information about the behaviour of the object function is taken from the function values calculated at earlier iterations.

3. Derivative methods, e.g. (quasi) Newton, Steepest Descent, Levenberg- Marquardt methods. They need at any point in parameter space the Jacobian. Sometimes explicit expressions for the Jacobian are available, but often it has to be calculated numerically. In the latter case extra function evaluations per iteration step are needed.

The ad-hoc methods come not into consideration here, because the problem is too intricate.

Remarks:

i) All methods work iteratively, so an initial quess has to be supplied. This starting vector .Ao should not differ too much from the exact solution, because otherwise all methods men-tioned above will converge to an undesired local minimum. An alternative method, which surmounts this rock, is simulated annealing. This approach is, however, extremely time consuming and completely inadequate in the present context.

ii)

The methods under b) and c) perform best if all variables are scaled to, say, the interval [0, 1] and if the objectfunction more or less varies in the same interval. In the present study, which has to be as general as possible, we do not apply any scaling. In many practical cases one has not the slightest idea about the order of magnitude of some parameters.

(10)

§5. All concentration curves measured

In this section we present a. solution method, which is applicable if all Yi are measured, i.e.

K

=

N. A basic ingredient is to fit the data. Yi( r;),j

=

1, ... , L by means of a smoothing spline. If L is not too small and the measurement errors are not systematic but randomly distributed, the splines will represent the real concentration curves sufficiently well. Splining has always a subjective element in it, because the bending stiffness of the spline can freely be chosen. In the procedure described in [2] and used throughout this report this parameter is set via a customer independent decision criterion. It is strongly recommended to check always the performance of the fit by visual inspection. The spline approximations provide a. continuous fit

ym(t)

to the data. for all

t

>

0. The parameter vector..\ could be estimated by minimizing the function F2 (..\) defined as

N p

(5.1)

F2(A)

=I: I:

Wij (R2)~; i=l i=l

with the residuals R2 given by

M

(5.2)

(R2)ii

ili(t;)- l:Aik(Ym(t;))Ak·

k=l

We note three important points:

i. In the r.h.s. of (5.2) we substitute the measured values

ym(t;)

instead of the calculated values

y(t;,

A).

If the measurements are not exact, this implies that minimizing (5.1) with (5.2) yields a zero order approximation. To improve this one could start an iteration procedure by putting the obtained..\ into IVP (2.1), integrating, subsituting the resulting

y(ti, ..\)

values into the r.h.s. of (3.2), and repeating the minimization procedure. In

practice, however, this iteration procedure is not recommendable, because the slope values

ym(t;)

will never be exactly known.

ii. The object function in (5.1) is a discretized version of the norm in (3.7) provided with weight factors Wii· Because the integration is replaced by the summation it is recommend-able to choose the times ti equidistant. There is no need to identify the time points

t;

in (5.2) with the measurement time points

r;.

The choice of the

t;

in (5.2) can be adjusted if appropriate in the minimization procedure.

iii. The parameters ..\ to be estimated often occur linearly in the residuals. Then, the

estima-tion problem can be calculated via. an explicit expression. If one wants to make explicit use of the linearity, one could apply specific methods as generalized least squares or sin-gular value decomposition. In practice, however, it is equally appropriate to minimize

F2

via. a standard procedure. This is extremely fast, because no integration is involved. In

(11)

§6. Missing data

We assume here that the concentrations of some reactants can not be measured, i.e. K

<

N.

We order the concentrations such that 1/i,

i

=

1, ... , k are measured and 'Jii,

i

=

k

+

1, ... , N,

not. Then, the ODE in (2.1) splits into two parts: (6.1a)

u

= A1(u, v).A

(6.1b)

v

=

A2(u, v).A

with u containing the measured components of y, thus u

=

(y

1 , ... , 1/K )T, and v the remaining

part, thus v

= ('JIK+t, ... ,'JINl·

The matrix A1 contains K rows of A and the matrix A2 the

( N - K) remaining ones. At

t

=

0 both u and v are given. The data are um(

r; ),

j

=

1, ... , L.

The idea is to apply the method in §5 to the system (6.1a), while the needed v(t) curves are obtained from numerical integration of (6.1b). Analogous to (3.1) and (3.2) the object function to be minimized is taken to be

K P

(6.2) Fa(.A) =

2.::

2.::

Wij (Ra)i;

i=l j=l

with Ra given by

M

(6.3) (R3)i:i

=

ui(t;)-

l:(At)ilc(um(t;), v(tJ,.A))-X~c.

lc=l

The time points

t;

can be chosen freely. The um(t) curves are obtained by splining the data. The procedure to evaluate F3(.A) for given

A

is as follows:

a. Substitute the

u,m( t)

curves into the r .h.s. of ( 6.1 b). Integrate ( 6.1 b) numerically starting with the known initial values for v. Denote the result by v(t,.A).

b. Evaluate F3 by sustituting the um(t) and v(t,.A) data into R3 •

The computation time of this approach is mainly determined by step a. In the following section we discuss some modifications of this procedure, which are inspired by its application in practice.

(12)

§7. Refinements

Application of the ideas in §6 motivated us to consider some modifications of the algorithm. We discuss here three modifications. Two of them (7a and 7b) have already been implemented and appeared to be improvements.

7a. Object functions

Contrary

toFt,

the object functions

F

2 and

Fa

are measures for the deviations between the

slopes of the calculated and measured concentration curves. Minimizing

F

2 or

Fa

might lead

to solution curves which increasingly deviate from the data for increasing

t.

As for

F2

in §5 (all concentrations measured) this objection appeares to be of no great importance. However,

if less information is available, the solution method in §6 appears to become sensitive for the detailed choice of the object function. Instead of F4 we therefore propose to use

K p

(7.1) F4(A)

=I: I:

Wi;(R4)t;

i=t .:i=l with the residuals R4 given by

(7.2)

in order to introduce the weighting function derived in (3.8)- (3.13).

7b. Boundary layers

Many systems, e.g. the one in the Appendix, show stiff behaviour. In Figure 1 boundary layers in the Y2, Ya and Ys curves are present for small

t.

These three curves start at zero for t

=

0 and jump to a certain level within a relatively very short time interval. This is usually the case in chemical reactions and gives rise to two difficulties in the current context. First, the numerical integration through the boundary layers consumes unacceptably much computer time. Second, no data are available in the boundary layers, which implies that the smoothing splines urn(

t)

are completely unreliable there. This renders the numerical integration of (6.1.b) using um(t) in fact impossible. To overcome this problem we propose to start the integration of (6.1b) not at t

=

0 but at t

=

t1 , with t1 a point in time beyond

all boundary layers. In practice, we use t1 = r1 , the first measurement time point. Because

for

t

~ ft the concentration curves are quite smooth, the reduction in computer time is considerably. The drawback is that the initial values v1 := v(t1 ) are not known in advance.

Together with the A components they have to be estimated now. Initial values v1 are readily

obtained from the initial estimates for A by integrating (2.1) once from

t

=

0 up

tot=

ft. In this modification the object function depends not only on A but also on v1, thus F

=

F(A, Vt)·

7c. Linearization

The more sophisticated optimization procedures require per iteration step the evaluation of the derivatives of the object function F with respect to the variables. Here, we propose ideas to perform these calculations in an efficient and accurate way. In the first instance, the parameter space is M dimensional with elements A. If we apply the ideas presented in 7b, not only A but also N - K initial values v1 have to be estimated. The corresponding parameter

(13)

space is then M

+

N- J( dimensional with elements(.\, v1)T. The ideas in this subsection are also applicable in that case. To be as general as possible we present the formulae in a form, which assumes that the ideas in §6 en §7b are applied. However, linearization is of course also applicable directly to the problems stated in §2. See, e.g., [8]. As for the notation, we write F1

(.\, v1j 6.\,6v1 ) for the derivative ofF at the point (.\, v1 )T in the direction of

6 := ( 6.\, 6v1l· In practice, 6 is chosen along the axes of the parameter space, so 6 contains in general zeros except at one position.

Optimization procedures based on derivative methods can be divided into two classes. Either

F' is calculated numerically by the routine itself, or the user has to provide a subroutine which yields F' directly. We deal with these two cases successivily.

(i) The most direct way of calculating the directional derivative F' is via the approximation

(€

<<

1)

(7.4) 'F'(~ A 1V1.~, ~ 1 UA,uVt ,....,. )...,F(.\+€8A,vt+€8v1)-F(.\,vl) •

(

It needs the evaluation of v(.\, v1) and v(.\+ fDA, v 1 +<:c5v1), so the number of numerical integrations is M

+

N- J(

+

1. In this notation, v(.\, v1 ) is the solution of the IVP

while v(.\ + <:hA, v1 + <:8v1 ) is the solution of the IVP

{

v

= A2(um, v)(.\+ fDA) (7.6)

v(tt) = Vt

+

€0Vt •

Note, that problem (7.6) is a regular perturbation of problem (7.5). In practice, either the vector field is perturbed via 8.\ (and c5v1

=

0) or the initial value is perturbed via

6v1 (and 6.\

=

0). Now we may apply the regular perturbation theorem [3, §11.6]. If only the vector field is perturbed, the differential quotient z := (v(.\ +fDA, Vt)- v(.\, vt))/£ satisfies for small £ the linear IVP

{

z

= J(um, v(.\, Vt))z+ A2(um, v(.\, Vt))6.X/II6.\II2 (7.7)

z(tt) = 0

.We note that this approach to calculate the derivatives also forms the basis of the work reported in [8].

The Jacobi matrix J is obtained by differentiating the vector field A2(um, v).\ with respect to v. This straightforward exercize could be done by hand or, in case of very complicated systems, via a formula manipulation package.

If only the initial value is perturbed, the difference z := (v(.\, v 1 + €6v1) -v(.\, vt))/£ obeys for small £ the linear IVP

(14)

(ii) If the optimization procedure needs a user supplied subroutine for

F',

the ideas under

(i)

could be elaborated once more. In that case F1 should be calculated via the exact

expression (7.9)

The derivative vector 8F/8v is quite easily obtained for the object functions

Ft, ... ,

F4

introduced above. It remains to determine v1

• From the considerations above it follows

that v' is the solution of either (7.7) or (7.8).

Let us discuss the (dis )ad vantages of the methods proposed in this subsection. A straight-forwardly applied standard procedure calculates

F

1 via (7.4) and integrates both (7.5) and

(7.6) numerically. This has several disadvantages:

- It is tricky to differentiate numerically via approximation (7.4) because of rounding errors. The increment may not be too small, but its value can often not be adjusted in standard routines,

- The accuracy in integrating (7.5) and (7.6) should be an order of magnitude higher than the accuracy in approximation (7.4). This implies the numerical integrator will take small steps, which costs a lot of computer time.

Both problems can be remedied by applying (7.9). Then, the accuracies of taking derivatives and integrating coincide more a less. A second improvement is reached by calculating V1 via

(7.7) and (7.8) instead of using (7.5) and (7.6). Equations (7.7) and (7.8) are linear and in numerical libraries routines are available which exploit this property to obtain more accurate results with less computer effort and thus time.

An intermediate approach might be the combination of (7.4) and (7.7) and (7.8). The com-bination of (7.9) with (7.7) and (7.8), however, is much more consistent.

(15)

§8. Simulations

The methods presented above have been partly implemented and tested using the model in the Appendix. We have chosen to perform the tests via simulations rather than via real data, though these are available, because this approach gives us the freedom to vary the circumstances in a systematic way.

Furthermore, it avoids the need of taking into account the effects of measurement errors, which would obscure the analysis. We generated data by choosing the parameters as mentioned in (A.4). Next, we tried to estimate (part of) the parameters while varying the number of missing concentrations. In all runs we used [2] as splining procedure and [5] as integration procedure. In simulation 1 we applied a no-derivative (viz.

[6]}

and a derivative method (viz. [4]) for the optimization.

Simulation 1. No missing data, 8 parameters.

In this simulation we used the data of all concentrations Yi, i

=

1, ... , 8. Then, the method of §5 is applicable, so object function

F

2(.\) is minimized. As starting vector we took the

real values in (A.4) divided by some factor. Because the problem is linear we may solve it directly. However, to keep comparisons with other simulation runs consistently, we applied the methods in [4) and

[6].

The no-derivative method in

[6]

performed quite badly. If we used as starting vector the (A.4) values divided by a factor of ten, this method found reliable results only if not more than two of the ki and/or Ki in (A.4) were assumed to be unknown. The derivative method in [4] always succeeded to determine the solution exactly, even if eigth parameters, i.e. ki,

i

=

1, ... ,5 and K1 , K3 and K6 in (A.4), are to be estimated. Because

no numerical integration is involved this takes hardly any computer time. As an alternative we also performed the same exercize with F4 in (7.1) instead of F2 in (5.1). This made the

performance considerably worse.

Simulation 2. Missing data, 8 and 5 parameters.

In this case we modified simulation 1 by assuming one concentration, namely Y1, being not measured. As starting values we used the figures in (A.4) divided by a factor of ten. First, eight parameters were assumed to be estimated. The method in §6, extended with the procedure in §7b, yielded no results for both object functions F2 and F4 , i.e. the minimization

procedure stopped after a few iterations. Decreasing the number of unknown parameters, so that only ki,

i

=

1, ... , 5 were estimated, did not improve this result. Choosing the starting values a factor of two closer to the real values also did not help .

.

Simulation 3. Missing data, 3 parameters.

We restricted the set of parameters to be estimated to (K1 , K3 , K6 ). As starting values we

again used the real values, given in (A.4), divided by a factor of ten. We gradually increased the number of missing data by taking (yt),(y11y2),(y11y2,y3 ) etc. to be not measured and

applied the methods in §6 and §7b. The use of

F

2 appeared to be not succesful, whereas

application of F4 led to the appropriate convergence. With (y1 ) and (yt, Y2) missing, the

exact solutions were obtained with 11 and 25 function evaluations resp. With (Yh Y2, Y3) and

(Yl,

Y2, Ys,

Y4)

missing the obtained solutions deviated slightly but appreciately from the real

parameter values. The numbers of function evaluations were 81 and 209 resp. In the case of four missing concentrations, the solution was

(K

1

,K

3

,K

6 ) = (2.50

10-

3,8.2210-3,8.93106

)

(16)

Because we applied the method in §7b, also the values of Y1(10), Y2(10), Y3(10), Y4(10) were estimated. The error in the initial value estimations appeared to be much larger than the error in the parameter estimations. However, if the system (2.1) is integrated from t = 0 to,

t =

T

using the estimated (K1 , K3 , K5 ) values, the agreement with the· experimental data is

splendid. In Fig. 1 the real (solid lines) and estimated (dotted lines) concentration curves are plotted for this case.

Varying the missing concentrations we found that data. for y8 are absolutely necessary to

obtain results. If y8 were assumed to be missing the procedure did not yield reliable results, even if the starting values were closer to the real values.

(17)

§9. Conclusions

i. If the concentrations of all reactants are measured, the method in §5 provides a simple and extremely efficient scheme to estimate the kinetical parameters.

ii. If some concentration curves are not known, one could apply the "classical" method described in §2. This approach takes awkwardly much computer time because it makes use of the numerical integration of the complete system (2.1 ), which is in general stiff. m. Much computer time can be saved using the alternative method described in §6 and in

§7b.

iv. The applicability of the methods in §6 and §7b is shown via simulations with the model in the Appendix. The results depend of course on this model, which is known to be quite intricate. We obtained good results but only if the number of unknown parameters is not too big(~ 3). As for the missing data it makes great difference which concentration curves are not measured. Missing data for (y1 , y2 , y3 , Y4) did not obstruct the satisfactory

estimation of the parameters, whereas missing y8 data appeared to be disastrous in this model.

v. In our simulations no-derivative optimization methods, e.g. as in [6], were quite unsuc-cesful even if applied to rather simple problems. The derivative method in [4] (modified Levenberg-Marquardt) showed satisfactory performance.

vi. The more sophisticated optimization procedures employ the derivatives 8F / 8>.. of the object function F(>..). Their fast and accurate evaluation is essential for the applicability of these methods. Commonly, and also in our simulations, these derivatives are numer-ically determined via (7.4). It is well-known that the choice of the increment in (7.4) is a tedious point. If this is too small, rounding errors creep in, whereas taking it too big also leads to inaccurate results. A complicating factor is that both F(>..) and F(>..

+

6>..) are calculated via numerical integration, the accuracy of which is hard to control. A possible strategy is to increase this accuracy as much as possible, but this enhances the

integration time considerably. ·

A promising idea to overcome this problem is presented in §7c. There, we point out that

8Ff8>.. can be calculated by solving a linear set of ODE's. One then uses (7.9) with (7.7) and (7.8) instead of (7.4). In this approach numerical differentation is avoided. Further-more, the remaining initial value problems are linear which allows the application of special algorithms, which are accurate and fast. The idea could be used in combina-tion with the "classical" method in §2. Examples concerning quite simple systems are worked out in [8]. For complicated system its application in the context of the methods presented in §6 seems to be more attractive and certainly deserves further research to determine its merits in practice.

(18)

References

1. K.K. Hong, Parameter Estimation in Chemical Reactions, final report of the postgraduate programme Mathematics for Industry, Eindhoven University of Technology, may 1990. 2. RENELA, Fortran subroutine described in THE-RC64650, inf PP-4.18 and in Computing

Centre Note 19 by C.J.J.M. van Ginneken of the Computing Centre of the EUT, march 1983.

3. R.M.M. Mattheij and J. Molenaar, Beginwaardeproblemen in theorie en praktijk, analyse, numerieke methode, modellen, Epsilon uitgave nr. 19, 1991, ISBN 90-5041-025-1.

4. LMDIF, Fortran subroutine, MINPACK library, Argonne National Laboratory, B.S. Gar-bow, K.E. Hillstrom, J.J. More, 1980.

5. LSODE, Fortran subroutine, ODEPACK, a systematized collection of ODE solvers, in Scientific Computing, R.S. Stepleman et al. (ed.s), North-Holland, Amsterdam 1983, pp. 55-64.

6. E04F DF, Fortran rubroutine, NAG workstation library, see "Algorithms for the solution of the non-linear least-squares problem", by P.E. Gill and W. Murray, SIAM J. Numer. Anal., 15, pp 977-992, 1978.

7. R. Fletcher, Practical Methods of Optimization, John Wiley & Sons, Chicester etc., 1987. 8. P.W. Hemker, Parameter estimation in non-linear differential equations, Mathematical

(19)

Appendix

For test purposes we consider the following reaction in which 8 reagents are involved. We denote the reagents by

Yi

and the corresponding concentrations by Yi·

Yt +2Ys Y3 +2Ys Y3+Y4 Y6 +2Ys y2 +-+ y3 +-+

y4

+-+ Ys

+

Y6 +-+

y7

+-+ 2Ys

The corresponding equations are in dimensionless form:

Yt

=

-r1

Y2

= -rs

+

f(ya, Ys)

iJa

= r1 - ra - ra

Y4

= ra- ra

ils

= r3

+

g(ya, Ys)

Y6 = rs- r4

iJ1

=

r4

ils

=

2(rs- r1 - ra- r4)/co with the reaction rates ri given by

rt = kt(YtY~ r2 = k2(YaYi r3 = k3(YaY4 r4 = k4(Y6Y~ rs

=

ks(Y2 -J(tY3) -/C2y4) -J(3Y5Y6) -J(4Y7) -J(sy~) and

co

a given constant.

All initial values Yi(O) vanish except for y1(0) and Ys(O). The terms j(y2, Ys) and g(y2, Ys) describe that Y2 and y5 may be supplied during the process. The explicit forms of

f

and g

are not of interest here. In the first instance, 10 parameters ki, i

=

1, ... , 5, and J(i, i

=

1, ... , 5 are to be estimated. After some time the system reaches equilibrium, from which it can be deduced that ](2 = ](4 = 0.

We use this system for simulations. To that end we set the parameters and initial values at certain values and integrate numerically. The solution curves Yi(

t),

i

=

1, ... , 8 are considered to represent measured data. We have not tried to introduce measurement errors by adding random noise, because the emphasis of this study is to focus on the optimization methods. In the simulations we use as "real" values of the parameters the set

kt = 6.67 101 *Co k2 = 1.04 105

*

c0 ks = 3.87 10-2

*

c0 k4 = 1.14 101

*

Co k5

=

5.30 10-3

*

c0 J(l

=

2.84 10-3 ](2

=

0 ](3

=

6. 70 10-3 ](4

=

0 ](5

=

8.91 106

(20)

3[3 ~ .§ ~

!

c 0 u

20 •a 61 Tlld 0.9 8.8 \ 8.7· 8.6 ~ ' .~ 8.5 ~ ~

...

c 0 u 9.3 9.2 8.1

8 28 48 68 TIJd 2a •a &a TIJd

••

59

---••

-

-... > ,' ' ' .§ ' '

;

30 ' ' ' ~ , , u , I c , 8 , 29 , , , , , , , 19 ' ' 8a 1BB 88 188 88 1BB ...

---·

f! -~

.

L ~

-g 8 ~ ~ ~ ~ g 8

,.

5 0 , . _

••

~ .§ ~ 30 ~ ~ 21! 10

3[3 2598 "' > -~

.

L

1

3E-3 a 8.8825 1[-3 5[-4 , ... ... ... 2a •a 61 8a 1BB TIJd 28 48 68 88 10a Tlld 21! 48 68 ee 191 Tlld

Referenties

GERELATEERDE DOCUMENTEN

The remarkable progress in understanding of various rates entering the kinetic equations describing the asymmetry generation along with considerable improvements of the

M.P. Hagenzieker; SWOV Institute for Road Safety Research, Leidschendam, The Netherlands Annex XII: Bibliography.. Assumptions used in road design M. Slop. SWOV Institute for Road

The L´ evy processes considered in Chapter 4 are the Generalized Hyperbolic process, the Variance Gamma process, the Normal Inverse Gaussian process and the Meixner process....

extensive human presence and interaction as chicks did not influence carcass attributes, meat quality or skin traits 34.. at slaughter age, but more importantly, it did not

This ap- proach was first explored in [ 19 ], where two Newton-type methods were proposed, and combines and extends ideas stemming from the literature on merit functions for

Om die proses verder te verstaan word die volgende vraag gevra: “Het die kultuur en invloede wat identiteit in die verlede gevorm het, verdwyn?” Is dit waar dat die tyd waarin

More specifically, we propose an adaptive formula for the Levenberg–Marquardt parameter and analyse the local convergence of the method under H¨ older metric subregularity of

Instead, the data support an alternative scaling flow for which the conductivity at the Dirac point increases logarithmically with sample size in the absence of intervalley scattering