Recomputing Causality Assignments on Lumped Process Models When Adding New Simpliﬁcation Assumptions

: This paper presents a new algorithm for the resolution of over-constrained lumped process systems, where partial differential equations of a continuous time and space model of the system are reduced into ordinary differential equations with a ﬁnite number of parameters and where the model equations outnumber the unknown model variables. Our proposal is aimed at the study and improvement of the algorithm proposed by Hangos-Szerkenyi-Tuza. This new algorithm improves the computational cost and solves some of the internal problems of the aforementioned algorithm in its original formulation. The proposed algorithm is based on parameter relaxation that can be modiﬁed easily. It retains the necessary information of the lumped process system to reduce the time cost after introducing changes during the system formulation. It also allows adjustment of the system formulations that change its differential index between simulations.


Introduction
There are well known reasons for modeling and simulating in engineering: testing and optimizing the design of a system before its construction, avoiding design errors, enhancing quality, reducing costs, evaluating performances, predicting responses, and so on [1]. Over the past 30 years, several modeling tools and environments have been developed [2], such as Modelica [3], Ascend [4], EcosimPro [5] and so on. Most of them are concerned with the model equations describing the system dynamics, which can be directly introduced by the user or automatically generated by the tool from some previous process description or schematic. Many physical systems are often modeled by lumped process models, particularly in the chemical process domain, where partial differential equations of a continuous time and space model are reduced into a finite set of differential-algebraic equations (DAEs) [6]. The general form of a DAE system has the following expression: f(x , x, y, t) = 0 g(x, y, t) = 0 , where the first equation is the differential part of the system and the second one is the algebraic part [7]. The x variables are fitted to a vector of states of n , where x are their time-derivatives and t, the time, is the independent variable. The algebraic variables y, which are not derived, belong to m . Thus, the functions f and g must be defined in some regions such as f:Ω f ⊆ 2n+m+1 → n and g:Ω g ⊆ n+m+1 → m . The set of equations in (1) must have n + m unknowns, x and y, with n + m equations in order to achieve a system with zero degrees of freedom for solvability. Modeling tools sort the equations and identify the sequence for how variables will be solved during simulation. The incidence matrix denotes which variables appear in which equation. For each variable, the matching solution determines the corresponding equation that solves this variable, that is, the causality. The Block Lower Triangular (BLT) transformation is the classical approach for this step [8] together with Tarjan's algorithm [9]. For semi-explicit index-1 DAE models, such as (2), causalization transforms equations into assignments. In the case of high-index DAE systems, index reduction, usually by the Pantelides' algorithm [10], must be performed to achieve an index-1 DAE or an explicit ODE model prior to matching [11]. However, dynamic simulation frequently results in numerical difficulties in the DAE resolution that can arise as a consequence of the differential index and stiffness of the DAE [12,13]. At the modeling level, an index problem can sometimes be avoided by using an alternative model formulation [14]. Therefore, the analysis of the solvability properties is basic for dynamic simulations.
x = f(x, y, t) g(x, y, t) = 0 In the modeling process, the DAE system is usually parameterized by the modeler with some values that are assumed to be known during the simulation. These are the design or input variables which influence the system behavior. The selection of the design variables depends highly on the problem or application area and it may greatly influence the mathematical properties of the model equation such as the differential index. Structural analysis methods can use the visual interface of modeling tools to analyze equation-variable patterns to choose the appropriate variables that have zero degrees of freedom or re-sort the DAEs to improve resolution methods [15,16]. In [17], the effect of modeling decisions on structural solvability are studied, and it was found that a change in the system specification can transform an index-1 model into a high-index system. Structural analysis methods traditionally use matrix representation and linear algebraic tools. More recent developments used a graph and combinatorial approach where the problem is addressed as a maximum weight assignment problem [18][19][20].
In the lumped process, modeling assumptions can be translated into additional mathematical relationships or constraints between model variables, equations or parameters. Adding a new modeling assumption to the system involves an additional equation, and consequently, the set of equations changes and becomes over-constrained. To recover the property of zero degrees of freedom and solve the system, two options can be chosen: removing any equation from the initial systems [21] or relaxing design variables or system parameters. As discussed later, the first option may not be a good alternative due to the logic contradictions of the initial system formulation. The second one converts the selected design variable into a new algebraic unknown to obtain a new system with zero degrees of freedom. The user must decide which design parameters are relaxed and thus, depending on this selection, different solutions can be obtained. This second method does not remove any equation from the original system, only the specification equation corresponding to the relaxed design parameter. In addition, assumption transformations may change the variable structure of the model and induce changes in the causality assignments, which need to be recomputed. Therefore, it is also important to develop algorithms that minimize the computational cost of introducing these assumptions in the initial phases of the modeling tools.
No study has been developed focused on examining the overall effect of modeling assumptions in all modeling processes. In [13], a formal representation of these assumptions is proposed to deduce basic properties of the resultant equations. The simplification assumptions, applied during the modeling phase, may seriously affect the differential index of lumped process models. In [21], adding new model assumptions and their corresponding effects on the differential index of a DAE system are studied by graph transformation. Two algorithms are proposed to follow the changes in the variable-equation assignments and the differential index. However, as will be shown later, Hangos' method may result in contradictions and inconsistent modeling in some cases and therefore, it needs to be improved. In addition, its computational cost can be reduced.
This work proposes such improvement of Hangos' method using the option of relaxing system parameters in the over-constrained system after adding new modeling assumptions. A polynomial time incremental algorithm is developed based on graph theory to analyze the change of differential index in the model simplification phase of the modeling process and reduce the time cost of these assumptions. The proposal can be implemented using either the Bellman-Ford or Dijkstra algorithms together with a representation of the information by minimum heaps or linked list [22,23]. According to that, the best algorithm has an almost linear time complexity if the number of assumption equations added to the system is small. The paper is structured as follows: a brief background on graph theory is given in Section 2. Section 3 describes Hangos' method and demonstrates its limitations by means of an example. The basis of the new proposed method is introduced in Section 4 using the same illustrative example. Then, the proposed algorithm is developed and discussed. Conclusions are summarized in Section 5.

Bipartite Graph Preliminzaries
Next, the minimal background on bipartite graphs used to perform the proposed algorithm is discussed. A bipartite graph, denoted as B(L, R, E), is an undirected graph where the vertex set is decomposed into two disjointed sets, L and R, such that every edge e ∈ E has a vertex in L and R. In this work, DAE systems are represented as bipartite graphs where the first set consists of equation vertices and the second one consists of variable vertices. A matching in a bipartite graph is a set M ⊆ E, such that for every vertex v there is an edge e ∈ M incident to v at most; no edge in a matching can share a node. The cardinality of matching M is the number of edges covered by M. A maximum matching is a matching of maximum cardinality. A perfect matching M is a matching in which for every vertex v there is only one edge of M incident to v. Using graph theory, the equation-variable assignment problem in a DAE system can be done executing a maximum matching algorithm. If the maximum matching association includes all variables and equations (a perfect matching), the system is structurally nonsingular. Given the equation system in (3), the corresponding bipartite graph is shown in Figure 1a where the edges which are part of a perfect matching are in bold. The form of the equations f i in (3) is irrelevant. For structural analysis, only the equation-variable relationship is considered. This optimal assignment is not unique. From this initial matching, an alternative representation, the directed graph, can be derived. In a directed graph C(U, F) induced by the bipartite graph B(L, R, E), the set of nodes is {u 1 , . . . , u n } and (u i , u j ) is a directed edge if and only if i = j and (l i , r j ) ∈ E [24]. The directed graph obtained from the bipartite graph of Figure 1b is also shown on the right.  The problems of finding maximum matching and perfect matching if one exists are two of the most fundamental algorithmic graph problems. Several time polynomial algorithms have been developed. The first developments run in time O(n 5/2 ) [25][26][27]. More recent proposals reduce the time complexity [24,28,29]. Most of them are based on a search for augmenting paths, such as the depth first search (DFS), the breadth first search (BFS) or a combination of both strategies. There also another approaches which employ a push-relabel strategy designed for maximum flow problems [30].
The proposed method is based on a search for augmenting paths. The following Theorems 1-3 establish the importance of the augmented graph in bipartite graphs [23,27,31], and allow the definition of the following Algorithm 1 to compute the matching of maximum weight onto a bipartite graph B(L, R, E) with weight function w [32].   First, a directed bipartite graph is calculated using the weight function z, and two new free vertices s and t are added. In step 3, the augmented path of minimum cost in G M can be computed by finding a shortest path from s to t in G M . There are several classic algorithms for this calculation with polynomial time [22,23]. The Bellman-Ford algorithm computes the shortest path in time O(V·E) which yields an overall running time of O(V 2 ·E). If there are no negative weights, the Dijkstra's algorithm can be used improving the time cost. Using an implementation with Fibonacci heaps, the overall cost is O(V·E + V 2 ·log(V)).

Hangos' Method for New Model Assumptions
In this section, Hangos' method for dealing with new model assumptions is briefly described [21]. Then, its limitations are illustrated by means of a process example. The Hangos method transforms a semi-explicit DAE in the form of (2) into an integral system of equations as follows: where p is the vector of design parameters (or design variables) and x 0 is the initial condition vector of x. They must be specified by the user. Note that the first equation in (2) has been changed to the two first equations in (4), and a new vector of variables u has been introduced. The parameters p and x 0 are specified as constants in the two last equations of system (4), that is, specification equations. For each new assumption equation added to the system, Hangos' method removes another one from the initial system (4) to obtain a new zero degrees of freedom system. Hangos removes algebraic equations, such as the third one in (4), or specification equations. The assumption equations fit into one of the two following types: a constant algebraic variable y i or a steady state differential variable, that is, u i = 0. In [21], two algorithms are introduced to obtain a new matching transformed system from the initial assignment of equation-variable when new assumption equations are included.

Hangos' Algorithms
Let B(L, R, E) be a bipartite graph corresponding to a physical system of equations and variables where L is the set of equations, R is the set of variables and E is the set of edges of the graph B. There is an edge from l ∈ L to r ∈ R if the equation l contains the variable r. Let F be a perfect matching to B. If new assumption equations are added to the system without adding new variables, the task is to find a new matching with the least amount of change in F. Hangos' method is based on writing the system equations in an integral form, reassigning the equation-variable associations by adding the design variables D and state variables S to the graph to achieve a new perfect matching. In [21], Hangos defines two algorithms, depending on the number of equations added to the system:

1.
Hangos' Algorithm 1: it is applied when only one assumption equation e is added. It has the following features: -The free vertexes on the new graph are 'e' and variables from D ∪ S. - The shortest path from 'e' to any variable of D ∪ S is found by a Breadth First Search (BFS). - The resulting matching is that having the greatest number of edges shared with F. - The computational cost is O(V + E), where V is the number of vertex in the graph.

2.
Hangos' Algorithm 2: it is used when several assumption equations Q are added, and its properties are: -The vertexes of the new graph become: The weight function w: E → is defined as follows: The purpose of this function is to provide higher weights to the edges of F to obtain a solution with the maximum number of edges shared with F. - The maximum matching of maximum weight for w has the greatest number of edges shared with F.

Problems of Hangos' Method
In order to analyze and show the problems that can arise using Hangos' method, the open evaporator of Figure 2 is considered. This example is also suggested by Hangos in [13,21]. The DAE system of this model is given by the system of equations in (6), where M is the mass of the fluid in the vessel, U is the fluid internal energy, T is the fluid temperature, E is the amount of fluid evaporated into the environment, Q e is the energy transferred by the fluid to the environment and P* is the fluid vapor pressure. The following variables are assumed to be constant: fluid enthalpy h L , liquid to gas enthalpy h LV , heat transfer coefficient k LV , specific heat C e , atmospheric pressure P 0 , atmospheric temperature T 0 and an invertible thermodynamic function H(.). The design variables are the heat flux Q transferred to the reactor, the inflow F to the evaporator and the outflow L; they must be specified by the modeler.
Equations f 1 and f 2 are the differential parts of the DAE and belong to the input/output mass and heat flows of the system. The equations f 3 , f 4 , f 5  To transform the DAE system of Equation (6) into an integral one, the following new variables are introduced: the overall mass flow of the reactor m = dM/dt, the overall energy flow u = dU/dt, and the initial values M 0 and U 0 for the mass and energy into the reactor, respectively. Therefore, the system becomes: Each variable is assigned to the equation in which it is solved. Therefore, the system causality is 11 -U, f 12 -M 0 , f 13 -U 0 }. All this information is collected into an equation-variable directed graph where every node is labeled as "equation-variable", as it is shown in Figure 3. An edge goes from node 1 to node 2, if node 2 computes a variable which also belongs to the equation of node 1.
Next, two examples are considered adding a new model assumption. These examples illustrate problems and contradictions that can arise with Hangos' method.

Example 1: Removing Initial Conditions
Suppose that M is constant in the integral system of Equation (7), that is, the following assumption equation is added to the system: Each time the user adds a new equation, the graph is enlarged with a node. However, the system has zero degrees of freedom if the number of equations is equal to the number of unknowns. Hangos' method proposes to remove equation f 12 of (7)   Hangos states that the system has a perfect matching; hence, the system differential index is one. Nevertheless, analyzing the new assignments, the following implications are obtained: 1.
From equation f 14 (8), m = 0, and consequently, from f 1 it is obtained that E is a function of F and L: E = F − L ⇒ E = E(F,L).

2.
Similarly, from equation f 3 with P 0 and k LV constant, P* is a function of F and L: From f 4 , and assuming H invertible: From f 5 with T 0 and q LV constant: Q e = q LV · (T − T 0 ) ⇒ Q e = Q e (T) ⇒ Q e = Q e (F, L) .

5.
By equation f 6 with c e and M constant: Finally, from f 2 with h F , h L , h LV constant: Analyzing the last equation, the left side is a function of F and L design variables, and the right side only depends on Q, which is another design variable specified separately. Thus, a contradiction appears because the addition of Equation (8) makes the system over-constrained. To solve this new system, the number of equations must be the same as the number of unknowns. Therefore, another equation should be removed. Conversely, considering this system in a strict way, the initial conditions M 0 and U 0 are considered variables. If this integral system expects to be compatible with a differential one, the time derivatives of these two variables must be zero. Removing the equation f 12 : M 0 = spec implies that the system from Example 1 does not satisfy M 0 = 0 and, therefore, this system is different from the original one.

Example 2: Removing an Arbitrary Specification Equation
Let us suppose that M is constant again; however, the equation f 7 is selected to be removed from Equation (7). This is equivalent to relax the design variable Q, which means that the design variable Q turns into an unknown algebraic variable. The changed equation-variable assignments are given by: f 14 -m, f 1 -E, f 3 -P*, f 4 -T, f 6 -U, f 11 -u, f 2 -Q, as shown in Figure 5. Sorting the equations and highlighting the variables that are computed by each equation in square brackets, the new system can be written as in Equation (9).
Note that the system has a perfect matching, so Hangos' method affirms that the system differential index is one. However, examining equation f 11 , the variable u is under integration. To isolate it, there are two alternatives: 1.
Derive equation f 11 . In this case, equations f 13 , f 6 , f 10 , f 12 , f 4 , f 3 , f 1 , f 8 and f 9 must also be derived. After that, a perfect matching is obtained and the index is two, contradicting the previous statement of index one.

2.
Solve integral f 11 by a quadrature method. In this case, it can be thought that no derivative is necessary. However, applying some rules such as the trapezoidal rule, the following results are obtained: To decrease errors, these two formulas are subtracted: Solving the above integral by the trapezoidal rule: and isolating the term u(t + h), the following expression is achieved: The last expression is the numeric derivative of U(t), which shows that the system solution could not converge. Thus, this second alternative is not a good general alternative.

Proposed Method by Relaxation of Design Variables
From the examples of the previous section, it can be concluded that Hangos' method is incomplete and needs to be improved. Modeling assumptions may induce changes in computational properties of the model and modifications in the variable-equation assignments. If the set of design variables is not changed properly, the index of the model can be increased or the zero degree of freedom will not be satisfied. In the next section, a modification of Hangos' method is proposed. This modification relaxes a design variable for each new assumption equation and restores a system with zero degrees of freedom. Two kinds of model assumptions are considered: steady state equations and algebraic equations. First, the basis of the new procedure is introduced and exemplified. The effect of the model assumptions is described as a graph transformation acting on the structure graph of the model. Then, a new algorithm is proposed based on a graph approach and the fundamentals of this proposal are developed.

Basis of the Proposed Method
This section discusses how Hangos' algorithm can be modified when a new algebraic or steady equation is introduced into the DAE system. The following notation is used: The lowercase letters (u, v, w, . . . ) represent system variables, which are placed to the right of the graph node. In particular, x and s are states and u, v and w denote algebraic variables. • An arrow marked with an asterisk indicates that there is a path from the first node to the last one.
To restore the property of zero degrees of freedom after adding a new equation, it is proposed to reduce in one the number of design variables specified by the user, that is relaxing a design variable by removing its corresponding specification equation. After this modification, two vertices of the associated graph remain without matching edges and a new assignment can be constructed starting from the previous one. This minimizes the analysis time, making the modeling tool even more responsive to the user. Then, a closest maximum assignment is performed. This assignment has the largest number of matching edges in common with the original full assignment.

Steady Equation Addition
When a state variable is assumed at steady state, a new equation d is added with the form d: x = 0. To relax a design variable, there are two possibilities: connecting the new equation d with the selected design variable I with a straight path, or linking the equation through a state variable s.
In the first case, connecting the new equation with the design variable with a straight path (Figure 6a), the new causality assignment is completed by combining every equation with a variable of the next one and combining the last equation with the design variable I, as shown in Figure 6b. The final assignment causality is shown in Figure 6c. Note that this new system has a differential index of one. For instance, let suppose again that the assumption equation f 14 in (8) is added to the original DAE system in (6) in the process example of Section 3.2. The original set of equation is preferred instead of the integral one proposed by Hangos in (7), in order to avoid possible inconsistency problems. The selected design variable to relax is F, so the specification equation f 9 is removed from (6). In this case, there is a direct path from f 14 to F. Figure 7 shows the causality assignment of the initial system in the corresponding equation-variable graph and the changed causality assignments of the modified system which are represented by dotted lines. The new causality assignment is {f 14 The second possibility consists of linking equation d indirectly with the selected design variable I through the state variable s, see Figure 8. There is a differential edge from node s to node c-s represented as a dashed arrow. In this case, an index reduction algorithm, such as Pantelides' algorithm [10], must be applied. This algorithm shows that the path from d to node s must be derived. Equations involved in this path must also be derived. The new system causality is shown in Figure 8c. Thus, in this situation, the differential index is higher than one.
Consider again the previous example of Section 3.2 and suppose now that the selected design variable to relax is Q. Then, the specification equation f 7 is removed. The path from equation f 14 to variable Q is shown in Figure 9. Note that there is a differential edge from node U to node f 2 -U . After applying Pantelides' algorithm to reduce the differential index of the system, the new causality assignment is {f 14  The new system has zero degrees of freedom, a perfect matching and a differential index of one. Thus, the original system has a differential index of two.

Algebraic Equation Addition
Suppose that an algebraic equation d is added to the DAE system with the form d: z = constant. In this situation, the same alternatives as those of previous subsections are found. Figure 10 shows a straight path from d to design variable I and the new causality assignment. Note that the final system has a perfect matching.
In the second case, equation d is linked with a design variable I by a path through a state variable s ( Figure 11). To obtain a perfect matching, it is necessary to apply an index reduction algorithm, such as Pantelides' algorithm, as shown in Figure 11c. Therefore, in this example, the differential index of the system is two.

General Case
The most general assumption appears when an equation with several algebraic variables is added. In this case, it is enough to find a link from any of the paths that begin at vertex d to the selected design variable I. As an example, let us suppose that the equation d is added with the form d: z = y. In Figure 12, two alternatives for finding the design variable that solves this problem are shown.

Proposed Algorithm
The suggested procedure is generalized in the proposed Algorithm 2, that relates algebraic and state variables in the assumption equations with design variables through straight paths or through the state variables of the initial system. This algorithm allows the addition of several equations simultaneously and tries to maintain the highest number of intact equation-graph associations of the initial system as possible, that is, the closest assignment. The algorithm can be combined with another index reduction algorithm to solve the causality assignment in one step. In comparison with the Hangos' Algorithm 2, discussed in Section 2.1, the proposal offers some improvements that optimize its computational cost and behavior. Since a maximum weight matching of a bipartite graph can be calculated in polynomial time [27], the closest assignment problem has the same time bound as valid.
A first improvement concerns the weight function w in (5) used by Hangos. It will be proven that there is a family of functions which obtain the same solution as the w in (5) and that they are simpler to compute.
By arranging all of the expressions to the left-hand side and gathering the expressions that are preceded by a or b, it follows: By (15), the following relation is deduced: This implies: Because a > b, the expression of the first bracket is a positive number. Therefore, the proposed result is obtained.
From Lemma 1, any w ab function can be selected to obtain a maximum matching of maximum weight M, which has the maximum number of element shared with F, considering the set of maximum matchings, that is: ∀ M ∈ {M /M is a maximum matching} ⇒ |M ∩ F| ≥ |M ∩ F|. Theorem 4. Let a, b, c and d ∈ , such that a > b and c > d. If we let B(L, R, E) be a bipartite graph with a matching F, then M is a maximum matching of maximum weight with w ab if and only if M is a maximum matching of maximum weight with w cd .
Proof. ⇒ It is assumed that M is a matching of maximum weight for w ab , and it will be proven that this is also true for w cd . Let M be any maximum matching, and by Lemma 1: By multiplying (19) by c − d, c > d, we obtain: (c − d)(|M ∩ F| − |M ∩ F|) ≥ 0, and therefore: Because M and M are maximum matchings, it follows: By arranging (21): By replacing (22) into (20), it follows: This implies that ⇐ This result is obtained in a similar manner.
Theorem 4 allows taking any w ab function as a weight function, with a > b. The simplest function of this family is the characteristic function of the matching F, which is w F : E → : As a consequence of Theorem 3, taking a matching of maximum weight in a run, another matching of maximum weight is in the next run. Therefore, if F is considered initially, then the number of iterations is drastically reduced.
The next proposition characterizes a maximum matching of maximum weight regarding to the initial matching F.

Lemma 2.
Let B(L, R, E) be a bipartite graph. Let F and M be two matching with |M| > |F|. If |M ∩ F| is the maximum among all matchings with cardinality |M|, then M ∆F has neither cycles nor alternating paths with even lengths.

Proof. (by contradiction)
Suppose for the sake of contradiction that M ∆ F has cycles or alternating paths with even lengths.
• All of the sets p i , d i and c i are disjointed.
• p i represents the augmented paths in F. There are k = |M| − |F| > 0 paths of this kind, and all of them have odd length. Moreover, the number of edges of p i in F + 1 is equal to number of edges of pi in M − F, that is: • d i represents an alternating path of even length. The number of edges of d i in F is equal to the number of edges of d i in M − F: • c i represents cycles, and therefore, their length is even. They fulfill the same equation as d i , that is: If M is a matching in B and its cardinality is |M | = |F| + k = |M|, then: In contrast, (3) According to our hypothesis, if |M ∩ F| is maximum among all the matchings of cardinality |M|, then: |M ∩ F| ≥ |M ∩ F|. By (29) and (30): Equation (31) can be simplified to: By replacing the Expressions (27) and (28) into (32): This expression implies that: |d i | = 0 ∀ i = 1, . . . , r and |c i | = 0 ∀ i = 1, . . . ., s, which is a contradiction.
Theorem 5. Let B(L, R, E) be a bipartite graph. Let F be a matching of B and w F , the characteristic function of F. Let M be a maximum matching of maximum weight for w F . If k = |M| − |F| > 0, then there are k augmented paths disjoined by vertexes {p 1 , . . . , p k } in F, such that: Proof.

1.
It is known by Lemma 1 that the cardinality of M ∩ F is the maximum among all the maximum matchings. After applying Lemma 2, the desired result is obtained.

2.
By (26) and (29) of Lemma 2, it follows that: 3. By the last expression: As k = |M| − |F|, then: The above theorem states the relation between the cardinals of the augmented path of F, F itself and the maximum matching of maximum weight M. Proof. Let {q 1 , . . . , q k } be any k vertex-disjoint augmenting paths of F. The matching M = F ∆q 1 ∆ . . . ∆q k has a cardinality of |M | = |F| + k = |M|. By applying Theorem 5 we deduce that: In contrast: By operating in the last expression, we deduce that Using the previous results, the Algorithm 2 is proposed. It is based on the construction of assignments using graph theory and the Algorithm 1 of Section 2, and is an improvement of Hangos' Algorithm 2. The algorithm expects as input a bipartite graph B(L, R, E) representing the system equations, the weight function w such as that in (25), and an initial matching F before adding modeling assumptions. The complexity to compute the maximum matching starting from an empty matching is O(V·E); however, it is reduced to O(E) when the graph is enlarged by adding only one new node and a previous maximum matching is already known [20]. This algorithm uses the information generated in previous steps of the simulation, so the initial matching F is not discarded, which is an improvement over algorithms that need to start from scratch (these kinds of algorithms have a cost O(V 2 + V·E)). This can make the modeling tool more responsive to user specification changes.
If Q is the number of assumption equations, the time complexity varies depending on the algorithm used in the computation of the path of the lowest cost in G M (step 3.1). When the Bellman-Ford algorithm is employed, the overall time cost is polynomial with O(Q·V·E), since it is performed one time for each new assumption added to the system. The cost of the algorithm can also be improved if it is noted that the edges of G M can take only values from the set {−1, 0, 1}. In [26], the authors proposed an algorithm based on Dijkstra's algorithm with linked list or minimal heaps that improves the time complexity of algorithm 1. This algorithm can be used if the weight of the edges is a positive integer in the range {0, 1, 2, . . . , W}, where W is the highest value of any edge. To apply this algorithm to the proposed one, it is necessary to reduce the set {−1, 0, 1} of the weight function in (25) to a set of a positive integers. The next theorem resolves this limitation. Theorem 6. Let G(V, E) be a graph with weight function w and let k be any real number. Then, M is a matching of maximum weight among all matchings of cardinality |M| for w ⇔ M, which is a matching of maximum weight among all the matchings of cardinality |M| for w 2 = w + k. According to (25), the weight function w of graph G M only takes on values in the set {−1, 0, 1}. By applying Theorem 6 for k = 1, the problem is reduced to another one where the values of the weight function belong to the set {0, 1, 2}. Thus, Dijkstra's algorithms can be applied to step 3.1 of the proposed algorithm. Using the Dijkstra's algorithm with minimal heaps, the time complexity is O(Q·(E + V·log(V))). When the Dijkstra's algorithm is implemented with linked list, the overall time cost is reduced to O(Q·(V + E)) [22]. In this last case, the algorithm is almost linear because Q is supposed to be small. Therefore, if the cardinality of Q is negligible against V + E, the cost is O(V + E).

Conclusions
Often the system under study is described by a system of DAEs. These equations are usually parameterized by the design variables, which are known during the simulation. If a new assumption equation is added to the system, it becomes over-constrained. There are several methods to reduce the new system to a well-constrained one, such as Hangos' method. This method basically consists of removing an equation from the original system to equal the number of equations and variables. However, this procedure requires that the new system contradicts the first one and, in some cases, the solution of the transformed system is far from one which describes the initial system. This paper presented an alternative to Hangos' method which preserves the original formulation of the system. The basis of this new method is to relax design variables, converting them into new algebraic variables. This procedure achieves a solution compatible with that of the initial system.
The example proposed in [21] has been studied with a simulation of a reactor with input and output flows. The differential system has been transformed into an integral one, as Hangos' method shows. Then, a new assumption equation has been added and a new system that is not related to the first one has been obtained. However, by applying the method suggested by the authors, an accurate system with the initial equations is obtained.
Also, an algorithm has been developed based on graph theory, starting from the initial matching of equation-variables. This algorithm has the followings characteristics: • It relates algebraic and state variables of the assumptions equations to the design variables of the system selected by the user.

•
It allows the addition of several assumption equations and obtaining a greater number of common edges with the starting system.

•
In the best case, depending on the implementation, the computational cost is O(Q(E + V)), where Q is the number of assumption equations, V the number of vertexes in the bipartite equation-variable graph and E is the number of edges in that graph. Therefore, if the number of assumption equations is small, the cost of the algorithm is linear.

•
The algorithm can be combined with an algorithm of index reduction.
Author Contributions: Antonio Belmonte and Juan Garrido conceived and design the method; Jorge E. Jiménez and Francisco Vázquez analyzed the method; and Francisco Vázquez and Juan Garrido wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.