Keywords: structural analysis, differential-algebraic equations (DAE), multi-mode systems, variable-structure models analyse structurelle, équations algébro-différentielles (DAE), systèmes multi-mode, modèles à structure variable

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 , x, u) = 0. This drastically facilitates modeling from first principles of the physics, as well as model reuse. In recent works [2, 3], we presented the mathematical theory needed to establish the development of compilers and 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 form amenable 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, as well as the latent equations (additional equations 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 possibly differentiated) is part of this structural analysis.

 by collecting all the needed background on structural analysis. The body of knowledge on structural analysis is large and scattered, which also motivated 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 briefly review the body of graph theory used in this context. We develop some extensions, for which we are not aware of any reference, namely the structural analysis of systems of equations with existential quantifiers.

For the structural analysis of DAE systems, we focus on John Pryce's Σ-method, that we both summarize and extend to non-square systems.

 are highlighted in this report.

Introduction

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 , x, u) = 0. This drastically facilitates modeling from first principles of the physics, as well as model reuse. In recent works [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], we presented the mathematical theory needed to establish the development of compilers and 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 form amenable 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, as well as the latent equations (additional equations 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 possibly differentiated) is part of this structural analysis. The body of knowledge on structural analysis is large and scattered, which motivated us to collect it in a single report. This report complements [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] by collecting all the needed background on structural analysis. The uses of these tools and methods in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] are highlighted in this report. Structural analysis of systems of equations: We first explain the primary meaning of the structural analysis of systems of equations, namely the study of their regularity or singularity in some generic sense. Structural information about a system of equations can be summed up by its incidence matrix; in particular, a square system yields a square matrix. If, after some permutation of row and columns, this matrix only has non-zero diagonal entries, then we say that this matrix, as well as the underlying system of equations, are structurally regular: this notion is repeatedly used in the Modelica community.

However, the most commonly used structural representation of a system of equations is actually the bipartite graph characterizing the pattern of the incidence matrix (i.e., the (i, j)-locations of its non-zero entries). Structural regularity or singularity, as well as derived notions, are then checked on 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 of graphs with sparse matrices is pervasive in linear algebra for high performance computing [START_REF] Hogben | Handbook of Linear Algebra[END_REF][START_REF] Hogben | Handbook of Linear Algebra[END_REF]. We briefly review the body of graph theory used in this context. We recall in particular the Dulmage-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. Structural analysis extends to systems of numerical equations, by considering the Jacobian of the considered system.

The link between numerical regularity and structural regularity belongs to the folklore of linear algebra; however, its understanding is critical for a correct use of structural analysis for systems of numerical equations. We were, however, unable to find a reference where related results were detailed; we propose such results here.

We also develop some useful extensions, for which we are not aware of any reference, particularly the structural analysis of systems of equations with existential quantifiers. Structural analysis of DAE systems: We focus on DAE and summarize John Pryce's Σ-method for structural analysis [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], 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 [START_REF] Campbell | The index of general nonlinear DAEs[END_REF], which we complement with its discrete-time counterpart, namely "difference arrays". These methods rely on the structural analysis of systems of equations with existential quantifiers.

Use of structural analysis on a toy example

This section is a summary of Sections 1.2.1 and 2 of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], to which the reader is referred for further details.

Inria

An ideal clutch

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. The dynamics of each shaft i = 1, 2, is modeled by ω 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 ω i denotes the time derivative of ω i . Depending on the value of the input Boolean variable γ, the clutch is either engaged (γ = t, the constant "true") or released (γ = f, the constant "false"). When the clutch is released, the two shafts rotate freely: no torque is applied to them (τ i = 0). When the clutch is engaged, it ensures a perfect join between the two shafts, forcing them to have the same angular velocity (ω 1 -ω 2 = 0) and opposite torques (τ 1 + τ 2 = 0). Here is the model:

               ω 1 = f 1 (ω 1 , τ 1 ) (e 1 ) ω 2 = f 2 (ω 2 , τ 2 )
(e 2 ) if γ then ω 1 -ω 2 = 0 (e 3 ) and τ 1 + τ 2 = 0 (e 4 ) if not γ then τ 1 = 0 (e 5 ) and τ 2 = 0 (e 6 )

When γ = t, equations (e 3 , e 4 ) are active and equations (e 5 , e 6 ) are disabled, and vice-versa when γ = f. If the clutch is initially released, then, at the instant of contact, the relative speed of the two rotating shafts jumps to zero; as a consequence, an impulse is expected on the torques.

Analyzing the two modes separately

This model yields an ODE system when the clutch is released (γ = f), and a DAE system of index 1 when the clutch is engaged (γ = t), due to equation (e 3 ) being active in this mode. The basic way of generating simulation code for the engaged mode (γ = t) is to add, to the model in this mode, the latent equation (e 3 ):

           ω 1 = f 1 (ω 1 , τ 1 ) (e 1 ) ω 2 = f 2 (ω 2 , τ 2 ) (e 2 ) ω 1 -ω 2 = 0 (e 3 ) ω 1 -ω 2 = 0 (e 3 ) τ 1 + τ 2 = 0 (e 4 ) (2)
By doing so, we can apply the following basic policy:

Regard System (2) as a system of equations in which the dependent variables are the leading variables of System (2), i.e., the torques τ i and the accelerations ω i , and consider the ω i 's as dummy derivatives, i.e., fresh variables not related to the velocities ω i ; ignore equation (e 3 ) and solve the four remaining equations for the dependent variables, using an algebraic equation solver.

(

) 3
This policy amounts to, first, bringing the DAE system back to an "ODE-like" form, then solving it as such. In this case, once the latent equation (e 3 ) was added, the Jacobian matrix of the system RT n°9387

(e 1 , e 2 , e 3 , e 4 ) with respect to its leading variables was invertible. This is the classical approach to DAE simulation via "index reduction"-having differentiated (e 3 ) once indicates that System (2) has index 1.

Basic tool 1

The current practice in DAE-based modeling languages is that adding latent equations is performed on the basis of a structural analysis, which is a semi-decision procedure using graph theoretic criteria, of much lower computational cost than that of studying the numerical invertibility of the Jacobian. We provide the bases for structural analysis of DAEs in Sections 3 and 4 of this report.

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, the mode changes are events at which restart values for states must be computed from states before the change, in one or several discrete-time transitions-note that such restart transitions are not explicitly specified in System (1).

Analyzing mode changes

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, by discretizing the real time line R ≥0 using an infinitesimal time step ∂. By doing so, the DAEs acting within 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 can be mathematically formalized by relying on Nonstandard Analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Cutland | Nonstandard analysis and its applications[END_REF], see [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] for details. In particular, System (1) is mapped to discrete-time by using the discrete-time line

T = def {n∂ | n = 0, 1, 2, . . . } (4)
and replacing in it all the derivatives by their first order explicit Euler schemes, that is:

ω 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 ≥0 is covered-infinite integers 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:

                 ω • 1 -ω 1 ∂ = f 1 (ω 1 , τ 1 ) (e 1 ) ω • 2 -ω 2 ∂ = f 2 (ω 2 , τ 2 ) (e 2 ) if γ then ω 1 -ω 2 = 0 (e 3 ) and τ 1 + τ 2 = 0 (e 4 ) if not γ then τ 1 = 0 (e 5 ) and τ 2 = 0 (e 6 ) (7)
Index reduction for the engaged mode can now be performed by shifting (e 3 ) and adding it to [START_REF] Campbell | The index of general nonlinear DAEs[END_REF], shown in red:

                     ω • 1 -ω 1 ∂ = f 1 (ω 1 , τ 1 ) (e 1 ) ω • 2 -ω 2 ∂ = f 2 (ω 2 , τ 2 ) (e 2 ) if γ then ω 1 -ω 2 = 0 (e 3 ) and ω • 1 -ω • 2 = 0 (e • 3 ) and τ 1 + τ 2 = 0 (e 4 ) if not γ then τ 1 = 0 (e 5 ) and τ 2 = 0 (e 6 ) (8)
Inria System (8) will be used to generate restart equations at impulsive mode change γ : f → t. At the considered instant, we have • γ = f and γ = t, see [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF] for the definition of shift operators. Unfolding System (8) at the two successive previous (shown in blue) and current (shown in black) instants yields, by taking the actual values for the guard at those instants into account:

previous          ω 1 -• ω 1 ∂ = f 1 ( • ω 1 , • τ 1 ) ( • e ∂ 1 ) ω 2 -• ω 2 ∂ = f 2 ( • ω 2 , • τ 2 ) ( • e ∂ 2 ) • τ 1 = 0 • τ 2 = 0 current                ω • 1 -ω 1 ∂ = f 1 (ω 1 , τ 1 ) ω • 2 -ω 2 ∂ = f 2 (ω 2 , τ 2 ) ω 1 -ω 2 = 0 (e 3 ) ω • 1 -ω • 2 = 0 τ 1 + τ 2 = 0 (9)
We regard System (9) as an algebraic system of equations with dependent variables being the leading variables of System (8) at the previous and current instants: 9) is singular since the following subsystem possesses five equations and only four dependent variables

• τ i , ω i ; τ i , ω • i for i = 1, 2. System (
ω 1 , ω 2 , • τ 1 , • τ 2 :              ω 1 -• ω 1 ∂ = f 1 ( • ω 1 , • τ 1 ) ( • e ∂ 1 ) ω 2 -• ω 2 ∂ = f 2 ( • ω 2 , • τ 2 ) ( • e ∂ 2 ) • τ 1 = 0 • τ 2 = 0 ω 1 -ω 2 = 0 (e 3 ) (10)
We propose to resolve the conflict in [START_REF] Dieudonné | Fondements de l'analyse moderne[END_REF] by applying the following causality principle:

Principle 1 (causality) What was done at the previous instant cannot be undone at the current instant.

This leads to removing, from subsystem [START_REF] Dieudonné | Fondements de l'analyse moderne[END_REF], the conflicting equation (e 3 ), thus getting the following problem for computing restart values at mode change γ : f → t:

solve for ω • 1 , ω • 2 , τ 1 , τ
               ω 1 , ω 2 , • τ 1 , • τ 2 set by previous instant ω • 1 -ω 1 ∂ = f 1 (ω 1 , τ 1 ) ω • 2 -ω 2 ∂ = f 2 (ω 2 , τ 2 ) ω • 1 -ω • 2 = 0 τ 1 + τ 2 = 0 (11)
Note that the consistency equation (e 3 ) : ω 1 -ω 2 = 0 has been removed from System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF], thus modifying 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 of the constraints in force in the new mode γ = t. Since our time step is infinitesimal, this takes zero real time. The move from System [START_REF] Cutland | Nonstandard analysis and its applications[END_REF] to System (11) must be automatized, hence we need:

Basic tool 2 We need a graph based algorithm for identifying the conflicting equations possibly occurring at a mode change. This is addressed in Section 3.2.3 of this report, where the Dulmage-Mendelsohn decomposition of a bipartite graph is presented. [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF] is not effective yet, since it involves, in the denominator of the left-hand side of the first two equations, the infinitesimal number ∂. Letting ∂ 0 in System (11) requires a different kind of technique, not related to structural analysis, to identify impulsive variables (here, the two torques) and to handle them properly. Details for these techniques are found in Sections 2.5.2 and 10 of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF].

Now, System

RT n°9387 3 Structural vs. numerical analysis of systems of equations

In this section, we develop the structural analysis for a system of algebraic equations, i.e., equations in a Euclidian space, with no dynamics (no time, no derivative), of the form

f i (y 1 , . . . , y k , x 1 , . . . , x n ) = 0, i = 1, . . . , m . ( 12
)
This system is rewritten as F (Y, X) = 0, where Y and X denote the vectors (y 1 , . . . , y k ) and (x 1 , . . . , x n ), respectively, and F is the vector (f 1 , . . . , f m ). This system has m equations, k free variables (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 1 at least. If the considered system ( 12) is square, i.e., if m = n, the Implicit Function Theorem (see, e.g., Theorem 10.2.2 in [START_REF] Dieudonné | Fondements de l'analyse moderne[END_REF]) states that, if (ν 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 , locally around ν 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 even considering it at all), by focusing instead on its structural nonsingularity, introduced in Section 3.1.

Structural nonsingularity of square matrices

Say that a subset of a Euclidian space is exceptional if it has zero Lebesgue measure. Say that a property P (x 1 , . . . , x k ) involving the real variables x 1 , . . . , x k holds almost everywhere if it holds for every x 1 , . . . , 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 {1, . . . , 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 the rows (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 of High Performance Computing. However, we were unable to find any proper reference in which it is clearly 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 1 B 2 ], where B 1 is square with 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 over some neighborhood.

Proof We consider the linear equation AX = Y , where Y ∈ R n has entries y 1 , . . . , y n , and X is the unknown vector with real entries x 1 , . . . , 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 we mean by "letting the non-zero entries of matrix A vary over a neighborhood." [START_REF] Hogben | Handbook of Linear Algebra[END_REF] 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 by lexicographic order, namely (i, j) < (i , j ) iff either i < i , or i = i ∧ j < j . Doing this defines an order preserving bijection ψ : [1, K] → P. The non-zero entries of matrix A are then collected in Inria the vector 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 = A = (a ij ) a ij = ξ ψ -1 (i,j) if (i, j) ∈ P, 0 otherwise , (14)
where the vector of R K , whose entries are the ξ k for 1≤k≤K, ranges over U . A set V of matrices obtained in this way is called a zero-pattern preserving neighborhood of A, or simply neighborhood of A when "zero-pattern preserving" is understood from the context. With this preliminary, we successively prove the two implications.

1 =⇒ 2: 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 X = Y holds almost everywhere when matrix A ranges over this neighborhood. With a suitable renumbering of the equations, we can assume that P is the identity matrix. Let V 1 be a zero-pattern preserving neighborhood of A such that the first diagonal term a 11 of A = PA remains non-zero when A varies over V 1 . Since a 11 = 0 holds when A varies over V 1 , one can perform, on A , the first step of a Gaussian elimination by using a 11 as a pivot. As x 1 is thus expressed in terms of x 2 , . . . , x m , y 1 , . . . , y n , the first equation may then be removed. This yields a reduced equation 2) , . . . , y σ(n) , and X 1 ∈ R m-1 collects the unknowns x 2 , . . . , x m . Matrix A 1 has diagonal entries equal to a 1 kk = a kk -a k1 /a 11 , for k = 2, . . . , n. Since a 11 is non-zero, the set of entries a ij of matrix A causing a 1 kk = 0 for some k is exceptional in R n×m since it requires the condition a k1 = a kk a 11 to hold. The corresponding set of entries a ij of matrix A is exceptional as well; we denote by Ξ 1 this exceptional subset of R n×m . Thus, we assume that A stays within

A 1 X 1 = Y 1 , where A 1 is an (n -1)×(m -1)-matrix, Y 1 ∈ R n-1 collects y σ(
V 1 \ Ξ 1 .
This ensures that the diagonal entries of matrix A are all non-zero. Thus, we can express x 2 in terms of x 3 , . . . , x m , y σ(2) . We are then left with a further reduced equation A 1 X 1 = Y 1 for which the same argument applies, namely: if matrix A stays within a certain zero-pattern preserving neighborhood V 2 of A, but does not belong to some exceptional set Ξ 2 , then the diagonal entries of A 1 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 1 , . . . , V n and Ξ 1 , Ξ 2 , . . . , Ξ n such that, for every matrix belonging to (

V 1 ∩• • •∩V n )\(Ξ 1 ∪Ξ 2 ∪• • •∪Ξ n ),
the values of x 1 , x 2 , . . . , 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 1 ∩ • • • ∩ V n is a zero-pattern preserving neighborhood of A and set Ξ 1 ∪ Ξ 2 ∪ • • • ∪ Ξ n is exceptional. 2 =⇒ 1:
We assume the existence of a zero-pattern preserving neighborhood V of A and an exceptional set Ξ such that A remains onto when varying over V \ Ξ. For each onto A , there exists at least one permutation matrix Q such that

A Q = B 1 B 2 (15)
where B 1 is a square invertible matrix. Neighborhood V decomposes as the union V = σ V σ , where V σ collects the matrices A for which decomposition [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF] holds with the matrix Q representing permutation σ. Since the number of permutations is finite, at least one of the V σ has non-empty interior, say, for σ = σ * . We can thus trim V to V σ * , which still is a zero-pattern preserving neighborhood 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 ∈ V . By Condition 2 of Lemma 1, the pre-image of 0 by the map V A → det(B 1 ), where det(M ) denotes the determinant of a square matrix M , is an exceptional set. Now, we have det(B 1 ) = n i=1 (-1) i+1 b i1 det(B i1 ), where B ij is the submatrix of B 1 obtained by erasing row i and column j. Since B 1 is almost everywhere invertible, there must be some i ∈ {1, . . . , n} such that b i1 = 0 and B i1 is almost everywhere invertible-otherwise we would have det(B 1 ) = 0 for any A ∈ W , where W is some neighborhood contained in V . Let P 1 be the permutation matrix exchanging rows 1 and i, so we replace B 1 by B 1 1 = P 1 B 1 . Now, since we have b 1 11 = b i1 = 0 and B 1 11 is almost everywhere invertible. This process is iterated m -1 times. For the last 1-row matrix,

RT n°9387

its diagonal entry is non-zero, otherwise the determinant of the corresponding B 1 matrix would be zero in a neighborhood of A. The wanted left permutation matrix is P = P m-1 . . .P 2 P 1 .

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-zero entries 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 some neighborhood.

Structural analysis of algebraic equations

In this section, we focus on systems of algebraic equations F (Y, X) = 0 of the form [START_REF] Hogben | Handbook of Linear Algebra[END_REF], 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 occurs from using more straightforward notations):

Notations 3 Lowercase letters (x, y, ...) denote scalar real variables; capitals (X, Y, ...) denote vectors of real variables; whenever convenient, we regard X as a set of variables and write x ∈ X to refer 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 .

As explained in the introduction of this section, the existence and uniqueness of solutions for system [START_REF] Hogben | Handbook of Linear Algebra[END_REF] 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 full development requires some background on graphs.

Background on graphs:

This short background is compiled from [START_REF] Berge | The theory of graphs and its applications[END_REF][START_REF] Hogben | Handbook of Linear Algebra[END_REF][START_REF] Hogben | Handbook of Linear Algebra[END_REF]. 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 no two edges in M are incident on a common vertex. Matching M is called complete if it covers all the vertices of G. An M-alternating path is a path of G whose edges are alternatively in M and not 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 of equations of the form [START_REF] Hogben | Handbook of Linear Algebra[END_REF]; its bipartite graph G F is the graph having F ∪X as its set of vertices and an edge (f i , x j ) if and only if variable x j occurs in function f i .

Square systems

We first consider square systems, in which the equations and the dependent variables are in equal numbers.

Definition 4 (structural nonsingularity) System F (Y, X) = 0 is called structurally nonsingular if its bipartite graph G F possesses a complete matching.

A structurally nonsingular system is necessarily square, i.e., with equations and dependent variables in equal numbers. Note that Condition 1 of Corollary 2 is a matrix reformulation of the graph theoretic condition stated in Definition 4. Thus, F is structurally nonsingular in the sense of Definition 4, if and only it its Jacobian J X F , evaluated at any pair (ν Y , ν X ) such that F (ν Y , ν X ) = 0, is structurally Inria nonsingular in the sense of Corollary 2. The link to numerical regularity is thus formalized in the following lemma:1 Lemma 5 Assume that F is of class C 1 . 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 ) remains
generically2 nonsingular when its non-zero coefficients vary over some neighborhood.

Property 2 is interpreted using the argument developed in [START_REF] Hogben | Handbook of Linear Algebra[END_REF][START_REF] Mattsson | Index reduction in Differential-Algebraic Equations using dummy derivatives[END_REF]: every n×n-matrix having the same non-zero pattern as J X F (ν Y , ν X ) identifies with a unique vector of R K , where K is the cardinal of this non-zero pattern. We denote by J X F ∈ R K the image of the Jacobian J X F (ν Y , ν X ) obtained via 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.

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 Each strongly connected component defines a block of F , consisting of the set of all equations that are vertices of this component. Blocks are partially ordered and we denote by F this order. Extending F to a total order (via topological sorting) yields an ordering of equations that puts the Jacobian matrix J X F in Block Triangular Form (BTF).

(f i , x j ) is directed from f i to x j if (f i , x j ) ∈ M, from x j to f

Non-square systems, Dulmage-Mendelsohn decomposition

In our development, we will also encounter non-square systems of algebraic equations. For general systems (with a number of variables not necessarily equal to the number of equations, n =m in [START_REF] Hogben | Handbook of Linear Algebra[END_REF]), call block a pair β = (f , x), consisting of a subsystem f =0 of system F =0 and the subset x of its dependent variables.

Definition 7 (Dulmage-Mendelsohn) For F = 0 a general system of algebraic equations, the Dulmage-Mendelsohn decomposition [START_REF] Pothen | Computing the block triangular form of a sparse matrix[END_REF] of G F yields the partition of system F = 0 into the following three blocks, given some matching of maximal cardinality for G F :

• block β O collects the variables and equations reachable via some alternating path from some unmatched equation; • block β U collects the variables and equations reachable via some alternating path from some unmatched variable; • block β E collects the remaining variables and equations.

Blocks β O , β U , β E are the Overdetermined, Underdetermined, and Enabled 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 on the 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 underdetermined blocks β O and β U are both empty.

A refined DM decomposition exists, which in addition refines block β E into its indecomposable block triangular 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 [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF]:

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 reference to the matching of maximal cardinality used for generating DM(F ). Then, denoting by F the remaining subsystem of F , we would still get that no overdetermined block exists in the Dulmage-Mendelsohn decomposition of F , and F ⊇ (F \ β O ) (the inclusion being strict in general). However, this policy depends on the particular matching. In contrast, the rule defined by Lemma 9 is independent from the matching. This Lemma is used in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] to generate restart code at mode changes, see Section 5.3.2 of this report; the independence with respect to the particular matching considered was crucial to ensure the determinism of code generation.

Practical use of structural analysis

We distinguish between two different uses of structural analysis. The classical use assumes that we already 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

Local use of structural analysis

If a square matrix is structurally singular, then it is singular. The converse is false: structural nonsingularity does not guarantee nonsingularity. Therefore, the practical use of Definition 4 is as follows: We first check if F = 0 is structurally nonsingular. If not, then we abort searching for a solution. Otherwise, we then proceed to computing a solution, which may or may not succeed depending 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 the beginning of this section. Therefore, the classical use of structural analysis is local: let F (Y, X) = 0 be 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 DAE solvers, since a close guess of the solution is known from the previous time step.

Non-local use of structural analysis

We are, however, also interested in invoking the structural nonsingularity of system F (Y, X) = 0 when no close guess is known. This is the situation encountered when handling mode changes of multimode DAE systems [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF]. We thus need another argument to justify the use of structural analysis in this case.

As a prerequisite, we recall basic definitions regarding smooth manifolds (see, e.g., [START_REF] Cartan | Formes Différentielles. Collection Méthodes[END_REF], Section 4.7). In the next lemma, we consider the system F (Y, X) = 0 for a fixed value of Y , thus Y is omitted.

Lemma 11 Let k, m, n ∈ N be such that m < n and set p = n -m. For S ⊂ R n and x * = (x * 1 , . . . , 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 such that 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 0 in R p and a homeomorphism ψ : U → S ∩ V , such that ψ, seen as a map with values in R n , is of class C k and has rank p at 0-i.e., ψ (0) has rank p.

Statement 1 means that S is, in a neighborhood of x * , the solution set of the system of equations f 1 (X) = 0, . . . , f m (X) = 0 in the n-tuple of dependent variables X, where F = (f 1 , . . . , f m ). Statement 2 expresses that ψ is a parameterization of S in a neighborhood of x * , of class C k and rank 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 study the case of square systems of equations (i.e., let m = n). Consider a system F (X) = 0, where X = (x 1 , . . . , x n ) is the n-tuple of dependent variables and

F = (f 1 , . . . , 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 more function in such a way that the resulting (n+1)-tuple Ḡ is a C k -diffeomorphism, from V to an open set W of R n+1 . The extended system Ḡ = 0 satisfies Property 1 of Lemma 11. By Property 2 of Lemma 11, (S × R) ∩ V has dimension 1, implying that S ∩ proj R n (V ) has dimension < 1. This analysis is summarized in the following lemma:

RT n°9387

Lemma 12 Consider the square system F (X) = 0, where X = (x 1 , . . . , x n ) is the n-tuple of dependent variables and F = (f 1 , . . . , f n ) is an n-tuple of functions R n → R of class C k , k≥1. Let S ⊆ R n be the solution set of this system. Assume that system F = 0 has solutions and let x * be such a solution. Assume further that the Jacobian matrix J X F (x) is nonsingular at x * . Then, there exists an open neighborhood V of x * in R n , such that S ∩ V has dimension < 1.

The implications for the non-local use of structural analysis are as follows. Consider the square system F (X) = 0, where X = (x 1 , . . . , x n ) is the n-tuple of dependent variables and F = (f 1 , . . . , 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 proceed to solving the system and the following cases can occur: 1. the system possesses no solution; 2. its solution set is nonempty, of dimension < 1; or 3. its solution set is nonempty, of dimension ≥ 1.

Based on Lemma 12, one among cases 1 or 2 will hold, generically, whereas case 3 will hold exceptionally. Thus, structurally, we know that the system is not underdetermined. This is our justification of the non-local use of structural analysis. Note that the existence of a solution is not guaranteed. Furthermore, unlike for the local use, uniqueness is not guaranteed either. Lemma 12 only states that there cannot be two arbitrarily close solutions of F = 0, as its conclusion S ∩ V = {x * } involves an open neighborhood V of x * . Of course, subclasses of systems for which existence and uniqueness are guaranteed are of interest.

Equations with existential quantifiers

To support the structural analysis for the method of differential arrays advocated by Campbell and Gear [START_REF] Campbell | The index of general nonlinear DAEs[END_REF] (see also Section 4.3 of this report), we will need to develop a structural analysis for the following class of (possibly non-square) algebraic systems of equations with existential quantifier:3

∃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 their elimination 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 the satisfaction of (19), given a value for tuple Y . At this point, a fundamental difficulty arises. Whereas (19) is well defined as an abstract relation, it cannot in general be represented by a projected system of smooth algebraic equations of the form G(X, Y ) = 0. Such a G can be associated to (19) only in subclasses of systems. 4 In general, 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 reformulate our requirement differently. Say that ν Y is consistent if the system of equations F (X, W, ν Y ) = 0 possesses a solution for (X, W ).

(20)

Consider the following property for the systen F (X, W, Y ) = 0: for every consistent ν Y ,

F (ν 1 X , ν 1 W , ν Y ) = 0 F (ν 2 X , ν 2 W , ν Y ) = 0 =⇒ ν 1 X = ν 2 X , (21)
expressing that X is independent of W given a consistent tuple of values for Y . To find structural criteria guaranteeing (21), we will consider the algebraic system F (X, W, Y ) = 0 with X, W, Y as dependent 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, almost everywhere when the non-zero coefficients of the Jacobian matrix J X,W,Y F vary over some neighborhood: 5

• consistent values for Y exist;

• condition (21) holds.

The structural nonsingularity of (19) can be checked by using the DM decomposition of F as follows. Let (β U , β E , β O ) = DM(F ) be the DM decomposition of F (X, W, Y ) = 0, with (X, W, Y ) as dependent variables. We further assume that the regular block β E is expressed in its block triangular form, see Lemma 6 and comments thereafter. Let B be the set of all indecomposable blocks of β E and be the partial order on B following Lemma 6. The following holds: 3. For every β ∈ B containing some variable from X, then, β contains no variable from W , and for every β ≺β and every directed edge (z , z)∈ G M F such that z ∈β , z∈β, and M is an arbitrary complete matching for G F , then z ∈ W .

Lemma
Conditions 1 and 2 speak by themselves. Regarding Condition 3, the intuition is the following. Every directed path of G M F , originating from W and terminating in X, of minimal length, must traverse 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. By condition (3), each indecomposable block β ∈ B involving a variable of X has the form G(Z, U ) = 0, where 1. Z is the n-tuple of dependent variables of G = (g 1 , . . . , g n ), 2. no variable of W belongs to Z, and 3. U ⊆ X∪Y is a subset of the dependent variables of the blocks β immediately preceding β.

Since G = 0 is structurally regular, fixing a value for U entirely determines all the dependent variables of block β. This proves the "if" part.

"Only if " part:

We prove it by contradiction. If condition 1 does not hold, then the existence of consistent values for Y is not structurally guaranteed. If condition 2 does not hold, then the variables of X that belong to β U are not determined, even for given values of W, Y . If condition 3 does 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 5 The precise meaning of this statement is the same as in Lemma 5, see the discussion therafter.

RT n°9387

1. Z is the n-tuple of dependent variables of G = (g 1 , . . . , g n ), and Z contains a variable x ∈ X, 2. no variable of W belongs to Z, and 3. some variable w ∈ W belongs to U , the set of the dependent variables of the blocks β immediately preceding β.

Hence, w influences x, structurally, which again violates (21). This finishes the proof.

Algorithm 1 ExistQuantifEqn

Require:

F ; return (b O , b U , F Σ , F Σ ) 1: (β U , β E , β O ) = DM(F ) 2:
if condition 1 of Lemma 14 holds then partition

β E = F Σ ∪ F Σ 7: else 8: b U ← f 9: else 10: b O ← f
We complement Lemma 14 with the algorithm ExistQuantifEqn (Algorithm 1), which requires a system 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 the value 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 next section, we recall the background on the structural analysis of (single-mode) DAE systems.

Structural analysis of DAE systems

In this section, we consider DAE systems of the form

f j (x i 's and derivatives) = 0 , (22)
where x 1 , . . . , x m denote the dependent signal variables (their valuations are trajectories) and f 1 = 0, . . . , f n = 0 denote the equations-with reference to notation [START_REF] Hogben | Handbook of Linear Algebra[END_REF], we omit the input signals collected in Y . By analogy with DAEs, we use the acronym dAE to mean difference Algebraic Equations, which define discrete-time dynamical systems of the form f j (x i 's and shifts) = 0 .

where x 1 , . . . , x m denote the dependent stream variables (their valuations are data streams, which can be regarded as real sequences), f 1 = 0, . . . , f n = 0 denote the equations, and, for x = {x n | n = 1, 2, . . . } a stream of variables, its successive shifts x •k , where k is a nonnegative integer, are defined by

x •k n = def x n+k ; we write for short x • = def x •1 . (24)
Also, the notation x k is adopted throughout this paper, instead of the more classical x (k) , for the k-th derivative of x.

Inria

John Pryce's Σ-method for square systems

The DAE systems we consider are "square", i.e., have the form (22), with m = n. Call leading variables of System (22) the d i -th derivatives x d i i for i = 1, . . . , n, where d i is the maximal differentiation degree of variable x i throughout f 1 = 0, . . . , f n = 0. The problem addressed by the structural analysis of DAE systems of the form (22) is the following. Regard (22) as a system of algebraic equations with the leading variables as unknowns. If this system is structurally nonsingular, then, given a value for all the x k i for i = 1, . . . , n and k = 0, . . . , d i -1, a unique value for the leading variables can be computed, structurally; hence, System (22) is "like an ODE". If this is not the case, finding additional latent equations by differentiating suitably selected equations from (22) will bring the system to an ODE-like form, while not changing its set of solutions. Performing this is known as index reduction.

Algorithms were proposed in the literature for doing it efficiently. Among them, the Pantelides' algorithm [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF] is the historical solution. We decided, however, to base our subsequent developments on the beautiful method proposed in [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF] by J. Pryce, called the Σ-method. The Σ-method also covers the construction of the block triangular form and addresses numerical issues, which we do not discuss here.

Weighted bipartite graphs: We consider System (22), which is entirely characterized by its set of 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 differentiation degree. 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 maximal differentiation degree of variable x in function f . This yields a weight for any matching M of G by the formula w(M) = (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 ≥ 0 for all f ∈ F . ( 25
)
Then, differentiating c f times each function f yields a DAE system F Σ = 0 having the following properties. 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 dx . Consequently, system F Σ = 0, now seen as a system of algebraic equations having x dx as dependent 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. 6Definition 15 For F a DAE system, the solution to Problem 1 yields the DAE system

F Σ = {f c f | f ∈F } together with the system of consistency constraints F Σ = {f k | f ∈F, 0≤k<c f }.
Knowing the offsets also allows transforming F into a system of index 1, by not performing the final round of differentiations. There are infinitely many solutions to Problem 1 with unknowns c f and d x , since, for example, adding the same ∈ N ≥0 to all c f 's and d x 's yields another solution. We thus seek for a smallest solution, elementwise. Hence, Problem 1 is the fundamental problem we must solve, and we seek a smallest solution for it. The beautiful idea of J. Pryce is to propose a linear program encoding Problem 1. As a preliminary step, we claim that a bruteforce encoding of Problem 1 is the following: Find nonnegative integers

ξ fx , d x , c f such that f :(f,x)∈E G ξ fx = 1 x:(f,x)∈E G ξ fx = 1 ξ fx ≥ 0    complete matching d x -c f -d fx ≥ 0 c f ≥ 0 encodes "≥" in (25) (f,x)∈E G ξ fx (d x -c f -d fx ) = 0 encodes "=" in (25) (26)
where, for the formulas having no summation: x, f , and (f, x) range over X, F , and E G respectively. We now justify our claim. Focus on the first block. Since the ξ's are nonnegative integers, they can only take values in {0, 1} and one defines a subgraph of graph G by only keeping edges (f, x) such that ξ fx = 1. The first two equations formalize that the chosen subset of edges is a matching, which is complete since all vertices of G are involved. The second block is a direct encoding of (25) if we ignore the additional statement "with equality iff". The latter is encoded by the last constraint (since having the sum equal to zero requires that all the terms be equal to zero). Constraint problem (26) does not account for our wish for a "smallest" solution: this will be handled separately. Following the characterization of solutions of linear programs via complementary slackness conditions, every solution of problem ( 26) is a solution of the following dual linear programs (LP), where the ξ fx and the c f are real:

primal : maximize (f,x)∈E G d fx ξ fx subject to f :(f,x)∈E G ξ fx = 1 and x:(f,x)∈E G ξ fx ≥ 1 and ξ fx ≥ 0 (27) dual : minimize x d x -f c f subject to d x -c f ≥ d fx and c f ≥ 0 (28)
where, for the formulas having no summation: x, f , and (f, x) range over X, F , and E G respectively. In these two problems, f ranges over F , x ranges over X, and (f, x) ranges over G. Also, we have relaxed the integer LP to a real LP, as all solutions to the integer LP are solutions of the real LP. Note that LP (27) encodes the search for a complete matching of maximum weight for G. By the principle of complementary slackness in linear programming, for respective optima of problems ( 27) and (28), ξ fx >0 if and only if

d x -c f = d fx , ( 29
)
which is exactly the last constraint of (26). Using this translation into linear programs ( 27) and (28), it is proved in [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], Thm 3.6, that, if a solution exists to Problem 1, then a unique elementwise smallest solution exists. Based on the above analysis, the following Algorithm 2 (FindOffsets) was proposed in [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF] for solving (27,28) and was proved to provide the smallest (real) solution for (28), which happens to be integer. The reason for using the special iterative algorithm for solving the dual LP (28) is that a standard LP-solving algorithm will return an arbitrary solution, not necessarily the smallest one. The following lemma, whic is repeatedly used in Section 11.1 of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] about the structural analysis of multimode DAE systems, is an obvious consequence of the linearity of problems ( 27) and (28):

Lemma 16 Let G be a given bipartite graph and let two families of weights

(d 1 fx ) (f,x)∈E G and (d 2 fx ) (f,x)∈E G be related by d 2 fx = M × d 1 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 2 x = M × d 1
x for every variable x, and

c 2 f = M × c 1 f for every function f . Inria 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 method for solving LP can be used.

2. Apply the following iteration until a fixpoint is reached (in finitely many steps), from the initial 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 complete matchings for G is non-empty. This criterion can be evaluated prior to computing the optimal offsets. This completes our background material on structural analysis.

The Σ-method for non-square systems

Structural analysis must be extended to non-square systems if we wish to allow for guards being evaluated at the current instant-in this case, enabled/disabled equations are progressively determined along with the progressive evaluation of guards, see Appendix A of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF]. Structural analysis is also a powerful tool for addressing problems other than just the simulation of DAE systems, e.g., the generation of failure indicators for system diagnosis. For both uses, we will need to address other cases than square DAE systems. We develop here an adaptation of the Σ-method to non-square DAE systems. In our development, we could simply copy the original constructions, lemmas, and proofs, of [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], 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 nodes and variable nodes in different numbers, hence the concept of complete matching does not apply. In this case, we consider instead the following weaker notion: we say that a matching M is equation-complete if it covers all the equation nodes.

(

No equation-complete matching exists for an overconstrained system F = 0, i.e., a system with less variables 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 we can complete this matching to find a feasible solution. Hence, the existence of at least one matching is 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 variable in 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 the entire system, namely d x .

Call F Σ the DAE system defined by H and denote by Z = {x dx | x∈X} the set of its leading variables. 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 we know that the overconstrained block β O is empty.

RT n°9387

Justification of the extension

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 computing the offsets?

Question 2 Is existence and uniqueness of the smallest solution still guaranteed in the non-square case?

Addressing Question 1: Step 2a of Algorithm 2 is not modified.

Step 2b, on the other hand, cannot be 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 the resulting algorithm.

Addressing Question 2: This question is answered by the following lemma:

Lemma 17 (non-square Σ-method) Algorithm 2 converges to a fixpoint if and only if H is an optimal solution of the primal problem (27). Let c * , d * be this fixpoint when it exists. Then c * , d * is an integer-valued optimal solution for the dual problem (28) and, for any optimal solution c † , d † for this dual problem, c † ≥ c * and d † ≥ d * hold.

The last statement expresses the uniqueness of the smallest optimal solution. Lemma 17 was stated and proved in [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], albeit for square DAE systems. Here, we will be carefully following the steps of that proof, while checking that they suit the non-square case.

Proof Following [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], let φ : c → φ(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 dual problem (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) = max (f,x)∈E G (d fx + c f ) x∈X = { 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 * be any fixpoint for φ (note that φ depends on H). We have d *

x ≥ d fx + c * f for any (f, x) ∈ E G , with equality if (f, x) ∈ H. Assume that c * ≥ 0 also holds; as a matter of fact, c * , d * is a feasible solution of the dual problem (28). Finally, H is an optimal solution of the primal problem, and c * , d * an optimal solution of the dual problem, because they satisfy the slackness conditions (29), which are known to characterize optimal solutions. By contraposition, if H is not optimal, then no fixpoint of φ exists.

It remains to be proved that c * ≥ 0 holds. By construction, φ is nondecreasing: c ≤ c implies φ(c ) ≤ φ(c). Denote by c 0 , d 0 , c 1 , d 1 . . . the alternating sequence of c's and d's obtained by looping over steps (a) and (b), alternatively. We have c k = φ(c k-1 ) for k > 0. If c 0 = 0, then d 0 ≥ 0 and c 1 f = d 0 x -d fx ≥ 0, where x is any variable such that (f, x) ∈ H (the value d 0 x -d fx is independent from this choice). As a consequence, c 1 ≥ 0. Since φ is nondecreasing, applying φ to both sides of c 1 ≥ c 0 = 0 yields c 2 ≥ c 1 , and so on, so that the algorithm yields a nondecreasing sequence, showing that c * ≥ 0.

Finally, if c † , d † is another dual-optimal solution for the offsets, then c † is a fixpoint of φ and thus c † = φ(c † ) ≥ φ(0), hence c † ≥ φ k (0) for every k, whence c † ≥ c * follows.

Link with the Pantelides algorithm

In this section, we closely follow Section 5.3 of [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], regarding the comparison of the Σ-method with Pantelides' method [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF] 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 finds a minimally structurally singular (MSS) subset of equations and differentiates each equation in this set, creating an augmented system. This is repeated until no more MSS are found. Pantelides writes the system as n + m equations F (X, Z) = 0, where X stands for the n variables whose derivatives appear, Y stands for the m algebraic variables, and Z = (X , Y ). Pantelides' construction applies to this reformulation. For our comparison, the Z coincides with the collection of leading variables x dx 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 smallest optimal offsets c * f , d * x satisfying

x d * x -f c * f = w(H * ) (the total weight of H * , defined as (f,x)∈H * d fx ). Let 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 ) as the set

L = {(f, x) | d fx = d x }
For some equation f , respectively some subset E of equations, consider

L f = {x | (f, x) ∈ L} , and L(E) = 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 the following results. This means that H is optimal for the primal problem. We have

d * x -c * f = d fx = d x for (f, x) ∈ H hence c * f = 0, d * x = d
x is the (unique) smallest optimal solution for the dual problem. ¬2a⇒ ¬2c: If 2a fails to hold, then the optimal solution H * for the primal problem contains a (f, x) with d fx < d x . Then, by the slackness condition (29):

d * x = d fx + c * f < d x + c * f . ( 32
)
RT n°9387

= d xo -d f xo = d * xo -d f xo ≥ c * f ≥ 0 whence c * f = 0.
Lemma 20 For E any MSS, we have ∀f ∈ E : c f > 0.

Proof Suppose the lemma is false, meaning that

E 0 = def {f ∈E | c f = 0}
L(E \ E 0 ) ∩ [∃f. ((E 0 ×X) ∩ H * )] = ∅ (35)
The set on the right-hand side of ∩ is contained in L(E) (by statement 1 of Lemma 19) and has cardinality at least equal to |E 0 |. By (35), we get

|L(E \ E 0 )| ≤ |L(E)| -|E 0 | < |E \ E 0 |, hence E \ E 0 is a SS subset strictly contained in E, thus contradicting the minimality of E.
Following these results, Pantelides' algorithm can be extended to non-square systems in pretty much the same way the Σ-method was extended. This non-square Pantelides' algorithm locates MSS 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 found by 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 as the Σ-method, in that the same offsets c i are found.

Inria

Differential and difference arrays

Differential arrays were proposed in [START_REF] Campbell | The index of general nonlinear DAEs[END_REF] 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 ≥0 , by

A k = def         F d dt F d 2 dt 2 F . . . d k dt k F         (36)
where d dt denotes the total time derivative. Thus, for k sufficiently large, array system A k = 0 should contain all the latent equations that must be added to F = 0 in order to reduce its index to 0. More precisely, the index reduced counterpart of F = 0 is the following system with existential quantifiers:

∃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 time derivatives d j dt j F by the j-th forward shift F •j defined in (24), when building array (36):

A k = def        F F • F •2 . . . F •k        (38)
The structural analysis of systems with existential quantifiers, of the form (19), was studied in Section 3.4. We can apply it to the structural analysis of arrays (36) or (38) to find the latent equations 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 that the rows of the array are successive shifts of the same dAE system. In turn, the method of Section 3.4 extends to time-varying discrete-time systems. We just need to replace, in (38), the shifted versions of F by the successive F (t 0 ), . . . , F (t k ) at successive discrete instants t 0 , . . . , t k . This is used in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] to develop a structural analysis of finite cascades of mode changes.

5 Summary of results of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF]

regarding multimode DAE systems

In this section, we briefly review the results of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] regarding multimode DAE systems simulation. We indicate where, in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], the developments of Sections 3 and 4 of this report are used.

RT n°9387

Multimode DAE and dAE systems definition

In [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], a mathematical definition for multimode DAE systems is provided, and we also define their discrete-time counterpart, namely multimode dAE systems, of the respective forms multimode DAE: if γ j (x i 's and derivatives) then f j (x i 's and derivatives) = 0 multimode dAE: if γ j (x i 's and shifts) then f j (x i 's and shifts) = 0 ,

where x •k , the k-shift of x, is defined in (24).

Nonstandard semantics

To allow for a uniform handling of modes and mode change events in the structural analysis of multimode DAE, we use a nonstandard semantics for DAEs and multimode DAEs. In this nonstandard semantics, derivative x 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 forward shift 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", which can be given a formal meaning in nonstandard analysis [START_REF] Cutland | Nonstandard analysis and its applications[END_REF][START_REF] Robinson | Nonstandard Analysis[END_REF].

This mapping allows for a discretization of multimode DAE systems with an infinitesimal error, which yields a faithful discrete-time approximation. Of course, there is no free lunch: the nonstandard semantics is not effective in that there is no computer that can execute it. Still, it is highly useful for the structural analysis, which is a symbolic analysis. The nonstandard semantics of a multimode DAE 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 the change and its value x + right after the change, encoded in the nonstandard semantics as x and x • .

Structural analysis of multimode DAE/dAE systems

The structural analyses of both DAE and dAE systems coincide, by replacing the differentiation operator x → x arising in DAEs, by the forward shift operator x → x • arising in dAEs. Pryce's Σ-method applies to both DAE and dAE; for the latter, the equation offsets c f indicate how many times each equation f must be shifted.

Handling continuous modes

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 Σ-method presented in Section 4. This is summarized in the following Basic tool 3 Standard structural analysis can be performed for each continuous mode, in order to add the needed latent equations. The continuous-time structural analysis (where latent equations are found by differentiation) is therefore used to generate code within continuous modes.

Inria

Handling mode changes

Basic tool 4 In order to prepare for the generation of restart code at mode changes, the discretetime structural analysis (where latent equations are found by shifting) based on the nonstandard semantics 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 for each continuous mode, we may have a conflict, at mode change events separating two successive continuous 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 generated by the structural analysis, which may be conflicting with the dynamics of the continuous mode before the change. Since one cannot regard this model as being incorrect, we need a policy for handling these conflicts.

Handling isolated mode changes separating two successive continuous modes: In [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], it is proposed to call causality to the rescue: for isolated mode changes, we give priority to the previous dynamics and postpone for a while the conflicting consistency equations generated by the new 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 by using the Dulmage-Mendelsohn decomposition (Definition 7 and associated Lemmas 8 and 9).

Conflicting consistency equations of the new continuous mode are postponed for a finite number of nonstandard instants. The resulting discrete-time system operates for finitely many nonstandard instants, and is the basis for recovering the restart values of the states in the new mode. This establishes the structural analysis of mode changes separating two successive continuous modes.

Handling finite cascades of mode changes: We also handle in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] finite cascades of successive mode changes separating two successive continuous modes. Since the dynamics varies over the successive instants of the cascade, time-invariance no longer holds, and, thus, the Σ-method does not apply. However, the method of difference arrays, closely derived from Campbell and Gear's differential arrays [START_REF] Campbell | The index of general nonlinear DAEs[END_REF], can be adapted to time-varying discrete-time systems, and invoked. See Section 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 leading to a new continuous mode, however, we form the associated difference array by stacking, one below the other, the dynamics that holds at each successive instant of the cascade. More precisely, let t be the 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 second event 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°9387

cascade are uniquely determined as functions of the states before the change. As for differential arrays, the extra variables brought by the successive shifting must be eliminated by existential quantification-see the comments regarding formula [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF] introducing differential arrays in [START_REF] Campbell | The index of general nonlinear DAEs[END_REF]. Finite cascades of mode changes are thus handled by using the results of Section 3.4 regarding equations with 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 3.4 regarding equations with existential quantifiers.

Fixpoint between guards and the equations they control: In a multimode DAE system, a logico-numerical fixpoint occurs when a variable x is determined by a different dynamics depending on the value of a guard that itself depends, directly or indirectly, on x. In [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], only multimode DAE systems having no logico-numerical fixpoint are supported. In such models, the above guard should rather depend on the left-limit x -of x, where x -(t) = lim s t x(s). Nevertheless, in Appendix A of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], we propose a method for handling multimode DAE systems with fixpoints but no cascades of length > 1. This method uses the following basic tool:

Basic tool 7 The Σ-method for non-square systems, presented in Section 4.2.

Generic form of the generated simulation code

For S a multimode DAE model having continuous modes separated by finite cascades of mode change events, the method proposed in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] returns the following:

Case of a structurally incorrect model: The compilation algorithm returns, for each mode and mode change event in which structural analysis fails, the overdetermined subsets of equations and the underdetermined subsets of variables.

Case of a structurally correct model, producing an mDAE interpreter : The compilation algorithm 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 from state variables, and expressions exp determining the values of guards as function of state and leading 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 one of 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 associate a 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 β 1 , . . . , β k , . . . , β K by setting λ(π) = K 1 λ(β k ). Every circuit of G has label f, which ensures that the code is not globally circular. Also, for any two blocks β 1 = β 2 leading to the same variable x in G, we have λ(β 1 ) ∧ λ(β 2 ) = f, which ensures that two different blocks never compete at determining the same variable x. This interpreter 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 of the Signal synchronous language [START_REF] Benveniste | Compositionality in dataflow synchronous languages: specification & distributed code generation[END_REF][START_REF] Benveniste | The synchronous languages 12 years later[END_REF].

The notion of mDAE interpreter is illustrated in Figure 5, for an RLDC2 circuit example with schematic shown in Figure 3 and model given in System (42), see Section 12 of [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] for details. This form of mdAE interpreter is used in the IsamDAE tool [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF], by representing all the label using BDD-related data structures.

0 = i 1 + i 2 + j 1 + j 2 (K 1 ) x 1 + w 1 = u 1 + v 1 (K 2 ) u 1 + v 1 = u 2 + v 2 (K 3 ) u 2 + v 2 = x 2 + w 2 (K 4 ) w 1 = L 1 • j 1 (L 1 ) w 2 = L 2 • j 2 (L 2 ) i 1 = C 1 • v 1 (C 1 ) i 2 = C 2 • v 2 (C 2 ) x 1 = R 1 • j 1 (R 1 ) x 2 = R 2 • j 2
(R 2 ) s 1 = if γ 1 then i 1 else -u 1 (S 1 ) s 2 = if γ 2 then i 2 else -u 2 (S 2 ) 0 = if γ 1 then u 1 else i 1 (Z 1 ) 0 = if γ 2 then u 2 else i 2 (Z 2 )

γ 1 = (s - 1 ≥ 0) (γ - 1 ) γ 2 = (s - 2 ≥ 0) (γ - 2 ) (42)

Conclusion

This report complements references [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] by collecting and providing all missing details related to structural analysis, for both algebraic and DAE systems, in a systematic way.

RT n°9387 p is a propositional formula defining in which modes the block is evaluated; Y is the set of free variables for reading; β is the set of equations of the block; X is the set of dependent variables for writing. Edges are labeled by a propositional formula, defining in which modes the dependency applies-notation "!g" means "not g", "g1 & g2" is the conjunction, and "g1 | g2" is the disjunction.

Inria

Figure 1 :

 1 Figure 1: An ideal clutch with two shafts.

  i otherwise. Denote by G M F the resulting directed graph. The following result holds ([11], Chapter 6.10): Lemma 6 The strongly connected components of G M F are independent of M.

Figure 2 : 1 F 2 F 1 Fand G M 2 F

 21212 Figure 2: Illustrating Lemma 6 (we only show the projection of strongly connected components on function nodes).

Inria

Lemma 18 1 . 1 2.

 11 If E is MSS with k elements, then |L(E)| = k -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 0. Proof For statement 1, if L(E) has fewer than k -1 elements, then removing any element from E yields a set E satisfying |L(E )| ≤ |L(E)| < k -1 = |E |. Therefore, E is SS, which contradicts the minimality of E. The equivalence 2a⇔2b is just Hall's Theorem applied to the sets L e . 2a⇒2c: If L has an equation-complete matching H, then (by definition of L) for (f, x)∈H, d fx = d x . In other words, d fx is maximal among all the d f x for f ranging over the set of all equations.

7

 7

7

  BDD means: Binary Decision Diagram, a data structure developed for model checking. Inria

Figure 3 :

 3 Figure 3: Schematics of the RLDC2 circuit.

Figure 4 :

 4 Figure 4: The RLDC2 model (n.b.: variable y denotes the time derivative of y).

Figure 5 :

 5 Figure5: Block dependency graph of the RLDC2 model, generated by IsamDAE. Vertices are labeled p : Y → β → X, where: p is a propositional formula defining in which modes the block is evaluated; Y is the set of free variables for reading; β is the set of equations of the block; X is the set of dependent variables for writing. Edges are labeled by a propositional formula, defining in which modes the dependency applies-notation "!g" means "not g", "g1 & g2" is the conjunction, and "g1 | g2" is the disjunction.

  On the other hand, since d * y , c * f is a dual-optimal solution, we have d * y ≥ d f y + c * f ≥ d f y for all y, whence In particular, d * y ≥ d y ; combining this result with (32) yields c f > d * x -d x ≥ 0, so that 2c is false.Lemma 19 Let H * be a primal-optimal solution as before. If c * fo = 0 is satisfied for some equationf o , then 1. 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 = f o such that (f, x o ) ∈ L. foxo = d xo, which proves statement 1. For statement 2, if (f, x o ) ∈ L for some f , then 0

	d * y ≥ d y .	(33)
	Proof We have
	d foxo = d * xo (since d * is optimal and c * fo = 0)
	≥ d xo (by (33))
	≥ d foxo (by (31))
	hence we get
	d * xo = d xo	(34)
	and d

  is nonempty. Let |E| = k and |E 0 | = l, so that 1≤l≤k. By statement 1 of Lemma 19, we have L(E 0 ) ≥ l. Since |L(E)| < |E| by the assumption that E is an MSS, we deduce l < k, hence E \ E 0 is nonempty. Pick any f ∈ E \ E 0 , then c f > 0. By statement 2 of Lemma 19, L f cannot contain any x that is reachable from E 0 via H

* :

This result is implicitly used in[START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF][START_REF] Mattsson | Index reduction in Differential-Algebraic Equations using dummy derivatives[END_REF][START_REF] Pothen | Computing the block triangular form of a sparse matrix[END_REF][START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], as well as in a major part of the literature on the structural analysis of DAE systems.

Generically means: outside a set of values of Lebesgue measure zero. RT n°9387

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-but expensive-using Gröbner bases. Inria

We should rather say the differentiation index as, once again, other notions of index exist for DAEs[START_REF] Campbell | The index of general nonlinear DAEs[END_REF] that are not relevant to our work. Also note that the standard definition of the differentiation index, used in Pryce's article[START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], slightly differs from the one that we adopted in this report for the sake of clarity.RT n°9387

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-0803