Nonlinear Control and Filtering
Nonlinear Control and Filtering
Gerasimos G. Rigatos
Nonlinear Control
and Filtering Using
Differential Flatness
Approaches
Applications to Electromechanical
Systems
Studies in Systems, Decision and Control
Volume 25
Series editor
Janusz Kacprzyk, Polish Academy of Sciences, Warsaw, Poland
e-mail: [email protected]
About this Series
The series “Studies in Systems, Decision and Control” (SSDC) covers both new
developments and advances, as well as the state of the art, in the various areas of
broadly perceived systems, decision making and control- quickly, up to date and
with a high quality. The intent is to cover the theory, applications, and perspectives
on the state of the art and future developments relevant to systems, decision
making, control, complex processes and related areas, as embedded in the fields of
engineering, computer science, physics, economics, social and life sciences, as well
as the paradigms and methodologies behind them. The series contains monographs,
textbooks, lecture notes and edited volumes in systems, decision making and
control spanning the areas of Cyber-Physical Systems, Autonomous Systems,
Sensor Networks, Control Systems, Energy Systems, Automotive Systems, Bio-
logical Systems, Vehicular Networking and Connected Vehicles, Aerospace Sys-
tems, Automation, Manufacturing, Smart Grids, Nonlinear Systems, Power
Systems, Robotics, Social Systems, Economic Systems and other. Of particular
value to both the contributors and the readership are the short publication timeframe
and the world-wide distribution and exposure which enable both a wide and rapid
dissemination of research output.
Nonlinear Control
and Filtering Using
Differential Flatness
Approaches
Applications to Electromechanical Systems
123
Gerasimos G. Rigatos
Unit of Industrial Automation
Industrial Systems Institute
26504 Rion Patras
Greece
This book analyzes the design of nonlinear filters and nonlinear adaptive controllers,
using exact linearization which is based on differential flatness theory and differ-
ential geometry methods. The obtained filters exhibit specific advantages as they
outperform in terms of accuracy of estimation and computation speed of other
nonlinear filters. The obtained adaptive controllers can be applied to a wider class of
nonlinear systems with unknown dynamics and can assure reliable functioning
of the control loop under uncertainty and under varying operating conditions.
Moreover, the book analyzes differential flatness theory-based control and filtering
methods for distributed parameter systems. The book presents a series of application
examples to confirm the efficiency of the proposed control and filtering schemes for
various electromechanical systems. These include:
(i) Industrial Robotics: neuro-fuzzy adaptive control of multi-DOF robots and
underactuated robotic manipulators, observer-based neuro-fuzzy control of
multi-DOF robotic manipulators, state estimation-based control of multi-
DOF robots and underactuated robotic manipulators, state estimation-based
control of robots and mechatronic systems under disturbances and model
uncertainties.
(ii) Mobile robotics and vehicles: state estimation-based control of autonomous
vehicles, control of cooperating vehicles with the use of nonlinear filtering,
distributed fault diagnosis for autonomous vehicles, velocity control of
four-wheel vehicles, active vehicle suspension control, control of various
types of unmanned vehicles, such as AGVs, UAVs, USVs and AUVs.
(iii) Electric Power Generation: state estimation-based control of the PMSG
(Permanent Magnet Synchronous Generator), state estimation-based control
of the DFIG (Doubly-fed Induction Generator), state estimation-based
control and synchronization of distributed PMSGs.
(iv) Electric Motors and Actuators: neuro-fuzzy adaptive control of the DC
motor, neuro-fuzzy adaptive control of the Induction motor, state estimation-
based control of the DC motor, state estimation-based control of the
Induction Motor.
vii
viii Foreword
ix
x Preface
transmission systems. To this end, differential flatness theory is used for adaptive
control of the DC motor, for adaptive control of the induction motor, for state
estimation-based control of the DC motor, for state estimation-based control of
asynchronous electric motors and finally for observer-based adaptive fuzzy control
of microactuators (MEMS).
In Chap. 9, differential flatness theory is used to solve nonlinear filtering and
control problems which appear in power electronics such as voltage source con-
verters (VSC) and inverters, when these devices are used for connecting various
types of power generation units (AC and DC) to the grid. In particular, the chapter
proposes differential flatness theory for state estimation-based control of three-
phase voltage source converters, for state estimation-based control of voltage
inverters (finding application to the connection of photovoltaics to the electricity
grid), and finally for decentralized control and synchronization of distributed
inverters which are used to connect distributed DC power units with the electricity
network.
In Chap. 10, differential flatness theory is applied to nonlinear filtering and
control methods for internal combustion engines. The presented methods are con-
cerned with robust control and filtering for valves of diesel engines, with filtering
flatness-based of turbocharged diesel engines, with embedded adaptive control of
turbocharged Diesel engines, with embedded control and filtering of spark-ignited
engines, with embedded adaptive control of spark ignited engines and finally with
embedded control and filtering of the air-fuel ratio in combustion engines.
In Chap. 11, differential flatness theory-based methods for nonlinear filtering and
nonlinear control are applied to chaotic dynamical systems. Flatness-based adaptive
fuzzy control is proposed first for chaotic dynamical systems and manages to
modify the behavior of such systems without any knowledge of their dynamic
model. Next, differential flatness theory is proposed for developing a chaotic
communication system. A differential flatness theory-based Kalman Filtering
approach is proposed for performing blind equalization of the chaotic communi-
cation channel and for the synchronization between the transmitter and the receiver.
In Chap. 12, it is proposed to use differential flatness theory in nonlinear filtering
and control problems of distributed parameter systems, that is, systems which are
described by partial differential equations (of the parabolic, hyperbolic or elliptic
type). Methods of pointwise control of nonlinear PDE dynamics are introduced,
while methods for state estimation in nonlinear PDE models are developed. Such
results are applicable to communication systems (transmission lines, optical fibers,
and electromagnetic waves propagation), to electronics (Josephson junctions), to
manufacturing (heat diffusion control in the gas-metal arc welding proces), in
structural engineering (dynamic analysis of buildings under seismic waves,
mechanical structures subjected to vibrations, pendulum chains), in sensor networks
(fault diagnosis in distributed sensors monitoring systems with PDE dynamics), etc.
Finally, in Chap. 13, it is demonstrated that differential flatness theory is in the
background of other control methods, such as backstepping control and optimal
control. First, the chapter shows that differential flatness theory is in the background
of backstepping control. It is shown that to implement flatness-based control, it is
xiv Preface
not always necessary to transform the system into the trivial (linear canonical) form,
but it may suffice to decompose it into a set of flat subsystems associated with the
rows of its state-space description. The proposed flatness-based control method can
solve efficiently several nonlinear control problems, met for instance in robotic and
power generation systems. Moreover, the chapter shows how differential flatness
theory can be used for implementing boundary control of distributed parameter
systems. As a case study the control of a nonlinear wave PDE is considered.
The manuscript is suitable for teaching nonlinear control and nonlinear filtering
methods at late undergraduate and at postgraduate level. Such a course can be part
of the curriculum of several university departments, such as Electrical Engineering,
Mechanical Engineering, Computer Science, Physics, etc. The proposed book
contains teaching material which can be also used in a supplementary manner to the
content of undergraduate nonlinear control courses. The book can also serve per-
fectly the needs of a postgraduate course on nonlinear control and nonlinear esti-
mation methods.
Moreover, the book can be a useful reference for researchers, academic tutors,
and engineers that come against complex nonlinear control and nonlinear dynamics
problems. The book has both theoretical and practical value. It provides efficient
solutions to control problems associated with several electromechanical systems
met in industrial production, in transportation systems, in electric power generation,
and in many other applications. The control and filtering methods analyzed in the
book are generic and applicable to classes of systems, therefore one can anticipate
the application of the book’s methods to various types of electromechanical sys-
tems. I hope the book to become a useful reference not only for the industrial
systems, and the robotics and control community, but also for researchers and
engineers in the fields of signal processing, computational intelligence, and
mechatronics.
The author of this monograph would like to thank researchers and university staff
who through their constructive review and comments have contributed to the
improvement of the structure and content of this manuscript.
xv
Contents
xvii
xviii Contents
References. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 701
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 731
Acronyms
AC Alternate Current
AGV Automatic Ground Vehicle
ARE Algebraic Riccati Equation
ARMA Auto Regressive Moving Average
ARMAX Auto Regressive Moving Average with noise
AUV Autonomous Underwater Vessel
BER Bit Error Rate
CRLB Cramer-Rao Lower Bound
DC Direct Current
DFIG Doubly-Fed Induction Generator
DFKF Derivative-Free nonlinear Kalman Filtering
DOF Degrees of Freedom
DPS Distributed Parameter Systems
EKF Extended Kalman Filtering
KF Kalman Filtering
MSE Mean Square Error
MMSE Minimum Mean Square Estimator
MIMO Multi-Input Multi-Output
ODE Ordinary Differential Equation
PDE Partial Differential Equation
PMSG Permanent Magnet Synchronous Generator
SI Spark Ignited
SISO Single-input Single Output
SNR Signal-to-Noise Ratio
UAV Unmanned Aerial Vehicle
UKF Unscented Kalman Filtering
USV Unmanned Surface Vessel
VSC Voltage Source Converter
xxix
Chapter 1
Nonlinear Dynamical Systems
and Global Linearizing Control Methods
1.1 Introduction
This chapter overviews the basics of systems theory which can be used in the
modeling of nonlinear dynamics. To understand the oscillatory behavior of nonlinear
dynamical systems, benchmark examples of oscillators are given. Moreover, the fol-
lowing properties are analyzed: phase diagram, isoclines, attractors, local stability,
bifurcations of fixed points, and chaos properties.
Next, the chapter outlines differential geometry and Lie algebra-based control
as a predecessor to differential flatness theory-based control. First, the differential
geometric approach and the Frobenius theorem are analyzed. Next, the concept of
input–output linearization is introduced and its association to transformation of nor-
mal forms is explained. Additionally, the concept of input-state linearization is pre-
sented and the stages of its implementation are explained. Necessary and sufficient
conditions for applying input-state linearization and input–output linearization are
provided.
The main features characterizing the stability of nonlinear dynamical systems are
defined as follows [238, 555]:
1. Finite escape time: It is the finite time within which the state-vector of the nonlinear
system converges to infinity.
2. Multiple isolated equilibria: A linear system can have only one equilibrium to
which converges the state vector of the system in steady-state. A nonlinear system
can have more than one isolated equilibria (fixed points). Depending on the initial
state of the system, in steady-state the state vector of the system can converge to one
of these equilibria.
3. Limit cycles: For a linear system to exhibit oscillations it must have eigenvalues on
the imaginary axis. The amplitude of the oscillations depends on initial conditions.
In nonlinear systems one may have oscillations of constant amplitude and frequency,
which do not depend on initial conditions. These types of oscillations are known as
limit cycles.
4. Subharmonic, harmonic, and almost periodic oscillations: A stable linear system
under periodic input produces an output of the same frequency. A nonlinear system,
under periodic excitation can generate oscillations with frequencies that are several
times smaller (subharmonic) or multiples of the frequency of the input (harmonic).
It may also generate almost periodic oscillations with frequencies that are not nec-
essarily multiples of a basis frequency (almost periodic oscillations).
5. Chaos: A nonlinear system in steady-state can exhibit behavior that is not character-
ized as equilibrium, periodic oscillation, or almost periodic oscillation. This behavior
is characterized as chaos. As time advances the behavior of the system changes in a
random-like manner, and this depends on the initial conditions. Although the dynamic
system is deterministic it exhibits randomness in the way it evolves in time.
6. Multiple modes of behavior: It is possible for the same dynamical system to
exhibit simultaneously more than one of the aforementioned characteristics (1)–(5).
Thus, a system without external excitation may exhibit simultaneously more than
one limit cycles. A system receiving a periodic external input may exhibit harmonic
or subharmonic oscillations, or an even more complex behavior in steady state which
depends on the amplitude and frequency of the excitation.
ml 2 θ̈ = mgsin(θ )l − klv ⇒
ml 2 θ̈ = mgsin(θ )l − kl 2 θ̇ ⇒ (1.1)
ml θ̈ = mgsin(θ ) − kl θ̇
By defining the state variables x1 = θ and x2 = θ̇ one has the state-space description
ẋ1 = x2
(1.2)
ẋ2 = gl sin(x1 ) − k
m x2
ẋ1 = ẋ2 = 0⇒
(1.3)
x2 = 0 and x1 = nπ n = 0, ±1, ±2, · · ·
ẋ1 = x2
(1.4)
ẋ2 = gl sin(x1 )
ẋ1 = x2
(1.5)
ẋ2 = gl sin(x1 ) − mk x2 + 1
ml 2
T
1 d Vc (t)
Vc = Ic (s)⇒Ic (s) = sC Vc (s)⇒Ic (t) = C (1.6)
sC dt
Fig. 1.2 Tunnel diode
circuit
4 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
d I L (t)
VL (s) = sLI L (s)⇒VL (t) = L (1.7)
dt
I L (s) = IC + I R (1.8)
The following state variables are defined x1 = VC and x2 = i L . From Eq. (1.8) it
holds that
I L − IC − I R = 0⇒IC = I L − I R ⇒IC = x2 − h(x1 ) (1.10)
d IL
I L R − E + VL + VC = 0⇒E = VL + VC + i L R⇒E = L + VC + I L R⇒
dt
1 1
E = L ẋ2 + VC − x2 R⇒ẋ2 = [−VC − Rx2 + E]⇒ẋ2 = [−x1 − Rx2 + u]
L L
(1.11)
Additionally, from Eq. (1.6) it holds that
d Vc (t) 1
= Ic (t) (1.12)
dt C
By replacing Eq. (1.10) into Eq. (1.12) one obtains
d Vc (t) 1
= [x2 − h(x1 )] (1.13)
dt C
From Eqs. (1.11) and (1.13) one has the description of the system is state-space form
as
ẋ1 = C1 [x2 − h(x1 )]
(1.14)
ẋ2 = L1 [−x1 − Rx2 + u]
The associated equilibrium is computed from the condition ẋ1 = 0 and ẋ2 = 0 which
gives
−h(x1 ) + x2 = 0
(1.15)
−x1 − Rx2 + u = 0
and consequently
E−x1
−h(x1 ) + R =0
E−x1 (1.16)
x2 = R
Therefore, finally the equilibrium point is computed from the solution of the relation
E 1
h(x1 ) = − x1 (1.17)
R R
1.2 Characteristics of the Dynamics of Nonlinear Systems 5
One can also consider a model with nonlinear spring dynamics given by
The combination of a hardening spring, a linear viscous damping term, and a periodic
external force F = Acos(ωt) results in the Duffing oscillator
m ẍ + c ẋ + kx + ka 2 x 3 = Acos(ωt) (1.20)
The combination of a linear spring, a linear viscous damping term, a dry friction
term, and a zero external force generates the following oscillator model:
where
⎧
⎪
⎨ μk mg·sign(ẋ) if |ẋ > 0|
η(x, ẋ) = −kx if |ẋ| = 0 and |x|≤μs mg/k (1.22)
⎪
⎩
−μs mgsign(x) if ẋ = 0 and |x| > μs mg/k
ẋ1 = x2
(1.23)
ẋ2 = − mk x1 − mc x2 − m η(x, ẋ)
1
6 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
ẋ1 = x2
(1.24)
ẋ2 = − mk x1 − mc x2 − μk g
ẋ1 = x2
(1.25)
ẋ2 = − mk x1 − mc x2 + μk g
The method of the isoclines consists of computing the slope (ratio) between f 2 and
f 1 for every point of the trajectory of the state vector (x1 , x2 ).
f 2 (x1 , x2 )
s(x) = (1.27)
f 1 (x1 , x2 )
The case s(x) = c describes a curve in the x1 − x2 plane along which the trajectories
ẋ1 = f 1 (x1 , x2 ) and ẋ2 = f 2 (x1 , x2 ) have a constant slope.
The curve s(x) = c is drawn in the x1 − x2 plane and along this curve one also
draws small linear segments of length c. The curve s(x) = c is known as an isocline.
The direction of these small linear segments is according to the sign of the ratio
f 2 (x1 , x2 )/ f 1 (x1 , x2 ).
Example 1:
The following simplified pendulum equation is considered:
ẋ1 = x2
(1.28)
ẋ2 = −sin(x1 )
f 2 (x1 , x2 ) sin(x2 )
s(x) = ⇒s(x) = − (1.29)
f 1 (x1 , x2 ) x2
Setting s(x) = c it holds that the isoclines are given by the relation
1.3 Computation of Isoclines 7
2 c=1/2
c=1/3
c=1/4
1
x2
0
−1
−2
−3
−4
−5
−5 0 5
x1
1
x2 = − sin(x1 ) (1.30)
c
For different values of c one has the following isoclines diagram depicted in Fig. 1.4.
Example 2:
The following oscillator model is considered, being free of friction and with state-
space equations:
ẋ1 = x2
(1.31)
ẋ2 = −0.5x2 − sin(x1 )
s(x) = −0.5x2x−sin(x 1)
= c⇒ − 0.5x2 − sin(x1 ) = cx2 ⇒
2 (1.32)
(0.5 + c)x2 = sin(x1 )⇒x2 = − 0.5+c1
sin(x1 )
For different values of parameter c the isoclines are depicted in Fig. 1.5.
8 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
2 c6=0
c7=−1/6
c8=−1/4
1
x2
0
−1
−2
−3
−4
−5
−5 0 5
x1
The basic features that are important in the study of nonlinear dynamics are (i)
equilibria (fixed points), (ii) limit cycles, (ii) phase diagrams, (iv) periodic orbits, and
(v) bifurcations of fixed points [238, 555, 584, 605]. The definition of these features
is given through examples in the case of nonlinear dynamical systems models.
One can consider the nonlinear model with the two state variables, V and η. The
dynamics of the nonlinear model can be written as
dV
dt = f (V, η)
dη (1.33)
dt = g(V, η)
The phase diagram consists of the points on the trajectories of the solution of the
associated differential equation, i.e., (V (tk ), η(tk )).
At a fixed point or equilibrium it holds that f (V R , η R ) = 0 and g(V R , η R ) = 0.
The closed trajectories are associated with periodic solutions. If there are closed
trajectories the ∃T > 0 such that (V (tk ), η(tk )) = (V (tk + T ), η(tk + T )).
1.4 Basic Features in the Study of Nonlinear Dynamics 9
ẋ = Ax + Bu (1.34)
where ⎛ ∂ f1 ∂ f1 ⎞
∂x ∂x· · · ∂∂ xfn1
⎜ ∂ f12 ∂ f22 ⎟
· · · ∂∂ xfn2 ⎟
⎜ ∂ x1 ∂ x2
A = ∇x f = ⎜ ⎟ |x=x0 (1.35)
⎝··· ··· ··· ···⎠
∂ fn ∂ fn ∂ fn
∂ x1 ∂ x2 · · · ∂ xn
and ⎛ ∂ f1 ∂ f1 ∂ f1 ⎞
∂u ∂u· · · ∂u
⎜ ∂ f21 ∂ f22∂ f2 ⎟
n
· · · ∂u
⎜ ∂u 1 ∂u 2 ⎟
B = ∇u f = ⎜ n ⎟ | x=x (1.36)
⎝··· ··· ··· ···⎠ 0
∂ fn ∂ fn ∂ fn
∂u 1 ∂u 2 · · · ∂u n
The eigenvalues of matrix A define the local stability features of the system:
Example 1:
Assume the nonlinear system
d2x dx
2
+ 2x + 2x 2 − 4x = 0 (1.37)
dt dt
It holds that ẋ = 0 if f (x) = 0 that is (x1∗,1 , x2∗,1 ) = (0, 0) and (x1∗,2 , x2∗,2 ) = (2, 0).
Round the first equilibrium (x1∗,1 , x2∗,1 ) = (0, 0) the system’s dynamics is written as
10 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
01 ẋ1 01 x1
A = ∇x f = or = (1.39)
40 ẋ2 40 x2
The eigenvalues of the system are λ1 = 2 and λ2 = −2. This means that the fixed
point (x1∗,1 , x2∗,1 ) = (0, 0) is an unstable one.
Next, the fixed point (x1∗,2 , x2∗,2 ) = (2, 0) is analyzed. The associated Jacobian
matrix is computed again. It holds that
0 1
A = ∇x f = (1.40)
−4 −4
The eigenvalues of the system are λ1 = −2 and λ2 = −2. Consequently, the fixed
point (x1∗,2 , x2∗,2 ) = (2, 0) is a stable one.
The Lyapunov method analyzes the stability of a dynamical system without the need
to compute explicitly the trajectories of the state vector x = [x1 , x2 , . . . , xn ]T .
Theorem: The system described by the relation ẋ = f (x) is asymptotically stable in
the vicinity of the equilibrium x0 = 0 if there is a function V (x) such that
(i) V (x) to be continuous and to have a continuous first order derivative at x0
(ii) V (x) > 0 if x=0 and V (0) = 0
(iii) V̇ (x) < 0, ∀x=0.
The Lyapunov function is usually chosen to be a quadratic (and thus positive) energy
function of the system, however, there in no systematic method to define it.
Assume now that ẋ = f (x) and x0 = 0 is the equilibrium. Then the system is
globally asymptotically stable if for every ε > 0, ∃δ(ε) > 0, such that if ||x(0)|| < δ
then ||x(t)|| < ε, ∀t≥0.
This means that if the state vector of the system starts in a disk of radius δ then
as time advances it will remain at the same disk, as shown in Fig. 1.6. Moreover, if
lim t→∞ ||x(t)|| = x0 = 0, then the system is globally asymptotically stable.
Example 1: Consider the system
ẋ1 = x2
(1.41)
ẋ2 = −x1 − x32
The equilibrium point is (x1 = 0, x2 = 0). It holds that V (x) > 0∀(x1 , x2 )=
(0, 0) and V (x) = 0 for (x1 , x2 ) = (0, 0). Moreover, it holds that
1.4 Basic Features in the Study of Nonlinear Dynamics 11
Therefore, the system is asymptotically stable and lim t→∞ (x1 , x2 ) = (0, 0).
Example 2: Consider the system
1 2
V (x) = x + x22 (1.45)
2 1
The equilibrium point is x1 = 0, x2 = 0. It holds that
V̇ (x) = x1 ẋ1 + 2x2 ẋ2 = −x12 (1 + 2x1 x22 + 2x2 (x13 x2 ))⇒
(1.46)
V̇ (x) = −x12 < 0 ∀ (x1 , x2 )=(0, 0)
Therefore, the system is asymptotically stable and lim t→∞ (x1 , x2 ) = (0, 0).
Local stability of the previously described nonlinear model can be studied around the
associated equilibria. Local linearization can be performed around equilibria. Using
the set of differential equations ẋ = h(x) that describe the nonlinear system and
performing Taylor series expansion, that is, ẋ = h(x)⇒ẋ = h(x0 )|x0 + ∇x h(x −
x0 ) + · · · .
12 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
kα 2 x 3
where f (x1 , x2 ) = x2 and g(x1 , x2 ) = − mc x2 − mk x1 − m 1 . The fixed points of
this model are computed from the condition ẋ1 = 0 and ẋ2 = 0, and are found to be
(x1∗ , x2∗ ) = (0, 0). The Jacobian matrix ∇x h = M is given as
∂f ∂f
J= ∂ x1 ∂ x2 (1.48)
∂g ∂g
∂ x1 ∂ x2
At the fixed point (x1∗ , x2∗ ) = (0, 0) one obtains det (λI − J ) = λ2 + mc λ + mk which
for c2 < 4 km has imaginary roots with negative real part. Consequently, the system
is locally stable at the operating point and exhibits attenuating oscillations.
The eigenvalues of matrix M define stability around fixed points (stable or unstable
fixed point). To this end, one has to find the roots of the associated characteristic
polynomial which is given by det (λI − J ) = 0 where I is the identity matrix.
ẋ = Ax (1.50)
The eigenvalues of matrix A define the system dynamics. Some terminology asso-
ciated with fixed points is as follows:
A fixed point for the system of Eq. (1.50) is called hyperbolic if none of the
eigenvalues of matrix A has zero real part. A hyperbolic fixed point is called a saddle
if some of the eigenvalues of matrix A have real parts greater than zero and the
rest of the eigenvalues have real parts less than zero. If all of the eigenvalues have
negative real parts then the hyperbolic fixed point is called a stable node or sink. If
all of the eigenvalues have positive real parts then the hyperbolic fixed point is called
1.5 Phase Diagrams and Equilibria of Nonlinear Models 13
x2
0
−2
−4
−6
−8
−10
−15 −10 −5 0 5 10 15
x1
an unstable node or source. If the eigenvalues are purely imaginary then one has an
elliptic fixed point which is said to be a center.
Case 1: Both eigenvalues of matrix A are real and unequal, that is, λ1 =λ1 =0. For
λ1 < 0 and λ2 < 0 the phase diagram for z 1 and z 2 is shown in Fig. 1.7.
In case λ2 is smaller than λ1 the term eλ2 t decays faster than eλ1 t .
For λ1 > 0 > λ2 the phase diagram of Fig. 1.8 is obtained:
In the latter case there are stable trajectories along eigenvector v1 and unstable
trajectories along eigenvector v2 of matrix A. The stability point (0, 0) is said to be
a saddle point.
When λ1 > λ2 > 0 then one has the phase diagrams of Fig. 1.9:
Case 2: Complex eigenvalues:
Typical phase diagrams in the case of stable complex eigenvalues are given in
Fig. 1.10.
Typical phase diagrams in the case of unstable complex eigenvalues are given in
Fig. 1.11.
Typical phase diagrams in the case of imaginary eigenvalues are given in Fig. 1.12.
Case 3: Matrix A has nonzero eigenvalues which are equal to each other. The asso-
ciated phase diagram is given in Fig. 1.13.
14 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
0.5
x2
0
−0.5
−1
−1.5
−2
−30 −20 −10 0 10 20 30
x1
1
x2
−1
−2
−3
−30 −20 −10 0 10 20 30
x1
1.5 Phase Diagrams and Equilibria of Nonlinear Models 15
x2
0
−5
−10
−15
−10 −5 0 5 10
x1
0.4
0.2
x2
−0.2
−0.4
−0.6
−0.8
−1
−8 −6 −4 −2 0 2 4 6 8
x1 x 10
9
16 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
10
x2
0
−5
−10
−15
−20
−25
−15 −10 −5 0 5 10 15
x1
2
x2
−2
−4
−6
−8
−10
−15 −10 −5 0 5 10 15
x1
1.5 Phase Diagrams and Equilibria of Nonlinear Models 17
A nonlinear system can have multiple equilibria as shown in the following example.
Consider for instance the model of a pendulum under friction
ẋ1 = x2
(1.51)
ẋ2 = − gl sin(x1 ) − K
m x2
The associated phase diagram is designed for different initial conditions and is given
in Fig. 1.14.
For the previous model of the nonlinear oscillator, local linearization around
equilibria with the use of Taylor series expansion enables analysis of the local stability
properties of nonlinear dynamical systems
ẋ1 = f 1 (x1 , x2 )
(1.52)
ẋ2 = f 2 (x1 , x2 )
which gives
2
x2
−2
−4
−6
−30 −20 −10 0 10 20 30
x1
18 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
where
∂ f 1 (x1 ,x2 ) ∂ f 1 (x1 ,x2 )
α11 = ∂ x1 |x1 = p1 ,x2 = p2 α12 = ∂ x2 |x1 = p1 ,x2 = p2
∂ f 2 (x1 ,x2 ) ∂ f 2 (x1 ,x2 ) (1.54)
α21 = ∂ x1 |x1 = p1 ,x2 = p2 α22 = ∂ x2 |x1 = p1 ,x2 = p2
By omitting the higher order terms one can approximate the initial nonlinear system
with its linearized equivalent
ẏ1 = α11 y1 + α12 y2
(1.57)
ẏ2 = α21 y1 + α22 y2
ẏ = Ay (1.58)
where
∂ f1 ∂ f1
α11 α12 ∂ x1 ∂ x2 ∂f
A= = ∂ f2 ∂ f2 |x= p = ∂ x |x= p (1.59)
α21 α22 ∂ x1 ∂ x2
ẋ = f (x)⇒
(1.60)
ẋ1 x2
=
ẋ2 −sin(x1 ) − 0.5x2
1.5 Phase Diagrams and Equilibria of Nonlinear Models 19
There are two equilibrium points (0, 0) and (π, 0). The linearization around the
equilibria gives
0 1 0 1
A1 = A2 =
−1 −0.5 1 −0.5
(1.62)
with eigenvalues with eigenvalues
Consequently, the equilibrium (0, 0) is a stable focus (matrix A1 has stable complex
eignevalues) and the equilibrium (π, 0) is a saddle point (matrix A2 has an unstable
eigenvalue).
A dynamical system is considered to exhibit limit cycles when it admits the nontrivial
periodic solution
x(t + T ) = x(t) ∀ t≥0 (1.63)
for some T > 0 (the trivial periodic solution is the one associated with x(t) =
constant). An example about the existence of limit cycles is examined in the case of
the Van der Pol oscillator. The state equations of the oscillator are
ẋ1 = x2
(1.64)
ẋ2 = −x1 + ε(1 − x!2 )x2
Next the phase diagram of the Van der Pol oscillator is designed for three different
values of parameter ε, namely ε = 0.2, ε = 1 and ε = 5.0. In Figs. 1.15 and 1.16,
it can be observed that in all cases there is a closed trajectory to which converge all
curves of the phase diagram that start from points far from it.
To study conditions under which dynamical systems exhibit limit cycles, a second
order autonomous nonlinear system is considered next, given by
ẋ1 f 1 (x1 , x2 )
= (1.65)
ẋ2 f 2 (x1 , x2 )
20 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
1.5
0.5
x2
0
−0.5
−1
−1.5
−2
−2.5
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1
1
x2
−1
−2
−3
−2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1
1.5 Phase Diagrams and Equilibria of Nonlinear Models 21
Next, the following theorem defines the appearance of limit cycles in the phase
diagram of dynamical systems [238, 555].
Theorem 1 (Poincaré-Bendixson): If a trajectory of the nonlinear system of Eq. (1.65)
remains in a finite region Ω, then one of the following is true: (i) the trajectory goes
to an equilibrium point, (ii) the trajectory tends to a limit cycle, (iii) the trajectory is
itself a limit cycle.
Moreover, the following theorem provides a sufficient condition for the nonexis-
tence of limit cycles:
Theorem 2 (Bendixson): For the nonlinear system of Eq. (1.65), no limit cycle can
exist in a region Ω of the phase plane in which ∂∂ xf11 + ∂∂ xf22 does not vanish and does
not change sign.
As the parameters of the nonlinear model of the system are changed, the stability of
the equilibrium point can also change and also the number of equilibria may vary.
Values of these parameters at which the locus of the equilibrium (as a function of the
parameter) changes and different branches appear are known as critical or bifurca-
tion values. The phenomenon of bifurcation has to do with quantitative changes in
parameters leading to qualitative changes in the system’s properties [94, 137, 194,
258, 287].
Another issue in bifurcation analysis has to do with the study of the segments of the
bifurcation branches in which the fixed points are no longer stable but either become
unstable or are associated with limit cycles. The latter case is called Hopf bifurcation
and the system’s Jacobian matrix has a pair of complex imaginary eigenvalues.
√
For μ > 0 the dynamical system has two fixed points located at ± μ. One fixed
point is stable and is associated with the upper branch of the parabola. The other
fixed point is unstable and is associated with the lower branch of the parabola. The
value μ = 0 is considered to be a bifurcation value and the point (x, μ) = (0, 0) is a
bifurcation point. This particular type of bifurcation where one branch is associated
with fixed points and the other branch is not associated to any fixed points is known
as saddle-node bifurcation.
In pitchfork bifurcations the number of fixed points varies with respect to the values
of the bifurcation parameter. The dynamical system ẋ = x(μ − x 2 ) is considered.
1.6 Bifurcations in Nonlinear Dynamics 23
The associated fixed points are found by the condition ẋ = 0. For μ < 0 there is
one fixed point at zero which is stable. For μ = 0 there is still one fixed point at
zero which is still stable. For μ > 0 there are three fixed points, one at x = 0, one at
√ √
x = + μ which is stable and one at x = − μ which is also stable. The associated
phase diagrams and fixed points are presented in Fig. 1.19.
The bifurcations diagram is given next. The diagram shows how the fixed points
of the dynamical system vary with respect to the values of parameter μ. In the above
case, it represents a parabola in the μ − x plane as shown in Fig. 1.20.
24 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
Bifurcation of the equilibrium point means that the locus of the equilibrium on the
phase plane changes due to variation in the system’s parameters. Equilibrium x ∗ is
a hyperbolic equilibrium point if the real parts of all its eigenvalues are nonzero. A
Hopf bifurcation appears when the hyperbolicity of the equilibrium point is lost due
to variation of the system parameters and the eigenvalues become purely imaginary.
By changing the values of the parameters at a Hopf bifurcation, an oscillatory solution
appears [373, 576].
The stages for finding Hopf bifurcations in nonlinear dynamical systems are
described next. The following autonomous differential equation is considered:
dx
= f 1 (x, λ) (1.66)
dt
where x is the state vector x∈R n and λ∈R m is the vector of the system parameters. In
Eq. (1.66) a point x ∗ satisfying f 1 (x ∗ ) = 0 is an equilibrium point. Therefore, from
the condition f 1 (x ∗ ) = 0 one obtains a set of equations which provide the equilibrium
point as function of the bifurcating parameter. The stability of the equilibrium point
can be evaluated by linearizing the system’s dynamic model around the equilibrium
point and by computing eigenvalues of the Jacobian matrix. The Jacobian matrix at
the equilibrium point can be written as
∂ f 1 (x)
J f1 (x ∗ ) = |x=x ∗ (1.67)
∂x
and the determinant of the Jacobian matrix provides the characteristic equation which
is given as
ẋ1 = x2
(1.69)
ẋ2 = −x1 + (m − x12 )x1
Setting ẋ1 = 0 and ẋ2 = 0 one obtains the system’s fixed points. For m≤1 the
system has the fixed point (x1∗ , x2∗ ) = (0, 0). For m > 1 the system has the fixed
√ √
points (x1∗ , x2∗ ) = (0, 0), (x1∗ , x2∗ ) = ( m − 1, 0), (x1∗ , x2∗ ) = (− m − 1, 0). The
system’s Jacobian is
1.6 Bifurcations in Nonlinear Dynamics 25
0.5
x2
0
−0.5
−1
−1.5
−2
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5
x1
0
x
−1
−2
−3
−4
1 2 3 4 5 6 7 8 9 10 11
m
26 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
0 1
J= (1.70)
−1 + m − 3x12 0
If m≤1, the eigenvalues of the Jacobian at the fixed point (x1∗ , x2∗ ) = (0, 0) will
be imaginary. Thus the system exhibits Hopf bifurcations, as shown in Fig. 1.21. If
m > 1 then the fixed point (x1∗ , x2∗ ) = (0, 0) is an unstable one. On the other hand,
√
if m > 1 the eigenvalues of the Jacobian at the fixed points (x1∗ , x2∗ ) = ( m − 1, 0)
√
and (x1∗ , x2∗ ) = (− m − 1, 0) will be imaginary and the system exhibits again Hopf
bifurcations.
The bifurcations diagram is given next. The diagram shows how the fixed points
of the dynamical system vary with respect to the values of parameter μ. In the above
case it represents a parabola in the μ − x plane as shown in Fig. 1.22.
Before proceeding to the analysis of the differential geometric approach for the
control of nonlinear systems the following definitions are given [216, 238]:
Vector Field: A mapping f : D→R n where D⊂R n is a domain, is said to be a vector
field on D. A vector field is a n-dimensional column.
Covector Field: The transpose of a vector field is said to be a covector field. A
covector field is an n-dimensional row.
Inner Product: If f and w are, respectively, a vector field and a covector field on D,
then the inner product < w, f > is defined as
n
< w, f >= w(x) f (x) = wi (x) f i (x) (1.71)
i=1
∂h ∂h ∂h ∂h
dh = =[ , ,..., ] (1.72)
∂x ∂ x1 ∂ x1 ∂ xn
Lie Derivatives: Let h : D→R and f : D→R n . The Lie derivative of h, with respect
to f or along f , is denoted as L f h and is defined as
∂h
Lfh = f (x) (1.73)
∂x
This is the meaning of computing the derivatives of h along the trajectories of the
system ẋ = f (x). The previous notation can be generalized towards the computation
1.7 Predecessors of Differential Flatness Theory 27
of Lie derivatives of higher order with respect to the same vector field, or with respect
to a new vector field. Thus, one has
∂L f h
L g L f h(x) = g(x) (1.74)
∂x
while it also holds
L 0f h(x) = h(x)
∂h n
L f h(x) = ∂h
∂ x1 f 1 +
1 ∂h 2
∂ x2 f 2 + ··· + ∂ xn f n
··· (1.75)
∂ L k−1
L kf h(x) = L f L k−1
f = ∂x
f
f (x)
Using the Lie derivative notation one can arrive at the definition of the relative degree
for nonlinear dynamical systems. Consider for example the following single-input
single-output nonlinear system
ẋ = f (x) + g(x)u
(1.76)
y = h(x)
Lie Bracket: Let f and g be two vector fields on D⊂R n . The Lie Bracket of f and
g is written as [ f, g] and is a third vector field which is defined as
∂g ∂f
[ f, g](x) = f − g (1.78)
∂x ∂x
where ∂∂gx and ∂∂ xf are Jacobian matrices. The definition of the Lie Brackets can be
repetitively used as follows:
ad 0f g(x) = g(x)
ad f g(x) = [ f, g](x)
(1.79)
···
ad kf g(x) = [ f, ad k−1
f g](x)
[r1 f 1 + r2 f 2 , g1 ] = r1 [ f 1 , g1 ] + r2 [ f 2 , g2 ]
(1.80)
[ f 1 , r1 g1 + r2 g2 ] = r1 [ f 1 , g1 ] + r2 [ f 1 , g2 ]
28 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
is a subspace on R n . To each point x∈R n a subspace Δ(x) is assigned. One can refer
to this assignment using
Δ = span( f 1 , f 2 , . . . , f n ) (1.83)
and this is called a distribution.This implies that, Δ is the collection of all vector
spaces Δ(x) for xin D. The dimension of Δ(x) which is denoted as
r
g(x) = ci (x) f i (x) (1.85)
i=1
or equivalently
r
[ f i , f j ](x) = αi jk f k (x) ∀ i, j (1.88)
k=1
Equivalently, involutivity means that if one forms the Lie bracket of any pairs of vector
fields from the set [ f 1 , f 2 , . . . , f m ], the produced vector field can be expressed as a
linear combination of the original set of vector fields. It is noted that checking if a
set of vector fields [ f 1 , . . . , f m ] is involutive is equivalent to checking if
and the covector fields dh j (x) are linearly independent for all x∈D. Equivalently,
one can write
Δ⊥ = span{dh 1 , . . . , dh n−r } (1.92)
Finally, a key result in differential geometry is Frobenius’ theorem which states that:
A nonsingular distribution is completely integrable if it is involutive.
The Frobenius theorem is important for the formal treatment of feedback linearization
of n-order nonlinear systems. It provides a necessary and sufficient condition for the
solvability of a special class of partial differential equations. Consider for example
the set of first-order partial differential equations
30 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
∂h ∂h ∂h
∂ x1 f 1 + ∂ x2 f 2 + ∂ x3 f 3 =0
∂h ∂h ∂h (1.93)
∂ g1 g1 + ∂ x2 g2 + ∂ x3 g3 =0
[ f, g] = α1 f + α2 g (1.94)
which means that the Lie bracket of f and g can be expressed as a linear combination
of f and g. This is the involutivity condition on the vector fields [ f, g]. Thus, the
Frobenius theorem states that the set of vector fields [ f, g] is completely integrable,
if and only if it is involutive.
Next, the Frobenius theorem is restated, according to the previous definitions of
complete integrability and involutivity.
Frobenius theorem: Let f 1 , f 2 , . . . , f m be a set of linearly independent vector fields.
The set is completely integrable, if and only if it is involutive.
ẋ = f (x) + g(x)u
(1.95)
y = h(x)
ẋ = f (x) + g(x)u
(1.96)
y = h(x)
Proof : The proof is based on induction. For j = 0 Eq. (1.98) it holds because it
becomes equivalent to the definition of the relative degree of the system of Eq. (1.96).
Next, it is assumed that Eq. (1.98) holds for some j and it is proven for j + 1. Using
the Jacobi identity one has
for any real-valued function λ and any vector fields f and β. Taking λ = L kf ,
j
β = ad f g one obtains
In the previous equation, the first term of the right-hand side is equal to 0, that is,
L f L ad i g L kf h(x) = 0 because j + k + 1≤r − 1 or j + k < r − 1.
f
Moreover, using that Eq. (1.98) holds for j one has also
0 if 0≤ j + k + 1 < r − 1
L ad i g L k+1
f h(k) = (1.101)
f (−1) j L g L rf−1 h(x)=0 if j + k + 1 = r − 1
where denotes a nonzero element. Thus, the matrix is nonsingular. This provides
a proof for the Lemma, because if any of the two matrices appearing on the left side
of the equation had rank less than n, their product should be a singular matrix.
The previous lemma proves that r ≤n. Next, the main result of input–output lin-
earization will be formulated: A single-input single-output system of relative degree
r < n can be locally transformed into the canonical (normal) form. This is given in
the following theorem:
Theorem: The system of Eq. (1.96) is considered and it is assumed that it has relative
degree r < n in D. Then, for every x0 ∈D, there exist a neighborhood N of x0 and
smooth functions φ1 (x), . . . , φn−r (x), such that
restricted to N , is a diffeomorphism on N .
Proof : The distribution Δ = span{g} is nonsingular, involutive, and has dimension
1. By Frobenius’ theorem Δ is completely integrable. Consequently, for every x0 ∈D
1.7 Predecessors of Differential Flatness Theory 33
Equivalently,
Δ⊥ = span{dφ1 , . . . , dφn−1 } (1.108)
on N1 . Since
L g L if h(x) = 0, for 0≤i≤r − 2 (1.109)
and dh 1 , . . . , d L rf−2 h(x) are linearly independent, one can use h, . . . , L rf−2 h as part
of these n − 1 functions. In particular they can be denoted as φn−r +1 , . . . , φn−1 .
Since L g L rf−1 h(x)=0, one has
If L g L f h(x) is again zero for all x∈Ωx , additional differentiations with respect to
time can be attempted
34 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
L g L i−1
f h(x) =0 (1.115)
for some x = x0 in Ωx . If the above relation indeed holds then one can compute the
stabilizing control input that is finally applied to the system
1
u= (−L rf h + v) (1.116)
L g L i−1
f h(x)
By applying the control input of Eq. (1.116) to the equivalent description of the
system’s dynamics one gets
y (r ) = v (1.118)
The number r of the differentiations required for the input u to appear is the relative-
degree of the system. One has r ≤n, where n is the system’s order. If r = n then
input–output linearization becomes equivalent to the input-state linearization of the
system.
When the relative degree r is defined and r < n, the nonlinear system of Eq. (1.95)
can be transformed to the normal form using the new state vector [y, ẏ, . . . , y (r −1) ].
Let
μ = [μ1 , μ2 , . . . , μr ]T = [y, ẏ, . . . , y (r −1) ]T (1.119)
ψ̇ = w(μ, ψ) (1.121)
y = μ1 (1.122)
such that Eqs. (1.120) and (1.121) are verified. To show that a transformation φ is a
diffeomorphism it suffices to show that its Jacobian is invertible, i.e., the gradients
∇μi and ∇ψ j are linearly independent.
First, it will be shown that the gradients ∇μi are linearly independent which
means that the components of μ can serve as state variables. This result states that
the system’s output and its first r − 1 derivatives can serve as new state variables.
Next, it will be shown that other n − r vector fields ψ j can be found to complete
the new state vector, which means that the gradients ∇μi and ∇ψ j are also linearly
independent.
The following Lemma holds about the linear independence of the gradients ∇μi
[498]:
Lemma: If the relative degree of the system of Eq. (1.95) is r in the region Ω, then
the gradients ∇μ1 , ∇μ2 , . . ., ∇μr are also linearly independent in Ω.
Proof : It is first noted that in terms of μ Eq. (1.115) can be written as
Next, the proof is provided for the case that r = 3, however, it can be generalized to
higher values of r . It is assumed that there exist smooth functions αi (x), such that
everywhere in Ω
a1 ∇μ1 + a2 ∇μ2 + a3 ∇μ3 = 0 (1.125)
Multiplying Eq. (1.127) by the Lie bracket ad f g, and using Jacobi’s identity, as well
as the definition of the relative degree one has
0 = a1 L ad f g μ1 + a2 L ad f g μ2 =
= α1 [L f L g h − L g L f h] + α2 [L f L g − L g L f ]L f h = (1.128)
−a2 L g L 2f h
which from the definition of the relative degree implies that a2 = 0, everywhere in
Ω. Replacing a2 = 0 in Eq. (1.127) implies that
a1 ∇μ1 = 0 (1.129)
Multiplying Eq. (1.129) by ad 2f g and using the definition of the relative degree, as
well as the following result from Jacobi’s identity
L ad 2 g h = L ad f ad f g h = L f L ad f g h − L ad f g L f h =
f
= L f [L f L g h − L g L f h] − [L f L g − L g L f ]L f h = (1.130)
= L 2f L g h − 2L f L g L f h + L g L 2f h
one gets
0 = a1 L ad 2 g μ1 = a1 [L 2f L g h − 2L f L g L f h + L g L 2f h] = a1 L g L 2f h (1.131)
f
which again from the definition of the relative degree implies that a1 = 0, everywhere
in Ω. This finally confirms that the gradients ∇μi are linearly independent. This
completes the Lemma’s proof .
Next, it is shown that there exist n − r more functions ψ j to complete the coor-
dinates transformation. The proof about this is as follows:
It is needed to show that the gradient vectors ∇ψ j can be found such that the
∇μi and ∇ψ j are all linearly independent. At x0 Eq. (1.124) indicates that the first
r − 1 vectors ∇μi , (i = 1, 2, . . . , r − 1) which have been shown to be linearly
independent, are all within an hyperplane orthogonal to g. Since the dimension of
such hyperplane is n − 1, one can find
n − r = (n − 1) − (r − 1) (1.132)
vectors in that hyperplane that are linearly independent of the ∇μi and linearly
independent of each other. These vectors can be called ∇ψ j ( j = 1, 2, . . . , n − r )
and by definition they verify
1.7 Predecessors of Differential Flatness Theory 37
∇ψ j g = 0 1≤ j≤n − r (1.133)
Through Frobenius theorem it can also be assured that there exist n − 1 linearly
independent gradient functions ∇h k k = 1, . . . , n − 1 which satisfy ∇h k g = 0.
Furthermore, from Eq. (1.124) it is assured that μr is not in the hyperplane orthogonal
to g. Thus the gradients ∇μi (i = 1, . . . , r ) and ∇ψ j ( j = 1, . . . , n − r ) are
all linearly independent. Thus the Jacobian matrix of the transformation given in
Eq. (1.123) is invertible.
By continuity, the Jacobian remains invertible in a neighborhood Ω1 of x0 .
Redefining Ω as the intersection of Ω1 and Ω, the transformation φ defines a dif-
feomorphism in Ω. Hence, in Ω this transformation is a state transformation which
brings the nonlinear system in the form of Eq. (1.120), with
while in Eq. (1.121) the input u does not appear, since the ψ j verifies
L g ψ j (x) = 0 ∀x ∈ Ω (1.135)
From a practical point of view, to find the vector field ψ that completes the trans-
formation into a normal form requires the nontrivial step of solving the set of partial
differential equations given in Eq. (1.133).
ẋ = f (x) + g(x)u
(1.139)
y = h(x)
The “only if” part follows from the definition of input-to-state linearization while the
“if part” follows from applying the change of variables which is given by z = T (x) =
[h(x), L f h(x), . . . , L n−1 T
f h(x)] . The previously given Lemma 2 confirms that the
Jacobian matrix [ ∂∂Tx ](x) is nonsingular for all x∈Dx . Hence, for each x0 ∈Dx , there
is a neighborhood N of x0 such that T (x) is restricted to N , is a diffeomorphism. For
the theorem’s proof it is also necessary to show that the existence of h(x) satisfying
Eq. (1.140) is equivalent to the previous stated conditions 1 and 2.
Necessity: Suppose there is h(x) satisfying Eq. (1.140) shows that rank(G) = n.
Then D is nonsingular and has dimension n − 1. From Eq. (1.98), with k = 0 and
r = n one has
Furthermore,
Since rank(G) = n and dh(x)=0, it must be true that L ad n−1 g h(x)=0. Using the
f
Jacobi identity, it can be verified that L g L n−1
f h(x) =0, which completes the proof of
the theorem.
ż = Az + Bv (1.146)
where ⎛ ⎞ ⎛ ⎞
0 1 0 ··· 0 0
⎜0 0 1 ··· 0⎟ ⎜0⎟
⎜ ⎟ ⎜ ⎟
⎜· · · ··· ··· ··· · · ·⎟ ⎜ ⎟
A=⎜ ⎟ B = ⎜ · · ·⎟ (1.147)
⎜· · · ··· ··· ··· ⎟
· · ·⎟ ⎜ · · ·⎟
⎜ ⎜ ⎟
⎝0 0 0 ··· 1 ⎠ ⎝0⎠
0 0 0 ··· 0 1
The new state vector z is called linearizing state and the control input of Eq. (1.146) is
the linearizing control law. The transformed system dynamics given in Eqs. (1.146)
and (1.147) is in the so-called linear canonical form, which is both controllable and
observable.
The canonical form of Eq. (1.146) stands for a special case of the input–output
linearization, where the output function leads to a relative degree n. This means that
if a system is input–output linearizable with relative degree n, it is also input-state
linearizable. On the other hand, is a system is input-state linearizable with a first new
state z 1 representing the output, the system is also input–output linearizable with
relative degree n.
Consequently, the relation between input-state linearization and input–output lin-
earization can be summarized as follows:
Lemma: An nth order nonlinear system is input-state linearizable, if and only if, there
exists a scalar function z 1 (x) such that the system’s input–output linearization with
z 1 (x) as output has relative degree n.
40 1 Nonlinear Dynamical Systems and Global Linearizing Control Methods
However, the above lemma provides no guidance about how to find the desirable
output function z 1 (x).
It should be specified which are the nonlinear systems of the form of Eq. (3.44) which
can be subjected to input-state linearization. The above question stands for one of the
most important issues in feedback linearization theory. The necessary and sufficient
conditions for input–output linearization of a dynamical system have been stated in
the theorem given in Sect. 1.7.5. Next, the proof of this theorem is revisited [498].
Theorem: The nonlinear system of Eq. (3.44) with f (x) and g(x) being smooth vector
fields, is input-state linearizable, if and only if, there exists a region Ω such that the
following conditions hold:
• the vector fields [g, ad f g, . . . , ad n−1
f g] are linearly independent in Ω
• the set [g, ad f g, . . . , ad n−2
f g] is involutive in Ω
Prior to proving the theorem some remarks about the aforementioned conditions are
given:
The first condition can be interpreted as a controllability condition for the nonlin-
ear system of Eq. (1.139). For linear systems, the vector fields [g, ad f g, . . . , ad n−1
f g]
become [b, Ab, . . . , An−1 b] and independence implies that the associated matrix is
a full-rank one which is also invertible. It can be also shown that if a system’s linear
approximations in a closed connected region Ω∈R are all controllable, then under
some mild smoothness assumptions, the system can be driven from any point in Ω
to any point in Ω (the inverse does not hold because, a nonlinear system can be
controllable whereas its linear approximation may not be).
The involutivity condition is less intuitive. It is trivially satisfied for linear systems
(which have constant vector fields) but does not hold in general in the nonlinear case.
Next, to proceed to the Theorem’s proof, the following Lemma is analyzed:
Lemma: Let z(x) be a smooth function in a region Ω. Then, in Ω, the set of equations
L g z = L g L f z = · · · = L g L kf z = 0 (1.148)
is equivalent to
L g z = L g L ad f gz = · · · = L g L ad k gz = 0 (1.149)
f
L ad f g z = L f L g z − L g L f z = 0 − 0 = 0 (1.150)
1.7 Predecessors of Differential Flatness Theory 41
L ad 2 g z = L 2f L g z − 2L f L g L f z + L g L 2f z = 0 − 0 + 0 = 0 (1.151)
f
Repeating the procedure, it can be shown by induction that Eq. (1.150) implies
Eq. (1.151) for any k. Finally, using Jacobi’s identity in an inverse manner it can
be shown that from Eq. (1.151) one can arrive at Eq. (1.150).
Proof of Theorem: First the necessity condition is proven. It is assumed that there
exist a state transformation z = z(x) and an input transformation u = α(x) + β(x)v,
such that z and v satisfy Eq. (1.146). The first row of Eq. (1.146) can be rewritten as
∂z 1
ż 1 = ( f + gu) = z 2 (1.152)
∂x
Continuing, in a similar manner with the rest of the elements of the state vector one
has the following differential equations:
∂z 1 ∂z 1
f + gu = z 2 (1.153)
∂x ∂x
∂z 2
∂x f + ∂z
∂ x gu = z 3
2
(1.154)
··· ··· ···
∂z n ∂z n
f + gu = v (1.155)
∂x ∂x
Since z 1 , . . . , z n are independent of u, while v is not, from the above equations it
can be concluded that
L f z i = z i+1 i = 1, 2, . . . , n − 1 (1.157)
The above equations on z i can be compressed into one single equation concerning
only z 1 . Using the Lemma proven above and Eq. (1.156) one gets
∇z 1 ad kf g = 0 k = 0, 1, 2, . . . , n − 2 (1.158)
Moreover, by continuing as in the proof of the previous Lemma it can be shown that
∇z 1 ad n−1
f g = (−1)n−1 L g z n (1.159)
From Eqs. (1.158)–(1.160) one can conclude that the state vector fields [g, ad f g,
. . . , ad n−1
f g] are linearly independent. Indeed, if for some number i(i≤n − 1) there
existed scalar functions α1 (x), . . . , αi−1 (x), such that
i−1
ad if g = αk ad kf g (1.161)
k=0
n−2
∇z 1 ad n−1
f g= αk ad kf g = 0 (1.163)
k=0
L g z 1 = L ad f g z 1 = · · · = L ad n−2 g z 1 = 0 (1.164)
f
L g z 1 = L g L f z 1 = · · · = L g L n−2
f z1 = 0 (1.165)
Consequently, if one uses as new state variables the elements of the vector z =
[z 1 , L f z 1 , . . . , L n−1
f z 1 ] , the first n − 1 state equations give
T
ż k = z k+1 k = 1, . . . , n − 1 (1.166)
ż n = L nf z 1 + L g L n−1
f z1u (1.167)
1.7 Predecessors of Differential Flatness Theory 43
L g L n−1
f z 1 = (−1)
n−1
L ad n−1 g z 1 (1.168)
f
∇z 1 [g, ad f g, . . . , ad n−1
f g] = 0 (1.170)
∇z 1 ad if = 0 i = 0, 1, . . . , n − 2
(1.172)
∇z 1 ad n−1
f =0
L nf z 1
α(x) = −
L g L n−1
f z1 (1.173)
β(x) = 1
L g L n−1
f z1
(n)
The transformed control input that makes the linearized system z 1 = v be a Hurwitz
stable one, i.e., to have all its poles in the left complex semi-plane is given by
(n)
v = z 1,d − K 1 (z 1(n−1) − z 1,d
(n−1)
) − K 2 (z 1(n−2) − z 1,d
(n−2)
) − ···−
(1.174)
−K n−1 (ż 1 − ż 1,d ) − K n (z 1 − z 1,d )
(ri )
m
yi = L rfi h i + L g j L rfi −1 h i u j (1.175)
j=1
In the above description the linearized dynamical system is a square one, that is, the
number of its outputs is equal to the number of control inputs. Therefore, through
the followed procedure one succeeds at the same time linearization and decoupling
of the system’s dynamics. Thus, finally the linearized system is written in the form
Ỹ = F̃ + G̃u⇒
(1.177)
Ỹ = v
1.7 Predecessors of Differential Flatness Theory 45
and the feedback control law that makes the elements of the output vector Ỹ converge
to the desirable setpoints is u = G̃ −1 (v− F̃). The above system is said to have relative
degree (r1 , . . . , rm ) and the scalar r = r1 + . . . + rm is the total relative degree.
The control input ṽ is chosen to place the poles of the input–output linearized system
to the left complex semiplane. The control input [u̇ 1 , u 2 ]T is given by
u̇ 1
= Ẽ −1 (ṽ − B̃(x, u 1 )) (1.181)
u2
2.1 Introduction
First, the chapter analyzes flatness-based control for lumped parameter systems,
that is systems which are described by ordinary differential equations. The chapter
overviews the definition and properties of differential flatness and presents basic
classes of differentially flat systems. It is explained that all dynamical systems which
satisfy differential flatness properties can be transformed through a change of vari-
ables into the linear canonical form. The first section of the chapter presents examples
of single-input dynamical systems which are written into the linear canonical form by
using the differential flatness theory diffeomorphism and the design of the associated
feedback control loop is explained. The case of MIMO differentially flat dynamical
system is also examined. It is shown that differentially flat systems which admit static
feedback linearization can be transformed into the linear canonical form. Moreover, it
is shown that for MIMO differentially flat systems, that admit only dynamic feedback
linearization, it is again possible to succeed transformation to the linear canonical
form and subsequently to design state feedback controllers.
Next, the chapter examines flatness-based control for distributed parameter sys-
tems, that is, systems which are described by partial differential equations. Unlike
control of lumped parameter systems, distributed parameter systems control has been
less investigated. Such systems are described by partial differential equations and
the associated boundary conditions and play a critical role in several engineering
problems, such as vibrating structures, flexible-link robots, waveguides and optical
fibers, heat conduction, etc. Differential flatness theory enables the solution of such
control problems. A flatness-based control method for distributed parameter systems
proposes the decomposition of the desirable trajectory into a series of a reference flat
output and its derivatives, and enables to compute control commands that succeed
trajectory tracking. One can also consider flatness-based control for PDE systems
which are based on the transformation of the PDE model into a finite differences
models and the associated state-space description in a canonical form. The chapter
Differential flatness theory and flatness-based control were introduced in the late
1980s by Michel Fliess and coresearchers and since then they keep on being devel-
oped and on providing efficient solutions to advanced control and state estimation
problems [153].
The definition of a differentially flat system is as follows: A system ẋ = f (x, u)
with state vector x ∈ R n , input vector u ∈ R m where f is a smooth vector field, is
differentially flat if there exists a vector y ∈ R m in the form
such that
where h, φ and α are smooth functions. This means that the new system’s description
is given by the m algebraic variables yi , i = 1, 2, . . . , m. The definition of the flat
output given above was y = h(x, u, u̇, . . . , u (r ) ). If the flat output is exclusively a
function of the state vector x then the system is a 0-flat one. However, there may be
a need to express the flat output as a function of not only the state vector x but also
as a function of the control u and of its derivatives. For instance, the latter holds in
the case of dynamic feedback linearization and in the application of the so-called
dynamic extension. This means that the state vector of the system is extended by
considering as additional state variables the control inputs and its derivatives.
Equation (2.2) shows that the state vector of the differentially flat system and its
control inputs can be expressed as function of the flat output and of the flat output’s
derivatives. The basic question that arises in the study of differential flatness is
whether, given the differential equations that describe the nonlinear system dynamics
ẋ = f (x, u), there exists a function y = h() given by y = h(x, u, u̇, . . . , u (r ) ), such
that the state vector of the system x and the control input u can be expressed as
functions of y and of its derivatives, as in Eq. (2.2).
This problem was initially set by D. Hilbert in 1912, for the second-order Monge’s
equation
d2 y dy dz
= F(x, y, z, , ) (2.3)
dx2 dx dx
2.2 Definition of Differentially Flat Systems 49
The time variable is t ∈ R, the state vector is x(t) ∈ R n with initial conditions x(0) =
x0 , and the input variable is u(t) ∈ R m . Next, the main principles of differentially
flat systems are given [465, 535]:
The finite dimensional system of Eq. (2.4) can be written in the general form of an
ordinary differential equation (ODE), i.e., Si (w, ẇ, ẅ, . . . , w(i) ), i = 1, 2, . . . , q.
The entity w is a generic notation for the system variables (these variables are,
for instance, the elements of the system’s state vector x(t) and the elements of
the control input u(t)) while w(i) , i = 1, 2, . . . , q are the associated derivatives.
Such a system is said to be differentially flat if there is a collection of m functions
y = (y1 , . . . , ym ) of the system variables and of their time derivatives, i.e., yi =
φ(w, ẇ, ẅ, . . . , w(αi ) ), i = 1, . . . , m satisfying the following two conditions [152,
340, 362, 364, 422]:
50 2 Differential Flatness Theory and Flatness-Based Control
1. There does not exist any differential relation of the form R(y, ẏ, . . . , y (β) ) = 0
which implies that the derivatives of the flat output are not coupled in the sense of an
ODE, or equivalently it can be said that the flat output is differentially independent.
2. All system variables (i.e., the elements of the system’s state vector w and the
control input) can be expressed using only the flat output y and its time derivatives
wi = ψi (y, ẏ, . . . , y (γi ) ), i = 1, . . . , s. An equivalent definition of differentially
flat systems is as follows:
Definition: The system ẋ = f (x, u), x ∈ R n , u ∈ R m is differentially flat if there
exist relations
h : R n × (R m )r +1 → R m ,
φ : (R m )r → R n and (2.5)
ψ : (R m )r +1 → R m
such that
y = h(x, u, u̇, . . . , u (r ) ),
x = φ(y, ẏ, . . . , y (r −1) ), and (2.6)
u = ψ(y, ẏ, . . . , y (r −1) , y (r ) ).
This means that all system dynamics can be expressed as a function of the flat output
and its derivatives; therefore, the state vector and the control input can be written as
Next, an example is given to explain the design of a differentially flat controller for
finite dimensional systems of known parameters.
Example 1: Flatness-based control for a nonlinear system of known parameters [263].
Consider the following model:
ẋ1 = x3 − x2 u
ẋ2 = −x2 + u (2.8)
ẋ3 = x2 − x1 + 2x2 (u − x2 )
x22
The flat output is chosen to be y1 = x1 + 2 . Thus one gets:
x2
y1 = x1 + 22
y2 = ẏ1 = (x3 − x2 u) + x2 (u − x2 ) = x3 − x22 (2.9)
y3 = ẏ2 = ÿ1 = x2 − x1 + 2x2 (u − x2 ) − 2x2 (u − x2 ) = −x1 + x2
(3)
v = ẏ3 = y1 = −x3 + x2 u − x2 + u = −x2 − x3 + u(1 + x2 )
It can be verified that property (1) holds, i.e., there does not exist any differential
relation of the form R(y, ẏ, . . . , y (β) ) = 0, and this implies that the derivatives of
2.2 Definition of Differentially Flat Systems 51
the flat output are not coupled. Moreover, it can be shown that property (2) also
holds i.e., the components w of the system (elements of the system’s state vector and
control input) can be expressed using only the flat output y and its time derivatives
wi = ψi (y, ẏ, . . . , y (γi ) ), i = i, . . . , s.
(3)
For instance to calculate x1 with respect to y1 , ẏ1 , ÿ1 and y1 the relation of ÿ1 is
used, i.e.,
x12 + 2x1 (1 + ÿ1 ) + ÿ12 − 2y1 = 0 (2.10)
√
√ solutions are derived, i.e., x1 = −(1+ ÿ1 − 1 + 2(y1 + ÿ1 ))
from which two possible
and x1 = −(1 + ÿ1 + 1 + 2(y1 + ÿ1 )). Keeping the biggest out of these two solu-
tions one obtains: √
x1 = −(1 + ÿ1 ) + 1 + 2(y1 + ÿ1 )
x2 = ÿ1 + x1
x3 = ẏ1 + ÿ12 + 2x1 ÿ1 + x12 (2.11)
(3)
y1 + ÿ12 + ÿ1 + ẏ1 +x1 +2x1 ÿ1 +x12
u= 1+x1 + ÿ1
The computation of the equivalent model of the system in the linear canonical form
is summarized as follows: By finding the derivatives of the flat output one gets a
set of equations which can be solved with respect to the state variables and the
control input of the initial state-space description of the system. First, the binomial
of variable x1 given in Eq. (2.10) is solved providing x1 as a function of the flat output
and its derivatives. Next, using the expression for x1 and Eq. (2.11), state variable
x2 is also written as a function of the flat output and its derivatives. Finally, using
the expressions for both x1 and x2 and Eq. (2.11), state variable x3 is written as a
function of the flat output and its derivatives. Thus one can finally express the state
vector elements and the control input as function of the flat output and its derivatives,
which completes the proof about differential flatness of the system.
From Eq. (2.11) it can be concluded that the initial system of Eq. (2.8) is indeed
differentially flat. Using the flat output and its derivatives, the system of Eq. (2.8) can
be written in Brunovsky (canonical) form:
⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
y 010 y1 0
d ⎝ 1⎠ ⎝
y2 = 0 0 1⎠ ⎝ y2 ⎠ + ⎝0⎠ v (2.12)
dt y 000 y3 1
3
where the new control input is v = f (x) + g(x)u. Therefore, a transformation of the
system into a linear equivalent description is obtained and then it is straightforward to
design a controller based on linear control theory. Thus, given the reference trajectory
[x1∗ , x2∗ , x3∗ ]T , one can find the transformed reference trajectory [y1∗ , ẏ1∗ , ÿ1∗ ]T and
select the appropriate control input v that succeeds tracking of the reference setpoints.
Knowing v, the control u of the initial system can be found. Knowing v the control
input that is actually applied to the system is u = g −1 (x)[v − f (x)].
52 2 Differential Flatness Theory and Flatness-Based Control
It is noted that for linear systems, the property of differential flatness is equivalent
to that of controllability. Next, two examples of differentially flat MIMO dynamical
systems are given. It is shown that the definition of the differential flat outputs also
enables to transform the system into the Brunovksy (canonical) form:
Example 2: Differential flatness of a nonlinear spring–damper–mass system which
consists of two masses.
The spring–damper–mass model is described in Fig. 2.1. The dynamic equations
of the model are given by [89]
where M1 and M2 are the masses of the system, x(t) = [x1 , ẋ1 , x2 , ẋ2 ]T is the state
vector which has as elements the positions and the velocities of the two masses, and
f K 1 (x) and f K 2 (x) are spring forces defined by the following equations
f K 1 (x) = K 10 + ΔK 1 x13
(2.14)
f K 2 (x) = K 20 + ΔK 2 (x2 − x1 )3
y1 = x1 y2 = x2 (2.16)
Obviously, it holds
x1 = y1 x2 = y2
(2.17)
ẋ1 = ẏ1 ẋ2 = ẏ2
Thus, it is observed that the state variables of the spring-damper-mass model can be
expressed as functions of the flat outputs and of the associated derivatives. Moreover,
the following relations can be obtained about the control inputs of the model
M1 ẍ1 =
−(K 10 x1 + ΔK 1 x13 ) − (b10 ẋ1 + Δb1 ẋ12 )+
(2.18)
+(K 20 (x2 − x1 ) + ΔK 2 (x2 − x1 )3 )+
+(b20 (ẋ2 − ẋ1 ) + Δb2 (ẋ2 − ẋ1 )2 ) + u 1 + 0.2u 2
M2 ẍ2 =
−(K 20 (x2 − x1 ) + ΔK 2 (x2 − x1 )3 ) − (b20 (ẋ2 − ẋ1 )+ (2.19)
Δb2 (ẋ2 − ẋ1 )2 + 0.25u 1 + u 2
Using the definition of the flat outputs in Eqs. (2.18) and (2.19), one obtains
M1 ÿ1 =
−(K 10 y1 + ΔK 1 y13 ) − (b10 ẏ1 + Δb1 ẏ12 )+
(2.20)
+(K 20 (y2 − y1 ) + ΔK 2 (y2 − y1 )3 )+
+(b20 ( ẏ2 − ẏ1 ) + Δb2 ( ẏ2 − ẏ1 )2 ) + u 1 + 0.2u 2
M2 ÿ2 =
−(K 20 (y2 − y1 ) + ΔK 2 (y2 − y1 )3 ) − (b20 ( ẏ2 − ẏ1 )+ (2.21)
+Δb2 ( ẏ2 − ẏ1 )2 + 0.25u 1 + u 2
or equivalently
Defining the matrix of the coefficients of the control inputs and its inverse as
1 0.2 1.0526 −0.2105
A= A−1 = (2.24)
0.25 1 −0.2632 1.0526
one obtains the following relation about the control inputs and the flat outputs defined
in Eq. (2.16)
u1 1.0526 −0.2105 q1
= (2.25)
u2 −0.2632 1.0526 q2
ÿ1 =
M1 {−(K 10 y1 + ΔK 1 y1 ) − (b10 ẏ1 + Δb1 ẏ1 )}+
1 3 2
(2.26)
M1 {(K 20 (y2 − y1 ) + ΔK 2 (y2 − y1 ) )+
1 3
ÿ2 =
M2 {−(K 20 (y2 − y1 ) + ΔK 2 (y2 − y1 ) )−
1 3
(2.27)
−(b20 ( ẏ2 − ẏ1 ) + Δb2 ( ẏ2 − ẏ1 ) )} + 0.25 M12 u 1 +
2 1
M2 u 2
with
f 1 (y1 , ẏ1 , y2 , ẏ2 ) = f 1 (x) =
M1 {−(K 10 y1 + ΔK 1 y1 ) − (b10 ẏ1 + Δb1 ẏ1 )}+
1 3 2
(2.30)
M1 {(K 20 (y2 − y1 ) + ΔK 2 (y2 − y1 ) )+
1 3
Thus, one obtains the following canonical form description for the nonlinear system
⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
ẋ1 0 1 0 0 x1 0 0
⎜ẋ2 ⎟ ⎜0 0 0 0⎟ ⎜ x 2 ⎟ ⎜1 0⎟
⎜ ⎟=⎜ ⎟⎜ ⎟ + ⎜ ⎟ v1 (2.34)
⎝ẋ3 ⎠ ⎝0 0 0 1⎠ ⎝x3 ⎠ ⎝0 0⎠ v2
ẋ4 0 0 0 0 x4 0 1
where
v1 = f 1 (x) + g11 (x)u 1 + g12 (x)u 2
(2.35)
v2 = f 2 (x) + g21 (x)u 1 + g22 (x)u 2
or equivalently
v = F(x) + G(x)u (2.36)
where
v1 u1 f 1 (x)
v= , u= ,F = ,
v2 u2 f 2 (x)
(2.37)
g (x) g12 (x)
G = 11
g21 (x) g22 (x)
Example 3: Differential flatness of the VTOL aircraft (vertically take-off and landing
aircraft).
The dynamic model of the vertically take-off and landing aircraft (Fig. 2.2) is
described by the following set of equations [561]
ẍ = u 1 sin(θ ) − εu 2 cos(θ )
z̈ = u 1 cos(θ ) + εu 2 sin(θ ) − 1 (2.38)
θ̈ = u 2
56 2 Differential Flatness Theory and Flatness-Based Control
ε ÿ1
x = y1 −
,
ÿ1 +( ÿ2 +1)2
2
z = y2 −
ε ÿ2 +1 , (2.39)
ÿ12 +( ÿ2 +1)2
θ = tan −1 ( ÿ2ÿ+1
1
)
Since all state variables are expressed as functions of the flat outputs and their deriv-
atives, one can also write the control inputs u 1 and u 2 as functions of the flat outputs
and their derivatives. This confirms that the model of the vertically take-off and
landing aircraft is differentially flat. The system can be written in the Brunovsky
form
(4) (4)
y1 = v1 , y2 = v2 (2.40)
The previous Brunovsky-form model of the vertically take-off and landing aircraft
can be also written using state-space equations
⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
ẏ1,1 0 1 0 0 0 0 0 0 y1,1 0 0
⎜ ẏ1,2 ⎟ ⎜0 0 1 0 0 0 0 0⎟ ⎜ y1,2 ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟
⎜ ẏ1,3 ⎟ ⎜0 0 0 1 0 0 0 0⎟ ⎜ ⎟ ⎜ 0⎟
⎜ ⎟ ⎜ ⎟ ⎜ y1,3 ⎟ ⎜0 ⎟
⎜ ẏ1,4 ⎟ ⎜0 0 0 0 0 0 0 ⎟ ⎜
0⎟ ⎜ y1,4 ⎟ ⎜ 0⎟
⎜ ⎟ ⎜ ⎟ + ⎜1 ⎟ v1
⎜ ẏ2,1 ⎟ = ⎜0 0 0 0 0 1 0 ⎟ ⎜
0⎟ ⎜ y2,1 ⎟ ⎜
⎟ 0⎟
(2.43)
⎜ ⎟ ⎜ ⎜0 ⎟ v2
⎜ ẏ2,2 ⎟ ⎜0 0 0 0 0 0 1 0⎟ ⎜ y2,2 ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟
⎝ ẏ2,3 ⎠ ⎝0 0 0 0 0 0 0 1⎠ ⎝ y2,3 ⎠ ⎝0 0⎠
ẏ2,4 0 0 0 0 0 0 0 0 y2,4 0 1
ẋ = f (x) x ∈ X ⊂ R n (2.44)
ẋ = f (x, u) (2.46)
58 2 Differential Flatness Theory and Flatness-Based Control
ẋ = f (x, u)
(2.47)
u̇ = v
In the latter case, in place of the state-space X one has the state-space X × U .
Assume now that the solution of ẋ = f (x, u) is the function t → (x(t, u(t))) which
has values in X × U , such that ẋ(t) = f (x(t), u(t)) ∀ t ≥ 0. Moreover, one can
consider the infinite function
relation
ξ̇ (t) = ( f (x(t), u(t)), u̇(t), ü(t), . . .) ∀ t ≥ 0 (2.49)
According to the above, the dynamical system ẋ = f (x, u) is defined by the space
X × U × Rm ∞ and an infinite number of vector fields F on that space. As in the case
If h() depends on a finite number of variables then the above sum becomes finite
too. For the above type of functions, one can define the so-called Fréchet topology
which implies computations on smooth functions which are defined in k copies of
R m , for k being sufficiently large.
After a change of variables, a differentially flat system is written in the canonical
Brunovsky form. Thus, it becomes equivalent to the trivial system Rm ∞ , F with
m
(1) (2)
coordinates y, y , y , . . . , and the vector fields
If the two systems are equivalent, the endogenous transformation takes the form
where ū is the abbreviated notation ū = (u, u (1) , . . . , u (k) ). Similarly, one can define
the inverse endogenous transformation
About vector fields F and G which are related to the derivatives of x and y, respec-
tively, one has
f (φ(y, v̄), α(y, v̄)) = Dφ(y, v̄) · ḡ(y, v̄) (2.59)
Example 1: The VTOL aircraft (Vertical Take-off and Landing Aircraft) studied in
Sect. 2.2.2 is revisited. State-space transformation for the VTOL model provides an
example of equivalence between dynamical systems descriptions.
The initial dynamic model of the VTOL aircraft is given by
ẍ = −u 1 sin(θ ) + εu 2 cos(θ )
z̈ = u 1 cos(θ ) + εu 2 sin(θ ) − 1 (2.63)
θ̇ = u 2
ÿ1 = −ξ sin(θ )
ÿ2 = ξ cos(θ ) − 1 (2.64)
θ̇ = u 2
where ξ and θ stand for the new control inputs. Indeed, choosing
and with the abbreviated notations given above one can define the functions of the
transformed state vector and of the transformed control input Y = ψ(X, Ū ) and
V = β(X, Ū ), as
⎛ ⎞
x − εsin(θ )
⎜ z + εcos(θ ) ⎟
⎟ and β(X, Ū ) = u 1 − εθ̇
2
ψ(X, U ) = ⎜
⎝ẋ − εθ̇ cos(θ )⎠ (2.66)
θ
ż − εθ̇ sin(θ )
2.3 Properties of Differentially Flat Systems 61
which finally enable to obtain the endogenous transformation Ψ . The inverse trans-
formation Φ is given by X = φ(Y, V̄ ) and U = α(Y, V̄ ) according to the following
⎛ ⎞
y1 + εsin(θ )
⎜ y1 − εcos(θ ) ⎟
⎜ ⎟
⎜ ẏ1 + εθ̇ cos(θ )⎟ ξ + εθ̇ 2
⎜
φ(Y, V̄ ) = ⎜ ⎟
⎟ and α(Y, V̄ ) = (2.67)
⎜ ẏ2 − εθ̇ sin(θ )⎟ θ̈
⎝ θ ⎠
θ̇
Continuing with the vertical take-off and landing aircraft (VTOL) and using the
dynamics of the aircraft which has been defined in Eq. (2.63), it is found that this
system admits as flat output
To express the state variables and the control input as functions of the flat output and
its derivatives (this enables to find also the elements of the inverse transformation of
Φ, that is X = φ( ȳ) and U = α(y)), the following implicit relations are used:
Solving the above set of equations with respect to x, z, and θ , one gets
ÿ1
x = y1 ±ε
One has simply to differentiate so as to obtain ẋ, ż, θ̇ , and u as functions of the deriv-
atives of the flat output y. Singularity of the system is avoided if ÿ12 + ( ÿ2 + 1)2 = 0.
Next, the system is transformed to the linear canonical form described in Eq. (2.43)
through dynamic feedback linearization. The feedback control consists of the fol-
lowing elements
ξ̇ = −v1 sin(θ ) + v2 cos(θ ) + ξ θ̇ 2
u 1 = ξ + εθ̇ 2 (2.71)
u 2 = − ξ (v1 cos(θ ) + v2 sin(θ ) + 2ξ̇ (θ̇))
1
62 2 Differential Flatness Theory and Flatness-Based Control
(4)
y1 = v1
(2.72)
y2(4) = v2
The only singularity which may appear in this feedback control loop is when ξ = 0
that is ÿ12 + ( ÿ2 + 1)2 = 0 (which practically means that the aircraft is in free fall).
Example 2: The underactuated vessel dynamics. The model of an autonomous hov-
ercraft provides another example about equivalence between an initial complicated
nonlinear description of its dynamics and the linear canonical (Brunovsky form)
[461].
The hovercraft model is obtained from the generic ship’s model, after setting
specific values for the elements of the inertia and Coriolis matrix and after reducing
the number of the available control inputs. The state-space equation of the nonlinear
underactuated hovercraft model (Fig. 2.3) is given by
ẋ = ucos(ψ) − vsin(ψ)
ẏ = usin(ψ) + vcos(ψ)
ψ̇ = r
(2.74)
u̇ = v · r + τu
v̇ = −u · r − βv
ṙ = τr
where x and y are the cartesian coordinates of the vessel, ψ is the orientation angle,
u is the surge velocity, v is the sway velocity, and r is the yaw rate. Coefficient β
is a function of elements of the inertia matrix and hydrodynamic damping matrix
of the vessel. The control inputs are the surge force τu and the yaw torque τr . The
hovercraft’s model is also written in the matrix form:
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
ẋ ucos(ψ) − vsin(ψ) 00
⎜ ẏ ⎟ ⎜ usin(ψ) + vcos(ψ) ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ψ̇ ⎟ ⎜ r ⎟ ⎜0 0⎟ τu
⎜ ⎟=⎜ ⎟+⎜ ⎟ (2.75)
⎜ u̇ ⎟ ⎜ vr ⎟ ⎜1 0⎟ τr
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ v̇ ⎠ ⎝ −ur − βv ⎠ ⎝0 0⎠
ṙ 0 01
It holds that
u ψ̇ + v̇ + βv = u · r − ur − βv + βv = 0
(2.82)
u̇ − vψ̇ + βu = vr + τu − vr + βu = τu + βu
one obtains
ÿ+β ẏ cos(ψ)0+sin(ψ)(τu +βu)
ẍ+β ẋ = cos(ψ)(τu +βu)−sin(ψ)0 →
(2.83)
ÿ+β ẏ
ẍ+β ẋ = tan(ψ) → ψ = atan −1 ( ẍ+β
ÿ+β ẏ
ẋ )
Thus, through Eq. (2.83) it is proven that the state variable ψ (heading angle of the
vessel) is a function of the flat output and of its derivatives.
From Eq. (2.81) one also has that
while using Eq. (2.84) and after intermediate computations one finally obtains
Dividing Eq. (2.86) with the square root of Eq. (2.84) one obtains
u(τu +βu)
√ẋ(ẍ+β ẋ)+ ẏ( ÿ+β ẏ) = (τu +βu) (2.87)
(ẍ+β ẋ) +( ÿ+β ẏ)
2 2
ẏ ẍ − ẋ ÿ
v= (2.91)
(ẍ + β ẋ)2 + ÿ + β ẏ)2
r = ψ̇ (2.92)
ÿ + β ẏ
ψ = atan −1 ( ) (2.93)
ẍ + β ẋ
which means that r is also a function of the flat output and of its derivatives. This
can be also confirmed analytically. Indeed from Eq. (2.93) it holds that
cos 2 (ψ)ψ̇ + sin 2 (ψ)ψ̇ (y (3) + β ψ̈)(ẍ + β ẋ) − ( ÿ + β ẏ)(x (3) + β ẍ)
=
cos 2 (ψ) (ẍ + β ẋ)2
(2.94)
which also gives
Finally, for the control input τr it holds that τr = ṙ and using Eq. (2.99) this implies
that τr is also a function of the flat output and of its derivatives. This can be also
shown analytically according to the following:
τr = ṙ ⇒ τr =
y (4) (ẍ+βx)−x (4) ( ÿ+β ẏ)+β(y (3) ẍ−x (3) ÿ)−β 2 (x (3) ẏ−y (3) ẋ)
[(ẍ+β ẋ)2 +( ÿ+β ẏ)2 ]· (2.102)
(3)
−2 [y (ẍ+β ẋ)−x 3 ( ÿ+β ẏ)−β 2 (ẍ ẏ− ÿ ẋ)]
[(ẍ+β ẋ)2 +( ÿ+β ẏ)2 ]2
·
·{(ẍ + β ẋ)(x + β ẍ) + ( ÿ + β ẏ)(y (3) + β ÿ)}
(3)
Through Eq. (2.102) it is confirmed that that all state variables and the control input
of the hovercraft’s model can be written as functions of the flat output and of its
derivatives. Consequently, the vessel’s model is a differential flat one.
Next, it will be shown that a flatness-based controller can be developed for the
hovercraft’s model. It has been shown that it holds
By differentiating once more with respect to time and after intermediate operations
one finally obtains
x (3) = τ̇u cos(ψ) − τu sin(ψ)r +
(2.103)
+β(−ur − βv)sin(ψ) + βvcos(ψ)r
As explained in Eq. (2.104), the state vector of the system is extended so as to include
as new state variables the control input τu and its first derivative τ̇u . The new state
variables are denoted as z 1 = τu and z˙1 = τ̇u . Using that the extended state-space
description of the system is
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
ẋ ucos(ψ) − vsin(ψ) 0 0
⎜ ẏ ⎟ ⎜ usin(ψ) + vcos(ψ) ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ ψ̇ ⎟ ⎜ r ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ u̇ ⎟ ⎜ vr + z 1 ⎟ ⎜0 0⎟
⎜ ⎟=⎜ ⎟+⎜ ⎟ τ̈u (2.104)
⎜ v̇ ⎟ ⎜ −ur − βv ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟ τr
⎜ ṙ ⎟ ⎜ 0 ⎟ ⎜0 1⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ż 1 ⎠ ⎝ z2 ⎠ ⎝0 0⎠
ż 2 0 1 0
It can be noticed that in the equations of the third-order derivatives for both x and
y only the control input τu and its derivative τ̇u appear, while the control input τr
is missing. Therefore, differentiation of x (3) once more with respect to time is per-
formed. This gives
68 2 Differential Flatness Theory and Flatness-Based Control
while after substituting the time derivative according to Eq. (2.74) and after regroup-
ing terms one obtains a description of the form
where
L ga L 3f y1 = cos(ψ) (2.110)
In a similar manner, differentiating once more with respect to time the expression
about y (3) one gets
while after substituting the time derivative according to Eq. (6.166) and after regroup-
ing terms one obtains a description of the form
(2.113)
2.3 Properties of Differentially Flat Systems 69
v1 = L 4f y1 + L ga L 3f y1 τ̈u + L gb L 3f y1 τr
(2.119)
v2 = L 4f y2 + L ga L 3f y2 τ̈u + L gb L 3f y2 τr
the following description for the input–output linearized hovercraft model is obtained
x (4) = v1
(2.120)
y (4) = v2
For the dynamics of the linearized equivalent model of the system, the following
new state variables can be defined
ż = Az + Bv
(2.122)
zm = C z
70 2 Differential Flatness Theory and Flatness-Based Control
or equivalently
⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
ż 1,1 0 1 0 0 0 0 0 0 z 1,1 0 0
⎜ż 1,2 ⎟ ⎜0 0 1 0 0 0 0 0⎟ ⎜z 1,2 ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟
⎜ż 1,3 ⎟ ⎜0 0 0 1 0 0 0 0⎟ ⎜ ⎟ ⎜ 0⎟
⎜ ⎟ ⎜ ⎟ ⎜z 1,3 ⎟ ⎜0 ⎟
⎜ż 1,4 ⎟ ⎜0 0 0 0 0 0 0 ⎟ ⎜
0⎟ ⎜z 1,4 ⎟ ⎜ 0⎟
⎜ ⎟ ⎜ ⎟ + ⎜1 ⎟ v1
⎜ż 2,1 ⎟ = ⎜0 0 0 0 0 1 0 ⎟ ⎜
0⎟ ⎜z 2,1 ⎟ ⎜
⎟ 0⎟
(2.123)
⎜ ⎟ ⎜ ⎜0 ⎟ v2
⎜ż 2,2 ⎟ ⎜0 0 0 0 0 0 1 0⎟ ⎜z 2,2 ⎟ ⎜0 0⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟
⎝ż 2,3 ⎠ ⎝0 0 0 0 0 0 0 1⎠ ⎝z 2,3 ⎠ ⎝0 0⎠
ż 2,4 0 0 0 0 0 0 0 0 z 2,4 0 1
Φμ : X × U × (R m+k )μ → Y × V × (R s )μ
(2.125)
(x, u, u 1 , . . . , u k+μ ) → (φ, α, α̇, . . . , α (μ) )
is as follows: Variable (X, τx ) comprising the state vector of the nonlinear system
and its successive derivatives with respect to time, as well as the Lie derivatives of
the mapping transformation along a vector field defined by the state vector elements
and their derivatives is called a manifold of jets or a diffiety.
Given two diffieties (M, CTM) and (N , CTN) it is said that a function Ψ from
M to N is Lie-Backlünd or a C-morphism, if the tangent (projection) function T Ψ
satisfies T Φ(CTM) ⊂ CTN. Moreover, if Ψ has an inverse function Φ, such that
T Ψ (C T N ) ⊂ CTM, then Ψ is a Lie-Backlünd isomorphism or a C isomorphism.
When such an isomorphism exists, the diffieties are called equivalent. Therefore,
an endogenous transformation is a special type of a Lie-Bakclünd transformation
which preserves the parametrization in time of the integral curves. It is also possible
to introduce a more general concept which is the orbital equivalence, about the Lie-
Bakclünd isomorphisms which preserve only the geometric locus of the integral
curves.
Definition: The dynamical system (M, F) where M: is the state vector and F: is a
vector field, is differentially flat in p (or respectively flat), if and only if it is equivalent
in p (respectively equivalent) with the trivial system (therefore it can be described
in the Brunovsky canonical form, as shown in the examples of Sect. 2.2.2).
The definition of equivalent systems can be used, and one can consider a differ-
entially flat system (X × U × Rm ∞ , F) associated with
By definition, the above system is equivalent to the trivial system Rn∞ , Fs , where the
endogenous transformation Ψ takes the form
In other terms, Ψ is the infinite extension of function h(). The inverse transformation
of Ψ is denoted as Φ and holds
Definition: Assume that {M, F} is a system and that Ψ is the transformation that
brings it to the canonical form. The first component of Ψ is considered to be the flat
output.
Another important property of differentially flat systems is that the dimension of
the flat output is equal to the number of the control inputs (and thus linearization
and decoupling is succeeded). However, it is noted that the flat output is not unique.
If y is a flat output, then γ (y) is also a flat output, provided that γ is a smooth
diffeomorphism. For single-input systems, the flat output is a scalar s = m = 1.
It is easy to show that one can pass from one flat output to another by using only
a static change of variables. The multivariable case is clearly more complicated.
For example, for s = m = 2 it is easy to see that if y1 , y2 is a flat output then
(y1 , y2 + γ (y1(r ) )) is also a flat output, provided that γ () is smooth function and
r ≥ 0. Since the flat output of a system is not unique, it is preferable to select the flat
output which leads to simple computation for the linearization of the system and the
design of the feedback controller.
Next, it will be shown how differential flatness can be used to solve the problem
of trajectory planning (setpoints definition) for dynamical systems control. Assume
the following nonlinear system in which the control input u does not implement
feedback control
ẋ = f (x, u), x ∈ R n u ∈ R m (2.132)
with flat output y = h(x, u, u̇, . . . , u (r ) ). In such a manner, the system’s trajec-
tories (x(t), u(t)) are written as a function of the flat output y and its derivatives
x = φ(y, ẏ, . . . , y (q) ) and u = α(y, ẏ, . . . , y (q) ). The trajectory planning problem
contains a start and an arrival state. It is assumed that the flat output is a vector with
elements yi which are written as
yi (t) = Ai j λ j (t) (2.133)
j
Next, it is assumed that the initial state of the system is x0 at t0 , while the final
state of the system is x f at t f . The coefficients Ai j should satisfy the following
yi (t0 ) = j Ai j λ j (t0 ) yi (t f ) = j Ai j λ j (t f )
··· ··· ··· (2.134)
(q) (q) (q) (q)
yi (t0 ) = j Ai j λ j (t0 ) yi (t f ) = j Ai j λ j (t f )
Without loss of generality it assumed that the vector of the flat output is of dimension
i = 1, that is y = [y1 ] (these results can be also generalized to the multidimensional
case). The following vectors are defined
(q)
ỹ0 = (y1 (τ0 ), . . . , y1 (τ0 ))
(q)
ỹ f = (y1 (τ f ), . . . , y1 (τ f )) (2.135)
ỹ = ( ỹ0 , ỹ f )
The elements of A should satisfy Eq. (2.136). The only condition for Eq. (2.136) to
have a solution, is that matrix Λ is a full rank one. This means that the space of
functions λ j has to be sufficiently rich. It can be concluded that the path planning
problem for the case of differentially flat systems ends up to linear algebra theory.
In this case, the objective is to find a trajectory going from point a to point b, which
will satisfy the constraints K (x, u, . . . , u (r ) ) ≤ 0 at each time instant. In the flat
coordinates, this consists in finding T > 0 and a function [0, T ]
t → y(t) with
(y, . . . , y (q) ) to be given at t = 0 and T and to verify that ∀t ∈ [0, T ] the constraint
K (y, . . . , y (v) )(t) ≤ 0 will hold, for a certain v and a specific function K , obtained
from k, φ, and α, where k is the constraint condition and φ, α are the previously
defined mappings ẋ = φ(y, v̄) and u = α(y, v̄). The difficulty of this problem
increases when q = v = 0.
Next, it is simply assumed that the initial state Y0 and the final state Y f are equi-
libria. It is also assumed that the motion from the initial to the final state satisfies the
following constraint: there exists a path [0, 1]
σ → Y (σ ) such that Y (0) and Y (1)
correspond to start and arrival equilibria , for all σ ∈ [0, 1], K (Y (σ ), 0, . . . , 0) < 0.
Therefore, there exists T > 0 and [0, T ]
t → y(t) which is a solution to the
problem. It suffices to take Y (η(t/T )) where T is sufficiently large, and with η to be
a smooth increasing function [0, 1]
s → η(s) ∈ [0, 1] with η(0) = 0, η(1) = 1,
and ddsηi = 0 for i = 1, . . . , max(q, v).
i
74 2 Differential Flatness Theory and Flatness-Based Control
The problem becomes more complicated in the case of trajectory planning for
systems subjected to nonholonomic constraints.
With the previous analysis it has been shown that the endogenous transformation
associated with the flat output y = h(x, ū) is defined everywhere, is smooth and
invertible in a manner that always enables state variable x and control input u to be
expressed as functions of the flat output and its derivatives
However, a point of singularity may exist in the area of the state-space into which
the control tries to bring the system’s state vector (y is no longer invertible and
Eq. (2.138)) cannot be solved. As function φ cannot be defined in such a point, the
previous computations become meaningless. A manner to circumvent this problem
and to avoid the singularity is to choose a trajectory for the system
such that the point of potential singularity is bypassed. This procedure is depicted in
the following example:
Example: The following system with differentially flat dynamics is considered:
where the flat output is defined as y = (x1 , x3 ). In case that u 1 = 0 it holds that
ẋ1 = ẏ1 = 0 and a singularity arises because it holds
φ ẏ2 ÿ2 ẏ1 − ÿ1 ẏ2
(y, ẏ, ÿ) (x , x , x , u , u ) = y1 , ẏ1 , y2 , ẏ1 , (2.141)
→ 1 2 3 1 2 ẏ 3
1
and the inverse transformation φ is not defined at the points of singularity. However,
if one considers trajectories of the following form t → y(t) := (σ (t), p(σ (t))),
with σ and π being smooth functions he obtains that
dp
ẏ2 (t) dσ (σ (t)) · σ̇ (t)
= (2.143)
ẏ1 (t) σ̇ (t)
2.3 Properties of Differentially Flat Systems 75
and
d2 p
ÿ2 ẏ1 − ÿ1 ẏ2 dσ 2
(σ (t)) · σ̇ 3 (t)
= (2.144)
ẏ13 σ̇ 3 (t)
Therefore, one can extend t → φ(y(t), ẏ(t), ÿ(t)) everywhere in time using
dp d2 p
t → σ (t), (σ (t)), p(σ (t)), σ̇ (t), (σ (t)) (2.145)
dσ dσ 2
Considering that these systems are equivalent, they have the same trajectories. A
question that arises is if a transition is possible from the system ẋ = f (x, u) to the
76 2 Differential Flatness Theory and Flatness-Based Control
ż = α(x, z, v) z ∈ Z ⊂ R q
(2.149)
u = κ(x, z, v)
such that
F(Φ(y, v̄)) = DΦ(y, v̄) · G(y, v̄) (2.153)
Next, one sets ȳ = (y, v, v(1) , . . . , v(μ) ) and w = v(μ+1) . For μ sufficiently large,
it holds that φ (respectively α) depends exclusively on ȳ (respectively on ȳ and w).
With these notations Φ is written as
where ḡ := (g, v(1) , . . . , v(k) ). Since Φ is invertible, and φ is of full rank it can thus
be completed by a function π to obtain a change of coordinates
u = κ(x, z, w)
(2.160)
ż = α(x, z, w)
The term endogenous has the meaning that variables z and w are not included in the
state vector of the system (case also known as “zero dynamics”). Using the above one
arrives at a more refined definition of equivalence between systems and differential
flatness.
78 2 Differential Flatness Theory and Flatness-Based Control
Theorem: Two dynamics ẋ = f (x, u) and ẏ = g(y, v) are equivalent, if and only if,
ẋ = f (x, u) can be transformed by endogenous feedback and change of coordinates.
u̇ = κ(x, z, w)
(2.163)
ż = α(x, z, w)
for μ sufficiently large (thus finally the system’s control input contains multiple inte-
gral terms). This shows clearly the properties which are preserved by equivalence.
The property which is maintained by introducing integral terms and change of coor-
dinates is controllability. An endogenous feedback is actually reversible to the effect
of multiple integrators.
The problem of trajectory tracking for the system ẋ = f (x, u) consists in finding a
controller u which will permit to the system state variables to track the reference set-
point t → (xr (t), u r (t)). It is convenient to add to the open loop control a correction
term Δu, which is a function of the state vector error Δx = x − xr . For differen-
tially flat systems, there is a systematic method to calculate Δu from the tracking
error Δx.
If the dynamics admits as flat output y = h(x, ū), one can use the corollary
that a differentially flat system is linearizable by (dynamic) feedback and change of
coordinates, into the linear canonical form y (μ+1) = w. For the linearized equivalent
(μ+1)
model, the control input is defined as yr − K Δ ỹ, where K is a gains vector and
Δy is a vector having as elements the flat output tracking error and its derivatives
up to order μ. By applying such a feedback control input, the tracking error dynamics
becomes
Δy (μ+1) = −K Δ ỹ (2.165)
2.3 Properties of Differentially Flat Systems 79
where yr (t) := (xr (t), ū r (t)) and ỹ := (y, ẏ, . . . , y μ ), and Δy represents the vector
of the flat output tracking error. Such a feedback control law assures asymptotic
tracking. There exists an invertible transformation
which connects the infinite vector fields F(x, ū) := ( f (x, u), u, u (1) , . . .) and
G( ȳ) := (y, y (1) , . . .). From the proof of the theorem on the equivalence of two
systems (theorem given in Sect. 2.3.3.1), the above means that for the state vector
one has
x = φ( ỹr (t) + Δ ỹ) =
= φ( ỹr (t)) + Rφ (yr (t), Δ ỹ)Δ ỹ = (2.167)
= xr (t) + Rφ (yr (t), Δ ỹ)Δ ỹ
Up to now, for the systems used for analyzing differential flatness properties and
flatness-based control it has been assumed that an exact dynamic model was avail-
able. Flatness-based control can be also applied to systems characterized by model
uncertainties and exogenous disturbances. Denoting the coefficients vector of the
system as θ ∈ R p , the generalized system dynamics of Eq. (2.170) is written as
In the generic case, the system’s coefficients θ are subject to uncertainty, i.e.,
where θ0 is the nominal value of the parameter vector. In Eq. (2.170) the vector
field f : R p × R n × R m → R n is smooth. Denoting the flat output y(t) ∈ R
and parameters vector θ , the definition of differential flatness given in Sect. 2.2.2 is
rewritten as [189, 190]
y = h(θ, x, u, u̇, . . . , u (r ) )
x = φ(θ, y, ẏ, . . . , y (r −1) ) (2.172)
u = ψ(θ, y, ẏ, . . . , y (r ) )
which means that for each admissible nominal trajectory y ∗ (t) there is a nominal
control input u ∗ . The initial condition of the desired trajectory of the flat output
t → y ∗ (t) is defined as the vector y0∗ = [y ∗ (0), ẏ ∗ (0), . . . , y ∗ (r −1) (0)]T . The sys-
tem is consistent with respect to the initial conditions x0 if x0 = φ(θ, y0∗ ). The
following theorem has been proven [189]:
Theorem: If the desired trajectory of the flat output is consistent with the initial con-
dition x0 and θ = θ0 , then when applying the nominal control input of Eq. (2.173)
the system becomes equivalent, by a change of coordinates, to a linear system in
Brunovksy (canonical) form.
where h = φ −1 with φ defined in Eq. (2.172), into the control normal (canonical)
form
ẏi (t) = yi+1 (t) i ∈ {1, . . . , n − 1}
(2.175)
ẏn (t) = α(θ, y(t), u(t))
where y (i) is the i − 1th order derivative of the flat output. From the second row of
Eq. (2.175), one can see that the nth order derivative of the flat output is finally written
as a nonlinear function of the parameters vector θ , of the flat output’s vector y(t) and
of the control input u(t). The system in canonical form described by Eq. (2.175) and
the system of Eqs. (2.170), (2.172), and (2.173) have the same solution t → y(t).
Defining the new control input v(t) = α(θ, y(t), u(t)), the Brunovsky form of
the transformed initial system is written as
Once the system has been written in the Brunovksy (canonical) form of Eq. (2.176),
a control input that can assure tracking of any desirable trajectory ẏn∗ is given by
v = ẏn∗ + K T ē (2.177)
Thus, in case that the exact values of the model’s parameters θ are known and these
are equal to the nominal values θ0 , and no integration of the flat output error is
nin the feedback loop, i.e. the feedback control signal is of the form K e =
used T
− i=1 kn+1−i ei (t), the closed-loop system error vector dynamics become
∗
where ei = yi − yi∗ = y (i−1) − y (i−1) . The dynamics of the flat output tracking
error can be also written in a matrix form
⎛ ⎞ ⎛ ⎞⎛ ⎞
ė1 0 1 0 ··· 0 e1
⎜ ė2 ⎟ ⎜ 0 0 1 ··· 0 ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟ ⎜ e2 ⎟
⎜· · ·⎟ = ⎜ · · · · · · ··· ··· ··· ⎟ ⎜ ⎟
⎜ ⎟ ⎜ ⎟ ⎜· · ·⎟ (2.180)
⎝· · ·⎠ ⎝ 0 0 0 · · · 1 ⎠ ⎝· · ·⎠
ėn −kn −kn−1 −kn−2 · · · −k1 en
In case of the error dynamics of Eq. (2.180) to succeed elimination of the tracking
error it suffices to select the controller gains k1 , k2 , . . . , kn so as the polynomial
e(n) + k1 e(n−1) + k2 e(n−2) + · · · + kn−1 ė + kn e = 0 to be a Hurwitz one.
When that the model’s parameters θ are not precisely known, or have deviated
from their nominal values θ0 , one can apply again a control input of the form described
in Eqs. (2.177) and (2.178), and for specific magnitude of the exogenous disturbances
or parametric uncertainties can expect again satisfactory tracking of the reference
trajectory. In such a case, the robustness features of the closed loop can be analyzed
with the use of results from interval polynomial stability theory (e.g. Kharitonov’s
theory) or other results on interval polynomial analysis (e.g. see [431]). The problem
of robustness of a flatness-based controller in case of parametric uncertainties or
exogenous perturbations has been studied in [189, 190, 433, 443].
the control input. Equivalently, in a 0-flat system its flat output depends exclusively
on the elements of the system’s state vector.
The next example shows that such an upper bound should depend, in the best
case, linearly to the dimension of the state vector. Consider the dynamical system
(a1 ) (a2 )
x1 = u 1 x1 = u 2 ẋ3 = u 1 u 2 (2.181)
It has been shown that dynamical systems which can be linearized with static state
feedback can be also transformed to the linear canonical (Brunovsky) form. Thus,
such systems are differentially flat, with the flat outputs to be exclusively functions
of the state vector x (0-flat system). It is noted that a differentially flat system is
not always a system which can be linearized with static state feedback. There exist
systems which can be linearized through dynamic feedback linearization, that is their
state vector is extended by including as additional state variables the control inputs
and their derivatives (see previous examples on the linearization of the VTOL and
the hovercraft model, given in Sect. 2.2.2). Such systems can be also written in the
linear canonical form, and thus are also differentially flat.
ẋ = f (x, u) x ∈ R n , u ∈ R (2.183)
can be linearized with static state feedback linearization and can be finally written
in the linear canonical form. Such systems are differentially flat.
that is systems in which the number of control inputs is by one less than the dimen-
sion of the system’s state vector, are 0-flat, as long as they are controllable (strongly
accessible in x). The problem becomes more complicated when the system is not an
affine-in-the-input one. However, the latter systems can be finally transformed into
affine-in-the-input ones too.
are flat if and only if, the generic rank of E k is equal to k + 2, for k = 0, . . . , n − 2,
where E 0 = { f 1 , f 2 } and E k+1 = {E k , [E k , E k ]}, k > 0. A dynamical system
without drift and with two control inputs is always 0-flat. Such a system can be
written in the chained form (which is a precursor of the linear canonical form),
through static state feedback and a change of variables. For example, for a system
of the form
ẋ1 = u 1 ẋ2 = u 2 ẋ3 = x2 u 1 · · · ẋn = xn−1 u 1 (2.186)
is a differentially flat one if it is controllable (strongly accessible for almost all x).
More precisely, such a system is 0-flat is n is an odd number and is a 1-flat one if n
is an even number.
The following theorem provides a sufficient condition for showing that dynamical
systems are not differentially flat [340]:
Theorem: The criterion of the regulated variety. Assume that the dynamical system
ẋ = f (x, u) is a differentially flat one. Then, the projection to the p-dimensional
space of its state-space equation p = f (x, u) is a subvariety regulated for all x.
This criterion signifies that the elimination of the control input u from the n
equations ẋ = f (x, u) leads to n − m equations of the form F(x, ẋ) = 0 with the
following property: For all (x, p) such that F(x, p) = 0 there exists α ∈ R n , α = 0
such that
∀λ ∈ R F(x, p + λa) = 0 (2.189)
Proof : Consider (x̄, ū, p̄) such that p̄ = f (x̄, ū). Generally, f is of the order m =
dim(u) with respect to u. If x is removed from the state-space description then one
arrives at n−m equations which are based exclusively on x and p, that is F(x, p) = 0.
Thus one has the constraint about the trajectory around (x̄, ū) : F(x, ẋ) = 0. If the
system is differentially flat it holds that
F(x̄, x̄ + φ y α ξ ) = 0 (2.192)
is not a differentially flat one. This is because the subvariety p3 = p12 + p23 is not a
regulated one. There does not exist an α ∈ R 3 , α = 0 such that
Actually, the cubic term implies that α2 = 0 and the quadratic term implies that
α1 = 0 and thus one has α3 = 0.
There are some nonflat systems which do not satisfy differential flatness properties.
Such systems cannot be linearized by state feedback (neither by static feedback
linearization nor by dynamic feedback linearization). For such cases, the notion
of Liouvillian systems has been introduced, which can be seen as an extension of
differentially flat systems.
Nonflat systems are algebraically characterized by an integer called the defect.
This integer is well defined by introducing the notion of maximal linearizing output.
Using an informal definition, a set of variables y = (y1 , . . . , ym ) is a maximal
linearizing output if the number of system variables which can be expressed as a
differential function of y is maximal. Therefore, the maximal linearizing output
characterizes the largest subsystem which satisfies the differential flatness property.
This largest subsystem is called a flat subsystem.
2.5 Classification of Types of Differentially Flat Systems 87
The defect represents the number of the system’s variables which do not belong
to the flat subsystem. A nonlinear system is said to be Liouvillian or integrable by
quadratures if the elements which do not belong to the flat subsystem can be obtained
by elementary integrations of the elements of the flat subsystem. In the following
examples, the notion of integrability by quadratures and that of Liouvillian systems
is explained [340]
Example 1:
This system is not linearizable by static or dynamic state feedback. The system is of
defect equal to 1 and is found to be Liouvillian (there is a state variable which does
not belong to the flat subsystem). A maximum linearizing output is given by y = x2
and x1 is the only variable that does not belong
to the flat subsystem. However, x1
can be obtained from the integration x1 = y + ẏ 2 , that is x1 is given by a pure
integral of a differential function of y.
This system is nilpotent and completely controllable since it is driftless and sat-
isfies the Lie algebra rank condition. However, the system is not flat and can-
not be transformed into the canonical form. Nevertheless, there is a subsystem
within it defined by the first four equations, which is differentially flat and con-
sequently the system’s defect is equal to 1. The maximal linearizing output of
the system is given by y = (y1 , y2 ), where y1 = x1 and y2 = (zx1 − x4 )/3,
with z to be defined by z = (x1 x2 − x3 )/2. Indeed z = ẏ2 / ẏ1 , thus x4 =
y1 z − 3y2 = α(y1 , ẏ1 , y2 , ẏ2 ), x3 = ẋ4 / ẏ1 = b(y1 , ẏ1 , ÿ1 , y2 , ẏ2 , ÿ2 ), x2 =
(2z + x3 )/y1 = c(y1 , ẏ1 , ÿ1 , y2 , ẏ2 , ÿ2 ), u 2 = ẋ2 = q(y1 , . . . , y1(3) , y2 , . . . , y2(3) ).
(3)
Since x5 = x3 u 2 = b(y1 , . . . , ÿ2 )q(y1 , . . . , y2 ) one gets that the above system
of Eq. (2.196) is a Liouvillian one.
are elaborated criteria that enable to determine if a system is differential flat [286].
These state that flatness is equivalent to the property of strong closedness of the ideal
1-form representing the differentials of all possible trivializations.
For a nonlinear dynamical system ẋ = f (x, u) (that is a system in the explicit
form) which is described by a set of ordinary differential equations, its implicit
description is obtained after eliminating the control input from the state-space model.
This takes the form:
F(x, ẋ) = 0 (2.197)
given by
while also the Lie derivative of function φ ∈ C ∞ (X, R) along vector field τX is
introduced
dφ n ( j+1) ∂φ
= L τX φ = x ( j)
(2.200)
dt i=1 j≥0 i ∂x i
( j) ( j) ( j) ( j+1)
Since it holds that dtd
xi = L τX xi = ẋi = xi , the Cartan vector field acts
on coordinates as a shift to the right. Variable (X, τX ), comprising the state vector
of the nonlinear system and its successive derivatives with respect to time, as well
as the Lie derivatives of the mapping transformation function φ along a vector field
defined by the state vector elements and their derivatives is called manifold of jets of
infinite order or diffiety.
Thus, in regular implicit control systems, one has no longer the description of the
system in the state-space form but uses a description based on the state vector and
its derivatives, as well as the derivatives of the mapping φ. Next, the definition of a
regular implicit control system is given:
Definition: A regular implicit control system is defined as a triple (X, τX , F) with
X = X × R∞ n , τ is the trivial Cartan field on X and F ∈ C ∞ (TX; R n−m ), satisfying
X
∂F
rank( ∂ x = n − m) in a suitable open dense subset of TX.
2.6 Elaborated Criteria for Checking Differential Flatness 89
The notion of the Lie-Backlünd equivalence for explicit control systems is now
extended to implicit control systems (systems described by manifolds of jets or
diffieties). Two regular implicit control systems are defined [286].
Definition 1: Let X, τX , F and Y, τY , G be two regular implicit systems. These are
Lie-Backlünd equivalent at (x̄0 , ȳ0 ) ∈ X0 × Y0 if and only if
(i) There exist neighborhoods X0 of x̄0 ∈ X0 and ȳ0 ∈ Y0 and a one-to-one map-
ping Φ = (φ0 , φ1 , . . .) ∈ C ∞ (Y0 ; X 0 ) with C ∞ (Y0 ; X 0 ) inverse Ψ , satisfying
Φ( ȳ0 ) = x̄0 and such that the restrictions of the trivial Cartan fields τY|y0 are
Φ-related, namely Φ ∗ τY|y0 = τX|x0
(ii) The C ∞ (Y0 |X0 ) inverse mapping Ψ = (ψ0 , ψ1 , . . .) is such that Ψ (x̄0 ) = ȳ0
and Ψ ∗ τX|x0 = τY|y0 .
The mappings Φ and Ψ are called mutually inverse Lie-Backlünd isomorphisms at
(x0 , y0 ).
It holds that two explicit control systems are Lie-Backlünd equivalent if and
only if the associated implicit control systems are Lie-Backlünd equivalent. This is
analytically stated as follows:
Two implicit systems (X, τX , F) and (Y, τY , G) are given and two vector fields f
and g compatible with (X, τX , F) and Y, τY , G, respectively are defined. The corre-
sponding explicit systems ẋ = f (x, u) and ẏ = g(y, v) are differentially equivalent
or Lie-Backlünd equivalent at the pair (x0 , u 0 , u̇ 0 ), . . . , and (y0 , v0 , v̇0 ), . . . (with u 0
such that ẋ0 = f (x0 , u 0 ) and v0 such that ẏ0 = g(x0 , v0 )), if and only if (X, τX , F)
and (Y, τY , G) are Lie-Backlünd equivalent in accordance to the previous Definition
1, at (x0 , ẋ0 , . . .) and (y0 , ẏ0 , . . .).
Next, it is pointed out once more that an explicit system is flat, if and only if, it
is equivalent to the trivial system (linear canonical Brunovsky form). This enables
to define also necessary and sufficient condition for implicit systems to be differen-
tially flat.
Definition 2: The implicit system (X, τX , F) is flat at x̄0 if and only if there exists
y0 ∈ R∞ m , such that (X, τ , F) is Lie-Backlünd equivalent at ( x̄ , ȳ ) ∈ X × R n
X 0 0 0 ∞
to the m-dimensional trivial implicit system (R∞m , τ , 0). In this case, the mutually
m
inverse Lie-Backlünd isomorphisms Φ and Ψ are called inverse trivializations.
Having denoted the implicit system as F(x, ẋ), its differential can be defined next
n
∂ Fi ∂ Fi
dF = dx j + d ẋ j i = 1, . . . , n − m (2.201)
j=1 ∂x j ∂ ẋ j
j ( j)
where it holds that φi (ψ̄) = xi . Then the following theorem holds about the
differential flatness of the implicit system
Theorem: The implicit system (X, τX , F) is flat at x̄0 , with x̄0 ∈ X0 if and only if
there exists ȳ0 ∈ R∞ m and a local Lie-Backlünd isomorphism Φ from R m to X
∞ 0
satisfying Φ( ȳ0 ) = x̄0 and such that Φ ∗ · dF = 0.
∂F ∂F d or d(φ 0 ) ∂φ 0 d j
P(F) = ∂x + ∂ ẋ dt P(φ 0 ) = j=0 ∂ y ( j) dt j
(2.203)
p , 0 p,q−p ) if p ≤ q
VMU = (Δ
Δq (2.204)
VMU = if p > q
0 p−q,q
Other definitions associated with the Smith decomposition of the differential matrix
P(F) of the implicit system F(x, ẋ) are as follows:
It holds that VP(F)U = (Δ, 0n−m,m ) and matrix U is denoted as R-Smith(P(F)).
Next one denotes
0
Û = U n−m,m (2.207)
Im
Next, it is considered that U is the right matrix in the Smith decomposition of the
implicit system’s differential matrix P(F), that is U ∈ R − Smith(P(F)). Then the
following m-dimensional vector 1-form ω is defined
⎛ ⎞ ⎛ ⎞
n j α
ω1 (x̄) j=1 α≥0 Q 1,a ( x̄)d x j |X0
⎜ ⎟
ω(x̄) = ⎝ · · · ⎠ = Q̃ 0 d x|X0 = ⎝ ··· ⎠ (2.211)
n
ωm (x̄) j=1 α≥0 Q
j
m,a ( x̄)d x α|
X0 j
Using the above one arrives at a criterion for defining flat outputs.
Definition: Variational flat outputs. The vector 1-form ω = Q̃ 0 d x|X0 , defined in
Eq. (2.211), using matrices Q 0 and Q̃ 0 defined in Eqs. (2.229) and (2.210), respec-
tively, is a flat output of the variational system.
The following relations are defined next:
Definition: It is said that the [ dt
d
]-ideal T generated by τ1 , . . . , τΓ is strongly closed,
if there exists a matrix M ∈ Um [ dτ d
] such that d(Mτ ) = 0.
Theorem: A necessary and sufficient condition for systems X, τX , F to be flat at
(x̄0 , ȳ0 ) (over the class of meromorphic functions) is that there exist U ∈ R −
Smith(P(F)) and Q ∈ L = Smith(Û ), with Û to be given by
0n−m.m
Û = U (2.212)
Im
92 2 Differential Flatness Theory and Flatness-Based Control
The condition given in Eq. (2.214) may be seen as a generalization in the framework
of manifolds of jets of infinite order of the well-known moving frame structure
equations. Furthermore, the general matrix solution μ = (μik )i,k=1,...,m of dω = μω,
with ω defined in Eq. (2.211) is given by
the integer ord(μ) being arbitrary but otherwise finite and satisfying or d(μ) ≥
j,k
or d(Γ ), the Γi,α,β being given by the following relation
j,k
and where vi,α,β ’s are meromorphic functions depending on successive derivatives
of x.
Moreover, the following hold:
The [ dt
d
]-ideal Ω generated by the components of the vector 1-form ω defined
in Eq. (2.211) is strongly closed in X0 , or equivalently the system (X, τX , F) is flat,
if and only if, there exists μ ∈ L 1 (Λ(X)m ), and two matrices M ∈ Mm,m [ dt d
] and
N ∈ Mm,m [ dt ], such that
d
The procedure for testing if a system is differentially flat and for computing the
associated flat outputs is:
1. Write the system in the implicit form F(x, ẋ) by eliminating the control input
from its state-space equations.
2. Compute the differential matrix of the system P(F).
3. Perform Smith decomposition of P(F) and retain the R-Smith matrix U .
4. Compute matrix Û from Eq. (2.228).
5. Compute matrices Q 0 and Q̃ 0 from Eqs. (2.210) and (2.231), respectively.
6. Compute the vector 1-form that is described in Eq. (2.211). The elements of the
ω vector are the flat outputs of the system.
If for some or d(μ), the algorithm produces an invertible M, a flat output is obtained
by integration of dy = Mω, which is possible since d(Mω) = 0. In the opposite
case, the system is not differentially flat.
ẋ = vcos(θ )
ẏ = vsin(θ ) (2.219)
θ̇ = vl tan(φ)
The system is written in implicit form F(x, ẋ) = 0. It has n = 3 states and m = 2
inputs, with n − m = 1 and thus it is equivalent to a single implicit equation which
is obtained
F(x, y, θ, ẋ, ẏ, θ̇ ) = ẋsin(θ ) − ẏcos(θ ) = 0 (2.220)
with ⎛ ⎞
001
U0 = ⎝0 1 0⎠ (2.224)
100
U1 = ⎝ 0 1 0 ⎠ (2.225)
0 0 1
finally gives
P(F)U = 1 0 0 (2.226)
and matrices Q̃ 0 , Q̂ 0 as
Q̃ 0 = Im 0m,n−m Q 0 (2.230)
Q̂ 0 = 0n−m In−m Q 0 (2.231)
one has ⎛ ⎞
0 1 0
Q0 = ⎝ 1 0 0⎠ (2.232)
sin(θ) d
L dt − cos(θ)
L dt
d
1
010
Q̃ 0 = (2.233)
100
which comes from using Eq. (2.220) about the implicit system dynamics. Next matrix
Q̃ 0 is defined as
Q̃ 0 = (Im , 0m,n−m )Q 0 (2.235)
The following nonlinear heat diffusion equation is considered, describing the spa-
tiotemporal variations of temperature in a surface
∂φ ∂ 2φ
= K 2 + f (φ) + u(x, t) (2.237)
∂t ∂x
Using the approximation for the partial derivative in the partial differential equation
given in Eq. (12.38) one has
∂φi K 2K K
= φi+1 − φi + φi−1 + f (φi ) + u(xi , t) (2.239)
∂t Δx 2 Δx 2 Δx 2
By considering the associated samples of φ given by φ0 , φ1 , . . . , φ N , φ N +1 one has
∂φ1
∂t = K
φ
Δx 2 2
− 2K
φ
+ Δx
Δx 2 1
K
2 φ0 + f (φ1 ) + u(x 1 , t)
∂φ2
∂t = K
φ
Δx 2 3
− 2K
φ
+ Δx 2 φ1 + f (φ2 ) + u(x2 , t)
Δx 2 2
K
∂φ3
∂t = K
φ
Δx 2 4
− 2K
φ
+ Δx
Δx 2 3
K
2 φ2 + f (φ3 ) + u(x 3 , t)
(2.240)
···
∂φ N −1
∂t = Δx
K
2 φ N − Δx 2 φ N −1 + Δx 2 φ N −2 + f (φ N −1 ) + u(x N −1 , t)
2K K
∂φ N
∂t = Δx 2 φ N +1 − Δx 2 φ N + Δx 2 φ N −1 + f (φ N ) + u(x N , t)
K 2K K
2.7 Distributed Parameter Systems and Their Transformation … 97
ẋ1 = K
x
Δx 2 2
− 2K
x
Δx 2 1
+ ΔxK
2 φ0 + f (x 1 ) + u(x 1 )
ẋ2 = K
x
Δx 2 3
− 2K
x
Δx 2 2
+ Δx 2 x1 + f (x2 ) + u(x2 )
K
ẋ3 = K
x
Δx 2 4
− 2K
x
Δx 2 3
+ ΔxK
2 x 2 + f (x 3 ) + u(x 3 )
(2.242)
···
ẋ N −1 = Δx K
2 x N − Δx 2 x N −1 + Δx 2 x N −2 + f (x N −1 ) + u(x N −1 )
2K K
ẋ N = Δx 2 φ N +1 − Δx 2 x N + Δx
K 2K K
2 x N −1 + f (x N ) + u(x N )
y1,i = xi
(2.243)
y2,i = ẋi
ẏ1,1 = Δx K
2 y1,2 − Δx 2 y1,1 + Δx 2 φ0 + f (y1,1 ) + u(y1,1 )
2K K
ẏ1,2 = Δx K
2 y1,3 − Δx 2 y1,2 + Δx 2 y1,1 + f (y1,2 ) + u(y1,2 )
2K K
···
(2.244)
···
ẏ1,N −1 = Δx 2 y1,N − Δx 2 y1,N −1 + Δx
K 2K K
2 y1,N −2 + f (y1,N −1 )
+u(y1,N −1 )
ẏ1,N = Δx K
2 φ N +1 − Δx 2 y1,N + Δx y1,N −1 + f (y1,N ) + u(y1,N )
2K K
The dynamical system described in Eq. (13.52) is a differentially flat one with flat
output defined as the vector ỹ = [y1,1 , y1,2 , . . . , y1,N ]. Indeed all state variables can
be written as functions of the flat output and its derivatives. Moreover, by defining
the new control inputs
v1 = Δx K
2 φ0 + f (y1,1 ) + u(y1,1 )
v2 = f (y1,2 ) + u(y1,2 )
v3 = f (y1,3 ) + u(y1,3 )
(2.245)
···
v N −1 = f (y1,N −1 ) + u(y1,N −1 )
v N = Δx K
2 φ N +1 + f (y1,N ) + u(y1,N )
98 2 Differential Flatness Theory and Flatness-Based Control
⎛ 2K ⎞
⎛ ⎞ − Δx 2 ΔxK
2 0 0 ··· 0 0 0 0 ⎛ ⎞
ẏ1,1 ⎜ K − 2K K
0 ··· 0 0 0 0 ⎟ y1,1
⎜ ẏ1,2 ⎟ ⎜ Δx 2 Δx 2 Δx 2 ⎟⎜
⎟ ⎜ y1,2 ⎟
⎜ ⎟ ⎜ K
− Δx 2 Δx 2 · · ·
2K K ⎟
⎜ ··· ⎟ = ⎜ 0 0 0 0 0 ⎟⎜
⎟⎜ ··· ⎟ ⎟+
Δx 2
⎜ ⎟ ⎜
⎜ · · · · · · ··· ··· ··· ··· ··· ··· ··· ⎟⎝
⎝ ẏ1,N −1 ⎠ ⎜ ⎟ y1,N −1 ⎠
⎝ 0 0 0 0 ··· 0 K
− Δx
2K K ⎠
ẏ1,N Δx 2 2 Δx 2 y1,N
0 0 0 0 ··· 0 0 K
Δx 2
− Δx
2K
2
⎛ ⎞⎛ ⎞
1 0 0 ··· 0 0 v1
⎜ 0 1 0 · · · 0 0 ⎟ ⎜ v2 ⎟
⎜ ⎟⎜ ⎟
⎜ 0 0 1 · · · 0 0 ⎟ ⎜ v3 ⎟
⎜
+⎜ ⎟ ⎜ ⎟
⎟⎜ ⎟
⎜· · · · · · · · · · · · ⎟ ⎜ · · · ⎟
⎝ 0 0 0 · · · 1 0 ⎠ ⎝v N −1 ⎠
0 0 0 ··· 0 1 vN
(2.246)
Assuming that all measurements from the set of points x j j ∈ [1, 2, . . . , m] are
available, the associated observation (measurement) equation becomes
⎛ ⎞
⎛ ⎞ y1,1
⎛ ⎞ 1 0 0 ··· 0 0 ⎜
z1 ⎟
⎜ ⎟ ⎜ y1,2 ⎟
⎜ z 2 ⎟ ⎜ 0 1 0 · · · 0 0 ⎟ ⎜ y1,3 ⎟
⎜ ⎟ = ⎜· · · · · · · · · · · ·⎟ ⎜ ⎟ (2.247)
⎝ · · ·⎠ ⎜ ⎟⎜ ⎟
⎝ 0 0 0 ··· 1 0 ⎠⎜ ··· ⎟
zm ⎝ y1,N −1 ⎠
0 0 0 ··· 0 1
y1,N
Thus, in matrix form one has the following state-space description of the system
ỹ˙ = A ỹ + Bv
(2.248)
z̃ = C ỹ
Moreover, denoting a = Dx K
2 and b = − Dx 2 , the initial description of the system
2K
The previously defined state vector of the PDE model is considered without the effect
of the external control input u(x, t).
∂φ ∂ 2φ
= K 2 + f (φ) (2.250)
∂t ∂x
The state vector is again Ỹ = [y1,1 , y1,2 , . . . , y1, i, . . . , y1,N −1 , y1,N ], where y1,1 =
V1 , y1,2 = V2 , . . . , y1,i = Vi , . . . , y1,N −1 = VN −1 and y1,N = VN . It will be shown
in a different manner that the state-space description of the nonlinear PDE dynamics,
which has as control input only the boundary condition φ0 is a differentially flat one.
One has
2K K K
ẏ1,1 = − 2 y1,1 + φ0 + f (y1,1 ) + y1,2 (2.251)
Δx Δx 2 Δx 2
2K K K
ẏ1,2 = − y1,2 + y1,1 + f (y1,2 ) + y1,3 (2.252)
Δx 2 Δx 2 Δx 2
2K K K
ẏ1,3 = − y1,3 + y1,2 + f (y1,3 ) + y1,4 (2.253)
Δx 2 Δx 2 Δx 2
··· ···
2K K K
ẏ1,i = − y1,i + y1,i−1 + f (y1,i ) + y1,i+1 (2.254)
Δx 2 Δx 2 Δx 2
··· ···
2K K K
ẏ1,N −1 = − y1,N −1 + y1,N −2 + f (y1,N −1 ) + y1,N (2.255)
Δx 2 Δx 2 Δx 2
100 2 Differential Flatness Theory and Flatness-Based Control
2K K K
ẏ1,N = − y1,N + y1,N −1 + f (y1,N ) + φ N +1 (2.256)
Δx 2 Δx 2 Δx 2
The flat output is considered to be the state variable y1,N , which is denoted as
y = y1,N . Next, it is shown that all state variables which stand also for virtual
control inputs of the system αi = y1,N −i , can be written as functions of the flat
output y = y1,N
From Eq. (13.62) and considering that φ N +1 is constant one obtains
y1,N −1 = α1 = 1
[ ẏ
K /Δx 2 1,N
+ Δx
2K
2 y1,N − 2Δx 2 φ N +1 − f (y1,N )]
K
(2.257)
⇒y1,N −1 = h 1 (y, ẏ, · · · )
y1,N −2 = α2 = 1
[ ẏ
K /Δx 2 1,N −1
+ 2K
y
Δx 2 1,N −1
− 2Δx K
2 y1,N − f (y1,N −1 )]
(2.258)
⇒y1,N −2 = h 2 (y, ẏ, · · · )
y1,2 = α N −2 = 1
[ ẏ
K /Δx 2 1,3
+ Δx2K
2 y1,3 − Δx 2 y1,4 − f (y1,3 )]
K
(2.260)
⇒y1,2 = h N −2 (y, ẏ, · · · )
y1,1 = α N −1 = 1
[ ẏ
K /Δx 2 1,2
+ Δx2K
2 y1,2 − Δx 2 y1,3 − f (y1,2 )]
K
(2.261)
⇒y1,1 = h N −1 (y, ẏ, · · · )
φ0 = α N = 1
[ ẏ
+ Δx
K /Δx 2 1,1
2K
2 y1,1 − Δx 2 y2 − f (y1,1 )]
K
(2.262)
⇒φ0 = h N (y, ẏ, · · · )
The above procedure confirms that all state variables of the model
and the control input which is the boundary condition φ0 can be written as functions of
the flat output y = y1,N and of the flat output’s derivatives. Consequently, differential
flatness for the heat PDE is shown.
Additionally, one can consider decomposition of the PDE state-space equation
into submodels, where at each submodel the virtual control input is αi = y1,N −i
ẏ1,i = − Δx
2K
2 y1,i +
K
y
Δx 2 1,i−1
+ f (y1,i ) + K
y
Δx 2 1,i+1 (2.264)
··· ···
and with local flat output y1,i one can confirm that all subsystems are differentially
flat.
Chapter 3
Nonlinear Adaptive Control Based
on Differential Flatness Theory
3.1 Introduction
3.2.1 Overview
where f s (x, t), gs (x, t) are nonlinear vector fields defining the system’s dynamics,
u denotes the control input and d̃ denotes additive input disturbances. Knowing
that the system of Eq. (3.1) is differentially flat, the next step is to try to write it
into a Brunovsky form. It has been shown that, in general, transformation into the
Brunovsky (canonical) form can be succeeded for systems that admit static feedback
linearization [340]. Single input differentially flat systems, admit static feedback
linearization, and can be transformed into the Brunovksy form.
The selected flat output is denoted by y. Then, as analyzed in Sect. 3.4.2, for the
state variables xi of the system of Eq. (3.1) it holds that
where v = f (x, t) + g(x, t)(u + d̃) is the control input for the linearized model, and
d̃ denotes additive input disturbances. Thus one can use
where f (x, t), g(x, t) are unknown nonlinear functions, while as mentioned above,
d̃ is an unknown additive disturbance. It is possible to make the system’s state vector
106 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
where ρ is the attenuation level and corresponds to the maximum singular value of
the transfer function G(s) of the linearized model associated to Eqs. (3.4) and (3.5).
The concept of H∞ control comes from the need to compensate for the worst
type of disturbances that may affect a nonlinear dynamical system. This was initially
formulated for linear dynamics described with the use of transfer functions. The
H∞ norm of a linear system with transfer function G(s), is denoted by ||G||∞ and
is defined as ||G||∞ = supω σmax [G( jω)] where σmax is the system’s maximum
singular value [407, 410, 413]. The term sup denotes the supremum or least upper
bound of the function σmax [G( j (ω)], and thus the H∞ norm of G(s) is the maximum
value of σmax [G( j (ω)] over all frequencies ω. H∞ norm has a physically meaningful
interpretation when considering the system y(s) = G(s)u(s). When this system is
driven with a unit sinusoidal input at a specific frequency, σmax |G( jω)| is the largest
possible output for the corresponding sinusoidal input. Thus, the H∞ norm is the
largest possible amplification over all frequencies of a sinusoidal input.
For the measurable state vector x of the system of Eqs. (3.4) and (3.5), and for
uncertain functions f (x, t) and g(x, t) an appropriate control law is
1 (n)
u= [y − fˆ(x, t) − K T e + u c ] (3.7)
ĝ(x, t) d
where the supervisory control term u c aims at the compensation of the approximation
error
as well as of the additive disturbance term d1 = g(x, t)d̃. The above relation
can be written in a state-equation form. The state vector is taken to be e T =
[e, ė, . . . , e(n−1) ], which after some operations yields
and K = [kn , kn−1 , . . . , k2 , k1 ]T . This is again a canonical form for the dynamics of
the tracking error. As explained above, the control signal u c is an auxiliary control
term, used for the compensation of d̃ and w, which can be selected according to H∞
control theory:
u c = − r1 B T Pe (3.12)
The approximation of functions f (x, t) and g(x, t) of Eq. (3.5) can be carried out
with neuro-fuzzy networks (Fig. 3.1). The estimation of f (x, t) and g(x, t) can be
written as [543, 544]:
It is assumed that the weights θ f and θg vary in the bounded areas Mθ f and Mθg
which are defined as
Mθ f = {θ f ∈ R h : ||θ f || ≤ m θ f },
(3.15)
Mθg = {θg ∈ R h : ||θg || ≤ m θg }
108 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
with m θ f and m θg positive constants. The values of θ f and θg that give optimal
approximation are:
where: (i) fˆ(x|θ ∗f ) is the approximation of f for the best estimation θ ∗f of the weights’
vector θ f , (ii) ĝ(x|θg∗ ) is the approximation of g for the best estimation θg∗ of the
weights’ vector θg .
The approximation error w can be decomposed into wa and wb , where
θ̃ f = θ f − θ ∗f
(3.19)
θ̃g = θg − θg∗ .
A difference between the neuro-fuzzy approximator given in Fig. 3.1 and a radial
basis functions (RBF) neural network is the normalization layer that appears between
the Gaussian basis functions layer and the weights output layer. After normalization,
3.3 Flatness-Based Adaptive Fuzzy Control for SISO Dynamical Systems 109
the sum (over the complete set of rules) of the fuzzy membership values of each
input pattern becomes equal to 1.
The adaptation law of the weights θ f and θg as well as of the supervisory control term
u c is derived by the requirement for negative derivative of the quadratic Lyapunov
function
1 T 1 T 1 T
V = e Pe + θ̃ θ̃ f + θ̃ θ̃g (3.20)
2 2γ1 f 2γ2 g
1 T ˙ 1 T ˙
V̇ = 21 ė T Pe + 21 e T P ė + γ1 θ̃ f θ̃ f + γ2 θ̃g θ̃g ⇒ (3.21)
Assumption 1: For given positive definite matrix Q and coefficients r and ρ there
exists a positive definite matrix P, which is the solution of the following matrix
equation:
T 2 1
(A − B K T ) P + P(A − B K T ) − PB( − 2 )B T P + Q = 0 (3.23)
r ρ
V̇ = − 21 e T Qe + 21 e T PB( r2 − 1
ρ2
)B T Pe + B T Pe(− r1 e T PB)+
(3.24)
+B T Pe(w + d1 ) + γ11 θ̃ Tf θ̃˙ f + γ12 θ̃gT θ̃˙g
It holds that
θ̃˙ f = θ̇ f − θ˙∗f = θ̇ f
(3.25)
θ̃˙g = θ̇g − θ˙g∗ = θ̇g
The update of θ f stands for a gradient algorithm on the cost function 21 ( f − fˆ)2 .
The update of θg is also of the gradient type, while u c implicitly tunes the adaptation
gain γ2 . Substituting Eqs. (3.26) and (3.27) into V̇ finally gives
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe + e T PB(w + d1 )−
−e T PB(θ f − θ ∗f )T φ(x) − e T PB(θg − θg∗ )T φ(x)u c ⇒
(3.29)
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe + e T PB(w + d1 ) + e T PBwα .
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe + e T PBw1 or equivalently,
(3.30)
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe + 21 e T PBw1 + 21 w1T B T Pe
Proof : The binomial (ρa − ρ1 b)2 ≥ 0 is considered. Expanding the left part of the
above inequality one gets
ρ2a2 + 1 2
ρ2
b − 2ab ≥ 0 ⇒ 21 ρ 2 a 2 + 2ρ1 2 b2 − ab ≥ 0 ⇒
(3.32)
ab − 1 2
2ρ 2
b ≤ 21 ρ 2 a 2 ⇒ 21 ab + 21 ab − 2ρ1 2 b2 ≤ 21 ρ 2 a 2
The following substitutions are carried out: a = w1 and b = e T PB and the previous
relation becomes
1 T
2 w1 B T Pe + 21 e T PBw1 − 1 T
2ρ 2
e PBB T Pe ≤ 21 ρ 2 w1T w1 (3.33)
The attenuation coefficient ρ can be chosen such that the right part of Eq. (3.34) is
always upper bounded by 0. For instance, it suffices at any iteration of the control
algorithm to have
− 21 e T Qe + 21 ρ 2 ||w||21 ≤0⇒
− 21 ||e||2Q + 21 ρ 2 ||w||21 ≤0⇒ (3.35)
||e||2
ρ 2 ≤ ||w||Q2
1
The performance of the proposed flatness-based adaptive fuzzy controller was first
tested in the benchmark problem of cart-pole balancing. The derivation of the state
equations of the cart-pole system stems from the solution of the associated Euler-
Lagrange equation. The model presented here follows [543, 544] and considers only
the dynamics of the pole (inverted pendulum). Denoting by θ the angle of the pole
112 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
(with respect to the vertical axis), m the mass of the pole, M the mass of the cart,
and l the length of the pole, one gets
where v = f (x, t) + g(x, t)u + d̃. Equation (3.41) describes a system in the
Brunovksy form. The flat output is taken to be y = x1 . It can be easily proven
that all state variables of the system of Eq. (3.41) and the associated control input can
be written as functions of the flat output y and its derivatives. Indeed, it holds that
x1 = y
x2 = ẏ (3.42)
u = g(y,
1
ẏ) [ ÿ − f (y, ẏ)]
3.3 Flatness-Based Adaptive Fuzzy Control for SISO Dynamical Systems 113
The following parameter values were chosen: m = 0.1 kg, M = 1.0 kgr, l = 0.5 m.
The basis functions used in the estimation of f (x, t) and g(x, t) were μ A j (x̂) =
x̂−c j 2
e( σ ) , j = 1, . . . , 5. Since there are two inputs x1 and ẋ1 and each fuzzy input
variable consists of 5 fuzzy sets there will be 25 fuzzy rules of the form:
(l)
The centers ci , i = 1, 2 and the variances v(l) of each rule are as follows:
(l) (l)
Rule c1 c2 v(l)
R (1) −0.3 −0.03 3
R (2) −0.3 −0.01 3
R (3) −0.3 0.00 3
R (4) −0.3 0.01 3
R (5) −0.3 0.03 3
R (6) −0.1 −0.03 3
··· ··· ··· ···
··· ··· ··· ···
R (25) 0.3 0.03 3
The estimation ĝ(x, t) was derived in a similar way. A measurable state vector was
assumed and used for the training of the neuro-fuzzy approximators. The positive
definite matrix P was found from the solution of the algebraic Riccati equation (3.23),
where Q was chosen to be a multiple of the identity matrix I2 .
Two different setpoints were studied: (i) a sinusoidal signal of amplitude 0.1 and
period T = 15 s, (ii) a seesaw setpoint of amplitude 0.15 and period T = 15 s.
At the beginning of the second half of the simulation time an additive sinusoidal
disturbance of amplitude A = 1.0 and period T = 7.5 s was applied to the system.
The approximations fˆ and ĝ were used in the control law, given by Eq. (3.7). The
position and velocity variations for the sinusoidal setpoint are depicted in Fig. 3.3a, b,
respectively.
The performance of the proposed flatness-based adaptive fuzzy control is also
tested in the tracking of a seesaw setpoint. The associated position and velocity vari-
ation are demonstrated in Fig. 3.4a, b, respectively. The control signal in the case of
tracking of a sinusoidal setpoint is shown in Fig. 3.5a, while the control signal when
tracking a seesaw setpoint is shown in Fig. 3.5b. Finally, the approximation of func-
tion g(x, t) in the case of tracking of a sinusoidal setpoint is shown in Fig. 3.6a (and
is marked as a dashed line), while when tracking a seesaw setpoint the approximated
function g(x, t) is shown in Fig. 3.6b.
114 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
1
0.5
position x1 of the pendulum
−0.5
−0.5
−1
−1
−1.5 −1.5
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t (sec) t (sec)
Fig. 3.3 a Tracking of a sinusoidal position setpoint (red line) by the angle of the pendulum.
b Tracking of a sinusoidal velocity setpoint (red line) by the angular velocity of the pendulum
1
1
position x1 of the pendulum
0.8
0.5
0.6
0.4 0
0.2 −0.5
0
−1
−0.2
−1.5
−0.4
−0.6 −2
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t (sec) t (sec)
Fig. 3.4 a Tracking of a seesaw position setpoint (red line) by the angle of the pendulum. b Tracking
of a seesaw velocity setpoint (red line) by the angular velocity of the pendulum
Comparing the proposed flatness-based adaptive fuzzy control method with other
adaptive fuzzy control approaches and with the analysis on neural adaptive control
methods given in the relevant bibliography (e.g., [167, 168, 205]), the following
can be noted: (i) the transformation of the initial nonlinear system into the lin-
earized Brunovksy (canonical) form does not require the computation of partial
derivatives or Lie derivatives, (ii) there is no need to make restrictive assumptions
about the number of truncated higher order terms in the linearization of the system’s
nonlinear model or about a bounded error in the linearization of the output of the
neural/fuzzy approximators, (iii) the number of adaptable parameters involved in the
training of the neural/fuzzy approximators remains small and there is no need to use
an excessive number of neural/fuzzy approximators to produce the control signal,
3.3 Flatness-Based Adaptive Fuzzy Control for SISO Dynamical Systems 115
(a) 50 (b) 60
50
40
40
30
30
20 20
10
10
0
−10
−10 −20
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t (sec) t (sec)
Fig. 3.5 a Control input of the pendulum when tracking a sinusoidal setpoint. b Control input of
the pendulum when tracking a seesaw setpoint
(a) 0 (b) 0
−0.1 −0.1
real vs estimated function g
−0.2 −0.2
−0.3 −0.3
−0.4 −0.4
−0.5 −0.5
−0.6 −0.6
−0.7 −0.7
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t (sec) t (sec)
Fig. 3.6 a Approximation of function g(x, t) of the pendulum model when tracking a sinusoidal
setpoint. b Approximation of function g(x, t) of the pendulum model when tracking a seesaw
setpoint
(iv) the tracking performance of the neuro-fuzzy control loop is evaluated with the use
of specific metrics (e.g., H∞ tracking performance), and (v) the proposed flatness-
based control method can also be extended to MIMO systems [426, 433].
It is also noted that the proposed flatness-based adaptive fuzzy control law assures
stability of the control loop and the asymptotic elimination of the tracking error.
Therefore, it can be ensured that θ̃g = θg − θg∗ →0 which means that ĝ(x)→g(x).
Provided that function g(x) does not become zero as long as x remains bounded, then
ĝ(x) will be also different from zero. For the cart-pole system described in Eqs. (3.39)
and (3.40) it holds that g(x) = 0 for x1 = θ = 0, however, for x1 ∈ (0, π ) this case is
116 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
not going to occur, even if ideally there is no function approximation error and ĝ(x)
coincides with g(x). In the generic case, to assure the avoidance of singularities in
the proposed control scheme one has to exclude singularity points from the reference
trajectory that the system’s state vector has to track.
Singularities may appear not only in the proposed adaptive fuzzy control scheme
but in all control systems that are based on static feedback linearization. For exam-
ple, linearization of the system through the use of a new control input of the form
v = f (x) + g(x)u means that u = g(x)−1 [v − f (x)], which does not exclude
the appearance of singularities. Therefore, singularities do not concern only the
proposed control method but the whole class of static feedback-based linearization
schemes. Modifications, as described in Chap. 2, can be introduced in the design of
the controller to prohibit the appearance of singularities, for example, a change of
coordinates that results in a new state-space representation that does not include any
points of singularity [340].
3.4.1 Overview
The results of Sect. 3.2 and can be extended toward adaptive fuzzy control based on
differential flatness theory for multi-input multi-output (MIMO) nonlinear dynamical
systems. There are numerous control problems for electromechanical systems with
MIMO nonlinear dynamics. The application of geometric techniques and mainly of
the feedback linearization method has enabled advancements in this research direc-
tion [428, 434, 524, 543, 544]. These techniques can be applied to nonlinear systems,
of uncertain or unknown dynamics. Based on the universal approximation theorem,
many stable adaptive fuzzy control schemes have been developed for unknown multi-
input multi-output (MIMO) dynamical systems [89, 91, 92, 289, 301, 501, 508, 526,
527]. The capability of neuro-fuzzy controllers to compensate for model parametric
uncertainties, external disturbances, as well as for incomplete measurement of the
system’s state vector has been analyzed in several studies [14, 110, 288, 484, 523,
524]. Moreover, results on neural and fuzzy control have been presented for MIMO
dynamical systems subjected to parametric changes and external perturbations, as
well as for systems in which model uncertainty may come from the inability to mea-
sure precisely the complete state vector. [101, 164, 302, 393, 527]. To improve the
robustness of these control schemes a supervisory control term is added to the control
signal, which can take the form of sliding mode control or H∞ control [92, 96, 291].
In this chapter an approach to the design of robust controllers for MIMO nonlinear
dynamical systems of unknown dynamics will be developed according to flatness-
based control theory. This approach extends the class of nonlinear systems to which
adaptive fuzzy control can be applied.
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 117
The transformation of MIMO nonlinear systems into the canonical (Brunovsky) from,
through the application of differential flatness theory, are explained first. Next, a new
control method for such systems are developed, also in accordance with differential
flatness theory. It is shown that the proposed control and filtering methods can be
efficiently applied to the model of the nonlinear MIMO dynamical systems.
The MIMO nonlinear dynamical system can be written as a state-space model
of the form ẋ = f (x, u). Moreover, the previous state-space description can be
transformed to the form of an affine in-the-input system by adding an integrator to
each input [55, 286]. In the latter case one obtains
m
ẋ = f (x) + i=1 gi (x)u i (3.44)
(ii) Lie Bracket: ad if g stands for a Lie Bracket which is defined recursively as
ad if g = [ f, ad i−1
f g] with ad f g = g and ad f g = [ f, g] = ∇g f − ∇fg.
0
with mj=1 v j = n, then y j = z 1, j for 1≤ j≤m are the 0-flat outputs which can be
written as functions of only the elements of the state vector x. In accordance with the
linearization criteria provided in Chap. 1, and to provide conditions for transforming
the system of Eq. (3.44) into the canonical form described in Eq. (3.45) the following
theorem has been stated [55]:
Theorem: For nonlinear systems described by Eq. (3.44) the following variables
are defined: (i) G 0 = span[g1 , . . . , gm ], (ii) G 1 = span[g1 , . . . , gm , ad f g1 , . . . ,
j
ad f gm ], . . . (k) G k = span{ad f gi for 0≤ j≤k, 1≤i≤m}. Then, the linearization
problem for the system of Eq. (3.44) can be solved if and only if: (1) The dimension
of G i , i = 1, . . . , k is constant for x ∈ X ⊆R n and for 1≤i≤n −1, (2) The dimension
of G n−1 if of order n, (3) The distribution G k is involutive for each 1≤k≤n − 2.
It is assumed now that after defining the flat outputs of the previous state-space
description of the nonlinear MIMO dynamics and after expressing the system’s state
variables and control inputs as functions of the flat output and of the associated
derivatives, the system can be transformed into the Brunovsky canonical form:
ẋ1 = x2
···
ẋr1 −1 = xr1
p
ẋr1 = f 1 (x) + j=1 g1 j (x)u j + d1
ẋr1 +1 = xr1 +2
··· (3.46)
ẋ p−1 = x p
p
ẋ p = f p (x) + j=1 g p j (x)u j + d p
y1 = x1
···
y p = xn−r p +1
y = [y1 , . . . , y p ]T is the output vector, f i are the drift functions and gi, j , i, j =
1, 2, . . . , p are smooth functions corresponding to the control input gains, while d j
is a variable associated to external disturbances. It holds that r1 + r2 + · · · + r p = n.
Having written the initial nonlinear system in the canonical (Brunovsky) form it
holds that
p
(r )
yi i = f i (x) + j=1 gi j (x)u j + d j (3.47)
Next, the following vectors and matrices can be defined: f (x) = [ f 1 (x), . . . ,
f n (x)]T , g(x) = [g1 (x), . . . , gn (x)]T , with gi (x) = [g1i (x), . . . , g pi (x)]T , A =
diag[A1 , . . . , A p ], and B = diag[B1 , . . . , B p ], C T = diag[C1 , . . . , C p ], d =
[d1 , . . . , d p ]T , where matrix A has the MIMO canonical form, i.e., with block-
diagonal elements
⎛ ⎞
0 1 ··· 0
⎜0 0 · · · 0 ⎟
⎜ ⎟
⎜ ⎟
Ai = ⎜ ... ... · · · ... ⎟
⎜ ⎟
⎝0 0 · · · 1 ⎠
(3.48)
0 0 · · · 0 r ×r
i i
BiT = 0 0 · · · 0 1 1×r
i
Ci = 1 0 · · · 0 0 1×r
i
ẋ = Ax + Bv + B d̃
(3.49)
y = Cx
where the control input is written as v = f (x) + g(x)u. The system of Eqs. (3.48)
and (3.49) is in controller and observer canonical form.
The reference setpoints for the system’s outputs y1 , . . . , y p are denoted as
y1m , . . . , y pm , thus for the associated tracking errors it holds that
e1 = y1 − y1m
e2 = y2 − y2m
(3.50)
···
e p = y p − y pm
The error vector of the outputs of the transformed MIMO system is denoted as
E 1 = [e1 , . . . , e p ]T
ym = [y1m , . . . , y pm ]T
(3.51)
···
ym(r ) = [y1m
(r )
, . . . , y (r ) T
pm ]
120 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
(r )
where yim denotes the r th order derivative of the ith reference output of the MIMO
dynamical system. Thus, one can also define the following vectors: (i) a vector
containing the state variables of the system and the associated derivatives, (ii) a
vector containing the reference outputs of the system and the associated derivatives
r −1 T
x = [x1 , . . . , x1r1 −1 , . . . , x p , . . . , x pp ] (3.52)
r1 −1 r −1
Ym = [y1m , . . . , y1m , . . . , y pm , . . . , y pm
p
]T (3.53)
while in a similar manner one can define a vector containing the tracking error of the
system’s outputs and the associated derivatives
r −1 T
e = Ym − x = [e1 , . . . , e1r1 −1 , . . . , e p , . . . , e pp ] (3.54)
It is assumed that matrix g(x) is a nonsingular one, i.e., g −1 (x) exists and is bounded
for all x ∈ Ux , where Ux ⊂R n is a compact set. In any case, the problem of singular-
ities in matrix g(x) can be handled by adding to ĝ(x) a small perturbation
ĝ which
results in avoidance of the singularity points and which is asymptotically eliminated
by the robustness of the control loop.
The objective of the adaptive fuzzy controller, denoted as u = u(x, e|θ ) is: (i)
all the signals involved in the controller’s design are bounded and it holds that
lim t→∞ e = 0, (ii) the H∞ tracking performance criterion is succeeded for a pre-
scribed attenuation level.
As in the case of adaptive fuzzy control for single-input dynamical systems
described in Sect. 3.2, successful tracking of the reference signal is denoted by the
H∞ criterion. This has already been formulated as [413, 524]:
T T
0 e T Qedt ≤ ρ 2 0 wd T wd dt (3.55)
where ρ is the attenuation level coefficient and denotes the aggregate effect of model
uncertainty and external perturbations.
The control signal of the MIMO nonlinear system which has been transformed into
the Brunovsky form as described by Eq. (3.49) contains the unknown nonlinear func-
tions f (x) and g(x), which can be approximated by
fˆ(x|θ f ) = Φ f (x)θ f
(3.56)
ĝ(x|θg ) = Φg (x)θg
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 121
where
with ξ if (x), ı = 1, . . . , n being the vector of kernel functions (e.g., normalized fuzzy
Gaussian membership functions), where
thus giving
⎛ ⎞
φ 1,1 1,2 1,N
f (x) φ f (x) · · · φ f (x)
⎜ 2,1 ⎟
⎜φ f (x) φ 2,2 2,N
f (x) · · · φ f (x)⎟
⎟
Φ f (x) = ⎜
⎜ ··· (3.59)
⎝ ··· ··· ··· ⎟ ⎠
φ n,1
f (x) φ n,2
f (x) · · · φ n,N
f (x)
with ξgi (x), ı = 1, . . . , N being the vector of kernel functions (e.g., normalized fuzzy
Gaussian membership functions), where
ξgi (x) = φgi,1 (x), φgi,2 (x), . . . , φgi,N (x) (3.62)
thus giving
⎛ ⎞
φg1,1 (x) φg1,2 (x) · · · φg1,N (x)
⎜ 2,1 2,2 2,N ⎟
⎜φ (x) φg (x) · · · φg (x)⎟
Φg (x) = ⎜ g ⎟ (3.63)
⎝ ··· ··· ··· ··· ⎠
φgn,1 (x) φgn,2 (x) · · · φgn,N (x)
It holds that
⎛ ⎞ ⎛ 1 p⎞
g1 g1 g12 · · · g1
⎜ g2 ⎟ ⎜ p⎟
⎟ ⎜ g2 g22 · · · g2 ⎟
1
g=⎜
⎝ · · ·⎠ = ⎜ ⎟ (3.67)
⎝· · · · · · · · · · · ·⎠
gn p
gn1 gn2 · · · gn
Using the above, one finally has the relation of Eq. (3.56), i.e., ĝ(x|θg ) = Φg (x)θg .
If the state variables of the system are available for measurement then a state-
feedback control law can be formulated as
(r )
u = ĝ −1 (x|θg )[− fˆ(x|θ f ) + ym + K T e + u c ] (3.68)
where fˆ(x|θ f ) and ĝ(x|θg ) are fuzzy models to approximate f (x) and g(x), respec-
tively. u c is a supervisory control term, e.g., H∞ control term which is used to com-
pensate for the effects of modeling inaccuracies and external disturbances. Using the
system’s state-space description of Eq. (3.49) the control term u c is defined as
u c = − r1 B T Pe (3.69)
Moreover, K T is the feedback gain matrix that assures that the characteristic poly-
nomial of matrix A − B K T will be a Hurwitz one.
Adaptive fuzzy control for MIMO nonlinear dynamical systems with the use
of differential flatness theory will be explained through an application example.
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 123
A 2-DOF rigid link robotic manipulator is considered. The dynamic model of the
robot is a differentially flat one. Indeed, it holds that
M11 M12 θ̈1 F1 (θ, θ̇ ) G 1 (θ ) T
+ + = 1 (3.70)
M21 M22 θ̈2 F2 (θ, θ̇ ) G 2 (θ ) T2
where M ∈ R 2×2 is the robot’s inertia matrix, F ∈ R 2×1 is the Coriolis and cen-
trifugal forces matrix and G ∈ R 2×1 is the gravitational effects matrix. The above
dynamic model is obtained with the use of the Euler-Lagrange method [160]. Equiv-
alently,
−1
θ̈1 M11 M12 F1 (θ, θ̇ )
=− −
θ̈2 M21 M22 F2 (θ, θ̇ )
−1
−1
(3.71)
M11 M12 G 1 (θ ) M11 M12 T1
− +
M21 M22 G 2 (θ ) M21 M22 T2
one obtains
θ̈1 N11 N12 F1 (θ, θ̇ )
=− −
θ̈2 N21 N22 F2 (θ, θ̇ )
(3.73)
N11 N12 G 1 (θ ) N11 N12 T1
− +
N21 N22 G 2 (θ ) N21 N22 T2
It holds that
where
It holds that
ẋ1 = x2
ẋ2 = f 1 (x) + g1 (x)u
(3.80)
ẋ3 = x4
ẋ4 = f 2 (x) + g2 (x)u
therefore, all system state variables can be written as functions of the flat output y
and its derivatives. The same holds for the control input u
x1 = [1 0]y T x2 = [1 0] ẏ T
(3.81)
x3 = [0 1]y T x4 = [0 1] ẏ T
Therefore, the considered robotic system is a differentially flat one. Next, taking
into account also the effects of additive disturbances to the joints of the robotic
manipulator the dynamic model becomes
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 125
where [u c1 u c2 ]T is a robust control term that is used for the compensation of the
model’s uncertainties as well as of the external disturbances and the vector of the
feedback gain is K iT = [k1i , k2i , . . . , kn−1
i , kni ].
Substituting Eq. (3.86) into Eq. (3.85) the closed-loop tracking error dynamics is
obtained
−1
ẍ1 f (x, t) g (x, t) ĝ1 (x, t)
= 1 + 1 ·
ẍ3 f 2 (x, t) g2 (x, t) ĝ2 (x, t)
(3.87)
ẍ1d fˆ1 (x, t) K 1T u c1 d̃1
− − e + +
ẍ3d fˆ2 (x, t) K 2T u c2 d̃2
Next, the following approximators of the unknown system dynamics are defined:
fˆ (x|θ f ) x ∈ R 4×1 fˆ1 (x|θ f ) ∈ R 1×1
fˆ(x) = ˆ1 (3.93)
f 2 (x|θ f ) x ∈ R 4×1 fˆ2 (x|θ f ) ∈ R 1×1
where l = 1, 2 and μ Ai (x) is the ith membership function of the antecedent (IF) part
j
of the lth fuzzy rule. Similarly, the following approximators of the unknown system
dynamics are defined:
ĝ (x|θg ) x ∈ R 4×1 ĝ1 (x|θg ) ∈ R 1×2
ĝ(x) = 1 (3.95)
ĝ2 (x|θg ) x ∈ R 4×1 ĝ2 (x|θg ) ∈ R 1×2
Mθ f = {θ f ∈ R h : ||θ f ||≤m θ f }
(3.97)
Mθg = {θg ∈ R h : ||θg ||≤m θg }
For the value of the approximation error defined in Eq. (3.90) that corresponds to the
optimal values of the weights vectors θ ∗f and θg∗ one has
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 127
w= f (x, t) − fˆ(x|θ ∗f ) +
−1 (3.98)
+ g(x, t) − ĝ(x|θg∗ ) ĝ(x, t) u
where
−1
wa = {[ fˆ(x, θ f ) − fˆ(x|θ ∗f )] + [ĝ(x, θg ) − ĝ(x|θg∗ )]} · ĝ(x, t) u (3.101)
θ̃ f = θ f − θ ∗f
(3.103)
θ̃g = θg − θg∗
Parameter γ1 is the learning rate used in the adaptation of the weights of the neu-
rofuzzy approximator for f (x), while parameter γ2 is the learning rate used in the
adaptation of the weights of the neurofuzzy approximation for g(x). It holds that
1 ˙T ˙T
V̇ = 21 ė T Pe + 21 e T P ė + γ1 θ̃ f θ̃ f + γ2 tr [θ̃g θ̃g ]
1 (3.105)
128 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
V̇ = 21 {e T (A − B K T )T + u cT B T +
+(w + d̃)T B T }Pe + 21 e T P{(A − B K T )e+ (3.109)
+Bu c + B(w + d̃)} + γ11 θ̃˙ Tf θ̃ f + γ12 tr [θ̃˙gT θ̃g ]
Assumption 1: For given positive definite matrix Q there exists a positive definite
matrix P, which is the solution of the following matrix equation:
T
(A − B K T ) P + P(A − B K T ) − PB( r2 − 1
ρ2
)B T P+Q=0 (3.111)
Substituting Eqs. (3.111) and (3.69) into V̇ yields after some operations
V̇ = 21 e T {−Q + PB( r2 − 1
ρ2
)B T P}e+
(3.112)
e T PB{− r1 B T Pe} + B T P(w + d̃) + γ11 θ̃˙ Tf θ̃ f + γ12 tr [θ̃˙gT θ̃g ]
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 129
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe + e T PB(w + d̃)+
(3.113)
1 ˙T ˙T
γ1 θ̃ f θ̃ f + γ2 tr [θ̃g θ̃g ]
1
θ̇ f = −γ1 Φ(x)T B T Pe
(3.115)
θ̇g = −γ2 Φ(x)T B T Peu T
where assuming N fuzzy rules and associated kernel functions the matrix dimensions
are θ f ∈ R N ×1 , θg ∈ R N ×2 , Φ(x) ∈ R 2×N , B ∈ R 4×2 , P ∈ R 4×4 and e ∈ R 4×1 .
Therefore, it holds that
or
Taking into account that u ∈ R 2×1 and e T PB(ĝ(x|θg ) − ĝ(x|θg∗ )) ∈ R 1×2 it holds
that
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe+
(3.122)
+e T PB(w + d̃) + e T P Bwα
w1 = w + d̃ + wα (3.123)
V̇ = − 21 e T Qe − 1 T
2ρ 2
e PBB T Pe + e T PBw1 (3.124)
Proof : The binomial (ρa − ρ1 b)2 ≥ 0 is considered. Expanding the left part of the
above inequality one gets
ρ2a2 + 1 2
ρ2
b − 2ab ≥ 0 ⇒ 21 ρ 2 a 2 + 2ρ1 2 b2 − ab ≥ 0 ⇒
(3.127)
ab − 1 2
2ρ 2
b ≤ 21 ρ 2 a 2 ⇒ 21 ab + 21 ab − 2ρ1 2 b2 ≤ 21 ρ 2 a 2
The following substitutions are carried out: a = w1 and b = e T PB and the previous
relation becomes
1 T
2 w1 B T Pe + 21 e T PBw1 − 1 T
2ρ 2
e PBB T Pe ≤ 21 ρ 2 w1T w1 (3.128)
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 131
The previous inequality is used in V̇ , and the right part of the associated inequality
is enforced
1 1
V̇ ≤ − e T Qe + ρ 2 w1T w1 (3.129)
2 2
As explained in the SISO systems control case, the attenuation coefficient ρ can be
chosen such that the right part of Eq. (3.129) is always upper bounded by 0. For
instance, it suffices at every iteration of the control algorithm to have
Again without knowledge of the uncertainties and disturbance term ||w1 || a suffi-
ciently small value of ρ can assure that the above inequality holds and thus that the
loop’s stability is assured. Actually, ρ should be given the least value which permits
to obtain a solution of the Riccati equation (3.111).
Equation (3.129) can be used to show that the H∞ performance criterion is satis-
fied. The integration of V̇ from 0 to T gives
T T T
0 V̇ (t)dt ≤ − 21 0 ||e||2 dt + 21 ρ 2 0 ||w1 ||2 dt ⇒
(3.131)
T T
2V (T ) + 0 ||e|| Q dt
2 ≤ 2V (0) + ρ 2 0 ||w1 ||2 dt
one gets
∞
0 ||e||2Q dt ≤ 2V (0) + ρ 2 Mw (3.133)
∞
Thus, the integral 0 ||e||2Q dt is bounded and according to Barbalat’s Lemma one
obtains lim t→∞ e(t) = 0.
In the presented adaptive fuzzy control approach there is no prior assumption
about bounded weights of the neuro-fuzzy approximators. It is the Lyapunov stability
analysis that assures that the tracking error of the control loop will follow the H∞
tracking performance criterion and that the weights values will remain bounded.
The only assumption needed for the design of the controller and for succeeding H∞
tracking performance for the control loop is that there exists a solution for a Riccati
equation associated to the linearized error dynamics of the differentially flat model.
This assumption is quite reasonable for several nonlinear systems, thus providing a
systematic approach to the design of reliable controllers for such systems.
132 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
It is also noted that the training of the neuro-fuzzy networks and the approximation
of the nonlinear system dynamics is efficient if the persistence of excitation condi-
tion holds. The problem of persistence of excitation in the training of neuro-fuzzy
approximators has been studied in several research articles [98, 179, 542]. According
to [179], a Persistence of Excitation condition on the input variables of the neuro-
fuzzy model holds if the input patterns belong to domains around the network nodes
centers. In the case of a neuro-fuzzy approximator that comprises Gaussian basis
functions centered on a regular mesh, these domains cover part of the input domain
volume. Therefore, sufficient coverage of the input space by basis functions makes
the Persistence of Excitation condition hold.
The linearized tracking error dynamics of Eq. (3.108) is considered again, i.e.,
ė = (A − B K T )e + Bu c + B(w + d̃)
The weight r determines how much the control signal should be penalized and the
weight ρ determines how much the disturbances influence should be rewarded in the
sense of a mini-max differential game. The H∞ control law is u(t) = − r1 B T Pe(t)
where P is the positive definite symmetric matrix derived from the algebraic Riccati
equation Eq. (3.111).
The parameter ρ in Eq. (3.134) is an indication of the closed-loop system robust-
ness. If the values of ρ > 0 are excessively decreased with respect to r , then the
solution of the Riccati equation is no longer a positive definite matrix. Consequently,
there is a lower bound ρmin of ρ for which the H∞ control problem has a solution.
The acceptable values of ρ lie in the interval [ρmin , ∞). If ρmin is found and used in
the design of the H∞ controller, the closed-loop system will have increased robust-
ness. Unlike this, if a value ρ > ρmin is used, then an admissible stabilizing H∞
controller will be derived but it will be a suboptimal one. The Hamiltonian matrix
A −( r2 − ρ12 )B B T
H= (3.135)
−Q −A T
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 133
provides a criterion for the existence of a solution of the Riccati equation Eq. (3.111).
A necessary condition for the solution of the algebraic Riccati equation to be a positive
semidefinite symmetric matrix is that H has no imaginary eigenvalues [132].
The performance of the proposed flatness-based adaptive fuzzy MIMO controller was
tested in the nonlinear dynamic model of the 2-DOF rigid-link robotic manipulator
(Fig. 3.7). The control loop is depicted in Fig. 3.8. The differentially flat model of the
robot and its transformation to the Brunovksy form has been analyzed in Sect. 3.4.
The flat outputs were taken to be the robot’s joint angles y1 = x1 and y2 = x3 . It has
been proven that all state variables of the robotic model and the associated control
inputs, i.e., the torques applied by the motors to the links’ joints can be written as
functions of the flat output [y1 , y2 ] and of the associated derivatives.
The state feedback gain was K ∈ R 2×4 . The basis functions used in the estimation
x̂−c j 2
of f i (x, t), i = 1, 2 and gi j (x, t), i = 1, 2, j = 1, 2 were μ A j (x̂) = e( σ ) , j =
1, . . . , 3. Since there are four inputs x1 , ẋ1 and x2 , ẋ2 and the associated definition
sets (universe of discourse) consist of 3 fuzzy sets, for the approximation of functions
f i (x, t) i = 1, 2, there will be 81 fuzzy rules of the form:
θ2
θ1 + θ2
θ1
X
O
134 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
Fig. 3.8 Flatness-based adaptive-fuzzy control loop for the robotic manipulator
In the considered fuzzy rule-base there are four input parameters in the antecedent
parts of the fuzzy rules, i.e., x1 = θ1 , x2 = θ̇1 , x3 = θ2 and x4 = θ̇2 . Each parameter
is partitioned into 3 fuzzy sets. Therefore, as mentioned before, by taking all possible
combinations between the fuzzy sets one has 34 = 81 fuzzy rules. The finer the par-
tition of the input variables into fuzzy sets is, the more accurate the approximation of
the nonlinear system dynamics by the neuro-fuzzy model is expected to be (although
some of the rules of the fuzzy rule base may not be sufficiently activated due to
little coverage of the associated region of the patterns space by input data). However,
considering a large number of fuzzy sets for each input variable induces the curse
of dimensionality, which means that there is an excessive and rather unnecessary
increase in the number of the adaptable parameters that constitute the neuro-fuzzy
model (Table 3.1).
The estimation of the control input gain functions ĝi j (x, t) i = 1, 2 was derived in
a similar way. The time step was taken to be 0.01 s. In the beginning of the training of
the neuro-fuzzy approximators their weights were initialized to zero. Moreover, the
elements of the robot’s state vector were also initialized to zero. The positive definite
matrix P ∈ R 4×4 stems from the solution of the algebraic Riccati equation (3.111),
for Q also positive definite.
Two different setpoints were studied: (i) a sinusoidal signal of variable amplitude
and frequency, (ii) a seesaw setpoint of amplitude 0.30 and period T = 20 s. The
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 135
approximations fˆ and ĝ were used in the derivation of the control law, given by
Eq. (3.115). To show the disturbance rejection capability of the proposed adaptive
fuzzy controller, at the beginning of the second half of the simulation time additive
sinusoidal disturbances of amplitude A = 0.5 and period T = 10 s were applied to
the robot’s joints. The external disturbances which appear as additive torques to the
robot’s joints can be due to a force F exerted on the manipulator’s end-effector (e.g.,
due to contact with a surface) and which through the relation T = J T F (J stands for
the Jacobian matrix of the robot’s kinematic model) generates torques on the joints
(Fig. 3.9).
In the simulation results that follow the position and velocity setpoints are noted
as continuous line while the position and velocity signals of the robot’s joints are
denoted as dashed lines. The position variation for the first joint of the robotic
2
20
1.5
state variable x2 of the robot
state variable x1 of the robot
1
10
0.5
0 0
−0.5
−10
−1
−1.5
−20
−2
−2.5 −30
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.9 a Tracking of sinusoidal position setpoint by joint 1 of the robot. b Tracking of sinusoidal
velocity setpoint by joint 1 of the robot
136 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
2
20
1.5
state variable x3 of the robot
0 0
−0.5
−10
−1
−1.5
−20
−2
−2.5 −30
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.10 a Tracking of sinusoidal position setpoint by joint 2 of the robot. b Tracking of sinusoidal
velocity setpoint by joint 2 of the robot
200 200
control signal of joint 1
100 100
0 0
−100 −100
−200 −200
−300 −300
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.11 a Tracking of a sinusoidal setpoint: Control input to joint 1. b Tracking of a sinusoidal
setpoint: Control input to joint 2
manipulator when tracking a sinusoidal setpoint is depicted in Fig. 3.9a, while the
velocity variation is shown in Fig. 3.9b. For the second joint of the 2-DOF robot the
tracking of the position setpoint is depicted in Fig. 3.10a while the tracking of the
velocity setpoint is shown in Fig. 3.10b. The control inputs (motor torques) applied
to the first and second joints of the robotic manipulator are shown in Fig. 3.11a, b,
respectively.
The performance of the proposed flatness-based adaptive fuzzy control is also
tested in the tracking of a seesaw setpoint. It is shown that after suitable tuning
of the feedback gain matrix K defined in Eq. (3.91) precise tracking of the seesaw
setpoint can be succeeded, with fast elimination of the tracking error. The controller
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 137
2
20
1.5
state variable x1 of the robot
0 0
−0.5
−10
−1
−1.5
−20
−2
−2.5 −30
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.12 a Tracking of a seesaw position setpoint by joint 1 of the robot. b Tracking of a seesaw
velocity setpoint by joint 1 of the robot
2
20
1.5
state variable x3 of the robot
1
10
0.5
0 0
−0.5
−10
−1
−1.5
−20
−2
−2.5 −30
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.13 a Tracking of a seesaw position setpoint by joint 2 of the robot. b Tracking of a seesaw
velocity setpoint by joint 2 of the robot
is also suitable for compensating the effects of the additive torques disturbance that
is applied to the robot’s joints.
The position variation of the first joint is demonstrated in Fig. 3.12a, while the
velocity variation is shown in Fig. 3.12b. Similarly, the tracking of the position
reference setpoint for the second joint is depicted in Fig. 3.13a, while the tracking of
the associated velocity setpoint is given in Fig. 3.13b. The control signal in the case
of tracking a seesaw setpoint by the two joints of the robotic manipulator is shown
in Fig. 3.14a, b, respectively.
138 3 Nonlinear Adaptive Control Based on Differential Flatness Theory
200 200
100 100
0 0
−100 −100
−200 −200
−300 −300
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.14 a Tracking of a seesaw setpoint: Control input to joint 1. b Tracking of a seesaw setpoint:
Control input to joint 2
1 1
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
t (sec) t (sec)
Fig. 3.15 a Approximation of function g22 (x, t) when tracking a sinusoidal setpoint. b Approxi-
mation of function g22 (x, t) when tracking a seesaw setpoint
Finally, the approximation of function g22 (x, t) in the case of tracking of a sinu-
soidal setpoint is shown in Fig. 3.15a (and is marked as a dashed line), while when
tracking a seesaw setpoint the approximated function g22 (x, t) is shown in Fig. 3.15b.
The RMSE (root mean square error) of the examined control loop is also calculated
(assuming the same parameters of the controller) in the case of tracking of previous
setpoints (a) sinusoidal setpoint of variable amplitude/frequency and (b) seesaw
setpoint. The results are summarized in Table 3.2. It can be observed that the examined
controller, apart from elimination of tracking error for the parameters of the robot’s
state vector, succeeds also in good transient performance. Taking into account that
no prior knowledge about the robot’s dynamics has been used by the controller, and
3.4 Flatness-Based Adaptive Fuzzy Control for MIMO Systems 139
at the first iterations of the control loop the controller undergoes self-tuning (learning
of values for specific gains) the transient characteristics of the control scheme can
be also deemed as quite satisfactory.
Chapter 4
Nonlinear Kalman Filtering Based
on Differential Flatness Theory
4.1 Introduction
A first type of uncertainties that flatness-based control has to cope with is the one
associated with parametric changes and external perturbations. This was analyzed
in Chap. 3, through flatness-based adaptive fuzzy control. A second type of uncer-
tainties in dynamical systems modeling and control has to do with the inability to
measure specific elements of the system’s state vector. To this end, the present chapter
introduces nonlinear Kalman Filtering based on differential flatness theory.
The chapter analyzes first the Derivative-free nonlinear Kalman Filter for SISO
and MIMO nonlinear dynamical systems. The considered nonlinear filtering scheme
which is based on differential flatness theory extends the class of systems to which
Kalman Filtering can be applied without the need for calculation of Jacobian matri-
ces. As explained, nonlinear systems satisfying the differential flatness property can
be written in the Brunovsky form via a transformation of their state variables and
control inputs. The Derivative-free nonlinear Kalman Filter consists of the Kalman
Filter recursion applied on the linearized equivalent of the system together with an
inverse transformation based on differential flatness theory that enables to retrieve
estimates for the state variables of the initial nonlinear model. In terms of accuracy
of the provided state estimation, the performance of the new filter is equivalent to
the one of the Unscented Kalman Filter and significantly improved comparing to
the Extended Kalman Filter. In terms of speed of computation, the Derivative-free
nonlinear Kalman Filter outperforms all other nonlinear estimation algorithms.
Aiming also at finding more efficient implementations of distributed nonlinear
filtering, in this chapter a derivative-free approach to Extended Information Filtering
is introduced and applied to state estimation-based control of nonlinear dynamical
systems, such as robots. First, the system undergoes a linearization transformation
and next state estimation is performed by applying the Kalman Filter recursion to the
linearized model. At a second level, the standard Information Filter is used to fuse the
state estimates obtained from local derivative-free Kalman filters running at the local
4.2.1 Overview
State estimation of nonlinear dynamical systems with the use of filters is a significant
topic in the area of control systems design, robotics, and mechatronics since it can
provide improved methods for sensorless control of electromechanical systems. For
nonlinear systems, and under the assumption of Gaussian noise, the Extended Kalman
Filter (EKF) is frequently applied for estimating the nonmeasurable state variables
through the processing of input and output sequences. This has been analyzed in [31,
51, 192, 229, 372, 417]. The Extended Kalman Filter is based on linearization of
the system using a first-order Taylor expansion [405, 408, 412, 422, 573]. Although
EKF is efficient in several estimation problems, it is characterized by cumulative
errors due to the local linearization assumption and this may affect the accuracy of
the state estimation or even risk the stability of the observer-based control loop.
Aiming at finding more efficient implementations of nonlinear Kalman Filtering,
in this chapter a derivative-free approach to Kalman filtering is introduced and applied
to state estimation-based control of a wide class of nonlinear electromechanical
systems. In the proposed derivative-free Kalman Filtering method the system is
first subject to a linearization transformation that is based on differential flatness
theory and next state estimation is performed by applying the standard Kalman
Filter recursion to the linearized model. Unlike EKF, the proposed method provides
estimates of the state vector of the nonlinear system without the need for derivatives
and Jacobian matrices calculation. By avoiding linearization approximations, the
proposed filtering method improves the accuracy of estimation of the system state
variables, and results in smooth control signal variations and in minimization of the
tracking error of the associated control loop. The chapter extends the results of [424,
429, 434] toward single-input single-output (SISO) and multi-input multi-output
(MIMO) nonlinear dynamical systems.
4.2 The Derivative-Free Nonlinear Kalman Filter 143
where x ∈ R m×1 is the system’s state vector, and z ∈ R p×1 is the system’s output.
Matrices A, B, and C can be time-varying and w(t), v(t) are uncorrelated white
Gaussian noises. The covariance matrix of the process noise w(t) is Q(t), while the
covariance matrix of the measurement noise is R(t). Then, the Kalman Filter is a
linear state observer which is given by
⎧
⎪ ˙
⎨x̂ = A x̂ + Bu + K [z − C x̂], x̂(t0 ) = 0
K (t) = PC T R −1 (4.2)
⎪
⎩
Ṗ = A P + P A T + Q − PC T R −1 CP
where x̂(t) is the optimal estimation of the state vector x(t) and P(t) is the covariance
matrix of the state vector estimation error with P(t0 ) = P0 . The Kalman Filter
consists of the system’s state equation plus a corrective term K [z−C x̂]. The selection
of gain K corresponds actually to the solution of an optimization problem. This is
expressed as the minimization of a quadratic cost functional and is performed through
the solution of a Riccati equation. In that case the observer’s gain K is calculated
by K = PC T R −1 considering an optimal control problem for the dual system
(A T , C T ), where the covariance matrix of the estimation error P is found by the
solution of a continuous-time Riccati equation of the form
Ṗ = AP + P A T + Q − PC T R −1 CP (4.3)
where matrices Q and R stand for the process and measurement noise covariance
matrices, respectively.
144 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
where the state x(k) is a m-vector, w(k) is a m-element process noise vector, and A
is a m × m real matrix. Moreover, the output measurement z(k) is a p-vector, C is
a p × m-matrix of real numbers, and v(k) is the measurement noise. It is assumed
that the process noise w(k) and the measurement noise v(k) are uncorrelated.
Now the problem of interest is to estimate the state x(k) based on the sequence
of output measurements z(1), z(2), . . . , z(k). The initial value of the state vector
x(0), and the initial value of the error covariance matrix P(0) is unknown and an
estimation of it is considered, i.e., x̂(0) is a guess of E[x(0)] and P̂(0) is a guess of
Cov[x(0)].
For the initialization of matrix P one can set P̂(0) = λI , with λ > 0. The
state vector x(k) has to be estimated taking into account x̂(0), P̂(0), and the output
measurements Z = [z(1), z(2), . . . , z(k)]T , i.e., x̂(k) = αn (x̂(0)), P̂(0), Z (k)).
This is a linear minimum mean squares estimation problem (LMMSE) formulated
as x̂(k + 1) = an+1 (x̂(k), z(k + 1)). The process and output noise are white and their
covariance matrices are given by: E[w(i)wT ( j)] = Qδ(i − j) and E[v(i)vT ( j)] =
Rδ(i − j).
Using the above, the discrete-time Kalman filter can be decomposed into two parts:
(i) time update (prediction stage), and (ii) measurement update (correction stage).
The first part employs an estimate of the state vector x(k) made before the output
measurement z(k) is available (a priori estimate). The second part estimates x(k)
after z(k) has become available (a posteriori estimate). The following measurement
processing stages are distinguished:
• When the set of measurements Z − = {z(1), . . . , z(k − 1)} is available. From
Z − an a priori estimation of x(k) is obtained which is denoted by x̂ − (k) and this
denotes the estimate of x(k) given Z − .
• When z(k) is available, the output measurements set becomes Z = {z(1), . . . , z(k)},
where x̂(k) and this denotes the estimate of x(k) given Z .
The associated estimation errors are defined by e− (k) = x(k) − x̂ − (k)= the a priori
error, and e(k) = x(k)− x̂(k)= the a posteriori error. The estimation error covariance
matrices associated with x̂ − (k) and x̂(k) are defined as P − (k) = Cov[e− (k)] =
E[e− (k)e− (k) ] and P(k) = Cov[e(k)] = E[e(k)e T (k)], respectively [229]. From
T
the definition of the trace of a matrix, the mean square error of the estimates can
be written as MSE(x̂ − (k)) = E[e− (k)e− (k) ] = tr(P − (k)) and MSE(x(k)) =
T
time update:
State estimation can be also performed for nonlinear dynamical systems using the
Extended Kalman Filter recursion. The following nonlinear state model is considered
[408, 422]:
x(k + 1) = φ(x(k)) + L(k)u(k) + w(k)
(4.7)
z(k) = γ (x(k)) + v(k)
where x ∈ R m×1 is the system’s state vector and z ∈ R p×1 is the system’s out-
put, while w(k) and v(k) are uncorrelated, Gaussian zero mean noise processes with
covariance matrices Q(k) and R(k), respectively. The operators φ(x) and γ (x) are the
vectors φ(x) = [φ1 (x), φ2 (x), . . ., φm (x)]T , and γ (x) = [γ1 (x), γ2 (x), . . . , γ p (x)]T ,
respectively. It is assumed that φ and γ are sufficiently smooth in x so that each
one has a valid series Taylor expansion. Following a linearization procedure, φ is
expanded into Taylor series about x̂:
where x̂ − (k) is the estimation of the state vector x(k) before measurement at the
kth instant to be received and x̂(k) is the updated estimation of the state vector after
measurement at the kth instant has been received. The Jacobian Jγ (x) is
⎛ ∂γ1 ∂γ1 ∂γ1 ⎞
∂ x1 ∂ x2 ··· ∂ xm
⎜ ∂γ ∂γ ⎟
⎜ 2 2 ··· ∂γ2 ⎟
∂γ ⎜ ∂ x1 ∂ x2 ∂ xm ⎟
Jγ (x) = |x=x̂ − (k) = ⎜
⎜ .. .. ..
⎟
.. ⎟ (4.11)
∂x ⎜ .
⎝ . . . ⎟
⎠
∂γ p ∂γ p ∂γ p
∂ x1 ∂ x2 ··· ∂ xm
Now, the EKF recursion is as follows: First, the time update is considered: by x̂(k)
the estimation of the state vector at time instant k is denoted. Given initial conditions
x̂ − (0) and P − (0) the recursion proceeds as:
• Measurement update. Acquire z(k) and compute:
K (k) = P − (k)JγT (x̂ − (k))·[Jγ (x̂ − (k))P − (k)JγT (x̂ − (k)) + R(k)]−1
x̂(k) = x̂ − (k) + K (k)[z(k) − γ (x̂ − (k))] (4.13)
P(k) = P − (k) − K (k)Jγ (x̂ − (k))P − (k)
for electric machines tasks have been given in [138] where an EKF-based fault diag-
nosis scheme for actuators has been studied, in [578] where the EKF has been pro-
posed for residual generation and diagnosis of incipient faults in DC electric motors,
and in [13] where the EKF has been used in fault diagnosis of an electrohydraulic
actuation system.
Next, the functioning of the Unscented Kalman Filter (UKF) will be explained.
The UKF belongs to a wider class of nonlinear estimators known as Sigma-Point
Kalman Filters (SPKF) [224, 225, 231]. Unlike EKF, in the Unscented Kalman Filter
no analytical Jacobians of the system equations need to be calculated. Some basic
operations performed in the UKF algorithm are summarized as follows:
(1) Denoting the current state mean as x̂, √ a set of 2n + 1 sigma points is taken
from√the columns of the n × n matrix (n + √ λ)Pxx as follows: x = x̂, x =
0 i
(m) λ (c) λ
W0 = (n+λ) W0 = (n+λ)+(1−α 2 +b)
(4.15)
(m) (c)
Wi = 1
2(n+λ) Wi = 2(n+λ)1
148 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
Measurement update: Obtain the new output measurement z(k) and compute the pre-
dicted mean ẑ(k) and covariance of the measurement Pzz (k), and the cross-covariance
of the state and measurement Pxz (k)
Then compute the filter gain K (k), the state mean x̂(k), and the covariance Pxx (k),
conditional to the measurement y(k)
The filter starts from the initial mean m(0) and covariance Px x (0). The stages of
state vector estimation with the use of the Unscented Kalman Filter algorithm are
depicted in Fig. 4.2.
4.2 The Derivative-Free Nonlinear Kalman Filter 149
Differential flatness theory will be used for arriving at a filtering method that will have
improved performance compared to the Extended Kalman Filter. It will be shown
that through a nonlinear transformation it is possible to design a state estimator for
a class of nonlinear systems, which can substitute for the Extended Kalman Filter.
The results will be towards derivative-free Kalman Filtering for nonlinear systems.
The following continuous-time nonlinear single-output system is considered [334,
335, 434].
p
ẋ = f (x) + q0 (x, u) + i=1 θi qi (x, u), or
ẋ = f (x) + q0 (x, u) + Q(x, u)θ x ∈ R n , u ∈ R m , θ ∈ R p (4.21)
z = h(x), z ∈ R
with ⎛ ⎞
0 1 0 ··· 0
⎜0 0 1 ··· 0⎟
⎜ ⎟
Ac = ⎜ . .. .. .. .. ⎟ (4.24)
⎝ .. . . . .⎠
0 0 0 ··· 0
Cc = 1 0 0 · · · 0 (4.25)
and ad f g = [ f, g] = ∇g f − ∇ f g.
(iii) [qi , adif g] = 0, 0 ≤ i ≤ p, 0 ≤ j ≤ n − 2 ∀ u ∈ R m .
(iv) the vector fields adif g, 0 ≤ i ≤ n − 1 are complete, in which g is the vector
field satisfying ⎛ ⎞ ⎛ ⎞
dh 0
⎜ .. ⎟ ⎜ .. ⎟
<⎝ . ⎠ , g >= ⎝ . ⎠ (4.26)
d(Ln−1
f h) 1
∂h
where dh has been defined as dh = ∂x = [ ∂∂h ∂h ∂h
x1 , ∂ x2 , . . . , ∂ xn ]. Then for every
parameter vector θ , the system
p
ζ̂ = Ac ζ̂ + ψ0 (z, u) + i=1 θi ψi (z, u) + K (z − Cc ζ̂ )
(4.27)
x̂ = T −1 (ζ̂ )
4.2 The Derivative-Free Nonlinear Kalman Filter 151
is an asymptotic observer for a suitable choice of K provided that the state x(t) is
bounded, with estimation error dynamics
⎛ ⎞
−k1 1 0 ··· 0
⎜ −k2 0 1 ··· 0⎟
⎜ ⎟
ė = (Ac − K Cc )e = ⎜
⎜ ··· ··· ··· ··· · · ·⎟
⎟e (4.28)
⎝−kn−1 0 0 ··· 1⎠
−kn 0 0 ··· 0
Since Eq. (4.27) provides an asymptotic observer for the initial nonlinear system of
Eq. (4.21) one can consider a case in which the observation error gain matrix K can
be provided by the Kalman Filter equations given initially in the continuous-time
KF formulation, or in discrete-time form by Eqs. (10.62) and (10.63). The following
single-input single-output nonlinear dynamical system is considered
where z = x is the system’s output, and f (x, t), g(x, t) are nonlinear functions. It
can be noticed that the system of Eq. (4.29) belongs to the general class of systems of
Eq. (4.21). Moreover, Eq. (4.29) is the canonical form one obtains after the application
of a differential flatness theory-based transformation. Assuming the transformation
ζi = x (i−1) , i = 1, . . . , n, and x (n) = f (x, t) + g(x, t)u(x, t) = v(ζ, t), i.e.,
ζ̇n = v(ζ, t), one obtains the linearized system of the form
ζ̇1 = ζ2
ζ̇2 = ζ3
··· ··· (4.30)
ζ̇n−1 = ζn
ζ̇n = v(ζ, t)
152 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
The system of Eqs. (4.31) and (4.32) has been written in the form of Eq. (4.23),
which means that Eq. (4.27) is the associated asymptotic observer. Therefore, the
observation gain K appearing in Eq. (4.27) can be found using either linear observer
design methods (in that case the elements of the observation error gain matrix K
have fixed values), or the recursive calculation of the continuous-time Kalman Filter
gain described in Sect. 4.2.2.2. If the discrete-time Kalman Filter is to be used then
one has to apply the recursive formulas of Eqs. (4.5) and (4.6) on the discrete-time
equivalent of Eqs. (4.31) and (4.32).
The model of the linear DC motor shown in Fig. 4.3 is described by the set of equations
L I˙ = −ke ω − R I + V
(4.33)
J ω̇ = ke I − kd ω − Γd
Next, control for a nonlinear DC motor model will be presented. Usually, the DC
motor model is considered to be linear by neglecting the effect of armature reaction
or by assuming that the compensating windings remove this effect. Introducing the
armature reaction leads to a nonlinear system. In that case the dynamic model of the
4.2 The Derivative-Free Nonlinear Kalman Filter 153
ẋ1 = x2
ẋ2 = k1 x2 + k2 x3 + k3 x32 + k4 T1 (4.37)
ẋ3 = k5 x2 + k6 x2 x3 + k7 x3 + k8 u
where T1 the load torque and u is the terminal voltage. From the second row of
Eq. (4.37) one obtains,
where, f¯(x) = k1 ẋ2 + k2 ẋ3 + 2k3 k5 x2 x3 + 2k3 k6 x2 x32 + 2k3 k7 x32 , and ḡ(x) =
2k3 k8 x3 . Taking k1 = 0 and considering the torque T1 as external disturbance, the
nonlinear DC motor model of Eq. (4.37) becomes
ẋ1 = x2
ẋ2 = k2 x3 + k3 x32 (4.40)
ẋ3 = k5 x2 + k6 x2 x3 + k7 x3 + k8 u
Then, setting the output to be y = x1 one can see that all state variables xi , i = 1, 2, 3
and the control input u can be expressed as functions of this output and its derivatives.
Indeed it holds
x1 = y
x2 = ẏ
−k2 + |k22 +4k3 ÿ| (4.41)
x3 = 2k3
u = ḡ(x)
1
[y (3) − f¯(x)]
The aforementioned system of Eq. (4.40) can be written in the Brunovsky (canonical)
form: ⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
ẏ1 010 y1 0
⎝ ẏ2 ⎠ = ⎝0 0 1⎠ ⎝ y2 ⎠ + ⎝0⎠ v (4.42)
ẏ3 000 y3 1
where v = f¯(x, t) + ḡ(x, t)u. Denoting the state vector of the transformed descrip-
tion of the system as y = [y1 , y2 , y)3]T , the measurement (observation) equation
becomes
z 1 = y1 ⇒
(4.43)
z 1 = [1 0 0]y
(r )
Finally, using a control input of the form u = ḡ(x̂, t)−1 [− f¯(x̂, t)+ ym + K cT e +u c ]
it is possible to make the electric motor’s angle track any desirable setpoint.
The dynamic model of the DC motor given in Eq. (4.39) is considered. The control
law described above is used to make the motor track the desirable trajectory. It is
assumed that only encoder measurements are available (angle θ ), thus the state vector
x used in the control law is not directly measured and has to be estimated through
filtering. First, state estimation-based control is implemented with the use of the
4.2 The Derivative-Free Nonlinear Kalman Filter 155
Extended Kalman Filter. The associated covariance matrices Q and R were taken to
be diagonal, both having diagonal elements equal to 10−3 . The feedback controller’s
gain appearing was taken to be K T = [60, 30, 20]. The initial value of the system’s
state vector was x(0) = [0, 0, 0]T while the initial value of the estimated state vector
was x̂(0) = [1, 1, 1]T . The estimation error covariance matrix P ∈ 3×3 and the
EKF gain K ∈ 3×1 were used in Eqs. (4.5) and (4.6). The sampling period was
taken to be Ts = 0.01 s. The obtained results for the case of a seesaw and a sinusoidal
setpoint are shown in Figs. 4.4, 4.5, 4.6 and 4.7. The reference setpoint in denoted
by the red line, the state vector variables of the motor are denoted by the dashed blue
line, while the estimated state vector elements are described by the green line.
Next, state estimation-based control was implemented with the use of the
Unscented Kalman Filter. As explained the UKF does not require linearization of
the system’s dynamical model and does not incorporate the calculation of Jacobians.
Since linearization errors are avoided, the Unscented Kalman Filter results in more
robust state estimation thus also making more accurate the tracking of the reference
setpoints by the associated control loop. It can be also noticed that the variations of
the control signal are quite smooth. The associated simulation results are depicted
in Figs. 4.8, 4.9, 4.10 and 4.11.
Finally, in the case of state estimation-based control with the use of the derivative-
free nonlinear Kalman Filter (DKF), the estimation of the state vector x̂ is obtained
by applying derivative-free Kalman Filtering to the linearized model of Eqs. (4.44)
and (4.45), according to Eqs. (4.5) and (4.6). Starting from Eq. (4.39) and defining
ζ1 = θ , ζ2 = θ̇ , and ζ3 = θ̈ the state-space equations of the motor are written as:
(a) 1.5
EKF − x1 (b) 1.5
EKF − x1
1 1
position x of the motor (rad)
0.5 0.5
0 0
1
−0.5 −0.5
−1 −1
−1.5 −1.5
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.4 Variable x1 of the state vector in state estimation with use of the Extended Kalman Filter,
a when tracking a seesaw setpoint, b when tracking a sinusoidal setpoint
156 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
(a) 3
EKF − x2 (b) 3
EKF − x2
2 2
1 1
0 0
2
−1 −1
−2 −2
−3 −3
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.5 Variable x2 of the state vector in state estimation with use of the Extended Kalman Filter,
a when tracking a seesaw setpoint, b when tracking a sinusoidal setpoint
(a) 6
EKF − x3 (b) 6
EKF − x3
acceleration x of the motor (rad/sec2)
acceleration x of the motor (rad/sec2)
4 4
2 2
0 0
−2 −2
3
−4 −4
−6 −6
−8 −8
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.6 Variable x3 of the state vector in estimation performed with use of the Extended Kalman
Filter, a when tracking a seesaw setpoint, b when tracking of a sinusoidal setpoint
⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
ζ̇1 010 ζ1 0
⎝ζ̇2 ⎠ = ⎝0 0 1⎠ ⎝ζ2 ⎠ + ⎝0⎠ v(ζ, t) (4.44)
ζ̇3 000 ζ3 1
z= 100 ζ (4.45)
where v(ζ, t) = f¯(x, t) + ḡ(x, t)u(x, t). To perform state estimation-based (sensor-
less) control only measurements of the rotor’s angle θ were used. As in the EKF case
the associated covariance matrices Q and R were taken to be diagonal, both having
4.2 The Derivative-Free Nonlinear Kalman Filter 157
(a)
control signal Vin of the motor (v.u.)
400
EKF − control input (b) 400 EKF − control input
200 200
100 100
0 0
−100 −100
−200 −200
−300 −300
−400 −400
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.7 Control signal of the Extended Kalman Filter-based control loop, a when tracking a seesaw
setpoint, b when tracking a sinusoidal setpoint
(a) 1.5
UKF − x1 (b) 1.5
UKF − x1
position x1 of the motor (rad)
1 1
position x of the motor (rad)
0.5 0.5
0 0
1
−0.5 −0.5
−1 −1
−1.5 −1.5
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.8 Variable x1 of the state vector in state estimation with use of the Unscented Kalman Filter,
a when tracking a seesaw setpoint, b when tracking a sinusoidal setpoint
(a)
velocity x2 of the motor (rad/sec)
3
UKF − x2 (b) 3
UKF − x2
1 1
0 0
−1 −1
−2 −2
−3 −3
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.9 Variable x2 of the state vector in state estimation with use of the Unscented Kalman Filter,
a when tracking a seesaw setpoint, b when tracking a sinusoidal setpoint
4 4
2 2
0 0
−2 −2
3
−4 −4
−6 −6
−8 −8
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.10 Variable x3 of the state vector in estimation was performed with use of the Unscented
Kalman Filter, a when tracking a seesaw setpoint, b when tracking of a sinusoidal setpoint
as the accuracy of tracking of the state estimation-based control loops was the Mean
Square Error (MSE). Alternatively, the Cramer–Rao Lower Bound (CRLB) could
have been considered [574].
It can be observed that comparing to EKF the derivative-free nonlinear Kalman
Filter results in smaller estimation errors for the system’s state variables and con-
sequently the accuracy of tracking of the reference setpoints is also improved. The
ranges of variation of the control signal are also kept small. Using the derivative-
free nonlinear Kalman Filter the performance of the state estimation-based control
is comparable to the one obtained when using the Unscented Kalman Filter for
4.2 The Derivative-Free Nonlinear Kalman Filter 159
200 200
100 100
0 0
−100 −100
−200 −200
−300 −300
−400 −400
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.11 Control signal of the Unscented Kalman Filter-based control loop, a when tracking a
seesaw setpoint, b when tracking a sinusoidal setpoint
(a) 1.5
DKF − x1 (b) 1.5
DKF − x1
1 1
position x of the motor (rad)
position x of the motor (rad)
0.5 0.5
0 0
1
1
−0.5 −0.5
−1 −1
−1.5 −1.5
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.12 Variable x1 of the state vector in state estimation with use of the derivative-free Kalman
Filter, a when tracking a seesaw setpoint, b when tracking a sinusoidal setpoint
estimating the nonmeasurable state vector parameters. It can be also noted that in
terms of computation cost the proposed derivative-free nonlinear Kalman Filtering
(DKF) approach has some advantages comparing to the Unscented Kalman Filter:
since the DKF follows the fast recursion of the standard Kalman Filter it does not
need to update sigma points at each iteration, to propagate the nonlinear estimation
through these sigma points and to calculate cross-covariance matrices.
The performance of the nonlinear filtering schemes was also tested in state
estimation-based control for the tracking of a setpoint of variable amplitude, fre-
quency, and phase. The associated results for the Extended Kalman Filter have
160 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
(a)
velocity x 2 of the motor (rad/sec)
3
DKF − x2 (b) 3
DKF − x2
1 1
0 0
−1 −1
−2 −2
−3 −3
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.13 Variable x2 of the state vector in state estimation with use of the derivative-free Kalman
Filter, a when tracking a seesaw setpoint, b when tracking a sinusoidal setpoint
(a) 6
DKF − x3 (b) 6
DKF − x3
acceleration x of the motor (rad/sec 2)
4 4
2 2
0 0
−2 −2
3
−4 −4
−6 −6
−8 −8
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.14 Variable x3 of the state vector in estimation was performed with use of the derivative-free
Kalman Filter, a when tracking a seesaw setpoint, b when tracking of a sinusoidal setpoint
been given in Figs. 4.16 and 4.17, for the Unscented Kalman Filter have been given
in Figs. 4.18 and 4.19, while for derivative-free Kalman filter have been given in
Figs. 4.20 and 4.21. Again the reference setpoint is denoted by the red line, the state
vector variables of the motor are denoted by the dashed blue line, while the estimated
state vector elements are described by the green line. The improved performance suc-
ceeded by the derivative-free nonlinear Kalman Filter (DKF) can be noticed once
more.
4.2 The Derivative-Free Nonlinear Kalman Filter 161
200 200
100 100
0 0
−100 −100
−200 −200
−300 −300
−400 −400
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.15 Control signal of the derivative-free Kalman Filter-based control loop, a when tracking
a seesaw setpoint, b when tracking a sinusoidal setpoint
(a) 1.5
EKF − x1 (b) EKF − x2
3
0.5 1
0 0
1
−0.5 −1
−1 −2
−1.5 −3
0 10 20 30 40 50 60 0 10 20 30 40 50 60
t (sec) t (sec)
Fig. 4.16 State estimation-based control with the use of the Extended Kalman Filter for tracking a
sinusoidal setpoint of variable amplitude, frequency, and phase: a state variable x1 , b state variable
x2
(a) 6
EKF − x3 (b) EKF − control input
400
acceleration x of the motor (rad/sec2)
4 300
200
2
100
0
0
−2
3
−100
−4
−200
−6 −300
−8 −400
0 10 20 30 40 50 60 0 10 20 30 40 50 60
t (sec) t (sec)
Fig. 4.17 State estimation-based control with the use of the Extended Kalman Filter for tracking a
sinusoidal setpoint of variable amplitude, frequency, and phase: a state variable x3 , b control input
Kalman Filtering approach and the transformation to the observer canonical form
make possible the application of smoothing algorithms (initially designed for linear
systems) that compensate for the effects of delayed sensor measurements in the
control loop. This can be particularly useful in the case of networked control and
fault diagnosis schemes where communication delays and packet drops can affect
the accuracy of the state estimation.
4.2 The Derivative-Free Nonlinear Kalman Filter 163
0.5 1
0 0
1
2
−0.5 −1
−1 −2
−1.5 −3
0 10 20 30 40 50 60 0 10 20 30 40 50 60
t (sec) t (sec)
Fig. 4.18 State estimation-based control with the use of the Unscented Kalman Filter for tracking a
sinusoidal setpoint of variable amplitude, frequency, and phase: a state variable x1 , b state variable
x2
4 300
200
2
100
0
0
−2
3
−100
−4
−200
−6 −300
−8 −400
0 10 20 30 40 50 60 0 10 20 30 40 50 60
t (sec) t (sec)
Fig. 4.19 State estimation-based control with the use of the Unscented Kalman Filter for tracking
a sinusoidal setpoint of variable amplitude, frequency, and phase: a state variable x3 , b control input
The results of Chap. 3 about transformation of MIMO nonlinear systems into the
canonical form will be generalized toward Derivative-free Kalman Filtering for
MIMO nonlinear dynamical systems. The proposed method for derivative-free
Kalman Filtering for MIMO nonlinear systems will be analyzed through an applica-
tion example. A 2-DOF rigid-link robotic manipulator is considered. The dynamic
164 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
(a) 1.5
DKF − x1 (b) 3
DKF − x2
0.5 1
0 0
1
2
−0.5 −1
−1 −2
−1.5 −3
0 10 20 30 40 50 60 0 10 20 30 40 50 60
t (sec) t (sec)
Fig. 4.20 State estimation-based control with the use of the derivative-free nonlinear Kalman Filter
for tracking a sinusoidal setpoint of variable amplitude, frequency, and phase: a state variable x1 ,
b state variable x2
4 300
200
2
100
0
0
−2
3
−100
−4
−200
−6 −300
−8 −400
0 10 20 30 40 50 60 0 10 20 30 40 50 60
t (sec) t (sec)
Fig. 4.21 State estimation-based control with the use of the derivative-free nonlinear Kalman Filter
for tracking a sinusoidal setpoint of variable amplitude, frequency, and phase: a state variable x3 ,
b control input
model of the robot has been analyzed in Chap. 3 and it has been shown that (i) it satis-
fies differential flatness properties, (ii) it can be transformed into the linear canonical
form.
As previously shown this robotic system is a differentially flat one. Next, consid-
ering also the effects of additive disturbances to the joints of the robotic manipulator
the dynamic model becomes
4.2 The Derivative-Free Nonlinear Kalman Filter 165
Moreover, considering that the model’s structure and parameters are known and
that the additive input disturbance d̃i , i = 1, 2, . . . are also known, the following
feedback linearizing control input is defined
−1
g1 (x̂, t) ẍ1d f 1 (x̂, t) K 1T u c1
u= − − e+ (4.48)
g2 (x̂, t) ẍ3d f 2 (x̂, t) K 2T u c2
Using the matrices of Eq. (4.49) one obtains the Brunovsky form of the MIMO robot
model dynamics
ẋ = Ax + Bv
(4.50)
y = Cx
Prior to performing discrete-time Kalman Filtering for the model of Eq. (4.50) matri-
ces A, B, and C of Eq. (4.49) are turned into discrete-time ones using common dis-
cretization methods. These discrete-time matrices are denoted as Ad , Bd , and Cd ,
respectively.
For the robotic model of Eqs. (4.50) and (4.51) state estimation can be performed
using the standard Kalman Filter recursion, as described in Eqs. (4.5) and (4.6).
166 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
The performance of the proposed derivative-free nonlinear Kalman Filter was tested
in the nonlinear state estimation-based control for a 2-DOF rigid-link robotic manip-
ulator (Fig. 3.7). The differentially flat model of the robot and its transformation to
the Brunovksy form has been analyzed in Sect. 4.2.5. The flat outputs were taken to
be the robot’s joint angles y1 = x1 and y2 = x3 . It has been proven that all state
variables of the robotic model and the associated control inputs, i.e., the torques
applied by the motors to the links’ joints can be written as functions of the flat output
[y1 , y2 ] and of the associated derivatives.
The position and velocity variations for the sinusoidal setpoint for the first joint of
the robotic manipulator are depicted in Fig. 4.22. For the second joint of the 2-DOF
robot the tracking of the position and velocity setpoints is depicted in Fig. 4.23. The
control inputs (motor torques) applied to the first and second joint of the robotic
manipulator are shown in Fig. 4.24.
The performance of the state estimation-based control has been also tested in the
tracking of a seesaw setpoint. The position and velocity variations of the first joint
are demonstrated in Fig. 4.25. Similarly, the tracking of the position and velocity
reference setpoints for the second joint are depicted in Fig. 4.26. The control signal
in the case of tracking of a tracking a seesaw setpoint by the two joints of the robotic
manipulator is shown in Fig. 4.27. The estimated state variables were denoted as
green line whereas the real state variables were denoted as blue lines.
1.5
1
state variable x of the robot
1
0.5
0.5
1
0 0
−0.5
−0.5
−1
−1
−1.5
−1.5 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.22 Derivative-free Kalman Filtering: a tracking of a sinusoidal position setpoint (dashed
line) by the angle of the first joint of the robot (continuous line), b tracking of a sinusoidal velocity
setpoint (dashed line) by the angular velocity of the first joint of the robot (continuous line)
4.2 The Derivative-Free Nonlinear Kalman Filter 167
1.5
1
1
0.5
0.5
0 0
−0.5
−0.5
−1
−1
−1.5
−1.5 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.23 Derivative-free Kalman Filtering: a tracking of a sinusoidal position setpoint (dashed
line) by the angle of the second joint of the robot (continuous line), b tracking of a sinusoidal
velocity setpoint (dashed line) by the angular velocity of the second joint of the robot (continuous
line)
(a) 60 (b) 26
24
55
22
control signal of joint 1
50
20
45 18
16
40
14
35
12
30 10
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.24 Tracking of a sinusoidal setpoint by the 2-DOF robotic manipulator with the use of the
derivative-free Kalman Filter: a control input (motor torque) applied to the first joint, b control
input (motor torque) applied to the second joint
Moreover, the derivative-free nonlinear Kalman Filter was tested against the
Extended Kalman Filter which is among the most commonly used estimation meth-
ods for nonlinear dynamical systems. In the latter case, the position and velocity
variations for the sinusoidal setpoint for the first joint of the robotic manipulator are
depicted in Fig. 4.28. For the second joint of the 2-DOF robot the tracking of the
168 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
0.5 1.5
0.4 1
0.3 0.5
0.2 0
0.1 −0.5
0 −1
−0.1 −1.5
−0.2 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.25 Derivative-free Kalman Filtering: a tracking of a seesaw position setpoint (dashed line)
by the angle of the first joint of the robot (continuous line), b tracking of a seesaw velocity setpoint
(dashed line) by the angular velocity of the first joint of the robot (continuous line)
0.5 1.5
state variable x of the robot
state variable x of the robot
0.4 1
0.3 0.5
4
3
0.2 0
0.1 −0.5
0 −1
−0.1 −1.5
−0.2 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.26 Derivative-free Kalman Filtering: a tracking of a seesaw position setpoint (dashed line)
by the angle of the second joint of the robot (continuous line), b tracking of a seesaw velocity
setpoint (dashed line) by the angular velocity of the second joint of the robot (continuous line)
position and velocity setpoints is depicted in Fig. 4.29. The control inputs (motor
torques) applied to the first and second joint of the robotic manipulator are shown in
Fig. 4.30.
The performance of the state estimation-based control, with the use of the
Extended Kalman Filter is also shown in the tracking of a seesaw setpoint. The
position and velocity variations of the first joint are demonstrated in Fig. 4.31. Simi-
larly, the tracking of the position and velocity reference setpoints for the second joint
4.2 The Derivative-Free Nonlinear Kalman Filter 169
(a) 60 (b) 30
55
25
50
control signal of joint 1
40 15
35
10
30
5
25
20 0
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.27 Tracking of a seesaw setpoint by the 2-DOF robotic manipulator with the use of the
derivative-free Kalman Filter: a control input (motor torque) applied to the first joint, b control
input (motor torque) applied to the second joint
1.5
1
state variable x of the robot
state variable x of the robot
1
0.5
0.5
2
1
0 0
−0.5
−0.5
−1
−1
−1.5
−1.5 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.28 Extended Kalman Filtering: a tracking of a sinusoidal position setpoint (dashed line) by
the angle of the first joint of the robot (continuous line), b tracking of a sinusoidal velocity setpoint
(dashed line) by the angular velocity of the first joint of the robot (continuous line)
are depicted in Fig. 4.32. The control signal in the case of tracking a seesaw setpoint
by the two joints of the robotic manipulator is shown in Fig. 4.33.
Comparing the estimation performed by the derivative-free MIMO nonlinear
Kalman Filter with the one performed by the Extended Kalman Filter it can be
noticed that the derivative-free filtering approach results in more accurate state
170 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
1.5
1
state variable x3 of the robot
0 0
−0.5
−0.5
−1
−1
−1.5
−1.5 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.29 Extended Kalman Filtering: a tracking of a sinusoidal position setpoint (dashed line)
by the angle of the second joint of the robot (continuous line), b tracking of a sinusoidal velocity
setpoint (dashed line) by the angular velocity of the second joint of the robot (continuous line)
(a) 60 (b) 30
59 29
58 28
control signal of joint 1
57 27
56 26
55 25
54 24
53 23
52 22
51 21
50 20
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.30 Tracking of a sinusoidal setpoint by the 2-DOF robotic manipulator with the use of the
Extended Kalman Filter: a control input (motor torque) applied to the first joint, b control input
(motor torque) applied to the second joint
estimates. The precision of the provided estimation is the one depicted in Tables 4.1
and 4.2 of Sect. 4.2.4 for the case of the DC motor. Moreover, comparing the associ-
ated state estimation-based control loop that was based on the derivative-free MIMO
nonlinear Kalman Filter to the control that was based on the Extended Kalman
4.2 The Derivative-Free Nonlinear Kalman Filter 171
0.6 1.5
state variable x of the robot
2
0
0.3
−0.5
0.2
−1
0.1 −1.5
0 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.31 Extended Kalman Filtering: a tracking of a seesaw position setpoint (dashed line) by
the angle of the first joint of the robot (continuous line), b tracking of a seesaw velocity setpoint
(dashed line) by the angular velocity of the first joint of the robot (continuous line)
0.6 1.5
state variable x of the robot
1
0.5
0.5
0.4
3
0
0.3
−0.5
0.2
−1
0.1 −1.5
0 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.32 Extended Kalman Filtering: a tracking of a seesaw position setpoint (dashed line) by
the angle of the second joint of the robot (continuous line), b tracking of a seesaw velocity setpoint
(dashed line) by the angular velocity of the second joint of the robot (continuous line)
Filter it was observed that the first control scheme was significantly more robust and
capable of tracking with better accuracy the desirable trajectories. These findings
show the suitability of the considered derivative-free MIMO nonlinear Kalman Fil-
ter for tracking, control, and fault diagnosis tasks.
172 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
(a) 60 (b) 35
50 30
control signal of joint 1
20
30
15
20
10
10
5
0 0
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.33 Tracking of a seesaw setpoint by the 2-DOF robotic manipulator with the Extended
Kalman Filter: a control input (motor torque) applied to the first joint, b control input (motor
torque) applied to the second joint
4.3.1 Overview
This part of the chapter extends the use of the Derivative-free nonlinear Kalman
Filter toward distributed state estimation-based control of nonlinear MIMO systems.
Up to now an established approach for performing distributed state estimation has
been the Extended Information Fitter [427]. In the Extended Information Filter the
local filters do not exchange raw measurements but send to an aggregation filter
their local information matrices (local inverse covariance matrices which can be also
associated to Fisher Information Matrices) and their associated local information
state vectors (products of the local information matrices with the local state vectors)
[264, 414]. The Extended Information Filter performs fusion of state estimates from
local distributed Extended Kalman Filters which in turn are based on the assump-
tion of linearization of the nonlinear system dynamics by first-order Taylor series
expansion and truncation of the higher order linearization terms [329, 369, 402].
Moreover, the Extended Kalman Filter requires the computation of Jacobians which
in the case of high order nonlinear dynamical systems can be a cumbersome proce-
dure. This approach introduces cumulative errors to the state estimation performed
by the local Extended Kalman Filter recursion which is finally transferred to the mas-
ter filter where the aggregate state estimate of the controlled system is computed.
Consequently, these local estimation errors may result in the deterioration of the
performance of the associated control loop or even risk its stability [415, 427].
4.3 The Derivative-Free Distributed Nonlinear Kalman Filter 173
Again, the discrete-time nonlinear system of Eqs. (4.46) and (4.47) is considered. The
Extended Information Filter (EIF) performs fusion of the local state vector estimates
which are provided by local Extended Kalman Filters, using the Information matrix
and the Information state vector [264]. The Information Matrix is the inverse of the
state vector covariance matrix, and can be also associated to the Fisher Information
matrix [414]. The Information state vector is the product between the Information
matrix and the local state vector estimate
174 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
The update equation for the Information Matrix and the Information state vector are
given by
where I (k) = JγT (k)R ( k)−1 Jγ (k) is the associated information matrix and, i(k) =
JγT (k)R ( k)−1 [(z(k) − γ (x(k))) + Jγ x̂ − (k)] is the information state contribution.
The predicted information state vector and Information matrix are obtained from
It is assumed that an observation vector z i (k) is available for the N different sensor
sites (vision nodes) i = 1, 2, . . . , N and each vision node observes the robot accord-
ing to the local observation model, expressed by z i (k) = γ (x(k)) + vi (k), i =
1, 2, . . . , N , where the local noise vector vi (k)∼N (0, R i ) is assumed to be white
Gaussian and uncorrelated between sensors. The variance of a composite observa-
tion noise vector vk is expressed in terms of the block diagonal matrix R(k) =
diag[R(k)1 , . . . , R N (k)]T . The information contribution can be expressed by a lin-
ear combination of each local information state contribution i i and the associated
information matrix I i at the ith sensor site
N iT −1 i −
i(k) = i=1 Jγ (k)R (k) [z (k) − γ (x(k)) + Jγ (k) x̂ (k)]
i i i
N i T −1 i
(4.55)
I (k) = i=1 Jγ (k)R (k) Jγ (k).
i
Thus, the update equations for fusing the local state estimates is
N
ŷ(k) = ŷ− (k) + iT −1 i −
i=1 Jγ (k)R (k) [z (k) − γ (x(k)) + Jγ (k) x̂ (k)]
i i i
N i T (4.56)
Y(k) = Y− (k) + i=1 Jγ (k)R i (k)−1 Jγi (k)
−1
Y − (k + 1) = P − (k + 1) = [Jφ (k)P(k)Jφ (k)T + Q(k)]−1
−1 −
(4.58)
and y − (k + 1) = P − (k + 1) x̂ (k + 1)
If the derivative-free nonlinear Kalman Filter is used in place of the Extended Kalman
Filter then in the EIF equations the following matrix substitutions should be per-
formed: Jφ (k)→Ad , Jγ (k)→Cd , where matrices Ad and Cd are the discrete-time
equivalents of matrices Ac and Cc which have been defined in Eq. (4.49). Matrices
Ad and Cd can be computed using established discretization methods, as already
explained in Chap. 3. Moreover, the covariance matrices P(k) and P − (k) are the
ones obtained from the linear Kalman Filter update equations given in Sect. 4.2.2.
The outputs of the local filters are treated as measurements which are fed into the
aggregation fusion filter (see Fig. 4.34) [264]. Then each local filter is expressed by
its respective error covariance and estimate in terms of information contributions and
is described by
−1
Pi −1 (k) = Pi− (k) + JγT (k)R(k)−1 Jγ (k)x̂i (k) =
= Pi (k)(Pi− (k)−1 x̂i− (k)) + JγT (k)R(k)−1 [z i (k) − γ i (x(k)) + Jγi (k)x̂i− (k)]
(4.59)
The global estimate and the associated error covariance for the aggregate fusion filter
can be rewritten in terms of the computed estimates and covariances from the local
filters using the relations
176 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
Fig. 4.34 Fusion of the distributed state estimates with the use of the Extended Information Filter
For the general case of N local filters i = 1, . . . , N , the distributed filtering archi-
tecture is described by the following equations
N −1 − P − (k)−1 ]
P(k)−1 = P − (k)−1 + i=1 [Pi (k) i
N (4.61)
x̂(k) = P(k)[P − (k)−1 x̂ − (k) + i=1 (Pi (k)−1 x̂i (k) − Pi− (k)−1 x̂i− (k))]
The global state update equation in the above distributed filter can be written in
terms of
the information state vector and of the information
N matrix, i.e., ŷ(k) =
ŷ − (k) + i=1
N
( ŷi (k) − ŷi− (k)), and Ŷ (k) = Ŷ − (k) + i=1 (Ŷi (k) − Ŷi− (k)). From
Eq. (4.61) it can be seen that if a local filter (processing station) fails, then the local
covariance matrices and the local state estimates provided by the rest of the filters
will enable an accurate computation of the monitored system’s state vector.
4.3 The Derivative-Free Distributed Nonlinear Kalman Filter 177
As mentioned above, for the system of Eqs. (4.46) and (4.47), state estimation is
possible by applying the standard Kalman Filter. The system is first turned into
discrete-time form using common discretization methods and then the recursion of
the linear Kalman Filter described in Eqs. (4.5) and (4.6) is applied.
If the derivative-free Kalman Filter is used in place of the Extended Kalman Filter
then in the EIF equations the following matrix substitutions should be performed:
Jφ (k)→Ad , Jγ (k)→Cd , where matrices Ad and Cd are the discrete-time equivalents
of matrices A and C which have been defined Eq. (4.49) and which appear also in
the measurement and time update of the standard Kalman Filter recursion. Matrices
Ad and Cd can be computed using established discretization methods. Moreover, the
covariance matrices P(k) and P − (k) are the ones obtained from the linear Kalman
Filter update equations given in Sect. 4.2.2.
The example to be presented describes the control of a planar robot with the use of
a position-based visual servo that comprises multiple fixed cameras. The chapter’s
approach relies on neither position nor velocity sensors, and directly sets the motor
control current using only visual feedback. Direct visual servoing is implemented
using a distributed filtering scheme which permits to fuse the estimates of the robot’s
state vector computed by local filters, each one associated to a camera in the cameras
network (see Fig. 4.35). The cameras’ network can be based on multiple cameras
connected to a computer with a frame grabber to form a vision node (see Figs. 4.35
and 4.36).
The master computer communicates video synchronization information over the
network to each vision node. Typical sources of measurement noise include charge-
coupled device (CCD) noise, analog-to-digital (A/D) noise, and finite word-length
effects. Under ideal conditions, the effective noise variance from these sources should
remain relatively constant. Occlusions can be also considered as a noise source.
Finally, communication delays and packet drops in the transmission of measure-
ments from the vision sensors to the information processing nodes induce additional
disturbances which should be compensated by the virtual servoing control loop.
178 4 Nonlinear Kalman Filtering Based on Differential Flatness Theory
Fig. 4.35 Distributed cameras network and distributed information processing units for visual
servoing
Fusion of the local state estimates which are provided by filters running on the vision
nodes can improve the accuracy and robustness of the performed state estimation,
thus also improving the performance of the robot’s control loop [126, 478, 506,
507, 594]. Under the assumption of Gaussian noise, a possible approach for fusing
the state estimates from the distributed local filters is the derivative-free distributed
nonlinear Kalman Filter. As explained in Sect. 4.3.2, this filter (that is also named
derivative-free Extended Information Filter) provides an aggregate state estimate by
weighting the state vectors produced by local Kalman Filters with the inverse of the
associated estimation error covariance matrices.
Visual servoing over the previously described cameras network is considered for
the nonlinear dynamic model of a two-link robotic manipulator. The robot can be
programmed to execute a manufacturing task, such as disassembly or welding [532].
The position of the robot’s end effector in the cartesian space (and consequently
the angle for the robotic link) is measured by the aforementioned m distributed
4.3 The Derivative-Free Distributed Nonlinear Kalman Filter 179
Fig. 4.36 Visual servoing based on fusion of state estimates provided by local derivative-free
nonlinear Kalman Filters
cameras. The proposed multi-camera based robotic control loop can be also useful
in other vision-based industrial robotic applications where the vision is occluded or
heavily disturbed by noise sources, e.g., cutting. In such applications there is need to
fuse measurements from multiple cameras so as to obtain redundancy in the visual
information and permit the robot to complete safely and within the specified accuracy
constraints its assigned tasks [365, 594].
1.5
1
state variable x1 of the robot
0 0
−0.5
−0.5
−1
−1
−1.5
−1.5 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.37 Derivative-free Kalman Filtering: a tracking of a sinusoidal position setpoint (dashed
line) by the angle of the first joint of the robot (continuous line), b tracking of a sinusoidal velocity
setpoint (dashed line) by the angular velocity of the first joint of the robot (continuous line)
1.5
1
state variable x of the robot
1
0.5
0.5
3
0 0
−0.5
−0.5
−1
−1
−1.5
−1.5 −2
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.38 Derivative-free Kalman Filtering: a tracking of a sinusoidal position setpoint (dashed
line) by the angle of the second joint of the robot (continuous line), b tracking of a sinusoidal
velocity setpoint (dashed line) by the angular velocity of the second joint of the robot (continuous
line)
4.3 The Derivative-Free Distributed Nonlinear Kalman Filter 181
(a) 60 (b) 26
24
55
22
control signal of joint 1
45 18
16
40
14
35
12
30 10
0 5 10 15 20 0 5 10 15 20
t (sec) t (sec)
Fig. 4.39 Tracking of a sinusoidal setpoint by the 2-DOF robotic manipulator with the use of the
derivative-free Kalman Filter: a control input (motor torque) applied to the first joint, b control
input (motor torque) applied to the second joint
The position and velocity estimates for the first joint of the robot when tracking a
sinusoidal setpoint are depicted in Fig. 4.37. For the second joint of the 2-DOF robot,
and in case of tracking again a sinusoidal setpoint, the position and angular velocity
estimates are depicted in Fig. 4.38. The control inputs (motor torques) applied to the
first and second joint of the robotic manipulator are shown in Fig. 4.39.
Chapter 5
Differential Flatness Theory
and Industrial Robotics
5.1 Overview
elements which depend on the system’s parameters. The nonlinear terms which
appear in the control inputs of the robot are approximated with the use of neuro-fuzzy
networks. Moreover, since only the system’s output is measurable the complete state
vector has to be reconstructed with the use of a state observer. It is shown that a suit-
able learning law can be defined for the aforementioned neuro-fuzzy approximators
so as to preserve the closed-loop system stability. With the use of Lyapunov stability
analysis, it is proven that the proposed observer-based adaptive neuro-fuzzy control
scheme results in H-inifnity tracking performance.
The chapter also shows that apart from adaptive control for the robotic manipula-
tors, it is possible to develop Kalman Filter-based control with the use of differential
flatness theory. The Derivative-free nonlinear Kalman Filter is used for developing
a robust controller which can be applied to underactuated MIMO robotic systems.
The control problem for underactuated robots is nontrivial and becomes further com-
plicated if the robot are subjected to model uncertainties and external disturbances.
Using differential flatness theory, it is shown that the model of a closed-chain 2-DOF
robotic manipulator can be transformed to linear canonical form. For the linearized
equivalent of the robotic system, it is shown that a state feedback controller can be
designed. Since certain elements of the state vector of the linearized system cannot
be measured directly, it is proposed to estimate them with the use of the differential
flatness theory-based method for state estimation (analyzed in Chap. 4) which is the
Derivative-free nonlinear Kalman Filter. Moreover, by redesigning the Kalman Filter
as a disturbance observer, it is shown that one can estimate simultaneously exter-
nal disturbances terms that affect the robotic model or disturbance terms which are
associated with parametric uncertainty. The efficiency of the proposed Kalman Filter-
based control scheme is checked in the case of a 2-DOF planar robotic manipulator
that has the structure of a closed-chain mechanism.
Another section of the chapter analyzes a distributed filtering scheme for robotic
visual servoing capable of compensating for modeling uncertainties and external
disturbances. The visual servoing robotic scheme that was introduced in Chap. 4 is
revisited. This time, the following disturbances are assumed to affect the control
loop: (i) Contact forces at the manipulator’s end effector, (ii) missing measurements
about the robot’s motion characteristics (the measurements transmission link used for
network-based control exhibits delays or data losses). First, it is assumed that during
the compliance task performed by the robot unknown or varying forces are exerted
on its end effector. The derivative-free Extended Information Filter (derivative-free
distributed nonlinear Kalman Filter) is suitably modified through the use of concepts
from disturbance observers theory to compensate for the effect of these perturbations.
Subsequently, the chapter examines the use of the Derivative-free nonlinear
Kalman Filter in the compensation of delayed measurements and data losses in
robotic visual servoing. An analysis is given about the problem of distributed non-
linear filtering over a communication/sensors network, and the use of the estimated
state vector in a control loop. Distributed filtering is now based on a derivative-free
implementation of nonlinear Kalman Filtering. It is shown how the proposed distrib-
uted filtering method can succeed compensation of random delays and packet drops
5.1 Overview 185
which may appear during the transmission of measurements and of state vector esti-
mates, thus assuring a reliable performance of the distributed filtering-based control
scheme. Evaluation tests confirm the improved performance of the Derivative-free
distributed nonlinear Kalman Filter against the Extended Information Filter.
5.2.1 Overview
differentially flat model. This condition holds for underactuated robotic systems,
thus providing a systematic approach to the design of reliable controllers for such
systems [426, 433].
The considered closed-chain 2-DOF robotic system depicted in Fig. 5.1 consists of
four bodies: (i) Bodies 1 and 2 are two sliders with masses m1 and m2 , respectively
[159, 597]. Body 3 is connected with a revolute joint to body 1 and has mass m3 ,
length l3 while its moment of inertia is I3 . Similarly body 4 is connected to body 2
with a revolute joint, has mass m4 , length l4 while the associated moment of inertia
is I4 . The motion of the system takes place in the 2D xy plane depicted in Fig. 5.1
while its dynamics is subjected to gravity.
The state variables for the robotic system are as follows q = [q1 , q2 , q3 , q4 ]T : q1
is the displacement of mass m1 along the x-axis, q3 is the turn angle of body 3 round
the revolute joint A. q2 is the displacement of m2 along the x-axis and q4 is the turn
angle of body 4 round the revolute joint B. The following geometric constraints hold:
l3 sin(q3 ) = l4 sin(q4 )
(5.1)
q1 + l3 cos(q3 ) = q2 + l4 cos(q4 )
The control inputs exerted on the robotic model are the horizontal force F1 that
causes the displacement of mass m1 along the x-axis, the torque T3 that causes the
rotation of the link with length l3 round the revolute joint A, the horizontal force
F2 that causes the displacement of mass m2 along the x-axis and the torque T4 that
causes the rotation of the link with length l4 round the revolute joint B. Thus, in the
most generic case the input vector can be u = [ f1 , T3 , f2 , T4 ]T . According the Euler-
Lagrange analysis, the dynamic model of the robotic manipulator is obtained [597]
where ⎛ ⎞
a11 a12 0 0
⎜a21 a22 0 0 ⎟
A(q) = ⎜
⎝ 0
⎟ (5.3)
0 a23 a24 ⎠
0 0 a43 a44
with a11 = m1 + m3 , a12 = a21 = −m3 lc3 sin(q3 ), a22 = m3 lc23 + I3 , a33 = m2 + m4 ,
a34 = a43 = −m4 lc4 sin(q4 ), a44 = m4 lc24 + I4 , and
⎛ ⎞
−m3 lc3 cos(q3 )q̇32
⎜ m3 glc cos(q3 ) ⎟
h(q, q̇) = ⎜ 3
⎝−m4 lc4 cos(q4 )q̇2 ⎠
⎟ (5.4)
4
m4 glc4 cos(q4 )
The closed-chain robotic system of Fig. 5.1, presented initially in Sect. 5.2.2 is con-
sidered again. This comprises four masses, out of which mass m1 and mass m2 slide
along the x-axis, while masses m3 and m4 associated with links l3 and l4 perform
both translational and rotational motion in the xy-plane.
For mass m1 one has that the potential energy is P1 = 0 while the kinetic energy
is K1 = 21 m1 q̇12 . For mass m2 one has again that the potential energy is P2 = 0 while
the kinetic energy is K2 = 21 m2 q̇22 .
For mass m3 the potential energy is given by
The kinetic energy of mass m3 (link l3 ) is due to both translational motion and
rotational motion. Denoting by I3 the moment of inertia for the rotation of the link
l3 , one has the kinetic energy that is associated with the rotational motion is K3,1 =
1 2
2 I3 q̇3 . The kinetic energy of link l3 that is associated with the translational motion
is denoted as K3,2 and is given by
Ċ3 is the velocity of the center of gravity of mass (link) 3 in cartesian coordinates.
Using that the position of the center of gravity of link 3 is the vector
Thus, using Eq. (5.8) the kinetic energy of link 3 that is due to its translational motion
is given by
Using K3,1 and K3,2 one has that the aggregate kinetic energy of mass (link) m3 is
The kinetic energy of mass m4 (link l4 ) is due to both translational motion and
rotational motion. Denoting by I4 the moment of inertia for the rotation of the link
l4 , one has the kinetic energy that is associated with the rotational motion is K4,1 =
1 2
2 I4 q̇4 . The kinetic energy of link l4 that is associated with the translational motion
is denoted as K4,2 and is given by
Ċ4 is the velocity of the center of gravity of mass (link) 4 in cartesian coordinates.
Using that the position of the center of gravity of link 4 is the vector
Thus, using Eq. (5.14) the kinetic energy of link 4 that is due to its translational
motion is given by
5.2 Adaptive Fuzzy Control of Underactuated MIMO Robots 189
Using K4,1 and K4,2 one has that the aggregate kinetic energy of mass (link) m4 is
L = 21 m1 q̇12 + 21 m2 q̇22 +
+ 21 I3 q̇32 + 21 m3 [q̇12 + lc23 q̇32 − 2lc3 sin(q3 )q̇1 q̇3 ]+ (5.17)
+ 21 I4 q̇42 + 2 m4 [q̇2
1 2 + lc24 q̇42 − 2lc4 sin(q4 )q̇2 q̇4 ]+
−m3 glc3 sin(q3 ) − m4 glc4 sin(q4 )
∂ ∂L
∂t ∂ q̇3 = I3 q̈3 + m3 lc23 q̈3 − m3 lc3 cos(q3 )q̇1 q̇3 − m3 lc3 sin(q3 )q̈1 (5.25)
Next, by using the previous equations and by applying the Euler-Lagrange principle
one gets
∂ ∂L ∂L
∂t ∂ q̇1 − ∂q1 = F1 ⇒ (5.30)
(m1 + m3 )q̈1 − m3 lc3 sin(q3)q̈3 − m3 lc3 cos(q3 )q̇32 = F1
∂ ∂L ∂L
− ∂q
∂t ∂ q̇2 = F2 ⇒
2 (5.31)
(m2 + m4 )q̈2 − m4 lc4 sin(q4)q̈4 − m4 lc4 cos(q4 )q̇42 = F2
∂ ∂L ∂L
− ∂q
∂t ∂ q̇3 3
= T1 ⇒
−m3 lc3 sin(q3 )q̈1 + (I3 + m3 lc23 )q̈3 − m3 lc3 cos(q3 )q̇1 q̇3 +
+m3 lc3 cos(q3 )q̇1 q̇3 + m3 glc3 cos(q3 ) = T1 ⇒ (5.32)
∂ ∂L ∂L
− ∂q
∂t ∂ q̇4 4
= T2 ⇒
−m4 lc4 sin(q4 )q̈2 + (I4 + m4 lc24 )q̈4 − m4 lc4 cos(q4 )q̇2 q̇4 +
+m4 lc4 cos(q4 )q̇2 q̇4 + m4 glc4 cos(q4 ) = T2 ⇒ (5.33)
4
m4 glc4 cos(q4 ) T2
(5.36)
Equation (5.36) can also take the compact matrix form
This is the description of the robotic mechanism which has been given in Eqs. (5.2),
(5.3) and (5.4) .
Next, the case in which l3 = l4 is examined. Moreover, it is considered that the mass
m2 is connected to a spring with elasticity k2 . Finally, it is assumed that the only
inputs applied to the robotic model are u1 = f1 and u2 = T3 . Then the dynamic
model of the robot becomes [597]
M1 + M2 −2M2 l3 sin(q3 )
A(ql ) = (5.39)
−2M2 l3 sin(q3 ) I3 + I4 + 4M2 l32 sin2 (q3 )
where
⎛ ⎞
q̇1
⎜ q̇3 ⎟
⎜ −k l (I +I )+2l M (k (π −q )sin(q )+(I +I )q̇2 cos(q )) ⎟
⎜ 2d 3 4 3 2 4 3 3 3 4 3 3 ⎟
f (x) = ⎜
⎜ M2 (I3 +I4 )+M1 (I3 +I4 +4M2 l32 sin2 (q3 )) ⎟
⎟ (5.42)
⎜ ⎟
⎝ ⎠
k4 (M1 +M2 )(π −q3 )+2I3 M1 sin(q3 )(k2 ld −2l3 M2 q̇3 cos(q3 ))
2
⎛ ⎞
0
⎜ 0 ⎟
⎜ ⎟
⎜ I3 +I4 +4M2 l32 sin2 (q3 ) ⎟
g1 (x) = ⎜ ⎟
⎜ (I3 +I4 )M2 +M1 (I3 +I4 +4M2 l3 sin (q3 )) ⎟
2 2 (5.43)
⎜ ⎟
⎝ ⎠
2M2 l3 sin(q3 )
(I3 +I4 )M2 +M1 (I3 +I4 +4M2 l32 sin2 (q3 ))
⎛ ⎞
0
⎜ 0 ⎟
⎜ 2M2 l3 sin(q3 ) ⎟
⎜ ⎟
g2 (x) = ⎜ (I3 +I4 )M2 +M1 (I3 +I4 +4M2 l32 sin2 (q3 )) ⎟ (5.44)
⎜ ⎟
⎝ ⎠
M1 +M2
(I3 +I4 )M2 +M1 (I3 +I4 +4M2 l32 sin2 (q3 ))
The robotic system is underactuated when only one of the control inputs is enabled.
It can be proven that when the only input is u1 = F1 then the robotic system is not
static feedback linearizable.
Next, the case in which the only control input is u2 = T3 is examined. When
k2 = 0 and k4 = 0 then the robotic model admits static feedback linearization.
Equivalently, with the use of the feedback linearization criteria stated in Chap. 1
this means that (i) the distribution D3 = [g(x), adf (x) g(x), adf2(x) g(x), adf3(x) g(x)]
has full rank and (ii) the vector fields D0 = [g(x)], D1 = [g(x), adf (x) g(x)] and
D2 = [g(x), adf (x) g(x), adf2(x) g(x)] are involutive.
After having proven that the closed-chain robotic manipulator admits static feedback
linearization, its transformation to a linear form will be attempted using suitably
chosen feedback linearizing outputs. The following linearizing output is defined first
Similarly, it holds
⎛ ⎞
f1
⎜f ⎟
∂z2 ∂z2 ∂z2 ∂z2 ⎜ 2 ⎟
z3 = Lf2 h1 = ∂q1 ∂q3 ∂ q̇1 ∂ q̇3 ⎝f ⎠ ⇒
3
f4 ⎛ ⎞
f1 (5.47)
⎜ f2 ⎟
⎜
z3 = Lf h1 = 0 −2M2 l3 cos(q3 )q̇3 (M1 + M2 ) −2M2 l3 sin(q3 ) ⎝ ⎠ ⇒
2 ⎟
f3
f4
z3 = Lf2 h1 = (M1 + M2 )f3 − 2M2 l3 sin(q3 )f4 − 2M2 l3 cos(q3 )q̇3 f2
It holds that
Therefore, it holds
or equivalently
z3 = (M1 + M2 )q̈1 − 2M2 l3 sin(q3 )q̈3
(5.50)
(M1 +M2 )2M2 l3 sin(q3 )
+ (I u2
3 +I4 )M2 +M1 (I3 +I4 +4M2 l3 sin (q3 ))
2 2
Therefore,
⎛ ⎞
f1
f ⎟
∂z4 ∂z4 ∂z4 ∂z4 ⎜
Lf h1 = ∂q1 ∂q3 ∂ q̇1 ∂ q̇3 ⎜
4 2⎟
⎝ f3 ⎠ ⇒
f4 ⎛ ⎞
f1
⎜f2 ⎟ (5.56)
Lf h1 = 0 2k2 l3 cos(q3 )q̇3 −k2 2k2 l3 sin(q3 ) ⎜
4 ⎟
⎝ f3 ⎠ ⇒
f4
Lf4 h1 = 2k2 l3 cos(q3 )q̇3 f2 − k2 f3 + 2k2 l3 sin(q3 )f4 ⇒
Lf4 h1 = 2k2 l3 cos(q3 )q̇32 − k2 f3 + 2k2 l3 sin(q3 )f4
5.2 Adaptive Fuzzy Control of Underactuated MIMO Robots 195
and that ż1 = z2 , ż2 = z3 , ż3 = z4 one has that the robotic model can be written in
the linear canonical (Brunovsky) form
⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
ż1 0 1 0 0 z1 0
⎜ż2 ⎟ ⎜0 0 1 0⎟ ⎜ z2 ⎟ ⎜ 0 ⎟
⎜ ⎟=⎜ ⎟⎜ ⎟ + ⎜ ⎟v (5.60)
⎝ż3 ⎠ ⎝0 0 0 1 ⎠ ⎝ z3 ⎠ ⎝ 0 ⎠
ż4 0 0 0 0 z4 1
v = ż4d − k1 (z4 − z4d ) − k2 (z3 − z3d ) − k3 (z2 − z2d ) − k4 (z1 − z1d ) (5.61)
Next, a linearized equivalent description for the robot’s kinematic chain will be
obtained with the use of differential flatness theory. The following flat output is
chosen
y = (M1 + M2 )q1 + 2M2 l3 cos(q3 ) (5.62)
196 5 Differential Flatness Theory and Industrial Robotics
It holds that
Consequently
y(3) = −k2 [q̇1 − 2l3 sin(q3 )q̇3 ] (5.66)
and, respectively
y(4) = −k2 q̈1 + 2k2 l3 cos(q3 )q̇32 + 2k2 l3 sin(q3 )q̈3 (5.67)
From Eq. (12.112) describing the flat output and from the equations of its higher
order derivatives one has a set of equations which can be solved with respect to the
state variables q1 , q3 , q̇1 and q̇3 . It holds that
k2 y−M2 ÿ+M2 k2 L
q1 = k 2 M1 (5.68)
2 q1 −k2 L
q3 = cos−1 ( ÿ+k2k 2 l3
) (5.69)
k2 ẏ+M2 y(3)
q̇1 = k 2 M1
(5.70)
Having expressed the elements of the state vector q as functions of the flat output and
its derivatives and knowing that q̈1 = f3 +ga3 u1 +gb3 u2 and q̈3 = f4 +ga4 u1 +gb4 u2 ,
it holds that the control inputs u1 and u2 can be also written as functions of the flat
output and its derivatives. Therefore, the robotic system stands for a differentially
flat model.
From the relation q̇ = f (q) + ga (q)u1 + gb (q)u2 and the associated relations about
f (q), ga q and gb (q) it holds
5.2 Adaptive Fuzzy Control of Underactuated MIMO Robots 197
−k2 ld (I3 +I4 )+2l3 M2 [k4 (π −q3 )sin(q3 )+(I3 +I4 )q̇32 cos(q3 )]
q̈1 = M2 (I3 +I4 )+M1 (I3 +I4 +4M2 l32 sin2 (q3 ))
+ (5.72)
2M2 l3 sin(q3 )
+ (I +I )M +M (I +I +4M l2 sin2 (q )) u2
3 4 2 1 3 4 2 3 3
k4 (M1 +M2 )(π −q3 )+2l3 M1 sin(q3 )[K2 ld −2l3 M2 q̇32 cos(q3 )]
q̈3 = M2 (I3 +I4 )+M1 (I3 +I4 +4M2 l32 sin2 (q3 ))
+ (5.73)
1 +M2
+ (I +I )M +M M u2
1 (I3 +I4 +4M2 l3 sin (q3 ))
2 2
3 4 2
Consequently, it holds
Considering that the linear displacement of the left part of the kinematic chain q1 is
measurable, and that the same holds for the turn angle q3 of joint A one has that the
flat output y = (M1 + M2 )q1 + 2M2 l3 cos(q3 ) is also a measurable variable.
Using the description of system in the linear canonical form, the appropriate
control law is
v = ż4d − k1 (z4 − z4d ) − k2 (z3 − z3d ) − k3 (z2 − z2d ) − k4 (z1 − z1d ) (5.77)
198 5 Differential Flatness Theory and Industrial Robotics
Moreover, using that the transformed control input that is applied to the robot is
v = fv + gv u2 , one can compute the real control input u2 = gv−1 (v − fv ).
After transforming the dynamic model of the robotic mechanism into the linear
canonical form, the design of a state feedback controller becomes possible. To cope
with the lack of knowledge about the model’s parameters adaptive fuzzy control will
be implemented. The design of the adaptive fuzzy control scheme for the robotic
model that is transformed to the canonical (Brunovsky) form has been explained
in Chap. 3. Moreover, the stages of the associated Lyapunov stability analysis have
been given, and it has been explained that the proposed adaptive fuzzy control method
succeeds H∞ tracking performance for the closed-loop.
The unknown dynamics of the robot is modeled with the use of neuro-fuzzy
approximators. These have as inputs the state variables of the robotic model xi , i =
1, . . . , 4. Thus one arrives at the approximation of the unknown dynamics of the
robot fˆ (x) and ĝ(x), which in turn are used for the computation of the control signal
u.
The efficiency of the proposed adaptive fuzzy control scheme for the underactuated
robotic manipulator was evaluated in the tracking of various setpoints for state vari-
able q1 i.e., the linear displacement of mass M1 of the mechanism and state variable
q3 , i.e., the rotation angle of joint A. As it can be observed in Figs. 5.2, 5.3 and 5.4,
the proposed control scheme enabled accurate tracking of the reference setpoints
by the robot’s state variables. Indicative results about the estimation of unknown
functions fv (x, t) and gv (x, t) in the robot’s dynamics by the neuro-fuzzy approxi-
mators is shown in Fig. 5.5. By including an additional term in the control loop that
was based on the disturbances estimation, it became possible to compensate for the
disturbances effects.
5.3 Observer-Based Adaptive Fuzzy Control of MIMO Robots 199
(a) 20
x1 − xd1 3 (b) 20 2.5
x2 − xd2
x1 − xd1
x2 − xd2
2.5
10 10 2
2
0 0 1.5
1.5
−10 1 −10 1
0 20 40 60 0 20 40 60 0 20 40 60 0 20 40 60
time (sec) time (sec) time (sec) time (sec)
6 0.6 6 0.3
x3 − xd3
x4 − xd4
x3 − xd3
x4 − xd4
4 0.4 4 0.2
2 0.2 2 0.1
0 0 0 0
−2 −0.2 −2 −0.1
0 20 40 60 0 20 40 60 0 20 40 60 0 20 40 60
time (sec) time (sec) time (sec) time (sec)
x2 − xd2
x2 − xd2
10 2.5
2 10
5 2
1.5 0
0 1.5
−5 1 −10 1
0 20 40 60 0 20 40 60 0 20 40 60 0 20 40 60
time (sec) time (sec) time (sec) time (sec)
6 0.3 6 0.6
x3 − xd3
x4 − xd4
x3 − xd3
x4 − xd4
4 0.2 4 0.4
2 0.1 2 0.2
0 0 0 0
−2 −0.1 −2 −0.2
0 20 40 60 0 20 40 60 0 20 40 60 0 20 40 60
time (sec) time (sec) time (sec) time (sec)
5.3.1 Overview
(a) 20
x1 − xd1 3 (b) 30 3
x2 − xd2
x1 − xd1
x2 − xd2
2.5 20 2.5
10
2 10 2
0
1.5 0 1.5
−10 1 −10 1
0 20 40 60 0 20 40 60 0 20 40 60 0 20 40 60
time (sec) time (sec) time (sec) time (sec)
6 x4 − xd4 0.6 6 0.6
x3 − xd3
x4 − xd4
x3 − xd3
4 0.4 4 0.4
2 0.2 2 0.2
0 0 0 0
−2 −0.2 −2 −0.2
0 20 40 60 0 20 40 60 0 20 40 60 0 20 40 60
time (sec) time (sec) time (sec) time (sec)
(a) 100 50
80 45
real vs estimated function g
real vs estimated function f
60 40
40 35
20 30
0 25
−20 20
−40 15
−60 10
−80 5
−100 0
0 10 20 30 40 50 60 0 10 20 30 40 50 60
time (sec) time (sec)
(b) 100 50
80 45
real vs estimated function g
real vs estimated function f
60 40
40 35
20 30
0 25
−20 20
−40 15
−60 10
−80 5
−100 0
0 10 20 30 40 50 60 0 10 20 30 40 50 60
time (sec) time (sec)
Fig. 5.5 Estimation of functions fv (x, t) and gv (x, t) of the robot dynamics by neuro-fuzzy networks
when tracking: a setpoint 1, b setpoint 2
5.3 Observer-Based Adaptive Fuzzy Control of MIMO Robots 201
errors and the external disturbances on the tracking error is attenuated to an arbitrary
desirable level. Differential flatness theory is used to transform the MIMO robotic
system into the canonical form. The resulting control inputs are shown to contain
nonlinear elements which depend on the system’s parameters. Since the parameters
of the system are unknown, then the nonlinear terms which appear in the control
inputs have to be approximated with the use of neuro-fuzzy networks. Moreover,
since only the system’s output is measurable the complete state vector has to be
reconstructed with the use of a state observer. In the current section, it is shown that
a suitable learning law can be defined for the aforementioned neuro-fuzzy approxi-
mators so as to preserve the closed-loop system stability. Lyapunov stability analysis
also proves that the proposed observer-based adaptive fuzzy control scheme results
in H∞ tracking performance, in accordance to the results of [282, 407, 410, 413].
For the design of the observer-based adaptive fuzzy controller, one has to solve
two Riccati equations, where the first one is associated with the controller and the
second one is associated with the observer. Parameters that affect the closed-loop
robustness are: (i) The feedback gain vector K, (ii) the observer’s gain vector Ko ,
and (iii) the positive definite matrices P1 and P2 which come from the solution of
the two algebraic Riccati equations and which weigh the above-mentioned observer
and controller terms. The proposed control architecture guarantees that, the output
of the closed-loop system will asymptotically track the desired trajectory and that
H∞ performance will be achieved.
In the case of a system with non-completely measurable state vector, the structure
and the stability properties of a MIMO adaptive fuzzy controller based on differential
flatness theory have been explained in Sect. 3.4. Now, the problem of controller is
more demanding because one cannot use the complete state vector of the system in the
feedback loop. With the application of differential flatness theory, in Chap. 3 it was
shown that the 2-DOF robotic model can be decomposed into the trivial subsystems
of Eq. (3.47), that is
(r ) p
yi i = fi (x) + j=1 gij (x)uj + dj (5.78)
As noted, the control of the robotic system already depicted in Fig. 3.7 described by
Eq. (3.70) and by its equivalent canonical form of Eq. (3.49) becomes more compli-
cated when the state vector x is not directly measurable and has to be reconstructed
through a state observer. The following definitions are used
202 5 Differential Flatness Theory and Industrial Robotics
Applying Eq. (5.80) to the nonlinear system described by Eq. (5.79), results into
(r)
y(r) = f (x) + g(x)ĝ−1 (x̂)[−fˆ (x̂) + ym − K T ê + uc ] + d⇒
(r)
= f (x) + [g(x) − ĝ(x̂) + ĝ(x̂)]ĝ (x̂)[−fˆ (x̂) + ym − K T ê + uc ] + d⇒
y(r) −1
(r)
y(r) = [f (x) − fˆ (x̂)] + [g(x) − ĝ(x̂)]u + ym − K T ê + uc + d
(5.81)
(r)
It holds that e = x − xm ⇒ y(r) = e(r) + ym . Substituting y(r) in the above equation
gives
(r) (r)
e(r) + ym = ym − K T ê + uc + [f (x) − fˆ (x̂)]+
(5.82)
+[g(x) − ĝ(x̂)]u + d
and equivalently
e1 = C T e (5.84)
ê1 = C T ê (5.86)
The feedback gain matrix is denoted as K ∈ Rn×p . The observation gain matrix is
denoted as Ko ∈ Rn×p and its elements are selected so as to assure the asymptotic
elimination of the observation error.
5.3 Observer-Based Adaptive Fuzzy Control of MIMO Robots 203
The dynamic model of the 2-DOF robotic manipulator was first analyzed in Chap. 3.
The considered 2-DOF robotic model consists of two rigid links which are rotated
by DC motors, as shown in Fig. 5.6. It has been proven that, the considered robotic
system is a differentially flat one.
Next, taking into account also the effects of additive disturbances to the joints of
the robotic manipulator the dynamic model becomes
where [uc1 uc2 ]T is a robust control term that is used for the compensation of the
model’s uncertainties as well as of the external disturbances and using the feedback
gain vector KiT = [k1i , k2i , . . . , kn−1
i , kni ], the closed-loop tracking error dynamics
becomes
Fig. 5.6 Rigid-link 2-DOF robotic manipulator and joints’ actuation with the use of DC motors
204 5 Differential Flatness Theory and Industrial Robotics
f1 (x, t) − fˆ1 (x, t)
ė = (A − BK T )e + Buc +B +
f2 (x, t) − fˆ2 (x, t)
(5.90)
g1 (x, t) − ĝ1 (x, t)
+ u + d̃
g2 (x, t) − ĝ2 (x, t)
When the estimated state vector x̂ is used in the feedback control loop, equivalently
to Eq. (10.182) one has
f1 (x, t) − fˆ1 (x̂, t)
ė = Ae − BK T ê + Buc + B +
f2 (x, t) − fˆ2 (x̂, t)
(5.92)
g1 (x, t) − ĝ1 (x̂, t)
+ u + d̃
g2 (x, t) − ĝ2 (x̂, t)
The associated state observer will be described again by Eqs. (5.85) and (5.86).
For the transformed state vector of the robot, the observation error is defined as
ẽ = e − ê = x − x̂. Substructing Eq. (5.85) from Eq. (5.83) as well as Eq. (5.86) from
Eq. (5.84) one gets
5.3 Observer-Based Adaptive Fuzzy Control of MIMO Robots 205
e1 − ê1 = C T (e − ê)
or equivalently
ẽ˙ = Aẽ + Buc + B{[f (x, t) − fˆ (x̂, t)] + [g(x, t) − ĝ(x̂, t)]u + d̃} − Ko C T ẽ
ẽ1 = C T ẽ
ẽ˙ = (A − Ko C T )ẽ + Buc + B{[f (x, t) − fˆ (x̂, t)] + [g(x, t) − ĝ(x̂, t)]u + d̃}
(5.95)
ẽ1 = C T ẽ (5.96)
ẽ1 = C T ẽ (5.98)
Next, the following approximators of the system’s unknown dynamics are defined
ˆ 4×1 fˆ (x̂|θ ) ∈ R1×1
ˆf (x̂) = f1 (x̂|θf ) x̂ ∈ R 1 f
(5.99)
fˆ2 (x̂|θf ) x̂ ∈ R4×1 fˆ2 (x̂|θf ) ∈ R1×1
where l = 1, 2, x̂ is the estimate of the state vector and μAi (x̂) is the ith membership
j
function of the antecedent (IF) part of the lth fuzzy rule. Similarly, the following
approximators of the unknown system dynamics are defined with the use of the
estimated state vector x̂
Using Eq. (5.93), the value of the approximation error that corresponds to the optimal
values of the weights vectors θf∗ and θg∗ one has
w = f (x, t) − fˆ (x̂|θf∗ ) + g(x, t) − ĝ(x̂|θg∗ ) u (5.104)