Getting around consistent initialization of DAE systems?

Computers chem. Engng Vol.2 I, No. 2, pp. 145-158, 1997 Copyright© 1996 ElsevierScienceLtd Printed in Great Britain.All rights reserved

Pergamon 0098-1354(95)00264-2

0098-1354/97 $15.00 + 0.00

GETTING AROUND CONSISTENT INITIALIZATION OF DAE SYSTEMS? A. KRONER',W. MARQUARDTb and E. D. GILLESa " Universit?it Stuttgart, lnstitutfur Systemdynamik und Regelungstechnik, Pfaffenwaldring 9, D-70569 Stuttgart, Germany b RWTH Aachen, Lehrstuhlfiir Prozefltechnik, Turmstr. 46. D-52064 Aachen, Germany (Received 29 September 1994; revised 6 October 1995) A b s t r a c t - - T h e initialization of DAE systems in general is a difficult numerical problem. If step responses in the vicinity of some steady state are studied, rigorous numerical initialization can be bypassed under certain restrictive conditions. For that purpose a sufficiently smooth input function approximating the discontinuous step is applied to the steady state and integration with some automatic integration code is started without initialization. This contribution analyzes the restrictions in applying this method successfully as well as the potential shortcomings of the computed approximate solutions. Copyright © 1996 Elsevier Science Ltd

INTRODUCTION

comprising n underlying ODEs (3) and q algebraic equations (4) can be obtained from (2) by symbolic manipulation according to a procedure discussed by Gear (1990), Bachmann et al. (1990) and Unger et al. (1995). The algebraic equations (4) are invariants to (3) and define a n - q dimensional manifold which the solution of the problem is restricted to. The right hand sides • and ~ may depend on first and higher order time derivatives of the inputs due to differentiations of the equations, To start numerical integration a set of initial values

A growing number of integration codes have become available (Brenan et al., 1989; Deuflhard et al., 1987; Hairer and Wanner, 1991) in the past for solving DAE systems

f(Lx,u)=0,

(i)

with f, i , x E R", af/0x being singular. The variables u R p are either constant parameters or sufficiently smooth input functions depending explicitly on time t. The major characteristic of DAE systems is their differential index ~, which we define, according to Brenan et al. (I 989), as the number of differentiations of all or some of the equations in (I), such that the consistency equations

f(x,x,u)=0 ~t f(£,x,u)=O = ~ Ofi + •

Of --~+

Of •

Ox

Ouu u

(2) d" Of x~, j~+ dt f(Lx,u)=0= ~ ~ ...

can be solved for £ as a function of x and t, at least locally. The corresponding extended system

t =~(x,u,it,...)

(3)

0 = q~(x,u,fi,...),

(4)

Current address: Linde AG, Dr.-Carl-von- Linde-Str. 6-14, D82049 Hrllrieglskreuth, Germany. * Author to whom all correspondence should be addressed. Phone: +49 241 806712, Fax: +49 241 8888-326, e-mail: marquardt @1fpt.rwth-aachen.de,

x(to) = Xo, ~l(to)= ~

(4)

must be provided. If initial values are chosen arbitrarily or such that they only satisfy equation (1) they are in general not a solution of the extended system (3), (4), i.e. they do not tie on the solution manifold and are therefore called inconsistent. Initial values which are a solution of the DAE systems (1) and satisfy the extended system are said to be consistent. Often, their determination turns out to be a nontrivial task. Inconsistent initial values always lead to discontinuities in the solution during the first integration step(s) (see Sincovec et al., 1981). This property can easily be understood from the analytical solution of linear constant coefficient DAE systems as presented in Brenan et al. (1989). Numerical solution with fixed stepsize integration using for example an implicit Euler method is feasible even with inconsistent initial conditions (Sincovec et al., 1981). This approach cannot be recommended, however, because of the lack of efficiency and error control. State of the art automatic DAE integration methods should rather be used employing step size and order adaptation to bound the local error and to minimize the computational effort at the same time. If inconsistent initial values are used to start

145

146

A. KRONERet al.

integration with standard automatic DAE integration Both problems seem to be overcome if not a step change codes like DASSL (Brenan et al., 1989), RADAU5 but a smoother forcing function continuous at t = to is (Hairer and Wanner, 199 I) or LIMEX (Deuflhard et al., applied to the steady state of the process model (see 1987) the various stepsize and order selection strategies Simeon et al. (1991), Yu (1991), or Rovaglio et al. normally fail due to this discontinuity. The local error in (1992) for examples). the first step contains a contribution which depends on In order to compute step responses of a differentialthe distance between the inconsistent initial values and algebraic process model one may first compute the the solution manifold. This error portion cannot be steady state xs directly by solving the algebraic equation controlled by conventional stepsize and order reduction system policies. The local error may either tend to a constant as 0=f(0,xs,us). (6) stepsize h ---*0, or it may grow infinitely if it depends on negative powers of the stepsize. for given u~ by means of an iterative method and next Hence, consistent initialization is crucial for a suc- apply a continuous forcing function which approximates cessful solution of dynamic process simulation problems the step reasonably well to start integration without and turns out to be a major hurdle in large-scale practical rigorous initialization. This approach to DAE initializaapplications. Compared to the work on solving DAE tion is subsequently referred to as the "direct method" systems with index one and higher much less effort has which is obviously only applicable if the dynamics in the been spent on computing consistent initial values. First, vicinity of a steady state needs to be investigated. Sincovec et al. (1981) investigated the numerical This method may sometimes be appropriate but solution together with consistent initialization of higher neither guarantees a successful start of the code nor a index linear constant coefficient DAE systems. Based on correct solution in general. This is due to the fact, that a structural considerations Pantelides (1988) determined successful start of automatic DAE integration codes the number of differentiations of each equation in (1) from a steady state requires in general the solution x(t), necessary to allow determination of consistent initial t -> to to be at least continuous in the first step. As values. Purely numerical approaches to initializing illustrated subsequently, this can often only be achieved nonlinear DAE systems are presented by Berzins et al. if forcing functions applied to inputs are higher order (1989), Leimkuhler et al. (1991), and Simeon et al. differentiable. Typically, the differentiability of some (1991). Kr6ner et al. (1992) enhanced the efficiency of input functions must be of a higher order than that of Leimkuhler's numerical approach incorporating struc- others. tural information about the DAE systems available from In this contribution, we will investigate conditions a structural analysis (Pantelides, 1988; Unger and under which the direct method for initialization of DAE Marquardt, 1991). systems is appropriate. We first analyze whether the The rigorous techniques for consistent initialization steady state solution as computed from (6) can serve as are in general difficult to apply. Even worse, since they consistent initial conditions for a DAE of arbitrary always require the solution of a set of nonlinear index. Then, the required degree of smoothness of the algebraic equations, they may fail due to convergence forcing functions is investigated in order to guarantee a problems (e.g. Majer et al., 1995). Due to these continuous solution in particular in the first time step. difficulties simpler approaches to initialization are of This is a necessary condition for starting integration with considerable interest in computational practice. most automatic DAE solvers. In addition, fundamental Rather than investigating process dynamics for arbi- limitations of smoothing are considered for index one trary initial conditions and time varying forcing func- DAE systems. The findings will be illustrated by a tions, step responses in the vicinity of a steady state are number of index one and index two examples. of prime importance in engineering analysis. Hence, the It should be noted that our results not only apply to process being at steady state initialization for step response computations but in ~ = 0 , x,=const, u~=const (5) analogous way to reinitialization after a discontinuity during transients. is perturbed by a stepwise change of one or more inputs u(t) = u~ + Au, t --> to. Obviously, such a step change in ARE STEADYSTATESALWAYSCONSISTENTINITIAL u leads to inconsistent initial values since neither the CONDITIONS? original problem (i) nor the extended system (3), (4) is fulfilled at to+. Further, as shown in (3), (4) a certain It has been taken for granted previously that steady degree of smoothness of the forcing functions is required states computed from (6) are always consistent initial in order not to give rise to distributions in the extended conditions even in higher index cases (e.g. Rovaglio et system resulting from time differentials of u which al., 1992). Though this seems to be correct for semicannot be handled rigorously by integration algorithms. explicit index one DAE systems it is not at all obvious in

Getting around consistent initialization of DAE systems? the general case, because x~ as computed from (6) must also satisfy the extended system (3), (4) at steady state in order to be consistent. Simple reasoning proves this conjecture. Assume a steady state solution (5) has been computed from (6). Now we insert (5) into the consistency equation system (2) to check whether it is also fulfilled. Since x, and u~ are constants, all higher order time derivatives of the state and input vector are zero, and we end up exactly with equation (6) which defines the steady state solution. As the extended system (3), (4) is obtained from the consistency equations (2) by symbolical manipulations it must also be satisfied by the steady state (5). Hence, the steady state (5) forms a set of consistent initial conditions for any DAE system in general. We illustrate our argumentation by means of two simple examples with linear constant coefficient DAE systems of index one and index two, both containing time dependent forcing functions uj and u2. Example 1: The first DAE system ±+~=ax+by+u,

(7)

O=cx+u,

(8)

147

Example 2: For the index two DAE system ~=ax+by+u~

(17)

O=cx+u2

(18)

a steady state solution is readily obtained with both forcing functions constant, cf. equations (12), and setting xs=0 in (17). Solving (17), (18) for x~ and y, we obtain 1

x ~ = - - u2.,

(19)

c

1

a

y , = - ~ u,..,+ bc u,_.~

(20)

~t,=O

(21)

.9.,.=0.

(22)

and

Using the more rigorous approach to compute a steady state solution, we first derive the corresponding extended system (3), (4) by differentiating equation (17) once and equation (18) twice. Solving for the time derivatives and the state variables we obtain

is of index one. The corresponding extended system i

x=-

- u.,,

(23)

C

1

~ t = - - t~_,

(9)

C

1 1

p=ax+by+ul + - it,_

a

1

y = - -b ul + bc u,.- bc u2,

(24)

(10)

C

1

O=cx+u2

is obtained by differentiating equation (8) once and solving for the state derivatives. With both input functions u~, u2 assumed constant Ul=Ul.~, U,.=U2..,,

1

C

l

(25)

C

3' = -

1

a

1

~it,+ ~cU2 - ~cti2.

(26)

(12)

i.e. u~ =0, t~2=O, we obtain the steady state solution x , = - - u,..,

, t = - - u,,

(11)

(13)

a

y . , = - ~ u,..,+ bc u2~

(14)

x,=0

(15)

y~=0

(16)

from the extended system (9)-(11). Following our theoretical argumentation we obtain the same result (13)-(16) if we solve the original DAE (7), (8) for X, and y, with both state variable derivatives set to zero and constant forcing functions u~.s, u~.,.

Since all forcing functions are constants (23)--(26) reduce to the steady state solution (I 9)-(22) as obtained directly from the original equations (17),(18). In summary, since a steady state is always consistent it can be used as an initial value at t = to if sufficiently smooth forcing functions are applied in order to maintain consistency and to guarantee continuous behavior of the state variables in the first integration step. The required smoothness of the forcing functions is discussed in more detail in the next section. DIFFERENTIABILITY OF FORCING FUNCTIONS

If an arbitrary time varying input forcing function u(t)---which is not continuous at time to., i.e. u(t0) U(to÷)---is applied to the DAE for t -> to+, it directly follows from (6) as well as from (2) that the solution at

148

A. KRONERel al.

to is no longer consistent. If we apply continuous but not continuously differentiable forcing functions such as ramps with U(to) = u(t0+) the original system (6) is still satisfied at to+. However, since the extended system (3), (4) also depends on input function derivatives u ~, i = I, .... v, it is not fulfilled any more and the steady state u(to), u~(t0) = 0, i= 1..... v, X(to), :~(to)=0 becomes inconsistent at to+ due to discontinuities in the derivatives of u(t) at

be obtained by structural analysis. Such an analysis which is based on the sparsity patterns (i.e. the nonzero entries) of the Jacobians o f f has been applied previously to analyze the properties of DAE systems (e.g. Pantelides, 1988). Unger et al. (1995) reported an algorithm ALGO based on the Jacobian patterns and structural calculus, which on termination provides the structural index v~, the structural number of dynamics degrees of freedom r , i.e. the number of arbitrary initial values and a structural representation of the underlying ODE. ALGO is applicable to any well-posed DAE system with a structurally regular matrix pencil

to÷: 0g~(x(tt0, u(to), fi(to+). . . . )=~(to+), 0~alt(x(to), U(to), U(to.) . . . . ).

(0, 0,)

To maintain consistency at to+ the forcing functions u~ and their derivatives u~~, i = 1..... p , j = 1 ..... d~, must be continuous at to, i.e. u~'( to+ ) = u~'( to+ ), j =O, 1. . . . . d~; i= 1. . . . . p.

~+a~

.

Instead of analyzing the autonomous system (Unger et al., 1995), the forced equation system (1) is considered in the algorithm ALGOU, an extension of ALGO, as described in Table 1. The inputs to ALGOU are the sparsity patterns of f with respect to x, ~, and u. In addition to ALGO's collection of the underlying ODE in Step 3, ALGOU also collects all underlying algebraic equations appearing in Step 6 in vector II~. The differentiation of the algebraic subset G~(x, IJ~) in Step 7 leads to new higher

Here, the integers d~ refer to the highest order differential of every input function u~ actually occuring in the extended system (3), (4). Symbolical derivation of the extended system (3), (4) would reveal the required order of differentiability of the forcing functions. However, this approach is cumbersome and impractical for large-scale systems. In most cases symbolic derivation is not necessary, since the order of the forcing function differentials occurring can

Table I. Algorithm ALGOU for index analysis, determination of the extended system and the degree of differentiability of input functions Step 0: Step I :

Initialization: i=0,n0=n,r=n,~o=X, Fo=L Y_t=O, U0=u, Y_~=O, II_~=O Determine: ri= rank -oFi - , where i~ e R" and set r---* r - (n~ - r,).

,

Step 2:

Step 3:

FIZ~(x,~',,z,,U,) ]

0 ~i ~0, where I~i",yi • Rr',~i= yir,g T ~ Rn,,~2~,zl • R~,,-r,}

Solve Flll(x,yi,zi,Ui)=O for yi='yi(x,zi,Ui) and append it to Yi_l= Yo,Yl,'-',Yi-I =Yi-i( x, ~ i, Ui-~) to obtain

~[i=[yTO,YTI,...,yiT-t,yS]T=yi(X,~Ii,~'i,Ui) • Step 4: Step 5:

y ,,x,

intoY:,(x, z 1 to o o,n

Ifr~ = ni then i", = ~, v, =iand STOP. / \ / \

Step 6: ~

3

Gi , (x, U~_~)]to obtain O= II~ = [Go(X, Uo), a,(x,U,),... G~.,(x, U~_,), G,~x, U~)I Step 7:

cgGi ~t + flGi 0Gi " Differentiate G,(x,U,) = 0 with respect to time: ~-Yi ~ ~ ~t~+ ~ Ui=0.

Step 8:

Substitute ~'i=Yi(x,~i,Ui) into the result of Step 7 to obtain Fi+l ' (x,zi,Oi,Ui) " " =0. Analyze Fi+ i '(x,zi, Ui,[~Ji) in order to find new higher derivatives of the input variables u not yet contained in U~and append them to U~to give U;+~and to obtain Fi. ,(x,7'i,Ui+ i)-- 0.

Step 9: Step I 0:

xi+l=zi,vli+l-~ni-ri, i ~ i 4- 1. Restart from Step 1.

Getting around consistent initialization of DAE systems? order differentials of the forcing functions not yet contained in vector U~. They are appended to U~ to give Ui+ ~. After termination of the algorithm the function vectors Y,. and II,._, correspond to the functions ~ and of the extended system (3), (4), respectively. The vector U,_ ~ contains all inputs together with all their derivatives ui°~, i = 1... p , j = 0,... di ": v occurring in Yv and II,._~. Hence, all members of U,._~ must be continuous at to+ in order to keep consistency after applying the forcing function. Unfortunately, ALGOU does not reveal the scalar equation in which the highest derivative of a certain input ui is occurring*. This is due to Step 3, where an inversion of the structural matrix ~

must be applied

in the algorithm. Since the structural inverse of a matrix is in general full, all inputs and their derivatives are considered to occur in all equations y~. Note that the algorithm developed by Bujakiewicz and van den Bosch (1994) for the determination of the perturbation index would render such detailed derivative information. To successfully (re)start variable order DAE integration algorithms it is almost always sufficient, that the solution x(t) is at least continuous across to. Hence, (3) need not be fulfilled at t0+, since/~(to+)~(to) does not prevent continuity of the solution x across to. Consistency is only required for (4) in order to prevent a discontinuity in x at to due to a discontinuity in a forcing function or its derivatives. Since the order of the derivatives in (3) is always one higher than the order occuring in (4) (see Table 1) the differentiability requirements on the input function are weakened. The required order of differentiability /x~ for every input function u~ to guarantee at least a continuous solutions for x(t) is given by /z~=max(0, d~- 1),

i=l...p.

(27)

The integer/zj refers to the highest derivative of input u~ which must be continuous. For index one problems, continuity of the solution can be guaranteed and the direct method to start DAE integration is not expected to cause problems if all forcing functions are continuous across to. We revisit Examples 1 and 2 to illustrate the theoretical considerations. Example 1 continued: Comparing the steady state solution (13)--(16) of the index one DAE system with the extended system (9)-(11) at t = to+, it is seen that a step perturbation of the steady state by means of either u or u2 and even a ramp perturbation in uz with u2(to÷)~0 causes the solution to become inconsistent at time to.. As

* This information is not required here but may be useful in other contexts.

149

ut occurs only undifferentiated in the extended system a nonzero first order time derivative in uj for t >- to+ such as in a ramp function does not lead to an inconsistency of the state at to+. Even though a ramp in u2 starting at to, causes inconsistency the solution for x and y remains continuous at to. Thus, automatic stepsize and error control strategies in numerical integration codes applied to index one DAE systems with input ramp functions are expected to allow almost always a successful start of the algorithm. If the structural algorithm of Table 1 is applied to this example instead of the analytical investigations above, it terminates after the first loop with the following result: structural index v~ = 1 structural degree of freedom r~. = 1 structural extended system

, I.g]-17,.(x,y,u,,i~,.)]

i

Ilo=Go=(go.l(x,uz) ) input function derivatives vector U~=[u~,u2,d2] with the highest order derivatives d r = [0,1]. Hence, consistency is guaranteed if u~ and li, are continuous across to whereas continuity of the solution only requires continuity in u, and u 2. The results from the structural analysis are identical with the findings from the analysis of the symbolically derived system. Example 2 continued: It is obvious that after a step change at to in either u~ or u 2 or both, (19)-(22) does not fulfill (23}--(26) anymore and hence, the state x~, 5~s,x , Ys, u,(to÷), u2(to÷) is inconsistent at to+. An analytical analysis of the extended system (23)--(26) reveals that the forcing functions must fulfill the conditions tiz(to÷)=O

(28)

ci~,(to. )+ ii2(to. )=O.

(29)

if the steady state solution at to can be used as a set of consistent initial values after a perturbation at to+.Applying ALGOU to the structural equivalent of equations (17), (18) the following result is obtained after two loops: structural index v~ = 2 structural degree of freedom r, = 0 structural extended system

150

A. KRONERet al.

input function derivatives vectorUT=[u~,u,_,t~,u2,ii2] This is an explicit DAE with the state variables molar and order of derivatives d r = [1, 2]. holdup N, temperature T, pressure p and the liquid If the structural representation of the extended system flowrate L leaving the condenser. The vapor feed flow obtained from ALGOU is compared to the analytical rate F, the feed temperature T~,, the cooling temperature extended system as given by (23)-(26) different depend- T,. as well as the volume V of the condenser are ency patterns are obtained. This is due to the formulation considered as time varying input functions in our of ALGOU where the underlying algebraic equations are numerical experiments which can be prescribed externot solved for state variables which are then substituted nally. in the underlying ODEs as done during the analytical To arrive at the corresponding extended system derivation. The results about the differentiability of the analytically we must differentiate equations (30), (31) input functions are however the same for the structural once and equations (32), (33) twice. Besides the original and the analytical analysis. equations (30)-(33) the extended system comprises three Hence, to guarantee continuous solutions for DAE additional equations system (17)-(18) which is necessary for successfully applying numerical solution algorithms the input function u, must be continuous--e.g, a ramp function--and u2 must be continuously differentiable at least once--e.g. a second order function in time--according to equation (r,.-N-C, + (24). These conditions are more restrictive than it may be +p~Z+ RT(L - F) (34) expected without such a detailed analysis. It is not sufficient to start integration from a consistent steady L=L(N,T,p,L,V,f',I/,F,F,T~,,~,,,ToT"c) (35) state with u~ and u 2 as ramps, i.e. continuous functions, according to

ul.s t<-to u/=] ui..~+/3it to.<-t<-tl l ui.s+/3itt t>-t~+

Bp (FC,(T~,,- T)+LAH+k,(T,.- T)) (36) P= NCp(T+ C) 2

i=1,2.

A discontinuity in the solution ofy due to the step in the first derivative of u2

[ 0 t<--to u 2 = / ~ to.<--t<--t,

t 0 t>--fi+ (see equation (24)) will arise at to+ as well as at t~ where the derivative of the forcing functions u~ and u2 switch from zero to/32 and back again. NUMERICAL ILLUSTRATION

To illustrate our discussion we first present results from numerical experiments computed with the model of a single component condenser according to Pantelides et al. (1988)*: /~= F - L,

(30)

NCeT"=FCt,(T~,, - T)+LAH+k~(T: - TL

(31)

B

0=p - Ae- ~Z?,

(32)

O=pV-NRT.

(33)

* Note that this model is thermodynamically not consistent. Despite this fact it is employed also here due to its widespread use in other papers for illustrative purposes.

For shorter notation equation (35) is not stated explicitly. It is obtained from (34) by differentiation, substitution of available time derivatives of state variables and solving for the unknown L. The condenser problem leads to an extended system of the seven equations (30)-(36). The above DAE system has a differential index of two and one dynamic degree of freedom. Details of this example are provided in Unger et al. (1995). From a structural analysis of equations (30)-(33) by means of ALGOU identical results are obtained. In particular, the input function derivatives vector is UrI=[F,T,.,,,ToV,P, Ti,,,L,~:,~ and the derivatives order vector is d r = [1,1,1,2]. The analysis reveals, that a start of the integration without rigorous initialization using automatic DAE codes is only feasible if the inputs F, T,., and T,, and the input derivative V are at least continuous. In addition, we observe from equations (30), (31), and (36) that the time derivatives of N, T, and p depend on the input functions. So we expect the solution of these state variables to be one order higher differentiable than the input functions F, T~,,, T,. and V themselves. The differential equation (35) for the liquid flowrate L depends on the second time derivative of the volume V and the first time derivatives of all other input functions. Hence, the solution of L is as differentiable as the input functions F, T,,, T,. but also as differentiable as the first derivative of the volume function V. We demonstrate the consequences of these properties

Getting around consistent initialization of DAE systems? during step response computations in two sets of experiments. Assuming all input functions and the volume constant we compute a consistent steady state solution first. This solution serves as starting values at t = 0 in all our experiments. At time to = 5 s we apply a step change or some continuous approximation in order to compute the step response at least approximately. Rigorous initialization using Leirnkuhlefs modified algorithm (Kr6ner et al., 1992) is compared to a direct start of the integration without prior initialization. Since automatic integration codes for index one DAE systems are not applicable in some of the cases studied, we use the implicit Euler method with fixed stepsize of Is and full Newton iteration in the corrector step to compute the simulation results. In the first simulation experiment we change the feed flowrate F(t) from F~ = 50 mol/s at t = to to F 2 = 30 mol/s for t -> to+ after the first five seconds have been simulated. We want to investigate the transient solution to a new steady state. Four alternatives are considered: a: If we introduce the step change in F, do not care about consistency, and start integration immediately, we obtain the solution shown as curve 'a' in Figs 1 and 2. As expected, the solution of the holdup N (Fig. 1) shows an abrupt change in the slope which is indeed one order more continuous than the jump in the feed flowrate F. The solution of T and p behave similarly and are therefore not shown. The flowrate L (line 'a' in Fig. 2) shows a discontinuous jump during the integration step after the step change in F. This step in L causes automatic index one DAE systems integrators like 550

1

151

DASSL, LIMEX and RADAU5 to quit due to a constant contribution (see curve 'b') to the local error in the first step after the step change. b: The same change in F is applied with a subsequent computation of consistent initial values according to Kr6ner et al. (1992). The results for N and L are shown as curves 'b' in Fig. 1 and Fig. 2, respectively. The degree of freedom is chosen as the holdup N which is assumed constant across the discontinuity. The step in the liquid flowrate L now occurs instantaneously due to the consistency computation immediately after the step in F at time t = 5s. If we chose T or p as the degree of freedom and assume it constant across the discontinuity at t = 5s, we obtain the same results. c: Consistent initialization of the DAE system can be bypassed if the discontinuous step change in the forcing function is replaced by the smoothing ramp function

F(t)=Ft t --

F(t)=F~+(F2-FO(--) tl --

to

to

F(t)=F,_

for

t<-to

for

to+<--t<--fi

for

t>tl+

over the time interval [t0,q]. In Fig. 1 and Fig. 2 the ramp from F, to F2 starts at to = 5s and ends at t, = 10s. This method transforms the problem from a single step in F into two steps in its first time derivative occurring at to and fi, respectively. With F now being continuous, we fulfill the requirements for the direct method and obtain from Eqs. (30)-(36) a continuously differentiable solution for the holdup N (curve 'c' in Fig. I) and a continuous solution for the liquid flowrate L (curve 'c' in

I

I

500 Z

450 o o

400

350

00

0

L

i

i

i

2

4

6

8

|

i

10 12 t i m e Is]

i

14

i

16

i

18

20

Fig. I. Initial transients of holdup N for change in feed flowrate F, (a) step change in forcing function without consistent initialization, (b) step change in forcing function with consistent initialization, (c) ramp function in forcing function without consistent initialization at starting and end point, (d) ramp function in forcing function with consistent initialization at starting and end point.

152

A. KRONERet al. f

50

48

I~\

'" , ,

b

c,d X,\

~ 46 o

•~- 44

42 40 0

I

I

I

I

I

I

I

I

I

2

4

6

8

10

12

14

16

18

20

time [s] Fig. 2. Initial transients of flowrate L for change in feed flowrate F, for legend see Fig. 1. Fig. 2) at both points to and tt. However, L jumps whenever F changes slope. An integration algorithm with stepsize and error control is likely to pass both points of discontinuity successfully after an appropriate number of stepsize and order reductions. d: Computing consistent initial values rigorously according to KrGner et al. (1992) at each time the ramp function changes slope, i.e. at t0. and t~+ with N fixed across the discontinuities does not give any difference in the solution as compared to case 'c'. The solution curves of cases 'c' and 'd' are identical for both, N in Fig. 1 and L in Fig. 2. For a change in F a ramp function would be adequate to circumvent the consistency problem numerically when starting from a steady state solution. However, we

obtain a significant difference in the transients to a new steady state. The full transients for N and L are shown in Figs 3 and 4, respectively. The difference between the step responses (cases 'a', 'b') and the ramp responses (cases 'c', 'd') is due to more moles being added in the interval [to,tt] in the ramp case than in the step case resulting in higher values of holdup N. Therefore, the response to a smoothed step should never be interpreted as the step response. Since the differences are largest initially, it may have a significant impact in control system performance studies, such as inventory control in muiticomponent distillation (Rovaglio et al., 1992). In the second simulation experiment we change the volume V of the condenser from Vt = 1 m3 at t = to to V2 = 0.8 m3 at t -> to+. If rigorous initialization ought to

550

~

500

~.

c,d

450 ~"\as ~ -,\

350

250

200 150 i

i

I

i

10

20

30

40

i

i

50 60 time [s]

'

70

80

100

Fig. 3. Full transients of holdup N for change in feed flowrate F, for legend see Fig. 1.

Getting around consistent initialization of DAE systems?

55

i

f

i

i

i

i

i

153

i

50

i ~

45

~o

40

•-" 35 30 5

i

0

i

10

i

20

i

30

40

50 60 time Is]

i

70

i

80

i

90

100

Fig. 4. Full transients of flowrate L for change in feed flowrate F, for legend see Fig. 1. be avoided a continuously differentiable smoothing function must be employed as revealed by the structural analysis. The equation defining the liquid flowrate L depends on the time derivative ~/ of the gas volume. Hence, the solution of L is expected to be continuously differentiable as many times as the first derivative of the volume function, ~'. The same ways as in the first experiment are considered to handle the change: a: Applying the step change in V to the condenser being at a steady state and starting integration immediately without rigorous initialization we get the solution curves 'a' in Fig. 5 and Fig. 6 for the holdup N and the temperature T respectively. The pressure transient is qualitatively similar. Figure 7 shows the evolution of the

540

Q

Z

liquid tlowrate L. The holdup N and temperature T show a step change during the first integration interval. The peak in the solution of L is due to its dependency on the time derivative ~/which in this case is a pulse. The first integration step not only includes the contribution of the system dynamics but also the discontinuous step change as a constant contribution to the difference of V at to and to + h. Therefore, the peak value of the pulse increases with decreasing stepsize. Common error control strategies which assume the error proportional to the change in the state variables during an integration step apply stepsize reductions if the error exceeds the threshold. Here, the error estimate as well as the solution update increase with decreasing

,

a•--.,

c,d

530 520 510

"%....% 500 490

480

0

. 2

.

4

.

. 6

.

8

.

. . . 10 12 time Is]

. . 14 16

18

20

Fig. 5. Initial transients of molar holdup N for change in condenser volume V,for legend see Fig. I.

154

A. KRONERet al.

505

500

/

495

/

/

/

c,d

/ /

/

/

/'

490

485

, 2

0

, 4

I

I

6

8

I

I

I

10 12 time Is]

14

I

I

16

18

20

Fig. 6. Initial transients of temperature T for change in condenser volume V, for legend see Fig. I.

stepsize causing the code to abort with stepsize control failure. The reason for this behavior is the approximation of the pulse function L. Consistently initializing the DAE system rigorously according to Krrner et al. (1992) after the step change in V with the dynamic degree of freedom N kept constant across the discontinuity results in the solution curves 'b' in Figs 5-7. Except N, all state variables show a step disturbance which fades out over time. The peak in the flowrate is eliminated because of the consistent initialization computation procedure. Introducing the change in V by a ramp function without consistent initialization at the initial and final time of the ramp does not lead to a satisfying situation either. This can be expected from the structural inves-

tigations since continuous differentiable forcing functions are required. Even though N and T (curves 'c' in Fig. 5 and Fig. 6, respectively) are continuous the flowrate L in Fig. 7, curve 'c', still shows discontinuous steps at both discontinuities in V. Thus, an automatic integration code is still likely to quit due to constant error contributions in the steps over the discontinuities (see curves 'd' in Figs 5-7). Solutions shown as curves 'd' in Figs 5-7 are obtained with consistent solutions rigorously computed after each discontinuity in V. The only difference to case 'c' is seen in the flowrate L (Fig. 7) where the steps are carried out instantaneously. There are differences among the transients 'a' and 'b', or 'c' and 'd', but they all lead to the same new steady

59 58 57 56 55 o "o °M

54 i" /"i~:'%

53

t

o~,,q

52

I ,/" i / I '"~f i',/ C

51 50 49 0

I

I

I

I

2

4

6

8

I

I

I0 12 time [s]

I

14

I

16

I

18

20

Fig. 7. Initial transients of flowrate L for change in condenser volume V, for legend see Fig. 1.

Getting around consistent initialization of DAE systems? state solution. Due to the pulse in the flowrate L of case 'a' initially a large amount of material is removed from the condenser leading to less holdup (Fig. 8) and a slightly smaller flowrate L (Fig. 9) for the rest of the transition time than in case 'b'. Cases 'c' and 'd' do not differ significantly in the transient after the smoothing interval. Compared to the step responses, a ramp change withdraws fewer holdup resulting in a still larger amount of material left during the transition time. Hence, inaccuracies of the solution are not confined to the first few steps of the simulation interval as might be conjectured. To achieve a solution without rigorous initialization all the state variables must be continuous across the

155

discontinuities. Therefore, according to our structural analysis, the forcing function must be continuously differentiable throughout the entire simulation interval. For a change in the volume V from V, to V2 a suitable smoothing function could be

V(t)=1 V(t)=V,+(V_VI)sin2( t-to -~zc) . tl__t-----~O V(t)= V_,

~~...c,d

520

~.,.a

Z o

500 480

o

460 440 420

I

0

10

I

20

I

30

I

40

I

I

50 60 time [s]

I

70

I

80

I

90

100

Fig. 8. Full transients of holdup N for change in condenser volume V.

60

L,..a

58 56

o

•-

54

b 52 50 I

0

10

I

20

I

30

I

40

t<-to

for

to+<--t<-4

for

t>-t,+

where to is the starting time and t~ the final time of the smoothing function. Figure 10 shows the initial transients of the flowrate L

540

o

for

I

I

50 60 time [s]

I

70

I

80

I

90

Fig. 9. Full transients of flowrate L for change in condenser volume V.

100

A. KRONERet al.

156

59 ,.., 58 rA0

57 56 ~. 55 =0

54

52

#e

51 50 49 0

, 2

, 4

I

I

6

8

I

10

I

I

12

14

I

I

16

18

time [s] Fig. 10. Initial transients of flowrate L for higher order change in condenser volume V,,(a) step change in forcing function without consistent initialization, (b) step change in forcing function with consistent initialization, (e) sin 2 smoothing function without consistent initialization. (curve 'e') for a change in V according to the above function. For comparison the transients caused by a step change in V with and without consistent initialization (eases ' a ' and 'b') are also shown. Now, the flowrate has also a continuous solution. However, transient functions like the one used here often cause the integration algorithms to drastically reduce the integration stepsize and hence, loose efficiency. In our case the stepsize for the Euler method had to be decreased to h = 0.1 sec in the transient phase to obtain correct results.

FUNDAMENTAL LIMITATIONS OF SMOOTHING

The numerical results of the last section show that the type of smooth approximation of the step response as well as the way of initializing the integration have significant impact on the accuracy of the transient solution. One may argue that these differences can be made arbitrarily small at least after the first integration steps if the smoothing time interval is made very small. Implicitly, this argument is based on the assumption that the state x(fi) at the end of the smoothing interval approaches the one obtained from a discontinuous step in the limit tt - to ~ 0. Example 1 continued: If we check this conjecture with Example 1 through integration of equations (9), (I0) in the interval [t0,tj] with smooth u~.,(t), i = 1,2, with ui.,(to) = uM0), and u~.,(fi) = u,(to+) we obtain in the limit t~ to --~ 0 the initial values

y(to÷) = y(to)+ 1 (u2(to.+) _ u2(to))

which are obviously independent of the type of smoothing function and of the length of the smoothing interval. This result is not true for general DAE systems. Brtill and Pallaske (1991) and Briill and Pallaske (1992) study this question for a class of linear implicit index one DAE systems. We review their very simple counter-example subsequently. Example 3: They introduce the following linear implicit DAE of index one: ±+zj~=0

(37)

y - ul =0

(38)

z - u,. = 0.

(39)

Let x(t0) = Y(t0) = Z(to) = ul(to) = u2(to) = 0 and consider unit step changes in Ul and u v The steps are replaced by bounded smoothing functions ui.,, i = 1,2, with ui.,(to) = 0 and uj.,(fi) = 1. For the sake of illustration we choose ul.,(t) = u,(t) and u2.,(t) = [u,.(t)]". From eq. (37)-(39) one obtains after differentiation and manipulation ,t= - uza~.

X(to÷) = - - u,.(to.,), C

(39)

Inserting the smooth forcing functions u~.,, i = 1,2, and integrating over the interval [t0,fi] results in

JL f" d

x(t,)=X(to)- n + l J,o ~ [u,(s)]"÷lds. I

(36)

C

(36) In the limit t, - to "-', 0 one obtains

(39)

Getting around consistent initialization of DAE systems? 1 x(to+)=X(to)- n+---1'

(39)

a value which is obviously depending on the choice of the smoothing functions. These analytical results can be confirmed by numerical simulation. Briiil and Pallaske (1991) prove that the state vector at the end of an arbitrarily short smoothing interval is unique and therefore independent of the type of smoothing only, if there exists an integral form for the left hand side of every differential equation of the DAE system. In Example 1, the integral form is ~ = ax + by + u~ with tp = x + y, whereas in Example 3 no integral form exists for eq. (37). More general conditions are given by these authors in a later paper (BrOil and Pailaske, 1992). DAE systems with no integral form are not at all uncommon in chemical engineering applications. They always arise if the reversible work is not neglected in an energy balance and if at the same time pressure and volume are time-varying input functions or state variables. The following example illustrates this statement. Example 4: Consider for example a closed vessel containing an ideal gas with forced volume u~(t) (by means of a piston) and heat input u2(t) (by means of an electrical heater). The model equations read as U+pf' 0 0 0

= u,, =V-u2, = p V - nRT, = U - nc,,T+ U,,.

Obviously, no integral form exists in this case. However, despite the similarity of these equations to those of Example 3, the difficulties discussed above do not exist here as confirmed by numerical experiments. From the results of BrOil and Pallaske (1991) it is obvious that even for the index one case smoothing of discontinuous forcing functions may be pointless since no unique state may be reached at the end of an arbitrarily short smoothing interval. This complication is a fundamental problem even with rigorous initialization procedures for DAE systems with discontinuities (BrOil and Pallaske, 1991; Briill and Pallaske, 1992; Majer et al., 1995). SUMMARY

In this article we raised the question whether consistent initialization by rigorous methods can be avoided at least in certain cases. In particular, we have been analyzing the frequently occuring problem of studying the transient behavior of a system if a step change in an input is applied to the system in steady state. This case gives rise to the so-called direct method of initialization, where a smoothed step function being continuous at

157

initial time is applied to the system in steady state in order to start integration successfully. We have shown that the steady state as computed from (6) is a consistent initial condition even for high index DAE systems. This fact is a precondition for the application of the direct method of initialization. Further, the forcing functions must display a certain degree of differentiability to prevent failing of automatic integrators. For example, for DAE systems with index one, input functions which are continuous at to+ do not guarantee consistency but at least continuity of the solution. The latter property usually allows error and stepsize control mechanisms to be passed successfully in the initial integration step. In general, an analysis of the extended system is required to investigate the order of differentiability of each single forcing function to maintain continuity of the solution at to. A general algorithm has been given for this purpose. An index two example illustrates convincingly that-even in case of a successful start of the integratiorv-significant errors may be introduced into the solution not only during the initial steps due to nonrigorous initialization. Often, the application of a smoothed step instead of a step in the forcing functions results in a poor approximation of the step response even if rigorous initialization is employed. Finally, it is shown that there are problems where these errors cannot be made arbitrarily small by choosing the smoothing time interval arbitrarily small. In summary, the direct method of initialization must be applied with great care. It may be considered to extend general purpose dynamic simulators by an algorithm to analyze the required and the actual degree ofdifferentiability of the forcing functions. If the applied forcing functions do not display the required degree of differentiability at to or at later time, an automatic smoothing of the forcing functions could be carried out. This action would increase robustness at the price of inaccuracy. Accuracy can only be maintained if robust initialization methods can be provided if discontinuities in the forcing functions or its derivatives occur. REFERENCES

Bachmann, R., L. BrOil, Th. Mrziglod and U. Pallaske, On methods for reducing the index of differential algebraic equations. Comput. Chem. Engng. 14, 1271--1273 (1990). Berzins, M., P. M. Dew and R. M. Furzeland, Developing software for time-dependent problems using the methods of lines and differential-algebraic integrators. Appl. Num. Math. 5, 375---397 (1989). Brenan, K. E., S. L. Campbell and L. R. Petzold. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. North-Holland, New York (1989). BrOil, L. and U. Pallaske. On consistent initialization of differential-algebraic equations with discontinuities. In: Proceedings of the Fourth European Conference on Mathematics in Industry (Hj. Wacker and W. Zulehner, eds.), Teubner, Stuttgart, pp 231-217 (1991).

158

A. KRONERet al.

BrUlI, L. and U. Pallaske, On differential-algebraic equations with discontinuities. 7_. angew. Math. Phys. 43, 319--327 (I 992). Bujakiewicz, P. and P. P. J. van den Bosch. Determination of perturbation index of a DAE with maximum weighted matching algorithm. Preprints IEEE/IFAC Symposium CACSD '94, Tucson, AZ (1994). Deuflhard, P., E. Hairer and J. Zugck, One-step and extrapolation methods for differential-algebraic systems. Numer, Math. 51,501---516 (1987). Gear, C. W. Differential algebraic equations, indices, and integral algebraic equations. SIAM J. Numer. Anal. 27, 1527--1534 (1990). Hairer, E. and G. Wanner. Solving Ordinary Differential Equations II-Stiff and Differential-Algebraic Problems. Springer, Berlin ( 1991). KrOner, A., W. Marquardt and E. D. Gilles. Computing consistent initial conditions for differential-algebraic equations. In: European Symposium on Computer Aided Process Engineering-I (R. Gani, ed.). Proceedings of ESCAPE-I, Elsinore, Denmark, May 24-28, 1992. Supplement to Com. put. Chem. Engng., 16, pp. S131-S138 (1992). Leimkuhler, B., L. R. Petzold and C. W. Gear, Approximation methods for the consistent initialization of differentialalgebraic equations. SIAM J. Numer. Anal. 28, 205--226 (1991). Majer, C., W. Marquardt and E. D. Gilles. Reinitialization of DAEs after discontinuities. Proc. ESCAPE-5, June 1995, Bled, Slovenia. Comput. Chem. Engng., 19, Suppl., $507512 (1995).

Pantelides, C. C. The consistent initialization of differentialalgebraic systems. SIAM J. Sci. Star. Comput. 9, 213--231 (1988). Pantelides, C. C., D. Gritsis, K. R. Morsion and R. W. H. Sargent, The mathematical modelling of transient systems using differential-algebraic equations. Comput. Chem. Engng. 12, 449---454 (1988). Rovaglio, M., T. Faravelli, G. Biardi, P. Gaffuri, and S. Soccol. Precise composition control of heterogeneous azeotropic distillation towers In: European Symposium on Computer Aided Process Engineering-I (R. Gani, ed.). Proceedings of ESCAPE-I, Elsinore, Denmark, May 24-28, 1992. Supplement to Comput. Chem. Engng., 16, pp. SI81-SI88 0992). Simeon, B., C. Filhrer and P. Rentrop, Differential-algebraic equations in vehicle system dynamics. Sun,. Math. Ind. 1, 1--37 (1991). Sincovec, R. F., A. M. Erisman, E. L. Yip and M. A. Epton, Analysis of descriptor systems using numerical algorithms. IEEE Trans. Aut. Control AC,26, 139--147 (1981). Unger, J. and W. Marquardt. Structural analysis of differentialalgebraic equation systems. In: Computer-Aided Process Engineering (L. Piugjaner and A. Espuna, eds.). Proceedings of COPE-91, Barcelona, Spain, October 14-16, 1991. Elsevier Sc. Publ., Amsterdam, pp. 241-246 (1991). Unger, J., A. KrOner and W. Marquardt, Structural analysis of differential-algebraic equation systems-theory and applications. Comput. chem. Engng. 19, 867--882 (1995). Yu, Q. Improving the computational efficiency of step inputs in adsorption processes. Comput. Chem. Engng. 17, 675--677 (199D.