Albert Benveniste

Benoît Caillaud

Mathias Malandain

From Hybrid Automata to DAE-based modeling

Keywords: Hybrid automata, multiphysics modeling, DAE, multimode DAE, structural analysis, index reduction

Tom Henzinger was among the co-founders of the paradigm of hybrid automata in 1992. Hybrid automata possess different locations, holding different ODE-based dynamics; exit conditions from a location trigger transitions, resulting in starting conditions for the next location. A large research activity was developed in the formal verification of hybrid automata; this paradigm still grounds popular commercial tools such as Stateflow for Simulink. However, modeling from first principles of physics requires a different approach: balance equations and conservation laws play a central role, and elementary physical components come with no prespecified input/output profile. All of this leads to grounding physical modeling on DAEs (Differential Algebraic Equations, of the form f (x ′ , x, v) = 0) instead of ODEs. DAE-based modeling, implemented for example in the Modelica language, allows for modularity and reuse of models. Unsurprisingly, DAE-based hybrid systems (also known as multimode DAE systems) emerge as the central paradigm in multiphysics modeling. Despite the growing popularity of modeling tools based on this paradigm, fundamental problems remain in the handling of multiple modes and mode changes-corresponding to multiple locations and transitions in hybrid automata. Deep symbolic analyses (grouped under the term "structural analysis" in the related community), grounded on solid foundations, are required to generate simulation code. This paper reviews the issues related to multimode DAE systems and proposes algorithms for their analysis. Computer science is instrumental in these works, with a lot to offer to the simulation scientific community.

Tom Henzinger and Hybrid Systems in computer science

According to the vision of Hybrid Systems by the computer science community [START_REF] Alur | Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems[END_REF], a hybrid system possesses control locations, characterized by invariants expressed as properties on continuous states; in each location, continuous-time dynamics are specified using Ordinary Differential Equations (ODEs); transitions between control locations are characterized by their pre-and post-conditions or, alternatively, by pre-conditions and a reset action. For the class of hybrid automata, a large body of knowledge, techniques and tools has been developed since the pioneering paper [START_REF] Alur | Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems[END_REF]. Tom Henzinger himself contributed extensively to this body of knowledge, with theory and tools [START_REF] Henzinger | HYTECH: The Cornell HYbrid TECHnology Tool[END_REF]. Interestingly enough, this vision still grounds commercial modeling tools, e.g., Stateflow, which highlights the depth of the 1990's vision of hybrid automata. Still, pushed by practical considerations, alternative approaches have emerged in the 1990's in the control engineering community, with ties going back to bond graphs [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF], a formalism that grounds the widely used industrial tool Simcenter Amesim.

2 Cyber-Physical Systems modeling: the need for DAEs As a running example, we consider the idealized clutch of Fig. 1, involving two rotating shafts with no motor or brake being connected. Each shaft possesses its individual dynamics, relating rotational velocity and torque. This clutch possesses two modes: the two shafts evolve freely when the clutch is released, and are coupled by a perfect contact when the clutch is engaged. At the mode change when the clutch gets engaged, a discontinuity occurs in the rotation velocities (they are generally different before engagement, and get equal in zero time as a result of the engagement). This causes an impulse in the torques. This part of the modeling task requires particular care, as we shall see.

Modeling from first principles of physics naturally leads to considering

acausal models [START_REF] Elmqvist | A structured model language for large continuous systems[END_REF] Modeling the clutch with hybrid state machines: Let us first consider the model of the clutch expressed using ODE-based hybrid state machines.

released location

ω ′ 1 = f 1 (ω 1 , τ 1 ) τ 1 = 0 ω ′ 2 = f 2 (ω 2 , τ 2 ) τ 2 = 0 γ:f→t restartF→T > < γ:t→f restartT→F ω ′ = f 12 (ω, τ ) τ = 0 engaged location (1)
State machine [START_REF] Alur | Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems[END_REF] possesses two locations: released (left) and engaged (right).

In the left-to-right transition, the label "γ : f → t" sitting above the transition indicates the event that triggers the transition; the label "restart f→t " sitting below the transition refers to the restart action resulting from the transition-its explicit formulation is given below. The ODE sitting in the released location is the model of two non-interacting shafts; for each one, the angular velocity ω i and torque τ i are related via an ODE, with the torque being zero when the clutch is released. The ODE sitting in the engaged location is clear as well: it is the model of a single shaft obtained by gluing the two shafts together; as such, the "engaged" model is global and its architecture does not reflect the physical architecture.

The two restart actions are now detailed. First, when the clutch gets released (t → f), the restart values for the two velocities are equal to the velocity of the global shaft just before mode change:

restart T→F : ω + 1 = ω + 2 = ω -. (2)
The symmetric restart action (f → t) is more involved. To derive it, we need to get closer to the physics, by specializing the two functions f i for i = 1, 2:

f i (ω i , τ i ) = 1 Ji (τ i -g i (ω i )) , (3)
where J i is the moment of inertia of shaft i, and g i (ω i ) is the friction. When the clutch gets engaged, the two shaft velocities are discontinuous, as they become equal in zero time. The resulting common value follows from the law of preservation of angular momentum, which writes:

(J 1 + J 2 )ω + = J 1 ω - 1 + J 2 ω - 2 .
Hence, the restart action when the clutch gets engaged is given by:

restart F→T : ω + = J 1 ω - 1 + J 2 ω - 2 J 1 + J 2 . (4)
Formulas ( 1)-( 4) fully specify the model of the clutch using ODE-based hybrid state machines. It can then easily be brought to the classical Alur-Henzinger hybrid automaton form [START_REF] Alur | Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems[END_REF]; this transformation has no impact on the issues discussed below.

DAE-based modeling: Generally, modeling in physics does not rely on hybrid automata. Let us consider electrical circuits as an illustration. First, the Kirchhoff laws are expressed as balance equations: the sum of currents at a node is zero; the sum of voltages along a loop is zero. Second, most components (e.g., resistors and capacitors) come with no prespecified input/output orientation.

A similar situation actually arises in mechanics or in thermodynamics. What happens if we follow a similar approach for the clutch? Getting the clutch model from first principles is straightforward and elegant, as shown in Fig. 1. The system and model architectures coincide: the overall model consists of a model of each individual shaft, completed with a model for the junction (the velocities coincide and the torques sum up to zero if the clutch is engaged, γ = t; the torques are zero if the clutch is released, γ = f). At an instant of mode change when the cluch gets engaged, the a priori different velocities of the two shafts will, in zero time, merge to a single identical velocity while the efforts on the shafts will climb up to impulses. The resulting common velocity will depend, in a determinate way, on the inertia and velocity of each shaft prior to engagement. Thus, the main characteristics of the model of Fig. 1 are the following:

1. The model has two modes: γ=f and γ=t yield different dynamics. 2. The model is acausal: In the released mode γ=f, the two torques are inputs and the two angular positions are outputs; in the engaged mode γ=t, both angular positions and torques are unknowns, with no input/output orientation being defined: the model is acausal in this mode.

ω1 ω2 τ1 τ2 J1 J2 shaft 1       ω ′ 1 = f1(ω1, τ1) (e1)       junction         if γ then ω1 -ω2 = 0 (e3) τ1 + τ2 = 0 (e4) else τ1 = 0 (e5) τ2 = 0 (e6)         shaft 2       ω ′ 2 = f2(ω2, τ2) (e2)      
Fig. 1. A simple clutch with two shafts and its model from first principles.

3. The mode transition is not explicitly specified: No equation directly addresses mode changes γ : t → f and γ : f → t. Still, restart conditions can be automatically inferred from the model of Fig. 1 at compile time, as we shall see.

Modularity in modeling:

What happens if we have a chain of clutches with N shafts? The DAE-based model of Fig. 1 just extends with more replicates. In contrast, the Hybrid Automaton model (1,2,4) must be globally rebuilt: the number of states is 2 N -1 ; the model in each state is modified if we glue one more shaft to the chain; the restart conditions are global and we must compute the equivalent moments of inertia for each chain of shafts with engaged clutches. The conclusion is clear: physical modeling with hybrid automata is not modular; multimode DAE-based modeling is stringently needed to achieve modularity.

Conclusion:

Acausal modeling makes the specification of large Cyber Physical Systems more modular, comfortable and elegant, and allows for better reusability. As a matter of fact, multimode DAE-based modeling is routinely used in focused areas such as electrical circuits [START_REF] Heemels | On the dynamic analysis of piecewise-linear networks[END_REF], contact mechanics [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Wiley[END_REF], or hydraulics. The modeling of large CPS, however, requires physics-agnostic modeling allowing for the combination of different physics together with the embedded control software. The Modelica language, which supports physics-agnostic multimode modeling since version 3.3 [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF][START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF], is a de facto standard in this area.

A new avenue of paradigms, issues, and questions

From the observations above, we identify new paradigms, issues and questions for the computer science community:

-How can we formally define the class of models that contains the clutch model shown in Fig. 1? -Is there room for a computer scientist's way of thinking in this new area? Issue 1 (acausal models) Rather than only considering ODEs (i.e., equations of the form y ′ =g(y, u) where y is the state and u is the input), acausal modeling relies on DAEs, i.e., equations of the form f (x ′ , x, v)=0 where both x and v are unknowns, with no specification of any input/output status. What are the consequences in terms of generating simulation code? Issue 1 was addressed by the DAE community by introducing the notion of index [START_REF] Campbell | The index of general nonlinear DAEs[END_REF][START_REF] Mattsson | Index reduction in differential-algebraic equations using dummy derivatives[END_REF]. An abstraction of the notion of index, in which regularity/singularity of the Jacobian is considered in a generic sense, was proposed as a way to better scale up [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF][START_REF] Pryce | A simple structural analysis method for DAEs[END_REF]. This solution, known as structural analysis of DAE systems, is extensively used in Modelica tools.

In addition to being naturally described by an acausal model, a system may have different modes for various reasons. Mechanical impacts where the bodies may remain in contact after a collision or an active grasping/docking, idealized electrical or hydraulic switching elements, are examples of situations requiring different sets of equations for their modeling. Another reason is the control of complex scenarios like repairing a satellite with robots, or switching between start-up, normal operation and shutdown during the operation of a power plant. Models of this kind are called multimode models. Hence the second issue for consideration:

Issue 2 (multiple modes) What is the consequence of considering multimode DAE models-in other words, DAE-based hybrid systems? Issue 2 was identified by authors, but is still incorrectly addressed in industrially used DAE-based modeling tools, as we shall see in the next Section. Authors have proposed dynamic approaches that perform structural analysis when needed [START_REF] Höger | Dynamic structural analysis for DAEs[END_REF][START_REF] Höger | Elaborate control: variable-structure modeling from an operational perspective[END_REF][START_REF] Elmqvist | Modiamath webpage[END_REF], at the cost of significant computational overheads at runtime. In addition, the compile time assistance in model debugging provided by structural analysis is lost with such approaches. On top of that, new difficulties can also arise at mode changes: Issue 3 (synthesis of restarts at mode changes) The designer specifies the model for each mode, but should not always be requested to specify the restart conditions resulting from a transition. How can these conditions be inferred? For example, restart equations for the clutch model (when entering the engaged mode) should be synthesized during compilation, not manually specified. For larger models, restart conditions may be very difficult to synthetize by hand, which makes Issue 3 of uttermost importance.

Last but definitely not least, we lack a reference semantics for multimode DAE systems in general:

Issue 4 (notion of solution) No notion of solution exists for general multimode DAE systems, but only partial notions on restricted subclasses of systems.

This came as a surprise to us when we started our study-the reader is referred to [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF] for a bibliographical discussion. The work of Stephan Trenn, relying on distributions, is an important contribution to the subject. However, Trenn himself points out in [START_REF] Trenn | Distributional Differential Algebraic Equations[END_REF][START_REF] Trenn | Regularity of distributional differential algebraic equations[END_REF] inherent difficulties to this approach, which indicates that distributions are not the ultimate answer to deal with impulsive variables in multimode DAE systems. Still, Trenn was able in [START_REF] Liberzon | Switched nonlinear differential algebraic equations: Solution theory, lyapunov functions, and stability[END_REF] to define complete solutions for a class of switched DAE systems in which each mode is in quasi-linear form: switching conditions are time-based, not state-based. The bottom line is that, unlike for ODE-based hybrid systems or DAE-based single-mode systems, we lack reference semantics. Moreover, spurious situations can arise, leading to infinitely fast mode switching [START_REF] Utkin | Sliding mode control in mechanical systems[END_REF]. The best we can hope is to prove that our approach boils down to known approaches for certain subclasses of systems: such an equivalence result was proved in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF] for the subclass of semi-linear DAE systems supported by Modiamath [START_REF] Elmqvist | Modiamath webpage[END_REF].

Contribution: In [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF], we addressed Issue 2 by proposing a scalable, compiletime, multimode structural analysis; in this paper, we present our approach to Issues 3 and 4. For this purpose, we need to give semantics to dynamics that combine discrete steps and continuous evolutions. This was already an issue for ODE-based hybrid systems; it becomes a steep challenge when moving to DAEs. The use of nonstandard analysis and hyperreals was already advocated in the study of ODE-based hybrid systems in the computer science community [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF][START_REF] Benveniste | Building a hybrid systems modeler on synchronous languages principles[END_REF]. While it was an "alternative approach" in that case, it turned out to become instrumental for the study of DAE-based hybrid systems.

The clutch example in existing physical modeling tools

The clutch model of Fig. 1 translates as the Modelica code of Fig. 2, which is accepted by the two Modelica tools we had the opportunity to test: OpenModelica 1.17.0 [START_REF] Fritzson | The OpenModelica Integrated Environment for Modeling, Simulation, and Model-Based Development[END_REF] and Dymola 2021 [START_REF]Dymola official webpage[END_REF]. However, none of these tools generates correct simulation code. Simulations fail precisely at the instant when the clutch switches from the released mode (g=false) to the engaged one (g=true), as evidenced by a division by zero exception resulting from the pivoting of a linear system of equations that becomes singular when g becomes equal to true.

We also tested the Mathematica Advanced Hybrid & DAE toolbox, a library advertised as supporting multimode DAE systems. With this toolbox, the simulation does not abort, but it returns incorrect results: the angular momentum of the system is not preserved, in contradiction with the physics; also, the state jump seems to be randomly selected: a small change in the velocities before the change results in totally different values for the restart.

model C l u t c h B a s i c parameter Real w01=1; parameter Real w02 = 1 . 5 ; parameter Real j 1 =1; parameter Real j 2 =2; parameter Real k1 = 0 . 0 1 ; parameter Real k2 = 0 . 0 1 2 5 ; parameter Real t 1 =5; parameter Real t 2 =7; Boolean g ( s t a r t=f a l s e ) ; Real w1 ( s t a r t=w01 , f i x e d=t r u e ) ; Real w2 ( s t a r t=w02 , f i x e d=t r u e ) ; Real f 1 ; Real f 2 ; e q u a t i o n g = ( time>=t 1 ) and ( time<=t 2 ) ; j 1 * der ( w1 ) = -k1 * w1 + f 1 ; j 2 * der ( w2 ) = -k2 * w2 + f 2 ; 0 = i f g then w1-w2 e l s e f 1 ; f 1 + f 2 = 0 ; end C l u t c h B a s i c ; Fig. 2. Modelica code for the idealized clutch. This is a faithful translation in the Modelica language of the model of Fig. 1 with fi as in [START_REF] Benveniste | Nonstandard semantics of hybrid systems modelers[END_REF]. The input guard γ (called g) takes the value t between t1 and t2, and f elsewhere.

Takeaway: Fundamental studies are clearly needed to ensure a correct handling of events of mode change by the Modelica tools. Smoothing "if then else" equations could help solving the problems arising at mode changes, but requires a delicate and definitely non-modular tuning, as this tuning depends on the different time scales arising in the system. We advocate that, as the tools reputedly support multimode DAE models, they should handle them correctly. Analysis of the clutch example: In Section 4, we perform the structural analysis of the model in each individual mode. This kind of symbolic analysis for acausal, DAE-based, models provided the grounds for the handling of mode changes that we propose in Section 5. The clutch example illustrates a major difficulty of this task, namely, the fact that restart conditions are, in general, not explicitly specified.

Structural analysis of each individual mode

In this section, we analyze separately the model for each mode of the clutch. In the released mode (γ = f in the model of Fig. 1), the two shafts are independent and one obtains the following two independent ODEs for ω 1 and ω 2 :

ω ′ 1 = f 1 (ω 1 , τ 1 ) (e 1 ) ω ′ 2 = f 2 (ω 2 , τ 2 ) (e 2 ) τ 1 = 0 (e 5 ) τ 2 = 0 (e 6 ) (5)
In the engaged mode, however (γ = t in the model of Fig. 1), the two velocities and torques are algebraically related:

ω ′ 1 = f 1 (ω 1 , τ 1 ) (e 1 ) ω ′ 2 = f 2 (ω 2 , τ 2 ) (e 2 ) ω 1 -ω 2 = 0 (e 3 ) τ 1 + τ 2 = 0 (e 4 ) (6)
Equation (e 3 ) relates the two velocities ω 1 and ω 2 that are otherwise subject to the ODE system (e 1 , e 2 ). This makes System (6) a DAE instead of an ODE.

In [START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF], torques τ 1 and τ 2 are not differentiated: they are algebraic variables; velocities ω 1 and ω 2 appear together with derivatives: they are state variables; ω ′ 1 , ω ′ 2 , τ 1 , τ 2 are the leading variables of System (6): they represent the unknowns. Suppose that the leading variables ω ′ 1 , ω ′ 2 , τ 1 , τ 2 are determined by given consistent values for the state variables ω 1 , ω 2 , i.e., values satisfying (e 3 ). Then, by using an ODE solver, one could perform an integration step and update the current velocities ω 1 , ω 2 using the computed values for their derivatives ω ′ 1 , ω ′ 2 . In this case, we say that the considered DAE is an "extended ODE" [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF].

It turns out that this condition does not hold for System [START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF] as is. To intuitively explain the issue, we move to sampled time T δ = def {0, δ, 2δ, 3δ, . . . }, by applying an explicit first-order Euler scheme with constant step size δ > 0 :

ω • 1 = ω 1 + δ f 1 (ω 1 , τ 1 ) (e δ 1 ) ω • 2 = ω 2 + δ f 2 (ω 2 , τ 2 ) (e δ 2
)

ω 1 -ω 2 = 0 (e 3 ) τ 1 + τ 2 = 0 (e 4 ) , (7)
where, for every t ∈ T δ ,

ω • (t) = def ω(t+δ) and • ω(t) = def ω(t-δ) (8)
respectively denote the forward and backward time shift operators by an amount of δ. Suppose we are given consistent values ω 1 = ω 2 and we wish to use System [START_REF] Benveniste | The synchronous languages 12 years later[END_REF], seen as a system of algebraic equations, to determine the value of the dependent variables (i.e., the unknowns)

τ 1 , τ 2 , ω • 1 , ω • 2 .
This attempt fails since we have only three equations e δ 1 , e δ 2 and e 4 to determine four unknowns τ 1 , τ 2 , ω • 1 and ω • 2 ; indeed, the omitted equation (e 3 ) does not involve any leading variable. Since System [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] is time invariant, if the system remains in the engaged mode for at least δ seconds, then the following shifted latent equation also holds:

ω • 1 -ω • 2 = 0 (e • 3 ) (9)
Replacing (e 3 ) with (e • 3 ) in System (7) yields a system with four equations and four dependent variables, which is nonsingular in a generic sense. One can now use System [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] augmented with equation (e •

3 ) to get an execution scheme for the engaged mode of the clutch. This is shown in Execution Scheme 1.

Execution Scheme 1 System (7)+Eq. [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF].

Require: consistent ω1 and ω2, i.e., satisfying (e3).

1:

Solve {e δ 1 , e δ 2 , e • 3 , e4} for (ω • 1 , ω • 2 , τ1, τ2) ▷ 4 equations, 4 unknowns 2: (ω1, ω2) ← (ω • 1 , ω • 2 ) ▷ update the states (ω1, ω2) 3: Tick ▷ move to next discrete step
Since the next values of the state variables satisfy (9) by construction, the consistency condition is met at the next iteration step.

Comment 1 (structural nonsingularity [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF]) The implicit assumption behind Line 1 in Execution Scheme 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 variables and/or parameters. 1We refer to this as structural nonsingularity in the sequel. □

Observe that the same analysis could be applied to the original continuous-time dynamics from System (6) by replacing (e 3 ) by the differentiated latent equation

ω ′ 1 -ω ′ 2 = 0 (e ′ 3 ) (10)
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 Execution Scheme 2, which parallels Execution Scheme 1.

Execution Scheme 2 System (6)+Eq. [START_REF] Campbell | The index of general nonlinear DAEs[END_REF].

Require: consistent ω1 and ω2, i.e., satisfying (e3).

1: Solve {e1, e2, e ′ 3 , e4} for (ω ′ 1 , ω ′ 2 , τ1, τ2) ▷ 4 equations, 4 unknowns 2: ODESolve (ω1, ω2) ▷ update the states (ω1, ω2) 3: Tick ▷ move to next discretization step Line 1 is identical for the two schemes and is assumed to give a unique solution, generically; it fails if one omits the latent equation (e ′ 3 ). Then, although getting the next values for ω 1 and ω 2 is easy in Execution Scheme 1, the situation changes in Execution Scheme 2: the derivatives ω ′ 1 , ω ′ 2 are first evaluated, then an ODE solver (here denoted by ODESolve) is used to update the state (ω 1 , ω 2 ) until the next time horizon.

Comparing Line 1 of Execution Scheme 1 and Line 1 of Execution Scheme 2 reveals the common principle driving the execution of dAE/DAE-based models:2 Principle 1 (dAE/DAE execution) Assume an algebraic equation system solver. The key task is to massage the source model in such a way that the system of equations determining the leading variables becomes regular (i.e., possesses, locally, a unique solution).

This transforms the original DAE system F (X ′ , X, Y ) = 0 into a form equivalent to an ODE X ′ = G(X, Y ); the same holds for the discrete-time counterpart.

Differentiation index and structural analysis of DAEs:

The replacement of (e 3 ) by (e ′ 3 ), which resulted in Execution Scheme 2, is known in the literature as index reduction [START_REF] Mattsson | Index reduction in differential-algebraic equations using dummy derivatives[END_REF]. It requires adding the (smallest set of) latent equations needed for the system in Line 1 of the execution scheme to become solvable and deterministic. The number of successive time shifts or differentiations needed for obtaining all the latent equations is called the differentiation index [START_REF] Campbell | The index of general nonlinear DAEs[END_REF]. Generally, for any given DAE model S, which is a system of numerical equations f i (the x j and their derivatives) = 0 , i = 1, . . . , m; j = 1, . . . , n , [START_REF]Dymola official webpage[END_REF] the structural analysis consists in producing Execution Scheme 2 associated to this model. By the implicit function theorem, the existence and uniqueness of solutions of System [START_REF]Dymola official webpage[END_REF] rely on the invertibility of the Jacobian matrix collecting the partial derivatives ∂f i /∂u j , where u j is the leading variable associated to the real variable x j .

An approximate analysis of System [START_REF]Dymola official webpage[END_REF] was developed since the 1980's that guarantees the regularity of this Jacobian matrix in a structural sense, taking only into account the occurrence/absence of a variable in an equation and its differentiation order: this data is stored into the weighted incidence graph G S of the system, which is a labeled bipartite graph in which edge (f i , x j ) is present with weight d ij iff variable x j occurs in f i with degree d ij [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF][START_REF] Mattsson | Index reduction in differential-algebraic equations using dummy derivatives[END_REF][START_REF] Pryce | A simple structural analysis method for DAEs[END_REF]. In [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], the problem of finding the offsets of each equation f i (i.e., the number of times it should be differentiated) was reduced to solving a certain linear program associated to the weighted incidence graph of System [START_REF]Dymola official webpage[END_REF]; this is known as the Σ-method. Details are found in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF].

Handling mode changes

Caveat: What is the semantics of a multimode DAE? The semantics of ODEs (i.e., the sets of runs or trajectories they define) has been known for a long time, as existence and uniqueness theorems for their solutions are well known in mathematics. The same holds for DAEs, although the corresponding developments are much more recent [START_REF] Campbell | The index of general nonlinear DAEs[END_REF]. In formal studies on hybrid systems, defining the semantics of a model is usually not an issue, hence, the scientific community can focus on difficult decidability issues for various classes of properties and models.

When entering the domain of multimode DAE systems, we did not expect semantics to be a problem. Many physical subdomains involve multimode DAE for their models. Electrical circuits with ideal diodes [START_REF] Heemels | On the dynamic analysis of piecewise-linear networks[END_REF], contact and friction mechanics [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Wiley[END_REF], and hydraulics with multiple-phase materials, are examples of this. For each of these domains, solutions of the corresponding multimode DAEs are well defined, albeit with reference to the dedicated formalism used to describe them-these formalisms differ among the various fields of physics.

In contrast, multiphysics modeling combined with the modeling of embedded digital control is beyond the scope of such dedicated modeling approaches. A notable exception stems from bond graphs [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF], where generic notions of effort, flow, junction, and gyrator generalize the notions of potential, current, junction, and motor in electrical circuits. This yields a generic multiphysics formalism well suited to DAE-based modeling. However, so-called switch bond graphs, extending bond graphs to multimode systems, are not well developed [START_REF] Broenink | Describing discontinuities in bond graphs[END_REF]. We are therefore faced with the following demanding problem: Problem 1 (Semantics). Given that there is no reference semantics in mathematics for general multimode DAE systems, find a sensible subclass of multimode DAE systems:

1. that is rich enough for supporting complex multiphysics systems modeling, 2. so that checking that a given system belongs to this subclass can be checked at compile time, and 3. for which a semantics is well defined.

An example of property that does not meet requirement 2 is "the class of systems whose solution trajectories are continuous with continuous derivatives". Our work [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF] provides a systematic answer to Problem 1. We illustrate our approach on the clutch example in the next section.

Structural analysis at mode changes

As discussed above, discontinuous velocities and impulsive torques are expected when the clutch gets engaged. To analyze this, we reuse the discrete-time approximation [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] with time step δ and the forward and backward shift operators introduced in (8). To capture the transition γ : f → t, we unfold the dynamics of the clutch over two successive instants, by adding the latent equation (e •  3 ) associated to mode γ=t:

previous instant γ = f          ω 1 -• ω 1 δ = f1( • ω1, • τ1) ( • e δ 1 ) ω 2 -• ω 2 δ = f2( • ω2, • τ2) ( • e δ 2 ) • τ1 = 0 • τ2 = 0 current instant γ = t                ω • 1 -ω 1 δ = f1(ω1, τ1) ω • 2 -ω 2 δ = f2(ω2, τ2) ω1 -ω2 = 0 (e3) ω • 1 -ω • 2 = 0 (e • 3 ) τ1 + τ2 = 0   subsystem ( • e δ 1 , • e δ 2 , e3) is conflicting   (12)
Subsystem ( • e δ 1 , • e δ 2 , e 3 ) possesses 3 equations, but only 2 dependent variables ω 1 , ω 2 as the remaining variables are all determined at the previous time step. This subsystem is thus structurally singular (over-specified). Should the clutch model be rejected? Not quite. This problem is only an artifact of time discretization.

We will get rid of the conflict in model ( 12) by removing some equations from the conflicting subsystem, at the instant of mode change. To respect causality, the past will not be undone, so that equations ( • e δ 1 , • e δ 2 ), which were handled at the previous instant, cannot be modified. However, we can resolve the conflict by erasing equation (e 3 ) in System [START_REF] Elmqvist | A structured model language for large continuous systems[END_REF]. Thus, at the instant of mode change, the consistency equation (e 3 ) from the new mode will not be satisfied. At the next instant, however, consistency equation ω 1 -ω 2 = 0 will be guaranteed by equation (e • 3 ) of ( 12), and it will then remain satisfied as long as the system stays in the engaged mode. Hence, erasing (e 3 ) in System [START_REF] Elmqvist | A structured model language for large continuous systems[END_REF] amounts to postponing by one time step the satisfaction of this state constraint.

Acceptance/rejection of a model: The analysis of the released → engaged mode change succeeds for the clutch example. It might have failed if erasing contradicting equations in the present and future of the mode change had yielded a structurally under-specified system, expressing that the model would be underspecified at this mode change-see the cup-and-ball example in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF].

Discussion: The analysis developed here meets the requirements of Problem 1: the acceptance/rejection of the submitted model relies on a compile-time structural analysis. For an accepted model, the δ-discretized scheme defines the semantics, which is well defined in a generic sense.

Working with a fixed step size discretization of time is, however, unacceptable, even for single-mode ODEs, as stiff dynamics, whose presence cannot be checked at compile time, can always occur. In the next section, we address this difficulty by taking for our time step an "infinitesimal" number, thus getting "infinitesimal" accuracy for our discretization. This is achieved by calling nonstandard analysis to the rescue.

Nonstandard analysis and hyperreals to the rescue

If DAE dynamics is approximated in discrete time, then the whole model becomes discrete-time. To avoid the problem of approximation error, our idea is to use an "infinitesimal" time step in the discrete-time approximation. This will yield an approximation that is accurate up to an infinitesimal error. This can be made rigorous by relying on nonstandard analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], which extends the set R of real numbers to a superset ⋆ R of hyperreals that includes infinitely large and infinitely small numbers. For the understanding of this paper, it is enough to know the following about nonstandard analysis.

A glimpse of nonstandard analysis: There exist infinitesimals, defined as hyperreals of ⋆ R that are smaller in absolute value than any nonzero real number. The usual arithmetic operations and relations are lifted to ⋆ R. Say that x ∈ ⋆ R is finite if there exists a real number y ∈ R + such that |x| ≤ y.

For every finite hyperreal x ∈ ⋆ R, there is a unique standard real number st(x) ∈ R such that st(x) -x is infinitesimal; st(x) is called the standard part (or standardization) of x. [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF] Standardizing functions or systems of equations, however, requires some care. One important issue is derivatives. For t → x(t) an R-valued signal (t ∈ R),

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 is infinitesimal; then, a = x ′ (t). [START_REF] Elmqvist | Modiamath webpage[END_REF] Nonstandard time: For ε > 0 an infinitesimal time step, we consider the time index set T ⊆ ⋆ R:

T = 0, ε, 2ε, 3ε, • • • = {nε | n ∈ ⋆ N} ( 15
)
where ⋆ N denotes the set of hyperintegers, consisting of all integers augmented with additional infinite numbers called nonstandard. The following features of T are important: (1) any finite real time t ∈ R is infinitesimally close to some element of T (hence, T covers R and can be used to index continuous-time dynamics); and (2) T is "discrete": every instant nε has a predecessor (n-1)ε (except for n = 0) and a successor (n+1)ε. Let x be a nonstandard signal indexed by T. The forward-and backward-shifted signals x • and • x are defined as in [START_REF] Broenink | Describing discontinuities in bond graphs[END_REF], with ε replacing δ.

Nonstandard semantics: Solutions of multi-mode DAE systems may be nondifferentiable or even non-continuous at events of mode change. To give a meaning to x ′ at any instant, we define it everywhere as

x ′ = def 1 ε (x • -x) , (16)
which agrees with ( 14) at the instants where x is differentiable. Consider again System [START_REF] Elmqvist | A structured model language for large continuous systems[END_REF] where conflicting equation (e 3 ) is removed, and focus on the dynamics in the current instant, where derivatives are expanded using [START_REF] Heemels | On the dynamic analysis of piecewise-linear networks[END_REF]:

   ω • 1 = ω 1 + ε f 1 (ω 1 , τ 1 ) ω • 2 = ω 2 + ε f 2 (ω 2 , τ 2 ) ω • 1 -ω • 2 = 0 and τ 1 + τ 2 = 0 (17)
We wish to use System [START_REF] Henzinger | HYTECH: The Cornell HYbrid TECHnology Tool[END_REF] by identifying current values for the states ω i with the left-limits ω - i i.e., the values of the velocities just before the mode change. From these values, we would then compute the restart values for the velocities ω + i = def ω • i , together with the torques τ i .

From hyperreals to effective code

Unfortunately, hyperreals are unknown to computers; hence, System (17) cannot be used as such, but needs to be standardized, by "washing out" the infinitesimal time step ε. As it is infinitesimal, it is tempting to get rid of of it in [START_REF] Henzinger | HYTECH: The Cornell HYbrid TECHnology Tool[END_REF] by simply setting ε = 0. Unfortunately, doing this leaves us with a structurally singular system, since the two torques are then involved in only one equation. 3This problem of structural singularity is in fact due to the existence of impulsive variables. To discover them in a systematic way, we perform an impulse analysis.

Impulse analysis: When the clutch is released, we have ω 1 -ω 2 ̸ =0, generically.

Since

ω • 1 -ω • 2 = 0 holds, (ω • 1 -ω • 2 )-(ω1-ω2) ε = f 1 (ω 1 , τ 1 )-f 2 (ω 2 , τ 2
) cannot be finite because, if it was, then the function ω 1 -ω 2 would be continuous, contradicting the assumption that ω 1 -ω 2 ̸ =0. Hence, the hyperreal 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 mode change, and because of equation (e 4 ) : τ 1 + τ 2 = 0, both torques are in fact infinite, i.e., are impulsive.

Eliminating impulsive variables: We now assume that each f i has the form (3). Erasing (e 3 ) from the black part of (12) yields the following system of equations, to be solved for ω

• 1 , ω • 2 , τ 1 , τ 2 at the instant when γ switches from f to t:        J 1 ω • 1 = J 1 ω 1 + ε(τ 1 -g 1 (ω 1 )) (e ε 1 ) J 2 ω • 2 = J 2 ω 2 + ε(τ 2 -g 2 (ω 2 )) (e ε 2 ) ω • 1 -ω • 2 = 0 (e • 3 ) τ 1 + τ 2 = 0 (e 4 ) (18)
We eliminate the impulsive variables from System (18), namely, the two torques.

Using (e 4 ) yields -τ

2 = τ 1 = def τ . Using (e • 3 ) by setting ω • = def ω • 1 = ω • 2
and adding the two equations (e ε 1 ) and (e ε 2 ) yields

(J 1 + J 2 )ω • = J 1 ω 1 + J 2 ω 2 -ε(g 1 (ω 1 ) + g 1 (ω 2 )) .
It is now legitimate to set ε to 0 to standardize the right-hand side of this equation. Finally, identifying st(ω i )=ω - i and st(ω • i )=ω + i synthesizes formula (4). Alternative approaches for the computation of the reset values, which do not require the elimination of impulsive variables, are presented in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF].

We thus recovered the restart condition derived from conservation laws of the mechanics in Section 2.1, by only using physics-agnostic symbolic reasoning. This addresses Issue 3 of Section 2.2, and can be made fully automatic [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF].

Discussion and perspectives

A new modeling paradigm: In this paper, we advocate the use of multimode DAE systems (or DAE-based hybrid systems) for physical modeling and, thus, CPS modeling. We argued for its stringent need to bring modularity and reuse in modeling. We illustrated on a toy example the difficulties in compiling such models to derive simulation code, and proposed effective solutions. Major lessons from this study are the following:

-Since DAEs, not ODEs, are involved, we must rely on structural analysis for both all modes and mode changes of the system. -In many cases, restart conditions at mode changes can be synthesized from a proper structural analysis. -The compilation steps detailed above provide a compile-time analysis for accepting or rejecting models, based on over/under-specifications of it, within modes or at mode changes. -An impulse analysis is also presented in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF]; this is a compile-time symbolic analysis identifying impulsive variables, which allows for the handling of models involving impulsive behaviors. -Hyperreals and nonstandard analysis were instrumental in doing this: they allow to unify the model of time, combining continuous-time dynamics within long modes and discrete-time steps for computing restarts at mode changes.

Our approach gives a semantics to multimode DAE systems possessing long modes and finite cascades of mode switches. Note that dense infinite cascades, such as those resulting from sliding mode control [START_REF] Utkin | Sliding mode control in mechanical systems[END_REF], are not covered. Our semantics was proved to coincide with the semantics provided by the physics for the restricted class of semilinear systems, to which the clutch belongs [START_REF] Benveniste | Multi-mode DAE models -challenges, theory and implementation[END_REF][START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF].

No comparison can be stated for our method in full generality, however, due to the lack of reference semantics for multimode multiphysics systems. A full justification of the use of nonstandard analysis and standardization was provided in [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF].

Making our approach systematic with algorithms that properly scale up: In the grounding paper [START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF], we develop a comprehensive theory handling, in full generality, multimode DAE systems belonging to the class described above. The whole procedure is automatic, and effective algorithms are specified and mathematically justified.

In this paper, we did not discuss algorithmic complexity. For the clutch example, structural analyses were needed for each of the two modes (released and engaged), and each of the two transitions (released → engaged and vice-versa). Since the number of modes can grow exponentially with the number of subsystems, the mode-centric modeling approach associating a DAE model to each different mode cannot scale up. This can be overcome in practice by adopting the dual equation-centric approach, in which each equation is labeled with a predicate characterizing the set of modes in which this equation is active. 4 This dual approach, combined with an efficient encoding of the multimode Σ-method using Binary Decision Diagrams (BDD) techniques, proved to scale up nicely on significant experiments [START_REF] Caillaud | Implicit structural analysis of multimode DAE systems[END_REF]. A large part of this body is implemented in our tool IsamDAE. 5 The examples coming with this tool already include thermodynamical, electrical and pneumatic models. Work is in progress toward developing modular structural analysis as a means to further scale up, and possibly achieve separate compilation of such models.

Grand challenges: Let us wrap up this paper by proposing a few grand challenges for computer scientists in this area.

Correct compilation of multimode DAE systems: With the current reference tools, physically correct models need to be nontrivially twisted and massaged to get simulated, which can require assistance by tool experts. Our work addresses this challenge. Scaling up: Large industrial models can involve up to hundreds of thousands of equations; examples are found in the energy sector. Such models are not fully supported by the existing tools. Progress is stringently needed on this issue, and should involve both high performance computing and compiletime analyses such as the ones developed above. Modularity and separate compilation: The first step of model compilation in current tools is a flattening of the model structure. This can yield prohibitive computational costs, particularly for large systems built from many instantiations of a few component classes. As such, separate compilation of components is stringently needed. Initialization of multimode DAE models: Specifying correct initial conditions is a widely acknowledged difficulty for designers of large DAE-based models, especially for CPS. Significant work was devoted to initialization for (singlemode) DAEs, but this remains a mostly open issue for multimode DAE models. General concept of solutions for multimode DAEs: This is a core challenge from a mathematical viewpoint, that we partly addressed.

Numerical singularity occurs when the Jacobian matrix of the system is singular; in this specific case, this can only occur when ∂f1(ω1, τ1)/∂τ1 + ∂f2(ω2, τ2)/∂τ2 = 0. Unlike structural singularity, numerical singularity depends on numerical values of both coefficients and variables.

dAE stands for difference Algebraic Equation (related to discrete time), while DAE stands for Differential Algebraic Equation (related to continuous time).

Actually, this is an example of erroneous standardization of systems of equations. In[START_REF] Benveniste | The mathematical foundations of physical systems modeling languages[END_REF][START_REF] Benveniste | Structural analysis of multimode DAE systems: summary of results[END_REF] it is proved that standardizing a system of equations by setting ε = 0 is correct if and only if doing so yields a structurally regular system of equations.

This equation-centric viewpoint was first used in the clock calculus and conditional dependency graph used in the compilation of the Signal synchronous language[START_REF] Benveniste | The synchronous languages 12 years later[END_REF].

To date, the following compilation tasks are implemented: multimode Σ-method, multimode initialization, and multimode Dulmage-Mendelsohn decomposition.

This work was supported by the FUI ModeliScale DOS0066450/00 French national grant (2018-2021) and the Inria IPL ModeliScale large scale initiative (2017-2021, https://team.inria.fr/modeliscale/).