Model Identication and H∞ Attitude Control

for Quadrotor MAV's

Ole Falkenberg1 , Jonas Witt1 , Ulf Pilz2 , Uwe Weltin1 and Herbert Werner2
1
Hamburg University of Technology, Institute for Reliability Engineering,
Eiÿendorfer Str. 40, 21073 Hamburg, Germany
2
Hamburg University of Technology, Institute of Control Systems,
Eiÿendorfer Str. 40, 21073 Hamburg, Germany
{[Link],[Link],[Link],weltin,[Link]}@[Link]

Abstract This paper presents the results of modelling, parameter iden-

tication and control of the rotational axes of a quadrotor robot. The
modelling is done in Newton-Euler Formalism and has been published be-
fore. Contrarily, our method uses a Grey-Box-based, iterative parameter
identication approach, the results of which can easily be reproduced and
oers great accuracy. By neglecting nonlinear and cross-coupling eects,
only three to four parameters have to be identied per axis, depending
on the order of the motor dynamics. Based on the achieved results we
were able to design an aggressive H∞ attitude controller, which shows
superior performance to the normal PID-like controllers. With an anti-
windup compensator based on Riccatiequations we are able to show
exceptional input disturbance rejection, even with disturbances saturat-
ing the engines.

Keywords: micro unmanned aerial vehicle (MAV), robust control, sys-

tem identication, quadrotor

1 Introduction
Quadrotor robots are a popular research platform e.g. for autonomous navigation
or multi-agent-control. Famous projects are located in the GRASP-laboratory at
the University of Pennsylvania [1], the Flying Machine Area at the ETH Zürich
[2] or the STARMAC IIproject at the University of Stanford [3]. These projects
rely on commercial quadrotors, except for the latter. The Hamburg University
of Technology has developed an own quadrotor platform for research purposes,
which shows exceptional research capabilities because of its insight to the low-
est levels of interest. These insights reveal the drawbacks of common quadrotor
controller designs. Especially the nonideal performance of the common onboard
MEMSsensors in position estimation has driven many projects to use external
camera systems, which provide great accuracy for the estimation of the quadro-
tors states, but also limits the operational area to the laboratory. Unlike those
projects we do not rely on external camera systems, but only onboard sensors,
onboard cameravision and GPS, see [4]. Our research interest cover visionbased
2 Model Identication and H∞ Attitude Control for Quadrotor MAV's

navigation, formation control and aggressive manoevering. This paper focuses on

linear parameter identication and attitude control, though. Attitude control is
the most inner control loop and is crucial for a good overall ight performance.
Therefore an advanced controller design scheme like H∞ -loopshaping should be
favoured over traditional PID-control.
Beside the choice of a controller design scheme, an accurate model of the
dynamics to be controlled is another crucial aspect for controller performance.
The quadrotors physical equations of motion are well known and were published
in e.g. [5] or [6]. But only few has been published about actual parameter iden-
tication based on real experiments. The reason for this absence is the unstable
system dynamics of the quadrotor, which makes open-loop identication non-
practical. Additionally, a bad signal to noise ratio and cross-coupling eects make
for a challenging setup. The noise sources are mainly vibration and sensor noise.
Some rst-order studies on aerodynamics have been done in [7] which have
been extended in [3]. Infact, aerodynamics can have a major inuence on attitude
control. Especially the nonlinear variation of the induced velocity νi in forward
ight, which eectively changes the relation of power P and thrust T , can lead to
undesired behaviour. In [3], blade apping is specied as another major source
of nonlinearity causing undesired forces and torques. In this paper we present an
advanced attitude control scheme, which is based on a linearized model, feasible
in most common ight situations. Though its validity is limited:
 Nonlinearities caused by aerodynamics increase during forward ight with
relative movement speed.
 During fast descent the quadrotor can enter the so-called vortex-ring-state,
which cannot be estimated with normal momentum theory.
 The so-called Ground Eect augments thrust, if the quadrotor operates close
to the ground (≈ one rotor radius, see [5]).
However, it is obvious that a linear approach has advantages due to very
powerful controller design and analysis tools. Using the linearized system dy-
namics we show that each quadrotor axis can be identied separately with a
Grey-Boxed-based identication. With that model we have designed an aggress-
ive H∞ based attitude controller with anti-windup, which shows superior per-
formance to a common PID-like controller. Particularly the robustness against
input disturbances is considered, because main error sources in attitude control
are aerodynamic torques and wind gusts.

2 Modelling
In various papers the quadrotor dynamics are derived, see [3], [6] or [5]. In this
paper we are concerned with the results of our identication experiments in order
to identify the dynamics that have an actual impact.

2.1 Quadrotor Dynamics

In gure 1, a schematic representation of a quadrotor is shown: Each motor pro-
duces a torque τ with respect to the center of gravity. By individually varying
 Model Identication and H∞ Attitude Control for Quadrotor MAV's 3

T0
τ0 φ , ωxb
Ω0
b

P
T3 Ti T1
i
τ3 τ1
Ω1
b b

Ω3 xb
yb

T2 l θ , ωyb
zb
τ2 ψ , ωzb
Ω2 b

mg

Figure 1. Quadrotor forces T , torques τ , rotational velocities Ω , Eulerangles (φ, θ, ψ)

and coordinate system

the produced thrust of the motors, the overall-torque will be nonzero, permit-
ting the quadrotor to rotate about all three body-xed axes with an angular
velocity ω . If we assume a diagonal inertia matrix J (which is a quite good ap-
proximation, regarding the fact that the quadrotor is mechanically symmetric)
and following the formalism of Newton-Euler, one gets the following model:

Jyy − Jzz b b 1
ω̇xb = ωy ωz + τx
Jxx Jxx
Jzz − Jxx b b 1
ω̇yb = ωx ωz + τy (1)
Jyy Jyy
Jxx − Jyy b b 1
ω̇zb = ωx ωy + JR Ω̇r + τz
Jzz Jzz

with

ΩR = (Ω0 + Ω1 + Ω2 + Ω3 )
τx = l(T3 − T1 )
(2)
τy = l(T0 − T2 )
τz = (τ0 + τ2 − τ3 − τ1 )

The superscript b denotes the body frame. Ti is the thrust and τi the torque of
the individual motors. The distance between the motors to the center of gravity
is denoted as l, see gure 1.

2.2 Thrust and Torque

Motor thrusts Ti and torques τi are a direct result of the motors performing
work on the air.
4 Model Identication and H∞ Attitude Control for Quadrotor MAV's

The input of the motor plant is a nondimensional number b ∈ [0..1] which is

sent via a digital interface to the motor controller, which actually performs the
sensorless commutation for the brushless DC-motors. The output is the angular
speed Ω . With the help of a simple experiment, one can identify the relation of
Ω and b to be
Ω(b) = a0 · ba1 + a2 , (3)
with scalar coecients a0 , a1 , a2 ∈ R. The relation between motor speed and
thrust/torque can be found by any introductory aerodynamic text to be

T = CT ·ρAR2 Ω 2 (4)
3
Q = CQ ·ρAR Ω , 2
(5)

in which CT and CQ denotes the nondimensional thrust and torque coecients,

A the rotor disc area and R the rotor radius.

2.3 Motor Dynamics

We use brushless DC-motors with very low inductance. So the dynamics are
approximated by a rst-order ODE of the form

Ω̇ = κ1 Ω + κ2 Vin . (6)

The coecients κ1 and κ2 are motor specic. The motor dynamics are then
modelled by the transfer-function
|p|
Gm = . (7)
s+p
Note that eq. (7) is normalized to 1, so it only covers the dynamics of the motor.
Note also that Q(Ω)|ΩH = Q(Ω(b)|bH )|Ω(bH ) = Q(Ω(bH )) + qΓ (b − bH ), where
q = 2CQ ρAR3 ΩH is the linearization of eq. (5) around ΩH and Γ = a0 a1 baH1 −1
is the linearization of eq. (3) around bH .

3 Model for controller design

Both quadrotor dynamics and force/torque generation are nonlinear, so we need
to linearize if we want to use linear controller design techniques. We rst have
a look at quadrotor dynamics in eq. (1): As one can see, the nonlinear coupling
terms only act, if the quadrotor rotates about at least two axis at a time. The
moments of inertia for x and y are approximately the same, so the term for z
cancels, while Jzz is slightly higher. Experiments show, that the nonlinear eects
of the left over terms for x and y can be neglected in closed-loop operation, see
[8]. Looking at the yaw-equation, one sees the additional linear term Jr Ω̇ . This
torque directly relates to the time derivative of angular momentum L

L̇ = M = Jr Ω̇, (8)
 Model Identication and H∞ Attitude Control for Quadrotor MAV's 5

with Jr being the combined inertia of the rotors and motors. That means when
the motors change their speed, the time derivative of angular momentum will
be nonzero, causing a torque about the z -axis. Once the speed change of the
motors has stopped, L will be constant again. This results in a very fast response
on yaw with a fast decaytime, comparable to the D-part of a PID-controller.
By neglecting the nonlinear coupling terms and dening (Φ, Φ̇)T with output
Φ = (φ, θ, ψ)T , we get the following transfer function:
   1
0 0

φ Jxx s2
1
θ  0 0   
     Jyy s2 
τx
Φ ψ   0 0 (s+z) 
=   Jzz s2   
 · τy (9)
 φ̇  =  1
0 0
Φ̇    Jxx s 
τz
 θ̇    0 1
0


Jyy s
ψ̇ 0 0 (s+z)
Jzz s

Now, let's have a look at the origin of the forces and torques: Equation (4)
represents the equation of thrust in hover which depends quadratically on the
angular rotor velocity Ω . By substituting eq. (3) in (4), (5) and linearizing about
an operating point bh one gets:
CT CT
T = Γ̃ + qΓ · b + · · · and (10)
CQ R CQ R

Q= Γ̃ + qΓ · b + · · · , (11)

with

Γ̃ = CQ ρAR3 (a0 baH1 + a2 ) [1 − 2(a0 baH1 + a2 )a0 a1 bH ] . (12)

By neglecting higher order terms, (10) and (11) become linear equations of
the form y = mx + c.
With that we can write our torques in form of τ = (τx , τy , τz )T = Kb+const.
The overall model used for identication is then:
 1 
Jxx s2 0 0
 0 1
0   |p0 | −sTd 
   
Jyy s2
(s+z) 
 Kx s+p e · bx
 0 0 0
Φ Jzz s2   K |p0 | e−sTd · b 
= 1  ·  y s+p0 y  (13)
Φ̇  Jxx s 0 0 
 1  K̃z |p1 | −sTd
Kz s+p1 e · bz
 0 Jyy s 0 
(s+z)
0 0 Jzz s

Note that we introduced an additional factor K̃z and that the engine pole for
yaw is dierent than for roll and pitch. This is a direct result of the identication
procedure: While xing the pole at the same value as roll/pitch, the optimization
process needs more variance to converge to a reasonable solution in form of
a additional zero and pole. These can then be approximated by changing the
6 Model Identication and H∞ Attitude Control for Quadrotor MAV's

xed pole and gain, like denoted in eq. (13). A possible explanation is, that
experiments for roll and pitch have been performed in a test rig, which only
permits movement about one axis. In contrast, the experiment for yaw has been
performed in free ight, which allows more vibration.

4 GreyBox Parameter Identication

Due to unstable system dynamics of the quadrotor, the parameter identication
has been performed in closed-loop. A Pseudorandom Binary Sequence (PRBS)
γ of full length is used to excite the system. To improve the signal-to-noise ratio
we averaged the experiments over ten periods, respectively. As stated before all
three axes of the quadrotor have been identied separately. For that purpose
we xed the quadrotor in a test rig, allowing rotation only about one axis. The
identication for the yawaxis has been performed in free ight, though, because
of the lack of an approriate test rig. Even though an operator is still required
during free ight experiments for safety reasons, the operator can minimize pos-
sible noise introduced by reference steps, by limiting his commands to roll and
pitch. Recall from eq. (1), that a reference step on roll/pitch inuences yaw
only through the nonlinear coupling term. Because of that it is also possible to
average over several periods of γ , as stated above. The generic scheme of the
identication process is depicted in gure 2. Note that the used PRBS-signal is
persistently exciting of order  4, which is sucient for our purpose.

γ (y, ẏ)T
G(s)

u
C(s)

Figure 2. Indirect closedloop identication setup

The controller C(s) used during identication is a simple stabilizing, hand-

tuned PD-controller with known parameters. Using equation (13), the identi-
cation signal γ , the controller C and output data (y, ẏ)T , we can use non-
linear optimization to estimate our parameter vector ∆φ,θ = (p, J, Td )T and
∆ψ = (p, z, J, Td )T respectively. Note that we excluded K from the parameter
vector, because it can be calculated apriori: The aerodynamic coecients CQ
and CT have been calculated using Blade Element Momentum Theory (BEMT),
[9]. Note also, that you do not need this apriori knowledge, even though you
would not be able to distinguish between K and J in that case. See table 1 and
2 for results.
Figure 3 compares the simulation results of the identied models to a valida-
tion data set. For both, the simulation and the validation experiment, the same
stabilizing controllers were used to close the unstable openloop.
 Model Identication and H∞ Attitude Control for Quadrotor MAV's 7

Table 1. Results of parameter estimation Table 2. A-priori knowledge

φ θ ψ l 0.23m
J 0.0181 0.0196 0.0273 R 0.127m
p -27.106 -26.963 -7.171 CT 0.0157
z   -2.434 CQ 0.0015
Td 0.024 0.028 0.016 a0 701.2 1s
a1 0.8069
a2 −34.15 1s
bH 0.6

40 measured 20 measured
simulated simulated
20 10
Angle[°]

Angle[°]
0 0

−20 −10

−40 −20
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time [sec] Time [sec]

400 40
Angular Velocity[°/s]

Angular Velocity[°/s]
measured
simulated 20
200
0
0
−20

−200 −40
0 2 4 6 8 10 12 0 2 4 6 8 10 12
Time [sec] Time [sec]

(a) pitch (b) yaw

40 measured
simulated
20
Angle[°]

−20

−40
0 2 4 6 8 10 12
Time [sec]

400
Angular Velocity[°/s]

measured
simulated
200

−200
0 2 4 6 8 10 12
Time [sec]

(c) roll

Figure 3. Simulation of the nal linear models for a validation data set for the pitch
(a), yaw(b) and roll(c) axis. Due to unstable system dynamics simulation has been
performed in closed-loop with the same stabilizing controller as used for identication.

5 Controller Design
We have designed three separate attitude controllers using H∞ -loop-shaping
with the derived models of the quadrotors rotational axes from section 4. Recall
8 Model Identication and H∞ Attitude Control for Quadrotor MAV's

that our models, which will be denoted as G(s) respectively from now on, are
single-input, multiple output (SIMO). Accordingly, our controllers C(s) will be
MISO. Our design goals are zero steady-state error, fast rise-time with little
overshoot and good input-disturbance rejection.
The generalized plant used for the controller design is plotted in gure 4(a).
To deal with input disturbances (e.g. wind gusts) we introduced an input d addi-
tionally to the reference input r. The prelters Vd and Vr are scalar coecients
to weight between the channels with
 
V r1 0
Vr = . (14)
0 Vr2

The shaping lters WK (s) and WS (s) are rst order transfer functions with
 
WS (s) 0
WS (s) = . (15)
0 0

The remaining lter WT (s) is only a static gain to reduce overshoot. The band-
width of the controller has already been dened by WK (s) and WS (s).

r
Vr
d
Vd
zK
WK (s)
zT
WT (s)
y zS
G(s) WS (s) r y
P
R u
I G(s) ẏ
D
P(s) C(s)
u e
C(s)

(a) Generalized plant for H∞  controller (b) Structure of modied PID

design

Figure 4.(a) Generalized plant for H∞  controller design, (b) Structure of the modi-
ed PID-controller.

For the evaluation of our controller performance, we have designed a modied

PID controller, commonly used for quadrotor attitude control, see e.g. [10]. Its
parameters are tuned by means of optimizing the simulated system response.
Its generic structure is shown in gure 4(b). The main modication is the direct
feedback of the time derivative ẏ of the output y , because it can be measured
directly, so there is no need for additional dierentiation. Thus the D-part is just
a static gain.
In gure 5 the results of a 10◦ reference step with the corresponding con-
trol inputs are plotted. It can be seen that both controllers have a rise time
of ≈ 200ms with little overshoot, even though the PID-controllers have worse
reference tracking. In gure 6 the results of an input disturbance of 10% of the
maximum possible torque is shown. Here, our controller shows superior perform-
ance with a maximum displacement of only 4◦ , whereas the PID goes up to a
 Model Identication and H∞ Attitude Control for Quadrotor MAV's 9

displacement of 15◦ and it takes about four seconds longer to compensate the
disturbance.

12 0.1
H∞−design
10 PID−design
0.05
8

roll control output

roll [degree]

6 reference
H∞−design 0
4 PID−design

2
−0.05
0

−2 −0.1
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time [seconds] time [seconds]

(a) roll (b) roll

12
H∞−design
0.08
10 PID−design
0.06

pitch control output

8
0.04
pitch [degree]

6 reference 0.02
H∞−design
4 0
PID−design
−0.02
2
−0.04
0
−0.06

−2 −0.08
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time [seconds] time [seconds]

(c) pitch (d) pitch

Figure 5. (left): Response to a 10◦ step for roll and pitch. (right): Corresponding
control output

6 Anti-Windup Synthesis
This section deals with the design of an anti-windup compensator that accounts
for input saturation due to actuator constraints. For this purpose, the existing
linear H∞ controller is augmented with an anti-windup compensation scheme
which becomes active only when saturation occurs. The augmentation technique
also prevents a re-design of the existing control algorithm while limiting degrad-
ation of the closed-loop performance during saturation periods. The idea of the
anti-windup synthesis using Riccati equations is taken from [11].

6.1 Design Procedure

Figure 7 shows the generic anti-windup scheme, where G(s) is the plant given
in section 4 and K(s) is the H∞ controller described in section 5 designed to
stabilize the nominal plant while fullling desired performance specications.
The anti-windup compensator is given as Θ(s) which only becomes active once
the actuators are saturated. Note that the compensator Θ(s) has one input
10 Model Identication and H∞ Attitude Control for Quadrotor MAV's

16 0.1
reference input disturbance
14 H∞−design H∞−design
0.05
12 PID−design PID−design
0

roll control output

10

roll [degree]
8 −0.05

6 −0.1
4
−0.15
2
−0.2
0

−2 −0.25
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time [seconds] time [seconds]

(a) roll (b) roll

0.1
reference input disturbance
15 H∞−design 0.05 H∞−design
PID−design PID−design

pitch control output

0
pitch [degree]

10
−0.05

−0.1
5
−0.15

−0.2
0
−0.25
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time [seconds] time [seconds]

(c) pitch (d) pitch

Figure 6. (left): Response to a input disturbance equivalent to 10% of maximum engine

thrust for roll, pitch. (right): Correspondig control output

ũ = u−um and two outputs ud and yd which act on the controller output and the
controller input, respectively. A parametrization of the anti-windup compensator

r u + lin u um
K(s) G(s) y
y lin −

ud
+ −

ũ
yd Θ(s)
+ +

Figure 7. Generic anti-windup scheme

which is completely independent of the controller K(s) is given in [12]. A state-

space representation for the anti-windup compensator is given as

  ẋaw = (A + BF ) xaw + B ũ
Θ1
Θ= ∼ ud = F xaw (16)
Θ2 

yd = (C + DF ) xaw + Dũ,
 Model Identication and H∞ Attitude Control for Quadrotor MAV's 11

where A, B, C, D are the state space matrices of the plant G(s), F is a free
parameter and A + BF must be Hurwitz. The synthesis procedure is based on
minimizing a scalar γ > 0 such that the induced 2-norm of the transfer function
from the controller output ulin to the compensator output yd is smaller than a
positive real scalar γ , or
kTyd ulin ki,2 ≤ γ. (17)

The details of the construction of the anti-windup compensator Θ(s) are given
in [11], Theorem 1.

6.2 Experimental Results

Experimental validation of the anti-windup compensator was carried out with

the help of a test rig, where a input disturbance of 30% of the maximal motor
power was injected to the quadrotor robot. A comparison between the H∞ 
controller, described in section 5, with and without an anti-windup compensator
was performed. Figure 8 shows the results of experiments in terms of the roll
angle φ over time.

30 3
φ with anti-windup u with anti-windup
φ usat with anti-windup
20 2.5 u
usat

10 2
Control Input u [%]
Roll angle φ [◦ ]

0 1.5

−10 1

−20 0.5

−30 0

−40 −0.5
0 0.5 1 1.5 0 0.5 1 1.5
Time [sec] Time [sec]

(a) Roll angle φ (b) Control input for φ

Figure 8. Roll angle φ and corresponding control input of an input disturbance rejec-
tion experiment

From gure 8(a) it can be observed that the anti-windup compensation

scheme is able to stabilize the quadrotor robot despite the input disturbance
and to bring the roll angle φ asymptotically to the origin, whereas without anti-
windup compensator unstable behavior occurs.
In gure 8(b), it is shown how the anti-windup compensator inuences the
control input such that the linear region is reached after an actuator satura-
tion event has occured. After recovery from saturation, the control input always
stays close below the saturation boundary which also indicates good disturbance
rejection properties.
12 Model Identication and H∞ Attitude Control for Quadrotor MAV's

7 Conclusion
In this paper we have presented our results for linear model identication of
the quadrotors rotational axes under real, non-ideal conditions. The derived
models include a rst-order transfer function for the motors, while aerodynamic
coecients have been estimated using BEMT. The experimental results show
a very good correlation with real data, which conrms the proposed GreyBox
approach. In addition to the model identication, we designed robust H∞ at-
titude controllers. They have superior performance to commonly used modied
PID controllers, particularly with regard to input disturbance rejections, but
also to rise-time and reference tracking. Just to deal with input saturations, we
augmented the attitude controllers with an anti-windup compensator based on
Riccati-equations. This integrated anti-windup compensator leads to a further
improvement of the performance, e.g. under harsh environmental conditions.
In the future we will account for nonlinear aerodynamic eects to especially
optimize the control during forward ight.

References
1. Michael, N., Mellinger, D., Lindsey, Q., Kumar, V.: The GRASP Multiple Micro-
UAV Testbed. Robotics & Automation Magazine, IEEE 17(3) (2010) 5665
2. Lupashin, S., Schollig, A., Hehn, M., D'Andrea, R.: The Flying Machine Arena as of
2010. In: Robotics and Automation (ICRA), 2011 IEEE International Conference
on, IEEE (2011) 29702971
3. Homann, G., Huang, H., Waslander, S., Tomlin, C.: Quadrotor Helicopter Flight
Dynamics and Control: Theory and Experiment. In: Proceedings of the AIAA
Guidance, Navigation, and Control Conference. Volume 4. (2007) 44
4. Witt, J., Annighöfer, B., Falkenberg, O., Weltin, U.: Design of a High Performance
Quad-Rotor Robot Based on a Layered Real-Time System Architecture. Intelligent
Robotics and Applications (2011) 312323
5. Bouabdallah, S.: Design and Control of Quadrotors with Application to Autonom-
ous Flying. Ecole Polytechnique Federale de Lausanne (2007)
6. Beard, R.: Quadrotor Dynamics and Control. (2008)
7. Fay, G.: Derivation of the Aerodynamic Forces for the Mesicopter. (2001)
8. Falkenberg, O.: Robuste Lageregelung und GPS/INS-Integration eines autonomen
Quadrokopters. Hamburg University of Technology (2010)
9. Leishman, J.: Principles of Helicopter Aerodynamics. Cambridge Univ Pr (2006)
10. Tayebi, A., McGilvray, S.: Attitude Stabilization of a VTOL Quadrotor Aircraft.
Control Systems Technology, IEEE Transactions on 14(3) (may 2006) 562  571
11. Sofrony, J., Turner, M., Postlethwaite, I.: Anti-Windup Synthesis Using Riccati
Equations. International Journal of Control 80(1) (2007) 112128
12. Weston, P., Postlethwaite, I.: Linear Conditioning for Systems Containing Satur-
ating Actuators. Automatica 36(9) (2000) 13471354
13. Pilz, U., Gropengieÿer, W., Walder, F., Witt, J., Werner, H.: Quadrocopter Localiz-
ation Using RTK-GPS and Vision-Based Trajectory Tracking. Intelligent Robotics
and Applications (2011) 1221
14. Ljung, L.: System Identication: Theory for the User. Pearson Education (1998)
15. Witt, J.: Approximate Model Predictive Control for Nonlinear Multivariable Sys-
tems. In: Model Predictive Control, InTech (2010) 141166