Transportation of Payload Using Multiple Quadrotors via Rigid Connection

Abstract

Due to the limited payload capability of an aerial robot, multiple quadrotors can be used to manipulate payloads in aerial transportation, construction, and assembly tasks. This paper focuses on the cooperative transportation of a payload rigidly attached to multiple quadrotor bodies. These quadrotors may have different orientations. The dynamics equation of a rigid body in 3-D space is derived to describe the motion of such a transportation system. Robust position and attitude controllers are designed to drive the system to the desired pose. To assign control signals for each quadrotor, the control command allocation method compatible with the case that partial or all quadrotors are in parallel planes is developed. Finally, experimental results are presented to validate the effectiveness of the proposed controllers and control command allocation methods. Different from classical works in this field, this paper can solve the dynamics modeling, controller design, and control command allocation problems for the transportation of a rigidly connected payload using a team of quadrotors with different orientations.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Due to the advancement in unmanned aerial vehicles and the capabilities of vertical take-off and landing, quadrotors have been implemented in various applications, such as wildfire mapping, search and rescue, payload delivery, and agricultural surveys [1, 2]. As stated in Ref. [3], the market size of the global UAV-assisted logistics and transportation will grow from 5.3 billion dollars in 2019 to about 11 billion dollars by 2026. Several research groups have developed notable algorithms and presented experimental results to prove the quadrotor’s ability in payload transportation [4–9]. However, an individual quadrotor usually has limited payload capability. One promising method to address this limitation is to transport a heavy payload using a team of quadrotors. Compared with transportation using a single quadrotor, cooperative transportation has to consider the collaboration and synchronization among quadrotors and the stability of the entire system [4, 10, 11]. Generally, there are mainly two carrying strategies, i.e., connecting the payload to the quadrotor bodies using cables or rigid fixtures. The cooperative transportation of a cable-suspended payload has been studied extensively [12–18]. For example, Sanalitro et al. designed a controller for aerial transportation using the minimum number of quadrotors and cables considering some system uncertainties [14]. Geng and Langelaan presented a centralized load-leading control method for the transportation of a slung payload using multiple rotorcraft to drive the payload to track the desired trajectory [15]. Jiang and Kumar focused on transportation using three aerial robots and proposed an analytic algorithm to solve the possible solution to the kinematics problem based on dialytic elimination [16]. Qian and Liu adopted Kane’s method to build the dynamics equation for multiple quadrotors carrying a slung payload and designed a cascade controller to drive the payload to follow the desired path [18]. Generally, it is convenient to attach payload using cables, and the long distances from quadrotors to the payload result in a negligible influence on the vehicles’ aerodynamics. However, the cable-suspended payload introduces additional dynamics. It is usually not possible to add more sensors and actuators to control the oscillations of the payload. Hence, the number of underactuated degrees of freedom increases, and the possible cable slackness will introduce system uncertainties, both of which pose a great challenge in its dynamics modeling and controller design.

To enable direct control of the payload’s position and attitude, some rigidly attached methods, such as the manipulators [19–21], graspers [22], and permanent electromagnets [23], have been proposed, and the resulting dynamics, navigation, and control problems have been studied by some scholars. For example, Lee et al. proposed a framework with estimation, control, and path planning for payload transportation using several aerial manipulators based on decentralized dynamics [19]. Mellinger et al. developed control algorithms for a team of quadrotors to grasp and transport a payload in 3-D space based on a dynamics model for the entire system [22]. Loianno and Kumar presented the basic dynamics models, estimation algorithms, and control methods of carrying a payload using multiple quadrotors via permanent electromagnets [23]. To reduce the system complexity, this work mainly focuses on the case with a simple rigidly attached method, such as the permanent magnets. In [22, 23], all $n$ quadrotors are placed in parallel planes, i.e., they are in a special configuration. Essentially, the entire system in such a case can be considered a classical multirotor aircraft with $4n$ propellers. The relative positions of these $4n$ propellers vary with different placements of $n$ quadrotors. Hence, the dynamics model of such a transportation system is the same as that of the classical quadrotor, but the command allocation method is different.

However, if the payload surface to be attached is not planar, all propellers of quadrotors may not be in parallel planes, such as the payloads with triangular and quadrangle crosssections shown in Figure 1. In such cases, the orientations of quadrotors are different, i.e., the quadrotors are in a more general configuration than those in [22, 23]. Furthermore, for quadrotors with different orientations, their thrust forces may have some horizontal components, i.e., the horizontal movement of the entire system can be achieved by just increasing or decreasing the thrusts of some quadrotors, which may provide a fast response in the horizontal movement of the entire system. In [24], Ritz and D’Andrea studied the cooperative transportation of a flexible ring using multiple quadrotors with different orientations and some assumptions, such as each quadrotor only providing two control inputs: a force in $z$ direction and a torque in roll direction in its body frame, to simplify the transportation problem. The physical coupling between quadrotors caused by the flexible payload is smaller than that in the case of the transportation of rigid bodies. Hence, essentially, the quadrotors are just considered actuators in [24]. A general framework for the dynamics analysis, controller design, and control command allocation method for the transportation of a rigid payload using multiple quadrotors in different configurations has not been studied yet. This work is aimed at providing a solution to such an open problem. Note that the relative orientations of the quadrotors in a general configuration are arbitrary under the condition that the lift force generated by quadrotors is larger than the gravity of the entire system. That is, the possible quadrotor configurations of the transportation system in this study have a wider range than those in [22, 23], and the systems in [22, 23] can be called a trivial case. In summary, the main contributions of this study are the following:

(i) The transportation system with quadrotors in various configurations is studied for the first time

(ii) A new dynamics modeling method and the corresponding control algorithm are presented in this work

(iii) A control command allocation method compatible with the case that some or all propellers are in parallel planes is developed

[figure(s) omitted; refer to PDF]

The remainder of this paper is organized as follows. Section 2 presents the dynamics model of the entire transportation system. The designed control law and control command allocation method are given in Section 3. Some experimental results are shown in Section 4. Conclusions are drawn in Section 5.

2. Transportation System Modeling

Suppose that, as shown in Figure 2, there are $n$ identical quadrotors attached to the payload rigidly. Assume that the payload will not affect the aerodynamic forces generated by quadrotors’ propellers. $n+2$ frames, i.e., the inertial frame $XYZ$, the body frame $xyz$ of the entire system, and the body frame $x_i{y}_i{z}_i$ of quadrotor $i$ for $i=1,2,\cdots ,n$, are built to describe the system. The inertial frame is the Earth-fixed West-South-Up frame. The origin of the frame $xyz$ is at the mass center $C$ of the whole system consisting of the payload and the $n$ quadrotors and the frame $xyz$ is parallel with the inertial frame at the initial time. The position vector of $C$ is denoted as ${\mathbf{r}}_c={x_c,y_c,z_c}^T$ in the inertial frame. The attitude of the entire system is described by the rotation matrix $\mathbf{R}\in SO3$ from the frame $xyz$ to the inertial frame and $\mathbf{R}0={\mathbf{I}}_3$ holds. Note that the Lie group $SO3=\mathbf{R}\in {ℝ}^{3\times 3}{\mathbf{R}}^T\mathbf{R}=\mathbf{R}{\mathbf{R}}^T={\mathbf{I}}_3,\det \mathbf{R}=1$ represents the group of $3\times 3$ orthogonal matrices with a determinant of one. In the body frame of the entire system, the position of the mass center of the $i$th quadrotor is represented by ${\mathbf{r}}_i={x_i,y_i,z_i}^T$. The rotation matrix from the body frame of the $i$th quadrotor to the frame $xyz$ is denoted by ${\mathbf{R}}_i\in SO3$. Since $n$ quadrotors are attached to the payload rigidly, both ${\mathbf{r}}_i$ and ${\mathbf{R}}_i$ are constant.

[figure(s) omitted; refer to PDF]

Essentially, the entire system shown in Figure 2 can be considered an unmanned aerial vehicle driven by 4$n$ propellers. Since all quadrotors are in a wide-range configuration, i.e., all propellers may not be in a plane, a three-dimensional control force vector may be generated. Therefore, the dynamics equations of the payload together with these $n$ quadrotors can be written as follows: $$\begin{array}{}\text{(1)}& {\dot{\mathbf{r}}}_c={\mathbf{v}}_c,\text{(2)}& m{\dot{\mathbf{v}}}_c=-mg{\mathbf{z}}_I+\mathbf{R}\mathbf{f}+{\mathbf{d}}_1,\text{(3)}& \dot{\mathbf{R}}=\mathbf{R}{\mathbf{\Omega }}^{\times },\text{(4)}& \mathbf{J}\dot{\mathbf{\Omega }}=-{\mathbf{\Omega }}^{\times }\mathbf{J}\mathbf{\Omega }+\tau +{\mathbf{d}}_2,\end{array}$$

where $m=m_p+n\times m_q$ is the total mass of the whole system, in which $m_p$ is the payload mass and $m_q$ is the quadrotor mass. $\mathbf{J}$ is the total inertia matrix of the entire system, whose calculation method based on the inertia matrices of the payload and quadrotor and the relative poses of quadrotors with respect to the frame $xyz$ can be found in Appendix A. ${\mathbf{v}}_c$ is the velocity vector of point $C$ in the inertial frame. $g$ is the constant gravitational acceleration. ${\mathbf{z}}_I={0,0,1}^T$. $\mathbf{\Omega }={{\omega }_1,{\omega }_2,{\omega }_3}^T$ is the angular velocity in the body frame $xyz$. Note that for any 3-dimensional vector $\mathbf{z}={z_1\hspace{1em}-z_2\hspace{1em}{z}_3}^T\in {ℝ}^3,$${\mathbf{z}}^{\times }$ is defined as the skew-symmetric matrix $0\hspace{1em}-z_3\hspace{1em}{z}_2;z_3\hspace{1em}0\hspace{1em}-z_1;-z_2\hspace{1em}{z}_1\hspace{1em}0$. ${\mathbf{d}}_1$ and ${\mathbf{d}}_2$ are unmodeled terms. $\mathbf{f}$ and $\tau$ are the total force and torque vectors provided by these $n$ quadrotors expressed in the system body frame $xyz$. Their detailed definitions are $$\begin{array}{}\text{(5)}& \begin{array}{c}\mathbf{f}\\ \tau \end{array}=\sum _{i=1}^n{\mathbf{A}}_i{\mathbf{u}}_i=\mathbf{A}\mathbf{u},\end{array}$$where $$\begin{array}{}\text{(6)}& {\mathbf{A}}_i=\begin{array}{cc}{\mathbf{R}}_i{\mathbf{z}}_{bi}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{r}}_i^{\times }{\mathbf{R}}_i{\mathbf{z}}_{bi}& {\mathbf{R}}_i\end{array}\in {ℝ}^{6\times 4},\end{array}$$${\mathbf{u}}_i={f_i,{\tau }_i^T}^T\in {ℝ}^4$ is the vector consisting of the thrust force and the 3-dimensional control torques of the $i$th quadrotor, $\mathbf{A}=\begin{array}{cccc}{\mathbf{A}}_1& {\mathbf{A}}_2& \cdots & {\mathbf{A}}_n\end{array}\in {ℝ}^{6\times 4n}$ and $\mathbf{u}={{\mathbf{u}}_1^T,{\mathbf{u}}_2^T,\cdots ,{\mathbf{u}}_n^T}^T\in {ℝ}^{4n}$. Note that ${\mathbf{z}}_{bi}={0,0,1}^T$ and ${\mathbf{u}}_i$ is described in the body frame of the $i$th quadrotor.

In [22, 23], it is assumed that all quadrotors are in parallel planes for the transportation system consisting of $n$ quadrotors carrying a payload rigidly. Consequently, the dynamics equations are the same as those of the classical quadrotor. The only difference from the classical quadrotor is that more propellers are included to generate the desired control command. It should be noted that the dynamics equations in [22, 23] is still underactuated; hence, the motion in $X$ and $Y$ directions is realized based on the attitude maneuver in roll and pitch directions. However, the quadrotors in the cooperative transportation system considered in this paper can be in any feasible configuration, i.e., all quadrotors may not be in a parallel plane. Consequently, the dynamics equation in Equation (2) is different from the classical quadrotor. The accelerations in $X$ and $Y$ directions may be generated by changing thrust forces of some quadrotors, which implies that maneuvers in $X$ and $Y$ directions are more agile than the classical quadrotor-like system.

The three-dimensional force in the frame $xyz$ provided by the $i$th quadrotor can be expressed as ${\mathbf{R}}_i{\mathbf{z}}_{bi}{f}_i$. Hence, the total force vector from these $n$ quadrotors expressed in the inertial frame is $\mathbf{R}·{\sum }_{i=1}^n{\mathbf{R}}_i{\mathbf{z}}_{bi}{f}_i$. Essentially, $\mathbf{R}·{\mathbf{R}}_i$, the rotation matrix from the body frame of the $i$th quadrotor to the inertial frame, can be written as follows: $$\begin{array}{}\text{(7)}& \mathbf{R}·{\mathbf{R}}_i=\begin{array}{ccc}{c}_{{\overline{\psi }}_i}{c}_{{\overline{\theta }}_i}& -s_{{\overline{\psi }}_i}{c}_{{\overline{\varphi }}_i}+c_{{\overline{\psi }}_i}{s}_{{\overline{\theta }}_i}{s}_{{\overline{\varphi }}_i}& s_{{\overline{\psi }}_i}{s}_{{\overline{\varphi }}_i}+c_{{\overline{\psi }}_i}{s}_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}\\ s_{{\overline{\psi }}_i}{c}_{{\overline{\theta }}_i}& c_{{\overline{\psi }}_i}{c}_{{\overline{\varphi }}_i}+s_{{\overline{\psi }}_i}{s}_{{\overline{\theta }}_i}{s}_{{\overline{\varphi }}_i}& -c_{{\overline{\psi }}_i}{s}_{{\overline{\varphi }}_i}+s_{{\overline{\psi }}_i}{s}_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}\\ -s_{{\overline{\theta }}_i}& c_{{\overline{\theta }}_i}{s}_{{\overline{\varphi }}_i}& c_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}\end{array},\end{array}$$where $s_p$ and $c_p$ represent $\sin p$ and $\cos p$ with $p={\overline{\varphi }}_i,{\overline{\theta }}_i,{\overline{\psi }}_i$. Note that ${\overline{\varphi }}_i$, ${\overline{\theta }}_i$, and ${\overline{\psi }}_i$ are the roll, pitch, and yaw angles associated with the rotation described by $\mathbf{R}·{\mathbf{R}}_i$.

Hence, the lift provided by these $n$ quadrotors can be written as ${\sum }_{i=1}^n{c}_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}{f}_i$. The forces acting on the system along $X$ and $Y$ axes are ${\sum }_{i=1}^n{s}_{{\overline{\psi }}_i}{s}_{{\overline{\varphi }}_i}+c_{{\overline{\psi }}_i}{s}_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}{f}_i$ and ${\sum }_{i=1}^n-c_{{\overline{\psi }}_i}{s}_{{\overline{\varphi }}_i}+s_{{\overline{\psi }}_i}{s}_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}{f}_i$. One easy way to check whether the whole system can take off successfully is to ensure that ${\sum }_{i=1}^n{c}_{{\overline{\theta }}_i}{c}_{{\overline{\varphi }}_i}{f}_i^{\max }$ is larger than the system gravity with the roll angle ${\overline{\varphi }}_i$ and pitch angle ${\overline{\theta }}_i$ in the initial configuration and the maximum thrust $f_i^{\max }$ of each quadrotor. Another important conclusion is that ${\overline{\varphi }}_i$ and ${\overline{\theta }}_i$ closer to zero will provide a larger lift force. However, in such a case, the forces acting on the system along $X$ and $Y$ axes will tend to zero.

3. Control System Design

This section is aimed at designing a controller with the following assumptions for the system governed by Equations (1)–(4) and presents a command allocation method to realize the desired control command using $n$ quadrotors.

Assumption 1.

The time derivative of the possible unmodeled terms ${\mathbf{d}}_1$ and ${\mathbf{d}}_2$ is bounded.

The unmodeled terms ${\mathbf{d}}_1$ and ${\mathbf{d}}_2$ are mainly from the uncertainty of the actuator model and the modeling error from the measuring errors of inertial parameters and relative pose between quadrotors. For a practical quadrotor transportation system, these modeling errors cannot change with an infinite speed. Hence, it is reasonable to make Assumption 1.

Assumption 2.

The control of the payload’s position and yaw angle is more significant in the transportation system. Suppose that the desired position vector and yaw angle of the entire system are ${\mathbf{r}}_d=x_d,y_d,z_d$ and ${\psi }_d$.

3.1. Position Control

Rewrite dynamics Equation (2) as $$\begin{array}{}\text{(8)}& m{\dot{\mathbf{v}}}_c=-mg{\mathbf{z}}_I+{\mathbf{f}}_d+{\mathbf{d}}_1,\end{array}$$where ${\mathbf{f}}_d=\mathbf{R}\mathbf{f}$ represents the reference control force command in the inertial frame, i.e., the control command for the system in Equation (2) can be written as $\mathbf{f}={\mathbf{R}}^T{\mathbf{f}}_d$. Consider the following linear extended state observer to estimate unmodeled terms ${\mathbf{d}}_1$: $$\begin{array}{c}\text{(9)}& {\dot{\widehat{\mathbf{q}}}}_1={\widehat{\mathbf{q}}}_2+3w_0{\mathbf{r}}_c-{\widehat{\mathbf{q}}}_1,{\dot{\widehat{\mathbf{q}}}}_2={\widehat{\mathbf{q}}}_3+3w_0^2{\mathbf{r}}_c-{\widehat{\mathbf{q}}}_1+\frac{{\mathbf{f}}_d}{m}-g{\mathbf{z}}_I,\\ {\dot{\widehat{\mathbf{q}}}}_3=w_0^3{\mathbf{r}}_c-{\widehat{\mathbf{q}}}_1,\end{array}$$where $w_0$ is a positive observer gain. Letting ${\tilde{\mathbf{q}}}_1={\mathbf{r}}_c-{\widehat{\mathbf{q}}}_1$, ${\tilde{\mathbf{q}}}_2={\dot{\mathbf{r}}}_c-{\widehat{\mathbf{q}}}_2$, and ${\tilde{\mathbf{q}}}_3={\mathbf{d}}_1/m-{\widehat{\mathbf{q}}}_3$ yields $$\begin{array}{c}\text{(10)}& {\dot{\stackrel{~}{\mathbf{q}}}}_1={\tilde{\mathbf{q}}}_2-3w_0{\tilde{\mathbf{q}}}_1,{\dot{\stackrel{~}{\mathbf{q}}}}_2={\tilde{\mathbf{q}}}_3-3w_0^2{\tilde{\mathbf{q}}}_1,\\ {\dot{\stackrel{~}{\mathbf{q}}}}_3=-w_0^3{\tilde{\mathbf{q}}}_1+\frac{{\dot{\mathbf{d}}}_1}{m}.\end{array}$$

The above equation can be rewritten as $$\begin{array}{}\text{(11)}& \dot{\stackrel{~}{\mathbf{q}}}=\overline{\mathbf{A}}\tilde{\mathbf{q}}+\overline{\mathbf{d}},\end{array}$$where $$\begin{array}{c}\text{(12)}& \tilde{\mathbf{q}}={{\tilde{\mathbf{q}}}_1^T,{\tilde{\mathbf{q}}}_2^T,{\tilde{\mathbf{q}}}_3^T}^T,\overline{\mathbf{d}}={\mathbf{0}}^T,{\mathbf{0}}^T,{\dot{\mathbf{d}}}_1^T/m,\\ \overline{\mathbf{A}}=\begin{array}{ccc}-3w_0{\mathbf{I}}_3& {\mathbf{I}}_3& \mathbf{0}\\ -3w_0^2{\mathbf{I}}_3& \mathbf{0}& {\mathbf{I}}_3\\ -w_0^3{\mathbf{I}}_3& \mathbf{0}& \mathbf{0}\end{array}.\end{array}$$

It is easy to verify all eigenvalues of $\mathbf{A}$ are negative with a positive constant $w_0$, which implies that $\dot{\stackrel{~}{\mathbf{q}}}=\overline{\mathbf{A}}\tilde{\mathbf{q}}$ is asymptotically stable. Therefore, with a finite $\overline{\mathbf{d}}$, $\tilde{\mathbf{q}}$ is bounded.

According to [25, 26], one has the following lemma.

Lemma 3 (see [25, 26]).

With Assumption 1 and ${\text{w}}_0>0$, ${\tilde{\mathbf{q}}}_1$, ${\tilde{\mathbf{q}}}_2$, and ${\tilde{\mathbf{q}}}_3$ are bounded and there are a constant ${\sigma }_{i,j}>0$ and a finite time $T_1$ such that ${\tilde{q}}_{i,j}\le {\sigma }_{i,j}$ for any $t>T_1$, where ${\tilde{q}}_{i,j}$ represents the $j$th element of ${\tilde{\mathbf{q}}}_i$ for $i=1,2,3$ and $j=1,2,3$ and ${\sigma }_{i,j}=O1/w_0^k$ for some positive integer $k$.

Essentially, ${\widehat{\mathbf{q}}}_1$, ${\widehat{\mathbf{q}}}_2$, and ${\widehat{\mathbf{q}}}_3$ are the estimates of ${\mathbf{r}}_c$, ${\dot{\mathbf{r}}}_c$, and ${\mathbf{d}}_1/m$. Based on the extended state observer, the position controller can be designed as $$\begin{array}{}\text{(13)}& {\mathbf{f}}_d=mg{\mathbf{z}}_I-{\mathbf{K}}_p+{\mathbf{K}}_i{\mathbf{r}}_c-{\mathbf{r}}_d-{\mathbf{K}}_i{\int }_0^t{\mathbf{r}}_c-{\mathbf{r}}_{d}dt-{\mathbf{K}}_d{\mathbf{v}}_c-{\dot{\mathbf{r}}}_d+m{\ddot{\mathbf{r}}}_d-m{\widehat{\mathbf{q}}}_3,\end{array}$$where ${\mathbf{K}}_p$, ${\mathbf{K}}_i$, and ${\mathbf{K}}_d$ are three diagonal positive definite matrices and ${\int }_0^t\mathbf{a}dt$ is defined as ${{\int }_0^t{a}_{1}dt,{\int }_0^t{a}_{2}dt,{\int }_0^t{a}_{3}dt}^T$ for a vector $\mathbf{a}={a_1,a_2,a_3}^T$.

Theorem 4.

Under the controller ${\mathbf{f}}_{\text{d}}$ in Equation (13), if these diagonal positive definite matrices ${\mathbf{K}}_{\text{p}}$, ${\mathbf{K}}_{\text{i}}$, and ${\mathbf{K}}_{\text{d}}$ are chosen such that all eigenvalues of $$\begin{array}{}\text{(14)}& \begin{array}{ccc}\mathbf{0}& {\mathbf{I}}_3& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{I}}_3\\ -\frac{{\mathbf{K}}_{\text{i}}}{m}& -\frac{{\mathbf{K}}_{\text{p}}+{\mathbf{K}}_i}{m}& -\frac{{\mathbf{K}}_{\text{d}}}{m}\end{array}\end{array}$$have negative real parts, the position tracking error is bounded.

Proof.

See Appendix B.

Remark 5.

As shown in the proof of Theorem 4, a large ${\omega }_0$ will result in a smaller tracking error. However, a too large ${\omega }_0$ will make the observer sensitive to noise [27]. Hence, ${\omega }_0$ should be chosen to balance the tracking performance and the noise tolerance.

Remark 6.

Theorem 4 presents the stability of the designed position controller, but beyond that, different control gains may result in different system behaviors. This remark is aimed at presenting some tips on how to choose the gains in the designed controller (13). As shown in Equation (B.1), the designed controller can be considered a classical PID controller for a double-integrator system with ${\mathbf{K}}_{\text{p}}+{\mathbf{K}}_{\text{i}}$, ${\mathbf{K}}_{\text{d}}$, and ${\mathbf{K}}_{\text{i}}$ as the proportional, derivative, and integral gains. Hence, the control gains ${\mathbf{K}}_{\text{p}}+{\mathbf{K}}_{\text{i}}$, ${\mathbf{K}}_{\text{d}}$, and ${\mathbf{K}}_{\text{i}}$ can be adjusted according to the PID controller tuning methods.

3.2. Attitude Control

In the last section, to realize the desired control force ${\mathbf{f}}_{\text{d}}$, the control force $\mathbf{f}$ can be simply chosen as $\mathbf{f}={\mathbf{R}}^T{\mathbf{f}}_{\text{d}}$. On the one hand, according to Assumption 2, the attitude maneuver, especially in roll and pitch directions, in aerial payload transportation is usually not significant. On the other hand, to improve the maneuverability in $X$ and $Y$ directions, the reference rotation matrix can be set as $$\begin{array}{}\text{(15)}& {\mathbf{R}}_d=\begin{array}{ccc}{c}_{{\psi }_d}{c}_{{\theta }_d}& -s_{{\psi }_d}{c}_{{\varphi }_d}+c_{{\psi }_d}{s}_{{\theta }_d}{s}_{{\varphi }_d}& s_{{\psi }_d}{s}_{{\varphi }_d}+c_{{\psi }_d}{s}_{{\theta }_d}{c}_{{\varphi }_d}\\ s_{{\psi }_d}{c}_{{\theta }_d}& c_{{\psi }_d}{c}_{{\varphi }_d}+s_{{\psi }_d}{s}_{{\theta }_d}{s}_{{\varphi }_d}& -c_{{\psi }_d}{s}_{{\varphi }_d}+s_{{\psi }_d}{s}_{{\theta }_d}{c}_{{\varphi }_d}\\ -s_{{\theta }_d}& c_{{\theta }_d}{s}_{{\varphi }_d}& c_{{\theta }_d}{c}_{{\varphi }_d}\end{array},\end{array}$$where ${\varphi }_d$, ${\theta }_d$, and ${\psi }_d$ are the desired roll, pitch, and yaw angles satisfying $s_{{\varphi }_d}=f_{d,1}{s}_{{\psi }_d}-f_{d,2}{c}_{{\psi }_d}/\parallel {\mathbf{f}}_d\parallel$ and $\tan {\theta }_d=f_{d,1}{c}_{{\psi }_d}+f_{d,2}{s}_{{\psi }_d}/f_{d,3}$. Note that $f_{d,i}$ represents the $i$th element of ${\mathbf{f}}_d$. For such a transportation system, it is reasonable to assume that both roll and pitch angles belong to $-\pi /2,\pi /2$. Hence, ${\varphi }_d$ and ${\theta }_d$ can be set as ${\varphi }_d=\arcsin f_{d,1}{s}_{{\psi }_d}-f_{d,2}{c}_{{\psi }_d}/\parallel {\mathbf{f}}_d\parallel$ and ${\theta }_d=\arctan f_{d,1}{c}_{{\psi }_d}+f_{d,2}{s}_{{\psi }_d}/f_{d,3}$, respectively. It should be noted that ${\psi }_d$ is the desired yaw angle defined according to the mission requirement. The method to define the desired attitude is the same as the one for a classical quadrotor. Hence, based on such a desired attitude, the proposed control method is compatible with the transportation mission with partial or all quadrotors in parallel planes.

To estimate the unmodeled attitude dynamics, the following observer is designed as $$\begin{array}{c}\text{(16)}& \dot{\widehat{\mathbf{\Omega }}}=2\Gamma \mathbf{\Omega }-\widehat{\mathbf{\Omega }}-{\mathbf{J}}^{-1}\mathbf{\Omega }\times \mathbf{J}\mathbf{\Omega }+{\mathbf{J}}^{-1}\tau +{\widehat{\mathbf{d}}}_2,{\dot{\widehat{\mathbf{d}}}}_2={\Gamma }^2\mathbf{\Omega }-\widehat{\mathbf{\Omega }},\end{array}$$where $\Gamma$ is a positive constant. Letting ${\mathbf{e}}_{\Omega }=\mathbf{\Omega }-\widehat{\mathbf{\Omega }}$ and ${\mathbf{e}}_{d2}={\mathbf{J}}^{-1}{\mathbf{d}}_2-{\widehat{\mathbf{d}}}_2$ yields $$\begin{array}{}\text{(17)}& {\dot{\mathbf{e}}}_{\Omega }=-2\Gamma {\mathbf{e}}_{\Omega }+{\mathbf{e}}_{d2},\text{(18)}& {\dot{\mathbf{e}}}_{d2}={\mathbf{J}}^{-1}{\dot{\mathbf{d}}}_2-{\Gamma }^2{\mathbf{e}}_{\Omega }.\end{array}$$

According to Theorem 2.2 in [28], one can have the following lemma.

Lemma 7.

With Assumption 1 and $\Gamma >0$, ${\mathbf{e}}_{\Omega }$ and ${\mathbf{e}}_{d2}$ are convergent in the sense that for any ${\sigma }_{\omega }\in 0,1$, there exists ${\epsilon }_{\sigma }\in$(0, 1) such that $e_{\Omega ,i}<{\sigma }_{\omega }$ and $e_{d2,i}<{\sigma }_{\omega }$ hold for any ${\epsilon }_1\in 0,{\epsilon }_{\sigma }$ and any $t\in T_{\epsilon },\infty$, where $T_{\epsilon }>0$ depends on ${\epsilon }_1$ and $e_{\Omega ,i}$, where $e_{\Omega ,i}$ and $e_{d2,i}$ are the $i$th elements of ${\mathbf{e}}_{\Omega }$ and ${\mathbf{e}}_{d2}$, respectively.

To track the desired attitude ${\mathbf{R}}_d$, a robust controller is designed as follows: $$\begin{array}{}\text{(19)}& \tau ={\mathbf{\Omega }}^{\times }\mathbf{J}\mathbf{\Omega }-k_R{\mathbf{e}}_R-k_{\Omega }{\mathbf{E}}_{\Omega }+\mathbf{J}{\dot{\mathbf{\Omega }}}_d-\mathbf{J}{\widehat{\mathbf{d}}}_2,\end{array}$$

where ${\mathbf{e}}_R={{\mathbf{R}}_d^T\mathbf{R}-{\mathbf{R}}^T{\mathbf{R}}_d}^{\vee }$ and ${\mathbf{E}}_{\Omega }=\mathbf{\Omega }-{\mathbf{\Omega }}_d$. Note that ${\dot{\mathbf{e}}}_R={-{\mathbf{\Omega }}_d^{\times }{\mathbf{R}}_d^T\mathbf{R}+{\mathbf{R}}_d^T\mathbf{R}{\mathbf{\Omega }}^{\times }+{\mathbf{\Omega }}^{\times }{\mathbf{R}}^T{\mathbf{R}}_d-{\mathbf{R}}^T{\mathbf{R}}_d{\mathbf{\Omega }}_d^{\times }}^{\vee }=\text{trace}{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{\Omega }-\text{trace}{\mathbf{R}}_d^T\mathbf{R}\mathbf{I}-{\mathbf{R}}_d^T\mathbf{R}{\mathbf{\Omega }}_d=\text{trace}{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d{\mathbf{E}}_{\Omega }+{\mathbf{e}}_R^{\times }{\mathbf{E}}_{\Omega }$. ${\widehat{\mathbf{d}}}_2$ is the estimate of the unknown term ${\mathbf{J}}^{-1}{\mathbf{d}}_2$ and updated according to the law in Equation (17).

Theorem 8.

If $0<\mu <\sqrt{k_R{\lambda }_{\min }\mathbf{J}/2}$ holds and $$\begin{array}{}\text{(20)}& \begin{array}{cc}{k}_R\mu {\lambda }_{\min }{\mathbf{J}}^{-1}& -\frac{k_{\Omega }\mu {\lambda }_{\max }{\mathbf{J}}^{-1}}{2}\\ -\frac{k_{\Omega }\mu {\lambda }_{\max }{\mathbf{J}}^{-1}}{2}& k_{\Omega }-2\mu \end{array}\end{array}$$

is positive definite, the attitude tracking error is bounded with the designed controller exponentially almost globally.

Proof.

See Appendix C.

Remark 9.

As shown in the proof of Theorem 8, the bounded stability of the system governed by Equation (C.2) is almost globally exponential, i.e., there may exist some undesired equilibria satisfying $\cos \Theta =-1$. According to [29], these undesired critical points are unstable. In a practical system, the ubiquitous noise can drive the system to leave these unstable equilibria.

3.3. Control Command Allocation

Suppose that the desired six-dimensional control command for the rigid body governed by Equations (1)–(4) is ${\mathbf{f},{\tau }^T}^T$. The rotation matrix ${\mathbf{R}}_i$ from the $i$th quadrotor’s body frame to the body frame of the entire system can be written as $$\begin{array}{}\text{(21)}& {\mathbf{R}}_i=\begin{array}{ccc}{c}_{{\psi }_i}{c}_{{\theta }_i}& -s_{{\psi }_i}{c}_{{\varphi }_i}+c_{{\psi }_i}{s}_{{\theta }_i}{s}_{{\varphi }_i}& s_{{\psi }_i}{s}_{{\varphi }_i}+c_{{\psi }_i}{s}_{{\theta }_i}{c}_{{\varphi }_i}\\ s_{{\psi }_i}{c}_{{\theta }_i}& c_{{\psi }_i}{c}_{{\varphi }_i}+s_{{\psi }_i}{s}_{{\theta }_i}{s}_{{\varphi }_i}& -c_{{\psi }_i}{s}_{{\varphi }_i}+s_{{\psi }_i}{s}_{{\theta }_i}{c}_{{\varphi }_i}\\ -s_{{\theta }_i}& c_{{\theta }_i}{s}_{{\varphi }_i}& c_{{\theta }_i}{c}_{{\varphi }_i}\end{array},\end{array}$$where ${\varphi }_i$, ${\theta }_i$, and ${\psi }_i$ are the roll, pitch, and yaw angles associated with the rotation ${\mathbf{R}}_i$.

To determine the control inputs for each quadrotor, one can minimize the following cost function $$\begin{array}{}\text{(22)}& J_c={\mathbf{u}}^T\mathbf{W}\mathbf{u}\end{array}$$with the following equation constraint $$\begin{array}{}\text{(23)}& \mathbf{A}\mathbf{u}={\mathbf{f},{\tau }^T}^T,\end{array}$$where $\mathbf{W}=\mathrm{diag}{\mathbf{w}}_i\in {ℝ}^{4n\times 4n}$ is the weight matrix, in which ${\mathbf{w}}_i=\mathrm{diag}{w}_{i,j}$ with $w_{i,j}>0$ for $i=1,2,\cdots ,n$ and $j=1,2,3,4$. Essentially, this is a constrained quadratic problem, whose solution is $$\begin{array}{}\text{(24)}& \mathbf{u}={\mathbf{K}}^{+}{f,{\tau }^T}^T,\end{array}$$where ${\mathbf{K}}^{+}={\mathbf{W}}^{-1}{\mathbf{A}}^T{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T}^{-1}$. Note that $\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T$ can be written as ${\sum }_{i=1}^n{\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T$. Some special cases are discussed below.

Case 1.

The weight matrix ${\mathbf{w}}_i$ is chosen the same for each quadrotor, and the weights in roll, pitch, and yaw directions are equal. Hence, ${\mathbf{w}}_i$ can be rewritten as $\mathrm{diag}{w}_1,w_2,w_2,w_2$. Since each quadrotor’s payload capability is limited, to improve the payload capability of the quadrotor group in a transportation mission, it is desirable that both ${\varphi }_i$ and ${\theta }_i$ are small enough such that $\cos {\varphi }_i\approx 1$, $\cos {\theta }_i\approx 1$, $\sin {\varphi }_i\approx {\varphi }_i$, $\sin {\theta }_i\approx {\theta }_i$, ${\varphi }_i^2\approx 0$, ${\theta }_i^2\approx 0$, and ${\varphi }_i{\theta }_i\approx 0$ hold. Hence, ${\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T$ can be written as $$\begin{array}{}\text{(25)}& {\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T=\begin{array}{ccc}{\mathbf{0}}_{2\times 2}& {\mathbf{h}}_1& {\mathbf{0}}_{2\times 1}\\ {\mathbf{h}}_1^T& {\mathbf{h}}_2& {\mathbf{h}}_3\\ {\mathbf{0}}_{1\times 2}& {\mathbf{h}}_3^T& h_4\end{array},\end{array}$$where $$\begin{array}{c}\text{(26)}& {\mathbf{h}}_1=\begin{array}{ccc}\frac{c_{{\psi }_i}{\theta }_i+s_{{\psi }_i}{\varphi }_i}{w_1}& \frac{y_i{c}_{{\psi }_i}{\theta }_i+s_{{\psi }_i}{\varphi }_i}{w_1}& -\frac{x_i{c}_{{\psi }_i}{\theta }_i+s_{{\psi }_i}{\varphi }_i}{w_1}\\ -\frac{c_{{\psi }_i}{\varphi }_i-s_{{\psi }_i}{\theta }_i}{w_1}& -\frac{y_i{c}_{{\psi }_i}{\varphi }_i-s_{{\psi }_i}{\theta }_i}{w_1}& \frac{x_i{c}_{{\psi }_i}{\varphi }_i-s_{{\psi }_i}{\theta }_i}{w_1}\end{array},{\mathbf{h}}_2=\begin{array}{ccc}{w_1}^{-1}& \frac{c_{{\psi }_i}{\varphi }_i{z}_i-s_{{\psi }_i}{\theta }_i{z}_i+y_i}{w_1}& \frac{c_{{\psi }_i}{\theta }_i{z}_i+s_{{\psi }_i}{\varphi }_i{z}_i-x_i}{w_1}\\ \frac{c_{{\psi }_i}{\varphi }_i{z}_i-s_{{\psi }_i}{\theta }_i{z}_i+y_i}{w_1}& \frac{2c_{{\psi }_i}{\varphi }_i{w}_2{y}_i{z}_i-2s_{{\psi }_i}{\theta }_i{w}_2{y}_i{z}_i+w_2{y_i}^2+w_1}{w_1{w}_2}& -\frac{c_{{\psi }_i}{\varphi }_i{x}_i{z}_i-c_{{\psi }_i}{\theta }_i{y}_i{z}_i-s_{{\psi }_i}{\varphi }_i{y}_i{z}_i-s_{{\psi }_i}{\theta }_i{x}_i{z}_i+x_i{y}_i}{w_1}\\ \frac{c_{{\psi }_i}{\theta }_i{z}_i+s_{{\psi }_i}{\varphi }_i{z}_i-x_i}{w_1}& -\frac{c_{{\psi }_i}{\varphi }_i{x}_i{z}_i-c_{{\psi }_i}{\theta }_i{y}_i{z}_i-s_{{\psi }_i}{\varphi }_i{y}_i{z}_i-s_{{\psi }_i}{\theta }_i{x}_i{z}_i+x_i{y}_i}{w_1}& -\frac{2c_{{\psi }_i}{\theta }_i{w}_2{x}_i{z}_i+2s_{{\psi }_i}{\varphi }_i{w}_2{x}_i{z}_i-w_2{x_i}^2-w_1}{w_1{w}_2}\end{array},\\ {\mathbf{h}}_3=\begin{array}{c}-\frac{c_{{\psi }_i}{\varphi }_i{x}_i+c_{{\psi }_i}{\theta }_i{y}_i+s_{{\psi }_i}{\varphi }_i{y}_i-s_{{\psi }_i}{\theta }_i{x}_i}{w_1}\\ -\frac{y_i{c}_{{\psi }_i}{\varphi }_i{x}_i+c_{{\psi }_i}{\theta }_i{y}_i+s_{{\psi }_i}{\varphi }_i{y}_i-s_{{\psi }_i}{\theta }_i{x}_i}{w_1}\\ \frac{x_i{c}_{{\psi }_i}{\varphi }_i{x}_i+c_{{\psi }_i}{\theta }_i{y}_i+s_{{\psi }_i}{\varphi }_i{y}_i-s_{{\psi }_i}{\theta }_i{x}_i}{w_1}\end{array},\\ h_4={w_2}^{-1}.\end{array}$$

It should be noted that the determinant of ${\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T$ is equal to zero, which implies that only one quadrotor cannot generate the six-dimensional control command based on the equipped four propellers. If there are four quadrotors in the team with ${\varphi }_1=0.1$ rad, ${\varphi }_2=-0.05$ rad, ${\varphi }_3={\varphi }_4=0$, ${\theta }_1={\theta }_2=0$, ${\theta }_3=0.1$ rad, ${\theta }_4=-0.05$ rad, ${\psi }_1={\psi }_2={\psi }_3={\psi }_4=0$, $x_1=-x_2=1$, $x_3=x_4=0$, $y_1=y_2=0$, $y_3=-y_4=1$, $z_1=z_2=z_3=z_4=0$, and $w_1=w_2=1$, the determinant of ${\sum }_{i=1}^4{\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T$ is 0.01, which implies that these four quadrotors can provide both three-dimensional forces and three-dimensional torques. However, if letting ${\theta }_3={\theta }_4=0$ and keeping all remaining parameters the same as those in the previous sentence, ${\sum }_{i=1}^4{\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T$ will be irreversible, i.e., the six-dimensional control command cannot be provided if two quadrotors are in a parallel plane for such a four-quadrotor team. Hence, to improve the compatibility of Equation (24), a new control command allocation method is defined as follows: $$\begin{array}{}\text{(27)}& \mathbf{u}=\begin{array}{cc}{\mathbf{W}}^{-1}{\mathbf{A}}^T{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T}^{-1}{{\mathbf{f}}^T,{\tau }^T}^T,\hfill & \text{if}\text{Det}\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T>\epsilon ,\hfill \\ {\mathbf{W}}^{-1}{\mathbf{A}}^T{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}{{\mathbf{f}}^T,{\tau }^T}^T,\hfill & \text{otherwise},\hfill \end{array}\end{array}$$where $\epsilon$ and $\delta$ are small positive constants.

Case 2.

If the $x_i{y}_i$ plane in the body frame $x_i{y}_i{z}_i$ is parallel to $xy$ plane in the frame $xyz$, ${\varphi }_i=0$ and ${\theta }_i=0$ hold for each quadrotor and $$\begin{array}{}\text{(28)}& {\mathbf{A}}_i=\begin{array}{c}{\mathbf{0}}_{2\times 4}\\ {\widehat{\mathbf{A}}}_i\end{array},\end{array}$$where ${\widehat{\mathbf{A}}}_i$ is $$\begin{array}{}\text{(29)}& {\widehat{\mathbf{A}}}_i=\begin{array}{cccc}1& 0& 0& 0\\ y_i& c_{{\psi }_i}& -s_{{\psi }_i}& 0\\ -x_i& s_{{\psi }_i}& c_{{\psi }_i}& 0\\ 0& 0& 0& 1\end{array}.\end{array}$$

Hence, ${\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T$ is $$\begin{array}{}\text{(30)}& {\mathbf{A}}_i{\mathbf{w}}_i^{-1}{\mathbf{A}}_i^T=\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& {w_1}^{-1}& \frac{y_i}{w_1}& -\frac{x_i}{w_1}& 0\\ 0& 0& \frac{y_i}{w_1}& \frac{w_2{y_i}^2+w_1}{w_1{w}_2}& -\frac{x_i{y}_i}{w_1}& 0\\ 0& 0& -\frac{x_i}{w_1}& -\frac{x_i{y}_i}{w_1}& \frac{w_2{x_i}^2+w_1}{w_1{w}_2}& 0\\ 0& 0& 0& 0& 0& {w_2}^{-1}\end{array}.\end{array}$$

The nonzero block matrix in the bottom right corner agrees with the matrix format in [22], where the aerial transportation using quadrotors in parallel planes is studied. In fact, the case in [22] is a special case of the general transportation configuration in Figure 2. If the assumption that the quadrotors are much heavier than the payload in [22] holds, $\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T$ would be $$\begin{array}{}\text{(31)}& \mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T=\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{n}{w_1}& 0& 0& 0\\ 0& 0& 0& \sum _{i=1}^n\frac{w_2{y_i}^2+w_1}{w_1{w}_2}& 0& 0\\ 0& 0& 0& 0& \sum _{i=1}^n\frac{w_2{x_i}^2+w_1}{w_1{w}_2}& 0\\ 0& 0& 0& 0& 0& \frac{n}{w_2}\end{array}.\end{array}$$

It is clear that the matrix in the above equation is singular. Based on Equation (27), ${\mathbf{A}}^T{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}$ can be written as $$\begin{array}{}\text{(32)}& \begin{array}{c}{\mathbf{A}}_1{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}\\ \cdots \\ {\mathbf{A}}_n{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}\end{array},\end{array}$$where $$\begin{array}{}\text{(33)}& {\mathbf{A}}_i{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}=\begin{array}{cc}{\mathbf{0}}_{4\times 2}& {\widehat{\mathbf{A}}}_i^T\end{array}{\mathrm{diag}0,0,\frac{n}{w_1},\sum _{i=1}^n\frac{w_2{y_i}^2+w_1}{w_1{w}_2},\sum _{i=1}^n\frac{w_2{x_i}^2+w_1}{w_1{w}_2},\frac{n}{w_2}+\delta {\mathbf{I}}_6}^{-1}\approx \begin{array}{cc}{\mathbf{0}}_{4\times 2}& {\widehat{\mathbf{A}}}_i^T{\mathrm{diag}\frac{n}{w_1},\sum _{i=1}^n\frac{w_2{y_i}^2+w_1}{w_1{w}_2},\sum _{i=1}^n\frac{w_2{x_i}^2+w_1}{w_1{w}_2},\frac{n}{w_2}}^{-1}\end{array}.\end{array}$$

Let $f_3$ be the third element of $\mathbf{f}$. ${\mathbf{A}}_i{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}{{\mathbf{f}}^T,{\tau }^T}^T$ can be written as $$\begin{array}{}\text{(34)}& {\widehat{\mathbf{A}}}_i^T{\mathrm{diag}\frac{n}{w_1},\sum _{i=1}^n\frac{w_2{y_i}^2+w_1}{w_1{w}_2},\sum _{i=1}^n\frac{w_2{x_i}^2+w_1}{w_1{w}_2},\frac{n}{w_2}}^{-1}{f_3,{\tau }^T}^T,\end{array}$$which agrees with the result from the method in [22].

Another interesting point is that both ${\mathbf{W}}^{-1}{\mathbf{A}}^T{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T}^{-1}$ and ${\mathbf{W}}^{-1}{\mathbf{A}}^T{\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T+\delta {\mathbf{I}}_6}^{-1}$ are only determined by the relative pose of quadrotors with respect to the payload. Hence, it can be calculated in advance in a transportation task. In other words, the control command for each quadrotor is obtained just from the product of a constant matrix and a six-dimensional control input vector.

In conclusion, the method to calculate the control command for each quadrotor presented in this study is illustrated by Figure 3. With the position vector of the mass center of the entire system and the desired transportation destination, Equation (13) will provide the value of ${\mathbf{f}}_d$, based on which the reference attitude of the entire system and the control force in Equation (2) can be obtained. With the attitude command from Equation (19), the control command allocation algorithm in Equation (27) provides a 4-dimensional control command for each quadrotor. It should be noted that one quadrotor can be set as a leader to run the controller and control command allocation algorithm and then send other control commands to the remaining quadrotors (followers).

[figure(s) omitted; refer to PDF]

4. Experimental Verification

As shown in Figure 4, the indoor test rig mainly consists of the OptiTrack motion capture system with 16 Flex 13 cameras, a work station, a router, two Quanser® QDrone quadrotors, and a payload to be transported. The QDrone quadrotor is a midsize quadrotor equipped with a powerful onboard Intel® Aero Compute Board, multiple high-resolution cameras, a BMI160 IMU, and the built-in WiFi. The motion capture system can measure the position and attitude of an object at 100 Hz. The real-time control of quadrotors is realized using Simulink from MathWorks and QUARC software [30] from Quanser Consulting Inc.

[figure(s) omitted; refer to PDF]

Magnets are used to connect the quadrotors and the board. As shown in the subfigure at the bottom of Figure 4, the right attachment point of quadrotor 1 has more magnets. Hence, the propeller plane of the first quadrotor is not parallel to the top surface of the board. Similarly, for quadrotor 2, the left connection point has more magnets than the right one. Mathematically, quadrotors 1 and 2 rotate $-0.102$ rad and 0.107 rad along the $x$ axis in the frame of the entire system, respectively. Hence, the thrust forces of both quadrotors contribute to the movement in $y$ direction. However, the movement along $x$ axis can only be realized by the attitude maneuver in the pitch direction. In the body frame of the entire system, the position vectors of these two quadrotors are ${\mathbf{r}}_1={\mathrm{0.022,0.56},-0.001}^T$ m and ${\mathbf{r}}_2={-0.022,-\mathrm{0.56,0.001}}^T$ m, respectively. Therefore, the system in Figure 4 can well verify the effectiveness of the proposed controller and control command allocation method since the movement in $x$ direction can be considered a special case with ${\theta }_i=0$ for $i=1,2$.

The detailed physical parameters of quadrotor and payload are given in Table 1. Since the payload mass is larger than the maximum payload capability of a single quadrotor, at least two quadrotors have to be used to transport such a payload. In tests, the mass center of the entire system is defined as the middle point of two quadrotors approximately because these two quadrotors are much heavier than the payload. Therefore, the workstation will send the measured and desired position and yaw angles of the entire system to the first quadrotor via WiFi at 100 Hz, and then, its onboard compute board will calculate the control command for both quadrotors based on the designed controllers (13) and (19) with roll and pitch angles measured by IMU and the control command allocation law (27) at a rate of 1 kHz. The last four elements of $\mathbf{u}$ from Equation (27) will be sent to quadrotor 2 using WiFi, which implies that the onboard compute board of the second quadrotor will just implement the received four-dimensional control command without calculating the control command by itself. In essence, a centralized leader-follower architecture is adopted here. Note that some markers are attached to the payload to record the payload position only. The parameters in position and attitude controllers and control command allocation algorithm are given in Table 2.

Table 1

System parameters.

Parameter	Value
Quadrotor mass	1.121 kg
Inertia matrix of each quadrotor	$\mathrm{diag}0.01,0.0082,0.0148$ kg·m2
Maximum payload capability of each quadrotor	$~$0.3 kg
Payload mass	0.453 kg
Payload dimensions in $x$, $y$, and $z$ directions	$$0.038\text{m}\times 1.219\text{m}\times 0.013\text{m}$$
Inertia matrix of payload	$\mathrm{diag}0.0561,0.00006,0.0561$ kg·m2
Constant gravitational acceleration	9.81 m/s2

Table 2

Parameters in controller and control command allocation algorithm.

Parameter	Value
$$w_0$$	1.5
$${\mathbf{K}}_p$$	$$\mathrm{diag}6.9,5.56,35$$
$${\mathbf{K}}_i$$	$$\mathrm{diag}0.492,0.588,2.8$$
$${\mathbf{K}}_d$$	$$\mathrm{diag}10.8,7.87,23.3$$
$$\Gamma$$	0.0005
$$k_R$$	3.9
$$k_{\Omega }$$	1.3
$$\epsilon$$	$$10^{-6}$$
$$\delta$$	$$10^{-9}$$
$${\mathbf{w}}_1$$	$$\mathrm{diag}0.5,10,1,1$$
$${\mathbf{w}}_2$$	$$\mathrm{diag}0.5,10,1,1$$

Two experimental tests are performed with different desired trajectories. The experiment video can be seen here: http://youtu.be/Dv1wOyhUnlE. In the first test, the system is expected to move 1 m along $X$ axis and keep the positions in $Y$ and $Z$ directions at 0 and 0.6 m. The desired yaw angle is set as 0 rad. To avoid large overshoot, the desired trajectory in $X$ direction is set as a stair-step signal shown as the blue dashed curve in the first subfigure in Figure 5, i.e., the system is expected to keep moving 0.1 m per second for 10 seconds in $X$ direction. As shown in Figure 5, the transportation system starts taking off at $t=0$. The entire system arrives at the hover status at 5 s. Due to the time delay of sending control command from the first quadrotor to the second one, and other unknown disturbances and uncertainties, a disturbance torque in roll direction is introduced in the take-off stage. Such a disturbance torque causes the system to move 13.2 cm in $Y$ direction. As shown in the second subfigure in Figure 5, the designed control method is able to drive the system back to the hover state. At about 7.4 s, the reference trajectory is sent to the first quadrotor. Finally, all tracking errors in $X$, $Y$, and $Z$ directions are within 5 cm, which validates the effectiveness of the proposed control method. Figure 6 gives the snapshots of this test.

[figure(s) omitted; refer to PDF]

In the second test, the system will only move 1.8 m in 10 seconds in $Y$ direction with the desired height 0.6 m and the desired zero yaw angle. Note that both thrust force and attitude maneuver of these two quadrotors can help finish the movement along $Y$ axis. As shown in Figure 7, the designed control and command allocation methods can complete the transportation task with a tracking error of 5 cm. The snapshots of the test in $Y$ direction are given in Figure 8.

[figure(s) omitted; refer to PDF]

In some sense, the simulation with 4 drones with different orientations is the easiest case for the designed controller because the matrix $\mathbf{A}{\mathbf{W}}^{-1}{\mathbf{A}}^T$ is nonsingular, which implies that the thrust forces of four quadrotors are sufficient to contribute to the movement in both $x$ and $y$ directions. However, two quadrotors are used to verify the developed methods in this work. In our experimental testbed, the thrust forces of both quadrotors contribute to the movement in $y$ direction. However, the movement along $x$ axis can only be realized by the attitude maneuver in the pitch direction. That is, the presented experimental result in $y$ direction can verify the effectiveness of the developed method for the case where the thrust forces can contribute to the translational movement and the result in $x$ direction can be used to verify the trivial case.

5. Conclusions

In this study, a general framework of the dynamics modeling, controller design, and control command allocation method is presented for the transportation of a payload attached to quadrotors rigidly. Since quadrotors may have different orientations, horizontal control forces can be obtained by adjusting the thrust forces of some quadrotors. Hence, different from the case for the classical quadrotors that the movement in the horizontal plane can only be realized based on attitude maneuvers in roll and pitch directions, the transportation system in this study can move in the horizontal plane by regulating some quadrotors’ thrust forces. To get a universal control system, a hierarchical controller and a modified control command allocation algorithm is adopted. Experimental results show the effectiveness of the proposed control system. This work is the first trial to solve the dynamics modeling, controller design, and control command allocation problems for the transportation of a rigidly connected payload using a team of quadrotors with different orientations. In the future, the dynamics and control of such a transportation system can be studied with the uncertainties and communication delays between quadrotors.

Acknowledgments

This work was supported by the Science and Technology on Space Intelligent Control Laboratory (Grant No. 2021-JCJQ-LB-010-17) and the National Natural Science Foundation of China under Grant 12102174.

Appendix

A. Detailed Expressions of the Total Inertia Matrix

Suppose that a new frame $x^{\prime }{y}^{\prime }{z}^{\prime }$ is obtained by translating the origin of $xyz$ to the point $C_1$ in Figure 2, where $C_1$ is the mass center of the payload. Hence, $x^{\prime }{y}^{\prime }{z}^{\prime }$ can be called the payload body frame. Denote the inertia matrices of the payload and quadrotors with respect to their body frames as ${\mathbf{J}}_p$, ${\mathbf{J}}_1$, $\cdots$, ${\mathbf{J}}_n$. According to [31], the inertia matrix of the quadrotor $i$ can be expressed in the frame $xyz$ as $$\begin{array}{}\text{(A.1)}& \mathbf{J}{\prime }_i={\mathbf{R}}_i{\mathbf{J}}_i{\mathbf{R}}_i^T-m_q{\mathbf{r}}_i^{\times }{\mathbf{r}}_i^{\times }.\end{array}$$

Similarly, in the frame $xyz$, the payload inertia matrix is $$\begin{array}{}\text{(A.2)}& \mathbf{J}{\prime }_p={\mathbf{J}}_p-m_p{\mathbf{r}}_p^{\times }{\mathbf{r}}_p^{\times },\end{array}$$

where ${\mathbf{r}}_p$ represents the vector from $C$ to $C_1$ in the frame $xyz$. Hence, the total inertia matrix of the entire system in the frame $xyz$ is $$\begin{array}{}\text{(A.3)}& \mathbf{J}=\mathbf{J}{\prime }_p+\sum _{i=1}^n\mathbf{J}{\prime }_i.\end{array}$$

B. Proof of Theorem 4

With the controller in Equation (13), the closed-loop system can be written as $$\begin{array}{}\text{(B.1)}& m{\ddot{\mathbf{e}}}_c=-{\mathbf{K}}_p+{\mathbf{K}}_i{\mathbf{e}}_c-{\mathbf{K}}_i{\int }_0^t{\mathbf{e}}_{c}dt-{\mathbf{K}}_d{\dot{\mathbf{e}}}_c+m{\tilde{\mathbf{q}}}_3,\end{array}$$where ${\mathbf{e}}_c={\mathbf{r}}_c-{\mathbf{r}}_d$. According to Lemma 3, ${\tilde{\mathbf{q}}}_3$ is bounded by a small constant $c_{\Delta }$ with a large $w_0$. Equation (B.1) can be rewritten as follows: $$\begin{array}{}\text{(B.2)}& \dot{\mathbf{q}}=\mathcal{A}\mathbf{q}+\mathbf{\Delta },\end{array}$$where $$\begin{array}{c}\text{(B.3)}& \mathbf{q}={{{\int }_0^t{\mathbf{e}}_{c}dt}^T,{\mathbf{e}}_c^T,{\dot{\mathbf{e}}}_c^T}^T,\mathcal{A}=\begin{array}{ccc}\mathbf{0}& {\mathbf{I}}_3& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{I}}_3\\ -\frac{{\mathbf{K}}_i}{m}& -\frac{{\mathbf{K}}_p+{\mathbf{K}}_i}{m}& -\frac{{\mathbf{K}}_d}{m}\end{array},\\ \mathbf{\Delta }={{\mathbf{0}}^T,{\mathbf{0}}^T,{\tilde{\mathbf{q}}}_3^T}^T.\end{array}$$

Since all eigenvalues of $\mathcal{A}$ have negative real parts, the following Lyapunov equation holds: $$\begin{array}{}\text{(B.4)}& \mathcal{P}\mathcal{A}+{\mathcal{A}}^T\mathcal{P}=-\mathbf{I},\end{array}$$

where the matrix $\mathcal{P}$ is positive definite. Choose the following Lyapunov function: $$\begin{array}{}\text{(B.5)}& V_p={\mathbf{q}}^T\mathcal{P}\mathbf{q},\end{array}$$whose time derivative is $$\begin{array}{}\text{(B.6)}& {\dot{V}}_p={\dot{\mathbf{q}}}^T\mathcal{P}\mathbf{q}+{\mathbf{q}}^T\mathcal{P}\dot{\mathbf{q}}=-{\mathit{q}}^T\mathit{q}+{\mathit{\Delta }}^T\mathcal{P}\mathit{q}+{\mathit{q}}^T\mathcal{P}\mathit{\Delta }\le -\parallel \mathit{q}{\parallel }_{\mathbf{2}}^{\mathbf{2}}+\mathbf{2}\parallel \mathit{\Delta }{\parallel }_{\mathbf{2}}\parallel \mathcal{P}{\parallel }_{\mathbf{2}}\parallel \mathit{q}{\parallel }_{\mathbf{2}}.\end{array}$$

In the case of $\parallel \mathbf{q}{\parallel }_2>2\parallel \mathbf{\Delta }{\parallel }_2\parallel \mathcal{P}{\parallel }_2$, ${\dot{V}}_p<0$ holds. Hence, $\mathbf{q}$ will converge to the set $\mathbf{q}\parallel \mathbf{q}{\parallel }_2\le 2\parallel \mathbf{\Delta }{\parallel }_2\parallel \mathcal{P}{\parallel }_2$. According to Lemma 3, by increasing $w_0$, one can get a smaller tracking error.

C. Proof of Theorem 8

The closed-loop attitude dynamics is $$\begin{array}{}\text{(C.1)}& \mathbf{J}{\dot{\mathbf{E}}}_{\Omega }=-k_R{\mathbf{e}}_R-k_{\Omega }{\mathbf{E}}_{\Omega }-\mathbf{J}{\widehat{\mathbf{d}}}_2+{\mathbf{d}}_2.\end{array}$$

Consider the system governed by the following equation: $$\begin{array}{}\text{(C.2)}& \mathbf{J}{\dot{\mathbf{E}}}_{\Omega }=-k_R{\mathbf{e}}_R-k_{\Omega }{\mathbf{E}}_{\Omega }.\end{array}$$

Choose the following Lyapunov function: $$\begin{array}{}\text{(C.3)}& V_a=\frac{k_R}{2}\parallel {\mathbf{E}}_R{\parallel }_F^2+\frac{1}{2}{\mathbf{E}}_{\Omega }^T\mathbf{J}{\mathbf{E}}_{\Omega }+\mu {\mathbf{e}}_R^T{\mathbf{E}}_{\Omega },\end{array}$$

where ${\mathbf{E}}_R=\mathbf{R}-{\mathbf{R}}_d$, $\mu$ is a constant satisfying $0<\mu <\sqrt{k_R{\lambda }_{\min }\mathbf{J}/2}$, and $\parallel ·{\parallel }_F$ represents the Frobenius norm of a matrix. $\parallel {\mathbf{E}}_R{\parallel }_F$ has the following property: $$\begin{array}{c}\text{(C.4)}& \parallel {\mathbf{E}}_R{\parallel }_F^2=\parallel \mathbf{R}-{\mathbf{R}}_d{\parallel }_F^2=\text{trace}{\mathbf{R}-{\mathbf{R}}_d}^T\mathbf{R}-{\mathbf{R}}_d=\text{trace}{\mathbf{R}-{\mathbf{R}}_d}^T\mathbf{R}{\mathbf{R}}^T\mathbf{R}-{\mathbf{R}}_d=\parallel \mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d{\parallel }_F^2,\frac{d}{dt}\parallel {\mathbf{E}}_R{\parallel }_F^2=\frac{d}{dt}\text{trace}2\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d-{\mathbf{R}}_d^T\mathbf{R}=2\text{trace}{\mathbf{\Omega }}^{\times }{\mathbf{R}}^T{\mathbf{R}}_d+{\mathbf{\Omega }}_d^{\times }{\mathbf{R}}_d^T\mathbf{R}=-2{\mathbf{\Omega }}^T{{\mathbf{R}}^T{\mathbf{R}}_d-{\mathbf{R}}_d^T\mathbf{R}}^{\vee }-2{\mathbf{\Omega }}_d^T{{\mathbf{R}}_d^T\mathbf{R}-{\mathbf{R}}^T{\mathbf{R}}_d}^{\vee }=2{\mathbf{\Omega }}^T{\mathbf{e}}_R-2{\mathbf{\Omega }}_d^T{\mathbf{e}}_R=2{\mathbf{E}}_{\Omega }^T{\mathbf{e}}_R.\end{array}$$

Using Rodrigues’ rotation formula, one has $$\begin{array}{c}\text{(C.5)}& \parallel {\mathbf{E}}_R{\parallel }_F^2=41-\cos \Theta ,\parallel {\mathbf{e}}_R{\parallel }_2=2\sin \Theta ,\end{array}$$where $\Theta$ is the rotation angle in Rodrigues’ rotation formula associated with the attitude motion described by ${\mathbf{R}}^T{\mathbf{R}}_d$. Hence, the following statement holds with the assumption that $\cos \Theta >-1$: $$\begin{array}{}\text{(C.6)}& \parallel {\mathbf{E}}_R{\parallel }_F^2\le \alpha \parallel {\mathbf{e}}_R{\parallel }_2^2,\end{array}$$where the constant $\alpha >1$ is large enough. Furthermore, it is straightforward to show that $$\begin{array}{}\text{(C.7)}& \parallel {\mathbf{e}}_R{\parallel }_2^2\le 2\parallel {\mathbf{E}}_R{\parallel }_F^2.\end{array}$$

Hence, $V_a\ge k_R/4\parallel {\mathbf{e}}_R{\parallel }_2^2+1/2{\mathbf{E}}_{\Omega }^T\mathbf{J}{\mathbf{E}}_{\Omega }+\mu {\mathbf{e}}_R^T{\mathbf{E}}_{\Omega }\ge k_R/4\parallel {\mathbf{e}}_R{\parallel }_2^2+1/2{\lambda }_{\min }\mathbf{J}\parallel {\mathbf{E}}_{\Omega }{\parallel }_2^2-\mu \parallel {\mathbf{e}}_R{\parallel }_2\parallel {\mathbf{E}}_{\Omega }{\parallel }_2={\mathbf{z}}^T{\mathbf{W}}_0\mathbf{z}$ holds, where $$\begin{array}{c}\text{(C.8)}& {\mathbf{W}}_0=\begin{array}{cc}\frac{k_R}{4}& -\frac{\mu }{2}\\ -\frac{\mu }{2}& \frac{1}{2}{\lambda }_{\min }\mathbf{J}\end{array},\mathbf{z}={\parallel {\mathbf{e}}_R{\parallel }_2,\parallel {\mathbf{E}}_{\Omega }{\parallel }_2}^T.\end{array}$$With $0<\mu <\sqrt{k_R{\lambda }_{\min }\mathbf{J}/2}$, ${\mathbf{W}}_0$ is positive definite and $V_a\ge 0$ holds. If and only if $\parallel {\mathbf{e}}_R{\parallel }_2=0$, $\parallel {\mathbf{E}}_{\Omega }{\parallel }_2=0$, and $\parallel {\mathbf{E}}_R{\parallel }_F=0$, $V_a$ equals to 0. Furthermore, $V_a\le k_R/2\alpha \parallel {\mathbf{e}}_R{\parallel }_2^2+1/2{\lambda }_{\max }\mathbf{J}\parallel {\mathbf{E}}_{\Omega }{\parallel }_2+\mu \parallel {\mathbf{e}}_R{\parallel }_2\parallel {\mathbf{E}}_{\Omega }{\parallel }_2={\mathbf{z}}^T{\mathbf{W}}_1\mathbf{z}$ with $$\begin{array}{}\text{(C.9)}& {\mathbf{W}}_1=\begin{array}{cc}\frac{k_R\alpha }{2}& \frac{\mu }{2}\\ \frac{\mu }{2}& \frac{1}{2}{\lambda }_{\max }\mathbf{J}\end{array}.\end{array}$$If $0<\mu <\sqrt{k_R{\lambda }_{\min }\mathbf{J}/2}$ holds, it is obvious that ${\mathbf{W}}_1$ is positive definite.

${\dot{V}}_a$ can be expressed as $$\begin{array}{}\text{(C.10)}& {\dot{V}}_a=k_R{\mathbf{E}}_{\Omega }^T{\mathbf{e}}_R+{\mathbf{E}}_{\Omega }^T-k_R{\mathbf{e}}_R-k_{\Omega }{\mathbf{E}}_{\Omega }+\mu {\text{trace}{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d{\mathbf{E}}_{\Omega }+{\mathbf{e}}_R\times {\mathbf{E}}_{\Omega }}^T{\mathbf{E}}_{\Omega }+\mu {\mathbf{e}}_R^T{\dot{\mathbf{E}}}_{\Omega }=-k_{\Omega }{\mathit{E}}_{\Omega }^T{\mathit{E}}_{\Omega }+\mu {\text{trace}{\mathit{R}}^T{\mathit{R}}_d\mathit{I}-{\mathit{R}}^T{\mathit{R}}_d{\mathit{E}}_{\Omega }}^T{\mathit{E}}_{\Omega }+\mu {\mathit{e}}_R^T{\mathit{J}}^{-1}-k_R{\mathit{e}}_R-k_{\Omega }{\mathit{E}}_{\Omega }.\end{array}$$

It is straightforward from Rodrigues’ rotation formula to show that $\parallel \text{trace}{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d{\parallel }_F=2$. Hence, one has $\parallel \text{trace}{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d{\parallel }_2\le \parallel \text{trace}{\mathbf{R}}^T{\mathbf{R}}_d\mathbf{I}-{\mathbf{R}}^T{\mathbf{R}}_d{\parallel }_F=2$. ${\dot{V}}_a$ satisfies the following inequality: $$\begin{array}{}\text{(C.11)}& {\dot{V}}_a\le -k_{\Omega }\parallel {\mathbf{E}}_{\Omega }{\parallel }_2^2+2\mu \parallel {\mathbf{E}}_{\Omega }{\parallel }_2^2-k_R\mu {\mathbf{e}}_R^T{\mathbf{J}}^{-1}{\mathbf{e}}_R-k_{\Omega }\mu {\mathbf{e}}_R^T{\mathbf{J}}^{-1}{\mathbf{E}}_{\Omega }\le -k_{\Omega }-2\mu \parallel {\mathit{E}}_{\Omega }{\parallel }_2^2-k_R\mu {\lambda }_{\min }{\mathit{J}}^{-1}\parallel {\mathit{e}}_R{\parallel }_2^2+k_{\Omega }\mu {\lambda }_{\max }{\mathit{J}}^{-1}\parallel {\mathit{e}}_R{\parallel }_2{\mathit{J}}^{-1}\parallel {\mathit{E}}_{\Omega }{\parallel }_2=-{\mathit{z}}^T{\mathit{W}}_2\mathit{z},\end{array}$$where $$\begin{array}{}\text{(C.12)}& {\mathbf{W}}_2=\begin{array}{cc}{k}_R\mu {\lambda }_{\min }{\mathbf{J}}^{-1}& -\frac{k_{\Omega }\mu {\lambda }_{\max }{\mathbf{J}}^{-1}}{2}\\ -\frac{k_{\Omega }\mu {\lambda }_{\max }{\mathbf{J}}^{-1}}{2}& k_{\Omega }-2\mu \end{array}.\end{array}$$There must exist a positive constant $0<\sigma \le {\lambda }_{\min }{\mathbf{W}}_2/{\lambda }_{\max }{\mathbf{W}}_1$ such that ${\mathbf{W}}_2-\sigma {\mathbf{W}}_1\ge 0$ holds. Therefore, one has $$\begin{array}{}\text{(C.13)}& {\dot{V}}_a\le -\sigma V_a.\end{array}$$

Therefore, the closed-loop system in Equation (B.3) is exponentially stable except at some points with $\cos \Theta =-1$.

As shown in Lemma 7, $-\mathbf{J}{\widehat{\mathbf{d}}}_2+{\mathbf{d}}_2$ is ultimately bounded. According to [32], the designed controller can drive the system to a small neighborhood of the desired attitude almost globally.

References

[1] J. Euchi, "Do drones have a realistic place in a pandemic fight for delivering medical supplies in healthcare systems problems?," Chinese Journal of Aeronautics, vol. 34 no. 2, pp. 182-190, DOI: 10.1016/j.cja.2020.06.006, 2021.

[2] D. Xilun, G. Pin, X. Kun, Y. Yushu, "A review of aerial manipulation of small-scale rotorcraft unmanned robotic systems," Chinese Journal of Aeronautics, vol. 32 no. 1, pp. 200-214, DOI: 10.1016/j.cja.2018.05.012, 2019.

[3] A. Gupta, T. Afrin, E. Scully, N. Yodo, "Advances of UAVs toward future transportation: The State-of-the-Art, challenges, and Opportunities," Transportation, vol. 1 no. 2, pp. 326-350, DOI: 10.3390/futuretransp1020019, 2021.

[4] D. K. Villa, A. S. Brandao, M. Sarcinelli-Filho, "A survey on load transportation using multirotor UAVs," Journal of Intelligent and Robotic Systems, vol. 98 no. 2, pp. 267-296, DOI: 10.1007/s10846-019-01088-w, 2020.

[5] I. Palunko, P. Cruz, R. Fierro, "Agile load transportation: safe and efficient load manipulation with aerial robots," IEEE Robotics and Automation Magazine, vol. 19 no. 3, pp. 69-79, DOI: 10.1109/MRA.2012.2205617, 2012.

[6] F. Ruggiero, V. Lippiello, A. Ollero, "Aerial manipulation: a literature review," IEEE Robotics and Automation Letters, vol. 3 no. 3, pp. 1957-1964, DOI: 10.1109/LRA.2018.2808541, 2018.

[7] T. Chen, J. Shan, "A novel cable-suspended quadrotor transportation system: from theory to experiment," Aerospace Science and Technology, vol. 104, article 105974,DOI: 10.1016/j.ast.2020.105974, 2020.

[8] L. Qian, H. H. Liu, "Path-following control of a quadrotor UAV with a cable-suspended payload under wind disturbances," IEEE Transactions on Industrial Electronics, vol. 67 no. 3, pp. 2021-2029, DOI: 10.1109/TIE.2019.2905811, 2020.

[9] P. J. Cruz, R. Fierro, "Cable-suspended load lifting by a quadrotor UAV: hybrid model, trajectory generation, and control," Autonomous Robots, vol. 41 no. 8, pp. 1629-1643, DOI: 10.1007/s10514-017-9632-2, 2017.

[10] A. Tagliabue, M. Kamel, R. Siegwart, J. Nieto, "Robust collaborative object transportation using multiple MAVs," The International Journal of Robotics Research, vol. 38 no. 9, pp. 1020-1044, DOI: 10.1177/0278364919854131, 2019.

[11] T. Chen, J. Shan, H. H. Liu, "Cooperative transportation of a flexible payload using two quadrotors," Journal of Guidance, Control, and Dynamics, vol. 44 no. 11, pp. 2099-2107, DOI: 10.2514/1.G005914, 2021.

[12] T. Chen, J. Shan, "Cooperative transportation of cable-suspended slender payload using two quadrotors," 2019 IEEE International Conference on Unmanned Systems (ICUS), pp. 432-437, 2019.

[13] N. Michael, J. Fink, V. Kumar, "Cooperative manipulation and transportation with aerial robots," Autonomous Robots, vol. 30 no. 1, pp. 73-86, DOI: 10.1007/s10514-010-9205-0, 2011.

[14] D. Sanalitro, H. J. Savino, M. Tognon, J. Cort’es, A. Franchi, "Full-pose manipulation control of a cable-suspended load with multiple UAVs under uncertainties," IEEE Robotics and Automation Letters, vol. 5 no. 2, pp. 2185-2191, DOI: 10.1109/LRA.2020.2969930, 2020.

[15] J. Geng, J. W. Langelaan, "Cooperative transport of a slung load using load-leading control," Journal of Guidance, Control, and Dynamics, vol. 43 no. 7, pp. 1313-1331, DOI: 10.2514/1.G004680, 2020.

[16] Q. Jiang, V. Kumar, "The inverse kinematics of cooperative transport with multiple aerial robots," IEEE Transactions on Robotics, vol. 29, pp. 136-145, 2013.

[17] F. A. Goodarzi, T. Lee, "Stabilization of a Rigid Body Payload With Multiple Cooperative Quadrotors," Journal of Dynamic Systems, Measurement, and Control, vol. 138 no. 12, article 121001,DOI: 10.1115/1.4033945, 2016.

[18] L. Qian, H. H. Liu, "Path following control of multiple quadrotors carrying a rigid-body slung payload," AIAA Scitech 2019 Forum, 2019.

[19] H. Lee, H. Kim, W. Kim, H. J. Kim, "An integrated framework for cooperative aerial manipulators in unknown environments," IEEE Robotics and Automation Letters, vol. 3 no. 3, pp. 2307-2314, DOI: 10.1109/LRA.2018.2807486, 2018.

[20] Y. Qi, J. Wang, J. Shan, "Aerial cooperative transporting and assembling control using multiple quadrotor–manipulator systems," International Journal of Systems Science, vol. 49 no. 3, pp. 662-676, DOI: 10.1080/00207721.2017.1412538, 2018.

[21] H. Lee, H. Kim, H. J. Kim, "Planning and control for collision-free cooperative aerial transportation," IEEE Transactions on Automation Science and Engineering, vol. 15 no. 1, pp. 189-201, DOI: 10.1109/TASE.2016.2605707, 2018.

[22] D. Mellinger, M. Shomin, N. Michael, V. Kumar, "Cooperative grasping and transport using multiple quadrotors," Distributed autonomous robotic systems, pp. 545-558, 2013.

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

[24] R. Ritz, R. D’Andrea, "Carrying a flexible payload with multiple flying vehicles," 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 3465-3471, 2013.

[25] Q. Zheng, L. Q. Gaol, Z. Gao, "On stability analysis of active disturbance rejection control for nonlinear time-varying plants with unknown dynamics," 2007 46th IEEE conference on decision and control, pp. 3501-3506, 2007.

[26] J. Yao, Z. Jiao, D. Ma, "Adaptive robust control of DC motors with extended state observer," IEEE Transactions on Industrial Electronics, vol. 61, pp. 3630-3637, 2014.

[27] T. Chen, J. Shan, H. Wen, "Distributed adaptive attitude control for networked underactuated flexible spacecraft," IEEE Transactions on Aerospace and Electronic Systems, vol. 55, pp. 215-225, 2019.

[28] B.-Z. Guo, Z. L. Zhao, "On the convergence of an extended state observer for nonlinear systems with uncertainty," Systems & Control Letters, vol. 60 no. 6, pp. 420-430, DOI: 10.1016/j.sysconle.2011.03.008, 2011.

[29] Z. Zheng, M. Shen, "Inertial vector measurements based attitude synchronization control for multiple spacecraft formation," Aerospace Science and Technology, vol. 93, article 105309,DOI: 10.1016/j.ast.2019.105309, 2019.

[30] QUARC, "real-time control software," 2021. https://www.quanser.com/products/quarc-real-time-control-software/ Accessed November 16, 2021

[31] F. P. Beer, E. R. Johnston, Vector Mechanics for Engineers, 2004.

[32] D. Angeli, L. Praly, "Stability robustness in the presence of exponentially unstable isolated equilibria," IEEE Transactions on Automatic Control, vol. 56 no. 7, pp. 1582-1592, DOI: 10.1109/TAC.2010.2091170, 2011.