Structural analysis for static and dynamic models

Mathematical and Computer Modelling 55 (2012) 1051–1067

Contents lists available at SciVerse ScienceDirect

Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm

Structural analysis for static and dynamic models R. de P. Soares a,∗ , Argimiro R. Secchi b a

Departamento de Engenharia Química, Escola de Engenharia, Universidade Federal do Rio Grande do Sul, Rua Engenheiro Luis Englert, s/n, Bairro Farroupilha, CEP 90040-040, Porto Alegre, RS, Brazil b

Programa de Engenharia Qumica, PEQ/COPPE/UFRJ, Universidade Federal do Rio de Janeiro, Caixa Postal 68.502 – CEP 21945-970, Rio de Janeiro, Brazil

article

info

Article history: Received 1 December 2009 Received in revised form 14 July 2011 Accepted 20 September 2011 Keywords: Structural analysis DAE Debugging Graph theory

abstract Modeling and simulation is an expanding area in several domains. In this work, the currently available methods for detecting inconsistencies in systems of equations coming from both static and dynamic models are briefly reviewed. For the dynamic case, a new algorithm is proposed. The proposed algorithm allows a detailed diagnosis of the dynamic degrees of freedom and valid initial conditions of differential-algebraic equation (DAE) systems. The algorithm is scalable for large problems and is a promising diagnosis tool to spread the usage of equation-oriented dynamic simulators. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction The growth of interest in the modeling and simulation of physical systems in varied domains is evident. One of these domains which is rich in both complex problems and tools is the modeling and simulation of chemical processes. The tools employed to solve the problems in this field are known as process simulators. The current process simulators may roughly be classified into two groups: modular and equation-oriented [1]. In modular tools, the models of process units are precoded in a programming language by a modeling expert and incorporated in a model library. The end-user selects the models from the library and connects them to form the plant model. The plant model is also known as the process flow diagram or simply flowsheet. Usually a graphical interface is available, allowing the graphical construction of the process flow diagram, as shown in Fig. 1. The incorporated engineering knowledge as well as the model structure are largely fixed and not accessible [2]. Frequently, the unit models contain a numerical algorithm for the convergence of its internal variables. Another algorithm is necessary for the convergence of the plant model. This algorithm needs to determine an appropriate set of tear streams in order to build a sequence of nested convergence loops [3]. In equation-oriented (EO) or equation-based implementations, the construction of the plant model is very similar and a graphical interface is also usually available. However, the unit models are written in some descriptive or modeling language and usually are opened for inspection and extension. These models share with the plant model only their equations, not their numerical solution. As a consequence, the implementation of component models is independent of any particular application or algorithm that may be used for their solution. Recognition of potential benefits of the EO technology has led to the development of several tools. Examples of EO implementations are SpeedUp [4], gPROMS [5], ASCEND [6], ABACUS II [7], and EMSO [8]. All these tools are focused on modeling and solving chemical process (chemical engineering) problems. A more general approach is proposed by the Modelica language, which is a new language for hierarchical object-oriented physical modeling that is being developed

∗

Corresponding author. Tel.: +55 51 33083528; fax: +55 51 33083277. E-mail addresses: [email protected] (R. de P. Soares), [email protected] (A.R. Secchi).

0895-7177/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2011.09.030

1052

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 1. Process flow diagram sample (ammonia synthesis process).

Fig. 2. Heat exchanger system with a bypass.

through an international effort [9]. Modelica is designed to give the user the possibility to combine electrical, mechanical, hydraulic, and thermodynamic component models within the same application model. There are some commercial and free simulation environments which implement the Modelica language specification; more details can be found at http://www.modelica.org. Although Modelica is designed for multidomain, to the best of our knowledge, its application in the simulation of chemical process problems is almost nonexistent. Despite this fact, Modelica-based and EO process simulators have to face the same challenge: implement flexible integrated environments with advanced debugging capabilities. Here, the term debug comes from an analogy with software development. We call debugging for any method which aids in detecting and removing problems of models. Methods for diagnosing and fixing ill-posed models for static systems are available in the literature. These methods are reviewed in Section 2. For dynamic models, a new algorithm is proposed in Section 3. The proposed method can convert a dynamic model into an under-constrained static model. This enables the analysis of dynamic degrees of freedom (valid initial conditions) with the well-established methods for the analysis of static models.

1.1. Difficulty when modeling with EO tools In general, modeling is a difficult activity. When the plant model is constructed by connecting the unit models, EO or otherwise, the result is a model with missing information. The user needs to find out which variables need to be specified in order to solve the model. In modular tools, the set of variables which can be specified is usually very limited and the task of closing the degrees of freedom sometimes is aided by some kind of a wizard. On the other hand, in EO tools, the user can make relatively arbitrary specifications with ease. This is one of the greatest strengths and at the same time one of the greatest weaknesses of such systems. To create a consistent model in EO tools involving even as few as 30 equations is difficult [10]. For instance, consider the plant model of a heat exchanger system with a bypass as shown in Fig. 2.

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1053

If the unit models contain only the mass balances, the plant model for Fig. 2 is simply: f1 : y1 = y2 + y3 f2 : y2 = y4 f3 : y3 = y5

(1)

f4 : y6 = y4 + y5 , where y1 , y2 , . . . , y6 are the mass flow rates. Without any additional specification, the system has 6 variables and only 4 equations. In order to simulate this system, the user needs to specify 2 variables (add 2 equations). For this small system, it is clear that the user can specify either y1 and y2 or y1 and y3 but neither y1 and y6 nor y3 and y5 . Actually, for such a small problem, it is easy to enumerate all valid specification options by hand. For larger models, this kind of analysis by manual inspection becomes impossible and both novice and experienced users can have great difficulty in accomplishing this step correctly. In Section 2, we review the methods that can, besides verifying whether a given specification set is feasible, automatically enumerate all valid specification sets. 1.2. Dynamic models and initial conditions Besides the problem of finding a valid set of specifications, for dynamic models, the consistency and completeness of initial conditions must be checked. For instance, consider the following system of equations: f1 (x′1 , x′2 ) = 0 f2 (x2 ) = 0,

(2)

where x′1 and x′2 are the time derivatives of the variables x1 and x2 (the symbol x is being used to distinguish variables with time derivatives from the algebraic ones y). The system (2) contains one differential equation and one algebraic equation. The presence of one differential equation suggests that the system has only one dynamic degree of freedom and requires initial condition for only one variable. Actually, the number of differential equations and the number of dynamic degrees of freedom match only for particular systems (see Section 3). Even though, assuming that the system (2) requires only one initial condition, which variable should be fixed? For this small system, it is clear that x2 cannot be supplied as initial condition, because it is already fixed by f2 (x2 ) = 0. In Section 3, a new algorithm for the analysis of differential-algebraic systems is proposed. Using this algorithm, all valid initial conditions can be enumerated and also the well-posedness of a given initial condition can be checked. 2. Static models A general system of nonlinear algebraic (NLA) equations can be represented by: F (y) = 0,

(3)

where F is the function vector and y is the variable vector. System (3) has m equations in terms of n variables and appears in the solution of steady-state simulations. The equations F represent balances, constitutive equations and specifications. The variables y are associated with the streams (e.g. flows, compositions, temperatures, pressures) and with unit blocks (e.g. sizes and operating conditions). Spatially distributed models are assumed to be discretized in (3). If all equations are linear, which is usually not the case, the system (3) can be reduced to: Ay = b.

(4)

The solution of (3) is usually obtained by a Newton-like method. Such methods iteratively solve problems like (4), where A is some approximation of the Jacobian Fy . This means that Fy needs to be invertible along the solution path. Otherwise, the numerical method will fail. Most codes for matrix inversion are able to detect numerical singularities before failing. Some more robust codes can check if the matrix is structurally singular1 before actually trying to solve the problem, e.g. UMFPACK [11]. However, the linear algebra codes cannot be used to diagnose the source of the problem because they stop the analysis when the singularity is detected. The early works on structural analysis of equation systems were concerned with partitioning large systems into small ones and with testing if a problem is solvable [12]. If a system of equations is solvable, there must be at least one possible complete assignment of variables to equations. This assignment can be accomplished by the well-known algorithms for block lower triangularization to make the structure of the Jacobian Fy become Block Lower Triangular (BLT). A basic step of a BLT partitioning is to permute the equations to make all diagonal elements of the structural Jacobian matrix non-zero. This

1 A property of a matrix is called structural if it is independent from the numerical values of the nonzero elements of the matrix.

1054

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 3. Graph for the NLA system (5).

step can also be seen as a procedure which assigns to each variable yi an equation Fj . According to [12], the first to analyze the structure of systems of equations in this way was Steward in 1962 [13]. In the early works, either graphs or incidence matrices were used to represent the systems. The assignment problem can be used to check if a problem is solvable but cannot help to fix ill-posed models. For this case, Mattsson [14] proposed to add fictitious variables and equations for unassigned equations and variables. Then the fictitious elements are collected in blocks and reported to the user with the aim of giving some hints about what is wrong. Morton and Collingwood [15] proposed an equation analyzer for static models of chemical processes. Sparse storage is used to represent the structure of the Jacobian and then a modification of the Duff algorithm [16] is used to produce the equation-variable assignment. Bunus and Fritzson [17] proposed a methodology for detecting and correcting over-constrained and under-constrained models. This methodology is based on the same previous ideas of equation-variable assignment, where graphs are used to represent the Jacobian matrix. The graph representation is advantageous, mainly because graphs can naturally represent sparse matrices in a very efficient way. A noticeable improvement in the Bunus and Fritzson method is to debug the systems by using the canonical decomposition proposed by Dulmage and Mendelsohn [18] (see Section 2.2) instead of adding fictitious variables or equations as previously proposed by Mattsson. In this section, graph based methods which are able to detect and report the source of matrix structural singularities are briefly reviewed. As structural nonsingularity is a necessary condition for the numerical solution of NLA problems, methods for structural analysis arise as an interesting debugging tool. Later in Section 3, the structural analysis of DAE systems, which is not as mature as the analysis for NLA systems, is visited. Some graph theory basics are given in the Appendix A. 2.1. NLA systems as graphs From the comments above, it is clear that when analyzing systems of equations, the equation-variable assignment is very important. For this combinatorial problem, the concept of bipartite graphs can be used. A graph G = (V = Ve ∪ Vv , E ) is called bipartite if V admits a partition in two classes such that every edge has its ends in different classes. In this work, the partition classes are called Ve (equation vertices) and Vv (variable vertices). Using this notation, the NLA system of Eqs. (5), as studied by [17], can be drawn as the bipartite graph shown in Fig. 3. f1 (y1 ) = 0 f2 (y1 , y2 ) = 0 f3 (y1 , y2 ) = 0 f4 (y2 , y3 , y4 ) = 0

(5)

f5 (y4 , y5 ) = 0 f6 (y3 , y4 , y5 ) = 0 f7 (y5 , y6 , y7 ) = 0. As can be seen in Fig. 3, the values or form of the equations in (5) are irrelevant; only the equation-variable relationship is considered. This is the essence of the structural analysis. 2.2. Debugging NLA systems When checking for defects in NLA systems, the first and easiest test is the degrees of freedom check. But systems with zero degrees of freedom still can be inconsistent. This is the case of system (5) with 7 variables and 7 equations. As cited before, the structural singularity of a NLA system can be checked by a structural block lower triangularization which is equivalent to an equation-variable assignment problem. In graph theory field, this is done by executing a maximum matching algorithm [19]. A maximum matching is an edge subset with as many edges as possible, given that no two edges can have a common end vertex. More details can be found in Appendix A. One maximum matching association for system (5) can be seen in Fig. 4. In this figure, the edges which are part of the matching are in bold and nodes not covered by the association are filled in gray. It should be noted that the optimal assignment or association in Fig. 4 is not unique.

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1055

Fig. 4. Maximum matching for system (5).

Fig. 5. DM decomposition for system (5).

Fig. 6. Directed graph for the determination of the over-constrained partition.

If a maximum matching association includes all variables and all equations (a perfect matching), then the system is structurally nonsingular. This is equivalent to get a zero-free diagonal in a structural lower block triangularization. Otherwise the system is structurally singular and any attempt to numerically solve the system will fail. Although useful for a structural singularity check, the maximum matching cannot be used as a tool for fixing the problem because the source of the problem is still obscured. Furthermore, the maximum matching assignment is not unique and a direct analysis of the result would not generate reliable results. One further step can be achieved by using the canonical decomposition by Dulmage and Mendelsohn (DM decomposition) [18]. This method canonically decomposes any maximum matching of a bipartite graph into three distinct parts: over-constrained, under-constrained, and well-constrained. The application of the DM decomposition to the graph shown in Fig. 4 produces Fig. 5. It is interesting to note that a bipartite graph can admit multiple different maximum matchings. Essentially, the DM decomposition takes into account all possible maximum matchings admitted by the graph. The over-constrained partition is determined by the set of equations exposed (not associated) by all possible maximum matchings. Similarly, the underconstrained partition is the set of variables exposed by all possible maximum matchings. The construction of these sets by enumerating all possible maximum matchings, although apparently complex, is a very simple task. In graph theory jargon, it is just a matter of finding alternating paths starting at the exposed nodes. Practically, we just need to start with one maximum matching M max and then replace the undirected edges by directed edges. For the determination of the over-constrained partition, all edges covered by M max are replaced by directed edges from the variables vertices Vv to the equations vertices Ve . Edges not covered by M max are replaced by directed edges in the opposite direction. The application of this procedure on the graph shown in Fig. 4 will produce the directed graph in Fig. 6. Starting from the equation nodes not covered by the matching, the over-constrained partition is formed by all equation nodes that can be reached by following the direction of the edges. For example in Fig. 6, following the directed edges starting at f1 , only the nodes f2 and f3 can be reached. Then these three nodes are the over-constrained partition, as depicted in Fig. 5. There is a symmetry in the determination of the under- and over-constrained partitions. In a similar way, one just need to invert the direction of the edges and start from the variable nodes not covered by the matching to get the under-constrained partition. The graph with the directions inverted can be seen in Fig. 7. From this figure, it is clear that the under-constrained partition involves only the variables y6 and y7 .

1056

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 7. Directed graph for the determination of the under-constrained partition.

Fig. 8. Directed graph for the under-constrained partition determination of the system (1).

Finally, the well-constrained partition consists of the nodes not included in any of the other two partitions. The procedure described above represents a very simple way to determine the system partitions according to the DM decomposition. This formulation of the problem also makes it clear that the under- and over-constrained partition determination are symmetric problems. Further, we consider the procedure by directional edges explained in this work more intuitive than the one presented in [17] which uses bidirectional edges. From Fig. 5, a debugging tool can conclude that one of the equations {f1 , f2 , f3 } needs to be removed and one additional equation involving y6 or y7 needs to be added. From the end-user perspective of an EO tool, one possible conclusion could be: y7 should be specified and f1 should be removed. Basically we consider two fixing options:

• Remove the extra equations from the over-constrained partition • Add extra equations or specifications to the under-constrained partition In [17], the authors present a comprehensive text about these fixing strategies. However, they consider the removal of variables when an under-constrained partition is found as a valid fixing option. There is no sense in arbitrarily removing a variable from a system and consequently removing variables from equations. This option only makes sense if the variable does not appear in any equation. If a model has an isolated variable it is indeed defective and this should be reported earlier to the unit model developer. Some additional remarks about the fixing strategies follow. 2.2.1. Under-constrained models If a maximum matching algorithm finishes with one or more free variables, then there is missing information in the model. This is usually the case when the end-user has just finished the plant model and then needs to add some specifications. If the DM decomposition is applied, it will reveal an under-constrained partition which will contain the variables not covered by the initial maximum matching and probably much more variables will be involved. For instance, consider again the plant model of a heat exchanger system with a bypass introduced in Section 1 (Fig. 2). Without any additional specification, the system has 6 variables and only 4 equations. A possible maximum matching for this system can be seen in Fig. 8. This figure has the edges directed as previously described for the determination of the under-constrained partition. Following the paths starting at the free variables y5 and y6 , all other variables can be reached. Apparently the DM decomposition could not help in this case. But this result could be reported in a more informative way by reporting separately the subgraphs induced by each free variable:

• Starting at y5 we can reach {y1 , y2 , y3 , y4 }; • Starting at y6 we can reach {y1 , y2 , y4 }. Now let us say that, given the above hints, the user adds a specification for the variable y5 . The plant model is analyzed again and the graph becomes the one shown in Fig. 9. Now the user still has to specify one of {y1 , y2 , y4 , y6 } in order to get a structurally non-singular problem. The complexity of the classical or primitive maximum matching algorithm when it starts from an empty matching is O(nm) (see Appendix A). But the complexity to compute a maximum matching for a graph which was enlarged by adding only one new node and a previous maximum matching is already known is just O(m). So a tool can be very responsive to the end-user as specifications, unit blocks or connections are added/removed even for large scale plant problems.

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1057

Fig. 9. Directed graph for the system (1) with a specification for y5 .

Fig. 10. Directed graph for the system (1) when over-specified.

As can be seen in this very simple heat exchanger problem, the number of fixing options when there is an underconstrained partition can be large. For instance, consider the model of an ammonia synthesis process as shown in Fig. 1 [3]. A static model for this process with 134 variables was constructed using EMSO. This model is found in file ammonia.mso available at http://www.enq.ufrgs.br/alsoc. If all specifications are supplied correctly, the maximum matching algorithm finishes with a perfect matching. But if only one specification is missing, for instance the process feed flow rate, then the under-constrained partition will involve 96 variables. This means that the well-constrained partition covers only about 30% of the variables. Unfortunately, the majority of models have a similar behavior: in the case of a non-empty under-constrained partition, the number of fixing options is quite large. To the best of our knowledge, this is not recognized in the literature and no doubt this is a weak point of the currently available debugging techniques. In order to reduce this deficiency, the fixing options could be ranked by heuristic rules that present the more meaningful options first. These rules are currently the object of study. 2.2.2. Over-constrained models Over-constrained models can be fixed in an analogous way by the analysis of the DM-decomposition results. The simpler case to solve is when the user over-specifies the system. This is a common issue in GUI/canvas style flowsheet modeling. For instance, if instead of adding only 2 specifications for the heat exchanger system (Fig. 2) we add specifications for y4 , y5 and y6 , being f4 , f5 and f6 , respectively. A possible maximum matching for this case is shown in Fig. 10. Following the directed edges in Fig. 10, we conclude that the over-constrained partition contains {f4 , f5 , f6 , f7 }. One of these equations should be removed. We know that f4 is a mass balance and should not be removed, but this information cannot be automatically detected only by the graph structure. In [17], the authors try to overcome this kind of problem by annotating the equations with different flexibility levels. Suppose that the mass balances should be annotated with a low flexibility level and should not be considered for removal if there is some other candidate with a higher flexibility. This transfers the problem to the unit model developer and clearly an automatic way to rank the candidate equations for removal is desirable. A very simple way to do this is to advise the end-user to first remove the equations which are specifications. This can drastically reduce the number of fixing options and will be effective if the unit models are not internally inconsistent. But if some unit model is internally defective, the developer will be probably prompted with several fixing options. Another simple approach which can be explored is the automatic association of the flexibility level with the number of edges of the equation. The lower the number of edges, the higher is the flexibility level. This is effective for small problems but its practical applicability still has to be studied. It should be noted that a system can have simultaneously an under-constrained partition and an over-constrained partition. Even in this case, there is no impact in the approaches presented above. The two defects can be analyzed and reported to the user simultaneously. 2.2.3. Other structural limitations Structural analysis is limited. Several structurally nonsingular problems have no numerical solutions. For instance, the result of the structural analysis for the heat exchanger problem (1) allows the user to specify simultaneously the input flow y1 and the output flow y6 . It is easy to see, by a plant global mass balance, that y1 and y6 are identical. The specification of these variables can produce a perfect matching, as can be seen in Fig. 11.

1058

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 11. Graph for the heat exchanger system (1) when wrongly specified.

There are other simple procedures that can ruin the structural analysis result. For instance, the pre-multiplication of a model F (x) = 0 by a full-rank random matrix R as RF (x) = 0. Each equation of the resulting system RF (x) = 0 will be given by the sum of all equations of the original system F (x) = 0 weighted by a row of R. This will produce little effect in the numerical solution of the system but will introduce several non-zero elements on the Jacobian matrix making any structural analysis result useless. 3. Dynamic models Differential-Algebraic Equation (DAE) systems arise naturally when dealing with dynamic simulation in EO tools. A general DAE system can be represented by: F (t , y, y′ ) = 0,

(6) ′

where t is the time and y are the derivatives of y with respect to t. When convenient, the variables vector y can be further partitioned into variables which actually contains its time derivative x and inputs u. The index of DAE systems is of great importance in the numerical classification of a DAE system [20]. Early works on the structural analysis of DAE systems addressed the problem of index and dynamic degrees of freedom determination [21,22]. In [21], the authors have demonstrated how to compute the structural index. Later, Pantelides [22] developed a graphbased algorithm for the analysis of DAE systems. Today, this algorithm is the most widely used structural technique for the analysis of DAE problems and it has been incorporated in several professional simulation tools like gPROMS [5] and Dymola [23]. The main objective of [22] was to determine the number of initial conditions required for consistent initialization of general DAE systems. Although the algorithm can discover the dynamic degrees of freedom (in a structural sense) there is no hint about which variables should be specified as initial conditions. Unfortunately, the DM decomposition cannot be applied to the resulting graph of Pantelides’ algorithm. Mattsson [14] tried to apply the well-known assignment methods for the structural analysis of NLA to DAE systems. The proposal was to consider any appearance of a derivative y′i as the appearance of the variable yi so the same method for NLA problems can be used. The usefulness of this technique is very limited, as no information about the index or about the dynamic degrees of freedom is obtained. In [24], the authors review the theory of DAE systems and a new structural algorithm is developed. This algorithm, which is conceptually different to the one proposed by Pantelides, is derived from a symbolical algorithm for index reduction proposed earlier by Gear [25]. The authors named their implementation ALGO. While Pantelides’ approach is fairly well suited for sparse systems (graph based), the approach in ALGO is more complex if handled in sparse matrix form and no sparse implementation is available. In [24], the authors have also shown how to examine whether a specified initial condition set allows consistent initialization. However, as noted in [24], neither ALGO nor Pantelides’ approach can terminate if the structural pencil (Fy and Fy′ structurally added, which is simply the union of the matrices) is structurally singular. This is a serious weakness of these algorithm; they stop their execution if a structural rank deficiency is encountered, disabling any diagnosis for this kind of ill-posed model. In [26], Murota showed how to determine the structural index by using maximum weight matchings. In [27], Pryce described another method for the structural analysis of DAE systems. The method can also be seen as a maximum weight assignment problem. A signature matrix is constructed based on the structure of the problem and then a highest-value transversal is solved. This is equivalent to the solution of a maximum weight assignment problem. This method is coupled with automatic differentiation to initialize and solve high-index DAE problems by Taylor series in the numerical package DAETS [28]. DAETS is competitive with standard DAE integrators only if very high-accuracy is required or for the solution of very-high index problems. Further, a sparse implementation is still not available limiting the use to problems with at most some hundreds of variables. In this work, a new algorithm for the analysis of DAE systems is described (Section 3.1). This algorithm is also based on the construction of assignments using graph theory. It can be seen as an extension of the classical algorithm proposed by Pantelides [22]. Thus, it is conceptually different from the approach in [24]. This new algorithm also produces results very similar to the algorithm by [27]. However, we see advantages in our graph-based implementation because it is naturally sparse and it handles problems with a singular pencil. Further, the DM decomposition can be directly applied to the resulting graph. This is a great achievement because it enables debugging for dynamic models at the same level as for the static models.

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1059

Fig. 12. Graph for system (2) with the second equation structurally differentiated.

3.1. New DAE analysis algorithm DAE systems can also be represented as bipartite graphs. But in the dynamic case, there are two new concepts: the derivatives of the variables are also considered and the equations can be differentiated by inserting new elements into the graph. In order to illustrate these concepts, consider again the system of Eqs. (2). Fig. 12 shows the graph representation for this system, but with the second equation differentiated with respect to time. As can be seen in this figure, the variables are classified as algebraic (nodes in gray) and differential. Edges for algebraic variables are also shown in gray. It should be noted that, although x1 does not appear in system (2), it is considered on the graph. It is also considered that the time derivative of f2 involves only x′2 and not x2 (this is called structurally linear differentiation in [24]). The motivation for this consideration is that, regardless of the form of the f2 function, its time derivative will always include x′2 , but only for some special cases will it also involve x2 . Further, this is a reasonable consideration because, as we will see later, the motivation to differentiate an equation (in our algorithm) is to determine a missing variable derivative. In this work, we propose Algorithm 1 for the analysis of DAE systems. This algorithm expects as input a graph G representing the system of equations and returns an enlarged graph (with possibly more equations and variables) and its maximum matching M. This algorithm and all other methods presented in this work were implemented in C++ and Java. The Java version is available as accompanying code at http://www.enq.ufrgs.br/enqlib. Algorithm 1 Pseudocode for the proposed DAE analysis algorithm. input G(Ve ∪ Vv , E ) output: G, M 1: for ve ∈ Ve do 2: if not Augment ( G, M , ve ) then 3: mark all colored vk ∈ Ve 4: uncolor all ve 5: if not Augment (G, M , ve ) then 6: return false 7: end if 8: diff all filled vk ∈ Ve 9: end if 10: end for 11: return true The algorithm can start from any arbitrary matching. If it is the first run, we start with an empty matching. Basically, the algorithm will try to find a match for all equations (line 2). The graph  G is defined as the current graph G without the algebraic variables. Naturally, in  G, all edges with one end in an algebraic variable are also excluded. Actually,  G includes only the highest order derivative of each variable. For instance, if the graph G contains x′′1 , both x1 and x′1 will not be considered in  G. If a match in  G is not possible, it tries to find a match in G, i.e. considering all variables (line 5). If the match in G including all variables is not possible, then the system is structurally singular. Algorithm 1 relies on the Augment algorithm (lines 2 and 5). The Augment algorithm is the same classical algorithm used for finding maximum cardinality matchings in bipartite graphs in Section 2 and described in more detail in Appendix A. If the matching could be augmented, Augment returns true, otherwise it returns false. When Augment returns false, the subset of the equations Ve reached by all alternating paths is colored. When this happens in line 2 of Algorithm 1, all colored equations are filled. Equations filled at that step are structurally differentiated with respect to time in line 9. When equations are differentiated, new equations are added to the system and possibly new variables. In order to clarify the application of Algorithm 1, consider again the system of Eqs. (2) and its graph in Fig. 12. For the first equation f1 , line 2 will find the match {f1 − x′1 }. For f2 , line 9 will fail because it cannot be matched in  G (when the algebraic variables are ignored). As a consequence, f2 will be filled at line 3. In line 5 (when all variables are considered), the match

1060

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 13. Graph for system (2) after the first three steps of the algorithm.

Fig. 14. Graph for system (2) after the analysis.

{f2 − x2 } is found. Fig. 13 shows the graph, the equation filled in gray, and the matching in bold at this point (gray edges are only used to distinguish the algebraic variables). At line 9, the new equation f2′ will be added. Finally, the match {f2′ − x′2 } is found at line 9 and the algorithm finishes. The final matching can be seen in Fig. 14. The main advantage of the new algorithm is that its final association also is a maximum matching and then suitable for a DM decomposition. For instance, the under-constrained partition will reveal all variables which can be supplied as initial conditions. Taking the system (2), the under-constrained partition will include only x1 , as it is the only variable not assigned to any equation in the final matching (Fig. 14). Using this information, an EO tool can tell the end-user that the only option for this model is to supply an initial value for x1 . All other variables {x′1 , x2 , x′2 } are discarded from the initial conditions candidates. The structural index is also computed by the algorithm; it is the maximum number of times that an equation was differentiated. For system (2), it is 1. This result matches with the differential index for the system while Pantelides’ algorithm computes a structural index of 2 for this simple system. The algorithms by Pryce and Murota both also correctly detect this system as an index-1 system. 3.1.1. Structurally singular pencil As mentioned before, neither ALGO nor Pantelides’ approach can terminate if the structural pencil (Fy and Fy′ structurally added, which is simply the union of the matrices) is structurally singular. This kind of model should be rejected because a structurally vanishing matrix pencil could not be solved numerically by the available methods [24]. Thus, this should be reported as a model defect, preferably with fixing hints. A simple model, taken from [22], with this kind of defect, is as follows: f1 (x, u1 , u2 ) = 0 f2 (x, x′ , y1 ) = 0 f3 (x, y2 ) = 0

(7)

f4 (y1 ) = 0 f5 (y2 ) = 0, where, y1 and y2 are the outputs, u1 and u2 are the inputs, and x is the state. As introduced in the original context [22], this problem is known as the uncontrollable DAE system. The specification of desired profiles for the output variables (f4 (y1 ) = 0 and f5 (y2 ) = 0) makes this an optimum control problem. The desired solution for this problem will be the profiles for the input variables u1 and u2 which satisfy the desired outputs. The graph in Fig. 15 represents the system (7). If Algorithm 1 is used to analyze the system (7), a structurally singular problem is found and the graph in Fig. 16 is produced. It should be noted that Pantelides’ algorithm cannot detect the structural singularity of this problem and will keep differentiating the same subset of equations ad infinitum. With the proposed algorithm, we can see that there is an over-constrained partition given by {f2 , f3 , f4 , f5 , f3′ , f5′ }. One of these equations should be removed in order to fix the problem, and this can be reported to the user (with the autogenerated

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1061

Fig. 15. Graph for the DAE system (7).

Fig. 16. Graph for the DAE system (7) after the analysis by Algorithm 1.

Fig. 17. Graph for the system (7) with u1 fixed instead of y2 after the analysis by Algorithm 1.

Fig. 18. Graph for the system (7) with u1 fixed instead of y2 after the analysis by Algorithm 1 with directed edges for the DM decomposition.

equations f3′ and f5′ filtered out). Using this valuable information, the user can conclude that this model does not support the specification of two optimum outputs but only one. The natural choice should be to replace one of the desired output profiles by the specification of one input. For instance, instead of y2 , we specify u1 , i.e., f5 (u1 ) = 0. If the model is fixed this way, the proposed algorithm will find a structurally well-posed model as shown in Fig. 17. Then this model can now be further analyzed with respect to the valid initial conditions. As can be seen, the variable u2 ended up exposed by the final association and the missing information must be supplied by the initial condition, as will be described next. 3.1.2. Initial condition candidates As can be seen in Fig. 17, when the system (7) has u1 fixed instead of y2 , the analysis finishes with one exposed variable: u2 . This means that only one initial condition should be specified to completely determine the initial condition. Because u2 is not covered by the maximum matching, it is a candidate for initial condition. But according to the DM decomposition — following the directed edges in Fig. 18 — it is easy to see that the under-constrained partition consists of {u2 , x, y2 , x′ , u′2 , y′2 }. Since the initial condition is usually not given in terms of differentials, the user should specify as initial condition one of { u2 , x , y 2 } .

1062

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 19. Graph for system (8).

Fig. 20. Graph for system (2) after the analysis.

Another simple problem is the example of a closed vessel containing an ideal gas with a heat input u1 (t ) (by means of an electrical heater) and a forced volume u2 (t ) (by means of a piston). In [29], this is the example 4 and was modeled as: f1 : U ′ + PV ′ = u1 (t ) f 2 : V = u2 ( t ) f3 : PV = nRT

(8)

f4 : U − U0 = ncv T , were U is the internal energy, P is the pressure, V is the volume, T is the temperature, and the other are constants. The graph for the system (8) can be seen in Fig. 19. Since this is a dynamic model, the end-user should supply initial conditions. In order to help the user, a simulator should inform the number of required initial conditions and preferably which variables can and which cannot be specified as initial conditions. If system (8) is analyzed by Algorithm 1, the graph shown in Fig. 20 is obtained. It ends with one free variable T , this means that only one initial condition should be specified to completely determine the initial point. Because T is not covered by maximum matching, it is a candidate for initial condition. According to the DM decomposition, Following the directed edges, it is easy to see that the under-constrained partition consists of {P , T , U , U ′ , P ′ , T ′ }. It is unusual to specify initial conditions in terms of variable derivatives (unless steady-state is given as the initial condition, y′ = 0). Hence the initial condition candidate set can be filtered out to just {P , T , U } and this can be reported to the user. 3.1.3. Large scale problems It should be noted that the differentiations executed by the proposed algorithm are structurally linear. This kind of differentiation can be implemented very efficiently (see the accompanying code) and does not depend on the actual form of the equations. Because the proposed algorithm is actually very similar to a maximum matching algorithm (a loop calling the Augment algorithm), it is expected to perform similarly — O(mn) (see Appendix A). In order to check if our implementation of Algorithm 1 is coded properly, we have tested how it performs for a large scale problem of a dynamic model for a distillation process. This model has mass and energy balances for each tray besides thermodynamics and hydrodynamics equations and was also constructed using EMSO. The graph for this model is also found in the accompanying code under the name ugav40.dot. For the case of the separation of isobutane from a mixture of 13 compounds in a 40-tray column, the number of variables is almost 4,000. The computational time required to analyze the dynamic model of this system with different numbers of trays can be seen in Table 1. The results shown in Table 1 were obtained in a Pentium M 1.70 GHz PC with 2 Mb of cache memory running Ubuntu Linux version 6.06. All analyzed cases are well-posed models — this test was not used to test the debugging power of the algorithm. As can be seen in the table, the performance is approximately quadratic as expected. Another good result is that the time required by the analysis is very acceptable for user interaction. With the current available hardware, this time is of the order of one second for systems with 10,000 variables. Moreover, the algorithm can be applied incrementally by adding new equations and variables as the user interacts with the modeling environment. Each time the user adds a unit block to the diagram, new variables and equations are added. If the user adds a specification, a new equation is added. Once equations and variables are added, the graph can be enlarged and a new assignment can be

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1063

Table 1 Time to analyze the dynamic model of a distillation column varying the number of trays. Trays

Variables

Time (s)

Time/N 2 · 109 (s)

20 40 80

2157 3877 7317

0.04 0.14 0.52

9.46 9.58 9.79

Fig. 21. Classical index-3 pendulum problem.

Fig. 22. Graph for the pendulum problem (9).

Fig. 23. Graph for the pendulum problem (9) after the analysis by Algorithm 1 and with the edges directed for the possible initial conditions determination.

constructed starting from the previous one. The graph based algorithms require no modification to start from a non-empty matching. This process can break up the analysis time, making the software even more responsive to the end-user. 3.1.4. High-index problems The proposed algorithm can be applied without modifications to analyze high-index DAE systems even with high-order derivatives. For instance, consider the classical index-3 pendulum problem [22,24] as depicted in Fig. 21. A possible formulation in Cartesian coordinates is: f1 : Tx = x′′ f2 : Ty − g = y′′

(9)

f3 : x + y = L , 2

2

2

where x and y represent the horizontal and vertical positions, T is the cable tension. The constants L and g are the cable length and the gravity acceleration. This system can also be modeled without the second order derivatives by adding additional variables, namely the velocities w = x′ and z = y′ . A graph for (9) can be seen in Fig. 22. When the system (9) is analyzed by Algorithm 1, the graph in Fig. 23 is obtained. The matching in this figure exposes only the variables y and y′ , hence we have 2 dynamic degrees of freedom for this problem. Once f3 is derived 3 times by the algorithm, the system has a structural index of 3.

1064

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 24. Linear electrical network which has the structural index higher than the differential index.

Fig. 25. Graph for the system (10) after the analysis by Algorithm 1.

Following the directed edges emanating from y and y′ , two variable subsets can be constructed: {y, x, T , y′′ , x′′ } and {y , x′ , T ′ , y′′′ x′′′ }, respectively. About the valid initial conditions, the user should supply two, one for each set. From the geometrical constraint x2 + y2 = L2 , it is clear that one cannot specify the initial condition simultaneously for x and y. From the sets above, it can be seen that this situation is automatically filtered out by the structural algorithm. It should be noted that the equations differentiated by the algorithm can also be used to generate an index-reduced system, actually an ODE is always obtained with the final graph, but index reduction is out of the scope of this paper. For more details in this regard, please see [30,31]. ′

3.1.5. When the structural index exceeds the index The fact that structural algorithms can underestimate the index is well known [22,24]. But only recently [32,26], it was found that some kinds of systems can have the structural index arbitrarily higher than the differential index. In [26], there is a mention of this problem as an embarrassing phenomenon. The same is observed with Algorithm 1 as well as any other currently known structural algorithm. As a practical example, consider the linear electrical network shown in Fig. 24, which consists of an independent voltage source v(t ), a resistor of nonzero resistance R, and a capacitor of nonzero capacitance C . As discussed in [32], if node 2 is chosen as the ground node, the equations of the modified node analysis (MNA) are: f1 : C (v1′ − v3′ ) = i f2 : C (v1′ − v3′ ) = v3 /R

(10)

f3 : v1 = v(t ), where vj is the node potential at node j and i is the current through the voltage source. The DAE (10) has differential index 1. If the proposed algorithm is applied to system (10), the graph shown in Fig. 25 is obtained. As can be seen in Fig. 25, equation f3 was differentiated twice, indicating a structural index of 2 which is higher than the differential index. The same is observed by any other structural algorithm available. The question is whether this fact impacts the debugging power of the proposed algorithm. The maximum matching found in Fig. 25 exposes only the variable v3 . Hence, only one initial condition should be specified for the problem, which is correct. Following the directed paths in the graph, the under-constrained partition consists of {v3 , v3′ , v3′′ , i, i′ }. As mentioned before, we could filter out the derivatives from this result to get just {v3 , i}. This is also a correct result; an initial condition for v1 will be redundant with the last model equation f3 : v1 = v(t ). Thus, even when the algorithm over-estimates the index, the debugging results are still correct. 4. Conclusions In this work, techniques which aid in the location and removal of inconsistencies of the models coming from equationoriented simulators were studied. For static models (NLA systems), mature methods were found in the literature and were briefly reviewed.

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1065

For the dynamic case (DAE systems), the algorithms available in the literature can be used to check if a given initial condition set is suitable for consistent initialization or not. However, these algorithms cannot terminate if the structural pencil is structurally singular, disabling any diagnosis for this case. In this context, a new algorithm for structural analysis of DAE systems was proposed. This algorithm terminates irrespective of the pencil being structurally singular or not. But the key advantage of the proposed algorithm is that the resulting graph can be analyzed in the light of the DM decomposition. Thus, it can be used to analyze DAE systems in a way very similar to NLA systems. All methods presented in this work were implemented in C++ and Java. The Java version is available as accompanying code at http://www.enq.ufrgs.br/enqlib. The major deficiency of the studied methods is the large number of fixing options when a singular model is found. In order to reduce this deficiency, heuristic rules for ranking the fixing options are being studied. Furthermore, these codes are being incorporated in the EMSO [8] process simulator. Appendix A. Graph theory basics A graph consists of a pair G = (V , E ) of sets satisfying E ⊆ [V ]2 . Thus the elements of E are pairs of the elements of V . The elements of V are the vertices (or nodes, or points) of the graph G, and the elements of E are the edges (or lines). In Fig. 26, a typical graph is drawn. How the nodes and edges are drawn is considered irrelevant: all that matters is the information about which pairs of vertices form an edge and which do not. A good source for graph theory concepts is [33]. In the graph in Fig. 26, the subset of edges P = {{3, 6}, {3, 7}, {4, 7}} form what is called a path. In this case, the path {6 − 3 − 7 − 4} links or connects the nodes 6 and 4. An important class of combinatorial optimization problem which can be solved using graph theory is the maximum matching problem. Examples of this class of problem are the associations of workers to tasks, classes to classrooms, and so on. This class of problem includes the structural analysis of system of equations, of interest in the present paper, where the equation-variable relationship is studied. In this context, the concept of bipartite graphs is very useful. A graph G = (V = Ve ∪ Vv , E ) is called bipartite if V admits a partition into two classes such that every edge has its ends in different classes. This means that there is no edge linking vertices in the same partition, or in other words, that the vertices on the same partition are not adjacent. In this work, the V partition classes are conveniently called Ve and Vv . The graph in Fig. 26 admits the properties of a bipartite graph, with the possible partition Ve = {1 . . . 4} and Vv = {5 . . . 8}. The maximum cardinality matching problem (or for the sake of brevity, simply maximum matching problem) involves finding a matching in a bipartite graph with as many independent edges as possible. A matching M is a subset of independent edges of a bipartite graph, where independent edges are incident to at most one edge. This means that no edge in a matching can share a node. The cardinality of a matching M is the number of edges matched or covered by M. In the same way, a vertex is said to be covered by M if it is an endpoint of an edge covered by M. Vertices not incident with any edge of M are called unmatched or exposed by M. The classical way of finding a matching in G with as many edges as possible starts with an initial arbitrarily matching M (usually an empty matching). For instance, considering the graph in Fig. 26, an arbitrary initial association could be M = {{1, 5}, {3, 7}}. This association is shown in Fig. 27. The bold lines are the covered edges. The path P = {6 − 3 − 7 − 4} alternates between matched and unmatched edges. Paths with this property are called alternating paths. Furthermore, an alternating path that ends in an unmatched vertex is called an augmenting path because it can be used to turn M in a larger matching. As can be seen in Fig. 27, P is an augmenting path. By dropping from the matching M the edges that are on an augmenting path P and adding to M the other edges of P, we get a new matching, with one more edge than M. This can be seen in Fig. 27. Applying this procedure to the path P with respect to the matching M (cardinality 2), we get a new matching M ′ = {{1, 5}, {3, 6}, {4, 7}} (cardinality 3). We should also note that paths linking two unmatched edges can also be used to augment the cardinality of a matching by one. It is clear that using augmenting paths will always lead, after a number of steps, to a maximum cardinality matching. When a maximum matching covers all nodes of a graph, it is called a perfect matching. If there is no augmenting path, then the current association is optimal and an association with more edges is not possible. This is the base of the theorem developed by Berge [34]. Theorem 1 ([34]). A matching M in a bipartite graph G is of maximum cardinality if and only if G contains no M augmenting path.

Fig. 26. The graph on V = {1 . . . 8} with edge set E = {{1, 5}, {3, 6}, {3, 7}, {4, 7}}.

1066

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

Fig. 27. Increasing the cardinality of M by the augmenting path P.

A depth-first search (DFS) can be used to find augmenting paths starting at a given exposed node, as described in Algorithm 2. This algorithm tries to discover a possible augmenting path. If found, the given matching M is enlarged by one. Algorithm 2 Pseudocode for the primitive augment matching algorithm for bipartite graphs. Augment(G = (Ve ∪ Vv , E ), M , ve ) 1: color ve 2: if exists {ve , vv } ∈ E and vv ∈ / M then 3: M ← M ∪ {ve , vv } 4: return true 5: end if 6: for all {ve , vv } ∈ E do 7: if exists {ve2 , vv } ∈ M and ve2 not colored then 8: if Augment(G = (Ve ∪ Vv , E ), M , ve2 ) then 9: M ← M ∪ {ve , vv } 10: return true 11: end if 12: end if 13: end for 14: return false Remarks about Algorithm 2:

• The algorithm returns true if it is possible to increase the cardinality of the current matching by including a given vertex ve . • A vertex ve can be included directly in the current matching when there is an edge {ve , vv } and vv is not covered by M (line 3); or by an augmenting path (line 9).

• It should be noted that the length of the augmenting path is not limited because the algorithm is recursive (line 8). • In line 1, the ve node is colored, later in line 7, there is a check which avoids any possible infinite recursion. • When the algorithm fails to find a possible augmenting path with an end in ve , all nodes which can be reached by alternating paths starting at the node ve are colored. A maximum matching can be obtained by the sequential application of Algorithm 2. This is the principle of Algorithm 3. The algorithm can start with an arbitrary matching, in a first run it is usual to start with an empty matching. Line 6 assures that all nodes are uncolored in each step because nodes are possibly colored in a previous augment trial. Algorithm 3 Pseudocode for the determination of a maximum matching M max . MaxMatching(G = (Ve ∪ Vv , E ), M ) 1: isPerfect ← true 2: for all ve ∈ Ve do 3: if not Augment(G = (Ve ∪ Vv , E ), M , ve ) then 4: isPerfect ← false 5: end if 6: uncolor Ve 7: end for 8: return isPerfect Algorithm 3 is also known as a classic or primitive algorithm [19]. This algorithm has a O(mn) complexity, where n is the number of vertices and m is the number of edges of G. The literature provides several modifications which can reduce the complexity of the primitive algorithm. For instance, in [35], the authors developed a method which the matching  √augments  in more than one edge per phase. This modification limits the complexity of the algorithm to O m n . In [19], the authors present a comprehensive review of algorithms for maximum matching problems.

R. de P. Soares, A.R. Secchi / Mathematical and Computer Modelling 55 (2012) 1051–1067

1067

References [1] J.F. Boston, H.I. Britt, M.T. Tayyabkhan, Tackling tougher tasks, Chemical Engineering Progress 89 (11) (1993) 38–49. [2] W. Marquardt, Trends in computer-aided process modeling, Computers & Chemical Engineering 20 (6) (1996) 591–609. [3] L.T. Biegler, I.E. Grossmann, A.W. Westerberg, Systematic Methods of Chemical Process Design. Prentice Hall International Series in Phisical and Chemical Engineering Sciences, 1997. [4] C.C. Pantelides, Speedup-recent advances in process simulation, Computers & Chemical Engineering 12 (7) (1988) 745–755. [5] M. Oh, C.C. Pantelides, A modelling and simulation language for combined lumped and distributed parameter systems, Computers & Chemical Engineering 20 (1996) 611–633. [6] B.A. Allan, A more reusable modeling system. Ph.D. Thesis, Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, PA, USA, 1997. [7] J.E. Tolsma, J. Clabaugh, P.I. Barton, ABACUSS II: Advanced modeling environment and embedded simulator. URL: http://yoric.mit.edu/abacuss2/ abacuss2.html, 1999. [8] R. Soares, A. Secchi, EMSO: A new environment for modelling, simulation and optimisation, in: Computer Aided Chemical Engineering, vol. 14, Elsevier, 2003, pp. 947–952. [9] Modelica, Modelica homepage — http://www.modelica.org. URL: http://www.modelica.org, 2010. [10] A.W. Westerberg, D. R. Benjamin, Thoughts on a future equation-oriented flowsheeting system, Computers and Chemical Engineering 9 (5) (1985) 517–526. [11] P. Amestoy, Enseeiht-Irit, T.A. Davis, I.S. Duff, Algorithm 837: Amd, an approximate minimum degree ordering algorithm, ACM Transactions on Mathematical Software 30 (3) (2004) 381–388. [12] R.W.H. Sargent, The decomposition of systems of procedures and algebraic equations, Lecture Notes in Mathematics 630 (1978) 158–178. [13] D.V. Steward, On an approach to techniques for the analysis of the structure of large systems of equations, SIAM Review 4 (4) (1962) 321–342. [14] S.E. Mattsson, Simulation of object-oriented continuous time models, Mathematics and Computers in Simulation 39 (5–6) (1995) 513–518. [15] W. Morton, C. Collingwood, An equation analyser for process models, Computers & Chemical Engineering 22 (4–5) (1998) 571–585. http://www.sciencedirect.com/science/article/B6TFT-3TGSGMN-B/1/8ee359ee60823022d1ce64a47fbb98ed. [16] I.S. Duff, On algorithms for obtaining a maximum transversal, ACM Transactions on Mathematical Software 7 (3) (1981) 315–330. [17] P. Bunus, P. Fritzson, Automated static analysis of equation-based components, Simulation — Transactions of the Society for Modeling Simulation International 80 (7–8) (2004) 321–345. [18] A.L. Dulmage, N.S. Mendelsohn, Coverings of bipartite graphs, Canadian Journal of Mathematics (10) (1958) 517–534. [19] H. Saip, C. Lucchesi, Matching algorithms for bipartite graph. Tech. Rep. DCC-03/93, Departamento de Ciência da Computação, Universidade Estudal de Campinas, Campinas, Brazil, 1993. URL: http://www.ic.unicamp.br/~reltech/1993/93-03.pdf. [20] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problem in Differential–Algebraic Equations, North-Holland, New York, 1989. [21] I.S. Duff, C.W. Gear, Computing the structural index, SIAM Journal on Algebraic and Discrete Methods 7 (4) (1986) 594–603. [22] C.C. Pantelides, The consistent initialization of differential–algebraic systems, SIAM Journal on Scientific Statistical Computing 9 (2) (1988) 213–231. [23] S.E. Mattsson, H. Olsson, H. Elmqvist, Dynamic selection of states in dymola. in: Modelica Workshop. 2000, pp. 61–67. [24] J. Unger, A. Kroner, W. Marquardt, Structural analysis of differential–algebraic equation systems–theory and applications, Computers & Chemical Engineering 19 (8) (1995) 867–882. [25] C.W., Jan. Gear, Differential–algebraic equation index transformation, SIAM Journal on Scientific Statistical Computing 9 (1) (1988) 39–47. [26] K. Murota, Matrices and Matroids for Systems Analysis, Springer-Verlag, 2000. [27] J.D. Pryce, A simple structural analysis method for DAEs, BIT Numerical Mathematics 41 (2) (2001) 364–394. URL: http://dx.doi.org/10.1023/A: 102199862479910.1023/A:1021998624799. [28] N.S. Nedialkov, J.D. Pryce, Solving differential–algebraic equations by taylor series (i): computing taylor coefficients, BIT Numerical Mathematics 45 (3) (2005) 561–591. URL: http://www.scopus.com/scopus/inward/record.url?eid=2-s2.0-28144454918&partnerID=40&rel=R8.0.0. [29] A. Kroner, W. Marquardt, E.D. Gilles, Getting around consistent initialization of DAE systems? Computers & Chemical Engineering 21 (2) (1997) 145–158. [30] R. Bachmann, L. Brüll, T. Mrziglod, U. Pallaske, On methods for reducing the index of differential–algebraic equations, Computers & Chemical Engineering 14 (11) (1990) 1271–1273. [31] R. Soares, A.R. Secchi, Direct initialisation and solution of high-index dae systems, in: Computer Aided Chemical Engineering, vol. 20, Elsevier, 2005, pp. 157–162. [32] G. Reißig, W.S. Martinson, P.I. Barton, Differential–algebraic equations of index 1 may have an arbitrarily high structural index, SIAM Journal on Scientific Statistical Computing 21 (6) (2000) 1987–1990. electronic. [33] R. Diestel, Graph Theory, 2nd ed. Springer-Verlag, New York, 2000, http://www.math.uni-hamburg.de/home/diestel/books/graph.theory/ GraphTheoryII.pdf. [34] C. Berge, Two theorems in graph theory, Proceedings of the National Academy of Sciences 43 (1957) 842–844. [35] J. Hopcroft, R. Karp, An n5/2 algorithm for maximum matching in bipartite graphs, SIAM Journal on Computing 2 (1973) 225–231.