Mechanical Design, Modelling and Control of a Novel Aerial
Manipulator
Alexandros Nikou, Georgios C. Gavridis and Kostas J. Kyriakopoulos
Abstract— In this paper a novel aerial manipulation system
is proposed. The mechanical structure of the system, the
number of thrusters and their geometry will be derived from
technical optimization problems. The aforementioned problems
are defined by taking into consideration the desired actuation
forces and torques applied to the end-effector of the system.
The framework of the proposed system is designed in a CAD
Package in order to evaluate the system parameter values.
Following this, the kinematic and dynamic models are developed
and an adaptive backstepping controller is designed aiming to
control the exact position and orientation of the end-effector in
the Cartesian space. Finally, the performance of the system is
demonstrated through a simulation study, where a manipulation
task scenario is investigated.
I. INTRODUCTION
Aerial manipulation is a new scientific field which has
been gaining significant research attention and a wide variety
of structures have been proposed in the last years. These
manipulation systems possess several features which have
lately brought them in the spotlight, with their objective
mainly oriented towards performing effectively complex
manipulating tasks in unstructured and dynamic environ-
ments. Having them include active manipulation as a major
functionality, would vastly broaden the applications of these
systems, as they move from mere passive observation and
sensing to interaction with the environment. Therefore, new
scientifically applicable horizons will be introduced related
to cooperative manipulation, surveillance, industrial inspec-
tions, inspection and maintenance of aerial power lines,
assisting people in rescue operations and constructing in
inaccessible sites by repairing and assembling. Naturally,
both designing and controlling aerial manipulators could be
considered as nontrivial engineering challenges.
The first theoretical and experimental results on aerial
robots interacting with the environment were developed
in [1], [2] using a ducted-fan prototype UAV within the
framework of AIRobots project. The design of a quadrotor
capable of applying force to a wall maintaining flight stability
was performed in [3]. In [4] experimental results with a small
helicopter with grasping capabilities were derived, along with
the stability proofs while grasping. Several grippers that al-
low quadrotors to grasp, pick up and transport payloads were
introduced in [5]. An implementation of indoor gripping
using a low-cost quadrotor has been introduced in [6]. The
Alexandros
Nikou,
Georgios
C.
Gavridis
and
Kostas
J.
Kyriakopoulos
are
with
the
Control
Systems
Lab,
Department
of
Mechanical
Engineering,
National
Technical
University
of
Athens,
9
Heroon
Polytechniou
Street,
Zografou
15780,
Greece.
Email:
{mcp12214,mc08042,kkyria}@mail.ntua.gr
authors in [7] addressed the problem of controlling multiple
quadrotor robots that cooperatively grasp and transport a
payload in three dimensions. Another significant work with
cooperative quadrotors throwing and catching a ball with a
net was performed in [8].
A dexterous holonomic hex-rotor platform equipped with
a six DoF end-effector that can resist any applied wrench was
proposed in [9]. A system for aerial manipulation, composed
of a helicopter platform and a fully actuated seven DoF
redundant robotic arm, has been introduced in [10]. Another
hex-rotor manipulator that consists of three pair of propellers
with a two-link manipulator aiming to trajectory tracking
control was studied in [11].
More recently, significant experiments using commercial
quadrotors equipped with external robotic arms have been
conducted in [12]–[14].
In this work, a completely novel aerial manipulator is
introduced. This could be considered as a small autonomous
aerial robot that interacts with the environment via an end-
effector by applying desired forces and torques in a 6 DoF
task space. The proposed system provides mechanical design
flexibility achieved through technical optimization problems.
The structural geometric distribution is the outcome of the
aforementioned problems with the main goals being oriented
towards low body volume, controllability of the system,
avoidance of possible aerodynamic interactions and effi-
ciency in performing desired manipulation tasks in dynamic
environments. The system is fully integrated as it is not a
commercial aerial robot equipped with an external robotic
arm, as many of the aerial manipulators mentioned above.
The optimal number of thrusters, their positions/orientations
and the optimal position of the end-effector on the body
structure are defined with respect to the modelling design
limitations. Taking all the above into account, the remaining
challenge is is to actually construct this novel aerial robot.
The rest of paper is organized as follows. In Section II a
functional description of the robot and the mechanical design
analysis is discussed. A mathematical model that captures
the proposed system dynamics and govern the behaviour of
the system is derived in Section III. Based on this highly
nonlinear model, an adaptive backstepping control law is
designed in Section IV. In Section V, simulation results are
presented in order to study the performance of the system.
Finally, the main conclusions are discussed in Section VI.
II. MECHANICAL DESIGN
The overall description of the Aerial Manipulator was
based on the idea of designing an aerial robot composed
arXiv:1502.07424v1  [cs.RO]  26 Feb 2015
of a set number n of similar thrusters and an end-effector, in
order to interact with objects in the environment. The exact
geometry of the structure will be the result of the analysis
of this section.
A. Principles of the Problem
Initially, we define the Body-Fixed frame and the End-
Effector frame as FB = {ˆxB, ˆyB, ˆzB}, FE = {ˆxE, ˆyE, ˆzE}.
These frames are attached to the rigid body of the aerial
manipulator as in Fig. 1. The vectors ri, re ∈R3 denote
the position of each thruster and the position of the end-
effector respectively with reference to the Body-Fixed frame.
The thruster orientations are given by the unit vectors ˆFi ∈
R3, i = 1, ..., n, the thrust forces are defined as λi and
the propulsion vectors are given by λi ˆFi. At this point,
it is assumed that the total system is considered to be a
rigid body and, without loss of generality, the End-Effector
frame and the Body frame have the same orientation. Thus,
the actuation force applied to the end-effector is Fact

B =
Fact

E ∈R3, where

B,

E denote the expressions to the
frames FB, FE respectively. The corresponding actuation
torque is obtained via the formula Tact

B = Tact

E +re ×fe
where Tact

E = r × fe is the torque produced by the
end-effector. The terms r, fe denote the displacement vector
(length of the lever arm) and the vector force that tends to
rotate a gripped object from the end-effector.
Fig. 1: Aerial Manipulator Frame Configuration System
B. Forces and Torques
The forces transmitted essentially through the end-effector
are written as
n
X
i=1
(λi ˆFi) + W = Fact

B
(1)
where W ∈R3 is the vector that corresponds to the total
weight of the system. By separating the weight as ws =
W −n w, where w ∈R3 is the weight of each thruster, (1)
is modified as
F λ + n w + ws = Fact

B
(2)
where λ = [λ1 · · · λn]τ ∈Rn and F = [ ˆF1 · · · ˆFn]τ ∈R3×n.
Similarly, the torque from each thruster is Ti = ri ×
(λi ˆFi) = λiS(ri) ˆFi. The well-known skew-symmetric ma-
trix S(·) ∈R3×3 is defined as a × b = S(a) b for the
cross-product × and any vectors a, b ∈R3. The torque due
to the weight is calculated as
TW = rG×W =
n
X
i=1
ri×w+rs×ws =
n
X
i=1
S(ri)w+S(rs)ws
where rG is the centre of gravity of the system and rs
is the centre of gravity of the system, when omitting the
mass of each thruster (mthr). The reaction torque of each
thruster is τi = µ (λ ˆFi) where µ is a coefficient that
represents the relationship between the thrust force and the
reaction torque [15]. Therefore, by combining all torques the
following equation holds
n
X
i=1
n
λi S(ri) ˆFi
o
+ µ
n
X
i=1

λi ˆFi

+
n
X
i=1
S(ri) w + S(rs) ws = Tact

B
(3)
Using the matrices
r =
r1
· · ·
rn
τ ∈R3×n, ¯F = Fact

B −n w −ws (4)
E(r, F) =

S(r1) ˆF1
· · ·
S(rn) ˆFn

∈R3×n
(5)
in (2),(3) we get





F λ = ¯F
E λ = Tact

B −µ ¯F −
n
X
i=1
S(ri) w −S(rs) ws
(6)
By defining the matrix D(r, F) =

F
E(r, F)

∈R6×n, from
the system (6) it is implied that
D(r, F) λ = WR
(7)
where the augmented wrench vector WR ∈R6 is given by
WR =


¯F
Tact

B −µ ¯F −
n
X
i=1
S(ri)w −S(rs)ws


(8)
C. Negative Thrust Forces
It is clear that when solving (7), the vector that corre-
sponds to the thrust force λ can obtain any value in R6.
However, the thrusters are optimally designed to produce
thrust force towards a specific direction, which we set to
correspond to the positive values of λi. In order to alleviate
the problem of negative λi, a conservative solution is adopted
in this analysis, which is based on the idea of introducing
one additional thruster. Thus, (7) is rewritten as
n
X
i=1
λi ti = WR
(9)
where
D(r, F) = [t1 t2 · · · tn] and ti =

ˆFi
S(ri) ˆFi

∈R6 (10)
for all i = 1, ..., n. The vector
ta = −
n
X
i=1
ti =


−
nP
i=1
ˆFi
−
nP
i=1
n
S(ri) ˆFi
o

=

ˆFa
S(ra) ˆFa

(11)
that corresponds to the additional thruster is introduced.
Using (11), the position vector ra and the direction ˆFa
of the new thruster should satisfy the equations
ˆFa =
−
nP
i=1
ˆFi , S( ˆFa) ra = −
nP
i=1
n
S( ˆFi) ri
o
.
If we assume that (7) results in some negative thrust
forces, then the set σN = {k : λk < 0, k = 1, ..., 6} denotes
the indexes for every negative thrust force and σP
=
{1, 2, 3, 4, 5, 6}−σN the corresponding set of positive thrust
forces. Observing that λk < 0 ⇔(−λk) > 0, ∀k ∈σN,
(9) can be separated in
X
i∈σP
λi ti +
X
k∈σN
λk tk = WR
⇔
X
i∈σP
λi ti +
X
k∈σN
(−λk) (−tk) = WR
(12)
Now, from (11) the following can be exported
ta = −
n
X
i=1
ti = −
X
i∈σP
ti −
X
j∈σN
tj
(13)
It is obvious that
−
X
j∈σN
tj = −tk −
X
j∈σN
j̸=k
tj, ∀k ∈σN
(14)
Combining (13), (14) we obtain
−tk = ta +
X
i∈σP
ti +
X
j∈σN
j̸=k
tj, ∀k ∈σN
(15)
By substituting (15) into (12) we have
X
i∈σP
λi ti +
X
k∈σN
(−λk)

ta +
X
i∈σP
ti +
X
j∈σN
j̸=k
tj

= WR
Defining
∆=
X
k∈σN
(−λk) > 0, Ek =
X
j∈σN
j̸=k
(−λj) > 0
(16)
and rearranging the terms we result in
X
i∈σP
λi + ∆


ti +
X
k∈σN
Ek tk + ∆ta = WR
(17)
From (17) the thruster redistribution among all thrusters after
adding the new thruster is provided. It has been proven
that the issue of negative thrust forces can be alleviated
with adding one extra thruster. This equation can be better
analysed in Fig. 2 in which the thrust redistribution algorithm
is depicted. The variables λ′, λ′
i, λ′
k denote the initial thrust
forces and the other variables the thrust force after the
redistribution, plus the additional thrust force λa. Thus, the
six thrust forces (not necessary all positive) are equivalent
to seven thrust forces, all positive with redistributed thrust
forces as in (17). By using the additional thruster, (4), (8)
are reformed into
¯F = Fact

B −(n + 1) w −ws
(18)
WR =


¯F
Tact

B −µ ¯F −
n+1
X
i=1
S(ri) w −S(rs) ws


(19)
Fig. 2: Thrust Force Equivalence
D. Aerodynamic Interaction
At this point, the aerodynamic interaction between the
operation thrusters is investigated. The aerodynamic effects
produced by each thruster, are based on experiments that took
place in Control Systems Lab NTUA on a 8 × 4.7SF APC
propeller accompanied with the Neu Motor NEU 1902/2Y -
2035 motor, which produces at 17550 rpm, a λmax = 28N
thrust force. The surface, corresponding to every thruster,
that approximates these effects is described by a third order
equation in (SI). By expressing this equation in the thruster
frame FTi = {ˆx′
i, ˆy′
i, ˆz′
i}, i = 1, ..., 6, a we get
−0.06 ≤x′
i ≤0.91
(20)
(y′
i)2 + (z′
i)2 ≤

−1.1(x′
i)3 + 1.56(x′
i)2 −0.3(x′
i) + 0.11
2
Hence, the aerodynamic effects of the air flow throughout the
rotor are extended from x = −0.06m to x = 0.91m. The x
axis shows the length of the aerodynamic effect of the exit
flow. In order to understand this, one should consider the
rotor/blade to be positioned at x = 0. On the other hand,
y axis shows the distance of the effect measured from the
rotation axis of the blade, where at position (x = 0, y =
0.102m) (SI) is the blade radius in approximation (because
of the existence of an offset).
An arbitrary point p = [x y z]τ expressed in FB and the
corresponding p′
i = [x′
i y′
i z′
i]τ expressed in FTi, can be
linked together by the equation p = {TR
FTi
FB (ri, ˆFi)} p′
i,
where TR
FTi
FB (ri, ˆFi) is the appropriate homogeneous frame
transformation corresponding to the translation and orien-
tation vectors (ri, ˆFi). By combining the previous coordi-
nates transformation equation with the constraints (20), a
set of constraints that can be described in matrix form as
G(ri, ˆFi, p) ≤0, is produced. The distance between two
such volumes i, j can be defined and evaluated via the
optimization problem (P1) of Table I.
E. Design Problem
Given a particular structure defined by the matrices
(r, F), for a set of required actuation forces and torques
(Fact

E, Tact

E) it is necessary to find the associate levels of
the thrust forces λi. Since WR ∈R6, in order for (7) to have
a solution for λ ∈Rn, the conditions {rank(D) = 6, n ≥6}
are required. The rank condition is adequate from a strict
mathematical perspective but from a practical point of view,
as (7) leads to the thrust forces values λ ∈Rn, the sought
solutions should not be very sensitive to small deviations.
This is partially achieved by using the condition number
κ(D) = σmax(D)/σmin(D) where σ(D) =
p
eig(DτD)
are the singular values of the matrix D, eig(·) denotes
the eigenvalues of a matrix and σmax(D), σmin(D) are the
maximum and minimum singular values of the matrix D
respectively. Thus, a low condition number κ(D) ≥1 is
required [16]. Although the condition number is bounded
to take feasible values (not equal to zero/infinity) when
σ(D) →0, the matrix D(r, F) might be ill-conditioned i.e.
det(D(r, F)) →0. Thus, σ(D) ≥ϵ1 > 0. Furthermore, to
avoid the fan interaction an other constraint is introduced
as dij(ri, ˆFi, rj, ˆFj) ≥ϵ2 > 0, ∀i, j = 1, 2, . . . , n, α. Note
that, similarly to (P1), the position re should be introduced
to the design problem as the intersection avoidance between
a sphere (with radius Re) that encloses the end-effector,
and the thrusters. This sphere, when expressed in the End-
Effector frame FE, is given by (x′
e)2 + (y′
e)2 + (z′
e)2 ≤R2
e.
Therefore, the constraint associated with the end-effector is
dei(re, ri) ≥Re > 0, ∀i = 1, 2, . . . , n, α. An optimization
is also required to minimize the volume of the system, by
using the norm J(r) = ∥r∥2. Taking all the above into
consideration, the design problem is essentially recast to the
optimization problem (P2) from Table I. The optimization
parameters are chosen as K = 5, ϵ1 = 10−3, ϵ2 =
10−2 m, Re = 10−2 m.
(P1)
dij(ri, ˆFi, rj, ˆFj) = min
pi,pj∥pi −pj∥
s.t.
G(ri, ˆFi, pi) ≤0
G(ri, ˆFi, pj) ≤0
(P2)
min
r,re, ˆ
F
J(r)
s.t.
σ(D) ≥ϵ1
dij ≥ϵ2, ∀i, j = 1, 2, . . . , n, α
dei ≥Re, ∀i = 1, 2, . . . , n, α
ˆFa = −
n
X
i=1
ˆFi
S( ˆFa) ra = −
n
X
i=1
S( ˆFi) ri
1 ≤κ(D) ≤K
TABLE I
OPTIMIZATION PROBLEMS
F. Solving the Optimization Problem
It should be noted that when solving the optimization
problem (P2), each time the inner problem (P1) should be
solved. There are 45 decision variables of the optimization
problem, which correspond to the seven position vectors
(ri) of the thrusters, the position vector (re) of the end-
effector and the direction vectors ( ˆFi) of the seven thrusters.
This issue, entails the necessity of solving 28 optimization
problems for each evaluation attempt of the outer problem
(P2).
The inner problem, that refers to the avoidance of the fan
interaction, is smooth but in terms of the outer problem
(P2) is nonsmooth and nonlinear. The objective function
and the constraints of the problem (P1) are continuous
and this problem, according to the inputs, has one and
only one global minimum. Using the appropriate rotation
and transformation matrices, the (P1) was solved by the
active-set strategy [17],[18]. On the other hand, the design
problem (P2) has nonsmooth, discontinuous and nonlinear
inequality constraints, but smooth objective function. Con-
sequently, a non-gradient-based methodology that searches
disjoint feasible regions, is utilized. For the pre-search of the
design space, a Latin Hypercube (LHS) [19] was chosen, in
order to ensure that the points are distributed throughout the
search space. The Latin Hypercube sampling is known to
provide better coverage than the simple random sampling
[20]. Following this, a Generalized Pattern Search (GPS)
direct search algorithm [21],[22] was used.
The thrust force (λ) and the momentum (Q) can be
calculated from [23],[24] as
(
Q = πρ CQ R5 Ω2
λ = πρ Cλ R4 Ω2
⇔
Q = CQ
Cλ
R λ
(21)
where the term CQ
Cλ R corresponds to the coefficient µ, R
is the radius of the rotor and ρ, Ωdenote the air density
and the rotational speed of the rotor respectively. Applying
a combination of the Blade Element Theory [24] and the
Momentum Theory [15], using the modified versions pro-
posed in [23] and invoking the experimental results extracted
by our lab on the APC propeller, it was calculated that
Cλ = 0.008, CQ = 0.0095, µ = 0.1473, R = 0.124m. By
solving the optimization problems, with the results depicted
in Table II, the matrix D(r, F) is full rank and using (7),
the thrust forces can be calculated as λ = D−1 WR. All the
constraints were satisfied and a low volume body structure
with condition number κ(D) = 3.36 resulted. The wrench
vector WR can be determined by substituting the desired
actuation forces/torques (Fact

E, Tact

E) in (19). The maxi-
mum thrust force and torque which can be applied from the
system are λmax = 28N and 3Nm respectively. The values
of the components, proposed for the Aerial Manipulator,
are the following: the motor and the propeller (0.12kg),
the frame (0.66kg), the battery (0.25kg) and the electronic
components (0.15kg). The total mass of the proposed system
is m = 1.90kg. Ultimately, the production of a carefully
studied framework (Fig. 3) was achieved by using the 3D
CAD Package SolidWorks. Using this Package, the system
parameter values of the Table II have been evaluated.
Fig. 3: Aerial Manipulator 3D Caption of the Framework
Param.
Description
Value
Units
m
Total Mass
1.90
kg
mthr
Thruster Mass
0.12
kg
IG
Moment of
Inertia Tensor

0.3488
0.0683 −0.0457
0.0683
0.1588
0.0144
−0.0457 0.0144
0.4081

kg m2
rG
Centre of
Gravity
Position
[0.0737 0.0083 −0.0781]τ
m
re
End-Effector
Position
[−0.23 0.015 0.23]τ
m
rs
Centre of Grav.
from (3)
[0.1267 −0.0052 −0.1900]τ
m
J(r)
Total Structure
Volume
1.80018
m3
g
Gravitational
Acceleration
9.81
m/s2
ri
r1 = [0.43 −0.15 −0.44]τ
m
r2 = [0.08 −0.22 −0.14]τ
r3 = [0.1 −0.9 −0.2]τ
Thruster
Positions
r4 = [−0.34 0.25 0.006]τ
r5 = [0.184 0.359 −0.254]τ
r6 = [−0.22 −0.44 −0.04]τ
r7 = [0.51 0.79 −0.06]τ
ˆFi
ˆF1 = [0.08 0.39 0.92]τ
ˆF2 = [−0.33 −0.90 0.29]τ
ˆF3 = [0.13 −0.87 −0.48]τ
Thruster
Orientations
ˆF4 = [0.56 0.08 0.82]τ
ˆF5 = [0.83 0.11 −0.55]τ
ˆF6 = [−0.66 0.57 −0.49]τ
ˆF7 = [−0.59 0.62 −0.51]τ
TABLE II
AERIAL MANIPULATOR PARAMETERS
III. MATHEMATICAL MODEL OF THE AERIAL
MANIPULATOR
In this section, the kinematic and dynamic equations of
motion in case there are no interaction forces and torques
from the environment applied to the end-effector are pre-
sented.
A. Kinematic Model
Fig. 1 shows the reference frames defined to derive the
kinematic and dynamic model of the proposed system. The
Earth-Fixed inertial frame is defined as FA = {ˆxA, ˆyA, ˆzA}
and it should be noted that the Body-Fixed frame’s origin
does not coincide with the centre of gravity G. The position
of FB relative to FA can be represented by p = [x y z]τ ∈
R3 and the corresponding orientation by the rotation angles
Θ = [φ θ ψ]τ
∈R3. The translational and rotational
kinematic equations of the moving rigid body are given (see
[25]) in matrix form by
˙ξ =
 ˙p
˙Θ

=

Jt(Θ)
O(3×3)
O(3×3)
Jr(Θ)
 
v
ω

(22)
where O(3×3) is the 3 × 3 zero matrix, v = [vx vy vz]τ ∈
R3, ω = [ωx ωy ωz]τ ∈R3 denote the translational velocity
and the angular velocity of FB relative to FA respectively,
both expressed in the Body-Fixed frame. The transformation
matrices Jt(Θ), Jr(Θ) ∈R3×3 are given by
Jt(Θ) =


cθcψ
sφsθcψ −sψcφ
sθcφcψ + sφsψ
sψcθ
sφsθsψ + cφcψ
sθsψcφ −sφcψ
−sθ
sφcθ
cφcθ

(23)
Jr(Θ) =


1
sφtθ
cφtθ
0
cφ
−sφ
0
sφ/cθ
cφ/cθ


(24)
The position of the end-effector with respect to FA is pe =
[xe ye ze]τ = p+Jt(Θ) re ∈R3. Its derivative is obtained as
˙pe = Jt(Θ) v −Jt(Θ) S(re) ω, using the formula ˙Jt(Θ) =
Jt(Θ)S(ω) from [26]. The Body-Fixed and the End-Effector
frame have the same orientation with reference to FA, as
mentioned in Section II, hence Θe = Θ. By combining the
last results the following kinematic equation holds
˙ξe =
 ˙pe
˙Θ

=
 Jt(Θ)
−Jt(Θ) S(re)
O(3×3)
Jr(Θ)

|
{z
}
J(ξe)

v
ω

(25)
where Θe is the orientation of FE relative to FA. The
Jacobian matrix of the system J(ξe) ∈R6×6 relates in
a straightforward way the linear velocity ˙pe and the rate
of change in the rotational angles ˙Θ of the end-effector
expressed in FA, with the Body-Fixed velocities v, ω.
B. Dynamic Model
The dynamic equations can be conveniently written with
respect to the Body-Fixed frame by using the Newton-Euler
formalism (the main concept is discussed extensively in [27],
[28]), as
M
 ˙v
˙ω

+ C(ν)
v
ω

=
F

B
T

B

(26)
where
M =

mI3
−mS(rG)
mS(rG)
IB

, M > 0,
˙M = 0
(27)
is the inertia matrix,
C =

mS(ω)
−mS(ω)S(rG)
mS(rG)S(ω)
−S(IBω)

, C = −Cτ
(28)
is the Coriolis-centripetal matrix, I3 is the 3 × 3 identity
matrix, m is the total mass of the system, IB is the inertia
tensor expressed in FB and ν = [vτ ωτ]τ ∈R6 is the
vector of the Body-Fixed velocities. The inertia tensor can be
written as IB = IG −mS(rG)S(rG) where IG is the inertia
tensor relative to the body’s centre of gravity. The vectors
F

B, T

B ∈R3 describe the forces and torques acting on
the system expressed in the Body-Fixed frame and can be
derived as
F

B
T

B

=

F λ −m g Jτ
t (Θ) e3
E λ + µ F λ −m g S(rG) Jτ
t (Θ) e3

=
F
E

λ
| {z }
propulsion
forces/torques
+
O(3×6)
µ F

λ
|
{z
}
reaction
torques
−m g

I3
S(rG)

Jτ
t (Θ)e3
|
{z
}
gravitational
forces/torques
(29)
where e3 = [0 0 1]τ. Combining (26), (29) and solving with
respect to [˙vτ ˙ωτ]τ we get
 ˙v
˙ω

= −M −1C(ν)
v
ω

+ M −1

F
E + µ F

λ
−m g M −1

I3
S(rG)

Jτ
t (Θ) e3
(30)
⇔˙ν = H(ν) + G(ξe)
|
{z
}
B(ξe,ν)
+N λ
(31)
where the matrices are defined as
H(ν) = −M −1 C(ν) ν, N = M −1

F
E + µ F

> 0 (32)
G(ξe) = −m g M −1

I3
S(rG)

Jτ
t (Θ)e3
(33)
B(ξe, ν) = H(ν) + G(ξe), B : R6 × R6 →R6
(34)
IV. NONLINEAR CONTROL OF THE AERIAL
MANIPULATOR
A manipulation task is usually given in terms of the
desired position and orientation of the end-effector. The
objective of this section is to design a controller for the aerial
manipulator ensuring that the position pe(t) and the orien-
tation Θ(t) of the end-effector track the desired Cartesian
trajectory ξdes(t) =

pτ
e,des(t) Θτ
des(t)
τ ∈R6 asymptotically
while all the closed loop signals remain bounded for all
t ≥0. Firstly, by using formulas (25), (31) the aerial
manipulator model, including the kinematics and dynamics,
can be written as
(S) :
( ˙ξe = J(ξe) ν
˙ν = B(ξe, ν) + N θ⋆
λ λ + d(ξe, ν, t)
(35)
where d : R6 × R6 × R+ →R6 represents the unmodelled
nonlinear dynamics and the environmental disturbances. The
unknown matrix θ⋆
λ
= diag{θ⋆
1, . . . , θ⋆
6} ∈R6×6 with
θ⋆
i
∈[θmin, θmax] = [0.1, 1], is introduced to model the
control actuation failures and the modeling errors among
the thrusters of the system, e.g. if θ⋆
i = 0.8 then the i−th
actuator has 20 % controller effectiveness reduction. The
control inputs of the system are the six independent thrust
forces λi(t), i = 1, . . . , 6 as mentioned in Section II. The
matrix N is full rank with low condition number which
constitutes a vital result of the control oriented optimization
from Section II.
The system (35) is highly nonlinear, cascaded and fully
actuated in the well-known strict feedback form, with vec-
tor relative degree 2. For such systems, the backstepping
controller design has proven to be successful [29], [30].
Due to the fact that the system is in the presence of the
uncertainties θ⋆
λ and the disturbances d(ξe, ν, t), a robust
adaptive controller will be designed in order to tackle them.
The aim is to study if the proposed system with the resulting
geometry from the optimization problems (P1), (P2), the
system specifications from Table II and the aforementioned
uncertainties/disturbances from (35), is capable to perform
specific trajectory tasks efficiently. In order to design the
controller of the system (35), the following assumptions are
required:
Assumption 1: The states of the system ξe, ν are available for
measurement ∀t ≥0 for the following control development.
Assumption 2: The desired trajectories ξdes are known and
bounded functions of time (ξdes ∈L∞) with known and
bounded derivatives ( ˙ξdes, ¨ξdes ∈L∞). Assumption 3: The
disturbance d(ξe, ν, t) = [d1(ξe, ν, t) · · · d6(ξe, ν, t)]τ is
unknown but bounded with |di(ξe, ν, t)| ≤∆i where ∆i
are unknown positive constants for all i = 1, . . . , 6 and
t ≥0. Assumption 4: It is assumed for all t ≥0 that
−π
2 < θ(t) < π
2 . This ensures that the Jacobian matrix is
nonsingular since det(J(ξe)) = 1/cθ. This assumption is
likewise utilized in [26], [31].
• Step 1: To begin with the backstepping controller design,
the position-orientation error of the end-effector is defined as
z1 = ξe −ξdes ∈R6. By differentiating it and using (25) we
get
˙z1 = J(ξe) ν −˙ξdes
(36)
We view ν as a control variable and we define a virtual
control law νdes ∈R6 for (36). The error signal representing
the difference between the virtual and the actual controls
is defined as z2 = ν −νdes ∈R6. Thus, in terms of
the new state variable, (36) can be rewritten as ˙z1
=
J(ξe) z2 + J(ξe) νdes −˙ξdes. Consider now the positive
definite and radially unbounded quadratic Lyapunov function
V1(z1) = 1
2∥z1∥2 = 1
2zτ
1z1. By differentiating it with respect
to time yields
˙V1 = zτ
1 ˙z1 = zτ
1
n
J(ξe) νdes −˙ξdes
o
+ zτ
1J(ξe)z2
(37)
The stabilization of z1 can be obtained by designing an
appropriate virtual control law
νdes = J−1(ξe)
n
˙ξdes −K1z1
o
(38)
where the matrix K1 ∈R6×6, K1 = Kτ
1 > 0 represents the
first controller gain to be designed. Hence, the time derivative
of V1 becomes ˙V1 = −zτ
1K1z1+zτ
1J(ξe)z2. The first term of
on the right-hand of this equation is negative and the second
term will be canceled in the next step.
• Step 2: For the second step we define the matrices of the
parameter estimation errors as ˜∆= [ ˜∆1 · · · ˜∆6]τ = [( ˆ∆1 −
∆1) · · · ( ˆ∆6−∆6)]τ and ˜θλ = diag{(ˆθ1−θ⋆
1), . . . , (ˆθ6−θ⋆
6)}
where ˆ∆i, ˆθi are the estimations of the unknown parameters
∆i, θ⋆
i respectively. The time derivative of the error z2 is
˙z2 = B(ξe, ν) + N θ⋆
λ λ + d(ξe, ν, t) −˙νdes. The Lyapunov
function candidate in this step is chosen as
V2(z1, z2, ˜∆, ˜θλ) = V1+1
2zτ
2z2+1
2
˜∆τΓ−1
∆˜∆+1
2tr(˜θτ
λΓ−1
θ ˜θλ)
where Γθ = Γτ
θ > 0, Γ∆= Γτ
∆> 0 are diagonal adaptation
gain matrices and tr(·) denotes the matrix trace. The time
derivative of V2(z1, z2, ˜∆, ˜θλ) is obtained as
˙V2 = −zτ
1K1z1 + zτ
2 {Jτ(ξe)z1 + B(ξe, ν) + N θ⋆
λ λ −˙νdes}
+ zτ
2d(ξe, ν, t) + ˜∆τ Γ−1
∆
˙ˆ∆+ tr(˜θτ
λΓ−1
θ
˙ˆθλ)
(39)
Using the −zτ
1K1z1 ≤−λmin(K1)∥z1∥2, zτ
2 d(ξe, ν, t) ≤
zτ
2sgn(z2)∆
and
adding
and
subtracting
the
terms
zτ
2N ˆθλλ,
zτ
2sgn(z2) ˆ∆in (39) the following inequality
holds
˙V2 ≤−λmin(K1)∥z1∥2 + zτ
2

Jτ(ξe)z1 + B(ξe, ν) −˙νdes
+ sgn(z2) ˆ∆+ N ˆθλλ
	
−zτ
2N ˜θλλ
−zτ
2sgn(z2) ˜∆+ ˜∆τ Γ−1
∆
˙ˆ∆+ tr(˜θτ
λΓ−1
θ
˙ˆθλ)
(40)
where λmin(K1) denotes the minimum eigenvalue of matrix
K1, sgn(z2) = diag{sgn(z2,1), ..., sgn(z2,6)} and sgn(·)
denotes the sign function. Rearranging the terms and using
the property aτb = tr(b aτ), ∀a, b ∈Rn we get
˙V2 ≤−λmin(K1)∥z1∥2 + zτ
2

Jτ(ξe)z1 + B(ξe, ν) −˙νdes
+ sgn(z2) ˆ∆+ N ˆθλ λ
	
+ ˜∆τ
Γ−1
∆
˙ˆ∆−sgn(z2)z2
	
+ tr{˜θτ
λ(Γ−1
θ
˙ˆθλ −N τz2 λτ)}
(41)
Given the form of ˙V2 from (41) the adaptive control law
and the corresponding parameter estimation update laws for
the nonlinear system (35) to be designed, are
λ(ξe, ν, ˆ∆, ˆθλ) = (ˆθλ)−1N −1
˙νdes −B(ξe, ν) −Jτ(ξe)z1
−sgn(z2) ˆ∆−K2z2}
(42)
˙ˆ∆= Γ∆{sgn(z2)z2 −σ ˆ∆}
(43)
˙ˆθλ = Γθ Proj(ˆθλ, N τz2 λτ)
(44)
where K2 = Kτ
2 > 0 is the second controller gain matrix,
σ is a strictly positive gain (σ-modification rule [32]) and
the projection operator Proj(·, ·) is the same as the one in
[33] with the parameter δ to be designed. By substituting
(42), (43), (44) into (41) and using the property ˜∆τ ˆ∆=
1
2∥˜∆∥2 + 1
2∥ˆ∆∥2 −1
2∥∆∥2 the following inequality holds
˙V2 ≤−λmin(K1)∥z1∥2 −λmin(K2)∥z2∥2+
tr
n
˜θτ
λ
h
Proj(ˆθλ, y) −y
io
|
{z
}
≤0, y = N τz2 λτ
−σ
2 ∥˜∆∥2 −σ
2 ∥ˆ∆∥2 + σ
2 ∥∆∥2
|
{z
}
≤−σ
2 ∥˜∆∥2+ σ
2 ∥∆∥2
The projection operator invoked from [33] contributes to the
negative semi-negativeness of the Lyapunov function since
by definition tr
n
˜θτ
λ
h
Proj(ˆθλ, y) −y
io
≤0, ∀y. Moreover,
it guarantees that if ˆθi(0) ∈[θmin, θmax] is chosen, then
ˆθi(t) ∈[θmin −δ, θmax +δ], ∀i = 1, ..., 6, ∀t ≥0 for suitable
δ > 0. The last result protects the term (ˆθλ)−1 in (42) from
singularity. By defining ¯w = σ
2 ∥∆∥2 > 0 we result in
˙V2 ≤−λmin(K1)∥z1∥2 −λmin(K2)∥z2∥2 −σ
2 ∥˜∆∥2 + ¯w
from which it follows that both errors z1, z2 and the param-
eter estimation ˜∆are uniformly ultimately bounded with re-
spect to the sets Ω1 =
n
z1 ∈R6 : ∥z1∥≤
p
¯w/λmin(K1)
o
,
Ω2
=
n
z2 ∈R6 : ∥z2∥≤
p
¯w/λmin(K2)
o
and Ω∆
=
n
˜∆∈R6 : ∥˜∆∥≤
p
2 ¯w/σ
o
.
Invoking
that
z1, z2
are
bounded and ξdes, νdes ∈L∞then ξe, ν ∈L∞. Since
˜∆, ∆, ˆθλ, θ⋆
λ, ˙νdes are bounded then ˆ∆, ˜θλ, λ ∈L∞. Based
on the above, it is proven that all closed loop signals remain
bounded.
One important issue associated with the controller de-
sign is the analytical form of the time derivative of
νdes,
which
can
be
obtained
from
(38)
as
˙νdes
=
J−1(ξe)
n
¨ξdes −˙J(ξe) νdes −K1 ˙z1
o
, and the time deriva-
tive of J(ξe), which can be calculated by using the ˙Jr(Θ) =
∂Jr
∂φ
˙φ + ∂Jr
∂θ
˙θ, ˙J(ξe) =
 ˙Jt(Θ)
−˙Jt(Θ)S(re)
O(3×3)
˙Jr(Θ)

.
0
5
10
15
20
25
30
−1
−0.5
0
0.5
1
1.5
2
P osition Errors
t [sec]
ex(t), ey(t), ez(t)[m]
 
 
xe −xdes
ye −ydes
ze −zdes
0
5
10
15
20
25
30
−2
−1
0
1
2
Orientation Errors
t [sec]
eφ(t), eθ(t), eψ(t)[rad]
 
 
φ −φdes
θ −θdes
ψ −ψdes
Fig. 4: Position and Orientation Errors
0
5
10
15
20
25
30
−5
0
5
10
15
20
Required T hrust F orces
t [sec]
λ(t) [Newton]
 
 
λ1(t)
λ2(t)
λ3(t)
λ4(t)
λ5(t)
λ6(t)
0
5
10
15
20
25
30
0
5
10
15
20
25
Required T hrust F orces After Redistribution
t [sec]
λ(t)[Newton]
 
 
λ1(t)
λ2(t)
λ3(t)
λ4(t)
λ5(t)
λ6(t)
λa(t)
λmax(t)
Fig. 5: Required Thrust Forces Using the Algorithm (17)
V. SIMULATION RESULTS
In this section, the results of a numerical simulation sce-
nario are presented in order to demonstrate the performance
of the proposed system. The dynamic model in (35) is uti-
lized with system parameters which are depicted in Table II.
15 % controller effectiveness reduction is chosen with θ⋆=
0.85 diag{1, 1, 1, 1, 1, 1}. The end-effector is forced to track
the trajectory pe,des(t) = [cos(0.5t) sin(0.5t) 1.5 + 0.3t]τ
with regulated orientation at Θdes =
 π
3
π
6 −π
4
τ with
reference to the Earth-Fixed frame. The initial conditions
of the system are set to pe(0) = re, p(0) = Θ(0) =
v(0) = ω(0) = 0(3×1). The disturbance is set as d(t) =
[0.5 0.4 sin(2t) 0.4 cos(t) 0.5 0.5 cos(0.8t) 0.6 sin(t)]τ rep-
resenting the unmodelled forces/torques and external distur-
bances. The initial values of the parameters estimations are
set to ˆθ1(0) = · · · = ˆθ6(0) = 0.7 and ˆ∆1(0) = · · · =
ˆ∆6(0) = 0. The controller gains are chosen as K1 =
diag{1, 1, 1, 0.3, 0.3, 0.3}, K2 = 8 diag{1, 1, 1, 1, 1, 1}. The
adaptation gains are selected as σ
=
1.5,
Γ∆
=
13 diag{1, 1, 1, 1, 1, 1}, Γθ = 0.1 diag{1, 1, 1, 1, 1, 1}. The
parameter of the projection operator is set to δ = 0.05. Fig. 4
shows the position and orientation tracking errors. The thrust
forces are provided in Fig. 5. This paper is accompanied by a
video demonstrating the simulation procedure of this Section.
Due to space limitations, a video with an additional Scenario
in better quality (HD) can be found at
https://www.youtube.com/watch?v=DXnzu6XOrXs
VI. CONCLUSIONS
Aerial robots physically interacting with the environment
could be very useful for many applications. In this paper,
we have presented the mechanical design of a novel aerial
manipulator which was the result of technical optimization
problems. A mathematical model for the kinematics and dy-
namics was derived in order to design an adaptive nonlinear
controller to study the system while performing manipulation
tasks. The simulation results illustrate the effectiveness of
the proposed system and the controller to achieve tracking
irrespectively of actuator failures, unmodelled dynamics and
external disturbances. Future work mainly involves the con-
struction of the aerial robot and the conduction of experimen-
tal trials for the proposed framework with the actual system,
in order to verify the theoretical results of this paper.
REFERENCES
[1] L. Marconi, R. Naldi, and L. Gentili, “Modelling and Control of a
Flying Robot Interacting with the Environment,” Automatica, vol. 47,
no. 12, pp. 2571–2583, 2011.
[2] A. Keemink, M. Fumagalli, S. Stramigioli, and R. Carloni, “Mechani-
cal Design of a Manipulation System for Unmanned Aerial Vehicles,”
IEEE International Conference on Robotics and Automation (ICRA),
pp. 3147–3152, 2012.
[3] A. Albers, S. Trautmann, T. Howard, T. Nguyen, M. Frietsch, and
C. Sauter, “Semi-Autonomous Flying Robot for Physical Interaction
with Environment,” IEEE Conference on Robotics Automation and
Mechatronics (RAM), pp. 441–446, 2010.
[4] P. Pounds, D. Bersak, and A. Dollar, “Grasping from the Air: Hov-
ering Capture and Load Stability,” IEEE International Conference on
Robotics and Automation (ICRA), pp. 2491–2498, 2011.
[5] D. Mellinger, Q. Lindsey, M. Shomin, and V. Kumar, “Design, Mod-
eling, Estimation and Control for Aerial Grasping and Manipulation,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
(IROS), pp. 2668–2673, 2011.
[6] V. Ghadiok, J. Goldin, and W. Ren, “Autonomous Indoor Aerial
Gripping Using A Quadrotor,” IEEE/RSJ International Conference on
Intelligent Robots and Systems (IROS), pp. 4645–4651, 2011.
[7] D. Mellinger, M. Shomin, N. M. Nathan, and V. Kumar, “Coopera-
tive Grasping and Transport Using Multiple Quadrotors,” Distributed
Autonomous Robotic Systems, vol. 83, pp. 545–558, 2013.
[8] R. Ritz, W. Muller, M. Hehn, and R. D’Andrea, “Cooperative
Quadrocopter Ball Throwing and Catching,” IEEE/RSJ International
Conference on Intelligent Robots and Systems (IROS), pp. 4972–4978,
2012.
[9] G. Jiang and R. Voyles, “Hexrotor UAV Platform Enabling Dextrous
Interaction with Structures-Flight Test,” 2013 IEEE International
Symposium on Safety, Security, and Rescue Robotics (SSRR), pp. 1–6,
2013.
[10] F. Huber, K. Kondak, K. Krieger, D. Sommer, M. Schwarzbach,
M. Laiacker, I. Kossyk, S. Parusel, S. Haddadin, and A. Albu-Schaffer,
“First Analysis and Experiments in Aerial Manipulation Using Fully
Actuated Redundant Robot Arm,” (IROS), 2013.
[11] M. Kobilarov, “Nonlinear Trajectory Control of Multi-Body Aerial
Manipulators,” Journal of Intelligent and Robotic Systems, vol. 73,
no. 1-4, pp. 679–692, 2014.
[12] S. Kim, S. Choi, and H. Kim, “Aerial Manipulation Using a Quadrotor
with a Two DOF Robotic Arm,” IEEE/RSJ International Conference
on Intelligent Robots and Systems (IROS), pp. 4990–4995, 2013.
[13] M. Orsag, C. Korpela, and P. Oh, “Modeling and Control of MM-
UAV: Mobile Manipulating Unmanned Aerial Vehicle,” Journal of
Intelligent and Robotic Systems, vol. 69, no. 1-4, pp. 227–240, 2013.
[14] A. Jimenez-Cano, J. Martin, G. Heredia, A. Ollero, and R. Cano,
“Control of an Aerial Robot with Multi-Link Arm for Assembly
Tasks,” (ICRA), pp. 4916–4921, 2013.
[15] G. Padfield, “Helicopter Flight Dynamics, The Theory and Application
of Flying Qualities and Simulation Modelling”.
AIAA education
series, American Institute of Aeronautics and Astronautics, 1996.
[16] L. N. Trefethen and D. Bau, “Numerical Linear Algebra”. SIAM,
1997.
[17] P. Gill, W. Murray, and M. Wright, “Numerical Linear Algebra and
Optimization”. Addison-Wesley Publishing Company, 1991.
[18] K. Murty, “Linear Complementarity, Linear and Nonlinear Program-
ming”. Sigma Series in Applied Mathematics, Berlin: Heldermann
Verlag, 1988.
[19] M. Stein, “Large Sample Properties of Simulations Using Latin
Hypercube Sampling,” Technometrics, vol. 29, no. 2, pp. 143–151,
1987.
[20] M. McKay, R. Beckman, and W. Conover, “A Comparison of Three
Methods for Selecting Values of Input Variables in the Analysis of
Output from a Computer Code,” Technometrics, vol. 42, no. 1, pp. 55–
61, 2000.
[21] C. Audet and J. Dennis, “Analysis of Generalized Pattern Searches,”
SIAM Journal on Optimization, vol. 13, no. 3, pp. 889–903, 2003.
[22] C. Audet and J. Dennis, “A Pattern Search Filter Method for Nonlinear
Programming without Derivatives,” SIAM Journal on Optimization,
vol. 14, no. 4, pp. 980–1010, 2004.
[23] V. Martinez, “Modelling of the Flight Dynamics of a Quadrotor
Helicopter,” Master Thesis, Cranfield University, 2007.
[24] R. Prouty, “Helicopter Performance, Stability and Control”.
R.E.
Krieger Publishing Company, 1995.
[25] G. Antonelli, “Underwater Robots”.
Springer Tracts in Advanced
Robotics, Springer International Publishing, 2013.
[26] T. Madani and A. Benallegue, “Backstepping Control for a Quadrotor
Helicopter,” IEEE/RSJ International Conference on Intelligent Robots
and Systems (IROS), pp. 3255–3260, 2006.
[27] T. Fossen, “Guidance and Control of Ocean Vehicles”. John Wiley
and Sons, 1994.
[28] D. Bernstein, “Geometry, Kinematics, Statics and Dynamics”. Uni-
versity of Michigan, 2013.
[29] M. Krsti´c, I. Kanellakopoulos, and P. Kokotovi´c, “Nonlinear and
Adaptive Control Design”. Wiley, 1995.
[30] J. Zhou and C. Wen, “Adaptive Backstepping Control of Un-
certain Systems: Nonsmooth Nonlinearities, Interactions Or Time-
Variations”. Springer, 2008.
[31] M. Huang, B. Xian, C. Diao, K. Yang, and Y. Feng, “Adaptive Track-
ing Control of Underactuated Quadrotor Unmanned Aerial Vehicles
via Backstepping,” American Control Conference (ACC), 2010.
[32] E. Lavretsky and K. Wise, “Robust and Adaptive Control: with
Aerospace Applications”. Springer, 2012.
[33] H. Khalil, “Adaptive Output Feedback Control of Nonlinear Systems
Represented by Input-Output Models,” IEEE Transactions on Auto-
matic Control, vol. 41, pp. 177–188, Feb 1996.