{'fraction_non_alphanumeric': 0.04588663906943638, 'fraction_numerical': 0.01400621874239047, 'mean_word_length': 4.193438556385126, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 99, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Mathematical Foundations of Physical Systems Modeling Languages

Modern modeling languages for general physical systems, such as Modelica, Amesim, or Simscape, rely on Differential Algebraic Equations (DAE), i.e., constraints of the form f(dot{x},x,u)=0. This drastically facilitates modeling from first principles of the physics and 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, DAE exhibit subtle issues because of the notion of differentiation index and related latent equations -- ODE are DAE 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 \emph{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 mode-dependent DAE based dynamics and state-dependent mode switching. Multimode DAE are much more difficult than DAE. The main difficulty is the handling 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 well sound criteria for accepting or rejecting models. For our mathematical development, we rely on nonstandard analysis, which allows us to cast hybrid systems dynamics to discrete time dynamics with infinitesimal step size, thus providing a uniform framework for handling both continuous dynamics and mode change events.


Overall motivations
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 [8]. 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 [8], 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 [2,45]; 2. Developing existence/uniqueness results and discretization schemes, for a given class of (single-mode or multimode) DAE systems [3]; 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 (1). 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: (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 and 4.
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 [40], 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
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 [1], 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: (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 ) (5) 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 (5). 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 (5) would consist in replacing, in (k 0 ), s by its left-limit s − (t) = def lim sup u t s(u). 3 Our 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  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 [1], 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: The system of Figure 8, where (ν 1 , ν 2 , ν 3 ) is replaced by (v 0 , . . . , v 4 ) in (6), has the form (1). 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 RR n°9334 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., [16] 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 [44,43] 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 [50], 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 [30]. However, it is not obvious how to generalize methods in this restricted area to general multimode DAEs.
The article by Mehrmann et al. [38] 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 [54], 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. [41]. Dynamically changing the structural analysis at runtime is also performed by Hoeger [31,32]. In [31], the author proposes a dynamic execution of John Pryce's Σ-method [48].
Elmqvist et al. [28,36] 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 [1] 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. [7,6] 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 [7], by significantly modifying the original Pantelides algorithm [40]. 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 [52] and his article [53], 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 [34] 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 [8]. 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 [8] 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 [8] 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 [8] when applied to the subclass of semi-linear systems, see Section 11.2. Our present work thus extends and significantly improves [8].

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 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 [49,23]. Impulsive values become normal citizens in nonstandard analysis. Continuous time and dynamics can be discretized using infinitesimal time steps, e.g., by setting 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 [8]). 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 [19]. 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 (3). 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 : In the engaged mode, however, γ holds true, and both the two velocities and the two torques are algebraically related: 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 (8), 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 (8). 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" [48].
It turns out that this condition does not hold for System (8) 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: 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 (9), 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 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.
update the states (ω 1 , ω 2 ) 3: Tick move to next discrete step Since the new values of the state variables satisfy (10) 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: 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.
The replacement of (e 3 ) by (e 3 ), which resulted in Exec. Sch. 2, is known in the literature as index reduction [37]. 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 [21], 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.
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 (3).
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 [49,35,5] 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 [35]: 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).
We can then consider, for a given positive infinitesimal ∂, the following time index set T ⊆ R: 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.
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 . 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 Definition 3 Substituting, in a multi-mode DAE system, every occurrence of x (t) by its expansion (15) yields a difference algebraic equation (dAE) system 9 that we call its nonstandard semantics.
For instance, the nonstandard semantics of System (3) is the following: 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 (16) 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 (16). 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): System (17) 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 • γ = 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:  (18) 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: The above reasoning may seem an overshoot for this mode change γ : t → f, since one could simply say: let the previous part of System (18) 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: Once more, we regard System (20) as an algebraic system of equations with dependent variables , the leading variables of System (17) at the previous and current instants.

RR n°9334
In contrast to the previous mode change, System (20) is no longer structurally nonsingular since the following subsystem possesses five equations and only four dependent variables ω 1 , ω 2 , • τ 1 , • τ 2 : 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 (21), the conflicting equation (e 3 ). This yields the following code for the mode change γ : f → t: Note that the consistency equation (e 3 ) : ω 1 − ω 2 = 0 has been removed from System (22), 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.

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 (13) 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 (7) in the released mode γ = f.
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.
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 ∂ 1 , e ∂ 2 , e • 3 , e 4 } to determine the leading variables (τ 1 , τ 2 , ω • 1 , ω • 2 ). In particular, from e ∂ i , we get 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 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: 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: 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.
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 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. 10 When A is singular (as in our case), the solution consists in eliminating the impulsive variables from System (25), namely, the two torques. This is performed by symbolic pivoting: 1. Using (e 4 ) yields −τ 2 = τ 1 = def τ ; 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.

Premultiplying the system of equations
by the row matrix b 2 (ω 2 ) b 1 (ω 1 ) yields the equation

Using in addition (e • 3 ) and setting
The symbolic-only rewriting used here actually alleviates any difficulty in standardizing equation (26): 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 : 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 (25). 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).

The Cup-and-Ball example
In this section, we study the example of Section 1.2.2 and address the objectives stated at the end of that section.

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 (5).
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 ∂: In the sequel, when referring to model (5), 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 (28).

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): 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 (29), where the modification is highlighted in red: This model is understood in the nonstandard semantics, meaning that the derivatives are expanded using (28). 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: 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: Applying again Principle (1) of causality, we erase, in model (30), equations (k 1 ) and (k • 1 ) at the instant of mode change • γ=f, γ=t. This yields the following system: 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: 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 (31) and (32). 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: In System (33), we expand second derivatives using (28). Note that substitution (28) implies that we can as well expand the second derivatives as functions of first derivatives: Since we will use System (34) to evaluate restart values for the velocities, we will use expansion (35) 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): Regarding (36), 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: 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 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 (35): By expanding x •2 = x • + ∂x • , the right-hand side of the last equation rewrites Using this expansion of L 2 −(x 2 +y 2 ) •2 , setting ∂ = 0 in (39) yields 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 (39). 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: and the restart velocity is evaluated by After similar substitutions for y, we get 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 (29). 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 (29). 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 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 [21], 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: 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 where 0 ≤ α < 1 is a fixed damping coefficient. This is reinterpreted in the nonstandard domain as yielding the following refined system for use at mode γ=t within the cascade γ:f→t→f: 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 (47), which leaves us with the following system: 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: 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 [15]. 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. 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 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 , 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. 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 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. 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 for which a complete matching is shown in red in Figure 15. This system is structurally nonsingular. Figure 12: The active system when γ = f. We show in red a structurally singular subsystem. 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 (28). 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 Figure 14: The active system when γ = t. We show in red a structurally singular subsystem. Compare with Figure 12. 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 γ : 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.
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 (28) 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. Figure 16: Restart at mode change γ : t → f.
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.
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 Figure 17: Restart at mode change γ : f → t.
mode change.

RR n°9334
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.
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 Figure 11 become identical due to the zero value for s. As a consequence, whatever next mode is, the consistency equations for it ( 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 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 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 [19,20], 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.
• 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.
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: 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 [24]) 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 [11]. 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 (53); 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 [40,37,46,48]): Lemma 13 Assume that F is C 1 . The following properties are equivalent: 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 11 Generically means: outside a set of values of Lebesgue measure zero.

RR n°9334
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 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. 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).
x 1 x 2 x 3 x 6 x 4 x 5 This lemma is illustrated in Figure 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 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 (53)), 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 [46] 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: 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., [22], 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: 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 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: 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 , 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). (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.

Definition 21 System
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: 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 partition β E = F Σ ∪ F Σ 7: 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.
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 [40] is the historical solution. We decided, however, to base our subsequent developments on the beautiful method proposed in [48] 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: 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. 16 Definition 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 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. primal : f :(f,x)∈G ξ fx = 1 and x:(f,x)∈G ξ fx ≥ 1 and ξ fx ≥ 0 (61) dual : 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 which is exactly the last constraint of (60). Using this translation into linear programs (61) and (62), it is proved in [48], 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 [48] 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}: 2. Apply the following iteration until a fixpoint is reached (in finitely many steps), from the initial values c f = 0: 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 where M is a fixed positive integer. Then, the offsets of the corresponding Σ-method are also related in the same way: 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 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 Having b = t returned by ExistsMatching() is a prerequisite to calling FindOffsets(), which is then guaranteed successful.
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.

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 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. For mdAE systems, the formalization of the second statement is straightforward: equation f =0 is enabled when For mDAE systems: equation f =0 is enabled when where γ − denotes the left-limsup of γ at the current instant, defined by: The left-limsup is always defined and coincides with the left limit when the latter exists. 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. where b is a long/transient-typed Boolean expression of predicates over X ( ) (respectively X (•) ).

Long versus transient modes
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 µ.

RR n°9334
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: 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. 17 Therefore, 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: 18 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 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 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): 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 ∂: In addition, we replace, in (76), the derivative d . 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: 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 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 [12,9,10,47], 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:

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,  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 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 (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: We denote by i u, u f, and u t the increasing changes of status for a given S-variable. 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 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.
1: Tick(σ) = def (σ • , ∆ • ), where: 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
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). (β U , β E , ∅) ← DM(K) with dependent variables set by σ 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
if ¬b then return Fail (σ) σ ← Eval(σ, H); increase σ 10: 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 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 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 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  Figure  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.
Detecting continuous modes 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.
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.
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 (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 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 Σ , 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 ).

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 [21]) 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. (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.

Definition 43 A time-varying dAE system over T is a pair
Example 44 Consider T = N, X = (u, v, w), and F is defined by t=0 : 0 = v • − 3u ; and ∀t > 0 : The solutions of this system are trajectories t → (u(t), v(t), w(t)) such that .
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: 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) In its Section 3.1 on 'standard indices', Campbell and Gear's landmark paper [21] 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: 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 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.

4:
A k ← (F 0 , . . . , F k ) 5: arrange the variables of A k as (Y k ∪X k ∪W k ) 6: if ¬b U then go to Line 2 10: else return (b ← t, F Σ , F Σ ); increase σ 11: else consider A defined in (88) 12: 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. 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
if ¬b then return Fail (σ) 9: else 10: σ ← Eval(σ, H); increase σ 11: 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 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. This set can be identified using structural algorithms.

Example 47
The statement expands, according to Definition 27, as whose nonstandard semantics is Since x • −x is not guaranteed infinitesimal by (e 1 ), x • −x is non-infinitesimal at the instant of zerocrossing (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 (46). 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 (46).
In the sequel, we denote by M the set of all possible modes of system S. 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.  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: The resulting formalism will be called mdDAE, or ∂AE, for short. We formalize this next.

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 where ∇f denotes the Jacobian of f . Then, differentiation and shift commute, meaning that: Proof Statement 1 follows immediately from identities (94). Then, the equalities 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 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.

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

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 . . . 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 and the definitions of B h (F, W ) and B t (F, W ) follow accordingly. In addition, we will write 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 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: (b) case ν = : define h j := h j−1 and 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 where ∇ is the Jacobian of .
Lemma 55 Applying map (104) to all the entries of B h (f, w), B h (f, w), and B t (f, w), returns respectively, where n = |w|. Doing the same for B(F, W ) defined in (103) returns where N is the vector of integers collecting the lengths of the entries of W . 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 [49,23]. 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 [35]. 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. [33] first proposed using nonstandard analysis to discuss the nature of time in hybrid systems. Bliudze and Krob [14,13] 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. [5] used extensively the nonstandard semantics of hybrid systems, with practical consequences for how to design a specification formalism for multimode ODE systems [15,4]. 22 See also [51,39,29] for similar studies using the axiomatic approach to nonstandard analysis. The following presentation is borrowed from [5].

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 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 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.
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 y i = x −1 i if x i = 0 and y i = 0 otherwise. This construction yields µ{i | y i x i = 1} = 1, whence [y i ][x i ] = 1 in X. 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

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 [u n ] ≤ [v n ] iff µ{n | u n >v n } = 0. 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 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 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.
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: 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 : 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 One can rewrite the positive infinitesimal quantity 1/|T | as ∂, so that 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, By (121), the following holds: 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 [35] 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 [35]. 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
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 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 .
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 (25) 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.
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
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.   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 (25): rewriting this equation as a 1 (ω 1 ) =

The rules of impulse analysis
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

RR n°9334
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: 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 and 22 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

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: The dependent variables are X • , V . It will be convenient to introduce the auxiliary variables so that the systems we consider take the following form: 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.
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 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 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 : prior information on impulse orders : 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 and 22, 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 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. 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: We then saturate System (134) using Rule (R8) of Figure 22 (added equations are in blue): RR n°9334 and we consider the prior information (132): 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: We then saturate System (138) using Rule (R8) of Figure 22 (added equations are in blue): We also consider the prior information (132): We now apply the rules of Figure 22 to System (139), then use (140):

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 (25), which we recall for convenience: 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 (25). The details were given in Section 2.5.2 and the resulting reduced system of equations is 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: 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 (25), replace in it the infinitesimal time step ∂ by a small real positive time step δ > 0, and expand the derivatives as before. This yields 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. 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 (54) 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 satisfying property (54) 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.
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).
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.

DAE of index zero
Invoking the implicit function theorem, one can also address DAEs of index 0, i.e., of the form: 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: 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 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 (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 B t where whereas B t (H, N ) and B h (H, N ) are constructed from (153) as in Notations 41.

Mode changes of mDAE
Assumption 6 (stated on page 74) 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: 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 (16).
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 (25), used to determine the restart values ω • 1 , ω • 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.
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: 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 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: as explained above; 2. The considered model is a reduction of the original one, in that, equivalently: when γ do τ 1 = NaN; τ 2 = NaN; ; 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: Let us investigate how the analysis of the clutch is modified when using the nonstandard interpretation for the derivative. The nonstandard semantics of the clutch model is: 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 ) 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: ) 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 Assuming again that the f i 's are linear in τ , i.e., of the form We can eliminate τ : 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 ω •(N +1) = def ω . 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 (27).

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).   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).
Step 2: Removing, from the set of dependent variables of G Σ and G a Σ , the variables involved in the contexts ∆ and ∆ a Let 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, 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 RR n°9334  and define  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 74), 8 (Page 84) and 9 (Page 87) 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: 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: 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:

RR n°9334
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:

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 [8] 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 [26]. As a consequence, if f has bounded variation, then Definition 84 (semilinear systems) Call semilinear an mDAE system S such that, for each mode µ, the active DAE system takes the restricted form 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.

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

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 [8], 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: 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 s (X s (u))X s (u)+B + (X(u)) du A n (X s (u))X n (u)du 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 ( [8]) The restart conditions X + are determined from the left-limits X − by using the following system of equations: 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: and applying the Σ-method to it yields the augmented system 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 * : Eliminating • X in (183) yields 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. 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 that remains structurally nonsingular when δ = 0.
Proof We successively prove the two statements.
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: (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 ) 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 γ : 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 ω • 1 − ω 1 , ω • 2 − ω 2 , τ 1 , τ 2 are solution of the following system of equalities and inequalities: Let us replace each unknown by its unique decomposition into a standard part and an infinitesimal part. Regarding velocities, we have: ω • i =ω • i +φ i , 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): RR n°9334 This standard system of equations admits a unique solution: 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 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 2: repeat 3: y ← x 4: u ← NSolve(system = F h , guess = u, tolerance = h) 5: 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 solutionx = (ω 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 .

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, 'onthe-fly' techniques for the structural analysis of such systems, as advocated in [17] and [27], 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 [18]) for implicit representations of the mode-dependent structure of an mDAE. This is described more thoroughly in [19,20]. 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: To avoid a logico-numerical fixpoint, we consider instead the system obtained by replacing s with its left-limit s − in equation (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.

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. 0 = i 1 + i 2 + j 1 + j 2 (K 1 )

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 and 34, 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 [19,20], 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 [19,20].

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 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 Modia-Math). 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. Logico-numerical 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   The structure of the code we generate is reminiscent of the so-called dummy derivatives method [37]: 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 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  Figure 35: 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". 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 [12,9,10]. However, a particular difficulty holds: no mathematical definition of solutions of an mDAE system exists in full generality, 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 [19], where our IsamDAE tool under development is described.
• 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.