Caracotsios

Prévia do material em texto

Computers & Chemical Engineering Vol. 9. No. 4. pp. 359-365. 1985 009s1354/85 $3.00 + .a0 
Printed I” Great Britain. Pcrgamon Press Ltd. 
SENSITIVITY ANALYSIS OF INITIAL VALUE PROBLEMS 
WITH MIXED ODES AND ALGEBRAIC EQUATIONS 
MAKIS CARACOTSIOS and WARREN E. STEWART? 
Department of Chemical Engineering, University of Wisconsin, Madison, WI 53706, U.S.A. 
(Received 2 May 1984; revision received 1 October 1984; receivedfor publication 23 October 1984) 
Abstract-An efficient method is described for sensitivity analysis of nonlinear initial value 
problems, which may include algebraic equations as well as ordinary differential equations. 
The linearity of the sensitivity equations is utilized to solve them directly via the local 
Jacobian of the state equations. The method is implemented with the implicit integrator DASSL 
and is demonstrated on a stiff industrial reaction model. 
Scope-With the rapid development of digital computers, increasingly realistic mathematical 
models are being used to investigate chemical phenomena. New mechanistic features, however, 
call for new physicochemical parameters whose values may not be accurately known. Conse- 
quently, there is an increasing need for parametric sensitivity analysis of proposed differential 
and algebraic models. 
Parametric sensitivity analysis is a very active research area. Extensive reviews can be found 
in Rabitz et al [l] and in Tilden et al [2]. Applications occur in every engineering and scientific 
discipline. Potential areas of application in chemical engineering include optimization, parame- 
ter estimation, model simplification, process sensitivity and multiplicity and experimental 
design. 
In this article we address the problem of numerical computation of sensitivity functions for 
systems of ordinary differential and algebraic equations. We develop a simple, efficient 
algorithm for this purpose by extension of a standard implicit integrator. 
Conclusions and Significance-An efficient method has been developed here for the calculation 
of parametric sensitivities for mixed systems of ordinary differential and algebraic equations. 
The method has been tested successfully on many problems, ranging from simple ordinary 
differential equations to large stiff systems of differential and algebraic equations. The scheme 
has been implemented for a robust implicit integrator (DASSL) and is applicable to any implicit 
integrator. 
PROBLEM STATEMENT 
Consider a dynamic system, described by the following 
set of differential and algebraic equations: 
Eu’(t) = f(t, u(t); 0) (la) 
u(t = to) = u,(8). (lb) 
Here u is an n-dimensional vector of state variables, 8 
is an m-dimensional vector of time-independent 
parameters and E is an (n, n) matrix of constant 
coefficients. Most frequently in chemical kinetics cal- 
culations, the matrix E assumes the form 
E= 
TAuthor to whom correspondence should be addressed. 
where I”’ is the identity matrix of order s. Ifs = n, the 
system in eqn (la) consists of purely differential 
equations. If 1 5 s < n, the system in eqn (la) consists 
of ordinary differential and algebraic equations. The 
latter case arises, for instance, in analysis of reaction 
schemes where equilibria give algebraic constraints on 
the concentrations. 
We define the (n, m) matrix W(t) of sensitivity 
functions as 
au(t) 
W(t): = ---g (3) 
This matrix satisfies a set of differential/algebraic 
equations that can be derived by partial differentiation 
of eqns (la) and (1 b) with respect to the parameter 
vector 9: 
EW’(t) - J(r)W(t) = $ f(t, u(t); 63) (da) 
au, (0) 
W(t = to) = 7 (4b) 
CACE 9:4-D 
359 
360 MAKIS CARACOTSIOS and WARREN E. STEWART 
where the matrix J(t) (shorthand for J(t, u(t); e)) is 
defined as 
J(t): = t f(t, u(t); e). 
Various properties of eqns (4a) and (4b) are described 
in Tomovic and Vukobratovic [3]. The most striking 
feature of these sensitivity equations is that they are 
linear, regardless of the linearity or nonlinearity of the 
state eqns (1 a) and (1 b). The problem studied here is 
the numerical computation of the matrix W(t) from 
eqns (4a) and (4b). 
LITERATURE REVIEW ANDTHEORETICALBACKCROUND 
Before describing the new algorithm, we review a 
few known facts about the solution of mixed systems of 
ordinary differential and algebraic equations. Several 
investigators [4-61 have considered this subject, and 
recently Petzold [7] published an algorithm called 
DASSL for the solution of such systems, 
Not all systems of differential/algebraic equations 
are solvable. The reader is encouraged to consult, in 
the literature, Petzold and Gear [8] and Campbell and, 
Petzold [9] on this peculiar feature of mixed systems. 
However, for the systems that we are considering, 
where the matrix E assumes the form shown in eqn 
(2), sufficient conditions are known [8] for the solva- 
bility of eqns (la) and (1 b). 
Let the function f(t, u(t); e) be continuously differ- 
entiable with respect to u(t). Now consider the Jaco- 
bian matrix J(t) defined by eqn (5) for the system 
shown in eqn (la). If we partition the Jacobian matrix 
according to the partition of E, i.e. 
J(t) == . . 
where J,,(t) is an (s, s) matrix, then the system in eqns 
(1 a) and (1 b) is solvable if 
det J**(f) # 0 for all t. (7) 
Under this condition, the solution obtained by a k-step 
backward differentiation formula algorithm with k < 
7 and fixed step size h converges to O(h’) if all initial 
values are correct. Further aspects of eqns (la) and 
(lb) and their numerical treatment are discussed in 
Petzold [lo]. 
I et us now review some of the methods used for the 
computation of the sensitivity matrix W(t). With one 
exception (Stewart and S$rensen [l l]), the known 
methods are for systems of ordinary differential equa- 
tions (ODES) only; that is, for systems with E = I’“‘. 
The available algorithms include the Fourier ampli- 
tude test [ 121, direct differential methods [ 131, 
Green’s function methods [ 141, the analytically inte- 
grated Magnus method [l] (a modification of the 
Green’s function method) and finite difference meth- 
ods. The Green’s function method and its variations 
exploit the fact that the sensitivity equations are linear 
inhomogeneous with time-varying coefficients; conse- 
quently they can be solved by first calculating the 
solution of the homogeneous part and then deter- 
mining the particular solution corresponding to each 
parameter. 
Several authors have proposed solving the sensitiv- 
ity equations by extending known solution algorithms 
for the state equations. This idea is based on the 
identity of the coefficient matrices in the sensitivity 
equations with those in the locally linearized form of 
the state equations on the locus u(t). Stewart and 
S$rensen [I I], Vemuri and Raefsky [ 151, Lojek [ 161 
and Hwang [ 171 have developed various aspects of this 
method; nevertheless the idea is still under develop- 
ment. 
In the present work, we fully exploit the similarity 
of the sensitivity and state equations by building the 
sensitivity analysis into a robust differential/algebraic 
equation solver. Then we illustrate the algorithm by 
solving a stiff industrial kinetics problem. 
MATHEMATICAL DEVELOPMENT OFTHESENSITIVITY 
ANALYSIS ALGORITHM 
One of the most important steps in developing the 
sensitivity analysis algorithm was the selection of the 
integrator. For mixed systems of differential and 
algebraic equations, there are several codes [4-71 
designed to perform the integration. These codes are 
primarily based on an idea developed by Gear [4]; 
specifically, the derivative u’(t,, ,) is approximated by 
a backward difference formula with adaptable order 
and step size, and the resulting system of nonlinear 
equations is solved for u(t,, ,) via a modified Newton 
scheme. 
We chose for this work the package DASSL devel- 
oped by Petzold [7]. This package is portable, robust 
and easy to use. We tested the code successfully on a 
wide variety of stiff problems,both differential and 
mixed, before adding the sensitivity analysis algo- 
rithm. 
First, consider the solution of state and sensitivity 
eqs (1 a), (1 b), (4a) and (4b) as a single system. In this 
approach, one needs the Jacobian matrix of the total 
system in eqns (la) and (4a). If we partition the 
sensitivity matrix W(t) into column vectors as 
W(r) = [w,(t) I W?(t) I I W,,(f)1 (8) 
where 
‘“(‘) i = 1 2 .m w,(t): = - 
dQ, ’ ’ ’ 
(9) 
Sensitivity analysis of initial value problems 361 
then the Jacobian matrix J*(t) [shorthand for J* (t, 
u(t), W(t); e)] of the total system in eqns (la) and 
(4a) is 
J*(t) = 
where 
J(t) 0 0 . . . . . -. 0 
j,(t) J(t) 0 - - - - - - - 0 
j,(t) 0 J(t) 0 
. . . 
. . . 
. . . 
. . . 
j,,,(t) 0 0 . - - . . - . J(t) 
(10) 
ii(t): = f$ aJ(r) w,(t) + F i= i 2 . .m. . 1 (11) 
“, 
The evaluation of J*(t) is a formidable calculation; 
however, it is a natural requirement of Newton itera- 
tion on the total equation system; eqns (la) and (4a). 
A simpler and quicker approach is to solve eqn (la) 
before eqn (4a) at each time step, as shown below; then 
the matrices j,(t) are not required. 
Let i”‘(t) be the local interpolant of u(t) obtained 
in r Newton iterations of a kth-order integrator within 
a given time step. Then the next iteration will give the 
interpolant i”‘(t) + AC”‘(t), which satisfies the fol- 
lowing linearized form of eqn (la): 
E(ii’(“(t) + Aii”“(t)) = f(t, tic’) (t); 9) 
+ J(t, ii’“(t); 8)AO”’ + O(hL). (12) 
Hence the correction Ai?” satisfies 
EAu”” - J”‘(t)A,P”’ = f(t, ii”‘(t); 8) 
- E@‘(t) + O(h’) (13) 
when the standard Newton method [with J”‘(t) 
updated for each iteration] is used. If A@)(t) con- 
verges toward zero with increasing r, then J”‘(t) 
converges toward J(t), and eqn (13) becomes formally 
similar to eqn (4a). Therefore, we can defer consider- 
ation of eqn (4a) until P”‘(t) has converged to i(t) at 
the current value of t. Then we can update the 
sensitivity solution W(t) directly by the use of eqn 
(4a), which has the same coefficients as eqn (13) but a 
different, now computable, right-hand function. More 
specifically, the corrections AC,(t) are computed in a 
single iteration by solving 
EAG;(t) - J(t) A%,(t) == -Eti;‘P’(t) + J(t) %-+‘)(t) 
in which GiP’(t) is the predicted value of 9,(t) via a 
kth-order predictor formula. 
The vectors Ei(t) can be calculated at the current t 
as follows: 
G,(t) = iv!@(t) + Ait, I’ = 1, 2,. . m. (15) 
On completion of the update, the local truncation error 
is tested and, if necessary, the step size h and approxi- 
mation order k are adjusted to achieve the specified 
accuracy for u(t) and W(t). 
Numerical tests show that stringent tolerances on 
u(t) normally lead to a good solution for W(t) as well, 
provided J(t) is updated before W(t) is computed. 
On the other hand, if the iteration matrix is only 
occasionally updated (as is usual in implicit integra- 
tors), then the local tolerances on W(t) are essential to 
control the calculation. 
COMPUTER IMPLEMENTATION 
The computer implementation of the sensitivity 
analysis algorithm was done as follows: 
The working arrays used by DASSL were modi- 
fied to provide storage allocation for the sensitiv- 
ity functions. 
An algorithm was written for the automatic 
formation of the sensitivity equations. 
The integrator in DASSL was properly extended 
to include the solution for the sensitivities. 
The following definitions were adopted: 
1. lteration matrix D(t): = cE - J(t), where c is a 
constant that depends on step-size history. 
2. Residual of the state equations: 
R,(t): = Eu’(t) - f(t, u(t); 0). 
3. Residual of the sensitivity equations: 
Ri(t): = E[w;@‘(t) - ~wi’~‘(t)] 
& f(t, u(t); 0) i= 1,2,...m. 
The calculation sequence is outlined in the adjoining 
flow chart. 
NUMERICAL EXAMPLE 
The algorithm was tested on a wide variety of 
problems, ranging from linear and nonlinear systems 
of differential equations to large systems of differen- 
tial and algebraic equations. We present here a batch 
reactor example given by the Dow Chemical Company 
[ 181. Figure 1 shows the proposed mechanism for the 
reaction system. The time-dependent concentrations 
are modeled by the following system of differential 
and algebraic equations: 
u;(t) = -k,u,(t)u,(t) 
u;(t) = -k,u2(t)ue(f) + k_,u,,(t) 
+ &f(t. ti(t);e) -t O(P) i = 1,2,. .m (14) 
“Vi - k+,(t) ug(t) 
362 MAKIS CARACOTSIOS and WARREN E. STEWART 
lnilialize arrays for 
new time step 
Compute coefficients for 
backward differentiation 
formulas 
Predict ~ltl,~~tl.~~tl.~‘(t) 
SOIVS D”‘(1) &‘* -6 
I _., .,0(t) 
I 
Convsrpsncs ? 
NO 
6 P 
FLOW CHART C 
FLOW CHART A 
Slow Kinetic React ions 
0 5 
13 Solve o(t) ii = -Fii (11 
Estimate local truncation error 
i 
for u(t).w(t) 
_ w 
%4- + 20M 
‘A- + BM 
M- + 4AB 
k1 - IoMBM- 
k-l 
k2 -‘ABM- 
k3 _ 
k-3 
ABM- 
Rapid Acid-Base React ions 
5MBMH 
KI _ 
MBM- + 
‘HA 
K2 ~ 
A- + 
3HABM K3 - ABM- + 
7HS 
H+ 
H+ 
Fig. I. Chemical reaction model for the numerical example 
[IS]. The chemical species are numbered as in eqn (16). FLOW CHART B 
Sensitivity analysis of initial value problems 363 
u;(t) = -k,u,(t)u,(t) + k~,u,(t) 
u;(t) = k,uz(r)u,(t) ~ k ,u,,(t) 
d(f) = -k,u,(t)u,(t) ~ k+q(f)ug(r) 
+ k-,u,o(t) + k A(Z) 
u,(t) = -0.0131 + u,(t) + u,(r) + U,(f) 
+ u,o(t) 
43(t) = 
Kzu,(t) 
K, + u,(r) 
u,(t) = 
K,u,(t) 
K, + u,(t) 
ho(t) = 
K, us(t) 
4 + u,(t) 
with initial conditions 
u,(O) = 1.5776 
u*(O) = 8.32 
u,(O) = 0 j = 3, 4, 5,9, 10 
&(O) = 0.013 1 
u,(O) = 0.5 i--K, + JK; + 4K,u,“‘} 
u,(O) = u,(O). 
(16) 
(17) 
The following values of rate and equilibrium constants 
were used [ 191: 
k, = 21.893 hr ’ Kg gmole ’ 
k_, = 2.14 E09 hr-’ 
k, = 32.318 hr ’ Kg gmole ’ 
k, = 21.893 hr-’ Kg gmole ’ 
k 3 = 1.07 E09 hr-’ 
K, = 7.65 E-18 gmole Kg ’ 
K, = 4.03 E-l I gmole Kg ’ 
K, = 5.32 E-18 gmole Kg ‘. (18) 
0 1 1’- tlrs3 4 5 6 
Fig. 2. Concentration functions u,(r) in gmole/kg. 
0.015 ~ 
AA- _ _..._.- 
6 /_.C.---- - 
,/ 
,/” 
...--R 
- / .\ 
0.00 5- ,/’ ‘l 
“‘..< 
--A. A- 
--+__ 
0.000, 
0.0 
, . I . ( , , I . , I I , I , I , I I 
0.5 I.0 1.5 2.0 
t-hrs 
Fig. 3. Concentration functions u,(t) in gmole/kg. 
kl 
1 
-0.66 
t-hrs 
Fig. 4. Elements of the sensitivity vector &s(t)/alnk,. The 
number I denotes the sensitivity coefficient u,(t) with respect 
to Ink,. Numbers not shown correspond to sensitivity func- 
tions indistinguishable from zero at all times. 
Table 1. Computational effort for the solution of the Dow problem 
State equations Sensitivity equations 
et = l.OE-6 c = l.OE-7 P = l.OE-6 e _ l.OE-7 
Time steps I91 265 185 283 
Function evaluations 418 587 377 571 
CPU set 4.0 5.5 18 27 
‘I’Local error test at each time step which requires that abs (local error) not exceed rtol . abs (u) + atol 
where rtol = atol = e. 
364 MAKIS CARACOTSIOS and WARREN E. STEWART 
The natural logarithms of these constants make up the 
parameter vector 8. Our goal is to estimate the sensi- 
tivity matrix W(t) == &i(t)/*. 
The combined system of state and sensitivity func- 
tions consists of 90 equations, 54 of which are differen- 
tial and 36 are algebraic. This problem presents a 
severe test for the DASSL integrator and the sensitiv- 
ity analysis algorithm. 
Table 1 summarizes the computational effort for 
the solution of the above problem on a VAX I l/780 
computer. All calculations were carried out in double 
precision. A mixed local truncation error control pro- 
vided by Petzold [7] was used. The tolerances for the 
sensitivities were equated with the tolerances of the 
state variables. The total reaction time considered was 
53 hr. 
Figures 2 and 3 show the evolution of the concentra- 
tion profiles as functions of time, while Figs 4 through 
11 show the dynamic behavior of the sensitivity func- 
tions. 
-0.4.0 -“.4~ IO 
t-hrs t-hrs 
Fig. 5. Elements of the sensitivity vector du(t)/dlnk_, Fig. 9. Elements of the sensitivity vector du(t)/dlnK,k2 
-0.75 I ’ , ’ I ’ I ’ 
0 2 t4hrs 6 6 IO 
Fig. 6. Elements of the sensitivity vector &~(t)/dlnk,. 
0.4 
k3 
Fig. 7. Elements of the sensitivity vector du(t)/~Ylnk,. 
Sensitivity plots, like those in Figs 4 to 11, can be of 
considerable use to the theoretician as well as to the 
experimentalist. From Figs 4 through I1 we see that 
all of the rate and equilibrium constants have compa- 
rable effects on the concentrations. and therefore we 
-0.4 I * , * , I , I 
0 2 
, I 
t!hrs 6 8 10 
Fig. 8. Elements of the sensitivity vector ~?u(t)/dlnk_, 
o.4iY-TT7+0 
t-hrs 
Fig. IO. Elements of the sensitivity vector du(t)/alnK,. 
0.4 
0.2- 
*0.0 
K3 dK, _ 
-0.2- 
-0.4.0 
t-hrs 
Fig. I I. Elements of the sensitivity vector du(t)/alnK,. 
Sensitivity analysis of initial value problems 365 
cannot eliminate any step from the proposed mecha- 
nism in Fig. 1. However, a close inspection of the 
results reveals the following linear relations: 
k au(t)=K_ au(t) 
-’ ak., ’ dK, 
k au(t) _ _K au(t) 
-3 ak-~, ’ aK, 
(19) 
(20) 
Thus one cannot estimate all parameters of the model 
from data on u(t) for the range of conditions shown. 
This information proved useful in fitting the model in 
eqns ( 16) and (17) to the published experimental data 
for this system [ 191. 
The software package for the present algorithm is 
available from the authors. 
AcknowledgPmenr~This work was supported by the 
National Science Foundation through Grant CPE-8308748 
and by the United States Army under Contract DAAG29- 
80-C-0041. We also express our appreciation to Linda Pet- 
zold for helpful discussions regarding the use of her integra- 
tion package. DASSL. 
I. 
2. 
3. 
4. 
5. 
REFERENCES 
H. Rabitz, M. Kramer and D. Dacol, Sensitivity analysis 
in chemical kinetics. .4nn. Rev. Phys. Chem. 34,419-461 
(1983). 
J. W. Tilden, V. Constanza, G. J. McRae and J. H. 
Seinfeld. In: Modeling of Chemical Reaction Systems, 
K. H. Ebert, P. Deuflhard and W. Jager, eds. Springer- 
Verlag, Berlin, pp. 69-91 (1981). 
R. Tomovic and M. Vukobratovic, General Sensitivity 
Theory. American Elsevier, New York (I 972). 
C. W. Gear, Simultaneous numerical solution of differ- 
ential-algebraic equations. lEEE Trans. Circuit Theory 
CT-18. 89-95 (1971). 
A. C. Hindmarsh, LSODE and LSODI, two initial value 
ordinary differential equation solvers. ACM SIGNUM 
Newsletter 15, IO-1 1 (1980). 
6. 
7. 
8. 
9. 
IO. 
II. 
12. 
13. 
14. 
15. 
16. 
17. 
18. 
19. 
J. W. Starrier, A numerical algorithm for the solution of 
implicit algebraic-differential systems of equations. 
Technical Report 3 18. Department of Mathematics and 
Statistics, University of New Mexico (1976). 
L. R. Petzold, A description of DASSL: a differential/ 
algebraic system solver. Sandia Tech. Rep. 82-8637 
(1982). 
L. R. Petzold and C. W. Gear, ODE methods for the 
solution of differential/algebraic systems. Sandia Tech. 
Rep. 82-8051 (1982). 
S. L. Campbell and L. R. Petzold, Canonical forms and 
solvable singular systems of differential equations. San- 
dia Tech. Rep. 82-8891 (1982). 
L. R. Petzolh. Differential/algebraic equations are not 
ODE’s, SIAM J. Sci. Siat. Comout. 3. 367-384 (19821. 
W. E. Stewart and J. P. S&ensen; Sensitivity and 
regression of multicomponent reactor models. Fourth 
International Symposium on Chemical Reaction Engi- 
neering, DECHEMA, Frankfurt, I-12 to I-20 (1976). 
R. I. Cukier, H. B. Levine and K. E. Shuler, Nonlinear 
sensitivity analysis for multiparameter model systems. J. 
Comput. Phys. 26, 1-42 (I 978). 
R. P. Dickinson and R. J. Gelinas, Sensitivity analysis of 
ordinary differential equations-a direct method. J. 
Compuf. Phys. 21, 1233143 (1976). 
J. T. Hwang, E. P. Dougherty, S. Rabitz and H. Rabitz, 
The Green’s function method of sensitivity analysis in 
chemical kinetics. J. Chem. Phys. 69, 5180-5191 
(1978). 
V. Vemuri and A. Raefsky, On a new approach to 
parameter estimation by the method of sensitivity func- 
tions. Inf. J. Syst. Sci. 10, 3955407 (1979). 
B. Lojek, Sensitivity analysis of nonlinear circuits. IEEE 
Proc. 1296(3), 85-88 (1982). 
J. T. Hwang. Sensitivity analysis in chemical kinetics by 
the method of polynomial approximations. In?. J. Chem. 
Kin&. IS,3955407 (1983). 
G. Blau, L. Kirkby and M. Marks, An industrial kinetics 
problem for testing nonlinear parameter estimation algo- 
rithms. Process Math Modeling Department, The Dow 
Chemical Company (1981). 
M. Caracotsios, W. E. Stewart and J. P. Sdrensen, 
Nonlinear parameter estimation for an industrial kinet- 
ics problem. In preparation.