# Suspended Load Path Tracking Control Strategy Using a Tilt-Rotor UAV

Authors: M. A. Santos; B. S. Rego; G. V. Raffo; A. Ferramosca
Journal: Journal of Advanced Transportation
 (2017)

---

## Abstract

This work proposes a control strategy to solve the path tracking problem of a suspended load carried by a tilt-rotor unmanned aerial vehicle (UAV). Initially, the equations of motion for the multibody mechanical system are derived from the load’s perspective by means of the Euler-Lagrange formulation, in which the load’s position and orientation are chosen as degrees of freedom. An unscented Kalman filter (UKF) is designed for nonlinear state estimation of all the system states, assuming that available information is provided by noisy sensors with different sampling rates that do not directly measure the load’s attitude. Furthermore, a model predictive control (MPC) strategy is proposed for path tracking of the suspended load with stabilization of the tilt-rotor UAV when parametric uncertainties and external disturbances affect the load, the rope’s length and total system mass vary during taking-off and landing, and the desired yaw angle changes throughout the trajectory. Finally, numerical experiments are presented to corroborate the good performance of the proposed strategy.


---

## Body

## 1. Introduction

Firstly used mainly for military purposes, recent advances in technology allowed large-scale production of unmanned aerial vehicles (UAVs) and consequently their civilian usage. Nowadays, UAVs receive great interest from engineers and researchers and have a wide range of applications, such as precision agriculture, cargo transportation, cinematographic filming, search-and-rescue missions, surveillance, fire inspection, and archeology. The most commonly found configurations of UAVs are the fixed-wing, helicopter, and quadrotor ones. Due to the vertical take-off and landing (VTOL) and hovering capabilities of helicopter and quadrotor UAVs, they are used in tasks that require high maneuverability in slow velocities. Thereby, fixed-wing UAVs are used in tasks that require improved forward flight, mostly in open environments.In the last decades, some hybrid aircrafts have drawn substantial attention due to their ability to perform VTOL, hover, and improved forward flights. The tilt-rotor configuration is among the most popular ones, being provided with fixed wings and rotary wings, and is capable of switching between helicopter and airplane flight modes by tilting its thrusters. Inspired by the capabilities of those aircrafts, recent researches are looking into the design of small-scale tilt-rotor UAVs [1–3]. Their hybrid characteristics offer advantages over fixed-wing and quadrotor UAVs, which also come with design and control challenges since they are complex underactuated mechanical systems with highly coupled dynamics.Moreover, an important application of unmanned aerial vehicles is the load transportation in risky and inaccessible zones, allowing dealing with rapid deployment of supplies in search-and-rescue missions [4], vertical replenishment of seaborne vessels [5], and safe landmine detection [6]. However, this kind of task is a challenging subject in terms of modeling and control. Since the payload is often connected to the UAV through a rope, the dynamic behavior of the system varies due to the load’s swing, which can destabilize the whole system if it is not well attenuated. Moreover, the suspended load by a rope adds unactuated degrees of freedom to the system, increasing its underactuation degree. Apart from the exposed, in order to accomplish the load transportation task, the knowledge on the load position is usually required. The problem of estimating the load position then arises, being mostly addressed through visual systems and state estimators [7, 8].Therefore, due to its hybrid capabilities, a tilt-rotor UAV becomes a promising platform for aerial load transportation, providing improved forward speed when compared with rotary-wing UAVs, which is a desired feature for missions requiring rapid deployment. Furthermore, missions that demand precise positioning of the load cannot be addressed by fixed-wing UAVs, since they are not able to perform hover flights. Thereby, the tilt-rotor UAV can be used in both scenarios.In the literature, control objectives for aerial load transportation include path tracking of the aircraft with reduced load’s swing [9–14], obstacle avoidance [15, 16], transportation by multiple aircrafts [17, 18], and path tracking of the suspended load [19–24]. Since large load swing may destabilize the aircraft, many works address the problem of swing-free aerial load transportation. In [9], a control strategy based on a cascade structure is proposed for load transportation using an unmanned helicopter. A delayed feedback controller is introduced in an outer loop to reduce load’s swing, while the inner-loop controller does not take into account the dynamics of the load. Open-loop approaches based on trajectory generation for load transportation using quadrotor UAVs are proposed in [10, 11], in order to obtain reduced load’s swing motion. In [13], a robust nonlinear control strategy is proposed in order to transport a suspended load by a quadrotor UAV along a predefined trajectory while avoiding load’s swing. The controller is designed based on nonlinear H∞ control theory and Lyapunov redesign taking into account the whole dynamics of the system. Some works have also used tilt-rotor UAVs for load transportation. In [12], a nonlinear cascade control strategy is proposed for path tracking of a tilt-rotor UAV with load’s swing improvement, which is composed of three levels of feedback linearization controllers. In [14], a model predictive control (MPC) strategy is studied for path tracking of a tilt-rotor UAV while carrying a load.When precise positioning of the load is of concern, swing-free motion may not be sufficient for the accomplishment of the task; then the problem of path tracking of the suspended load is addressed. For instance, in [19], the authors propose a model-free open-loop approach based on trajectory generation by reinforcement learning for path tracking of the suspended load using a quadrotor UAV. A nonlinear control strategy based on cascade structure and system decoupling is proposed in [20], also for path tracking of a suspended load using a quadrotor UAV. In [22], a nonlinear solution to the suspended load path tracking problem using a quadrotor UAV is presented, in which the authors assume the quadrotor as a system actuated by total thrust and orientation. The control problem is casted into the framework described in [21], which deals with path tracking of underactuated systems driven by directed thrust and angular velocity with a double-integrator structure. By using tilt-rotor UAVs, in [23, 24], the path tracking problem of the suspended load is solved through the design of control and state estimation strategies based on linearized, time-invariant state-space equations, which did not allow yaw angle tracking, neither the occurrence of changes in the load’s mass and rope’s length.This work intends to solve the problem of suspended load path tracking considering a realistic scenario in which the rope’s length and total system mass are not constant during take-off and landing maneuvers, the desired yaw angle varies during the mission execution, and the sensors embedded in the UAV have different sampling rates and do not measure directly the load’s position and orientation. Due to the tilt-rotor UAV highly nonlinear dynamic behavior, the unscented Kalman filter is chosen as state estimation strategy since it is not based on model linearization. Further, to perform path tracking while dealing with yaw regulation and the rope’s length variation without relying only on controller’s robustness, the MPC strategy is selected because its model-based nature allows dealing with those problems in a simple fashion way. Therefore, the main contributions of this work are as follows: (i) a detailed modeling in which the multibody system’s equations of motion are obtained from the load’s point of view, yielding an input-affine state-space representation with the load’s position and orientation as state variables; (ii) a nonlinear estimation strategy dealing with multirate noisy sensors that do not provide directly information about the load’s pose; and (iii) a control strategy that is able to cope with the path tracking of the suspended load considering yaw angle regulation, to incorporate the rope’s length variation to the control problem formulation, and to keep the UAV stabilized while rejecting external disturbances and parametric uncertainties.The remainder of this paper is organized as follows: Section2 develops the equations of motion of the tilt-rotor UAV with suspended load from the load’s perspective, yielding a state-space representation with the load’s position and orientation represented by state variables; Section 3 presents the unscented Kalman Filter design based on the model and measurement equations given by available sensors; Section 4 describes the model predictive control strategy designed to achieve path tracking of the suspended load based on the linearized time-varying model; Section 5 presents the model and design parameters, besides the simulation scenario used to validate the proposed control strategy; in Section 6, the numerical experiments used to demonstrate the good performance of the controller and filter are presented and discussed; finally, Section 7 concludes the work.


## 2. Tilt-Rotor UAV with Suspended Load Modeling

This section deals with the development of the equations of motion of the tilt-rotor UAV with suspended load, from the load’s perspective. The dynamic equations, assuming a multibody system, are obtained through the Euler-Lagrange formulation, in which the coupling between the aircraft and the load is considered naturally.


### 2.1. System Description

The computer-aided design (CAD) model of the tilt-rotor UAV with suspended load is shown in Figure1. The system is regarded as a multibody mechanical system composed of four rigid bodies: (i) the aircraft’s main body, composed of Acrylonitrile Butadiene Styrene (ABS) structure, landing gear, batteries, and electronics; (ii) the right thruster group, composed of the right thruster and its tilting mechanism (a revolute joint); (iii) the left thruster group, composed of the left thruster and its tilting mechanism; and (iv) the suspended load group, comprising the load and the rope. The system is actuated through the aircraft’s thrusters and tilting mechanisms. For modeling purposes, the following assumptions are made:(A1)
The rope is rigid and has negligible mass.(A2)
The rope is connected to the aircraft’s geometric center.(A3)
The main body’s center of mass does not coincide with the aircraft’s geometric center.(A4)
The thrusters groups’ centers of mass are located at their respective tilting axes.Figure 1
The tilt-rotor UAV with suspended load (CAD model).Regarding the assumptions stated above, it is important to describe the main reasons of assumption (A3), which are as follows: (i) the center of mass is vertically displaced in order to improve the pitch moment; and (ii) the center of mass is displaced with respect to thex-axis of the geometric center, allowing obtaining non-null equilibrium points for the angular positions of the tilting mechanisms and pitch angle. This mechanical feature improves the controllability of the system in hover flight, allowing the projection of the thrusts generated by propellers on the x-axis without the need to tilt the thrusters’ group.


### 2.2. Kinematics from the Load’s Perspective

The load’s perspective approach consists in formulating the system’s kinematics regarding the suspended load as a free rigid body, with the aircraft as a multibody mechanical system rigidly coupled to it. Six reference frames are defined, as shown in Figure2: (i) the inertial reference frame, I; (ii) the aircraft’s geometric center frame, B; (iii) the main body center of mass frame, C1; (iv) the right thruster group center of mass frame, C2; (v) the left thruster group center of mass frame, C3; and (vi) the suspended load group center of mass frame, L. The load’s position with respect to the inertial frame I is denoted by ξ≜xyzT. The displacement vector from L to B corresponds to the rope and is expressed in L by dBL≜00lT, with l as the rope’s length. The displacement vectors from B to Ci are model parameters of the tilt-rotor UAV and are denoted by dCiB, expressed in B, with i∈1,2,3.Figure 2
Kinematic definitions, input forces, and torques.The orientation of the load with respect toI is parametrized by Euler angles, η≜ϕθψT, using the ZYX convention about local axes. The associated rotation matrix is given by (1)RLI≜Rz,ψRy,θRx,ϕ=cψcθcψsθsϕ-sψcϕcψsθcϕ+sψsϕsψcθsψsθsϕ+cψcϕsψsθcϕ-cψsϕ-sθcθsϕcθcϕ.On the other hand, the orientation of the aircraft’s geometric center frame with respect toL, corresponding to the orientation of the UAV with respect to the rope, is parametrized by two angles, γ≜γ1γ2T, such that (2)RBL≜Rx,-γ1Ry,-γ2=cγ20-sγ2sγ1sγ2cγ1sγ1cγ2cγ1sγ2-sγ1cγ1cγ2.The reference framesB and C1 are parallel to each other and attached to the same rigid body; thus the relative orientation is null, that is, RC1B=I3×3 (in this work In×n is an identity matrix with dimension n, 0n×m is n by m zeros matrix, and 1n×m is n by m ones matrix). Furthermore, the orientations of the thrusters’ groups with respect to B are described by(3)RC2B≜Rx,-βRy,αR=cαR0sαR-sβsαRcβsβcαR-cβsαR-sβcβcαR,RC3B≜Rx,βRy,αL=cαL0sαLsβsαLcβ-sβcαL-cβsαLsβcβcαL,where αR and αL are the tilting angles of the right and left thrusters, respectively, and β is a fixed inclination angle of the thrusters towards the aircraft geometric center, designed to improve the aircraft’s controllability [25]. The angular velocities of the system are given by ωILL=Wηη˙, ωLBB=Qγ˙, ωBC1C1=03×1, ωBC2C2=ayα˙R, and ωBC3C3=ayα˙L, where, for instance, ωILL denotes the angular velocity of L with respect to I, expressed in L, with(4)Wη≜10-sθ0cϕsϕcθ0-sϕcϕcθ,Q≜-cγ200-1sγ20,and ay≜010T.From the rigid transformations of the system, the forward kinematics of points that belong to each rigid body are given by(5)pLI=ξ+RLIpLL,(6)piI=ξ+RLIdBL+RLIRBLdCiB+RLIRBLRCiBpiCi,where pL is the position of a point that belongs to the suspended load body and pi belongs to the rigid body with attached frame Ci. Taking the time derivatives of (5) and (6) and making use of several properties of skew-symmetric matrices [26] yield(7)p˙LI=ξ˙+RLISpLLTωILL,(8)p˙iI=ξ˙+RLISdBLT+RLIRBLSdCiBTRBLTξ˙++RLIRBLRCiBSpiCiTRBLRCiBTωILLξ˙++RLIRBLSdCiBT+RLIRBLRCiBSpiCiTRCiBTωLBBξ˙++RLIRBLRCiBSpiCiTωBCiCi,where S· denotes an operator that maps a vector to a skew-symmetric matrix [26].The generalized coordinates of the system are chosen according to their degrees of freedom. For the equations of motion describing explicitly the time evolution of the load’s position and orientation, these are included in the generalized coordinates, which are chosen as(9)q≜ξTηTγTαRαLT∈R10.Note that, due to the chosen perspective for the system kinematics, the position and orientation of the aircraft with respect toI are not degrees of freedom of the system. Thus, they are not included in (9), and consequently their time evolution will not be described explicitly by the developed equations of motion.


### 2.3. Kinetic and Potential Energies

Initially, in order to derive the equations of motion through the Euler-Lagrange formulation, the kinetic and potential energies of each body of the mechanical system must be obtained. For eachith rigid body, these energies can be computed through the volume integrals [27](10)Ki=12∫Viρip˙iITp˙iIdVi,(11)Ui=-∫Viρig^TpiIdVi=-mig^ToCiI,respectively, where ρi corresponds to its density, Vi to its volume, and mi=∫ViρiVi‍ to its mass, g^≜00-g^zT denotes the gravitational acceleration vector expressed in I, and oCiI is the position vector obtained from the forward kinematics of the origin of Ci.The quadratic termsp˙LITp˙LI and p˙iITp˙iI,i∈1,2,3, are obtained using (7) and (8), respectively. The kinetic energy of the load KL and the ith body of the aircraft Ki are obtained through (10). The total kinetic energy of the system is then computed through K=KL+∑i=13Ki‍. Writing the result in the quadratic form K=1/2q˙TMqq˙, and by defining the inertia tensors IL=-∫VLρLSpLL2dVL‍ and Ii=-∫ViρiSpiCi2dVi‍, and taking into account the parallel axis theorem [28], yielding Ji≜-miSdCiB2+RCiBIiRCiBT and Di≜-miSdBL2+RBLJiRBLT, we have that the inertia matrix Mq∈R10×10 is given by(12)Mq=mL+mI3×3M12-RLIRBLSdmQ03×103×1∗M22M23WηTRBLRC2BI2ayWηTRBLRC3BI3ay∗∗QTJQQTRC2BI2ayQTRC3BI3ay∗∗∗ayTI2ay0∗∗∗∗ayTI3ay,with ∗ denoting terms that are deduced by symmetry, and (13)M12=-mRLISdBLWη-RLIRBLSdmRBLTWη,M22=WηTIL+D-SdBLRBLSdmRBLTM22=WηT2-RBLSdmRBLTSdBLWη,M23=WηT-SdBLRBLSdm+RBLJQ,where m≜∑i=13mi‍, J≜∑i=13Ji‍, D≜∑i=13Di‍, and dm≜∑i=13midCiB‍. Note that, from the inertia matrix (12), the four body dynamics are completely coupled, allowing to consider this interaction in the control law design and avoid the need of cascade control structures.The Coriolis and centripetal forces matrix,Cq,q˙∈R10×10, can be calculated via Christoffel symbols of the first kind [26]. The element from its kth row and jth column is computed through(14)Ckj=∑i=11012∂Mkj∂qi+∂Mki∂qj-∂Mij∂qkq˙i,where M is an element of the inertia matrix and k,j∈1,2,…,10.The forward kinematics of each body’s center of mass is obtained by makingpLL=03×1 in (5) and piCi=03×1 in (6). Thereafter, the potential energies of the load and of each body of the aircraft are obtained using (11). The total potential energy of the system is then computed by U=UL+∑i=13Ui‍=-g^TmL+mξ+mRLIdBL+RLIRBLdm.Then, the gravitational force vector is given by(15)gq=∂U∂q∈R10.


### 2.4. Generalized Forces

This section obtains the contributions to generalized forces of all nonconservative forces and torques that actuate on the tilt-rotor UAV with suspended load. Letf∈R3 and τ∈R3 denote a nonconservative force and a nonconservative torque, respectively, p∈R3 denote the point of application of f, and F be a reference frame rigidly attached to the body to which τ is applied. According to [29], the contributions of f and τ to the generalized forces can be computed through the following mappings:(16)ϑf=JpTfI∈Rn,(17)ϑτ=WFTτI∈Rn,where Jp≜∂p˙I/∂q˙∈R3×n and WF≜∂ωIFI/∂q˙∈R3×n.The thrust forces generated by the aircraft’s propellers, denoted byfR and fL, and the torques generated by the servomotors composing the tilting mechanisms, denoted by ταR and ταL, correspond to the input forces and torques of the system. Expressed in their respective thrusters’ frames (see Figure 2), they are given by fRC2=azfR, fLC3=azfL, ταRC2=ayταR, and ταLC3=ayταL, where az≜001T. In the inertial reference frame, these vectors are expressed as fRI=RC2IfRC2, fLI=RC3IfLC3, ταRI=RC2IταRC2, and ταLI=RC3IταLC3, with RCiI≜RLIRBLRCiB.This work assumes that the thrust forces are applied to the centers of mass of the respective thrusters’ groups, which correspond the origins ofC2 and C3. Making piCi=oCiCi=03×1 in (6) to obtain JoC2=∂o˙C2I/∂q˙ and JoC3=∂o˙C3I/∂q˙ yields(18)o˙CiI=I3×3RLISdBLT+RLIRBLSdCiBTRBLTWηRLIRBLSdCiBTQ03×103×1︸JoCiq˙.Then, using (16) and (18) yields(19)ϑfR=JoC2TfRI=RLIRBLRC2BazWηTSdBLRBLRC2Baz+WηTRBLSdC2BRC2BazQTSdC2BRC2Baz00fR,(20)ϑfL=JoC3TfLI=RLIRBLRC3BazWηTSdBLRBLRC3Baz+WηTRBLSdC3BRC3BazQTSdC3BRC3Baz00fL.The servomotors’ torques are applied to the respective thrusters’ bodies, and opposite torques due to reaction are also applied to the aircraft’s main body. These pairs of torques are mapped to generalized forces through (17). From the addition of angular velocities [26], we have(21)ωIBI=ωILI+ωLBI=RLIωILL+RLIRBLωLBBωIBI=03×3RLIWηRLIRBLQ03×103×1q˙,(22)ωIC2I=ωILI+ωLBI+ωBC2IωIC2I=03×3RLIWηRLIRBLQRLIRBLRC2Bay03×1q˙,(23)ωIC3I=ωILI+ωLBI+ωBC3IωIC3I=03×3RLIWηRLIRBLQ03×1RLIRBLRC3Bayq˙.Comparing (21), (22), and (23) to ωIBI=WBq˙, ωIC2I=WC2q˙, and ωIC3I=WC3q˙, respectively, and using (17) lead to(24)ϑταR=WC2TταRI+WBT-ταRI=03×1T03×1T02×1T10TταR,(25)ϑταL=WC3TταLI+WBT-ταLI=03×1T03×1T02×1T01TταL.This work also takes into account drag torques generated by the propellers. These are reaction torques applied to the thrusters’ bodies, due to the blades’ acceleration and drag [30]. Assuming steady-state for the angular velocity of the blades, the drag torques are given in the thrusters’ reference frames by τdrag,RC2=λRkτ/bfRC2, τdrag,LC3=λLkτ/bfLC3, where kτ and b are parameters obtained experimentally, and λR and λL are given according to the direction of rotation of the corresponding propeller: if counter-clockwise, 1; if clockwise, -1. In the inertial reference frame, we then have τdrag,RI=RC2Iτdrag,RC2 and τdrag,LI=RC3Iτdrag,LC3.As the drag torques are applied to the thrusters’ bodies, through (17), (22), and (23), we have(26)ϑτdrag,R=WC2Tτdrag,RI=λRkτb03×1WηTRBLRC2BazQTRC2Baz00fR,(27)ϑτdrag,L=WC3Tτdrag,LI=λLkτb03×1WηTRBLRC3BazQTRC3Baz00fL.Finally, the total mapping of the system inputs to generalized forces is obtained by summing up the contributions of the thrust forces, servomotor torques, and drag torques. Then, from (19)-(20) and (24)–(27),(28)ϑin=ϑfR+ϑfL+ϑταR+ϑταL+ϑτdrag,R+ϑτdrag,L=RLIRBLrRRLIRBLrL03×103×1WηTΛRrRWηTΛLrL03×103×1QTΓRrRQTΓLrL02×102×100100001fRfLταRταL≜Linqu,where(29)rR≜RC2Baz,ΛR≜SdBLRBL+RBLSdC2B+λRkτbRBL,ΓR≜SdC2B+λRkτbI3×3,rL≜RC3Baz,ΛL≜SdBLRBL+RBLSdC3B+λLkτbRBL,ΓL≜SdC3B+λLkτbI3×3.Note that, although in this current work no aerodynamic surfaces are considered in the tilt-rotor UAV (see Figure1), the system modeling developed here is general enough to describe the dynamics of any tilt-rotor carrying a suspended load. The only constraint is that the aircraft needs to be seen as a multibody system with similar frames definitions as those shown in Figure 2. Further, in presence of aerodynamics surfaces (e.g., wings and horizontal and vertical stabilizers), the lift and drag forces generated by them can be added to the model in a straightforward manner, by including these terms in (28), which will allow dealing with both helicopter and airplane flight modes.In addition, viscous friction is taken into account at the revolute joints of the tilting mechanisms and also at the point of connection between the rope and the tilt-rotor UAV. It is assumed that the friction torques are mapped to generalized forces as(30)ϑfr=-Lfrq˙,where Lfr≜diag0,0,0,0,0,0,μγ,μγ,μα,μα with μγ and μα being constant parameters.External disturbances applied to the suspended load are also considered, which may represent wind gusts affecting the system. Defining these disturbances in the inertial reference frame as the force vectord≜dxdydzT∈R3, and assuming that is applied to the load’s center of mass, it can be mapped to generalized forces through (16), yielding(31)ϑdb=∂o˙LI∂q˙Td=I3×303×303×203×103×1Td≜Ldbd,where o˙LI=ξ˙ is obtained by making pLL=03×1 in (7).


### 2.5. Equations of Motion

From (12), (14), and (15), the equations of motion of the tilt-rotor UAV with suspended load can be written in the Euler-Lagrange formulation as [27](32)Mqq¨+Cq,q˙q˙+gq=ϑ,where ϑ is the total generalized forces vector, obtained by summing up the contributions from the input forces and torques (28), viscous friction torques (30), and external disturbances (31). Thus, substituting ϑ=ϑin+ϑfr+ϑdb in (32) yields(33)Mqq¨+Cq,q˙+Lfrq˙+gq=Linqu+Ldbd.Finally, by defining the state vector(34)x≜qTq˙TT∈R20,and recalling the input vector defined in (28), the dynamic equations (33) can be written in the state-space representation(35)x˙=φx,u,d=q˙Mq-1-Cq,q˙+Lfrq˙-gq+Linqu+Ldbd,which is nonlinear and highly coupled. Since the load’s position and orientation are among the generalized coordinates (9), they are represented by the state variables (34), and, consequently, the load’s behavior is described explicitly by (35). On the other hand, the aircraft’s position and orientation are described only with respect to the load, thus appearing in (35) only implicitly.


## 2.1. System Description

The computer-aided design (CAD) model of the tilt-rotor UAV with suspended load is shown in Figure1. The system is regarded as a multibody mechanical system composed of four rigid bodies: (i) the aircraft’s main body, composed of Acrylonitrile Butadiene Styrene (ABS) structure, landing gear, batteries, and electronics; (ii) the right thruster group, composed of the right thruster and its tilting mechanism (a revolute joint); (iii) the left thruster group, composed of the left thruster and its tilting mechanism; and (iv) the suspended load group, comprising the load and the rope. The system is actuated through the aircraft’s thrusters and tilting mechanisms. For modeling purposes, the following assumptions are made:(A1)
The rope is rigid and has negligible mass.(A2)
The rope is connected to the aircraft’s geometric center.(A3)
The main body’s center of mass does not coincide with the aircraft’s geometric center.(A4)
The thrusters groups’ centers of mass are located at their respective tilting axes.Figure 1
The tilt-rotor UAV with suspended load (CAD model).Regarding the assumptions stated above, it is important to describe the main reasons of assumption (A3), which are as follows: (i) the center of mass is vertically displaced in order to improve the pitch moment; and (ii) the center of mass is displaced with respect to thex-axis of the geometric center, allowing obtaining non-null equilibrium points for the angular positions of the tilting mechanisms and pitch angle. This mechanical feature improves the controllability of the system in hover flight, allowing the projection of the thrusts generated by propellers on the x-axis without the need to tilt the thrusters’ group.


## 2.2. Kinematics from the Load’s Perspective

The load’s perspective approach consists in formulating the system’s kinematics regarding the suspended load as a free rigid body, with the aircraft as a multibody mechanical system rigidly coupled to it. Six reference frames are defined, as shown in Figure2: (i) the inertial reference frame, I; (ii) the aircraft’s geometric center frame, B; (iii) the main body center of mass frame, C1; (iv) the right thruster group center of mass frame, C2; (v) the left thruster group center of mass frame, C3; and (vi) the suspended load group center of mass frame, L. The load’s position with respect to the inertial frame I is denoted by ξ≜xyzT. The displacement vector from L to B corresponds to the rope and is expressed in L by dBL≜00lT, with l as the rope’s length. The displacement vectors from B to Ci are model parameters of the tilt-rotor UAV and are denoted by dCiB, expressed in B, with i∈1,2,3.Figure 2
Kinematic definitions, input forces, and torques.The orientation of the load with respect toI is parametrized by Euler angles, η≜ϕθψT, using the ZYX convention about local axes. The associated rotation matrix is given by (1)RLI≜Rz,ψRy,θRx,ϕ=cψcθcψsθsϕ-sψcϕcψsθcϕ+sψsϕsψcθsψsθsϕ+cψcϕsψsθcϕ-cψsϕ-sθcθsϕcθcϕ.On the other hand, the orientation of the aircraft’s geometric center frame with respect toL, corresponding to the orientation of the UAV with respect to the rope, is parametrized by two angles, γ≜γ1γ2T, such that (2)RBL≜Rx,-γ1Ry,-γ2=cγ20-sγ2sγ1sγ2cγ1sγ1cγ2cγ1sγ2-sγ1cγ1cγ2.The reference framesB and C1 are parallel to each other and attached to the same rigid body; thus the relative orientation is null, that is, RC1B=I3×3 (in this work In×n is an identity matrix with dimension n, 0n×m is n by m zeros matrix, and 1n×m is n by m ones matrix). Furthermore, the orientations of the thrusters’ groups with respect to B are described by(3)RC2B≜Rx,-βRy,αR=cαR0sαR-sβsαRcβsβcαR-cβsαR-sβcβcαR,RC3B≜Rx,βRy,αL=cαL0sαLsβsαLcβ-sβcαL-cβsαLsβcβcαL,where αR and αL are the tilting angles of the right and left thrusters, respectively, and β is a fixed inclination angle of the thrusters towards the aircraft geometric center, designed to improve the aircraft’s controllability [25]. The angular velocities of the system are given by ωILL=Wηη˙, ωLBB=Qγ˙, ωBC1C1=03×1, ωBC2C2=ayα˙R, and ωBC3C3=ayα˙L, where, for instance, ωILL denotes the angular velocity of L with respect to I, expressed in L, with(4)Wη≜10-sθ0cϕsϕcθ0-sϕcϕcθ,Q≜-cγ200-1sγ20,and ay≜010T.From the rigid transformations of the system, the forward kinematics of points that belong to each rigid body are given by(5)pLI=ξ+RLIpLL,(6)piI=ξ+RLIdBL+RLIRBLdCiB+RLIRBLRCiBpiCi,where pL is the position of a point that belongs to the suspended load body and pi belongs to the rigid body with attached frame Ci. Taking the time derivatives of (5) and (6) and making use of several properties of skew-symmetric matrices [26] yield(7)p˙LI=ξ˙+RLISpLLTωILL,(8)p˙iI=ξ˙+RLISdBLT+RLIRBLSdCiBTRBLTξ˙++RLIRBLRCiBSpiCiTRBLRCiBTωILLξ˙++RLIRBLSdCiBT+RLIRBLRCiBSpiCiTRCiBTωLBBξ˙++RLIRBLRCiBSpiCiTωBCiCi,where S· denotes an operator that maps a vector to a skew-symmetric matrix [26].The generalized coordinates of the system are chosen according to their degrees of freedom. For the equations of motion describing explicitly the time evolution of the load’s position and orientation, these are included in the generalized coordinates, which are chosen as(9)q≜ξTηTγTαRαLT∈R10.Note that, due to the chosen perspective for the system kinematics, the position and orientation of the aircraft with respect toI are not degrees of freedom of the system. Thus, they are not included in (9), and consequently their time evolution will not be described explicitly by the developed equations of motion.


## 2.3. Kinetic and Potential Energies

Initially, in order to derive the equations of motion through the Euler-Lagrange formulation, the kinetic and potential energies of each body of the mechanical system must be obtained. For eachith rigid body, these energies can be computed through the volume integrals [27](10)Ki=12∫Viρip˙iITp˙iIdVi,(11)Ui=-∫Viρig^TpiIdVi=-mig^ToCiI,respectively, where ρi corresponds to its density, Vi to its volume, and mi=∫ViρiVi‍ to its mass, g^≜00-g^zT denotes the gravitational acceleration vector expressed in I, and oCiI is the position vector obtained from the forward kinematics of the origin of Ci.The quadratic termsp˙LITp˙LI and p˙iITp˙iI,i∈1,2,3, are obtained using (7) and (8), respectively. The kinetic energy of the load KL and the ith body of the aircraft Ki are obtained through (10). The total kinetic energy of the system is then computed through K=KL+∑i=13Ki‍. Writing the result in the quadratic form K=1/2q˙TMqq˙, and by defining the inertia tensors IL=-∫VLρLSpLL2dVL‍ and Ii=-∫ViρiSpiCi2dVi‍, and taking into account the parallel axis theorem [28], yielding Ji≜-miSdCiB2+RCiBIiRCiBT and Di≜-miSdBL2+RBLJiRBLT, we have that the inertia matrix Mq∈R10×10 is given by(12)Mq=mL+mI3×3M12-RLIRBLSdmQ03×103×1∗M22M23WηTRBLRC2BI2ayWηTRBLRC3BI3ay∗∗QTJQQTRC2BI2ayQTRC3BI3ay∗∗∗ayTI2ay0∗∗∗∗ayTI3ay,with ∗ denoting terms that are deduced by symmetry, and (13)M12=-mRLISdBLWη-RLIRBLSdmRBLTWη,M22=WηTIL+D-SdBLRBLSdmRBLTM22=WηT2-RBLSdmRBLTSdBLWη,M23=WηT-SdBLRBLSdm+RBLJQ,where m≜∑i=13mi‍, J≜∑i=13Ji‍, D≜∑i=13Di‍, and dm≜∑i=13midCiB‍. Note that, from the inertia matrix (12), the four body dynamics are completely coupled, allowing to consider this interaction in the control law design and avoid the need of cascade control structures.The Coriolis and centripetal forces matrix,Cq,q˙∈R10×10, can be calculated via Christoffel symbols of the first kind [26]. The element from its kth row and jth column is computed through(14)Ckj=∑i=11012∂Mkj∂qi+∂Mki∂qj-∂Mij∂qkq˙i,where M is an element of the inertia matrix and k,j∈1,2,…,10.The forward kinematics of each body’s center of mass is obtained by makingpLL=03×1 in (5) and piCi=03×1 in (6). Thereafter, the potential energies of the load and of each body of the aircraft are obtained using (11). The total potential energy of the system is then computed by U=UL+∑i=13Ui‍=-g^TmL+mξ+mRLIdBL+RLIRBLdm.Then, the gravitational force vector is given by(15)gq=∂U∂q∈R10.


## 2.4. Generalized Forces

This section obtains the contributions to generalized forces of all nonconservative forces and torques that actuate on the tilt-rotor UAV with suspended load. Letf∈R3 and τ∈R3 denote a nonconservative force and a nonconservative torque, respectively, p∈R3 denote the point of application of f, and F be a reference frame rigidly attached to the body to which τ is applied. According to [29], the contributions of f and τ to the generalized forces can be computed through the following mappings:(16)ϑf=JpTfI∈Rn,(17)ϑτ=WFTτI∈Rn,where Jp≜∂p˙I/∂q˙∈R3×n and WF≜∂ωIFI/∂q˙∈R3×n.The thrust forces generated by the aircraft’s propellers, denoted byfR and fL, and the torques generated by the servomotors composing the tilting mechanisms, denoted by ταR and ταL, correspond to the input forces and torques of the system. Expressed in their respective thrusters’ frames (see Figure 2), they are given by fRC2=azfR, fLC3=azfL, ταRC2=ayταR, and ταLC3=ayταL, where az≜001T. In the inertial reference frame, these vectors are expressed as fRI=RC2IfRC2, fLI=RC3IfLC3, ταRI=RC2IταRC2, and ταLI=RC3IταLC3, with RCiI≜RLIRBLRCiB.This work assumes that the thrust forces are applied to the centers of mass of the respective thrusters’ groups, which correspond the origins ofC2 and C3. Making piCi=oCiCi=03×1 in (6) to obtain JoC2=∂o˙C2I/∂q˙ and JoC3=∂o˙C3I/∂q˙ yields(18)o˙CiI=I3×3RLISdBLT+RLIRBLSdCiBTRBLTWηRLIRBLSdCiBTQ03×103×1︸JoCiq˙.Then, using (16) and (18) yields(19)ϑfR=JoC2TfRI=RLIRBLRC2BazWηTSdBLRBLRC2Baz+WηTRBLSdC2BRC2BazQTSdC2BRC2Baz00fR,(20)ϑfL=JoC3TfLI=RLIRBLRC3BazWηTSdBLRBLRC3Baz+WηTRBLSdC3BRC3BazQTSdC3BRC3Baz00fL.The servomotors’ torques are applied to the respective thrusters’ bodies, and opposite torques due to reaction are also applied to the aircraft’s main body. These pairs of torques are mapped to generalized forces through (17). From the addition of angular velocities [26], we have(21)ωIBI=ωILI+ωLBI=RLIωILL+RLIRBLωLBBωIBI=03×3RLIWηRLIRBLQ03×103×1q˙,(22)ωIC2I=ωILI+ωLBI+ωBC2IωIC2I=03×3RLIWηRLIRBLQRLIRBLRC2Bay03×1q˙,(23)ωIC3I=ωILI+ωLBI+ωBC3IωIC3I=03×3RLIWηRLIRBLQ03×1RLIRBLRC3Bayq˙.Comparing (21), (22), and (23) to ωIBI=WBq˙, ωIC2I=WC2q˙, and ωIC3I=WC3q˙, respectively, and using (17) lead to(24)ϑταR=WC2TταRI+WBT-ταRI=03×1T03×1T02×1T10TταR,(25)ϑταL=WC3TταLI+WBT-ταLI=03×1T03×1T02×1T01TταL.This work also takes into account drag torques generated by the propellers. These are reaction torques applied to the thrusters’ bodies, due to the blades’ acceleration and drag [30]. Assuming steady-state for the angular velocity of the blades, the drag torques are given in the thrusters’ reference frames by τdrag,RC2=λRkτ/bfRC2, τdrag,LC3=λLkτ/bfLC3, where kτ and b are parameters obtained experimentally, and λR and λL are given according to the direction of rotation of the corresponding propeller: if counter-clockwise, 1; if clockwise, -1. In the inertial reference frame, we then have τdrag,RI=RC2Iτdrag,RC2 and τdrag,LI=RC3Iτdrag,LC3.As the drag torques are applied to the thrusters’ bodies, through (17), (22), and (23), we have(26)ϑτdrag,R=WC2Tτdrag,RI=λRkτb03×1WηTRBLRC2BazQTRC2Baz00fR,(27)ϑτdrag,L=WC3Tτdrag,LI=λLkτb03×1WηTRBLRC3BazQTRC3Baz00fL.Finally, the total mapping of the system inputs to generalized forces is obtained by summing up the contributions of the thrust forces, servomotor torques, and drag torques. Then, from (19)-(20) and (24)–(27),(28)ϑin=ϑfR+ϑfL+ϑταR+ϑταL+ϑτdrag,R+ϑτdrag,L=RLIRBLrRRLIRBLrL03×103×1WηTΛRrRWηTΛLrL03×103×1QTΓRrRQTΓLrL02×102×100100001fRfLταRταL≜Linqu,where(29)rR≜RC2Baz,ΛR≜SdBLRBL+RBLSdC2B+λRkτbRBL,ΓR≜SdC2B+λRkτbI3×3,rL≜RC3Baz,ΛL≜SdBLRBL+RBLSdC3B+λLkτbRBL,ΓL≜SdC3B+λLkτbI3×3.Note that, although in this current work no aerodynamic surfaces are considered in the tilt-rotor UAV (see Figure1), the system modeling developed here is general enough to describe the dynamics of any tilt-rotor carrying a suspended load. The only constraint is that the aircraft needs to be seen as a multibody system with similar frames definitions as those shown in Figure 2. Further, in presence of aerodynamics surfaces (e.g., wings and horizontal and vertical stabilizers), the lift and drag forces generated by them can be added to the model in a straightforward manner, by including these terms in (28), which will allow dealing with both helicopter and airplane flight modes.In addition, viscous friction is taken into account at the revolute joints of the tilting mechanisms and also at the point of connection between the rope and the tilt-rotor UAV. It is assumed that the friction torques are mapped to generalized forces as(30)ϑfr=-Lfrq˙,where Lfr≜diag0,0,0,0,0,0,μγ,μγ,μα,μα with μγ and μα being constant parameters.External disturbances applied to the suspended load are also considered, which may represent wind gusts affecting the system. Defining these disturbances in the inertial reference frame as the force vectord≜dxdydzT∈R3, and assuming that is applied to the load’s center of mass, it can be mapped to generalized forces through (16), yielding(31)ϑdb=∂o˙LI∂q˙Td=I3×303×303×203×103×1Td≜Ldbd,where o˙LI=ξ˙ is obtained by making pLL=03×1 in (7).


## 2.5. Equations of Motion

From (12), (14), and (15), the equations of motion of the tilt-rotor UAV with suspended load can be written in the Euler-Lagrange formulation as [27](32)Mqq¨+Cq,q˙q˙+gq=ϑ,where ϑ is the total generalized forces vector, obtained by summing up the contributions from the input forces and torques (28), viscous friction torques (30), and external disturbances (31). Thus, substituting ϑ=ϑin+ϑfr+ϑdb in (32) yields(33)Mqq¨+Cq,q˙+Lfrq˙+gq=Linqu+Ldbd.Finally, by defining the state vector(34)x≜qTq˙TT∈R20,and recalling the input vector defined in (28), the dynamic equations (33) can be written in the state-space representation(35)x˙=φx,u,d=q˙Mq-1-Cq,q˙+Lfrq˙-gq+Linqu+Ldbd,which is nonlinear and highly coupled. Since the load’s position and orientation are among the generalized coordinates (9), they are represented by the state variables (34), and, consequently, the load’s behavior is described explicitly by (35). On the other hand, the aircraft’s position and orientation are described only with respect to the load, thus appearing in (35) only implicitly.


## 3. Nonlinear State Estimation

This section presents a state estimation strategy to cope with the problem of predicting the states related to the load’s pose, using the information of noisy sensors with different sampling rates, in order to gather the information necessary to build the state vector (34). The following sensors are assumed to be available: (i) a Global Positioning System (GPS) equipment to measure the x and y positions of the UAV; (ii) a barometer to measure the UAV’s altitude; (iii) an Inertial Measurement Unit (IMU) that provides the orientation and angular velocity of the UAV; (iv) a camera that gives the position of the load with respect to the UAV; and (v) embedded sensors at the servomotors to provide the tilting angles and their time derivatives.The measurement model is highly nonlinear and therefore some well-known estimation techniques based on linearized transformations may have limited performance. Thereby, this paper considers the unscented Kalman filter strategy that estimates the mean and covariance of the likelihood distribution by means of the Unscented Transform (UT) and then uses the Kalman Filter (KF) equations to compute the posterior distribution [31, 32].


### 3.1. Measurement Equation

The relation between the variables measured by the embedded sensors and the state vectorx at time instant k is given by the measurement equation(36)yk=πxk+υk, where yk is the measured vector, πxk corresponds to a nonlinear mapping to be obtained, and υk is the measurement noise.LetξB≜xByBzBT be the position of the aircraft with respect to I (see Figure 3). Then, through forward kinematics, the following holds:(37)ξBξ,η=ξ+RLIdBL.Figure 3
Tilt-rotor UAV pose and load’s position measured by the available sensors.Additionally, letηB≜ϕBθBψBT be the aircraft attitude with respect to I parametrized by Euler angles using the local roll-pitch-yaw convention. Therefore, (38)RBI≜Rz,ψBRy,θBRx,ϕB=cψBcθBcψBsθBsϕB-sψBcϕBcψBsθBcϕB+sψBsϕBsψBcθBsψBsθBsϕB+cψBcϕBsψBsθBcϕB-cψBsϕB-sθBcθBsϕBcθBcϕB.Recalling that the relationRBI=RLIRBL holds and considering θB≠±π/2, it is possible to state that(39)ϕBη,γ=arctan⁡RLIRBL32RLIRBL33,θBη,γ=arcsin⁡-RLIRBL31,ψBη,γ=arctan⁡RLIRBL21RLIRBL11.Furthermore, the angular velocity provided by the IMU is given by(40)ωIBBη,γ,η˙,γ˙=ωILB+ωLBB=RBLTWηη˙+Qγ˙,where the matrices Wη and Q were defined in (4).LetdBLB denote the displacement vector from B to L, expressed in B, which is the measurement provided by the camera (see Figure 3). Then,(41)dBLBγ=-dLBB=-RBLTdLBL≜-RBLTdBL.Finally, (37), (39), (40), and (41) along with the system’s states αR, αL, α˙R, α˙L, and the measurement noise νk can be grouped into the nonlinear measurement equation(42)yk=πxk+νk≜ξBξ,ηϕBη,γθBη,γψBη,γωIBBη,γ,η˙,γ˙dBLBγαRαLα˙Rα˙L+νk=ξ+RLIdBLarctan⁡RLIRBL32RLIRBL33arcsin⁡-RLIRBL31arctan⁡RLIRBL21RLIRBL11RBLTWηη˙+Qγ˙-RBLTdBLαRαLα˙Rα˙L+νk.This work assumes synchronized sensors with sampling rates equal to10Ts, 2Ts, and Ts, respectively, for the GPS, camera, and remaining sensors, Ts being the control sampling time. Since not all information is available at time instant k, the dimension of the vector yk, as well as the transformation π·, will change with time. Thus, from now on, the variables that have different data acquisition rates will be denoted with the subscript ·Ts.


### 3.2. Unscented Kalman Filter

In the unscented Kalman filter (UKF), a fixed number of sigma points are chosen deterministically to capture the mean and covariance of the prior distribution. These points are then propagated through a nonlinear transformation to estimate the posterior distribution [32].Consider the discrete representation of (35) and the measurement equation (42)(43)xk=φxk-1,uk-1,dk-1+wk-1,(44)yk,Ts=πTsxk+νk,Ts,where wk is the process noise. Also, let ·^ denote estimated variables, (·)m∣n denote information at time instant m given measurements up to instant n, and E[·] be the expected value operator. This work assumes that measures of dk are always available, and the mean x^0∣0=E[x0] and covariance P0∣0x=E[(x0-x^0∣0)(x0-x^0∣0)T] are known. Furthermore, the process and measurement noises are assumed to have zero mean and covariances E[wkwkT]=Qk and E[νk,{Ts}νk,{Ts}T]=Rk,{Ts}, respectively, and the cross-covariance between them is assumed to be null, that is, E[wiνj,{Ts}T]=0∀i,j.In order to describe the prior statistics using UT, onlynsp=2n sigma points are necessary, where n is the dimension of the state vector x [31]. The sigma points X1,…,Xnsp must satisfy(45)x^=∑j=1nspγjXj,Px=∑j=1nspγjXj-x^Xj-x^T,where γj are weights defined as γj=1/2n subject to ∑j=1nspγj‍=1.The sigma point matrixX≜X1X2⋯Xnsp is chosen as(46)X=x^11×nsp+nPx1/2-Px1/2,where ·1/2 is the Cholesky square root.The likelihood distribution statistics can be obtained by propagating the sigma pointsX1,…,Xnsp through the nonlinear measurement equation (42), yielding(47)Yj,Ts=πTsXj,j=1,…,nspsuch that(48)y^Ts=∑j=1nspγjYj,Ts,PTsy=∑j=1nspγjYj,Ts-y^TsYj,Ts-y^TsT.Algorithm. The UKF algorithm is formed by a forecast and a data assimilation step [32]. The operations presented below must be performed at each measurement step k=1,2,3,….(i) Forecast Step
(1)
Form the sigma points using (46),(49)Xk-1∣k-1=x^k-1∣k-111×nsp+nPk-1∣k-1x1/2-Pk-1∣k-1x1/2.
(2)
Propagate the sigma points through the nonlinear transformation (43) (50)Xj,k∣k-1=φXj,k-1∣k-1,uk-1,dk-1,∀j=1,…,nsp.
(3)
Compute the predicted mean and covariance using (45) together with the process noise covariance (51)x^k∣k-1=∑j=1nspγjXj,k∣k-1,Pk∣k-1x=∑j=1nspγjXj,k∣k-1-x^k∣k-1Xj,k∣k-1-x^k∣k-1T+Qk-1.
(4)
Form the sigma point matrix with the estimation ofx^k∣k-1(52)Xk∣k-1=x^k∣k-111×nsp+nPk∣k-1x1/2-Pk∣k-1x1/2.
(5)
Propagate the sigma pointsXj,k∣k-1 through the measurement model (44)(53)Yj,k∣k-1,Ts=πTsXj,k∣k-1,∀j=1,…,nsp.
(6)
Compute the measurement statistics using (48) together with the measurement noise covariance and the cross-covariance of the states and measurements (54)y^k∣k-1,Ts=∑j=1nspγjYj,k∣k-1,Ts,Pk∣k-1,Tsy=∑j=1nspγjYj,k∣k-1,Ts-y^k∣k-1,TsYj,k∣k-1,Ts-y^k∣k-1,TsT+Rk,Ts,Pk∣k-1,Tsxy=∑j=1nspγjXj,k∣k-1-x^k∣k-1Yj,k∣k-1,Ts-y^k∣k-1,TsT.(ii) Data Assimilation Step
(1)
Compute the filter gainKk,Ts and the innovations υk,Ts(55)Kk,Ts=Pk∣k-1,TsxyPk∣k-1,Tsy-1,υk,Ts=yk,Ts-y^k∣k-1,Ts.
(2)
Compute the corrected mean and covariance conditional on the measurement information(56)x^k∣k=x^k∣k-1+Kk,Tsυk,Ts,Pk∣kx=Pk∣k-1x-Kk,TsPk∣k-1,TsyKk,TsT.Observe that the UKF algorithm uses at all time instantk all available information from the sensors to estimate the posterior distribution. For those time instants in which the measurement vector is full, that is, it has all sensors information, the estimation of x^k∣k will be more accurate and the covariance Pk∣kx reduced. On the other hand, when neither GPS nor camera information is available, the data assimilation step will have a less important role on the estimation algorithm, and x^k∣k will be less accurate and Pk∣kx bigger.Finally, the state vector estimated in this section,x^k∣k, will be used by the control strategy presented in the next section to perform the state feedback control.


## 3.1. Measurement Equation

The relation between the variables measured by the embedded sensors and the state vectorx at time instant k is given by the measurement equation(36)yk=πxk+υk, where yk is the measured vector, πxk corresponds to a nonlinear mapping to be obtained, and υk is the measurement noise.LetξB≜xByBzBT be the position of the aircraft with respect to I (see Figure 3). Then, through forward kinematics, the following holds:(37)ξBξ,η=ξ+RLIdBL.Figure 3
Tilt-rotor UAV pose and load’s position measured by the available sensors.Additionally, letηB≜ϕBθBψBT be the aircraft attitude with respect to I parametrized by Euler angles using the local roll-pitch-yaw convention. Therefore, (38)RBI≜Rz,ψBRy,θBRx,ϕB=cψBcθBcψBsθBsϕB-sψBcϕBcψBsθBcϕB+sψBsϕBsψBcθBsψBsθBsϕB+cψBcϕBsψBsθBcϕB-cψBsϕB-sθBcθBsϕBcθBcϕB.Recalling that the relationRBI=RLIRBL holds and considering θB≠±π/2, it is possible to state that(39)ϕBη,γ=arctan⁡RLIRBL32RLIRBL33,θBη,γ=arcsin⁡-RLIRBL31,ψBη,γ=arctan⁡RLIRBL21RLIRBL11.Furthermore, the angular velocity provided by the IMU is given by(40)ωIBBη,γ,η˙,γ˙=ωILB+ωLBB=RBLTWηη˙+Qγ˙,where the matrices Wη and Q were defined in (4).LetdBLB denote the displacement vector from B to L, expressed in B, which is the measurement provided by the camera (see Figure 3). Then,(41)dBLBγ=-dLBB=-RBLTdLBL≜-RBLTdBL.Finally, (37), (39), (40), and (41) along with the system’s states αR, αL, α˙R, α˙L, and the measurement noise νk can be grouped into the nonlinear measurement equation(42)yk=πxk+νk≜ξBξ,ηϕBη,γθBη,γψBη,γωIBBη,γ,η˙,γ˙dBLBγαRαLα˙Rα˙L+νk=ξ+RLIdBLarctan⁡RLIRBL32RLIRBL33arcsin⁡-RLIRBL31arctan⁡RLIRBL21RLIRBL11RBLTWηη˙+Qγ˙-RBLTdBLαRαLα˙Rα˙L+νk.This work assumes synchronized sensors with sampling rates equal to10Ts, 2Ts, and Ts, respectively, for the GPS, camera, and remaining sensors, Ts being the control sampling time. Since not all information is available at time instant k, the dimension of the vector yk, as well as the transformation π·, will change with time. Thus, from now on, the variables that have different data acquisition rates will be denoted with the subscript ·Ts.


## 3.2. Unscented Kalman Filter

In the unscented Kalman filter (UKF), a fixed number of sigma points are chosen deterministically to capture the mean and covariance of the prior distribution. These points are then propagated through a nonlinear transformation to estimate the posterior distribution [32].Consider the discrete representation of (35) and the measurement equation (42)(43)xk=φxk-1,uk-1,dk-1+wk-1,(44)yk,Ts=πTsxk+νk,Ts,where wk is the process noise. Also, let ·^ denote estimated variables, (·)m∣n denote information at time instant m given measurements up to instant n, and E[·] be the expected value operator. This work assumes that measures of dk are always available, and the mean x^0∣0=E[x0] and covariance P0∣0x=E[(x0-x^0∣0)(x0-x^0∣0)T] are known. Furthermore, the process and measurement noises are assumed to have zero mean and covariances E[wkwkT]=Qk and E[νk,{Ts}νk,{Ts}T]=Rk,{Ts}, respectively, and the cross-covariance between them is assumed to be null, that is, E[wiνj,{Ts}T]=0∀i,j.In order to describe the prior statistics using UT, onlynsp=2n sigma points are necessary, where n is the dimension of the state vector x [31]. The sigma points X1,…,Xnsp must satisfy(45)x^=∑j=1nspγjXj,Px=∑j=1nspγjXj-x^Xj-x^T,where γj are weights defined as γj=1/2n subject to ∑j=1nspγj‍=1.The sigma point matrixX≜X1X2⋯Xnsp is chosen as(46)X=x^11×nsp+nPx1/2-Px1/2,where ·1/2 is the Cholesky square root.The likelihood distribution statistics can be obtained by propagating the sigma pointsX1,…,Xnsp through the nonlinear measurement equation (42), yielding(47)Yj,Ts=πTsXj,j=1,…,nspsuch that(48)y^Ts=∑j=1nspγjYj,Ts,PTsy=∑j=1nspγjYj,Ts-y^TsYj,Ts-y^TsT.Algorithm. The UKF algorithm is formed by a forecast and a data assimilation step [32]. The operations presented below must be performed at each measurement step k=1,2,3,….(i) Forecast Step
(1)
Form the sigma points using (46),(49)Xk-1∣k-1=x^k-1∣k-111×nsp+nPk-1∣k-1x1/2-Pk-1∣k-1x1/2.
(2)
Propagate the sigma points through the nonlinear transformation (43) (50)Xj,k∣k-1=φXj,k-1∣k-1,uk-1,dk-1,∀j=1,…,nsp.
(3)
Compute the predicted mean and covariance using (45) together with the process noise covariance (51)x^k∣k-1=∑j=1nspγjXj,k∣k-1,Pk∣k-1x=∑j=1nspγjXj,k∣k-1-x^k∣k-1Xj,k∣k-1-x^k∣k-1T+Qk-1.
(4)
Form the sigma point matrix with the estimation ofx^k∣k-1(52)Xk∣k-1=x^k∣k-111×nsp+nPk∣k-1x1/2-Pk∣k-1x1/2.
(5)
Propagate the sigma pointsXj,k∣k-1 through the measurement model (44)(53)Yj,k∣k-1,Ts=πTsXj,k∣k-1,∀j=1,…,nsp.
(6)
Compute the measurement statistics using (48) together with the measurement noise covariance and the cross-covariance of the states and measurements (54)y^k∣k-1,Ts=∑j=1nspγjYj,k∣k-1,Ts,Pk∣k-1,Tsy=∑j=1nspγjYj,k∣k-1,Ts-y^k∣k-1,TsYj,k∣k-1,Ts-y^k∣k-1,TsT+Rk,Ts,Pk∣k-1,Tsxy=∑j=1nspγjXj,k∣k-1-x^k∣k-1Yj,k∣k-1,Ts-y^k∣k-1,TsT.(ii) Data Assimilation Step
(1)
Compute the filter gainKk,Ts and the innovations υk,Ts(55)Kk,Ts=Pk∣k-1,TsxyPk∣k-1,Tsy-1,υk,Ts=yk,Ts-y^k∣k-1,Ts.
(2)
Compute the corrected mean and covariance conditional on the measurement information(56)x^k∣k=x^k∣k-1+Kk,Tsυk,Ts,Pk∣kx=Pk∣k-1x-Kk,TsPk∣k-1,TsyKk,TsT.Observe that the UKF algorithm uses at all time instantk all available information from the sensors to estimate the posterior distribution. For those time instants in which the measurement vector is full, that is, it has all sensors information, the estimation of x^k∣k will be more accurate and the covariance Pk∣kx reduced. On the other hand, when neither GPS nor camera information is available, the data assimilation step will have a less important role on the estimation algorithm, and x^k∣k will be less accurate and Pk∣kx bigger.Finally, the state vector estimated in this section,x^k∣k, will be used by the control strategy presented in the next section to perform the state feedback control.


## 4. Suspended Load Path Tracking Controller

This section describes the model predictive controller design for performing the suspended load path tracking. The main objectives of the control system are as follows: ensure closed-loop stability, reject constant external disturbances and parametric uncertainties, and satisfy constraints on state deviations and control inputs.Aiming at an improved path tracking control and recalling the system’s underactuated behavior, the load’s positionξ=xyzT and its yaw angle ψ are chosen to be regulated, while the other degrees of freedom will be only stabilized.This paper works with the incremental MPC framework [33], for which the prediction model is obtained using a discrete linear time-varying state-space model that is not affine in the parameters. On one hand, the use of this kind of linear system improves the aircraft’s nonlinear dynamics representation, when comparing to time-invariant linearized models, and, particularly for the proposed application, allows yaw movements regulation and copes with the rope’s length variation during the take-off and landing. On the other hand, because the time-varying linearized model used in this work is non-affine in the time-varying parameters, it is not possible to obtain, for instance, a polytopic representation for the linearized model and, then, use some well-know techniques to improve the MPC robustness.


### 4.1. Linearized Error Dynamics

Seeking to obtain the discrete linear time-varying state-space model to construct the prediction model, the equations of motion (35) must be linearized around a time-varying trajectory. Additionally, due to limited computational resources, this process needs to be done with most of the physical parameters numerically evaluated. However, it is possible to let some physical parameters be variables in a way that they will appear in the linearized Jacobians after finishing the linearization process.Letxtr and utr denote trajectory values and, for linearization purposes, consider d=03×1 (see (31)). This work assumes that the desired trajectory is feasible, that is,(57)x˙tr=φxtr,utr,d.Then, linearizing the state-space equations (35) around these trajectories, through first-order expansion in Taylor series, yields(58)Δx˙=AtΔx+BtΔu,where Δx≜x-xtr, Δu≜u-utr, and(59)At=∂φx,u,d∂xx=xtru=utr∈R20×20,Bt=∂φx,u,d∂ux=xtru=utr∈R20×4.In this work, the trajectory values forx and u are given by(60)xtr=qtrTq˙trTT,(61)utr=Linqtr+Mqtrq¨tr+Cqtr,q˙tr+Lfrq˙tr+gqtr,where Linqtr+ denotes the left pseudoinverse of Linqtr and qtr, q˙tr, and q¨tr are provided reference signals with qtr≜xtrtytrtztrtϕeqθeqψtrtγ1eqγ2eqαReqαLeqT, in which ·eq is the state’s equilibrium value. Notice that utr, since it is computed using a left pseudoinverse, will be an exact solution to the dynamic equations (33) only if the desired trajectory is feasible.Therefore, by linearizing the system using (59) with the trajectories defined in (60) and (61) added to the rope’s length lt as a time-varying parameter, the linearized Jacobians are Aζt and Bζt, where ζ(t)≜[xtr(t)ytr(t)ztr(t)ψtr(t)x˙tr(t)y˙tr(t)z˙tr(t)ψ˙tr(t)x¨tr(t)y¨tr(t)z¨tr(t)ψ¨tr(t)l(t)]T is the vector of time-varying parameters.To improve the trajectory tracking of the regulated variables and provide constant disturbance and parametric uncertainties rejection, the state vectorΔx is augmented with integral actions [34], yielding(62)Δx¯≜Δx∫ξ-ξtr∫ψ-ψtr∈R24,whose dynamics are given by(63)Δx¯˙=Aζt020×410000001000000100000000104×1404×4︸A¯ζtΔx¯+Bζt04×4︸B¯ζtΔu.To obtain a discrete prediction model using the incremental form and, thereafter, improve performance with input integrators added to the closed-loop system [33], it will be necessary to map the model (63) from the continuous-time to the discrete-time domain, which yields (64)Δx¯k+1=A¯ζkΔx¯k+B¯ζkΔuk,being the matrices A¯ζk and B¯ζk obtained after discretizing the model using a zero-order hold with sampling time Ts (in this work, variables in continuous-time and discrete-time domain are differentiated by the time variable t and the sampling variable k).Finally, choosing the control incrementδuk≜Δuk-Δuk-1 to be the control input, the extended discrete linearized system can be rewritten in the incremental form as(65)Δx¯k+1Δuk︸ΔXk+1=A¯ζkB¯ζk04×24I4×4︸AζkΔx¯kΔuk-1︸ΔXk+B¯ζkI4×4︸Bζkδuk.


### 4.2. Prediction Model

The state-space model (65) gives the one-step ahead prediction and can be used recursively to obtain the prediction model considering a prediction horizon Np and control horizon Nc. Thus, considering the case where Nc<Np and assuming δuk+i=04×1,∀i≥Nc, the Np-step ahead prediction yields(66)ΔXk+Np=∏l=1NpAζk+Np-lΔXkΔXk+Np+∏l=1Np-1Aζk+Np-lBζkδuk+⋯ΔXk+Np+∏l=1Np-NcAζk+Np-lBζk+Nc-1δuk+Nc-1.Thereby, definingΔX→≜ΔXk+1T⋯ΔXk+NpTT and δu→≜δukT⋯δuk+Nc-1TT, it is possible to write the prediction model as(67)ΔX→=PΔXk+Hδu→,where the matrices P∈R28·Np×28 and H∈R28·Np×4·Nc are given by (68)P=∏l=11Aζk+1-lT∏l=12Aζk+2-lT⋯∏l=1NpAζk+Np-lTT,H=Bζk028×4⋯028×4∏l=12-1Aζk+2-lBζkBζk+1⋯028×4⋮⋮⋱⋮∏l=1Np-1Aζk+Np-lBζk∏l=1Np-2Aζk+Np-lBζk+1⋯∏l=1Np-NcAζk+Np-lBζk+Nc-1.


### 4.3. Optimization Problem

Consider the standard quadratic cost function(69)J=∑i=1NpWk+i-ΔXk+iQ2+∑j=0Nc-1δuk+jR2,where Wk+i=Xk+itr-Xktr for i=1,…,Np, Xk+itr being the reference trajectory at the instant k+i. Moreover, Q∈R28×28 and R∈R4×4 are, respectively, weighting matrices of states error and control effort.The cost function (69) can be written in the matrix form by means of the prediction model (67) as(70)J=Hδu→+PΔXk-W→TΩQHδu→+PΔXk-W→+δuT→ΩRδu→,where ΩQ≜blkdiagQ,…,Q∈R28·Np×28·Np and ΩR≜blkdiagR,…,R∈R4·Nc×4·Nc are block diagonal matrices, and W→≜Wk+1T⋯Wk+NpTT.Finally, (70) can be rewritten in the canonical quadratic form as(71)J=12δuT→Λδu→+fTδu→+f0,where f0=PΔXk-W→TΩQPΔXk-W→, fT=2PΔX-W→TΩQH, and Λ=2HTΩQH+ΩR.Adding constraints on the objective variableδu→, the most general optimization problem must be solved:(72)minδu→12δuT→Λδu→+fTδu→+f0s.t.Mδu→≤N.The constraints considered above can be used to limit the control signal amplitude avoiding saturation in the actuators and to limit the maximum state error. In both cases, the constraints must be mapped to the amplitude of the control incrementδu as in [14](73)M≜C2-C2H-H,N≜C1umax-Δuk-1-uktr→-C1umin-Δuk-1+uktr→ΔXmax-PΔXk-ΔXmin+PΔXk,where (·)max and (·)min are, respectively, the maximum and minimum values allowed for the variable, and the matrices C1 and C2 are defined as (74)C1=I4×4I4×4⋮I4×4,C2=I4×404×4⋯04×4I4×4I4×4⋯04×4⋮⋮⋱⋮I4×4I4×4⋯I4×4.


### 4.4. Applied Control Signal

Using the relationsΔuk=uk-uktr and δuk=Δuk-Δuk-1, the control signal at the time instant k can be written as(75)uk=uktr+Δuk-1+δuk,where δuk results from the optimization problem (72) and uktr is the feedforward term (61).


## 4.1. Linearized Error Dynamics

Seeking to obtain the discrete linear time-varying state-space model to construct the prediction model, the equations of motion (35) must be linearized around a time-varying trajectory. Additionally, due to limited computational resources, this process needs to be done with most of the physical parameters numerically evaluated. However, it is possible to let some physical parameters be variables in a way that they will appear in the linearized Jacobians after finishing the linearization process.Letxtr and utr denote trajectory values and, for linearization purposes, consider d=03×1 (see (31)). This work assumes that the desired trajectory is feasible, that is,(57)x˙tr=φxtr,utr,d.Then, linearizing the state-space equations (35) around these trajectories, through first-order expansion in Taylor series, yields(58)Δx˙=AtΔx+BtΔu,where Δx≜x-xtr, Δu≜u-utr, and(59)At=∂φx,u,d∂xx=xtru=utr∈R20×20,Bt=∂φx,u,d∂ux=xtru=utr∈R20×4.In this work, the trajectory values forx and u are given by(60)xtr=qtrTq˙trTT,(61)utr=Linqtr+Mqtrq¨tr+Cqtr,q˙tr+Lfrq˙tr+gqtr,where Linqtr+ denotes the left pseudoinverse of Linqtr and qtr, q˙tr, and q¨tr are provided reference signals with qtr≜xtrtytrtztrtϕeqθeqψtrtγ1eqγ2eqαReqαLeqT, in which ·eq is the state’s equilibrium value. Notice that utr, since it is computed using a left pseudoinverse, will be an exact solution to the dynamic equations (33) only if the desired trajectory is feasible.Therefore, by linearizing the system using (59) with the trajectories defined in (60) and (61) added to the rope’s length lt as a time-varying parameter, the linearized Jacobians are Aζt and Bζt, where ζ(t)≜[xtr(t)ytr(t)ztr(t)ψtr(t)x˙tr(t)y˙tr(t)z˙tr(t)ψ˙tr(t)x¨tr(t)y¨tr(t)z¨tr(t)ψ¨tr(t)l(t)]T is the vector of time-varying parameters.To improve the trajectory tracking of the regulated variables and provide constant disturbance and parametric uncertainties rejection, the state vectorΔx is augmented with integral actions [34], yielding(62)Δx¯≜Δx∫ξ-ξtr∫ψ-ψtr∈R24,whose dynamics are given by(63)Δx¯˙=Aζt020×410000001000000100000000104×1404×4︸A¯ζtΔx¯+Bζt04×4︸B¯ζtΔu.To obtain a discrete prediction model using the incremental form and, thereafter, improve performance with input integrators added to the closed-loop system [33], it will be necessary to map the model (63) from the continuous-time to the discrete-time domain, which yields (64)Δx¯k+1=A¯ζkΔx¯k+B¯ζkΔuk,being the matrices A¯ζk and B¯ζk obtained after discretizing the model using a zero-order hold with sampling time Ts (in this work, variables in continuous-time and discrete-time domain are differentiated by the time variable t and the sampling variable k).Finally, choosing the control incrementδuk≜Δuk-Δuk-1 to be the control input, the extended discrete linearized system can be rewritten in the incremental form as(65)Δx¯k+1Δuk︸ΔXk+1=A¯ζkB¯ζk04×24I4×4︸AζkΔx¯kΔuk-1︸ΔXk+B¯ζkI4×4︸Bζkδuk.


## 4.2. Prediction Model

The state-space model (65) gives the one-step ahead prediction and can be used recursively to obtain the prediction model considering a prediction horizon Np and control horizon Nc. Thus, considering the case where Nc<Np and assuming δuk+i=04×1,∀i≥Nc, the Np-step ahead prediction yields(66)ΔXk+Np=∏l=1NpAζk+Np-lΔXkΔXk+Np+∏l=1Np-1Aζk+Np-lBζkδuk+⋯ΔXk+Np+∏l=1Np-NcAζk+Np-lBζk+Nc-1δuk+Nc-1.Thereby, definingΔX→≜ΔXk+1T⋯ΔXk+NpTT and δu→≜δukT⋯δuk+Nc-1TT, it is possible to write the prediction model as(67)ΔX→=PΔXk+Hδu→,where the matrices P∈R28·Np×28 and H∈R28·Np×4·Nc are given by (68)P=∏l=11Aζk+1-lT∏l=12Aζk+2-lT⋯∏l=1NpAζk+Np-lTT,H=Bζk028×4⋯028×4∏l=12-1Aζk+2-lBζkBζk+1⋯028×4⋮⋮⋱⋮∏l=1Np-1Aζk+Np-lBζk∏l=1Np-2Aζk+Np-lBζk+1⋯∏l=1Np-NcAζk+Np-lBζk+Nc-1.


## 4.3. Optimization Problem

Consider the standard quadratic cost function(69)J=∑i=1NpWk+i-ΔXk+iQ2+∑j=0Nc-1δuk+jR2,where Wk+i=Xk+itr-Xktr for i=1,…,Np, Xk+itr being the reference trajectory at the instant k+i. Moreover, Q∈R28×28 and R∈R4×4 are, respectively, weighting matrices of states error and control effort.The cost function (69) can be written in the matrix form by means of the prediction model (67) as(70)J=Hδu→+PΔXk-W→TΩQHδu→+PΔXk-W→+δuT→ΩRδu→,where ΩQ≜blkdiagQ,…,Q∈R28·Np×28·Np and ΩR≜blkdiagR,…,R∈R4·Nc×4·Nc are block diagonal matrices, and W→≜Wk+1T⋯Wk+NpTT.Finally, (70) can be rewritten in the canonical quadratic form as(71)J=12δuT→Λδu→+fTδu→+f0,where f0=PΔXk-W→TΩQPΔXk-W→, fT=2PΔX-W→TΩQH, and Λ=2HTΩQH+ΩR.Adding constraints on the objective variableδu→, the most general optimization problem must be solved:(72)minδu→12δuT→Λδu→+fTδu→+f0s.t.Mδu→≤N.The constraints considered above can be used to limit the control signal amplitude avoiding saturation in the actuators and to limit the maximum state error. In both cases, the constraints must be mapped to the amplitude of the control incrementδu as in [14](73)M≜C2-C2H-H,N≜C1umax-Δuk-1-uktr→-C1umin-Δuk-1+uktr→ΔXmax-PΔXk-ΔXmin+PΔXk,where (·)max and (·)min are, respectively, the maximum and minimum values allowed for the variable, and the matrices C1 and C2 are defined as (74)C1=I4×4I4×4⋮I4×4,C2=I4×404×4⋯04×4I4×4I4×4⋯04×4⋮⋮⋱⋮I4×4I4×4⋯I4×4.


## 4.4. Applied Control Signal

Using the relationsΔuk=uk-uktr and δuk=Δuk-Δuk-1, the control signal at the time instant k can be written as(75)uk=uktr+Δuk-1+δuk,where δuk results from the optimization problem (72) and uktr is the feedforward term (61).


## 5. Simulation Scenario

This section describes the simulation scenario used to evaluate the performance of the proposed control and state estimation strategy. The implemented control structure is presented in Figure4.Figure 4
The proposed control and state estimation strategy.


### 5.1. Model and Design Parameters

The model parameters of the tilt-rotor UAV with suspended load are shown in Table1. Mass, inertia, and displacement parameters of the aircraft were obtained from CAD model, designed in Solidworks® software. The parameters related to the suspended load, as well as kτ and b, are the same considered in [12]. The gravitational acceleration is assumed constant, and λR and λL are given according to the following: the right propeller of the aircraft rotates counter-clockwise and the left one rotates clockwise. From those data, and assuming ψ=0, the following equilibrium points was obtained to construct the vector qtr used in (60) and (61)(76)ϕeq=0,θeq=0,γ1eq=0.0001563,γ2eq=0.0287134,αReq=0.0288375,αLeq=0.0283718,fReq=9.7903455,fLeq=9.8252665,ταReq=0,ταLeq=0.Table 1
Model parameters of the tilt-rotor UAV with suspended load.
Parameter
Value
m
L
0.05000 Kg
m1
1.70249 Kg
m
2, m3
0.13973 Kg
lmax
1 m
d
B
L
0
0
l
T m
d
C
1
B
-
0.00433
0.00060
-
0.04559
T m
d
C
2
B
0.00002
-
0.27761
0.05493
T m
d
C
3
B
0.00077
0.27761
0.05493
T m
I
L
2.645
⋅
1
0
-
6
⋅
I
3
×
3 Kg⋅m2
I
1
3697.66749
0.36342
-
9.51029
∗
840.10403
0.61804
∗
∗
3865.05354
⋅
1
0
-
6 Kg⋅m2
I
2
441.68245
0
0
∗
441.67985
-
1.07006
∗
∗
0.64418
⋅
1
0
-
6 Kg⋅m2
I
3
441.68245
0
0
∗
441.67985
1.07006
∗
∗
0.64418
⋅
1
0
-
6 Kg⋅m2
g^
0
0
-
9.81
T m/s2
kτ
1.7
⋅
1
0
-
7 N⋅m⋅s2
b
9.5
⋅
1
0
-
6 N⋅s2
(λR,λL)
(1,-1)
β
5
°
μ
γ, μα
0.005N⋅m/(rad/s)For simplicity, the sensors’ noise was assumed to have a Gaussian probability distribution, and the measurement error is defined as three times the standard deviation. However, such assumption is not required by the UKF to estimate the system states (34). Additionally, this work considers known disturbances affecting the aircraft, which may represent the presence of wind gusts, and all sensors are synchronized. Table 2 shows the sensors’ measurement error and their sampling time. The relations between the measured variables and those presented in (42) are as follows: ξB≡s1,s2,s3, ηB≡s4,s5,s6, ωIBB≡s7,s8,s9, dBLB≡s10,s11,s12, αR≡s13, αL≡s14, α˙R≡s15, and α˙L≡s16.Table 2
Sensors parameters for the considered scenario.
Sensor
Measurement error
Sampling time
GPS
{
s
1
,
s
2
}
±0.15 m
120 ms
Barometer
{
s
3
}
±0.51 m
12 ms
IMU
{
s
4
,
s
5
,
s
6
}
±2.618⋅10-3 rad
12 ms
{
s
7
,
s
8
,
s
9
}
±16.558⋅10-3 rad/s
Camera
{
s
10
,
s
11
}
±0.005 m
24 ms
{
s
12
}
±0.02 m
Servos
{
s
13
,
s
14
}
±5.67⋅10-3 rad
12 ms
{
s
15
,
s
16
}
±0.50772 rad/sFor time instantsk multiples of 120 ms, the measurement vector is full, being yk,{Ts} = [s1s2s3s4s5s6s7s8s9s10s11s12s13s14s15s16]T, while for those multiples of 24 ms the vector does not have GPS information, and then yk,Ts=s3s4s5s6s7s8s9s10s11s12s13s14s15s16T. Finally, for time instants multiples of 12 ms, the measurement vector does not have neither GPS nor camera information; thus yk,Ts=s3s4s5s6s7s8s9s13s14s15s16T.The initial state vector is considered to be precisely known, that is,x^0∣0=x0, and the covariance matrices used in the filter algorithm are given by (77)P0∣0x=0.001·I20×20,Qk=diag0.0001·12×1,0.001,0.00001·13×1,0.05·12×1,mmm0.001·12×1,0.01·16×1,0.05·12×1,0.1·12×1,and Rk,Ts being a diagonal matrix composed by the sensor standard deviation, that is, the measurement error divided by three, regarding the sensors information available at the time instant k.For discretization of the linearized Jacobians,A¯ζt and B¯ζt, in (63), it is assumed that Ts = 12 ms. The prediction and control horizons, chosen considering the trade-off between good performance and small computational cost, are Np=100 and Nc=10. Furthermore, the saturation level of the tilt-rotor UAV actuators and the maximum state error allowed, used in (73), are(78)Δx=-1,1,Δy=-1,1,Δz=-1,1,Δϕ=-0.8,0.8,Δθ=-0.8,0.8,Δψ=-0.8,0.8,Δγ1=-0.8,0.8,Δγ2=-0.8,0.8,ΔαR=-0.8,0.8,ΔαL=-0.8,0.8,fR=0,15,fL=0,15,ταR=-2,2,ταL=-2,2,where the error limitations were chosen regarding measurement error and disturbances effect on the system. The actuators’ bounds are due to physical constraints.The Bryson’s rule [35] was used as starting point to synthesize the MPC controller’s weighting matrices, which are given by(79)Q=diag4022,4022,2022,5π/22,5π/22,10π2,10π/22,10π/22,0.1π/22,0.1π/22,1022,1022,522,1π/32,1π/32,1π/42,53π2,Q=diag253π2,0.13π2,0.13π2,40,40,40,20,40fReq-152,40fLeq-152,20ταReq-22,20ταLeq-22,R=diag200fReq-152,200fLeq-152,1000ταReq-22,1000ταLeq-22.


### 5.2. Desired Trajectory

To explore the controller capabilities, the proposed trajectory to be tracked by the suspended load is composed of several interconnected paths, which are described by polynomial and/or sinusoidal functions. It starts with vertical take-off, followed by straight line tracking with changes in direction, and ends with vertical landing. Along with the desired position, yaw’s movements are specified in order to always have the aircraft performing the trajectory head-on (see Figure5). Moreover, to evaluate the disturbance compensation of the proposed strategy, external forces are applied to the suspended load. Figure 6 shows the disturbance profile for the desired trajectory, which may represent sustained wind gusts affecting the load. The magnitude of the disturbances may seem low at a first glance; however the mass of the load is very small (see Table 1).Figure 5
Three-dimensional trajectories of the suspended load (blue), the tilt-rotor UAV (red), and the desired trajectory (dashed black). At the beginning and at the end of the path, the rope’s length is reduced, while throughout the trajectory the rope reaches its maximum value.Figure 6
Profile of the disturbance forces applied to the load and expressed in the inertial frameI.


## 5.1. Model and Design Parameters

The model parameters of the tilt-rotor UAV with suspended load are shown in Table1. Mass, inertia, and displacement parameters of the aircraft were obtained from CAD model, designed in Solidworks® software. The parameters related to the suspended load, as well as kτ and b, are the same considered in [12]. The gravitational acceleration is assumed constant, and λR and λL are given according to the following: the right propeller of the aircraft rotates counter-clockwise and the left one rotates clockwise. From those data, and assuming ψ=0, the following equilibrium points was obtained to construct the vector qtr used in (60) and (61)(76)ϕeq=0,θeq=0,γ1eq=0.0001563,γ2eq=0.0287134,αReq=0.0288375,αLeq=0.0283718,fReq=9.7903455,fLeq=9.8252665,ταReq=0,ταLeq=0.Table 1
Model parameters of the tilt-rotor UAV with suspended load.
Parameter
Value
m
L
0.05000 Kg
m1
1.70249 Kg
m
2, m3
0.13973 Kg
lmax
1 m
d
B
L
0
0
l
T m
d
C
1
B
-
0.00433
0.00060
-
0.04559
T m
d
C
2
B
0.00002
-
0.27761
0.05493
T m
d
C
3
B
0.00077
0.27761
0.05493
T m
I
L
2.645
⋅
1
0
-
6
⋅
I
3
×
3 Kg⋅m2
I
1
3697.66749
0.36342
-
9.51029
∗
840.10403
0.61804
∗
∗
3865.05354
⋅
1
0
-
6 Kg⋅m2
I
2
441.68245
0
0
∗
441.67985
-
1.07006
∗
∗
0.64418
⋅
1
0
-
6 Kg⋅m2
I
3
441.68245
0
0
∗
441.67985
1.07006
∗
∗
0.64418
⋅
1
0
-
6 Kg⋅m2
g^
0
0
-
9.81
T m/s2
kτ
1.7
⋅
1
0
-
7 N⋅m⋅s2
b
9.5
⋅
1
0
-
6 N⋅s2
(λR,λL)
(1,-1)
β
5
°
μ
γ, μα
0.005N⋅m/(rad/s)For simplicity, the sensors’ noise was assumed to have a Gaussian probability distribution, and the measurement error is defined as three times the standard deviation. However, such assumption is not required by the UKF to estimate the system states (34). Additionally, this work considers known disturbances affecting the aircraft, which may represent the presence of wind gusts, and all sensors are synchronized. Table 2 shows the sensors’ measurement error and their sampling time. The relations between the measured variables and those presented in (42) are as follows: ξB≡s1,s2,s3, ηB≡s4,s5,s6, ωIBB≡s7,s8,s9, dBLB≡s10,s11,s12, αR≡s13, αL≡s14, α˙R≡s15, and α˙L≡s16.Table 2
Sensors parameters for the considered scenario.
Sensor
Measurement error
Sampling time
GPS
{
s
1
,
s
2
}
±0.15 m
120 ms
Barometer
{
s
3
}
±0.51 m
12 ms
IMU
{
s
4
,
s
5
,
s
6
}
±2.618⋅10-3 rad
12 ms
{
s
7
,
s
8
,
s
9
}
±16.558⋅10-3 rad/s
Camera
{
s
10
,
s
11
}
±0.005 m
24 ms
{
s
12
}
±0.02 m
Servos
{
s
13
,
s
14
}
±5.67⋅10-3 rad
12 ms
{
s
15
,
s
16
}
±0.50772 rad/sFor time instantsk multiples of 120 ms, the measurement vector is full, being yk,{Ts} = [s1s2s3s4s5s6s7s8s9s10s11s12s13s14s15s16]T, while for those multiples of 24 ms the vector does not have GPS information, and then yk,Ts=s3s4s5s6s7s8s9s10s11s12s13s14s15s16T. Finally, for time instants multiples of 12 ms, the measurement vector does not have neither GPS nor camera information; thus yk,Ts=s3s4s5s6s7s8s9s13s14s15s16T.The initial state vector is considered to be precisely known, that is,x^0∣0=x0, and the covariance matrices used in the filter algorithm are given by (77)P0∣0x=0.001·I20×20,Qk=diag0.0001·12×1,0.001,0.00001·13×1,0.05·12×1,mmm0.001·12×1,0.01·16×1,0.05·12×1,0.1·12×1,and Rk,Ts being a diagonal matrix composed by the sensor standard deviation, that is, the measurement error divided by three, regarding the sensors information available at the time instant k.For discretization of the linearized Jacobians,A¯ζt and B¯ζt, in (63), it is assumed that Ts = 12 ms. The prediction and control horizons, chosen considering the trade-off between good performance and small computational cost, are Np=100 and Nc=10. Furthermore, the saturation level of the tilt-rotor UAV actuators and the maximum state error allowed, used in (73), are(78)Δx=-1,1,Δy=-1,1,Δz=-1,1,Δϕ=-0.8,0.8,Δθ=-0.8,0.8,Δψ=-0.8,0.8,Δγ1=-0.8,0.8,Δγ2=-0.8,0.8,ΔαR=-0.8,0.8,ΔαL=-0.8,0.8,fR=0,15,fL=0,15,ταR=-2,2,ταL=-2,2,where the error limitations were chosen regarding measurement error and disturbances effect on the system. The actuators’ bounds are due to physical constraints.The Bryson’s rule [35] was used as starting point to synthesize the MPC controller’s weighting matrices, which are given by(79)Q=diag4022,4022,2022,5π/22,5π/22,10π2,10π/22,10π/22,0.1π/22,0.1π/22,1022,1022,522,1π/32,1π/32,1π/42,53π2,Q=diag253π2,0.13π2,0.13π2,40,40,40,20,40fReq-152,40fLeq-152,20ταReq-22,20ταLeq-22,R=diag200fReq-152,200fLeq-152,1000ταReq-22,1000ταLeq-22.


## 5.2. Desired Trajectory

To explore the controller capabilities, the proposed trajectory to be tracked by the suspended load is composed of several interconnected paths, which are described by polynomial and/or sinusoidal functions. It starts with vertical take-off, followed by straight line tracking with changes in direction, and ends with vertical landing. Along with the desired position, yaw’s movements are specified in order to always have the aircraft performing the trajectory head-on (see Figure5). Moreover, to evaluate the disturbance compensation of the proposed strategy, external forces are applied to the suspended load. Figure 6 shows the disturbance profile for the desired trajectory, which may represent sustained wind gusts affecting the load. The magnitude of the disturbances may seem low at a first glance; however the mass of the load is very small (see Table 1).Figure 5
Three-dimensional trajectories of the suspended load (blue), the tilt-rotor UAV (red), and the desired trajectory (dashed black). At the beginning and at the end of the path, the rope’s length is reduced, while throughout the trajectory the rope reaches its maximum value.Figure 6
Profile of the disturbance forces applied to the load and expressed in the inertial frameI.


## 6. Numerical Results and Discussion

This section presents the numerical results obtained with the proposed control and state estimation strategy when subjected to the simulation scenario described before. The simulations have been carried out using the MATLAB/Simulink® environment. A detailed analysis about the performance of the control system when solving the path tracking problem of a suspended load carried by a tilt-rotor UAV is provided.The trajectories performed by the tilt-rotor UAV and the suspended load are shown in Figures5 and 7. The tracking error is illustrated in Figure 8. The path tracking was performed successfully, from take-off to landing, and throughout the different paths that compose the trajectory, including yaw angle regulation. These results demonstrate the joint performance of the designed MPC controller and the UKF estimator, using the adopted control structure.Figure 7
Time evolution of the regulated variables (x,y,z,ψ) and their desired trajectories.Figure 8
Time evolution of the tracking errors. The oscillatory behavior around zero is due to the presence of external disturbances and noisy measurements.During the vertical take-off, the aircraft starts to fly while the load remains in the ground. Only when the distance between the aircraft and the floor is greater than the rope’s length, the tilt-rotor UAV starts to carry the load. Vice versa, in the landing maneuver, the aircraft flies free of load once it has touched the ground. Figure9 highlights this behavior in the first 5 and last 10 seconds of simulation. At the beginning of the simulation, the rope’s length increases until it reaches its maximum value; then, the aircraft starts to carry the suspended load and the total mass increases. Likewise, during the landing, last 10 seconds, the load touches the ground and the rope’s length decreases until the aircraft landing finishes. Also, Figure 9 shows the total mass reduction when the load touches the ground. The controller was capable of dealing with this problem due to the incorporation of lt in the time-varying parameters vector ζt and the model-based nature of predictive controllers. The load’s mass was not incorporated in ζt for two main reasons: (i) in order to consider the mass variation as a time-varying parameter in the same way that lt was considered, the load’s mass needs to be estimated or informed to the controller before the flight starts; (ii) it is reasonable to assume that the relation between the UAV’s mass and the load’s mass is small enough to be rejected by the controller as parametric uncertainty when the actual load’s mass is different from the one considered in Table 1.Figure 9
Analysis of the time evolution of the load and UAV altitude in take-off and landing maneuvers together with the rope’s length and the total mass variation.Figure10 shows the time evolution of the remaining degrees of freedom, which are kept stable as the trajectory is performed. Since the aircraft’s behavior is described implicitly by those variables, one can conclude that the UAV was stabilized along the trajectory. It is important to note that the designed MPC controller was able to stabilize the aircraft without the need of a cascade control structure, since the proposed modeling considered all the dynamic couplings between the rigid bodies.Figure 10
Time evolution of the remaining degrees of freedom. The variations onϕ, θ, γ1, and γ2 are due to changes on their equilibrium points when external disturbances affect the system.Figure11 shows the actuator signals generated by the MPC controller. Despite being noisy, the control signals did not reach the actuators’ limits. In an experimental setup, the inherent inertial properties of the aircraft actuators, which are commonly brushless DC motors and servomotors, would attenuate such noise naturally.Figure 11
Applied thrusts and torques to the tilt-rotor UAV by the propellers and servomotors, respectively.Figures12 and 13 present, respectively, the estimation error for the generalized coordinates and for their time derivatives, along with the confidence limits (i.e., three times the standard deviation). In all cases, the estimation error is close to zero and is kept inside the confidence limits, demonstrating the performance of the unscented Kalman filter. Recalling the nonlinearities of the system (35) and the measurement equation (42), the UKF was then able to estimate the means and covariances of the posterior distributions in a consistent manner, despite the sensors’ different sampling rates.Figure 12
Time evolution of the estimation error of the generalized coordinates, along with the confidence limits. Patterns are highlighted to illustrate the estimation behavior in presence of different sampling rates.Figure 13
Time evolution of the estimation error of the generalized coordinates derivatives, along with the confidence limits. Patterns are highlighted to illustrate the estimation behavior in presence of different sampling rates.Some patterns arose in the confidence limits due to the greater sampling rates of the GPS and the camera, whose measurements are available only every 120 ms and 24 ms for performing the data assimilation step, respectively. The other sensors’ data are available every 12 ms, which is also the controller sampling time. These patterns are expected since the estimation is more accurate and has smaller confidence limits every time that more data are available to be used in the data assimilation step. This result illustrates that, despite having nonlinear dynamics and sensors with different sampling rates, the unscented Kalman filter is able to recover the state vector from the information provided only by the UAV’s embedded sensors.Based on the presented results, the proposed modeling, control and estimation strategies were demonstrated to be appropriate to solve the problem of path tracking control of the suspended load. A next research step consists in validating the proposed approach considering an experimental setup.


## 7. Conclusions

This paper proposed a detailed model and the design of a control and state estimation strategy to solve the path tracking control problem of a suspended load using a tilt-rotor UAV when it is operating in the helicopter flight-mode.A modeling approach was presented, in which the kinematics of the system were formulated from the load’s perspective, being the load’s position and orientation chosen as degrees of freedom of the multibody mechanical system. The UAV’s position and orientation were described only with respect to the load. By using the Euler-Lagrange formulation, the equations of motion were obtained, while taking into account the dynamic coupling between the aircraft and the load, the existence of viscous friction at the suspension point and at the tilting mechanisms, and also external forces affecting the load. The nonlinear state-space representation of the system was obtained, with the load’s position and orientation as state variables.Assuming that the UAV’s embedded sensors provide noisy information with different sampling rates, an unscented Kalman filter was proposed for nonlinear state estimation of all the state variables, based on the model attained and measurement equations developed from kinematic concepts. Moreover, based on linearized time-varying state-space equations augmented with integral actions, a model predictive control strategy was designed for path tracking of the suspended load with stabilization of the tilt-rotor UAV. The proposed model predictive controller allows yaw angle tracking, take-off and landing maneuvers, situations where the rope’s length and the total system mass are not constant. The feedback control loop was performed using the estimation provided by the UKF.The proposed strategy in this work was evaluated through numerical experiments in MATLAB/Simulink environment. The trajectory performed by the load comprised vertical take-off, straight line following with direction changing, and vertical landing, with disturbances being applied to the load and also under total mass and rod’s length variation. The unscented Kalman filter was able to estimate the entire state vector from the information provided by the sensors, while the MPC controller was able to perform the control task using the provided estimation, without saturating the aircraft actuators. The presented results demonstrated the good performance of the designed MPC controller and the unscented Kalman filter for path tracking of the load using a tilt-rotor UAV, in this challenging, not yet addressed scenario.One of the biggest appeals of tilt-rotor UAVs is their ability to take-off and land vertically, perform hover flights, and achieve improved forward velocities when in airplane-like mode. This work dealt only with helicopter-like flights, in which the thrusters’ groups motions are limited to small tilting angles. However, as stated before, the inclusion of aerodynamic surfaces into the proposed model is straightforward. Further, the proposed estimation and control strategy is general enough to cope with both flight modes, including the transition between them.


### 7.1. Future Works and Limitations

The main limitation of the proposed approach is the computational cost associated with the estimation and control algorithms. Since one of the goals of this work was to perform yaw angle tracking, the use of a nonlinear estimator and a controller based on a linearized time-varying model became necessary. As discussed in Section4, this model is not affine in the parameters; that is, it is not possible to obtain a polytopic representation for the linearized model. This fact prevents the use of techniques that could allow the prediction horizon reduction, which is directed related to the controller computational cost. Regarding the estimator, the requirement to perform yaw motion on flight, which makes the linearized model unsuitable far from the chosen operational points, together with the inability to derive Jacobian and Hessian matrices from the model derived in Section 2 using numerical computation software, made the choice of a filter with higher computational cost necessary.Future works will propose ways to reduce the computational cost of the proposed controller allowing implementing it in embedded systems. Approaches that use piecewise linearized models or models identified in subspace models could be a possible way to address this issue. Also, studying computational efficient ways to implement UKF algorithms will be a future work. Once the problems related to computational cost have been solved and the prototype is finished, future works will present the aircraft aerodynamic design, that is currently being developed, together with simulation and experimental results for load transportation tasks.


## 7.1. Future Works and Limitations

The main limitation of the proposed approach is the computational cost associated with the estimation and control algorithms. Since one of the goals of this work was to perform yaw angle tracking, the use of a nonlinear estimator and a controller based on a linearized time-varying model became necessary. As discussed in Section4, this model is not affine in the parameters; that is, it is not possible to obtain a polytopic representation for the linearized model. This fact prevents the use of techniques that could allow the prediction horizon reduction, which is directed related to the controller computational cost. Regarding the estimator, the requirement to perform yaw motion on flight, which makes the linearized model unsuitable far from the chosen operational points, together with the inability to derive Jacobian and Hessian matrices from the model derived in Section 2 using numerical computation software, made the choice of a filter with higher computational cost necessary.Future works will propose ways to reduce the computational cost of the proposed controller allowing implementing it in embedded systems. Approaches that use piecewise linearized models or models identified in subspace models could be a possible way to address this issue. Also, studying computational efficient ways to implement UKF algorithms will be a future work. Once the problems related to computational cost have been solved and the prototype is finished, future works will present the aircraft aerodynamic design, that is currently being developed, together with simulation and experimental results for load transportation tasks.



---

*Source: 9095324-2017-10-19.xml*