In recent years, the application of unmanned aerial vehicles (UAVs) has expanded significantly across various fields, including military operations, agriculture, scientific research, and disaster management. Among these, multirotor drones have gained prominence due to their excellent hovering capabilities, relatively simple control logic, and substantial load-carrying capacity. However, traditional single-drone payload transport systems are limited in their carrying capacity, making multi-drone collaboration an emerging trend. This paper addresses the challenge of cooperative payload suspension using multiple multirotor drones, focusing on modeling and control strategies to enhance stability and efficiency.
The conventional approach to modeling rotor-load systems often employs the Euler-Lagrange method, which requires introducing multiple auxiliary variables and complicates the modeling process. In contrast, we utilize the Udwadia-Kalaba (U-K) equation to handle constraints in the multirotor drone system, providing an explicit solution for constrained dynamics. By establishing kinematic and dynamic models of the multirotor drone, we apply the U-K equation to decompose the effects of the payload into horizontal and vertical directions, treating them as step inputs for computational simplicity. Our control strategy employs a cascade PID control model to ensure robustness and stability under cooperative conditions. Simulations conducted in Matlab/Simulink demonstrate that the proposed method achieves rapid and accurate formation stability, optimizing the collaborative control of multirotor drones.
Introduction to Multi-Rotor Drone Systems
Multirotor drones, particularly quadrotors, are characterized by their ability to perform vertical take-off and landing (VTOL) and maintain stable hover. The X-frame and cross-frame configurations are common, with the X-type offering better balance due to its centered gravity. In this study, we adopt the X-frame multirotor drone for its superior stability in collaborative tasks. The coordinate system follows the right-hand rule, with the body frame and spatial frame related through rotation matrices.
The dynamics of a multirotor drone involve multiple models: the power unit model, control efficiency model, and rigid-body model. The rigid-body model comprises kinematic and dynamic components. The kinematic model describes the position and velocity in the spatial frame, while the dynamic model accounts for forces and moments. The rotation matrix $$ R_e^b $$ transforms vectors from the body frame to the spatial frame and is given by:
$$ R_e^b = \begin{bmatrix} \cos\theta \cos\psi & \cos\psi \sin\theta \sin\phi – \sin\psi \cos\phi & \cos\psi \sin\theta \cos\phi + \sin\psi \sin\phi \\ \cos\theta \sin\psi & \sin\phi \sin\theta \sin\phi + \cos\psi \cos\phi & \sin\psi \sin\theta \cos\phi – \cos\psi \sin\phi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi \cos\theta \end{bmatrix} $$
where $$ \psi $$, $$ \theta $$, and $$ \phi $$ represent the yaw, pitch, and roll angles, respectively. The position kinematics are described by $$ \dot{P}^e = V^e $$, and the attitude kinematics by $$ \dot{\Theta}^e = W \cdot \omega^b $$, where $$ W $$ is a parameter matrix and $$ \omega^b $$ is the angular velocity in the body frame.
The dynamics of the multirotor drone are derived from Newton’s second law. The position dynamics are:
$$ m \dot{V}^e = mg – f^e $$
where $$ f^e = R_e^b f^b $$ is the lift force in the spatial frame. Simplifying, we get:
$$ \dot{V}^e = g – R_e^b \frac{f^b}{m} $$
The attitude dynamics, ignoring gyroscopic moments due to their minimal impact, are:
$$ \tau = J \dot{\omega}^b $$
where $$ J $$ is the moment of inertia. The control efficiency model relates the rotor speeds to the forces and moments:
$$ \begin{bmatrix} f^b \\ \tau_x \\ \tau_y \\ \tau_z \end{bmatrix} = \begin{bmatrix} c_T & c_T & c_T & c_T \\ -\frac{\sqrt{2}}{2} d c_T & \frac{\sqrt{2}}{2} d c_T & \frac{\sqrt{2}}{2} d c_T & -\frac{\sqrt{2}}{2} d c_T \\ \frac{\sqrt{2}}{2} d c_T & -\frac{\sqrt{2}}{2} d c_T & \frac{\sqrt{2}}{2} d c_T & -\frac{\sqrt{2}}{2} d c_T \\ c_M & c_M & -c_M & -c_M \end{bmatrix} \begin{bmatrix} T_1^2 \\ T_2^2 \\ T_3^2 \\ T_4^2 \end{bmatrix} $$
Here, $$ d $$ is the fuselage radius, $$ c_T $$ and $$ c_M $$ are the thrust and moment coefficients, and $$ T_i $$ are the rotor speeds.
Udwadia-Kalaba Equation for Constrained Systems
The U-K equation provides a straightforward method for modeling constrained systems without introducing auxiliary variables. For a system with $$ n $$ degrees of freedom, the unconstrained dynamics are:
$$ M \ddot{q}_u = F(q_u, \dot{q}_u, t) $$
where $$ M $$ is the symmetric positive definite mass matrix, $$ q_u $$ is the unconstrained generalized coordinate, and $$ F $$ is the generalized force. When constraints are applied, such as those from payload suspension, the system becomes:
$$ M \ddot{q} = F + F_c $$
where $$ F_c $$ is the constraint force. The U-K equation expresses $$ F_c $$ explicitly:
$$ F_c = M^{1/2} (V M^{-1/2})^+ (W – V \dot{q}_u) $$
Here, $$ V $$ and $$ W $$ are matrices derived from the constraints, and $$ ^+ $$ denotes the Moore-Penrose pseudoinverse. The constrained acceleration is:
$$ \ddot{q} = \ddot{q}_u + M^{-1/2} (V M^{-1/2})^+ (W – V \dot{q}_u) $$
For a multirotor drone system carrying a payload via $$ n $$ cables, the constraint ensures the distance between each drone and the payload remains constant. Defining $$ L_i = \eta_i – \eta_l $$ as the vector from the payload to drone $$ i $$, the constraint is:
$$ g_i = \| L_i \|^2 – d_i^2 = 0 $$
where $$ d_i $$ is the length of cable $$ i $$. The matrices $$ V $$ and $$ W $$ are constructed as:
$$ V_{j,1 \times 3(n+1)} = 2 L_j^T [ \mathbf{0}_{3 \times 3(j-1)} \quad I \quad \mathbf{0}_{3 \times 3(n-j)} \quad -I ] $$
$$ W_{j,1 \times 1} = -2 \dot{L}_j^T \dot{L}_j $$
The mass matrix is $$ M = \text{diag}(m_1, m_2, \dots, m_n, m_l) $$, where $$ m_i $$ are the drone masses and $$ m_l $$ is the payload mass.
Payload Equivalence and Decomposition
To simplify the model, the payload’s effect on the multirotor drone is decomposed into horizontal and vertical components. During takeoff, the tension in the cables causes a sudden vertical force and a horizontal force pulling the drones toward the formation center. The horizontal force can be equated to a negative pitch or roll angle, depending on the drone’s position in the formation.
For a square formation, the horizontal equations are:
$$ f_T \cos \beta = (m + \frac{1}{4} m_l) g $$
$$ f_T \sin \beta = \frac{1}{4} m_l g $$
where $$ f_T $$ is the total thrust after payload attachment, and $$ \beta $$ is the angle between the thrust vector and the vertical. This equivalence allows us to treat the payload influence as a step input in the control system.
Cascade PID Controller Design
The multirotor drone is an underactuated system, with four inputs (total thrust and three moments) controlling six outputs (position and orientation). We design a cascade PID controller to manage this complexity. The outer loop regulates position, while the inner loop controls attitude, ensuring rapid response and disturbance rejection.
The attitude controller uses cascade PID for angle (outer loop) and angular velocity (inner loop). The angular velocity loop outputs the desired moments:
$$ \tau_d = K_{\omega p} e_\omega + K_{\omega i} \int e_\omega \, dt + K_{\omega d} \dot{e}_\omega $$
where $$ e_\omega $$ is the angular velocity error, and $$ K_{\omega p} $$, $$ K_{\omega i} $$, $$ K_{\omega d} $$ are PID gains.
The position controller handles horizontal and vertical channels separately. For horizontal control, the desired pitch and roll angles are computed from the position error:
$$ \Theta_h = -g^{-1} A_\psi^{-1} \dot{V}_h $$
where $$ A_\psi = \begin{bmatrix} \sin\psi & \cos\psi \\ -\cos\psi & \sin\psi \end{bmatrix} $$. The vertical controller calculates the desired thrust:
$$ m \dot{V}_z = -f_b + mg $$
$$ \dot{V}_z = K_{pz} e_z + K_{iz} \int e_z \, dt + K_{dz} \dot{e}_z $$
with $$ e_z $$ being the height error.
Formation Strategy for Cooperative Transport
We employ a square formation for the multirotor drones, with each drone at a vertex and the payload at the center. This symmetric configuration ensures balanced force distribution. The cable length $$ l $$ is fixed to maintain a constant angle $$ \alpha $$ between the drones and the payload, satisfying:
$$ l \tan \alpha = \Delta x = \Delta y $$
A virtual leader point at the formation center generates trajectory commands for the group. The trajectory is planned using minimum snap (MS) methods, with waypoints and zero initial/final velocities and accelerations. The trajectory is represented as a piecewise polynomial:
$$ p(t) = \begin{cases} [1, t, t^2, \dots, t^n] \cdot p_1 & t_0 \leq t < t_1 \\ [1, t, t^2, \dots, t^n] \cdot p_2 & t_1 \leq t < t_2 \\ \vdots & \vdots \\ [1, t, t^2, \dots, t^n] \cdot p_k & t_{k-1} \leq t < t_k \end{cases} $$
This approach ensures smooth and stable motion for the multirotor drone formation.
Stability Analysis
The stability of the multirotor drone system is analyzed using the state-space representation. The state vector is:
$$ X = [\dot{x}, \dot{y}, \dot{z}, x, y, z, \dot{\phi}, \dot{\theta}, \dot{\psi}, \phi, \theta, \psi, g]^T $$
and the input vector is $$ U = [f_b, \tau_x, \tau_y, \tau_z]^T $$. The transfer function $$ G(s) = C(sI – A)^{-1}B + D $$ is derived for each channel. For the vertical channel, with simplified dynamics ($$ \theta \approx \phi \approx 0 $$), the transfer function is:
$$ G_f(s) = \frac{1.19}{s^2 (0.02s + 1)} $$
With cascade PID control, the closed-loop transfer function becomes:
$$ \Phi_f(s) = \frac{3.78s^2 + 21.4s + 24.2}{0.02s^4 + s^3 + 5.67s^2 + 32.1s + 37.8} $$
The Routh-Hurwitz criterion confirms stability, as all coefficients in the first column of the Routh array are positive. Similar analyses apply to other channels, ensuring overall system stability for the multirotor drone formation.
Simulation Results and Discussion
We simulate a scenario with four multirotor drones in a square formation, starting at positions (0.5, 0, 0), (0, 0.5, 0), (0, -0.5, 0), and (-0.5, 0, 0). The leader point moves from (0, 0, 0) to (5, 5, 5), with waypoints at (2, 2) and (4, 4). The total trajectory time is 10 seconds. The drones first ascend to the desired height, then move horizontally to maintain formation.
The position trajectories show that the multirotor drones achieve the target within 10 seconds while preserving the square formation. Velocity profiles converge to zero, indicating stable hover at the destination. The attitude angles adjust to compensate for payload effects: for example, drone 1 settles at a pitch of -10°, while drone 2 has a roll of 10°. These adjustments generate horizontal forces to counter payload pull.
The following table summarizes key performance metrics for the multirotor drone formation:
| Drone | Final Position (x, y, z) | Settling Time (s) | Max Attitude Angle (°) |
|---|---|---|---|
| UAV1 | (5.5, 5, 5) | 9.8 | -10 (pitch) |
| UAV2 | (5, 5.5, 5) | 9.7 | 10 (roll) |
| UAV3 | (5, 4.5, 5) | 9.9 | -10 (roll) |
| UAV4 | (4.5, 5, 5) | 9.8 | 10 (pitch) |
The Euler angle variations demonstrate the controller’s ability to handle disturbances. The cascade PID ensures that attitude changes are smooth and convergent. Compared to prior work using LQR controllers, our approach with U-K equations and cascade PID achieves faster stabilization (within 10 seconds vs. 15 seconds), enhancing both stability and rapidity for multirotor drone systems.
Conclusion and Future Work
This paper presents a comprehensive framework for collaborative payload suspension using multirotor drones. By integrating U-K equations for constraint handling and cascade PID for control, we achieve stable and efficient formation transport. The square formation and trajectory planning ensure balanced load distribution and smooth motion. Simulations validate the system’s performance, with rapid convergence and minimal oscillation.
Future work will focus on refining the model to include more complex dynamics, such as aerodynamic effects and variable cable lengths. Incorporating consensus-based algorithms could improve formation robustness. Experimental validation with physical multirotor drones is also planned to verify practical applicability. Overall, this research contributes to advancing multi-drone cooperation, enabling larger payload capacities and broader applications for multirotor drone technologies.