# Backstepping passivity-based trajectory tracking control of frictional oscillators: continuous and discrete time analyses

Bernard Brogliato, Félix Miranda-Villatoro, Aya Younes

► To cite this version:

Bernard Brogliato, Félix Miranda-Villatoro, Aya Younes. Backstepping passivity-based trajectory tracking control of frictional oscillators: continuous and discrete time analyses. 2024. hal-04842497

HAL Id: hal-04842497

https://inria.hal.science/hal-04842497v1

Preprint submitted on 17 Dec 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Backstepping passivity-based trajectory tracking control of frictional oscillators: continuous and discrete time analyses

Bernard Brogliato and Félix Miranda-Villatoro and Aya Younes*

December 16, 2024

Passivity-based control, backstepping control, LMI, maximal monotone operator, trajectory tracking, set-valued Coulomb's friction, superpotential, proximal-point algorithms, sliding mode control, finite-time convergence.

## Abstract

This article is largely concerned with the design of trajectory tracking controllers for a class of frictional oscillators, represented by differential inclusions with maximal monotone operators. A backstepping passivity-based control strategy is proposed, which relies on set-valued friction and its monotonicity properties (hence the goal is not to compensate for friction but to take advantage of it). Both the continuous-time design and its discrete-time counterpart are analysed. Robustness with respect to uncertainties in the friction coefficient is carefully analyzed in both continuous and discrete time settings. Details about the set-valued controller implementation are provided. Numerical simulations illustrate the theoretical developments.

# 1 Introduction

## 1.1 Frictional oscillators

This article focuses on the control of one-degree-of-freedom oscillators, a class of systems widely studied in the Applied Mathematics, Nonlinear Dynamics and in the Bifurcation and Chaos scientific communities, see, e.g., [1, Chapter 9] [2, 3, 4, 5]. It may also be seen as a simplified Burridge-Knopoff model [6]. The behaviour of the stick-slip oscillators with various friction models has been studied see, e.g., [2, 7, 8, 9]. The existence and stability analysis of periodic solutions of a periodically excited frictional oscillator was originally

<footnote>*Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, 38000 Grenoble, France.</footnote>
presented by [4, 5], see also [10, 3]. Externally excited frictional oscillators, as studied by
[11, 7, 8], may exhibit chaotic behaviour. They undergo rich bifurcational behaviour, see,
e.g., [12, 13, 2, 8]. Stick-slip vibrations, which are self-sustained oscillations induced by dry
friction are described in [11, 10]. Control techniques have been analyzed and implemented,
such as active feedback control to superimpose frictional force in [14], and control strategies
based on delayed time feedback in [15, 16].

1.2 Friction Modeling

The analysis and the control of systems with friction have witnessed a large number of articles
in both Automatic Control and the Nonlinear Dynamics. Roughly speaking, friction models
of interest for Control may be classified as follows: static models (which may be set-valued or
single-valued, with constant or with varying friction coefficient), and dynamic models which
incorporate pre-sliding, bristle, hysteretic effects, etc. Static models comprise Coulomb's
friction with constant coefficient or varying coefficient (e.g., Stribeck [17], Dieterich-Ruina
[18, Equation (21)]), which are set-valued at zero relative tangential velocity, as well as
regularized models (for planar friction, where the set-valued signum function is replaced by
a single-valued saturation, or a sigmoid function). Set-valuedness at zero tangential velocity
allows to correctly handle sticking modes. It is thus crucial in multibody applications.
Dynamic models with internal state are numerous (LuGre, Leuven, Bliman-Sorine, Dahl,
elastoplastic etc [19]). They are better suited for modeling micro stick-slip effects, though
Coulomb with Stribeck coefficient of friction also provides good results in this respect [19].
Dynamic models also present severe drawbacks [20, 21]. Coulomb friction (with constant
coreffunctions) has limitations [18, 8], but it can also provide quite good results in
many instances, see, e.g., [22, 23, 8]. It is sometimes not an easy task to choose between a
regularized (around zero tangential velocity) model and a set-valued model, see, e.g., [24]. It
is inferred that conducting preliminary research work using the set-valued Coulomb's model
(with constant coefficient of friction) is a reasonable path for trajectory tracking, especially
if robustness with respect to uncertain coefficients is shown.

1.3 Control with Friction

The subject of friction in controlled systems is wide. It has mainly concerned friction com-
pensation at low tangential velocity, see, e.g., [25] and references therein, or stability/sta-
bilization of equilibrium points or of sets of equilibria [17, 26]. To the best of the authors'
knowledge, the problem that is tackled in this article is original, because it is based on using
friction to achieve trajectory tracking with a backstepping passivity-based input.
<figure><img src="image_2.png" /><figcaption>Figure 1: Frictional oscillator with Coulomb’s friction.</figcaption></figure>

## 1.4 Controlled Frictional Oscillator’s Dynamics

Consider the frictional oscillator system depicted in Fig. 1. The objective is to perform trajectory tracking for the mass with position $q$. Neglecting the frictional effects on the belt dynamics so that the pulleys-belt’s velocity can be considered not influenced by the friction with the mass, is a common assumption in the Nonlinear Dynamics literature [2, 3, 4, 5]. If friction on the pulley-belt dynamics is not neglected, because it does not hold that $\frac{r\mu mg}{I} \ll 1$, the following set-valued dynamics is obtained:

$$m\ddot{q}(t) = -kq(t) - f\dot{q}(t) + u_1(t) - \mu mg\lambda_t(t) \quad (1a)$$

$$I\ddot{\theta}(t) = \tau(t) + r\mu mg\lambda_t(t) \quad (1b)$$

$$\lambda_t(t) \in \operatorname{sgn}(\dot{q}(t) - r\dot{\theta}(t)) \quad (1c)$$

where $\tau$ is the control torque applied on the pulleys, $I$ is the inertia momentum of each pulley, $\dot{\theta} = \omega$, $\mu > 0$ is the friction coefficient, $r$ is their radius, $g$ is the gravity acceleration, $m$ kg the mass weight, $k \ge 0$ is a stiffness, $f \ge 0$ is a viscous friction coefficient, $-\mu mg\lambda_t$ is the tangential force exerted by the belt on the mass (thus, $\lambda_t$ is a selection of the set-valued signum function), $q$ is the mass’ position, $\dot{q}$ its velocity, $\ddot{q}$ its acceleration, $\theta$ is the pulley’s angle, $\dot{\theta}$ its angular velocity, $\ddot{\theta}$ its angular acceleration, $\operatorname{sgn}(x) = 1$ if $x > 0$, $-1$ if $x < 0$, $[-1, 1]$ if $x = 0$. It may also be assumed that the spring is acted upon by the input $u_1$ (which is a displacement instead of a force) instead of the mass [8]. Then the term $-kq$ in (1a) is replaced by $-k(q + u_1)$. Provided that $k$ is known, this does not change the control problem tackled in this article.
## 1.5 Backstepping Passivity-based Tracking Control

Let $u_1 = u_1(q, \dot{q}, t)$. Then a remarkable feature of (1) is that the only link between (1a) and (1b) is made through the selection $\lambda_t$ in (1c). The basic idea in this article is to use a backstepping-like strategy to give $\lambda_t$ a desired value which allows to control (1a), using a suitable $\tau$ in (1b). To that aim we first recall how to achieve trajectory tracking when $r\dot{\theta}$ is considered as a control input $u_2$, using $u = (u_1, u_2)^\top$. This is based on the results in [27]. Then a fictitious input $\vartheta_2$ is used, and the error input $\tilde{u}_2$ is forced to be finite-time stable. The maximal monotonicity of the signum set-valued function is used to show stability, in a similar way as in [28, 27]. It is noteworthy that this can be seen as an extension of [29, 30] where backstepping is applied to mechanical systems where the coupling between both (actuated, and non-actuated) dynamics is a force that derives from the elasticity potential energy, which is a quadratic function of the positions. In (1) the coupling $\lambda_t$ derives from a nonsmooth potential $x \mapsto |x|$ with corresponding set-valued force $\partial|x| = \text{sgn}(x)$, where $\partial$ is the subdifferential of Convex Analysis [31]. This is named a superpotential or pseudopotential [32].

To the best of authors' knowledge, this work is the first application of backstepping to
mechanical systems with superpotentials set-valued interconnections.

## 1.6 Article's Outline

This article is organized as follows: section 2 recalls the results obtained when the belt's velocity is considered as the control input, as a consequence of the material in [27]; section 3 is dedicated to present the backstepping controller and its stability analysis; section 4 studies robustness of the backstepping scheme with respect to uncertain friction coefficients; the discrete-time controller, its implementation, stability and robustness are analysed in detail in section 5; numerical simulations illustrate the theoretical findings in section 6; conclusions end the article in section 7.

**Notation and Definitions** A quadruple $(A, B, C, D)$ is passive if there exists $P = P^\top \succcurlyeq 0$ such that $$ \begin{pmatrix} -A^\top P - PA & PB - C^\top \\ B^\top P - C & D + D^\top \end{pmatrix} \succcurlyeq 0 $$ (this is called the passivity matrix inequality PMI$(A, B, C, D$)). Equivalently, the system $\dot{x} = Ax + Bu$, $y = Cx + Du$ is passive if there exists a storage function $V(x) = \frac{1}{2}x^\top Px$ such that for any $t_2 \ge t_1$, $V(x(t_2)) - V(x(t_1)) \le \int_{t_1}^{t_2} u(t)^\top y(t)dt$. The system is strictly state passive if the PMI holds with $-A^\top P - PA - \varepsilon P$ for some $\varepsilon > 0$ (then we speak of the SSPMI$(A, B, C, D, \varepsilon$)), and $P \succcurlyeq 0$. It is strongly passive if the PMI$(A, B, C, D)$ is satisfied with strict inequalities. Let $M \in \mathbb{R}^{n \times m}$, then $M_{i, \bullet}$ is its $i$-th row and $M_{\bullet, i}$ its $i$-th column. $M \in \mathbb{R}^{n \times n}$ is a P-matrix if all its principal minors are positive. Let $K \subseteq \mathbb{R}^n$ be a closed nonempty convex set. Its normal cone at $x \in K$ is $N_K(x) = \{z \in \mathbb{R}^n | z^\top(y-x) \le 0 \ \forall y \in K\}$. Let $\mathcal{M}: \mathbb{R}^n \to \mathbb{R}^n$ be a maximal
monotone operator. Its resolvent of degree $\alpha$ is $J_{\mathcal{M}}^{\alpha} = (I_d + \alpha \mathcal{M})^{-1}$, $\alpha > 0$. It is Lipschitz continuous, single-valued and non-expansive. The $L^{\infty}$-norm of a signal $u: \mathbb{R}_+ \to \mathbb{R}$ is given by $\|u\|_{L^{\infty}} = \ess \sup \{u(t)|t \in [0, +\infty)\}$. Let $K \subseteq \mathbb{R}^n$ be a set, the orthogonal projection of $x \in \mathbb{R}^n$ onto $K$ is denoted by $\text{proj}_K(x)$. AC stands for absolutely continuous. GUUB stands for globally uniformly ultimately bounded.

## 2 Direct Controller Design ($u_2 = r\dot{\theta}$)

First let us assume that the pulley's angular velocity can be given an arbitrary value, hence the controlled system reduces to (1a), (1c) with control input $u = (u_1, u_2)^{\top}$, $u_2 = r\dot{\theta}$. The objective is to track trajectories $(q_d(t), \dot{q}_d(t))^{\top}$. Following the passivity-based strategy in [27], desired trajectories are generated by desired dynamics. The goal is not to compensate for $\lambda_t$ in (1a) using a sliding-mode input $u_1$, but to use $u_2$ through friction coupled to the maximal monotonicity of the signum multifunction (similar to the maximal monotonicity of the normal cone mapping in [27]). To that aim we set:

$$u = Ke + u_d(t) \qquad (2)$$

with $e = (\tilde{q}, \dot{\tilde{q}})^{\top}$, $\tilde{q} = q - q_d(t)$, $K = \begin{pmatrix} k_{11} & k_{12} \\ k_{21} & k_{22} \end{pmatrix}$, and

$$\left\{ m\ddot{q}_d(t) = -kq_d(t) - f\dot{q}_d(t) + u_{1d}(t) - \mu mg\lambda_{t,d}(t) \right. \qquad (3a)$$

$$\left. \lambda_{t,d}(t) \in \text{sgn}(\dot{q}_d(t) - u_{2d}(t)) \right\} \qquad (3b)$$

is the desired dynamics which generates the desired trajectories (i.e., $u_d(t)$ has to be chosen suitably to generate desired solutions, see [33, 34] for efficient solutions to compute periodic solutions). Let us define the following matrices: $A_{cl} = \begin{pmatrix} 0 & 1 \\ \frac{k_{11}-k}{m} & \frac{k_{12}-f}{m} \end{pmatrix}$, $B_{cl} = \begin{pmatrix} 0 & 0 \\ -1 & 0 \end{pmatrix}$, $C_{cl} = \begin{pmatrix} k_{21} & -1+k_{22} \\ 0 & 0 \end{pmatrix}$, $D_{cl} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$, $\lambda = (\lambda_1 \ \lambda_2)^{\top}$, $\lambda_d = (\lambda_{1d} \ \lambda_{2d})^{\top}$. The closed-loop error system (1a), (1c), (2), (3) writes as the linear complementarity system (LCS):

$$\left\{ \begin{array}{l} \dot{e} = A_{cl}e + 2\mu g B_{cl}(\lambda - \lambda_d(t)) \\ 0 \le \lambda \perp w = C_{cl}e + D_{cl}\lambda + \begin{pmatrix} -\dot{q}_d(t) + u_{2d}(t) \\ 1 \end{pmatrix} \ge 0 \\ 0 \le \lambda_d \perp w_d = D_{cl}\lambda_d + \begin{pmatrix} -\dot{q}_d(t) + u_{2d}(t) \\ 1 \end{pmatrix} \ge 0, \end{array} \right. \qquad (4)$$

where the representation of the signum set-valued function in a complementarity framework is used [35, 36], see Appendix A. Hence (4) is equivalent to the differential inclusion:

$$\dot{e}(t) \in A_{cl}e(t) + \mu g E_{cl} \text{sgn}(\dot{q}_d(t) - u_{2d}(t)) - \mu g E_{cl} \text{sgn}(-[C_{cl}]_{1,\bullet}e(t) + \dot{q}_d(t) - u_{2d}(t)), \qquad (5)$$
where $\lambda_1 = \frac{\lambda_t+1}{2} \in [0, 1]$, and $E_{cl} = (0 \ 1)^T$. We have $\lambda_1 = \frac{\lambda_t+1}{2}$, with $\lambda_t$ in (1c). It is noteworthy that $\mu$ does not play any role in the complementarity constraints, which equivalently represent (1c) and (3b). Still following the developments in [27] (which deal with LCS control) the following holds true:

**Proposition 1** *Consider the error dynamics in (4). Assume that the quadruple $(A_{cl}, B_{cl}, C_{cl}, D_{cl})$ is strictly state passive. Then $e(t) \to 0$ as $t \to +\infty$ exponentially fast for any initial condition.*

The SSPMI $(A_{cl}, B_{cl}, C_{cl}, D_{cl}, \epsilon)$ with unknowns $K$ and $P$ [27, Section 3, Appendix C] is used to calculate the gains. Numerical examples prove that the PMI is solvable, with $-2p_{12}\mu g = k_{21}$ and $-2p_{22}\mu g = k_{22} - 1 < 0$. The proof uses the monotonicity of the set-valued part, and is led similarly as [27, Proposition 3.8]. It does not use the boundedness of the signum function to compensate for it. The oscillator block diagram is depicted in Fig. 2. In addition, some robustness results can be obtained, which are a particular case of [27, section 4] since the matrix $P$ is computed independently of the friction coefficient $\mu$.

<figure><img src="image_3.png" /><figcaption>Figure 2: Block diagram for (1a), (1c) (2) (two-channel set-valued system). Same applies to (3) with $u_{1d}$ and $u_{2d}$.</figcaption></figure>

**Remark 1** A summary of LCS well-posedness is proposed in [27, section 2]. Passivity is used to prove that (4) possesses AC solutions with uniqueness. The desired dynamics always has an AC solution imposing mild regularity conditions on $u_{1d}$ and $u_{2d}$ (relying for instance on Filippov's differential inclusions). Uniqueness may be obtained using LCS well-posedness tools.

When $k > 0$ and $f > 0$, it is possible to solve the problem with $k_{11} = k_{12} = 0$ (provided that the PMI has a solution). Then $u_1$ reduces to a feedforward input $u_1(q_d, \dot{q}_d, u_{1d})$, and the control problem is that of a robot-object system [37].
# 3 Backstepping Controller Design with Pulley-Belt's Dynamics

Consider the mass dynamics in (1a). For the sake of notation simplicity we assume without loss of generality that $k = 0$ and $f = 0$ (or, we just replace formally $k_{11} - k$ by $k_{11}$ and $k_{12} - f$ by $k_{12}$).

## 3.1 Controller Design

Define $\vartheta_2 \triangleq k_{21}\tilde{q} + k_{22}\dot{\tilde{q}} + u_{2d}$, then (1a) is equivalently rewritten as:

$$ \ddot{q}(t) \in \frac{u_1(t)}{m} - \mu g \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t) - \tilde{u}_2(t)) \quad (6) $$

with $\tilde{u}_2 \triangleq r\dot{\theta} - \vartheta_2$. Here, $\vartheta_2$ plays the role of $u_2$ in section 2. Thus using (3) and $u_1(q, \dot{q}, t)$ from (2), the error dynamics is:

$$ \ddot{q}(t) = \frac{1}{m}(k_{11}\tilde{q}(t) + k_{12}\dot{\tilde{q}}(t)) - \mu g (\lambda_t(t) - \lambda_{t,d}), \quad (7) $$

where $\lambda_{t,d}(t) \in \operatorname{sgn}(\dot{q}_d(t) - u_{2d}(t))$ and $\lambda_t(t) \in \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t) - \tilde{u}_2(t))$. It is seen that (7) and (1b) (1c) has the required structure for backstepping control. The objective is to design the controller $\tau(\cdot)$ such that $\tilde{u}_2 \to 0$, so that the control problem reduces to that of section 2 (and (7) becomes equal to (5)). Thus $\dot{\tilde{u}}_2 = r\ddot{\theta} - \dot{\vartheta}_2$, i.e.:

$$ \dot{\tilde{u}}_2(t) \in \frac{r}{I}\tau(t) + \frac{r^2\mu g}{I} \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t) - \tilde{u}_2(t)) - k_{21}\dot{\tilde{q}}(t) - k_{22}\ddot{\tilde{q}}(t) - \dot{u}_{2d}(t), \quad (8) $$

where $u_{2d}(\cdot)$ is assumed to be differentiable almost everywhere. Inserting (6) into (8) yields the differential inclusion:

$$ \begin{aligned} \dot{\tilde{u}}_2(t) \in \frac{r}{I}\tau(t) & - \frac{k_{22}k_{11}}{m}\tilde{q}(t) - \left(k_{21} + \frac{k_{22}k_{12}}{m}\right)\dot{\tilde{q}}(t) + \left(\frac{r^2m}{I} + k_{22}\right)\mu g \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t) - \tilde{u}_2(t)) \\ & - \frac{k_{22}}{m}u_{1d}(t) - \dot{u}_{2d}(t) + k_{22}\ddot{q}_d(t). \end{aligned} \quad (9) $$

Consequently, it is possible to rewrite (9) as:

$$ \begin{aligned} \dot{\tilde{u}}_2(t) \in \frac{r}{I}\tau(t) & - \Lambda(q(t), \dot{q}(t), q_d(t), \dot{q}_d(t), t) \\ & + \left(\frac{r^2m}{I} + k_{22}\right) \mu g \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t) - \tilde{u}_2(t)), \end{aligned} \quad (10) $$
where $\Lambda(q, \dot{q}, q_d, \dot{q}_d, t) = k_{21}(\dot{q}(t) - \dot{q}_d(t)) + k_{22}\frac{u_1(t)}{m} - k_{22}\ddot{q}_d(t) + \dot{u}_{2d}(t)$, with $u_1 = k_{11}(q - q_d) + k_{12}(\dot{q} - \dot{q}_d) + u_{1d}$. The next step is to design an input $\tau(\cdot)$ such that $\tilde{u}_2$ converges to zero. For this purpose, let us choose $\tau(t)$ such that:

$$ \frac{r}{I}\tau(t) = \Lambda(q, \dot{q}, q_d, \dot{q}_d, t) + \tau_s(t), \quad (11) $$

with

$$ \tau_s \in -k_{sl} \operatorname{sgn}(\tilde{u}_2), \quad (12) $$

and $k_{sl} > \mu g |\frac{r^2 m}{I} + k_{22}|$ (recall that $\tilde{u}_2(q, \dot{q}, \dot{\theta}, q_d, \dot{q}_d, u_{2d})$ is calculable from state measurement). The closed-loop dynamics for $\tilde{u}_2$ is the differential inclusion:

$$ \begin{aligned} \dot{\tilde{u}}_2(t) &\in -k_{sl} \operatorname{sgn}(\tilde{u}_2(t)) + \mu g \left( \frac{r^2 m}{I} + k_{22} \right) \lambda_t(t), \\ \lambda_t(t) &\in \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t) - \tilde{u}_2(t)). \end{aligned} \quad (13) $$

Thus in this context of sliding mode control, $\lambda_t$ which is a (at least Lebesgue measurable) selection of the set-valued signum function, is considered as a bounded disturbance.

**Remark 2** The closed-loop (1b) (11) (12) reads as $I\ddot{\theta} \in \frac{I}{r}\Lambda(q, \dot{q}, t) - \frac{I}{r}k_{sl}\operatorname{sgn}(r\dot{\theta} - k_{21}\tilde{q} - k_{22}\dot{\tilde{q}} - u_{2d}(t)) + r\mu mg\lambda_t$. The proposed strategy is thus quite different from classical first-order sliding-mode control $\tau_{smc}$ of the $(\theta, \dot{\theta})$-dynamics, which would aim at compensating $r\mu mg\lambda_t$ in (1b), treating it as a disturbance, and then designing a suitable desired $\dot{\theta}_d$ to be used as a fictitious control (recall that $u_2 = r\dot{\theta}$). Indeed, letting $\sigma = \dot{\tilde{\theta}} + a\tilde{\theta}$, $\tilde{\theta} = \theta - \theta_d(t)$, $a > 0$, we get $\tau_{smc} = I\ddot{\theta}_d - Ia\dot{\tilde{\theta}} - \alpha \operatorname{sgn}(\sigma)$, $\alpha > \frac{r\mu mg}{I}$, so that $\dot{\sigma} \in -\alpha \operatorname{sgn}(\sigma) + \frac{r\mu mg}{I}\lambda_t$: $\sigma$ vanishes in finite time, and both $\tilde{\theta}$ and $\dot{\tilde{\theta}} \to 0$ exponentially globally. But this implies that $\ddot{\theta}_d$ depends on $\lambda_t$, which we do not assume to be available. Moreover the exponential convergence of $\dot{\theta}$ to a suitable $\dot{\theta}_d$ implies the presence of a disturbance inside $\lambda_t$ because $\tilde{u}_2 = r\dot{\tilde{\theta}}$ vanishes asymptotically only.

Notice that the controller $\tau$ can be computed online with measurements of $q, \dot{q}, \dot{\theta}$, since $\tilde{u}_2 = r\dot{\theta} - \vartheta_2(q, \dot{q}, u_{2d}(t))$. The dynamics for $\tilde{u}_2$ in (13) is a consequence of the choice of $\tau(q, \dot{q}, \dot{\theta}, t)$ in (11) (12), it doesn't need to be implemented.

## 3.2 Stability Analysis

The main stability result is stated.

**Proposition 2** Consider the error dynamics in (7) with bounded initial data and functions $q_d(\cdot), \dot{q}_d(\cdot), \ddot{q}_d(\cdot), u_{1d}(\cdot), u_{2d}(\cdot), \dot{u}_{2d}(\cdot)$. Let K be chosen such that $(A_{cl}, B_{cl}, C_{cl}, D_{cl})$ is strictly state passive, and $k_{sl} > \mu g |\frac{r^2 m}{I} + k_{22}|$ in (12). Then $e = (\tilde{q}, \dot{\tilde{q}})^\top$ is bounded and converges exponentially fast to zero. Moreover $\dot{\theta}$ and $\ddot{\theta}$ are bounded on $[0, +\infty)$.
**PROOF:** The proof is performed in three steps. First, the focus is on demonstrating the finite time convergence of $\tilde{u}_2$ in (13) to zero. Then, the stability of (7) is analysed when $\tilde{u}_2 = 0$. The third step demonstrates the boundedness of the variables.

• First step: Consider the dynamics of $\dot{\tilde{u}}_2(t)$ in (13). The objective is to show that $\tilde{u}_2 = 0$ for all $t \ge t^*$ where $t^* < \infty$, under the sliding mode control (12), with the sliding surface $\tilde{u}_2 = 0$. Let us take the Lyapunov function $V(\tilde{u}_2(t)) = \frac{1}{2}\tilde{u}_2^2(t)$, then along trajectories of (13):

$$ 
\begin{aligned}
\dot{V} &= \tilde{u}_2 \dot{\tilde{u}}_2 = \tilde{u}_2 \left( \tau_s + \left( \frac{r^2 \mu mg}{I} + k_{22} \mu g \right) \lambda_t \right) \\
&= -k_{sl} \tilde{u}_2 \operatorname{sgn}(\tilde{u}_2) + \tilde{u}_2 \left( \frac{r^2 \mu mg}{I} + k_{22} \mu g \right) \lambda_t
\end{aligned}
\quad (14) $$

Given that $\lambda_t \in [-1, 1]$ and $\tilde{u}_2 \operatorname{sgn}(\tilde{u}_2) = |\tilde{u}_2|$:

$$ 
\begin{aligned}
\dot{V} &= -k_{sl}|\tilde{u}_2| + \tilde{u}_2 \left( \frac{r^2 \mu mg}{I} + k_{22} \mu g \right) \lambda_t \\
&\le -k_{sl}|\tilde{u}_2| + |\tilde{u}_2| \left| \frac{r^2 \mu mg}{I} + k_{22} \mu g \right| \\
&\le |\tilde{u}_2| \left( \left| \frac{r^2 \mu mg}{I} + k_{22} \mu g \right| - k_{sl} \right)
\end{aligned}
\quad (15) $$

Considering that $k_{sl} > \left| \frac{r^2 \mu mg}{I} + k_{22} \mu g \right|$, let us define $-a \triangleq \left| \frac{r^2 \mu mg}{I} + k_{22} \mu g \right| - k_{sl} < 0$ with $a > 0$. Then: $\dot{V} \le -a|\tilde{u}_2| \le -a\sqrt{2}V^{\frac{1}{2}}$. Thus, according to [38, Theorem 4.2], the equilibrium point $\tilde{u}_2^* = 0$ of the dynamics in (13) is globally finite-time-stable. Let $t^*$ to be the finite time required for $\tilde{u}_2$ to reach zero (i.e., the solution $\tilde{u}_2(t)$ of (13) is defined on $[0, t^*]$ and $\lim_{t \to t^*} \tilde{u}_2(t) = 0$). Applying [38, Theorem 4.2] yields $t^* = \frac{2V^{\frac{1}{2}}(\tilde{u}_2(0))}{a\sqrt{2}}$.

• Second step: In this step, the stability of (7) is analysed for all $t \in [t^*, \infty)$ when $\tilde{u}_2(t) = 0$. The closed-loop system (6) is reduced to the following on $[t^*, \infty)$:

$$ 
\begin{aligned}
\ddot{q}(t) \in & \frac{k_{11}}{m} \tilde{q}(t) + \frac{k_{12}}{m} \dot{\tilde{q}}(t) + \frac{u_{1d}(t)}{m} \\
& - \mu g \operatorname{sgn} (\dot{q}(t) - k_{21} \tilde{q}(t) - k_{22} \dot{\tilde{q}}(t) - u_{2d}(t))
\end{aligned}
\quad (16) $$

Consequently, the error dynamics in (7) is reduced on the interval $[t^*, \infty)$ and is written as:

$$ 
\ddot{q}(t) = \frac{1}{m}(k_{11}\tilde{q}(t) + k_{12}\dot{\tilde{q}}(t)) + \mu g [\lambda_{t,d}(t) - \lambda_t(t)] \quad (17) $$

where $\lambda_{t,d}(t) \in \operatorname{sgn}(\dot{q}_d(t) - u_{2d}(t))$ and $\lambda_t(t) \in \operatorname{sgn}(\dot{q}(t) - \vartheta_2(t))$. Given that $\vartheta_2(t) = k_{21}(q(t) - q_d(t)) + k_{22}(\dot{q}(t) - \dot{q}_d(t)) + u_{2d}(t)$, the error dynamics in (17) are equivalent to the error dynamics in (4) or (5). Hence a similar result as in [27, Proposition 3.8] is applied in this step (it can be used to prove Proposition 1). Let us recall the proof for the sake of readability. Let $V(e) = \frac{1}{2}e^\top Pe$, where $P = P^\top > 0$ is a solution of the SSPMI($A_{cl}, B_{cl}, C_{cl}, D_{cl}, \varepsilon$). Then
along the closed-loop trajectories:

$$
\begin{align}
\dot{V}(t) &= \frac{1}{2} e^{\top} (A_{cl}^{\top} P + P A_{cl}) e + \mu g e^{\top} C_{cl}^{\top} (\lambda - \lambda_d) \nonumber \\
&= \frac{1}{2} e^{\top} (A_{cl}^{\top} P + P A_{cl}) e \nonumber \\
&\quad + \mu g (w - w_d - D_{cl} (\lambda - \lambda_d))^{\top} (\lambda - \lambda_d) \nonumber \\
&= \frac{1}{2} e^{\top} (A_{cl}^{\top} P + P A_{cl}) e + \mu g (w - w_d)^{\top} (\lambda - \lambda_d) \nonumber \\
&\le \frac{1}{2} e^{\top} (A_{cl}^{\top} P + P A_{cl}) e \le -\varepsilon e^{\top} P e. \tag{18}
\end{align}
$$

The fact that $D_{cl} = -D_{cl}^{\top} \Rightarrow PB_{cl} = C_{cl}^{\top}$, the monotonicity for each $t$ of the normal cone mapping $\mathcal{N}_{S(t)}(C_{cl}e + D_{cl}\lambda)$, $S(t) \triangleq \{\xi \in \mathbb{R}^m | \xi + Cx_d(t) + Fu_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} - C_{cl}x_d(t) \ge 0\}$, and the fact that $w(t) - w_d(t) = C_{cl}e(t) + D_{cl}(\lambda(t) - \lambda_d(t))$, are used. Using (17) it follows also that $\lambda(t) - \lambda_d(t) \to 0$.

• Third step: This step is dedicated to show the boundedness of variables. Let us analyze the dynamics in (1b). According to the analysis in section 2 and given that $\tilde{u}_2 = 0$ for $t \ge t^*$, it follows that the variables $q(\cdot)$, $\dot{q}(\cdot)$ are bounded on $[t^*, +\infty)$. Assuming that $\dot{u}_{2d}(\cdot)$, $\dot{q}_d(\cdot)$ and $\ddot{q}_{2d}(\cdot)$ are bounded, then $\ddot{q}(\cdot)$ is also bounded. Therefore, on $[t^*, +\infty)$, and given that $r\ddot{\theta} = k_{21}(\dot{q} - \dot{q}_d) + k_{22}(\ddot{q} - \ddot{q}_d) + \dot{u}_{2d}$, it follows that $\dot{\theta}(\cdot)$ and consequently $\tau(\cdot)$ are bounded. Similarly, $\dot{\theta}(\cdot)$ is bounded on $[t^*, +\infty)$. Since $\theta$ is defined modulo $2\pi$, its boundedness is a secondary issue. The boundedness of all variables on $[0, t^*)$ is ensured from the fact that there does not exist any finite escape in the closed-loop system, since all the right-hand sides are sums of linear single-valued terms and bounded set-valued terms. $\square$

**Remark 3** It follows from (18) and as a consequence of the maximal monotonicity [39, Section 4.3.3] that the system (4) is incrementally passive with supply rate $(\lambda - \lambda_d)^{\top}(w - w_d)$ and storage function $\frac{1}{2}e^{\top}Pe = \frac{1}{2}(x - x_d)^{\top}P(x - x_d)$, $x = (q, \dot{q})^{\top}$, $x_d = (q_d, \dot{q}_d)^{\top}$. Similarly to section 2, it is noticed that the backstepping scheme is still valid when $\mu = \mu(\dot{q}, \dot{\theta}, t) > 0$, provided the coefficient is known to design the desired system to obtain the error dynamics (7).

# 4 Robustness Analysis: Unknown Friction Coefficient

$\mu$

In this section it is assumed that the friction coefficient $\mu(\dot{q}, \dot{\theta}, t) > 0$ is an unknown function (which is denoted by $\mu(t)$ in the sequel to simplify notation). It is also assumed that the nominal value $\mu_0(t)$ is uniformly bounded and an upper bound on its uncertainty $\Delta\mu(t)$
(defined below) is known. Thus, in this new setting, the pulley-belt dynamics is still given by (1), but the desired dynamics is given by

$$
\begin{equation}
\begin{cases}
m\ddot{q}_d(t) = -kq_d(t) - f\dot{q}_d(t) + u_{1d}(t) \\
\qquad - \mu_0(t)mg\lambda_{t,d}(t) & (19a) \\
\lambda_{t,d}(t) \in \text{sgn}(\dot{q}_d(t) - u_{2d}(t)), & (19b)
\end{cases}
\end{equation}
$$

where $\mu_0(t) > 0$ for all $t \ge 0$. In other words, a time-varying coefficient is allowed in the desired dynamics, since setting $\mu_0(\dot{q}_d, \dot{\theta}, t)$ would destroy the independency of (19) with respect to the plant's dynamics (recall that $\dot{\theta}_d(t)$ is not defined as an arbitrary exogenous signal).

Let us analyse the behaviour of the closed-loop system (4) when $\Delta\mu(t) \stackrel{\Delta}{=} \mu(t) - \mu_0(t) \neq 0$. The error dynamics is thus described by

$$
\dot{e}(t) = A_{cl}e(t) - \mu_0(t)gE_{cl}(\lambda_t(t) - \lambda_{td}(t)) - \Delta\mu(t)gE_{cl}\lambda_{td}(t), \quad (20)
$$

where $E_{cl}$ is as in (5), $\lambda_t$ and $\lambda_{t,d}$ are as in (1c) and (19b), respectively. It is worth to notice that the last term in (20) is globally bounded, by assumption. Thus, robustness of the backstepping design, with respect to the friction coefficient, is expected.

**Proposition 3** Consider the system (1) and the desired dynamics (19), with $u_1 = k_{11}(q - q_d) + k_{21}(\dot{q} - \dot{q}_d) + u_{1d}$ and $\tau$ from the backstepping design, i.e.,

$$
\frac{r}{I}\tau \in \Lambda(q, \dot{q}, q_d, \dot{q}_d, t) - k_{sl} \operatorname{sgn}(\tilde{u}_2), \quad (21)
$$

where

$$
k_{sl} \geq g \left| \frac{r^2 m}{I} + k_{22} \right| \|\mu\|_{L^{\infty}},
$$

$$
\Lambda(q, \dot{q}, q_d, \dot{q}_d, t) = k_{21}(\dot{q} - q_d) - \frac{k_{22}}{m}u_1 - \dot{u}_{2d} + k_{22}\ddot{q}_d, \quad \tilde{u}_2 = r\dot{\theta} - \vartheta_2, \quad \text{and} \quad \vartheta_2 = k_{21}(q - q_d) + k_{22}(\dot{q} - \dot{q}_d) + u_{2d}.
$$
Let SSPMI($A_{cl}$, $B_{cl}$, $C_{cl}$, $D_{cl}$, $\epsilon$) have a solution $P = P^\top > 0$. Then the origin of (20) is GUUB in a ball with radius $\rho \sqrt{\frac{\lambda_{\max}(P)}{\lambda_{\min}(P)}}$, where $\rho \stackrel{\Delta}{=} \frac{g \| [C_{cl}]_1 \cdot \|}{\epsilon \lambda_{\min}(P)} \|\Delta\mu\|_{L^\infty}$, and all signals are bounded.

**PROOF:** The proof follows the same steps as in the previous section. First, notice that the dynamics of $\tilde{u}_2$ in (13) remains unchanged for the case of unknown $\mu$. Consequently, the finite-time convergence of $\tilde{u}_2$ to zero remains unchanged. For the second step, i.e., the analysis of the error dynamics on the interval $[t^*, \infty)$, let us consider the Lyapunov function candidate $V(e) = \frac{1}{2}e^\top Pe$, where $P = P^\top > 0$ is a solution of the SSPMI($A_{cl}$, $B_{cl}$, $C_{cl}$, $D_{cl}$, $\epsilon$). The time derivative of $V$ along trajectories of (20) is given by

$$
\begin{equation}
\begin{split}
\dot{V}(e) ={}& \frac{1}{2} e^\top (A_{cl}^\top P + P A_{cl}) e - \mu(t) g e^\top P E_{cl} (\lambda_t - \lambda_{td}) \\
& - \Delta \mu(t) g e^\top P E_{cl} \lambda_{td}.
\end{split}
\tag{22}
\end{equation}
$$
Using the fact that $PE_{cl} = -[C_{cl}^{\top}]_{\bullet,1} = \begin{pmatrix} -k_{21} \\ 1 - k_{22} \end{pmatrix}$, and the maximal monotonicity of the
sign multifunction, it follows that

$$
\begin{aligned}
\dot{V}(e) &\leq -\epsilon e^{\top} P e - \Delta\mu(t) g e^{\top} P E_{cl} \lambda_{td} \\
&\leq - (\epsilon \lambda_{\min}(P) \|e\| - g \|\Delta\mu\|_{L^{\infty}} \|[C_{cl}]_{1,\bullet}\|) \|e\|.
\end{aligned}
$$

Hence, whenever $\|e\| > \rho$, then $\dot{V}(e) < 0$. It follows from classical arguments, see, e.g., [40], that the error converges towards the smallest level set of $V$ containing the ball $\rho B$. Hence, the error is guaranteed to converge to a ball of radius $\rho \sqrt{\frac{\lambda_{\max}(P)}{\lambda_{\min}(P)}}$. Boundedness of all signals holds using the same arguments as in Proposition 2. $\square$

# 5 The Backward-Euler Discrete-time Implementation of the Backstepping Controller

It is of interest to analyse the discrete-time implementation of the controller (11) (12), which is set-valued and thus prone to numerical chattering when badly implemented. Backward Euler discretizations drastically decrease the digital chattering [41, 42, 43, 44, 45, 46, 47, 48, 49].

## 5.1 Plant and Controller Discretization

To that aim let us start with a backward (implicit) Euler discretization of the dynamics in (1):

$$ \dot{q}_{k+1} = \dot{q}_k + \frac{h}{m} u_{1,k} - h \mu g \lambda_{t,k+1} \quad (23a) $$

$$ q_{k+1} = q_k + h \dot{q}_k \quad (23b) $$

$$ I \dot{\theta}_{k+1} = I \dot{\theta}_k + h \tau_k + h r \mu m g \lambda_{t,k+1} \quad (23c) $$

$$ \theta_{k+1} = \theta_k + h \dot{\theta}_k \quad (23d) $$

$$ \lambda_{t,k+1} \in \text{sgn}(\dot{q}_{k+1} - r\dot{\theta}_{k+1}) \quad (23e) $$

with $h > 0$ the timestep. It is assumed in this section that $\mu$ is known. The scheme in (23) is named implicit because it is implicit in the set-valued frictional terms (23e). This allows us to get numerical results avoiding chattering at sticking modes [50]. In other words, backward Euler schemes correctly approximate the contact force, while forward Euler schemes do not. It is therefore postulated that (23) is a good time-discretization of the plant's dynamics for the analysis (see Appendix B). The discretization of (3) is given by (assuming $f = 0$ and $k = 0$):

$$
\left\{
\begin{array}{l}
\dot{q}_{d,k+1} = \dot{q}_{d,k} + \frac{h}{m} u_{1,d,k} - h \mu g \lambda_{t,d,k+1} \\
q_{d,k+1} = q_{d,k} + h \dot{q}_{d,k} \\
\lambda_{t,d,k+1} \in \text{sgn}(\dot{q}_{d,k+1} - u_{2,d,k+1})
\end{array}
\right.
\quad (24)
$$
Let us define $\tilde{u}_{2,k} = r\dot{\theta}_k - \vartheta_{2,k}$, where $\vartheta_{2,k} = k_{21}(q_k - q_{d,k}) + k_{22}(\dot{q}_k - \dot{q}_{d,k}) + u_{2d,k}$. After some manipulations we obtain:

$$ \tilde{u}_{2,k+1} = \tilde{u}_{2,k} + h \frac{r}{I} \tau_k - h \Lambda_d(\dot{q}_k, \dot{q}_{d,k}, k) + h \mu g \left( \frac{r^2 m}{I} + k_{22} \right) \lambda_{t,k+1} \quad (25) $$

where

$$ \Lambda_d(\dot{q}_k, \dot{q}_{d,k}, k) = k_{21}(\dot{q}_k - \dot{q}_{d,k}) - k_{22} \left( \frac{\dot{q}_{d,k+1} - \dot{q}_{d,k}}{h} \right) + \frac{k_{22}}{m} u_{1,k} + \frac{u_{2d,k+1} - u_{2d,k}}{h} \quad (26) $$

Mimicking the continuous-time case, let us consider the control input $\tau_k$ as

$$ \frac{r}{I} \tau_k = \Lambda_d(\dot{q}_k, \dot{q}_{d,k}, k) - k_{sl} \lambda_{3,k+1}, \quad (27a) $$

$$ \lambda_{3,k+1} \in \operatorname{sgn}(\tilde{u}_{2,k+1}). \quad (27b) $$

This controller is an implicit (or backward) Euler discrete-time input, due to (27b). Indeed the substitution of (27) into (25) yields the GE:

$$ \tilde{u}_{2,k+1} = \tilde{u}_{2,k} - h k_{sl} \lambda_{3,k+1} + h \mu g \left( \frac{r^2 m}{I} + k_{22} \right) \lambda_{t,k+1} \quad (28a) $$

$$ \lambda_{3,k+1} \in \operatorname{sgn}(\tilde{u}_{2,k+1}) \quad (28b) $$

The next step is to prove that the controller in (27) can be computed at $t = t_k$, in a unique way (i.e., it is well-posed).

## 5.2 Well-posedness of (23), (27), (28)

The system in (23), (27), (28) is a difference inclusion, which is a GE with unknowns $\dot{q}_{k+1}$, $\dot{\theta}_{k+1}$, $\lambda_{t,k+1}$, $\tilde{u}_{2,k+1}$, $\lambda_{3,k+1}$. The inclusion in (23e) is equivalently written as: $\dot{q}_{k+1} - r\dot{\theta}_{k+1} \in \mathcal{N}_{[-1,1]}(\lambda_{t,k+1})$. Using (23a), (23c), and (27), then:

$$ \dot{q}_k + h \frac{u_{1,k}}{m} - r \dot{\theta}_k - h \Lambda_d(q_k, \dot{q}_k, t_k) + h k_{sl} \lambda_{3,k+1} - h \mu g \left(1 + \frac{r^2 m}{I}\right) \lambda_{t,k+1} \in \mathcal{N}_{[-1,1]}(\lambda_{t,k+1}) \quad (29) $$

with:

$$ u_{1,k} = k_{11}(q_k - q_{d,k}) + k_{12}(\dot{q}_k - \dot{q}_{d,k}) + u_{1d,k}. \quad (30) $$
In a similar way, it follows from (28b) that $\tilde{u}_{2,k+1} \in \mathcal{N}_{[-1,1]}(\lambda_{3,k+1})$. Thus, the substitution of (28a), leads us to

$$r\dot{\theta}_k - \vartheta_{2,k} - h k_{sl} \lambda_{3,k+1} + h \mu g \left( k_{22} + \frac{r^2 m}{I} \right) \lambda_{t,k+1} \in \mathcal{N}_{[-1,1]}(\lambda_{3,k+1}) \quad (31)$$

Combining (29) and (31), the following GE with unknown $(\lambda_{t,k+1}, \lambda_{3,k+1})$ is obtained:

$$H(q_k, \dot{q}_k, \dot{\theta}_k, t_k) - M_h \begin{pmatrix} \lambda_{t,k+1} \\ \lambda_{3,k+1} \end{pmatrix} \in \mathcal{N}_{[-1,1]^2} \begin{pmatrix} \lambda_{t,k+1} \\ \lambda_{3,k+1} \end{pmatrix} \quad (32)$$

with $H = (H_1, H_2)^\top$, $H_1 = \dot{q}_k + h \frac{u_{1,k}}{m} - r \dot{\theta}_k - h \Lambda_d(\dot{q}_k, \dot{q}_{d,k}, k)$, $H_2 = r \dot{\theta}_k - k_{21}(q_k - q_{d,k}) - k_{22}(\dot{q}_k - \dot{q}_{d,k} - u_{2d,k})$, and $M_h = h \begin{pmatrix} \mu g \left(1 + \frac{r^2 m}{I}\right) & -k_{sl} \\ -\mu g \left(k_{22} + \frac{r^2 m}{I}\right) & k_{sl} \end{pmatrix}$.

**Lemma 1** Assume that $k_{22} < 1$, $k_{sl} > 0$, and $\mu > 0$. Then the GE in (32) has a unique solution $(\lambda_{t,k+1}, \lambda_{3,k+1})$ for any data $q_k, \dot{q}_k, \dot{\theta}_k, t_k, k \ge 0$.

**PROOF:** The GE in (32) can be written as an affine variational inequality AVI(K, H, $M_h$) [51] where $K = [-1, 1]^2$ is a bounded polyhedral set. It has a solution according to [51, Corollary 2.2.5]. In order to show the uniqueness of the solution to (32), it is sufficient to prove that the matrix $M_h$ is P-matrix [52, Theorem 3]. It is sufficient to check that the second principal minor $k_{sl}\mu g(1 - k_{22}) > 0$, which holds for all $k_{22} < 1$. Thus, the matrix $M_h$ is a P-matrix if $k_{22} < 1$. $\square$

Therefore the controller $\tau_k$ in (27a) to be applied on $[t_k, t_{k+1})$, can be computed with measurements of $q, \dot{q}$ and $\dot{\theta}$ at $t = t_k$.

Existence and uniqueness of the multipliers as shown in Lemma 1, prove that the following is true.

**Corollary 1** The difference inclusion in (23) (27) (28) is well-posed with unique solution $q_{k+1}, \dot{q}_{k+1}, \theta_{k+1}, \dot{\theta}_{k+1}$.

## 5.3 Controller Calculation

There is no unique way to compute the controller by solving the GE. Let us indicate one of them. First it is possible to rewrite the GE (32) as a linear complementarity problem (LCP). Let $\bar{\lambda}_{k+1} \triangleq (\lambda_{t,k+1} \ \lambda_{3,k+1})^\top$ and let us define $\mathcal{C}$ as the set of $\bar{\lambda}_{k+1}$ such that $\bar{\lambda}_{k+1} \in [-1, 1]^2$, then:

$$\mathcal{C} = \{ \bar{\lambda}_{k+1} \in \mathbb{R}^2 \mid R\bar{\lambda}_{k+1} + r \ge 0 \}. \quad (33)$$
$R = \begin{pmatrix} 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & -1 \end{pmatrix}^\top$ and $r = (1 1 1 1)^\top$. Let us define the normal cone to $C$ at $\bar{\lambda}_{k+1}$ as follows:

$$ \begin{aligned} \mathcal{N}_C(\bar{\lambda}_{k+1}) &= \{z \in \mathbb{R}^2 \mid z_1 = -\gamma_1 + \gamma_2 \text{ and } z_2 = -\gamma_3 + \gamma_4, \gamma_i \ge 0 \text{ for } i \in \{1, 2, 3, 4\}\} \\ &= \{z \in \mathbb{R}^2 \mid z = -R^\top \gamma, 0 \le \gamma \perp R \bar{\lambda}_{k+1} + r \ge 0\} \end{aligned} \quad (34) $$

where $\gamma = (\gamma_1, \gamma_2, \gamma_3, \gamma_4)^\top$ is a vector of slack variables. Then, the GE in (32) is equivalently written as the mixed LCP (MLCP):

$$ \left\{ \begin{array}{l} M_h \bar{\lambda}_{k+1} - H(q_k, \dot{q}_k, t_k) = R^\top \gamma \\ 0 \le \gamma \perp R \bar{\lambda}_{k+1} + r \ge 0. \end{array} \right. \quad (35) $$

Given that $M_h$ is P-matrix, the following LCP is derived:

$$ 0 \le \gamma \perp R M_h^{-1} R^\top \gamma + R M_h^{-1} H(q_k, \dot{q}_k, t_k) + r \ge 0. \quad (36) $$

By solving (36), the variables $\lambda_{t,k+1}$ and $\lambda_{3,k+1}$ are obtained. The matrix $R M_h^{-1} R^\top$ is not a P-matrix, since it has low-rank 2 provided that $k_{22} < 1$. However it is noteworthy that the 4 constraints defining $C$ in (33) cannot be activated simultaneously. Therefore an enumeration procedure can be used to compute the solution, which exists with uniqueness. The controller $\tau_k$ in (27a) to be applied on $[t_k, t_{k+1})$, can be computed with measurements of $q, \dot{q}$ and $\dot{\theta}$ at $t = t_k$.

Another way to compute the multipliers $\lambda_{t,k+1}$ and $\lambda_{3,k+1}$ in (32), consists in using proximal splitting methods from convex optimization [53] (another way, not presented here, consists of solving a linear complementarity problem). The following proposition, states the conditions for convergence of the Douglas-Rachford splitting, when applied to the GE (32).

**Proposition 4** Assume that $k_{22} < 1$, $k_{sl} > 0$, and $\mu > 0$. Let $z_0 \in \mathbb{R}^2$ be arbitrary, and consider the following iteration in $s \in \mathbb{N}$:

$$ \hat{\lambda}_s = \operatorname{Proj}_{[-1,1]^2}(z_s), \quad (37a) $$

$$ z_{s+1} = z_s - \hat{\lambda}_s + (I + \delta M_h)^{-1} (2\hat{\lambda}_s - z_s + \delta H_k), \quad (37b) $$

where $M_h$ and $H_k = H(q_k, \dot{q}_k, \dot{\theta}_k, t_k)$ are as in (32). If

$$ \beta - \sqrt{\beta^2 - (\mu\alpha_2)^2} < k_{sl} < \beta + \sqrt{\beta^2 - (\mu\alpha_2)^2}, \quad (38) $$

where $\beta = \mu(2\alpha_1 - \alpha_2)$. Then, for any $\delta > 0$, $\lim_{s \to \infty} \hat{\lambda}_s = [\lambda_{t,k+1}, \lambda_{3,k+1}]^\top$, solution of (32).
**PROOF:** The proof relies on verifying that the map $A(\hat{\lambda}) = M_h\hat{\lambda} - H_k$ is strongly monotone and $B(\hat{\lambda}) = N_{[-1,1]^2}(\hat{\lambda})$ is maximal monotone and then applying [54, Theorem 25.6]. Since $B = \partial\Psi_{[-1,1]^2}$, it is maximal monotone [54, Theorem 20.40]. Regarding $A$, it is maximal strongly monotone if and only if

$$ \frac{1}{2}(M_h + M_h^\top) = h \left[ -\frac{\mu\alpha_1}{k_{sl} + \mu\alpha_2} - \frac{k_{sl} + \mu\alpha_2}{k_{sl}^2} \right] > 0. $$

This last condition holds, if and only if, $\mu\alpha_1 > 0$, $k_{sl} > 0$ and (38) holds. Therefore, under the assumptions of Lemma 1, the GE: $0 \in A(\hat{\lambda}) + B(\hat{\lambda})$ has a unique solution and [54, Theorem 25.6] guarantees the convergence of $\hat{\lambda}_s$ towards it. $\square$

Under the assumptions of Proposition 4 the multiplier $\lambda_{3,k+1}$ appearing in (27) can be approximated, at each time step, by the iteration (37). As demonstrated in section 6, such an approximation is sufficiently accurate not to destroy the closed-loop properties.

## 5.4 Study of (28)

The input in (27) is the counterpart of (11)-(12), while (28) is the counterpart of (13). Thus, similar to the continuous time case, $\tilde{u}_{2,k}$ is finite-time stable.

**Lemma 2** Assume that $k_{22} < 1$, $k_{sl} > 0$, and $\mu > 0$. The GE (28) has always a unique solution $\tilde{u}_{2,k+1}$. Moreover for $k_{sl} > \mu g |\frac{r^2m}{I} + k_{22}|$, the sequence $\{\tilde{u}_{2,k}\}_{k \ge 0}$ converges to zero in a finite number of steps.

The proof is identical to the proofs that can be found in [41, 42] and it is omitted. If a forward Euler discretization is chosen, then the sequence $\{\tilde{u}_{2,k}\}_{k \ge 0}$ does not converge to zero. Consequently the multiplier $\lambda_{3,k}$ in (27b) keeps on oscillating between 1 and -1, and chattering is present in the controller in (27a), see section 6.

## 5.5 Closed-loop Extended Proximal-Point Algorithm

Resolvants of maximal monotone operators and proximal-point algorithms are well-known in Optimization and maximal monotone operators theory [54]. They also appear in backward discretization of sliding-mode systems [44, 55]. Interestingly the following holds.

**Proposition 5** Consider (23), (24), (27), (28), (30). Then:

$$ \tilde{u}_{2,k+1} = J_{\text{sgn}}^{hk_{sl}}(\tilde{u}_{2,k} + h\alpha\lambda_{t,k+1}(q_k, \dot{q}_k, \dot{\theta}_k)), \quad (39) $$

where $\alpha \triangleq \mu g \left(\frac{r^2m}{I} + k_{22}\right)$, and

$$ \begin{aligned} \dot{q}_{k+1} = & \frac{1}{1-k_{22}}(J_{\text{sgn}}^{(1-k_{22})h\mu g}((1-k_{22})(1+\frac{h}{m}k_{12})\dot{q}_k \\ & +\zeta_k(q_k, \dot{q}_k, \dot{\theta}_k, t_k) + \beta_k(q_k, \dot{q}_k, \dot{\theta}_k, t_k)) \\ & -\beta_k(q_k, \dot{q}_k, \dot{\theta}_k, t_k)) \end{aligned} \quad (40) $$
where $\zeta_k \stackrel{\Delta}{=} \frac{h}{m}k_{11}(q_k - q_{d,k}) + \dot{q}_{d,k} - \dot{q}_{d,k+1} + \frac{h}{m}u_{1d,k}$, and $\beta_k \stackrel{\Delta}{=} -k_{21}(q_k + h\dot{q}_k - q_{d,k+1}) - u_{2d,k+1}$, $1 - k_{22} > 0$, and

$$r\dot{\theta}_{k+1} - \vartheta_{2,k+1} = J_{\text{sgn}}^{hk_{sl}}(r\dot{\theta}_k - \vartheta_{2,k} + \gamma_k(q_k, \dot{q}_k, \dot{\theta}_k, t_k)), \quad (41)$$

$$\text{with } \gamma_k \stackrel{\Delta}{=} \vartheta_{2,k} - \vartheta_{2,k+1} + h\Lambda(q_k, \dot{q}_k, t_k) + h \frac{r^2 \mu g}{I} \lambda_{t,k+1}.$$

## PROOF:

**Proof of (39):** From Lemma 1, $\lambda_{t,k+1} = \lambda_{t,k+1}(q_k, \dot{q}_k, \dot{\theta}_k)$. Then (39) is a direct consequence of (28) and the definition of the resolvant.

**Proof of (40):** Using (23a) and (24) yields:

$$\dot{q}_{k+1} - \dot{q}_{d,k+1} = \left(1 + \frac{h}{m}k_{11}\right)(\dot{q}_k - \dot{q}_{d,k}) + \frac{h}{m}k_{12}(q_k - q_{d,k}) + \dot{q}_{d,k} - \dot{q}_{d,k+1} + \frac{h}{m}u_{1d,k} - h\mu g \lambda_{t,k+1}. \quad (42)$$

Assume that $\tilde{u}_{2,k} = 0$ for all $k > k^*$ (from Lemma 2 such a $k^* < +\infty$ exists), then $r\dot{\theta}_k = \vartheta_{2,k}$ and:

$$\dot{q}_{k+1} - \dot{q}_{d,k+1} \in \left(1 + \frac{h}{m}k_{12}\right)(\dot{q}_k - \dot{q}_{d,k}) + \frac{h}{m}k_{11}(q_k - q_{d,k}) + \dot{q}_{d,k} - \dot{q}_{d,k+1} + \frac{h}{m}u_{1d,k} - h\mu g \operatorname{sgn}(\dot{q}_{k+1} - k_{21}(q_{k+1} - q_{d,k+1}) - k_{22}(\dot{q}_{k+1} - \dot{q}_{d,k+1}) - u_{2d,k+1}) \quad (43)$$

Then (43) is rewritten equivalently as:

$$\dot{q}_{k+1} - \dot{q}_{d,k+1} \in \left(1 + \frac{h}{m}k_{12}\right)(\dot{q}_k - \dot{q}_{d,k}) + \zeta_k - h\mu g \operatorname{sgn}\left((1 - k_{22})(\dot{q}_{k+1} - \dot{q}_{d,k+1}) + \beta_k\right). \quad (44)$$

Let us denote $\dot{\tilde{q}}_k \stackrel{\Delta}{=} \dot{q}_k - \dot{q}_{d,k}$. Then (44) is rewritten equivalently as:

$$\begin{align}
\dot{\tilde{q}}_{k+1} &\in \left(1 + \frac{h}{m}k_{12}\right) \dot{\tilde{q}}_k + \zeta_k - h\mu g \operatorname{sgn}\left((1 - k_{22})\dot{\tilde{q}}_{k+1} + \beta_k\right) \nonumber \\
&\Downarrow \nonumber \\
&\quad (1 - k_{22})\dot{\tilde{q}}_{k+1} + \beta_k \in \beta_k + (1 - k_{22}) \left(1 + \frac{h}{m}k_{12}\right) \dot{\tilde{q}}_k \nonumber \\
&\quad + \zeta_k - (1 - k_{22})h\mu g \operatorname{sgn}\left((1 - k_{22})\dot{\tilde{q}}_{k+1} + \beta_k\right) \tag{45} \\
&\Downarrow \nonumber \\
&\quad (1 - k_{22})\dot{\tilde{q}}_{k+1} + \beta_k = J_{\operatorname{sgn}}^{(1-k_{22})h\mu g} \left((1 - k_{22}) \left(1 + \frac{h}{m}k_{12}\right) \dot{\tilde{q}}_k + \zeta_k + \beta_k\right), \nonumber
\end{align}$$

where it is used that $1 - k_{22} > 0$ as in Lemmata 1 and 2. Equivalence with (40) follows. By definition $r\dot{\theta}_k = \tilde{u}_{2,k} + \vartheta_{2,k}$. Therefore the above reasoning can be redone, introducing $\tilde{u}_{2,k}$ in (43), in $\operatorname{sgn}(\dot{q}_{k+1} - k_{21}(q_{k+1} - x_{1d,k+1}) - k_{22}(\dot{q}_{k+1} - x_{2d,k+1}) - u_{2d,k+1} + \tilde{u}_{2,k+1})$ and adjusting $\beta_k$ accordingly.

**Proof of (41):** We have $\vartheta_{2,k+1} = \vartheta_{2,k+1}(q_{k+1}, \dot{q}_{k+1}, t_k) = \vartheta_{2,k+1}(q_k, \dot{q}_k, \dot{\theta}_k, t_k)$. Hence $\gamma_k =$
$\gamma_k(q_k, \dot{q}_k, \dot{\theta}_k, t_k)$. Consider now (23c). Inserting (27) into (23c) yields:

$$
\begin{align}
& r\dot{\theta}_{k+1} \in r\dot{\theta}_k + h\Lambda(q_k, \dot{q}_k, t_k) \nonumber \\
& \qquad -hk_{sl} \operatorname{sgn}(-\vartheta_{2,k+1} + r\dot{\theta}_{k+1}) + h\frac{r^2\mu mg}{I}\lambda_{t,k+1} \nonumber \\
& \Uparrow \\
& r\dot{\theta}_{k+1} - \vartheta_{2,k+1} \in r\dot{\theta}_k - \vartheta_{2,k} + \vartheta_{2,k} + h\Lambda(q_k, \dot{q}_k, t_k) \nonumber \\
& \qquad -\vartheta_{2,k+1} - hk_{sl} \operatorname{sgn}(-\vartheta_{2,k+1} + r\dot{\theta}_{k+1}) + h\frac{r^2\mu mg}{I}\lambda_{t,k+1} \nonumber \\
& \Uparrow \\
& r\dot{\theta}_{k+1} - \vartheta_{2,k+1} \in r\dot{\theta}_k - \vartheta_{2,k} + \gamma_k - hk_{sl} \operatorname{sgn}(r\dot{\theta}_{k+1} - \vartheta_{2,k+1}). \tag{46}
\end{align}
$$

This shows that (41) holds. $\square$

It is noteworthy that (39), (40) and (41) are is similar to a proximal point algorithm [54]. However the presence of "perturbations" introduces significant differences. They may be named *extended proximal-point algorithm*. The behaviour of (39) is analysed in Lemma 2. The behaviour of (40) and (41) is a consequence of the analysis in section 5.6. *This allows us to extend the results on classical proximal-point algorithms.*

## 5.6 Stability Analysis

The stability analysis of (43) is inspired from its continuous-time counterpart. Using (23)-(24) the discrete-time error dynamics are given by:

$$
\left\{
\begin{aligned}
\dot{\tilde{q}}_{k+1} &= \dot{\tilde{q}}_k + h\mu g(\lambda_{t,d,k+1} - \lambda_{t,k+1}) \\
&\quad + \frac{h}{m}(u_{1,k} - u_{1,d,k}) && (47a) \\
u_{1,k} &= k_{11}\tilde{q}_k + k_{12}\dot{\tilde{q}}_k + u_{1,d,k} && (47b) \\
\lambda_{t,d,k+1} &\in \operatorname{sgn}(\dot{q}_{d,k+1} - u_{2,d,k}) && (47c) \\
\lambda_{t,k+1} &\in \operatorname{sgn}(\dot{q}_{k+1} - r\dot{\theta}_{k+1}) = \operatorname{sgn}(\dot{q}_{k+1} - \tilde{u}_{2,k+1} \\
&\quad - \vartheta_{2,k+1}) && (47d) \\
\vartheta_{2,k} &= k_{21}\tilde{q}_k + k_{22}\dot{\tilde{q}}_k + u_{2,d,k} && (47e) \\
\tilde{q}_{k+1} &= \tilde{q}_k + h\dot{\tilde{q}}_k. && (47f)
\end{aligned}
\right.
$$

with $\tilde{q}_k = q_k - q_{d,k}$, $\dot{\tilde{q}}_k = \dot{q}_k - \dot{q}_{d,k}$. Thus, similarly to the continuous-time design, we now focus on the stability of the zero solution of the discrete error dynamics. Setting $e_k = [\tilde{q}_k, \dot{\tilde{q}}_k]^\top$, it follows from (30) (47) that the discrete time error dynamics is given by:

$$
\left\{
\begin{aligned}
e_{k+1} &= (I + hA_{cl})e_k - h\mu g E_{cl}(\lambda_{t,k+1} - \lambda_{td,k+1}) && (48a) \\
\lambda_{t,k+1} &\in \operatorname{sgn}(\dot{q}_{k+1} - r\dot{\theta}_{k+1}) && (48b) \\
\lambda_{t,d,k+1} &\in \operatorname{sgn}(\dot{q}_{d,k+1} - u_{2d,k+1}) && (48c)
\end{aligned}
\right.
$$
with the matrices $A_{cl}$ and $E_{cl}$ as in (5) and $\dot{\theta}_{k+1}$ satisfies (23c). Recalling that $r\dot{\theta}_{k+1} = \tilde{u}_{2,k+1} + \vartheta_{2,k+1}$, (with $\vartheta_{2,k}$ defined after (24)) and that $\tilde{u}_{2,k+1}$ becomes zero after a finite number of steps $k^*$ (Lemma 2), it is sufficient to study the difference inclusion (48) for $k \ge k^*$, so that the error dynamics becomes equal to

$$e_{k+1} = (I + hA_{cl})e_k - h\mu g E_{cl}(\lambda_{t,k+1} - \lambda_{t,d,k+1}) \quad (49a)$$

$$\lambda_{t,k+1} \in \mathbf{sgn}(C_d e_{k+1} + \dot{q}_{d,k+1} - u_{2d,k+1}) \quad (49b)$$

$$\lambda_{t,d,k+1} \in \mathbf{sgn}(\dot{q}_{d,k+1} - u_{2d,k+1}) \quad (49c)$$

where $C_d = -[C_{cl}]_{1,\bullet}$ with $C_{cl}$ as defined in Section 2. The following proposition states the asymptotic stability of the discrete error dynamics.

**Proposition 6** Let all the assumptions of Lemma 2 hold. Consider the difference inclusion in (49), with bounded initial data. Assume that $(A_{cl}, B_{cl}, C_{cl}, D_{cl})$ is strictly state passive, and that $P = P^\top > 0$ is a solution to the SSPLMI$(A_{cl}, B_{cl}, C_{cl}, D_{cl}, \varepsilon)$. Thus, for $h$ sufficiently small, such that

$$\varepsilon P - hA_{cl}^\top PA_{cl} > 0, \quad (50)$$

the origin of the error dynamics (49) is globally asymptotically stable.

**PROOF:** It is easy to see that any $P = P^\top > 0$ solution to the SSPLMI$(A_{cl}, B_{cl}, C_{cl}, D_{cl}, \varepsilon)$ is also a solution to the SSPLMI$(A_{cl}, E_{cl}, C_d, 0, \varepsilon)$, with $C_d = (-k_{21}, 1 - k_{22})^\top$. In particular: i) $A_{cl}^\top P + PA_{cl} \le -\varepsilon P$ and ii) $PE_{cl} = C_d^\top$. Keeping this two facts in mind, consider the candidate Lyapunov function $V(e_k) = \frac{1}{2}e_k^\top Pe_k$, and let $\Delta V = V(e_{k+1}) - V(e_k)$. It follows from (49) that

$$\begin{aligned} \Delta V &= \frac{1}{2}e_{k+1}^\top Pe_{k+1} - \frac{1}{2}e_k^\top Pe_k \\ &= \frac{1}{2}e_k^\top (I + hA_{cl})P(I + hA_{cl})e_k - \frac{1}{2}e_k^\top Pe_k \\ &\quad - h\mu g \Delta\lambda_{t,k+1}^\top E_{cl}^\top P(I + hA_{cl})e_k + \frac{(h\mu g)^2}{2} \Delta\lambda_{t,k+1} E_{cl}^\top P E_{cl} \Delta\lambda_{t,k+1} \\ &= \frac{1}{2}e_k^\top (I + hA_{cl})P(I + hA_{cl})e_k - \frac{1}{2}e_k^\top Pe_k - h\mu g \Delta\lambda_{t,k+1}^\top E_{cl}^\top P e_{k+1} \\ &\quad - \frac{(h\mu g)^2}{2} \Delta\lambda_{t,k+1} E_{cl}^\top P E_{cl} \Delta\lambda_{t,k+1}, \end{aligned} \quad (51)$$

where $\Delta\lambda_{t,k+1} = \lambda_{t,k+1} - \lambda_{t,d,k+1}$. Using $E_{cl}^\top P = C_d$ and the monotonicity of the sign multifunction it follows from (49b) and (49c) that $\Delta\lambda_{t,k+1}^\top E_{cl}^\top P e_{k+1} \ge 0$ and

$$\begin{aligned} \Delta V &\le \frac{1}{2}e_k^\top (I + hA_{cl})P(I + hA_{cl})e_k - \frac{1}{2}e_k^\top Pe_k \\ &= \frac{h}{2}e_k^\top (A_{cl}^\top P + PA_{cl} + hA_{cl}^\top PA_{cl})e_k \\ &\le -\frac{h}{2}e_k^\top (\varepsilon P - hA_{cl}^\top PA_{cl})e_k, \end{aligned} \quad (52)$$
and the conclusion follows in view of (50). □

**Remark 4** The variable $\dot{\theta}_k$ is uniformly bounded as well under the conditions in Lemma 2 and in Proposition 6. Indeed we have $r\dot{\theta}_k = \tilde{u}_{2,k} + \vartheta_{2,k} = \tilde{u}_{2,k} + k_{21}(q_k - q_{d,k}) + k_{22}(\dot{q}_k - \dot{q}_{d,k}) + u_{2d,k}$. Lemma 2 and Proposition 6 allow us to conclude. After a finite number of steps it is true that $r\dot{\theta}_k = \tilde{u}_{2,k} + \vartheta_{2,k} = k_{21}(q_k - q_{d,k}) + k_{22}(\dot{q}_k - \dot{q}_{d,k}) + u_{2d,k}$.

## 5.7 Case with Unknown $\mu$

Similarly to the continuous-time case, in this section we allow for $\mu$ to be time-dependent, i.e., $\mu_k = \mu(q_k, \dot{q}_k, t_k)$. The explicit expression for $\mu_k$ is assumed unknown, however, an upper bound is assumed available, that is $|\mu_k| \le \bar{\mu}$ for all $k \in \mathbb{N}$ and some known $0 \le \bar{\mu} < \infty$. In this case, it is not possible to consider (27b) with this selection $\lambda_{3,k+1}$. Indeed using (23a) and (23c), it follows that $\tilde{u}_{2,k+1} = r\dot{\theta}_{k+1} - \vartheta_{2,k+1}$ depends on $\mu$ which is unknown. Another way to see this is to consider that $M_h$ in (32) depends on $\mu$, hence $\lambda_{3,k+1}$ in (27) used to compute the controller, depends on $\mu$. Therefore the reasoning in section 5.2 cannot be applied directly to the current case.

### 5.7.1 Controller Calculation

To overcome this issue, the control input $\tau_k$ is modified in the following way

$$ \frac{r}{I} \tau_k = \Lambda_d(\dot{q}_k, \dot{q}_{d,k}, k) - k_{sl} \lambda_{3,k+1} \quad (53a) $$

$$ \lambda_{3,k+1} \in \operatorname{sgn}\left(H_{2,k} - h k_{sl} \lambda_{3,k+1} + h \mu_{0,k} \alpha_2 \tilde{\lambda}_{t,k+1}\right) \quad (53b) $$

$$ \tilde{\lambda}_{t,k+1} \in \operatorname{sgn}\left(H_{1,k} + h k_{sl} \lambda_{3,k+1} - h \mu_{0,k} \alpha_1 \tilde{\lambda}_{t,k+1}\right) \quad (53c) $$

where $\mu_{0,k} \ge 0$ denotes the nominal (known) value of the friction coefficient and, like above, $\alpha_1 = g\left(\frac{r^2m}{I} + 1\right)$ and $\alpha_2 = g\left(\frac{r^2m}{I} + k_{22}\right)$. The arguments of the signum multifunction in (53b) and (53c) replace the unknown $\lambda_{t,k+1}$ with its nominal value denoted $\tilde{\lambda}_{t,k+1}$. Note that (53) is independent of $\mu_k$. Moreover, it is well-defined for all $k \in \mathbb{N}$. Indeed, the subsystem (53b)-(53c) can be rewritten equivalently as:

$$ H(q_k, \dot{q}_k, \dot{\theta}_k, t_k) - h \begin{pmatrix} \mu_{0,k} \alpha_1 & -k_{sl} \\ -\mu_{0,k} \alpha_2 & k_{sl} \end{pmatrix} \begin{pmatrix} \tilde{\lambda}_{t,k+1} \\ \lambda_{3,k+1} \end{pmatrix} \in \mathcal{N}_{[-1,1]^2} \begin{pmatrix} \tilde{\lambda}_{t,k+1} \\ \lambda_{3,k+1} \end{pmatrix} \quad (54) $$

with same matrices as in (32). Thus the well-posedness of the GE (54) is established, for each $k \in \mathbb{N}$, under the assumptions of Lemma 1. Hence, the control input $\tau_k$ in (53) is uniquely defined for all $k \in \mathbb{N}$. It is worth to remark that if $\mu_{0,k} = 0$, then the solution to the GE (54) is not unique when $H_{1,k} + h k_{sl} \lambda_{3,k+1} = 0$. Nevertheless, the multiplier (selection)
$\lambda_{3,k+1}$, and consequently the input $\tau_k$, is always uniquely defined. Indeed, if $\mu_{0,k} = 0$ then (53b) becomes

$$ \lambda_{3,k+1} \in \operatorname{sgn}(H_{2,k} - h k_{sl} \lambda_{3,k+1}) \qquad (55) $$

whose unique solution is

$$ \lambda_{3,k+1} = \operatorname{Proj}_{[-1,1]} \left( \frac{H_{2,k}}{h k_{sl}} \right), \qquad (56) $$

such that $\tau_k$ is uniquely determined by (53a). Now, the new control input (53) transforms the dynamics of $\tilde{u}_{2,k+1}$ in (25) into

$$ \tilde{u}_{2,k+1} = \varrho_{k+1} + h \alpha_2 (\mu_k \lambda_{t,k+1} - \mu_{0,k} \tilde{\lambda}_{t,k+1}) \qquad (57a) $$

$$ \varrho_{k+1} = H_{2,k} - h k_{sl} \lambda_{3,k+1} + h \mu_{0,k} \alpha_2 \tilde{\lambda}_{t,k+1} \qquad (57b) $$

$$ \lambda_{3,k+1} \in \operatorname{sgn}(\varrho_{k+1}) \qquad (57c) $$

$$ \tilde{\lambda}_{t,k+1} \in \operatorname{sgn}(H_{1,k} + h k_{sl} \lambda_{3,k+1} - h \mu_{0,k} \alpha_1 \tilde{\lambda}_{t,k+1}) \qquad (57d) $$

To see explicitly the difference between $\lambda_{t,k+1}$ and $\tilde{\lambda}_{t,k+1}$, we recall using (23a) (23c) that $\dot{q}_{k+1} - r \dot{\theta}_{k+1} = H_{1,k} + h k_{sl} \lambda_{3,k+1} - h \alpha_1 \mu_k \lambda_{t,k+1}$. Hence, under the assumption that $\mu_k > 0$ for all $k \in \mathbb{N}$, it follows from (23e) and [54, Example 23.4] that

$$ \lambda_{t,k+1} = \operatorname{Proj}_{[-1,1]} \left( \frac{H_{1,k} + h k_{sl} \lambda_{3,k+1}}{h \alpha_1 \mu_k} \right), \qquad (58) $$

whereas, assuming that $\mu_{0,k} > 0$ for all $k \in \mathbb{N}$, (53c) is equivalent to

$$ \tilde{\lambda}_{t,k+1} = \operatorname{Proj}_{[-1,1]} \left( \frac{H_{1,k} + h k_{sl} \lambda_{3,k+1}}{h \alpha_1 \mu_{0,k}} \right), \qquad (59) $$

It is thus clear that when there is no uncertainty in the friction coefficient (i.e., $\mu_k = \mu_{0,k}$ for all $k \in \mathbb{N}$), then $\lambda_{t,k+1} = \tilde{\lambda}_{t,k+1}$, implying that $\tilde{u}_{2,k+1} = \varrho_{k+1}$ and the scheme (57) reduces to (28). Notice also that the existence and uniqueness of $\lambda_{3,k+1}$ implies that of $\lambda_{t,k+1}$ whenever $\mu_k \neq 0$. On the other hand, if $\mu_k = 0$, from (23e) we have $\lambda_{t,k+1} \in \operatorname{sgn}(H_{1,k} + h k_{sl} \lambda_{3,k+1}) \subseteq [-1, 1]$, but in this case, the exact value of $\lambda_{t,k+1}$ becomes irrelevant, since it plays no role in the dynamics (23) when $\mu_k = 0$. We have thus proved the following.

**Corollary 2** For any two sequences of uniformly bounded, non-negative friction coefficients, $\{\mu_k\}_{k \in \mathbb{N}}$ and $\{\mu_{0,k}\}_{k \in \mathbb{N}}$, the controller $\tau_k$ is unique for all $k \in \mathbb{N}$ and the difference inclusion (23), (53), (57) possesses a unique solution $\{q_k\}_{k \in \mathbb{N}}$, $\{\dot{q}_k\}_{k \in \mathbb{N}}$, $\{\theta_k\}_{k \in \mathbb{N}}$, $\{\dot{\theta}_k\}_{k \in \mathbb{N}}$.

## 5.7.2 Stability Analysis

In what follows we focus on the stability properties of the zero solution of the closed-loop error dynamics. We start with the study of the dynamics of $\tilde{u}_{2,k+1}$ in (57).
**Proposition 7** Consider the scheme (57) and let $0 \le \bar{\mu}_0 < \infty$ be such that $\mu_{0,k} \le \bar{\mu}_0$ for all $k \in \mathbb{N}$. If

$$k_{sl} > \left| \frac{r^2 m}{I} + k_{22} \right| \left( \bar{\mu}_0 + \sup_{k \in \mathbb{N}} |\mu_k - \mu_{0,k}| \right), \quad (60)$$

then there is $0 < k^* < \infty$ such that $\varrho_{k+1} = 0$ for all $k \ge k^*$. In particular

$$|\tilde{u}_{2,k+1}| \le h \left| \frac{r^2 m}{I} + k_{22} \right| \sup_{k \in \mathbb{N}} |\mu_k - \mu_{0,k}|, \quad (61)$$

for all $k \ge k^*$.

**PROOF:** It follows from (57a)-(57b) that $\varrho_{k+1}$ satisfies

$$\varrho_{k+1} = \varrho_k - h k_{sl} \lambda_{3,k+1} + h \mu_{0,k} \alpha_2 \tilde{\lambda}_{t,k+1} + h \alpha_2 (\mu_{k-1} \lambda_{t,k} - \mu_{0,k-1} \tilde{\lambda}_{t,k}), \quad (62)$$

since $H_{2,k} = r\dot{\theta}_k - \vartheta_{2,k} = \tilde{u}_{2,k}$. Taking the candidate Lyapunov function $V(\varrho_k) = \frac{1}{2}|\varrho_k|^2$, it follows that

$$\begin{aligned} \Delta V &= \varrho_{k+1}^2 - \frac{1}{2}(\varrho_{k+1}^2 + \varrho_k^2) \\ &= \varrho_{k+1} \left( \varrho_k - h k_{sl} \lambda_{3,k+1} + h \mu_{0,k} \alpha_2 \tilde{\lambda}_{t,k+1} \right. \\ & \qquad \left. - \frac{1}{2} (\varrho_{k+1}^2 + \varrho_k^2) \right) \\ & \le -\frac{1}{2} (\varrho_{k+1} - \varrho_k)^2 - h \left( k_{sl} - |\alpha_2|(\bar{\mu}_0 - |\mu_{k-1} \lambda_{t,k} - \mu_{0,k-1} \tilde{\lambda}_{t,k}|) \right) |\varrho_{k+1}|. \end{aligned} \quad (63)$$

Now, from (58)-(59), simple computations show that

$$\mu_k \lambda_{t,k+1} = \operatorname{Proj}_{[-\mu_k, \mu_k]} (a_k), \quad (64)$$

$$\mu_{0,k} \tilde{\lambda}_{t,k+1} = \operatorname{Proj}_{[-\mu_{0,k}, \mu_{0,k}]} (a_k), \quad (65)$$

where $a_k = \frac{H_{1,k-1} + h k_{sl} \lambda_{3,k}}{h \alpha_1}$. Hence,

$$\begin{aligned} & |\mu_{k-1} \lambda_{t,k} - \mu_{0,k-1} \tilde{\lambda}_{t,k}| \\ &= \left| \operatorname{Proj}_{[-\mu_{k-1}, \mu_{k-1}]} (a_k) - \operatorname{Proj}_{[-\mu_{0,k-1}, \mu_{0,k-1}]} (a_k) \right| \\ & \le |\mu_{k-1} - \mu_{0,k-1}| \le \sup_{k \in \mathbb{N}} |\mu_k - \mu_{0,k}|. \end{aligned} \quad (66)$$

The conclusion follows by noticing that (60) and (66) imply in (63) that $\Delta V \le -\varepsilon|\varrho_{k+1}|$, for some $\varepsilon > 0$ sufficiently small, and the finite-time stability of the zero solution associated with $\varrho_{k+1}$ follows. $\square$
Proposition 7 shows that $\tilde{u}_{2,k+1}$ stays uniformly bounded for all $k \in \mathbb{N}$, with an upper bound proportional to the difference between the friction coefficients. Let us show that the error is ultimately bounded. The associated error dynamics is now given by

$$e_{k+1} = (I + hA_{cl})e_k - hgE_{cl}(\mu_k\lambda_{t,k+1} - \mu_{0,k}\lambda_{t,d,k+1}) \quad (67a)$$

$$\lambda_{t,k+1} \in \text{sgn}(C_d e_{k+1} + \dot{q}_{d,k+1} - u_{2d,k+1} - \tilde{u}_{2,k+1}) \quad (67b)$$

$$\lambda_{t,d,k+1} \in \text{sgn}(\dot{q}_{d,k+1} - u_{2d,k+1}), \quad (67c)$$

with $\tilde{u}_{2,k+1}$ satisfying (57) and $\dot{q}_{d,k+1}$, $\lambda_{t,d,k+1}$ given by (24) with friction coefficient $\mu_{0,k}$.

**Corollary 3** Consider the difference inclusion (57) and (67), with bounded initial data and bounded friction coefficients ($0 \le \mu_k \le \bar{\mu}$ and $0 \le \mu_{0,k} \le \bar{\mu}_0$ for all $k \in \mathbb{N}$). Assume that $P = P^\top > 0$ is a solution to the SSPLMI($A_{cl}, B_{cl}, C_{cl}, D_{cl}, \varepsilon$). If $k_{sl} > 0$ satisfies (60) and $h > 0$ is sufficiently small such that (50) holds, then the origin of (67) is GUUB in a ball with radius $\rho\sqrt{\frac{\lambda_{\max}(P)}{\lambda_{\min}(P)}}$ where

$$\rho = \frac{\sqrt{g \|\Delta\mu\|_{l^{\infty}}}}{\lambda_{\min}(\varepsilon P - hA_{cl}^{\top}PA_{cl})} \Gamma, \quad (68)$$

$$\Gamma = \xi + \sqrt{\xi^2 + 2h\lambda_{\min}(\varepsilon P - hA_{cl}^{\top}PA_{cl})\delta}, \quad (69)$$

$$\xi = \|C_d(I + hA_{cl})\| \sqrt{g \|\Delta\mu\|_{l^{\infty}}}, \quad \text{and} \quad \delta = 2\bar{\mu}\alpha_2 + g(\bar{\mu} + \bar{\mu}_0) \|E_{cl}\|.$$

**PROOF:** Similarly to the proof of Proposition 6, consider the function $V(e_k) = \frac{1}{2}e_k^\top Pe_k$. This yields:

$$\begin{align}
\Delta V &\le -\frac{h}{2}e_k^\top(\varepsilon P - hA_{cl}^\top PA_{cl})e_k - hg(\mu_k\lambda_{t,k+1} - \mu_{0,k}\lambda_{td,k+1})C_d e_{k+1} \nonumber \\
&\quad - \frac{(hg)^2}{2}(\mu_k\lambda_{t,k+1} - \mu_{0,k}\lambda_{td,k+1})^2 E_{cl}^\top P E_{cl} \nonumber \\
&\le -\frac{h}{2}e_k^\top(\varepsilon P - hA_{cl}^\top PA_{cl})e_k + hg\mu_k(\lambda_{t,k+1} - \lambda_{td,k+1})\tilde{u}_{2,k+1} \nonumber \\
&\quad - hg(\mu_k - \mu_{0,k})\lambda_{t,k+1}C_d e_{k+1} \nonumber \\
&\le -\frac{h}{2}e_k^\top(\varepsilon P - hA_{cl}^\top PA_{cl})e_k \nonumber \\
&\quad + hg\bar{\mu}|\lambda_{t,k+1} - \lambda_{td,k+1}|\|\tilde{u}_{2,k+1}\| + hg|\mu_k - \mu_{0,k}|(\|C_d(I + hA_{cl})e_k\| \nonumber \\
&\quad + hg(\bar{\mu} + \bar{\mu}_0)\|E_{cl}\|), \tag{70}
\end{align}$$
where we have used the monotonicity property of the sign multifunction to get the second inequality. Since, by assumption, $k_{sl}$ satisfies (60), it follows from (61) in Proposition 7 that, after a finite number of steps, the difference $\Delta V$ satisfies

$$ \Delta V \leq -\frac{h}{2} e_k^\top (\varepsilon P - h A_{cl}^\top P A_{cl}) e_k \\ + h g \sup_{k \in \mathbb{N}} |\mu_k - \mu_{0,k}| \left( 2h\bar{\mu}\alpha_2 \right. \\ \left. + \|C_d(I + hA_{cl})e_k\| + h g(\bar{\mu} + \bar{\mu}_0) \|E_{cl}\| \right). \quad (71) $$

Simple computations show that the right-hand side of (71) is negative whenever $\|e_k\|^2 > \rho$. The ultimate boundedness of the error $e_k$ follows in view of [56, Corollary 4.11.2]. □

# 6 Numerical Simulations

Consider the tracking problem for the plant (1) and (3), with the following parameters for both systems: $m = 1$[Kg], $k = 0$[N/m], $f = 0$[(N/m)/s], $r = 0.25$[m], $g = 9.81$ [m/s²], $I = \frac{mr^2}{2}$ [Kg·m²]. The nominal friction coefficient of the desired system is set to the constant value $\mu_0 = 1.0$, whereas the friction coefficient of the plant is set to

$$ \mu(t) = \max \left\{ \sin(5t) \cos(t) + \cos(\pi t) \cos(t/\sqrt{7}), 0.1 \right\}. \quad (72) $$

The inputs for the desired dynamics are set as

$$ u_{1,d}(t) = 1.5 \cos(2\pi t) \quad \text{and} \quad u_{2,d}(t) = \sqrt{2} \sin(\sqrt{3}t). \quad (73) $$

Two different discrete-time controllers are considered to address the tracking problem. The first one is the controller (53) which is designed with a backward Euler discretization of (21). As second controller, we consider the forward Euler discretization of (21)

$$ \frac{r}{I} \tau_{\text{exp},k} = \Lambda_d(\dot{q}_k, \dot{q}_{d,k}, t_k) - k_{sl} \operatorname{sgn}(r\dot{\theta}_k - \vartheta_{2,k}), \quad (74) $$

with $\Lambda_d(\dot{q}_k, \dot{q}_{d,k}, t_k)$ as in (26) and $\vartheta_{2,k} = k_{21}(q_k - q_{d,k}) + k_{22}(\dot{q}_k - \dot{q}_{d,k}) + u_{2,k}$. In what follows the subindex 'exp' is used to indicate the variables obtained from closing the loop with (74). The resulting closed-loops are simulated considering a sampling time of $h = 10^{-2}$ seconds. Regarding the parameters of each controller, the gain matrix $K$ is computed from the SSPLMI($A_{cl}, B_{cl}, C_{cl}, D_{cl}, \varepsilon$) with $\varepsilon = 3.0$ and with $P$ satisfying (50) yielding

$$ P = \begin{bmatrix} 45.42 & 8.52 \\ 8.52 & 2.9 \end{bmatrix}, \quad K = \begin{bmatrix} -21.506 & -6.995 \\ -8.521 & -1.902 \end{bmatrix}. \quad (75) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Positions and velocities in (1) and (3) with the control laws (53) and (74).</figcaption></figure>

For the controller (53), the computation of the multipliers $\tilde{\lambda}_{t,k+1}$ and $\lambda_{3,k+1}$ is done using Proposition 4. Thus, with the above parameters, and taking into account the constraint imposed by (60) in Proposition 7, the gain $k_{sl}$ is constrained to the set

$$1.56 < k_{sl} < 115.806. \qquad (76)$$

In the following we set $k_{sl} = 15.0$ for both control laws. Figures 3 through 8 show the time evolution of closed-loop variables. Notably the proposed strategy (53) shows no chattering in the input, leading to small tracking errors, contrary to the forward Euler controller (74), where the input chattering propagates to the state $\dot{q}_{exp}$, affecting also the friction force $g\mu(t)\lambda_{t,exp}(t)$. This chattering reduction is a landmark of backward Euler schemes [41, 42], allowing better precision in the closed-loop. It is noteworthy that the splitting approach of Proposition 4 converges to a solution of (54) with a precision of the order of $10^{-15}$, taking in average 25 iterations at each time step. Pulleys' angular velocity remains bounded (Fig. 7), while the relative velocity shows a stick/slip behaviour (Fig. 8). Due to lack of space the issue of optimal choice of control gains in Proposition 3 and Corollary 3 is not tackled. It can be formulated as an optimization problem under constraints.

# 7 Conclusions

This article is largely concerned with trajectory tracking control of frictional oscillators, a class of set-valued nonsmooth dynamical systems which has received much attention in the Nonlinear Dynamics scientific community. It is shown how backstepping and passivity can be combined to achieve trajectory tracking when set-valued friction is considered, thus
<figure><img src="image_5.png" /><figcaption>Figure 4: (Upper part) Position errors $|q_{\text{exp}}(t) - q_d(t)|$ and $|q(t) - q_d(t)|$, with controllers (74) and (53), respectively. (Lower part) Velocity errors $|\dot{q}_{\text{exp}}(t) - \dot{q}_d(t)|$ and $|\dot{q}(t) - \dot{q}_d(t)|$, with controllers (74) and (53), respectively.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: (Upper part) Control inputs (74) (gray line) and (53) (blue line). (Lower part) Backstepping variables $\tilde{u}_{2,\text{exp}}(t)$ and $\tilde{u}_2(t)$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Friction forces $f_1(t) = g\mu(t)\lambda_{\text{exp},t}(t)$ (gray line), $f_1(t) = g\mu(t)\lambda_{t,d}(t)$ (blue line), and $f_{1,d}(t) = g\mu_0\lambda_{t,d}(t)$ (red line).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Angular velocities $\dot{\theta}_{\text{exp}}$ and $\dot{\theta}$ with inputs (74) and (53), respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Relative velocities $\dot{q}_{\text{rel,exp}}(t) = \dot{q}_{\text{exp}}(t) - r\dot{\theta}_{\text{exp}}(t)$ and $\dot{q}_{\text{rel}}(t) = \dot{q}(t) - r\dot{\theta}(t)$ with controllers (74) and (53), respectively.</figcaption></figure>

extending backstepping to a class of maximal monotone interactions. Both the continuous and the discrete-time cases are analysed. Robustness is studied carefully. It is noteworthy that the results in this paper, extend naturally when the signum multifunction is replaced by a bounded maximal monotone mapping $\lambda_t \in \mathcal{M}(\dot{q} - r\dot{\theta})$. The case of hypomonotone friction (like Stribeck effects) is worth investigating: the addition of a linear feedback (using $u_1$) from both the mass and the pulley's velocities may be necessary to recover the monotone case. It would also be interesting to investigate whether or not the mass controller $u_1(q_d, \dot{q}_d, \ddot{q}_d)$ could be generated through $\lambda_t$, without any direct action $u_1$. Finally extension towards higher dimensional systems (with several masses on the belt, or stacked masses) could be tackled.

## A Relay Systems and LCS

Consider the dynamics in (1a), (1c) with $r\dot{\theta} = u_2$. Using [35] this is equivalently rewritten as:

$$ \left\{ \begin{array}{l} \dot{x} = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} x + \begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix} u + \begin{pmatrix} 0 \\ \mu g \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix} \lambda \\ 0 \le \lambda \perp w = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \lambda \\ \qquad + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x \ge 0, \end{array} \right. \quad (77) $$
with $x = (q, \dot{q})^\top$, $\lambda = (\lambda_1, \lambda_2)^\top$, $u = (u_1, u_2)^\top$. From the complementarity constraints: $\lambda_1 \in [0, 1]$. Same applies to (3). As alluded to above, uncertainties in $\mu$ do not affect multipliers $\lambda$. The LCS in (77) is compactly written as $\dot{x}(t) = Ax(t) + Eu(t) + E' + B\lambda(t)$, $0 \le \lambda \perp Cx(t) + D\lambda(t) + Fu(t) + F' \ge 0$. In (4) we have: $A_{cl} = A + EK$, $2\mu gB_{cl} = B$, $C_{cl} = C + FK$, $D_{cl} = D$. The multiplier $\lambda_1$ satisfies $2\mu g\lambda_1 \in \text{sgn}(\dot{q} - u_2)$. Hence $\lambda_1 \in [-\frac{1}{2\mu g}, \frac{1}{2\mu g}]$. The LCS form (77) shows that the controller acts in both the differential and the complementarity parts of the dynamics, so that it is natural to follow a passivity-based control strategy as in [27].

## B Analysis of Convergence of Discrete-time Solutions to Continuous-time Limits

This short section is a preliminary study of the convergence of piecewise-linear approximations (constructed from the discrete-time system's variables). Let us consider a countable sequence of positive timesteps $\{h_k\}_{k \ge 0}$, $t_{k+1} - t_k = h_k > 0$, and let $T > 0$ such that $T = \sum_{k=0}^n h_k$ for some $n \in \mathbb{N}$. Let us analyze the convergence of the continuous piecewise-linear functions:

$$ \bar{e}_n(t) \triangleq (e_{k+1} - e_k) \frac{t - t_k}{h_k} + e_k, \quad t \in [t_k, t_{k+1}), \quad (78) $$

defined on $[0, T]$, and $\{e_k\}_{k \ge 0}$ is a trajectory of the difference inclusion (48) (or of (49)). Thus we study the convergence of the sequences $\{\bar{e}_n\}_{n \in \mathbb{N}}$ on $[0, T]$ as $n \to +\infty$. It is clear from (78) and Proposition 6 that the sequence is uniformly bounded on $[0, T]$. Secondly consider $||\bar{e}_n(\tau_1) - \bar{e}_n(\tau_2)||$ for some $n \in \mathbb{N}$, $\tau_1, \tau_2 \in [0, T]$. From continuity of the functions $\bar{e}_n(\cdot)$, for any $\varepsilon > 0$ there exists $\delta > 0$ (that does not depend on $n$) such that $||\bar{e}_n(\tau_1) - \bar{e}_n(\tau_2)|| < \varepsilon$ for any $|\tau_1 - \tau_2| < \delta$. Thus the sequence $\{\bar{e}_n(\cdot)\}_{n \in \mathbb{N}}$ is uniformly equicontinuous. By the Arzelà-Ascoli Theorem, it follows that there exists a subsequence $\{\bar{e}_{n_i}\}_{i \in \mathbb{N}}$ which converges uniformly to a continuous limit $\bar{e}(\cdot)$ on $[0, T]$. Consider now the functions:

$$ \dot{\bar{e}}_n(t) = \frac{e_{k+1} - e_k}{h_k}, \quad t \in [t_k, t_{k+1}), \quad (79) $$

defined on $[0, T]$, which are the almost-everywhere derivatives of $\bar{e}_n(\cdot)$. From Proposition 6 it follows that $\dot{\bar{e}}_n(t)$ is bounded, hence the sequence $\{\dot{\bar{e}}_n(\cdot)\}_{n \in \mathbb{N}}$ is uniformly bounded on $[0, T]$. From the Banach-Alaoglu Theorem, this sequence converges weakly* towards a bounded limit $\dot{\bar{e}}(\cdot)$, i.e., for all $\varphi \in \mathcal{L}_1([0, T], \mathbb{R}^2)$: $\lim_{n \to +\infty} \int_0^T (\dot{\bar{e}}_n(t) - \dot{\bar{e}}(t))^{\top} \varphi(t) dt = 0$. Let $\bar{\lambda}_{t,n}^*(t) = \lambda_{t,k+1}$, and $\bar{\lambda}_{d,t,n}^*(t) = \lambda_{d,t,k+1}$, for all $t \in [t_k, t_{k+1})$. For identical reasons (both multipliers are bounded) it follows that $\{\bar{\lambda}_{t,n}^*(t)\}_{n \in \mathbb{N}}$, $\{\bar{\lambda}_{t,d,n}^*(t)\}_{n \in \mathbb{N}}$ converge weakly* to bounded limits $\bar{\lambda}_t(\cdot)$ and $\bar{\lambda}_{t,d}(\cdot)$, respectively, on $[0, T]$. Notice that same results hold for the sequences constructed from $\{\dot{\theta}_k\}_{k \in \mathbb{N}}$, i.e., $\bar{\theta}_n(t) \triangleq (\dot{\theta}_{k+1} - \dot{\theta}_k) \frac{t-t_k}{h_k} + \dot{\theta}_k$, $t \in [t_k, t_{k+1})$, and $\dot{\bar{\theta}}_n(t) = \frac{\dot{\theta}_{k+1} - \dot{\theta}_k}{h_k}$, $t \in [t_k, t_{k+1})$, see Remark 4.
The next step is to prove that limits are solutions of the closed-loop system (1) (11) (12), equivalently of (47) or (49). From (49) the following holds for all $t \in [t_k, t_{k+1})$:

$$ \dot{\bar{e}}_n(t) = \frac{e_{k+1} - e_k}{h_k} = A_{cl}\bar{e}_n^*(t) - \mu g E_{cl}(\bar{\lambda}_{t,n}^*(t) - \bar{\lambda}_{d,t,n}^*(t)), \quad (80) $$

where $\bar{e}_n^*(t) = e_k$, for all $t \in [t_k, t_{k+1})$. Mimicking [44, section 4.5], we have that $||\bar{e}_n - \bar{e}_n^*||_{L_2([0,T])} \le \frac{C^2T(\sup_k h_k)^2}{3}$ where $C$ is an upperbound of $\dot{\bar{e}}_n$. It is deduced that $\{\bar{e}_n^*(\cdot)\}_{n \ge 0}$ converges to $\bar{e}(\cdot)$ strongly in $L_2([0,T])$. Now $\lambda_{t,k} \in \text{sgn}(C_d e_k + \dot{q}_{d,k} - u_{2d,k})$, see (67b). Thus $\bar{\lambda}_{t,n}^*(t) \in \text{sgn}(C_d \bar{e}_n^*(t) + \dot{q}_{d,n}(t) - u_{2d,n}(t))$ for all $t \in [t_k, t_{k+1})$, where both sequences $\{\dot{q}_{d,n}(\cdot)\}_{n \ge 0}$ and $\{u_{2d,n}(\cdot)\}_{n \ge 0}$ are constructed as in (78). Desired trajectories and inputs are uniformly bounded, and we can assume that they satisfy regularity properties such that both sequences converge uniformly to continuous $q_d(\cdot)$ and $u_{2d}(\cdot)$ on $[0, T]$. Recalling that for a compact set $[0, T]: L_\infty([0, T]) \subset L_2([0, T]) \subset L_1([0, T])$ and that the sequences $\{\bar{\lambda}_{t,n}^*\}_{n \in \mathbb{N}}$ and $\{\bar{\lambda}_{d,t,n}^*\}_{n \in \mathbb{N}}$ converge in the weak$^*$ topology of $L_\infty([0, T])$, it follows that indeed both sequences also converge in the weak topology of $L_2([0, T])$ towards $\bar{\lambda}_t$ and $\bar{\lambda}_{d,t}$, respectively. Now, using [54, Proposition 20.33] we deduce that $\bar{\lambda}_t \in \text{sgn}(C_d \bar{e} + \dot{q}_d - u_{2d})$ and $\bar{\lambda}_{d,t} \in \text{sgn}(\dot{q}_d - u_{2d})$. Taking the limits of both sides in (80) it is inferred that the triplet of limit functions $(\bar{e}, \bar{\lambda}_t, \bar{\lambda}_{d,t})$ are solution of (49).

This shows that the approximate discrete-time solutions obtained using the plant discretization in (23), are "close" to the continuous-time differential inclusion.

## References

[1] M. Kunze, Non-smooth Dynamical Systems, Vol. 1744 of Lecture Notes in Mathematics, Springer-Verlag, 2000.

[2] R. I. Leine, H. Nijmeijer, Dynamics and Bifurcations of Non-Smooth Mechanical Systems, no. 18 in Lecture Notes in Applied and Computational Mechanics, Springer-Verlag, Berlin Heidelberg, 2004.

[3] G. Csernák, G. Stépán, On the periodic response of a harmonically excited dry friction oscillator, Journal of Sound and Vibration 295 (3) (2006) 649–658.

[4] J. P. Den Hartog, Forced vibrations with combined Coulomb and viscous friction, Transactions of the American Society of Mechanical Engineers 53 (2) (1930) 107–115.

[5] S. W. Shaw, On the dynamic response of a system with dry friction, Journal of Sound and Vibration 108 (1986) 305–325.

[6] B. Erickson, B. Birnir, D. Lavallée, Periodicity, chaos and localization in a Burridge-Knopoff model of an earthquake with rate-and-state friction, Geophysical Journal International 187 (2011) 178–198.
[7] U. Andreaus, P. Casini, Dynamics of friction oscillators excited by a moving base and/or driving force, Journal of Sound and Vibration 245 (2001) 685-699.

[8] N. Hinrichs, M. Oestreich, K. Popp, On the modelling of friction oscillators, J. Sound and Vibration 216 (3) (1998) 435-459.

[9] C. S. Liu, W. T. Chang, Frictional behaviour of a belt-driven and periodically excited oscillator, Journal of Sound and Vibration 258 (2) (2002) 247-268.

[10] K. Popp, P. Stelter, Stick-slip vibrations and chaos, Phil. Trans.: Physical Sciences and Eng. 332 (1624) (1990) 89-105.

[11] J. J. Thomsen, A. Fidlin, Analytical approximations for stick-slip vibration amplitudes, International Journal of Non-Linear Mechanics 38 (2003) 389-403.

[12] M. Oestreich, N. Hinrichs, K. Popp, Bifurcation and stability analysis for a non-smooth friction oscillator, Archive of Applied Mechanics 66 (5) (1996) 301-314.

[13] G. Csernak, G. Licsko, Asymmetric and chaotic responses of dry friction oscillators with different static and kinetic coefficients of friction, Meccanica 56 (2021) 2401-2414.

[14] M. A. Heckl, I. D. Abrahams, Active control of friction-driven oscillations, J. Sound and Vibration 193 (1) (1996) 416-426.

[15] J. Das, A. K. Mallik, Control of friction driven oscillation by time-delayed state feedback, Journal of Sound and Vibration 297 (3) (2006) 578-594.

[16] Z. Li, Q. Wang, H. P. Gao, Control of friction oscillator by lyapunov redesign based on delayed state feedback, Acta Mechanica Sinica 25 (2009) 257-264.

[17] M. Barreau, S. Tarbouriech, F. Gouaisbaut, Lyapunov stability analysis of a mass spring system subject to friction, Systems and Control Letters 150 (2021) 104910.

[18] A. Cabboi, L. Marino, A. Cicirello, A comparative study between Amontons-Coulomb and Dieterich-Ruina friction laws for the cyclic response of a single degree of freedom system, European J. Mechanics - A/Solids 96 (2022) 104737.

[19] Y. F. Liu, J. Li, Z. M. Zhang, X. H. Hu, W. J. Zhang, Experimental comparison of five friction models on the same test-bed of the micro stick-slip motion system, Mechanical Sciences 6 (2015) 15-28.

[20] G. Rill, T. Schaeffer, M. Schuderer, Lugre or not Lugre, Multibody System Dynamics 60 (2024) 191-218.
[21] M. Schuderer, G. Rill, T. Schaeffer, C. Schulz, Friction modeling from a practical point of view, Multibody System Dynamics https://doi.org/10.1007/s11044-024-09978-0 (2024).

[22] L. Marino, A. Cicirello, Experimental investigation of a single-degree-of-freedom system with Coulomb friction, Nonlinear Dynamics 99 (2020) 1781–1799.

[23] Z. Li, H. Wei, C. Liu, Y. He, G. Liu, H. Zhang, W. Li, An improved iterative approach with a comprehensive friction model for identifying dynamic parameters of collaborative robots, Robotica 42 (5) (2024) 1500–1522.

[24] P. Vigué, C. Vergez, S. Karkar, B. Cochelin, Regularized friction and continuation: Comparison with Coulomb's law, Journal of Sound and Vibration 389 (2017) 350–363.

[25] L. Marton, B. Lantos, Modeling, identification, and compensation of stick-slip friction, IEEE Transactions on Industrial Electronics 54 (1) (2007) 511–521.

[26] R. Beerens, H. Nijmeijer, W. P. M. H. Heemels, N. van de Wouw, Set-point control of motion systems with uncertain set-valued Stribeck friction, IFAC-PapersOnLine 50 (1) (2017) 2965–2970, 20th IFAC World Congress.

[27] A. Younes, F. A. Miranda-Villatoro, B. Brogliato, Trajectory tracking in linear complementarity systems with and without state jumps: A passivity approach, Nonlinear Analysis: Hybrid Systems 54 (2024) 101520.

[28] B. Brogliato, W. P. M. H. Heemels, Observer design for Lur'e systems with multivalued mappings: A passivity approach, IEEE Trans. Automatic Control 54 (8) (2009) 1996–2001.

[29] R. Lozano, B. Brogliato, Adaptive control of robot manipulators with flexible joints, IEEE Transactions on Automatic Control 37 (2) (1992) 174–181.

[30] B. Brogliato, R. Ortega, R. Lozano, Globally stable nonlinear controllers for flexible joint manipulators: a comparative study, Automatica 31 (7) (1995) 941–956.

[31] R. T. Rockafellar, Convex Analysis, Princeton University Press, New Jersey, 1970.

[32] J. J. Moreau, La notion de surpotentiel et les liaisons unilatérales en élastostatique, C. R. Acad. Sci. Paris Sér. A-B (1968) A954–A957.

[33] V. Sessa, L. Iannelli, F. Vasca, V. Acary, A complementarity approach for the computation of periodic oscillations in piecewise linear systems, Nonlinear Dynamics 85 (2016) 1255–1273.
[34] M. Legrand, C. Pierre, A compact, equality-based weighted residual formulation for periodic solutions of systems undergoing frictional occurrences, Journal of Structural Dynamics 2 (2024) 144–170.

[35] M. K. Camlibel, Complementarity Methods in the Analysis of Piecewise Linear Dynamical Systems, Ph.D. thesis, Katholieke Univ. Brabant, Tilburg, NL, https://pure.uvt.nl/ws/portalfiles/portal/419687/86166.pdf (2001).

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

[37] B. Brogliato, Modeling, analysis and control of robot object nonsmooth underactuated Lagrangian systems: A tutorial overview and perspectives, Annual Reviews in Control 55 (2023) 297–337.

[38] S. P. Bhat, D. S. Bernstein, Finite-time stability of continuous autonomous systems, SIAM Journal on Control and Optimization 38 (3) (2000) 751–766.

[39] B. Brogliato, R. Lozano, B. Maschke, O. Egeland, Dissipative Systems Analysis and Control, 3rd Edition, Communications and Control Eng., Springer Nature Switzerland AG, 2020.

[40] H. K. Khalil, Nonlinear Systems, 3rd Edition, Prentice Hall, Upper Saddle River, NJ, 2002.

[41] V. Acary, B. Brogliato, Y. Orlov, Chattering-free digital sliding-mode control with state observer and disturbance rejection, IEEE Transactions on Automatic Control 57 (5) (2012) 1087–1101.

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

[43] F. Miranda-Villatoro, B. Brogliato, F. Castanos, Multivalued robust tracking control of Lagrange systems: Continuous and discrete-time algorithms, IEEE Transactions on Automatic Control 62 (9) (2017) 4436–4450.

[44] F. Miranda-Villatoro, B. Brogliato, F. Castanos, Set-valued sliding-mode control of uncertain linear systems: continuous and discrete-time analysis, SIAM Journal on Control and Optimization 56 (3) (2018) 1756–1793.

[45] B. Andritsch, L. Watermann, S. Koch, M. Reichhartinger, J. Reger, M. Horn, Modified implicit discretization of the super-twisting controller, IEEE Transactions on Automatic Control 69 (8) (2024) 5620–5626.
[46] X. Xiong, R. Kikuuwe, S. Kamal, S. Jin, Implicit-Euler implementation of super-twisting observer and twisting controller for second-order systems, IEEE Transactions on Circuits and Systems II: Express Briefs 67 (11) (2020) 2607-2611.

[47] R. Kikuuwe, S. Yasukouchi, H. Fujimoto, M. Yamamoto, Proxy-based sliding mode control: A safer extension of PID position control, IEEE Transactions on Robotics 26 (4) (2010) 670-683.

[48] R. Kikuuwe, Y. Yamamoto, B. Brogliato, Implicit implementation of nonsmooth controllers to nonsmooth actuators, IEEE Transactions on Automatic Control 67 (9) (2022) 4645-4657.

[49] R. Seeber, B. Andritsch, Proper implicit discretization of the super-twisting controller without and with actuator saturation, Automatica 173 (2025) 112027.

[50] V. Acary, B. Brogliato, Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics, Vol. 35 of Lectures Notes in Applied and Computational Mechanics, Springer Verlag, 2008.

[51] F. Facchinei, J. S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer New York, NY, 2003.

[52] K. Addi, B. Brogliato, D. Goeleven, A qualitative mathematical analysis of a class of variational inequalities via semi-complementarity problems. Applications in electronics, Mathematical Programming 126 (2011) 31-67.

[53] L. Condat, D. Kitahara, A. Contreras, A. Hirabayashi, Proximal splitting algorithms for convex optimization: a tour of recent advances, with new twists, SIAM Review 65 (2) (2023) 375-435.

[54] H. H. Bauschke, P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces, Canadian Mathematical Society-Société mathématique du Canada, Springer Science+Business Media, 2011.

[55] F. A. Miranda-Villatoro, F. Castanos, B. Brogliato, When proximal-point algorithms meet set-valued systems. An optimization point of view of discrete-time sliding modes, Preprint hal-04626631, Univ. Grenoble-Alpes, Inria, https://inria.hal.science/hal-04626631v1/document (2024).

[56] V. Lakshmikantham, D. Trigiante, Theory of Difference Equations: Numerical Methods and Applications, Academic Press Inc., 1988.