Jerk Control of Floating Base Systems with Contact-Stable Parametrised Force Feedback

TL;DR

This paper introduces a jerk-based control framework for floating base systems that incorporates contact stability constraints and force feedback, transforming constrained optimization into an unconstrained problem for smoother control.

Contribution

It presents a novel jerk control architecture with a one-to-one contact stability manifold correspondence, enabling unconstrained optimization and improved force feedback integration.

Findings

Successfully validated with simulations and experiments on iCub robot.

Achieved contact-stable control with reduced discontinuities.

Effectively incorporated FT sensor feedback into control loop.

Abstract

Nonlinear controllers for floating base systems in contact with the environment are often framed as quadratic programming (QP) optimization problems. Common drawbacks of such QP based controllers are: the control input often experiences discontinuities; no force feedback from Force/Torque (FT) sensors installed on the robot is taken into account. This paper attempts to address these limitations using jerk based control architectures. The proposed controllers assume the rate-of-change of the joint torques as control input, and exploit the system position, velocity, accelerations, and contact wrenches as measurable quantities. The key ingredient of the presented approach is a one-to-one correspondence between free variables and an inner approximation of the manifold defined by the contact stability constraints. More precisely, the proposed correspondence covers completely the contact stability manifold except for the so-called friction cone, for which there exists a unique correspondence for more than 90% of its elements. The correspondence allows us to transform the underlying constrained optimisation problem into one that is unconstrained. Then, we propose a jerk control framework that exploits the proposed correspondence and uses FT measurements in the control loop. Furthermore, we present Lyapunov stable controllers for the system momentum in the jerk control framework. The approach is validated with simulations and experiments using the iCub humanoid robot.

Tables

Table 1. TABLE I : Robustness tests on the real iCub while performing highly dynamic movements.

Robustness w.r.t. modeling errors on real robot
Error on dynamic model $[\%]$	Success rate $\left[\mathrm{trials}\right]$
Overest. $5\%$	$5/5$
Overest. $7.5\%$	$3/5$
Overest. $10\%$	$1/5$
Underest. $5\%$	$1/5$

Equations

$$M(q)\dot{\nu }+C(q,\nu )\nu +G(q)=B\tau +\sum _{k=1}^{n_c}{J}_{C_k}^T{f}^k,$$

$$J(q)\nu =$$

$$J\dot{\nu }+\dot{J}\nu =0.$$

$$M\dot{\nu }+h$$

$$J\dot{\nu }+\dot{J}\nu$$

$$f_z>f_z^{min}\ge 0,$$

$$\sqrt[]{f_x^2+f_y^2}<{\mu }_c{f}_z,$$

$$y_c^{min}<\frac{M_x}{f_z}<y_c^{max},$$

$$x_c^{min}<-\frac{M_y}{f_z}<x_c^{max},$$

$$\left|\frac{M_z}{f_z}\right|<{\mu }_z.$$

$$\underset{y=(\dot{\nu },f,\tau )}{\mathrm{minimize}}$$

$$f^{\ast }=f_d-K_i{\int }_0^t(f_m-f_d)ds,\text{vspace}-0.08cm$$

$$\underset{y}{\mathrm{minimize}}$$

$$\underset{\xi }{\mathrm{minimize}}$$

$$\varphi (\xi )\text{mathchar}58=\left[\begin{array}{c}{\mu }_c\frac{\tanh ({\xi }_1)(e^{{\xi }_3}+f_z^{min})}{\sqrt[]{1+{\tanh }^2({\xi }_2)}}\\ {\mu }_c\frac{\tanh ({\xi }_2)(e^{{\xi }_3}+f_z^{min})}{\sqrt[]{1+{\tanh }^2({\xi }_1)}}\\ e^{{\xi }_3}+f_z^{min}\\ ({\delta }_y\tanh ({\xi }_4)+{\delta }_{y_0})(e^{{\xi }_3}+f_z^{min})\\ ({\delta }_x\tanh ({\xi }_5)+{\delta }_{x_0})(e^{{\xi }_3}+f_z^{min})\\ {\mu }_z\tanh ({\xi }_6)(e^{{\xi }_3}+f_z^{min})\end{array}\right],$$

$${\delta }_x\text{mathchar}58=\frac{x_c^{max}-x_c^{min}}{2},{\delta }_{x_0}\text{mathchar}58=-\frac{x_c^{min}+x_c^{max}}{2}$$

$${\delta }_y\text{mathchar}58=\frac{y_c^{max}-y_c^{min}}{2},{\delta }_{y_0}\text{mathchar}58=\frac{y_c^{max}+y_c^{min}}{2},$$

$$\Phi (\xi )$$

$$f$$

$$\dot{f}$$

$$\underset{y=(\dot{\nu },\xi ,\tau )}{\mathrm{minimize}}$$

$$D$$

$$\beta$$

$$\underset{\dot{y}=(\ddot{\nu },\dot{f},\dot{\tau })}{\mathrm{minimize}}$$

$$\underset{u=(\ddot{\nu },\dot{\xi },\dot{\tau })}{\mathrm{minimize}}$$

$$P\text{mathchar}58=\left(\begin{array}{ccc}{\mathbf{I}}_{n+6}& 0_{n+6,6n_c}& 0_{n+6,n}\\ 0_{6n_c,n+6}& \Phi (\xi )& 0_{6n_c,n}\\ 0_{n,n+6}& 0_{n,6n_c}& {\mathbf{I}}_n\end{array}\right).$$

$$\dot{H}$$

$$A$$

$$A_k$$

$$\ddot{H}$$

$$\dot{A}$$

$${\dot{A}}_k$$

$${\ddot{H}}^{\ast }={\ddot{H}}_d-K_d\dot{\tilde{H}}-K_p\tilde{H},$$

$$\ddot{H}(\dot{\xi })={\ddot{H}}^{\ast }$$

$${\dot{H}}_d=A{f}_e-mg{e}_3;$$

$${\ddot{H}}^{\ast }={\ddot{H}}_d-K_d\dot{\tilde{H}}-K_p\tilde{H}-K_i{\int }_0^t\tilde{H}dt,$$

Full text

I Introduction

Nonlinear feedback control of fixed-based (e.g. manipulators) and floating-base (e.g. humanoids) robots is not new to the Control community [1, 2, 3, 4]. Feedback-linearisation, robust control, and adaptive laws are only few examples of the large variety of control methods developed for steering these systems towards desired quantities. When fixed and floating-base systems make contact with the environment, the robot control has to deal with regulating also the forces and torques towards values that ensure a desired interaction. This paper contributes towards the stabilisation of floating-base systems in contact with the environment by proposing controllers that ensure contact stability (see, e.g., [5]) while including force feedback in the control laws. The proposed approach uses the rate-of-change of the joint torques as control input, and for this reason it is here referred to as jerk control.

Force control strategies for fixed-base systems can be roughly divided into two categories: direct and indirect force control methods [1, 6]. Indirect methods achieve compliance without explicitly closing a feedback loop on the measured contact forces. In particular, impedance control is a common objective for indirect techniques whose goal is often that of achieving a desired dynamic behavior of the robot’s end-effector. The control design requires to include force and torque measurements as feedforward terms to achieve full feedback linearization of the end-effector dynamics. If no force measurement is available, a simplified stiffness control can still be applied [1]. On the other hand, direct force control methods include explicit feedback from the measured interaction wrenches, usually related to the force error [6, 7]. An example of these techniques is the hybrid position/force control, which is often applied when the environment is rigid and the end-effector has to continuously maintain a contact. The rigid environment assumption enables the decomposition of the end-effector dynamics into a constrained and a free direction [8]. Along the constrained directions, feasible desired forces are exerted. Here, additional feedback terms are added to ensure convergence in presence of external disturbances and unmodeled dynamics. Although the desired force satisfies the contact stability constraints, the commanded force, which includes the feedback from Force/Torque (FT) sensors, may instantaneously violate the contact constraints.

The recent research effort on humanoid robots gave impetus to the force control of floating-base systems [9, 10, 11, 12]. These systems are often underactuated, namely, the number of control inputs is fewer than the system’s degrees of freedom [13]. The lack of actuation is usually circumvented by means of the contacts between the robot and the environment, but this requires close attention to the forces the robot exerts at the contact locations. If not regulated appropriately, uncontrolled contact forces may break a contact, which makes the robot control critical [10, 11]. The contacts between a floating-base system and the environment are often assumed to be rigid and planar, although compliant contacts and uneven ground are also considered in the literature [14, 15, 16, 17]. Furthermore, all contacts are also unilateral, being the robot not physically attached to ground and in general able to make and break contacts. Contact activation and deactivation occur continuously, e.g. in the case of humanoid robot walking, and can be addressed with the design of a proper state machine that plans references for a balancing/walking controller [18, 19]. From the control design perspective, a common strategy for floating-base systems is based on the so called stack-of-task approach [20]. These strategies usually consider several control objectives organized in a hierarchical or weighted prioritization. Often, the high-priority task is the stabilization of the robot momentum [21]: its objective is the control of the robot center-of-mass and angular momentum by means of contact forces while guaranteeing stable zero-dynamics [22]. Quadratic programming (QP) solvers can be used to monitor contact forces to ensure both robot and contact stability [23].

Analogously to direct and indirect force controllers, QP based force control of floating base systems usually suffers from the following limitations: $i)$ optimal control inputs (i.e., the joint torques) may be discontinuous in certain cases, e.g. during the switching from two feet and one foot balancing of humanoid robots. In this case, the robot control may become critical. Although further (and often numerous) constraints can be added to the QP solver to enforce continuity, the effectiveness of this approach is not always satisfactory in practice; $ii)$ force feedback from FT sensors is missing in the control action.

This paper proposes a control approach that attempts to address these two main limitations of QP based controllers for floating base systems. The key ingredients of our approach are: $a)$ propose an invertible one-to-one mapping between a set of constraint-free variables and an inner approximation of the contact stability manifold. $b)$ propose controllers that use the wrench mapping and exploit the rate-of-change of the joint torques as control input; $c)$ extend the proposed controllers to the case when the joint torques (and not their rate-of-change) is the available control input. The proposed approach exploits a relative degree augmentation of the underlying system, i.e., the system state is composed of the system position, velocity, and acceleration. For this reason, the proposed approach is referred to as jerk control, which also incorporates force feedback from FT sensors. Furthermore, we present control laws that stabilise a desired robot momentum using joint-torques as input and having guaranteed Lyapunov stability properties. The proposed approach is validated with simulations and experiments using the humanoid robot iCub balancing on rigid contacts.

The paper is organized as follows. Section II recalls the notation, the robot modeling, the contact stability constraints, and introduces the problem statement. Section III presents a novel contact-stable parametrization of the contact wrenches. Section IV introduces the main ideas behind jerk control using the contact-stable parametrisation. Section V presents Lyapunov stable jerk controllers when the control objective is the stabilisation of the robot momentum. Section VI presents validations of the approach on the humanoid robot iCub. Conclusions and perspectives conclude the paper.

II Background

II-A Notation

• •

$\mathcal{I}$ is an inertial frame, with its $z$ axis pointing against the gravity, $\mathcal{B}$ is a frame attached to the robot’s base link.

• •

The constant $m$ represents the total mass of the robot, and $g$ is the norm of the gravitational acceleration.

• •

Given a matrix $A\in\mathbb{R}^{p\times n}$, we denote with $A^{\dagger}\in\mathbb{R}^{n\times p}$ its Moore-Penrose pseudoinverse.

• •

$e_{i}\in\mathbb{R}^{6}$ is the canonical vector, consisting of all zeros but the $i$-th component that is equal to one.

• •

Let $S(x)\in\mathbb{R}^{3\times 3}$ be the skew-symmetric matrix such that $S(x)y{=}x{\times}y$, with $\times$ the cross product operator in $\mathbb{R}^{3}$.

II-B Robot Modeling

The robot is modeled as a multi-body system composed of $n+1$ rigid bodies, called links, connected by $n$ joints with one degree of freedom each. We also assume that none of the links has an a priori constant position-and-orientation with respect to an inertial frame, i.e. the system is free floating.

The robot configuration space is the Lie group $\mathbb{Q}=\mathbb{R}^{3}\times\mathbb{SO}(3)\times\mathbb{R}^{n}$ and it is characterized by the pose (i.e. position-and-orientation) of a base frame attached to the robot’s base link, and the joint positions. An element $q\in\mathbb{Q}$ can be defined as the following triplet: $q=(^{\mathcal{I}}o_{\mathcal{B}},{}^{\mathcal{I}}R_{\mathcal{B}},s)$ where ${}^{\mathcal{I}}o_{\mathcal{B}}\in\mathbb{R}^{3}$ denotes the origin of the base frame with respect to the inertial frame, ${}^{\mathcal{I}}R_{\mathcal{B}}$ is the rotation matrix representing the orientation of the base frame, and $s\in\mathbb{R}^{n}$ is the joint configuration characterizing the robot posture. The velocity of the multi-body system can be characterized by the algebra of $\mathbb{Q}$. We here choose to represent the velocity of the multi-body system by the set $\mathbb{V}=\mathbb{R}^{3}\times\mathbb{R}^{3}\times\mathbb{R}^{n}$, where an element of $\mathbb{V}$ is the vector $\nu=(^{\mathcal{I}}\dot{o}_{\mathcal{B}},\,^{\mathcal{I}}\omega_{\mathcal{B}},\,\dot{s})=(\mathrm{v}_{\mathcal{B}},\dot{s})$, and ${}^{\mathcal{I}}\omega_{\mathcal{B}}$ is the angular velocity of the base frame expressed w.r.t. the inertial frame, i.e. ${}^{\mathcal{I}}\dot{R}_{\mathcal{B}}=S(^{\mathcal{I}}\omega_{\mathcal{B}})^{\mathcal{I}}R_{\mathcal{B}}$. A more detailed description of the configuration space and its algebra is provided in [24].

We assume that the robot interacts with the environment by exchanging $n_{c}$ distinct wrenches. The equations of motion of the multi-body system can be described by applying the Euler-Lagrange formalism [25, Chapter 13.5]:

[TABLE]

where $M\in\mathbb{R}^{n+6\times n+6}$ is the mass matrix, $C\in\mathbb{R}^{n+6\times n+6}$ is the matrix accounting for Coriolis and centrifugal effects, $G$ $\in\mathbb{R}^{n+6}$ is the gravity term, $B=[0_{n\times 6},\,\,1_{n}]^{\top}$ is a selector of the actuated DoFs, $\tau\in\mathbb{R}^{n}$ is a vector representing the internal actuation torques, and $f^{k}\in\mathbb{R}^{6}$ represents an external wrench applied by the environment to the link of the $k$-th contact. The Jacobian $J_{\mathcal{C}_{k}}=J_{\mathcal{C}_{k}}(q)$ is the map between the robot’s velocity $\nu$ and the linear and angular velocity of the $k$-th contact link.

Lastly, it is assumed that a set of holonomic constraints acts on Eq. (1) of the form $c(q)=0$: they may represent, for instance, a frame having a constant position-and-orientation w.r.t. the inertial frame. Hence, we represent the holonomic constraints on links in rigid contact with the environment as $J_{\mathcal{C}_{k}}(q)\nu=0$. The holonomic constraints associated with all the rigid contacts can be then compactly represented as:

[TABLE]

By differentiating the kinematic constraints Eq. (2), one has:

[TABLE]

By combining the system dynamics (1) and the constraint equations (3), one obtains the following set of equations:

[TABLE]

where $h\mathrel{\mathop{\mathchar 58\relax}}=C(q,\nu)\nu+G(q)$, and $f\mathrel{\mathop{\mathchar 58\relax}}=[f^{1};...;f^{n_{c}}]\in\mathbb{R}^{6n_{c}}$ is the vector of contact wrenches making Eq. (3) satisfied.

II-C Contact Stability Constraints

Often, the holonomic constraints acting on the system represent a robot flat surface in complete rigid contact with the environment. To maintain the constraints stable, then, conditions on the contact forces and moments shall be met. More precisely, let $f^{k}\in\mathbb{R}^{6}$ denote the vector of the $k$-th contact wrench associated with the $k$-th active rigid contact of a planar robot surface, namely $f^{k}=[f_{x},\hskip 2.84526ptf_{y},\hskip 2.84526ptf_{z},\hskip 2.84526ptM_{x},\hskip 2.84526ptM_{y},\hskip 2.84526ptM_{z}]^{\top}$. Then, the contact stability conditions to preserve the planar unilateral contact can be formulated as follows:

[TABLE]

Being the constraints unilateral, condition (6a) imposes that the force normal to the contact is greater than a value $f_{z}^{min}$, which must be greater than or equal to zero. Eq. (7a) limits the magnitude of the forces parallel to the contact surface not to overcome the static friction characterised by the coefficient $\mu_{c}$. Conditions (8a)-(9a) constrain the local Center of Pressure – see, e.g., [5, Appendix B, Eq. (A7)] – to remain inside the contact surface, which is assumed to be rectangular of dimensions $x_{c}^{min},x_{c}^{max},y_{c}^{min},y_{c}^{max}$, calculated w.r.t. the contact reference frame $\mathcal{C}$ and defined as shown in Fig. 1. Eq. (10a) imposes no foot rotation along the axis normal to the contact surface, and $\mu_{z}$ is the torsional friction coefficient.

II-D Problem Statement

A common control approach for system (4a) usually considers several control objectives organized in a hierarchical or weighted prioritization [20, 22]. More precisely, let $a^{*}$ be a desired acceleration that the system should achieve. Then, a single priority111When several priorities are defined into the optimisation problem, higher priority tasks can be defined as constraints of (11). stack-of-task can be represented by the following optimisation problem:

[TABLE]

with $A$ a proper projection matrix, and $\mathcal{K}$ the manifold given by the constraints (6a). The above optimisation problem is usually framed as a Quadratic Programming (QP) one, and its solutions may suffer from the following limitations:

The solution may be discontinuous, e.g. at contact switching or after reference trajectory sharp variations; 2. 2.

The closed-loop dynamics does not include any feedback term from the measured contact wrenches $f_{m}$.

Limitation $1)$ is often addressed by approximating the continuity property with a set of inequality constraints to be added to (11), but the effectiveness of this approach is often unsatisfactory from the experimental standpoint [26]. Limitation $2)$ is the most critical one, since FT sensor information are not used in the optimal control law $\tau$ that solves (11), thus potentially wasting important feedback information at the control level.

Let us observe that Limitation $2)$ may be attenuated when desired force tasks are added to the problem (11) [27]. For instance, if we aim to achieve a desired force $f_{d}$, then the force task can be achieved by adding equality constraints in the form $f=f_{d}$ to the problem (11). At this point, one may attempt at using the FT measurements by replacing $f_{d}$ with

[TABLE]

where $f_{m}$ is the measured force, and $f=f^{*}$ being the equality constraint. The main limitation of this approach is that this equality constraint may require $f$ to violate the constraint $f\in\mathcal{K}$. Putting the desired force as part of the cost function of (11) may be an option, but this alters the priorities that the force task has over the acceleration one.

What follows presents an alternative, theoretically sound approach that aims at addressing the above limitations $1)$, and $2)$ of classical QP based stack-of-task approaches for the control of floating base systems in contact with the environment.

III A Contact-Stable Wrench Parametrization

Parametrisations can be an effective way to transform constrained optimisation problems into unconstrained ones [28]. Consider, for instance, the following optimisation problem:

[TABLE]

with $\text{cost}(\cdot)\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}{\rightarrow}\mathbb{R}$, and $y\in\mathbb{R}$. If there exists a solution to (12), the process of seeking for this solution is equivalent to solving the following problem:

[TABLE]

with $\xi\in\mathbb{R}$. For the specific case (12), it is trivial to find a parametrisation ensuring $y>0$. Note, however, that the mapping $y=e^{\xi}$ is one-to-one, and its gradient is always invertible, namely $\frac{\partial}{\partial\xi}(e^{\xi})=e^{\xi}\neq 0\ \forall\xi$. These two additional properties are of particular importance for the numerical stability of solvers addressing the problem (13).

Next Lemma proposes a wrench parametrisation that may be used to remove the constraint $f\in\mathcal{K}$ into the problem (11).

Lemma 1**.**

Consider the parametrization $f^{k}{=}\phi(\xi)\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{6}\rightarrow\mathcal{K^{\prime}}$ defined by

[TABLE]

with $f_{z}^{min}\geq 0$ the minimum value of the vertical force $f_{z}$,

[TABLE]

and $x_{c}^{max},y_{c}^{max}$ and $x_{c}^{min},y_{c}^{min}$ the contact surface dimensions as described in Fig. 1.

Then, the following properties hold:

The contact constraints (6a) are always satisfied, namely, $\mathcal{K^{\prime}}\subset\mathcal{K}$ or, equivalently, $\phi(\xi)\in\mathcal{K}\ \ \forall\xi\in\mathbb{R}^{6}$. 2. 2.

The function $\phi(\xi)\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{6}\rightarrow\mathcal{K^{\prime}}$ is a bijection, namely, a one-to-one correspondence from $\mathbb{R}^{6}$ to $\mathcal{K^{\prime}}$. 3. 3.

*The gradient of the function $\phi(\xi)$, i.e. *

[TABLE]

is an invertible matrix $\forall\hskip 1.42262pt\xi\in\mathbb{R}^{6}$.

The proof is in the Appendix -A. Lemma 1 shows that (14) generates wrenches that always belong to the contact stability manifold defined by (6a). Furthermore, it shows that $\phi(\cdot)$ is a one-to-one correspondence between a set of free parameters $\xi$ and the manifold $\mathcal{K^{\prime}}$, i.e. the image of (14). Clearly, one may find other functions for which the contact constraints (6a) are always satisfied. The proposed function $\phi(\cdot)$ in (14), however, has an image that approximates the friction cones (7a) similarly to what one obtains using a set of linear inequalities.

More precisely, when (14) is substituted into (7a), instead of a cone we obtain a set that much resembles an octagonal flipped pyramid. In fact, it is easy to verify that $\sqrt[]{f_{x}^{2}+f_{y}^{2}}=\mu_{c}f_{z}$ is satisfied eight times when $f^{k}{=}\phi(\xi)$ and either, or both, $\tanh^{2}(\cdot)=1$. By computing the ratio between the areas of an octagon and that of a circle, we conclude that (14) covers more than 90 % of the set given by (7a). For a visual representation of this fact, see Figure 2 that shows a typical approximation of (7a) given by the function (14). Let us also observe that the gradient of this function is invertible for any value of the parameter $\xi$. This property will be of pivotal importance in Sections IV and V when designing stable controllers for system (4a).

The parametrization (14) can be easily extended in case of $n_{c}$ distinct contact wrenches. In this case, define:

[TABLE]

where $\Phi=\text{blkdiag}(\Phi_{1},...\hskip 1.42262pt,\Phi_{n_{c}})\in\mathbb{R}^{6n_{c}\times 6n_{c}}$ and $\xi=[\xi^{1};...\hskip 1.42262pt;\xi^{n_{c}}]\in\mathbb{R}^{6n_{c}}$. It is then straightforward to verify that the properties described in Lemma 1 are retained even in case of multiple contact wrenches.

IV Jerk Control

This section proposes control laws that exploit the contact wrench parametrisation (14) and attempt to address the limitations – listed in Section II-D – of the classical torque-based controllers framed as stack-of-tasks optimisation problems.

IV-A *

Jerk control with parametrised contact wrenches*

The wrench parametrisation (14) can be used to remove the constraint $f\in\mathcal{K}$ from the optimisation problem (11). This process would lead to the following formulation:

[TABLE]

with $g(y)\mathrel{\mathop{\mathchar 58\relax}}=[\dot{\nu};\phi(\xi);\tau]$. The main drawbacks of the above approach are: $i)$ the optimisation problem (20) can no longer be casted in a QP being the parametrisation $\phi(\xi)$ nonlinear; we would then need nonlinear – and often slower than QPs – optimisers to solve (20); $ii)$ the limitations 1), 2) listed in Section II-D are not addressed.

To include feedback terms into the control laws, the contact wrenches, or accelerations, shall become part of the system state. In the language of Automatic Control, we shall then proceed with augmenting the relative degree of the output (or task) that one wants to stabilise [29].

More precisely, assume that: $hp{-}i)$ the control objective is the stabilisation of a desired jerk $\dot{a}^{*}$; $hp{-}ii)$ the joint torque rate-of-change $\dot{\tau}$ can be considered as a control input; $hp{-}iii)$ both the joint torques $\tau$ and the contact forces $f$ are measurable quantities, so the robot acceleration $\dot{\nu}$ – if not measurable – can be obtained from (4a). Now, define

[TABLE]

As a consequence of $hp{-}iii)$, one has a measurement of the vector $y=(\dot{\nu},f,\tau)$, while the variable $\dot{y}$ can be used as a search variable. Then, control laws for $\dot{\tau}$ that contain feedback information from the FT sensors are obtained as an outcome of the following optimisation problem:

[TABLE]

The solutions to the above problem are continuous in $\tau$ (even if $\dot{\tau}$ is discontinuous) and contain FT measurement feedback from the vector $y$. One of the main difficulties when solving (23a) is given by (24a). Since the variable $f$ no longer is a search variable, in fact, one cannot instantaneously choose values of the contact wrenches such that $f\in\mathcal{K}$. One may attempt to make (24a) satisfied by regulating appropriately the variable $\dot{f}$, which influences the wrench $f$ at the next time step.

A possibility to make (24a) satisfied is to use the parametrisation in Lemma 1: the gradient of the parametrisation automatically enforces the fact that $f(t)\in\mathcal{K}\ \forall t$. More precisely, in view of (17), one has $\dot{y}=(\ddot{\nu},\Phi(\xi)\dot{\xi},\dot{\tau})$, which leads to the following optimisation problem:

[TABLE]

with $P$ defined as:

[TABLE]

In order to be solved at each time instant, the optimisation problem (25a) requires the variable $\xi$. This variable may be retrieved from either time integration of $\dot{\xi}$, or by inverting the relationship $f=\phi(\xi)$: being the parametrisation a one-to-one correspondence (see Lemma 1), there exists a unique $\xi$ for any value of the contact wrench $f$ provided that it belongs to $\mathcal{K^{\prime}}\subset\mathcal{K}$. The latter allows us to inject further information from the FT sensor measurements into the optimal control laws $u$ solving the optimisation problem (25a).

Note also that the matrix $P$ is invertible thanks to the property $3)$ of Lemma 1. The invertibility of $P$ clearly plays a pivotal role when solving the optimisation problem (25a).

IV-B On the modeling and control requirements for jerk control

The optimal value $\dot{\tau}$ solving (25a) may be sent directly to the robot low-level control system if it allows to set desired rate-of-changes of joint torques. This may be feasible, for instance, when the low-level robot control exploits the model between the joint torques $\tau$ and the motor currents $i$, e.g. $\tau=k_{\tau}i$. More precisely, the motor currents are usually subject to electrical dynamics of the kind $\tfrac{d}{dt}{i}=k_{i}v$, where $v$ is often the motor voltages to be applied to the motors – namely, the real control input. Then, it is straightforward to express the optimisation problem (25a) so that the search variable $u$ contains $v$. Let us observe, however, that this control architecture in general requires high-frequency control loops (e.g. $5-20$ KHz) for generating the motor voltages $v$: these loops have to compute inverse dynamics within a short control cycle. If the control loops are not fast enough, sampling effects may be preponderant phenomena that render the assumption $\tfrac{d}{dt}{i}~{}=~{}k_{i}v$ not representative of the underlying physical dynamics. In this case, the associated control strategy resulting from (25a) may prove to be ineffective.

Another necessary requirement for achieving jerk control is the calculation of the terms $\dot{A}$, $\dot{D}$ and $\dot{\beta}$ to solve the optimisation problem (25a). These terms in general depend on the robot configuration space $q$, velocity $\nu$, and accelerations $\ddot{\nu}$, and need the derivatives of the system inverse dynamics. Besides numerical approximations for computing these terms, existing libraries nowadays provide users with the support of automatic differentiation and/or directly derivatives of inverse dynamics [30, 31]. If some of the terms in (25a) are not available, one may attempt setting them equal to zero and tune the feedback control gains in $\dot{a}^{\ast}$ so as to achieve robustness against them. However, we present below a jerk control architecture that overcomes the above modeling and control limitations of the mere application of (25a).

V Momentum-based Jerk control

This section proposes control laws that can be obtained from the problem (25a) when it is explicitly solved and extended for a two layer stack-of-task. These laws can also be shown to possess stability properties. Interestingly, the architecture presented below does not need the feedforward terms that depend on the inverse dynamics derivatives required by (25a). This is achieved by loosing the continuity property of $\tau$ but retaining the continuity of the contact wrenches $f$.

More precisely, we assume that: i) the highest priority task is the stabilisation of a desired robot centroidal momentum [24, 32]; ii) the lower priority task aims at stabilising the robot posture to regulate the system zero dynamics [29].

Let us recall that the momentum rate-of-change equals the summation of all the external wrenches acting on the robot. In a multi-contact scenario, the external wrenches reduce to the contact wrenches plus the gravity force:

[TABLE]

where $H\in\mathbb{R}^{6}$ is the robot’s momentum, $A_{k}\in\mathbb{R}^{6\times 6}$ is the matrix mapping the $k$-th contact wrench to the momentum dynamics, ${}^{\mathcal{I}}o_{\mathcal{C}_{k}}\in\mathbb{R}^{3}$ is the origin of the frame associated with the $k$-th contact, and ${}^{\mathcal{I}}o_{CoM}\in\mathbb{R}^{3}$ is the CoM position.

Recall that the rate-of-change of the robot momentum (27a) is related to the system accelerations (e.g. acceleration of the system center of mass). So, to derive jerk-based control laws, we need to differentiate (27a) w.r.t. time, which writes:

[TABLE]

Note that Eq. (30) is linear w.r.t. $\dot{\xi}$ and optimisation problems similar to (25a) may be laid down. In particular, to obtain the robot momentum stability, one may: $i)$ consider $\dot{\xi}$ as control input – or search variable – of the momentum acceleration (30); $ii)$ apply feedback linearization to (30) in order to impose a momentum acceleration $\ddot{H}^{*}$ of the form:

[TABLE]

where $K_{d},K_{p}\in\mathbb{R}^{6\times 6}$ are symmetric and positive definite matrices, $H_{d}\in\mathbb{R}^{6}$ is the reference momentum, and $\tilde{H}=H-H_{d}$ is the momentum error. Observe that it is always possible to find $\dot{\xi}$ such that

[TABLE]

because of the item $3)$ of Lemma 1. More precisely, the gradient $\Phi$ being always invertible ensures that the matrix $A\Phi$ in (30) is full rank $\forall\ \xi$. Consequently, $\dot{\xi}$ has full control authority on the momentum acceleration for any value of $\xi$. Clearly, one can impose (32) as long as $\xi$ remains bounded.

The stability properties of the control laws ensuring (32) are then presented below.

Lemma 2**.**

Assume that:

• •

the robot makes at least one rigid contact with the environment, i.e. $n_{c}\geq 1$;

• •

the desired momentum $H_{d}$ is a feasible system equilibrium such that there exists $f_{e}(t)\in\mathcal{K^{\prime}}$ satisfying

[TABLE]

• •

the variable $\dot{\xi}$ is chosen so as (32) is satisfied, with $\ddot{H}$ given by (30) and $\ddot{H}^{*}$ by (31).

Then:

the equilibrium point $(\tilde{H},\dot{\tilde{H}})=(0,0)$ is locally asymptotically stable if the robot makes one rigid contact with the environment, i.e. $n_{c}=1$; 2. 2.

the equilibrium point $(\tilde{H},\dot{\tilde{H}})=(0,0)$ is globally asymptotically stable if $\xi$ is bounded, namely if there exists a constant $c\in\mathbb{R}^{+}$ such that $|\xi(t)|<c\ \forall t$.

The proof is in the Appendix -B. Lemma 2 shows that as long as we satisfy Eq. (32), the system trajectories converge towards the desired values. The possibility of satisfying Eq. (32) is inherently related to the boundedness of $\xi$, which guarantees that the matrix $A\Phi(\xi)$ in (30) remains of full rank. More precisely, the item $1)$ shows that the boundeness of $\xi$ is achieved locally to the equilibrium point when the number of contacts is equal to one. When the number of contacts is greater than one, there is a redundancy for $\dot{\xi}$ that solve (32), and this redundancy should be chosen so as $\xi$ is always bounded. In this case, in fact, the item $2)$ of Lemma 2 shows that the system trajectories globally converge towards the desired values. We show in the next section a possible choice for the redundancy of $\dot{\xi}$ that proved to work effectively both in simulation and in real experiments. Let us remark, however, that the aim of a global bounded $\xi$ cannot be in general achieved, and so global asymptotic stability. For instance, think of a humanoid robot standing on one foot and starting with an high velocity of its center of mass: there is a limit for this velocity that would cause the contact to break, namely, the variable $\xi$ growing indefinitely. However, we can monitor the adverse conditions, where contacts are about to break, by looking at the overall norm of the variable $\xi$. This is, in the authors’ opinion, an interesting result of the proposed approach.

Often, the control objective of the center-of-mass trajectory tracking is framed as momentum control. In this case, one is tempted to use the approach presented above with

[TABLE]

where the (linear momentum) integral correction terms can be replaced by the position errors between the center of mass trajectory and its desired values. The resulting third order system (32)-(33), however, is in general very sensitive to gain tuning, as not all possible combinations of the gain matrices guarantee stability of the associated closed-loop system. This limitation affects the controller’s performances when applied to the real robot, where phenomena as modeling errors, measurements noise and external disturbances further limit the control gain choice.

V-A Momentum-based jerk control with integral terms

We propose a control algorithm alternative to pure feedback linearization with the goal of facilitating the gain tuning of the closed-loop system dynamics. In particular, consider as control objective the stabilization of $(I,\tilde{H},\zeta)$ towards the reference values $(0,0,0)$, with $I$ the integral of the momentum error, $\tilde{H}$ the momentum error, and $\zeta$ an exogenous state that will be used to prove Lyapunov stability, i.e.

[TABLE]

whose dynamics write:

[TABLE]

The dynamics (37a) are obtained by taking the derivative of $(I,\tilde{H},\zeta)$, and substituting their corresponding definitions in (34a)-(36a). Then, the following result holds.

Lemma 3**.**

Assume that:

• •

the robot makes at least one rigid contact with the environment, i.e. $n_{c}\geq 1$;

• •

the desired momentum $H_{d}$ is a feasible system equilibrium such that there exists $f_{e}(t)\in\mathcal{K^{\prime}}$ satisfying

[TABLE]

Choose,

[TABLE]

with $K_{o},K_{p},K_{d}{\in}\mathbb{R}^{6\times 6}$ symmetric positive definite matrices,

[TABLE]

the projector in the null space of $A\Phi$, and $\dot{\xi}_{0}$ a free variable of proper dimension. Then:

the equilibrium point $(I,\tilde{H},\zeta)=(0,0,0)$ is locally asymptotically stable if the robot makes one rigid contact with the environment, i.e. $n_{c}=1$; 2. 2.

the equilibrium point $(I,\tilde{H},\zeta)=(0,0,0)$ is globally asymptotically stable if $\xi$ is bounded, namely if there exists a constant $c\in\mathbb{R}^{+}$ such that $|\xi(t)|<c\ \forall t$.

The proof is in the Appendix -C. Lemma 3 shows that the additional integral terms do not break the stability properties achieved in Lemma 2, and that one is left with free positive definite gains $K_{o},\ K_{p},\ K_{d}$. Let us recall that the constraints (6a) remain satisfied while ensuring the stability properties of the associated closed-loop system, and such a claim cannot usually be made in classical stack-of-task approaches (11).

The control law (40) contains both feedforward and feedback terms that depend on the measured contact wrenches. It makes use, in fact, of (27a) for computing $\dot{H}$, which depends on the measured contact wrenches. In case of a single contact, there exists a unique control input $\dot{\xi}$ that satisfies (40), and the null space of the matrix $A\Phi$ is empty, i.e. $N_{A\Phi}=0$. In case of multiple contacts ($n_{c}>1$), instead, infinite control inputs satisfy (40). We solve the associated redundancy using the free variable $\dot{\xi}_{0}$ to minimize the norm of the robot joint torques. The computation of $\dot{\xi}_{0}$ is detailed in Appendix -D.

Let us remark again the importance of the invertibility of the gradient $\Phi$ – see Lemma 1. This property guarantees that the matrix $A\Phi$ in (40) is full rank, so $\dot{\xi}$ has full control authority on the momentum acceleration for any value of $\xi$.

V-B Computation of $f$, $\dot{H}(f)$ and $\Phi(\xi)$

The control input (40) requires: the contact wrenches $f$; the momentum derivative $\dot{H}=\dot{H}(f)$; and the associated variable $\xi$ such that $f=\phi(\xi)$. The contact wrenches can be measured/estimated using the measurements from 6-axis FT sensors installed on the robot. Once the wrenches $f$ are retrieved, we can compute the momentum rate of change via (27a) The associated $\xi$ can be computed by applying the parametrisation inverse mapping, namely $\xi=\phi^{-1}(f)$. The inverse mapping exists provided that the measured contact wrenches remain inside the set $\mathcal{K^{\prime}}$. If the measured wrenches do not belong to $\mathcal{K^{\prime}}$ (because, e.g., measurements noise, external unmodeled disturbances, etc.), a saturation shall be applied in the calculation of the inverse mapping so that the control input ${\xi}$ always remains finite.

V-C Computation of the joint torques to realize $\dot{\xi}$

To realize a $\dot{\xi}$, e.g. the law in (40), we have to choose the real control input of the system properly. We assume in this section that the control input is the joint torque $\tau$, so we cannot impose a desired $\dot{\tau}$ instantaneously.

A possibility for computing the joint torques is to find $\dot{\tau}$ realising $\dot{\xi}$, and then perform time-integration of $\dot{\tau}$ to obtain $\tau$. This procedure, however, requires to compute some derivatives of the inverse dynamics, which may not be available in practice.

For this reason, we follow here another route for finding the joint torques $\tau$ attempting to realise $\dot{\xi}$. First, we find the instantaneous relationship between the joint torques $\tau$ and the contact wrenches $f$. This relationship can be found, for instance, by substituting the state accelerations $\dot{\nu}$ from (4a) into the constraints (5a), which leads to:

[TABLE]

with $\Lambda=JM^{-1}B$. Then, we proceed as follows:

• •

Integrate the control input $\dot{\xi}$ to obtain $\xi$. The initial conditions for the integrator can be calculated by measuring the initial contact forces $f(0)$ and by applying the parametrization inverse mapping, i.e. $\xi(0)=\phi^{-1}(f(0))$;

• •

Apply the parametrization direct mapping to evaluate the wrenches $f$ from $\xi$, i.e. $f=\phi(\xi)$. By doing so, note that $f$ always satisfy the contact stability constraints;

• •

Retrieve the input torques $\tau$ from (41), which write

[TABLE]

where $N_{\Lambda}=(1_{n}-{\Lambda}^{\dagger}\Lambda)$ is the projector in the null space of $\Lambda$, and $\tau_{0}$ is a free variable, that can be chosen as in [22] to guarantee the stability of the system’s zero-dynamics.

VI Simulations and Experimental Results

VI-A Simulation Environment

The modeling and control framework presented in Sec. V is tested on the 23-DoFs iCub humanoid robot [33], both on the real robot and in simulations using Gazebo [34]. The controller is implemented in Simulink, and runs at a frequency of $100~{}\text{[Hz]}$. An advantage of using the Simulink-Gazebo simulator consists in the possibility to test directly on the real robot the same control software used in simulation.

Gazebo offers different physic engines to integrate the system’s dynamics. We chose the Open Dynamics Engine (ODE), that uses a fixed step semi-implicit Euler integration scheme, with a maximum simulation time step of $1$ [ms].

On the real iCub, the Simulink controller runs on an external PC and provides reference joint torques to an inner joint torque control loop, that runs on board the robot at $1000~{}\text{[Hz]}$. At the moment, iCub is not endowed with joint torque sensors. The measured joint torques are achieved by combining the FT sensor information, the joint encoders, IMU data and the robot model. The robot is equipped with 6 FT sensors: two located in the robot’s upper arms, two in the robot thighs, and two in the robot’s feet [33]. During experiments on the real robot, we verified that the desired controller step time of $0.01\ \rm[s]$ is respected most of the times. More precisely, statistics show that the desired step time is respected around $97.9\%$ of iterations for the jerk control. This is comparable with the performances of the classical momentum-based QP control [22], which meets the desired step time around $98.2\%$ of iterations.

The choice of running the high-level Simulink controller with a frequency of $100~{}\text{[Hz]}$ is due to: a CPU limitation of the current iCub; the limited bandwidth of the current FT sensors at the robot’s feet. However, we are working on robot hardware upgrades that will allow us to run the Simulink controllers with a standard frequency of $1000~{}\text{[Hz]}$.

VI-B Robustness analysis for jerk control

A preliminary analysis of the robot balancing behavior with the control law (40)–(42) indicates that the proposed jerk control strategy may be particularly sensitive to bias errors on the estimated momentum rate of change $\dot{H}$, which is used as feedback in Eq. (40). The momentum rate of change is estimated as detailed in Sec. V-B. The sources of this bias may be errors in the robot dynamic parameters, or low FT sensor accuracy. To enforce the closed loop system robustness w.r.t. errors in the estimation of $\dot{H}$, we modified Eq. (40) by adding a regularization term as follows:

[TABLE]

where $k_{e}>0$ is a positive scalar, $\xi$ the integral of $\dot{\xi}$ and $\xi_{d}$ is obtained by applying the parametrization inverse mapping on the set of wrenches satisfying the desired momentum rate of change, i.e. $\dot{H}_{d}(f_{d})$. In case of multiple solutions, the one ensuring minimum norm of $f_{d}$ is chosen.

It is important to point out that the regularization term $-k_{e}(\xi-\xi_{d})$ is only necessary in case of errors when estimating the momentum rate-of-change $\dot{H}$. To support this statement, Figure 3 shows the linear and angular momentum error norm during balancing simulations when $7.5\%$ modeling errors are added to the parameters $M,m,C$ and $G$. Note that the robot falls after few seconds if no regularization is added to Eq. (40) (orange line). When the regularization term is added (red line) or the errors on $\dot{H}$ are removed (blue line), stability is retained. In this stability retained case, let us remark that only the estimation of $\dot{H}$ is evaluated correctly, while the model errors are kept present for all the other calculations.

We also carried out robustness tests on the real iCub during a contact switching scenario [22]. The robot starts balancing on two feet, then it switches to balance on the left foot via a finite state machine, and performs highly dynamic movements on the left foot. Finally, the robot returns back to two feet balancing. Results are reported in Table I: the robot succeeded to conclude the demo $60\%$ of times in case of parameters overestimation by $7.5\%$. Dealing with parameters underestimation, instead, seems to be a more challenging task for the controller (43) despite the presence of the regularization term.

When mounted on the robot, the FT sensors accuracy is affected by several phenomena such as temperature, internal stresses and vibrations. More specifically, we observed that even after FT sensor fine calibration the linear forces measurements still have an offset of $\pm 2.5N$, and this offset biases the estimation of $\dot{H}$.

Figure 4 shows the behavior of the linear and angular momentum error norm during several two feet balancing simulations and experiments. On the real iCub, the robot falls after few seconds when $k_{e}=0$, as pointed out by the green line. When adding the regularization term in Eq. (43), stability is retained and the momentum error does not diverge (purple line). On the other hand, the blue line is obtained in simulation with perfect estimation of external forces, and with $k_{e}=0$. In this case, the momentum error does not diverge, thus showing again that the regularization term is not needed when $\dot{H}$ is properly estimated. In simulation, we injected a constant offset of amplitude $2.5N$ to the ”measured” $f_{x}$ component of one of the two contact wrenches. Results correspond to the orange line in Figure 4: stability is no longer retained and the robot falls after few seconds. With the regularization term, the previous balancing performances are restored (red line).

VI-C Comparison with a momentum-based QP controller

We compared the performance of the momentum-based jerk controller (43) with a classical momentum-based QP controller that solves the optimization problem (11) on the real iCub during the contact switching scenario introduced in Section VI-B. Both controllers have been fine tuned for the specific demo. The goal is to show that the momentum-based jerk control guarantees performances that are comparable with a controller already available in the literature. Also, the momentum-based jerk control provides smoother references to the torque controller, as the desired contact wrenches $f^{*}$ are always continuous.

Figure 5 depicts the norm of the left foot (input) contact forces and moments for both the momentum-based jerk control and the momentum-based QP control. Results have been achieved by running 10 experiments for each control strategy. The solid lines represents the average values, while the transparent regions are the variance over the 10 experiments. The orange background represents the instants at which the robot is balancing on two feet, while the white background is when the robot is balancing on the left foot.

The momentum-based jerk controller helps provide smoother references to the torque controller during transitions. In Figure 6, we compared the norms of the linear and angular momentum error. During transition from two feet to left foot balancing, a peak of error is present when the momentum-QP control is used. The peak is caused by the sharp change in the input forces. When jerk control is used (and smoother force input is required), the peak error is reduced by  $90\%$. During highly dynamic movements (white background), jerk and momentum-based QP control show similar tracking performances.

Both controllers use Eq. (42) to generate the input torques. Figure 7 verifies the boundedness of the system’s zero dynamics. In both cases, the zero dynamics does not diverge. Convergence to zero of the joint position error is not necessarily expected as the controllers implement strict task priorities, and the postural task is chosen as the lowest priority task. For further details, a video of the experiment is attached to this paper.

VI-D Disturbance rejection

To evaluate the robustness of the momentum-based jerk controller against unmodeled external force perturbations, we performed the following experiment. The robot balances on two feet and a person pushes and pulls continuously the robot’s upper body. The applied external force is unmodeled, so it is treated as a disturbance by the momentum-based jerk controller. Figure 8 shows the momentum rate of change error norm and the momentum error norm during interaction. Despite the high peaks of errors when the external force is applied, the controller is still able to retain stability. When the force is removed, in fact, the momentum error and its rate of change converge to a value. Exact convergence to zero of the momentum derivative error on the real iCub is difficult to obtain because of the low sensitivity of the FT sensors, so a compromise is achieved by properly tuning the corresponding feedback gains. A video describing the experiment is attached to the paper.

VII Conclusions and Future Work

In this paper, we addressed some common limitations of force and torque controllers for floating base systems based on Quadratic Programming. More specifically, we removed inequality constraints from the optimization problem by designing an invertible, one-to-one mapping that parametrises the contact wrenches into a new set of unconstrained variables. This parametrization guarantees that the contact wrenches always satisfy the contact stability constraints. Based on this mapping, we designed a jerk control framework for floating base systems. We then analyzed a specific use case of the jerk controller, namely a momentum-based jerk control architecture for balancing a 23 DoFs humanoid robot. The controller has been validated both in simulation and on the real iCub, and compared with a classical momentum-based controller. Sensitivity to errors in the momentum rate of change estimation is identified as a drawback of the approach, as it may affect stability and convergence of the closed loop dynamics. A solution for increasing robustness of the controller w.r.t. biases on momentum estimation is presented.

A limitation is that the proposed jerk control architecture does not take into account joint position and torque limits. A future work may involve the integration of these limits in the control framework by extending the approach presented in [28] to the case of floating base robots.

In this paper, the jerk control architecture implements a momentum task and a postural task for balancing. Further development will be done to extend the jerk control to more complex control objectives, as for example humanoid robot walking.

-A Proof of Lemma 1

Proof of 1): when $f^{k}=\phi(\xi^{k})$, the constraint on the vertical force Eq. (6a) is given by:

[TABLE]

which is satisfied for all $\xi_{3}$. Now, substitute the parametrization (14) in the remaining stability constraints (7a)–(10a):

[TABLE]

The vertical force $f_{z}$ is greater than zero and can be removed from (44a), thus leading to the following set of inequalities:

[TABLE]

where also the coefficients $\mu_{c}$ and $\mu_{z}$ have been removed from Eq. (44a) and (47a), respectively. It is straightforward to verify that the constraint (51a) is verified for all finite $\xi_{6}$. Also, direct calculations on Eq. (49a)-(50a) show that:

[TABLE]

Using the facts that $y_{c}^{min}<y_{c}^{max}$, $x_{c}^{min}<x_{c}^{max}$ and $|\tanh(\cdot)|<1$, one easily shows that the above quantities are always satisfied for all $\xi_{4},\xi_{5}$. Concerning the remaining constraint Eq. (48a), note that the argument of the square root can be rearranged as follows:

[TABLE]

So, the constraint (48a) is satisfied if and only if

[TABLE]

Since $|\tanh(\cdot)|<1$, then $\tanh^{4}(\cdot)\leq\tanh^{2}(\cdot)$. As a consequence, one easily finds a minorant of the left hand side of (52), namely

[TABLE]

Since $|\tanh(\cdot)|<1$, then (-A) is (strictly) greater than zero, which renders (52) satisfied.

Proof of 2): First, it follows from (14) that any value of $\xi$ generates a unique $f^{k}\in\mathcal{K^{\prime}}$. Now, assume that $f^{k}\in\mathcal{K^{\prime}}$. Then, we have to show that there exists a unique $\xi=\phi^{-1}(f^{k})$. It is straightforward to compute the inverse mapping of the vertical force and moments parametrization:

[TABLE]

Since $\mathcal{K^{\prime}}\subset\mathcal{K}$, then $f^{k}$ satisfies (6a). Using this fact and the definitions (15a) we conclude that the solutions $(\xi_{3},\xi_{4},\xi_{5},\xi_{6})$ exist. Furthermore, since the above equations are composed of one-to-one correspondences (hyperbolic tangent, logarithm), then these solutions are also unique.

For what concerns the tangential forces $f_{x}$ and $f_{y}$, let us recall the expressions of the $f_{x}$ and $f_{y}$ parametrizations:

[TABLE]

An easy way to compute the inverse mapping is to raise to the square Eq. (54a), which gives:

[TABLE]

In the resulting equations, the square hyperbolic tangents $\tanh^{2}(\xi_{1})$, $\tanh^{2}(\xi_{2})$ appear linearly. Therefore they can be easily computed through matrix inversion:

[TABLE]

Since $f^{k}\in\mathcal{K^{\prime}}$, then the right hand side of (58) is necessarily smaller than one and there exists at least one $\xi$ satisfying (58). Resolving Eq. (58) w.r.t. $\xi_{1},\xi_{2}$ gives two possible solutions, namely $\xi_{1(2)}=\pm\text{atanh}\left(\sqrt{\tanh^{2}(\xi_{1(2)})}\right)$. However, only one of the two solutions satisfies the parametrization (54a): in fact, the terms $f_{z},\mu_{c}$, and $\sqrt{{1+\tanh^{2}(\xi_{1(2)})}}$ on the right-hand side of (54a) are always positive. So, the sign of $\xi_{1}$ ($\xi_{2}$) must correspond to the sign of $f_{x}$ ($f_{y}$), leading to the unique solution:

[TABLE]

Remark: it is possible to verify that if $f_{x}$ and $f_{y}$ belong to $\mathcal{K^{\prime}}$, then the matrix inversion in Eq. (58) can always be performed. In fact, singularities arise when

[TABLE]

Substituting Eq. (56a) in the expression of the determinant allows to verify that the condition $\mu_{c}^{4}f_{z}^{4}=f_{x}^{2}f_{y}^{2}$ never occurs for any $\xi_{1},\xi_{2}$.

Proof of 3): let $\Phi_{k}\in\mathbb{R}^{6\times 6}$ denote the gradient of $f^{k}=\phi(\xi)$. Then, straightforward calculations show that $\Phi_{k}$ is given by:

[TABLE]

Applying the Laplace’s formula for the calculation of the determinant of $\Phi_{k}$ leads to:

[TABLE]

where one has:

[TABLE]

It can be verified that $F_{33},F_{44},F_{55}$ and $F_{66}$ are always different from zero for any finite $\xi_{3},\xi_{4},\xi_{5},\xi_{6}$. We are then left to evaluate the determinant of $\begin{bmatrix}F_{11}&F_{12}\\ F_{21}&F_{22}\end{bmatrix}$, which is $F_{11}F_{22}-F_{12}F_{21}$. After noting that the product $F_{11}F_{22}$ is contained in the expression of $F_{12}F_{21}$, one is left with:

[TABLE]

which is non-zero for any finite $\xi_{1},\xi_{2},\xi_{3}$. Hence, matrix $\Phi_{k}$ is invertible for any finite $\xi$.

-B Proof of Lemma 2

Proof of 1): first, observe that Lyapunov stability of the equilibrium point $(\tilde{H},\dot{\tilde{H}})=(0,0)$ follows from (32) with $\ddot{H}^{*}$ given by (31). By definition, Lyapunov stability implies that there exist initial conditions $(\tilde{H},\dot{\tilde{H}})(0)$ such that the system trajectories $(\tilde{H},\dot{\tilde{H}})(t)$ remain as close as we wish to the point $(0,0)$. We are then left to show that the variable $\xi$ remains bounded while the system trajectories evolve close to the equilibrium point: the boundedness of $\xi$ would imply that the matrix $A\Phi(\xi)$ in (30) remains of full rank, thus ensuring that the equality (32) can be satisfied $\forall t$.

Now, since the desired momentum $H_{d}$ is a feasible system equilibrium, then there exists a bounded $\xi_{e}(t)$ such that

[TABLE]

Furthermore, the behaviour of system

[TABLE]

close to the equlibrium point can be obtained via linearisation techniques, which in view of (59) yields

[TABLE]

with $\tilde{\xi}=\xi-\xi_{e}$. Since it is assumed that the robot makes a single contact with the environment, namely $n_{c}=1$, then the matrix $A\in\mathbb{R}^{6\times 6}$. From (29a) one easily verifies that the matrix $A$ is always invertible, with $A^{-1}$ a bounded matrix. Then, one has

[TABLE]

with $c_{a}>0$. As a consequence of the Lyapunov stability of $(\tilde{H},\dot{\tilde{H}})=(0,0)$, there exist initial conditions $(\tilde{H},\dot{\tilde{H}})(0)$ for which the variable $\tilde{\xi}$ stays bounded and as close as we wish to zero, thus rendering the linearisation (60) consistent $\forall\ t$. Finally, since $\xi_{e}$ is bounded, then $\xi$ is also bounded, and the equality (32) can be satisfied $\forall t$. This latter fact, along with $\ddot{H}^{*}$ given by (31), then implies local asymptotic stability since (32) is satisfied $\forall\ t$.

Proof of 2): since it is assumed that the variable $\xi$ is bounded $\forall t$, then (32) can be satisfied $\forall\ t$ because the matrix $A\Phi(\xi)$ in (30) remains of full rank. Then, the closed loop dynamics is $\ddot{\tilde{H}}=-K_{d}\dot{\tilde{H}}-K_{p}{\tilde{H}}\ \forall t$, which implies global asymptotic stability of the equilibrium point $(\tilde{H},\dot{\tilde{H}})=(0,0)$.

-C Proof of Lemma 3

Proof of 1): define the following Lyapunov function candidate:

[TABLE]

Note that $V=0\iff(I,\tilde{H},\zeta)=(0,0,0)$. Compute the Lyapunov function derivative $\dot{V}$:

[TABLE]

It is easy to prove that $\dot{V}\leq 0$ when $\dot{\zeta}+K_{o}^{-1}\tilde{H}=-\zeta$. Then, in view of $\dot{\zeta}$ given by (39a) and the definition of $\zeta$ one has:

[TABLE]

and after a rearrangement, Eq. (62) leads to the definition of the control input $\dot{\xi}^{*}$ as in Eq. (40), which gives

[TABLE]

This result implies stability of the equilibrium point and the boundedness of system’s trajectories. Furthermore, as long as Eq. (62) holds, the closed loop dynamics is given by

[TABLE]

and Eq. (37a)-(38a). The system is therefore autonomous, and the convergence of $\tilde{H},\zeta$ and $\dot{\tilde{H}},\dot{\zeta}$ to zero can be proved via the La Salle’s theorem. Convergence to zero of $I$ can be proven by computing Eq. (38a) on the invariant manifold.

Analogously to the proof of Lemma 2, we are left to show that in a neighborhood of the equilibrium point, the variable $\xi$ remains bounded, which guarantees that (62) holds. Note that since the desired momentum $H_{d}$ is a feasible equilibrium, then there exists a bounded $\xi_{e}(t)$ such that

[TABLE]

Consider now the right hand side of (38a), namely

[TABLE]

and linearize $Af-mge_{3}-\dot{H}_{d}$ close to the equilibrium trajectory $\xi_{e}$, which yields

[TABLE]

with $\bar{K}=(-K_{p},\ -K_{d},\ 1_{6})$. Since it is assumed that the robot makes a single contact with the environment, namely $n_{c}=1$, then the matrix $A\in\mathbb{R}^{6\times 6}$. From (29a) one easily verifies that the matrix $A$ is always invertible, with $A^{-1}$ a bounded matrix. Then, one has

[TABLE]

with $c_{b}>0$. As a consequence of the Lyapunov stability of $(I,\tilde{H},\zeta)=(0,0,0)$, there exist initial conditions $(I,\tilde{H},\zeta)(0)$ for which the variable $\tilde{\xi}$ stays bounded and as close as we wish to zero, thus rendering the linearisation (65) consistent $\forall\ t$. Finally, since $\xi_{e}$ is bounded, then $\xi$ is also bounded, and the equality (62) can be satisfied $\forall t$.

Proof of 2): it is assumed that the $\xi$ is bounded $\forall t$; then (62) is satisfied $\forall\ t$ because the matrix $A\Phi(\xi)$ in (62) remains of full rank. In view of the radial unboundedness of (61), and of the expression of (63) that is valid $\forall t$, then global asymptotic stability follows from the same La Salle arguments above.

-D Computation of $\dot{\xi}_{0}$:

The calculation of $\dot{\xi}_{0}$ is carried out by substituting (40) in the expression of (42). Fist, rewrite (42) as:

[TABLE]

Then, consider the following Lyapunov function:

[TABLE]

Substitute now $\dot{\xi}$ given by (40) into $\dot{V}$, which leads to:

[TABLE]

$\dot{\xi}_{1}{\mathrel{\mathop{\mathchar 58\relax}}=}{(A\Phi)}^{\dagger}\,[\ddot{H}_{d}{-}(K_{d}{+}1_{6})\dot{\tilde{H}}{-}(K_{d}{+}K^{-1}_{o}{+}K_{p})\tilde{H}{-}K_{p}I{-}\dot{A}f]$. A solution that minimizes the joint torque norm is to impose:

[TABLE]

with $K_{\tau}$ a symmetric and positive definite matrix. When the equivalence (68) is satisfied, the Lyapunov derivative Eq. (67) becomes $\dot{V}=-\tau^{\top}K_{\tau}\tau\leq 0$ and the input joint torques converge to zero. However, this is not the case as the rank of the matrix $(\Theta\Phi N_{A\Phi})$ that multiplies the free variable $\dot{\xi}_{0}$ is lower than the dimension of the joint torques vector $\tau\in\mathbb{R}^{n}$. Nevertheless, we compute the closest solution to Eq. (68), that leads to the following expression of $\dot{\xi}_{0}$:

[TABLE]