# Equilibria, Stability, and Sensitivity for the Aerial Suspended Beam Robotic System Subject to Parameter Uncertainty

Chiara Gabellieri, Marco Tognon, Dario Sanalitro, Antonio Franchi

► To cite this version:

Chiara Gabellieri, Marco Tognon, Dario Sanalitro, Antonio Franchi. Equilibria, Stability, and Sensitivity for the Aerial Suspended Beam Robotic System Subject to Parameter Uncertainty. IEEE Transactions on Robotics, 2023, 39 (5), pp.3977 - 3993. 10.1109/TRO.2023.3279033. hal-04161342

HAL Id: hal-04161342

https://laas.hal.science/hal-04161342

Submitted on 13 Jul 2023

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.

Public Domain
# Equilibria, Stability, and Sensitivity for the Aerial Suspended Beam Robotic System subject to Parameter Uncertainty

C. Gabellieri¹ IEEE Member, M. Tognon² IEEE Member, D. Sanalitro³,⁴, A. Franchi¹,³,⁵, IEEE Fellow

**Abstract**—This work studies how parametric uncertainties affect the cooperative manipulation of a cable-suspended beam-shaped load by means of two aerial robots not explicitly communicating with each other. In particular, the work sheds light on the impact of the uncertain knowledge of the model parameters available to an established communication-less force-based controller. First, we find the closed-loop equilibrium configurations in the presence of the aforementioned uncertainties, and then we study their stability. Hence, we show the fundamental role played in the robustness of the load attitude control by the internal force induced in the manipulated object by non-vertical cables. Furthermore, we formally study the sensitivity of the attitude error to such parametric variations, and we provide a method to act on the load position error in the presence of the uncertainties. Eventually, we validate the results through an extensive set of numerical tests in a realistic simulation environment including underactuated aerial vehicles and sagging-prone cables, and through hardware experiments.

**Index Terms**—Cooperative Aerial Manipulation; Aerial Systems: Mechanics and Control; Motion Control

## I. INTRODUCTION

It is nowadays acknowledged that the interest in Unmanned Aerial Vehicles (UAVs) is becoming wider and wider by virtue of their ability to embrace an ample set of applications. A very recent and popular topic in aerial robotics is physical interaction using aerial manipulators [1], [2], [3] for applications such as contact-based inspection, assembly, human assistance, etc. To solve these challenges, aerial platforms are endowed with physical interaction tools, such as cables [4] or more complex robotic arms [5].

Researchers have considered taking advantage of the co-operation between multiple robots to enhance the overall payload and manipulate large objects [6], [7], [8], [9]. Different methods have been developed to tackle multi-robot aerial manipulation. In [10] and [11] the authors use passive manipulation tools to solve the cooperative aerial transportation of rigid and elastic objects, respectively. Multiple flying arms are instead used in [12], [13]. Cables have been often

considered in multi-robot manipulation scenarios, because, in addition to being lightweight and low-cost, they also mitigate the coupling between the system dynamics and the robots' attitude, which can simplify the control problem, especially when using underactuated aerial platforms.

### A. Related Works

The problem of manipulating a cable-suspended load through a team of aerial vehicles has been studied, e.g., in [14], [15], [16], [17], [18], [19]. In [20], a robust pose controller for a cable-suspended load manipulated by multiple UAVs is presented. Stability is ensured through gain tuning given a bound on the uncertainties affecting the kinematic parameters. Formation control to transport the payload with a focus on robustness is described in [21], where modeling uncertainties are also taken into account.

A standard system that has attracted substantial interest in the research community is composed of two aerial vehicles manipulating a beam-like load through cables [22], [23], [24], [17], [25], [26], [27]. Such standard configuration is of interest for several real-world applications, especially in the construction field, where we find columns, wooden pillars, iron beams for cement walls, scaffolds, pipes, pieces of roofs, and other beam-like building elements. Two is the minimum number of aerial robots allowing to control both the position and attitude of a cable-suspended beam-like load [23]. While three aerial robots allow controlling the entire pose of a generic rigid body [28], using more than two robots for a beam-like load it is arguably not the optimal solution in most of the cases because of the increased complexity of the system without being necessary for the control of the load.

In [23] and [29], the authors propose a method for the transportation of a cable-suspended beam load by two aerial vehicles that have access to the state of the load; they consider rigid and elastic cables, respectively. In [24], centralized and decentralized model predictive control is proposed for a system of two UAVs manipulating a beam load through cables.

Decentralized algorithms as [30] are more robust and scalable with respect to (w.r.t.) the number of robots. However, decentralized communication-less approaches have been also intensively studied in the literature [31] because communication delays and packet losses are among the principal causes undermining the performance and stability of the system in real implementations, and because the hardware and software complexity can be reduced by confining explicit communication. In [26], a method relying on visual feedback is presented.

This work has been partially funded by the European Union's Horizon Europe research and innovation program [ID: 101059875] Flyflic and Horizon 2020 research and innovation programme [ID 871479] AERIAL-CORE.

¹Robotics and Mechatronics group, EEMCS faculty, University of Twente, Enschede, The Netherlands ²Univ Rennes, CNRS, Inria, IRISA, F-35000 Rennes, France ³LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France. ⁴Electrical electronic and computer engineering department

at University of Catania, 95125, Catania, Italy ⁵Department of Computer, Control and Management Engineering, Sapienza University of Rome, 00185 Rome, Italy. Email: c.gabellieri@utwente.nl, a.franchi@utwente.nl
As an alternative to vision, a force-based method that uses admittance controllers and a leader-follower scheme is typically used to address communication-less aerial manipulation of cable-suspended objects [25], [32], [33], [34]. The leader robot guides the system following a predetermined trajectory, while the second robot, which carries a portion of the load weight, follows its lead by sensing the cable force variations.

A primary goal of [25] is to keep the cables always vertical during transportation, meaning that no internal force is induced in the object. The authors in [32] extend the results of [25] towards the N-robot case and provides a method for tuning the gains of the admittance controllers to improve the robustness against unmodeled dynamics and parametric uncertainties. In both works experiments are shown in which, however, the altitude of the robots is set to a predetermined reference, implying either a centralized vertical movement coordination or restricted vertical motion.

For such a popular class of communication-less, admittance-controlled, and leader-follower schemes, the formal analysis of the closed-loop system equilibrium configurations and their stability was presented for the first time in our previous work [33]. There, we showed that inducing an internal force on the load through non-vertical cables is required for full-pose regulation, especially to prevent arbitrary vertical movements of the robots that would interfere with the regulation of the load pitch and center position. In [34], we considered N robots, empirically showing through extensive simulations the effect of changing the number of leader robots on the stability and robustness against disturbances. Both works tackle only the ideal case where perfect knowledge of the system parameters is available to the admittance controllers of each aerial robot.

Despite being of primary interest, it has been unclear until now if and how in-practice-unavoidable uncertainties impact the pose regulation in the aforementioned control framework. Such a gap is filled in this work by introducing uncertainties on those system parameters used in the control action.

Adaptive control laws have been proposed in the literature for the system in question, however, they are based on different assumptions than those used in this work. For instance, the full state has been considered available for feedback in [31], or the robots rigidly attached to the object [35], [36]. Other works assume the load mass is the sole uncertain parameter and only focus on the translational velocity regulation [13], or rely on a communication network [37].

In this work, we have found that the internal force induced by non-vertical cables plays a fundamental role in enabling task execution, especially in realistic conditions characterized by uncertainties. The importance of this is masked when vertical movements of the robots are prevented or the leader-follower approach is used solely to regulate the load motion in the horizontal plane, as in [25], [32]. On the other hand, the role of the internal force is crucial if the admittance-based communication-less approach is applied in the full 3D space and, hence, communication-less full-pose regulation is sought.

While some loads may be damaged by internal forces, this is easily prevented in practice by enclosing the loads in suitable cases. Also, internal forces require additional control effort, which is justified by the benefits in terms of convergence and

robustness of the load pose control, as will be clear in the following. Indeed, internal forces have been often proposed also in the robotic grasping literature as a tool to make the grip on the object robust thanks to friction [38].

## B. Contributions of the Work and Outline of the Paper

The contribution of this work is showing the effects of parametric uncertainties on the static regulation of the load pose when the usual approach [25], [32], [33], [34] based on admittance-controlled leader-follower aerial robots is used for manipulating a cable-suspended beam in the absence of explicit communication. In this approach, each robot knows only its own state and the force in its cable, retrievable from the robot's state using an external force observer.

Note that, unlike in [29], it is not feasible to assume all robots have knowledge of the object state. This is because the object state is based on the state of all robots, but data exchange among them is not considered in our scenario.

Throughout the manuscript, we show that it is best to avoid the intuitive idea of having the cables vertical when performing the manipulation with force-based methods in real scenarios, i.e., when uncertainties are present. We address the problem aiming at a mathematically-sound point of view. We point out that we restrict the analysis to pose regulation, hence, to quasi-static motion of the load. Tracking of more aggressive trajectories is left to future work. With the above context in mind, the key contributions can be succinctly summarized:

* after formally studying the equilibrium points of the closed-loop system in the presence of uncertain parameters, their stability is proved using Lyapunov's theory;

* the impact of an internal force induced by non-vertical cables on the robustness of the load pose control is formally studied;

* the effect of the internal force of diminishing the sensitivity of the load attitude error to parametric uncertainty or variations is shown;

* a method for correcting the load position inaccuracy induced by the uncertainties is also presented;

* we present extensive numerical results and hardware experiments supporting the claims conveyed by the theoretical analysis;

* last but not least, this work generalizes the system model by considering a generic position of the center of mass (CoM) of the load rather than assuming it to be in the middle of the cable anchoring points as done in [33].

The paper is organized as follows. Sec II contains some background useful to better understand the results of the work. In Secs III to V, we present the three main contributions of the work: Sec III contains the derivation of the equilibrium points, and Sec IV their stability analysis; Sec V highlights the role of the internal forces in the load error robustness and sensitivity to parametric variations. The results of the simulations and experiments are presented in Sec VI and Sec VII, respectively. Conclusive discussions are in Sec VIII.

## II. BACKGROUND

In this section, we provide the background needed to understand the contribution of the work. Specifically, we quickly
<table><caption>TABLE I: Notation—General symbols and reference frames</caption><tr><td><b>I</b><sub>i</sub></td><td>i × i identity matrix</td></tr><tr><td><b>e</b><sub>i</sub></td><td>i-th column of <b>I</b><sub>3</sub></td></tr><tr><td><b>S</b>(★)</td><td>skew operator</td></tr><tr><td><b>Ker</b>(★)</td><td>nullspace of ★</td></tr><tr><td><b>diag</b>(★)</td><td>diagonal matrix</td></tr><tr><td><b>★</b><sup>⊥</sup></td><td>transpose of ★</td></tr><tr><td><b>||★||</b></td><td>2-norm of ★</td></tr><tr><td><b>̂★</b></td><td>desired value of ★</td></tr><tr><td><b>★</b><sup>Eq</sup></td><td>value of ★ at the equilibrium</td></tr><tr><td><b>★</b><sup>̂</sup></td><td>uncertain value of ★</td></tr><tr><td><b>Δ</b><sub>*</sub></td><td>★ - ★̂</td></tr><tr><td><b>F</b><sub>W</sub></td><td>inertial reference frame</td></tr><tr><td><b>F</b><sub>L</sub></td><td>load reference frame</td></tr><tr><td><b>F</b><sub>Ri</sub></td><td>i-th robot reference frame</td></tr><tr><td>{<b>O</b><sub>A</sub>, <b>x</b><sub>A</sub>, <b>y</b><sub>A</sub>, <b>z</b><sub>A</sub>}</td><td>Origin, X-, Y-, and Z- axis of <b>F</b><sub>A</sub></td></tr><tr><td><b>L</b><sub>*</sub></td><td>★ expressed in <b>F</b><sub>L</sub></td></tr><tr><td><b>★</b></td><td>time derivative of ★</td></tr></table>

<table><caption>TABLE II: Notation—System variables</caption><tr><td><b>g</b></td><td>gravity acceleration</td></tr><tr><td><b>p</b><sub>L</sub>, <b>ṗ</b><sub>L</sub>, <b>ṗ̃</sub><sub>L</sub></td><td>load position, velocity, acceleration</td></tr><tr><td><b>R</b><sub>L</sub></td><td>rotation of <b>F</b><sub>L</sub> w.r.t. <b>F</b><sub>W</sub></td></tr><tr><td><b>ω</b><sub>L</sub></td><td>load angular velocity</td></tr><tr><td><b>q</b><sub>L</sub></td><td>set composed of (<b>p</b><sub>L</sub>, <b>R</b><sub>L</sub>)</td></tr><tr><td><b>v</b><sub>L</sub></td><td>[<b>p</b><sub>L</sub><sup>⊥</sup><b>ω</b><sub>L</sub><sup>⊥</sup>]<sup>⊥</sup></td></tr><tr><td><b>ψ</b>, θ</td><td>yaw and pitch of the load</td></tr><tr><td><b>B</b><sub>i</sub></td><td>i-th cable attaching point on load</td></tr><tr><td><b>b</b><sub>1</sub></td><td>position of point <b>B</b><sub>i</sub></td></tr><tr><td><b>m</b><sub>L</sub>, <b>J</b><sub>L</sub></td><td>load mass and rotational inertia</td></tr><tr><td><b>L</b></td><td>load length, with <i>l</i> = 1/L</td></tr><tr><td><b>M</b><sub>L</sub></td><td>diag(<b>m</b><sub>L</sub><b>I</b><sub>3</sub>, <b>J</b><sub>L</sub>)</td></tr><tr><td><b>g</b><sub>L</sub></td><td>gravity terms in load dynamics</td></tr><tr><td><b>c</b><sub>L</sub></td><td>Coriolis terms in load dynamics</td></tr><tr><td><b>G</b></td><td>load grasp matrix</td></tr><tr><td><b>p</b><sub>Ri</sub>, <b>ṗ</b><sub>Ri</sub>, <b>ṗ̃</sub><sub>Ri</sub></td><td>i-th robot position, velocity, acceleration</td></tr><tr><td><b>R</b><sub>Ri</sub></td><td>rotation of <b>F</b><sub>Ri</sub> w.r.t. <b>F</b><sub>W</sub></td></tr><tr><td><b>q</b><sub>R</sub></td><td>[<b>p</b><sub>Ri</sub><sup>⊥</sup> <b>p</b><sub>R2</sub><sup>⊥</sup>]<sup>⊥</sup></td></tr><tr><td><b>k</b><sub>i</sub></td><td>stiffness of the i-th cable</td></tr><tr><td><b>l</b><sub>0i</sub></td><td>rest length of the i-th cable</td></tr><tr><td><b>l</b><sub>i</sub></td><td><b>p</b><sub>Ri</sub> = <b>b</b><sub>i</sub>, vector along the i-th cable</td></tr><tr><td><b>f</b><sub>i</sub></td><td>force of the i-th cable on the load</td></tr><tr><td><b>f</b></td><td>[<b>f</b><sub>1</sub><sup>⊥</sup> <b>f</b><sub>2</sub><sup>⊥</sup>]<sup>⊥</sup></td></tr><tr><td><b>q</b></td><td>set (<b>q</b><sub>R</sub>, <b>q</b><sub>R</sub>)</td></tr><tr><td><b>v</b></td><td>[<b>v</b><sub>R</sub><sup>⊥</sup> <b>v</b><sub>L</sub><sup>⊥</sup>]<sup>⊥</sup></td></tr><tr><td><b>x</b></td><td>(<b>q</b>, <b>v</b>)</td></tr><tr><td><b>u</b><sub>Ri</sub></td><td>control input of robot i-th</td></tr><tr><td><b>M</b><sub>Ai</sub></td><td>i-th robot control gain: apparent inertia, <b>M</b><sub>A</sub> diag(<b>M</b><sub>A1</sub>, <b>M</b><sub>A2</sub>)</td></tr><tr><td><b>B</b><sub>Ai</sub></td><td>i-th robot control gain: apparent damping, <b>B</b><sub>A</sub> diag(<b>B</b><sub>A1</sub>, <b>B</b><sub>A2</sub>)</td></tr><tr><td><b>K</b><sub>Ai</sub></td><td>i-th robot control gain: virtual spring stiffness, <b>K</b><sub>A</sub> diag(<b>K</b><sub>A1</sub>, <b>K</b><sub>A2</sub>)</td></tr><tr><td><b>π</b><sub>Ai</sub></td><td>feedforward term of i-th robot's control, <b>π</b><sub>A</sub>=[<b>π</b><sub>A1</sub> <b>π</b><sub>A2</sub>]<sup>⊥</sup></td></tr><tr><td><b>t</b><sub>L</sub></td><td>load internal force</td></tr><tr><td><b>ξ</b></td><td>(<b>b</b><sub>1</sub><b>m</b><sub>L</sub> - <b>b</b><sub>1</sub><b>ṁ</b><sub>L</sub>)</td></tr><tr><td><b>m</b>(<b>q</b>, <b>v</b>, <b>π</b><sub>A</sub>)</td><td>closed-loop dynamics</td></tr><tr><td><b>Q</b>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>)</td><td>{<b>q</b> satisfying Theorem 2}</td></tr><tr><td><b>Q</b><sup>+</sup>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>)</td><td><b>Q</b>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>) where <b>R</b><sub>L</sub><sup>eq</sup><b>e</b><sub>1</sub> has same sign of (17)</td></tr><tr><td><b>R</b><sub>L</sub><sup>eq+</sup></td><td><b>R</b><sub>L</sub><sup>eq</sup> in <b>q</b> ∈ <b>Q</b><sup>+</sup>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>)</td></tr><tr><td><b>R</b><sub>L</sub><sup>eq-</sup></td><td><b>R</b><sub>L</sub><sup>eq+</sup> <b>R</b><sub>zL</sub>(<b>π</b>)</td></tr><tr><td><b>Q</b><sup>-</sup>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>)</td><td><b>q</b> ∈ <b>Q</b>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>) s.t. <b>R</b><sub>L</sub><sup>eq</sup> = <b>R</b><sub>L</sub><sup>eq-</sup></td></tr><tr><td><b>Q</b><sub>1</sub>(0, <b>q̃</b><sub>L</sub>)</td><td><b>q</b> ∈ <b>Q</b>(0, <b>q̃</b><sub>L</sub>) s.t. (<b>R</b><sub>L</sub><sup>eq</sup><b>e</b><sub>1</sub>)<sup>⊥</sup><b>e</b><sub>3</sub> = +1</td></tr><tr><td><b>Q</b><sub>2</sub>(0, <b>q̃</b><sub>L</sub>)</td><td><b>q</b> ∈ <b>Q</b>(0, <b>q̃</b><sub>L</sub>) s.t. (<b>R</b><sub>L</sub><sup>eq</sup><b>e</b><sub>1</sub>)<sup>⊥</sup><b>e</b><sub>3</sub> = -1</td></tr><tr><td><b>Q</b>(0, <b>q̃</b><sub>L</sub>)</td><td>{<b>x</b> : <b>q</b> ∈ <b>Q</b>(0, <b>q̃</b><sub>L</sub>), <b>v</b> = 0}</td></tr><tr><td><b>Q</b><sub>i</sub>(0, <b>q̃</b><sub>L</sub>)</td><td>{<b>x</b> : <b>q</b> ∈ <b>Q</b>(0, <b>q̃</b><sub>L</sub>), <b>v</b> = 0}</td></tr><tr><td><b>Q</b>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>)<sup>±</sup></td><td>{<b>x</b> : <b>q</b> ∈ <b>Q</b>(<b>t</b><sub>L</sub>, <b>q̃</b><sub>L</sub>)<sup>±</sup>, <b>v</b> = 0}</td></tr><tr><td><b>e</b><sub>R</sub><sub>L</sub>, <b>e</b><sub>pL</sub></td><td>load attitude and position errors</td></tr></table>

recall the system's main variables, the dynamics equations, and the already established findings. For the sake of readability, the notation is also summarized in Tables I and II.

The considered system, sketched in Fig. 1, is the typical rigid beam load attached to two aerial robots through cables.

1) Load Model: The beam-like load has mass $m_L \in \mathbb{R}_{>0}$ and positive-definite rotational inertia $J_L \in \mathbb{R}^{3 \times 3}$. The frame $\mathcal{F}_L = \{O_L, x_L, y_L, z_L\}$, where $O_L$ coincides with the load CoM, is rigidly attached to the load. The inertial frame is denoted with $\mathcal{F}_W = \{O_W, x_W, y_W, z_W\}$ where $z_W$ is oriented in the direction opposite to the gravity. The position and orientation of $\mathcal{F}_L$ w.r.t. $\mathcal{F}_W$, defined by the vector$^1 W p_L \in \mathbb{R}^3$ and the rotation matrix $R_L$, respectively, describe the full configuration of the load. We recall that the rotation along the axis that passes between the two cable anchoring points is not controllable by the robots. Being the load a beam, only the yaw angle, $\psi$, and pitch angle, $\theta$, are used to describe its attitude. The usual equations of a rigid body subject to gravity and contact forces describe the dynamics of the load as

$$ \dot{v}_L = M_L^{-1} (-c_L(v_L) - g_L + G(q_L)f), \quad (1) $$

where $q_L = (p_L, R_L)$; $v_L = [\hat{p}_L^⊥ \hat{\omega}_L^⊥]^⊥$ with $^⊥ \omega_L \in \mathbb{R}^3$ the angular velocity of $\mathcal{F}_L$ w.r.t. $\mathcal{F}_W$ expressed in $\mathcal{F}_L$; $M_L = \text{diag}(m_L I_3, J_L)$ with $I_3 \in \mathbb{R}^{3 \times 3}$ the identity matrix; $g_L = [m_L g_e^⊥ 0]^⊥$, where g is the gravitational acceleration and $e_i$ is the canonical unit vector with a 1 in the i-th entry. Coriolis and centrifugal terms are given by

$$ c_L = \begin{bmatrix} 0 \\ S(\omega_L)J_L\omega_L \end{bmatrix} $$

where $S(★)$ is the skew operator$^2$, and the grasp matrix is

$$ G = \begin{bmatrix} I_3 & I_3 \\ S(^L b_1) R_L^⊥ & S(^L b_2) R_L^⊥ \end{bmatrix}. $$

The load is suspended by two cables from two anchoring points, $B_i$ with $i=1,2$, for which the position w.r.t. $\mathcal{F}_L$ is described by the vector $^L b_i \in \mathbb{R}^3$. Each cable exerts on the load a force $f_i$ such that $f = [f_1^⊥ f_2^⊥]^⊥$ in (1). By simple kinematics, the position of $B_i$ w.r.t. $\mathcal{F}_W$ is then given by $b_i = p_L + R_L^⊥ b_i$. Since we are considering a beam-like load, the object CoM is aligned with the two anchoring points of the cables. Without loss of generality, we assume that $^L b_1 = [b_1 \ 0 \ 0]^⊥$ and $^L b_2 = [-b_2 \ 0 \ 0]^⊥$, where $b_i \in \mathbb{R}_{>0}$, for $i=1,2$. We also define the beam's length $L = b_1 + b_2$.

2) Robot Model: We define a frame $\mathcal{F}_{Ri} = \{O_{Ri}, x_{Ri}, y_{Ri}, z_{Ri}\}$ rigidly attached to the i-th robot and centered in its CoM. The i-th cable is attached to the i-th aerial vehicle at the point $O_{Ri}$, which allows decoupling the robot's attitude dynamics from the rest [25], [29]. $\mathcal{F}_{Ri}$ is used to describe the position and rotation of the vehicle w.r.t. $\mathcal{F}_W$, denoted by the vector $p_{Ri} \in \mathbb{R}^3$, and the rotation matrix $R_{Ri} \in SO(3)$, respectively.

The use of recent controllers for unidirectional- and multidirectional-thrust vehicles [39], [40] and disturbance observers for aerial vehicles has been experimentally proven to

<footnote>$^1$The left superscript indicates the reference frame. From now on, $\mathcal{F}_W$ is considered as a reference frame when the superscript is omitted.</footnote>

<footnote>$^2$Given $x \in \mathbb{R}^3$, $S(x) \in \mathbb{R}^{3 \times 3}$ is such that $S(x)y = x \times y$ for all $y \in \mathbb{R}^3$.</footnote>
<figure><img src="image_1.png" /><figcaption>Fig. 1: Representation of the system and its major variables. The two aerial vehicles do not need to be necessarily quadrotors since the analysis and control design is valid for general aerial vehicles.</figcaption></figure>

result in negligible tracking errors even in the presence of external disturbances.

Consequently, due to the time-scale separation between the fast attitude dynamics and the slow translational dynamics [41], the closed-loop translational dynamics of the robot under the influence of the position controller effectively behaves like that of a double integrator $\ddot{p}_{Ri} = u_{Ri}$, where $u_{Ri}$ is a virtual input. In other words, it is safe to assume that the aerial robots together with a sufficiently accurate position controller can track any desired $C^2$ trajectory with negligible error in the domain of interest [42], independently from external disturbances [33]. In this work, we follow such experimentally validated common practice for the theoretical derivations contained in Secs III, IV, and V (see, e.g., the experiments on cooperative load transport in [25], [32]), and we utilize underactuated quadrotors for both the numerical and experimental validations in Sec VI and Sec VII.

### 3) Cable Model
Cable-to-robot and cable-to-load connections are modeled as passive and mass-negligible rotational joints. Besides, the i-th cable is represented as a unilateral spring along its principal direction, which is a frequently adopted model [43], [33], [44], [45]. As commonly done in the state of the art, the cables' mass and inertia are assumed negligible in comparison to the robots' and load's. Its parameters are the constant elastic coefficient $k_i \in \mathbb{R}_{>0}$ and the constant rest length denoted by $l_{0i}$.

The attitude of the i-th cable w.r.t. $\mathcal{F}_W$ is expressed by the normalized vector $^3l_i/||l_i|| \in \mathbb{S}^2$, where $l_i = p_{Ri} - b_i$. The force acting on the load at $B_i$, given a certain length $||l_i||$ of the cable, is given by the simplified Hooke's law:

$$f_i = \|f_i\| \frac{l_i}{\|l_i\|}, \quad \|f_i\| = \begin{cases} k_i(\|l_i\| - l_{0i}) & \text{if } \|l_i\| - l_{0i} > 0 \\ 0 & \text{otherwise} \end{cases}. \qquad (2)$$

The force produced on the other hand of the cable, i.e., on the i-th robot at $O_{Ri}$, is equal to $-f_i$.

### 4) Controller
We recall that to regulate the pose of the manipulated load to a desired configuration $\bar{q}_L = (\bar{p}_L, \bar{R}_L)$, an admittance controller is used on the robots [33]:

$$u_{Ri} = M_{Ai}^{-1} (-B_{Ai} \dot{p}_{Ri} - K_{Ai} p_{Ri} - f_i + \pi_{Ai}), \quad (3)$$

where the positive-definite symmetric matrices $M_{Ai}, B_{Ai}, K_{Ai} \in \mathbb{R}^{3 \times 3}$ are, respectively, the virtual inertia

$$^3\mathbb{S}^2 = \{v \in \mathbb{R}^3 \mid \|v\| = 1\}$$

<figure><img src="image_2.png" /><figcaption>Fig. 2: Schematic representation of the overall system including both physical and control blocks. The input of the admittance controller is affected by the uncertainty of the system parameters.</figcaption></figure>

of the robot, and the damping and stiffness coefficients of a virtual spring-damper system that links the robot and a desired reference frame; $\pi_{Ai} \in \mathbb{R}^3$ is an additional forcing input that is properly set to steer the load to the desired configuration.

Remark 1. One can notice that (3) requires only local information, i.e. the robot's state ($p_{Ri}, \dot{p}_{Ri}$), which can be retrieved with standard onboard sensors like IMU, GPS, and cameras; the force applied by the cable $f_i$, which can be directly measured by an onboard force sensor or estimated by a sufficiently precise model-based observer as done in [46], [25]. Therefore, the described method is decentralized and does not require explicit communication between the robots.

### 5) Closed-loop Model
From equations (1) and (3), the closed-loop system dynamics can be written as $\dot{v} = m(q, v, \pi_A)$ where

$$m(q, v, \pi_A) = \begin{bmatrix} M_A^{-1}(-B_A v_R - K_A q_R - f + \pi_A) \\ M_L^{-1}(-c_L(v_L) - g_L + Gf) \end{bmatrix}, \quad (4)$$

with $q_R = [\dot{p}_{R1} \ \dot{p}_{R2}]^\top$, $q = (q_R, q_L)$, $v_R = [\dot{p}_{R1}^\top \ \dot{p}_{R2}^\top]^\top$ and $v = [v_R^\top \ v_L^\top]^\top$; $\pi_A = [\pi_{A1}^\top \ \pi_{A2}^\top]^\top$. Furthermore $M_A = \text{diag}(M_{A1}, M_{A2})$, $B_A = \text{diag}(B_{A1}, B_{A2})$ and $K_A = \text{diag}(K_{A1}, K_{A2})$.

In order to coordinate the motion of the robots in a decentralized way, a leader-follower approach is used. In this way, only the designated leader will have active control over the position of the load. On the other hand, the other robot will follow, partially sustaining the weight of the load and contributing to the control of the load attitude. Choosing without loss of generality, robot 1 as the leader and robot 2 as the follower, the leader-follower approach is achieved as previously proposed in [25], [32], [33], [34] by setting $K_{A1} \neq 0$ and $K_{A2} = 0$.

In the following, we present Theorem 1, from [33], along with two definitions that will help the reader comprehend the contribution of this work.

**Definition 1** (Equilibrium configuration). $q$ is an equilibrium configuration, indicated as $\bar{q}$, if $\exists \pi_A \text{ s.t. } 0 = m(q, 0, \pi_A)$, i.e, if the corresponding zero-velocity state is a forced equilibrium for the system (4) for a certain forcing input $\pi_A$.

**Definition 2** (Load internal force). For the considered system, the load internal force is defined as

$$t_L := \frac{1}{2} \dot{f}^\top [I_3 - I_3]^\top R_L e_1, \quad (5)$$
where $[\mathbf{I}_3 - \mathbf{I}_3]^\top \mathbf{R}_L e_1 \in \text{Ker}(\mathbf{G})$. We have that

* if $t_L > 0$ the internal force causes a tension in the load;

* if $t_L < 0$ the internal force causes a compression.

The following result, proven in [33], provides the expression of the forcing input $\pi_A$ and the robot configurations $q_R$ for which, given a desired load configuration $\bar{q}_L$, $q = (q_R, \bar{q}_L)$ is an equilibrium configuration of the system.

**Theorem 1** (equilibrium inverse problem, provided in [33], reported here for completeness). Consider the closed-loop system (4) and assume that the load is at a given desired configuration $\bar{q}_L = (\bar{p}_L, \bar{R}_L)$. For each internal force $t_L \in \mathbb{R}$, there exists a unique constant value of the forcing input $\pi_A = \bar{\pi}_A$ (and a unique position of the robots $q_R = \bar{q}_R$) such that $\bar{q} = (\bar{q}_L, \bar{q}_R)$ is an equilibrium of the system.

In particular $\bar{\pi}_A$ and $\bar{q}_R = [\bar{p}_{R1} \bar{p}_{R2}]^\top$ are given by

$$ \bar{\pi}_A(\bar{q}_L, t_L) = \mathbf{K}_A \bar{q}_R + \bar{f}(\bar{q}_L, t_L) \quad (6) $$

$$ \bar{p}_{Ri}(\bar{q}_L, t_L) = \bar{p}_L + \bar{R}_L^L b_i + \left( \frac{\|\bar{f}_i\|}{k_i} + l_{0i} \right) \frac{\bar{f}_i}{\|\bar{f}_i\|}, \quad (7) $$

for $i = 1, 2$, where

$$ \bar{f}(\bar{q}_L, t_L) = \begin{bmatrix} \bar{f}_1 \\ \bar{f}_2 \end{bmatrix} = \begin{bmatrix} \frac{b_2 m_L g}{L} \\ \frac{b_1 m_L g}{L} \end{bmatrix} \begin{bmatrix} I_3 \\ I_3 \end{bmatrix} e_3 + t_L \begin{bmatrix} I_3 \\ -I_3 \end{bmatrix} \bar{R}_L e_1. \quad (8) $$

From (6), we can see that the forcing input is made up of two parts: one that depends on the robots' positions computed from the load equilibrium configuration according to kinematic relations, and the other one that depends on the equilibrium forces. The equilibrium forces are composed, according to (8), by one term that compensates the gravity and one term that produces an internal force on the load whose intensity is $t_L$. [33] confirms that, if $\bar{\pi}_A$ is exactly applied to the closed-loop system (4), $\bar{q}_L$ is an isolated load equilibrium configuration if $t_L \neq 0$, which is asymptotically stable if $t_L > 0$ and unstable if $t_L < 0$. Instead, $\bar{q}_L$ belongs to a continuum of equilibrium points containing any possible attitude of the load if $t_L = 0$. In the remainder, Sec. III-Sec. V contain the main theoretical contributions of the work.

### III. EQUILIBRIA UNDER UNCERTAINTY

In this section, the uncertainties are introduced and the equilibrium configurations of the system subject to those uncertainties are derived. Note that, in reality, $\bar{\pi}_A$ in (6) cannot be applied exactly because of parametric uncertainties. Instead, one can apply only a version of $\bar{\pi}_A$, denoted with $\hat{\pi}_A$, computed using the nominal, uncertain values of the system parameters, (see Figure 2 for a schematic representation of the control scheme with the nominal forcing input). In the following, if not differently stated, we consider the general case in which a whole set of uncertainties are present. These uncertainties affect the control law (6) and, in turn, affect the system equilibrium configurations. The uncertainties are the following:

* $m_L$ is unknown, but only its nominal value $\hat{m}_L$ is available for the control design. We define the corresponding uncertainty as $\Delta_m = m_L - \hat{m}_L$;

* $b_1$ is unknown, but only its nominal value $\hat{b}_1$ is available. The corresponding uncertainty, affecting the load CoM position, is $\Delta_b = b_1 - \hat{b}_1$;

* L is unknown, but only its nominal value $\hat{L}$ is available, and we define $\Delta_L = \frac{1}{L} - \frac{1}{\hat{L}} = \ell - \hat{\ell}$ and $\Delta_L = L - \hat{L}$;

* the model of the cable i-th is inexact. Therefore, the nominal length $l_{0i}$ and stiffness $k_i$ are unknown, but their nominal values $\hat{l}_{0i}$ and $\hat{k}_i$ are available for the control design. We define the uncertainties $\Delta_{ki} = k_i - \hat{k}_i$, $\Delta_{l0i} = l_{0i} - \hat{l}_{0i}$.

Note that the nominal value of $b_2$, $\hat{b}_2$, depends on the previously defined quantities according to the relationship $\hat{b}_2 = \hat{L} - \hat{b}_1$. However, for convenience, we also define $\Delta_{b2} = b_2 - \hat{b}_2$.

We shall now study the system's equilibrium configurations when $\hat{\pi}_A$ is applied.

**Theorem 2** (equilibrium direct problem). Given a desired load configuration $\bar{q}_L = (\bar{p}_L, \bar{R}_L)$ and the internal force $t_L \in \mathbb{R}$, assume that the forcing input $\hat{\pi}_A$ is computed from (6) and is applied to the closed-loop system (4). Then, the equilibrium configurations are all and only the ones satisfying the following conditions:

$$ \mathbf{p}_{R1} = \hat{\mathbf{p}}_{R1} - \mathbf{K}_{A1}^{-1} (\Delta_m g e_3) := \mathbf{p}_{R1}^{eq} \quad (9) $$

$$ \mathbf{R}_L := \mathbf{R}_L^{eq} \quad \text{s.t.} \quad \mathbf{S}(e_1) \mathbf{R}_L^{eq} \top \left[ \left( b_1 m_L - \frac{\hat{b}_1 \hat{m}_L L}{\hat{L}} \right) g e_3 + L t_L \bar{R}_L e_1 \right] = \mathbf{0} \quad (10) $$

$$ \mathbf{f}_1 = m_L g e_3 - \frac{\hat{m}_L \hat{b}_1 g}{\hat{L}} e_3 + t_L \bar{R}_L e_1 := \mathbf{f}_1^{eq} \quad (11) $$

$$ \mathbf{f}_2 = \frac{\hat{b}_1 \hat{m}_L g}{\hat{L}} e_3 - t_L \bar{R}_L e_1 = \hat{\mathbf{f}}_2 := \mathbf{f}_2^{eq} \quad (12) $$

$$ \mathbf{p}_L = \mathbf{p}_{R1}^{eq} - \mathbf{R}_L^{eqL} b_1 - \left( \frac{\|\mathbf{f}_1^{eq}\|}{k_1} + l_{01} \right) \frac{\mathbf{f}_1^{eq}}{\|\mathbf{f}_1^{eq}\|} := \mathbf{p}_L^{eq}, \quad (13) $$

where $\hat{\mathbf{p}}_{R1}$ indicates the reference position of the leader robot computed as in (7), namely starting from $\bar{p}_L$, $\bar{R}_L$, but using the uncertain parameters.

**Proof.** $\hat{\pi}_A$ is defined according to (6), where (8) becomes

$$ \hat{\mathbf{f}}(\bar{q}_L, t_L) = \begin{bmatrix} \hat{\mathbf{f}}_1 \\ \hat{\mathbf{f}}_2 \end{bmatrix} = \begin{bmatrix} \frac{(\hat{L} - \hat{b}_1) \hat{m}_L g}{\hat{L}} \\ \frac{\hat{b}_1 \hat{m}_L g}{\hat{L}} \end{bmatrix} \begin{bmatrix} I_3 \\ I_3 \end{bmatrix} e_3 + t_L \begin{bmatrix} I_3 \\ -I_3 \end{bmatrix} \bar{R}_L e_1. \quad (14) $$

The control (3) is

$$ \mathbf{u}_{Ri} = \mathbf{M}_{Ai}^{-1} (-\mathbf{B}_{Ai} \hat{\mathbf{p}}_{Ri} - \mathbf{K}_{Ai} \mathbf{p}_{Ri}) - \mathbf{f}_i + \hat{\pi}_{Ai}. $$

Consider the equilibrium condition

$$ \mathbf{0} = m(\mathbf{q}, 0, \hat{\pi}_A). \quad (15) $$

Equation (12) is obtained by substituting the last three lines of (14) into (15) and solving the equilibrium condition for the follower robot. Then, (12) can be substituted into the load translational equilibrium (lines 7, 8, and 9 of (15)) to retrieve (11). (9) results from the first three lines of (15) using (11). Finally, (10) can be obtained using (11) and (12) in the last
<figure><img src="image_3.png" /><figcaption>Fig. 3: Representation of the different equilibrium configurations of the system depending on $t_L$ and $\xi$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4: Load attitude in the asymptotically stable equilibrium points $\mathcal{R}^+(t_L, \bar{q}_L)$ when $t_L \neq 0$. In dashed green, the desired configuration, and in grey the system with the actual attitude at the equilibrium. The leader robot is red. In $\mathcal{R}^+(t_L, \bar{q}_L)$, the attitude at the equilibrium is always such that $R_L^{eq}e_1$ is aligned with $t_L \bar{R}_L e_1$. However, for $t_L < 0$, this means the system is 'flipped' compared to the desired configuration.</figcaption></figure>

three lines of (15). Equation (13) is obtained applying the analogous of (42). □

**Definition 3.** Given a desired load configuration $\bar{q}_L = (\bar{p}_L, \bar{R}_L)$, internal force $t_L \in \mathbb{R}$, and forcing input $\pi_A = \hat{\pi}_A(\bar{q}_L, t_L)$, we define the set of equilibrium configurations as $\mathcal{Q}(t_L, \bar{q}_L) = \{q \text{ s.t. conditions of Theorem 2 are satisfied}\}$

From Theorem 2, we can distinguish between two scenarios:

**Scenario 1:** If $t_L = 0$, condition (10) implies that the attitude of the load is such that $R_L^{eq}e_1$ is aligned to $e_3$, and conditions (11) and (12) imply that both cables are vertical. In other words, the load at the equilibrium is, irrespective of the parametric uncertainties, aligned with the vertical direction; even an infinitesimal parametric uncertainty would lead the load to this undesired configuration in which the vertical load is aligned with the two vertical cables. Such a configuration is clearly not realizable. Note that the position error of the system at the equilibrium still depends on the parametric uncertainties (see condition (9)). One can express the alignment between $R_L^{eq}e_1$ and $e_3$ as $(R_L^{eq}e_1)^\top e_3 = \pm 1$. By convention, let us indicate with $\mathcal{Q}_1(0, \bar{q}_L)$ the system equilibrium configuration in which $(R_L^{eq}e_1)^\top e_3 = +1$ holds, namely the one in which the leader robot is above and the follower robot below, and with $\mathcal{Q}_2(0, \bar{q}_L)$ the other equilibrium configuration. Fig 3a illustrates the aforementioned equilibrium configurations. Note also that there is an additional possibility. With simple ma-

nipulation, remembering that $\hat{b}_1 = b_1 + \Delta_b$, $\hat{b}_2 = b_2 + \Delta_{b2}$, $\hat{L} = L + \Delta_L$, and defining $\xi$ as follows, the term in (10) becomes:

$$ \xi := \left( b_1 m_L - \frac{\hat{b}_1 \hat{m}_L L}{\hat{L}} \right) = \Delta_m b_1 + \frac{\Delta_{b2}}{\hat{L}} \hat{m}_L b_1 - \frac{\Delta_{b1}}{\hat{L}} \hat{m}_L b_2. \quad (16) $$

If $\xi = 0$, (10) is verified for every value of $R_L$, and hence the equilibrium configurations $\mathcal{Q}(0, \bar{q}_L)$ are infinite and such that the attitude of the load at equilibrium is arbitrary. This happens in the special case in which the parameters of the system are exactly known (this situation is the one we analyzed in [33] and which we can now see as a special case with $\xi = 0$). Indeed, $\xi = 0$ is verified also if the cable parameters are the sole uncertain ones, as it will be also more deeply discussed in the following. See 3b for a schematic representation of the mentioned equilibrium configurations.

**Scenario 2:** If $t_L \neq 0$, condition (10) holds when the vectors $R_L^{eq}e_1$ and

$$ (\xi g e_3 + L t_L \bar{R}_L e_1) \quad (17) $$

are aligned. Similar to before, this condition holds in two possible cases: when the vectors are aligned and point in the same direction, or when they are aligned but point in opposite directions. Let us indicate with $R_L^{eq+}$ the attitude of the load for which condition (10) holds and the two vectors $R_L^{eq}e_1$ and
(17) point in the same direction. We indicate the corresponding load equilibrium configuration as $\mathcal{Q}^+(t_L, \bar{q}_L)$. In the other case, when the two aforementioned vectors point in opposite directions, at the equilibrium one has that $R_L^{eq-} = R_L^{eq+}R_{z_L}(\pi)$; we indicate the corresponding equilibrium configuration as $\mathcal{Q}^-(t_L, \bar{q}_L)$. Depending on the sign of $t_L$ in $\hat{\pi}_A$, the forces in the cables place the load under tension in one equilibrium configuration and under compression in the other. Figure 3c represents these equilibrium configurations.

Remark 2. Under the hypothesis that $\bar{\theta} \neq \pi/2 + k\pi$, with $k \in \mathbb{N}$ and $t_L \neq 0$, as shown by (10), at the equilibrium

$$ \psi = \bar{\psi} + k\pi \quad (18) $$

$$ \tan \theta = \tan \bar{\theta} + \frac{-\xi g}{L t_L \cos \bar{\theta}}. \quad (19) $$

In other words, the uncertainties have no effect on the yaw angle at equilibrium. $\psi$ may differ from $\bar{\psi}$ by $\pi$ because, as already discussed, both $\mathcal{Q}^+(t_L, \bar{q}_L)$ and $\mathcal{Q}^-(t_L, \bar{q}_L)$ are equilibrium configurations. Moreover, (19) tells us that not only is the attitude error proportional to the amount of uncertainty but also that, as $t_L$ decreases, the load at the equilibrium becomes increasingly vertical. Eventually, for $t_L = 0$ and uncertain parameters ($\xi \neq 0$), (10) leads to $e_1 \times R_L^\top e_3 = 0$. Namely, as previously observed, the load at the equilibrium is aligned with the vertical direction and the two cables are vertical despite the value of $\xi \neq 0$. In other words, if $t_L = 0$ the load attitude error is unaffected by the parametric uncertainties: the load will reach the same, clearly undesired, configuration regardless of the smallest $\xi \neq 0$.

In the remainder of this section, we briefly analyze the effects of each uncertain parameter on the final equilibrium.

### A. Uncertainty on the load mass $m_L$

In this subsection, we only discuss uncertainty in the load's mass, while the other parameters are assumed to be perfectly known. Equations (9)-(12) become:

$$ \mathbf{p}_{R1}^{eq} = \hat{\mathbf{p}}_{R1} - \mathbf{K}_{A1}^{-1} \Delta_m g e_3 \quad (20) $$

$$ b_1 \mathbf{S}(e_1) \mathbf{R}_L^{eq\top} g \Delta_m e_3 + t_L L \mathbf{S}(e_1) \mathbf{R}_L^{eq\top} \bar{\mathbf{R}}_L e_1 = 0 \quad (21) $$

$$ \mathbf{f}_1^{eq} = m_L g e_3 - \frac{b_1 \hat{m}_L g}{L} e_3 + t_L \bar{\mathbf{R}}_L e_1 = \hat{\mathbf{f}}_1 + \Delta_m g e_3 \quad (22) $$

$$ \mathbf{f}_2^{eq} = \frac{b_1 \hat{m}_L g}{L} e_3 - t_L \bar{\mathbf{R}}_L e_1 = \hat{\mathbf{f}}_2. \quad (23) $$

The position of the load CoM at the equilibrium is different from $\bar{p}_L$ and can be computed from (13) using (20)-(22).

It is worth noting that the leader robot can detect a mismatch between the known commanded $\hat{\mathbf{f}}_1$ and the actual force $\mathbf{f}_1^{eq}$ measured at steady state. Such a discrepancy solely depends on $\Delta_m$, according to (22). Thus, the leader robot can compute $\Delta_m$ and, by knowing the nominal value $\hat{m}_L$, retrieve the actual value of the load mass $m_L$, which can be used to adjust its own reference force and position. However, note that in a communication-less setup it is impossible for both robots to know the correct parameter value based simply on their own state. In fact, according to (23), the follower robot has no mismatch between the equilibrium and the force reference value.

### B. Uncertainty on the load length, L, or CoM position, $b_1$

Uncertainties on one of these two parameters have similar effects. In one case, $\hat{b}_1 \neq b_1$, namely the load CoM is aligned to the cables attachment points on the load at an uncertain position but L is exactly known; in the other case, $b_1$ is exactly known but L is not. In both cases, at the equilibrium, the following conditions hold:

$$ \mathbf{p}_R^{eq} = \hat{\mathbf{p}}_{R1} \quad (24) $$

$$ \mathbf{S}(e_1) \mathbf{R}_L^{eq\top} (\mathrm{Lt}_L \bar{\mathbf{R}}_L e_1 + y m_L g e_3) = 0 \quad (25) $$

$$ \mathbf{f}_1^{eq} = \hat{\mathbf{f}}_1 \quad (26) $$

$$ \mathbf{f}_2^{eq} = \hat{\mathbf{f}}_2, \quad (27) $$

where $y = \Delta_b$ in one case, and $y = b_1 L \Delta_\ell$ in the other. $\hat{\mathbf{p}}_{R1}$, $\hat{\mathbf{f}}_1$, and $\hat{\mathbf{f}}_2$ are computed from (7) and (8), where the corresponding uncertain parameter is used in place of the real one.

In this case, the leader robot position and both cable forces at the equilibrium coincide with the respective reference values available to the robots (see (24)-(27)). Consequently, it is not possible for any of the robots to estimate the uncertain parameter at the equilibrium based on the local information they possess.

### C. Uncertainty on the cable length $l_{0i}$ or stiffness $k_i$

Consider an uncertainty on the parameters of the i-th cable such that the rest length is $l_{0i} \neq \hat{l}_{0i}$ and the stiffness is $k_i \neq \hat{k}_i$. At the equilibrium, $\mathbf{f}_i^{eq} = \hat{\mathbf{f}}_i$,

$$ \mathbf{R}_L^{eq} = \bar{\mathbf{R}}_L. \quad (28) $$

$$ \mathbf{p}_{R1}^{eq} = \hat{\mathbf{p}}_{R1} = \bar{\mathbf{p}}_L + \bar{\mathbf{R}}_L^L \mathbf{b}_1 + \left( \frac{\|\bar{\mathbf{f}}_1\|}{\hat{k}_1} + \hat{l}_{01} \right) \frac{\bar{\mathbf{f}}_1}{\|\bar{\mathbf{f}}_1\|}, \quad (29) $$

and the value of $\mathbf{p}_L$ at the equilibrium is

$$ \mathbf{p}_L^{eq} = \hat{\mathbf{p}}_{R1} - \bar{\mathbf{R}}_L^L \mathbf{b}_1 - \left( \frac{\|\bar{\mathbf{f}}_1\|}{k_1} + l_{01} \right) \frac{\bar{\mathbf{f}}_1}{\|\bar{\mathbf{f}}_1\|} \quad (30) $$

We highlight that knowledge about the cable properties is required only when computing the reference position of the leader robot, according to (7). What is more, only the information about $k_1$ and $l_{01}$ is required. We conclude that knowledge of $l_{02}$ and $k_2$, is not necessary to stabilize the load at a desired pose. Moreover, $l_{01}$ and $k_1$ have no effect on the load attitude at equilibrium but they do influence the load position. This is evident by simply substituting (29) into (30) with $l_{01} \neq \hat{l}_{01}$ and $k_1 \neq \hat{k}_1$. Note that, since the robots' forces and the leader robot's position at the equilibrium coincide with the reference values available to the robots themselves, they are unaware of the load pose error induced by this uncertainty.

## IV. STABILITY ANALYSIS

In this section, we shall analyze the stability of the equilibrium configurations discovered in Sec III. First, being $\mathbf{x} = (\mathbf{q}, \mathbf{v})$ the state of the system, we define the following equilibrium states (subspaces of the state space):

$$ \mathcal{X}(0, \bar{q}_L) = \{ \mathbf{x} : \mathbf{q} \in \mathcal{Q}(0, \bar{q}_L), \mathbf{v} = 0 \}, $$
• $\mathcal{X}_1(0, \bar{q}_L) = \{x : q \in \mathcal{Q}_1(0, \bar{q}_L), v = 0\}$,

• $\mathcal{X}_2(0, \bar{q}_L) = \{x : q \in \mathcal{Q}_2(0, \bar{q}_L), v = 0\}$,

• $\mathcal{X}^+(t_L, \bar{q}_L) = \{x : q \in \mathcal{Q}^+(t_L, \bar{q}_L), v = 0\}$,

• $\mathcal{X}^-(t_L, \bar{q}_L) = \{x : q \in \mathcal{Q}^-(t_L, \bar{q}_L), v = 0\}$.

**Theorem 3.** Let us consider a desired load configuration $\bar{q}_L$. For the system (4), let the constant forcing input be $\hat{\pi}_A$. Then,

• $\mathcal{X}_1(0, \bar{q}_L)$ is asymptotically stable if $\xi > 0$ and unstable if $\xi < 0$;

• $\mathcal{X}_2(0, \bar{q}_L)$ is asymptotically stable if $\xi < 0$ and unstable if $\xi > 0$.

• $\mathcal{X}(0, \bar{q}_L)$ is a set of marginally stable equilibrium points if $\xi = 0$.

• $\mathcal{X}^+(t_L, \bar{q}_L)$ is asymptotically stable

• $\mathcal{X}^-(t_L, \bar{q}_L)$ is unstable.

*Proof.* Consider the following Lyapunov candidate function:

$$V(x) = \frac{1}{2}(v_R^\top M_A v_R + e_R^\top K_A e_R + v_L^\top M_L v_L + \\ + k_1 (\|l_1\| - l_{01})^2 + k_2 (\|l_2\| - l_{02})^2 - l_1^\top f_1^{eq} + \\ - l_2^\top f_2^{eq} + V_0 + V_R(x), \quad (31)$$

where the robot position error is $e_R = p_R - p_R^{eq}$, $V_0$ is constant, and $V_R(x)$ is an additional term explained in the following. Function (31) is composed of standard positive definite quadratic terms equal to zero in the equilibrium points and by two terms of the form $\frac{1}{2}k_i(\|l_i\| - l_{0i})^2 - l_i^\top f_i^{eq}$, call them $V_i(x)$: these are linked to the elastic energy of the cables and have a minimum at the equilibrium as well. A detailed proof of the former point can be found in [33]. The proof first shows that $V_i(x)$ is radially unbounded, i.e., $\lim_{\|x\| \to \infty} V_i(x) = \infty$. Then, based on this result and Theorem 1.15 of [47], the term has a global minimum. Finally, it has been shown that the global minimum of $V_i(x)$ corresponds to the considered equilibrium [33].

We define the value of $V_i(x)$ at the equilibrium (its minimum value) as $-V_0$, and we cancel it in (31) so that its value at the equilibrium is zero.

Let us start considering $\mathcal{X}_1(0, \bar{q}_L)$ and $\xi > 0$. In this case, we set $V_R(x) = \xi g(1 - e_3^\top R_L e_1)$. With this choice, (31) is zero in $\mathcal{X}_1(0, \bar{q}_L)$ because also the term $1 - e_3^\top R_L e_1$ is zero in $\mathcal{X}_1(0, \bar{q}_L)$ by definition (load aligned with the vertical with $R_L e_1$ and $e_3$ pointing in the same direction) and positive elsewhere (the scalar product $e_3^\top R_L e_1 \le 1$ because $e_3$ and $R_L e_1$ have both unit norm).

Studying the sign of the time derivative of (31), using (4), (2), and (8), we obtain $\dot{V}(x) = -v_R^\top B_A v_R$, which is clearly negative semidefinite. In particular, let us define $\mathcal{E} = \{x : \dot{V}(x) = 0\}$. In this case, we have $\mathcal{E} = \{x : v_R = 0, \omega_L = 0\}$.

Since $\dot{V}(x)$ is only negative semidefinite, we rely on LaSalle's invariance principle to complete the proof: one can easily verify from (4) that the largest invariant set in $\mathcal{E}$ is $\mathcal{X}_1(0, \bar{q}_L)$.

Analogous reasoning can be used when $\xi < 0$. The computation of $\dot{V}$ does not change, and it is, thus, negative semidefinite. However, $\mathcal{X}_1(0, \bar{q}_L)$ is a set of accumulation for the points where $V(x) < 0$ if $\xi < 0$. To see this, consider

$v = 0$ and all quantities at the equilibrium apart from $R_L$, which is such that $e_3^\top R_L e_1 = 1 - \epsilon$, with $\epsilon > 0$ arbitrarily small, meaning that $R_L$ is arbitrarily close to $R_L^{eq}$. Under these conditions, $V(x) = g\xi\epsilon < 0$. All conditions of Chetaev's theorem (the formulation of both this and La Salle's invariance principle can be found, e.g., in [48]) are satisfied. Hence, we can conclude that $\mathcal{X}^+(t_L, \bar{q}_L)$ is unstable. To show that $\mathcal{X}_2(0, \bar{q}_L)$ is asymptotically stable if $\xi < 0$, we set $V_R(x) = -\xi g(1 + e_3^\top R_L e_1)$, which is zero in $\mathcal{X}_2(0, \bar{q}_L)$ (when $e_3^\top R_L e_1 = -1$ by definition), and positive elsewhere. The same Lyapunov candidate function is used to show that $\mathcal{X}_2(0, \bar{q}_L)$ is unstable if $\xi > 0$. The reasoning is exactly dual to the previous case, hence it is here omitted for the sake of space.

Consider now $\xi = 0$. We set $V_R(x) = 0$, so that (31) is zero in $\mathcal{X}(0, \bar{q}_L)$ and positive elsewhere. Moreover, $\dot{V}$ is negative semidefinite as before. One can easily show that the largest invariant set is $\mathcal{X}(0, \bar{q}_L)$. We can thus say that the system state converges to a state $x \in \mathcal{X}(0, \bar{q}_L)$, which is, however, composed of a continuum of equilibrium points. Hence they are only marginally stable.

Finally, we study the stability of the equilibrium points when $t_L \ne 0$. In this case, we set $V_R(x) = -(\xi g e_3 + t_L L \bar{R}_L e_1)^\top R_L e_1 + V_0'$, with $V_0' = (\xi g e_3 + t_L L \bar{R}_L e_1)^\top R_L^{eq} e_1$. Clearly, $V_R(x)$ and hence $V(x)$ are zero at the equilibrium. Moreover, $V_R(x)$ is positive elsewhere by definition of $\mathcal{X}^+(t_L, \bar{q}_L)$, $(\xi g e_3 + t_L L \bar{R}_L e_1)$ and $R_L e_1$ are aligned and point in the same direction when $R_L = R_L^{eq}$, so that $V_R(x)$ has its minimum in $\mathcal{X}^+(t_L, \bar{q}_L)$. Moreover, $\dot{V}(x) = -v_R^\top B_A v_R$, and the application of LaSalle's invariance principle leads to the conclusion that $\mathcal{X}^+(t_L, \bar{q}_L)$ is an asymptotically stable equilibrium point, similarly to before. To show the instability of $\mathcal{X}^-(t_L, \bar{q}_L)$, we use the same choice for $V_R(x)$. However, since $\xi g e_3 + t_L L \bar{R}_L e_1$ and $R_L e_1$ are anti-parallel in $\mathcal{X}^-(t_L, \bar{q}_L)$, $V_R(x)$ is still zero at the equilibrium but negative when $R_L$ is arbitrarily close to $R_L^{eq}$. $\mathcal{X}^-(t_L, \bar{q}_L)$ is a point of accumulation for the points in which $\dot{V}(x)$ is negative, while $\dot{V}(x)$ remains negative semi-definite. For Chetaev's theorem, we conclude that $\mathcal{X}^-(t_L, \bar{q}_L)$ is unstable. $\square$

It is important to highlight that, as shown in Figure 4, for $t_L > 0$, $\mathcal{X}^+(t_L, \bar{q}_L)$ corresponds to a configuration of the system in which $R_L^{eq}$ (irrespective of the sign of $\xi$) is the closest condition to $\bar{R}_L$, namely to the desired attitude, with a displacement due to the parametric uncertainty. Instead, for $t_L < 0$, is $\mathcal{X}^-(t_L, \bar{q}_L)$ the equilibrium point in which the configuration of the load is the closest to the desired one. We can say that these configurations are the most desirable equilibrium configuration of the load in the presence of parametric uncertainties. As stated in Theorem 3, $t_L > 0$ stabilizes the most desirable equilibrium configuration of the load, which is, instead, unstable if $t_L < 0$.

V. THE ROLE OF THE INTERNAL FORCES ON THE LOAD ERROR CAUSED BY PARAMETRIC UNCERTAINTIES

In this section, we provide a formal analysis of the role that the internal force plays in determining the load pose at the
equilibrium in the presence of parametric uncertainties. We shall consider the simultaneous presence of all the uncertainties listed in Sec III. We start considering the load attitude.

### A. Load attitude error

**Theorem 4.** The load attitude error at the equilibrium $e_{R_L}$, is inversely proportional to the intensity of a positive internal force $t_L$. Furthermore, defining

$$e_{R_L} = \| \mathbf{R}_L e_1 \times \bar{\mathbf{R}}_L e_1 \|^2, \quad (32)$$

the error sensitivity w.r.t. $\Delta_m, \Delta_b, \Delta_{k_i}, \Delta_{l_{0i}}, \Delta_l$, defined as $\frac{\partial e_{R_L}}{\partial \Delta_m}, \frac{\partial e_{R_L}}{\partial \Delta_{k_i}}, \frac{\partial e_{R_L}}{\partial \Delta_{l_{0i}}}, \frac{\partial e_{R_L}}{\partial \Delta_l}$, respectively, is given by:

$$\frac{\partial e_{R_L}}{\partial \Delta_m} = \frac{-2\hat{b}_1 \hat{l} g^2 \alpha \cos \theta^2}{t_L^2 L^2} \quad (33)$$

$$\frac{\partial e_{R_L}}{\partial \Delta_b} = \frac{-2\hat{m}_L \hat{l} g^2 \alpha \cos \theta^2}{t_L^2 L^2} \quad (34)$$

$$\frac{\partial e_{R_L}}{\partial \Delta_{k_i}} = \frac{\partial e_{R_L}}{\partial \Delta_{l_{0i}}} = 0 \quad (35)$$

$$\frac{\partial e_{R_L}}{\partial \Delta_l} = \frac{-2\hat{b}_1 \hat{m}_L g^2 \alpha \cos \theta^2}{t_L^2 L^2} \quad (36)$$

where

$$\alpha := (b_1 - \Delta b)(m_L - \Delta m)(l - \Delta_l) - m_L b_1.$$

*Proof.* With a positive internal force, at the equilibrium, (18) with $k=0$ holds. Thus, the quantity $|\theta^{eq} - \bar{\theta}|$ is a viable indicator of the attitude error. According to (19) and for monotonicity of the tan() function, the difference between $\bar{\theta}$ and $\theta^{eq}$ varies with the quantity $\xi g/(L t_L \cos \bar{\theta})$, which is inversely proportional to $t_L$. Hence, also the error is. Rewrite now (10) in $\mathcal{F}_W$ as:

$$\mathbf{R}_L^{eq} e_1 \times \left[ \left( b_1 m_L - \frac{\hat{b}_1 \hat{m}_L L}{\hat{L}} \right) g e_3 + L t_L \bar{\mathbf{R}}_L e_1 \right] = \mathbf{0}. \quad (37)$$

Define also:

$$\frac{\hat{b}_1 \hat{m}_L L - b_1 m_L \hat{L}}{t_L L \hat{L}} (\mathbf{R}_L^{eq} e_1 \times g e_3) := \mathbf{x}. \quad (38)$$

Thus, from (37), we have that $\mathbf{R}_L^{eq} e_1 \times \bar{\mathbf{R}}_L e_1 = \mathbf{x}$ and, from (32), that $e_{R_L} = \mathbf{x}^\top \mathbf{x}$. Regarding the sensitivity, we show the proof for (33) only, because the other cases follow the exactly same analysis. We can write the sensitivity as:

$$\frac{\partial e_{R_L}}{\partial \Delta_m} = 2\mathbf{x}^\top \frac{\partial \mathbf{x}}{\partial \Delta_m} = 2\left[ \frac{1}{t_L L} \mathbf{R}_L e_1 \times (\alpha) g e_3 \right]^\top \left[ \frac{1}{t_L L} \mathbf{R}_L e_1 \times (\Delta b - b_1) g e_3 \right] \quad (39)$$

Eventually, (39) can be rewritten as (33) by remembering that, given three vectors $\mathbf{a}, \mathbf{b}$, and $\mathbf{c}$

$$(\mathbf{a} \times \mathbf{b})^\top (\mathbf{a} \times \mathbf{c}) = |\mathbf{a}|^2 (\mathbf{b}^\top \mathbf{c}) - (\mathbf{a}^\top \mathbf{b})(\mathbf{a}^\top \mathbf{c}).$$

Note that we are considering $t_L \neq 0$ by assumption. □

**Remark 3.** The definition in (32) is a suitable metric for the attitude error. If we consider the equilibrium point $\mathcal{X}^+(t_L, \bar{q}_L)$,

namely the one in which the displacement between $\mathbf{R}_L^{eq} e_1$ and $\bar{\mathbf{R}}_L e_1$ is the smallest, and hence our desired equilibrium point. Firstly, $\mathbf{R}_L e_1$ is enough to describe the entire attitude of the beam-like load. Secondly, $e_{R_L}$ is zero when $\mathbf{R}_L^{eq} = \bar{\mathbf{R}}_L$ and increases with the displacement between the two vectors $\mathbf{R}_L^{eq} e_1$ and $\bar{\mathbf{R}}_L e_1$, at least locally (for displacements smaller than $\pm \pi/2$).

Moreover, Theorem 4 shows that increasing the intensity of the internal force $t_L$ not only makes the attitude error smaller in presence of parametric uncertainties, but it also makes the error more robust to variations of such uncertainties.

This last aspect may be of particular practical interest: as a matter of fact, parametric uncertainty variations take place every time the actual physical parameters of the system change. A possible real-world scenario is the transportation of objects that are slightly different from each other, e.g., in mass and length. One may want to transport the objects without changing every time the controller parameters for the sake of time, thus dealing with varying parametric uncertainties. Especially interesting, as also highlighted in [49], is the variation affecting the CoM position, which may change online when transporting moving masses, i.e. containers of liquids, or boxes with smaller objects free to move inside. The previous analysis suggests that in all these cases having a larger value of $t_L$ is of uttermost benefit, resulting in an error less sensitive to the aforementioned parametric variations.

### B. Load position error

Differently from what happens to the load attitude error, the load position error at the equilibrium does not necessarily decrease when $t_L$ increases. While to claim a positive statement a comprehensive proof is needed, as we did in Sec V-A, to deny a positive statement, as we do in this section, a counterexample is enough. First, it is easy to see that, when only $l_{01}$ is uncertain, the load position error at the equilibrium, $e_{p_L} := \mathbf{p}_L^{eq} - \bar{\mathbf{p}}_L$ is

$$e_{p_L} = \Delta_l \frac{b_2 m g e_3 + t_L \bar{\mathbf{R}}_L e_1}{\| b_2 m g e_3 + t_L \bar{\mathbf{R}}_L e_1 \|}. \quad (40)$$

Eq. (40) suggests that $e_{p_L}$ is equal to a unit vector multiplied by $\Delta_{l_{01}}$, thus, its module is independent of the value of $t_L$. Moreover, in the next section, we provide two numerical examples showing that, depending on the specific combination and values of uncertainties, $e_{p_L}$ may even have a non-monotonic evolution for increasing values of $t_L$, with an initial increase or decay.

We show, however, that the load position error at the equilibrium can be corrected, ideally to zero, without altering the leader-follower architecture nor requiring direct communication between the robots.

We recall that, due to parametric uncertainties, the reference position given to the leader robot is

$$\hat{\mathbf{p}}_{R1} = \bar{\mathbf{p}}_L + \bar{\mathbf{R}}_L^L \hat{\mathbf{b}}_1 + \left( \frac{\| \hat{\mathbf{f}}_1 \|}{\hat{k}_1} + \hat{l}_{01} \right) \frac{\hat{\mathbf{f}}_1}{\| \hat{\mathbf{f}}_1 \|}. \quad (41)$$
<figure><img src="image_5.png" /><figcaption>Fig. 5: Sagging cable model in the simulated scenario.</figcaption></figure>

By using kinematics, (9), and (41), the load position at the equilibrium is

$$ 
\begin{aligned}
\boldsymbol{p}_{L}^{eq} &= \boldsymbol{p}_{R1}^{eq} - \left( \frac{\|\boldsymbol{f}_{1}^{eq}\|}{k_{1}} + l_{01} \right) \frac{\boldsymbol{f}_{1}^{eq}}{\|\boldsymbol{f}_{1}^{eq}\|} - \boldsymbol{R}_{L}^{eq} \boldsymbol{b}_{1} \\
&= \bar{\boldsymbol{p}}_{L} + \bar{\boldsymbol{R}}_{L}^{L} \hat{\boldsymbol{b}}_{1} + \left( \frac{\|\hat{\boldsymbol{f}}_{1}\|}{\hat{k}_{1}} + \hat{l}_{01} \right) \frac{\hat{\boldsymbol{f}}_{1}}{\|\hat{\boldsymbol{f}}_{1}\|} - K_{A}^{-1} \Delta_{m} g e_{3} + \\
&\quad - \left( \frac{\|\boldsymbol{f}_{1}^{eq}\|}{k_{1}} + l_{01} \right) \frac{\boldsymbol{f}_{1}^{eq}}{\|\boldsymbol{f}_{1}^{eq}\|} - \boldsymbol{R}_{L}^{eq} \boldsymbol{b}_{1}.
\end{aligned}
\quad (42) $$

From (42), we have an expression for $e_{pL}$. Now, if the leader robot knows the load position, it can recognize that, at steady state, $\bar{\boldsymbol{p}}_L \neq 0$ holds, and it can adjust its position reference to $^2\hat{\boldsymbol{p}}_{R1}$ accordingly, with

$$ ^2\hat{\boldsymbol{p}}_{R1} = \hat{\boldsymbol{p}}_{R1} - e_{pL}. \quad (43) $$

In this way, there will be a new equilibrium in which the leader robot position is

$$ \boldsymbol{p}_{R1}^{eq} = ^2\hat{\boldsymbol{p}}_{R1} - K_{A}^{-1} \Delta_{m} g e_{3} \quad (44) $$

and thus (42) becomes

$$ \boldsymbol{p}_{L}^{eq} = \bar{\boldsymbol{p}}_{L} \quad (45) $$

It is important to highlight that the leader robot position only influences the load position and not the attitude at the equilibrium, which depends only upon the reference forces computed based on $\bar{\boldsymbol{R}}_L$. Indeed, by evaluating (4) at the equilibrium, the last three rows are

$$ S(^{L}\boldsymbol{b}_{1}) \boldsymbol{R}_{L}^{eq \top} \boldsymbol{f}_{1}^{eq} + S(^{L}\boldsymbol{b}_{2}) \boldsymbol{R}_{L}^{eq \top} \boldsymbol{f}_{2}^{eq} = \boldsymbol{0}, $$

which becomes, substituting $\boldsymbol{f}_i^{eq}$, eq. (10). Hence, the leader robot can correct the load position error, while the internal force independently acts decreasing the attitude error. Because the load's position can be steered relying solely on the leader robot, unlike the load's attitude which is determined by the cooperative actions of both robots, the control approach can maintain its distributed nature. However, to correct the load position error, the leader robot must have access to the load position. This implies that the leader robot would require additional sensors, such as cameras to accomplish this task.

## VI. NUMERICAL VALIDATION

Extensive numerical simulations have been carried out using a URDF description of the system and ODE physics engine in Gazebo. We avoided validating the theoretical results on the same equations used to derive them. The main differences

between the control model used to derive a fully satisfactory theoretical analysis and the complex simulation model used to study the applicability of the theoretical results in the real world are in the following.

* Under-actuated quadrotors have been deliberately preferred for the validation since they represent the worst-case in terms of the validity of some of the assumptions made in the theoretical analyses. Validation using fully-actuated aerial robots would have seemed, instead, limiting.

* The cables are subject to sagging, which is obtained by using a series of several links interconnected by passive universal joints, as can be seen in Figure 5.

* In the validation, there is no guarantee of perfect trajectory tracking as assumed in the theory but a standard position controller [50] is implemented for each robot.

* The wrench observer proposed in [51] is used to estimate the force applied by the cable on the robot. The observer introduces noisy and delayed measurements when compared to ideal force measurement.

The control software has been implemented in Matlab-Simulink using the Generator of Modules GenoM<fnref n="4" />. The interface between Matlab and Gazebo is also managed by a Gazebo-genom3 plugin<fnref n="5" />. All phases of a physical experiment, starting with takeoff, are replicated in the simulated environment using a state machine, ensuring that the results are as realistic as possible. After the takeoff, the two robots lift the load, and the admittance controller is activated right after.

The robot models are two quadrotors weighing 1.03 kg and having a maximum thrust for each propeller of 6 N. They are equipped with two light cables of length 1 m and attached to a bar. The bar is a one-meter-long link with a mass of 0.5 kg.

1) Case of $t_L > 0$: Figure 6 contains the average load attitude error at a steady state in a total of 40 simulations. The average is computed in a 2-second time window. In all those simulations, $\bar{\boldsymbol{p}}_L = [1 \ 1 \ 1]^\top$ m, the load desired yaw is $\bar{\boldsymbol{\psi}} = \frac{\pi}{8}$ radians, and the desired pitch $\bar{\boldsymbol{\theta}} = -\frac{\pi}{12}$ radians. In each of the three plots in Figure 6, for four different values of the internal force, $t_L = \{0.5, 0.75, 1, 1.25\}$N, four different simulation results are displayed for each relative error equal to 0%, 5%, 10%, 15% on a specific uncertain parameter considered separately from the others. Specifically, Figure 6a, considers the uncertainty on $m_L$, Figure 6b on $L$, and Figure 6c on $b_1$. Even in the absence of uncertainties, small errors of less than 2.5 degrees in the bar's attitude control can be found. This can be due to minor tracking errors or possible biases in the wrench observer, which estimations are unbiased as soon as the robot takes off. These considerations ignore the external forces applied by the loose cables at the startup phase. From all the three figures one can appreciate the beneficial effect of larger values of $t_L$ on the attitude error: for the same value of the uncertainty, the attitude error decreases if the $t_L$ increases. Moreover, the plots show that for every value of $t_L$, increasing the uncertainty on one parameter increases the attitude error, as expected, but, especially, the increase

<footnote n="4">⁴https://git.openrobots.org/projects/genom3</footnote>

<footnote n="5">⁵https://git.openrobots.org/projects/mrsim-gazebo</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 6: Each point in the plots is a value of the average attitude error at steady state in a simulation with $t_L$ as indicated in the legend, and parametric uncertainty as indicated in the x-axis of the corresponding plot. A total of 40 simulation results are known in these plots.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7: Simulations for cable parameter uncertainties and $t_L = 1$N. Dotted lines of the same color indicate the corresponding desired quantities.</figcaption></figure>

is smaller for high values of $t_L$ (this can be seen by the slope of the lines in the plots). These results confirm the theoretical findings collected in Theorem 4. Figure 7a and 7b provide validation of the theoretical results on the effect of the uncertainties affecting the cable parameters. The leader robot position reference is not given as a step, but the robot follows a 5-th order polynomial trajectory to reach the desired position. An error of 15% is considered to affect the length of the leader and follower robot's cable in Figure 7a and 7b, respectively. In both cases, $t_L = 1$N. Note that the displayed time starts after the admittance controller activation. The reader can appreciate how uncertainties on the leader robot's cable model only cause $\dot{\boldsymbol{p}}_L^{eq} \neq \bar{\boldsymbol{p}}_L$, while $\boldsymbol{R}_L^{eq} = \bar{\boldsymbol{R}}_L$, and how the follower robot's cable parameters are not needed to control the load pose, as expected from (30) and (28).

2) Case of $t_L < 0$: Figure 7a and 7b provide validation of the theoretical results on the effect of the uncertainties affecting the cable parameters. The leader robot position reference is not given as a step, but the robot follows a 5-th order polynomial trajectory to reach the desired position. An error of 15% is considered to affect the length of the leader and follower robot's cable in Figure 7a and 7b, respectively.

In both cases, $t_L = 1$N. Note that the displayed time starts after the admittance controller activation. The reader can appreciate how uncertainties on the leader robot's cable model only cause $\dot{\boldsymbol{p}}_L^{eq} \neq \bar{\boldsymbol{p}}_L$, while $\boldsymbol{R}_L^{eq} = \bar{\boldsymbol{R}}_L$, and how the follower robot's cable parameters are not needed to control the load pose, as expected from (30) and (28).

3) Case of $t_L < 0$: Here we show the behavior of the system with $t_L < 0$. The unstable nature of the desired configuration with no parametric uncertainties was shown in [33]. The simulations of the realistic system, in accordance with Theorem 3, show that $\mathcal{X}^-(t_L, \bar{q}_L)$ is unstable when the parametric uncertainties are considered. We report the results of two simulations with $\bar{p}_{R1} = [0 \ 0 \ 1]^\top$ m, $\bar{\psi} = \frac{\pi}{8}$ radians, and $\bar{\theta} = 0$ radians (desired horizontal bar). We simulate an uncertainty of 5% both on $m_L$ and $L$ such that (i) $\xi > 0$ (we chose $\hat{m}_L < m_L$ and $\hat{L} > L$) and (ii) $\xi < 0$ (thanks to $\hat{m}_L > m_L$ and $\hat{L} < L$). In both cases, we obtained that, in accordance to Theorem 3, the system converges to $\mathcal{X}^+(t_L, \bar{q}_L)$. Since $t_L < 0$, that means, as reported in Figure 9, that we have $\psi^{eq} = \bar{\psi} - \pi$, while $\theta^{eq}$ varies according to the sign of $\xi$, as expected (see Figure 4). The cable forces were observed to be as desired, except for the vertical component of $\hat{\boldsymbol{f}}_1$, as expected due to
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Simulations for cable parameter uncertainties and $t_L = -1$N. Dotted lines of the same color indicate the corresponding desired quantities.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Frames from two simulated scenarios for $t_L < 0$. Top: $\xi < 0$; bottom: $\xi > 0$. The leader has a red arm and the follower a blue arm. However, to facilitate the distinction among them, a red letter 'L' indicated the leader and a blue 'F' the follower.</figcaption></figure>

$\Delta_m \neq 0$ according to (22). Figure 9 shows the behavior of the system in the described cases through screenshots of the Gazebo environment, and Figure 8 shows the evolution of the main quantities during the simulated tasks.

4) Case of $t_L = 0$: When it comes to the case in which $\xi \neq 0$ and $t_L = 0$, clearly, the sole equilibrium configurations are not really attainable: all elements of the system are supposed to be aligned vertically, one on top of the others (see Figure 3a). When simulating such condition in Gazebo, we found that numerical issues arise as the system approaches the expected configuration in which the link that models the load and those that model the cables are vertically aligned. Despite the practical irrelevance of the considered case, with the objective of demonstrating the validity of the theoretical results, simulations have been carried out also for this case, using the Matlab-Simulink simulator used in [33].

In that simulator, the cables are modeled as mass-less extensible elements and the force is directly retrieved by the model of the cable without resorting to a wrench observer.

<figure><img src="image_11.png" /><figcaption>Fig. 10: Superposition of different instants of simulations with $\xi \neq 0$ and $t_L = 0$. The solid image is the final equilibrium. The grey cylinder is the load, and the red cross is the leader robot. The red cylinder is the desired (identical to the initial) pose of the load.</figcaption></figure>

Nevertheless, underactuated quadrotors are still considered, as well as the same trajectory controller. The results of two simulations can be found in Figure 10. Even though the load
<figure><img src="image_12.png" /><figcaption>Fig. 11: $||e_{pL}||$ for different values of $t_L$ in two cases in which different values of the uncertainties are considered on two parameters, $m$ and $l_{01}$. $e_{pL}$ does not always decrease when $t_L$ increases.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12: Simulation results for $t_L = 1$ N and a 5% error on each parameter. Around Time = 41 s, signed by a red vertical line, the leader robot reads the load position, and corrects its own reference position in order to zero the load position error. Dotted lines of the same color indicate the corresponding desired quantities.</figcaption></figure>

has been initialized in the desired configuration, with position $\bar{p}_L = [0 \ 0 \ 1]^{\top}$ and the same desired yaw and pitch as before, it moves to the vertical equilibrium, with the leader on top when $\xi > 0$, and the follower on top when $\xi < 0$, as explained by the stability analysis in Sec. IV.

5) Position Error: First, we provide in Figure 11 two examples of the different behavior of $e_{pL}$ when $t_L$ increases and different values of the uncertainties are present. This fully supports the finding that the load position error at the equilibrium does not necessarily decrease when $t_L$ is increased.

Anyway, as we have seen from the theory, it is possible to correct the error of the load position by acting solely on the leader robot reference position. This, in turn, does not affect the regulation of the load attitude. In Figure 12, we report the results of a Gazebo simulation in which the initial and desired load pose are as in Sec. VI-1, $t_L = 1$ N, and an error equal to 5% of the nominal value is considered on each uncertain parameter. After 41 s, the leader robot corrects its reference position based on the position of the load according to (43). The results show that, consequently, the load is steered to the desired position when the new equilibrium is reached. On the other hand, as expected from the theory (see Eq. (10)

<figure><img src="image_14.png" /><figcaption>Fig. 13: Picture from the experiment showing: one of the motion capture system cameras 1); two quadrotors 2); two cables 3) connecting the robots to the carbon fiber load 4). The markers on the bar track the pose of the object for validation purposes.</figcaption></figure>

and (11)), due to the inaccurate knowledge of the system parameters, the value of the pitch angle and the leader robot's cable force at the equilibrium do not match the desired values. Also, as predicted, one can observe in Figure 12 that their values are not affected by the change in the leader robot position.

## VII. EXPERIMENTAL VALIDATION

### A. Experimental Setup

1) Hardware: The system is made of a 2-meter-long carbon fiber bar carried by two UAVs by means of two cables that connect the robots at the bar's end. Each cable is 1 m long, the bar weighs 0.300 kg and each UAV weighs 1.03 kg. The cable anchoring points are installed on the robots' underside at a distance $d = [0 \ 0 \ -d]^{\top}$ from their CoM, where $d = 0.15$ cm. Such a geometrical configuration changes the process by which the leader reference position $\hat{p}_{R1}$ is generated as it is explained in the Appendix. In addition, the aerial vehicles have an onboard PC, four ESCs (Electronic Speed Controllers) that control the propeller speed in closed-loop [52], and a flight controller [52].

2) Software: The control architecture runs in part onboard and in part on a desktop PC. A state-of-the-art UKF-based state estimation, which fuses Motion Capture measurements at 120 Hz with the IMU measurements at 1 kHz, and a geometric control are carried out as part of the onboard task at 1 kHz. The admittance filter and wrench observer are implemented in Matlab/Simulink and run on the desktop PC. Wi-fi is used for command and data transfer between the desktop PC and the onboard computers at 100 Hz.

A picture taken from the experiments and highlighting the main setup components is in Fig 13.

3) Experimental Results: Two main sets of experiments were carried out: one in which the controllers use as accurate as possible values of the system parameters; one in which the controllers use values that differ by 10% from the accurate corresponding value. We refer to the former case as 'without uncertainty', and to the latter as 'with uncertainty'.

For each case, we performed 3 tests in which the same manipulation task is carried out for 3 values of $t_L$, equal to 0 N,
<figure><img src="image_15.png" /><figcaption>Fig. 14: Evolution of the pitch (top) and yaw (bottom) angle errors during three experiments with no parameter uncertainties for different values of $t_L$. As expected, the yaw angle converges to the desired value while for the pitch to do the same, $t_L > 0$ is needed. On the left is the boxplot of the average attitude error over the last 20 seconds of the experiments.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15: Evolution of the pitch (top) and yaw (bottom) angle errors during three experiments with 10% uncertainty on each parameter for different values of $t_L$. As expected, the yaw angle converges to the desired value while the pitch angle increases for $t_L = 0$ as the load becomes more and more vertical. For $t_L > 0$, the pitch becomes closer and closer to the desired value. On the left is the boxplot of the average attitude error over the last 20 seconds of the experiments.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16: Experimental results: $t_L = 3$ N and 10% error on each parameter. At Time=100 s, signed by a vertical red stripe, the leader robot corrects its reference position to zero the load position error.</figcaption></figure>

1.5 N, and 3 N. The task execution starts with initial steps in which the load, from position $\mathbf{p}_L(0) = [0 \ 0 \ 0]^\top$ m and zero yaw and pitch angles, is lifted by the robots through simple upwards motions; hence, the proposed controller is activated and the robots try to bring the load to $\bar{\mathbf{p}}_L = [0.5 \ 0 \ 1.5]^\top$ m with $\bar{\psi} = 11.5$ deg and $\bar{\theta} = -6.9$ deg.

The case with no uncertainties is depicted in Figure 14 for all three experiments. The evolution of the attitude error of the load is displayed, in the form of quantities $\theta - \bar{\theta}$ and $\psi - \bar{\psi}$. As expected from (18), the yaw angle at the equilibrium coincides with the desired value. Instead, the pitch angle converges to an arbitrary value when $t_L = 0$, in this case with an error around 18 deg, when $t_L = 0$. When a positive internal force is applied, the attitude error is reduced up to about 2.5 deg.

Figure 15 shows the three tests in the case with uncertainties. Again, in accordance with (18), there is no error on the load yaw angle at the equilibrium. The equilibrium configuration could not be reached with $t_L = 0$. This is also the case in general since so would imply that all bodies are vertically aligned. However, the reader can clearly appreciate the increasing pitch angle evolution, in line with (10). As

soon as $t_L > 0$, the pitch angle at the equilibrium approaches the desired value, as expected by (19). Box plots of the equilibrium error between 130 s and 150 s of the task execution are displayed in Figure 14 and 15 for both cases, i.e. with and without uncertainties, respectively.

The evolution of all the most relevant quantities can be appreciated in Figure 16 for another task execution with uncertainties and $t_L = 3$ N. Furthermore, in that experiment, the leader robot corrects its own reference position at Time=100 s according to (43) and, as a consequence, the load equilibrium position is also adjusted. This validates the load position correction method involving solely the leader robot.

VIII. CONCLUSIONS

This work concerns the decentralized cooperative manipulation of a cable-suspended load by two aerial robots in the absence of direct communication. The robots are controlled with a leader-follower scheme achieved through an admittance controller on each robot. The controllers make use of system parameters that are subject to uncertainty. The equilibrium points and their stability were formally studied. The theory demonstrates how an internal force that stretches the load longitudinally, generated by non-vertically operated cables, is beneficial in terms of stability of the load pose control as well as robustness to the uncertainties. The theoretical results were validated through realistic numerical simulations and experiments.

In the future, an extension to non-beam loads with uncertain parameters will be formally addressed. In this case, $N > 2$ robots will be considered as two robots would not be able to control the full pose of the cable-suspended object: rotations around the line connecting the two cables attaching points on the object would not be controlled [53]. Experimental tests outdoors could be valuable to assess the robustness of the method in windy conditions and when relying on outdoor state-estimation techniques. Investigating how the approach could benefit from limited communication between the robots, e.g., low-frequency communication, is an interesting future
direction. Exploring the possibility of communication-less trajectory tracking and of adaptive laws is left as future work.

APPENDIX

Underneath the robots, at a distance *d* from the CoM (*d* = [0 0 0.15]†cm), the cable anchoring points are attached. A reference position taking into account such a displacement can be provided to the leader robot at the equilibrium according to

$$ \hat{\mathbf{p}}_{R1}^{r} = \hat{\mathbf{p}}_{R1} - \mathbf{R}_{R1}^{r} d \quad (46) $$

where $\mathbf{R}_{R1}^r$ is the leader robot rotational matrix at the equilibrium. This matrix is computed from the leader robot's equilibrium condition under the assumption that the thrust is aligned with the external forces

$$ \mathbf{R}_{R1}^{r} e_{3} = \frac{m_{R} g e_{3} + \hat{\mathbf{f}}_{1}}{\|m_{R} g e_{3} + \hat{\mathbf{f}}_{1}\|} := \begin{bmatrix} A_{1} \\ A_{2} \\ A_{3} \end{bmatrix}. $$

Hence, assuming the yaw is controlled to zero, $\mathbf{R}_{R1}^r = \mathbf{R}_Y(\theta_{R1}^r)\mathbf{R}_X(\phi_{R1}^r)$ with $\theta_{R1}^r = \operatorname{atan}(A1/A3)$ and $\phi_{R1}^r = \operatorname{asin}(-A2)$.

REFERENCES

[1] A. Ollero Baturone, M. Tognon, A. Suárez Fernández-Miranda, D. Lee, and A. Franchi, "Past, present, and future of aerial robotic manipulators," 2021.

[2] F. Ruggiero, V. Lippiello, and A. Ollero, "Aerial manipulation: A literature review," *IEEE Robotics and Automation Letters*, vol. 3, no. 3, pp. 1957-1964, 2018.

[3] H. B. Khamsch, F. Janabi-Sharifi, and A. Abdessameud, "Aerial manipulation—a literature survey," *Robotics and Autonomous Systems*, vol. 107, pp. 221-235, 2018.

[4] M. Tognon and A. Franchi, *Theory and Applications for Control of Aerial Robots in Physical Interaction Through Tethers*. Springer Nature, 2020, vol. 140.

[5] M. Tognon, B. Yüksel, G. Buondonno, and A. Franchi, "Dynamic decentralized control for protocentric aerial manipulators," in *2017 IEEE Int. Conf. on Robotics and Automation*, Singapore, May 2017, pp. 6375-6380.

[6] G. Skorobogatov, C. Barrado, and E. Salamí, "Multiple uav systems: A survey," *Unmanned Systems*, vol. 8, no. 02, pp. 149-169, 2020.

[7] I. Maza, K. Kondak, M. Bernard, and A. Ollero, "Multi-UAV cooperation and control for load transportation and deployment," *Journal of Intelligent & Robotics Systems*, vol. 57, no. 1-4, pp. 417-449, 2010.

[8] A. Mohiuddin, T. Tarek, Y. Zweiri, and D. Gan, "A survey of single and multi-uav aerial manipulation," *Unmanned Systems*, vol. 8, no. 02, pp. 119-147, 2020.

[9] G. Loianno and V. Kumar, "Cooperative transportation using small quadrotors using monocular vision and inertial sensing," *IEEE Robotics and Automation Letters*, vol. 3, no. 2, pp. 680-687, 2017.

[10] H.-N. Nguyen, S. Park, and D. J. Lee, "Aerial tool operation system using quadrotors as rotating thrust generators," in *2015 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*, Hamburg, Germany, Oct. 2015, pp. 1285-1291.

[11] R. Ritz and R. D'Andrea, "Carrying a flexible payload with multiple flying vehicles," in *2013 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*, 2013, pp. 3465-3471.

[12] F. Caccavale, G. Giglio, G. Muscio, and F. Pierri, "Cooperative impedance control for multiple uavs with a robotic arm," in *2015 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*, 2015, pp. 2366-2371.

[13] S. Thapa, H. Bai, and J. Acosta, "Cooperative aerial load transport with force control," *IFAC-PapersOnLine*, vol. 51, no. 12, pp. 38-43, 2018.

[14] K. Sreenath and V. Kumar, "Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots," in *Robotics: Science and Systems*, Berlin, Germany, June 2013.

[15] C. Masone, H. H. Bült Hoff, and P. Stegagno, "Cooperative transportation of a payload using quadrotors: A reconfigurable cable-driven parallel robot," in *2016 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*, Oct 2016, pp. 1623-1630.

[16] M. Manubens, D. Devaurs, L. Ros, and J. Cortés, "Motion planning for 6-D manipulation with aerial towed-cable systems," in *2013 Robotics: Science and Systems*, Berlin, Germany, May 2013.

[17] A. Mohiuddin, Y. Zweiri, T. Taha, and D. Gan, "Energy distribution in dual-uav collaborative transportation through load sharing," *Journal of Mechanisms and Robotics*, 04 2020.

[18] T. Lee, "Geometric control of quadrotor uavs transporting a cable-suspended rigid body," *IEEE Transactions on Control Systems Technology*, vol. 26, no. 1, pp. 255-264, 2017.

[19] G. Li, R. Ge, and G. Loianno, "Cooperative transportation of cable suspended payloads with mavs using monocular vision and inertial sensing," *IEEE Robotics and Automation Letters*, vol. 6, no. 3, pp. 5316-5323, 2021.

[20] D. Sanalitro, H. J. Savino, M. Tognon, J. Cortés, and A. Franchi, "Full- pose manipulation control of a cable-suspended load with multiple uavs under uncertainties," *IEEE Robotics and Automation Letters*, vol. 5, no. 2, pp. 2185-2191, 2020.

[21] F. Rossomando, C. Rosales, J. Gimenez, L. Salinas, C. Soria, M. Sarcinelli-Filho, and R. Carelli, "Aerial load transportation with multiple quadrotors based on a kinematic controller and a neural smc dynamic compensation," *Journal of Intelligent & Robotic Systems*, vol. 100, no. 2, pp. 519-530, 2020.

[22] V. Spurny, M. Petrlík, V. Vonasek, and M. Saska, "Cooperative transport of large objects by a pair of unmanned aerial systems using sampling-based motion planning," in *2019 24th IEEE International Conference on Emerging Technologies and Factory Automation (ETFA)*. IEEE, 2019, pp. 955-962.

[23] P. O. Pereira and D. V. Dimarogonas, "Pose stabilization of a bar tethered to two aerial vehicles," *Automatica*, vol. 112, p. 108695, 2020.

[24] R. C. Sundin, P. Roque, and D. V. Dimarogonas, "Decentralized model predictive control for equilibrium-based collaborative uav bar transportation," in *39th IEEE International Conference on Robotics and Automation (Accepted)*, 2022.

[25] A. Tagliabue, M. Kamel, S. Verling, R. Siegwart, and J. Nieto, "Collaborative transportation using MAVs via passive force control," in *2017 IEEE Int. Conf. on Robotics and Automation*, Singapore, 2016, pp. 5766-5773.

[26] M. Gassner, T. Cieslewski, and D. Scaramuzza, "Dynamic collaboration without communication: Vision-based cable-suspended load transport with two quadrotors," in *2017 IEEE Int. Conf. on Robotics and Automation*, Singapore, May 2017, pp. 5196-5202.

[27] D. K. D. Villa, A. S. Brandão, R. Carelli, and M. Sarcinelli-Filho, "Cooperative load transportation with two quadrotors using adaptive control," *IEEE Access*, vol. 9, pp. 129 148-129 160, 2021.

[28] J. Fink, N. Michael, S. Kim, and V. Kumar, "Planning and control for cooperative manipulation and transportation with aerial robots," in *14th Int. Symp. on Robotics Research*, Lucerne, Switzerland, Sep. 2009.

[29] J. Goodman and L. Colombo, "Geometric control of two quadrotors carrying a rigid rod with elastic cables," *Journal of Nonlinear Science*, vol. 32, no. 5, pp. 1-31, 2022.

[30] D. Mellinger, M. Shomin, N. Michael, and V. Kumar, "Cooperative grasping and transport using multiple quadrotors," in *Int. Symp. on Distributed Autonomous Robotic Systems*, 2013, pp. 545-558.

[31] Z. Wang and M. Schwager, "Force-amplifying n-robot transport system (force-ants) for cooperative planar manipulation without communication," *The International Journal of Robotics Research*, vol. 35, no. 13, pp. 1564-1586, 2016.

[32] A. Tagliabue, M. Kamel, R. Siegwart, and J. Nieto, "Robust collaborative object transportation using multiple mavs," *The International Journal of Robotics Research*, vol. 38, no. 9, pp. 1020-1044, 2019.

[33] M. Tognon, C. Gabellieri, L. Pallottino, and A. Franchi, "Aerial co- manipulation with cables: The role of internal force for equilibria, stability, and passivity," *IEEE Robotics and Automation Letters, Special Issue on Aerial Manipulation*, vol. 3, no. 3, pp. 2577-2583, 2018.

[34] C. Gabellieri, M. Tognon, D. Sanalitro, L. Pallottino, and A. Franchi, "A study on force-based collaboration in swarms," *Swarm Intelligence*, vol. 14, no. 1, pp. 57-82, 2020.

[35] H. Lee, H. Kim, and H. J. Kim, "Planning and control for collision-free cooperative aerial transportation," *IEEE Transactions on Automation Science and Engineering*, vol. 15, no. 1, pp. 189-201, 2016.

[36] H. Lee, H. Kim, W. Kim, and H. J. Kim, "An integrated framework for cooperative aerial manipulators in unknown environments," *IEEE Robotics and Automation Letters*, vol. 3, no. 3, pp. 2307-2314, 2018.
[37] S. Thapa, H. Bai, and J. Á. Acosta, "Cooperative aerial manipulation with decentralized adaptive force-consensus control," *Journal of Intelligent & Robotic Systems*, vol. 97, no. 1, pp. 171–183, 2020.

[38] A. Bicchi, "On the problem of decomposing grasp and manipulation forces in multiple whole-limb manipulation," *Robotics and Autonomous Systems*, vol. 13, no. 2, pp. 127–147, 1994.

[39] M. Kamel, T. Stastny, K. Alexis, and R. Siegwart, "Model predictive control for trajectory tracking of unmanned aerial vehicles using robot operating system," in *Robot operating system (ROS)*. Springer, 2017, pp. 3–39.

[40] M. Ryll, D. Bicego, and A. Franchi, "Modeling and control of FAST-Hex: a fully-actuated by synchronized-tilting hexarotor," in *2016 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems*, Daejeon, South Korea, Oct. 2016, pp. 1689–1694.

[41] K. Nonami, F. Kendoul, S. Suzuki, W. Wang, and D. Nakazawa, *Autonomous flying robots: unmanned aerial vehicles and micro aerial vehicles - Chapter 12*. Springer Science & Business Media, 2010.

[42] F. Ruggiero, V. Lippiello, and A. Ollero, "Aerial manipulation: A literature review," *IEEE Robotics and Automation Letters*, vol. 3, no. 3, pp. 1957–1964, 2018.

[43] J. R. Goodman, J. S. Cely, T. Beckers, and L. J. Colombo, "Geometric control for load transportation with quadrotor uavs by elastic cables," *arXiv preprint arXiv:2111.00777*, 2021.

[44] C. Bisig, J. B. Montejo, M. R. Verbryke, A. Sathyan, and O. Ma, "Genetic fuzzy systems for decentralized, multi-uav cargo handling," in *AIAA Scitech 2020 Forum*, 2020, p. 1117.

[45] A. Yiğit, M. A. Perozo, L. Cuvillon, S. Durand, and J. Gangloff, "Novel omnidirectional aerial manipulator with elastic suspension: Dynamic control and experimental performance assessment," *IEEE Robotics and Automation Letters*, vol. 6, no. 2, pp. 612–619, 2021.

[46] M. Ryll, G. Muscio, F. Pierri, E. Cataldi, G. Antonelli, F. Caccavale, and A. Franchi, "6D physical interaction with a fully actuated aerial robot," in *2017 IEEE Int. Conf. on Robotics and Automation*, Singapore, May 2017, pp. 5190–5195.

[47] R. Horst, P. M. Pardalos, and N. V. Thoai, *Introduction to global optimization*. Springer Science & Business Media, 2000.

[48] H. K. Khalil and J. W. Grizzle, *Nonlinear systems*. Prentice hall Upper Saddle River, NJ, 2002, vol. 3.

[49] A. S. Aghdam, M. B. Menhaj, F. Barazandeh, and F. Abdollahi, "Cooperative load transport with movable load center of mass using multiple quadrotor uavs," in *2016 4th International Conference on Control, Instrumentation, and Automation (ICCIA)*. IEEE, 2016, pp. 23–27.

[50] T. Lee, M. Leoky, and N. H. McClamroch, "Geometric tracking control of a quadrotor UAV on SE(3)," in *49th IEEE Conf. on Decision and Control*, Atlanta, GA, Dec. 2010, pp. 5420–5425.

[51] M. Ryll, G. Muscio, F. Pierri, E. Cataldi, G. Antonelli, F. Caccavale, D. Bicego, and A. Franchi, "6D interaction control with aerial robots: The flying end-effector paradigm," *The International Journal of Robotics Research*, vol. 38, no. 9, pp. 1045–1062, 2019.

[52] A. Franchi and A. Mallet, "Adaptive closed-loop speed control of BLDC motors with applications to multi-rotor aerial vehicles," in *2017 IEEE Int. Conf. on Robotics and Automation*, Singapore, May 2017, pp. 5203–5208.

[53] N. Michael, J. Fink, and V. Kumar, "Cooperative manipulation and transportation with aerial robots," *Autonomous Robots*, vol. 30, no. 1, pp. 73–86, 2011.

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

**Marco Tognon** (S'15-M'18) from November 2022 is a tenuered researcher at INRIA Rennes, France, in the Rainbow team. Before he was a postdoctoral researcher in the Autonomous System Lab (ASL) at ETH Zurich, Switzerland. He received the Ph.D. degree in Robotics in 2018, from INSA Toulouse, France, developing his thesis at LAAS-CNRS. His thesis has been awarded with three prizes. He received the M.Sc. degree in automation engineering in 2014, from the University of Padua, Italy, with a master thesis carried out at MPI for Biological

Cybernetics, Tübingen, Germany. From 2018 to 2020 he has been postdoctoral researcher at LAAS-CNRS, Toulouse, France. He served as Area Chair for RSS2023, and Associate Editor for ICRA and IROS for the years 2019-2020. His current research interests include robotics, aerial physical interaction, multi-robot systems, and human-robot physical interaction.

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

**Dario Sanalitro** received the M.Sc. degree in automation engineering and control of complex systems from the University of Catania, Catania (Italy) in 2017. In 2022, he received the Ph.D. degree in robotics from the National Institute for Applied Sciences, Toulouse, France. He has been a Postdoctoral Researcher at LAAS-CNRS, Toulouse, France in 2022 and he is currently an Assistant Professor at University of Catania, Catania, Italy. His research interests are robotics, aerial physical interactions, multirobot systems.

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

**Chiara Gabellieri** is a researcher in the Robotics and Mechatronics group at the University of Twente. She received the Marie Skłodowska-Curie postdoctoral fellowship with the Flyflic project in 2022 and she leads a Work Package in the European coordination and support project AeroSTREAM. From the University of Pisa, she received her Ph.D. in Information Engineering in 2021, her MSc. in Robotics and Automation Engineering with honors in 2017, and her BSc. in Bioengineering in 2014. She was at LAAS-CNRS in Toulouse, France, from

Nov. 2017 to May 2018 and a visiting Ph.D. student at the German Aerospace Center (DLR) in Oberpfaffenhofen from Nov. 2019 to May 2020. She served as Associate Editor for ICRA in 2022-23 and for IROS in 2023. Her current research interests include environmental robotics and aerial manipulation.

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

**Antonio Franchi** (F'23-SM'16-M'07) received the Ph.D. degree in system engineering from the Sapienza University of Rome, Rome, Italy in 2010, and the HDR degree in Science, from the National Polytechnic Institute of Toulouse in 2016. From 2014 to 2019 he was a tenured CNRS researcher at LAAS-CNRS, Toulouse, France. From 2019 to 2021 he was an Associate Professor and then since 2022 he has been a Full Professor in aerial robotics control, both in the Robotics and Mechatronics group at the University of Twente. Enschede, The Netherlands.

Since 2023 he is also a Full Professor in the Department of Computer, Control and Management Engineering at Sapienza University of Rome. He co-authored more than 160 publications in peer-reviewed international journals, books, and conferences on design, control, and estimation for robotic systems applied to multi-robot systems and aerial robots.