- Email: [email protected]
Structural solvability analysis of dynamic process models
Computers and Chemical Engineering 25 (2001) 1633– 1646 www.elsevier.com/locate/compchemeng
Structural solvability analysis of dynamic process models Adrien Leitold a, Katalin M. Hangos b,* a
b
Department of Mathematics, Uni6ersity of Veszpre´m, PO Box 158, H-8201 Veszpre´m, Hungary Systems and Control Research Laboratory, Computer and Automation Institute, Hungarian Academy of Science, PO Box 63, H-1518 Budapest, Hungary Received 6 March 2000; received in revised form 17 July 2001; accepted 17 July 2001
Abstract The variable structure of dynamic process models is represented by a directed graph termed as the representation graph for the purpose of solvability analysis in this paper. Structural solvability analysis, the determination of the structural differential index and the structural decomposition of the differential– algebraic equations (DAE) model set can be performed using the representation graph. The characteristic features of the representation graph for both index 1 and high index semi-explicit DAE models are presented. Based on the above a novel index reduction procedure for high index models is proposed. The notions and methods are illustrated on simple process examples. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Process models; DAE models; Differential index; Solvability; Structural analysis
1. Introduction Lumped process models are in the form of differential – algebraic equations (DAE), which are sometimes difficult to solve numerically due to index and stiffness problems. Therefore the solvability analysis and index reduction of higher index DAE process models are of primary importance for dynamic simulation. The intelligent front-end of dynamic simulators, the so-called computer-aided modeling tools usually provide a visual, interactive means of analyzing the equation –variable patterns in order to choose appropriate variables to satisfy the degrees of freedom or to investigate equation ordering to enhance solution methods. Structural analysis methods can be found in various forms in ModKit (2000), Model.la (Stephanopoulos, Henning, & Leone (1990)) and ICAS (Krogh Jensen, 1998) to mention only a few. It is an important related question how high index process models arise. The effect of the modeling assumptions on the structural solvability has also been investigated in details (Moe, 1995). It is known that the * Corresponding author. E-mail addresses: [email protected] [email protected] (K.M. Hangos).
(A.
Leitold),
change in the specification may transform an index 1 model to a higher index one (Hangos & Cameron, 2001). It is also known that algebraic transformations, that is transformations which only change the algebraic form of a DAE model may change its structural solvability properties, namely its decomposition in a rather drastic way (Leitold & Hangos, 1998). Various methods and techniques have been proposed for structural solvability analysis based on both algebraic (see Unger, Kroner, & Marquardt, 1995; Duff, 1972) and graphical combinatorial (Murota, 1987; Yajima & Tsunekawa, 1982; Sargent & Westerberg, 1964) techniques. Marquardt et al. defined a structural approach to the characterization of DAEs, which is derived from symbolical algorithm (Unger et al., 1995). In addition to this, two systematically different structural approaches were applied to support the development of low index models. Moe et al. also developed two index reduction algorithms, which make the originally high index models manageable by commercially available process simulators (Moe, 1995). The aim of this paper is to propose graph-theoretical methods for index analysis of DAE process models and to develop index reduction techniques for high index models based thereon. The graph representation of DAE process models offers computationally cheap
0098-1354/01/$ - see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 2 7 - X
1634
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
methods for analyzing the differential index and the structural solvability decompositions, which can effectively be used for large systems with several hundred variables (Murota, 1987; Yajima & Tsunekawa, 1982; Moe, 1995). The paper is organized as follows. We start with basic notions on structural solvability analysis and decomposition in the next section. Thereafter, the characteristic features of the representation graph belonging to lumped dynamic process models are described and the analysis and index reduction methods are proposed. Finally conclusions are drawn.
2. Basic notions The basic notions on structural solvability are described here for both algebraic and DAE models. Related questions, such as the representation of DAE models and the handling of model specifications are also discussed. From the viewpoint of the solvability, the most important characteristic feature of process models is the differential index, which determines the difficulty of their solution. Therefore, the basic notions on differential index of DAE models are also discussed in this section. The detailed solvability analysis of the model equations gives information on the computation order and on the size of the nonlinear sub-sets of equations to be solved iteratively. Decompositions of the original DAE model are used for the detailed solvability analysis. The structural solvability of a DAE model and its decompositions can be investigated on the so-called representation graph (Murota, 1987). From the viewpoint of structural solvability two types of model decompositions can be defined and constructed using the representation graph. The goal of these decompositions is to obtain the sub-systems of a model and to determine a successful way for solving the model equations.
2.1. Structural sol6ability and model decomposition As a first step, we consider a system of linear or non-linear algebraic equations in its so-called standard form (Murota, 1987). yi = fi (x, u),
i= 1, …, M
(1)
uk = gk (x, u),
k =1, …, K
(2)
where xj, ( j=1, …, N) and uk, (k =1, …, K) are unknowns, yi, (i=1, …, M) are known parameters, fi, (i= 1, …, M) and gk, (k =1, …, K) are assumed to be sufficiently smooth real-valued functions. We assume that the partial derivatives of functions fi and gk are elements of a field F that is an extension of the rational number field Q. The system of equations above is
structurally sol6able, if the Jacobian matrix J(x, u) as a matrix over F is non-singular, where J(x, u)=
n n
J( f, x) J(g, u)
n n n
J( f, u) , J(g, u)−Ik
(3)
(fi (fi , J( f, u)= (4) (xj (ul (gk (gk J(g, x)= , J(g, u)= (5) (xj (ul and Ik is the k× k identity matrix. It is possible to consider the standard form with no loss of generality in such a way that the variables in left-hand side of equations appear only once (Murota, 1987). It is important to note that the standard form in itself puts no restriction on the form of an algebraic set of equations because it allows both explicit and implicit equations to be present. There are computationally easy algorithms (Murota, 1987) to transform a given set of algebraic equations into standard form. Simple examples of how this is done will be given later on. Consider a system of equations in the standard form. We construct a directed graph to represent the structure of the set of equations in the following way. The vertex-set corresponding to unknowns and parameters is partitioned as X U Y, where X= {x1, …, xN }, U= {u1, …, uK } and Y= {y1, …, yM }. The functional dependence corresponding to the equations is expressed by arcs coming into yi or uk, respectively, from those xj and ul, which appear on the right-hand side. This graph is called the representation graph of the system of equations. It is denoted by G= (V, A, X, Y), where V is the vertex-set, A is the arc-set, X is the entrance and Y is the exit of the graph. The representation graph satisfies the following properties. 1. Each vertex xX has no in-coming arcs (these are initial vertices); 2. each vertex yY has no out-going arcs (these are terminal vertices). The representation graph may be regarded as the flow of information in the system as it expresses the existence of functional dependence among variables. It is important to note that a one-to-one correspondence can be established between the representation graph of a semi-explicit DAE system and its variable– equation structure graph. The occurrence matrix of the representation graph is exactly the variable–equation structure graph of the underlying semi-explicit DAE system when we take into account the assignment between the equations and their left-hand side variables required by the standard form. This way an equivalent algebraic characterization of the representation graph is obtained. In fact a lot of important graph algorithm uses this equivalence and operates on the occurrence matrix of the graph (van Leeuwen, 1990). J( f, x)=
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Eample 1. equations.
Consider
the
following
system
of
y1 = f1(x1, u1, u3)
(6)
y2 = f2(u1, u3)
(7)
y3 = f3(x2, u3, u4)
(8)
u1 = g1(x1, u2)
(9)
u2 = g2(x1, x2, u3)
(10)
u3 = g3(u1)
(11)
u4 = g4(x2, x3, u3)
(12)
The representation graph of the system above is shown in Fig. 1. There is a generality assumption (abbreviated by GA) on the functions fi, (i=1, …, M) and gk, (k=1, …, K) on which the graphical methods of decomposition based which is stated as follows. GA. Let D be the collection of all partial derivatives of fi and gk. The non-vanishing elements of D are algebraically independent over the rational number field Q. This assumption means, that the non-vanishing partial derivatives of fi and gk are so general, that they do not satisfy any polynomial relation with integer coefficients. A Menger-type linking (Murota, 1987) from X to Y is a set of pair-wise vertex-disjoint directed paths from a vertex in X to a vertex in Y. The size of a linking is the number of directed paths from X to Y contained in the linking. In case X = Y , (M= N), a linking of size X is called a complete linking. The graphical condition of the structural solvability is then the following:
Fig. 1. Representation graph of model in Example 1.
1635
Linkage theorem. ( Murota, 1987) Assume GA is true. A system of equations in the standard form is structurally sol6able if and only if there exists on the representation graph a Menger-type complete linking from X to Y. A system of standard form equations is conditionally structurally sol6able if there exists a Menger-type complete linking from X to Y on the representation graph. Here we do not take into account the generality assumption GA. Obviously the conditional structural solvability is a necessary condition of structural solvability. Example 1. (continued) As shown in Fig. 1 there exists a Menger-type complete linking {x1 y1, x2 u2 ul u3 y2, x3 u4 y3} on the representation graph, so the system described in Example 1 is structurally solvable under GA. There are two kinds of decomposition of a graph by Menger-type linking: the so-called L-decomposition and the M-decomposition (Murota, 1987). The L-decomposition is constructed in the following way. Consider a graph G= (V, A, X, Y) with entrance X and exit Y; (XY=¥) such that a Menger-type complete linking exists. Fix a Menger-type complete linking from X to Y. Let be G% a graph, that we construct by adding an arc (y, x) to G for each pair of vertices xX and yY which are linked by the linking. The strong components of G% determine a partition of V and there exists a partial order of that components. The L-decomposition is uniquely determined independently of the choice of the Menger-type complete linking. For a system of equations, the L-decomposition disintegrates the system into structurally solvable subsystems with a hierarchical structure among them. Example 1. (continued) In Fig. 1 the L-decomposition of the system of Example 1 is shown by dashed lines. It can be seen that the representation graph is decomposed into two parts. The first subsystem V1 of this decomposition shows, for example, that we can obtain the unknowns x1, x2, u1, u2 and u3 from parameters y1 and y2 by solving Eqs. (6), (7), (9)–(11). In short the L-decomposition shows which x-unknowns can be obtained from which y parameters. It is important to note that the determination of the strong components in a graph is a polynomial type problem, which can also be done using the occurrence matrix of the graph (van Leeuwen, 1990). Therefore, the L-components can also be found if the equivalent algebraic characterization is given only or if the use of the algebraic description is more convenient. The M-decomposition of a graph leads to a finer decomposition than the L-decomposition. The essence of the decomposition method is the following. There exists a network associated to the representation graph G= (V, A, Y, U) of a system of equations in the standard form. The decomposition of that associated
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1636
network by minimum cuts (so called min-cut decomposition) determines a partition of V and a hierarchical structure of components. Some properties of M-decomposition are the following. 1. The M-components of a system are not necessarily disjoint. A vertex belonging to two components is called a connector. 2. The M-decomposition leads to the finest decomposition, when we do not discriminate between xj and uk. Example 1. (continued) To illustrate above we can see the M-components of the system of Example 1 in Fig. 2.
2.2. Representation of dynamic process models We can adapt the graphical techniques described in Section 2.1 to DAE-system as well (Yajima & Tsunekawa, 1982). An ordinary differential equation of a DAE system can be described by the following equation. x%= f(x1, ..., xn )
(13)
The variable x of the above differential equation is usually determined using a numerical integration method.
&
&
x = x% dt = x%
Fig. 3. Representation of an ordinary differential equation by dynamic graph.
namic graph there are direct arcs attached from the previous static graph to succeeding static graph corresponding to the method for solving ordinary differential equations. In case of an explicit method for solving ordinary differential equations the value of a differential variable at a given time is computed using the corresponding differential value and its value at previous time, for example. Therefore, if there is an ordinary differential equation x% = f(x1, … xn ) in the model, we can rewrite it into the following standard form.
&
x= x%
(15)
x%= f(x1, …, xn )
(16)
(14)
In DAE systems there are two types of variables. Differential 6ariables are the variables with their time derivative present in the model. Variables, which do not have their time derivative present are called algebraic 6ariables. The derivative x% is called deri6ati6e (velocity) 6ariable. A system of equations including also differential equations, can be represented by dynamic graph. A dynamic graph is a sequence of static graphs corresponding to each time step of integration. On a dy-
The representation graph corresponding to Eqs. (15) and (16) is shown in Fig. 3. It can be seen that there are vertices on the static graphs referring to the differential variable x, the derivative variable x% and the variables x1, … xn occurring in Eq. (16). Direct arcs attached to static graphs at time step t and t+ 1 correspond to the applied Euler or other explicit one-step method for solving ordinary differential equations.
Fig. 2. M-decomposition of model in Example 1.
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
2.3. Standard form and the reduced graph The form of the representation graph for DAE system models suitable for structural solvability analysis is developed in steps performed sequentially. The first step for investigation of a dynamic process model by the graph theoretical technique described in Section 2.1 is to rewrite the model into the standard form. The transformation of ordinary differential equations into standard form is given in Section 2.2. The second step is the assignment of types to 6ertices in the representation graph (Iri, Tsunekawa, & Yajima, 1972). The important types of vertices corresponding to the model specification are the following. S (set) -type 6ariablesThese represent variables, which are assigned to the specified given values. These variables require no computation and are present in the specification associated to the process model. The type S is usually assigned to initial vertices of the representation graph.In case of dynamic representation graph we can assume an explicit Euler method (or another one-step explicit method) for solving the differential equations. Therefore, the differential variables will be labeled by type S because their initial value can be obtained from the initial conditions then their values can be calculated step by step by numerical integration. In order to distinguish the label S assigned to differential variables from the label S of other variables with specified given values the label of differential variables is denoted by S*. Labels S and S* are treated the same way in the rest of this paper. G (gi6en) -type 6ariablesA variable assigned to a specific value of a left hand side is a G-type variable. Unlike in case of S-type variables the values of the right hand side variables will be suitably adjusted so as to preserve the equality of the two sides. The type G is usually assigned to terminal vertices of representation graph. A non-terminal vertex assigned to type G can be split into two copies one of which being initial and the other which being terminal vertices. The initial vertex is labeled as type S and the terminal as type G. As a result, vertices of type G form a subset of the terminal vertices. According to the representation graph the value of every variable, which have incoming arcs only from vertices labeled by type S can be calculated by simple substitution into the corresponding equation. These variables become secondarily labeled by type S. After this we can calculate the value of all variables, which have incoming arcs only from vertices primarily and secondarily labeled by type S by simple substitution. These variables will be tertiarily labeled by type S and this process can be repeated if necessary.
1637
Omitting all vertices labeled primarily, secondarily, etc. by type S and all arcs starting from them from the representation graph we obtain the reduced graph. The classification of vertices of a reduced graph is as follows. All initial vertices form the unknown variable set X; all terminal vertices labeled by type G constitute the known variable (parameter) set Y; all other vertices constitute the unknown variable set U.
2.4. Differential algebraic equations and the differential index The degree of difficulty to solve a DAE-system F(z(t), z%(t), t)= 0 can be characterized by the differential index of the model. The definition of the differential index (Brenan, Campbell, & Petzold, 1989) is as follows. The minimum number of times that all or part of a DAE-system F(z(t), z%(t), t)= 0 must be differentiated with respect to time in order to determine z% as a continuous function of z, t is the differential index of the DAE-system. Numerical solution of DAEs includes both initialization and integration. To solve DAEs successfully the initial values must be consistent. In case of ODEs the initial values can be chosen independently hence the so called dynamic degree of freedom is equal to the number of differential equations. In contrast to this the initial values of DAEs can be constrained by the algebraic equations, so the initialization of DAEs can be cumbersome and the problems related to initialization increase with higher index values. Difficulty of numerical integration of DAEs also increase with higher index values (Hairer, Lubich, & Roche, 1989). Dynamic process models can be described by semiexplicit DAEs as follows (Hangos & Cameron, 2001). z%1 = f(z1, z2, t), 0 = g(z1, z2, t)
z1(t0)= z10
(17) (18)
According to the definition of the differential index a semi-explicit DAE has index 1 if and only if one differentiation is sufficient to express z %2 as a continuous function of z1, z2 and t. One differentiation is sufficient only if the Jacobian matrix gz 2 =
n (gi (z2j
is non-singular. Structural non-singularity of gz 2 i.e. the full structural rank of gz 2 results in a model, which is structurally semi-explicit index 1. The initial values of the differential variables in the semi-explicit index one models may be chosen independently.
1638
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
index of dynamic process models. We demonstrate the methods and principles first on an introductory example and then present the general treatment.
3.1.1. An introductory example Let us consider the following example. Example 2. Consider a liquid tank system with one inlet stream and one exit stream as is shown in Fig. 4. Let the vessel be perfectly stirred. Heat is transferred to the liquid using a heater. The flow rate and the enthalpy of inlet stream as well as the flow rate of the outlet stream and the heat transfer rate are functions of time. The model equations of the above system are the following. Fig. 4. A simple liquid system.
If gz 2 is singular then the differential index is greater than 1. Structural singularity of gz 2 is a sufficient condition for differential index being \ 1 numerically. In practice most of dynamic process models having differential index \1 are semi-explicit index 2 models, but there are models having arbitrarily large differential index (Moe, 1995). The initial values of the differential variables in the higher index semi-explicit models cannot be chosen independently. High index models can either be simulated directly by an integrator, which tackles high index DAEs or be transformed to semi-explicit index 1 and then integrated. The principles of the reported index reduction algorithms (Moe, 1995) are the following. Removing algebraic equations from the model, that is to transform sets of DAEs to ODEs; removing differential equations from the model by replacing some of the differential equations of the original model by algebraic equations; increasing the number of equations and variables, that is introducing into the model a new algebraic variable by replacing one of the derivatives by a so called dummy derivative. 3. Investigation of structural solvability of dynamic process models The basic notions and techniques of structural solvability analysis is applied to lumped dynamic process models in this section. The effect of model specification on the solvability is analyzed and an index reduction method is proposed based thereon.
3.1. Relationship between differential index and representation graph of dynamic process models There are several examples in literature, which show the effect of model specification on the differential
M% = −L +F
(19)
U%= − L · hL + F · hF + Q
(20)
hL = f1(TL, p)
(21)
hF = f2(TF, pF)
(22)
U= M · uL
(23)
uL = f3(hL, p)
(24)
L=f4(M)
(25)
Here M denotes the mass; U the internal energy; Q denotes heat transfer rate; F and L are inlet and outlet mass flow rate, respectively; hF and hL are specific enthalpy of inlet and outlet flow, respectively; uL is specific internal energy; TF temperature of inlet flow; pF pressure of inlet flow; TL, temperature in the vessel; p is atmospheric pressure and f1, f2 and f3, are given functions. The standard form of this model is obtained the way as it is described in Sections 2.2 and 2.3, which is as follows.
& &
M= M%
(26)
U= U%
(27)
M%= − L+F
(28)
U%= − L · hL + F · hF + Q
(29)
hL = f1(TL, p)
(30)
hF = f2(TF, pF)
(31)
U M
(32)
u*L = f3(hL, p)
(33)
uL =
s= uL − u*, L
s= 0
(34)
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1639
Fig. 5. Representation graph of model with Specification 1 in Example 2.
L =f4(M) (35) Introduction of new variables u *, L s and the equation s =uL −u *, L s = 0 into the original model guarantees that all variables on the left-hand side of the model are different. In Eq. (34) the ‘satelite’ note s =0 refers to the fact that the value of variable s is set to 0, i.e. the variable vertex s is labeled by type G. The representation graph of the liquid tank system is shown in Fig. 5. We indicate differential variables on the static graph at time step t+1 only, but the structure of this graph is the same as the structure of static graph at time step t. Furthermore, consider the following specification. Specification 1. Given the properties of the inlet flow F, TF, pF and the heat transfer rate Q as functions of time; the initial values of mass and internal energy M0, U0 and the pressure p are constants. To be calculated the mass, the internal energy and temperature of the liquid in the vessel and the outlet flow rate M, U, TL, L as functions of time. It can be shown that the differential index of this semi-explicit DAE system is equal to 1. Dynamic degree of freedom is 2, so initial values of mass M and internal energy U can be chosen independently. The types of variables assigned corresponding to Specification 1 are shown in Fig. 5. The reduced graph according to Specification 1 is shown in Fig. 6. It can be seen that there exists a Menger-type complete linking on the reduced graph hence the model is conditionally structurally solvable.
Next let us consider the following specification. Specification 2. Given F, TF, pF, TL as a function of time; M0, U0, p constants. To be calculated the mass, the internal energy of the liquid in the vessel, the heat transfer rate and the outlet flow rate M, U, Q, L as functions of time. It is easy to show that the model with Specification 2 has semi-explicit index 2, because the unknown variable Q is not present in the algebraic equations, so two differentiations are needed to express Q%. The dynamic degree of freedom is 1, so initial values of mass and internal energy cannot be chosen independently.
Fig. 6. Reduced graph of model with Specification 1 in Example 2.
1640
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Fig. 7. Representation graph of model with Specification 2 in Example 2.
The types of variables, which are assigned corresponding to Specification 2 are shown in Fig. 7. Based on the structure of the representation graph on Fig. 7 the following structural solvability properties can be determined. 1. The standard form model is not structurally solvable (there is no Menger-type complete linking on the graph). 2. There are an o6erspecified part and an underspecified part on the representation graph. The overspecified part indicates the fact that the initial values of the model cannot be chosen independently. The underspecified part Q U% indicates, that Q cannot be calculated from the algebraic equations, so the index of the model is \ 1.
3.1.2. General treatment Inspired by the introductory example the following general statements can be established. The L-decomposition of a reduced static graph results in two types of L-components. Definition. An L-component is essential if it contains at least one vertex labeled by G (the corresponding variable belongs to the known variable set Y). An M-component is essential if it contains at least two vertices. If an L-component is not essential then it contains one vertex only and in this case the corresponding variable can be calculated with simple substitution. Vertices belonging to a Menger-type linking on the reduced graph are always in essential L-components.
Theorem. Consider a dynamic process model M in the standard form according to Sections 2.2 and 2.3 described by a semi-explicit DAE system. Assume that the generality assumption GA is true. Then the differential index of M is equal to 1 if and only if there exists a Menger-type complete linking on the reduced graph at any time step t. Proof. Consider a semi-explicit DAE system in the following form. z%1 = f(z1, z2, t) 0 = g(z1, z2, t)
z1(t0)= z10
(36) (37)
The differential index of this model is equal to 1 if and only if gz 2 is non-singular. The structure of dynamic graphs described in Section 2.2 yields that the derivative variables are missing from the reduced graph or they are not in the essential L-component of reduced graph. Assuming an explicit method for solving ordinary differential equations all differential variables can be labeled by type S* on the static graph at any time step t. Hence the variables belonging to z1 do not appear on the reduced graph. Therefore, all unknown variables of essential L-components belong to z2. So the fact that gz 2 is non-singular is equivalent with the existence of the Menger-type complete linking on the reduced graph. Consequence. Under the assumptions of Theorem above, the following statement holds. If the differential index of M is greater than 1 then there is no Menger-
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1641
type complete linking on the static graph at any time step t.
hF = f2(TF, pF)
(45)
U= M · uL
(46)
The properties of a static graph of a dynamic model, which has differential index \ 1 are as follows. 1. The fact that the initial values of integral variables cannot be chosen independently results in an overspecified part on the graph. This situation can be easy shown by assignment of types to vertices corresponding to the model specification. There is an overspecified part on the graph if a vertex labeled by type G can also be labeled secondarily or tertiarily or etc. by type S. 2. Non-singularity of gz 2 results in an underspecified part on the graph. In this part those algebraic variables appear, which cannot be calculated from algebraic equations and those differential variables, which we want to calculate from them.
uL = f3(hL, p)
(47)
L=f4(M)
(48)
3.2. Index reduction on the representation graph There is a natural method for index reduction suggested by the structure of the representation graph of semi-explicit DAE-systems. Example2. (continued) Consider the Example 2 with Specification 2. According to Section 3.1 this model has differential index 2. Modify the standard form of the model as follows. Instead of equation U = U% introduce the following equation. 1 U%= (U t + h − U t) h where U t and U t + h denote the value of internal energy U at time step t and t + h, respectively, and h is the step length. This equation corresponds to numerical deri6ation so if we get the derivative U% by numerical derivation we can calculate Q from it. In order to simplify notation, let us scale time locally such that h=1 holds. Then the new equation replacing U = U% is the following in the scaled time. U% = U t + 1 − U t
(38)
According to the above transformation of the standard form let the modified standard form of model of Example 2 be the following.
&
M = M%
(39)
M%= − L+ F
(40)
U%= − L · hL +F · hF +Q U%*=U
t+1
−U
s = U% − U%*, hL = f1(TL, p)
t
s =0
(41) (42) (43) (44)
During the modification the internal energy U turns from a differential variable into an algebraic variable and U% turns from a derivative variable into an algebraic one. An important consequence of the modification is that there is no integration for variable U in the modified model, so it is not necessary to gi6e initial 6alue for the 6ariable U. The representation graph of the modified model can be seen in Fig. 8. The static graph at time step t+ 1 shows vertices belonging to variables only, which are necessary to calculate Q in time step t. The reduced graph and the L- and M-components of this part of calculation is shown in Fig. 9. It can be seen that there exists Menger-type complete linking on the reduced graph, so the modified model is conditionally structurally solvable. Inspired by the above example we generalize the idea for a wide class of semi-explicit DAE systems and formulate a general index reduction algorithm. Consider a dynamic representation graph belonging to the standard form of a model M described by a semi-explicit DAE-system. Assume that an underspecified and overspecified part can be found on its static graphs belonging to any time step t. Moreover, the overspecified part contains a differential variable while the underspecified part contains a derivative variable. Assume that the underspecified part contains only one input vertex, i.e. it is under specified by once (this is the most usual form of the models with differential index larger than 1).
3.2.1. Index reduction algorithm 1. Let us form the following variable sets. I0 is the set of the differential variables belonging to the overspecified subgraph, D0 is the set of the derivative variables referring to the differential variables of set I0 I1 is the set of differential variables from which directed paths lead to the derivative variables in the set D0, D1 is the set of derivative variables referring to the differential variables of set I1,… Ik is the set of differential variables from which directed paths lead to the derivative variables in the set Dk − 1, Dk is the set of derivative variables referring to the differential variables of set Ik …
1642
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Fig. 8. Representation graph of modified model in Example 2.
2. Let n be the smallest natural number for which the set Dn contains some deri6ati6e 6ariables of the underspecified subgraph. Then the differential index of the model is wd =n+ 2
(49)
If there is no such number n then the model is not structurally solvable. 3. Based on the sets Ik, Dk (k =0, …, n) let us form a finite series of differential and derivative variables as follows: x0, x %0, x1, x%1 …, xn, x %n, where xk Ik (k=0, …, n), there is a directed path from xk to x − %k l (k =1, …, n). 4. Using the series of variables constructed in item 3 let us modify the model as in Example 2 with Specification 2 before. The modification sub-steps can be seen in Fig. 10. This modification essentially means one index reduction step of the model. As a result of transformations the differential index of the model decreases by one in every iteration step and there is an index 1 model obtained at the end of this iterative index reduction process.
Remarks. There are some important remarks about the above general index reduction algorithm, which are as follows.
1. In case of index 2 models the algorithm presented in Fig. 10 includes only one iteration step. Then the above mentioned differential variable x0 can be determined using the simple condition as follows. let x0 be on the overspecified subgraph, and let x %0 be on the underspecified subgraph. 2. If we do not take into account the fulfillment of the generality assumption GA defined in the Linkage theorem then only the structural differential index and the conditional structural solvability of the model modified by the index reduction algorithm can be investigated. 3. The steps in the algorithm can all be performed using the equivalent algebraic occurrence matrix description of the representation graph augmented
Fig. 9. Reduced graph of modified model in Example 2.
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1643
dT K = K1(Tin − T)− K2 · K3 · C exp − 4 dt T − K5(T−Tc)
Fig. 10. The index reduction algorithm.
with the vertex labels appearing as the symbolic values of the entries in the occurrence matrix. Moreover, every step can be executed in polynomial time including the partitioning of the graph into L-components, therefore, the index reduction algorithm is a polynomial time algorithm. These properties enable to use the algorithm for large problems with several hundred variables, as well. 4. The idea of using the ‘reverse information flow’ for computations, which is the essence of the index reduction algorithm, has appeared in a pioneering work by Smith (1985) but in another context.
(50)
K4 dC = K1(Cin − C)− K3 · C exp − dt T
(51)
Cin = Cin(t)
(52)
Tin = Tin(t)
(53)
Tc = Tc(t)
(54)
where C is the reactant concentration in the reactor; Cin is the reactant concentration in the feed; T stands for the temperature in the reactor; Tin the temperature in the feed; Tc the temperature of cooling water and Ki (i= 1, …, 5) are constants. Specification 1. Given Cin, Tin and Tc, as functions of time; K1, K2, K3, K4, K5 and C0, T0 (initial values) constants. To be calculated C, T, as functions of time. The differential index of DAE system above with Specification 1 is equal to 1 while the substitution of Eqs. (52)–(54) into the differential equations results in an ODE system of which the differential index is equal to 0. Assume now that the modeling goal is changed and the modeler now wants to calculate the values of cooling water Tc(t) which makes the reactant concentration C satisfy the specification C(t)= b(t). The specification for this case is as follows. Specification 2. Given Cin, Tin, and C, as functions of time; K1, K2, K3, K4, K5 and T0 (initial values) constants. To be calculated T, Tc, as functions of time. It can be shown that the differential index of this system with Specification 2 is equal to 3 (Moe, 1995). Now let us apply the method proposed in Section 3.2 for the determination of the structural differential index and for the index reduction. First we need the standard form of the model, which is as follows.
& &
3.3. A further example
T = T%
(55)
In this section there is a further example to demonstrate the proposed index reduction algorithm on a simple high index dynamic process model.
C= C%
(56)
Example 3. Consider a well-stirred reactor equipped with a cooling jacket. A first order exothermic chemical reaction takes place in the reactor. We assume that the volume of the liquid phase mixture is constant in the reactor and the physico-chemical properties are constants. The system can then be described by the following model equations (Moe, 1995).
(57)
T%= K1(Tin − T)− K2 · K3 · C exp −
K4 − K5(T−Tc) T
C%= K1(Cin − C)− K3 · C exp −
K4 T
(58)
The dynamic representation graph belonging to the standard form can be seen in Fig. 11. The assignment of labels to variables is selected according to Specification 2.
1644
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
It can be shown that there is no Menger type complete linking on the representation graph, i.e. the standard form is not structurally solvable. There is an underspecified and an overspecified subgraph on the static graphs. As it is seen from the standard form equations, the reason of the overspecification is that the variable C has to be determined by integration (the label S* refers to that) but the value of C is given at every time step from elsewhere (the label S refers to that). The variable sets defined in step 1 of the index reduction algorithm are determined from the representation graph as follows. I0 = {C}
(59)
D0 ={C%}
(60)
I1 = {T, C}
(61)
D1 ={T%, C%}
(62)
Since the variable T% is present in the underspecified subgraph, n= 1 and wd =n + 2 = 3. The variable series, which can be used for the model transformation is C, C%, T, T%. Therefore, the index reduction can be done in two iteration steps. As first step we replace the equation C = C% for variable C with the equation C%= C t + 1 −C t
(63)
in the model and cancel the label S* of C. The representation graph of the resulted model after the first iteration step can be seen in Fig. 12. In the second iteration step we replace the equation T= T% for variable T with the equation
T% = T t + 1 − T t
(64)
and cancel the label S* of T. The representation graph of the modified model can be seen in Fig. 13 while its reduced graph is exhibited in Fig. 14. There is now a Menger-type complete linking on the reduced graph, i.e. after the two iteration steps the model is conditionally structurally solvable. The following two facts indicate that the differential index of the original model is equal to 3. The first one is that two equations of x = x% type have to be dropped from the original model. The other is that for the calculation of the values of all variables of the modified model at time step t the values of some variables at time steps t+1 and t+ 2 have to be used, as it can be seen in Fig. 13.
4. Conclusion In this paper a graph-theoretical graphical technique is presented, which is suitable to compute the structural differential index of lumped dynamic process models as well as to determine the structural decomposition of the model equations. The characteristic features of the representation graph of index 1 and high index semi-explicit DAE models are also described. Based on the structural solvability properties a novel index reduction algorithm is proposed for both index 2 and higher index models. The developed algorithm can also be used for investigating the structural differential index and the structural solvability of these models.
Fig. 11. The representation graph of the standard form model in Example 3.
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1645
Fig. 12. The representation graph of Example 3 after the first iteration step
Fig. 13. The representation graph of Example 3 after the second iteration step.
It is important to note that the applied model transformation method for index reduction, i.e. to re-formulate ODEs with ‘reverse flow’ of information as difference equations, cannot only be used in the framework of structural analysis using a representation graph but also in connection with other analysis methods based on algebraic matrix representation, too.
Acknowledgements This research was supported by the Hungarian National Research Fund through grant no. T-026575, which is gratefully acknowledged. We also thank the helpful comments and criticism of one of the anonymous reviewers, which inspired us to improve the paper substantially.
1646
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Fig. 14. The reduced graph of Example 3 after the second iteration step.
References Brenan, K. E., Campbell, S. L., & Petzold, L. R. (1989). Numerical solution of initial 6alue problems in differential –algebraic equations. New York: North-Holland. Duff, I. (1972). Analysis of sparse systems. D.Phil, Oxford University. Hairer, E., Lubich, C., & Roche, M. (1989). The numerical solution of differential– algebraic systems by Runge –Kutta methods. Berlin: Springer.
Hangos, K. M., & Cameron, I. T. (2001). Process modelling and model analysis. New York: Academic Press. Iri, M., Tsunekawa, J., Yajima, K. (1972). The Graphical Techniques Used for a Chemical Process Simulator ‘JUSE GIFS’, Information Processing 71 (Proc. IFIP Congr. 71), vol. 2. (Appl.), 1150 –1155. Krogh Jensen, A. (1998). Generation of Problem Specific Simulation Models Within an Integrated Computer Aided System. Ph.D. thesis, Danish Technical University. van Leeuwen, J. (1990). Handbook of theoretical computer science, Algorithms and complexity, vol. A. Amsterdam: Elsevier-MIT Press. Leitold, A., Hangos, K.M. (1998). On algebraic model transformations. Proceedings of third IEEE workshop on computer-intensi6e methods in control and data processing (pp. 133 – 138), CMP’98 Prague, Czech Republic. ModKit (2000). Computer aided process modeling (ModKit) http:// www.lfpt.rwth-aachen.de/research/modeling/modkit.html. Moe, H. I. (1995). Dynamic process simulation, Studies on modeling and index reduction. Ph.D. Thesis, University of Trondheim. Murota, K. (1987). Systems analysis by graphs and matroids. Berlin: Springer. Sargent, R. W. H., & Westerberg, A. W. (1964). Speedup in chemical engineering design. Transactions of Institute of Chemical Engineers, 42, 190 – 197. Smith, G. J. (1985). Dynamic simulation of chemical processes. Ph.D. Thesis, Cambridge University. Stephanopoulos, G., Henning, G., & Leone, H. (1990). MODEL.LA a modeling language for process engineering — 1 the formal framework. Computers and Chemical Engineering, 8, 813–846. Unger, J., Kroner, A., & Marquardt, W. (1995). Structural analysis of differential-algebraic equation systems —theory and application. Computers and Chemical Engineering, 19 (8), 867 – 882. Yajima, K., & Tsunekawa, J. (1982). Graphical techniques used for a dynamic chemical process simulation. System modeling and optimization proceedings 10th IFIP conference New York City, 1981. In R. F. Denrick, & F. Kozin, Lecture notes in control and information sciences, vol. 38 (pp. 826 – 833). Berlin: Springer.
Structural solvability analysis of dynamic process models Adrien Leitold a, Katalin M. Hangos b,* a
b
Department of Mathematics, Uni6ersity of Veszpre´m, PO Box 158, H-8201 Veszpre´m, Hungary Systems and Control Research Laboratory, Computer and Automation Institute, Hungarian Academy of Science, PO Box 63, H-1518 Budapest, Hungary Received 6 March 2000; received in revised form 17 July 2001; accepted 17 July 2001
Abstract The variable structure of dynamic process models is represented by a directed graph termed as the representation graph for the purpose of solvability analysis in this paper. Structural solvability analysis, the determination of the structural differential index and the structural decomposition of the differential– algebraic equations (DAE) model set can be performed using the representation graph. The characteristic features of the representation graph for both index 1 and high index semi-explicit DAE models are presented. Based on the above a novel index reduction procedure for high index models is proposed. The notions and methods are illustrated on simple process examples. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Process models; DAE models; Differential index; Solvability; Structural analysis
1. Introduction Lumped process models are in the form of differential – algebraic equations (DAE), which are sometimes difficult to solve numerically due to index and stiffness problems. Therefore the solvability analysis and index reduction of higher index DAE process models are of primary importance for dynamic simulation. The intelligent front-end of dynamic simulators, the so-called computer-aided modeling tools usually provide a visual, interactive means of analyzing the equation –variable patterns in order to choose appropriate variables to satisfy the degrees of freedom or to investigate equation ordering to enhance solution methods. Structural analysis methods can be found in various forms in ModKit (2000), Model.la (Stephanopoulos, Henning, & Leone (1990)) and ICAS (Krogh Jensen, 1998) to mention only a few. It is an important related question how high index process models arise. The effect of the modeling assumptions on the structural solvability has also been investigated in details (Moe, 1995). It is known that the * Corresponding author. E-mail addresses: [email protected] [email protected] (K.M. Hangos).
(A.
Leitold),
change in the specification may transform an index 1 model to a higher index one (Hangos & Cameron, 2001). It is also known that algebraic transformations, that is transformations which only change the algebraic form of a DAE model may change its structural solvability properties, namely its decomposition in a rather drastic way (Leitold & Hangos, 1998). Various methods and techniques have been proposed for structural solvability analysis based on both algebraic (see Unger, Kroner, & Marquardt, 1995; Duff, 1972) and graphical combinatorial (Murota, 1987; Yajima & Tsunekawa, 1982; Sargent & Westerberg, 1964) techniques. Marquardt et al. defined a structural approach to the characterization of DAEs, which is derived from symbolical algorithm (Unger et al., 1995). In addition to this, two systematically different structural approaches were applied to support the development of low index models. Moe et al. also developed two index reduction algorithms, which make the originally high index models manageable by commercially available process simulators (Moe, 1995). The aim of this paper is to propose graph-theoretical methods for index analysis of DAE process models and to develop index reduction techniques for high index models based thereon. The graph representation of DAE process models offers computationally cheap
0098-1354/01/$ - see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 2 7 - X
1634
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
methods for analyzing the differential index and the structural solvability decompositions, which can effectively be used for large systems with several hundred variables (Murota, 1987; Yajima & Tsunekawa, 1982; Moe, 1995). The paper is organized as follows. We start with basic notions on structural solvability analysis and decomposition in the next section. Thereafter, the characteristic features of the representation graph belonging to lumped dynamic process models are described and the analysis and index reduction methods are proposed. Finally conclusions are drawn.
2. Basic notions The basic notions on structural solvability are described here for both algebraic and DAE models. Related questions, such as the representation of DAE models and the handling of model specifications are also discussed. From the viewpoint of the solvability, the most important characteristic feature of process models is the differential index, which determines the difficulty of their solution. Therefore, the basic notions on differential index of DAE models are also discussed in this section. The detailed solvability analysis of the model equations gives information on the computation order and on the size of the nonlinear sub-sets of equations to be solved iteratively. Decompositions of the original DAE model are used for the detailed solvability analysis. The structural solvability of a DAE model and its decompositions can be investigated on the so-called representation graph (Murota, 1987). From the viewpoint of structural solvability two types of model decompositions can be defined and constructed using the representation graph. The goal of these decompositions is to obtain the sub-systems of a model and to determine a successful way for solving the model equations.
2.1. Structural sol6ability and model decomposition As a first step, we consider a system of linear or non-linear algebraic equations in its so-called standard form (Murota, 1987). yi = fi (x, u),
i= 1, …, M
(1)
uk = gk (x, u),
k =1, …, K
(2)
where xj, ( j=1, …, N) and uk, (k =1, …, K) are unknowns, yi, (i=1, …, M) are known parameters, fi, (i= 1, …, M) and gk, (k =1, …, K) are assumed to be sufficiently smooth real-valued functions. We assume that the partial derivatives of functions fi and gk are elements of a field F that is an extension of the rational number field Q. The system of equations above is
structurally sol6able, if the Jacobian matrix J(x, u) as a matrix over F is non-singular, where J(x, u)=
n n
J( f, x) J(g, u)
n n n
J( f, u) , J(g, u)−Ik
(3)
(fi (fi , J( f, u)= (4) (xj (ul (gk (gk J(g, x)= , J(g, u)= (5) (xj (ul and Ik is the k× k identity matrix. It is possible to consider the standard form with no loss of generality in such a way that the variables in left-hand side of equations appear only once (Murota, 1987). It is important to note that the standard form in itself puts no restriction on the form of an algebraic set of equations because it allows both explicit and implicit equations to be present. There are computationally easy algorithms (Murota, 1987) to transform a given set of algebraic equations into standard form. Simple examples of how this is done will be given later on. Consider a system of equations in the standard form. We construct a directed graph to represent the structure of the set of equations in the following way. The vertex-set corresponding to unknowns and parameters is partitioned as X U Y, where X= {x1, …, xN }, U= {u1, …, uK } and Y= {y1, …, yM }. The functional dependence corresponding to the equations is expressed by arcs coming into yi or uk, respectively, from those xj and ul, which appear on the right-hand side. This graph is called the representation graph of the system of equations. It is denoted by G= (V, A, X, Y), where V is the vertex-set, A is the arc-set, X is the entrance and Y is the exit of the graph. The representation graph satisfies the following properties. 1. Each vertex xX has no in-coming arcs (these are initial vertices); 2. each vertex yY has no out-going arcs (these are terminal vertices). The representation graph may be regarded as the flow of information in the system as it expresses the existence of functional dependence among variables. It is important to note that a one-to-one correspondence can be established between the representation graph of a semi-explicit DAE system and its variable– equation structure graph. The occurrence matrix of the representation graph is exactly the variable–equation structure graph of the underlying semi-explicit DAE system when we take into account the assignment between the equations and their left-hand side variables required by the standard form. This way an equivalent algebraic characterization of the representation graph is obtained. In fact a lot of important graph algorithm uses this equivalence and operates on the occurrence matrix of the graph (van Leeuwen, 1990). J( f, x)=
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Eample 1. equations.
Consider
the
following
system
of
y1 = f1(x1, u1, u3)
(6)
y2 = f2(u1, u3)
(7)
y3 = f3(x2, u3, u4)
(8)
u1 = g1(x1, u2)
(9)
u2 = g2(x1, x2, u3)
(10)
u3 = g3(u1)
(11)
u4 = g4(x2, x3, u3)
(12)
The representation graph of the system above is shown in Fig. 1. There is a generality assumption (abbreviated by GA) on the functions fi, (i=1, …, M) and gk, (k=1, …, K) on which the graphical methods of decomposition based which is stated as follows. GA. Let D be the collection of all partial derivatives of fi and gk. The non-vanishing elements of D are algebraically independent over the rational number field Q. This assumption means, that the non-vanishing partial derivatives of fi and gk are so general, that they do not satisfy any polynomial relation with integer coefficients. A Menger-type linking (Murota, 1987) from X to Y is a set of pair-wise vertex-disjoint directed paths from a vertex in X to a vertex in Y. The size of a linking is the number of directed paths from X to Y contained in the linking. In case X = Y , (M= N), a linking of size X is called a complete linking. The graphical condition of the structural solvability is then the following:
Fig. 1. Representation graph of model in Example 1.
1635
Linkage theorem. ( Murota, 1987) Assume GA is true. A system of equations in the standard form is structurally sol6able if and only if there exists on the representation graph a Menger-type complete linking from X to Y. A system of standard form equations is conditionally structurally sol6able if there exists a Menger-type complete linking from X to Y on the representation graph. Here we do not take into account the generality assumption GA. Obviously the conditional structural solvability is a necessary condition of structural solvability. Example 1. (continued) As shown in Fig. 1 there exists a Menger-type complete linking {x1 y1, x2 u2 ul u3 y2, x3 u4 y3} on the representation graph, so the system described in Example 1 is structurally solvable under GA. There are two kinds of decomposition of a graph by Menger-type linking: the so-called L-decomposition and the M-decomposition (Murota, 1987). The L-decomposition is constructed in the following way. Consider a graph G= (V, A, X, Y) with entrance X and exit Y; (XY=¥) such that a Menger-type complete linking exists. Fix a Menger-type complete linking from X to Y. Let be G% a graph, that we construct by adding an arc (y, x) to G for each pair of vertices xX and yY which are linked by the linking. The strong components of G% determine a partition of V and there exists a partial order of that components. The L-decomposition is uniquely determined independently of the choice of the Menger-type complete linking. For a system of equations, the L-decomposition disintegrates the system into structurally solvable subsystems with a hierarchical structure among them. Example 1. (continued) In Fig. 1 the L-decomposition of the system of Example 1 is shown by dashed lines. It can be seen that the representation graph is decomposed into two parts. The first subsystem V1 of this decomposition shows, for example, that we can obtain the unknowns x1, x2, u1, u2 and u3 from parameters y1 and y2 by solving Eqs. (6), (7), (9)–(11). In short the L-decomposition shows which x-unknowns can be obtained from which y parameters. It is important to note that the determination of the strong components in a graph is a polynomial type problem, which can also be done using the occurrence matrix of the graph (van Leeuwen, 1990). Therefore, the L-components can also be found if the equivalent algebraic characterization is given only or if the use of the algebraic description is more convenient. The M-decomposition of a graph leads to a finer decomposition than the L-decomposition. The essence of the decomposition method is the following. There exists a network associated to the representation graph G= (V, A, Y, U) of a system of equations in the standard form. The decomposition of that associated
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1636
network by minimum cuts (so called min-cut decomposition) determines a partition of V and a hierarchical structure of components. Some properties of M-decomposition are the following. 1. The M-components of a system are not necessarily disjoint. A vertex belonging to two components is called a connector. 2. The M-decomposition leads to the finest decomposition, when we do not discriminate between xj and uk. Example 1. (continued) To illustrate above we can see the M-components of the system of Example 1 in Fig. 2.
2.2. Representation of dynamic process models We can adapt the graphical techniques described in Section 2.1 to DAE-system as well (Yajima & Tsunekawa, 1982). An ordinary differential equation of a DAE system can be described by the following equation. x%= f(x1, ..., xn )
(13)
The variable x of the above differential equation is usually determined using a numerical integration method.
&
&
x = x% dt = x%
Fig. 3. Representation of an ordinary differential equation by dynamic graph.
namic graph there are direct arcs attached from the previous static graph to succeeding static graph corresponding to the method for solving ordinary differential equations. In case of an explicit method for solving ordinary differential equations the value of a differential variable at a given time is computed using the corresponding differential value and its value at previous time, for example. Therefore, if there is an ordinary differential equation x% = f(x1, … xn ) in the model, we can rewrite it into the following standard form.
&
x= x%
(15)
x%= f(x1, …, xn )
(16)
(14)
In DAE systems there are two types of variables. Differential 6ariables are the variables with their time derivative present in the model. Variables, which do not have their time derivative present are called algebraic 6ariables. The derivative x% is called deri6ati6e (velocity) 6ariable. A system of equations including also differential equations, can be represented by dynamic graph. A dynamic graph is a sequence of static graphs corresponding to each time step of integration. On a dy-
The representation graph corresponding to Eqs. (15) and (16) is shown in Fig. 3. It can be seen that there are vertices on the static graphs referring to the differential variable x, the derivative variable x% and the variables x1, … xn occurring in Eq. (16). Direct arcs attached to static graphs at time step t and t+ 1 correspond to the applied Euler or other explicit one-step method for solving ordinary differential equations.
Fig. 2. M-decomposition of model in Example 1.
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
2.3. Standard form and the reduced graph The form of the representation graph for DAE system models suitable for structural solvability analysis is developed in steps performed sequentially. The first step for investigation of a dynamic process model by the graph theoretical technique described in Section 2.1 is to rewrite the model into the standard form. The transformation of ordinary differential equations into standard form is given in Section 2.2. The second step is the assignment of types to 6ertices in the representation graph (Iri, Tsunekawa, & Yajima, 1972). The important types of vertices corresponding to the model specification are the following. S (set) -type 6ariablesThese represent variables, which are assigned to the specified given values. These variables require no computation and are present in the specification associated to the process model. The type S is usually assigned to initial vertices of the representation graph.In case of dynamic representation graph we can assume an explicit Euler method (or another one-step explicit method) for solving the differential equations. Therefore, the differential variables will be labeled by type S because their initial value can be obtained from the initial conditions then their values can be calculated step by step by numerical integration. In order to distinguish the label S assigned to differential variables from the label S of other variables with specified given values the label of differential variables is denoted by S*. Labels S and S* are treated the same way in the rest of this paper. G (gi6en) -type 6ariablesA variable assigned to a specific value of a left hand side is a G-type variable. Unlike in case of S-type variables the values of the right hand side variables will be suitably adjusted so as to preserve the equality of the two sides. The type G is usually assigned to terminal vertices of representation graph. A non-terminal vertex assigned to type G can be split into two copies one of which being initial and the other which being terminal vertices. The initial vertex is labeled as type S and the terminal as type G. As a result, vertices of type G form a subset of the terminal vertices. According to the representation graph the value of every variable, which have incoming arcs only from vertices labeled by type S can be calculated by simple substitution into the corresponding equation. These variables become secondarily labeled by type S. After this we can calculate the value of all variables, which have incoming arcs only from vertices primarily and secondarily labeled by type S by simple substitution. These variables will be tertiarily labeled by type S and this process can be repeated if necessary.
1637
Omitting all vertices labeled primarily, secondarily, etc. by type S and all arcs starting from them from the representation graph we obtain the reduced graph. The classification of vertices of a reduced graph is as follows. All initial vertices form the unknown variable set X; all terminal vertices labeled by type G constitute the known variable (parameter) set Y; all other vertices constitute the unknown variable set U.
2.4. Differential algebraic equations and the differential index The degree of difficulty to solve a DAE-system F(z(t), z%(t), t)= 0 can be characterized by the differential index of the model. The definition of the differential index (Brenan, Campbell, & Petzold, 1989) is as follows. The minimum number of times that all or part of a DAE-system F(z(t), z%(t), t)= 0 must be differentiated with respect to time in order to determine z% as a continuous function of z, t is the differential index of the DAE-system. Numerical solution of DAEs includes both initialization and integration. To solve DAEs successfully the initial values must be consistent. In case of ODEs the initial values can be chosen independently hence the so called dynamic degree of freedom is equal to the number of differential equations. In contrast to this the initial values of DAEs can be constrained by the algebraic equations, so the initialization of DAEs can be cumbersome and the problems related to initialization increase with higher index values. Difficulty of numerical integration of DAEs also increase with higher index values (Hairer, Lubich, & Roche, 1989). Dynamic process models can be described by semiexplicit DAEs as follows (Hangos & Cameron, 2001). z%1 = f(z1, z2, t), 0 = g(z1, z2, t)
z1(t0)= z10
(17) (18)
According to the definition of the differential index a semi-explicit DAE has index 1 if and only if one differentiation is sufficient to express z %2 as a continuous function of z1, z2 and t. One differentiation is sufficient only if the Jacobian matrix gz 2 =
n (gi (z2j
is non-singular. Structural non-singularity of gz 2 i.e. the full structural rank of gz 2 results in a model, which is structurally semi-explicit index 1. The initial values of the differential variables in the semi-explicit index one models may be chosen independently.
1638
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
index of dynamic process models. We demonstrate the methods and principles first on an introductory example and then present the general treatment.
3.1.1. An introductory example Let us consider the following example. Example 2. Consider a liquid tank system with one inlet stream and one exit stream as is shown in Fig. 4. Let the vessel be perfectly stirred. Heat is transferred to the liquid using a heater. The flow rate and the enthalpy of inlet stream as well as the flow rate of the outlet stream and the heat transfer rate are functions of time. The model equations of the above system are the following. Fig. 4. A simple liquid system.
If gz 2 is singular then the differential index is greater than 1. Structural singularity of gz 2 is a sufficient condition for differential index being \ 1 numerically. In practice most of dynamic process models having differential index \1 are semi-explicit index 2 models, but there are models having arbitrarily large differential index (Moe, 1995). The initial values of the differential variables in the higher index semi-explicit models cannot be chosen independently. High index models can either be simulated directly by an integrator, which tackles high index DAEs or be transformed to semi-explicit index 1 and then integrated. The principles of the reported index reduction algorithms (Moe, 1995) are the following. Removing algebraic equations from the model, that is to transform sets of DAEs to ODEs; removing differential equations from the model by replacing some of the differential equations of the original model by algebraic equations; increasing the number of equations and variables, that is introducing into the model a new algebraic variable by replacing one of the derivatives by a so called dummy derivative. 3. Investigation of structural solvability of dynamic process models The basic notions and techniques of structural solvability analysis is applied to lumped dynamic process models in this section. The effect of model specification on the solvability is analyzed and an index reduction method is proposed based thereon.
3.1. Relationship between differential index and representation graph of dynamic process models There are several examples in literature, which show the effect of model specification on the differential
M% = −L +F
(19)
U%= − L · hL + F · hF + Q
(20)
hL = f1(TL, p)
(21)
hF = f2(TF, pF)
(22)
U= M · uL
(23)
uL = f3(hL, p)
(24)
L=f4(M)
(25)
Here M denotes the mass; U the internal energy; Q denotes heat transfer rate; F and L are inlet and outlet mass flow rate, respectively; hF and hL are specific enthalpy of inlet and outlet flow, respectively; uL is specific internal energy; TF temperature of inlet flow; pF pressure of inlet flow; TL, temperature in the vessel; p is atmospheric pressure and f1, f2 and f3, are given functions. The standard form of this model is obtained the way as it is described in Sections 2.2 and 2.3, which is as follows.
& &
M= M%
(26)
U= U%
(27)
M%= − L+F
(28)
U%= − L · hL + F · hF + Q
(29)
hL = f1(TL, p)
(30)
hF = f2(TF, pF)
(31)
U M
(32)
u*L = f3(hL, p)
(33)
uL =
s= uL − u*, L
s= 0
(34)
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1639
Fig. 5. Representation graph of model with Specification 1 in Example 2.
L =f4(M) (35) Introduction of new variables u *, L s and the equation s =uL −u *, L s = 0 into the original model guarantees that all variables on the left-hand side of the model are different. In Eq. (34) the ‘satelite’ note s =0 refers to the fact that the value of variable s is set to 0, i.e. the variable vertex s is labeled by type G. The representation graph of the liquid tank system is shown in Fig. 5. We indicate differential variables on the static graph at time step t+1 only, but the structure of this graph is the same as the structure of static graph at time step t. Furthermore, consider the following specification. Specification 1. Given the properties of the inlet flow F, TF, pF and the heat transfer rate Q as functions of time; the initial values of mass and internal energy M0, U0 and the pressure p are constants. To be calculated the mass, the internal energy and temperature of the liquid in the vessel and the outlet flow rate M, U, TL, L as functions of time. It can be shown that the differential index of this semi-explicit DAE system is equal to 1. Dynamic degree of freedom is 2, so initial values of mass M and internal energy U can be chosen independently. The types of variables assigned corresponding to Specification 1 are shown in Fig. 5. The reduced graph according to Specification 1 is shown in Fig. 6. It can be seen that there exists a Menger-type complete linking on the reduced graph hence the model is conditionally structurally solvable.
Next let us consider the following specification. Specification 2. Given F, TF, pF, TL as a function of time; M0, U0, p constants. To be calculated the mass, the internal energy of the liquid in the vessel, the heat transfer rate and the outlet flow rate M, U, Q, L as functions of time. It is easy to show that the model with Specification 2 has semi-explicit index 2, because the unknown variable Q is not present in the algebraic equations, so two differentiations are needed to express Q%. The dynamic degree of freedom is 1, so initial values of mass and internal energy cannot be chosen independently.
Fig. 6. Reduced graph of model with Specification 1 in Example 2.
1640
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Fig. 7. Representation graph of model with Specification 2 in Example 2.
The types of variables, which are assigned corresponding to Specification 2 are shown in Fig. 7. Based on the structure of the representation graph on Fig. 7 the following structural solvability properties can be determined. 1. The standard form model is not structurally solvable (there is no Menger-type complete linking on the graph). 2. There are an o6erspecified part and an underspecified part on the representation graph. The overspecified part indicates the fact that the initial values of the model cannot be chosen independently. The underspecified part Q U% indicates, that Q cannot be calculated from the algebraic equations, so the index of the model is \ 1.
3.1.2. General treatment Inspired by the introductory example the following general statements can be established. The L-decomposition of a reduced static graph results in two types of L-components. Definition. An L-component is essential if it contains at least one vertex labeled by G (the corresponding variable belongs to the known variable set Y). An M-component is essential if it contains at least two vertices. If an L-component is not essential then it contains one vertex only and in this case the corresponding variable can be calculated with simple substitution. Vertices belonging to a Menger-type linking on the reduced graph are always in essential L-components.
Theorem. Consider a dynamic process model M in the standard form according to Sections 2.2 and 2.3 described by a semi-explicit DAE system. Assume that the generality assumption GA is true. Then the differential index of M is equal to 1 if and only if there exists a Menger-type complete linking on the reduced graph at any time step t. Proof. Consider a semi-explicit DAE system in the following form. z%1 = f(z1, z2, t) 0 = g(z1, z2, t)
z1(t0)= z10
(36) (37)
The differential index of this model is equal to 1 if and only if gz 2 is non-singular. The structure of dynamic graphs described in Section 2.2 yields that the derivative variables are missing from the reduced graph or they are not in the essential L-component of reduced graph. Assuming an explicit method for solving ordinary differential equations all differential variables can be labeled by type S* on the static graph at any time step t. Hence the variables belonging to z1 do not appear on the reduced graph. Therefore, all unknown variables of essential L-components belong to z2. So the fact that gz 2 is non-singular is equivalent with the existence of the Menger-type complete linking on the reduced graph. Consequence. Under the assumptions of Theorem above, the following statement holds. If the differential index of M is greater than 1 then there is no Menger-
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1641
type complete linking on the static graph at any time step t.
hF = f2(TF, pF)
(45)
U= M · uL
(46)
The properties of a static graph of a dynamic model, which has differential index \ 1 are as follows. 1. The fact that the initial values of integral variables cannot be chosen independently results in an overspecified part on the graph. This situation can be easy shown by assignment of types to vertices corresponding to the model specification. There is an overspecified part on the graph if a vertex labeled by type G can also be labeled secondarily or tertiarily or etc. by type S. 2. Non-singularity of gz 2 results in an underspecified part on the graph. In this part those algebraic variables appear, which cannot be calculated from algebraic equations and those differential variables, which we want to calculate from them.
uL = f3(hL, p)
(47)
L=f4(M)
(48)
3.2. Index reduction on the representation graph There is a natural method for index reduction suggested by the structure of the representation graph of semi-explicit DAE-systems. Example2. (continued) Consider the Example 2 with Specification 2. According to Section 3.1 this model has differential index 2. Modify the standard form of the model as follows. Instead of equation U = U% introduce the following equation. 1 U%= (U t + h − U t) h where U t and U t + h denote the value of internal energy U at time step t and t + h, respectively, and h is the step length. This equation corresponds to numerical deri6ation so if we get the derivative U% by numerical derivation we can calculate Q from it. In order to simplify notation, let us scale time locally such that h=1 holds. Then the new equation replacing U = U% is the following in the scaled time. U% = U t + 1 − U t
(38)
According to the above transformation of the standard form let the modified standard form of model of Example 2 be the following.
&
M = M%
(39)
M%= − L+ F
(40)
U%= − L · hL +F · hF +Q U%*=U
t+1
−U
s = U% − U%*, hL = f1(TL, p)
t
s =0
(41) (42) (43) (44)
During the modification the internal energy U turns from a differential variable into an algebraic variable and U% turns from a derivative variable into an algebraic one. An important consequence of the modification is that there is no integration for variable U in the modified model, so it is not necessary to gi6e initial 6alue for the 6ariable U. The representation graph of the modified model can be seen in Fig. 8. The static graph at time step t+ 1 shows vertices belonging to variables only, which are necessary to calculate Q in time step t. The reduced graph and the L- and M-components of this part of calculation is shown in Fig. 9. It can be seen that there exists Menger-type complete linking on the reduced graph, so the modified model is conditionally structurally solvable. Inspired by the above example we generalize the idea for a wide class of semi-explicit DAE systems and formulate a general index reduction algorithm. Consider a dynamic representation graph belonging to the standard form of a model M described by a semi-explicit DAE-system. Assume that an underspecified and overspecified part can be found on its static graphs belonging to any time step t. Moreover, the overspecified part contains a differential variable while the underspecified part contains a derivative variable. Assume that the underspecified part contains only one input vertex, i.e. it is under specified by once (this is the most usual form of the models with differential index larger than 1).
3.2.1. Index reduction algorithm 1. Let us form the following variable sets. I0 is the set of the differential variables belonging to the overspecified subgraph, D0 is the set of the derivative variables referring to the differential variables of set I0 I1 is the set of differential variables from which directed paths lead to the derivative variables in the set D0, D1 is the set of derivative variables referring to the differential variables of set I1,… Ik is the set of differential variables from which directed paths lead to the derivative variables in the set Dk − 1, Dk is the set of derivative variables referring to the differential variables of set Ik …
1642
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Fig. 8. Representation graph of modified model in Example 2.
2. Let n be the smallest natural number for which the set Dn contains some deri6ati6e 6ariables of the underspecified subgraph. Then the differential index of the model is wd =n+ 2
(49)
If there is no such number n then the model is not structurally solvable. 3. Based on the sets Ik, Dk (k =0, …, n) let us form a finite series of differential and derivative variables as follows: x0, x %0, x1, x%1 …, xn, x %n, where xk Ik (k=0, …, n), there is a directed path from xk to x − %k l (k =1, …, n). 4. Using the series of variables constructed in item 3 let us modify the model as in Example 2 with Specification 2 before. The modification sub-steps can be seen in Fig. 10. This modification essentially means one index reduction step of the model. As a result of transformations the differential index of the model decreases by one in every iteration step and there is an index 1 model obtained at the end of this iterative index reduction process.
Remarks. There are some important remarks about the above general index reduction algorithm, which are as follows.
1. In case of index 2 models the algorithm presented in Fig. 10 includes only one iteration step. Then the above mentioned differential variable x0 can be determined using the simple condition as follows. let x0 be on the overspecified subgraph, and let x %0 be on the underspecified subgraph. 2. If we do not take into account the fulfillment of the generality assumption GA defined in the Linkage theorem then only the structural differential index and the conditional structural solvability of the model modified by the index reduction algorithm can be investigated. 3. The steps in the algorithm can all be performed using the equivalent algebraic occurrence matrix description of the representation graph augmented
Fig. 9. Reduced graph of modified model in Example 2.
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1643
dT K = K1(Tin − T)− K2 · K3 · C exp − 4 dt T − K5(T−Tc)
Fig. 10. The index reduction algorithm.
with the vertex labels appearing as the symbolic values of the entries in the occurrence matrix. Moreover, every step can be executed in polynomial time including the partitioning of the graph into L-components, therefore, the index reduction algorithm is a polynomial time algorithm. These properties enable to use the algorithm for large problems with several hundred variables, as well. 4. The idea of using the ‘reverse information flow’ for computations, which is the essence of the index reduction algorithm, has appeared in a pioneering work by Smith (1985) but in another context.
(50)
K4 dC = K1(Cin − C)− K3 · C exp − dt T
(51)
Cin = Cin(t)
(52)
Tin = Tin(t)
(53)
Tc = Tc(t)
(54)
where C is the reactant concentration in the reactor; Cin is the reactant concentration in the feed; T stands for the temperature in the reactor; Tin the temperature in the feed; Tc the temperature of cooling water and Ki (i= 1, …, 5) are constants. Specification 1. Given Cin, Tin and Tc, as functions of time; K1, K2, K3, K4, K5 and C0, T0 (initial values) constants. To be calculated C, T, as functions of time. The differential index of DAE system above with Specification 1 is equal to 1 while the substitution of Eqs. (52)–(54) into the differential equations results in an ODE system of which the differential index is equal to 0. Assume now that the modeling goal is changed and the modeler now wants to calculate the values of cooling water Tc(t) which makes the reactant concentration C satisfy the specification C(t)= b(t). The specification for this case is as follows. Specification 2. Given Cin, Tin, and C, as functions of time; K1, K2, K3, K4, K5 and T0 (initial values) constants. To be calculated T, Tc, as functions of time. It can be shown that the differential index of this system with Specification 2 is equal to 3 (Moe, 1995). Now let us apply the method proposed in Section 3.2 for the determination of the structural differential index and for the index reduction. First we need the standard form of the model, which is as follows.
& &
3.3. A further example
T = T%
(55)
In this section there is a further example to demonstrate the proposed index reduction algorithm on a simple high index dynamic process model.
C= C%
(56)
Example 3. Consider a well-stirred reactor equipped with a cooling jacket. A first order exothermic chemical reaction takes place in the reactor. We assume that the volume of the liquid phase mixture is constant in the reactor and the physico-chemical properties are constants. The system can then be described by the following model equations (Moe, 1995).
(57)
T%= K1(Tin − T)− K2 · K3 · C exp −
K4 − K5(T−Tc) T
C%= K1(Cin − C)− K3 · C exp −
K4 T
(58)
The dynamic representation graph belonging to the standard form can be seen in Fig. 11. The assignment of labels to variables is selected according to Specification 2.
1644
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
It can be shown that there is no Menger type complete linking on the representation graph, i.e. the standard form is not structurally solvable. There is an underspecified and an overspecified subgraph on the static graphs. As it is seen from the standard form equations, the reason of the overspecification is that the variable C has to be determined by integration (the label S* refers to that) but the value of C is given at every time step from elsewhere (the label S refers to that). The variable sets defined in step 1 of the index reduction algorithm are determined from the representation graph as follows. I0 = {C}
(59)
D0 ={C%}
(60)
I1 = {T, C}
(61)
D1 ={T%, C%}
(62)
Since the variable T% is present in the underspecified subgraph, n= 1 and wd =n + 2 = 3. The variable series, which can be used for the model transformation is C, C%, T, T%. Therefore, the index reduction can be done in two iteration steps. As first step we replace the equation C = C% for variable C with the equation C%= C t + 1 −C t
(63)
in the model and cancel the label S* of C. The representation graph of the resulted model after the first iteration step can be seen in Fig. 12. In the second iteration step we replace the equation T= T% for variable T with the equation
T% = T t + 1 − T t
(64)
and cancel the label S* of T. The representation graph of the modified model can be seen in Fig. 13 while its reduced graph is exhibited in Fig. 14. There is now a Menger-type complete linking on the reduced graph, i.e. after the two iteration steps the model is conditionally structurally solvable. The following two facts indicate that the differential index of the original model is equal to 3. The first one is that two equations of x = x% type have to be dropped from the original model. The other is that for the calculation of the values of all variables of the modified model at time step t the values of some variables at time steps t+1 and t+ 2 have to be used, as it can be seen in Fig. 13.
4. Conclusion In this paper a graph-theoretical graphical technique is presented, which is suitable to compute the structural differential index of lumped dynamic process models as well as to determine the structural decomposition of the model equations. The characteristic features of the representation graph of index 1 and high index semi-explicit DAE models are also described. Based on the structural solvability properties a novel index reduction algorithm is proposed for both index 2 and higher index models. The developed algorithm can also be used for investigating the structural differential index and the structural solvability of these models.
Fig. 11. The representation graph of the standard form model in Example 3.
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
1645
Fig. 12. The representation graph of Example 3 after the first iteration step
Fig. 13. The representation graph of Example 3 after the second iteration step.
It is important to note that the applied model transformation method for index reduction, i.e. to re-formulate ODEs with ‘reverse flow’ of information as difference equations, cannot only be used in the framework of structural analysis using a representation graph but also in connection with other analysis methods based on algebraic matrix representation, too.
Acknowledgements This research was supported by the Hungarian National Research Fund through grant no. T-026575, which is gratefully acknowledged. We also thank the helpful comments and criticism of one of the anonymous reviewers, which inspired us to improve the paper substantially.
1646
A. Leitold, K.M. Hangos / Computers and Chemical Engineering 25 (2001) 1633–1646
Fig. 14. The reduced graph of Example 3 after the second iteration step.
References Brenan, K. E., Campbell, S. L., & Petzold, L. R. (1989). Numerical solution of initial 6alue problems in differential –algebraic equations. New York: North-Holland. Duff, I. (1972). Analysis of sparse systems. D.Phil, Oxford University. Hairer, E., Lubich, C., & Roche, M. (1989). The numerical solution of differential– algebraic systems by Runge –Kutta methods. Berlin: Springer.
Hangos, K. M., & Cameron, I. T. (2001). Process modelling and model analysis. New York: Academic Press. Iri, M., Tsunekawa, J., Yajima, K. (1972). The Graphical Techniques Used for a Chemical Process Simulator ‘JUSE GIFS’, Information Processing 71 (Proc. IFIP Congr. 71), vol. 2. (Appl.), 1150 –1155. Krogh Jensen, A. (1998). Generation of Problem Specific Simulation Models Within an Integrated Computer Aided System. Ph.D. thesis, Danish Technical University. van Leeuwen, J. (1990). Handbook of theoretical computer science, Algorithms and complexity, vol. A. Amsterdam: Elsevier-MIT Press. Leitold, A., Hangos, K.M. (1998). On algebraic model transformations. Proceedings of third IEEE workshop on computer-intensi6e methods in control and data processing (pp. 133 – 138), CMP’98 Prague, Czech Republic. ModKit (2000). Computer aided process modeling (ModKit) http:// www.lfpt.rwth-aachen.de/research/modeling/modkit.html. Moe, H. I. (1995). Dynamic process simulation, Studies on modeling and index reduction. Ph.D. Thesis, University of Trondheim. Murota, K. (1987). Systems analysis by graphs and matroids. Berlin: Springer. Sargent, R. W. H., & Westerberg, A. W. (1964). Speedup in chemical engineering design. Transactions of Institute of Chemical Engineers, 42, 190 – 197. Smith, G. J. (1985). Dynamic simulation of chemical processes. Ph.D. Thesis, Cambridge University. Stephanopoulos, G., Henning, G., & Leone, H. (1990). MODEL.LA a modeling language for process engineering — 1 the formal framework. Computers and Chemical Engineering, 8, 813–846. Unger, J., Kroner, A., & Marquardt, W. (1995). Structural analysis of differential-algebraic equation systems —theory and application. Computers and Chemical Engineering, 19 (8), 867 – 882. Yajima, K., & Tsunekawa, J. (1982). Graphical techniques used for a dynamic chemical process simulation. System modeling and optimization proceedings 10th IFIP conference New York City, 1981. In R. F. Denrick, & F. Kozin, Lecture notes in control and information sciences, vol. 38 (pp. 826 – 833). Berlin: Springer.