Keywords: structural analysis, differential-algebraic equations (DAE), multi-mode systems, variable-structure models, nonstandard analysis

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 the reuse of models. In this paper, we develop the mathematical theory needed to establish the development of compilers and tools for DAE-based physical modeling languages on solid mathematical bases. Unlike Ordinary Differential Equations (ODEs, of the form x = g(x, u)), DAEs exhibit subtle issues because of the notion of differentiation index and related latent equations-ODEs are DAEs of index zero, for which no latent equation needs to be considered. Prior to generating execution code and calling solvers, the compilation of such languages requires a nontrivial structural analysis step that reduces the differentiation index to a level acceptable by DAE solvers. The models supported by tools of the Modelica class involve multiple modes, with modedependent DAE-based dynamics and state-dependent mode switching. However, multimode DAEs are much more difficult to handle than DAEs, especially because of the events of mode change. Unfortunately, the large literature devoted to the mathematical analysis of DAEs does not cover the multimode case, typically saying nothing about mode changes. This lack of foundations causes numerous difficulties to the existing modeling tools. Some models are well handled, others are not, with no clear boundary between the two classes. In this paper, we develop a comprehensive mathematical approach supporting compilation and code generation for this class of languages. Its core is the structural analysis of multimode DAE systems. As a byproduct of this structural analysis, we propose sound criteria for accepting or rejecting multimode models. Our mathematical development relies on nonstandard analysis, which allows us to cast hybrid system dynamics to discrete-time dynamics with infinitesimal step size, thus providing a uniform framework for handling both continuous dynamics and mode change events.

Fondements mathématiques des langages de modélisation de systèmes physiques

Résumé : Les langages modernes de modélisation de systèmes physiques, tels que Modelica, Amesim, ou Simscape, s'appuient sur des Équations Algébro-Différentielles (DAE), qui sont des contraintes de la forme f (x , x, u) = 0. Cette approche rend naturelle la modélisation directe à partir des principes de la physique, ainsi que la réutilisation de modèles. Dans cet article, nous développons les bases mathématiques qui fondent ces langages.

Par rapport aux Équations Différentielles Ordinaires (ODE, de la forme x = g(x, u)), les DAE sont plus compliquées, en raison des notions d'index et d'équations latentes-les ODE sont des DAE d'index zéro, sans équations latentes. Avant toute génération de code, une analyse structurelle des modèles est requise, afin de ramener l'index à des valeurs acceptables pour les solveurs de DAE (deux ou trois maximum).

Les modèles acceptés par les langages de la classe Modelica sont des DAE multi-modes, avec une DAE différente dans chaque mode, tandis que les changements de mode sont provoqués par des conditions portant sur l'état et/ou le temps. Les DAE multi-mode sont beaucoup plus difficiles à traiter que les DAE (mono-mode). La difficulté principale est le traitement des changements de mode et des conditions de redémarrage associées. Ce point n'est pas abordé dans la littérature mathématique sur les DAE multi-mode: rien n'est dit sur comment les changements de mode doivent être traités. Cette situation est source de problèmes importants pour tous les outils existants. Certains modèles sont correctement traités, d'autres non, sans qu'on sache bien pourquoi-et les cas problématiques ne sont pas pathologiques, mais se rencontrent naturellement.

Dans ce travail, nous développons une approche cohérente et complète pour la compilation des DAE multi-mode. Au coeur de celle-ci se trouve une analyse structurelle multi-mode, qui couvre à la fois les dynamiques dans les divers modes, et les changements de mode. Notre approche traite des modèles présentant des valeurs impulsives pour certaines variables, lors des changements de mode. Notre approche permet d'accepter ou de rejeter des modèles sur des critères bien fondés (modèle sur-ou sous-déterminé). Ces travaux eussent été impossibles sans le recours à l'analyse non-standard. Grâce à elle, nous pouvons réinterpréter les dynamiques en temps continu comme du temps discret à pas infinitésimal, ce qui fournit une vision uniforme de toute la trajectoire-en mode continu ou en changement de mode.

Mots-clés : analyse structurelle, équations algébro-différentielles (DAE), systèmes multi-mode, modèles à structure variable, analyse non-standard 1 Introduction

Overall motivations

Multimode DAE systems constitute the mathematical framework supporting the physical modeling of systems possessing different modes. Each mode exhibits a different dynamics, captured by Differential Algebraic Equations (DAEs). Multimode DAE systems are the underlying framework of 'object-oriented' modeling languages such as Modelica, Amesim, or Simscape. Corresponding models can be represented as systems of guarded equations of the form if γ j (x i 's and derivatives) then f j (x i 's and derivatives) = 0 [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] where x i , i = 1, . . . , n denote the system variables, and, for j = 1, . . . , m, γ j (. . . ) is a Boolean combination of predicates guarding the differential or algebraic equation f j (. . . ) = 0. The meaning is that, if γ j has the value t (the constant true), then equation f j (. . . ) = 0 has to hold; otherwise, it is discarded. In particular, when all the predicates have the value t, one obtains a single-mode DAE, that is, a classical DAE defined by the set of equations f j (x i 's and derivatives) = 0 [START_REF] Alur | Hybrid automata: an algorithmic approach to the specification and verification of hybrid systems[END_REF] where i=1, . . . , n and j=1, . . . , m. When all f j 's have the special form x j -g j (x 1 , . . . , x n ), one recovers the Ordinary Differential Equations (ODE) x j =g j (x 1 , . . . , x n ). DAEs are a strict generalization of ODEs, where the so-called state variables x 1 , . . . , x n are implicitly related to their time derivatives x 1 , . . . , x n . This modeling framework is fully compositional, since systems of systems of equations of the form (1) are just systems of equations of the form [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF], with no restriction. DAE systems (having a single mode) are well understood. With comparison to ODE systems, a new difficulty arises with the notion of differentiation index, introduced in the late 1980's [START_REF] Petzold | Differential algebraic equations are not ODEs[END_REF][START_REF] Campbell | The index of general nonlinear DAEs[END_REF]. For simplicity, in this introduction, we discuss it for the particular case in which System (2) involves all derivatives x i but no higher-order derivative. The problem addressed by this notion of differentiation index relates to the Jacobian matrix J associated to System [START_REF] Alur | Hybrid automata: an algorithmic approach to the specification and verification of hybrid systems[END_REF], defined by J ij = ∂f j /∂x i . If this Jacobian matrix is invertible (requiring, in particular, m=n), then all the derivatives x i are uniquely determined as functions of the x i , so that System (2) is equivalent to an ODE system. If, however, matrix J is singular, then System (2) as such is no longer equivalent to an ODE system. This situation can occur, even for a square DAE system possessing a unique solution for any given initial conditions. This situation makes it difficult to design DAE solvers that work for any kind of DAE. Now, f j =0 in (2) implies d dt f j =0, revealing that such latent equations come implicitly with System [START_REF] Alur | Hybrid automata: an algorithmic approach to the specification and verification of hybrid systems[END_REF]. Adding latent equations to the DAE system does not change the system solutions, but it can bring additional constraints on highest-order derivatives of the original system states-on the other hand, this can introduce "spurious" derivatives (of higher order than in the original system), which are to be eliminated. Adding more latent equations while eliminating spurious derivatives eventually leads to highest-order derivatives of the original system states being uniquely determined as functions of the lower order state derivatives. The differentiation index [START_REF] Petzold | Differential algebraic equations are not ODEs[END_REF][START_REF] Campbell | The index of general nonlinear DAEs[END_REF] is the smallest integer k such that it is enough to differentiate with respect to time every equation f j =0 up to order at most k for the above situation to occur. DAE solvers exist for DAE of index 1, 2, or 3 [START_REF] Mattsson | Index reduction in Differential-Algebraic Equations using dummy derivatives[END_REF]. An efficient and popular approach, called the dummy derivatives method [START_REF] Mattsson | Index reduction in Differential-Algebraic Equations using dummy derivatives[END_REF], consists in reducing the index of the DAE system to 1, producing in addition a form in which the algebraic constraints are all preserved by the solver-this no longer holds if the index is reduced to zero, leading to an ODE. Structural analysis methods have been proposed [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF][START_REF] Pryce | A simple structural analysis method for DAEs[END_REF] which perform index reduction by only exploiting the bipartite graph associated to the different equations and variables of the system. Structural analyses scale up much better and provide valid results outside exceptional values for the coefficients arising in the system equations [START_REF] Campbell | The index of general nonlinear DAEs[END_REF]; structural analysis with dummy derivatives is implemented in most Modelica tools. Efforts are still ongoing to improve the efficiency and range of applicability of these methods, but one can say that DAE systems are now reasonably well understood.

Inria

In contrast, the handling of mode changes is much less mature. No structural analysis exists for mode changes. Due to discontinuities in system trajectories, the notion of differentiation index does not apply to multimode DAE systems. The notion of solutions of such systems is not even well understood, except for some subclasses of models, such as semi-linear systems [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF]. Still, modeling languages exist that support multimode DAE systems, e.g., Modelica, VHDL-AMS, and the proprietary languages Amesim and Simscape. However, the lack of mathematical understanding of mode changes can result in a spurious handling of some physically meaningful models. Indeed, with the exception of the recent work [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF], the class of "safe" models (well supported by the considered tool) is never characterized, thus leading to a "try-and-see" methodology. This situation motivates our work.

Establishing the mathematical foundations for compilers 1 of multimode DAE systems raises non-classical difficulties. To substantiate this claim, let us compare the following three subject matters:

1. Developing formal verification for a given class of hybrid systems [START_REF] Alur | Hybrid automata: an algorithmic approach to the specification and verification of hybrid systems[END_REF][START_REF] Platzer | Logical Foundations of Cyber-Physical Systems[END_REF]; 2. Developing existence/uniqueness results and discretization schemes, for a given class of (single-mode or multimode) DAE systems [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF];

3. Developing the mathematical foundations for multiphysics/multimode DAE systems modeling languages (our focus).

For subject 1, restrictions on the class are stated, under which the proposed algorithms are proved correct, and their complexity is analyzed. Similarly, for subject 2, assumptions are stated on which discretization schemes are proved correct, and convergence rates can be given. In both cases, assumptions are formulated and it is the responsibility of the user to check their validity when performing verification or applying discretization schemes.

In contrast, for subject 3, one cannot expect a compiler to check the existence/uniqueness of solutions of a given multimode DAE system. Worse, we cannot expect the user to check this as part of his/her model design activity. In fact, the compiler must handle any submitted model, and it has the responsibility for accepting or rejecting a model on the sole basis of syntactic or symbolic (but never numerical) analyses. This is a demanding task that we claim is not well addressed today.

Some illustration examples

To illustrate the above discussion, we now review two specific examples of multimode DAE systems.

An ideal clutch

This clutch is depicted in Figure 1. It is a simple, idealized clutch involving two rotating shafts where no motor or brake are connected. More precisely, we assume that this system is closed, with no interaction other than explicitly specified. We provide hereafter in (3) a model for it, complying with the general form [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. The dynamics of each shaft i 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. In addition, the model yields an ODE system when the clutch is released, and a DAE system of index 1 when the clutch is engaged, due to equation (e 3 ) being active in this mode. Note that the condition τ 1 + τ 2 = 0 is an invariant of the system in both modes, which reflects the preservation of the angular momentum due to the assumption that the system is closed. Model (3) is one possible specification. One could have instead put τ 1 + τ 2 = 0 outside the scope of any guard, thus making explicit that this is an invariant-this is actually the form adopted for the Modelica model of Figure 2; then, mode γ = f only states τ 1 = 0. Both models are equivalent. It turns out that our approach yields the same generated code for both source codes.

model ClutchBasic to parameter Real w01=1; to parameter Real w02=1.5; to parameter Real j1=1; to parameter Real j2=2; to parameter Real k1=0.01; to parameter Real k2=0.0125; to parameter Real t1=5; to parameter Real t2=7; to Real t(start=0, fixed=true); to Boolean g(start=false); to Real w1(start = w01, fixed=true); to Real w2(start= w02, fixed=true); to Real f1; Real f2; equation to der(t) = 1; to g = (t >= t1) and (t <= t2); to j1*der(w1) = -k1*w1 + f1; to j2*der(w2) = -k2*w2 + f2; to 0 = if g then w1-w2 else f1; to f1 + f2 = 0; end ClutchBasic; The clutch in Modelica Figure 2 details the Modelica model of the Ideal Clutch system, with two shafts interconnected by an idealized clutch. It is a faithful translation in the Modelica language of the mDAE (3), with the exception that the two differential equations have been linearized. Also, the trajectory of the input guard γ (here called g) has been fully defined: it takes the value t between instants t 1 and t 2 , and f elsewhere. The model is deemed to be structurally nonsingular by the two Modelica tools we had the opportunity to test: OpenModelica and Dymola. However, none of these tools generates correct simulation code from this model. Indeed, simulations fail precisely at the instant when the clutch switches from the uncoupled mode (g=false) to the coupled one (g=true). This is evidenced by a division by zero exception, as shown in Figures 3 and4.

The root cause of this exception is that none of these tools performs a multimode structural analysis. Instead, the structure of the model is assumed to be invariant, and the structural analysis implemented in these tools is a Dummy Derivatives method, which is proved to be correct only on single-mode DAE systems; it is possibly correct on multimode systems, but under the stringent assumption that the model structure (including, but not limited to, its differentiation index) is independent of the mode. The consequence is that the structural analysis methods implemented in these tools do not detect that the differentiation index jumps from 0 to 1 when the shafts are coupled, and that the structure is not invariant. The division by zero results from the pivoting of a linear system of equations that becomes singular when g becomes equal to true.

RR n°9334

The clutch in Mathematica DAE systems can be solved with the Mathematica "Advanced Hybrid & DAE" toolbox. 2 This library is also advertised as supporting multimode DAE systems.

Figure 5 shows the solution of the Ideal Clutch system, computed by Mathematica. It can be seen that the angular momentum of the mechanical system is not preserved, in contradiction with the physics. Also, a small change in the velocities before the change results in totally different values for restart, although the system is not chaotic by itself. A likely explanation is that Mathematica regards the mode switching as a consistent initialization problem, with one degree of freedom. Indeed, any choice of angular velocities satisfying the constraint ω 1 = ω 2 is a consistent initial state for the DAE in the engaged mode. The conclusion of these two experiments is that, clearly, some fundamental study is needed to ensure a correct handling of events of mode change by the Modelica tools. Smoothing the 'if then else' equation could help solving the above problem, but requires a delicate and definitely non-modular tuning, as it depends on the different time scales arising in the system. We believe that, as the tools reputedly support multimode DAE models, they should handle them correctly. Our objective For this example, our objective is to be able to compile System (3) and produce correct simulations for it, without any need for smoothing the mode change. This objective is addressed in Section 2. We sketch here a multimode extension of the popular example of the pendulum in Cartesian coordinates [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF], namely the planar Cup-and-Ball game, whose 3D version is illustrated by Figure 6. A ball, modeled by a point mass, is attached to one end of a rope, while the other end of the rope is fixed, to the origin of the plane in the model. The ball is subject to the unilateral constraint set by the rope, but moves freely while the distance between the ball and the origin is less than its actual length. Again, we assume that the system is closed and subject to no interaction other than specified. Using Cartesian coordinates x and y, the equations of the planar Cup-and-Ball game are the following:

A Cup-and-Ball game

           0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) 0 ≤ L 2 -(x 2 +y 2 ) (κ 1 ) 0 ≤ λ (κ 2 ) 0 = L 2 -(x 2 +y 2 ) × λ (κ 3 ) (4)
where the unknowns (also called dependent variables) are the position (x, y) of the ball in Cartesian coordinates and the rope tension λ.

The subsystem (κ 1 , κ 2 , κ 3 ) expresses that the distance of the ball from the origin is less than or equal to L, the tension is nonnegative, and one cannot have a distance less than L and a nonzero tension at the same time. This is known as a complementarity condition, written as 0 ≤ L 2 -(x 2 +y 2 ) ⊥ λ ≥ 0 in the nonsmooth systems literature [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF], and is an adequate modeling of ideal valves, diodes, and contact in mechanics. Constraints κ 1 and κ 2 are unilateral, which does not belong to our framework of guarded equations. Therefore, using a technique communicated to us by H. Elmqvist and M. Otter, we redefine the graph of this complementarity condition as a parametric curve, represented by the following three equations:

s = if γ then -λ else L 2 -(x 2 +y 2 ) 0 = if γ then L 2 -(x 2 +y 2 ) else λ γ = [s ≤ 0]
which allows us to rewrite the complementarity condition (κ 1 , κ 2 , κ 3 ) as follows:

                   0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) γ = [s ≤ 0] (k 0 ) if γ then 0 = L 2 -(x 2 +y 2 )
(k 1 ) and 0 = λ + s (k 2 ) if not γ then 0 = λ (k 3 ) and 0 = (L 2 -(x 2 +y 2 )) -s (k 4 ) [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF] Unsurprinsingly (considering the reasons for the wrong handling of the clutch example), the Modelica tools also fail to handle this model correctly.

Some new issues emerge from Example [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF]. First, subsystem (κ 1 , κ 2 , κ 3 ) of ( 4) leaves the impact law at mode change insufficiently specified: it could be fully elastic, fully inelastic, or inbetween. Second, (5) exhibits a logico-numerical fixpoint equation in s, which we regard as problematic. While the Jacobian matrix yields a regularity criterion for systems of smooth algebraic equations, logico-numerical systems of equations are unfriendly: there is no simple criterion for the existence and/or uniqueness of solutions. Note that a sensible modification of model [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF] would consist in replacing, in (k 0 ), s by its left-limit s -(t) = def lim sup u t s(u). 3Our objectives For this example, our objectives are to address the above difficulties at compile time, and to generate correct simulation code once structural problems are fixed. These objectives are addressed in Section 3.

A Westinghouse air brake

The second example is a simplified model of the railway air brake patented by George Westinghouse in 1868. A railway car is equipped with the pneumatic system whose schematics is shown in Figure 7. It has two pneumatic ports that are interconnected when several cars are coupled to form a train; in this simple example, only one car is considered, with the boundary conditions shown (control is an input pressure). As for the clutch example, this program is not correctly executed by the Modelica 3.3 engine. Using the smoothed version of a valve from the Modelica library of predefined components solves the problem. This, however, requires a careful tuning in accordance with the time scales of the whole system.

A simplified model corresponding to the schematics of Figure 7 is given in Figure 8, where we have taken the boundary conditions into account -the equations are simplified with reference to hydraulics, but the model structure (states and index) are correct. The model of the ideal valve is given by the subsystem (ν 1 , ν 2 , ν 3 ). Inequalities (ν 1 , ν 2 ) express nonnegativity conditions satisfied by the pressure difference and the flow through the valve. Equality (ν 3 ) expresses that, either the valve is open (no pressure drop), or it is closed (no air pass). This is known as a complementarity condition, written 0 ≤ p r -p t ⊥ f v ≥ 0

                                           0 = f b -f v -f cl + f l (m 1 ) 0 = f v -f ch -f l -f t (m 2 ) 0 = f b -F 1 .ϕ(p b -p r ) (f 1 ) 0 = F 1 .ϕ(p bn -p r ) (f 2 ) 0 = f l -F 2 .ϕ(p t -p r ) (f 3 ) 0 = f t -ρ.V
P0 .p t (r) 0 = ρ.S.(x .p r + (x -L).p r ) + P 0 .f cl (p 1 ) 0 = ρ.S.(x .S.p t + x.p t ) -P 0 .f ch (p 2 ) 0 = S.(p t -p r ) -b (l 1 ) 0 = K.x -b (l 2 ) 0 ≤ p r -p t (ν 1 ) 0 ≤ f v (ν 2 ) 0 = f v × (p r -p t ) (ν 3 )

Figure 8: Model corresponding to the schematics of Figure 7; p b , the control, and f bn determined by the boundary condition f bn =0, are known.

in the nonsmooth systems literature [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF], and is the adequate modeling of systems such as ideal valve, diode, and contact in mechanics. The two constraints (ν 1 , ν 2 ) are unilateral, which does not belong to our framework of guarded equations. Therefore, using a technique communicated to us by H. Elmqvist and M. Otter, we rewrite the complementarity condition (ν 1 , ν 2 , ν 3 ) by using an auxiliary variable s as follows:

           γ = [s ≤ 0] (v 0 ) if γ then 0 = p r -p t (v 1 ) and 0 = f v + s (v 2 ) if not γ then 0 = f v (v 3 ) and 0 = p r -p t -s (v 4 ) (6)
The system of Figure 8, where (ν 1 , ν 2 , ν 3 ) is replaced by (v 0 , . . . , v 4 ) in [START_REF] Benveniste | Structural Analysis of Multi-Mode DAE Systems[END_REF], has the form [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF]. In this model, the subsystem collecting equations (v 0 , . . . , v 4 ) constitutes a logico-numerical fixpoint equation, in that γ depends on s and the equations determining s are guarded by γ. Given that no Inria solver exists that reliably supports this kind of fixpoint equation, we think it should be rejected at compile time, with a proper diagnosis. Still, our model is physically meaningful.

Our objectives: For this example, our objectives are to be able to reject, with a proper diagnosis, the model of Figure 8 (using (6) for the valve), and to correct it so that the corrected model simulates as expected. These objectives are addressed in Section 4.

Related work

There exists a huge litterature devoted to the multi-physics DAEs (having a single mode); see, e.g., [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF] as a basic reference. Many references related to the Modelica language can also be found on the website of the Modelica organization. 4 We focus this review on works addressing multimode DAE systems with an emphasis on the compilation and code generation, for related models.

A large body of work was developed in the domain of multi-body systems with contacts where impulses can occur; see for instance [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Pfeiffer | On non-smooth multibody dynamics[END_REF] for an overview of this domain. The current research focuses on the handling of contacts between many bodies using time-stepping methods. Contact equations are usually expressed on the velocity or acceleration, and since the constraints on position level are not explicitly taken into account, a drift typically occurs. An exception is [START_REF] Schoeder | Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF], where constraints on position level are enforced. It is not obvious how the specialized multimode methods for multi-body systems can be generalized to any type of multimode DAE. For electrical circuits with idealized switches, like ideal diodes, impulses can also occur. Again, a large literature is available, concentrating mostly on specialized electrical circuits, such as piecewise-linear networks with idealized switches [START_REF] Heemels | On the dynamic analysis of piecewise-linear networks[END_REF]. However, it is not obvious how to generalize methods in this restricted area to general multimode DAEs.

The article by Mehrmann et al. [START_REF] Mehrmann | Hybrid systems of differential-algebraic equations -analysis and numerical solution[END_REF] contains interesting results regarding numerical techniques to detect chattering between modes. It however assumes that consistent reset values are explicitly given for each mode. Such an assumption does not hold in general, particularly for models derived from first principles of the physics; the clutch and cup-and-ball examples are good illustrations of this.

In the PhD thesis of Zimmer [START_REF] Zimmer | Equation-Based Modeling of Variable-Structure Systems[END_REF], variable-structure multi-domain DAEs are defined thanks to an ad hoc modeling language, and a runtime interpreter is used that processes the equations at runtime, when the structure and/or the index changes. Limitations of this work are that impulsive behavior is not supported and that the user has to explicity define the transfer of variable values from one mode to another, which is not practical for large models.

Describing variable-structure systems with causal state machines is discussed by Pepper et al. [START_REF] Pepper | A compositional semantics for modelica-style variable-structure modeling[END_REF]. Dynamically changing the structural analysis at runtime is also performed by Hoeger [START_REF] Höger | Dynamic structural analysis for daes[END_REF][START_REF] Höger | Elaborate control: variable-structure modeling from an operational perspective[END_REF]. In [START_REF] Höger | Dynamic structural analysis for daes[END_REF], the author proposes a dynamic execution of John Pryce's Σ-method [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF].

Elmqvist et al. [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF][START_REF] Mattsson | Multi-Mode DAE Systems with Varying Index[END_REF] propose a high-level description of multimode models as an extension to the synchronous Modelica 3.3 state machines, by using continuous-time state machines having continuous-time models as 'states'. Besides ODEs as used in hybrid automata, acausal DAE models with physical connectors can also be a 'state' of a state machine. Such a state machine is mapped into a form such that the resulting equations can be processed by standard symbolic algorithms supported by Modelica tools. The major restrictions of this approach are that mode changes with impulsive behavior are not supported and that not all types of multimode systems can be handled due to the static code generation.

Correct restart at mode changes is clearly a central issue in multi-mode DAE systems. Techniques from nonsmooth dynamical systems [START_REF] Vincent | Numerical Methods for Nonsmooth Dynamical Systems[END_REF] address this issue by performing a warm restart at mode changes using time-stepping discretization schemes: not only the numerical scheme does not interrupt at mode changes, but these mode changes are not even detected. This is very elegant and effective. However, this relies on the handling of complementarity conditions (see Section 1.2.2), thus, it does not apply to general multimode DAE systems but only to a certain subclass.

We observe that, in general, the issue of correct restart at mode changes is not addressed in the above mentioned literature. It is, however, considered in the few references to follow.

Benveniste et al. [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF][START_REF] Benveniste | Structural Analysis of Multi-Mode DAE Systems[END_REF] considered this issue, as well as the problem of varying structure and index, from a fundamental point of view, by relying on nonstandard analysis to capture continuous-time dynamics and mode change events in a unified framework. A first structural analysis algorithm was presented in [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF], by significantly modifying the original Pantelides algorithm [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF]. A proof-of-concept mockup named SunDAE was developed implementing this approach. This first attempt suffers from some issues: the proposed structural analysis does not boil down trivially to the Pantelides algorithm in the case of single-mode systems; it involves nondeterministic decisions, an unwanted feature for the mathematical foundation of compilers; and its mathematical study is incomplete.

The work of Stephan Trenn is an important contribution to the subject. In his PhD thesis [START_REF] Trenn | Distributional Differential Algebraic Equations[END_REF] and his article [START_REF] Trenn | Regularity of distributional differential algebraic equations[END_REF], he points out the difficulty in defining piecewise smooth distributions: there is no intrinsic definition for their 'Dirac' part, as several definitions comply with the requirements for it. This indicates that distributions are not the ultimate answer to deal with impulsive variables in multimode DAE systems. Still, Trenn was able in [START_REF] Liberzon | Switched nonlinear differential algebraic equations: Solution theory, lyapunov functions, and stability[END_REF] to define complete solutions for a class of switched DAE systems in which each mode is in quasi-linear form: switching conditions are time-based, not state-based. A main difference with our approach is that Trenn's solution is complete for all variables and relies on piecewise smooth distributions, whereas our works on restart conditions consist in eliminating impulsive variables while computing restart values for state variables. (Beyond the existence of solutions, the work of Trenn et al. contains additional interesting results that are not relevant to our study.)

An important step forward was done in [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF]. The interesting subclass of multimode DAE systems that is now called semi-linear (see Section 11.2) was first identified. Semi-linear multimode DAE systems can exhibit impulsive variables at mode changes. They generalize the 'semi-linear systems' proposed by Trenn in the sense that switching conditions are no longer restricted to time-based ones, instead including state-based switching conditions. The analysis and discretization schemes proposed in [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF] are mathematically sound. Building on this work, Martin Otter has developed the ModiaMath tool for semi-linear multimode DAE systems. Hence, it makes sense to compare the schemes proposed in [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF] to the ones we develop in this paper, for general mDAE systems. It turns out that our general approach coincides with the schemes proposed in [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF] when applied to the subclass of semi-linear systems, see Section 11.2. Our present work thus extends and significantly improves [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF].

Contributions

In this work, we develop a systematic approach that addresses, as particular cases, the objectives stated for the examples of Sections 1.2.1 and 1.2.2. Its main contributions can be detailed as follows.

Structural analysis of mode changes using nonstandard analysis As our main contribution, we extend the notion of structural analysis to mode change events. More precisely, we develop a structural analysis for multimode DAE systems that is valid at any time, that is, for both continuous dynamics and mode changes.

Doing so raises a number of difficulties. First, mode changes occur under various kinds of dynamics; there can be isolated events or finite cascades thereof, but Zeno situations can also occur, i.e., events of mode changes can accumulate and form an infinite sequence bounded by some finite instant. As a result, sliding modes can emerge out of mode changes, and the modeling and simulation tools can only identify such situations at runtime. Second, impulses may occur at mode changes for certain variables, which again is most of the time discovered by tools at runtime. Compilation, however, is a task performed prior to generating simulation code; hence, it encounters serious difficulties in both the nature of time at mode changes and the domain of variables, due to impulsive behaviors.

The main objective of structural analysis, akin to compilation in general, is precisely to keep Inria away from numerical considerations, by focusing only on the structure (or the syntax) of the system of equations. This is captured by the bipartite graph associated to the considered DAE system (also referred to as 'incidence graph' or 'incidence matrix' or 'Σ-matrix' depending on the context and authors). The basic principle supporting structural analysis consists in: (1) abstracting an equation f (x, y, z) = 0 as the bipartite graph having f, x, y, z as vertices and (f, x), (f, y), (f, z) as edges, and a system of equations as the union of the above bipartite graphs; and (2) exploiting the resulting graph to generate efficient simulation code. This means focusing on syntax and symbols, while abstracting away numerical aspects. We thus wish to symbolically manipulate both normal and impulsive values for variables, as well as both continuous time and events (in finite or infinite cascades), and possibly stranger structures for time.

We see no vehicle for supporting all of this, other than nonstandard analysis. Nonstandard analysis was proposed by Abraham Robinson in the 1960s to allow the explicit manipulation of 'infinitesimals' and 'infinities' in analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Cutland | Nonstandard analysis and its applications[END_REF]. Impulsive values become normal citizens in nonstandard analysis. Continuous time and dynamics can be discretized using infinitesimal time steps, e.g., by setting x (t) ≈ x(t+∂)-x(t) ∂ with ∂ infinitesimal, thus providing perfect fidelity up to infinitesimal errors. This discretized time can then support the various kinds of time needed when modeling mode changes. Since techniques from structural analysis translate almost verbatim to discrete-time dynamics, having dynamical systems indexed by discrete (nonstandard) time in a uniform setting paves the way toward structural analysis for multimode DAE systems in their generality. Using this approach, we are able to extend the structural analysis, from single-mode to multimode DAE systems, by encompassing both modes and mode changes.

Generating code for the restart at mode changes The structural analysis developed hereafter works in the nonstandard analysis domain. As such, it cannot be executed on any real-life computer. A further stage is needed to produce executable code, by mapping the nonstandard execution schemes to ordinary world (standard) schemes; this is called standardization. This allows us to produce mathematically justified schemes for the restart at mode changes.

Rejecting and accepting programs on a clear basis Our structural analysis is precise enough to properly identify models that are structurally over-or under-specified at mode change events, or models that contain a fixpoint equation involving a subset of variables and some guard. We reject such models, with a proper diagnosis explaining the rejection. In turn, mode-dependent index/state/dynamics are not reasons for rejection per se. Our compilation technique is able to handle such cases.

The paper is organized as follows This report is an extended version of a paper under submission. With reference to this submission, it contains additional material, which the reader may skip for a first reading. In the following list of sections, we indicate in blue these additional sections.

The clutch example is informally developed in Section 2; we show how to map this model to the nonstandard domain, how the structural analysis works (in particular, the handling of mode changes is emphazised), and how standardization is performed to generate actual simulation code. The Cup-and-Ball example is developed in Section 3, following the same guidelines; based on the case of elastic impact, it brings forth the need for considering transient modes, i.e., modes that last for zero time. Section 4 develops the Westinghouse air brake example; we address the objectives stated in the introduction and we advocate for the interest of adding assertions, restricting the allowed sequences of modes in the system. The insight gained from studying these examples leads to the approach described in Section 5.

The extension of structural analysis to multimode DAE systems is then developed. In Section 6, the basics of structural analysis for algebraic systems of equations are recalled. Then, we review J. Pryce's Σ-method for the structural analysis of DAE systems. We extend the Σ-method to the multimode case, including the handling of mode change events, in Section 7, for the subclass of systems involving only long modes-which excludes the Cup-and-Ball example with elastic impact.

Motivated by the latter case, we identify the need for handling transient modes, in which the system spends zero time; this is addressed in Section 7.4. In Section 8 we widen our scope by adding to our framework, in addition to derivatives, the left-and right-limits of a signal, which is a useful primitive operator provided by modeling languages such as Modelica or Simulink.

Section 9 recalls the needed background on nonstandard analysis, mostly related to differential calculus. In Section 10, we complement this background with additional results related to structural analysis and impulse analysis, and we apply this machinery to the task of standardization.

Three important results are developed in Section 11.

• In Section 11.1, we prove that modifying the nonstandard expansion of the derivative (there are infinitely many possibilities) does not change the final generated code, thus showing that our approach is intrinsic.

• We prove in Section 11.2 that our approach does compute correct solutions for semi-linear multimode DAE systems, a nontrivial subclass possibly involving impulsive behaviors and containing, as a subclass, multi-body mechanics (see [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF]). We would be happy to prove that the simulation code we generate does compute the solution of any multimode DAE system; this, however, is beyond what is doable since no notion of solution is known for multimode DAE systems in general.

• We propose in Section 11.3 a numerical scheme for approximating restart values for nonimpulsive variables at mode changes. This alleviates the need for performing standardization when generating code in practice; the proof of this result, however, uses standardization as a mathematical tool.

Finally, we draw in Section 12 how our approach can be implemented in a tool, and we illustrate this on the RLDC2 circuit example. This example, provided to us by Sven-Erik Mattsson, is problematic for the existing Modelica tools, as it has mode-dependent index and structure. The structural analysis of the corresponding DAE system in all its modes is successfully handled by our IsamDAE tool [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF]. This tool, currently under development, will support multimode DAE structural analysis following our theory. Its current version does not support impulsive variables at mode changes-the RLDC2 circuit example, however, does not exhibit this difficulty.

In Section 7, we assumed that the effect of guards is delayed by one nonstandard time shift. This is a necessary condition for being able to support both long and transient modes. In Apprendix A, we propose an alternative execution scheme that is constructive in that it can handle models in which the effect of guards is not delayed by one nonstandard time shift. In turn, only models with long modes are supported.

The ideal clutch

In this section, we introduce our approach by discussing in detail the ideal clutch example presented in Section 1.2.1. Its model was given in [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF]. We first analyze separately the model for each mode of the clutch (Section 2.1). Then, we discuss the difficulties arising when handling mode changes (Section 2.2). Finally, we propose a global comprehensive analysis in Sections 2.3 -2.5. For convenience, we present an intuitive and informal introduction to nonstandard analysis in Section 2.3.

Separate Analysis of Each Mode

In the released mode, when γ is false in System (3), the two shafts are independent and one obtains the following two independent ODEs for ω 1 and ω 2 :

ω 1 = f 1 (ω 1 , τ 1 ) (e 1 ) ω 2 = f 2 (ω 2 , τ 2 ) (e 2 ) τ 1 = 0 (e 5 ) τ 2 = 0 (e 6 ) (7)

Inria

In the engaged mode, however, γ holds true, and both the two velocities and the two torques are algebraically related:

ω 1 = f 1 (ω 1 , τ 1 ) (e 1 ) ω 2 = f 2 (ω 2 , τ 2 ) (e 2 ) ω 1 -ω 2 = 0 (e 3 ) τ 1 + τ 2 = 0 (e 4 ) (8)
Equation (e 4 ), that relates the two torques, is not an issue by itself, as it can be rewritten as τ 2 = -τ 1 and used to eliminate τ 2 . Equation (e 3 ), however, relates the two velocities ω 1 and ω 2 that are otherwise subject to the ODE system (e 1 , e 2 ), which makes System (8) a DAE instead of an ODE.

In System [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF], the torques τ 1 and τ 2 are not differentiated: we call them algebraic variables. The derivatives of the velocities ω 1 and ω 2 appear: they are state variables. The ω 1 , ω 2 , τ 1 , τ 2 are called the leading variables of System [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF]. Now, suppose one is able to uniquely determine the leading variables ω 1 , ω 2 , τ 1 , τ 2 given consistent values for the state variables ω 1 , ω 2 , i.e., values satisfying (e 3 ). Then, by using an ODE solver, one could perform an integration step and update the current velocities ω 1 , ω 2 using the computed values for their derivatives ω 1 , ω 2 . In this case, we say that the considered DAE is an "extended ODE" [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF].

It turns out that this condition does not hold for System [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF] as is. To intuitively explain the issue, we move to discrete time by applying an explicit first-order Euler scheme with constant step size δ > 0:

ω • 1 = ω 1 + δ • f 1 (ω 1 , τ 1 ) (e δ 1 ) ω • 2 = ω 2 + δ • f 2 (ω 2 , τ 2 ) (e δ 2 )
ω 1 -ω 2 = 0 (e 3 ) τ 1 + τ 2 = 0 (e 4 ) (9)
where ω • (t) = def ω(t+δ) denotes the forward time shift operator by an amount of δ. Suppose we are given consistent values ω 1 = ω 2 and we wish to use System [START_REF] Benveniste | Compositionality in dataflow synchronous languages: Specification and distributed code generation[END_REF], seen as a system of algebraic equations, to determine the value of the dependent variables (i.e., the unknowns)

τ 1 , τ 2 , ω • 1 , ω • 2 .
This attempt fails since we have only three equations e δ 1 , e δ 2 and e 4 to determine four unknowns τ 1 , τ 2 , ω • 1 and ω • 2 ; indeed, the omitted equation (e 3 ) does not involve any dependent variable. However, since System (9) is time invariant, and assuming that the system remains in the engaged mode for at least δ seconds, there exists an additional latent equation on the set of variables, namely

ω • 1 -ω • 2 = 0 (e • 3 ) (10)
obtained by shifting (e 3 ) forward. Now, replacing (e 3 ) with (e • 3 ) in System (9) yields a system with four equations and four dependent variables; moreover, this system is, in fact, nonsingular in a generic sense. One can now use System (9) along with equation (e •

3 ) to get an execution scheme for the engaged mode of the clutch. This is shown in Exec. Sch. 1 below.

Execution Scheme 1 System (9)+Eq. [START_REF] Benveniste | The synchronous languages 12 years later[END_REF]. Require: consistent ω 1 and ω 2 , i.e., satisfying (e 3 ).

1: Solve {e δ 1 , e δ 2 , e • 3 , e 4 } for (ω • 1 , ω • 2 , τ 1 , τ 2 ) 4 equations, 4 unknowns 2: (ω 1 , ω 2 ) ← (ω • 1 , ω • 2 ) update the states (ω 1 , ω 2 ) 3: Tick move to next discrete step
Since the new values of the state variables satisfy [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] by construction, the consistency condition is met at the next iteration step (should the system remain in the same mode).

Comment 1 (structural nonsingularity) The implicit assumption behind Line 1 in Exec. Sch. 1 is that solving {e δ 1 , e δ 2 , e • 3 , e 4 } always returns a unique set of values. In our example, this is true in a 'generic' or 'structural' sense, meaning that it holds for all but exceptional values for the parameters in the considered equations (through the unspecified functions f 1 and f 2 ). We will repeatedly use structural nonsingularity (also called structural regularity) in the sequel. A formalization of structural reasoning and analysis is developed in Section 6.

Observe that the same analysis could be applied to the original continuous-time dynamics from System (8) by augmenting it with the following latent equation:

ω 1 -ω 2 = 0 (e 3 ) (11)
obtained by differentiating (e 3 ); as a matter of fact, since (e 3 ) holds at any instant, (e 3 ) follows as long as the solution is smooth enough for the derivatives ω 1 and ω 2 to be defined. The resulting execution scheme is given in Exec. Sch. 2, that parallels Exec. Sch. 1.

Execution Scheme 2 System (8)+Eq. [START_REF] Berge | The theory of graphs and its applications[END_REF].

Require: consistent ω 1 and ω 2 , i.e., satisfying (e 3 ). update the states (ω 1 , ω 2 ) 3: Tick move to next discretization step Line 1 is identical for the two schemes and is assumed to give a unique solution, generically; it fails if one omits the latent equation (e 3 ). Then, although getting the next values for the ω 1 and ω 2 is easy in Exec. Sch. 1, the situation changes in Exec. Sch. 2: the derivatives (ω 1 , ω 2 ) are first evaluated, and then an ODE solver (here denoted by ODESolve) is used to update the state (ω 1 , ω 2 ) for the next discretization step. Note that, when considering an exact mathematical solution, if ω 1 -ω 2 = 0 initially holds and ω 1 -ω 2 = 0, then the linear constraint (e 3 ) will be satisfied at any positive time.

The replacement of (e 3 ) by (e 3 ), which resulted in Exec. Sch. 2, is known in the literature as index reduction [START_REF] Mattsson | Index reduction in Differential-Algebraic Equations using dummy derivatives[END_REF]. It requires adding the (smallest set of) latent equations needed for the system in Line 1 of the execution scheme to become solvable and deterministic. The number of successive time differentiations needed for obtaining all the latent equations is called the differentiation index [START_REF] Campbell | The index of general nonlinear DAEs[END_REF], which we simply refer to as the index in the sequel. 5 Finally, observe that both execution schemes rely on an algebraic equation system solver to evaluate the leading variables as a function of state variables.

We conclude this part by briefly discussing the initialization problem. Unlike for ODE systems, initialization is far from trivial for DAE systems. For the clutch example, if one considers System (8) as a standalone DAE, the initialization is performed as indicated in Exec. Sch. 3.

Execution Scheme 3 Initialization of System (8)+Eq. [START_REF] Benveniste | The synchronous languages 12 years later[END_REF].

1: (ω 1 , ω 2 ) ← Solve{e 3 } for (ω 1 , ω 2 )
2 unknowns for only 1 equation

The system for solving possesses two unknowns for only one equation, which leaves us with one degree of freedom; in mathematical terms, the set of all initial values for the 2-tuple of variables is a manifold of dimension 1, generically or structurally. For example, one may freely fix the initial common rotation speed so that (e 3 ) is satisfied.

Mode Transitions

In an attempt to reduce the full clutch model to the analysis of the DAE of each mode, one hopes that the handling of a mode change boils down to applying the initialization given in Exec. Sch. 3. If one was to handle restarts at mode changes like initializations, it would make the clutch system nondeterministic, due to the extra degree of freedom of Exec. Sch. 3. In contrast, the physics tells us that the state of the system is entirely determined when the clutch gets engaged. This comforts the intuition that resets at mode changes strongly differ from mere initializations.

However, inferring by hand the reset values for rotation velocities when the clutch gets engaged is definitely non-trivial. Furthermore, these values depend on the whole system model, so that the task of determining them becomes complex if external components are added. It is therefore highly desirable, for this example, to let the compiler infer these reset values from model [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF].

This problem involves a mix of continuous-time and discrete-time dynamics, the latter consisting in one or several computation steps to restart the dynamics of the next mode. As a result, it will be convenient to cast everything in discrete time thanks to nonstandard analysis. We develop this next.

Nonstandard Semantics

Nonstandard analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF][START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF] extends the set R of real numbers into a superset R of hyperreals (also called nonstandard reals) that includes infinite sets of infinitely large numbers and infinitely small numbers. The key properties of hyperreals, needed for the informal discussion of the clutch example, are the following [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF]: 6• There exist infinitesimals, defined as hyperreals that are smaller in absolute value than any real number: an infinitesimal ∂ ∈ R is such that |∂| < a for any positive a ∈ R. For x, y two hyperreals, write x ≈ y if x -y is an infinitesimal.

• All relations, operators, and propositional formulas that are valid over R are also valid over R; for example, R is a totally ordered set. The arithmetic operations +, ×, etc. can be lifted to R. We say that a hyperreal x is finite if there exists some standard finite positive real number a such that |x| < a.

• For every finite hyperreal x ∈ R, there exists a unique standard real number st(x) ∈ R such that st(x) ≈ x, and st(x) is called the standard part (or standardization) of x. As we shall see, standardizing more complex objects, such as functions or systems of equations, requires some care.

• Every real function lifts in a systematic way to a hyperreal function (regardless of its continuity properties). 7 This allows us to write f (x) where f is a real function and x is a nonstandard number.

• Let t → x(t), t ∈ R be an R-valued (standard) signal. Then:

x is continuous at instant t ∈ R if and only if, for any infinitesimal

∂ ∈ R, one has x(t + ∂) ≈ x(t); (12)
x is differentiable at instant t ∈ R if and only if there exists a ∈ R such that, for any infinitesimal ∂ ∈ R, x(t+∂)-x(t) ∂ ≈ a. In this case, a = x (t). [START_REF] Bliudze | Un cadre formel pour l'étude des systèmes industriels complexes: un exemple basé sur l'infrastructure de l'UMTS[END_REF] We can then consider, for a given positive infinitesimal ∂, the following time index set T ⊆ R:

T = 0, ∂, 2∂, 3∂, • • • = {n∂ | n ∈ N} ( 14
)
where N denotes the set of hyperintegers, consisting of all natural integers augmented with additional infinite numbers called nonstandard. For convenience, the elements of T will be generically denoted by the symbol t, the explicit expression n∂, or even simply by t when the possible confusion with (usual) real time does not matter. The important features of T are:8

• Any finite positive real time t ∈ R, t ≥ 0 is infinitesimally close to some element of T (informally, T covers the set R + of nonnegative real numbers and can be used to index continuous-time dynamics); and

• T is 'discrete': every instant n∂ has a predecessor (n-1)∂ and a successor (n+1)∂, except 0 that has no predecessor.

Definition 2 (Forward and Backward Shifts) Let x be a nonstandard signal indexed by T. We define the forward shifted signal x • through x • (n∂) = def x((n+1)∂) and, for f (X) a function of the tuple X of signals, we set (f (X)) • = def f (X • ) where the forward shift X → X • applies pointwise to all the components of the tuple. It will also be convenient to consider the backward shift • x defined by • x((n+1)∂) = def x(n∂), implying that an initial value for • x(0) must be provided.

For example, f • (x, y)(t) = def f (x • , y • )(t) = f (x(t+∂), y(t+∂)).
Observe that the definition of the forward shift depends on the infinitesimal ∂. Thanks to (13), using forward shifts allows us to represent, up to an infinitesimal, the derivative x of a signal by its first-order explicit Euler approximation 1 ∂ (x • -x). Solutions of multimode DAE systems may, however, be non-differentiable and even noncontinuous at events of mode change. To give a meaning to x at any instant, we decide to define it everywhere as the nonstandard first-order Euler increment

x = def 1 ∂ (x • -x) . ( 15
)
Definition 3 Substituting, in a multi-mode DAE system, every occurrence of x (t) by its expansion (15) yields a difference algebraic equation (dAE) system9 that we call its nonstandard semantics.

For instance, the nonstandard semantics of System (3) is the following:

                 ω • 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 )

The state variables are ω 1 , ω 2 whereas the leading variables are now τ 1 , τ 2 , ω • 1 , ω • 2 , in both modes γ = f and γ = t. Reproducing the reasoning of Section 2.1, one obtains, for each mode, a discrete system very much like the explicit Euler scheme of Section 2.1, except that the step size is now infinitesimal and that the variables are all nonstandard. The added value of System [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF] with respect to the explicit Euler scheme of Section 2.1 is twofold: first, it is exact up to infinitesimals within each mode; second, it will allow us to carefully analyze what happens at events of modes change.

Nonstandard structural analysis

We now develop a structural analysis for the nonstandard multimode dAE system [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF]. As a preliminary step, we perform the structural analysis in each mode and combine the two dynamics with the due guards (the added latent equation is 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 ) (17)
System [START_REF] Broman | Modelyze: a gradually typed host language for embedding equation-based modeling languages[END_REF] coincides with the known structural analysis in the interior of each of the two modes γ = t and γ = f. What about the instants of mode change?

Mode change γ : t → f At the considered instant, we have

             ω 1 -• ω 1 ∂ = f 1 ( • ω 1 , • τ 1 ) ω 2 -• ω 2 ∂ = f 2 ( • ω 2 , • τ 2 ) • ω 1 -• ω 2 = 0 ( • e 3 ) ω 1 -ω 2 = 0 • τ 1 + • τ 2 = 0 current            ω • 1 -ω 1 ∂ = f 1 (ω 1 , τ 1 ) ω • 2 -ω 2 ∂ = f 2 (ω 2 , τ 2 ) τ 1 = 0 τ 2 = 0 (18)
We assume that values are given for the state variables at the previous instant, namely • ω 1 and • ω 2 , and we regard System (18) as an algebraic system of equations with dependent variables • τ i , ω i ; τ i , ω • i for i = 1, 2, i.e., the leading variables of System (17) at the previous and current instants. Note that equation ( • e 3 ) does not involve any dependent variable: it is actually a fact that is inherited from the instant t -2∂, where t is the current instant; we can thus ignore this equation in System [START_REF] Bryant | Graph-based algorithms for Boolean function manipulation[END_REF]. Seen this way, System [START_REF] Bryant | Graph-based algorithms for Boolean function manipulation[END_REF] is structurally nonsingular and we can solve it as two successive blocks, corresponding to the previous instant followed by the current instant: the latter is our desired code for the mode change. This yields the following code for the mode change γ : t → f:

               ω 1 , ω 2 set by previous instant ω • 1 -ω 1 ∂ = f 1 (ω 1 , τ 1 ) ω • 2 -ω 2 ∂ = f 2 (ω 2 , τ 2 ) τ 1 = 0 τ 2 = 0 (19)
The above reasoning may seem an overshoot for this mode change γ : t → f, since one could simply say: let the previous part of System [START_REF] Bryant | Graph-based algorithms for Boolean function manipulation[END_REF] fix ω 1 , ω 2 and, then, let the current part of System (18) determine τ i , ω • i for i = 1, 2. However, the same line of reasoning will allow us to address the other, more difficult, mode change as well.

Mode change γ : f → t At the considered instant, we have • γ = f and γ = t. We proceed as for the other case. Unfolding System (17) 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 (20)
Once more, we regard System (20) as an algebraic system of equations with dependent variables

• τ i , ω i ; τ i , ω • i for i = 1, 2, i.e.
, the leading variables of System (17) at the previous and current instants.

In contrast to the previous mode change, System [START_REF] Caillaud | Implicit Structural Analysis of Multimode DAE Systems[END_REF] is no longer structurally nonsingular since the following subsystem possesses five equations and only four dependent variables

ω 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 ) (21)
This subsystem exhibits a conflict. We propose to resolve it by applying the following causality principle:

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

Applying Principle 1 leads to removing, from subsystem [START_REF] Campbell | The index of general nonlinear DAEs[END_REF], the conflicting equation (e 3 ). This yields the following code for the mode change γ : f → t:

               ω 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 (22)
Note that the consistency equation (e 3 ) : ω 1 -ω 2 = 0 has been removed from System [START_REF] Cartan | Formes Différentielles. Collection Méthodes[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 nonstandard instant the satisfaction of some of the constraints in force in the new mode γ = t. Since our time step ∂ is infinitesimal, this takes zero standard time.

The corresponding nonstandard execution scheme is summarized in Exec. Sch. 4. We use variable ∆ to encode the context, that is, the set of equations known to be satisfied by the state variables as a result of having executed the previous instant. For instance, e 3 ∈ ∆ corresponds to the mode change γ : f → t, expressing that (e 3 ) is not provably true if the guard was false at the previous instant. At each tick, that is, every time the scheme moves to a new instant, the context gets updated to account for the equations satisfied by the new state. The procedure 'Solve' solves the (algebraic) system to determine the values of the leading variables.

Execution Scheme 4 for Nonstandard System [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF]. Comments are written in blue. Require: ω 1 and ω 2 .

1: if γ (engaged mode) then

2:
if e 3 / ∈ ∆ (mode change) then 3:

(ω • 1 , ω • 2 ) ← Solve e ∂ 1 , e ∂ 2 , e • 3 , e 4
4:

Tick: ∆ ← ∆ ∪ {e 3 } 5: else 6: (τ 1 , τ 2 , ω • 1 , ω • 2 ) ← Solve e ∂ 1 ,
(τ 1 , τ 2 , ω • 1 , ω • 2 ) ← Solve e ∂ 1 ,

Standardization

Exec. Sch. 4 cannot be run in its present form, since it involves nonstandard reals. To recover executable code over the real numbers, a supplementary standardization step is needed. Recall that any finite nonstandard real x has a unique standard part st(x) ∈ R such that x ≈ st(x). The standardization procedure aims at recovering the standard parts of the leading variables from their nonstandard version. We distinguish two cases: continuous evolutions within each mode, assuming that the sojourn time in each mode has (standard) positive duration, and discrete evolutions at events of mode change.

Within continuous modes

Let x : t → x(t), t ∈ [s, p), denote the real continuous and differentiable solution in a given mode (assuming it exists). In the sequel, we use the notation o(ξ) to denote any (unspecified) nonstandard real number ζ such that ζ/ξ is infinitesimal. Using [START_REF] Bliudze | Un cadre formel pour l'étude des systèmes industriels complexes: un exemple basé sur l'infrastructure de l'UMTS[END_REF] and this notation, equations (e ∂ 1 ) and (e ∂ 2 ) rewrite as ω i = f i (ω i , τ i ) + o(1), which can be shown to standardize as the differential equations (e 1 ) and (e 2 ), respectively. This suffices to recover the dynamics [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF] in the released mode γ = f.

In the engaged mode γ = t, however, we need to handle (e • 3 ) :

ω • 1 -ω • 2 = 0. The nonstandard characterization of derivatives writes ω • i = ω i + ∂.ω i + o(∂)
, where ω i and ω i are both standard. Since ω 1 -ω 2 = 0 and ω

• 1 -ω • 2 = 0 both hold, we inherit ω 1 -ω 2 = o(1)
, which implies ω 1 -ω 2 = 0 since the ω i are standard. We thus recover the dynamics (8) augmented with the latent equation [START_REF] Berge | The theory of graphs and its applications[END_REF] in the engaged mode γ = t.

To summarize, within continuous modes, we perform structural analysis as usual for DAE systems in continuous time: nothing new happens.

At the instants of mode change

Suppose we have an event of mode change at time t, meaning that γ(t) = γ(t -∂). Our aim is to use the next values ω • i (t) = ω i (t + ∂) to initialize the state variables for the next mode with the value ω + i (t) = ω • i (t). Thus, unlike for continuous time dynamics, the occurrence of the infinitesimal time step ∂ in the expression ω i (t + ∂) is not an issue: it just points to the supplementary instant at which restart is performed.

However, the equations defining ω • i as functions of the ω i 's involve the hyperreal ∂ in space: standardization must be performed to get rid of it.

Transition γ : t → f It is easy, as the new mode has an ODE dynamics whose state variables are initialized with corresponding exit values when leaving the previous mode.

ω • i -ω i ∂ = f i (ω i , τ i ), i = 1, 2. ( 23
)
At this point, we must identify possible impulsive variables, which is the objective of the following impulse analysis.

Impulse analysis. Before engaging the clutch, we must assume

ω 1 -ω 2 = 0. Since ω • 1 -ω • 2 = 0 holds, (ω • 1 -ω • 2 )-(ω 1 -ω 2 ) ∂ = f 1 (ω 1 , τ 1 ) -f 2 (ω 2 , τ 2
) cannot be a finite nonstandard real because, if it was, then the function t → ω 1 (t) -ω 2 (t) would be right-continuous, in contradiction with the above assumption. Hence, the nonstandard real f 1 (ω 1 , τ 1 ) -f 2 (ω 2 , τ 2 ) is necessarily infinite. However, we assumed continuous functions f i and finite state (ω 1 , ω 2 ). Thus, one of the torques τ i must be infinite at t, and because of equation (e 4 ) : τ 1 + τ 2 = 0, both torques are in fact infinite, i.e., are impulsive.

Eliminating impulsive variables. It remains to compute the (standard) restart values for the state variables. To this end, we assume that the f i 's are linear in the torques, i.e., each f i has the following form:

f i (ω i , τ i ) = a i (ω i ) + b i (ω i )τ i , (24)
where b 1 and b 2 are the inverse moments of inertia of the rotating masses and a 1 and a 2 are damping factors divided by the corresponding moments of inertia. This yields the following system of equations, to be solved for ω • 1 , ω • 2 , τ 1 , τ 2 at the instant when γ switches from f to t:

       ω • 1 = ω 1 + ∂.(a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ) (e ∂ 1 ) ω • 2 = ω 2 + ∂.(a 2 (ω 2 ) + b 2 (ω 2 )τ 2 ) (e ∂ 2 ) ω • 1 -ω • 2 = 0 (e • 3 ) τ 1 + τ 2 = 0 (e 4 ) (25)
It is tempting to standardize System (25) by simply setting ∂ = 0 in it. Unfortunately, doing this leaves us with a structurally singular system, since the two torques are then involved in only one equation.

Comment 4 (standardizing functions vs. standardizing equations) To get a better insight into this problem, let us replace for a while the infinitesimal ∂ by a 'small' (standard) time step δ and rewrite System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF] as

F (ω • 1 , ω • 2 , τ 1 , τ 2 , ω 1 , ω 1 , δ) = 0
, where F =0 stacks the four equations of (25) into one vector equation. Then, the function [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF], the standardization of

δ → F (ω • 1 , ω • 2 , τ 1 , τ 2 , ω 1 , ω 1 , δ) is continuous in δ at δ=0. Consequently, by
F (ω • 1 , ω • 2 , τ 1 , τ 2 , ω 1 , ω 1 , ∂), seen as a function, is obtained by setting ∂=0 in it.
This, however, becomes incorrect if, instead of the function F , we are interested in the equation F =0. More precisely, in our example, let Y be the vector that collects the states ω 1 , ω 2 , and X be the vector that collects the dependent variables ω • 1 , ω • 2 , τ 1 , τ 2 ; then, the function F defined by System (25) is linear in δ, i.e., it has the form Y → X, where X is solution of the algebraic equation

[A -δ×B(Y )]X = C(δ, Y ), in which A =     1 0 0 0 0 1 0 0 1 -1 0 0 0 0 1 1     ; B = B(Y ) =     0 0 b 1 (ω 1 ) 0 0 0 0 b 2 (ω 2 ) 0 0 0 0 0 0 0 0    
and vector C(δ, Y ) is appropriately fixed; from now on, we omit Y in the writing. For every δ > 0, the matrix A-δ×B is invertible and we have X=(A-δ×B) -1 C.

If A was nonsingular, we could simply write X = (I-δ×D) -1 A -1 C, where D = A -1 B. Expanding the inverse around δ = 0 would yield (I-δ×D) -1 = I+δ×D+o(δ), implying X = (I+δ×D)A -1 C+o(δ). Therefore, the solution X(∂) of equation (A-∂×B(Y ))X=C(∂, Y ) would standardize as the solution X of the standard equation AX = C(0, Y ). All of this, however, requires that A be nonsingular, which is wrong in our case. 10When A is singular (as in our case), the solution consists in eliminating the impulsive variables from System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF], namely, the two torques. This is performed by symbolic pivoting:

1. Using (e 4 ) yields -τ 2 = τ 1 = def τ ; 2. Premultiplying the system of equations

ω • 1 = ω 1 + ∂.(a 1 (ω 1 ) + b 1 (ω 1 )τ ) (e ∂ 1 ) ω • 2 = ω 2 + ∂.(a 2 (ω 2 ) -b 2 (ω 2 )τ ) (e ∂ 2 )
by the row matrix b 2 (ω 2 ) b 1 (ω 1 ) yields the equation

b 2 (ω 2 )ω • 1 + b 1 (ω 1 )ω • 2 = b 2 (ω 2 )(ω 1 + ∂.a 1 (ω 1 )) + b 1 (ω 1 )(ω 2 + ∂.a 2 (ω 2 )) .

Using in addition (e •

3 ) and setting ω

• = def ω • 1 = ω • 2 finally yields ω • = b 2 (ω 2 )ω 1 + b 1 (ω 1 )ω 2 b 1 (ω 1 ) + b 2 (ω 2 ) + ∂ a 1 (ω 1 )b 2 (ω 2 ) + a 2 (ω 2 )b 1 (ω 1 ) b 1 (ω 1 ) + b 2 (ω 2 ) . (26)
The symbolic-only rewriting used here actually alleviates any difficulty in standardizing equation [START_REF] Dunford | Linear Operators, Part I, General Theory[END_REF]: this can now be done by just setting ∂ = 0 in its right-hand side. This yields, by identifying st(ω i ) = ω - i and st(ω

• i ) = ω + i : ω + 1 = ω + 2 = b 2 (ω - 2 )ω - 1 + b 1 (ω - 1 )ω - 2 b 1 (ω - 1 ) + b 2 (ω - 2 ) , (27)
that is, the reset value for the two velocities is the weighted arithmetic mean of ω - 1 and ω - 2 . Eq. ( 27) provides us with the reset values for the positions in the engaged mode, which is enough to restart the simulation in this mode. The actual impulsive values for the torques are useless and can be discarded. Note, however, that their sum explicitly remains equal to 0 during the instant of mode change. As a result, the momentum is preserved by the mode change, so that the obtained solution is satisfactory from a physical point of view.

As a final observation, instead of computing the exact standard part of ω • i , one could instead attempt to approximate it numerically by substituting ∂ with a small (non-infinitesimal) step size δ in System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]. This alternative approach is developed in Section 11.3. It has the advantage of not requiring the elimination of the impulsive variables.

We thus achieved the objectives stated at the end of Section 1.2.1. Figure 9 shows a simulation of the clutch model where the resets are computed following our approach developed hereabove. As one may expect, the reset value sits between the two values of ω - 1 and ω - 2 when γ : f → t (at t = 5s), and the transition is continuous at the second reset (at t = 10s).

Rejecting the submitted model

We first explain why we reject the model proposed in Section 1.2.2, and how we do so. Let us consider again System [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF].

Nonstandard semantics Following Section 2.3 for the clutch, the derivatives x , y , x , y are replaced using the explicit first-order Euler scheme with an infinitesimal time step ∂:

x ← x • -x ∂ ; x ← x •2 -2x • + x ∂ 2 . ( 28
)
In the sequel, when referring to model [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF], substitution (28) will be implicit.

Structural analysis x, y, x • , y • are the state variables, the value of which is known from previous nonstandard instants. The system has guard γ. In contrast to the clutch example, γ is no longer an input, but is rather evaluated based on system variables using equation (k 0 ). For both modes, the leading variables are x •2 , y •2 , λ; the first two of these arise from the second derivatives x , y by using substitution [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF].

Comment 5 (logico-numerical fixpoints) We could consider the direct solving of System (5), seen as a mixed logico-numerical system of equations with dependent variables x •2 , y •2 , λ, γ. However, this system mixes guarded numerical equations with the evaluation of a predicate. A possible solution would consist in performing a relaxation, by iteratively updating the numerical variables based on the previous value for the guards, and then re-evaluating the guard based on the updated values of the numerical variables, hoping for a fixpoint to occur. Such fixpoint equation, however, can have zero, one, several, or infinitely many solutions. No characterization exists that could serve as a basis for a (graph-based) structural analysis. We thus decided not to try to solve such mixed logico-numerical systems.

As a consequence, we are unable to evaluate guard γ, so that the mode the system is in cannot be determined. Should we give up? Not yet. It could be that, based on information obtained by solving the subsystem consisting of the equations that are always enabled, guard γ can be evaluated. In this model, this subsystem is G = {(e 1 ), (e 2 )}, with dependent variables x •2 , y •2 , λ. Unfortunately, subsystem G is underdetermined: it cannot be solved, and the evaluation of guard γ is still impossible. We thus reject model ( 5), because we are unable to evaluate guard γ on the basis of the subsystem that is enabled prior to evaluating this guard.

Comment 6 (forbidding logico-numerical fixpoints based on syntax) The reader may wonder why we did not simply reject the model right from the beginning, on the sole basis of the existence of the logico-numerical fixpoint equation. The rationale is that the method discussed above allows, in principle, more models to get accepted than just the ones without such equations. This design issue will be addressed more thoroughly in Comment 32 of Section 7.1.2.

Correcting the model

To break the fixpoint equation defining γ, the idea is to replace the predicate s ≤ 0 defining the guard γ by its left-limit: γ = [s -≤ 0], where s -(t) = def lim u t s(u):

                   0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) γ = [s -≤ 0] (k 0 ) if γ then 0 = L 2 -(x 2 +y 2 ) (k 1 ) and 0 = λ + s (k 2 ) if not γ then 0 = λ (k 3 ) and 0 = (L 2 -(x 2 +y 2 )) -s (k 4 )

Inria

We map this model to nonstandard semantics by stating that predicate s ≤ 0 defines the value of the guard at the next nonstandard instant-therefore, an initial condition for guard γ is required; we initialize the system in free motion mode γ(0) = f. This yields the corrected model [START_REF] Hasuo | Metamathematics for systems design -comprehensive transfer of formal methods techniques to cyber-physical systems[END_REF], where the modification is highlighted in red:

                   0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) γ • = [s ≤ 0]; γ(0) = f (k 0 ) if γ then 0 = L 2 -(x 2 +y 2 ) (k 1 ) and 0 = λ + s (k 2 ) if not γ then 0 = λ (k 3 ) and 0 = (L 2 -(x 2 +y 2 )) -s (k 4 ) ( 29
)
This model is understood in the nonstandard semantics, meaning that the derivatives are expanded using [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF]. Therefore, the leading variables in all modes are λ, s, x •2 , y •2 .

Nonstandard structural analysis

This section focuses on the core difficulty, namely the mode change γ : f → t, i.e., the event when the rope gets straight. Due to equation (k 1 ), the mode γ = t (where the rope is straight) requires index reduction. We thus augment model ( 29) with the two latent equations shown in red:

                           0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) γ • = [s ≤ 0]; γ(0) = f (k 0 ) if γ then 0 = L 2 -(x 2 +y 2 ) (k 1 ) and 0 = L 2 -(x 2 +y 2 ) • (k • 1 ) and 0 = L 2 -(x 2 +y 2 ) •2 (k •2 1 ) and 0 = λ + s (k 2 ) if not γ then 0 = λ (k 3 ) and 0 = (L 2 -(x 2 +y 2 )) -s (k 4 ) (30)
In this model, the two equations (k 1 ) and (k • 1 ) are in conflict with equations from the previous two instants, shown in blue in the following subsystem:

                     0 = x-2 • x+ •2 x ∂ 2 + •2 λ •2 x ( •2 e 1 ) 0 = y-2 • y+ •2 y ∂ 2 + •2 λ •2 y + g ( •2 e 2 ) 0 = x • -2x+ • x ∂ 2 + • λ • x ( • e 1 ) 0 = y • -2y+ • y ∂ 2 + • λ • y + g ( • e 2 ) 0 = L 2 -(x 2 +y 2 ) (k 1 ) 0 = L 2 -(x 2 +y 2 ) • (k • 1 )
Applying again Principle (1) of causality, we erase, in model [START_REF] Heemels | On the dynamic analysis of piecewise-linear networks[END_REF], equations (k 1 ) and (k • 1 ) at the instant of mode change • γ=f, γ=t. This yields the following system: at

• γ=f γ=t :

       0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) 0 = L 2 -(x 2 +y 2 ) •2 (k •2 1 ) 0 = λ + s (k 2 ) (31)
It uniquely determines all the leading variables from the state variables x, y and x • , y • . In turn, equations (k 1 ) and (k • 1 ), which were erased from this model, are not satisfied.

Also, we erase, in model ( 30), the only equation (k 1 ) at the next instant, i.e., when •2 γ=f, • γ=t, γ=t. This yields the following system:

at   •2 γ=f • γ=t γ=t   :            0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) 0 = L 2 -(x 2 +y 2 ) • (k • 1 ) 0 = L 2 -(x 2 +y 2 ) •2 (k •2 1 ) 0 = λ + s (k 2 ) (32)
This uniquely determines all the leading variables, given values for the state variables x, y and x • , y • . Note that (k • 1 ) is a consistency equation that is satisfied by the state variables x • , y • . In turn, equation (k 1 ), which was erased from this model, is not satisfied. At subsequent instants, equation erasure is no longer needed. This completes the nonstandard structural analysis of the mode change γ : f → t, i.e., when the rope gets straight.

Standardization

Code generation for restarts consists in standardizing nonstandard systems [START_REF] Höger | Dynamic structural analysis for daes[END_REF] and [START_REF] Höger | Elaborate control: variable-structure modeling from an operational perspective[END_REF]. Remember, from our development of the clutch example, that standardizing systems of equations requires more care than standardizing numbers, due to impulsive behaviors and singularity issues that result. The simplest method from a conceptual point of view is to eliminate impulsive variables from the restart system, as they are of no use for restarting the new mode.

We focus on the standardization of the mode change γ : f → t, i.e., when the rope gets straight. Our task is to standardize systems ( 31) and ( 32), by targeting discrete-time dynamics, for the two successive instants composing the restart phase. This will provide us with restart values for positions and velocities.

Due to the expansion of derivatives in equations (e 1 , e 2 , e • 1 , e • 2 ), tensions λ and λ • are both impulsive, hence so are s and s • by (k 2 , k • 2 ). We eliminate the impulsive variables by ignoring (k 2 , k •

2 ), combining (e 1 ) and (e 2 ) to eliminate λ, and (e • 1 ) and (e • 2 ) to eliminate λ • . This yields:

at • γ=f γ=t : 0 = y x + gx -x y 0 = L 2 -(x 2 +y 2 ) •2 (33) at   •2 γ=f • γ=t γ=t   :    0 = y x + gx -x y 0 = L 2 -(x 2 +y 2 ) • 0 = L 2 -(x 2 +y 2 ) •2 (34)
In System [START_REF] Iwasaki | Modeling time in hybrid systems: How fast is "instantaneous[END_REF], we expand second derivatives using [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF]. Note that substitution [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF] implies that we can as well expand the second derivatives as functions of first derivatives:

x ← x • -x ∂ . (35)
Since we will use System [START_REF] Liberzon | Switched nonlinear differential algebraic equations: Solution theory, lyapunov functions, and stability[END_REF] to evaluate restart values for the velocities, we will use expansion [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF] in it. Consequently, System (33) has dependent variables x •2 , y •2 , whereas System (34) has dependent variables x • , y • . We are now ready to standardize the two systems.

System (33) to define restart positions We expand second derivatives using (28):

0 = (y •2 -2y • + y)x -(x •2 -2x • + x)y + ∂ 2 gx 0 = L 2 -(x 2 +y 2 ) •2 (36)
Regarding [START_REF] Mattsson | Multi-Mode DAE Systems with Varying Index[END_REF], one can safely set ∂ = 0 to standardize it, since the resulting system is still structurally regular-see Comment 4 regarding the clutch example. The standardization of ( 36) is thus:

0 = (y •2 -2y • + y)x -(x •2 -2x • + x)y 0 = L 2 -(x 2 +y 2 ) •2 (37)

Inria

In contrast, had we directly set ∂ = 0 in System (31) (without eliminating impulsive variable λ), we would get a structurally singular system, an incorrect standardization.

In System (37), we can interpret x and x • as the left-limit x -of state variable x in previous mode, and x •2 as the restart value x + for the new mode. This yields

0 = (y + -y -)x --(x + -x -)y - 0 = L 2 -(x 2 +y 2 ) + (38)
which determines the restart values for positions. Note that the constraint that the rope is straight is satisfied. Furthermore, if we remember that the rope is straight at mode change (this is the mode change condition), then 0 = L 2 -(x 2 +y 2 ) -also holds, which ensures that x + = x -, y + = y -is the unique solution of ( 38): positions are continuous.

System (34) to define restart velocities We expand second derivatives using [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF]:

   0 = (y • -y )x -(x • -x )y + ∂.gx 0 = L 2 -(x 2 +y 2 ) • 0 = L 2 -(x 2 +y 2 ) •2 (39)
By expanding x •2 = x • + ∂x • , the right-hand side of the last equation rewrites

L 2 -(x 2 +y 2 ) •2 = L 2 -(x 2 +y 2 ) • + 2∂(x • x • + y • y • ) + ∂ 2 (x • ) 2 + (y • ) 2
= 0 (by second equation of ( 39))

+ 2∂(x • x • + y • y • ) + O(∂ 2 ) (40)
Using this expansion of

L 2 -(x 2 +y 2 ) •2 , setting ∂ = 0 in (39) yields 0 = (y • -y )x -(x • -x )y 0 = x • x • + y • y • (41)
where the dependent variables are now x • , y • , whereas other variables are state variables whose values are determined by previous time steps. Note that System (41) is structurally regular, hence it is the correct standardization of System [START_REF] Nakamura | A nonstandard functional programming language[END_REF]. We are now ready to get restart values for velocities.

In System (41), we perform the following substitutions, recalling that superscripts -and + denote left-and right-limits:

x = x -; x • = x + ; x = x - (42)
and the restart velocity is evaluated by

x + = x • . ( 43
)
After similar substitutions for y, we get

0 = (y + -y -)x --(x + -x -)y - 0 = x + x + + y + y + (44)
System (44) determines x + and y + , which are the velocities for restart. The second equation guarantees that the velocity will be tangent to the constraint. With ( 38) and ( 44), we determine the restart conditions for positions and velocities. Invariants from the physics are satisfied. Our reasoning so far produces a behavior in which the two modes (free motion and straight rope) gently alternate; the system always stays in one mode for some positive period of time before switching to the other mode. This indeed amounts to assuming that the impact is totally inelastic at mode change, an assumption that was not explicit at all in [START_REF] Hasuo | Metamathematics for systems design -comprehensive transfer of formal methods techniques to cyber-physical systems[END_REF]. So, what happened?

Handling transient modes

In fact, the straight rope mode was implicitly assumed to last for at least three nonstandard successive instants, since we allowed ourselves to shift (k 1 ) twice. Now, let us instead assume elastic impact. Hence, the cascade of mode changes γ : f → t → f is possible in three successive nonstandard instants, reflecting that the straight rope mode is transient (it is left immediately after being reached). We wish to adapt our reasoning to handle this transient mode.

Consider again model [START_REF] Hasuo | Metamathematics for systems design -comprehensive transfer of formal methods techniques to cyber-physical systems[END_REF]. We regard the instant of the cascade when γ = t occurs as the current instant. We cannot add latent equations by simply shifting (k 1 ), since these shifted versions are not active in the mode γ = f. Instead, our line of reasoning goes as follows. Set

S(t) = {(e 1 ), (e 2 ), (k 1 ), (k 2 )} , S(f) = {(e 1 ), (e 2 ), (k 3 ), (k 4 )} .
Systems S • (t) and S • (f) are obtained by shifting the equations constituting S(t) and S(f); systems S •2 (t) and S •2 (f) are defined similarly. Now, the idea is to consider the differentiation array originally proposed by Campbell and Gear [START_REF] Campbell | The index of general nonlinear DAEs[END_REF], except that we take into account the trajectory t, f, f, . . . for guard γ. Using shifting instead of differentiation yields the following array of shifted systems, where size n is a parameter:

A n (S) = def        S(t) S • (f) S •2 (f) . . . S •n (f)        (45)
We search for an n such that A n (S) provides us with the due latent equations. Thus, the question is: can we find latent equations for S(t) by shifting appropriate equations from S(f)? Unfortunately, the answer is no: although shifting twice (k 4 ) in System (29) produces one more equation involving the leading variables x •2 , y •2 , this equation also involves the new variable s •2 , which keeps the augmented system underdetermined; shifting other equations fails as well. Therefore, the structural analysis rejects this model as being underdetermined at transient mode γ = t. The user is then asked to provide one more equation. For example, he/she could specify an impact law for the velocity y by providing the equation

(y ) + = -(1 -α)(y ) -,
where 0 ≤ α < 1 is a fixed damping coefficient. This is reinterpreted in the nonstandard domain as

y • = -(1 -α)y , (46)
yielding the following refined system for use at mode γ=t within the cascade γ:f→t→f:

           0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) 0 = y • + (1 -α)y (τ 1 ) 0 = L 2 -(x 2 +y 2 ) (k 1 ) 0 = λ + s (k 2 ) (47)
We proceed again with the structural analysis. Variables x, y are the states, so that their values are set by the previous instants. Current equation (k 1 ) creates a conflict with the past. Hence, we discard it from System [START_REF] Potop-Butucaru | Compiling Esterel[END_REF], which leaves us with the following system:

       0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) 0 = y • + (1 -α)y (τ 1 ) 0 = λ + s (k 2 ) ( 48
)
Inria Model ( 48) is structurally nonsingular, recalling that y and y • can be interchanged for the structural analysis. This refined model is therefore accepted. Standardization proceeds accordingly, with no new difficulty.

Consequences for the modeling language

Through the Cup-and-Ball example, we demonstrated the need for the following user-given information: is the current mode long or transient? We argue that this information cannot be derived from inspecting the system model.

In the clutch example, guard γ was a system input, and we implicitly assumed that each of the two modes was long. Had we composed the clutch model with a 'feedback' determining γ as a predicate over system variables, this implicit assumption would have become questionable. This situation indeed occurs with the Cup-and-Ball example, where the elastic/inelastic impact assumption is a side information that cannot be deduced from the equation syntax. This leads to the following important observation:

Comment 7 (Long / Transient) Long / Transient is an information regarding modes, that cannot be found by inspecting the equations of the system. It must be inferred from understanding the system physics. Therefore, it has to be provided as an extra information by the model designer.

The natural way of addressing Comment 7 is to provide a different syntax for specifying long modes on the one hand, and events corresponding to transient modes on the other hand. In contrast, mode changes separating two successive long modes need not be specified. A candidate syntax is:

long mode if predicate then S transient mode when predicate then S (49)
The 'when' statement of Modelica, which we reuse here, points to the event when 'predicate' switches from false to true. For this 'when' statement, we could furthermore restrict 'predicate' to be a zero-crossing condition, by which a real-valued expression crosses zero from below [START_REF] Bourke | Zélus: A synchronous language with ODEs[END_REF]. We devote the 'if' statement to long-lasting modes specified by 'predicate'.

To illustrate this, our syntax offering if-equations and when-equations allows specifying the Cup-and-Ball example with elastic impact as shown in Figure 10.

               0 = x + λx (e 1 ) 0 = y + λy + g (e 2 ) γ = [s -≤ 0]; γ(0) = f (k 0 ) when γ then y = -αy - (τ 1 ) if not γ then 0 = λ (k 3 ) and 0 = (L 2 -(x 2 +y 2 )) -s (k 4 )
Figure 10: Model of the Cup-and-Ball example with elastic impact specified by equation (τ 1 ), highlighted in red. This equation replaces the two equations (k 1 ), (k 2 ) specifying the straight rope mode in case of inelastic impact.

The Westinghouse air brake example

In this section we develop the example of Section 1.2.3 and we address the objectives stated at the end of that section.

We reject the model proposed in Section 1.2.3 for the very same reasons as for the Cup-and-Ball example, namely the existence of a logico-numerical fixpoint equation defining the guard. The corrected system is shown in Figure 11. The two systems differ in equation (v 0 ), where the predicate

0 = f b -f v -f cl + f l (m 1 ) 0 = f v -f ch -f l -f t (m 2 ) 0 = f b -F 1 .ϕ(p b -p r ) (f 1 ) 0 = F 1 .ϕ(p bn -p r ) (f 2 ) 0 = f l -F 2 .ϕ(p t -p r ) (f 3 ) 0 = f t -ρ.V P 0 .p t (r) 0 = ρ.S.(x .p r + (x -L).p r ) + P 0 .f cl (p 1 ) 0 = ρ.S.(x .S.p t + x.p t ) -P 0 .f ch (p 2 ) 0 = S.(p t -p r ) -b (l 1 ) 0 = K.x -b (l 2 ) γ • = [s ≤ 0]; γ(0) = f (v 0 ) if γ then 0 = p r -p t (v 1 ) and 0 = f v + s (v 2 ) if not γ then 0 = f v (v 3 ) and 0 = p r -p t -s (v 4 )
Figure 11: The system of of Figure 8 corrected. The correction is highlighted in red.

s < 0 sets the value of guard γ for the next nonstandard instant. This models the fact that γ is evaluated based on • s, which is the nonstandard semantics of the left-limit of s. Note that an initial condition for γ is needed in this corrected version. This model is interpreted in the nonstandard semantics, meaning that substitution (28) applies.

Nonstandard structural analysis

The value of guard γ is known from the previous instant, together with the value for the two state variables x, p r , p t . For both modes γ = t/f, the leading variables for evaluation are the 12

variables x • , p • r , p • t , f b , f v , f cl , f l , f ch , f t , p bn , b
, s, hence the active system is square in both modes. We successively explore the two long modes γ = t/f for the guard. Then we investigate the mode changes.

Long mode γ = f. The active system is G f = G ∪ {(v 3 ), (v 4 )}, see Figure 12.

The subsystem {(l 1 ), (l 2 )} involves two equations but only one leading variable: b. It is structurally singular. Following the technique of index reduction used in the clutch example, we augment G f with the latent equations (l • 1 ), (l • 2 ). The result is shown in Figure 13. Note that b is now a state variable. Equations (l 1 ), (l 2 ) are consistency conditions. They are satisfied as a result of executing the previous instant. We are thus left with the system

G Σ f = G ∪ {(l • 1 ), (l • 2 )}.
Can we check its regularity by using structural (graph based) methods?

A complete matching is shown in red in Figure 13, i.e., a one-to-one assignment of leading variables to equations. The intent is that each variable is evaluated using its assigned equation, by pivoting. For example, f cl is evaluated using equation (m 1 ). The heart of structural analysis for algebraic systems of equations is that, if a complete matching exists, then the system is regular in a generic sense, i.e., outside exceptional values for the non-zero coefficients of the equations. See Section 6 for a formal development.

The system G Σ f is thus structurally nonsingular. Hence, the system of Figure 12 was of index 1 and adding the latent equations reduced the index to 0.

Long mode γ=t. The active system is G t = G ∪ {(v 1 ), (v 2 )}, shown in Figure 14. The subsystem {(l 1 ), (l 2 ), (v 1 )} is structurally singular since it has 3 equations and only one dependent variable b. We thus shift it forward. The result is shown in Figure 15.

Equations (l 1 ), (l 2 ), (v 1 ) are consistency conditions. They are satisfied as a result of executing the previous instant. We are thus left with the index reduced system

G Σ t = G ∪ {(l • 1 ), (l • 2 ), (v • 1 )
}, for which a complete matching is shown in red in Figure 15. This system is structurally nonsingular.

Inria 0 = f b -f v -f cl + f l (m 1 ) 0 = f v -f ch -f l -f t (m 2 ) 0 = f b -F 1 .ϕ(p b -p r ) (f 1 ) 0 = F 1 .ϕ(p bn -p r ) (f 2 ) 0 = f l -F 2 .ϕ(p t -p r ) (f 3 ) 0 = f t -ρ.V P 0 .p t (r) 0 = ρ.S.(x .p r + (x -L).p r ) + P 0 .f cl (p 1 ) 0 = ρ.S.(x .S.p t + x.p t ) -P 0 .f ch (p 2 ) 0 = S.(p t -p r ) -b (l 1 ) 0 = K.x -b (l 2 ) 0 = f v (v 3 ) 0 = p r -p t -s (v 4 )
Figure 12: The active system when γ = f. We show in red a structurally singular subsystem.

0 = f b -f v -f cl + f l (m 1 ) 0 = f v -f ch -f l -f t (m 2 ) 0 = f b -F 1 .ϕ(p b -p r ) (f 1 ) 0 = F 1 .ϕ(p bn -p r ) (f 2 ) 0 = f l -F 2 .ϕ(p t -p r ) (f 3 ) 0 = f t -ρ.V P 0 .p t (r) 0 = ρ.S.(x .p r + (x -L).p r ) + P 0 .f cl (p 1 ) 0 = ρ.S.(x .S.p t + x.p t ) -P 0 .f ch (p 2 ) 0 = S.(p • t -p • r ) -b • (l • 1 ) 0 = K.x • -b • (l • 2 ) 0 = f v (v 3 ) 0 = p r -p t -s (v 4 ) consistency equations: 0 = S.(p t -p r ) -b (l 1 ) 0 = K.x -b (l 2 )
Figure 13: The system of Figure 12 in which latent equations are found by shifting (l 1 ), (l 2 ). Recall that x and x • are related by [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF]. The assignment of a variable to each equation, highlighted in red, defines a complete matching.

The system of Figure 14 was of index 1. The set of latent equations differs from that of the other mode, however.

So far we have completed the study of the long modes. It remains to investigate the mode changes.

Mode change γ : t → f. We reproduce the reasoning following System (20) for the clutch example. Joining together the systems at previous and current instant may result in a conflict between the consistency equations of the current instant and selected equations from the previous instant. This does not occur here, however: consistency equations (l 1 ), (l 2 ) in Figure 13 coincide with the latent equations (l • 1 ), (l • 2 ) associated with the previous instant in Figure 15 (the set of latent equations of mode γ = t contains the set of latent equations of mode γ = f).

Mode change γ : f → t. Unlike for the previous case, joining together the systems at previous and current instant actually results in a conflict between the consistency equations of the current instant and the dynamics at the previous instant. More precisely, the consistency equation (v 1 ) : 0=p r -p t , which holds in current mode γ = t, is not a consequence of executing the previous instant (in mode

0 = f b -f v -f cl + f l (m 1 ) 0 = f v -f ch -f l -f t (m 2 ) 0 = f b -F 1 .ϕ(p b -p r ) (f 1 ) 0 = F 1 .ϕ(p bn -p r ) (f 2 ) 0 = f l -F 2 .ϕ(p t -p r ) (f 3 ) 0 = f t -ρ.V P 0 .p t (r) 0 = ρ.S.(x .p r + (x -L).p r ) + P 0 .f cl (p 1 ) 0 = ρ.S.(x .S.p t + x.p t ) -P 0 .f ch (p 2 ) 0 = S.(p t -p r ) -b (l 1 ) 0 = K.x -b (l 2 ) 0 = p r -p t (v 1 ) 0 = f v + s (v 2 )
Figure 14: The active system when γ = t. We show in red a structurally singular subsystem.

Compare with Figure 12.

0 = f b -f v -f cl + f l (m 1 ) 0 = f v -f ch -f l -f t (m 2 ) 0 = f b -F 1 .ϕ(p b -p r ) (f 1 ) 0 = F 1 .ϕ(p bn -p r ) (f 2 ) 0 = f l -F 2 .ϕ(p t -p r ) (f 3 ) 0 = f t -ρ.V P 0 .p t (r) 0 = ρ.S.(x .p r + (x -L).p r ) + P 0 .f cl (p 1 ) 0 = ρ.S.(x .S.p t + x.p t ) -P 0 .f ch (p 2 ) 0 = S.(p • t -p • r ) -b • (l • 1 ) 0 = K.x • -b • (l • 2 ) 0 = p • r -p • t (v • 1 ) 0 = f v + s (v 2 ) consistency equations: 0 = S.(p t -p r ) -b (l 1 ) 0 = K.x -b (l 2 ) 0 = p r -p t (v 1 )
Figure 15: Adding latent equations to the system of Figure 14.

γ = f). As for the clutch example (Section 2.4 and particularly ( 1)),

we remove, at the instant of mode change γ :

f → t, the consistency equation (v 1 ) : 0=p r -p t . ( 50
)
This completes the structural analysis.

Standardization

To derive the actual executable code we must return to the standard domain by performing standardization. The task is very much like for the clutch example and we only provide a brief summary:

Within long modes γ = f or t: Here, time and dynamics are continuous: We standardize x • -x ∂ as the derivative x . Also, we standardize equation (v 0 ) as γ=[s -≤0], where s -denotes the left limit of s, i.e., the value of s before the mode change in case a mode change occurs. Doing this removes all the occurrences of ∂ in both modes. The resulting system is standard.

Inria

For mode changes: Here we target discrete time: x • is interpreted as the next value for x and raises no difficulty. In contrast, the occurrence of ∂ in the Euler schemes x • -x ∂ , where it acts in space, is now the difficulty.

Mode change γ : t → f. The restart conditions are solutions of the system of Figure 13 in which the consistency equations are ignored. This system coincides with the nonstandard dynamics of the long mode γ = f. However, we do not seek for a strandardization in continuous time, but rather in discrete time, which is different: derivatives are expanded using [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF] and we solve the system for the discrete time leading variables p • t , p • r , x • , etc. This nonstandard system is shown in Figure 16, where the infinitesimal ∂ is highlighted in blue. We need to standardize this system. We leave to the reader as an exercise that setting ∂ = 0 in this system makes it structurally nonsingular (no complete matching can be found). By Comment 4 of Section 2.5.2, this implies that setting ∂ = 0 in this system does not lead to a correct standardization.

0 = f b -fv -f cl + f l (m1) 0 = fv -f ch -f l -ft (m2) 0 = f b -F1.ϕ(p b -pr) (f1) 0 = F1.ϕ(p bn -pr) (f2) 0 = f l -F2.ϕ(pt -pr) (f3) 0 = ∂.ft -ρ.V P 0 .(p • t -pt) (r) 0 = ρ.S.((x • -x).pr + (x -L).(p • r -pr)) + ∂.P0.f cl (p1) 0 = ρ.S.((x • -x).S.pt + x.(p • t -pt)) -∂.P0.f ch (p2) 0 = S.(p • t -p • r ) -b • (l • 1 ) 0 = K.x • -b • (l • 2 ) 0 = fv (v3) 0 = pr -pt -s (v4)
The problem originates from equations (r), (p 1 ), (p 2 ). Focus on (r): the state variable x is finite before and after the mode change and x • -x = 0 possibly holds; then, ∂ infinitesimal requires that f t is possibly infinite, i.e., impulsive. The same holds for f cl and f ch .

Following the approach of Section 2.5.2, we must eliminate, by algebraic manipulations, the possibly impulsive variables f t , f cl , f ch . This is easy since the latter enter linearly in equations (r), (p 1 ), (p 2 ).

Setting ∂ = 0 in the resulting system leaves us with a structurally singular system that sets the restart values for the system state variables p r , p t , b, x. This yields the desired standardization, which is the code for correct restart at this mode change. Details are omitted.

Mode change γ : f → t. We proceed similarly. The restart equations at mode change γ : f → t are given in Figure 17, where the infinitesimal ∂ is highlighted in blue. We conclude as for the other Comment 8 (fake impulsive variables) The reader acquainted with the physics of this system may be quite surprised by our reasoning above. We found that, for the two mode changes, algebraic variables f t , f cl , and f ch are possibly impulsive. However, the physics of the Westinghouse air brake does not tell so. There seems to be something wrong here.

0 = f b -fv -f cl + f l (m1) 0 = fv -f ch -f l -ft (m2) 0 = f b -F1.ϕ(p b -pr) (f1) 0 = F1.ϕ(p bn -pr) (f2) 0 = f l -F2.ϕ(pt -pr) (f3) 0 = ∂.ft -ρ.V P 0 .(p • t -pt) (r) 0 = ρ.S.((x • -x).pr + (x -L).(p • r -pr)) + ∂.P0.f cl (p1) 0 = ρ.S.((x • -x).S.pt + x.(p • t -pt)) -∂.P0.f ch (p2) 0 = S.(p • t -p • r ) -b • (l • 1 ) 0 = K.x • -b • (l • 2 ) 0 = p • r -p • t (v • 1 ) 0 = fv + s (v2)
In the next section we investigate how the right solution can be found at compile time.

Assertions

Focus on the system for restart at mode change γ : f → t, see Figure 17 and assume for a while that the trajectories of system variables remain continuous, which applies in particular to s. Now, the predicate defining the guard γ tells us that the value of s at both mode changes is zero. As a consequence, at the instant of both mode changes γ : f → t and t → f, the two subsystems {(v 1 ), (v 2 )} and {(v 3 ), (v 4 )} of the system of Figure 11 become identical due to the zero value for s.

As a consequence, whatever next mode is, the consistency equations for it ({(v 1 ), (v 2 )} or {(v 3 ), (v 4 )}) will be satisfied. Therefore, no action is required at mode changes other than commuting from code of Figure 13 to code of Figure 15 or vice-versa. Since our compiler is only equipped with our graph based structural analysis, it will not be able to perform the above reasoning. We can, however, overcome this difficulty by giving to the user the possibility to manually refine her model by adding guarded assertions of the form

when γ assert β(X) (51)
where γ is some guard and β(X) is a predicate in the numerical variables X of the system. For example, in the Westinghouse example, the user could state the following assertion, which follows from the fact that s = 0 holds at mode changes: when γ -and not γ assert s = 0 (52) Assertion ( 52) can be used by the compiler in deducing that no specific restart action is required other than commuting between the code of the two modes, which, in turn, implies that no state variable is impulsive.

Comment 9 (are assertions needed for correctness?) Not quite. In fact, the approach followed in Section 4.2 for deriving the restart conditions at mode changes is still correct. It is just incomplete in that it does not allow to compute the value for the algebraic variables that were erroneously considered impulsive. This is indeed harmless. Still, our approach of Section 4.2 was an overshoot and adding a proper assertion is preferred. Also, assertions can be very useful in narrowing the domain of reachable modes, which may improve efficiency of the compilation.

Assertions are provided in our IsamDAE tool under development [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF][START_REF] Caillaud | Implicit Structural Analysis of Multimode DAE Systems[END_REF], see also Section 12. We will not discuss them any further in this paper, however.

Issues and systematic approach

From our discussion of the examples, we see that the key issue is the handling of mode changes, and particularly the reconciliation of the dynamics just before the change, and the consistency conditions set by the new mode. How to properly address this in the original model of continuous real-time was an open question; mapping the system dynamics to the nonstandard domain clarified this task.

In the rest of this paper, we extend our method, from the toy examples to a systematic approach for mDAE system compilation. This approach, illustrated in Figure 18, decomposes into several steps.

Inria

• First, the mDAE system is transformed into a system of multi-mode difference Algebraic Equations (mdAE) using the nonstandard Euler expansion of derivatives x ← x • -x ∂ . The resulting mdAE system is the nonstandard semantics of the original model.

• Second, the structural analysis of this nonstandard semantics is performed, resulting in a new mdAE system where latent equations are made explicit and conflicts at mode changes are solved.

• Finally, the standardization step recovers both the smooth dynamics in each mode, and the possibly discontinuous/impulsive restart conditions at mode changes. The resulting code can then be executed with the help of a solver. The intuitions developed on the examples are formalized in the sequel of the paper. We first recall the background on structural analysis and we extend it to our context. We then focus on standardization.

Background on Structural Analysis

We recall in this section the basics of structural analysis for algebraic systems of equations and provide in passing a few new lemmas for subsequent use. We then present the Σ-method, which is the structural analysis method we adopt for DAE systems.

Notations 10 In this paper, N = {0, 1, 2, . . . } is the set of nonnegative natural integers and R + = [0, ∞) is the set of nonnegative real numbers.

We will sometimes 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 11 Lowercase letters (x, y, etc.) denote scalar real variables; capitals (X, Y , etc.) 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, etc.) and vectors of functions (F, G, etc.). A value for a variable x is generically denoted by ν x , and similarly for X and ν X .

Structural analysis of algebraic equations

In this section we focus on systems of algebraic equations, that is, equations involving no time and no dynamics.

Structural nonsingularity of algebraic equations

Consider a system of smooth algebraic equations:

f i (y 1 , . . . , y k , x 1 , . . . , x n ) = 0, i = 1, . . . , m (53)
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 ). The system has m equations, k input variables 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 C 1 at least. If the system 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. ( 53) 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 considering J X F by focusing on structural nonsingularity, a weaker but much cheaper criterion previously introduced for handling sparse matrices in high performance computing.

Structural analysis

We begin with some background on graphs, for which a basic reference is [START_REF] Berge | The theory of graphs and its applications[END_REF]. Given a graph G = (V, E), where V and E 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 say simply "alternating path" when M is understood from the context. A vertex is said to be matched if there is an edge in M that is incident on the vertex, and unmatched otherwise.

Let F (Y, X)=0 be a system of algebraic equations of the form [START_REF] Trenn | Regularity of distributional differential algebraic equations[END_REF]; its bipartite graph G F is the graph having F ∪X as 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 equations and dependent variables are in equal numbers.

Definition 12 (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 variables in equal numbers. The link to numerical regularity is formalized in the following lemma (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]):

Lemma 13 Assume that F is 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 generi-
cally 11 nonsingular when its nonzero coefficients vary over some neighborhood.

Let us detail Property 2. We can represent graph G F by the n×n-incidence matrix M F = (m ij ), having a 1 at (i, j) if (f i , x j ) is an edge of G F , and a 0 otherwise. We order the entries {a ij | m ij = 1} of Jacobian J X F (ν Y , ν X ), for example by lexicographic order over the pair (i, j). In this way, every n×n-matrix whose pattern of nonzero coefficients is M F identifies with a unique vector of R K , where K is the number of 1's in M F . 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.

(54)

In the sequel, the so defined sets U and V will be denoted by

U F (ν Y , ν X ) and V F (ν Y , ν X ) (55)
or simply U F and V F when no confusion can result. Practical use of structural nonsingularity. 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 12 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 . 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 G M
F the resulting directed graph. The following result holds ([25], Chapter 6.10):

Lemma 14 The strongly connected components of G M F are independent of M. 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 1 f 2 f 3 f 4 f 5 f 6 x 1 x 2 x 3 x 6 x 4 x 5 f 1 f 2 f 3 f 4 f 5 f 6 f 1 f 2 f 3 f 4 f 5 f 6 f 1 f 2 f 3 f 4 f 5 f 6 x 1 x 2 x 3
x 6

x 4

x 5 Nonsquare systems In our development, we will also encounter nonsquare systems of algebraic equations. For general systems (with a number of variables not equal to the number of equations, n =m in [START_REF] Trenn | Regularity of distributional differential algebraic equations[END_REF]), call block a pair β = (f , x) collecting a subset f of the set of f i 's and a subset x of variables x j involved in it.

Definition 15 (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 15 is meaningful:

Lemma 16
1. The triple (β O , β E , β U ) defined by the Dulmage-Mendelsohn decomposition does not depend on the particular maximal matching used for its definition. We thus write

(β U , β E , β O ) = DM(F ).
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 14 to β E in order to get its BTF.

The following obvious lemma will be useful in the sequel. Comment 18 Suppose that, instead of removing from F all equations belonging to β O , we only remove all unmatched equations with reference to the matching of maximal cardinality used for generating DM(F ), and call F the remaining subsystem of F . Then, we still get that no overdetermined block exists in the Dulmage-Mendelsohn decomposition of F and we have F ⊇ (F \ β O ), the inclusion being strict in general. However, this policy depends on the particular matching. In contrast, our policy is canonical since the DM decomposition is independent from the matching. Since we will rely on Lemma 17 for our compilation scheme in Section 7, adopting the finer policy would introduce a risk of nondeterminism due to the arbitrary choice of the particular complete matching selected. This is the reason for sticking with the policy stated in Lemma 17.

Local versus non-local use of structural analysis

Structural analysis relies on the Implicit Function Theorem for its justification, as we have seen. 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 given values ν 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. 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, see our examples of Sections 2 and 3. 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 , and thus we omit Y .

Lemma 19 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:

Inria 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 19 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 19, 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 19. By Property 2 of Lemma 19, (S × R) ∩ V has dimension 1, implying that S ∩ proj R n (V ) has dimension < 1. This analysis is summarized in the following lemma: Lemma 20 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 20, 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. 12 Of course, subclasses of systems for which existence and uniqueness are guaranteed are of interest-the illustration examples we develop in this paper have this property.

Equations with existential quantifiers

To our knowledge, the results of this section are new. For the handling of mode change events, we will also need to develop a structural analysis for the following class of (possibly nonsquare) algebraic systems of equations with existential quantifier:

∃W : F (X, W, Y )=0 , (56)
where (X, W ) collects the dependent variables of system F (X, W, Y )=0. Our aim is to find structural conditions ensuring that (56) defines a partial function Y → X, meaning that the value of tuple X is uniquely defined by the satisfaction of (56), given a value for tuple Y . At this point, a fundamental difficulty arises. Whereas (56) 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 (56) only in subclasses of systems. 13 In general, no extension of the Implicit Function Theorem exists for systems of the form (56); thus, one cannot apply as such the arguments developed in Section 6.1.3. 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 ). We thus 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 , (57)
expressing that X is independent of W given a consistent tuple of values for Y . To find structural criteria guaranteeing (57), 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 21 System (56) is structurally nonsingular if the following two conditions hold, almost everywhere when the nonzero coefficients of the Jacobian matrix J X,W,Y F vary over some neighborhood: 14 consistent values for Y exist and condition (57) holds.

The structural nonsingularity of (56) 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 14 and comments thereafter. Let B be the set of all indecomposable blocks of β E and be the partial order on B following Lemma 14. The following holds:

Lemma 22 System (56) 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 β ≺β 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 .

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 "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 (57) 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 (57) 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 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 (57). This finishes the proof.

We complement Lemma 22 with the algorithm ExistQuantifEqn (Algorithm 5), which requires a system F of the form (56). If condition 1 of Lemma 22 fails to be satisfied, then b O ← f is returned, indicating overdetermination. If conditions 2 or 3 of Lemma 22 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.

The Σ-method for DAE systems

The DAE systems we consider are "square", i.e., have the form f j (the x i 's and their derivatives) = 0 (58)

where x 1 , . . . , x n denote the dependent variables and f 1 = 0, . . . , f n = 0 denote the equations. Call leading variables of System (58) the d i -th derivatives 15 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 Algorithm 5 ExistQuantifEqn Require:

F ; return (b O , b U , F Σ , F Σ ) 1: (β U , β E , β O ) = DM(F ) 2: if condition 1 of Lemma 22 holds then 3: b O ← t; 4:
if 2 and 3 of Lemma 22 hold then partition b O ← f structural analysis of DAE systems of the form (58) is the following. Regard (58) 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 (58) is like an ODE. If this is not the case, finding additional latent equations by differentiating suitably selected equations from (58) will bring the system to an ODE-like form, while not changing its set of solutions. Performing this is known as index reduction.

β E = F Σ ∪ F Σ 7:
In our simple examples, finding latent equations was easy. In general, this is difficult and 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 (58), 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 (58) the bipartite graph G having F ∪ X as set of vertices and having an edge (f, x) 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)∈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 23 Find a complete matching M for G and 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 with equality if (f, x)∈M c f ≥ 0 (59)
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 12. Hence, F Σ =0 is like an ODE. The integer k = max f ∈F c f is called the index of the system. 16Definition 24 For F a DAE system, the solution to Problem 23 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 23 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 23 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 23. As a preliminary step, we claim that a bruteforce encoding of Problem 23 is the following: Find integers ξ fx , d x , c f such that

f :(f,x)∈G ξ fx = 1 x:(f,x)∈G ξ fx = 1 ξ fx ≥ 0    complete matching d x -c f -d fx ≥ 0 c f ≥ 0 encodes "≥" in (59) (f,x)∈G ξ fx (d x -c f -d fx ) = 0 encodes "=" in (59) (60)
We now justify our claim. Focus on the first block. Since the ξ's are integers, they can only take values in {0, 1} and one can define a bipartite graph by stating that edge (f, x) exists iff ξ fx = 1.

The first two equations formalize that this graph is a matching, which is complete since all equations are involved. The second block is a direct encoding of (59) 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 (60) 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 (60) is a solution of the following dual linear programs (LP), where the ξ fx are real:

primal : maximize (f,x)∈G d fx ξ fx subject to f :(f,x)∈G ξ fx = 1 and x:(f,x)∈G ξ fx ≥ 1 and ξ fx ≥ 0 (61) dual : minimize x d x -f c f subject to d x -c f ≥ d fx and c f ≥ 0 (62)
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 the two possess identical solutions in this case. Note that LP (61) 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 (61) and (62), ξ fx >0 if and only if

d x -c f =d fx , (63)
which is exactly the last constraint of (60). Using this translation into linear programs (61) and (62), 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 23, then a unique elementwise smallest solution exists. Based on the above analysis, the following algorithm was proposed in [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF] for solving (61,62) and was proved to provide the smallest solution for (62). Let G be a bipartite graph with weights {d fx |(f, x)∈G}:

1. Solve LP (61), 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) ∈ G} ; (b) ∀f : c f ← d x -d fx where (f, x) ∈ M .
The reason for using the special iterative algorithm for solving the dual LP (62) is that a standard LP-solving algorithm will return an arbitrary solution, not necessarily the smallest one. The following lemma, that will be crucially put to use in Section 11.1, is an obvious consequence of the linearity of problems ( 61) and (62):

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

(d 1 fx ) (f,x)∈G and (d 2 fx ) (f,x)∈G be related by d 2 fx = M × d 1 fx for every (f, x) ∈ 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 . A necessary and sufficient condition for Problem 23 to have a solution is that the set of complete matchings for G is non-empty. We thus consider the function

b ← ExistsMatching(G) , (64)
which returns a Boolean b such that b = t iff there exists a complete matching for the bipartite graph G (regardless of its associated weights). In this case, step 1 of the above algorithm is guaranteed successful (since the set of complete matchings for G is nonempty and finite). In the sequel, we encode the above algorithm as the function

(c, d) ← FindOffsets(G, D) , where    c : F → N d : X → N D : F ×X → N (65)
Having b = t returned by ExistsMatching() is a prerequisite to calling FindOffsets(), which is then guaranteed successful.

Algorithm 6 SigmaMethod

Require: F ; return (b, F Σ , F Σ ) 1: b ← ExistsMatching(G). 2: if b then 3: (c, d) ← FindOffsets(G, D) 4: (F Σ , F Σ ) ← LatentEqns(c, F )
We encapsulate these two functions as Algorithm 6 (SigmaMethod). It requires a DAE system F =0. G is the weighted bipartite graph of F with weights D. SigmaMethod returns a Boolean b and, if b = t, a pair (F Σ , F Σ ) defining an index-0 system F Σ =0 and a (possibly underdetermined) system F Σ =0 of consistency constraints following Definition 24. In Line 4, F Σ collects the f c f , for f ∈ F , and F Σ collects the f k , for f ∈ F and 0 ≤ k < c f . This completes our background material on structural analysis.

Structural analysis of multimode DAE systems

In this section, we extend structural analysis, from (single-mode) DAE systems, to multimode DAE systems. As a prerequisite, we precisely define the class of multimode DAE systems considered. Since the structural analysis will actually be applied to the nonstandard semantics, we also need to define multimode difference Algebraic Systems (dAE ), operating in discrete time.

Inria

Defining multimode DAE/dAE systems

In this section, the class of systems of multimode Differential/difference Algebraic Equations considered in this paper is formally defined. By analogy with DAE's, we use the acronym dAE to mean difference algebraic equations. Notation x • denotes the shift of x, as per Definition 2.

Notations 26 Let X be an underlying set of variables. For x∈X and m∈N, the m-differentiation and m-shift of x are denoted by x m and x •m , respectively. Let X m and X •m denote the set of all x m and x •m , for x ranging over the set X of variables. Let

X ( ) = def m∈N X m and X (•) = def m∈N X •m . (66)
For X⊂X , sets X m , X •m , X ( ) , X (•) are defined in a similar way.

The following assumption will be in force in the remainder of this paper:

Assumption 1 All the functions and predicates over X ( ) are assumed to involve only finitely many of the x m for x∈X and m∈N. The same holds regarding X (•) .

Syntax and meaning

From the discussion in Section 3.6 regarding the Cup-and-Ball example, we know we need to type the modes as either long or transient, since the two situations require different index reduction techniques. This justifies the consideration of two different kinds of guarded equations. | when γ then f =0 (69) mDAE S is a finite set of s; s is either a guard evaluation g (67), or a guarded equation e, of two different kinds (68) or (69); f is a smooth scalar function over X ( ) , and b is a Boolean expression of predicates over X ( ) . The same holds for the syntax of mdAE, with X (•) replacing X ( ) .

Notations 28 Let e be a guarded equation. Its body f =0 and guard γ are denoted by f (e) and γ(e). For S an mDAE or mdAE system, we denote by X(S), E(S), and Γ(S) its sets of numerical variables, guarded equations, and guards, respectively. The mention of S is omitted when it is clear from the context. Finally, for e : if/when γ then f =0 a guarded equation following Definition 27, we define

e • /e : if/when γ then f • /f =0 ( 70
)
where f • follows Definition 2, and, for E a set of equations, Notations 26 are adapted so that E m , E •m , E ( ) , E (•) can be considered. We insist that, in (70), the guard γ is left unchanged.

Meaning The meanings of if-equation (68) and when-equation (69) differ:

• In the if-equation (68), equation f =0 is enabled if and only if the guard γ holds. Otherwise, this equation is disabled.

• In the when-equation (69), equation f =0 is enabled exactly at the events when guard γ switches from f to t. Otherwise, this equation is disabled.

For mdAE systems, the formalization of the second statement is straightforward: equation f =0 is enabled when γ ↑ holds, where γ ↑ = def (not • γ) and γ .

For mDAE systems: equation f =0 is enabled when γ ↑ holds, where γ ↑ = def (not γ -) and γ (72

)
where γ -denotes the left-limsup of γ at the current instant, defined by:

γ -(t) = def lim sup s t γ(s) = if ∃s n t : γ(s n )=t then t else f
The left-limsup is always defined and coincides with the left limit when the latter exists.

Comment 29 Why bother with the two types of equations, since when-equations are mapped to if-equations by mapping γ to γ ↑ ? The motivation was given in Section 3.6. We need to type the modes as either long or transient, since the two situations require different index reduction techniques. This is precisely what our when-equations guarantee, since, by construction, switching from f to t occurs in a single nonstandard instant. Hence, by construction, when-equations are active at transient modes only. Actually, this distinction between if-equations and when-equations already exists in object-oriented DAE-based languages such as Modelica.

Regarding if-equations, we state the following assumption:

Assumption 2 In if-equations, guard γ holds true for a (standard) positive duration.

It is the responsibility of the programmer to make sure that if-equations are properly used so that Assumption 2 holds.

Long versus transient modes

For our theoretical development and the presentation of algorithms, it will be convenient to simplify our syntax by replacing when-equations by if-equations using (71) and Assumption 2. Our syntax of Definition 27 with arbitrary guards is thus replaced with the following (equivalent) syntax:

Definition 30 (mDAE/mdAE revisited) Multimode DAE systems (mDAEs), respectively multimode dAE systems (mdAEs), are defined by the following syntax:

S ::= s | S; s s ::= g | e g ::= γ = b e ::= if γ then f =0 ( 73
)
where b is a long/transient-typed Boolean expression of predicates over X ( ) (respectively X (•) ).

From now on, and unless otherwise specified, we will use the modified syntax of Definition 30. The notion of mode is then easily introduced.

Definition 31 (mode) For S a mDAE or mdAE, a mode is a valuation, in {f, t}, of all of its guards γ ∈ Γ(S). The mode is called transient if at least one guard typed transient takes the value t in it; otherwise, it is called long. Modes will be generically denoted by the symbol µ.

Inria

A mode enables a subset of the equations f =0 and disables the other ones. Nonstandard semantics. An mDAE S is transformed to a (nonstandard) mdAE S through the following syntactic transformation:

S = def S x is replaced with x • -x ∂ (74)
Since we perform the structural analysis on the nonstandard semantics, we must develop a structural analysis for mdAE systems, operating in discrete time. Therefore, in the rest of this section, we work in nonstandard semantics.

Comment 32 (Progressive evaluation of the modes) It may be that some guards have their value unknown when the execution of the current instant starts; this was, for example, the case with the Cup-and-Ball example in its first form given in Section 1.2.2. We can use the known enabled equations to evaluate part of the variables, hoping that some of the yet unknown guards get their value defined. This yields more equations getting enabled/disabled, allowing for more guards being evaluated. If this process leads to the evaluation of all guards, then our execution scheme succeeds; otherwise, we report a failure due to a logico-numerical fixpoint equation, as was the case in Section 1.2.2. This progressive evaluation of guards becomes problematic when both long and transient modes are considered. Indeed, different structural analyses are used for each case (see Section 3.3 for the long mode case of the Cup-and-Ball example, versus Section 3.5 for the transient mode case). Unfortunately, we do not know which one should be used if the mode is only partially known, since the long/transient qualification is only known when all guards have been evaluated. Consequently, we have to choose between the following two incompatible features for our compilation method:

1. either we support the progressive evaluation of modes, but are unable to handle transient modes;

2. or we support both long and transient modes, but we do not support the progressive evaluation of modes.

Since many physically meaningful models involve transient modes, our opinion is that the second option is much more useful in practice, hence it is adopted in our work.

Comment 6, regarding the Cup-and-Ball example, discussed what we miss by favoring the second option. Indeed, delaying the effect of guards was the appropriate correction to remove logiconumerical fixpoint problems in the Cup-and-Ball example. Remember that the so introduced delay is infinitesimal, so this does not really matter. 17Therefore, throughout the rest of this work, Definition 30 is modified to restrict ourselves by the syntax to systems in which all guards can be evaluated at the initialization of each nonstandard instant: 18Definition 33 (mdAE modified) Multimode dAE systems (mdAEs), are defined by the following syntax:

S ::= s | S; s s ::= g | e g ::= γ • = b (75) e ::= if γ then f =0
where b is a long/transient-typed Boolean expression of predicates over X (•) .

Notations 34 Referring to (75), we write b(γ) to refer to the expression defining γ • .

From now on, by mdAE, we will mean a system of the form given in Definition 33.

Structural analysis: intuition

In this section, we give an intuition about how to extend the structural analysis from dAE systems to mdAE systems. To simplify our explanation, we consider the restricted case in which both modes (previous and current) are long. Also, it will be convenient to assume, for the two modes before and after the change, that we have two index-1 models, of the form

S i : X i = F i (X i , Y i ) 0 = H i (X i ) , for i = 0, 1 ,
where indices 0 and 1 refer to the system before and after the change. Observe that system state X i may collect different variables before and after change, and so does the algebraic tuple Y i . Of course, the dynamics F i , H i are subject to changes as well. Let J i = def J X i H i denote the Jacobian of H i w.r.t. X i . Hence, the system augmented with its latent equations

S iΣ :    X i = F i (X i , Y i ) 0 = H i (X i ) 0 = J i (X i )X i
, for i = 0, 1 , seen as a system of algebraic equations with X i , Y i as dependent variables, is structurally nonsingular. Equivalently, we may replace, in this system, the differentiation of the algebraic constraint, d dt H i (X i ), by its shifting H i (X • i ) (their bipartite graphs are equal):

S iΣ :    X i = F i (X i , Y i ) 0 = H i (X i ) 0 = H i (X • i ) , for i = 0, 1 . (76)
Form (76) can be seen as the result of applying Pantelides' or Pryce's structural analysis in each mode, before and after change.

Comment 35 (key idea) Our key idea is to handle the mode change as a task of reconciling the possible conflict between the predictions generated by the previous mode and the consistency conditions associated to the new mode. Indeed, at the instant of mode change, there is a possible conflict between:

• on the one hand, the value for X 0 predicted by model S 0Σ which sets the value of the derivative X 0 just before the change;

• on the other hand, the consistency constraints on the pair X 1 , Y 1 resulting from model S 1Σ after the change.

To perform this reconciliation, we move to the nonstandard semantics by syntactically replacing, in (76), the derivatives x by their first-order Euler approximations with infinitesimal time step ∂:

x (t) ← is replaced by x(t + ∂) -x(t) ∂ = x • (t) -x(t) ∂ .
In addition, we replace, in (76), the derivative d dt H i (X i ) = J i (X i )X i by the time shift H i (X • i ). This substitution is legitimate, since the bipartite graph of the system is left unchanged by performing this replacement.

Performing the above substitutions in (76) and sampling it on the nonstandard discrete time basis T={0, ∂, 2∂, . . . } yields a discrete-time two-modes system. At the instant of mode change, we inherit the following system of algebraic equations:

S = • S 0Σ previous instant S 1Σ current instant (77) Inria where • S 0Σ :    X 0 -• X 0 ∂ = F 0 ( • X 0 , • Y 0 ) 0 = H 0 ( • X 0 ) 0 = H 0 (X 0 ) , S 1Σ :    X • 1 -X 1 ∂ = F 1 (X 1 , Y 1 ) (e 0 ) 0 = H 1 (X 1 ) (e 1 ) 0 = H 1 (X • 1 ) (e • 1 ) (78)
and ∪ refers to the union of sets of equations. For each state variable x that is shared between the two systems, i.e., x ∈ X 0 ∩ X 1 , the value for x must be identical in • S 0Σ and S 1Σ . To avoid duplicates of a same equation while taking the union of previous and current systems in S, we identify previous and current equations that are identical in their syntax, if any. The dependent variables of S are the X 0 inherited from S 0Σ and the X 1 , Y 1 , X • 1 inherited from S 1Σ . When seen as sets of variables, X 0 on the one hand, and X 1 , Y 1 , X • 1 on the other hand, possess in general a nonempty intersection. Being the result of a structural analysis, system S 1Σ is:

• structurally nonsingular, when seen as a system of algebraic equations in its leading variables Y 1 , X • 1 , for given values for X 1 ; and

• possibly structurally underdetermined, when seen as a system of algebraic equations in all its variables X 1 , Y 1 , X • 1 (ensuring consistent initialization).

Focus on S 1Σ in decomposition (77). Denote by Φ the subsystem of S 1Σ collecting all its equations that also belong to • S 0Σ , and are thus already solved.

The remaining system S 1Σ \ Φ must then be solved for its dependent variables, which consist of the subset of

X 1 ∪ Y 1 ∪ X •
1 collecting all variables not yet determined by • S 0Σ . To prepare for this, we submit S 1Σ \ Φ to the Dulmage-Mendelsohn (DM) decomposition, see Definition 15:

(β U , β E , β O ) ← DM(S 1Σ \ Φ). If β U = β O = ∅, then S 1Σ \ Φ is
structurally nonsingular by Lemma 16 and we solve it. Otherwise, the following can occur:

Case β U = ∅. This points to an underspecified subsystem at the instant of mode change. This information is returned to the designer for a correcting action.

Case β U = ∅ and β O = ∅. This points to a conflict between • S 0Σ and S 1Σ and is handled differently. We solve this conflict by applying Principle 1 of causality. Accordingly, the equations belonging to β O are removed from the consistency equations (e 1 ) of the new system S 1Σ . We thus keep the system S ↓ 1Σ = def {(e 0 ), (e 1 ) \ β O , (e • 1 )}. By Lemma 17, since the system {(e 0 ), (e • 1 )} is structurally nonsingular, performing the DM decomposition on system S ↓ 1Σ returns a triple with an empty overdetermined block. Hence, the so obtained reduced system S ↓ 1Σ is structurally regular and ready to be solved.

Comment 36 If the two systems S iΣ , i = 0, 1 are identical (there is no mode change), then X 1 = X 0 , Y 1 = Y 0 , F 1 = F 0 , and H 1 = H 0 , so that Φ consists of the consistency subsystem 0 = H 1 (X 1 ), see the definition of S 1Σ in (78). Hence, S 1Σ \ Φ coincides in this case with the structurally nonsingular system obtained by performing index reduction (see Definition 24), and we are left with the usual structural analysis for DAE.

The reconciliation may take several nonstandard instants following the mode change, where less and less equations from S 1Σ are erased when applying DM, until we reach the steady regime. The interpretation of the erasure of some equations from S 1Σ is that their satisfaction gets postponed for a few instants. Recall that we are in nonstandard semantics, so that all of this will occur in a total zero duration, in standard real time.

In the rest of this section, we formalize this intuition and generalize it.

Structural Analysis of long modes

In this section we develop the structural analysis of mdAE systems for the subclass of systems possessing only long modes. Within each long mode, we simply reuse the existing structural analysis based on the Σ-method. The only difficulty sits in the mode changes, for which we use conflict reconciliation as explained above. Throughout this section, we work in the nonstandard semantics as given by Equation (74).

Constructive Semantics

We reuse the techniques of constructive semantics, first introduced in the context of reactive synchronous programming languages [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF][START_REF] Benveniste | Compositionality in dataflow synchronous languages: Specification and distributed code generation[END_REF][START_REF] Benveniste | The synchronous languages 12 years later[END_REF][START_REF] Potop-Butucaru | Compiling Esterel[END_REF], for the purpose of grounding compilation on solid mathematical foundations. For synchronous languages, the execution of a program proceeds through successive reactions, by which discrete time progresses from the current instant to the next one. The compilation of a program consists in generating the code for executing a reaction. This task is formalized through the notion of constructive semantics, which relies on the following two pillars:

1. The specification of the set of atomic actions, which are effective, non-interruptible, state transformations. For synchronous languages, atomic actions consist of: (a) the evaluation of a single expression, and (b) control flow actions. Executing an instance of an atomic action is referred to as performing a micro-step.

2. The causality analysis of the reaction, which is a partial order on the set of micro-steps, abstracting the dependencies between them. A scheduling of the micro-steps is correct if and only if it complies with the causality analysis.

The constructive semantics consists in decomposing a reaction as the symbolic execution of a sequence of micro-steps subject to causality constraints. Intermediate stages of this execution are represented by the 'status' of each variable, belonging to the abstract alphabet {not evaluated/evaluated}. The assignment of a status to all variables is called a configuration. Such executions are monotonic, in that the so produced sequence of configurations is increasing with respect to the product order derived from the order not evaluated < evaluated for each variable. The constructive semantics is complete if it terminates by having all variables evaluated. Failure occurs if an incomplete sequence cannot be extended. This approach is suited to mdAE systems, since it will provide the formal support for correctly chaining the different actions composing the structural analysis, and describing the resulting increase in knowledge on all involved equations, variables and guards. For mdAE systems, the atomic actions, however, consist of: (a) the evaluation of a guard;

(b) manipulations of systems of equations: adding latent equations or deleting conflicting equations at mode changes; and

(c) solving a block of algebraic equations.

Actions (a) and (b) participate in the structural analysis, whereas (c) is the runtime solving of systems of equations to evaluate the leading variables. To keep this evaluation symbolic, we abstract away the values of numerical variables and regard the values of guards as oracles, so that all the possible valuations must be explored.

Statuses and Contexts

The different statuses of the guards, variables and equations of an mdAE are the following: an evaluated guard may be true or false; a variable may be not evaluated, or evaluated; an equation may be not evaluated, disabled (its guard is false), or solved. This is summed up in Table 1. There is no need to consider the status not evaluated for a guard, since guards can be evaluated at the initialization of each nonstandard instant, by Definition 33. Unlike for single-mode DAEs, the set of equations defining the dynamics of an mdAE is modedependent. This information should therefore be traced as part of the constructive semantics of mdAE systems. To capture this, we tag as irrelevant in the considered mode all the variables, Inria i u f t guard irrelevant evaluated to f evaluated to t variable irrelevant not evaluated evaluated equation irrelevant not evaluated disabled solved Table 1: The partially ordered domain D = {i, u, f, t} and its interpretation for guards, variables, and equations.

guards, and equations that are not used to define the dynamics of the mdAE system in this mode. Formally, the domain of statuses is the partially ordered set (D, <) defined by

D = {i, u, f, t} and i < u < f, t . (79)
The interpretation of this domain for variables, equations and guards, is summarized in Table 1 and detailed next:

• The minimal element i is used to represent the fact that a variable, a guard, or an equation is irrelevant.

• Value u means that a variable or equation has not been evaluated yet. At the beginning of a time step, only state variables are known, and all other variables are set to u, reflecting that their numerical values are not known yet.

• Maximal element t has different meanings, depending on whether it applies to a variable, a guard or an equation. In the case of a variable, it means that the numerical value of the variable has been computed, whatever it could be. For a guard, it means that the guard has been evaluated to true. For an equation, it means that the equation has been solved.

• Maximal element f also has different meanings, depending on whether it applies to a guard or an equation. In the context of a guard, it means that the guard has been evaluated to false. When it applies to an equation, it means that the equation is disabled, i.e., its guard is false. This value does not apply to variables.

We call status the assignment, to each guard, variable, and equation, of a value over D. It is formally defined as follows:

Definition 37 (status) Let S be an mdAE system. Its set V of S-variables 19 is

V = def Γ ∪ X (•) ∪ E (•) (80)
(using Notations 26 and 28). A status σ of S is a valuation in D of the S-variables, that is, a mapping σ : V → D. A status is called finite if it is equal to i, except for a finite number of S-variables. The set of statuses is partially ordered by the product order:

σ 1 ≤ σ 2 if and only if for all v ∈ V, σ 1 (v) ≤ σ 2 (v).
We denote by i u, u f, and u t the increasing changes of status for a given S-variable. Enabled Equations and Leading Variables. Let σ be a status. Equation e •m is enabled in σ (respectively disabled in σ) if and only if σ(γ) = t (respectively σ(γ) = f), where γ is the guard of e •m , see (70). Denote by Enab(σ) (respectively Disab(σ))

the set of equations that are enabled (respectively disabled) in σ. Recall that, for any finite status σ, these sets are finite. The set Enab(σ) defines the DAE system that is enabled at status σ. The leading variables in status σ are defined, with reference to this enabled DAE system, following the introduction of Section 6.2. We define the set

Undef (σ) = def {v ∈ V | σ(v) ≤ u}
collecting the S-variables that are either irrelevant or not evaluated in status σ. We are now ready to introduce micro-steps, and runs as finite sequences of micro-steps.

Definition 38 (micro-step) A micro-step is a transformation, via some atomic action, of a status σ to a status σ by updating the values of a finite subset of S-variables, from i to u, or from u to t or f.

A run is a finite sequence of micro-steps. It is complete when every equation is either irrelevant, disabled, or solved, and all the enabled leading variables are evaluated. Formally (the reader is referred to Definition 37 for the notions used):

Definition 39 (runs) Given a finite initial status σ 0 , a run of mdAE system S is a finite increasing sequence of statuses

σ 0 < σ 1 < • • • < σ k < σ k+1 < • • • < σ end (81)
such that, for every k < end, σ k+1 is obtained from σ k by performing a micro-step. The final status σ end is called complete if, in σ end , all equations e i have either the value t or f and no leading variable has the value u. A run σ 0 < . . . < σ end of mdAE system S is complete if and only if its final status σ end is complete.

Handling conflicts between past and present. To handle the conflicts between the previous and current instants when a mode change occurs, we will need to keep track of the possible sources of conflict between previous and current instants. If equation e •(m+1) , m≥0, was solved at the previous instant (meaning that σ(e •(m+1) ) = t for the final status of that instant), then e •m is a possible source of conflict at the current instant. The set of all such equations is called the context at the current instant. All the variables involved in the context are already evaluated at the initial status of the current run (value t).

If, however, equation e •m ∈ ∆ is also enabled at status σ of the current instant, then e •m is no longer a source of conflict, but rather a satisfied consistency equation. We call it a fact at σ. Formally: Definition 40 (contexts and facts) Let S be an mdAE system. Its context ∆ at the current instant, and set of facts Φ(σ) at the current status, are defined as follows:

∆ = e •m | the final status of e •(m+1) at previous instant is t Φ(σ) = ∆ ∩ Enab(σ)
By abuse of notation, we will also denote by ∆ and Φ the sets of functions f (e), for e ranging over ∆ and Φ.

To prepare for the algorithm computing a run for given initial status σ 0 and context ∆, we now collect a few useful auxiliary algorithms.

Auxiliary algorithms

Tick algorithm (Algorithm 7) When a run is complete, the system can proceed to the next time step by executing the algorithm Tick. Tick requires a complete (final) status σ and returns a pair (σ • , ∆ • ) consisting of an initial status σ • and a context ∆ • , for use as initial conditions for the next run.

Performing Tick has the effect of shifting backward defined variables, and setting all other S-variables v ∈ V to either u, if v is relevant to the mdAE S, or i (irrelevant) otherwise. Guards are set either to i or to the value σ(b(γ)) ∈ {f, t}, using Notations 34 and the definition (75) for γ • . Thus, the value of all guards at the next instant is set. The new context is defined to be the set of equations that are known to be satisfied in the next instant. Tick is not increasing and it does not have to be, since it applies when moving to the next time step.

Inria

Algorithm 7 Tick

Require: σ, a complete status; return pair (σ • , ∆ • ).

1: Tick(σ) = def (σ • , ∆ • ), where:

σ • (γ) = if γ is involved in S then σ(b(γ)) else i σ • x •m = if σ x •m+1 = t then t else if x •m is a variable of S then u else i σ • e •m = if e •m is a variable of S then u else i ∆ • = e •m e∈E, m∈N : σ e •m+1 =t
IndexReduc algorithm (Algorithm 8) Given a system G with its set of leading variables determined by the status σ, the IndexReduc algorithm is a renaming of the Σ-method (Algorithm 6), together with a refined interpretation of its result. It is understood that we use the discrete-time adaptation of it, in which forward shift replaces differentiation but everything else remains unchanged.

Algorithm 8 IndexReduc

Require:

(G, σ); return (σ, b, G Σ , G Σ ) 1: (b, G Σ , G Σ ) ← SigmaMethod(G) with dependent variables set by σ; increase σ 2: if b then return (G Σ , G Σ )
Recall that this algorithm requires a square system G=0 and returns a success/failure flag b, as well as, in case of success, a pair consisting of G Σ =0, a structurally nonsingular system determining some leading variables, and the system G Σ =0 collecting the consistency equations. As a consequence, the status of every latent equation added by the Σ-method is updated as i u.

SolveConflict algorithm (Algorithm 9) This algorithm requires a system K=0 and a status σ for all variables involved in K. System K is submitted to the DM decomposition (Line 1); in doing so, the dependent variables are those with status u. Possible conflicts are collected in β O . If the latter is nonempty, corresponding equations are deleted (Line 3). The status is updated as u f for every removed equation. The so reduced system is again submitted to the DM decomposition (Line 4), which, by Lemma 17, will return an empty overdetermined block: conflicts are removed. If block β U is nonempty, then system K is underdetermined. Therefore, b = f is returned (Line 5). Otherwise, b = t and H = β E are returned (Line 6).

Algorithm 9 SolveConflict

Require: (K, σ); return (σ, b, H) 1: (β U , β E , β O ) ← DM(K) with dependent variables set by σ 2: if β O = ∅ then 3: K ← K \ β O ; increase σ 4:
(β U , β E , ∅) ← DM(K) with dependent variables set by σ

5: b ← [β U = ∅] 6: if b then H ← β E Eval algorithm σ ← Eval(H, σ) (82)
This algorithm requires a status σ and a structurally regular system H, enabled at σ. It evaluates H for its dependent variables, which sets the values of the leading variables. The corresponding change in status is u t for both the variables that are evaluated and the equations that are solved.

Executing a nonstandard instant

Algorithm 10 ExecRun Require: (σ, ∆); return (Fail , (σ

• , ∆ • )) 1: Φ ← {e ∈ ∆ | σ(γ(e)) = t} 2: G ← [Enab(σ) ∩ Undef (σ)] 3: (σ, b, G Σ , G Σ ) ← IndexReduc(G, σ); increase σ 4: if ¬b then return Fail (σ) 5: else 6: (σ, b, H) ← SolveConflict((G Σ ∪ G Σ ) \ Φ, σ); increase σ 7:
if ¬b then return Fail (σ)

8: else 9:
σ ← Eval(σ, H); increase σ 10:

(σ • , ∆ • ) ← Tick(σ)
ExecRun requires a finite status σ and a finite context ∆ to serve as initial conditions for the run. It returns, either a documented Fail, or a final value for status σ, as well as values for σ • and ∆ • to serve as initial conditions for the next instant. A line-by-line description of ExecRun is given next.

Line 1 Not all the equations belonging to the context are active in the current instant. The active ones are collected in the set Φ of facts.

Line 2 System G is set to the enabled guarded equations that are not evaluated yet in status σ.

Line 3 We submit G to index reduction (IndexReduc algorithm, Algorithm 8).

Line 4 If IndexReduc returns b=f, then ExecRun returns Fail (σ) and stops.

Line 5 Otherwise, ExecRun returns a pair consisting of G Σ =0, a structurally nonsingular system determining some leading variables, and the system G Σ =0 collecting consistency equations. In this case, the algorithm can further progress.

Line 6

Taking the set Φ of facts into account, the possible conflict is solved using the SolveConflict algorithm (Algorithm 9). When applying SolveConflict, the dependent variables of (G Σ ∪G Σ )\Φ include all the variables (both leading and non-leading) of the current system that have status σ = u, i.e., whose value was not set by executing the previous nonstandard instant.

Line 7 If SolveConflict fails, ExecRun returns Fail (σ) and stops.

Line 9 Otherwise, SolveConflict returns a structurally regular system H, which is solved using Eval, and the so reached status σ is complete.

Line 10 Since σ is complete, Tick (Algorithm 7) is performed, so that the values of guards, the status of all S-variables and the context are known for the next nonstandard instant.

Important properties

Evaluating guards Since ExecRun (Algorithm 10) performs a symbolic interpretation of the model, all the possible modes (valuations for the guards) must be hypothesized and explored at

Inria Line 10. 20 The program is accepted if and only if, for all possible modes, no Fail is returned by ExecRun. In case of success, all reachable pairs (σ, ∆) of statuses and contexts have been explored.

ω 1 , ω 2 start γ, ω 1 , ω 2 , e 3 , e 4 γ, ω 1 , ω 2 , τ 1 , τ 2 , ω • 1 , ω • 2 , e ∂ 1 ,
ExecRun as a labeled automaton We can see ExecRun as a labeled finite state automaton having nodes labeled with the different encountered status-context pairs, and transitions labeled with the micro-steps. In Figure 20, we depict this graph for the clutch example. The Tick transitions indicate a move to the next nonstandard instant.

Rejecting models Firstly, a system exhibiting logico-numerical fixpoint equations is rejected via synctactic checks, since it does not agree with Definition 33. Then, a failure is found in the following cases:

• SigmaMethod (Algorithm 6) itself fails, indicating that no matching is found for the bipartite graph G G of the running dAE system G=0, a failure criterion for the Σ-method. Failure message indicates either over-or under-determination.

• The conflict between the predictions from the past and the consistency conditions from the current mode cannot be solved, i.e., β U = ∅ at Line 5 of SolveConflict, revealing some missing restart equation(s) for the new mode. Solving conflicts at mode changes takes only finitely many nonstandard instants This follows from the fact that, at a mode change, the maximal shifting degree of context ∆ is finite and the maximal shifting degree of the conflicting subsystem decreases by one at each subsequent time step. See Section 3.3 on the Cup-and-Ball example as an illustration of this.

Detecting continuous modes

Determinism of the algorithm All the algorithms called throughout the execution of Algorithm 10 possess a unique return (at Line 3 and Line 6). Since Algorithm 10 involves no choice by itself, it is deterministic. Comment 18 following Lemma 17 finds its justification here: our design choice in Lemma 17 was essential in guaranteeing the determinism of our compilation method.

Generic Blocks within a long mode

The results of this section will be used in Section 10, and more specifically in Section 10.4. The successive application of Line 3 and Line 6 of Algorithm 10 yields blocks of a specific form, that is investigated next.

Notations 41 For f : R K → R a smooth numerical function and n ≥ 0 an integer, set

B • (f, n) = def B h • (f, n) B t • (f, n)
where

B h • (f, n) = def      f f • . . . f •(n-1)      and B t • (f, n) = def f •n . ( 83
) B • (f, n) is called an (f, n)-block, B h • (f, n) is its head and B t • (f, n) is its tail. By convention, B • (f, n) is empty for n < 0 and B h • (f, n) is empty for n = 0. For {f 1 , . . . , f k , . . . , f K } a tuple of K smooth functions in K dependent variables, define F =    f 1 . . . f K   , N =    n 1 . . . n K   , B • (F, N )=    B • (f 1 , n 1 ) . . . B • (f K , n K )    ( 84
)
and the definitions of B h • (F, N ) and B t • (F, N ) follow accordingly. Lemma 42 Assume that the system is within a long mode. Then, successively applying Line 3 and Line 6 of Algorithm 10 to the system G=0 returns, either 'Fail', or a structurally nonsingular system H=B t

• (G, N ) of the form (103).

Proof Consider Line 3 of Algorithm 10. Since the system is in a long mode, SigmaMethod returns, either Fail, or (G

Σ , ḠΣ ) = (B t • (G, N ), B h • (G, N ))
, where G is the subset of enabled equations at status σ.

Move to Line 6. Since the system is in a long mode, we have ḠΣ =Φ, expressing that consistency conditions are guaranteed by previous nonstandard instants. Hence, the Dulmage-Mendelsohn decomposition performed at Line 6 simply returns β E =G Σ =B t

• (G, N ).

Inria

Structural Analysis: general case

So far, we restricted ourselves to systems having only long modes. In this section, we consider the general case, where transient modes can also occur, possibly in cascades. The Cup-and-Ball example of Section 3 can be kept in mind while reading this section. Our main task is to revisit index reduction. So far, this was performed by invoking the Σ-method, which requires the new mode to last for long enough to allow for an unlimited number of time shifts to be applied. If the current mode is transient, this no longer holds and the continuation of the system dynamics is determined by its successive future modes-see the analysis of the Cup-and-Ball example in Section 3.5.

Whereas the Σ-method applies to time-invariant systems only, the method of differentiation arrays (see Section 3.1 of [START_REF] Campbell | The index of general nonlinear DAEs[END_REF]) can be adapted to time-varying systems.

Index reduction with Difference Arrays

The index reduction algorithm proposed in this section applies to time-varying DAE systems. It is, however, less computationally efficient than the Σ-method (when the latter can be applied).

Difference arrays for time-varying systems

In the following, T denotes a discrete-time index set possessing a minimal element, denoted by 0.

Definition 43 A time-varying dAE system over T is a pair (X, F), where X is a tuple of numerical variables taking its values in an Euclidian space D, and F : T → E maps every instant t ∈ T to a dAE system F t = 0 whose set of dependent variables is contained in X. A solution of (X, F) is a trajectory x : T → D such that, for every t ∈ T, x(t) satisfies F t = 0.

Example 44 Consider T = N, X = (u, v, w), and F is defined by t=0 : 0 = v • -3u ; and ∀t > 0 :

   0 = u • -v 0 = v -uw 0 = w • + u .
The solutions of this system are trajectories t → (u(t), v(t), w(t)) such that 0 = v(1) -3u(0) ; ∀t > 0 :

   0 = u(t + 1) -v(t) 0 = v(t) -u(t)w(t) 0 = w(t + 1) + u(t) .
Note that solutions are only partially specified at t = 0 since w is not involved in the initial conditions.

Given a time-varying dAE system (X, F), we consider, for each k ∈ T, the following dAE array:

A k = def      F 0 F 1 . . . F k      . ( 85
)
This array allows us to perform the 'index reduction' of F 0 , by adding 'latent equations' taken from the future systems F 1 , F 2 , etc., generally different from F • 0 , F •2 0 , etc. To this end, the set of all variables of array A k decomposes as Y ∪ X ∪ W , where:

Y collects the variables in F 0 that participate in the consistency; X collects the leading variables of F 0 ; W collects all the remaining variables of the array.

(86)

RR n°9334

In its Section 3.1 on 'standard indices', Campbell and Gear's landmark paper [START_REF] Campbell | The index of general nonlinear DAEs[END_REF] states that, for time-invariant systems (F t ≡ F 0 ), the index is the smallest k such that ∃W : A k (X, W, Y )=0 uniquely defines X as a function of Y , provided that Y is consistent.

Therefore, we consider the following structural translation of the above problem (compare with Problem 23 in Section 6.2):

Problem 45 Find the smallest integer k such that the system ∃W : A k (X, W, Y )=0 is structurally nonsingular in the sense of Definition 21.

For time-invariant systems (F t ≡ F for all t), Problem 45 is monotonic in that, if this problem has a solution for k, then it also has a solution for k > k. However, this no longer holds for time-varying systems, since it can be that system F k+1 = 0 is in conflict with A k .

Handling finite cascades of transient modes

Now, we specialize the structural analysis of difference arrays to the handling of finite cascades of transient modes. That is, we assume the following finite chain of modes:

µ 0 , µ 1 . . . , µ K finite cascade of transient modes , µ ∞ final long mode (87)
We adopt the following principles in developing our structural analysis:

• For the final long mode µ ∞ , it is legitimate to postpone the satisfaction of some of the consistency equations. That is, we reuse the conflict solving approach of the algorithm SolveConflict (Algorithm 9).

• For the transient modes µ 0 , µ 1 . . . , µ K , however, this is not legitimate, so that we require the corresponding equations to be satisfied and refuse to erase any of them.

We associate to the chain (87) the array

A = def        F 0 F 1 . . . F K F ∞,Σ ∪ F ∞,Σ       
with a decomposition of its variables as Y ∪ X ∪ W following (86).

(88)

The last row of this array originates from applying index reduction on the long mode µ ∞ .

Based on the results of Section 6.1.4, in particular Lemma 22, the structural analysis for timevarying systems is described in Algorithm 11 (DiffArray). This algorithm requires as inputs an array A as in (88) and a status; it returns, either a Fail information (b = f), or a pair (F Σ , F Σ ), where F Σ =0 is a structurally nonsingular system extracted from A and determining the leading variables of F 0 , and F Σ =0 is a solvable system of consistency equations. The algorithm proceeds according to two phases.

In a first phase (Line 1-Line 10), the subchain of transient modes is explored, with increasing array size k = 0, . . . , K. Using status σ, the variables of A k are sorted at Line 5 following (86). At Line 6, array A k is submitted to ExistQuantifEqn (Algorithm 5), which returns the two Boolean flags b O , indicating whether the overdetermined block is empty, and b U , indicating whether the underdetermined block is empty. If b O = f, which indicates an overdetermined system, a failure indication is returned at Line 7. If b O = t and b U = f is returned, the array is increased (Line 3) and reprocessed the same way, until k = K is reached. Otherwise, the first phase returns the pair (F Σ , F Σ ). The status is updated as i u for each latent equation added, and the algorithm successfully terminates.

Inria

Algorithm 11 DiffArray Require: (A, σ) as in (88); return (σ, b, F Σ , F Σ )

1: k = -1; 2: if k < K then 3: k ← k+1 4: A k ← (F 0 , . . . , F k ) 5:
arrange the variables of

A k as (Y k ∪X k ∪W k ) 6: (b O , b U , F Σ , F Σ ) ← ExistQuantifEqn(A k ) 7: if ¬b O then return b ← f 8: else 9:
if ¬b U then go to Line 2

A ← A where F ∞,Σ ← (F ∞,Σ \ A K ) 13: (β U , β E , β O ) = DM(A) 14: if β O = ∅ then 15: A ← A where F ∞,Σ ← F ∞,Σ \ β O ; 16:
go to Line 13 partition

β E = F Σ ∪ F Σ 20: return (b ← t, F Σ , F Σ ); increase σ 21:
else return b ← f

If k=K is reached with no success, we move to the second phase, starting with the removal of duplicates between A K and F ∞,Σ (Line 12). Then, we apply the DM decomposition at Line 13 to identify possible conflicts and underdetermined parts. We handle conflicts as in the SolveConflict algorithm (Algorithm 9), by deleting conflicting equations from F ∞,Σ (Line 15) and returning to Line 13; this loop is visited only once thanks to Lemma 17. When reaching Line 17, we know that condition 1 of Lemma 22 is satisfied, so that we can conclude by using the successful case of Lemma 22 at Lines 18 and subsequent. The status is updated as i u for each latent equation added (Line 20).

Executing a nonstandard instant

We first revisit the algorithm IndexReduc. If the current mode is long, then IndexReduc is implemented as Algorithm 8. Otherwise, we will implement it by hypothesizing a future mode trajectory and calling the algorithm DiffArray (Algorithm 11). IndexReduc, revisited, is specified as Algorithm 12.

Algorithm 12 IndexReduc, revisited Require:

(G, σ, type); return (σ, b, G Σ , G Σ ) 1: if type=long then 2:
(b, G Σ , G Σ ) ← SigmaMethod(G) with dependent variables set by σ; increase σ 3: else 4:

hypothesize F = (G 0 , G 1 , . . . ) with G 0 =G; let A be the corresponding array

5: (b, G Σ , G Σ ) ← DiffArray(A, σ); increase σ 6: if b then return (G Σ , G Σ )
In this algorithm, Line 4 requires hypothesizing a continuation for the current mode. This, of course, is problematic if all modes can be visited in this continuation, particularly transient modes.

It is therefore essential to be able, at compile time, to restrict the number of possible continuations for the current mode. This question is investigated in Section 7.4.3.

Algorithm 13 ExecRun , revisited Require: (σ, ∆); return (Fail , σ, (σ

• , ∆ • )) 1: type ← EvalModeType(σ) 2: Φ ← {e ∈ ∆ | σ(γ(e)) = t} 3: G ← [Enab(σ) ∩ Undef (σ)] 4: (σ, b, G Σ , G Σ ) ← IndexReduc(G,
(σ, b, H) ← SolveConflict((G Σ ∪ G Σ ) \ Φ, σ); increase σ 8:
if ¬b then return Fail (σ)

(σ • , ∆ • ) ← Tick(σ)
Having revisited index reduction, we can now revisit ExecRun, see Algorithm 13. Only two lines are modified with reference to the original algorithm ExecRun: Line 1 We add the evaluation of the type (long or transient) of the current mode.

Line 4

We call instead the IndexReduc corresponding to Algorithm 12, which selects the appropriate method for index reduction, depending on the type of the current mode.

Continuations of transient modes

We need to investigate this issue for transient modes only. Therefore, we begin with assumptions regarding transient modes. Recall that we follow Notations 26, hence, by numerical variable of an mDAE system S, we mean a variable of the form x m (m-th derivative of x) for some algebraic or state variable x of the considered system. The degree of x m is m. The same holds regarding mdAE system S, with x •m (m-shifted version of x) replacing x m . We repeatedly use these notations in the sequel.

Assumptions regarding transient modes According to Definition 31, a mode µ is transient if at least one of its t-valued guards is transient. Guards are Boolean expressions of predicates over numerical variables. We will restrict ourselves to the class of nonstandard mdAE models S with set X of numerical variables, satisfying the following assumption regarding transient predicates over numerical variables:

Assumption 3 Let γ 1 , γ 2 ∈ Γ(S)
be any two distinct transient predicates over numerical variables of mdAE model S. Then, for ν 1 X and ν 2 X any two valuations of the tuple X such that ν

1 X -ν 2 X is infinitesimal, γ 1 (ν 1 X ) = γ 2 (ν 2 X ) = t cannot occur.
Discussion. Assumption 3 formalizes that two different transient predicates never take the value t simultaneously, "by chance". For example, for x, y ∈ X, the two predicates γ 1 = when x≥0 and γ 2 = when y≤5 specify the events when x becomes ≥0 and when y becomes ≤5. The assumption says that these two events will never occur simultaneously. For this case, Assumption 3 is very reasonable. Of course, if γ 1 = when x 3 ≥8 and γ 2 = when x≥2, the two guards are considered different (their syntax differ) but generate simultaneous events. It is not the duty of the compiler to discover if two guards happen to be mathematically identical whereas they differ in their syntax. It is rather the responsibility of the programmer to use identical syntax for different occurrences of a same predicate.

Inria

Using Assumption 3 we can in particular bound the number of transient modes that can occur at a given instant t ∈ T. For γ o ∈ Γ(S) a transient predicate, let

Γ(S, γ o ) ⊆ Γ(S) (89)
be the subset of all guards γ of S that are Boolean expressions of predicates such that γ = t requires γ o = t. An example would be γ = γ o ∧ γ . The following obvious lemma holds:

Lemma 46 Under Assumption 3, the number of transient modes that can possibly occur at a given instant t ∈ T is bounded by ,γo) where γ o ranges over the set of all transient predicate guards of S, and |A| denotes the cardinal of set A.

M (S) = max γo 2 Γ(S
In practice, this bound is generally much smaller than the number of possible modes for S, which is bounded by |2 Γ(S) |.

Assumptions regarding continuations of transient modes Let µ 0 = µ(t) be the current mode at time t ∈ T, and let µ 0 , µ 1 , µ 2 , . . . be any continuation for µ 0 at instants t, t

finite cascade : µ 0 , µ 1 , . . . , µ m-1 , µ ∞ (90) infinite cascade : µ 0 , µ 1 , . . . , µ m-1 , . . . (91)
where the final mode µ ∞ in (90) is long (and thus repeated), and all other listed modes, including the ". . . ", are transient. We will state assumptions that statically rule out (91) and limit the number of possible continuations in (90). For µ a transient mode, let X(µ) ⊆ X be the subset of numerical variables x such that x • -x is possibly non-infinitesimal at this mode.

In other words, X(µ) is the set of variables that are reset by this transient mode. We then define, for M a set of transient modes:

X(M ) = µ∈M X(µ) .
This set can be identified using structural algorithms.

Example 47

The statement

x = f (x, u) reset 0 when x ≥ 1
expands, according to Definition 27, as

when x ≥ 1 then x + = 0 else x = f (x, u)
whose nonstandard semantics is

when x ≥ 1 then x • = 0 (e 1 ) else x • -x ∂ = f (x, u) (e 2 ) ( 92
)
Since x • -x is not guaranteed infinitesimal by (e 1 ), x • -x is non-infinitesimal at the instant of zero-crossing (i.e., the first time when x ≥ 1 occurs). Thus, X(µ zc ) = {x}, where µ zc is the mode when the zero-crossing occurs. This can be found at compile time, based on the syntax of System (92).

Example 48 Consider the Cup-and-Ball example with elastic impact, see Section 3.5, in its original specification prior to adding the restart condition [START_REF] Pothen | Computing the block triangular form of a sparse matrix[END_REF]. Let µ imp be the mode when the elastic impact occurs. Then, X(µ imp ) = ∅. Thus, we know that the only continuation is: µ imp , µ free , in which µ free is the mode in which the rope is not straight and, thus, free motion occurs. The same holds after adding restart condition [START_REF] Pothen | Computing the block triangular form of a sparse matrix[END_REF].

In the sequel, we denote by M the set of all possible modes of system S. For Y ⊂ X a subset of variables, let M (Y ) ⊆ M be the set of all transient modes that can possibly be visited as a consequence of making noninfinitesimal changes to the variables of Y -this set can be identified at compile time using structural algorithms.

In bounding the possible continuations of a transient mode µ 0 , we consider the following chain of sets of possibly reachable transient modes:

Continuations (µ 0 ) = M 0 , M 1 , M 2 , . . .
where :

M 0 = {µ 0 } and ∀k > 0 : M k = M (X(M k-1 ))
Once again thanks to structural algorithms, the following assumption can be checked at compile time on the considered system S: Assumption 4 For every transient mode µ, the chain X k , k ≥ 0 in Continuations (µ) converges to ∅ in finitely many steps. Assumption 5 holds if, when a transient mode is reached due to a predicate on positions x, y, . . . , then only velocities x , y , . . . or accelerations x , y , . . . can be reset. This kind of property is often encountered in classes of physical systems, e.g., contact mechanics.

Systems with shifts and differentiations

So far, our formalism of multi-mode DAE systems handles only guarded DAEs. For more elaborate examples, we need to support also left-and right-limits, to be able to explicitly define restart or reset values for the dynamics of a given mode. With reference to model class (1), our extension takes the form if γ j (the x i their derivatives and left-limits) do f j (the x i their derivatives and right-limits) = 0 (93)

Right-limits are typically used in defining restart values at mode changes. Since we will map the dynamics in the nonstandard domain, we will represent the left-and right-limits x -and x + through the backward and forward shifts • x and x • , see Definition 2. We thus consider a syntax offering derivative and forward shift •, both related by the following identities in the nonstandard semantics:

x = 1 ∂ (x • -x) x • = x + ∂x (94)
The resulting formalism will be called mdDAE, or ∂AE, for short. We formalize this next.

Inria

Preliminaries

In the nonstandard semantics, shift and derivative are related via identities (94). In these preliminaries we investigate additional relations and properties that follow from identities (94).

Lemma 50

1. The two operators • and commute: for every variable x, we have

x • = x • .
2. For f : R k → R a C 1 -function and X a vector variable of dimension k, define

f • (X) = def [f (X)] • = def f (X • ) (f (X)) = def ∇f (X).X ( 95
)
where ∇f denotes the Jacobian of f . Then, differentiation and shift commute, meaning that:

[f (X)] • = [f (X)] • .
Proof Statement 1 follows immediately from identities (94). Then, the equalities

[f (X)] • = [f • (X)] = ∇f (X • ).X • (by Statement 1) = [∇f (X).X ] • = [f (X)] • prove Statement 2.
We consider the alphabet {•, } related through identities (94), representing the shift and the derivative operators, respectively. Successive differentiations or shifts of a variable x are denoted by x w for some word w ∈ {•, } * where * denotes the usual Kleene closure. We denote by •n and n the concatenation of n successive • and , respectively. In x w , word w is read from left to right, e.g., in x • , x is first differentiated, and then x is shifted. We denote by X w the set of all x w for x ranging over the set X of variables. For Σ ⊆ {•, } * , we set

X (Σ) = def w∈Σ x w and {X} = def X ({•, } * ) (96)
Definition 51 A ∂AE system is a finite set of guarded equations of the form:

e ::= if γ do f =0 (97) γ • ::= bexp ( 98
)
where f is a smooth scalar function over {X}, the shifted guard γ • is defined by bexp, which is a typed boolean expression of predicates over {X}, and the guard γ is typed long/transient.

The reader is invited to compare this definition with Definition 33. The notations of Definition 33 and the explanations that follow it extend to ∂AE. In a ∂AE model we may find x, x , x • , x • , x • and so on. Using (94) and Lemma 50 repeatedly, x • can be expressed as a function of x, x • , x •2 , x •3 or as a function of x, x , x , x 3 . Consequently, for w ∈ { , •} * a word over the alphabet { , •}, it is natural to define the degree of x w with respect to x as being |w|, the length of word w. For X, Σ ⊆ {•, } * , and X (Σ) as in (96), let E be a system of algebraic equations over the set X (Σ) .

We denote by (•) E the system obtained by expressing any monomial x w as a function of x, x • , x •2 , x •3 , . . . , by repeatedly using the first identity of (94) and Lemma 50.

System (•) E involves only shifts and no differentiations. Vice-versa, recalling (95):

For E of the shifted form E = F •K , let ( ) E be obtained by repeatedly applying, to E, until all shifts get eliminated from it, the map G(X • ) → G(X) + ∇G(X).X , where ∇G denotes the Jacobian of G.

(100)

Note the difference in defining (•) E and ( ) E.

Shifting versus differentiating

The whole Sections 7 and 7.3, where the structural analysis of mdAE systems was developed, extends verbatim to ∂AE systems, up to the following changes:

• Replace everywhere the superscript •m, for m∈N, by a word w∈{ , •} * of length |w|=m; define, for e w :: if γ do f w =0, the operators

ForwardShift : e w• = def if γ do f w• = 0 Differentiate : e w = def if γ do f w = 0 (101)
• When invoking the Σ-method (via IndexReduc) in Algorithm 10 (respectively Algorithm 13), we can either use differentiation or forward shift at Line 3 (respectively Line 4). The question is: does this choice influence the outcome of the Dulmage-Mendelsohn decomposition at Line 1 of SolveConflict (Algorithm 9)? Other steps of Algorithm 10 or Algorithm 13 are for sure not affected.

The alternative "shift or differentiate" raises the following question: what is the effect, on the structural analysis, of chosing a particular alternative each time we reach Line 3 of Algorithm 10 or Line 4 of Algorithm 13 ? We address this issue next.

Invariance results

We use the notion of weighted bipartite graph associated to a ∂AE system F =0. We refer the reader to Section 6.2 for this notion associated with DAE systems, and we extend it to ∂AE systems by deducing weights from the length |w| of words over the alphabet {•, } * . Thanks to the two-way correspondence (94), the following holds:

Lemma 52 For F =0 any system of algebraic equations over X (Σ) : F , (•) F , and ( ) F (when the latter can be considered) possess identical weighted bipartite graphs. As a consequence, applying Algorithm 6 (the Σ-method) to F , (•) F , or ( ) F yields the same result in terms of offsets and outcome of the DM decomposition.

Recall the following definition: An execution tree for Algorithm 10 is a labeled tree T representing all its possible runs. Nodes of T are labeled by one of the 10 action labels of Algorithm 10 and each maximal path of T describes a possible run of Algorithm 10 through the sequence of action labels sitting at the nodes of this path. In this section, we investigate the effect of selecting one of the two alternatives in (101), each time index reduction is performed.

Performing (102) results in corresponding variations of the context ∆ at the next instant. We thus consider the set T of all executions trees obtained by performing (102) in Algorithm 10. From Lemma 52 we immediately deduce:

Lemma 53 T is a singleton consisting of an execution tree that we denote by T S . In addition, the contribution, to the bipartite graph, by the context ∆, is independent from these choices. The same holds for Algorithm 13.

Lemma 53 means that Algorithms 10 and 13 will visit the same sequence of action labels, irrespectively of all the choices made (either shifting or differentiating). The systems of equations produced by Line 3 of Algorithm 10 and Line 4 of Algorithm 13 may differ, but possess identical weighted bipartite graphs. From now on we focus on Algorithm 13, of which Algorithm 10 is a specialization.

Inria

Executing a nonstandard instant, revisited

In Algorithm 13, forward shifting is used when performing index reduction. In Section 8, however, we have seen that, in Algorithm 13 forward shifting and differentiation can be substituted for one another. This leads us to consider the following strategy, whose motivation is to make standardization easier (see the standardisation of the Cup-and-Ball example in Section 3.4):

Strategy 54 In Algorithm 13, we priorize shifts or differentiations depending on the following cases:

1. The considered nonstandard instant belongs to a cascade of events: we priorize shifts;

2. The considered nonstandard instant is in a continuous mode: we priorize differentiation;

3. We are entering a continuous mode: we priorize differentiation unless it brings leading derivatives.

The justification of point 2 is clear: we should reproduce the usual index reduction. In a cascade of mode changes occurring at successive nonstandard instants, we are in a discrete time dynamics and thus point 1 follows. Finally, point 3 is justified by the fact that, if x n is the leading derivative for x in the new mode, we need to know restart values for the lower derivatives x, x , . . . , x (n-1) , given by the values for the variables x • , x • , . . . , x (n-1)• .

Generic Blocks

Despite all the different choices (either or •) give rise to the same execution tree T S by Lemma 53, they may result in different standardizations. We thus want to refine our understanding of the effect of choices (102) when applying Algorithm 13. What we are really interested in is what happens at Line 10 where enabled blocks are solved: it may matter whether shift or differentiation was used at each particular traverse of Line 4 of Algorithm 13. To this end we investigate the effect, on a given e ∈ E S , of the choices (102) performed in Algorithm 13, where E S is the set of all the guarded equations of a given ∂AE system S. The successive application of Line 4 and Line 7 of Algorithm 13 yields blocks of a specific form that we investigate next.

Notations: Let w ∈ { , •} * be a word of length n, let w[0], . . . , w[n -1] be the successive prefixes of w, starting from the empty prefix. For f : R K → R a smooth numerical function and w as above, set

B(f, w) = def B h (f, w) B t (f, w) where B h (f, w) = def      f (w[0]) f (w[1]) . . . f (w[n-1])      and B t (f, w) = def f w B(f, w) is called an (f, w)-block, B h (f, w
) is its head and B t (f, w) is its tail. For {f 1 , . . . , f k , . . . , f K } a tuple of K smooth functions in K dependent variables, w 1 , . . . , w k , . . . , w K ∈ { , •} * a set of K words of respective lengths n 1 , . . . , n k , . . . , n K , define

F =    f 1 . . . f K   , W =    w 1 . . . w K   , B(F, W )=    B(f 1 , w 1 ) . . . B(f K , w K )    (103)
and the definitions of B h (F, W ) and B t (F, W ) follow accordingly. In addition, we will write

B h • (f, n) = def B h (f, •n) and B h (f, n) = def B h (f, n)
Exchanging ↔ •: In this paragraph we detail such exchanges for the above identified blocks. Let f : R m → R be a function and let f (X)=0 be the associated equation with m dependent variables collected in the vector X. Consider w∈{ , •} * , and let

f w → f ( n) , f (•n) , where n = |w| , (104)
be the map defined by the following algorithm:

1. Initialize h 0 = g 0 := f w and v 0 := w;

2. For j = 1, . . . , n:

(a) Decompose v j-1 := v j .ν where ν ∈ { , •};

(b) case ν = : define h j := h j-1 and

g j := g j-1 x ← x • -x ∂ , (105)
where x ranges over the set of derivatives involved in g j-1 ; case ν = • : h j-1 has the form • (Y ) for some function h and Y its set of dependent variables; define g j := g j-1 and

h j (Y ) := ∇ (Y ).Y , ( 106
)
where ∇ is the Jacobian of .

Return f

( n) , f (•n) = (h n , g n ).
Lemma 55 Applying map (104) to all the entries of B h (f, w), B h (f, w), and B t (f, w), returns

B h (f, n) and B t (f, n) , B h • (f, n) and B t • (f, n) , (107)
respectively, where n = |w|. Doing the same for B(F, W ) defined in (103) returns

B (F, N ) and B • (F, N ) ( 108
)
where N is the vector of integers collecting the lengths of the entries of W .

Comment 56 By regarding all x , x , etc. as dummy variables,21 every structurally nonsingular block B(F, W ) = 0 defines a dAE system over nonstandard time basis T. This also holds for B (F, N ) = 0 and B • (F, N ) = 0. On the other hand, by not involving shifts, B (F, N ) = 0 defines a (continuous time, standard) DAE system. This is not true for B(F, W ) = 0 in general.

This finishes our study of structural analysis for mDAE or mdAE systems. Referring to Figure 18 in Section 5, we are ready to move to the two downgoing arrows, consisting in the standardization of the nonstandard model resulting from the structural analysis, together with the impulse analysis. The rest of the paper is devoted to this.

Background on nonstandard analysis

The development of the clutch example in Section 2 consisted of three main steps: (1) the mapping of the original clutch model to its nonstandard semantics, (2) the symbolic execution of this nonstandard semantics (Section 2.4), and (3) the back-standardization of this symbolic execution (Section 2.5). All of this resulted in the actual code for the simulation results reported in Figure 9.

Formally justifying back-standardization requires a non-trivial use of nonstandard analysis, way beyond the quick and superficial introduction we gave in Section 2.3. In this and the next sections, we provide the mathematical background supporting back-standardization procedures in full generality. Nonstandard analysis was proposed by Abraham Robinson in the 1960s to allow for the explicit manipulation of "infinitesimals" in analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Cutland | Nonstandard analysis and its applications[END_REF]. Robinson's approach was axiomatic, by adding three new axioms to the basic Zermelo-Fraenkel (ZFC) framework. An alternative presentation was later proposed by Lindstrøm [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF]. Its interest is that it does not require any fancy axiomatic material but only makes use of Zorn's lemma, equivalent to the axiom of choice in the ZF set theory. The proposed construction bears some resemblance to the construction of R as the set of equivalence classes of Cauchy sequences in Q modulo the equivalence relation

∼ defined by (u n ) ∼ (v n ) iff lim n→∞ (u n -v n ) = 0.
The important point for us is that nonstandard analysis allows the use of the nonstandard discretization of continuous dynamics "as if" it was operational and with infinitesimal discretization error. Iwasaki et al. [START_REF] Iwasaki | Modeling time in hybrid systems: How fast is "instantaneous[END_REF] first proposed using nonstandard analysis to discuss the nature of time in hybrid systems. Bliudze and Krob [START_REF] Bliudze | Modelling of complex systems: Systems as dataflow machines[END_REF][START_REF] Bliudze | Un cadre formel pour l'étude des systèmes industriels complexes: un exemple basé sur l'infrastructure de l'UMTS[END_REF] have also used nonstandard analysis as a mathematical support for defining a system theory for hybrid systems. In their study of mathematical foundations of hybrid systems modelers of the ODE class, Benveniste et al. [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF] used extensively the nonstandard semantics of hybrid systems, with practical consequences for how to design a specification formalism for multimode ODE systems [START_REF] Bourke | Zélus: A synchronous language with ODEs[END_REF][START_REF] Benveniste | A Type-based Analysis of Causality Loops in Hybrid Systems Modelers[END_REF]. 22 See also [START_REF] Suenaga | Hyperstream processing systems: nonstandard modeling of continuous-time signals[END_REF][START_REF] Nakamura | A nonstandard functional programming language[END_REF][START_REF] Hasuo | Metamathematics for systems design -comprehensive transfer of formal methods techniques to cyber-physical systems[END_REF] for similar studies using the axiomatic approach to nonstandard analysis. The following presentation is borrowed from [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF].

Intuitive introduction

We begin with an intuitive introduction to the construction of the nonstandard reals. The goal is to augment R ∪ {±∞} by adding, to each x in the set, a set of elements that are "infinitesimally close" to it. We will call the resulting set R. Another requirement is that all operations and relations defined on R should extend to R.

A first idea is to represent such additional numbers as convergent real sequences. For example, elements infinitesimally close to the real number zero are the sequences u n = 1/n, v n = 1/ √ n and w n = 1/n 2 . Observe that the above three sequences can be ordered: v n > u n > w n > 0 where 0 denotes the constant zero sequence. Of course, infinitely large elements (close to +∞) can also be considered, e.g., sequences x u = n, y n = √ n, and z n = n 2 . Unfortunately, this way of defining R does not yield a total order since two sequences converging to zero cannot always be compared: if u n and u n are two such sequences, the three sets {n | u n > u n }, {n | u n = u n }, and {n | u n < u n } may even all be infinite. The beautiful idea of Lindstrøm is to enforce that exactly one of the above sets is important and the other two can be ignored.

The key step in Lindstrøm's construction consists in fixing once and for all a finitely additive positive measure µ over the set N of integers with the following properties:

23 µ : 2 N → {0, 1} µ(X) = 0 whenever X is finite µ(N) = 1
Once µ is fixed, one can compare any two sequences u n and u n as follows. Exactly one of the three sets {n | u n > u n }, {n | u n = u n }, or {n | u n < u n } has µ-measure 1 whereas the other two must

Nonstandard reals and integers

The above general construction applies to X = R and I = N. The result, denoted by R, is a field (according to the transfer principle). By the same principle, R is totally ordered by Proof To prove this, let

[u n ] ≤ [v n ] iff µ{n | u n >v n } = 0.
x = sup{u ∈ R | [u] ≤ [x n ]}. Since [x n ] is finite,
x exists and we only need to show that [x n ]-x is infinitesimal. If not, then there exists y ∈ R, y>0 such that either

[x] < [x n ] -[y] or [x] > [x n ] + [y]
, which both contradict the definition of x. The uniqueness of x is clear, thus we can define st([

x n ]) = x.
It is also of interest to apply the general construction (111) to X = I = N, which results in the set N of nonstandard natural numbers; Z is defined similarly, from X = Z and I = N. The nonstandard set N differs from N by the addition of infinite natural numbers, which are equivalence classes of sequences of integers whose essential limit is +∞.

Internal functions and sets

Any sequence (f n ) of functions f n : R → R pointwise defines a function [f n ] : R → R by setting

[f n ]([x n ]) = def [f n (x n )] . (112)
To justify this definition based on a representant of the equivalence class, note that [x n ] = [x n ] means that the set of integers n such that x n = x n is neglectible, and we have f n (x n ) = f n (x n ) for every n outside of this set. A function R → R obtained in this way is called internal. Properties of, and operations on, ordinary functions extend pointwise to internal functions of R → R. The same notions apply to sets by considering their characteristic functions.

An internal set A = [A n ] is called hyperfinite if µ{n | A n finite} = 1; the cardinal |A| of A is defined as [ |A n | ],
where |A n | denotes the cardinal of the finite set A n in the usual sense. For A a bounded set of hyperreals, we define its shadow Sh(A) ⊆ R as follows:

Sh(A) = {st(x) ∈ R | ∃y∈A such that y ≈ x} . (113)
The nonstandard lifting of f : R → R is the internal function f = [f, f, f, . . . ] . Let f : R → R be a standard function that is differentiable at x ∈ R, and f = [f, f, f, . . . ] its nonstandard lifting. Then, we have, for every x ∈ R :

f (x) = st f (y) -f (x) y -x (114)
for every y ∈ R such that y ≈ [x]-the proof is similar to that of Lemma 59. Now, consider an infinite number K ∈ N and the set

T = 0, 1 K , 2 K , 3 K , . . . K -1 K , 1 . (115) By definition, if K = [K n ], then T = [T n ] with T n = 0, 1 K n , 2 K n , 3 K n , . . . K n -1 K n , 1 hence |T | = [ |T n | ] = [K n + 1] = K + 1.
Note that the shadow of T (see ( 113)) is the standard interval [0, 1]: for any given real number a in [0, 1], aK K standardizes to a. In our nonstandard semantics, we will be using the following time basis:

T = def {k∂ | k∈ N} (116)
where ∂ > 0 is some fixed infinitesimal time. This definition is reminiscent of Equation ( 14), in Section 2.3. The shadow of T is R + , the set of nonnegative real numbers, and that, informally speaking, T is both "discrete" (every element has a previous and a next element, except for 0) and "continuous" (as the shadow of T is R + ).

Notations 60 Generic elements of T will be denoted by the special symbol t or simply t when no confusion can result.

Integrals

Now, consider an internal function f = [f n ] and a hyperfinite set

A = [A n ]. The sum of f over A can be defined by a∈A f (a) = def a∈An f n (a) .
If T is as above, and f : R → R is a standard function, we obtain

t∈T 1 |T | f (t) = t∈Tn 1 |T n | f (t) . (117)
Now, the continuity of f implies the convergence of the Riemann sums: t∈Tn

1 |Tn| f (t) → 1 0 f (t)dt. Hence, 1 0 f (t)dt = st t∈T 1 |T | f (t) . (118)
Under the same assumptions, for any t ∈ [0, 1],

t 0 f (u)du = st   u∈T,u≤t 1 |T | f (u)   . (119)

ODE

Consider the following ODE:

x = f (x, t), x(0) = x 0 . (120)
Assume (120) possesses a solution x : [0, 1] → R such that the function t → f (x(t), t) is continuous. Rewriting (120) in its equivalent integral form x(t) = x 0 + t 0 f (x(u), u)du and using (119) yields

x(t) = st( x(t)) where x(t) = def x 0 + u∈T,u≤t 1 |T | f (x(u), u) . (121)

Inria

One can rewrite the positive infinitesimal quantity 1/|T | as ∂, so that

T = {t k = k∂ | k = 0, . . . , |T |} .
Then, after substituting in (121), one gets that the piecewise-constant right-continuous function x(t), t∈ R, 0≤t≤1 satisfies the following difference equation, for k=0, . . . , |T |-1,

x(t k+1 ) = x(t k ) + ∂ . f ( x(t k ), t k ) x(t 0 ) = x 0 (122)
By (121), the following holds:

Theorem 61 The solution x(t), t ∈ R, of ODE (120), and x(t), t ∈ R + , the piecewise constant interpolation of the trajectory of system (122), are related by x = st( x).

The same line of arguments applies for systems of ODEs possessing a unique continuous solution, i.e., x and f take their values in R m in (120). Formula (122) can be seen as a nonstandard semantics for ODE (120). This semantics seems to depend on the choice of infinitesimal step parameter ∂. Property (121), though, expresses that all of these nonstandard semantics are equivalent from the standard viewpoint, regardless of the choice made for ∂.

The paper [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF] goes beyond our short exposure by including a direct proof of the Peano theorem on the existence of solutions of (120) when f is continuous and bounded.

Infinitesimal calculus

We conclude this background material with a summary of the section "Infinitesimal calculus" from [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF]. We formulate the essential characterizations of properties of the infinitesimal calculus in nonstandard terms. f denotes a standard function and f = [f, f, f, . . . ] is its nonstandard lifting. The reader is referred to that text for the proofs.

Theorem 62

1. The function f : R → R is continuous at a ∈ R if and only if f (x) ≈ f (a) for all x ≈ a.
2. The function f : R → R is uniformly continuous on set A ⊆ R if and only if f (x) ≈ f (y) for all x, y ∈ A such that x ≈ y.

The function

f : R → R is differentiable at a ∈ R if and only if there exists a number b ∈ R such that f (x) -f (a) x -a ≈ b , for all x ∈ R, x ≈ a, x = a.
Moreover, if such a number b exists, then it equals f (a).

The toolkit supporting standardization

By building on the background of Section 9, we develop here the specific material we need for our developments. First, we develop the compile-time impulse analysis. Second, we formally justify our use of structural analysis in the nonstandard domain. Third, referring to Comment 4 in support of the clutch example, we formally justify our method for standardizing systems of equations. We consider systems of equations of the form

0 = H(X , X • , V, X) , (124)
To identify impulsive and non-impulsive variables in a systematic way, we develop the impulse analysis, which consists in abstracting hyperreals with their "magnitude order" compared to the infinitesimal ∂. Since our original DAE models are all standard, only the occurrence of the infinitesimal ∂ will have a role in this matter. The impulse analysis will be useful as a preparation step for standardization.

Definition 64 (impulse order and analysis)

1. Given a system of equations such as (124) or (125), say that a dependent variable x has impulse order (or simply order) o ∈ R, if the solution of the considered system is such that x∂ o is provably a finite non-zero (standard) real number. We denote by [[x]] the impulse order of x. By convention, the constant 0 has impulse order -∞.

Say that

x is impulsive if [[x]] > 0.
3. The impulse analysis of a system of equations such as (124) or (125) is the system of constraints satisfied by the impulse orders of the dependent variables of the system. An expression e is 0, a nonzero (standard) real constant c, the infinitesimal ∂, a variable x, the monomial e c , a sum, or a product. Figure 21 describes the syntax of a mini-language specifying such systems of equations. The left column of Figure 22 gives the rules for mapping expressions to their corresponding impulse orders. All the rules are self-explanatory, with the exception of Rule (R6), where the inequality deserves some explanation. The sum e 1 +e 2 , the dominant terms in the expansion of the e i 's as power series over ∂ may compensate each other for their impulse orders. For an example of this, see equation (e ∂ 1 ) in System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]: rewriting this equation as a 1 (ω 1 ) =

(R1) [[0]] = -∞ (R2) [[c]] = 0 (R3) [[∂]] = -1 (R4) [[e c ]] = c[[e]] (R5) [[e 1 × e 2 ]] = [[e 1 ]] + [[e 2 ]] (R6) [[e 1 + e 2 ]] ≤ max{[[e 1 ]], [[e 2 ]]} E e = e [[E]] [[e]] = [[e ]] (
ω • 1 -ω 1 ∂
-b 1 (ω 1 )τ 1 , we see a case of strict inequality for (R6) since a 1 (ω 1 ) has order zero, whereas it is equal to the difference of two terms of order one, see Section 10.1.4 for details.

We will use Rule (R6) in the following way, thereby reinforcing it. Consider an equation e : z = x + y .

We can rewrite e in the following equivalent ways: 0 = x + y -z , x = z -y, or y = z -x. To each of them we apply the max rule. This yields, for the impulse analysis of equation e, the following system of constraints:

impulse analysis of e :

[[z]] ≤ max{[[x]], [[y]]} ; [[0]] ≤ max{[[x]], [[y]], [[z]]} [[x]] ≤ max{[[z]], [[y]]} ; [[y]] ≤ max{[[x]], [[z]]} (126)
Note that the constraint [[0]] ≤ . . . is vacuously satisfied since [[0]] = -∞. Also note that, among the three nontrivial inequalities of (126), at least two must be saturated. We will use impulse analysis (126) for handling sums of terms. This reinforcement of the max rule is formalized by Rule (R8) of Figure 22, which mechanizes the association, to equation e of (126), of its different rewritings.

Using the rules of Figures 21 and22 in the numerical expressions, we map any system of rational equations of the form (124) or (125) into a system of constraints over impulse orders.

To cover functions beyond polynomials, we need to extend R∪{-∞} with +∞. In this extension, we take the convention that -∞ + ∞ = -∞, justified by the equality 0×x = 0 for any nonstandard x. For functions f (x) = ∞ k=0 a k x k that can be represented as absolutely converging power series, we then get

[[f (x)]] = [[ ∞ k=0 a k x k ]] = [[x]]. sup(A), where A = {k | a k = 0} (127)
is the support of the series and sup(A) is the supremum of set A. In particular, if [[x]] > 0 and if the support of the series is infinite, we get [[f (x)]] = +∞.

Impulse analysis of systems of equations for restarts

Here we particularize the impulse analysis to systems of equations of the form (124), where the only reason for ∂ to occur is the expansion of derivatives using the Euler scheme:

0 = H X • -X ∂ , X • , V, X
The dependent variables are X • , V . It will be convenient to introduce the auxiliary variables

U = def X • -X ,
so that the systems we consider take the following form:

dependent variables X • , V, U in : 0 = H U ∂ , X • , V, X U = X • -X (128)
The following condition for System (128) can be assumed, based on physical considerations (restart values for an ODE or a DAE cannot be impulsive in a physically meaningful model):

Assumption 7 Since X is a state, both X (a known value) and X • must be finite.

First, the impulse orders [[X]] are all known, from previous nonstandard instants. Next, from Assumption 7 we deduce the inequalities

[[X • ]] ≤ 0 and [[U ]] ≤ 0 . ( 129
)
The impulse orders [[V ]] are a priori unknown. We have, however, more prior information, thanks to the structural analysis. As part of both the simpler algorithm ExecRun (Algorithm 10) and revisited algorithm ExecRun (Algorithm 13), the auxiliary algorithm SolveConflict (Algorithm 9) is called. At the considered mode change, we thus know which consistency equation(s) of the new mode was/were conflicting with the dynamics of the previous mode. Formally, call G=0 the subsystem collecting Inria all the equations that were erased by algorithm SolveConflict at Line 6 of Algorithm 10 or Line 7 of Algorithm 13. As a result, G=0 is no longer satisfied at the considered mode change, and thus, G defines a tuple R of finite nonzero variables called residuals, by setting

R = G . ( 130
)
An example of residual in the clutch is r = ω 1 -ω 2 , which is both finite and nonzero at mode change γ : f → t; this information was found by the structural analysis. Finally, the system of equations that we need to solve collects all the above items, namely:

dependent variables X • , V, U, R in :      0 = H U ∂ , X • , V, X U = X • -X R = G (131)
prior information on impulse orders :

         [[ 1 ∂ ]] = 1 [[X • ]] ≤ 0 [[U ]] ≤ 0 [[R]] = 0 (132)
Decomposing algebraic dependent variables into impulsive and non-impulsive ones In the following, we assume that the vector functions H or H can be represented as absolutely convergent power series in their arguments. We can then apply the rules of Figures 21 and22, complemented with (127), to derive the system of equations that impulse orders must satisfy. In particular, this allows us to partition the set V of algebraic variables as

V = W i ∪ W (133)
where W i collects the impulsive variables, having impulse order > 0, and W collects the nonimpulsive variables. In Section 11.3, we will propose a numerical scheme for computing restarts that only requires the splitting of algebraic dependent variables into impulsive and non-impulsive ones.

For now, let us develop the impulse analysis for the two transitions γ : t → f and γ : f → t of System ( 16) assuming linear f i as in Eq. [START_REF] Dieudonné | Fondements de l'analyse moderne[END_REF]. We develop this example here based on the rules of Figures 21 and22. To better illustrate this mechanical reasoning, we forbid ourselves algebraic manipulations that otherwise could help simplifying a manual analysis. We repeat here the corresponding system of equations with its expansion (128)-there is no residual for this mode change:

               u 1 ∂ = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 (e ∂ 1 ) u 2 ∂ = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 (e ∂ 2 ) τ 1 = 0 (e 5 ) τ 2 = 0 (e 6 ) u 1 = ω • 1 -ω 1 u 2 = ω • 2 -ω 2 (134)
We then saturate System (134) using Rule (R8) of Figure 22 (added equations are in blue):

               u 1 ∂ = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ; a 1 (ω 1 ) = u 1 ∂ -b 1 (ω 1 )τ 1 ; b 1 (ω 1 )τ 1 = u 1 ∂ -a 1 (ω 1 ) u 2 ∂ = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 ; a 2 (ω 2 ) = u 2 ∂ -b 2 (ω 2 )τ 2 ; b 2 (ω 2 )τ 2 = u 2 ∂ -a 2 (ω 2 ) τ 1 = 0 τ 2 = 0 u 1 = ω • 1 -ω 1 ; u 1 + ω 1 = ω • 1 ; u 1 -ω • 1 = -ω 1 u 2 = ω • 2 -ω 2 ; u 2 + ω 2 = ω • 2 ; u 2 -ω • 2 = -ω 2 (135) RR n°9334
and we consider the prior information (132):

[[ω • i ]] ≤ 0 ; [[u i ]] ≤ 0 ; [[ 1 ∂ ]] = 1 (136)
We now apply the rules of Figure 22 to System (135), then use (136): We repeat here the corresponding system of equations with its expansion (128), and we include the residual r = def ω 1 -ω 2 (from structural analysis, we know it is nonzero), highlighted in red:

               1 + [[u 1 ]] ≤ 0 ; 0 ≤ 1 + [[u 1 ]] ; vacuous 1 + [[u 2 ]] ≤ 0 ; 0 ≤ 1 + [[u 2 ]] ; vacuous [[τ 1 ]] = -∞ [[τ 2 ]] = -∞ [[u 1 ]]
                   u 1 ∂ = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 (e ∂ 1 ) u 2 ∂ = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 (e ∂ 2 ) ω • 1 -ω • 2 = 0 (e • 3 ) τ 1 + τ 2 = 0 (e 4 ) u 1 = ω • 1 -ω 1 u 2 = ω • 2 -ω 2 r = ω 1 -ω 2 (138)
We then saturate System (138) using Rule (R8) of Figure 22 (added equations are in blue):

                   u 1 ∂ = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ; a 1 (ω 1 ) = u 1 ∂ -b 1 (ω 1 )τ 1 ; b 1 (ω 1 )τ 1 = u 1 ∂ -a 1 (ω 1 ) u 2 ∂ = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 ; a 2 (ω 2 ) = u 2 ∂ -b 2 (ω 2 )τ 2 ; b 2 (ω 2 )τ 2 = u 2 ∂ -a 2 (ω 2 ) ω • 1 -ω • 2 = 0 ; ω • 1 = ω • 2 τ 1 + τ 2 = 0 ; τ 1 = -τ 2 u 1 = ω • 1 -ω 1 ; u 1 + ω 1 = ω • 1 ; u 1 -ω • 1 = -ω 1 u 2 = ω • 2 -ω 2 ; u 2 + ω 2 = ω • 2 ; u 2 -ω • 2 = -ω 2 r = ω 1 -ω 2 (139)
We also consider the prior information (132):

[[ω • i ]] ≤ 0 ; [[u i ]] ≤ 0 ; [[r]] = 0 ; [[ 1 ∂ ]] = 1 . ( 140
)
We now apply the rules of Figure 22 to System (139), then use (140):

                       1 ≤ [[τ 1 ]] ; vacuous ; [[τ 1 ]] ≤ 1 1 ≤ [[τ 2 ]] ; vacuous ; [[τ 2 ]] ≤ 1 vacuous ; [[ω • 1 ]] = [[ω • 2 ]] vacuous ; [[τ 1 ]] = [[τ 2 ]] vacuous ; vacuous ; 0 ≤ max{[[u 1 ]], [[ω • 1 ]]} vacuous ; vacuous ; 0 ≤ max{[[u 2 ]], [[ω • 2 ]]} 0 = [[r]] = [[ω 1 ]] = [[ω 2 ]] [[ω • i ]] ≤ 0 , [[u i ]] ≤ 0 , i = 1, 2 (141)

Using impulse analysis in code generation

Code generation for restarts consists in standardizing System (124) or System (125). Recall that standardizing systems of equations requires more care than standardizing numbers, due to impulsive behaviors and singularity issues that result. We can exploit the impulse analysis through three different approaches, which we illustrate by using the clutch example and particularly the mode change γ : f → t when the clutch gets engaged. The system of equations for the corresponding restart is [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF], which we recall for convenience:

       ω • 1 = ω 1 + ∂.(a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ) (e ∂ 1 ) ω • 2 = ω 2 + ∂.(a 2 (ω 2 ) + b 2 (ω 2 )τ 2 ) (e ∂ 2 ) ω • 1 -ω • 2 = 0 (e • 3 ) τ 1 + τ 2 = 0 (e 4 )
As we have shown, the two torques τ 1 and τ 2 are impulsive of order 1.

Eliminating impulsive variables When this is practical, the simplest method is to eliminate impulsive variables from the restart system, namely the two torques in System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]. The details were given in Section 2.5.2 and the resulting reduced system of equations is

ω • 1 = ω • 2 = b 2 (ω 2 )ω 1 + b 1 (ω 1 )ω 2 b 1 (ω 1 ) + b 2 (ω 2 ) + ∂ a 1 (ω 1 )b 2 (ω 2 ) + a 2 (ω 2 )b 1 (ω 1 ) b 1 (ω 1 ) + b 2 (ω 2 ) , ( 142
)
whose standardization is simply achieved by substituting ∂ ← 0. No numerical information is provided on the torques with this method. This is a satisfactory solution when elimination of impulsive variables is practical. In our example, they entered linearly in the restart system, so that elimination was straightforward. When this is not the case, elimination becomes costly or even impossible. We thus need to look for alternatives.

Rescaling impulsive variables Since [[τ i ]] = 1, we apply the change of variable τ i = ∂ τi , which transforms System (138) into:

       ω • 1 -ω 1 = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 (e ∂ 1 ) ω • 2 -ω 2 = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 (e ∂ 2 ) ω • 1 -ω • 2 = 0 (e • 3 ) τ1 + τ2 = 0 (e 4 ) (143)
which yields again the solution (142) for ω • 1 = ω • 2 . However, the rescaled value τ1 can now be computed. Rescaling impulsive variables is simpler than eliminating them. Unfortunately, this method does not work in full generality since impulse orders can be infinite as we have shown in (127). The last method addresses such cases, at the price of a possibly poor numerical conditioning.

Bruteforce solving of the restart system Consider once more System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF], replace in it the infinitesimal time step ∂ by a small real positive time step δ > 0, and expand the derivatives as before. This yields

       ω • 1 = ω 1 + δ.(a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ) (e δ 1 ) ω • 2 = ω 2 + δ.(a 2 (ω 2 ) + b 2 (ω 2 )τ 2 ) (e δ 2 ) ω • 1 -ω • 2 = 0 (e • 3 ) τ 1 + τ 2 = 0 (e 4 ) (144)
Then, it will be proved in Theorem 91 of Section 11.3 that solving System (144) for its dependent variables and then discarding the values found for the impulsive variables yields a converging approximation for the states ω + 1 and ω + 2 at restart. Of course, numerical conditioning for System (144) is likely to be less good than for System (143), so that rescaling is recommended when impulse orders are finite.

Qualitative information obtained from impulse analysis The dynamical profile of impulsive variables yields interesting information about the nature of the impulsive behaviors. For the clutch example, we get the following profile for [[τ i ]] : . . . , 0, 1, 0, . . . , where the 1 occurs at the event γ : f → t. This indicates a Dirac-like impulse. As a second example (not developed here), the impulse analysis of the Cup-and-Ball example with inelastic impact yields the following profile for the impulse order of the tension [[λ]] : . . . , 0, 1, 0, . . . , where the 1 occurs at the event following the mode change γ : f → t, i.e., when System (34) gets executed to evaluate restart values for velocities. This again indicates a Dirac impulse for the tension. Note that having a profile . . . , 0, k, 0, . . . for an impulsive variable, where k is a positive integer, still indicates a Dirac behavior. The amplitude does not matter: only the profile matters. Developing a systematic analysis of this kind requires further study and is left for future investigations.

Nonstandard Structural Analysis

The motivation for this section was developed in Section 2.5. The following two questions will be addressed:

1. We apply structural analysis to the nonstandard semantics. Is this justified? In the standard setting, the structural analysis was motivated by the quest for low cost criteria for almost everywhere nonsingularity of systems of algebraic equations. Can we have a similar argument here?

2. We need to formalize and generalize the Comment 4 of Section 2.5.2, where we discussed the difference between standardizing functions and standardizing equations.

We begin with question 1. Algorithm 13 returns at its Line 10 systems of algebraic equations of the form H(Z, X, ∂) = 0 introduced in (125), where tuple Z collects the dependent variables and tuple X collects the state variables. We justify the call of the atomic action Eval(σ, H) by the fact that "H is structurally nonsingular". But the justification of structural analysis was developed in Section 6.1.1, i.e., in a standard setting. We thus need a justification for taking the liberty of using it for nonstandard semantics as well. Formally, we consider nonstandard systems of equations of the form

H (Z, X, ∂) =0,
µ(A) = 1, where A = {n | H (z n , x n , δ n ) =0} ( 146
)
and µ is the finitely additive measure over N introduced in (109). Now, assume that H is structurally nonsingular according to Definition 65. Then, using notations (55) of Section 6.1.1, for every n ∈ A, there exist U n = def U Hn (z n , x n , δ n ) and V n = def V Hn (z n , x n , δ n ) satisfying property [START_REF] Zimmer | Equation-Based Modeling of Variable-Structure Systems[END_REF] with respect to the system H (Z n , X n , δ n ) =0. Define the internal sets U = [U n ] and V = [V n ]. We have λ(V \ U) = 0, where λ is the internal set function λ = [λ, λ, λ, . . . ]. For every n ∈ A, every J n ∈ V n yields an invertible matrix, implying that J = def [J n ] also yields an invertible matrix. As a result, we have exhibited two sets U and V Inria satisfying property [START_REF] Zimmer | Equation-Based Modeling of Variable-Structure Systems[END_REF] with respect to the system H (Z, X, ∂) = 0. Definition 65 thus formally justifies our shameless use of the structural analysis in nonstandard domains. Note that Definition 65 tells nothing about how U is-it could be infinitely small. Neither does it require uniqueness of the solutions of the considered equations.

We now move to question 2 and investigate how equations should be standardized.

Theorem 67 (standardizing equations) Consider a square system of the form (145), and assume that H is of class C 1 . The following two properties hold:

1. If the standard pair (z, x) satisfies H(z, x, 0)=0 and the Jacobian J Z H at the point (z, x, 0) is nonsingular, then, for every nonstandard pair ( x, ∂) such that ( x, ∂) ≈ (x, 0), there exists z such that st( z)=z and H ( z, x, ∂) =0.

2. Vice-versa, let ( z, x, ∂) be a finite nonstandard triple satisfying H ( z, x, ∂) = 0, and let (z, x, 0) = st( z, x, ∂) be its standard part. Then H(z, x, 0) = 0 holds.

Under these assumptions, the standardization of the equation H ( z, x, ∂) =0 is independent from the particular value for ∂ we have chosen (provided that it is infinitesimal).

Corollary 68 We can standardize H(Z, X, ∂)=0 as the system H(Z, X, 0)=0, provided that the latter is structurally nonsingular.

The nonsingularity condition is needed for statement 1 to hold. We proceed to the proof of Theorem 67.

Proof We successively prove the two statements. For Statement 1, by the Implicit Function Theorem applied to the standard function H, there exists a C 1 -class function (x , δ) → G(x , δ) such that H(G(x , δ), x , δ)=0 holds for (x , δ) ranging over a neighborhood V of (x, 0). Let For Statement 2, let the finite triple ( z, x, ∂) satisfy H ( z, x, ∂) = 0 and let (z, x, 0)=st( z, x, ∂). Since (z, x, 0)≈( z, x, ∂) and H is of class C 1 , we can apply Statement 2 of Theorem 62 to derive H(z, x, 0) ≈ H ( z, x, ∂) = 0, whence H(z, x, 0) = 0 follows since H(z, x, 0) is standard.

DAE of index zero

Invoking the implicit function theorem, one can also address DAEs of index 0, i.e., of the form:

F (X , X) = 0 ( 148
)
where the Jacobian ∇ V F (V, X) is invertible in a neighborhood of the solution of the equation F (V, X) = 0, where V is the vector of dependent variables.

Theorem 69 Let F (X , X) = 0 be a structurally nonsingular DAE system of index 0, and let F = [F, F, F, . . . ] be the nonstandard lifting of F . Consider the following nonstandard transition system having X • as dependent variables:

0 = F X • -X ∂ , X , X(0) = X 0 (149)
Then, system (149) possesses a unique solution and we denote by X(t), t∈ R + its piecewise constant interpolation. Furthermore, the solution X of DAE system F (X , X)=0, X(0)=X 0 satisfies X = st( X).

Proof The theorem is proved by componentwise induction on N , with the help of the Lemma 71 to follow.

Lemma 71 Consider a structurally nonsingular nonstandard dAE system of the form

H(X • , Y • , X, Y ) = def   G(X • , Y • , X, Y ) g(X) g(X • )   = 0 (154)
with leading variables X • , Y • , and consider

H(X • , Y • , X, Y ) = def   G(X • , Y • , X, Y ) g(X) (Jg)(X).X   = 0 (155)
where Jg(X) is the Jacobian of the function g at X. Let ( X(t), Y (t)), t ∈ R + , and ( X(t), Y (t)), t ∈ R + , be the piecewise constant interpolations of the solutions of nonstandard dAE systems (154) and (155), respectively. Assume that we know a standard DAE block

G(X , Y , X, Y ) = 0 such that the solution (X(t), Y (t)) t∈R of H(X , Y , X, Y ) = def   G(X , Y , X, Y ) g(X) Jg(X).X   = 0 (156)
exists, is unique, and satisfies (X, Y )=st X, Y . Then, we also have (X, Y )=st( X, Y ) .

Proof By statement 3 of Theorem 62, we have

g(X • ) = g(X) + ∂. (Jg)(X). X + ε where ε is infinitesimal. Hence H rewrites H(X • , Y • , X, Y ) =   G(X • , Y • , X, Y ) g(X) g(X) + ∂. (Jg)(X). (X + ε)   Since H(X • , Y • , X, Y ) = 0 implies in particular g(X) = 0, (154) rewrites   G(X • , Y • , X, Y ) g(X) (Jg)(X). (X + ε)   = 0
and removing the infinitesimal ε does not change the standardization of this dAE system. We conclude by invoking the assumption of the lemma.

Mode changes of mDAE

Assumption 6 (stated on page 73) is in force throughout this section. We study the standardization of mode change events. Having successfully performed the nonstandard structural analysis, the system at mode changes can be put in the generic form (124), and we use the decomposition (133) of the algebraic variables into impulsive and non-impulsive ones. For convenience, we repeat these formulas here:

(124) : X = X • -X ∂ in 0 = H(X , X • , V, X) (133) : V = W i ∪ W
Note that, in (124), the infinitesimal parameter ∂ is involved in both time (by being the time step of T) and space (since derivatives x in the original standard model are expanded as x • -x ∂ ). Concrete instances of this are found in the clutch example, e.g., in System [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF].

Compared with the previous section that dealt with long modes, the novelty, here, is that we take a different target domain for the standardization, namely the domain of discrete-time transition systems-standard transition systems indexed by N. Clearly, we cannot expect mapping to this domain the entire trajectory of System (124). We only aim at standardizing the finite prefix of it needed to specify the restart conditions for the new mode. The "•" shift operator defines the discrete time of this dynamics, hence we do not need to eliminate ∂ from the time domain. In contrast, we will have to eliminate the occurrences of the infinitesimal parameter ∂ acting in space. Consider again the clutch example and in particular System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF], used to determine the restart values ω • 1 , ω The standardization of (124) relies on Theorem 67. Unfortunately, it is not true in general that substituting 0 for ∂ in (124) yields a structurally nonsingular system with respect to the dependent variables Z = def (X • , V ). We thus need further work to be able to invoke Theorem 67. This will consist in eliminating, from System (124), the impulsive variables. For the example of the clutch, the corresponding analysis was developed in Section 2.5.2. In the following, we propose a systematic approach.

To this end, we also use the rewriting of (124), i.e., we use both:

(124) : 0 = H(X , X • , V, X)

(125) : 0 = H(Z, X, ∂)
where we recall that the form (125) is obtained by expanding X as X • -X ∂ in (124)-the same notational convention is used in the assumption below. Using decomposition (133) of algebraic variables into impulsive and non-impulsive ones: V = W i ∪ W , and writing

Y = def (X • , W ) ,

we assume the following about (124):

Assumption 8 There exists a reduced system K(X , X • , W, X) = 0 with dependent variables Y , satisfying the following two conditions:

1. System K(Y, X, 0)=0 is structurally nonsingular, where K(Y, X, ∂) is derived from K(X , X • , W, X)
as explained above;

2. The considered model is a reduction of the original one, in that, equivalently:

K(X , X • , W, X) = 0 ⇐⇒ ∃W i .H(X , X • , V, X) = 0 K(Y, X, ∂) = 0 ⇐⇒ ∃W i .H(Z, X, ∂) = 0 Inria mode ¬γ : index 0 τ 1 = 0; τ 2 = 0; ω 1 = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ; ω 2 = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 start mode γ : index 1 τ 1 = (a 2 (ω 2 ) -a 1 (ω 1 ))/(b 1 (ω 1 ) + b 2 (ω 2 )); τ 2 = -τ 1 ; ω 1 = a 1 (ω 1 ) + b 1 (ω 1 )τ 1 ; ω 2 = a 2 (ω 2 ) + b 2 (ω 2 )τ 2 ; constraint ω 1 -ω 2 = 0 when γ do τ 1 = NaN; τ 2 = NaN; ω + 1 = b2(ω - 2 )ω - 1 +b1(ω - 1 )ω - 2 b1(ω - 1 )+b2(ω - 2 )
;

ω + 2 = ω + 1 done when ¬γ do τ 1 = 0; τ 2 = 0; ω + 1 = ω - 1 ; ω + 2 = ω - 2 done
Figure 23: Standardization of ExecRun for the ideal clutch (Figure 20). Blocks have been standardized and then symbolically pivoted. x -is the previous value of state variable x, i.e., the left limit of x when exiting a mode; x + is the value for restart. Long modes are colored blue; non-impulsive (resp. impulsive) state-jumps are colored black (resp. red). The dynamics in mode ¬γ is defined by an ODE system, while in mode γ, it is defined by an over-determined index-1 DAE system consisting of an ODE system coupled to an algebraic constraint. In the transition from mode ¬γ to mode γ, variables τ 1 and τ 2 are impulsive, and their standardization is undefined. This explains why they are set to NaN (Not a Number).

In words, we can eliminate, from the original model, the impulsive variables (belonging to W i ), and the result is structurally nonsingular for ∂ = 0. We constructed such a reduced system in Section 2.5.2 by eliminating the torques for the released-to-engaged mode change in the clutch example.

Theorem 72 Assumptions 6 and 8 are in force. Then, standardizing System (124) is performed by solving, for Y , the structurally nonsingular system K(Y, X, 0) = 0. Note that this standardization is independent from the particular infinitesimal time step ∂.

Note that the latter system is standard.

Proof Theorem 72 is a consequence of Theorem 67.

Theorem 72 allows us to compute the restart values for the state variables that are not impulsive, which is sufficient to restart the dynamics in the coming long mode (impulsive variables at mode change cannot be initial conditions for state variables in the new mode).

We have shown in Section 2.5.2 that Assumption 8 holds for the clutch example. In general, the method would consist in, first, computing the decomposition (133), and, then, eliminating impulsive variables. Figure 23 shows the effective code resulting from the standardization of the labeled automaton of Figure 20.

To get practical, Assumption 8 relies on computer algebraic techniques for eliminating variables from systems of equations. Algebraic elimination is easy for systems in which impulsive variables enter linearly. If they enter polynomially, Gröbner bases can be used, but at a high computational cost. Effective elimination techniques may not exist in general, or may be too costly. Hence, the alternative based on the numerical scheme developed in Section 11.3 is of great practical interest.

Main results

This section collects results that are essential in grounding our approach. First, we prove that our approach is intrinsic, in that the actual code we generate does not depend on the particular nonstandard scheme used, instead of the explicit first-order Euler scheme, to expand the derivative. Then, we compare our approach with known schemes, for subclasses for which they are known to exist. Finally, we propose a numerical scheme to compute restart values at mode changes, in full generality, for programs accepted by our structural analysis. This numerical scheme does not require eliminating impulsive variables using computer algebra.

What if we change the interpretation of the derivative?

So far, an essential step in our reasoning was the interpretation (74) for the derivative ω . However, several nonstandard interpretations of the derivative can be equally considered. What happens if we take an interpretation different from (74)? We first discuss this on the clutch example, then provide a general result.

The clutch example

If the system is within a long continuous mode, the following alternative expansion exists for the derivative of ω: for n ∈ N finite but otherwise arbitrary, ω ≈ (ω • -ω) •n ∂ . More generally, for a=(α P , . . . , α N ) such that P ≤N , N ≥ 0, α N = 0, and N n=P α n = 1, one has:

ω ≈ 1 ∂ N n=P α n (ω • -ω) •n . (157)
Let us investigate how the analysis of the clutch is modified when using the nonstandard interpretation

ω a = def 1 ∂ N n=P α n (ω • -ω) •n (158)
for the derivative. The nonstandard semantics of the clutch model is:

                 ω a1 = f 1 (ω 1 , τ 1 ) (e ∂ 1 ) ω a2 = 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 ) (159)
Algorithm 13 tells us that we must shift forward (e 3 ) by N +1 and keep other equations unchanged; other latent equations (e •n 3 ) for n ≤ N are erased when solving the conflict with the previous mode:

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

Using τ 1 = -τ 2 = def τ , the resetting equations when entering the mode γ = t are

     ω a1 = f 1 (ω 1 , τ ) ω a2 = f 2 (ω 2 , -τ ) ω •(N +1) 1 -ω •(N +1) 2 = 0 (161)
Inria -Let G=0 and • G=0 be the DAE system at the current (t ≥ t o ) and previous (t < t o ) modes, respectively; we denote by (c f , d x ) and (c - g , d - y ) the offsets returned by the Σ-method for the current and previous mode, respectively; -Let (G Σ , G Σ ) and (G a Σ ,(G) a Σ ) be the returns of the Σ-method when applied to G=0 and G a =0;

-Let ( • G Σ , • G Σ ) and ( • G a Σ , • (
G) a Σ ) be the returns of the Σ-method when applied to • G=0 and • G a =0; -We denote by ∆ and ∆ a , and by Φ and Φ a , the contexts and sets of facts constructed by ExecRun if expansions (164) and (165) are used.

Line 6 takes as inputs the systems (G

Σ ∪ G Σ ) \ Φ and (G a Σ ∪(G) a Σ ) \ Φ a .
When applying SolveConflict, the dependent variables for (G Σ ∪ G Σ ) \ Φ are all the variables of this system, both leading and non-leading, whose values were not already set by executing the previous nonstandard instant. The same holds for the dependent variables of (G a Σ ∪(G) a Σ ) \ Φ a . The following result holds:

Theorem 74 Set F = def (G Σ ∪ G Σ ) \ Φ with dependent variables defined by context ∆, and
F a = def (G a Σ ∪(G) a Σ )
\ Φ a with dependent variables defined by context ∆ a . Then, the mappings

ψ f : F a → F defined by ψ f f •M (c f -m) = f •(c f -m) ψ x : X(F a ) → X(F ) defined by ψ x x •M (dx-m) = x •(dx-m) preserve the Dulmage-Mendelsohn decompositions of G F and G F a , i.e., g ∈ β a O/E/U if and only if ψ f (g) ∈ β O/E/U , and y ∈ β a O/E/U if and only if ψ x (y) ∈ β O/E/U .
The proof of this theorem decomposes into a series of lemmas. Let f =0 range over the set of equations of system G=0, and let n = c f be the equation offset computed for f when applying the Σ-method to G=0.

Lemma 75

1. The following formulas hold:

G Σ ∪ G Σ = f ∈G B(f, c f ) and G a Σ ∪ G a Σ = f ∈G B(f, M c f ) (a) ∆ = f ∈ • G B(f, c - f -1) and ∆ a = f ∈ • G B(f, M c - f -1) (b) Φ = f ∈ • G∩G B(f, k f ) and Φ a = f ∈ • G∩G B(f, K f ) (c)
where

k f = min(c - f -1, c f ) and K f = min(M c - f -1, M c f ) (d) (167) 2. If c - f ≤ c f , then k f = c - f -1 and K f = M c - f -1; If c - f > c f , then k f = c f and K f = M c f .
Proof For the first statement, (167-a) is a direct consequence of Lemma 73 and Lemma 42 from Section 7.

> 1, c - f ≤ c f is equivalent to c - f -1 M ≤ c f , or M c - f -1 ≤ M c f . The case c - f ≤ c f follows.
The other case is proved similarly.

In the following lemma, we use Notations 41 and recall that X(G) denotes the set of variables involved in the set G of functions. We also introduce the following notations.

Notations 76 For G=0 a DAE system, we extend by convention its offsets (c f ) f ∈G and (d x ) x∈X(G) to arbitrary f and x by setting

c f = d x = -∞ if f ∈ G and x ∈ X(G).
For -∞ ≤ m and 0 ≤ n, and ζ a variable or a function, we set

B(ζ, m, n) = def      ζ • m ζ •( m+1) . . . ζ •n     
, where m = max(m, 0)

with the convention that B(ζ, m, n) is empty if m > n. Lemma 77 Set F = def (G Σ ∪ G Σ ) \ Φ and F a = def (G a Σ ∪(G) a Σ ) \ Φ a .
Then: 1. F and F a rewrite as follows:

F = f ∈G B(f, k f +1, c f ) and F a = f ∈G B(f, K f +1, M c f ) . (168)
2. Let Z and Z a be the sets of dependent variables of F =0 and F a =0. Then:

Z = x∈X(G) B(x, d - x , d x ) and Z a = x∈X(G) B(x, M d - x , M d x ) . ( 169
) Proof Obvious. If c - f ≤ c f , then k f +1 = c - f and K f +1 = M c - f ; thus, B(f, k f +1, c f ) and B(f, K f +1, c f ) are both nonempty. If c - f > c f , then k f +1 = c f +1 and K f +1 = M c f +1; thus, B(f, k f +1, c f ) and B(f, K f +1, c f )
are both empty. Hence, (168) rewrites:

F = f ∈G:c - f ≤c f B(f, c - f , c f ) and F a = f ∈G:c - f ≤c f B(f, M c - f , M c f ) . (170)
Similarly, (169) rewrites:

Z = x∈X(G):d - x ≤dx B(x, d - x , d x ) and Z a = x∈X(G):d - x ≤dx B(x, M d - x , M d x ) . (171)
In words, F a and its set of dependent variables are obtained, from F and its set of dependent variables, via an M -stretching. From formulas (170) and (171), we deduce the following result:

Lemma 78 Let f ∈ G and x ∈ X(G). Then, (f •k , x •l ) is an edge of G F if and only if (f •M k , x •M l ) is an edge of G F a .
The situation is illustrated in Figures 24,25, and 26, for the clutch and two variants of it, for N = 1 and N = 2.

Procedure for constructing matchings of maximal cardinality for F and F a

We want to design a procedure for determining such matchings so that the two constructions parallel each other. This will allow us to compare the results.

Step 1: starting from G Σ and G a Σ By construction, G Σ and G a Σ are both structurally nonsingular and determine the leading variables of the new mode, given a consistent valuation of other variables. They both admit a complete matching.

Let N 1 be a complete matching for G Σ . Its edges have the form

f •c f , x •dx for appropriate pairs (f, x) ∈ G × X(G). Then, N a 1 = def f •M c f , x •M dx f •c f , x •dx ∈ N 1 is a complete matching for G a Σ . (172)
This parallel construction is illustrated in Figure 24.

Inria 24, the added edges correspond to consistency equations and are colored in blue. Subsequently applying the Dulmage-Mendelsohn decomposition returns β O = ∅ on both. Actually, the two matchings are complete. As a result, the mode change is fully determined, despite the new mode for the clutch having underdetermined consistency conditions (because the common rotation velocity is free).

Step 2: Removing, from the set of dependent variables of G Σ and G a Σ , the variables involved in the contexts ∆ and ∆ a Let

Y = def x ∈ X(G) x •dx ∈ X(∆) and Y a = def x ∈ X(G) x •M dx ∈ X(∆ a ) . Then, Y = Y a holds. Denote by G 2 (resp. G a 2 ) the system G Σ (resp. G a Σ )
where the variables involved in the context ∆ (resp. ∆ a ) have been removed. The bipartite graphs of G 2 and G a 2 , denoted by G 2 and G a 2 , are isomorphic. Hence: If N 2 denotes a matching of maximal cardinality of G 2 , then,

N a 2 = def f •M c f , x •M dx f •c f , x •dx ∈ N 2 is a matching of maximal cardinality for G a 2 . (173)
In particular, their Dulmage-Mendelsohn decompositions are isomorphic too. An example of this is illustrated in Figure 26. We start from N 2 and N a 2 , perform their respective Dulmage-Mendelsohn decompositions, and remove, from both bipartite graphs, their respective overdetermined blocks β O and β a O . The resulting bipartite graphs are denoted by G 2 and G a 2 : they are also isomorphic. Hence: Let N 2 be a matching of maximal cardinality for G 2 . The formula

N a 2 = def f •M c f , x •M dx f •c f , x •dx ∈ N 2 defines a matching of maximal cardinality for G a 2 . ( 174
)
RR n°9334 Step 3: Adding to G Σ the consistency subsystem G Σ and subtracting Φ; performing the corresponding changes to G a Σ We start from G 2 , N 2 and G a 2 , N a 2 . We do not need to remove equations belonging to Φ and Φ a , since such equations, if any, have already been removed at Step 2, when moving from N 2 to N 2 . Therefore, the only final modification we need to perform is adding the consistency equations. The following result holds:

Lemma 79 1. Let f ∈ G such that c - f ≤ c f and let (f •c f , x •dx ) ∈ N 2 . For every m such that 0 < m ≤ c f -c - f , the edge (f •(c f -m) , x •(dx-m)
) is an augmenting path of N 2 with respect to system F .

Let

f ∈ G such that c - f ≤ c f and let (f •M c f , x •M dx ) ∈ N a 2 . For every m such that 0 < m ≤ M (c f -c - f ), the edge (f •(M c f -m) , x •(M dx-m)
) is an augmenting path of N a 2 with respect to system F a .

Proof Consider Statement 1. If c f ≤ c - f for every f ∈ G, then matching N 2 is of maximal cardinality for F , so that no augmenting path exists. From now on, we assume c f > c - f for some f ∈ G. Select any such f and let (f

•c f , x •dx ) be the edge of N 2 involving f . Then, f •(c f -1) and x •(dx-1) are involved in F , but x •(dx-1) is not a vertex of N 2 . Thus, for every f ∈ G such that c f > c - f , then (f •(c f -1) , x •(dx-1)
) is an augmenting path of N 2 in F . The same reasoning applies for (f

•(c f -m) , x •(dx-m) ), provided that 0 < m ≤ c f -c -
f . This proves Statement 1. The same reasoning applies for proving Satement 2.

Lemma 80

1. The following formulas define matchings of maximal cardinality for F and F a , respectively:

M = f •(c f -m) , x •(dx-m) f •c f , x •dx ∈ N 2 , 0 < m ≤ c f -c - f M a = f •M (c f -m) , x •M (dx-m) f •c f , x •dx ∈ N 2 , 0 < m ≤ c f -c - f 2.
Consider the surjective mapping ψ : M a → M, defined by

ψ f •M (c f -m) , x •M (dx-m) = f •(c f -m) , x •(dx-m)
Inria and define

ψ f : F a → F by ψ f (f •M (c f -m) ) = f •(c f -m) , ψ x : X(F a ) → X(F ) by ψ x (x •M (dx-m) ) = x •(dx-m) .
Then

(f •(c f -m) , x •(dx-m) ) ∈ M if and only if (f •M (c f -m) , x •M (dx-m) ) ∈ M a . Let g = f •(c f -m)
be a vertex of G F associated with a function. Then, g is unmatched in M if and only if ψ f (g) is unmatched in M a . We conclude by invoking Lemma 78: the mappings ψ x and ψ f preserve the overdetermined blocks. A similar reasoning holds regarding the preservation of underdetermined blocks, by starting from unmatched vertices of variable type.

This finishes the proof of Theorem 74.

Comparing the standardizations

The independence of the standardization with respect to the particular semantics of the derivative is immediate within long (continuous) modes. We thus focus on mode changes, where we must target a standard discrete-time dynamical system indexed by N. Throughout this section, Assumptions 6 (Page 73), 8 (Page 83) and 9 (Page 86) are in force. We want to compare the respective standardizations of the returns of Algorithm 10 (ExecRun), when using the two expansions (164) and (165) for the derivatives. We first adapt the generic system (124), computing the restart at a mode change:

expand X as 1 ∂ N M α n (X • -X) •n in 0 = H(X , X • , V, X) , (175)
also written 0 = H(Z, X, ∂) .

In (175), X collects the state variables, while V collects the algebraic dependent variables, decomposed into impulsive and non-impulsive variables: V = W i ∪ W . Since N ≥ 0, the dependent variables of system (175) are Z = def (X •(N +1) , V ), and we set Y = def (X •(N +1) , W ); in words, Y collects the non-impulsive dependent variables. In (175), the functions F and G are standard; so is the forward shift x → x • , since it will be mapped to a corresponding forward shift in the discrete time index N. Hence, the infinitesimal time step ∂, arising in the expansion of the derivative in F , is the only reason for System (175) to be nonstandard.

Assumption 8 states that impulsive variables can be eliminated from system (175), resulting in the following reduced system:

K(X , X • , W, X) = 0 ⇐⇒ ∃W i .H(X , X • , V, X) = 0 . ( 177
)
Recall that ∂ enters System (177)-left, i.e., the left-hand part of the equivalence given by Equation (177), via the expansion used for X in (175). As part of Assumption 8:

System (177)-left remains structurally nonsingular when ∂ = 0.

In (177), (X •(N +1) , W ) standardizes as the right-limit (X + , W + ) whereas all the X •n , n = P, . . . , 0 standardize as the left-limit X -. Therefore, in the expansion used for X in (175), all terms cancel out but the last one. Hence, System (177) standardizes as:

• Case N > 0: substitute ∂ ← 0 in K α N ∂ (X + -X -), X -, W, X = 0; • Case N = 0: substitute ∂ ← 0 in K α N ∂ (X + -X -), X + , W, X = 0.
The difference lies in the X -vs. X + , highlighted in red. These substitutions are both legitimate, thanks to (178). In both cases, the result does not depend on α N , since α N = 0. Hence, we get the following result:

Theorem 81 Under Assumptions 6, 8, and 9, the standardization of the restart conditions does not depend on the expansion (165) for the derivative.

Comment 82 So far, Theorem 81 applies to mDAE systems possessing only long modes. We could extend this theorem in several directions. First, we could include in our modeling language leftor right-limits: x -or x + . Then, our reasoning would apply with no change provided that x -and x + are encoded as • x or x • if expansion (164) is used, and

A correctness result

Our approach, consisting of the structural analysis in the nonstandard domain (Algorithm 10 and its refined version Algorithm 13) followed by standardization, either fails, in which case it returns proper diagnosis, or succeeds and produces executable simulation code. The natural question is:

Question 83 Does this approach compute the solution of the given mDAE system?

Once again, answering Question 83 in its full generality is out of reach as, to the best of our knowledge, there is no general mathematical definition of the solution of an mDAE system. Pathologies may exist in the complement of the long 'continuous' modes: possibly unbounded cascades of events, chattering or sliding modes, Zeno behaviors for events of mode changes or transient modes, and more. Even when the set of events of mode changes is 'gentle' (a finite or diverging sequence of instants), the characterization of the impulses may not be available. Still, as mentioned in the Introduction, notions of solutions of mDAE systems exist for some restricted classes.

In this section, we provide an answer to Question 83 for a restricted subclass of mDAE systems, by reproducing the results of [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF] and giving their detailed proofs. For the description of the class of systems considered, we use the notion of function of bounded variation: a function f : R → R has bounded variation if it is the primitive of a Lebesgue integrable function [START_REF] Dunford | Linear Operators, Part I, General Theory[END_REF]. As a consequence, if f has bounded variation, then lim

h 0 t+h t f (s)ds = 0 . (179)
Definition 84 (semilinear systems) Call semilinear an mDAE system S such that, for each mode µ, the active DAE system takes the restricted form

0 = A(X s )X + B µ (X) 0 = C µ (X) (180)
where:

1. X s collects the smooth elements of X, i.e., those being continuous and of bounded variation around each instant of mode change. Other elements of X might be discontinuous.

2. Both matrix A(.) and the mode-dependent vectors B µ (.) and C µ (.) are continuous functions of their arguments, and remain bounded in a neighborhood of each instant of mode change.

The matrix

J = A(X s ) J X C µ (X)
is regular, where J X C µ (X) denotes the Jacobian matrix of C µ with respect to X, evaluated at X.

Inria

Comment 85 The special form (180) and Condition 3 for the Jacobian matrix express that mode dynamics have been processed using the Σ-method: differentiating the algebraic constraint once yields a regular linear system determining X . Then, Conditions 1 and 2 together make the impulse analysis easier.

Comment 86 In [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF], Section 4.3, it is shown that multi-body systems with contacts yield semilinear mDAE systems.

For semilinear mDAE systems, the solution at an instant t * of mode change can be determined as follows. After proper reordering of the entries of X, we decompose vector X into its smooth and nonsmooth parts and decompose matrix A accordingly:

X = X s X n and A = A s A n .
Since X s is smooth at any instant, by Condition 1, matrix A X s (t * ) is known prior to computing the mode change. In contrast, X n (u)du is a Dirac measure at mode change.

To evaluate the restart values for all elements of X, we integrate the system dynamics over the interval [t * -ε, t * + ε] for ε > 0 sufficiently small, taking into account the dynamics before and after the mode change. In the following calculations, µ -and µ + denote the mode prior and after the mode change, and B + and C + stand for B µ + and C µ + , respectively.

[t * -ε,t * ) A(X s (u))X (u)+B -(X(u)) du + A s (X s (u))X s (u)+B + (X(u)) du

+ [t * -ε,t * +ε]
A n (X s (u))X n (u)du ≈ 0 + 0 + A n (X s (t * )) X n (t * +ε) -X n (t * -ε) ≈ A(X s (t * )) X(t * +ε) -X(t * -ε)

The evaluation to zero of the first two integrals follows from (179). For the third integral, we use the fact that X n (u)du is approximately a Dirac measure at t * in the indicated interval. This reasoning leads to the following result determining the restart conditions at t * :

Theorem 87 ( [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF]) The restart conditions X + are determined from the left-limits X -by using the following system of equations:

0 = A(X - s )(X + -X -) 0 = C + (X + ) ( 181
)
which is locally regular thanks to Condition 3 of Definition 84.

Whether we can solve system (181) depends on the particular system at hand. It remains to compare scheme (181) with our approach. The nonstandard semantics of (180) is:

0 = A(X s )(X • -X) + ∂ × B µ (X) 0 = C µ (X)
and applying the Σ-method to it yields the augmented system

   0 = A(X s )(X • -X) + ∂ × B µ (X) 0 = C µ (X) 0 = C µ (X • ) . ( 182
)
We briefly recall the corresponding notations for convenience: in (185), X collects the state variables and Z = def (X • , V ) collects the dependent variables. The algebraic variables V can be split into impulsive and non-impulsive variables: V = W i ∪ W , and we set Y = def (X • , W ). Assumption 8 states that impulsive variables can be eliminated from System (186), resulting in the reduced system K(Y, X, δ)=0 ,

that remains structurally nonsingular when δ = 0.

Theorem 91 The following statements hold:

1. Let Y * be a solution of (187). There exists a selection map δ → (Y (δ), W i (δ)), from R >0 to the set of solutions of the system (186), such that lim δ 0 Y (δ)=Y * .

2. Conversely, let (Y (δ), W i (δ)) be a solution of (186) such that Y * = def lim δ 0 Y (δ) exists; then, Y * is a solution of (187).

Proof We successively prove the two statements.

Statement 1. By the Implicit Function Theorem applied to the standard function K, there exists a function ( X, δ) → Y ( X, δ) of class C 1 such that Y (X, 0)=Y * and K(Y ( X, δ), X, δ)=0 holds for ( X, δ) ranging over a neighborhood of (X, 0). Then, by Assumption 8, ∃W i .H(Z, X, δ)=0 ⇐⇒ K(Y, X, δ)=0 holds, which concludes the proof of this statement.

Statement 2. This is trivial and does not use Assumption 8. Since K is smooth, the two conditions (1) K(Y (δ), X, δ)=0 for every δ > 0, and (2) Y (δ) converges to Y * , imply K(Y * , X, 0)=0.

Corollary 92 Let Assumption 8 be in force and let the systems (186) and (187) possess a unique solution for Z and Y , respectively. Then, H(Z(δ), X, δ)=0 implies K(lim δ 0 Y (δ), X, 0)=0.

For its practical application, this scheme requires the prior knowledge of the decomposition V = W i ∪ W of algebraic variables into impulsive and non-impulsive ones: we solve System (186) for Z and simply discard the values of W i from this solution. The drawback of this scheme is that system (186) will be likely ill-conditioned.

Thus, in practice, we would recommend a rescaling of the impulsive variables based on the impulsion order: for every impulsive variable w or order µ = def [[w]], we perform the substitution w ← ν × δ -µ , where ν is the new variable replacing w-such a rescaling works only for variables of finite impulsion order, see Section 10.1.

Experimental results: mode changes of a nonsemilinear system

To exercise a numerical method computing state-jumps at mode changes, when some variables are impulsive, let us consider the following multimode DAE system:

               ω 1 = a 1 ω 1 + b 1 τ 3 1 (e 1 ) ω 2 = a 2 ω 2 + b 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 )

(188

)
This system is identical to the clutch model (3), with the only difference that torque τ 1 appears non-linearly in equation (e 1 ). Therefore, this system is not semilinear, and the direct approach based on linear-algebraic techniques used in Section 11.2 does not apply. Nevertheless, we shall demonstrate on this example that our general approach applies. We focus on the mode change

We now show that solving System (189) with a small standard positive parameter instead of the nonstandard ∂ yields a numerical scheme converging to the solution (194), if we disregard impulsive variables. Algorithm 14 describes our numerical method, whereas Figure 27 x ← π(u); h ← θh 6: until |x -y| ∞ ≤ def n_error(eps):

(n,u,q) = approx(eps,1e-2,0.5) return (q[0] -5) Experimental results on this method have been obtained using Sage. The nonlinear solver used is the least squares solver of the Scipy's optimize library. Figure 27 gives, on the left, the number of iterations required, and on the right, the absolute error η = |x -x| ∞ between the computed solution and the exact solution x = (ω 2 , ω 2 ). The scaling factor is θ = 0.5, the initial time step h 0 = 10 -2 and u 0 = (ω 1 , ω 2 , 0, 0). It can be observed that the method converges in not more than 37 iterations, down to = 10 -9 . It should also be remarked that, experimentally, the absolute error is bounded : η < 2 .

12 Toward a tool supporting our approach

The general approach we presented sets formidable challenges for a tool development. We identify three important requirements for such a tool.

• First, it should implement the theory: structural analysis must be adjusted to every different mode and mode change. Both aspects are included in what we call a 'multimode structural analysis'.

• Second, the described approach aims at the compilation of mDAE systems. As such, 'on-the-fly' techniques for the structural analysis of such systems, as advocated in [START_REF] Broman | Modelyze: a gradually typed host language for embedding equation-based modeling languages[END_REF] and [START_REF] Elmqvist | Systems modeling and programming in a unified environment based on Julia[END_REF], cannot be taken into consideration. As soon as a model is structurally unsound in any given mode or at any given mode change, it should be rejected at compile time. For accepted models, the output of the tool should enable the generation of efficient simulation code.

• Third, enumerating all modes and mode changes to perform multimode structural analysis is bound to collapse. An mDAE compiler should be able to handle models in a clever way to avoid such enumeration.

We believe that we achieved an important first step towards alleviating these challenges with our tool IsamDAE under development. IsamDAE implements the structural analysis of all modes of an mDAE without any enumeration of the modes, instead relying on Binary Decision Diagrams (BDD, see for instance [START_REF] Bryant | Graph-based algorithms for Boolean function manipulation[END_REF]) for implicit representations of the mode-dependent structure of an mDAE. This is described more thoroughly in [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF][START_REF] Caillaud | Implicit Structural Analysis of Multimode DAE Systems[END_REF].

In order to illustrate how a multimode model is currently handled by IsamDAE, this Section focuses on the so-called RLDC2 model, a simple electronic circuit shown in Figure 28. This example, provided to us by Sven-Erik Mattsson, is actually a difficult one to handle for the existing Modelica tools because of its mode-dependent index and structure, even though it exhibits no impulsive behavior. Indeed, we show once more how this model is mishandled by two leading Modelica tools, namely, Dymola and OpenModelica, and hint at reasons for this fact, thereby justifying our work once again. We then detail its handling by the IsamDAE tool, in order to provide an insight about why this tool could be expanded into a complete multimode structural analysis suite.

The RLDC2 model

The circuit consists in two RLC circuits interconnected in parallel and in which two diodes have been inserted. The diodes are considered to be ideal, meaning that they are not ruled by the customary Shockley law i = I e u/U -1 , but rather by the complementarity condition 0 ≤ i ⊥ -u ≥ 0, meaning that, at any moment, both i and -u are nonnegative and iu = 0 (Section 1.2.2 already mentions complementarity conditions). Figure 29 illustrates these two diode laws.

As in Section 1.2.2, we redefine the graph of this complementarity condition as a parametric curve, represented by the following three equations:

s = if γ then i else -u (S) 0 = if γ then u else i (Z) γ = (s ≥ 0) (G)
Inria To avoid a logico-numerical fixpoint, we consider instead the system obtained by replacing s with its left-limit s -in equation (G):

s = if γ then i else -u (S) 0 = if γ then u else i (Z) γ = (s -≥ 0) (G -)
Under the assumption that u and i are continuous functions of time (which turns out to be a valid assumption for the RLDC2 circuit), s(t) = s -(t) holds at every instant t. This implies that systems (S), (Z), (G) and (S), (Z), (G -) are equivalent. Using this encoding of the complementarity condition, the RLDC2 circuit is modeled as the mDAE shown in Figure 30, with two Boolean variables γ 1 and γ 2 and four multimode equations (S 1 ), (Z 1 ), (S 2 ) and (Z 2 ), incident to a set of variables that depends on the mode of the system (i.e., the values of Boolean variables γ 1 and γ 2 ). Note that this model belongs to the class introduced in Definition 33. 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 ) (195)

Handling the RLDC2 model with Modelica tools

The Modelica model of the RLDC2 circuit is given in Figure 31; it is a direct translation of the model above. The authors had the opportunity to test it on two implementations of the Modelica language: OpenModelica v1.12 26 and Dymola 2019. 27 Identical results are obtained with both tools: the model is deemed nonsingular at compile time, yet a runtime error immediately occurs (at time t = 0) during simulation, as shown in Figure 32.

The issue is that the structural analysis implemented in these tools treats the model as a single-mode DAE, disregarding the mode-dependent variability of the incidence relations. The consequence is that, despite the model being deemed structurally nonsingular by these tools, blocks of equations that are structurally singular in some modes are produced. This immediately leads to runtime errors: both tools attempt a pivoting of the Jacobian matrix by an element that is equal to zero in some modes.

Structural analysis of the RLDC2 model

The mDAE system of Figure 30 has two modes of nonzero index, when diodes are either both passing or both open. Having performed index reduction on each of them, we show the corresponding block dependency graphs in Figures 33 and34, respectively-see Lemma 14 and Figure 19. These two graphs significantly differ, highlighting the fact that a mode-dependent structural analysis is necessary.

Figure 35 shows a graph generated by the structural analysis performed by our software Isam-DAE [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF][START_REF] Caillaud | Implicit Structural Analysis of Multimode DAE Systems[END_REF], for the mDAE system of Figure 30. This graph does not consist of one graph per mode, but, rather, of a set of edges labeled by a propositional formula characterizing the set of modes in which the considered branch is involved. This data structure prevents the combinatorial explosion arising in mDAE systems composed of a large number of mDAE subsystems, when modes are enumerated.

The block triangular form of the system can then be found at runtime for any given mode, for instance by evaluating the block dependency graph in this mode and performing a topological sort on the resulting graph. For the RLDC2 example, when evaluated in the mode in which both diodes are passing (resp. blocking), the graph of Figure 35 yields the dependency graph of Figure 33 (resp. Figure 34); the same holds in the remaining modes. All equation blocks can be turned into efficient simulation code at compile time, so that the computational overhead due to a mode switching at runtime is minimized.

Perspectives

As shown on the RLDC2 example above, our tool IsamDAE addresses the issue of performing the structural analysis of all modes of an mDAE while avoiding any explicit enumeration of these modes. Although the example detailed above is, for clarity purposes, a small one, IsamDAE is provably able to handle models with a very large number of modes, as detailed in [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF][START_REF] Caillaud | Implicit Structural Analysis of Multimode DAE Systems[END_REF].

As such, the tool already provides the data structures and key algorithms required for addressing the second and third challenges invoked at the beginning of this section. This encourages us to implement, in the same tool, the structural analysis of mode changes, in the same 'all-at-once' fashion as for the modes themselves. This is indeed the main perspective of our work, as we believe it will be a crucial step towards the design of a mathematically sound compiler for multimode DAE systems.

Conclusion

We have proposed sound foundations for the compilation and code generation for physical systems modeling languages relying on DAE for their continuous dynamics; major instances are Modelica and Simscape. As its main contribution, our theory explicitly considers events of mode change and submitting her/his model.

A second requirement was the ability to handle the conflicts that may result, at mode changes, between the dynamics of the previous mode and the consistency conditions implied by the new mode.

We do not see how these two requirements could be addressed without relying on nonstandard analysis. Mapping real time to a discrete time with infinitesimal step size allowed us to cast both (DAE-based) continuous dynamics, and the restart conditions at mode changes, to a uniform framework. This allowed us to properly solve the conflicts mentioned above, between the dynamics of the previous mode and the consistency conditions implied by the new mode. To make all of this doable, we developed a toolbox of new results on nonstandard analysis, mainly related to the needs of structural analysis. We proved that the result of our compilation remains independent from the particular scheme we used to map derivatives to the nonstandard domain.

From our experience in using nonstandard analysis, we can state that this is an extremely convenient tool to establish the mathematical foundations for the compilation of continuous-time systems modeling tools, for both the Simulink and Modelica classes.

A desirable objective of a work like ours is to prove that every execution scheme produced by the approach actually generates solutions of the given source mDAE system. As a reference, such soundness analyses were performed for all synchronous languages [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF][START_REF] Benveniste | Compositionality in dataflow synchronous languages: Specification and distributed code generation[END_REF][START_REF] Benveniste | The synchronous languages 12 years later[END_REF]. However, a particular difficulty holds: no mathematical definition of solutions of an mDAE system exists in full generality, , where: p is a propositional formula defining in which modes the block is evaluated; R is the set of variables to read; B is the set of equations of the block; W is the set of variables to write. Edges are labeled by a propositional formula, defining in which modes the dependency applies-notation "!g" means "not g".

Inria

that would serve as a reference. So far, we were only able to prove such a correctness result for the subclass of "semi-linear multimode systems", that may involve impulsive behaviors at mode changes.

In this work, we paid no specific attention to the efficiency of our structural analysis algorithm, and did not discuss how to generate scalable simulation code. We also only considered flat guards, whereas nested guards (naturally occurring, for instance, in the case of nested "if-then-else" constructs) are very useful in a practical formalism: flattening a model with nested guards is not desirable for getting efficient code. Work is in progress to address these issues by relying on the technique of Binary Decision Diagrams (BDD), and first results are reported in [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF], where our IsamDAE tool under development is described.

A Appendix: constructive execution of an instant in systems having only long modes

In this appendix, we briefly develop the execution of a run, for mDAE systems according to Definition 27, i.e., without delaying the effect of guards by one nonstandard instant. Since the effect of guards may not be postponed, the current mode is not known from evaluating the previous instant. We still refuse to solve mixed logico-numerical fixpoint equations. Hence we must instead identify the current mode progressively and constructively, according to an iterative scheme of the following form:

1. Enter the new nonstandard instant; the status of a subset of the S-variables is known from previous instant;

2. Using the subset of evaluated guards, select the corresponding set of enabled equations and submit them to index reduction, conflict solving, and then evaluate them;

3. This yields more enabled guards; evaluate them and return to Step 2.

This iteration is performed until, either all S-variables are evaluated to t, f, or i, or no additional guard gets enabled at Step 3-a failure. Since knowing if the current mode is long or transient modifies the way index reduction is performed, and since knowing this requires having evaluated all the guards, we are unable to process general mDAE systems, and we must restrict ourselves to systems having only long modes.

To execute a nonstandard instant through a sequence of microsteps, the algorithm ExecRunProgress (Algorithm 15) is iterated until a failure is reported or a Tick atomic action occurs. In case of success, the pair (σ, ∆) produced by Tick produces the initial status and the context for starting the next instant. ExecRunProgress requires a finite coherent status σ and a finite context ∆; it returns Algorithm 15 ExecRunProgress; compare with Algorithm 10 Require: (σ, ∆); return (Fail , updated (σ, ∆))

1: if Success(σ) then 2:

(σ, ∆) ← Tick(σ) Line 2 If σ is complete, then Tick (Algorithm 7) is performed and the instant is completed.

Line 3 Not all the equations belonging to the context are active in the current instant. A system Φ collecting all the facts is introduced and is initialized to the subset of the context ∆ consisting of all its enabled equations, taking into account the guards that are already known to be true at initialization of the run.

Line 5

The system G is set to the enabled guarded equations not evaluated yet, in status σ.

Line 6 We submit G to the index reduction (IndexReduc atomic action, Algorithm 8). This returns, Inria

Figure 1 :

 1 Figure 1: An ideal clutch with two shafts.

Figure 2 :

 2 Figure 2: Modelica code for the idealized clutch.

Figure 3 :

 3 Figure 3: Division by zero exception occuring when simulating the Ideal Clutch Modelica model.

Figure 4 :

 4 Figure 4: Trajectory of the Ideal Clutch Modelica model: it stops at t = 5s, when the exception occurs.

Figure 5 :

 5 Figure 5: Solution of the Ideal Clutch model, computed by Mathematica.

Figure 6 :

 6 Figure 6: The Cup-and-Ball game.

Figure 7 :

 7 Figure 7: Westinghouse brake for one car: schematics.p b : connector pressure (a control) p bn : next connector pressure p r : railcar pressure p t : reservoir pressure f b : mass flow brake connector f bn : mass flow next brake connectorf bn = 0 f v : mass flow one way valve f l : leakage mass flow f cl : mass flow low press side of cylinder f ch : mass flow high press side of cylinder f t : mass flow reservoir b : brake force x : position of the piston

Figure 9 :

 9 Figure 9: Simulation of the clutch model with resets, using our method. Mode change f → t occurs at t = 5s and mode change t → f occurs at t = 10s.

Figure 16 :

 16 Figure 16: Restart at mode change γ : t → f.

Figure 17 :

 17 Figure 17: Restart at mode change γ : f → t.

Figure 18 :

 18 Figure 18: Our approach to the compilation of mDAE systems.

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

 191212 Figure 19: Illustrating Lemma 14 (we only show the projection of strongly connected components on function nodes).This lemma is illustrated inFigure 19, inspired by lecture notes by J-Y. L'Excellent and Bora Uçar, 2010. In this figure, a same bipartite graph G is shown twice (left and right), with the equations sitting on the left-hand side in green, and the variables on the right-hand side in red. Two complete matchings M 1 (left) and M 2 (right) are shown in thick blue, with other edges of the bipartite graph being black. The restriction, to the equation vertices, of the two directed graphs G M 1 F (left) and G M 2 F (right) are shown on both sides. Although the two directed graphs G M 1 F and G M 2 F differ, the resulting block structures (encircled in dashed brown) are identical.

Lemma 17

 17 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 .

Definition 27 (

 27 mDAE/mdAE syntax) Multimode DAE systems (mDAEs), respectively multimode dAE systems (mdAEs), are defined by the following syntax: S ::= s | S; s s ::= g | e g ::= γ = b (67) e ::= |if γ then f =0 (68)

Figure 20 :

 20 Figure 20: ExecRun as a labeled automaton for the Simple Clutch. Notations: For all statuses (shown in boxes), v (respectively v) means v = t (respectively v = f), and not mentioning v means v = u. e means that e f belongs to context ∆. The assignment ∆ ← e 3 refers to Line 10 of Algorithm 10. Blue (respectively black) transitions belong to a continuous-time (respectively discrete-time) dynamics. The red transition is impulsive. A semicolon is the sequential composition of computations, and the + symbol defines blocks of equations.

  Elementary cycles of ExecRun capture runs with stationary valuations of the guards and define the continuous dynamics in each mode. Other runs capture mode changes and their reset actions. During runs with constant valuations of the guards, we have (G Σ ∪ G Σ ) \ Φ = G Σ . Hence, Line 1 returns β U = β O = ∅ and, thus, β E = G Σ : we recover the usual structural analysis for DAE systems.

Lemma 49 Assumption 5

 495 Assumption 5 implies Assumption 4, where Assumption 5 states: For every subset Y of reset variables, d o (X(M (Y ))) > d o (Y ) where d o (Z) is the smallest (shifting or differentiation) degree among all variables of Z.

10. 1 . 1 Figure 21 :

 1121 Figures 21 and 22 display the rules defining the translation of a system of equations of the form (124) or (125) into its impulse analysis, for the restricted class where only rational expressions are involved. e ::= 0 | c | ∂ | x | e c | e + e | e × e E ::= e = e | E and E

Figure 22 :

 22 Figure 22: Rules: The left column displays the impulse order of the primitive expressions. Rule (R7) indicates that [[e]] = [[e ]] is an equation of the impulse analysis [[E]] if e = e is an equation of E; Rule (R8) indicates that, if E involves the equation x = y + e but not the equation y = x -e, then we augment E with the latter, i.e., we saturate E with the rule x = y + e =⇒ y = x -e.

10. 1 . 3

 13 Example: mode change γ : t → f of the clutch

  which proves [[τ 1 ]] = [[τ 2 ]] = 1, expressing that the two torques are impulsive. Note, on the other hand, that (141) does not fully determine the impulse orders of ω • i and u i , since one may, e.g., have [[ω • i ]] < 0 and [[u i ]] = 0. Inria 10.1.

  ( x, ∂) = ([x n ], [δ n ]) be such that st( x, ∂)=(x, 0). Then, by definition of the standard part, there exists B ⊆ N with µ(B) = 1 such that, for every n ∈ B, (x n , δ n ) ∈ V . The sequence (G(x n , δ n )) n∈B satisfies H(G(x n , δ n ), x n , δ n )=0, for every n ∈ B. (147) We extend the sequence (G(x n , δ n )) n∈B to N by assigning arbitrary values to indices n ∈ N \ B and consider the nonstandard z = def [G(x n , δ n )]. By (147), we get H ( z, x, ∂) =0.

  3.5. (167-b) and (167-c) follow from Definition 40 and Lemma 73. For the second statement, we note that, since M

1 Figure 24 :

 124 Figure 24: The clutch is an example in which c - f = c f holds. Maximal cardinality matchings for N = 1 (left) and N = 2 (right) are given. Subsequently applying the Dulmage-Mendelsohn decomposition returns β O = {e 3 } on the left, and β O = {e 3 , e • 3 } on the right.

1 Figure 25 :

 125 Figure 25: An example with c - f < c f : a variation on the clutch. The previous mode does not involve ω • 1 , hence ω 1 remains a dependent variable when entering the new mode. Maximal cardinality matchings for N = 1 (left) and N = 2 (right) are shown. With reference with the example of Figure 24, the added edges correspond to consistency equations and are colored in blue. Subsequently applying the Dulmage-Mendelsohn decomposition returns β O = ∅ on both. Actually, the two matchings are complete. As a result, the mode change is fully determined, despite the new mode for the clutch having underdetermined consistency conditions (because the common rotation velocity is free).

1 Figure 26 :

 126 Figure 26: An example with c - f > c f : another variation on the clutch. The previous mode determines ω •2 2 , hence ω • 2 is no longer a dependent variable when entering the new mode. Maximal cardinality matchings for N = 1 (left) and N = 2 (right) are shown. Subsequently applying the Dulmage-Mendelsohn decomposition on the left-hand graph returns β O = G F , hence Line 6 of Algorithm 10 returns H = ∅. No equation of the new mode is satisfied at the first instant and we move to the next instant. The same holds for the right-hand graph.

  [t * ,t * +ε] A(X s (u))X (u)+B + (X(u)) du = [t * -ε,t * ) A s (X s (u))X s (u)+B -(X(u)) du + [t * ,t * +ε]

  displays the results. Algorithm 14 Numerical method for approximating the state-jump of non-semilinear system (188) for mode change γ : f → t; F h (u) denotes the value of the right-hand side of System (189) for a 4-tuple u of values for the unknowns and ∂ replaced by h. Require: Difference equations F h (u) : R×R N →R N , u 0 ∈R N initial guess, π : R N →R M projection on the set of non-impulsive variables, h 0 >0 initial time step, >0 tolerance bound, 0<θ<1 scaling factor ; return u estimated state-jump 1: u ← u 0 ; x ← π(u 0 ); h ← h 0 NSolve(system = F h , guess = u, tolerance = h) 5:

Figure 27 :

 27 Figure 27: Numerical method for computing the state-jump of non-semilinear system (188): (a) number of iterations as a function of ; (b) error η as a function of .

Figure 28 :

 28 Figure 28: Schematics of the RLDC2 circuit.

Figure 29 :

 29 Figure 29: Shockley law (dashed, magenta) vs. ideal complementarity condition (solid, green).

Figure 30 :

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

Figure 33 :

 33 Figure 33: Block dependency graph for the RLDC2 example with both diodes passing (i.e., γ 1 = γ 2 = T).

Figure 34 :

 34 Figure 34: Block dependency graph for the RLDC2 example with both diodes blocking (i.e., γ 1 = γ 2 = F).

Figure 35 :

 35 Figure35: Block dependency graph of the RLDC2 model, generated by IsamDAE. Vertices are labeled p : R → B → W , where: p is a propositional formula defining in which modes the block is evaluated; R is the set of variables to read; B is the set of equations of the block; W is the set of variables to write. Edges are labeled by a propositional formula, defining in which modes the dependency applies-notation "!g" means "not g".

3 :ΦLine 1

 31 ← {e ∈ ∆ | σ(γ(e)) = t} 4: else 5:G ← [Enab(σ) ∩ Undef (σ)] 6: (σ, b, G Σ , G Σ ) ← IndexReduc(G, σ); increase σ 7: if b = f then return Fail (σ) b, H) ← SolveConflict((G Σ ∪ G Σ ) \ Φ, σ); increase σ 10: if b=f then return Fail (σ) 11:else (σ, Φ) ← Eval(∆, Φ, σ, H); increase σ updated σ and ∆, or documented failure information. Its details are commented next. Function Success(σ) checks if status σ is complete according to Definition 39.

1 :

 1 Solve {e 1 , e 2 , e 3 , e 4 } for (ω 1 , ω 2 , τ 1 , τ 2 ) 4 equations, 4 unknowns 2: ODESolve (ω 1 , ω 2 )

  • γ = t and γ = f, where • γ denotes the backward shift of γ following Definition 2. Unfolding System (17) 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

  Observe that Exec. Sch. 4 would work without changes if the guard γ was a predicate on the state variables ω 1 and ω 2 .

		e ∂ 2 , e 5 , e 6
	10:	Tick: ∆ ← ∆ \ {e 3 }
		Inria

  Transition γ : f → t This one is more involved. As established in Exec. Sch. 4, in order to compute the reset values, we use the system of 4 equations {e ∂

	variables (τ 1 , τ 2 , ω • 1 , ω • 2 ). In particular, from e ∂ i , we get	1 , e ∂ 2 , e • 3 , e 4 } to determine the leading

  e ∂ 2 , e 3 , e 4 , e 5 , e 6

		e 5 ; e 6 ;
		e ∂ 1 ; e ∂ 2
	γ; e 3 ; e 4	Tick
	Tick	γ; e 5 ; e 6
	γ, ω 1 , ω 2 ,
	τ 1 , τ 2 , ω • 1 , ω • 2 , e ∂ 1 , e ∂ 2 , e 3 , e 4 , e 5 , e 6		γ, ω 1 , ω 2 , e 5 , e 6	e ∂ 1 + e ∂ 2 + e • 3 + e 4
	e 5 ; e 6 ; e ∂ 1 ; e ∂ 2	e 3 , e 4 γ, ω 1 , ω 2 ,		γ, ω 1 , ω 2 , τ 1 , τ 2 , ω • 1 , ω • e ∂ 1 , e ∂ 3 , 2 , e • 2 ,
				e 4 , e 5 , e 6
		γ; e 3 ; e 4		Tick
		ω 1 , ω 2 , e 3
		Tick		γ; e 5 ; e 6 ; ∆ ← e 3
		γ, ω 1 , ω 2 ,
		τ 1 , τ 2 , ω • 1 , ω • 2 , e ∂ 1 , e ∂ 2 , e 3 , e • 3 , e 4 , e 5 , e 6	e ∂ 1 + e ∂ 2 + e • 3 + e 4	γ, ω 1 , ω 2 , e 3 , e 5 , e 6

  Call infinitesimal any nonstandard real number whose absolute value is smaller than any positive real number: [ |x n | ] ≤ [ε] for any ε ∈ R, ε > 0. For x and y two nonstandard real numbers, write x ≈ y if x -y is infinitesimal. Call infinite any nonstandard real number whose absolute value is larger than any positive real number: [ |x n | ] ≥ [K] for any finite K ∈ R. Call finite a nonstandard real number that is not infinite. Call bounded a subset A ⊆ R such that there exists a finite standard positive real K such that |x| ≤ K for every x ∈ A. Lemma 59 (standard part) For any finite [x n ] ∈ R, there exists a unique standard real number x ∈ R, such that [x] -[x n ] is infinitesimal. We call x the standard part of [x n ] and denote it by st([x n ]). Infinite nonstandard reals have no standard part in R.

  where X∈ R m , Z collects the dependent (nonstandard) variables, H=[H, H, H, . . . ] is the nonstandard lifting of the standard function H, and ∂=[δ n ] is infinitesimal.Definition 65 Say that H (Z, X, ∂) = 0 is structurally nonsingular if and only if so is the standard equation H (Z, X, δ) = 0 in the sense of Section 6.1.1, where δ > 0 is a standard step size.Lemma 66 Lemma 13 extends to nonstandard systems of the form (145).Proof Select z = [z n ] and x = [x n ] satisfying H ( z, x, ∂) = 0. By definition, we have [H (z n , x n , δ n )] ∼ [0, 0, . . . ],where relation ∼ was defined in (110). Equivalently,

	(145)

  • 2 and the current torques τ 1 , τ 2 from the current values for the states ω 1 , ω 2 . The occurrence of ∂ in the right-hand sides of equations (e ∂ 1 ), (e ∂ 2 ) is the difficulty, since ∂ is nonstandardsee Comment 4 regarding this.

  , ψ preserves the Dulmage-Mendelsohn decompositions of G F and G F a , meaning that g ∈ β a O/E/U if and only if ψ f (g) ∈ β O/E/U , and y ∈ β a O/E/U if and only if ψ x (y) ∈ β O/E/U .

	Proof Statement 1 is a direct consequence of Lemma 79 and (174). We thus focus on Statement 2. By
	Statement 1, we can see that

  •(M +1) x or x •(M +1) if expansion (165) is used. Removing Assumption 9 could be handled similarly, by redefining transient modes as modes of nonstandard duration M ∂ if expansion (165) is used, instead of ∂ if expansion (164) is used.

RR n°9334

By compilation, we mean here the suite of symbolic analyses and transformations that must be performed prior to generating simulation code.RR n°9334

https://www.wolfram.com/mathematica/new-in-9/advanced-hybrid-and-differential-algebraic-equations/ Inria

The left-limit is available in Simulink under the 'state port' construct, whereas the pre keyword in Modelica serves the same purpose.RR n°9334

https://www.modelica.org/publications RR n°9334

There are indeed a great number of possible definitions for "the index" of a DAE, serving various purposes, see[START_REF] Campbell | The index of general nonlinear DAEs[END_REF] for instance. This, however, is not relevant to the present work. Inria

The skeptical reader is gently referred to Section 9 for mathematical details.

The result of this lifting may not be intuitive in general, see Section 9 for mathematical details.

In French, one uses to qualify this as: avoir le beurre et l'argent du beurre.RR n°9334

Throughout this paper, we consistently use, in acronyms, the letters "D" and "d" to refer to "Differential" and "difference", respectively. Inria

This informal reasoning will be formalized in Section 10.2, in the context of nonstandard analysis. Inria

The Cup-and-Ball exampleIn this section, we study the example of Section 1.2.2 and address the objectives stated at the end of that section.RR n°9334

Generically means: outside a set of values of Lebesgue measure zero. Inria

Lemma 20 only states that there cannot be two arbitrarily close solutions of F = 0, as its conclusionS ∩ V = {x * } involves an open neighborhood V of x * . RR n°9334

Examples are linear systems for which elimination is easy, and polynomial systems for which it is doable-but expensive-using Gröbner bases.

The precise meaning of this statement is the same as in Lemma 13, see the discussion therafter. Inria

The notation x k is adopted throughout this paper, instead of the more classical x (k) , for the k-th derivative of x. See Notations 26.RR n°9334

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. Inria

This is unlike Synchonous Languages for discrete time systems, in which the restriction resulting from adopting Definition 33 would be significant.

With this new syntax, the original forms for the Cup-and-Ball and Westinghouse examples (given in the introduction) could not be expressed.RR n°9334

The prefix "S" is reminiscent of constructive Semantics.RR n°9334

Assertions may be used to restrict the set of possible modes, by stating some physical knowledge the programmer may have, e.g., 'this sequence of modes is not possible'. This is a minor change to our approach, not developed here.RR n°9334

This means that we regard x , x , etc. as different variable names and forget about the fact that they are related via differentiation. The same term holds whith shifting. Inria

Corresponding to the Simulink class of modeling languages, with no DAE.

Observe that, as a consequence, µ cannot be sigma-additive (in contrast to probability measures or Radon measures) in that it is not true that µ( n An) = n µ(An) holds for an infinite denumerable sequence An of pairwise disjoint subsets of N.

More generally, the Transfer Principle states that every first-order formula is true over X if and only if it is true over X. Recall that a first-order formula is a formula involving quantifiers over variables but not functions. Inria

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

Acknowledgements

The authors wish to thank a number of people for this work. First, Hilding Elmqvist and Martin Otter introduced them to the subject in the early 2010's; our collaboration then established for a couple of years, leading to both a first version of our approach and the handling of semi-linear mDAE systems; since then, Hilding and Martin have focused on their grand child of Modelica called Modia. We had extensive fruitful discussions with John Pryce on both structural analysis of DAEs, including the Σ-method}, and physical modeling in general; this helped us bring out the need for an explicit distinction between long and transient modes. Khalil Ghorbal participated to the first version of this approach and further interacted with the authors for this revised approach. Finally, Vincent Acary introduced us to nonsmooth systems techniques and solvers and we had frequent exchanges regarding the notions of solution for mDAE systems, revealing to us how strange these can be.

have µ-measure 0. Thus, we say that

respectively. Indeed, the same trick works for many other relations and operations on nonstandard real numbers, as we shall see. We now proceed with a more formal presentation.

Nonstandard domains

For I an arbitrary set, a filter F over I is a family of subsets of I such that: (1) the empty set does not belong to F; (2) P, Q ∈ F implies P ∩ Q ∈ F; and (3) P ∈ F and P ⊂ Q ⊆ I implies Q ∈ F. Consequently, F cannot contain both a set P and its complement P c . A filter that contains one of the two for any subset P ⊆ I is called an ultra-filter. At this point we recall Zorn's lemma, known to be equivalent to the axiom of choice:

Lemma 57 (Zorn's lemma) Any partially ordered set (X, ≤) such that any chain in X possesses an upper bound has a maximal element.

A filter F over I is an ultra-filter if and only if it is maximal with respect to set inclusion. By Zorn's lemma, any filter F over I can be extended to an ultra-filter over I. Now, if I is infinite, the family of sets F = {P ⊆ I | P c is finite} is a free filter, meaning it contains no finite set. It can thus be extended to a free ultra-filter over I. Hence: Lemma 58 Any infinite set has a free ultra-filter.

Every free ultra-filter F over I uniquely defines, by setting µ(P ) = 1 if P ∈ F and 0 otherwise, a finitely additive measure 24 µ : 2 I → {0, 1}, which satisfies µ(I) = 1 and, if P is finite, then µ(P ) = 0.

(109)

Now, fix an infinite set I and a finitely additive measure µ over I as above. Let X be a set and consider the Cartesian product X I = (x i ) i∈I . Define

if and only if µ{i ∈ I | x i = x i } = 0. Relation ∼ is an equivalence relation whose equivalence classes are denoted by [x i ] and we define

X is naturally embedded into X by mapping every x ∈ X to the constant tuple such that x i = x for every i ∈ I; we denote it by [x] or, simply, by x, if no confusion can result. Any algebraic structure over X (group, ring, field) carries over to X by almost pointwise extension. In particular, if [x i ] = 0, meaning that µ{i | x i = 0} = 0 we can define its inverse [x i ] -1 by taking

The existence of an inverse for any non-zero element of a ring is indeed stated by the formula: ∀x (x = 0 ∨ ∃y (xy = 1)). 25 where Z = def (X • , V ) collects the dependent variables. Such systems are also written as

after expanding X as X • -X

∂

. In System (124), V collects the algebraic variables, X collects the state variables, and X • -X ∂ is the nonstandard semantics of X . Writing (124) involves both X and X • but not ∂. As (125) is obtained from (124) by performing the expansion, it involves both X • and ∂ but not X .

As an illustration, with reference to the Cup-and-Ball example, System (30) is of the form (124)-we repeat it here for convenience:

Its counterpart following form (125) is obtained by expanding the second derivatives by using Euler scheme:

Comment 63 System (124) is the generic form of systems we need to solve at Line 10 of Algorithm 13, when we eliminate higher order derivatives and shifts by introducing auxiliary variables.

We consider the following assumptions, which hold for the nonstandard semantics of any standard mDAE system: Assumption 6

1. The functions H(•) and H(•) are standard (they result from applying Algorithm 13 to our given standard mDAE model);

2. The state variables X are standard and finite (their values were assigned by the standardized dynamics of the previous mode).

Impulse Analysis

We refer the reader to Section 2.5 about the standardization of the clutch model at mode changes, and particularly System [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF] in it. As a prerequisite to standardization, we had to identify which variables could become impulsive at mode changes (the two torques in this example), and then, we had to eliminate them.

Inria

Proof Recall that the time basis of nonstandard semantics is T = {k×∂ | k∈ N}. Denoting by t k the generic element of T, transition system (149) rewrites:

Using the Implicit Function Theorem, let G be the standard function such that X = G(X) solves (148) for X , and set G = [G, G, G, . . . ]. Then, (150) rewrites

Theorem 69 follows from Theorem 61 applied to (151).

In the next two sections we investigate the standardization of the blocks returned by Algorithm 13, as identified in Section 7.3.5. Since, in this section, we pay attention to the status standard/nonstandard, we will make this status explicit in our development whenever needed.

Long modes of mDAE

We use the notations and results of Section 7.3.5. More precisely, we reuse notations (103) with replacing •, and we invoke Lemma 42. We consider a system B (H, N ) = (B t (H, N ), B h (H, N )) returned by Algorithm 13 for long modes. By construction of Algorithm 13, the algebraic system G Σ = 0, where G Σ = def B t (H, N ), is structurally nonsingular and admits ḠΣ = 0, where ḠΣ = def B h (H, N ), as associated consistency conditions.

At this point, we need to relate the above blocks to their standard continuous-time counterparts. The reduced index systems (G Σ , ḠΣ ) returned by Algorithm 8 (IndexReduc) in case of success, coincide with the two systems returned by the Σ-method, when replacing differentiation by shifting (see the emphasized sentence when introducing this algorithm). In our study of standardization, we will also need to consider the original outcome of the Σ-method (using differentiation); we denote this pair by (G c Σ , Ḡc Σ ), where the superscript c refers to continuous time. Continuous-time systems (G c Σ , Ḡc Σ ) define the leading equations and consistency conditions of a (standard) DAE system that is generically regular. In the following theorem, we denote by X = (X(t)) t∈R + the solution of this DAE system.

Since we are concerned with standardization, we carefully distinguish between H, a standard vector function, and H = [H, H, . . . ], its nonstandard lifting.

Theorem 70 The dAE system B t

• ( H, N )=0, with consistency conditions B h • ( H, N )=0, possesses a unique solution. Its piecewise constant interpolation, denoted by X(t), t∈ R + , satisfies X = st( X). Note that this standardization is independent from the particular infinitesimal time step ∂.

Theorem 70 formalizes the practical informal rule that for long modes, the standardization of

where

whereas B t (H, N ) and B h (H, N ) are constructed from (153) as in Notations 41.

Inria

Assuming again that the f i 's are linear in τ , i.e., of the form

We can zero the term in ∂ in the first equation: doing so keeps the system structurally nonsingular, hence it is safe. It then remains to standardize

and, finally, solve it for the dependent variable ω

. For n = P, . . . , N , ω •n i standardizes as the left-limit ω - i for i = 1, 2, whereas ω •(N +1) standardizes as the right-limit ω + . Hence, (163) standardizes as

which does not depend on α N since α N = 0. Finally, we get

, which coincides with [START_REF] Elmqvist | Systems modeling and programming in a unified environment based on Julia[END_REF].

A general result

In this section, we consider again the generic nonstandard interpretation (158) for the derivative. We investigate under which assumptions the result of our whole compilation suite is independent from the particular way derivatives are expanded when mapping the considered mDAE system to its nonstandard semantics. We consider the following assumption on S, thus characterizing a restricted class of ∂AE:

Assumption 9 mDAE-system S involves only long modes, no transient modes.

Let us compare, under Assumption 9 and Comment 63, the results of Algorithm 10 on S for the following two nonstandard semantics for the derivative:

and

where a = (α P , . . . , α N ) satisfies P ≤N , N ≥0, α N =0, and N n=P α n = 1, whereas α P -1 =0 by convention. Note that nothing prevents P from being negative. For the two formulas, we have given both the direct and inverse forms. Equations (164) are obtained by taking P =N =0 and α N =1 in (165).

Comparing the structural analyses

Let S be an mDAE and let H S be the system H returned by Algorithm 10 at Line 6, i.e., the outcome of index reduction, for both a long mode and a mode change event.

We will consider the following problem: how should H S be modified to become the outcome of the index reduction of S, when derivatives are expanded using (165) instead of the first-order Euler expansion (164)? We denote by S a and S the corresponding two nonstandard semantics for S. Regarding this question, two lines of Algorithm 10 are important: Line 3 and Line 6. We investigate them successively.

Line 3 of Algorithm 10 This consists in applying the Σ-method to the original mDAE system S, within a given mode; let G=0 be the dynamics of S in this mode. The outcome of this algorithm only depends on the weighted bipartite graph G associated to G. In the original Σ-method, the weight of an edge (f, x) of graph G is equal to the highest differentiation degree of x in f . Similarly, if G is a dAE system (in discrete time), the weight of an edge (f, x) of G is equal to the highest shifting degree of x in f . Having recalled this, we can now compare the two expansions:

• Using expansion (164):

1. We begin with DAE system G and its bipartite graph G;

2. Replacing DAE system G by dAE system G using expansion (164) yields G, that is identical to G.

• Using expansion (165):

1. We begin with DAE system G and its bipartite graph G;

2. Replacing DAE system G by dAE system G using expansion (165) yields G a , in which the vertices and edges are exactly those in G and all weights have been multiplied by

The case N = 0 is the base case of the forward first-order Euler scheme, so that we do not need to consider it. Therefore, in the sequel, we assume N > 0, i.e., M > 1.

Using the notations of (65), let (c f ) f ∈F and (c a f ) f ∈F be the offsets associated to G and G a . Then, by Lemma 25, one gets: Lemma 73 The Σ-method succeeds on G if and only if it succeeds on G a , and, in this case, c a f = M × c f for every f ∈F and d a x = M × d x for every x∈X(F ).

In particular, each guarded equation in S a corresponds to an equation in S shifted M times. We denote by

the corresponding two returns of the Σ-method.

Line 6 of Algorithm 10 If the current instant sits within a long mode, then this line just erases the consistency conditions G Σ and(G) a Σ . So, we are left with G Σ if we use expansion (164), and G a Σ if we use expansion (165). Now, dividing the infinitesimal step ∂ by M in G Σ yields G a Σ . Theorem 70 tells us that the standardizations of the two systems G Σ and G a Σ yield the same DAE system. As such, the difficulty is the handling of mode changes by Line 6, when the two different expansions are used. Is the handling of conflicts preserved, when the two different expansions are used? More precisely, for t o ∈ T an instant of mode change:

Inria

The reasoning performed at Section 7.2 applies to semilinear mDAE systems. Following (77), we unfold (182) at the two successive instants • t * and t * :

Equations (e -) and (e + 2 ) are in conflict: we discard the latter. This finally yields

In System ( 184), (e -) is a context, whereas (e + 1 ) and (e + 3 ) determine X • as a function of X. Setting ∂ ← 0 in (e + 1 , e + 3 ) yields a structurally nonsingular system by Condition 3 of Definition 84. Hence, by Theorem 67, setting ∂ ← 0 in (e + 1 , e + 3 ) performs the standardization: the scheme of Theorem 87 is recovered.

Theorem 88 (correctness result) For semilinear mDAE systems, our approach computes the correct solution.

Comment 89 Our approach computes solutions for any mDAE system whose structural analysis (Algorithm 10 or Algorithm 13) succeeds. This is way beyond the class of semilinear mDAE systems. It would be desirable to compare the so obtained scheme with schemes known for more dedicated physics.

Comment 90 Hilding Elmqvist and Martin Otter have developed the ModiaMath tool for semilinear mDAE systems by effectively implementing scheme (181).

Numerical scheme for restarts

In this section, we focus on the system of equations defining restarts, in its form (125), which we recall for convenience:

Our aim is to justify our claim that, if we remain 'blind' with respect to impulsive variables, the solution of this system is well approximated by the standard system

where the infinitesimal ∂ has been replaced by a small standard parameter δ > 0. This is essentially the proposed numerical scheme. Throughout this section, Assumptions 6 and 8 of Section 10.5 are in force.

Inria γ : f → t. First, we perform impulse analysis following Section 10.1. Based on this impulse analysis, we apply a change of variable that bypasses the need for eliminating impulsive variables. We then standardize the resulting algebraic system of equations and we prove that it possesses a unique solution.

Using the structural analysis method detailed in Section 2.4, the state-jump is solution of the following system of algebraic equations, where ω • 1 , ω • 2 , τ 1 and τ 2 are the unknowns:

Applying the method of Section 10.1, the impulse orders of the dependent variables ω

τ 2 are solution of the following system of equalities and inequalities:

. This impulse analysis tells us that if a solution to System (189) exists, then unknowns τ 1 and τ 2 are of order O(∂ -1

3 ). The following substitution helps solving System (189):

With this substitution, System (189) becomes:

Let us replace each unknown by its unique decomposition into a standard part and an infinitesimal part. Regarding velocities, we have:

, where φ i ≈ 0, for i = 1, 2. Next, since the rescaled variables τi are non-impulsive, the same substitution is legitimate by Lemma 59: τi = τi + ψ i , with ψ i ≈ 0, for i = 1, 2. System (191) then becomes:

Zeroing the infinitesimals in these equations yields the following nonsingular system of equations, which, by Theorem 67, is then the standardization of System (192):

This standard system of equations admits a unique solution: properly handles any multimode system having at most finite cascades of mode changes-sliding modes are not supported by our method.

A key step was the extension of structural analysis from single-mode to multimode systems. In this extension, the handling of mode changes was the main difficulty. We illustrated our approach on small examples that are not properly supported by existing tools (with the exception of ModiaMath). These examples exhibit mode-dependent index, which is not well supported by existing techniques.

For the clutch example, the most natural model consists in specifying the dynamics for the two modes "released" and "engaged"; then, the compiler automatically generates the restart conditions at mode changes. The "Westinghouse air brake" example motivated the consideration of assertions-not studied here but included in our tool under development. The Cup-and-Ball example is interesting in that, when the rope gets straight, inelastic or elastic impact can be hypothesized. This example pointed the issue of long vs. transient modes (modes having zero duration). Depending on whether the straight rope mode was assumed long (inelastic impact) or transient (elastic impact), the original model was accepted or rejected. For the latter, underspecification was detected and this information was returned to the user, who could refine his/her model by specifying an impact law. The so enriched model was accepted and code properly generated at mode change.

Generally, our method allows us to reject models on the basis of under/over-specification. Logiconumerical fixpoint equations (cases in which numerical variables determine the value of a Boolean condition that, in turn, guards the equation allowing to evaluate the same numerical variables) are prohibited by setting restrictions on the language syntax.

The structure of the code we generate is reminiscent of the so-called dummy derivatives method [START_REF] Mattsson | Index reduction in Differential-Algebraic Equations using dummy derivatives[END_REF]: for both the continuous dynamics and the handling of transient modes or mode change events, we call an algebraic equation system solver for computing derivatives or next state values, as a function of the current state, while meeting algebraic constraints.

A first requirement on our theory was that it should support the analysis of any model related to any kind of physics-we aim at being physics agnostic. This support should hold regardless of the particular properties of the model, such as conditions ensuring existence and uniqueness of solutions. The motivation for this was that the user is not expected to check such conditions before Inria • either b = f (Line 7), in which case ExecRun returns Fail (σ) and stops, or • (Line 8) a pair consisting of G Σ , a structurally nonsingular system determining some leading variables, and the system G Σ collecting consistency equations. In this case, the algorithm can further progress.

Line 9 With the additional help of the set Φ of facts, the possible conflict is solved using the SolveConflict atomic action (Algorithm 9).

Line 10 SolveConflict fails if no structurally regular system survives as a result of applying Line 9.

In this case, ExecRun returns Fail (σ) and stops.

Line 11 Otherwise, SolveConflict returns a structurally regular system H, which is solved using Eval (82), which updates the status σ by giving a value to more S-variables. As a result, the set Φ of facts is updated.