# Multivalued Hamiltonian Systems with Multivalued Dissipation: Analysis of the Backward-Euler Discretisation

Fernando Castaños, Felix A Miranda-Villatoro, Bernard Brogliato

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

## To cite this version:

Fernando Castaños, Felix A Miranda-Villatoro, Bernard Brogliato. Multivalued Hamiltonian Systems with Multivalued Dissipation: Analysis of the Backward-Euler Discretisation. 2024. hal-04625231

HAL Id: hal-04625231

https://inria.hal.science/hal-04625231

Preprint submitted on 26 Jun 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Multivalued Hamiltonian Systems with Sliding Motions: Analysis of the Backward-Euler Discretisation

Fernando Castaños, Felix A Miranda-Villatoro, Bernard Brogliato

► To cite this version:

Fernando Castaños, Felix A Miranda-Villatoro, Bernard Brogliato. Multivalued Hamiltonian Systems with Sliding Motions: Analysis of the Backward-Euler Discretisation. Automatica, 2025, 182, pp.112567. 10.1016/j.automatica.2025.112567. hal-04625231v2

HAL Id: hal-04625231

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

Submitted on 5 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Multivalued Hamiltonian Systems with Multivalued Dissipation: Analysis of the Backward-Euler Discretisation*

Fernando Castaños<sup>b</sup>, Felix A. Miranda-Villatoro<sup>a</sup>, Bernard Brogliato<sup>a</sup>

## Abstract

This article is mainly concerned with the time-discretisation of multivalued Hamiltonian systems with multivalued dissipation, a special class of differential inclusions. Two classes of set-valued Hamiltonian systems are considered, depending on the dissipation function being position or momentum dependent. The backward-Euler discretisation is analyzed in both cases: well-posedness of the generalized equation obtained after discretisation is proved, then stability of fixed-points is tackled. The well-known sliding-mode twisting and super-twisting algorithms, as well as an example from Contact Mechanics, illustrate the theoretical developments.

## Key words:

Hamiltonian systems, differential inclusions, backward Euler time-discretisation, Lyapunov stability, twisting algorithm, super-twisting algorithm.

## 1 Introduction

The Hamiltonian formalism encompasses a large class of systems with diverse applications such as mechanics [3], circuit theory [29], thermodynamics [14], economy [32], and economy [36]. In the classical formulation, a Hamiltonian system is characterized by two objects: a two-form on a smooth manifold (called a symplectic form) and a smooth function (called a Hamiltonian function) [3]. The main property of Hamiltonian systems is the simultaneous preservation of the symplectic form and the Hamiltonian function along the system trajectories [17]. For electrical or mechanical systems, the Hamiltonian function corresponds to the energy stored in the system, so preservation of the Hamiltonian amounts to energy conservation.

The Hamiltonian approach has been generalised in several directions. An important one in the control literature is the port-Hamiltonian formalism, where Hamilto-

nian systems are endowed with dissipative elements and external ports so that they can interact with the environment and, e.g., a controller [28,?]. Port-Hamiltonian systems are not necessarily energy-conserving, but satisfy the more general property of passivity and are amenable to passivity-based control [39]. Another interesting generalisation consists in allowing systems to be described by differential inclusions in place of differential equations. There are several reasons for doing so, the main one being the possibility to treat nondifferentiable Hamiltonian functions [13]. This paper considers systems with possibly nonsmooth Hamiltonian functions and nonsmooth Rayleigh dissipation potentials (systems with additional external ports will be addressed in future contributions). Our main motivation is that this class of systems contains several higher-order sliding-mode controllers, such as the twisting and super-twisting algorithms [23], which have been difficult to classify in the literature. The primary interest in such controllers is the finite-time stability of the closed-loop system.

An important issue concerning dynamical systems is time discretisation. In the context of pure Hamiltonian systems, research has focused on the development of

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

<footnote><sup>b</sup> Automatic Control Department, Cinvestav-IPN, 2508 Av. IPN, 07360, Mexico City, Mexico</footnote>

<footnote>* This paper was not presented at any IFAC meeting.
Email addresses: felix.miranda-villatoro@inria.fr
(Felix A. Miranda-Villatoro),
bernard.brogliato@inria.fr (Bernard Brogliato).</footnote>
# Multivalued Hamiltonian Systems
with Sliding Motions:
Analysis of the Backward-Euler Discretisation ⋆

Fernando Castaños<sup>b</sup>, Félix Miranda-Villatoro<sup>a</sup>, Bernard Brogliato<sup>a</sup>

**Abstract**

This article is mainly concerned with the time-discretisation of multivalued Hamiltonian systems with multivalued dissipation, a special class of differential inclusions. Two classes of set-valued Hamiltonian systems are considered, depending on whether the dissipation function is position- or momentum-dependent. The backward-Euler discretisation is analysed in both cases: the well-posedness of the generalised equation obtained after discretisation is proved, and then *finite-time stability* of fixed points is tackled. The well-known twisting and super-twisting sliding-mode algorithms, as well as an example from Contact Mechanics and dynamical optimisation, illustrate the theoretical developments.

**Key words:**

Hamiltonian systems, differential inclusions, backward-Euler time-discretisation, Lyapunov stability, twisting algorithm, super-twisting algorithm, sliding mode, dynamical optimisation, robust control.

## 1 Introduction

The Hamiltonian formalism encompasses a large class of systems with diverse applications such as mechanics [2], circuit theory [7], thermodynamics [13], ecology [26], and economy [31]. In the classical formulation, a Hamiltonian system is characterised by two objects: a two-form on a smooth manifold (called a symplectic form) and a smooth function (called a Hamiltonian function) [2]. The main property of Hamiltonian systems is the simultaneous preservation of the symplectic form and the Hamiltonian function along the system trajectories [15]. For electrical or mechanical systems, the Hamiltonian function corresponds to the energy stored in the system, so preservation of the Hamiltonian amounts to energy conservation. The Hamiltonian approach has been generalised in several directions. An important one in the control literature is the port-Hamiltonian formalism, where Hamiltonian systems are endowed with dissipative elements and external ports so that they can interact with the environment and, e.g., a controller.

Port-Hamiltonian systems are not necessarily energy-conserving, but satisfy the more general property of passivity and are amenable to passivity-based control [33]. Another interesting generalisation consists in allowing systems to be described by differential inclusions (DI) in place of differential equations. There are several reasons for doing so, the main one being the possibility to treat non-differentiable Hamiltonian functions [11]. This paper considers systems with possibly nonsmooth Hamiltonian functions and nonsmooth Rayleigh dissipation potentials. Our main motivation is that this class of systems contains some higher-order sliding-mode controllers, such as the second-order nested, twisting, and super-twisting algorithms [18]. The primary interests in such controllers are the finite-time stability and the robustness of the closed-loop system.

An important issue concerning dynamical systems is time discretisation. In the context of pure Hamiltonian systems, research has focused on the development of discrete-time maps that approximate the Hamiltonian flow while ensuring that the symplectic form and a perturbed Hamiltonian function are invariant along dissipative elements and external ports so that they can interact with the environment and, e.g., a controller.

Email addresses: fernando.castanos@cinvestav.mx (Fernando Castaños), felix.miranda-villatoro@inria.fr, (Félix Miranda-Villatoro), bernard.brogliato@inria.fr (Bernard Brogliato).

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

<footnote><sup>b</sup> Automatic Control Department, Cinvestav-IPN, 2508 Av. IPN, 07360, Mexico City, Mexico</footnote>

<footnote>* This paper was not presented at any IFAC meeting.</footnote>
discrete-time maps that approximate the Hamiltonian flow while ensuring that the symplectic form and a perturbed Hamiltonian function are invariant along discrete orbits [17]. In the context of port-Hamiltonian systems, the usual objective is to preserve either passivity [37, 20, 16, 2, 31, 27] or the closely related Dirac structure [38].

In the framework of this work, multivalued Hamiltonian systems result from closing the loop with multivalued controls. The objective is to derive well-posed discretisation schemes that preserve finite-time stability. The practical incentive is that the corresponding discrete-time control law is immediately revealed once an appropriate discrete-time model is obtained.

The following section contains the minimal background. The main contributions of this paper are distributed in Sections 3 and 4. These are:

* A unified framework for proving the finite-time stability of the twisting and super-twisting controllers. The proposed framework, multivalued Hamiltonian systems with dissipation, can be used to study other systems and creates the possibility for devising new sliding-mode algorithms (Section 3).

* A criterion for establishing finite-time stability in discrete time (Proposition 2 in Section 4).

* Sufficient conditions for the backward-Euler discretisation of multivalued Hamiltonian systems to inherit finite-time convergence (Theorem 3 in Subsection 4.1 and Theorem 5 in Subsection 4.3). The twisting and super-twisting algorithms satisfy such conditions.

* A splitting method, inspired by the Douglas-Rachford method in Optimisation, for approximating the proposed discrete-time maps (and hence the corresponding discrete-time controllers) of multivalued Hamiltonian systems with intricate Hamiltonian or dissipation functions (Theorem 4 in Subsection 4.2). The method is applied to the twisting algorithm and shown to preserve finite-time stability, even though it is an approximation.

Finally, Section 5 contains some conclusions.

## 2 Background

In this section we provide the necessary material to recall a useful invariance theorem for differential inclusions and level sets of nonsmooth functions [4]. The following preliminary results and definitions are taken from classical references [33, 34, 18, 6]. We consider convex functions $V: \mathbb{R}^n \to \mathbb{R}$. Note that, in contrast with the above-mentioned references, we exclude the value $+\infty$ from the range of the functions we consider. Thus, convex functions are automatically continuous and proper (indeed,

their domains are equal to $\mathbb{R}^n$). The convex subdifferential of $V$ at the point $x$ is given by

$$ \partial V(x) := \{\eta \in \mathbb{R}^n \mid \langle \eta, \xi - x \rangle \le V(\xi) - V(x), \\ \text{for all } \xi \in \mathbb{R}^n \}. \quad (1) $$

By [18, Remark 4.1.7], $\partial V$ also has domain $\mathbb{R}^n$ and is locally bounded. The graph of a set-valued map $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^m$ is the set

$$ \text{gph } \mathbf{F} := \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^m \mid y \in \mathbf{F}(x)\}. \quad (2) $$

The set of zeros of $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$ is denoted by

$$ \text{Zero } \mathbf{F} := \mathbf{F}^{-1}(0) = \{\xi \in \mathbb{R}^n \mid 0 \in \mathbf{F}(\xi)\}. \quad (3) $$

Let $r \in [1, +\infty)$, the $r$-norm of a vector $x \in \mathbb{R}^n$ is denoted as

$$ \|x\|_r = \left( \sum_{i=1}^{n} |x_i|^r \right)^{\frac{1}{r}}. \quad (4) $$

By convention, we set $\|x\|_{\infty} = \max_{i \in 1, \dots, n} |x_i|$. The open ball, with respect to the $r$-norm, centered at the origin and having radius 1, is denoted by

$$ \mathcal{B}_r = \{\xi \in \mathbb{R}^n \mid \|\xi\|_r < 1\}. \quad (5) $$

To avoid cumbersome notation, $\| \cdot \|_2$ and $\mathcal{B}_2$ are denoted simply as $\| \cdot \|$ and $\mathcal{B}$, respectively. The closure of a set $S$ is denoted $\mathrm{cl} S$.

$\mathbf{F}$ is said to be outer semicontinuous (osc) at $x^* \in \mathrm{dom} \mathbf{F}$ if, for any $y \notin \mathbf{F}(x^*)$, there are neighborhoods of $x^*$ and $y$ (denoted respectively as $\mathcal{N}(x^*)$ and $\mathcal{N}(y)$) such that $\mathcal{N}(x^*) \cap \mathbf{F}^{-1}(\mathcal{N}(y)) = \emptyset$. When $\mathbf{F}$ is closed-valued and locally bounded, outer semicontinuity at $x^*$ is equivalent to the following (see e.g., [11, Proposition 2.5.24]): For every $\varepsilon > 0$, there exists a neighborhood $\mathcal{N}(x^*)$ such that $x \in \mathcal{N}(x^*)$ implies that

$$ \mathbf{F}(x) \subset \mathbf{F}(x^*) + \varepsilon \mathcal{B}. \quad (6) $$

A set-valued operator $\mathbf{F}$ is said to be monotone (respectively, strongly monotone) if, for any two pairs $(x_i, y_i)$, $x_i \in \mathbb{R}^n$, $y_i \in \mathbf{F}(x_i)$, $i = 1, 2$,

$$ \langle x_1 - x_2, y_1 - y_2 \rangle \ge 0 \quad (7) $$

(respectively, if

$$ \langle x_1 - x_2, y_1 - y_2 \rangle \ge \gamma \|x_1 - x_2\|^2, \quad (8) $$

for some $\gamma > 0$). The operator $\mathbf{F}$ is maximal monotone if it is monotone and its graph is not strictly contained inside the graph of any other monotone operator. It is
crete orbits [15]. In the context of port-Hamiltonian systems, the usual objective is to preserve either passivity [17,1,23] or the closely related Dirac structure [32]. In the framework of this work, multivalued Hamiltonian systems may result from closing the loop with multivalued laws like those arising in sliding-mode control. The objective is to derive well-posed discretisation schemes that preserve finite-time stability. The practical incentive is that the corresponding discrete-time control law is immediately revealed once an appropriate discrete-time model is obtained.

Section 2 contains the minimal background. The main contributions are distributed in Sections 3 and 4:

* A unified framework for proving the finite-time stability of the second-order nested, twisting and super-twisting controllers. The proposed framework, multivalued Hamiltonian systems with dissipation, can be used to study other systems and creates the possibility for devising new sliding-mode algorithms (Section 3).

* A criterion for establishing finite-time stability in discrete time (Proposition 2 in Section 4).

* Sufficient conditions for the backward-Euler discretisation of multivalued Hamiltonian systems to inherit finite-time convergence (Theorem 3 in Section 4.1 and Theorem 5 in Section 4.3). The twisting and super-twisting algorithms satisfy such conditions.

* A splitting method, inspired by the Douglas-Rachford method in Optimisation, for approximating the proposed discrete-time maps (and hence the corresponding discrete-time controls) of multivalued Hamiltonian systems with intricate Hamiltonian or dissipation functions (Theorem 4 in Section 4.2). The method is applied to the twisting scheme: it preserves finite-time stability, even though it is an approximation.

Finally, Section 5 presents some open problems, and conclusions are in Section 6.

## 2 Background

The following preliminary results and definitions are taken from classical references [3,29,16,5].

**Notation.** Let $r \in [1, +\infty)$. The r-norm of a vector $x \in \mathbb{R}^n$ is denoted as $\|x\|_r = (\sum_{i=1}^n |x_i|^r)^{\frac{1}{r}}$. The open ball, with respect to the r-norm, centred at the origin and having radius 1, is written as $B_r = \{\xi \in \mathbb{R}^n | \|\xi\|_r < 1\}$. To avoid cumbersome notation, $\| \cdot \|_2$ and $B_2$ are written simply as $\| \cdot \|$ and $B$, respectively. We use $I_n$ for the $n \times n$-dimensional identity matrix. The closure of a set $S$ is denoted $\text{cl } S$. We denote the distance from $x \in \mathbb{R}^n$ to the set $S \subset \mathbb{R}^n$ by $\text{Dist}(x, S) = \inf_{\xi \in S} \|x - \xi\|$, and the projection of $x \in \mathbb{R}^n$ onto the closed set $S \subset \mathbb{R}^n$ by

$$ \text{Proj}(x, S) = \arg \min_{\xi \in S} \|x - \xi\|. $$

Continuous differentiability is denoted as $C^1$. The gradient of a $C^1$ function $f: \mathbb{R}^n \to \mathbb{R}$ at $x \in \mathbb{R}^n$ is denoted by $\nabla f(x)$, and the signed bracket by $[x]^\gamma = |x|^\gamma \text{sgn}(x)$, $\gamma \ge 0$, $x \in \mathbb{R}$.

## 2.1 Convex Analysis and Set-valued Functions

We consider convex functions $V: \mathbb{R}^n \to \mathbb{R}$. Note that we exclude the value $+\infty$ from the range of functions we consider. Thus, convex functions are automatically proper (recall that $V$ is proper if $V(x) < +\infty$ for some $x \in \mathbb{R}^n$) and continuous [29, Corollary 2.36] The convex subdifferential of $V$ at the point $x$ is given by

$$ \partial V(x) = \{\eta \in \mathbb{R}^n | \langle \eta, \xi - x \rangle \le V(\xi) - V(x) \text{ for all } \xi \in \mathbb{R}^n\}. $$

By [16, Remark 4.1.7], $\partial V$ also has domain $\mathbb{R}^n$ and is locally bounded. Note that, whenever $V$ is $C^1$, $\partial V(x)$ reduces to $\{\nabla V(x)\}$. The graph of a set-valued map $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$ is the set

$$ \text{gph} \mathbf{F} = \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^n | y \in \mathbf{F}(x)\}. $$

The set of zeros of $\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$ is denoted by

$$ \text{Zero} \mathbf{F} = \mathbf{F}^{-1}(0) = \{\xi \in \mathbb{R}^n | 0 \in \mathbf{F}(\xi)\}. \quad (1) $$

$\mathbf{F}$ is said to be outer semicontinuous (osc) at $x^* \in \text{dom} \mathbf{F}$ if, for any $y \notin \mathbf{F}(x^*)$, there are neighbourhoods of $x^*$ and $y$ (denoted respectively as $\mathcal{N}(x^*)$ and $\mathcal{N}(y)$) such that $\mathcal{N}(x^*) \cap \mathbf{F}^{-1}(\mathcal{N}(y)) = \emptyset$. When $\mathbf{F}$ is closed-valued and locally bounded, outer semicontinuity at $x^*$ is equivalent to the following (see e.g., [10, Proposition 2.5.24]): For every $\varepsilon > 0$, there exists a neighbourhood $\mathcal{N}(x^*)$ such that $x \in \mathcal{N}(x^*)$ implies that $\mathbf{F}(x) \subset \mathbf{F}(x^*) + \varepsilon\mathcal{B}$.

$\mathbf{F}$ is said to be monotone (respectively, strongly monotone) if, for any two pairs $(x_i, y_i)$, $x_i \in \mathbb{R}^n$, $y_i \in \mathbf{F}(x_i)$, $i = 1, 2$, $\langle x_1 - x_2, y_1 - y_2 \rangle \ge 0$ (respectively, if $\langle x_1 - x_2, y_1 - y_2 \rangle \ge \gamma \|x_1 - x_2\|^2$, for some $\gamma > 0$). The operator $\mathbf{F}$ is maximal monotone if it is monotone and its graph is not strictly contained inside the graph of any other monotone operator. In finite dimensions, maximal monotone operators are osc [10, Proposition 4.2.1].

The proximal map associated with a convex function $V: \mathbb{R}^n \to \mathbb{R}$ is

$$ \text{Prox}_V(x) = \arg \min_{w \in \mathbb{R}^n} \left\{ V(w) + \frac{1}{2} \|w - x\|^2 \right\}. $$

Equivalently, $p_x = \text{Prox}_V(x)$ if, and only if, $x \in (\text{I}_d + \partial V)(p_x)$, where $\text{I}_d$ denotes the identity map. In other words, $\text{Prox}_V(x) = (\text{I}_d + \partial V)^{-1}(x)$. Recall that, by the maximal monotonicity of $V$, $\text{Prox}_V$ is a well-defined single-valued operator [5, Corollary 23.10]
well known that, in finite dimensions, maximal mono-
tone operators are osc [11, Proposition 4.2.1].

The proximal map associated with a convex function
$V : \mathbb{R}^n \to \mathbb{R}$ is

$$
\operatorname{Prox}_V(x) = \arg \min_{w \in \mathbb{R}^n} \left\{ V(w) + \frac{1}{2} \|w - x\|^2 \right\}. \quad (9)
$$

Equivalently, $p_x = \text{Prox}_V(x)$ if, and only if, $x \in (\text{I}_d + \partial V)(p_x)$, where $\text{I}_d$ denotes the identity map. In other words, $\text{Prox}_V(x) = (\text{I}_d + \partial V)^{-1}(x)$.

Consider a differential inclusion

$$
\dot{x}(t) \in \mathbf{F}(x(t)), \quad x(0) = x_0. \tag{10}
$$

We assume that the multivalued vector field **F** is such that, for each $x_0 \in \mathbb{R}^n$, there exists a unique absolutely continuous solution $x(\cdot)$. Recall that the set-valued derivative of a convex function $V : \mathbb{R}^n \to \mathbb{R}$ along **F** is defined as [4]

$$
\mathcal{L}_{\mathbf{F}} V(x) := \{a \in \mathbb{R} \mid \exists \xi \in \mathbf{F}(x) \text{ such that} \\
\qquad \langle \eta, \xi \rangle = a, \text{ for all } \eta \in \partial V(x)\}. \quad (11)
$$

**Lemma 1.** [4] Let $x(\cdot)$ be a solution of (10) and let $V : \mathbb{R}^n \to \mathbb{R}$ be a convex function. Then $\frac{d}{dt}V(x(t))$ exists almost everywhere and $\frac{d}{dt}V(x(t)) \in \mathcal{L}_{\mathbf{F}}V(x(t))$ almost everywhere.

**Theorem 1.** [4] Let $V : \mathbb{R}^n \to \mathbb{R}$ be a convex positive-definite function. If, for all $x \in \mathbb{R}^n$, one has either $\max \mathcal{L}_{\mathbf{F}}V(x) \le 0$ or $\mathcal{L}_{\mathbf{F}}V(x) = \emptyset$, then (10) is Lyapunov stable at $x = 0$.

To simplify the exposition, we follow the convention
$\max \emptyset = -\infty$. A set $\Omega$ is said to be invariant with respect to (10) if for every initial condition $x_0 \in \Omega$, the
maximal solution $x(\cdot)$ of (10) lays in $\Omega$.

**Theorem 2.** [4] Let $V : \mathbb{R}^n \to \mathbb{R}$ be a function satisfying the conditions of Theorem 1 and let us assume that, for some $l > 0$, the set $\{\xi \in \mathbb{R}^n \mid V(\xi) \le l\}$ is bounded. Define the set

$$
\mathcal{Z}_{\mathbf{F}} V := \text{Zero} \mathcal{L}_{\mathbf{F}} V, \quad (12)
$$

let $x_0 \in L_1$, and let $\mathcal{M}$ be the largest invariant subset of $\mathcal{Z}_{\mathbf{F}}V \cap L_1$. Then,

$$
\mathrm{Dist}(x(t), \mathcal{M}) \to 0 \quad \text{as } t \to \infty. \qquad (13)
$$

3 Multivalued Hamiltonian systems
with dissipation

Let us consider a Hamiltonian system with generalized
coordinate $q \in \mathbb{R}^n$ and generalized momentum $p \in \mathbb{R}^n$.
For a given convex function $H: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$, $(q, p) \mapsto
H(q, p)$, the subdifferentials of $H$ with respect to the

first and second arguments are denoted as ∂_q H and ∂_p H,
respectively. We consider differential inclusions of the
form

$$
\begin{equation}
\begin{aligned}
\dot{q} &\in \partial_p H(q, p) - \partial_q g(q, p) \\
\dot{p} &\in -\partial_q H(q, p) - \partial_p g(q, p),
\end{aligned}
\tag{14}
\end{equation}
$$

where $g: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$ is also convex. Let us define
the multivalued Hamiltonian vector field $\mathbf{H}: \mathbb{R}^n \times \mathbb{R}^n \to
\mathbb{R}^n \times \mathbb{R}^n$ as

$$
\mathbf{H}(q, p) = \begin{pmatrix} \partial_p H(q, p) \\ -\partial_q H(q, p) \end{pmatrix}. \quad (15)
$$

Notice that, if $H$ is continuously differentiable and $g$ is a constant function, the multivalued Hamiltonian vector field reduces to a classic one. In analogy with single-valued Hamiltonian vector fields, multivalued Hamiltonian vector fields preserve the Hamiltonian function $H$. **Proposition 1.** *Consider a convex function $H : \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$. The function $H$ remains constant along the integral curves of the vector field $\mathbf{H}$.*

*Proof.* The set-valued derivative of $H$ along the trajectories of $\mathbf{H}$ is

$$
\begin{equation}
\begin{aligned}
\mathcal{L}_{\mathbf{H}} H(q, p) = {}& \{a \in \mathbb{R} \mid \exists (\xi_1, \xi_2) \in \partial H(q, p), \\
& \text{such that } \langle \eta_1, \xi_2 \rangle + \langle \eta_2, -\xi_1 \rangle = a, \\
& \text{for all } (\eta_1, \eta_2) \in \partial H(q, p) \}.
\end{aligned}
\tag{16}
\end{equation}
$$

The condition $\langle \eta_1, \xi_2 \rangle + \langle \eta_2, -\xi_1 \rangle = a$ must hold, in particular, for $\eta_1 = \xi_1$ and $\eta_2 = \xi_2$, in which case we have

$$
\langle \xi_1, \xi_2 \rangle + \langle \xi_2, -\xi_1 \rangle = a = 0, \quad (17)
$$

so either $\mathcal{L}_{\mathbf{H}}H(q,p) = \emptyset$ or $\mathcal{L}_{\mathbf{H}}H(q,p) = \{0\}$. By Lemma 1, $H$ remains constant along the integral curves of $\mathbf{H}$. $\square$

**Assumption 1.** The vector field $-\partial g$ in (14) is dissipative. That is,

$$
\max \mathcal{L}_{-\partial g} H(q, p) \le 0 \quad (18)
$$

for all $(q,p) \in \mathbb{R}^n \times \mathbb{R}^n$.

This assumption is common, e.g., in classical mechan-
ics. Indeed, in the Lagrangian and Hamiltonian frame-
works, the dissipative forces are typically included in the
form of *Rayleigh dissipation potential* functions. These
are typically convex functions depending on the second
argument only, *i.e.*, such that $\partial_q g = \{0\}$. The *twisting*
sliding-mode algorithm [23] also falls in this category, as
shown below.

The dissipativity of $-\partial g$ implies that $\max \mathcal{L}_{(\mathbf{H}-\partial g)}H \le 0$, so $H$ is a Lyapunov function if it is also positive definite. Asymptotic stability can then be established using the Invariance Principle formulated in Theorem 2.
## 2.2 Stability of Differential Inclusions

Consider a differential inclusion (DI)

$$ \dot{x}(t) \in \mathbf{F}(x(t)), \quad t \text{ - a.e.} \quad x(0) = x_0. \quad (2) $$

where a.e. means almost everywhere. We assume that the vector field $\mathbf{F}$ is such that, for each $x_0 \in \mathbb{R}^n$, there exists a (non-necessarily unique) absolutely continuous solution $x: I \to \mathbb{R}^n$, where $I$ is a real interval containing the origin. A solution is said to be maximal if there are no solutions $\bar{x}: \bar{I} \to \mathbb{R}^n$ such that $\bar{x}|_I = x$ with $I$ strictly contained in $\bar{I}$. To alleviate the notation, we will henceforth omit the time argument. Stability of DI's is surveyed in [12]. Recall that the set-valued derivative of a convex function $V: \mathbb{R}^n \to \mathbb{R}$ along $\mathbf{F}$ is defined as [3]

$$ \mathcal{L}_{\mathbf{F}}V(x) = \{a \in \mathbb{R} \mid \exists \xi \in \mathbf{F}(x) \text{ such that} \\ \langle \eta, \xi \rangle = a, \text{ for all } \eta \in \partial V(x)\}. $$

Whenever $V$ is $C^1$ and $\mathbf{F}$ single valued, the set-valued derivative reduces to the usual Lie derivative. In the sequel, we refer to $\mathcal{L}_{\mathbf{F}}V(x)$ as the derivative of $V$ along $\mathbf{F}$. Instead of specifying a vector field $\mathbf{F}$ in the derivative of $V$, we may specify a given DI. For example, we may write $\mathcal{L}_{(2)}V(x)$ in place of $\mathcal{L}_{\mathbf{F}}V(x)$.

**Lemma 1** [3] Let $x(\cdot)$ be a solution of (2) and let $V: \mathbb{R}^n \to \mathbb{R}$ be a convex function. Then $\frac{d}{dt}V(x(t))$ exists a.e. and $\frac{d}{dt}V(x(t)) \in \mathcal{L}_{\mathbf{F}}V(x(t))$ a.e..

**Theorem 1** [3] Let $V: \mathbb{R}^n \to \mathbb{R}$ be a convex positive-definite function. If, for all $x \in \mathbb{R}^n$, one has either $\max \mathcal{L}_{\mathbf{F}}V(x) \le 0$ or $\mathcal{L}_{\mathbf{F}}V(x) = \emptyset$, then (2) is Lyapunov stable at $x = 0$.

We follow the convention $\max \emptyset = -\infty$. A set $\Omega$ is said to be weakly invariant with respect to (2) if, for every initial condition $x_0 \in \Omega$, a maximal solution $x$ of (2) lies in $\Omega$.

**Theorem 2** [3] Let $V: \mathbb{R}^n \to \mathbb{R}$ be a function satisfying the conditions of Theorem 1 and let us assume that, for some $l > 0$, the set $L_l = \{\xi \in \mathbb{R}^n | V(\xi) \le l\}$ is bounded. Define the set

$$ Z_{\mathbf{F}}V = \text{Zero } \mathcal{L}_{\mathbf{F}}V, $$

let $x_0 \in L_l$, and let $\mathcal{M}$ be the largest weakly invariant subset of $Z_{\mathbf{F}}V \cap L_l$. Then, $\text{Dist}(x(t), \mathcal{M}) \to 0$ as $t \to \infty$.

## 2.3 Single-valued Hamiltonian systems

Locally, Hamiltonian systems have the canonical representation [15]

$$ \dot{x} = J \nabla H(x), \quad (3) $$

where

$$ x = \begin{pmatrix} q \\ p \end{pmatrix} \in \mathbb{R}^{2n} \quad \text{and} \quad J = \begin{pmatrix} 0 & I_n \\ -I_n & 0 \end{pmatrix}. \quad (4) $$

Here, $q \in \mathbb{R}^n$ is the generalised coordinate and $p \in \mathbb{R}^n$ is the generalised momentum. Note that the use of the gradient implicitly assumes that $H$ is at least $C^1$. An easily verifiable property is that $H$ remains constant along the system trajectories. Since for mechanical and electrical systems $H$ equals the energy stored in the system, the latter property amounts to energy conservation. In the context of passivity-based control, it has proved useful to generalise (3) to [33]

$$ \dot{x} = (J - R(x))\nabla H(x), \quad (5) $$

where $R(x) = R(x)^\top \succ 0$ is the damping matrix. In particular, the Interconnection and Damping Assignment technique [28] yields feedback laws such that the closed-loop system has the form (5).<fnref n="1" /> Note that the vector field in (5) is readily decomposed in two parts: $J\nabla H$ and $-R\nabla H$. By the skew-symmetry of $J$, we have $\mathcal{L}_{J\nabla H}H(x) = 0$, so we say that $J\nabla H$ is conservative. By the positive semi-definiteness of $R$, we have $\mathcal{L}_{-R\nabla H}H(x) \le 0$ and we will say that $-R\nabla H$ is dissipative. By the linearity of the derivative, we have

$$ \frac{d}{dt} H(x) = \mathcal{L}_{(J-R)\nabla H} H(x) \le 0. \quad (6) $$

The main interest of (6) is that, if $H$ is chosen positive definite, then it serves as a Lyapunov function and the stability of the origin is assured. Asymptotic stability is usually concluded at a later stage using LaSalle's Invariance Principle or Barbalat's Lemma.

# 3 Multivalued Hamiltonian systems

In this section, we extend the class of Hamiltonian systems by including functions $H$ that are convex but not necessarily differentiable. As a motivation, we will present several systems belonging to the extended class.

## 3.1 Autonomous systems with dissipation

An natural extension of (5) is

$$ \dot{x} \in (J - R(x))\partial H(x), \quad (7) $$

where $x$ and $J$ are as in (4), $R(x) = R(x)^\top \succ 0$, and $H$ is a convex function.

<footnote n="1">¹ Actually, it is possible to consider interconnection matrices $J$ that are state-depend, but we will not do so here.</footnote>
**Example 1 (Twisting Algorithm).** Consider a differential inclusion of the form (14), where $n = 1$,

$$H(q,p) = \gamma_1|q| + \frac{1}{2}p^2, \quad g(q,p) = \gamma_2|p|, \quad (19)$$

and

$$\gamma_1 > \gamma_2 > 0. \quad (20)$$

These parameters give the inclusion

$$\begin{aligned} \dot{q} &= p \\ \dot{p} &\in -\gamma_1 \operatorname{sgn}(q) - \gamma_2 \operatorname{sgn}(p), \end{aligned} \quad (21)$$

which is the twisting algorithm (we have used $\partial|\cdot| = \operatorname{sgn}(\cdot)$, the set-valued signum function). The Hamiltonian $H$ is positive definite and its set-valued derivative $H$ along (21) is

$$\mathcal{L}_{(\mathbf{H}-\partial g)} H(q,p) = \begin{cases} a \in \mathbb{R} & |\exists \xi_2 \in -\gamma_1 \operatorname{sgn}(q) - \gamma_2 \operatorname{sgn}(p) \\ \text{such that } \langle \eta_1, p \rangle + \langle p, \xi_2 \rangle = a, \\ \text{for all } \eta_1 \in \gamma_1 \operatorname{sgn}(q) \end{cases}. \quad (22)$$

Note that $\mathcal{L}_{(\mathbf{H}-\partial g)}H(q, 0) = 0$ for any $q \in \mathbb{R}$, that $\mathcal{L}_{(\mathbf{H}-\partial g)}H(q, p) = -\gamma_2|p|$ for any $p \in \mathbb{R}$ if $q \neq 0$, and that $\mathcal{L}_{(\mathbf{H}-\partial g)}H(q, p) = \emptyset$ otherwise. This can be compactly written as

$$\mathcal{L}_{(\mathbf{H}-\partial g)} H(q,p) = \begin{cases} \emptyset & \text{if } q = 0 \text{ and } p \neq 0 \\ -\gamma_2|p| & \text{otherwise} \end{cases}, \quad (23)$$

from which Lyapunov stability can be readily concluded. To establish asymptotic stability, we begin by computing

$$Z_{(\mathbf{H}-\partial g)} H = \text{Zero } \mathcal{L}_{(\mathbf{H}-\partial g)} H = \mathbb{R} \times \{0\}. \quad (24)$$

Clearly, any element $(q, p)$ of an invariant subset of $\mathbb{R} \times \{0\}$ must also be contained in the set

$$S = \{(q, p) \in \mathbb{R}^2 \mid 0 \in -\gamma_1 \operatorname{sgn}(q) - \gamma_2 \operatorname{sgn}(p)\}. \quad (25)$$

In other words, the set $\mathcal{M}$ in Theorem 2 is a subset of $(\mathbb{R} \times \{0\}) \cap S$. It can be readily seen that

$$(\mathbb{R} \times \{0\}) \cap S = \{(q, 0) \in \mathbb{R}^2 \mid 0 \in -\gamma_1 \operatorname{sgn}(q) + [-\gamma_2, \gamma_2]\}, \quad (26)$$

and that the condition (20) ensures that $(\mathbb{R} \times \{0\}) \cap S = \{(0, 0)\}$. This establishes the asymptotic stability of the origin. Finite-time convergence can be established by verifying that (28) is homogeneous with weights $(2, 1)$ and negative degree $-1$ [5, Cor. 5.4].

Before presenting the following example, allow us to define the signed bracket as $[x]^\gamma = |x|^\gamma \operatorname{sgn}(x)$, $\gamma \ge 0$, $x \in \mathbb{R}$.

**Example 2 (Super-Twisting Algorithm).** Consider a differential inclusion of the form (14), where $n = 1$,

$$H(q,p) = \gamma_2|q| + \frac{1}{2}p^2, \quad g(q,p) = \gamma_1\frac{2}{3}|q|^{3/2}, \quad (27)$$

$\gamma_1 > 0$, and $\gamma_2 > 0$. These parameters give the differential inclusion

$$\begin{aligned} \dot{q} &= -\gamma_1 \lfloor q \rfloor^{1/2} + p, \\ \dot{p} &\in -\gamma_2 \operatorname{sgn}(q) \end{aligned}, \quad (28)$$

which corresponds to the super-twisting algorithm [25,22,24]. The only minimum of $H$ is the origin. Following the same reasoning as in Example 1, we can see that the set-valued derivative of $H$ along (28) is

$$\mathcal{L}_{(\mathbf{H}-\partial g)} H(q,p) = \begin{cases} \emptyset & \text{if } q = 0 \text{ and } p \neq 0 \\ -\gamma_1|q|^{1/2} & \text{otherwise} \end{cases}, \quad (29)$$

from which Lyapunov stability follows. Asymptotic stability of the origin can be easily confirmed by noting that $Z_{(\mathbf{H}-\partial g)}H = \{(0, 0)\}$. Finite-time convergence can be established by verifying that (28) is homogeneous with weights $(2, 1)$ and negative degree $-1$ [5, Corollary 5.4].

**Example 3 (First-order sliding-mode control with actuator dynamics).** Consider a simple integrator

$$\dot{q} = p \quad (30a)$$

with $p$ a control input. It is well-known that the control law $p \in -\gamma \operatorname{sgn}(q)$ renders the system finite-time stable. Suppose that the actuator has parasitic dynamics, so that the control law cannot be enforced exactly. Instead, we have

$$\dot{p} \in -\frac{1}{\tau} (p + \gamma \operatorname{sgn}(q)), \quad 0 < \tau \ll 1. \quad (30b)$$

The differential inclusion (30) is of the form (14) with $H(q,p) = \frac{\gamma}{\tau}|q| + \frac{1}{2}p^2$ and $g(q,p) = \frac{1}{\tau}p^2$. Similar computations as those performed in Example 1 show that the system is asymptotically stable. However, the system is not homogeneous, so finite-time stability cannot be claimed. $\triangle$

**Example 4 (Mechanical system with Coulomb's friction).** Let us consider a multibody system with $n$ degrees of freedom and $m$ contact points with 3-dimensional friction. It is assumed that the Lagrange dynamics read as [8, Section 5.3.2]

$$\begin{aligned} M(q)\ddot{q} + C(q, \dot{q})\dot{q} + G(q) &= \sum_{i=1}^{m} H_{t,i}(q)\lambda_{t,i} + \nabla h_n(q)\lambda_n, \\ \lambda_{t,i} &\in -\mu_i |\lambda_{n,i}| \partial_{\dot{q}} \|H_{t,i}^{\top}(q)\dot{q}\|, \\ h_n(q) &= 0 \end{aligned}, \quad (31)$$
**Proposition 1** If $H: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$ is a convex positive-definite function, the origin of (7) is Lyapunov stable.

**PROOF.** The derivative of $H$ along the trajectories of (7) is

$$ \mathcal{L}_{(7)}H(x) = \{a \in \mathbb{R} \mid \exists \xi \in \partial H(x) \text{ such that} \\ \langle \eta, (J - R(x))\xi \rangle = a \text{ for all } \eta \in \partial H(x)\}. \quad (8) $$

The condition $\langle \eta, (J - R(x))\xi \rangle = a$ must hold, in particular, for $\eta = \xi$, in which case we have $\langle \xi, (J - R(x))\xi \rangle = a \le 0$ so $\max \mathcal{L}_{(7)}H(x) \le 0$. By Theorem 1, the origin is Lyapunov stable.

The main motivation for studying the stability of (7) is that some sliding-mode control systems have a mathematical structure that resembles that of mechanical systems, so that an energy function can be easily inferred from their structure. As with mechanical systems, energy functions may serve as Lyapunov functions. This is the case of the Lyapunov function presented in [27,12] for the twisting algorithm, and the one recalled in [24] for the super-twisting algorithm. However, such Lyapunov functions are weak, so additional considerations are required to conclude asymptotic stability.

**Example 1 (Super-Twisting Algorithm)** Consider a simple integrator $\dot{q} = u$ with $u$ a control input. The super-twisting scheme [20,18] is

$$ u = -\gamma_1 [q]^{1/2} + p, $$

where $\dot{p} \in -\gamma_2 \operatorname{sgn}(q)$ and $\gamma_1, \gamma_2 > 0$. This gives the undisturbed closed-loop

$$ \begin{aligned} \dot{q} &= -\gamma_1 [q]^{1/2} + p. \\ \dot{p} &= -\gamma_2 \operatorname{sgn}(q). \end{aligned} \quad (9) $$

It can be written as a Hamiltonian system of the form (7) with $n = 1$, and

$$ H(x) = \gamma_2|q| + \frac{1}{2}p^2, \quad R(x) = \begin{pmatrix} \frac{\gamma_1}{\gamma_2}|q|^{1/2} & 0 \\ 0 & 0 \end{pmatrix}. $$

We can conclude stability of the origin using Proposition 1. To establish asymptotic stability, we explicitly compute the derivative of $H$. For $q \ne 0$, $\partial H$ is single-valued, so we have

$$ \mathcal{L}_{(7)}H(x) = \gamma_2 \operatorname{sgn} q \left( -\gamma_1 [q]^{1/2} + p \right) \\ \qquad + p (-\gamma_2 \operatorname{sgn}(q)) = -\gamma_1 \gamma_2 |q|^{1/2}. $$

For $q = 0$ we have

$$ \mathcal{L}_{(7)}H(x) = \{a \in \mathbb{R} \mid \exists \xi_1 \in -\gamma_1 \operatorname{sgn}(q) \text{ such that} \\ \langle \eta_1, p \rangle + \langle p, \xi_1 \rangle = a \text{ for all } \eta_1 \in \gamma_1 \operatorname{sgn}(q)\}. $$

Note that $\mathcal{L}_{(7)}H(0) = \{0\}$ and $\mathcal{L}_{(7)}H(x) = \emptyset$ for $q = 0$, $p \ne 0$. This can be compactly written as

$$ \mathcal{L}_{(7)}H(x) = \begin{cases} \emptyset & \text{if } q = 0 \text{ and } p \ne 0 \\ \{-\gamma_1 \gamma_2 |q|^{1/2}\} & \text{otherwise} \end{cases}. $$

Asymptotic stability of the origin is confirmed by noting that $\mathcal{Z}_{(7)}H = \{(0,0)\}$. Finite-time convergence is established by verifying that (9) is homogeneous with weights (2, 1) and negative degree $-1$ [4, Corollary 5.4]. $\triangle$

**Example 2 (Dynamical Optimisation and Nesterov Gradient Acceleration)** Consider the system

$$ \ddot{z} + (1 - \beta)\dot{z} + \frac{1}{L}\nabla f(z + \beta\dot{z}) = 0, \quad L > 0. \quad (10) $$

It has been shown in [25] that, for $\beta = \frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$, where $\kappa$ is the condition number associated with $f$, a semi-implicit Euler discretisation of (10) leads to the celebrated accelerated gradient method of Nesterov. Notably, it is shown in [25] that, by taking the change of coordinates $q = z$, $p = \dot{z}$, system (10) can be written as a Hamiltonian system plus an extra term. Here, we consider the change of variables $q = z + \beta \dot{z}$, $p = \beta \dot{z}$, and the case where $f: \mathbb{R}^n \to \mathbb{R}$ is convex but not necessarily differentiable. In such case, system (10) is rewritten in the form (7) with

$$ \begin{aligned} H(x) &= \frac{\beta}{L} f(q) + \frac{\beta^2 - \beta + 1}{2\beta} \|p\|^2 \\ R(x) &= \begin{pmatrix} 1 & 0 \\ 0 & \frac{\beta(1-\beta)}{\beta^2 - \beta + 1} \end{pmatrix} \end{aligned} $$

$$ \text{and } 0 < \beta < 1. \quad \triangle $$

Another possible extension of (5) is

$$ \dot{x} \in J \partial H(x) + \mathbf{D}(x), \quad (11) $$

where $\mathbf{D}$ is a (possibly set-valued) dissipative vector field satisfying

$$ \max \mathcal{L}_{\mathbf{D}}H(x) \le 0 \quad (12) $$

for all $x \in \mathbb{R}^{2n}$.

**Example 3 (Mechanical Systems with Coulomb's Friction and Twisting Control)**

Consider an actuated planar rotating homogeneous disc with one frictional contact point (see Fig. 1). The disc
where $q \in C_{\text{conf}} \subseteq \mathbb{R}^n$, $C_{\text{conf}}$ is the configuration space of the system when all $m$ constraints are removed, $M(q) = M(q)^\top > 0$ is the inertia matrix, $C(q, \dot{q})$ contains Coriolis and centripetal torques, $G(q)$ is the potential force that derives from a smooth potential $P(q)$ (gravity, elasticity), $\mu_i > 0$ are the friction coefficients, $\lambda_{n,i}$ are the normal components of the contact forces at contact $i$, $H_{t,i}(q)$ and $h_n(q) : C_{\text{conf}} \to \mathbb{R}^m$ are defined from the local kinematics at contact points, such that the relative tangential velocity at contact $i$ satisfies $v_{t,i} = H_{t,i}^\top(q)\dot{q}$, and the relative normal velocity is $v_{n,i} = \nabla h_{n,i}^\top(q)\dot{q}$ [8, Chapters 4, 5]. Proceeding as usual with the generalized momentum $p = M(q)\dot{q}$, Equation (31) is rewritten in a Hamiltonian formalism:

$$ 
\begin{aligned}
\dot{q} &= \frac{\partial H}{\partial p}(q, p) \\
\dot{p} &\in -\frac{\partial H}{\partial q}(q, p) - \sum_{i=1}^{m} \mu_i |\lambda_{n,i}| H_{t,i}(q) \partial_{\dot{q}} \|H_{t,i}^{\top}(q) M^{-1}(q)p\| + \nabla h_n(q)\lambda_n \\
h_n(q) &= 0
\end{aligned}
\quad (32)
$$

where $H(q,p) = \frac{1}{2}p^\top M^{-1}(q)p + P(q)$. We may rewrite the dissipative contact force as

$$ 
\begin{aligned}
\mu_i |\lambda_{n,i}| M(q) M^{-1}(q) H_{t,i}(q) \partial_{\dot{q}} \| H_{t,i}^{\top}(q) M^{-1}(q) p \| \\
&= \mu_i |\lambda_{n,i}| M(q) \partial_p g_i(q, p),
\end{aligned}
\quad (33)
$$

with $g_i(q, \cdot) = \| \cdot \| \circ H_{t,i}^\top(q) M^{-1}(q)$, i.e., $g_i(q, p) = \|H_{t,i}^\top(q) M^{-1}(q)p\|$. Notice that $\dot{q}^\top \nabla h_n(q) \lambda_n = 0$, while $\dot{q}^\top M(q) \partial_p g_i(q, p) = p^\top \partial_p g_i(q, p) \ge 0$. Also, it is well-known that bilaterally constrained systems such as (31) may possess solutions with discontinuous velocities [8, section 5.6.4], [15,21]. Therefore (32) does not fit (14) in general. However in some cases it has absolutely continuous solutions and it fits with the above dissipative multivalued Hamiltonian framework. This is for instance the case for a planar disc sliding on a ground with persistent contact, or n stacked blocks with horizontal motion as in Fig. 4. In both cases $M(q) = M$ and $H_t(q) = H_t$ are constant matrices, and $\lambda_n = \lambda_n(t)$ is known, hence $g_i(q, p) = g_i(p)$. Thus, in a reduced set of coordinates Equation (32) boils down to

$$ 
\begin{aligned}
\dot{q} &= \frac{\partial H}{\partial p}(q, p) \\
\dot{p} &\in -\frac{\partial H}{\partial q}(q, p) - \sum_{i=1}^{m} \mu_i |\lambda_{n,i}(t)| M \partial g_i(p)
\end{aligned}
\quad (34)
$$

where the term $\nabla h_n(q) \lambda_n$ is no longer present, because it is irrelevant for the considered set of reduced coordinates. Notice that $q$ and $p$ in (34) are not the same as $q$ and $p$ in (32), however we keep the same notation (see the next examples).

For the homogeneous disc system with one frictional contact point we have $q = (x, \theta)^\top$ (the disc center has coordinates $(x, y)$, its rotation angle is $\theta$, and it is acted upon by a linear spring horizontally). The kinetic energy is $T(\dot{x}, \dot{\theta}) = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}I\dot{\theta}^2$ with $I = \frac{mr^2}{2}$, while the potential energy is $P(x) = \frac{1}{2}k(x-l)^2$ with $l > 0$ the unloaded linear spring's length, and $h_n(y) = y - r = 0$ (hence the vertical coordinate $y$ of the center of mass can be eliminated from the dynamics). Thus $p = (m\dot{x}, I\dot{\theta})^\top$, $H(q,p) = \frac{1}{2m}p_1^2 + \frac{1}{2I}p_2^2 + \frac{1}{2}k(x-l)^2$, $g_i(q,p) = g_i(p) = |r\dot{\theta} + \dot{x}| = |r p_2 + \frac{p_1}{m}|$, $i = 1, 2$, $\lambda_n = mg$, $H_{t,i} = (1, r)^\top$. The disc's dynamics is in the set-valued Hamiltonian formalism

$$ 
\begin{aligned}
\dot{x} &= \frac{1}{m} p_1 \\
\dot{\theta} &= \frac{1}{I} p_2 \\
\dot{p}_1 &\in -k(x-l) - \mu m^2 \mathbf{g} \partial_{p_1} \left| \left( \frac{1}{m}, \frac{r}{I} \right) p \right| \\
\dot{p}_2 &\in -\mu m \mathbf{g} I \partial_{p_2} \left| \left( \frac{1}{m}, \frac{r}{I} \right) p \right|.
\end{aligned}
\quad (35)
$$

Let us now examine the dynamics of n stacked blocks with mass $m_i$, $i = 1, \dots, n$, moving horizontally. The generalized coordinates and momentum are $q = (q_1, \dots, q_n)^\top$ and $p = (p_1, \dots, p_n)^\top$, respectively, with $p_i = m_i \dot{q}_i$. Let $H(q,p) = \frac{1}{2}p^\top M^{-1}p + P(q)$, where the potential energy is $P(q) = \frac{1}{2}(q-L)^\top K(q-L)$, $L = (l_1, \dots, l_n)^\top$, $l_i > 0$ are the springs' unloaded linear lengths, $K = K^\top \gtrsim 0$ is the stiffness matrix, and $M = \text{diag}(m_i)$. Coulomb friction is acting between each block and between the first block and the ground. Hence $\lambda_{n,i} = \sum_{j=i}^n m_j \mathbf{g}$. Then:

$$ 
\begin{aligned}
\dot{q} &= M^{-1}p \\
\dot{p} &= -K(q-L) + H_t \lambda_t \\
\lambda_{t,i} &\in -\mu_i \left( \sum_{j=i}^{n} m_j \right) \mathbf{g} \operatorname{sgn} \left( \frac{p_i}{m_i} - \frac{p_{i-1}}{m_{i-1}} \right)
\end{aligned}
\quad (36)
$$

where $p_0 = 0$ (the ground velocity is null),

$$ 
\begin{aligned}
H_{t,1} &= (1, 0, \dots, 0)^\top \\
H_{t,i} &= (0, \dots, 0, -1, 1, 0, \dots, 0)^\top \\
H_{t,n} &= (0, \dots, 0, -1, 1)^\top,
\end{aligned}
\quad (37)
$$

$i = 2, \dots, n-1$ and $H_t = (H_{t,1}, \dots, H_{t,n})$. Therefore,

$$ 
\begin{aligned}
\dot{p} &\in -K(q-L) \\
&\quad - \sum_{i=1}^{n} \mu_i \left( \sum_{j=i}^{n} m_j \right) \mathbf{g} M M^{-1} H_{t,i} \partial_i \left| H_{t,i}^{\top} M^{-1} p \right|.
\end{aligned}
\quad (38)
$$
centre has coordinates $(x, y)$ and rotation angle $\theta$. We choose the generalised coordinate vector $q = (q_1, q_2)^\top = (x, \theta)^\top$. The kinetic energy is $T(\dot{x}, \dot{\theta}) = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}I\dot{\theta}^2$ with $I = \frac{mr^2}{2}$. The disc is horizontally acted upon by a linear spring with stiffness $k > 0$, so the potential energy is $P(x) = \frac{1}{2}k(x - l)^2$ with $l > 0$ the unloaded linear spring's length. Thus, the generalised momentum is $p = (m\dot{x}, I\dot{\theta})^\top$ and the Hamiltonian is

$$H(q,p) = \frac{1}{2m}p_1^2 + \frac{1}{2I}p_2^2 + \frac{1}{2}k(q_1 - l)^2.$$

To take Coulomb's friction into account, let $g_1(p) = |r^{p_2/I} + \frac{p_1}{m}|$. The normal force at the contact point is $\lambda_n = mg$, so the disc's dynamics is given by

$$\begin{aligned}
\dot{q}_1 &= \frac{1}{m}p_1, \quad \dot{q}_2 = \frac{1}{I}p_2 \\
\dot{p}_1 &\in -k(q_1 - l) - \mu m^2 g \partial_{p_1} g_1(p) + u_1 \quad (13) \\
\dot{p}_2 &\in -\mu mgI \partial_{p_2} g_1(p) + u_2.
\end{aligned}$$

The control u can be used to modify the potential and the dissipation as

$$P(q) = \frac{1}{2}k(q_1 - l)^2 + \gamma_1|q_1| + \gamma_2|q_2|$$

and $g_{tw}(p) = \gamma_3|p_1| + \gamma_4|p_2|$, $\gamma_i > 0$, $i \in \{1, 4\}$, with the purpose of including a twisting control term. Then, the p-dynamics in (13) becomes

$$\begin{aligned}
\dot{p}_1 &\in -k(q_1 - l) - \gamma_1 \mathbf{sgn}(q_1) - \gamma_3 \mathbf{sgn}(p_1) \\
&\quad - \mu m^2 \mathbf{g} \partial_{p_1} g_1(p) \\
\dot{p}_2 &\in -\gamma_2 \mathbf{sgn}(q_2) - \gamma_4 \mathbf{sgn}(p_2) - \mu mgI \partial_{p_2} g_1(p) \quad (14)
\end{aligned}$$

The system is rewritten compactly as (11) with $H(q,p)$ as above, $\mathbf{D}(q,p) = -S\partial g_1(p) - \partial g_{tw}(p)$, with $S = \text{diag}(\mu m^2 \mathbf{g}, \mu mgI)$. Let us calculate the derivative

$$\begin{aligned}
\mathcal{L}_{(11)}H(q,p) = \{&a \in \mathbb{R} \mid \exists \xi \in (J\partial H + \mathbf{D})(p,q) \\
&\text{such that } \langle \eta, \xi \rangle = a, \\
&\text{for all } \eta \in (\gamma_1 \mathbf{sgn}(q_1) + k(q_1 - l), \gamma_2 \mathbf{sgn}(q_2), \frac{p_1}{m}, \frac{p_2}{I})^\top \} \\
= &\left\{ \emptyset, \text{ if } q_1 = 0 \text{ or } q_2 = 0; \right. \\
&\left. -\frac{\gamma_3}{m}|p_1| - \frac{\gamma_4}{I}|p_2| - \mu mgI \frac{r}{I} p_2 + \frac{p_1}{m}, \text{ otherwise.} \right\}
\end{aligned}$$

It is inferred that $\max \mathcal{L}_{(11)}H(q,0) = \{0\}$, and that $\max \mathcal{L}_{(11)}H(q,p) \le 0$ for all $p$ and $q$. Let us study now the set of zeros. We have

$$\mathcal{Z}_{(11)}H = \text{Zero } \mathcal{L}_{(11)}H = \mathbb{R}^2 \times \{0\} \times \{0\} \subset \mathbb{R}^4$$

if $\gamma_3 > 0$ or $\gamma_4 > 0$. If $\gamma_3 = \gamma_4 = 0$ then $\mathcal{L}_{J\partial H+\mathbf{D}}H = \mathbb{R}^2 \times \{p_1\} \times \{-\frac{I}{mr}p_1\}$, $p_1 \in \mathbb{R}$. Let $\gamma_3 > 0$ or $\gamma_4 > 0$ and let us compute the largest invariant set inside $\mathbb{R}^2 \times \{0\} \times$

$\{0\}$. Using (13) with $p_1 = p_2 = 0$ boils down to analysing the existence of $\xi_1 \in \mathbf{sgn}(q_1)$ and $\xi_2 \in \mathbf{sgn}(q_2)$ such that

$$k(q_1 - l) + \gamma_1 \xi_1 \in \gamma_3[-1, 1] + \mu mg[-1, 1] \quad (15a)$$

$$\gamma_2 \xi_2 \in \gamma_4[-1, 1] + \mu mg[-1, 1]. \quad (15b)$$

Assume that $\gamma_2 > \gamma_4 + \mu mg$, then $q_2 = 0$ is the only solution of (15b). If $k = 0$ and $\gamma_1 > \gamma_3 + \mu mg$, then $q_1 = 0$ is the only solution of (15a). If $k > 0$, $q_1 = 0$ is an equilibrium if and only if $kl \in (\gamma_1 + \gamma_3 + \mu mg)[-1, 1]$, $q_1 > 0$ is an equilibrium only if $\gamma_1 < \gamma_3 + \mu mg + kl$, $q_1 < 0$ is an equilibrium only if $\gamma_1 < \gamma_3 + \mu mg - kl$.

Let us note that other mechanical systems could be cast into such set-valued Hamiltonian dynamics, like the stacked blocks in Fig. 1. Following [6, section 5.5.2], this could be extended to Lagrange dynamics with set-valued friction at $m \ge 1$ contacts.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Left: Sliding disc. Right: stacked blocks.</figcaption></figure>

$\Delta$

In the following section, we present some examples of multivalued systems subject to perturbations.

## 3.2 Perturbed systems with dissipation

Our main interest is in the time discretisation of systems of the forms (7) and (11). In this paper, we are primarily focused on preserving properties like asymptotic or finite-time stability and leave robustness for future contributions. Although discrete-time perturbations will not be considered here, we will provide some examples that illustrate how continuous-time perturbations can be addressed within the proposed Hamiltonian framework. We will now consider perturbed versions of (7) and (11), that is,

$$\dot{x} \in (J - R(x))\partial H(x) + \mathbf{W}(x) \quad (16)$$

and

$$\dot{x} \in J\partial H(x) + \mathbf{D}(x) + \mathbf{W}(x). \quad (17)$$

Here, $\mathbf{W}$ is of the form $\mathbf{W}(x) = w(x) \cdot \mathcal{W}$, where $w$ is a single-valued vector field and the interval $\mathcal{W} = [-C, C]$ with $C \ge 0$ models (possibly time-varying) uncertainty.
That is,

$$ \dot{p} \in -K(q-L) - \sum_{i=1}^{n} \mu_i \left( \sum_{j=i}^{n} m_j \right) g M \partial g_i(p), \quad (39) $$

$$ \text{where } g_i(p) = (|\cdot| \circ H_{t,i}^{\top} M^{-1})(p) = |H_{t,i}^{\top} M^{-1} p|. \quad \triangle $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Left: Sliding disc. Right: stacked blocks.</figcaption></figure>

# 4 Time discretisation

A vast literature exists on the time discretisation of classical Hamiltonian systems, i.e., single-valued and without dissipation, see, e.g., [17] and references therein. It has been found that discrete-time models are well-posed and best retain the continuous-time properties when discretised using sophisticated combinations of backward-forward schemes [17]. In the multivalued scenario, however, well-posedness becomes critical, and many of such schemes are not implementable. In this section we show that, under reasonable assumptions, pure backward-Euler schemes for multivalued Hamiltonian systems are well-posed and inherit finite-time stability.

Motivated by Examples 1 and 2, we consider the time discretisation of (14) for two classes of dissipation potentials, one that depends on $p$ only and is analogous to mechanical friction, and one that depends on $q$ alone and has been disregarded in the literature. We will show that backward-Euler discretisation schemes preserve the dissipativity property of multivalued systems with dissipation, inasmuch as the Hamiltonian function is decreasing along the orbits of the discrete scheme. Thus, when positive definite, the Hamiltonian serves as a Lyapunov function for proving stability.

Finite-time convergence is proved using a simple argument that is particularly well suited for backward-Euler schemes.

**Proposition 2.** Consider the discrete-time system $x_{k+1} = T(x_k)$ and let

$$ \text{Fix } T := \{\xi \in \mathbb{R}^n \mid T(\xi) = \xi\} \quad (40) $$

be the set of fixed points of $T$. Suppose that Fix $T$ is a globally asymptotically stable set and suppose that there

exist $m \in \mathbb{N}$ and a uniform neighborhood $U$ of Fix $T$ such that

$$ U \subseteq T^{-m}(\text{Fix } T). \quad (41) $$

Then, Fix $T$ is globally finite-time stable.

**Proof.** If Fix $T$ is globally asymptotically stable, there exists a finite $k^* \in \mathbb{N}$ such that $x_k \in U$ for all $k \ge k^*$. Thus, the inclusion (41) implies that Fix $T$ is reached after $k^* + m$ steps. □

## 4.1 Momentum-dependent dissipation

Motivated by mechanical friction, we consider first the case in which dissipation depends on the momenta alone.

**Assumption 2.** The functions $H$ and $g$ in (14) satisfy the following:

(1) The Hamiltonian function is of the form $H(q,p) = P(q) + \frac{1}{2}\|p\|^2$ for some convex function $P: \mathbb{R}^n \to \mathbb{R}$.

(2) The function $g$ depends only on $p$, i.e., $g(q,p) = \tilde{g}(p)$ for some convex function $\tilde{g}$.

**Lemma 2.** Consider a set-valued Hamiltonian system (14) satisfying Assumption 2. For each current state $(q_k, p_k)$ there exists a unique future state

$$ (q_{k+1}, p_{k+1}) = T(q_k, p_k) \quad (42) $$

satisfying the backward-Euler discretisation scheme

$$ q_{k+1} \in q_k + h \partial_p H(q_{k+1}, p_{k+1}) \quad (43a) $$

$$ p_{k+1} \in p_k - h \partial_q H(q_{k+1}, p_{k+1}) - h \partial \tilde{g}(p_{k+1}). \quad (43b) $$

**Proof.** Taking into account Assumption 2, system (43) takes the particular form

$$ q_{k+1} = q_k + h p_{k+1} \quad (44a) $$

$$ p_{k+1} \in p_k - h \partial P(q_{k+1}) - h \partial \tilde{g}(p_{k+1}) \quad (44b) $$

The substitution of (44a) into (44b) yields

$$ p_{k+1} \in p_k - h \partial \tilde{G}_k(p_{k+1}), \quad (45) $$

where $\tilde{G}_k$ is the step-dependent map $p \mapsto \frac{1}{h}P(q_k + hp) + \tilde{g}(p)$. Since $\tilde{G}_k$ is maximal monotone for each $k$, the future momentum $p_{k+1}$ is uniquely determined by the expression

$$ p_{k+1} = \text{Prox}_{h \tilde{G}_k} (p_k). \quad (46) $$

The future position $q_{k+1}$ is now uniquely determined via (44a), and the conclusion follows. □

A usual strategy for showing the asymptotic stability of continuous-time dissipative Hamiltonian systems consists of three steps: i) Use $H$ as a Lyapunov function to establish stability, ii) use Barbalat's lemma to show that
### Example 4 (Nested Algorithm with Disturbances)

Consider the double integrator

$$ \begin{aligned} \dot{z}_1 &= z_2 \\ \dot{z}_2 &= \delta(t, z) + u \end{aligned} $$

where $z \in \mathbb{R}^2$ is the state, $u \in \mathbb{R}$ is the control and $\delta(t, z) \in \mathbb{R}$ is a time- and state-dependent perturbation. Assume that $|\delta(t, z)| \le C$ for all $z \in \mathbb{R}^2$, $t \ge 0$, and some $C \ge 0$. Following, e.g., [19], we can recast the perturbed system as the autonomous DI

$$ \begin{aligned} \dot{z}_1 &= z_2 \\ \dot{z}_2 &\in W + u \end{aligned} \qquad (18) $$

with $W = [-C, C]$. It is known that the second-order nested controller [19]

$$ u \in -\gamma \operatorname{sgn}(\lceil z_1 \rceil^{1/2} + z_2) $$

stabilises the origin in finite time for $\gamma$ high enough. The change of coordinates $q = z_1 + \lceil z_2 \rceil^2$ and $p = z_2$ takes the system to the form (16) with

$$ \begin{aligned} R(x) &= \begin{pmatrix} -2|p| & 0 \\ 0 & 0 \end{pmatrix}, & H(x) &= \gamma|q| + \frac{1}{2}p^2, \\ \text{and} & \mathbf{W}(x) = \left(2|p|, 1\right)^{\top} \mathbf{W}. \end{aligned} $$

For $q \ne 0$ we have

$$ \mathcal{L}_{(16)} H(x) = -|p| \left(2\gamma^2 - (2\gamma \operatorname{sgn} q + \operatorname{sgn} p) \cdot \mathbf{W}\right) $$

and $\max \mathcal{L}_{(16)} H(x) \le -|p| (2\gamma^2 - (2\gamma + 1)C)$. The inequality readily suggests that the gain satisfies

$$ \gamma = \frac{C + \sqrt{C^2 + 2C + 2\alpha}}{2}, \quad \alpha > 0, $$

in which case we have $\max \mathcal{L}_{(16)} H(x) \le -\alpha|p|$. For $q = 0$ and $\gamma \ge C + \frac{1}{2}$ we have $\mathcal{L}_{(16)} H(x) = \{-\frac{1}{2}|p|\}$. We have $\max \mathcal{L}_{(16)} H(x) \le -|p| \min \{\frac{1}{2}, \alpha\}$ whenever

$$ \gamma \ge \max \left\{ \frac{C + \sqrt{C^2 + 2C + 2\alpha}}{2}, C + \frac{1}{2} \right\}. \qquad (19) $$

Thus the origin is robustly stable. The largest invariant set contained in $\mathbb{R} \times \{0\}$ is the origin. This implies that stability is robustly asymptotic.

**Remark 1** While the finite-time stability of the nested controller has been established for sufficiently large $\gamma$, the explicit sufficient condition (19) is, to the best of our knowledge, new.

### Example 5 (Twisting Algorithm with Disturbances)

Consider the perturbed double integrator (18) and the twisting controller

$$ u \in -\gamma_1 \operatorname{sgn}(q) - \gamma_2 \operatorname{sgn}(p), \quad \gamma_1, \gamma_2 > 0. $$

It is a DI of the form (17) with $n = 1$,

$$ H(x) = \gamma_1|q| + \frac{1}{2}p^2, \quad \mathbf{D}(x) = -\partial g(q, p), $$

$g(q, p) = \gamma_2|p|$, and $w(x) = (0, 1)^\top$. We have $\mathcal{L}_{(17)} H(x) = \emptyset$ if $q = 0$ and $p \ne 0$ and $\mathcal{L}_{(17)} H(x) = -|p|(\gamma_2 + \mathcal{W})$ otherwise. Lyapunov stability can be concluded after setting $\gamma_2 > C$.

To establish asymptotic stability, we begin by computing $\mathcal{Z}_{(17)} H = \text{Zero } \mathcal{L}_{(17)} H = \mathbb{R} \times \{0\}$. Clearly, any element $x$ of an invariant subset of $\mathbb{R} \times \{0\}$ must also be contained in the set

$$ \mathcal{S} = \left\{ (q, p) \in \mathbb{R}^2 \mid 0 \in -\gamma_1 \operatorname{sgn}(q) - \gamma_2 \operatorname{sgn}(p) + \mathcal{W} \right\}. $$

In other words, the set $\mathcal{M}$ in Theorem 2 is a subset of $(\mathbb{R} \times \{0\}) \cap \mathcal{S}$. It can be seen that

$$ (\mathbb{R} \times \{0\}) \cap \mathcal{S} = \left\{ (q, 0) \in \mathbb{R}^2 \mid 0 \in -\gamma_1 \operatorname{sgn}(q) + \bar{\mathcal{W}} \right\} $$

with $\bar{\mathcal{W}} = [-\gamma_2 - C, \gamma_2 + C]$. The condition $\gamma_1 > \gamma_2 + C$ ensures that $(\mathbb{R} \times \{0\}) \cap \mathcal{S} = \{(0, 0)\}$. This establishes the asymptotic stability of the origin. Finite-time convergence can be concluded after verifying that the closed-loop is homogeneous with weights $(2, 1)$ and negative degree $-1$ [4, Corollary 5.4]. $\triangle$

Let us finish this section by emphasising some differences between the models (7) and (11), and pointing out some important caveats.

**Remark 2** In general,

$$ (J - R)\partial H(x) \subseteq J\partial H(x) - R\partial H(x). \qquad (20) $$

The practical significance of the previous remark is the following. Given that $\mathbf{D} = -R\partial H$ satisfies (12), a naive observation would be that (7) is a particular case of (11). However, this is not true in general, the main reason being that the inclusion in (20) may be strict. Thus, systems (7) and (11) need not be the same, even if $\mathbf{D} = -R\partial H$. Take, for simplicity, the unperturbed version
the dissipation function $\mathcal{L}_{\mathbf{H}-\partial\tilde{g}}H$ converges to zero, and iii) show (or assume) that $\mathcal{L}_{\mathbf{H}-\partial\tilde{g}}H = 0$ implies that the state converges to zero. We will follow a similar approach in the discrete-time setting. The following lemma will be useful in accomplishing the discrete-time counterparts of steps i) and ii).

**Lemma 3.** Consider a convex function $V: \mathbb{R}^n \to \mathbb{R}$ and a sequence $\{(v_k, w_k)\}_{k \in \mathbb{N}}$ such that $(v_k, w_k) \subset \text{gph} \partial V$ for all $k \in \mathbb{N}$. If $\lim_{k \to \infty} v_k = 0$ then

$$ \lim_{k \to \infty} \text{Dist}(\partial V(0), w_k) = 0. \quad (47) $$

*Proof.* Assume, for the sake of contradiction, that $\text{Dist}(\partial V(0), w_k) \to 0$ as $k \to \infty$. That is, that there exists $\varepsilon > 0$ such that, for any $k^* \in \mathbb{N}$, there is a $k > k^*$ such that $\text{Dist}(\partial V(0), w_k) > \varepsilon$. In particular, $w_k \notin \partial V(0)$ for some $k > k^*$. Thus, as $\partial V$ is outer semi-continuous at zero [11, Proposition 4.2.1-ii)], it follows that there exist neighborhoods of $w_k$ and zero (denoted, respectively, as $\mathcal{N}(w_k)$ and $\mathcal{N}(0)$) such that

$$ \mathcal{N}(0) \cap (\partial V)^{-1}(\mathcal{N}(w_k)) = \emptyset, \quad (48) $$

a contradiction, since $v_k \to 0$ as $k \to \infty$ and $v_k \in (\partial V^{-1})(\mathcal{N}(w_k))$. Therefore, (47) holds. □

The following assumption is used to establish the finite-time convergence to $(q, p) = (0, 0)$ for the scheme (44). It is a generalization of inequality (20).

**Assumption 3.** The subdifferentials $\partial P$ and $\partial \tilde{g}$ in Assumption 2 satisfy

$$ 0 \in -\partial \tilde{g}(0) \subset \text{int} \, \partial P(0), \quad (49) $$

where $\text{int} \, \partial P(0)$ is the interior of $\partial P(0)$. In addition, for any pair $(p, \xi) \in \text{gph} \partial \tilde{g}$ there exists $m \ge 1$ and $\alpha > 0$ such that

$$ \langle p, \xi \rangle \ge \alpha \|p\|^m. \quad (50) $$

Since $0 \in \partial \tilde{g}(0)$, condition (50) holds under diverse situations. For instance, if $\tilde{g}$ is $\mu$-strongly convex [18], then it holds with $m = 2$ and $\alpha = \mu$. Alternatively, if $\tilde{g}$ is convex and $0 \in \text{int} \, \partial \tilde{g}(0)$, then there is $\alpha > 0$ such that (50) holds with $m = 1$. Indeed, if $0 \in \text{int} \, \partial \tilde{g}(0)$, then there exists $\alpha > 0$, sufficiently small, such that $\alpha \mathcal{B} \in \partial \tilde{g}(0)$. Thus, it follows from the maximal monotonicity of $\partial \tilde{g}$ that, for any $(p, \xi) \in \text{gph} \partial \tilde{g}$,

$$ \langle p - 0, \xi - \alpha b \rangle \ge 0 \text{ for all } b \in \mathcal{B}. \quad (51) $$

Consequently,

$$ \langle p, \xi \rangle \ge \alpha \sup_{b \in \mathcal{B}} \langle b, p \rangle = \alpha \|p\|. \quad (52) $$

Also functions of the form $\tilde{g}(p) = \|p\|_r^{1+\beta}$ where $\beta > 0$ satisfy condition (50). Notice that in this case $\tilde{g}$ is neither strongly monotone for all cases of $\beta > 0$ nor $0 \in \text{int} \, \partial \tilde{g}(0)$, since $\tilde{g}$ is continuously differentiable. Nevertheless, the inequality (50) holds since, for any $(p, \xi) \in \text{gph} \partial \tilde{g}$,

$$ \begin{aligned} \langle p, \xi \rangle &= (1 + \beta) \|p\|_r^{\beta} \langle p, \tilde{\xi} \rangle \\ &= (1 + \beta) \|p\|_r^{\beta} (\|p\|_r + \Psi_{\mathcal{B}_s}(\xi)) \\ &= (1 + \beta) \|p\|_r^{1+\beta} \\ &\ge \alpha \|p\|^{1+\beta}, \end{aligned} \quad (53) $$

where we have used the fact that $\xi = (1 + \beta)\|p\|_r^{\beta}\tilde{\xi}$, $\tilde{\xi} \in \partial\|p\|_r$ in the first equality and in the second equality the fact that, for any $(p, \xi) \in \text{gph} \partial V$, $\langle p, \xi \rangle = V(p) + V^*(\xi)$, see e.g., [6, Proposition 16.9].

Note that, by maximal monotonicity, the condition $0 \in \text{int} \, \partial P(0)$ implies that

$$ \text{Zero} \, \partial P = \{0\}. \quad (54) $$

Indeed, for the sake of contradiction, let us assume there is $v \ne 0$ such that $0 \in \partial P(v)$. Thus, using Assumption 3, there exists $\varepsilon > 0$ such that $\varepsilon b \in \partial P(0)$ for all $b \in \mathcal{B}$. Now, the maximal monotonicity of $\partial P$ implies that

$$ 0 \le \langle 0 - \varepsilon b, v - 0 \rangle = -\varepsilon \langle b, v \rangle \text{ for all } b \in \mathcal{B}, \quad (55) $$

which leads to $v = 0$, a contradiction. Therefore, (54) holds. This in turn ensures that the origin is the only equilibrium and that $P$ is positive definite.

**Theorem 3.** Consider the backward-Euler discretisation (44) of an inclusion (14) satisfying Assumptions 2 and 3. Then, the origin is globally finite-time stable.

*Proof.* Take $H$ as a Lyapunov function candidate. By the definition of the convex subdifferential we have

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, q_{k+1} - q_k \rangle + \langle \eta_{2,k+1}, p_{k+1} - p_k \rangle \quad (56) $$

for any $\eta_{1,k} \in \partial P(q_k)$ and $\eta_{2,k} = p_k$. It follows from the discrete dynamics (44) that

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, h p_{k+1} \rangle + \langle p_{k+1}, -h \xi_{1,k+1} - h \xi_{2,k+1} \rangle, \quad (57) $$

where $\xi_{1,k} \in \partial P(q_k)$ and $\xi_{2,k} \in \partial \tilde{g}(p_k)$. Inequality (57) holds in particular for $\eta_{1,k} = \xi_{1,k}$, which implies that

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le -h \langle p_{k+1}, \xi_{2,k+1} \rangle. \quad (58) $$

Because of Assumption 3, we have

$$ \langle p_k, \xi_{2,k} \rangle \ge \gamma \alpha \|p_k\|^m, \quad (59) $$
($C=0$) of Example 4. Note that the vector fields

$$ (J - R)\partial H(x) = \begin{pmatrix} -2|p| & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \gamma \mathbf{sgn}(q) \\ p \end{pmatrix}, \\ J\partial H(x) = \begin{pmatrix} p \\ -\gamma \mathbf{sgn}(q) \end{pmatrix}, \quad R\partial H = \begin{pmatrix} 2\gamma|p| \mathbf{sgn}(q) \\ 0 \end{pmatrix} \tag{21} $$

satisfy (20) in a strict way. Moreover, at any point $q=0$, $p \neq 0$, we have $f \in J\partial H$ and $f \in R\partial H$ with $f(0,p) = (p,0)^\top$ (recall that $\gamma > \frac{1}{2}$). Thus, the trivial vector field $f-f$ is a selection of $J\partial H - R\partial H$ and we conclude that the whole $p$-axis is an equilibrium set. On the other hand, it has already been shown in Example 4 that the only equilibrium of $(J-R)\partial H$ is the origin.

**Remark 3** The set-valued derivative is sub-linear:

$$ \mathcal{L}_{F_a} H(x) + \mathcal{L}_{F_b} H(x) \subseteq \mathcal{L}_{F_a + F_b} H(x). \tag{22} $$

Since $\max \mathcal{L}_{J\partial H} H(x) = 0$, one maybe tempted to conclude that (12) is sufficient for asserting the stability of the origin of (11). This is not the case either, as the inclusion in (22) can also be strict. Thus, the condition (12) is only necessary. To conclude stability, we may impose the additional linearity condition

$$ \mathcal{L}_{J\partial H} H(x) + \mathcal{L}_{D} H(x) = \mathcal{L}_{J\partial H+D} H(x). \tag{23} $$

As an example, take again (21) and recall that $D = -R\partial H$. We have $\mathcal{L}_{J\partial H} H(x) + \mathcal{L}_{D} H(x) = \emptyset$ if $q=0$ and $p \neq 0$, and $\mathcal{L}_{J\partial H} H(x) + \mathcal{L}_{D} H(x) = \{-2\gamma^2|p|\}$ otherwise. On the other hand, we have

$$ \mathcal{L}_{J\partial H+D} H(x) = \begin{cases} [-\gamma|p|, \gamma|p|] & \text{if } q=0 \\ \{-2\gamma^2|p|\} & \text{otherwise.} \end{cases} $$

This problem does not occur with the twisting or super-twisting algorithms, as (23) is satisfied in those cases.

# 4 Time discretisation

Let us pass now to the analysis of time-discretisation methods of set-valued Hamiltonian systems like the ones described in the foregoing section. A vast literature exists on the time discretisation of classical Hamiltonian systems, i.e., single-valued and without dissipation, see, e.g., [15] and references therein. It has been found that discrete-time models are well-posed and best retain the continuous-time properties when discretised using sophisticated combinations of backward-forward schemes [15]. In the multivalued scenario, however, well-posedness becomes critical, and many of such schemes

are not implementable. In this section we show that, under reasonable assumptions, pure backward-Euler schemes for multivalued Hamiltonian systems are well-posed and inherit finite-time stability. In addition, such methods are very efficient in suppressing chattering at both input and output [8]. Because of Remarks 2 and 3, care has to be taken in regards with the type of chosen dissipation. Our choice is stated in the following.

**Assumption 1** The continuous-time system is of the form (11) with

$$ H(x) = P(q) + \frac{1}{2} \|p\|^2, \quad D(x) = -\partial g(q, p), $$

and with $P$ and $g$ convex functions.

Our motivation for the assumption is that Examples 1-3 and 5 can be written in the form (11) with $D(x) = -\partial g(q,p)$ (for Example 3 this is obtained after a coordinate change $\bar{q} = S^{-\frac{1}{2}}q$ and $\bar{p} = S^{-\frac{1}{2}}p$, for Example 1, we have $g(q,p) = g(q) = \frac{2\gamma_1}{3}|q|^{\frac{3}{2}}$). This type of dissipation is common in the Lagrangian and Hamiltonian frameworks, where the dissipative forces are typically included in the form of Rayleigh dissipation [7, Definition 6.12]. These are typically convex functions depending on $p$ only. All the examples have the Hamiltonian function

$$ H(x) = P(q) + \frac{1}{2} p^\top M p $$

with $M = M^\top > 0$ constant. Only Examples 3 and 2 have $M \neq I_n$, and the results that follow can be easily accommodated to such case. We show that backward-Euler discretisation schemes preserve the dissipativity property of multivalued systems with dissipation, inasmuch as the Hamiltonian is decreasing along the orbits of the discrete scheme. Thus, when positive definite, the Hamiltonian serves as a Lyapunov function for stability.

Note that Assumption 1 leads to the backward-Euler method

$$ \begin{aligned} q_{k+1} &= q_k + h p_{k+1} - h \partial_q g(q_{k+1}, p_{k+1}) \\ p_{k+1} &= p_k - h \partial P(q_{k+1}) - h \partial_p g(q_{k+1}, p_{k+1}). \end{aligned} \tag{24} $$

It is not clear in general if, for any $(q_k, p_k)$, the future state $(q_{k+1}, p_{k+1})$ can be uniquely found. Except for Example 2, the potential function $g$ depends either on $q$ or on $p$ alone. Thus, we consider the time discretisation for two dissipation classes and show that the method (24) is well-defined under mild additional conditions.

Finite-time convergence is proved using a simple argument that is particularly well suited for backward-Euler schemes.

**Proposition 2** Consider the discrete-time map $x_{k+1} = T(x_k)$ and let $\text{Fix} T = \{\xi \in \mathbb{R}^n | T(\xi) = \xi\}$ be the set
for all $k \in \mathbb{N}$ and for some $\gamma > 0$, $\alpha > 0$, $m \ge 1$, so that $H(q_{k+1}, p_{k+1}) \le H(q_k, p_k)$. This establishes the stability of the origin. By summing both sides of (58) we obtain

$$ 
\begin{aligned}
& \sum_{k=0}^{N} (H(q_{k+1}, p_{k+1}) - H(q_k, p_k)) \\
&= H(q_{N+1}, p_{N+1}) - H(q_0, p_0) \\
&\le -h \sum_{k=0}^{N} \langle p_{k+1}, \xi_{2,k+1} \rangle. \quad (60)
\end{aligned}
$$

Since the sequence $\{H(q_N, p_N)\}_{N \in \mathbb{N}}$ is decreasing and bounded from below, we conclude that it converges to a limit $\bar{H} \ge \min_{q,p} H(q,p)$. Thus,

$$ h \sum_{k=0}^{\infty} \langle p_{k+1}, \xi_{2,k+1} \rangle \le H(q_0, p_0) - \bar{H}. \quad (61) $$

The fact that $\langle p_k, \xi_{2,k} \rangle$ is non-negative implies that $\lim_{k \to \infty} \langle p_k, \xi_{2,k} \rangle = 0$, while the bound (59) ensures that $\lim_{k \to \infty} p_k = 0$. It follows from Lemma 3 that $\lim_{k \to \infty} \text{Dist}(\partial \tilde{g}(0), \xi_{2,k}) = 0$ but, according to (44b), we have

$$ \xi_{1,k} + \xi_{2,k} \to 0 \quad \text{as} \quad k \to \infty, \quad (62) $$

so $\lim_{k \to \infty} \text{Dist}(\partial \tilde{g}(0), -\xi_{1,k}) = 0$ and, by (49), there exist $\varepsilon > 0$ and $k^* \in \mathbb{N}$ such that, for all $k \ge k^*$, $\xi_{1,k} + b_k \in \partial P(0)$ for any $b_k \in \varepsilon\mathcal{B}$. From $\xi_{1,k} \in \partial P(q_{k+1})$ and the monotonicity property of subdifferentials we know that

$$ \langle b_k, q_k \rangle \ge 0, \quad (63) $$

for any $b_k \in \varepsilon\mathcal{B}$ and all $k \ge k^*$. Therefore, $q_k = 0$ for all $k \ge k^*$, and hence $p_{k+1} = 0$ for all $k \ge k^*$ in view of (44). $\square$

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

**Example 5.** Examples 1 and 3 satisfy Assumptions 2 and 3. Consider Example 4. It satisfies Assumption 2 if $M = M^\top > 0$ is constant (at the price of changing the norm accordingly). The linear elasticity potential $P(q)$ considered both in (35) and (36) does not satisfy (49) (in both cases it is possible to change $q$ to $q' := q - L$). It is however possible to shape the potential by introducing a set-valued controller of the form $u = -G \partial \sum_{j=1}^n |q_j|$, $G \in \mathbb{R}^{n \times n}$ a suitable control gain (thus the disc has one input for $x$ and one for $\theta$, while each mass is controlled in the system of stacked blocks). The potential $P(q)$ is changed to $P(q) + G \sum_{j=1}^n |q_j|$. Then (50) can be satisfied also. The obtained closed-loop dynamics is similar to a twisting scheme and may be seen as an extension of [1].

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

## 4.2 Douglas-Rachford-like Splitting Discretisation

Lemma 2 and Theorem 3 show that, when the dissipation function depends only on the momenta $p$, the

discrete-time system (44) is well-posed and finite-time stable. However, an explicit expression for the proximal map (46) is, in general, hard to obtain. Nevertheless, approximations can be computed via splitting schemes as follows. Let us start by considering the following system:

$$ q_{k+1} = q_k + hp_{k+1} \quad (64a) $$

$$ p_{k+1} = p_k - h\xi_{1,k+1} - h\xi_{2,k+1} \quad (64b) $$

$$ \xi_{1,k+1} \in \partial P(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \quad (64c) $$

$$ \xi_{2,k+1} \in \partial \tilde{g}(p_{k+1}). \quad (64d) $$

It can be seen as a dynamic extension of (44). Indeed, the particular value of the $\xi_{2,k+1} \in \partial \tilde{g}(p_{k+1})$ now depends implicitly on $\xi_{2,k}$, effectively making $\xi_2$ a new state variable. By the outer semicontinuity of $\partial P$ we know that, for any $\varepsilon > 0$, there exists $h > 0$ sufficiently small such that

$$ \partial P(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \subset \partial P(q_{k+1}) + \varepsilon \mathcal{B}. \quad (65) $$

If $\partial \tilde{g}$ is uniformly bounded, then for all $k \in \mathbb{N}$, $\|\xi_{2,k}\| \le M$ for some $M > 0$. Thus, $q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k}) \in \{q_{k+1}\} + 2h^2M \cdot \mathcal{B} = \{q_{k+1}\} + O(h^2) \cdot \mathcal{B}$. In this sense, the scheme (64) approximates (44). We now address its well-posedness.

**Lemma 4.** Consider a difference equation (64) satisfying Assumption 2. For each current state $(q_k, p_k, \xi_{2,k})$ there exists a unique future state

$$ (q_{k+1}, p_{k+1}, \xi_{2,k+1}) = T(q_k, p_k, \xi_{2,k}) \quad (66) $$

Furthermore, the selections $\xi_{1,k+1}$, $\xi_{2,k+1}$ are determined as

$$ \xi_{1,k+1} = \frac{1}{h^2} (\mathrm{Id} - \mathrm{Prox}_{h^2 P}) (q_k + hp_k - h^2\xi_{2,k}) \quad (67a) $$

$$ \xi_{2,k+1} = \frac{1}{h} (\mathrm{Id} - \mathrm{Prox}_{h\tilde{g}}) (p_k - h\xi_{1,k+1}). \quad (67b) $$

*Proof.* Define $\theta_{k+1} := q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})$. It follows from (64a)-(64c) that

$$ 
\begin{aligned}
\theta_{k+1} + h^2 \xi_{1,k+1} &= q_k + hp_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k}) + h^2 \xi_{1,k+1} \\
&= q_k + hp_k - h^2 \xi_{2,k}
\end{aligned}
\quad (68) $$

Since $\xi_{1,k+1} \in \partial P(\theta_{k+1})$, we have

$$ \theta_{k+1} = \mathrm{Prox}_{h^2 P} (q_k + hp_k - h^2 \xi_{2,k}). \quad (69) $$

Substituting (69) into (68) and solving for $\xi_{1,k+1}$ gives (67a), while (64b) gives

$$ p_{k+1} + h\xi_{2,k+1} = p_k - h\xi_{1,k+1}, \quad (70) $$

Hence, (64d) implies that

$$ p_{k+1} = \mathrm{Prox}_{h\tilde{g}}(p_k - h\xi_{1,k+1}). \quad (71) $$
of fixed points of T. Suppose that Fix T is a globally asymptotically stable set and suppose that there exist $m \in \mathbb{N}$ and a uniform neighbourhood U of Fix T such that

$$U \subseteq T^{-m}(\text{Fix } T). \qquad (25)$$

Then, Fix T is globally finite-time stable.

**PROOF.** If Fix T is globally asymptotically stable, there exists a finite $k^* \in \mathbb{N}$ such that $x_k \in U$ for all $k \ge k^*$. Thus, the inclusion (25) implies that Fix T is reached after $k^* + m$ steps. ■

## 4.1 Momentum-dependent dissipation

Motivated by Coulomb's friction, we consider first the case in which dissipation depends on the momenta alone. In the sequel, $h > 0$ denotes the timestep (or the sampling period), and $f_k = f(t_k)$, $t_k = kh$, $k \in \mathbb{N}$.

**Assumption 2** The system satisfies Assumption 1 with $g$ depending only on $p$, i.e., $g(q, p) = \tilde{g}(p)$ for some convex function $\tilde{g}$.

**Lemma 2** Consider a Hamiltonian system satisfying Assumption 2. For each current state $(q_k, p_k)$ there exists a unique future state

$$(q_{k+1}, p_{k+1}) = T(q_k, p_k) \qquad (26)$$

satisfying the backward-Euler discretisation scheme

$$q_{k+1} \in q_k + h\partial_p H(q_{k+1}, p_{k+1}) \qquad (27a)$$

$$p_{k+1} \in p_k - h\partial_q H(q_{k+1}, p_{k+1}) - h\partial\tilde{g}(p_{k+1}). \qquad (27b)$$

**PROOF.** Taking into account the assumptions on the system, it takes the particular form

$$q_{k+1} = q_k + hp_{k+1} \qquad (28a)$$

$$p_{k+1} \in p_k - h\partial P(q_{k+1}) - h\partial\tilde{g}(p_{k+1}). \qquad (28b)$$

The substitution of (28a) into (28b) yields $p_{k+1} \in p_k - h\partial\tilde{G}_k(p_{k+1})$, where $\tilde{G}_k$ is the step-dependent map $p \mapsto \frac{1}{h}P(q_k + hp) + \tilde{g}(p)$. Since $\tilde{G}_k$ is maximal monotone for each $k$, the future momentum $p_{k+1}$ is uniquely determined by the expression

$$p_{k+1} = \text{Prox}_{h\tilde{G}_k}(p_k). \qquad (29)$$

The future position $q_{k+1}$ is now uniquely determined via (28a) and the conclusion follows. ■

As mentioned above, a usual strategy for showing the asymptotic stability of continuous-time dissipative

Hamiltonian systems consists of three steps: i) Use $H$ as a Lyapunov function to establish stability, ii) use Barbalat's lemma to show that the dissipation function $\mathcal{L}_{H-\partial\tilde{g}}H$ converges to zero, and iii) show (or assume) that $\mathcal{L}_{H-\partial\tilde{g}}H = 0$ implies that the state converges to zero. We will follow a similar approach in the discrete-time setting. The next lemma is useful in accomplishing the discrete-time counterparts of steps i) and ii).

**Lemma 3** Consider a convex function $V: \mathbb{R}^n \to \mathbb{R}$ and a sequence $\{(v_k, w_k)\}_{k \in \mathbb{N}}$ such that $(v_k, w_k) \subset \text{gph} \partial V$ for all $k \in \mathbb{N}$. If $\lim_{k \to \infty} v_k = 0$ then

$$\lim_{k \to \infty} \text{Dist}(w_k, \partial V(0)) = 0. \qquad (30)$$

**PROOF.** Assume, for the sake of contradiction, that $\text{Dist}(w_k, \partial V(0)) \to 0$ as $k \to \infty$. That is, that there exists $\varepsilon > 0$ such that, for any $k^* \in \mathbb{N}$, there is a $k > k^*$ such that $\text{Dist}(w_k, \partial V(0)) > \varepsilon$. In particular, $w_k \notin \partial V(0)$ for some $k > k^*$. However, as $\partial V$ is osc at zero [10, Proposition 4.2.1-ii)], it follows that there exist neighbourhoods of $w_k$ and zero (denoted, respectively, as $\mathcal{N}(w_k)$ and $\mathcal{N}(0)$) such that $\mathcal{N}(0) \cap (\partial V)^{-1}(\mathcal{N}(w_k)) = \emptyset$, a contradiction, since $v_k \to 0$ as $k \to \infty$ and $v_k \in (\partial V^{-1})(\mathcal{N}(w_k))$. Therefore, (30) holds. ■

The following assumption is used to establish the finite-time convergence to $(q, p) = (0, 0)$ for the scheme (28). It is a generalisation of inequality $\gamma_1 > \gamma_2 > 0$ required for the unperturbed twisting controller.

**Assumption 3** The subdifferentials $\partial P$ and $\partial\tilde{g}$ in Assumption 2 satisfy

$$0 \in -\partial\tilde{g}(0) \subset \text{int} \, \partial P(0), \qquad (31)$$

where $\text{int} \, \partial P(0)$ is the interior of $\partial P(0)$. In addition, for any pair $(p, \xi) \in \text{gph} \, \partial\tilde{g}$ there exists $m \ge 1$ and $\alpha > 0$ such that

$$\langle p, \xi \rangle \ge \alpha \|p\|^m. \qquad (32)$$

Condition (31) implies the set-valuedness of $\partial P$ at 0. Since $0 \in \partial\tilde{g}(0)$, condition (32) holds under diverse situations. For instance, if $\tilde{g}$ is $\mu$-strongly convex [16], then it holds with $m = 2$ and $\alpha = \mu$. Alternatively, if $\tilde{g}$ is convex and $0 \in \text{int} \, \partial\tilde{g}(0)$, then there is $\alpha > 0$ such that (32) holds with $m = 1$. Indeed, if $0 \in \text{int} \, \partial\tilde{g}(0)$, then there exists $\alpha > 0$, sufficiently small, such that $\alpha \in \partial\tilde{g}(0)$. Thus, it follows from the maximal monotonicity of $\partial\tilde{g}$ that, for any $(p, \xi) \in \text{gph} \, \partial\tilde{g}$, $\langle p-0, \xi-\alpha b \rangle \ge 0$ for all $b \in \mathcal{B}$. Consequently, $\langle p, \xi \rangle \ge \alpha \sup_{b \in \mathcal{B}} \langle b, p \rangle = \alpha \|p\|$. Functions $\tilde{g}(p) = \|p\|_{\mathcal{B}}^{1+\beta}$ where $\beta > 0$ satisfy condition (32). Notice that in this case neither is $\tilde{g}$ strongly monotone for all $\beta > 0$ nor $0 \in \text{int} \, \partial\tilde{g}(0)$, since $\tilde{g}$ is $\mathcal{C}^1$. Nevertheless,
Substituting (71) into (70) and solving for $\xi_{2,k+1}$ gives (67b). Thus, $p_{k+1}$, and hence $q_{k+1}$, depend on data available at time $t_k$ so that (64) is well-defined. This concludes the proof. □

**Lemma 4** shows the explicit expression for the selection values $\xi_{1,k+1}$ and $\xi_{2,k+1}$. Notice that $\xi_{2,k+1}$ depends on its past values, so that the full state of the proposed approximation (64) lies in $\mathbb{R}^{3n}$.

**Remark 1.** As introduced in [30] for the case of multi-variable super-twisting-like schemes, the approximation of the selection values in (67) is reminiscent of the celebrated Douglas-Rachford splitting [26, ?]. Indeed, set

$$ \vartheta_{k+1} = p_{k+1} + h\xi_{2,k+1} \quad (72) $$

and note that it follows from (64d) that $p_{k+1} = \text{Prox}_{h\tilde{g}}(\vartheta_{k+1})$, so that after simple computations we arrive at the expressions

$$ \vartheta_{k+1} = \frac{1}{h} \left( \text{Prox}_{h^2 P} \left( q_k + h(2 \text{Prox}_{h\tilde{g}} - \text{I}_d)(\vartheta_k) \right) - q_k \right) + (\text{I}_d - \text{Prox}_{h\tilde{g}})(\vartheta_k) \quad (73a) $$

$$ p_{k+1} = \text{Prox}_{h\tilde{g}}(\vartheta_{k+1}), \quad (73b) $$

which define the Douglas-Rachford splitting associated with the generalized equation

$$ 0 \in h\partial P(q + hp) + h\partial\tilde{g}(p), \quad (74) $$

see e.g., [35, Section 2.7].

In what follows the stability properties of the set

$$ \Gamma := \left\{ \begin{bmatrix} q \\ p \\ \xi_2 \end{bmatrix} \in \mathbb{R}^{3n} \mid q = p = 0, \xi_2 \in \mathbb{R}^n \right\} \quad (75) $$

are studied.

**Theorem 4.** Consider a discrete-time system (64) satisfying Assumption 2. Strengthen Assumption 3 to

$$ \mu\mathcal{B} \subset -\partial\tilde{g}(0) \subset \text{int} \, \partial P(0) \quad (76) $$

for some $\mu > 0$. In addition, assume that $\partial\tilde{g}$ is uniformly bounded. Then, there exists $h^* > 0$, sufficiently small, such that for all $0 < h < h^*$ the set $\Gamma$ in (75) is globally asymptotically stable.

**Proof.** Let us consider the function $H$ introduced in Assumption 2. Using again the definition of convex subdifferential, we can readily show that

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, q_{k+1} - q_k \rangle + \langle \eta_{2,k+1}, p_{k+1} - p_k \rangle \quad (77) $$

for any $\eta_{1,k} \in \partial P(q_k)$ and $\eta_{2,k} = p_k$. It follows from (64c), (65), and the outer semicontinuity of $\partial P$ that, for any $0 < \varepsilon_1 < \mu$, there exists $h > 0$ sufficiently small such that

$$ \xi_{1,k+1} \in \partial P(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \subset \partial P(q_{k+1}) + \varepsilon_1 \mathcal{B}. \quad (78) $$

It is inferred that, for any $k \in \mathbb{N}$, there always exists $b_k \in \mathcal{B}$ such that

$$ \xi_{1,k} + \varepsilon_1 b_k \in \partial P(q_k). \quad (79) $$

Hence, taking $\eta_{1,k} = \xi_{1,k} + \varepsilon_1 b_k$ and $\eta_{2,k} = p_k$ in (77) leads us to

$$ \begin{aligned} H(q_{k+1}, p_{k+1}) - H(q_k, p_k) &\le \langle \xi_{1,k+1} + \varepsilon_1 b_{k+1}, q_{k+1} - q_k \rangle \\ &\quad + \langle p_{k+1}, p_{k+1} - p_k \rangle \\ &= h \langle \xi_{1,k+1} + \varepsilon_1 b_{k+1}, p_{k+1} \rangle - h \langle p_{k+1}, \xi_{1,k+1} + \xi_{2,k+1} \rangle \\ &\le h \varepsilon_1 \|p_{k+1}\| - h \langle p_{k+1}, \xi_{2,k+1} \rangle. \end{aligned} \quad (80) $$

It follows from (76) that

$$ \langle p_k, \xi_{2,k} \rangle > \mu \|p_k\|. \quad (81) $$

To see this, note that $-\mu b_k \in \partial\tilde{g}(0)$ for any $b_k \in \mathcal{B}$ and that, by the maximal monotonicity of $\partial\tilde{g}$, we have

$$ \langle p_k, \xi_{2,k} - \mu b_k \rangle \ge 0. \quad (82) $$

In other words, $\langle p_k, \xi_{2,k} \rangle \ge \mu \langle p_k, b_k \rangle$ for all $b_k \in \mathcal{B}$. Therefore,

$$ \langle p_k, \xi_{2,k} \rangle \ge \mu \sup_{b_k \in \mathcal{B}} \langle p_k, b_k \rangle, \quad (83) $$

which implies (81). Hence, (80) can be further bounded as

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le -h(\mu - \varepsilon_1) \|p_{k+1}\| \quad (84) $$

Since by assumption $0 < \varepsilon_1 < \mu$, it follows that $h(\mu - \varepsilon_1)\|p_k\| \ge 0$ for all $k \in \mathbb{N}$ and stability of the set $\Gamma$ is concluded. To show asymptotic convergence, we apply the same arguments as in the proof of Theorem 3, leading us to

$$ \lim_{k \to \infty} p_k = 0 \quad \text{and} \quad \lim_{k \to \infty} \text{Dist}(\partial\tilde{g}(0), -\xi_{1,k}) = 0. \quad (85) $$

By following the same steps as in the proof of Theorem 3, we can conclude from (64c), (76), and (85) that there exists $k^* \in \mathbb{N}$ such that for all $k > k^*$,

$$ q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k}) = 0. \quad (86) $$

Now, the substitution of (64a) into (86) yields

$$ \liminf_{k \to \infty} (q_k + h^2 \xi_{2,k+1}) = \liminf_{k \to \infty} (-hp_{k+1} + h^2 \xi_{2,k}) \quad (87) $$
the inequality (32) holds since, for any $(p, \xi) \in \text{gph} \partial \tilde{g}$,

$$ \langle p, \xi \rangle = (1+\beta)\|p\|_r^\beta \langle p, \tilde{\xi} \rangle = (1+\beta)\|p\|_r^\beta (\|p\|_r + \Psi_{\mathcal{B}_s}(\xi)) \\ = (1+\beta)\|p\|_r^{1+\beta} \ge \alpha \|p\|^{1+\beta}, $$

where we have used $\xi = (1+\beta)\|p\|_r^\beta \tilde{\xi}$, $\tilde{\xi} \in \partial \|p\|_r$, in the first equality and the fact that, for any $(p, \xi) \in \text{gph} \partial V$, $\langle p, \xi \rangle = V(p) + V^*(\xi)$ [5, Proposition 16.9] in the second one. Note that, by maximal monotonicity, the condition $0 \in \text{int} \partial P(0)$ implies that

$$ \text{Zero} \partial P = \{0\}. \qquad (33) $$

Indeed, for the sake of contradiction, let us assume there is $v \ne 0$ such that $0 \in \partial P(v)$. Using Assumption 3, there exists $\varepsilon > 0$ such that $\varepsilon b \in \partial P(0)$ for all $b \in \mathcal{B}$. Now, maximal monotonicity of $\partial P$ implies that $0 \le \langle 0 - \varepsilon b, v - 0 \rangle = -\varepsilon \langle b, v \rangle$ for all $b \in \mathcal{B}$, which leads to $v = 0$, a contradiction. Therefore, (33) holds. This in turn ensures that the origin is the only equilibrium and that $P$ is positive definite.

**Theorem 3** Consider the backward-Euler discretisation (28) of an inclusion satisfying Assumptions 2 and 3. The origin is globally finite-time stable.

**PROOF.** Take $H$ as a candidate Lyapunov function. By definition of the convex subdifferential we have

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, q_{k+1} - q_k \rangle + \langle \eta_{2,k+1}, p_{k+1} - p_k \rangle $$

for any $\eta_{1,k+1} \in \partial P(q_{k+1})$ and $\eta_{2,k+1} = p_{k+1}$. It follows from (28) that

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, h p_{k+1} \rangle + \langle p_{k+1}, -h \xi_{1,k+1} - h \xi_{2,k+1} \rangle, \quad (34) $$

where $\xi_{1,k+1} \in \partial P(q_{k+1})$ and $\xi_{2,k+1} \in \partial \tilde{g}(p_{k+1})$. Inequality (34) holds in particular for $\eta_{1,k+1} = \xi_{1,k+1}$, which implies that

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le -h \langle p_{k+1}, \xi_{2,k+1} \rangle. \quad (35) $$

From Assumption 3, we have

$$ \langle p_k, \xi_{2,k} \rangle \ge \gamma \alpha \|p_k\|^m \quad (36) $$

for all $k \in \mathbb{N}$ and for some $\gamma > 0$, $\alpha > 0$, $m \ge 1$, so that $H(q_{k+1}, p_{k+1}) \le H(q_k, p_k)$. This establishes the

stability of the origin. Summing both sides of (35) yields

$$ \sum_{k=0}^{N} (H(q_{k+1}, p_{k+1}) - H(q_k, p_k)) = H(q_{N+1}, p_{N+1}) - H(q_0, p_0) \le -h \sum_{k=0}^{N} \langle p_{k+1}, \xi_{2,k+1} \rangle. $$

Since the sequence $\{H(q_N, p_N)\}_{N \in \mathbb{N}}$ is decreasing and bounded from below, we conclude that it converges to a limit $\bar{H} \ge \min_{q,p} H(q,p)$. Thus,

$$ h \sum_{k=0}^{\infty} \langle p_{k+1}, \xi_{2,k+1} \rangle \le H(q_0, p_0) - \bar{H}. $$

Inequality $\langle p_{k+1}, \xi_{2,k+1} \rangle \ge 0$ implies $\lim_{k \to \infty} \langle p_k, \xi_{2,k} \rangle = 0$, while the bound (36) ensures that $\lim_{k \to \infty} p_k = 0$. It follows from Lemma 3 that $\lim_{k \to \infty} \text{Dist}(\xi_{2,k}, \partial \tilde{g}(0)) = 0$ but, according to (28b), we have $\xi_{1,k} + \xi_{2,k} \to 0$ as $k \to \infty$, so $\lim_{k \to \infty} \text{Dist}(-\xi_{1,k}, \partial \tilde{g}(0)) = 0$. By (31), there exist $\varepsilon > 0$ and $k^* \in \mathbb{N}$ such that, for all $k \ge k^*$, $\xi_{1,k+1} + b_{k+1} \in \partial P(0)$ for any $b_{k+1} \in \varepsilon \mathcal{B}$. From $\xi_{1,k+1} \in \partial P(q_{k+1})$ and the monotonicity of subdifferentials we know that $\langle b_{k+1}, q_{k+1} \rangle \ge 0$, for any $b_{k+1} \in \varepsilon \mathcal{B}$ and all $k \ge k^*$. Therefore, $q_{k+1} = 0$ for all $k \ge k^*$, hence $p_{k+2} = 0$ for all $k \ge k^*$ in view of (28). ■

## 4.2 Douglas-Rachford-like Splitting Discretisation

Lemma 2 and Theorem 3 show that, when the dissipation function depends only on the momenta $p$, the discrete-time system (28) is well-posed and finite-time stable. However, an explicit expression for the proximal map (29) is, in general, hard to obtain. Nevertheless, approximations can be computed via splitting schemes. Let us start by considering the following system:

$$ q_{k+1} = q_k + h p_{k+1} \quad (37a) $$

$$ p_{k+1} = p_k - h \xi_{1,k+1} - h \xi_{2,k+1} \quad (37b) $$

$$ \xi_{1,k+1} \in \partial P(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \quad (37c) $$

$$ \xi_{2,k+1} \in \partial \tilde{g}(p_{k+1}). \quad (37d) $$

This can be seen as a dynamic extension of (28). Indeed, the particular value of $\xi_{2,k+1} \in \partial \tilde{g}(p_{k+1})$ now depends implicitly on $\xi_{2,k}$, making $\xi_2$ a new state variable. Since $\partial P$ is osc, we know that, for any $\varepsilon > 0$, there exists $h > 0$ sufficiently small such that

$$ \partial P(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \subset \partial P(q_{k+1}) + \varepsilon \mathcal{B}. \quad (38) $$

If $\partial \tilde{g}$ is uniformly bounded, then for all $k \in \mathbb{N}$, $\|\xi_{2,k}\| \le M$ for some $M > 0$. Thus, $q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k}) \in \{q_{k+1}\} + 2h^2 M \cdot \mathcal{B} = \{q_{k+1}\} + \mathcal{O}(h^2) \cdot \mathcal{B}$. In this sense, the scheme (37) approximates (28). We now address its well-posedness.
Notice that both lim inf in (87) are well-defined, as their arguments are bounded. Moreover, it follows from (85) and (64a) that $\lim_{k \to \infty} q_k$ exists, as $\{q_k\}_{k \in \mathbb{N}}$ is a Cauchy sequence. Thus, (87) becomes

$$ 
\begin{aligned}
\lim_{k \to \infty} q_k + h^2 \liminf_{k \to \infty} \xi_{2,k+1} &= \lim_{k \to \infty} -hp_{k+1} + h^2 \liminf_{k \to \infty} \xi_{2,k} \\
&= h^2 \liminf_{k \to \infty} \xi_{2,k}
\end{aligned}
\quad (88) $$

and the convergence of $\{q_k\}_{k \in \mathbb{N}}$ towards zero follows. The proof is complete. □

Let us illustrate the above with the twisting algorithm. In addition it will be proved for this particular (yet important) case, that the finite-time stability is obtained. **Example 6 (Twisting Algorithm).** The twisting algorithm (21) has a dissipation function that depends only on the momenta $p$. Indeed, it satisfies Assumption 2 with $P(q) = \gamma_1|q|$ and $\tilde{g}(p) = \gamma_2|p|$. Thus, its discretised version

$$ q_{k+1} = q_k + hp_{k+1} \quad (89a) $$

$$ p_{k+1} = p_k - h\xi_{1,k+1} - h\xi_{2,k+1} \quad (89b) $$

$$ \xi_{1,k+1} \in \gamma_1 \operatorname{sgn}(q_{k+1}) \quad (89c) $$

$$ \xi_{2,k+1} \in \gamma_2 \operatorname{sgn}(p_{k+1}) \quad (89d) $$

is globally finite-time stable whenever $\gamma_1 > \gamma_2 > 0$, in view of Theorem 3. However, an explicit expression for the proximal map (46) (and for the selection values $\xi_{1,k}$ and $\xi_{2,k}$) is, in general, hard to obtain. Nevertheless, the approximation proposed in (64) can be computed. Indeed, simple computations show that

$$ \operatorname{Prox}_{h^2P}(s) = s - h^2\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{s}{h^2\gamma_1}\right), \quad (90a) $$

$$ \operatorname{Prox}_{h\tilde{g}}(s) = s - h\gamma_2 \operatorname{Proj}\left([-1, 1]; \frac{s}{h\gamma_2}\right), \quad (90b) $$

so that (64) becomes

$$ q_{k+1} = q_k + hp_{k+1} \quad (91a) $$

$$ p_{k+1} = p_k - h\xi_{1,k+1} - h\xi_{2,k+1} \quad (91b) $$

$$ \xi_{1,k+1} = \gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{q_k + hp_k - h^2\xi_{2,k}}{h^2\gamma_1}\right) \quad (91c) $$

$$ \xi_{2,k+1} = \gamma_2 \operatorname{Proj}\left([-1, 1]; \frac{p_k - h\xi_{1,k+1}}{h\gamma_2}\right) \quad (91d) $$

Like the origin of (89), the origin (91) is also finite-time stable. Indeed, from (64) simple computations show that $\Gamma = \operatorname{Fix} T$, so that $T^{-1}(\Gamma)$ is computed from the implicit

<figure><img src="image_6.png" /><figcaption>Fig. 2. Discrete orbits of backward-Euler implementation of the twisting algorithm (44) with $P(\cdot) = \gamma_1|\cdot|$, $\tilde{g} = \gamma_2|\cdot|$, $\gamma_1 = 4$, $\gamma_2 = 2$, and $h = 0.3s$. The red line denotes the set $T^{-1}(0)$, whereas the light blue region denotes the set $T^{-2}(0)$.</figcaption></figure>

model (64) as

$$ 
\begin{aligned}
T^{-1}(\Gamma) &= \left\{ \begin{bmatrix} q \\ p \\ \xi_2 \end{bmatrix} \in \mathbb{R}^3 \mid q = 0, |p - h\xi_1| \le h\gamma_2, \right. \\
& \qquad \left. |\xi_1| \le \gamma_1 \right\} \\
&= \{0\} \times [-h(\gamma_1 + \gamma_2), h(\gamma_1 + \gamma_2)] \times \mathbb{R}, \quad (92)
\end{aligned}
 $$

whereas the second pre-image of $\Gamma$ is given as

$$ 
\begin{aligned}
T^{-2}(\Gamma) &= \left\{ \begin{bmatrix} q \\ p \\ \xi_2 \end{bmatrix} \in \mathbb{R}^3 \mid 0 = q + h\bar{p}, \right. \\
& \qquad \left. \bar{p} \in p - h\xi_1 - h\gamma_2 \operatorname{sgn}(\bar{p}), \quad |\bar{p}| \le h(\gamma_1 + \gamma_2), \right. \\
& \qquad \left. \xi_1 = \gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{q + hp - h^2\xi_2}{h^2\gamma_1}\right) \right\} \\
&= \left\{ \begin{bmatrix} q \\ p \\ \xi_2 \end{bmatrix} \in \mathbb{R}^3 \mid q = -h \operatorname{Prox}_{h\tilde{g}}(p - h\xi_1), \right. \\
& \qquad \left. |q| \le h^2(\gamma_1 + \gamma_2), \right. \\
& \qquad \left. \xi_1 = \gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{q + hp - h^2\xi_2}{h^2\gamma_1}\right) \right\} \quad (93)
\end{aligned}
 $$

Hence, $\Gamma$ is finite-time stable in view of Proposition 2. Notice that the above algorithm is different from the discrete-time twisting studied in [19], which is studied with a ZOH discretisation of the plant, hence yielding a closed-loop different from (89). △

**Example 7 (Multivariable twisting algorithm).** Con-
**Lemma 4** Consider a difference inclusion (37) satisfying Assumption 2. For each current state $(q_k, p_k, \xi_{2,k})$ there exists a unique future state

$$ (q_{k+1}, p_{k+1}, \xi_{2,k+1}) = T(q_k, p_k, \xi_{2,k}) \quad (39) $$

Furthermore, the unique selections $\xi_{1,k+1}$, $\xi_{2,k+1}$ are determined as

$$ \xi_{1,k+1} = \frac{1}{h^2} (I_d - \text{Prox}_{h^2 P}) (q_k + hp_k - h^2\xi_{2,k}) \quad (40a) $$

$$ \xi_{2,k+1} = \frac{1}{h} (I_d - \text{Prox}_{h\tilde{g}}) (p_k - h\xi_{1,k+1}). \quad (40b) $$

**PROOF.** Define $\theta_{k+1} = q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})$. It follows from (37a)-(37c) that

$$ \begin{aligned} \theta_{k+1} + h^2 \xi_{1,k+1} &= q_k + hp_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k}) \\ & \quad + h^2 \xi_{1,k+1} \\ &= q_k + hp_k - h^2 \xi_{2,k} \end{aligned} \quad (41) $$

Since $\xi_{1,k+1} \in \partial P(\theta_{k+1})$, we have

$$ \theta_{k+1} = \text{Prox}_{h^2 P} (q_k + hp_k - h^2 \xi_{2,k}). \quad (42) $$

Substituting (42) into (41) and solving for $\xi_{1,k+1}$ gives (40a), while (37b) gives

$$ p_{k+1} + h\xi_{2,k+1} = p_k - h\xi_{1,k+1}. \quad (43) $$

Hence, (37d) implies that

$$ p_{k+1} = \text{Prox}_{h\tilde{g}}(p_k - h\xi_{1,k+1}). \quad (44) $$

Substituting (44) into (43) and solving for $\xi_{2,k+1}$ gives (40b). Thus, $p_{k+1}$, and hence $q_{k+1}$, depend on data available at time $t_k$ so that (37) is well-defined. This concludes the proof. ■

Lemma 4 shows the explicit expression for the selection values $\xi_{1,k+1}$ and $\xi_{2,k+1}$. Notice that $\xi_{2,k+1}$ depends on its past values, so that the full state of the proposed approximation (37) belongs to $\mathbb{R}^{3n}$.

**Remark 4** As introduced in [22] for the case of multi-variable super-twisting-like schemes, the approximation of the selection values in (40) is reminiscent of the celebrated Douglas-Rachford splitting [21,14]. Indeed, set $\vartheta_{k+1} = p_{k+1} + h\xi_{2,k+1}$ and note that it follows from (37d) that $p_{k+1} = \text{Prox}_{h\tilde{g}}(\vartheta_{k+1})$, so that after simple computations we arrive at the expressions

$$ \vartheta_{k+1} = \frac{1}{h} \left( \text{Prox}_{h^2 P} \left( q_k + h(2 \text{Prox}_{h\tilde{g}} - I_d)(\vartheta_k) \right) - q_k \right) + (I_d - \text{Prox}_{h\tilde{g}})(\vartheta_k) $$

$$ p_{k+1} = \text{Prox}_{h\tilde{g}}(\vartheta_{k+1}), $$

which define the Douglas-Rachford splitting associated with the generalised equation $0 \in h\partial P(q + hp) + h\partial\tilde{g}(p)$, see, e.g., [30, Section 2.7].

In what follows, the stability properties of the set

$$ \Gamma = \{(q, p, \xi_2) \in \mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R}^n \mid q = p = 0, \xi_2 \in \mathbb{R}^n\} \quad (45) $$

are studied. Assumption 3 is strengthened.

**Theorem 4** Consider (37), the backward-Euler discretisation of a system satisfying

$$ \mu\mathcal{B} \subset -\partial\tilde{g}(0) \subset \int \partial P(0) \quad (46) $$

for some $\mu > 0$. Assume $\partial\tilde{g}$ is uniformly bounded. Then, there exists $h^* > 0$ such that, for all $0 < h < h^*$, the set $\Gamma$ in (45) is globally asymptotically stable.

**PROOF.** Let us consider the function $H$ introduced in Assumption 2. Using again the definition of convex subdifferential, we can show that

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \leq \langle \eta_{1,k+1}, q_{k+1} - q_k \rangle + \langle \eta_{2,k+1}, p_{k+1} - p_k \rangle \quad (47) $$

for any $\eta_{1,k+1} \in \partial P(q_{k+1})$ and $\eta_{2,k+1} = p_{k+1}$. It follows from (37c), (38), and the fact that $\partial P$ is osc, that, for any $0 < \varepsilon_1 < \mu$, there exists $h > 0$ sufficiently small such that

$$ \xi_{1,k+1} \in \partial P(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \subset \partial P(q_{k+1}) + \varepsilon_1 \mathcal{B}. $$

It is inferred that, for any $k \in \mathbb{N}$, there always exists $b_{k+1} \in \mathcal{B}$ such that $\xi_{1,k+1} + \varepsilon_1 b_{k+1} \in \partial P(q_{k+1})$. Hence, taking $\eta_{1,k+1} = \xi_{1,k+1} + \varepsilon_1 b_{k+1}$ and $\eta_{2,k+1} = p_{k+1}$ in (47) leads us to

$$ \begin{aligned} H(q_{k+1}, p_{k+1}) - H(q_k, p_k) &\leq \langle \xi_{1,k+1} + \varepsilon_1 b_{k+1}, q_{k+1} - q_k \rangle \\ & \quad + \langle p_{k+1}, p_{k+1} - p_k \rangle \\ &= h \langle \xi_{1,k+1} + \varepsilon_1 b_{k+1}, p_{k+1} \rangle - h \langle p_{k+1}, \xi_{1,k+1} + \xi_{2,k+1} \rangle \\ &\leq h \varepsilon_1 \|p_{k+1}\| - h \langle p_{k+1}, \xi_{2,k+1} \rangle. \end{aligned} \quad (48) $$

It follows from (46) that

$$ \langle p_{k+1}, \xi_{2,k+1} \rangle > \mu \|p_{k+1}\|. \quad (49) $$

To see this, note that $-\mu b_k \in \partial\tilde{g}(0)$ for any $b_k \in \mathcal{B}$ and that, by the maximal monotonicity of $\partial\tilde{g}$, we have $\langle p_k, \xi_{2,k} - \mu b_k \rangle \geq 0$. In other words, $\langle p_{k+1}, \xi_{2,k+1} \rangle \geq \mu \langle p_{k+1}, b_{k+1} \rangle$ for all $b_{k+1} \in \mathcal{B}$. Therefore, $\langle p_{k+1}, \xi_{2,k+1} \rangle \geq \mu \sup_{b_{k+1} \in \mathcal{B}} \langle p_{k+1}, b_{k+1} \rangle$, which implies (49). Hence, (48) is further bounded as

$$ H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \leq -h(\mu - \varepsilon_1) \|p_{k+1}\|. $$
<figure><img src="image_7.png" /><figcaption>Fig. 3. Discrete orbits of the twisting algorithm with splitting (64) with $P(\cdot) = \gamma_1|\cdot|$, $\tilde{g} = \gamma_2|\cdot|$, $\gamma_1 = 4$, $\gamma_2 = 2$, and $h = 0.3s$. The red line denotes the projection onto the plane $q-p$ of the set $T^{-1}(\Gamma)$, whereas the light blue region denotes the projection of the set $T^{-2}(\Gamma)$.</figcaption></figure>

sider the multivariable twisting scheme

$$q_{k+1} = q_k + hp_{k+1} \quad (94a)$$

$$p_{k+1} = p_k - h\xi_{1,k+1} - h\xi_{2,k+1} \quad (94b)$$

$$\xi_{1,k+1} \in \partial P(q_{k+1}) \quad (94c)$$

$$\xi_{2,k+1} \in \partial \tilde{g}(p_{k+1}), \quad (94d)$$

where $q_k \in \mathbb{R}^n$, $p_k \in \mathbb{R}^n$ for all $k \in \mathbb{N}$, $P(\cdot) = \gamma_1\|\cdot\|_{r_1}$, and $\tilde{g}(\cdot) = \gamma_2\|\cdot\|_{r_2}$. It is easy to verify that, for the case of $n=1$, (94) reduces to the classical twisting algorithm. The finite-time stability of (94) is guaranteed by Theorem 3 whenever Assumption 2 and 3 hold. Clearly, Assumption 2 holds for the choice of maps $P$ and $\tilde{g}$ considered above, whereas for the specific case of (94), (49) becomes

$$0 \in \gamma_2 \operatorname{cl} \mathcal{B}_{s_2} \subset \gamma_1 \mathcal{B}_{s_1}, \quad (95)$$

where $1/r_i + 1/s_i = 1$, for $i \in \{1, 2\}$. Indeed, let $f(x) = \|x\|_r$ for $r \in [1, +\infty]$. Hence,

$$\begin{aligned} \partial f(0) &= \{\xi \in \mathbb{R}^n \mid \langle \xi, \eta \rangle \le \|\eta\|_r \text{ for all } \eta \in \mathbb{R}^n\} \\ &= \left\{ \xi \in \mathbb{R}^n \mid \sup_{\eta \in \mathbb{R}^n} \{\langle \xi, \eta \rangle - \|\eta\|_r\} \le 0 \right\} \\ &= \{\xi \in \mathbb{R}^n \mid f^*(\xi) \le 0\} \\ &= \operatorname{cl} \mathcal{B}_s \end{aligned} \quad (96)$$

The last equality follows from the fact that, for $f(\cdot) = \|\cdot\|_r$, $f^*(\cdot) = \Psi_{\operatorname{cl}\mathcal{B}_s}(\cdot)$, where $1/r + 1/s = 1$ and $\Psi_C$ denotes the classical indicator function of the set C. If, for instance, $r_1 = r_2$, then $s_1 = s_2$ and (95) holds whenever $\gamma_1 > \gamma_2$, which is the same condition as (20) for the finite-time stability of the conventional twisting algorithm. Thus, if (95) holds, then the origin is finite-time stable in view of Theorem 3. Figure 4 depicts the time evolution of the state trajectories of (94) for the case of $h = 0.1$, $n = 3$, $r_1 = 1$, $r_2 = 2$, $\gamma_1 = 5$ and $\gamma_2 = 3$. It is shown that all trajectories converge to zero in finite time.

<figure><img src="image_8.png" /><figcaption>Fig. 4. Time evolution of state trajectories of (44) for several initial conditions. In this case $h = 0.1$, $n = 3$, $P(\cdot) = \gamma_1\|\cdot\|_1$, $\tilde{g}(\cdot) = \|\cdot\|_2$, $\gamma_1 = 5$ and $\gamma_2 = 3$. For this simulation the state $p_{k+1}$ is computed via (46) using the solver CVXOPT at each iteration.</figcaption></figure>

Now, let us consider the splitting scheme (64) applied to the multivariable twisting algorithm. In this case Lemma 4 provides the explicit expressions for the control selections $\xi_{1,k+1}$ and $\xi_{2,k+1}$, which for the case of $P(\cdot) = \gamma_1\|\cdot\|_{r_1}$ and $\tilde{g}(\cdot) = \gamma_2\|\cdot\|_{r_2}$, they become

$$\xi_{1,k+1} = \gamma_1 \operatorname{Proj} \left( \mathcal{B}_{s_1}; \frac{1}{h^2\gamma_1} (q_k + hp_k - h^2\xi_{2,k}) \right) \quad (97a)$$

$$\xi_{2,k+1} = \gamma_2 \operatorname{Proj} \left( \mathcal{B}_{s_2}; \frac{1}{h\gamma_2} (p_k - h\xi_{1,k+1}) \right) \quad (97b)$$

where we have used Moreau's decomposition, see e.g., [6, Theorem 14.3] and the fact that $f^*(\cdot) = \Psi_{\mathcal{B}_s}(\cdot)$ for $f = \|\cdot\|_r$ with $1/r + 1/s = 1$. Figure 5 shows the time evolution of the state trajectories with the same set of parameters as for the previous case. It is shown that all trajectories converge to zero in finite time. △

## 4.3 Position-dependent dissipation

As mentioned before, the model (14) with $\partial_q g = \{0\}$ encompasses many mechanical systems with set-valued frictional forces. Inspired by the super-twisting example, we propose a different, somewhat complementary set of assumptions.

**Assumption 4.** The functions $H$ and $g$ in (14) satisfy the following:

(1) The Hamiltonian function is of the form $H(q, p) = P(q) + \frac{1}{2}\|p\|^2$ for some convex function $P: \mathbb{R}^n \to \mathbb{R}$.

(2) The subdifferential of the dissipation potential $g$ is single valued and $g$ depends on $q$ only, i.e., $g(q, p) = \hat{g}(q)$ for some convex function $\hat{g}$.
Since by assumption $0 < \varepsilon_1 < \mu$, it follows that $h(\mu - \varepsilon_1)\|p_{k+1}\| \ge 0$ for all $k \in \mathbb{N}$. Stability of the set $\Gamma$ is concluded. To show asymptotic convergence, we apply same arguments as in the proof of Theorem 3 (cf. last paragraph), leading to

$$ \lim_{k \to \infty} p_k = 0 \quad \text{and} \quad \lim_{k \to \infty} \text{Dist}(-\xi_{1,k}, \partial \tilde{g}(0)) = 0. \quad (50) $$

Continuing with the same steps as in the proof of Theorem 3, we can conclude from (37c), (46), and (50) that there exists $k^* \in \mathbb{N}$ such that, for all $k > k^*$,

$$ q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k}) = 0. \quad (51) $$

Now, the substitution of (37a) into (51) yields

$$ \liminf_{k \to \infty} (q_k + h^2 \xi_{2,k+1}) = \liminf_{k \to \infty} (-hp_{k+1} + h^2 \xi_{2,k}), \quad (52) $$

where the $\liminf$ is taken component-wise. Notice that both $\liminf$ in (52) are well-defined, as their arguments are bounded. Moreover, it follows from (50) and (37a) that $\lim_{k \to \infty} q_k$ exists, as $\{q_k\}_{k \in \mathbb{N}}$ is a Cauchy sequence. Thus, (52) becomes

$$ \lim_{k \to \infty} q_k + h^2 \liminf_{k \to \infty} \xi_{2,k+1} = \lim_{k \to \infty} -hp_{k+1} + h^2 \liminf_{k \to \infty} \xi_{2,k} \\ = h^2 \liminf_{k \to \infty} \xi_{2,k} $$

and the convergence of $\{q_k\}_{k \in \mathbb{N}}$ towards zero follows. ■

Let us illustrate the above with the twisting algorithm. In addition it will be proved for this particular (yet important) case, that the finite-time stability is obtained.

### Example 6 (Multivariable twisting algorithm)

Consider the multivariable twisting scheme

$$ q_{k+1} = q_k + hp_{k+1} \quad (53a) $$

$$ p_{k+1} = p_k - h\xi_{1,k+1} - h\xi_{2,k+1} \quad (53b) $$

$$ \xi_{1,k+1} \in \partial P(q_{k+1}) \quad (53c) $$

$$ \xi_{2,k+1} \in \partial \tilde{g}(p_{k+1}), \quad (53d) $$

where $q_k \in \mathbb{R}^n$, $p_k \in \mathbb{R}^n$ for all $k \in \mathbb{N}$, $P(\cdot) = \gamma_1 \| \cdot \|_{r_1}$, and $\tilde{g}(\cdot) = \gamma_2 \| \cdot \|_{r_2}$. It is easy to verify that, for the case of $n=1$, (53) reduces to the classical twisting algorithm. The finite-time stability of (53) is guaranteed by Theorem 3 whenever Assumption 2 and 3 hold. Clearly, Assumption 2 holds for the choice of maps $P$ and $\tilde{g}$ considered above, whereas for the specific case of (53), (31) becomes

$$ 0 \in \gamma_2 \mathrm{cl} \mathcal{B}_{s_2} \subset \gamma_1 \mathcal{B}_{s_1}, \quad (54) $$

where $1/r_i + 1/s_i = 1$, for $i \in \{1, 2\}$. Indeed, let $f(x) =$

$\|x\|_r$ for $r \in [1, +\infty]$. Hence,

$$ \begin{aligned} \partial f(0) &= \{\xi \in \mathbb{R}^n \mid \langle \xi, \eta \rangle \le \|\eta\|_r \text{ for all } \eta \in \mathbb{R}^n\} \\ &= \{\xi \in \mathbb{R}^n \mid \sup_{\eta \in \mathbb{R}^n} \{\langle \xi, \eta \rangle - \|\eta\|_r\} \le 0\} \\ &= \{\xi \in \mathbb{R}^n \mid f^*(\xi) \le 0\} = \mathrm{cl} \mathcal{B}_s \end{aligned} $$

The last equality follows from the fact that, for $f(\cdot) = \| \cdot \|_r$, $f^*(\cdot) = \Psi_{\mathrm{cl} \mathcal{B}_s}(\cdot)$, where $1/r + 1/s = 1$ and $\Psi_C$ denotes the classical indicator function of the set $C$. If, for instance, $r_1 = r_2$, then $s_1 = s_2$ and (54) holds whenever $\gamma_1 > \gamma_2$, which is the same condition as $\gamma_1 > \gamma_2 > 0$ for the finite-time stability of the conventional twisting algorithm. Thus, if (54) holds, then the origin is finite-time stable in view of Theorem 3. Figure 2 (top) depicts the time evolution of the Hamiltonian function of (53) for $n=3$, $r_1=1$, $r_2=2$, $\gamma_1=7$, $\gamma_2=5$, and several values of the step-size. It can be verified that in all cases the full state converges to zero in finite time.

Now, let us consider the splitting scheme (37) applied to the multivariable twisting algorithm. In this case Lemma 4 provides the explicit expressions for the control selections $\xi_{1,k+1}$ and $\xi_{2,k+1}$. For the case of $P(\cdot) = \gamma_1 \| \cdot \|_{r_1}$ and $\tilde{g}(\cdot) = \gamma_2 \| \cdot \|_{r_2}$, they become

$$ \xi_{1,k+1} = \gamma_1 \mathrm{Proj} \left( \frac{1}{h^2\gamma_1} (q_k + hp_k - h^2\xi_{2,k}), \mathcal{B}_{s_1} \right) \quad (55a) $$

$$ \xi_{2,k+1} = \gamma_2 \mathrm{Proj} \left( \frac{1}{h\gamma_2} (p_k - h\xi_{1,k+1}), \mathcal{B}_{s_2} \right) \quad (55b) $$

where Moreau's decomposition was used [5, Theorem 14.3] and the fact that $f^*(\cdot) = \Psi_{\mathcal{B}_s}(\cdot)$ for $f = \| \cdot \|_r$ with $1/r + 1/s = 1$. Figure 2 (middle) shows the time evolution of the state trajectories with the same set of parameters. It is inferred that the full state converges to zero in finite time for values of $h$ sufficiently small. Finally, we also consider the case when the forward-Euler discretisation is used for the computation of the control $u_k = -\xi_{1,k+1} - \xi_{2,k+1}$, that is

$$ \xi_{1,k+1}^{\top} = \gamma_1 \left( \frac{q_{1,k}}{|q_{1,k}|} \frac{q_{2,k}}{|q_{2,k}|} \frac{q_{3,k}}{|q_{3,k}|} \right)^{\top} \quad (56a) $$

$$ \xi_{2,k+1} = \frac{\gamma_2}{\|p_k\|_2} p_k. \quad (56b) $$

Figure 2 (bottom) shows the associated Hamiltonians for different values of the step-size: the full state converges to zero only when $h$ is zero. Chattering is visible.

Δ

### 4.3 Position-dependent dissipation

As mentioned before, Hamiltonian systems with potential dissipation functions depending on $p$ alone encom-
<figure><img src="image_9.png" /><figcaption>Fig. 5. Time evolution of state trajectories of (64) for several initial conditions. In this case $h = 0.1$, $n = 3$, $P(\cdot) = \gamma_1 \|\cdot\|_1$, $\tilde{g}(\cdot) = \gamma_2 \|\cdot\|_2$, $\gamma_1 = 5$ and $\gamma_2 = 3$. The explicit expressions for the values of selections $\xi_{1,k+1}$ and $\xi_{2,k+1}$ are given in (97).</figcaption></figure>

**Lemma 5.** Consider the differential inclusion (14) and let $H$ and $g$ satisfy Assumption 4. For each current state $(q_k, p_k)$, there exists a unique future state $(q_{k+1}, p_{k+1})$ satisfying the backward-Euler discretisation scheme

$$q_{k+1} \in q_k + h \partial_p H(q_{k+1}, p_{k+1}) - h \partial \hat{g}(q_{k+1}) \quad (98a)$$

$$p_{k+1} \in p_k - h \partial_q H(q_{k+1}, p_{k+1}). \quad (98b)$$

*Proof.* By taking into account the assumptions and substituting (98b) into (98a), we can rewrite (98) as

$$q_{k+1} \in q_k + hp_k - h\partial\hat{g}(q_{k+1}) - h^2\xi_{2,k+1} \quad (99a)$$

$$p_{k+1} = p_k - h\xi_{2,k+1} \quad (99b)$$

with $\xi_{2,k} \in \partial P(q_k)$ the potential force. Equation (99a) is equivalent to the inclusion

$$q_{k+1} \in q_k + hp_k - h\partial G(q_{k+1}), \quad (100)$$

where $G$ is a new function defined by $G(q) = \hat{g}(q) + hP(q)$. By the maximal monotony of $G$, we can uniquely compute $q_{k+1}$ as

$$q_{k+1} = \text{Prox}_{hG}(q_k + hp_k). \quad (101)$$

Having $q_{k+1}$ fixed, we can uniquely determine the potential force $\xi_{2,k+1} \in \partial P(q_{k+1})$ from (99a)

$$\xi_{2,k+1} = \frac{p_k - \partial \hat{g}(q_{k+1})}{h} - \frac{q_{k+1} - q_k}{h^2}. \quad (102)$$

Finally, we can uniquely compute $p_{k+1}$ from (99b). □

**Theorem 5.** Consider the backward-Euler discretisation (98) of an inclusion (14) satisfying Assumptions 1 and 4. Suppose, further, that

$$\text{Zero } \partial P = 0 \quad (103)$$

and

$$\langle \eta, \partial \hat{g}(q) \rangle \ge \|q\|^m \quad (104)$$

for some $m > 0$ and all $\eta \in \partial P(q)$. Then, the origin is globally asymptotically stable.

Condition (103) implies that $H$ has a unique minimum and that it takes place at $(0,0)$. The latter is without loss of generality up to a simple change of coordinates. The condition (104) ensures that the dissipation is zero if and only if $q = 0$.

*Proof.* By the definition of the subdifferential, we have

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, q_{k+1} - q_k \rangle + \langle \eta_{2,k+1}, p_{k+1} - p_k \rangle \quad (105)$$

for all $\eta_{1,k} \in \partial P(q_k)$ and $\eta_{2,k} = p_k$. We have $q_{k+1} - q_k = hp_{k+1} - h\partial \hat{g}(q_{k+1})$ and it follows from Lemma 5 that, for every $k \in \mathbb{N}$, there exists a unique vector $\xi_{1,k} \in \partial P(q_k)$ such that $p_{k+1} - p_k = h\xi_{1,k+1}$. Thus,

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, hp_{k+1} - h\partial \hat{g}(q_{k+1}) \rangle + \langle p_{k+1}, h\xi_{1,k+1} \rangle. \quad (106)$$

This is true, in particular, for $\eta_{1,k} = \xi_{1,k}$, which implies that

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le -h \langle \xi_{1,k+1}, \partial \hat{g}(q_{k+1}) \rangle. \quad (107)$$

By (104), we have

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le 0, \quad (108)$$

which shows that $H$ is a Lyapunov function for the set of its minima.

By summing both sides of (107) we obtain

$$\sum_{k=0}^{N} (H(q_{k+1}, p_{k+1}) - H(q_k, p_k)) = H(q_{N+1}, p_{N+1}) - H(q_0, p_0) \le -h \sum_{k=0}^{N} \langle \xi_{1,k+1}, \partial \hat{g}(q_{k+1}) \rangle, \quad (109)$$

Since the sequence $\{H(q_N, p_N)\}_{N \in \mathbb{N}}$ is decreasing and bounded from below, we conclude that it converges to a
<figure><img src="image_3.png" /><figcaption>Fig. 2. Hamiltonians for three multivariable twisting algorithms. (Top) $H_a$ for closed-loop (53). (Middle) $H_b$ for closed-loop (53a)-(53b) with control (55). (Bottom) $H_c$ for closed-loop (53a)-(53b) with control (56).</figcaption></figure>

pass many mechanical systems with set-valued friction. Inspired by the super-twisting example, we propose a different, somewhat complementary set of assumptions.

**Assumption 4** The system satisfies Assumption 1 with $g$ depending only on $q$, i.e., $g(q, p) = \hat{g}(q)$ for some convex $C^1$ function $\hat{g}$.

**Lemma 5** Consider continuous-time system satisfying Assumption 4. For each current state $(q_k, p_k)$, there exists a unique future state $(q_{k+1}, p_{k+1})$ satisfying the backward-Euler discretisation scheme

$$q_{k+1} \in q_k + h \partial_p H(q_{k+1}, p_{k+1}) - h \partial \hat{g}(q_{k+1}) \quad (57a)$$

$$p_{k+1} \in p_k - h \partial_q H(q_{k+1}, p_{k+1}). \quad (57b)$$

**PROOF.** By taking into account the assumptions and substituting (57b) into (57a), we can rewrite (57) as

$$q_{k+1} \in q_k + hp_k - h\partial\hat{g}(q_{k+1}) - h^2\xi_{2,k+1} \quad (58a)$$

$$p_{k+1} = p_k - h\xi_{2,k+1} \quad (58b)$$

with $\xi_{2,k+1} \in \partial P(q_{k+1})$ the potential force. Equation (58a) is equivalent to the inclusion $q_{k+1} \in q_k + hp_k - h\partial G(q_{k+1})$, where $G$ is a new function defined by $G(q) = \hat{g}(q) + hP(q)$. By the maximal monotonicity of $G$, we can uniquely compute $q_{k+1}$ as

$$q_{k+1} = \text{Prox}_{hG}(q_k + hp_k) \quad (59a)$$

Having $q_{k+1}$ fixed, we can uniquely determine the potential force $\xi_{2,k+1} \in \partial P(q_{k+1})$ from (58a),

$$\xi_{2,k+1} = \frac{p_k - \partial \hat{g}(q_{k+1})}{h} - \frac{q_{k+1} - q_k}{h^2}.$$

Finally, we can uniquely compute $p_{k+1}$ from (58b) as

$$p_{k+1} = \frac{1}{h} (I_d + h \nabla \hat{g})(\text{Prox}_{hG}(q_k + hp_k)) - q_k). \quad (59b)$$

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

**Remark 5** Assumption 4 is sufficient to guarantee the well-posedness of the closed-loop (57). It can be relaxed as follows: Let $D_1$ and $D_2$ the (zero Lebesgue-measure) sets where the functions $\hat{g}$ and $P$ are non-differentiable. Then, if $D_1 \cap D_2 = \emptyset$ then the closed-loop is well-posed, see [22, Theorem 3].

**Remark 6** Similar to the case with momentum-dependent dissipation, splitting methods can be used for the explicit computation of the selections $\xi_{1,k+1} \in \partial \hat{g}(q_{k+1})$ and $\xi_{2,k+1} \in \partial P(q_{k+1})$ in (58), by considering the system

$$q_{k+1} = q_k + hp_k - h\xi_{1,k+1} - h^2\xi_{2,k+1} \quad (60a)$$

$$p_{k+1} = p_k - h\xi_{2,k+1} \quad (60b)$$

$$\xi_{1,k+1} \in \partial \hat{g}(q_{k+1} + h^2(\xi_{2,k+1} - \xi_{2,k})) \quad (60c)$$

$$\xi_{2,k+1} \in \partial P(q_{k+1}) \quad (60d)$$

The computations are similar to those shown in the proof of Lemma 4 and are omitted here because of space limitations.

**Theorem 5** Consider the backward-Euler discretisation (57) of a continuous-time system satisfying Assumption 4. Suppose, further, that

$$0 \in \text{int} \partial P(0) \quad (61)$$

and

$$\langle \eta, \partial \hat{g}(q) \rangle \ge \|q\|^a \quad (62)$$

for some $a > 0$ and all $\eta \in \partial P(q)$. Then, the origin is globally finite-time stable.

Condition (61) implies that $H$ has a unique minimum at $(0, 0)$. The latter is without loss of generality up to a simple change of coordinates. The condition (62) ensures that the dissipation is zero if and only if $q = 0$.

**PROOF.** From the definition of the subdifferential,

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \le \langle \eta_{1,k+1}, q_{k+1} - q_k \rangle + \langle \eta_{2,k+1}, p_{k+1} - p_k \rangle$$
limit $\bar{H} \ge \min_{q,p} H(q,p)$. Thus,

$$h \sum_{k=0}^{\infty} \langle \xi_{1,k+1}, \partial \hat{g}(q_{k+1}) \rangle = H(q_0, p_0) - \bar{H}, \quad (110)$$

which implies that $\lim_{k \to \infty} \langle \xi_{1,k}, \hat{g}(q_k) \rangle = 0$. It follows from (104) that

$$\lim_{k \to \infty} q_k = 0. \quad (111)$$

Finally, we can see from (98a) that (111) implies that

$$\lim_{k \to \infty} p_k = 0. \quad (112)$$

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

**Remark 2.** The inequality in (109) is a dissipation inequality showing that the system with input $\xi_{2k} = \partial \hat{g}(q_k)$, output $-\xi_{1,k}$, is passive with respect to the supply rate $\langle -h\xi_{1,k}, \xi_{2,k} \rangle$ and storage function $H(q_{k+1}, p_{k+1})$ [9, section 3.15].

**Example 8** (Super-Twisting Algorithm). The super-twisting algorithm (28) satisfies Assumption 4 with

$$P(q) = \gamma_2 |q| \quad \text{and} \quad \hat{g}(q) = \gamma_1 \frac{2}{3} |q|^{3/2}. \quad (113)$$

Indeed, the potential function $P$ has a unique minimum at $q=0$, and we have $\partial \hat{g}(0) = 0$ and

$$\partial \hat{g}(q) \cdot \partial P(q) = \gamma_1 \gamma_2 |q|^{1/2}, \quad (114)$$

that is, it satisfies (104). The latter also implies Assumption 1. Thus, the iteration satisfying

$$q_{k+1} = q_k + hp_{k+1} - h\gamma_1 [q_{k+1}]^{1/2} \quad (115a)$$

$$p_{k+1} \in p_k - h\gamma_2 \operatorname{sgn}(q_{k+1}) \quad (115b)$$

represents the origin globally asymptotically stable. The explicit expression associated with (115a) was obtained in [10] as

$$q_{k+1} = q_k + hp_k - h\gamma_1 \xi_{1,k+1} - h^2 \gamma_2 \xi_{2,k+1} \quad (116a)$$

$$p_{k+1} = p_k - h\gamma_2 \xi_{2,k+1}, \quad (116b)$$

where

$$\xi_{1,k+1} = \tilde{\beta}(q_k + hp_k) \operatorname{Proj} \left( [-1, 1]; \frac{q_k + hp_k}{h^2 \gamma_2} \right) = [q_{k+1}]^{\frac{1}{2}} \quad (117a)$$

$$\xi_{2,k+1} = \operatorname{Proj} \left( [-1, 1]; \frac{q_k + hp_k}{h^2 \gamma_2} \right) \in \operatorname{sgn}(q_{k+1}) \quad (117b)$$

and $\beta: \mathbb{R} \to \mathbb{R}_+$ is

$$\tilde{\beta}(s) = -\frac{h\gamma_1}{2} + \frac{1}{2}\sqrt{h^2\gamma_1^2 + 4 \max\{0, |s| - h^2\gamma_2\}}. \quad (117c)$$

<figure><img src="image_11.png" /><figcaption>Fig. 6. Discrete orbits of the backward-Euler implementation of the super-twisting algorithm (115) with $\gamma_1 = 4$, $\gamma_2 = 2$, and $h = 0.3$. The pre-image $T^{-1}(0)$ (red) does not contain a neighborhood of the origin, however, $T^{-2}(0)$ (light blue) does. Thus, discrete-time orbits (black) reach the origin in finite time.</figcaption></figure>

To establish finite-time convergence, we will begin by computing $T^{-1}(0, 0)$. This is easily accomplished directly in the implicit form of $T$. Setting $q_{k+1} = 0$ and $p_{k+1} = 0$ in (115) gives

$$T^{-1}(0, 0) = 0 \times [-h\gamma_2, h\gamma_2]. \quad (118)$$

The set does not contain a neighborhood of the origin. We proceed to compute $T^{-1}(0 \times [-h\gamma_2, h\gamma_2])$. Setting $q_{k+1} = 0$ in (115) yields the conditions

$$q_k = -hp_{k+1} \quad \text{and} \quad p_k \in p_{k+1} + [-h\gamma_2, h\gamma_2] \quad (119)$$

with $p_{k+1} \in [-h\gamma_2, h\gamma_2]$. Thus, the second pre-image of the origin is

$$T^{-2}(0, 0) = \{(q, p) \in \mathbb{R}^2 \mid |q + hp| \le h^2\gamma_2, |q| \le h^2\gamma_2\}. \quad (120)$$

## 5 Conclusions

Gradient and Hamiltonian systems play important and somewhat complementary roles in Dynamical Systems Theory [7]. Hamiltonian systems with Rayleigh dissipation incorporate both classes of systems, as the dissipative vector field effectively behaves as a gradient vector field. First-order sliding-mode control systems can be framed into the generalised class of multivalued gradient systems, also called subgradient systems [12]. Unfortunately, higher-order sliding-mode control systems cannot. However, we have showed that at least two second-order sliding-mode control systems can be recast as multivalued Hamiltonian systems with dissipation, providing further insight into higher-order sliding motions.

It is well-known that the properties of continuous-time subgradient systems are best preserved under backward-Euler discretisation schemes [12]. We have shown here
for all $\eta_{1,k+1} \in \partial P(q_{k+1})$ and $\eta_{2,k+1} = p_{k+1}$. We have $q_{k+1} - q_k = hp_{k+1} - h\partial\hat{g}(q_{k+1})$ and it follows from Lemma 5 that, for every $k \in \mathbb{N}$, there exists a unique vector $\xi_{1,k} \in \partial P(q_k)$ such that $p_{k+1} - p_k = h\xi_{1,k+1}$. Thus,

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \leq \\ \langle \eta_{1,k+1}, hp_{k+1} - h\partial\hat{g}(q_{k+1}) \rangle + \langle p_{k+1}, h\xi_{1,k+1} \rangle.$$

This is true, in particular, for $\eta_{1,k+1} = \xi_{1,k+1}$, which implies

$$H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \leq -h\langle\xi_{1,k+1}, \partial\hat{g}(q_{k+1})\rangle. \qquad (63)$$

By (62), we have $H(q_{k+1}, p_{k+1}) - H(q_k, p_k) \leq 0$, which shows that $H$ is a Lyapunov function for the set of its minima. By summing both sides of (63) we obtain

$$\sum_{k=0}^{N} (H(q_{k+1}, p_{k+1}) - H(q_k, p_k)) = \\ H(q_{N+1}, p_{N+1}) - H(q_0, p_0) \leq -h \sum_{k=0}^{N} \langle \xi_{1,k+1}, \partial \hat{g}(q_{k+1}) \rangle, \qquad (64)$$

Since the sequence $\{H(q_N, p_N)\}_{N \in \mathbb{N}}$ is decreasing and bounded from below, we conclude that it converges to a limit $\bar{H} \geq \min_{q,p} H(q,p)$. Thus,

$$h \sum_{k=0}^{\infty} \langle \xi_{1,k+1}, \partial \hat{g}(q_{k+1}) \rangle = H(q_0, p_0) - \bar{H},$$

which implies that $\lim_{k \to \infty} \langle \xi_{1,k}, \hat{g}(q_k) \rangle = 0$. It follows from (62) that

$$\lim_{k \to \infty} q_k = 0. \qquad (65)$$

From (57a), we see that (65) implies $\lim_{k \to \infty} p_k = 0$. We now show that the convergence is indeed in finite time. To that end, we compute the first two pre-images of the origin under the map $T$ defining the dynamics in (59) as follows:

$$T^{-1}(0) = \{(q, p) \in \mathbb{R}^{2n} \mid 0 = \text{Prox}_{hG}(q + hp), q = 0\} \\ = \{(q, p) \in \mathbb{R}^{2n} \mid q = 0, p \in h\partial P(0)\} \qquad (66)$$

and

$$T^{-2}(0) = \{(q, p) \in \mathbb{R}^{2n} \mid 0 = \text{Prox}_{hG}(q + hp), \\ -q \in h^2 \partial P(0)\} \\ = \{(q, p) \in \mathbb{R}^{2n} \mid q + hp \in h^2 \partial P(0), \\ -q \in h^2 \partial P(0)\}.$$

From (61), there exists $\varepsilon > 0$ such that $T^{-2}(0) \supset \varepsilon\mathcal{B} \times \varepsilon\mathcal{B}$. Thus, the origin is finite-time stable by Proposition 2. ■

**Remark 7** The inequality in (64) is a dissipation inequality: the system with input $\xi_{2,k+1} = \partial\hat{g}(q_{k+1})$, output $-\xi_{1,k+1}$, is passive w.r.t. the supply rate $\langle -h\xi_{1,k+1}, \xi_{2,k+1} \rangle$ and storage function $H(q_k, p_k)$ [7, section 3.15].

**Example 7 (Super-Twisting Algorithm)** The algorithm satisfies Assumption 4 with $P(q) = \gamma_2|q|$ and $\hat{g}(q) = \gamma_1 \frac{2}{3}|q|^{3/2}$. Indeed, the potential function $P$ has a unique minimum at $q = 0$, and we have $\partial\hat{g}(0) = 0$ and $\partial\hat{g}(q) \cdot \partial P(q) = \gamma_1\gamma_2|q|^{1/2}$, that is, it satisfies (62). The latter also implies (12). Thus, the iteration satisfying

$$q_{k+1} = q_k + hp_{k+1} - h\gamma_1 \lfloor q_{k+1} \rfloor^{1/2} \qquad (67a)$$

$$p_{k+1} \in p_k - h\gamma_2 \operatorname{sgn}(q_{k+1}) \qquad (67b)$$

renders the origin finite-time stable. The explicit expression associated with (67a) is equal to [9]

$$q_{k+1} = q_k + hp_k - h\gamma_1\xi_{1,k+1} - h^2\gamma_2\xi_{2,k+1} \qquad (68a)$$

$$p_{k+1} = p_k - h\gamma_2\xi_{2,k+1}, \qquad (68b)$$

where

$$\xi_{1,k+1} = \tilde{\beta}(q_k + hp_k) \operatorname{Proj} \left( \frac{q_k + hp_k}{h^2\gamma_2}, [-1, 1] \right) \\ = \lfloor q_{k+1} \rfloor^{\frac{1}{2}} \qquad (69a)$$

$$\xi_{2,k+1} = \operatorname{Proj} \left( \frac{q_k + hp_k}{h^2\gamma_2}, [-1, 1] \right) \in \operatorname{sgn}(q_{k+1}) \qquad (69b)$$

and $\tilde{\beta} : \mathbb{R} \to \mathbb{R}_+$ is

$$\tilde{\beta}(s) = -\frac{h\gamma_1}{2} + \frac{1}{2}\sqrt{h^2\gamma_1^2 + 4\max\{0, |s| - h^2\gamma_2\}}. \qquad (69c)$$

The discrete-time controller is identified from (68a) as $u_k = p_k - \gamma_1\xi_{1,k+1} - h\gamma_2\xi_{2,k+1}$, where $p_k, \xi_{1,k+1}$, and $\xi_{2,k+1}$ are given by (68b)-(69b). △

The results in Example 7 are similar to those in [9, Corollaries 1, 2], however the proof is quite different.

## 5 Open issues

* Relaxing Assumptions 2 and 4: In the general case of models (7) and (11), the set-valued nature of both the Hamiltonian and the dissipation function poses difficulties regarding the well-posedness of the discretised closed-loop. Consider for instance the model (7) with $R(x) = R(x)^\top > 0$ a constant diagonal matrix. In the planar case, the backward-Euler discretisation gives the implicit recursion

$$x_{k+1} \in x_k - h \begin{pmatrix} R_{1,1} & -1 \\ 1 & R_{2,2} \end{pmatrix} \partial H(x_{k+1}). \qquad (70)$$
that the same holds for multivalued Hamiltonian systems with dissipation, as finite-time stability is also well-preserved. Finally, splitting schemes that have been traditionally applied to subgradient systems can be effectively adapted to the multivalued Hamiltonian framework.

Future research could focus on extensions towards port-Hamiltonian systems and Hamiltonian systems with disturbances and uncertainties.

## References

[1] J. Alvarez, I. Orlov, and L. Acho. An invariance principle for discontinuous dynamic systems with application to a Coulomb friction oscillator. *Journal of Dynamic Systems, Measurement, and Control*, 122(4):687–690, 02 2000.

[2] S. Aoues, M. Di Loreto, D. Eberard, and W. Marquis-Favre. Hamiltonian systems discrete-time approximation: Losslessness, passivity and composability. *Systems & Control Letters*, 110:9–14, 2017.

[3] V.I. Arnold. *Mathematical Methods of Classical Mechanics*. Springer-Verlag, New York, 1989.

[4] A. Bacciotti and F. Ceragioli. Stability and stabilization of discontinuous systems and nonsmooth Lyapunov functions. *ESAIM: COCV*, 4:361–376, 1999.

[5] A. Bacciotti and L. Rosier. *Liapunov Functions and Stability in Control Theory*. Springer-Verlag, The Netherlands, 2005.

[6] H. H. Bauschke and P. L. Combettes. *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*. CMS Books in Mathematics. Springer, 2011.

[7] H.W. Broer and F. Takens. Chapter 1 - preliminaries of dynamical systems theory. volume 3 of *Handbook of Dynamical Systems*, pages 1-42. Elsevier Science, 2010.

[8] B. Brogliato. *Nonsmooth Mechanics. Models, Dynamics, and Control*. Communications and Control Eng. Springer Int. Pub. Switzerland, 3rd edition, 2016.

[9] B. Brogliato, B. Maschke, R. Lozano, and O. Egeland. *Dissipative Systems Analysis and Control: Theory and Applications*. Communications and Control Eng. Springer, London, 3rd edition, 2020.

[10] B. Brogliato, A. Polyakov, and D. Efimov. The implicit discretization of the supertwisting sliding-mode control algorithm. *IEEE Transactions on Automatic Control*, 65(8):3707–3713, August 2020.

[11] R. S. Burachik and A. N. Iusem. *Set-valued Mappings and Enlargements of Monotone Operators*. Springer, 2008.

[12] F. Castaños, F. Miranda Villatoro, and B. Brogliato. When proximal-point algorithms meet set-valued systems. an optimization point of view of discrete-time sliding modes. Research Report hal-04362282, INRIA, Grenoble, France, 2023. https://dumas.ccsd.cnrs.fr/LJK-MAD-TRIPOP/hal-04362282v1.

[13] F.H. Clarke. Periodic solutions to Hamiltonian inclusions. *Journal of Differential Equations*, 40(1):1–6, 1981.

[14] D. Eberard, B. Maschke, and A.J. van der Schaft. An extension of Hamiltonian systems to the thermodynamic phase space: Towards a geometry of nonreversible processes. *Rep. on Mathematical Physics*, 60:175–198, October 2007.

[15] M. Frémond. *Collisions*. Facolta di Ingegneria di Roma Tor Vergata, Roma, Italy, 2007.

[16] L. Gören-Sümer and Y. Yalçin. Gradient based discrete-time modeling and control of Hamiltonian systems. *IFAC Proceedings Volumes*, 41(2):212–217, 2008. 17th IFAC World Congress.

[17] E. Hairer, C. Lubich, and G. Wanner. *Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations*. Springer, second edition, 2006.

[18] J.B. Hiriart-Urruty and C. Lemaréchal. *Fundamentals of Convex Analysis*. Grundlehren Text Editions. Springer-Verlag Berlin Heidelberg, 2001.

[19] O. Huber, V. Acary, and B. Brogliato. Lyapunov stability analysis of the implicit discrete-time twisting control algorithm. *IEEE Transactions on Automatic Control*, 65(6):2619–2626, June 2020.

[20] D. S. Laila and A. Astolfi. Construction of discrete-time models for port-controlled Hamiltonian systems with applications. *Systems & Control Letters*, 55(8):673–680, 2006. New Trends in Nonlinear Control.

[21] X.A. Le. *Dynamics of Mechanical Systems with Coulomb Friction*. Foundations of Engineering Mechanics. Springer, Berlin, 2003.

[22] A. Levant. Higher-order sliding modes, differentiation and output-feedback control. *Int. Journal of Control*, 76:924–941, 1993.

[23] A. Levant. Sliding order and sliding accuracy in sliding mode control. *Int. J. Control*, 58:1247–1263, 1993.

[24] A. Levant. Higher-order sliding modes, differentiation and output-feedback control. *International Journal of Control*, 76(9-10):924–941, 2003.

[25] L.V. Levantovsky. Second order sliding algorithms. Their realization. *Dynamics of Heterogeneous Systems*, pages 32–43, 1985. Moscow: Institute for System Studies (in Russian).

[26] P. L. Lions and B. Mercier. Splitting algorithms for the sum of two nonlinear operators. *SIAM Journal of Numerical Analysis*, 16(6):964–979, 1979.

[27] A. Macchelli. Control design for a class of discrete-time port-Hamiltonian systems. *IEEE Transactions on Automatic Control*, 68(12):8224–8231, 2023.

[28] B. Maschke and A.J. van der Schaft. Port-controlled Hamiltonian systems: modelling origins and system-theoretic properties. In *Proc. IFAC Symposium on Nonlinear Control Systems*, pages 359 – 365, Bordeaux, France, June 1992.

[29] B. Maschke, A.J. van der Schaft, and P.C. Breedveld. An intrinsic Hamiltonian formulation of the dynamics of LC-circuits. *IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications*, 42:73–82, feb 1995.

[30] F. Miranda-Villatoro. A variational approach to the design of discrete-time super-twisting-like algorithms. https://inria.hal.science/hal-04511213v1, March 2024.

[31] S. Monaco, D. Normand-Cyrot, M. Mattioni, and A. Moreschini. Nonlinear Hamiltonian systems under sampling. *IEEE Transactions on Automatic Control*, 67(9):4598–4613, 2022.

[32] Y. Nutku. Hamiltonian structure of the Lotka-Volterra equations. *Physics Letters A*, 145(1):27–28, 1990.

[33] R.T. Rockafellar. *Convex Analysis*. Princeton University Press, 1970.

[34] R.T. Rockafellar and J.B. Wets. *Variational Analysis*, volume 317 of *Grundlehren der mathematischen Wissenschaften*. Springer Verlag, 3rd printing edition, 1998.

[35] E.K. Ryu and W. Yin. *Large-Scale Convex Optimization: Algorithms & Analyses via Monotone Operators*. Cambridge University Press, 2022.
Since the difference $R - J$ is non-symmetric, there is no direct way of computing the next state $x_{k+1}$ in the general case. Indeed, since the operator $(R - J)\partial H$ need not be monotone (recall that the product of monotone operators is not necessarily monotone), the existence and uniqueness of $x_{k+1}$ is not guaranteed. Nevertheless, the are cases where the structure of $\partial H$ can be exploited, as is the case for the backward-Euler discretisation of Example 2:

$$q_{k+1} = q_k + h \frac{\beta^2 - \beta + 1}{\beta} p_{k+1} - h \frac{\beta}{L} \xi_{k+1} \quad (71a)$$

$$p_{k+1} = p_k - h \frac{\beta}{L} \xi_{k+1} - h(1-\beta)p_{k+1} \quad (71b)$$

$$\xi_{k+1} \in \partial f(q_{k+1}). \quad (71c)$$

After simple manipulations and computations, it follows that (71) is equivalent to

$$q_{k+1} = \text{Prox}_{h\alpha_1 H}(q_k + h\alpha_2 p_k)$$

$$p_{k+1} = \frac{1}{1+h(1-\beta)} p_k - h \frac{\beta}{L} (I_d - \text{Prox}_{h\alpha_1 H})(q_k + h\alpha_2 p_k),$$

where $\alpha_1 = \frac{\beta}{L} + \frac{h(\beta^2 - \beta + 1)}{L(1+h(1-\beta))}$ and $\alpha_2 = \frac{\beta^2 - \beta + 1}{\beta(1+h(1-\beta))}$.

* Different discretisations: Backward-Euler is efficient in solving several issues like digital chattering, existence of a discrete-time sliding-surface, preservation of finite-time convergence, even in the case of multi-sliding surfaces [8,22]. Could interesting features be added to these methods without destroying such properties? Semi-implicit and $\theta$-methods could be tackled, especially in cases with $J(x)$ and $R(x)$ as in (7) or (16).

* Port-Hamiltonian systems: study of interconnected backward-Euler discretised set-valued Hamiltonian systems (preservation of finite-time stability).

## 6 Conclusions

This article is largely concerned with set-valued Hamiltonian systems with set-valued dissipation that exhibit sliding modes. They make particular differential inclusions. Applications can be found in sliding-mode controllers (twisting and super-twisting algorithms, nested controllers), mechanical systems with multi-valued Coulomb's friction, dynamical optimisation (Nesterov gradient acceleration). The first part is dedicated to detailed analysis of the finite-time stability in continuous-time. The second part is devoted to the analysis in discrete-time, when the backward-Euler (or, implicit-Euler) method is applied. Finite-time stability is studied, as well as implementation issues using splitting algorithms. Numerical simulations illustrate theoretical findings.

## References

[1] S. Aoues, M. Di Loreto, D. Eberard, and W. Marquis-Favre. Hamiltonian systems discrete-time approximation: Losslessness, passivity and composability. Systems & Control Letters, 110:9–14, 2017.

[2] V.I. Arnold. Mathematical Methods of Classical Mechanics. Springer-Verlag, New York, 1989.

[3] A. Bacciotti and F. Ceragioli. Stability and stabilization of discontinuous systems and nonsmooth Lyapunov functions. ESAIM: COCV, 4:361–376, 1999.

[4] A. Bacciotti and L. Rosier. Liapunov Functions and Stability in Control Theory. Springer-Verlag, The Netherlands, 2005.

[5] H. H. Bauschke and P. L. Combettes. Convex Analysis and Monotone Operator Theory in Hilbert Spaces. CMS Books in Mathematics. Springer, 2011.

[6] B. Brogliato. Nonsmooth Mechanics. Models, Dynamics, and Control. Communications and Control Eng. Springer Int. Pub. Switzerland, 3rd edition, 2016.

[7] B. Brogliato, B. Maschke, R. Lozano, and O. Egeland. Dissipative Systems Analysis and Control: Theory and Applications. Springer, London, 3rd edition, 2020.

[8] B. Brogliato and A. Polyakov. Digital implementation of sliding-mode control via the implicit method: A tutorial. Int. J. Robust and Nonlinear Contr., 31(9):3528–3586, 2021.

[9] B. Brogliato, A. Polyakov, and D. Efimov. The implicit discretization of the supertwisting sliding-mode control algorithm. IEEE Trans. Aut. Contr., 65(8):3707–3713, August 2020.

[10] R. S. Burachik and A. N. Iusem. Set-valued Mappings and Enlargements of Monotone Operators. Springer, 2008.

[11] F.H. Clarke. Periodic solutions to Hamiltonian inclusions. Journal of Differential Equations, 40(1):1–6, 1981.

[12] J. Cortes. Discontinuous dynamical systems. IEEE Control Systems Magazine, 28(3):36–73, 2008.

[13] D. Eberard, B. Maschke, and A.J. van der Schaft. An extension of Hamiltonian systems to the thermodynamic phase space: Towards a geometry of nonreversible processes. Rep. on Mathematical Physics, 60:175–198, October 2007.

[14] J. Eckstein and D. P. Bertsekas. On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators. Math. Prog., 55:293–318, 1992.

[15] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration: Structure-preserving algorithms for ordinary differential equations. Springer, second edition, 2006.

[16] J.B. Hiriart-Urruty and C. Lemaréchal. Fundamentals of Convex Analysis. Grundlehren Text Editions. Springer-Verlag Berlin Heidelberg, 2001.

[17] D. S. Laila and A. Astolfi. Construction of discrete-time models for port-controlled Hamiltonian systems with applications. Systems & Control Letters, 55(8):673–680, 2006.

[18] A. Levant. Sliding order and sliding accuracy in sliding mode control. Int. Journal of Control, 58:1247–1263, 1993.

[19] A. Levant and A. Michael. Adjustment of high-order sliding-mode controllers. Int. J. Robust Nonlinear Control, 19:1657–1672, 2009.

[20] L.V. Levantovsky. Second order sliding algorithms. Their realization. Dynamics of Heterogeneous Systems, pages 32–43, 1985. Moscow: Institute for System Studies (in Russian).

[21] P.L. Lions and B. Mercier. Splitting algorithms for the sum of two nonlinear operators. SIAM Journal of Numerical Analysis, 16(6):964–979, 1979.
[36] S.P. Sethi. *Optimal Control Theory: Applications to Management Science and Economics*. Springer New York, NY.

[37] S. Stramigioli, C. Secchi, A.J. van der Schaft, and C. Fantuzzi. Sampled data systems passivity and discrete port-hamiltonian systems. *IEEE Transactions on Robotics*, 21(4):574–587, 2005.

[38] V. Talasila, J. Clemente-Gallardo, and A.J. van der Schaft. Discrete port-hamiltonian systems. *Systems & Control Letters*, 55(6):478–486, 2006.

[39] A.J. van der Schaft. *L<sub>2</sub>-Gain and Passivity Techniques in Nonlinear Control*. Springer-Verlag, London, 2017.
[22] F.A. Miranda-Villatoro. A variational approach to the design of multivariable discrete-time super-twisting-like algorithms. *IEEE Transactions on Automatic Control*, pages 1-12, 2024.

[23] S. Monaco, D. Normand-Cyrot, M. Mattioni, and A. Moreschini. Nonlinear Hamiltonian systems under sampling. *IEEE Trans. Aut. Contr.*, 67(9):4598-4613, 2022.

[24] J.A. Moreno. *Lyapunov Approach for Analysis and Design of Second Order Sliding Mode Algorithms*, pages 113-149. Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.

[25] M. Muehlebach and M. I. Jordan. A dynamical systems perspective on Nesterov acceleration. In *Proc. 36th Int. Conference on Machine Learning*, pages 4656-4662, 2019.

[26] Y. Nutku. Hamiltonian structure of the Lotka-Volterra equations. *Physics Letters A*, 145(1):27-28, 1990.

[27] Y. Orlov. *Discontinuous Systems, Lyapunov Analysis and Robust Synthesis under Uncertainty Conditions*. Springer-Verlag, London, 2009.

[28] R. Ortega, A.J. van der Schaft, B. Maschke, and G. Escobar. Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems. *Automatica*, 38:585–596, April 2002.

[29] R.T. Rockafellar and J.B. Wets. *Variational Analysis*, volume 317 of *Grundlehren der mathematischen Wissenschaften*. Springer Verlag, 3rd printing edition, 1998.

[30] E.K. Ryu and W. Yin. *Large-Scale Convex Optimization: Algorithms & Analyses via Monotone Operators*. Cambridge University Press, 2022.

[31] S.P. Sethi. *Optimal Control Theory: Applications to Management Science and Economics*. Springer New York, NY, 2019.

[32] V. Talasila, J. Clemente-Gallardo, and A.J. van der Schaft. Discrete port-hamiltonian systems. *Systems & Control Letters*, 55(6):478-486, 2006.

[33] A.J. van der Schaft. *L<sub>2</sub>-Gain and Passivity Techniques in Nonlinear Control*. Springer-Verlag, London, 2017.