# The Time-Freezing Reformulation for Numerical Optimal Control of Complementarity Lagrangian Systems with State Jumps - Extended Version

Armin Nurkanović, Sebastian Albrecht, Bernard Brogliato, Moritz Diehl

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

## To cite this version:

Armin Nurkanović, Sebastian Albrecht, Bernard Brogliato, Moritz Diehl. The Time-Freezing Reformulation for Numerical Optimal Control of Complementarity Lagrangian Systems with State Jumps - Extended Version. Automatica, 2023, 158, pp.111295. 10.1016/j.automatica.2023.111295. hal-03427800v2

HAL Id: hal-03427800

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

Submitted on 1 Aug 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The Time-Freezing Reformulation for Numerical Optimal Control of Complementarity Lagrangian Systems with State Jumps - Extended Version *

Armin Nurkanović<sup>a</sup>, Sebastian Albrecht<sup>c</sup>, Bernard Brogliato<sup>d</sup>, Moritz Diehl<sup>a,b</sup>,

## Abstract

This paper introduces a novel time-freezing reformulation and numerical methods for optimal control of complementarity Lagrangian systems (CLS) with state jumps. We cover the difficult case when the system evolves on the boundary of the dynamic's feasible set after the state jump. In nonsmooth mechanics, this corresponds to inelastic impacts. The main idea of the time-freezing reformulation is to introduce a clock state and an auxiliary dynamical system whose trajectory endpoints satisfy the state jump law. When the auxiliary system is active, the clock state is not evolving, hence by taking only the parts of the trajectory when the clock state was active, we can recover the original solution. The resulting time-freezing system is a Filippov system that has jump discontinuities only in the first time derivative instead of the trajectory itself. This enables one to use the recently proposed Finite Elements with Switch Detection [30], which makes high accuracy numerical optimal control of CLS with impacts and friction possible. We detail how to recover the solution of the original system and show how to select appropriate auxiliary dynamics. The theoretical findings are illustrated on a nontrivial numerical optimal control example of a hopping one-legged robot.

**Key words:** nonsmooth and discontinuous problems, modeling for control optimization, numerical algorithms, algorithms and software

## 1 Introduction

Complementarity Lagrangian Systems (CLS) model the dynamics of rigid bodies with friction and impact. They are indispensable in modern robotic control applications, as any complex task requires exploiting contacts and friction [7,18,19,32,35,38]. Although many mature and efficient simulation methods [2,37] and software [19] ex-

ist, solving optimal control problems remains very difficult and is often heuristic-based. In this paper, we introduce a time-freezing reformulation that transforms the CLS with state jumps into a Piecewise Smooth System (PSS), which is treated as a Filippov system [14]. The resulting system has discontinuities only in its vector field, but not in the solution itself anymore. Our main motivation is to formulate optimal control problems that we can solve with high accuracy and to avoid convergence to spurious solutions, as we discuss below in more detail.

From the theoretical side, it is desirable to know whether solutions to optimal control problems involving CLS or Filippov systems exist. Currently, to the best of the authors' knowledge, there are no general results for these problems, but only for some special cases. This is an intriguing and challenging problem but beyond the scope of this paper. We will assume that optimal solutions do exist. Under this assumption, Pontryagin-type condi-

<footnote><sup>a</sup>Department of Microsystems Engineering (IMTEK), University of Freiburg, Germany</footnote>

<footnote><sup>b</sup>Department of Mathematics, University of Freiburg, Germany</footnote>

<footnote><sup>c</sup>Siemens Technology, Munich, Germany</footnote>

<footnote><sup>d</sup>Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, Grenoble, France</footnote>

<footnote>* This research was supported by the German Federal Ministry of Education and Research (BMBF) via the funded Kopernikus project: SynErgie (03SFK3U0), DFG via Research Unit FOR 2401 and project 424107692 and by the EU via ELO-X 953348. Corresponding author A. Nurkanović.

*Email address: armin.nurkanovic@imtek.uni-freiburg.de (Armin Nurkanović), sebastian.albrecht@siemens.com (Sebastian Albrecht), bernard.brogliato@inria.fr (Bernard Brogliato), moritz.diehl@imtek.uni-freiburg.de (Moritz Diehl)</footnote>
tions for dynamic complementarity systems (with abso-
lutely continuous solutions) are provided in [16,39]. Nec-
essary optimality conditions for a controlled first-order
sweeping process [8] are derived in [12]. Existence results
for differential inclusions with Lipschitz properties are
derived in [11,24] and with upper semi-continuity prop-
erties in [10].

From the computational side, time-stepping methods for CLS and differential inclusions have in general first-order accuracy [2,35]. Despite this, they are often used in direct optimal control approaches [27,32,38,39]. However, it was first noticed by Stewart and Anitescu [36] that for ODE with a discontinuous right-hand side (r.h.s.), contrary to smooth ODE, the numerical sensitivities are wrong independent of the step size. This often impairs the progress of the optimizer and results in spurious local minima [27]. Even smoothing provides serious pitfalls besides the introduced stiffness, since the numerical sensitivities are only correct if the step size is smaller than the smoothing parameter [36]. This requires large optimization problems even for moderate accuracy. Alternatively, high accuracy event-based methods [23,34] are, due to external switch detecting algorithms and internal logical statements, very difficult to incorporate into direct optimal control [23]. In this paper, we overcome these difficulties by applying the recently introduced Finite Elements with Switch Detection (FESD) method [30] to the time-freezing Filippov system. This method automatically detects the switching events without an external routine and delivers correct numerical sensitivities [30]. It was demonstrated in several benchmarks [31,28,30] that FESD outperforms standard discretization and mixed-integer approaches while providing several orders of magnitude more accurate for the same computational effort.

The time-freezing reformulation was first introduced in [29] for CLS with partially elastic impacts. In this paper, we extend these ideas to the inelastic case. Time-freezing was later independently introduced by Halm and Posa [18]. They treat the inherent nonuniqueness in the simulation of CLS with multiple frictional impacts in a stochastic way. The multiple outcomes are achieved by randomly varying parameters of the auxiliary dynamics and by considering a set of possible simulation solutions. In their work, time-freezing is used to obtain a differential inclusion with a bounded right-hand side, to be able to apply standard solution existence results. In this paper, our goal is to transform an Optimal Control Problem (OCP) with CLS into an equivalent OCP subject to a Filippov system, which we can solve with high accuracy via FESD[30].

Other reformulation approaches use coordinate transfor-
mations [22,42], penalization/smoothing of the comple-
mentary [36,38] and compliant contact models [7, Chap-
ter 2]. The former have the advantage that they are ex-
act, i.e., we can recover the solution of the original sys-

tem. However, they are usually limited to special cases,
namely, partially elastic impacts and a single scalar con-
straint. Examples are the Zhuravlev-Ivanov transforma-
tion [42],[7, Sec 1.4.3.] and the gluing function approach
[22] in the hybrid systems formalism.

1.1 Contributions

This paper extends the time-freezing reformulation from
[29], which exactly transforms a CLS with a single uni-
lateral constraint into a piecewise smooth system. We
discuss how to select auxiliary dynamics and formalize
the relationship between the time-freezing system and
CLS. Additionally, we present an extension to handle
state jumps in tangential directions resulting from fric-
tion. We demonstrate how to apply the time-freezing
system in optimal control and show that its solutions
are also optimal for the initial OCP with a CLS. More-
over, we show how to reformulate the time-freezing sys-
tem into a dynamic complementarity system, which fa-
cilitates the application of the FESD method for direct
optimal control. Furthermore, we introduce time trans-
formations and constraints to achieve equidistant con-
trol grid discretization and desired final time despite the
nonsmooth clock state. The discretized optimal control
problems result in mathematical programs with comple-
mentarity constraints. They are solved via a homotopy
approach. Thereby, we solve only a few smooth Nonlin-
ear Programs (NLP) and recover a highly accurate non-
smooth solution approximation with state jumps. The
theoretical considerations and efficacy of the proposed
numerical methods are demonstrated on a challenging
OCP example where the dynamic trajectory of a robot
hopping over holes is computed. All methods and exam-
ples from this paper, including a fully-automated refor-
mulation of the CLS into a PSS, are implemented in the
open-source tool NOSNOC [1,28].

1.2 Outline

The remainder of this paper is structured as follows.
Section 2 gives an introduction to piecewise smooth sys-
tems, their embedding into Filippov's framework, and
complementarity Lagrangian systems with state jumps.
This is followed by Section 3 where the time-freezing re-
formulation is discussed in detail. In Section 4 we ex-
tend these ideas to the frictional impact case. Section
5 discusses numerical optimal control for time-freezing
systems. Section 6 provides a numerical example of a
robotics optimal control problem. We conclude and list
some future research directions in Section 7.

1.3 Notation

Time derivatives of a function $x(t)$ w.r.t. to the phys-
ical time $t$ are compactly denoted by $\dot{x}(t) := \frac{dx(t)}{dt}$,
and of a function $y(\tau)$ w.r.t. to the numerical time
$\tau$ by $y'(\tau) := \frac{dy(\tau)}{d\tau}$. For the left and the right lim-
its, we use the notation $x(t_s^+) = \lim_{t \to t_s, t > t_s} x(t)$ and
$x(t_s^-) = \lim_{t \to t_s, t < t_s} x(t)$, respectively. For ease of notation, when clear from the context we drop the $t$, $\tau$, or $x$-dependencies. All vector inequalities are to be understood element-wise. The complementarity conditions for two vectors $a, b \in \mathbb{R}^n$ read as $0 \le a \perp b \ge 0$, where $a \perp b$ means $a^\top b = 0$. The matrix $I_n \in \mathbb{R}^{n \times n}$ is the identity matrix, and $0_{m,n} \in \mathbb{R}^{m \times n}$ is a matrix whose entries are all zeros. The concatenation of two column vectors $a \in \mathbb{R}^m, b \in \mathbb{R}^n$ is denoted by $(a, b) := [a^\top, b^\top]^\top$. The concatenation of several column vectors is defined analogously. A vector with all ones is denoted by $e = (1, 1, \dots, 1) \in \mathbb{R}^n$, and its dimension is clear from the context. The closure of a set $X$ is denoted by $\bar{X}$, its boundary by $\partial X$.

The set-valued sign function is defined as

$$ \operatorname{sign}(x) = \begin{cases} \{1\}, & x > 0, \\ [-1, 1], & x = 0, \\ \{-1\}, & x < 0. \end{cases} $$

The vector-valued version $\operatorname{Sign}: \mathbb{R}^n \to \mathbb{R}^n$ is defined as $\operatorname{Sign}(x) = (\operatorname{sign}(x_1), \dots, \operatorname{sign}(x_n))$. Key symbols and definitions used in this paper are summarized in Table 1.

## 2 Nonsmooth differential equations

In this section, we define piecewise smooth systems (PSS), Filippov's notion of solutions for the PSS [14], and complementarity Lagrangian systems (CLS).

### 2.1 Piecewise smooth systems

We regard piecewise smooth systems of the following form:

$$ \dot{x} = f_i(x, u), \text{ if } x \in R_i \subset \mathbb{R}^{n_x}, i \in \mathcal{I} := \{1, \dots, n_f\}, \quad (1) $$

with regions $R_i \subset \mathbb{R}^{n_x}$ and associated dynamics $f_i(\cdot)$, which are at least twice continuously differentiable functions on an open neighborhood of $\bar{R}_i$. The control function $u \in \mathbb{R}^{n_u}$ is assumed to be given and can be obtained, e.g., by solving an optimal control problem. The right-hand side of (1) is in general discontinuous in $x$. We assume that the sets $R_i$ are disjoint, nonempty, connected, and open. They have piecewise smooth boundaries $\partial R_i$. Moreover, it is assumed that $\bigcup_{i \in \mathcal{I}} \bar{R}_i = \mathbb{R}^n$ and that $\mathbb{R}^n \setminus \bigcup_{i \in \mathcal{I}} R_i$ is a set of measure zero.

The ODE (1) is not properly defined on the boundaries $\partial R_i$. To have a meaningful solution concept for the PSS (1) we regard its Filippov extension [14]. The ODE (1) is replaced by a differential inclusion whose r.h.s. is a convex and bounded set. Due to the assumed structure

<table><caption>Table 1
Key symbols used throughout this paper.</caption><thead><tr><th>Symbol</th><th>Meaning and reference</th></tr></thead><tbody><tr><td>$n_x$</td><td>dimension of the differential state x, Sec. 2</td></tr><tr><td>$n_f$</td><td>number of regions/modes in the PSS, Eq.(1)</td></tr><tr><td>$n_q$</td><td>dimension of the states q and v, Sec. 2.2</td></tr><tr><td>$n_u$</td><td>dimension of the control function, Sec. 2.2</td></tr><tr><td>$n_y$</td><td>dimension of the time-freezing state, Sec. 3.2</td></tr><tr><td>$n_t$</td><td>dimension of the tangent space, Sec. 4.1</td></tr><tr><td>$x$</td><td>differential state, Sec. 2.</td></tr><tr><td>$\theta$</td><td>Filippov's convex multipliers, Eq. (2)</td></tr><tr><td>$q$</td><td>position state, Eq. (3)</td></tr><tr><td>$v$</td><td>velocity state, Eq. (3)</td></tr><tr><td>$v_t$</td><td>tangential velocity at contact points, Sec. 4.2</td></tr><tr><td>$u$</td><td>control function, Eq. (3)</td></tr><tr><td>$\lambda_n$</td><td>Lagrange multiplier, norm. contact force, Eq. (3)</td></tr><tr><td>$\lambda_t$</td><td>Lagrange multiplier, friction force, Eq. (20)</td></tr><tr><td>$y$</td><td>extended state of time-freezing system, Sec. 3.2</td></tr><tr><td>$t$</td><td>physical time, Sec. 3.2</td></tr><tr><td>$\tau$</td><td>numerical time, Sec. 3.2</td></tr><tr><td>$a_n$</td><td>constant of the auxiliary dynamics, Prop. 3</td></tr><tr><td>$a_t$</td><td>constant of the auxiliary dynamics, Eq. (24)</td></tr><tr><td>$\lambda^n, \lambda^p$</td><td>dual variables of the linear program, Eq. (31)</td></tr><tr><td>$\alpha$</td><td>primal variables of the linear program, Eq. (30)</td></tr><tr><td>$\mu$</td><td>coefficient of friction, Eq. (20f)</td></tr><tr><td>$\epsilon_t$</td><td>relaxation parameter in friction model, Eq. (26)</td></tr><tr><td>$s$</td><td>speed of time control variable, Eq. (33)</td></tr><tr><td>$f_i(x, u)$</td><td>modes of the PSS system, Eq. (1)</td></tr><tr><td>$f_v(q, v, u)$</td><td>vector field of the velocity state, Eq. (3b)</td></tr><tr><td>$n(q)$</td><td>normal to the CLS constraint surface, Eq. (3b)</td></tr><tr><td>$M(q)$</td><td>inertia matrix, Eq. (3b)</td></tr><tr><td>$f_c(q)$</td><td>constraint function in the CLS, Eq. (3c)</td></tr><tr><td>$f_{ODE}(x, u)$</td><td>unconstrained dynamics when $f_c(q) > 0$, Eq. (6)</td></tr><tr><td>$D(q)$</td><td>Delassus' matrix/scalar, Eq. (11a)</td></tr><tr><td>$\varphi(x, u)$</td><td>determines if contact persists, Eq (11b)</td></tr><tr><td>$f_{DAE}(x, u)$</td><td>dynamics equivalent to the DAE (7) (constrained dynamics when $f_c(q) = 0$), Eq. (13)</td></tr><tr><td>$c_i(y)$</td><td>time-freezing PSS switching functions, Sec. 3.2, Eq. (14) and Sec. 4.3, Eq. (25)</td></tr><tr><td>$f_{aux,n}(y)$</td><td>auxiliary dynamics, Def. 1 and Prop. 3</td></tr><tr><td>$\gamma(x, u)$</td><td>time-rescaling factor of the time-freezing sliding mode, Eq. (17)</td></tr><tr><td>$B(q)$</td><td>matrix whose columns span the tangent space at contact points, Eq. (20)</td></tr><tr><td>$b_j(q)$</td><td>j-th column of $B(q)$, Eq. (20)</td></tr><tr><td>$f_{aux,t}^-(y)$</td><td>auxiliary dynamics for tan. directions, Eq. (24)</td></tr><tr><td>$f_{Slip}(x, u)$</td><td>dynamics for slipping motion in contact phases, Eq. (22)</td></tr><tr><td>$f_{Stick}(x, u)$</td><td>dynamics for sticking motion in contact phases, Eq. (23)</td></tr><tr><td>$\tilde{D}(q)$</td><td>generalization of $D(q)$, Sec. 4.2</td></tr><tr><td>$\tilde{\varphi}(x, u)$</td><td>generalization of $\varphi(x, u)$, Sec. 4.2</td></tr><tr><td>$g_F(\theta, \alpha)$</td><td>expression for relating $\theta$ and $\alpha$, Eq. (32)</td></tr><tr><td>$\Psi(x(T))$</td><td>terminal cost, Eq. (27)</td></tr><tr><td>$g(x, u)$</td><td>path and terminal constraints, Eq. (27)</td></tr><tr><td>$r(x)$</td><td>terminal constraints, Eq. (27)</td></tr><tr><td>$R_i$</td><td>regions of the PSS, Eq. (1)</td></tr><tr><td>$\mathcal{I}$</td><td>index set for PSS modes, Eq. (1)</td></tr><tr><td>$F_F(x, u)$</td><td>Filippov set, Eq. (2)</td></tr><tr><td>$\Sigma$</td><td>switching surface, Sec. 3.2</td></tr><tr><td>$F_{TF}(y, u)$</td><td>Filippov set for the time-freezing system, Def. 2 and 6</td></tr><tr><td>$Q$</td><td>region where all auxiliary dynamics are defined, Eq. (25)</td></tr></tbody></table>
of the sets $R_i$, if $\dot{x}$ exists, functions $\theta_i(\cdot)$ that serve as convex multipliers, can be introduced, and the Filippov differential inclusion for (1) reads as [14,34]:

$$ \dot{x} \in F_F(x, u) := \left\{ \begin{array}{l} \sum_{i \in \mathcal{I}} f_i(x, u) \theta_i \mid \sum_{i \in \mathcal{I}} \theta_i = 1, \theta_i \ge 0, \\ 0 = \theta_i \text{ if } x \notin \bar{R}_i, \forall i \in \mathcal{I} \end{array} \right\}. \quad (2) $$

Note that in the interior of the regions $R_i$ the Filippov set $F_F(x, u)$ is equal to $\{f_i(x, u)\}$ and on the boundary between regions it is a convex combination of the neighboring vector fields. Sufficient conditions for the existence and uniqueness of solutions are given in [14]. We assume that the boundaries of the regions $\partial R_i$ are defined by the zero-level sets of, at least twice continuously differentiable, scalar functions $c_i(x) = 0$. These functions are called switching functions. The evolution of $x(\cdot)$ on region boundaries $\partial R_i$, are called sliding modes. The dynamics of sliding modes are implicitly defined by differential algebraic equations, since the corresponding constraints $c_i(x) = 0$ must hold [14].

## 2.2 Complementarity Lagrangian systems

This paper regards complementarity Lagrangian systems (CLS) with a single unilateral constraint. A CLS with inelastic impacts reads as

$$ \dot{q} = v, \quad (3a) $$

$$ \dot{v} = f_v(q, v, u) + M(q)^{-1}n(q)\lambda_n, \quad (3b) $$

$$ 0 \le \lambda_n \perp f_c(q) \ge 0, \quad (3c) $$

$$ 0 = n(q(t_s))^{\top} v(t_s^{+}), \quad (3d) $$

$$ \text{if } f_c(q(t_s)) = 0 \text{ and } n(q(t_s))^{\top} v(t_s^{-}) < 0, $$

where $x := (q, v)$ and the states $q \in \mathbb{R}^{n_q}$ and $v \in \mathbb{R}^{n_q}$ correspond to the position and velocity of a rigid body, respectively. The variable $\lambda_n \in \mathbb{R}$ is the Lagrange multiplier and is physically interpreted as the normal contact force. The function $f_c(q) \in \mathbb{R}$ is the signed distance between a contact point of a rigid body and an obstacle or another rigid body, and $n(q) := \nabla_q f_c(q)$. The matrix $M(q)$ is the inertia matrix and the function $f_v(q, v, u)$ corresponds to the total acceleration of the rigid body, i.e., it collects all internal and external forces (except the contact forces), multiplied by the inverse of $M(q)$. The functions $f_v: \mathbb{R}^{n_q} \times \mathbb{R}^{n_q} \times \mathbb{R}^{n_u} \to \mathbb{R}^{n_q}$, $M(q): \mathbb{R}^{n_q} \to \mathbb{R}^{n_q \times n_q}$, $f_c: \mathbb{R}^{n_q} \to \mathbb{R}$ are assumed to be at least twice continuously differentiable and the matrix $M(q)$ is assumed to be symmetric positive definite. The complementarity condition (3c) states that: either the system is in contact and there is a reaction force ($f_c(q) = 0, \lambda_n \ge 0$) or there is no contact and no reaction force ($f_c(q) > 0, \lambda_n = 0$). When the body makes contact, the negative normal velocity $n(q)^{\top}v$ must jump to zero. This is modeled via Eq. (3d).

The function $u \in \mathbb{R}^{n_u}$ in Eq. (3b) is the control function. In this paper, we aim to find the control function $u(t)$ by solving the Optimal Control Problem (OCP):

$$ \min_{x(\cdot), \lambda_n(\cdot), u(\cdot)} \Psi(x(T)) \quad (4a) $$

$$ \text{s.t.} \quad x(0) = \bar{x}_0, \quad (4b) $$

$$ \text{Eq. (3), } t \in [0, T] \quad (4c) $$

$$ 0 \le g(x(t), u(t)), t \in [0, T], \quad (4d) $$

$$ 0 \le r(x(T)), \quad (4e) $$

where $\Psi: \mathbb{R}^{n_x} \to \mathbb{R}$ is the terminal cost of the OCP, $\bar{x}_0$ is a given initial value. The functions $g: \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \to \mathbb{R}^{n_g}$ and $r: \mathbb{R}^{n_x} \to \mathbb{R}^{n_r}$ are the path and terminal constraints, respectively. The CLS dynamics in Eq. (4c) (resp. Eq. (3)) make this OCP nonsmooth and nonconvex.

The nonsmooth dynamics and discontinuous velocity state make this OCP difficult to solve numerically. In the subsequent sections, we develop the time-freezing reformulation such that we can reformulate the problem (4) into an OCP subject to a Filippov system, for which more efficient numerical methods are available [30]. In this paper, we focus on the case of a single unilateral constraint $f_c(q)$. Extensions for multiple and simultaneous impacts will be studied in future work. The extension depends on the chosen impact model [7].

Next, we introduce a guiding example on which we will illustrate the main ideas behind the time-freezing reformulation throughout the paper.

**Example 1 (Guiding example)** Consider a frictionless point mass in two dimensions above a horizontal table. The mass is $m = 1$ kg and $g = 9.81 \text{ m/s}^2$ is the gravitational acceleration. Denote by $q := (q_1, q_2)$ and $v := (v_1, v_2)$ its position and velocity, respectively, and let $\lambda_n$ be the normal contact force. The dynamics are given by the CLS:

$$ \dot{q} = v, \quad (5a) $$

$$ m\dot{v} = \begin{bmatrix} 0 \\ -mg \end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} \lambda_n + \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}, \quad (5b) $$

$$ 0 \le \lambda_n \perp q_2 \ge 0, \quad (5c) $$

$$ v_2(t_s^+) = 0, \text{ if } q_2(t_s) = 0 \text{ and } v_2(t_s^-) < 0. \quad (5d) $$

where $u = (u_1, u_2) \in \mathbb{R}^2$ is an externally chosen thrust force, which shall be found, e.g., by solving an optimal control problem.

# 3 The time-freezing reformulation

This section develops the time-freezing reformulation that enables one to transform the CLS (3) with inelastic impacts into a PSS of the form of (1). For simplicity
and ease of exposition, we first focus on the case without friction. Extensions with frictional impacts are given in Section 4. We start by investigating different possible modes of the CLS (3). Afterward, the time-freezing reformulation with its needed ingredients is introduced. The section finishes by formally relating the two regarded systems.

## 3.1 The different modes of the CLS

For the CLS (3) we can distinguish two modes of operation:

(i) the system is not in contact (unconstrained case, free flight), i.e., $f_c(q) > 0$ which implies $\lambda_n = 0$, (ii) the system is in contact, i.e., $f_c(q) = 0$ and $\lambda_n \ge 0$. In the first case, the system evolves according to the ODE:

$$ \dot{q} = v, \quad \dot{v} = f_v(q, v, u). \qquad (6) $$

We write this ODE compactly as $\dot{x} = f_{ODE}(x, u) := (v, f_v(q, v, u))$.

We call an active-set change from $\lambda_n(t_s^-) = 0$, $f_c(q(t_s^-)) \ge 0$ to $\lambda_n(t_s^+) \ge 0$, $f_c(q(t_s^+)) = 0$, which triggers a state jump, an *impact*.

After an impact, it holds that $0 = n(q(t_s))^{\top}v(t_s^{\perp})$. Subsequently, the system evolves according to a differential algebraic equation (DAE) of index 3:

$$ \dot{q} = v, \qquad (7a) $$

$$ \dot{v} = f_v(q, v, u) + M(q)^{-1}n(q)\lambda_n, \qquad (7b) $$

$$ 0 = f_c(q). \qquad (7c) $$

Note that $f_c(q(t))$ needs to be differentiated twice w.r.t. to time until $\lambda_n(t)$ appears explicitly. The next question to be answered is: will the system stay in contact (dynamics defined by (7) with $f_c(q) = 0$) or will the contact break (dynamics defined by (6) with $f_c(q) > 0$)? The answer can be found by looking at the contact Linear Complementarity Problem (LCP) [7, Section 5.1.2]. Under our standing assumptions, during contact on some time interval $[t_1, t_2]$ the consistent initialization conditions hold

$$ 0 = f_c(q(t)), \quad 0 = \frac{d}{dt} f_c(q(t)) = \nabla_q f_c(q(t))^{\top} v(t). \qquad (8) $$

Consequently, $\lambda_n(t) \ge 0, t \in [t_1, t_2]$. Due to the continuity of $q(t)$, $f_c(q(t))$ and $\frac{d}{dt}f_c(q(t))$, for contact breaking (i.e., $f_c(q)$ becomes strictly positive) it is required that $\frac{d^2}{dt^2}f_c(q(t)) \ge 0$ for $t \in [t_2, t_2 + \hat{\epsilon})$, for some $\hat{\epsilon} > 0$.

Therefore, from (3c) we deduce that

$$ 0 \le \frac{d^2}{dt^2} f_c(q(t)) \perp \lambda_n(t) \ge 0, \quad t \in [t_1, t_2]. \qquad (9) $$

Then, by computing $\frac{d^2}{dt^2}f_c(q(t))$ and using the r.h.s. of (7b), we obtain the contact LCP in $\lambda_n(t)$:

$$ 0 \le D(q)\lambda_n + \varphi(x, u) \perp \lambda_n \ge 0, \qquad (10) $$

with

$$ D(q) = \nabla_q f_c(q)^{\top} M(q)^{-1} \nabla_q f_c(q), \qquad (11a) $$

$$ \varphi(x, u) = \nabla_q f_c(q)^{\top} f_v(q, v, u) + \nabla_q (\nabla_q f_c(q)^{\top} v)^{\top} v, \qquad (11b) $$

where $D(q) > 0$ [7].

The solution map of the LCP (10) is given by

$$ \lambda_n = \max(0, -D(q)^{-1}\varphi(x, u)). \qquad (12) $$

From the last equation we deduce that contact breaking or sticking depends on the sign of the function $\varphi(x, u)$.

In the case of $\varphi(x, u) \le 0$ from Eq. (10) and (12) it follows that $\lambda_n(t) \ge 0$ and $\frac{d^2}{dt^2}f_c(q(t)) = 0$. Therefore, we have a persistent contact and the system evolves according to the DAE (7). Using index reduction and the solution map (12) we can derive an ODE that is equivalent to the DAE (7):

$$ \dot{q} = v, \qquad (13a) $$

$$ \dot{v} = f_v(q, v, u) - M(q)^{-1}n(q)D(q)^{-1}\varphi(x, u). \qquad (13b) $$

We compactly denote this ODE by $\dot{x} = f_{DAE}(x, u)$.

In the second case, $\varphi(x, u) > 0$ implies $\lambda_n(t) = 0$ and $\frac{d^2}{dt^2}f_c(q(t)) > 0$, therefore the contact breaks and the system evolves according to the ODE (6).

To summarize, if the system switches from the ODE mode in (6) to the DAE mode in (13), a state jump must occur, except if the active-set changes happen with $n(q(t_s^-))^{\top}v(t_s^-) = 0$. Now the system evolves on the boundary of the feasible set with $f_c(q) = 0$ according to the DAE (7), or equivalently according to the ODE defined by (13). On the other hand, if we switch from DAE to ODE mode, we have a continuous transition without state jumps, i.e., contact breaking occurs.

## 3.2 Main ideas and auxiliary dynamics

The arguments above reveal that the CLS (3) switches between an ODE and a DAE of index 3. This already
bears similarity to a PSS, but the main obstacle to completing this transition are the state jumps. Note that large parts of the state space, namely $f_c(q) < 0$, are prohibited for the solution trajectories of the CLS.

The time-freezing reformulation is based on the following two main ideas [29]. First, we relax the constraint and allow $f_c(q) < 0$. We define an auxiliary dynamical system in this infeasible region whose trajectory endpoints satisfy the state jump law (3d) on some finite time interval. Second, we introduce a clock state $t(\tau)$ that stops counting (i.e., $t'(\tau) = 0$), when the auxiliary ODE is active (for $f_c(q) < 0$). By taking the pieces of the trajectory when the clock state was active, one can recover the solution of the original system with discontinuous trajectories. Note that the time-freezing system has no discontinuities in its solution, but only in its r.h.s.. The extended state of the time-freezing system reads as $y := (x, t) \in \mathbb{R}^{n_y}$, $n_y = n_x + 1$. The time of the time-freezing system $\tau$ is called numerical time. The intervals with $t'(\tau) > 0$ are referred to as physical time and those with $t'(\tau) = 0$ as virtual time. The properties of the auxiliary dynamics are summarized in the following definition. Later a constructive way to find such systems is given.

**Definition 1 (Auxiliary dynamics)** An auxiliary dynamical system $y'(\tau) = f_{aux,n}(y(\tau))$ satisfies for every initial value $y(\tau_s) = (q_s, v_s, t_s)$, with $f_c(q_s) = 0$ and $n(q_s)^\top v_s < 0$, for every well-defined and finite time interval $(\tau_s, \tau_r)$, with the length $\tau_{jump} = \tau_r - \tau_s$, the following properties: (i) $f_c(q(\tau)) \le 0$, $t'(\tau) = 0$ $\forall \tau \in (\tau_s, \tau_r)$, (ii) $n(q(\tau_r))^\top v(\tau_r) = 0$, and (iii) $f_c(q(\tau_r)) = 0$.

To construct the time-freezing system we take several steps. First, observe that the post-impact velocity (3d) is equal to the total time derivative of the constraint $f_c(q) = 0$, i.e., $\frac{d}{dt}f_c(q) = n(q)^\top v = 0$. We choose these functions as switching functions, i.e., $c_1(y) = f_c(q)$ and $c_2(y) = n(q)^\top v$, and define the following regions:

$$R_1 = \{y \in \mathbb{R}^{n_y} | c_1(y) > 0\} \cup \{y \in \mathbb{R}^{n_y} | c_1(y) < 0, c_2(y) > 0\}, \\ R_2 = \{y \in \mathbb{R}^{n_y} | c_1(y) < 0, c_2(y) < 0\}. \qquad (14)$$

Second, we associate with the region $R_1$ the unconstrained dynamics $y' = (f_{ODE}(x, u), 1)$ and with $R_2$ the auxiliary dynamics from Definition 1. Note that the unconstrained dynamics is Eq. (6) augmented by the clock state dynamics $t' = 1$. Moreover, the control functions $u(\cdot)$ does not influence the dynamics for $y \in R_2$, i.e., whenever the time is time frozen the control can take any value.

We formally define the time-freezing system and its Filippov extension in the next definition.

**Definition 2 (Time-freezing system)** Let $\tau \in \mathbb{R}$ be the numerical time and $y(\tau) := (x(\tau), t(\tau)) \in \mathbb{R}^{n_y}$ the

<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration of a phase plot of the time-freezing system from Definition 2. The red and yellow shaded areas are infeasible for the CLS (3). The trajectories of the auxiliary dynamics (the blue dashed line) flow in the red-shaded area.</figcaption></figure>

differential states and $u(\tau) \in \mathbb{R}^{n_u}$ a given control function. The time-freezing PSS is a PSS of the form of Eq. (1), defined by the regions $R_1$ and $R_2$ in (14) with $f_1(y, u) = (f_{ODE}(x, u), 1)$ and $f_2(y) = f_{aux,n}(y)$. The corresponding Filippov system, which we call the time-freezing system, is defined as

$$y' \in F_{TF}(y, u) := \left\{ \begin{array}{l} \theta_1 f_1(y, u) + \theta_2 f_2(y) \mid e^{\top} \theta = 1, \theta \ge 0 \\ \theta_i = 0, \text{ if } y \notin \bar{R}_i, i = 1, 2 \end{array} \right\}, \qquad (15)$$

with $\theta = (\theta_1, \theta_2)$. It is assumed that appropriate dynamics $f_{aux,n}(y)$ exist.

The phase plot of the time-freezing system is depicted in Fig. 1. Note that the region $R_1$ consists of $R_1^a = \{y \in \mathbb{R}^{n_y} | f_c(q) > 0\}$ (green area), that corresponds to the feasible set of the unconstrained dynamics (6) and the set $R_1^b := \{y \in \mathbb{R}^{n_y} | c_1(y) < 0, c_2(y) > 0\}$ (yellow area). The solution trajectories never flow in $R_1^b$ and the system should not be initialized in there. However, as we show later, it is crucial for sliding modes and contact breaking. Region $R_2$ (red shaded area) contains the auxiliary dynamics that mimic the state jump.

To make further use of Definition 2 we must specify how to select an appropriate auxiliary ODE. The next proposition provides a constructive way of selecting the auxiliary ODE from Definition 1 for any smooth scalar constraint $f_c(q) = 0$.

**Proposition 3 (Auxiliary dynamics)** Suppose that $y(\tau_s) = (q_s, v_s, t_s)$ is given such that $f_c(q_s) = 0$ and
$n(q_s)^\top v_s \le 0$ holds. Then the ODE given by

$$y' = f_{\text{aux,n}}(y) := \begin{bmatrix} \mathbf{0}_{n_q, 1} \\ M(q)^{-1}n(q)a_n \\ 0 \end{bmatrix} \quad (16)$$

with $a_n > 0$ is an auxiliary dynamical system from Definition 1 with $\tau_{\text{jump}} = -\frac{n(q_s)^\top v_s}{D(q_s)a_n}$.

**PROOF.** According to (16) we have $q'(\tau) = \mathbf{0}_{n_q, 1}$, $\forall \tau \ge \tau_s$, which implies $q(\tau) = q_s$ and $f_c(q(\tau)) = 0$, $\forall \tau \ge \tau_s$. This means also that $M(q(\tau)) = M(q_s)$, $\forall \tau \ge \tau_s$. Second, regard the dynamics of $v' = M(q)^{-1}n(q)a_n = M(q_s)^{-1}n(q_s)a_n$ and rewrite this equation in integral form. By multiplying it from the left by $n(q_s)^\top$ we obtain:

$$n(q_s)^\top v(\tau) = n(q_s)^\top v_s + n(q_s)^\top M(q_s)^{-1}n(q_s)a_n(\tau - \tau_s).$$

Since the first term on the r.h.s. is negative and the second strictly positive, we deduce that $n(q(\tau_r))^{\top} v(\tau_r) = 0$ and $f_c(q(\tau_r)) = 0$ with $\tau_{\text{jump}} = \tau_r - \tau_s = -\frac{n(q_s)^\top v_s}{D(q_s)a_n}$. Hence, all conditions from Definition 1 are satisfied and the proof is complete. $\square$

Next, we discuss which mode of the time-freezing system matches the persistent contact dynamics of the CLS (3). We observe that the set $\Sigma := \{y | c_1(y) = f_c(q) = 0, c_2(y) = n^\top v = 0\}$ is defined by the same equations as the consistent initialization conditions (8) (but now in $\mathbb{R}^{n_y}$ instead of $\mathbb{R}^{n_x}$ due to the clock state). A sliding mode of the time-freezing PSS evolves on $\Sigma$ just as the solution of the persistent contact DAE (7). Therefore, its dynamics should match the dynamics of the DAE (7). We detail in the next two subsections that this is indeed the case.

To summarize, for $f_c(q) > 0$ (which is a subset of $R_1$) the time-freezing system and CLS have the same dynamics. In $R_2$ the auxiliary dynamics mimics the state jump and the sliding mode $y \in \Sigma$ should match the dynamics of the CLS in contact mode. To illustrate the developments so far, we derive a time-freezing PSS for our guiding example.

**Example 2 (Guiding example as time-freezing PSS)**
The state space is $y = (q, v, t) \in \mathbb{R}^5$ and the switching functions read as $c_1(y) = y_2$ and $c_2(y) = v_2$. The two PSS regions are $R_1 = \{y \in \mathbb{R}^5 | q_2 > 0\} \cup \{y \in \mathbb{R}^5 | q_2 < 0, v_2 > 0\}$ and $R_2 = \{y \in \mathbb{R}^5 | q_2 < 0, v_2 < 0\}$, cf. Fig. 1. The PSS dynamics are given by $f_1 = (v_1, v_2, u_1, -g + u_2, 1)$. The constraint normal is $n = (0, 1)$, so that by applying Eq. (16) we find that $f_2 = (0, 0, 0, a_n, 0)$.

### 3.3 Persistent contact and sliding mode

Depending on the sign of the function $\varphi(x, u)$, a solution $y$ initialized at $\Sigma$ should either stay at $\Sigma$ (sliding mode, persistent contact) or leave it (contact breaking). In this subsection, we study the case when a solution of the time-freezing system satisfies the conditions $y(\tau) \in \Sigma$ and $\varphi(x(\tau), u(\tau)) \le 0$ for some $\tau \in [\tau_1, \tau_2]$ (persistent contact). During contact, the CLS system satisfies the consistent initialization (8) which corresponds to $\Sigma$ without the clock state. It is desired that under these conditions $y(\tau)$ stays on $\Sigma$ and that the corresponding sliding mode dynamics match the DAE dynamics (13). For a solution to stay in the sliding mode, the surface $\Sigma$ must be stable, i.e., all neighboring vector fields point toward $\Sigma$. Since $\nabla c_1(y)^\top f_1(y, u) = 0$, $\nabla c_1(y)^\top f_2(y) = 0$ and $\nabla c_2(y)^\top f_1(y, u) = \varphi(x, u) \le 0$, $\nabla c_2(y)^\top f_2(y) = D(q)a_n > 0$, we see that this is indeed the case, cf. Fig. 1. We show next that the sliding mode of the time-freezing system is unique and that it matches the dynamics of the DAE of index 3 after the state jump, as required.

**Theorem 4 (Unique sliding mode)**
Regard the time-freezing system from Definition 2 with the auxiliary dynamics from Proposition 3. Let $y(\tau)$ be a solution of this system with $y(0) \in \Sigma$ and $\tau \in [0, \tau_f]$. Suppose that $\varphi(x(\tau), u(\tau)) \le 0$ for all $\tau \in [0, \tau_f]$ (persistent contact), then the following statements are true:

(i) the convex multipliers $\theta_1, \theta_2 \ge 0$ in Eq. (15) are unique,

(ii) the dynamics of the sliding mode are given by $y' = \gamma(x, u)(f_{\text{DAE}}(x, u), 1)$, where $\gamma(x, u) \in (0, 1]$ is a time-rescaling factor given by

$$\gamma(x, u) := \frac{D(q)a_n}{D(q)a_n - \varphi(x, u)}. \quad (17)$$

**PROOF.** To compute the convex multipliers $\theta_1$ and $\theta_2$ we use Definition 2 and the fact that $y \in \Sigma$. This results in the conditions:

$$c_1(y) = 0, c_2(y) = 0, \theta_1 + \theta_2 = 1.$$

Since $\theta$ does not explicitly appear in the first two conditions, we differentiate them w.r.t. to $\tau$ and use Eq. (15). We have two unknowns and three conditions, hence the system is over-determined. However, we have by assumption that $c_2(y) = \nabla_q c_1(y)^\top v = 0$ and by direct evaluation, we conclude that $\frac{dc_1(y)}{d\tau} = 0$ is satisfied for every $\theta_1$ and $\theta_2$. The conditions that are left are $\frac{dc_2(y)}{d\tau} = 0$ and $\theta_1 + \theta_2 = 1$. Since $\nabla_q c_2(y) = \nabla_q (\nabla_q c_1(y)^\top v) = \nabla_q c_1(y)$ and $\frac{dc_2(y)}{dt} = 0$ we obtain from $\frac{dc_2(y)}{d\tau} = 0$ the following
equations

$$
\begin{align*}
0 = \frac{dc_2(y)}{d\tau} &= \theta_1 [\nabla_q c_2(y)^\top \nabla_v c_2(y)^\top] [f_v(q, v, u)] \\
&\quad + \theta_2 [\nabla_q c_2(y)^\top \nabla_v c_2(y)^\top] [M(q)^{-1} \nabla_q c_1(y) a_n], \\
0 &= \theta_1 \underbrace{\nabla_q (\nabla_q f_c(q)^\top v)^\top v + \nabla_q f_c(q)^\top f_v(q, v, u)}_{=\varphi(x,u)<0} \\
&\quad + \theta_2 \underbrace{\nabla_q f_c(q)^\top M(q)^{-1} \nabla_q f_c(q) a_n}_{=D(q)>0}.
\end{align*}
$$

Thus we obtain a system linear in θ

$$
\begin{bmatrix} \varphi(x, u) & D(q)a_n \\ 1 & 1 \end{bmatrix} \begin{bmatrix} \theta_1 \\ \theta_2 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix},
$$

and by solving it we have that $\theta_1 = \frac{D(q)a_n}{D(q)a_n - \varphi(x,u)} = \gamma(x, u)$ and $\theta_2 = \frac{-\varphi(x,u)}{D(q)a_n - \varphi(x,u)}$. Since $D(q)a_n - \varphi(x, u) > 0$ and $\varphi(x, u) \le 0$, we have always unique $\theta_1, \theta_2 \ge 0$. This completes the first part of the proof.

For the second part, we evaluate

$$
y' = \theta_1 f_1(y, u) + \theta_2 f_2(y) = \gamma(x, u) \begin{bmatrix} v \\ f_v(q, v, u) \\ 1 \end{bmatrix} + \frac{-\varphi(x, u)}{D(q)a_n - \varphi(x, u)} \begin{bmatrix} \mathbf{0}_{n_q, 1} \\ M(q)^{-1} \nabla f_c(q) a_n \\ 0 \end{bmatrix} \cdot D(q)^{-1} D(q).
$$

In the second term we use that $\lambda_n = -D(q)\varphi(x, u)$ (cf. Eq. (12)) and the expression for $\gamma(x, u)$ in Eq. (17). By comparing the last expression to Eq. (13) we obtain $y' = \gamma(x, u)(f_{\text{DAE}}(x, u), 1)$. This completes the proof. $\square$

This theorem shows that the sliding mode of the time-
freezing system on Σ is unique and equal to the dynamics
of the CLS in the persistent contact mode given by Eq.
(13), but slowed down by the factor γ(x, u). Note that
for larger values of a_n the factor γ(x, u) comes closer to
one, which reduces the slow-down, cf. Example 3.

However, in the function $t(\tau) \to x(t(\tau)))$ the solution
with a speed of time of one is recovered. We briefly dis-
cuss the intuition behind the time slow down by $\gamma(x, u)$.
To achieve the sliding mode on $\Sigma$, the vector fields from
$R_1$ and $R_2$ (the yellow part in Fig. 1) must push toward
$\Sigma$. The resulting dynamics is a convex combination of
the two vector fields, and since the speed of time in $R_1$
is one and in $R_2$ zero we obtain a slow down equal to
$\theta_1 = \gamma(x, u)$. Moreover, the vector field $f_{ODE}(\cdot)$ in the
yellow area "stops" the trajectory coming from $R_2$ and

<figure><img src="image_4.png" /><figcaption>Fig. 2. Illustration of a phase plot of the time-freezing system from Definition 2 with $\varphi(x, u) > 0$. Compared to Fig. 1, the vector field in $R_1$ is changed and $\Sigma$ is not stable anymore, thus leaving $\Sigma$ into $R_1$ is possible.</figcaption></figure>

thus enables the sliding mode. This shows the signifi-
cance of having the vector field ($f_{ODE}(\cdot)$, 1) in the yellow
area $R_1^b$, even though by construction the solution never
flows there.

3.4 Contact breaking

It is left to study the case when $y(\tau_1) \in \Sigma$ but
$\varphi(x(\tau), u(\tau)) > 0$ for $\tau \in [\tau_1, \tau_2]$. In the CLS for
$\varphi(x, u) > 0$ the contact breaks, cf. Sec 3.1. In the
time-freezing system, we expect the trajectory to leave
the sliding mode from $\Sigma$. Since $\nabla c_1(y)^\top f_1(y, u) = 0$,
$\nabla c_1(y)^\top f_2(y) = 0$ and $\nabla c_2(y)^\top f_1(y, u) = \varphi(x, u) > 0$,
$\nabla c_2(y)^\top f_2(y) = D(q)a_n > 0$, the surface $\Sigma$ is not sta-
ble anymore. This scenario is illustrated in Fig. 2. We
conclude that under these conditions $y$ leaves $\Sigma$ and en-
ters $R_1$ (into the green region with $f_c(q) > 0$). In this
case $\theta = (1,0)$, hence $y' = f_1(y, u) = (f_{ODE}(x, u), 1)$,
which matches the unconstrained CLS dynamics (6) aug-
mented by the clock state.

Effectively, the time-freezing system switches between
the DAE and ODE modes, just as the CLS and the state
jump is performed by the auxiliary dynamics while the
time is frozen. For $\varphi(x, u) \le 0$ (persistent contact) the
solution of the time-freezing system stays on $\Sigma$, just as
the solution of the CLS. It leaves the sliding mode when
$\varphi(x, u) > 0$, which corresponds to contact breaking in
the CLS. This relationship is formalized in the next sub-
section. To illustrate the developments of the last two
subsections, we revisit our guiding example and provide
a simulation that encompasses all effects discussed so far.

**Example 3 (Speed of time and sliding modes)** Let us consider the time-freezing PSS from Example 2. We choose $a_n = g$. It follows that $D(q) = 1$ and $\varphi(x, u) =$
$-g + u_2$. We choose a control function

$$u(t) = \begin{cases} (7, 0), & t < 1, \\ (7, 2g(t(\tau) - 1)), & t \ge 1. \end{cases}$$

Let us make a simulation of the time-freezing system with $y(0) = (0, 1, 0, 0, 0)$ for $\tau \in [0, 3.5]$. The result is depicted in Fig. 3. The particle hits the ground, slides horizontally on it, and lifts off when the control force $u_2(t)$ is stronger than gravity, cf. the top plot. We see that when the particle hits the ground, the time is frozen and the auxiliary ODE is active (red strips). The vertical velocity $v_2$ becomes zero with the rate $a_n$. The system is then a sliding mode with a time slow down factor of $\gamma(x, u) = \frac{a_n}{a_n - (-g + u_2)} = \frac{g}{g - (-g + 0)} = 0.5$ for $t < 1$ ($\tau < 2$), cf. bottom left plot. At $t = 1$, which corresponds to $\tau = 2$ the vertical control force becomes nonzero and $\gamma(x, u) = \frac{g}{2g + u_2}$ grows. For $\tau > 2.8$ we have $\varphi(x(\tau), u(\tau)) = -g + u_2 > 0$ and $\Sigma$ is not stable anymore. The particle lifts off and the contact breaks. Note that the solution of the time-freezing system $y(\cdot)$ is continuous in numerical time $\tau$ (middle left plot) and discontinuous in physical time $t$ (middle right plot).

## 3.5 Solution relationship

We formalize now how to recover the solution of the initial value problem corresponding to the CLS (3) from the solution of the time-freezing system from Definition 2.

**Theorem 5 (Solution relationship)** Regard the initial value problems corresponding to: i) the time-freezing system in Eq. (15) with a given $y(0) = (q_0, v_0, 0) \in \mathbb{R}^{n_y}$ and $f_c(q_0) \ge 0$ on a time interval $[0, \tau_f]$, and ii) the CLS from Eq. (3) with the initial value $x(0) = (q_0, v_0) \in \mathbb{R}^{n_x}$ on a time interval $[0, t_f] := [0, t(\tau_f)]$, with $f_c(q(t_f)) \ge 0$

<figure><img src="image_5.png" /><figcaption>Fig. 3. Trajectories of the time-freezing system corresponding to the example CLS. The top plot shows the position of the particle. The middle left plot shows the continuous velocities $v_1$ and $v_2$ in numerical time $\tau$ and the middle right plot shows the discontinuous velocities $v_1$ and $v_2$ in physical time $t$. The bottom left plot shows the speed of time $\frac{dt}{d\tau}$ and the bottom right plot the control function $u(t)$.</figcaption></figure>

and $n(q(t_f))^{\top}v(t_f) \ge 0$. Suppose the following assumptions hold:

(a) the auxiliary dynamics $f_{aux,n}(y)$ from Proposition 3 is used in the time-freezing system in Definition 2,

(b) there is at most one time point $t_s = t(\tau_s)$ where $f_c(q(t_s)) = 0$ and $n(q(t_s))^{\top}v(t_s^{-}) < 0$ on the time interval $[0, t_f]$,

Then, the solutions to the two problems are related as follows:

(1) for $t \ne t_s$:

$$\begin{align} x(t(\tau)) &= Ry(\tau), \quad \text{with } R = [I_{n_x} \quad \mathbf{0}_{n_x,1}], \tag{18a} \\ \lambda_n(t(\tau)) &= \begin{cases} -D(q(t(\tau)))\varphi(x(t(\tau))), & \text{if } y(\tau) \in \Sigma, \\ 0, & \text{otherwise.} \end{cases} \tag{18b} \end{align}$$

(2) for $t = t_s$:

$$\lim_{\epsilon \to 0} \lim_{\tau \to \tau_s} \int_{t_s - \epsilon}^{t_s + \epsilon} \lambda_n(t) dt = \int_{\tau_s}^{\tau_r} a_n d\tau. \quad (19)$$

**PROOF.** See Appendix A. □

Equivalence for several subsequent impacts is trivially obtained by sequentially applying the argument of the last theorem. Time-freezing enables one to make the state jump in "slow motion". By plotting the state as a function of physical time we make the "slow" transition "infinitely fast" and recover the discontinuity in time. More formally, this is encapsulated in (19), which shows that the integral of a Dirac impulse $\lambda_n(t)$ is the same as the integral of the $v$-state of the auxiliary ODE over a finite time interval $[\tau_s, \tau_r]$ of nonzero length. To simulate a time-freezing system with Zeno's effect, the numerical time horizon would have to be infinitely long, as every state jump requires $\tau_{jump} > 0$. In this paper, we assume to have a finite number of impacts. In practical robotics applications, one is usually interested in solutions with a finite number of impacts.

## 3.6 Possible extensions

In this paper, we consider a single unilateral constraint. To extend the same ideas for multiple and/or simultaneous impacts one must take care of the chosen impact model [7]. The extension can be made in several ways, e.g., time continues to flow when the first normal velocity component reaches zero or when all of them reach zero. Some multiple impact models suffer from nonuniqueness of solutions [7,26,35], and how to proceed is a modeling decision. In the special case when the constraints, e.g., $f_{c,i}(q)$ and $f_{c,j}(q)$ are orthogonal in the kinetic metric, i.e., $f_{c,i}(q)^{\top}M(q)^{-1}f_{c,j}(q) = 0$, the impacts can be
treated independently [7]. In this case, one would take for every constraint the auxiliary dynamics from Prop. 3 and the time would continue to flow when all normal velocity components reach zero. Stochastic approaches with multiple outcomes are also possible, e.g., within the binary collision law as shown in [26] or a stochastic version of Routh's impact model [18].

# 4 Frictional impact

If friction is present at the contact point, frictional impulses cause state jumps in the tangential directions. This section extends the time-freezing reformulation for CLS with frictional impacts. Appropriate auxiliary dynamics for the state jumps in the tangential directions are introduced. The time-freezing system covers both stick and slip motions.

## 4.1 CLS with frictional impacts

The extension of the CLS (3) with Coulomb friction reads as:

$$ \dot{q} = v, \qquad (20a) $$

$$ \dot{v} = f_v(q, v, u) + M(q)^{-1}(n(q)\lambda_n + B(q)\lambda_t), \quad (20b) $$

$$ 0 \le \lambda_n \perp f_c(q) \ge 0, \qquad (20c) $$

$$ 0 = n(q(t_s))^\top v(t_s^+), \qquad (20d) $$

$$ \text{if } f_c(q(t_s)) = 0 \text{ and } n(q(t_s))^\top v(t_s^-) < 0, $$

$$ \lambda_t \in \arg \min_{\tilde{\lambda}_t \in \mathbb{R}^{n_t}} -v^\top B(q) \tilde{\lambda}_t \qquad (20e) $$

$$ \text{s.t.} \quad \|\tilde{\lambda}_t\|_2 \le \mu \lambda_n. \qquad (20f) $$

Compared to the CLS (3), the model is extended by the term $M(q)^{-1}B(q)\lambda_t$ in the r.h.s. of velocity dynamics (friction force) and Eq. (20e)-(20f) (friction model). The matrix $B(q) = [b_1(q), \dots, b_{n_t}(q)] \in \mathbb{R}^{n_q \times n_t}$ spans the tangent space at contact points $\{q \in \mathbb{R}^{n_q} | f_c(q) = 0\}$ and $\lambda_t \in \mathbb{R}^{n_t}$ is the friction force. For planar contacts, we have $n_t = 1$, and for 3D contacts $n_t = 2$. Denote the tangential velocity during contact by $v_t := B(q)^\top v \in \mathbb{R}^{n_t}$. We assume that the vectors $b_1(q), b_2(q)$ and $n(q)$ are orthogonal in the kinetic metric, i.e., $b_1(q)^\top M(q)^{-1}n(q) = 0$, $b_2(q)^\top M(q)^{-1}n(q) = 0$ and $b_1(q)^\top M(q)^{-1}b_2(q) = 0$. While restrictive, this assumption simplifies computations and maintains brevity in the exposition. Following a similar approach without it leads to the same conclusions but with more complex equations. The convex optimization problem (20e)-(20f) is the maximum dissipation principle [25] and $\mu > 0$ is the coefficient of friction. This model expresses that the dissipation of the kinetic energy between two objects in contact is maximized. When impacts occur, the impulsive $\lambda_n$ results in an impulsive $\lambda_t$ via Eq. (20e)-(20f) and thus we have state jumps in the velocity $v$ in the tangential directions $b_j(q), j = 1, \dots, n_t$.

## 4.2 Stick and slip dynamics of the CLS

For a given $\lambda_n$ the solution map of the convex optimization problem (20e)-(20f) is given by [25]:

$$ \lambda_t \in \begin{cases} \{-\mu\lambda_n \frac{v_t}{\|v_t\|_2}\}, & \text{if } \|v_t\|_2 > 0, \\ \{\tilde{\lambda}_t | \|\tilde{\lambda}_t\|_2 \le \mu\lambda_n\}, & \text{if } \|v_t\|_2 = 0. \end{cases} \qquad (21) $$

When the system is in contact, it can be in slipping motion, i.e., it has nonzero tangential velocity $v_t \ne 0$, or in sticking motion with $v_t = 0$. Similar to Section 3.1, we derive equivalent ODEs which model the stick and slip dynamics during contact phases. If the system is in slip motion and $\lambda_n > 0$, it follows from (21) that $\lambda_t = -\mu\lambda_n \frac{v_t}{\|v_t\|_2}$. Similar to Eq. (13) we obtain

$$ \dot{q} = v, \qquad (22a) $$

$$ \dot{v} = f_v(q, v, u) - M(q)^{-1}D(q)^{-1}\varphi(x, u) \qquad (22b) $$

$$ (n(q) - B(q)\mu \frac{v_t}{\|v_t\|_2}). $$

The r.h.s. of this ODE is compactly denoted by $f_{\text{Slip}}(x, u)$. If the system is in sticking motion we have $v_t = B(q)^\top v = 0$ and $n^\top v = 0$. By differentiating these equations w.r.t. time we can explicitly compute the multipliers $\lambda := (\lambda_n, \lambda_t)$ by

$$ \lambda = -\tilde{D}(q)^{-1}\tilde{\varphi}(x, u), $$

$$ \tilde{D}(q) = [n(q) B(q)]^\top M(q)^{-1} [n(q) B(q)], $$

$$ \tilde{\varphi}(x, u) := [n(q) B(q)]^\top f_v(q, v, u) + \nabla_q([n(q) B(q)]^\top v)^\top v. $$

The ODE describing the sticking motion reads as

$$ \dot{q} = v, \qquad (23a) $$

$$ \dot{v} = f_v(q, v, u) - M(q)^{-1} [n(q) B(q)] \tilde{D}(q)^{-1} \tilde{\varphi}(x, u). \qquad (23b) $$

Its r.h.s. is compactly denoted by $f_{\text{Stick}}(x, u)$. The transition from the stick to the slip mode occurs when the other tangential forces are greater than the maximal friction force $\lambda_t$, cf. [7, Chapter 5].

## 4.3 Time-freezing for frictional impacts in the 2D case

In the planar case we have $n_t = 1$ and Eq. (21) simplifies to $\lambda_t \in -\mu\lambda_n\text{sign}(v_t)$. Denote the single column of $B(q)$ by $b(q)$. Depending on the sign of $v_t$, we define an auxiliary dynamical system to mimic the state jump in the tangential direction $b(q)$. For the $n(q)$-direction we use the dynamics from Proposition 3. For the tangential direction and $b(q)^\top v < 0$ we define the tangential auxil-
iary dynamics analogously:

$$y' = f_{\text{aux,t}}^{-}(y) := \begin{bmatrix} \mathbf{0}_{n_q,1} \\ M(q)^{-1}b(q)a_t \\ 0 \end{bmatrix}. \quad (24)$$

To account for the sign of the tangential velocity, for $b(q)^{\top}v > 0$ we use $y' = f_{\text{aux,t}}^{+}(y) := -f_{\text{aux,t}}^{-}(y)$. Depending on the sign of $v_t$, one of these ODE is active for the same numerical time interval of the length $\tau_{\text{jump}}$ as $y' = f_{\text{aux,n}}(y)$. Furthermore, we know from Eq. (19) in Theorem 5 that the impulse bringing the normal velocity $n(q)^{\top}v < 0$ to zero after an impact is proportional to $a_n\tau_{\text{jump}}$. Thus, by settings $a_t = \mu a_n$, we conclude that the integrals of the auxiliary dynamics satisfy the maximum dissipation principle, i.e., $a_t\tau_{\text{jump}} = -\mu a_n\tau_{\text{jump}}\text{sign}(b(q)^{\top}v)$.

State jumps in both the normal and tangential directions are treated simultaneously with different auxiliary dynamics. They should be active whenever $y \in Q := \{y \in \mathbb{R}^{n_y} | c_1(y) < 0, c_2(y) < 0\}$, cf. Fig. 1. To treat different signs of the tangential velocity we introduce the switching function $c_3(y) = b(q)^{\top}v$. Hence, we have in total $n_f = 3$ regions, one for the unconstrained dynamics and two to mimic the state jumps. We extend the definition of the regions in Eq. (14) as follows:

$$R_1 = \{y \in \mathbb{R}^{n_y} | c_1(y) > 0\} \cup \{y \in \mathbb{R}^{n_y} | c_1(y) < 0, c_2(y) > 0\}, \quad (25a)$$

$$R_2 = Q \cap \{y \in \mathbb{R}^{n_y} | c_3(y) > 0\}, \quad (25b)$$

$$R_3 = Q \cap \{y \in \mathbb{R}^{n_y} | c_3(y) < 0\}. \quad (25c)$$

The sum of the corresponding auxiliary dynamics accounts for the simultaneous state jumps, i.e., $f_2(y) = f_{\text{aux,n}}(y) + f_{\text{aux,t}}^{+}(y)$, $f_3(y) = f_{\text{aux,n}}(y) + f_{\text{aux,t}}^{-}(y)$. The time-freezing system for the CLS (20) is given in the next definition.

**Definition 6 (Time-freezing system with friction)**
Let $\tau \in \mathbb{R}$ be the numerical time, $y(\tau) := (x(\tau), t(\tau)) \in \mathbb{R}^{n_y}$ the differential states and $u(\tau) \in \mathbb{R}^{n_u}$ a given control function. The time-freezing PSS is defined by the regions $R_i$, $i = 1, \dots, 3$, in Eq. (25) with

$$\begin{aligned} f_1(y, u) &= (f_{\text{ODE}}(x, u), 1), \\ f_2(y) &= f_{\text{aux,n}}(y) + f_{\text{aux,t}}^{+}(y), \\ f_3(y) &= f_{\text{aux,n}}(y) + f_{\text{aux,t}}^{-}(y). \end{aligned}$$

The corresponding Filippov system, which we call the time-freezing system, is denoted by $y' \in F_{\text{TF}}(y, u)$. The set $F_{\text{TF}}(y, u)$ is defined as in (2). It is assumed that appropriate dynamics $f_{\text{aux,n}}(y)$, $f_{\text{aux,t}}^{+}(y)$ and $f_{\text{aux,t}}^{-}(y)$ exist.

As in the frictionless case, we are interested in the relation of the CLS in contact mode and the corresponding

sliding mode of the time-freezing system on $\Sigma$.

**Theorem 7 (Slip-stick sliding mode)** Suppose that the auxiliary dynamics from Proposition 3 and Eq. (24) are used in the time-freezing system from Definition 6. Let $y(\tau)$ be a solution of this system with $y(0) \in \Sigma$ and $\tau \in [0, \tau_f]$. Suppose that $\varphi(x(\tau), u(\tau)) \le 0$ for all $\tau \in [0, \tau_f]$ (persistent contact), then the following statements are true:

(i) If $v_t \ne 0$ (slip motion), then the sliding mode dynamics are given by $y' = \gamma(x, u)(f_{\text{Slip}}(x, u), 1)$.

(ii) If $v_t = 0$ (stick motion), then the sliding mode dynamics are given by $y' = \gamma(x, u)(f_{\text{Stick}}(x, u), 1)$,

where $\gamma(x, u) \in (0, 1]$ is a time-rescaling factor defined in Eq. (17).

**PROOF.** See Appendix B. $\square$

This result generalizes Theorem 4 and one can see that the sliding mode dynamics match the slip or stick dynamics of the CLS. Finally, we show that CLS with planar contacts, friction, and impacts are equivalent to Filippov systems.

**Theorem 8 (Solution relationship)**

Regard the initial value problems corresponding to: i) the time-freezing system in Definition 6 with a given $y(0) = (q_0, v_0, 0) \in \mathbb{R}^{n_y}$ and $f_c(q_0) \ge 0$ on a time interval $[0, \tau_f]$, ii) the CLS from Eq. (20) with the initial value $x(0) = (q_0, v_0) \in \mathbb{R}^{n_x}$ on a time interval $[0, t_f] := [0, t(\tau_f)]$, with $f_c(q(t_f)) \ge 0$ and $n(q(t_f))^{\top}v(t_f) \ge 0$. Suppose the following assumptions hold:

(a) the auxiliary dynamics $f_{\text{aux,n}}(y)$ from Proposition 3 and $f_{\text{aux,t}}^{-}(y)$, $f_{\text{aux,t}}^{+}(y)$ from Eq. (24) are used in the time-freezing system in Definition 6,

(b) there is at most one time point $t_s = t(\tau_s)$ where $f_c(q(t_s)) = 0$ and $n(q(t_s))^{\top}v(t_s^-) < 0$ on the time interval $[0, t_f]$,

Then, the solutions to the two problems are related as follows:

(1) For $t \ne t_s$:

$$\begin{aligned} x(t(\tau)) &= Ry(\tau), \quad \text{with } R = [I_{n_x} \quad \mathbf{0}_{n_x,1}], \\ \lambda(t(\tau)) &= \begin{cases} \lambda_{\text{Slip}}(t(\tau)), & \text{if } y \in \Sigma, v_t \ne 0, \\ \lambda_{\text{Stick}}(t(\tau)), & \text{if } y \in \Sigma, v_t = 0, \\ 0, & \text{otherwise.} \end{cases} \end{aligned}$$

with

$$\lambda_{\text{Slip}} = -D(q)^{-1}\varphi(x, u)(1, -\mu\text{sign}(v_t))$$

$$\lambda_{\text{Stick}} = -\tilde{D}(q)^{-1}\tilde{\varphi}(x, u).$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. The exact friction force Eq. (21) (left) and its approximation in Eq. (26) (right).</figcaption></figure>

(2) For $t = t_s$:

$$ \lim_{\epsilon \to 0} \int_{t_s - \epsilon}^{t_s + \epsilon} \lambda_n(t) dt = \int_{\tau_s}^{\tau_r} a_n d\tau, \\ \lim_{\epsilon \to 0} \int_{t_s - \epsilon}^{t_s + \epsilon} |\lambda_n(t)| dt = \int_{\tau_s}^{\tau_r} \mu a_n d\tau, $$

**PROOF.** Theorem 7 is applied for $y \in \Sigma$, otherwise the proof follows similar lines as the proof of Theorem 5. The absolute value in the last equation accounts for all signs of $v_t$. $\square$

## 4.4 Time-freezing for frictional impacts in the 3D case

This case is more difficult since we cannot easily treat different directions of $v_t$ with different auxiliary dynamics as in the planar case. The solution map of (21) depends discontinuously on $\|v_t\|_2$, hence, we must take it as a switching function. The set $\|v_t\|_2 = 0$ has no interior and we cannot use the Filippov extension from Eq. (2), which assumed regions $R_i$ with nonempty interior. More general definitions without multipliers $\theta$ that could treat this case exist [14], but they are not computationally useful in our case, as we see in the next section. To alleviate this difficulty we propose an approximation for (21):

$$ \lambda_t = \begin{cases} -\mu\lambda_n \frac{v_t}{\|v_t\|_2}, & \text{if } \|v_t\|_2 > \epsilon_t, \\ v_t, & \text{if } \|v_t\|_2 < \epsilon_t, \end{cases} \quad (26) $$

with a small parameter $\epsilon_t > 0$. This expression is exact for $\|v_t\|_2 > \epsilon_t$, thus we can make it arbitrarily accurate. For $\|v_t\|_2 < \epsilon_t$, the vector field drives the tangential velocity towards $\|v_t\|_2 = \epsilon_t$, see Fig. 4. In a Filippov setting, a convex combination of the two cases in (26) keeps the velocity at $\|v_t\|_2 = \epsilon_t$. Hence, in sticking mode, we have a velocity drift of $\epsilon_t$. Now by taking $c_3(y) = \|v_t\|_2 - \epsilon_t$ as a switching function, we can define regions with nonempty interiors and the corresponding auxiliary dynamics. The auxiliary dynamics mimicking

<figure><img src="image_7.png" /><figcaption>Fig. 5. Trajectories of the time-freezing system from Ex. 4.</figcaption></figure>

the behavior of (26) read as

$$ \begin{aligned} f_{\text{aux},t}^{+}(y) &= \begin{bmatrix} \mathbf{0}_{n_q,1} \\ -M(q)^{-1}B(q)a_t \frac{v_t}{\|v_t\|} \\ 0 \end{bmatrix}, \\ f_{\text{aux},t}^{-}(y) &= \begin{bmatrix} \mathbf{0}_{n_q,1} \\ M(q)^{-1}B(q)v_t \\ 0 \end{bmatrix}. \end{aligned} $$

The regions for the time-freezing system are defined as in (25) and the matching time-freezing system is defined analogously to Definition 6. Furthermore, one could derive stick-slip dynamics corresponding to the solution map approximation (26) and relate it to the time-freezing system by following similar lines as in Theorems 7 and 8, but we omit the details here. We conclude this section by revisiting Example 3, but now with adding friction.

**Example 4 (Frictional impact)** The time-freezing system model from Example 2 is extended by adding friction with a coefficient $\mu = 0.6$. In the planar case, we have $n_t = 1$ and the tangent at the contact point is $b(q) = [1, 0]$. We have the switching functions $c(y) = (q_2, v_2, v_1)$. Following Eq. (25), the regions of the time-freezing system are $R_1 = \{y \mid q_2 > 0\} \cup \{y \mid q_2 < 0, v_2 > 0\}$, $R_2 = \{y \mid q_2 < 0, v_2 < 0, v_1 > 0\}$ and $R_3 = \{y \mid q_2 < 0, v_2 < 0, v_1 < 0\}$. The dynamics of the PSS are $f_1 = (v_1, v_2, u_1, -g + u_2, 1)$, $f_2 = (0, 0, -\mu a_n, a_n, 0)$ and $f_3 = (0, 0, \mu a_n, a_n, 0)$. The results of the simulation are depicted in Fig. 5. Note that due to friction there is now also a state jump in the tangential velocity $v_1$, cf. middle plots. Afterward, the acceleration of $v_1$ is during contact phases smaller due to the friction force. However, the tangential acceleration is increasing over time as the normal contact force becomes weaker because of $u_2$. At $\tau = 2.8$ the particle lifts off as in the previous example.
# 5 The application of time-freezing in optimal control problems

This section regards OCP formulations with time-freezing systems and numerical methods to solve them. We derive an equivalent OCP, now subject to the time-freezing system. The section concludes with discussing the numerical methods and software for solving OCPs subject to time-freezing systems.

## 5.1 Continuous-time OCP with a CLS

We regard a modification of the OCP (4), where we consider the CLS with friction (20). This continuous-time OCP read as:

$$ \min_{x(\cdot), \lambda(\cdot), u(\cdot)} \Psi(x(T)) \quad (27a) $$

$$ \text{s.t.} \quad x(0) = \bar{x}_0, \quad (27b) $$

$$ \text{Eq. (20), } t \in [0, T] \quad (27c) $$

$$ 0 \le g(x(t), u(t)), \quad t \in [0, T], \quad (27d) $$

$$ 0 \le r(x(T)), \quad (27e) $$

The functions $g: \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \to \mathbb{R}^{n_g}$ and $r: \mathbb{R}^{n_x} \to \mathbb{R}^{n_r}$ are the path and terminal constraints, respectively. W.l.o.g. we only consider a terminal cost term here and remind the reader that the integral of a running cost $L: \mathbb{R}^{n_x} \times \mathbb{R}^{n_u} \to \mathbb{R}$ over $[0, T]$ can be treated via a terminal cost term by introducing a quadrature state $\ell(t)$

$$ \frac{d}{dt} \ell(t) = L(x(t), u(t)), \quad t \in [0, T], \quad \ell(0) = 0, \quad (28) $$

and adding $\ell(T)$ to the objective.

## 5.2 Continuous-time OCP with a time-freezing system

Based on the results from the previous sections, we derive now an OCP subject to a time-freezing system. The new OCP is regarded in numerical time $\tau \in [0, \tilde{T}]$. We take four steps in this transformation: (1) we modify the quadrature state in Eq. (28) so that the cost integrated over both numerical and physical time remains unchanged; (2) we reformulate the time-freezing system into an equivalent dynamic complementarity system to make it possible to apply FESD; (3) we introduce a time-transformation to ensure that the terminal physical time $t(\tilde{T})$ matches the true control horizon of Eq.(27), i.e., $t(\tilde{T}) = T$; (4) we express the remaining constraints in terms of numerical time.

We start with adapting the objective. This is achieved

by replacing the quadrature state (28) by:

$$ \frac{d}{d\tau} \ell(\tau) = \begin{cases} L(x(\tau), u(\tau)), & \text{if } y \in R_1, \\ 0, & \text{otherwise.} \end{cases} \quad (29) $$

When the time is frozen the cost integral is zero and there are no contributions to the overall objective, i.e., the cost is unchanged when the time is frozen.

Next, we rewrite the time-freezing system from Definition 6 as an equivalent dynamic complementarity system. This enables the application of the FESD method to the time-freezing system. To achieve this, we use the set-valued step function $\alpha(x) = 0.5(1 + \text{Sign}(x))$. It is well-known that Filippov multipliers $\theta$ can be expressed as products of set-valued step functions $\alpha_i$ [3,13]. We define $c(y) := (c_1(y), c_2(y), c_3(y))$. The set-valued step function $\alpha(\cdot)$ can be expressed as the solution map of a parametric linear program [3]:

$$ \alpha(c(y)) \in \arg \min_{\tilde{\alpha} \in \mathbb{R}^3} -c(y)\tilde{\alpha} \quad \text{s.t.} \ 0 \le \tilde{\alpha} \le e. \quad (30) $$

Using the Karush-Kuhn-Tucker (KKT) conditions of (30) we can derive from $y' \in F_{TF}(y, u)$ in Definition 6 the equivalent dynamic complementarity system:

$$ y' = F(y, u) \theta, \quad (31a) $$

$$ 0 = g_F(\theta, \alpha), \quad (31b) $$

$$ 0 = c(y) - \lambda^P + \lambda^n, \quad (31c) $$

$$ 0 \le \alpha \perp \lambda^n \ge 0, \quad (31d) $$

$$ 0 \le e - \alpha \perp \lambda^P \ge 0. \quad (31e) $$

The matrix $F(y, u) = [f_1(y, u), \dots, f_{n_f}(y)] \in \mathbb{R}^{n_y \times n_f}$ collects the modes of the PSS and $\theta = (\theta_1, \dots, \theta_{n_f})$. The last three lines are the KKT conditions of (30), where $\lambda^n, \lambda^P \in \mathbb{R}^3$ are the Lagrange multipliers for the lower and upper bounds in (30), respectively. We group all algebraic variables of the DCS in the vector $z = (\theta, \alpha, \lambda^P, \lambda^n)$. The function $g_F$ relates the Filippov multipliers $\theta$ with the evaluations of the step functions $\alpha$:

$$ g_F(\theta, \alpha) := \begin{bmatrix} \theta_1 - \alpha_1 + (1 - \alpha_1)\alpha_2, \\ \theta_2 - (1 - \alpha_1)(1 - \alpha_2)(1 - \alpha_3) \\ \theta_3 - (1 - \alpha_1)(1 - \alpha_2)(\alpha_3) \end{bmatrix}. \quad (32) $$

These expressions correspond to the signs of $c_j(y)$ in the definitions of $R_i$, e.g., $c_1(y) > 0, c_2(y) < 0$ results in $\alpha_1(1 - \alpha_2)$ [3,13].

The time-freezing system evolves over $\tau \in [0, \tilde{T}]$. During state jumps the physical time evolution is stopped. As a consequence, we have that $t(\tilde{T}) < T$, i.e., the terminal physical time in the time-freezing problem does not match the desired time $T$. To resolve this, we introduce a time-transformation variable $s(\tau) \in \mathbb{R}$ and impose the terminal constraint on the clock state $t(\tilde{T}) = T$. The
function $s(\cdot)$ can be interpreted as virtual control that controls the physical time in the numerical time. Consequently, we obtain $t' = s$ and with $s > 1$, which speeds up the physical time and allows us to catch up and reach the desired time $T$, cf. Example 5. Such time transformations are very common in optimal control when one wants to optimize over the terminal time [33].

It is left to impose the path (27d) and terminal constraints (27e) in numerical time for $x(\tau)$ and $u(\tau)$. Finally, the OCP subject to the time-freezing system reads as:

$$ \min_{\substack{y(\cdot), z(\cdot), \\ u(\cdot), s(\cdot)}} \Psi(x(\tilde{T})) \qquad (33a) $$

$$ \text{s.t.} \quad x(0) = \bar{x}_0, \quad t(0) = 0, \qquad (33b) $$

$$ y'(\tau) = s(\tau)F(y(\tau), u(\tau))\theta(\tau), \tau \in [0, \tilde{T}], \qquad (33c) $$

$$ 0 = g_F(\theta(\tau), \alpha(\tau)), \tau \in [0, \tilde{T}], \qquad (33d) $$

$$ 0 = c(y(\tau)) - \lambda^p(\tau) + \lambda^n(\tau), \tau \in [0, \tilde{T}], \qquad (33e) $$

$$ 0 \le \alpha(\tau) \perp \lambda^n(\tau) \ge 0, \tau \in [0, \tilde{T}], \qquad (33f) $$

$$ 0 \le e - \alpha(\tau) \perp \lambda^p(\tau) \ge 0, \tau \in [0, \tilde{T}], \qquad (33g) $$

$$ 0 \le g(x(\tau), u(\tau)), \quad \tau \in [0, \tilde{T}], \qquad (33h) $$

$$ 0 \le r(x(\tilde{T})), \qquad (33i) $$

$$ t(\tilde{T}) = T. \qquad (33j) $$

It is important to note, that when the time is frozen ($t' = 0$) the control $u(\tau)$ does not influence $x(\tau)$, since the auxiliary dynamics does not depend on the control, cf. Eq (16). Therefore, one could even omit the path constraints whenever $t' = 0$, but we keep it for notational simplicity. Additionally, the integral of the stage cost remains unchanged, since $\frac{d}{d\tau}\ell(\tau) = 0$ in this case, cf. (29).

Next, we show that the optimal controls obtained by solving the initial OCP (27), with appropriate modifications, are also optimal for (33). Let $u^*(t), t \in [0, T]$ be an optimal control of (27). We construct an $\tilde{u}^*(\tau)$, $\tau \in [0, \tilde{T}]$ as follows. It can be seen that when the physical time is evolving ($t' > 0$), we can find the inverse function $t^{-1}(\cdot)$ to find the corresponding numerical time $\tau$. We construct a control function for the time-freezing system:

$$ \tilde{u}^*(\tau) = \begin{cases} u(t^{-1}(t(\tau))), & \text{for } t(\tau)' > 0 \\ \hat{u}(\tau), & \text{for } t(\tau)' = 0, \end{cases} \qquad (34) $$

where $\hat{u}(\tau)$ is any function such that $g(x(\tau), \hat{u}(\tau)) \ge 0$ holds, whenever $t'(\tau) = 0$. Recall that $\hat{u}(\tau)$ does not change the objective nor it changes $x(\tau)$, its only purpose is to extend $u(t)$ to intervals when the time is frozen. For example, we can choose a constant value that does not violate the path constraints. With (34) we can extend $u^*(t)$ for the time interval where the physical time is frozen. Conversely, given an optimal control $\tilde{u}^*(\tau)$ of Eq. (33), then we expect $u^*(t(\tau))$ to be optimal for (27).

**Theorem 9** Let $\tilde{u}^*(\tau)$, $\tau \in [0, \tilde{T}]$ be an optimal control obtained by solving the OCP (33). Then $u^*(t) = \tilde{u}^*(t(\tau))$, $t \in [0, T]$ is an optimal control of the OCP (27). Conversely, let $u^*(t)$, $t \in [0, T]$ be an optimal control of the OCP (27), then the control function $\tilde{u}^*(\tau)$, $\tau \in [0, \tilde{T}]$ obtained via Eq. (34) is optimal for (33).

**PROOF.** For a fixed control function $u(\tau)$ and $s(\tau)$ such that $t(\tilde{T}) = T$, the time-freezing system (33c)-(33g) and the CLS (20) with $u(t(\tau))$ are equivalent in the sense of Theorem 8. Thus, a feasible $y(\tau)$ in (33) results in a feasible $x(t)$ in (27). Due to equation (29), both OCPs have the same objective value. Consequently, given a $\tilde{u}(\tau) = \tilde{u}^*(\tau) + \delta\tilde{u}(\tau)$ that improves the objective (33a), the corresponding $\tilde{u}(t(\tau))$ would also improve (27a). Conversely, for every modified $u(t(\tau)) = u^*(t(\tau)) + \delta u(t(\tau))$ that improves the objective (27a), we can construct an appropriate control function $u(\tau)$ via (34) that improves (33a). Thus, $u^*(t) = \tilde{u}^*(t(\tau))$ is optimal for (27). The converse is proved by similar arguments. $\square$

### 5.3 Discrete-time OCP with the time-freezing system

In principle, one can discretize the OCP (27) by using any time-stepping integration method for CLS [2,7,35] e.g., the Stewart-Trinkle method [37]. Such an approach for direct optimal control was used in [32]. As discussed in Section 1, standard time-stepping methods for CLS with friction (20) have at best first-order accuracy [2,35]. Moreover, the numerical sensitivities obtained from such a discretization are always wrong and the NLP solvers converge to spurious solutions [36,41]. Therefore, for a moderately accurate solution usually a large computational effort is needed. Standard time-stepping methods for PSS encounter the same difficulties as methods for CLS. These fundamental limitations motivate the derivation of the OCP formulation with a time-freezing system (33) since we can use the recently introduced FESD method that overcomes these difficulties [30]. It discretizes dynamic complementarity systems that are equivalent to PSS (such as Eq. (31)) and ensures automatic switch detection, higher-order integration accuracy, and correct numerical sensitivities [28,30]. In conclusion, the fundamental limitations of standard direct optimal control methods are resolved by combining time-freezing and FESD. This enables one to find a more accurate solution approximation for the continuous-time OCP (27) by solving (33).

We proceed by introducing the discrete-time version of the OCP (33) with a multiple shooting-type discretization [6]. The numerical time horizon $[0, \tilde{T}]$ is split into $N$ control intervals $[\tau_k, \tau_{k+1}]$ of equal length [30]. The controls are assumed to be constant over every interval, i.e., $u(\tau) = u_k, \tau \in [\tau_k, \tau_{k+1}]$, $k = 0, \dots, N-1$, and $y_k = (x_k, t_k) \in \mathbb{R}^{n_y}$ is the discrete-time approximation of the time-freezing state, i.e., $x_k \approx x(\tau_k) \quad t_k \approx t(\tau_k)$.
The vectors $z_k$ collect all algebraic and internal integration variables for the $k$-th control interval. The vector $w := (y_0, z_0, u_0, s_0, \dots, y_{N-1}, z_{N-1}, u_{N-1}, s_{N-1}, y_N)$ groups all optimization variables.

Our goal is to have an equidistant control grid, as this is typically required in feedback control applications. It is important to note that, due to intervals with frozen physical time evolution ($t' = 0$), an equidistant grid in numerical time $\{\tau_0, \dots, \tau_N\}$ does not imply an equidistant grid in physical time $\{t_0, \dots, t_N\}$. To address this issue, we propose to use a piecewise constant discretization of the time-transformation variable $s(\tau)$, i.e., we have $s_k \in \mathbb{R}$, $k = 0, \dots, N-1$. Additionally, we introduce the constraints $t_k = k \frac{T}{N}$, $k = 0, \dots, N$, cf. Eq. (35f) below. It is worth noting that for $k = N$, we have the discrete-time versions of the terminal clock constraint (33j). The steps above result in an equidistant control discretization grid in physical time, i.e., $u(t) = u_k$ for $t \in [t_k, t_{k+1}]$ with $t_0 = 0$ and $t_k = t_{k-1} + T/N$. This is further illustrated in Example 5.

The discretization of (33) reads as:

$$ \min_{w} \Psi(x_N) \qquad (35a) $$

$$ \text{s.t.} \quad x_0 = \bar{x}_0, \qquad (35b) $$

$$ y_{k+1} = \Phi_f(y_k, z_k, u_k, s_k), \quad k = 0, \dots, N-1, \qquad (35c) $$

$$ 0 = \Phi_{\text{int}}(y_k, z_k, u_k), \quad k = 0, \dots, N-1, \qquad (35d) $$

$$ 0 \le \Phi_{c,1}(z_k) \perp \Phi_{c,2}(z_k) \ge 0, \quad k = 0, \dots, N-1, \qquad (35e) $$

$$ t_k = k \frac{T}{N}, \quad k = 0, \dots, N, \qquad (35f) $$

$$ 1 \le s_k \le \bar{s}, \quad k = 0, \dots, N-1, \qquad (35g) $$

$$ 0 \le g(x_k, u_k), \quad k = 0, \dots, N-1, \qquad (35h) $$

$$ 0 \le r(x_N). \qquad (35i) $$

It is common in direct optimal control to write discretization method equations in a compact discrete-time system manner [33, Chapter 8], as we do here in Eq. (35c)-(35e). The function $\Phi_f: \mathbb{R}^{n_y} \times \mathbb{R}^{n_z} \times \mathbb{R}^{n_u} \times \mathbb{R} \to \mathbb{R}^{n_x}$ is the discrete-time state transition map which approximates $y(\tau)$. The function $\Phi_{\text{int}}: \mathbb{R}^{n_y} \times \mathbb{R}^{n_z} \times \mathbb{R}^{n_u} \to \mathbb{R}^{n_{\Phi}}$ collects all internal computations of the underlying integration scheme. The constraints (35e) arise from the discretization of the complementarity conditions (33f)-(33g). These functions are obtained via the FESD discretization, for more details see [30, Section 3.2.4]. The constraint (35g) bounds $s_k$, where $\bar{s}$ is its maximal value that has to be sufficiently large to ensure feasibility of (35f).

The FESD method for PSS is implemented in the CasADi [4] based open-source tool NOSNOC [1, 28]. Note that the NLP (35) is a mathematical program with complementarity constraints. They are degenerate nonsmooth NLP which are solved in NOSNOC with a homotopy approach, cf. [28, Section IV.B]. The advantage of the homotopy approach is that only a (finite) sequence

<figure><img src="image_8.png" /><figcaption>Fig. 6. Solution to the guiding optimal control example.</figcaption></figure>

of related, but smooth NLP is solved. Under some regularity assumptions, the solution of the last NLP is a solution of the initial nonsmooth NLP [5, 17]. The main drawback of the homotopy approach is that it sometimes requires some tuning of the homotopy parameters.

**Example 5 (OCP Example)** We solve an OCP of the form of (27) with our guiding example. The initial value is unchanged, i.e., $y_0 = (0, 1, 0, 0, 0)$. The particle should reach at $t(\bar{T}) = 2$, with $\bar{T} = 2$, the position $q(T) = (3, 0)$ with zero terminal velocity $v(T) = (0, 0)$. We bound the horizontal thrust force $|u_1| \le 10$ and set for simplicity $u_2 = 0$. The ball should reach the goal with minimum control effort, which is modeled with the stage cost $L(x, u) = u_1^2$. We take $N = 20$ control intervals and discretize the equivalent time-freezing OCP with a third-order FESD-Radau II-A scheme with three integration steps on every control interval [28]. The solution is depicted in Fig. 6. We can see that maximum force is applied before the impact since there is still no friction and the motion is cheaper. After the impact, a smaller control force is applied just to reach the target. Note that the $s_k$ (yellow line in the bottom right plot) is higher during the control interval when the state jump happens (to catch up the frozen time) and $s_k = 2$ during contact to compensate for the slow down due to $\gamma(x, u)$. The resulting speed of time is always one (bottom left plot), except when the state jump happens where a speed-up is needed to compensate for the frozen time. This ensures an equidistant control grid (as intended with the constraint (35f)) and $t(\bar{T}) = 2$ as desired.

# 6 Numerical optimal control of a jumping robot

We consider a hopping robot that must jump over three holes to reach a desired target. This example showcases all theoretical developments of Sections 3, 4, and the use of numerical methods described in Section 5. Thereby, an OCP formulation for synthesizing dynamic motions of the single-legged 2D robot Capler [9, 21] is derived. The robot is described by four degrees of freedom $q = (q_x, q_z, \phi_{\text{knee}}, \phi_{\text{hip}})$. Here, $(q_x, q_z)$ are the coordinates of
the robot's base at the hip and $\phi_{\text{knee}}$, $\phi_{\text{hip}}$ are the angles of the hip and knee, respectively, cf. left plot in Fig. 7. It is actuated by two direct-drive motors at the hip and knee joints. The robot's dynamics are compactly described by the CLS in the form of (20). The torques of the two motors $u(t) = (u_{\text{knee}}(t), u_{\text{hip}}(t))$ are the control variables. A detailed derivation of the model equations and all parameters for the robot can be found in [15, Appendix A].

Denote by $p_{\text{foot}}(q) = (p_{\text{foot},x}(q), p_{\text{foot},z}(q))$ and $p_{\text{knee}}(q) = (p_{\text{knee},x}(q), p_{\text{knee},z}(q))$ the kinematic position of the robot's foot and knee, respectively. For the unilateral constraint function we take $f_c(q) = p_{\text{foot},z}(q)$. For a planar robot, we need just one tangent, i.e., $b(q) = \nabla_q p_{\text{foot},x}(q)$ and the friction model is exact. The coefficient of friction is $\mu = 0.8$ and the auxiliary ODE constant is $a_n = 200$.

The objective of the OCP is to minimize the integral of the squared control torques, i.e., we have the stage cost $L(x, u) = u(\tau)^{\top}u(\tau)$. The robots should reach a given target position $q_{\text{target}} = (3, 0.4, 0, 0)$ starting from the initial position $q_0 = (0, 0.4, 0, 0)$ with zero velocity $v_0 = \mathbf{0}_{4,1}$. The initial value is $y_0 = (q_0, v_0, 0)$.

The prediction horizon is $T = 2.5$ s. We add the following constraints on the states and kinematic positions:

$$ 
\begin{aligned}
-0.05e &\le (q_x(t), p_{\text{foot},x}(q(t)), p_{\text{knee},x}(q(t))), \\
0.2 &\le q_z(t) \le 0.55, \\
-\frac{3\pi}{8} &\le \phi_{\text{hip}}(t) \le \frac{3\pi}{8}, \\
-\frac{\pi}{2} &\le \phi_{\text{knee}}(t) \le \frac{\pi}{2}, \\
0.05 &\le p_{\text{knee},z}(q(t)), \\
-0.005 &\le p_{\text{foot},z}(q(t)) \le 0.2, \quad t \in [0, T].
\end{aligned}
$$

Their goal is twofold. On one hand, they should avoid unnatural and too extensive bending of the joints. On the other hand, they serve as guiding constraints during the early phases of the homotopy procedure. In the early iterations, the physics is relaxed and we want to prohibit the optimizer to go to undesired regions. The control bounds read as

$$ -60e \le u(t) \le 60e, \quad t \in [0, T]. $$

On the way to the target, the robot must overcome three holes in the ground. Instead of using very complicated expressions for $f_c(q)$ we model the holes as regions that the robot should not enter. This is achieved by constraints inside the OCP requiring that $p_{\text{foot}}$ is outside $n_e = 3$ ellipsoids:

$$ \left(\frac{p_{\text{foot},x} - x_{c,k}}{a_k}\right)^2 + \left(\frac{p_{\text{foot},z} - z_{c,k}}{b_k}\right)^2 \ge 1, k = 1, \dots, n_e. $$

By appropriately picking $a_k, b_k, x_{c,k}$ and $z_{c,k}$ the desired shapes are trivially selected. In our example, we pick $z_{c,k} = 0$, $a_k = 0.5$ (width of the hole), $b_k = 0.1$ (kept low, should not enforce unnecessarily high jumps). For the centers of the holes, we pick $x_{c,1} = 0.5$, $x_{c,2} = 1.5$, $x_{c,3} = 2.5$. We collect all path constraints (for the holes, on the kinematics, control, and state bounds) into the function $g(x, u) \ge 0$.

**Remark 10** Note that the constraints $g(x, u) \ge 0$ cannot become active if the corresponding normal velocity is nonzero, as opposed to activating a constraint $f_c(q) \ge 0$, since no state jump law is associated with path constraints in an OCP. This is one of the main differences between constraints that are part of the dynamics (equipped with a state jump law) and path constraints in the OCP.

We have now all ingredients to formulate an OCP of the form of Eq. (27). NOSNOC automatically reformulates the CLS into a time-freezing system and discretizes the OCP, such that we obtain a discrete-time problem of the form of (33). The resulting mathematical program with complementarity constraints is solved in a homotopy procedure with IPOPT [40] equipped with the MA57 linear solver [20]. The source code for this example is available in NOSNOC's repository [1]. The OCP is discretized with a FESD Radau-IIA scheme of order 5 [30]. We consider $N = 20$ control intervals with 3 intermediate integration steps on every interval.

For the initialization of the differential states, we take $y_0$ at every discretization node. All discrete-time control variables are initialized with zero. Hence, no information about the order, number, or timing of the nonsmooth transitions and jumps is provided. Treating the contact dynamics directly in the OCP and thus implicitly discovering all nonsmooth transitions is in the robotics community called contact implicit optimization [9,38]. The results of the optimization are shown in the right plot of Fig. 7. The approach finds an intuitive dynamic movement by solving only smooth NLP, without providing any hints about the order and number of nonsmooth transitions. The optimal torques are depicted in Fig. 8.

## 7 Conclusion and outlook

This paper introduced a novel time-freezing reformulation for transforming complementarity Lagrangian systems (CLS) with inelastic impacts and friction into piecewise smooth systems. This reduces the level of non-smoothness significantly. We prove solution equivalence under mild conditions and derive constructive ways to select the auxiliary ODE. Moreover, we show that we can obtain a solution to an optimal control problem with a CLS by solving the simpler OCP with time-freezing systems. We derive a reformulation of the time-freezing Filippov system into a dynamic complementarity system that allows the use of the high-accuracy Finite Elements
<figure><img src="image_9.png" /><figcaption>Fig. 7. Illustration of the robot kinematics (left), several frames of the solution of the discretized OCP (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. The optimal control input $u(t)$ in physical time $t$ obtained by solving the discretization of the optimal control problem.</figcaption></figure>

with Switch Detection (FESD) method for numerical optimal control [30]. To the best of our knowledge, this is the first reformulation that enables one to treat CLS with inelastic impacts as Filippov systems. The practicality of the discussed methods is demonstrated in an OCP considering a one-legged robot with frictional impact. All methods from this paper are implemented in the open-source package NOSNOC [28].

In further work, we aim to extend the ideas to multiple and simultaneous impacts, e.g., based on Moreau's impact law in second-order sweeping processes [7]. Moreover, to make this approach more practical, good initialization strategies and more sophisticated homotopy procedures for the mathematical programs with complementarity constraints would be useful.

## Acknowledgements

We thank Jan Carius from ETH Zürich in Switzerland, for providing details for the Capler robot.

## A Proof of Theorem 5

**PROOF.** The idea of the proof is to consider the different modes, in which the CLS and time-freezing system can be and to compare the solutions to establish the result of the theorem. A solution of the initial value problem given by the time-freezing system in Eq. (15) with $y(0) = y_0$ is denoted by $y_{\text{sol}}(\tau; y_0)$ for $\tau \in [0, \hat{\tau}]$. Similarly, for the CLS in Eq. (3) and $x(0) = x_0$ for $t(\tau) \in [0, t(\hat{\tau})]$ we use $x_{\text{sol}}(t(\tau); x_0)$. We must distinguish all possible cases, hence we split the proof into several parts.

**Part I (Unconstrained case).** Regard the case $f_c(q(\tau)) > 0, \tau \in [0, \hat{\tau}]$. This means that $y \in R_1$, $y' = f_1(y, u) = (f_{\text{ODE}}(x, u), 1)$, $\tau \in [0, \hat{\tau}]$. It holds that $t(\tau) = \int_0^\tau ds = \tau$ and by setting $\hat{\tau} = \tau_f$ we have $t(\tau_f) = t_f$. Note that $Ry' = Rf_1(y, u)$, is equivalent to $x' = f_{\text{ODE}}(x, u)$. Since $[0, \tau_f] = [0, t_f]$, this ODE has the same solution as $\dot{x} = f_{\text{ODE}}(x, u)$, therefore relation (18a) holds for $t \in [0, t_f]$. This means that $f_c(q(t)) > 0$ and $\lambda_n(t) = 0$ for $t \in [0, t_f]$. For the time-freezing system this means that $y(\tau) \notin \Sigma$ for $\tau \in [0, \tau_f]$, hence equation (18b) is also satisfied.

**Part II (Sliding mode/persistent contact).** Regard the case $f_c(q(0)) = 0$ and $n(q(0))^{\top}v(0) = 0$, i.e., $y(0) \in \Sigma$. Assume that $\varphi(x(\tau), u(\tau)) < 0$, $\tau \in [0, \tau_f]$. This means that $y(\tau) \in \Sigma$, $\tau \in [0, \tau_f]$, cf. Section 3.3. From Theorem 4 we have that $y' = \gamma(x, u)(f_{\text{DAE}}(x, u), 1)$ for $\tau \in [0, \tau_f]$ and $t'(\tau) = \gamma(x, u) > 0$, $\tau \in [0, \tau_f]$, thus $t(\tau_f) = t_f > 0$. On one hand, from $\frac{dy(t(\tau))}{d\tau} = \frac{dy}{dt}\frac{dt}{d\tau}$ we have that $R\frac{dy}{dt} = R\frac{1}{\gamma(x,u)}\gamma(x,u)(f_{\text{DAE}}(x,u),1) = f_{\text{DAE}}(x,u)$. On the other hand, for the CLS we have for $f_c(q(0)) = 0, n(q(0))^{\top}v(0) = 0, z(t) \ge 0$, $t \in [0, t_f]$. Consequently, the CLS reduces to the ODE $\dot{x} = f_{\text{DAE}}(x, u)$. Similar to part I, we conclude that (18a) holds. Since $f_c(q(t(\tau))) = 0$ and $n(q((\tau)))^{\top}v((\tau)) = 0$ for $\tau \in [0, \tau_f]$, the relation for $\lambda_n(t(\tau))$ in Eq. (18b) follows from Eq. (12) and (13).

**Part III (Leaving sliding mode).** Now we consider a similar scenario as in part II, with $y_0 \in \Sigma$, $\varphi(x(\tau), u(\tau)) \le 0$ with $\tau \in [0, \tau_e)$, $\tau_e < \tau_f$ (sliding mode) and $\varphi(x(\tau), u(\tau)) > 0$ for $\tau \in [\tau_e, \tau_f]$ (leaving sliding mode). Relations (18a) and (18b) hold for $\tau \in [0, \tau_e)$ by the same arguments as in part II. For $\tau \ge \tau_e$, following the arguments in Section 3.4, $y(\tau)$ leaves $\Sigma$ and $y(\tau) \in R_1$ for $\tau \in [0, \tau_e)$. We can apply the arguments of part I and establish the result of the theorem.

**Part IV (State jump).** This part regards the case of $\tau_s \in [0, \tau_f]$, i.e., $t_s \in [0, t_f]$. For $\tau \in [0, \tau_s)$ and $t \in [0, t(\tau_s^-))$ we can apply Part I of the proof by simply setting $\hat{\tau} = \tau_s$ and deduce that (18a) and (18b) hold. For $\tau = \tau_s$ we have $f_c(q(\tau_s)) = 0$ and $n(q(\tau_s))^{\top}v(\tau_s) < 0$. Consequently, $y \in R_2$ and $y' = f_{\text{aux,n}}(y)$. The assumption $f_c(q(\tau_f)) \ge 0$ and $n(q(t_f))^{\top}v(t_f) \ge 0$ ensures that the time evolution of $y'(\tau) = f_{\text{aux,n}}(y(\tau))$ is finished in $[\tau_s, \tau_f]$, i.e., $\tau_r \le \tau_f$. From the proof of Proposition 3 we know that by construction $q(\tau) = q(\tau_s) =: q_s$, $\tau \in [\tau_s, \tau_r]$. Consequently, $f_c(q(\tau)) = 0$, $\tau \in [\tau_s, \tau_r]$. For $v(\tau)$, from (16) we obtain
that:

$$v(\tau_r) = v(\tau_s) + \int_{\tau_s}^{\tau_r} M(q(\tau))^{-1} n(q(\tau)) a_n d\tau. \quad (A.1)$$

Multiplying both sides with $n(q_s)^\top$ from the left and noting that $M^{-1}(q(\tau))n(q(\tau))$ is constant since $q(\tau) = q_s$, $\tau \in [\tau_s, \tau_r]$, we have

$$\underbrace{n(q_s)^\top v(\tau_r)}_{=0} - \underbrace{n(q_s)^\top v(\tau_s)}_{=D(q_s)} = \underbrace{n(q_s)^\top M(q_s)^{-1} n(q_s)}_{\underbrace{\int_{\tau_s}^{\tau_r} a_n d\tau,}_{=:\Lambda_1}} \\ \Lambda_1 = - \frac{n(q_s)^\top v(\tau_s)}{D(q_s)} > 0.$$

Next, we look at the post-impact states of the CLS and compare it to the solution of the time-freezing system. Since in CLS, $v(t)$ is a function of bounded variation [7], we have that $q(t)$ is a continuous function. Thus, $q(t_s^+) = q(t_s^-) = q(t_s)$. Furthermore, notice that $q(t_s) = q_s$ which implies $n(q(t_s))^\top M^{-1}(q(t_s))n(q(t_s)) = D(q_s)$. Examining,

$$v(t_s^+) = v(t_s^-) + \lim_{\epsilon \to 0} \underbrace{\int_{t_s-\epsilon}^{t_s+\epsilon} f_v(q(t), v(t)) dt}_{=0} \quad (A.2) \\ + \lim_{\epsilon \to 0} \int_{t_s-\epsilon}^{t_s+\epsilon} M(q(t))^{-1} n(q(t)) \lambda_n(t) dt,$$

and multiplying both sides with $n(q_s)^\top$ from the left, introducing $\Lambda_2 := \lim_{\epsilon \to 0} \int_{t_s-\epsilon}^{t_s+\epsilon} \lambda_n(t) dt$, we conclude that

$$\Lambda_2 = - \frac{n(q(t_s))^\top v(t_s^-)}{D(q_s)} = \Lambda_1. \quad (A.3)$$

By comparing (A.1) and (A.2), due to the last relation we conclude that $v(\tau_r) = v(t_s^+) =: v_s$. Furthermore, by Proposition 3 we have $f_c(q(\tau_r)) = 0$ and $n(q(\tau_r))^\top v(\tau_r) = 0$. Since $t'(\tau) = 0$ with $\tau \in [\tau_s, \tau_r]$, it follows that $t(\tau_r) = t(\tau_s) = t_s$. Consequently,

$$f_c(q(t_s)) = f_c(q(\tau_r)) = 0, \quad (A.4a)$$

$$n(q(t_s))^\top v(t_s^+) = n(q(\tau_r))^\top v(\tau_r) = 0. \quad (A.4b)$$

Let $y_s := (q_s, v_s, t_s)$. Note that $y_{sol}(\tau - \tau_r, y_s) = y(\tau, y_0)$ for $\tau \in [\tau_r, \tau_f]$. Likewise, $x_{sol}(t - t_s, x_s) = x(t, x_0)$ for $t \in (t_s, t_f]$, with $x_s = Ry_s$. The two initial value problems are initialized with the same initial condition. Since (A.4) holds we can apply Theorem 4 for $y \in \Sigma$. Therefore, by using the arguments of parts II or III (depending on $\varphi(x(\tau), u(\tau)))$, we deduce that (18a) and (18b) hold on $[\tau_r, \tau_f]$. Additionally, for $\tau \in (\tau_s, \tau_r]$ we have $t = t_s$ and Eq. (19) follows directly from (A.3).

Part V (Summary). Parts I-IV cover all possible modes of the CLS and the time-freezing system: evolution according to $f_{ODE}$ (Part I), evolution on $\Sigma$ according to $f_{DAE}$ without leaving it (Part II), leaving $\Sigma$ and continuing to evolve according to $f_{ODE}$ (Part III), and the state jump (Part IV). To regard any other possible sequence of mode on $[0, \tau_f]$, the time interval is simply split into sub-intervals with the different mode, and we apply subsequently the arguments from Parts I-IV to verify that (18a) and (18b) hold for $t \neq t_s$ and (19) for $t = t_s$. This completes the proof. □

## B Proof of Theorem 7

**PROOF.** For the first part of the proof, we assume that $b(q)^\top v > 0$. This means that $y \notin R_3$ and it follows that $\theta_3 = 0$. From the conditions $c_2(y) = 0$ and $\theta_1 + \theta_2 = 1$ we can compute $\theta_1$ and $\theta_2$. By using the fact that $\nabla n(q)^\top M(q)^{-1}b(q) = 0$ and following similar lines as in the proof of Theorem 4 we can compute that $y' = (v, f_v(q, v, u) - M(q)^{-1}D(q)^{-1}\varphi(x, u)(n(q) - \mu b(q)))$. By similar reasoning for $b(q)^\top v < 0$ we obtain that $y' = (v, f_v(q, v, u) - M(q)^{-1}D(q)^{-1}\varphi(x, u)(n(q) + \mu b(q)))$. Therefore, it holds for $v_t \neq 0$ that $y' = \gamma(x, u)(f_{Slip}(x, u), 1)$, whereby the sign(·) in the r.h.s. of $f_{Slip}(x, u)$ accounts for the sign of $v_t$.

In the second part, we have $v_t = b(q)^\top v = 0$. Together with the assumption that $y \in \Sigma$ it follows that $y \in \partial R_i$, for all $i = 1, 2, 3$, hence no $\theta_i$ can be set to be zero a priori. Since the vectors $b(q)$ and $n(q)$ are orthogonal in the kinetic metric, it can be seen that the matrix $\tilde{D}(q)$ is a diagonal matrix. We denote the first and second entries on its diagonal by $D(q)$ (cf. (11a)) and $D_t(q)$, respectively. Similarly, the first and second components of the vector $\tilde{\varphi}(x, u)$ are denoted by $\varphi(x, u)$ and $\varphi_t(x, u)$, respectively. This, we obtain that $\lambda_n = -D(q)^{-1}\varphi(x, u)$ and $\lambda_t = D_t(q)^{-1}\varphi_t(x, u)$ holds.

From the condition $\frac{dc_2(y)}{d\tau} = 0$ and the definition of the time-freezing system we compute that

$$\varphi(x, u)\theta_1 + D(q)a_n(\theta_2 + \theta_3) = 0.$$

Using $\sum_{i=1}^3 \theta_i = 1$ and the last equation we can compute that:

$$\theta_1 = \frac{D(q)a_n}{D(q)a_n - \varphi(x, u)} = \gamma(x, u), \quad (B.1a)$$

$$\theta_2 + \theta_3 = \frac{-\varphi(x, u)}{D(q)a_n - \varphi(x, u)} \quad (B.1b)$$

Next, we use the condition $\frac{dc_3(y)}{d\tau} = \frac{d(b(q)^\top v)}{d\tau} = 0$ and compute:

$$\varphi_t(x, u)\theta_1 + \mu D_t(q)a_n(-\theta_2 + \theta_3) = 0,$$
and by using (B.1a) we establish the relation:

$$-\theta_2 + \theta_3 = -\gamma(x, u) \frac{\varphi_t(x, u)}{\mu D_t(q) a_n} \quad (B.2)$$

Next, we compute the sliding mode vector field of the time-freezing system $y' = \sum_{i=1}^{3} \theta_i f_i$. By rearranging the terms we obtain

$$y' = \theta_1 \begin{bmatrix} v \\ f_v \\ 1 \end{bmatrix} + (\theta_2 + \theta_3) \begin{bmatrix} \mathbf{0}_{n_q, 1} \\ M(q)^{-1} n(q) a_n \\ 0 \end{bmatrix} + (-\theta_2 + \theta_3) \begin{bmatrix} \mathbf{0}_{n_q, 1} \\ M(q)^{-1} b(q) \mu a_n \\ 0 \end{bmatrix}$$

Now, by multiplying the second term by $D(q)^{-1}D(q)$, plugging in the expressions for the sums (B.1b) and (B.2) and by comparing it to the r.h.s. of Eq. (23) we can conclude that $y' = \gamma(x, u)(f_{\text{Stick}}(x, u), 1)$. This completes the proof. □

## References

[1] NOSNOC. MATLAB: https://github.com/nurkanovic/nosnoc, python: https://github.com/FreyJo/nosnoc_python.

[2] Vincent Acary and Bernard Brogliato. *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*. Springer Science & Business Media, 2008.

[3] Vincent Acary, Hidde De Jong, and Bernard Brogliato. Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems. *Physica D: Nonlinear Phenomena*, 269:103–119, 2014.

[4] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – a software framework for nonlinear optimization and optimal control. *Mathematical Programming Computation*, 11(1):1–36, 2019.

[5] Mihai Anitescu, Paul Tseng, and Stephen J. Wright. Elastic-mode algorithms for mathematical programs with equilibrium constraints: global convergence and stationarity properties. *Mathematical Programming*, 110(2):337–371, 2007.

[6] Hans G. Bock and K. J. Plitt. A multiple shooting algorithm for direct solution of optimal control problems. In *Proceedings of the IFAC World Congress*, pages 242–247. Pergamon Press, 1984.

[7] Bernard Brogliato. *Nonsmooth Mechanics: Models, Dynamics and Control*. Springer International Publishing Switzerland, 2016. Third edition.

[8] Bernard Brogliato and Aneel Tanwani. Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability. *SIAM Review*, 62(1):3–129, 2020.

[9] Jan Carius, René Ranftl, Vladlen Koltun, and Marco Hutter. Trajectory optimization with implicit hard contacts. *IEEE Robotics and Automation Letters*, 3(4):3316–3323, 2018.

[10] Arrigo Cellina and António Ornelas. Existence of solutions to differential inclusions and to time optimal control problems in the autonomous case. *SIAM journal on control and optimization*, 42(1):260–265, 2003.

[11] Frank H. Clarke. Optimal solutions to differential inclusions. *Journal of Optimization Theory and Applications*, 19(3):469–478, 1976.

[12] Giovanni Colombo, Boris S. Mordukhovich, and Dao Nguyen. Optimization of a perturbed sweeping process by constrained discontinuous controls. *SIAM Journal on Control and Optimization*, 58(4):2678–2709, 2020.

[13] Luca Dieci and Luciano Lopez. Sliding motion on discontinuity surfaces of high co-dimension. a construction for selecting a filippov vector field. *Numerische Mathematik*, 117(4):779–811, 2011.

[14] Aleksei F. Filippov. *Differential Equations with Discontinuous Righthand Sides*. Dordrecht, The Netherlands: Kluwer, 1988.

[15] Christian Gehring. Operational space control of single legged hopping. Master’s thesis, Eidgenössische Technische Hochschule Zürich, Autonomous Systems Lab, 2011.

[16] Lei Guo and Jane J. Ye. Necessary optimality conditions for optimal control problems with equilibrium constraints. *SIAM Journal on Control and Optimization*, 54(5):2710–2733, 2016.

[17] Jonas Hall, Armin Nurkanović, Florian Messerer, and Moritz Diehl. A sequential convex programming approach to solving quadratic programs and optimal control problems with linear complementarity constraints. *IEEE Control Systems Letters*, 6:536–541, 2022.

[18] Mathew Halm and Michael Posa. Set-valued rigid body dynamics for simultaneous frictional impact. *arXiv preprint arXiv:2103.15714*, 2021.

[19] Taylor A. Howell, Simon Le Cleac'h, J. Zico Kolter, Mac Schwager, and Zachary Manchester. Dojo: A differentiable simulator for robotics. *arXiv preprint arXiv:2203.00806*, 2022.

[20] HSL. A collection of Fortran codes for large scale scientific computation. http://www.hsl rl.ac.uk, 2011.

[21] Jemin Hwangbo, Vassilios Tsounis, Hendrik Kolvenbach, and Marco Hutter. Cable-driven actuation for highly dynamic robotic systems. In *2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, pages 8543–8550. IEEE, 2018.

[22] Jisu Kim, Hansung Cho, Artem Shamsuarov, Hyungbo Shim, and Jin H. Seo. State estimation strategy without jump detection for hybrid systems using gluing function. In *53rd IEEE International Conference on Decision and Control (CDC)*, pages 139–144. IEEE, 2014.

[23] Christian Kirches. A Numerical Method for Nonlinear Robust Optimal Control with Implicit Discontinuities and an Application to Powertrain Oscillations. Diploma thesis, University of Heidelberg, October 2006.

[24] Boris S. Mordukhovich and Yuan Tian. Implicit euler approximation and optimization of one-sided lipschitzian differential inclusions. *Nonlinear analysis and optimization*, 659:165–188, 2016.

[25] Jean Jacques Moreau. Application of convex analysis to some problems of dry friction. In *Trends in applications of pure mathematics to mechanics*, pages 263–280. Pitman, 1977.

[26] Ngoc Son Nguyen and Bernard Brogliato. Comparisons of multiple-impact laws for multibody systems: Moreau’s law, binary impacts, and the LZB approach. In *Advanced Topics in Nonsmooth Dynamics*, pages 1–45. Springer, 2018.
[27] Armin Nurkanović, Sebastian Albrecht, and Moritz Diehl. Limits of MPCC formulations in direct optimal control with nonsmooth differential equations. In *2020 European Control Conference (ECC)*, pages 2015–2020, 2020.

[28] Armin Nurkanović and Moritz Diehl. NOSNOC: A software package for numerical optimal control of nonsmooth systems. *IEEE Control Systems Letters (L-CSS)*, 6:3110–3115, 2022.

[29] Armin Nurkanović, Tommaso Sartor, Sebastian Albrecht, and Moritz Diehl. A time-freezing approach for numerical optimal control of nonsmooth differential equations with state jumps. *IEEE Control Systems Letters (L-CSS)*, 5(2):439–444, 2021.

[30] Armin Nurkanović, Mario Sperl, Sebastian Albrecht, and Moritz Diehl. Finite elements with switch detection for direct optimal control of nonsmooth systems. *arXiv preprint arXiv:2205.05337*, 2022.

[31] Armin Nurkanović and Moritz Diehl. Continuous optimization for control of hybrid systems with hysteresis via time-freezing. *IEEE Control Systems Letters*, 6:3182–3187, 2022.

[32] Michael Posa, Cecilia Cantu, and Russ Tedrake. A direct method for trajectory optimization of rigid bodies through contact. *The International Journal of Robotics Research*, 33(1):69–81, 2014.

[33] James B. Rawlings, David Q. Mayne, and Moritz Diehl. *Model Predictive Control: Theory, Computation, and Design*. Nob Hill, 2nd edition, 2017.

[34] David E. Stewart. A high accuracy method for solving ODEs with discontinuous right-hand side. *Numerische Mathematik*, 58(1):299–328, 1990.

[35] David E. Stewart. Rigid-body dynamics with friction and impact. *SIAM Review*, 42(1):3–39, 2000.

[36] David E. Stewart and Mihai Anitescu. Optimal control of systems with discontinuous differential equations. *Numerische Mathematik*, 114(4):653–695, 2010.

[37] David E. Stewart and Jeffrey C. Trinkle. An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and coulomb friction. *International Journal for Numerical Methods in Engineering*, 39(15):2673–2691, 1996.

[38] Yuval Tassa, Tom Erez, and Emanuel Todorov. Synthesis and stabilization of complex behaviors through online trajectory optimization. In *2012 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, pages 4906–4913. IEEE, 2012.

[39] Alexandre Vieira, Bernard Brogliato, and Christophe Prieur. Quadratic optimal control of linear complementarity systems: First-order necessary conditions and numerical analysis. *IEEE Transactions on Automatic Control*, 65(6):2743–2750, 2020.

[40] Andreas Wächter and Lorenz T. Biegler. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming*, 106(1):25–57, 2006.

[41] Yaofeng Desmond Zhong, Jiequn Han, and Georgia Olympia Brikis. Differentiable physics simulations with contacts: Do they have correct gradients wrt position, velocity and control? In *ICML 2022 2nd AI for Science Workshop*.

[42] Viktor F. Zhuravlev. Equations of motion of mechanical systems with ideal one-sided constraints. *Prikladnaia Matematika i Mekhanika*, 42:781–788, 1978.