# Structural Methods for handling mode changes in multimode DAE systems

Albert Benveniste, Benoit Caillaud, Yahao Chen, Khalil Ghorbal, Mathias Malandain

► To cite this version:

Albert Benveniste, Benoit Caillaud, Yahao Chen, Khalil Ghorbal, Mathias Malandain. Structural Methods for handling mode changes in multimode DAE systems. RR-9603, Inria. 2026, pp.61. hal-05385651v2

HAL Id: hal-05385651

https://inria.hal.science/hal-05385651v2

Submitted on 8 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Inria

Structural Methods for handling
mode changes in multimode
DAE systems

Albert Benveniste, Benoit Caillaud, Yahao Chen, Khalil Ghorbal, Mathias Malandain

RESEARCH
REPORT

Nº 9603

December 2025

Project-Team Hycomes

ISSN 0249-6399 ISRN INRIA/RR-9603-FR+ENG
<figure><img src="image_2.png" /></figure>
Inria

Structural Methods for handling mode
changes in multimode DAE systems

Albert Benveniste*, Benoit Caillaud, Yahao Chen, Khalil Ghorbal, Mathias Malandain

Project-Team Hycomes

Research Report n° 9603 — version 2† — initial version December 2025 — revised version January 2026 — 55 pages

RESEARCH CENTRE
Centre Inria de l'Université de Rennes

Campus universitaire de Beaulieu
Avenue du Général Leclerc
35042 Rennes Cedex

<footnote>* All authors are with Inria center at Rennes University, Campus de Beaulieu, 35042 Rennes cedex, France; email: surname.name@inria.fr</footnote>

<footnote>† The handling of cascaded mode changes was added.</footnote>
**Abstract:** Hybrid systems are an important concept in Cyber-Physical Systems modeling, for which multiphysics modeling from first principles and the reuse of models from libraries are key. To achieve this, DAEs must be used to specify the dynamics in each discrete state (or *mode* in our context). This led to the development of DAE-based equational languages supporting multiple modes, of which Modelica is a popular standard. Mode switching can be time- or state-based. Impulsive behaviors can occur at mode changes. While mode changes are well understood in particular physics (e.g., contact mechanics), this is not the case in physics-agnostic paradigms such as Modelica. This situation causes difficulties for the compilation of programs, often requiring users to manually "smooth out" mode changes. In this paper, we propose a novel approach for the hot restart at mode changes in such paradigms. We propose a mathematical meaning for hot restarts (such a mathematical meaning does not exist in general), as well as a combined *structural-and-impulse analysis* for mode changes, generating the hot restart even in the presence of impulses. Our algorithm detects at compile time if the mode change is insufficiently specified, in which case it returns diagnostics information to the user.

**About this revision:** In a first version of this report, we restricted ourselves to the handling of mode changes separating *long* modes (lasting for a positive duration). In this revised version, we extend our approach to handling cascades of *transient* modes (of zero duration) such as, e.g., occurring in complex impacts in multibody mechanics.

**Key-words:** structural analysis, differential-algebraic equations (DAE), multi-mode systems, switched system
# Analyse Structurelle des changements de mode
dans les DAE multimodes

**Résumé** : La modélisation des systèmes cyber-physiques repose sur une modélisation à partir des principes de la physique, et en réutilisant au maximum des modèles prédéfinis issus d'une bibliothèque. Cela exige le recours aux Equations Différentielles Algébriques (DAE) admettant plusieurs modes (une DAE commutée, ou une DAE hybride). Le standard de modélisation est le langage Modelica. Les changements de mode peuvent être déclenchés de manière externe, ou par des conditions portant sur les états. Ces changements de mode sont connus et traités à l'intérieur de physiques particulières (mécanique avec contacts). Il en va autrement dans un cadre multi-physique général, qui est, pourtant, celui de Modelica et d'autres langages de modélisation multi-physique. Dans ce papier, nous proposons une approche nouvelle pour le redémarrage à chaud suite à un changement de mode. Noter qu'il n'existe pas de définition mathématique de ce qu'est une solution dans notre cadre général. Notre méthode utilise une analyse structurelle doublée d'un calcul symbolique des comportements impulsifs. Notre méthode s'applique lors de la phase de compilation et permet de détecter, avant toute simulation, si le modèle soumis est éventuellement insuffisamment spécifié.

**Concernant cette version** : Dans une première version de ce rapport, nous traitions seulement de changements de mode séparant deux modes *longs* (de durée strictement positive). Dane cette version-ci, **nous étendons notre approche aux cascades of modes transitoires** (de durée nulle) telle qu'on en rencontre en mécanique multi-corps.

**Mots-clés** : analyse structurelle, équations algébro-différentielles (DAE), systèmes multi-mode, changement de mode à chaud
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>A cup-and-ball example</td><td>5</td></tr><tr><td>3</td><td>Background</td><td>10</td></tr><tr><td>3.1</td><td>Structural analysis of systems of equations</td><td>11</td></tr><tr><td>3.2</td><td>The Σ-method for DAEs</td><td>12</td></tr><tr><td>3.3</td><td>Implicit function theorem</td><td>13</td></tr><tr><td>4</td><td>Problem setting: hot restart of mode changes</td><td>13</td></tr><tr><td>5</td><td>Mode change arrays</td><td>14</td></tr><tr><td>5.1</td><td>Shifting and differentiating</td><td>14</td></tr><tr><td>5.2</td><td>Mode change arrays and facts</td><td>16</td></tr><tr><td>5.3</td><td>Structural Implicit Function Theorem</td><td>18</td></tr><tr><td>6</td><td>Rescaling analysis of mode change arrays</td><td>19</td></tr><tr><td>7</td><td>Main theorems</td><td>22</td></tr><tr><td>7.1</td><td>Correctness of our approach</td><td>22</td></tr><tr><td>7.2</td><td>Solving the rescaling problem</td><td>23</td></tr><tr><td>7.3</td><td>Characterizing the solutions</td><td>24</td></tr><tr><td>8</td><td>Proof of correctness</td><td>26</td></tr><tr><td>9</td><td>Addendum: solving nondeterminism and handling cascaded mode changes</td><td>29</td></tr><tr><td>9.1</td><td>Revisiting hot restart</td><td>29</td></tr><tr><td>9.2</td><td>Revisiting mode change arrays</td><td>30</td></tr><tr><td>9.3</td><td>Reusing our previous analysis</td><td>31</td></tr><tr><td>9.4</td><td>Cup-and-ball with elastic impact</td><td>31</td></tr><tr><td>9.5</td><td>On the modeling power of our approach</td><td>35</td></tr><tr><td>10</td><td>Conclusion and perspectives</td><td>36</td></tr><tr><td>A</td><td>Proof of Corollary 1</td><td>41</td></tr><tr><td>B</td><td>Proof of Lemma 8</td><td>41</td></tr><tr><td>C</td><td>Proof of Theorem 2</td><td>41</td></tr><tr><td>D</td><td>Examples having no good solution</td><td>45</td></tr><tr><td>D.1</td><td>The cup-and-ball with K>1</td><td>45</td></tr><tr><td>D.2</td><td>The cup-and-ball with exogenous mode change</td><td>47</td></tr><tr><td>D.3</td><td>A “strange” cup-and-ball</td><td>49</td></tr><tr><td>E</td><td>Nonlinear systems: a clutch</td><td>50</td></tr><tr><td>F</td><td>Comparing with Trenn et al.</td><td>52</td></tr><tr><td>F.1</td><td>A linear DAE system</td><td>52</td></tr><tr><td>F.2</td><td>Further investigating the clutch</td><td>53</td></tr></table>
# 1 Introduction

## Motivation and related work

In this paper, we consider multimode DAE systems, i.e., collections of finitely many (possibly nonlinear) equations of the form

$$ \text{if } b \text{ then } f(\text{the } x_i\text{'s and their derivatives}) = 0, \qquad (1) $$

where $x_i, i=1, \dots, n$ are (time-dependent) numerical variables, $f$ is a sufficiently differentiable numerical function,<fnref n="1" /> and $b$ is a Boolean expression of predicates over system variables and time. Boolean $b$ guards equation $f = 0$, meaning that this equation is enabled when $b$ evaluates to $T$ (the constant “true”), disabled otherwise. We assume that each predicate takes the form $g(x^-) \ge 0$, where $g$ is a smooth function of system variables, and $x^-(t) = \lim_{s \to t} x(s)$ is the left-limit of $x$ at instant $t$. Function $g$ is often called a zero-crossing function: its successive switches from negative to positive values trigger the onsets of the predicate.

Multimode DAE systems constitute a generic form for the hybrid systems specified by Modelica [30, 26], a popular DAE-based modeling language dedicated to the modeling of multiphysical systems.<fnref n="2" /> Modelica is by itself physics-agnostic, which allows the user to include in the model a specification of the software code controlling the system.<fnref n="3" /> It is important to establish such modeling languages on solid mathematical grounds.

Multimode DAE systems, however, exhibit significant difficulties. First, there is no reference mathematical definition of the solution of a multimode DAE in general, but only for specific subclasses of physics, e.g., contact mechanics [8]. Second, modes can be numerous—the number of modes typically grows exponentially with the number of subsystems. Due to these difficulties, existing tools do not implement mode-aware compilation: some models, although clearly valid, fail to get correctly simulated. Such pathological models are by no means exceptional, nor are they difficult to exhibit [6]. This work addresses the first difficulty by:

1. giving a mathematically sound definition of mode change events in multimode DAEs, and

2. proposing an effective algorithm computing the hot restart after a mode change.

**Related work:** An extensive review of the literature on multimode DAEs is available in [2]; in this brief discussion, we only collect its major findings. For selected physics, e.g., multi-body systems with contacts and electrical circuits with idealized switches and diodes, dedicated methods are proposed to handle possible impulses [33, 32, 36, 21, 1]. It is not clear if the special methods used in these areas extend to general models. Mehrmann et al. [20] propose numerical techniques to handle chattering between modes. In Zimmer's PhD thesis [41], multimode DAEs are considered with varying structure and index; however, impulsive behaviors are not supported. Both references assume that consistent reset values are explicitly given for each mode: this assumption is not suited to a compositional framework where one wants to assemble predefined physical components. Elmqvist et al. [17, 27] propose a high level description of multimode models as an extension to the Modelica 3.3 state machines. However, mode changes with impulsive behavior are not supported and not all types of multimode systems can be handled, as mentioned above.

In our previous work [2], we proposed for the first time an alternative approach for handling multimode DAEs. A restricted class of multimode DAEs, possibly involving impulsive behaviors, was considered (since then known as semi-linear multimode DAEs [3]). For this subclass, two alternative approaches were investigated and shown to be equivalent: a Gear-Gupta-Leimkuhler

<footnote n="1">¹What “sufficiently” means will be made precise in Section 3.2.</footnote>

<footnote n="2">²Modelica is the most well-known instance of such languages; Amesim and Simscape are other languages used in the industry.</footnote>

<footnote n="3">³This is in contrast to physics-oriented modeling methodologies such as bond graphs [37] or port-Hamiltonian modeling [40].</footnote>
method [18], implemented in the Julia package Modia, and a novel approach relying on nonstandard analysis [25]. The latter approach is the root of the one presented in this paper.

**Zooming in:** Key references for our work are the important contributions by Stephan Trenn and coworkers [38, 39, 24, 22, 13]. The pioneering work [38] investigates the inconsistent DAE initialization problem: a regular linear DAE $Ex'=Ax+f$ is considered for $t \ge 0$, with an inconsistent initialization $x_{(-\infty,0)}$, meaning that the left limit $\lim_{s \to 0} x(s)$ does not satisfy the consistency constraints imposed by the DAE. This possibly causes impulsive behavior for some variables at instant 0. The classical theory of distributions over the timeline $\mathbb{R}$ is inadequate to handle this problem, as the value of a distribution at a given instant is undefined in general. The use of piecewise-smooth distributions was proposed in [38, Chapter 2] to give an explicit solution for the inconsistent initialization of a linear DAE, by which the behavior of the impulsive variables is described as a linear combination of Dirac derivatives. A key pillar in computing this solution is the quasi-Weierstrass form for a linear DAE, which consists in constructing a state basis and recombining equations, such that the impulsive order (differentiation degree of the Dirac function $\delta$) of each variable is well identified, and pure constraints are separated from the ODE part.

An interesting step forward is [22], in which the inconsistent initialization problem for DAE $Ex'=Ax+g(x)+f$ is investigated, where $g$ is a smooth nonlinear function. The authors assume that the nonlinear part $g(x)$ remains “foreign to impulsions” (see [22] for a formalization)—the reason being that, for $x$ impulsive and $g$ nonlinear, $g(x)$ is undefined in general ($g(\delta)$ being an example). Under this assumption, an appropriate time-invariant change of coordinates and equations allows the authors to statically decompose the DAE system into three subsystems: 1) a nonlinear ODE, 2) a nonlinear static constraint, and 3) a linear DAE carrying the impulsive part of the system. Using this decomposition and reusing the background from linear DAEs, the authors propose an explicit solution to the inconsistent initialization problem. These results comply with the linear control systems vision: states are internal and are only defined up to a change of basis.

## Objectives and approach

Our work differs from the works of Trenn *et al.* in a number of aspects.

First, a key difference is our overall objective, namely: *to provide mathematical soundness for the compilation of DAE-based modeling languages such as Modelica*. In Modelica models, parameters occurring in models are entered right before launching a simulation, i.e., after the simulation code was compiled. Therefore, compilation involves the structure of the model, not the actual values of its parameters. In particular, a key transformation of a DAE model is the *structural analysis*, mainly consisting of the *index reduction*, which transforms the DAE into an ODE-like system by suitably differentiating selected equations. The technique used in Modelica tools to perform index reduction uses the model structure, not the numerical details of the model. Our work complies with this philosophy:

1. *Changes of state basis are prohibited*^4—such changes are little relevant in the context of physical modeling (unlike in black-box modeling).

2. Whereas in classical linear systems theory, regularity (of matrices or linear pencils) plays a central role, *we give up numerical regularity and replace it by structural regularity* (also called structural nonsingularity).

A square matrix is *structurally nonsingular* if all of its diagonal entries can be made nonzero by pre- and post-multiplying it by permutation matrices. A static system of nonlinear equations is *structurally nonsingular* if so is its Jacobian around its solution. Structural nonsingularity holds if and only if a one-to-one mapping exists between equations and dependent variables, i.e., a *perfect*

<footnote>^4This prohibits numerical steps such as the construction of the above mentioned quasi-Weierstrass form.</footnote>
*matching*. Structural nonsingularity is necessary, and generically sufficient,<fnref n="5" /> for exact regularity of a static system. Structural regularity is checked on the incidence graph of the system (an abstraction of its structure), relying on graph-based algorithms that can scale up much better than numerical ones: those are used in high-performance computing, and in all the compilers of DAE-based modeling tools [9]. Moreover, as structural methods do not rely on parameter values, they allow us to identify, prior to simulation, if the model is over-, under-, or well-determined, providing useful information to model designers.

The second difference is that we do not describe the full trajectories of all the variables, including impulsion events. Instead, *we detect impulsive behaviors, without identifying their exact nature,<fnref n="6" /> and rescale impulsive variables to make them non-impulsive*. Our rescaling algorithm bears similarities with J. Pryce's Σ-method for the structural analysis of (single-mode) DAEs [35].

## Contribution

Based on the above approach, our contributions are the following:

* A general definition of *hot restart* (Problem 1);

* A system of *rescaling equations* (Problem 2);

* A procedure to generate hot restarts (Procedure 1);

* A proof of correctness of Procedure 1 (Theorem 1);

* Bounds for design parameters of Problem 2 (Theorem 2);

* An effective algorithm for solving rescaling equations.

Procedure 1 generates the system of equations by which restart values for the states are determined from the values just before the mode change. Uniqueness of the restart system of equations is guaranteed, when it exists, by Theorem 1, and its non-existence typically corresponds to a lack of determinism, expressing that the model is insufficiently specified for hot restart.

In our previous work [3], the same problem was addressed, by building on top of nonstandard analysis [25]. The resulting algorithm, however, was not completely specified and was difficult to analyze. The method introduced in the present paper fixes both issues. In particular:

* It handles general nonlinear systems; mode changes can be state-based, not only time-based.<fnref n="7" /> Yet, hot restart succeeds only if impulsive variables are involved *linearly* in the system model—this restriction mirrors the one formulated in [22], where it was formulated in Definition 6 as an extra and complex “condition ($G_p$)”. In our case, there is no need for a separate check: our algorithm discovers by itself if this condition is satisfied or violated.

* It is physics-agnostic, yet is able to reveal hidden physical invariants: in the cup-and-ball example of Section 2, the hot restart we generate preserves angular momentum, although no such law was explicitly stated in the model.

* By being graph-based, it has the potential to scale up to very large systems, unlike previous approaches.

This report is organized as follows. Section 2 introduces our approach by means of an illustrative example. This example is simple; yet, state-of-the-art DAE-based modeling tools fail to simulate it correctly—they actually crash when reaching mode change events. Background material is recalled

<footnote n="5">⁵Generically means that the matrix remains nonsingular almost everywhere when its non-zero entries vary over some neighborhood. See [5] for a short tutorial on structural methods.</footnote>

<footnote n="6">⁶That is, without describing them as a suitable linear combination of Dirac measures. Recall that the occurrence of $g(\delta)$, for $g$ nonlinear, prevents us from providing a full definition of solutions in this case.</footnote>

<footnote n="7">⁷Consequently, all multimode DAE models are supported by our compilation method.</footnote>
in Section 3: Sections 3.1 and 3.2 focus on the structural analysis of both static systems of equations and (single-mode) DAE systems, and Section 3.3 states a slight reformulation of the implicit function theorem. The core problem addressed by this paper, namely, the hot restart problem, is stated in Section 4. Sections 5 and 6 develop our key result, namely the rescaling analysis. Main theorems regarding our approach are collected in Section 7.1, and details on the resulting algorithm are developed in Section 7.2. Finally, a mathematical characterization of the hot restarts generated by our method is presented in Section 7.3. Main proofs are collected in Section 8.

Additional proofs are collected in Appendices A, B, and C. Additional examples illustrating the features of our approach are presented in Appendices D and E. Finally, a comparison with the approach by Trenn et al. is discussed in Appendix F.

In a first version of this report, we only considered mode changes between two successive “long modes” (lasting for a positive duration with the same DAE dynamics), with no additional hot restart constraint. The case of finite cascades of successive transient modes (of zero duration) is now supported in Section 9 of this revised version, with limited objectives. This is achieved by considering additional hot restart constraints.

## 2 A cup-and-ball example

We develop a cup-and-ball game example to give an intuitive presentation of our approach; missing background is introduced informally and will be developed in Section 3.1. This example illustrates the main challenges to be addressed: mode changes are state-based, and they involve impulsive behaviors. As a matter of fact, major tools fail to simulate this model.

A ball, modeled by a point mass, is attached to one end of a rope, while the other end of the rope is fixed to the origin of the plane in the model. The ball is subject to the unilateral constraint set by the rope, but moves freely while the distance between the ball and the origin is less than its length. A model for a 2D version of this example is:

<figure><img src="image_3.png" /></figure>

$$ 
\left\{ 
\begin{array}{ll}
0 = x'' + \lambda x & (e_1) \\
0 = y'' + \lambda y + g & (e_2) \\
0 \le L^2 - (x^2+y^2) & (\kappa_1) \\
0 \le \lambda & (\kappa_2) \\
0 = [L^2 - (x^2+y^2)] \times \lambda & (\kappa_3)
\end{array}
\right.
\qquad (2) 
$$

where the dependent variables are the position $(x, y)$ of the ball in Cartesian coordinates and the rope tension $\lambda$.

Subsystem $(\kappa_1, \kappa_2, \kappa_3)$ expresses that the tension is nonnegative, the distance between the ball and the origin is less than or equal to $L$, and one cannot have a nonzero tension and a distance less than $L$ at the same time: (2) is a complementarity system, and such systems are key in non-smooth mechanics [32]. This model is not of the form (1) yet because of unilateral constraints $(\kappa_1, \kappa_2)$. Using the technique presented in [26], we redefine the graph of $(\kappa_1, \kappa_2, \kappa_3)$ as a parametric curve, represented by the following three equations:

$$ 
\left\{ 
\begin{array}{l}
\gamma = [s \le 0] \\
0 = \text{if } \gamma \text{ then } L^2 - (x^2 + y^2) \text{ else } \lambda \\
s = \text{if } \gamma \text{ then } -\lambda \text{ else } L^2 - (x^2 + y^2)
\end{array}
\right.
\qquad (3) 
$$

The model now constitutes a logico-numerical fixpoint equation with dependent variables $x'', y'', \lambda, \gamma, s$. Such equation can have zero, one, or infinitely many solutions. No characterization exists that could serve as a basis for a structural analysis. We thus decide to refuse solving such mixed logico-numerical systems.

To break the fixpoint equation defining $\gamma$, we choose to base guards on left-limits of signals. This yields the modified model (4), where the modification is highlighted in red. For convenience,
we also grouped the equations that are only active in modes $\gamma = T$ and $\gamma = F$, respectively:

$$ 
\left\{
\begin{array}{ll}
0 = x'' + \lambda x & (e_1) \\
0 = y'' + \lambda y + g & (e_2) \\
\gamma = [s' \le 0]; \gamma(0) = F & (k_0) \\
\text{if } \gamma \text{ then } 0 = L^2 - (x^2+y^2) & (k_1) \\
\text{and } 0 = \lambda + s & (k_2) \\
\text{else } 0 = \lambda & (k_3) \\
\text{and } 0 = (L^2 - (x^2+y^2)) - s & (k_4)
\end{array}
\right.
\qquad (4)
$$

where “else” is a shorthand for “if not $\gamma$ then”. Let us assume that the ball is initially under free motion ($\gamma = F$). Eventually, the rope will get straight ($\gamma : F \to T$), with two possible continuations, depending on whether the impact is assumed elastic (the ball bounces inward when the rope gets straight) or inelastic (the rope remains straight).

* If the impact is inelastic, then the radial velocity of the ball becomes zero, hence the rope gets straight ($\gamma = T$) until the gravity pushes the ball back to free motion ($\gamma = F$).

* If the impact is elastic, then the radial velocity changes its sign (the simplest model is that it gets opposite and equal in magnitude) and the ball bounces back to free motion ($\gamma = F$); in this case, the ball spends zero time in $\gamma = T$ mode, we call it *transient*.

The respective timing views are the following:

$$ 
\begin{array}{l}
\text{inelastic : } \underbrace{\text{free motion}}_{\text{long}} \text{ straight rope } \underbrace{\text{free motion}}_{\text{long}} \\
\text{elastic : } \underbrace{\text{free motion}}_{\text{long}} \text{ straight rope } \underbrace{\text{free motion}}_{\text{transient}}
\end{array}
\qquad (5)
$$

Note that subsystem $(\kappa_1, \kappa_2, \kappa_3)$ leaves the impact law at mode change insufficiently specified: the model does not state whether the impact is elastic or inelastic. We focus here on the inelastic case, where only long modes occur. An impulsive behavior is expected at the mode change from free motion to straight rope, $\gamma : F \to T$; we will illustrate our method on this mode change.

We first consider each mode separately in (4). Under the mode $\gamma = F$, free motion occurs, which yields an ODE. In contrast, under $\gamma = T$, the rope is straight, thus enabling the algebraic constraint $(k_1)$: this yields a DAE. Index reduction is performed [12, 31], which consists here in adding to the model the two *latent equations* $(k'_1, k''_1)$ obtained by successive differentiations of $(k_1)$. The resulting model splits into two parts [31]:

$$ 
\begin{array}{ll}
\text{leading equations :} & \left\{
\begin{array}{ll}
0 = x'' + \lambda x & (e_1) \\
0 = y'' + \lambda y & (e_2) \\
0 = xx'' + x'^2 + y'^2 + yy'' & (k''_1)
\end{array}
\right. \\
\text{consistency equations :} & \left\{
\begin{array}{ll}
0 = L^2 - (x^2+y^2) & (k_1) \\
0 = xx' + yy' & (k'_1)
\end{array}
\right.
\end{array}
\qquad (6)
$$

In (6) and in the sequel, we highlight in blue an injective matching of equations to variables. It is one-to-one for the leading equations, showing structural nonsingularity — see the discussion on structural nonsingularity in the introduction, and Section 3.1 for a formal presentation. Indeed, the leading equations generically determine the leading variables $x'', y'', \lambda$, while the consistency equation constraints the state variables $x, y, x', y'$ at initialization, leaving two degrees of freedom.

Focus next on mode change $\gamma : F \to T$, and let $t$ be the instant when it occurs. We wish to identify how the previous mode $\gamma = F$ influences the restart of the new mode $\gamma = T$.

**Approach 1** Our approach unfolds as follows. In a neighborhood of the mode change instant:
1. Each derivative is discretized using a simple forward Euler scheme with time step $\epsilon$, i.e., any derivative $x'$ is identified with $\epsilon^{-1}(x(t+\epsilon)-x(t))$; we call the added equation $x' = \epsilon^{-1}(x(t+\epsilon)-x(t))$ an Euler identity;

2. Using this interpretation of derivatives, we stack finitely many successive transitions of the resulting discrete-time system around the mode change instant, thus building an array of equations called *mode change array* — the last transition stacked in this array will serve to initialize the next mode. The identification and compensation of impulsive behaviors is performed on this array, leading to its *rescaling*;

3. Finally, we let $\epsilon:=0$ in the rescaled array.<fnref n="8" /> As a result, with the expansion of $x'(t)$ as $\epsilon^{-1}(x(t+\epsilon)-x(t))$, we retrieve the derivative of $x$ at $t$. The resulting *restart system* relates the latest state variables of this system, to the right limits of the trajectories of the previous mode; it is used for hot restart. □

The approximations errors introduced by enforcing the Euler identity in step 1 are handled in step 3; we prove later that, under certain conditions, the limit exists when $\epsilon \searrow 0$ and is indeed obtained by setting $\epsilon:=0$. To be properly used for the hot restart of the new mode, this restart system should satisfy the following requirements:

**Requirement 1 (for the restart system)**

1. It should be deterministic, i.e., it should uniquely determine the restart state variables of the new mode;

2. The restart state variables should satisfy the consistency conditions of the new mode.

Requirement 1.1 means that the mode change is determined. Requirement 1.2 expresses that a hot restart should at least satisfy the conditions for being a consistent start.

Performing Step 1 of Approach 1 results in a discretization of the dynamics of the clutch in a neighborhood of the mode change. We focus on the neighborhood of the mode change, and we detail Step 2. To simplify the writing, we use the following notation using the instant $t$ of mode change as a reference: for any variable $x$, let

$$ \bullet x \equiv x(t-\epsilon), x \equiv x(t) \text{ and } x \bullet \equiv x(t+\epsilon). $$

Using this notation, we display in Fig. 1 the dynamics of the new mode at instant $t$, as completed following (6). We regard it as a static system of equations $A_0$ and call it *mode change array*. Variables $x, y, x', y'$ are fully determined by the previous mode, due to the Euler identities $x = \bullet x + \epsilon \times \bullet x'$ and $x' = \bullet x' + \epsilon \times \bullet x''$ (and similarly for $y$): we call them *past variables*. The dependent variables of $A_0$ are $x'', y'', \lambda$ (past variables are excluded).

**Finding 1 (facts)** The equation ($k_1$) is highlighted in green to indicate that it is satisfied up to $O(\epsilon)$, despite it involving only past variables; the reason for this is that the mode change was detected at instant $\bullet t$ by the *zero-crossing event* $0=\bullet(x^2+y^2-L^2)$, hence, $x^2+y^2-L^2=O(\epsilon)$. We call *fact* an equation that involves past variables only and is satisfied up to an $O(\epsilon)$. Being nearly satisfied and involving no dependent variable, *facts can be ignored*.<fnref n="9" /> □

This is in contrast to the next equation ($k'_1$): it also involves only past variables but it is not satisfied. We decide to disable it at this instant, which amounts to postponing its consideration for a while. We are left with the three equations ($e_1, e_2, k''_1$), which determine the leading variables $x'', y'', \lambda$. Hence, array $A_0$ uniquely determines the leading variables from the values of past variables: Requirement 1.1 is met. In contrast, Requirement 1.2 is not met since consistency equation ($k'_1$) was violated.

<footnote n="8">⁸Throughout this paper, symbol $:==$ means the assignment of a value to a variable or parameter.</footnote>

<footnote n="9">⁹The reader is kindly advised against inferring anything about the human world from this statement.</footnote>
$$A_0: \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \gamma=T, t \\ 0 = y'' + \lambda y + g & (e_2) & \\ 0 = L^2 - (x^2+y^2) & (k_1) & \text{fact} \\ 0 = xx' + yy' & (k'_1) & \text{disabled} \\ 0 = xx'' + x'^2 + y'^2 + yy'' & (k''_1) & \end{array} \right. $$

Figure 1: Cup-and-ball example: Mode change array $A_0$. In the last column we point the facts and the disabled conflicting equations; black equations are enabled. Triple $(e_1, e_2, k''_1)$ is structurally nonsingular, with a one-to-one matching $\mathcal{M} = \{(e_1, x''), (e_2, \lambda), (k''_1, y'')\}$ between enabled equations and variables.

$$A_1: \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \text{instant } t \\ 0 = y'' + \lambda y + g & (e_2) & \\ 0 = L^2 - (x^2+y^2) & (k_1) & \text{fact} \\ 0 = xx' + yy' & (k'_1) & \text{disabled} \\ 0 = xx'' + x'^2 + y'^2 + yy'' & (k''_1) & \text{disabled} \\ 0 = x'' - \varepsilon^{-1}(x'^{\bullet} - x') & (\mathcal{E}_{x''}) & \\ 0 = y'' - \varepsilon^{-1}(y'^{\bullet} - y') & (\mathcal{E}_{y''}) & \\ 0 = (x'' + \lambda x)^{\bullet} & (e_1^{\bullet}) & \text{instant } t+\varepsilon \\ 0 = (y'' + \lambda y + g)^{\bullet} & (e_2^{\bullet}) & \\ 0 = (L^2 - (x^2+y^2))^{\bullet} & (k_1^{\bullet}) & \text{fact} \\ 0 = (xx' + yy')^{\bullet} & (k'_1^{\bullet}) & \\ 0 = (xx'' + x'^2 + y'^2 + yy'')^{\bullet} & (k''_1^{\bullet}) & \end{array} \right. $$

Figure 2: Cup-and-ball example: Mode change array $A_1$. Facts and disabled conflicting equations are pointed. Euler identities $(\mathcal{E}_{x''})$, $(\mathcal{E}_{y''})$, relating second derivatives to shifts of first derivatives, were added. The enabled subsystem (in black) is structurally nonsingular, with a one-to-one matching $\mathcal{M}$ between equations and variables, highlighted in blue.

$$A_1: \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (f_1) \\ 0 = y'' + \lambda y + g & (f_2) \\ 0 = (xx' + yy')^{\bullet} & (f_3) \\ 0 = x'' - \varepsilon^{-1}(x'^{\bullet} - x') & (f_4) \\ 0 = y'' - \varepsilon^{-1}(y'^{\bullet} - y') & (f_5) \end{array} \right. $$

Figure 3: Cup-and-ball example: The black subsystem of Fig. 2, in which leading equations of instant $t^{\bullet} = t+\varepsilon$ were omitted, and equations renumbered for convenience.

Let us extend array $A_0$ by one more instant, by considering array $A_1$, shown in Fig. 2. Euler identities $(\mathcal{E}_{x''})$ and $(\mathcal{E}_{y''})$ were added, to make explicit the relations between $x''$ and $x'^{\bullet}$, and $y''$ and $y'^{\bullet}$. Variables $x, y, x', y'$ are still determined by the previous mode. Equation $(k_1)$ remains a fact, and so does $(k'_1)$, for the same reasons. Let us disable the conflicting equations $(k'_1, k''_1)$. The remaining system is structurally nonsingular, as evidenced by the perfect matching highlighted in blue. Array $A_1$ determines the values of all its dependent variables knowing the past variables: Requirement 1.1 is met. In addition, Requirement 1.2 is now satisfied since all the consistency equations at instant $t^{\bullet}$ are enabled.

**Finding 2** We solved the hot restart of the new mode, in discretized setting, i.e. Step 2 of
Approach 1 was successful when using $A_1$.

<figure><img src="image_4.png" /></figure>

What happens when performing Step 3, i.e., when letting $\epsilon:=0$? Here, we only consider the black subsystem of Fig. 2, which yields the system shown in Fig. 3, where the same perfect matching is highlighted in blue.

**Finding 3 (trying $\epsilon:=0$)** Letting $\epsilon:=0$ in $A_1$ causes trouble, due to the occurrence of $\epsilon^{-1}$ in equations ($f_4, f_5$). Multiplying these two equations by $\epsilon$ does not solve the problem either, since letting $\epsilon:=0$ in this case erases $y''$ from equation ($f_5$), which makes $A_1$ singular. □

<figure><img src="image_5.png" /></figure>

Only one action can get rid of this difficulty: properly identifying impulsive variables and rescaling them.

**Notations 1** In the sequel, we use the convention that ($f_i$) denotes an equation of array $A_1$, whereas $f_i$ shall denote the corresponding function defined by its right-hand side. □

<figure><img src="image_6.png" /></figure>

We propose a simple approach to identify and quantify impulsive behaviors:<fnref n="10" />

1. Say that variable $x$ has rescaling offset $k$, written $\mu_x=k$, if its value in the solution of system $A_1$ is $0(\epsilon^{-k})$; (7)

2. $\mu_x=k$ results in the rescaling $x^\downarrow = \text{def } \epsilon^k \times x$.

Focus on equation ($f_5$). Since $y'$ and $y'^\bullet$ are both non impulsive, equation ($f_5$) in array $A_1$ implies $\mu_{y''}=1+\max(\mu_{y'}, \mu_{y'^\bullet})=1$. In addition, we also associate, to function $f_5$, a rescaling offset indicating the magnitude order w.r.t. $\epsilon^{-1}$ of $f_5$ in the neighborhood of a solution of this equation. That is:

$$ \mu_{f_5} = \max(\mu_{y''}, 1+\max(\mu_{y'}, \mu_{y'^\bullet})) = 1. $$

Performing this kind of reasoning for all the equations of array $A_1$ yields the following rescaling analysis:

$$ \left\{ \begin{array}{l} \mu_{x''} = \mu_{f_1} = \max(\mu_{x''}, \mu_\lambda) \\ \mu_\lambda = \mu_{f_2} = \max(\mu_{y''}, \mu_\lambda) \\ \mu_{y'^\bullet} = \mu_{f_3} = \max(\mu_{x'^\bullet}, \mu_{y'^\bullet}) \\ 1 + \mu_{x'^\bullet} = \mu_{f_4} = \max(\mu_{x''}, 1 + \mu_{x'^\bullet}) \\ \mu_{y''} = \mu_{f_5} = \max(\mu_{y''}, 1 + \mu_{y'^\bullet}) \end{array} \right. $$

We search for solutions of the rescaling analysis, satisfying the following requirement, expressing that state variables for restart should not be impulsive:

**Requirement 2 (on states)** State variables of the last instant possess rescaling offset zero.

The only solution making restart variables non-impulsive is the following:

$$ \mu_{x''} = \mu_{\lambda} = \mu_{y''} = 1 = \mu_{f_5} = \mu_{f_2} = \mu_{f_1} \qquad (8) $$

whereas other variables and functions have rescaling offset zero. Rescaling variables and equations using (7), Step 2, yields:

$$ \begin{aligned} x''^\downarrow &= \text{def } \epsilon \times x'' = x'^\bullet - x' \quad ; \quad f_1^\downarrow &= \text{def } \epsilon \times f_1 \\ y''^\downarrow &= \text{def } \epsilon \times y'' = y'^\bullet - y' \quad ; \quad f_5^\downarrow &= \text{def } \epsilon \times f_5 \\ \lambda^\downarrow &= \text{def } \epsilon \times \lambda \quad ; \quad f_2^\downarrow &= \text{def } \epsilon \times f_2 \end{aligned} \qquad (9) $$

By using (9), we have

$$ A_1^\downarrow : \left\{ \begin{array}{l} 0 = x''^\downarrow + \lambda^\downarrow x \\ 0 = y''^\downarrow + \lambda^\downarrow y \\ 0 = (xx' + yy')^\bullet \\ 0 = x''^\downarrow - (x'^\bullet - x') \\ 0 = y''^\downarrow - (y'^\bullet - y') \end{array} \right. \qquad (10) $$

<footnote n="10">¹⁰What follows is an informal definition; the formalization will be given in Section 6.</footnote>
**Finding 4 (retrying ε := 0)** At this point, letting ε := 0 in A↓₁ preserves regularity, structurally— actually, ε already disappeared from A↓₁. □

We are now ready to conclude. New positions were determined by the previous mode: $x^{\bullet} = x$ and $y^{\bullet} = y$, expressing that positions are continuous at mode change. The restart system is finally obtained by renaming, in (10):

1. the variables set by the previous mode $x, y, x', y'$, by the left-limits at mode change $x^-, y^-, x'^-, y'^-$;

2. the tail variables $x^{\bullet}, y^{\bullet}, x'^{\bullet}, y'^{\bullet}$, by the restart values for the new mode $x^{+}, y^{+}, x'^{+}, y'^{+}$; whereas

3. other variables are auxiliary and are not renamed.

The second renaming action is sensible only if the following requirement is satisfied by rescaled array (10):

**Requirement 3** Performing rescaling shall bring, in array $A_1^↓$, no variable attached to instants later than $t+\epsilon$.

This requirement is indeed satisfied by rescaled array (10). Finally, performing the above legitimate renaming yields the

$$ \text{restart system: } \left\{ \begin{array}{l} 0 = x^{+} - x^{-} \\ 0 = y^{+} - y^{-} \\ 0 = x'^{+} - x'^{-} + \lambda^{\downarrow} x^{-} \\ 0 = y'^{+} - y'^{-} + \lambda^{\downarrow} y^{-} \\ 0 = x^{-} x'^{+} + y^{-} y'^{+} \end{array} \right. \qquad (11) $$

Note that eliminating $\lambda^↓$ from (11) yields

$$ x'^{+}y^{+} - y'^{+}x^{+} = x'^{-}y^{-} - y'^{-}x^{-}, \qquad (12) $$

which expresses the preservation of angular momentum, a physical invariant that was not explicitly specified in the original modeling: our approach “discovered” it. Also note that our solution does not give a meaning to variables $x'', y'', \lambda$ beyond their status of being impulsive. Only the rescaled variable $\lambda^↓$ is well defined. This approach is formalized and generalized in Sections 4 and 5; before that, we shall proceed with background material.

## 3 Background

In this section, we recall our background, consisting of the structural analysis of DAE systems. Throughout this section we consider DAE systems

$$ S : f_i(\text{the } x_j\text{'s and their derivatives}) = 0 \qquad (13) $$

where $x_1, \dots, x_p$ are the variables and $f_1=0, \dots, f_n=0$ are the equations. The functions $f_i: \mathbb{R}^p \to \mathbb{R}$ are smooth, i.e., of class $C^m$ for sufficiently large integer $m$. System (13) is called *static* (also called *algebraic*) if no derivatives are involved.

* Call leading variables of System (13) the $d_j$-th derivatives$^{11}$ $x_j'^{d_j}$ for $j=1, \dots, p$, where $d_j$ is the maximal differentiation degree of variable $x_j$ through $f_1=0, \dots, f_n=0$.

* Remaining variables $x_j'^{m}$ for $j = 1, \dots, m$ and $0 \le m < d_j$ are called state variables.

* For $f=0$ an equation and $x$ a variable of $S$, we denote by $\sigma_{f,x}$ the maximal differentiation degree of $x$ in $f$; by convention, $\sigma_{f,x} = -\infty$ if $x$ does not occur in $f$.

<footnote n="11">¹¹The notation $x'^{k}$ is adopted throughout this paper, instead of $x^{(k)}$, for the k-th derivative of $x$.</footnote>
System (13) is denoted by

$$S = (F, X), \quad \text{where} \quad F = \{f_i \mid i = 1, \dots, n\}, \quad \text{and} \quad X = \{x_j \mid j = 1, \dots, p\} \qquad (14)$$

In the following, we regard $X$ and $F$ as sets of variables and functions. In particular, the conjunction of systems $S_1$ and $S_2$ is $(F_1 \cup F_2, X_1 \cup X_2)$, which we write also $S_1 \cup S_2$.

**Convention 1** By abuse of language, we will identify the set of functions $F$ and the set of equations $F=0$ it defines. Hence, we will often call “equation” an $f \in F$. □

**Comment 1** We can alternatively regard $S$ in (13) as a system of static equations with the leading variables as dependent variables (unknowns). The state variables are free variables, however subject to certain consistency conditions, as we shall see in the Section 3.2 on the $\Sigma$-method. □

## 3.1 Structural analysis of systems of equations

In this section, we consider the subcase of *static systems*, i.e., systems of the form (13) involving no derivative. We will need to consider additional *free variables* collected in $Y$, whose value is set by some environment in the system $F(X, Y)=0$. By contrast, we call *dependent variables* the variables collected in $X$, meaning that $F(X, Y)=0$ is to be solved for $X$.

The generic term of *structural analysis* refers to any analysis of a system of equations that relies on its incidence graph only—such analyses are therefore much cheaper than numerical ones. To $S = (F, X)$ (since $Y$ is irrelevant here), we associate its *incidence graph*, which is a nondirected bipartite graph $G_S = (F \cup X, E)$, where $(f, x) \in F \times X$ is an edge of $G_S$ if and only if variable $x$ occurs in equation $f$.

A matching is a subset $M \subseteq E$ involving at most once each equation and variable of $S$. Matching $M$ is perfect (we also say *complete*) if all equations and all variables of $S$ are involved (implying $p=n$). Matching $M$ is *variable-complete* (respectively *equation-complete*) if all variables (respectively equations) are involved in it. A vertex of $G_S$ is called *unmatched* in $M$ if it does not occur in $M$.

**Definition 1** Say that system $S$ is structurally nonsingular (or structurally regular) if its incidence graph $G_S$ possesses a perfect matching $M$.

Structural nonsingularity is a necessary (and generically sufficient) condition for the existence and uniqueness of solutions, see [5]. This condition is widely used in all software dealing with very large, but sparse, systems of equations. Structural analysis only deals with the incidence graph of a system of equations. Hence, in structural analysis, the following rule is enforced:

if variable $x$ occurs with differentiation degree $m$ in equation $f=0$, then it occurs with differentiation degree $m+1$ in equation $f'=0$. (15)

For $S$ a (possibly nonsquare, i.e., $p \neq n$) system of equations, the *Dulmage-Mendelsohn (DM) decomposition* [34, 3, 4] of $G_S$ uniquely partitions $S$ into three pairwise disjoint subsystems:

$$S = S^o \cup S^r \cup S^u, \qquad (16)$$

where $S^r$ is *structurally regular*, whereas $S^o$ and $S^u$ are the structurally *over-determined* and *under-determined* parts of $S$. System $S$ is structurally regular if and only if its Dulmage-Mendelsohn decomposition yields $S^o = S^u = \emptyset$.

**Notations 2** If $M \subseteq F \times X$ is a matching, the notation $f \in M$ means that there exists $x \in X$ such that $(f, x) \in M$; the corresponding meaning holds for the notation $x \in M$. □
## 3.2 The Σ-method for DAEs

Here we consider square DAE systems, meaning that $p=n$ in (13). We recall the so-called *index reduction*, by which suitably differentiating the different equations of a DAE system makes it ODE-like. Suppose we have a solution to the following problem:

Find a perfect matching $\mathcal{M}$ for $\mathcal{G}_S$ and integer valued equation offsets $\{c_f | f \in F\}$ and variable offsets $\{d_x | x \in X\}$, satisfying the following conditions:

$$ \begin{array}{l} d_x - c_f \geq \sigma_{f,x} \quad \text{with equality if } (f, x) \in \mathcal{M} \\ c_f \geq 0 \end{array} $$

Following Comment 1, we denote by $S_\downarrow$ the static system collecting the leading equations, i.e., such that:

$$ S_{\downarrow} = (F_{\downarrow}, X_{\downarrow}) \quad \text{where} \quad \begin{cases} F_{\downarrow} = \{f^{\prime c_f} | f \in F\} \\ X_{\downarrow} = \{x^{\prime d_x} | x \in X\} \end{cases} \qquad (18) $$

and we call it the *index reduced system*. System $S_\downarrow$ has the following properties: inequality $d_x \geq c_f + \sigma_{f,x}$ holds for each $x \in X$, and $d_x = c_f + \sigma_{f,x}$ holds for the unique $f$ such that $(f, x) \in \mathcal{M}$. Consequently, the set of pairs

$$ \mathcal{M}_{\downarrow} =_{\text{def}} \left\{ (f^{\prime c_f}, x^{\prime d_x}) \mid (f, x) \in \mathcal{M} \right\} \qquad (19) $$

is a perfect matching for the system $S_\downarrow$, seen as a system of static equations having $x^{\prime d_x}$ as dependent variables, hence this system is structurally nonsingular, see Definition 1. Hence, having reduced the index of $S$ makes it ODE-like: $S_\downarrow$ (generically) uniquely determines the values of leading variables, knowing the values of state variables (belonging to $X_\uparrow$).<fnref n="12" /> As announced in footnote 1, the consideration of latent equations governs how many times each $f \in F$ should be differentiable.

The static system

$$ S_{\uparrow} = (F_{\uparrow}, X_{\uparrow}) \quad \text{where} \quad \begin{cases} F_{\uparrow} = \{f^{\prime l} | f \in F, 0 \le l < c_f\} \\ X_{\uparrow} = \{x^{\prime k} | x \in X, 0 \le k < d_x\} \end{cases} \qquad (20) $$

defines the *consistency conditions*. System $S_\uparrow$ sets the constraints that must be satisfied by any consistent initial condition for the considered DAE [31]. Assuming that the Σ-method has a solution, $S_\uparrow$ possesses no overdetermined subsystem; hence, consistent valuations for the state variables exist, and each of them determines a value for the leading variables. Actually,

$$ \mathcal{M}_{\uparrow} =_{\text{def}} \left\{ (f^{\prime m}, x^{\prime d}) \mid \begin{array}{l} f \in F \\ 0 \le m < c_f \\ d = d_x - c_f + m \end{array} \right\} \qquad (21) $$

defines an equation-complete matching for $S_\uparrow$.

Finally, the static system

$$ S_{\uparrow} = (F_{\uparrow}, X_{\uparrow}) =_{\text{def}} (F_{\downarrow} \cup F_{\uparrow}, X_{\downarrow} \cup X_{\uparrow}) \qquad (22) $$

is called the *completion* of DAE $S$. Referring to the cup-and-ball example, (6) shows the decomposition $S_\downarrow = S_\downarrow \cup S_\uparrow$.

**Comment 2** Systems $S_\downarrow$, $S_\uparrow$, and $S_\uparrow$ are all static, meaning that $x, x', \dots, x'^k, \dots$ are seen as different variables, not as successive derivatives of the same variable $x$—they are called *dummy derivatives* in the DAE literature [28]. □

<footnote n="12">¹²The very elegant Σ-method, proposed by J. Pryce in 2001 [35] translates Problem (17) into a pair of primal/dual Linear Programs, of which the primal determines $\mathcal{M}$ and the dual determines the offsets. This proves that variable and equation offsets are independent of the particular choice for $\mathcal{M}$. Furthermore, it is proven that offsets can be selected pointwise minimal, and that such a minimal solution is unique. Index reduction is performed by all major DAE based modeling tools, by using the Σ-method or the original method by C. Pantelides [31]. See [3, 4] for omitted details.</footnote>
### 3.3 Implicit function theorem

This section reports useful results on the structural analysis of static systems of equations involving a small parameter. More precisely, we consider a system

$$F(X, Z, \varepsilon) = 0, \qquad (23)$$

where $F$ is a finite set of $C^1$-functions, $X$ and $Z$ are the sets of dependent and free variables, and $\varepsilon>0$ is a small parameter. Recall the following result, which is a rephrazing of the Implicit Function Theorem (see, e.g., Theorem 10.2.2 in [15]):

**Proposition 1 (Implicit Function Theorem)** Assume that the valuation $(x, z)$ for the pair $(X, Z)$ satisfies $F(x, z, 0) = 0$ and the Jacobian matrix $\partial F/\partial X$ at $(x, z, 0)$ is nonsingular. Then:

1. There exists a neighborhood of $(z, 0)$ in which, for every pair $(\hat{z}, \varepsilon)$, there exists $\hat{x}$ such that $F(\hat{x}, \hat{z}, \varepsilon) = 0$.

2. For any sequence $(z_n, \varepsilon_n)$ converging to $(z, 0)$, the solution $x_n$ of system $F(X, z_n, \varepsilon_n)=0$ converges to $x$.

## 4 Problem setting: hot restart of mode changes

In this section we formalize the hot restart problem, that is, how to properly generate deterministic restart conditions for a new mode, knowing the previous mode. We consider mode changes in isolation. Cascades of transient modes (of zero duration, as in the cup-and-ball with elastic impact) and Zeno phenomena (in which events of mode change accumulate in a finite duration of time) are not addressed.

**Definition 2 (mode change)** A mode change is defined as a triple $(S^-, t_*, S)$, where $S^- = (F^-, X^-)$ and $S = (F, X)$ are square<fnref n="13" /> DAE systems of the form (13) and $t_* \in \mathbb{R}_+$. $S^-$ is the previous mode, $t_*$ is the instant of mode change, and $S$ is the new mode.

Without loss of generality, we can assume that:

1. $F^-$ is index-reduced, i.e., $F^- = F_{\downarrow}^-$; and

2. $F$ is completed with its latent equations, i.e., $F = F_{\uparrow}$.

Recall that $X_{\downarrow}^-$ (resp., $X_{\uparrow}$) denotes the set collecting the leading and state variables of the previous (resp., new) mode. Then, we consider the subset $X^+ \subseteq X_{\uparrow}$ of state variables of the new mode. The hot restart consists in finding values for all of them.

As announced before, mode changes can be time- or state-based. More precisely, the following assumption holds throughout this work regarding mode changes:

**Assumption 1 (zero-crossing)** Mode change $(S^-, t_*, S)$ is caused by a zero-crossing, i.e., the crossing of zero from below by a $C^1$-function $g(X_{\downarrow}^-)$. Formally, setting $g(t) =_{\text{def}} g(X_{\downarrow}^-(t))$, there exists some duration $\delta>0$ such that, if the dynamics of $S^-$ is enforced everywhere, then, $g(t) < 0$ holds for $t \in (t_*-\delta, t_*)$ and $g(t) > 0$ holds for $t \in (t_*, t_*+\delta)$.

**Problem 1 (hot restart)** For $(S^-, t_*, S)$ a mode change, construct its hot restart system, which is a system of equations $R(X_{\downarrow}^-, X^+) = 0$, relating the left limits of variables and their derivatives at the mode change to the state variables for the new mode, satisfying the following conditions:

1. System $R(X_{\downarrow}^-, X^+)$ should be deterministic, meaning that the values for $X^+$ should be uniquely determined from the values for $X_{\uparrow}$;

<footnote n="13">¹³With as many equations as dependent variables.</footnote>
2. The restart should be consistent, meaning that the values for $X^+$ should satisfy the consistency conditions $F_\uparrow$;

3. After proper rescaling to compensate for possible impulsive behaviors, every invariant dynamics (i.e., common to previous and new mode) should be satisfied at the mode change.

Condition 1 reflects Requirement 1.1; condition 2 reflects Requirement 1.2. Condition 3 is partly informal, as “impulsive” and “rescaling” are not formally defined yet. This condition will be formalized and studied in Section 7.3. It is vacuously true if no invariant dynamics exists; otherwise, it indicates that hot restart is a strict refinement of cold start.

Our solution for hot restart will actually take the form of a system of equations

$$ \bar{R}(X_{\downarrow}^{-}, Z, X^{+}) = 0 $$

where $Z$ collects auxiliary variables. Eliminating $Z$ from $\bar{R}$ will yield a system

$$ R(X_{\downarrow}^{-}, X^{+}) = 0 $$

satisfying Requirements 1, 2, and 3. All of this was illustrated in the cup-and-ball example by (11,12), where $Z = \{\lambda^{\downarrow}\}$.

# 5 Mode change arrays

In this section, we introduce mode change arrays, as the main data structure of our approach. We first formalize the discrete time setting we introduced in Approach 1. For a fixed time step $\varepsilon > 0$, we use the time line

$$ \mathbb{T} =_{\text{def}} \{n\varepsilon \mid n \in \mathbb{N}\}, \qquad (24) $$

in a neighborhood of the mode change instant, where $\mathbb{N}$ denotes the set of nonnegative integers. Assuming $\varepsilon$ small enough and with reference to Assumption 1, we redefine the instant of detection of the mode change as being the first instant belonging to $\mathbb{T} \cap (t_{*}-\delta, t_{*}+\delta)$, such that the zero-crossing function $g(X_{\downarrow}^{-})$ is positive. To ensure that mode change detection remains causal, the new mode begins at instant $t_{*}+\varepsilon$; by abuse of notation, we redefine $t_{*}$ to be this instant. To summarize,

in our discretized setting, the zero-crossing is detected at instant $t_{*}-\varepsilon$ and the new mode is effective at instant $t_{*}$ (see Fig. 4).

<figure><img src="image_7.png" /><figcaption>Figure 4: Zero-crossing in discrete time. The detection instant is in red and the first instant of the new mode is in blue.</figcaption></figure>

## 5.1 Shifting and differentiating

Since we are interested in the system behavior around the instant of mode change, we take $t_{*}$ as time reference:
**Notations 3 (variables and functions)**

1. We assume an underlying set $\Xi$ of base variables. For $x \in \Xi$ and integers $k \in \mathbb{Z}, m \in \mathbb{N}$:

$$
\begin{align}
x &\text{ shall denote } x(t_*) \nonumber \\
x^{\bullet k} &\text{ shall denote } x(t_*+k\epsilon) \nonumber \\
x'^m &\text{ shall denote } x'^m(t_*) \nonumber \\
x'^{m\bullet k} &\text{ shall denote } x'^m(t_*+k\epsilon) \tag{26}
\end{align}
$$

and we denote by $\Xi'^{\bullet}$ the set of all $x'^{m\bullet k}$. Notations (26) extend to $\Xi'^{\bullet}$.

2. We also assume an underlying set $\mathbb{F}$ of base functions. For $f \in \mathbb{F}$ we define $f'^{m\bullet k}$ by $f'^{m\bullet k}(X) =_{\text{def}} f(X'^{m\bullet k})$ and $\mathbb{F}'^{\bullet}$ follows accordingly.

3. For $S$ a static system of equations over $\Xi'^{\bullet}$, we denote by $\mathcal{X}_S$ the set of variables occurring in $S$. When $S = (\{f\}, X)$, we simply write $\mathcal{X}_f$. $\square$

Following step 1 of Approach 1, derivatives are discretized by enforcing, for every $x \in \Xi'^{\bullet}$, the Euler identity

$$
x' = \epsilon^{-1}(x^{\bullet} - x). \tag{27}
$$

**Notations 4** For $y \in \Xi$ and $g \in \mathbb{F}$, we define:

$$
\begin{align*}
x &= y'^{m\bullet k} \in \Xi'^{\bullet} \quad \Rightarrow \quad (y_x, m_x, k_x) =_{\text{def}} (y, m, k) \\
f &= g'^{m\bullet k} \in \mathbb{F}'^{\bullet} \quad \Rightarrow \quad (g_f, m_f, k_f) =_{\text{def}} (g, m, k) \quad \Box
\end{align*}
$$

**Definition 3** Let $\sim \subseteq \Xi'^{\bullet} \times \Xi'^{\bullet}$ be defined by:

$$
x_1 \sim x_2 \quad \text{iff} \quad y_{x_1} = y_{x_2} \text{ and } [[x_1]] = [[x_2]],
$$

where $[[x]] = m_x+k_x$ is the total degree of $x$. Define the partial order $\leq$ on $\Xi'^{\bullet}$ by $x_1 \leq x_2$ iff $y_{x_1} = y_{x_2}$ and $[[x_1]] \leq [[x_2]]$, and write $\prec =_{\text{def}} \leq \wedge \neq$.

As mentioned in Section 3,

DAEs are seen as systems of static equations over base variables, their derivatives,
and their shifts, i.e., as systems of static equations over $\Xi'^{\bullet}$.
(28)

We now investigate some consequences of identity (27). Repeatedly applying this identity yields the
Euler identity

$$
y^m = \epsilon^{-n} \sum_{i=0}^{n} \binom{n}{i} (-1)^i y^{\bullet(n-i)}, \quad \text{where} \quad \binom{n}{i} = \frac{n!}{i!(n-i)!} \tag{29}
$$

is the binomial coefficient.

Consider two variables $x \sim z$. We have $x=y'^{m_x} \bullet k_x$ and $z=y'^{m_z} \bullet k_z$, where $m_x+k_x=m_z+k_z$. If, in addition, $x \neq z$ holds, then two cases can occur: either $m_x > m_z$, or $m_x < m_z$. We assume $m_x > m_z$, set $n = m_x - m_z = k_z - k_x > 0$, and apply (29) with this pair $(y, n)$. Differentiating $m_z$-times and shifting $k_x$-times the two sides of (29) yields the (generalized) Euler identity

$$
\begin{align}
0 &= x - \epsilon^{-n} \sum_{i=0}^{n} \binom{n}{i} (-1)^i z^{\bullet(-i)}, \text{ denoted by} \notag \\
0 &= \mathcal{E}_{xz}^{n}(x, z, U), \text{ where } U = (z^{\bullet(-i)})_{i=1, \dots, n} \tag{30}
\end{align}
$$

where tuple $U$ collects the variables that are $\prec x$.
**Convention 2** We write $\mathcal{E}_{xz}^n$ when the arguments need not be mentioned, or even $\mathcal{E}_{xz}$ when only $x$ and $z$ matter. □

We summarize the above discussion in the following lemma, where Notations 4 are used:

**Lemma 1** Any two different variables $x \sim z$ such that $m_x > m_z$ are related by Euler identity $\mathcal{E}_{xz}$.

**Definition 4** ($\sim$-closed system) Say that static system $S=(F, X)$, where $X \subseteq \Xi'^{\bullet}$, is $\sim$-closed if it includes all the Euler identities $\mathcal{E}_{xz}$ for any two $x, z \in \mathcal{X}_S$ such that $x \sim z$ and $m_x > m_z$. For any such system $S = (F, X)$, its $\sim$-closure, denoted by $\tilde{S}$, is defined by:

$$ \tilde{S} = \bigcap \left\{ T \mid \begin{array}{l} T \supseteq S \\ \forall x, z \in \mathcal{X}_T, \begin{bmatrix} x \sim z, x \neq z \\ m_x > m_z \end{bmatrix} \Rightarrow \mathcal{E}_{xz} \in T \end{array} \right\} \quad (31) $$

where $T$ ranges over the set of all extensions of system $S$.

The following result justifies the above definition:

**Lemma 2** $\tilde{S}$ always exists and is unique. It can be obtained by performing the following recursion, until fixpoint:

1. Initialization: $F_0 \leftarrow F, X_0 \leftarrow X$;

2. While $X_n \supset X_{n-1}$, update $(F_{n+1}, X_{n+1}) \leftarrow (F_n, X_n)$ by

* adding to $F_n$ every Euler identity $\mathcal{E}_{xz}$ such that $x, z \in X_n$, $x \sim z$, $x \neq z$, and $m_x > m_z$;

* adding to $X_n$ every $v$ involved in $\mathcal{E}_{xz}$.

*Proof* The recursion terminates in finitely many steps, since the total degree of the added variables is strictly decreasing. The limit of the recursion satisfies the fixpoint equation (31). □

**Convention 3** Unless the joint consideration of $S$ and its $\sim$-closure $\tilde{S}$ is needed, superscript “$\sim$” shall be omitted. □

## 5.2 Mode change arrays and facts

In this section, we start addressing Step 2 of Approach 1, for mode change $(S^{-}, t_{*}, S)$ following Definition 2. We use the notations from Section 3.2: for $F$ a DAE system, $F_{\downarrow}$, $F_{\uparrow}$, and $F_{\updownarrow}$ refer to the “index-reduced”, “consistency”, and “completion” of $F$ once index reduction has been performed. When considering mode changes, we use notations (25,26).

**Notations 5** Let $S = (F, X)$ be any DAE system, and let $t_{*}$ be a time origin. Then, for $k \in \mathbb{Z}$, $(F^{\bullet k}, X^{\bullet k})$ denotes the dynamics of the mode $S$ snapshot at instant $t_{*}+k\varepsilon$. □

For example, if function $f: x+y'$ belongs to $F$, then $f$ refers to $x(t_{*})+y'(t_{*})$, and $f^{\bullet k}$ refers to $x(t_{*}+k\varepsilon)+y'(t_{*}+k\varepsilon)$.

**Convention 4** Throughout this section, arrays are implicitly $\sim$-closed by completing them with the Euler identities (30) following Lemmas 1 and 2, and Convention 3 applies. □

To the new mode $S=(F, X)$ we associate the infinite array $A$ of equations obtained by stacking successive shifts of the completion $F_{\updownarrow}$:

$$ A =_{\text{def}} \bigcup_{k \ge 0} F_{\updownarrow}^{\bullet k} = \{ f^{\prime m \bullet k} \mid f \in F, 0 \le m \le c_f, k \ge 0 \}. \quad (32) $$
The set of past variables of A is its subset of variables whose value is determined by the previous mode:

$$ \mathcal{X}^{-} =_{\text{def}} \{ x \in \mathcal{X}_A \mid \exists x^{-} \in \bigcup_{k<0} (X^{-})^{\bullet k} \text{ s.t. } x^{-} \sim x \}. \quad (33) $$

The previous mode $S^{-}$ acts on the new mode $S$ only through the set $\mathcal{X}^{-}$ of past variables. To illustrate past variables, if $y'^{m} \in X^{-}$, then $x = y'^{m} \bullet (-k) \sim y^{\bullet (m-k)} \in \mathcal{X}^{-}$ whenever $0 < k \le m$.

Referring to (33), since $x^{-} \sim x$, $x$ is computed from $x^{-}$ via the Euler identity relating them (Lemma 1). Hence, the set of dependent variables of A is

$$ \mathcal{X}_{A}^{\text{dep}} =_{\text{def}} \mathcal{X}_{A} \setminus \mathcal{X}^{-}, \quad (34) $$

where $\mathcal{X}_A$ is the set of all the variables involved in A (both leading and state variables). The reason for taking the dependent variables of A as in (34) is that we follow the causality principle that the past cannot be undone. Variables that were set by the previous mode cannot be modified by the new mode, thus should not be considered as dependent variables of A. Besides this, all the variables involved in A need to be evaluated when performing the mode change (including the state variables).

At mode changes, additional equations hold, due to zero-crossings (Assumption 1):

**Lemma 3** If the mode change is caused by a zero-crossing of $g(X^{-})$, where the tuple $X^{-}$ of variables possesses continuous trajectories in the previous mode, then $g(X^{-})=0$ holds, up to an $O(\varepsilon)$, at instants $t_{*}, \dots, t_{*}+\ell\varepsilon$, where $\ell$ is the largest integer such that $(X^{-})^{\bullet\ell} \subseteq \mathcal{X}^{-}$.

*Proof* The continuity of trajectories of $X^{-}$ within the previous mode and the smoothness of function $g$, together imply that the value of $g(X_{\downarrow}^{-})$ changes by at most $O(\varepsilon)$ during a finite number of discrete time steps. $\square$

**Definition 5 (facts)** For $g(X_{\downarrow}^{-})$ a zero-crossing function:

1. Call root fact any equation that is entailed, up to $O(\varepsilon)$, by $g(X^{-}(t_{*}))=0$; the set $\Phi^{-}$ of all root facts is independent from the new mode;

2. With $g$ a root fact and $\ell$ as in Lemma 3, call fact any equation of A that is entailed, up to $O(\varepsilon)$, by the system

$$ g(X^{-})=0, \dots, g(X^{-})^{\bullet\ell}=0. $$

The set of all facts is denoted by $\Phi_A$.

Since facts hold true up to an $O(\varepsilon)$ and do not involve dependent variables, we set the following convention:

**Convention 5** In the sequel, facts are removed from A, i.e., A is replaced by $A \setminus \Phi_A$. $\square$

We continue addressing Step 2 of Approach 1 by introducing finite arrays, based on the notion of height. A height for mode $S = (F, X)$ denotes a function$^{14}$

$$ K : F \to \mathbb{N}, \text{ denoted } f \mapsto K_f. \quad (35) $$

In the following definition, A denotes the infinite array (32), and its set $\mathcal{X}_A^{\text{dep}}$ of dependent variables is defined by (34):

<footnote>$^{14}$Reasons for considering nonconstant heights are stated in Comment 4.</footnote>
**Definition 6 (mode change array)** Any height K for mode S=(F, X) defines the K-mode change array

$$A_K = \{ f'^{m} \bullet k \mid f \in F, 0 \le m \le c_f, 0 \le k \le K_f \} \subset A. \qquad (36)$$

The set of dependent variables of $A_K$ is $\mathcal{X}_{A_K}^{\text{dep}} = \mathcal{X}_A^{\text{dep}} \cap \mathcal{X}_{A_K}$. Define the tail of $A_K$ and its set of dependent variables by

$$A_K^{\text{tail}} = \left\{ f'^{m} \bullet K_f \mid \begin{array}{l} f \in F, f, x \in \mathcal{M} \\ 0 \le m \le c_f \\ d = \sigma_{f,x} + m \end{array} \right\}, \quad \mathcal{X}_{A_K}^{\text{tail}} = \left\{ x'^{d} \bullet K_f \mid \begin{array}{l} f \in F, (f, x) \in \mathcal{M} \\ 0 \le m \le c_f \\ d = \sigma_{f,x} + m \end{array} \right\} \\ \text{and} \\ \mathcal{X}_{A_K}^{\text{head}} = \text{def} \mathcal{X}_{A_K} \setminus \mathcal{X}_{A_K}^{\text{tail}} \qquad (37)$$

where $\mathcal{M}$ is the matching found by solving the $\Sigma$-method for S. We define the subset of equations that must be enabled by

$$A_K^{\text{ena}} = A_K^{\text{tail}} \cup (A_K \cap \bigcup_{k \ge 0} (F^{-})^{\bullet k}) \qquad (38)$$

and we set $A_K^{\text{dis}} = \text{def} A_K \setminus A_K^{\text{ena}}$.

The rationale for (38) is that 1) following Requirement 1, we want the tail to be enabled, and 2) the invariant equations (belonging to both modes) should never be disabled.

**Lemma 4** For any given base variable $x$ and $0 \le d \le d_x$, the set $\{z'^{m} \bullet k \in \mathcal{X}_{A_K}^{\text{tail}} | z=x, m=d\}$ possesses cardinality 1.

*Proof* Variable $x$ uniquely determines $f$ by the condition $(f, x) \in \mathcal{M}$. Hence, the height $K_f$ is uniquely defined. The lemma follows immediately. □

**Example 1** Fig. 2 shows the $\sim$-closed array $A_1$ (see Definition 4) associated to the mode change $\gamma: F \to T$ for the cup-and-ball. Disabled equations (in red) are included; $K=1$ is a constant; $\mathcal{X}^- = \{x, y, x', y', x^\bullet, y^\bullet\}$; and $\mathcal{X}_{A_1}$ collects all the variables involved in blocks occurring at $t$ and $t+\varepsilon$. $A_1^{\text{tail}}$ is the block of equations attached to instant $t+\varepsilon$ and $\mathcal{X}_{A_1}^{\text{tail}} = \{x'', y'', \lambda, x', y'\}^\bullet$. Note that $\mathcal{X}_{A_1}^{\text{tail}}$ has cardinality 5, whereas the cardinality of $A_1^{\text{tail}}$ is 4. Adding, to $A_1^{\text{tail}}$, Euler identity $(\mathcal{E}_{x''})$ makes it structurally nonsingular. □

## 5.3 Structural Implicit Function Theorem

We now provide the main mathematical tool justifying step 3 of Approach 1, when we set $\varepsilon$ to zero. This is achieved through a structural interpretation of Proposition 1 (Implicit Function Theorem).

We assume a reference instant $t_*$, a time step $\varepsilon$, and the resulting set $\Xi'^{\bullet}$ of variables, associated with them following Notations 3. On top of this, we consider a static system $S = (F, X \cup Z)$, $X \cup Z \subseteq \Xi'^{\bullet}$, where $X$ and $Z$ are the dependent and free variables of $S$. We assume that $S$ is structurally nonsingular and $\sim$-closed, see Definition 4.15 System $S$ involves dependent and free variables, as well as the parameter $\varepsilon$, due to the occurrence of Euler identities. $S$ is thus of the form $F(X, Z, \varepsilon)=0$, to which Proposition 1 applies. As its main smoothness condition, Proposition 1 requires that $\partial F / \partial X$ shall be invertible at a solution $(x, z, 0)$ of $F$.

As now usual in our approach, we weaken invertibility to generic invertibility (equivalently, structural nonsingularity). Accordingly, Proposition 1 weakens as follows:

**Proposition 2 (Structural Implicit Function Theorem)** If setting $\varepsilon$ to zero in the system $F(X, Z, \varepsilon) = 0$ yields a structurally nonsingular system, then:

<footnote>15An example is the array $A_1$ shown in Fig. 3.</footnote>
1. System $F(X, Z, 0) = 0$ structurally determines $X$ as a function of $Z$; and

2. For a same valuation of $Z$, the so obtained valuation of $X$ is the limit, when $\varepsilon \searrow 0$, of the valuation of $X$ defined by system $F(X, Z, \varepsilon) = 0$.

Statement 1 is equivalent to assuming that $F(X, Z, 0) = 0$ is structurally nonsingular, hence only Statement 2 has added value. Proposition 2 induces the following rule governing the setting of $\varepsilon$ to zero in the system $F(X, Z, \varepsilon)=0$:

**Rule 1** It is legitimate to set $\varepsilon$ to zero in system $F(X, Z, \varepsilon) = 0$, provided that the resulting system remains structurally nonsingular.

We will use Rule 1 as a validity criterion when performing Step 3 of Approach 1.

# 6 Rescaling analysis of mode change arrays

In this section we develop a rescaling analysis for mode change arrays. Conventions 3 and 4 are in force. This section focuses on exploiting array A introduced in (32), recalling that facts were removed from it (Convention 5).

We first specify more precisely the class of “model functions”, on top of which A is built.

**Definition 7 (model functions)** The set $\mathbb{F}'$ of model functions is defined as the set of all smooth functions defined on the set of variables $\Xi'$, augmented with the set of all Euler identities relating the variables belonging to $\Xi'$.

**Lemma 5** $\mathbb{F}'$ contains all the Euler identities, as well as the closure of $\mathbb{F}$ under time differentiation $f \mapsto f'$.

*Proof* The first statement is immediate. The second statement follows from the chain rule for differentiation. □

We now zoom on the features of model functions that we will be using in the rescaling analysis:

* For $f \in \mathbb{F}'$ a smooth function defined on $\Xi'$, we distinguish the subset $\mathcal{L}_f$ of its variables that enter linearly:

$$f = \sum_{x \in \mathcal{L}_f} x \times f_x, \qquad (39)$$

where each $f_x$ is a smooth nonlinear function of a finite subset of variables belonging to $\Xi' \setminus \mathcal{L}_f$.

* Euler identity (30), where $x \sim z$ and $n = \text{def } m_x - m_z > 0$, has the following form, where $a_i$ are constants:

$$\mathcal{E}_{xz} = x - \varepsilon^{-n} \sum_{i=0}^{n} a_i \times z^{\bullet(-i)}. \qquad (40)$$

* Classes (39) and (40) of functions are subsumed by the following larger class collecting all functions of the form

$$f = \sum_{x \in \mathcal{L}_f} \varepsilon^{-n_x} \times x \times f_x, \qquad (41)$$

where each $n_x$ is a nonnegative integer, and each $f_x$ is a smooth nonlinear function of a finite subset of variables belonging to $\Xi' \setminus \mathcal{L}_f$.
In the sequel, we use the form (41) for our model functions.

**Definition 8 (rescaling offset)** *Call rescaling offset a function $\mu : \mathbb{F}' \to \mathbb{N} \cup \{+\infty\}$, denoted $f \mapsto \mu_f$, satisfying the following axioms, where we write $\mu_x$ when $f \equiv x$, and $\mathcal{X}_f$ is defined in Notations 3:*

$$ \forall n \in \mathbb{N} \Rightarrow \mu_{\varepsilon^{-n}} = n \qquad (42) $$

$$ f \text{ of the form } (41) \Rightarrow \mu_f = \max_{x \in \mathcal{L}_f} (n_x + \mu_x + \mu_{f_x}) \qquad (43) $$

$$ f \text{ non-linear} \Rightarrow \mu_f = \text{if } \mu_{\mathcal{X}_f} = 0 \text{ then } 0 \text{ else } +\infty \qquad (44) $$

where $\mu_{\mathcal{X}_f} = 0$ means $\mu_x = 0$ for every $x \in \mathcal{X}_f$.

**Comment 3** *Axioms (42-44) are supported by the following intuition. The rescaling offset aims to represent the “impulse order”, namely $\mu_x = n$ captures that $x$ is $O(\varepsilon^{-n})$. Axiom (42) expresses that intent, and axiom (43) is immediate considering (41). Axiom (44) formalizes that, if $f$ is nonlinear, we do not know to assign an “impulse order” to it, unless it is known that all the variables involved in $f$ are not impulsive. □*

Say that a function is *impulsive* if its rescaling offset is positive. The *rescaling* of function $f$ is defined by

$$ f^{\downarrow} =_{\text{def}} \text{if } \mu_f < \infty \text{ then } \varepsilon^{\mu_f} f \text{ else "undefined"}. \qquad (45) $$

Of course, the intent of rescaling is that, if defined, $f^{\downarrow}$ will be non impulsive:

**Lemma 6** *For every $f$ such that $\mu_f < \infty$: we have $\mu_{f^{\downarrow}} = 0$ and $(f^{\downarrow})^{\downarrow} = f^{\downarrow}$.*

*Proof* Using (45), we get $f = \varepsilon^{-\mu_f} \times f^{\downarrow}$. Hence, by axiom (43), we get $\mu_f = \mu_f + \mu_{f^{\downarrow}}$, showing that $\mu_{f^{\downarrow}} = 0$. □

We are given a mode change $(S^{-}, t_{*}, S)$. To formalize Approach 1, we consider the following *rescaling problem*, where we use Notations 2. For $K$ a height, let $\mathcal{M}_K$ be a matching of $(A_K, \mathcal{X}_{A_K}^{\text{dep}})$, and let $(f, x_f) \in \mathcal{M}_K$, where $f$ has the form (41). Define

$$ f_{\downarrow} \mathcal{M}_K \qquad (46) $$

by restricting, in (41), the sum to the subset of its terms involving the variable $x_f$ (note that $x_f$ may or may not belong to $\mathcal{L}_f$). We denote this subset by

$$ \mathcal{L}_{f \downarrow \mathcal{M}_K} . $$

For example, referring to array $A_1$ of Fig. 2 and equation $(k_1^{\bullet})$ in it, this equation is matched with variable $y^{\bullet}$ in matching $\mathcal{M}_1$ highlighted in blue, hence $(k_1^{\bullet})_{\downarrow \mathcal{M}_1} = (yy')^{\bullet}$.

**Problem 2 (rescaling problem)** *For a given pair $(A, \mathcal{X}^{-})$ following (32, 33), find:*

1. a height $K$ following (35);

2. a variable-complete matching $\mathcal{M}_K$ for $A_K$, satisfying (see (38) for the definition of $A_K^{\text{ena}}$)

$$ \forall f \in A_K^{\text{ena}} \Rightarrow f \in \mathcal{M}_K \qquad (47) $$

and such that the following system of equations over rescaling offsets has a solution $\mu = (\mu_x)_{x \in X}$ that is finite for every $x$:

$$ \forall f \in \mathcal{M}_K \Rightarrow \mu_f = \mu_{(f_{\downarrow} \mathcal{M}_K)}, \quad \text{where } f_{\downarrow} \mathcal{M} \text{ was defined in (46).} \qquad (48) $$
Triple $(K, \mathcal{M}_K, \mu)$ is a solution of the rescaling problem.

Expanding equation (48) by using Axiom (43) yields

$$ \max_{x \in \mathcal{L}_f} (n_x + \mu_x + \mu_{f_x}) = \max_{x \in \mathcal{L}_{f \downarrow \mathcal{M}_K}} (n_x + \mu_x + \mu_{f_x}) \quad (49) $$

which, by applying Axiom (44) to $\mu_{f_x}$, reveals that (48) is indeed to be solved for $(\mu_x)_{x \in X}$.

**Example 2** These notions are illustrated on the cup-and-ball example (Section 2). In Fig. 1, $A_0^{\text{ena}}$ collects equations $(e_1, e_2)$ since they are common to both modes, and $(k'_1)$ since it is a consistency equation ot the last instant of the array. Thus, $A_0^{\text{dis}} = \emptyset$; hence, disabling $(k'_1)$ violates this requirement. In Fig. 2, $A_1^{\text{tail}}$ collects all the equations attached to the last instant $t+\varepsilon$; $A_1^{\text{head}} = A_1 \setminus A_1^{\text{tail}}$; $A_1^{\text{dis}}$ collects equations $(k'_1, k''_1)$, since they belong to the new mode only and do not sit in the last instant; $A_1^{\text{ena}}$ is the complement. The solution of the corresponding rescaling problem consists of the matching highlighted in blue, and the rescaling offsets listed in (8). $\square$

Using a solution $(K, \mathcal{M}_K, \mu)$, we can apply the following procedure to (hopefully) generate the system of equations characterizing the hot restart. In this procedure, we make time step $\varepsilon$ explicit, since we let $\varepsilon := 0$ in its step 3.

**Procedure 1** Let $(K, \mathcal{M}_K, \mu)$ be a solution of Problem 2, and let $A_\varepsilon$ be the subarray of $A_K$ collecting the equations matched in $\mathcal{M}_K$. Perform the following:

1. erase from $A_\varepsilon$ all equations not matched in $\mathcal{M}_K$;

2. rescale the remaining equations and variables using (45); the resulting rescaled array is denoted by $A_\varepsilon^\perp$;

3. set $\varepsilon := 0$ in $A_\varepsilon^\perp$, and denote the result by $A_0^\perp$;

4. rename as follows all the variables occurring in $A_0^\perp$:

$$ \begin{align} x \in \mathcal{X}^- &\leftarrow x^- \nonumber \\ x'^{m \cdot k} \in \mathcal{X}_{A_K}^{\text{tail}} &\leftarrow (x'^m)^+ \quad \text{(restart variables)} \tag{50} \\ x'^{m \cdot k} \in \mathcal{X}_{A_K}^{\text{head}} &\leftarrow x'^{m \cdot k} \quad \text{(auxiliary variables)} \nonumber \\ x'^{m \cdot k} \notin \mathcal{X}_{A_K} &\leftarrow \text{undefined} \nonumber \end{align} $$

Let $A_\star$ be the array obtained from $A_\varepsilon$ by performing the above sequence of operations: we call it the restart system.

For the cup-and-ball, this procedure was illustrated in (9,10,11).

The restart system has a unique solution under certain conditions, as will be seen in Theorem 1 and its proof. These conditions define good solutions of the rescaling problem:

**Definition 9 (good solution)** Triple $(K, \mathcal{M}_K, \mu)$ is a good solution of the rescaling problem if the following holds:

$$ \forall f \in \mathcal{M}_K \quad \Rightarrow \quad \mu_f < \infty \quad (51) $$

$$ \forall x \in \mathcal{X}_{A_K}^{\text{tail}} \quad \Rightarrow \quad \mu_x = 0 \quad (52) $$

$$ \forall x = y'^{m \cdot k} \in \mathcal{X}_{A_K}^{\text{head}} \quad \Rightarrow \quad y'^{(m-n) \cdot (k+n)} \in \mathcal{X}_{A_K} \quad (53) $$

$$ \text{where } n =_{\text{def}} \min(\mu_x, m) \quad \square $$

Condition (51) expresses that the enabled array must be rescalable. Condition (52) expresses that restart variables should not be impulsive. Condition (53) will ensure that the last case of renaming (50) never occurs.
**Comment 4** Note that, if the new mode S is the composition of two noninteracting DAE systems, then Problem 2 also decomposes into two independent rescaling problems. Requiring heights to be constant would then be inadequate. □

The solution (8) for the rescaling problem of the cup-and-ball was good. These results are further illustrated by a series of examples in Appendices D and E. Differences with the approach by Trenn et al. are illustrated in Appendix F.

# 7 Main theorems

## 7.1 Correctness of our approach

**Theorem 1** The hot restart system obtained by solving Problem 2 and then applying Procedure 1, satisfies the following properties:

1. Given $(K, \mathcal{M}_K)$ satisfying (47), system (49) possesses a solution $\mu$, generically.

2. If solution $(K, \mathcal{M}_K, \mu)$ is good, then:

(a) rescaled array $A_{\varepsilon}^{\downarrow}$ involves only variables belonging to $\mathcal{X}_{A_K}$;

(b) $A_0^{\downarrow}$ is structurally nonsingular;

(c) $A_0^{\downarrow}$ is the unique limit of $A_{\varepsilon}^{\downarrow}$ when $\varepsilon \searrow 0$; and

(d) restart system $A_*$ uniquely determines the states $(y^m)^+$ for the hot restart of the new mode.

3. If good solutions $(K, \mathcal{M}_K, \mu)$ exist, then height $K$ and rescaling offsets $\mu$ are unique; the only possible source of non uniqueness is indeed the matching $\mathcal{M}_K$.

In Statement 1, “generically” refers to the structural analysis setup. More precisely, we extend $\mu$ to take values in $\mathbb{R}_{\ge 0}$ and exponents of $\varepsilon^{-1}$ to belong to $\mathbb{R}_{\ge 0}$. Problem 2 extends to this situation. Then, “generically” means: for almost every value for the exponents of $\varepsilon^{-1}$ in some neighborhood of their nominal values.

*Proof* See Section 8. □

Theorem 1 allows us to reject a model at compile time if no good solution is found, but gives little feedback in this case, regarding how to correct the model. The following corollary of Theorem 1 improves the feedback to the designer by identifying the subset of variables that are insufficiently determined. Notations of Theorem 1 are used.

Consider the array $A_{\varepsilon}$ and assume that goodness conditions (51) and (52) are satisfied. Let $\mathcal{Y}$ be the subset of variables $x$ violating the goodness condition (53). Remove, from $A_{\varepsilon}$, the equations that are matched with a variable belonging to $\mathcal{Y}$, and keep $\mathcal{X}_{A_{\varepsilon}}^{\text{dep}}$ unchanged. The resulting system $(A_{\varepsilon}^{\mathcal{Y}}, \mathcal{X}_{A_{\varepsilon}}^{\text{dep}})$ is now rectangular, with more dependent variables than equations. Its Dulmage-Mendelsohn decomposition (16) yields $A_{\varepsilon}^{\mathcal{Y}} = A_{\varepsilon}^{r} \oplus A_{\varepsilon}^{u}$ (the overdetermined subsystem is empty). Let $\mathcal{M}_K^r$ be a perfect matching for $(A_{\varepsilon}^{r}, \mathcal{X}^r)$.

**Corollary 1** Consider the objects of Theorem 1 and perform the substitutions

$$A_{\varepsilon} \leftarrow A_{\varepsilon}^{r}, \mathcal{X}_{A}^{\text{dep}} \leftarrow \mathcal{X}^{r}, \mathcal{M}_{K} \leftarrow \mathcal{M}_{K}^{r} \qquad (54)$$

whereas the rescaling offsets $\mu_x$ and $\mu_f$ are kept unchanged for $x \in \mathcal{X}^r$ and $f \in A_{\varepsilon}^r$. Let $A_*^r$ be the array obtained by performing, mutatis mutandis, the operations listed in Procedure 1. Then, $A_*^r$ determines the values of state variables belonging to $\mathcal{X}^r$ for a hot restart of the new mode. Other state variables may be undetermined.
*Proof* See Appendix A.

<figure><img src="image_8.png" /></figure>

In Appendix D.2 we analyze a slight modification of the cup-and-ball game, by which the mode change is no longer state based, but rather an external input. For this example, the rescaling problem possesses no good solution. However, applying Corollary 1 shows that restart positions are determined, but restart velocities are not.

So far nothing tells us how large the height $K$ for the array should get explored. The next result proposes bounds for it:

**Theorem 2** The exist two heights $K_* \le K^*$, depending on the equation offsets of the new mode $S$ and the set $\mathcal{X}^-$ of past variables, such that:

1. To find a variable-complete matching $\mathcal{M}_K$ satisfying requirement (47), it is enough to select $K \ge K^*$, whereas no such matching can be expected unless $K \ge K_*$;

2. Let $K$ be the variable height of a successful solution of Problem 2. Then, $K \le K^*$ actually holds.

*Proof* See Appendix C.

<figure><img src="image_9.png" /></figure>

For the cup-and-ball example, we have $K_* = K^* = 1$. Although not true in general, the equality $K_* = K^*$ is expected to hold in practice. Theorem 2 is illustrated in Appendix D.1.

## 7.2 Solving the rescaling problem

Recall our functions $f$ belong to class (41). However, since our actual equations will be of the form either (39) or (40), the following assumptions can be formulated on class (41):

**Assumption 2** The functions $f$ of class (41) that we consider, satisfy the following conditions:

$$ \exists x \in \mathcal{L}_f \text{ s.t. } n_x > 0 \quad \text{implies} \quad \mathcal{L}_f = X \qquad (55) $$

$$ x_f \in \mathcal{L}_f \quad \text{implies} \quad \mathcal{L}_{f \downarrow \mathcal{M}_K} = \{x_f\} \qquad (56) $$

**Justification** The left hand side of (55) holds only if $f$ is of class (40), i.e., a Euler identity; in this case, all the variables are linearly involved, whence the right hand side. Then, (56) holds for every function $f$ of class (41).

<figure><img src="image_10.png" /></figure>

**Theorem 3** Assumption 2 is in force.

1. Finding a good solution for (49) decomposes as the following two steps, where

$$ \mathcal{M}_{\mathcal{L}} = \text{def} \{ (f, x_f) \in \mathcal{M}_K \mid x_f \in \mathcal{L}_f \} $$

(a) Find a finite minimal nonnegative solution for $\mu_{\mathcal{L}} = \text{def} \{ \mu_{x_f} \mid (f, x_f) \in \mathcal{M}_{\mathcal{L}} \}$, of the system

$$ \forall f \in \mathcal{M}_{\mathcal{L}} : \max_{x \in \mathcal{L}_{f \downarrow \mathcal{M}}} (n_x + \mu_x) = \max_{x \in \mathcal{L}_f} (n_x + \mu_x) \qquad (57) $$

(b) Extend this solution to a good solution for (49) by checking the satisfiability of the following property:

$$ \forall f \in \mathcal{M}_K, \forall x \notin \mathcal{L}_f \Rightarrow \mu_x = 0 \qquad (58) $$

2. Rescaling equation (57) simplifies as

$$ \forall f \in \mathcal{M}_{\mathcal{L}} : n_{x_f} + \mu_{x_f} = \max_{x \in \mathcal{L}_f} (n_x + \mu_x) \qquad (59) $$
Proof We first prove Statement 1. Consider, as part of the rescaling problem, the system of equations (49). Suppose it possesses a good solution. By (51) we look for a solution ensuring that all the variables and enabled equations possess a finite rescaling offset. If such a solution $\mu$ exists, then:

1. For $(f, x_f) \in \mathcal{M}_L$, (57) holds, and

2. Property (58) holds.

In addition, if a solution is good, then any smaller solution is also good. By Statement 3 of Theorem 1, there is at most one good solution. Hence the good solution, if any, must be a minimal solution of (57). This proves Statement 1.

Focus on Statement 2. By Property (56) of Assumption 2, the max sitting on the left hand side of (57) boils down to the singleton $x = x_f$. Thus (57) rewrites as (59). □

System (59) is of max/+ type. Yet, it belongs to a subclass amenable of much more efficient algorithms. More precisely, straightforward manipulations bring (59) to a set of inequalities of the form $\mu_x - \mu_y \le D_{xy}$, where $x, y \in X \cup \{0\}$, the unknowns are $\mu_x$ taking values in $\mathbb{Z} =_{\text{def}} \mathbb{Z} \cup \{\pm\infty\}$, $0$ is a distinguished element such that $\mu_0 = 0$, and $D_{xy}$ is a matrix of entries belonging to $\mathbb{Z}$. Such systems of inequality constraints are known as Difference Bound Matrices (DBM) [16, 29]. DBM can also be used for the dual problem of the $\Sigma$-method, for which they provide a modular (compositional) algorithm, see [10].

## 7.3 Characterizing the solutions

In this section, we complete the proof that our approach solves the hot restart problem, by formalizing point 3 of Problem 1 and proving that our solution solves it.

**Background regarding rings and ideals:** Our reference is the textbook [14]. The set of solutions of a (static) system of smooth equations is best studied by using the notions of ring $\mathbb{F}'$ of real-valued functions $f$, and of ideal of this ring. For $S = (F, X)$ such a system, we denote by $V_S$ its set of solutions, i.e., the set of all valuations $\nu_X \in \text{dom}(X)$ such that $F(\nu_X) = 0$. If $S = S_1 \cup S_2$, then $V_S = V_{S_1} \cap V_{S_2}$, where the intersection symbol first requires to take the inverse projections of $V_{S_1}$ and $V_{S_2}$, from $X_1$ and $X_2$ to $X_1 \cup X_2$. The set of all functions $f$ such that $F(\nu_X) = 0$ implies $f(\nu_X) = 0$ is an ideal, i.e., is closed under addition and multiplication by an arbitrary element of $\mathbb{F}'$, we denote this ideal by $\mathcal{I}_S$. If $S = S_1 \cup S_2$, then $\mathcal{I}_S = \mathcal{I}_{S_1} + \mathcal{I}_{S_2} =_{\text{def}} \{f_1 + f_2 \mid f_i \in \mathcal{I}_{S_i}\}$.

### Analyzing the cup-and-ball example

Consider the cup-and-ball example. The subsystem common to the two modes is

$$S^{\pm} : \begin{cases} 0 = x'' + \lambda x \\ 0 = y'' + \lambda y + g \end{cases}$$

The algebraic variable (the tension $\lambda$) plays no role in hot restart. Eliminating it from $S^{\pm}$ yields the equation

$$f_{\pm} = 0, \text{ where } f_{\pm} = y''x - x''y + gx, \qquad (60)$$

which belongs to the ideal generated by $S^{\pm}$. However, equation (60) exhibits impulses at the mode change, and we have $\mu_{f_{\pm}} = \mu_{x''} = \mu_{y''} = 1$. Therefore:

$$x'' \downarrow = \varepsilon x'' = x' - x'$$

$$y'' \downarrow = \varepsilon y'' = y' - y'$$

$$f_{\pm} = \varepsilon f_{\pm} = (y' - y)x - (x' - x)y + 0(\varepsilon)$$
which, by letting $\epsilon := 0$, yields

$$0 = (y'_{\bullet} - y')x - (x'_{\bullet} - x')y. \qquad (61)$$

Focus on our solution (11) for the restart of the cup-and-ball, which yields (12) by eliminating $\lambda^{\downarrow}$ from it. Observe that (61) identifies with (12) after renaming (50) is performed.

### Generalizing

We now generalize the analysis of the cup-and-ball example. Following Section 6, we consider the ring $\mathbb{F}'_{\bullet}$ of functions of the variables belonging to $\Xi'_{\bullet}$ augmented with $\{\epsilon^{-n} | n \ge 0\}$. We also consider the ring $\mathbb{M} =_{\text{def}} (\mathbb{N} \cup \{+\infty\}, \max, +)$ of rescaling offsets equipped with the $\max/+-$ algebra.

From now on, we consider $(S^{-}, t_{*}, S)$ a mode change, and we assume that Problem 2 possesses a solution $(K, \mathcal{M}_{K}, \mu)$ satisfying goodness conditions (51–53). By axioms (42–44), the following result easily follows regarding $\mu$:

**Lemma 7** $\mu : \mathbb{F}'_{\bullet} \to \mathbb{M}$ is a ring homomorphism.

Consider the map $\Lambda : \mathbb{F}'_{\bullet} \to \mathbb{F}'_{\bullet}$ consisting in setting $\epsilon:=0$ in $f^{\downarrow}$:

$$\Lambda(f) =_{\text{def}} f^{\downarrow}[\epsilon:=0]. \qquad (62)$$

**Lemma 8** Writing $\mu_{i}$ instead of $\mu_{f_{i}}$, map $\Lambda$ satisfies

$$\Lambda(f_1+f_2) = \sum_{i:\mu_i=\max(\mu_1,\mu_2)} \Lambda(f_i)$$

$$\Lambda(f_1 \times f_2) = \Lambda(f_1) \times \Lambda(f_2)$$

*Proof* See Appendix B. $\square$

The following theorem shows that our solution meets Requirement 3 of Problem 1 (by the way, this requirement is now formalized):

**Theorem 4** Every function belonging to the ideal spanned by the subsystem that is common to both modes, i.e., $S^{\pm} = S^{-} \cap S$, is mapped via $\Lambda$ to a function that belongs to the ideal generated by our solution for restart.

*Proof* Convention 4 is in force throughout this proof. Let $A_{\epsilon}$ be the mode change array associated to the considered solution $(K, \mathcal{M}_{K}, \mu)$, and let $A_{\epsilon}^{\mathcal{M}}$ be the subarray collecting the enabled equations (i.e., matched in $\mathcal{M}_{K}$). Rephrasing Theorem 2 by using map $\Lambda$, and seeing arrays as sets of equations, we have:

$$\begin{aligned} A_{0}^{\downarrow} &= \Lambda(A_{\epsilon}^{\mathcal{M}}) \text{ is structurally nonsingular, and} \\ A_{\star} &= A_{0}^{\downarrow} \text{ followed by renaming (50).} \end{aligned} \qquad (63)$$

Define the ideals

$$\mathcal{I}_{\epsilon} =_{\text{def}} \text{span}(A_{\epsilon}^{\mathcal{M}}) \quad \text{and} \quad \mathcal{I}_{0} =_{\text{def}} \text{span}(A_{0}^{\downarrow}). \qquad (64)$$

Using Lemma 8 yields

$$\forall f \in \mathcal{I}_{\epsilon} \Rightarrow \Lambda(f) \in \mathcal{I}_{0}. \qquad (65)$$

Let $A_{\epsilon}^{\pm}$ be the subarray of $A_{K}^{\mathcal{M}}$ collecting shifted versions of the equations involved in both modes. By (47), invariant equations cannot be disabled, i.e., $A_{\epsilon}^{\pm} \subseteq A_{\epsilon}^{\text{ena}}$. Hence equations belonging to $A_{\epsilon}^{\pm}$ are all matched in $\mathcal{M}_{K}$. Therefore, $A_{\epsilon}^{\pm} \subseteq \mathcal{I}_{\epsilon}$. Hence, by (65), $\Lambda(A_{\epsilon}^{\pm}) \subseteq \mathcal{I}_{0}$, which proves the theorem. $\square$

Array $A_{0}^{\downarrow}$ involves the following categories of variables:
1. free variables belonging to $\mathcal{X}^{-}$; they are non-impulsive;

2. dependent state variables sitting in the tail, i.e., of the form $x^{m} \bullet K_{f}$ for $(f, x) \in \mathcal{M}$, the solution of Problem (17) and $0 \le m < c_{f}$; they are also non-impulsive;

3. other variables, possibly rescaled and/or auxiliary.

Variables of categories 1 and 2 are physically interpretable (left- and right-limits of state variables at the mode change). In contrast, variables of category 3 are not: it makes sense to eliminate them. The resulting system of equations is included in ideal $\mathcal{I}_{0}$. Such an elimination can be performed using numerical methods for linear systems, and Gröbner bases and elimination theory [14] for polynomial systems. For the latter class, large systems are totally beyond reach. Even for linear systems, elimination is of order $n^{3}$, where $n$ is the dimension of the mode change array.

# 8 Proof of correctness

In this section we present the proof of the main result, namely the correctness of our approach formalized by Theorem 1. Other proofs are deferred to appendices. Throughout this section, Conventions 3, 4, and 5 are in force. We will use the following result regarding rescaling offsets:

**Lemma 9** Let $x, z$ be as in Euler identity (30), and let $\mu$ be a rescaling offset. Then, $\mu_{x} = n + \max_{0 \le i \le n} \mu_{(z \bullet (-i))}$.

*Proof* The lemma is a direct consequence of identity (30) and Axiom (43). □

We now proceed with the proof of Theorem 1.

**Proof of Statement 1 of Theorem 1:** Let $\chi$ be the Heaviside function and let $\chi^{\eta}$ be a smoothing of it, such that $0 \le \chi^{\eta} \le 1$ and $\chi^{\eta} = \chi$ outside the interval $[0, \eta]$, with $\eta$ being a small positive parameter. Expand $\max(x, y) = x\chi(x - y) + y\chi(y - x)$ and define, inductively:

$$ \begin{aligned} \max^{\eta}(x, y) &= \text{def} \quad x \chi^{\eta}(x - y) + y \chi^{\eta}(y - x) \\ \max^{\eta}(x, \text{list}) &= \text{def} \quad yx \chi^{\eta}(x - \max^{\eta}(\text{list})) + y \chi^{\eta}(\max^{\eta}(\text{list}) - x) \end{aligned} $$

Then, operator $\max^{\eta}$ is smooth in its arguments. Next, with reference to Axiom (43), setting

$$ \tilde{\mu}_{f} = \text{def} \max_{x \in \mathcal{L}_{f}}^{\eta} \tilde{\mu}_{x} - \max_{x \in \mathcal{L}_{f} \downarrow \mathcal{M}_{K}}^{\eta} \tilde{\mu}_{x}, \quad \text{where} \quad \tilde{\mu}_{x} = \text{def} n_{x} + \mu_{x} + \mu_{f_{x}}, $$

the mapping

$$ (f, x) \mapsto (\tilde{\mu}_{f}, \mu_{x}) $$

maps the variable-complete matching $\mathcal{M}_{K}$ for system $A_{\varepsilon}=0$ to a variable-complete matching $\mathcal{M}_{\mu}$ for the system

$$ \forall f \in A_{\varepsilon} : \tilde{\mu}_{f} = 0 \qquad (66) $$

whose dependent variables are the $(\mu_{x})_{x \in X}$. Hence, system (66) is structurally nonsingular, let $\mu_{x}^{\eta}$ denote its generically unique solution.

Let $\eta_{n} \searrow 0$. The sequence $n \mapsto (\mu_{x}^{\eta_{n}})_{x \in X}$ is bounded. It thus possesses converging subsequences. By abuse of notation we rename $n \mapsto (\mu_{x}^{\eta_{n}})_{x \in X}$ one of them, and we denote by $(\mu_{x}^{*})_{x \in X}$ its limit. It is a solution of system (49). This proves the existence of solutions to this system. Uniqueness, however, cannot be guaranteed at this point.
**Proof of Statement 2 of Theorem 1:** We decompose it into its successive substatements:

a) Array $A_{\epsilon}^{\downarrow}$ involves only variables belonging to $\mathcal{X}_{A_{\epsilon}}$;

b) Array $A_{0}^{\downarrow}$ is structurally nonsingular;

c) $A_{0}^{\downarrow}$ is the unique limit of $A_{\epsilon}^{\downarrow}$ when $\epsilon \searrow 0$; and

d) The restart system $A_{\star}$ uniquely determines the states $(y^m)^{+}$ for the hot restart of the new mode. $\square$

*Proof of Substatement a):* This substatement deals with rescaling offsets of variables. Let $x$ be as in Euler identity (30), and select $z$ such that $n =_{\text{def}} m_x - m_z = \min(\mu_x, m_x)$. Using (30) yields:

$$x^{\downarrow} = \varepsilon^{\mu_x} x = \varepsilon^{\mu_x - n} \underbrace{\sum_{i=0}^{n} \binom{n}{i} (-1)^i z^{\bullet(-i)}}_{z_*} .$$

By Lemma 9, the following two cases occur:

* Case $\mu_x \le m$: we have

- $x^{\downarrow} = z_{*}$ and, in the expansion of $z_{*}$, $\varepsilon$ does not occur and no variable $z_i$ is impulsive;

- by condition (53), all the variables $z_i$ occurring in the expansion of $z$ sit within $\mathcal{X}_A$.

* Case $\mu_x > m$: we have

- $x^{\downarrow} = z_{*}^{\downarrow} = \varepsilon^{\mu_x - m} z_{*}$ is not impulsive;

- for $1 \le i \le n$, $\varepsilon^{\mu_x - m} z_i = \varepsilon^{\mu_i} z_i^{\downarrow}$ for $\mu_i = \mu_x - m + \mu_{z_i}$;

- by condition (53), all the variables occurring in the expansion of $z$ sit within $\mathcal{X}_A$.

This analysis is summarized by the following statement:

In the two cases, $x^{\downarrow}$ expands as a linear combination of terms $\varepsilon^{\mu_i} z_i^{\downarrow}$, where $\mu_i \ge 0$ and variable $z_i^{\downarrow}$ sits within array $A_{\varepsilon}$. (67)

This proves Substatement a).

*Proof of Substatement b), the crux of the proof (guided by Rule 1):* Let

$$f = \sum_{x \in \mathcal{L}_f} \varepsilon^{-n_x} x f_x$$

be as in (41) and such that $\mu_f < \infty$. Using axioms (42-44), we get $\mu_{f_x} = 0$, implying $f_x^{\downarrow} = f_x$. Setting $\tilde{\mu}_x =_{\text{def}} n_x + \mu_x + \mu_{f_x} = n_x + \mu_x$, rescaling (45) for $f$ yields:

$$\begin{align*}
f^{\downarrow} &= \varepsilon^{\mu_f} \sum_{x \in \mathcal{L}_f} \varepsilon^{-n_x} x f_x \\
&= \varepsilon^{\mu_f} \sum_{x \in \mathcal{L}_f} (\varepsilon^{-n_x} (\varepsilon^{-\mu_x} x^{\downarrow}) f_x) \\
&= \varepsilon^{\mu_f} \sum_{x \in \mathcal{L}_f} (\varepsilon^{-\tilde{\mu}_x} x^{\downarrow} f_x) \\
&= \sum_{x \in \mathcal{L}_f} (\varepsilon^{\mu_f - \tilde{\mu}_x} x^{\downarrow} f_x).
\end{align*}$$

Using Axiom (43) and equation (49), we get $\mu_f - \tilde{\mu}_x \ge 0$, with equality for the term in the sum involving the unique variable $x$ such that $(f, x) \in \mathcal{M}_K$. This shows that

$$\begin{align}
&\text{when setting } \varepsilon := 0 \text{ in the rescaled array } A_{\varepsilon}^{\downarrow}, \text{ the presence of term } x^{\downarrow} f_x \tag{68} \\
&\text{remains guaranteed for } (f, x) \in \mathcal{M}_K. \nonumber
\end{align}$$
Hence, rescaling $A_\varepsilon$ and then setting $\varepsilon := 0$ ensures that $M_0^\downarrow$ is a variable-complete matching for $A_0^\downarrow$, where $M_0^\downarrow =_{\text{def}} \{(f^\downarrow, x^\downarrow) | (f, x) \in \mathcal{M}_K\}$. This proves Substatement (b).

*Proof of Substatement c):* It is a direct consequence of Substatement b), by Proposition 2 in Section 5.3.

*Proof of Substatement d):* Substatements a,b,c) together imply the following, for array $A_0^\downarrow$:

(i) it is structurally nonsingular;

(ii) it involves only variables belonging to $X_{A_\varepsilon}$;

(iii) when constructing it, no equation belonging to the tail of $A_\varepsilon$ was disabled.

Then, $A_\star$ is obtained from $A_0^\downarrow$ via renaming (50). By (ii) and Lemma 4, the renaming (50) is total and injective. By (iii), all the consistency constraints of the new mode are satisfied by the solution of $A_\star = 0$. Since the renaming (50) is injective by (i), $A_\star$ is also structurally nonsingular. This proves Substatement (d) and finishes the proof of Statement 2 of Theorem 1.

**Proof of Statement 3 of Theorem 1:** We refine Statement 3 as the following lemma:

**Lemma 10**

1. For every variable $x \in X_{A_\varepsilon}^{\text{dep}}$, its rescaling offset $\mu_x$ is independent from the particular choice for the pair $(K, \mathcal{M}_K)$ in the solution of Problem 2 satisfying goodness Conditions (51-53);

2. For $f \in A^{\text{ena}}$, its rescaling offset $\mu_f$ is independent from the particular choice for the pair $(K, \mathcal{M}_K)$ in the solution of Problem 2 satisfying goodness Conditions (51-53);

3. The height function $K: f \mapsto K_f$, for $f \in F$, is unique.

We prove the three statements successively. For the first two, since the variable height $K$ of the array $A_\varepsilon$ is not fixed, we extend rescaling offsets beyond array $A_\varepsilon$, by setting $\mu_x = 0$ for $x \notin X_{A_\varepsilon}^{\text{dep}}$ and $\mu_f = 0$ for $f \notin A_\varepsilon$.

Consider Statement 1 of Lemma 10. With the above convention, denote by $\mu_x^{K, \mathcal{M}}$ the rescaling offset of $x$ for a pair $(K, \mathcal{M}_K)$ being part of a solution of Problem 2 satisfying goodness Conditions (51-53). Consider two such different pairs $(K_1, \mathcal{M}_1)$ and $(K_2, \mathcal{M}_2)$. Suppose $\mu_1 =_{\text{def}} \mu_x^{K_1, \mathcal{M}_1} < \mu_2 =_{\text{def}} \mu_x^{K_2, \mathcal{M}_2}$ for two such pairs $(K_i, \mathcal{M}_i)$, $i=1,2$. Let $x_i^\downarrow = \varepsilon^{\mu_i} x$, $i=1,2$ be the rescaling of $x$ based on the first and second solution. Then, $x_2^\downarrow = \varepsilon^{\mu_2} x = \varepsilon^{\mu_2 - \mu_1} x_1^\downarrow$. Therefore, setting $\varepsilon := 0$ will zero $x_2^\downarrow$ since $\mu_2 - \mu_1 > 0$, thus erasing $x^\downarrow$ from the matching $\mathcal{M}_2$. This contradicts Statement 2 of Theorem 1.

Statement 2 of Lemma 10 is proved by using a similar argument. However, the condition that $f \in A^{\text{ena}}$ is used to ensure that $f \in \mathcal{M}$ whatever the matching selected in the considered solution of Problem 2. Let $\mu_f^{K, \mathcal{M}}$ be the rescaling offset of $f$ for a pair $(K, \mathcal{M}_K)$ being part of a solution of Problem 2 satisfying goodness Conditions (51-53). Suppose $\mu_1 =_{\text{def}} \mu_f^{K_1, \mathcal{M}_1} < \mu_2 =_{\text{def}} \mu_f^{K_2, \mathcal{M}_2}$ for two such pairs $(K_i, \mathcal{M}_{\varepsilon i})$, $i=1,2$. Let $f_1^\downarrow = \varepsilon^{\mu_1} f$ be the rescaling of $f$ based on the first solution. Then, using the same notation, $f_2^\downarrow = \varepsilon^{\mu_2} f = \varepsilon^{\mu_2 - \mu_1} f_1^\downarrow$. Therefore, setting $\varepsilon := 0$ will zero $f_2^\downarrow$ since $\mu_2 - \mu_1 > 0$, thus making the corresponding equation a tautology. This, again, contradicts Statement 2 of Theorem 1.

We also prove Statement 3 of Lemma 10 by contradiction. If $K_1$ and $K_2$ are two variable heights, seen as functions, then $K_1 \neq K_2$ iff either $K_1 \neq \max(K_1, K_2)$ holds, or the opposite holds. Thus, it is enough to investigate one of the latter cases. Let $(K_1, \mathcal{M}_1)$ and $(K_2, \mathcal{M}_2)$ be two solutions such that $K_1 \leq K_2$, $K_1 \neq K_2$, and denote by $\mu_x^1$ and $\mu_x^2$ the rescaling offset of variable $x$ for solutions 1 and 2. Extend $(K_1, \mathcal{M}_1)$ until $K_2$ by adding, to $\mathcal{M}_1$, the set of pairs

$$ \{(f^{\downarrow c_f \bullet k}, x^{\downarrow d_x \bullet k}) \mid (f, x) \in \mathcal{M}, K_{1,f} < k \leq K_{2,f}\}. $$
Denote by $(K_2, \overline{M}_1)$ this extension, and, as before, we extend rescaling offsets by setting

$$ \mu_x^1 = 0 \text{ if } x \text{ sits between } K_1 \text{ and } K_2. \qquad (69) $$

By statement 1 of Lemma 10, $\mu_x$ is the same for both solutions:

$$ \mu_x^1 = \mu_x^2, \text{ for every dependent variable } x \text{ of } A_{K_2}. \qquad (70) $$

Consider the tail of the mode change array of height $K_2$. The equations of this tail must all be matched in $M_2$. Let $f$ be an equation of the new mode with equation offset $c_f>0$ and let $x$ be the variable matched with $f$ in the matching $M$ found in Problem (17). Then, $f^{(c_f-1)\bullet K_{2,f}}$ is a consistency equation sitting in the tail, hence it is enabled and matched in $M_2$ with the state variable $x^{(d_x-1)\bullet K_{2,f}}$, where $d_x$ is the variable offset of $x$ in Problem (17). Consider the variable $x^{d_x\bullet(K_{2,f}-1)}$ sitting in the previous block. We have $x^{d_x\bullet(K_{2,f}-1)} \sim x^{(d_x-1)\bullet K_{2,f}}$, implying that $x^{d_x\bullet(K_{2,f}-1)}$ cannot be matched with an equation in $M_2$. Therefore, it must be matched in $M_2$ with the Euler identity

$$ 0 = \varepsilon \times x^{d_x \bullet (K_{2,f}-1)} - (x^{(d_x-1) \bullet K_{2,f}} - x^{(d_x-1) \bullet (K_{2,f}-1)}). $$

Setting $z =_{\text{def}} x^{(d_x-1)\bullet(K_{2,f}-1)}$, this Euler identity rewrites

$$ 0 = \varepsilon \times z' - (z \bullet - z), $$

which yields the rescaling equation

$$ \mu_{z'}^2 = 1 + \max(\mu_z^2, \mu_z^2) \geq 1. $$

This contradicts (70) since, by (69), $\mu_{z'}^1 = 0$ should hold because $K_{2,f} - 1 \geq K_{1,f}$. This completes the proof of Lemma 10, and thus proves Statement 3 of Theorem 1.

# 9 Addendum: solving nondeterminism and handling cascaded mode changes

So far we considered mode changes between two successive “long modes” (lasting for a positive duration with the same DAE dynamics). Such mode change may not be sufficiently determined, in which case either the structural analysis of the array or Corollary 1 may provide appropriate diagnosis. In this section we consider active determinization, by which constraints are added to the mode change to determinize it. By doing so, we will also encompass cascades of transient (of zero-duration) modes, with limited objectives, however. We will illustrate this with the cup-and-ball example (4) in the elastic impact case (5).

Most of the previous material is reused with little change, if any. Pointers to **reused** material will be given, and emphasis will be put in **red** on what is new.

## 9.1 Revisiting hot restart

**Definition 10 (mode change revisited)** A mode change is defined as a tuple $(S^-, t_*, T, S)$, where $S^-$ is the previous long mode, $S$ is the new long mode, $T$ is the (optional) hot restart constraint, and $t_* \in \mathbb{R}_+$ is the instant when the cascaded mode change occurs. We assume that:

* $S^-= (F^-, X^-)$ and $S=(F, X)$ are square DAE systems; $S^-$ is index-reduced, i.e., $S^- = S_{\downarrow}^-$; $S$ is completed with its latent equations, i.e., $S = S_{\uparrow}$; and

* $T=(F_T, X_T^- \cup Z \cup X_T)$ is a static system of equations where

$$ X_T^- \subset \bigcup_{k \geq 0} X^{-ik} \quad \text{and} \quad X_T \subset \bigcup_{k \geq 0} X^{ik}, $$

and $Z$ collect auxiliary variables. Sets $X_T^-$ and $Z \cup X_T$ collect the free and dependent variables of $T$, respectively. We require $T$ not to be structurally overdetermined.
Static system T relates variables of the previous long mode to variables of the new long mode, thus contributing to the determination of the hot restart.

Recalling that $X_\uparrow$ denotes the set collecting the leading and state variables of the new long mode, we consider the subset $X^+ \subseteq X_\uparrow$ of state variables of the new mode. The hot restart consists in finding values for all the variables belonging to $X^+$.

**Reused 1** The following material is reused with no change:

* Problem 1 of hot restart, with the exception that “previous/new mode” is replaced by “previous/new long mode”;

* Assumption 1 regarding how the instant $t_*$ of mode change is defined. □

**Syntax for hot restart:** To specify a pair $(T, S)$ collecting the new long mode and its restart transient mode, we propose the following syntax:<fnref n="16" />

$$ \begin{array}{l} \text{if } \gamma \text{ then } S \\ \text{and when } \gamma \text{ then } T \end{array} \qquad (71) $$

## 9.2 Revisiting mode change arrays

**Reused 2** The following material is reused with no change:

* The introduction of Section 5 introducing our model of discrete time for the mode change;

* Section 5.1 regarding the toolbox relating shifting and differentiating. □

The construction of arrays, however, must be revisited, so we detail it carefully. We use the notations from Section 3.2: for $F$ a DAE system, $F_\downarrow$, $F_\uparrow$, and $F_\uparrow$ refer to the “index-reduced”, “consistency”, and “completion” of $F$ once index reduction has been performed. When considering mode changes, we use notations (25,26) for variables and functions, and Notations 5 for shifting transition relations. Convention 4 regarding $\sim$-completion is again in force.

**Mode change array:** To the new long mode $S=(F, X)$ we again associate following (32) the infinite system $A$ of equations obtained by stacking successive shifts of the completion $F_\uparrow$. The set of past variables of $A$ is its subset of variables whose value is determined by the previous mode:

$$ \mathcal{X}^{-} = \{x \in \mathcal{X}_{A} \mid \exists x^{-} \in \bigcup_{k<0} (X^{-})^{\bullet k} \text{ s.t. } x^{-} \sim x\}. \qquad (72) $$

The set of dependent variables of $A$ is

$$ \mathcal{X}_A^{\text{dep}} =_{\text{def}} \mathcal{X}_A \setminus \mathcal{X}^{-}, \qquad (73) $$

where $\mathcal{X}_A$ is the set of all the variables involved in $A$ (both leading and state variables).

**Reused 3** The following material is reused with no change:

* Lemma 3 and Definition 5 regarding the definition of facts;

* Convention 5 still applies, i.e., facts are removed from $A$. □

We address Step 2 of Approach 1 by introducing finite arrays, based on the notion of height. A height for mode $S = (F, X)$ denotes a function

$$ K : F \to \mathbb{N}, \text{ denoted } f \mapsto K_f. \qquad (74) $$

In the following definition, $A$ denotes the infinite array (32), and its set $\mathcal{X}_A^{\text{dep}}$ of dependent variables is defined by (73). The changes with reference to Definition 6 are pointed in red:

<footnote n="16">¹⁶This syntax conforms to Modelica, in which “when $\gamma$” denotes the onset of condition $\gamma$.</footnote>
**Definition 11 (mode change array)** Any height K for mode S=(F, X) defines the K-mode change array

$$A_K = \{ f'^m \bullet k \mid f \in F, 0 \le m \le c_f, 0 \le k \le K_f \} \cup T_K, \text{ where} \quad (75)$$

$T_K$ is equal to T in which each $x'^d \in X_T^-$ is replaced by $x'^d \bullet (-1)$ and each $x'^d \in X_T$ is replaced by $x'^d \bullet K_f$, for $(f, x) \in \mathcal{M}$, $\mathcal{M}$ being the matching found by solving the $\Sigma$-method for S.

The set of dependent variables of $A_K$ is $\mathcal{X}_{A_K}^{\text{dep}} = \mathcal{X}_A^{\text{dep}} \cap \mathcal{X}_{A_K}$. Define the tail of $A_K$ and its set of dependent variables by

$$ 
\begin{aligned}
A_K^{\text{tail}} &= \left\{ f'^m \bullet K_f \mid \begin{array}{l} f \in F \\ 0 \le m \le c_f \end{array} \right\} \cup T_K \\
\mathcal{X}_{A_K}^{\text{tail}} &= \left\{ x'^d \bullet K_f \mid \begin{array}{l} f \in F, (f, x) \in \mathcal{M} \\ 0 \le m \le c_f \\ d = \sigma_{f,x} + m \end{array} \right\} \\
\mathcal{X}_{A_K}^{\text{head}} &= \mathcal{X}_{A_K} \setminus \mathcal{X}_{A_K}^{\text{tail}}
\end{aligned} 
\quad (76) $$

We define the subset of equations that must be enabled by

$$A_K^{\text{ena}} = A_K^{\text{tail}} \cup (A_K \cap \bigcup_{k \ge 0} (F^-)^{\bullet k}) \quad (77)$$

and we set $A_K^{\text{dis}} =_{\text{def}} A_K \setminus A_K^{\text{ena}}$.

**Reused 4** The following material is reused with no change:

* Lemma 4 guaranteeing the soundness of the definition of $T_K$ in (75).

☐

## 9.3 Reusing our previous analysis

**Reused 5** The following material is reused with no change:

* Structural implicit function theorem (Proposition 2 and Rule 1);

* Section 6 on the rescaling analysis of mode change arrays;

* Section 7 collecting the main theorems, and their associated proofs.

☐

Indeed, once mode change arrays are properly adapted (see the modifications in Definition 11), the above listed material can be reused with no change (except for how a mode change is named), since it only takes the mode change array as its input.

## 9.4 Cup-and-ball with elastic impact

We discuss here the cup-and-ball example of Section 2 in the case of elastic impact, see (5).

**A specification where transience is made explicit:** The model (4) is modified as follows:

$$ 
\left\{
\begin{array}{ll}
0 = x'' + \lambda x & (e_1) \\
0 = y'' + \lambda y + g & (e_2) \\
\gamma = [s^- \le 0]; \gamma(0) = F & \\
\text{when } \gamma \text{ then } & 0 = L^2 - (x^2+y^2) \quad (k_0) \quad \text{(hot restart)} \\
& 0 = \lambda + s \quad (k_1) \quad \text{(hot restart)} \\
\text{and} & 0 = \lambda \quad (k_3) \\
& 0 = (L^2 - (x^2+y^2)) - s \quad (k_4)
\end{array}
\right.
\quad (78) $$
Model (78) is obtained from (4) by replacing the long guard “if $\gamma$” by the transient guard “when $\gamma$”. We consider the cascade $\gamma: F \to T \to F$, where restart occurs. We regard the transient mode $\gamma=T$ as contributing, by being the “T” player, to the hot restart of the new long mode $\gamma=F$ following Definition 10. Following the construction (11), the resulting mode change array $A_0$ is

$$A_0 = \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \gamma=F, t_* \\ 0 = y'' + \lambda y + g & (e_2) & \\ 0 = \lambda & (k_3) & \\ 0 = (L^2 - (x^2+y^2)) - s & (k_4) & \\ 0 = L^2 - (x^2+y^2) & (k_0) & \text{(hot restart) fact} \\ 0 = \lambda + s & (k_1) & \text{(hot restart)} \end{array} \right. \qquad (79)$$

Array $A_0$ consists of two blocks. The first block collects the dynamics of the new long mode $\gamma=F$, snapshot at $t_*$. The second block collects the transient dynamics $\gamma=T$ at the same instant $t_*$. Note that no latent equation was added: it makes no sense to add them since a transient mode takes zero duration whereas differentiations tightly relate to time shifts as we have seen in Section 5.1.

Equation (k0) is a fact, since it is the zero-crossing triggering the mode change. The dependent variables of this array are $\lambda, s, x'', y'',$ whereas $x, y, x', y'$ are the past variables. $A_0$ is thus structurally singular with 4 dependent variables and 5 equations. Actually, the subsystem $(k_3, k_4, k_1)$, which involves only the two dependent variables $\lambda, s$, is conflicting. Hence, one of these equations must be disabled, which contradicts the definition (77) for $A_0^{\text{ena}}$. We thus need to consider $A_1$:

$$A_1 = \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \gamma=F, t_* \\ 0 = y'' + \lambda y + g & (e_2) & \\ 0 = \lambda & (k_3) & \text{conflict} \\ 0 = (L^2 - (x^2+y^2)) - s & (k_4) & \\ 0 = \varepsilon \times y'' - (y' - y) & (\mathcal{E}_{y''y'}) & \\ 0 = \varepsilon \times x'' - (x' - x) & (\mathcal{E}_{x''x'}) & \\ 0 = x'' + \lambda \cdot x & (e_1) & \gamma=F, t_* + \varepsilon \\ 0 = y'' + \lambda \cdot y & (e_2) & \\ 0 = \lambda & (k_3) & \text{conflict} \\ 0 = (L^2 - (x^2+y^2) - s) & (k_4) & \\ 0 = L^2 - (x^2+y^2) & (k_0) & \text{(hot restart) fact} \\ 0 = \lambda + s & (k_1) & \text{(hot restart)} \end{array} \right. \qquad (80)$$

We pinpoint in red the conflicting equations. They belong to $A_1^{\text{dis}}$, so we disable them. Then, we show in blue an equation-complete matching. Unfortunately, the dependent variable $y'$ is not matched. Further augmenting the array height will not solve this problem. The hot restart specified by model (78) was insufficiently determined. To cope with this, we propose adding the following fully elastic Newton impact law as an additional hot restart constraint following Definition 10:

$$\left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \\ 0 = y'' + \lambda y + g & (e_2) & \\ \gamma = [s^{-} \le 0]; \gamma(0) = F & & \\ \text{when } \gamma \text{ then } & 0 = L^2 - (x^2+y^2) & (k_0) \quad \text{(hot restart)} \\ & 0 = \lambda + s & (k_1) \quad \text{(hot restart)} \\ & 0 = (xx' + yy') + (xx' + yy')^{-} & (\mathcal{N}) \quad \text{(additional hot restart)} \\ \text{and} & 0 = \lambda & (k_3) \\ \text{and} & 0 = (L^2 - (x^2+y^2)) - s & (k_4) \end{array} \right. \qquad (81)$$

With reference to Definition 10 of cascaded mode change and its pair $(T, S)$, equations $(k_0), (k_1), (\mathcal{N})$ constitute the T-part whereas the else branch augmented with the always active equations $(e_1), (e_2)$,
constitutes the S-part. We consider the cascade $\gamma : F \to T \to F$, where restart occurs. Its mode change array $A_1$ with matching $\mathcal{M}$ highlighted in blue is

$$A_1 = \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \gamma=F, t_* \\ 0 = y'' + \lambda y + g & (e_2) & \\ 0 = \lambda & (k_3) & \text{conflict} \\ 0 = (L^2 - (x^2+y^2)) - s & (k_4) & \\ 0 = \varepsilon \times y'' - (y'\cdot - y') & (\mathcal{E}_{y''y'\cdot}) & \\ 0 = \varepsilon \times x'' - (x'\cdot - x') & (\mathcal{E}_{x''x'\cdot}) & \\ 0 = x''\cdot + \lambda\cdot x\cdot & (e_1) & \gamma=F, t_*+\varepsilon \\ 0 = y''\cdot + \lambda\cdot y\cdot + g & (e_2) & \\ 0 = \lambda\cdot & (k_3) & \text{conflict} \\ 0 = (L^2 - (x^2+y^2)\cdot) - s\cdot & (k_4) & \\ 0 = L^2 - (x^2+y^2)\cdot & (k_0) & \text{fact} \\ 0 = \lambda\cdot + s\cdot & (k_1) & \\ 0 = (xx'+yy')\cdot + (xx'+yy')\cdot(-1) & (\mathcal{N}) & \end{array} \right. \quad (82)$$

The three equations constituting the “when” branch sit at the end of the array, with the due replacements following (36). Note that $(k_0)$ is a fact. The four equations constituting the “else” branch occur at both instants $t_*$ and $t_*+\varepsilon$. The needed Euler identities are included. The two equations $(k_3)$, $(k_3^\bullet)$ are conflicting with the rest of the array: since disabling them is allowed, they get disabled. A perfect matching is shown, showing structural nonsingularity of the black part of the array. The useful subarray (obtained by discarding leading equations of the last block) is:

$$A_1 = \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \gamma=F, t_* \\ 0 = y'' + \lambda y + g & (e_2) & \\ 0 = (L^2 - (x^2+y^2)) - s & (k_4) & \\ 0 = \varepsilon \times y'' - (y'\cdot - y') & (\mathcal{E}_{y''y'\cdot}) & \\ 0 = \varepsilon \times x'' - (x'\cdot - x') & (\mathcal{E}_{x''x'\cdot}) & \\ 0 = (xx'+yy')\cdot + (xx'+yy')\cdot(-1) & (\mathcal{N}) & \end{array} \right. \quad (83)$$

The corresponding rescaling analysis is:

$$\begin{align} \mu_{x''} &= \mu_{e_1} &&= \max(\mu_{x''}, \mu_{\lambda}) \nonumber \\ \mu_{\lambda} &= \mu_{e_2} &&= \max(\mu_{y''}, \mu_{\lambda}) \nonumber \\ \mu_s &= \mu_{k_4} \nonumber \\ \mu_{y''} &= \mu_{\varepsilon_{y''y'\cdot}} &&= \max(\mu_{y''}, 1 + \mu_{y'\cdot}) \nonumber \\ 1 + \mu_{x'\cdot} &= \mu_{\varepsilon_{x''x'\cdot}} &&= \max(\mu_{x''}, 1 + \mu_{x'\cdot}) \nonumber \\ \mu_{y'\cdot} &= \mu_{\mathcal{N}} &&= \max(\mu_{y'\cdot}, \mu_{x'\cdot}) \tag{84} \end{align}$$

Its good solution is:

$$\mu_{x''} = \mu_{\lambda} = \mu_{y''} = 1 = \mu_{\varepsilon_{y''y'\cdot}} = \mu_{\varepsilon_{x''x'\cdot}} = \mu_{e_2} = \mu_{e_1}$$

whereas other variables and functions have rescaling offset zero. Rescaling variables and equations using (7) yields

$$\begin{align*} x''\downarrow &= \text{def } \varepsilon \times x'' = x'\cdot - x' \quad ; \qquad e_1^\downarrow &= \text{def } \varepsilon \times e_1 \\ y''\downarrow &= \text{def } \varepsilon \times y'' = y'\cdot - y' \quad ; \qquad \mathcal{E}_{y''y'\cdot}^\downarrow &= \text{def } \varepsilon \times \mathcal{E}_{y''y'\cdot} \\ \lambda\downarrow &= \text{def } \varepsilon \times \lambda \quad ; \qquad e_2^\downarrow &= \text{def } \varepsilon \times e_2 \\ & \qquad \mathcal{E}_{x''x'\cdot}^\downarrow &= \text{def } \varepsilon \times \mathcal{E}_{x''x'\cdot} \end{align*}$$
The rescaled array $A_1^↓$ is then:

$$A_1^↓: \left\{ \begin{array}{l} 0 = x''^↓ + \lambda^↓ x \\ 0 = y''^↓ + \lambda^↓ y \\ 0 = (L^2 - (x^2+y^2)) - s \\ 0 = x''^↓ - (x'^{\bullet} - x') \\ 0 = y''^↓ - (y'^{\bullet} - y') \\ 0 = (x^{\bullet}x'^{\bullet} + y^{\bullet}y'^{\bullet}) + (xx' + yy')^{\bullet}(-1) \end{array} \right.$$

and the resulting (expected) restart system is obtained by renaming as in (50):

$$\text{restart :} \quad \left\{ \begin{array}{l} 0 = x''^↓ + \lambda^↓ x^- \\ 0 = y''^↓ + \lambda^↓ y^- \\ 0 = (L^2 - (x^2+y^2)^-) - s^+ \\ 0 = x''^↓ - (x'^+ - x'^-) \\ 0 = y''^↓ - (y'^+ - y'^-) \\ 0 = (x^- x'^+ + y^- y'^+) + (xx' + yy')^- \end{array} \right. \qquad (85)$$

**An inconsistent specification:** By guarding the mode $\gamma=T$ with when only, model (81) directly specifies that the mode $\gamma=T$ is transient. One could, alternatively, consider the following model:

$$\left\{
\begin{array}{ll}
0 = x'' + \lambda x & (e_1) \\
0 = y'' + \lambda y + g & (e_2) \\
\gamma = [s^- \le 0]; \gamma(0) = F & \\
\text{if } \gamma \text{ then } 0 = L^2 - (x^2+y^2) & (k_0) \\
\text{and } 0 = \lambda + s & (k_1) \\
\text{and when } \gamma \text{ then } 0 = (xx' + yy') + (xx' + yy')^- & (\mathcal{N}) \\
\text{else } 0 = \lambda & (k_3) \\
\text{and } 0 = (L^2 - (x^2+y^2)) - s & (k_4)
\end{array}
\right.
\qquad (86)$$

With reference to model (81), in (86) equations $(k_0)$, $(k_1)$ were shifted, from the $T$ part of the new mode to its $S$ part, see Definition 10. Since mode $\gamma=T$ is possibly long, we must consider the mode change $\gamma : F \to T$. Also, following (32), latent equations must be added to (82) when forming the array:

$$A_1 = \left\{
\begin{array}{lll}
0 = x'' + \lambda x & (e_1) & \gamma=T, t_* \\
0 = y'' + \lambda y + g & (e_2) & \\
0 = L^2 - (x^2+y^2) & (k_4) & \text{fact} \\
0 = xx' + yy' & (k'_4) & \text{conflict} \\
0 = xx'' + x'^2 + y'^2 + yy'' & (k''_4) & \text{conflict} \\
0 = \varepsilon \times y'' - (y'^{\bullet} - y') & (\mathcal{E}_{y''y'^{\bullet}}) & \\
0 = \varepsilon \times x'' - (x'^{\bullet} - x') & (\mathcal{E}_{x''x'^{\bullet}}) & \\
0 = L^2 - (x^2+y^2)^{\bullet} & (k'_4) & \text{fact} \\
0 = (xx' + yy')^{\bullet} & (k''_4) & \text{conflict} \\
0 = (xx' + yy')^{\bullet} + (xx' + yy')^{\bullet}(-1) & (\mathcal{N}) & t_* + \varepsilon
\end{array}
\right.
\qquad (87)$$

Note that equations $(k'_4)$ and $(\mathcal{N})$ are in conflict. Solving this conflict requires disabling one of them, which is forbidden since both belong to the set of must-be-enabled equations by (77). Augmenting the array will not change this. The conclusion is that model (86) is inconsistent due to the conjunction of “if $\gamma$” and “when $\gamma$”. This is indeed expected from physical standpoint: one cannot have both elastic impact $(\mathcal{N})$ and mode $\gamma=T$ being long, meaning inelastic impact.
**A redundant specification of the impact law:** What happens, however, if we explicitly add as (N) a law modeling inelastic impact?

$$
\left\{
\begin{array}{lll}
0 = x'' + \lambda x & (e_1) & \\
0 = y'' + \lambda y + g & (e_2) & \\
\gamma = [s^{-} \le 0]; \gamma(0) = F & & \\
\text{if } \gamma \text{ then } & 0 = L^2 - (x^2+y^2) & (k_0) \\
\text{and } & 0 = \lambda + s & (k_1) \\
\text{and when } \gamma \text{ then } & 0 = (xx'+yy') & (\mathcal{N}) \\
\text{else } & 0 = \lambda & (k_3) \\
\text{and } & 0 = (L^2 - (x^2+y^2)) - s & (k_4)
\end{array}
\right.
\qquad (88)
$$

Then array $A_1$ of (87) is modified as follows in its last two equations:

$$
A_1 =
\left\{
\begin{array}{lll}
0 = x'' + \lambda x & (e_1) & \gamma=T, t_* \\
0 = y'' + \lambda y + g & (e_2) & \\
0 = L^2 - (x^2+y^2) & (k_4) & \text{fact} \\
0 = xx' + yy' & (k'_4) & \text{conflict} \\
0 = xx'' + x'^2 + y'^2 + yy'' & (k''_4) & \text{conflict} \\
0 = \varepsilon \times y'' - (y'^{-} - y') & (\varepsilon_{y''} y'^{-}) & \\
0 = \varepsilon \times x'' - (x'^{-} - x') & (\varepsilon_{x''} x'^{-}) & \\
0 = L^2 - (x^2+y^2)^{-} & (k'_4) & \text{fact} \\
0 = (xx' + yy')^{-} & (k''_4) & \\
0 = (xx' + yy')^{-} & (\mathcal{N}) & t_*+\varepsilon
\end{array}
\right.
\qquad (89)
$$

The situation is now different: equations $(k_4')$ and $(\mathcal{N})$ are syntactically identical, showing redundancy. So we could fuse them by discarding $(\mathcal{N})$, which brings us back to the inelastic impact case handled in Section 2. This is a clear and easy case.

However, depending on how the impact law is formulated, it could be redundant with the first latent equation, but still different from it in its syntax. Checking redundancy could be performed at run time: $(\mathcal{N})$ is discarded at compile time and evaluated at run time to check that $(k_4')$ and $(\mathcal{N})$ take equal values. Performing the same checking on array (87) would violate redundancy checking.

9.5 On the modeling power of our approach

To understand the power and limits of our approach, we consider the Newton cradle displayed in Fig. 5.

We assume elastic impact. The transient modes are the impacts. Long modes inherit an index two from the pendulum in Cartesian coordinates. In performing the modeling of this device, the difficult part is the model of the impacts, see, e.g., the interesting discussion in Section 4 of [23]. The case of two balls is easy since the two tangent velocities resulting from the impact are determined by using two laws, namely the conservation of momentum and the conservation of kinetic energy. The resulting model can then be handled as the cup-and-ball in the elastic impact case. The modeling is more subtle if a larger number of balls is considered. Using only the above two laws results in nondeterminism. The analysis must be complemented by considering the potential energy associated to deformations and its propagation through waves. To circumvent such difficulties, several authors determinize their model by using smooth but stiff versions of the contact—the resulting simulations highly vary with how this smoothing is performed [23].

The impact modeling problem is indeed by itself a subject matter [19]. Clearly, our approach provides no help for discovering the impact model. We only offer means to specify the impact model once it is known. Admittedly, this impairs modularity of modeling, since one can hardly expect being able to specify a single pendulum with its general context, and combine such models to get a model for the Newton cradle.
<figure><img src="image_11.png" /><figcaption>Figure 5: Newton cradle. Source: https://fr.wikipedia.org/wiki/Pendule_de_Newton.</figcaption></figure>

# 10 Conclusion and perspectives

We contributed to the mathematical soundness of DAE-based modeling languages such as Modelica, by formally solving the hot restart problem for a multimode DAE. We targeted modeling languages that are physics agnostic. Accordingly, our approach is physics agnostic. Despite being physics agnostic, our solution for the hot restart preserves (possibly latent) invariants.

We handle impulsive behaviors with no restriction. Our approach generates hot restart for general nonlinear DAE provided that the non-polynomial part does not interfere with the impulsive part—this condition is checked by our algorithm.<fnref n="17" />

Our approach was first developed for mode changes involving two successive long modes. We extended it to handle the hot restart for *finite cascades of transient mode changes*, such as occurring in the cup-and-ball with elastic impact.<fnref n="18" /> This was achieved by introducing additional hot restart constraints, and then reusing with little changes the analysis developed in Sections 4–8 for mode changes separating long modes.

Our approach builds on the principles of *structural analysis* abstracting numerical equations as their incidence graphs. Our algorithm for synthesizing the hot restart bears similarities with the Σ-method for the index reduction of (single mode) DAE systems.

So far we solved the hot restart problem for a mode change considered in isolation. However, a multimode DAE system has a large number of modes, its number of mode changes is even larger. Enumerating mode changes is the beyond reach. One way of overcoming this consists in reducing the set of mode changes by specifying assertions in the form of invariants. A complementary way is to reuse the techniques proposed in [6] for the all-modes-at-once structural analysis of multimode DAE systems.

We proved in [11, 7, 10] that a notion of structural interface can be associated to any DAE

<footnote n="17">¹⁷The class of multimode DAE for which hot restart is generated can be widened. By allowing for rational offsets, our rescaling analysis can be extended to handle the partition polynomial/nonpolynomial (instead of linear/nonlinear). Variables inheriting a non-integer offset can still be handled, provided that they do not need to be integrated or differentiated using Euler identities. It is not clear how useful this widening would be in practice.</footnote>

<footnote n="18">¹⁸This mode change takes the form free-motion → straight-rope → free-motion, with zero time spent in the “straight-rope” transient mode.</footnote>
system having more dependent variables than equations — the mathematical model for a model class
in Modelica. With this notion of structural interface, index reduction can be performed in a modular
way (i.e., at the class level) instead of globally — as performed in DAE based tools today. This
allows an impressive scaling-up for the compilation of Modelica models [10], by handling models with
millions of equations. From the above mentioned similarity, we expect that the structural interface
can be extended to include mode change information. This is in contrast with the computational
cost of methods using the Quasi-Weierstrass decomposition, not to mention the use of elimination
which becomes totally prohibitive when polynomial systems are considered. The way forward is
clear: modular methods must be developed for the hot restart of mode changes.
## References

[1] Mario Barela. *A Complementarity Approach to Modeling Dynamic Electric Circuits*. PhD thesis, University of Iowa, 2016.

[2] Albert Benveniste, Benoît Caillaud, Hilding Elmqvist, Khalil Ghorbal, Martin Otter, and Marc Pouzet. Multi-mode DAE models - challenges, theory and implementation. In Bernhard Steffen and Gerhard J. Woeginger, editors, *Comp. Software Sci. – SOTA Persp.*, volume 10000 of *Lecture Notes in Comp. Sci.*, pages 283–310. Springer, 2019.

[3] Albert Benveniste, Benoît Caillaud, and Mathias Malandain. The mathematical foundations of physical systems modeling languages. *Annual Reviews in Control*, 50:72–118, 2020.

[4] Albert Benveniste, Benoît Caillaud, and Mathias Malandain. The Mathematical Foundations of Physical Systems Modeling Languages. Research Report RR-9334, Inria, April 2020.

[5] Albert Benveniste, Benoît Caillaud, and Mathias Malandain. Structural Analysis of Multimode DAE Systems: summary of results. Research Report RR-9387, Inria Rennes – Bretagne Atlantique, January 2021.

[6] Albert Benveniste, Benoît Caillaud, Mathias Malandain, and Joan Thibault. Algorithms for the structural analysis of multimode Modelica models. *Electronics*, 11(17), 2022.

[7] Albert Benveniste, Benoît Caillaud, Mathias Malandain, and Joan Thibault. Towards the separate compilation of Modelica: Modularity and interfaces for the index reduction of incomplete DAE systems. In *Proc. 15th Int. Modelica Conf.*, volume 204, page 10, Aachen, Germany, October 2023.

[8] Bernard Brogliato. *Nonsmooth Mechanics: Models, Dynamics and Control*. Communications and Control Engineering. Springer London, 2012.

[9] Peter Bunus and Peter Fritzson. Methods for structural analysis and debugging of modelica models. *Proc. 2nd Int. Modelica Conf.*, 05 2002.

[10] Benoît Caillaud, Albert Benveniste, and Mathias Malandain. Benchmarking the Modular Structural Analysis Algorithm. In *2025 - 16th International Modelica & FMI Conference*, pages 1–14, Lucerne, Switzerland, September 2025. Ulf Christian Müller and Dirk Zimmer, Linköping University Press.

[11] Benoît Caillaud, Mathias Malandain, and Joan Thibault. Implicit structural analysis of multimode DAE systems. In *23rd ACM Int. Con. on Hybrid Syst.: Comp. and Control (HSCC 2020)*, Sydney, Australia, April 2020.

[12] Stephen L. Campbell and C. William Gear. The index of general nonlinear DAEs. *Numerische Mathematik*, 72:173–196, 1995.

[13] Yahao Chen and Stephan Trenn. On impulse-free solutions and stability of switched nonlinear differential-algebraic equations. *Automatica*, 156(111208):1–14, 2023.

[14] David A. Cox, John Little, and Donal O’Shea. *Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra*. Springer Publishing Company, Incorporated, 3rd edition, 2010.

[15] Jean Dieudonné. *Fondements de l'analyse moderne*. Gauthier-Villars, 1963.

[16] David L. Dill. Timing assumptions and verification of finite-state concurrent systems. In Joseph Sifakis, editor, *Automatic Verification Methods for Finite State Systems, International Workshop, Grenoble, France, June 12-14, 1989, Proceedings*, volume 407 of *Lecture Notes in Computer Science*, pages 197–212. Springer, 1989.
[17] Hilding Elmqvist, Sven Erik Mattsson, and Martin Otter. Modelica extensions for multi-mode DAE systems. In Hubertus Tummescheit and Karl-Erik Arzen, editors, *Proc. 10th Int. Modelica Conf.*, Lund, Sweden, September 2014. Modelica Association.

[18] Charles William Gear, Ben Leimkuhler, and Gopal K Gupta. Automatic integration of Euler-Lagrange equations with constraints. *J. Comput. Appl. Math.*, 12:77-90, 1985.

[19] G. Gilardi and I. Sharf. Literature survey of contact dynamics modelling. *Mechanism and Machine Theory*, 37(10):1213-1239, 2002.

[20] Peter Hamann and Volker Mehrmann. Numerical solution of hybrid systems of differential-algebraic equations. *Comput. Methods Appl. Mech. Eng.*, 197(6):693-705, 2008.

[21] Maurice Heemels, Kanat Camlibel, and J. M. Schumacher. On the dynamic analysis of piecewise-linear networks. *IEEE Trans. Circuits Syst. I: Fundam. Theory Appl.*, 49(3):315-327, Mar 2002.

[22] Rukhsana Kausar and Stephan Trenn. Impulses in structured nonlinear switched DAEs. In *Proc. 56th IEEE Conf. Decis. Control*, pages 3181 – 3186, Melbourne, Australia, 2017.

[23] Edward A. Lee. Fundamental limits of cyber-physical systems modeling. *ACM Trans. Cyber Phys. Syst.*, 1(1):3:1-3:26, 2016.

[24] Daniel Liberzon and Stephan Trenn. Switched nonlinear differential algebraic equations: Solution theory, Lyapunov functions, and stability. *Automatica*, 48(5):954-963, 2012.

[25] Ton Lindstrøm. An invitation to nonstandard analysis. In Nigel J. Cutland, editor, *Nonstandard Analysis and its Applications*, pages 1-105. Cambridge Univ. Press, 1988.

[26] Sven Erik Mattsson, Martin Otter, and Hilding Elmqvist. Modelica hybrid modeling and efficient simulation. In *Proc. 38th IEEE Conf. on Decis. Control*, pages 3502-3507, 1999.

[27] Sven Erik Mattsson, Martin Otter, and Hilding Elmqvist. Multi-Mode DAE systems with varying index. In Hilding Elmqvist and Peter Fritzson, editors, *Proc. 11th Int. Modelica Conf.*, Versailles, France, September 2015. Modelica Association.

[28] Sven Erik Mattsson and Gustaf Söderlind. Index reduction in differential-algebraic equations using dummy derivatives. *SIAM J. Sci. Comput.*, 14(3):677-692, 1993.

[29] Antoine Miné. A new numerical abstract domain based on difference-bound matrices. In Olivier Danvy and Andrzej Filinski, editors, *Programs as Data Objects, Second Symposium, PADO 2001, Aarhus, Denmark, May 21-23, 2001, Proceedings*, volume 2053 of *Lecture Notes in Computer Science*, pages 155-172. Springer, 2001.

[30] Martin Otter and Dirk Zimmer, editors. *Proceedings of 9th International Modelica Conference*. Munich, Germany, September 2012.

[31] Constantinos C. Pantelides. The consistent initialization of differential-algebraic systems. *SIAM J. Sci. Stat. Comput.*, 9(2):213-231, 1988.

[32] Friedrich Pfeiffer. On non-smooth multibody dynamics. *Proc. Inst. Mech. Eng., Part K: J. Multi-body Dyn.*, 226(2):147-177, 2012.

[33] Friedrich Pfeiffer and Christoph Glocker. *Multibody Dynamics with Unilateral Contacts*. Wiley, 2008.

[34] Alex Pothen and Chin-Ju Fan. Computing the block triangular form of a sparse matrix. *ACM Trans. Math. Softw.*, 16(4):303-324, 1990.
[35] John D. Pryce. A simple structural analysis method for DAEs. *BIT Numer. Math.*, 41(2):364-394, 2001.

[36] Svenja Schoeder, Heinz Ulbrich, and Thorsten Schindler. Discussion of the Gear-Gupta-Leimkuhler methods for impacting mechanical systems. In *Multibody Syst. Dyn.*, volume 31, pages 477-495, 2013.

[37] Jean U. Thoma. *Introduction to Bond Graphs and Their Applications*. Pergamon Int. Lib. Sci., Technol., Eng. Soc. Stud. Pergamon Press, 1975.

[38] Stephan Trenn. *Distributional Differential Algebraic Equations*. PhD thesis, Technischen Universität Ilmenau, 2009.

[39] Stephan Trenn. Regularity of distributional differential algebraic equations. *Math. Control Signals Syst.*, 21(3):229-264, 2009.

[40] Arjan van der Schaft and Dimitri Jeltsema. *Port-Hamiltonian systems theory: An introductory overview*. *Foundations and Trends in Systems and Control*, 1(2-3):173-378, 2014.

[41] Dirk Zimmer. *Equation-Based Modeling of Variable-Structure Systems*. PhD thesis, ETH Zürich, 2010.
# Appendices

In these appendices, missing proofs are collected and additional illustrative examples are developed.

## A Proof of Corollary 1

In the proof of Statement 2 of Theorem 1, goodness Conditions (51-53) play a role only in the proof of Substatement d). Substitutions (54), however, restore the validity of Conditions (51-53) for this substatement.

## B Proof of Lemma 8

Writing $\mu_i$ instead of $\mu_{f_i}$ for $i=1, 2$, and using Lemma 7, we have, by setting $\tilde{\mu}_i =_{\text{def}} \max(\mu_1, \mu_2) - \mu_i$:

$$ 
\begin{aligned}
\Lambda(f_1+f_2) &= (f_1+f_2)^{\downarrow}[\varepsilon := 0] \\
&= (\varepsilon^{\mu(f_1+f_2)}(f_1+f_2))[\varepsilon := 0] \\
&= (\varepsilon^{\max(\mu_1,\mu_2)}(f_1+f_2))[\varepsilon := 0] \\
&= (\varepsilon^{\tilde{\mu}_1}\varepsilon^{\mu_1}f_1 + \varepsilon^{\tilde{\mu}_2}\varepsilon^{\mu_2}f_2)[\varepsilon := 0] \\
&= \sum_{i:\mu_i=\max(\mu_1,\mu_2)} \Lambda(f_i),
\end{aligned}
\quad (90) $$

and

$$ 
\begin{aligned}
\Lambda(f_1 \times f_2) &= (f_1 \times f_2)^{\downarrow}[\varepsilon := 0] \\
&= (\varepsilon^{\mu(f_1 \times f_2)}(f_1 \times f_2))[\varepsilon := 0] \\
&= (\varepsilon^{(\mu_1+\mu_2)}(f_1 \times f_2))[\varepsilon := 0] \\
&= \varepsilon^{\mu_1} f_1[\varepsilon := 0] \times \varepsilon^{\mu_2} f_2[\varepsilon := 0] \\
&= \Lambda(f_1) \times \Lambda(f_2).
\end{aligned}
\quad (91) $$

The lemma is proved.

## C Proof of Theorem 2

This proof will require a more in-depth analysis of the structure of matchings $M_\varepsilon$, which in turn requires an extended analysis of Euler identities and their role in matchings. As a consequence, we will need to handle $\sim$-closures with care, thus we will indicate them explicitly (alternative case of Convention 3).

The mode change arrays $A_\varepsilon$ considered in Problem 2 are $\sim$-closed. To further study their matchings, it will be convenient to handle equivalence classes of variables modulo $\sim$. Recalling Definition 3, $\hat{x}$ denotes the equivalence class of $x$ modulo $\sim$. Relation $\sim$ is stable under both differentiation and shifting. Hence we can define

$$ \hat{x}' =_{\text{def}} \hat{x}' \quad \text{and} \quad \hat{x}^{\bullet} =_{\text{def}} \hat{x}^{\bullet} $$

Let $S = (F, X)$, $X \subseteq \Xi'^{\bullet}$ be a system, and set $\hat{X} =_{\text{def}} \{\hat{x} \mid x \in X\}$. The incidence graph $G_S = (F \cup X, E)$ of $S$ was introduced in Section 3.1. We will also consider the $\sim$-quotient of $G_S$, which is a bipartite graph $\hat{G}_S = (F \cup \hat{X}, \hat{E})$ over the pair $(F, \hat{X})$, defined as follows:

$$ (f, \hat{x}) \in \hat{E} \quad \text{iff} \quad \exists x \in \hat{x} : (f, x) \in E. \quad (92) $$

The whole apparatus of bipartite graphs applies, which includes in particular Section 3.1 regarding matchings. Consider both

* the $\sim$-closure $\tilde{S}$ of $S$ and its incidence graph $G_{\tilde{S}}$, and
• the quotient graph $\hat{\mathcal{G}}_S$ defined in (92).

The following result holds:

**Lemma 11** For $S, \tilde{S}$, and $\hat{\mathcal{G}}_S$ given as above, there exists a total injective function $\Gamma$, mapping any variable-complete (respectively equation-complete) matching $\tilde{\mathcal{M}}$ for $\hat{\mathcal{G}}_S$, to a variable-complete (respectively equation-complete) matching $\tilde{\mathcal{M}} = \Gamma(\tilde{\mathcal{M}})$ for $\tilde{S}$.

*Proof* By definition of $\hat{\mathcal{G}}_S$, for every $(f, \hat{x}) \in \tilde{\mathcal{M}}$, we can select a variable $x \in \hat{x}$ such that $(f, x) \in \mathcal{G}_S$, we denote by $\chi(\hat{x})$ this selection. Perform the following finite recursion:

1. Initialization: $\mathcal{M}_0 = \{(f, x) | x = \chi(\hat{x})\}, X_0 = X$;

2. While $X_n \supset X_{n-1}$, update $(\mathcal{M}_{n+1}, X_{n+1}) \leftarrow (\mathcal{M}_n, X_n)$ by

• adding the pair $(\mathcal{E}_{xz}, z)$ to $\mathcal{M}_n$, for every $x, z \in X_n$ such that $x \in \mathcal{M}_n, z \sim x, z \notin \mathcal{M}_n$ (we use Notations 2);

• adding to $X_n$ every variable $u$ occurring in some above added $\mathcal{E}_{xz}$ and such that $u \notin \mathcal{M}_n, u \neq z$.

Let $(\mathcal{M}_N, X_N)$ be the fixpoint. We set $\tilde{\mathcal{M}} = \Gamma(\tilde{\mathcal{M}}) =_{\text{def}} \mathcal{M}_N$ and we consider the map $\Gamma: \tilde{\mathcal{M}} \mapsto \tilde{\mathcal{M}}$. If $\tilde{\mathcal{M}}$ is variable-complete, then so is $\tilde{\mathcal{M}}$ by construction. Assume next that $\tilde{\mathcal{M}}$ is equation-complete. Then, the above recursion only adds equations that occur in the increasing sequence of matchings. Thus, $\tilde{\mathcal{M}}$ is equation-complete as well. $\square$

**Lemma 11** allows us to reduce the search for a matching over the $\sim$-closure of a system, to the (simpler) search for a matching over its $\sim$-quotient. Let us apply this to the array $(A_K, \mathcal{X}_{A_K})$, where $A_K$ was defined in (36). In particular,

$$ \hat{\Psi}(K, \mathcal{M}) =_{\text{def}} \left\{ (f^m, \hat{x}^d) \middle| \begin{array}{l} (f, x) \in \mathcal{M} \\ 0 \le m < c_f \\ d = d_x - c_f + m \end{array} \right\} \cup \left\{ (f'^{c_f} \bullet k, \hat{x}^{d_x} \bullet k) \middle| \begin{array}{l} (f, x) \in \mathcal{M} \\ 0 \le k \le K_f \end{array} \right\} \quad (93) \\ =_{\text{def}} \tilde{\mathcal{M}}_{\uparrow} \cup (\tilde{\mathcal{M}}_K)_{\downarrow} $$

defines a matching $\hat{\Psi}(K, \mathcal{M}) =_{\text{def}} \tilde{\mathcal{M}}_{\uparrow} \cup (\tilde{\mathcal{M}}_K)_{\downarrow}$, where

• $\tilde{\mathcal{M}}_{\uparrow}$ identifies with $\mathcal{M}_{\uparrow}$, the equation-complete matching (21) associated to consistency conditions in the $\Sigma$-method;

• $(\tilde{\mathcal{M}}_K)_{\downarrow}$ is a variable-complete matching for the quotient incidence graph of the pair $((A_K)_{\downarrow}, (\hat{\mathcal{X}}_{A_K})_{\downarrow})$, where

$$ (A_K)_{\downarrow} =_{\text{def}} \{ f'^{c_f} \bullet k \mid f \in F, 0 \le k \le K_f \}, \text{ and} $$

$$ (\hat{\mathcal{X}}_{A_K})_{\downarrow} =_{\text{def}} \{ \hat{x}^{d_x} \bullet k \mid x \in X, (f, x) \in \mathcal{M}, 0 \le k \le K_f \} $$

collect all the leading equations and the equivalence classes of all the leading variables of the array.

By Lemma 11

$$ \mathcal{M}_{\varepsilon} =_{\text{def}} \Gamma(\hat{\Psi}(K, \mathcal{M})) \quad (94) $$

yields a matching for the $\sim$-closure $A_{\varepsilon}$ of array $A_K$, and the images by $\Gamma$ of $\tilde{\mathcal{M}}_{\uparrow}$ and $(\tilde{\mathcal{M}}_K)_{\downarrow}$ are equation-complete and variable-complete, respectively.

We now consider requirement (47) stating, in particular, that all the equations of the tail block should be enabled by $\mathcal{M}_{\varepsilon}$.
**Lemma 12** To ensure that the above matching $\mathcal{M}_\varepsilon$ for $(A_\varepsilon, \chi_{A_\varepsilon})$ satisfies requirement (47), it is enough to select a variable height $K$ such that

$$ \forall f \in F \Rightarrow K_f \ge c_f. \qquad (95) $$

*Proof* The reader is referred to Problem (17) (index reduction) applied to the new mode $(F, X)$. The lemma is vacuously true if $c_f = 0$ for every $f$, since no consistency equation exists in this case and $K = 0$ is fine. So, we assume that $c_f > 0$ holds for some $f$. Let $f \in F$ be such that $c_f > 0$, and consider

$$ g=f^{\prime m \bullet K_f} \in A_{\varepsilon}^{\text {tail }} \text { with } 0 \leq m< c_f. \qquad (96) $$

Since consistency equations $f^{\prime m \bullet k}, 0 \le k \le K_f$ are all unmatched in $(\widehat{\mathcal{M}}_K)_\downarrow$, so is equation $g$. By (95), $m < c_f \le K_f$ holds. Consider

$$ h =_{\text{def}} f^{\prime c_f \bullet k}, \quad \text{where} \quad k =_{\text{def}} K_f - c_f + m. \qquad (97) $$

Since $c_f+k = c_f+K_f-c_f+m = K_f+m$, we have $h \sim g$. Equation $h$ is matched in $(\widehat{\mathcal{M}}_K)_\downarrow$ since it is a shifted leading equation; it is actually matched with $x^{\prime d_x \bullet k}$. Setting $d=d_x-c_f+m$, we have $d+K_f = d_x-c_f+m+K_f = d_x+k$, hence $x^{\prime d_x \bullet k} \sim x^{\prime d \bullet K_f}$. By Lemma 1, variables $x^{\prime d_x \bullet k}$ and $x^{\prime d \bullet K_f}$ are related by Euler identity $0=\mathcal{E}_{uv}$, where $u=x^{\prime d \bullet K_f}$ and $v=x^{\prime d_x \bullet k}$, and this Euler identity belongs to $A_\varepsilon$ by Lemma 2 since $A_\varepsilon$ is the $\sim$-closure of $A_K$.

The introduction of map $\chi$ is the key step of the proof: Consider the subarray $A_g$ of $A_\varepsilon$:

$$ A_g =_{\text{def}} \left\{ \begin{array}{l} 0 = h(v, V) \\ 0 = \mathcal{E}_{uv}(u, v, U) \end{array} \right. \qquad (98) $$

where $u = x^{\prime d \bullet K_f}, v = x^{\prime d_x \bullet k}, U$ collects variables $\prec u$, and the matchings are indicated in blue. We associate to $A_g$ the array $\chi(A_g)$ defined by

$$ \chi(A_g) =_{\text{def}} \left\{ \begin{array}{l} 0 = g(u, W) \\ 0 = \mathcal{E}_{uv}(u, v, U) \end{array} \right. \qquad (99) $$

Since this change is only local, the remaining part of the matching is not modified. Performing the mapping $A_g \mapsto \chi(A_g)$ for every consistency equation $g$ belonging to the tail of $A_K$, yields a matching satisfying (47). The lemma is proved. $\square$

**Taking past variables into account:** Lemma 12 does not take past variables $\mathcal{X}^-$ into account. The effect of $\mathcal{X}^-$ is twofold. First, some equations become facts, and thus get removed from the array. Second, the past variables are no longer dependent variables, and thus cannot participate in matchings. We investigate these two effects in detail.

**Taking facts into account:** If $K$ is the variable height of the considered mode change array, consistency equations $f^{\prime m \bullet K_f}$ that are facts are removed from the array. To account for this, we modify Lemma 12 by replacing Condition (95) by

$$ \forall f \in F \Rightarrow K_f \ge \max \left\{ c_f - m \mid \begin{array}{l} 0 \le m \le c_f \\ f^{\prime m \bullet K_f} \notin \Phi_{A_\varepsilon} \end{array} \right\} \qquad (100) $$

No other change is needed in the proof.

Lemma 12 with its extension (100) is illustrated in Fig. 2, displaying array $A_1$ for the cup-and-ball example. The exchange map $\chi$ is exemplified by the substitution of $(k_1'')$ for $(k_1')$ when disabling equations.
**Removing past variables from dependent variables:** Consider $g$ and $h$ defined in (96) and (97). We know that $m+K_f = c_f+k$, and $k = m+K_f-c_f \ge m$ holds by (100). Equation $h$ belongs to the $k$-th instant $(F_{\downarrow}^{\bullet k}, X_{\downarrow}^{\bullet k})$ in array $A_K$. Therefore, in any possible choice for $M_\epsilon$, equation $g$ is paired with some variable belonging to $X_{\downarrow}^{\bullet k}$. This leads to considering the set $X_{\downarrow}^g$ collecting all the variables occurring in $g$ that could possibly be paired with $g$ in some matching $M_\epsilon$:

$$X_{\downarrow}^g = \text{def} \left\{ x^{\prime d \bullet K_f} \text{ occurs in } g \mid \begin{array}{l} x^{\prime d_x \bullet k} \in X_{\downarrow}^{\bullet k} \\ d_x + k = d + K_f \end{array} \right\} \quad (101)$$

Let $K$ be a variable height for the array satisfying condition (100). For $f \in F$ such that $c_f>0$, and $0 \le m < c_f$, we consider the following properties, where $g=f^{\prime m \bullet K_f}$:

$$P_*(K, f, m) = \text{def} \text{ either } f^{\prime m \bullet K_f} \text{ is a fact, or } X_{\downarrow}^g \setminus \mathcal{X}^- \neq \emptyset$$

$$P^*(K, f, m) = \text{def} \text{ either } f^{\prime m \bullet K_f} \text{ is a fact, or } X_{\downarrow}^g \cap \mathcal{X}^- = \emptyset$$

so that $P^*(K, f, m) \Rightarrow P_*(K, f, m)$, and we set

$$P_*(K) = \text{def} \quad \bigwedge \{P_*(K, f, m) \mid f \in F, 0 \le m < c_f\}$$

$$P^*(K) = \text{def} \quad \bigwedge \{P^*(K, f, m) \mid f \in F, 0 \le m < c_f\}$$

and, finally:

$$K_* = \text{def} \quad \min\{K \mid K \models (100) \text{ and } K \models P_*(K)\}$$

$$K^* = \text{def} \quad \min\{K \mid K \models (100) \text{ and } K \models P^*(K)\}$$

so that $K_* \le K^*$.

**Lemma 13** With these notations, to find a variable-complete matching $M_\epsilon$ satisfying requirement (47):

1. it is enough to select $K \ge K^*$, whereas

2. no such matching can be expected unless $K \ge K_*$.

*Proof* Let $K$ satisfy (100), and assume the existence of a variable-complete matching $M_\epsilon$ satisfying requirement (47). Assume that $f^{\prime m \bullet K_f}$ is not a fact, and there exists a variable $x^{\prime d \bullet K_f}$ that is paired with $f^{\prime m \bullet K_f}$ in the matching $\hat{\Psi}(K, M)$ defined in (93). The set $\mathcal{X}^-$ of past variables forbids such a pairing if and only if it contains $x^{\prime d \bullet K_f}$. We deduce that

$$\text{there exists a matching } \mathcal{M}_{\downarrow} \text{ following (19), such that} \quad (f^{\prime m \bullet K_f}, x^{\prime d \bullet K_f}) \in \hat{\Psi}(K, M) \text{ and } x^{\prime d \bullet K_f} \notin \mathcal{X}^{-}. \quad (102)$$

This leads to considering the property

$$P(K, f, m) = \text{def} \quad \text{either } f^{\prime m \bullet K_f} \text{ is a fact, or (102) holds.}$$

Property $P(f, m)$ explores all the alternatives in choosing a candidate variable $x^{\prime d \bullet K_f}$, for pairing with $f^{\prime m \bullet K_f}$ in some matching $\hat{\Psi}(K, M)$ considered in (93).

Condition (102) is too complicated, as it involves exploring the set of all the matchings $M$ solution of Problem (17). We thus consider the following properties:

$$P_*(K, f, m) = \text{def} \quad \text{either } f^{\prime m \bullet K_f} \text{ is a fact, or } X_{\downarrow}^g \setminus \mathcal{X}^- \neq \emptyset.$$

$$P^*(K, f, m) = \text{def} \quad \text{either } f^{\prime m \bullet K_f} \text{ is a fact, or } X_{\downarrow}^g \cap \mathcal{X}^- = \emptyset.$$

where $g = f^{\prime m \bullet K_f}$.
Then, the following implications hold:

$$ \mathcal{P}^*(K, f, m) \Rightarrow \mathcal{P}(K, f, m) \Rightarrow \mathcal{P}_*(K, f, m), $$

which proves statements 1 and 2 of the lemma.

<figure><img src="image_12.png" /></figure>

The following result states that there is no point going beyond $K^*$ in searching for a good solution of Problem 2:

**Lemma 14** Let $K$ be the variable height of a good solution of Problem 2. Then, $K \le K^*$ holds.

*Proof* Let $(f, m)$ be a pair such that: 1) $m = c_f - 1$, and 2) $f'^{m} \bullet K_f$ is not a fact and is paired with $x'^{(d_x-1) \bullet K_f}$ in the matching $M_\epsilon$. Then, by (99), the Euler identity

$$ \mathcal{E}_{yz} : \quad \varepsilon \times \underbrace{x'^{(d_x-1) \bullet (K_f-1)}}_{y} = \underbrace{x'^{(d_x-1) \bullet K_f}}_{z} - \underbrace{x'^{(d_x-1) \bullet (K_f-1)}}_{u} $$

needs to be added to the array. Since $\mu_z = 0$ by goodness condition (51), this Euler identity induces the following rescaling equation:

$$ \mu_y = \mu_{\mathcal{E}_{yz}} \ge 1 + \mu_u. \qquad (103) $$

The sequel of the proof is by contradiction: We will prove:

$$ \text{if } K > K^*, \text{ we can find a pair } (f, m) \text{ as above, such that } \mu_u > 0. \qquad (104) $$

Property (104) would imply $\mu_y > 1$ by (103). In this case, goodness condition (53) would get violated for $y = x'^{d_x} \bullet (K_f - 1)$, since $K_f - 1 + \min(d_x, \mu_y) > K_f$.

It thus remains to prove (104): Since $K > K^*$, we can find a pair $(f, m)$ as above, such that $X_{\downarrow}^h \cap \mathcal{X}^- = \emptyset$, where $h = f'^{m} \bullet (K_f - 1)$, and $u \in X_{\downarrow}^h$. Since $u$ is not a past variable, it must be paired in the matching $M_\epsilon$, namely with the Euler identity

$$ \mathcal{E}_{uv} : \quad \varepsilon \times \underbrace{x'^{(d_x-1) \bullet (K_f-1)}}_{u} = \underbrace{x'^{(d_x-2) \bullet K_f}}_{v} - \underbrace{x'^{(d_x-2) \bullet (K_f-1)}}_{w} $$

which induces the following rescaling equation:

$$ \mu_u = \mu_{\mathcal{E}_{uv}} \ge 1 + \mu_w > 0, $$

proving (104). The proof is now complete.

<figure><img src="image_13.png" /></figure>

Lemmas 13 and 14 together prove Theorem 2 and make the heights $K_*$ and $K^*$ precise.

Theorem 2 was illustrated for the cup-and-ball example by the results $K_* = K^* = 1$. We illustrate in Appendix D.1 what happens if we still insist taking $K > 1$ for this example.

# D Examples having no good solution

## D.1 The cup-and-ball with $K > 1$

In Section 2, we successfully studied the cup-and-ball example with a minimal mode change array $A$ of height $K=1$. In this section we study what happens if we select a non-minimal array. Note that we know by Statement 3 of Theorem 1 that no solution should be found. Nevertheless, we like to investigate how the failure to find a solution occurs. Hence, we consider the choice $K=2$. We will show that Problem 2 possesses no solution with these choices. The same holds for any $K>1$. Thus, $K1$ is the only choice for the height of $A$.
Here follow the details. The mode change array with $K=2$ is the following:

$$A_2 = \left\{ \begin{array}{lll} 0 = x'' + \lambda x & (f_1) & t_* \\ 0 = y'' + \lambda y + g & (f_2) & \dots \\ 0 = L^2 - (x^2+y^2) & (k_1) & \dots \\ 0 = xx' + yy' & (k'_1) & \dots \\ 0 = xx'' + x'^2 + y'^2 + yy'' & (k''_1) & \dots \\ 0 = (x'' + \lambda x)\bullet & (f_3) & t_*+\epsilon \\ 0 = (y'' + \lambda y + g)\bullet & (f_4) & \dots \\ 0 = (L^2 - (x^2+y^2))\bullet & (k'_1) & \dots \\ 0 = (xx' + yy')\bullet & (k'_1) & \dots \\ 0 = (xx'' + x'^2 + y'^2 + yy'')\bullet & (k''_1) & \dots \\ 0 = (x'' + \lambda x)\bullet^2 & (f_5) & t_*+2\epsilon \\ 0 = (y'' + \lambda y + g)\bullet^2 & (f_6) & \dots \\ 0 = (L^2 - (x^2+y^2))\bullet^2 & (f_7) & \dots \\ 0 = (xx' + yy')\bullet^2 & (f_8) & \dots \\ 0 = (xx'' + x'^2 + y'^2 + yy'')\bullet^2 & (f_9) & \dots \\ 0 = x'' - \epsilon^{-2}(x^{\bullet 2} - 2x^{\bullet} + x) & (f_{10}) & \text{Euler ids.} \\ 0 = y'' - \epsilon^{-2}(y^{\bullet 2} - 2y^{\bullet} + y) & (f_{11}) & \dots \\ 0 = x''\bullet - \epsilon^{-1}(x'^{\bullet 2} - x'^{\bullet}) & (f_{12}) & \dots \\ 0 = y''\bullet - \epsilon^{-1}(y'^{\bullet 2} - y'^{\bullet}) & (f_{13}) & \dots \\ 0 = x'\bullet - \epsilon^{-1}(x^{\bullet 2} - x^{\bullet}) & (f_{14}) & \dots \\ 0 = y'\bullet - \epsilon^{-1}(y'^{\bullet 2} - y'^{\bullet}) & (f_{15}) & \dots \end{array} \right.$$

$\mathcal{X}^-$ is the same as for $K=1$. Facts and conflicts are pointed. The subsystem in black is structurally nonsingular, with a perfect matching $M$ highlighted in blue. In the right most column we indicate the origin of each equation: for example, $t_*+\epsilon$ indicates that the corresponding equation originates from the 1-shifted discretized dynamics.

The corresponding rescaling calculus is:

$$\begin{align*} \mu_{x''} &= \mu_{f_1} \ge \mu_{\lambda} \\ \mu_{\lambda} &= \mu_{f_2} \ge \mu_{y''} \\ \mu_{x''\bullet} &= \mu_{f_3} \ge \mu_{\lambda\bullet} \\ \mu_{\lambda\bullet} &= \mu_{f_4} \ge \mu_{y''\bullet} \\ \mu_{x''\bullet 2} &= \mu_{f_5} \ge \mu_{\lambda\bullet 2} \\ \mu_{\lambda\bullet 2} &= \mu_{f_6} \ge \mu_{y''\bullet 2} \\ \mu_{f_7} &= \mu_{x\bullet 2} = \mu_{y\bullet 2} = 0 \\ \mu_{y\bullet 2} + \mu_{y'\bullet 2} &= \mu_{f_8} \ge \mu_{x\bullet 2} + \mu_{x'\bullet 2} \\ \mu_{y\bullet 2} + \mu_{y''\bullet 2} &= \mu_{f_9} \ge \max(\mu_{x\bullet 2} + \mu_{x''\bullet 2}, 2\mu_{x'\bullet 2}, 2\mu_{y'\bullet 2}) \\ 2 + \mu_{x\bullet 2} &= \mu_{f_{10}} \ge \mu_{x''} \\ \mu_{y''} &= \mu_{f_{11}} \ge 2 + \mu_{y\bullet 2} \\ 1 + \mu_{x'\bullet} &= \mu_{f_{12}} \ge \mu_{x''\bullet} \\ \mu_{y''\bullet} &= \mu_{f_{13}} \ge 1 + \mu_{y'\bullet} \\ \mu_{x'\bullet} &= \mu_{f_{14}} \ge 1 + \mu_{x\bullet 2} \\ \mu_{y'\bullet} &= \mu_{f_{15}} \ge 1 + \mu_{y\bullet 2} \end{align*}$$

Note that the rescaling equation for $f_7$ corresponds to rescaling equation (58) for nonlinear functions. By goodness condition (52), the following equations and variables possess a zero rescaling offset: $f_7, x, y, x^{\bullet 2}, y^{\bullet 2}, x'^{\bullet 2}, y'^{\bullet 2}, x''^{\bullet 2}, y''^{\bullet 2}, \lambda^{\bullet 2}$. The solution of the rescaling calculus is the following, where we list the variables and equations with rescaling offsets 2, 1, 0:

$$2 : f_1, f_2, f_3, f_4, f_{10}, f_{11}, f_{12}, f_{13}, x'', y'', \lambda, x''', y''', \lambda'$$

$$1 : f_{14}, f_{15}, x', y'$$

0 : other equations and variables

The rescaling offsets violating goodness condition (53) were highlighted in red. Hence, Problem 2 possesses no solution for $K=2$.
Still, suppose that we insist producing the restart system following Procedure 1. The useful part of the rescaled system $A_2^↓$ is the following (rescaled Euler identities have been directly used to expand impulsive derivatives, and are no longer shown):

$$A_2^↓ = \begin{cases} 0 = x^{\bullet 2} - 2x^{\bullet} + x + \lambda^{\downarrow}x \\ 0 = y^{\bullet 2} - 2y^{\bullet} + y + \lambda^{\downarrow}y \\ 0 = x^{\bullet 3} - 2x^{\bullet 2} + x^{\bullet} + \lambda^{\bullet \downarrow}x^{\bullet} \\ 0 = y^{\bullet 3} - 2y^{\bullet 2} + y^{\bullet} + \lambda^{\bullet \downarrow}y^{\bullet} \\ 0 = (L^2 - (x^2 + y^2))^{\bullet 2} \\ 0 = (xx' + yy')^{\bullet 2} \end{cases}$$

The restart system follows by performing renaming (50). The array involves variables beyond the tail instant $t_* + 2\epsilon$, namely $x^{\bullet 3}$, $y^{\bullet 3}$. Hence, the renaming (50) is no longer bijective, which kills the structural nonsingularity of the restart system:

$$\text{restart system} = \begin{cases} 0 = x^{+} - x^{-} + \lambda^{\downarrow}x^{-} \\ 0 = y^{+} - y^{-} + \lambda^{\downarrow}y^{-} \\ 0 = x^{+} + x^{-} + \lambda^{\bullet \downarrow}x^{-} \\ 0 = y^{+} + y^{-} + \lambda^{\bullet \downarrow}y^{-} \\ 0 = L^2 - ((x^{+})^2 + (y^{+})^2) \\ 0 = x^{+}x'^{+} + y^{+}y'^{+} \end{cases}$$

The resulting restart system possesses 6 equations and 6 dependent variables: $x^+$, $x'^+$, $y^+$, $y'^+$, $\lambda^↓$, $\lambda^{\bullet \downarrow}$. Nevertheless it is structurally singular: $x'^+$ and the third equation are both unmatched. However, adding the equation $\lambda^↓ = \lambda^{\bullet \downarrow}$ has the following effects:

* The first and third equations become identical; so we can discard the third equation: the resulting system is no longer structurally conflicting;

* Still, variable $x'^+$ remains unmatched.

The conclusion is that the additional non-structural post-processing we applied did not provide any further progress toward finding a solution to the hot restart.

## D.2 The cup-and-ball with exogenous mode change

In this section we modify the cup-and-ball example, by making the control of mode changes external (with reference to the original model (4), we removed equation ($k_0$)):

$$\left\{ \begin{array}{ll} 0 = x'' + \lambda x & (e_1) \\ 0 = y'' + \lambda y + g & (e_2) \\ \text{if } \gamma \text{ then } 0 = L^2 - (x^2 + y^2) & (k_1) \\ \text{and } 0 = \lambda + s & (k_2) \\ \text{if not } \gamma \text{ then } 0 = \lambda & (k_3) \\ \text{and } 0 = (L^2 - (x^2 + y^2)) - s & (k_4) \end{array} \right.$$
We focus again on the mode change $\gamma : F \to T$ and we show the mode change array with $K = 2$:

$$
\left\{
\begin{array}{l@{\qquad}l@{\qquad}l}
0 = x'' + \lambda x & (e_1) & \gamma=T, t_* \\
0 = y'' + \lambda y + g & (e_2) & \dots \\
0 = L^2 - (x^2+y^2) & (k_1) & \dots \\
0 = xx' + yy' & (k'_1) & \dots \\
0 = xx'' + x'^2 + y'^2 + yy'' & (k''_1) & \dots \\
0 = (x'' + \lambda x)\bullet & (e'_1) & \gamma=T, t'_* \\
0 = (y'' + \lambda y)\bullet + g & (e'_2) & \dots \\
0 = (L^2 - (x^2+y^2))\bullet & (k'_1) & \dots \\
0 = (xx' + yy')\bullet & (k'_1) & \dots \\
0 = (xx'' + x'^2 + y'^2 + yy'')\bullet & (k''_1) & \dots \\
0 = (L^2 - (x^2+y^2))\bullet^2 & (k_1^2) & \gamma=T, t_*^2 \\
0 = (xx' + yy')\bullet^2 & (k'_1^2) & \dots
\end{array}
\right.
\tag{105}
$$

Observe the effect of the change in the model. No fact occurs. Instead, we have two more conflicting equations, namely $(k_1)$ and $(k'_1)$. This causes the need for a longer array compared to the one in Fig. 2 for the original cup-and-ball. In (105), the conflicts were pointed. Corresponding equations are disabled. A perfect matching modulo $\sim$ is highlighted in blue. We make it explicit by adding the needed Euler identities:

$$
\left\{
\begin{array}{l@{\qquad}l}
0 = x'' + \lambda x & (e_1 : 1) \\
0 = y'' + \lambda y + g & (e_2 : 2) \\
0 = (x'' + \lambda x)\bullet & (e_1^\bullet : 3) \\
0 = (y'' + \lambda y)\bullet + g & (e_2^\bullet : 4) \\
0 = (L^2 - (x^2+y^2))\bullet^2 & (k_1^\bullet^2 : 5) \\
0 = (xx' + yy')\bullet^2 & (k'_1^\bullet^2 : 6) \\
0 = \varepsilon^2 \times x'' - (x^\bullet^2 - 2x^\bullet + x) & (\mathcal{E}_{x'',x^\bullet^2} : 7) \\
0 = \varepsilon^2 \times y'' - (y^\bullet^2 - 2y^\bullet + y) & (\mathcal{E}_{y'',y^\bullet^2} : 8) \\
0 = \varepsilon \times x'^\bullet - (x^\bullet^2 - x^\bullet) & (\mathcal{E}_{x'^\bullet,x^\bullet^2} : 9) \\
0 = \varepsilon \times y'^\bullet - (y^\bullet^2 - y^\bullet) & (\mathcal{E}_{y'^\bullet,y^\bullet^2} : 10) \\
0 = \varepsilon \times x''^\bullet - (x'^\bullet^2 - x'^\bullet) & (\mathcal{E}_{x''^\bullet,x'^\bullet^2} : 11) \\
0 = \varepsilon \times y''^\bullet - (y'^\bullet^2 - y'^\bullet) & (\mathcal{E}_{y''^\bullet,y'^\bullet^2} : 12)
\end{array}
\right.
\tag{106}
$$

The resulting rescaling calculus is:

$$
\begin{align*}
\mu_{x''} &= \mu_{f_1} \ge \mu_{\lambda} \\
\mu_{\lambda} &= \mu_{f_2} \ge \mu_{y''} \\
\mu_{x''\bullet} &= \mu_{f_3} \ge \mu_{\lambda\bullet} \\
\mu_{\lambda\bullet} &= \mu_{f_4} \ge \mu_{y''\bullet} \\
\mu_{f_5} &= \mu_{y^\bullet 2} = \mu_{x^\bullet 2} = 0 \\
\mu_{y^\bullet 2} + \mu_{y^\bullet 2} &= \mu_{f_6} \ge \mu_{x'^\bullet 2} + \mu_{x^\bullet 2} \\
2 + \mu_{x^\bullet 2} &= \mu_{f_7} \ge \mu_{x''} \\
2 + \mu_{y^\bullet 2} &= \mu_{f_8} \ge \mu_{y''} \\
1 + \mu_{x'^\bullet} &= \mu_{f_9} \ge \mu_{x''\bullet} \\
\mu_{y''\bullet} &= \mu_{f_{10}} \ge 1 + \mu_{y'^\bullet} \\
\mu_{x'^\bullet} &= \mu_{f_{11}} \ge 1 + \mu_{x^\bullet 2} \\
\mu_{y^\bullet \cdot} &= \mu_{f_{12}} \ge 1 + \mu_{y^\bullet 2}
\end{align*}
$$

The rescaling equation for $f_5$ corresponds to rescaling equation (58) for nonlinear functions. The solution of the rescaling calculus is the following, where we list the variables and equations with rescaling offsets 2, 1, 0:

$2 : f_1, f_2, f_3, f_4, f_7, f_8, f_9, f_{10}, x'', y'', \lambda, x''\bullet, y''\bullet, \lambda\bullet$

$1 : f_{11}, f_{12}, x^{\prime\bullet}, y^{\prime\bullet}$

$0 : \text{other equations and variables}$

The rescaling offsets violating goodness condition (53) are highlighted in red. This solution violates goodness condition (53), hence Problem 2 possesses no good solution. The considered mode change is insufficiently determined.
Since goodness conditions (51,52) are satisfied, we can be more precise, however, by applying Corollary 1. The rescaling offsets violating goodness condition (53) are highlighted in red. The involved variables are $x''\bullet$ and $y''\bullet$. Remove the equations that are matched with these variables, namely $e_1^\bullet$ and $\mathcal{E}_{y''\bullet, y\bullet2}$. Then, a Dulmage-Mendelsohn decomposition of the system of (106), we denote it by $A_\epsilon$ yields the following result:

$$A_{\epsilon}^{o} = \emptyset$$

$$A_{\epsilon}^{r} = ((e_1, e_2, k_1^{\bullet 2}, \mathcal{E}_{x'',x\bullet2}, \mathcal{E}_{y'',y\bullet2}), (x'', x^{\bullet 2}, y'', y^{\bullet 2}, \lambda))$$

$$A_{\epsilon}^{u} = \text{other equations and variables}$$

Rescaling of the regular subsystem $A_\epsilon^r$ yields

$$\left\{ \begin{array}{l} 0 = x^{\bullet 2} - 2x^{\bullet} + (1 + \lambda^{\downarrow})x \\ 0 = y^{\bullet 2} - 2y^{\bullet} + (1 + \lambda^{\downarrow})y \\ 0 = (y^{\bullet 2})^2 + (y^{\bullet 2})^2 \end{array} \right.$$

which maps to the restart system

$$\text{restart system} = \left\{ \begin{array}{l} 0 = x^{+} - x^{-} + \lambda^{\downarrow}x^{-} \\ 0 = y^{+} - y^{-} + \lambda^{\downarrow}y^{-} \\ 0 = (y^{+})^{2} + (y^{+})^{2} \end{array} \right. \qquad (107)$$

Conclusion: the positions are determined by (107), whereas the velocities remain undetermined.

### D.3 A “strange” cup-and-ball

With this example we aim at analyzing difficult cases related to facts. The following “strange” cup-and-ball is a different modification of the cup-and-ball example, by which the zero-crossing function defining the mode change was modified, which modifies the facts. Its model is the following, where the modification is highlighted in red:

$$\left\{ \begin{array}{ll} 0 = x'' + \lambda x & (e_1) \\ 0 = y'' + \lambda y + g & (e_2) \\ \gamma = [s'^- \le 0]; \gamma(0) = \text{F} & (k_0) \\ \text{if } \gamma \text{ then } 0 = L^2 - (x^2 + y^2) & (k_1) \\ \text{and } 0 = \lambda + s & (k_2) \\ \text{if not } \gamma \text{ then } 0 = \lambda & (k_3) \\ \text{and } 0 = (L^2 - (x^2 + y^2)) - s & (k_4) \end{array} \right.$$

The dynamics in modes $\gamma = \text{F}$ and $\gamma = \text{T}$ are not modified. However, in mode $\gamma = \text{F}$, we have

$$s' = (L^2 - (x^2 + y^2))' = 2(xx' + yy').$$

Consequently, the candidate fact is now $xx'+yy' = 0$ (compared to $L^2-(x^2+y^2) = 0$ for the cup-and-ball). We show next the mode change array $A_2$:

$$\left\{ \begin{array}{lll} 0 = x'' + \lambda x & (e_1) & \gamma=\text{T}, t_* \\ 0 = y'' + \lambda y + g & (e_2) & \dots \\ 0 = L^2 - (x^2 + y^2) & (k_1) & \dots \\ 0 = xx' + yy' & (k'_1) & \dots \\ 0 = xx'' + x'^2 + y'^2 + yy'' & (k''_1) & \dots \\ 0 = (x'' + \lambda x)\bullet & (e_1^\bullet) & \gamma=\text{T}, t_*^\bullet \\ 0 = (y'' + \lambda y)\bullet + g & (e_2^\bullet) & \dots \\ 0 = (L^2 - (x^2 + y^2))\bullet & (k'_1^\bullet) & \dots \\ 0 = (xx' + yy')\bullet & (k'_1^\bullet) & \dots \\ 0 = (xx'' + x'^2 + y'^2 + yy'')\bullet & (k''_1^\bullet) & \dots \\ 0 = (L^2 - (x^2 + y^2))\bullet^2 & (k'_1^\bullet^2) & \gamma=\text{T}, t_*^\bullet^2 \\ 0 = (xx' + yy')\bullet^2 & (k'_1^\bullet^2) & \dots \end{array} \right. \qquad (108)$$
The fact is shown in green and the conflicts are shown in red. The dependent variables are

$$ \lambda, x^{\bullet 2} \sim x'^{\bullet} \sim x'', y^{\bullet 2} \sim y'^{\bullet} \sim y'', \lambda^{\bullet}, x'^{\bullet 2} \sim x''^{\bullet}, y'^{\bullet 2} \sim y''^{\bullet} $$

The black subsystems of (108) and (105) (cup-and-ball with exogenous mode change) coincide. Consequently, the analyses of the two examples coincide as well.

## E Nonlinear systems: a clutch

<figure><img src="image_14.png" /><figcaption>Figure 6: A simple clutch with two shafts and its model from first principles.</figcaption></figure>

The clutch and its model from first principles are shown on Fig. 6. We separately model each shaft: rotation velocities are $\omega_i, i=1,2$ and torques are $\tau_i, i=1,2$. In the model of their junction, the Boolean $\gamma$ is an input. The junction exhibits two modes: *engaged* $\gamma=T$, and *released* $\gamma=F$. We focus on the mode change from released to engaged $\gamma:F \to T$. The physics suggests the following: 1) The mode change $\gamma:F \to T$ should be determined; hot restart differs from cold start, which possesses one degree of freedom (the common rotation velocity $\omega_1 = \omega_2$ is not determined); 2) The torques will experience an impulsion. Physically meaningful modeling leads to take $f_1, f_2$ linear with respect to torques. We will explore the consequences of choosing $f_1, f_2$ nonlinear. Here follows the mode change array $A_1$ for this example:

$$ A_1: \left\{ \begin{array}{ll} \begin{array}{l} \omega'_1 = f_1(\omega_1, \tau_1) \\ \omega'_2 = f_2(\omega_2, \tau_2) \\ \omega_1 = \omega_2 \\ \omega'_1 = \omega'_2 \\ \tau_1 + \tau_2 = 0 \end{array} & (e_1) \quad \gamma=T, t \\ \\ \begin{array}{l} \omega'_1 = \varepsilon^{-1}(\omega_1^\bullet - \omega_1) \\ \omega'_2 = \varepsilon^{-1}(\omega_2^\bullet - \omega_2) \end{array} & (\mathcal{E}_1) \\ \\ \begin{array}{l} \omega'_1^\bullet = f_1(\omega_1^\bullet, \tau_1^\bullet) \\ \omega'_2^\bullet = f_2(\omega_2^\bullet, \tau_2^\bullet) \\ \omega_1^\bullet = \omega_2^\bullet \\ \omega'_1^\bullet = \omega'_2^\bullet \\ \tau_1^\bullet + \tau_2^\bullet = 0 \end{array} & (e_2) \quad \gamma=T, t+\varepsilon \\ \\ \begin{array}{l} \omega'_1^\bullet = f_1(\omega_1^\bullet, \tau_1^\bullet) \\ \omega'_2^\bullet = f_2(\omega_2^\bullet, \tau_2^\bullet) \\ \omega_1^\bullet = \omega_2^\bullet \\ \omega'_1^\bullet = \omega'_2^\bullet \\ \tau_1^\bullet + \tau_2^\bullet = 0 \end{array} & (e_3) \quad \gamma=T, t+\varepsilon \\ \\ \begin{array}{l} \omega'_1^\bullet = f_1(\omega_1^\bullet, \tau_1^\bullet) \\ \omega'_2^\bullet = f_2(\omega_2^\bullet, \tau_2^\bullet) \\ \omega_1^\bullet = \omega_2^\bullet \\ \omega'_1^\bullet = \omega'_2^\bullet \\ \tau_1^\bullet + \tau_2^\bullet = 0 \end{array} & (e_4) \quad \gamma=T, t+\varepsilon \end{array} \right. \quad (109) $$

Mode change array $A_1$ showing the dynamics around the mode change $\gamma:F \to T$ occurring at instant $t$, by stacking the discrete time dynamics at instants $t$, and $t+\varepsilon$. Disabled equations are highlighted
in red. The enabled subsystem (in black) is structurally nonsingular, as evidenced by the one-to-one matching between variables and equations highlighted in blue. We have $\mathcal{X}^- = \{\omega_1, \omega_2\}$. Since we focus on state variables for restart, we are not interested in leading equations of instant $t_*+1 = t_*^\bullet$. By omitting equations of this instant, the enabled part of the mode change array $A_1$ in (109) is:

$$ 
\left\{ 
\begin{array}{l} 
J_1\omega_1' = f_1(\tau_1, \omega_1) \quad (f_1) \\ 
J_2\omega_2' = f_2(\tau_2, \omega_2) \quad (f_2) \\ 
\tau_1 + \tau_2 = 0 \quad (f_3) \\ 
\omega_1^\bullet = \omega_2^\bullet \quad (f_4) \\ 
\omega_1' = \varepsilon^{-1} \times (\omega_1^\bullet - \omega_1) \quad (f_5) \\ 
\omega_2' = \varepsilon^{-1} \times (\omega_2^\bullet - \omega_2) \quad (f_6) 
\end{array} 
\right. 
$$

So far, no particular form was assumed for $f_i$, $i=1, 2$.

### Rescaling calculus with $f_i$ linear in $\tau_i$

Assume $f_i(\tau_i, \omega_i) = \tau_i - g_i(\omega_i)$ is linear with respect to the torque. The resulting rescaling calculus is:

$$ 
\left\{ 
\begin{array}{l} 
\mu_{\omega_1'} = \max(\mu_{\omega_1'}, \mu_{\tau_1}) \\ 
\mu_{\tau_2} = \max(\mu_{\omega_2'}, \mu_{\tau_2}) \\ 
\mu_{\tau_1} = \max(\mu_{\tau_1}, \mu_{\tau_2}) \\ 
\mu_{\omega_2^\bullet} = \max(\mu_{\omega_2^\bullet}, \mu_{\omega_1^\bullet}) \\ 
\mu_{\omega_1^\bullet} + 1 = \max(\mu_{\omega_1^\bullet} + 1, \mu_{\omega_1'}) \\ 
\mu_{\omega_2^\bullet} = \max(\mu_{\omega_2^\bullet}, \mu_{\omega_2^\bullet} + 1) 
\end{array} 
\right. 
\qquad (110) 
$$

The following solution of this rescaling analysis yields a good solution for Problem 2:

$$ 
\begin{array}{l} 
1 : \omega_1', \omega_2', \tau_1, \tau_2; f_1, f_2, f_3, f_5, f_6 \\ 
0 : \omega_1^\bullet, \omega_2^\bullet; f_4 
\end{array} 
\qquad (111) 
$$

Expanding $\omega_i^{\prime\downarrow} = \omega_i^\bullet - \omega_i$, $i = 1, 2$. The rescaled system is

$$ 
\left\{ 
\begin{array}{l} 
J_1(\omega_1^\bullet - \omega_1) = \tau_1^{\downarrow} \\ 
J_2(\omega_2^\bullet - \omega_2) = \tau_2^{\downarrow} \\ 
\tau_1^{\downarrow} + \tau_2^{\downarrow} = 0 \\ 
\omega_1^\bullet = \omega_2^\bullet 
\end{array} 
\right. 
$$

which yields, by setting $\omega^+ = \text{def} \omega_1^\bullet = \omega_2^\bullet$, the restart system:

$$ (J_1 + J_2)\omega^+ = J_1\omega_1^- + J_2\omega_2^-, $$

reflecting the preservation of angular momentum. This again illustrates Theorem 4. Note that having $g_i$ nonlinear was not an obstacle, since the rotation velocities are not impulsive at the mode change (they are past variables).

### Rescaling calculus with $f_i$ nonlinear in $\tau_i$

Assume now that $f_i$ is nonlinear in its arguments. Rescaling calculus (110) modifies as follows:

$$ 
\left\{ 
\begin{array}{l} 
\mu_{\omega_1'} = \mu_{f_1} = \text{if } \max(\mu_{\tau_1}, \mu_{\omega_1}) = 0 \quad \text{then } 0 \quad \text{else } \infty \\ 
\mu_{\tau_2} = \mu_{f_2} = \text{if } \max(\mu_{\tau_2}, \mu_{\omega_2}) = 0 \quad \text{then } 0 \quad \text{else } \infty \\ 
\mu_{\tau_1} = \max(\mu_{\tau_1}, \mu_{\tau_2}) \\ 
\mu_{\omega_2^\bullet} = \max(\mu_{\omega_2^\bullet}, \mu_{\omega_1^\bullet}) \\ 
\mu_{\omega_1^\bullet} + 1 = \max(\mu_{\omega_1^\bullet} + 1, \mu_{\omega_1'}) \\ 
\mu_{\omega_2^\bullet} = \max(\mu_{\omega_2^\bullet}, \mu_{\omega_2^\bullet} + 1) 
\end{array} 
\right. 
\qquad (112) 
$$
The last equation of this rescaling calculus prevents the existence of a good solution to Problem 2 since it prohibits $\mu_{\omega_2'} = 0$. Our approach does not solve the restart problem in this case.

### What if the junction is nonlinear in velocities

We keep $f_i = \tau_i - g(\omega_i)$, $i=1, 2$ but we change the junction model by changing $(e_3)$ to $\omega_2 = h(\omega_1)$, where $h$ is nonlinear. The rescaling calculus is now

$$ 
\left\{ 
\begin{array}{l}
\mu_{\omega_1'} = \max(\mu_{\omega_1'}, \mu_{\tau_1}) \\
\mu_{\tau_2} = \max(\mu_{\omega_2'}, \mu_{\tau_2}) \\
\mu_{\tau_1} = \max(\mu_{\tau_1}, \mu_{\tau_2}) \\
\\
\mu_{\omega_2'} = \max(\mu_{\omega_2'}, \text{if } \mu_{\omega_1'} = 0 \text{ then } 0 \text{ else } \infty) \\
\\
\mu_{\omega_1'} + 1 = \max(\mu_{\omega_1'} + 1, \mu_{\omega_1'}) \\
\mu_{\omega_2'} = \max(\mu_{\omega_2'}, \mu_{\omega_2'} + 1)
\end{array}
\right. 
\qquad (113) 
$$

for which (111) is still a good solution. Making the junction model non-polynomial added no difficulty with reference to the original clutch.

## F Comparing with Trenn et al.

### F.1 A linear DAE system

Consider the following linear model

$$ 
\left\{ 
\begin{array}{lll}
0 = v_1' - i & (f_1) \\
0 = v_2' - i + v_R & (f_2) \\
0 = i' - v_R & (f_3) \\
0 = v_1 + u & (f_4)
\end{array}
\right. 
$$

The dependent variables of this DAE are $i, v_1, v_2, v_R$ and $u$ is an input with $u, u', u''$ known at $t = 0^-$. There is a unique perfect matching (shown in blue) and its total weight is 1. The equations defining the offsets are (with $c_i =_{\text{def}} c_{f_i}$)

$$ 
\begin{array}{lclcl}
d_i & = & c_1 & c_4 - c_1 & \geq & 1 \\
d_{v_2} & = & 1 + c_2 & c_1 - c_2 & \geq & 0 \\
d_{v_R} & = & c_3 & c_3 - c_2 & \geq & 0 \\
d_{v_1} & = & c_4 & c_1 - c_3 & \geq & 1
\end{array} 
$$

The minimal solution is

$$ 
\begin{array}{l}
c_1 = 1, c_2 = c_3 = 0, c_4 = 2 \\
d_i = 1, d_{v_1} = 2, d_{v_2} = 1, d_{v_R} = 0
\end{array} 
$$

The index reduced system is

$$ 
\left\{ 
\begin{array}{l}
0 = v_1' - i \qquad (f_1) \\
0 = v_1'' - i' \qquad (f_1') \\
0 = v_2' - i + v_R \qquad (f_2) \\
0 = i' - v_R \qquad (f_3) \\
0 = v_1 + u \qquad (f_4) \\
0 = v_1' + u' \qquad (f_4') \\
0 = v_1'' + u'' \qquad (f_4'')
\end{array}
\right. 
$$

and we show in blue a perfect matching for the leading variables $i', v_2', v_R, v_1''$. The state variables are $i, v_2, v_1, v_1'$, $u, u'$, subject to the consistency constraints $(f_1, f_4)$, and $u''$ is a free (input) variable.
The mode change array $A_1$ completed with Euler identities, in which we omit the leading equations of the last block, together with the resulting rescaling calculus, is the following:

$$
\left\{
\begin{array}{lll}
0 = v_1' - i & \mu_i \ge \mu_{v_1'} & (f_1) \\
0 = v_1'' - i' & \text{conflict} & \\
0 = i' - v_R & \mu_{v_R} \ge \mu_{i'} & (f_3) \\
0 = v_1 + u & \mu_{v_1} \ge 0 & (f_4) \\
0 = v_1' + u' & \text{conflict} & \\
0 = v_1'' + u'' & \text{conflict} & \\
\\
0 = v_1^{\bullet} - i^{\bullet} & \mu_{i^{\bullet}} \ge \mu_{v_1^{\bullet}} & (f_5) \\
0 = v_1^{\bullet} + u^{\bullet} & \mu_{v_1^{\bullet}} \ge 0 & (f_6) \\
0 = v_1^{\bullet} + u^{\bullet} & \mu_{v_1^{\bullet}} \ge \mu_{u^{\bullet}} & (f_7) \\
\\
0 = \varepsilon \times i' - (i^{\bullet} - i) & \mu_{i^{\bullet}} \ge 1 + \max(\mu_{i^{\bullet}}, \mu_i) & (f_8) \\
0 = \varepsilon \times v_1' - (v_1^{\bullet} - v_1) & \mu_{v_1^{\bullet}} \ge 1 + \max(\mu_{v_1^{\bullet}}, \mu_{v_1}) & (f_9) \\
0 = \varepsilon \times v_1'' - (v_1^{\bullet} - v_1') & \mu_{v_1''} \ge 1 + \max(\mu_{v_1^{\bullet}}, \mu_{v_1'}) & (f_{10}) \\
0 = \varepsilon \times u'' - (u^{\bullet} - u') & \mu_{u''} \ge 1 + \max(\mu_{u^{\bullet}}, \mu_{u'}) & (f_{12})
\end{array}
\right.
$$

The values of •u, •u', •u'' are given. Hence, the set of past variables is

$$
\mathcal{X}^{-} = \{ u = \bullet u + \varepsilon \bullet u', u' = \bullet u' + \varepsilon \bullet u'', u^{\bullet} = \bullet u + 2\varepsilon \bullet u' + \varepsilon^2 \bullet u'' \}
$$

The solution of the rescaling calculus is the following, where we list the variables and equations
with rescaling offsets 2, 1, 0:

2 : $v_R, i', v_1''$; $f_{10}, f_8, f_3$

1 : $u'', i, v_1'$; $f_9, f_1, f_{12}$

0 : other variables and equations

Goodness condition (53) is violated by variable $v_1''$, highlighted in red. So this mode change is
nondetermined according to our approach. Equation $f_{10}$ is matched with $v_1''$. Hence, we remove it,
and perform the Dulmage-Mendelsohn decomposition of the so reduced mode change array $A_1$. The
following variables are reachable from $v_1''$ by following an alternating path: $v_1, v_1', i, v_1^\bullet, i', \dot{i}^\bullet, v_1^\bullet, u^\bullet$;
no variable is determined.

Our method fails finding a solution to hot restart, unlike Trenn's method. The reason is that we do not allow for changes in the state basis, which is precisely performed by the Weierstrass decomposition [38, 39, 24].

F.2 Further investigating the clutch

Being structural, our approach does not recombine states. We illustrate some consequences of this
on a linear version of the clutch, namely: $f_i(\omega_i, \tau_i) = \tau_i - a_i\omega_i$, where $a_i$ are constants, already
studied in Appendix E.

Now, reconsider the same example, but with the change of variable $(\omega_i, \tau_i) \mapsto (w_i, \sigma_i)$, where:

$$
w_1 = \omega_1 - \tau_1, w_2 = \omega_2 - \tau_2, \sigma_1 = \omega_1 + \tau_1, \sigma_2 = \omega_2 + \tau_2.
$$

or, equivalently,

$$
2\omega_1 = w_1 + \sigma_1, 2\omega_2 = w_2 + \sigma_2, 2\tau_1 = \sigma_1 - w_1, 2\tau_2 = \sigma_2 - w_2.
$$

The dynamics in the F mode rewrites

$$
\left\{
\begin{array}{l}
0 = w'_1 + \sigma'_1 - (\sigma_1 - w_1 - a(\sigma_1 + w_1)) \\
0 = w'_2 + \sigma'_2 - (\sigma_2 - w_2 - a(\sigma_2 + w_2)) \\
0 = w_1 - \sigma_1 \\
0 = w_2 - \sigma_2
\end{array}
\right.
$$
It is now a DAE, structurally! The index reduction yields

$$
\begin{cases}
0 = w_1' + \sigma_1' - (\sigma_1 - w_1 - a(\sigma_1 + w_1)) \\
0 = w_2' + \sigma_2' - (\sigma_2 - w_2 - a(\sigma_2 + w_2)) \\
0 = w_1 - \sigma_1 \\
0 = w_2 - \sigma_2 \\
0 = w_1' - \sigma_1' \\
0 = w_2' - \sigma_2'
\end{cases}
$$

Hence, for the mode change, $\mathcal{X}^{-} = \{w_1, w_2, \sigma_1, \sigma_2\}$. The mode change array for $\gamma: F \to T$ rewrites

$$
\left\{
\begin{aligned}
0 &= w_1' + \sigma_1' - (\sigma_1 - w_1 - a(\sigma_1 + w_1)) \\
0 &= w_2' + \sigma_2' - (\sigma_2 - w_2 - a(\sigma_2 + w_2)) \\
0 &= (w_1 + \sigma_1) - (w_2 + \sigma_2) \\
0 &= (\sigma_1 + \sigma_2) - (w_1 + w_2) \\
0 &= (w_1' + \sigma_1') - (w_2' + \sigma_2') \\
0 &= (\sigma_1' + \sigma_2') - (w_1' + w_2') \\
0 &= w_1'^{\bullet} + \sigma_1'^{\bullet} - (\sigma_1^{\bullet} - w_1^{\bullet} - a(\sigma_1^{\bullet} + w_1^{\bullet})) \\
0 &= w_2'^{\bullet} + \sigma_2'^{\bullet} - (\sigma_2^{\bullet} - w_2^{\bullet} - a(\sigma_2^{\bullet} + w_2^{\bullet})) \\
0 &= (w_1^{\bullet} + \sigma_1^{\bullet}) - (w_2^{\bullet} + \sigma_2^{\bullet}) \\
0 &= (\sigma_1^{\bullet} + \sigma_2^{\bullet}) - (w_1^{\bullet} + w_2^{\bullet}) \\
0 &= (w_1'^{\bullet} + \sigma_1'^{\bullet}) - (w_2'^{\bullet} + \sigma_2'^{\bullet}) \\
0 &= (\sigma_1'^{\bullet} + \sigma_2'^{\bullet}) - (w_1'^{\bullet} + w_2'^{\bullet})
\end{aligned}
\right.
$$

Keeping only the active equations, removing the leading equations of the last block, and adding the
Euler identities yields

$$
\left\{
\begin{aligned}
0 &= w_1' + \sigma_1' - (\sigma_1 - w_1 - a(\sigma_1 + w_1)) \\
0 &= w_2' + \sigma_2' - (\sigma_2 - w_2 - a(\sigma_2 + w_2)) \\
0 &= (w_1^{\bullet} + \sigma_1^{\bullet}) - (w_2^{\bullet} + \sigma_2^{\bullet}) \\
0 &= (\sigma_1^{\bullet} + \sigma_2^{\bullet}) - (w_1^{\bullet} + w_2^{\bullet}) \\
0 &= \varepsilon \times w_1' - (w_1^{\bullet} - w_1) \\
0 &= \varepsilon \times w_2' - (w_2^{\bullet} - w_2) \\
0 &= \varepsilon \times \sigma_1' - (\sigma_1^{\bullet} - \sigma_1) \\
0 &= \varepsilon \times \sigma_2' - (\sigma_2^{\bullet} - \sigma_2)
\end{aligned}
\right.
$$

The rescaling calculus is

$$
\left\{
\begin{array}{l@{\qquad}l@{\qquad}l}
\mu_{w_1'} = \mu_{f_1} & \geq & \mu_{\sigma_1'} \\
\mu_{w_2'} = \mu_{f_2} & \geq & \mu_{\sigma_2'} \\
\mu_{\sigma_1^\bullet} = \mu_{f_3} & \geq & \max(\mu_{w_1^\bullet}, \mu_{w_2^\bullet}, \mu_{\sigma_2^\bullet}) \\
\mu_{\sigma_2^\bullet} = \mu_{f_4} & \geq & \max(\mu_{w_1^\bullet}, \mu_{w_2^\bullet}, \mu_{\sigma_1^\bullet}) \\
1 + \mu_{w_1^\bullet} = \mu_{f_5} & \geq & \mu_{w_1'} \\
1 + \mu_{w_2^\bullet} = \mu_{f_6} & \geq & \mu_{w_2'} \\
\mu_{\sigma_1'} = \mu_{f_7} & \geq & 1 + \mu_{\sigma_1^\bullet} \\
\mu_{\sigma_2'} = \mu_{f_8} & \geq & 1 + \mu_{\sigma_2^\bullet}
\end{array}
\right.
$$

Its solution is the following, where we list the variables and equations with rescaling offsets 2, 1, 0:

$1 : \sigma_1', \sigma_2', w_1', w_2'; f_1, f_2, f_7, f_8$

$0 : \text{other variables and equations}$

This solution satisfies the Conditions (51–53). Hence, we can rescale the array, thus obtaining:

$$
\left\{
\begin{aligned}
0 &= (w_1^{\bullet} - w_1) + (\sigma_1^{\bullet} - \sigma_1) \\
0 &= (w_2^{\bullet} - w_2) + (\sigma_2^{\bullet} - \sigma_2) \\
0 &= (w_1^{\bullet} + \sigma_1^{\bullet}) - (w_2^{\bullet} + \sigma_2^{\bullet}) \\
0 &= (\sigma_1^{\bullet} + \sigma_2^{\bullet}) - (w_1^{\bullet} + w_2^{\bullet})
\end{aligned}
\right.
$$
At this point we have a problem: the above system is structurally nonsingular; nevertheless, it is numerically singular. The first equation yields $w_1^\bullet = w_1 + \sigma_1 - \sigma_1^\bullet$ and replacing $w_1^\bullet$ by $w_1 + \sigma_1 - \sigma_1^\bullet$ in the third equation cancels $\sigma_1^\bullet$, which cannot be a pivot.

Maybe the change of variables was responsible for this; we modify its coefficients:

$$ \omega_1 = w_1 + 2\sigma_1, \omega_2 = w_2 + \sigma_2, \tau_1 = \sigma_1 - 3w_1, \tau_2 = 2\sigma_2 - w_2. $$

The structural analysis is not modified. So we end up with the rescaled array

$$ \left\{ \begin{array}{l} 0 = (w_1^\bullet - w_1) + 2(\sigma_1^\bullet - \sigma_1) \\ 0 = (w_2^\bullet - w_2) + (\sigma_2^\bullet - \sigma_2) \\ 0 = (w_1^\bullet + 2\sigma_1^\bullet) - (w_2^\bullet + \sigma_2^\bullet) \\ 0 = (\sigma_1^\bullet + 2\sigma_2^\bullet) - (3w_1^\bullet + w_2^\bullet) \end{array} \right. $$

which suffers from the same problem.

The reason for this pitfall is that our change of state basis resulted in combining non-impulsive variables with impulsive ones. This would not be a trouble to the method developed by Trenn et al. [38, 39, 24] for linear systems, since the Weierstrass decomposition that is applied first, recovers the basis in which clean separation is made between non-impulsive and impulsive variables. We do not do this.
Inria

RESEARCH CENTRE
Centre Inria de l'Université de Rennes

Campus universitaire de Beaulieu
Avenue du Général Leclerc
35042 Rennes Cedex

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

ISSN 0249-0803