# A Load Transportation Nonlinear Control Strategy Using a Tilt-Rotor UAV

Authors: Guilherme V. Raffo; Marcelino M. de Almeida
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

This paper proposes a nonlinear control strategy to solve the trajectory tracking problem of a tilt-rotor Unmanned Aerial Vehicle (UAV) when transporting a suspended load. For the present study, the aim of the control system is to track a desired trajectory of the aircraft with load’s swing-free, even in the presence of external disturbances, parametric uncertainties, unmodeled dynamics, and noisy position measurements with lower sampling frequency than the controller. The whole system modeling is obtained through the Euler-Lagrange formulation considering the dynamics of the tilt-rotor UAV coupled to the suspended load. As for the nonlinear control strategy, an inner-loop control is designed based on input-output feedback linearization combined with the dynamic extension approach to stabilize the attitude and altitude of the UAV assuming nonlinearities, while an outer-loop control law is designed for guiding the aircraft with reduced load swing. The linearized dynamics are controlled using linear mixedH
2
/
H
∞ controllers with pole placement constraints. The solution is compared to two simpler control systems: the first one considers the load as a disturbance to the system but does not avoid its swing; the second one is a previous academic result with a three-level cascade strategy. Finally, in order to deal with the problem of position estimation in presence of unknown disturbances and noisy measurements with low sampling frequency, a Linear Kalman Filter with Unknown Inputs is designed for estimating both the aircraft’s translational position and translational disturbances. Simulation results are carried out to corroborate the proposed control strategy.


---

## Body

## 1. Introduction

Research in Unmanned Aerial Vehicles (UAVs) has gained much attention in the last years, mainly for its variety of applications. Some examples of uses for UAVs might be listed as cargo transportation and delivery, surveillance, field recognition, cave exploration, cinematographic filming, military purposes, 3D mapping, search and rescue, and wildlife research, among many others. The present work focuses on studying the problem of suspended load transportation. This kind of task is extremely important in some missions such as search and rescue, surface exploration, military applications, and personal assistance, among others. Piloted aircraft usually require an experienced pilot to transport suspended payloads to a destination while avoiding any accident with the load and the aircraft itself. Fully autonomous UAVs, on the other hand, are required to use more sophisticated control laws so as to achieve similar (or even better) performances.The most common UAVs that are studied in academia are helicopters, quadrotors, and fixed-wing airplanes. On one hand, rotorcrafts have the advantage over airplanes for performing Vertical Take-Off and Landing (VTOL), while on the other hand, airplanes are able to obtain higher forward flights with greater range and autonomy. Aiming to combine the advantages of both kinds of aircraft, present research on UAVs is increasingly interested in the tilt-rotor UAV, a hybrid copter-plane aircraft. The tilt-rotor is a type of aircraft that combines the vertical lift capacity of helicopters with the range, autonomy, and speeds of fixed-wing airplanes. For missions of search and rescue, for example, tilt-rotors might stand out since it can reach disaster zones faster than rotorcrafts, while also being able to hover over some position of interest.This work proposes a robust nonlinear control strategy for the tilt-rotor UAV in the helicopter flight-mode with the further requirement that it should track a desired trajectory carrying a suspended load. Furthermore, the control system should maintain both the aircraft and the load stable even in the presence of external disturbances, parametric uncertainties, unmodeled dynamics, and noisy position measurements with lower sampling frequency than the controller.The most notorious research works assessing control of tilt-rotor UAVs started being published after 2005. In [1] a back-stepping strategy was applied to a tilt-rotor with two degrees of freedom on each rotor. The use of two degrees of freedom was later abandoned due to difficulties on practical implementation. Reference [2] was able to experimentally maintain a tilt-rotor in hovering using nonlinear control in the vicinity of the equilibrium point. In [3] an adaptation to the previous solution was performed by including coupling gyroscopic body effects on the system’s model used for control design. Reference [4] explored the use of gain-scheduling to control a tilt-rotor’s roll and pitch by choosing a vast number of linearization points, with results being presented in [5]. A nonlinear control of a tilt-rotor UAV was proposed in [6] by means of cascade control using feedback linearization over a numerical model of the aircraft obtained on wind-tunnel tests. The work of [7] presented an approach to control a tilt-rotor UAV using Fuzzy Logic Control. In [8], a Model Predictive Control (MPC) was designed for the attitude of the aircraft. Reference [9] derived a simplified Euler-Lagrange model for the tilt-rotor UAV and used it to design a back-stepping control strategy. In [10], linear H
∞ and mixed H
2
/
H
∞ controllers were designed based on LMI (Linear Matrix Inequality) approach for trajectory tracking of a tilt-rotor UAV in helicopter flight-mode. In order to cover a large range of forward velocity of a tilt-rotor UAV, a robust adaptive mixing control strategy was proposed in [11].Regarding the load transportation control problem, many research works are found in the literature. Nonlinear controllers were introduced in [12, 13] for stabilization of suspended loads in crane operations. Quadrotor UAV was used in [14] to transport a suspended load from one desired point to another applying a machine learning approach to avoid load swing. Reference [15] also used a quadrotor UAV to stabilize the swing of a suspended load with unknown mass by combining a Proportional-Derivative controller with Retrospective Cost Adaptive Control. In [16], trajectories are generated to a quadrotor UAV so that a suspended load passes through a desired trajectory.Some works also address the problem of load transportation using a tilt-rotor UAV. In [17, 18], MPC strategies based on linearized models around desired trajectories were proposed for trajectory tracking of a tilt-rotor UAV with reduced load’s swing. In [19], the trajectory tracking problem of the suspended load was solved through the design of control and state estimation strategies based on linearized, time-invariant state-space equations but did not allow yaw angle tracking, or the occurrence of changes in the load’s mass and rope’s length. This latter work was improved in [20], where an MPC based on a linear time-varying model was designed to perform trajectory 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. Nevertheless, the above control strategies, developed to solve the load transportation problem using a tilt-rotor UAV, are based on linearized models, which limit the domain of attraction of the closed-loop control system. In order to improve that, a nonlinear cascade control strategy was proposed in [21, 22] for trajectory tracking of a tilt-rotor UAV with load’s swing-free. Although this control strategy enlarged the domain of attraction, this nonlinear solution used a three-level cascade strategy, which might not be very attractive from the control point of view, since outer-loops are capable of destabilizing inner-loops if the design is not properly tuned.This paper improves the results presented in [21, 22], in which the position of the aircraft was assumed to be perfectly known, by proposing a two-level cascade strategy, reducing computational costs and attaining a solution whose nonlinear controllers are simpler to tune due to lower number of cascaded loops. Each level of the cascade system executes a control law through the method of input-output feedback linearization (IOFL), which constantly linearizes the system so that linear control design techniques can be applied. In summary, the work in [21, 22] uses two loops to control altitude and attitude, while this work achieves the same in a single loop by using a dynamic extension approach for controlling these variables by actuating on their snaps. The external loop presented in this paper performs trajectory tracking for translational motion while reducing the load’s swing. Besides, the proposed control strategy is designed based on a detailed whole-body dynamic model of the tilt-rotor UAV using the Euler-Lagrange formulation.Some model simplifications are assumed in the control design, neglecting some dynamic cross coupling between generalized coordinates. These assumptions can be close to reality if the system’s angular velocities and generalized acceleration are not very high, which is acceptable when the tilt-rotor UAV operates in the helicopter flight-mode. In order to deal with those neglected terms, mixedH
2
/
H
∞ controllers with pole placement constraints are designed based on the linearized dynamics with the addition of integral terms, featuring robustness against unmodeled dynamics and constant disturbance rejection, while guaranteeing satisfactory time response.In order to solve the problem of position and speed estimation in presence of unknown disturbances and noisy measurements with low sampling frequency, a Linear Kalman Filter with Unknown Inputs (LKFUI) is designed for estimating the aircraft’s translational position and speed and the corresponding disturbances. It assumes that the position is actually measured by a positioning system (e.g., GPS, vision system) equipment with sampling timeT
s
>
τ
s, where τ
s is the controller sampling time. The estimator evaluates the position of the aircraft when no new measurements are available from the sensor, also taking into account its measurement uncertainty. The estimator also considers that the aircraft’s motion may be affected by disturbances (e.g., wind), which are not usually measured. The designed estimator is based on the formulation proposed by [23], which presents a generalized approach to Kalman Filters calledGain-Constrained Kalman Filtering in which both LKFUI and LKF are particular solutions.The proposed control system performance is evaluated through simulation results and compared to a simpler controller, whose design considers the load as a disturbance to the system but does not avoid its swing, and to the three-level cascade strategy proposed in [22].The remainder of the paper is structured as follows: Section2 presents the whole-body modeling of the tilt-rotor UAV with suspended load; Section 3 presents the proposed control structure describing the design of the aforementioned nonlinear control system with the linear mixed H
2
/
H
∞ controllers. In addition, an estimation model for the system is presented in order to develop the LKFUI algorithm; Section 4 shows simulation results; and, finally, Section 5 concludes this paper and provides suggestions of future works.Notation. Let i
→
≜
1
0
0
T, j
→
≜
0
1
0
T, k
→
≜
0
0
1
T be the unit vectors in the direction of the x-, y-, and z-axes, respectively, of a proper reference frame A; R
A
B represents the rotation matrix from frame A to frame B; R
i
,
φ expresses the rotation matrix of an angle φ around axis i; the vector ω
B
A
B
(
t
)
∈
R
3 denotes an angular velocity of frame B with respect to frame A represented in frame B; S
(
p
) is the skew symmetric matrix related to the vector p
∈
R
3. Let r
i denote the relative degree of a system’s output h
i; L
f
h
i is the Lie derivative operation. Let the state vectors x, x
k denote the value of the states at sampling instant k, x
^
k being their estimated values; y
k contains the values measured at instant k, while y
^
k represents the estimation of the same vector; x
^
k
∣
k
-
1 denotes the fact that x at instant k has been predicted using measured information up to instant k
-
1. Let T
s and τ
s be two different sampling times with T
s
>
τ
s and T
s
/
τ
s an integer; x
^
i
,
k (i
=
0,1
,
2
,
…
,
T
s
/
τ
s) denotes the estimated states of the system i
·
τ
s seconds after the instant given at x
^
k
∣
k, where index i represents increments in the controller’s cycles, while index k represents sampling instants of sensors. A complementary list of nomenclature is provided in the Section “Nomenclature.”


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

This section presents the whole-body dynamic modeling of the tilt-rotor UAV with suspended load using the Euler-Lagrange formulation. The multibody system is composed of four rigid bodies (see Figure1): the main body (carbon-fiber structure, landing gear, battery, and electronic devices); two thrusters groups, one on each side of the aircraft (tiltable mechanisms with rotors), connected to the main body by a revolute joint; and the suspended load, which is assumed to be attached to the main body via a rigid rod with negligible mass. The following coordinate frames used in the system modeling are defined as (see Figure 1) a fixed inertial frame I, a moving frame B rigidly attached to the main body, a frame C
1 rigidly attached to the main body’s center of mass, frames C
2 and C
3 rigidly attached to the rotation axes of the right and left tiltable mechanisms, respectively, and a frame C
4 rigidly attached to the suspended load’s center of mass.Figure 1
Tilt-rotor UAV frames and variables definition.


### 2.1. Kinematics

The UAV’s attitude with respect to the frameB is described by Euler angles around the local axes with the roll-pitch-yaw convention η
=
ϕ
θ
ψ
T. In addition, R
B
I represents the rotation matrix from frame B to I. The translation between the origins of frames I and B is represented by ξ
=
x
I
y
I
z
I
T, while d
i
B
=
d
x
i
B
d
y
i
B
d
z
i
B
T is the translation between the origins of frames B and C
i, for i
=
1,2
,
3,4. In addition, note that d
1
B, d
2
B, and d
3
B are all constants, while d
4
B varies due to the suspended load degrees of freedom (DOF). Also, this model assumes that d
x
2
B
=
d
x
3
B
=
0; that is, the rotors are aligned with frame B in the x direction (see Figure 1).The vectord
4
B can be calculated through a parametrization that considers thesuspended load by a cable as asimple pendulum, from which a weight (the load) is connected to the aircraft via a massless rigid rod of length l, and two degrees of freedom represented by γ
1 and γ
2 (rotations around x
B and y
B, respectively). The Forward Kinematic Model of the pendulum subsystem with respect to the aircraft’s body is given by d
4
B
=
R
C
4
B
(
-
l
k
→
), where R
C
4
B
=
R
y
,
γ
2
R
x
,
γ
1 is the load’s attitude with respect to frame B.In addition, frameC
1, which is rigidly attached to the main body’s center of mass, is assumed parallel to frame B, leading to R
C
1
B
=
I
3
×
3. As for the rotors’ attitudes, the rotation matrices are obtained by R
C
2
B
=
R
x
,
-
β
R
y
,
α
R and R
C
3
B
=
R
x
,
β
R
y
,
α
L, α
R, α
L being the right and left tilting angles, respectively. Moreover, β is a small fixed tilt angle of the rotors toward the origin of B and dependent of the aircraft’s mechanical design (see Figure 1). The introduction of β improves the controllability for this kind of VTOL aircraft, although the thrust efficiency is decreased [24].Finally, the position of a point rigidly attached to frameC
i with respect to the inertial frame I is obtained as(1)
p
i
I
=
R
B
I
R
C
i
B
p
i
C
i
+
d
i
B
+
ξ
,
i
=
1,2
,
3,4
.The generalized coordinates vector is given byq
=
ξ
T
η
T
α
T
γ
T
T
∈
R
10, where α
=
α
R
α
L
T and γ
=
γ
1
γ
2
T.


### 2.2. Equations of Motion

This section derives the Euler-Lagrange equations of motion for the tilt-rotor UAV with suspended load. In these equations the exogenous forces are separated into dissipative viscous friction forces and unknown external disturbances. Thus, the equations of motion are written as(2)
M
q
q
¨
+
C
q
,
q
˙
q
˙
+
G
q
=
F
q
+
F
e
x
t
+
F
d
r
a
g
,where M
(
q
)
∈
R
10
×
10 is the inertia matrix, C
(
q
,
q
˙
)
∈
R
10
×
10 is the Coriolis and centripetal forces matrix obtained using the Christoffel symbols of the first kind, G
(
q
)
∈
R
10 is the gravitational force vector, F
(
q
)
∈
R
10 is the independent generalized input force vector, F
e
x
t
∈
R
10 is the external disturbances vector, and F
d
r
a
g
∈
R
10 is the generalizeddrag force vector. The drag forces are assumed to be proportional to the generalized velocities and are given by F
d
r
a
g
=
-
μ
q
˙, with μ
∈
R
10
×
10 being a diagonal constant matrix. Consequently, (2) can be rewritten as(3)
M
q
q
¨
+
C
q
,
q
˙
+
μ
q
˙
+
G
q
=
F
q
+
F
e
x
t
.The inertia matrix can be found by calculating the system’s kinetic energy and expressing it in the formK
=
1
/
2
q
˙
T
M
(
q
)
q
˙. The kinetic energy of the whole system is given by the sum of the individual kinetic energies of each body K
=
∑
i
=
1
4
‍
K
i, where the kinetic energy of the i
t
h body can be obtained from the volume integral [25] (4)
K
i
q
=
1
2
∫
V
i
ρ
i
v
i
I
T
v
i
I
d
V
i
,with ρ
i the mass density at body i. The vector v
i
I is the velocity of a point of body i with respect to frame I and is given by the time derivative of (1):(5)
p
˙
i
I
=
R
˙
B
I
R
C
i
B
p
i
C
i
+
d
i
B
+
R
B
I
R
˙
C
i
B
p
i
C
i
+
R
C
i
B
p
˙
i
C
i
+
d
˙
i
B
+
ξ
˙
.Points p
i
C
i are rigidly attached to their respective frames, leading to p
˙
i
C
i
=
0
3
×
1 for i
=
1,2
,
3,4. Translations d
i
B for i
=
1,2
,
3 are constant, resulting in d
˙
1
B
=
d
˙
2
B
=
d
˙
3
B
=
0
3
×
1. Moreover, since the body’s frame B is fixed with respect to its center of mass’ frame C
1, then R
˙
C
1
B
=
0
3
×
3.With use of properties for skew symmetric matrices (S
(
p
)
q
=
S
(
q
)
T
p, S
(
a
p
+
b
q
)
=
a
S
(
p
)
+
b
S
(
q
), and S
(
R
p
)
=
R
S
(
p
)
R
T; besides, R
˙
B
A
(
t
)
=
R
B
A
(
t
)
S
(
ω
B
A
B
(
t
)
).), equation (5) can be rewritten for each body in the form(6)
p
˙
1
I
=
R
B
I
R
C
1
B
S
p
1
C
1
T
R
C
1
B
T
ω
B
I
B
+
R
B
I
S
d
1
B
T
ω
B
I
B
+
ξ
˙
,
(7)
p
˙
i
I
=
R
B
I
R
C
i
B
S
p
i
C
i
T
R
C
i
B
T
ω
B
I
B
+
R
B
I
S
d
i
B
T
ω
B
I
B
+
R
B
I
R
C
i
B
S
p
i
C
i
T
ω
C
i
B
C
i
+
ξ
˙
;
i
=
2,3
,
(8)
p
˙
4
I
=
R
B
I
R
C
4
B
S
p
4
C
4
T
R
C
4
B
T
ω
B
I
B
+
R
B
I
S
d
4
B
T
ω
B
I
B
+
R
B
I
R
C
4
B
S
p
4
C
4
T
ω
C
4
B
C
4
+
R
B
I
d
˙
4
B
+
ξ
˙
.Assuming that all the system’s bodies are symmetric and each frameC
i coincides with the center of mass of the i
t
h body, then ∫
‍
V
i
ρ
i
p
i
C
i
d
V
i
=
0
3
×
1. Using this property and substituting (6)–(8) into (4), the kinetic energies of the system’s bodies are found to be K
1
=
X
1
′, K
2
=
X
2
′
+
Y
2
′, K
3
=
X
3
′
+
Y
3
′, and K
4
=
X
4
′
+
Y
4
′
+
Z
4
′, where X
i
′, Y
i
′, and Z
i
′ are given by(9)
X
i
′
=
1
2
m
i
ξ
˙
T
ξ
˙
-
m
i
ξ
˙
T
R
B
I
S
d
i
B
ω
B
I
B
+
1
2
ω
B
I
B
T
J
i
ω
B
I
B
,
(10)
Y
i
′
=
ω
B
I
B
T
R
C
i
B
I
i
ω
C
i
B
C
i
+
1
2
ω
C
i
B
C
i
T
I
i
ω
C
i
B
C
i
,
(11)
Z
i
′
=
ω
B
I
B
T
m
i
S
d
i
B
d
˙
i
B
+
1
2
d
˙
i
B
T
m
i
d
˙
i
B
+
ξ
˙
T
m
i
R
B
I
d
˙
i
B
,with m
i the mass of the i
t
h body, I
i
≜
∫
‍
S
(
p
i
C
i
)
T
S
(
p
i
C
i
)
d
m the inertia tensor of the i
t
h body with respect to frame C
i, and J
i
≜
R
C
i
B
I
i
(
R
C
i
B
)
T
+
m
i
S
(
d
i
B
)
T
S
(
d
i
B
) the inertia tensor of body i expressed in frame B (Steiner’s theorem [26]).In order to represent the kinetic energy as a function of the generalized coordinates, the following mappings are applied(12)
ω
B
I
B
=
1
0
-
s
θ
0
c
ϕ
s
ϕ
c
θ
0
-
s
ϕ
c
ϕ
c
θ
ϕ
˙
θ
˙
ψ
˙
=
W
η
η
˙
,
(13)
ω
C
4
B
C
4
=
1
0
0
c
γ
1
0
-
s
γ
1
γ
˙
1
γ
˙
2
=
W
γ
γ
˙
,
(14)
d
˙
4
B
=
l
s
γ
1
s
γ
2
-
l
c
γ
1
c
γ
2
l
c
γ
1
0
l
s
γ
1
c
γ
2
l
c
γ
1
s
γ
2
γ
˙
1
γ
˙
2
=
L
γ
γ
˙
,
(15)
ω
C
2
B
C
2
=
α
˙
R
j
→
,
(16)
ω
C
3
B
C
3
=
α
˙
L
j
→
,with c
x
≜
cos
⁡
x and s
x
≜
sin
⁡
x.Then, the inertia matrixM
(
q
) can be written as(17)
M
q
=
m
I
3
×
3
R
B
I
H
W
η
0
3
×
1
0
3
×
1
m
4
R
B
I
L
γ
∗
W
η
T
J
W
η
m
23
m
24
m
25
∗
∗
j
→
T
I
2
j
→
0
0
1
×
2
∗
∗
∗
j
→
T
I
3
j
→
0
1
×
2
∗
∗
∗
∗
m
55
,where the ∗ terms indicate symmetry with respect to the main diagonal, m
=
∑
‍
m
i, J
=
∑
‍
J
i, H
=
-
S
∑
‍
m
i
d
i
B, m
23
=
W
η
T
R
C
2
B
I
2
j
→, m
24
=
W
η
T
R
C
3
B
I
3
j
→, m
25
=
W
η
T
R
C
4
B
I
4
W
γ
+
m
4
W
η
T
S
(
d
4
B
)
L
γ, and m
55
=
m
4
L
γ
T
L
γ
+
W
γ
T
I
4
W
γ.The gravitational force vectorG
(
q
) is obtained by calculating G
(
q
)
=
∂
P
/
∂
q, where P
=
∑
i
=
1
4
‍
P
i is the sum of the potential energies of the individual bodies given by P
i
=
∫
V
i
‍
ρ
i
(
g
I
)
T
p
i
I
d
V
i, with g
I
=
-
g
z
k
→ the gravity vector with respect to the inertial frame. Assuming again that all the system’s bodies are symmetric and each frame C
i coincides with the center of mass of the i
t
h body, then ∫
‍
V
i
ρ
i
p
i
C
i
d
V
i
=
0
3
×
1 and the potential energy of the whole system is given by(18)
P
=
-
g
I
T
R
B
I
∑
i
=
1
4
m
i
d
i
B
+
m
ξ
.The generalized force vectorF
(
q
)
=
(
T
ξ
I
)
T
(
τ
η
I
)
T
τ
α
T
τ
γ
T
T is composed of all translational forces T
i and rotational torques τ
k performing actuation on the system’s generalized coordinates. In order to calculate T
ξ
I, it is necessary to decompose the force provided by each propeller along frame B
(19)
F
R
B
=
R
C
2
B
k
→
f
R
=
r
R
f
R
,
(20)
F
L
B
=
R
C
3
B
k
→
f
L
=
r
L
f
L
,where f
R and f
L are the right and left propeller thrusts, respectively. By defining F
B
=
F
R
B
+
F
L
B, the translational forces expressed in the inertial frame are then given by(21)
T
ξ
I
=
R
B
I
F
B
=
R
B
I
r
R
R
B
I
r
L
f
R
f
L
.The aircraft’s rotational torques are obtained by adding the torque generated by the thrust of the propellers to the torque caused by the drag of the propellers. The drag torque generated by each propeller is assumed in steady-state and given byτ
d
r
a
g
=
k
τ
/
b
f, where k
τ and b are estimated aerodynamic constants, and f is the vertical thrust of the given propeller. Thus, the main body’s rotational torques expressed in the inertial frame are written as(22)
τ
η
I
=
W
η
T
-
c
α
R
c
β
d
y
-
k
τ
b
s
α
R
c
α
L
c
β
d
y
+
k
τ
b
s
α
L
s
α
R
d
z
+
k
τ
b
c
α
R
s
β
s
α
L
d
z
+
k
τ
b
c
α
L
s
β
s
α
R
d
y
+
k
τ
b
c
α
R
c
β
-
s
α
L
d
y
-
k
τ
b
c
α
L
c
β
f
R
f
L
=
W
η
T
τ
R
τ
L
f
R
f
L
,where d
y
=
d
y
2
B
=
d
y
3
B and d
z
=
d
z
2
B
=
d
z
3
B.Finally, definingτ
α
=
τ
α
R
τ
α
L
T, the input generalized force vector can be expressed in an input-affine form(23)
F
q
=
R
B
I
r
R
R
B
I
r
L
0
3
×
2
W
η
T
τ
R
W
η
T
τ
L
0
3
×
2
0
2
×
1
0
2
×
1
I
2
×
2
0
2
×
1
0
2
×
1
0
2
×
2
f
R
f
L
τ
α
=
B
q
Γ
,where Γ
=
f
R
f
L
τ
α
R
τ
α
L
T is the system’s input vector.Remark 1.
By considering the dynamics of the tiltable mechanisms into the tilt-rotor UAV modeling, the equations of motion become affine in the control input, in which the torques applied to tiltable mechanisms,τ
α
R and τ
α
L, are used to manipulate their angular positions. In opposite to this assumption, some works [1–3, 9] consider these angular positions as control inputs, which lead to high nonlinear coupling in the input mapping, besides neglecting the mentioned dynamics, which have considerable time scale.


### 2.3. State-Space Representation of the System

The state-space representation of the Euler-Lagrange’s dynamic equations (3), assuming the state vector x
(
t
)
=
x
1
⋯
x
20
T
=
q
T
q
˙
T
T, can be written as(24)
x
˙
=
q
˙
q
¨
=
q
˙
M
-
1
B
q
Γ
+
F
e
x
t
-
C
q
,
q
˙
+
μ
q
˙
-
G
q
,
=
F
x
,
Γ
,
F
e
x
t
,
=
f
x
+
g
u
x
Γ
+
g
d
x
F
e
x
t
.


## 2.1. Kinematics

The UAV’s attitude with respect to the frameB is described by Euler angles around the local axes with the roll-pitch-yaw convention η
=
ϕ
θ
ψ
T. In addition, R
B
I represents the rotation matrix from frame B to I. The translation between the origins of frames I and B is represented by ξ
=
x
I
y
I
z
I
T, while d
i
B
=
d
x
i
B
d
y
i
B
d
z
i
B
T is the translation between the origins of frames B and C
i, for i
=
1,2
,
3,4. In addition, note that d
1
B, d
2
B, and d
3
B are all constants, while d
4
B varies due to the suspended load degrees of freedom (DOF). Also, this model assumes that d
x
2
B
=
d
x
3
B
=
0; that is, the rotors are aligned with frame B in the x direction (see Figure 1).The vectord
4
B can be calculated through a parametrization that considers thesuspended load by a cable as asimple pendulum, from which a weight (the load) is connected to the aircraft via a massless rigid rod of length l, and two degrees of freedom represented by γ
1 and γ
2 (rotations around x
B and y
B, respectively). The Forward Kinematic Model of the pendulum subsystem with respect to the aircraft’s body is given by d
4
B
=
R
C
4
B
(
-
l
k
→
), where R
C
4
B
=
R
y
,
γ
2
R
x
,
γ
1 is the load’s attitude with respect to frame B.In addition, frameC
1, which is rigidly attached to the main body’s center of mass, is assumed parallel to frame B, leading to R
C
1
B
=
I
3
×
3. As for the rotors’ attitudes, the rotation matrices are obtained by R
C
2
B
=
R
x
,
-
β
R
y
,
α
R and R
C
3
B
=
R
x
,
β
R
y
,
α
L, α
R, α
L being the right and left tilting angles, respectively. Moreover, β is a small fixed tilt angle of the rotors toward the origin of B and dependent of the aircraft’s mechanical design (see Figure 1). The introduction of β improves the controllability for this kind of VTOL aircraft, although the thrust efficiency is decreased [24].Finally, the position of a point rigidly attached to frameC
i with respect to the inertial frame I is obtained as(1)
p
i
I
=
R
B
I
R
C
i
B
p
i
C
i
+
d
i
B
+
ξ
,
i
=
1,2
,
3,4
.The generalized coordinates vector is given byq
=
ξ
T
η
T
α
T
γ
T
T
∈
R
10, where α
=
α
R
α
L
T and γ
=
γ
1
γ
2
T.


## 2.2. Equations of Motion

This section derives the Euler-Lagrange equations of motion for the tilt-rotor UAV with suspended load. In these equations the exogenous forces are separated into dissipative viscous friction forces and unknown external disturbances. Thus, the equations of motion are written as(2)
M
q
q
¨
+
C
q
,
q
˙
q
˙
+
G
q
=
F
q
+
F
e
x
t
+
F
d
r
a
g
,where M
(
q
)
∈
R
10
×
10 is the inertia matrix, C
(
q
,
q
˙
)
∈
R
10
×
10 is the Coriolis and centripetal forces matrix obtained using the Christoffel symbols of the first kind, G
(
q
)
∈
R
10 is the gravitational force vector, F
(
q
)
∈
R
10 is the independent generalized input force vector, F
e
x
t
∈
R
10 is the external disturbances vector, and F
d
r
a
g
∈
R
10 is the generalizeddrag force vector. The drag forces are assumed to be proportional to the generalized velocities and are given by F
d
r
a
g
=
-
μ
q
˙, with μ
∈
R
10
×
10 being a diagonal constant matrix. Consequently, (2) can be rewritten as(3)
M
q
q
¨
+
C
q
,
q
˙
+
μ
q
˙
+
G
q
=
F
q
+
F
e
x
t
.The inertia matrix can be found by calculating the system’s kinetic energy and expressing it in the formK
=
1
/
2
q
˙
T
M
(
q
)
q
˙. The kinetic energy of the whole system is given by the sum of the individual kinetic energies of each body K
=
∑
i
=
1
4
‍
K
i, where the kinetic energy of the i
t
h body can be obtained from the volume integral [25] (4)
K
i
q
=
1
2
∫
V
i
ρ
i
v
i
I
T
v
i
I
d
V
i
,with ρ
i the mass density at body i. The vector v
i
I is the velocity of a point of body i with respect to frame I and is given by the time derivative of (1):(5)
p
˙
i
I
=
R
˙
B
I
R
C
i
B
p
i
C
i
+
d
i
B
+
R
B
I
R
˙
C
i
B
p
i
C
i
+
R
C
i
B
p
˙
i
C
i
+
d
˙
i
B
+
ξ
˙
.Points p
i
C
i are rigidly attached to their respective frames, leading to p
˙
i
C
i
=
0
3
×
1 for i
=
1,2
,
3,4. Translations d
i
B for i
=
1,2
,
3 are constant, resulting in d
˙
1
B
=
d
˙
2
B
=
d
˙
3
B
=
0
3
×
1. Moreover, since the body’s frame B is fixed with respect to its center of mass’ frame C
1, then R
˙
C
1
B
=
0
3
×
3.With use of properties for skew symmetric matrices (S
(
p
)
q
=
S
(
q
)
T
p, S
(
a
p
+
b
q
)
=
a
S
(
p
)
+
b
S
(
q
), and S
(
R
p
)
=
R
S
(
p
)
R
T; besides, R
˙
B
A
(
t
)
=
R
B
A
(
t
)
S
(
ω
B
A
B
(
t
)
).), equation (5) can be rewritten for each body in the form(6)
p
˙
1
I
=
R
B
I
R
C
1
B
S
p
1
C
1
T
R
C
1
B
T
ω
B
I
B
+
R
B
I
S
d
1
B
T
ω
B
I
B
+
ξ
˙
,
(7)
p
˙
i
I
=
R
B
I
R
C
i
B
S
p
i
C
i
T
R
C
i
B
T
ω
B
I
B
+
R
B
I
S
d
i
B
T
ω
B
I
B
+
R
B
I
R
C
i
B
S
p
i
C
i
T
ω
C
i
B
C
i
+
ξ
˙
;
i
=
2,3
,
(8)
p
˙
4
I
=
R
B
I
R
C
4
B
S
p
4
C
4
T
R
C
4
B
T
ω
B
I
B
+
R
B
I
S
d
4
B
T
ω
B
I
B
+
R
B
I
R
C
4
B
S
p
4
C
4
T
ω
C
4
B
C
4
+
R
B
I
d
˙
4
B
+
ξ
˙
.Assuming that all the system’s bodies are symmetric and each frameC
i coincides with the center of mass of the i
t
h body, then ∫
‍
V
i
ρ
i
p
i
C
i
d
V
i
=
0
3
×
1. Using this property and substituting (6)–(8) into (4), the kinetic energies of the system’s bodies are found to be K
1
=
X
1
′, K
2
=
X
2
′
+
Y
2
′, K
3
=
X
3
′
+
Y
3
′, and K
4
=
X
4
′
+
Y
4
′
+
Z
4
′, where X
i
′, Y
i
′, and Z
i
′ are given by(9)
X
i
′
=
1
2
m
i
ξ
˙
T
ξ
˙
-
m
i
ξ
˙
T
R
B
I
S
d
i
B
ω
B
I
B
+
1
2
ω
B
I
B
T
J
i
ω
B
I
B
,
(10)
Y
i
′
=
ω
B
I
B
T
R
C
i
B
I
i
ω
C
i
B
C
i
+
1
2
ω
C
i
B
C
i
T
I
i
ω
C
i
B
C
i
,
(11)
Z
i
′
=
ω
B
I
B
T
m
i
S
d
i
B
d
˙
i
B
+
1
2
d
˙
i
B
T
m
i
d
˙
i
B
+
ξ
˙
T
m
i
R
B
I
d
˙
i
B
,with m
i the mass of the i
t
h body, I
i
≜
∫
‍
S
(
p
i
C
i
)
T
S
(
p
i
C
i
)
d
m the inertia tensor of the i
t
h body with respect to frame C
i, and J
i
≜
R
C
i
B
I
i
(
R
C
i
B
)
T
+
m
i
S
(
d
i
B
)
T
S
(
d
i
B
) the inertia tensor of body i expressed in frame B (Steiner’s theorem [26]).In order to represent the kinetic energy as a function of the generalized coordinates, the following mappings are applied(12)
ω
B
I
B
=
1
0
-
s
θ
0
c
ϕ
s
ϕ
c
θ
0
-
s
ϕ
c
ϕ
c
θ
ϕ
˙
θ
˙
ψ
˙
=
W
η
η
˙
,
(13)
ω
C
4
B
C
4
=
1
0
0
c
γ
1
0
-
s
γ
1
γ
˙
1
γ
˙
2
=
W
γ
γ
˙
,
(14)
d
˙
4
B
=
l
s
γ
1
s
γ
2
-
l
c
γ
1
c
γ
2
l
c
γ
1
0
l
s
γ
1
c
γ
2
l
c
γ
1
s
γ
2
γ
˙
1
γ
˙
2
=
L
γ
γ
˙
,
(15)
ω
C
2
B
C
2
=
α
˙
R
j
→
,
(16)
ω
C
3
B
C
3
=
α
˙
L
j
→
,with c
x
≜
cos
⁡
x and s
x
≜
sin
⁡
x.Then, the inertia matrixM
(
q
) can be written as(17)
M
q
=
m
I
3
×
3
R
B
I
H
W
η
0
3
×
1
0
3
×
1
m
4
R
B
I
L
γ
∗
W
η
T
J
W
η
m
23
m
24
m
25
∗
∗
j
→
T
I
2
j
→
0
0
1
×
2
∗
∗
∗
j
→
T
I
3
j
→
0
1
×
2
∗
∗
∗
∗
m
55
,where the ∗ terms indicate symmetry with respect to the main diagonal, m
=
∑
‍
m
i, J
=
∑
‍
J
i, H
=
-
S
∑
‍
m
i
d
i
B, m
23
=
W
η
T
R
C
2
B
I
2
j
→, m
24
=
W
η
T
R
C
3
B
I
3
j
→, m
25
=
W
η
T
R
C
4
B
I
4
W
γ
+
m
4
W
η
T
S
(
d
4
B
)
L
γ, and m
55
=
m
4
L
γ
T
L
γ
+
W
γ
T
I
4
W
γ.The gravitational force vectorG
(
q
) is obtained by calculating G
(
q
)
=
∂
P
/
∂
q, where P
=
∑
i
=
1
4
‍
P
i is the sum of the potential energies of the individual bodies given by P
i
=
∫
V
i
‍
ρ
i
(
g
I
)
T
p
i
I
d
V
i, with g
I
=
-
g
z
k
→ the gravity vector with respect to the inertial frame. Assuming again that all the system’s bodies are symmetric and each frame C
i coincides with the center of mass of the i
t
h body, then ∫
‍
V
i
ρ
i
p
i
C
i
d
V
i
=
0
3
×
1 and the potential energy of the whole system is given by(18)
P
=
-
g
I
T
R
B
I
∑
i
=
1
4
m
i
d
i
B
+
m
ξ
.The generalized force vectorF
(
q
)
=
(
T
ξ
I
)
T
(
τ
η
I
)
T
τ
α
T
τ
γ
T
T is composed of all translational forces T
i and rotational torques τ
k performing actuation on the system’s generalized coordinates. In order to calculate T
ξ
I, it is necessary to decompose the force provided by each propeller along frame B
(19)
F
R
B
=
R
C
2
B
k
→
f
R
=
r
R
f
R
,
(20)
F
L
B
=
R
C
3
B
k
→
f
L
=
r
L
f
L
,where f
R and f
L are the right and left propeller thrusts, respectively. By defining F
B
=
F
R
B
+
F
L
B, the translational forces expressed in the inertial frame are then given by(21)
T
ξ
I
=
R
B
I
F
B
=
R
B
I
r
R
R
B
I
r
L
f
R
f
L
.The aircraft’s rotational torques are obtained by adding the torque generated by the thrust of the propellers to the torque caused by the drag of the propellers. The drag torque generated by each propeller is assumed in steady-state and given byτ
d
r
a
g
=
k
τ
/
b
f, where k
τ and b are estimated aerodynamic constants, and f is the vertical thrust of the given propeller. Thus, the main body’s rotational torques expressed in the inertial frame are written as(22)
τ
η
I
=
W
η
T
-
c
α
R
c
β
d
y
-
k
τ
b
s
α
R
c
α
L
c
β
d
y
+
k
τ
b
s
α
L
s
α
R
d
z
+
k
τ
b
c
α
R
s
β
s
α
L
d
z
+
k
τ
b
c
α
L
s
β
s
α
R
d
y
+
k
τ
b
c
α
R
c
β
-
s
α
L
d
y
-
k
τ
b
c
α
L
c
β
f
R
f
L
=
W
η
T
τ
R
τ
L
f
R
f
L
,where d
y
=
d
y
2
B
=
d
y
3
B and d
z
=
d
z
2
B
=
d
z
3
B.Finally, definingτ
α
=
τ
α
R
τ
α
L
T, the input generalized force vector can be expressed in an input-affine form(23)
F
q
=
R
B
I
r
R
R
B
I
r
L
0
3
×
2
W
η
T
τ
R
W
η
T
τ
L
0
3
×
2
0
2
×
1
0
2
×
1
I
2
×
2
0
2
×
1
0
2
×
1
0
2
×
2
f
R
f
L
τ
α
=
B
q
Γ
,where Γ
=
f
R
f
L
τ
α
R
τ
α
L
T is the system’s input vector.Remark 1.
By considering the dynamics of the tiltable mechanisms into the tilt-rotor UAV modeling, the equations of motion become affine in the control input, in which the torques applied to tiltable mechanisms,τ
α
R and τ
α
L, are used to manipulate their angular positions. In opposite to this assumption, some works [1–3, 9] consider these angular positions as control inputs, which lead to high nonlinear coupling in the input mapping, besides neglecting the mentioned dynamics, which have considerable time scale.


## 2.3. State-Space Representation of the System

The state-space representation of the Euler-Lagrange’s dynamic equations (3), assuming the state vector x
(
t
)
=
x
1
⋯
x
20
T
=
q
T
q
˙
T
T, can be written as(24)
x
˙
=
q
˙
q
¨
=
q
˙
M
-
1
B
q
Γ
+
F
e
x
t
-
C
q
,
q
˙
+
μ
q
˙
-
G
q
,
=
F
x
,
Γ
,
F
e
x
t
,
=
f
x
+
g
u
x
Γ
+
g
d
x
F
e
x
t
.


## 3. Control Strategy Design

This section presents the proposed robust nonlinear control strategy to solve the load transportation problem using a tilt-rotor UAV. The goal of the control design is to perform trajectory tracking of coordinatesx, y, z, and ψ, while ensuring stability of the remaining generalized coordinates, which includes the stabilization of the load. The strategy is based on a cascade scheme as shown in Figure 2. The inner-loop control law actuates in Γ and is designed by using the input-output feedback linearization (IOFL) method combined with the dynamic extension technique in order to control altitude, z, and attitude, η
=
ϕ
θ
ψ
T, also ensuring stability of the tilting angles α
R and α
L. The outer-loop actuates on the references of ϕ and θ in order to ensure trajectory tracking of the planar position x
-
y, whereas an additional control term is derived to stabilize the load angles (γ
1 and γ
2) with reduced swing.Figure 2
Nonlinear feedback linearization cascade control strategy.When designing linear controllers for the linearized blocks, the performances of coordinatesx, y, z, ϕ, θ, and ψ are achieved through PID-like controllers, ensuring trajectory tracking in presence of constant disturbances, unmodeled dynamics, and parametric uncertainties. As for γ
1 and γ
2, their derivatives with respect to the aircraft are controlled by using PI-like controllers.Furthermore, for control design purposes, some hypotheses are considered on the equations of motion (3) of the tilt-rotor UAV with suspended load, from which some DOF are assumed uncoupled and considered as unmodeled dynamics. These dynamics decoupling assumptions can be close to reality if the system’s angular velocities (the three first columns of the Coriolis Matrix have only zero terms, implying that translational velocity does not affect the dynamics of the remaining generalized coordinates) and generalized acceleration are small. Nevertheless, these kind of uncertainties must be compensated by the control system. Therefore, in order to compensate them and, additionally, provide good transient response, the linear controllers’ gains are obtained by the mixed H
2
/
H
∞ synthesis with pole placement constraints. The following assumptions are considered:(A1) Letd
1
B
=
0
3
×
1 (center of mass of the main body is in its geometric center), m
2
=
m
3 (mass of both rotors are the same), and d
2
B
=
0
d
y
2
B
0
T
=
-
d
3
B; then H
=
-
S
(
m
4
d
4
B
). In addition, assuming that m
4
≪
m
=
∑
‍
m
i, then the inertia matrix’s term m
12
=
R
B
I
H
W
η
≈
0
3
×
3.
(A2) The coupling between the altitude dynamics (z-dynamics) and the load motion is assumed to be negligible; that is, m
15
≈
i
→
j
→
0
3
×
1
m
15.
(A3) The coupling between the attitude dynamics and the tilting angles dynamics is assumed to be negligible; that is,m
23
≈
0
3
×
1 and m
24
≈
0
3
×
1.
(A4) The coupling between the attitude dynamics and the load motion is also assumed to be negligible; that is,m
25
≈
0
3
×
2.From these assumptions, (3) can be split into two uncoupled equations of motion: (i) one representing the dynamics for q
1
=
z
η
T
α
T
T, and (ii) another representing q
2
=
x
y
γ
T
T:(25)
M
q
1
q
¨
1
+
C
q
1
q
˙
1
+
G
q
1
=
B
q
1
Γ
-
μ
q
1
q
˙
1
+
δ
q
1
,
(26)
M
q
2
q
¨
2
+
C
q
2
q
˙
2
+
G
q
2
=
B
q
2
Γ
-
μ
q
2
q
˙
2
+
δ
q
2
,where M
q
i
=
I
q
i
M
I
q
i
T, C
q
i
=
I
q
i
C
I
q
i
T, G
q
i
=
I
q
i
G
q
i, B
q
i
=
I
q
i
B, μ
q
i
=
I
q
i
μ
I
q
i
T, and δ
q
i is a vector including unmodeled dynamics and the unknown external disturbances F
e
x
t, i
∈
{
1,2
}, with(27)
I
q
1
=
0
6
×
2
I
6
×
6
0
6
×
2
,
(28)
I
q
2
=
I
2
×
2
0
2
×
6
0
2
×
2
0
2
×
2
0
2
×
6
I
2
×
2
.The state-space representations of (25)-(26) will be denoted as x
˙
q
i
=
f
q
i
(
x
q
i
)
+
g
u
q
i
(
x
q
i
)
Γ
+
g
d
q
i
(
x
q
i
)
δ
q
i.Additionally, in order to deal with a more realistic scenario, it is assumed that the translational variables are not available during all controller’s sampling time. Therefore, the problem of state estimation for thetranslational position and speed of the aircraft arises, in which it is assumed that the remaining state-space variables are known. It considers that the position is actually measured by a positioning system (e.g., GPS, vision system) equipment with sampling time T
s, while the controller has a sampling time τ
s, with τ
s
<
T
s and assuming T
s
/
τ
s
∈
N
+. Moreover, velocity is assumed to be measured throughout a sensor based in visual odometry with same sampling frequency as the positioning system. Therefore, the estimator must evaluate the position and speed of the aircraft when no new measurements are available from the sensors, also taking into account their measurement uncertainty. In order to be more realistic, the estimator must also consider that the aircraft’s motion may be affected by disturbances (e.g., wind gusts), which are not usually measured.State estimation algorithms typically rely on knowing all the inputs that affect a system; otherwise the estimates become biased due to the unknown input parameters. In order to estimate states with unknown inputs, the techniqueLinear Kalman Filter with Unknown Inputs (LKFUI) is used. Unlike classicLinear Kalman Filter (LKF) [27], the LKFUI [28] can deal with problems involving unknown inputs such as disturbances.


### 3.1. Inner-Loop Control Design

In this section the proposed IOFL control law is designed for controlling system (25) considering the systems’ outputs h
q
1
x
q
1
=
z
ϕ
θ
ψ
T and the system’s inputs Γ
=
f
R
f
L
τ
α
R
τ
α
L
T.Recall that the relative degreer
i of a system’s output h
i is exactly the number of times one has to differentiate h
i in order to have at least one component of the vector input Γ explicitly appearing. It is computed by r
i
=
(
inf
⁡
k
,
∃
j
,
1
≤
j
≤
4
,
L
g
j
L
f
k
-
1
h
i
), where L
f
h
i is the Lie derivative operation and g
j is the j
t
h column of g
u
q
1
(
x
q
i
). It is possible to verify that r
i
=
2
,
∀
i
∈
{
1,2
,
3,4
} (this statement can be trivially understood due to the fact that z
¨, ϕ
¨, θ
¨, and ψ
¨ are all function of at least one input Γ
i for arbitrary α
R and α
L), which implies that this set of inputs/outputs are not fully feedback linearizable, since r
=
∑
‍
r
i
=
8 and the system has twelve state variables (q
1 and q
˙
1).Since it is only possible to regulate on a desired value the number of outputs equal to the number of control inputs, and the time-derivativesz
¨, ϕ
¨, θ
¨, and ψ
¨ are all affected only by the control inputs f
R and f
L (see (23)), and none by the torques τ
R and τ
L, it is required to delay the appearance of f
R and f
L to higher order time-derivatives of z, ϕ, θ, and ψ in order to allow the other control inputs to appear. Therefore, the present problem can be solved by using the dynamic extension technique, in which the state vector of the system is extended by setting f
R and f
L equal to the output of two double integrators and, then, the state variables x
13
=
f
R, x
14
=
f
˙
R, x
15
=
f
L, and x
16
=
f
˙
L are included. The new system’s inputs are defined as Γ
¯
=
f
¨
R
f
¨
L
τ
α
R
τ
α
L
T. Thus, the new state vector is x
¯
q
1
=
q
1
T
q
˙
1
T
x
13
x
14
x
15
x
16
T, and the state-space equations can be rewritten as(29)
x
¯
˙
q
1
=
f
¯
q
1
x
¯
q
1
+
∑
i
=
1
4
g
¯
i
x
¯
q
1
Γ
¯
i
+
g
¯
d
q
1
x
¯
q
1
δ
q
1
,
y
F
L
=
h
q
1
x
q
1
=
z
ϕ
θ
ψ
T
,where f
¯
q
1
x
¯
q
1
=
q
˙
1
T
q
¨
1
T
x
14
0
x
16
0
T, q
¨
1
=
M
q
1
-
1
(
-
[
C
q
1
+
μ
q
1
]
-
G
q
1
)
+
B
1
x
13
+
B
2
x
15, Γ
¯
i is the i
t
h element of the system’s input vector Γ
¯, g
¯
i
(
x
¯
q
1
) is given by(30)
g
¯
1
=
0
13
×
1
1
0
2
×
1
,
(31)
g
¯
2
=
0
15
×
1
1
,
(32)
g
¯
3
=
0
6
×
1
B
3
0
4
×
1
,
(33)
g
¯
4
=
0
6
×
1
B
4
0
4
×
1
,and B
i is the i
t
h column of B
q
1, i
∈
{
1,2
,
3,4
}. Note that the control inputs f
R and f
L of system (24) are not anymore control inputs for the extended system (29), being now state variables. Additionally, it improves the capability of the inner-loop controller to deal with external disturbances.Computing the relative degreer
i for the new system, it is now possible to verify that r
i
=
4
,
∀
i
∈
{
1,2
,
3,4
}. This only holds provided the condition that the terms of the inertia matrix m
23
=
m
24
=
0. If this was not true, then the coupling between η and α would imply that it is physically possible to actuate on attitude angles by Newton’s action-reaction when applying torques to rotate the propellers. Even though this is mathematically attainable, this is not the form of actuation desired for the conception of a tilt-rotor aircraft. Therefore, the system’s relative degree is r
=
∑
‍
r
i
=
16, which means that the system is now fully feedback linearizable. Thus, in order to apply IOFL techniques, the following needs to be calculated:(34)
Δ
x
¯
q
1
∈
R
4
×
4
:
Δ
i
j
x
¯
q
1
=
L
g
¯
j
L
f
¯
3
h
i
x
¯
q
1
,
(35)
b
x
¯
q
1
∈
R
4
×
1
:
b
i
x
¯
q
1
=
L
f
¯
4
h
i
x
¯
q
1
,where Δ
i
j
(
x
¯
q
1
) is the element of the i
t
h row and j
t
h column of Δ
(
x
¯
q
1
) and b
i
(
x
¯
q
1
) is the i
t
h row of b
(
x
¯
q
1
).Then, the following control law may be applied(36)
Γ
¯
=
Δ
-
1
v
q
1
-
b
,where v
q
1
=
v
z
v
ϕ
v
θ
v
ψ
T is the additional control input. The matrix Δ
(
x
¯
q
1
) has full-rank on the neighborhood of the system’s equilibrium point, implying that(37)
h
…
.
q
1
x
¯
q
1
=
v
q
1
+
π
1
δ
q
1
,with π
1
(
δ
q
1
) a function of the residual terms due to model simplifications, unmodeled dynamics, and unknown external disturbancesTherefore, it is then possible to apply linear techniques to obtain the additional control law in order to regulate the outputsz, ϕ, θ, and ψ. These variables are controlled by Proportional-Integral-Derivative- (PID-) like controllers with feed-forward term. This choice, together with a desired faster closed-loop time response, allows assuming the steady-state behavior of the inner-loop as unitary gain at the moment of tuning the outer-loop controllers. Moreover, the variables z, ϕ, θ, and ψ are required to track desired trajectories when constant external disturbances and modeling errors affect the system. Thus, the linear control laws compounding the additional control input vector, v
q
1, are designed as(38)
v
z
=
z
…
.
R
e
f
+
K
d
d
d
z
e
⃛
z
+
K
d
d
z
e
¨
z
+
K
d
z
e
˙
z
+
K
p
z
e
z
+
K
I
z
∫
e
z
d
t
,
(39)
v
ϕ
=
ϕ
…
.
R
e
f
+
K
d
d
d
ϕ
e
⃛
ϕ
+
K
d
d
ϕ
e
¨
ϕ
+
K
d
ϕ
e
˙
ϕ
+
K
p
ϕ
e
ϕ
+
K
I
ϕ
∫
e
ϕ
d
t
,
(40)
v
θ
=
θ
…
.
R
e
f
+
K
d
d
d
θ
e
⃛
θ
+
K
d
d
θ
e
¨
θ
+
K
d
θ
e
˙
θ
+
K
p
θ
e
θ
+
K
I
θ
∫
e
θ
d
t
,
(41)
v
ψ
=
ψ
…
.
R
e
f
+
K
d
d
d
ψ
e
⃛
ψ
+
K
d
d
ψ
e
¨
ψ
+
K
d
ψ
e
˙
ψ
+
K
p
ψ
e
ψ
+
K
I
ψ
∫
e
ψ
d
t
,where e
i
=
(
i
-
i
R
e
f
) with i
∈
{
z
,
ϕ
,
θ
,
ψ
}, and z
…
.
R
e
f, ϕ
…
.
R
e
f, θ
…
.
R
e
f, and ψ
…
.
R
e
f are feed-forward terms due to the time-varying reference.The above control laws will be synthesized using the mixedH
2
/
H
∞ control technique with pole placement constraints in order to ensure reference tracking for the variables z, ϕ, θ, and ψ, while α
R and α
L are stabilized.


### 3.2. Outer-Loop Control Design

The outer-loop IOFL controller is designed in order to control the subsystem (26) by actuating on the references of roll, ϕ
R
e
f, and pitch, θ
R
e
f, angles. Moreover, an additional control term is considered to reduce the load swing.From (23), the generalized control input vector of subsystem (26) can be rewritten as B
q
2
Γ
=
T
x
I
T
y
I
0
0
T. The variables T
x
I and T
y
I are, respectively, the translational forces along x and y expressed in the inertial frame and are obtained from(42)
T
x
I
T
y
I
=
i
→
T
j
→
T
R
B
I
F
R
B
+
F
L
B
=
r
11
r
12
r
13
r
21
r
22
r
23
F
R
B
+
F
L
B
,where r
i
j is the element of R
B
I located at the i
t
h line and j
t
h column. From F
R
B
+
F
L
B
=
f
x
B
f
y
B
f
z
B
T, the actuation approach in this part of the design is performed by changing the projection of f
z
B along x
I and y
I. In other words, the roll (ϕ) and pitch (θ) angles need to be changed so as to obtain the desired projections of f
z
B along x
I and y
I. Besides, the projections of f
x
B and f
y
B along the inertial frame are assumed as known disturbances which should be compensated by the controller.First, in order to perform the feedback linearization approach and design an additional linear control input, from (26) and (42), the vector of dynamics to be canceled n
q
2
=
n
x
n
y
n
γ
1
n
γ
2
T is defined as(43)
n
q
2
=
C
q
2
+
μ
q
2
q
˙
2
+
G
q
2
-
F
q
2
′
,where F
q
2
′
=
r
11
f
x
B
+
r
12
f
y
B
r
21
f
x
B
+
r
22
f
y
B
0
0
T. Consequently, subsystem (26) can be rewritten as(44)
M
q
2
q
¨
2
+
n
x
n
y
n
γ
1
n
γ
2
︸
n
q
2
=
s
ψ
s
ϕ
+
c
ψ
s
θ
c
ϕ
f
z
B
s
ψ
s
θ
c
ϕ
-
c
ψ
s
ϕ
f
z
B
0
0
+
δ
q
2
.By inspecting (44), it is not trivial to obtain a nonlinear state-space representation affine in the control inputs ϕ
R
e
f and θ
R
e
f. To solve this problem, which will allow applying the IOFL technique, and assuming that vector of dynamics n
q
2 is known, one can define the intermediary input values Υ
1 and Υ
2 and state the following nonlinear set of equations:(45)
s
ψ
s
ϕ
+
c
ψ
s
θ
c
ϕ
f
z
B
-
n
x
=
Υ
1
,
s
ψ
s
θ
c
ϕ
-
c
ψ
s
ϕ
f
z
B
-
n
y
=
Υ
2
.As the term f
z
B is obtained from the inner-loop controller, which has a faster settling time than the outer one, for the outer-loop control design, this variable is assumed as a time-varying parameter. Furthermore, since the roll (ϕ) and pitch (θ) angles are controlled variables of the inner-loop, they can be written with respect to their error and reference as e
ϕ
=
(
ϕ
-
ϕ
R
e
f
) and e
θ
=
(
θ
-
θ
R
e
f
). Therefore, being the inner-loop controller in charge of carrying these errors to zero, and assuming that the yaw (ψ) angle is measurable, (45) can be written with respect to ϕ
R
e
f and θ
R
e
f, whose solutions are given by(46)
ϕ
R
e
f
=
arcsin
⁡
σ
s
ψ
Υ
1
+
n
x
-
c
ψ
Υ
2
+
n
y
f
z
B
,
θ
R
e
f
=
arcsin
⁡
σ
c
ψ
Υ
1
+
n
x
+
s
ψ
Υ
2
+
n
y
f
z
B
c
ϕ
R
e
f
,where σ
(
p
) is a saturation function that avoids the a
r
c
s
i
n operation from becoming undefined and is given by σ
(
p
)
=
min
⁡
1
,
max
⁡
-
1
,
p (the function min
⁡
a
,
b returns the minimum value between a and b and max
⁡
c
,
d returns the maximum value between c and d).By replacing the transformation (45) in (44), (26) is reduced to(47)
M
q
2
q
¨
2
+
n
¯
q
2
=
Υ
q
2
+
δ
q
2
,where n
¯
q
2
=
0
0
n
γ
1
n
γ
2
T and Υ
q
2
=
Υ
1
Υ
2
0
0
T. Since (47) is affine in the intermediary input vector Υ
q
2, IOFL techniques can be used to linearize this subsystem.Then, by choosing the outputsh
q
2
x
q
2
=
x
y
T, matrices Δ
(
x
q
2
) and b
(
x
q
2
) can be computed so as to apply a control law similar to the one in (36) leading to(48)
h
¨
q
2
x
q
2
=
v
q
2
+
π
2
δ
q
2
,where v
q
2
=
v
x
v
y
T. Consequently, it is possible to control the aircraft’s motion along x
I and y
I, while γ
1 and γ
2 compose the internal dynamics.Therefore, linear controllers are designed taking into account the motion of the load expressed in the inertial frame. First, the vectord
B
C
4
I
=
R
B
I
d
4
B
=
d
x
I
d
y
I
d
z
I
T is defined as the translation between frame B and the load expressed in frame I. Then, the reference for the load variation is given as d
˙
x
I
d
˙
x
I
T
=
0
0
T. Note that d
4
B is the vector that goes from the origin of frame B to the origin of frame C
4. Hence, d
˙
4
B is the relative velocity between the origins of frames B and C
4. That is, the objective is to avoid load swinging speed regardless of its position with respect to the aircraft. In order to accomplish this task, PID-like controllers for the x
y motion with additional PI-like control terms to regulate the load’s speed with respect to the tilt-rotor are applied:(49)
v
x
=
x
¨
R
e
f
+
K
d
x
e
˙
x
+
K
p
x
e
x
+
K
I
x
∫
e
x
d
t
︸
x
-
m
o
t
i
o
n
+
K
p
d
x
d
˙
x
I
+
K
I
d
x
∫
d
˙
x
I
d
t
︸
l
o
a
d
s
w
i
n
g
,
(50)
v
y
=
y
¨
R
e
f
+
K
d
y
e
˙
y
+
K
p
y
e
y
+
K
I
y
∫
e
y
d
t
︸
y
-
m
o
t
i
o
n
+
K
p
d
y
d
˙
y
I
+
K
I
d
y
∫
d
˙
y
I
d
t
︸
l
o
a
d
s
w
i
n
g
,where x
¨
R
e
f and y
¨
R
e
f are feed-forward terms due to the time-varying reference.When tuning parameters of controllers (49) and (50), it is again important to bear in mind that this strategy assumed null coupling between x, y, γ
1, and γ
2 and the remaining variables, which is incorporated into the residual terms, π
2
(
δ
q
2
). Again, in order to synthesize the controllers’ gains providing good transient response and disturbance attenuation, the mixed H
2
/
H
∞ synthesis with pole placement constraints is employed.


### 3.3. Robust MixedH
2
/
H
∞ Control Design

After applying the input-output feedback linearization procedure to both inner and outer loops, the additional control inputs (38)–(41) and (49)-(50) are synthesized using the mixed H
2
/
H
∞ technique with pole placement constraints [11]. These linearized dynamics can be written as(51)
e
˙
i
=
A
i
e
i
+
B
u
i
u
~
i
+
B
π
i
π
i
,
z
i
=
C
z
i
e
i
+
D
u
z
i
u
~
i
+
D
π
z
i
π
i
,
u
~
i
=
K
i
e
i
,where e
i
=
[
∫
‍
e
i
d
t
e
i
e
˙
i
e
¨
i
e
⃛
i
]
T and K
i
=
[
K
I
i
K
p
i
K
d
i
K
d
d
i
K
d
d
d
i
] with i
∈
{
z
,
ϕ
,
θ
,
ψ
} for the inner-loop, e
i
=
[
∫
‍
e
i
d
t
e
i
e
˙
i
]
T and K
i
=
[
K
I
i
K
p
i
K
d
i
] with i
∈
{
x
,
y
} for the x
y-motion, and e
i
=
[
∫
‍
d
˙
i
I
d
t
d
˙
i
I
]
T and K
i
=
[
K
I
d
i
K
p
d
i
] with i
∈
{
γ
1
,
γ
2
} for the term in charge of load’s swing. Matrices A
i, B
u
i, B
π
i are defined with proper dimensions, and C
z
i, D
u
z
i, D
π
z
i are weighting matrices.The optimalH
∞ controller asymptotically stabilizes the linear system and minimizes the H
∞ norm, H
π
i
z
i
(
s
)
∞
=
s
u
p
(
σ
H
π
i
z
i
(
s
)
), where H
π
i
z
i
(
s
) is the transfer function between the external disturbance π
i and the cost variable z
i. Applying Parseval’s theorem to the H
∞ norm, defined in the frequency domain, it results in z
i
2
≤
H
π
i
z
i
s
∞
π
i
2, in which minimizing the H
∞ norm minimizes the system disturbance effect. Solving the H
∞ optimization problem, the linear system is asymptotically stable with H
∞ norm bounded by H
π
i
z
i
(
s
)
∞
<
γ
∗, where γ
∗ is the optimal value.Considering again system (51) with D
π
z
i
=
0, H
2 control problem leads to an optimal controller that minimizes the H
2 norm h
i
(
t
)
2
=
∫
0
∞
‍
h
i
2
(
t
)
d
t for h
i
(
t
)
=
L
-
1
{
H
π
i
z
i
(
s
)
}
. The minimization of this norm implicates in the optimal transient response of the system. Therefore, the mixed H
2
/
H
∞ control approach is formulated by finding an optimal H
∞ norm and using γ
^
>
γ
∗ to solve the H
2 problem maintaining the H
∞ constraints, providing the advantage of a small H
∞ norm and a better transient response.In addition, to ensure pole placement constraints, the closed-loop poles must be within a defined convex region in the complex plane (D-stability). An LMI region is defined as a subset D of the complex plane C, where D is given by D
=
s
∈
C
:
L
+
s
M
+
s
∗
M
T
<
0, where L
=
L
T and M are matrices that define region D, s
=
σ
+
j
ω, and s
∗ is the complex conjugate of s [29]. In this work the following regions are used: (R
1) R
e
(
s
)
<
-
α, which is defined by s
+
s
∗
<
-
2
α, given that s
+
s
∗
=
σ
+
j
ω
+
σ
-
j
ω
=
2
σ. Therefore, it is possible to rewrite this expression in the form 2
α
+
s
+
s
∗
<
0 and then L
=
2
α and M
=
1; (R
2) is a disk centered at (
c
,
0
) with radius r, in which L
=
-
r
-
c
-
c
-
r and M
=
0
1
0
0; and (R
3) is a cone defined by I
m
s
<
t
a
n
φ
R
e
s, from which the matrices L and M are given as L
=
0
0
0
0 and M
=
s
i
n
(
φ
)
c
o
s
(
φ
)
-
c
o
s
(
φ
)
s
i
n
(
φ
). According to [29], a closed-loop state matrix A
+
B
K is D-stable if and only if there exists a real definite positive symmetric matrix Q
∈
R
n
×
n such that L
⊗
Q
+
M
⊗
V
+
M
T
⊗
V
T
<
0, where ⊗ is the Kronecker product operator, and V
=
(
A
Q
+
B
Y
) with K
=
Y
Q
-
1.Thus, the robust mixedH
2
/
H
∞ control problem with pole placement constraints is formulated for each linearized system (51) through the following optimization problem subject to linear matrix inequalities assuming γ
=
γ
^:(52)
m
i
n
t
r
N
i
,
s
.
t
.
Q
i
>
0
(53)
H
2
c
o
n
s
t
r
a
i
n
t
s
N
i
C
z
i
+
D
u
z
i
K
i
Q
i
Q
i
C
z
i
+
D
u
z
i
K
i
T
Q
i
>
0
,
Q
i
A
i
T
+
A
i
Q
i
+
Y
i
T
B
u
i
T
+
B
u
i
Y
i
+
B
π
i
B
π
i
T
<
0
,
(54)
H
∞
c
o
n
s
t
r
a
i
n
t
s
Υ
B
π
i
Q
i
C
z
i
T
+
Y
i
T
D
u
z
i
T
∗
-
γ
I
n
w
i
D
π
z
i
T
∗
∗
-
γ
I
n
z
i
<
0
,
(55)
R
e
g
i
o
n
R
1
⟶
2
α
Q
+
A
Q
+
Q
A
T
+
B
Y
+
Y
T
B
T
<
0
,
(56)
R
e
g
i
o
n
R
2
⟶
-
r
Q
-
c
Q
+
A
Q
+
B
Y
-
c
Q
+
Q
A
T
+
Y
T
B
T
-
r
Q
<
0
,
(57)
R
e
g
i
o
n
R
3
⟶
s
φ
A
Q
+
Q
A
T
+
B
Y
+
Y
T
B
T
c
φ
A
Q
-
Q
A
T
+
B
Y
-
Y
T
B
T
c
φ
-
A
Q
+
Q
A
T
-
B
Y
+
Y
T
B
T
s
φ
A
Q
+
Q
A
T
+
B
Y
+
Y
T
B
T
<
0
,where Υ
=
A
i
Q
i
+
Q
i
A
i
T
+
B
u
i
Y
i
+
Y
i
T
B
u
i
T, Y
i
=
K
i
Q
i, and the ∗ term indicates symmetry.Remark 2.
RegionR
1 guarantees that the system is faster than a minimum requirement. Region R
2, in turn, avoids the poles from being allocated too far away from the origin of the plane, avoiding high gains in the controller. Region R
3 limits the system’s maximum percentage overshoot.


### 3.4. State Estimation Using the Linear Kalman Filter with Unknown Inputs

As stated previously, the system has two different sampling times (controller and positioning system). Then, in order to design the state estimator, two state-space equations are formulated distinguishing both sampling times. For the controller’s sampling timeτ
s, we consider a stochastic linear discrete-time dynamic system of the form(58)
x
i
,
k
=
A
x
i
-
1
,
k
+
B
u
i
-
1
,
k
+
G
d
i
-
1
,
k
+
w
i
-
1
,
k
,
y
k
=
C
x
k
+
v
k
,where A
∈
R
n
×
n, B
∈
R
n
×
p, G
∈
R
n
×
s, and C
∈
R
m
×
n are known matrices, in which n is the number of states of the system, p is the number of known inputs, s is the number of unknown inputs, and m is the number of outputs. It is assumed that for all i
,
k
≥
1 the input u
i
-
1
,
k
∈
R
p is known. The process noise w
i
-
1
,
k
∈
R
n is assumed to be white, Gaussian, zero-mean, and mutually independent with known diagonal covariance matrix Q. Besides, the outputs are measured by the sensors, in which the measurement noise v
k
∈
R
m is also assumed to be white, Gaussian, zero-mean, and mutually independent with known diagonal covariance matrix R
k. x
i
,
k
∈
R
n is the state vector, which should be estimated. The initial state vector x
0,0
∈
R
n is assumed to have Gaussian distribution with initial estimate x
^
0,0 and error covariance P
0,0
x
x
≜
E
[
(
x
0
-
x
^
0,0
)
(
x
0
-
x
^
0,0
)
T
], where E
[
·
] denotes expected value.Since the model in (58) uses sampling time τ
s, these equations can be used iteratively to predict the system’s states for sampling time T
s. Thus, in order to predict x
^
k
∣
k
-
1 given x
^
k
-
1
∣
k
-
1 and all the inputs u
i
,
k
-
1 and d
i
,
k
-
1 for i
=
1
,
…
,
T
s
/
τ
s, this prediction can be calculated by finding x
^
i
=
T
s
/
τ
s
,
k. Therefore, defining h
=
T
s
/
τ
s, when estimating the states at an instant k given information up to time k
-
1 the following state-space equations are used:(59)
x
^
k
∣
k
-
1
=
A
T
s
x
^
k
-
1
∣
k
-
1
+
B
T
s
u
→
k
-
1
+
G
T
s
d
→
k
-
1
y
^
k
=
C
x
^
k
∣
k
-
1
,where A
T
s
=
A
h, B
T
s
=
A
h
-
1
B
⋯
B, G
T
s
=
A
h
-
1
G
⋯
G, u
→
k
-
1
=
u
1
,
k
-
1
T
⋯
u
h
,
k
-
1
T
T, and d
→
k
-
1
=
d
1
,
k
-
1
T
⋯
d
h
,
k
-
1
T
T; y
^
k
∈
R
m is the estimated output of the system and is calculated at each time instant k.In order to obtain the system model for state estimation (position and velocity estimation), only the translational dynamics (ξ
¨) are considered from the equations of motion of the tilt-rotor UAV with suspended load (3). Besides, since the generalized coordinates ξ are external variables for the mechanical system (i.e., they do not shape the inertia matrix of the UAV), the Coriolis and centripetal matrix does not depend on these variables and their time-derivatives. Therefore, it can be assumed that all known exerted forces along an axis are grouped into a single component (inertia, Coriolis, and centripetal forces due to existing coupling between the translational dynamics and the remaining generalized coordinates, the gravity force, and the generalized control input forces), compounding the vector T
I
=
T
x
I
T
y
I
T
z
I
T, which is function of the aircraft and load orientation (η, η
˙, α, α
˙, γ, and γ
˙). For instance, on the x-axis it is possible to use the expression x
¨
=
T
x
I
/
m
+
T
d
x
I
/
m, where T
x
I is the known translational force along the axis x expressed in the inertial frame, assuming that all angular positions and velocities are measured with sampling time τ
s, T
d
x
I is the unknown net disturbance force along the same axis, and m is the body’s mass. Note that, for the translational dynamics, the state-space system is linear.Since attitude values are obtained by measurements, it is fair to say thatT
x
I, T
y
I, and T
z
I are not precisely known. Therefore, choosing the state vector x
=
x
y
z
x
˙
y
˙
z
˙
T, with known input vector T
I
=
T
x
I
T
y
I
T
z
I
T and unknown input vector T
d
I
=
T
d
x
I
T
d
y
I
T
d
z
I
T, in which the uncertainties of T
I along with uncertainties of the estimated disturbance force and modeling errors can be expressed by the process noise vector w, the translational dynamics can be written in the stochastic discrete-time domain with the system’s sampling time τ
s through the following state-space representation:(60)
x
i
+
1
=
I
3
×
3
τ
s
I
3
×
3
0
3
×
3
I
3
×
3
x
i
+
τ
s
2
2
m
I
3
×
3
τ
s
m
I
3
×
3
T
i
I
+
τ
s
2
2
m
I
3
×
3
τ
s
m
I
3
×
3
T
d
i
I
+
w
i
,where w
i
=
w
1
w
2
w
3
w
4
w
5
w
6
T is the process noise vector at instant i. Thus, the matrices A, B, and G of (58) are function of the sampling time τ
s and can be extracted from (60). Further, matrices A
T
s, B
T
s, and G
T
s are easily obtained from (60) and, by assumption, C
=
I
6
×
6.In addition, the unknown inputs are estimated once at each instant of the sensors’ measurement and are denoted byd
k
-
1
∣
k. When using (59) to estimate x
^
k
∣
k
-
1, the vector of unknown inputs d
→
k
-
1 possesses its terms equal to d
k
-
1
∣
k for all i
=
1
,
…
,
h. Therefore, G
T
s
d
→
k
-
1
=
∑
j
=
0
h
-
1
A
j
G
d
k
-
1
∣
k, in which(61)
A
j
=
I
3
×
3
j
·
τ
s
I
3
×
3
0
3
×
3
I
3
×
3
,
(62)
A
j
·
G
=
2
j
+
1
τ
s
2
2
m
I
3
×
3
τ
s
m
I
3
×
3
.Thus, ∑
j
=
0
h
-
1
‍
A
j
G can be obtained by calculating the summations ∑
j
=
0
h
-
1
‍
2
j
+
1
τ
s
2
/
2
m
I
3
×
3 and ∑
j
=
0
h
-
1
‍
τ
s
/
m
I
3
×
3 separately. Besides, noting that ∑
j
=
0
h
-
1
‍
(
2
j
+
1
) is the sum of the first h odd numbers, it is given by ∑
j
=
0
h
-
1
‍
2
j
+
1
τ
s
2
/
2
m
=
(
h
τ
s
)
2
/
2
m. Otherwise, since h
=
T
s
/
τ
s, then h
τ
s
=
T
s and ∑
j
=
0
h
-
1
‍
A
j
G
=
G
¯, where G
¯ is equal to G except that τ
s is substituted by T
s. Therefore, it can be said that G
T
s
d
→
k
-
1
=
G
¯
d
k
-
1
∣
k.Now, considering models (58) and (59), the LKFUI will use a three-step algorithm, which is described below.Initialization. Initial estimates are assumed for x
^
0
∣
0, d
^
0
∣
0 and their uncertainties. Every time that the predictor assimilates new measurements, x
^
i
=
0
,
k is initialized as(63)
x
^
i
=
0
,
k
=
x
^
k
∣
k
.Prediction
(i)
Calculatex
^
i
,
k using the controller’s sampling frequency τ
s, for i
=
1,2
,
…
,
T
s
/
τ
s
(64)
x
^
i
,
k
=
A
x
^
i
-
1
∣
k
+
B
u
i
-
1
∣
k
+
G
d
^
k
-
1
∣
k
.
(ii)
Assignx
^
k
+
1
∣
k
=
x
^
i
=
T
s
/
τ
s
,
k.Correction
(i)
Calculatey
^
k
+
1
∣
k given x
^
k
+
1
∣
k:(65)
y
^
k
+
1
∣
k
=
C
x
^
k
+
1
∣
k
.
(ii)
Calculatex
^
k
+
1
∣
k
+
1 given x
^
k
+
1
∣
k, y
k
+
1 and their respective covariance:(66)
x
^
k
+
1
∣
k
+
1
=
x
^
k
+
1
∣
k
+
L
k
+
1
y
k
+
1
-
y
^
k
+
1
∣
k
.
(iii)
Calculate the estimation of the disturbanced
^
k
∣
k
+
1 as follows:(67)
d
^
k
∣
k
+
1
=
G
¯
k
T
G
¯
k
-
1
G
¯
k
T
L
k
+
1
y
k
+
1
-
y
^
k
+
1
∣
k
.
(iv)
Incrementk and then go back to thePrediction step.The filter gainL
k
∈
R
n
×
m must be chosen such that it minimizes the cost function J
k
(
L
k
)
≜
E
[
(
x
k
-
x
^
k
∣
k
)
T
(
x
k
-
x
^
k
∣
k
)
] subject to the constraint L
k
E
k
=
F
k, where E
k and F
k solve the problem of the LKFUI when they are given by E
k
=
C
G
¯ and F
k
=
G
¯ [23].Consider theforecast error
e
k
∣
k
-
1
≜
x
k
-
x
^
k
∣
k
-
1, theinnovation
ν
k
∣
k
-
1
≜
y
k
-
y
^
k
∣
k
-
1, thedata-assimilation error
e
k
∣
k
≜
x
k
-
x
^
k
∣
k, theforecast error covariance
P
k
∣
k
-
1
x
x
≜
E
[
e
k
∣
k
-
1
e
k
∣
k
-
1
T
], theinnovation covariance
P
k
∣
k
-
1
y
y
≜
E
[
ν
k
∣
k
-
1
ν
k
∣
k
-
1
T
], and thecross covariance
P
k
∣
k
-
1
x
y
≜
E
[
e
k
∣
k
-
1
ν
k
∣
k
-
1
T
]. From the equations of filter presented in the LKFUI algorithm, thedata-assimilation error covariance
P
k
∣
k
x
x is given by [23](68)
P
k
∣
k
x
x
≜
E
e
k
∣
k
e
k
∣
k
T
=
P
k
∣
k
-
1
x
x
-
L
k
P
k
∣
k
-
1
x
y
T
-
P
k
∣
k
-
1
x
y
L
k
T
+
L
k
P
k
∣
k
-
1
y
y
L
k
T
,where P
k
∣
k
-
1
x
x
=
A
T
s
P
k
-
1
∣
k
-
1
x
x
A
T
s
T
+
Q, P
k
∣
k
-
1
y
y
=
C
P
k
∣
k
-
1
x
x
C
T
+
R
k, and P
k
∣
k
-
1
x
y
=
P
k
∣
k
-
1
x
x
C
T.For convenience, the following terms are defined asE
k
L
≜
(
E
k
T
E
K
)
-
1
E
k
T, Ω
k
⊥
≜
I
m
×
m
-
Ω
k, K
k
≜
P
k
∣
k
-
1
x
y
(
P
k
∣
k
-
1
y
y
)
-
1, and Ω
k
≜
E
k
[
E
k
T
(
P
k
∣
k
-
1
y
y
)
-
1
E
k
]
-
1
E
k
T
(
P
k
∣
k
-
1
y
y
)
-
1.The gainL
k that solves the optimization problem is given by [23](69)
L
k
=
K
k
Ω
k
⊥
+
F
k
E
k
L
Ω
k
.Thedata-assimilation error covariance is then given by the Riccati equation (70)
P
k
∣
k
x
x
=
P
k
∣
k
-
1
x
x
-
P
k
∣
k
-
1
x
y
P
k
∣
k
-
1
y
y
-
1
P
k
∣
k
-
1
x
y
T
+
F
k
E
k
L
Ω
k
P
k
∣
k
-
1
y
y
F
k
E
k
L
Ω
k
T
-
Δ
1
,
k
+
Δ
1
,
k
T
+
P
k
∣
k
-
1
x
y
P
k
∣
k
-
1
y
y
-
1
Ω
k
P
k
∣
k
-
1
y
y
P
k
∣
k
-
1
x
y
P
k
∣
k
-
1
y
y
-
1
Ω
k
T
,where Δ
1
,
k
≜
P
k
∣
k
-
1
x
y
(
P
k
∣
k
-
1
y
y
)
-
1
Ω
k
P
k
∣
k
-
1
y
y
Ω
k
T
(
F
k
E
k
L
)
T.


## 3.1. Inner-Loop Control Design

In this section the proposed IOFL control law is designed for controlling system (25) considering the systems’ outputs h
q
1
x
q
1
=
z
ϕ
θ
ψ
T and the system’s inputs Γ
=
f
R
f
L
τ
α
R
τ
α
L
T.Recall that the relative degreer
i of a system’s output h
i is exactly the number of times one has to differentiate h
i in order to have at least one component of the vector input Γ explicitly appearing. It is computed by r
i
=
(
inf
⁡
k
,
∃
j
,
1
≤
j
≤
4
,
L
g
j
L
f
k
-
1
h
i
), where L
f
h
i is the Lie derivative operation and g
j is the j
t
h column of g
u
q
1
(
x
q
i
). It is possible to verify that r
i
=
2
,
∀
i
∈
{
1,2
,
3,4
} (this statement can be trivially understood due to the fact that z
¨, ϕ
¨, θ
¨, and ψ
¨ are all function of at least one input Γ
i for arbitrary α
R and α
L), which implies that this set of inputs/outputs are not fully feedback linearizable, since r
=
∑
‍
r
i
=
8 and the system has twelve state variables (q
1 and q
˙
1).Since it is only possible to regulate on a desired value the number of outputs equal to the number of control inputs, and the time-derivativesz
¨, ϕ
¨, θ
¨, and ψ
¨ are all affected only by the control inputs f
R and f
L (see (23)), and none by the torques τ
R and τ
L, it is required to delay the appearance of f
R and f
L to higher order time-derivatives of z, ϕ, θ, and ψ in order to allow the other control inputs to appear. Therefore, the present problem can be solved by using the dynamic extension technique, in which the state vector of the system is extended by setting f
R and f
L equal to the output of two double integrators and, then, the state variables x
13
=
f
R, x
14
=
f
˙
R, x
15
=
f
L, and x
16
=
f
˙
L are included. The new system’s inputs are defined as Γ
¯
=
f
¨
R
f
¨
L
τ
α
R
τ
α
L
T. Thus, the new state vector is x
¯
q
1
=
q
1
T
q
˙
1
T
x
13
x
14
x
15
x
16
T, and the state-space equations can be rewritten as(29)
x
¯
˙
q
1
=
f
¯
q
1
x
¯
q
1
+
∑
i
=
1
4
g
¯
i
x
¯
q
1
Γ
¯
i
+
g
¯
d
q
1
x
¯
q
1
δ
q
1
,
y
F
L
=
h
q
1
x
q
1
=
z
ϕ
θ
ψ
T
,where f
¯
q
1
x
¯
q
1
=
q
˙
1
T
q
¨
1
T
x
14
0
x
16
0
T, q
¨
1
=
M
q
1
-
1
(
-
[
C
q
1
+
μ
q
1
]
-
G
q
1
)
+
B
1
x
13
+
B
2
x
15, Γ
¯
i is the i
t
h element of the system’s input vector Γ
¯, g
¯
i
(
x
¯
q
1
) is given by(30)
g
¯
1
=
0
13
×
1
1
0
2
×
1
,
(31)
g
¯
2
=
0
15
×
1
1
,
(32)
g
¯
3
=
0
6
×
1
B
3
0
4
×
1
,
(33)
g
¯
4
=
0
6
×
1
B
4
0
4
×
1
,and B
i is the i
t
h column of B
q
1, i
∈
{
1,2
,
3,4
}. Note that the control inputs f
R and f
L of system (24) are not anymore control inputs for the extended system (29), being now state variables. Additionally, it improves the capability of the inner-loop controller to deal with external disturbances.Computing the relative degreer
i for the new system, it is now possible to verify that r
i
=
4
,
∀
i
∈
{
1,2
,
3,4
}. This only holds provided the condition that the terms of the inertia matrix m
23
=
m
24
=
0. If this was not true, then the coupling between η and α would imply that it is physically possible to actuate on attitude angles by Newton’s action-reaction when applying torques to rotate the propellers. Even though this is mathematically attainable, this is not the form of actuation desired for the conception of a tilt-rotor aircraft. Therefore, the system’s relative degree is r
=
∑
‍
r
i
=
16, which means that the system is now fully feedback linearizable. Thus, in order to apply IOFL techniques, the following needs to be calculated:(34)
Δ
x
¯
q
1
∈
R
4
×
4
:
Δ
i
j
x
¯
q
1
=
L
g
¯
j
L
f
¯
3
h
i
x
¯
q
1
,
(35)
b
x
¯
q
1
∈
R
4
×
1
:
b
i
x
¯
q
1
=
L
f
¯
4
h
i
x
¯
q
1
,where Δ
i
j
(
x
¯
q
1
) is the element of the i
t
h row and j
t
h column of Δ
(
x
¯
q
1
) and b
i
(
x
¯
q
1
) is the i
t
h row of b
(
x
¯
q
1
).Then, the following control law may be applied(36)
Γ
¯
=
Δ
-
1
v
q
1
-
b
,where v
q
1
=
v
z
v
ϕ
v
θ
v
ψ
T is the additional control input. The matrix Δ
(
x
¯
q
1
) has full-rank on the neighborhood of the system’s equilibrium point, implying that(37)
h
…
.
q
1
x
¯
q
1
=
v
q
1
+
π
1
δ
q
1
,with π
1
(
δ
q
1
) a function of the residual terms due to model simplifications, unmodeled dynamics, and unknown external disturbancesTherefore, it is then possible to apply linear techniques to obtain the additional control law in order to regulate the outputsz, ϕ, θ, and ψ. These variables are controlled by Proportional-Integral-Derivative- (PID-) like controllers with feed-forward term. This choice, together with a desired faster closed-loop time response, allows assuming the steady-state behavior of the inner-loop as unitary gain at the moment of tuning the outer-loop controllers. Moreover, the variables z, ϕ, θ, and ψ are required to track desired trajectories when constant external disturbances and modeling errors affect the system. Thus, the linear control laws compounding the additional control input vector, v
q
1, are designed as(38)
v
z
=
z
…
.
R
e
f
+
K
d
d
d
z
e
⃛
z
+
K
d
d
z
e
¨
z
+
K
d
z
e
˙
z
+
K
p
z
e
z
+
K
I
z
∫
e
z
d
t
,
(39)
v
ϕ
=
ϕ
…
.
R
e
f
+
K
d
d
d
ϕ
e
⃛
ϕ
+
K
d
d
ϕ
e
¨
ϕ
+
K
d
ϕ
e
˙
ϕ
+
K
p
ϕ
e
ϕ
+
K
I
ϕ
∫
e
ϕ
d
t
,
(40)
v
θ
=
θ
…
.
R
e
f
+
K
d
d
d
θ
e
⃛
θ
+
K
d
d
θ
e
¨
θ
+
K
d
θ
e
˙
θ
+
K
p
θ
e
θ
+
K
I
θ
∫
e
θ
d
t
,
(41)
v
ψ
=
ψ
…
.
R
e
f
+
K
d
d
d
ψ
e
⃛
ψ
+
K
d
d
ψ
e
¨
ψ
+
K
d
ψ
e
˙
ψ
+
K
p
ψ
e
ψ
+
K
I
ψ
∫
e
ψ
d
t
,where e
i
=
(
i
-
i
R
e
f
) with i
∈
{
z
,
ϕ
,
θ
,
ψ
}, and z
…
.
R
e
f, ϕ
…
.
R
e
f, θ
…
.
R
e
f, and ψ
…
.
R
e
f are feed-forward terms due to the time-varying reference.The above control laws will be synthesized using the mixedH
2
/
H
∞ control technique with pole placement constraints in order to ensure reference tracking for the variables z, ϕ, θ, and ψ, while α
R and α
L are stabilized.


## 3.2. Outer-Loop Control Design

The outer-loop IOFL controller is designed in order to control the subsystem (26) by actuating on the references of roll, ϕ
R
e
f, and pitch, θ
R
e
f, angles. Moreover, an additional control term is considered to reduce the load swing.From (23), the generalized control input vector of subsystem (26) can be rewritten as B
q
2
Γ
=
T
x
I
T
y
I
0
0
T. The variables T
x
I and T
y
I are, respectively, the translational forces along x and y expressed in the inertial frame and are obtained from(42)
T
x
I
T
y
I
=
i
→
T
j
→
T
R
B
I
F
R
B
+
F
L
B
=
r
11
r
12
r
13
r
21
r
22
r
23
F
R
B
+
F
L
B
,where r
i
j is the element of R
B
I located at the i
t
h line and j
t
h column. From F
R
B
+
F
L
B
=
f
x
B
f
y
B
f
z
B
T, the actuation approach in this part of the design is performed by changing the projection of f
z
B along x
I and y
I. In other words, the roll (ϕ) and pitch (θ) angles need to be changed so as to obtain the desired projections of f
z
B along x
I and y
I. Besides, the projections of f
x
B and f
y
B along the inertial frame are assumed as known disturbances which should be compensated by the controller.First, in order to perform the feedback linearization approach and design an additional linear control input, from (26) and (42), the vector of dynamics to be canceled n
q
2
=
n
x
n
y
n
γ
1
n
γ
2
T is defined as(43)
n
q
2
=
C
q
2
+
μ
q
2
q
˙
2
+
G
q
2
-
F
q
2
′
,where F
q
2
′
=
r
11
f
x
B
+
r
12
f
y
B
r
21
f
x
B
+
r
22
f
y
B
0
0
T. Consequently, subsystem (26) can be rewritten as(44)
M
q
2
q
¨
2
+
n
x
n
y
n
γ
1
n
γ
2
︸
n
q
2
=
s
ψ
s
ϕ
+
c
ψ
s
θ
c
ϕ
f
z
B
s
ψ
s
θ
c
ϕ
-
c
ψ
s
ϕ
f
z
B
0
0
+
δ
q
2
.By inspecting (44), it is not trivial to obtain a nonlinear state-space representation affine in the control inputs ϕ
R
e
f and θ
R
e
f. To solve this problem, which will allow applying the IOFL technique, and assuming that vector of dynamics n
q
2 is known, one can define the intermediary input values Υ
1 and Υ
2 and state the following nonlinear set of equations:(45)
s
ψ
s
ϕ
+
c
ψ
s
θ
c
ϕ
f
z
B
-
n
x
=
Υ
1
,
s
ψ
s
θ
c
ϕ
-
c
ψ
s
ϕ
f
z
B
-
n
y
=
Υ
2
.As the term f
z
B is obtained from the inner-loop controller, which has a faster settling time than the outer one, for the outer-loop control design, this variable is assumed as a time-varying parameter. Furthermore, since the roll (ϕ) and pitch (θ) angles are controlled variables of the inner-loop, they can be written with respect to their error and reference as e
ϕ
=
(
ϕ
-
ϕ
R
e
f
) and e
θ
=
(
θ
-
θ
R
e
f
). Therefore, being the inner-loop controller in charge of carrying these errors to zero, and assuming that the yaw (ψ) angle is measurable, (45) can be written with respect to ϕ
R
e
f and θ
R
e
f, whose solutions are given by(46)
ϕ
R
e
f
=
arcsin
⁡
σ
s
ψ
Υ
1
+
n
x
-
c
ψ
Υ
2
+
n
y
f
z
B
,
θ
R
e
f
=
arcsin
⁡
σ
c
ψ
Υ
1
+
n
x
+
s
ψ
Υ
2
+
n
y
f
z
B
c
ϕ
R
e
f
,where σ
(
p
) is a saturation function that avoids the a
r
c
s
i
n operation from becoming undefined and is given by σ
(
p
)
=
min
⁡
1
,
max
⁡
-
1
,
p (the function min
⁡
a
,
b returns the minimum value between a and b and max
⁡
c
,
d returns the maximum value between c and d).By replacing the transformation (45) in (44), (26) is reduced to(47)
M
q
2
q
¨
2
+
n
¯
q
2
=
Υ
q
2
+
δ
q
2
,where n
¯
q
2
=
0
0
n
γ
1
n
γ
2
T and Υ
q
2
=
Υ
1
Υ
2
0
0
T. Since (47) is affine in the intermediary input vector Υ
q
2, IOFL techniques can be used to linearize this subsystem.Then, by choosing the outputsh
q
2
x
q
2
=
x
y
T, matrices Δ
(
x
q
2
) and b
(
x
q
2
) can be computed so as to apply a control law similar to the one in (36) leading to(48)
h
¨
q
2
x
q
2
=
v
q
2
+
π
2
δ
q
2
,where v
q
2
=
v
x
v
y
T. Consequently, it is possible to control the aircraft’s motion along x
I and y
I, while γ
1 and γ
2 compose the internal dynamics.Therefore, linear controllers are designed taking into account the motion of the load expressed in the inertial frame. First, the vectord
B
C
4
I
=
R
B
I
d
4
B
=
d
x
I
d
y
I
d
z
I
T is defined as the translation between frame B and the load expressed in frame I. Then, the reference for the load variation is given as d
˙
x
I
d
˙
x
I
T
=
0
0
T. Note that d
4
B is the vector that goes from the origin of frame B to the origin of frame C
4. Hence, d
˙
4
B is the relative velocity between the origins of frames B and C
4. That is, the objective is to avoid load swinging speed regardless of its position with respect to the aircraft. In order to accomplish this task, PID-like controllers for the x
y motion with additional PI-like control terms to regulate the load’s speed with respect to the tilt-rotor are applied:(49)
v
x
=
x
¨
R
e
f
+
K
d
x
e
˙
x
+
K
p
x
e
x
+
K
I
x
∫
e
x
d
t
︸
x
-
m
o
t
i
o
n
+
K
p
d
x
d
˙
x
I
+
K
I
d
x
∫
d
˙
x
I
d
t
︸
l
o
a
d
s
w
i
n
g
,
(50)
v
y
=
y
¨
R
e
f
+
K
d
y
e
˙
y
+
K
p
y
e
y
+
K
I
y
∫
e
y
d
t
︸
y
-
m
o
t
i
o
n
+
K
p
d
y
d
˙
y
I
+
K
I
d
y
∫
d
˙
y
I
d
t
︸
l
o
a
d
s
w
i
n
g
,where x
¨
R
e
f and y
¨
R
e
f are feed-forward terms due to the time-varying reference.When tuning parameters of controllers (49) and (50), it is again important to bear in mind that this strategy assumed null coupling between x, y, γ
1, and γ
2 and the remaining variables, which is incorporated into the residual terms, π
2
(
δ
q
2
). Again, in order to synthesize the controllers’ gains providing good transient response and disturbance attenuation, the mixed H
2
/
H
∞ synthesis with pole placement constraints is employed.


## 3.3. Robust MixedH
2
/
H
∞ Control Design

After applying the input-output feedback linearization procedure to both inner and outer loops, the additional control inputs (38)–(41) and (49)-(50) are synthesized using the mixed H
2
/
H
∞ technique with pole placement constraints [11]. These linearized dynamics can be written as(51)
e
˙
i
=
A
i
e
i
+
B
u
i
u
~
i
+
B
π
i
π
i
,
z
i
=
C
z
i
e
i
+
D
u
z
i
u
~
i
+
D
π
z
i
π
i
,
u
~
i
=
K
i
e
i
,where e
i
=
[
∫
‍
e
i
d
t
e
i
e
˙
i
e
¨
i
e
⃛
i
]
T and K
i
=
[
K
I
i
K
p
i
K
d
i
K
d
d
i
K
d
d
d
i
] with i
∈
{
z
,
ϕ
,
θ
,
ψ
} for the inner-loop, e
i
=
[
∫
‍
e
i
d
t
e
i
e
˙
i
]
T and K
i
=
[
K
I
i
K
p
i
K
d
i
] with i
∈
{
x
,
y
} for the x
y-motion, and e
i
=
[
∫
‍
d
˙
i
I
d
t
d
˙
i
I
]
T and K
i
=
[
K
I
d
i
K
p
d
i
] with i
∈
{
γ
1
,
γ
2
} for the term in charge of load’s swing. Matrices A
i, B
u
i, B
π
i are defined with proper dimensions, and C
z
i, D
u
z
i, D
π
z
i are weighting matrices.The optimalH
∞ controller asymptotically stabilizes the linear system and minimizes the H
∞ norm, H
π
i
z
i
(
s
)
∞
=
s
u
p
(
σ
H
π
i
z
i
(
s
)
), where H
π
i
z
i
(
s
) is the transfer function between the external disturbance π
i and the cost variable z
i. Applying Parseval’s theorem to the H
∞ norm, defined in the frequency domain, it results in z
i
2
≤
H
π
i
z
i
s
∞
π
i
2, in which minimizing the H
∞ norm minimizes the system disturbance effect. Solving the H
∞ optimization problem, the linear system is asymptotically stable with H
∞ norm bounded by H
π
i
z
i
(
s
)
∞
<
γ
∗, where γ
∗ is the optimal value.Considering again system (51) with D
π
z
i
=
0, H
2 control problem leads to an optimal controller that minimizes the H
2 norm h
i
(
t
)
2
=
∫
0
∞
‍
h
i
2
(
t
)
d
t for h
i
(
t
)
=
L
-
1
{
H
π
i
z
i
(
s
)
}
. The minimization of this norm implicates in the optimal transient response of the system. Therefore, the mixed H
2
/
H
∞ control approach is formulated by finding an optimal H
∞ norm and using γ
^
>
γ
∗ to solve the H
2 problem maintaining the H
∞ constraints, providing the advantage of a small H
∞ norm and a better transient response.In addition, to ensure pole placement constraints, the closed-loop poles must be within a defined convex region in the complex plane (D-stability). An LMI region is defined as a subset D of the complex plane C, where D is given by D
=
s
∈
C
:
L
+
s
M
+
s
∗
M
T
<
0, where L
=
L
T and M are matrices that define region D, s
=
σ
+
j
ω, and s
∗ is the complex conjugate of s [29]. In this work the following regions are used: (R
1) R
e
(
s
)
<
-
α, which is defined by s
+
s
∗
<
-
2
α, given that s
+
s
∗
=
σ
+
j
ω
+
σ
-
j
ω
=
2
σ. Therefore, it is possible to rewrite this expression in the form 2
α
+
s
+
s
∗
<
0 and then L
=
2
α and M
=
1; (R
2) is a disk centered at (
c
,
0
) with radius r, in which L
=
-
r
-
c
-
c
-
r and M
=
0
1
0
0; and (R
3) is a cone defined by I
m
s
<
t
a
n
φ
R
e
s, from which the matrices L and M are given as L
=
0
0
0
0 and M
=
s
i
n
(
φ
)
c
o
s
(
φ
)
-
c
o
s
(
φ
)
s
i
n
(
φ
). According to [29], a closed-loop state matrix A
+
B
K is D-stable if and only if there exists a real definite positive symmetric matrix Q
∈
R
n
×
n such that L
⊗
Q
+
M
⊗
V
+
M
T
⊗
V
T
<
0, where ⊗ is the Kronecker product operator, and V
=
(
A
Q
+
B
Y
) with K
=
Y
Q
-
1.Thus, the robust mixedH
2
/
H
∞ control problem with pole placement constraints is formulated for each linearized system (51) through the following optimization problem subject to linear matrix inequalities assuming γ
=
γ
^:(52)
m
i
n
t
r
N
i
,
s
.
t
.
Q
i
>
0
(53)
H
2
c
o
n
s
t
r
a
i
n
t
s
N
i
C
z
i
+
D
u
z
i
K
i
Q
i
Q
i
C
z
i
+
D
u
z
i
K
i
T
Q
i
>
0
,
Q
i
A
i
T
+
A
i
Q
i
+
Y
i
T
B
u
i
T
+
B
u
i
Y
i
+
B
π
i
B
π
i
T
<
0
,
(54)
H
∞
c
o
n
s
t
r
a
i
n
t
s
Υ
B
π
i
Q
i
C
z
i
T
+
Y
i
T
D
u
z
i
T
∗
-
γ
I
n
w
i
D
π
z
i
T
∗
∗
-
γ
I
n
z
i
<
0
,
(55)
R
e
g
i
o
n
R
1
⟶
2
α
Q
+
A
Q
+
Q
A
T
+
B
Y
+
Y
T
B
T
<
0
,
(56)
R
e
g
i
o
n
R
2
⟶
-
r
Q
-
c
Q
+
A
Q
+
B
Y
-
c
Q
+
Q
A
T
+
Y
T
B
T
-
r
Q
<
0
,
(57)
R
e
g
i
o
n
R
3
⟶
s
φ
A
Q
+
Q
A
T
+
B
Y
+
Y
T
B
T
c
φ
A
Q
-
Q
A
T
+
B
Y
-
Y
T
B
T
c
φ
-
A
Q
+
Q
A
T
-
B
Y
+
Y
T
B
T
s
φ
A
Q
+
Q
A
T
+
B
Y
+
Y
T
B
T
<
0
,where Υ
=
A
i
Q
i
+
Q
i
A
i
T
+
B
u
i
Y
i
+
Y
i
T
B
u
i
T, Y
i
=
K
i
Q
i, and the ∗ term indicates symmetry.Remark 2.
RegionR
1 guarantees that the system is faster than a minimum requirement. Region R
2, in turn, avoids the poles from being allocated too far away from the origin of the plane, avoiding high gains in the controller. Region R
3 limits the system’s maximum percentage overshoot.


## 3.4. State Estimation Using the Linear Kalman Filter with Unknown Inputs

As stated previously, the system has two different sampling times (controller and positioning system). Then, in order to design the state estimator, two state-space equations are formulated distinguishing both sampling times. For the controller’s sampling timeτ
s, we consider a stochastic linear discrete-time dynamic system of the form(58)
x
i
,
k
=
A
x
i
-
1
,
k
+
B
u
i
-
1
,
k
+
G
d
i
-
1
,
k
+
w
i
-
1
,
k
,
y
k
=
C
x
k
+
v
k
,where A
∈
R
n
×
n, B
∈
R
n
×
p, G
∈
R
n
×
s, and C
∈
R
m
×
n are known matrices, in which n is the number of states of the system, p is the number of known inputs, s is the number of unknown inputs, and m is the number of outputs. It is assumed that for all i
,
k
≥
1 the input u
i
-
1
,
k
∈
R
p is known. The process noise w
i
-
1
,
k
∈
R
n is assumed to be white, Gaussian, zero-mean, and mutually independent with known diagonal covariance matrix Q. Besides, the outputs are measured by the sensors, in which the measurement noise v
k
∈
R
m is also assumed to be white, Gaussian, zero-mean, and mutually independent with known diagonal covariance matrix R
k. x
i
,
k
∈
R
n is the state vector, which should be estimated. The initial state vector x
0,0
∈
R
n is assumed to have Gaussian distribution with initial estimate x
^
0,0 and error covariance P
0,0
x
x
≜
E
[
(
x
0
-
x
^
0,0
)
(
x
0
-
x
^
0,0
)
T
], where E
[
·
] denotes expected value.Since the model in (58) uses sampling time τ
s, these equations can be used iteratively to predict the system’s states for sampling time T
s. Thus, in order to predict x
^
k
∣
k
-
1 given x
^
k
-
1
∣
k
-
1 and all the inputs u
i
,
k
-
1 and d
i
,
k
-
1 for i
=
1
,
…
,
T
s
/
τ
s, this prediction can be calculated by finding x
^
i
=
T
s
/
τ
s
,
k. Therefore, defining h
=
T
s
/
τ
s, when estimating the states at an instant k given information up to time k
-
1 the following state-space equations are used:(59)
x
^
k
∣
k
-
1
=
A
T
s
x
^
k
-
1
∣
k
-
1
+
B
T
s
u
→
k
-
1
+
G
T
s
d
→
k
-
1
y
^
k
=
C
x
^
k
∣
k
-
1
,where A
T
s
=
A
h, B
T
s
=
A
h
-
1
B
⋯
B, G
T
s
=
A
h
-
1
G
⋯
G, u
→
k
-
1
=
u
1
,
k
-
1
T
⋯
u
h
,
k
-
1
T
T, and d
→
k
-
1
=
d
1
,
k
-
1
T
⋯
d
h
,
k
-
1
T
T; y
^
k
∈
R
m is the estimated output of the system and is calculated at each time instant k.In order to obtain the system model for state estimation (position and velocity estimation), only the translational dynamics (ξ
¨) are considered from the equations of motion of the tilt-rotor UAV with suspended load (3). Besides, since the generalized coordinates ξ are external variables for the mechanical system (i.e., they do not shape the inertia matrix of the UAV), the Coriolis and centripetal matrix does not depend on these variables and their time-derivatives. Therefore, it can be assumed that all known exerted forces along an axis are grouped into a single component (inertia, Coriolis, and centripetal forces due to existing coupling between the translational dynamics and the remaining generalized coordinates, the gravity force, and the generalized control input forces), compounding the vector T
I
=
T
x
I
T
y
I
T
z
I
T, which is function of the aircraft and load orientation (η, η
˙, α, α
˙, γ, and γ
˙). For instance, on the x-axis it is possible to use the expression x
¨
=
T
x
I
/
m
+
T
d
x
I
/
m, where T
x
I is the known translational force along the axis x expressed in the inertial frame, assuming that all angular positions and velocities are measured with sampling time τ
s, T
d
x
I is the unknown net disturbance force along the same axis, and m is the body’s mass. Note that, for the translational dynamics, the state-space system is linear.Since attitude values are obtained by measurements, it is fair to say thatT
x
I, T
y
I, and T
z
I are not precisely known. Therefore, choosing the state vector x
=
x
y
z
x
˙
y
˙
z
˙
T, with known input vector T
I
=
T
x
I
T
y
I
T
z
I
T and unknown input vector T
d
I
=
T
d
x
I
T
d
y
I
T
d
z
I
T, in which the uncertainties of T
I along with uncertainties of the estimated disturbance force and modeling errors can be expressed by the process noise vector w, the translational dynamics can be written in the stochastic discrete-time domain with the system’s sampling time τ
s through the following state-space representation:(60)
x
i
+
1
=
I
3
×
3
τ
s
I
3
×
3
0
3
×
3
I
3
×
3
x
i
+
τ
s
2
2
m
I
3
×
3
τ
s
m
I
3
×
3
T
i
I
+
τ
s
2
2
m
I
3
×
3
τ
s
m
I
3
×
3
T
d
i
I
+
w
i
,where w
i
=
w
1
w
2
w
3
w
4
w
5
w
6
T is the process noise vector at instant i. Thus, the matrices A, B, and G of (58) are function of the sampling time τ
s and can be extracted from (60). Further, matrices A
T
s, B
T
s, and G
T
s are easily obtained from (60) and, by assumption, C
=
I
6
×
6.In addition, the unknown inputs are estimated once at each instant of the sensors’ measurement and are denoted byd
k
-
1
∣
k. When using (59) to estimate x
^
k
∣
k
-
1, the vector of unknown inputs d
→
k
-
1 possesses its terms equal to d
k
-
1
∣
k for all i
=
1
,
…
,
h. Therefore, G
T
s
d
→
k
-
1
=
∑
j
=
0
h
-
1
A
j
G
d
k
-
1
∣
k, in which(61)
A
j
=
I
3
×
3
j
·
τ
s
I
3
×
3
0
3
×
3
I
3
×
3
,
(62)
A
j
·
G
=
2
j
+
1
τ
s
2
2
m
I
3
×
3
τ
s
m
I
3
×
3
.Thus, ∑
j
=
0
h
-
1
‍
A
j
G can be obtained by calculating the summations ∑
j
=
0
h
-
1
‍
2
j
+
1
τ
s
2
/
2
m
I
3
×
3 and ∑
j
=
0
h
-
1
‍
τ
s
/
m
I
3
×
3 separately. Besides, noting that ∑
j
=
0
h
-
1
‍
(
2
j
+
1
) is the sum of the first h odd numbers, it is given by ∑
j
=
0
h
-
1
‍
2
j
+
1
τ
s
2
/
2
m
=
(
h
τ
s
)
2
/
2
m. Otherwise, since h
=
T
s
/
τ
s, then h
τ
s
=
T
s and ∑
j
=
0
h
-
1
‍
A
j
G
=
G
¯, where G
¯ is equal to G except that τ
s is substituted by T
s. Therefore, it can be said that G
T
s
d
→
k
-
1
=
G
¯
d
k
-
1
∣
k.Now, considering models (58) and (59), the LKFUI will use a three-step algorithm, which is described below.Initialization. Initial estimates are assumed for x
^
0
∣
0, d
^
0
∣
0 and their uncertainties. Every time that the predictor assimilates new measurements, x
^
i
=
0
,
k is initialized as(63)
x
^
i
=
0
,
k
=
x
^
k
∣
k
.Prediction
(i)
Calculatex
^
i
,
k using the controller’s sampling frequency τ
s, for i
=
1,2
,
…
,
T
s
/
τ
s
(64)
x
^
i
,
k
=
A
x
^
i
-
1
∣
k
+
B
u
i
-
1
∣
k
+
G
d
^
k
-
1
∣
k
.
(ii)
Assignx
^
k
+
1
∣
k
=
x
^
i
=
T
s
/
τ
s
,
k.Correction
(i)
Calculatey
^
k
+
1
∣
k given x
^
k
+
1
∣
k:(65)
y
^
k
+
1
∣
k
=
C
x
^
k
+
1
∣
k
.
(ii)
Calculatex
^
k
+
1
∣
k
+
1 given x
^
k
+
1
∣
k, y
k
+
1 and their respective covariance:(66)
x
^
k
+
1
∣
k
+
1
=
x
^
k
+
1
∣
k
+
L
k
+
1
y
k
+
1
-
y
^
k
+
1
∣
k
.
(iii)
Calculate the estimation of the disturbanced
^
k
∣
k
+
1 as follows:(67)
d
^
k
∣
k
+
1
=
G
¯
k
T
G
¯
k
-
1
G
¯
k
T
L
k
+
1
y
k
+
1
-
y
^
k
+
1
∣
k
.
(iv)
Incrementk and then go back to thePrediction step.The filter gainL
k
∈
R
n
×
m must be chosen such that it minimizes the cost function J
k
(
L
k
)
≜
E
[
(
x
k
-
x
^
k
∣
k
)
T
(
x
k
-
x
^
k
∣
k
)
] subject to the constraint L
k
E
k
=
F
k, where E
k and F
k solve the problem of the LKFUI when they are given by E
k
=
C
G
¯ and F
k
=
G
¯ [23].Consider theforecast error
e
k
∣
k
-
1
≜
x
k
-
x
^
k
∣
k
-
1, theinnovation
ν
k
∣
k
-
1
≜
y
k
-
y
^
k
∣
k
-
1, thedata-assimilation error
e
k
∣
k
≜
x
k
-
x
^
k
∣
k, theforecast error covariance
P
k
∣
k
-
1
x
x
≜
E
[
e
k
∣
k
-
1
e
k
∣
k
-
1
T
], theinnovation covariance
P
k
∣
k
-
1
y
y
≜
E
[
ν
k
∣
k
-
1
ν
k
∣
k
-
1
T
], and thecross covariance
P
k
∣
k
-
1
x
y
≜
E
[
e
k
∣
k
-
1
ν
k
∣
k
-
1
T
]. From the equations of filter presented in the LKFUI algorithm, thedata-assimilation error covariance
P
k
∣
k
x
x is given by [23](68)
P
k
∣
k
x
x
≜
E
e
k
∣
k
e
k
∣
k
T
=
P
k
∣
k
-
1
x
x
-
L
k
P
k
∣
k
-
1
x
y
T
-
P
k
∣
k
-
1
x
y
L
k
T
+
L
k
P
k
∣
k
-
1
y
y
L
k
T
,where P
k
∣
k
-
1
x
x
=
A
T
s
P
k
-
1
∣
k
-
1
x
x
A
T
s
T
+
Q, P
k
∣
k
-
1
y
y
=
C
P
k
∣
k
-
1
x
x
C
T
+
R
k, and P
k
∣
k
-
1
x
y
=
P
k
∣
k
-
1
x
x
C
T.For convenience, the following terms are defined asE
k
L
≜
(
E
k
T
E
K
)
-
1
E
k
T, Ω
k
⊥
≜
I
m
×
m
-
Ω
k, K
k
≜
P
k
∣
k
-
1
x
y
(
P
k
∣
k
-
1
y
y
)
-
1, and Ω
k
≜
E
k
[
E
k
T
(
P
k
∣
k
-
1
y
y
)
-
1
E
k
]
-
1
E
k
T
(
P
k
∣
k
-
1
y
y
)
-
1.The gainL
k that solves the optimization problem is given by [23](69)
L
k
=
K
k
Ω
k
⊥
+
F
k
E
k
L
Ω
k
.Thedata-assimilation error covariance is then given by the Riccati equation (70)
P
k
∣
k
x
x
=
P
k
∣
k
-
1
x
x
-
P
k
∣
k
-
1
x
y
P
k
∣
k
-
1
y
y
-
1
P
k
∣
k
-
1
x
y
T
+
F
k
E
k
L
Ω
k
P
k
∣
k
-
1
y
y
F
k
E
k
L
Ω
k
T
-
Δ
1
,
k
+
Δ
1
,
k
T
+
P
k
∣
k
-
1
x
y
P
k
∣
k
-
1
y
y
-
1
Ω
k
P
k
∣
k
-
1
y
y
P
k
∣
k
-
1
x
y
P
k
∣
k
-
1
y
y
-
1
Ω
k
T
,where Δ
1
,
k
≜
P
k
∣
k
-
1
x
y
(
P
k
∣
k
-
1
y
y
)
-
1
Ω
k
P
k
∣
k
-
1
y
y
Ω
k
T
(
F
k
E
k
L
)
T.


## 4. Simulations Results and Analysis

This section shows simulation results carried out to analyze the performance of the proposed nonlinear control strategy with swing-free load and two control loops (NLLS-2L), which was designed for the tilt-rotor UAV with suspended load in the presence of unknown disturbances, parametric uncertainties, unmodeled dynamics, and noisy measurements. Simulations were performed using the UAV model parameters of Table1. The simulations considered that the model’s masses m
i and inertia tensors I
i, for i
=
1,2
,
3,4, had all uncertainties ranging from -
30
% to 30
% of their nominal values. The thrusters’ inputs were bounded to work between the range from 0 N to 17 N, and the rotors’ input torques were bounded between -
1 N·m and 1 N·m (these values were chosen so as to match commercial specifications for those actuators).Table 1
System parameters.
Variable
Value
Variable
Value (Kg⋅m2)
β
5∘
I
1
x
x
0.01902947
l
0.5 m
I
1
y
y
0.00881577
g
9.81 m/s2
I
1
z
z
0.01747731
m
1
1.402 Kg
I
1
x
y
0.00002074
m
2
,
m
3
0.1566 Kg
I
1
x
z
−0.00087669
m
4
0.05 Kg
I
1
y
z
0.00000808
d
1
(6.72, 0.342, −78.9) mm
I
2
x
x
,
I
3
x
x
0.00004223
d
2
(0, −247, 123) mm
I
2
y
y
,
I
3
y
y
0.00004096
d
3
(0, 247, 123) mm
I
2
z
z
,
I
3
z
z
0.00002658
k
d
0.0179 m
I
4
x
x
,
I
4
y
y
0.00000264
μ
γ
0.005 N⋅s/m
I
4
z
z
0.00000264The simulated system assumes the position and speed sensors with sampling period equal toT
s
=
0.1 s and the controller’s sampling rate is equal to τ
s
=
0.01 s. The simulation tool used a fixed-step simulation with period τ
s and the measures pass through a zero-order-holder with period T
s. Furthermore, positioning system measurements are assumed to have Gaussian distribution with accuracy of 3
σ
=
±
0.15 m (this value is equivalent to the specifications of the Novatel OEMStar GPS receiver, which is used in our project) and a speed sensor with accuracy of 3
σ
=
±
0.12
m
/
s (specification of the DJI Guidance system), where σ is the standard deviation. The translational forces T
x
I, T
y
I, and T
z
I are estimated using the tilt-rotor model and the applied control input, and Gaussian noise with 3
σ
=
±
1
N is added to their values. The matrix Q is given by(71)
Q
=
B
T
s
B
T
s
T
.In order to evaluate the performance of the control system in damping the load swing, the NLLS-2L controller is compared with another control strategy that only assumes the load as a disturbance to the system. This controller is named here as nonlinear controller for trajectory tracking with two control loops (NLPT-2L), which makes use of the same controller designed for the inner-loop, while for the outer-loop the linear control laws (49) and (50) are replaced simply by PID controllers for the motion of x and y given by (72)
v
x
=
x
¨
R
e
f
+
K
d
x
e
˙
x
+
K
p
x
e
x
+
K
I
x
∫
e
x
d
t
,
(73)
v
y
=
y
¨
R
e
f
+
K
d
y
e
˙
y
+
K
p
y
e
y
+
K
I
y
∫
e
y
d
t
.By using the NLPT-2L controller and linearizing the system in the neighborhood of the equilibrium point it is possible to verify that γ
1 and γ
2 are stable if there are drag terms multiplying γ
˙
1 and γ
˙
2 in matrix μ. If, otherwise, these terms are null, then the load swings uninterruptedly (γ
1 and γ
2 are marginally stable). Since drag forces are present in mechanical systems at atmospheric conditions, it is fair to state that these terms are nonnull, implying stability of the whole system in a neighborhood of the equilibrium point.In addition, a third control strategy is also compared with the proposed one, which makes use of the nonlinear cascade controller with three loops (NLLS-3L) presented in [22] combined with the LKFUI described in Section 3.4. For the NLLS-3L strategy, the linear controllers of the outermost loop and altitude are synthesized with the same parameters of the NLLS-2L. However, since the attitude and tiltable dynamics are controlled by two control loops in the NLLS-3L strategy presenting high coupling, their gains are obtained similarly with the tuning used in [22].The closed-loop requirements used for the linear mixedH
2
/
H
∞ controllers synthesis with pole placement constraints are presented in Table 2, for the LMI regions (55)–(57), and Table 3, for the weighting matrices of H
2 and H
∞ norms. From these requirements, the mixing H
2
/
H
∞ optimization problem described in Section 3.3 was solved using the Yalmip toolbox [30] (version R20170626) with the Sedumi solver, yielding the controllers gains for each linearized dynamic.Table 2
Closed-loop pole placement requirements for the NLLS-2L controller.
Variable
α
c
r
φ
x
-
1
0
-
3
π
/
36
y
-
1
0
-
3
π
/
36
z
-
1.1
0
-
80
π
/
36
ϕ
-
7
-
7
-
110
π
/
4
θ
-
7
-
7
-
110
π
/
4
ψ
-
3
-
3
-
20
π
/
10
d
˙
x
I
-
1.5
0
-
2
π
/
20
d
˙
y
I
-
1.5
0
-
2
π
/
20Table 3
Weighting matrices ofH
2 and H
∞ norms for the NLLS-2L controller.
Variable
C
z
D
u
z
D
π
z
x
y
d
i
a
g
0.1
1.0
1
/
(
2
)
2
0.5
0.1
0.1
T
[
0.1
0.1
1
]
T
z
ψ
d
i
a
g
(
0.1
1.0
1
/
(
2
)
2
1
/
(
0.2
)
2
1
/
(
0.3
)
2
)
[
0.5
0.1
0.1
0.1
0.1
]
T
[
0.1
0.1
0.1
0.1
1
]
T
ϕ
θ
d
i
a
g
(
1500
1000
1
/
(
0.5
)
2
1
/
(
0.7
)
2
1
/
(
0.3
)
2
)
[
0.05
0.001
0.1
0.1
1.3
]
T
[
0.3
0.01
0.01
1.0
1.0
]
T
d
˙
x
I
d
˙
y
I
d
i
a
g
(
0.1
1.0
)
[
0.5
0.1
]
T
[
0.1
0.1
]
TFor the sake of comparison, the simulations of NLLS-2L, NLPT-2L, and NLLS-3L took place with exactly the same noise in function of time (same source). The external disturbance profiles applied for the three axis are shown in Figure3. Figure 4 shows the position estimation error in the direction x for the ten first seconds of all simulations. It is possible to see that the estimation for NLLS-2L was not much different than the ones for NLPT-2L and NLLS-3L. It can also be seen that LKFUI was able to converge the position estimation, with a priori estimation error bounded to be ±
3.81
c
m with 99.8
% of confidence or, likewise, bounded to be ±
2.54
c
m with 95
% of confidence. Figure 5 shows the disturbance estimation in the z direction for the NLLS-2L simulation, which is estimated with σ
=
0.92N as standard deviation. The position estimation for y and z and the disturbance estimation for d
x and d
y present similar profile of estimation convergence with respect to the estimation of x and d
z, respectively.Figure 3
System disturbances for the simulation.Figure 4
e
x
(
i
∣
k
) in function of time for NLLS-3L, NLLS-2L, and NLPT-2L simulations.Figure 5
Disturbance estimation in function of time for the NLLS-2L simulation.Figure6 shows the results of the simulated trajectories. The aircraft started on the point x
=
0
,
y
=
0 and its initial values for γ
1 and γ
2 were both π
/
10 and -
π
/
10, respectively. Both NLLS-2L and NLLS-3L strategies reduced the load’s swing, but the controller with three loops presented an overall worse trajectory tracking. On the other hand, the NLLS-2L presented a trajectory tracking error for x, y, and z close to the one obtained using NLPT-2L, as it can be observed in Figure 7. This figure also shows that γ
1 and γ
2 presented greatly reduced swing when using the NLLS controllers, with both two and three loops. The time evolution of the remaining generalized coordinates is shown in Figure 8, showing that they converge to operating points without destabilizing throughout the simulation. Figure 9 shows the system’s control inputs.Figure 6
Trajectory tracking of the aircraft for NLPT-2L, NLLS-2L, and NLLS-3L.Figure 7
Tracking error and load swing on NLPT-2L, NLLS-2L, and NLLS-3L.Figure 8
Body and tiltable mechanism angles on NLPT-2L, NLLS-2L, and NLLS-3L.Figure 9
Inputs of the system using NLPT-2L, NLLS-2L, and NLLS-3L.Table4 presents the root mean square error (RMSE) of the aircraft’s trajectory with respect to its reference and the RMSE of the load’s translational velocity with respect to the aircraft for the control strategies. The NLLS-2L controller presented 85
% of the root mean square error on the direction x with respect to the NLPT-2L controller and 68
% in the direction y. The tracking errors for z and ψ were similar on both controllers. The RMSE for the derivatives of the load angles γ
1 and γ
2 presented an improvement of approximately 136
% for γ
1 and 47
% for γ
2 when comparing NLLS-2L with NLPT-2L.Table 4
Root-mean-square-error comparison between NLPT-2L, NLLS-2L, and NLLS-3L.
NLPT-2L
NLLS-2L
NLLS-3L
R
M
S
E
N
L
P
T
-
2
L
/
R
M
S
E
N
L
L
S
-
2
L
R
M
S
E
N
L
L
S
-
3
L
/
R
M
S
E
N
L
L
S
-
2
L
R
M
S
E
x
(
m
)
5
,
58
⋅
1
0
-
2
6
,
54
⋅
1
0
-
2
1
,
09
⋅
1
0
-
1
0.85
1
,
67
R
M
S
E
y
(
m
)
3
,
96
⋅
1
0
-
2
5
,
78
⋅
1
0
-
2
1
,
04
⋅
1
0
-
1
0
,
68
1
,
79
R
M
S
E
z
(
m
)
6
,
37
⋅
1
0
-
2
6
,
35
⋅
1
0
-
2
4
,
40
⋅
1
0
-
2
1
,
00
0
,
69
R
M
S
E
ψ
(
r
a
d
)
2
,
02
⋅
1
0
-
1
2
,
02
⋅
1
0
-
1
1
,
23
⋅
1
0
-
1
0
,
99
0
,
60
R
M
S
E
d
˙
x
I
(
m
/
s
)
3
,
45
⋅
1
0
-
1
1
,
46
⋅
1
0
-
1
2
,
48
⋅
1
0
-
1
2
,
36
1
,
70
R
M
S
E
d
˙
y
I
(
m
/
s
)
2
,
31
⋅
1
0
-
1
1
,
56
⋅
1
0
-
1
2
,
31
⋅
1
0
-
1
1
,
47
1
,
47The NLLS-2L controller presented a trajectory tracking improvement of67
% of the root mean square error on the direction x with respect to the NLLS-3L controller and 79
% in the direction y. The tracking errors for z and ψ were better in NLLS-3L, with an improvement of 31
% and 40
%, respectively. The RMSE for the derivatives of the load angles γ
1 and γ
2 presented an improvement of approximately 70
% for γ
1 and 47
% for γ
2 when comparing NLLS-2L with NLLS-3L.Table5 presents the Integrated Absolute Variation of the Control signal (IAVU) index for all controllers. This index evaluates the control effort and is given by(74)
I
A
V
U
i
=
∫
0
t
f
d
u
i
d
t
d
t
,
i
∈
f
R
,
f
L
,
τ
α
R
,
τ
α
L
.Table 5
IAVU index comparison between NLPT-2L, NLLS-2L, and NLLS-3L.
NLPT-2L
NLLS-2L
NLLS-3L
I
A
V
U
N
L
P
T
-
2
L
/
I
A
V
U
N
L
L
S
-
2
L
I
A
V
U
N
L
L
S
-
3
L
/
I
A
V
U
N
L
L
S
-
2
L
I
A
V
U
f
R
2
,
359
⋅
1
0
4
2
,
349
⋅
1
0
4
1
,
091
⋅
1
0
5
1.004
4.642
I
A
V
U
f
L
2
,
420
⋅
1
0
4
2
,
379
⋅
1
0
4
9
,
450
⋅
1
0
4
1.017
3.972
I
A
V
U
τ
α
R
1
,
401
⋅
1
0
3
1
,
403
⋅
1
0
3
2
,
127
⋅
1
0
3
0.998
1.515
I
A
V
U
τ
α
L
1
,
592
⋅
1
0
3
1
,
601
⋅
1
0
3
2
,
765
⋅
1
0
3
0.994
1.727Table5 indicates that the control efforts in NLPT-2L and NLLS-2L are quite similar. In addition, the two-loop strategies reduced the control effort over 300
% for f
R and f
L when compared with the three-loop strategy. In addition, the control effort for the tilting torques is greatly reduced in the two-loop strategy with respect to the three-loop one.


## 5. Conclusion

This paper provided a solution for the nonlinear control problem of the tilt-rotor UAV with suspended load. The proposed control strategy was able to stabilize the whole system while tracking a desired trajectory in the presence of external disturbances, parametric uncertainties, unmodeled dynamics, and noisy measurements with lower sampling frequency than the controller.The state-space model that was used for the LKFUI algorithm simplifies the dynamics of the tilt-rotor, since the translational generalized coordinates (x, y, and z) are external variables. This way, there is no need to use a nonlinear state estimation algorithm.Simulation results were presented where the LKFUI algorithm was able to estimate the disturbances and use its information along with the aircraft translational forces so as to estimate the system’s position while the position and speed sensors do not provide new measurements. At the moment that the sensors provide a new measurement, a correction step is taken to merge the information from the measurement along with the previous estimation. Consequently, the LKFUI was able to reduce estimation error variability from±
0.15m (positioning system’s measurement error) to approximately ±
0.0381m.Three nonlinear control strategies were compared while performing the task of trajectory tracking: NLPT-2L, NLLS-2L, and NLLS-3L. NLLS strategies greatly reduce the load’s swing with respect to NLPT-2L (up to147
% improvement in the RMSE sense) at the cost of worsening its tracking error for x and y when subject to stronger changes in reference (worsening up to 32
% in the RMSE sense). However, even though NLPT-2L performed better trajectory tracking than NLLS-2L, it is not advised to be used since too much swing of the load, especially for heavy loads, might destabilize the aircraft or also cause the load to collide with the environment.In terms of control effort, NLLS-2L and NLPT-2L were quite similar. On the other hand, these strategies largely reduced the control effort when compared with NLLS-3L. This is a compelling feature that motivates the abandoning of the three-loop strategy in favor of NLLS-2L.An interesting future work would be to reduce the simplifications taken and also study the domain of attraction of the designed controllers, since they were only analyzed in the neighborhood of the equilibrium point. Another important step is to implement the solutions in a real tilt-rotor UAV to perform load transportation. This is expected to be performed in the future by the researchers of the ProVant project, whose research is presently carried out in a partnership between the Brazilian universities UFMG and UFSC.



---

*Source: 1467040-2018-06-27.xml*