# Analysis of a Dry Friction Force Law for the Covariant Optimal Control of Mechanical Systems with Revolute Joints

Juan Antonio Rojas-Quintero, François Dubois, Hedy César

Ramírez-De-Ávila, Eusebio Bugarin, Bruno Sánchez-García, Nohe R

Cazarez-Castro

► To cite this version:

Juan Antonio Rojas-Quintero, François Dubois, Hedy César Ramírez-De-Ávila, Eusebio Bugarin, Bruno Sánchez-García, et al.. Analysis of a Dry Friction Force Law for the Covariant Optimal Control of Mechanical Systems with Revolute Joints. Mathematics, 2024, 12 (20), pp.3239. 10.3390/math12203239. hal-05335037

HAL Id: hal-05335037

https://hal.science/hal-05335037v1

Submitted on 28 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Analysis of a Dry Friction Force Law for the Covariant Optimal Control of Mechanical Systems with Revolute Joints

Juan Antonio Rojas-Quintero <fnref n="1" />\*, François Dubois <fnref n="2" />,<fnref n="3" />\*, Hedy César Ramírez-de-Ávila <fnref n="4" />\*, Eusebio Bugarin <fnref n="5" />\*,
Bruno Sánchez-García <fnref n="4" />\* and Nohe R. Cazarez-Castro <fnref n="4" />,<fnref n="6" />\*

**Abstract:** This contribution shows a geometric optimal control procedure to solve the trajectory generation problem for the navigation (generic motion) of mechanical systems with revolute joints. The mechanical system is analyzed as a nonlinear Lagrangian system affected by dry friction at the joint level. Rayleigh's dissipation function is used to model this dissipative effect of joint-level friction, and regarded as a potential. Rayleigh's potential is an invariant scalar quantity from which friction forces derive and are represented by a smooth model that approaches the traditional Coulomb's law in our proposal. For the optimal control procedure, an invariant cost function is formed with the motion equations and a Riemannian metric. The goal is to minimize the consumed energy per unit time of the system. Covariant control equations are obtained by applying Pontryagin's principle, and time-integrated using a Finite Elements Method-based solver. The obtained solution is an optimal trajectory that is then applied to a mechanical system using a proportional-derivative plus feedforward controller to guarantee the trajectory tracking control problem. Simulations and experiments confirm that including joint-level friction forces at the modeling stage of the optimal control procedure increases performance, compared with scenarios where the friction is not taken into account, or when friction compensation is performed at the feedback level during motion control.

**Keywords:** Lagrangian systems; optimal control; Rayleigh's dissipation function; friction force laws; nonlinear mechanical systems; Riemannian geometry; Pontryagin's maximum principle

**MSC:** 49S05; 51P05; 53A35; 70E60

## 1. Introduction

Friction is a resistive and non-conservative force. It arises when there is contact between two surfaces. Therefore, it is present in most realistic systems because mechanical interfaces are present. Being resistive, it systematically opposes the motion, and because it is a non-conservative force, some of the mechanical energy that contributes to the system motion transforms into heat, noise, and material deformation and wear in some extreme cases. These energy losses result in motion control performance decline as vibration propagates and the system is dragged down in a sometimes hardly predictable manner.

This study involves an actuated nonlinear mechanical system in which friction plays a significant role during motion control. Friction drags down the system, introduces

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

check for
updates

**Citation:** Rojas-Quintero, J.A.; Dubois, F.; Ramírez-de-Ávila, H.C.; Bugarin, E.; Sánchez-García, B.; Cazarez-Castro, N.R. Analysis of a Dry Friction Force Law for the Covariant Optimal Control of Mechanical Systems with Revolute Joints. Mathematics **2024**, *12*, 3239. https://doi.org/10.3390/math12203239

Academic Editor: Andrea Scozzari

Received: 17 September 2024

Revised: 12 October 2024

Accepted: 14 October 2024

Published: 16 October 2024

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

**Copyright:** © 2024 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

<footnote n="1">¹ CONAHCYT—Tecnológico Nacional de México, I. T. Ensenada, Ensenada 22780, BC, Mexico</footnote>

<footnote n="2">² Laboratoire de Mathématiques d'Orsay, Université Paris-Saclay, 91400 Orsay, France;
francois.dubois@universite-paris-saclay.fr</footnote>

<footnote n="3">³ Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et Métiers,
75141 Paris, France</footnote>

<footnote n="4">⁴ Tecnológico Nacional de México, I. T. Tijuana, Tijuana 22414, BC, Mexico;
d22210006@tectijuana.edu.mx (H.C.R.-d.-Á.); d24210005@tectijuana.edu.mx (B.S.-G.);
nohe@ieee.org (N.R.C.-C.)</footnote>

<footnote n="5">⁵ Tecnológico Nacional de México, I. T. Ensenada, Ensenada 22780, BC, Mexico; ebugarin@ite.edu.mx</footnote>

<footnote n="6">⁶ Facultad de Ingeniería Química, Universidad Michoacana de San Nicolas de Hidalgo,
Morelia 58030, MC, Mexico</footnote>

<footnote>* Correspondence: jarojas@conahcyt.mx</footnote>
undesirable stick-slip motions, and results in positioning errors. As a result, appropriately modeling and predicting the effects of friction on dynamical systems is highly desirable. A set of reviews covering a large spectrum of issues in modeling, numerical simulation, and control of dynamical systems subject to friction effects can be found in [1–7]. These provide a detailed overview of the main problems and advances in the field.

Our focus is on the analysis, numerical and experimental, of a smooth friction model that approximates friction forces at the joints of a mechanical system when used in the framework of the indirect optimal control (OC) approach proposed in [8,9]. This OC method is a Riemannian formulation of Pontryagin's Maximum Principle (PMP). It is covariant and holds interesting properties when compared to more traditional formulations [9]: equation stiffness is reduced, leading to shorter computing times, and resulting motions require less power consumption and less motion amplitude. This methodology is a geometric description of optimal control that can be applied to different types of multibody systems such as robotic manipulators and mechatronic systems of the sort. Note that OC problems often lead to solving a stiff set of nonlinear ordinary differential equations (ODEs) [10–13]. Therefore, adding friction forces into the motion equations shall make the task even more challenging. Indeed, classical friction models usually showcase discontinuity at null velocity [2]. Therefore, special attention is required to solve the ODE arising from this OC methodology.

We model friction at the system joints by invoking Coulomb's law, one of the most fundamental friction models. More precise and complete friction models exist. The Dahl and LuGre models [7] are two examples that are widely used for friction compensation (see, e.g., [3]). Nevertheless, Coulomb's model is considered to be a suitable first step [14] because it has led to successful practical applications. However, Coulomb's model does feature a discontinuity at zero velocity (see Figure 1a) because it is usually defined using the signum function. Non-smooth friction approximations such as the Dahl and LuGre models involve an internal variable and an additional restriction. The numerical treatment of the classical Coulomb friction model requires specific procedures deal with the discontinuity at zero velocity. We approach Coulomb friction by a smooth approximation that replaces the signum function, with the continuously differentiable hyperbolic tangent function. This smooths out the discontinuity at zero velocity [5,7] (see Figure 1b) and avoids the need for a specific procedure or auxiliary equations. Being an approximation, one should expect limitations, but our experimental observations prove that it could be a good starting point.

Direct optimization procedures, (where the dynamics are first discretized and then optimized) for OC, routinely deal with friction by making use of complementarity formulations of Coulomb's model as in [15–17] for example. However, this is usually carried out to deal with friction forces that come from contact or impact situations with external objects in the system. Our approach focuses on joint-level frictional forces (that are internal to the system) within an indirect method (where we first optimize the dynamics and then discretize them). Direct optimization procedures fall out of the scope of this work.

Regarding indirect OC methods, joint-level friction effects (internal to the system) have previously been considered in robotics for example. However, either viscous friction is exclusively considered [18–21], or friction dynamics are linearized through the Linear Quadratic Regulator framework [18,22]. Indeed, joint-level Coulomb's friction models are rarely associated with the nonlinear OC of dynamical systems in the modeling phase of joint motion. One of those few examples can be found in [23] where an indirect optimal control method based on Pontryagin's Maximum Principle (PMP) was used to optimally control the motions of an industrial robot while previously identifying a friction model for its joints. Note that even though the identified friction model combines Coulomb's and viscous friction effects, the authors report that it is continuous and specific to a particular robot model. Another example can be found in [24], where Dynamic Programming is used, but it is unspecified if a traditional Coulomb model or a continuous approximation was selected. The work by Kim and Ha [25] is perhaps most closely related to our work. There, an explicit solution is given for the time-optimal control of a one-degree-of-freedom (DOF)
mechanical system. Our work differs in that actuating forces at the joint level are being minimized. Either way, no mention is usually made about the numerical complexities that arise when dealing with an internal, joint-level Coulomb's friction force in conjunction with OC-related equations. It seems that such friction effects have been routinely regarded as a purely control issue where its forces require compensation at the feedback level [3].

<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Classical Coulomb friction model. (b) Coulomb-tanh alternative model. The alternative Coulomb-tanh model ensures continuity at null velocity (represented by the $\dot{q}$ variable) and approaches the classical model as $k$ increases.</figcaption></figure>

Consequently, there is motivation to manage friction at the modeling level. In this way, the optimal control process shall lead to trajectories that inherently minimize friction effects during motion. We hypothesize that optimal control performance increases by including joint-level friction into the motion equations of the system. This hypothesis will be experimentally verified. We shall see that this is still a difficult task because even if the selected friction model (that we refer to as “Coulomb-tanh”) is continuous, equations become stiffer when the slope coefficient within the hyperbolic tangent function increases. Consequently, the used ODE solver is revealed to be rapidly affected by this coefficient when set up to operate using traditional algorithms such as shooting methods and finite differences. The optimal control equations are solved using a commercial ODE solver set up to operate with a Finite Element Method (FEM) that discretizes over the time dimension, to overcome such limitations. Additionally, we will experimentally demonstrate that the hyperbolic tangent slope coefficient does not require reaching excessively high values to obtain acceptable performance in our OC process that considers internal joint-level friction effects.

The goal of the present work is not to propose a numerical method designed to tackle our optimal control problem. Instead, we present some first steps toward incorporating joint-level friction forces into the system model, in a way that is compatible with a Riemannian formulation, and such that these forces are directly taken into account by the indirect optimal control method. Therefore, we analyze the impact and the benefits of incorporating joint-level friction forces into the model. We compare our results with a scenario where these forces are compensated at the feedback level. The main contributions of our work consist of showing the following:

* The comprehensive integration of joint- level friction forces into the motion equations through Rayleigh's dissipation function;

* Numerical evidence that the smooth friction approximation converges toward the classical Coulomb friction within the optimal control procedure;

* Numerical confirmation that control performance is improved by including friction forces in the modeling stage, compared to a feedback-based friction compensation;

* Experimental validation of the optimal control procedure (with satisfactory results even in an open-loop control).
The remainder of the manuscript is organized as follows. First, we present the OC methodology, starting with the Riemannian metric selection in Section 2. Joint-level friction will be incorporated into the model by ensuring that the covariance of motion equations is preserved. We will show that friction models for rotational joints can be derived from Rayleigh's dissipation function, regarded as a scalar potential. Since joint-level friction forces are included in the system dynamics, these are managed by the optimal control process at the motion planning stage. Section 3 is devoted to the numerical simulation analysis of our situation. The optimal path planning of a simplified nonlinear Lagrangian system considering joint friction will be taken as the case scenario. An example where the control variable is determined with the Riemannian OC method will expose the numerical difficulties that arise when including friction forces. The stiffness of the optimal control equations will be briefly assessed. Then, further simulations will show that the proposed procedure is compatible with multi-DOF systems. We will then realize the optimally controlled motion on a simple experimental setup built for this study. Section 4 presents these results. We will show that the selected friction model leads to a performance increase of the optimal control. Finally, we will discuss the obtained results in this work that we deem to be a first step toward more complex scenarios. We will also discuss perspectives to treat more precise friction models on more complex systems.

## 2. Riemannian Optimal Control

Let us consider a multibody mechanical system described by a finite number of $n$ configuration parameters $q^i$ where $i \in [1, n]$. In multibody systems such as these, a joint supports another one on the chain. An actuator ensures rotational joint motion by exerting torque. $Q$ designates the set of states $q = \{q^i\}$. In the remainder of the manuscript, we shall rely on tensor notation and apply Einstein's summation convention on repeated subscript-superscript indices [26–28]. This is a convention that is used to ease notation in tensor calculus. For example,

$$a_{kl}^j b_j^k = \sum_{j,k} a_{kl}^j b_j^k,$$

and the live index of the resulting quantity is $\ell$.

### 2.1. Riemannian Metric

System dynamics is derived from the kinetic energy $K$ of the system (as in [26]). It is a quadratic, strictly convex form of $\dot{q}^i$. An overdot denotes the time rate of change of the variable. $K$ coefficients define a positive-definite, configuration-dependent Hessian $M(q)$, called mass tensor. Nonlinear regular functions of the state $q \in Q$ compose these tensor components. The kinetic energy of the system is expressed as

$$K(q, \dot{q}) = \frac{1}{2} M_{kl}(q) \dot{q}^k \dot{q}^l, \qquad (1)$$

where the mass tensor $M(q)$ defines a metric in a Riemannian manifold [26,29–31] and Einstein's summation convention is implied over repeated subscript-superscript indices. The choice of a metric tensor is fundamental for calculating lengths of vectors and curves [28]. In the current setting, the selected metric stands above and so the following methodology lies in the framework of Riemannian geometry.

### 2.2. Riemannian Manifold

With the selected Riemannian metric, $Q$ has a Riemannian structure and the following manifold can be established.

(a) Covariant space derivation $\partial_j \equiv \frac{\partial}{\partial q^j}$.

(b) Contravariant space derivation $\partial^j$ such that $\langle \partial^j, \partial_k \rangle = \delta_k^j$.

(c) Components of the inverse mass tensor $M^{-1} = M^{jl}$ such that $M_{ij}M^{jl} = \delta_i^j$.
(d) The connection $\Gamma_{ik}^{\ell} = \frac{1}{2} M^{il} (\partial_i M_{lk} + \partial_k M_{li} - \partial_{\ell} M_{ik})$ is symmetric with respect to its lower indices: $\Gamma_{ki}^{j} = \Gamma_{ik}^{j}$. It defines the differentials $de_j = \Gamma_{jk}^{\ell} dq^k e_{\ell}$ and $de^j = -\Gamma_{k\ell}^{j} dq^k e^{\ell}$.

(e) The relationship between covariant components $\varphi_j$ and contravariant components $\varphi^k$ is established through the metric: $\varphi_j = M_{jk} \varphi^k$, $\varphi^k = M^{jk} \varphi_j$.

(f) Covariant differential of a vector field $\varphi = \varphi^i e_j$, $d\varphi = (\partial_\ell \varphi^j + \Gamma_{lk}^j \varphi^k) dq^\ell e_j$.

(g) Covariant differential of a covector field $\varphi = \varphi^i e^\ell$, $d\varphi = (\partial_k \varphi^\ell - \Gamma_{k\ell}^k \varphi^j) dq^k e^\ell$.

(h) Covariant derivative of a scalar field: $dV = \partial_\ell V dq^\ell = \langle \nabla V, dq^i e_j \rangle$ and $\nabla V = \partial_\ell V e^\ell$.

(i) Covariant derivative of a covector field $\varphi = \varphi^i e^\ell$, $d\varphi = \langle \nabla \varphi, dq^i e_j \rangle$, and $\nabla \varphi = (\partial_k \varphi^\ell - \Gamma_{k\ell}^j \varphi^j) e^k e^\ell$.

## 2.3. System Dynamics Subject to Dissipative Effects

We shall consider the configuration-dependent gravitational potential $V = V(q)$. It models gravity actions and is defined as a product between the system mass, the local gravitational field intensity, and the center of mass height. The system Lagrangian is defined as $L(q, \dot{q}) = K(q, \dot{q}) - V(q)$. When a forcing control $u_i$ is exerted at the joint, the Euler-Lagrange equations

$$ \frac{d}{dt} \frac{\partial L}{\partial \dot{q}^i} - \frac{\partial L}{\partial q^i} = u_i $$

are satisfied (see [31,32]).

The equations above describe the Lagrangian system dynamics (c.f. [33]) without considering energy dissipation. However, motion equations of most realistic multibody mechanical systems are not conservative. Such systems are generally affected by forces that are velocity-dependent and, ultimately, dissipate energy. These forces oppose motion and may be grouped in a non-conservative generalized forces term that derives from a velocity-dependent scalar potential $F \equiv F(\dot{q})$. This potential corresponds to Rayleigh's dissipation function ([32], ch. 1). Therefore, when accounting for dissipative effects at each joint, the Euler-Lagrange equations become

$$ \frac{d}{dt} \frac{\partial L}{\partial \dot{q}^i} - \frac{\partial L}{\partial q^i} = u_i - \frac{\partial F}{\partial \dot{q}^i}. \qquad (2) $$

Equation (2) is prescribed by some authors (see [34] (ch. 6) or [35] (ch. 4), for example) and is a generalization after the works of J. W. Strutt (otherwise known as Baron or Lord Rayleigh) [36], perhaps best advocated by [37]. Note that Strutt only considered the case of viscous friction that is linear in the system velocities. However, Rayleigh's dissipation function can be used to model Coulomb's friction (see [37]) and ensures covariance preservation of motion equations when regarded as an invariant potential function. Indeed, for an invariant $F(\dot{q})$, the derivative $\frac{\partial F}{\partial \dot{q}^i}$ is covariant.

Here, we propose to approximate Coulomb's friction at the system joints, using the Coulomb-tanh friction force model invoked in [4], so that

$$ \frac{\partial F}{\partial \dot{q}^i} = f_{ij} \tanh(k \dot{q}^j), \qquad (3) $$

where $k$ adjusts the transition slope between negative and positive velocity; $f_{ij}$, such that $f_{ij} = 0$ for $i \neq j$, are dry friction coefficients at the joints; and $\dot{q}^j$ are the configuration parameters' time derivatives.

The Coulomb-tanh model of Equation (3) is an approximation of the traditional Coulomb's model for dry friction, which is usually represented through the signum function [2,3,5]. Note that we wish to focus on Coulomb's friction, and evaluate how usable
the proposed model from Equation (3) is. Therefore, viscous friction is not treated in this section. As such, when k is large enough, Equation (3) approaches the traditional Coulomb friction model (see Figure 1). The Coulomb-tanh model is nonlinear but should be easier to simulate than the more classical model, as suggested in [4,5]. This is because any conventional solver is able to deal with the motion Equation (2), and includes the Coulomb-tanh model of Equation (3). The numerical implementation is straightforward.

Upon inserting Equations (1) and (3) into Equation (2), the covariant second order equations of motion are found:

$$M_{ij}(\ddot{q}^i + \Gamma^j_{\ell m} \dot{q}^\ell \dot{q}^m) + \nabla_i V + f_{ij} \tanh(k\dot{q}^j) = u_i. \quad (4)$$

The system motion is thus governed by Equation (4). By recalling that $u^j = M^{ij}u_i$, the contravariant components of the torques $u$ can also be formulated:

$$\ddot{q}^i + \Gamma^i_{j\ell} \dot{q}^j \dot{q}^\ell + \nabla^i V + M^{ij} f_{j\ell} \tanh(k\dot{q}^\ell) = u^i. \quad (5)$$

## 2.4. Optimal Control

Consider a dynamical system at an initial state $x_0 = x(0)$. The idea is to take it to a final state $x_1 = x(T)$ in a period $T$, such that the control action minimizes the convex, invariant, and Riemannian cost functional

$$J(u) = \int_{0}^{T} \gamma(u(q, \dot{q}, \ddot{q})) dt = \int_{0}^{T} \frac{1}{2} M_{ij} u^i u^j dt = \int_{0}^{T} \frac{1}{2} u_i u^i dt. \quad (6)$$

Being endowed with the metric $M$ in this Riemannian setting, the above cost functional contains the scalar product $\langle u, u \rangle \equiv \frac{1}{2} u_i u^i$. Torques $u$ are defined covariantly ($u_i$) by Equation (4) and contravariantly ($u^i$) by Equation (5), and belong to the set $U$ of controls $u \equiv \{u^i\}$. These are functions of the states $q \equiv \{q^i\}$, and their first and second derivatives. States and controls are assumed to belong to a restricted class of smooth and differentiable functions in the domain $[0, T]$. Controls are admissible if they satisfy the imposed boundary conditions [38]. The integrand $\gamma = \gamma(u)$ is a convex and invariant running cost function that confers a Riemannian structure to our problem and is assumed to be continuously differentiable in $t \in [0, T]$, $q \in Q$, and $u \in U$. The fact that it is convex eases the minimization process of Equation (6).

The selected cost functional described by Equation (6) is regarded as the performance index. For systems with revolute joints, this performance index has the units of power (N m s⁻¹ ≡ kg m² s⁻³ ≡ W). Minimizing (6) amounts to lower energy consumption across time.

The control action can be minimized by applying Pontryagin's Maximum Principle (PMP) [38–40], which introduces Lagrange multipliers, regarded as adjoint variables. With these new variables, we form a Hamiltonian function $H$ of the states and the adjoint variables such that

$$H(\dot{x}, \lambda, \gamma) = \lambda \dot{x} - \gamma.$$

The states $x = \{q, \zeta\} = \{q, \dot{q}\}$ and the adjoint variables $\lambda = \{p, \zeta\}$ are continuous and $\gamma$ is the running cost function from Equation (6).

An important physical interpretation can be given to the adjoint variable $\zeta$ by carrying out the analysis in Proposition 2 of [9]: it is equal to the control variable $u$. This short analysis can be carried out in the following few lines. For instance, using the cost functional (6), $H(\dot{x}, \lambda, \gamma)$ takes the following form:

$$H(\dot{x}, \lambda, \gamma) = \xi_j \left[ u^j - \Gamma^j_{\ell m} \zeta^\ell \zeta^m - M^{i\ell} \partial_\ell V - M^{i\ell} f_{\ell m} \tanh(k\zeta^m) \right] + p_j \zeta^j - \frac{1}{2} M_{j\ell} u^j u^\ell. \quad (7)$$
We then apply PMP by canceling the partial derivative of the above Hamiltonian with respect to the control variable $u^i$: $\frac{\partial H}{\partial u^i} = 0$. The outcome is that

$$ \xi_j = M_{j\ell} u^\ell = u_j. \qquad (8) $$

Therefore, $\xi$ is exactly equal to the selected control variable. This property results from the Riemannian structure of the Hamiltonian function.

By combining condition (8) with Equation (7), the optimal Hamiltonian function is obtained:

$$ H = p_j \xi^j + \xi_j \left[ -\Gamma^j_{\ell m} \xi^\ell \xi^m - M^{j\ell} \partial_\ell V - M^{j\ell} f_{\ell m} \tanh(k \xi^m) \right] + \frac{1}{2} M^{j\ell} \xi_j \xi_\ell. $$

Motion and adjoint equations derive from this Hamiltonian in a symplectic manner, so that states and adjoint variables are, respectively, described by the following first-order ODE:

$$ \dot{q}^j = \frac{\partial H}{\partial p_j}, \quad \dot{\xi}^j = \frac{\partial H}{\partial \xi_j}, \quad \dot{p}_j(t) = -\frac{\partial H}{\partial q^j}, \quad \dot{\xi}_j(t) = -\frac{\partial H}{\partial \xi^j}. \qquad (9) $$

Note that the first two equations of (9) describe the system dynamics as a first-order system. The last two equations of (9) constrain the control variable behavior (see Proposition 3 of [9]).

In summary, the optimization procedure consists of solving the first-order system (9), subject to a set of four boundary conditions on the states. In our case, boundary conditions will be applied to initial and final positions and velocities.

### 3. Optimal Control Simulation of a Nonlinear Lagrangian System Considering Friction

Let us consider the actuated mechanical system depicted in Figure 2. It can be assimilated to a nonlinear pendulum, actuated at its center of rotation O through a rotating motor. Following from Equation (4), its motion can be approximated through the equation

$$ M(\ddot{q} + \omega^2 \sin(q)) + f_c \tanh(k\xi) = u, \qquad (10) $$

where M is the system mass tensor, and the $\omega$ parameter embeds gravity effects, system mass, and link length properties. Parameters $f_c$ and $k$ are, respectively, the Coulomb (or dry) friction coefficient and the hyperbolic tangent slope coefficient.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Actuated mechanical system. Motion is generated around the center of rotation O by an actuator delivering a forcing torque u. The generalized coordinate q is the angular position. C is the position of the center of mass. Vector g indicates the gravity action direction.</figcaption></figure>

We will now present the numerical analysis of this system optimal motion simulation. Our computations were performed on Wolfram's Mathematica® scientific computing software.
ware (version 12.3) [41]. It provides the differential equations solving tool NDSolveValue, which can be set up to operate with a wide range of numerical methods. The results presented in this section were obtained with this solver.

### 3.1. Optimal Control Simulation Taking Friction into Account

To optimally control the motions of this system, we follow the methodology presented in [9,42]. As such, to minimize the cost functional of the system torques (6), we describe the system dynamics and optimal controls by Equation (9). These become

$$ egin{cases} \dot{q} = \zeta \\ \dot{\zeta} = -\omega^2 \sin(q) - \frac{f_c \tanh(k\zeta) - \zeta}{M} \\ \dot{p} = \omega^2 \zeta \cos(q) \\ \dot{\zeta} = -p + \frac{f_c}{M} k \zeta \sech^2(k\zeta). \end{cases} \quad (11) $$

This scenario does not correspond to a typical friction compensation problem. It is an optimal motion planning procedure, where friction is managed at the modeling level. Therefore, the optimization process provides trajectories that minimize joint-level friction effects during motion.

It has been established through numerical trials that this type of ODE system can be difficult to solve by conventional numerical methods because it is subject to numerical stiffness (see [9,13]). Regarding stiffness, it may be safe to assume that the friction model will only exacerbate the issue [14]. Recall that we use a friction model that approaches Coulomb's friction when the parameter $k$ is sufficiently large (see Figure 1). Thus, when fixing the rest of the parameters, we can expect stiffness to increase with a growing value of $k$. Therefore, we wish to reach the highest value of $k$, for which Equation (11) can still be solved.

In the following numerical experiments of optimally controlled motion, we will solve the first-order system (11) subject to the parameters and boundary values of Table 1. Other boundary conditions were used in our trials; however, the following results are presented without loss of generality. Different boundary conditions change the trajectory but not the resulting motion properties.

<table><caption>**Table 1.** Simulation parameters for the optimal motion simulation of a one-DOF system subject to joint-level friction.</caption><thead><tr><th>Final Time</th><th colspan="4">Boundary Conditions: Initial and Final States</th><th>System Mass</th><th>Pulsation</th><th>Dry Friction Coefficient</th></tr><tr><th><i>T</i></th><th><i>q</i>(0)</th><th>ζ(0)</th><th><i>q</i>(<i>T</i>)</th><th>ζ(<i>T</i>)</th><th><i>M</i></th><th>ω</th><th><i>f</i><sub>c</sub></th></tr></thead><tbody><tr><td>2 s</td><td>0 rad</td><td>0 rad s<sup>-1</sup></td><td>0.8 rad</td><td>0 rad s<sup>-1</sup></td><td>1.384 kg m<sup>2</sup></td><td>2.609 s<sup>-1</sup></td><td>13.683 N m</td></tr></tbody></table>

We use Mathematica's NDSolveValue ODE solver set up to operate with an FEM option on a mono-dimensional, uniform mesh. The fixed step size is chosen to be of $h = 10^{-4}$ s. The integration method is of order 2 and uses 20 digits of precision. Other NDSolveValue options covering finite differences methods, stiffness-specific procedures, and shooting methods were tried but were unsuccessful in solving Equation (11) subject to the boundary values of Table 1. The idea is to seek the largest possible value of $k$ so that the solver still converges to a solution.

Unless the FEM option of NDSolveValue is used, solutions can be obtained only up to $k = 0.5$ s rad$^{-1}$, which does not allow for much of an analysis. Therefore, using the FEM option, the value of $k$ was gradually increased until NDSolveValue failed to return a solution that met the boundary values of Table 1. Note that when simulating stiff systems, step sizes should be chosen so that these are smaller than the smoothing parameter's inverse (see e.g., [43]). For our chosen smooth friction model, the smoothing parameter is
$k$. In our case, the step size is much smaller than the inverse of the smoothing parameter:
$h = 10^{-4} \text{ s} \ll 1/k$.

During the optimally controlled motion simulation of the one-DOF system depicted in Figure 2, the maximum value of $k$ for which NDSolveValue could find a solution was $k = 55 \text{ s rad}^{-1}$, even when further diminishing the step size.

Figure 3 shows the solutions obtained with four different values of $k$. Note that solutions start reaching a limit curve already for $k = 10 \text{ s rad}^{-1}$ because the corresponding optimal position solution is almost the same as the one obtained for $k = 55 \text{ s rad}^{-1}$. Indeed, those two corresponding position curves are virtually indistinguishable in Figure 3a. The difference in the corresponding joint torque solutions is a little more perceptible without being overwhelming (see Figure 3b). The torque solution tends to a limit curve in which we expect sharp corners (see Figure 3b). For practical purposes, the solution for $k = 10 \text{ s rad}^{-1}$ might suffice.

The NDSolveValue solver faces limitations when the hyperbolic tangent function coefficient $k$, of the chosen friction model (3), reaches high values. The specialized literature of friction modeling and compensation often mentions numerical stiffness to justify issues such as this one, i.e., the difficulty of obtaining solutions with ODE solvers [2]. Classical methods to estimate and measure stiffness generally involve the evaluation of eigenvalues, Jacobians, Lipschitz constants, and norms [44,45]. However, in our specific optimal control example, it is easy to see that stiffness increases just as $k$ increases. Consider Figure 4a: it shows the evolution of the maximum slope at each torque solution obtained for growing values of $k$. Figure 3b shows that the maximum slope increases with $k$. This evolution is linear but quite steep: the slope of the fitted model for this evolution is of $104.21 \text{ N m s}^{-2}$. Figure 4b shows a zoom-in around the corners of the torque curves where the maximum slope for the torque curve when $k = 55 \text{ s rad}^{-1}$ is displayed. This maximum slope is of $5725.1 \text{ N m s}^{-1}$, which is very steep. Therefore, we can expect most numerical methods to fail in solving the system (11).

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Optimal motion of the actuated mechanical system subject to dissipative effects modeled by the Coulomb-tanh friction model. Solutions to the system (11) subject to the boundary values of Table 1, obtained with Mathematica's NDSolveValue ODE solver, set up to operate with an FEM. The solver finds a solution for $k \le 55 \text{ s rad}^{-1}$: (a) optimal joint torque solutions for different values of $k$; (b) zoom-in at the corners. Note that solutions change very little for $k > 10 \text{ s rad}^{-1}$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Optimal motion of the actuated mechanical system subject to dissipative effects modeled by the Coulomb–tanh friction model. Analysis of the optimal torque solution evolution with a growing value of *k*: (a) maximum slope $\dot{\zeta}_{max}$ of torque curves’ evolution as *k* increases; (b) zoom-in around the sharp corners of torque curves for different values of *k*. The maximum slope evolution is linear with a steep rate of change. The maximum slope tends to infinity as *k* reaches very high values. Convergence toward the classical Coulomb’s model is revealed.</figcaption></figure>

In an extreme case where *k* could reach a much higher value, the torque curve in Figure 3b would showcase very sharp corners with a maximum slope that tends to infinity. However, considering the maximum torque curve slope, we can conclude that such a value for *k* already approximates Coulomb’s friction very well in the current situation. Of course, this needs to be verified experimentally. Therefore, the experimental validation of this hypothesis will be presented in the next section.

Now, consider Figure 5. This figure shows a short analysis of the friction force evolution during optimal motion simulations as the value of *k* grows. Friction force stabilizes around a limit curve for $k \ge 10 \text{ s rad}^{-1}$. Again, as shown by Figure 5a, the only notable difference between the curve for $k = 10 \text{ s rad}^{-1}$ and the curve for $k = 55 \text{ s rad}^{-1}$ is around the corners, which are already acute for the latter case. Note that a plateau is present in the friction force curves. In this plateau, the force corresponds to Coulomb’s friction coefficient $f_c = 13.683 \text{ N m}$. This behavior is typical of Coulomb’s friction force law, so we can conclude that the optimal control procedure gives a good representation of the Coulomb friction force even for small values of *k*. Additionally, consider Figure 5b, which shows the evolution of the ratio between the maximum friction force $f_{max}$ and the maximum torque $\zeta_{max}$ that was recorded during optimal motion simulations for growing values of *k*. It is interesting to notice that as *k* grows, $f_{max}$ stabilizes around a value of exactly half the maximum torque for each corresponding motion.

Let us remark that other popular solvers could be used. For example, Matlab’s bvp4c and bvp5c solvers can be invoked to solve the two-point boundary value problem of equations (11) subject to the boundary values of Table 1. However, even with a step size of $h = 10^{-4} \text{ s}$ (as in our example), usable solutions could not be obtained for $k \ge 5$.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Optimal motion of the actuated mechanical system subject to dissipative effects modeled by the Coulomb–tanh friction force law. Analysis of the friction force evolution during optimal motion with a growing value of *k*: (a) friction torque *f* as *k* increases for the optimal motion; (b) ratio between maximum friction force *f*<sub>max</sub> and maximum optimal torque ζ<sub>max</sub> for growing values of *k*. Note that for *k* ≥ 10, friction force stabilizes at the value of the Coulomb friction coefficient *f*<sub>c</sub> = 13.683 N m. Subfigure (b) shows the evolution of the ratio between the maximum friction force *f*<sub>max</sub> and maximum optimal torque ζ<sub>max</sub> for growing values of *k*. The maximum friction force is about half the maximum torque recorded during optimal motion. Notice that friction torques converge toward the traditional Coulomb’s law.</figcaption></figure>

3.2. *Trajectory Duration Influence*

The previous subsection analysis was conducted with a fixed trajectory duration of
T = 2 s. However, it was observed that T does have an impact not only on the solver
performance (equations become stiffer as T increases) but also on motion characteristics (as
observed in [9,13]).

In particular, the analyzed motion where T = 2 s (see previous subsection) showcased
only one velocity sign change in the absence of friction (see Figure 3a). One could then
wonder what would happen if there were more velocity sign changes when omitting
friction in the motion equations. Would this also affect the optimal motion where friction
has been directly taken into account? We have conducted a series of numerical tests for
four different trajectory durations. Let us remark that the results and conclusions presented
in the previous subsection still apply when the trajectory duration changes. Therefore, we
will focus on other aspects from here on.

Consider Table 2 which lists some results of our numerical trials. The first observation is that even when there are several changes in the velocity sign during optimal motion without friction, there are no sign changes in the case where friction is considered, regardless of trajectory duration. The second observation is that for T > 3 s, the optimal control problem becomes considerably stiffer (numerically speaking) and so the *k* coefficient cannot reach a very high value. In that regard, a coefficient of *k* = 14 s rad⁻¹ is still representative of Coulomb friction, as was shown in Section 3.1, but is lower than for the other cases. Interestingly, *k* reached a lower value for T = 1 s or T = 2 s than for T = 3 s. This is counterintuitive and surprising because the analysis carried out in [9] showed that stiffness tends to increase with higher values of T.

Let us now focus on Figures 6 and 7 where optimal solutions for position and torque
are displayed for four different trajectory durations. It is striking that optimal solutions, in
the case where friction is considered, are practically the same for every value of T. Figure 6
shows the optimal positions for the case without friction (when k = 0 s rad⁻¹). Notice
that there are various velocity sign changes for T > 1 s. Intuitively, there would at least
be one velocity sign change somewhere in the optimal position for T = 4 s but this does
not happen. Instead, the system performs only one swing to get to the final position for every duration. The torque curves (Figure 7) confirm that the optimal solution is the same for every duration (with the corresponding amount of zero values before motion starts). Maximum torque reaches the same value in every case. It is twice the value of the Coulomb friction coefficient $f_c$. This observation may not be reproduced when the system has more degrees of freedom but at least we can observe that motion is effectively damped.

<table><caption>**Table 2.** Optimal motion characteristics for four different values of trajectory duration *T*. An optimal velocity solution may showcase multiple sign changes when friction is not taken into account. When considering friction within the model, no velocity sign changes are recorded. The highest value of *k* (*k*<sub>max</sub>) is attained for *T* = 3 s.</caption><thead><tr><th rowspan="2">T</th><th colspan="2">Optimal Velocity Sign Changes</th><th rowspan="2">k<sub>max</sub></th></tr><tr><th>No Friction (k = 0)</th><th>With Friction (k = k<sub>max</sub>)</th></tr></thead><tbody><tr><td>1 s</td><td>0</td><td>0</td><td>48 s rad<sup>-1</sup></td></tr><tr><td>2 s</td><td>1</td><td>0</td><td>55 s rad<sup>-1</sup></td></tr><tr><td>3 s</td><td>2</td><td>0</td><td>70 s rad<sup>-1</sup></td></tr><tr><td>4 s</td><td>3</td><td>0</td><td>14 s rad<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Optimal positions for four trajectories with growing duration *T*, without and with friction. Curves calculated without friction are shown along those calculated by taking friction into account for the maximum value of *k* attained in each case (*k*<sub>max</sub>). Notice how optimal positions that consider friction tend to the same limit curves in each case, regardless of the trajectory duration.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7. Optimal torques for four trajectories with growing duration T, without and with friction. Each curve shows: optimal torques without considering friction ($\xi_{k=0}$); optimal torques considering friction ($\xi_{k=k_{max}}$); optimal torques without considering friction with added friction compensation ($\xi_{comp}$). Notice how solutions that consider friction ($\xi_{k=k_{max}}$) tend to the same limit curves in each case, regardless of the trajectory duration.</figcaption></figure>

### 3.3. Optimal Control Performance

Optimal control performance is now assessed. Let us recall that the goal of including friction forces into the system model is for the optimal control procedure to take them directly into account so that friction effects are inherently minimized. We aim to show the difference between our optimal control that considers friction and a more traditional approach where friction is compensated at the feedback control level. The latter procedure was presented in [8,13] and did not consider friction in the system model. This case would correspond to the above solutions for $k = 0$ s rad$^{-1}$. Friction forces were then added at the feedback control level. This process can be simulated by creating a friction-compensated torque variable

$$ \xi_{\text{comp}} = \xi_{k=0} + f_c \tanh(k_{\text{max}}\xi_{k=0}) \quad (12) $$

where we shall recall that $\xi_{k=0}$ is the optimal torque in the absence of friction and $\xi_{comp}$ is the friction-compensated torque; $f_c$ is the Coulomb friction coefficient; $k_{max}$ is the highest possible value of $k$; and $\xi_{k=0}$ is the optimal velocity in the absence of friction. Other variables may arise in feedback control, depending on the chosen controller. However, (12) gives a good approximation of what to expect.

Figure 7 shows the $\xi_{comp}$ variable evolution for four different trajectory durations T. The difference with our proposed solution ($\xi_{k=k_{max}}$) is quite noticeable. Obviously, in the added friction compensation case, the curve is an exaggeration of the optimal solution that does not consider friction ($\xi_{k=0}$) because we are simply adding a force that has the value of the Coulomb friction coefficient $f_c$. Consequently, the compensated solution suffers the
same number of transitions between positive and negative values as the optimal motion that does not consider friction (see the second column of Table 2), resulting in an oscillatory torque evolution.

To assess the performance increase by our proposed method, we analyze the behavior of the cost functional (6) that is regarded as the performance index. We compare the value of $J(\xi_{\text{comp}})$ against the values of $J(\xi_k)$ for different values of $k$. Table 3 shows the evaluation of $J(\xi_{\text{comp}})$ for the case where friction compensation is performed on the optimal solution that does not consider friction, in the second column. $J(\xi_k)$ was evaluated for other values of $k$, in the optimal case that takes friction into account, after the third column of Table 3.

Let us focus on the second and third columns of Table 3: $J(\xi_{\text{comp}})$ grows as $T$ increases, while $J(\xi_{k=0})$ drops its value. This was already observed in [13]. Most importantly, the performance increase becomes evident with our proposal: while $J(\xi_{\text{comp}})$ grows as $T$ increases, $J(\xi_{k=k_{\text{max}}})$ seems to settle around the value of 78 N m s⁻¹, regardless of the trajectory duration. This property shall be experimentally verified in the next section. Finally, the last column of Table 3 shows the performance increase that should be obtained with our proposal, which grows with the trajectory duration.

<table><caption>**Table 3.** Performance indicators for the one-DOF system. Cost functional $J(\xi)$ evaluation for: optimal torques without considering friction with added friction compensation ($\xi_{\text{comp}}$); optimal torques with various values of the $k$ coefficient $\xi_{k=+}$. The cost function is regarded as the performance index of our optimal control method: lower values indicate better performance. Optimal torques $\xi_{k=0}$ require friction compensation.</caption><thead><tr><th rowspan="2">T (s)</th><th colspan="6">J(ξ) = ∫<sub>0</sub><sup>T</sup> M<sup>ij</sup>ξ<sub>i</sub>ξ<sub>j</sub>, Measured in N m s<sup>-1</sup></th><th rowspan="2">Performance Increase</th></tr><tr><th>J(ξ<sub>comp</sub>)</th><th>J(ξ<sub>k=0</sub>)</th><th>J(ξ<sub>k=14</sub>)</th><th>J(ξ<sub>k=48</sub>)</th><th>J(ξ<sub>k=55</sub>)</th><th>J(ξ<sub>k=70</sub>)</th></tr></thead><tbody><tr><td>1</td><td>104.24</td><td>4.99</td><td>77.27</td><td>78.24</td><td>–</td><td>–</td><td>24.94%</td></tr><tr><td>2</td><td>168.77</td><td>2.77</td><td>77.27</td><td>78.24</td><td>78.29</td><td>–</td><td>53.64%</td></tr><tr><td>3</td><td>235.33</td><td>1.99</td><td>77.27</td><td>78.24</td><td>–</td><td>78.36</td><td>66.75%</td></tr><tr><td>4</td><td>276.25</td><td>1.51</td><td>77.27</td><td>–</td><td>–</td><td>–</td><td>71.68%</td></tr></tbody></table>

### 3.4. Optimal Control Simulation of a Two-DOF System Considering Friction

This section illustrates that the proposed procedure can be used on multi-DOF systems. The following simulations were achieved on a two-DOF mechanical system that can be regarded as an actuated double pendulum. The parameters for this system are reported in [13] and recalled in Figure 8, which depicts the system. For the following simulations, the parameters are indicated by Table 4. These are used in Equations (4) or (5) to describe the system dynamics. The optimal control procedure exposed in Section 2.4 is then followed for this simulation. In this example, the resulting first-order system (9) is subject to the boundary values from Table 4.

<table><caption>**Table 4.** Simulation parameters for the optimal motion simulation of a two-DOF system subject to joint-level friction.</caption><thead><tr><th>Final Time</th><th colspan="5">Boundary Conditions: Initial and Final States</th><th>Dry Friction Coefficients</th></tr><tr><th>T</th><th>(q<sup>1</sup>(0), q<sup>2</sup>(0))</th><th>(ζ<sup>1</sup>(0), ζ<sup>2</sup>(0))</th><th>(q<sup>1</sup>(T), q<sup>2</sup>(T))</th><th>(ζ<sup>1</sup>(T), ζ<sup>2</sup>(T))</th><th>(f<sub>11</sub>, f<sub>12</sub>)</th><th></th></tr></thead><tbody><tr><td>1.5 s</td><td>(0, 0) rad</td><td>(0, 0) rad s<sup>-1</sup></td><td>(0.8, 1) rad</td><td>(0, 0) rad s<sup>-1</sup></td><td>(13.683, 2.018) N m</td><td></td></tr></tbody></table>

Figure 9 presents the obtained simulation results for the optimal control of the two-DOF system depicted by Figure 8. Naturally, numerical stiffness increases in this configuration. The maximum value that was attained for $k = 20$ s rad⁻¹ for this two-DOF
configuration. However, the system seems to be reaching limit curves at $k = 10 \text{ s rad}^{-1}$. This confirms the tendency also observed for the one-DOF case (see Figures 3 and 5a).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Representation of the two-DOF mechanical system. Motion is generated around the centers of rotation ($O_1, O_2$) by actuators delivering the respective forcing torques ($u_1, u_2$). The generalized coordinates ($q^1, q^2$) are the angular positions. Parameters ($m_1, m_2, I_1, I_2$) = (23.902 kg, 3.880 kg, 1.266 kg m$^2$, 0.093 kg m$^2$) are the respective masses and inertias of each link. ($C_1, C_2$) are the positions of the centers of mass. Vector $g$ indicates the gravity action direction.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Optimal motion of a two-DOF mechanical system subject to dissipative effects modeled by the Coulomb–tanh friction model. Solutions to the system (9) subject to the boundary conditions of Table 4, obtained with Mathematica's NDSolveValue ODE solver, set up to operate with an FEM. Solutions are shown for growing values of $k$ up to $k \le 20 \text{ s rad}^{-1}$. Note that solutions change very little for $k > 10 \text{ s rad}^{-1}$.</figcaption></figure>
Finally, let us consider Figure 10. It shows the optimal control solution on the idealized case where friction is absent ($\tilde{\xi}_{k=0}$). Then, $\tilde{\xi}_{k=0}^{\text{comp}}$ simulates the traditional case where friction is managed at feedback-level, by simple addition, as in Equation (17). In this added friction compensation case, curves $\tilde{\xi}_{k=0}^{\text{comp}}$ become an exaggeration of $\tilde{\xi}_{k=0}$ because a force is simply added afterward. Conversely, our proposed solution is represented by the $\tilde{\xi}_{k=20}$ curves. Being smoother, these are easier on the actuators of the system. Most importantly, our proposed solution shall lead to a better performance, as compared to the traditional feedback-level friction compensation case. This is shown by Table 5 where the performance increase is shown for this simulation. Note that such performance further increases for longer trajectories, as shown by Table 3.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Optimal torques for the trajectory of the two-DOF system. The curve $\tilde{\xi}_{k=0}$ is the optimal torque when friction is not considered. The curve $\tilde{\xi}_{k=0}^{\text{comp}}$ represents $\tilde{\xi}_{k=0}$ with added friction compensation (simulating feedback control-level friction compensation). The curve $\tilde{\xi}_{k=20}$ is the proposed optimal solution when friction is directly managed by the optimal control procedure. The proposed solution $\tilde{\xi}_{k=20}$ is easier on the actuators of the system.</figcaption></figure>

<table><caption>**Table 5.** Performance indicators for the two-DOF system. Cost functional $J(\xi)$ evaluation for: optimal torques with feedback-level friction compensation ($\tilde{\xi}_{k=0}^{\text{comp}}$); optimal torques with two values of the k coefficient $\tilde{\xi}_{k=1}$.. The cost function is regarded as the performance index of our optimal control method: lower values indicate better performance. Optimal torques $\tilde{\xi}_{k=0}$ require friction compensation.</caption><thead><tr><th rowspan="2">T (s)</th><th colspan="4">J(ξ) = ∫<sub>0</sub><sup>T</sup> Mijξ<sub>i</sub>ξ<sub>j</sub> (N m s<sup>-1</sup>)</th><th rowspan="2">Performance Increase</th></tr><tr><th>J(ξ<sub>k=0</sub><sup>comp</sup>)</th><th>J(ξ<sub>k=0</sub>)</th><th>J(ξ<sub>k=11</sub>)</th><th>J(ξ<sub>k=20</sub>)</th></tr></thead><tbody><tr><td>1.5</td><td>312.09</td><td>39.88</td><td></td><td>293.55</td><td>5.94%</td></tr><tr><td>3</td><td>391.90</td><td>17.26</td><td>348.65</td><td>–</td><td>11.04%</td></tr></tbody></table>

The simulations in this section show that our proposed optimal control procedure is compatible with multi-DOF systems and also result in increased performance.

## 4. Experimental Validation

We have conducted a series of experiments on a one-DOF mechanical system with a rotating actuator. The goal is to validate the simulation results presented in the previous section. We begin by briefly describing our experimental setup and the parameter identification procedure. The obtained results are then presented.
## 4.1. Experimental Setup

To validate our simulation results, we perform some experiments on an actuated single DOF mechanical system. A direct current (DC) motor is directly mounted on the axis so that the system operates as a direct drive joint. It is a McMillan Electric Company C3350B3006 brushed DC motor to which a metal disc holding a pendulum arm is directly attached at its center of rotation (see Figure 11). A 4096 pulses per round optical encoder was installed on the motor axis at the rear part of the motor. The result is a resolution of $1.5 \times 10^{-3}$ rad ($8.79 \times 10^{-2}\circ$). Velocity is estimated by Euler backward differentiation of the pendulum position. Then, the average of 3 consecutive velocity values is taken. A round, homogeneous, and rigid ABS plastic part was placed on the rim of the metallic disc. This slightly increases joint-level friction.

<figure><img src="image_14.png" /><figcaption>Figure 11. One-DOF experimental system.</figcaption></figure>

System control is achieved through Matlab Simulink (version 7.0) and QuaRC software (version 1.1), operated in real-time by a Windows desktop computer linked to a Quanser Q4 data acquisition card. QuaRC ensures communication between Matlab's Simulink environment and the Quanser Q4 device. The real-time system was configured with a sampling period of $1 \times 10^{-3}$ s. The encoder is connected to the Quanser Q4 card. Data are then sent to the UPM-2405 power supply that powers up the DC motor. An in-house current meter was designed to drive the motor in torque mode.

## 4.2. System Parameters

Our procedure to estimate the system parameters, which shall be used in Equation (10), requires that the torque mode control operates accurately. In this sense, we need to identify some electrical parameters first. It is known that the electric model of a DC motor, neglecting the armature inductance, is given by (see e.g., [46])

$$v = k_b \dot{q} + R_a i_a, \qquad (13)$$

where $v$ is the armature voltage, $R_a$ its resistance, $i_a$ its current, $k_b$ its counter-electromotive force constant, and $\dot{q}$ its angular velocity. It is also known that the motor outputs a torque $u$, which is proportional to $i_a$,

$$u = k_a i_a, \qquad (14)$$

where $k_a$ is the motor torque constant. If we consider that losses are minimal during motor operation, then $k_a = k_b$. Therefore, if we divide (13) by $i_a$ on both sides,

$$\frac{v}{i_a} = k_a \frac{\dot{q}}{i_a} + R_a. \qquad (15)$$
This means that we can estimate $k_a$ and $R_a$ by measuring $v, i_a$, and $\dot{q}$ since Equation (15) can be regarded as a linear equation where $v/i_a$ is the dependent variable and $\dot{q}/i_a$ is the independent variable. This experiment gave the values $k_a = 0.1762 \text{ N m A}^{-1}$ and $R_a = 5.54 \text{ }\Omega$. Per (14), it suffices to know $k_a$ and supply an adequate $i_a$ to generate the required torque. Therefore, to achieve the torque mode operation of the DC motor, a proportional-integral (PI) controller is designed to regulate the current, using gain values of $k_p = 5$ and $k_i = 5000$.

Next, the identification tool included in Simulink was used for the parametric identification of the actuated mechanical system (10). The “Nonlinear Least Squares” method was selected within the “Reflective Trust Region” algorithm. An open-loop experiment was then conducted by using the input torque signal

$$u(t) = 0.3(1 - e^{-1.8t}) + 0.1 \sin(26t + 0.08) + 0.1 \sin(12t + 0.34). \quad (16)$$

A hyperbolic tangent slope parameter of $k = 100 \text{ s rad}^{-1}$ was used during the parameter identification. The identified values are indicated in Table 6. Figure 12 shows the simulated joint position time evolution compared to the experimental joint position result. Both curves are acceptably close to one another.

<table><caption>**Table 6.** Identified parameters of the experimental system of Figure 11.</caption><thead><tr><th>System Mass</th><th>Pulsation</th><th>Dry Friction Coefficient</th></tr><tr><th><i>M</i></th><th><i>ω</i></th><th><i>f</i><sub>c</sub></th></tr></thead><tbody><tr><td>1.367 × 10<sup>-2</sup> kg m<sup>2</sup></td><td>5.591 × 10<sup>-1</sup> s<sup>-1</sup></td><td>1.075 × 10<sup>-1</sup> N m</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of position time evolution: simulation versus experimental result after identifying system parameters. Input signal (16) and the identified system parameters of Table 6 were used on the motion Equation (10) for this experiment. The measured position is acceptably close to the simulated one.</figcaption></figure>

### 4.3. *Experimental Results Analysis*

For this experimental exercise, we followed the optimal control procedure exposed in Section 3.1. That is, we solved the nonlinear system (11) subject to the boundary conditions of Table 4, with a fixed value of $T = 2$ s, but using the system parameters given in Table 6. The maximum value of $k$ that could be attained using these parameters was of $k = 56 \text{ s rad}^{-1}$ with which we obtained the optimal solutions $q_{k=56}$ and $\bar{q}_{k=56}$. These shall be used as references for the main experiments.

Before proceeding, let us clarify two nomenclature details.

1. A tilde over a variable indicates an experimental open-loop control result. For example, $\bar{q}_{k=56}$ is the open-loop position resulting from taking $q_{k=56}$ as a reference.

2. A bar over a variable indicates an experimental closed-loop control result. For example, $\bar{q}_{k=56}$ is the closed-loop position when taking $q_{k=56}$ as a reference for feedback.

We will begin by presenting the results obtained when controlling the system in an open loop. Then, we will show the results obtained in a closed loop with friction
compensation at the feedback level. Finally, we will show the results obtained when the
optimal control procedure managed friction.

4.3.1. Open-Loop Control

This first open-loop experiment is important to verify the system response when controlled in an open loop. If the system reaches the goal position within acceptable bounds, this means that the selected model is adequate because there is no controller feedback in this experiment.

Equation (11) is solved subject to the boundary conditions of Table 1 on a $T = 2$ s duration, fixing the system parameters as in Table 6. This is carried out for several values of the $k$ slope parameter. Then, the corresponding torques $u_k^-=$. are directly fed to the system without feedback. Figure 13 shows the tendency of open-loop control positions $q_k^-=$. as $k$ grows, compared with the best optimal reference attained at $k = 56$. Notice that without added friction compensation at the feedback level, the system cannot move and remains at the initial position ($\tilde{q}_{k=0}$ curve of Figure 13). Reaching the end goal position is impossible in this manner. Therefore, when $k = 0$, friction must be compensated at the feedback level in a closed loop to reach the end goal position.

However, notice how for $k = 1$ there is a substantial improvement even if the goal end position is still unattained. Then, a limit curve seems to be attained for $k = 10$. Curves cease to significantly improve for $k \ge 10$. For this reason, Figure 13 also shows a zoom-in around the end goal position so that one can better appreciate what is happening: as $k$ grows, open-loop position curves seem to settle around some region near the end goal position without attaining it. Nevertheless, this experiment proves that a very high value of $k$ is not required to properly approach Coulomb's friction force because good results are already obtained for $k = 10 \text{ s rad}^{-1}$ with little deviation from the optimal reference trajectory.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Open-loop positions $\tilde{q}_{k=\infty}$. for various values of $k$ compared against the best optimal reference $q_{k=56}$. Notice that open-loop positions are very close to the optimal reference for $k \ge 10$. Zooming near the end of the trajectory reveals that the attained end positions converge around a region close to the desired end goal position for $k \ge 10$.</figcaption></figure>

4.3.2. Classical Closed-Loop Control with Feedback-Level Friction Compensation

This experiment shows a more classical procedure to deal with friction at the joint level. In this experiment, friction is not considered in Equation (11), meaning that $k = 0 \text{ s rad}^{-1}$. The idea is to show that the optimal control reference is heavily changed when dealing with friction in this manner.
The proportional-derivative (PD) plus feedforward controller found in [46] was used in our experiments. In this particular case, we added friction compensation at the feedback level, so that

$$ \bar{\xi}_{k=0}^{\text{comp}} = k_p(q_{k=0} - q) + k_v(\dot{q}_{k=0} - \dot{q}) + \xi_{k=0} + f_c \tanh(100\dot{q}), \quad (17) $$

because we have fixed the hyperbolic tangent slope parameter at $k = 100 \text{ s rad}^{-1}$ for this experiment. In the above control law, $q$ and $\dot{q}$ are, respectively, the measured position and velocity of the system at each instant.

The closed-loop curves in Figure 14 were obtained by using $k_p = 6$ and $k_v = 0.005$. Notice how, for practical purposes, the position tracking goal is achieved. The maximum recorded positioning error was of $3.68 \times 10^{-2} \text{ rad}$ (2.11°) at $t = 1.77 \text{ s}$. This means that the controller is performing acceptably well as can be seen in Figure 14a where $\bar{q}_{k=0}^{\text{comp}}$ tracks $q_{k=0}$ with minor deviations during the trajectory. Now, consider Figure 14b where the reference optimal torque $\bar{\xi}_{k=0}^{\text{comp}}$ was modified according to Equation (17) in order to produce $\bar{\xi}_{k=0}^{\text{comp}}$. Consequently, some spikes were recorded along the trajectory. This is the effect of adding the Coulomb-tanh friction model to the control law in (17), which occurs when velocity is very close to zero. Therefore, $\bar{\xi}_{k=0}^{\text{comp}}$ is very different from the optimal reference $\bar{\xi}_{k=0}$, reaching maximum values that are more than five times higher than the maximum reference ones. This means that the system response is far from the optimal reference.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Comparison between optimal reference variables and experimental results for the case where friction is compensated at feedback level: (a) position and (b) torques.</figcaption></figure>

The results from this experiment motivate our procedure. In the following experiment, we will use the optimal reference for which $k \neq 0$ in the Coulomb-tanh friction model. The results will then be compared with those obtained in this experiment.

### 4.3.3. Closed-Loop Control with Friction Managed by the Optimal Control Procedure

Our final experiments use the same previous conditions and parameters. However, friction is now directly taken into account by the optimal control procedure. Therefore, in (10), parameter $k$ takes the highest attained value for this system: $k = 56 \text{ s rad}^{-1}$.

In this scenario, the controller no longer compensates friction at the feedback level, so now,

$$ \bar{\xi}_{k=k_{\text{max}}} = k_p(q_{k=k_{\text{max}}} - q) + k_v(\dot{q}_{k=k_{\text{max}}} - \dot{q}) + \xi_{k=k_{\text{max}}}, \quad (18) $$
where $q$ and $\dot{q}$ correspond to the measured system position and velocity at each instant. Note that for this experiment, the same gain values as in the previous experiment were used ($k_p = 6$ and $k_v = 0.005$).

Figure 15 shows the curves obtained for this experiment: system position and required torque for the trajectory. Notice how position tracking was properly achieved using the closed-loop control procedure (18) since curves $q_{k=56}$ and $\bar{q}_{k=56}$ are indistinguishable (see Figure 15a). Indeed, the results are now better than those obtained in the scenario where friction is compensated at feedback level (see Section 4.3.2): maximum error occurred at $t = 1.818$ s with a deviation of $1.67 \times 10^{-2}$ rad (0.96°). Figure 15a clearly shows when the controller compensates deviations by either increasing or decreasing the input torque value with respect to the optimal reference $\bar{\zeta}_{k=56}$. Notice how by the end of the experiment, around $t = 1.9$ s, the closed-loop torque $\bar{\zeta}_{k=56}$ settles at the desired values. These small deviations are the consequence of unmodeled dynamics and uncertainties. We hypothesize that a more complete friction model would further reduce these deviations.

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Comparison between optimal reference variables and experimental results for the case where friction is directly taken into account by the optimal control procedure: (a) positions and (b) torques.</figcaption></figure>

Recall that the maximum torque resulted in being twice as big as the friction force (which is equal to the friction coefficient $f_c$, see Figure 5b) in Section 3.1's simulations. During these experiment simulations, however, we recorded that the ratio $\frac{\max(\bar{\zeta}_{k=56})}{f_c} \approx 3$ and in the closed-loop experiment, $\frac{\max(\bar{\zeta}_{k=56})}{f_c} \approx 3.5$. Interestingly, some relationship might exist between these values (see Section 3.1), but we cannot conclude on a specific tendency here.

We must also remark that the open-loop position $\bar{q}_{k=56}$ is very close to the reference $q_{k=56}$ (see Figure 15a). This was observed in Figure 13. More importantly, open-loop control quality is radically increased by our approach.

### 4.3.4. Performance Increase with the Proposed Optimal Control Procedure

Consider Figure 16 which shows the torque with feedback-compensated friction $\bar{\zeta}_{k=0}^{\text{comp}}$ and the closed-loop torque $\bar{\zeta}_{k=56}$. Notice how these curves generally behave as predicted by simulations in Section 3.3 (with added signal noise). Some differences with the simulations of Section 3.3 reveal unmodeled uncertainties and friction effects. Indeed, we focused on approximating Coulomb friction. Nevertheless, the resulting closed-loop torques were globally anticipated through simulations. We can observe that $\bar{\zeta}_{k=0}^{\text{comp}}$ and $\bar{\zeta}_{k=56}$ reach similar maximum values. However, the system exerts less force across the trajectory since only one swing is required to reach the end goal position.
Table 7 shows the performance increase due to our optimal control procedure: the performance index value is lower. The ratio between $J(\bar{\xi}_{k=0}^{\text{comp}})$ and $J(\bar{\xi}_{k=k_{\text{max}}})$ indicates the performance increase. We must remark that the moment of inertia is smaller in our experimental platform than in our simulations (see Section 3). This factor explains why the figures in Table 7 are lower than those in Table 3. Additionally, Table 3 shows that $J(\bar{\xi}_{k=0}^{\text{comp}})$ grows as $T$ increases; Table 7 further verifies this result experimentally. Therefore, as $T$ grows, the performance gap increases between both situations just as we predicted in Section 3. Finally, we must retain the following: $J(\bar{\xi}_{k=0}^{\text{comp}}) \gg J(\bar{\xi}_{k=0})$, whereas $J(\bar{\xi}_{k=k_{\text{max}}}) \approx J(\bar{\xi}_{k=k_{\text{max}}})$, but most importantly, $J(\bar{\xi}_{k=k_{\text{max}}}) < J(\bar{\xi}_{k=0}^{\text{comp}})$.

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Closed-loop torques for two scenarios: friction compensation at feedback level $(\bar{\xi}_{k=0}^{\text{comp}})$, and friction taken into account by the optimal control procedure $(\bar{\xi}_{k=56})$. The solution $\bar{\xi}_{k=56}$ is the better one because it is much closer to its optimal reference and is more cost-effective.</figcaption></figure>

<table><caption>**Table 7.** Cost functional $J(\xi)$ evaluation for two scenarios: feedback level friction compensation $J(\bar{\xi}_{k=0}^{\text{comp}})$, and optimal control procedure that takes friction into account $J(\bar{\xi}_{k=k_{\text{max}}})$ compared with their optimal reference counterparts. Lower values indicate better performance.</caption><thead><tr><th rowspan="2">T (s)</th><th rowspan="2">k<sub>max</sub></th><th colspan="3">J(ξ) = ∫<sub>0</sub><sup>T</sup> M<sup>ij</sup>ξ<sub>i</sub>ξ<sub>j</sub> (N m s<sup>-1</sup>)</th><th rowspan="2">Performance Increase</th></tr><tr><th>J(ξ̃<sub>k=0</sub><sup>comp</sup>)</th><th>J(ξ̃<sub>k=k<sub>max</sub></sub>)</th><th>J(ξ̃<sub>k=k<sub>max</sub></sub>)</th></tr></thead><tbody><tr><td>1.5</td><td>50</td><td><b>1.7189</b></td><td>1.1787</td><td><b>1.0960</b></td><td>36.24%</td></tr><tr><td>2</td><td>56</td><td><b>1.9779</b></td><td>1.1789</td><td><b>1.1585</b></td><td>41.42%</td></tr><tr><td>2.5</td><td>26</td><td><b>1.8706</b></td><td>1.1736</td><td><b>1.2076</b></td><td>35.44%</td></tr><tr><td>3</td><td>23</td><td><b>2.4580</b></td><td>1.1717</td><td><b>1.2218</b></td><td>50.29%</td></tr></tbody></table>

## 5. Discussion

The selected smooth dry friction law (3) is an approximation of Coulomb's friction force law. Therefore, some limitations are observed. For instance, the Coulomb-tanh model supposes motion at very low velocities, so stiction cannot be represented in this manner. Effects such as Stribeck's friction, or the break-away force, are non-existent in our numerical trials. Computationally speaking, the selected dry friction approximation does increase numerical stiffness when the slope coefficient grows (see Figure 4). As a result, most numerical methods fail to integrate the optimal control problem equations in high-stiffness situations. Note that acceptable results are obtained for a relatively small value of the slope coefficient though (see Figures 3, 5 and 9). However, our contribution is an exploration into how to deal with friction within the indirect optimal control method proposed in [8,9,42], which could be useful for the control of dynamical systems. A different friction law can be
selected without changing the optimal control methodology, but an appropriate numerical
strategy should be required.

Regarding the selected geometric optimal control procedure, it has been compared
with more classical formulations in [8,9], where some numerical advantages are described.
In this first step, we focused on joint-level Coulomb's friction because it is the most pre-
dominant dissipative effect during motion. Indeed, most models generalize Coulomb's
friction with added effects. These at least represent Coulomb's friction when simplified to
the extreme. This is the case for the widely used Dahl and LuGre models for example [2,5].

A straightforward enhancement of the selected friction model consists of adding a viscous friction term such as $f_v \dot{q}$, where $f_v$ denotes the viscous friction coefficient. However, our focus was on dry friction to assess if the selected Coulomb–tanh friction law is actually usable at the optimal planning level. Hopefully, our simulations (Section 3) and experiments (Section 4) convince the reader that Coulomb’s friction can be successfully modeled using the approximate Coulomb–tanh friction law. Note that even if the friction model was considerably enhanced, a precise parameter identification method is still required to increase performance on a real system. It has been pointed out that joint-level friction identification is affected by time [6]. Therefore, the identified parameters vary from trial to trial, even if every other factor remains unchanged. Other ambient conditions have also been reported to negatively affect the parameter identification procedure. Indeed, it is almost impossible to obtain the same friction characteristics in two separate experiments because factors such as temperature and lubricant distribution within the actuator and joint have an (often non-negligible) impact on friction conditions [47]. These problems have been identified a long time ago but remain open. In our case, these issues were minimized by conducting our experiments immediately after identifying the system parameters.

Our experimental observations prove that the simple addition of joint-level Coulomb
friction force at the modeling stage of an optimal control procedure significantly enhances
performance. Our focus is on the proper manner of including friction forces in motion
equations. Their addition led to interesting numerical challenges. For instance, traditional
methods included in most commercial solvers fail to deal with this problem considering
realistic parameters and conditions.

A FEM approach did provide acceptable results. In future developments, we may
adapt symplectic integrators to treat our frictional optimal control problem. These inte-
grators usually give good results for the time integration of dynamical systems in terms
of accuracy, long-term convergence rate preservation, and stability (see [48,49] for exam-
ple). Therefore, one perspective is to formulate our optimal control procedure as a fully
discrete method based on symplectic variational integrators that the authors are currently
developing in this sense [50,51]. On a different tone, to enhance the modeling and simula-
tion approach, subdifferentials [52] or differential inclusions [14,43] may also need to be
considered to deal with discontinuous friction models.

**6. Conclusions**

The proposed optimal control procedure manages joint-level friction by including it in
the system motion equations through Rayleigh’s dissipation function. The procedure is
suitable for mechanical systems with revolute joints. The covariance of the system motion
equations was preserved. The optimality condition is obtained by applying Pontryagin’s
Maximum Principle. The result is that one adjoint parameter is equal to the control variable.
This optimality condition prevails with the added joint-level dissipative effects. It confers
practical and physical interpretation to the adjoint variable.

The smooth Coulomb–tanh friction model was selected to approximate the classical
Coulomb’s friction law. Simulations were conducted on an actuated mechanical system
affected by joint-level friction (see Figure 2).

Simulations showed the following:

1. Optimal solutions quickly approach limit curves for $k \ge 10 \text{ s rad}^{-1}$ (see Figure 3).
2. Stiffness of the optimal control Equation (11) radically increases as $k$ grows (see Figure 4).

3. Coulomb's friction force is well-represented for values of $k \ge 10 \text{ s rad}^{-1}$ (see Figures 5 and 9c,f).

4. Control performance is increased by our approach when compared to the case where friction is compensated at feedback level (see Tables 3 and 5, and Figure 7).

Experiments were then conducted on an actuated pendulum specifically built for this work (see Figure 11). System parameters were identified following the procedure detailed in Section 4.2. Analysis of the experimental results showed that:

1. Open-loop control quality is radically increased with our optimal control procedure that manages friction (see Figure 13).

2. Feedback level friction compensation results in a closed-loop torque $\tilde{\xi}_{k=0}^{\text{comp}}$ that is heavily modified with respect to the optimal reference $\xi_{k=0}$ (see Figures 7 and 14).

3. Conversely, closed-loop system torque $\tilde{\xi}_{k=k_{\text{max}}}$ is very close to the optimal reference $\xi_{k=k_{\text{max}}}$ when joint-level friction is directly taken into account by the optimal control procedure (see Figure 15).

4. Closed-loop control quality is also increased with our proposal (see Sections 4.3.2 and 4.3.3).

5. Optimal control performance increases by taking joint-level dry friction into account when compared to a feedback level friction compensation scenario (see Figure 16 and Table 7).

Simulations and experiments show that incorporating frictional forces into the modeling phase of the optimal control procedure significantly enhances performance compared to a scenario where friction compensation is applied at the feedback level during motion control. Moreover, it was verified that the Coulomb-tanh friction force law is suitable and effectively converges the traditional Coulomb's friction law. The proposed methodology, involving the Coulomb-tanh model, is a way to obtain acceptable results, avoiding complementarity formulations and differential inclusions in the optimal control procedure.

**Author Contributions:** Conceptualization, J.A.R.-Q., F.D. and H.C.R.-d.-Á.; methodology, J.A.R.-Q., F.D., E.B. and N.R.C.-C.; software, J.A.R.-Q., H.C.R.-d.-Á. and B.S.-G.; validation, J.A.R.-Q., F.D. and E.B.; formal analysis, J.A.R.-Q. and F.D.; resources, J.A.R.-Q., E.B. and N.R.C.-C.; data curation, H.C.R.-d.-Á., E.B. and B.S.-G.; writing, J.A.R.-Q., F.D. and E.B.; funding acquisition, J.A.R.-Q., E.B. and N.R.C.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by TecNM research grants 14980.22-P and 18009.23-P.

**Data Availability Statement:** The raw data supporting the conclusions of this article will be made available by the authors on reasonable request.

**Acknowledgments:** This research was concluded while the corresponding author held a one-year visiting position at the Laboratoire des Sciences du Numérique de Nantes (LS2N) of the Institut Mines-Télécom Atlantique (IMTA) of Nantes, France (October 2023 to October 2024); under the authorization of CONAHCYT. The corresponding author acknowledges the support received from IMTA and CONAHCYT. This work was partially developed in the framework of the international network "Red Internacional de Control y Cómputo Aplicados" supported by TecNM, Mexico.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## Abbreviations

The following abbreviations are used in this manuscript:

OC Optimal control

PMP Pontryagin's Maximum Principle

ODE Ordinary differential equations
DOF Degree-of-freedom

FEM Finite Element Method

## References

1. Armstrong-Hélouvry, B.; Dupont, P.; De Wit, C.C. A survey of models, analysis tools and compensation methods for the control of machines with friction. Automatica **1994**, *30*, 1083–1138. [CrossRef]

2. Olsson, H.; Åström, K.; Canudas de Wit, C.; Gäfvert, M.; Lischinsky, P. Friction Models and Friction Compensation. Eur. J. Control **1998**, *4*, 176–195. [CrossRef]

3. Bona, B.; Indri, M. Friction Compensation in Robotics: An Overview. In Proceedings of the 44th IEEE Conference on Decision and Control, Seville, Spain, 15 December 2005; pp. 4360–4367. [CrossRef]

4. Andersson, S.; Söderberg, A.; Björklund, S. Friction models for sliding dry, boundary and mixed lubricated contacts. Tribol. Int. **2007**, *40*, 580–587. NORDTRIB 2004. [CrossRef]

5. Marques, F.; Flores, P.; Pimenta Claro, J.C.; Lankarani, H.M. A survey and comparison of several friction force models for dynamic analysis of multibody mechanical systems. Nonlinear Dyn. **2016**, *86*, 1407–1443. [CrossRef]

6. Khan, Z.A.; Chacko, V.; Nazir, H. A review of friction models in interacting joints for durability design. Friction **2017**, *5*, 1–22. [CrossRef]

7. Marques, F.; Flores, P.; Claro, J.C.P.; Lankarani, H.M. Modeling and analysis of friction including rolling effects in multibody dynamics: A review. Multibody Syst. Dyn. **2019**, *45*, 223–244. [CrossRef]

8. Rojas-Quintero, J.A.; Rojas-Estrada, J.A.; Villalobos-Chin, J.; Santibañez, V.; Bugarin, E. Optimal controller applied to robotic systems using covariant control equations. Int. J. Control **2022**, *95*, 1576–1589. [CrossRef]

9. Rojas-Quintero, J.A.; Dubois, F.; Ramírez-de Ávila, H.C. Riemannian Formulation of Pontryagin's Maximum Principle for the Optimal Control of Robotic Manipulators. Mathematics **2022**, *10*, 1117. [CrossRef]

10. Grancharova, A.; Johansen, T.A. Survey of Explicit Approaches to Constrained Optimal Control. In Switching and Learning in Feedback Systems: European Summer School on Multi-Agent Control, Maynooth, Ireland, 8–10 September 2003; Revised Lectures and Selected Papers; Murray-Smith, R., Shorten, R., Eds.; Springer: Berlin/Heidelberg, Germany, 2005; pp. 47–97. [CrossRef]

11. Rath, G.; Harker, M.; Zaev, E. Direct numerical solution of stiff ODE systems in optimal control. In Proceedings of the 2017 6th Mediterranean Conference on Embedded Computing (MECO), Bar, Montenegro, 11–15 June 2017; pp. 1–5. [CrossRef]

12. Almuslimani, I.; Vilmart, G. Explicit Stabilized Integrators for Stiff Optimal Control Problems. SIAM J. Sci. Comput. **2021**, *43*, A721–A743. [CrossRef]

13. Rojas-Quintero, J.A.; Villalobos-Chin, J.; Santibanez, V. Optimal Control of Robotic Systems Using Finite Elements for Time Integration of Covariant Control Equations. IEEE Access **2021**, *9*, 104980–105001. [CrossRef]

14. Stewart, D.E. Rigid-Body Dynamics with Friction and Impact. SIAM Rev. **2000**, *42*, 3–39. [CrossRef]

15. Tassa, Y.; Erez, T.; Todorov, E. Synthesis and stabilization of complex behaviors through online trajectory optimization. In Proceedings of the 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, Vilamoura-Algarve, Portugal, 7–12 October 2012; pp. 4906–4913. [CrossRef]

16. Posa, M.; Cantu, C.; Tedrake, R. A direct method for trajectory optimization of rigid bodies through contact. Int. J. Robot. Res. **2014**, *33*, 69–81. [CrossRef]

17. Carius, J.; Ranftl, R.; Koltun, V.; Hutter, M. Trajectory Optimization With Implicit Hard Contacts. IEEE Robot. Autom. Lett. **2018**, *3*, 3316–3323. [CrossRef]

18. Lin, F.; Brandt, R. An optimal control approach to robust control of robot manipulators. IEEE Trans. Robot. Autom. **1998**, *14*, 69–77. [CrossRef]

19. Braun, D.J.; Petit, F.; Huber, F.; Haddadin, S.; van der Smagt, P.; Albu-Schäffer, A.; Vijayakumar, S. Robots Driven by Compliant Actuators: Optimal Control Under Actuation Constraints. IEEE Trans. Robot. **2013**, *29*, 1085–1101. [CrossRef]

20. Pan, H.; Xin, M. Nonlinear robust and optimal control of robot manipulators. Nonlinear Dyn. **2014**, *76*, 237–254. [CrossRef]

21. Asgari, M.; Nikoobin, A. Analysis of Optimal Dynamic Manipulation for Robotic Manipulator Based on Pontryagin's Minimum Principle. Arab. J. Sci. Eng. **2020**, *45*, 9159–9169. [CrossRef]

22. Chignoli, M.; Wensing, P.M. Variational-Based Optimal Control of Underactuated Balancing for Dynamic Quadrupeds. IEEE Access **2020**, *8*, 49785–49797. [CrossRef]

23. Gaudin, H.; Bessonnet, G. From identification to motion optimization of a planar manipulator. Robotica **1995**, *13*, 123–132. [CrossRef]

24. Kaserer, D.; Gattringer, H.; Müller, A. Nearly Optimal Path Following With Jerk and Torque Rate Limits Using Dynamic Programming. IEEE Trans. Robot. **2019**, *35*, 521–528. [CrossRef]

25. Kim, T.H.; Ha, I.J. Time-optimal control of a single-DOF mechanical system with friction. IEEE Trans. Autom. Control **2001**, *46*, 751–755. [CrossRef]

26. Brillouin, L. Les Tenseurs en Mécanique et en Élasticité; Masson et Cie: Paris, France, 1938.

27. Lovelock, D.; Rund, H. Tensors, Differential Forms, and Variational Principles; Dover Books on Mathematics; Dover Publications: Mineola, NY, USA, 1989.

28. Grinfeld, P. Introduction to Tensor Analysis and the Calculus of Moving Surfaces; Springer: New York, NY, USA, 2013. [CrossRef]
29. Gu, Y.L. Modeling and simplification for dynamic systems with testing procedures and metric decomposition. In Proceedings of the 1991 IEEE International Conference on Systems, Man, and Cybernetics, Charlottesville, VA, USA, 13–16 October 1991; Volume 1, pp. 487–492. [CrossRef]

30. Spong, M. Remarks on robot dynamics: Canonical transformations and Riemannian geometry. In Proceedings of the 1992 IEEE International Conference on Robotics and Automation, Nice, France, 12–14 May 1992; Volume 1, pp. 554–559. [CrossRef]

31. Žefran, M.; Bullo, F. *Robotics and Automation Handbook*; CRC Press: Boca Raton, FL, USA, 2005; Chapter Lagrangian Dynamics.

32. Goldstein, H.; Poole, C.; Safko, J. *Classical Mechanics*; Addison-Wesley Series in Physics; Addison Wesley: Boston, MA, USA, 2002.

33. Arnold, V.I. *Mathematical Methods of Classical Mechanics*, 2nd ed.; Springer: New York, NY, USA, 1989.

34. Brogliato, B.; Lozano, R.; Maschke, B.; Egeland, O. *Dissipative Systems Analysis and Control*, 3rd ed.; Communications and Control Engineering; Springer: Cham, Switzerland, 2020. [CrossRef]

35. van der Schaft, A. *L2-Gain and Passivity Techniques in Nonlinear Control*, 3rd ed.; Communications and Control Engineering; Springer: Cham, Switzerland, 2017. [CrossRef]

36. Strutt, J.W. Some General Theorems relating to Vibrations. *Proc. Lond. Math. Soc.* **1871**, s1-4, 357–368. [CrossRef]

37. Minguzzi, E. Rayleigh's dissipation function at work. *Eur. J. Phys.* **2015**, 36, 035014. [CrossRef]

38. Mesterton-Gibbons, M. *A Primer on the Calculus of Variations and Optimal Control Theory*; American Mathematical Society: Providence, RI, USA, 2009; Volume 50. [CrossRef]

39. Pontryagin, L.S.; Boltyanskii, V.G.; Gamkrelidze, R.V.; Mishchenko, E.F. *The Mathematical Theory of Optimal Processes*; Interscience Publishers (Division of John Wiley & Sons, Inc.): New York, NY, USA, 1962.

40. Liberzon, D. *Calculus of Variations and Optimal Control Theory: A Concise Introduction*; Princeton University Press: Princeton, NJ, USA, 2012.

41. Wolfram Research, Inc. *Mathematica*, Version 12.3; Wolfram Research Inc.: Champaign, IL, USA, 2021.

42. Dubois, F.; Fortuné, D.; Rojas Quintero, J.A.; Vallée, C. Pontryagin Calculus in Riemannian Geometry. In *Proceedings of the Geometric Science of Information*; Nielsen, F., Barbaresco, F., Eds.; Springer: Cham, Switzerland, 2015; pp. 541–549. [CrossRef]

43. Stewart, D.E.; Anitescu, M. Optimal control of systems with discontinuous differential equations. *Numer. Math.* **2010**, *114*, 653–695. [CrossRef]

44. Shampine, L. Measuring Stiffness. *Appl. Numer. Math.* **1985**, *1*, 107–119. [CrossRef]

45. Spijker, M. Stiffness in numerical initial-value problems. *J. Comput. Appl. Math.* **1996**, *72*, 393–406. [CrossRef]

46. Kelly, R.; Santibáñez Davila, V.; Loría Perez, J.A. *Control of Robot Manipulators in Joint Space*; Advanced Textbooks in Control and Signal Processing; Springer: London, UK, 2005.

47. Grotjahn, M.; Daemi, M.; Heimann, B. Friction and rigid body identification of robot dynamics. *Int. J. Solids Struct.* **2001**, *38*, 1889–1902. [CrossRef]

48. Terze, Z.; Müller, A.; Zlatar, D. An Angular Momentum and Energy Conserving Lie-Group Integration Scheme for Rigid Body Rotational Dynamics Originating From Störmer–Verlet Algorithm. *J. Comput. Nonlinear Dyn.* **2015**, *10*, 051005. [CrossRef]

49. Terze, Z.; Müller, A.; Zlatar, D. Singularity-free time integration of rotational quaternions using non-redundant ordinary differential equations. *Multibody Syst. Dyn.* **2016**, *38*, 201–225. [CrossRef]

50. Dubois, F.; Antonio Rojas-Quintero, J. Simpson's Quadrature for a Nonlinear Variational Symplectic Scheme. In *Proceedings of the Finite Volumes for Complex Applications X—Volume 2, Hyperbolic and Related Problems*; Franck, E., Fuhrmann, J., Michel-Dansac, V., Navoret, L., Eds.; Springer: Cham, Switzerland, 2023; pp. 83–92. [CrossRef]

51. Rojas-Quintero, J.A.; Dubois, F.; Cabrera-Díaz, J.G. Simpson's Variational Integrator for Systems with Quadratic Lagrangians. *Axioms* **2024**, *13*, 255. [CrossRef]

52. Moreau, J.J. Evolution problem associated with a moving convex set in a Hilbert space. *J. Differ. Equ.* **1977**, *26*, 347–374. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.