Structural Analysis of Multimode DAE Systems: summary of results
II SS N - I S RN I NR I A / RR -- -- F R + E N G RESEARCHREPORTN° 9387
January 2021Project-Team Hycomes
Structural Analysis ofMultimode DAE Systems:summary of results
Albert Benveniste, Benoit Caillaud, Mathias Malandain a r X i v : . [ c s . P L ] J a n ESEARCH CENTRERENNES – BRETAGNE ATLANTIQUE
Campus universitaire de Beaulieu35042 Rennes Cedex
Structural Analysis of Multimode DAESystems: summary of results
Albert Benveniste ∗ , Benoit Caillaud, Mathias Malandain Project-Team HycomesResearch Report n ° This work was supported by the FUI ModeliScale DOS0066450/00 French national grant (2018-2021) and theInria IPL ModeliScale large scale initiative (2017-2021, https://team.inria.fr/modeliscale/ ). ∗ All authors are with Inria-Rennes, Campus de Beaulieu, 35042 Rennes cedex, France; email: [email protected] bstract:
Modern modeling languages for general physical systems, such as Modelica,Amesim, or Simscape, rely on Differential Algebraic Equations (DAEs), i.e., constraintsof the form f ( x (cid:48) , x, u ) = 0. This drastically facilitates modeling from first principles ofthe physics, as well as model reuse. In recent works [2, 3], we presented the mathematicaltheory needed to establish the development of compilers and tools for DAE-based physicalmodeling languages on solid mathematical grounds.At the core of this analysis sits the so-called structural analysis , whose purpose, atcompile time, is to either identify under- and overspecified subsystems (if any), or torewrite the model in a form amenable of existing DAE solvers, including the handling ofmode change events. The notion of “structure” collects, for each mode and mode changeevent, the variables and equations involved, as well as the latent equations (additionalequations redundant with the system), needed to prepare the code submitted to the solver.The notion of DAE index (the minimal number of times any equation has to be possiblydifferentiated) is part of this structural analysis.This report complements [2, 3] by collecting all the needed background on structuralanalysis. The body of knowledge on structural analysis is large and scattered, which alsomotivated us to collect it in a single report.We first explain the primary meaning of structural analysis of systems of equations,namely the study of their regularity or singularity in some generic sense. We then brieflyreview the body of graph theory used in this context. We develop some extensions, forwhich we are not aware of any reference, namely the structural analysis of systems ofequations with existential quantifiers.For the structural analysis of DAE systems, we focus on John Pryce’s Σ-method, thatwe both summarize and extend to non-square systems.The uses of these tools and methods in [2, 3] are highlighted in this report. Key-words: structural analysis, differential-algebraic equations (DAE), multi-modesystems, variable-structure models nalyse Structurelle des Syst`emes de DAE multimodes:recueil des r´esultats de base
R´esum´e :
Les langages modernes de mod´elisation de syst`emes physiques, tels que Modelica,Amesim ou Simscape, s’appuient sur des ´Equations Alg´ebro-Diff´erentielles (DAE, de l’anglais
Differential-Algebraic Equations ), c’est-`a-dire des contraintes de la forme f ( x (cid:48) , x, u ) = 0. Lamod´elisation `a partir des premiers principes de la physique, ainsi que la r´eutilisation de mod`eles,sont facilit´ees par cette approche. Dans des travaux r´ecents [2, 3], nous pr´esentons la th´eoriemath´ematique requise pour le d´eveloppement de compilateurs et d’outils pour les langages demod´elisation `a base de DAE sur des bases math´ematiques solides.Cette analyse s’appuie sur l’ analyse structurelle , dont le but, `a la compilation d’un mod`ele, estsoit d’identifier d’´eventuels sous-syst`emes sous- et sur-d´etermin´es, soit de r´e´ecrire le mod`ele en vuede son traitement par des solveurs de DAE existants, en prenant en compte les ´ev´enements dechangements de mode. La notion de ”structure” rassemble, pour chaque mode et chaque changementde mode, les ´equations et les variables impliqu´ees, ainsi que les ´equations latentes (des ´equationssuppl´ementaires redondantes), requises pour construire le code soumis au solveur num´erique. Lanotion d’ index d’un syst`eme de DAE, li´e aux nombres de diff´erentiations successives devant ˆetreappliqu´ees `a ses ´equations, est partie int´egrante de cette analyse structurelle.Le pr´esent rapport vient en compl´ement de [2, 3], en rassemblant toutes les connaissancesrequises en analyse structurelle. Ce travail a ´egalement ´et´e motiv´e par le fait que ce vaste corpus deconnaissances est ´eparpill´e dans la litt´erature.Nous expliquons d’abord le sens premier de l’analyse structurelle de syst`emes d’´equations, `asavoir, l’´etude de leur r´egularit´e ou singularit´e dans un sens g´en´erique. Nous passons ensuite enrevue les ´el´ements de th´eorie des graphes utilis´es dans ce contexte. Nous d´eveloppons ´egalement uneextension qui est, `a notre connaissance, in´edite dans la litt´erature sur le sujet : l’analyse structurellede syst`emes d’´equations contenant des quantificateurs existentiels.Pour l’analyse structurelle de syst`emes de DAE, nous nous focalisons sur la Σ-m´ethode de J.Pryce, qui est r´esum´ee puis ´etendue `a des syst`emes non-carr´es.Les utilisations de ces outils et m´ethodes dans [2, 3] sont mises en ´evidence dans ce rapport. Mots-cl´es : analyse structurelle, ´equations alg´ebro-diff´erentielles (DAE), syst`emes multi-mode,mod`eles `a structure variable tructural Analysis of Multimode DAE Systems: summary of results Contents
RT n ° A. Benveniste, B. Caillaud, M. Malandain
Modern modeling languages for general physical systems, such as Modelica, Amesim, or Simscape,rely on Differential Algebraic Equations (DAEs), i.e., constraints of the form f ( x (cid:48) , x, u ) = 0. Thisdrastically facilitates modeling from first principles of the physics, as well as model reuse. In recentworks [2, 3], we presented the mathematical theory needed to establish the development of compilersand tools for DAE-based physical modeling languages on solid mathematical grounds.At the core of this analysis sits the so-called structural analysis , whose purpose, at compile time,is to either identify under- and overspecified subsystems (if any), or to rewrite the model in a formamenable of existing DAE solvers, including the handling of mode change events. The notion of“structure” collects, for each mode and mode change event, the variables and equations involved, aswell as the latent equations (additional equations redundant with the system), needed to preparethe code submitted to the solver. The notion of DAE index (the minimal number of times anyequation has to be possibly differentiated) is part of this structural analysis. The body of knowledgeon structural analysis is large and scattered, which motivated us to collect it in a single report. Thisreport complements [2, 3] by collecting all the needed background on structural analysis. The usesof these tools and methods in [2, 3] are highlighted in this report. Structural analysis of systems of equations:
We first explain the primary meaning of the structuralanalysis of systems of equations, namely the study of their regularity or singularity in some genericsense. Structural information about a system of equations can be summed up by its incidencematrix; in particular, a square system yields a square matrix. If, after some permutation of row andcolumns, this matrix only has non-zero diagonal entries, then we say that this matrix, as well asthe underlying system of equations, are structurally regular: this notion is repeatedly used in theModelica community.However, the most commonly used structural representation of a system of equations is actuallythe bipartite graph characterizing the pattern of the incidence matrix (i.e., the ( i, j )-locations of itsnon-zero entries). Structural regularity or singularity, as well as derived notions, are then checkedon this bipartite graph. As a matter of fact, the systems of equations under study are often sparse(i.e., their incidence matrices have a small proportion of non-zero entries), and the association ofgraphs with sparse matrices is pervasive in linear algebra for high performance computing [12, 13].We briefly review the body of graph theory used in this context. We recall in particular theDulmage-Mendelsohn decomposition of a bipartite graph, which is the basis for structuring a general(non-square) matrix into its overdetermined, regular, and underdetermined blocks, structurally.This body of knowledge in linear algebra is referred to as structural analysis of matrices. Structuralanalysis extends to systems of numerical equations, by considering the Jacobian of the consideredsystem.The link between numerical regularity and structural regularity belongs to the folklore of linearalgebra; however, its understanding is critical for a correct use of structural analysis for systemsof numerical equations. We were, however, unable to find a reference where related results weredetailed; we propose such results here.We also develop some useful extensions, for which we are not aware of any reference, particularlythe structural analysis of systems of equations with existential quantifiers.
Structural analysis of DAE systems:
We focus on DAE and summarize John Pryce’s Σ-method forstructural analysis [17], and we present an extension of this method to non-square DAE systems.We also include the structural form of Campbell and Gear’s method of differential arrays [7], whichwe complement with its discrete-time counterpart, namely “difference arrays”. These methods relyon the structural analysis of systems of equations with existential quantifiers.
This section is a summary of Sections 1.2.1 and 2 of [3], to which the reader is referred for furtherdetails.
Inria tructural Analysis of Multimode DAE Systems: summary of results Our example is an idealized clutch involving two rotating shafts with no motor or brake connected(Figure 1). We assume this system to be closed, with no interaction other than explicitly specified.Figure 1: An ideal clutch with two shafts.The dynamics of each shaft i = 1 ,
2, is modeled by ω (cid:48) i = f i ( ω i , τ i ) for some, yet unspecified,function f i , where ω i is the angular velocity, τ i is the torque applied to shaft i , and ω (cid:48) i denotes thetime derivative of ω i . Depending on the value of the input Boolean variable γ , the clutch is eitherengaged ( γ = t , the constant “true”) or released ( γ = f , the constant “false”). When the clutch isreleased, the two shafts rotate freely: no torque is applied to them ( τ i = 0). When the clutch isengaged, it ensures a perfect join between the two shafts, forcing them to have the same angularvelocity ( ω − ω = 0) and opposite torques ( τ + τ = 0). Here is the model: ω (cid:48) = f ( ω , τ ) ( e ) ω (cid:48) = f ( ω , τ ) ( e ) if γ then ω − ω = 0 ( e ) and τ + τ = 0 ( e ) if not γ then τ = 0 ( e ) and τ = 0 ( e ) (1)When γ = t , equations ( e , e ) are active and equations ( e , e ) are disabled, and vice-versa when γ = f . If the clutch is initially released, then, at the instant of contact, the relative speed of the tworotating shafts jumps to zero; as a consequence, an impulse is expected on the torques. This model yields an ODE system when the clutch is released ( γ = f ), and a DAE system of index1 when the clutch is engaged ( γ = t ), due to equation ( e ) being active in this mode. The basicway of generating simulation code for the engaged mode ( γ = t ) is to add, to the model in thismode, the latent equation ( e (cid:48) ): ω (cid:48) = f ( ω , τ ) ( e ) ω (cid:48) = f ( ω , τ ) ( e ) ω − ω = 0 ( e ) ω (cid:48) − ω (cid:48) = 0 ( e (cid:48) ) τ + τ = 0 ( e ) (2)By doing so, we can apply the following basic policy:Regard System (2) as a system of equations in which the dependent variables arethe leading variables of System (2), i.e., the torques τ i and the accelerations ω (cid:48) i ,and consider the ω (cid:48) i ’s as dummy derivatives, i.e., fresh variables not related to thevelocities ω i ; ignore equation ( e ) and solve the four remaining equations for thedependent variables, using an algebraic equation solver. (3)This policy amounts to, first, bringing the DAE system back to an “ODE-like” form, then solving itas such. In this case, once the latent equation ( e (cid:48) ) was added, the Jacobian matrix of the system RT n ° A. Benveniste, B. Caillaud, M. Malandain ( e , e , e (cid:48) , e ) with respect to its leading variables was invertible. This is the classical approach toDAE simulation via “index reduction”—having differentiated ( e ) once indicates that System (2)has index 1. Basic tool 1
The current practice in DAE-based modeling languages is that adding latent equationsis performed on the basis of a structural analysis, which is a semi-decision procedure using graphtheoretic criteria, of much lower computational cost than that of studying the numerical invertibilityof the Jacobian. We provide the bases for structural analysis of DAEs in Sections and of thisreport. The analysis of the system above in both its modes is classical. The handling of mode changes,however, is not, due to the occurrence of impulsive torques in the “released → engaged” transition.Also, the nature of time differs between the study of modes and that of mode changes. Each mode,“released” or “engaged”, evolves in continuous-time according to a DAE system. In contrast, themode changes are events at which restart values for states must be computed from states beforethe change, in one or several discrete-time transitions—note that such restart transitions are notexplicitly specified in System (1). Let us focus on the “released → engaged” transition, where torques are impulsive.To unify the views of time, we map the different kinds of time to a same discrete-time line, bydiscretizing the real time line R ≥ using an infinitesimal time step ∂ . By doing so, the DAEs actingwithin the two “released” or “engaged” modes, get approximated with an infinitesimal error. By“infinitesimal”, we mean: “smaller in absolute value than any non-zero real number”. All of this canbe mathematically formalized by relying on Nonstandard Analysis [18, 9], see [2, 3] for details. Inparticular, System (1) is mapped to discrete-time by using the discrete-time line T = def { n∂ | n = 0 , , , . . . } (4)and replacing in it all the derivatives by their first order explicit Euler schemes, that is: ω (cid:48) i ( t ) is replaced by ω i ( t + ∂ ) − ω i ( t ) ∂ . (5)In (4), integer n must become large enough so that the whole time line R ≥ is covered—infiniteintegers must be used, which is also formalized in Nonstandard Analysis. Using the shift operators • x ( t ) = def x ( t − ∂ ) and x • ( t ) = def x ( t + ∂ ) , (6)we can rewrite System (1) as follows: ω • − ω ∂ = f ( ω , τ ) ( e ) ω • − ω ∂ = f ( ω , τ ) ( e ) if γ then ω − ω = 0 ( e ) and τ + τ = 0 ( e ) if not γ then τ = 0 ( e ) and τ = 0 ( e ) (7)Index reduction for the engaged mode can now be performed by shifting ( e ) and adding it to (7),shown in red: ω • − ω ∂ = f ( ω , τ ) ( e ) ω • − ω ∂ = f ( ω , τ ) ( e ) if γ then ω − ω = 0 ( e ) and ω • − ω • = 0 ( e • ) and τ + τ = 0 ( e ) if not γ then τ = 0 ( e ) and τ = 0 ( e ) (8) Inria tructural Analysis of Multimode DAE Systems: summary of results γ : f → t . At theconsidered instant, we have • γ = f and γ = t , see (6) for the definition of shift operators. UnfoldingSystem (8) at the two successive previous (shown in blue) and current (shown in black) instantsyields, by taking the actual values for the guard at those instants into account:previous ω − • ω ∂ = f ( • ω , • τ ) ( • e ∂ ) ω − • ω ∂ = f ( • ω , • τ ) ( • e ∂ ) • τ = 0 • τ = 0current ω • − ω ∂ = f ( ω , τ ) ω • − ω ∂ = f ( ω , τ ) ω − ω = 0 ( e ) ω • − ω • = 0 τ + τ = 0 (9)We regard System (9) as an algebraic system of equations with dependent variables being the leadingvariables of System (8) at the previous and current instants: • τ i , ω i ; τ i , ω • i for i = 1 ,
2. System (9) issingular since the following subsystem possesses five equations and only four dependent variables ω , ω , • τ , • τ : ω − • ω ∂ = f ( • ω , • τ ) ( • e ∂ ) ω − • ω ∂ = f ( • ω , • τ ) ( • e ∂ ) • τ = 0 • τ = 0 ω − ω = 0 ( e ) (10)We propose to resolve the conflict in (10) by applying the following causality principle: Principle 1 (causality)
What was done at the previous instant cannot be undone at the currentinstant.
This leads to removing, from subsystem (10), the conflicting equation ( e ), thus getting the followingproblem for computing restart values at mode change γ : f → t :solve for ω • , ω • , τ , τ the follow-ing system, and then, use ω • , ω • as restart values for the velocities: ω , ω , • τ , • τ set by previous instant ω • − ω ∂ = f ( ω , τ ) ω • − ω ∂ = f ( ω , τ ) ω • − ω • = 0 τ + τ = 0 (11)Note that the consistency equation ( e ) : ω − ω = 0 has been removed from System (11), thusmodifying the original model. However, this removal occurs only at mode change events γ : f → t .What we have done amounts to delaying by one infinitesimal time step the satisfaction of some ofthe constraints in force in the new mode γ = t . Since our time step is infinitesimal, this takes zeroreal time. The move from System (9) to System (11) must be automatized, hence we need: Basic tool 2
We need a graph based algorithm for identifying the conflicting equations possiblyoccurring at a mode change. This is addressed in Section . . of this report, where the Dulmage-Mendelsohn decomposition of a bipartite graph is presented. Now, System (11) is not effective yet, since it involves, in the denominator of the left-hand side ofthe first two equations, the infinitesimal number ∂ . Letting ∂ (cid:38) RT n ° A. Benveniste, B. Caillaud, M. Malandain
In this section, we develop the structural analysis for a system of algebraic equations, i.e., equationsin a Euclidian space, with no dynamics (no time, no derivative), of the form f i ( y , . . . , y k , x , . . . , x n ) = 0 , i = 1 , . . . , m . (12)This system is rewritten as F ( Y, X ) = 0, where Y and X denote the vectors ( y , . . . , y k ) and( x , . . . , x n ), respectively, and F is the vector ( f , . . . , f m ). This system has m equations, k freevariables (whose values are given) collected in vector Y , and n dependent variables (or unknowns)collected in vector X . Throughout this section , we assume that the f i ’s are all of class C at least.If the considered system (12) is square, i.e., if m = n , the Implicit Function Theorem (see, e.g.,Theorem 10 . . ν Y , ν X ) ∈ R k + n is a value for the pair ( Y, X ) such that F ( ν Y , ν X ) = 0 and the Jacobian matrix of F with respect to X evaluated at ( ν Y , ν X ) is nonsingular,then there exists, in an open neighborhood U of ν Y , a unique vector of functions G such that F ( v, G ( v )) = 0 for all v ∈ U . In words, Eq. (12) uniquely determines X as a function of Y , locallyaround ν Y .Denote by J X F the above-mentioned Jacobian matrix. Solving F = 0 for X , given a value ν Y for Y , requires forming J X F ( ν Y ) as well as inverting it. One could avoid inverting J X F (or evenconsidering it at all), by focusing instead on its structural nonsingularity, introduced in Section 3.1. Say that a subset of a Euclidian space is exceptional if it has zero Lebesgue measure. Say that aproperty P ( x , . . . , x k ) involving the real variables x , . . . , x k holds almost everywhere if it holds forevery x , . . . , x k outside an exceptional subset of R k .Square matrix P of size n is a permutation matrix if and only if there exists a permutation σ of the set { , . . . , n } such that p ij = 1 if j = σ ( i ) and p ij = 0 otherwise. Pre-multiplication(respectively, post-multiplication) of a matrix A by a permutation matrix results in permuting therows (resp., the columns) of A .These preliminary definitions and properties make it possible to prove the following lemma.This result, dealing in particular with the invertibility of sparse matrices, is part of the folklore ofHigh Performance Computing. However, we were unable to find any proper reference in which it isclearly stated (let alone proven). Lemma 1
Let A be an n × m -matrix with m ≥ n . The following two properties are equivalent:1. There exist two permutation matrices P and Q , such that PAQ = [ B B ] , where B is squarewith non-zero entries on its main diagonal—we say that A is structurally onto ;2. Matrix A remains almost everywhere onto (surjective) when its non-zero entries vary oversome neighborhood.Proof We consider the linear equation AX = Y , where Y ∈ R n has entries y , . . . , y n , and X isthe unknown vector with real entries x , . . . , x m . With a suitable renumbering of the coordinates of X , we can assume that Q is the identity matrix. For this proof, we will need to formalize what wemean by “letting the non-zero entries of matrix A vary over a neighborhood.” (13)To this end, we consider the non-zero pattern of matrix A , i.e., the subset P ⊆ [1 , n ] × [1 , m ]collecting the pairs ( i, j ) such that a ij is non-zero; let K be the cardinal of P . Order the set P bylexicographic order, namely ( i, j ) < ( i (cid:48) , j (cid:48) ) iff either i < i (cid:48) , or i = i (cid:48) ∧ j < j (cid:48) . Doing this defines anorder preserving bijection ψ : [1 , K ] → P . The non-zero entries of matrix A are then collected in Inria tructural Analysis of Multimode DAE Systems: summary of results A of R K whose components are the a ψ ( k ) , k = 1 , . . . , K . To every neighborhood U of A ,we can thus associate the set of ( n × m )-matrices V = (cid:8) A (cid:48) = ( a (cid:48) ij ) (cid:12)(cid:12) a (cid:48) ij = ξ ψ − ( i,j ) if ( i, j ) ∈ P , (cid:9) , (14)where the vector of R K , whose entries are the ξ k for 1 ≤ k ≤ K , ranges over U . A set V of matricesobtained in this way is called a zero-pattern preserving neighborhood of A , or simply neighborhoodof A when “zero-pattern preserving” is understood from the context. With this preliminary, wesuccessively prove the two implications. = ⇒ We need to prove the existence of a zero-pattern preserving neighborhood of A such that,for every Y , there exists an X such that A (cid:48) X = Y holds almost everywhere when matrix A (cid:48) rangesover this neighborhood. With a suitable renumbering of the equations, we can assume that P is theidentity matrix.Let V be a zero-pattern preserving neighborhood of A such that the first diagonal term a (cid:48) of A (cid:48) = PA (cid:48) remains non-zero when A (cid:48) varies over V . Since a (cid:48) (cid:54) = 0 holds when A (cid:48) varies over V ,one can perform, on A (cid:48) , the first step of a Gaussian elimination by using a (cid:48) as a pivot. As x is thus expressed in terms of x , . . . , x m , y , . . . , y n , the first equation may then be removed. Thisyields a reduced equation A (cid:48) X = Y , where A (cid:48) is an ( n − × ( m − Y ∈ R n − collects y σ (2) , . . . , y σ ( n ) , and X ∈ R m − collects the unknowns x , . . . , x m . Matrix A (cid:48) has diagonal entriesequal to a (cid:48) kk = a (cid:48) kk − a (cid:48) k /a (cid:48) , for k = 2 , . . . , n . Since a (cid:48) is non-zero, the set of entries a (cid:48) ij of matrix A (cid:48) causing a (cid:48) kk = 0 for some k is exceptional in R n × m since it requires the condition a (cid:48) k = a (cid:48) kk a (cid:48) tohold. The corresponding set of entries a (cid:48) ij of matrix A (cid:48) is exceptional as well; we denote by Ξ thisexceptional subset of R n × m . Thus, we assume that A (cid:48) stays within V \ Ξ .This ensures that the diagonal entries of matrix A (cid:48) are all non-zero. Thus, we can express x interms of x , . . . , x m , y σ (2) . We are then left with a further reduced equation A (cid:48) X (cid:48) = Y (cid:48) for whichthe same argument applies, namely: if matrix A (cid:48) stays within a certain zero-pattern preservingneighborhood V of A , but does not belong to some exceptional set Ξ , then the diagonal entries of A (cid:48) are all non-zero. And so on, by reducing the number of equations and unknowns.By iterating this process, we prove the existence of sets of matrices V , . . . , V n and Ξ , Ξ , . . . , Ξ n such that, for every matrix belonging to ( V ∩· · ·∩ V n ) \ (Ξ ∪ Ξ ∪· · ·∪ Ξ n ), the values of x , x , . . . , x n are uniquely determined as functions of the coordinates of Y and the other variables x n +1 , . . . , x m ,by pivoting. This concludes the proof of this implication, since V ∩ · · · ∩ V n is a zero-patternpreserving neighborhood of A and set Ξ ∪ Ξ ∪ · · · ∪ Ξ n is exceptional. = ⇒ We assume the existence of a zero-pattern preserving neighborhood V of A and anexceptional set Ξ such that A (cid:48) remains onto when varying over V \ Ξ. For each onto A (cid:48) , there existsat least one permutation matrix Q such that A (cid:48) Q = (cid:2) B B (cid:3) (15)where B is a square invertible matrix. Neighborhood V decomposes as the union V = (cid:83) σ V σ , where V σ collects the matrices A (cid:48) for which decomposition (15) holds with the matrix Q representingpermutation σ . Since the number of permutations is finite, at least one of the V σ has non-emptyinterior, say, for σ = σ ∗ . We can thus trim V to V σ ∗ , which still is a zero-pattern preservingneighborhood of A , and we call it again V for convenience.Hence, there exists a permutation matrix Q such that decomposition (15) holds for every A (cid:48) ∈ V . By Condition of Lemma 1, the pre-image of 0 by the map V (cid:51) A (cid:48) (cid:55)→ det( B ), wheredet( M ) denotes the determinant of a square matrix M , is an exceptional set. Now, we havedet( B ) = (cid:80) ni =1 ( − i +1 b i det( B i ), where B ij is the submatrix of B obtained by erasing row i and column j . Since B is almost everywhere invertible, there must be some i ∈ { , . . . , n } suchthat b i (cid:54) = 0 and B i is almost everywhere invertible—otherwise we would have det( B ) = 0 forany A (cid:48) ∈ W , where W is some neighborhood contained in V . Let P be the permutation matrixexchanging rows 1 and i , so we replace B by B = P B . Now, since we have b = b i (cid:54) = 0 and B is almost everywhere invertible. This process is iterated m − RT n ° A. Benveniste, B. Caillaud, M. Malandain its diagonal entry is non-zero, otherwise the determinant of the corresponding B matrix would bezero in a neighborhood of A . The wanted left permutation matrix is P = P m − . . .P P . (cid:3) Lemma 1 specializes to the following simpler result:
Corollary 2
Let A be an n × n -matrix. The following two properties are equivalent:1. There exist two permutation matrices P and Q , such that PAQ = B , where B has non-zeroentries on its main diagonal—we say that A is structurally nonsingular ;2. Matrix A remains almost everywhere invertible when its non-zero entries vary over someneighborhood. In this section, we focus on systems of algebraic equations F ( Y, X ) = 0 of the form (12), that is,equations involving no time and no dynamics. We will need to handle variables, their valuations,and vectors thereof. To this end, the following conventions will be used (unless no confusion occursfrom using more straightforward notations):
Notations 3
Lowercase letters ( x, y, ... ) denote scalar real variables; capitals (
X, Y, ... ) denotevectors of real variables; whenever convenient, we regard X as a set of variables and write x ∈ X torefer to an entry of X . We adopt similar conventions for functions ( f, g ,...) and vectors of functions( F, G ,...). A value for a variable x is generically denoted by ν x , and similarly for X and ν X . (cid:3) As explained in the introduction of this section, the existence and uniqueness of solutions forsystem (12) relates to the invertibility of its Jacobian J X F . The structural analysis of system (12)is built on top of the structural nonsingularity of its Jacobian, introduced in Section 3.1. Its fulldevelopment requires some background on graphs. This short background is compiled from [5, 12, 13]. Given a graph G = ( V, E ), where V and E ⊆ V × V are the sets of vertices and edges, a matching M of G is a set of edges of G such that notwo edges in M are incident on a common vertex. Matching M is called complete if it covers allthe vertices of G . An M -alternating path is a path of G whose edges are alternatively in M andnot in M —we simply say “alternating path” when M is understood from the context. A vertex is matched in M if it is an endpoint of an edge in M , and unmatched otherwise.A bipartite graph is a graph G = ( L ∪ R, E ), where E ⊆ L × R . Let F ( Y, X ) = 0 be a system ofequations of the form (12); its bipartite graph G F is the graph having F ∪ X as its set of vertices andan edge ( f i , x j ) if and only if variable x j occurs in function f i . We first consider square systems, in which the equations and the dependent variables are in equalnumbers.
Definition 4 (structural nonsingularity)
System F ( Y, X ) = 0 is called structurally nonsingu-lar if its bipartite graph G F possesses a complete matching. A structurally nonsingular system is necessarily square, i.e., with equations and dependent variables inequal numbers. Note that Condition 1 of Corollary 2 is a matrix reformulation of the graph theoreticcondition stated in Definition 4. Thus, F is structurally nonsingular in the sense of Definition 4, ifand only it its Jacobian J X F , evaluated at any pair ( ν Y , ν X ) such that F ( ν Y , ν X ) = 0, is structurally Inria tructural Analysis of Multimode DAE Systems: summary of results Lemma 5
Assume that F is of class C . The following properties are equivalent:1. System F ( Y, X ) = 0 is structurally nonsingular;2. For every ( ν Y , ν X ) satisfying F ( ν Y , ν X ) = 0 , the Jacobian matrix J X F ( ν Y , ν X ) remainsgenerically nonsingular when its non-zero coefficients vary over some neighborhood. Property 2 is interpreted using the argument developed in (13,14): every n × n -matrix having thesame non-zero pattern as J X F ( ν Y , ν X ) identifies with a unique vector of R K , where K is the cardinalof this non-zero pattern. We denote by J X F ∈ R K the image of the Jacobian J X F ( ν Y , ν X ) obtainedvia this correspondence. Property 2 says:There exist an open neighborhood U of J X F in R K , and a subset V ⊆ U such that: ( i ) the set U \ V has zero Lebesgue measure, and( ii ) every J ∈ V yields a regular matrix. (16)In the sequel, the so defined sets U and V will be denoted by U F ( ν Y , ν X ) and V F ( ν Y , ν X ) (17)or simply U F and V F when no confusion can result. (cid:3) Let F ( Y, X ) = 0 be structurally nonsingular and let M be a complete matching for it. Using M , graph G F can be directed as follows: edge ( f i , x j ) is directed from f i to x j if ( f i , x j ) ∈ M , from x j to f i otherwise. Denote by (cid:126) G M F the resulting directed graph. The following result holds ([11],Chapter 6.10): Lemma 6
The strongly connected components of (cid:126) G M F are independent of M . Each strongly connected component defines a block of F , consisting of the set of all equationsthat are vertices of this component. Blocks are partially ordered and we denote by (cid:22) F this order.Extending (cid:22) F to a total order (via topological sorting) yields an ordering of equations that putsthe Jacobian matrix J X F in Block Triangular Form (BTF).Figure 2: Illustrating Lemma 6 (we only show the projection of strongly connected components onfunction nodes).This lemma is illustrated in Figure 2, inspired by lecture notes by J-Y. L’Excellent and BoraU¸car, 2010. In this figure, a same bipartite graph G is shown twice (left and right), with the This result is implicitly used in [15, 14, 16, 17], as well as in a major part of the literature on the structuralanalysis of DAE systems. Generically means: outside a set of values of Lebesgue measure zero.RT n ° A. Benveniste, B. Caillaud, M. Malandain equations sitting on the left-hand side in green, and the variables on the right-hand side in red.Two complete matchings M (left) and M (right) are shown in thick blue, with other edges of thebipartite graph being black. The restriction, to the equation vertices, of the two directed graphs (cid:126) G M F (left) and (cid:126) G M F (right) are shown on both sides. Although the two directed graphs (cid:126) G M F and (cid:126) G M F differ, the resulting block structures (encircled in dashed brown) are identical. In our development, we will also encounter non-square systems of algebraic equations. For generalsystems (with a number of variables not necessarily equal to the number of equations, n (cid:54) = m in (12)),call block a pair β = ( f , x ), consisting of a subsystem f =0 of system F =0 and the subset x of itsdependent variables. Definition 7 (Dulmage-Mendelsohn)
For F = 0 a general system of algebraic equations, the Dulmage-Mendelsohn decomposition [16] of G F yields the partition of system F = 0 into the followingthree blocks, given some matching of maximal cardinality for G F : • block β O collects the variables and equations reachable via some alternating path from someunmatched equation; • block β U collects the variables and equations reachable via some alternating path from someunmatched variable; • block β E collects the remaining variables and equations.Blocks β O , β U , β E are the O verdetermined, U nderdetermined, and E nabled parts of F . Statement 1 of the following lemma ensures that Definition 7 is meaningful:
Lemma 8
1. The triple ( β U , β E , β O ) defined by the Dulmage-Mendelsohn decomposition does not depend onthe particular maximal matching used for its definition. We thus write DM( F ) = ( β U , β E , β O ) . (18)
2. System F is structurally nonsingular if and only if the overdetermined and underdeterminedblocks β O and β U are both empty. A refined DM decomposition exists, which in addition refines block β E into its indecomposable blocktriangular form. This is a direct consequence of applying Lemma 6 to β E in order to get its BTF.The following lemma is useful in the structural analysis of multimode DAE systems [2, 3]: Lemma 9
Let
DM( F ) = ( β U , β E , β O ) be the Dulmage-Mendelsohn decomposition of F = 0 . Then,no overdetermined block exists in the Dulmage-Mendelsohn decomposition of F \ β O . Comment 10
We could refine this result by only removing unmatched equations with referenceto the matching of maximal cardinality used for generating DM( F ). Then, denoting by F (cid:48) theremaining subsystem of F , we would still get that no overdetermined block exists in the Dulmage-Mendelsohn decomposition of F (cid:48) , and F (cid:48) ⊇ ( F \ β O ) (the inclusion being strict in general). However,this policy depends on the particular matching. In contrast, the rule defined by Lemma 9 isindependent from the matching. This Lemma is used in [2, 3] to generate restart code at modechanges, see Section 5.3.2 of this report; the independence with respect to the particular matchingconsidered was crucial to ensure the determinism of code generation. We distinguish between two different uses of structural analysis. The classical use assumes that wealready have a close guess of a solution, we call it the local use.
By contrast, the other use is called global.
We explain both here.
Inria tructural Analysis of Multimode DAE Systems: summary of results If a square matrix is structurally singular, then it is singular. The converse is false: structuralnonsingularity does not guarantee nonsingularity. Therefore, the practical use of Definition 4 is asfollows: We first check if F = 0 is structurally nonsingular. If not, then we abort searching fora solution. Otherwise, we then proceed to computing a solution, which may or may not succeeddepending on the actual numerical regularity of the Jacobian matrix J X F .Structural analysis relies on the Implicit Function Theorem for its justification, as stated at thebeginning of this section. Therefore, the classical use of structural analysis is local: let F ( Y, X ) = 0be a structurally nonsingular system of equations with dependent variables X , and let ν X satisfy F ( ν Y , ν X ) = 0 for a given value ν Y for Y ; then, the system remains (generically) nonsingular if ν varies in a neighborhood of ν Y . This is the situation encountered in the running of ODE or DAEsolvers, since a close guess of the solution is known from the previous time step. We are, however, also interested in invoking the structural nonsingularity of system F ( Y, X ) = 0when no close guess is known. This is the situation encountered when handling mode changes ofmultimode DAE systems [2, 3]. We thus need another argument to justify the use of structuralanalysis in this case.As a prerequisite, we recall basic definitions regarding smooth manifolds (see, e.g., [8], Section4.7). In the next lemma, we consider the system F ( Y, X ) = 0 for a fixed value of Y , thus Y isomitted. Lemma 11
Let k, m, n ∈ N be such that m < n and set p = n − m . For S ⊂ R n and x ∗ =( x ∗ , . . . , x ∗ n ) ∈ S , the following properties are equivalent:1. There exists an open neighborhood V of x ∗ and a C k -diffeomorphism F : V → W ⊂ R n suchthat F ( x ∗ ) = 0 and F ( S ∩ V ) is the intersection of W with the subspace of R n defined by the m equations w p +1 = 0 , . . . , w n = 0 , where w i denote the coordinates of points belonging to W .2. There exists an open neighborhood V of x ∗ , an open neighborhood U of in R p and ahomeomorphism ψ : U → S ∩ V , such that ψ , seen as a map with values in R n , is of class C k and has rank p at —i.e., ψ (cid:48) (0) has rank p . Statement 1 means that S is, in a neighborhood of x ∗ , the solution set of the system of equations f ( X ) = 0 , . . . , f m ( X ) = 0 in the n -tuple of dependent variables X , where F = ( f , . . . , f m ).Statement 2 expresses that ψ is a parameterization of S in a neighborhood of x ∗ , of class C k andrank p —meaning that S has dimension p in a neighborhood of x ∗ .In order to apply Lemma 11 to the non-local use of structural analysis, we have to studythe case of square systems of equations (i.e., let m = n ). Consider a system F ( X ) = 0, where X = ( x , . . . , x n ) is the n -tuple of dependent variables and F = ( f , . . . , f n ) is an n -tuple of functions R n → R of class C k . Let S ⊆ R n be the solution set of this system, that we assume is non-empty.To be able to apply Lemma 11, we augment X with one extra variable z , i.e., we set Z = def X ∪ { z } and we now regard our formerly square system F ( X ) = 0 as an augmented system G ( Z ) = 0 where G ( X, z ) = def F ( X ). Extended system G possesses n equations and n +1 dependent variables, hence, p = 1, and its solution set is equal to S × R . Let x ∗ ∈ S be a solution of F = 0. Then, we have G ( x ∗ , z ∗ ) = 0 for any z ∗ ∈ R . Assume further that the Jacobian matrix J X F ( x ) is nonsingular at x ∗ . Then, it remains nonsingular in a neighborhood of x ∗ . Hence, the Jacobian matrix J X G ( x, z )has rank n in a neighborhood V of ( x ∗ , z ∗ ) in R n +1 . We can thus extend G by adding one morefunction in such a way that the resulting ( n +1)-tuple ¯ G is a C k -diffeomorphism, from V to an openset W of R n +1 . The extended system ¯ G = 0 satisfies Property 1 of Lemma 11. By Property 2 ofLemma 11, ( S × R ) ∩ V has dimension 1, implying that S ∩ proj R n ( V ) has dimension <
1. Thisanalysis is summarized in the following lemma:
RT n ° A. Benveniste, B. Caillaud, M. Malandain
Lemma 12
Consider the square system F ( X ) = 0 , where X = ( x , . . . , x n ) is the n -tuple ofdependent variables and F = ( f , . . . , f n ) is an n -tuple of functions R n → R of class C k , k ≥ . Let S ⊆ R n be the solution set of this system. Assume that system F = 0 has solutions and let x ∗ besuch a solution. Assume further that the Jacobian matrix J X F ( x ) is nonsingular at x ∗ . Then, thereexists an open neighborhood V of x ∗ in R n , such that S ∩ V has dimension < . The implications for the non-local use of structural analysis are as follows. Consider the squaresystem F ( X ) = 0, where X = ( x , . . . , x n ) is the n -tuple of dependent variables and F = ( f , . . . , f n )is an n -tuple of functions R n → R of class C k . We first check the structural nonsingularity of F = 0.If structural nonsingularity does not hold, then we abort solving the system. Otherwise, we proceedto solving the system and the following cases can occur:1. the system possesses no solution;2. its solution set is nonempty, of dimension <
1; or3. its solution set is nonempty, of dimension ≥ F = 0, as its conclusion S ∩ V = { x ∗ } involvesan open neighborhood V of x ∗ . Of course, subclasses of systems for which existence and uniquenessare guaranteed are of interest. To support the structural analysis for the method of differential arrays advocated by Campbell andGear [7] (see also Section 4.3 of this report), we will need to develop a structural analysis for thefollowing class of (possibly non-square) algebraic systems of equations with existential quantifier: ∃ W : F ( X, W, Y ) = 0 , (19)where ( X, W ) collects the dependent variables of system F ( X, W, Y ) = 0. Y and X are as before,and the additional tuple W collects supplementary variables, of no interest to us, whence theirelimination by existential quantification. Our aim is to find structural conditions ensuring that (19)defines a partial function Y → X , meaning that the value of tuple X is uniquely defined by thesatisfaction of (19), given a value for tuple Y .At this point, a fundamental difficulty arises. Whereas (19) is well defined as an abstractrelation, it cannot in general be represented by a projected system of smooth algebraic equationsof the form G ( X, Y ) = 0. Such a G can be associated to (19) only in subclasses of systems. Ingeneral, no extension of the Implicit Function Theorem exists for systems of the form (19); thus,one cannot apply as such the arguments developed in Section 3.3.2. Therefore, we will reformulateour requirement differently. Say that ν Y is consistent if the system of equations F ( X, W, ν Y ) = 0possesses a solution for ( X, W ). (20)Consider the following property for the systen F ( X, W, Y ) = 0: for every consistent ν Y , F ( ν X , ν W , ν Y ) = 0 F ( ν X , ν W , ν Y ) = 0 (cid:41) = ⇒ ν X = ν X , (21) To our knowledge, the results of this section are not part of the folklore of sparse matrix algebra. Examples are linear systems for which elimination is easy, and polynomial systems for which it is doable—butexpensive—using Gr¨obner bases. Inria tructural Analysis of Multimode DAE Systems: summary of results X is independent of W given a consistent tuple of values for Y . To find structuralcriteria guaranteeing (21), we will consider the algebraic system F ( X, W, Y ) = 0 with
X, W, Y asdependent variables (i.e., values for the entries of Y are no longer seen as being given). Definition 13
System (19) is structurally nonsingular if the following two conditions hold, al-most everywhere when the non-zero coefficients of the Jacobian matrix J X,W,Y F vary over someneighborhood: • consistent values for Y exist; • condition (21) holds. The structural nonsingularity of (19) can be checked by using the DM decomposition of F asfollows. Let ( β U , β E , β O ) = DM( F ) be the DM decomposition of F ( X, W, Y ) = 0, with (
X, W, Y ) asdependent variables. We further assume that the regular block β E is expressed in its block triangularform, see Lemma 6 and comments thereafter. Let B be the set of all indecomposable blocks of β E and (cid:22) be the partial order on B following Lemma 6. The following holds: Lemma 14
System (19) is structurally nonsingular if and only if:1. Block β O is empty;2. Block β U involves no variable belonging to X ;3. For every β ∈ B containing some variable from X , then, β contains no variable from W ,and for every β (cid:48) ≺ β and every directed edge ( z (cid:48) , z ) ∈ (cid:126) G M F such that z (cid:48) ∈ β (cid:48) , z ∈ β , and M is anarbitrary complete matching for G F , then z (cid:48) (cid:54)∈ W . Conditions 1 and 2 speak by themselves. Regarding Condition 3, the intuition is the following.Every directed path of (cid:126) G M F , originating from W and terminating in X , of minimal length, musttraverse Y . Consequently, W influences X “through” Y only. Proof
We successively prove the if and then the only if parts. “If ” part:
By Condition 1, there exist consistent values for Y . By condition 2, no variable of X belongs to the underdetermined part of DM( F ). It remains to show that condition (21) holds. Bycondition (3), each indecomposable block β ∈ B involving a variable of X has the form G ( Z, U ) = 0,where1. Z is the n -tuple of dependent variables of G = ( g , . . . , g n ),2. no variable of W belongs to Z , and3. U ⊆ X ∪ Y is a subset of the dependent variables of the blocks β (cid:48) immediately preceding β .Since G = 0 is structurally regular, fixing a value for U entirely determines all the dependentvariables of block β . This proves the “if” part. “Only if ” part: We prove it by contradiction. If condition 1 does not hold, then the existenceof consistent values for Y is not structurally guaranteed. If condition 2 does not hold, then thevariables of X that belong to β U are not determined, even for given values of W, Y . If condition 3does not hold, then two cases can occur: • Some x ∈ X and w ∈ W are involved in a same indecomposable block β ∈ B . Then, condition(21) will not hold for x due to the structural coupling with w through block β . • There exists an indecomposable block β ∈ B of the form G ( Z, U ) = 0, where The precise meaning of this statement is the same as in Lemma 5, see the discussion therafter.RT n ° A. Benveniste, B. Caillaud, M. Malandain Z is the n -tuple of dependent variables of G = ( g , . . . , g n ), and Z contains a variable x ∈ X ,2. no variable of W belongs to Z , and3. some variable w ∈ W belongs to U , the set of the dependent variables of the blocks β (cid:48) immediately preceding β .Hence, w influences x , structurally, which again violates (21). This finishes the proof. (cid:3) Algorithm 1
ExistQuantifEqn
Require: F ; return ( b O , b U , F Σ , F Σ ) ( β U , β E , β O ) = DM( F ) if condition 1 of Lemma 14 holds then b O ← t ; if then b U ← t ; partition β E = F Σ ∪ F Σ else b U ← f else b O ← f We complement Lemma 14 with the algorithm
ExistQuantifEqn (Algorithm 1), which requires asystem F of the form (19). If condition 1 of Lemma 14 fails to be satisfied, then b O ← f is returned,indicating overdetermination. If conditions 2 or 3 of Lemma 14 fail to be satisfied, then b U ← f is returned, indicating underdetermination. Otherwise, ExistQuantifEqn succeeds and returns thevalue t for both Booleans, together with the decomposition F Σ ∪ F Σ of β E . In this decomposition: • Subsystem F Σ collects the indecomposable blocks involving variables belonging to X , so that F Σ determines X as a function of ν Y when ν Y is consistent; • Subsystem F Σ collects the consistency conditions, whose dependent variables belong to W ∪ Y .Our background on the structural analysis of algebraic equations is now complete. In the nextsection, we recall the background on the structural analysis of (single-mode) DAE systems. In this section, we consider DAE systems of the form f j ( x i ’s and derivatives) = 0 , (22)where x , . . . , x m denote the dependent signal variables (their valuations are trajectories) and f = 0 , . . . , f n = 0 denote the equations—with reference to notation (12), we omit the input signalscollected in Y . By analogy with DAEs, we use the acronym dAE to mean difference AlgebraicEquations, which define discrete-time dynamical systems of the form f j ( x i ’s and shifts) = 0 . (23)where x , . . . , x m denote the dependent stream variables (their valuations are data streams, whichcan be regarded as real sequences), f = 0 , . . . , f n = 0 denote the equations, and, for x = { x n | n = 1 , , . . . } a stream of variables, its successive shifts x • k , where k is a nonnegative integer, aredefined by x • kn = def x n + k ; we write for short x • = def x • . (24)Also, the notation x (cid:48) k is adopted throughout this paper, instead of the more classical x ( k ) , for the k -th derivative of x . Inria tructural Analysis of Multimode DAE Systems: summary of results Σ -method for square systems The DAE systems we consider are “square”, i.e., have the form (22), with m = n . Call leadingvariables of System (22) the d i -th derivatives x (cid:48) d i i for i = 1 , . . . , n , where d i is the maximaldifferentiation degree of variable x i throughout f = 0 , . . . , f n = 0. The problem addressed by thestructural analysis of DAE systems of the form (22) is the following. Regard (22) as a system ofalgebraic equations with the leading variables as unknowns. If this system is structurally nonsingular,then, given a value for all the x (cid:48) ki for i = 1 , . . . , n and k = 0 , . . . , d i −
1, a unique value for the leadingvariables can be computed, structurally; hence, System (22) is “like an ODE”. If this is not thecase, finding additional latent equations by differentiating suitably selected equations from (22) willbring the system to an ODE-like form, while not changing its set of solutions. Performing this isknown as index reduction.
Algorithms were proposed in the literature for doing it efficiently. Among them, the Pantelides’algorithm [15] is the historical solution. We decided, however, to base our subsequent developmentson the beautiful method proposed in [17] by J. Pryce, called the Σ-method. The Σ -method alsocovers the construction of the block triangular form and addresses numerical issues, which we donot discuss here.
Weighted bipartite graphs:
We consider System (22), which is entirely characterized by its setof dependent variables X (whose generic element is denoted by x ) and its set of equations F = 0(whose generic element is written f = 0). We attach to (22) the bipartite graph G = ( F ∪ X, E G )having an edge ( f, x ) ∈ E G if and only if x occurs in function f , regardless of its differentiationdegree. Recall that a matching is complete iff it involves all equations of F .So far, G is agnostic with respect to differentiations. To account for this, we further equip G with weights : to each edge ( f, x ) ∈ E G is associated a nonnegative integer d fx , equal to the maximaldifferentiation degree of variable x in function f . This yields a weight for any matching M of G bythe formula w ( M ) = (cid:80) ( f,x ) ∈M d fx . Suppose we have a solution to the following problem: Problem 1
Find a complete matching M for G and nonnegative integer offsets { c f | f ∈ F } and { d x | x ∈ X } , satisfying the following conditions, where the d fx are the weights as before: d x − c f ≥ d fx for all ( f, x ) ∈ E G , with equality if ( f, x ) ∈M c f ≥ for all f ∈ F . (25)Then, differentiating c f times each function f yields a DAE system F Σ = 0 having the followingproperties. Inequality d x ≥ c f + d fx holds for each x ∈ X , and d x = c f + d fx holds for the unique f such that ( f, x ) ∈M . Hence, the leading variables of DAE system F Σ = 0 are the d x -th derivatives x (cid:48) d x . Consequently, system F Σ = 0, now seen as a system of algebraic equations having x (cid:48) d x asdependent variables, is structurally nonsingular by Definition 4. Hence, F Σ = 0 is “like an ODE”.The integer k = max f ∈ F c f is called the index of the system. Definition 15
For F a DAE system, the solution to Problem yields the DAE system F Σ = { f (cid:48) c f | f ∈ F } together with the system of consistency constraints F Σ = { f (cid:48) k | f ∈ F, ≤ k
Lemma 16
Let G be a given bipartite graph and let two families of weights ( d fx ) ( f,x ) ∈ E G and ( d fx ) ( f,x ) ∈ E G be related by d fx = M × d fx for every ( f, x ) ∈ E G , where M is a fixed positive integer.Then, the offsets of the corresponding Σ -method are also related in the same way: d x = M × d x forevery variable x , and c f = M × c f for every function f . Inria tructural Analysis of Multimode DAE Systems: summary of results Algorithm 2
FindOffsets ( G is a bipartite graph with weights { d fx | ( f, x ) ∈ E G } )1. Solve LP (27), which gives a complete matching M of maximum weight for G ; any methodfor solving LP can be used.2. Apply the following iteration until a fixpoint is reached (in finitely many steps), from theinitial values c f = 0:(a) ∀ x : d x ← max { d fx + c f | ( f, x ) ∈ E G } ;(b) ∀ f : c f ← d x − d fx where ( f, x ) ∈ M .A necessary and sufficient condition for Problem 1 to have a solution is that the set of completematchings for G is non-empty. This criterion can be evaluated prior to computing the optimaloffsets. This completes our background material on structural analysis. Σ -method for non-square systems Structural analysis must be extended to non-square systems if we wish to allow for guards be-ing evaluated at the current instant—in this case, enabled/disabled equations are progressivelydetermined along with the progressive evaluation of guards, see Appendix A of [3]. Structuralanalysis is also a powerful tool for addressing problems other than just the simulation of DAEsystems, e.g., the generation of failure indicators for system diagnosis. For both uses, we will needto address other cases than square DAE systems. We develop here an adaptation of the Σ-methodto non-square DAE systems. In our development, we could simply copy the original constructions,lemmas, and proofs, of [17], while performing marginal adjustments to handle non-square systems.This is presented next.If F = 0 is a non-square DAE system, then its weighted bipartite graph G has equation nodesand variable nodes in different numbers, hence the concept of complete matching does not apply. Inthis case, we consider instead the following weaker notion:we say that a matching M is equation-complete if it covers all the equation nodes. (30)No equation-complete matching exists for an overconstrained system F = 0, i.e., a system with lessvariables than equations. This being said, we still consider the primal/dual problems (27) and (28).Again, if G possesses an equation-complete matching, then problem (27) has a solution, since wecan complete this matching to find a feasible solution. Hence, the existence of at least one matchingis still used, except that this matching has to be equation-complete instead of complete.Next, consider respective optima of (27) and (28). They still satisfy the slackness conditions (29),which yields the following results. • An optimal solution of (27) defines a subgraph
H ⊆ G by: ( f, x ) ∈ H if and only if ξ fx = 1.Note that H is not a matching, as one equation may be associated to more than one variablein this graph. In contrast, each variable is associated to exactly one equation. • By the slackness conditions (29), for each equation f = 0 and each x such that ( f, x ) ∈ H , x occurs in f with differentiation degree equal to the maximum differentiation degree in theentire system, namely d x .Call F Σ the DAE system defined by H and denote by Z = { x d x | x ∈ X } the set of its leadingvariables. F Σ is then regarded as an algebraic system of equations with Z as dependent variables,and we apply the Dulmage-Mendelsohn (DM) decomposition to it. This yields ( β U , β E , ∅ ), as weknow that the overconstrained block β O is empty. RT n ° A. Benveniste, B. Caillaud, M. Malandain
To justify our above adaptation of Pryce’s Σ-method, the following questions must be answered:
Question 1
How to adapt Step 2 of Algorithm 2, that is, the fixpoint iteration used for computingthe offsets?
Question 2
Is existence and uniqueness of the smallest solution still guaranteed in the non-squarecase?
Addressing Question : Step 2a of Algorithm 2 is not modified. Step 2b, on the other hand, cannotbe kept as is, since H is no longer a matching and there may be more than one x such that ( f, x ) ∈H .Step 2 of Algorithm 2 is thus modified as follows:(a) ∀ x : d x ← max { d fx + c f | ( f, x ) ∈ E G } ;(b) ∀ f : c f ← max { d x − d fx | ( f, x ) ∈ H} .We call Algorithm 2 (cid:48) the resulting algorithm. Addressing Question : This question is answered by the following lemma:
Lemma 17 (non-square Σ -method) Algorithm 2 (cid:48) converges to a fixpoint if and only if H is anoptimal solution of the primal problem (27) . Let c ∗ , d ∗ be this fixpoint when it exists. Then c ∗ , d ∗ isan integer-valued optimal solution for the dual problem (28) and, for any optimal solution c † , d † forthis dual problem, c † ≥ c ∗ and d † ≥ d ∗ hold. The last statement expresses the uniqueness of the smallest optimal solution. Lemma 17 was statedand proved in [17], albeit for square DAE systems. Here, we will be carefully following the steps ofthat proof, while checking that they suit the non-square case.
Proof
Following [17], let φ : c (cid:55)→ φ ( c )be the mapping corresponding to the successive application of steps (a) and (b) to c .Assume that H is an optimal solution of (27) and c , d is an optimal solution of its dualproblem (28). By slackness conditions (29), one has d x = c f + d fx for ( f, x ) ∈ H . Starting from c ,step (a) yields d ( c ) = (cid:8) max ( f,x ) ∈ E G ( d fx + c f ) (cid:12)(cid:12) x ∈ X (cid:9) = { d fx + c f where ( f, x ) ∈ H | x ∈ X } so that, in turn, applying step (b) yields c : hence, φ ( c ) = c (actually, in step (b), all the d x − d fx for ( f, x ) ∈ H are equal). In other words, c is a fixpoint of φ .Conversely, let H ⊆ G be any feasible solution of (27), not necessarily optimal, and let c ∗ , d ∗ beany fixpoint for φ (note that φ depends on H ). We have d ∗ x ≥ d fx + c ∗ f for any ( f, x ) ∈ E G , withequality if ( f, x ) ∈ H . Assume that c ∗ ≥ c ∗ , d ∗ is a feasible solutionof the dual problem (28). Finally, H is an optimal solution of the primal problem, and c ∗ , d ∗ anoptimal solution of the dual problem, because they satisfy the slackness conditions (29), which areknown to characterize optimal solutions. By contraposition, if H is not optimal, then no fixpoint of φ exists.It remains to be proved that c ∗ ≥ φ is nondecreasing: c (cid:48) ≤ c implies φ ( c (cid:48) ) ≤ φ ( c ). Denote by c , d , c , d . . . the alternating sequence of c ’s and d ’s obtained by loopingover steps (a) and (b), alternatively. We have c k = φ ( c k − ) for k >
0. If c = 0, then d ≥ c f = d x − d fx ≥
0, where x is any variable such that ( f, x ) ∈ H (the value d x − d fx is independentfrom this choice). As a consequence, c ≥
0. Since φ is nondecreasing, applying φ to both sidesof c ≥ c = 0 yields c ≥ c , and so on, so that the algorithm yields a nondecreasing sequence,showing that c ∗ ≥ c † , d † is another dual-optimal solution for the offsets, then c † is a fixpoint of φ andthus c † = φ ( c † ) ≥ φ (0), hence c † ≥ φ k (0) for every k , whence c † ≥ c ∗ follows. (cid:3) Inria tructural Analysis of Multimode DAE Systems: summary of results In this section, we closely follow Section 5.3 of [17], regarding the comparison of the Σ-method withPantelides’ method [15] for the square case, and adapt it to the non-square case.Pantelides considers DAE systems with degree at most 1 for all variables. His construction findsa minimally structurally singular (MSS) subset of equations and differentiates each equation in thisset, creating an augmented system. This is repeated until no more MSS are found. Pantelides writesthe system as n + m equations F ( X, Z ) = 0, where X stands for the n variables whose derivativesappear, Y stands for the m algebraic variables, and Z = ( X (cid:48) , Y ). Pantelides’ construction appliesto this reformulation. For our comparison, the Z coincides with the collection of leading variables x d x where d x = max f d fx is the leading differentiation degree of variable x in the system.We assume that the primal AP (27) has an optimal solution H ∗ , and that there exist smallestoptimal offsets c ∗ f , d ∗ x satisfying (cid:80) x d ∗ x − (cid:80) f c ∗ f = w ( H ∗ ) (the total weight of H ∗ , defined as (cid:80) ( f,x ) ∈H ∗ d fx ). Let (cid:98) d x = def max f d fx (31)and define the leading derivative pattern L of the weighted bipartite graph G (with weights d fx ) asthe set L = { ( f, x ) | d fx = (cid:98) d x } For some equation f , respectively some subset E of equations, consider L f = { x | ( f, x ) ∈ L } , and L ( E ) = (cid:91) f ∈ E L f . Following Pantelides, we say that E is structurally singular (SS) if | L ( E ) | < | E | . E is minimally SS(MSS) if it is SS and has no proper SS subset. The link with the Σ-method is established in thefollowing results. Lemma 18
1. If E is MSS with k elements, then | L ( E ) | = k −
2. The following three statements are equivalent:(a) L has an equation-complete matching, see (30) ;(b) There are no SS subsets;(c) All c i are equal to .Proof For statement 1, if L ( E ) has fewer than k − E yields a set E (cid:48) satisfying | L ( E (cid:48) ) | ≤ | L ( E ) | < k − | E (cid:48) | . Therefore, E (cid:48) is SS, which contradictsthe minimality of E .The equivalence 2a ⇔
2b is just Hall’s Theorem applied to the sets L e .2a ⇒ L has an equation-complete matching H , then (by definition of L ) for ( f, x ) ∈H , d fx = (cid:98) d x . In other words, d fx is maximal among all the d (cid:98) fx for (cid:98) f ranging over the set of all equations.This means that H is optimal for the primal problem. We have d ∗ x − c ∗ f = d fx = (cid:98) d x for ( f, x ) ∈ H hence c ∗ f = 0 , d ∗ x = (cid:98) d x is the (unique) smallest optimal solution for the dual problem. ¬ ⇒ ¬ H ∗ for the primal problem contains a( f, x ) with d fx < (cid:98) d x . Then, by the slackness condition (29): d ∗ x = d fx + c ∗ f < (cid:98) d x + c ∗ f . (32) RT n ° A. Benveniste, B. Caillaud, M. Malandain
On the other hand, since d ∗ y , c ∗ f is a dual-optimal solution, we have d ∗ y ≥ d fy + c ∗ f ≥ d fy for all y ,whence d ∗ y ≥ (cid:98) d y . (33)In particular, d ∗ y ≥ (cid:98) d y ; combining this result with (32) yields c f > d ∗ x − (cid:98) d x ≥
0, so that 2c is false. (cid:3)
Lemma 19
Let H ∗ be a primal-optimal solution as before. If c ∗ f o = 0 is satisfied for some equation f o , then1. Let x o be such that ( f o , x o ) ∈H ∗ , then ( f o , x o ) ∈ L ;2. c f = 0 holds for any f (cid:54) = f o such that ( f, x o ) ∈ L .Proof We have d f o x o = d ∗ x o (since d ∗ is optimal and c ∗ f o = 0) ≥ (cid:98) d x o (by (33)) ≥ d f o x o (by (31))hence we get d ∗ x o = (cid:98) d x o (34)and d f o x o = (cid:98) d x o , which proves statement 1. For statement 2, if ( f, x o ) ∈ L for some f , then0 = (cid:98) d x o − d fx o = d ∗ x o − d fx o ≥ c ∗ f ≥ c ∗ f = 0. (cid:3) Lemma 20
For E any MSS, we have ∀ f ∈ E : c f > .Proof Suppose the lemma is false, meaning that E = def { f ∈ E | c f = 0 } is nonempty. Let | E | = k and | E | = l , so that 1 ≤ l ≤ k . By statement 1 of Lemma 19, we have L ( E ) ≥ l . Since | L ( E ) | < | E | by the assumption that E is an MSS, we deduce l < k , hence E \ E is nonempty. Pick any f ∈ E \ E , then c f >
0. By statement 2 of Lemma 19, L f cannot contain any x that is reachablefrom E via H ∗ : L ( E \ E ) ∩ [ ∃ f. (( E × X ) ∩ H ∗ )] = ∅ (35)The set on the right-hand side of ∩ is contained in L ( E ) (by statement 1 of Lemma 19) and hascardinality at least equal to | E | . By (35), we get | L ( E \ E ) | ≤ | L ( E ) | − | E | < | E \ E | , hence E \ E is a SS subset strictly contained in E , thus contradicting the minimality of E . (cid:3) Following these results, Pantelides’ algorithm can be extended to non-square systems in prettymuch the same way the Σ-method was extended. This non-square Pantelides’ algorithm locatesMSS subsets and differentiates them, until there is no such subset left. Note that, when Pantelides’algorithm locates a MSS subset I and differentiates it, this amounts to modifying the offsets foundby the Σ-method in the following way: • the d j remain unchanged • the c i decrease by 1 if i ∈ I , remain unchanged otherwise.From this and the preceding results, one gets: Theorem 21 (equivalence with Pantelides)
Pantelides’ algorithm gives the same results asthe Σ -method, in that the same offsets c i are found. Inria tructural Analysis of Multimode DAE Systems: summary of results Differential arrays were proposed in [7] as a tool for finding the latent equations of a DAE system.To F = 0 a DAE system, we associate its k -th differential array A k defined, for any k ∈ N ≥ , by A k = def F ddt F d dt F ... d k dt k F (36)where ddt denotes the total time derivative. Thus, for k sufficiently large, array system A k = 0should contain all the latent equations that must be added to F = 0 in order to reduce its index to0. More precisely, the index reduced counterpart of F = 0 is the following system with existentialquantifiers: ∃ W : A k ( X, W, Y ) = 0 (37)where: • Y collects the free variables of system F = 0; • X collects the dependent variables of system F = 0; • W collects other supplementary variables involved in array A k . Difference arrays are defined similarly for discrete-time dAE systems, by replacing the j -th timederivatives d j dt j F by the j -th forward shift F • j defined in (24), when building array (36): A k = def FF • F • ... F • k (38)The structural analysis of systems with existential quantifiers, of the form (19), was studied inSection 3.4. We can apply it to the structural analysis of arrays (36) or (38) to find the latentequations of a DAE system, and, thus, its index.The resulting algorithm is much less efficient than Pryce’s Σ-method or Pantelides algorithm.The reason is that, unlike the above methods, the one of Section 3.4 does not exploit the fact thatthe rows of the array are successive shifts of the same dAE system. In turn, the method of Section 3.4extends to time-varying discrete-time systems. We just need to replace, in (38), the shifted versionsof F by the successive F ( t ) , . . . , F ( t k ) at successive discrete instants t , . . . , t k . This is used in [2, 3]to develop a structural analysis of finite cascades of mode changes. In this section, we briefly review the results of [2, 3] regarding multimode DAE systems simulation.We indicate where, in [2, 3], the developments of Sections 3 and 4 of this report are used.
RT n ° A. Benveniste, B. Caillaud, M. Malandain
In [2, 3], a mathematical definition for multimode DAE systems is provided, and we also definetheir discrete-time counterpart, namely multimode dAE systems, of the respective formsmultimode DAE: if γ j ( x i ’s and derivatives) then f j ( x i ’s and derivatives) = 0multimode dAE: if γ j ( x i ’s and shifts) then f j ( x i ’s and shifts) = 0 , (39)where x • k , the k -shift of x , is defined in (24). To allow for a uniform handling of modes and mode change events in the structural analysisof multimode DAE, we use a nonstandard semantics for DAEs and multimode DAEs. In thisnonstandard semantics,derivative x (cid:48) is interpreted as its first-order forward Euler scheme x • − x∂ , (40)where continuous-time is discretized by using an infinitesimal positive time step ∂ and the forwardshift x • is defined in this discrete-time. Here, as already stated at the beginning of Section 2.3,“infinitesimal” means “smaller than any positive real number”, whichcan be given a formal meaning in nonstandard analysis [9, 18].This mapping allows for a discretization of multimode DAE systems with an infinitesimal error, whichyields a faithful discrete-time approximation. Of course, there is no free lunch: the nonstandardsemantics is not effective in that there is no computer that can execute it. Still, it is highly usefulfor the structural analysis, which is a symbolic analysis. The nonstandard semantics of a multimodeDAE systems is a multimode dAE system, in which • occurrences of increments x • − x∂ captures the continuous-time dynamics of the considered DAE(for this reason, we call continuous the modes in which such dynamics occurs), whereas • mode changes are represented by transitions relating the value x − of a state just before thechange and its value x + right after the change, encoded in the nonstandard semantics as x and x • . The structural analyses of both DAE and dAE systems coincide, by replacing the differentiationoperator x (cid:55)→ x (cid:48) arising in DAEs, by the forward shift operator x (cid:55)→ x • arising in dAEs. Pryce’sΣ-method applies to both DAE and dAE; for the latter, the equation offsets c f indicate how manytimes each equation f must be shifted. Since structural analyses mirror each other in continuous and discrete-times through the mapping(40), the structural analysis of continuous modes is performed as usual, with the help of the Σ-methodpresented in Section 4. This is summarized in the following
Basic tool 3
Standard structural analysis can be performed for each continuous mode, in order toadd the needed latent equations. The continuous-time structural analysis (where latent equations arefound by differentiation) is therefore used to generate code within continuous modes.
Inria tructural Analysis of Multimode DAE Systems: summary of results In order to prepare for the generation of restart code at mode changes, the discrete-time structural analysis (where latent equations are found by shifting) based on the nonstandardsemantics is applied to the continuous modes before and after the change.
Basic tools 3 and 4 rely on J. Pryce’s Σ-method. Once latent equations have been added foreach continuous mode, we may have a conflict, at mode change events separating two successivecontinuous modes, between: • the next values of states predicted by the dynamics in the previous mode; and • the consistency equations, generated by the structural analysis of the new mode.Is this pathological? Not quite. Such a situation can occur even for physically meaningful models.If the new continuous mode has positive index, then, nontrivial consistency equations are generatedby the structural analysis, which may be conflicting with the dynamics of the continuous modebefore the change. Since one cannot regard this model as being incorrect, we need a policy forhandling these conflicts. Handling isolated mode changes separating two successive continuous modes:
In [2, 3],it is proposed to call causality to the rescue: for isolated mode changes, we give priority to theprevious dynamics and postpone for a while the conflicting consistency equations generated by thenew continuous mode. Here, “for a while” means “for a finite number of infinitesimal time steps”,which amounts to zero time in real duration.
Basic tool 5
Identifying the conflicting consistency equations at mode changes is performed byusing the Dulmage-Mendelsohn decomposition (Definition and associated Lemmas and ). Conflicting consistency equations of the new continuous mode are postponed for a finite number ofnonstandard instants. The resulting discrete-time system operates for finitely many nonstandardinstants, and is the basis for recovering the restart values of the states in the new mode. Thisestablishes the structural analysis of mode changes separating two successive continuous modes.
Handling finite cascades of mode changes:
We also handle in [2, 3] finite cascades of successivemode changes separating two successive continuous modes. Since the dynamics varies over thesuccessive instants of the cascade, time-invariance no longer holds, and, thus, the Σ-method doesnot apply. However, the method of difference arrays, closely derived from Campbell and Gear’s differential arrays [7], can be adapted to time-varying discrete-time systems, and invoked. SeeSection 4.3 and formula (38) for a short introduction to difference arrays.For time-invariant dynamics, difference arrays are obtained by stacking, one below the other,the dAE dynamics and its successive shifted versions. For a finite cascade of mode changes leadingto a new continuous mode, however, we form the associated difference array by stacking, one belowthe other, the dynamics that holds at each successive instant of the cascade. More precisely, let t bethe current instant in the nonstandard semantics, and let ∂ be the infinitesimal time step. Then,the time-varying array associated to the cascade of mode changes is A k ( t ) = def F ( t ) F ( t + ∂ ) F ( t + 2 ∂ )... F ( t + k∂ ) (41)where F ( t ) is the dynamics at the first event of the cascade, F ( t + ∂ ) is the dynamics at the secondevent of the cascade, and so on until the instant t + k∂ , which sits within the new continuous mode.Enough rows should be added to the array, so that the leading variables of the first event of the RT n ° A. Benveniste, B. Caillaud, M. Malandain cascade are uniquely determined as functions of the states before the change. As for differentialarrays, the extra variables brought by the successive shifting must be eliminated by existentialquantification—see the comments regarding formula (11) introducing differential arrays in [7]. Finitecascades of mode changes are thus handled by using the results of Section 3.4 regarding equationswith existential quantifiers:
Basic tool 6
Finite cascades of mode changes are handled by forming the difference array (41) and analyzing it by using the method of Section . regarding equations with existential quantifiers. Fixpoint between guards and the equations they control:
In a multimode DAE system, alogico-numerical fixpoint occurs when a variable x is determined by a different dynamics dependingon the value of a guard that itself depends, directly or indirectly, on x . In [2, 3], only multimode DAEsystems having no logico-numerical fixpoint are supported. In such models, the above guard shouldrather depend on the left-limit x − of x , where x − ( t ) = lim s (cid:37) t x ( s ). Nevertheless, in Appendix Aof [3], we propose a method for handling multimode DAE systems with fixpoints but no cascades oflength >
1. This method uses the following basic tool:
Basic tool 7
The Σ -method for non-square systems, presented in Section 4.2. For S a multimode DAE model having continuous modes separated by finite cascades of modechange events, the method proposed in [3] returns the following: Case of a structurally incorrect model:
The compilation algorithm returns, for each modeand mode change event in which structural analysis fails, the overdetermined subsets of equationsand the underdetermined subsets of variables.
Case of a structurally correct model, producing an mDAE interpreter : The compilationalgorithm returns an mDAE interpreter, which is a labeled bipartite graph G whose vertices are: • state and leading variables x of the system, as well as guards γ ; • structurally nonsingular blocks β of algebraic equations determining leading variables fromstate variables, and expressions exp determining the values of guards as function of state andleading variables.A branch x → β exists in G if x is a free variable (input) of block β and a branch β → x exists in G if x is a dependent variable (output) of block β . A branch x → exp exists in G if exp has x as oneof its arguments, and a branch exp → γ exists in G if γ is computed using expression exp .Not every block is involved in a given mode or mode change. Thus, to each block β we associatea label λ ( β ), which is a property characterizing the set of all modes or mode changes in which β is active. A label is then assigned to every path π traversing blocks β , . . . , β k , . . . , β K by setting λ ( π ) = (cid:86) K λ ( β k ). Every circuit of G has label f , which ensures that the code is not globally circular.Also, for any two blocks β (cid:54) = β leading to the same variable x in G , we have λ ( β ) ∧ λ ( β ) = f ,which ensures that two different blocks never compete at determining the same variable x . Thisinterpreter encodes all the different schedulings that are valid for the different modes of the system.Note that modes are not enumerated, which is essential in order for this technique to scale up.This technique is reminiscent of the conditional dependency graphs used in the compilation ofthe Signal synchronous language [1, 4].The notion of mDAE interpreter is illustrated in Figure 5, for an RLDC2 circuit example withschematic shown in Figure 3 and model given in System (42), see Section 12 of [3] for details. Thisform of mdAE interpreter is used in the IsamDAE tool [6], by representing all the label usingBDD-related data structures. BDD means: Binary Decision Diagram, a data structure developed for model checking. Inria tructural Analysis of Multimode DAE Systems: summary of results i + i + j + j ( K ) x + w = u + v ( K ) u + v = u + v ( K ) u + v = x + w ( K ) w = L · j (cid:48) ( L ) w = L · j (cid:48) ( L ) i = C · v (cid:48) ( C ) i = C · v (cid:48) ( C ) x = R · j ( R ) x = R · j ( R ) s = if γ then i else − u ( S ) s = if γ then i else − u ( S )0 = if γ then u else i ( Z )0 = if γ then u else i ( Z ) γ = ( s − ≥
0) ( γ − ) γ = ( s − ≥
0) ( γ − ) (42)Figure 4: The RLDC2 model ( n.b. : variable y (cid:48) denotes the time derivative of y ). This report complements references [2, 3] by collecting and providing all missing details related tostructural analysis, for both algebraic and DAE systems, in a systematic way.
RT n ° A. Benveniste, B. Caillaud, M. Malandain g1 & !g2 :Z2 ->i2 !g1 | !g2 :i2 ->C2 ->v2'g1 & !g2g1 & !g2 :i2 j1 j2 ->K1 ->i1g1 & !g2!g1 & g2 :Z2 ->u2 g1 | g2 :u2 v2 x2 ->K4 ->w2!g1 & g2!g1 & g2 :u2 v1 v2 ->K3 ->u1!g1 & g2!g1 & !g2 :Z2' ->i2' !g1 & !g2 :i1' i2' v1 v2 x1 x2 -> K1' K2 K3 K4 L1 L2 ->w2 j1' u2 u1 w1 j2'!g1 & !g2g1 & g2 :Z2' ->u2' g1 & g2 :j1 j2 u1' u2' ->K1 K3' C1 C2 ->v2' i1 v1' i2g1 & g2 g2 :i2 ->S2 ->s2!g2 :u2 ->S2 ->s2!g1 & g2 :Z1 ->i1 !g1 | !g2 :i1 ->C1 ->v1'!g1 & g2!g1 & g2 :i1 j1 j2 ->K1 ->i2!g1 & g2g1 & !g2 :Z1 ->u1 g1 & !g2 :u1 v1 v2 ->K3 ->u2g1 & !g2 g1 | g2 :u1 v1 x1 ->K2 ->w1g1 & !g2!g1 & !g2 :Z1' ->i1' !g1 & !g2g1 & g2 :Z1' ->u1' g1 & g2 g1 :i1 ->S1 ->s1!g1 :u1 ->S1 ->s1 g1 | g2 :w2 ->L2 ->j2'g1 | g2 :w1 ->L1 ->j1'j2 ->R2 ->x2 g1 | g2!g1 & !g2j1 ->R1 ->x1 g1 | g2!g1 & !g2 g1 | g2!g1 & g2!g1 & g2g1 & !g2g1 & !g2 g1 | g2g1 & !g2g1 & !g2!g1 & g2!g1 & g2g1 & g2g1 & g2!g1 & !g2!g1 & !g2
Figure 5: Block dependency graph of the RLDC2 model, generated by IsamDAE. Vertices arelabeled p : Y → β → X , where: p is a propositional formula defining in which modes the block isevaluated; Y is the set of free variables for reading; β is the set of equations of the block; X is theset of dependent variables for writing. Edges are labeled by a propositional formula, defining inwhich modes the dependency applies—notation “!g” means “not g”, “g1 & g2” is the conjunction,and “g1 | g2” is the disjunction. Inria tructural Analysis of Multimode DAE Systems: summary of results References [1] A. Benveniste, B. Caillaud, and P. L. Guernic. Compositionality in dataflow synchronouslanguages: specification & distributed code generation.
Information and Computatin , 163:125–171, 2000.[2] A. Benveniste, B. Caillaud, and M. Malandain. The mathematical foundations of physicalsystems modeling languages.
Annual Reviews in Control , 2020.[3] A. Benveniste, B. Caillaud, and M. Malandain. The mathematical foundations of physicalsystems modeling languages.
CoRR , abs/2008.05166, 2020.[4] A. Benveniste, P. Caspi, S. Edwards, N. Halbwachs, P. Le Guernic, and R. de Simone. Thesynchronous languages 12 years later.
Proceedings of the IEEE , 91(1), Jan. 2003.[5] C. Berge.
The theory of graphs and its applications . Wiley, 1962.[6] B. Caillaud, M. Malandain, and J. Thibault. Implicit structural analysis of multimode DAEsystems. In , Sydney, Australia, April 2020. to appear.[7] S. L. Campbell and C. W. Gear. The index of general nonlinear DAEs.
Numer. Math. ,72:173–196, 1995.[8] H. Cartan.
Formes Diff´erentielles . Collection M´ethodes. Hermann, 1967.[9] N. Cutland.
Nonstandard analysis and its applications . Cambridge Univ. Press, 1988.[10] J. Dieudonn´e.
Fondements de l’analyse moderne . Gauthier-Villars, 1965.[11] I. S. Duff, A. M. Erisman, and J. K. Reid.
Direct Methods for Sparse Matrices . NumericalMathematics and Scientific Computation. Oxford University Press, 1986.[12] L. Hogben, editor.
Handbook of Linear Algebra . CRC Press, Boca Raton, FL, USA, 2006.[13] L. Hogben.
Handbook of Linear Algebra . 2nd edition, 2014.[14] S. E. Mattsson and G. S¨oderlind. Index reduction in Differential-Algebraic Equations usingdummy derivatives.
Siam J. Sci. Comput. , 14(3):677–692, 1993.[15] C. Pantelides. The consistent initialization of differential-algebraic systems.
SIAM J. Sci. Stat.Comput. , 9(2):213–231, 1988.[16] A. Pothen and C. Fan. Computing the block triangular form of a sparse matrix.
ACM Trans.Math. Softw. , 16(4):303–324, 1990.[17] J. D. Pryce. A simple structural analysis method for DAEs.
BIT , 41(2):364–394, 2001.[18] A. Robinson.
Nonstandard Analysis . Princeton Landmarks in Mathematics, 1996. ISBN0-691-04490-2.