Structural analysis of the computational properties of QP-DAE systems

European Symposium on Computer-Aided Process Engineering - 14 A. Barbosa-P6voa and H. Matos (Editors) 9 2004 Elsevier B.V. All rights reserved.

1063

Structural Analysis of the Computational Properties of QP-DAE systems Gordon D. Ingram l, Adrien Leitold 2 and Katalin M. Hangos 3. Division of Chemical Engineering, University of Queensland 4072, Australia 2 Department of Mathematics, University of Veszpr6m, POB 158, H-8201, Veszpr6m, Hungary 3 Systems and Control Research Laboratory, Computer and Automation Institute, Hungarian Academy of Sciences, POB 63, H-1518, Budapest, Hungary

Abstract Many lumped dynamic models of process systems can be written without approximation in quasi-polynomial differential-algebraic equation (QP-DAE) form. We propose nonlinearity indices for QP-DAEs based on the system invariants. There are global (whole model) and local versions of the indices. The local ones are based on the implicit subsets (essential L-components) of the model. The indices are assessed against simulations of a vaporiser model with adjustable nonlinearity.

Keywords: computational properties, GLV models, decompositions, nonlinearity

1. Introduction Lumped dynamic process models, which have no spatial distribution of their variables, comprise sets of differential-algebraic equations (DAEs). Their structure can be represented by a 'dynamic representation graph'. In earlier work we used this graph for model analysis, including structural solvability, differential index assessment and model decomposition (Leitold and Hangos, 2001). The quasi-polynomial (QP) formalism can potentially extend the current range of computational properties, such as those above, to measures of nonlinearity and interdependence of the system components (Dfaz-Sierra et al., 1999). The QP formalism has a great advantage: it can capture many kinds of nonlinearity without approximation in a strongly structured universal form (HermindezBermejo and Fair6n, 1995). QP systems have both graph and matrix representations; techniques from both these disciplines can be applied. In this paper, we extend the earlier structural analysis work to assess the nonlinearity of DAE process models in QP form. The QP approach allows the use not only of structural information for model analysis, but also the model's constituent mathematical functions and parameter values.

2. Quasi-polynomial DAEs A set of semi-explicit ordinary differential-algebraic equations is a quasi-polynomial differential-algebraic equation (QP-DAE) system if it can be written in the form:

* Author to whom correspondence should be adressed: [email protected]

1064

jCi

=

/~iXi "l- X i

m A o ;i

E

j =1

"k

-"

+i +

n

n+i)j j =l

7Bj(n+k) "~ '

X i(O)

--

X/0,

i = 1. . . . .

n

(1)

k =1

m

0

Bjk

k =1

d

"*k ~ =I

"k

i=1 ..... d ,

'

m>n+d

k =I

where x and z are vectors of the differential and algebraic variables, respectively. The real matrices AE O, i = 1 ..... n; zi > 0, i = 1..... d. Significantly, many terms in process models naturally have QP form. 2.1 Form invariance and equivalence classes The m quasi-monomials of the system, which are the products

qJ =

",k Bj~ k=l

" kBj(.+~) '

j = 1, " " " ~ m

k=l

that appear in Eq. (1), play an important role in analysing QP-DAEs. Magyar et al. (2003) show that QP-DAEs form equivalence classes jointly with non-minimal QPODEs. The members of a class are related through quasi-monomial transformations: n+d X i

:

HX

Cit' ,

i=1 .....

n+d.

k=l

The vector X contains both the differential and algebraic variables, X is its transformed counterpart, and C~ 9~(n+a)• is an arbitrary invertible matrix. A QP-DAE equivalence class can be characterised by the set of transformation-invariant quantities { m, B/],, BA }.

3. Structural Analysis of DAEs The solvability properties of DAlE models can be investigated with the dynamic representation graph (Leitold and Hangos, 2001). This graph is a sequence of static graphs corresponding to each time step of the numerical method used to integrate the differential equations. The internal structure of each static graph is the same. The graph vertices correspond to the model's variables and given parameters, while directed arcs express their functional dependence, corresponding to the model equations. Directed arcs also join the static graphs at different times to indicate the integration method. Our structural analysis procedure determines the most important solvability properties of a model: structural solvability, differential index and model decomposition. In outline, the steps of the structural analysis method are: 9 Transformation of the model into standard form. 9 Construction of the dynamic representation graph. 9 Type declaration of the vertices according to the model specification. 9 Construction of a reduced graph of the model. 9 Investigation of the solvability properties using the reduced graph.

1065 The reduced graph represents the implicit part of the model. Notably, L-decomposition can be performed on the reduced graph to reveal the structurally solvable subsystems of the model, named the L-components, with a hierarchy between them. The essential Lcomponents are those subsets of the model that cannot be solved by simple substitution.

4. Nonlinearity Measures for QP-DAEs Nonlinearity measures in general may be helpful for classifying a model's behaviour and predicting the difficulty of its solution. Dfaz-Sierra et al. (1999) proposed a graphbased measure for the nonlinearity and connectivity of QP-ODEs. We have taken inspiration from their work to formulate alternative nonlinearity measures for QPDAlEs. First, it seems reasonable to characterise the nonlinearity of a QP-DAE system by global indices, vc, obtained from norms of the invariants of the complete model:

V~ =

{m,IIBtAIII.},

118tAlL --

m i=lmax(~[ .=

(n[A])o]m

)'

(2)

where Ilzll denotes the infinity norm of some matrix Z and [Z] indicates a structural version of Z formed by replacing its non-zero elements with ones. We wish also to explore local versions, VL, of the global indices in Eq. (2) that correspond to the implicit parts of the model. Each essential L-component Li of a QP-DAE system can be characterised by the number of quasi-monomials it contains, mLi, and the local norm

liB,,

[AL/][[~. The matrices ALi and BLi a r e found from the quasi-monomials in the QP form of the equations corresponding to component Li. We assess the nonlinearity measures vs and VLthrough a case study.

5. Case Study: Computational Analysis of a Vaporiser Model Table 1 presents an adaptation of the single component partial vaporiser model of Ponton and Gawthrop (1991). Different nonlinearities can be introduced by selecting the equation of state (case El: ideal gas law, or E2: van der Waals equation) and the enthalpy function (case HI: constant heat capacity, or H2: heat capacity varies linearly with temperature). The variables of the system in Table 1 may be classified as follows:

Differential variables: Mv, ME, Uv, UL Constants~parameters: R, g, M, C~v,flv, ~, ~ , fiE, ~L, a, b, A *, B*, C*, PL, k, U, VT, A, f Algebraic variables: F, V, L, TF, Tv, TL, E, Q, QE, hF, hv, hL, hE, P, P*, Po, Vv, VL, 17v Specification of the set of design variables {F, V, TF, Q, Po } along with the above set of parameters ensures that the system has zero degrees of freedom.

5.1 QP form of vaporiser model The QP-DAE form of the vaporiser model appears in Table 2. In accordance with the auxiliary variables procedure for QP systems (Hermindez-Bermejo and Fair6n, 1995), new variables (y, w and u) and their respective equations were introduced to transform the non-QP nonlinearities in the original model into QP form.

1066 Table 1. Equations of the vaporiser model. M v = E-V lJv = E h E - V h v

+QE

U v = Mvh v E

= kA(e*-

P)

(1.1)

M E = F-E-L

(1.3)

/-)L = F h F - E h E - L I k

(1.5)

U L = MLhL

(1.6)

(1.7)

QE = U A ( T L - T v )

(1.8)

(1.9)

VL = M L / p L

(1.10)

L = f ~L

(1.12)

VV

-- V T - V L

r

= VvM/Mv

(1.11)

p* = 1 0 A'-B'/(TL+C*)

(1.13)

Case E1 P : RTv / Vv

(1.2) +Q-QE

4P+gMLIA-Po

(1.4)

Case E2 (1.14)

Case H 1

P = R T v / ( V v - b ) - a / V v2

(1.14')

Case H2

h v = otv + flvTv

(1.15)

hv = a' v + flvTv + rvTv~

(1.15')

h E = a L + flLTL

(1.16)

hE = aL + flLr~ + rLT/

(1.16')

hE = 6gv + flvrL

(1.17)

hE = o'v + flvrL + YvrL2

(1.17')

h F = OtL + flLTF

(1.18)

hF = a'L + flLTF + rLT~

(1.18')

5.2 Computational analysis L-decomposition of the vaporiser model reveals that there are two essential Lcomponents (Fig. 1) in all four cases considered. The variables Uv*, UL*, Sl and s2 appear when writing the model in standard form (refer to Leitold and Hangos, 2001). The vaporiser model in Table 1 was coded in MATLAB 5.3 and solved using the DAE solver ODE15S. Each combination of options (El/E2 and H1/H2) was simulated for the same initial conditions and period. Table 3 presents the values of the parameters and design variables used in the simulations. Table 4 shows some computational properties measured from the simulations, namely execution times and the number of floating point operations (flops), along with our proposed nonlinearity measures. The execution times and flop counts reflect an intuitive assessment of the difficulty, that is, the nonlinearity, of the problems. Both m and IIB[A]II- provide some indication of computational nonlinearity, but is unclear which is the more useful. The L-components are also a little too small and simple to draw any definitive conclusion about their value for this kind of analysis. Note that only those nonlinearities that appear in the implicit (iterated) part of the model will be reflected in the local indices. Hence they could not capture the case El/E2 differences. We also note that the L-components themselves depend on the algebraic form of the model - they may change if one applies an algebraic (equivalence) transformation to the model (Leitold and Hangos, 1998). Therefore, it is natural to expect that the nonlinearity measures will also change with algebraic transformations. Further investigations are needed on much more complicated models that can be adjusted to display a wider range of nonlinearity.

1067 T a b l e 2. The vaporiser model in Q P - D A E form.

IVI v = M v ( M v ' E - V M v '

(2.1)

)

/~L = M L ( F M L ' - M L ' E - M L

]L)

(2.2)

(Iv = Uv (Uv'Eh~ - v Uv'hv + UV'OE)

(2.3)

0 L = UL(FU;'hF - U ~ E h E - U ~ L h L + QU~ ~ - U ~ Q E )

(2.4)

0 = U v-Mvhv

(2.5)

0 = U L-MLh L

(2.6)

0 = E-kAP*+kAP

(2.7)

0 = QE - U A T L + UATv

(2.8)

(2.9)

0 = VL - p~l ML

o

=

0

=

-V~+Vv+VL

(2.10) (2.11)

Vv - M M v ' V v

0 = E - fzpL P - (f2pLg/A) M L + f2pLPo

(2.12)

0 = P*-y

(2.13)

Case E1

Case E2

0 = P - Rv,,~:' y-

(2.14)

u

0 = PCv - b P - R T v + a V v l - abCv 2 (2.14')

Case H 1

Case H2

-av+hv-~vrv-rvrv

~

0 = - a v + h v - f l v Tv

(2.15)

o =

0 = - a L + h L - flL TL

(2.16)

0 = - - a L + h L - - f l L T L - - y L T2

(2.16')

0 = - a v + h E - f l v TL

(2.17)

o = - ~ v +hE--C~VrL--rVr~

(2.17')

0 = - ( a L + flLTF) + h F

(2.18)

0

(2.18')

All cases:

--(O~'L

=

"[-/~LTF +

YLTF2 )

+

(2.15')

hF

0 = - C * + w - TL

(2.19)

Case H1

= y(B* ln(10)fltr { F M~'hF w-2 - M~'EhE w-2 - M~'LhL w-2 + Q M ~ w -2 - m -L" ~~gEw -2 _ FM;2UL w-2 + M~ZULEw -2 + M~ZULLw -z })

(2.20)

Case H2

= y(B*ln(lO){F M~lhvw-Zu -1 - M~lEhEw-2u -1 - M~'LhLw-Zu -' + a M~'w-2u -~

}

- M~'QEW-2U -' -- F M~2VLw-Zu -' + ML2ULEW-2U -I + ML2ULLw-2u -1 ) (2.20') 0 = - flL + U -- 27/L TL

Tv

(1.15)

-"

or (1.15")

TL

(1.16) or (1.16")

(2.21')

hv

(1.5)

~:

hr. ~:

(1.6)

Uv* I) ~

sl | ~ 9 J

Ut."

s2 ~] I)o J

~

L.1

Figure 1. Essential L-components of the vaporiser model. The numbers in parentheses refer to the corresponding equations in Table 1; < G > is a vertex type declaration.

1068 Table 3. Values of parameters and design variables used in vaporiser simulations. ,OL A* B* C*

790 kg.m-3 9.42448 1312.253 K -32.445 K

Case HI

M 58.08kg.kmol 1 a 1.568x106mns2.kgmol2 b 0.1010 m2s.kgmol1 k 0.000125 kg.s-l.m2.pa1 U 1000W.m-2K-1 Ctv 1.853x105J.kg1 flv CtL -5.930X10s J.kg x fit. Ctv 3.438x10 s J.kg 1 flv O:L -4.687X105J.kg1 fit. ,,,

Case H2

F 3.5 kg.s 1 1.57m3 V 2.5 kg.s 1 0.785m2 5.69x10 4 m2 TF 298.15 K Q 1.547x106W Po 120000Pa 1329 J.kglK l 2171 J.kglK "1 331.1 J.kglK 1 1.570 J.kglK 2 1337 J.kglK 1 ~L 1.388J.kglK -2

VT A f

Table 4. Computational measures from simulations and the nonlinearity indices. Case

Execution

Hops

El, H1 El, H2 E2, H1 E2, H2

time (s) 2.27 2.35 2.33 2.57

3307581 3363403 3390089 3664666

Complete model m 44 47 46 49

IIB[A]II18 20 18 20

Component Lx mL~ 3 4 3 4

IlBLI[AL~]II2 6 2 6

Component mL2 3 4 3 4

IIBL2[AL2]II2 6 2 6

6. Conclusions We proposed two indices to measure the nonlinearity of DAE systems that can be written in quasi-polynomial form. Both global (whole model) and local (L-component) versions of the indices were explored through a case study on a vaporiser. The indices broadly reflected model nonlinearity, but we need to tackle other problems of more widely adjustable nonlinearity to discriminate between them. The QP approach is useful for model analysis because both structural and quantitative measures can be developed.

References Dfaz-Sierra, R., B. Hernhndez-Bermejo and V. Fair6n, 1999, Math. Biosci. 156, 229. Hern~.ndez-Bermejo, B. and V. Fair6n, 1995, Phys. Lett. A 206, 31. Leitold, A. and K.M. Hangos, 1998, In CPM'98:3 rd IEEE European workshop on computerintensive methods in control and data processing, Prague, Czech Republic, 133. Leitold, A. and K.M Hangos, 2001, Computers chem. Engng 25, 1633. Magyar, A., G. Ingram, B. Pongfftcz, G. Szederk6nyi and K.M. Hangos, 2003, On some properties of quasi-polynomial ordinary differential equations and differential-algebraic equations: Report SCL-010/2003, MTA SZTAKI, Budapest (http://daedalus.scl.sztaki.hu). Ponton, J.W. and P.J. Gawthrop, 1991, Computers chem. Engng 15, 803.

Acknowledgments We thank Prof Ian Cameron (University of Queensland) for discussions on the vaporiser model. Hungarian National Research Fund grant T042710 supported part of the research.