Shared Control of an Aerial Cooperative Transportation System with a Cable-suspended Payload

This paper presents a novel bilateral shared framework for a cooperative aerial transportation and manipulation system composed by a team of micro aerial vehicles with a cable-suspended payload. The human operator is in charge of steering the payload and he/she can also change online the desired shape of the formation of robots. At the same time, an obstacle avoidance algorithm is in charge of avoiding collisions with the static environment. The signals from the user and from the obstacle avoidance are blended together in the trajectory generation module, by means of a tracking controller and a filter called dynamic input boundary (DIB). The DIB filters out the directions of motions that would bring the system too close to singularities, according to a suitable metric. The loop with the user is finally closed with a force feedback that is informative of the mismatch between the operator’s commands and the trajectory of the payload. This feedback intuitively increases the user’s awareness of obstacles or configurations of the system that are close to singularities. The proposed framework is validated by means of realistic hardware-in-the-loop simulations with a person operating the system via a force-feedback haptic interface.


Introduction
The topic of aerial transportation and manipulation has recently gained a lot of consideration in research [39,41,50,54], both thanks to the technological maturity of unmanned aerial vehicles (UAVs) and to the broad range of applications that could benefit from it. Evidence of this is Carlo Masone carlo.masone@iit.it Paolo Stegagno pstegagno@uri.edu 1 the success of several research projects, such as ARCAS 1 , AEROARMS 2 and FlyCrane.
Developing a (semi) autonomous aerial transportation and/or manipulation system requires solving several challenges. Firstly, it is necessary to devise ways for one or more UAVs to interact with objects while considering their limitations in terms of payload capacity and flight time duration. Additionally, the complexity of such an interconnected system requires advanced models and control strategies. Various technical solutions have been proposed to tackle this problem. The first and perhaps most intuitive approach is to equip a single, large sized UAV with a serial arm [23], thus granting it the capability to grasp and move objects. However, this solution has limited applicability as it cannot be used in narrow spaces or indoor. Hence, recent research has been focusing primarily on using micro aerial vehicles (MAVs), multi-rotor systems in particular. Researchers have successfully demonstrated the use of ad-hoc lightweight manipulators mounted on MAVs [22,28,55,64], however the additional mass of these robotic arms drastically reduces the available payload of these vehicles thus compromising their usefulness for transportation.
The limited payload capacity of MAVs can be overcome by using the combined effort of several of them to cooperatively transport a single load [28]. One solution in this regard is to have a rigid formation of MAVs that are physically attached to the load via a special gripper or similar mechanism [32,38,40]. Another approach is to equip the MAVs with a special tooltip and let the UAVs grasp the object by exerting contrasting forces, as in a flying hand [16,47]. However, this last method requires the object to have planar surfaces for the grasp and it allows only for very limited changes in the shape of the formation of robots during operations. A flexible floating transportation structure can be achieved by replacing the rigid links that connect the MAVs to the load with cables (see Fig 1). This idea has some attractive properties: 1. cables are very light thus, allowing to carry heavier objects; 2. the formation of robots is flexible, provided that spherical joints are used to attach the cables; 3. cables-supported transportation does not rely on contact forces to hold the object thus preventing the risk of "losing grasp" of the load; 4. with enough cables/MAVs it is possible to control the full 6-DoF pose of the load.
The scientific community has recently shown a considerable interest in studying these flying cable-suspended systems, focusing in particular on modeling and developing suitable planning and control strategies. However, the question of how to bring a human operator in the control loop to maneuver the system, which is important for the applicability of this system, has only been marginally addressed. As explained in numerous publications [14,26,42,57,61], the concept of human-in-the-loop is pervasive in aerial robotics because UAVs lack the autonomy and intelligence to cope with unpredictable/unplanned situations and to perform complex tasks in possibly challenging environments. The presence of a human operator is a way to endow the system with advanced intelligence and decision making capabilities, until the gap with higher levels of autonomy [48,59] is reduced. Moreover, as mentioned by Dietrich et al. [6], obtaining outdoor autonomous flight authorizations from legal authorities can be problematic. Therefore, including a human operator is also a practical way to simplify the deployment of UAVs in real applications. Despite the advances made in recent years in terms of the autonomy of aerial robots, there is evidence that reliance on a human operator is still fundamental in real-world applications. Skorobogatov et al. [59] published a study in 2020 on the utilization of multiple UAVs in different applications showing that fully autonomy is not yet viable in the real-world. For the specific case of a cable-suspended cooperative aerial transportation system, to the best of our knowledge, so far there is only one study from Prajapati et al. [51] that has focused on including a human operator in-the-loop. In [51], a person acts as a pilot who steers a single MAV while the others move according to a leader-follower paradigm. Yet, this is done without haptic feedback and with very limited autonomy.
In the spirit of assisting the human operator in the task and increasing the autonomy of the system, we propose a novel shared control framework for an aerial cooperative transportation system with cables (see Fig. 1). The proposed framework is designed to let a human operator steer the pose of the carried load with the assistance of a collision avoidance algorithm. The interface also lets the user modify online the shape of the robots' formation. A novel filter, called dynamic input boundary (DIB), is introduced to blend the commands from the user and from the collision avoidance module. Lastly, a force feedback interface is designed to close the loop with the user in order to increase his/her situation awareness. The proposed method leverages the model and kinematic decomposition derived in [35] to separate the motion of the load from the changes in the MAVs' formation. The rest of the paper, is organized as follows. In Section 2 we provide a review of previous works related to our research. In Section 3 we recall the model of the cable-suspended aerial transportation system from [35], with a description of the kinematics and dynamics. Section 4 provides a detailed description of the novel shared control architecture. In Section 5 we present the results of simulations performed with an haptic joystick that demonstrate our solution. Finally, Section 6 concludes the paper with some observations regarding future directions of research beyond this work.

Related Works
In the last decade, the problem of cooperative manipulation of a payload suspended via cables by unmanned multirotor vehicles has been an active topic of research. In [43,44] the authors describe and solve the inverse and forward kinematics of such a system to determine the configurations of the quadrotors that correspond to static equilibria of the payload while respecting the constraint of non-negative cable tensions. The desired motion of the payload is then achieved by navigating through successive static equilibria, assuming that the trajectory is sufficiently slow. This approach is extended in [11] by a planning strategy that, given a series of desired waypoints for the payload, generates a sequence of optimal robot configurations corresponding to static equilibria of the load. Another strategy for cooperative aerial transportation of a cable suspended load is described and experimentally validated with unmanned helicopters in [2,37], but it is based on the assumption that the load is a point mass, i.e., its orientation is not controlled. Full 6-dimensional aerial manipulation of cable-suspended load is addressed in [33] with the FlyCrane, an architecture for aerial towing that uses three quadrotors, each one of them connected to the load via two cables. The authors then adopt a cost-based motion planning that includes constraints on the cable tensions as well as on the thrust of the robots. All the papers mentioned so far either base their results on a quasi-static model of the system or (partially) neglect the dynamics of the payload. The complete dynamics of the payload has been considered in [60] where the authors first prove that the system is differentially flat and then demonstrate with an experiment that considering the full dynamics in the controller improves the tracking accuracy. This approach requires high order derivatives of the load and the specification of the tensions of some cables up to the fourth derivative. However, it is not easy to predict how changes to these tensions affect the trajectories of the MAVs.
Several papers have proposed distributed controllers for this problem. In [34] it is presented a distributed controller based on an incremental nonlinear dynamic inversion that does not require global positioning but rather relative measurements. In [62] the authors propose a masterslave paradigm where the master agent is commanded remotely while the single slave complies with an admittance controller. This approach does not rely on communication nor on the knowledge of the relative positions between agents and load. A leader-follower approach is also used in [15], however this is framed in a more general setting not limited to aerial robots. The role of internal forces in a communicationless cable-suspended aerial transportation system is investigated in [63], albeit only in the case of two robots. The problem of limited communication has also been explored in [52], with a focus on trajectory planning. There the authors combine a global motion planner that generates offline a sequence of waypoints with a local motion planner planner that yields online references for the vehicles without requiring continuous communication.
Robustness is the focus of [29,53,56], which propose controllers that account for perturbations and modeling uncertainties.
More recently, researchers have highlighted the similarity between these aerial cooperative transportation systems (ACTS) and cable driven parallel robots (CDPR). This parallelism was made explicit in [35], which first describes an ACTS as a CDPR with reconfigurable winches and fixed cable lengths, and then introduces a geometric decomposition of the kinematics to provide an interpretation of the internal motions of the system. The authors of [7,9] analyze the wrench capability of an ACTS, expanding previous studies on wrench feasible workspaces for CDPRs by accounting also for the thrust limits of the MAVs propellers. This wrench analysis is combined in [8] with a modified version of a controller developed for discretely reconfigurable CDPRs. The parallelism between ACTS and CDPR is even more pronounced in [30], where the MAVs of the transportation system are equipped with actuated winches.
One facet of cable-suspended cooperative aerial systems that has been only marginally addressed is how to enable a single human operator to maneuver the transported object in an intuitive way. Intuitively, it would be infeasible for a single person to independently steer all the robots at once. The complexity of such a system, where the object to be controlled is actuated by several independent robots via flexible connections, makes it ideal for a shared control strategy that splits the responsibility of the task between one person and an autonomous subsystem. Shared control has been widely researched in aerial robotics [12], because it represents a promising way to accelerate the adoption of UAVs in real world applications where the presence of a human supervisor is often mandatory. Jiang et al. [21] demonstrate an application of shared control in the case of a single UAV, using a hysteresis switch as a blending function. In [26] the authors introduce a collaborative force-feedback scheme in which a person tele-operates a formation of MAVs that is modelled as a deformable objects. Human commands and other commands generated by artificial potentials are summed together to generate the reference velocities to be tracked by the robot. A different shared control architecture for a formation of MAVs is proposed in [13], where the commands of the two subsystems are separated in orthogonal spaces so that they can never interfere. The shared control paradigm is lifted to the planning level in [36] where a human operator is tasked to manipulate the planned path for a single MAV and his/her commands are blended with other assistive commands via a null-space projector. These examples of shared control are however limited to UAVs not physically interacting with objects. In the more challenging case of cable-suspended cooperative aerial transportation, which has emerged only recently, the investigation of how to include the human operator in the task is still in its infancy. So far this problem has only been studied in [51] as a classical teleoperation problem, where the human operator acts as a pilot that maneuvers a single robot of the formation in a leader-follower paradigm, without haptic feedback nor other forms of assistance.
With respect to the state of the art in research we propose a shared control framework tailored for an ACTS with a cable-suspended load. The main characteristics of this architecture are: 1. It contains an obstacle avoidance term to assist the pilot. 2. It exploits the internal motions of the system [35] to separate both the operator's and the autonomously generated commands in two terms: commands that steer the load and commands that modify the formation of MAVs. 3. It introduces the dynamic input boundary (DIB), a novel mechanism for dynamically blending the inputs for the payload based on a measure of the distance of the system from singular configurations. The inputs that affect the shape of the formation are instead mixed by a closed-loop control like approach. 4. A force-feedback is rendered on the haptic device based on the mismatch between the operator's command and the executed motion.
The proposed framework is validated with hardware-in-theloop simulations using an haptic joystick.

System Model
In this paper we consider an aerial cooperative transportation system composed of: 1. a rigid body (payload), tasked to follow a trajectory in a n-dimensional task space, with n ≤ 6 3 ; 2. m ≥ n quadrotors that are connected to the payload via cables (one cable per UAV) 4 .
In this system, as illustrated in Fig. 1, the i-th cable is attached at one end to a point B i on the payload (onboard connection) and at the other end it is attached to a point A i on the quadrotor (moving anchor). Both ends behave like spherical joints, i.e., they are free to rotate. Following the state of the art [7-9, 35, 56] we make two assumptions: The anchor A i is coincident with the center of mass of the i-th quadrotor.

Assumption 2
The cables are massless, inextensible and always taut (straight lines and constant cable length).
With this setting, we introduce the model of the system as described in [35].

Kinematics
In this section we study the differential kinematics of the system and derive a mapping between the motion of the payload and the motion of the anchors (i.e., of the UAVs). This mapping is used in our framework (Section 4.2) to translate the desired trajectory of the payload into a reference trajectory for the quadrotors. Let us represent the pose of the payload with the vector x ν = [p T ν T ] T ∈ SE(3) where p ∈ R 3 is the position and ν = (φ θ ψ) T ∈ SO(3) is the orientation (roll-pitch-yaw angles) of a reference frame F P fixed on the load w.r.t. an inertial world frame F W (see Section 2). The orientation is alternatively represented as the rotation matrix W R P from F P to F W . Velocity and acceleration of the payload are denoted aṡ where ω andω are the angular velocity and acceleration of the rigid body in F W 5 . The position of the end-point of the i-th cable, i.e., the moving anchor A i , is described in F W by the vector a i ∈ R 3 (see Fig. 2). The position vectors of all the anchors can be stacked together in the more compact representation χ = a T 1 · · · a T m T ∈ R 3m . Our goal is to find a mapping that relates x ν ,ẋ andẍ to χ,χ and χ. These are the signals that define the trajectories of the and it is factorized into the unit vector n i ∈ S 2 (cable direction in F W ) and the scalar ρ i > 0 (cable length).

Equation 1
can be written for all the cables in the following compact form where ρ = ρ 1 · · · ρ m T ∈ R m is the vector of cable lengths; χ = a T 1 · · · a T m T ∈ R 3m is the previously defined vector of anchor positions; T is the vector of stacked positions of the connection points B i in F W , with ⊗ being the Kronecker product and 1 m×1 being a vector with m entries all equal to 1; By differentiating Eq. 2 twice we havė Nρ =χ −ξ where we used Assumption 2 to imposeρ =ρ = 0 m×1 . Since the cables directions n 1 , . . . , n m are unit vectors, their derivatives are given by pure rotations. Indicating with ω 1 , . . . ω m the angular velocities and withω 1 , . . . ,ω m the angular accelerations of the single cables, we can expressṄ andN aṡ where W,Ẇ ∈ R 3m×3m are the following block diagonal matrices with [•] × being the skew-symmetric matrix operator.

Remark 1 (Forward Kinematics)
Given the payload pose, the forward kinematics Eq. 2 cannot be uniquely solved because in order to determine the position χ of the anchors it is also necessary to specify the direction of the cables N.
In general, Eq. 2 can be solved either with an optimization algorithm that implicitly chooses N according to some metric or by assigning N. In the architecture presented in Section 4 we only need to solve the forward kinematics to generate the reference trajectory, therefore we can assign the initial value N and then integrate it with the signalsṄ andN that are obtained from Eqs. 6 and 7 by assigning the matrices W andẆ in order to achieve a desired behaviour.
Equations 2, 4 and 5 provide a relation between motion of the anchors (χ,χ andχ) and motion of the load (implicitly expressed by ξ ,ξ andξ ). However, from Eq. 4 and 5 it is not immediately evident how the motion of the anchors relates to the motion of the payload. This explicit relation was provided in [35] through a geometric decomposition of the trajectories Eqs. 4 and 5 of the anchors. Here, we recall that result: Proposition 1 For the cable suspended payload with fully movable anchors and fixed cable lengths, the decompositioṅ χ =χ R +χ K andχ =χ R +χ K of the anchors velocity and accelerations into a component in the m-dimensional range space R(N) of N (denoted by • R ) and a component in the 2m-dimensional null space K(N) where Proof The proof is reported in [35].
This result shows that, given the current state of the kinematic system, the motion of the payload uniquely defines the motion of the anchors on R(N), but not on K(N). The motion of the anchors on K(N) requires to assign also the derivativesṄ andN of the cables configuration. The interpretation of this fact is that only the motion of the anchors on R(N), i.e., along the cables instantaneously affects the motion of the payload whereas the 2m degrees of freedom of the anchors on K(N) can be used to assign internal motions that change the configuration of cables without moving the payload. We will leverage this geometric decomposition as a foundation to build our shared control architecture discussed in Section 4.

Remark 2 (Apparent acceleration)
The term NṄ TṄ ρ in Eqs. 11 and 12 is the apparent acceleration due to the rotation of the cables represented byṄ . Indeed, this term recalls the well known expression of centrifugal acceleration in a rotating frame withṄ as the angular velocity of the frame and Nρ as the position from the center of the frame.

Remark 3 (Connection to CDPRs)
The model presented here is similar to that of a classic cable driven parallel robot (CDPR). In a CDPR with fixed anchors, each cable is associated with only one degree of freedom, the cable length. The relation between motion of the payload and variations in the cable length is given by (see [25]) Clearly, for this system there is only one degree of freedom per cable (cables length) which can only cause an instantaneous motion directed along the direction of the corresponding cable. Indeed, the kinematic relation Eq. 15 resembles Eq. 9 and in fact it is J ≡ N T J R , where the term N T essentially restricts the degrees of χ from three per cable to one per cable.

Dynamics
In this section we describe the dynamics of the system composed of payload and UAVs. We will model these components separately and show their interaction through the cable forces. This model is used by the controller to determine the required cable tensions for the desired payload dynamics and, consequently, the control inputs for the UAVs.

Payload
We consider the payload as a single rigid body with fixed mass and subject to gravity and to the forces exerted by the cables, and we describe it using the classical model commonly adopted to represent the end-effector in cabledriven parallel robots [24,25,45]. The derivation of this model can obtained using the Newton-Euler or Euler-Lagrange methods, as shown in textbooks [10,58], and can be found in previous works in literature [1]. 6 Let us indicate with m P the mass of the payload, with P J P the 3 × 3 inertia matrix w.r.t. F P , and with c P = [c x c y c z ] T ∈ R 3 the displacement vector from the origin of F P to the center of mass of the payload expressed in F W 7 . Additionally, we denote with C I C the inertia matrix of the payload w.r.t. a frame of reference attached that is oriented as F P , but placed on the center of mass of the payload. 8 Using this notation, the dynamics of the rigid payload pulled by the cables, as shown in [24,25,45], is (16) where -B P (x ν )ẍ is dependent on the acceleration and represents the inertia of the system, with -C P (x ν ,ẋ)ẋ describes the centrifugal and Coriolis effects, with g P (x ν ) is the gravitational effects that are dependent only on the configuration of the system, i.e., -J R (N) T Nt is the wrench exerted on the payload by the cables, being t = [t 1 , . . . , t m ] ∈ R m the vector 6 The derivation in [1] is for the more general case of cables with non-negligible mass and variable length. However, these terms can be discarded to fit the case-study considered in this manuscript. 7 This vector is typically expressed in F P where it is constant, i.e., P c P . Moving to F W is straightforward, i.e., c P = W R P P c P . 8 The inertia matrix C I C is transformed to world coordinates as W R P C I C P R W . of tensions of the cables. Substituting the definitions Eqs. 3 and 13 in the matrix N T J R (N), we see that this takes the form which corresponds to the wrench matrix that maps the cable tensions to wrenches applied on the end-effector in cable-driven parallel robots [18].
Note that in the case in which the origin of the frame F P is coincident with the center of mass of the payload, the vector c p is identically null and thus the expressions Eqs. 17 to 20 are simplified.
One important remark about the model Eq. 16 is that the cables cannot push the payload, but only pull it. Namely, the tensions t 1 , . . . , t m must be positive. More in general, following the classic formulation from cable-robots [24], we assume the following constraint where the upper bound t depends on the physical properties of the cables and actuators, and the lower bound t is chosen high enough to prevent cable slackness in accordance to Assumption 2.

Quadrotor UAV
To model the (underactuated) 6 DoF dynamics of the i-th UAV connected to the cable we consider a north-west-up body frame F Q i that is attached to the center of mass A i of the robot. The dynamics of the i-th UAV depends on its mass m i and on the diagonal inertia matrix J i w.r.t. the body frame, and has the well known form [26] where -Q i ω Q i is the angular velocity of the UAV in body frame; -W R Q i is the rotation matrix of the body frame w.r.t. F W and e 3 = [0 0 1] T ; -−t i n i is the reaction force applied by the cable on the UAV's center mass; -J i is the inertia matrix of the UAV w.r.t. F Q i ; -τ i ∈ R is the thrust control input in body frame and ζ i ∈ R 3 is the attitude torque control input.
Remark 4 (Reaction from the cable) We can divide the translational dynamics of the i-th UAV into two components, parallel and orthogonal to the cable direction n i . With the same approach used in Section 3.1, applying the projector operators π R = n i n T i and π K = (I 3 − n i n T i ) to Eqs. 23 and 24 gives Equation 25 shows that the cable tension enters the translational dynamics of the UAV only along the cable direction. Namely, only the motion of the quadrotor along the cable applies a wrench on the load, in accordance with the spatial decomposition of the kinematics. This effect is related to Assumption 1. In practice, the anchor may A i may not be coincident with the center of mass of the i-th UAV and in this case the reaction from the cable would also apply a torque on the UAV. In the results presented in Section 5 we demonstrate the effectiveness controller derived using the model Eq. 24 even when Assumption 1 is not fulfilled.

Shared Control Architecture
The shared control architecture that we have designed for the aerial cooperative transportation system is based on the model described in Section 3 and it follows three principles: -There is a single human operator who is in charge of steering the load, i.e., commanding translations and rotations. This commands are given without the person having to directly control any of the MAVs. However, we want to give the person the ability to modify the shape of the formation, if needed. Different reshaping motions are possible but, for the sake of keeping the user interface simple, we consider a single reshaping function that coordinately rotates the cables upwards (formation shrinking) or downwards (formation blossoming) relatively to the load. Even though this form of reshaping is simple, in practice it can be used to make the formation of robots more compact when traversing a narrow passage, or more open in order to increase the cable forces for better stability. -Automation assists the human operator in the form of an obstacle avoidance module based on artificial potentials. This module is meant to help particularly in the phases that require precise navigation and positioning in narrow passages or amidst static objects. This choice is motivated by the consideration that such a transportation system is not expected to operate with high dynamics and to avoid fast moving objects. -A force feedback increases the situational awareness of the operator by encoding how the system is tracking her/his commands. For this purpose, the input device used by the operator must be actuated to be able to exert forces.
The shared control architecture that we propose is composed of four subsystems (see Fig. 3): -Inputs Provider -this subsystem produces the inputs that eventually generate the trajectories for the payload and MAVs. Following the criteria previously established, the inputs come from two sources: i) the human operator, who steers the load and can change the desired shape of the formation, and ii) an obstacle avoidance algorithm, that is responsible for avoiding collisions with the static environment. -Reference Trajectory Generator -this subsystem blends the inputs from the user and from the collision avoidance, and generates the reference trajectories for the payload and for the UAVs. This mapping is based on the kinematic model from Section 3.1 and leverages the geometric decomposition presented therein. -Closed Loop Controller -this subsystem generates the control inputs for the UAVs with the goal of achieving accurate trajectory tracking for both the load and the robots. The controller is based on the dynamic model from Section 3.2. -Haptic Controller -this subsystem closes the loop with the user by rendering on the input device a force feedback that encodes the mismatch between her/his input and how they are executed by the system. This feedback is not only based on the motion of the payload, but it also encodes how the MAVs are tracking the desired formation.

Inputs Provider
There are two sources that, independently from each other, provide the inputs for the trajectory generation subsystem.
Here we describe these sources and the signals they output.

Human Inputs
The primary role of the user is to steer the payload. For this purpose, we map the axes of an actuated input device (see e.g., Fig. 9b) to the signalẋ ν,hmn = [ṗ T hmnν T hmn ] T ∈ R 6 , whereṗ hmn is a translational velocity andν hmn is a rotational velocity expressed as roll-pitch-yaw rates. We map the degrees of freedom (DoF) of the input device to a velocity because this does not limit the workspace and it is easier to use than an acceleration mapping. Additionally, it is more intuitive for a person to control the orientation of a system in terms of roll-pitch-yaw angles rather than commanding angular velocities. Nevertheless, the conversion from the roll-pitch-yaw ratesν hmn to an angular velocity in ω hmn in F W is trivial, therefore when needed we will use the signalẋ The user is given the freedom to control the vertical angles of the cables w.r.t. the payload (elevation angles), thus effectively changing the shape of the formation of MAVs. Formally, let us express the direction of the generic i-th cable in the frame F P by using polar coordinates, i.e., where P η i and P α i are the elevation and azimuth angles respectively. Going from this representation in the local frame F P to the world frame is a matter of a rotation, i.e., P n i = P R W n i . The second input given by the human is then an angular rate Pη hmn that is broadcast to all the cables, i.e., Pη hmn,1 = . . . = Pη hmn,m = Pη hmn .

Collision Avoidance Input
For the collision avoidance input we assume that the load/cables/MAVs system is approximated by a finite number of simple shapes (see Fig. 4a) and that there is an obstacle provider that supplies the distances and relative directions of these shapes w.r.t. to any nearby obstacle. We do not investigate how this information is estimated because this is not the focus of the manuscript and there is a rich body of research on this topic, e.g., [17,49]. Example of the approximation of the system as a collection of simple shapes. The payload is approximated by the green sphere, the −i-th cable is approximated by the blue spheres and the i-th UAV is approximated by the red sphere. b) Example of a distance based potential function that is zero when δ jk ≥ δ max and approaches infinity when δ jk approaches δ min from the right side Using distance-based artificial potentials we generate the obstacle avoidance inputs: a translational velocityṗ obs that acts on the payload, and a set of angular velocities ω obs,1 , . . . , ω obs,m that act on the cables and UAVs. Formally, these signals are computed aṡ where -P, C i and O are the set of shapes approximating payload, cables (including the UAVs) and obstacles, respectively; -δ jk ∈ R ≥0 and d jk ∈ S 2 are the distance and and unit direction from the shape S j to the obstacle O k ; -U : R + → R + is the artificial potential (see Fig. 4b for an example); -l j is the distance from the center of the shape S j to the connection point B i of the i-th cable (see Fig. 4a).
Intuitively, the reactive term Eq. 28 acts on the payload, changing its translation speed to avoid collisions with the environment. On the other hand, the reactive terms Eq. 28 act individually on each cable/MAV, modifying the corresponding orientation and thus changing the overall configuration of the cables, without instantaneously affecting the payload.

Reference Trajectory Generator
The role of this subsystem is to take both the human and the collision avoidance inputs and blend them to generate the reference trajectories for the payload and for the quadrotors.
As represented in Fig. 5, there are two generators for these trajectories: one for the internal motions and one for the motions that produce the desired trajectory of the payload.

Internal Motions Generator
The The term W obs is directly obtained by plugging the signals ω obs,1 , . . . , ω obs,m in Eq. 8. Similarly, the termẆ obs is obtained from the signalsω obs,1 , . . . ,ω obs,m which are not readily available but can be computed through numerical differentiation and filtering. The two terms W hmn andẆ hmn based on the user's input are less straightforward. Firstly, the signal Pη hmn,i is filtered to produce the smooth reference signals P η des,i , Pη des,i and Pη des,i . Filtering a possibly discontinuous and unbounded signal to produce a smooth and bounded trajectory is a common procedure often referred to as input shaping. Many such methods exist in literature, and we refer the interested reader to Appendix A for a brief description of the specific method we adopted in our experiments. For uniformity of notation we introduce a similar notation for the desired azimuth angles, even though they are not modified by the user. Therefore we have P α des,i = const and Pα des,i = Pα des,i = 0. Lastly, using well known changes of coordinates we express the desired azimuth and elevation angles, velocities and accelerations from the frame F P to the world frame F W . These new signals are denoted as η des,i ,η des,i andη des,i , for the elevation, and α des,i , α des,i andα des,i for the azimuth. Note that the even though the azimuth angles P α des,i in F P are constant, the angles α des,i expressed in F W will change as the orientation of the payload changes.
With this setting, the signals ω hmn,i andω hmn,i , with i = 1, . . . , m are computed as the output of a configuration tracking controller

Singularities
In Section 4.2.1 it was discussed how the human operators can command internal motions, i.e., motions that only change the configuration of the cables without directly affecting the motion of the payload. However, it must be noted that these internal motions may have an indirect effect on the payload as well. This has to do with the concept of singular configurations: Definition 1 (Singularity) A singularity for the aerial transportation system is a configuration N of the cables in which the dimension of the set of feasible wrenches The reduction in the set of feasible wrenches mentioned in the definition of singularity may be due to two different factors: -The cables are oriented in such a way as to make the wrench matrix N T J R (N) singular. For example, if all the cables are vertical, no horizontal force can be applied on the payload. -Even when the wrench matrix N T J R (N) is not singular, singularities can occur due to constraint Eq. 22, which limits the range of viable cable tensions. Figure 6 shows an example of a singular configuration in which N T J R (N) is full rank.
It is important to prevent the system from reaching a singular configuration because this would limit the wrench that can be exerted on the payload via the cable forces. In particular, the internal motions may lead to a singularity in two different ways: The elevation command from the user drives the system to a singularity Depending on the initial configuration of the cables, the user can potentially steer the system to a singularity with the elevation rate command Pη hmn . To illustrate this possibility, consider for simplicity the example illustrated in Fig. 7. The system starts in a configuration with the elevation angles of the cables all set to the same value, specifically η i = 40 • (Fig. 7 left). The user, commanding Pη hmn could drive these angles to be η i ≈ 0 • (i.e., the cables are near horizontal, see Fig. 7 bottom right). This is a singular configuration because it is not possible to exert vertical forces on the payload. Similarly, the operator could steer the elevation angles to become η i ≈ 90 • (i.e., the cables are near vertical, see Fig. 7 top right), in which case it would not be possible to exert lateral forces on the payload. In both cases, the wrench matrix N T J R (N) would lose rank. However, this problem can be prevented by limiting the range of elevation angles that can be controlled by the user, so that the manifold of reachable configurations does not contain singularities. With respect to the example shown in Fig. 7, it would be reasonable to impose 0 < η i ≤ η i ≤ η i < π/2, with η i and η i suitably chosen. This is the solution adopted in the results presented in Section 5 and it is implemented via the input shaping method described in Appendix A.
The obstacle avoidance acting on the cables drives the system to a singularity As the operator steers the system near obstacles, the obstacle avoidance reaction on the cables may lead the system to a singularity. To understand how this could happen, consider the example illustrated in Fig. 8, which depicts a transportation system with 4 robots and the task is to control the payload's position (n=3). Assume that the operator is steering the payload towards a narrow gap in a wall (Fig. 8a). The configuration of the cables in Fig. 8a is the one desired by the user. As the robots move, the cables closer to the wall, i.e., the ones pulling the object, are pushed away from the obstacle by the artificial forces (Fig. 8b). As a consequence, the actual force that can be exerted in the direction of the wall (red arrow) is reduced w.r.t. the desired configuration (the one shown in Fig. 8a). If the human pilot keeps steering the object to the right eventually the system reaches a singularity in which the cables pulling the object become vertical (Fig. 8c). In this situation, even though the wrench matrix J R (N) T N has full rank, it is not possible to exert on the payload forces towards the right because there are no cables that can pull in that direction. Preventing this kind of situation where the obstacle avoidance causes to reach a singularity is not trivial. On one hand, limiting the action of the obstacle avoidance with a maximum deviation of a cable's orientation from its desired state would hamper the ability to avoid collisions. On the other hand, although this mechanism is caused by the obstacle avoidance, it is the operator's commands to the payload that push the system close to the static obstacles in the environment. Therefore, preventing this phenomenon requires a solution that takes into consideration both the operator's commands to the payload and the obstacle avoidance action on the cables. We present such a solution in the next section, in the form of a filter that blends the inputs from the operator and obstacle avoidance.  (N) is full rank, the system is at a singularity. The black arrows indicate the direction of the forces exerted on the payload. There is no force that can pull the payload to the right, because the two quadrotor on the right are aligned vertically Fig. 7 Example of how the operator command Pη hmn could drive the system to a singular configuration. Left) The system starts with all the elevation angles η i = 40 • . Top-Right) The operator drives the cables to become nearly vertical, i.e., η i = 90 • . In this case it is not possible to exert lateral forces on the payload. Bottom-Right) The operator drives the cables to become nearly horizontal, i.e., η i = 0 • . In this case it is not possible to exert vertical forces on the payload

Payload Trajectory Generator
The payload trajectory generator takes the signalsẋ hmn anḋ p obs and generates as output the reference trajectory for the payload. Its core component, as shown in Fig. 5, is the dynamic input boundary (DIB), a system that blends the two inputs and filters out their contribution along the directions where the achievable force (or torque) is too small w.r.t. the one expected according to the desired configuration N des . This heuristic stems from the observation made in the previous section, that these directions are those where the system is close to a singular configuration (see example in Fig. 8). Hence, by blocking the commands along those directions the DIB prevents the user from driving the system to a singularity. The second purpose of the DIB is to allow the haptic feedback to capture changes in the cables' orientation caused by the obstacles. This point will be discussed in Section 4.4.
In order to formally describe the DIB, let us define the direction of the translational and rotational inputs as We also introduce two functions that are based on the dynamic model Eq. 16 and describe the maximum force and torque The red arrow indicates the maximum force that can be exerted on the payload towards right. a) The user is steering the payload to the right, towards a passage between some obstacles. b) The obstacle avoidance starts acting by pushing away the two closest MAVs. c) The two MAVs on the right have been pushed so much that their corresponding cables are almost vertical (near singular configuration). In this state it is not possible to exert on the payload a force oriented towards the right that can be exerted on the payload along a certain direction d ∈ R 3 and in a certain configuration N ∈ S 2 m , i.e.
With these definitions, the DIB's translational and rotational outputs are defined aṡ and f (d trn , N des ) = 0 p hmn +ṗ obs otherwise (33) and τ (d rot , N des ) = 0 ω hmn otherwise (34) with trn , rot ∈ (0, 1). Together, Eqs. 33 and 34 implement a novel mechanism that blends the commands given to the payload by the human operator and the obstacle avoidance. This mechanism takes into account the changes of configuration caused by the obstacle avoidance to filter out motions that could lead to singularities and which are identified by a reduction in the attainable force or torque.
Finally, the piecewise continuous output of the DIB is filtered with an input shaping method to produce the smooth reference signals

Closed Loop Controller
The controller takes the reference trajectories (x ν,ref ,

ẍ ref ) and (χ ref ,χ ref ,χ ref )
and computes the commands for the UAVs. To tackle this problem, we adopt the dual-space control approach with tension distribution presented in [35]. This controller is composed by two elements: i) a task space controller that computes the wrench to be applied on the payload and maps it to a feasible set of cable tensions; i) a set of joint space controllers that compute the inputs for the quadrotors subject to the required forces from the cables.

Task space controller
The task space controller is implemented as a classic closedloop inverse dynamics control, where K 1 , K 2 are suitable diagonal matrices of positive gains. The vector f r ∈ R 6 is the desired wrench applied to the payload. Following the dynamic model derived in Section 3.2.1, applying this desired wrench on the payload means finding a vector of cable tensions t r such that Provided that the system is not in a singularity, thanks to the DIB, and that the required wrench is feasible, then the required vector of cable tensions t r is obtained by inverting Eq. 36. In practice this assumption can be guaranteed by restricting the payload and cables configurations to the socalled Wrench-Feasible-Workspace (WFW) [4], i.e., the set of configurations in which, for any wrench in a desired set there is a tensions vector that solves Eq. 36 and satisfies Eq. 22. Furthermore, provided that the required set of wrenches contains a neighborhood of the origin, the wrench matrix J T R N has full rank [19]. If the system is not redundant w.r.t. the task, i.e., m = n, then if the desired wrench is feasible there is only a single solution to the problem Eq. 36. However, in practice it is much more common to assume that the number of cables is greater than the dimension of the task space, i.e., m > n, because: i) using more robots allows to increase the weight that can be carried, and ii) having more cables gives more flexibility to change the configuration of the system and extends the wrench feasible workspace. In the redundant case, the solution to relation Eq. 36 is not unique but it lies in the m − n dimensional null space of the wrench matrix J T R N. One specific solution would be the one obtained by the pseudo-inverse of the wrench matrix, i.e., Solution Eq. 37 corresponds to the vector of cable tensions with the minimum Euclidean norm, but this is not guaranteed to satisfy the constraints Eq. 22 even if the wrench f r is feasible. In fact, some of the cable tensions found by Eq. 37 may be even negative. In the redundant case, instead of using Eq. 37 we look for a solution of the following optimization problem The parameter α is used to modulate the kind of solution that we are looking for. When α = 1 we try to find a solution to Eq. 36 that is as close as possible to the maximum cable tension, which could lead to a more rigid formation, whereas when α = 0 we want the solution to be as close as possible to t, which could help reducing the effort of the MAVs. In our simulations we always set α = 0.5, which means that we look for the solution to Eq. 36 that is furthest away from the tension limits. It is important to remark that the solution of this problem exists if the wrench w r is feasible. As stated beforehand, in practice this is verified if the range of configurations achievable by the cables are chosen so that the Wrench-Feasible-Workspace contains the range of wrenches expected for the application. In practice, if a feasible solution is not found we keep the last valid vector of cable tensions, yet this situation never occurred in our simulations.
Since the optimization problem Eqs. 38 to 42 is well studied in the literature of cable-driven parallel robots (under the name tension distribution) and since this specific problem is not the focus of this framework, we do not go into the details of how to solve it. One of the many methods existing in literature can be used, e.g., [24,46]. Nevertheless, the interested reader can find in Appendix B a few more details on how in practice we reformulate the optimization problem Eqs. 38 to 42 so that it can be solved more efficiently.

Joint Space Controller
We implement the tracking controller for quadrotors that is detailed in [26]. This controller has an inner/outer loop structure. The slower outer loop position controller computes the thrust input and determines the desired roll and pitch commands as The faster inner loop attitude controller computes the torque input as where ν i,r are the desired roll-pitch-yaw angles, k d , k p are positive gains and E(ν) is the well known matrix that gives the mappingν = E(ν) Q i ω Q i . The stability of the controller is proven in [26] for near-hovering configurations.

Haptic Controller
The operator to commandẋ hmn = [ṗ T hmn ω T hmn ] T using an actuated input device, i.e., a joystick that can exert forces through its end-effector on the hand of the operator. Assuming gravity compensation, the model of such an actuated input device is where q ∈ R 6 is the configuration vector representing the 3D position and orientation of the tip of its end effector, B H D (q) ∈ R 6×6 is the positive-definite/symmetric inertia matrix, C H D (q,q) ∈ R 6×6 is the Coriolis matrix, and τ H D , f H D ∈ R 6 are the control and human wrenches, respectively. The configuration of the input device is mapped to the commandẋ h by direct scalinġ where λ H D ∈ R 6×6 is a diagonal matrix of positive scaling factors.
The purpose of the haptic feedback interface is to render the command force τ H D in order to provide an information of how the payload is following the operator's command. Namely, we desire the force feedback to be representative of the mismatch between the executed velocity and the one specified by the operator (ẋ hmn ), i.e., with k 1 ≥ 0 and k 2 > 0. The command force is composed by a damping term (it can be made very small, to improve transparency) and by the information term. It is well known from the literature that the implementation of such a control action can be particularly challenging due to problems such as delays and packet losses. It is not our interest to study the implementation of Eq. 47, noting that there are several state of the art approaches to do it. We used the passive setposition modulation algorithm which ensures passivity of the closed teleoperation system and we refer the interested reader to [27] for further details. We must point out that the physical device used in our simulations (see Fig. 9b) only has three actuated degrees of freedom. Due to this practical limitation we allowed the operator to control alternatively either the position or orientation of the payload, receiving back the corresponding feedback.
Assuming that the payload tracks accurately the reference trajectory, the main information conveyed by the feedback Eq. 47 is the mismatch between the user commands and the tracked reference. In order to give an interpretation to this mismatch we must discuss its causes. From the design of the command interface in Section 4.2, there are three causes for the mismatch: 1. The commands given by the human are passed through shaping filters that cause slower transients in the reference trajectory and such difference is captured by the force τ H D . For example, if the operator commands a sudden change of direction to the load he/she would feel a resistance on the input device that will disappear as the reference velocity catches up with the commanḋ x hmn . This kind of reaction provides a very intuitive and simple means to inform the pilot of the limits of the system and can be interpreted as a 'virtual inertia'. 2. Another source for the mismatch is the obstacle reactionṗ obs that pushes the payload away from nearby obstacles with an intensity that increases monotonically as the distance decreases (see Fig. 4b). The operator would therefore feel this reaction as a resistance on the input device whenever he/she is driving the payload towards a close enough obstacle. 3. The DIB can be a cause of the mismatch. This mechanism, defined by Eqs. 33 and 34, cancels motion commands that would lead the system in a singular configuration. To understand its effect in terms of force feedback we refer the reader once again to the example illustrated in Fig. 8. The operator is commanding the system between some obstacles and whenever the reconfiguration caused by the obstacles exceeds a predefined limit the input from the user is cancelled. When this happens the translational mismatch is exactly −ṗ hmn and the operator feels an opposing force, like a virtual wall, pushing away from the obstacle. This force informs the user that this direction of motion is infeasible because it leads the system to a singularity. Namely, even though the force feedback is only computed on the mismatch in the payload trajectory, the presence of the DIB allows it to be informative also of the obstacle avoidance reactions applied on the MAVs/cables that bring the system close to a singularity. Note that the piecewise continuous velocityẋ DI B is passed through an input shaping filter that yields a smooth reference. Therefore the feedback caused by the DIB is continuous. The effect of the DIB is similar to the dynamic kinesthetic boundary (DKB) proposed in [20]. Both solutions allow normal operations in absence of obstacles and limit the inputs of the user when approaching the environment, but there are some conceptual differences between the two approaches. On one hand, the DKB projects the environment into hard boundaries of the input device workspace and, as a consequence, it limits the inputs of the user. On the other hand, the DIB uses an opposite perspective because it cancels out the inputs and, as a consequence, it produces a feedback on the actuated input device. Moreover the action of the DIB is not based directly on the distance from the obstacles but rather on the configuration of the system.

Simulations
The proposed shared control architecture has been tested in simulation. All the parts of our shared control framework have been implemented in Simulink, using standard block and Matlab functions. On the other hand, the physical system is not simulated by numerically evaluating the kinematic and dynamic models derived in Section 3. Rather, we used SimScape's Multibody simulation environment for 3D mechanical systems. In SimScape, mechanical systems are created by connecting together masses with joints and other elements. At each time step in Simulink simulation, SimScape Multibody formulates and solves the equations of motion for the complete mechanical system by applying the external forces and torques (including gravity) and imposing the constraints originated from the connection of its parts. The cooperative transportation system in SimScape consists of three main elements: -Payload: It is implemented as a single rigid body, namely a box with uniform density and weighing 1.5 kg. In the controller we use the exact values of mass and inertia matrix of the payload that are measured using SimScape's inertia sensor. -MAVs: The physical model of the robots is available from MATLAB Central file exchange 9 [5] and it is created in CAD (see Fig. 9a). As for the payload, the mass and inertia matrix of the MAV were measured using SimScape's inertia sensor. However, in the controller we utilize the exact value fore the mass (0.7368 kg) but not for inertia matrix. There are two main reasons for this choice. Firstly, in practice it is not hard to measure the mass of a quadrotor using a scale, but estimating its inertia matrix is less trivial and more prone to errors. Secondly, the frame F Q i that is used in Section 3.2.2 to model the dynamics of the MAV is supposed to be co-planar with the propellers and coincident with the principal frame of inertia, however this is not true in the SimScape model. In practice, the measured principal inertia matrix of the UAV is J p i = diag(0.005095, 0.005217, 0.009955) whereas the one we used in the controller is J p i = diag(0.005131, 0.005185, 0.009955), where small off-diagonal terms have been set to zero. -Cables: SimScape provides out-of-the-box an element that models an ideal cable, without mass and inextensible, that is fixed at two points by means of spherical joints. These two points are one on the payload and the other on the bottom of the MAV (anchor). Note that the anchor does not coincide with the center of mass of the MAV (see Fig. 9a) as it is hypothesized in Assumption 1, thus making the simulation more realistic.
The loop with the shared control algorithm in Simulink is closed by inserting SimScape sensors that measure the state of the robots and payload, with Gaussian noise added to the measurements. Concerning the controller, we empirically set the minimum and maximum tensions that can be exerted by each cable to 0.1 N and 10 N, respectively. With this setup we perform two sets of simulations. In the first simulation we focus on assessing the behaviour of the proposed shared control with the human-in-the-loop. The input device used by the operator is a SensAble (Geomagic) Phantom Omni haptic device (see Fig. 9b). This device has six continuous axes, but only three are actuated (the wrist axes are passive), hence we performed two different simulations with the human operator in-the-loop. In the first one, the three actuated axes are used by the operator to control the position of the payload (signalṗ hmn ) while its desired orientation remains fixed (i.e.ν hmn ≡ 0 3×1 ). In the second one, the desired position of the payload is fixed (i.e. p hmn ≡ 0 3×1 ) and the three actuated axes are mapped to rotations of the payload (signalν hmn ). In both simulations the operator is allowed to change the configuration of the cables by using the two buttons on the pen tooltip that are mapped to the fixed rotation rates Pη hmn = ±0.1 rad/s. In the second simulation we investigate the robustness of the system to external forces that could represent wind or other disturbances acting on the system. In this second simulation the operator gives no commands and the payload is tasked to remain in a fixed pose in order to better evaluate the impact of the disturbances. We encourage the reader to watch the videos of the simulations before reading the detailed discussion of the results presented hereinafter.

Translation Commands
This simulation mimics a transportation task in which the operator has to pilot the payload from the initial position p = (15, 0, 4) [m] to a goal location defined as a box centered in p = (−15, 5, 1) [m] while keeping the orientation fixed to ν = (0 • , 0 • , 0 • ). In order to reach the goal the payload must be carried through a narrow passage between two walls. This non trivial setting allows to demonstrate the combined effect of user's commands, obstacle avoidance and DIB. Figure 10 shows four snapshots from the simulation. Starting from the initial configuration (Fig. 10a) the user drives the payload towards the narrow passage where the obstacle avoidance starts changing the formation (Fig. 10b). Eventually, the DIB stops the motion of the system and the operator uses the reconfiguration command Pη hmn to group the robots closer, enter the passage (Fig. 10c) and finally reach the goal (Fig. 10d).
Let us now inspect how the reference trajectory generation worked during the execution of the task. The inputs for the payload from the userṗ hmn and from the obstacle avoidanceṗ obs are shown in Fig. 11a and Fig. 11b, respectively. The plots show that the user commands translations along all the three Cartesian axes. Note also thatṗ obs is not null only at around 40s, which corresponds to the time when the payload is transported over the transversal obstacle in the middle of the passage as shown in Fig. 10c. The total inpuṫ p hmn +ṗ obs (Fig. 11c), is then passed through the DIB. The DIB operates by blocking the inputs when the condition f (d trn , N ref )/f (d trn , N des ) ≤ trn = 0.92 is verified. We see in Fig. 11d that the ratio f (d trn , N ref )/f (d trn , N des ) (blue line) remains always above trn (red line) except for a short interval approximately between 17s and 21s. This is the moment when the user first tries unsuccessfully to push the payload through the passage, because the formation is to large to pass in the gap and it is blocked by the collision avoidance. During this interval the DIB filters the inputs and the signalṗ DI B is set to zero (see Fig. 11e). Finally, the input shaping method filtersṗ DI B producing the smoother reference velocityṗ ref (see Fig. 11f). The other side of the trajectory generation concerns the internal motions of the system. Figure 11a shows the elevation command Pη hmn given by the user. We see that the user increases the elevation angle at around 20s, when the system is blocked at the entrance of the passage (situation shown in Fig. 10b). By doing so, the footprint of the formation reduces to the point that it fits the passage, the obstacle avoidance effect vanishes and the DIB filtering action stops. Afterwards, when the formation successfully crosses the passage (around 55s), the user decreases again the elevation angles. Besides the user input, the internal motions of the formation are determined also by the obstacle avoidance acting on the cables and MAVs. To understand how this obstacle avoidance affected the system we look at the elevation and azimuth angles of the formation (see Fig. 11h and 11i). From Fig. 11i we see that when the DIB is active (between 17s and 20s) the reference azimuth angles for two cables deviate from the user given desired values. This deviation is the effect of the obstacle avoidance that pushes two of the MAVs away from the walls in the environment (Fig. 10b). Similarly, we observe deviations in the azimuth and elevation angles of four cables during the traversing of the narrow passage (Fig. 10c) due to the obstacle avoidance reacting to the lateral walls.
In terms of task execution, the load tracks the the reference position accurately (see Fig. 12a) while keeping the orientation error below 0.3 • in roll, pitch and yaw (see Fig. 12b). We report also the tracking performance of one of the quadrotors (they are all comparable). Without loss of generality we select the first UAV and Fig. 12c shows that the reference trajectory is followed accurately. Finally, we report the attitude of the UAVs during the execution of the task in Fig. 12d to 12f. Now we inspect the force feedback produced by the haptic controller to close the loop with the user. This force feedback, shown in Fig. 13a, is informative of the mismatchṗ hmn −ṗ. In particular, by looking at the signals in the trajectory generation pipeline (Fig. 11), we can observe that the force feedback has three main sources (or equivalently gives three main pieces of information): 1. The DIB, which generates an intense cue between 17s and 21s. This mechanism allows the user to perceive the obstacles that deform the formation 'too much'. In particular, we can see that while the user is pushing in the negative x direction (see Fig. 11a), the force is opposite. This also demonstrates how the DIB, besides impeding the motion of the system towards singularities, further helps preventing reaching them by giving a force feedback that pushes the user away from the obstacles and, consequently, away from the singular configuration. 2. The inputṗ obs , which generates a spike at around 40s.
This informs the user of the obstacles too close to the payload. In particular, the obstacle is below the payload and the force feedback is then in the opposite direction. 3. The input shaping, which produces smaller cues because the filtering ofṗ DI B introduces slower transients in the Finally, we can verify that the additional objectives concerning distance from the obstacles and cable forces have been satisfied throughout the execution of the task. Firstly, the minimum distance between the cables/robots and the obstacles is always above the minimum value 0.3m (see Fig. 13b). This value was arbitrarily chosen for the simulation, but it can be increased/reduced by modifying the artificial potential functions that generate the inputs of the obstacle avoidance. Secondly, the forces measured at the cables always satisfy the constraint 0.1N ≤ t i ≤ 10N (see Fig. 13c).

Rotational commands
This simulation mimics a manipulation task in which the operator is asked to rotate the payload from the initial orientation ν = 0 3×1 to the goal orientation ν ≈ (12 • , 12 • , 44 • ) while keeping its position fixed to p = (0, 0 4) [m]. The task is not trivial due the presence  Figure 14 shows four snapshots from the simulation. Starting from the initial pose (Fig. 14a) commands a counter-clockwise yaw rotation to reach the goal yaw, however the execution of this command drives some cables/robots towards the columns. (Fig. 14b). Eventually, the DIB stops the motion of the system and the operator uses the reconfiguration command Pη hmn to group the robots tighter and reprise the yaw rotation (Fig. 10c). Finally, the user also drives the payload to the target roll and pitch orientation (see Fig. 14d).
Let us now inspect how the reference trajectory generation worked during the execution of the task. The inputs for   Fig. 11a. We can observe that the user uses all the commands, i.e., roll, pitch and yaw rotations. Unlike the previous simulation, there is no obstacle avoidance input on the payload since in our architecture the obstacle avoidance does not act directly on the orientation of the payload. The rotational input passes through the DIB which blocks the commands when the condition τ (d rot , N ref )/τ (d rot , N des ) ≤ 0.88 is verified. We see in Fig. 15b that the ratio τ (d rot , N ref )/τ (d rot , N des ) (blue line) remains always above rot (red line) except for a short interval approximately between 3s and 5s. This moment corresponds to the initial yaw rotation, when the formation is too large and the obstacle avoidance stops some MAVs are from colliding with the columns Fig. 14b. Additionally, we observe that when the user gives no commands, i.e., the direction d rot is not defined, the ratio τ (d rot , N ref )/τ (d rot , N des ) is set to 1. As a result of the DIB activation the rotational trajectory is set to zero. In Fig. 15c we show directly the signalν ref produced by passing the output of the DIB through the input shaping method, which removes the discontinuity that would be otherwise present when the DIB is first activated. For what concerns the internal motions, we show the operator's command Pη hmn in Fig. 15d. We see that the user increases the elevation angle at around 5s, when the obstacle avoidance stops the motion of the MAVs that are too close to the columns (see Fig. 14b), thus resulting in the activation of the DIB. By squeezing the formation, the user manages to free the MAVs that were blocked by the columns, thus resulting in the deactivation of the DIB and the reprisal of the task. Afterwards, the user opens up the formation again at around 20s. We can now look at the plots of the azimuth and elevation angles of the cables, in Fig. 15e and 15f, to get a better understanding of how the shape of the formation is modified by the user command Pη hmn and by the obstacle avoidance terms ω obs,i acting on each cable/MAV. We see that the reference elevation angles follow accurately the desired values specified by the user. On the other hand, the reference azimuth angles for some of the cables deviate from their desired values at around 3s. This deviation is due to the obstacle avoidance that prevents the MAVs/cables to collide with the columns when the user is first commanding the yaw rotation (Fig. 14b). The mismatch of the azimuths vanishes once the user squeezes the formation.
In terms of task execution, the payload remains in a fixed position (see Fig. 16a) while tracking the the reference orientation (see Fig. 16b). We report also the tracking performance of one of the quadrotors (they are all comparable). Without loss of generality we select the first UAV and Fig. 16c shows that the reference trajectory is followed accurately. Finally, we report the attitude of the UAVs during the execution of the task in Fig. 16d to 16f. Now we inspect the force feedback produced by the haptic controller. This force feedback, shown in Fig. 17a, is informative of the mismatchν hmn −ν Since the obstacle avoidance does not modify directly the orientation of the payload, in this case we see that the feedback has two main causes: 1. The DIB, which generates an intense cue between 3s and 5s. Just like in the previous simulation, this cue pushes the user away from the obstacles and, consequently, away from the singular configuration. 2. The input shaping, which produces larger peaks at around 30s. These peaks are due to the fact that, unlike the translational case, for rotational commands the implemented input shaping imposes a saturation on the maximum roll and pitch angles (more details in Appendix A).
Finally, we can verify that the additional objectives concerning distance from the obstacles and cable forces have been satisfied throughout the execution of the task. Firstly, the minimum distance between the cables/robots and the obstacles is always above the minimum value 0.3m (see Fig. 17b). This value was arbitrarily chosen for the simulation, but it can be increased/reduced by modifying the artificial potential functions that generate the inputs of the obstacle avoidance. Secondly, the forces measured at the cables always satisfy the constraint 0.1N ≤ t i ≤ 10N (see Fig. 17c).

Robustness to external disturbances
The previous simulations illustrated the behaviour of the system with a human-in-the-loop. Now, to demonstrate the robustness of the system, we assess its behaviour in presence of external forces that are applied to the payload and MAVs to mimic the effect of wind or other disturbances. For this purpose, we consider a hovering task without any obstacles nor inputs from the human operator and focus on analyzing the deviations from the desired state induced by the external forces. Since in SimScape the cables are modelled as ideal connections (they have fixed length and are massless), the physical environment cannot simulate vibrations nor compression/elongation of the cables. Despite these limitations, the physical engine  can simulate lateral oscillations induced by wind, which is a major concern for such a system. To increase the difficulty of this test we consider external forces with different profiles (step and sinusoidal) and we apply these forces either to the payload alone or both to the payload and robots. This choice follows the intuition that a disturbance that is not applied uniformly to all the parts of the system (payload and MAVs) may induce more pronounced oscillations of the cables. The simulation is divided in four parts, in which different disturbances are considered: D1. a step force directed along the X in F W and with magnitude 5N is applied to the payload (see Fig. 18a); D2. a step force directed along the X in F W and with magnitude 5N is applied to the payload and a step force directed along the X in F W and with magnitude 2N is applied to every MAV (see Fig. 18b For what concerns the state of the MAVs we only plot the quantities from one robot, because drawing all the signals would make the graphs not readable. In particular, we show the results from the MAV that displays the largest deviations in its roll and pitch angles. Figure 19a and 19d show the forces applied on the payload and on the MAV, respectively. As mentioned earlier, the sinusoidal disturbances have different frequencies and phases, as it is noticeable in the part of the graphs pertaining to D4. The effect of the disturbances on the payload are visible in Fig. 19b and 19c. We can see that despite the strong forces the payload is stable, showing both in position and orientation bounded deviations from the desired configuration. Moreover, when the disturbance are removed the payload converges back to the desired state. From the plots we can appreciate that the payload remains stable in all sequences D1-D4. Curiously, comparing the configuration of the payload in D1 and D2 we observe that the position error is smaller in D1, whereas the orientation error is smaller in D2. This can be attributed to the fact that in D2 the MAVs are subjected to a force with the same orientation as the one applied to the payload, hence the whole system translates but the cables do not change orientation significantly. Indeed, this is also visible in Fig. 18, where in the case of D2 the cables are  (Fig. 18b), whereas in D1 the orientation of the cables changes more (Fig. 18a). This is confirmed by the plots of the azimuth and elevation angles of the cables during the simulation, in Fig. 20a and 20b. Finally, when sinusoidal forces are applied to system we can observe bounded oscillations in the position and orientation of the payload. As expected, the oscillations are more pronounced in D4, when we applied different sinusoidal forces to the payload and MAVs. The position and orientation of the selected MAV are plotted in Fig. 19e and 19f. We can observe that even when the disturbances are applied to the payload alone (D1 and D3) the MAV, being physically connected to it, is still affected. In particular, its attitude is modified quite significantly in order to compensate for the disturbances (see Fig. 19f), with the pitch angle that reaches a maximum tilt of -35 degrees from the nominal value of -15 degrees.

Conclusion
In this paper we have presented a novel solution for integrating a human operator in the control loop of a cooperative aerial transportation system with a cable-suspended payload. This is a problem that has not been addressed before and, given the complexity of the system, it needs a specialized solution that is manageable by a single person. Our framework is, to the best of our knowledge, the first shared control solution tailored for the cooperative aerial transportation system. The main feature of this paradigm is that it relieves the user from having to monitor the individual robots and allows him/her to focus on manipulating the payload. The obstacle avoidance is an important element to achieve this goal, because it allows the user to be less precise in his/her controls while guaranteeing the safety of the system. Similarly, we have provided a new mechanism for coping online with the possible singularities that are characteristic of this transportation system. The distinctive characteristic of this mechanism is the fact that it is embedded in the inputs blending of the shared control architecture and it requires no extra effort or additional training to the user. All the elements of the proposed architecture are also integrated by a bilateral interface that provides a force feedback that informs the user of how the system is executing the commands. Specifically, the force feedback encodes information about i) how well the payload is tracking the commands, ii)presence of obstacles, and ii)closeness to a singular configuration. In other words, through the force feedback the user becomes immediately aware of the limits of the system and of the automatic modifications injected by shared architecture, which otherwise may not be immediately transparent from the user perspective. Lastly, we have also formulated the architecture so as to exploit the internal motions of the system, which is one of the big advantages in using a system of robots connected by cables. We have demonstrated via simulations that the framework allows a single person to control such a transportation system even in complex environments and with difficult manipulation tasks.
In the future we plan to assess the framework with a user study and further validate it with experiments. Additionally, we believe that the proposed architecture can be extended in several direction. Firstly, the current controller is based on a model inversion, which is not robust to modeling uncertainties and perturbations. We plan to investigate a different strategy that does not require exact knowledge of the parameters of the system, and reformulate the controller in a distributed fashion with the feedback action based on the measurements/estimations of the cable forces and orientation. We also plan to move the shared control architecture towards higher levels of autonomy by fully automatizing the internal motions of the system so that they optimize some desired criteria. In particular, we want to investigate a predictive solution to control the shape of the formation so that it optimizes a combined cost that combines wrench capability, effort of the quadrotors and distance from the obstacles. Consider first the DIB filtering, restricting our discussion to the translational for the sake of simplicity. The tracking controller is defined as where k r > 0 and sat(•) is a function that imposes the limits where v(k) =p LP +k r (ṗ LP −ṗ ref ) p LP +k r (ṗ LP −ṗ ref ) is the direction of the unsaturated command Eq. 48.
In regards to Pη LP , the controller must additionally satisfy a limit on the maximum angle, i.e. P η min ≤ P η des ≤ P η max . For this purpose, a standard bang-bang strategy is used, i.e. decreased to m−n. Given the solution λ of Eq. 39, the optimal cable tension is simply t r = J T R N † w r + Nλ . In the simulations presented in Section 5 we solve the tension distribution problem using the reduced formulation Eqs. 51 and 53. In particular we impose Lagrange optimality conditions which, being the problem convex, ensure that a candidate point is a global minima.

Pη
Funding Open access funding provided by Istituto Italiano di Tecnologia within the CRUI-CARE Agreement.

Conflict of Interests
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.