Albert Benveniste

Benoît Caillaud

Hilding Elmqvist

Khalil Ghorbal

Martin Otter

Marc Pouzet

Multi-Mode DAE Models Challenges, Theory and Implementation

Keywords: Multi-mode systems, Cyber-Physical Systems, CPS, Modia, Modelica, differential algebraic equations, DAE, differential index, structural analysis, operational semantics, constructive semantics, nonstandard analysis

Our objective is to model and simulate Cyber-Physical Systems (CPS) such as robots, vehicles, and power plants. The structure of CPS models may change during simulation due to the desired operation, due to failure situations or due to changes in physical conditions. Corresponding models are called multi-mode. We are interested in multidomain, component-oriented modeling as performed, for example, with the modeling language Modelica that leads naturally to Differential Algebraic Equations (DAEs). This paper is thus about multi-mode DAE systems. In particular, new methods are discussed to overcome one key problem that was only solved for specific subclasses of systems before: How to switch from one mode to another one when the number of equations may change and variables may exhibit impulsive behavior? An evaluation is performed both with the experimental modeling and simulation system Modia, a domain specific language extension of the programming language Julia, and with SunDAE, a novel structural analysis library for multi-mode DAE systems.

Introduction

Modeling with block diagrams described with Ordinary Differential Equations in state space form (ODE) has become a key pillar in the development of Cyber-Physical Systems (CPS). Block diagrams, however, suffer from a lack of modularity and reuse that is best illustrated in Fig. 1. This figure shows two models of the same system consisting of a simple model of an electrical motor and of the rotational inertia of the motor. On the left hand side a (Modelica) schematic/component diagram of the system is shown that connects physical components through non-directed interactions resulting from the first principles of physics. On the right hand side the same model is shown as a block diagram in which input/output oriented blocks are connected with a directed wiring manually specified by the designer. Adding one more physical component is straightforward in the schematic, whereas it may need a complete redesign in the block diagram.

This should not come as a surprise, as the first principles of physics naturally lead to considering acausal models such as the one in Fig. 1-left. Consider, for example, the case of electric circuits. So-called circuit laws such as Kirchhoff laws, are naturally expressed as balance equations: the algebraic sum of currents in a network of conductors meeting at a point is zero; or, the sum of all the voltages around a loop is equal to zero. Similarly, some components (such as, e.g., resistors or capacitors) come with no input/output prespecified orientation. A same circuit can be assigned different input/output status for its variables, depending on which ones are declared as sources. The same situation arises in mechanics or in thermodynamics. Engineers interested in multi-physics modeling have identified this fact since the 70's by proposing bond graphs [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF][START_REF] Karnopp | System dynamics: a unified approach[END_REF], in which electric circuits, mechanical systems, and thermodynamical systems, are abstracted to a common framework manipulating efforts, flows, and junctions. To summarize, for systems made of a large number of interconnected components, getting the usual input/output state space model (such as in Simulink block diagrams) becomes intractable, whereas it remains manageable if acausal models are supported.

In addition to being naturally described by an acausal model, a system may have different modes for different reasons. One reason is the physics: 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, changing operations of a system like a power plant with start-up, normal operation and shut down. Finally, safety needs may lead to so-called degraded modes in which the system dynamics may become very different (some components, sensors or actuators getting down for example). Models of this kind are called multi-mode models.

Acausal multi-mode models present, however, subtle difficulties at mode changes. Consider for instance the case of an idealized clutch model, illustrated in Fig. 2. The clutch possesses two modes: released (the two shafts rotate freely) Fig. 2. A simple clutch with two shafts.

and engaged (the two shafts are in contact). Deriving acausal models from the first principles is simple and rather elegant for each single mode alone. The intuition tells us that, when the clutch gets engaged, the different velocities of the two shafts will, in zero time, merge to a unique identical angular velocity under impulsive torques. One also expects that the resulting common velocity is entirely determined and depends only on the individual inertia and velocity of each shaft, prior to engagement. Determining manually the restart conditions is required if one resorts to using oriented block-diagrams like in Simulink. We argue that such task is difficult for this simple example and becomes even impossible for real-life models. Those restart conditions should rather be synthesized automatically by the compiler and made available at run time to simulate the model. This work advocates therefore the use of acausal modelling languages since they 1. make the specification of large (single mode) Cyber Physical Systems more modular, more elegant, and with better reuse; 2. are indispensable to synthesize proper restart conditions at mode changes, for models where this cannot be done manually -the idealized clutch model being a simple example.

It should also be clear that classical state space input/output formalisms (such as Simulink block diagrams) do not address the above challenges. Neither does the formalism of hybrid automata [START_REF] Alur | Hybrid automata: an algorithmic approach to the specification and verification of hybrid systems[END_REF] used as the model of hybrid systems by several verification tools such as Flow* or SpaceEx. Causal input/output state space models are based on Ordinary Differential Equations (ODE), of the form ẋ = f (x, u) where x is the state and u the input, whereas acausal models à la Modelica are based on Differential Algebraic Equations (DAE) of the form f ( ẋ, x, v) = 0 where x is a state and v an algebraic variable, for which no input or output status can be a priori specified.

In this paper we propose a theory to support the compilation of multi-mode DAE based models in a systematic and mathematically sound way. We want to reject some models that are in some sense spurious. And we want to synthesize (whenever needed, not in all cases) the computation of the restart conditions at mode changes. All of this will be performed at compile time, through a special static analysis called the structural analysis, which is an important topic of this paper. This structural analysis is an essential preprocessing step prior to generating simulation code. We will illustrate our purpose by an example from electro-mechanics that includes a clutch as one of its subsystems. The specification, analysis, compilation, and simulation of this model will be illustrated using the recently proposed experimental modeling language called Modia.

State-of-the-art and Related Work

There exists an extensive body of work related to the numerical solution of single-mode DAEs, especially [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF], as well as to the modeling and simulation of single-mode multi-domain models leading to DAEs5 .

Mechanical systems featuring impulses are known since the 18th century. A large literature about multi-body systems with contacts and impulses exists, for example [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Pfeiffer | On non-smooth multibody dynamics[END_REF]. The current research focuses on simulating contacts between many bodies using time stepping methods. Contact equations are usually described approximately on velocity or acceleration level (if contact occurs, the relative velocity or acceleration between the relevant bodies is constrained to be zero). Since the constraints on position level are not explicitly taken into account, typically a drift occurs that is usually not relevant for systems with many contacting bodies but is not acceptable for general multi-mode DAEs. An exception is [START_REF] Schoeder | Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF], where constraints on position level are enforced. It is not obvious how the specialized multi-mode methods for multi-body systems can be generalized to any type of multi-mode DAE.

Also for electrical circuits with idealized switches, like ideal diodes, impulses can occur. Again a large literature is available concentrating mostly on specialized electrical circuits such as piecewise-linear networks with idealized switches [START_REF] Heemels | On the dynamic analysis of piecewise-linear networks[END_REF][START_REF] Trenn | Distributional Differential Algebraic Equations[END_REF][START_REF] Barela | A complementarity approach to modeling dynamic electric circuits[END_REF]. Again, it is not obvious how to generalize methods in this specialized area to general multi-mode DAEs.

The paper of Mehrmann et al. [START_REF] Mehrmann | Hybrid systems of differential-algebraic equations -analysis and numerical solution[END_REF] contains interesting results regarding numerical techniques to detect chattering between modes. It, however, assumes that consistent reset values are explicitly given for each mode. Such an assumption does not hold in general, especially for a compositional framework where one wants to assemble pre-defined physical components.

In the PhD-thesis of Zimmer [START_REF] Zimmer | Equation-Based Modeling of Variable-Structure Systems[END_REF] variable-structure multi-domain DAEs are defined with a special modeling language and a run-time interpreter is used that processes the DAE equations at run-time, when the structure and/or the index is changing. Limitations of this work are that impulsive behavior is not supported and that the user has to define the transfer of variable values from one mode to the next mode explicitly, which is not practical for large models.

Describing variable structure systems with causal state machines is discussed by Pepper et al. [START_REF] Pepper | A Compositional Semantics for Modelica-style Variable-structure Modeling[END_REF].

Elmqvist et al. [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF][START_REF] Mattsson | Multi-Mode DAE Systems with Varying Index[END_REF] propose a high level description of multi-mode 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, also acausal DAE models with physical connectors can be a "state" of a state machine. Such a state machine is mapped into a form so 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 multi-mode systems can be handled due to the static code generation.

Benveniste et al. [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF][START_REF] Benveniste | Structural Analysis of Multi-Mode DAE Systems[END_REF] tackle the problem of variable structure, varying index DAEs from a fundamental point of view by using a low level description of DAEs with a few language elements only and a precise mathematical description of the semantics based on non-standard analysis. A proof-of-concept mockup, SunDAE, was developed implementing this approach.

The Modia Language

Modelica [START_REF] Modelica | A Unified Object-Oriented Language for Systems Modeling[END_REF] is a state-of-the-art modeling language for multi-domain modeling. Recently, an experimental language which is similar to Modelica and called Modia [START_REF] Elmqvist | Systems Modeling and Programming in a Unified Environment based on Julia[END_REF][START_REF] Elmqvist | Innovations for Future Modelica[END_REF] has been designed and implemented. Modia is a domain specific language extension of Julia [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF] by means of structured macros, that is, the Julia parser is used to parse Modia models. The Modia language elements will be introduced through a small set of examples.

Elementary Modia Constructs A model of a rotating inertia can be defined as follows (# start a comment until end of the line): @model Inertia begin J = 1 w0 = 0.0 w = Float(start = w0) flange = Flange() @equations begin w = flange.w J*der(w) = flange.tau end end

The @model macro has one section for declarations and one section for equations and connections. The construct J = 1 defines the parameter J with default value of 1. A variable w with a start value of w0 is introduced by w = Float(start = w0). The construct flange = Flange() calls the constructor for Flange which is defined as follows: @model Flange begin w = Float() # Angular velocity tau = Float(flow=true)

# Torque end This model consists just of two variable declarations with type Float modeling the interaction when the Inertia is mechanically coupled via a shaft to other mechanical components. The attribute flow=true tells that if several shafts are connected to flange, the internal and external torques are all summed to zero.

The equation J*der(w) = flange.tau is Euler's equation for the rotational motion of the inertia. The operator der() denotes time derivative.

A model with two rigidly connected inertias can be defined as shown below:

@model TwoInertias begin inertia1 = Inertia(J=0.1) inertia2 = Inertia(J=0.4, w0=10.0) @equations begin connect(inertia1.flange, inertia2.flange) end end Note that parameters are changed in the Inertia constructor call. The connect() primitive models a physical coupling. In this case the semantics is: Note, that inertia2 has initial angular velocity of 10 rad/s, whereas inertia1 has the default of 0. The result of simulation shows a constant angular velocity for both inertias of 8 rad/s. Since the initial angular velocities are different there will be a Dirac impulse in the torque to make them equal. The resulting joint angular velocity is according to theory weighted as w = (J1 * w1+J2 * w2)/(J1+J2) = 8.

Advanced Modia Constructs Modia has many other features, such as inheritance, size and type inference, type declarations, component redeclarations, time events, state event, and nested simulations. More information can be found in [START_REF] Elmqvist | Systems Modeling and Programming in a Unified Environment based on Julia[END_REF][START_REF] Elmqvist | Innovations for Future Modelica[END_REF]. It is possible to declare variables with specific SI units. As an example, an electric Pin connector can be defined as shown below using predefined variable constructors Voltage() and Current() with associated SI units (Volt and Ampere):

@model Pin begin v = Voltage() i = Current(flow=true) end
Modia Multi-mode Modeling: Running Example We model in Modia a multimode DAE system containing the idealized clutch introduced informally in Section 1. The model has an electric motor, two load inertia and a clutch, see Fig. When the clutch is engaged, the angular velocities of the two flanges are the same and the torques sum to zero. When not engaged, the torques are zero and there are no constraints on the angular velocities. The definition of the input engaged is made in the surrounding environment of the clutch:

clutch.engaged = time < 100 || time >= 300
The DAE index is changing depending on the state of the clutch. Furthermore, Dirac impulses occur initially and when the mode changes to engaged. The techniques needed to analyze and simulate such models are described in subsequent sections.

The results of a simulation are shown in Figure 4. The initial angular velocity of inertia1 is 0 rad/s and 10 rad/s for inertia2. The capacitor is initially uncharged. The upper two partly overlapping curves are the angular velocities of inertia1 and inertia2. When the clutch is disengaged, the angular velocity of inertia2 is constant. The lower curve shows the voltage over the capacitor. Since the clutch is engaged at initialization, Dirac impulses occur at time=0 s. As a result, after initialization the angular velocities of the two inertia are identical. The common angular velocity at time=0 s is not 8 rad/s as in the example TwoInertias, but 6.4 rad/s. The reason is that the capacitor acts in the

Simulating a Restricted Class of Multi-Mode DAEs

In this section, we propose a method to simulate a restricted class of multimode DAE systems, encompassing in particular our example of Figure 3 and its simulation result in Figure 4.

Throughout this paper, we use the classical notations ẋ and ẍ to denote the first and second time-derivatives of a function of continuous time x : R → R. When such a notation is not appropriate for readability reasons, we sometimes write instead x and x . The notation for higher order derivatives is indicated when needed. We write vectors and matrices in boldface. The transpose of a matrix G is written G T .

Problem setting

The goal is to simulate the following class of multi-mode DAE systems:

if γ 1 ( ẋ- γ , x - γ , t) then f 1 ( ẋ1 , x 1 , t) = 0 elseif γ 2 ( ẋ- γ , x - γ , t) then f 2 ( ẋ2 , x 2 , t) = 0 . . . else f m ( ẋm , x m , t) = 0 (1)
where vector x γ of length n γ collects variables that are present in all modes and in the predicates γ i (. . . ), vector x i of length n i collects all variables that are used in mode i (x γ is a subset of x i ), and

-γ i ∈ R nγ × R nγ × R → Boolean, f i ∈ R ni × R ni × R → R ni .
Hence the considered system switches between different modes, each described by a DAE having n i equations and n i unknowns. The switching conditions are predicates denoted by γ i and depend on time and/or on the left limit of the DAE variables assuming they exist and are finite, that is

x - γ (t) = def lim ε 0 x γ (t -ε) and ẋ- γ (t) = def lim ε 0 ẋγ (t -ε) .
When changing from mode i to mode j it is assumed that the initial conditions x - γ are known from the previous mode. A concrete example will be given in section 4.3. we also refer to [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF] for a more generic treatment. When changing from one mode to another the set of equations that govern the dynamics of the system may change. Thus, some variables may experience discontinuities or even impulses. Throughout this section we assume that the system spends a strictly positive duration in each mode in such a way that instants of mode changes are cleanly isolated.

Note that the class (1) of multi-mode systems is not compositional as the predicates x γ are global. We restricted ourselves to this class in order to simplify the exposure below. This restriction is relaxed in the actual implementation in Modia. For instance, @model Clutch in section 3 has a local definition of the predicate to engage or disengage the clutch and there is no global if-clause. In section 5, where an alternative approach is discussed, such restrictions are not present.

Handling Mode Changes

In a mode i, an initial value problem f i ( ẋi , x i , t) = 0 has to be solved. It is well-known that only special classes of DAEs can be directly numerically integrated [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF]. For this reason, general DAEs might be first transformed in to a DAE class where reliable integration methods exist. For nonlinear systems such a transformation means that equations of the DAE might need to be (analytically) differentiated. The best integration methods exist for ODEs ẋ = f (x, t). However, transforming a DAE to this form requires in general to solve nonlinear algebraic equation systems at every model evaluation and therefore the benefit of ODE integration methods might get lost. For this reason, another approach is used in this section where the transformation to a numerically solvable form is performed without solving algebraic equation systems. In Assumption 1 we define the target DAE class for this transformation: Assumption 1 For each mode i of System (1), the DAE f i ( ẋi , x i , t) = 0 has the following form (we omit the index i for better readability)

f d ( ẋ, x, t) = 0 f c (x, t) = 0 (2)
where

J =   ∂f d ∂ ẋ ∂f c ∂x   is regular (that is, invertible). (3)
In (2), subscripts " d " and " c " are reminiscent of "dynamics" and "constraint", respectively. Condition (3) means that System (2) has differentiation index one. 6Transforming general DAEs to System (2,3). If the model equations in a given mode are not in the form (2,3), they are transformed to it. This transformation is non-trivial and it is beyond the scope of this paper to explain the details. Only a short overview is given here: Typically, the Pantelides algorithm [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF], or a variant like [START_REF] Pryce | A simple structural analysis method for DAEs[END_REF], is used to determine which equations of the original DAE must be differentiated in order that the highest derivative variables can be uniquely determined from the highest derivative equations. By differentiating the corresponding equations analytically it is then possible to transform to ODE form. Hereby, algebraic equation systems might need to be solved. In [START_REF] Otter | Transformation of Differential Algebraic Array Equations to Index One Form[END_REF] a new algorithm is proposed to transform every DAE that can be treated with the Pantelides algorithm to the form of Assumption 1 without solving algebraic equation systems. This algorithm is a generalization of [START_REF] Gear | Automatic integration of eulerlagrange equations with constraints[END_REF] that was developed for single-mode multi-body systems. In section 4.3 it is applied to multi-mode multi-body systems.

Simulating System (2,3). A number of methods exist for solving system (2,3) numerically. In particular, fixed or variable step-size BDF (Backward Differentiation Formula) methods can be used [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF]. In all cases, consistent initial conditions have to be provided that fulfill (2) and the so-called latent equations

∂f c ∂x (x, t) ẋ + ∂f c ∂t (x, t) = 0 (4)
obtained by differentiating the constraint f c (x, t) = 0, see [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF][START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF]. Integration methods, including BDF, assume that x is smooth. Therefore, standard integration methods cannot be applied if x is discontinuous at a mode change.

Computing Restart Values. We now propose a scheme for computing restart conditions. More precisely, let t ev be an instant of mode change of System (1), meaning that one of the guards γ i changes its value. The objective is to compute consistent restart values

x + = x(t + ev ) = def lim s tev x(s) and ẋ+ = ẋ(t + ev ) = def lim s tev ẋ(s) (5)
for the index one system (2,3), given values

x -= x(t - ev ) = def lim s tev x(s) (6)
just prior to entering the new mode. Note that, since x -are variable values in the previous mode, they will not, in general, be consistent for the new mode, that is the second equation of System (2) might not be satisfied by x -. Still, these equations must be satisfied by the (yet unknown) x + . Since x(t) may be discontinuous at t ev , the derivative of x(t) may be a Dirac impulse.

To derive our method for computing the restart values, we first restrict the class of systems with the following additional assumption. We then discuss the general case.

Assumption 2 Assume DAE System (2) has the following special structure:

0 = A(x s , t) ẋ + b(x, t) (= f d ( ẋ, x, t)) 0 = f c (x, t) (7)
where x s collects the smooth elements of x, that is those being continuous and of bounded variation around the mode change. 7 Other elements of x might be discontinuous. Furthermore, we assume that x, A(. . . ), b(. . . ) are continuous functions of their arguments and remain bounded around the instant of mode change.

Note that Assumption 2 does not forbid that the triple defining the dynamics (7) actually varies with the mode i, that is has the form (A i , b i , f c,i ). Since elements of x may be discontinuous, ẋ may have Dirac impulses. In [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF] it is therefore assumed that ẋ having potentially Dirac impulses appear linearly and the linear factors are continuous functions without discontinuities. Under Assumption 2, the mathematical solution of the restart problem can be determined as follows: In a first step, since x s is smooth at the instant t ev of mode change, the matrix A(x s , t ev ) is immediately known once a change has been detected at t ev . To evaluate all the elements of x right after t ev we proceed by integrating (7) over [t ev -ε, t ev + ε] using the well known properties of Lebesgue integrals and Dirac measure. To this end, we decompose x =

x s x ns and A = A s A ns where x ns collects the nonsmooth entries of x that may experience discontinuities at the instant of mode change and A is decomposed accordingly. The following approximations hold, where the integrals

tev +ε tev -ε mean [tev -ε,tev +ε] : tev +ε tev -ε f d ( ẋ, x, t)dt = tev +ε tev -ε (A s (x s , t) ẋs + b(x, t)) dt
≈0 by footnote 7

+ tev +ε tev -ε A ns (x s , t) ẋns dt ẋns dt is a Dirac at tev ≈ A ns (x s , t ev ) (x ns (t ev + ε) -x ns (t ev -ε)) ≈ A(x s , t ev ) (x(t ev + ε) -x(t ev -ε))
In the above two approximations, the first one follows from the property of the Dirac measure and the fact that t → A ns (x s (t), t) is continuous, and the second one follows from the fact that A s (x s , t ev ) (x s (t ev + ε) -x s (t ev -ε)) ≈ 0. This leads to proposing the following scheme for computing the restart values x + from x -at the instant t ev of mode change, cf. Eqns ( 5) and ( 6):

0 = A(x s , t ev )(x + -x -) 0 = f c (x + , t ev ) (8)
Since [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF] may not hold at t - ev , and x + -x -need not to be small, the Implicit Function Theorem cannot be invoked to argue about the existence and uniqueness of a solution around x -. Since the Jacobian ∂f ∂x + of ( 8) is regular due to (3) and Assumption 2, (8) can be numerically solved with a Newton-type method. Fortunately, in many practical cases (such as a mechanical impact) the nonlinear part of f c (x, t) is often satisfied by x -and then only a linear equation system needs to be solved for x + , in which case a unique solution exists, see also section 4.3. The physical interpretation of ( 8) is that we look for restart values that are consistent for the new mode (second equation) and meet the integral of f d (. . . ) across the mode change (first equation).

In a modeling language such as Modelica or Modia, the model code has the form (2), not the special form [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF]. It is therefore desirable to compute restart values only with the form (2), without being forced to generate special code that reveals the details of the equations as in [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF]. On the basis of Assumption 1, we can nevertheless consider the following two systems of equations, where x + are the unknowns and h > 0 is small:

mixed explicit/implicit Euler implicit Euler 0 = hf d x + -x - h , x -, t ev 0 = f c (x + , t ev + h) 0 = hf d x + -x - h , x + , t ev + h 0 = f c (x + , t ev + h) (9)
Schemes ( 9) only require Assumption 1. Under Assumption 2, Schemes (9) take the form mixed explicit/implicit Euler implicit Euler

0 = A(x s , t ev )(x + -x -) +h • b(x -, t ev ) 0 = f c (x + , t ev + h) 0 = A(x s , t ev + h)(x + -x -) +h • b(x + , t ev + h) 0 = f c (x + , t ev + h) (10)
where we have used the fact that x - s ≈ x + s since x s is smooth. These schemes reduce to [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF] for h ≈ 0.

Thus, if Assumption 2 holds, the form (2) for the system model at mode i can be used to numerically compute restart values for this mode: At the instant of the mode change integrate (2) with either mixed explicit/implicit Euler or implicit Euler schemes, from x -over a small step-size h, thereby scaling the dynamic part of the model, f d (..), with h. The solution of the nonlinear equation system (9) converges to the solution of the system [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF].

Discussing the General case. Schemes (9) can be applied to any index one system. Such a brute force use without Assumption 2, however, raises questions that we review now:

-Relaxing the continuity assumption on A(. . . ), so that A(x, t) might have discontinuous elements: As a result the first contribution in the decomposition of the integral tev +ε tev f d ( ẋ, x, t)dt is no longer negligible and it is unclear, from the literature, whether a well-defined meaning to this term exists.

-Removing the linearity assumption on ẋ: The solution is no longer a Dirac and it is again unclear whether a well-defined solution exists. It is not even clear that h is the proper scaling factor for f d . Note that the linearity assumption quite often holds for physical system models, since balance equations in physics are linear in their derivatives.

To summarize, there is some evidence that precautions must be taken when relaxing Assumption 2.

Completing consistent restart. After the consistent restart values x + have been computed, the corresponding consistent restart values of ẋ+ can be determined with (4) and the first equation of (2), by solving the following non-linear equation system in the unknowns ẋ+ (in case Assumption 2 holds, this is a linear equation system with a regular Jacobian, so a unique solution exists):

f d ( ẋ+ , x + , t ev ) = 0 ∂f c ∂x (x + , t ev ) ẋ+ + ∂f c ∂t (x + , t ev ) = 0 (11)
Addressing initialization. Note that the combined use of ( 9) and ( 11) provides a method for consistent initialization: The initial mode i and guesses for x - i,init must be provided by the modeler. Afterwards, consistent values ( ẋ+ i,init , x + i,init ) for the initialization are computed with [START_REF] Brenan | Numerical Solution of Initial Value Problems in Differential-Algebraic Equations[END_REF][START_REF] Elmqvist | Systems Modeling and Programming in a Unified Environment based on Julia[END_REF].

A class of multi-mode multi-body systems satisfying Assumptions 1 and 2

In this section we exhibit a practically useful class of systems satisfying Assumptions 1 and 2. We consider a multi-body system whose model at each mode i has the following structure:

q = v (m 1 ) M (q, t) v + G T i (q, t)λ i = h(q, v, t) (m 2 ) 0 = g i (q, t) (m 3 ) (12)
where M = M T is positive definite, and G i = ∂g i ∂q has full row rank. [START_REF] Elmqvist | Modelica extensions for multi-mode DAE systems[END_REF] Eq. ( 12) describes a multi-body system with generalized position coordinates q, generalized velocity coordinates v and generalized constraint forces λ i due to the constraints (m 3 ) of the ith mode. Both the constraints (m 3 ) and the constraint forces G T i (q, t)λ i can vary with the mode i. In particular the constraints can also be completely removed. M , however, remains invariant through the different modes. Whenever a set of constraints is changing at a new mode, impulses might occur. In case an impulse is due to an impact, it is assumed that the impact is completely inelastic. We now show that Eq. ( 12) can be put to a form satisfying Assumptions 1 and 2.

Equation ( 12) is first transformed to an index two DAE with the method of Gear, Gupta and Leimkuhler [START_REF] Gear | Automatic integration of eulerlagrange equations with constraints[END_REF]:

q = v -G T i (q, t)µ i (m a 1 ) M (q, t) v + G T i (q, t)λ i = h(q, v, t) (m 2 ) 0 = g i (q, t) (m 3 ) 0 = G i (q, t)v + g (1) i (q, t) ( ṁ3 ) (14) with g (1) i (q, t) = ∂g i ∂t . (15)
In Eq. ( 14) the constraint equation (m 3 ) is differentiated once and new auxiliary variables µ i are introduced so that the number of equations and unknowns remain the same. ( 14) has the same solution as [START_REF] Elmqvist | Innovations for Future Modelica[END_REF] because µ i = 0 holds as we explain now: Inserting (m a 1 ) in ( ṁ3 ) yields:

0 = G i ( q -G T i µ i ) + G (1) i (16)
Subtracting the derivative of (m 3 ) from ( 16) results in equation 0 = G i G T i µ i . Since G i has full row rank according to assumption (13), G i G T i is regular and therefore µ i = 0.

As proposed by Gear [START_REF] Gear | Differential-algebraic equation index transformations[END_REF], by using the substitutions (λ int,i is the integral of λ i and µ int,i is the integral of µ i ):

λint,i = λ i , μint,i = µ i (17)
the DAE ( 14) which is index 2 in the variables q, v, λ i , µ i is transformed to the following DAE which is index 1 in the variables q, v, λ int,i , µ int,i :

q = v -G T i (q, t) μint,i (m b 1 ) M (q, t) v + G T i (q, t) λint,i = h(q, v, t) (m b 2 ) 0 = g i (q, t) (m 3 ) 0 = G i (q, t)v + g (1)
i (q, t) ( ṁ3 ) [START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF] (18) has equation structure [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF], with

x =     q v λ int,i µ int,i     A = I 0 0 G T i (q, t) 0 M (q, t) G T i (q, t) 0 b = -v -h(q, v, t) f c = g i (q, t) G i (q, t)v + g (1)
i (q, t) .

The Jacobian (3) of ( 18) is (P is a permutation matrix to exchange (m 3 ) and ( ṁ3 ) in order that the regularity of the Jacobian is at once visible):

J =   ∂f d ∂ ẋ ∂f c ∂x   = P     I 0 0 G T i 0 M G T i 0 0 G i 0 0 G i 0 0 0     (20)
As required by Assumption 1, this Jacobian is regular because M is positive definite and G i has full row rank. With reference to Assumption 2, it remains to show that it is legitimate to take x s = q, which amounts to requiring that q is continuous and has bounded variation around the instant of mode change. Before an impact occurs, the normal distance d j to a constraint surface j is defined as d j = g j (q, t) and d j > 0 when the multi-body system is not in contact with the constraint surface. Contact occurs when d j = 0 and therefore the constraint g j (q, t) = 0 at the time instant of the impact and q is continuous. In the fully inelastic impact case, q has bounded variation because the multibody system remains in contact with the contact surface after the impact and it is assumed that g j (q, t) is smooth. The scheme above can be easily extended to elastic impacts. In such a case, q has bounded variation provided instants of mode changes do not form a Zeno sequence.

Example: Ideal Clutch With Motor

We will show in this section how the example of Fig. 3 can be simulated with the method developed in the previous sections. The modular description with Modia is mapped to a set of about 25 equations that are first pre-processed (such as performing alias elimination). To avoid overloading the development with unnecessary details, we only show the results after this pre-processing is performed (these equations can also be easily derived manually from the circuit). In the following model, C, k, J, R are constants; u 0 (t), γ(t) are given time functions; γ = T if the clutch is released and F if it is engaged:

R(i 1 + i 2 ) + u = u 0 (sum of voltages in left loop) C u = i 1 (capacitor) kω 1 = u (emf) τ 1 = ki 2 (emf) J 1 ω1 = τ 1 -τ 2 (inertia1) J 2 ω2 = τ 2 (inertia2) if γ then τ 2 else ω 1 -ω 2 = 0 (idealized clutch) (21)
All the following analysis could be performed with [START_REF] Modelica | A Unified Object-Oriented Language for Systems Modeling[END_REF]. To concentrate on the essential details, the equations are further simplified by using the following relationships from ( 21), as well as the substitution τint,2 = τ 2

τ 1 := J 1 ω1 + τ 2 i 2 := τ 1 /k = (J 1 ω1 + τint,2 )/k i 1 := C u ( 22
)
and the equation system ( 21) can be simplified to the following four equations:

J 1 ω1 + τint,2 + kC u -k(u 0 -u)/R = 0 J 2 ω2 -τint,2 = 0 if γ then τint,2 else ω 1 -ω 2 = 0 kω 1 -u = 0 (23)
With x = [u ; ω 1 ; ω 2 ; τ int,2 ] T , ( 23) is in both modes an index one DAE of the form [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF]:

γ = t (mode 1) : A 1 =   kC J 1 0 +1 0 0 J 2 -1 0 0 0 +1   b 1 =   -k(u 0 -u)/R 0 0   f c,1 = kω 1 -u γ = f (mode 2) : A 2 = kC J 1 0 +1 0 0 J 2 -1 b 2 = -k(u 0 -u)/R 0 f c,2 = kω 1 -u ω 1 -ω 2 (24)
With ( 8), or alternatively with ( 9) and h ≈ 0, restart values can be computed:

γ = f → t (assuming x -satisfies mode 2 with A 2 , b 2 , f c,2 ) : ω + 1 = ω - 1 , ω + 2 = ω - 2 , u + = u - γ = t → f (assuming x -satisfies mode 1 with A 1 , b 1 , f c,1 ) : ω + 1 = J 1 ω - 1 + J 2 ω - 2 + kCu - J 1 + J 2 + k 2 C , ω + 2 = ω + 1 , u + = kω + 1 ( 25
)
When the clutch is disengaging, the variables are continuous at the mode change.

When the clutch is engaging, the variables are discontinuous at the mode change and Dirac impulses occur.

Implementation of Multi-mode Features in Modia

The implementation of Modia is described in [START_REF] Elmqvist | Innovations for Future Modelica[END_REF]. An extension of Modia to support a restricted class of multi-mode DAE systems along the lines developed in this section is currently under development. Since completely different sets of equations can be present in different modes, the implementation is based on just-in-time symbolic transformations and code generation of residue functions.

A dictionary from a Boolean vector of mode flags to functions is used as a cache to avoid symbolic transformations and code generation in case the same modes have been active before and a corresponding function is already available to calculate the residues. The final values of all variables from one simulation is extracted when a mode change event happens. These are inputs to one very short implicit Euler step with the residue function for the newly enabled modes in order to correctly simulate possible impulses. The new values of the state vector is used to start a new simulation until the next mode change.

Structural Analysis of Multi-Mode DAE Systems

In Section 4 we proposed an approach to analyze and simulate multi-mode DAE systems based on a generalization of DAE theory.

In this section we propose an alternative approach, more computer science oriented and detailed in [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF] (as well as its companion technical report [START_REF] Benveniste | Structural Analysis of Multi-Mode DAE Systems[END_REF]), which works for general multi-mode systems and uses a small number of principles. The key ideas are as follows.

1. Mode changes may result in discontinuous jumps and, therefore, resets must be performed in discrete computation steps. Hence, we first map the original system to discrete time by using a first order Euler scheme. This brings to discrete time both the reset actions and the dynamics within each mode, hence the principles of index analysis uniformly apply, albeit in discrete time. Also, a new principle of mode causality is invoked. 2. Mapping the dynamics to discrete time results in approximations. No approximation, however, results if we interpret the Euler scheme via nonstandard analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF] by using an infinitesimal time step. The analysis is then performed over the nonstandard reals. A final standardization step is applied to recover an effective numerical scheme.

To keep the exposure simple, we develop this on the example of ideal clutch with motor, see Section 4.4. More precisely, we consider again [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF] where, for indeed satisfied and can thus be seen as a fact for both modes. The same holds for (e 5 ) in the engaged mode. In turn these equations cannot be used, when determining the leading variables from the state variables. In this case, for the two systems (26), we have 4 variables but only 3 and 2 equations, respectively: one cannot determine the leading variables as functions of the state variables by using System [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]. Since the considered models are time-invariant, every solution has also satisfy the equations obtained by shifting forward any equation of the model. Shifting forward equations that do not bring variables that are shifted more than originally present in the system, yields so-called latent equations. For our two models, we add the latent equations, highlighted in blue:

(e1) : 0 = J1(ω • 1 -ω1) + kC(u • -u) 0 = -δ.(τ -k(u0 -u)/R) (e2) : 0 = J2(ω • 2 -ω2) -δ.τ (e3) : 0 = kω1 -u (e • 3 ) : 0 = kω • 1 -u • (e4) : 0 = τ (e5)                   
Eq. ( 29)-released and

                   (e1) : 0 = J1(ω • 1 -ω1) + kC(u • -u) 0 = -δ.(τ -k(u0 -u)/R) (e2) : 0 = J2(ω • 2 -ω2) -δ.τ (e3) : 0 = kω1 -u (e • 3 ) : 0 = kω • 1 -u • (e5) : 0 = ω1 -ω2 (e • 5 ) : 0 = ω • 1 -ω • 2 Eq. ( 29
)-engaged (29)
The associated incidence graphs are augmented accordingly and we remove the consistency equations:

e 1 --ω • 1 , u • , τ e 2 --ω • 2 , τ e • 3 --ω • 1 , u • e 4 --τ        Eq. (29)-released and        e 1 --ω • 1 , u • , τ e 2 --ω • 2 , τ e • 3 --ω • 1 , u • e • 5 --ω • 1 , ω • 2 Eq. ( 29
)-engaged (30)
In the two incidence graphs of (30), we show in red a pairing function e, i.e., a bijection, from the set of variables, to the set of equations, such that (e(x), x) is a edge of G for every leading variables x. This pairing function defines an orientation for G as follows: for each edge (e(x), x) ∈ G, we set e(x) → x and y → e(x) for every y = x such that (e(x), y) is a edge of G. The minimal cycles of the resulting directed graph form blocks of equations to be solved for their assigned variables. The blocks are partially ordered by the directed graph.

According to [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF], if a pairing function exists for its incidence graph, the system of equations possesses, in a generic sense, a unique solution for its leading variables, assuming consistent values for the state variables. "Generic" here means that the statement holds outside some exceptional numerical values when the non-zero coefficients of the Jacobian of this system of equations vary in a neighborhood. The above reasoning applies in continuous time where the forward shift is substituted back using differentiation. As far as structural analysis is concerned, we can freely exchange continuous and discrete time using the correspondence ω ↔ ω • -ω δ .

(e • 3 , e

Case ω 1 = ω 2 does not follow from previous time step. This arises if the system is engaged at the current instant t, but was released at the immediate previous time step, t -δ, i.e., t is an instant of mode change γ : t → f. This yields a new situation, not seen in Section 5.1. The engaged mode requires the consistency equation ω 1 = ω 2 , whereas ω 1 (t) and ω 2 (t) were both evaluated at the previous time step t -δ, at which ω • 1 (t -δ) = ω • 2 (t -δ) was not enforced. As a consequence, equation (e 5 ) of System ( 31) is enabled. Unfortunately, (e 5 ) possesses no dependent variable and the values of the state variables ω 1 and ω 2 were set at previous time step t -δ, with no guarantee that ω 1 (t) = ω 2 (t) would result. System [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF] is thus overdetermined at time t. A first idea would be to reject this kind of model. This would be unfortunate as our original model ( 26) was natural for our electromechanical system. To overcome this issue, we invoke the following Principle 1 (Mode Causality Principle) A guard must be evaluated before the equation it controls.

This principle leads to shifting forward equation (e 5 ) in model (29-engaged). We expect this modification to be very mild since it consists in delaying the satisfaction of the constraint by a small amount of time. Performing this yields:

assuming : (e 3 ) : 0 = kω 1 -u if γ = f do (29-engaged \ {(e 3 ), (e 5 )}) (34) Systems ( 33) and (34) possess identical right hand sides but were obtained by different reasoning-the fact that identical right hand sides were obtained is incidental to this example. System [START_REF] Trenn | Distributional Differential Algebraic Equations[END_REF][START_REF] Zimmer | Equation-Based Modeling of Variable-Structure Systems[END_REF]34) replaces the original System [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF]. As a final remark, observe that the same analysis would work without changes if the guard γ was a predicate in the state variables u, ω 1 , ω 2 .

The corresponding complete execution scheme is shown in Fig. 5. In this figure, boxes are the states of the execution scheme and their content specifies the configuration of guard, variables, and equations. For the guard and the variables:

v (resp. v) means v = t (resp. v = f).
For equations, e (resp. e) means that e is active (resp. disabled) and e means that the body of e is assumed from previous time step. Not mentioning a variable or an equation in a box means that this variable is not evaluated yet and this equation is not solved yet; for shifted equations added by the algorithm, however, we mention them underlined. to be satisfied. LE(e • ) indicates that we add e • as a latent equation. Blue (resp. black) transitions belong to a continuous-time (resp. discrete-time) dynamics. The red transition is impulsive. A semicolon is the sequential composition of computations, and the + sign denotes enabled blocks of equations, ready to be solved. The following comments are in order.

1. Observe first the parallel between the models sitting in the boxes of the diagram of Fig. 5 on the one hand, and the mixed explicit/implicit scheme (9) on the other hand. For this comparison, variables with superscript " + " in (9) correspond to shifted variables in Fig. 5 and variables with superscript " -" in (9) correspond to non shifted variables in Fig. 5. 2. Our development relies on a small set of principles:

-We map the continuous time multi-mode System (26) to discrete time System (31) by mapping derivatives ẋ to their explicit Euler scheme

x • -x δ .
-Our massaging of the equations only depends on the values taken by the guards, and the assumption regarding the satisfaction/violation of the consistency equations (here (e 3 ), and (e 5 ) for the engaged mode).

Continuous-time dynamics are colored blue; non-impulsive (resp. impulsive) resets are colored black (resp. red). The dynamics in each mode are 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 γ, variable τ is impulsive, and its value is not computed-it is set to NaN (Not a Number). Let us explain how the code of Fig. 6 is deduced from the execution scheme of Fig. 5.

Nonstandard Analysis First, we select the step size δ to be infinitesimal in the sense of nonstandard analysis [START_REF] Robinson | Nonstandard Analysis[END_REF][START_REF] Lindstrøm | An invitation to nonstandard analysis[END_REF], of which we informally recall the background we need.

In nonstandard analysis, the set R of real numbers is extended with infinite numbers, which are larger than any real number, and infinitesimal numbers, which are smaller (in absolute value) than any nonzero real number. The resulting set R is an extension of R that keeps its basic properties. In particular, it is an ordered field. We will be writing x ≈ y if x -y is an infinitesimal. Any finite element z ∈ R has a standard part, defined as the unique real number st(z) such that z ≈ st(z). Any element of R is called standard.

If x(t) is a standard differentiable function of time, its derivative ẋ(t) satisfies the property that, for any infinitesimal nonzero time step δ,

ẋ(t) ≈ x(t + δ) -x(t) δ . (35)
That is, taking an explicit first order Euler scheme with infinitesimal step size yields an exact match for an ODE, up to an infinitesimal error. Formally, one says that this Euler scheme standardizes as the solution of the ODE.

DAE model For Each Mode Each (blue) mode in Fig. 6 corresponds to a blue cycle in Fig. 5. For instance, the equation

(e 2 ) : 0 = J 2 (ω • 2 -ω 2 ) -δ.τ becomes (e 2 ) : 0 = J 2 ω • 2 -ω2 δ -τ ,
which, by Eq. ( 35), standardizes as the ODE (e 2 ) : 0 = J 2 ω2 -τ . The reader is referred to [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF] for details about standardization.

Computing The Reset Values As in the development of Section 4, one key contribution of our approach is the reset code for the mode transitions. Let us now explain how this part of Fig. 6 is derived from Fig. 5. Here the reasoning is different since we do not target a continuous time dynamics, but rather a finite sequence of discrete time steps implementing the reset actions. Let's focus on the transition γ : t → f shown in red in Fig. 6, which originates from the path sitting at the top, from the right-to the left blue cycle in Fig. 5. The corresponding dynamics is (34) and we target a discrete time dynamics involving the forward shift • and no differentiation. Hence, in Eq. (34), equations (e • 3 ), (e • 5 ) have the convenient form. In turn, one can no longer interpret (e 1 ) or (e 2 ) as differential equations. We must rather regard them as difference equations. Since they involve the infinitesimal parameter δ, standardization must be performed with care. Indeed, due to the equation (e • 5 ) of (34) and since (e 5 ) is not assumed, the velocities ω i experience a discrete jump. By (e 1 ) and (e 2 ) and since δ is infinitesimal, we infer that τ must be impulsive. This propagates throughout the different equations of the system. For impulses, the rescaling τ = def δ.τ of the system variables is applied, guided by the equations and the incidence graph (30-right) of the system, see the impulse analysis developed in [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF]. Using this rescaling, the block of equations of (34) becomes (e 1 ) : 0

= J 1 (ω • 1 -ω 1 ) + kC(u • -u) -τ + δ.k(u 0 -u)/R (e 2 ) : 0 = J 2 (ω • 2 -ω 2 ) -τ (e • 3 ) : 0 = kω • 1 -u • (e • 5 ) : 0 = ω • 1 -ω • 2 (36)
The term shown in blue involves state variables and is multiplied by the infinitesimal δ. Zeroing this blue term in System (36) leaves us with a structurally regular system that we can solve for its dependent variables τ , u

+ i = def ω • 1 = ω • 2 .
We recover in particular the formulas (25) from Section 4.4.

Constructive Semantics

The essential step in getting the final code of Fig. 6 was the construction of the state machine of Fig. 5. This state machine is called the Constructive Semantics of the original System [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF]. The notion of constructive semantics was first introduced in the context of reactive synchronous programming languages [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF][START_REF] Benveniste | Compositionality in dataflow synchronous languages: Specification and distributed code generation[END_REF][START_REF] Benveniste | The synchronous languages 12 years later[END_REF], where it played an important role in grounding compilation on solid mathematical foundations. Essentially, a constructive semantics for a discrete time dynamical system consists of:

1. A specification of the set of atomic actions, which are effective, non interruptible, state transformation operations. Executing an atomic action is often referred to as performing a micro-step; 2. A specification of the correct scheduling of the set of micro-steps constituting a reaction, by which discrete time progresses, from the current instant to the next one.

The effect of atomic actions is to propagate knowledge regarding the statuses (not evaluated, evaluated ) and values of variables. The computation of the constructive semantics of a synchronous program may succeed, and then the execution code is generated. Or it may fail, and then the program is rejected. The decision success/failure is formally sound, see [START_REF] Berry | Constructive semantics of Esterel: From theory to practice[END_REF][START_REF] Benveniste | Compositionality in dataflow synchronous languages: Specification and distributed code generation[END_REF][START_REF] Benveniste | The synchronous languages 12 years later[END_REF]. For synchronous languages, atomic actions are restricted to either (i) the evaluation of a single expression, or (ii) control flow operations.

In contrast, the class of atomic actions for multi-mode difference algebraic equations systems comprises: (i) Evaluating a guard; (ii) Solving a block of numerical equations; (iii) Equation management operations, such as shifting an equation, adding a latent equation, or changing the status involved/not involved of an equation at a given mode.

In [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF] we develop in detail the constructive semantics for multi-mode DAE systems. This allows us to formally define which model can be compiled and then simulated, and which cannot. Models that are under-or overdetermined are rejected. In addition, models which are not handled by the Principle 1 of mode causality, are rejected as well. This approach is implemented in the SunDAE proof of concept tool.

The standardization of the constructive semantics remains open in part although its main principles have been clarified. The Standardization, however, requires symbolic processing related to computer algebra.

6 Challenges in DAE Based Modeling Languages DAE based modeling languages are essential to the design of CPS. The development of such languages raise a number of challenges.

The correct simulation of mode changes and the need for resetting state variables is a first-mostly open-difficulty, particularly when impulses occur. In this paper, we proposed two different approaches for addressing this issue. The first approach relies on a transformation of the system model to a special index one form, followed by the application of a new formula for resetting the state variables. Arguments supporting this formula were given under additional assumptions on the model. The second approach builds on the use of nonstandard analysis, combined with the heritage of synchronous languages in computer science, particularly on the concept of constructive semantics. The classes of accepted/rejected models are well defined and simulation code is always produced for accepted models. In turn, the physical interpretation is understood only for restricted classes of models.

A particular difficulty of DAE based modeling languages is the need for sophisticated symbolic preprocessing of the model, called structural analysis. The Pantelides algorithm for computing the latent equations of a (single-mode) DAE model gave birth to a large body of literature since 1988. Our paper shows that structural analyses are also essential to handle mode changes.

Structural analyses play also a central role in scaling-up to huge models involving millions of equations. The community of High Performance Computing has provided important related contributions, for single-mode DAE systems. Yet, modularity in compilation and simulation remains open.

Overall, we see as a grand challenge the development of a DAE based modeling language and tool with the following features: DAE models with a very large number of modes are supported; accepted/rejected models are formally characterized; huge models are supported and can be handled in a modular way. We see the new language Modia and the SunDAE library as good starting points for tackling such challenges.

Fig. 1 .

 1 Fig. 1. DAE (left) vs. ODE (right) based modeling.

  inertia1.flange.w = inertia2.flange.w inertia1.flange.tau + inertia2.flange.tau = 0

3 .Fig. 3 .

 33 Fig. 3. Schematics of MotorAndLoad model.

Fig. 4 .

 4 Fig. 4. Simulation results of the MotorAndLoad model.

• 5 )

 5 is(29-engaged), in which (e 3 , e 5 ) are taken as assumptions: assuming :(e 3 ) : 0 = kω 1 -u (e 5 ) : 0 = ω 1 -ω 2 if γ = f do (29-engaged \ {(e 3 ), (e 5 )})

Fig. 5 .

 5 Fig. 5. Structural analysis of System (31): state machine describing the execution scheme of one time step of System (31).

  The transitions of the state machine indicate the actions performed when moving to the next state. FS(e) indicates that e is shifted. PR(e) indicates that e is known

	ω1, ω2, u, τ, ω • 1 , ω • 2 , u • , γ, e1, e2, e3, e4, e5, e • 3 , e • 5	ω1, ω2, u, γ, e3, e4, e5, e • 3 , e • 5	ω1, ω2, u, γ, e3, e4
	ω1, ω2, u, e3, e5	start	ω1, ω2	ω1, ω2, u, e3
		ω1, ω2, u,
		γ, e3, e4, e5
	ω1, ω2, u,
	γ, e3, e4, e5, e3

  • , ω • 1 , ω • 2 . Solving (36) yields in particular the reset values for the state variables u + = def u • , and ω

See the extensive literature available in https://www.modelica.org/publications.

A DAE f ( ẋ, x, t) = 0 has differentiation index n if one or more equations must be differentiated n-times until the equations can be algebraically transformed to an ODE form ẋ = f (x, t), see for example[START_REF] Campbell | The index of general nonlinear DAEs[END_REF].

A function f : R → R is said to have bounded variation if it is the primitive of a Lebesgue integrable function[START_REF] Dunford | Linear Operators, Part I, General Theory[END_REF]. As a consequence, t+h t ḟ (s)ds → 0 when h → 0.

Acknowledgements: The authors want to thank Toivo Henningsson, Lund, Sweden for the collaboration regarding Julia and the design of Modia, as well as the implementation of instantiation and flattening.

simplicity, we substitute τint,2 by τ :

We first analyze separately the model for each mode of the clutch.

Separate analysis of each mode, in discrete time

We begin by providing the model for each mode. We highlight in blue the equations that are unique to the considered mode. Other equations are shared. In the "released" mode, the two shafts are independent. In the "engaged" mode, the velocities of the two shafts are algebraically related.

Following key idea 1, we replace derivatives by their first order explicit Euler scheme, in discrete time with constant step size δ > 0. Let

be the forward time shift operator by an amount of δ. System (26) expands as:

System (26)-engaged [START_REF] Robinson | Nonstandard Analysis[END_REF] Let us first focus on System (26)-released. The state variables are u, ω 1 , ω 2 and their respective values are known initially. The current step must determine the values of the leading variables τ, u

Towards this, we first form the incidence graph G of each system, which is a nondirected bipartite graph having as vertices: the four leading variables, plus the two systems of equations {(e 1 ), (e 2 ), (e 3 ), (e 4 )} and {(e 1 ), (e 2 ), (e 3 ), (e 5 )}. An edge from an equation to a leading variable exists if and only if this variable is involved in that equation. Incidence graphs for models (26-released) and (26engaged) are:

Observe that, for both models, equation (e 3 ) involves no leading variable: it is a consistency equation, i.e., a constraint that must be satisfied as a result of the execution of previous time steps. Once initialization is performed, (e 3 ) is

Global Discrete-Time Analysis

In this section, we handle each of the two modes as well as the mode changes in a uniform way. Let's consider now the discretized version of our two modes System [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]:

Notice how model [START_REF] Thoma | Introduction to Bond Graphs and Their Applications[END_REF] encompasses both the released and engaged modes, as well as the mode changes, in the same uniform setting of discrete time systems. This calls for using similar reasoning to determine the new values of the leading variables regardless of whether the system is evolving in continuous mode or is at an event of mode change. We now detail the how and when the different involved equations are used.

1. Using (e 1 ), . . . , (e 3 ): Those three equations are active in all modes. Equation (e 3 ) is useless, so we are left with 2 equations and 4 leading variables and no subset of them can be evaluated using the 3 available equations only. We thus have to evaluate the guard γ in order to know which equation among (e 4 ) or (e 5 ) is active. We successively analyze the two cases below.

2-released. Case γ = t: Equation (e 4 ) is enabled and equation (e 5 ) is disabled. The reasoning proceeds exactly as for getting the model (29-released). The difference is that we take the consistency equation (e 3 ) as an assumption, since its satisfaction results from the execution of the previous time step. The resulting model is thus: assuming : (e 3 ) : 0 = kω 1 -u if γ = t do (29-released \ {(e 3 )}) [START_REF] Trenn | Distributional Differential Algebraic Equations[END_REF] where (29-released \ {(e 3 )}) means that equation (e 3 ) is removed from system 29-released. Note that model (29-released) applies both within the "released" mode and at the instant of mode change γ : f → t.

2-engaged. Case γ = f: Equation (e 4 ) is disabled and equation (e 5 ) is enabled. With reference to system (29-engaged), an important difference occurs: assuming that the consistency equation (e 5 ) results from having executed the previous time step is not always valid (equation (e 5 ) is guarded and can thus be disabled). Consequently the following two sub-cases need to be considered:

Case ω 1 = ω 2 follows from previous time step. This corresponds to the case in which the system was already in mode "engaged" at the previous time step. The separate analysis developed for the "engaged" mode in Section 5.1 applies with no change. The model, as augmented with the two latent equations Otherwise, we make no distinction between instants of mode changes and other instants: our treatment is uniform. -Our massaging of the equations has two objectives: finding latent equations if needed, and shifting forward equations when required by the principle of mode causality (Principle 1). This set of principles is small, clean, and powerful enough to encompass general systems with a structural analysis alike the one we developed here for the clutch example. See [START_REF] Benveniste | Structural analysis of multi-mode DAE systems[END_REF] for a presentation of this approach for general multi-mode DAE systems.

Effective Simulation Code mode

Fig. 6. Actual simulation code for the clutch.

So far the execution scheme of Fig. 5 is not satisfactory, as it uses an explicit first order Euler scheme. From the execution scheme of Fig. 5, however, the actual simulation code shown on Fig. 6 is derived. Observe that it consists of a DAE model of index zero or one for each mode, and code for resetting the state variables at mode changes.

The model for each mode can be simulated, e.g., using a BDF method as proposed in Section 4.2. Regarding the resets at mode changes, ω - i is the previous value of state variable ω i , which is the left limit of ω i when exiting a mode. Similarly, ω + i is the reset value for state variable ω i when entering the new mode.


domain: spi phys

text:
Fabrizio Di Donfrancesco
email: fabrizio.didonfrancesco@onera.fr

Antoine Placzek
email: antoine.placzek@onera.fr

Christophe Blondeau
email: christophe.blondeau@onera.fr

Jean-Camille Chassaing
email: jean-camille.chassaing@sorbonne-universite.fr

REDUCED ORDER TIME SPECTRAL MODEL FOR THE ALE FORMULATION OF THE COMPRESSIBLE NAVIER-STOKES EQUATIONS

Keywords: Reduced Order Models, Proper Orthogonal Decomposition, Time Spectral Method, Discrete Empirical Interpolation Method

In this paper we address the construction of a time and a frequency domain Reduced Order Model (ROM) for the ALE formulation of the Navier-Stokes equations. The wide variety of fluid dynamics ROMs reported in literature shares the aim of reducing the dimensionality of dynamical systems by performing a projection of the governing equations onto a basis that may be often constructed via the Proper Orthogonal Decomposition (POD). Unfortunately, the applicability of this kind of ROMs in the time domain still remains problematic as non-linear effects induce numerical stability problems, especially for aeroelastic applications. Nevertheless, in many cases, only the periodic steady state is of concern. In this context, spectral methods like the Time Spectral Method (TSM) avoid stability problems by expressing the conservative variables as Fourier series in time with spatially varying coefficients to be computed. In this work, we develop a POD based ROM which aims to reduce the computational cost of the TSM. Numerical tests have been carried out in order to highlight the potentiality of the proposed technique.

INTRODUCTION

Dynamical systems are the basic framework for modelling and control of a large variety of complex systems of scientific interest or industrial value. The aeronautic and aerospace fields are important examples in which some complex physical phenomena, such as fluid-structure interaction, shock interaction, flow separation, limit-cycle oscillations and so on, can be investigated with the aforementioned dynamical systems. Interactions between fluid and moving (deforming) structures are among the most important issues of aircraft design and numerical simulation has been one of the few available means for studying this kind of problems. However, the growing need for an improved accuracy attainable with a high fidelity (or full order) model leads inevitably to an unaffordable computational cost whose reduction represents one of the main motivation of this work. The projection-based reduced order models (ROMs) are physics based methods that can potentially yield very high speedups. These approaches involve IFASD-2019-134 the resolution of the dynamic equations of the system by projecting them on a suitably chosen low-dimensional subspace. The effective dimension reducibility for these methods is usually limited to problems with linear or multi-linear terms. In fact, when a general non-linearity is present, the cost to evaluate the projected non-linear function still depends on the dimension of the original system, resulting in a ROM simulation time of the same order as the original system. This is known as the lifting bottleneck and various approaches exist to tackle this issue. Masked projection approaches, such as DEIM [START_REF] Chaturantabut | Nonlinear model Reduction via Discrete Empirical Interpolation[END_REF], QDEIM [START_REF] Drmac | A New Selection Operator for the Discrete Empirical Interpolation Method -improved a priori error bound and extensions[END_REF], MPE [START_REF] Astrid | Missing Point Estimation in Models Described by Proper Orthogonal Decomposition[END_REF] approximate a given general non-linear term by interpolating (or least square fitting) it on the basis of a certain number of judiciously selected mesh points. Various approaches can be distinguished by the different greedy selection algorithms even though they all minimize the representation error associated to the non-linear terms. Considering aeroelastic dynamical systems solved on a deforming mesh, additional difficulties arise since, on one side, the mesh evolution must be taken into account in the definition of the POD basis and the non-linear operator and, on the other side, the involved highly non-linear responses severely limit the use of approximation techniques like masked projections. The first point has been discussed by Antonnen et al. [START_REF] Anttonen | POD-Based reduced-order models with deforming grids[END_REF] who introduced an index based formulation for the POD. Among the several studies entailed by this POD formulation, the one of Freno et al. [START_REF] Freno | A proper orthogonal decomposition method for nonlinear flows with deforming meshes[END_REF] addresses subsonic and transonic flows around airfoils subject to forced oscillations. They approximated the flow using basis functions which, although not explicitly functions of time, depend on parameters associated with flow unsteadiness. However, the hyper reduction techniques and the long time stability for this kind of flows are not addressed which still remains a limiting factor. Nevertheless, for the simulation of time periodic solutions for internal or external flows only the periodic steady state is of concern. Spectral-methods exploit the time periodicity nature of the flow by expressing the conservative variables as Fourier series in time with spatially varying coefficients. Spectral computations can be performed in either the time domain, the frequency domain, or a combination of the two. Several formulations have been proposed in the literature. The time-spectral method (TSM) introduced by Gopinath and Jameson [START_REF] Gopinath | Time Spectral Method for Periodic Unsteady Computations over Two-and Three-Dimensional Bodies[END_REF][START_REF] Gopinath | Application of the Time Spectral Method to Periodic Unsteady Vortex Shedding[END_REF] is similar to the harmonic balance method of Hall et al. [START_REF] Hall | Computation of Unsteady Nonlinear Flows in Cascades Using a Harmonic Balance Technique[END_REF]. However, the time-spectral method is solved completely in the time domain. This paper presents a reduced order formulation for the TSM (RO-TSM). First, a projectionbased ROM is solved in the time domain. We show with a simple-test case that this kind of ROM fails when long term stability or hyper-reduction is envisaged. Then, the TSM is implemented as presented in [START_REF] Blondeau | A Modular Implementation Of The Time Spectral Method For Aeroelastic Analysis And Optimization On Structured Meshes[END_REF] and represents the full order model (FOM) to reduce. Provided a POD basis for the TSM solution, we propose a non-linear reduction via least-squares minimization which allows to reduce the number of unknowns and shows globally better efficiency and convergence properties. Numerical test cases are investigated to highlight the potentiality of the proposed formulation.

THEORETICAL FRAMEWORK

Full Order Model

The governing equations considered in the present paper for an aeroelastic CFD problem are represented by the Arbitrary Lagrangian-Eulerian (ALE) formulation of the Navier-Stokes equations. We consider the general problem of a laminar or turbulent compressible flow in an entrained and possibly deformable control volume Ω(t) of boundary Σ(t). In the semi-discrete form the ALE formulation for the basic cell of the computational domain reads:

d dt Ω(t) W dΩ = - 6 i=1 i (t) F c (W , s)•ndΣ- 6 i=1 i (t) F d (W )•ndΣ = -R Ω (W , t) (1)
IFASD-2019-134

where W is the vector of the conservatives variables ρ, ρU and ρE. F c and F d are, respectively, the convective and diffusive fluxes and Σ i represents the i-th face of the hexahedral cell considered in the structured mesh. The arbitrary motion of the computational mesh is taken into account by the definition of c = Us as the convective velocity, where s is the mesh velocity with respect the spatial domain. The flow equations are solved by a finite-volume method with ONERA in-house software elsA [START_REF] Cambier | The Onera elsA CFD software: Input from research and feedback from industry[END_REF] (ONERA-Airbus-Safran property). Once spatially discretized, we can write:

dV(Ω(t))W Ω dt + R Ω (W Ω , t) = 0 (2)
where W Ω is the average of the conservative variables W in the cell Ω, V(Ω(t)) is the volume of the related cell and R Ω is a non-linear residual operator that computes the flux balance for each cell. For a constant cell volume, we obtain the system of equations:

V dW Ω dt + R Ω (W Ω , t) = 0 (3)
that represents the FOM. For sake of clarity, from now on, the term W Ω is referred to as W .

Projection based ROM in time domain

The first step of construction of the ROM consists in approximating the conservative field by a base solution W 0 (steady, time-average flow or the initial snapshot) and a linear combination of solution modes Φ i :

W (t) ≈ W 0 + Φ a(t) (4)
In this work, the solution basis is derived by retaining a small number of the orthogonal modes computed via the Proper Orthogonal Decomposition (POD). Sirovich [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, Parts I-III[END_REF] introduced the socalled "method of snapshots" for computing a POD basis. Assuming that the snapshots are centered with respect to a base solution W 0 (steady or time-average flow) and collected in a matrix [ W ] ∈ R Nx×Nt with N x the number of degrees of freedom and N t the number of collected snapshots, the POD is equivalent to a singular value decomposition (SVD) of the matrix [ W ], so that:

[ W ] = [Φ][Σ][V] T (5)
where the matrix [Φ] ∈ R Nx×Nr is an orthonormal matrix that contains the POD mode vectors, with N r = rank([ W ]) ≤ min(N x , N t ). The diagonal matrix [Σ] ∈ R Nr×Nr contains the singular values of [ W ] listed in order of decreasing magnitude. If we define the diagonal element of [Σ] as σ i the rate of "energy" captured by the first N m modes is given by

E Nm = Nm i=1 σ 2 i / Nr i=1 σ 2 i .
This quantity is an indicator of the energy neglected by retaining only the first N m POD modes. Then, by substituting eq.( 4) in eq.( 2) and using the Galerkin projection the following system of N m ODEs is obtained:

da dt = Φ T R Ω (W 0 + Φ a(t)) (6)
It should be noted that, for the compressible Navier-Stokes equations, the operator R Ω is nonlinear and cannot be expressed explicitly in terms of the coordinates a(t) unless specific approximations or formulations are introduced to obtain a multi-linear form. The computation of the entire non-linear term, at each time step, would lead to a computational time for the ROM of the same order of magnitude of the FOM affecting the performance of the reduced order model.

To tackle this issue, masked projection approaches aim to approximate the non-linear term as:

R Ω (W , t) ≈ Ψ c(t) (7)
IFASD-2019-134

where, analogously to eq.( 5), a POD of the non-linear term snapshots is performed to obtain the basis Ψ of size N x × N l , where N l is the number of retained modes after truncating the POD basis of the non-linear term. The coefficients c j (t) represent the related amplitudes. The system of eq.( 7) is overdetermined because it involves N x equations and N l unknowns. At this stage a so-called mask matrix is introduced to select a certain number N f ≥ N l and N f N x of lines of eq.( 7) so as to obtain a smaller system to solve. More specifically, if e j is the j-th column of the identity operator of size N x × N x , the mask-matrix is the column-orthogonal matrix P = [e 1 , . . . , e N f ] ∈ R Nx×N f so that the following masked subsystem is obtained:

P T R Ω (W , t) ≈ (P T Ψ) c(t) (8)
In particular, when N f = N l the matrix P T Ψ of size N l × N f can be inverted provided that it is positive definite, otherwise a least-square fitting is required. The lines selected by applying the masked projection correspond to points (or cells) in the discretized domain so that the evaluation of the non-linear term at the interpolation points depends only on a limited number of neighboring grid points (or cells) depending on the order of the adopted space discretization scheme. Different masked projection approaches (DEIM [START_REF] Chaturantabut | Nonlinear model Reduction via Discrete Empirical Interpolation[END_REF], QDEIM [START_REF] Drmac | A New Selection Operator for the Discrete Empirical Interpolation Method -improved a priori error bound and extensions[END_REF] , MPE [START_REF] Astrid | Missing Point Estimation in Models Described by Proper Orthogonal Decomposition[END_REF]) provide different greedy algorithms to compute the mask-matrix P . It is worth noting that the additional approximation resulting from masked projection accumulates an error over each time step of the time integration affecting the long term stability of this kind of ROMs.

Time Spectral Method

Provided that W is periodic with period T = 2π/ω, the Fourier series of eq.( 3) is obtained:

∞ k=-∞ (ikωV Ŵ k + Rk ) exp(ikωt) = 0 (9)
where Ŵ k and Rk are the Fourier coefficients of W and R Ω corresponding to mode k. The complex exponential family forming an orthogonal basis, the only way for eq.( 9) to be true is that the weight of every mode k is zero. The time spectral method (TSM) [START_REF] Sicot | Block-Jacobi Implicit Algorithms for the Time Spectral Method[END_REF] uses a Inverse Discrete Fourier Transform (IDFT) to cast back in the time domain a set of 2N + 1 equations from eq.( 9). The IDFT induces linear relations between Fouriers coefficients Ŵ k and a uniform sampling of W within the period:

W n = N k=-N Ŵ k exp(iωn∆t) , 0 ≤ n ≤ 2N + 1 (10)
with W n = W (n∆t) and ∆t = T /(2N + 1). This leads to a time discretization with a new time operator D t as follows:

R Ω (W n ) + VD t W n = 0 , 0 ≤ n ≤ 2N + 1 (11)
These steady equations correspond to 2N + 1 instants equally spaced within the period. The new time operator connects all time levels and can be expressed analytically as

D t W n = N i=-N d i W n+i (12)
with

d i = π T (-1) i+1 csc πi 2N +1 , i = 0 0 , i = 0 (13)
IFASD-2019-134

The resolution of eq.( 11) can be challenging and pseudo-time stepping techniques are used to time march the equations to the steady-state solutions. In this work, explicit and fully implicit TSM implementations approaches are used. The reader interested to the mathematical derivation of the these approaches is referred to [START_REF] Blondeau | A Modular Implementation Of The Time Spectral Method For Aeroelastic Analysis And Optimization On Structured Meshes[END_REF].

Reduced Order Time Spectral Method

The eq.( 11) can be challenging to solve because of stability issues which are more penalizing with an increasing number of harmonics N , that is when the diagonal dominance of the system is reduced [START_REF] Sicot | Block-Jacobi Implicit Algorithms for the Time Spectral Method[END_REF]. In this section we introduce a Reduced Order Time Spectral Method (RO-TSM) which implies a lower computational cost and exhibits better convergence properties than the full order TSM (FO-TSM).

For the construction of the RO-TSM, analogously to eq.( 5), a POD of the snapshots W n is performed without any snapshot matrix centering in order to obtain the basis Φ ∈ R Nx×(2N +1) . The obtained POD modes span a reduced subspace that contains the TSM solutions. In such a case, no energy truncation of the POD basis is done. Then the solution is approximated as:

W n ≈ Φ a n , 0 ≤ n < 2N + 1 (14)
by substituting eq.( 14) in eq.( 11) the following coupled system of (2N + 1) × N x steady equations involving (2N + 1) × (2N + 1) unknowns is obtained:

R = R Ω (Φ a n ) + VD t Φ a n = 0 , 0 ≤ n < 2N + 1 (15)
A steady-state solution is typically sought using a pseudo-time stepping technique. This approach is a common choice when starting from a poor initial guess (such as uniform flow). However, when a good starting guess is available it can be more efficient to apply the Newton's method directly. Applying Newton's method to solve the fully coupled system [START_REF] Xie | Data-Driven Filtered Reduced Order Modeling of Fluid Flows[END_REF] with all instance amplitudes gathered in vector a = (a 0 , a 1 , . . . , a 2N +1 ) results in the following iterations. For k = 1 . . . K

A k ∆a k = -R(a k ) (16)
with K determined by satisfaction of a convergence criterion. The resulting Jacobian matrix is defined as

A k = J k Φ D =       J f k 0 Vd 1 I . . . Vd -1 I Vd -1 I J f k 1 . . . . . . . . . . . . . . . Vd 1 I Vd 1 I . . . Vd -1 I J f k (2N )            Φ 0 . . . 0 0 Φ . . . . . . . . . . . . . . . 0 0 . . . 0 Φ      (17)
where the Jacobian J f j = ∂R(a j )/∂(Φa j ). The system of equations ( 16) can be solved in a least-square sense by using the Gauss-Newton algorithm:

∆a k = ∆a k ∈ R (2N +1)×(2N +1) : A k ∆a k + R(a) k 2 2 is minimized (18)
that is:

∆a k = -(Φ T D J k T J k Φ D ) -1 (Φ T D J k T ) R(a) k (19)
It can be noted that the previous resolution of the Gauss-Newton algorithm is equivalent to the resolution of the Petrov-Galerkin projection of the system in eq.( 16) onto the subspace L = (J k Φ D ).

IFASD-2019-134

The presented method is beneficial for two important aspects. First, the solution is sought in a subspace that is a priori suitable as a result of the POD. Second, the matrices involved in the numerical resolution have smaller dimensions. If we consider the Jacobian matrix A k , for the FO-TSM it involves (2N + 1)N x columns, whereas it reduces to only (2N + 1)(2N + 1) for the RO-TSM.

Definition of the initial guess

As mentioned above, the success of the presented RO-TSM depends on the quality of the starting guess a 0 . In this work, the FO-TSM initial guess is set as the uniform flow at the free-stream conditions W f s . If we define W n as the exact solution of system [START_REF] Xie | Data-Driven Filtered Reduced Order Modeling of Fluid Flows[END_REF], the quality of the RO-TSM initial guess is represented by the error:

e = W n -W f s = e + e ⊥ = W n -Φ(Φ T W f s ) e -(I -ΦΦ T )W f s e ⊥ (20)
where e ∈ span(Φ) and e ⊥ ∈ span(Φ ⊥ ). By construction of the POD basis, the solution W n ∈ span(Φ). As a consequence, if we define the initial guess as the projection of the snapshot W f s onto the POD basis Φ

a 0 n = Φ T W f s (21)
the contribution of e ⊥ is eliminated from the first iteration onwards, resulting in a more suitable initial guess for the Newton algorithm. The mentioned procedure is graphically depicted in Fig. 1.

NUMERICAL EXPERIMENTS

The projection based ROM resulting from the formulation of eq.( 6) is used to model the inviscid flow around an oscillating NACA0012 airfoil under a prescribed harmonic pitch motion. The motion is sinusoidal with an angle of attack defined by the function

α(t) = α 0 + α sin(ωt) (22)
IFASD-2019-134

where α 0 = 0.016 • is the mean angle of attack and α = 2.51 • is the amplitude of the pitch motion. The flow parameters are those of the reference CT5 [START_REF] Landon | NACA 0012 Oscillatory and Transient Pitching[END_REF] except for the Mach number set to M = 0.5. Under these conditions the reduced frequency is κ = ωc/2U ∞ = 0.125.

Projection based ROM in the time domain

A finite volume solver for the Euler equations in used as FOM. We employ a sampling of 1000 snapshots over a cycle of oscillation to construct the POD basis for the flow solution and for the non-linear residual term. The energy distribution of the solution POD modes is plotted in Fig. 2. The value of E Nm = Nm i=1 σ i / Nr i=1 σ i for N m = 10 is already over 0.9999. The first 2 density POD modes are depicted in Fig. 3. The ROM is integrated using the DEIM and the QDEIM as masked projection techniques. The test case is apparently not challenging, considering the moderate reduced frequency κ and the regular structures of the POD modes. Nevertheless, using from 5 to 40 non-linear POD modes and the DEIM/QDEIM as masked projection, the ROM has always diverged before the end of the first cycle. Then, the ROM is integrated without using any kind of masked projection.

In such a case, the ROM is integrated successfully for the oscillation cycle but the efficiency is compromised. However, the ROM always diverges when the long term stability (beyond the time sampling interval) is investigated. In Fig. 4 the evolution of the amplitudes related to first 3 (and most energetic) modes is shown. It can be seen that after the first oscillation cycle the ROM diverges after few iterations.

Figure 4: Comparison of the first 3 modal coordinates a i computed as the solution of the projection based ROM and of the reference modal coordinates for the test case defined in reference CT5 [START_REF] Landon | NACA 0012 Oscillatory and Transient Pitching[END_REF].

IFASD-2019-134

Discussion

Freno et al. [START_REF] Freno | A proper orthogonal decomposition method for nonlinear flows with deforming meshes[END_REF] already investigated this kind of ROMs by using dynamic average and dynamic basis functions to model compressible flow using a deforming mesh. These dynamic functions vary continuously with respect to parameters associated with the flow unsteadiness and mesh deformation, and they are optimal, subject to the prescribed form. However, during the integration of the ROM the entire non-linear term is projected without providing any type of hyper-reduction. This would naturally lead to more stable ROMs, but it still implies a highly costly ROM. Also, it is well known that this kind of ROM run into problems of convergence when long-term stability is investigated [START_REF] Sirisup | A spectral viscosity method for correcting the long-term behavior of POD models[END_REF][START_REF] Xie | Data-Driven Filtered Reduced Order Modeling of Fluid Flows[END_REF]. The main issue is related to the approximation error that accumulates over the time integration leading to divergence. To overcome this problem, different solutions like the addition of artificial viscosity or the use of data-driven closure terms have been investigated. Anyway, the stability of fully non-linear ROMs is still an open question.

Classical unsteady techniques can be avoided when considering flows periodic in time. Methods like the TSM seek the periodic state of the solution by approximating the solution as Fourier series in time. These methods are intrinsically immune to the long-term stability problems mentioned above.

RO-TSM

To illustrates the potentialities of the proposed RO-TSM a more challenging test case is studied.

We have selected the test case noted as dynamic index 55 reported by the AGARD group in [START_REF] Davis | NACA 64A010 (NASA Ames model) oscillatory pitching[END_REF]. It corresponds to an oscillating NACA64A010 airfoil under the prescribed pitch motion of eq.( 22), where α 0 = -0.22 • and α = 1.01 • . The Mach number is 0.796 and the excitation frequency 34.4 Hz. Under these conditions the reduced frequency is κ = ωc/2U ∞ = 0.41.

In this case, the FOM is represented by the TSM. The reference solution is obtained by a finite volume solver for the Euler equations. We remind that the number N of harmonics of eq.( 10) and so the number (2N + 1) of coupled steady equations to solve is a parameter of the TSM solver. In Fig. 5 it can be noted that to reproduce correctly the lift coefficient (on the left) 1 harmonic is sufficient. However this choice implies a significant error on the drag coefficient prediction (on the right). Finally, the non-linear effects involved in this test case require a minimum number of harmonics equal to 3 in order to correctly predict the aerodynamic coefficients.

Figure 5: Comparison of the aerodynamic force coefficients between the finite volume solver and the TSM solver.

IFASD-2019-134

1 harmonic test case

Firstly, the RO-TSM is constructed by taking into account only 1 harmonic for the TSM resolution. For this test case an explicit TSM (with a CFL = 1) is used. Then, as already presented in Section 2.4, a POD basis is constructed using the 3 snapshots corresponding to the 3 time instances of the TSM. The 3 density POD modes are depicted in Fig. 6.

POM 1 POM 2 POM 3 In Fig. 7 the rate of convergence of the full order and the reduced order TSM is shown. Firstly, it is obvious that the RO-TSM presents a much faster and more regular rate of convergence (note the different abscissa). The zoomed part of the plot emphasizes the differences between the initial guesses. Indeed, as explained in Section 2.4.1, the projection of the free-stream state vector onto the POD basis provides a different and more suitable initial guess for the RO-TSM. The RO-TSM requires only about 10 iterations to attain a residual comparable to the one reached by the TSM in 1000 iterations. Moreover, each iteration of the RO-TSM required a smaller system to be solved resulting, globally, in a lower computational cost. We remind also that the POD basis is built from the TSM snapshots associated to a given residual norm. This norm is the convergence criterion that monitors the cost of the TSM. Since the RO-TSM seeks a solution in the vector space spanned by the TSM snapshots, it is not possible to reach a better numerical accuracy than the one the TSM simulation. Finally, in Fig. 8 the absolute difference between the density fields related to the RO-TSM and the TSM is shown as a proof of the high accuracy attained by the presented model.

3 harmonics

As previously introduced, the non-linear effects of the presented test case require at least 3 harmonics for the TSM computation to reach an acceptable accuracy on the aerodynamic drag coefficient prediction. The choice of 3 harmonics implies 7 steady coupled equations to be solved which affects the diagonal dominance of the TSM eq.( 11) and therefore the global complexity of the TSM resolution. This is shown by the convergence issues occurring when an explicit TSM solver is employed. Then, a fully implicit TSM solver (with a block-Jacobi implementation [START_REF] Blondeau | A Modular Implementation Of The Time Spectral Method For Aeroelastic Analysis And Optimization On Structured Meshes[END_REF]) is employed with a CFL number for the pseudo-time integration equal to 2.5.

Once the FO-TSM computation performed, the POD basis is constructed with the resulting 7 solution field snapshots. The density POD modes are depicted in Fig. 9. there are no significant differences in the stability of the model between the cases with 1 and the 3 harmonics. This is not the case for the FO-TSM, which has implied implicitation and a low CFL condition. Instead, provided a POD basis for the involved time instances, the same RO-TSM reaches a monotonic convergence. In Fig. 11 the absolute difference between the density

CONCLUSION

This paper presented the development of a Reduced Order Time Spectral Model approach. This method was motivated by the stability shortcomings related to the time dependent projectionbased ROMs when nonlinear flows are investigated. In many aeroelastic cases, periodicity is a characteristic of the flow. In this context, spectral methods represent a suitable alternative to the classical unsteady techniques. We formulated a RO-TSM which presents two main features. First, the solution is sought onto a subspace that is a priori suitable as a result of the POD. Second, the computational complexity of the resolution is globally reduced. This method is validated for the inviscid flow around an oscillating NACA64A010 airfoil under prescribed pitch motion at a reduced frequency κ = ωc/2U ∞ = 0.41. Even when considering 3 harmonics for the Fourier series approximation, the model converges monotonically and shows very accurate results.

Future works will extend the present method for a higher number of harmonics for different and more complex test cases. Then this RO-TSM will be extended to take into account parametric changes. In particular, to model an off-reference condition it will be necessary to interpolate the POD basis. To this end, interpolation techniques like the interpolation on a tangent space to a Grassmann manifold will be investigated [START_REF] Farhat | An Interpolation Method for the Adaptation of Reduced-Order Models to Parameter Changes and Its Application to Aeroelasticty[END_REF].

Figure 1 :

 1 Figure 1: A graphical description of the initialization of the RO-TSM by projecting onto the POD subspace.

2 iFigure 2 :POM 2 Figure 3 :

 2223 Figure 2: Proper orthogonal eigenvalue spectrum.

Figure 6 :

 6 Figure 6: Density fields of the 3 proper orthogonal modes (POMs) composing the POD basis for the RO-TSM for N = 1.

Figure 7 :Figure 8 :

 78 Figure 7: Comparison of the Euclidean norm of the density residual for the full order and the reduced order TSM.Note the different abscissa related to the two approaches.

Figure 9 :

 9 Figure 9: Density fields of the 7 proper orthogonal modes (POMs) composing the POD bases for the RO-TSM for N = 3.

Figure 10 :

 10 Figure 10: The Euclidean norm of the density residual for the 7 steady state fields of the reduced order TSM.

Figure 11 :

 11 Figure 11: On the first line, the density fields at the seven time instances of the TSM. On the second line, the absolute difference between the RO-TSM and the TSM density fields.

COPYRIGHT STATEMENT

The authors confirm that they, and/or their company or organization, hold copyright on all of the original material included in this paper. The authors also confirm that they have obtained permission, from the copyright holder of any third party material included in this paper, to publish it as part of their paper. The authors confirm that they give permission, or have obtained permission from the copyright holder of this paper, for the publication and distribution of this paper as part of the IFASD-2019 proceedings or as individual off-prints from the proceedings.