Automatic structural characterization of DAE systems

European Symposium on Computer Aided Process Engineering - 11 R. Gani and S.B. Jorgensen (Editors) 9 2001 Elsevier Science B.V. All rights reserved.

123

Automatic Structural Characterization of DAE Systems E.F.Costa Jr. b, R.C.Vieira b, A.R.Secchi a and E.C.Biscaia Jr. b a Departamento

de Engenharia Quimica- UFRGS - Porto Alegre - RS - Brazil

b Programa de Engenharia Quimica- COPPEAJFRJ - Caixa Postal 68.502 CEP: 21945 - 970 Rio de Janeiro, RJ - Brazil - [email protected] The characterization of the DAE system is unavoidable for the construction of generalpurpose robust codes for numerical integration of DAEs. Such information is needed for index reduction and consistent initialization of DAE systems. In the present contribution, it is proposed the utilization of a graph theoretical approach associated with automatic differentiation tools for the automatic generation and resolution of the consistency system. All differentiation required is performed exactly and at an affordable computational cost, and the index reduction prior to numerical resolution is performed automatically. The resulting code from this project has been incorporated to the numerical integration code DASSLC. Several examples illustrate the potentiality of the approach. Encouraging results were achieved on both traditional benchmarks for initialization procedures and more realistic problems. 1. INTRODUCTION When solving DAEs one must concem about the index of the system and about the consistency of the initial conditions. The index of a DAE is the number of times that all or part of the system must be differentiated with respect to time in order to convert it into an explicit set of first order ODEs. Index 1 systems may be solved with modified ODE methods, while higher index systems (systems with index 2 or greater) require specific methods. Generally, higher index systems are rewritten in an index 1 formulation and solved with standard integration codes like DASSL [ 1] or DASSLC [2], written in FORTRAN and C. During the index reduction, some extra algebraic equations are obtained which generally correspond to derivatives of the original algebraic equations. Those hidden algebraic equations along with the original DAEs compose the extended system. Consistent initial values must satisfy not only the DAE system but also the underlying extended system [3]. Several rigorous methods for initialization have been proposed in the literature, and they all depend to some extent on the characterization of the DAE system. In other words, in order to perform the numerical integration of a DAE with standard integration codes, it is necessary to: (i)perform the characterization of the DAE, in order to determine its index and, if required, the set of equations to be differentiated; (ii)perform the differentiation; (iii)check for feasibility of the initial states arbitrarily set; (iv)solve the consistency system for a consistent initial state. In this contribution it is discussed an extension to the numerical integration code DASSLC. The code developed aims the preprocessing of the DAE in order to reduce its index and/or

124 determine consistent initial conditions. In w the DAE characterization issue is addressed. Existing code and algorithms are discussed and their potentiality and limitations are shown. In w the approach proposed in this contribution is detailed. Some numerical examples presented in w show the potentiality of the proposed approach and explain some of the main features of the code. In w some perspectives of the research in the area and the next developments of the code are discussed. 2. CHARACTERIZATION OF DAE SYSTEMS Consistent initial conditions for DAEs must satisfy not only the original equations in the system but also their differentials in respect to time up to some order. Whether or not this additional requirement actually imposes extra constraints on the initial values depends on the particular problem. Pantelides [4] derived a criterion for determining whether differentiation of a subset of the DAE system provides further constraints to be satisfied by initial conditions. A graphtheoretical algorithm was proposed to locate those subsets which need to be differentiated. Unger et al. [5] implemented two systematically different structural approaches as tools for the analysis of DAEs, in the codes PALG and ALGO. In structural computations, only hard zeros (0) and nonzeros (*) are distinguished, and a DAE system is merely represented by the pattern of its jacobian. Since the computation of the index and of the number of degrees of freedom require (repeated) differentiation, this operation has to be defined in the structural sense. Unger et al. [5] defined the structural differentiation to be of the type linear if the derivative of f(z) depends only on z'. On the contrary, if the derivative of f(z) is a fimction of z and z' the differentiation is said to be of the type highly nonlinear. Generally, neither way of structural differentiation yields the correct pattern of the Jacobian of the differentiated function. The linear and highly nonlinear definitions yield the patterns with minimum and maximum nonzeros entries, respectively. The structural computation as implemented in PALG and ALGO codes present two main drawbacks: (i) as the structural differentiation yields the wrong pattern for the derived equations, it impacts on the feasibility check of an arbitrary initial set; (ii) a purely structural gaussian elimination can not identify a set of equations that is structurally regular but actually singular. A pure numerical approach for the initialization of DAEs with known index v was proposed in [6]. The authors construct the consistency equations with all total differentials of the DAE with respect to time up to order v and extra user specifications. Then a family of forward finite difference approximations of these equations is constructed and the resulting system is solved in a least square sense. The solution of this rank-deficient algebraic problem poses a nontrivial numerical task, and existence and uniqueness of the solution can only be proven for linear constant coefficients DAEs. If only the necessary differentiation was performed in the previous algorithm, the resulting non-linear system would have a full rank Jacobian, and could be in principle solved exactly. Kr6ner et al. [7] proposed an algorithm that consists of three steps: (i) a structural analysis of the problem employing PALG/ALGO codes, in order to determine the index of the DAE, the minimum set of equations to be differentiated and the number of degrees of freedom; (ii) the numerical differentiation of these equations as in Leimkuhler's algorithm; (iii) the solution of the non-linear algebraic system. It must be stressed that steps (i) and (ii) are not exact, and

125 hence the resulting algebraic system may present severe numerical problems. The authors pointed out that the solution of the algebraic system was the most difficult step of the algorithm, and recommended restrictive damping strategies if there are steep initial gradients. Murata and Biscaia Jr. [8] developed the symbolical code INDEX for characterization of DAEs based on symbolical computation and MAPLE software. As output the code identifies general systems of index one or higher and Hessemberg forms of index two and three, besides checking for linearity and solvability of DAEs with certain properties. As an auxiliary task INDEX supplies the patterns of the jacobian matrices to ALGO/PALG structural analysis codes. In spite of the ever known tendency of computer algebra systems to use much memory space and computing time, the authors report success on the characterization of several examples. However, the code can not be coupled with a numerical integration software developed in a standard programming language (as FORTRAN or C), and most of all presents severe limitations on the size of the problems it can handle. 3. THE EXTENDED DASSLC CODE In this contribution, it is discussed the characterization/initialization module recently added to DASSLC. Based on Pantelides' algorithm, it uses the proposed graph-theoretical approach to identify the minimum set of equations to be differentiated in order to generate an index one formulation of the original DAE. Through AD code ADOL-C [9], the differentiation is performed automatically at an affordable computational cost, and, most important of all, the differentiation technique employed does not incur on truncation error. The structural pattern of the jacobian matrices required are computed via numerical perturbation of the original equations and of the exact extended system. The user must supply the DAE ..system F(z,z')=0 according to DASSLC standards. The code will apply Pantelides' algorithm, construct the extended system and inform the number r of degrees of freedom of the system. The user is prompted to give r arbitrary initial conditions, on which the code will perform a structural check for feasibility. It must be stressed that structural feasibility is a necessary but not sufficient condition for feasibility. That comes from the fact that structural computation can not determine if a structurally regular matrix is not in fact a singular matrix, as previously discussed. In the PALG/ALGO codes, structurally feasibility was neither necessary nor sufficient, as the pattern of the derived equations was computed in an approximated sense. In the present code, when an initial arbitrary set is found to be structurally feasible, the code will carry out the resolution of the extended system by a Newton type method. The extended system was constructed via automatic differentiation, and hence is exact. The initial conditions obtained from its resolution can be fed directly to the numerical integrator DASSLC. If there is a failure in any of these stages the user is warned and prompted to supply another arbitrary set. 4. NUMERICAL EXAMPLES 4.1. Condenser The model of a single component condenser presented by Pantelides et al. [10] is frequently adopted in literature for illustrative purposes. Although being a contrived example with non rigorous modeling, the widespread use of this model has transformed it into a kind of "benchmark" for initialization procedures [4,5,11 ].

126

iV = F - L PV=NRT

(la) (lb)

N C p l ~ = FCp (Tin- T)+ LAH + UAs(Tc- T) P = A exp(- B/(T + C))

(lc) (ld)

The index 2 formulation shown in Equations (1) was analyzed. The automatic index reduction was successfully performed, and the number of degrees of freedom (r = 1) was correctly determined. Equations (1 b) and (1 d) have been automatically differentiated in order to uncover the extended system. All variables are feasible arbitrary conditions, and the Newton code was capable to perform the numerical resolution of the extended system. 4.2. The Unavoidable Pendulum

The well known index 3 pendulum problem [4,5,8] is used to illustrate the superiority of the developed code over reported implementations of Pantelides' algorithm. This system presents two degrees of freedom (r=2). =G

(2a)

j, = Vy

(2b)

~x = ~,x

(2c)

~,y = ~ y - 1

(2d)

x 2 + y2 = 1

(2e)

PALG and ALGO codes, due to the inexact differentiation used, can not determine whether the set [Vx, Vy] is feasible or not. Depending on the type of the differentiation used, the codes reached to conflicting solutions. The code reported herein has been capable to correctly determine r and identify all feasible sets. Equations (2a) and (2b) have been differentiated once, and (2e) twice. Additionally, sets leading to singular consistency systems and/or impossible initial values have been eliminated by the resolution of the extended system with the Newton type algorithm implemented. 4.3. Batch Distillation Column

The following example illustrate the utilization of the developed code on a more realistic problem. The start-up of a batch distillation column [12] is typically a high index DAE system with highly nonlinear state equations and equilibrium relationships. Additionally, it can easily become a large scale system as mass and energy balances must be solved for every plate and there are several components involved. In this contribution a column with np plates and nc components is simulated, what leads to system with approximately (np+2).(2 nc+6) equations.

127

13/0 = V/+I Yi+l,i. + Li- ni-l'J

(3a)

Li -ni'j ~ -- V~ Yi,j

N, = ~--,~c ,,,j

(3b)

N i h i = Vi+1 Hi+ ~ + Li_ ~ hi_ 1 - Z i h i -I'~ H i + Qi

(3c)

yi,j = x i j K j ( T ~ ) ,

(3d)

Ni hi

tic = Zj=ll"li,j hj(Ti)

H i = Z~c=lYi,j Hj(T~) ~.c i=l Yi,j = 1 where:

(3e) (3f) i=0, 1..... np+l

j=l,2,...,nc

Hp,i(T)=~:=oe,,kZk~

(3g)

Ky(T):T~3k=oaj,kTk~

ej-Z~clYi,jep,i(rj) In the start-up of the column, the fbllowing relations between liquid and vapor fluxes hold: Vi+ 1 = L i

i: 0,1,..., np

(4)

In previous equations, x and y represents molar fractions on liquid and vapor phases, h and H are enthalpies on liquid and vapor phases, T is temperature, i represents plate and j represents component. All parameter used can be found in [12, pp.211-212]. With the extended DASSLC code discussed herein, it is possible to integrate the high index system directly. All index reductions and manipulation in order to achieve consistent initialization are performed automatically. The user should only give the original DAEs according to DASSLC standards, determine which are the arbitrarily set initial variables and provide initial values for those variables. Characterization results for an example with 5 components and 12 plates are shown in Table 1. It should be emphasized that other sets could be determined arbitrarily. Table 1: Characterization results for example 4.3. nc=5, np=12 Index 2 Equations (3b), (3d) (3e) and (3g) must be differentiated once Degrees of freedom: 70 [ nc (np+2) ] Feasible arbitrary initial set: moles of component i on plate j (nij) 5. CONCLUSIONS

It is presented the extension of the code DASSLC that permits its utilization for the numerical integration of high index problems. The characterization/initialization module added preprocesses the original DAE systems and converts it into an index one problem with consistent initialization. Based on Pantelides' criteria and on the Automatic Differentiation

128 tool ADOL-C, the exact and affordable computation of derivatives required is one of the key features. With the exact calculation of derivatives, structural feasibility became a necessary condition for feasibility, and this criteria is used to perform the first search for consistent arbitrary values. A Newton-type method is then used to eliminate singular but structurally regular subsets and to compute the consistent initial values. The condenser and the pendulum problems have been analyzed. These "benchmarks" have shown that the results achieved with the proposed approach are either equivalent to of better than previously reported characterization results. Additionally, the batch distillation column example has shown that the methodology proposed is specially suited to the highly nonlinear large scale models often faced in the simulation of chemical engineering processes. REFERENCES

1. 2.

3. 4. 5. 6. 7. 8. 9. 10. 11. 12.

Petzold, L. DASSL code, Computing and Mathematics Research Division, Lawrence Livermore National Laboratory, L316, PO Box 808, Livermore, CA 94559 (1989) Secchi., A.R. "DASSLC: User's Manual, a Differential-Algebraic System Solver", Technical Report at Chemical Engineering Department, UFRGS, Porto Alegre, RS, pp.49 (1992), http://www.enq.ufrgs.br/enqlib/numeric.DASSLC Brenan, K., Campbell, S. and Petzold, L. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, New York, Elsevier Sc. Publ. Co (1989) Pantelides, C.C. SIAMJ. Sci. Stat. Comp. 9 (1988), 213-231 Unger, J., Kr6ner, A. and Marquardt, W. Comp. Chem. Eng. 19 (1995), 867-882 B.Leimkuler, L.Petzold and C.Gear, SIAMJ. Num. Anal. 28 (1991), 205-226 Kr6ner, A., Marquardt, W. and Gilles, E. Comp. Chem. Eng. 16 (1992), S 131 Murata, V. and Biscala, E. JR., Comp. Chem. Engng. 21 (1997) $829-$834 Griewank, A., Juedes, D., Mitev, H., Utke, J., Vogel, O. and Walther, A., ACM TOMS 22 (1996) 131-167 - Algor. 755 Pantelides, C.C., Gritsis, D., Morison, K. and Sargent, R.W.S. Comp. Chem. Engng. 12(1988), 449-454 Kr6ner, A., Marquardt, W. and Gilles, E. Comp. Chem. Eng. 21 (1997) 145-158 Holland, C. and Liapis, A.I. Computer Methods for Solving Dynamic Separation Problems, McGraw Hill Publishers, 1983