**Remark 5.** The above control strategy has been proposed in [21, 31]. For instance, equation (94b) is exactly [21, Equation (9)].

### 3.2.1 Stability Analysis and Sliding Mode

Let us briefly recall the main stability properties of (94). It is proved in [21, Proposition 2] that {$\tilde{s}_k$}$_{k \in \mathbb{N}}$ solution of (94b) converges to zero in a finite number $k_{\min} < +\infty$ of steps (cf. Proposition 8 in "Sidebar VII" as well), while {$s_k$}$_{k \in \mathbb{N}}$ solution of (94a) is bounded. During the sliding mode ($\tilde{s}_k = 0$ for $k \ge k_{\min}$) we have $s_{k+1} = \tilde{C}\tilde{B}\delta_k$, since $s_k + u_{sv}(x_k) = 0$ is equivalent to

$$u_{sv}(x_k) = -s_k = -\tilde{C}\tilde{B}\delta_{k-1}. \quad (98)$$

As a fundamental result, the control (21) is approximated by (98): in discrete time, the implicit method allows us to design an input which compensates for the perturbation with a one-step delay, in a similar way to (39). Once again the set-valued controller does not depend on the gain $\tilde{K}$ during the sliding phase.

### 3.2.2 Control-Input Calculation

Using Proposition 6 in "Sidebar VI", equation (94), and mapping inversion (see "Sidebar III"), we can see that

$$\tilde{K} (\tilde{K}^{-1} s_k - (-\tilde{K}^{-1} u_{sv}(x_k)) \in \mathbf{N}_{\mathcal{B}_{\infty}} (-\tilde{K} u_{sv}(x_k)) \quad (99)$$

is equivalent to $-u_{sv}(x_k) = \tilde{K}(I_d + \tilde{K}^{-1} \mathbf{N}_{\mathcal{B}_{\infty}})^{-1}(\tilde{K}^{-1} s_k)$, which in turn is equivalent to

$$-u_{sv}(x_k) = \tilde{K} \operatorname{Proj}_{\tilde{K}} (\mathcal{B}_{\infty}; \tilde{K}^{-1} s_k). \quad (100)$$

Again, the projection can be computed by solving a quadratic program or, similarly to (89), conveniently choose $\tilde{K} = \gamma I$ so that

$$(u_{sv}(x_k))_i = -\gamma \min \left\{ \frac{|(s_k)_i|}{\gamma}, 1 \right\} \operatorname{sgn}((s_k)_i), \quad (101)$$

see "Sidebar VI" for $\mathcal{B}_{\infty}$, (45a) and the associated piece of code for the computation of the projection in a general setting. The controller in (101) is in the minimum-operator formalism [61, 48, 50, 62, 63].

**Remark 6.** Dead-beat controller. The linear controllers in (17) and in (92) have different structures. Suppose that $u_{\text{lin}}(x_k)$ is computed as $u_{\text{lin}}(x_k) = -(\tilde{C}\tilde{B})^{-1}\tilde{C}\tilde{A}x_k$, which would be a naive implementation of (17). Then,

$$\tilde{C}x_{k+1} = \tilde{C}\tilde{A}x_k + \tilde{C}\tilde{B}(u_{\text{lin}}(x_k) + \tilde{\delta}_k) + v_k = v_k + \tilde{C}\tilde{B}\tilde{\delta}_k, \quad (102)$$

which is different from (93). Applying the above $u_{sv}(x_k)$ yields $\tilde{s}_{k+1} \in -\tilde{K} \operatorname{Sgn}(\tilde{s}_{k+1})$, which is equivalent to $\tilde{s}'_{k+1} = (I_d + \partial f)^{-1}(0)$. Since the operator $(I_d + \partial f)^{-1}$ is single-valued and its graph contains (0, 0), it follows that $\tilde{s}_{k+1} = 0$. Therefore, this algorithm is a one-step

dead-beat controller. The implicit algorithms introduced in [30, 31, 21] are not of the dead-beat type. It is noteworthy that a major discrepancy exists between (92) and the dead-beat input in the framework of discretization of continuous systems. Consider for instance the exact ZOH discretization with constant sampling time $h > 0$, where $\tilde{A} = e^{Ah}$, $\tilde{B} = e^{Ah} \int_0^h e^{-A\tau} B d\tau$. Then, as $h \to 0$ we have $(\tilde{C}\tilde{B})^{-1}\tilde{C}(I_n - \tilde{A}) = \mathcal{O}(1)$ while $(\tilde{C}\tilde{B})^{-1}\tilde{C}\tilde{A} = \mathcal{O}(h^{-1})$. Clearly, the dead-beat input grows unbounded as $h \to 0$ and cannot be a good candidate for convergence towards its continuous-time counterpart. This shows that such dead-beat controller is well-suited to purely discrete-time systems, but not for discretized continuous-time systems.

## 3.3 Recapitulation

Allow us to summarize some of the connections between sliding-mode control, passivity, and optimization.

### 3.3.1 Link with Optimization

It is interesting to note that operators of the form $(I_d + M)^{-1}$ are ubiquitous in the above discrete-time algorithms: they are the well-known resolvent mappings, very close to Yosida approximations (see (23) and (26) in "Sidebar V"). They have been exhibited in (10), (11), (73), (86), (96), and (100). They are a direct consequence of the presence of a maximal monotone mapping $M$ on the right-hand side of the differential inclusions in (22), (36). These operators are resolvents of mappings of the form $M = \partial f$ for some proper convex lower semi-continuous function $f$. This shows the strong connection between the implicit discretization and Optimization (see also (86)). Thus, a common feature between model-predictive control (MPC) and set-valued sliding-mode implicit discretization is that an optimization problem has to be solved to compute the controller. This is formalized in "Sidebar VIII", where the implicitly discretized scheme is reinterpreted as a one-step MPC. The explicit discretization (that is, choosing $u_{sv}(x_k) \in -\tilde{K} \operatorname{Sgn}(s_k)$) yields operators of the form $(I_d - \partial f)$, which are neither monotone nor possess interesting properties, neither in terms of stability nor robustness.

The above analysis is closely related to the analysis of proximal-point algorithms (Definition 12 in "Sidebar VII"). However, a peculiarity of discrete sliding-mode control is the notion of a robust proximal-point algorithm, formulated in (57) in "Sidebar VII" and shown to converge in a finite number of steps, given sufficiently large gains. It is worth emphasizing that (67) has the form $u_{sv}(x_k) \in -\partial \|\tilde{K} \tilde{s}_{k+1}\|_p$ with $p = 1$. Just as with the continuous-time case, other norms can be used (or other proper LSC convex functions for that matter). For arbitrary $p$, it suffices to replace $\mathcal{B}_{\infty}$ in (86) by $\mathcal{B}_{p*}$. For
Proximal-point algorithms and proximal operators are fundamental tools in Optimization. The big names and founding fathers of this field are Jean Jacques Moreau [S30, S31, S32], Kôsacu Yosida [S33], R. Tyrell Rockafellar [S18, 54], and Hedy Attouch [S34]. Let us recall that this article deals with finite-dimensional systems only. As seen in (11), (13), (74), (86), (96), and (97), proximal operators and proximal-point algorithms are important tools for the implicit discretization of subgradient systems. While the convergence properties of proximal-point algorithms have been thoroughly studied in the literature, the robustness of such algorithms in the perturbed scenarios that naturally arise in Control is apparently unknown in Optimization. Given a maximal monotone operator $M: \mathbb{R}^n \to \mathbb{R}^n$ (see "Sidebar V"), its set of zeros is denoted by zero $M = M^{-1}(0) = \{\xi \in \mathbb{R}^n | 0 \in M(\xi)\}$.

**Definition 12.** Proximal-Point Algorithm [S28, p. 345]. Assume that $M: \mathbb{R}^n \to \mathbb{R}^n$ is a maximal monotone operator satisfying zero $M \neq \emptyset$ and let $x_0 \in \mathbb{R}^n$ and $\gamma > 0$. The proximal-point algorithm associated with $M$ is defined as

$$x_{k+1} = \mathcal{J}_{\gamma M}(x_k), \quad (49)$$

for $k \in \mathbb{N}$, where $\mathcal{J}_{\gamma M}(\cdot)$ is the resolvent of $M$.

Theorem 1 in [S35] establishes the convergence of (49) towards zero $M$. Consider a proper, lower semicontinuous, and convex function $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$. The operator

$$\mathrm{Prox}_f = (I_d + \partial f)^{-1} \quad (50)$$

is called the proximal mapping [S28, Definition 12.23]. By definition, when $M = \partial f$ the algorithm (49) is equivalent to any of the following iterations:

$$x_{k+1} = \mathrm{Prox}_{\gamma f}(x_k) \quad (51)$$

$$x_{k+1} - x_k \in -\gamma \partial f(x_{k+1}) \quad (52)$$

$$x_{k+1} = \arg \min_{\xi \in \mathbb{R}^n} \left\{ f(\xi) + \frac{1}{2\gamma} \|\xi - x_k\|^2 \right\}. \quad (53)$$

The form (51) gives the algorithm its name. It can be readily seen from (50) that

$$\text{zero } \mathrm{Prox}_f = \partial f(0). \quad (54)$$

Also remark that, if $f = \Psi_C$ (see Definition 10) with $C \subseteq \mathbb{R}^n$ a closed convex nonempty set, then,

$$\mathrm{Prox}_f = \mathrm{Proj}(K; \cdot) = (I_d + \mathbf{N}_K)^{-1}. \quad (55)$$

When implicitly or semi-implicitly discretized, several higher-order differentiators and controllers involve the resolvent of the subdifferential of the proper convex LSC function $f: \mathbb{R} \to \mathbb{R}$, $x \mapsto \frac{1}{2}\gamma x^2 + \Psi_{[-1,1]}(x)$, $\gamma > 0$, that is, $\partial f(x) = \gamma x + \mathbf{N}_{[-1,1]}(x)$, see "Sidebar VI". This boils down to calculating the inverse of the set-valued function $\partial f' : x \mapsto (1 + \gamma)x + \mathbf{N}_{[-1,1]}(x)$, which is a saturation,

$$\mathcal{J}_{\partial f'}(y) = \begin{cases} -1 & \text{if } y \le -1 - \gamma \\ \frac{1}{1+\gamma}y & \text{if } y \in [-1-\gamma, 1+\gamma] \\ 1 & \text{if } y \ge 1+\gamma \end{cases}. \quad (56)$$

**Remark 7.** It can be readily seen from (52) that the proximal-point algorithm corresponds to the backward (implicit) Euler discretization of the gradient system $\dot{x}(t) \in -\frac{\gamma}{h} \partial f(x(t))$, which is well-known to be more efficient and have better approximation properties than its forward (explicit) Euler counterpart [55].

Implicit discrete-time sliding-mode controllers and differentiators will frequently yield perturbed iterations of the form

$$x_{k+1} = \mathcal{J}_{\gamma M}(x_k + w_k). \quad (57)$$

The next Proposition is the discrete-time counterpart of Proposition 2.
Proximal-point algorithms and proximal operators are fundamental tools in Optimization. The big names and founding fathers of this field are Jean Jacques Moreau [S30, S31, S32], Kôsacu Yosida [S33], R. Tyrell Rockafellar [S18, 54], and Hedy Attouch [S34]. Let us recall that this article deals with finite-dimensional systems only. As seen in (11), (13), (74), (86), (96), and (97), proximal operators and proximal-point algorithms are important tools for the implicit discretization of subgradient systems. While the convergence properties of proximal-point algorithms have been thoroughly studied in the literature, the robustness of such algorithms in the perturbed scenarios that naturally arise in Control is apparently unknown in Optimization. Given a maximal monotone operator $M: \mathbb{R}^n \to \mathbb{R}^n$ (see "Sidebar V"), its set of zeros is denoted by zero $M = M^{-1}(0) = \{\xi \in \mathbb{R}^n | 0 \in M(\xi)\}$.

**Definition 12.** Proximal-Point Algorithm [S28, p. 345]. Assume that $M: \mathbb{R}^n \to \mathbb{R}^n$ is a maximal monotone operator satisfying zero $M \neq \emptyset$ and let $x_0 \in \mathbb{R}^n$ and $\gamma > 0$. The proximal-point algorithm associated with $M$ is defined as

$$x_{k+1} = \mathcal{J}_{\gamma M}(x_k), \quad (49)$$

for $k \in \mathbb{N}$, where $\mathcal{J}_{\gamma M}(\cdot)$ is the resolvent of $M$.

Theorem 1 in [S35] establishes the convergence of (49) towards zero $M$. Consider a proper, lower semicontinuous, and convex function $f: \mathbb{R}^n \to \mathbb{R} \cup \{+\infty\}$. The operator

$$\mathrm{Prox}_f = (I_d + \partial f)^{-1} \quad (50)$$

is called the proximal mapping [S28, Definition 12.23]. By definition, when $M = \partial f$ the algorithm (49) is equivalent to any of the following iterations:

$$x_{k+1} = \mathrm{Prox}_{\gamma f}(x_k) \quad (51)$$

$$x_{k+1} - x_k \in -\gamma \partial f(x_{k+1}) \quad (52)$$

$$x_{k+1} = \arg \min_{\xi \in \mathbb{R}^n} \left\{ f(\xi) + \frac{1}{2\gamma} |\xi - x_k|^2 \right\}. \quad (53)$$

The form (51) gives the algorithm its name. It can be readily seen from (50) that

$$\text{zero } \mathrm{Prox}_f = \partial f(0). \quad (54)$$

Also remark that, if $f = \Psi_C$ (see Definition 10) with $C \subseteq \mathbb{R}^n$ a closed convex nonempty set, then,

$$\mathrm{Prox}_f = \mathrm{Proj}(K; \cdot) = (I_d + \mathbf{N}_K)^{-1}. \quad (55)$$

When implicitly or semi-implicitly discretized, several higher-order differentiators and controllers involve the resolvent of the subdifferential of the proper convex LSC function $f: \mathbb{R} \to \mathbb{R}$, $x \mapsto \frac{1}{2}\gamma x^2 + \Psi_{[-1,1]}(x)$, $\gamma > 0$, that is, $\partial f(x) = \gamma x + \mathbf{N}_{[-1,1]}(x)$, see "Sidebar VI". This boils down to calculating the inverse of the set-valued function $\partial f' : x \mapsto (1 + \gamma)x + \mathbf{N}_{[-1,1]}(x)$, which is a saturation,

$$\mathcal{J}_{\partial f'}(y) = \begin{cases} -1 & \text{if } y \le -1 - \gamma \\ \frac{1}{1+\gamma}y & \text{if } y \in [-1-\gamma, 1+\gamma] \\ 1 & \text{if } y \ge 1+\gamma \end{cases}. \quad (56)$$

**Remark 7.** It can be readily seen from (52) that the proximal-point algorithm corresponds to the backward (implicit) Euler discretization of the gradient system $\dot{x}(t) \in -\frac{\gamma}{h}\partial f(x(t))$, which is well-known to be more efficient and have better approximation properties than its forward (explicit) Euler counterpart [55].

Implicit discrete-time sliding-mode controllers and differentiators will frequently yield perturbed iterations of the form

$$x_{k+1} = \mathcal{J}_{\gamma M}(x_k + w_k). \quad (57)$$

The next Proposition is the discrete-time counterpart of Proposition 2.
**Proposition 8.** Robust proximal-point algorithm. Consider the perturbed proximal-point algorithm (57) with $\gamma > 0$ and $\mathbf{M}: \mathbb{R}^n \to \mathbb{R}^n$ a maximal monotone operator such that

$$ (\gamma \mathbf{M})^{-1} \left( \bigcup_{k \in \mathbb{N}} w_k \right) = \{0\}. \qquad (58) $$

The origin $x = 0$ is globally asymptotically stable. If, moreover,

$$ (\gamma \mathbf{M})^{-1} \left( \bigcup_{k \in \mathbb{N}} w_k + \varepsilon \mathcal{B} \right) = \{0\} \qquad (59) $$

for some $\varepsilon > 0$, then it is globally finite-time stable.

*Proof.* From the definition of Yosida approximation (26), it follows that $x_{k+1} = x_k + w_k - \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k)$. Thus

$$ \begin{aligned} \|x_{k+1}\|_2^2 &= x_{k+1}^\top (x_k + w_k - \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k)) \\ &\le \frac{1}{2} \|x_{k+1}\|_2^2 + \frac{1}{2} \|x_k\|_2^2 \\ &\quad - \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k)^\top (\gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k), \end{aligned} \qquad (60) $$

where we used the fact that $\langle x, y \rangle \le \frac{1}{2} \|x\|_2^2 + \frac{1}{2} \|y\|_2^2$. Setting the Lyapunov function candidate as $V_{k+1} = \frac{1}{2} \|x_{k+1}\|_2^2$, it is obtained

$$ V_{k+1} - V_k \le -\mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k)^\top (\gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k). \qquad (61) $$

In what follows, it is shown that the right-hand side of (61) is negative definite. To that end, we point out that assumption (59) implies that for any $b \in \mathcal{B}$ and any $k \in \mathbb{N}$, $\varepsilon b - w_k \in \gamma \mathbf{M}(0)$. Thus, it follows from (28) and the maximal monotonicity of $\gamma \mathbf{M}$, that for all $k \in \mathbb{N}$

$$ \langle \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k), \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - \varepsilon b - w_k \rangle \ge 0. \qquad (62) $$

After rearranging terms, we can see that

$$ \langle \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k), \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k \rangle \ge \varepsilon \sup_{b \in \mathcal{B}} \langle b, \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k) \rangle. \qquad (63) $$

That is,

$$ \begin{aligned} \langle \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k), \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k \rangle &\ge \varepsilon \| \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k) \| \\ &\ge \varepsilon \| x_{k+1} \|. \end{aligned} \qquad (64) $$

Therefore, $V_{k+1} - V_k \le -\varepsilon \|x_{k+1}\|$ and the global asymptotic stability of the origin follows. Now, for the finite-time convergence, notice that (57) is equivalent to

$$ x_k + w_k \in (I_d + \gamma \mathbf{M})(x_{k+1}), \qquad (65) $$

whereas (59) and the asymptotic stability of the origin imply that there is $k^* \in \mathbb{N}$ and $\varepsilon_2 > 0$ such that for all $k \ge k^*$ and any $b \in \mathcal{B}$,

$$ w_k + x_k + \varepsilon_2 b \in (I_d + \gamma \mathbf{M})(0). \qquad (66) $$

Thus, (65)-(66) and the maximal monotonicity of $I_d + \gamma \mathbf{M}$ imply that for all $b \in \mathcal{B}$

$$ \varepsilon_2 \langle b, x_{k+1} \rangle \ge 0. \qquad (67) $$

It thus follows that $x_{k+1} = 0$ for all $k \ge k^*$. $\square$

<figure><img src="image_7.png" /></figure>
**Proposition 8.** Robust proximal-point algorithm. Consider the perturbed proximal-point algorithm (57) with $\gamma > 0$ and $\mathbf{M}: \mathbb{R}^n \to \mathbb{R}^n$ a maximal monotone operator such that

$$ (\gamma \mathbf{M})^{-1} \left( \bigcup_{k \in \mathbb{N}} w_k \right) = \{0\}. \qquad (58) $$

The origin $x = 0$ is globally asymptotically stable. If, moreover,

$$ (\gamma \mathbf{M})^{-1} \left( \bigcup_{k \in \mathbb{N}} w_k + \varepsilon \mathcal{B} \right) = \{0\} \qquad (59) $$

for some $\varepsilon > 0$, then it is globally finite-time stable.

*Proof.* From the definition of Yosida approximation (26), it follows that $x_{k+1} = x_k + w_k - \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k)$. Thus

$$ \begin{aligned} \|x_{k+1}\|_2^2 &= x_{k+1}^\top (x_k + w_k - \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k)) \\ &\le \frac{1}{2} \|x_{k+1}\|_2^2 + \frac{1}{2} \|x_k\|_2^2 \\ &\quad - \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k)^\top (\gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k), \end{aligned} \qquad (60) $$

where we used the fact that $\langle x, y \rangle \le \frac{1}{2} \|x\|_2^2 + \frac{1}{2} \|y\|_2^2$. Setting the Lyapunov function candidate as $V_{k+1} = \frac{1}{2} \|x_{k+1}\|_2^2$, it is obtained

$$ V_{k+1} - V_k \le - \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k)^\top (\gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k). \qquad (61) $$

In what follows, it is shown that the right-hand side of (61) is negative definite. To that end, we point out that assumption (59) implies that for any $b \in \mathcal{B}$ and any $k \in \mathbb{N}$, $\varepsilon b - w_k \in \gamma \mathbf{M}(0)$. Thus, it follows from (28) and the maximal monotonicity of $\gamma \mathbf{M}$, that for all $k \in \mathbb{N}$

$$ \langle \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k), \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - \varepsilon b - w_k \rangle \ge 0. \qquad (62) $$

After rearranging terms, we can see that

$$ \langle \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k), \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k \rangle \ge \varepsilon \sup_{b \in \mathcal{B}} \langle b, \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k) \rangle. \qquad (63) $$

That is,

$$ \begin{aligned} \langle \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k), \gamma \mathcal{Y}_{\gamma \mathbf{M}}(x_k + w_k) - w_k \rangle &\ge \varepsilon \| \mathcal{J}_{\gamma \mathbf{M}}(x_k + w_k) \| \\ &\ge \varepsilon \| x_{k+1} \|. \end{aligned} \qquad (64) $$

Therefore, $V_{k+1} - V_k \le -\varepsilon \|x_{k+1}\|$ and the global asymptotic stability of the origin follows. Now, for the finite-time convergence, notice that (57) is equivalent to

$$ x_k + w_k \in (I_d + \gamma \mathbf{M})(x_{k+1}), \qquad (65) $$

whereas (59) and the asymptotic stability of the origin imply that there is $k^* \in \mathbb{N}$ and $\varepsilon_2 > 0$ such that for all $k \ge k^*$ and any $b \in \mathcal{B}$,

$$ w_k + x_k + \varepsilon_2 b \in (I_d + \gamma \mathbf{M})(0). \qquad (66) $$

Thus, (65)-(66) and the maximal monotonicity of $I_d + \gamma \mathbf{M}$ imply that for all $b \in \mathcal{B}$

$$ \varepsilon_2 \langle b, x_{k+1} \rangle \ge 0. \qquad (67) $$

It thus follows that $x_{k+1} = 0$ for all $k \ge k^*$.

□
## References

[S30] J.J. Moreau, Fonctions convexes duales et points proximaux dans un espace hilbertien, *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 255:2897-2899, 1962.

[S31] J.J. Moreau, Propriétés des applications “prox”, *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 256:1069-1071, 1963.

[S32] J.J. Moreau, Proximité et dualité dans un espace hilbertien, *Bulletin de la Soc. Math. France*, tome 93:273-299, 1965.

[S33] K. Yosida, Functional Analysis, Springer Berlin, Heidelberg, 1974.

[S34] H. Attouch, G. Buttazzo and G. Michaille, Variational Analysis in Sobolev and BV Spaces, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2014.

[S35] E.K. Ryu, W. Yin, Large-scale convex optimization via monotone operators, Cambridge University Press, 2023.

example, the implicit discretization of the unit control is

$$ -u_{sv}(x_k) = \tilde{K} \operatorname{Proj}_M (B_2; \tilde{K}^{-1}\xi(x_k)). \quad (103) $$

Again, it is possible to choose $\tilde{K}$ such that $M = \gamma^2 I$ and use (43) to compute an explicit form for the projection.

### 3.3.2 Link with Passivity

Both control design methodologies described above follow the celebrated passivity plus zero-state detectability or observability results [68, 75, 67] (see [72] for the discrete-time version). Regarding the direct passivity-based design, Theorem 4 shows that a triplet $(V(x_k), s_{k+1}, v_k + \tilde{\delta}_k)$ can be chosen, defining a passifiable system. Hence it is necessary that the system has relative degree zero and is of minimum phase. The system is actually rendered passive with $u_{\text{lin}}(x_k)$, so that

$$ \begin{aligned} x_{k+1} &= \tilde{A}x_k + \tilde{B}u_{\text{lin}}(x_k) + \tilde{B}(v_k + \tilde{\delta}_k) \\ s_{k+1} &= \tilde{C}x_{k+1} \end{aligned} \quad (104) $$

happens to be zero-state detectable (ZSD) [72, Definition 2.2]. In this setting, the discretization with an implicit scheme certainly turns out to be necessary for assuring passivity (a fact that was not noticed in previous works on the topic). The next step is in fact an extension of [72, Theorem 2.6], using a set-valued static output feedback controller. The principle is the same for the discrete-time classical controller design although, in most of the classical literature, the ZSD property is assumed a priori rather than shown like it was done in Theorem 4. A notion that is closely related to feedback passivity is that of control Lyapunov functions. The connection between sliding variables and control-Lyapunov functions has been explicitly recognized and exploited, for example, in [76,

71], although it is worth mentioning that the idea already appears in its early stages in [77, 78].

### 3.3.3 Closed-Loop Properties

It is well known that, in continuous time, a sliding-mode controller gives rise to a differential inclusion. For a properly designed controller, closed-loop solutions converge to a sliding surface $\{\xi \in \mathbb{R}^n | \sigma(\xi) = 0\}$ in finite time. Moreover, once it has reached the sliding surface, a solution stays on the surface and becomes independent from any perturbations matched by the control. The latter is a remarkable property that can only be explained by the set-valued nature of the closed-loop vector field. A suitably discretized controller should approximately replicate these properties. This is the case, for example, of the implicitly-defined control law (67), since the condition $\sigma(x_k, u_k) = 0$ is attained in finite time and maintained there after, regardless of any bounded matched perturbation.

The discretization $u_{sv}(x_k) \in -\tilde{K} \operatorname{Sgn}(\tilde{K}\tilde{s}_k)$ produces digital chattering and loses the robustness and stability properties of its continuous-time counterpart, no matter how small $h > 0$ is. The implicit discretization, on the other hand, maintains these properties, even for relatively large $h$ (a fact verified experimentally). At the same time, implicit discretizations are consistent, in the sense that discrete-time solutions and inputs converge to their continuous-time counterparts as $h \downarrow 0$ [43, 21]. When written explicitly, the implicit (or, backward-Euler) discretization takes the form of a projection function of the state on a set, which in some simple instances can be computed using saturation functions. Mark, however, that when using the implicit discretization, correct parameters for the saturation functions are automatically provided by
## References

[S30] J.J. Moreau, Fonctions convexes duales et points proximaux dans un espace hilbertien, *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 255:2897-2899, 1962.

[S31] J.J. Moreau, Propriétés des applications “prox”, *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, 256:1069-1071, 1963.

[S32] J.J. Moreau, Proximité et dualité dans un espace hilbertien, *Bulletin de la Soc. Math. France*, tome 93:273-299, 1965.

[S33] K. Yosida, Functional Analysis, Springer Berlin, Heidelberg, 1974.

[S34] H. Attouch, G. Buttazzo and G. Michaille, Variational Analysis in Sobolev and BV Spaces, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2014.

[S35] E.K. Ryu, W. Yin, Large-scale convex optimization via monotone operators, Cambridge University Press, 2023.

example, the implicit discretization of the unit control is

$$-u_{sv}(x_k) = \tilde{K} \operatorname{Proj}_M (B_2; \tilde{K}^{-1}\xi(x_k)). \quad (103)$$

Again, it is possible to choose $\tilde{K}$ such that $M = \gamma^2 I$ and use (43) to compute an explicit form for the projection.

### 3.3.2 Link with Passivity

Both control design methodologies described above follow the celebrated passivity plus zero-state detectability or observability results [68, 75, 67] (see [72] for the discrete-time version). Regarding the direct passivity-based design, Theorem 4 shows that a triplet $(V(x_k), s_{k+1}, v_k + \delta_k)$ can be chosen, defining a passifiable system. Hence it is necessary that the system has relative degree zero and is of minimum phase. The system is actually rendered passive with $u_{\text{lin}}(x_k)$, so that

$$\begin{align}
x_{k+1} &= \tilde{A}x_k + \tilde{B}u_{\text{lin}}(x_k) + \tilde{B}(v_k + \tilde{\delta}_k) \tag{104} \\
s_{k+1} &= \tilde{C}x_{k+1} \nonumber
\end{align}$$

happens to be zero-state detectable (ZSD) [72, Definition 2.2]. In this setting, the discretization with an implicit scheme certainly turns out to be necessary for assuring passivity (a fact that was not noticed in previous works on the topic). The next step is in fact an extension of [72, Theorem 2.6], using a set-valued static output feedback controller. The principle is the same for the discrete-time classical controller design although, in most of the classical literature, the ZSD property is assumed a priori rather than shown like it was done in Theorem 4. A notion that is closely related to feedback passivity is that of control Lyapunov functions. The connection between sliding variables and control-Lyapunov functions has been explicitly recognized and exploited, for example, in [76,

71], although it is worth mentioning that the idea already appears in its early stages in [77, 78].

### 3.3.3 Closed-Loop Properties

It is well known that, in continuous time, a sliding-mode controller gives rise to a differential inclusion. For a properly designed controller, closed-loop solutions converge to a sliding surface $\{\xi \in \mathbb{R}^n | \sigma(\xi) = 0\}$ in finite time. Moreover, once it has reached the sliding surface, a solution stays on the surface and becomes independent from any perturbations matched by the control. The latter is a remarkable property that can only be explained by the set-valued nature of the closed-loop vector field. A suitably discretized controller should approximately replicate these properties. This is the case, for example, of the implicitly-defined control law (67), since the condition $\sigma(x_k, u_k) = 0$ is attained in finite time and maintained there after, regardless of any bounded matched perturbation.

The discretization $u_{sv}(x_k) \in -\tilde{K} \operatorname{Sgn}(\tilde{K}\tilde{s}_k)$ produces digital chattering and loses the robustness and stability properties of its continuous-time counterpart, no matter how small $h > 0$ is. The implicit discretization, on the other hand, maintains these properties, even for relatively large $h$ (a fact verified experimentally). At the same time, implicit discretizations are consistent, in the sense that discrete-time solutions and inputs converge to their continuous-time counterparts as $h \downarrow 0$ [43, 21]. When written explicitly, the implicit (or, backward-Euler) discretization takes the form of a projection function of the state on a set, which in some simple instances can be computed using saturation functions. Mark, however, that when using the implicit discretization, correct parameters for the saturation functions are automatically provided by
Discretization with a backward-Euler scheme turns out to be necessary for
assuring passivity.

the method: no trial-and-error tuning is required.

The properties described in the previous two para-
graphs can be formulated using asymptotic notation.
Consider, for example, the nominal case. If the state is
inside a neighborhood of order h² of the sliding surface,
||s<sub>k</sub>|| = O(h²), then at the following step we only have
||s<sub>k+1</sub>|| = O(h) and the state exits the neighborhood [21,
Lemma 10]. If, on the other hand, the discretization is im-
plicit, the sliding variable remains inside a neighborhood
of the same order, ||s<sub>k+1</sub>|| = O(h²) [21, Lemma 11].

The implicit discretization method results in a con-
trolled system that is also robust with respect to uncer-
tainty in the system parameters [42, 43]. With the appro-
priate choice of the maximal monotone operators, it is also
possible to obtain robustness with respect to unmatched
external perturbations as well [79, 44]. The explicit com-
putation of the control law can be carried out in several
ways, depending on the complexity of the problem. In
some cases, it is possible to obtain an analytic expression
(for example, as in (100)). In more complex scenarios, it
is always possible to resort to the numerical machinery
described in section "6".

3.3.4 On the Nature of Sliding Mode Control

Contrarily to a widely spread idea, sliding-mode control
is not intrinsically related to discontinuous and infinitely-
fast switching inputs (bang-bang-like controllers). The
implicit discretization shows this fact, as it allows us to
correctly approximate not only the closed-loop system's
output *s* but, most importantly, the set-valued input.
Two fundamental notions are prominent: the selection of
a set-valued map (Definition 3 in "Sidebar III") and the
calculation of the controller with optimization tools, as
detailed below. See Figure 8 below, for a typical illustra-
tion of discontinuous and selections of set-valued control
inputs.

4 HIGHER-ORDER SET-VALUED CONTROLLERS AND DIFFERENTIATORS

The goal of this section is to show how the foregoing de-
velopments extend to “modern” sliding-mode algorithms
dedicated to control, state observation, or differentiation.
Most importantly, higher-order schemes, introduced by
Arie Levant and co-authors [9, 8, 10], have witnessed
considerable attention in the past twenty years. Roughly
speaking, first-order schemes permit the designer to guar-
antee that a sliding variable (denoted *s* in section “2”)

vanishes in finite-time. Then the origin is attained asymp-
totically while the state trajectories evolve on the slid-
ing surface. Higher-order schemes of order *r* (or, order-*r*-
schemes) guarantee the finite-time convergence of the slid-
ing variable *s* and its derivatives up to the (*r* − 1)th order,
to zero. This is named an *r*-sliding surface and the *r*-sliding mode is defined by {*x* ∈ Rⁿ | *s*(*x*) = *š*(*x*) = · · · = *s*(*r*−¹)(*x*) = 0}. From a practical point of view, higher-order sliding-mode method offers several advantages: 1) in the extreme case *r* = *n* > 1, the origin is attained in finite time, a feature not possible with first-order schemes, 2) by adding integrators at the plant’s input, it is possible to steer the plant using arbitrarily smooth controls, 3) given a signal with bounded (*r* + 1)th time derivative, it is possible to recover the first *r* derivatives.

The main principle used in the discretization schemes presented in previous sections is structure preservation, that is, the preservation of passivity and maximal monotonicity. Structure preservation presents several challenges in the case of *r*-sliding mode control. First, a necessary condition for establishing an *r*-sliding mode on *s* is that the sliding output has relative degree *r*. Thus, for *r* > 1 the plant cannot be feedback equivalent to a passive system. Second, most higher-order sliding-mode control algorithms are not defined by maximal monotone differential inclusions. Finally, the main underlying structure is homogeneity [80], a property that is not meaningful in discrete time. However, we will show in this section that maximal monotonicity still plays an essential role in the backward-Euler discretization of such schemes.

The first part of this section is dedicated to the
super-twisting algorithm [10], which is certainly the
most widespread modern sliding-mode scheme. Then,
other higher-order differentiators and controllers are pre-
sented. The section ends with a short paragraph about
homogeneity-based schemes and on how proximal-point
algorithms appear in this setting.

4.1 Super-Twisting Algorithm

The super-twisting is an order-2 algorithm. It has proven
useful for control, differentiation and state observation.
In the following, all three application fields are presented.

4.1.1 Differentiator

Consider a signal $f: \mathbb{R} \rightarrow \mathbb{R}$ and suppose that we wish
to differentiate it in real time. A common solution in
Discretization with a backward-Euler scheme turns out to be necessary for
assuring passivity.

the method: no trial-and-error tuning is required.

The properties described in the previous two para-
graphs can be formulated using asymptotic notation.
Consider, for example, the nominal case. If the state is
inside a neighborhood of order h² of the sliding surface,
||s<sub>k</sub>|| = O(h²), then at the following step we only have
||s<sub>k+1</sub>|| = O(h) and the state exits the neighborhood [21,
Lemma 10]. If, on the other hand, the discretization is im-
plicit, the sliding variable remains inside a neighborhood
of the same order, ||s<sub>k+1</sub>|| = O(h²) [21, Lemma 11].

The implicit discretization method results in a con-
trolled system that is also robust with respect to uncer-
tainty in the system parameters [42, 43]. With the appro-
priate choice of the maximal monotone operators, it is also
possible to obtain robustness with respect to unmatched
external perturbations as well [79, 44]. The explicit com-
putation of the control law can be carried out in several
ways, depending on the complexity of the problem. In
some cases, it is possible to obtain an analytic expression
(for example, as in (100)). In more complex scenarios, it
is always possible to resort to the numerical machinery
described in section "6".

3.3.4 On the Nature of Sliding Mode Control

Contrarily to a widely spread idea, sliding-mode control
is not intrinsically related to discontinuous and infinitely-
fast switching inputs (bang-bang-like controllers). The
implicit discretization shows this fact, as it allows us to
correctly approximate not only the closed-loop system's
output *s* but, most importantly, the set-valued input.
Two fundamental notions are prominent: the selection of
a set-valued map (Definition 3 in "Sidebar III") and the
calculation of the controller with optimization tools, as
detailed below. See Figure 8 below, for a typical illustra-
tion of discontinuous and selections of set-valued control
inputs.

4 HIGHER-ORDER SET-VALUED CONTROLLERS AND DIFFERENTIATORS

The goal of this section is to show how the foregoing de-
velopments extend to “modern” sliding-mode algorithms
dedicated to control, state observation, or differentiation.
Most importantly, higher-order schemes, introduced by
Arie Levant and co-authors [9, 8, 10], have witnessed
considerable attention in the past twenty years. Roughly
speaking, first-order schemes permit the designer to guar-
antee that a sliding variable (denoted *s* in section “2”)

vanishes in finite-time. Then the origin is attained asymp-
totically while the state trajectories evolve on the slid-
ing surface. Higher-order schemes of order *r* (or, order-*r*-
schemes) guarantee the finite-time convergence of the slid-
ing variable *s* and its derivatives up to the (*r* − 1)th order,
to zero. This is named an *r*-sliding surface and the *r*-sliding mode is defined by {*x* ∈ Rⁿ | *s*(*x*) = *š*(*x*) = · · · = *s*(*r*−¹)(*x*) = 0}. From a practical point of view, higher-order sliding-mode method offers several advantages: 1) in the extreme case *r* = *n* > 1, the origin is attained in finite time, a feature not possible with first-order schemes, 2) by adding integrators at the plant’s input, it is possible to steer the plant using arbitrarily smooth controls, 3) given a signal with bounded (*r* + 1)th time derivative, it is possible to recover the first *r* derivatives.

The main principle used in the discretization schemes presented in previous sections is structure preservation, that is, the preservation of passivity and maximal monotonicity. Structure preservation presents several challenges in the case of *r*-sliding mode control. First, a necessary condition for establishing an *r*-sliding mode on *s* is that the sliding output has relative degree *r*. Thus, for *r* > 1 the plant cannot be feedback equivalent to a passive system. Second, most higher-order sliding-mode control algorithms are not defined by maximal monotone differential inclusions. Finally, the main underlying structure is homogeneity [80], a property that is not meaningful in discrete time. However, we will show in this section that maximal monotonicity still plays an essential role in the backward-Euler discretization of such schemes.

The first part of this section is dedicated to the
super-twisting algorithm [10], which is certainly the
most widespread modern sliding-mode scheme. Then,
other higher-order differentiators and controllers are pre-
sented. The section ends with a short paragraph about
homogeneity-based schemes and on how proximal-point
algorithms appear in this setting.

4.1 Super-Twisting Algorithm

The super-twisting is an order-2 algorithm. It has proven
useful for control, differentiation and state observation.
In the following, all three application fields are presented.

4.1.1 Differentiator

Consider a signal $f: \mathbb{R} \rightarrow \mathbb{R}$ and suppose that we wish
to differentiate it in real time. A common solution in
## Sidebar VIII: Implicit Sliding-Mode and Model Predictive Control (MPC)

Implicit discrete-time sliding-mode control and differentiation require the solution of an optimization problem at each time-step. Model predictive control requires the solution of an optimization problem on longer time intervals. Besides this feature, there exists a closer and more formal link between implicit sliding modes and model predictive control. Consider an optimization problem over an horizon event of one step,

$$ \min_{(s_{k+1}, u_k)} f(s_{k+1} - h\delta_k) + \frac{h}{2} \|u_k\|^2, \\ \text{such that } s_{k+1} = s_k + hu_k + h\delta_k \qquad (68) $$

where the function $f: \mathbb{R}^m \to \mathbb{R}$ is proper, convex and LSC. In this problem, at time $t_k = kh$, $s_k$ is treated as a known parameter, whereas $\delta_k$ is an unknown parameter and the target is to find the optimal controller value $u_k$ and the optimal state value $s_{k+1}$ that minimize the given cost function. The Lagrangian associated with (68) is given as,

$$ L(s_{k+1}, u_k, \lambda_{k+1}) = f(s_{k+1} - h\delta_k) + \frac{h}{2} \|u_k\|^2 - \lambda_{k+1}^\top (s_{k+1} - s_k - hu_k - h\delta_k). \quad (69) $$

The Karush-Kuhn-Tucker conditions associated with (68) are

i) $0 \in \partial_{s_{k+1}} L(s_{k+1}^*, u_k^*)$,

ii) $0 = \nabla_{u_k} L(s_{k+1}^*, u_k^*)$, and

iii) $s_{k+1}^* = s_k + hu_k^* + h\delta_k$.

From KKT-i) one sees that $0 \in \partial f(s_{k+1}^* - h\delta_k) - \lambda_{k+1}^* \Leftrightarrow \lambda_{k+1}^* \in \partial f(s_{k+1}^* - h\delta_k)$, whereas KKT-ii) yields $0 = hu_k^* + h\lambda_{k+1}^* \Leftrightarrow \lambda_{k+1}^* = -u_k^*$. Hence, the combination of both relations together with condition KKT-iii) yields $h\lambda_{k+1}^* = -hu_k^* = s_k + h\delta_k - s_{k+1}^* \in h\partial f(s_{k+1}^* - h\delta_k)$. Consequently,

$$ s_k \in (I + h\partial f)(s_{k+1}^* - h\delta_k) \Leftrightarrow s_{k+1}^* - h\delta_k = \text{Prox}_{h f}(s_k), \qquad (70) $$

and the substitution of (70) back into KKT-iii) leads to

$$ s_{k+1}^* - h\delta_k = s_k + hu_k^* = \text{Prox}_{h f}(s_k), \qquad (71) $$

equivalently,

$$ u_k^* = -\frac{1}{h}(s_k - \text{Prox}_{h f}(s_k)). \qquad (72) $$

Setting the dummy variable $\tilde{s}_{k+1}$ as $\tilde{s}_{k+1} = s_{k+1}^* - h\delta_k$ yields

$$ \begin{aligned} s_{k+1}^* &= \tilde{s}_{k+1} + h\delta_k \\ \tilde{s}_{k+1} &= s_k + hu_k^*, \\ -u_k^* &\in \partial f(\tilde{s}_{k+1}) \end{aligned} \qquad (73) $$

which are the same equations as the implicit discretization with $f(\cdot) = \gamma \|\cdot\|$, $\gamma > 0$.
## Sidebar VIII: Implicit Sliding-Mode and Model Predictive Control (MPC)

Implicit discrete-time sliding-mode control and differentiation require the solution of an optimization problem at each time-step. Model predictive control requires the solution of an optimization problem on longer time intervals. Besides this feature, there exists a closer and more formal link between implicit sliding modes and model predictive control. Consider an optimization problem over an horizon event of one step,

$$ \min_{(s_{k+1}, u_k)} f(s_{k+1} - h\delta_k) + \frac{h}{2} \|u_k\|^2, \\ \text{such that } s_{k+1} = s_k + hu_k + h\delta_k \qquad (68) $$

where the function $f: \mathbb{R}^m \to \mathbb{R}$ is proper, convex and LSC. In this problem, at time $t_k = kh$, $s_k$ is treated as a known parameter, whereas $\delta_k$ is an unknown parameter and the target is to find the optimal controller value $u_k$ and the optimal state value $s_{k+1}$ that minimize the given cost function. The Lagrangian associated with (68) is given as,

$$ L(s_{k+1}, u_k, \lambda_{k+1}) = f(s_{k+1} - h\delta_k) + \frac{h}{2} \|u_k\|^2 - \lambda_{k+1}^\top (s_{k+1} - s_k - hu_k - h\delta_k). \quad (69) $$

The Karush-Kuhn-Tucker conditions associated with (68) are

i) $0 \in \partial_{s_{k+1}} L(s_{k+1}^*, u_k^*)$,

ii) $0 = \nabla_{u_k} L(s_{k+1}^*, u_k^*)$, and

iii) $s_{k+1}^* = s_k + hu_k^* + h\delta_k$.

From KKT-i) one sees that $0 \in \partial f(s_{k+1}^* - h\delta_k) - \lambda_{k+1}^* \Leftrightarrow \lambda_{k+1}^* \in \partial f(s_{k+1}^* - h\delta_k)$, whereas KKT-ii) yields $0 = hu_k^* + h\lambda_{k+1}^* \Leftrightarrow \lambda_{k+1}^* = -u_k^*$. Hence, the combination of both relations together with condition KKT-iii) yields $h\lambda_{k+1}^* = -hu_k^* = s_k + h\delta_k - s_{k+1}^* \in h\partial f(s_{k+1}^* - h\delta_k)$. Consequently,

$$ s_k \in (I + h\partial f)(s_{k+1}^* - h\delta_k) \Leftrightarrow s_{k+1}^* - h\delta_k = \text{Prox}_{h f}(s_k), \qquad (70) $$

and the substitution of (70) back into KKT-iii) leads to

$$ s_{k+1}^* - h\delta_k = s_k + hu_k^* = \text{Prox}_{h f}(s_k), \qquad (71) $$

equivalently,

$$ u_k^* = -\frac{1}{h}(s_k - \text{Prox}_{h f}(s_k)). \qquad (72) $$

Setting the dummy variable $\tilde{s}_{k+1}$ as $\tilde{s}_{k+1} = s_{k+1}^* - h\delta_k$ yields

$$ \begin{aligned} s_{k+1}^* &= \tilde{s}_{k+1} + h\delta_k \\ \tilde{s}_{k+1} &= s_k + hu_k^*, \\ -u_k^* &\in \partial f(\tilde{s}_{k+1}) \end{aligned} \qquad (73) $$

which are the same equations as the implicit discretization with $f(\cdot) = \gamma \|\cdot\|$, $\gamma > 0$.
Projected dynamical systems (PDS) make a class of nonsmooth systems which, as seen below, is closely related to sliding-mode systems. They were introduced in [S36, S37]. Recently, they have witnessed a strong interest from the automatic-control scientific community, see, for example, [S38, S39]. PDS may be expressed under different equivalent forms [S40, S41] (see definitions 8 and 9 in “Sidebar VI”). Some examples of such equivalent dynamics are:

$$ \dot{x} \in -f(x) - g(t) - \mathbf{N}_c(x) \quad \text{with} \quad \dot{x} = -f(x) - g(t) - \operatorname{Proj}[\mathbf{N}_c(x); -f(x) - g(t)] \quad (74a) $$

$$ \dot{x} \in -f(x) - g(t) - \mathbf{N}_{\mathbf{T}_c(x)}(\dot{x}) \quad (74b) $$

$$ \dot{x} = (I_d + \mathbf{N}_{\mathbf{T}_c(x)})^{-1}(-f(x) - g(t)) \quad (74c) $$

$$ \dot{x} = \operatorname{Proj}[\mathbf{T}_c(x); -f(x) - g(t)] \quad (74d) $$

see [S40, S41] for details. The second condition in (74a) means that solutions are slow (with $-\dot{x}(t)$ of minimal norm inside $f(x) + g(t) + \mathbf{N}_c(x)$). The resolvent of the normal cone to the tangent cone appears in (74c). Use was made of Fact 5 in “Sidebar VI”, (50), and the comment just after (54) in “Sidebar VII” to obtain (74d). The projection is a particular instance of Prox, see (55), showing that PDS belong to the class of proximal systems.

### 3.3.5 Discrete-time PDS

Let us now discretize (74a) as

$$ x_{k+1} \in x_k - h(f(x_k) - g_k) - \mathbf{N}_c(x_{k+1}) \quad (75) $$

$h > 0$. This is rewritten equivalently as

$$ x_{k+1} = (I_d + \mathbf{N}_c)^{-1}(x_k - h(f(x_k) + g_k)) \quad (76) $$

and as

$$ x_{k+1} = \operatorname{Proj}[C; x_k - h(f(x_k) + g_k)], \quad (77) $$

where the same tools as above have again been used to pass from (74c) to (74d). It is seen once again that resolvents are ubiquitous in PDS. The system (77) furnishes a convenient way to obtain solutions whenever the projection onto the set $C$ can be computed. This is the case when $C$ is nonempty convex polyhedral (48). Indeed, in this case the normal cone can be expressed as $\mathbf{N}_C = \{\xi \in \mathbb{R}^n | \xi = G^\top \lambda, 0 \le \lambda \perp -G\xi + b \ge 0\}$. In other words, the normal cone is generated by the outwards normals to the active constraints at $x$. Thus, using (75), the PDS is equivalently rewritten as a complementarity system

$$ \begin{gathered} x_{k+1} = x_k - h(f(x_k) - g_k) - G^\top \lambda_{k+1} \quad (78) \\ 0 \le \lambda_{k+1} \perp w_{k+1} = -Gx_{k+1} + b \ge 0. \end{gathered} $$

The discrete-time PDS in (78) is a mixed linear complementarity problem [S42] with unknowns $x_{k+1}$ and $\lambda_{k+1}$. It gives the linear complementarity problem (LCP)

$$ 0 \le \lambda_{k+1} \perp w_{k+1} = -G(x_k - h(f(x_k) - g_k)) + GG^\top \lambda_{k+1} + b \ge 0, \quad (79) $$

Obviously, $GG^\top > 0$ if and only if, $G$ has row-rank $m$ (which implies that $m \le n$), otherwise $GG^\top \ge 0$ and it is also a copositive-plus matrix. Using [S42, Theorem 3.8.6], if $b \in (\operatorname{SOL}(GG^\top, 0))^*$, where $\operatorname{SOL}(M, 0)$ denotes the set of solutions to $0 \le \lambda \perp M\lambda \ge 0$ and $(\operatorname{SOL}(M, 0))^*$ is its dual cone, then the LCP in (79) has a solution. Assuming a Slater condition (that is, that there exists $\lambda$ such that $GG^\top\lambda - G(x_k - h(f(x_k) - g_k)) + b > 0$) then by [S42, Theorem 5.3.9] the LCP in (79) has a solution. Also if $\lambda_1$ and $\lambda_2$ are any two solutions, then $GG^\top(\lambda_1 - \lambda_2) = 0$ [S42, Theorem 3.1.7]. Therefore $\lambda_1 - \lambda_2 \in \operatorname{Ker}(GG^\top) = \operatorname{Ker}(G^\top)$ [65, Theorem 3.5.3]. It is easily deduced that, if

$$ x_{k+1}^i = x_k^i + h(f(x_k^i) - g_k) - G^\top \lambda_{k+1}^i, \quad (80) $$

$i = 1, 2$, then $x_k^1 = x_k^2$ implies $x_{k+1}^1 = x_{k+1}^2$, $k \ge 0$. As expected, the projection onto $C$ is unique. A quite interesting fact is that efficient algorithms can be used to solve the LCP in (79) [S42, 59], at each time-step $k \ge 1$.
Projected dynamical systems (PDS) make a class of nonsmooth systems which, as seen below, is closely related to sliding-mode systems. They were introduced in [S36, S37]. Recently, they have witnessed a strong interest from the automatic-control scientific community, see, for example, [S38, S39]. PDS may be expressed under different equivalent forms [S40, S41] (see definitions 8 and 9 in “Sidebar VI”). Some examples of such equivalent dynamics are:

$$ \dot{x} \in -f(x) - g(t) - \mathbf{N}_c(x) \quad \text{with} \quad \dot{x} = -f(x) - g(t) - \operatorname{Proj}[\mathbf{N}_c(x); -f(x) - g(t)] \quad (74a) $$

$$ \dot{x} \in -f(x) - g(t) - \mathbf{N}_{\mathbf{T}_c(x)}(\dot{x}) \quad (74b) $$

$$ \dot{x} = (I_d + \mathbf{N}_{\mathbf{T}_c(x)})^{-1}(-f(x) - g(t)) \quad (74c) $$

$$ \dot{x} = \operatorname{Proj}[\mathbf{T}_c(x); -f(x) - g(t)] \quad (74d) $$

see [S40, S41] for details. The second condition in (74a) means that solutions are slow (with $-\dot{x}(t)$ of minimal norm inside $f(x) + g(t) + \mathbf{N}_c(x)$). The resolvent of the normal cone to the tangent cone appears in (74c). Use was made of Fact 5 in “Sidebar VI”, (50), and the comment just after (54) in “Sidebar VII” to obtain (74d). The projection is a particular instance of Prox, see (55), showing that PDS belong to the class of proximal systems.

### 3.3.5 Discrete-time PDS

Let us now discretize (74a) as

$$ x_{k+1} \in x_k - h(f(x_k) - g_k) - \mathbf{N}_c(x_{k+1}) \quad (75) $$

$h > 0$. This is rewritten equivalently as

$$ x_{k+1} = (I_d + \mathbf{N}_c)^{-1}(x_k - h(f(x_k) + g_k)) \quad (76) $$

and as

$$ x_{k+1} = \operatorname{Proj}[C; x_k - h(f(x_k) + g_k)], \quad (77) $$

where the same tools as above have again been used to pass from (74c) to (74d). It is seen once again that resolvents are ubiquitous in PDS. The system (77) furnishes a convenient way to obtain solutions whenever the projection onto the set $C$ can be computed. This is the case when $C$ is nonempty convex polyhedral (48). Indeed, in this case the normal cone can be expressed as $\mathbf{N}_c = \{\xi \in \mathbb{R}^n | \xi = G^\top \lambda, 0 \le \lambda \perp -G\xi + b \ge 0\}$. In other words, the normal cone is generated by the outwards normals to the active constraints at $x$. Thus, using (75), the PDS is equivalently rewritten as a complementarity system

$$ \begin{gathered} x_{k+1} = x_k - h(f(x_k) - g_k) - G^\top \lambda_{k+1} \quad (78) \\ 0 \le \lambda_{k+1} \perp w_{k+1} = -Gx_{k+1} + b \ge 0. \end{gathered} $$

The discrete-time PDS in (78) is a mixed linear complementarity problem [S42] with unknowns $x_{k+1}$ and $\lambda_{k+1}$. It gives the linear complementarity problem (LCP)

$$ 0 \le \lambda_{k+1} \perp w_{k+1} = -G(x_k - h(f(x_k) - g_k)) + GG^\top \lambda_{k+1} + b \ge 0, \quad (79) $$

Obviously, $GG^\top > 0$ if and only if, $G$ has row-rank $m$ (which implies that $m \le n$), otherwise $GG^\top \ge 0$ and it is also a copositive-plus matrix. Using [S42, Theorem 3.8.6], if $b \in (\operatorname{SOL}(GG^\top, 0))^*$, where $\operatorname{SOL}(M, 0)$ denotes the set of solutions to $0 \le \lambda \perp M\lambda \ge 0$ and $(\operatorname{SOL}(M, 0))^*$ is its dual cone, then the LCP in (79) has a solution. Assuming a Slater condition (that is, that there exists $\lambda$ such that $GG^\top \lambda - G(x_k - h(f(x_k) - g_k)) + b > 0$) then by [S42, Theorem 5.3.9] the LCP in (79) has a solution. Also if $\lambda_1$ and $\lambda_2$ are any two solutions, then $GG^\top(\lambda_1 - \lambda_2) = 0$ [S42, Theorem 3.1.7]. Therefore $\lambda_1 - \lambda_2 \in \operatorname{Ker}(GG^\top) = \operatorname{Ker}(G^\top)$ [65, Theorem 3.5.3]. It is easily deduced that, if

$$ x_{k+1}^i = x_k^i + h(f(x_k^i) - g_k) - G^\top \lambda_{k+1}^i, \quad (80) $$

$i = 1, 2$, then $x_k^1 = x_k^2$ implies $x_{k+1}^1 = x_{k+1}^2$, $k \ge 0$. As expected, the projection onto $C$ is unique. A quite interesting fact is that efficient algorithms can be used to solve the LCP in (79) [S42, 59], at each time-step $k \ge 1$.
### 3.3.6 Sliding Modes in PDS

Some kind of sliding modes can occur in PDS when the single-valued vector field is such that the trajectory evolves persistently on the boundary of $C$. In this case there exists an element of $\lambda(x, t) \in -N_C(x)$ (a selection, see Definition 3 in “Sidebar III”) which “compensates” for the part of $-f(x) - g(t)$ which tends to “push” the trajectory outside $C$ (in a way quite similar to what occurs in Contact Mechanics where the interaction contact force balances the other forces which otherwise would make the system leave the admissible domain). It is noteworthy that the set-valued controlled systems in [79, Equation (6)] can be transformed into a PDS after the classical state space transformation (see [S41, section 3.4]). It guarantees robust output tracking.

## References

[S36] C. Henry, Differential equations with discontinuous right-hand side for planning procedures, *J. Econom. Theory*, 4:545-551, 1972.

[S37] C. Henry, An existence theorem for a class of differential equations with multivalued right-hand side, *J. Math. Anal. Appl.*, 41:179-186, 1973.

[S38] D. Gadjov, L. Pavel, A passivity-based approach to Nash equilibrium seeking over networks, *IEEE Transactions on Automatic Control*, 64(3):1077-1092, March 2019.

[S39] L. Pavel, Dissipativity theory in game theory: On the role of dissipativity and passivity in Nash equilibrium seeking, *IEEE Control Systems Magazine*, 42(3):150-164, June 2022.

[S40] B. Brogliato, A. Daniliidis, C. Lemaréchal, V. Acary, On the equivalence between complementarity systems, projected systems and differential inclusions, *Systems Control Letters*, 55:45-51, 2006

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

[S42] R.W. Cottle, J.S. Pang, R.E. Stone, The Linear Complementarity Problem, Academic Press, Computer Science and Science Computing, 1992.
## Sidebar (cont.)

### 3.3.6 Sliding Modes in PDS

Some kind of sliding modes can occur in PDS when the single-valued vector field is such that the trajectory evolves persistently on the boundary of $C$. In this case there exists an element of $\lambda(x, t) \in -N_C(x)$ (a selection, see Definition 3 in “Sidebar III”) which “compensates” for the part of $-f(x) - g(t)$ which tends to “push” the trajectory outside $C$ (in a way quite similar to what occurs in Contact Mechanics where the interaction contact force balances the other forces which otherwise would make the system leave the admissible domain). It is noteworthy that the set-valued controlled systems in [79, Equation (6)] can be transformed into a PDS after the classical state space transformation (see [S41, section 3.4]). It guarantees robust output tracking.

## References

[S36] C. Henry, Differential equations with discontinuous right-hand side for planning procedures, *J. Econom. Theory*, 4:545-551, 1972.

[S37] C. Henry, An existence theorem for a class of differential equations with multivalued right-hand side, *J. Math. Anal. Appl.*, 41:179-186, 1973.

[S38] D. Gadjov, L. Pavel, A passivity-based approach to Nash equilibrium seeking over networks, *IEEE Transactions on Automatic Control*, 64(3):1077-1092, March 2019.

[S39] L. Pavel, Dissipativity theory in game theory: On the role of dissipativity and passivity in Nash equilibrium seeking, *IEEE Control Systems Magazine*, 42(3):150-164, June 2022.

[S40] B. Brogliato, A. Daniliidis, C. Lemaréchal, V. Acary, On the equivalence between complementarity systems, projected systems and differential inclusions, *Systems Control Letters*, 55:45-51, 2006

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

[S42] R.W. Cottle, J.S. Pang, R.E. Stone, The Linear Complementarity Problem, Academic Press, Computer Science and Science Computing, 1992.
practice is to replace $\dot{f}$ by its *dirty derivative*, $\hat{g}$, which is obtained by passing $f$ through a linear filter of the form

$$ 
\begin{aligned}
\dot{z}(t) &= -\gamma(z(t) + \gamma f(t)) \\
\hat{g}(t) &= z(t) + \gamma f(t)
\end{aligned}
\qquad (105) $$

with $\gamma \gg 1$. It can be easily shown that, as $\gamma$ approaches $+\infty$, the transfer function of the filter converges to the transfer function of a differentiator. This high-gain strategy benefits from simplicity, but it has two important drawbacks: noise amplification and phase lag. If $\dot{f}$ is Lipschitz continuous with Lipschitz constant $L > 0$, then a super-twisting differentiation algorithm can be used. The advantages of such scheme are: 1) exact differentiation with finite gains in the absence of noise, 2) limited noise amplification, 3) no phase lag.

We will make use of the operator $[w]^\alpha := |w|^\alpha \operatorname{sgn}(w)$, $\alpha \ge 0$, $w \in \mathbb{R}$. Notice that $[w]^\alpha = \partial\zeta_\alpha(w)$ with

$$ \zeta_\alpha(w) = \frac{1}{\alpha+1} |w|^{\alpha+1}. \qquad (106) $$

The convexity of $\zeta_\alpha$ immediately guarantees the maximal monotonicity of $[\cdot]^\alpha$, see Theorem 3 in "Sidebar V". Notice also that $[\cdot]^\alpha$ is single-valued for $\alpha > 0$, but multivalued at the origin when $\alpha = 0$. The super-twisting differentiation algorithm [6] reads as:

$$ 
\begin{aligned}
\dot{z}_0(t) &= z_1(t) - \gamma_0 L^{\frac{1}{2}} [s_0(t)]^{\frac{1}{2}}, \\
\dot{z}_1(t) &\in -\gamma_1 L \operatorname{sgn}(s_0(t))
\end{aligned}
\qquad (107) $$

where $s_0(t) := z_0(t) - f(t)$ and $\gamma_0, \gamma_1$ are positive gains. The differentiator properties are best appreciated under the coordinate transformation $\sigma_0(t) := \frac{1}{L}(z_0 - f(t))$ and $\sigma_1(t) := \frac{1}{L}(z_1(t) - \dot{f}(t))$. In these coordinates, system (107) is rewritten in recursive form as [81]

$$ 
\begin{aligned}
\dot{\sigma}_0(t) &= \sigma_1(t) - \gamma_0 [\sigma_0(t)]^{\frac{1}{2}} \\
\dot{\sigma}_1(t) &\in -\gamma_1 \operatorname{sgn}(\sigma_0(t)) - \ddot{f}(t)
\end{aligned}
\qquad (108) $$

If $\gamma_0$ and $\gamma_1$ are appropriately chosen, then the origin $(\sigma_0, \sigma_1) = (0, 0)$ is reached in finite time, and thus the condition $s_0(t) = \dot{s}_0(t) = 0$ is attained and maintained after a finite time [82]. Hence, the differentiator attains a 2-sliding surface. On the sliding surface $s_0(t) \equiv 0$ we have $z_1(t) \equiv \dot{f}(t)$, which shows that the signal's derivative is available in finite-time.

The backward-Euler (or, implicit-Euler) discretization of (107) is

$$ s_{0,k+1} = s_{0,k} + h z_{1,k+1} - h \gamma_0 L^{\frac{1}{2}} [s_{0,k+1}]^{\frac{1}{2}} - \Delta f_k \qquad (109a) $$

$$ z_{1,k+1} \in z_{1,k} - h\gamma_1 L \operatorname{sgn}(s_{0,k+1}), \qquad (109b) $$

where $\Delta f_k := f_{k+1} - f_k$. Our purpose now is to derive an explicit representation of (109). Substituting (109b) into (109a) and solving for $s_{0,k+1}$ yields

$$ s_{0,k+1} = (I_d + h M_h)^{-1} (s_{0,k} + h z_{1,k} - \Delta f_k), \qquad (110) $$

where

$$ M_h(w) = \gamma_0 L^{\frac{1}{2}} [w]^{\frac{1}{2}} + h \gamma_1 L \operatorname{sgn}(w) \qquad (111) $$

is a maximal monotone operator (being the sum of two maximal monotone operators, both with domain $\mathbb{R}$). The iteration (110) is compactly written as the resolvent-based algorithm

$$ s_{0,k+1} = \mathcal{J}_h M_h (s_{0,k} + h z_{1,k} - \Delta f_k) \qquad (112) $$

(see (57) and Proposition 8 in "Sidebar VII"), where the resolvent is defined in (23) in "Sidebar V" (compare with the classical algorithm (96) and the general iteration (57)). It is interesting to note that choosing an exact discretization as in [83] and neglecting unmeasurable terms also yields a form like (112) [83, Equation (24)]. Together with (10), (73), and (96), this shows the ubiquity of resolvent-based algorithms in the implicit discretization of sliding-mode systems. The algorithms provided in [84, 83] also correspond to the resolvent in (112). As seen later in this article, there is no unique way to compute the resolvents numerically.

With an explicit expression for $s_{0,k+1}$ available in (112), we can readily solve for the state $z_{1,k+1}$ for all $h > 0$ in (109a),

$$ z_{1,k+1} = \frac{1}{h} \left( s_{0,k+1} - s_{0,k} + h \gamma_0 L^{\frac{1}{2}} [s_{0,k+1}]^{\frac{1}{2}} + \Delta f_k \right). \qquad (113) $$

Therefore, the pair of equations (112), (113) forms a non-linear difference equation with state vector $(s_{0,k} \ z_{1,k})^\top$ and allows us to advance the discretized super-twisting differentiator from step $k$ to step $k+1$. Notice in passing that, having computed the future state $z_{1,k+1}$ in (113), the equation $z_{1,k+1} = z_{1,k} - h\gamma_1 L \lambda_{k+1}$ can be used to compute a selection $\lambda_{k+1} \in \operatorname{sgn}(s_{0,k+1})$ such that (109b) holds.

**Remark 8.** In real-time applications it is possible to substitute $\Delta f_k$ by $\Delta f_{k-1}$ in (109a), and hence in (113), to render the differentiator non-anticipative.

Semi-implicit discretization is also possible. Let us replace (109) by

$$ s_{0,k+1} \in s_{0,k} + h z_{1,k+1} - h \gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} \operatorname{sgn}(s_{0,k+1}) - \Delta f_k \qquad (114a) $$

$$ z_{1,k+1} \in z_{1,k} - h\gamma_1 L \operatorname{sgn}(s_{0,k+1}), \qquad (114b) $$

where the implicit terms have been restricted to the set-valued functions. Solving again for $s_{0,k+1}$ leads now to

$$ 
\begin{aligned}
& s_{0,k+1} \\
&= (I_d + (h\gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} + h^2 \gamma_1 L) \operatorname{sgn})^{-1} (s_{0,k} + h z_{1,k} - \Delta f_k).
\end{aligned}
\qquad (115) $$

Compared to (111), the maximal monotone operator is modified to

$$ M_{k,h}(w) = (\gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} + h \gamma_1 L) \operatorname{sgn}(w). \qquad (116) $$
practice is to replace $\dot{f}$ by its *dirty derivative*, $\hat{g}$, which is obtained by passing $f$ through a linear filter of the form

$$ 
\begin{aligned}
\dot{z}(t) &= -\gamma(z(t) + \gamma f(t)) \\
\hat{g}(t) &= z(t) + \gamma f(t)
\end{aligned}
\qquad (105) $$

with $\gamma \gg 1$. It can be easily shown that, as $\gamma$ approaches $+\infty$, the transfer function of the filter converges to the transfer function of a differentiator. This high-gain strategy benefits from simplicity, but it has two important drawbacks: noise amplification and phase lag. If $\dot{f}$ is Lipschitz continuous with Lipschitz constant $L > 0$, then a super-twisting differentiation algorithm can be used. The advantages of such scheme are: 1) exact differentiation with finite gains in the absence of noise, 2) limited noise amplification, 3) no phase lag.

We will make use of the operator $[w]^\alpha := |w|^\alpha \operatorname{sgn}(w)$, $\alpha \ge 0$, $w \in \mathbb{R}$. Notice that $[w]^\alpha = \partial\zeta_\alpha(w)$ with

$$ \zeta_\alpha(w) = \frac{1}{\alpha+1} |w|^{\alpha+1}. \qquad (106) $$

The convexity of $\zeta_\alpha$ immediately guarantees the maximal monotonicity of $[\cdot]^\alpha$, see Theorem 3 in "Sidebar V". Notice also that $[\cdot]^\alpha$ is single-valued for $\alpha > 0$, but multivalued at the origin when $\alpha = 0$. The super-twisting differentiation algorithm [6] reads as:

$$ 
\begin{aligned}
\dot{z}_0(t) &= z_1(t) - \gamma_0 L^{\frac{1}{2}} [s_0(t)]^{\frac{1}{2}}, \\
\dot{z}_1(t) &\in -\gamma_1 L \operatorname{sgn}(s_0(t))
\end{aligned}
\qquad (107) $$

where $s_0(t) := z_0(t) - f(t)$ and $\gamma_0, \gamma_1$ are positive gains. The differentiator properties are best appreciated under the coordinate transformation $\sigma_0(t) := \frac{1}{L}(z_0 - f(t))$ and $\sigma_1(t) := \frac{1}{L}(z_1(t) - \dot{f}(t))$. In these coordinates, system (107) is rewritten in recursive form as [81]

$$ 
\begin{aligned}
\dot{\sigma}_0(t) &= \sigma_1(t) - \gamma_0 [\sigma_0(t)]^{\frac{1}{2}} \\
\dot{\sigma}_1(t) &\in -\gamma_1 \operatorname{sgn}(\sigma_0(t)) - \ddot{f}(t)
\end{aligned}
\qquad (108) $$

If $\gamma_0$ and $\gamma_1$ are appropriately chosen, then the origin $(\sigma_0, \sigma_1) = (0,0)$ is reached in finite time, and thus the condition $s_0(t) = \dot{s}_0(t) = 0$ is attained and maintained after a finite time [82]. Hence, the differentiator attains a 2-sliding surface. On the sliding surface $s_0(t) \equiv 0$ we have $z_1(t) \equiv \dot{f}(t)$, which shows that the signal's derivative is available in finite-time.

The backward-Euler (or, implicit-Euler) discretization of (107) is

$$ s_{0,k+1} = s_{0,k} + h z_{1,k+1} - h \gamma_0 L^{\frac{1}{2}} [s_{0,k+1}]^{\frac{1}{2}} - \Delta f_k \qquad (109a) $$

$$ z_{1,k+1} \in z_{1,k} - h\gamma_1 L \operatorname{sgn}(s_{0,k+1}), \qquad (109b) $$

where $\Delta f_k := f_{k+1} - f_k$. Our purpose now is to derive an explicit representation of (109). Substituting (109b) into (109a) and solving for $s_{0,k+1}$ yields

$$ s_{0,k+1} = (I_d + h \mathbf{M}_h)^{-1} (s_{0,k} + h z_{1,k} - \Delta f_k), \qquad (110) $$

where

$$ \mathbf{M}_h(w) = \gamma_0 L^{\frac{1}{2}} [w]^{\frac{1}{2}} + h \gamma_1 L \operatorname{sgn}(w) \qquad (111) $$

is a maximal monotone operator (being the sum of two maximal monotone operators, both with domain $\mathbb{R}$). The iteration (110) is compactly written as the resolvent-based algorithm

$$ s_{0,k+1} = \mathcal{J}_h \mathbf{M}_h (s_{0,k} + h z_{1,k} - \Delta f_k) \qquad (112) $$

(see (57) and Proposition 8 in "Sidebar VII"), where the resolvent is defined in (23) in "Sidebar V" (compare with the classical algorithm (96) and the general iteration (57)). It is interesting to note that choosing an exact discretization as in [83] and neglecting unmeasurable terms also yields a form like (112) [83, Equation (24)]. Together with (10), (73), and (96), this shows the ubiquity of resolvent-based algorithms in the implicit discretization of sliding-mode systems. The algorithms provided in [84, 83] also correspond to the resolvent in (112). As seen later in this article, there is no unique way to compute the resolvents numerically.

With an explicit expression for $s_{0,k+1}$ available in (112), we can readily solve for the state $z_{1,k+1}$ for all $h > 0$ in (109a),

$$ z_{1,k+1} = \frac{1}{h} \left( s_{0,k+1} - s_{0,k} + h \gamma_0 L^{\frac{1}{2}} [s_{0,k+1}]^{\frac{1}{2}} + \Delta f_k \right). \qquad (113) $$

Therefore, the pair of equations (112), (113) forms a non-linear difference equation with state vector $(s_{0,k} \ z_{1,k})^\top$ and allows us to advance the discretized super-twisting differentiator from step $k$ to step $k+1$. Notice in passing that, having computed the future state $z_{1,k+1}$ in (113), the equation $z_{1,k+1} = z_{1,k} - h\gamma_1 L \lambda_{k+1}$ can be used to compute a selection $\lambda_{k+1} \in \operatorname{sgn}(s_{0,k+1})$ such that (109b) holds.

**Remark 8.** In real-time applications it is possible to substitute $\Delta f_k$ by $\Delta f_{k-1}$ in (109a), and hence in (113), to render the differentiator non-anticipative.

Semi-implicit discretization is also possible. Let us replace (109) by

$$ s_{0,k+1} \in s_{0,k} + h z_{1,k+1} - h \gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} \operatorname{sgn}(s_{0,k+1}) - \Delta f_k \qquad (114a) $$

$$ z_{1,k+1} \in z_{1,k} - h\gamma_1 L \operatorname{sgn}(s_{0,k+1}), \qquad (114b) $$

where the implicit terms have been restricted to the set-valued functions. Solving again for $s_{0,k+1}$ leads now to

$$ 
\begin{aligned}
& s_{0,k+1} \\
&= (I_d + (h\gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} + h^2\gamma_1 L) \operatorname{sgn})^{-1} (s_{0,k} + h z_{1,k} - \Delta f_k).
\end{aligned}
\qquad (115) $$

Compared to (111), the maximal monotone operator is modified to

$$ \mathbf{M}_{k,h}(w) = (\gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} + h \gamma_1 L) \operatorname{sgn}(w). \qquad (116) $$
Computing $z_{1,k+1}$ is somewhat more complicated than before. First let us define the selection $\lambda_{k+1} \in \mathbf{sgn}(s_{0,k+1})$ and recall that it holds if, and only if, $s_{0,k+1} \in \mathbf{N}_{[-1,1]}(\lambda_{k+1})$ (see Definition 3 in "Sidebar III", Fact 7 and (37) in "Sidebar VI"). By inverting the operators in (114), the system in the new coordinates can be rewritten as

$$s_{0,k} + h z_{1,k+1} - h \gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} \lambda_{k+1} - \Delta f_k \in \mathbf{N}_{[-1,1]}(\lambda_{k+1}) \quad (117a)$$

$$z_{1,k+1} = z_{1,k} - h\gamma_1 L \lambda_{k+1} \quad (117b)$$

and then solved for $\lambda_{k+1}$. To do so, first substitute (117b) in (117a) to obtain

$$\lambda_{k+1} \in \frac{1}{h^2 \gamma_1 L} \left( s_{0,k} + h z_{1,k} - h \gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} \lambda_{k+1} - \Delta f_k - \mathbf{N}_{[-1,1]}(\lambda_{k+1}) \right). \quad (118)$$

The inclusion is equivalent to

$$\lambda_{k+1} = \mathcal{J}_{M_{k,h}} \left( \frac{s_{0,k} + h z_{1,k} - \Delta f_k}{h^2 \gamma_1 L} \right), \quad (119)$$

$$\mathbf{M}_{k,h}(w) = \frac{\gamma_0}{h \gamma_1 L^{\frac{1}{2}}} |s_{0,k}|^{\frac{1}{2}} w + \mathbf{N}_{[-1,1]}(w). \quad (120)$$

Finally, $z_{1,k+1}$ is recovered from (117b) as

$$z_{1,k+1} = z_{1,k} - h\gamma_1 L \mathcal{J}_{M_{k,h}} \left( \frac{s_{0,k} + h z_{1,k} - \Delta f_k}{h^2 \gamma_1 L} \right). \quad (121)$$

The resolvent $\mathcal{J}_{M_{k,h}}(\cdot)$ can be calculated as indicated in "Sidebar VII", after (55). It happens that resolvents of this kind are basic tools which are to be found in all the higher-order algorithms in this section.

### 4.1.2 Controller

When the super-twisting algorithm is used for control (as opposed to differentiation), a perturbed integrator $\dot{x}_1(t) = u(t) + \varphi(t)$ is considered. This dynamics may stem from some transformation of the plant's dynamics, and $x_1$ may represent a sliding variable (similar to $s$ in the first order schemes), hence it may be a combination of the plant's state variables [85]. The proposed controller is given by

$$\begin{aligned} u(t) &= \nu(t) - \gamma_0 L^{\frac{1}{2}} [x_1(t)]^{\frac{1}{2}}. \\ \dot{\nu}(t) &\in -\gamma_1 L \operatorname{sgn}(x_1(t)) \end{aligned} \quad (122)$$

It is apparent that $u$ in (122) is a kind of nonlinear nonsmooth set-valued proportional-integral input. The main interest in this controller stems from the fact that, contrarily to the first-order case, the applied control, $u$, is a continuous function of the state. By making the time-varying change of variable $x_2(t) = \nu(t) + \varphi(t)$, one obtains the closed-loop system

$$\begin{aligned} \dot{x}_1(t) &= x_2(t) - \gamma_0 L^{\frac{1}{2}} [x_1(t)]^{\frac{1}{2}} \\ \dot{x}_2(t) &\in -\gamma_1 L \operatorname{sgn}(x_1(t)) + \delta(t) \end{aligned} \quad (123)$$

with $\delta(t) = \dot{\varphi}(t)$. The closed-loop dynamics (123) is the same as the differentiator's dynamics (108). Finite-time stability of the origin $x_1 = x_2 = 0$ is guaranteed under similar assumptions on the gains [82], hence $x_1 = \dot{x}_1 = 0$ defines the 2-sliding surface. Notice that, similarly to (21) and (39), the controller compensates exactly for the perturbation on the sliding surface $\{0\} \subset \mathbb{R}^2$, along which $u(t) = \nu(t) = -\varphi(t)$.

The objective now is to compute a discrete-time version of the controller, similar to the explicit expressions (12) and (13). A discrete-time model for the plant is

$$x_{1,k+1} = x_{1,k} + h u_k + h \bar{\varphi}_k \quad (124a)$$

$$\varphi_{k+1} = \varphi_k + h \bar{\delta}_k \quad (124b)$$

Here, $\bar{\varphi}_k$ and $\bar{\delta}_k$ are discretized versions of $\varphi(t)$ and $\delta(t)$ (see [34] for details). The implicit discretization can now be achieved by following the same steps as in (109). Let us introduce, for instance, a nominal virtual variable $\tilde{x}_{1,k} := x_{1,k} - h\bar{\varphi}_{k-1}$ as done in (62) and (93), see [34, Equations (8), (9)]. According to (124a), the nominal variable evolves as

$$\tilde{x}_{1,k+1} = x_{1,k} + h u_k, \quad (125)$$

while the nominal backward-Euler discretization of (122) is

$$\begin{aligned} u_k &\in -\gamma_0 L^{\frac{1}{2}} |\tilde{x}_{1,k+1}|^{\frac{1}{2}} \mathbf{sgn}(\tilde{x}_{1,k+1}) + \nu_{k+1} \\ \nu_{k+1} &= \nu_k - h \gamma_1 L \mathbf{sgn}(\tilde{x}_{1,k+1}). \end{aligned} \quad (126)$$

Substitution of (126) in (125) gives

$$\begin{aligned} \tilde{x}_{1,k+1} &= (I_d + h \mathbf{M}_h)^{-1} (x_{1,k} + h \nu_k) \\ &= \mathcal{J}_h \mathbf{M}_h (x_{1,k} + h \nu_k), \end{aligned} \quad (127)$$

with $\mathbf{M}_h$ as in (111). The explicit expression for the control in (126) is [34, 86]

$$\begin{aligned} u(x_{1,k}, \nu_k) &= -\gamma_0 L^{\frac{1}{2}} \beta_k \operatorname{sgn}(x_{1,k} + h\nu_k) \\ &\quad + \nu_k - h\gamma_1 L \operatorname{Proj}([-1, 1]; \frac{x_{1,k} + h\nu_k}{h^2 \gamma_0 L}) \end{aligned} \quad (128)$$

with

$$\beta_k = \frac{-h \gamma_0 L^{\frac{1}{2}}}{2} + \sqrt{\frac{h^2 \gamma_0^2 L}{4} + \rho_k} \quad (129)$$

$$\rho_k = \max\{0, |x_{1,k} + h\nu_k| - h^2 \gamma_1 L\}, \quad (130)$$

while

$$\nu_{k+1} = \nu_k - h\gamma_1 \operatorname{Proj}([-1, 1]; \frac{x_{0,k} + h\nu_k}{h\mu_k}). \quad (131)$$

Alternatively, we may rewrite the control as

$$u(x_{1,k}, \nu_k) = \frac{1}{h}(I_d + h \mathbf{M}_h)^{-1}(x_{1,k} + h \nu_k) - \frac{x_{1,k}}{h}. \quad (132)$$

It is noteworthy that this boils down once again to the calculation of a resolvent and of a projection (the saturation
Computing $z_{1,k+1}$ is somewhat more complicated than before. First let us define the selection $\lambda_{k+1} \in \mathbf{sgn}(s_{0,k+1})$ and recall that it holds if, and only if, $s_{0,k+1} \in \mathbf{N}_{[-1,1]}(\lambda_{k+1})$ (see Definition 3 in “Sidebar III”, Fact 7 and (37) in “Sidebar VI”). By inverting the operators in (114), the system in the new coordinates can be rewritten as

$$s_{0,k} + h z_{1,k+1} - h \gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} \lambda_{k+1} - \Delta f_k \in \mathbf{N}_{[-1,1]}(\lambda_{k+1}) \quad (117a)$$

$$z_{1,k+1} = z_{1,k} - h\gamma_1 L \lambda_{k+1} \quad \text{and (13). A discrete-time model for the plant is} \quad (117b)$$

and then solved for $\lambda_{k+1}$. To do so, first substitute (117b) in (117a) to obtain

$$\lambda_{k+1} \in \frac{1}{h^2 \gamma_1 L} \left( s_{0,k} + h z_{1,k} - h \gamma_0 L^{\frac{1}{2}} |s_{0,k}|^{\frac{1}{2}} \lambda_{k+1} - \Delta f_k - \mathbf{N}_{[-1,1]}(\lambda_{k+1}) \right). \quad (118)$$

The inclusion is equivalent to

$$\lambda_{k+1} = \mathcal{J}_{M_{k,h}} \left( \frac{s_{0,k} + h z_{1,k} - \Delta f_k}{h^2 \gamma_1 L} \right), \quad (119)$$

$$\mathbf{M}_{k,h}(w) = \frac{\gamma_0}{h \gamma_1 L^{\frac{1}{2}}} |s_{0,k}|^{\frac{1}{2}} w + \mathbf{N}_{[-1,1]}(w). \quad (120)$$

Finally, $z_{1,k+1}$ is recovered from (117b) as

$$z_{1,k+1} = z_{1,k} - h\gamma_1 L \mathcal{J}_{M_{k,h}} \left( \frac{s_{0,k} + h z_{1,k} - \Delta f_k}{h^2 \gamma_1 L} \right). \quad (121)$$

The resolvent $\mathcal{J}_{M_{k,h}}(\cdot)$ can be calculated as indicated in “Sidebar VII”, after (55). It happens that resolvents of this kind are basic tools which are to be found in all the higher-order algorithms in this section.

### 4.1.2 Controller

When the super-twisting algorithm is used for control (as opposed to differentiation), a perturbed integrator $\dot{x}_1(t) = u(t) + \varphi(t)$ is considered. This dynamics may stem from some transformation of the plant's dynamics, and $x_1$ may represent a sliding variable (similar to $s$ in the first order schemes), hence it may be a combination of the plant's state variables [85]. The proposed controller is given by

$$\begin{aligned} u(t) &= \nu(t) - \gamma_0 L^{\frac{1}{2}} [x_1(t)]^{\frac{1}{2}}. \\ \dot{\nu}(t) &\in -\gamma_1 L \operatorname{sgn}(x_1(t)) \end{aligned} \quad (122)$$

It is apparent that $u$ in (122) is a kind of nonlinear nonsmooth set-valued proportional-integral input. The main interest in this controller stems from the fact that, contrarily to the first-order case, the applied control, $u$, is a continuous function of the state. By making the time-varying change of variable $x_2(t) = \nu(t) + \varphi(t)$, one obtains the closed-loop system

$$\begin{aligned} \dot{x}_1(t) &= x_2(t) - \gamma_0 L^{\frac{1}{2}} [x_1(t)]^{\frac{1}{2}} \\ \dot{x}_2(t) &\in -\gamma_1 L \operatorname{sgn}(x_1(t)) + \delta(t) \end{aligned} \quad (123)$$

with $\delta(t) = \dot{\varphi}(t)$. The closed-loop dynamics (123) is the same as the differentiator's dynamics (108). Finite-time stability of the origin $x_1 = x_2 = 0$ is guaranteed under similar assumptions on the gains [82], hence $x_1 = \dot{x}_1 = 0$ defines the 2-sliding surface. Notice that, similarly to (21) and (39), the controller compensates exactly for the perturbation on the sliding surface $\{0\} \subset \mathbb{R}^2$, along which $u(t) = \nu(t) = -\varphi(t)$.

The objective now is to compute a discrete-time version of the controller, similar to the explicit expressions (12) and (13). A discrete-time model for the plant is

$$x_{1,k+1} = x_{1,k} + h u_k + h \bar{\varphi}_k \quad (124a)$$

$$\varphi_{k+1} = \varphi_k + h \bar{\delta}_k \quad (124b)$$

Here, $\bar{\varphi}_k$ and $\bar{\delta}_k$ are discretized versions of $\varphi(t)$ and $\delta(t)$ (see [34] for details). The implicit discretization can now be achieved by following the same steps as in (109). Let us introduce, for instance, a nominal virtual variable $\tilde{x}_{1,k} := x_{1,k} - h\bar{\varphi}_{k-1}$ as done in (62) and (93), see [34, Equations (8), (9)]. According to (124a), the nominal variable evolves as

$$\tilde{x}_{1,k+1} = x_{1,k} + h u_k, \quad (125)$$

while the nominal backward-Euler discretization of (122) is

$$\begin{aligned} u_k &\in -\gamma_0 L^{\frac{1}{2}} |\tilde{x}_{1,k+1}|^{\frac{1}{2}} \mathbf{sgn}(\tilde{x}_{1,k+1}) + \nu_{k+1} \\ \nu_{k+1} &= \nu_k - h \gamma_1 L \mathbf{sgn}(\tilde{x}_{1,k+1}). \end{aligned} \quad (126)$$

Substitution of (126) in (125) gives

$$\begin{aligned} \tilde{x}_{1,k+1} &= (I_d + h \mathbf{M}_h)^{-1} (x_{1,k} + h \nu_k) \\ &= \mathcal{J}_h \mathbf{M}_h (x_{1,k} + h \nu_k), \end{aligned} \quad (127)$$

with $\mathbf{M}_h$ as in (111). The explicit expression for the control in (126) is [34, 86]

$$\begin{aligned} u(x_{1,k}, \nu_k) &= -\gamma_0 L^{\frac{1}{2}} \beta_k \operatorname{sgn}(x_{1,k} + h\nu_k) \\ &\quad + \nu_k - h\gamma_1 L \operatorname{Proj}([-1, 1]; \frac{x_{1,k} + h\nu_k}{h^2 \gamma_0 L}) \end{aligned} \quad (128)$$

with

$$\beta_k = \frac{-h \gamma_0 L^{\frac{1}{2}}}{2} + \sqrt{\frac{h^2 \gamma_0^2 L}{4} + \rho_k} \quad (129)$$

$$\rho_k = \max\{0, |x_{1,k} + h\nu_k| - h^2 \gamma_1 L\}, \quad (130)$$

while

$$\nu_{k+1} = \nu_k - h\gamma_1 \operatorname{Proj}([-1, 1]; \frac{x_{0,k} + h\nu_k}{h\mu_k}). \quad (131)$$

Alternatively, we may rewrite the control as

$$u(x_{1,k}, \nu_k) = \frac{1}{h}(I_d + h \mathbf{M}_h)^{-1}(x_{1,k} + h \nu_k) - \frac{x_{1,k}}{h}. \quad (132)$$

It is noteworthy that this boils down once again to the calculation of a resolvent and of a projection (the saturation
function). That is, equation (127) defines a resolvent-based algorithm similar to (112) and (114).

As above, a semi-implicit method can be applied [35].
It yields

$$
\tilde{x}_{1,k+1} = (\mathbf{I}_d + h \mathbf{M}_{k,h})^{-1} (x_{1,k} + h\nu_k) \quad (133)
$$

but now with

$$
\mathbf{M}_{k,h}(w) = (\gamma_0 L^{\frac{1}{2}} |x_{1,k}|^{\frac{1}{2}} + h\gamma_1 L) \operatorname{sgn}(w), \quad (134)
$$

similar to (116). Another choice is

$$
\mathbf{M}_{k,h}(w) = (\gamma_0 L^{\frac{1}{2}} |\tilde{x}_{1,k}|^{\frac{1}{2}} + h\gamma_1 L) \operatorname{sgn}(w). \quad (135)
$$

Which alternative should be chosen is an open issue.

Note that, during the sliding motion, $\tilde{x}_{1,k+1} = \tilde{x}_{1,k} = 0$, and the inclusion

$$
\nu_k + \bar{\varphi}_k \in \frac{1}{h} (\mathbf{I}_d + h \mathbf{M}_h)(0) = h\gamma_0 L[-1, 1] \quad (136)
$$

follows from (127). Using (128), (130) we obtain

$$
u(x_{1,k}, \nu_k) = -\frac{x_{1,k}}{h} = -\frac{\tilde{x}_{1,k} + h\bar{\varphi}_{k-1}}{h} = -\bar{\varphi}_{k-1}. \quad (137)
$$

It is inferred that, similarly to first order controllers, see
(98), the implicit super-twisting scheme compensates for
the perturbation with one-step delay.

4.1.3 Multivariable super-twisting

The above extends to the multivariable super-twisting
[87], with $x_1, x_2 \in \mathbb{R}^n$. One replaces (106) by

$$
\zeta_{\alpha}(w) = \frac{1}{\alpha + 1} \|w\|_{2}^{\alpha + 1}, \quad (138)
$$

so that instead of $|x_1|^{\frac{1}{2}} \operatorname{sgn}(x_1)$ one has

$$
\|x_1\|_2^{\frac{1}{2}} \partial \|x_1\|_2 = \|x_1\|_2^{\frac{1}{2}} \mathbf{N}_{\mathbf{B}_2}^{-1}(x_1) \quad (139)
$$

(see (41) in “Sidebar VI”), and **sgn**(*x*₁) by ∂‖*x*₁‖₂ (or the
subdifferential of any other norm), in (122). This yields
the resolvent J<sub>h</sub> **M**<sub>h</sub>, to be computed at each step with

$$
\mathbf{M}_h(w) = \alpha_1 \|w\|_2^{\frac{1}{2}} \partial \|w\|_2 + h\alpha_2 \partial \|w\|_2, \quad (140)
$$

$\alpha_1 > 0, \alpha_2 > 0$. Note that $\mathbf{M}_h = \partial(\alpha_1 \zeta_{1/2} + h\alpha_2 \zeta_0)$.
It follows from the convexity of $\zeta_\alpha$ and Theorem 3 in
“Sidebar V” that $\mathbf{M}_h$ is maximal monotone.

4.1.4 Generalized Observer

The equivalence between uniform observability and the
ability to recover the complete state by differentiating the
system output is well known. In view of (107), it is not
surprising that the super-twisting algorithm can be used
to observe planar systems. Recall that a planar system is

uniformly observable if and only if there exists a coordi-
nate transformation such that the system can be written
as [88]

$$
\begin{align}
\dot{x}_1 &= f_1(x_1, u) + x_2 \\
\dot{x}_2 &= f_2(x_1, x_2, u) , \tag{141} \\
y &= x_1
\end{align}
$$

from where it is clear that the remaining state $x_2$ can be
recovered by differentiating the output $y$. Let us consider
the class of set-valued super-twisting observers presented
in [89, Equations (12.3)–(12.5)]

$$
\dot{x}_1 = f_1(\hat{x}_1, u) + \hat{x}_2 - l_1 \gamma \phi_1(e_1) \quad (142a)
$$

$$
\dot{x}_2 \in f_2(x_1, \hat{x}_2, u) - l_2 \gamma^2 \phi_2(e_1) \quad (142b)
$$

$$
\phi_1(e_1) = \mu_1 [e_1]^{\frac{1}{2}} + \mu_2 [e_1]^q \quad (142c)
$$

$$
\phi_2(e_1) = \frac{\mu_1^2}{2} \operatorname{sgn}(e_1) + \mu_1 \mu_2 [e_1]^{q-\frac{1}{2}} + \mu_2^2 [e_1]^{2q-1}, \quad (142\text{d})
$$

where $\mu_1 \ge 0$, $\mu_2 \ge 0$, $\mu_1\mu_2 \ne 0$, $q \ge \frac{1}{2}$, $l_1 > 0$, $l_2 > 0$,
$\gamma > 0$, and $e_1 = \hat{x}_1 - x_1$ is the state estimation error. An
implicit discretization of (142) reads

$$
\hat{x}_{1,k+1} = \hat{x}_{1,k} + h f_1(\hat{x}_{1,k+1}, u_k) + h \hat{x}_{2,k+1} - h l_1 \gamma \phi_1(e_{1,k+1}) \quad (143a)
$$

$$
\hat{x}_{2,k+1} \in \hat{x}_{2,k} + h f_2(\hat{x}_{1,k+1}, \hat{x}_{2,k+1}, u_k) - h l_2 \gamma^2 \phi_2(e_{1,k+1}) \quad (143b)
$$

with $e_{1,k+1} = \hat{x}_{1,k+1} - y_k$. Noting that $\phi_2(e_{1,k+1})$ is a set-valued function of $e_{1,k+1}$, it is assumed that (143b) can be solved to get $\hat{x}_{2,k+1} \in \Phi(e_{1,k+1}, y_k, \hat{x}_{2,k}, u_k, h\phi_2(e_{1,k+1}))$ for some function $\Phi$. Inserting the latter in (143a) yields

$$
\begin{equation}
\begin{split}
e_{1,k+1} \in y_{k-1} - y_k + e_{1,k} &+ h f_1(e_{1,k+1}, y_k, u_k) \\
&+ h \Phi(e_{1,k+1}, y_k, \hat{x}_{2,k}, u_k, h \phi_2(e_{1,k+1})) \\
&- h l_1 \gamma \phi_1(e_{1,k+1}).
\end{split}
\tag{144}
\end{equation}
$$

It is inferred that

$$
\begin{equation}
\begin{split}
e_{1,k+1} \in (\mathbf{I}_d - h f_1(\cdot, y_k, u_k) &+ h \Phi(\cdot, y_k, \hat{x}_{2,k}, u_k, h \phi_2(\cdot)) \\
&+ h l_1 \gamma \phi_1(\cdot))^{-1} (y_{k-1} - y_k + e_{1,k}).
\end{split}
\tag{145}
\end{equation}
$$

Then, provided that the operators $x \mapsto -f_1(x, y_k, u_k) + l_1 \gamma \phi_1(x)$ and $x \mapsto \Phi(x, y_k, \hat{x}_{2,k}, u_k, h\phi_2(x))$ are maximal monotone (notice that $\phi_1$ and $\phi_2$ are strongly monotone for $q > \frac{1}{2}$), the inclusion in (145) is an equality and the right-hand side is a resolvent. Depending on the nonlinearities, the closed form of this resolvent may be impossible to obtain, in which case the numerical solutions described below are necessary. The observed variable $\hat{x}_{2,k+1}$ is obtained from (143a) and (145). Of course, semi-implicit methods can also be applied so that the computations are easier and the assumptions may be relaxed.
function). That is, equation (127) defines a resolvent-based algorithm similar to (112) and (114).

As above, a semi-implicit method can be applied [35].
It yields

$$
\tilde{x}_{1,k+1} = (\mathbf{I}_d + h \mathbf{M}_{k,h})^{-1} (x_{1,k} + h\nu_k) \quad (133)
$$

but now with

$$
\mathbf{M}_{k,h}(w) = (\gamma_0 L^{\frac{1}{2}} |x_{1,k}|^{\frac{1}{2}} + h\gamma_1 L) \operatorname{sgn}(w), \quad (134)
$$

similar to (116). Another choice is

$$
\mathbf{M}_{k,h}(w) = (\gamma_0 L^{\frac{1}{2}} |\tilde{x}_{1,k}|^{\frac{1}{2}} + h\gamma_1 L) \operatorname{sgn}(w). \quad (135)
$$

Which alternative should be chosen is an open issue.

Note that, during the sliding motion, $\tilde{x}_{1,k+1} = \tilde{x}_{1,k} = 0$, and the inclusion

$$
\nu_k + \bar{\varphi}_k \in \frac{1}{h}(\mathbf{I}_d + h \mathbf{M}_h)(0) = h\gamma_0 L[-1, 1] \quad (136)
$$

follows from (127). Using (128), (130) we obtain

$$
u(x_{1,k}, \nu_k) = -\frac{x_{1,k}}{h} = -\frac{\tilde{x}_{1,k} + h\bar{\varphi}_{k-1}}{h} = -\bar{\varphi}_{k-1}. \quad (137)
$$

It is inferred that, similarly to first order controllers, see
(98), the implicit super-twisting scheme compensates for
the perturbation with one-step delay.

4.1.3 Multivariable super-twisting

The above extends to the multivariable super-twisting
[87], with $x_1, x_2 \in \mathbb{R}^n$. One replaces (106) by

$$
\zeta_{\alpha}(w) = \frac{1}{\alpha + 1} \|w\|_{2}^{\alpha + 1}, \quad (138)
$$

so that instead of $|x_1|^{\frac{1}{2}} \operatorname{sgn}(x_1)$ one has

$$
\|x_1\|_2^{\frac{1}{2}} \partial \|x_1\|_2 = \|x_1\|_2^{\frac{1}{2}} \mathbf{N}_{\mathbf{B}_2}^{-1}(x_1) \quad (139)
$$

(see (41) in “Sidebar VI”), and **sgn**(*x*₁) by ∂‖*x*₁‖₂ (or the
subdifferential of any other norm), in (122). This yields
the resolvent J<sub>h</sub> **M**<sub>h</sub>, to be computed at each step with

$$
\mathbf{M}_h(w) = \alpha_1 \|w\|_2^{\frac{1}{2}} \partial \|w\|_2 + h\alpha_2 \partial \|w\|_2, \quad (140)
$$

$\alpha_1 > 0, \alpha_2 > 0$. Note that $\mathbf{M}_h = \partial(\alpha_1 \zeta_{1/2} + h\alpha_2 \zeta_0)$.
It follows from the convexity of $\zeta_\alpha$ and Theorem 3 in
“Sidebar V” that $\mathbf{M}_h$ is maximal monotone.

4.1.4 Generalized Observer

The equivalence between uniform observability and the
ability to recover the complete state by differentiating the
system output is well known. In view of (107), it is not
surprising that the super-twisting algorithm can be used
to observe planar systems. Recall that a planar system is

uniformly observable if and only if there exists a coordi-
nate transformation such that the system can be written
as [88]

$$
\begin{align}
\dot{x}_1 &= f_1(x_1, u) + x_2 \\
\dot{x}_2 &= f_2(x_1, x_2, u) , \tag{141} \\
y &= x_1
\end{align}
$$

from where it is clear that the remaining state $x_2$ can be
recovered by differentiating the output $y$. Let us consider
the class of set-valued super-twisting observers presented
in [89, Equations (12.3)–(12.5)]

$$
\dot{x}_1 = f_1(\hat{x}_1, u) + \hat{x}_2 - l_1 \gamma \phi_1(e_1) \quad (142a)
$$

$$
\dot{x}_2 \in f_2(x_1, \hat{x}_2, u) - l_2 \gamma^2 \phi_2(e_1) \quad (142b)
$$

$$
\phi_1(e_1) = \mu_1 [e_1]^{\frac{1}{2}} + \mu_2 [e_1]^q \quad (142c)
$$

$$
\phi_2(e_1) = \frac{\mu_1^2}{2} \operatorname{sgn}(e_1) + \mu_1 \mu_2 [e_1]^{q-\frac{1}{2}} + \mu_2^2 [e_1]^{2q-1}, \quad (142\text{d})
$$

where $\mu_1 \ge 0$, $\mu_2 \ge 0$, $\mu_1\mu_2 \ne 0$, $q \ge \frac{1}{2}$, $l_1 > 0$, $l_2 > 0$, $\gamma > 0$, and $e_1 = \hat{x}_1 - x_1$ is the state estimation error. An implicit discretization of (142) reads

$$
\hat{x}_{1,k+1} = \hat{x}_{1,k} + h f_1(\hat{x}_{1,k+1}, u_k) + h \hat{x}_{2,k+1} - h l_1 \gamma \phi_1(e_{1,k+1}) \quad (143a)
$$

$$
\hat{x}_{2,k+1} \in \hat{x}_{2,k} + h f_2(\hat{x}_{1,k+1}, \hat{x}_{2,k+1}, u_k) - h l_2 \gamma^2 \phi_2(e_{1,k+1}) \quad (143b)
$$

with $e_{1,k+1} = \hat{x}_{1,k+1} - y_k$. Noting that $\phi_2(e_{1,k+1})$ is a set-valued function of $e_{1,k+1}$, it is assumed that (143b) can be solved to get $\hat{x}_{2,k+1} \in \Phi(e_{1,k+1}, y_k, \hat{x}_{2,k}, u_k, h\phi_2(e_{1,k+1}))$ for some function $\Phi$. Inserting the latter in (143a) yields

$$
\begin{equation}
\begin{split}
e_{1,k+1} \in y_{k-1} - y_k + e_{1,k} &+ h f_1(e_{1,k+1}, y_k, u_k) \\
&+ h \Phi(e_{1,k+1}, y_k, \hat{x}_{2,k}, u_k, h \phi_2(e_{1,k+1})) \\
&- h l_1 \gamma \phi_1(e_{1,k+1}).
\end{split}
\tag{144}
\end{equation}
$$

It is inferred that

$$
\begin{equation}
\begin{split}
e_{1,k+1} \in (\mathbf{I}_d - h f_1(\cdot, y_k, u_k) + h \Phi(\cdot, y_k, \hat{x}_{2,k}, u_k, h \phi_2(\cdot)) \\
\qquad \qquad \qquad \qquad + h l_1 \gamma \phi_1(\cdot))^{-1} (y_{k-1} - y_k + e_{1,k}).
\end{split}
\tag{145}
\end{equation}
$$

Then, provided that the operators $x \mapsto -f_1(x, y_k, u_k) + l_1 \gamma \phi_1(x)$ and $x \mapsto \Phi(x, y_k, \hat{x}_{2,k}, u_k, h\phi_2(x))$ are maximal monotone (notice that $\phi_1$ and $\phi_2$ are strongly monotone for $q > \frac{1}{2}$), the inclusion in (145) is an equality and the right-hand side is a resolvent. Depending on the nonlinearities, the closed form of this resolvent may be impossible to obtain, in which case the numerical solutions described below are necessary. The observed variable $\hat{x}_{2,k+1}$ is obtained from (143a) and (145). Of course, semi-implicit methods can also be applied so that the computations are easier and the assumptions may be relaxed.
## 4.2 Higher-Order Differentiators

It is sometimes necessary to compute higher-order derivatives of $f$. A naive approach is to cascade several differentiators of the form (107) (which may provide reasonably good results for the second derivative [90]). A more efficient solution is to directly construct a differentiator of order $r-1$, $r \ge 2$, by enforcing an r-sliding mode. The arbitrary-order super-twisting differentiator (AO-STD) has the dynamics in non-recursive form [7, 81]

$$ 
\begin{aligned}
\dot{z}_i(t) &= z_{i+1}(t) - \gamma_i L \frac{\Gamma}{\Gamma+1} \left[ s_0(t) \right] \frac{\Gamma}{\Gamma+1}, \quad i = 0, \dots, r-1 \\
\dot{z}_r(t) &\in -\gamma_r L \operatorname{sgn}(s_0(t))
\end{aligned}
\qquad (146) $$

with $s_0(t) = z_0(t) - f(t)$, which is the counterpart of (107), while the counterpart of (108) writes as the recursive form

$$ 
\begin{aligned}
\dot{\sigma}_i(t) &= \sigma_{i+1}(t) - \gamma_i \left[ \sigma_0(t) \right] \frac{\Gamma}{\Gamma+1}, \quad i = 0, \dots, r-1, \\
\dot{\sigma}_r(t) &\in -\gamma_r \operatorname{sgn}(\sigma_0(t)) + f^{(r+1)}(t)
\end{aligned}
\qquad (147) $$

where $\sigma_i(t) := \frac{1}{L}(z_i(t) - f^{(i)}(t))$. Accuracy and robustness of these differentiators have been studied by exploiting the homogeneity of (147) [81, 49]. As shown in [84] and [83, 91], their discrete-time versions lend themselves to transformations which yield the same structures as in (127) or (112). In particular it is shown in [84, Table 2] that “classical” higher-order differentiators known as URED (uniform robust exact differentiator) [92], AO-STD as above [7, 81], HDD (homogeneous discrete-time differentiator), GHDD (generalized homogeneous discrete-time differentiator) [93], FDFF (first-order differentiator with first-order sliding-mode filtering) [5, 94], AO-FDFF (arbitrary-order differentiator with first-order sliding-mode filtering) [94, 95] yield, when implicitly or semi-implicitly discretized, the equation

$$ s_{0,k+1} = \mathcal{J}_h \mathbf{M}_h (-b_k), \qquad (148) $$

where $\mathbf{M}_h = \mathbf{M}_{\text{sing},h} + \mathbf{M}_{\text{set},h}$ with $\mathbf{M}_{\text{set},h}$ set-valued maximal monotone, $\mathbf{M}_{\text{sing},h}$ single-valued monotone, and $b_k = b(s_{0,k}, z_{1,k}, \dots, z_{n,k}, \Delta f_k)$, $\Delta f_k = f_{k+1} - f_k$.

Let us provide the implicit URED operators, which consist of the super-twisting operators with additional higher-order terms [84, Table 2, Equation (39)]. Its defining operators are

$$ 
\begin{aligned}
\mathbf{M}_{\text{sing},h}(w) &= \gamma_0 L^{\frac{1}{2}} \left( \left[ w \right]^{\frac{1}{2}} + \mu \left[ w \right]^{\frac{3}{2}} \right) + \frac{3}{2} h \mu^2 \gamma_1 L \left[ w \right]^2 \\
&\quad + 2h \gamma_1 L \mu w
\end{aligned}
\qquad (149a) $$

$$ \mathbf{M}_{\text{set},h}(w) = \frac{1}{2} h \gamma_1 L \operatorname{sgn}(w) \qquad (149b) $$

and the argument is

$$ b(s_{0,k}, z_{1,k}, \Delta f_k) = -s_{0,k} - h z_{1,k} + \Delta f_k, \qquad (149c) $$

where $\lambda_0, \lambda_1, \mu$ are parameters. Similarly to the super-twisting differentiator, the URED has an additional state

variable $x_1$, see, for example, [84, equation (6)]. The URED counterpart of (113) takes the form

$$ z_{1,k+1} = \frac{1}{h} (s_{0,k+1} - s_{0,k} + \Delta f_k + h \gamma_0 L^{\frac{1}{2}} \left[ s_{0,k+1} \right]^{\frac{1}{2}} + h L^{\frac{1}{2}} \mu \left[ s_{0,k+1} \right]^{\frac{3}{2}}). \qquad (150) $$

The semi-implicit version of this algorithm is derived as follows. First, $s_{0,k+1} = \mathcal{J}_h \mathbf{M}_{k,h} (-b_k)$ is computed with

$$ \mathbf{M}_{k,h}(w) = \left( \gamma_0 L^{\frac{1}{2}} \left( |s_{0,k}|^{\frac{1}{2}} + \mu |s_{0,k}|^{\frac{3}{2}} \right) + \frac{3}{2} h \mu^2 \gamma_1 L s_{0,k}^2 \right. \\ \left. + h \gamma_L |s_{0,k}| + \frac{1}{2} h \gamma_1 L \right) \operatorname{sgn}(w), \qquad (151) $$

which is the URED counterpart of (115). Let $\lambda_{k+1} \in \operatorname{sgn}(s_{0,k+1})$ be a selection of the set-valued sign function. Solving for $s_{0,k+1} \in \mathbf{N}_{[-1,1]}(\lambda_{k+1})$ gives

$$ s_{0,k} + h x_{1,k} - \Delta f_k - h \gamma_0 L^{\frac{1}{2}} \left( |s_{0,k}|^{\frac{1}{2}} + \mu |s_{0,k}|^{\frac{3}{2}} \right) \lambda_{k+1} - h^2 \gamma_1 L \left( 2\mu |s_{0,k}| + \frac{3}{2} \mu^2 s_{0,k}^2 + \frac{1}{2} \right) \lambda_{k+1} \in \mathbf{N}_{[-1,1]}(\lambda_{k+1}), \qquad (152) $$

so that $\lambda_{k+1} = \mathcal{J}_{\mathbf{M}_{k,h}} \left( -\frac{2b_k}{h^2 \gamma_1 L} \right)$ with

$$ \mathbf{M}_{k,h} = \frac{2}{h \lambda_1 L^{\frac{1}{2}}} \left( \gamma_0 \left( |s_{0,k}|^{\frac{1}{2}} + \mu |s_{0,k}|^{\frac{3}{2}} \right) \right. \\ \left. + h \gamma_1 L^{\frac{1}{2}} \left( 2\mu |s_{0,k}| + \frac{3}{2} \mu^2 s_{0,k}^2 + \frac{1}{2} \right) \right) \mathbf{I}_d + \mathbf{N}_{[-1,1]} \qquad (153) $$

Finally, we can recover

$$ z_{1,k+1} = z_{1,k} + h \gamma_1 L \left( 2\mu |s_{0,k}| + \frac{3}{2} \mu^2 s_{0,k}^2 + \frac{1}{2} \right) \lambda_{k+1}, \qquad (154) $$

which is the URED counterpart of (121). Thus, (112), (113), and (115), (121), and (148), (150), and (151), (154) make nonlinear difference equations with state vector $(s_{0,k}, z_{1,k})^T$. Some of these may be named higher-order proximal-point algorithms since inserting, for example, (113) into (112) yields a difference equation with $s_{0,k}$ and $s_{0,k-1}$ in its right-hand side. From the point of view of the resolvent calculation at step $k$, all terms depending only on measurements at $k$ are considered as being constants.

Similar developments can be carried out for the AO-STD in both implicit and semi-implicit cases [84, sections 5.1, 5.8]. For instance, the AO-STD of order 3 yields $s_{0,k+1} = \mathcal{J}_h \mathbf{M}_{k,h} (-b_k)$ with

$$ \mathbf{M}_{k,h}(w) = \sum_{l=0}^{2} h^l \gamma_l L^{\frac{l+1}{4}} \left[ w \right]^{\frac{3-l}{4}} + h^3 \gamma_3 L \operatorname{sgn}(w) \qquad (155) $$
## 4.2 Higher-Order Differentiators

It is sometimes necessary to compute higher-order derivatives of $f$. A naive approach is to cascade several differentiators of the form (107) (which may provide reasonably good results for the second derivative [90]). A more efficient solution is to directly construct a differentiator of order $r-1$, $r \ge 2$, by enforcing an r-sliding mode. The arbitrary-order super-twisting differentiator (AO-STD) has the dynamics in non-recursive form [7, 81]

$$ \begin{aligned} \dot{z}_i(t) &= z_{i+1}(t) - \gamma_i L \frac{\Gamma}{\Gamma+1} [s_0(t)]^{\frac{r-i}{r+1}}, \quad i = 0, \dots, r-1 \\ \dot{z}_r(t) &\in -\gamma_r L \operatorname{sgn}(s_0(t)) \end{aligned} \quad (146) $$

with $s_0(t) = z_0(t) - f(t)$, which is the counterpart of (107), while the counterpart of (108) writes as the recursive form

$$ \begin{aligned} \dot{\sigma}_i(t) &= \sigma_{i+1}(t) - \gamma_i [\sigma_0(t)]^{\frac{r-i}{r+1}}, \quad i = 0, \dots, r-1, \\ \dot{\sigma}_r(t) &\in -\gamma_r \operatorname{sgn}(\sigma_0(t)) + f^{(r+1)}(t) \end{aligned} \quad (147) $$

where $\sigma_i(t) := \frac{1}{L}(z_i(t) - f^{(i)}(t))$. Accuracy and robustness of these differentiators have been studied by exploiting the homogeneity of (147) [81, 49]. As shown in [84] and [83, 91], their discrete-time versions lend themselves to transformations which yield the same structures as in (127) or (112). In particular it is shown in [84, Table 2] that “classical” higher-order differentiators known as URED (uniform robust exact differentiator) [92], AO-STD as above [7, 81], HDD (homogeneous discrete-time differentiator), GHDD (generalized homogeneous discrete-time differentiator) [93], FDFF (first-order differentiator with first-order sliding-mode filtering) [5, 94], AO-FDFF (arbitrary-order differentiator with first-order sliding-mode filtering) [94, 95] yield, when implicitly or semi-implicitly discretized, the equation

$$ s_{0,k+1} = \mathcal{J}_h \mathbf{M}_h (-b_k), \quad (148) $$

where $\mathbf{M}_h = \mathbf{M}_{\text{sing},h} + \mathbf{M}_{\text{set},h}$ with $\mathbf{M}_{\text{set},h}$ set-valued maximal monotone, $\mathbf{M}_{\text{sing},h}$ single-valued monotone, and $b_k = b(s_{0,k}, z_{1,k}, \dots, z_{n,k}, \Delta f_k)$, $\Delta f_k = f_{k+1} - f_k$.

Let us provide the implicit URED operators, which consist of the super-twisting operators with additional higher-order terms [84, Table 2, Equation (39)]. Its defining operators are

$$ \mathbf{M}_{\text{sing},h}(w) = \gamma_0 L^{\frac{1}{2}} \left( [w]^{\frac{1}{2}} + \mu [w]^{\frac{3}{2}} \right) + \frac{3}{2} h \mu^2 \gamma_1 L [w]^2 + 2h \gamma_1 L \mu w \quad (149a) $$

$$ \mathbf{M}_{\text{set},h}(w) = \frac{1}{2} h \gamma_1 L \operatorname{sgn}(w) \quad (149b) $$

and the argument is

$$ b(s_{0,k}, z_{1,k}, \Delta f_k) = -s_{0,k} - h z_{1,k} + \Delta f_k, \quad (149c) $$

where $\lambda_0, \lambda_1, \mu$ are parameters. Similarly to the super-twisting differentiator, the URED has an additional state

variable $x_1$, see, for example, [84, equation (6)]. The URED counterpart of (113) takes the form

$$ z_{1,k+1} = \frac{1}{h} (s_{0,k+1} - s_{0,k} + \Delta f_k + h \gamma_0 L^{\frac{1}{2}} [s_{0,k+1}]^{\frac{1}{2}} + h L^{\frac{1}{2}} \mu [s_{0,k+1}]^{\frac{3}{2}}). \quad (150) $$

The semi-implicit version of this algorithm is derived as follows. First, $s_{0,k+1} = \mathcal{J}_h \mathbf{M}_{k,h} (-b_k)$ is computed with

$$ \mathbf{M}_{k,h}(w) = \left( \gamma_0 L^{\frac{1}{2}} \left( |s_{0,k}|^{\frac{1}{2}} + \mu |s_{0,k}|^{\frac{3}{2}} \right) + \frac{3}{2} h \mu^2 \gamma_1 L s_{0,k}^2 + h \gamma_L |s_{0,k}| + \frac{1}{2} h \gamma_1 L \right) \operatorname{sgn}(w), \quad (151) $$

which is the URED counterpart of (115). Let $\lambda_{k+1} \in \operatorname{sgn}(s_{0,k+1})$ be a selection of the set-valued sign function. Solving for $s_{0,k+1} \in \mathbf{N}_{[-1,1]}(\lambda_{k+1})$ gives

$$ s_{0,k} + h x_{1,k} - \Delta f_k - h \gamma_0 L^{\frac{1}{2}} \left( |s_{0,k}|^{\frac{1}{2}} + \mu |s_{0,k}|^{\frac{3}{2}} \right) \lambda_{k+1} - h^2 \gamma_1 L \left( 2\mu |s_{0,k}| + \frac{3}{2} \mu^2 s_{0,k}^2 + \frac{1}{2} \right) \lambda_{k+1} \in \mathbf{N}_{[-1,1]}(\lambda_{k+1}), \quad (152) $$

so that $\lambda_{k+1} = \mathcal{J}_{\mathbf{M}_{k,h}} \left(-\frac{2b_k}{h^2 \gamma_1 L}\right)$ with

$$ \mathbf{M}_{k,h} = \frac{2}{h \lambda_1 L^{\frac{1}{2}}} \left( \gamma_0 \left( |s_{0,k}|^{\frac{1}{2}} + \mu |s_{0,k}|^{\frac{3}{2}} \right) + h \gamma_1 L^{\frac{1}{2}} \left( 2\mu |s_{0,k}| + \frac{3}{2} \mu^2 s_{0,k}^2 + \frac{1}{2} \right) \right) \mathbf{I}_d + \mathbf{N}_{[-1,1]} \quad (153) $$

Finally, we can recover

$$ z_{1,k+1} = z_{1,k} + h \gamma_1 L \left( 2\mu |s_{0,k}| + \frac{3}{2} \mu^2 s_{0,k}^2 + \frac{1}{2} \right) \lambda_{k+1}, \quad (154) $$

which is the URED counterpart of (121). Thus, (112), (113), and (115), (121), and (148), (150), and (151), (154) make nonlinear difference equations with state vector $(s_{0,k}, z_{1,k})^T$. Some of these may be named higher-order proximal-point algorithms since inserting, for example, (113) into (112) yields a difference equation with $s_{0,k}$ and $s_{0,k-1}$ in its right-hand side. From the point of view of the resolvent calculation at step $k$, all terms depending only on measurements at $k$ are considered as being constants.

Similar developments can be carried out for the AO-STD in both implicit and semi-implicit cases [84, sections 5.1, 5.8]. For instance, the AO-STD of order 3 yields $s_{0,k+1} = \mathcal{J}_h \mathbf{M}_{k,h} (-b_k)$ with

$$ \mathbf{M}_{k,h}(w) = \sum_{l=0}^{2} h^l \gamma_l L^{\frac{l+1}{4}} [w]^{\frac{3-l}{4}} + h^3 \gamma_3 L \operatorname{sgn}(w) \quad (155) $$
for the implicit method, and

$$ \mathbf{M}_{k,h}(w) = \left( \sum_{l=0}^{2} h^l \gamma_l L^{\frac{l+1}{4}} |s_{0,k}|^{\frac{3-l}{4}} + h^3 \gamma_3 L \right) \operatorname{sgn}(w) \qquad (156) $$

for the semi-implicit method. The remaining state variables $z_{1,k+1}$, $z_{2,k+1}$, $z_{3,k+1}$ (which are the higher-order derivatives of $f$) are calculated similarly as for the super-twisting and the URED, yielding other higher-order proximal-point algorithms.

## 4.3 Relay polynomial controllers

As mentioned above, an interesting feature of higher-order sliding modes is the possibility for the state to attain the origin in finite time. There exist several classes of algorithms achieving such objective (see, e.g., the quasi-continuous or the nested controllers [96]). The following algorithm is one of the simplest. Consider a perturbed chain of integrators

$$ \begin{aligned} \dot{x}_i(t) &= x_{i+1}(t), \quad i = 1, \dots, n-1, \\ \dot{x}_n(t) &= u(t) + \delta(t). \end{aligned} \qquad (157) $$

The relay polynomial controller $u(t) = \hat{u}(x(t))$ with

$$ \dot{u}(x) \in -\gamma_n \operatorname{sgn} \left( [x_n]^n + \gamma_{n-1} [x_{n-1}]^{\frac{n}{2}} + \dots + \gamma_1 x_1 \right) \qquad (158) $$

establishes an n-sliding mode at the origin in finite time, provided that the perturbation $\delta(t)$ is uniformly bounded and that the gains $\gamma_i > 0$ are appropriately chosen [97, 76]. The forward-Euler discretization of the open-loop system is

$$ \begin{aligned} x_{i,k+1} &= x_{i,k} + h x_{i+1,k}, \quad i = 1, \dots, n-1, \\ x_{n,k+1} &= x_{n,k} + h u_k + \delta_k. \end{aligned} \qquad (159) $$

The backward-Euler discretization of (158) takes the implicit form $u_k \in -\gamma_n \mathbf{M}_{k,h}(u_k)$ with

$$ \mathbf{M}_{k,h}(w) = \operatorname{sgn} \left( [x_{n,k} + hw]^n + \gamma_{n-1} [x_{n-1,k} + hx_{n,k}]^{\frac{n}{2}} + \dots + \gamma_1 (x_{1,k} + hx_{2,k}) \right). \qquad (160) $$

It is not difficult to see that, for every $x_k$, $\mathbf{M}_{k,h}$ is maximal monotone, hence a unique solution $u_k(x_k) = (I_d + \gamma_n \mathbf{M}_{k,h})^{-1}(0)$ exists. In fact, it is possible to write an analytic expression for the solution. Note that $u_k \in -\gamma_n \mathbf{M}_{k,h}(u_k)$ is equivalent to

$$ u_k \in -\gamma_n \operatorname{sgn}(u_k - \bar{u}(x_k)), \qquad (161) $$

where

$$ \bar{u}(x_k) = -\frac{1}{h} \left( \begin{aligned} & \left[ \gamma_{n-1} [x_{n-1,k} + hx_{n,k}]^{\frac{\alpha}{2}} + \dots \right]^{\frac{n}{2}} \\ & + \gamma_1 [x_{1,k} + hx_{2,k}]^{\frac{\alpha}{n}} \end{aligned} \right) \qquad (162) $$

The control action is thus $u_k = \bar{u}(x_k) - \mathcal{J}_{\gamma_n} \operatorname{sgn}(\bar{u}(x_k))$, that is,

$$ u_k = \min\{|\bar{u}(x_k)|, \gamma_n\} \operatorname{sgn}(\bar{u}(x_k)), \qquad (163) $$

which once again is a form of minimum-operator formalism as (13c). Consider now the ZOH discrete model of the chain of integrators,

$$ x_{i,k+1} = x_{i,k} + \sum_{j=1}^{n-i} \frac{h^j}{j!} x_{i+j,k} + \frac{h^{n+1-i}}{(n+1-i)!} (u_k + \bar{\delta}_{i,k}), \qquad (164) $$

where $i = 1, \dots, n$ and $\bar{\delta}_i$ are (possibly unmatched) disturbances resulting from $\delta(t)$. Consider also the inclusion $u_k \in -\gamma_n \mathbf{M}_{k,h}(u_k)$, where

$$ \mathbf{M}_{k,h}(w) = \operatorname{sgn} \left( [x_{n,k} + hw]^n + \gamma_{n-1} \left[ x_{n-1,k} + h x_{n,k} + \frac{h^2}{2!} w \right]^{\frac{n}{2}} + \dots + \gamma_1 \left( x_{1,k} + \sum_{j=1}^{n-1} \frac{h^j}{j!} x_{1+j,k} + \frac{h^n}{n!} w \right) \right). \qquad (165) $$

For each fixed $x_k$, the operator $\mathbf{M}_{k,h}$ is again maximal monotone, so the resolvent $(I_d + \mathbf{M}_{k,h})^{-1}$ has domain $\mathbb{R}$ and is single-valued. While an analytic expression is no longer available, it is still possible to compute it numerically by using the techniques described in Section “6”.

## 4.4 Stability Analysis

The closed-loop systems (112), (113), or (114), or (148), (150), or its semi-implicit version in (151), (154), are higher-order proximal-point algorithms which do not fit within the robust proximal-point algorithm in (57) in “Sidebar VII”. Boundedness of the closed-loop state with (128) is proved only in the absence of perturbations in [34]. As shown in [86], the algorithm is not robust with respect to an unbounded perturbation. A modification to (128) is proposed in [86] that enhances robustness. Some results concerning higher-order differentiators of the implicit AO-STD type can be found in [84]. Some properties of the implicit higher-order differentiators (sliding surface invariance, accuracy, finite convergence) are shown in [84, 91]. The stability of the semi-implicit super-twisting controller is analysed in [35], see also the semi-implicit discretization [31, Equation (60)] of the super-twisting velocity observer [98].

## 4.5 Homogeneity and maximal monotonicity

Homogeneity is a powerful tool for control design [99]. It plays an important role in all of the higher-order sliding-mode algorithms presented above. While homogeneity has no clear discrete-time counterpart, it is possible to construct consistent discretizations by keeping track of the maximal monotone operators that appear
for the implicit method, and

$$ \mathbf{M}_{k,h}(w) = \left( \sum_{l=0}^{2} h^l \gamma_l L^{\frac{l+1}{4}} |s_{0,k}|^{\frac{3-l}{4}} + h^3 \gamma_3 L \right) \operatorname{sgn}(w) \quad (156) $$

for the semi-implicit method. The remaining state variables $z_{1,k+1}$, $z_{2,k+1}$, $z_{3,k+1}$ (which are the higher-order derivatives of $f$) are calculated similarly as for the super-twisting and the URED, yielding other higher-order proximal-point algorithms.

## 4.3 Relay polynomial controllers

As mentioned above, an interesting feature of higher-order sliding modes is the possibility for the state to attain the origin in finite time. There exist several classes of algorithms achieving such objective (see, e.g., the quasi-continuous or the nested controllers [96]). The following algorithm is one of the simplest. Consider a perturbed chain of integrators

$$ \begin{aligned} \dot{x}_i(t) &= x_{i+1}(t), \quad i = 1, \dots, n-1, \\ \dot{x}_n(t) &= u(t) + \delta(t). \end{aligned} \quad (157) $$

The relay polynomial controller $u(t) = \hat{u}(x(t))$ with

$$ \dot{u}(x) \in -\gamma_n \operatorname{sgn} \left( \left[x_n\right]^n + \gamma_{n-1} \left[x_{n-1}\right]^{\frac{n}{2}} + \dots + \gamma_1 x_1 \right) \quad (158) $$

establishes an n-sliding mode at the origin in finite time, provided that the perturbation $\delta(t)$ is uniformly bounded and that the gains $\gamma_i > 0$ are appropriately chosen [97, 76]. The forward-Euler discretization of the open-loop system is

$$ \begin{aligned} x_{i,k+1} &= x_{i,k} + h x_{i+1,k}, \quad i = 1, \dots, n-1, \\ x_{n,k+1} &= x_{n,k} + h u_k + \delta_k. \end{aligned} \quad (159) $$

The backward-Euler discretization of (158) takes the implicit form $u_k \in -\gamma_n \mathbf{M}_{k,h}(u_k)$ with

$$ \mathbf{M}_{k,h}(w) = \operatorname{sgn} \left( \left[x_{n,k} + hw\right]^n + \gamma_{n-1} \left[x_{n-1,k} + hx_{n,k}\right]^{\frac{n}{2}} + \dots + \gamma_1 (x_{1,k} + hx_{2,k}) \right). \quad (160) $$

It is not difficult to see that, for every $x_k$, $\mathbf{M}_{k,h}$ is maximal monotone, hence a unique solution $u_k(x_k) = (I_d + \gamma_n \mathbf{M}_{k,h})^{-1}(0)$ exists. In fact, it is possible to write an analytic expression for the solution. Note that $u_k \in -\gamma_n \mathbf{M}_{k,h}(u_k)$ is equivalent to

$$ u_k \in -\gamma_n \operatorname{sgn}(u_k - \bar{u}(x_k)), \quad (161) $$

where

$$ \bar{u}(x_k) = -\frac{1}{h} \left( \left[ \gamma_{n-1} \left[x_{n-1,k} + hx_{n,k}\right]^{\frac{\alpha}{2}} + \dots + \gamma_1 \left[x_{1,k} + hx_{2,k}\right]^{\frac{\alpha}{n}} \right]^{\frac{1}{\alpha}} + x_{n,k} \right). \quad (162) $$

The control action is thus $u_k = \bar{u}(x_k) - \mathcal{J}_{\gamma_n} \operatorname{sgn}(\bar{u}(x_k))$, that is,

$$ u_k = \min\{|\bar{u}(x_k)|, \gamma_n\} \operatorname{sgn}(\bar{u}(x_k)), \quad (163) $$

which once again is a form of minimum-operator formalism as (13c). Consider now the ZOH discrete model of the chain of integrators,

$$ x_{i,k+1} = x_{i,k} + \sum_{j=1}^{n-i} \frac{h^j}{j!} x_{i+j,k} + \frac{h^{n+1-i}}{(n+1-i)!} (u_k + \bar{\delta}_{i,k}), \quad (164) $$

where $i = 1, \dots, n$ and $\bar{\delta}_i$ are (possibly unmatched) disturbances resulting from $\delta(t)$. Consider also the inclusion $u_k \in -\gamma_n \mathbf{M}_{k,h}(u_k)$, where

$$ \begin{aligned} \mathbf{M}_{k,h}(w) = \operatorname{sgn} & \left( \left[x_{n,k} + hw\right]^n \right. \\ & + \gamma_{n-1} \left[ x_{n-1,k} + h x_{n,k} + \frac{h^2}{2!} w \right]^{\frac{n}{2}} + \\ & \left. \dots + \gamma_1 \left( x_{1,k} + \sum_{j=1}^{n-1} \frac{h^j}{j!} x_{1+j,k} + \frac{h^n}{n!} w \right) \right). \end{aligned} \quad (165) $$

For each fixed $x_k$, the operator $\mathbf{M}_{k,h}$ is again maximal monotone, so the resolvent $(I_d + \mathbf{M}_{k,h})^{-1}$ has domain $\mathbb{R}$ and is single-valued. While an analytic expression is no longer available, it is still possible to compute it numerically by using the techniques described in Section “6”.

## 4.4 Stability Analysis

The closed-loop systems (112), (113), or (114), or (148), (150), or its semi-implicit version in (151), (154), are higher-order proximal-point algorithms which do not fit within the robust proximal-point algorithm in (57) in “Sidebar VII”. Boundedness of the closed-loop state with (128) is proved only in the absence of perturbations in [34]. As shown in [86], the algorithm is not robust with respect to an unbounded perturbation. A modification to (128) is proposed in [86] that enhances robustness. Some results concerning higher-order differentiators of the implicit AO-STD type can be found in [84]. Some properties of the implicit higher-order differentiators (sliding surface invariance, accuracy, finite convergence) are shown in [84, 91]. The stability of the semi-implicit super-twisting controller is analysed in [35], see also the semi-implicit discretization [31, Equation (60)] of the super-twisting velocity observer [98].

## 4.5 Homogeneity and maximal monotonicity

Homogeneity is a powerful tool for control design [99]. It plays an important role in all of the higher-order sliding-mode algorithms presented above. While homogeneity has no clear discrete-time counterpart, it is possible to construct consistent discretizations by keeping track of the maximal monotone operators that appear
Backward-Euler discretization of sliding-mode systems gives rise to new variants of proximal-point algorithms.

in the homogeneous vector fields. Consistent discretization of finite-time convergent systems have been presented in [19, Definitions 2.1, 2.2]. As an example, let us consider the systems analysed in [19, Section 6.2]. The first analysis concerns the discretization of a generalized homogeneous closed-loop system with homogeneity degree $d = 1$ (see [99] for definitions). After a suitable transformation into another homogeneous system, the consistent discrete-time dynamics

$$x_{k+1} = x_k - h\|x_k\|(I_n - \tilde{A})x_{k+1} \quad (166)$$

is obtained [19, Section 6.2, Equation (6.8)]. Similar to the discrete dynamics presented in the subsections above, this dynamics is given by the resolvent of a maximal monotone operator. Indeed, these are equivalent to

$$x_{k+1} = (I_d + h\|x_k\|(I_n - \tilde{A}))^{-1}(x_k), \quad k \ge 0, \quad (167)$$

with $\tilde{A}^\top P + P\tilde{A} = 0$ for some $P = P^\top \succ 0$ [19, Definition 2.1]. Let us define $R = R^\top \succ 0$, $R^2 = P$, then we obtain $R^{-1}(\tilde{A} - I_n)^\top R + R(\tilde{A} - I_n)R^{-1} = -2I_n \prec 0$, from where it follows that $R(I_n - \tilde{A})R^{-1} \succ 0$. Using the new state variable $\bar{x}_k = Rx_k$, the system dynamics takes the form

$$\bar{x}_{k+1} = \bar{x}_k - h\|R^{-1}\bar{x}_k\| R(I_n - \tilde{A})R^{-1}\bar{x}_{k+1}, \quad (168)$$

which is equivalent to

$$\bar{x}_{k+1} = (I_d + h\|R^{-1}\bar{x}_k\| R(I_n - \tilde{A})R^{-1})^{-1}(\bar{x}_k), \quad k \ge 0. \quad (169)$$

Therefore, at each step $k$ the scheme is updated by calculating the resolvent $\mathcal{J}_{M_k}$ of the maximal monotone operator $M_k: w \mapsto \|R^{-1}\bar{x}_k\| R(I_n - \tilde{A})R^{-1}w$. The proximal algorithm that is obtained is $\bar{x}_{k+1} = \mathcal{J}_{M_k}(\bar{x}_k)$, it may be named a time-varying proximal-point algorithm.

A second case, with homogeneity degree $d = -1$, is treated in [19, Section 6.2, Equation (6.4)]. It yields the consistent discretization

$$x_{k+1} \in x_k - h(I_n - \tilde{A})\tilde{F}(x_{k+1}), \quad (170)$$

which is equivalent to

$$x_{k+1} = (I_d + h(I_n - \tilde{A})\tilde{F})^{-1}(x_k) \quad (171)$$

with

$$\tilde{F}(x) = \begin{cases} \frac{x}{\|x\|_P} & \text{if } x \ne 0 \\ B_P & \text{if } x = 0, \end{cases} \quad (172)$$

and $\|x\|_P = \sqrt{x^\top P x}$. Using again the state variable $\bar{x}_k = Rx_k$ one obtains

$$\begin{aligned} \bar{x}_{k+1} &\in \bar{x}_k - hR(I_n - \tilde{A})R^{-1}\tilde{F}(x_{k+1}) \\ &= \bar{x}_k - hR(I_n - \tilde{A})R^{-1}\partial\|\bar{x}_{k+1}\|_2. \end{aligned} \quad (173)$$

For the same reasons given in the first case, $d = 1$, it follows that $\Lambda = R(I_n - \tilde{A})R^{-1} \succ 0$. However $\Lambda \ne \Lambda^\top$ in general. Now, let us proceed as in the proof of [19, Theorem 4.1]. For all $\bar{x}_k \ne 0$,

$$\begin{aligned} \bar{x}_k^\top \bar{x}_k &= \bar{x}_{k+1}^\top \bar{x}_{k+1} + 2h \bar{x}_{k+1}^\top \Lambda \partial \| \bar{x}_{k+1} \|_2 \\ &\quad + h^2 (\Lambda \partial \| z_{k+1} \|_2)^\top (\Lambda \partial \| z_{k+1} \|_2) \end{aligned} \quad (174)$$

$$\begin{aligned} &= \bar{x}_{k+1}^\top \bar{x}_{k+1} + 2h \frac{\bar{x}_{k+1}^\top \Lambda \bar{x}_{k+1}}{\|\bar{x}_{k+1}\|_2} \\ &\quad + h^2 (\Lambda \partial \| \bar{x}_{k+1} \|_2)^\top \Lambda \partial \| \bar{x}_{k+1} \|_2 \end{aligned} \quad (175)$$

$$> \bar{x}_{k+1}^\top \bar{x}_{k+1} + h^2 (\Lambda \partial \| \bar{x}_{k+1} \|_2)^\top \Lambda \partial \| \bar{x}_{k+1} \|_2 \quad (176)$$

$$= \bar{x}_{k+1}^\top \bar{x}_{k+1} + h^2 \frac{\bar{x}_{k+1}^\top \Lambda^\top \Lambda \bar{x}_{k+1}}{\|\bar{x}_{k+1}\|_2^2}. \quad (177)$$

Since $\Lambda^\top \Lambda \succ 0$,

$$\inf_{\xi \ne 0} h^2 \frac{\xi^\top \Lambda^\top \Lambda \xi}{\|\xi\|_2^2} > 0, \quad (178)$$

for all $h > 0$, from where the convergence of $\bar{x}_k$ to zero in a finite number of steps is inferred. Thus, in spite of the fact that neither the operator inside (171) nor inside $\bar{x}_{k+1} = (I_d + h\Lambda\partial\|\cdot\|_2)^{-1}(\bar{x}_k)$ are maximal monotone in general, finite-time convergence of this algorithm is guaranteed. Moreover the calculation is made in [19, Section 6.2]. As a side note, mark that (173) is equivalent to

$$\bar{x}_{k+1} \in \bar{x}_k - h(\Lambda_{\text{sym}} + \Lambda_{\text{sk}})\partial\|\bar{x}_{k+1}\|_2, \quad (179)$$

with $\Lambda_{\text{sym}} = \frac{1}{2}(\Lambda + \Lambda^\top)$, $\Lambda_{\text{sk}} = \frac{1}{2}(\Lambda - \Lambda^\top)$. Therefore, Equation (173) is interpreted as a discrete-time Hamiltonian system

$$\bar{x}_{k+1} \in \bar{x}_k + h(J - R)\partial H(\bar{x}_{k+1}) \quad (180)$$

with $J = -\Lambda_{\text{sk}}$, $R = \Lambda_{\text{sym}} \succ 0$, and Hamiltonian function $H(x) = \|x\|_2$. This is the implicit discretization of the Hamiltonian dynamics

$$\dot{x}(t) \in (J - R)\partial\|x(t)\|_2. \quad (181)$$

Thus, Hamiltonian proximal-point algorithms with resolvents $\mathcal{J}_{hM} = (I_d - h(J - R)\partial\|\cdot\|)^{-1}$, $M = -(J - R)\partial\|\cdot\|$ may be introduced. More details on the discretization of Hamiltonian systems can be found in [100].
Backward-Euler discretization of sliding-mode systems gives rise to new variants of proximal-point algorithms.

in the homogeneous vector fields. Consistent discretization of finite-time convergent systems have been presented in [19, Definitions 2.1, 2.2]. As an example, let us consider the systems analysed in [19, Section 6.2]. The first analysis concerns the discretization of a generalized homogeneous closed-loop system with homogeneity degree $d = 1$ (see [99] for definitions). After a suitable transformation into another homogeneous system, the consistent discrete-time dynamics

$$x_{k+1} = x_k - h\|x_k\|(I_n - \tilde{A})x_{k+1} \quad (166)$$

is obtained [19, Section 6.2, Equation (6.8)]. Similar to the discrete dynamics presented in the subsections above, this dynamics is given by the resolvent of a maximal monotone operator. Indeed, these are equivalent to

$$x_{k+1} = (I_d + h\|x_k\|(I_n - \tilde{A}))^{-1}(x_k), \quad k \ge 0, \quad (167)$$

with $\tilde{A}^\top P + P\tilde{A} = 0$ for some $P = P^\top \succ 0$ [19, Definition 2.1]. Let us define $R = R^\top \succ 0$, $R^2 = P$, then we obtain $R^{-1}(\tilde{A} - I_n)^\top R + R(\tilde{A} - I_n)R^{-1} = -2I_n \prec 0$, from where it follows that $R(I_n - \tilde{A})R^{-1} \succ 0$. Using the new state variable $\bar{x}_k = Rx_k$, the system dynamics takes the form

$$\bar{x}_{k+1} = \bar{x}_k - h\|R^{-1}\bar{x}_k\| R(I_n - \tilde{A})R^{-1}\bar{x}_{k+1}, \quad (168)$$

which is equivalent to

$$\bar{x}_{k+1} = (I_d + h\|R^{-1}\bar{x}_k\| R(I_n - \tilde{A})R^{-1})^{-1}(\bar{x}_k), \quad k \ge 0. \quad (169)$$

Therefore, at each step $k$ the scheme is updated by calculating the resolvent $\mathcal{J}_{M_k}$ of the maximal monotone operator $M_k: w \mapsto \|R^{-1}\bar{x}_k\| R(I_n - \tilde{A})R^{-1}w$. The proximal algorithm that is obtained is $\bar{x}_{k+1} = \mathcal{J}_{M_k}(\bar{x}_k)$, it may be named a time-varying proximal-point algorithm.

A second case, with homogeneity degree $d = -1$, is treated in [19, Section 6.2, Equation (6.4)]. It yields the consistent discretization

$$x_{k+1} \in x_k - h(I_n - \tilde{A})\tilde{F}(x_{k+1}), \quad (170)$$

which is equivalent to

$$x_{k+1} = (I_d + h(I_n - \tilde{A})\tilde{F})^{-1}(x_k) \quad (171)$$

with

$$\tilde{F}(x) = \begin{cases} \frac{x}{\|x\|_P} & \text{if } x \neq 0 \\ B_P & \text{if } x = 0, \end{cases} \quad (172)$$

and $\|x\|_P = \sqrt{x^\top P x}$. Using again the state variable $\bar{x}_k = Rx_k$ one obtains

$$\begin{aligned} \bar{x}_{k+1} &\in \bar{x}_k - hR(I_n - \tilde{A})R^{-1}\tilde{F}(x_{k+1}) \\ &= \bar{x}_k - hR(I_n - \tilde{A})R^{-1}\partial\|\bar{x}_{k+1}\|_2. \end{aligned} \quad (173)$$

For the same reasons given in the first case, $d = 1$, it follows that $\Lambda = R(I_n - \tilde{A})R^{-1} \succ 0$. However $\Lambda \neq \Lambda^\top$ in general. Now, let us proceed as in the proof of [19, Theorem 4.1]. For all $\bar{x}_k \neq 0$,

$$\begin{aligned} \bar{x}_k^\top \bar{x}_k &= \bar{x}_{k+1}^\top \bar{x}_{k+1} + 2h \bar{x}_{k+1}^\top \Lambda \partial \|\bar{x}_{k+1}\|_2 \\ &\quad + h^2 (\Lambda \partial \|z_{k+1}\|_2)^\top (\Lambda \partial \|z_{k+1}\|_2) \end{aligned} \quad (174)$$

$$\begin{aligned} &= \bar{x}_{k+1}^\top \bar{x}_{k+1} + 2h \frac{\bar{x}_{k+1}^\top \Lambda \bar{x}_{k+1}}{\|\bar{x}_{k+1}\|_2} \\ &\quad + h^2 (\Lambda \partial \|\bar{x}_{k+1}\|_2)^\top \Lambda \partial \|\bar{x}_{k+1}\|_2 \end{aligned} \quad (175)$$

$$> \bar{x}_{k+1}^\top \bar{x}_{k+1} + h^2 (\Lambda \partial \|\bar{x}_{k+1}\|_2)^\top \Lambda \partial \|\bar{x}_{k+1}\|_2 \quad (176)$$

$$= \bar{x}_{k+1}^\top \bar{x}_{k+1} + h^2 \frac{\bar{x}_{k+1}^\top \Lambda^\top \Lambda \bar{x}_{k+1}}{\|\bar{x}_{k+1}\|_2^2}. \quad (177)$$

Since $\Lambda^\top \Lambda \succ 0$,

$$\inf_{\xi \neq 0} h^2 \frac{\xi^\top \Lambda^\top \Lambda \xi}{\|\xi\|_2^2} > 0, \quad (178)$$

for all $h > 0$, from where the convergence of $\bar{x}_k$ to zero in a finite number of steps is inferred. Thus, in spite of the fact that neither the operator inside (171) nor inside $\bar{x}_{k+1} = (I_d + h\Lambda\partial\|\cdot\|_2)^{-1}(\bar{x}_k)$ are maximal monotone in general, finite-time convergence of this algorithm is guaranteed. Moreover the calculation is made in [19, Section 6.2]. As a side note, mark that (173) is equivalent to

$$\bar{x}_{k+1} \in \bar{x}_k - h(\Lambda_{\text{sym}} + \Lambda_{\text{sk}})\partial\|\bar{x}_{k+1}\|_2, \quad (179)$$

with $\Lambda_{\text{sym}} = \frac{1}{2}(\Lambda + \Lambda^\top)$, $\Lambda_{\text{sk}} = \frac{1}{2}(\Lambda - \Lambda^\top)$. Therefore, Equation (173) is interpreted as a discrete-time Hamiltonian system

$$\bar{x}_{k+1} \in \bar{x}_k + h(J - R)\partial H(\bar{x}_{k+1}) \quad (180)$$

with $J = -\Lambda_{\text{sk}}$, $R = \Lambda_{\text{sym}} \succ 0$, and Hamiltonian function $H(x) = \|x\|_2$. This is the implicit discretization of the Hamiltonian dynamics

$$\dot{x}(t) \in (J - R)\partial \|x(t)\|_2. \quad (181)$$

Thus, Hamiltonian proximal-point algorithms with resolvents $\mathcal{J}_{hM} = (I_d - h(J - R)\partial\|\cdot\|)^{-1}$, $M = -(J - R)\partial\|\cdot\|$ may be introduced. More details on the discretization of Hamiltonian systems can be found in [100].
# 5 FURTHER IMPLICIT ALGORITHMS

## 5.1 Other Classes of Set-Valued Algorithms

The implicit discretization for first-order sliding-mode controllers, super-twisting controller, observer and differentiator, and higher-order differentiators, has been studied in the foregoing sections. The implicit discretization of other classes of algorithms has been analysed in the literature: terminal sliding-mode control for first and second order systems [41, 101], twisting controller [102], proxy-based sliding-mode [32, 94, 103], homogeneous differentiators and observers [38, 39], super-twisting velocity observer [31], other types of differentiators [94, 95, 104, 105, 106, 107], the differentiator proposed in [5] discretized in [84], sliding-mode control of nonsmooth actuators [36], homogeneous systems and controllers [19, 108], nested controllers for second-order systems with unmatched perturbations [44], first-order controllers applied to Lagrangian nonlinear systems with parameter uncertainties [42] and to linear hyperbolic infinite dimensional systems [109], fixed-time or prescribed-time schemes [19, 110, 111, 101], adaptive sliding-mode algorithms [112, 113, 114], sliding-mode control with reaching law [115], fractional-order differentiators [116, 117], multivariable super-twisting [118], modified super-twisting using minimum-operator-based algorithms (similar to (13c)) [119].

The articles mentioned above study implicit or semi-implicit schemes. The solvability of the one-step nonsmooth problem (do the generalized equations possess a solution to advance the scheme from step k to step k + 1?) is usually analyzed. Some also analyze the boundedness/stability properties, and only a few of them analyse the discrete solutions convergence towards the solutions of the continuous-time closed-loop. The analysis made in foregoing sections focuses on the first step, which is necessary for the implementation.

## 5.2 Limitations and Modifications of Implicit Algorithms

The fact that implicit and semi-implicit discretization methods supersede the explicit one is widely accepted and validated (see section "7"). Nevertheless, open issues have emerged which show that room exists for significant improvements. The implicit method has been studied in [102] for the twisting scheme and in [34] for the super-twisting one. However, the straightforward implicit discretization of the twisting scheme (i.e., the implicit or backward-Euler emulation) does not guarantee finite-time convergence to the origin [120]: a modification of the emulated implicit method is proposed in [102] which secures it. As noticed in [86], the implicit emulation of the super-twisting scheme (127), (128) does not yield acceptable robustness (a fact also noticed experimentally in [121], see also [122] for the influence of the plant discretization). Again, a suitable modification of

the implicit emulation allows us to enhance it. Also, implicit emulation discretization of some higher-order differentiators may still undergo numerical chattering, and suitable modifications have to be made to avoid it [123]. Finally, finite-time stability may be lost after an implicit emulation discretization has been performed, and may be recovered if a suitable state transformation is applied before implicitly discretizing [19].

# 6 COMPUTATIONAL ISSUES

The analysis of closed-loop systems with implicit discrete-time controllers and the computation of output differentiators, relies on computing proximal maps to convex functions, or more generally, resolvents of maximal monotone operators (in order to solve the associated generalized equations). Figure 5 depicts the steps to follow in order to compute the resolvent of the sign multifunction. However, as seen in foregoing sections, some algorithms yield nontrivial resolvents. In general, finding closed-form expressions for such proximal maps is a hard task, so that the computation of the control law has to rely on numerical methods from convex optimization. Moreover, the numerical method used to calculate the resolvent can influence significantly the closed-loop behavior in some applications [40, 83, 90, 124].

## 6.1 Solving Generalized Equations with Proximal Operators

The computation of the variable of interest (control input or differentiator output) is reduced to solving either the generalized equation

$$ \begin{gathered} \zeta_{k+1} = \varphi(x_k) + \rho\eta(x_k) \\ -\eta(x_k) \in \mathbf{M}(\zeta_{k+1}) \end{gathered} \qquad (182a) $$

with unknown $\zeta_{k+1}$, or the generalized equation

$$ \begin{gathered} -\eta(x_k) = \frac{1}{\rho}\varphi(x_k) - \frac{1}{\rho}\zeta_{k+1}, \\ \zeta_{k+1} \in \mathbf{M}^{-1}(-\eta(x_k)) \end{gathered} \qquad (182b) $$

with unknown $\eta(x_k)$, where $\varphi: \mathbb{R}^n \to \mathbb{R}^m$ is a known function of the state, $\rho > 0$, and $\mathbf{M}: \mathbb{R}^m \to \mathbb{R}^m$ is maximal monotone.

As an example, note that equation (71) has this form with $\zeta_{k+1} = \tilde{s}'_{k+1}$, $\eta(x_k) = u_{sv}(x_k)$,

$$ \mathbf{M} = \partial \|\cdot\|_1 \circ \tilde{K}(\tilde{B}^\top \tilde{P}^{-1} \tilde{B})^{\frac{1}{2}}, \qquad (183) $$

and $\varphi(x_k) = \tilde{L}x_k$. Likewise, equation (95) with $\zeta_{k+1} = \tilde{s}'_{k+1}$, $\eta(x_k) = u_{sv}(x_k)$, $\mathbf{M} = \partial \|\cdot\|_1 \circ \tilde{K}^{\frac{1}{2}}$, and $\varphi(x_k) = \tilde{K}^{-\frac{1}{2}} s_k$ belongs to this class. Equations (124)–(130) also have this form with $\zeta_{k+1} = \tilde{x}_{0,k+1}$, $\eta(x_k) = u_k(x_{0,k}, \nu_k)$, $\mathbf{M}$ as in (111), $\varphi(x_k) = x_{0,k} + h\nu_k$. Equations (148)–(149) are an instance of (182) with $\zeta_{k+1} = s_{0,k+1}$, $\mathbf{M} = \mathbf{M}_{\text{sing}} + \mathbf{M}_{\text{set},h}$, $\varphi(x_k) = -b_k$. The unknown $\zeta_{k+1}$
# 5 FURTHER IMPLICIT ALGORITHMS

## 5.1 Other Classes of Set-Valued Algorithms

The implicit discretization for first-order sliding-mode controllers, super-twisting controller, observer and differentiator, and higher-order differentiators, has been studied in the foregoing sections. The implicit discretization of other classes of algorithms has been analysed in the literature: terminal sliding-mode control for first and second order systems [41, 101], twisting controller [102], proxy-based sliding-mode [32, 94, 103], homogeneous differentiators and observers [38, 39], super-twisting velocity observer [31], other types of differentiators [94, 95, 104, 105, 106, 107], the differentiator proposed in [5] discretized in [84], sliding-mode control of nonsmooth actuators [36], homogeneous systems and controllers [19, 108], nested controllers for second-order systems with unmatched perturbations [44], first-order controllers applied to Lagrangian nonlinear systems with parameter uncertainties [42] and to linear hyperbolic infinite dimensional systems [109], fixed-time or prescribed-time schemes [19, 110, 111, 101], adaptive sliding-mode algorithms [112, 113, 114], sliding-mode control with reaching law [115], fractional-order differentiators [116, 117], multivariable super-twisting [118], modified super-twisting using minimum-operator-based algorithms (similar to (13c)) [119].

The articles mentioned above study implicit or semi-implicit schemes. The solvability of the one-step nonsmooth problem (do the generalized equations possess a solution to advance the scheme from step $k$ to step $k+1$?) is usually analyzed. Some also analyze the boundedness/stability properties, and only a few of them analyse the discrete solutions convergence towards the solutions of the continuous-time closed-loop. The analysis made in foregoing sections focuses on the first step, which is necessary for the implementation.

## 5.2 Limitations and Modifications of Implicit Algorithms

The fact that implicit and semi-implicit discretization methods supersede the explicit one is widely accepted and validated (see section "7"). Nevertheless, open issues have emerged which show that room exists for significant improvements. The implicit method has been studied in [102] for the twisting scheme and in [34] for the super-twisting one. However, the straightforward implicit discretization of the twisting scheme (i.e., the implicit or backward-Euler emulation) does not guarantee finite-time convergence to the origin [120]: a modification of the emulated implicit method is proposed in [102] which secures it. As noticed in [86], the implicit emulation of the super-twisting scheme (127), (128) does not yield acceptable robustness (a fact also noticed experimentally in [121], see also [122] for the influence of the plant discretization). Again, a suitable modification of

the implicit emulation allows us to enhance it. Also, implicit emulation discretization of some higher-order differentiators may still undergo numerical chattering, and suitable modifications have to be made to avoid it [123]. Finally, finite-time stability may be lost after an implicit emulation discretization has been performed, and may be recovered if a suitable state transformation is applied before implicitly discretizing [19].

# 6 COMPUTATIONAL ISSUES

The analysis of closed-loop systems with implicit discrete-time controllers and the computation of output differentiators, relies on computing proximal maps to convex functions, or more generally, resolvents of maximal monotone operators (in order to solve the associated generalized equations). Figure 5 depicts the steps to follow in order to compute the resolvent of the sign multifunction. However, as seen in foregoing sections, some algorithms yield nontrivial resolvents. In general, finding closed-form expressions for such proximal maps is a hard task, so that the computation of the control law has to rely on numerical methods from convex optimization. Moreover, the numerical method used to calculate the resolvent can influence significantly the closed-loop behavior in some applications [40, 83, 90, 124].

## 6.1 Solving Generalized Equations with Proximal Operators

The computation of the variable of interest (control input or differentiator output) is reduced to solving either the generalized equation

$$ \begin{gathered} \zeta_{k+1} = \varphi(x_k) + \rho\eta(x_k) \\ -\eta(x_k) \in \mathbf{M}(\zeta_{k+1}) \end{gathered} \tag{182a} $$

with unknown $\zeta_{k+1}$, or the generalized equation

$$ \begin{gathered} -\eta(x_k) = \frac{1}{\rho}\varphi(x_k) - \frac{1}{\rho}\zeta_{k+1}, \qquad (182b) \\ \zeta_{k+1} \in \mathbf{M}^{-1}(-\eta(x_k)) \end{gathered} $$

with unknown $\eta(x_k)$, where $\varphi: \mathbb{R}^n \to \mathbb{R}^m$ is a known function of the state, $\rho > 0$, and $\mathbf{M}: \mathbb{R}^m \to \mathbb{R}^m$ is maximal monotone.

As an example, note that equation (71) has this form with $\zeta_{k+1} = \tilde{s}'_{k+1}$, $\eta(x_k) = u_{sv}(x_k)$,

$$ \mathbf{M} = \partial \|\cdot\|_1 \circ \tilde{K}(\tilde{B}^\top \tilde{P}^{-1} \tilde{B})^{\frac{1}{2}}, \tag{183} $$

and $\varphi(x_k) = \tilde{L}x_k$. Likewise, equation (95) with $\zeta_{k+1} = \tilde{s}'_{k+1}$, $\eta(x_k) = u_{sv}(x_k)$, $\mathbf{M} = \partial \|\cdot\|_1 \circ \tilde{K}^{\frac{1}{2}}$, and $\varphi(x_k) = \tilde{K}^{-\frac{1}{2}} s_k$ belongs to this class. Equations (124)–(130) also have this form with $\zeta_{k+1} = \tilde{x}_{0,k+1}$, $\eta(x_k) = u_k(x_{0,k}, \nu_k)$, $\mathbf{M}$ as in (111), $\varphi(x_k) = x_{0,k} + h\nu_k$. Equations (148)–(149) are an instance of (182) with $\zeta_{k+1} = s_{0,k+1}$, $\mathbf{M} = \mathbf{M}_{\text{sing}} + \mathbf{M}_{\text{set},h}$, $\varphi(x_k) = -b_k$. The unknown $\zeta_{k+1}$
Maximal monotone operators and their resolvents provide a powerful and unified tool for implicit algorithms computation.

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

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

<figure><img src="image_10.png" /><figcaption>Figure 5: Graphical computation of resolvent and Yosida approximation (lower left and right figures, respectively) for the sgn multifunction in the scalar case. Note that, at each step, a maximal monotone operator is obtained.</figcaption></figure>

for the generalized equation in (182a) corresponds to a “virtual”, “nominal” or “unperturbed” sliding variable, which can be used to determine precisely whether or not a sliding mode is occurring in discrete time. The generalized equation in (182b) has the unknown $\eta(x_k)$ and can be interpreted as the dual, or reciprocal of (182a). From (182a) the selection $\eta(x_k)$ is explicitly given as

$$-\eta(x_k) = \frac{1}{\rho} (I_d - \mathcal{J}_{\rho M}) (\varphi(x_k)) = \mathcal{Y}_{\rho M}(\varphi(x_k)) \quad (184)$$

(see (50) in “Sidebar VII”), whereas (182b) gives

$$-\eta(x_k) = \mathcal{J}_{\rho} M^{-1} \left( \frac{1}{\rho} \varphi(x_k) \right). \quad (185)$$

The expressions (184)–(185) yield the same control law. The main difference relies on the proximal map to be computed. That is, for a certain $M$, it could be easier to compute first the inverse map $M^{-1}$ and then (185), while for other cases (184) might be a simpler choice. In any case, a convex optimization problem has to be solved at each time-step.

For $M = \partial \|\cdot\|_p$, the calculation of the control inputs can be performed according to (39)–(41) in “Sidebar VI”, so that

$$-\eta(x_k) = \mathrm{Proj} \left( \mathcal{B}_{p^*}; \frac{1}{\rho} \varphi(x_k) \right). \quad (186)$$

When $p=1$ then $p^* = \infty$ and, according to (44), the control can be computed component-wise as

$$\eta(x_k)_i = - \mathbf{sgn}(\varphi(x_k)_i) \min \left\{ \frac{1}{\rho} |\varphi(x_k)_i|, 1 \right\}, \quad i = 1, \dots, n,$$

which is a minimum-operator formalism, whereas if $p = p^* = 2$ then, according to (43), it is inferred that

$$\eta(x_k) = \begin{cases} -\frac{1}{\rho}\varphi(x_k) & \text{if } \|\varphi(x_k)\| \le \rho, \\ -\frac{1}{\|\varphi(x_k)\|_2}\varphi(x_k) & \text{otherwise.} \end{cases}$$

For sliding-mode control, it is also common to find maps $f$ given as the composition of norm functions with linear maps. In this case, the projectors change accordingly by considering instead a weighted norm. Indeed, if $M = \partial(\|\cdot\|_p \circ R)$, where $R$ is non-singular, then it follows from [S28, Proposition 23.23] that

$$\begin{align}
\mathcal{J}_{\rho M}(\varphi(x_k)) & = (I_d - R^{\top} \circ (RR^{\top} + (\partial \rho \|\cdot\|_p)^{-1})^{-1} \circ R) (\varphi(x_k)) \nonumber \\
& = (I_d - R^{\top} \circ (RR^{\top} + N_{\rho \mathcal{B}_{p^*}})^{-1} \circ R) (\varphi(x_k)) \nonumber \\
& = \varphi(x_k) - R^{\top} \mathrm{Proj}^{RR^{\top}} (\rho \mathcal{B}_{p^*}; R^{-\top} \varphi(x_k)) \tag{187}
\end{align}$$

where we have made use of (39) to get the second equality and (45c) to get the last equality, both in the sidebar “Sidebar VI”.

In the more general case in which $M = \partial(g \circ R)$ with $g$ a proper, convex, LSC function, we have

$$\mathcal{J}_{\rho} M^{-1} \left( \frac{1}{\rho} \varphi(x_k) \right) = R^{\top} \theta, \quad (188)$$

where

$$\theta = \arg \min_{\xi \in \mathbb{R}^n} \left\{ g^*(\xi) + \frac{\rho}{2} \left( \xi - \frac{1}{\rho} R^{-\top} \varphi(x_k) \right) RR^{\top} \left( \xi - \frac{1}{\rho} R^{-\top} \varphi(x_k) \right) \right\}. \quad (189)$$

The following proposition will be useful for computing proximal maps associated with higher-order sliding-mode controllers and differentiators, as well as their associated splittings.
Maximal monotone operators and their resolvents provide a powerful and unified tool for implicit algorithms computation.

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

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

<figure><img src="image_8.png" /><figcaption>Figure 5: Graphical computation of resolvent and Yosida approximation (lower left and right figures, respectively) for the sgn multifunction in the scalar case. Note that, at each step, a maximal monotone operator is obtained.</figcaption></figure>

for the generalized equation in (182a) corresponds to a “virtual”, “nominal” or “unperturbed” sliding variable, which can be used to determine precisely whether or not a sliding mode is occurring in discrete time. The generalized equation in (182b) has the unknown $\eta(x_k)$ and can be interpreted as the dual, or reciprocal of (182a). From (182a) the selection $\eta(x_k)$ is explicitly given as

$$-\eta(x_k) = \frac{1}{\rho} (I_d - \mathcal{J}_{\rho M}) (\varphi(x_k)) = \mathcal{Y}_{\rho M}(\varphi(x_k)) \quad (184)$$

(see (50) in “Sidebar VII”), whereas (182b) gives

$$-\eta(x_k) = \mathcal{J}_{\rho} M^{-1} \left( \frac{1}{\rho} \varphi(x_k) \right). \quad (185)$$

The expressions (184)–(185) yield the same control law. The main difference relies on the proximal map to be computed. That is, for a certain $M$, it could be easier to compute first the inverse map $M^{-1}$ and then (185), while for other cases (184) might be a simpler choice. In any case, a convex optimization problem has to be solved at each time-step.

For $M = \partial \|\cdot\|_p$, the calculation of the control inputs can be performed according to (39)–(41) in “Sidebar VI”, so that

$$-\eta(x_k) = \mathrm{Proj} \left( \mathcal{B}_{p^*}; \frac{1}{\rho} \varphi(x_k) \right). \quad (186)$$

When $p = 1$ then $p^* = \infty$ and, according to (44), the control can be computed component-wise as

$$\eta(x_k)_i = - \mathbf{sgn}(\varphi(x_k)_i) \min \left\{ \frac{1}{\rho} |\varphi(x_k)_i|, 1 \right\}, \quad i = 1, \dots, n,$$

which is a minimum-operator formalism, whereas if $p = p^* = 2$ then, according to (43), it is inferred that

$$\eta(x_k) = \begin{cases} -\frac{1}{\rho}\varphi(x_k) & \text{if } \|\varphi(x_k)\| \le \rho, \\ -\frac{1}{\|\varphi(x_k)\|_2}\varphi(x_k) & \text{otherwise.} \end{cases}$$

For sliding-mode control, it is also common to find maps $f$ given as the composition of norm functions with linear maps. In this case, the projectors change accordingly by considering instead a weighted norm. Indeed, if $M = \partial(\|\cdot\|_p \circ R)$, where $R$ is non-singular, then it follows from [S28, Proposition 23.23] that

$$\begin{aligned} \mathcal{J}_{\rho M}(\varphi(x_k)) &= (I_d - R^{\top} \circ (RR^{\top} + (\partial \rho \|\cdot\|_p)^{-1})^{-1} \circ R) (\varphi(x_k)) \\ &= (I_d - R^{\top} \circ (RR^{\top} + N_{\rho \mathcal{B}_{p^*}})^{-1} \circ R) (\varphi(x_k)) \\ &= \varphi(x_k) - R^{\top} \mathrm{Proj}^{RR^{\top}} (\rho \mathcal{B}_{p^*}; R^{-\top} \varphi(x_k)) \end{aligned} \quad (187)$$

where we have made use of (39) to get the second equality and (45c) to get the last equality, both in the sidebar “Sidebar VI”.

In the more general case in which $M = \partial(g \circ R)$ with $g$ a proper, convex, LSC function, we have

$$\mathcal{J}_{\rho} M^{-1} \left( \frac{1}{\rho} \varphi(x_k) \right) = R^{\top} \theta, \quad (188)$$

where

$$\theta = \arg \min_{\xi \in \mathbb{R}^n} \left\{ g^*(\xi) + \frac{\rho}{2} \left( \xi - \frac{1}{\rho} R^{-\top} \varphi(x_k) \right) RR^{\top} \left( \xi - \frac{1}{\rho} R^{-\top} \varphi(x_k) \right) \right\}. \quad (189)$$

The following proposition will be useful for computing proximal maps associated with higher-order sliding-mode controllers and differentiators, as well as their associated splittings.
**Proposition 9.** Let $\mathbf{M}: \mathbb{R} \rightrightarrows \mathbb{R}$ be the maximal monotone map

$$ \mathbf{M}(x) = \left( \sum_{i=0}^{N} a_i |x|^{\frac{p_i}{q_i}} \right) \mathbf{sgn}(x) + b_0, \quad (190) $$

where $a_i \ge 0$ and $\frac{p_i}{q_i}$ is a nonnegative rational number for all $i \in \{1, ..., N\}$ and $p_0 = 0$. Then, for $\rho > 0$ the resolvent $\mathcal{J}_{\rho\mathbf{M}}$ has the form

$$ \mathcal{J}_{\rho\mathbf{M}}(x) = \beta_{\mathbf{M}}(x, b_0, \rho)^C \mathbf{sgn}(x - \rho b_0), \quad (191) $$

where $C$ is the least common denominator of the set of fractions $\{\frac{p_i}{q_i} | i = 1, ..., N\}$ and $\beta_{\mathbf{M}}: \mathbb{R} \times \mathbb{R} \times \mathbb{R}_{++} \to \mathbb{R}_{+}$ is the unique nonnegative root of the polynomial

$$ \beta^C + \rho \sum_{i=1}^{N} a_i \beta^{r_i} - \max\{0, |x - \rho b_0| - \rho a_0\} \quad (192) $$

with indeterminate $\beta$ and powers $r_i = C^{\frac{p_i}{q_i}} \in \mathbb{Z}_+$.

*Proof.* From the definition of resolvent in (23), it follows that $y = \mathcal{J}_{\rho\mathbf{M}}(x)$ if, and only if, $x \in y + \rho\mathbf{M}(y)$. That is,

$$ x - \rho b_0 \in \left( |y| + \rho a_0 + \rho \sum_{i=1}^{N} a_i |y|^{\frac{p_i}{q_i}} \right) \mathbf{sgn}(y). \quad (193) $$

Since all the coefficients $a_i$ are nonnegative, it follows from (193) that $y = 0$ if, and only if, $|x - \rho b_0| \le \rho a_0$. Hence, if $|x - \rho b_0| > \rho a_0$, then $y \ne 0$ and from (193) it follows that $\mathbf{sgn}(x - \rho b_0) = \mathbf{sgn}(y)$, so that (193) becomes single-valued as

$$ |y| + \rho \sum_{i=1}^{N} a_i |y|^{\frac{p_i}{q_i}} - (|x - \rho b_0| - \rho a_0) = 0. \quad (194) $$

Let $C > 0$ be the least common denominator of the set $\{\frac{p_i}{q_i} | i = 1, ..., N\}$. The change of variables $|y| = \beta^C$, where $\beta > 0$, transforms (194) into a root-finding problem for the polynomial

$$ \beta^C + \rho \sum_{i=1}^{N} a_i \beta^{r_i} - (|x - \rho b_0| - \rho a_0), \quad (195) $$

where $r_i = C^{\frac{p_i}{q_i}}$. Hence, in order to solve (194) let us look for the positive root of (195). Note that, the non-negativity of all $a_i$, $i \in \{1, ..., N\}$ together with Descartes' rule of signs [125, Theorem 2.23] imply the uniqueness of such positive root of (195) in the cases when $|x - \rho b_0| > \rho a_0$. The polynomial (192) includes both cases above. Indeed, if $|x - \rho b_0| \le \rho a_0$, then the unique non-negative root of (192) is $\beta = 0$, implying $y = 0$. Finally, it follows that

$$ \mathcal{J}_{\rho\mathbf{M}}(x) = y = |y| \mathbf{sgn}(y) = \beta_{\mathbf{M}}(x, b_0, \rho)^C \mathbf{sgn}(x - \rho b_0). \quad (196) $$

The proof is complete. □

Instances of Proposition 9 have been independently developed in [84] and [83] in the context of higher-order differentiators, since (191) is the unique solution of the generalized equation

$$ 0 \in y + \sum_{i=0}^{N} a_i |y|^{\frac{p_i}{q_i}} \mathbf{sgn}(y) + b_0 - x. \quad (197) $$

## 6.2 Computation of proximal maps via splittings

For complex strategies such as higher-order sliding-mode control and differentiation, the associated monotone operator is composed as the addition of several simpler operators. In this situation, the splitting algorithms in “Sidebar X” can be used to compute the associated proximal map in an iterative way.

**Example 1.** The splitting (89) can be used to compute the URED proximal map in (148), (149) as follows. Consider the optimization problem (84) with $A = 1$ and

$$ f(\xi) = \frac{2}{3}h\gamma_0 L^{\frac{1}{2}}|\xi|^{\frac{3}{2}} + \frac{1}{2}h^2\gamma_1 L|\xi|, \quad (198a) $$

$$ g(\xi) = \frac{1}{2}h^2\mu^2\gamma_1 L|\xi|^3, \quad (198b) $$

$$ h(\xi) = \frac{2}{5}h\gamma_0 L^{\frac{1}{2}}\mu|\xi|^{\frac{5}{2}} + h^2\gamma_1 L\mu\xi^2 + \frac{1}{2}(\xi + b_k)^2. \quad (198c) $$

With the given parameters, the solution of (84) indeed gives the proximal map $\mathcal{J}_{h\mathbf{M}_h}(-b_k)$, where $\mathbf{M}_h$ is as in (149). Also, note that $\nabla h(\cdot)$ is Lipschitz continuous in any bounded set of $\mathbb{R}$, since

$$ \frac{\partial^2}{\partial \xi^2} h(\xi) = \frac{3}{2}h\gamma_0 L^{\frac{1}{2}}\mu|\xi|^{\frac{1}{2}} + (1 + 2h^2\gamma_1 L\mu). \quad (199) $$

Notice that, with such decomposition, each map appearing in (89) is proximable. Indeed, it follows from Proposition 9 that

$$ \text{Prox}_{\alpha_1 f}(x) = \mathcal{J}_{\alpha_1 \partial f}(x) = \beta_{\partial f}(x, \alpha_1)^2 \mathbf{sgn}(x), \quad (200) $$

where

$$ \beta_{\partial f}(x, \alpha_1) = \frac{-\alpha_1 h \gamma_0 L^{\frac{1}{2}} + \sqrt{\left(\alpha_1 h \gamma_0 L^{\frac{1}{2}}\right)^2 + 4\rho(x)}}{2}, \quad (201) $$

where $\rho(x) = \max\{0, |x| - \frac{\alpha_1 h^2 \gamma_1 L}{2}\}$. Now, in order to compute the proximal map associated with $g^*(\cdot)$ in (89), we make use of Moreau's decomposition (see, e.g., [S28, Theorem 14.3] in "Sidebar V"),

$$ x = \text{Prox}_{\alpha_2 g^*}(x) + \alpha_2 \text{Prox}_{\frac{1}{\alpha_2} g}\left(\frac{x}{\alpha_2}\right). \quad (202) $$

Thus, let us compute $\text{Prox}_{\frac{1}{\alpha_2} g}\left(\frac{x}{\alpha_2}\right)$ using once again Proposition 9, so that

$$ \text{Prox}_{\frac{1}{\alpha_2} g}\left(\frac{x}{\alpha_2}\right) = \beta_{\partial g}\left(\frac{x}{\alpha_2}, \frac{1}{\alpha_2}\right) \mathbf{sgn}(x), \quad (203) $$
**Proposition 9.** Let $\mathbf{M}: \mathbb{R} \rightrightarrows \mathbb{R}$ be the maximal monotone map

$$ \mathbf{M}(x) = \left( \sum_{i=0}^{N} a_i |x|^{\frac{p_i}{q_i}} \right) \mathbf{sgn}(x) + b_0, \quad (190) $$

where $a_i \ge 0$ and $\frac{p_i}{q_i}$ is a nonnegative rational number for all $i \in \{1, ..., N\}$ and $p_0 = 0$. Then, for $\rho > 0$ the resolvent $\mathcal{J}_{\rho \mathbf{M}}$ has the form

$$ \mathcal{J}_{\rho \mathbf{M}}(x) = \beta_{\mathbf{M}}(x, b_0, \rho)^C \mathbf{sgn}(x - \rho b_0), \quad (191) $$

where $C$ is the least common denominator of the set of fractions $\{\frac{p_i}{q_i} | i = 1, ..., N\}$ and $\beta_{\mathbf{M}}: \mathbb{R} \times \mathbb{R} \times \mathbb{R}_{++} \to \mathbb{R}_{+}$ is the unique nonnegative root of the polynomial

$$ \beta^C + \rho \sum_{i=1}^{N} a_i \beta^{r_i} - \max\{0, |x - \rho b_0| - \rho a_0\} \quad (192) $$

with indeterminate $\beta$ and powers $r_i = C^{\frac{p_i}{q_i}} \in \mathbb{Z}_{+}$.

*Proof.* From the definition of resolvent in (23), it follows that $y = \mathcal{J}_{\rho \mathbf{M}}(x)$ if, and only if, $x \in y + \rho \mathbf{M}(y)$. That is,

$$ x - \rho b_0 \in \left( |y| + \rho a_0 + \rho \sum_{i=1}^{N} a_i |y|^{\frac{p_i}{q_i}} \right) \mathbf{sgn}(y). \quad (193) $$

Since all the coefficients $a_i$ are nonnegative, it follows from (193) that $y = 0$ if, and only if, $|x - \rho b_0| \le \rho a_0$. Hence, if $|x - \rho b_0| > \rho a_0$, then $y \ne 0$ and from (193) it follows that $\mathbf{sgn}(x - \rho b_0) = \mathbf{sgn}(y)$, so that (193) becomes single-valued as

$$ |y| + \rho \sum_{i=1}^{N} a_i |y|^{\frac{p_i}{q_i}} - (|x - \rho b_0| - \rho a_0) = 0. \quad (194) $$

Let $C > 0$ be the least common denominator of the set $\{\frac{p_i}{q_i} | i = 1, ..., N\}$. The change of variables $|y| = \beta^C$, where $\beta > 0$, transforms (194) into a root-finding problem for the polynomial

$$ \beta^C + \rho \sum_{i=1}^{N} a_i \beta^{r_i} - (|x - \rho b_0| - \rho a_0), \quad (195) $$

where $r_i = C^{\frac{p_i}{q_i}}$. Hence, in order to solve (194) let us look for the positive root of (195). Note that, the non-negativity of all $a_i$, $i \in \{1, ..., N\}$ together with Descartes' rule of signs [125, Theorem 2.23] imply the uniqueness of such positive root of (195) in the cases when $|x - \rho b_0| > \rho a_0$. The polynomial (192) includes both cases above. Indeed, if $|x - \rho b_0| \le \rho a_0$, then the unique non-negative root of (192) is $\beta = 0$, implying $y = 0$. Finally, it follows that

$$ \mathcal{J}_{\rho \mathbf{M}}(x) = y = |y| \mathbf{sgn}(y) = \beta_{\mathbf{M}}(x, b_0, \rho)^C \mathbf{sgn}(x - \rho b_0). \quad (196) $$

The proof is complete. $\square$

Instances of Proposition 9 have been independently developed in [84] and [83] in the context of higher-order differentiators, since (191) is the unique solution of the generalized equation

$$ 0 \in y + \sum_{i=0}^{N} a_i |y|^{\frac{p_i}{q_i}} \mathbf{sgn}(y) + b_0 - x. \quad (197) $$

## 6.2 Computation of proximal maps via splittings

For complex strategies such as higher-order sliding-mode control and differentiation, the associated monotone operator is composed as the addition of several simpler operators. In this situation, the splitting algorithms in “Sidebar X” can be used to compute the associated proximal map in an iterative way.

**Example 1.** The splitting (89) can be used to compute the URED proximal map in (148), (149) as follows. Consider the optimization problem (84) with $A = 1$ and

$$ f(\xi) = \frac{2}{3}h\gamma_0 L^{\frac{1}{2}}|\xi|^{\frac{3}{2}} + \frac{1}{2}h^2\gamma_1 L|\xi|, \quad (198a) $$

$$ g(\xi) = \frac{1}{2}h^2\mu^2\gamma_1 L|\xi|^3, \quad (198b) $$

$$ h(\xi) = \frac{2}{5}h\gamma_0 L^{\frac{1}{2}}\mu|\xi|^{\frac{5}{2}} + h^2\gamma_1 L\mu\xi^2 + \frac{1}{2}(\xi + b_k)^2. \quad (198c) $$

With the given parameters, the solution of (84) indeed gives the proximal map $\mathcal{J}_{h \mathbf{M}_h}(-b_k)$, where $\mathbf{M}_h$ is as in (149). Also, note that $\nabla h(\cdot)$ is Lipschitz continuous in any bounded set of $\mathbb{R}$, since

$$ \frac{\partial^2}{\partial \xi^2} h(\xi) = \frac{3}{2}h\gamma_0 L^{\frac{1}{2}}\mu|\xi|^{\frac{1}{2}} + (1 + 2h^2\gamma_1 L\mu). \quad (199) $$

Notice that, with such decomposition, each map appearing in (89) is proximable. Indeed, it follows from Proposition 9 that

$$ \mathrm{Prox}_{\alpha_1 f}(x) = \mathcal{J}_{\alpha_1 \partial f}(x) = \beta_{\partial f}(x, \alpha_1)^2 \mathbf{sgn}(x), \quad (200) $$

where

$$ \beta_{\partial f}(x, \alpha_1) = \frac{-\alpha_1 h \gamma_0 L^{\frac{1}{2}} + \sqrt{(\alpha_1 h \gamma_0 L^{\frac{1}{2}})^2 + 4\rho(x)}}{2}, \quad (201) $$

where $\rho(x) = \max\{0, |x| - \frac{\alpha_1 h^2 \gamma_1 L}{2}\}$. Now, in order to compute the proximal map associated with $g^*(\cdot)$ in (89), we make use of Moreau's decomposition (see, e.g., [S28, Theorem 14.3] in "Sidebar V"),

$$ x = \mathrm{Prox}_{\alpha_2 g^*}(x) + \alpha_2 \mathrm{Prox}_{\frac{1}{\alpha_2} g}\left(\frac{x}{\alpha_2}\right). \quad (202) $$

Thus, let us compute $\mathrm{Prox}_{\frac{1}{\alpha_2} g}\left(\frac{x}{\alpha_2}\right)$ using once again Proposition 9, so that

$$ \mathrm{Prox}_{\frac{1}{\alpha_2} g}\left(\frac{x}{\alpha_2}\right) = \beta_{\partial g}\left(\frac{x}{\alpha_2}, \frac{1}{\alpha_2}\right) \mathbf{sgn}(x), \quad (203) $$
where

$$ \beta_{\partial g}(x, \rho) = \frac{-1 + \sqrt{1 + 6h^2\mu^2\gamma_1 L\rho|x|}}{3h^2\mu^2\gamma_1 L\rho}. \quad (204) $$

Finally, it follows from (202) and (203) that

$$ \begin{aligned} \operatorname{Prox}_{\alpha_2 g^*}(x) &= x - \alpha_2 \operatorname{Prox}_{\frac{1}{\alpha_2} g} \left( \frac{x}{\alpha_2} \right) \\ &= x - \alpha_2 \beta_{\partial g}(x, \rho) \operatorname{sgn}(x). \end{aligned} \quad (205) $$

Therefore, by using the expressions (200) and (205) in the iteration (89) and by selecting the parameters $\alpha_1$ and $\alpha_2$ sufficiently small such that (90) holds, we have that $\xi_{j+1} \to \mathcal{J}_h M_h(-b_k)$ as $j \uparrow +\infty$.

It is worth mentioning that, for problems with more than two maximal monotone terms, parallel algorithms as those described in [126, S35] can be used.

**Example 2.** The Douglas-Rachford splitting can also be used to compute the proximal map in (148), (149). This time a splitting of the following form is considered:

$$ f(\xi) = \frac{2}{5}h\gamma_0 L^{\frac{1}{2}}\mu|\xi|^{\frac{5}{2}} + \frac{2}{3}h\gamma_0 L^{\frac{1}{2}}|\xi|^{\frac{3}{2}} + \frac{1}{2}h^2\gamma_1 L|\xi|, \quad (206) $$

$$ g(\xi) = \frac{1}{2}h^2\mu^2\gamma_1 L|\xi|^3 + h^2\gamma_1 L\mu\xi^2 + \frac{1}{2}(\xi + b_k)^2. \quad (207) $$

From Proposition 9, the expressions

$$ \operatorname{Prox}_{\alpha f}(x) = \beta_{\partial f}(x, \alpha)^2 \operatorname{sgn}(x), \quad (208) $$

$$ \operatorname{Prox}_{\alpha g}(x) = \beta_{\partial g}(x, b_k, \alpha) \operatorname{sgn}(x - \alpha b_k), \quad (209) $$

are obtained, where $\beta_{\partial f}(x, \alpha)$ is the unique positive root of the cubic polynomial

$$ \alpha h \lambda_0 L^{\frac{1}{2}} \mu \beta^3 + \beta^2 + \alpha \lambda_0 L^{\frac{1}{2}} \beta - \max \left\{ 0, |x| - \frac{\alpha h^2 \lambda_1 L}{2} \right\}, \quad (210) $$

that is,

$$ \beta_{\partial g}(x, b_k, \alpha) = \frac{-2(1 + h^2 \mu \lambda_1 L)}{3h^2 \mu^2 \lambda_1 L} + \frac{\sqrt{4(1 + h^2 \mu \lambda_1 L)^2 + 6h^2 \mu^2 \lambda_1 L |x - b_k|}}{3h^2 \mu^2 \lambda_1 L}. \quad (211) $$

Hence, by taking the Douglas-Rachford iteration (94) with $\operatorname{Prox}_{\alpha f}$ and $\operatorname{Prox}_{\alpha g}$ as in (208)-(211), we have that $\zeta_j$ converges to the resolvent (148)-(149) as $j \uparrow +\infty$.

## 6.3 Relaxation

One of the main objectives of sliding-mode control concerns the finite-time stability of the sliding variable, whose discrete-time dynamics are of the form

$$ s_{k+1} = s_k + \rho u_{sv}(x_k) + \delta_k. \quad (212) $$

Notice that, as the disturbance $\delta_k$ is assumed unknown, the finite-time stability of $s_k$ cannot be achieved, so that a

<figure><img src="image_11.png" /><figcaption>Figure 6: Closed-loop system. The state of the continuous-time dynamics is sampled every $h$ seconds. Such information is used for computing the control input $u_k$ which passes through a classical zero-order-hold (ZOH) mechanism to obtain the feedback input signal $u(t) = u_k$ for $kh < t \le (k+1)h$.</figcaption></figure>

new variable $\tilde{s}_k$ is introduced. The form of such “virtual” sliding variable is largely open. One option consists in considering a copy of the nominal behavior as in (182). One way of accelerating the convergence of $\tilde{s}_k$ towards the origin consists in setting the control law as

$$ -u_{sv}(x_k) \in \partial f(\tilde{s}_{k+1}) \quad (213) $$

$$ \tilde{s}_{k+1} = s_k + \mu u_{sv}(x_k), \quad (214) $$

where $\mu > 0$, so that the closed-loop is given as

$$ \tilde{s}_{k+1} = \operatorname{Prox}_{\mu f}(s_k) \quad (215a) $$

$$ s_{k+1} = s_k + \frac{\rho}{\mu}(\tilde{s}_{k+1} - s_k) + \delta_k. \quad (215b) $$

The iteration (215) is known as a relaxed iteration in the literature of proximal-point algorithms [127]. In the cases when $\delta_k \equiv 0$, it has been shown that (215) converges towards a minimum of $f(\cdot)$ [S28, Corollary 5.16]. Moreover, for certain values of $\frac{\rho}{\mu} \in (0, 2)$, the algorithm (215) has a faster convergence with respect to the case $\rho = \mu$, see, e.g., [127].

## 6.4 Numerical Experiments

Taking advantage of the fact that the super-twisting controller can be calculated in the closed form (128), it is used as a benchmark to test algorithms which approximate implicit or semi-implicit discrete-time inputs. It appears that there are several ways to implement a given controller or differentiator, and it is certainly crucial to launch a research effort to tackle such issues. Let us consider a perturbed integrator in feedback with a discrete-time controller as shown in Figure 6. The target amounts to driving the state $x_0(t)$ towards the origin in the presence of the persistent disturbance $\delta(t) = \cos(\pi t) \sin(5\sqrt{\frac{t}{3}})$. To this end, let us consider the super-twisting algorithm, taking in to account the two main approaches for its discrete-time implementation, i.e., emulation of the continuous-time controller and discrete-time design. For the emulation part, i) the explicit (forward) Euler discretization,
where

$$ \beta_{\partial g}(x, \rho) = \frac{-1 + \sqrt{1 + 6h^2\mu^2\gamma_1 L\rho|x|}}{3h^2\mu^2\gamma_1 L\rho}. \quad (204) $$

Finally, it follows from (202) and (203) that

$$ \begin{aligned} \operatorname{Prox}_{\alpha_2 g^*}(x) &= x - \alpha_2 \operatorname{Prox}_{\frac{1}{\alpha_2} g} \left( \frac{x}{\alpha_2} \right) \\ &= x - \alpha_2 \beta_{\partial g}(x, \rho) \operatorname{sgn}(x). \end{aligned} \quad (205) $$

Therefore, by using the expressions (200) and (205) in the iteration (89) and by selecting the parameters $\alpha_1$ and $\alpha_2$ sufficiently small such that (90) holds, we have that $\xi_{j+1} \to \mathcal{J}_h M_h(-b_k)$ as $j \uparrow +\infty$.

It is worth mentioning that, for problems with more than two maximal monotone terms, parallel algorithms as those described in [126, S35] can be used.

**Example 2.** The Douglas-Rachford splitting can also be used to compute the proximal map in (148), (149). This time a splitting of the following form is considered:

$$ f(\xi) = \frac{2}{5}h\gamma_0 L^{\frac{1}{2}}\mu|\xi|^{\frac{5}{2}} + \frac{2}{3}h\gamma_0 L^{\frac{1}{2}}|\xi|^{\frac{3}{2}} + \frac{1}{2}h^2\gamma_1 L|\xi|, \quad (206) $$

$$ g(\xi) = \frac{1}{2}h^2\mu^2\gamma_1 L|\xi|^3 + h^2\gamma_1 L\mu\xi^2 + \frac{1}{2}(\xi + b_k)^2. \quad (207) $$

From Proposition 9, the expressions

$$ \operatorname{Prox}_{\alpha f}(x) = \beta_{\partial f}(x, \alpha)^2 \operatorname{sgn}(x), \quad (208) $$

$$ \operatorname{Prox}_{\alpha g}(x) = \beta_{\partial g}(x, b_k, \alpha) \operatorname{sgn}(x - \alpha b_k), \quad (209) $$

are obtained, where $\beta_{\partial f}(x, \alpha)$ is the unique positive root of the cubic polynomial

$$ \alpha h \lambda_0 L^{\frac{1}{2}} \mu \beta^3 + \beta^2 + \alpha \lambda_0 L^{\frac{1}{2}} \beta - \max \left\{ 0, |x| - \frac{\alpha h^2 \lambda_1 L}{2} \right\}, \quad (210) $$

that is,

$$ \beta_{\partial g}(x, b_k, \alpha) = \frac{-2(1 + h^2 \mu \lambda_1 L)}{3h^2 \mu^2 \lambda_1 L} + \frac{\sqrt{4(1 + h^2 \mu \lambda_1 L)^2 + 6h^2 \mu^2 \lambda_1 L |x - b_k|}}{3h^2 \mu^2 \lambda_1 L}. \quad (211) $$

Hence, by taking the Douglas-Rachford iteration (94) with $\operatorname{Prox}_{\alpha f}$ and $\operatorname{Prox}_{\alpha g}$ as in (208)-(211), we have that $\zeta_j$ converges to the resolvent (148)-(149) as $j \uparrow +\infty$.

## 6.3 Relaxation

One of the main objectives of sliding-mode control concerns the finite-time stability of the sliding variable, whose discrete-time dynamics are of the form

$$ s_{k+1} = s_k + \rho u_{sv}(x_k) + \delta_k. \quad (212) $$

Notice that, as the disturbance $\delta_k$ is assumed unknown, the finite-time stability of $s_k$ cannot be achieved, so that a

<figure><img src="image_9.png" /><figcaption>Figure 6: Closed-loop system. The state of the continuous-time dynamics is sampled every $h$ seconds. Such information is used for computing the control input $u_k$ which passes through a classical zero-order-hold (ZOH) mechanism to obtain the feedback input signal $u(t) = u_k$ for $kh < t \le (k+1)h$.</figcaption></figure>

new variable $\tilde{s}_k$ is introduced. The form of such “virtual” sliding variable is largely open. One option consists in considering a copy of the nominal behavior as in (182). One way of accelerating the convergence of $\tilde{s}_k$ towards the origin consists in setting the control law as

$$ -u_{sv}(x_k) \in \partial f(\tilde{s}_{k+1}) \quad (213) $$

$$ \tilde{s}_{k+1} = s_k + \mu u_{sv}(x_k), \quad (214) $$

where $\mu > 0$, so that the closed-loop is given as

$$ \tilde{s}_{k+1} = \operatorname{Prox}_{\mu f}(s_k) \quad (215a) $$

$$ s_{k+1} = s_k + \frac{\rho}{\mu}(\tilde{s}_{k+1} - s_k) + \delta_k. \quad (215b) $$

The iteration (215) is known as a relaxed iteration in the literature of proximal-point algorithms [127]. In the cases when $\delta_k \equiv 0$, it has been shown that (215) converges towards a minimum of $f(\cdot)$ [S28, Corollary 5.16]. Moreover, for certain values of $\frac{\rho}{\mu} \in (0, 2)$, the algorithm (215) has a faster convergence with respect to the case $\rho = \mu$, see, e.g., [127].

## 6.4 Numerical Experiments

Taking advantage of the fact that the super-twisting controller can be calculated in the closed form (128), it is used as a benchmark to test algorithms which approximate implicit or semi-implicit discrete-time inputs. It appears that there are several ways to implement a given controller or differentiator, and it is certainly crucial to launch a research effort to tackle such issues. Let us consider a perturbed integrator in feedback with a discrete-time controller as shown in Figure 6. The target amounts to driving the state $x_0(t)$ towards the origin in the presence of the persistent disturbance $\delta(t) = \cos(\pi t) \sin(5\sqrt{\frac{t}{3}})$. To this end, let us consider the super-twisting algorithm, taking in to account the two main approaches for its discrete-time implementation, i.e., emulation of the continuous-time controller and discrete-time design. For the emulation part, i) the explicit (forward) Euler discretization,
ii) the implicit (backward) Euler discretization in [34] and
iii) the semi-implicit discretization in [35], are considered.
In each case the control law takes the form

$$
\begin{align}
u_k &= -\gamma_0 L^{\frac{1}{2}} |x_{0,k}|^{\frac{1}{2}} \operatorname{sgn}(x_{0,k}) + \nu_k \tag{216} \\
\nu_{k+1} &= \nu_k - h\gamma_1 L \operatorname{sgn}(x_{0,k}), \nonumber
\end{align}
$$

and

$$
\begin{align*}
u_k &= -\gamma_0 L^{\frac{1}{2}} \beta_k \operatorname{sgn}(x_{0,k} + h\nu_k) + \nu_k \\
    &\quad - h\gamma_1 L \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h^2\gamma_1 L}\right)
\end{align*}
$$

$$
\begin{align}
\beta_k &= \frac{-h\gamma_0 L^{\frac{1}{2}}}{2} + \sqrt{\frac{h^2\gamma_0^2 L}{4} + \max\{0, |x_{0,k} + h\nu_k| - h^2\gamma_1 L\}} \tag{217} \\
\nu_{k+1} &= \nu_k - h\gamma_1 L \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h^2\gamma_1 L}\right), \nonumber
\end{align}
$$

and

$$
\begin{align}
u_k &= \nu_k - \mu_k \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h\mu_k}\right) \notag \\
\nu_{k+1} &= \nu_k - h\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h\mu_k}\right) \tag{218} \\
\mu_k &= \gamma_0 L^{\frac{1}{2}} |x_{0,k}|^{\frac{1}{2}} + h\gamma_1 L. \notag
\end{align}
$$

For the second approach, which relies on the direct design in discrete-time, the following controller, which is proposed in [86], is considered:

$$
\begin{align}
u_k &= -\gamma_0 \beta_k \operatorname{sgn}(x_{1,k}) - 2h\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{x_{1,k}}{h^2\gamma_1}\right) + \nu_k \nonumber \\
\beta_k &= -\frac{h\gamma_0}{2} + \frac{1}{2}\sqrt{(h\gamma_0)^2 + 4\max\{0, |x_{1,k}| - h^2\gamma_1\}} \nonumber \\
\nu_{k+1} &= \nu_k - h\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{x_{1,k}}{h^2\gamma_1}\right) \tag{219}
\end{align}
$$

The controller (219) is in explicit form, though strongly
motivated from the implicit controller obtained from the
backward-Euler discretization (217). In addition, the
splitting schemes (89) and (94) are considered, applied
to the implicit super-twisting algorithm in (110). That
is, for solving the minimization problem

$$
\min_{\xi \in \mathbb{R}} \frac{2}{3} h \gamma_0 L^{\frac{1}{2}} |\xi|^{\frac{3}{2}} + h^2 \gamma_1 L |\xi| + \frac{1}{2} (\xi - (x_{0,k} + h \nu_k))^2 . \quad (220)
$$

As a last scheme to compare, the problem (220) is solved numerically at each time step by using the interior point solver *Clarabel* [128]. The listings 1 and 2 show the code implementation in Python 3 used for the Douglas-Rachford and Condat-Vũ splittings, respectively.

For the simulations, the initial condition $x_0(0) = 5$ is
taken and several discretization steps $100\mu s \le h \le 1s$ are
considered. For each simulation a final time of $t_f = 20s$
is chosen. Figure 7 displays the $L_2$-norm, $L_\infty$-norm, and
the step-to-step variation SVar of the state and control

signals for different sampling times. The latter is defined
for a finite sequence {$\xi_k$} of N samples as

$$
SVar(\{\xi_k\}) = \sum_{k=0}^{N-2} |\xi_{k+1} - \xi_k|. \quad (221)
$$

For large sampling times, 0.5s < h < 1s, all the con-
trollers lead to similar values on the L₂-norm of the state
x₀. Notice that, in terms of the L₂-norm of the control in-
put and the step-to-step variations, the controllers (219)
and (216) are less performing. Indeed, for the explicit
controller (216), SVar[0,20]({uₖ}) ≈ 500 units, regardless
of the sampling time. For smaller sampling times, i.e.,
2ms < h < 20ms, all controllers perform in a similar
way in terms of the L₂-norm of the state and the input,
even though the step-to-step variation of the state with
the controller (216) is also higher that in the other cases.
Finally, when h < 2ms, all the controllers show the same
performance, except for the numerical approach, which
shows a sudden increment in the step-to-step variation
of the input u (that at some point even surpasses that
of the explicit controller (216)). Such increment in the
step-to-step variation, whilst maintaining the L₂-norm
unchanged, is a sign of the presence of chattering in the
input.

<caption>Listing 1: Douglas-Rachford splitting for the super-twisting controller (217)</caption>

import numpy as np

def bracket(x, a):
    return np.sign(x) * np.abs(x) ** a

def J_A0(alpha, x, h_c, g0):
    beta_A = (-alpha * h_c * g0 + np.sqrt((alpha * h_c * g0) ** 2 + 4.0 * np.abs(x))) / 2.0
    return (beta_A ** 2) * np.sign(x)

def J_A1(x, h_c, g1):
    return 0 if np.abs(x) <= (h_c ** 2) * g1 else v - (h_c ** 2) * g1 * np.sign(x)

def sta_dgl_rfd(h_c, g0, g1, x_0_k, nu, x_0_tilde, w,
    N):
    r = x_0_k + h_c * nu[-1]
    z = w[-1]
    for i in range(N):
        R_A1 = 2.0 * J_A1(z, h_c, g1) - z
        z = 0.5 * (z - R_A1) + J_A0(0.5, 0.5 * (r + R_A1), h_c, g0)
        x_0_next_tilde = J_A1(z, h_c, g1)
        u_k = -(x_0_k - x_0_next_tilde) / h_c
        nu_next = u_k + g0 * bracket(x_0_next_tilde, 0.5)
        x_0_tilde.append(x_0_next_tilde)
        nu.append(nu_next)
        w.append(z)
    return u_k

<caption>Listing 2: Condat-Vũ splitting for the super-twisting con-
troller (217)</caption>
ii) the implicit (backward) Euler discretization in [34] and
iii) the semi-implicit discretization in [35], are considered.
In each case the control law takes the form

$$
\begin{align}
u_k &= -\gamma_0 L^{\frac{1}{2}} |x_{0,k}|^{\frac{1}{2}} \operatorname{sgn}(x_{0,k}) + \nu_k \tag{216} \\
\nu_{k+1} &= \nu_k - h\gamma_1 L \operatorname{sgn}(x_{0,k}), \nonumber
\end{align}
$$

and

$$
\begin{align*}
u_k &= -\gamma_0 L^{\frac{1}{2}} \beta_k \operatorname{sgn}(x_{0,k} + h\nu_k) + \nu_k \\
    &\quad - h\gamma_1 L \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h^2\gamma_1 L}\right)
\end{align*}
$$

$$
\begin{align*}
\beta_k &= \frac{-h\gamma_0 L^{\frac{1}{2}}}{2} + \sqrt{\frac{h^2\gamma_0^2 L}{4} + \max\{0, |x_{0,k} + h\nu_k| - h^2\gamma_1 L\}} \quad \text{way in terms of the } L_2\text{-norm of the state and the input,} \\
\nu_{k+1} &= \nu_k - h\gamma_1 L \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h^2\gamma_1 L}\right), \tag{217}
\end{align*}
$$

and

$$
\begin{align}
u_k &= \nu_k - \mu_k \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h\mu_k}\right) \notag \\
\nu_{k+1} &= \nu_k - h\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{x_{0,k} + h\nu_k}{h\mu_k}\right) \tag{218} \\
\mu_k &= \gamma_0 L^{\frac{1}{2}} |x_{0,k}|^{\frac{1}{2}} + h\gamma_1 L. \notag
\end{align}
$$

For the second approach, which relies on the direct design in discrete-time, the following controller, which is proposed in [86], is considered:

$$
\begin{align}
u_k &= -\gamma_0 \beta_k \operatorname{sgn}(x_{1,k}) - 2h\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{x_{1,k}}{h^2\gamma_1}\right) + \nu_k \nonumber \\
\beta_k &= -\frac{h\gamma_0}{2} + \frac{1}{2}\sqrt{(h\gamma_0)^2 + 4\max\{0, |x_{1,k}| - h^2\gamma_1\}} \nonumber \\
\nu_{k+1} &= \nu_k - h\gamma_1 \operatorname{Proj}\left([-1, 1]; \frac{x_{1,k}}{h^2\gamma_1}\right) \tag{219}
\end{align}
$$

The controller (219) is in explicit form, though strongly
motivated from the implicit controller obtained from the
backward-Euler discretization (217). In addition, the
splitting schemes (89) and (94) are considered, applied
to the implicit super-twisting algorithm in (110). That
is, for solving the minimization problem

$$
\min_{\xi \in \mathbb{R}} \frac{2}{3} h \gamma_0 L^{\frac{1}{2}} |\xi|^{\frac{3}{2}} + h^2 \gamma_1 L |\xi| + \frac{1}{2} (\xi - (x_{0,k} + h \nu_k))^2 . \quad (220)
$$

As a last scheme to compare, the problem (220) is solved numerically at each time step by using the interior point solver *Clarabel* [128]. The listings 1 and 2 show the code implementation in Python 3 used for the Douglas-Rachford and Condat-Vũ splittings, respectively.

For the simulations, the initial condition $x_0(0) = 5$ is
taken and several discretization steps $100\mu s \le h \le 1s$ are
considered. For each simulation a final time of $t_f = 20s$
is chosen. Figure 7 displays the $L_2$-norm, $L_\infty$-norm, and
the step-to-step variation SVar of the state and control

signals for different sampling times. The latter is defined
for a finite sequence {$\xi_k$} of N samples as

$$
SVar(\{\xi_k\}) = \sum_{k=0}^{N-2} |\xi_{k+1} - \xi_k|. \quad (221)
$$

For large sampling times, 0.5s < h < 1s, all the con-
trollers lead to similar values on the L₂-norm of the state
x₀. Notice that, in terms of the L₂-norm of the control in-
put and the step-to-step variations, the controllers (219)
and (216) are less performing. Indeed, for the explicit
controller (216), SVar[0,20]({uₖ}) ≈ 500 units, regardless
of the sampling time. For smaller sampling times, i.e.,
2ms < h < 20ms, all controllers perform in a similar
way in terms of the L₂-norm of the state and the input,
even though the step-to-step variation of the state with
the controller (216) is also higher that in the other cases.
Finally, when h < 2ms, all the controllers show the same
performance, except for the numerical approach, which
shows a sudden increment in the step-to-step variation
of the input u (that at some point even surpasses that
of the explicit controller (216)). Such increment in the
step-to-step variation, whilst maintaining the L₂-norm
unchanged, is a sign of the presence of chattering in the
input.

<caption>Listing 1: Douglas-Rachford splitting for the super-twisting controller (217)</caption>

import numpy as np

def bracket(x, a):
    return np.sign(x) * np.abs(x) ** a

def J_A0(alpha, x, h_c, g0):
    beta_A = (-alpha * h_c * g0 +
               np.sqrt((alpha * h_c * g0) ** 2 +
                         4.0 * np.abs(x))) / 2.0
    return (beta_A ** 2) * np.sign(x)

def J_A1(x, h_c, g1):
    return 0 if np.abs(x) <= (h_c ** 2) * g1 else
           v - (h_c ** 2) * g1 * np.sign(x)

def sta_dgl_rfd(h_c, g0, g1, x_0_k, nu, x_0_tilde, w,
                N):
    r = x_0_k + h_c * nu [-1]
    z = w [-1]
    for i in range(N):
        R_A1 = 2.0 * J_A1(z, h_c, g1) - z
        z = 0.5 * (z - R_A1) + J_A0(0.5, 0.5 * (r +
                                                R_A1), h_c, g0)
        x_0_next_tilde = J_A1(z, h_c, g1)
        u_k = -(x_0_k - x_0_next_tilde) / h_c
        nu_next = u_k + g0 * bracket(x_0_next_tilde, 0.5)
        x_0_tilde.append(x_0_next_tilde)
        nu.append(nu_next)
        w.append(z)
    return u_k

<caption>Listing 2: Condat-Vũ splitting for the super-twisting con-
troller (217)</caption>
<figure><img src="image_12.png" /><figcaption>Figure 7: Norms and total variations of state and control signals for several discrete-time super-twisting controllers. In terms of the $L_2$-norm of the state $x_0$ and the control input $u$, the implicit controller (217) (and its associated splittings); the semi-implicit controller (218); and the numerical solver for (220), have a similar behavior for all sampling times. In terms of the step-to-step variation (221) of the control input sequence, the explicit controller shows the larger values. It is noteworthy that for sampling times $h < 2\text{ms}$, the numerical solver presents an increment in the variation of the controller $u$, indicating the presence of numerical chattering. In terms of the $\|x_0\|_{L_\infty(15,20)}$, the controller (219) shows the best performance for $h \le 200\text{ms}$.</figcaption></figure>

import numpy as np

def bracket(x, a):
    return np.sign(x) * np.abs(x) ** a

def J_B0(x, h_c, g0):
    beta = (-h_c*g0 + np/sqrt((h_c*g0)**2 +
    4.0*np.abs(x)))/2.0
    return (beta)**2 * np sign(x)

def J_B1(x, a):
    return x if np abs(x) <= a else a * np sign(x)

def sta_cdt_vu(h_c, g0, g1, x_0_k, nu, x_0_tilde, mu,
               N, a1, a2):
    z = x_0_tilde[-1]
    m = mu[-1]
    for i in range(N):
        z_next = J_B0((1 - a1)*z + a1*(x_0_k +
        h_c*nu[-1] - m, h_c, g0)
        m = J_B1(m + a2*(2.0*z_next - z), g1+h_c**2)
        z = z_next
    x_0_next_tilde = z
    u_k = -(x_0_k - x_0_next_tilde)/h_c
    nu_next = u_k + g0*bracket(x_0_next_tilde, 0.5)
    x_0_tilde.append(x_0_next_tilde)
    nu.append(nu_next)
    mu.append(m)
    return u_k

# 7 EXPERIMENTAL RESULTS AND VALIDATIONS

## 7.1 A Quick Review

Implicit and semi-implicit control and differentiation algorithms have been experimentally tested, mostly on laboratory setups. It is noteworthy that the fundamental properties like insensitivity with respect to the gain during sliding phases (stemming from (77), (98), (137)), moderate deterioration of the performance with increasing sampling times, significant chattering alleviation at both outputs and inputs, are validated experimentally, see, e.g., [21, 33, 130, 131]. It is also often observed that explicit and implicit algorithms behave similarly for very small sampling times, though no clear explanation of this fact has been given yet. The classical first-order sliding-mode controller [33, 130, 131], the twisting controller [130], several higher order and homogeneous differentiators have been implemented on an electromechanical system [90, 124, 38] and compared. Homogeneous differentiators have been tested on an RLC circuit [132], differentiators have been implemented on DC-DC buck converters in [133]. Various discretizations of the classical first-order sliding-mode control [21], a switching input upgrading linear control [134], several higher-order differentiators [90, 124, 135] have been tested and compared in closed-loop on the inverted pendulum and the rotary inverted pendulum (including cascaded first-order differentiators). Set-valued admittance controllers [136, 107, 137, 138, 139], first-order differentiators [94],
<figure><img src="image_10.png" /><figcaption>Figure 7: Norms and total variations of state and control signals for several discrete-time super-twisting controllers. In terms of the $L_2$-norm of the state $x_0$ and the control input $u$, the implicit controller (217) (and its associated splittings); the semi-implicit controller (218); and the numerical solver for (220), have a similar behavior for all sampling times. In terms of the step-to-step variation (221) of the control input sequence, the explicit controller shows the larger values. It is noteworthy that for sampling times $h < 2\text{ms}$, the numerical solver presents an increment in the variation of the controller $u$, indicating the presence of numerical chattering. In terms of the $\|x_0\|_{L_\infty(15,20)}$, the controller (219) shows the best performance for $h \le 200\text{ms}$.</figcaption></figure>

```python
import numpy as np

def bracket(x, a):
    return np.sign(x) * np.abs(x) ** a

def J_B0(x, h_c, g0):
    beta = (-h_c*g0 + np/sqrt((h_c*g0)**2 +
    4.0*np.abs(x)))/2.0
    return (beta**2) * np.sign(x)

def J_B1(x, a):
    return x if np.abs(x) <= a else a * np.sign(x)

def sta_cdt_vu(h_c, g0, g1, x_0_k, nu, x_0_tilde, mu,
               N, a1, a2):
    z = x_0_tilde[-1]
    m = mu[-1]
    for i in range(N):
        z_next = J_B0((1 - a1)*z + a1*(x_0_k +
        h_c*nu[-1] - m, h_c, g0)
        m = J_B1(m + a2*(2.0*z_next - z), g1+h_c**2)
        z = z_next
    x_0_next_tilde = z
    u_k = -(x_0_k - x_0_next_tilde)/h_c
    nu_next = u_k + g0*bracket(x_0_next_tilde, 0.5)
    x_0_tilde.append(x_0_next_tilde)
    nu.append(nu_next)
    mu.append(m)
    return u_k
```

# 7 EXPERIMENTAL RESULTS AND VALIDATIONS

## 7.1 A Quick Review

Implicit and semi-implicit control and differentiation algorithms have been experimentally tested, mostly on laboratory setups. It is noteworthy that the fundamental properties like insensitivity with respect to the gain during sliding phases (stemming from (77), (98), (137)), moderate deterioration of the performance with increasing sampling times, significant chattering alleviation at both outputs and inputs, are validated experimentally, see, e.g., [21, 33, 130, 131]. It is also often observed that explicit and implicit algorithms behave similarly for very small sampling times, though no clear explanation of this fact has been given yet. The classical first-order sliding-mode controller [33, 130, 131], the twisting controller [130], several higher order and homogeneous differentiators have been implemented on an electromechanical system [90, 124, 38] and compared. Homogeneous differentiators have been tested on an RLC circuit [132], differentiators have been implemented on DC-DC buck converters in [133]. Various discretizations of the classical first-order sliding-mode control [21], a switching input upgrading linear control [134], several higher-order differentiators [90, 124, 135] have been tested and compared in closed-loop on the inverted pendulum and the rotary inverted pendulum (including cascaded first-order differentiators). Set-valued admittance controllers [136, 107, 137, 138, 139], first-order differentiators [94],
Consider the following problem. Given two maximal monotone operators $M_1, M_2: \mathbb{R}^n \to \mathbb{R}^n$, find $\xi \in \mathbb{R}^n$ such that

$$0 \in M_1(\xi) + M_2(\xi). \qquad (81)$$

Under constraint-qualification assumptions guaranteeing that $M_1 + M_2$ is maximal monotone, it follows from (23) in "Sidebar V" that problem (81) is equivalent to finding $\xi \in \mathbb{R}^n$ such that

$$\xi = \mathcal{J}_{\rho(M_1 + M_2)}(\xi) \qquad (82)$$

for any $\rho > 0$. The formulation (82) allows for a search of a solution to (81) in an iterative way. Indeed, it follows from the contents of "Sidebar VII" that the fixed-point iteration

$$\xi_{k+1} = \mathcal{J}_{\rho(M_1 + M_2)}(\xi_k) \qquad (83)$$

converges towards a solution of (82). However, an explicit expression for the resolvent map of the sum $M_1 + M_2$ in (83) may not be available. To solve such issue, proximal splitting schemes are applied to compute alternative iterations that converge to a solution of (81) and that only involve the individual resolvents of $M_1(\cdot)$ and $M_2(\cdot)$ at each iteration. In what follows, two popular splitting schemes are discussed: the Condat-Vū and the Douglas-Rachford splittings.

### 6.5 Condat-Vū splitting

The Condat-Vū splitting was independently proposed in [S43] and [S44] for solving the optimization problem

$$\min_{\xi \in \mathbb{R}^n} f(\xi) + g(A\xi) + h(\xi), \qquad (84)$$

where $f: \mathbb{R}^n \to \mathbb{R}$, $g: \mathbb{R}^m \to \mathbb{R}$, and $h: \mathbb{R}^n \to \mathbb{R}$ are proper, convex, LSC functions. The function $h$ is continuously differentiable and $A \in \mathbb{R}^{m \times n}$. Recalling that the proximal map is itself the solution of a convex optimization problem (see (51) and (53) in "Sidebar VII"), the Condat-Vū splitting can be used for approximating the values of elaborated proximal maps such as those appearing in higher-order sliding-mode controllers, observers and differentiators (see for instance (127) and (148)). The splitting relies on the fact that $g(A\xi) = \sup_{\mu \in \mathbb{R}^m} \{\langle A\xi, \mu \rangle - g^*(\mu)\}$ (Fenchel-Moreau theorem [S28, Theorem 13.32]) so that, as pointed out in [S43], the primal problem (84) can be rewritten as

$$\min_{\xi \in \mathbb{R}^n} \sup_{\mu \in \mathbb{R}^m} \mathcal{L}(\xi, \mu), \qquad (85)$$

where $\mathcal{L}: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ is of the form

$$\mathcal{L}(\xi, \mu) = f(\xi) + h(\xi) + \langle \mu, A\xi \rangle - g^*(\mu). \qquad (86)$$

Thus, the saddle points of $\mathcal{L}(\cdot)$ are characterized by the following generalized equations:

$$0 \in \nabla h(\xi) + A^{\top}\mu + \partial f(\xi), \qquad (87a)$$

$$0 \in -A\xi + \partial g^*(\mu). \qquad (87b)$$

Proceeding in a similar way as above, (87) is equivalent to the following fixed-point conditions:

$$\xi = \text{Prox}_{\alpha_1 f} (\xi - \alpha_1 (\nabla h(\xi) + A^{\top} \mu)), \qquad (88a)$$

$$\mu = \text{Prox}_{\alpha_2 g^*} (\mu + \alpha_2 A\xi), \qquad (88b)$$

for any $\alpha_1, \alpha_2 > 0$. It is well-known that, if $\xi^*$ and $\mu^*$ are solutions of the primal problem (84) and its associated dual problem, respectively, then they also satisfy (88), see e.g., [S28, Proposition 19.18]. The converse is true if, in addition, a constraint qualification such as $0 \in \text{rint}(A \text{dom } f - \text{dom } g)$ holds. Under such total duality assumption, the following iteration was proposed in [S43, S44]:

$$\xi_{j+1} = \text{Prox}_{\alpha_1 f} (\xi_j - \alpha_1 (\nabla h(\xi_j) + A^{\top} \mu_j)), \qquad (89a)$$

$$\mu_{j+1} = \text{Prox}_{\alpha_2 g^*} (\mu_j + \alpha_2 A(2\xi_{j+1} - \xi_j)). \qquad (89b)$$
Consider the following problem. Given two maximal monotone operators $M_1, M_2: \mathbb{R}^n \to \mathbb{R}^n$, find $\xi \in \mathbb{R}^n$ such that

$$0 \in M_1(\xi) + M_2(\xi). \qquad (81)$$

Under constraint-qualification assumptions guaranteeing that $M_1 + M_2$ is maximal monotone, it follows from (23) in “Sidebar V” that problem (81) is equivalent to finding $\xi \in \mathbb{R}^n$ such that

$$\xi = \mathcal{J}_{\rho(M_1 + M_2)}(\xi) \qquad (82)$$

for any $\rho > 0$. The formulation (82) allows for a search of a solution to (81) in an iterative way. Indeed, it follows from the contents of “Sidebar VII” that the fixed-point iteration

$$\xi_{k+1} = \mathcal{J}_{\rho(M_1 + M_2)}(\xi_k) \qquad (83)$$

converges towards a solution of (82). However, an explicit expression for the resolvent map of the sum $M_1 + M_2$ in (83) may not be available. To solve such issue, proximal splitting schemes are applied to compute alternative iterations that converge to a solution of (81) and that only involve the individual resolvents of $M_1(\cdot)$ and $M_2(\cdot)$ at each iteration. In what follows, two popular splitting schemes are discussed: the Condat-Vū and the Douglas-Rachford splittings.

## 6.5 Condat-Vū splitting

The Condat-Vū splitting was independently proposed in [S43] and [S44] for solving the optimization problem

$$\min_{\xi \in \mathbb{R}^n} f(\xi) + g(A\xi) + h(\xi), \qquad (84)$$

where $f: \mathbb{R}^n \to \mathbb{R}$, $g: \mathbb{R}^m \to \mathbb{R}$, and $h: \mathbb{R}^n \to \mathbb{R}$ are proper, convex, LSC functions. The function $h$ is continuously differentiable and $A \in \mathbb{R}^{m \times n}$. Recalling that the proximal map is itself the solution of a convex optimization problem (see (51) and (53) in “Sidebar VII”), the Condat-Vū splitting can be used for approximating the values of elaborated proximal maps such as those appearing in higher-order sliding-mode controllers, observers and differentiators (see for instance (127) and (148)). The splitting relies on the fact that $g(A\xi) = \sup_{\mu \in \mathbb{R}^m} \{\langle A\xi, \mu \rangle - g^*(\mu)\}$ (Fenchel-Moreau theorem [S28, Theorem 13.32]) so that, as pointed out in [S43], the primal problem (84) can be rewritten as

$$\min_{\xi \in \mathbb{R}^n} \sup_{\mu \in \mathbb{R}^m} \mathcal{L}(\xi, \mu), \qquad (85)$$

where $\mathcal{L}: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ is of the form

$$\mathcal{L}(\xi, \mu) = f(\xi) + h(\xi) + \langle \mu, A\xi \rangle - g^*(\mu). \qquad (86)$$

Thus, the saddle points of $\mathcal{L}(\cdot)$ are characterized by the following generalized equations:

$$0 \in \nabla h(\xi) + A^{\top}\mu + \partial f(\xi), \qquad (87a)$$

$$0 \in -A\xi + \partial g^*(\mu). \qquad (87b)$$

Proceeding in a similar way as above, (87) is equivalent to the following fixed-point conditions:

$$\xi = \text{Prox}_{\alpha_1 f} (\xi - \alpha_1 (\nabla h(\xi) + A^{\top}\mu)), \qquad (88a)$$

$$\mu = \text{Prox}_{\alpha_2 g^*} (\mu + \alpha_2 A\xi), \qquad (88b)$$

for any $\alpha_1, \alpha_2 > 0$. It is well-known that, if $\xi^*$ and $\mu^*$ are solutions of the primal problem (84) and its associated dual problem, respectively, then they also satisfy (88), see e.g., [S28, Proposition 19.18]. The converse is true if, in addition, a constraint qualification such as $0 \in \text{rint}(A \text{dom } f - \text{dom } g)$ holds. Under such total duality assumption, the following iteration was proposed in [S43, S44]:

$$\xi_{j+1} = \text{Prox}_{\alpha_1 f} (\xi_j - \alpha_1 (\nabla h(\xi_j) + A^{\top}\mu_j)), \qquad (89a)$$

$$\mu_{j+1} = \text{Prox}_{\alpha_2 g^*} (\mu_j + \alpha_2 A(2\xi_{j+1} - \xi_j)). \qquad (89b)$$
## Sidebar (cont.)

Thus, if the gradient of $h(\cdot)$ is Lipschitz continuous with constant $\ell$, and

$$ \frac{\alpha_1 \ell}{2} + \alpha_1 \alpha_2 \lambda_{\max}(A^\top A) < 1, \qquad (90) $$

then the iteration (89) converges towards solutions of (88), implying that the sequence $\{\xi_k\}_{k \in \mathbb{N}}$ converges to a solution of (84).

## 6.6 Douglas-Rachford splitting

Similar to the previous approach, the Douglas-Rachford splitting can be used to approximate proximal maps composed by the sum of two maximal monotone operators. Let us start considering the problem of finding $\xi \in \mathbb{R}^n$ such that

$$ 0 \in \partial f(\xi) + \partial g(\xi), \qquad (91) $$

where, as before, $f(\cdot)$ and $g(\cdot)$ are proper, convex, LSC functions. Thus, as pointed out in [S35, p. 44] the following equivalence is true:

$$ 0 \in \partial f(\xi) + \partial g(\xi) \qquad \Leftrightarrow z = (\mathcal{R}_{\alpha \partial f} \circ \mathcal{R}_{\alpha \partial g})(z), \text{ and } \xi = \text{Prox}_{\alpha g}(z), \qquad (92) $$

where $\mathcal{R}_{\alpha \partial f}(\cdot) = (2 \text{Prox}_{\alpha f} - \text{I}_d)(\cdot)$ is the reflected resolvent. Now, since the operator $\mathcal{R}_{\alpha \partial f} \circ \mathcal{R}_{\alpha \partial g}$ is simply non-expansive, a Picard iteration using (92) may not converge. To overcome such issue, the Douglas-Rachford splitting uses the following Krasnosel'ski-Mann iteration

$$ z_{j+1} = \frac{1}{2}z_j + \frac{1}{2}(\mathcal{R}_{\alpha \partial f} \circ \mathcal{R}_{\alpha \partial g})(z_j) \\ = \text{Prox}_{\alpha f} \circ (2 \text{Prox}_{\alpha g} - \text{I}_d)(z_j) + (\text{I}_d - \text{Prox}_{\alpha g})(z_j), \qquad (93) $$

for which its convergence is guaranteed for any $\alpha > 0$ [129, Theorem 1]. Thus, the iteration

$$ z_{j+1} = \text{Prox}_{\alpha f} \circ (2 \text{Prox}_{\alpha g} - \text{I}_d)(z_j) + (\text{I}_d - \text{Prox}_{\alpha g})(z_j) \qquad (94a) $$

$$ \xi_{j+1} = \text{Prox}_{\alpha g}(z_j) \qquad (94b) $$

converges and $\xi_{j+1}$ approaches a solution of (91) as $j \uparrow +\infty$.

## References

[S43] L. Condat, A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms. *Journal of Optimization Theory and Applications*, 158:460-479, 2013.

[S44] C.B. Vu, A splitting algorithm for dual monotone inclusions involving cocoercive operators. *Advances in Computational Mathematics*, 38:667-681,2013.
## Sidebar (cont.)

Thus, if the gradient of $h(\cdot)$ is Lipschitz continuous with constant $\ell$, and

$$ \frac{\alpha_1 \ell}{2} + \alpha_1 \alpha_2 \lambda_{\max}(A^\top A) < 1, \qquad (90) $$

then the iteration (89) converges towards solutions of (88), implying that the sequence $\{\xi_k\}_{k \in \mathbb{N}}$ converges to a solution of (84).

## 6.6 Douglas-Rachford splitting

Similar to the previous approach, the Douglas-Rachford splitting can be used to approximate proximal maps composed by the sum of two maximal monotone operators. Let us start considering the problem of finding $\xi \in \mathbb{R}^n$ such that

$$ 0 \in \partial f(\xi) + \partial g(\xi), \qquad (91) $$

where, as before, $f(\cdot)$ and $g(\cdot)$ are proper, convex, LSC functions. Thus, as pointed out in [S35, p. 44] the following equivalence is true:

$$ 0 \in \partial f(\xi) + \partial g(\xi) $$

$$ \Leftrightarrow z = (\mathcal{R}_{\alpha \partial f} \circ \mathcal{R}_{\alpha \partial g})(z), \text{ and } \xi = \text{Prox}_{\alpha g}(z), \qquad (92) $$

where $\mathcal{R}_{\alpha \partial f}(\cdot) = (2 \text{Prox}_{\alpha f} - \text{I}_d)(\cdot)$ is the reflected resolvent. Now, since the operator $\mathcal{R}_{\alpha \partial f} \circ \mathcal{R}_{\alpha \partial g}$ is simply non-expansive, a Picard iteration using (92) may not converge. To overcome such issue, the Douglas-Rachford splitting uses the following Krasnosel'ski-Mann iteration

$$ \begin{aligned} z_{j+1} &= \frac{1}{2}z_j + \frac{1}{2}(\mathcal{R}_{\alpha \partial f} \circ \mathcal{R}_{\alpha \partial g})(z_j) \\ &= \text{Prox}_{\alpha f} \circ (2 \text{Prox}_{\alpha g} - \text{I}_d)(z_j) + (\text{I}_d - \text{Prox}_{\alpha g})(z_j), \end{aligned} \qquad (93) $$

for which its convergence is guaranteed for any $\alpha > 0$ [129, Theorem 1]. Thus, the iteration

$$ z_{j+1} = \text{Prox}_{\alpha f} \circ (2 \text{Prox}_{\alpha g} - \text{I}_d)(z_j) + (\text{I}_d - \text{Prox}_{\alpha g})(z_j) \qquad (94a) $$

$$ \xi_{j+1} = \text{Prox}_{\alpha g}(z_j) \qquad (94b) $$

converges and $\xi_{j+1}$ approaches a solution of (91) as $j \uparrow +\infty$.

## References

[S43] L. Condat, A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms. *Journal of Optimization Theory and Applications*, 158:460–479, 2013.

[S44] C.B. Vu, A splitting algorithm for dual monotone inclusions involving cocoercive operators. *Advances in Computational Mathematics*, 38:667–681, 2013.
have been tested on robotic systems. The super-twisting observer and sliding mode observers have been validated in silico on pancreas and diabetes control, and glucose regulation with insulin infusion in [140, 141, 142]. A parabolic differentiator is implemented on a direct current servomotor equipped with an optical encoder in [106]. A set-valued anti-sway controller is applied on an overhead crane in [143]. A super-twisting observer is applied to a magneto-rheological clutch in [144, 145], while a super-twisting controller semi-implicitly discretized is applied to an electromechanical system with time-delay in [146]. Various controllers (super-twisting, first-order) are applied on a floating wind turbine in [147, 148, 149], with extensive comparisons with other types of controllers. The super-twisting controller is applied to functional electrical stimulation of index finger's metacarpophalangeal in [150]. A first-order sliding-mode set-valued input is designed for the stabilization of a nonsmooth model of an hydraulic actuator in [36] and applied to an industrial commercial excavator [151, 152, 153]. A toolbox for automatic gain-tuning of higher-order differentiators implemented with the backward-Euler method as those introduced in section 4.2 is proposed in [154, 84] (see also [155] for a toolbox dedicated to differentiators with the forward-Euler discretization).

## 7.2 Typical Experimental Results

Let us present a few experimental results which have been obtained in [120, 156, 21, 33, 130, 157, 131] on two setups: an electropneumatic system, and an inverted pendulum. Two controllers are tested: the classical first-order SMC and the twisting algorithm. They nicely illustrate some important properties of backward-Euler schemes. The aim is to illustrate the main specific features without presenting details which are available in the cited articles. The first property is digital chattering alleviation and magnitude decrease in the input, see Fig. 8 for first-order SMC applied to the electropneumatic setup [33] (see [158] for a video showing the setup). The explicit controller (Fig. 8 top) switches between the two extremes of the sign function at each sampling step: this is a high-frequency bang-bang signal. By contrast, the values of the implicit controller are selected from the graph of the sign set-valued map at each sampling step. Consequently, when the sliding surface has been reached, the amplitude of the explicit controller depends on the feedback gain, whereas the amplitude of the implicit controller is independent of such gain and it only depends on the disturbance (in this case a square force signal, see Fig. 8 bottom): we are therefore exactly in a situation where (98) (or (77) for the passivity-based first-order SMC) hold. Fig. 9 and 10 depict experiments on the same electropneumatic system, using the twisting algorithm for a trajectory tracking task [120, 130, 102]. It is seen in Fig. 9 that the explicit implementation performance deteriorates a lot when $h$ increases, while the implicit one keeps

a very good level of performance. In Fig. 10 the gain insensitivity property (as a consequence of (21) and (98)) is shown for three different sampling times $h$, and a large range of control gains (equivalent to $K$ in (19)). Very similar results on gain insensitivity are obtained in [21] (first-order SMC and inverted pendulum) and [33] (first-order SMC and electropneumatic system). Fig. 11 and 12 show the evolution of the system's output for two different sampling times (first-order SMC, electropneumatic system [33]). Explicit and implicit controllers produce similar outputs for very small $h$ (Fig. 11 and 12, top figures). However it is striking to notice that although the explicit method displays large chattering for $h = 15$ms, the implicit one remains with very low chattering for both $h = 2$ms and $h = 15$ms. This confirms the results depicted in Fig. 9. This also demonstrates that both methods yield similar outputs for very small $h$, however the inputs are still quite different one from each other (they always look like in Fig. 8). Quite similar results are obtained with a first-order SMC applied to a six-component thrust generator installed on the head of a floating wind turbine [148, Fig. 7-10], where codimension $n \ge 2$ sliding surfaces are considered.

## 8 CONCLUSIONS AND PERSPECTIVES

The implicit (or, backward-Euler) discretization of set-valued sliding mode systems has received significant attention during the last decade, since it was shown in [31] that it applies to systems with unknown matched disturbances and is an efficient way to significantly alleviate the digital chattering at both the output and the input. This article aims to show the close link between implicit/semi-implicit discretization of set-valued sliding-mode control and differentiation, and optimization. In particular, proximal-point algorithms using the resolvents of maximal monotone operators are shown to be at the core of the analysis and of the implementation of implicit methods. The crucial role played by passivity in the design of implicit discretization of first-order sliding-mode algorithms is highlighted. Robust, time-varying, Hamiltonian and higher-order versions of the proximal-point algorithm are introduced. Computational issues, which are crucial for real-time implementations, are analyzed. Splitting and relaxation algorithms, which provide systematic ways to solve the generalized equations associated with the one-step problems, are reviewed. Some numerical experiments illustrate the presentation.

It is the authors' belief that the analysis of this class of algorithms is still in its infancy, but existing results and experiments show its great potential. The implicit discretizations of infinite-dimensional systems (maximal monotonicity being well suited to the study of some classes of partial differential equations), output feedback control, fixed-time and prescribed-time convergent systems, are still largely open problems. The problem's complexity (which can stem from the sliding surface codi-