# Transportation of Payload Using Multiple Quadrotors via Rigid Connection

Authors: Ti Chen; Jinjun Shan; Hugh H. T. Liu
Journal: International Journal of Aerospace Engineering
 (2022)

---

## 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.


---

## Body

## 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 Figure1. 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 developedFigure 1
Payload examples.The remainder of this paper is organized as follows. Section2 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 Figure2, 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 xiyizi of quadrotor i for i=1,2,⋯,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 rc=xc,yc,zcT in the inertial frame. The attitude of the entire system is described by the rotation matrix R∈SO3 from the frame xyz to the inertial frame and R0=I3 holds. Note that the Lie group SO3=R∈ℝ3×3RTR=RRT=I3,detR=1 represents the group of 3×3 orthogonal matrices with a determinant of one. In the body frame of the entire system, the position of the mass center of the ith quadrotor is represented by ri=xi,yi,ziT. The rotation matrix from the body frame of the ith quadrotor to the frame xyz is denoted by Ri∈SO3. Since n quadrotors are attached to the payload rigidly, both ri and Ri are constant.Figure 2
The cooperative transportation system.Essentially, the entire system shown in Figure2 can be considered an unmanned aerial vehicle driven by 4n 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:
(1)ṙc=vc,(2)mv̇c=−mgzI+Rf+d1,(3)Ṙ=RΩ×,(4)JΩ̇=−Ω×JΩ+τ+d2,wherem=mp+n×mq is the total mass of the whole system, in which mp is the payload mass and mq is the quadrotor mass. 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. vc is the velocity vector of point C in the inertial frame. g is the constant gravitational acceleration. zI=0,0,1T. Ω=ω1,ω2,ω3T is the angular velocity in the body frame xyz. Note that for any 3-dimensional vector z=z1−z2z3T∈ℝ3,z× is defined as the skew-symmetric matrix 0−z3z2;z30−z1;−z2z10. d1 and d2 are unmodeled terms. f and τ are the total force and torque vectors provided by these n quadrotors expressed in the system body frame xyz. Their detailed definitions are
(5)fτ=∑i=1nAiui=Au,where
(6)Ai=Rizbi03×3ri×RizbiRi∈ℝ6×4,ui=fi,τiTT∈ℝ4 is the vector consisting of the thrust force and the 3-dimensional control torques of the ith quadrotor, A=A1A2⋯An∈ℝ6×4n and u=u1T,u2T,⋯,unTT∈ℝ4n. Note that zbi=0,0,1T and ui is described in the body frame of the ith 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 framexyz provided by the ith quadrotor can be expressed as Rizbifi. Hence, the total force vector from these n quadrotors expressed in the inertial frame is R·∑i=1nRizbifi. Essentially, R·Ri, the rotation matrix from the body frame of the ith quadrotor to the inertial frame, can be written as follows:
(7)R·Ri=cψ¯icθ¯i−sψ¯icϕ¯i+cψ¯isθ¯isϕ¯isψ¯isϕ¯i+cψ¯isθ¯icϕ¯isψ¯icθ¯icψ¯icϕ¯i+sψ¯isθ¯isϕ¯i−cψ¯isϕ¯i+sψ¯isθ¯icϕ¯i−sθ¯icθ¯isϕ¯icθ¯icϕ¯i,where sp and cp represent sinp and cosp with p=ϕ¯i,θ¯i,ψ¯i. Note that ϕ¯i, θ¯i, and ψ¯i are the roll, pitch, and yaw angles associated with the rotation described by R·Ri.Hence, the lift provided by thesen quadrotors can be written as ∑i=1ncθ¯icϕ¯ifi. The forces acting on the system along X and Y axes are ∑i=1nsψ¯isϕ¯i+cψ¯isθ¯icϕ¯ifi and ∑i=1n−cψ¯isϕ¯i+sψ¯isθ¯icϕ¯ifi. One easy way to check whether the whole system can take off successfully is to ensure that ∑i=1ncθ¯icϕ¯ifimax is larger than the system gravity with the roll angle ϕ¯i and pitch angle θ¯i in the initial configuration and the maximum thrust fimax of each quadrotor. Another important conclusion is that ϕ¯i and θ¯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 termsd1 and d2 is bounded.The unmodeled termsd1 and d2 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 arerd=xd,yd,zd and ψd.


### 3.1. Position Control

Rewrite dynamics Equation (2) as
(8)mv̇c=−mgzI+fd+d1,where fd=Rf 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 f=RTfd. Consider the following linear extended state observer to estimate unmodeled terms d1:
(9)q^̇1=q^2+3w0rc−q^1,q^̇2=q^3+3w02rc−q^1+fdm−gzI,q^̇3=w03rc−q^1,where w0 is a positive observer gain. Letting q~1=rc−q^1, q~2=ṙc−q^2, and q~3=d1/m−q^3 yields
(10)q~̇1=q~2−3w0q~1,q~̇2=q~3−3w02q~1,q~̇3=−w03q~1+ḋ1m.The above equation can be rewritten as(11)q~̇=A¯q~+d¯,where
(12)q~=q~1T,q~2T,q~3TT,d¯=0T,0T,ḋ1T/m,A¯=−3w0I3I30−3w02I30I3−w03I300.It is easy to verify all eigenvalues ofA are negative with a positive constant w0, which implies that q~̇=A¯q~ is asymptotically stable. Therefore, with a finite d¯, q~ is bounded.According to [25, 26], one has the following lemma.Lemma 3 (see [25, 26]).
With Assumption1 and w0>0, q~1, q~2, and q~3 are bounded and there are a constant σi,j>0 and a finite time T1 such that q~i,j≤σi,j for any t>T1, where q~i,j represents the jth element of q~i for i=1,2,3 and j=1,2,3 and σi,j=O1/w0k for some positive integer k.Essentially,q^1, q^2, and q^3 are the estimates of rc, ṙc, and d1/m. Based on the extended state observer, the position controller can be designed as
(13)fd=mgzI−Kp+Kirc−rd−Ki∫0trc−rddt−Kdvc−ṙd+mr¨d−mq^3,where Kp, Ki, and Kd are three diagonal positive definite matrices and ∫0tadt is defined as ∫0ta1dt,∫0ta2dt,∫0ta3dtT for a vector a=a1,a2,a3T.Theorem 4.
Under the controllerfd in Equation (13), if these diagonal positive definite matrices Kp, Ki, and Kd are chosen such that all eigenvalues of
(14)0I3000I3−Kim−Kp+Kim−Kdmhave negative real parts, the position tracking error is bounded.Proof.
See AppendixB.Remark 5.
As shown in the proof of Theorem4, a large ω0 will result in a smaller tracking error. However, a too large ω0 will make the observer sensitive to noise [27]. Hence, ω0 should be chosen to balance the tracking performance and the noise tolerance.Remark 6.
Theorem4 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 Kp+Ki, Kd, and Ki as the proportional, derivative, and integral gains. Hence, the control gains Kp+Ki, Kd, and Ki can be adjusted according to the PID controller tuning methods.


### 3.2. Attitude Control

In the last section, to realize the desired control forcefd, the control force f can be simply chosen as f=RTfd. 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
(15)Rd=cψdcθd−sψdcϕd+cψdsθdsϕdsψdsϕd+cψdsθdcϕdsψdcθdcψdcϕd+sψdsθdsϕd−cψdsϕd+sψdsθdcϕd−sθdcθdsϕdcθdcϕd,where ϕd, θd, and ψd are the desired roll, pitch, and yaw angles satisfying sϕd=fd,1sψd−fd,2cψd/∥fd∥ and tanθd=fd,1cψd+fd,2sψd/fd,3. Note that fd,i represents the ith element of fd. For such a transportation system, it is reasonable to assume that both roll and pitch angles belong to −π/2,π/2. Hence, ϕd and θd can be set as ϕd=arcsinfd,1sψd−fd,2cψd/∥fd∥ and θd=arctanfd,1cψd+fd,2sψd/fd,3, respectively. It should be noted that ψ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(16)Ω^̇=2ΓΩ−Ω^−J−1Ω×JΩ+J−1τ+d^2,d^̇2=Γ2Ω−Ω^,where Γ is a positive constant. Letting eΩ=Ω−Ω^ and ed2=J−1d2−d^2 yields
(17)ėΩ=−2ΓeΩ+ed2,(18)ėd2=J−1ḋ2−Γ2eΩ.According to Theorem 2.2 in [28], one can have the following lemma.Lemma 7.
With Assumption1 and Γ>0, eΩ and ed2 are convergent in the sense that for any σω∈0,1, there exists εσ∈(0, 1) such that eΩ,i<σω and ed2,i<σω hold for any ε1∈0,εσ and any t∈Tε,∞, where Tε>0 depends on ε1 and eΩ,i, where eΩ,i and ed2,i are the ith elements of eΩ and ed2, respectively.To track the desired attitudeRd, a robust controller is designed as follows:
(19)τ=Ω×JΩ−kReR−kΩEΩ+JΩ̇d−Jd^2,whereeR=RdTR−RTRd∨ and EΩ=Ω−Ωd. Note that ėR=−Ωd×RdTR+RdTRΩ×+Ω×RTRd−RTRdΩd×∨=traceRTRdI−RTRdΩ−traceRdTRI−RdTRΩd=traceRTRdI−RTRdEΩ+eR×EΩ. d^2 is the estimate of the unknown term J−1d2 and updated according to the law in Equation (17).Theorem 8.
If0<μ<kRλminJ/2 holds and
(20)kRμλminJ−1−kΩμλmaxJ−12−kΩμλmaxJ−12kΩ−2μ
is positive definite, the attitude tracking error is bounded with the designed controller exponentially almost globally.Proof.
See AppendixC.Remark 9.
As shown in the proof of Theorem8, 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Θ=−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 f,τTT. The rotation matrix Ri from the ith quadrotor’s body frame to the body frame of the entire system can be written as
(21)Ri=cψicθi−sψicϕi+cψisθisϕisψisϕi+cψisθicϕisψicθicψicϕi+sψisθisϕi−cψisϕi+sψisθicϕi−sθicθisϕicθicϕi,where ϕi, θi, and ψi are the roll, pitch, and yaw angles associated with the rotation Ri.To determine the control inputs for each quadrotor, one can minimize the following cost function(22)Jc=uTWuwith the following equation constraint
(23)Au=f,τTT,where W=diagwi∈ℝ4n×4n is the weight matrix, in which wi=diagwi,j with wi,j>0 for i=1,2,⋯,n and j=1,2,3,4. Essentially, this is a constrained quadratic problem, whose solution is
(24)u=K+f,τTT,where K+=W−1ATAW−1AT−1. Note that AW−1AT can be written as ∑i=1nAiwi−1AiT. Some special cases are discussed below.Case 1.
The weight matrixwi is chosen the same for each quadrotor, and the weights in roll, pitch, and yaw directions are equal. Hence, wi can be rewritten as diagw1,w2,w2,w2. 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 ϕi and θi are small enough such that cosϕi≈1, cosθi≈1, sinϕi≈ϕi, sinθi≈θi, ϕi2≈0, θi2≈0, and ϕiθi≈0 hold. Hence, Aiwi−1AiT can be written as
(25)Aiwi−1AiT=02×2h102×1h1Th2h301×2h3Th4,where
(26)h1=cψiθi+sψiϕiw1yicψiθi+sψiϕiw1−xicψiθi+sψiϕiw1−cψiϕi−sψiθiw1−yicψiϕi−sψiθiw1xicψiϕi−sψiθiw1,h2=w1−1cψiϕizi−sψiθizi+yiw1cψiθizi+sψiϕizi−xiw1cψiϕizi−sψiθizi+yiw12cψiϕiw2yizi−2sψiθiw2yizi+w2yi2+w1w1w2−cψiϕixizi−cψiθiyizi−sψiϕiyizi−sψiθixizi+xiyiw1cψiθizi+sψiϕizi−xiw1−cψiϕixizi−cψiθiyizi−sψiϕiyizi−sψiθixizi+xiyiw1−2cψiθiw2xizi+2sψiϕiw2xizi−w2xi2−w1w1w2,h3=−cψiϕixi+cψiθiyi+sψiϕiyi−sψiθixiw1−yicψiϕixi+cψiθiyi+sψiϕiyi−sψiθixiw1xicψiϕixi+cψiθiyi+sψiϕiyi−sψiθixiw1,h4=w2−1.It should be noted that the determinant ofAiwi−1AiT 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 ϕ1=0.1 rad, ϕ2=−0.05 rad, ϕ3=ϕ4=0, θ1=θ2=0, θ3=0.1 rad, θ4=−0.05 rad, ψ1=ψ2=ψ3=ψ4=0, x1=−x2=1, x3=x4=0, y1=y2=0, y3=−y4=1, z1=z2=z3=z4=0, and w1=w2=1, the determinant of ∑i=14Aiwi−1AiT is 0.01, which implies that these four quadrotors can provide both three-dimensional forces and three-dimensional torques. However, if letting θ3=θ4=0 and keeping all remaining parameters the same as those in the previous sentence, ∑i=14Aiwi−1AiT 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:
(27)u=W−1ATAW−1AT−1fT,τTT,ifDetAW−1AT>ε,W−1ATAW−1AT+δI6−1fT,τTT,otherwise,where ε and δ are small positive constants.Case 2.
If thexiyi plane in the body frame xiyizi is parallel to xy plane in the frame xyz, ϕi=0 and θi=0 hold for each quadrotor and
(28)Ai=02×4A^i,where A^i is
(29)A^i=1000yicψi−sψi0−xisψicψi00001.
Hence,Aiwi−1AiT is
(30)Aiwi−1AiT=00000000000000w1−1yiw1−xiw1000yiw1w2yi2+w1w1w2−xiyiw1000−xiw1−xiyiw1w2xi2+w1w1w2000000w2−1.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, AW−1AT would be
(31)AW−1AT=00000000000000nw1000000∑i=1nw2yi2+w1w1w2000000∑i=1nw2xi2+w1w1w2000000nw2.It is clear that the matrix in the above equation is singular. Based on Equation (27), ATAW−1AT+δI6−1 can be written as
(32)A1AW−1AT+δI6−1⋯AnAW−1AT+δI6−1,where
(33)AiAW−1AT+δI6−1=04×2A^iTdiag0,0,nw1,∑i=1nw2yi2+w1w1w2,∑i=1nw2xi2+w1w1w2,nw2+δI6−1≈04×2A^iTdiagnw1,∑i=1nw2yi2+w1w1w2,∑i=1nw2xi2+w1w1w2,nw2−1.Letf3 be the third element of f. AiAW−1AT+δI6−1fT,τTT can be written as
(34)A^iTdiagnw1,∑i=1nw2yi2+w1w1w2,∑i=1nw2xi2+w1w1w2,nw2−1f3,τTT,which agrees with the result from the method in [22].Another interesting point is that bothW−1ATAW−1AT−1 and W−1ATAW−1AT+δI6−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 Figure3. With the position vector of the mass center of the entire system and the desired transportation destination, Equation (13) will provide the value of fd, 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 3
Control system diagram.


## 3.1. Position Control

Rewrite dynamics Equation (2) as
(8)mv̇c=−mgzI+fd+d1,where fd=Rf 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 f=RTfd. Consider the following linear extended state observer to estimate unmodeled terms d1:
(9)q^̇1=q^2+3w0rc−q^1,q^̇2=q^3+3w02rc−q^1+fdm−gzI,q^̇3=w03rc−q^1,where w0 is a positive observer gain. Letting q~1=rc−q^1, q~2=ṙc−q^2, and q~3=d1/m−q^3 yields
(10)q~̇1=q~2−3w0q~1,q~̇2=q~3−3w02q~1,q~̇3=−w03q~1+ḋ1m.The above equation can be rewritten as(11)q~̇=A¯q~+d¯,where
(12)q~=q~1T,q~2T,q~3TT,d¯=0T,0T,ḋ1T/m,A¯=−3w0I3I30−3w02I30I3−w03I300.It is easy to verify all eigenvalues ofA are negative with a positive constant w0, which implies that q~̇=A¯q~ is asymptotically stable. Therefore, with a finite d¯, q~ is bounded.According to [25, 26], one has the following lemma.Lemma 3 (see [25, 26]).
With Assumption1 and w0>0, q~1, q~2, and q~3 are bounded and there are a constant σi,j>0 and a finite time T1 such that q~i,j≤σi,j for any t>T1, where q~i,j represents the jth element of q~i for i=1,2,3 and j=1,2,3 and σi,j=O1/w0k for some positive integer k.Essentially,q^1, q^2, and q^3 are the estimates of rc, ṙc, and d1/m. Based on the extended state observer, the position controller can be designed as
(13)fd=mgzI−Kp+Kirc−rd−Ki∫0trc−rddt−Kdvc−ṙd+mr¨d−mq^3,where Kp, Ki, and Kd are three diagonal positive definite matrices and ∫0tadt is defined as ∫0ta1dt,∫0ta2dt,∫0ta3dtT for a vector a=a1,a2,a3T.Theorem 4.
Under the controllerfd in Equation (13), if these diagonal positive definite matrices Kp, Ki, and Kd are chosen such that all eigenvalues of
(14)0I3000I3−Kim−Kp+Kim−Kdmhave negative real parts, the position tracking error is bounded.Proof.
See AppendixB.Remark 5.
As shown in the proof of Theorem4, a large ω0 will result in a smaller tracking error. However, a too large ω0 will make the observer sensitive to noise [27]. Hence, ω0 should be chosen to balance the tracking performance and the noise tolerance.Remark 6.
Theorem4 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 Kp+Ki, Kd, and Ki as the proportional, derivative, and integral gains. Hence, the control gains Kp+Ki, Kd, and Ki can be adjusted according to the PID controller tuning methods.


## 3.2. Attitude Control

In the last section, to realize the desired control forcefd, the control force f can be simply chosen as f=RTfd. 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
(15)Rd=cψdcθd−sψdcϕd+cψdsθdsϕdsψdsϕd+cψdsθdcϕdsψdcθdcψdcϕd+sψdsθdsϕd−cψdsϕd+sψdsθdcϕd−sθdcθdsϕdcθdcϕd,where ϕd, θd, and ψd are the desired roll, pitch, and yaw angles satisfying sϕd=fd,1sψd−fd,2cψd/∥fd∥ and tanθd=fd,1cψd+fd,2sψd/fd,3. Note that fd,i represents the ith element of fd. For such a transportation system, it is reasonable to assume that both roll and pitch angles belong to −π/2,π/2. Hence, ϕd and θd can be set as ϕd=arcsinfd,1sψd−fd,2cψd/∥fd∥ and θd=arctanfd,1cψd+fd,2sψd/fd,3, respectively. It should be noted that ψ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(16)Ω^̇=2ΓΩ−Ω^−J−1Ω×JΩ+J−1τ+d^2,d^̇2=Γ2Ω−Ω^,where Γ is a positive constant. Letting eΩ=Ω−Ω^ and ed2=J−1d2−d^2 yields
(17)ėΩ=−2ΓeΩ+ed2,(18)ėd2=J−1ḋ2−Γ2eΩ.According to Theorem 2.2 in [28], one can have the following lemma.Lemma 7.
With Assumption1 and Γ>0, eΩ and ed2 are convergent in the sense that for any σω∈0,1, there exists εσ∈(0, 1) such that eΩ,i<σω and ed2,i<σω hold for any ε1∈0,εσ and any t∈Tε,∞, where Tε>0 depends on ε1 and eΩ,i, where eΩ,i and ed2,i are the ith elements of eΩ and ed2, respectively.To track the desired attitudeRd, a robust controller is designed as follows:
(19)τ=Ω×JΩ−kReR−kΩEΩ+JΩ̇d−Jd^2,whereeR=RdTR−RTRd∨ and EΩ=Ω−Ωd. Note that ėR=−Ωd×RdTR+RdTRΩ×+Ω×RTRd−RTRdΩd×∨=traceRTRdI−RTRdΩ−traceRdTRI−RdTRΩd=traceRTRdI−RTRdEΩ+eR×EΩ. d^2 is the estimate of the unknown term J−1d2 and updated according to the law in Equation (17).Theorem 8.
If0<μ<kRλminJ/2 holds and
(20)kRμλminJ−1−kΩμλmaxJ−12−kΩμλmaxJ−12kΩ−2μ
is positive definite, the attitude tracking error is bounded with the designed controller exponentially almost globally.Proof.
See AppendixC.Remark 9.
As shown in the proof of Theorem8, 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Θ=−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 f,τTT. The rotation matrix Ri from the ith quadrotor’s body frame to the body frame of the entire system can be written as
(21)Ri=cψicθi−sψicϕi+cψisθisϕisψisϕi+cψisθicϕisψicθicψicϕi+sψisθisϕi−cψisϕi+sψisθicϕi−sθicθisϕicθicϕi,where ϕi, θi, and ψi are the roll, pitch, and yaw angles associated with the rotation Ri.To determine the control inputs for each quadrotor, one can minimize the following cost function(22)Jc=uTWuwith the following equation constraint
(23)Au=f,τTT,where W=diagwi∈ℝ4n×4n is the weight matrix, in which wi=diagwi,j with wi,j>0 for i=1,2,⋯,n and j=1,2,3,4. Essentially, this is a constrained quadratic problem, whose solution is
(24)u=K+f,τTT,where K+=W−1ATAW−1AT−1. Note that AW−1AT can be written as ∑i=1nAiwi−1AiT. Some special cases are discussed below.Case 1.
The weight matrixwi is chosen the same for each quadrotor, and the weights in roll, pitch, and yaw directions are equal. Hence, wi can be rewritten as diagw1,w2,w2,w2. 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 ϕi and θi are small enough such that cosϕi≈1, cosθi≈1, sinϕi≈ϕi, sinθi≈θi, ϕi2≈0, θi2≈0, and ϕiθi≈0 hold. Hence, Aiwi−1AiT can be written as
(25)Aiwi−1AiT=02×2h102×1h1Th2h301×2h3Th4,where
(26)h1=cψiθi+sψiϕiw1yicψiθi+sψiϕiw1−xicψiθi+sψiϕiw1−cψiϕi−sψiθiw1−yicψiϕi−sψiθiw1xicψiϕi−sψiθiw1,h2=w1−1cψiϕizi−sψiθizi+yiw1cψiθizi+sψiϕizi−xiw1cψiϕizi−sψiθizi+yiw12cψiϕiw2yizi−2sψiθiw2yizi+w2yi2+w1w1w2−cψiϕixizi−cψiθiyizi−sψiϕiyizi−sψiθixizi+xiyiw1cψiθizi+sψiϕizi−xiw1−cψiϕixizi−cψiθiyizi−sψiϕiyizi−sψiθixizi+xiyiw1−2cψiθiw2xizi+2sψiϕiw2xizi−w2xi2−w1w1w2,h3=−cψiϕixi+cψiθiyi+sψiϕiyi−sψiθixiw1−yicψiϕixi+cψiθiyi+sψiϕiyi−sψiθixiw1xicψiϕixi+cψiθiyi+sψiϕiyi−sψiθixiw1,h4=w2−1.It should be noted that the determinant ofAiwi−1AiT 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 ϕ1=0.1 rad, ϕ2=−0.05 rad, ϕ3=ϕ4=0, θ1=θ2=0, θ3=0.1 rad, θ4=−0.05 rad, ψ1=ψ2=ψ3=ψ4=0, x1=−x2=1, x3=x4=0, y1=y2=0, y3=−y4=1, z1=z2=z3=z4=0, and w1=w2=1, the determinant of ∑i=14Aiwi−1AiT is 0.01, which implies that these four quadrotors can provide both three-dimensional forces and three-dimensional torques. However, if letting θ3=θ4=0 and keeping all remaining parameters the same as those in the previous sentence, ∑i=14Aiwi−1AiT 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:
(27)u=W−1ATAW−1AT−1fT,τTT,ifDetAW−1AT>ε,W−1ATAW−1AT+δI6−1fT,τTT,otherwise,where ε and δ are small positive constants.Case 2.
If thexiyi plane in the body frame xiyizi is parallel to xy plane in the frame xyz, ϕi=0 and θi=0 hold for each quadrotor and
(28)Ai=02×4A^i,where A^i is
(29)A^i=1000yicψi−sψi0−xisψicψi00001.
Hence,Aiwi−1AiT is
(30)Aiwi−1AiT=00000000000000w1−1yiw1−xiw1000yiw1w2yi2+w1w1w2−xiyiw1000−xiw1−xiyiw1w2xi2+w1w1w2000000w2−1.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, AW−1AT would be
(31)AW−1AT=00000000000000nw1000000∑i=1nw2yi2+w1w1w2000000∑i=1nw2xi2+w1w1w2000000nw2.It is clear that the matrix in the above equation is singular. Based on Equation (27), ATAW−1AT+δI6−1 can be written as
(32)A1AW−1AT+δI6−1⋯AnAW−1AT+δI6−1,where
(33)AiAW−1AT+δI6−1=04×2A^iTdiag0,0,nw1,∑i=1nw2yi2+w1w1w2,∑i=1nw2xi2+w1w1w2,nw2+δI6−1≈04×2A^iTdiagnw1,∑i=1nw2yi2+w1w1w2,∑i=1nw2xi2+w1w1w2,nw2−1.Letf3 be the third element of f. AiAW−1AT+δI6−1fT,τTT can be written as
(34)A^iTdiagnw1,∑i=1nw2yi2+w1w1w2,∑i=1nw2xi2+w1w1w2,nw2−1f3,τTT,which agrees with the result from the method in [22].Another interesting point is that bothW−1ATAW−1AT−1 and W−1ATAW−1AT+δI6−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 Figure3. With the position vector of the mass center of the entire system and the desired transportation destination, Equation (13) will provide the value of fd, 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 3
Control system diagram.


## 4. Experimental Verification

As shown in Figure4, 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 4
Experimental system.Magnets are used to connect the quadrotors and the board. As shown in the subfigure at the bottom of Figure4, 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 r1=0.022,0.56,−0.001T m and r2=−0.022,−0.56,0.001T 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 θi=0 for i=1,2.The detailed physical parameters of quadrotor and payload are given in Table1. 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 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.
ParameterValueQuadrotor mass1.121 kgInertia matrix of each quadrotordiag0.01,0.0082,0.0148 kg·m2Maximum payload capability of each quadrotor~0.3 kgPayload mass0.453 kgPayload dimensions inx, y, and z directions0.038m×1.219m×0.013mInertia matrix of payloaddiag0.0561,0.00006,0.0561 kg·m2Constant gravitational acceleration9.81 m/s2Table 2
Parameters in controller and control command allocation algorithm.
ParameterValuew01.5Kpdiag6.9,5.56,35Kidiag0.492,0.588,2.8Kddiag10.8,7.87,23.3Γ0.0005kR3.9kΩ1.3ε10−6δ10−9w1diag0.5,10,1,1w2diag0.5,10,1,1Two 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 5
The position of the mass center of the entire system in the test inX direction.Figure 6
The snapshots of the test inX direction.In the second test, the system will only move 1.8 m in 10 seconds inY 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 7
The position of the mass center of the entire system in the test inY direction.Figure 8
The snapshots of the test inY direction.In some sense, the simulation with 4 drones with different orientations is the easiest case for the designed controller because the matrixAW−1AT 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.



---

*Source: 2486561-2022-04-22.xml*