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

HAL
open science

# Passivity-based Trajectory Tracking Control in Frictional Oscillators with Set-valued Friction

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

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

## To cite this version:

Aya Younes, Félix Miranda-Villatoro, Bernard Brogliato. Passivity-based Trajectory Tracking Control in Frictional Oscillators with Set-valued Friction. Nonlinear Analysis: Hybrid Systems, 2026, 60, pp.article 101672. 10.1016/j.nahs.2025.101672. hal-04810905v2

HAL Id: hal-04810905

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

Submitted on 22 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Passivity-based Trajectory Tracking Control in Frictional Oscillators with
Set-valued Friction *

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

<sup>a</sup>Univ. Grenoble-Alpes, Inria, CNRS, LJK, Grenoble INP, Grenoble, 38000, France

**Abstract**

This article is largely concerned with the trajectory tracking control of frictional oscillators, which are nonsmooth nonlinear dynamical systems. The trajectory tracking problem, which is studied under a passivity-based controller addresses three main cases: the nominal case with known friction coefficient, uncertain friction coefficient, and when the Coulomb friction model is enhanced by including Stribeck effects. Monotonicity (or hypomonotonicity) of the friction model is crucial for the stability analysis of the tracking error. It can be relaxed to hypomonotonicity to handle Stribeck model. The framework of linear complementarity systems is used for the analysis. The case of a two-mass system is tackled as an extension of the standard one-mass oscillator. Theoretical results are supported by numerical simulations.

**Keywords:** Linear complementarity systems, passivity, Lyapunov stability, mechanical oscillators, Coulomb friction, set-valued systems, nonsmooth dynamics, trajectory tracking, Stribeck effects, maximal monotonicity.

**1. Introduction**

Frictional oscillators are low-dimension mechanical systems where friction plays a prominent role due to its non-linear and nonsmooth effects. They have gained significant attention in the nonlinear dynamics, bifurcation and chaos theory scientific community due to their applications in mechanical engineering (like energy dissipation mechanisms to damp vibrations). One-degree-of-freedom oscillators have been widely studied in the applied mathematics and in the bifurcation and chaos scientific communities, see, e.g., [1, Chapter 9] [2, 3, 4, 5, 6, 7, 8, 9]. These oscillators may also be linked to a simplified Burridge-Knopoff model in seismology, highlighting their practical applications (see [10]). Many other applications exist in turbine blade joints, electric motor drives, wheel rail coupling of mass transit systems, and brake systems [11]. Many authors studied the behaviour of the stick-slip oscillators with various frictional laws see, e.g., [2, 3, 4, 12, 13, 14, 15, 16, 17, 18]. Stability analysis of periodic solutions for periodically excited frictional oscillator were first presented by [6, 7]. Later studies, see, e.g., [16, 19, 20, 15, 5], further investigated the conditions for the existence and the stability of periodic solutions. Externally excited frictional oscillators, as studied in [17, 18, 12, 13], may exhibit chaotic behaviour. These systems often undergo rich bifurcation behaviour, which is explored in [21, 22, 3, 4, 23, 24, 13, 19, 20]. Stick-slip vibrations, which are self-sustained oscillations induced by dry friction are described in [18, 16, 25]. Various control techniques have been developed for the stability analysis of frictional oscillators, including active feedback control to superimpose frictional force in [26], and other control strategies based on delayed position feedback [27, 28, 29, 30, 31, 32]. One-degree-of-freedom oscillator may be used as simplified models of feed drive assembly for grinding machines [33, Figure 3], and for the control of machine tools [34, 35, 36]. Extensions to nonlinear oscillators [37] and coupled masses [38, 39, 40, 41, 11] have been considered.

*The work of the first author was supported by the French National Research Agency in the framework of the "Investissements d'avenir" program "ANR-15-IDEX-02" and the LabEx PERSYVAL "ANR-11-LABX-0025-01"

Email addresses: aya.younes@inria.fr (Aya Younes), felix.miranda-villatoro@inria.fr (Félix Miranda-Villatoro),
bernard.brogliato@inria.fr (Bernard Brogliato)
## 1.1. Friction Modeling

Friction is the main source of difficulty in the analysis of frictional oscillators, making its modeling an important step. Given that there is a plethora of friction models, analysis and control of systems with friction have witnessed a significant research interest over the past thirty years both in automatic control and nonlinear dynamics. Friction models of interest for control may be classified into static models and dynamic models. Static friction includes set-valued or single-valued models. Set-valued models such as Coulomb's friction with constant coefficient or varying coefficient (e.g., Stribeck [42, 43], Dieterich-Ruina [44, Equation (21)]), allow to correctly handle sticking modes at zero relative tangential velocity and are thus essential in multibody applications, contrarily to regularized models which are mainly applied to planar friction where the signum function is replaced by a saturation, or a sigmoid function. Dynamic models with internal state are numerous (LuGre, Leuven, Dahl, Bliman-Sorine, elastoplastic etc [43]). Dynamic models are better for modeling micro stick-slip effects and pre-sliding motions, though Coulomb with Stribeck coefficient of friction can provide good results in this respect [43]. However, dynamic models also present severe drawbacks [45, 46]. Coulomb friction (with constant coefficient) has limitations [44], but it is able to give quite good results in practice [47, 48, 49, 50, 51], since they allow to better capture physical phenomena like sticking modes [52]. This is the reason why the approach we choose in this article is to design the control loop using Coulomb's model, then analyse robustness with respect to uncertain coefficients, then add some well-known effects like Stribeck's effects and investigate how the proposed approach extends. It is sometimes not easy task to choose between a regularized (around zero tangential velocity) model and the set-valued model [53]. See [54, 55, 56] for surveys and comparisons of friction models for simulation, identification, stability and control.

## 1.2. Control of Systems with Friction

Several classes of control strategies have been developed for controlling systems with friction. A popular approach in the control and robotics literature is to use model-based friction compensator in the controller. The compensation method (which, roughly speaking, consists of annihilating frictional forces using a model $\hat{F}_{fric}$ of the friction force $F_{fric}$, then applying some controller designed for frictionless systems) can be implemented through feedback control which is based on the actual system states or through feedforward control which is based on the desired system states, or using internal-state observers and identification schemes, see [57, 58, 47, 59, 60, 61, 62, 63, 64, 65, 66]. See, e.g., [67, 68] for applications in robotics, and [65] for a survey on model-based and model-free (neural networks, iterative learning, fuzzy logic) methods. Viscous and/or dry friction observers are proposed in [69, 70], and are used with P, PD or PI controllers, with experimental validations. However, since it is challenging to model friction exactly, the compensation approach has drawbacks due to error $\tilde{F}_{fric} = \hat{F}_{fric} - F_{fric} \neq 0$: overcompensating friction leads to limit cycles and undercompensating friction leads to steady state errors, see [71] for an analysis of the compensation methods. An adaptive friction compensation control is implemented without a complete model of friction in [68, 72]. Another class of control strategies, not based on friction compensation with an as-accurate-as-possible friction model $\hat{F}_{fric}$, is proposed in [73, 74, 42, 75, 76, 51]. Roughly speaking, such methods do not compensate-then-control, but they finely analyse stability using some friction model $\hat{F}_{fric}$ (considered to be an acceptable model in this case). An impulsive control is designed in [76]. PID controllers are studied in [75, 66]: they guarantee Lyapunov asymptotic stability of a set of sticking fixed points. Stabilization of the origin is tackled in [74, 77, 78], using PID control with switching integral action. A one-degree-of-freedom system with Stribeck friction is studied in [42]. Local stability of attractors and of sliding equilibria with basins of attraction is analysed. In [51], stability properties of a one-degree-of-freedom system with a PI feedback are analysed (the closed-loop system is a relay system with relative degree 1). Experimental validations are shown in [51]. In [79] a robust sliding-mode twisting sliding-mode set-valued controller is used to asymptotically stabilize the origin: friction is considered as a bounded unknown disturbance. As alluded to above, another class of feedback controllers consists of proportional position inputs with time delays [27, 28, 29, 30, 31, 80, 32]. They are mainly used to create regions of locally stable equilibria, hence avoiding unwanted vibrations.

## 1.3. Article's contributions and organization

In this article, the friction is modeled by a set-valued signum function, which represents planar Coulomb's friction with constant normal contact force. This function is quite amenable for complementarity [81, 82], allowing the study of frictional oscillators within the framework of linear complementarity systems. This connection to complementarity
systems motivates us to apply some theoretical results from our previous work on trajectory tracking in linear complementarity systems [83] and providing original extensions in this paper. To the best of our knowledge, trajectory tracking for frictional oscillators was not addressed before. This work introduces a feedback control design approach that incorporates friction, where the controller is embedded inside the signum function. The control strategy does not aim at compensating for friction, then applying some suitable controller to the frictionless system. On the contrary friction is pivotal to the feedback control design, which uses Coulomb's friction specific features in a suitable way to leverage the linear collocated control action. Then robust trajectory tracking is possible.

The main contribution of this article is to tackle the trajectory tracking control problem in frictional oscillators with set-valued friction, considering different cases and challenges. The stability result derived from our previous work in [83] is applied to the nominal case in this article, with a new result concerning the complementarity variables. Incremental passivity is used in the stability proof. In the case of uncertainties, the controller design and stability analysis apply the theoretical results in [83]. In addition, this article introduces a new approach for the control design based on rendering the hypomonotone friction model when Stribeck effects are incorporated, maximal monotone. The control strategy is extended to design a trajectory tracking controller for higher dimension systems. The center of mass of a frictional oscillator with two masses, as well as Lagrangian systems with several frictional contacts, are studied.

This article is organised as follows: section 2 provides useful notations and definitions; Section 3 introduces the dynamics of the simple frictional oscillator within the framework of linear complementarity systems, it addresses well-posedness issues and the characterization of equilibrium points of the error dynamics. The analysis of the desired trajectories is presented in Section 4, followed by the passivity-based control design and the closed-loop error stability analysis for the known parameters case in Section 5. Section 6 is dedicated to study robustness in the presence parametric uncertainties. Section 7 investigates the stability analysis with the Stribeck friction model. Section 8 is dedicated to the tracking control of higher dimensional systems with multiple frictional contacts. To illustrate theoretical results, numerical simulations are provided throughout the article. Conclusions end the article in Section 9, and several useful tools and results are provided in Appendix A through Appendix D.

## 2. Preliminaries

This section introduces the notation and definitions necessary for the analyses made in this article. The notation $\langle \cdot, \cdot \rangle$ denotes the standard Euclidean inner product in $\mathbb{R}^n$. The elements of a matrix $M \in \mathbb{R}^{n \times m}$ are denoted $M_{ij}$. $M_{\bullet, j}$ is its $j$th column, and $M_{i, \bullet}$ is its $i$th row. A matrix $M \in \mathbb{R}^{n \times n}$, possibly nonsymmetric, is said to be positive definite, $M > 0$ (resp. positive semi-definite, $M \succeq 0$) if $x^\top M x > 0$ for all $x \neq 0 \in \mathbb{R}^n$ (resp. $x^\top M x \geq 0$). $M$ is a P-matrix if all its principal minors are positive; a positive definite matrix is a P-matrix. The minimum and the maximum eigenvalues of $M$ are denoted $\lambda_{\min}(M)$ and $\lambda_{\max}(M)$. The maximum singular value of $M \in \mathbb{R}^{n \times m}$ is denoted as $\sigma_{\max}(M)$, $\mathrm{Im}(M)$ denotes its image, and $\mathrm{Ker}(M)$ its null space. The singular values are ordered as $\sigma_{\max}(M) = \sigma_1(M) \geq \sigma_2(M) > \dots > \sigma_r(M) > 0$, where $r = \mathrm{rank}(M)$. Let $M = M^\top$, its eigenvalues are real and ordered as $\lambda_{\max}(M) = \lambda_1(M) \geq \lambda_2(M) \geq \dots \geq \lambda_n(M) = \lambda_{\min}(M)$. $M^\dagger$ denotes the Moore-Penrose pseudo-inverse of $M$. The range of the function $f(\cdot)$ is $\mathrm{Im}(f)$, its domain is $\mathrm{dom}(f) = \{x | f(x) < +\infty\}$. Right and left limits of the function $f(\cdot)$ at $t$ are denoted $f(t^+)$ and $f(t^-)$. Let $S \subseteq \mathbb{R}^m$ be a closed convex set containing 0. Its polar cone $S^\circ = \{x \in \mathbb{R}^n | x^\top v \leq 0, \forall v \in S\} = -S^\star$ where $S^\star$ is the dual cone. AC stands for absolute continuity. Let $f: \mathbb{R}^m \to \mathbb{R}^n$, its Jacobian matrix is $\frac{\partial f}{\partial x} \in \mathbb{R}^{m \times n}$. Let us recall some useful definitions [84, 85, 86, 87, 88].

**Definition 1.** Let $S \subset \mathbb{R}^n$ be a closed, non-empty and convex set. The normal cone to $S$ at $x \in S$ is given by

$$ N_S(x) = \{g \in \mathbb{R}^n \mid \langle g, z - x \rangle \leq 0, \forall z \in S\}. \quad (1) $$

If $x \notin S$, then $N_S(x) = \emptyset$.

The domain of a set-valued mapping $T: \mathbb{R}^n \rightrightarrows \mathbb{R}^n$ is $\mathrm{Dom}(T) = \{x \in \mathbb{R}^n \mid T(x) \neq \emptyset\}$, whereas $\mathrm{Gph} T := \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^n \mid y \in T(x)\}$ denotes the graph of $T$.

**Definition 2.** A multivalued mapping $T: \mathrm{Dom}(T) \subseteq \mathbb{R}^n \rightrightarrows \mathbb{R}^n$ is monotone if for all $u_1, u_2 \in \mathrm{Dom}(T)$, for all $v_1 \in T(u_1)$, $v_2 \in T(u_2)$

$$ \langle u_1 - u_2, v_1 - v_2 \rangle \geq 0. $$
It is hypomonotone if there exists a real $k > 0$ such that: $\langle u_1 - u_2, v_1 - v_2 \rangle \ge -k\|u_1 - u_2\|^2$. In addition, $T$ is maximal monotone if it is monotone and its graph cannot be enlarged without destroying monotonicity.

**Definition 3.** Let $f: \mathbb{R}^n \to \mathbb{R}$ be a proper continuous convex function. A vector $g \in \mathbb{R}^n$ is called a subgradient of $f$ at a point $\bar{x}$ if, for all $x \in \mathbb{R}^n$: $f(x) \ge f(\bar{x}) + g^\top(x - \bar{x})$. The subdifferential of $f$ is the set of all subgradients and is denoted by $\partial f(\bar{x})$.

**Definition 4.** A linear complementarity system (LCS) with inputs is a nonsmooth nonlinear dynamical system defined as:

$$ \begin{cases} (a) \quad \dot{x}(t) = Ax(t) + B\lambda(t) + Eu(t) \quad (\text{dt-almost everywhere}) \\ (b) \quad 0 \le \lambda(t) \perp w(t) = Cx(t) + D\lambda(t) + Fu(t) \ge 0 \end{cases} \quad (2) $$

where $x(t) \in \mathbb{R}^n$, $\lambda(t) \in \mathbb{R}^m$ and $w(t) \in \mathbb{R}^m$, for all $t \ge 0$, A, B, C, D, E and F are constant matrices of appropriate dimensions, $u: \mathbb{R}_+ \to \mathbb{R}^p$ is an exogenous signal or a control input. The constraint in (2) (b) is the short form of the following three conditions: i) $\lambda(t) \in \mathbb{R}_+^m$, ii) $w(t) = Cx(t) + D\lambda(t) + Fu(t) \in \mathbb{R}_+^m$, and iii) $\langle \lambda(t), w(t) \rangle = 0$. Such constraints describe a linear complementarity problem (LCP) in $\lambda$, denoted LCP(D, Cx + Fu).

The LCP (2) (b) states that for $i = 1, \dots, k$, either $\lambda_i = 0$ or $w_i = 0$. The set of constraints corresponding to the indices where $w_i = 0$ is the set of active constraints and it varies as the system evolves over time. Thus, the LCS in (2) can transition from one operation mode to another. A state $x^*$ is an equilibrium point of (2) if and only if there exist $\lambda^*$ and $w^* \in \mathbb{R}^m$ such that the mixed LCP (MLCP)

$$ \begin{cases} 0 = Ax^* + B\lambda^* + Eu(t) \\ 0 \le \lambda^* \perp w^* = Cx^* + D\lambda^* + Fu(t) \ge 0 \end{cases} \quad (3) $$

holds for all $t \ge 0$. It is noteworthy that $\lambda^*$ may be time-varying.

State jumps, or discontinuities in the solution $x(\cdot)$ occur in LCS (2) at initial time or at times of discontinuity in $Fu(\cdot)$. The formulation of state jumps is based on dissipation principles and is given in [84, Lemma 2.3], [89, 90, 91, 92]. Let us introduce the jump rule briefly.

**Definition 5.** [92, Definitions 6, 7] The set of admissible state jumps for LCS in (2) is defined as:

$$ \mathcal{K} \stackrel{\Delta}{=} \{z \in \mathbb{R}^n \mid Cz + Fu(t^+) \in Q_D^*\} \quad (4) $$

where $Q_D = \{z \in \mathbb{R}^m \mid 0 \le z \perp Dz \ge 0\}$ and $Q_D^*$ is its dual cone. If $(A, B, C, D)$ is passive with $P = P^\top > 0$ and if $0 \in \mathcal{K}$, then the jumps dissipate energy (i.e., $V(x(t^+)) \le V(x(t^-))$) [92, Lemma 3] and the state jump rule is given by:

$$ x(t^+) = \underset{x \in \mathcal{K}}{\operatorname{argmin}} \frac{1}{2} (x - x(t^-))^\top P (x - x(t^-)). \quad (5) $$

**Definition 6.** A quadruple $(A, B, C, D)$ is said to be passive, or dissipative with respect to the supply rate $\lambda^\top w$, if there exists a non-negative function $V: \mathbb{R}^n \to \mathbb{R}_+$, called a storage function, such that for all $t_0 \le t_1$ and all functions $(\lambda, x, w) \in \mathcal{L}_2([t_0, t_1]; \mathbb{R}^m \times \mathbb{R}^m)$ such that $\dot{x}(t) = Ax(t) + B\lambda(t)$, $w(t) = Cx(t) + D\lambda(t)$, the following inequality holds:

$$ V(x(t_0)) + \int_{t_0}^{t_1} \lambda^\top(t)w(t)dt \ge V(x(t_1)) \quad (6) $$

The inequality (6) is called the dissipation inequality. Equivalently, the following linear matrix inequality (LMI)

$$ \begin{pmatrix} A^\top P + PA & PB - C^\top \\ B^\top P - C & -D - D^\top \end{pmatrix} \preceq 0 \quad (7) $$

has a solution $P = P^\top \succeq 0$. Then $V(x) = \frac{1}{2}x^\top P x$. The quadruple $(A, B, C, D)$ is said strictly state passive if the LMI (7) holds with $P = P^\top > 0$ and $A^\top P + PA + \epsilon P$ for some constant $\epsilon > 0$. In this case, LMI (7) is named SSPMI(A, B, C, D, $\epsilon$). It is said strongly passive if the LMI (7) holds with strict inequality and $P = P^\top > 0$.
### 3. Controlled Frictional Oscillators

<figure><img src="image_4.png" /><figcaption>Figure 1: Simple frictional oscillator with Coulomb's friction.</figcaption></figure>

Consider the simple frictional oscillator system depicted in Figure 1a. The objective is to perform trajectory tracking for the mass $m$, using as inputs both $u_1$ and the belt's velocity through friction exerted by the belt on the mass. Assuming that the belt can be given arbitrary velocity, the dynamics of the oscillator are given by:

$$ \ddot{q}(t) \in \frac{1}{m}u_1(t) - \mu g \operatorname{sgn}(\dot{q}(t) - u_2(t)) \quad (8) $$

with $u = (u_1, u_2)^\top$, $u_2 = r\omega$ is considered as a control input, $\mu > 0$ is the friction coefficient, $\omega$ rad/s is the pulleys' angular velocity, $r$ is their radius, $g$ is the gravity acceleration, $m$ in kg is the mass weight. The sign multifunction sgn: $R \Rightarrow [-1, 1]$ is defined as sgn($q$) = +1, if $q > 0$; sgn($q$) = -1, if $q < 0$; and sgn(0) = [-1, 1]. Notice that the dynamics (8) can be rewritten as $\ddot{q}(t) = \frac{1}{m}u_1(t) - \mu g \lambda_1(t)$, with $\lambda_1(t) \in \operatorname{sgn}(\dot{q}(t) - u_2(t))$. If $u_1 = 0$ or if $u_1 = u_1(q, \dot{q}, t)$, the dynamics (8) fit within the robot-object class of systems [93] where the robot's dynamics (the pulley-belt system) is neglected: the mass $m$ is the object which can be controlled only through the contact force multiplier $\lambda_1(\cdot)$. Neglecting the friction 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, 6, 7, 32, 94], see also references therein. In this setting, we may link the mass to a rigid wall with a linear spring-dashpot model as shown in Figure 1b, yielding the dynamics:

$$ \ddot{q}(t) \in \frac{1}{m}u_1(t) - \frac{k}{m}(q(t) - l) - \frac{f}{m}\dot{q}(t) - \mu g \operatorname{sgn}(\dot{q}(t) - u_2(t)) \quad (9) $$

where $k > 0$ is the stiffness, $f > 0$ is the viscous friction coefficient, and $l$ is the natural length of the spring. Under these conditions, the analysis made in the sequel shows that $u_1$ can be chosen as a feedforward controller that depends only on exogenous signals. It is also noteworthy that $u_1$ could be chosen as the displacement imposed to one tip of the spring, changing $k(q - l)$ to $k(q + u_1 - l)$ in (9). Provided that $k$ is known, this does not modify the analysis made in the next sections. Notice that rotational systems like clutch systems have the same dynamics as (8), with displacements replaced by angles, see Fig. 2. Its dynamics reads as: $I\dot{\theta}(t) \in -k(\theta(t) - \theta_0) - f\dot{\theta}(t) + u_1(t) - \mu |F_n| \operatorname{sgn}(\dot{\theta}(t) - \omega(t))$, where $I$ is the moment of inertia of the left body, $F_n$ is a normal force exerted on the right body to keep contact closed, $k \ge 0$ is the torsional stiffness of a torsional spring acting on the left shaft, $f\dot{\theta}$ is a viscous friction at the revolute joint, $u_1$ is a control torque applied to the left shaft, $\omega = u_2$ is the angular velocity of the right shaft. Coulomb's model may be a first approximation of the frictional process between both discs. See also [95] for a rotating conveyor-belt analysis which shares similar dynamics [95, Equation (1)]. Two stacked blocks sliding on each other, where one of them is velocity-controlled and whose dynamics is not influenced by the other block, also fit with same framework.

**Remark 3.1.** One important assumption that is considered here, is that the angular velocity of the pulley $\omega$ can be set arbitrarily. In this way, $u_2 = r\omega$ is an input to be specified below, (since $r$ is assumed constant). Thus, with this assumption in place, the dynamics of the actuator in the pulley (usually, an electrical motor) is not considered, simplifying the design of the closed-loop. Taking into account the actuator dynamics would require an extra control-loop, (based for instance in backstepping control), for the implementation of the controller $u_2$, which is not considered
<figure><img src="image_5.png" /><figcaption>Figure 2: A "clutch" oscillator.</figcaption></figure>

here due to the space limitation. Notice also that considering $\omega$ as an independent variable is a common assumption in the nonlinear dynamics literature.

**Remark 3.2.** There certainly exists several ways to achieve trajectory tracking for (8). In this work, we do not aim to compensate for friction using the control input $u_1$, a control strategy adopted in many works (see section 1.2). Such compensation is also achievable using for instance a sliding-mode controller, where the friction term $\mu g \lambda_t$ would be considered as a bounded disturbance. Our strategy is based on using both inputs $u_1$ and $u_2$ (hence this is a two-channel system as in Fig. 9), relying on the monotonicity of friction coupled to passivity of a suitable closed-loop operator for stability analysis and using a collocated input $u_1$. Then the proposed strategy may hold with a simple feedforward input $u_1(t)$ when $f > 0$ and $k > 0$.

### 3.1. The Complementarity System Framework

The dynamics (8) are equivalently rewritten in a complementarity systems framework as:

$$ \left\{ \begin{array}{l} \dot{x}(t) = \underbrace{\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}}_{\triangleq A} x(t) + \underbrace{\begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix}}_{\triangleq E} u(t) + \underbrace{\begin{pmatrix} 0 \\ \mu g \end{pmatrix}}_{\triangleq E'} + \underbrace{\begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix}}_{\triangleq B} \lambda(t) \\ 0 \le \lambda(t) \perp w(t) = \underbrace{\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}}_{\triangleq F} u(t) + \underbrace{\begin{pmatrix} 0 \\ 1 \end{pmatrix}}_{\triangleq D} + \underbrace{\begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}}_{\triangleq C} \lambda(t) + \underbrace{\begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix}}_{\triangleq C} x(t) \ge 0, \end{array} \right. \quad (10) $$

with $x = (q, \dot{q})^\top$, $\lambda = (\lambda_1, \lambda_2)^\top$, $u = (u_1, u_2)^\top$. The complementarity conditions in (10) stem from the equivalent representation of the set-valued signum function in [81]. It is verified that the quadruple $(A, B, C, D)$ is passive with storage matrix $P = \begin{pmatrix} p_{11} & 0 \\ 0 & \frac{1}{2g} \end{pmatrix}$. There are two peculiarities to system (10): a) $D \ge 0$ and $D + D^\top = 0$, b) there are constant terms $\begin{pmatrix} 0 \\ 1 \end{pmatrix} \notin \text{Im}(F)$ and $E' = \begin{pmatrix} 0 \\ -\mu g \end{pmatrix} \in \text{Im}(E)$, hence it does not fit with (2). Solutions to (10) are AC under some conditions on the control input, see below. The first step is the computation of fixed points of (10) in closed-loop, in order to determine whether or not a feedback controller of the form:

$$ u(x, \lambda, t) = K[x - x_d(t)] + G[\lambda - \lambda_d(t)] + u_d(t). \quad (11) $$

where $x_d$, $\lambda_d$, $u_d$ are defined later, and with feedback gains $K \in \mathbb{R}^{2 \times 2}$ and $G \in \mathbb{R}^{2 \times 2}$, is able to modify the set of equilibria. Notice that $\lambda_1 = \frac{\lambda_1+1}{2}$ where $\mu g \lambda_1$ is the tangent interaction force between the mass and the belt, i.e., $\lambda_1 \in \text{sgn}(\dot{q} - u_2)$. The multiplier $\lambda_2$ is a slack variable has no clear mechanical meaning. It is inferred that only
$\lambda_1$ is available for feedback. Moreover, we avoid the use of $\lambda_1$ in $u_1$, for otherwise the tangential force could be compensated for directly: this is not a desirable control strategy, see Remark 3.2. Most importantly, as alluded to above, our analysis applies to the case when the mass is attached to a rigid wall by a linear spring-dashpot system as in Figure 1b, while $u_1$ is to be interpreted as a feedforward controller, that depends only on $x_d(t)$ and $u_d(t)$. In such a setting, the tracking control problem is that of a mass+spring-dashpot system, with feedforward control $u_1$ and feedback control $u_2$ through the tangential contact force which is a maximal monotone operator. These observations constrain the input matrix $G$ to satisfy $G = \begin{pmatrix} 0 & 0 \\ g_{21} & 0 \end{pmatrix}$. It is also noticed that a feedback using $\lambda_2$ is useless as it is unable to modify the closed-loop matrix $D_{cl}$ to get strong passivity ($\Rightarrow D_{cl} > 0$) instead of strict state passivity, because of the structure of $F$. In general, it is even better to assume $G = 0$. However, the subsequent analysis shows that $g_{21} > 0$ allows to regularize the set-valued friction and to get robustness that is otherwise impossible to obtain. Inserting (11) to (10), with the above restrictions on $G$, gives the closed-loop system dynamics:

$$ \left\{ \begin{array}{l} (a) \ \dot{x}(t) = \underbrace{\begin{pmatrix} 0 & 1 \\ \frac{k_{11}}{m} & \frac{k_{12}}{m} \end{pmatrix}}_{\triangleq A_{cl}} x(t) + \underbrace{\begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix}}_{\triangleq B_{cl}} \lambda(t) - \begin{pmatrix} \frac{k_{11}}{m} & \frac{k_{12}}{m} \end{pmatrix} x_d(t) + \begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ \mu g \end{pmatrix} \\ \\ (b) \ 0 \le \lambda(t) \perp w(t) = \underbrace{\begin{pmatrix} k_{21} & -1+k_{22} \\ 0 & 0 \end{pmatrix}}_{\triangleq C_{cl}} x(t) + \underbrace{\begin{pmatrix} g_{21} & 1 \\ -1 & 0 \end{pmatrix}}_{\triangleq D_{cl}} \lambda(t) - \begin{pmatrix} k_{21} & k_{22} \\ 0 & 0 \end{pmatrix} x_d(t) \\ \qquad - \begin{pmatrix} g_{21} & 0 \\ 0 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} u_{2d}(t) \\ 1 \end{pmatrix} \ge 0 \end{array} \right. \quad (12) $$

**Remark 3.3.** No physical parameters appear in the complementarity conditions in (10) or (12), because the complementarity merely represents the set-valued signum function. This will have consequences on the robustness analysis (see section 6). Again it is noteworthy that even a feedback from $\lambda_1$ cannot guarantee strong passivity of the closed-loop, because at best $D_{cl} \ge 0$. But strict state passivity is possible. One question which arises is how much $g_{21} \ne 0$ can improve the robustness with respect to uncertainties in $\mu$. This issue is tackled in section 6.

Following the approach in [83], let us choose the system which generates the desired trajectories as:

$$ \left\{ \begin{array}{l} \dot{x}_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} x_d(t) + \begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ \mu g \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix} \lambda_d(t) \\ 0 \le \lambda_d(t) \perp w_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) \ge 0, \end{array} \right. \quad (13) $$

where $x_d = (q_d, \dot{q}_d)^\top = (x_{1d}, x_{2d})^\top$. This LCS is equivalent to its differential inclusion form:

$$ \ddot{q}_d(t) \in \frac{1}{m} u_{1d}(t) - \mu g \lambda_{t,d}(t), \quad \lambda_{t,d}(t) \in \operatorname{sgn}(\dot{q}_d(t) - u_{2d}(t)) \quad (14) $$

This gives rise to the closed-loop error dynamics $e = x - x_d$:

$$ \left\{ \begin{array}{l} \dot{e}(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix} K e(t) + \begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix} + (\lambda(t) - \lambda_d(t)) \\ 0 \le \lambda_d(t) \perp w_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) \ge 0 \\ 0 \le \lambda(t) \perp w(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} K + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} e(t) - \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} G \lambda_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} \\ \qquad + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} G \lambda(t) + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) + \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) \ge 0, \end{array} \right. \quad (15) $$
that is:

$$
\left\{
\begin{array}{l}
\dot{e}(t) = \overbrace{\begin{pmatrix} 0 & 1 \\ \frac{k_{11}}{m} & \frac{k_{12}}{m} \end{pmatrix}}^{A_{cl}} e(t) + \overbrace{\begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix}}^{B_{cl}} (\lambda(t) - \lambda_d(t)) \\
0 \le \lambda_d(t) \perp w_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) \ge 0 \\
0 \le \lambda(t) \perp w(t) = \overbrace{\begin{pmatrix} k_{21} & k_{22} - 1 \\ 0 & 0 \end{pmatrix}}^{C_{cl}} e(t) - \begin{pmatrix} g_{21} & 0 \\ 0 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \overbrace{\begin{pmatrix} g_{21} & 1 \\ -1 & 0 \end{pmatrix}}^{D_{cl}} \lambda(t) \\
\qquad + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) + \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) \ge 0
\end{array}
\right.
\qquad (16)
$$

with $G = \begin{pmatrix} 0 & 0 \\ g_{21} & 0 \end{pmatrix}$. The closed-loop LCP matrix is $D_{cl}$.

**Remark 3.4.** It is important to note that rank($B$) = rank($B_{cl}$) = 1. So, the passification conditions on the rank of the input matrix [88, Theorem 3.61] are not met in this context. It is also noteworthy that $A_{cl} + A_{cl}^{\top}$ cannot be negative definite. Hence the LMI in (7) applied to $(A_{cl}, B_{cl}, C_{cl}, D_{cl})$ can be satisfied only for $P \neq I_2$. More generally, assume that there is no control inside the signum multifunction, i.e., $k_{21} = k_{22} = 0$, and assume also that $g_{21} = 0$. Then $D_{cl} + D_{cl}^{\top} = 0$, hence necessarily $PB_{cl} = C_{cl}^{\top}$ which holds for all $P = \text{diag}(p_{11}, \frac{1}{2\mu g})$, $p_{11} > 0$. With such $P$, the matrix $PA_{cl} + A_{cl}^{\top}P$ cannot be negative definite. It is therefore necessary to consider a controller $u_2$ acting inside the signum multifunction (or, more generally, inside a maximal monotone set-valued mapping) to render the passivity LMI (7) solvable with strict state passivity ($\epsilon > 0$). It is inferred that the trajectory tracking problem tackled in this paper, can be solved only by considering a control acting through friction between the belt and the mass.

## 3.2. Well-posedness of the Closed-loop Dynamics

The well-posedness of (12) and (13) has to be guaranteed as a prerequisite to stability analysis. Let us state the following:

**Assumption 3.5.** The quadruple $(A_{cl}, B_{cl}, C_{cl}, D_{cl})$ in (12) is strictly state passive.

The satisfaction of Assumption 3.5 is equivalent to the existence of gain matrices $K$ and $G$ such that the following SSPMI with unknowns $P$, $K$ and $G$:

$$ M \triangleq \begin{pmatrix} (A + EK)^{\top} P + P(A + EK) & P(B + EG) - (C + FK)^{\top} \\ (B + EG)^{\top} P - (C + FK) & -(D + FG) - (D + FG)^{\top} \end{pmatrix} \preceq \begin{pmatrix} -\epsilon P & 0 \\ 0 & 0 \end{pmatrix}, \quad (17) $$

has a solution with $P = P^{\top} > 0$. This implies that $g_{21} \ge 0$. As noted in Remark 3.4, it is necessary to consider $u_2$ to ensure the solvability of the LMI (17) when $G = 0$. As examples will show, strict state passivity can hold with $g_{21} = 0$ in certain cases, which may be interesting in practice. Let us examine the well-posedness of the closed-loop system in (12) and of the desired dynamics in (13). First, let us verify that no state jump can occur in the closed-loop system.

Indeed, let us consider (5), we obtain: $Q_{D_{cl}} = \begin{pmatrix} 0 \\ R_+ \end{pmatrix}$, $Q_{D_{cl}}^{*} = \begin{pmatrix} R \\ R_+ \end{pmatrix}$, and $K_{cl} = R^2$. Thus, $x(t^+) = x(t^-)$ and there is no state jump in (12), even at times where $\lambda_d(\cdot)$ may be discontinuous. Therefore, solutions to (12), if they exist, must be continuous. Let us deal with the well-posedness of (13), equivalently of (14). Assume that $u_d(\cdot)$ is AC, then posing $z = \dot{q}_d - u_{2d}$ it follows that $\dot{z} = \ddot{q}_d - \dot{u}_{2d}$ where $\dot{u}_{2d}$ is the almost-everywhere derivative of $u_{2d}$. This gives the DI:

$$ \dot{z}(t) \in -\dot{u}_{2d}(t) + \frac{u_{1d}(t)}{m} - \mu g \operatorname{sgn}(z(t)). \quad (18) $$

Using for instance [96], it is deduced that the DI in (18) has a unique Lipschitz continuous solution on $[0, T)$ for any $T > 0$, for any $z(0) = z_0 \in \mathbb{R}$ and any $u_d(\cdot)$ with $\dot{u}_{1d}(\cdot)$ and $\dot{u}_{2d}(\cdot)$ bounded in the extended (or local) $L_2$ norm. Thus, $\dot{q}_d(\cdot)$ is AC (Lipschitzness implies AC), $\ddot{q}_d(\cdot)$ is defined almost everywhere, and $q_d(t) = q_d(0) + \int_0^t \dot{q}_d(s)ds$ is continuously differentiable. However, the multiplier $\lambda_d(\cdot)$ can undergo jumps at some junction times. Therefore, the
closed-loop plant (12) is an LCS with exogenous signal $F_{cl}(t) \triangleq -\begin{pmatrix} k_{21} & k_{22} \\ 0 & 0 \end{pmatrix} x_d(t) - \begin{pmatrix} g_{21} & 0 \\ 0 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} u_{2d}(t) \\ 1 \end{pmatrix}$ entering the complementarity constraints, where no state jump can occur.

**Remark 3.6.** Let us consider the DI in (8). Let us explain how a feedback using the variable $\lambda_1$ (equivalently $\lambda_t$) may be implemented inside the signum set-valued function, without direct measurement of the tangent contact force (which is unrealistic from a practical point of view), but only from position and velocity measurements. Let $u_2 = u_2(q, \dot{q}, \lambda_t)$, then $\lambda_t \in \text{sgn}(\dot{q} - u_2(q, \dot{q}, \lambda_t)) \Leftrightarrow \dot{q} - u_2(q, \dot{q}, \lambda_t) \in N_{[-1,1]}(\lambda_t)$, where the inversion of maximal monotone set-valued mapping is used. This is a generalized equation with unknown $\lambda_t$, whose well-posedness can be characterized. Indeed, $u_2(q, \dot{q}, \lambda_t) = k_{21}(q - q_d) + k_{22}(\dot{q} - \dot{q}_d) + g_{21}(\lambda_1 - \lambda_{1d}) + u_{2d}$, with $\lambda_1(q, \dot{q}) = \frac{\lambda_1(q, \dot{q})+1}{2}$. This yields when $g_{21} > 0$:

$$ \begin{gathered} -\dot{q} + k_{21}(q - q_d) + k_{22}(\dot{q} - \dot{q}_d) + g_{21}(\frac{\lambda_1+1}{2} - \lambda_{1d}) + u_{2d} \in N_{[-1,1]}(\lambda_t) \\ \lambda_t(q, \dot{q}) = \text{proj}([-1, 1]; 2\frac{\dot{q}-k_{21}(q-q_d)-k_{22}(\dot{q}-\dot{q}_d)+g_{21}\lambda_{1d}-\frac{g_{21}}{2}-u_{2d}}{g_{21}}), \end{gathered} \tag{19} $$

see Fact 1. Therefore, such feedback regularizes the set-valued signum function, as is known in other contexts [82, Remark 2.8] [97]. It is also possible to prove that the LCP in (12) (b) has a unique solution $\lambda(t)$ for any $C_{cl}x(t) + F_{cl}(t)$ and $g_{21} > 0$. In spite of the fact that $D_{cl} \supseteq 0$, the result holds due to the special structure of the LCP (especially the second line) and rank($D_{cl}$) = 2. Also, the result concerns only $\lambda_1 = \frac{\lambda_1+1}{2}$, not $\lambda_2$. It is noteworthy that the multiplier $\lambda_t$ in (19) depends only on known control parameters, mass state, and desired signals. This is expected since the LCP in (12) does not involve mechanical parameters. Taking $g_{21} > 0$ means that the multiplier $\lambda_1$, solution to the LCP in (12) (b), can be explicitly calculated. So there is no need for a contact force measurement in this case, measurement of $x$ is sufficient to solve the LCP in (12) (b). Note that $\lambda_1$ does not depend on $\mu$ nor on the normal contact force, and can therefore be computed even if these variables are subjected to uncertainties.

Let us continue the well-posedness study. When $g_{21} = 0$, it follows from [98, Theorems 11 and 26] that the closed-loop system (12) has unique AC solutions. This can be shown using passivity of the quadruple $(A_{cl}, B_{cl}, C_{cl}, D_{cl})$ in Assumption 3.5, the fact that $\text{Im}(D_{cl} + N_{R_+^2}) - F_{cl}(t) = \text{Im}(D_{cl}) + R_-^2 - F_{cl}(t) = R^2$ (here a result from [99] is used with the fact that $N_{R_+^2}(\cdot) = N_{R_-^2}^{-1}(\cdot)$). This together with the fact that $F_{cl}(\cdot)$ is AC allows us to guarantee that [98, Theorem 26 (ii)] holds true. When $g_{21} > 0$ (which will be imposed in some cases), then it follows from (19) that $\lambda_1(\cdot)$ (which is the only component of $\lambda(\cdot)$ that intervenes in (12) (a)) is a Lipschitz continuous function of the state and exogenous terms. Therefore (12) (a) boils down to an ordinary differential equation with Lipschitz continuous right-hand side, existence and uniqueness of locally AC solutions follows from general results [100, Proposition C. 3.8].

The conclusion of this section is that imposing mild conditions on $u_d(\cdot)$ together with passivity as in Assumption 3.5, guarantees the well-posedness of the closed-loop, the desired, and the error dynamics, with AC states $x(\cdot)$ and $x_d(\cdot)$ (hence the error dynamics (16) is well-posed also).

### 3.3. Equilibrium Points of the Error System

Let us assume that Assumption 3.5 holds. As pointed out above this implies that $g_{21} \ge 0$. The equilibria of the LCS (15) are characterized as follows:

**Proposition 3.7.** Consider the error dynamics in (15) or (16). Let Assumption 3.5 hold. Then, the unique equilibrium point is $e^* = (0, 0)^\top$ and the multipliers at equilibrium are time-varying such that $\lambda^*(t) = (\lambda_{1d}(t), |x_{2d}(t) - u_{2d}(t)|)^\top$.

The proof is in Appendix D.

## 4. Analysis of the Desired Trajectories Generator

### 4.1. Periodic Trajectories in Frictional Oscillators: a Short Review

In practice, it is certainly interesting to understand which desired trajectories can be designed, in spite of the fact that this is not the main topic of this article. Let us consider the case when $u_{2d} \ne 0$ which is studied in [2, 3, 4, 12] when $u_{2d}$ is constant [15] and when $u_{2d}(t)$ is time-varying. The stick-slip mechanical system considered in [2, 3, 4] is
a simple mass-spring system on a moving belt with constant speed. However, the system presented in Figure 1a differs because it includes an external force $u_{1d}(t)$ acting on the mass and allows for time-varying belt speed $u_{2d}(t)$. In [2], the authors introduce the shooting method as a periodic-solution-finder to address the computational challenges of solving stiff differential equations. To tackle this issue, the authors present the switch model as an alternate friction model to simulate stick-slip vibrations by solving a set of non-stiff differential equations. Furthermore, a time-dependent friction model is studied. The authors in [4] study bifurcations of periodic solutions in mechanical systems with dry friction and discuss the existence of infinitely unstable periodic solutions through repulsion sliding mode in stick/slip mechanical systems.

The paper [101] proposes two numerical methods for the computation of periodic solutions in dynamical systems with maximal monotone right-hand side, including LCS and systems with dry friction. It provides formal justifications to ensure the consistency of the proposed numerical schemes. The results in [101] could be applied to the frictional oscillator we work with, which is represented as an LCS in (13) and, equivalently, as a differential inclusion with maximal monotonicity properties in (14). A recent paper [102] introduces a new method for computing periodic solutions of systems with frictional occurrences. The approach formulates all equations, including friction, as equalities and allows for an exact Coulomb’s friction law without regularizing the friction force. This equality-based approach is implemented through a highly compact weighted residual formulation and it represents accurately multiple sticking and sliding phases. Given that Coulomb’s friction is represented within the complementarity framework in (10), the method developed in [102] is applicable to design periodic desired trajectories. However, this is not investigated further in this article.

## 4.2. Numerical Example of Periodic Desired Trajectory

This simulation involves implementing a suitable desired input $u_d$ and observe the resulting desired trajectory. The simulation is done with the INRIA software package siconos<fnref n="1" />[103, 104] dedicated to the simulation of this class of nonsmooth systems. In order to design a periodic desired trajectory $x_d$, the frictional oscillator is excited by periodic desired input $u_d = (u_{1d}, u_{2d})^\top$ where $u_{1d}$ and $u_{2d}$ have the same period $T$ but different amplitudes, as shown in Figure 3. Take $x_d(0) = (1, 0)^\top$ and the time step $h = 0.01$s.

<figure><img src="image_6.png" /><figcaption>Figure 3: Numerical simulation of (13) with $u_{1d} = 0.8\mu g \cos 2t$, $u_{2d} = 0.5\mu g \cos 2t$</figcaption></figure>

The simulation shows a periodic desired trajectory with period $T = \frac{2\pi}{2}$ s which is the period of $u_d$. Interestingly enough the frictional oscillator exhibits both sticking and slipping behaviours, as demonstrated by the relative velocity $V_{rel} = x_{2d} - u_{2d}$ plot.

<footnote n="1">¹https://nonsmooth.gricad-pages.univ-grenoble-alpes.fr/siconos/index.html</footnote>
### 4.3. Arbitrary Desired Trajectories

Generating desired trajectories using (13) has the advantage that $q_d(\cdot)$ is a "natural" trajectory of the system when subjected to an external excitation (like a harmonic signal, widely used in the nonlinear dynamics literature). It has the advantage of generating trajectories with an arbitrary number of stick/slip transitions (even chaotic oscillator systems could be tracked), and the drawback of restricting the class of tracked trajectories. Instead of constructing (13), we may define desired trajectories $q_d(\cdot)$, such that $\dot{q}_d(\cdot)$ and $\ddot{q}_d(\cdot)$ are defined almost everywhere. Let $u_1 = m\ddot{q}_d + k_{11}\tilde{q} + k_{12}\tilde{\dot{q}}_d$. Then from (8) we obtain:

$$m\ddot{\tilde{q}} - k_{12}\ddot{\tilde{q}} - k_{11}\tilde{q} \in -\mu mg \operatorname{sgn}(\dot{q} - u_2). \quad (20)$$

Let $u_2 = \dot{q}_d + k_{21}\tilde{q} + k_{22}\ddot{\tilde{q}}$. Then we get $m\ddot{\tilde{q}} - k_{12}\ddot{\tilde{q}} - k_{11}\tilde{q} \in -\mu mg \operatorname{sgn}(-k_{21}\tilde{q} + (1 - k_{22})\ddot{\tilde{q}})$. This is rewritten equivalently as $\dot{x} = A'_{cl}x + B'_{cl}\lambda$, with $\lambda \in \operatorname{sgn}(C'_{cl}x)$, $A'_{cl} = A_{cl} = \begin{pmatrix} 0 & 1 \\ \frac{k_{11}}{m} & \frac{k_{12}}{m} \end{pmatrix}$, $B'_{cl} = \begin{pmatrix} 0 \\ -\mu g \end{pmatrix}$, $C'_{cl} = (-k_{21} \quad 1 - k_{22})$. Stability analysis relies on the SSPMI $(A'_{cl}, B'_{cl}, C'_{cl}, \epsilon)$. Then the analysis can be led in a similar way as in section 5.1, using the fact that $x^T C'_{cl} \operatorname{Sgn}(C'_{cl}x) \ge 0$ since $(0,0)$ belongs to the graph of $\operatorname{Sgn}(\cdot)$. Notice, however, that the use of a generator as in (13) and incremental passivity (see Remark 5.2), allow us to consider set-valued nonlinearities which do not satisfy this property. Notice that it is quite possible to generate trajectories $q_d(\cdot)$ from (13), then use them as in (20), hence combining the advantages of both approaches. An important point is that the subsequent stability analyses and controller parameters calculations, are quite similar in both cases (see Remark 5.2).

## 5. Trajectory Tracking of Frictional Oscillator in Nominal Case

In this section, the trajectory tracking of the frictional oscillator is solved without uncertainties on the friction coefficient. Numerical simulations are presented to illustrate the theoretical result.

### 5.1. Closed-loop Error Stability Analysis

The stability analysis led in [83], based on the LCS formulations for the nominal case is applied in this section. Assume that all trajectories are AC (see section 3.2). Then, the following results hold. The proof of the first part mimics that of [83, Proposition 3.8] and is recalled for the sake of completeness. The second part is original.

**Proposition 5.1.** Consider the error dynamics (16) with bounded initial conditions, and let Assumption 3.5 hold, with constant $\epsilon > 0$. Then: (i) $e(\cdot)$ converges globally exponentially fast to zero; (ii) When $g_{21} > 0$, the multiplier error $(\lambda_1 - \lambda_{1d})(\cdot)$ converges globally exponentially fast to zero also.

*Proof.* (i) Consider the Lyapunov candidate function $V(e(t)) = e^T P e$, with $P = P^T > 0$ the solution of the closed-loop passivity matrix inequality. Then, along the closed-loop error dynamics trajectories:

$$\dot{V}(e(t)) = e^T (A_{cl}^T P + P A_{cl}) e(t) + 2e(t)^T P B_{cl} (\lambda(t) - \lambda_d(t))$$

Starting from (10) and (13), the following is obtained:

$$\left\{ \begin{array}{l} w(t) = C_{cl}x(t) + D_{cl}\lambda(t) + Cx_d(t) + Fu_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} - C_{cl}x_d(t) \\ w_d(t) = C_{cl}x_d(t) + D_{cl}\lambda_d(t) + Cx_d(t) + Fu_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} - C_{cl}x_d(t). \end{array} \right.$$

Therefore, it is inferred that:

$$\left\{ \begin{array}{l} \lambda(t) \in \mathcal{N}_{S(t)}(C_{cl}x(t) + D_{cl}\lambda(t)) \\ \lambda_d(t) \in \mathcal{N}_{S(t)}(C_{cl}x_d(t) + D_{cl}\lambda_d(t)), \end{array} \right. \quad (21)$$

with $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\}$. From the monotonicity of the normal cone, it follows that:

$$(w(t) - w_d(t))^T (\lambda(t) - \lambda_d(t)) \le 0 \quad (22)$$
Given that $w(t) - w_d(t) = C_{cl}e(t) + D_{cl}(\lambda(t) - \lambda_d(t))$, $\dot{V}(e(t))$ is written in matrix form as:

$$ 
\begin{aligned}
\dot{V}(t) &= \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top M \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top \begin{pmatrix} 0 & (C + FK)^\top \\ C + FK & D + FG + (D + FG)^\top \end{pmatrix} \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} \\
&= \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top M \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + 2\Delta\lambda^\top\Delta w \le \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top M \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} \le -\epsilon e(t)^\top P e(t),
\end{aligned} 
\quad (23) $$

where $M$ is the matrix defined in (17) and satisfying the strict state passivity LMI. Hence, $e^* = 0$ is globally exponentially stable. (ii) This follows from the projection form of $\lambda_t$ in (19), i.e.:

$$ \lambda_t(q, \dot{q}) = \operatorname{proj} \left( [-1, 1]; \frac{2}{g_{21}} \left( (\dot{q} - \dot{q}_d) - k_{21}(q - q_d) - k_{22}(\dot{q} - \dot{q}_d) + g_{21}\lambda_{1d} - \frac{g_{21}}{2} + \dot{q}_d - u_{2d} \right) \right) $$

Similarly, $\lambda_{t,d}$ is written as follows

$$ \lambda_{t,d} \in \operatorname{sgn}(\dot{q}_d - u_{2d}) = \operatorname{sgn} \left( \frac{2}{g_{21}} (\dot{q}_d - u_{2d}) \right) \Leftrightarrow \lambda_{t,d} = \operatorname{proj} \left( [-1, 1]; \frac{2}{g_{21}} (\dot{q}_d - u_{2d}) + \lambda_{t,d} \right) $$

Recall that $x \triangleq (q, \dot{q})^\top$ and $x_d \triangleq (q_d, \dot{q}_d)^\top$. Given that $\lambda_t = 2\lambda_1 - 1$ and $\lambda_{t,d} = 2\lambda_{1,d} - 1$, then:

$$ 
\begin{aligned}
|\lambda_1 - \lambda_{1d}| &= \frac{1}{2} |\lambda_t - \lambda_{t,d}| \\
&= \frac{1}{2} \left| \operatorname{proj}([-1, 1]; \frac{2}{g_{21}} (-[C_{cl}]_{1,\bullet}e + \lambda_{t,d} + x_{2d} - u_{2d})) - \operatorname{proj}([-1, 1]; \frac{2}{g_{21}} (\lambda_{t,d} + x_{2d} - u_{2d})) \right| \\
&\le \frac{1}{g_{21}} \|C_{cl}e\|,
\end{aligned} 
$$

where we have used the Lipschitz continuity property of the projection operator ([105, Proposition 4.8]) to obtain the last inequality. Hence, $\lim_{t \to \infty} |\lambda_1(t) - \lambda_{1d}(t)| \le \frac{1}{g_{21}} \lim_{t \to \infty} \|C_{cl}e(t)\| = 0$. $\square$

**Remark 5.2.** It is noteworthy that inequality (22) involves some incremental passivity of the closed-loop system, with $\Delta\lambda$ and $\Delta w$ as the incremental input/output pair, and $\Delta\lambda^\top\Delta w$ the supply rate. In this framework, (23) is the infinitesimal dissipation inequality. Notice also that both SSPMI($A_{cl}$, $B_{cl}$, $C_{cl}$, $D_{cl}$, $\epsilon$) and SSPMI($A'_{cl}$, $B'_{cl}$, $C'_{cl}$, $\epsilon$) imply that there exists a solution $P$ such that $PB_{cl} = C_{cl}^\top$ or $PB'_{cl} = C'_{cl}^\top$ with $P \ne \operatorname{diag}(p_{11}, p_{12})$, because otherwise $PA_{cl} + A_{cl}^\top P$ cannot be negative definite. It is inferred that both approaches for the desired trajectories generation (sections 4.1 and 4.3), yield similar calculations.

## 5.2. Numerical Simulations

Take $m = 1$ kg, $g = 9.8m/s^2$, and $\mu = 0.5$. Let us take $G = 0$ and check if there exist $K$ and $P$ such that the closed-loop system's quadruple ($A_{cl} = A + EK$, $B_{cl} = B$, $C_{cl} = C + FK$, $D_{cl} = D$) is strictly state passive, which means equivalently solving the SSPMI in (17). The SSPMI is transformed into an LMI (see Appendix A.1) and solved using MOSEK 9.3.14 solver [106], giving the solution:

$$ K = \begin{pmatrix} -1.49 & -1.17 \\ -5.46 & -12.1 \end{pmatrix} \quad \text{and} \quad P = \begin{pmatrix} 1.34 & 0.557 \\ 0.557 & 1.34 \end{pmatrix}, \quad (24) $$

with $\epsilon = 0.01$. Notice that when $G = 0$ in (11), the passivity LMI implies that $PB_{cl} = C_{cl}^\top$ since $D + D^\top = 0$ when $g_{21} = 0$, which yields $-2p_{12}\mu g = k_{21}$ and $-2p_{22}\mu g = k_{22} - 1 < 0$. Figure 4 shows the numerical simulation of the error dynamical system in (16) with external input $u = (u_1, u_2)^\top = K(x - x_d) + u_d$. Take $x(0) = (0, -1)^\top$, $x_d(0) = (1, 0)^\top$ and the time step $h = 0.01$.
<figure><img src="image_7.png" /><figcaption>Figure 4: Simulation of the error dynamics (16), with $u_{1d} = u_{2d} = 0.5m\mu g \cos 2t$ and $G = 0$.</figcaption></figure>

The numerical simulation demonstrates that tracking of the desired trajectory $x_d$ is achieved without the requirement of an additional feedback from $\lambda_1$ (i.e., $G = 0$). Additionally, it is interesting to observe that the complementarity variables of the closed-loop system ($\lambda$ and $w$) converge to the complementarity variables of the desired system ($\lambda_d$ and $w_d$). The controller $u$ shows a periodic behaviour in the steady state response (i.e., $e = x - x_d = 0$), where $u$ converges to the periodic desired input $u_d$.

## 6. Parametric Robustness Analysis

This section is dedicated to analyse the robustness of the above scheme, which is designed assuming the strict state passivity of the closed-loop system. As shown in [83], strong passivity of $(A_{cl}, B_{cl}, C_{cl}D_{cl})$ in (12) ensures robustness against a large class of uncertainties. However, strong passivity of this quadruple cannot be achieved for (12). The analysis focuses on the case when uncertainties $\Delta\mu$ are present in the friction coefficient $\mu$, which is supposed to be known only with a nominal value $\mu_0 > 0$. Then, it is crucial that the multiplier $\lambda_1$ is used in the controller, with $g_{21} > 0$. The analysis in this section relies on theoretical results from [83], specifically Propositions B.1, B.3 and Lemma B.2 which are presented in Appendix B.

### 6.1. Robustness analysis with bounded and time-varying uncertainties $\Delta\mu(\dot{q}_{rel}) = \mu(\dot{q}_{rel}) - \mu_0$

Let us recall that $\dot{q}_{rel} = \dot{q} - u_2$. Here, we consider that the nominal friction coefficient $\mu_0$ (used in the desired dynamics and the nominal model of the plant) and the "real" plant's friction coefficient $\mu$ are not the same. It could even be that $\mu_0$ is chosen constant while $\mu = \mu(t, \dot{q})$, then $\Delta\mu(\dot{q}_{rel}) = \mu(\dot{q}_{rel}) - \mu_0$. From (10) the uncertainties satisfy:

$$ \Delta A = 0, \Delta C = 0, \Delta D = 0, \Delta E = 0, \Delta F = 0, \Delta B = \begin{pmatrix} 0 & 0 \\ -2g\Delta\mu & 0 \end{pmatrix}, \Delta E' = \begin{pmatrix} 0 \\ g\Delta\mu \end{pmatrix}. \quad (25) $$

Compared to [83], the uncertainties are much simpler to handle. This will allow us to analyse the robustness in detail. The closed-loop system in (12) cannot be made strongly passive (see Remark 3.3). Thus, it is not possible to apply directly the results of robustness analysis for LCS from [83]. Instead, let us relax the conditions of strong passivity by applying the results derived in [83, section 4.3] (see Appendix B) to address the robustness analysis in frictional oscillators. The desired system dynamics in (13) are written as:

$$ \begin{cases} \dot{x}_d(t) = A_0x_d(t) + B_0\lambda_d(t) + E_0u_d(t) + E'_0 \\ 0 \le \lambda_d(t) \perp w_d(t) = C_0x_d(t) + D_0\lambda_d(t) + F_0u_d(t) + F'_0 \ge 0 \end{cases} $$
where $A_0 = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$, $B_0 = \begin{pmatrix} 0 & 0 \\ -2\mu_0 g & 0 \end{pmatrix}$, $C_0 = \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix}$, $D_0 = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$, $E_0 = \begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix}$, $F_0 = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$, $E'_0 = \begin{pmatrix} 0 \\ \mu_0 g \end{pmatrix}$ and $F'_0 = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$. The dynamics of the closed-loop system in (12) with parametric uncertainties are written as:

$$ \begin{cases} \dot{x}(t) = (A_0 + E_0 K_0)x(t) + (B_0 + \Delta B)\lambda(t) - E_0 K_0 x_d(t) + E_0 u_d(t) + E'_0 + \Delta E' \\ 0 \le \lambda(t) \perp w(t) = (C_0 + F_0 K_0)x(t) + (D_0 + F_0 G_0)\lambda(t) - F_0 K_0 x_d(t) - F_0 G_0 \lambda_d(t) + F_0 u_d(t) + F'_0 \ge 0 \end{cases} \quad (26) $$

where $K_0$ and $G_0$ are computed solving the nominal system matrix inequality as in (17), see the proof of Proposition 6.1. The closed-loop error dynamics ($e = x - x_d$) are:

$$ \left\{ \begin{aligned} \dot{e}(t) &= A_{cl}e(t) + B_{cl}\Delta\lambda(t) + \Delta B\lambda(t) + \Delta E' \\ \Delta w(t) &= w(t) - w_d(t) = C_{cl}e(t) + D_{cl}\Delta\lambda(t) \\ 0 &\le \lambda(t) \perp w(t) \ge 0 \quad \text{and} \quad 0 \le \lambda_d(t) \perp w_d(t) \ge 0 \end{aligned} \right. \quad (27) $$

where $\Delta\lambda(t) = \lambda(t) - \lambda_d(t)$ and $A_{cl} = A_0 + E_0 K_0$, $B_{cl} = B_0$, $C_{cl} = C_0 + F_0 K_0$ and $D_{cl} = D_0 + F_0 G_0$. The stability result is stated as follows.

**Proposition 6.1.** Consider the error dynamics in (27). Assume that:

* (i) The uncertainties are bounded such that $\Delta B^\top \hat{\Lambda}_B^{-1} \Delta B \preceq \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}$ and $\Delta E'^\top \hat{\Lambda}_1^{-1} \Delta E' \preceq I_m$, given the matrices $\hat{\Lambda}_k = \hat{\Lambda}_k^\top > 0$ where $k \in \{1, B\}$.

* (ii) There exists a matrix $P_0 = P_0^\top > 0$ such that the quadruple $(A_{cl}, B_0, C_{cl}, D_{cl})$, corresponding to the nominal closed-loop system in (12) with $\mu = \mu_0$, is strictly state passive with $\epsilon > 0$.

* (iii) The control gain $G = \begin{pmatrix} 0 & 0 \\ g_{21} & 0 \end{pmatrix}$ is chosen such that $g_{21} > \frac{1}{2}$.

* (iv) The following inequality holds: $\sigma_{\max}(-P_0(\hat{\Lambda}_1 + \hat{\Lambda}_B)P_0) + \sigma_{\max}((-P_0B_0 + C_{cl}^\top)(-B_0^\top P_0 + C_{cl})) \le \epsilon \lambda_{\max}(P_0)$ for $P_0 = P_0^\top > 0$ and $\epsilon > 0$.

Then, the solution of the error dynamics in (27) is globally uniformly ultimately bounded (GUUB), satisfying: $\|e\| \le \sqrt{\frac{\lambda_{\max}(P_0)\lambda_{\min}^{-1}(\hat{\Lambda}_1)\rho}{\lambda_{\min}^2(P_0)e'}}$, where $\rho = \frac{\det(\hat{\Lambda}_B)}{\hat{\Lambda}_{B11}}$ defines the bound on $\Delta E$.

For the sake of notation simplicity, the matrices $M_0$ and $\Delta M_0$ are not explicitly stated in the theorem. These matrices are defined in the following proof.

*Proof.* The proof is divided into two main parts. In the first part, a suitable storage function is derived for the error dynamics in (27) to serve as a Lyapunov function. Then, the second part presents explicit conditions to prove that the Lyapunov function decreases.

*Construction of Lyapunov function.* The storage function $V(t) = e^\top P_0 e$ is associated with the SSPMI$(A_{cl}, B_{cl}, C_{cl}, D_{cl}, \epsilon)$ transformed into LMI (see Appendix A.1) given by: $M_0 \preceq \begin{pmatrix} -\epsilon P_0 & 0 \\ 0 & 0 \end{pmatrix}$ where $M_0$ is defined as:

$$ M_0 \triangleq \begin{pmatrix} A_{cl}^\top P_0 + P_0 A_{cl} & P_0 B_{cl} - C_{cl}^\top \\ B_{cl}^\top P_0 - C_{cl} & -D_{cl} - D_{cl}^\top \end{pmatrix} \quad (28) $$

Let us define $p(t, \lambda) \triangleq \Delta B \lambda(t) + \Delta E'$. The rate of change of the storage function is:

$$ \dot{V} = e^\top(A_{cl}^\top P_0 + P_0 A_{cl})e + 2e^\top P_0 B_{cl} \Delta \lambda + 2e^\top P_0 p(t, \lambda) \pm 2 \Delta \lambda^\top \Delta w \quad (29) $$
The added and subtracted term is to recover the matrix $M_0$ in (28). In matrix form, it is represented as:

$$2\Delta\lambda^\top\Delta w = \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top \begin{pmatrix} 0 & C_{cl}^\top \\ C_{cl} & D_{cl} + D_{cl}^\top \end{pmatrix} \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}$$

Then,

$$\dot{V} = \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top M_0 \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + 2e^\top P_0 p(t, \lambda) + 2\Delta\lambda^\top \Delta w \quad (30)$$

Let us substitute the expression of $p(t, \lambda)$ in (30), then:

$$\dot{V} = \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top M_0 \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top \begin{pmatrix} 0 & P_0 \Delta B \\ \Delta B^\top P_0 & 0 \end{pmatrix} \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + 2e^\top P_0 \Delta E' + 2\Delta\lambda^\top \Delta w$$

By considering monotonicity of the mapping $-\lambda \mapsto w$ (or $-\lambda_d \mapsto w_d$) in (27), the following inequality is obtained:

$$\dot{V} \leq - \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top (-M_0 + \Delta M_0) \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + 2e^\top P_0 \Delta E'$$

with $\Delta M_0 \triangleq \begin{pmatrix} 0 & -P_0 \Delta B \\ -\Delta B^\top P_0 & 0 \end{pmatrix}$. For any $\hat{\Lambda}_1^\top = \hat{\Lambda}_1 > 0$, the following holds:

$$|2e^\top P_0 \Delta E'| \leq e^\top P_0 \hat{\Lambda}_1 P_0 e + \Delta E'^\top \hat{\Lambda}_1^{-1} \Delta E'$$

Thus,

$$\dot{V} \leq - \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top \left( -M_0 + \Delta M_0 - \begin{pmatrix} P_0 \hat{\Lambda}_1 P_0 & 0 \\ 0 & 0 \end{pmatrix} \right) \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix} + \Delta E'^\top \hat{\Lambda}_1^{-1} \Delta E' \quad (31)$$

The expression $\begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}^\top \Delta M_0 \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}$ is equivalently written as: $-2e^\top P_0 \Delta B \Delta \lambda$. Take $\hat{\Lambda}_B = \begin{pmatrix} d & e \\ e & f \end{pmatrix} > 0$ for some $d, e, f \in \mathbb{R}$. Let us assume that:

$$\Delta B^\top \hat{\Lambda}_B^{-1} \Delta B = \frac{1}{df - e^2} \begin{pmatrix} 4g^2 \Delta \mu^2 d & 0 \\ 0 & 0 \end{pmatrix} \asymp \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix} \quad (32)$$

Let $z \triangleq \begin{pmatrix} e \\ \Delta\lambda \end{pmatrix}$, then

$$z^\top \left[ -M_0 + \Delta M_0 - \begin{pmatrix} P_0 \hat{\Lambda}_1 P_0 & 0 \\ 0 & 0 \end{pmatrix} \right] z \geq z^\top \begin{pmatrix} -(M_0)_{11} - P_0(\hat{\Lambda}_1 + \hat{\Lambda}_B)P_0 & -(M_0)_{12} \\ -(M_0)_{21} & -(M_0)_{22} - \hat{I} \end{pmatrix} z$$

where $\hat{I} \triangleq \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}$. Let us define:

$$\begin{align}
\hat{L} &\triangleq -M_0 + \Delta M_0 - \begin{pmatrix} P_0 \hat{\Lambda}_1 P_0 & 0 \\ 0 & 0 \end{pmatrix} - \begin{pmatrix} \epsilon' P_0 & 0 \\ 0 & 0 \end{pmatrix} \nonumber \\
&\triangleq \begin{pmatrix} -A_{cl}^\top P_0 - P_0 A_{cl} - P_0(\hat{\Lambda}_1 + \hat{\Lambda}_B)P_0 - \epsilon' P_0 & -P_0 B_{cl} + C_{cl}^\top \\ -B_{cl}^\top P_0 + C_{cl} & D_{cl} + D_{cl}^\top - \hat{I} \end{pmatrix} \geq 0 \tag{33}
\end{align}$$

The term $\epsilon' P_0$ is introduced for bounding the error as will be shown later. For stability analysis, it is required to prove that $\dot{V} \leq 0$ holds in a certain domain to be defined. This means that it is necessary to show that the matrix inequality $\hat{L} \geq 0$ holds for $\epsilon' > 0$, as detailed below.
Conditions for $\hat{L} \ge 0$ in (33) to hold true. The conditions presented below are derived from proving positive semi-definiteness, relaxing strong passivity into strict state passivity (see Appendix B). Take

$$ \hat{Q} \triangleq \underbrace{-A_{cl}^{\top}P_{0} - P_{0}A_{cl} - \epsilon'P_{0}}_{\hat{Q}_{0}} + \underbrace{(-P_{0}(\hat{\Lambda}_{1} + \hat{\Lambda}_{B})P_{0})}_{\Delta\hat{Q}}, $$

$\hat{S} \triangleq \underbrace{-P_{0}B_{cl} + C_{cl}^{\top}}_{\hat{S}_{0}}$ and $\hat{R} \triangleq \underbrace{D_{cl} + D_{cl}^{\top}}_{\hat{R}_{0}} + \underbrace{(-\hat{I})}_{\Delta\hat{R}}$, where $P_{0} = \begin{pmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{pmatrix}$ such that $P_{0} = P_{0}^{\top} > 0$, is the solution of $M_{0} \preceq \begin{pmatrix} -\epsilon P_{0} & 0 \\ 0 & 0 \end{pmatrix}$ with $M_{0}$ defined in (28). Recall that $\hat{\Lambda}_{1} = \begin{pmatrix} a & b \\ b & c \end{pmatrix} > 0$ for some $a, b, c \in \mathbb{R}$. The matrices $A_{cl}, B_{cl}, C_{cl}$ and $D_{cl}$ are defined in (12). Let us explicitly write the matrices $\hat{Q} = \hat{Q}_{0} + \Delta\hat{Q}, \hat{R} = \hat{R}_{0} + \Delta\hat{R}$ and $\hat{S} = \hat{S}_{0}$:

$$ \hat{Q}_{0} = \begin{pmatrix} -2p_{12} \frac{k_{11}}{m} - \epsilon'p_{11} & -p_{11} - \frac{k_{12}}{m}p_{12} - \frac{k_{11}}{m}p_{22} - \epsilon'p_{12} \\ -p_{11} - \frac{k_{12}}{m}p_{12} - \frac{k_{11}}{m}p_{22} - \epsilon'p_{12} & -2p_{12} - 2 \frac{k_{12}}{m}p_{22} - \epsilon'p_{22} \end{pmatrix}, $$

$$ -\Delta\hat{Q} = \begin{pmatrix} p_{11}^{2}(a+d) + p_{12}^{2}(c+f) + 2p_{11}p_{12}(b+e) & (p_{12}^{2} + p_{11}p_{22})(b+e) + p_{11}p_{12}(a+d) + p_{12}p_{22}(c+f) \\ (p_{12}^{2} + p_{11}p_{22})(b+e) + p_{11}p_{12}(a+d) + p_{12}p_{22}(c+f) & p_{12}^{2}(a+d) + p_{22}^{2}(c+f) + 2p_{12}p_{22}(b+e) \end{pmatrix}, $$

$$ \hat{S} = \hat{S}_{0} = \begin{pmatrix} 2p_{12}\mu g + k_{21} & 0 \\ 2p_{22}\mu g - 1 + k_{22} & 0 \end{pmatrix}, \quad \hat{R}_{0} = \begin{pmatrix} 2g_{21} & 0 \\ 0 & 0 \end{pmatrix}, \quad \text{and} \quad \Delta\hat{R} = \begin{pmatrix} -1 & 0 \\ 0 & 0 \end{pmatrix}, $$

where the constants $a, b, c, d, e$ and $f$ are all defined above. By considering the conditions for positive semi-definiteness in Definition A.1, the matrix $\hat{L} \triangleq \begin{pmatrix} \hat{Q}_{0} + \Delta\hat{Q} & \hat{S}_{0} \\ \hat{S}_{0}^{\top} & \hat{R}_{0} + \Delta\hat{R} \end{pmatrix} \ge 0$ if and only if:

1. $\hat{R} \ge 0 \Leftrightarrow \begin{pmatrix} 2g_{21} - 1 & 0 \\ 0 & 0 \end{pmatrix} \ge 0$ which is satisfied for all $g_{21}$ such that $g_{21} \ge \frac{1}{2}$.

2. Im($\hat{S}^{\top}$) $\subseteq$ Im($\hat{R}$) and this condition holds as a consequence of the following:

$$ \begin{aligned} \operatorname{Im}(\hat{S}_{0}^{\top}) &= \operatorname{Im}(\hat{S}^{\top}) = \left\{ w \in \mathbb{R}^{2} \mid w = \begin{pmatrix} (2\mu g p_{12} + k_{12})x + (2\mu g p_{22} - 1 + k_{22})y \\ 0 \end{pmatrix} \text{ for some } x, y \in \mathbb{R} \right\} \\ &\subseteq \operatorname{Im}(\hat{R}) = \left\{ w \in \mathbb{R}^{2} \mid w = \begin{pmatrix} (2g_{21} - 1)t \\ 0 \end{pmatrix} \text{ for some } t \in \mathbb{R} \right\} \end{aligned} $$

3. $\hat{Q} \ge \hat{S}\hat{R}^{\dagger}\hat{S}^{\top}$. Let us follow the conditions presented in Proposition B.3 in order to check the validity of this item (i.e., $Q \ge SR^{\dagger}S^{\top}$), then: (i) $\hat{R} \ge 0$ for $g_{21} \ge \frac{1}{2}$, (ii) $\operatorname{Im}(\Delta\hat{R}) = \left\{ w \in \mathbb{R}^{2} \mid w = \begin{pmatrix} -y \\ 0 \end{pmatrix} \text{ for some } y \in \mathbb{R} \right\} \subseteq \operatorname{Im}(R_{0}) = \left\{ w \in \mathbb{R}^{2} \mid w = \begin{pmatrix} 2g_{21}x \\ 0 \end{pmatrix} \text{ for some } x \in \mathbb{R} \right\}$ for any $x \in \mathbb{R}$, (iii) $\operatorname{rank}(R_{0}) = 1$, $\sigma_{\max}(\Delta\hat{R}) = 1 < \sigma_{1}(R_{0}) = 2g_{21}$ which holds for $g_{21} > \frac{1}{2}$, (iv) $\sigma_{\max}(\operatorname{OT}) + \sigma_{\max}(\Delta\hat{Q}) < \epsilon\lambda_{\max}(P_{0})$ where $\operatorname{OT}$ is defined in (B.3) and in this case, where $\Delta\hat{S} = 0$, it is written as:

$$ \operatorname{OT} = -\hat{S}_{0}\hat{R}_{0}^{\dagger}\Delta\hat{R}\hat{R}_{0}^{\dagger}\hat{S}_{0}^{\top} $$

with $\hat{R}_{0}^{\dagger} = \begin{pmatrix} \frac{1}{2g_{21}} & 0 \\ 0 & 0 \end{pmatrix}$. Some manipulations yield:

$$ \begin{aligned} \operatorname{OT} &= \frac{1}{4g_{21}^{2}} \begin{pmatrix} (2p_{12}\mu g + k_{21})^{2} & (2p_{12}\mu g + k_{21})(2p_{22}\mu g - 1 + k_{22}) \\ (2p_{22}\mu g - 1 + k_{22})(2p_{12}\mu g + k_{21}) & (2p_{22}\mu g - 1 + k_{22})^{2} \end{pmatrix} \\ &= \frac{1}{4g_{21}^{2}} \hat{S}_{0}\hat{S}_{0}^{\top} \end{aligned} \quad (34) $$
According to [107, Fact 5.11.5] and given that *OT* is positive semi-definite, $\sigma_i(OT) = \lambda_i(OT)$ for $i \in \{1, 2\}$. Then, $\sigma_{\max}(OT) = \frac{1}{4g_{21}^2}((2p_{12}\mu g + k_{21})^2 + (2p_{22}\mu g - 1 + k_{22})^2) > 0$. Given that $P_0 = P_0^\top > 0$, then $\lambda_{\max}(P_0) > 0$. Note that the matrices $\Delta\hat{Q}$ and $-(\hat{\Lambda}_1 + \hat{\Lambda}_B)$ are congruent, so they have the same number and signs of eigenvalues [108, Theroem 4.5.8]. This means that $\lambda_i(\Delta\hat{Q}) < 0$ for $i \in \{1, 2\}$. Assume that $\hat{\Lambda}_1 = \alpha I_2$ and $\hat{\Lambda}_B = \beta I_2$, $\alpha > 0$, $\beta > 0$, then $\Delta\hat{Q} = -(\alpha + \beta)P_0^2$ which is symmetric (i.e., $\Delta\hat{Q} = \Delta\hat{Q}^\top$). As $P_0 = P_0^\top > 0$, it follows that $\sigma_i(P_0) = \lambda_i(P_0)$. Thus, $\sigma_{\max}(\Delta\hat{Q}) = (\alpha + \beta)\lambda_{\max}(P_0)^2 - \epsilon\lambda_{\max}(P_0) + \sigma_{\max}(OT) < 0$. In order to find the conditions such that item (iv) in Proposition B.3 is satisfied, let us solve the inequality with the variable $\hat{\theta}$ as follows:

$$ (\alpha + \beta)\hat{\theta}^2 - \epsilon\hat{\theta} + \sigma_{\max}(OT) < 0 \quad (35) $$

The roots are

$$ \hat{\theta}_1 = \frac{1}{2(\alpha + \beta)} \left( \epsilon - \sqrt{\epsilon^2 - 4(\alpha + \beta)\sigma_{\max}(OT)} \right) \quad \text{and} \quad \hat{\theta}_2 = \frac{1}{2(\alpha + \beta)} \left( \epsilon + \sqrt{\epsilon^2 - 4(\alpha + \beta)\sigma_{\max}(OT)} \right) $$

Based on the solution of the quadratic inequality (35), and following the condition of item (iv) in Proposition B.3, $\lambda_{\max}(P_0)$ must satisfy the condition

$$ \hat{\theta}_1 < \lambda_{\max}(P_0) < \hat{\theta}_2 \quad (36) $$

where the roots $\hat{\theta}_1$ and $\hat{\theta}_2$ are real and positive if the discriminant: $\epsilon^2 - 4(\alpha + \beta)\sigma_{\max}(OT) \ge 0$, given that $\alpha + \beta > 0$. It is possible to find $\lambda_{\max}(P_0)$ satisfying (36) by properly constraining $\alpha$ and $\beta$, consequently constraining the uncertainty $\Delta\mu$. Additionally, $\epsilon$ can be chosen to determine the control gain $K_0$, which determines $\sigma_{\max}(OT)$. Thus, this item (i.e., $\hat{Q} \succcurlyeq \hat{S}\hat{R}^\dagger\hat{S}^\top$) is satisfied such that all the conditions in Proposition B.3 holds as explained above.

Therefore, the matrix inequality $\hat{L} \succcurlyeq 0$ in (33) holds. Given that $\frac{4\Delta\mu^2g^2d}{df - e^2} \le 1$ derived from the bound on $\Delta B$ in (32) where $\hat{\Lambda}_B = \begin{pmatrix} d & e \\ e & f \end{pmatrix} > 0$, let us assume that $\|\Delta E'\|^2 \le \rho$ where $\rho = \frac{df - e^2}{4d}$. Thus, the inequality in (31) is written as:

$$ \begin{aligned} \dot{V} & \le -\left(\frac{e}{\Delta\lambda}\right)^\top \left(-M_0 + \Delta M_0 - \begin{pmatrix} P_0 \hat{\Lambda}_1 P_0 & 0 \\ 0 & 0 \end{pmatrix}\right) \left(\frac{e}{\Delta\lambda}\right) + \Delta E'^\top \hat{\Lambda}_1^{-1} \Delta E' \\ & \le -\epsilon' e^\top P_0 e + E'^\top \hat{\Lambda}_1^{-1} \Delta E' \\ & \le -\epsilon' \lambda_{\min}(P_0) \|e\|^2 + \lambda_{\min}^{-1}(\hat{\Lambda}_1) \|\Delta E'\|^2 \le -\epsilon' \lambda_{\min}(P_0) \|e\|^2 + \lambda_{\min}^{-1}(\hat{\Lambda}_1) \rho. \end{aligned} \quad (37) $$

So, $\dot{V} < 0$ outside the ball $B_r(0) \subset \mathbb{R}^n$, with $r = \frac{\Delta\lambda^{-1}(\hat{\Lambda}_1)\rho}{\lambda_{\min}(P_0)\epsilon}$. Thus, the solution of (27) is GUUB [109, 110] by relaxing strong passivity and, according to [83, Proposition 4.5], the ultimate bound is given by:

$$ \|e\| \le \sqrt{\frac{\lambda_{\max}(P_0)\lambda_{\min}^{-1}(\hat{\Lambda}_1)\rho}{\lambda_{\min}^2(P_0)\epsilon'}} \quad (38) $$

□

**Remark 6.2.** It is important to note that taking $g_{21} = 0$ within the framework of bounded uncertainties in this section is not possible. This is due to the condition in item 1 ($\hat{R} \succcurlyeq 0$) to ensure that the matrix inequality $\hat{L} \succcurlyeq 0$ in (33) holds. Furthermore, it is necessary to select $g_{21} \ge \frac{1}{2}$ as imposed by item 1 in the above analysis. It is also worth noting that nowhere is used the fact that the set-valued part in right-hand side of (8) or (9) (equivalently, the multiplier $\lambda_1$) is bounded. The feedback controller structure respects the passivity-based structure as designed in [83].
## 6.2. Numerical simulations

Let us check on an example if there exists $P_0$, $K_0$ and $G_0$ such that the SSPMI:

$$M_0 \lesssim \begin{pmatrix} -\epsilon P_0 & 0 \\ 0 & 0 \end{pmatrix}, \qquad (39)$$

holds where $M_0$ is defined in (28). For numerical purposes the BMI in (39) is transformed into an LMI (see Appendix A.1). Take $\mu_0 = 0.5, m = 1, g = 9.8$, and $\epsilon = 10$. The software package MOSEK 9.3.14 [106] computes a solution given by:

$$K_0 = \begin{pmatrix} -70.79 & -11.95 \\ -29.89 & -4.22 \end{pmatrix}, \quad G_0 = \begin{pmatrix} 0 & 0 \\ 0.67 & 0 \end{pmatrix} \quad \text{and} \quad P_0 = \begin{pmatrix} 31.8 & 3.05 \\ 3.05 & 0.53 \end{pmatrix} \qquad (40)$$

In order to find the admissible $\Delta\mu$, let us solve the LMI of uncertainties $\hat{L} \succ 0$ in (33) with the unknowns $\hat{\Lambda}_1$ and $\hat{\Lambda}_B$ given the values of $K_0$, $P_0$ and $G_0$ from the strict passivity of the nominal system ($A_0+E_0K_0, B_0, C_0+F_0K_0, D_0+F_0G_0$). The solution is:

$$\hat{\Lambda}_1 = \begin{pmatrix} 0.16 & 0 \\ 0 & 9.2 \end{pmatrix} \quad \text{and} \quad \hat{\Lambda}_B = \begin{pmatrix} 0.14 & 0 \\ 0 & 0.567 \end{pmatrix} \qquad (41)$$

with the constraint trace($\hat{\Lambda}_1$) $\le 10$ and without constraining $\hat{\Lambda}_B$. The admissible value of the upperbound on the uncertainty is $\Delta\mu \le 0.038$. It is possible to increase the admissible $\Delta\mu$ for the same $\epsilon$ by increasing the values of $\hat{\Lambda}_1$ or $\hat{\Lambda}_B$ (i.e., increasing the trace of the matrices $\hat{\Lambda}_1$ and $\hat{\Lambda}_B$), maintaining the BMI $\hat{L} \succ 0$. Let us increase the value of $\hat{\Lambda}_B$ consequently $\Delta\mu$ such that trace($\hat{\Lambda}_B$) = 20, the numerical solution is:

$$\hat{\Lambda}_1 = \begin{pmatrix} 0.032 & 0 \\ 0 & 0.86 \end{pmatrix} \quad \text{and} \quad \hat{\Lambda}_B = \begin{pmatrix} 0.023 & 0 \\ 0 & 19.9 \end{pmatrix} \qquad (42)$$

The admissible value of the uncertainty is bounded such that $\Delta\mu \le 0.22$. The following figure illustrates the robustness of the closed-loop error system (16) through simulations, considering perturbation on the parameter $\mu_0$ by time-varying and bounded uncertainty $\Delta\mu$. Take $x_d(0) = (1, 0)^T$, $x(0) = (0, -1)^T$, and the time step $h = 0.01s$ (simulations are still obtained using the time-stepping method available in siconos [103, 104]).

<figure><img src="image_8.png" /><figcaption>Figure 5: Closed-loop error system with control gains from (40) and uncertainty $\Delta\mu(t) = 0.2 \sin \sqrt{3}t$.</figcaption></figure>

In Figure 5, the value of the uncertainty $\Delta\mu$ is bounded such that $\Delta\mu < 0.22$. The numerical simulation shows that the $L_\infty$ norms of the error $e = x - x_d$ are $\|e_1\|_{[5,15],\infty} = 0.028$ and $\|e_2\|_{[5,15],\infty} = 0.154$. It is also noticeable that the synchronization between the complementarity variables $\lambda$ and $\lambda_d$, as well as between $w$ and $w_d$ is lost in the presence of uncertainties.
# 7. Controller Design and Analysis with Stribeck Effects

As usual in control theory, it is interesting to analyse whether or not refining the model yields better results than designing a robust feedback controller. The foregoing study guarantees some robustness with respect to time-varying bounded uncertainties in the friction coefficient. Now, let us assume that the Coulomb’s model is enhanced to include Stribeck effects. The reason of this assumption is that on one hand, Stribeck effects often represent in a quite acceptable way the frictional effects, on the other hand it will allow us to propose different control approach, using its hypomonotonicity (see Definition 2). Actually, any other modification of Coulomb’s set-valued model which satisfies a hypomonotonicity constraint, would fit with the control framework in this section. This new approach consists of "monotonifying" a hypomonotone friction model, so that it fits within the foregoing framework. As noted in the introduction, Stribeck effects often play a major role [57, 111, 64, 74, 42, 38, 33].

## 7.1. The Stribeck Model and its Properties

One example of a physical model that may be considered as a bounded uncertainty, is the Stribeck model of sliding friction [42], which states that $\mu = \mu(\dot{q}_{rel})$, where $\dot{q}_{rel} \triangleq \dot{q} - u_2$ is the relative tangential velocity. It can be expressed as [112, Equation (17)]:

$$F_{\text{str}}(\dot{q}_{\text{rel}}(t)) = -F_n(\mu_s - \mu_c) \left( \exp\left(-\frac{\dot{q}_{\text{rel}}^2(t)}{v_s^2}\right) - 1 \right) \text{sgn}(\dot{q}_{\text{rel}}(t)) - \underbrace{\mu_s F_n \text{sgn}(\dot{q}_{\text{rel}}(t))}_{F_{\text{diff}}(\dot{q}_{\text{rel}})} \quad (43)$$

where $\mu_s > \mu_c > 0$, $v_s$ are parameters, $F_n = mg > 0$ is the constant normal contact force. The total force is therefore the sum of a continuously differentiable, Lipschitz continuous, hypomonotone decreasing term $F_{diff}(\dot{x})$, and a set-valued maximal monotone term $F_{sv}(\dot{x})$. Figure 6 shows the graphical representations of different frictional forces, including Coulomb friction as $F_{sv}$ and Stribeck friction $F_{str}$ in (43) with respect to the relative velocity $\dot{q}_{rel}$. The behaviour of the function $F_{diff}$ is depicted with variations in its parameter $v_s$. Take $F_n = mg = 9.8$, $\mu_s = 0.5$ and $\mu_c = 0.4$.

<figure><img src="image_9.png" /><figcaption>Figure 6: Different frictional forces as a function of relative velocity $\dot{q}_{rel}$</figcaption></figure>
The plot in Figure 6a shows the Coulomb friction force with respect to the relative velocity $\dot{q}_{rel}$. It varies inside $F_{sv} \in [-F_n\mu_s, F_n\mu_s]$. In Figure 6b, the function $F_{diff}$ is depicted for two parameters $v_s$, and $F_{diff} \in [-F_n(\mu_s - \mu_c), F_n(\mu_s - \mu_c)]$. The plot in Figure 6c displays the Stribeck friction force $F_{str}$ without viscous friction as in (43) for two different values of the parameter $v_s$ which causes faster decay towards the value $F_n(\mu_s - \mu_c)$ as it decreases. Let us analyze the behaviour of the function $F_{diff}(\dot{q}_{rel})$. Denote $\dot{q}_{rel}$ as $x$ for simplicity, then the function $F_{diff}(x)$ is expressed as follows:

$$F_{diff}(x(t)) = F_n(\mu_s - \mu_c) \left( \exp\left(-\frac{x^2(t)}{v_s^2}\right) - 1 \right) \operatorname{sgn}(x(t)).$$

The derivative of $F_{diff}(x)$ with respect to $x$ is given for $x \neq 0$ by:

$$\frac{d}{dx}F_{diff}(x) = F_n(\mu_s - \mu_c)\left(-\frac{2x}{v_s^2}\exp\left(-\frac{x^2}{v_s^2}\right)\right)\operatorname{sgn}(x)$$

For $x > 0$: $\frac{d}{dx}F_{diff}(x) = -F_n(\mu_s - \mu_c)\frac{2x}{v_s^2}\exp\left(-\frac{x^2}{v_s^2}\right)$. For $x < 0$: $\frac{d}{dx}F_{diff}(x) = F_n(\mu_s - \mu_c)\frac{2x}{v_s^2}\exp\left(-\frac{x^2}{v_s^2}\right)$. Then, $\lim_{x \to 0^-} \frac{d}{dx}F_{diff}(x) = \lim_{x \to 0^+} \frac{d}{dx}F_{diff}(x) = 0$. Thus, the function $\frac{d}{dx}F_{diff}(\cdot)$ is continuous at 0, and the function $F_{diff}(\cdot)$ is continuously differentiable everywhere, with bounded derivatives (bounded in absolute value). It follows that it is a hypomonotone single-valued mapping. Let us write the dynamics in (8) with Stribeck friction presented in (43), then:

$$\ddot{q}(t) \in \frac{1}{m}u_1(t) - \mu g \operatorname{sgn}(\dot{q}_{rel}(t)) - F_{diff}(\dot{q}_{rel}(t)) \pm \alpha \dot{q}_{rel}(t) \quad (44)$$

where the linear term $\alpha \dot{q}_{rel}$ is added and subtracted. Thus, the dynamics in (44) is written as follows:

$$\ddot{q}(t) \in \frac{1}{m}u'_1(t) - \mu g \operatorname{sgn}(\dot{q}_{rel}(t)) - (F_{diff} + \alpha I_d)(\dot{q}_{rel}(t)) \quad (45)$$

where $u_1 = u'_1 - m\alpha\dot{q}_{rel}$ and $I_d$ is the identity mapping. It is required to find a suitable gain $\alpha > 0$ such that the mapping $\dot{q}_{rel} \mapsto (F_{diff} + \alpha I_d)(\dot{q}_{rel})$ is monotone (it is maximal since it is single-valued). For this purpose, it is required to find the maximum slope of the function $F_{diff}(\cdot)$. The second derivative of $F_{diff}(x)$ is:

$$\frac{d^2}{dx^2}F_{diff}(x) = \frac{2}{v_s^4}F_n(\mu_s - \mu_c)\left((2x^2 - v_s^2)\exp\left(-\frac{x^2}{v_s^2}\right)\right)\operatorname{sgn}(x). \quad (46)$$

Let $m_{max}$ be the maximum slope of $F_{diff}(x)$ which is the absolute value of $\frac{d}{dx}F_{diff}(x)$ at $x = \frac{v_s}{\sqrt{2}}$ and is given by:

$$m_{max} = \frac{2}{v_s \sqrt{2}} F_n(\mu_s - \mu_c) \exp\left(-\frac{1}{2}\right). \quad (47)$$

Thus, the monotonicity of the mapping $\dot{q}_{rel} \mapsto (F_{diff} + \alpha I_d)(\dot{q}_{rel})$ is guaranteed by choosing the gain $\alpha$ such that $\alpha > m_{max}$. Take

$$\alpha = \sqrt{2} \exp(\frac{1}{2}) m_{max} = \frac{2}{v_s} F_n(\mu_s - \mu_c) \quad (48)$$

in the following.

**Remark 7.1.** If viscous friction is included in the force of friction given in (43), then the dynamics in (45) is rewritten as:

$$\ddot{q}(t) \in \frac{1}{m}u'_1(t) - \mu g \operatorname{sgn}(\dot{q}_{rel}(t)) - (F_{diff} + k_v + \alpha I_d)(\dot{q}_{rel}(t))$$

where $k_v > 0$ is the coefficient of viscous friction. Given that $k_v > 0$, the term $k_v\dot{q}_{rel}$ increases the overall damping by $k_v$, thus enhancing the dissipativity of the system which is crucial for our analysis. The main effect of adding viscous friction is that the value of the compensating term $\alpha$ changes such that $\alpha > \frac{2}{v_s\sqrt{2}}F_n(\mu_s - \mu_c)\exp(-\frac{1}{2}) + k_v$. Therefore, it is important to note that the following analytical approach applies to the Stribeck effect in (43) with viscous friction.
## 7.2. Controller Calculation

Let a suitable $\alpha$ be calculated. The next step consists of applying the foregoing passivity-based (or monotonicity-based) technique to calculate $u = (u'_1, u_2)^\top = K(x - x_d) + u_d$ with $x \stackrel{\Delta}{=} (q, \dot{q})^\top$ as independent input, then we have $u_1 = u'_1 - m\alpha(\dot{q} - u_2)$. Let us write the dynamics in (45) according to the expression of $F_{diff}$ in (43) as:

$$ 
\begin{aligned}
\ddot{q}(t) & \in \frac{1}{m} u'_1(t) - g\lambda_t (\mu_s + (\mu_s - \mu_c) \left( \exp\left(\frac{-(\dot{q}(t)-u_2(t))^2}{v_s^2}\right) - 1 \right)) - \alpha(\dot{q}(t) - u_2(t)) \\
& = \frac{1}{m} u'_1(t) - g\mu_s\lambda_t - F_{diff,\alpha}(\dot{q} - u_2),
\end{aligned} 
\quad (49) $$

where $\lambda_t \in \text{sgn}(\dot{q} - u_2)$ and

$$ F_{diff,\alpha}(\dot{q}_{rel}(t)) \stackrel{\Delta}{=} (F_{diff} + \alpha I_d)(\dot{q}_{rel}(t)) = F_n(\mu_s - \mu_c) \left( \left( \exp\left(\frac{-(\dot{q}_{rel}(t))^2}{v_s^2}\right) - 1 \right) + \frac{2}{v_s} \right) (\dot{q}_{rel}(t)) $$

is maximal monotone, where $\alpha$ is given in (48). Let us keep the plant dynamics as in (45) and write the desired dynamics in (14) with Stribeck friction in (43) as:

$$ 
\begin{aligned}
\ddot{q}_d(t) & \in \frac{1}{m} u_{1d}(t) - \mu_s g \text{sgn}(\dot{q}_{rel,d}(t)) - F_{diff}(\dot{q}_{rel,d}(t)) \pm \alpha \dot{q}_{rel,d}(t) \\
& \in \frac{1}{m} u'_{1d}(t) - \mu_s g \text{sgn}(\dot{q}_{rel,d}(t)) - (F_{diff} + \alpha I_d)(\dot{q}_{rel,d}(t))
\end{aligned} 
\quad (50) $$

where $u_{1d}(t) = u'_{1d}(t) - m\alpha\dot{q}_{rel,d}(t)$ and $\dot{q}_{rel,d}(t) \stackrel{\Delta}{=} \dot{q}_d(t) - u_{2d}(t)$. The mapping $\dot{q}_{rel,d}(t) \mapsto F_{diff,\alpha}(\dot{q}_{rel,d}(t)) \stackrel{\Delta}{=} (F_{diff} + \alpha I_d)(\dot{q}_{rel,d}(t))$ is maximal monotone.

Therefore, the Stribeck model is incorporated into the desired dynamics. In view of the philosophy followed in this section, this is a logical step.

## 7.3. Error Dynamics Stability Analysis

Using (50) and (45), the error $\ddot{q} = \ddot{q} - \ddot{q}_d$ is written as:

$$ 
\begin{aligned}
\ddot{q}(t) & \in \frac{1}{m}(u_1(t) - u_{1d}(t)) - \mu_s g \text{sgn}(\dot{q}_{rel}(t)) + \mu_s g \text{sgn}(\dot{q}_{rel,d}(t)) - F_{diff,\alpha}(\dot{q}_{rel}(t)) + F_{diff,\alpha}(\dot{q}_{rel,d}(t)) \\
& \in \frac{1}{m}(u'_1(t) - m\alpha(\dot{q}_{rel}(t) - \dot{q}_{rel,d}(t)) - u'_{1d}(t)) - \mu_s g(\text{sgn}(\dot{q}_{rel}(t)) - \text{sgn}(\dot{q}_{rel,d}(t))) \\
& \quad -(F_{diff,\alpha}(\dot{q}_{rel}(t)) - F_{diff,\alpha}(\dot{q}_{rel,d}(t))) \\
& \in \frac{1}{m}(k_{11}(q(t) - q_d(t)) + k_{12}(\dot{q}(t) - \dot{q}_d(t)) - m\alpha(\dot{q}(t) - k_{21}(q(t) - q_d(t)) - k_{22}(\dot{q}(t) - \dot{q}_d(t)) - \dot{q}_d(t))) \\
& \quad -\mu_s g(\text{sgn}(\dot{q}_{rel}(t)) - \text{sgn}(\dot{q}_{rel,d}(t))) - (F_{diff,\alpha}(\dot{q}_{rel}(t)) - F_{diff,\alpha}(\dot{q}_{rel,d}(t)))
\end{aligned} 
\quad (51) $$

with $\tilde{q} = q - q_d$ and $e = (\tilde{q}, \dot{q})^\top$. The error dynamics in (51) is equivalently written as:

$$ \dot{e}(t) \in \underbrace{\left(A_{cl} + \begin{pmatrix} 0 \\ \alpha \end{pmatrix} (C_{cl})_{1,\bullet}\right)}_{\stackrel{\Delta}{=} \hat{A}_{cl}} e(t) - \begin{pmatrix} 0 \\ \mu_s g \end{pmatrix} (\text{sgn}(\dot{q}_{rel}(t)) - \text{sgn}(\dot{q}_{rel,d}(t))) - \left(F_{diff,\alpha}(\dot{q}_{rel}(t)) - F_{diff,\alpha}(\dot{q}_{rel,d}(t))\right) \quad (52) $$

where $\lambda_{1d}(t) \in \text{sgn}(\dot{q}_d(t) - u_{2d}(t))$, $\lambda_1(t) \in \text{sgn}(-\langle C_{cl} \rangle_{1,\bullet} e(t) + \dot{q}_d(t) - u_{2d}(t))$.

**Proposition 7.2.** Consider the error dynamics in (52) with bounded initial conditions, and assume that the quadruple $(\hat{A}_{cl}, B_{cl}, C_{cl}, D_{cl})$ of the closed loop system is strictly state passive with $\epsilon > 0$. Then, the equilibrium point $e^* = 0$ is globally exponentially stable.
*Proof.* Consider the error dynamics in (52). Assume that the new closed-loop quadruple ($\hat{A}_{cl}$, $B_{cl}$, $C_{cl}$, $D_{cl}$) is strictly state passive. This means that the following SSPMI:

$$ \begin{pmatrix} \hat{A}_{cl}^{\top} \hat{P} + \hat{P} \hat{A}_{cl} & \hat{P} B_{cl} - C_{cl}^{\top} \\ B_{cl}^{\top} \hat{P} - C_{cl} & -D_{cl} - D_{cl}^{\top} \end{pmatrix} \preceq \begin{pmatrix} -\epsilon \hat{P} & 0 \\ 0 & 0 \end{pmatrix} \quad (53) $$

has a solution with $\hat{P} = \hat{P}^{\top} > 0$. Note that this SSPMI is transformed into an LMI (see Appendix A.1) to be solved numerically. Consider the storage function $V(e(t)) = e^{\top}(t) \hat{P} e(t)$. The rate of change of the storage function is:

$$ \dot{V}(e(t)) = e^{\top}(\hat{A}_{cl}^{\top} \hat{P} + \hat{P} \hat{A}_{cl})e + 2e^{\top} \hat{P} \begin{pmatrix} 0 \\ \mu_{sg} \end{pmatrix} [\operatorname{sgn}(\dot{q}_{rel,d}(t)) - \operatorname{sgn}(\dot{q}_{rel}(t))] - 2e^{\top} \hat{P} \begin{pmatrix} 0 \\ F_{diff,\alpha}(\dot{q}_{rel,d}) - F_{diff,\alpha}(\dot{q}_{rel}) \end{pmatrix} $$

Take $\hat{P} = \begin{pmatrix} \hat{p}_{11} & \hat{p}_{12} \\ \hat{p}_{21} & \hat{p}_{22} \end{pmatrix}$. Given that $\hat{P} B_{cl} = C_{cl}$ holds and that $(C_{cl})_{1,\bullet} e(t) = u_2(t) - \dot{q}(t) - u_{2d}(t) + \dot{q}_d(t)$, then the rate of change of the storage function is written as follows:

$$ \begin{aligned} \dot{V}(e) &= e^{\top}(\hat{A}_{cl}^{\top} \hat{P} + \hat{P} \hat{A}_{cl})e - e^{\top}(\hat{P} B_{cl})_{\bullet,1} [\operatorname{sgn}(\dot{q}_{rel,d}(t)) - \operatorname{sgn}(\dot{q}_{rel}(t))] - \frac{1}{\mu_{sg}} e^{\top}(\hat{P} B_{cl})_{\bullet,1} [F_{diff,\alpha}(\dot{q}_{rel,d}) - F_{diff,\alpha}(\dot{q}_{rel})] \\ &= e^{\top}(\hat{A}_{cl}^{\top} \hat{P} + \hat{P} \hat{A}_{cl})e + e^{\top}(C_{cl}^{\top})_{\bullet,1} [\operatorname{sgn}(\dot{q}_{rel}(t)) - \operatorname{sgn}(\dot{q}_{rel,d}(t))] + \frac{1}{\mu_{sg}} e^{\top}(C_{cl}^{\top})_{\bullet,1} [F_{diff,\alpha}(\dot{q}_{rel}) - F_{diff,\alpha}(\dot{q}_{rel,d})] \\ &= e^{\top}(\hat{A}_{cl}^{\top} \hat{P} + \hat{P} \hat{A}_{cl})e - [\dot{q}_{rel}(t) - \dot{q}_{rel,d}(t)]^{\top} [\operatorname{sgn}(\dot{q}_{rel}(t)) - \operatorname{sgn}(\dot{q}_{rel,d}(t))] \\ &\quad - [\dot{q}_{rel}(t) - \dot{q}_{rel,d}(t)]^{\top} [F_{diff,\alpha}(\dot{q}_{rel}) - F_{diff,\alpha}(\dot{q}_{rel,d})] \\ &\le e^{\top}(\hat{A}_{cl}^{\top} \hat{P} + \hat{P} \hat{A}_{cl})e \le -\epsilon e(t)^{\top} \hat{P} e(t) \end{aligned} \quad (54) $$

Therefore, the equilibrium point of the error dynamics in (52), $e^* = 0$, is globally exponentially stable since $\|e(t)\|^2 \le \frac{V(0)}{\lambda_{\min}(\hat{P})} \exp(-\epsilon t)$. $\square$

**Remark 7.3.** The above tracking controller aims at improving the closed-loop accuracy when the friction effects follow Stribeck's model and the parameters are correctly identified. Other control strategies are possible and not tackled here. For instance, the term $F_{diff}(\dot{q}_{rel})$ is uniformly bounded and thus could be compensated for with a simple first-order sliding-mode input, hence dispensing the designer for the knowledge of parameters $\mu_c$ and $v_s$ (see (43)). This would however be at the price of applying an input $u_1$ which is no longer a feedforward input (in case $k > 0$ and $f > 0$). It would also allow us to design another desired system, as in (50) without the nonlinear term $F_{diff}(\dot{q}_{rel,d}(t))$. For the sake of briefness this is not developed further.

**Remark 7.4.** Starting from (49), and taking into account the maximal monotonicity of the mapping $F_{diff,\alpha}(\cdot)$, the arbitrary desired trajectory design in section 4.3 can be applied. Same comments about pros and cons of both approaches can be made in case of Stribeck effects. This shows in passing that the proposed control strategy applies to a much broader class of set-valued mappings, than mere signum multifunction.

### 7.4. Numerical simulation

Recall that $x \triangleq (q, \dot{q})$, $x_d \triangleq (q_d, \dot{q}_d)$, and $e \triangleq (\tilde{q}, \dot{\tilde{q}})$. The desired dynamics in (50), the closed-loop dynamics in (45) and the error dynamics in (52) can be equivalently expressed in the form of LCS as explained in section 3.1. The following numerical simulation is obtained using the LCS time-stepping scheme available in siconos [103, 104]. Take $m = 1$kg, $F_n = mg = 9.8m/s^2$, $v_s = 0.2$, $\mu_s = 0.5$ and $\mu_c = 0.4$. Let us check if there exists a control gain $K$ such that the closed-loop system's quadruple $(A + EK, B, C + FK, D)$ is strictly state passive. The LMI derived from the BMI of strict state passivity in (17) has a solution which is given by:

$$ K = \begin{pmatrix} -70.79 & -11.95 \\ -23.91 & -3.18 \end{pmatrix} \quad \text{and} \quad P = \begin{pmatrix} 31.8 & 3.05 \\ 3.05 & 0.53 \end{pmatrix} $$
Consider the value of the linear term $\alpha$ in (48), and the controller $u_1 = k_{11}(x_1 - x_{1d}) + k_{12}(x_2 - x_{2d}) + u_{1d} - m\alpha(x_2 - u_2)$ and $u_2 = k_{21}(x_1 - x_{1d}) + k_{22}(x_2 - x_{2d}) + u_{2d}$. The initial conditions are $x_d(0) = (-1, -2)^\top$ and $x(0) = (1, 3)^\top$, with desired inputs $u_{1d} = 0.7\mu g \cos 5t$ and $u_{2d} = 0.5\mu g \sin 3t$, and a time step $h = 0.01s$.

<figure><img src="image_10.png" /><figcaption>Figure 7: The desired, closed-loop and error systems of the LCS in (45), complementarity variables $(\lambda, w)$ and input $u = (u_1, u_2)^\top$.</figcaption></figure>

The numerical simulation shows the desired trajectories generated by (50). The numerical simulation in Figure 7 shows that the tracking error $e$ converges exponentially to zero, which confirms the result stated in Proposition 7.2 when the friction model includes Stribeck effects. It can be observed that the plant's complementarity variables $(\lambda, w)$ converge to the desired complementarity variables $(\lambda_d, w_d)$.

# 8. Extensions towards Systems with n Degrees of Freedom

It is of interest to investigate how the foregoing material may extend to more complex mechanical systems with frictional bilateral contacts.

## 8.1. Two-mass Oscillator: Center of Gravity Tracking Control

This section is dedicated to analyze an extension of the above oscillator, when several (here, two) masses are sliding on the belt [38].

### 8.1.1. Center of Mass' Dynamics

Let us now consider the extension of (9) as depicted in Fig. 8:

$$ \left\{ \begin{array}{l} (a) \quad \ddot{q}_1(t) \in \frac{u_1(t)}{m_1} - \frac{k}{m_1}(q_1 - q_2) - \frac{f}{m_1}(\dot{q}_1 - \dot{q}_2) \\ \qquad \qquad -\mu_1 g \operatorname{sgn}(\dot{q}_1(t) - r\dot{\theta}(t)) \\ (b) \quad \ddot{q}_2(t) \in \frac{u_2(t)}{m_2} + \frac{k}{m_2}(q_1 - q_2) + \frac{f}{m_2}(\dot{q}_1 - \dot{q}_2) \\ \qquad \qquad -\mu_2 g \operatorname{sgn}(\dot{q}_2(t) - r\dot{\theta}(t)) \\ (c) \quad I\ddot{\theta}(t) \in \tau(t) + r\mu_1 m_1 g \operatorname{sgn}(\dot{q}_1(t) - r\dot{\theta}(t)) \\ \qquad \qquad + r\mu_2 m_2 g \operatorname{sgn}(\dot{q}_2(t) - r\dot{\theta}(t)), \end{array} \right. \qquad (55) $$
<figure><img src="image_11.png" /><figcaption>Figure 8: Frictional oscillator with two masses.</figcaption></figure>

This yields:

$$ egin{cases} (m_1 + m_2)\ddot{q}_c(t) \in u_1(t) + u_2(t) - \mu_1 m_1 g \\ \operatorname{sgn}(\dot{q}_1(t) - r\dot{\theta}(t)) - \mu_2 m_2 g \operatorname{sgn}(\dot{q}_2(t) - r\dot{\theta}(t)) \\ I\ddot{\theta}(t) \in \tau(t) + r\mu_1 m_1 g \operatorname{sgn}(\dot{q}_1(t) - r\dot{\theta}(t)) \\ +r\mu_2 m_2 g \operatorname{sgn}(\dot{q}_2(t) - r\dot{\theta}(t)), \end{cases} \quad (56) $$

with $q_c \triangleq \frac{m_1 q_1 + m_2 q_2}{m_1 + m_2}$ the masses' center of gravity position. After some manipulations we obtain:

$$ \left\{ \begin{array}{l} (m_1 + m_2)\ddot{q}_c(t) \in u_1(t) + u_2(t) - \mathcal{F}_t(-u_3) \\ I\ddot{\theta}(t) \in \tau(t) + r\mathcal{F}_t(-u_3) \end{array} \right. \quad (57) $$

where the mapping $-u_3 \mapsto \mathcal{F}_t(-u_3) = \mu_1 m_1 g \operatorname{sgn}(\dot{q}_1 - u_3) + \mu_2 m_2 g \operatorname{sgn}(\dot{q}_2 - u_3)$, $u_3 = r\dot{\theta}$, is maximal monotone [105, Corollary 24.4]. It can be computed as a set-valued piecewise-constant mapping:

$$ \left\{ \begin{array}{ll} -\mu_1 m_1 g - \mu_2 m_2 g & \text{if } u_3 > \dot{q}_2(t) > \dot{q}_1(t) \\ \mu_1 m_1 g + \mu_2 m_2 g & \text{if } u_3 < \dot{q}_1(t) < \dot{q}_2(t) \\ -\mu_1 m_1 g + \mu_2 m_2 g & \text{if } \dot{q}_1(t) < u_3 < \dot{q}_2(t) \\ [-\mu_1 m_1 g - \mu_2 m_2 g, -\mu_1 m_1 g + \mu_2 m_2 g] & \text{if } u_3 = \dot{q}_2(t) \\ [-\mu_1 m_1 g + \mu_2 m_2 g, \mu_1 m_1 g + \mu_2 m_2 g] & \text{if } u_3 = \dot{q}_1(t) \end{array} \right. \quad (58) $$

and similarly if $\dot{q}_1(t) > \dot{q}_2(t)$. For each fixed $t$ there exists a convex proper continuous function $f_t(\cdot)$ such that $\mathcal{F}_t(-u_3) = \partial f_t(-u_3)$, where $\partial$ is the subdifferentiation of Convex Analysis (see Definition 3). The analogy between (9) and (57) (58) is clear. Therefore, the same control strategy can be applied to (57) (58) for trajectory tracking control of the center of mass position and velocity. To this aim only one input $u_1$ or $u_2$ is sufficient. Extension to $n \ge 3$ aligned blocks follows the same lines. One significant discrepancy compared with the one-mass oscillator, is that the desired dynamics when designed as above, depends on $\dot{q}_1$ and $\dot{q}_2$ since $\mathcal{F}_t$ does. Indeed a natural choice is:

$$ (m_1 + m_2)\ddot{q}_{cd}(t) \in u_{1d}(t) - \mathcal{F}_t(-u_{3d}(t)), \quad (59) $$

where $q_{cd}: \mathbb{R}_+ \to \mathbb{R}$ is the desired center of mass motion. This allows us to redo the analysis in the previous sections. However, it is noteworthy that contrarily to (13), the dynamics (59) is not autonomous, thus the design of $q_{cd}(\cdot)$ is not possible from it. The problem is more a synchronization problem between (59) and (57) than a trajectory control problem. To cope with this issue we may design the desired system as follows, where the dependence of $\mathcal{F}_t$ on $\dot{q}_1$ and $\dot{q}_2$ is replaced by dependence on signals $\dot{q}_{1d}$ and $\dot{q}_{2d}$. To that end let us rewrite (56) (equivalently (57) (58)) in its LCS
form as (time argument is dropped):

$$
\begin{align}
\dot{x}_c &= \underbrace{\begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}}_{\triangleq A_c} x_c + \underbrace{\begin{pmatrix} 0 \\ \frac{1}{m_1+m_2} \end{pmatrix}}_{\triangleq E_c} \begin{pmatrix} 0 \\ 0 \end{pmatrix} u_c + \underbrace{\begin{pmatrix} 0 & 0 \\ \frac{-2\mu_1gm_1}{m_1+m_2} & 0 \end{pmatrix}}_{\triangleq B_c} \begin{pmatrix} 0 \\ \frac{-2\mu_2m_2g}{m_1+m_2} \end{pmatrix} \lambda_c + \underbrace{\begin{pmatrix} 0 \\ \frac{(\mu_1m_1+\mu_2m_2)g}{m_1+m_2} \end{pmatrix}}_{\triangleq E'_c} \nonumber \\
0 \le \lambda_c \perp w_c &= \underbrace{\begin{pmatrix} 0 & 1 \\ 0 & 0 \\ 0 & 1 \\ 0 & 0 \end{pmatrix}}_{\triangleq F_c} u_c + \underbrace{\begin{pmatrix} 0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -1 & 0 \end{pmatrix}}_{\triangleq D_c} \lambda_c + \underbrace{\begin{pmatrix} 0 & \frac{m_2}{m_1} \\ 0 & 0 \\ \frac{m_1}{m_2} & 0 \\ 0 & 0 \end{pmatrix}}_{\triangleq H_c} \begin{pmatrix} \dot{q}_1 \\ \dot{q}_2 \end{pmatrix} + F'_c + \underbrace{\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \end{pmatrix}}_{\triangleq C_c} \begin{pmatrix} -\frac{m_1+m_2}{m_1} \\ \frac{-m_1+m_2}{m_2} \\ 0 \\ 0 \end{pmatrix} x_c \ge 0, \tag{60}
\end{align}
$$

with $x_c = (q_c \dot{q}_c)^\top$, $\lambda_c = (\lambda_1 \ \lambda_2 \ \lambda_3 \ \lambda_4)^\top$, $u_c = (u_{c1} \ u_{c2})^\top = (u_1 + u_2 \ u_3)^\top$, $F'_c = (0 \ 1 \ 0 \ 1)^\top$. The dependence of $\mathcal{F}_1$ on velocities is visible through $H_c \dot{q}_i = (\dot{q}_1 \ \dot{q}_2)^\top$. The idea is to define the desired dynamics using the complementarity constraint:

$$0 \le \lambda_{cd} \perp w_{cd} = F_c u_{cd} + D_c \lambda_{cd} + C_c x_{cd} + F'_c + H_c \dot{q}_d \ge 0, \quad (61)$$

while the dynamical part is set as:

$$\dot{x}_{cd} = A_c x_{cd} + E_c u_{cd} + E'_c + B_c \lambda_{cd}. \quad (62)$$

Clearly $\dot{q}_{1d}$, $\dot{q}_{2d}$ and $\dot{q}_{cd}$ are not independent variables.

**Remark 8.1.** The above developments extend to the case where m masses slide on the belt. Then a similar dynamics as in (56) holds for the gravity center of the m masses. The mapping $\mathcal{F}_1$ is then the sum of m maximal monotone mappings with same domains, hence it is again maximal monotone.

### 8.1.2. Closed-loop Error Dynamics and Stability Analysis

The error dynamics with $u_c = K(x_c - x_{cd}) + u_{cd}$ is given by:

$$
\left\{
\begin{array}{l}
\dot{\tilde{x}}_c(t) = (A_c + E_c K)\tilde{x}_c(t) + B_c(\lambda_c(t) - \lambda_{cd}(t)) \\
\quad 0 \le \lambda_{cd}(t) \perp w_{cd}(t) = F_c u_{cd}(t) + D_c \lambda_{cd}(t) + C_c x_{cd}(t) \\
\qquad + F'_c + H_c \dot{q}_d(t) \ge 0 \\
\\
\quad 0 \le \lambda_c(t) \perp w_c(t) = (C_c + F_c K)\tilde{x}_c(t) + F_c u_{cd}(t) \\
\qquad + C_c x_{cd} + D_c \lambda_c(t) + F'_c + H_c \dot{q}(t) \ge 0
\end{array}
\right.
\qquad (63)
$$

**Proposition 8.2.** where $\tilde{x}_c = x_c - x_{cd}$. Assume that the quadruple $(A_c + E_c K, B_c, C_c + F_c K, D_c)$ is strictly state passive, and that $\dot{q}_1$ is bounded. Then $(\tilde{q}_c, \dot{q}_c)^\top$ is globally uniformly ultimately bounded, and $\dot{q}_2$ is bounded.

*Proof.* Let $V_c(\tilde{x}_c) = \frac{1}{2} \tilde{x}_c^\top P_c \tilde{x}_c$, where $P_c = P_c^\top > 0$ is a solution of SSPMI$(A_c + E_c K, B_c, C_c + F_c K, D_c, \epsilon)$. We obtain along trajectories of (63):

$$
\begin{align}
\dot{V}_c(t) &= \tilde{x}_c^\top (P_c (A_c + E_c K) + (A_c + E_c K)^\top P_c) \tilde{x}_c + \tilde{x}_c^\top P_c B_c (\lambda_c - \lambda_{cd}) \nonumber \\
&\le -\epsilon \tilde{x}_c^\top P_c \tilde{x}_c + \tilde{x}_c^\top (C_c + F_c K)(\lambda_c - \lambda_{cd}) \nonumber \\
&\le -\epsilon \tilde{x}_c^\top P_c \tilde{x}_c + (w_c - w_{cd} + D_c(\lambda_{cd} - \lambda_c) - H_c \dot{q})^\top (\lambda_c - \lambda_{cd}) \nonumber \\
&\le -\epsilon \tilde{x}_c^\top P_c \tilde{x}_c + (w_c - w_{cd} - H_c \dot{q})^\top (\lambda_c - \lambda_{cd}) \le -\epsilon \tilde{x}_c^\top P_c \tilde{x}_c - \dot{q}^\top H_c^\top (\lambda_c - \lambda_{cd}) \tag{64}
\end{align}
$$

From (60) and (56) it is inferred that $H_c^\top (\lambda_c - \lambda_{cd})$ is bounded, since it involves only $\lambda_1$ and $\lambda_3$ which are selections of signum functions. We also have $\dot{q}_c = \frac{(m_1 m_2) \dot{q}}{m_1+m_2}$, $\tilde{x}_c = (\tilde{q}_c \ \dot{q}_c)^\top$. Thus:

$$
\begin{align}
\dot{V}_c(t) &\le -\epsilon \tilde{x}_c^\top P_c \tilde{x}_c + \sqrt{\mu_1^2 m_2^2 + \mu_2^2 m_1^2} g \|\dot{q}\| \nonumber \\
&\le -\epsilon_1(\epsilon) \tilde{q}_c^2 - \epsilon_2(\epsilon) \tilde{q}_c^2 + \epsilon_3(|\dot{q}_1| + |\dot{q}_2|) \nonumber \\
&\le -\epsilon_1(\epsilon) \tilde{q}_c^2 - \epsilon_2(\epsilon) \tilde{q}_c^2 + \epsilon_4 |\dot{q}_1| + \epsilon_5 |\dot{q}_c| \nonumber \\
&\le -\epsilon_1(\epsilon) \tilde{q}_c^2 - \epsilon_2(\epsilon) |\dot{q}_c| (|\dot{q}_c| - \frac{\epsilon_4}{\epsilon_2(\epsilon)} + \epsilon_4 |\dot{q}_1|_{\max}) \tag{65}
\end{align}
$$

where $\epsilon_i$ are positive constants, $\epsilon_i(\epsilon) \to +\infty$ as $\epsilon \to +\infty$, $i = 1, 2$.

□
**Remark 8.3.** It appears clearly in (64) why setting $H_c \dot{q}$ instead of $H_c \dot{q}_d$ in (61) simplifies the analysis. The price to pay is to prevent the generation of unwanted desired trajectories.

The inputs $u_{c1}$ and $u_3$ can be used to design the above controller $u_c$. The input $u_1$ can be used to guarantee a bounded $\dot{q}_1$ in (55) (a). Then $u_2 = u_{c1} - u_1$. It is noteworthy that $f$ and $k$ in Figure 8 do not play any role in the above analysis and could be set to zero. However it plays a role in the boundedness of the first mass' state. Assume that $k=0$ and $k=0$. Then (55) (a) is rewritten as:

$$ \ddot{q}_1(t) \in \frac{u_1(t)}{m_1} - \mu_1 g \operatorname{sgn}(\dot{q}_1(t) - r\dot{\theta}(t)) \quad (66) $$

The boundedness of $\dot{q}_1$ is guaranteed by any controller of the form $u_1 = -k_1(q_1 - q_{1d}) - f_1(\dot{q}_1 - \dot{q}_{1d}) + \ddot{q}_{1d}$. Similarly as in the one-mass system, the spring-dashpot terms can be realized passively by attaching the mass $m_1$ to the wall (see Figure 1b)), so that $u_1$ is a feedforward input $u_1(q_{1d}, \dot{q}_{1d}, \ddot{q}_{1d})$. Thus we arrived to the following.

**Proposition 8.4.** Consider the two-mass system in (55), with $k=0$, $f=0$, and mass $m_1$ attached to the wall with a spring-dashpot with stiffness $k_1 > 0$, damping $f_1 > 0$. Let $u_1 = k_1 q_{1d} + f_1 \dot{q}_{1d} + \ddot{q}_{1d}$, $u_c = (u_{c1}, u_{c2})^\top = K(x_c - x_{cd}) + u_{cd}$, $u_2 = u_{c1} - u_1$, $u_{c2} = u_3$, and assume that the quadruple $(A_c + E_c K, B_c, C_c + F_c K, D_c)$ is strictly state passive with constant $\epsilon > 0$. Then the state of the error dynamics in (63) is globally uniformly ultimately bounded in a ball whose radius $r(\epsilon)$ strictly decreases as $\epsilon$ strictly increases. Moreover both $(q_1, \dot{q}_1)^\top$ and $(q_2, \dot{q}_2)^\top$ are bounded.

The stiffness $k_1$ and the damping coefficient $f_1$ are attached between the wall and mass $m_1$, similarly as $k$ and $f$ in Fig. 1b.

**Remark 8.5.** All systems in (8), (9), (10), (44), (60) and (57) (58) belong to set-valued Lur'e dynamical systems [88, 84] with controllers. A block-diagram representation is as in Fig. 9. In this article the set-valued feedback may be the signum function (8), (9), or an operator defined from complementarity constraints (10), or the multifunction in (58). The closed-loop error systems (16), (27), (52), (63) can also be recast in the set-valued Lur'e class of systems (without control inputs). It is conjectured that the material in this paper can be extended to more general classes of maximal monotone set-valued feedback mappings.

<figure><img src="image_12.png" /><figcaption>Figure 9: Block-diagram of the control system.</figcaption></figure>

## 8.2. Lagrangian Systems with Frictional Bilateral Contacts

The foregoing oscillators have linear one-degree-of-freedom mass dynamics. Let us investigate a more general case of an $n$-degree-of-freedom system subjected to $m$ frictional bilateral contacts, where it is assumed that the frictional elements can be given arbitrary velocities (like the belt for frictional oscillators). We may consider two main classes of systems in Lagrangian form, which are extensions of (8) (an LCS formulation as in (10) is also possible). One with linear dynamics:

$$ M\ddot{q}_1 + K_v\dot{q}_1 + K_p q_1 \in u_1 - H_1[\mu_i|\lambda_{n,i}] \operatorname{Sgn}(H_1^\top \dot{q}_1 + H_2^\top \dot{q}_2), \quad (67) $$
where $K_v > 0$ and $K_p > 0$ are damping and stiffness matrices, and one with nonlinear dynamics:

$$M(q_1)\ddot{q}_1 + C(q_1, \dot{q}_1)\dot{q}_1 + G(q_1) \in u_1 - H_1(q)[\mu_i|\lambda_{n,i}] \text{ Sgn}(H_1(q)^\top \dot{q}_1 + H_2(q)^\top \dot{q}_2), \quad (68)$$

where: $q_1(t) \in \mathbf{R}^{n_1}$, $[\mu_i|\lambda_{n,i}] = \text{diag}(\mu_i|\lambda_{n,i}) \in \mathbf{R}^{m \times m}$, $\mu_i \ge 0$, $\lambda_{n,i}$ is the normal contact force at contact $i$, Sgn($x$) = (sgn($x_1$), ..., sgn($x_m$))$^\top$, $H_1(q) \in \mathbf{R}^{n_1 \times m}$, $H_2(q) \in \mathbf{R}^{n_2 \times m}$, $q_2 \in \mathbf{R}^{n_2}$ represents the positions of frictional elements (like the belt angle in the foregoing oscillators). Smooth nonlinearities are inertia matrix $M(q_1)$, Coriolis and centrifugal forces $C(q_1, \dot{q}_1)\dot{q}_1$, and forces that derive from a smooth potential $G(q_1)$ [88]. The form of the generalized friction force in (68) stems from virtual work principle [82, section 5.5] [113, 114, 37], where the vector of local relative tangential velocities satisfies $v_t = H_1(q)^\top \dot{q}_1 + H_2(q)^\top \dot{q}_2$. Mimicking the above control strategy we may consider that $\dot{q}_2$ is a control input $u_2$. The normal contact forces $\lambda_{n,i}$ in (68) are associated with $m$ bilateral holonomic constraints $h(q_1, q_2) = 0$. The dynamics in (67) and (68) are therefore reduced order dynamics, just as the mass dynamics in (8) or (10) is with respect to the mass plus pulley-belt dynamics.

The DI in (67) is equivalently rewritten as:

$$\begin{align}
\dot{x} &= \begin{pmatrix} 0 & I \\ M^{-1}K_p & M^{-1}K_v \end{pmatrix} x + \begin{pmatrix} 0 \\ M^{-1}H_1[\mu_i|\lambda_i] \end{pmatrix} e + \begin{pmatrix} 0 & 0 \\ M^{-1} & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ -2M^{-1}H_1[\mu_i|\lambda_i] & 0 \end{pmatrix} \begin{pmatrix} \lambda_1 \\ \lambda_2 \end{pmatrix} \tag{69} \\
0 &\le \begin{pmatrix} \lambda_1 \\ \lambda_2 \end{pmatrix} + \begin{pmatrix} -(0 & H_1^\top) \\ 0 & 0 \end{pmatrix} x + \begin{pmatrix} 0 & H_2^\top \\ 0 & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} + \begin{pmatrix} 0 \\ e \end{pmatrix} + \begin{pmatrix} 0 & I \\ -I & 0 \end{pmatrix} \begin{pmatrix} \lambda_1 \\ \lambda_2 \end{pmatrix} \ge 0, \nonumber
\end{align}$$

which is the counterpart of (10). Clearly a desired dynamics generator can be designed for (67) in a similar way as (13), still following [81] (see also [84]). Then same steps as in section 5 can be taken, inspired from [83], to design feedback gains from suitable SSPMI. The major assumption to be made is that there exists $u_2$ such that $H_1^\top \dot{q}_1 + H_2^\top u_2 = K_{21}\tilde{q}_1 + K_{22}\dot{\tilde{q}}_1 + u_{2d}(t)$, if a linear input as in (11) is sought. This implies conditions on $H_1$ and $H_2$. The counterpart of (14) is given as:

$$M\ddot{q}_{1d} + K_v\dot{q}_{1d} + K_p q_{1d} \in u_{1d} - H_1[\mu_i|\lambda_{n,i}] \text{ Sgn}(H_1^\top \dot{q}_{1d} + H_2^\top u_{2d}). \quad (70)$$

The case of (68) requires special attention because it may not be possible to define a generator as (13) and a linear controller as (11) to achieve tracking, because of smooth nonlinearities. The approach in section 4.3 may be a solution to the tracking problem in this case, as outlined in the next result.

**Proposition 8.6.** Consider (68). Let $q_{1d}(\cdot)$ be a desired trajectory, such that $\dot{q}_{1d}(\cdot)$ and $\ddot{q}_{1d}(\cdot)$ exist almost everywhere. Assume that there exists a radially unbounded positive definite function $V(\tilde{q}_1, \dot{\tilde{q}}_1)$ and a collocated feedback $u_1(q_1, \dot{q}_1, q_{1d}, \dot{q}_{1d}, \ddot{q}_{1d})$ such that $\frac{\partial V}{\partial \tilde{q}_1} \dot{\tilde{q}}_1 + \frac{\partial V}{\partial \dot{\tilde{q}}_1} M^{-1}(q_1)(-C(q_1, \dot{q}_1)\dot{q}_1 - G(q_1) + u_1(q, \dot{q}, q_{1d}, \dot{q}_{1d}, \ddot{q}_{1d}) - \ddot{q}_{1d}) \le -\gamma_1(\|\tilde{q}_1\|) - \gamma_2(\|\dot{\tilde{q}}_1\|)$, with $\gamma_1$ and $\gamma_2$ functions of class $\mathcal{K}$. Assume further that there exists $u_2(q_1, \dot{q}_1, q_{1d}, \dot{q}_{1d}, \ddot{q}_{1d})$ such that $H_1(q)^\top \dot{q}_1 + H_2(q)^\top u_2 = H_1(q)^\top M^{-1}(q_1) \frac{\partial V}{\partial \dot{\tilde{q}}_1}$. Then the closed-loop system is globally asymptotically Lyapunov stable.

*Proof.* We have the following:

$$\begin{align}
\dot{V}(\tilde{q}_1, \dot{\tilde{q}}_1) &= \frac{\partial V}{\partial \tilde{q}_1} \dot{\tilde{q}}_1 + \frac{\partial V}{\partial \dot{\tilde{q}}_1} (\ddot{q}_1 - \ddot{q}_{1d}) \nonumber \\
&= \frac{\partial V}{\partial \tilde{q}_1} \dot{\tilde{q}}_1 - \frac{\partial V}{\partial \dot{\tilde{q}}_1} \ddot{q}_{1d} \nonumber \\
&\quad + \frac{\partial V}{\partial \dot{\tilde{q}}_1} M^{-1}(q_1) [-C(q_1, \dot{q}_1)\dot{q}_1 - G(q_1) + u_1(q, \dot{q}, q_{1d}, \dot{q}_{1d}, \ddot{q}_{1d})] - H_1(q)[\mu_i|\lambda_{n,i}] \text{ Sgn}(H_1(q)^\top \dot{q}_1 + H_2(q)^\top \dot{q}_2)] \nonumber \\
&\le -\frac{\partial V}{\partial \dot{\tilde{q}}_1} M^{-1}(q_1)H_1(q)[\mu_i|\lambda_{n,i}] \text{ Sgn}(H_1(q)^\top \dot{q}_1 + H_2(q)^\top \dot{q}_2) \nonumber \\
&= -\frac{\partial V}{\partial \dot{\tilde{q}}_1} M^{-1}(q_1)H_1(q)[\mu_i|\lambda_{n,i}] \text{ Sgn}(H_1(q)^\top M^{-1}(q_1) \frac{\partial V}{\partial \dot{\tilde{q}}_1}^\top). \tag{71}
\end{align}$$

The mapping $x \mapsto [\mu_i|\lambda_{n,i}] \text{ Sgn}(x)$ is maximal monotone for each $t \ge 0$ (recall that $\lambda_{n,i} = \lambda_{n,i}(t)$ in general). Let us use the fact that $x^\top \text{Sgn}(x) \ge 0$ for any $x \in \mathbf{R}^m$, since $(0, 0)$ belongs to the graph of Sgn($\cdot$) and this is a maximal monotone mapping, and the fact that the mapping $x \mapsto [\mu_i|\lambda_{rmni}] \text{ Sgn}(x)$ is maximal monotone. It follows from [86, Lemma 12.43] that $\dot{V}(\tilde{q}_1, \dot{\tilde{q}}_1) \le -\gamma_1(\|\tilde{q}_1\|) - \gamma_2(\|\dot{\tilde{q}}_1\|)$. $\square$
## 9. Conclusions

This article deals with trajectory tracking control of frictional oscillators, which are low-dimensional strongly nonsmooth and set-valued systems made of a mass sliding on a moving belt. Passivity (feedback passification) and maximal monotonicity of the friction model are at the core of the stability analysis, which relies on incremental passivity. The overall objective is not to compensate for friction before setting a tracking controller designed for frictionless systems, but to use its monotonicity for control purposes. This is, to the best of authors' knowledge, an original approach. The problem is solved with 2-dimensional Coulomb's friction, both when the friction coefficient is known, and when uncertainties are present where robustness is carefully studied. It is extended to the case of varying friction coefficient (Stribeck's model), where the hypomonotonicity of the model is used. Extensions are proposed to multibody systems with multiple frictional contacts. The case of two masses moving on a belt is studied, where the objective is to control the masses' center of gravity. How the approach extends to Lagrangian systems is also investigated. Numerical simulations illustrate the theoretical results.

# Appendix A. Useful Mathematical Results

## Appendix A.1. Transformation of BMI to an LMI

The bilinear matrix inequality BMI in (17) can be reformulated to an LMI [115, Example 5.5 p.136]. The left-hand and the right-hand sides of (17) are multiplied by $\begin{pmatrix} Q & 0 \\ 0 & I \end{pmatrix}$, where $Q \triangleq P^{-1}$, and let $N \triangleq KQ$. It follows that

$$ \begin{pmatrix} QA^\top + AQ + N^\top E^\top + EN + \epsilon Q & B + EG - QC^\top - N^\top F^\top \\ B^\top + G^\top E^\top - CQ - FN & -D - FG - (D + FG)^\top \end{pmatrix} \preceq 0 \quad (A.1) $$

This LMI is a feasible problem and it is solved numerically in the new variables $Q = Q^\top > 0$, $G$, and $N$ where $K = NQ^{-1}$.

## Appendix A.2. Conditions for Positive Semi-Definiteness

**Lemma A.1.** [107, Proposition 8.2.4, Remark after Fact 6.4.5] Consider the matrix $M = \begin{pmatrix} Q & S \\ S^\top & R \end{pmatrix}$. Assume that $Q = Q^\top$ and $R = R^\top$. Then, the matrix $M$ is positive semi-definite (i.e. $M \succcurlyeq 0$) if and only if:

1. $R \succcurlyeq 0$,

2. $SR^\dagger R = S \Leftrightarrow \mathrm{Im}(S^\top) \subseteq \mathrm{Im}(R)$,

3. $Q \succcurlyeq SR^\dagger S^\top$.

**Fact 1.** [116] Let $R^{m \times m} \ni M = M^\top > 0$, $x \in R^m$, $q \in R^m$ and $K \subseteq R^m$ be a closed nonempty convex set. Then:

$$ M(x-q) \in -N_K(x) \Leftrightarrow x = \mathrm{proj}_M[K; q] \Leftrightarrow x = \underset{z \in K}{\operatorname{argmin}} \frac{1}{2}(z-q)^\top M(z-q), \quad (A.2) $$

where $\mathrm{proj}_M$ is the orthogonal projection in the metric defined by $M$. Moreover such a $x$ is unique.

# Appendix B. Material on Relaxing Strong Passivity to Strict State Passivity

As noted earlier, the closed-loop system in (12) cannot be made strongly passive (see Remark 3.3). This motivates us to apply the theoretical results on relaxing the strong passivity property from [83, section 4.3]. Let us recall some of these results, which mainly rely on the positive semi-definiteness conditions presented in Definition A.1.

**Proposition B.1.** [83, Proposition 4.4] Let Assumption 3.5 holds for the nominal system. Assume that $x_d: R_+ \to R^n$ is AC, $\lambda_d(u_d, x_d)$ is bounded, and that the closed-loop system in (26) is well-posed. Then, the matrix inequality

$$ -M_0 + \Delta M_0 - \begin{pmatrix} \epsilon' P_0 & 0 \\ 0 & 0 \end{pmatrix} - \begin{pmatrix} P_0 \Lambda_1 P_0 & 0 \\ 0 & \Lambda_2 \end{pmatrix} \succcurlyeq 0 \quad (B.1) $$

where $M_0$ in (28) corresponds to the strict state passivity LMI with $\epsilon > 0$ and $\Delta M_0$ is defined in section 6, holds for $\epsilon > \epsilon' > 0$ if and only if:
1. $R \stackrel{\Delta}{=} (D_0 + F_0G_0) + (D_0 + F_0G_0)^\top + (\Delta D + \Delta FG_0) + (\Delta D + \Delta FG_0)^\top - \Lambda_2 \ge 0,$

2. $Q \stackrel{\Delta}{=} -(A_0 + E_0K_0 + \frac{\epsilon'}{2}I_n)^\top P_0 - P_0(A_0 + E_0K_0 + \frac{\epsilon'}{2}I_n) - (\Delta A + \Delta EK_0)^\top P_0 - P_0(\Delta A + \Delta EK_0) - P_0\Lambda_1P_0 \ge 0$, with $-(A_0 + E_0K_0)^\top P_0 - P_0(A_0 + E_0K_0) \ge \epsilon P_0$, $\epsilon > \epsilon' > 0$,

3. $\mathrm{Im}(S^\top) \subseteq \mathrm{Im}(R)$, with $S \stackrel{\Delta}{=} -P_0(B_0 + E_0G_0) + (C_0 + F_0K_0)^\top - P_0(\Delta B + \Delta EG_0) + (\Delta C + \Delta FK_0)^\top,$

4. $Q \ge SR^\dagger S^\top$.

The following Lemma gives sufficient conditions for item 3 in Proposition B.1 to be satisfied.

**Lemma B.2.** [83, Lemma 4.6] Let Assumption 3.5 hold for the nominal system. Assume that $\mathrm{Im}(\Delta S^\top) \subseteq \mathrm{Im}(S_0^\top)$, and that $\mathrm{Im}(R) = \mathrm{Im}(R_0)$, then $\mathrm{Im}(S^\top) \subseteq \mathrm{Im}(R)$.

The proposition below provides conditions on the uncertainty matrices to ensure that the inequality in item 4 holds when the nominal system is strictly state passive.

**Proposition B.3.** [83, Proposition 4.7] Let Assumption 3.5 hold for the quadruple $(A_0, B_0, C_0, D_0)$ and define the matrix $OT$ such that:

$$OT = S_0 R_0^\dagger \Delta S^\top - S_0 R_0^\dagger \Delta R_0^\dagger S_0^\top - S_0 R_0^\dagger \Delta R_0^\dagger \Delta S^\top + \Delta S R_0^\dagger S_0^\top + \Delta S R_0^\dagger \Delta S^\top - \Delta S R_0^\dagger \Delta R_0^\dagger S_0^\top - \Delta S R_0^\dagger \Delta R_0^\dagger \Delta S^\top.$$

Assume that: (i) $R \ge 0$, (ii) $\mathrm{Im}(\Delta R) \subseteq \mathrm{Im}(R_0)$, (iii) $\mathrm{rank}(R_0) = r$, $R_0 \ne 0$, $\sigma_{\max}(\Delta R) < \sigma_r(R_0)$, (iv) $\sigma_{\max}(OT) + \sigma_{\max}(O((R_0^\dagger \Delta R)^2)) + \sigma_{\max}(\Delta Q) < \epsilon\lambda_{\max}(P_0)$. Then, $Q \ge SR^\dagger S^\top$.

## Appendix C. Stick/Slip Behaviour of the Desired LCS in (13)

Recall the desired system in (13):

$$\begin{cases} \dot{x}_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} x_d(t) + \begin{pmatrix} 0 & 0 \\ \frac{1}{m} & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ \mu g \end{pmatrix} + \begin{pmatrix} 0 & 0 \\ -2\mu g & 0 \end{pmatrix} \lambda_d(t) \\ 0 \le \lambda_d(t) \perp w_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) \ge 0, \end{cases}$$

where $x_d = (q_d, \dot{q}_d)^\top = (x_{1d}, x_{2d})^\top$. According to the complementarity problem in (13), the multiplier $\lambda_{1d} \in [0, 1]$. Let us determine the value of the desired input $u_{1d}$ for the different stick/slip behaviour by considering the following analysis of the desired system in (13):

* If $\lambda_{1d} \in ]0, 1[$, $1[ \Rightarrow x_{2d} = 0 \Rightarrow \dot{x}_{2d} = 0 \Rightarrow u_{1d} \in ] -m\mu g, m\mu g[$. The system is in sticking mode.

* If $\lambda_{1d} = 0 \Rightarrow x_{2d} \le 0$. When $u_{1d} = -m\mu g$, then $x_{2d} = 0$ and the system is in sticking mode. But, when $u_{1d}(t) < -m\mu g \ \forall t$, then the solution $x_{2d} < 0$ (i.e., sliding mode) and it diverges.

* If $\lambda_{1d} = 1 \Rightarrow x_{2d} = \lambda_{2d} \ge 0$. When $u_{1d} = m\mu g$, then $x_{2d} = 0$ and the system behaves in sticking mode. When $u_{1d}$ is chosen such that $u_{1d}(t) > m\mu g \ \forall t$, then the solution $x_{2d} > 0$, it behaves in sliding mode and it is unbounded.

It is noteworthy that, in the cases when $\lambda_{1d} = 1$ or $\lambda_{1d} = 0$ (i.e., $u_{1d} = m\mu g$ or $u_{1d} = -m\mu g$ respectively), and the initial velocity of the mass $x_{2d}(0)$ is chosen such that it is in the direction of the external force $u_{1d}$ (i.e., $u_{1d} < 0$ and $x_{2d}(0) < 0$). Then the trajectory $x_{1d}$ diverges and it is sliding mode at a constant velocity.

If $u_{1d}(t)$ is chosen as a periodic, time-varying, and bounded signal in the form of $u_{1d} = A\cos(wt)$ where $|A| > m\mu g$, then the system will behave in both sticking and sliding modes or non-sticking mode depending on the amplitude A.
## Appendix D. Equilibrium Points of the Closed-loop Error System in (16)

In order to find the equilibrium points of the error dynamics, the LCS in (16) is written as follows:

$$ egin{cases} 0 = e_2^* \\ 0 = \frac{k_{11}}{m} e_1^* - 2\mu g(\lambda_1^*(t) - \lambda_{1d}(t)) \\ 0 \le \lambda_d(t) \perp w_d(t) = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) + \begin{pmatrix} 0 \\ 1 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 & -1 \\ 0 & 0 \end{pmatrix} x_d(t) \ge 0 \\ 0 \le \lambda^*(t) \perp w^*(t) = \begin{pmatrix} k_{21} & k_{22} - 1 \\ 0 & 0 \end{pmatrix} e^*(t) + \begin{pmatrix} g_{21} & 1 \\ -1 & 0 \end{pmatrix} \lambda^*(t) - \begin{pmatrix} g_{21} & 0 \\ 0 & 0 \end{pmatrix} \lambda_d(t) + \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix} u_d(t) \ge 0 \end{cases} \quad (D.1) $$

The problem in (D.1) is a mixed LCP with unknowns $e_1^*, e_2^*, \lambda_1^*, \lambda_2^*$. The equilibrium point of the error dynamics is given by $e_1^*(t) = \frac{2m\mu g}{k_{11}}(\lambda_1^*(t) - \lambda_{1d}(t))$. It is dependent on the complementarity variables $\lambda_{1d}$ and $\lambda_1^*$. In the above analysis, the LCP of the desired system is solved which gives the values of $\lambda_d$ and $x_d$. Now, it is required to solve the LCP of the closed-loop system (i.e., $0 \le \lambda^* \perp w^* \ge 0$) in order to find $\lambda^*$ and then calculate the equilibrium point $e_1^*$. Let us substitute $e^*$ in terms of $\lambda^*$ and $\lambda_d$, then the LCP to be solved is:

$$ 0 \le \begin{pmatrix} \lambda_1^*(t) \\ \lambda_2^*(t) \end{pmatrix} \perp \begin{pmatrix} w_1^*(t) \\ w_2^*(t) \end{pmatrix} = \begin{pmatrix} 2m\mu g \frac{k_{21}}{k_{11}} + g_{21} & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} \lambda_1^*(t) \\ \lambda_2^*(t) \end{pmatrix} + \begin{pmatrix} -2m\mu g \frac{k_{21}}{k_{11}} - g_{21} & 0 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \lambda_{1d}(t) \\ \lambda_{2d}(t) \end{pmatrix} + \begin{pmatrix} u_{2d}(t) - x_{2d}(t) \\ 1 \end{pmatrix} \ge 0 \quad (D.2) $$

It is important to note that the system can be at equilibrium, while the multipliers are time-varying. We can observe the new term $2m\mu g \frac{k_{21}}{k_{11}} + g_{21}$, that appears in the matrix of the multiplier $\lambda^*(t)$ in (D.2). This term plays a significant role in the subsequent analysis of the equilibrium point, so let us study the sign of this term. Let Assumption 3.5 holds for the closed-loop system in (15), then it is necessary to ensure that the following holds:

$$ -A_{cl}^{\top}P - PA_{cl} - \epsilon P \ge 0 \quad (D.3) $$

where $P = P^\top = \begin{pmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{pmatrix} > 0$ and $\epsilon > 0$. Hence, the matrix inequality in (D.3) is written as:

$$ \begin{pmatrix} -\frac{2k_{11}}{m}p_{12} & -p_{11} - \frac{k_{11}}{m}p_{22} - \frac{k_{12}}{m}p_{12} \\ -p_{11} - \frac{k_{11}}{m}p_{22} - \frac{k_{12}}{m}p_{12} & -2p_{12} - \frac{2k_{12}}{m}p_{22} \end{pmatrix} - \epsilon \begin{pmatrix} p_{11} & p_{12} \\ p_{12} & p_{22} \end{pmatrix} \ge 0 \quad (D.4) $$

According to the equality implied by passivity, $PB_{cl} = C_{cl}^\top$ since $D_{cl} + D_{cl}^\top = 0$, we have $-2p_{12}\mu g = k_{21}$ and $-2p_{22}\mu g = k_{22} - 1$. Let us substitute the value of $p_{12}$ in (D.4), then

$$ -\frac{2k_{11}}{m}p_{12} - \epsilon p_{11} \ge 0 \Leftrightarrow \frac{k_{11}k_{21}}{m\mu g} - \epsilon p_{11} \ge 0 $$

Given that $p_{11} > 0$, then $\frac{k_{11}k_{21}}{m\mu g}$ must be positive. This means that $k_{11}$ and $k_{21}$ have same signs. So, the term $2m\mu g \frac{k_{21}}{k_{11}}$ and according to Assumption 3.5 that imposes strict passivity on the closed-loop system, $\ge 0$. Hence, the term $2m\mu g \frac{k_{21}}{k_{11}} + g_{21}$ which is multiplied by $\lambda_1^*(t)$ is non-negative. In order to study the well-posedness of the closed-loop system's LCP in (D.2), let us consider different cases as above:

case 1: $\lambda_1^* \in ]0, 1[$. In this case, $w_1^* = (2m\mu g \frac{k_{21}}{k_{11}} + g_{21})(\lambda_1^* - \lambda_{1d}) = 0$. Hence, $\lambda_1^* = \lambda_{1d}$. The solutions are given by:
$$ \begin{cases} \lambda^* = (\lambda_{1d}, 0)^\top \\ w^* = (0, 1 - \lambda_{1d})^\top \end{cases} $$
case 2: $\lambda_1^* = 0$. In this case, $w_1^* = -(2m\mu g \frac{k_{21}}{k_{11}} + g_{21})\lambda_{1d} + u_{2d} - x_{2d} \ge 0$ which holds only if $u_{2d} - x_{2d} \ge 0$. If $u_{2d} - x_{2d} > 0$, then the solution of the desired LCP in (D.1) is:

$$ \begin{cases} \lambda_d = (0,0)^\top \\ w_d = (u_{2d} - x_{2d}, 1)^\top \end{cases} $$

If $u_{2d} - x_{2d} = 0$, then the solution is $\lambda_{1d} = 0$. Thus, in this case, $\lambda_1^* = \lambda_{1d} = 0$. The unique solution is given by:
$$ \begin{cases} \lambda^* = (0,0)^\top \\ w^* = (u_{2d} - x_{2d}, 1)^\top \end{cases} $$

case 3: $\lambda_1^* = 1$. In this case, $w_1^* = (2m\mu g \frac{k_{21}}{k_{11}} + g_{21})(1 - \lambda_{1d}) + \lambda_2^* + u_{2d} - x_{2d} = 0 \Rightarrow \lambda_2^* = -(2m\mu g \frac{k_{21}}{k_{11}} + g_{21})(1 - \lambda_{1d}) + x_{2d} - u_{2d} \ge 0$ which holds only if $x_{2d} - u_{2d} \ge 0$. If $x_{2d} - u_{2d} > 0$, then the solution of the desired system's LCP in (D.1) is:

$$ \begin{cases} \lambda_d = (1, x_{2d} - u_{2d})^\top \\ w_d = (0,0)^\top \end{cases} $$

If $x_{2d} - u_{2d} = 0$, then the inequality $-(2m\mu g \frac{k_{21}}{k_{11}} + g_{21})(1 - \lambda_{1d}) \ge 0$ holds if and only if $\lambda_{1d} = 0$. So, $\lambda_1^* = \lambda_{1d} = 1$. The solution is unique and it is given by:
$$ \begin{cases} \lambda^* = (1, x_{2d} - u_{2d})^\top \\ w^* = (0,0)^\top \end{cases} $$

Therefore, in all the possible cases of $\lambda_1^*$, it is observed that $\lambda_1^* = \lambda_{1d}$. It is noteworthy that in all cases, we have $\lambda^* = \lambda_d$ and $w^* = w_d$. Given that $e_1^* = \frac{2m\mu g}{k_{11}}(\lambda_1^*(t) - \lambda_{1d}(t))$ and $\lambda^*(t) = \lambda_d(t)$, then $e^* = (e_1^*, e_2^*)^\top = (0,0)^\top$. This is in accordance with the fact that Proposition 5.1 guarantees global exponential stability (hence uniqueness of the origin of the error dynamics).

## References

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

[2] R. I. Leine, D. H. van Campen, A. de Kraker, L. van den Steen, Stick-slip vibrations induced by alternate friction models, Nonlinear Dynamics 16 (1) (1998) 41–54.

[3] R. I. Leine, D. H. van Campen, B. Vrande, Bifurcations in nonlinear discontinuous systems, Nonlinear Dynamics 23 (2000) 105–164.

[4] R. I. Leine, D. H. van Campen, Discontinuous bifurcations of periodic solutions, Mathematical and Computer Modelling 36 (2002) 259–273.

[5] 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.

[6] 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.

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

[8] M. di Bernardo, C. Budd, A. Champneys, P. Kowalczyk, Piecewise-smooth Dynamical Systems. Theory and Applications, no. 163 in Applied Mathematical Sciences, Springer-Verlag, London, 2008.

[9] 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.

[10] B. A. Erickson, B. Birnir, D. Lavallee, Periodicity, chaos and localization in a Burridge–Knopoff model of an earthquake with rate-and-state friction, Geophysical Journal International 187 (2011) 178–198.

[11] M. Pascal, Dynamics of coupled oscillators excited by dry friction, Nonlinear Dynamics 3 (2008) 20–26.

[12] 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.

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

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

[15] A. C. J. Luo, B. C. Gegg, Periodic motions in a periodically forced oscillator moving on an oscillating belt with dry friction, Journal of Computational and Nonlinear Dynamics 1 (3) (2006) 212–220.

[16] K. Popp, P. Stelter, Stick-slip vibrations and chaos, Philosophical Transactions: Physical Sciences and Engineering 332 (1624) (1990) 89–105.

[17] J. Thomsen, Using fast vibrations to quench friction-induced oscillations, Journal of Sound and Vibration 228 (5) (1999) 1079–1102.

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

[19] A. J. Luo, M. Patel, Bifurcation and stability of periodic motions in a periodically forced oscillator with multiple discontinuities, Journal of Computational and Nonlinear Dynamics 4 (1) (2008) 011011.
[20] Y. B. Kim, S. T. Noah, Stability and bifurcation analysis of oscillators with piecewise-linear characteristics: A general approach, Journal of Applied Mechanics 58 (2) (1991) 545-553.

[21] 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.

[22] G. Csernák, G. Lískó, Asymmetric and chaotic responses of dry friction oscillators with different static and kinetic coefficients of friction, Meccanica 56 (06 2021).

[23] B. C. Gegg, A. C. Luo, S. C. Suh, Grazing bifurcations of a harmonically excited oscillator moving on a time-varying translation belt, Nonlinear Analysis: Real World Applications 9 (5) (2008) 2156-2174.

[24] K. Mehran, B. Zahawi, D. Giaouris, J. Wang, Stability analysis of sliding-grazing phenomenon in dry-friction oscillator using Takagi-Sugeno fuzzy approach, Journal of Computational and Nonlinear Dynamics 10 (6) (2015) 061010.

[25] Y. Li, Z. Feng, Bifurcation and chaos in friction-induced vibration, Communications in Nonlinear Science and Numerical Simulation 9 (6) (2004) 633-647.

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

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

[28] S. Chatterjee, Time-delayed feedback control of friction-induced instability, International Journal of Non-Linear Mechanics 42 (9) (2007) 1127-1143.

[29] 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.

[30] A. Saha, B. Bhattacharya, P. Wahi, A comparative study on the control of friction-driven oscillations by time-delayed feedback, Nonlinear Dynamics 60 (2010) 15-37.

[31] A. Saha, S. S. Pandey, B. Bhattacharya, P. Wahi, Analysis and control of friction-induced oscillations in a continuous system, Journal of Vibration and Control 18 (3) (2011) 467-480.

[32] A. Saha, P. Wahi, An analytical study of time-delayed control of friction-induced vibrations in a system with a dynamic friction model, International Journal of Non-Linear Mechanics 63 (2014) 60-70.

[33] Z. Iklodi, P. Piiroinen, O. Franco, X. Beudaert, Z. Dombovari, Stick-slip oscillations in the low feed linear motion of a grinding machine due to dry friction and backlash, Int. J. Non-Linear Mechanics 168 (2025) 104940.

[34] O. Zirn, Machine tool analysis modelling, simulation and control of machine tool manipulators, Habilitation thesis, ETH Zürich, https://doi.org/10.3929/ethz-a-005825192 (2008).

[35] O. Franco, X. Beudaert, K. Erkorkmaz, J. Munoa, Influence of guideway friction on the cutting point receptance in machine tools, CIRP Annals 71 (1) (2022) 361-364.

[36] O. Franco, X. Beudaert, K. Erkorkmaz, J. Munoa, Feed drive control and non-linear friction interaction effect on machining chatter stability prediction, Journal of Manufacturing and Materials Processing 8 (5) (2024).

[37] R. Leine, B. Brogliato, H. Nijmeijer, Periodic motion and bifurcations induced by the Painlevé paradox, European Journal of Mechanics - A/Solids 21 (5) (2002) 869-896.

[38] T. S. Amer, M. O. Shaker, H. A. Dahab, Analysis of the stick-slip behavior of coupled oscillators with dry friction, Journal of Low Frequency Noise, Vibration and Active Control 42 (3) (2023) 1251-1275.

[39] S. Khizgiyayev, Self-excited oscillations of a two-mass oscillator with dry "stick-slip" friction, Journal of Applied Mathematics and Mechanics 71 (6) (2007) 905-913.

[40] J. Awrejcewicz, P. Olejnik, Friction pair modeling by a 2-dof system: Numerical and experimental investigations, International Journal of Bifurcation and Chaos 15 (06) (2005) 1931-1944.

[41] J. Awrejcewicz, L. Dzyubak, C. Grebogi, Estimation of chaotic and regular (stick-slip and slip-slip) oscillations exhibited by coupled oscillators with dry friction, Nonlinear Dynamics 42 (2005) 383-394.

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

[43] 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.

[44] 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 Journal of Mechanics - A/Solids 96 (2022) 104737.

[45] G. Rill, T. Schaeffer, M. Schuderer, LuGre or not LuGre, Multibody System Dynamics 60 (2) (2024) 191-218.

[46] M. Schuderer, G. Rill, T. Schaeffer, C. Schulz, Friction modeling from a practical point of view, Multibody System Dynamics (2024) 1-18.

[47] S. Tafazoli, C. de Silva, P. Lawrence, Tracking control of an electrohydraulic manipulator in the presence of friction, IEEE Transactions on Control Systems Technology 6 (3) (1998) 401-411.

[48] K. C. Lin, Observer-based tension feedback control with friction and inertia compensation, IEEE Transactions on Control Systems Technology 11 (1) (2003) 109-118.

[49] J. C. L. Barreto S., A. G. S. Conceição, C. E. T. Dórea, L. Martinez, E. R. de Pieri, Design and implementation of model-predictive control with friction compensation on an omnidirectional mobile robot, IEEE/ASME Transactions on Mechatronics 19 (2) (2014) 467-476.

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

[51] M. Ruderman, On convergence analysis of feedback control with integral action and discontinuous relay perturbation, Communications in Nonlinear Science and Numerical Simulation 145 (2025) 108698.

[52] R. Oprea, C. Tudorache, M. Spiroiu, S. Arsene, C. Craciun, Comparative study of continuous versus discontinuous numerical models for railway vehicles suspensions with dry friction, Applied Sciences 15 (21) (2025).

[53] 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.
[54] E. Berger, Friction modeling for dynamic system simulation, Applied Mechanics Reviews 55 (6) (2002) 535-577.

[55] F. Al-Bender, J. Swevers, Characterization of friction force dynamics, IEEE Control Systems Magazine 28 (6) (2008) 64-81.

[56] B. Bekesi, M. Antali, G. Csernak, Phase portraits and bifurcations induced by static and dynamic friction models, Nonlinear Dynamics 113 (13) (2025) 15863-15899.

[57] A. Kayihan, F. J. Doyle, Friction compensation for a process control valve, Control Engineering Practice 8 (7) (2000) 799-812.

[58] B. Armstrong-Hélouvry, Control of Machines with Friction, 1st Edition, Springer New York, NY, 1991.

[59] P. Vedagarbha, D. Dawson, M. Feemster, Tracking control of mechanical systems in the presence of nonlinear dynamic friction effects, IEEE Transactions on Control Systems Technology 7 (4) (1999) 446-456.

[60] S. Southward, C. Radcliffe, C. MacCluer, Robust nonlinear stick-slip friction compensation, Journal of Dynamic Systems Measurement and Control 113 (12 1991).

[61] E. Papadopoulos, G. Chasparis, Analysis and model-based control of servomechanisms with friction, in: IEEE/RSJ International Conference on Intelligent Robots and Systems, Vol. 3, 2002, pp. 2109-2114 vol.3.

[62] N. Mallon, N. van de Wouw, D. Putra, H. Nijmeijer, Friction compensation in a controlled one-link robot using a reduced-order observer, IEEE Transactions on Control Systems Technology 14 (2) (2006) 374-383.

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

[64] 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.

[65] S. Huang, W. Liang, K. K. Tan, Intelligent friction compensation: A review, IEEE/ASME Transactions on Mechatronics 24 (4) (2019) 1763-1774.

[66] A. Bisoffi, R. Beerens, W. Heemels, H. Nijmeijer, N. van de Wouw, L. Zaccarian, To stick or to slip: A reset PID control perspective on positioning systems with friction, Annual Reviews in Control 49 (2020) 37-63.

[67] C. Makkar, G. Hu, W. G. Sawyer, W. E. Dixon, Lyapunov-based tracking control in the presence of uncertain nonlinear parameterizable friction, IEEE Transactions on Automatic Control 52 (10) (2007) 1988-1994.

[68] E. Panteley, R. Ortega, M. Gäfvert, An adaptive friction compensator for global tracking in robot manipulators, Systems & Control Letters 33 (5) (1998) 307-313.

[69] M. Ruderman, Tracking control of motor drives using feedforward friction observer, IEEE Transactions on Industrial Electronics 61 (7) (2014) 3727-3735.

[70] M. Ruderman, M. Iwasaki, Observer of nonlinear friction dynamics for motion control, IEEE Transactions on Industrial Electronics 62 (9) (2015) 5941-5949.

[71] D. Putra, H. Nijmeijer, N. van de Wouw, Analysis of undercompensation and overcompensation of friction in 1-dof mechanical systems, Automatica 43 (8) (2007) 1387-1394.

[72] A. Taware, G. Tao, N. Pradhan, C. Teolis, Friction compensation for a sandwich dynamic system, Automatica 39 (3) (2003) 481-488.

[73] J. C. A. de Bruin, A. Doris, N. van de Wouw, W. P. M. H. Heemels, H. Nijmeijer, Control of mechanical motion systems with non-collocation of actuation and friction: A Popov criterion approach for input-to-state stability and set-valued nonlinearities, Automatica 45 (2) (2009) 405-415.

[74] R. Beerens, A. Bisoffi, L. Zaccarian, W. P. M. H. Heemels, H. Nijmeijer, N. van de Wouw, Hybrid PID control for transient performance improvement of motion systems with friction, in: 2018 Annual American Control Conference (ACC), 2018, pp. 539-544.

[75] A. Bisoffi, M. Da Lio, A. Teel, L. Zaccarian, Global asymptotic stability of a PID control system with Coulomb friction, IEEE Transactions on Automatic Control 63 (8) (2018) 2654-2661.

[76] N. van de Wouw, R. I. Leine, Robust impulsive control of motion systems with uncertain friction, International Journal of Robust and Nonlinear Control 22 (4) (2011) 369-397.

[77] R. Beerens, A. Bisoffi, L. Zaccarian, W. Heemels, H. Nijmeijer, N. van de Wouw, Reset integral control for improved settling of PID-based motion systems with friction, Automatica 107 (2019) 483-492.

[78] R. Beerens, A. Bisoffi, L. Zaccarian, H. Nijmeijer, M. Heemels, N. van de Wouw, Reset PID design for motion systems with Stribeck friction, IEEE Transactions on Control Systems Technology 30 (1) (2022) 294-310.

[79] J. Alvarez, I. Orlov, L. Acho, An invariance principle for discontinuous dynamic systems with application to a Coulomb friction oscillator, ASME J. of Dynamic Syst., Meas., and Control 122 (4) (2000) 687-690.

[80] A. Saha, P. Wahi, Delayed feedback for controlling the nature of bifurcations in friction-induced vibrations, Journal of Sound and Vibration 330 (25) (2011) 6070-6087.

[81] K. Camlibel, Complementarity methods in the analysis of piecewise linear dynamical systems, Ph.D. thesis, Katholieke Universiteit Brabant, Tilburg, NL, iSBN 90 5668 079 X (2001).

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

[83] A. Younes, F. 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.

[84] 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.

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

[86] R. T. Rockafellar, R. J. B. Wets, Variational Analysis, Vol. 317 of Grundlehren der mathematischen Wissenschaften, Springer Verlag, Berlin Heidelberg, 2009, corrected Third Printing.

[87] R. W. Cottle, J. S. Pang, R. E. Stone, The Linear Complementarity Problem, Academic Press, 1992.

[88] B. Brogliato, R. Lozano, B. Maschke, O. Egeland, Dissipative Systems Analysis and Control, 3rd Edition, Communications and Control Eng., Springer Nature Switzerland AG, 2020, erratum/addendum at https://inria.hal.science/hal-03364924v3/document.

[89] B. Brogliato, L. Thibault, Existence and uniqueness of solutions for non-autonomous complementarity systems, Journal of Convex Analysis
17 (3-4) (2010) 961–990.

[90] M. K. Camlibel, L. Iannelli, F. Vasca, Passivity and complementarity, Mathematical Programming Ser. A 145 (2014) 531–563.

[91] M. K. Camlibel, W. P. M. H. Heemels, J. M. Schumacher, On linear passive complementarity systems, European Journal of Control 58 (2002) 220–237.

[92] S. Greenhalgh, V. Acary, B. Brogliato, On preserving dissipativity properties of linear complementarity dynamical systems with the theta-method, Numerische Mathematik 125 (4) (2013) 601–637.

[93] 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.

[94] Z. Li, H. Ouyang, Y. Gu, S. Martelli, S. Yang, H. Wei, W. Wang, R. Wei, Non-stationary friction-induced vibration with multiple contact points, Nonlinear Dynamics 111 (2023) 9889–9917.

[95] N. Han, Y. Gao, H. Zhang, Z. Li, Equilibrium states for a rotating mass-spring dry friction model on circular conveyor belt, Nonlinear Dynamics 113 (25) (2025) 34645–34676.

[96] J. Bastien, M. Schatzman, C. Lamarque, Study of an elastoplastic model with an infinite number of internal degrees of freedom, European Journal of Mechanics - A/Solids 21 (2) (2002) 199–222.

[97] 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.

[98] M. K. Camlibel, L. Iannelli, A. Tanwani, Convergence of proximal solutions for evolution inclusions with time-dependent maximal monotone operators, Mathematical Programming Ser. A 145 (2021) 1017–1059.

[99] Q. H. Pham, B. Brogliato, Correction to "Analysis of the implicit Euler time-discretization of a class of descriptor-variable linear cone complementarity systems", J. Convex Analysis 29/2 (2022) 481–517, Journal of Convex Analysis 31 (3) (2024) 1035–1037.

[100] E. Sontag, Mathematical Control Theory: Deterministic Finite Dimensional Systems., Texts in Applied Mathematics, Springer, New York, 2013.

[101] W. P. M. H. Heemels, V. Sessa, F. Vasca, M. K. Camlibel, Computation of periodic solutions in maximal monotone dynamical systems with guaranteed consistency, Nonlinear Analysis: Hybrid Systems 24 (2017) 100–114.

[102] 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.

[103] V. Acary, O. Bonnefon, M. Brémond, O. Huber, F. Pérignon, S. Sinclair, An introduction to Siconos, Technical Report RT-0340, INRIA, https://inria.hal.science/inria-00162911v3/file/RR-0340-v3.pdf (2019).

[104] V. Acary, B. Brogliato, Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics, Vol. 35, Springer Verlag, 2008.

[105] 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.

[106] MOSEK ApS, MOSEK Optimizer API for Python (2023). URL https://docs.mosek.com/10.1/pythonapi.pdf

[107] D. S. Bernstein, Matrix Mathematics. Theory, Facts and Formulas, 2nd Edition, Princeton Univ. Press, 2009.

[108] R. A. Horn, C. R. Johnson, Matrix Analysis, 2nd Edition, Cambridge University Press, 2013.

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

[110] M. Corless, G. Leitmann, Continuous state feedback guaranteeing uniform ultimate boundedness for uncertain dynamic systems, IEEE Transactions on Automatic Control 26 (5) (1981) 1139–1144.

[111] L. Márton, B. Lantos, Control of mechanical systems with stribeck friction and backlash, Systems and Control Letters 58 (2) (2009) 141–147.

[112] B. Brogliato, Dissipative dynamical systems with set-valued feedback loops: Well-posed set-valued Lur'e dynamical systems, IEEE Control Systems 42 (3) (2022) 93–114.

[113] A. Blumentals, B. Brogliato, F. Bertails-Descoubes, The contact problem in Lagrangian systems subject to bilateral and unilateral constraints, with or without sliding Coulomb's friction: a tutorial, Multibody System Dynamics 38 (1) (2016) 43–76.

[114] 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.

[115] C. Scherer, S. Weiland, Linear Matrix Inequalities in Control, 2015, https://www.imng.uni-stuttgart.de/mst/files/LectureNotes.pdf.

[116] F. Fachinei, J. S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Volume I, Operations Research and Financial Engineering, Springer, New York, 2003.