# Agile and Cooperative Aerial Manipulation of a Cable-Suspended Load

Sihao Sun, Xuerui Wang, Dario Sanalitro, Antonio Franchi, Marco Tognon, Javier Alonso-Mora

► To cite this version:

Sihao Sun, Xuerui Wang, Dario Sanalitro, Antonio Franchi, Marco Tognon, et al. Agile and Cooperative Aerial Manipulation of a Cable-Suspended Load. Science Robotics, 2025, 10 (107), pp.1-46. 10.1126/scirobotics.adu8015. hal-05360829

HAL Id: hal-05360829

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

Submitted on 12 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Agile and Cooperative Aerial Manipulation of a
Cable-Suspended Load

Sihao Sun<fnref n="1" />*, Xuerui Wang<fnref n="1" />, Dario Sanalitro<fnref n="2" />,

Antonio Franchi<fnref n="3" />,<fnref n="4" />, Marco Tognon<fnref n="5" />, Javier Alonso-Mora<fnref n="1" />

Quadrotors can carry slung loads to hard-to-reach locations at high speed. Since a single quadrotor has limited payload capacities, using a team of quadrotors to collaboratively manipulate the full pose of a heavy object is a scalable and promising solution. However, existing control algorithms for multi-lifting systems only enable low-speed and low-acceleration operations due to the complex dynamic coupling between quadrotors and the load, limiting their use in time-critical missions such as search and rescue. In this work, we present a solution to substantially enhance the agility of cable-suspended multi-lifting systems. Unlike traditional cascaded solutions, we introduce a trajectory-based framework that solves the whole-body kinodynamic motion planning problem online, accounting for the dynamic coupling effects and constraints between the quadrotors and the load. The planned trajectory is provided to the quadrotors as a reference in a receding-horizon fashion and is tracked by an onboard controller that observes and compensates for the cable tension. Real-world experiments demonstrate that our framework can achieve at least eight times greater acceleration than state-of-the-art methods to follow agile trajectories. Our method can even perform complex maneuvers such as flying through narrow passages at high speed. Additionally, it exhibits high robustness against load uncertainties, wind disturbances, and does not require adding any sensors to the load, demonstrating strong practicality.

<footnote n="1">¹Delft University of Technology, the Netherlands.</footnote>

<footnote n="2">²Sorbonne Université, France.</footnote>

<footnote n="3">³University of Twente, the Netherlands.</footnote>

<footnote n="4">⁴Sapienza University of Rome, Italy</footnote>

<footnote n="5">⁵Univ Rennes, CNRS, Inria, IRISA, France</footnote>

<footnote>*Corresponding author. Email: sihao.sun@outlook.com</footnote>
# Video

A video of the experiment can be found at https://youtu.be/FBWN-rTK1YU

# Introduction

Quadrotors stand out for their unparalleled agility, speed, and mobility compared to other robotic systems. This unique capability has made them highly suitable for lifting and transporting objects to hard-to-reach locations at high speed (1, 2). However, the payload capacity of a single quadrotor is limited, prompting the exploration of utilizing multiple quadrotors in collaboration to transport (position control) and even manipulate (full pose control) heavy objects, resulting in a multi-lifting system (3–5). This strategy has great potential in a wide range of applications requiring heavy object manipulation, such as construction, disaster relief, and agriculture, as well as space exploration missions on Mars and Titan, where aerial vehicles have very limited resources and payload capacity (6, 7). Among the various manipulation mechanisms, the cable-suspended solution stands out for its simplicity and low weight (8–17). By connecting each quadrotor to a different location on the load through cables, a team of three quadrotors, or more, can change the full pose of the cable-suspended load by adjusting their positions, eliminating the need for additional mechanisms like robotic manipulators.

However, existing cooperative autonomous flight algorithms can only achieve pose control of a cable-suspended object at low speed and low acceleration, greatly limiting its performance and endurance in time-critical missions. The main challenge lies in addressing the complex dynamic coupling and kinematic constraints between the robots, cables, and the load. Early works typically resort to a quasi-static assumption to neglect the dynamic coupling effects (18–21), and only consider the kinematic constraints to determine the position and the path of quadrotors to reach the target pose of the load. Despite being simple, failing to account for dynamic coupling leads to undesired swinging motions, and cannot guarantee a safe load distribution on each quadrotor.

To account for dynamic coupling effects, recent works employ a force-based framework that employs the full dynamic model of the cable-suspended multi-lifting system. Given the pose of the reference load, the methods in this framework calculate a desired wrench (force and torque) that acts on the load through an outer-loop controller, for example, inverse-dynamics control (22), nonlinear model-predictive control (NMPC) (12, 23), and geometric control (14, 16, 24). Then the commanded wrench is allocated to each cable for their desired tension and directions through the Moore-Penrose inverse of the allocation matrix, which is determined by the connection points on the load (14). Some works further exploit the system redundancy in the null-space of the allocation matrix, offering capabilities for secondary tasks such as equal force distribution (15) and obstacle avoidance (12, 16) while retaining the collective wrench on the load. Once the required tension and
<figure><img src="image_2.png" /><figcaption>**Figure 1: Snapshot of the real-world experiments.** We propose an approach to control a cable-suspended load using multiple quadrotors with high agility. (A) Our approach enables agile full-pose control of a cable-suspended load. (B-D) It enables the quadrotors to dynamically control the load pose and fly through a narrow passage and a horizontally oriented gap. A summary of the experiments is highlighted in Movie 1.</figcaption></figure>
cable directions are determined, a mid-level controller calculates the command thrust and attitude for each quadrotor, which are then executed by an inner-loop attitude controller. However, despite considering the dynamic coupling effects, force-based methods are still far from fully exploiting the high agility of the cable-suspended system. In fact, only low-speed (under 1.5 m/s) and low acceleration (under 0.5 m/s²) flights have been successfully demonstrated in real-world experiments for load pose control using force-based approaches (12, 16, 24). Given that a loaded quadrotor with a thrust-to-weight ratio of 1.5 can easily reach an acceleration of over 4 m/s², existing solutions still largely compromise the inherent agility of quadrotors, making the cable-suspended multi-lifting system far from being able to operate in time-critical missions.

Here, we identify three major challenges obstructing the existing methods to achieve high agility in reality. First, the aforementioned force-based methods typically employ a cascaded control structure, which assumes that the load dynamics are substantially slower than those of the quadrotor. This assumption fails during agile flights, where the load needs to change its pose rapidly. With a cascaded control structure, the outer-loop commands can easily exceed the bandwidth of the inner loops, leading to instability, particularly in the presence of communication delays and actuator dynamics, which are often overlooked in simulation studies.

The second challenge is the high reliance on an accurate dynamic model, which is difficult to obtain. The mismatch of the model, especially the mass and inertia of the payload, leads to an error in the thrust command sent to each quadrotor, ultimately causing tracking error and even instability.

The third challenge is the reliance on high-frequency load and cable measurements for closed-loop control, requiring additional sensors to be installed onto the load, such as reflective markers for a motion capture system (15, 16, 24); or installing additional sensors on the quadrotors, such as downward-facing cameras (11), cable tension sensors, and cable direction sensors (25, 26). These methods inherently suffer from sensor noise and latency and typically require nontrivial engineering efforts for installation and calibration, making them largely impractical for day-to-day real-world operations.

## Trajectory-based framework

In this article, we propose a trajectory-based framework to address the above challenges. Our framework separates the controller into two submodules: an online kinodynamic motion planner and onboard trajectory tracking controllers. The kinodynamic motion planner considers the whole-body dynamics of the cable-suspended multi-lifting system, including the force-coupling effects, to generate dynamically feasible trajectories to each quadrotor in a receding-horizon fashion. Then, a trajectory tracking controller is deployed onboard each quadrotor to generate the rotor-speed-level commands to follow the online-generated trajectories while considering the effect of the cable forces.

Specifically, we formulate the kinodynamic motion planner into a finite-time optimal control
problem (OCP) which can be effectively solved within tens of milliseconds to generate predicted trajectories with a horizon of 2 s. The OCP formulates safety-related constraints as path constraints, including thrust limitations, cable tautness, inter-quadrotor collision avoidance, and obstacle avoidance. As the planner also takes into account the bandwidth and actuation constraints of the inner loop, the assumption of the time-scale separation principle required by existing solutions can be circumvented. The generated trajectories include the full state of quadrotors along the horizon, hence our method allows the planner to run at a considerably lower frequency (≤10Hz) than the outer-loop controllers of existing works (≥100Hz). This makes our method substantially more robust against the delay and noise on the state estimate of the load.

We deploy an estimator based on an extended Kalman filter (EKF) leveraging the load-cable dynamic model, quadrotor position and velocity estimates (generally available from an onboard state estimator), and accelerometers on quadrotors to provide satisfactory state estimates of the load and cables for the planner, achieving high-accuracy closed-loop tracking that outperforms state-of-the-art methods. The onboard trajectory tracking controller employs the incremental nonlinear dynamic inversion (INDI) technique (27-29) and leverages the differential-flatness property of quadrotors (30) to follow the reference trajectories and instantly compensates for the forces from cables using measurements from the inertial measurement unit (IMU). The mismatch in the planned cable tension that stems from the possible mismatch of the load inertia model is thereby effectively compensated for by the trajectory tracking controller, which eventually ensures high robustness against model uncertainties.

In the remainder of this article, we study the performance of the proposed trajectory-based framework in real-world experiments. The results reveal that the cable-suspended multi-lifting system controlled by our framework can achieve superior agility in pose control and trajectory following at high speeds (over 5 m/s) and accelerations (over 8 m/s²). It can even rapidly change configurations to avoid obstacles and fly through narrow passages dynamically (Fig. 1). Our method also shows robustness against load model uncertainties, external wind disturbances, and quadrotor state estimation errors. Moreover, the experiments were conducted without adding any sensors to the load to measure its pose, enhancing practicality in day-to-day real-world operations. The results and methods are summarized in Movie 1.

# Results

## Experimental Setup

We tested our algorithm through real-world experiments using three quadrotors to manipulate a 1.4 kg payload. Each quadrotor, weighing 0.6 kg, experienced a substantial additional force due to the payload. Without loss of generality, we set the cable length to 1 m for all quadrotors. These
<table><caption>**Table 1: Position tracking result.** Position root-mean-square-error (RMSE) in tracking references with different levels of agility. All reference trajectories were with a figure-eight shape. Our method substantially outperformed the two baseline methods (Geometric (14), and NMPC (12)), especially in tracking agile trajectories. The baseline methods were tested in a simulation environment, whereas our method was tested in both simulation and real-world experiments.</caption><thead><tr><th rowspan="2">Name of Ref.<br>trajectories</th><th rowspan="2">vel<sub>max</sub><br>[m/s]</th><th rowspan="2">acc<sub>max</sub><br>[m/s<sup>2</sup>]</th><th rowspan="2">jerk<sub>max</sub><br>[m/s<sup>3</sup>]</th><th>Geometric (14)</th><th>NMPC (12)</th><th>Ours</th><th>Ours (real-world)</th></tr><tr><th>[m]</th><th>[m]</th><th>[m]</th><th>[m]</th></tr></thead><tbody><tr><td>Slow</td><td>1</td><td>0.5</td><td>0.25</td><td>0.032</td><td>0.036</td><td><strong>0.031</strong></td><td>0.102</td></tr><tr><td>Medium</td><td>2</td><td>2</td><td>2</td><td>0.135</td><td>0.159</td><td><strong>0.067</strong></td><td>0.093</td></tr><tr><td>Medium Plus</td><td>2</td><td>4</td><td>8</td><td>Crash</td><td>Crash</td><td><strong>0.062</strong></td><td>0.117</td></tr><tr><td>Fast</td><td>5</td><td>8</td><td>16</td><td>Crash</td><td>Crash</td><td><strong>0.152</strong></td><td>0.197</td></tr></tbody></table>

cables were attached to three distinct points on the rigid-body payload to enable pose control, with the other end connected to each quadrotor 0.03 m below its center of gravity (CoG). The quadrotors were modified from the Agilicious open-source hardware platform (31), and each operates its onboard algorithms using a Raspberry Pi 5 mini PC. The centralized planner for our algorithm ran on a laptop at 10 Hz, sending commands to each quadrotor via WiFi.

We used motion capture systems to measure the poses of quadrotors at 100 Hz. These measurements were fused with onboard IMUs through an EKF to obtain state estimates of each quadrotor. On the other hand, the state of the load for closed-loop control was estimated from the quadrotor states, agnostic to the sensors and state estimation algorithms onboard each quadrotor. This also offers high practicality since no sensors are required to be attached to the load. The effect of quadrotor state estimation error, typically seen in field operations without a motion capture system, is analyzed in the Section **Robustness against quadrotor state estimation error**. A snapshot of the experimental setup is provided in Fig. S3.

## Agile Pose Control

To demonstrate that our method could control the cable-suspended multi-lifting system to achieve high agility, we tested its performance in tracking figure-eight trajectories with various levels of agility (increasing velocities, accelerations, and jerks), listed in Table 1. The algebraic expressions of the reference trajectories are given in Table S1. At the same time, the reference heading also varied over time with a constant yawing rate of 0.25 rad/s. We present the results of our method obtained from real-world experiments. We also present the results of our method compared with two baseline methods obtained in a simulation environment. It should be noted that the parameters of the quadrotors and loads in these experiments were kept consistent between different approaches to ensure a fair comparison.
<figure><img src="image_3.png" /><figcaption>**Figure 2: Performance in tracking the reference Fast.** (A) Experiment comparing our method against two baseline methods to follow the reference trajectory Fast, a figure-eight trajectory with a maximum speed of 5 m/s and a maximum acceleration of 8 m/s². The detailed expression of the reference is given in Table S1. (i) Top view of the flight path of the CoG of the load. (ii-iii) Time history of the root-mean-square error of the load position and attitude tracking error of the load. We used axis-angle representation for the attitude error. (B) Experiment comparing our methods with tightened thrust limits and without, while tracking the reference Fast. (i) Top view of the flight path of the CoG of the load. Once the maximum thrust was limited, the reference trajectory became dynamically infeasible for the system to follow accurately (red). (ii-iv) The commanded collective thrust of the three quadrotors with the reduced thrust limits (black dashed lines).</figcaption></figure>
We selected two representative state-of-the-art methods as the baseline: geometric control (14, 24) and NMPC (12) that have been successfully demonstrated in real-world experiments. These two force-based approaches both employ a conventional cascaded structure, namely using an outer-loop controller to generate the desired collective load wrench through a geometric controller or NMPC and distribute it to each quadrotor through an inner-loop controller. Table 1 lists the tracking error in these reference trajectories. Both baseline approaches could follow trajectories with relatively low agility (up until $v_{max}$ = 2 m/s, $a_{max}$ = 2 m/s²). However, they started to fail in following the reference Medium Plus, which involves higher peaks in acceleration and jerk, requiring rapid changes in the cables' directions to produce a fast time-varying wrench on the load. Our method, by contrast, avoids using the cascaded structure employed by the baseline methods and thus can consequently allow fast variation of load pose and cable directions. Therefore, it still successfully followed the reference Medium Plus and even the reference Fast, which has substantially larger accelerations and jerks. A video recording of the comparison in simulation environments is provided in Movie S1.

Fig. 2A presents the path and pose error while tracking the trajectory Fast that has a maximum velocity of 5 m/s and a maximum acceleration of 8 m/s². The reference velocity and acceleration started from zero and gradually reached their maximum values. As the reference velocity increased, both baseline methods failed to track the reference. By contrast, our method succeeded in tracking the reference trajectory with a position-tracking root-mean-square-error (RMSE) of 0.197 m, and an attitude-tracking RMSE of 12.9 deg, in real-world experiments. The high closed-loop tracking accuracy came from the combined efforts of our controller and estimator. The time history of pose reference, estimate, and ground truth is presented in Fig. S1.

Our algorithm considers dynamic coupling and thrust limits to prevent overloading the quadrotors. In another experiment, we limited the maximum thrust of two quadrotors on the same side of the load from 20 N to 11 N, and that of the third quadrotor on the opposite side to 15 N, since it needed to carry more lift. Then we let the system track the reference Fast. Consequently, the reference trajectory became dynamically infeasible for the thrust-limited multi-lifting system to follow precisely, leading to larger tracking error (0.363 m on position; 18.9 deg on attitude). The tracking result is presented in Fig. 2B. Despite these thrust limitations, our method still enabled the multi-lifting system to follow the reference trajectory and avoid instability. Interestingly, our controller modulated trajectory curvature around turns to lower the required acceleration. Throughout this process, our method ensured that the commanded thrust of each quadrotor was maintained within the reduced thrust limits. Additionally, the variation in the collective thrust of each quadrotor was notably reduced with a tightened thrust limit, demonstrating that our method automatically adjusted the level of agility to match the capabilities of the quadrotors.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Obstacle avoidance through dynamic motion. Both tasks were provided with a line segment reference that originally intersected the obstacles. (A) Task 1: Flight through a narrow passage between two walls. (i) Top view of the load center and three quadrotors with predicted trajectories at *t* = 1.5 s. (ii-iii) Velocity and acceleration profiles. (iv) Distances between quadrotors. (v) Snapshot of the experiment when the multi-lifting system flew through the narrow passage. (vi) Load inclination during traversal, defined as the angle between the load-fixed z-axis and the world-frame z-axis. (B) Task 2: Flight through a horizontally oriented narrow gap. (i) Side view of the trajectory and predicted trajectories at *t* = 1.5 s. (ii) Snapshot of the experiment when the multi-lifting system flew through the horizontally oriented gap. (iii-iv) Velocity and acceleration profiles. (v) Cable inclinations during traversal, defined as the angle between cable directions and the gravity.</figcaption></figure>
Obstacle Avoidance

Our algorithm enabled high-speed obstacle avoidance without designing obstacle-free trajectories
in advance, which required hundreds of seconds to generate for a multi-lifting system composed
of more than three quadrotors (16). Instead, we surrounded the obstacles with predefined no-fly
zones and formulated them as second-order inequality constraints in the OCP. This way, the online-
generated reference trajectories for the quadrotors ensured that the quadrotors and the load avoided
the no-fly zones, thereby preventing collisions with the obstacles.

We demonstrated the obstacle avoidance capability of our algorithm in two previously un-
explored challenging tasks. In both tasks, the multi-lifting system had to navigate through gaps
smaller than its original configuration size by leveraging its kinematic redundancy to reconfigure
and squeeze through the narrow passage. In the second task, traversal was performed dynamically,
exploiting the momentum gained at high speed, since the system could not counteract gravity with
its configuration at the moment of traversal if it attempted to fly statically. A video of the obstacle
avoidance experiments is provided in Movie S2.

Flying through a narrow passage

In the first scenario, the multi-lifting system was commanded to fly through a narrow passage with a width of 0.8 m. The width of the entire system in hovering condition was approximately 1.4 m, which was greater than the size of the gap. We first commanded the quadrotors to carry the load to hover at an initial position. Then we set a target at 6 m away from the initial position along the y-axis of the inertial frame. A minimum-snap reference trajectory (30) of the load was generated starting from the initial position to the target. However, this reference trajectory intersected with the obstacle. Without an obstacle avoidance mechanism, the system would have flown directly toward the wall and crashed. Note that we used a constant orientation reference, with the load frame aligned with the world frame.

To guide the system through the opening, we defined two vertical cylinders as no-fly zones, each with a radius of 1.5 m, encompassing the real obstacles. These two no-fly zones created a gap of 0.2 m for the system to pass through, ensuring a 0.3 m clearance from the real obstacles. We selected several reference points on the load and on each quadrotor. The algorithm then ensured that none of these reference points entered the no-fly zones. Specifically, the reference points in this experiment were the center of each quadrotor and the four edges of the payload. At the same time, the error between the actual and the reference pose was minimized in the cost function, encouraging the system to continue moving toward the final target pose.

Fig. 3A presents the experimental data, illustrating the maneuvering process. Our proposed algorithm generated predicted trajectories for both the load and the quadrotors at 10 Hz, allowing them to fly through the gap while adhering to the system kinodynamic model. The planner auto-
matically exploited the redundancy of the system to change the cable directions. Since the width of the load (0.54 m) was greater than the gap between the two no-fly zones (0.2 m), the quadrotors managed to steer the load at a steep inclination of approximately 70 deg to squeeze through the gap. The distances between quadrotors were also included as constraints in the optimization problem. Hence, their distances were kept greater than a safe margin (0.8 m) throughout the traversal.

Despite successfully avoiding the obstacles, the speed of the maneuver was not compromised. The system reached a top speed of 4 m/s during the fly-through maneuver, with a peak acceleration of over 5 m/s². The load passed through the gap within 1.2 s from the start of the maneuver and eventually stabilized at the target pose after successfully completing the traversal.

## Flying through a horizontally oriented narrow gap

In the second task, the quadrotors were commanded to carry the load through a horizontally oriented gap with a height of 0.6 m, while the height of the multi-lifting system in hover was around 1.2 m. The experimental data of this fly-through maneuver is presented in Fig. 3B. In this case, we defined two horizontal cylinders as no-fly zones, each with a radius of 4 m, ensuring that the obstacles were encompassed within the no-fly zones with a minimum safety margin of 0.2 m. The load was initially controlled to hover. Next, a target pose behind the gap was sent to the algorithm, which generated a minimum-snap reference trajectory. This trajectory, however, intersected with one of the obstacles.

Since the vertical size of the gap between the two no-fly zones was only 0.2 m, which was much smaller than the system's hovering height of approximately 1.2 m when all cables were nearly vertical, our algorithm controlled the quadrotors to spread out and stretch the cables, reducing the overall height of the system to enable it to pass through the gap. During this process, the cable directions changed rapidly from nearly vertical to almost horizontal within 1.2 s. At the moment of traversal, when the cables were nearly horizontal, the vertical components of the cable tensions could not compensate for the gravity of the load. Therefore, the algorithm induced dynamic motions and took advantage of the momentum of the system to complete the fly-through. In this maneuver, the load reached a maximum velocity of 4 m/s and a peak acceleration of over 7 m/s², generating the momentum necessary for a successful traversal.

## Robustness

In this section, we investigate our method's robustness against load model uncertainties, external wind disturbances, and quadrotor state estimation errors. These uncertainties are commonly seen in field operations.
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Test under load model uncertainties and communication delays. (A) Tracking performance of our method versus two baseline methods under various load model mismatches and communication delays while tracking reference Slow, Medium, and Fast as defined in Table S1. The baseline methods failed to follow the reference Fast even without mismatches, whereas our method remained robust. Each box corresponds to one run and summarizes the error at 4500 reference points: median (center line), 25th-75th percentiles (box), and whiskers extending to the minimum and maximum non-outlier values; outliers are defined as points lying beyond 1.5 interquartile range (IQR) from the box edges. (i–iii) Position error is in meters. (iv–vi) Attitude error is in degrees, calculated through axis-angle representations. (B) Real-world experiment where a 0.6 kg basketball was placed onto the 1.4 kg basket-shaped load and introduced a considerable inertia model mismatch of the load. Our method ran without knowing the presence of the basketball. (i) Top view of the path of the load CoG with and without a sloshing load. (ii–iii) Time history of the position and attitude tracking error. (iv) A snapshot of the experiment.</figcaption></figure>
## Robustness against load model uncertainties

Being robust to uncertainties in the load model is a desirable feature, as it is often impractical to obtain an accurate load model during real-world operations. In Fig. 4A, we compared our method with the two baseline controllers in the presence of various types of model uncertainties on the load in a well-controlled simulation environment that quantified the model mismatch. The results showed that the baseline controllers were more sensitive to the model mismatch, especially in attitude control. By contrast, our approach sustained over 50% mass and inertial mismatch in all the tests. These types of load model mismatch, commonly seen in practice, did not degrade the tracking performance, including the fastest reference trajectory given in Table 1. Only the case with a 10% bias of CoG led to an increase in attitude tracking error of about 5 deg. We also conducted simulations under step inputs of pose commands, which led to the same conclusion. The results are shown in Fig. S5.

We further conducted a challenging real-world experiment, where we placed a 0.6 kg basketball into the original basket-shaped load to introduce sloshing motion (Fig. 4B). This led to a mass mismatch of 43%, given that the mass of the original load was 1.4 kg. The motion of the basketball during flight also caused a pronounced time-varying CoG and inertia of the load if the basketball and the original load were considered as a single unit. We did not modify any parameters in the algorithm; the presence of the basketball was entirely unknown to our method. Despite that, we commanded the multi-lifting system to follow the trajectory Fast. Since the sloshing inevitably introduced additional swaying motions, particularly during dynamic maneuvers, the tracking error with the sloshing load was slightly larger (0.225 m vs. 0.197 m for position RMSE and 18.9 deg vs. 12.9 deg for attitude). Nevertheless, our algorithm managed to control the multi-lifting system to follow the reference Fast with an unknown sloshing load, which the baseline methods could not achieve even with a perfect model (see Table 1).

## Robustness against wind disturbance

We evaluated the performance of our method in both simulation and real-world experiments under windy conditions. In the simulation, we compared our method against the two baseline methods in the presence of various levels of wind. The system was commanded to hover at a target position. A horizontal wind was then ramped from zero to the designated speed over five seconds and persisted until the end of the simulation. We recorded the position error of the load once the system became stable. We employed the quadrotor drag model introduced in (29) with parameters identified in real-world experiments. For the load, we employed a second-order drag model with a reference surface area of 0.05 m² and a drag coefficient of 1.05 (value for a cube (32)). Fig. 5A presents the simulation results with wind up to 15 m/s (greater than the maximum wind resistance of most commercial quadrotors), showing that our method notably outperformed the baselines. Note that
the wind effect on the cables was neglected.

We also conducted real-world experiments under a wind field of around 5 m/s generated by a 1.5 m diameter fan. We first commanded the multi-lifting system (consisting of three and four quadrotors, respectively) to follow a straight line across the wind field at only 0.3 m/s, exposing the system directly to the wind for a duration of 5 seconds. We then evaluated higher-speed flight at 2 m/s by tracking a curved trajectory, which allowed us to assess performance under more dynamic conditions. In addition, we introduced the system into the wind field from an initially windless environment, creating a wind-gust-like scenario that further tested its ability to handle sudden changes in airflow.

Fig. 5D-F compares the top view of the trajectories under windy and windless conditions. Our framework enabled the multi-lifting system to operate at a moderate wind speed of 5 m/s. The disturbances acting on the quadrotors were effectively compensated for by the onboard flight controller. Since the aerodynamic model was not considered by the planner, the disturbance acting on the load led to greater tracking error compared to the case without wind disturbance (0.048 m vs. 0.055 m for position RMSE with three quadrotors; 0.048 m vs. 0.070 m with four quadrotors). Such tracking errors could be further reduced in future work by identifying and integrating a wind-effect model into the centralized planner. Videos of the above experiments are provided in Movie S4.

## Robustness against quadrotor state estimation error

Despite having no sensors on the load, our experiments still required state estimates from the quadrotors. Outside the lab environment without motion capture systems, greater estimation errors could occur due to imperfections in the quadrotor estimation algorithms and sensor noise from GPS, barometer, IMU, or onboard cameras. The quadrotor state estimation errors could deteriorate the closed-loop control performance of our framework. To test the sensitivity to quadrotor state estimation error, we deliberately added noise to the quadrotor position, attitude (Euler angle representation), and velocities. We selected random-walk noise on the quadrotor states, as it captured both the drift and stochastic fluctuations in the states. The level of random-walk noise was quantified by its standard deviation after 50 seconds starting from zero. We defined eight noise levels, each including 50 simulations where the system was commanded to follow the reference Fast.

Fig. 6A presents the Monte Carlo simulation results at each noise level, including the load position tracking RMSE and the success rate (i.e., no crash occurred). The closed-loop control performance gradually degraded as the noise level increased, which also reduced the success rate. Despite that, the success rate remained over 95% even under noise level 4, where the standard deviations of position, attitude, and velocities were respectively 0.07 m, 7 deg, and 0.07 m/s after 50 seconds of flight.

A real-world experiment was also conducted by adding random-walk pose and velocity noise to the original state estimator output of each quadrotor. Fig. 6B-E presents the estimation errors
<figure><img src="image_7.png" /><figcaption>**Figure 5: Test under wind disturbances.** (A) Comparison of position error between our approach and the baseline methods at different wind speeds in simulation environments. (B-C) Snapshot of experiments with three and four quadrotors, respectively, under windy conditions generated by a 1.5 m diameter fan. (D-E) Real-world experimental data from three and four quadrotors carrying a load to follow a straight line at a speed of 0.3 m/s in a 5 m/s wind field. (F) Real-world experimental data with four quadrotors carrying the payload and flying over the wind field while following a curved trajectory at a speed of 2 m/s. The videos of the experiments are provided in Movie S4.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Test under quadrotor state estimation errors. (A) Monte Carlo evaluation under varying levels of state estimation noise (50 runs per level), performed while tracking the Fast reference using noisy position, attitude, and velocity measurements. The red curve shows the success rate. The box plots summarize load position RMSE from the successful runs only: median (center line), 25th-75th percentiles (box), and whiskers extending to the minimum and maximum values. (B) Real-world flight test result under quadrotor state estimation error (level 3). We introduce random-walk noise on the original quadrotor state estimator that uses a motion capture system. (C-E) Time history of the quadrotor position, attitude, and velocity noise introduced in this flight test.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Computational load and scalability. (A) The CPU time to solve each OCP to follow the reference trajectory Fast with three quadrotors in real-world experiments. (B) Box plots of the CPU time of planning a trajectory with different numbers of quadrotors in a setpoint control task. Each box corresponds to data from 200 simulation steps: median (center line), 25th–75th percentiles (box), and whiskers extending to the minimum and maximum non-outlier values; outliers are defined as points lying beyond 1.5 IQR from the box edges. (C) Simulation result in a setpoint tracking task involving nine quadrotors. (i) 3D plot of the load and the CoG of the nine quadrotors. (ii-iii) Time history of the load pose (solid lines) compared with the reference pose (dashed lines). (iv) Time history of the distances between quadrotors (black solid lines) and the minimum distance allowed by our algorithm (red dashed line).</figcaption></figure>

of all the quadrotors and the tracking performance compared to the case without added errors. The position errors of the quadrotors exceeded 0.1 m over 90 seconds of flight, the attitude errors exceeded 10 deg, and the velocity errors exceeded 0.1 m/s, which was considerably less accurate than a commercial quadrotor operating in the field. Consequently, the position tracking error of the load increased from 0.197 m to 0.278 m, and the attitude tracking error of the load increased from 12.9 deg to 19.4 deg. Nevertheless, our method successfully controlled the multi-lifting system to follow the reference Fast without crashing. A video of the flights in simulation and in the real world under quadrotor state estimation error is provided in Movie S5.

## Computational Load and Scalability

Our method utilized a centralized structure, aiming to reach optimality in coordinating all the agents to manipulate the payload. The centralized structure, particularly the planner, posed a challenge to running in real time when multiple quadrotors were involved. The real-world experiments demonstrated that our method successfully coordinated three units in real time. Fig. 7A presents the
CPU time of the laptop running the planner (Intel Core i7-13700H) while tracking the trajectory Fast. The average CPU time consumed by the planner was 15.3 ms. Since the planner ran at 10 Hz (i.e., every 100 ms), our algorithm only took 15.3% of the computational budget.

We further explored the potential of scaling up our algorithm to include more units. First, we successfully scaled up our method to four units in real-world experiments, as shown in Movie S3. We also conducted simulations of setpoint tracking tasks with larger numbers of quadrotors. Without loss of generality, we assumed the load to be centrosymmetric, with cables connected to points evenly distributed around the load center at equal angular intervals depending on the number of quadrotors involved. The mass of the load was scaled proportionally with the number of quadrotors. We set the minimum distance constraint between quadrotors to be 1.6 times the distance between the corresponding contact points on the load, ensuring that constraint violations were possible if it was not imposed.

As a result, Fig. 7B shows that CPU time grew exponentially with the number of quadrotors included in the system. With our hardware setup, the planner supported nine units at 10 Hz without exceeding the compute budget. Fig. 7C presents the simulation results with nine quadrotors, including pose tracking, inter-quadrotor distances, and a 3D illustration. We believe that with a tailored optimizer for this particular problem, along with additional software and hardware optimizations, our method could be scaled up to substantially more units.

## Discussion

Our experiments have shown that the proposed trajectory-based framework can substantially enhance the agility, robustness, and practicality of cable-suspended multi-lifting systems compared to the state-of-the-art. To consider the dynamic coupling effects, the state-of-the-art framework requires a cascaded structure to coordinate and control multiple quadrotors to collaborate. This conventional structure is built upon the principle of time-scale separation, assuming that quadrotors can instantly generate a resultant wrench on the load requested by an outer-loop controller, which limits the outer-loop gains to prevent instability and makes the tuning process tedious and task-specific (16). Consequently, agility (high gain) and safety (low gain) are considered contradictory in the traditional framework.

In contrast, our method does not require the cascaded structure and addresses this issue by solving an online kinodynamic planning problem that considers the whole-body dynamics of the load-multi-quadrotor system. The solution generates the states of the system over a future horizon, offering a predictive capability that allows for the inclusion of safety constraints during agile motion, rather than simply limiting gains in the traditional framework. This key component of our solution enables precise pose control, trajectory tracking, and obstacle avoidance at high speeds (over 5 m/s) and accelerations (over 8 m/s²). Although the kinodynamic motion planning accounts
for the whole-body dynamics, we have demonstrated that this problem is solvable on a midrange CPU in just a few milliseconds, enabling fast online generation to adjust to disturbances and avoid obstacles, and even with the potential to scale up to nine units with the current hardware and software setup.

We send the predicted trajectories to each quadrotor instead of a single reference point, offering two major advantages. First, deploying a reference sampler together with a robust trajectory tracking controller on the quadrotors makes our framework substantially more robust to load model uncertainties, communication delays, and external wind disturbances. This enables our algorithm to safely control the cable-suspended multi-lifting system, even when handling an unknown-mass sloshing load that brings over 40% mismatch on the load mass model, or under moderate wind breeze at 5 m/s. Second, sending a trajectory allows us to run the planner at over ten times lower frequency than traditional controllers, avoiding reliance on high-frequency measurements from sensors installed on the load.

We see several opportunities for future work. Beyond the quadrotor multi-lifting problem, the trajectory-based framework has the potential to be applied to a wider range of robotic collaboration challenges, particularly those involving dynamic coupling, agility, and safety constraints. Although our method guarantees high accuracy in the presence of load mismatch, underestimating the load mass and inertia can lead to violations of the maximum thrust constraints. This problem can be alleviated by estimating the inertial properties online (33), or using the constrained tightening technique employed by a robust nonlinear optimal control framework, such as robust MPC (34, 35). This requires a pre-estimation of the uncertainties and provides a more conservative, yet safer reference for quadrotors to avoid any violation of their dynamical constraints. Another opportunity is combining our method with onboard perception algorithms in a GPS-denied environment. This requires aligning coordinate frames across quadrotors, which is still a challenge for multi-agent perception algorithms, especially in agile flights. To support future works in this direction, we provide a preliminary analysis in the Supplementary Discussion to demonstrate the performance of our algorithm in the presence of misalignment between quadrotor coordinate frames.

Overall, our work paves the way for future aerial manipulation systems with substantially higher resilience, versatility, and agility to perform complex collaborative tasks in day-to-day operations, from search and rescue to precision delivery in difficult terrains.

# Materials and Methods

An overview of the method is shown in Fig. 8. The proposed framework incorporates an optimization-based kinodynamic motion planner that generates real-time reference trajectories for the quadrotors. It also includes a time-based sampler and an INDI-based trajectory tracking controller onboard each quadrotor. In addition, the framework employs a centralized EKF to estimate the load pose and
cable directions from the quadrotors' position, velocity, and IMU measurements. All modules are model-based and rely on the dynamic model of the cable-suspended multi-lifting system. The following sections provide a detailed description of each module.

# Modeling of cable-suspended multi-lifting systems

## Load-cable dynamic model

The load-cable dynamic model describes the 6-DoF motion of the load and the motions of all the cables attached to the load. Specifically, we employed the following definition of the state of the load-cable dynamic model:

$$ \mathbf{x} = [\mathbf{p}, \mathbf{v}, \mathbf{q}, \omega, s_1, \mathbf{r}_1, \dot{\mathbf{r}}_1, \ddot{\mathbf{r}}_1, t_1, \dot{t}_1, \dots, s_n, \mathbf{r}_n, \dot{\mathbf{r}}_n, \ddot{\mathbf{r}}_n, t_n, \dot{t}_n]^\top, \quad (1) $$

where $n$ is the number of quadrotors, $\mathbf{p} \in \mathbb{R}^3$, $\mathbf{v} \in \mathbb{R}^3$ are positions and velocities of the load, $\mathbf{q} \in \mathbb{S}^3$ is the unit quaternion describing the load attitude, $\omega \in \mathbb{R}^3$ is the load angular velocity expressed in the load-fixed coordinate frame $\mathcal{F}_L$. Here the subscript $i$ indicates variables of the cable connected to the $i$-th quadrotor, where $s_i \in \mathbb{S}^2$ is the cable direction pointing from the quadrotor to the load, $\mathbf{r}_i \in \mathbb{R}^3$ is the cable angular velocity, $t_i \in \mathbb{R}_{\ge 0}$ is the cable tension. An illustration of the reference frames and some symbols defined above can be found in Fig. S4.

We also adopted the following dynamic equations for the load:

$$ \begin{aligned} \dot{\mathbf{p}} &= \mathbf{v}, & \dot{\mathbf{v}} &= -\sum_{i=1}^{n} t_i s_i / m + \mathbf{g}, \\ \dot{\mathbf{q}} &= \frac{1}{2} \mathbf{\Lambda}(\mathbf{q}) \begin{bmatrix} 0 \\ \omega \end{bmatrix}, \end{aligned} \quad (2) $$

$$ \mathbf{J}\dot{\mathbf{\omega}} = -\mathbf{\omega} \times \mathbf{J}\mathbf{\omega} + \sum_{i=1}^{n} t_i (\mathbf{R}(\mathbf{q})^{\top} \mathbf{s}_i \times \mathbf{\rho}_i), $$

where $\mathbf{J} \in \mathbb{R}^{3 \times 3}$ is the load inertia, $m$ is the load mass, $\mathbf{\rho}_i \in \mathbb{R}^3$ is the displacement of the $i$-th attachment point, expressed in the load frame, and $\mathbf{g}$ is the constant gravity vector. $\mathbf{\Lambda}(\mathbf{q})$ represents the quaternion multiplication, $\mathbf{R}(\mathbf{q}) \in \text{SO}(3)$ is the rotation matrix of the unit quaternion $\mathbf{q}$. To ensure smooth quadrotor reference trajectories up to the jerk level, we employed the cable kinematic model, with the third-order derivative of cable angular velocity and the second-order derivative of cable thrust treated as bounded inputs, yielding

$$ \dot{\mathbf{s}}_i = \mathbf{r}_i \times \mathbf{s}_i, \quad \ddot{\mathbf{r}}_i = \gamma_i, \quad \ddot{\lambda}_i = \lambda_i, \quad \text{for } i = \{1, \dots, n\}, \quad (3) $$

where $\gamma_i \in \mathbb{R}^3$ is the snap of cable directions, and $\lambda_i \in \mathbb{R}$ is the second-order derivative of cable tensions. In the Supplementary Methods, we proved that the generated quadrotor trajectories are smooth up to the jerk level and also lead to a smooth angular velocity reference with states defined in the load-cable dynamic model in Equation 1, as long as the $\gamma_i$ and $\lambda_i$ are bounded.
<figure><img src="image_10.png" /><figcaption>**Figure 8: Method overview.** Our framework includes a kinodynamic motion planner solving an OCP online at 10 Hz to generate receding-horizon reference trajectories of quadrotors given load reference pose and predefined no-fly zones. The OCP utilizes the whole-body dynamics of the system, including the quadrotor model and the load-cable model. The load's pose, twist, and cable directions are obtained from an EKF-based estimator. The remaining elements in the initial state of the OCP, namely the derivatives of the cable directions and tensions, are obtained by resampling the previously generated predicted trajectory to avoid oscillatory motion of the quadrotor when a new reference arrives. The load state estimator fuses the load-cable model, and the quadrotors' position, velocity, and IMU measurements to obtain estimates of load pose, twist, and cable directions. It is initialized through an iterative Kabsch-Umeyama algorithm given the initial quadrotor states. Onboard each quadrotor, a time-based sampler samples the received receding-horizon reference trajectory using the current timestamp to generate a single reference point, which is tracked by a trajectory tracking controller based on the INDI technique that regards the cable tensions as external disturbances, and compensates for them using the IMU measurements.</figcaption></figure>
## Quadrotor dynamic model

For the i-th quadrotor, we described the state space as $x_i = [p_i, v_i, q_i, \omega_i]^T$, which corresponds to the CoG position $p_i \in \mathbb{R}^3$, velocity $v_i \in \mathbb{R}^3$, unit quaternion rotation $q_i \in \mathbb{S}^3$, and angular velocity expressed in the quadrotor-fixed coordinate frame $\omega_i \in \mathbb{R}^3$. We employed the following rigid-body dynamics to derive the quadrotor equations of motion (36):

$$ 
\begin{aligned}
\dot{\mathbf{p}}_i &= \mathbf{v}_i, \quad \dot{\mathbf{v}}_i = (T_i z_i + t_i s_i + \mathbf{f}_{a,i})/m_i + \mathbf{g}, \\
\dot{\mathbf{q}}_i &= \frac{1}{2} \mathbf{\Lambda}(\mathbf{q}_i) \begin{bmatrix} 0 \\ \omega_i \end{bmatrix},
\end{aligned}
\quad (4) $$

$$ \mathbf{J}_i \dot{\omega}_i = -\omega_i \times \mathbf{J}_i \omega + \boldsymbol{\tau}_i + \boldsymbol{\tau}_{a,i}, $$

where $m_i$ and $\mathbf{J}_i \in \mathbb{R}^{3 \times 3}$ are respectively the mass and the inertia matrix of the i-th quadrotor; $T_i \in \mathbb{R}_{\ge 0}$ is the collective thrust; $z_i \in \mathbb{S}^2$ is the thrust direction aligning with the z-axis of the quadrotor body-fixed frame $\mathcal{F}_i$; $\mathbf{f}_{a,i} \in \mathbb{R}^3$ and $\boldsymbol{\tau}_{a,i} \in \mathbb{R}^3$ are the aerodynamic drag force and torque; $\boldsymbol{\tau}_i \in \mathbb{R}^3$ is the control torque generated by the rotors.

## Kinematic constraints

We assumed that the cables' tautness would be maintained by our algorithm throughout the operation, even during agile motions. On the i-th quadrotor, the position of the cable contact point in the inertial frame is denoted as $\mathbf{p}_i$. Then, the following kinematic constraint between the i-th quadrotor and load-cable dynamics holds:

$$ \mathbf{p}_i = \mathbf{p} + \mathbf{R}(\mathbf{q})\rho_i - l_i s_i, \quad (5) $$

where $l_i$ is the length of the cable.

## Online kinodynamic motion planner

### Finite-time optimal control problem

Our framework includes a centralized kinodynamic motion planner that generates smooth reference trajectories of all quadrotors in a receding-horizon fashion while considering the dynamic coupling between load and quadrotors. Specifically, the planner is a discretized finite-time OCP, solved by a
multi-shooting method (37)

$$
\begin{equation}
\begin{aligned}
\min J = & \sum_{k=0}^{N-1} \left( \|\mathbf{x}_k - \mathbf{x}_{k,\text{ref}}\|_Q^2 + \|\mathbf{u}_k - \mathbf{u}_{k,\text{ref}}\|_R^2 \right) \\
& + \|\mathbf{x}_N - \mathbf{x}_{N,\text{ref}}\|_P^2
\end{aligned}
\tag{6}
\end{equation}
$$

$$
\begin{array}{ll}
\text{subject to} & \mathbf{x}_0 = \mathbf{x}_{\text{init}}, \quad \mathbf{x}_{k+1} = f(\mathbf{x}_k, \mathbf{u}_k), \\
& \mathbf{h}(\mathbf{x}_{k+1}, \mathbf{u}_k) \le \mathbf{0}, \quad k \in \{0, \dots, N\},
\end{array}
$$

where the state equation uses the load-cable dynamic model (Equation 2 and 3), and the input
is $\mathbf{u} = [\gamma_1, \lambda_1, \dots, \gamma_n, \lambda_n]$. The quadrotor dynamics, albeit not explicitly included in the state
equation for the reason of numerical efficiency, are included in the path constraints, particularly to
avoid overloading quadrotors, and to perform obstacle avoidance.

The cost function is in a standard quadratic form to minimize the load pose tracking error and control effort for smoothness. The reference states used in the cost function of the OCP are precomputed based on a polynomial load pose reference, with the remaining states derived using the flatness property of the cable-suspended multi-lifting system (10). It is worth noting that in our experiments, the load reference did not consider avoiding obstacles. Instead, we left the planner to decide and generate commands for the quadrotors to carry the load to avoid the obstacles. In other words, the planner had the flexibility to deviate from the reference position or adjust the configurations to satisfy the obstacle avoidance constraints. Naturally, our method could also follow load references generated by a higher-level offline planner (e.g., (17)), other than the simple polynomial reference.

We discretized the horizon into $N = 20$ non-equidistant segments, with intervals linearly increasing along the horizon. Hence, it ensured higher fidelity of the predicted trajectory in the near future while extending the horizon length without increasing the number of discretization nodes. This OCP was subsequently solved through the sequential quadratic programming (SQP) algorithm in a real-time iteration (RTI) scheme (38), implemented using ACADOS toolkit (39). The solution of the OCP was the optimal input $\mathbf{u}_k^*$ and load-cable state $\mathbf{x}_k^*$ along the horizon

$$
\begin{align}
U^* &= [\mathbf{u}_0^*, \mathbf{u}_2^*, \dots, \mathbf{u}_{N-1}^*], \notag \\
X^* &= [\mathbf{x}_1^*, \mathbf{x}_2^*, \dots, \mathbf{x}_N^*] = \pi(U^*, \mathbf{x}_{\text{init}}). \tag{7}
\end{align}
$$

Once the optimal state sequence $\mathbf{X}^*$ was obtained, we converted it to the position, velocity, acceleration, and jerk of the quadrotor through kinematic constraints (Equation 5) and its derivatives. It is worth noting that the headings of quadrotors, defined as the rotation angle around $z_i$, did not affect the thrust directions or the motion of the load. Therefore, we avoided explicitly setting the heading reference for the quadrotors but let the quadrotors maintain a zero yaw rate instead.
The OCP was then solved at a fixed frequency to generate a new trajectory online. The initial state $\mathbf{x}_{\text{init}}$ in OCP was provided partially by the load-cable state estimator described in the Section Load-cable state estimator, and partially by resampling the trajectory from the latest OCP solution. Specifically, we used the estimated load pose, twist, and cable directions to renew $\mathbf{x}_{\text{init}}$ to make sure the trajectories employ the up-to-date state of the load for closed-loop control. The other states (cable rate, cable tensions, and their higher-order derivatives) were directly estimated by resampling on the previously generated trajectory. We observed that this treatment has two benefits. First, it ensured smooth transitions between consecutive reference trajectories, avoiding any abrupt and jerky maneuvers by the quadrotors. It also avoided numerical differentiation of state estimator values, which is usually impractical due to the high requirements for accuracy and smoothness in the estimator output.

Path constraints

We included several path constraints $\mathbf{h}(\mathbf{x}) \le \mathbf{0}$ in the OCP to ensure safety. To avoid overloading each quadrotor, we included the thrust constraint:

$$T_{i,\min} \le T_i(\mathbf{x}) \le T_{i,\max}, \qquad (8)$$

where $T_i(\mathbf{x})$ is the thrust of each quadrotor as a function of the state of the load-cable dynamic model. Specifically, $T_i(\mathbf{x})$ was obtained through quadrotor dynamics (Equation 4),

$$T_i(\mathbf{x}) = \| (\dot{\mathbf{v}}_i(\mathbf{x}) - \mathbf{g}) m_i - t_i s_i - \mathbf{f}_{a,i} \|, \qquad (9)$$

where $\dot{\mathbf{v}}_i(\mathbf{x})$ was calculated from the second-order derivative of the kinematic constraints (Equation 5).

We assumed that the cables would remain taut throughout the operation. Therefore, a cable tension constraint is included.

$$0 < t_{\min} \le t_i \le t_{\max}, \qquad (10)$$

To avoid inter-quadrotor collisions, the minimum distance constraints were also provided for every pair of quadrotors indexed by i and j

$$0 < d_{\min} \le \| \mathbf{p}_i(\mathbf{x}) - \mathbf{p}_j(\mathbf{x}) \|, \qquad (11)$$

where $d_{\min}$ is the predefined minimum distance, $\mathbf{p}_i(\mathbf{x})$ and $\mathbf{p}_j(\mathbf{x})$ are the positions of the i-th and the j-th quadrotor.

We also established several control points on the system to ensure it avoided obstacles. Without loss of generality, we used the CoG of each quadrotor together with the attaching points on the load. For each obstacle, and each control point denoted by $\mathbf{p}_c(\mathbf{x})$, the following constraint ensured that none of the control points entered the no-fly zone encompassing the obstacle

$$d_{o,\min}^2 \le (\mathbf{p}_c(\mathbf{x}) - \mathbf{p}_o)^{\top} \mathbf{C} (\mathbf{p}_c(\mathbf{x}) - \mathbf{p}_o), \qquad (12)$$
where $C \in \mathbb{R}^{3 \times 3}$ is a diagonal matrix controlling the shape of the no-fly zone, $\boldsymbol{p}_o$ is the center of the no-fly zone, $d_{o,min}$ is the safe distance from the control points to the center. The position and shape of the obstacle can be determined either offline or detected online. In the case when the obstacle is detected online, we can set up the problem to support a large number of obstacles and set the inactive obstacle constraints with a zero radius or place them far from the current positions. Once a new obstacle is detected, we activate one of the reserved inequality constraints by adjusting its parameters to match the size and location of the detected obstacle.

Lastly, to ensure the bounded input to the OCP, the control input constraint $\boldsymbol{u}_{\min} \le \boldsymbol{u} \le \boldsymbol{u}_{\max}$ was imposed in the experiments.

It is worth noting that all path constraints were inequality constraints and handled using slack variables to ensure problem feasibility.

## Load-cable state estimator

To update $\boldsymbol{x}_{\text{init}}$ in the OCP of the planner, the load pose, twist, and cable directions must be estimated in real-time. Instead of relying on additional sensors such as downward-facing cameras (11) or adding motion capture markers on the load (24) in the state-of-the-art approaches, our proposed estimator utilizes only the quadrotors' states and the load-cable dynamics. This eliminates the need for any hardware modifications.

We chose EKF to solve this state estimation problem thanks to its simplicity and computational efficiency. We omit the detailed steps in EKF and only describe the selections of states, measurements, models, and initialization.

The state vector of EKF comprised pose and twist of the load, as well as the positions and velocities of all quadrotors, namely $\hat{\boldsymbol{x}} = [\boldsymbol{p}, \boldsymbol{v}, \boldsymbol{q}, \omega, \boldsymbol{p}_1, \boldsymbol{v}_1, ..., \boldsymbol{p}_n, \boldsymbol{v}_n]^T$. State prediction was performed using the load dynamics and quadrotor dynamics (Equations 2 and 4). The cable directions to solve the load dynamic were obtained through the kinematic constraint (Equation 5). The cable forces in these equations were estimated through a spring-damper model, i.e., $t_i = k_{\text{stiff}}d_i + k_{\text{damp}}\dot{d}_i$ where $d_i$ is the distance between the position of the i-th quadrotor and its connection point on the load, namely $d_i = \|\boldsymbol{p}_i - \boldsymbol{R}(\boldsymbol{q})\boldsymbol{\rho}_i - \boldsymbol{p}\|$; $k_{\text{stiff}}$ and $k_{\text{damp}}$ are positive coefficients.

The EKF took the cables' directions, together with the quadrotors' positions and velocities as measurements, namely $\tilde{\boldsymbol{y}} = [\tilde{\boldsymbol{s}}_1, \tilde{\boldsymbol{p}}_1, \tilde{\boldsymbol{v}}_1, ..., \tilde{\boldsymbol{s}}_n, \tilde{\boldsymbol{p}}_n, \tilde{\boldsymbol{v}}_n]^T$. The quadrotor positions and velocities, and their covariances, were obtained directly from their onboard state estimators. The cable directions were obtained indirectly from the accelerometer sensor that is commonly available on a quadrotor drone. Since the accelerometer directly measures the specific force (the mass-normalized force excluding gravity), it captures the combined forces, including the cable tension, aerodynamic drag, wind force, and rotor thrusts. For each quadrotor, we identified a collective thrust model $\tilde{T}_i$, and a
drag model $\tilde{f}_{a,i}$ with the following expressions:

$$ \bar{T}_i = \sum_{j=1}^{4} c_j \omega_{j,i}^2, \quad \bar{f}_{a,i} = \mathbf{R}(\mathbf{q}_i) \mathbf{D}_a \mathbf{R}(\mathbf{q}_i)^{\top} \mathbf{v}_i, \qquad (13) $$

where $c_j$ is the thrust coefficient of the rotors, $\omega_{j,i}$ is the rotor speed, $\mathbf{D}_a \in \mathbb{R}^{3 \times 3}$ is the aerodynamic coefficient matrix (40). According to the quadrotor dynamics (Equation 4), subtracting them from the accelerometer readings provides the force vector from cables. Then the cable directions were approximated by

$$ \tilde{s}_i = (m_i \mathbf{a}_i - \bar{T}_i \mathbf{z}_i - \tilde{\mathbf{f}}_{a,i}) / \|m_i \mathbf{a}_i - \bar{T}_i \mathbf{z}_i - \tilde{\mathbf{f}}_{a,i}\|, \qquad (14) $$

where $\mathbf{a}_i = \dot{\mathbf{v}}_i - \mathbf{g}$ is the unbiased accelerometer measurement.

We observed that the covariance matrix of the cable direction was challenging to determine directly from the accelerometer properties, as it depends on the accuracy of $\bar{T}_i$, $\tilde{\mathbf{f}}_{a,i}$, and is also affected by wind disturbance. Therefore, we tuned this covariance experimentally. We observed that underestimating the covariance could result in a noisy load pose estimate, whereas overestimating it led to an excessively free motion in the load estimate.

We initialized the EKF with a first-order guess of load pose and twist, and cable directions. We assumed a static initial load state, namely a zero twist. As for the load pose and cable directions, we propose an algorithm to provide a guess iteratively through the Kabsch-Umeyama algorithm (41). The details of the algorithm are provided in Algorithm S1.

## Trajectory tracking controller on quadrotors

In our setup, every 100 ms, the most recently generated reference trajectories by the planner were sent to each quadrotor, and then followed by a differential-flatness-based trajectory tracking controller deployed onboard, modified from (29). Since the trajectory tracking controller operated at a higher frequency (300 Hz) than the interval between nodes in the reference trajectory, a time-based sampler was implemented to generate high-frequency reference states by linearly interpolating between the discretized nodes of the reference trajectory. The sampler continued to sample along the reference trajectory until a new reference was received.

The onboard trajectory tracking controller then computed the thrust command, including magnitude $T_{i,des}$ and direction $z_{i,des}$, through the following PD controller:

$$ T_{i,des} z_{i,des} / m_i = \mathbf{K}_p (\mathbf{p}_{i,ref} - \mathbf{p}_i) + \mathbf{K}_v (\mathbf{v}_{i,ref} - \mathbf{v}_i) + \dot{\mathbf{v}}_{i,ref} + \mathbf{f}_{\text{ext}} / m_i, \qquad (15) $$

where $\mathbf{K}_p \in \mathbb{R}^{3 \times 3}$ and $\mathbf{K}_v \in \mathbb{R}^{3 \times 3}$ are positive definite gain matrices; $\mathbf{f}_{\text{ext}}$ are external forces represents external forces on the quadrotor, excluding thrust and gravity, namely cable tension, aerodynamic drag, and wind. We estimated external forces using the accelerometer on the quadrotor
through the relationship $f_{\text{ext}} = m_i a_{i,\text{filtered}} - f_{i,\text{filtered}}$, where $a_{i,\text{filtered}}$ is the unbiased and low-pass-filtered accelerometer measurement, $f_{i,\text{filtered}}$ is the current collective thrust denoised with the same filter. Then, we used a tilt-prioritized attitude controller (42) to generate the angular acceleration command $\alpha_{i,\text{des}}$ from the desired attitude $z_{i,\text{des}}$, the reference jerk, and the zero yaw rate reference.

The angular acceleration and force commands were subsequently allocated to rotor speed commands through an INDI inner-loop controller, which is a sensor-based adaptive controller robust against external torque disturbances such as aerodynamic torque, motor differences, quadrotor CoG bias, etc. Finally, the rotor speed commands generated by INDI were sent to the electronic speed controllers (ESCs) through the DShot protocol. We refer interested readers to (31) for further details about the hardware implementations. We also provide key equations for INDI in the Supplementary Methods.
References and Notes

1. D. Hanover, P. Foehn, S. Sun, E. Kaufmann, D. Scaramuzza, Performance, precision, and payloads: Adaptive nonlinear mpc for quadrotors. *IEEE Robotics and Automation Letters* **7** (2), 690–697 (2021).

2. A. Saviolo, G. Loianno, Learning quadrotor dynamics for precise, safe, and agile flight control. *Annual Reviews in Control* **55**, 45–60 (2023).

3. G. Loianno, V. Kumar, Cooperative transportation using small quadrotors using monocular vision and inertial sensing. *IEEE Robotics and Automation Letters* **3** (2), 680–687 (2017).

4. A. Tagliabue, M. Kamel, S. Verling, R. Siegwart, J. Nieto, Collaborative transportation using mavs via passive force control, in *IEEE International Conference on Robotics and Automation* (2017), pp. 5766–5773.

5. A. Tagliabue, M. Kamel, R. Siegwart, J. Nieto, Robust collaborative object transportation using multiple MAVs. *The International Journal of Robotics Research* **38** (9), 1020–1044 (2019).

6. J. Balaram, M. Aung, M. P. Golombek, The ingenuity helicopter on the perseverance rover. *Space Science Reviews* **217** (4), 56 (2021).

7. R. D. Lorenz, et al., Dragonfly: A rotorcraft lander concept for scientific exploration at Titan. *Johns Hopkins APL Technical Digest* **34** (3), 14 (2018).

8. P. Foehn, D. Falanga, N. Kuppuswamy, R. Tedrake, D. Scaramuzza, Fast trajectory optimization for agile quadrotor maneuvers with a cable-suspended payload, in *Robotics: Science and Systems* (2017).

9. F. Panetsos, G. C. Karras, K. J. Kyriakopoulos, GP-based NMPC for aerial transportation of suspended loads. *IEEE Robotics and Automation Letters* **10** (1), 524–531 (2025).

10. K. Sreenath, V. Kumar, Dynamics, control and planning for cooperative manipulation of payloads suspended by cables from multiple quadrotor robots, in *Robotics: Science and Systems* (2013).

11. G. Li, R. Ge, G. Loianno, Cooperative transportation of cable suspended payloads with mavs using monocular vision and inertial sensing. *IEEE Robotics and Automation Letters* **6** (3), 5316–5323 (2021).

12. G. Li, G. Loianno, Nonlinear model predictive control for cooperative transportation and manipulation of cable suspended payloads with multiple quadrotors, in *IEEE/RSJ International Conference on Intelligent Robots and Systems* (2023), pp. 5034–5041.
13. T. Lee, K. Sreenath, V. Kumar, Geometric control of cooperating multiple quadrotor UAVs with a suspended payload, in *IEEE Conference on Decision and Control* (2013), pp. 5510–5515.

14. T. Lee, Geometric control of quadrotor UAVs transporting a cable-suspended rigid body. *IEEE Transactions on Control Systems Technology* **26** (1), 255–264 (2017).

15. J. Geng, P. Singla, J. W. Langelaan, Load-distribution-based trajectory planning and control for a multilift system. *Journal of Aerospace Information Systems* **19** (5), 366–381 (2022).

16. K. Wahba, W. Hönig, Efficient optimization-based cable force allocation for geometric control of a multirotor team transporting a payload. *IEEE Robotics and Automation Letters* **9** (4), 3688–3695 (2024).

17. K. Wahba, J. Ortiz-Haro, M. Toussaint, W. Hönig, Kinodynamic motion planning for a team of multirotors transporting a cable-suspended payload in cluttered environments. *arXiv preprint arXiv:2310.03394* (2023).

18. N. Michael, J. Fink, V. Kumar, Cooperative manipulation and transportation with aerial robots. *Autonomous Robots* **30**, 73–86 (2011).

19. M. Manubens, D. Devaurs, L. Ros, J. Cortés, Motion planning for 6-D manipulation with aerial towed-cable systems, in *Robotics: Science and Systems* (2013).

20. E. Fresk, G. Nikolakopoulos, Full quaternion based attitude control for a quadrotor, in *European Control Conference* (2013), pp. 3864–3869.

21. D. Sanalitro, H. J. Savino, M. Tognon, J. Cortés, A. Franchi, Full-pose manipulation control of a cable-suspended load with multiple UAVs under uncertainties. *IEEE Robotics and Automation Letters* **5** (2), 2185–2191 (2020).

22. C. Masone, P. Stegagno, Shared control of an aerial cooperative transportation system with a cable-suspended payload. *Journal of Intelligent & Robotic Systems* **103** (3), 40 (2021).

23. S. Sun, A. Franchi, Nonlinear MPC for full-pose manipulation of a cable-suspended load using multiple UAVs, in *2023 International Conference on Unmanned Aircraft Systems (ICUAS)* (IEEE) (2023), pp. 969–975.

24. G. Li, X. Liu, G. Loianno, RotorTM: A flexible simulator for aerial transportation and manipulation. *IEEE Transactions on Robotics* **40**, 831–850 (2024).

25. M. Bernard, K. Kondak, G. Hommel, Load transportation system based on autonomous small size helicopters. *The Aeronautical Journal* **114** (1153), 191–198 (2010).
26. M. Bernard, K. Kondak, I. Maza, A. Ollero, Autonomous transportation and deployment with aerial robots for search and rescue missions. *Journal of Field Robotics* **28** (6), 914–931 (2011).

27. E. J. Smeur, Q. Chu, G. C. De Croon, Adaptive incremental nonlinear dynamic inversion for attitude control of micro air vehicles. *Journal of Guidance, Control, and Dynamics* **39** (3), 450–461 (2016).

28. E. Tal, S. Karaman, Accurate tracking of aggressive quadrotor trajectories using incremental nonlinear dynamic inversion and differential flatness. *IEEE Transactions on Control Systems Technology* **29** (3), 1203–1218 (2020).

29. S. Sun, A. Romero, P. Foehn, E. Kaufmann, D. Scaramuzza, A comparative study of nonlinear mpc and differential-flatness-based control for quadrotor agile flight. *IEEE Transactions on Robotics* **38** (6), 3357–3373 (2022).

30. D. Mellinger, V. Kumar, Minimum snap trajectory generation and control for quadrotors, in *IEEE International Conference on Robotics and Automation* (2011), pp. 2520–2525.

31. P. Foehn, et al., Agilicious: Open-source and open-hardware agile quadrotor for vision-based flight. *Science robotics* **7** (67), eabl6259 (2022).

32. S. F. Hoerner, *Fluid-dynamic drag: practical information on aerodynamic drag and hydrodynamic resistance*, 2 ed. (1965).

33. A. Petitti, et al., Inertial estimation and energy-efficient control of a cable-suspended load with a team of UAVs, in *IEEE International Conference on Unmanned Aircraft Systems* (2020), pp. 158–165.

34. A. Richards, J. How, Robust stable model predictive control with constraint tightening, in *2006 American Control Conference (IEEE)* (2006), pp. 6–pp.

35. J. Köhler, M. A. Müller, F. Allgöwer, A novel constraint tightening approach for nonlinear robust model predictive control, in *IEEE American Control Conference* (2018), pp. 728–734.

36. P. Foehn, A. Romero, D. Scaramuzza, Time-optimal planning for quadrotor waypoint flight. *Science robotics* **6** (56), eabh1221 (2021).

37. H. G. Bock, K.-J. Plitt, A multiple shooting algorithm for direct solution of optimal control problems. *IFAC Proceedings Volumes* **17** (2), 1603–1608 (1984).

38. M. Diehl, et al., Real-time optimization and nonlinear model predictive control of processes governed by differential-algebraic equations. *Journal of Process Control* **12** (4), 577–585 (2002).
39. R. Verschueren, et al., Acados—a modular open-source framework for fast embedded optimal control. *Mathematical Programming Computation* **14** (1), 147–183 (2022).

40. M. Faessler, A. Franchi, D. Scaramuzza, Differential flatness of quadrotor dynamics subject to rotor drag for accurate tracking of high-speed trajectories. *IEEE Robotics and Automation Letters* **3** (2), 620–626 (2017).

41. S. Umeyama, Least-squares estimation of transformation parameters between two point patterns. *IEEE Transactions on Pattern Analysis & Machine Intelligence* **13** (04), 376–380 (1991).

42. D. Brescianini, R. D'Andrea, Tilt-prioritized quadrocopter attitude control. *IEEE Transactions on Control Systems Technology* **28** (2), 376–387 (2018).

43. D. Zou, P. Tan, W. Yu, Collaborative visual SLAM for multiple agents: A brief survey. *Virtual Reality & Intelligent Hardware* **1** (5), 461–482 (2019).

44. Y. Tian, et al., Kimera-multi: Robust, distributed, dense metric-semantic slam for multi-robot systems. *IEEE Transactions on Robotics* **38** (4) (2022).

45. T. Cieslewski, *Decentralized multi-agent visual SLAM*, Ph.D. thesis, University of Zurich (2021).

46. M. Peterson, P. Lusk, A. Avila, J. How, TCAFF: Temporal consistency for robot frame alignment. *arXiv preprint arXiv:2405.05210* (2025).

## Acknowledgments

The authors would like to thank Kseniia Khomenko, Maurits Pfaff, Dr. Gianluca Corsini, Riccardo Belletti, Dr. Dennis Benders, and Shantnav Agarwal for their support during the experiments. We also extend our gratitude to Dr. Yunlong Song, Dr. Guanrui Li, Prof. Guido de Croon, and Prof. Davide Scaramuzza for their valuable discussions and feedback.

**Funding:** This work was supported by the Dutch Research Council (NWO) through the Veni Talent Programme (Grant No. 20256, Accurate Aerial Manipulation), by the French National Research Agency (ANR) under Grant ANR-24-CE33-5799 (MATES), and by the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 871479 (AERIAL-CORE).
**Author contributions:** S.S. acquired funding, formulated the main ideas, developed the algorithm, conducted experiments, and wrote the manuscript. X.W. developed the algorithm, conducted experiments, and revised the manuscript. D.S helped with the algorithm implementation and revised the manuscript. A.F, M.T, J.A formulated ideas, acquired funding, and revised the manuscript.

**Competing interests:** There are no competing interests to declare.

**Data and materials availability:** All (other) data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials. The data for this study have been deposited in the database Dryad (https://doi.org/10.5061/dryad.n2z34tn6w).

**Supplementary materials**

Supplementary Methods

Supplementary Discussion

Figures S1 to S5

Table S1

Algorithm S1

Movies S1 to S5
# Supplementary Materials for
Agile and Cooperative Aerial Manipulation of a
Cable-Suspended Load

Sihao Sun*, Xuerui Wang, Dario Sanalitro,
Antonio Franchi, Marco Tognon, Javier Alonso-Mora

*Corresponding author. Email: sihao.sun@outlook.com

This PDF file includes:

Supplementary Methods

Supplementary Discussion

Figures S1 to S5

Table S1

Algorithm S1

Captions for Movies S1 to S5

Other Supplementary Materials for this manuscript:

Movies S1 to S5
# Supplementary Methods

## Proof of Smoothness of Quadrotor Trajectories

**Proposition 1:** When the i-th cable remains taut, the trajectory of i-th quadrotor, denoted as $\mathbf{p}_i(t)$, is $C^3$ smooth if $\lambda_i$ and $\gamma_i$ defined in Equation 3 are bounded.

**Proof:** For the i-th quadrotor, the kinematic constraint (Equation 5) holds when the corresponding cable is taut. Then we take the 3rd-order derivative of Equation 5 to obtain the jerk of the quadrotor

$$ \ddot{\mathbf{v}}_i = \ddot{\mathbf{v}} + \mathbf{R}(\mathbf{q}) \left\{ \omega \times \left[ \dot{\omega} \times \boldsymbol{\rho}_i + \omega \times (\omega \times \boldsymbol{\rho}_i) \right] + \ddot{\omega} \times \boldsymbol{\rho}_i + \dot{\omega} \times (\omega \times \boldsymbol{\rho}_i) + \omega \times (\dot{\omega} \times \boldsymbol{\rho}_i) \right\} - l_i \left\{ \ddot{\mathbf{r}}_i \times \mathbf{s}_i + 2\dot{\mathbf{r}}_i \times (\mathbf{r}_i \times \mathbf{s}_i) + \mathbf{r}_i \times (\dot{\mathbf{r}}_i \times \mathbf{s}_i) + \mathbf{r}_i \times \left[ \mathbf{r}_i \times (\mathbf{r}_i \times \mathbf{s}_i) \right] \right\}, \quad (S1) $$

where $\dot{\omega}$ is given in the load dynamics (Equation 2); $\ddot{\mathbf{v}}$ and $\ddot{\omega}$ are obtained by taking the derivative of both sides of Equation 2:

$$ \ddot{\mathbf{v}} = -\frac{1}{m} \sum_{i=1}^{n} \left[ \dot{t}_i \mathbf{s}_i + t_i (\mathbf{r}_i \times \mathbf{s}_i) \right], \quad (S2) $$

$$ \ddot{\omega} = \mathbf{J}^{-1} \left\{ -\dot{\omega} \times \mathbf{J}\omega - \omega \times \mathbf{J}\dot{\omega} + \sum_{i=1}^{n} \left[ \dot{t}_i \mathbf{R}^{\top} \mathbf{s}_i + t_i \left( -\omega \times \mathbf{R}^{\top} \mathbf{s}_i + \mathbf{R}^{\top} (\mathbf{r}_i \times \mathbf{s}_i) \right) \right] \times \boldsymbol{\rho}_i \right\}, \quad (S3) $$

Judging from the load-cable dynamics (Equation 2 and 3), the continuousness of $\ddot{\mathbf{v}}_i$ is determined by the highest-order states $\dot{t}_i$ and $\ddot{\mathbf{r}}_i$. Therefore, $\ddot{\mathbf{v}}_i$ is continuous; namely $\mathbf{p}_i$ is $C^3$-smooth, when $\ddot{t}_i = \gamma_i$ and $\ddot{\mathbf{r}}_i = \lambda_i$ are bounded. □

The OCP of the planner takes $\gamma_i$ and $\lambda_i$ as inputs, which can be bounded by setting input constraints. Hence, the generated trajectories of all quadrotors are smooth up jerk as long as the cable tautness is guaranteed. One step further, once the reference jerk is continuous, we can also obtain a smooth angular velocity reference.

**Proposition 2:** The angular velocity of the i-th quadrotor expressed in the inertial frame, denoted by $\omega_i^T \in \mathbb{R}^3$ is $C^0$-smooth if $\lambda_i$ and $\gamma_i$ defined in Equation 3 are bounded, and aerodynamic drag $f_{a,i}$ is at least $C^1$-smooth.

**Proof:** To obtain the angular velocity reference of each quadrotor, we need to revisit the translational dynamic equation of the i-th quadrotor

$$ \dot{\mathbf{v}}_i = (T_i z_i + t_i \mathbf{s}_i + \mathbf{f}_{a,i}) / m_i + \mathbf{g} \quad (S4) $$
Taking the derivative of both sides of Equation S4, we have

$$ 
\begin{aligned}
\mathbf{h}_i &\triangleq \omega_i^T \times z_i \\
&= \left[ m_i \ddot{v}_i - \dot{T}_i z_i - \dot{t}_i s_i - t_i (\mathbf{r}_i \times \mathbf{s}_i) - \dot{\mathbf{f}}_{a,i} \right] / T_i
\end{aligned} 
\quad (S5) 
$$

Since the yaw rate references are zero for all quadrotors, $\omega_i^T$ is perpendicular to $z_i$. Then we can obtain the expression of $\omega_i^T$ by

$$ \omega_i^T = z_i \times \mathbf{h}_i = \frac{1}{T_i} z_i \times \left[ \ddot{v}_i - \dot{t}_i s_i - t_i (\mathbf{r}_i \times \mathbf{s}_i) - \dot{\mathbf{f}}_{a,i} \right] \quad (S6) $$

According to Proposition 1 and Equation 9, $\ddot{v}_i$ is $C^0$-smooth and $T_i$ is at least $C^1$ smooth. When $\gamma_i$ is bounded, $\dot{t}_i$ is also $C^0$-smooth. Hence angular velocity $\omega_i^T$ is also $C^0$-smooth when $\mathbf{f}_{a,i}$ is at least $C^1$-smooth. □

If we use Equation S6 to generate the angular velocity reference of each quadrotor, its smoothness is guaranteed through Proposition 2, if we use a smooth drag model (In this work, we assume zero drag for simplicity). The smooth angular velocity reference, used as feed-forward terms by the trajectory tracking controller onboard the quadrotor, guarantees smooth quadrotor behavior, which is particularly crucial during dynamic motions.
Incremental Nonlinear Dynamic Inversion (INDI) Low-Level Controller

The INDI low-level controller of each quadrotor generates rotor speed commands, using the collective thrust command $T_{des}$ and angular acceleration command $\alpha_{des}$. Here, we summarize the key equations of the INDI controller introduced in our previous work (29). In the following context, we denote rotor speed commands as $\mathbf{u}_c \in \mathbb{R}^4$, and rotor speed measurement as $\mathbf{u}_m \in \mathbb{R}^4$. Note that the following equations apply to a single quadrotor. Hence we omit the subscript $i$ for readability.

The INDI low-level controller employs the following model that maps the rotor speeds to the collective thrust $T \in \mathbb{R}_{\ge 0}$ and body torque $\tau \in \mathbb{R}^3$

$$ \begin{bmatrix} T \\ \tau \end{bmatrix} = G_1 \mathbf{u}_m^{\circ 2} + G_2 \dot{\mathbf{u}}_m \quad (S7) $$

where $G_1$ and $G_2$ are the control effectiveness matrices with respect to the rotor speeds. Specifically, $G_1$ depends on the shape and size of the quadrotor and aerodynamic coefficients of the propellers. $G_2$ captures the inertial yawing torque due to the acceleration and deceleration of the rotors, which is a function of the moment of inertia of the rotors.

Therefore, once the desired collective thrust and torque $[T_{des}, \tau_{des}]^\top$ is computed, INDI numerically solves the following equation to obtain rotor speed command $\mathbf{u}_c$

$$ \begin{bmatrix} T_{des} \\ \tau_{des} \end{bmatrix} = G_1 \mathbf{u}_c^{\circ 2} + \Delta t^{-1} G_2 (\mathbf{u}_c - \mathbf{u}_{c,k-1}) \quad (S8) $$

where $\Delta t$ is the sampling interval of the controller; $\mathbf{u}_{c,k-1}$ is the last computed rotor speed command $\mathbf{u}_c$. And $\mathbf{u}_{c,0} = \mathbf{u}_m$ for initialization.

Unlike conventional dynamic inversion, the INDI low-level controller defines the desired body torque in the following incremental form

$$ \tau_{des} = \tau_f + J (\alpha_{des} - \dot{\omega}_f) \quad (S9) $$

where $\dot{\omega}_f$ is the angular acceleration obtained by numerically differentiating the filtered gyroscope measurement from the quadrotor. $\tau_f$ is the filtered body torque, which can be calculated using rotor speed measurements leveraging Equation S7, yielding

$$ \tau_f = [G_1 \mathbf{u}_f^{\circ 2} + \Delta t^{-1} G_2 (\mathbf{u}_f - \mathbf{u}_{f,k-1})]_{k=4} \quad (S10) $$

where $\mathbf{u}_f$ is the low-pass filtered rotor speed measurements. Note that the cutoff frequency of the filter for $\mathbf{u}_f$ and $\omega_f$ is the same to synchronize the delay introduced by the low-pass filter on these two measurements. The INDI low-level controller leverages the sensor measurements to effectively capture and compensate for the external torques that are not captured in Equation S7, such as aerodynamic torque, CoG bias of quadrotors, etc.
# Supplementary Discussions

## Load Pose Estimation Performance

In the above experiments (all real-world experiments shown in this paper), we need information on load pose and twist to achieve dynamic and accurate trajectory tracking. In the state-of-the-art method that includes real-world experiments, additional sensors are required for the load pose estimation. The most commonly used approach is attaching reflective markers on the load to measure its pose from the motion capture system (12, 16, 21), or resorting to additional downward-facing cameras and attaching additional circular tags on the load (11). However, it is impractical to attach these sensors in the field for day-to-day operations.

In comparison, our algorithm does not need to put any additional sensors on the load, nor does it make any modifications to the quadrotors. In the experiments, we demonstrated that by simply leveraging the IMU on each quadrotor and the dynamic model of the multi-lifting system can provide a sufficiently accurate load pose and twist estimate as well as the cable states to achieve agile pose control. Fig. S1 presents the comparison between the ground truth pose of the load and the estimated pose while tracking reference Fast. Despite the large acceleration of the motion, over 45 degrees of inclination, and continuous yawing motion, our method provided sufficiently high estimation accuracy to achieve closed-loop trajectory tracking. The position estimation RMSE was 0.136 m and the attitude estimation RMSE was 7.5 deg even with the highly dynamic motion of the system.
<figure><img src="image_11.png" /><figcaption>**Figure S1: Load pose estimation and tracking result.** Without adding additional sensors to the load or force sensors to quadrotors, our algorithm could accurately estimate the load's pose for an accurate closed-loop control. **Top:** time history of the position and attitude (Euler angles) estimate in comparison with the ground truth from the motion capture system, as well as the reference load pose while tracking the figure-eight trajectory Fast; **Bottom:** time history of the estimated, the ground truth, and the reference load pose while tracking a setpoint.</figcaption></figure>
## Sensitivity to Quadrotor Coordinate Frame Misalignment

We conducted Monte Carlo simulations to examine the effect of misalignment among ground-fixed reference frames (e.g., odometry frames) used by different quadrotors. This issue arises when no global positioning sensors (e.g., motion capture or GPS) are available, such as when using Visual Inertial Odometry (VIO) for state estimation in GPS-denied environments. In this simulation, we assumed that the initial odometry frames of all quadrotors were well-calibrated and aligned with the inertial frame. We then introduced the transformation between the odometry frame and the inertial frame for each quadrotor, in the form of a random-walk process. In this way, we simultaneously simulated the misalignment of reference frames among quadrotors and the pose drift typically observed in visual (-inertial) odometry algorithms.

It is worth noting that extensive research has been conducted on aligning estimated reference frames among multiple robots in the context of multi-robot VIO (e.g., (43–46)). A standard approach is to align these coordinate frames in real-time through place recognition to match landmarks seen by different cameras (quadrotors in our case) and by estimating the relative poses between coordinate frames (45). The accuracy of relative pose estimation depends on the quality frontend / backend in the VIO algorithms and the quality of the established map. A typical example alignment frequency is 1 Hz (46).

Since we were unable to simulate all possible approaches in our study, we only performed a worst-case and best-case scenario analysis. In the worst case, we ran simulations without any real-time alignment. In the best case, we ran simulations and precisely aligned the quadrotors’ odometry frames with the inertial frame at 1 Hz. In both cases, we let the system follow the reference trajectory Fast.

<figure><img src="image_12.png" /><figcaption>**Figure S2: Simulation result under misalignment between quadrotor reference frames.** The misalignment was induced by random-walk noise at different levels, on the transform between quadrotor VIO frames and the inertial frame. (A) The worst case: there was no frame alignment mechanism among quadrotors. (B) The best case: the frames were precisely aligned at 1 Hz.</figcaption></figure>
Fig. S2 presents the results of the two cases. It shows that, in the worst case, without any alignment between quadrotors, our method still maintains a success rate of 90% at the 2nd noise level, where the standard deviation of coordinate frame drift is 0.03 m on position and 3 deg on attitude after 50 seconds. We noticed that the attitude misalignment between the coordinate frames is the main cause of failure, as it introduced large position state estimation errors of quadrotors when they follow a trajectory that is far from the origin of the frames. On the other hand, we demonstrated that with an ideal alignment running at 1 Hz, the performance is almost unaffected by the drift of coordinate frames until it reaches noise level 4. These results have demonstrated the strong potential of our method to combine with multi-robot VIO algorithms and deploy in a GPS-denied environment.
<figure><img src="image_13.png" /><figcaption>**Figure S3: Experimental setup.** A snapshot of our experiment, together with the parameters of the load, the quadrotors, and the cables.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure S4: Definition of reference frames and symbols.** $F_I, F_L, F_i$ respectively denote the inertial frame, load-fixed frame, and the *i*-th quadrotor-fixed frame, where the *x*, *y*, *z* are marked in red, green, blue colors respectively.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure S5: Step response under load model uncertainties and communication delay.** Simulation result comparing the tracking performance between our method and the two baseline methods (NMPC (12) and Geometric (14)) under various types of model mismatch on the load, as well as communication delay between the centralized planner and quadrotors. We sent a load reference point at 2 m away along y-axis, and an attitude command of -90 deg, -20 deg, and 30 deg, respectively, on yaw, pitch, and roll. Our method clearly outperforms the two baseline methods in the presence of load model mismatch and communication delays.</figcaption></figure>
**Algorithm S1** Iterative Kabsch–Umeyama algorithm to initialize states of the EKF

**Input:** $n, \mathbf{p}_i, \rho_i, l_i$ for $i \in \{1, 2, \dots, n\}$

Define tolerance and maximum steps for iteration: $tol_{\text{pos}}$, $tol_{\text{att}}$, $iter_{\text{max}}$

Define average cable connection points $\bar{\rho} = \sum_{i=1}^{n} \rho_i/n, \mathbf{L} = [\rho_1 - \bar{\rho}, \rho_2 - \bar{\rho}, \dots, \rho_n - \bar{\rho}]$

Define initial load pose $\mathbf{p} = [0, 0, 0]^{\top}, \mathbf{R} = \mathbf{I}_3$

Initialize the last load pose $\mathbf{p}_{\text{last}} = \infty, \mathbf{R}_{\text{last}} = \mathbf{O}_3$

Initial guess of cable directions $s_i = [0, 0, -1]^{\top}$ for $i \in \{1, 2, \dots, n\}$

**for** $k = 1, \dots, iter_{\text{max}}$ **do**

**for** $i = 1, \dots, n$ **do**

$$c_i = \mathbf{p}_i + s_i l_i$$

$$\bar{c} = \sum_{i=1}^{n} c_i / n$$

$$C = [c_1 - \bar{c}, c_2 - \bar{c}, \dots, c_n - \bar{c}]$$

$$[U, V] = \text{SVD}(LC^{\top})$$

$$\mathbf{R} = \mathbf{V} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \text{sign}(\det(V\mathbf{U}^{\top})) \end{bmatrix} \mathbf{U}^{\top} \qquad \triangleright \text{ Estimated load attitude}$$

**for** $i = 1, 2, \dots, n$ **do**

$$\tilde{p}_i = c_i - \mathbf{R} \rho_i$$

$$\mathbf{p} = \sum_{i=1}^{n} \tilde{p}_i / n$$

► Estimated load position

**for** $i = 1, 2, \dots, n$ **do**

$$s_i = (\mathbf{R} \rho_i + \mathbf{p} - \mathbf{p}_i) / \| \mathbf{R} \rho_i + \mathbf{p} - \mathbf{p}_i \|$$

► Estimated cable direction

**if** $\|\mathbf{p} - \mathbf{p}_{\text{last}}\| < tol_{\text{pos}}$ **and** $\|\det(\mathbf{R}^{-1}\mathbf{R}_{\text{last}}) - 1\| < tol_{\text{att}}$ **then**

**Break**

$$\mathbf{p}_{\text{last}} = \mathbf{p}, \mathbf{R}_{\text{last}} = \mathbf{R}$$

**Return** $\mathbf{p}, \mathbf{q}(\mathbf{R}), s_1, s_2, \dots, s_n$
<table><caption>Table S1: Algebraic expression of load reference trajectories</caption><thead><tr><th>Name</th><th>p<sub>x</sub></th><th>p<sub>y</sub></th><th>p<sub>z</sub></th></tr></thead><tbody><tr><td>Slow</td><td>2.5 cos (0.25t)</td><td>2 sin (0.5t)</td><td>1.0</td></tr><tr><td>Medium</td><td>2.5 cos (0.5t)</td><td>2 sin (t)</td><td>1.0</td></tr><tr><td>Medium Plus</td><td>cos (t)</td><td>sin (2t)</td><td>1.0</td></tr><tr><td>Fast</td><td>2.5 cos (t)</td><td>2 sin (2t)</td><td>1.0</td></tr></tbody></table>
**Caption for Movie S1. Video comparing our method with baseline methods in simulation.**
The video shows the 6-DoF simulation result of our method and the two baseline methods in tracking the two reference trajectories (Fast and Slow). The video is associated with the result given in Table 1.

**Caption for Movie S2. Video of flight for obstacle avoidance.** The video shows real-world experiments performing two obstacle avoidance tasks: one through a narrow passage between two walls, and the other through a horizontally oriented gap.

**Caption for Movie S3. Video of flight with four quadrotors.** The video shows the real-world experiment of our method scaling up to a case with four quadrotors following a reference trajectory dynamically.

**Caption for Movie S4. Video of flight in windy conditions.** The video shows that our method effectively controls the system to follow trajectories under moderate wind disturbances.

**Caption for Movie S5. Video showing results with large quadrotor state estimation errors.**
The video shows real-world experiment and simulation results under different levels of quadrotor state estimation errors.