HAL Id: inria-00142588
https://hal.inria.fr/inria-00142588
Submitted on 19 Apr 2007
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
A revised look at numerical dierentiation with an
application to nonlinear feedback control
Mamadou Mboup, Cédric Join, Michel Fliess
To cite this version:
Mamadou Mboup, Cédric Join, Michel Fliess. A revised look at numerical dierentiation with an ap-
plication to nonlinear feedback control. The 15th Mediterrean Conference on Control and Automation
- MED’2007, Jun 2007, Athènes, Greece. �inria-00142588�
A revised look at numerical differentiation with an application to
nonlinear feedback control
Mamadou MBOUP
∗†
, C´edric JOIN
∗‡
, Michel FLIESS
∗§
∗
Projet ALIEN, INRIA-Futurs
†
UFR Math´ematiques et Informatique, Universit´e Ren´e Descartes (Paris 5),
45 rue des Saints-P`eres, 75270 Paris cedex 06, France
mboup@math-info.univ- paris5.fr
‡
CRAN (CNRS-UMR 7039), Universit´e Henri Poincar´e (Nancy I), BP 239,
54506 Vandœuvre-l`es-Nancy, France
cedric.join@cran.uhp-nancy.fr
§
´
Equipe MAX, LIX (CNRS-UMR 7161),
´
Ecole polytechnique,
91128 Palaiseau, France
Michel.FLiess@polytechnique.edu
Abstract—We are presenting new and efficient methods for
numerical differentiation, i.e., for estimating derivatives of a
noisy time signal. They are illustrated, via convincing numerical
simulations, by the analysis of an academic signal and by the
feedback control of a nonlinear system.
I. INTRODUCTION
A. Numerical differentiation
Numerical differentiation, i.e., the derivatives estimation
of noisy time signals, is an important but difficult ill-posed
theoretical problem. It has attracted a lot of attention in many
fields of engineering and applied mathematics (see, e.g., in
the recent control literature [2], [3], [4], [14], [15], [16], [17],
[18], [19], [23], and the references therein). Our purpose here
is to improve a new approach which started in [8], [13], and
in [7], [9], [10], for solving various questions in control and
in signal and image processing. Let us briefly describe our
differentiators which are obtained via iterated time integrals
of the noisy signal.
B. Short summary of our approach
Start with the first degree polynomial time function
p
1
(t) = a
0
+ a
1
t, t ≥ 0, a
0
, a
1
∈ R. Rewrite thanks to
classic operational calculus (cf. [25]) p
1
as P
1
=
a
0
s
+
a
1
s
2
.
Multiply both sides by s
2
:
s
2
P
1
= a
0
s + a
1
(1)
Take the derivative of both sides with respect to s , which
corresponds in the time domain to the multiplication by −t
(cf. [25]):
s
2
dP
1
ds
+ 2sP
1
= a
0
(2)
The coefficients a
0
, a
1
are obtained via the triangular system
of equations (1)-(2). We get rid of the time derivatives, i.e., of
sP
1
, s
2
P
1
, and s
2
dP
1
ds
, by multiplying both sides of equations
(1)-(2) by s
−n
, n ≥ 2. The corresponding iterated time
integrals are low pass filters which attenuate the corrupting
noises, which are viewed as highly fluctuating phenomena
(see [6] for more details). A quite short time window is
sufficient for obtaining accurate values of a
0
, a
1
.
The extension to polynomial functions of higher degree is
straightforward. For derivatives estimates up to some finite
order of a given smooth function f : [0 , +∞) → R, take
a suitable truncated Taylor expansion around a given time
instant t
0
, and apply the previous computations. Utilizing
sliding time windows permit to estimate derivatives of vari-
ous orders at any sampled time instant.
C. Difficulties and improvements
The above method becomes more and more ill-conditioned
for higher order truncation of the Taylor expansion. This
is a major impediment for obtaining good estimates for
higher order derivatives in a noisy setting. Using elimination
techniques and Jacobi orthogonal polynomials (see, e.g., [1],
[24]), we propose here individual and independent deriv-
atives estimators for each given order. A judicious choice
of the point at which the derivatives are estimated in each
sliding time window permits to take advantage of the extra-
modeling capability afforded by a higher order truncation.
D. Organization of our paper
Section II discusses the mathematical foundations of our
differentiators. The first illustration in Section III applies the
above techniques for estimating the derivative of a noisy
academic signal. The second illustration in Section IV, which
is borrowed from [5], deals with the nonlinear feedback con-
trol of a DC motor joined to an inverted pendulum through
a torsional spring. We provide in both cases convincing
computer simulations
1
.
1
Interested readers may obtain the corresponding computer programs
from one of the authors (cedric.join@cran.uhp-nancy.fr). Ref-
erence [11], from which the second example is taken, contains many more
applications to various topics in nonlinear control. Most useful discussions
and comparisons may be found in [21] where an interesting concrete case-
study is analyzed.
II. DERIVATIVES ESTIMATION
Let y(t) = x(t) + n(t) be a noisy observation on a finite
time interval of a real-valued signal x(t), the derivatives of
which we want to estimate. Assume that x(t) is analytic on
this time interval. Consider without any loss of generality
the convergent Taylor expansion x(t) =
P
i>0
c
i
t
i
i!
at t = 0.
The truncated Taylor expansion x
N
(t) =
P
N
i=0
c
i
t
i
i!
satisfies
the differential equation
d
N +1
dt
N +1
x
N
(t) = 0 . It reads in the
operational domain as
s
N+1
ˆx
N
(s) = s
N
x
N
(0) + s
N−1
˙x
N
(0) . . . + x
(N)
N
(0) (3)
where ˆx
N
(s) is the operational analog of x
N
(t). To ease the
notation, we subsequently ignore the argument s.
A. Simultaneous estimation
Replace x
N
(t) by the noisy observed signal y(t). Then the
estimates ˜x
(i)
N
(0) of the derivatives at the origin x
(i)
N
(0)
△
=
d
i
dt
i
x(t)
t=0
are directly obtained from the linear triangular
system of equations (see [8], [13])
s
−ν
d
m
ds
m
n
˜x
N
(0)s
N
+ . . . + ˜x
(N−1)
N
(0)s + ˜x
(N)
N
(0)
o
=
s
−ν
d
m
ds
m
n
s
N+1
ˆy
o
(4)
m = 0, . . . , N , where ν > N + 2 ensures strict properness.
To obtain the numerical estimates, it suffices to express
(4) back in the time domain, using the classical rules of
operational calculus ([25]). Denote by T the estimation time.
We end up with the following closed form expression
P
ν
(T )
˜x
N
(0)
˜
˙x
N
(0)
.
.
.
˜x
(N)
N
(0)
=
Z
T
0
Q
ν
(τ)y(τ)dτ (5)
where the nonzero entries of the triangular matrix P
ν
(T ) are
given, for i = 0, . . . , N , j = 0, . . . , N − i, by
{P
ν
(T )}
ij
=
(N − j)!
(N − i − j)!
T
ν−N+i+j−1
(ν − N + i + j − 1)!
and
{Q
ν
(τ)}
i
=
i
X
ℓ=0
q
i,ℓ
(T − τ)
ν−N−2−ℓ
τ
i−ℓ
with
q
i,ℓ
=
i
ℓ
(N + 1)!
(N + 1 − ℓ)!
(−1)
i−ℓ
(ν − N − 2 − ℓ)!
Finally, for each estimation time interval of length T , I
T
t
+
=
[t, t + T ], we obtain the derivatives estimates at time t by
replacing y(t) in (5) by y(t+τ ). These estimates are however
not causal. To obtain causal estimates, i.e., the estimates at
time t, based on the signal observation in I
T
t
−
= [t − T, t], it
suffices to replace y(t + τ) by −y(t − τ ), τ ∈ [0, T ].
B. Individual estimation
The matrix P
ν
(T ) in (5) is in general ill-conditioned,
and yields therefore poor estimates especially in a noisy
setting. A solution to this problem is to obtain an independent
estimator for each order of derivation
2
. Reconsider (3).
Examine for 0 6 n 6 N the n
th
order derivative. Annihilate
the remaining coefficients x
(j)
N
(0), j 6= n by multiplying by
linear differential operators of the form
Π
N,n
κ
=
d
n+κ
ds
n+κ
·
1
s
·
d
N−n
ds
N−n
κ > 0
It yields the following estimator for x
(n)
(0)
˜x
(n)
N
(0)
s
ν+n+κ+1
=
(−1)
n+κ
(n + κ)!(N − n)!
1
s
ν
Π
N,n
κ
s
N+1
ˆx
(6)
which is strictly proper whenever ν is of the form ν =
N + 1 + µ, µ > 0. We obtain a family of strictly proper
estimators which is parametrized by κ, µ and N . Write
˜x
(n)
0
(κ, µ; N ) the corresponding estimator. If N = n the
differential operator Π
N,n
κ
reduces to Π
n
κ
=
d
n+κ
ds
n+κ
·
1
s
: we will
then use the simplified notation ˜x
(n)
0
(κ, µ). The following
result is straightforward:
Lemma 1: For any N > n and µ > 0, ˜x
(n)
0
(κ, µ; N ) in
(6) belongs to the set
F = span
Q
{˜x
(n)
0
(κ
ℓ
, µ
ℓ
), ℓ = 0, . . . , min (n + κ, N − n)} (7)
where κ
ℓ
= κ + N − n − ℓ and µ
ℓ
= µ + ℓ.
Proof: Set q = N − n and p = n+κ. The proof follows
by direct inspection, upon writing Π
N,n
κ
(s
N+1
ˆx) in the form
Π
N,n
κ
(s
N+1
ˆx) =
q
X
i=0
q
i
!
(q + 1)!
(q + 1 − i)!
d
p
ds
p
n
s
q−i
(s
n
ˆx)
(q−i)
o
=
q
X
i=0
min (p,q−i)
X
j=i
a
i,j
s
q−j
Π
n
κ+q−j
{s
n+1
ˆx}
(8)
where
a
i,j
=
q
i
p
j − i
(q + 1)!
(q + 1 − i)(q − j)!
This lemma shows that an n
th
-order truncated Taylor expan-
sion is appropriate for estimating the n
th
-order derivative.
C. Least squares interpretation
3
A common way for estimating the derivatives of a signal
is to resort to a least squares polynomial fitting on an interval
and then take the derivatives of the resulting polynomial
function. The estimators derived here rely however on a
different approach: the derivatives are estimated pointwise.
This depature is furthermore apparent with the developments
of the preceding subsection. Nonetheless, a least squares
interpretation may be attached to our approach, as shown
below.
Start with the estimation of the first order derivative.
2
The system (5) being triangular, a closed-form expression for the
estimator of x
(i)
(t) may be derived from it. The corresponding solutions
would however exhibit the same sensitivity to noise perturbations.
3
The authors would like to thank A. Sedoglavic for bringing this question
to their attention.
a) N = 1: With ν = µ + 2, equation (6) becomes
1
s
µ+κ+4
˜
˙x
0
(κ, µ) =
(−1)
κ+1
(κ + 1)!
ˆx
(κ+1)
s
µ+1
+ (κ + 1)
ˆx
(κ)
s
µ+2
(9)
It reads in the time domain:
˜
˙x
0
(κ, µ) =
µ + 2
T
µ + κ + 3
κ + 1
Z
1
0
p(τ)τ
κ
(1 − τ)
µ
y(T τ )dτ
(10)
where p(τ) = (µ+κ+2)τ −(κ+1), and T is the estimation
time (the estimation interval is I
T
0
+
). We replace of course
the signal x by its noisy observation y.
Consider now the Jacobi orthogonal polynomials (cf.
[24], [1]) {P
κ,µ
i
(t)}
i>0
, associated to the weight function
w
κ,µ
(t) = t
κ+1
(1 − t)
µ+1
on the interval [0, 1].
Lemma 2: The first order derivative estimate, given in
equation (10), reads as
˜
˙x
0
(κ, µ) =
1
kP
κ,µ
0
(t)k
2
Z
1
0
P
κ,µ
0
(τ)w
κ,µ
(τ) ˙y(T τ )dτ
△
=
hP
κ,µ
0
(τ), ˙y(T τ)i
κ,µ
kP
κ,µ
0
(τ)k
2
(11)
Proof: Observe that P
κ,µ
0
(t) = 1 and p(τ)τ
κ
(1−τ )
µ
=
−
d
dτ
{τ
κ+1
(1 − τ )
µ+1
}. Integration by parts shows that
the integral in (10) reduces to T hP
κ,µ
0
(τ), ˙y(T τ)i
κ,µ
. The
equality k P
κ,µ
0
k
2
=
(µ+1)!(κ+1)!
(κ+µ+3)!
completes the proof.
This estimate of the first order derivative appears as the
orthogonal projection of the unobserved signal derivatives,
on P
κ,µ
0
(t). Expanding ˙x(T τ ), τ ∈ [0, 1] in the basis of the
Jacobi polynomials
˙x(T τ) = a
0
P
κ,µ
0
(τ) + a
1
P
κ,µ
1
(τ) + a
2
P
κ,µ
2
(τ) + · · · (12)
shows that
˜
˙x
0
(κ, µ), which coincides with a
0
, actually turns
out to correspond to an estimate of ˙x(T τ
0
) for some τ
0
> 0.
Using a first order approximation of ˙x(T τ)
˙x(T τ
0
) ≈ a
0
P
κ,µ
0
(τ
0
) + a
1
P
κ,µ
1
(τ
0
) = a
0
allows one to identify τ
0
as the solution of P
κ,µ
1
(τ
0
) = 0.
This value of τ
0
, given by
τ
0
=
κ + 2
µ + κ + 4
(13)
is “experimentally” confirmed by the numerical simulations
below. The resulting derivative estimation is thus subject to
a time delay.
b) N = 2: It allows to avoid such a delay, which may
not be tolerable in real-time processing. Equation (8) yields
˜
˙x
0
(κ, µ; 2)
s
µ+κ+5
=
(−1)
κ+1
(κ + 1)!
ˆx
(κ+2)
s
µ+1
+ (κ + 2)
ˆx
(κ+1)
s
µ+2
+ (κ + 3)
ˆx
(κ+1)
s
µ+2
+ (κ + 1)
ˆx
(κ)
s
µ+3
(14)
It reads in the time domain after some algebraic manipula-
tions:
˜
˙x
0
(κ, µ; 2) =
hP
κ,µ
0
(τ), ˙y(T τ)i
κ,µ
kP
κ,µ
0
k
2
+ P
κ,µ
1
(0)
hP
κ,µ
1
(τ), ˙y(T τ)i
κ,µ
kP
κ,µ
1
k
2
(15)
where P
κ,µ
1
(0) = −τ
0
= −
κ+2
µ+κ+4
. We therefore deduce
that
˜
˙x
0
(κ, µ; 2) corresponds to an estimate of the first order
derivative ˙x (t) at t = 0. The estimation of ˙x(t) from a second
order truncation of the Taylor expansion is therefore delay
free.
We now show how this interpretation can be exploited to
obtain better estimates. According to Lemma 1 and using
(14), it is easy to verify that the relation
˜
˙x
0
(κ, µ; 2) = λ
0
˜
˙x
0
(κ, µ + 1) + λ
1
˜
˙x
0
(κ + 1, µ) (16)
holds for λ
0
= κ + 3 and λ
1
= −(κ + 2). Let us now
extend the set F in (7), by allowing the coefficients of the
linear combinations therein to be real, rather than rational. In
particular, given a point τ
1
∈ [0, 1], one may always choose
the λ
i
’s such that (16) becomes
˜
˙x
τ
1
(κ, µ; 2) =
hP
κ,µ
0
(τ), ˙y(T τ)i
κ,µ
kP
κ,µ
0
k
2
+ P
κ,µ
1
(τ
1
)
hP
κ,µ
1
(τ), ˙y(T τ)i
κ,µ
kP
κ,µ
1
k
2
(17)
In this case,
˜
˙x
τ
1
(κ, µ; 2) will represent the estimate of
˙x(T τ
1
), obtained from the truncated Taylor expansion
x
2
(t) = x(τ)+ ˙x(τ)(t−τ)+
¨x(τ )
2
(t− τ )
2
= x
2
(0)+ [ ˙x(τ)−
τ ¨x(τ)]t +
¨x(τ )
2
t
2
at τ = T τ
1
. A direct verification will show
that the values for λ
0
and λ
1
, associated to τ
1
are given by
λ
0
= (κ + 3) − (µ + κ + 5)τ
1
and λ
1
= 1 − λ
0
.
Let us now turn to the question pertaining to the selection
of a value for τ
1
. Expansion (12) shows that a good choice
for τ
1
is given by the smallest (resp. largest) root of the
polynomial P
κ,µ
2
when the estimation interval is I
T
t
+
(resp.
I
T
t
−
). Indeed, choosing τ
1
as a zero of P
κ,µ
2
annihilates the
contribution of the orthogonal projection of the signal on
P
κ,µ
2
in the estimation error. On the other hand, recall that
τ
1
(resp. 1 − τ
1
for I
T
t
−
) represents a delay in the estimation.
Choosing the smallest (resp. largest) root thus translates to
a lower delay. Note also that the delay τ
1
is smaller than τ
0
,
since the zeros of P
κ,µ
2
and P
κ,µ
1
interlace.
The above analysis may be easily generalized to the n
th
order derivative estimation, n ≥ 2.
III. FIRST ILLUSTRATION: DERIVATIVE OF A NOISY
SIGNAL
Let y(t) = x(t)+n(t), 0 ≤ t ≤ 5, be a noisy measurement
of the signal
x(t) = tanh(t − 1) + e
−t/1.2
sin(6t + π)
The noise level, measured by the signal to noise ratio
in dB, i.e., SNR = 10 log
10
P
|y(t
i
)|
2
P
|n(t
i
)|
2
, corresponds to
SN R = 25 dB (see Figure 1). In all the subsequent
numerical simulations, the integrals are computed via the
classical trapezoidal rule.
Begin with the first order derivative and N = 1. The
estimates
˜
˙x
0
(κ, µ), obtained from (10) with κ = µ = 2,
are displayed in Figure 2 (solid line). It corresponds to
the estimation results in the successive intervals I
T
t+
, with
T = 60T
s
and for t = iT
s
, i = 0, . . . , ⌊
5−T
T
s
⌋. The exact
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
0.0
0.5
1.0
1.5
−2.0
−1.5
−1.0
−0.5
Fig. 1. Noisy observation signal, SNR = 25dB.
τ
0
T
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
−6
−4
−2
0
2
4
6
Fig. 2. Estimation of the signal derivative: N = 1.
derivative of the noise-free signal is also displayed (dashed
line) in order to gauge the estimation accuracy. For this,
we have introduced a shift corresponding to the delay τ
0
in (13). Observe how this predicted value of τ
0
fits with
the experiment. Of course with N = 1, the truncated
Taylor series model is linear, resulting in poor estimates
on the intervals where the signal’s dynamic is strong. For
high signal-to-noise ratio, the estimates may be improved
by reducing the estimation time T . Alternatively, one may
consider a richer signal model, e.g., with N = 2. This is the
case of the next experiment. We now consider the sliding
windows I
T
t−
, with T = 110T
s
. The estimates
˜
˙x
0
(κ, µ; 2),
based on (14), (15), are plotted (solid line) in Figure 3 below,
for κ = µ = 0.
There is no estimation delay, as expected. However, the
performance significantly degrades as compared to the pre-
ceding results although the signal model is more precise. If
we now relax this delay-free constraint, it becomes possible
to take advantage of the more flexible second order model
for the signal.
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
−6
−4
−2
0
2
4
6
Fig. 3. Estimation of the signal derivative: N = 2, no delay
This is illustrated in the following simulation (see Figure
4), where we keep the same settings for T , κ and µ. The
solid line curve in Figure 4 represents the estimates obtained
τ T
2
τ
1
T
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
−6
−4
−2
0
2
4
6
Fig. 4. Estimation of the signal derivative: N = 2, with delay τ
1
.
from (16)-(17) where τ
1
is the largest root of P
κ,µ
2
. Using the
same idea with N = 3, we obtain the estimates
˜
˙x
0
(κ, µ; 3),
plotted with dashed line. The associated estimator reads as
˜
˙x
0
(κ, µ; 3) =
2
X
i=0
λ
i
˜
˙x
0
(κ + 2 − i, µ + i)
where the real coefficients λ
i
, i = 0, 1, 2 are chosen so that
˜
˙x
0
(κ, µ; 3) corresponds to ˙x(τ
2
T ), with τ
2
being a root of
P
κ,µ
3
. If κ = µ, which is the case here, then τ
2
= 0.5 is a
common root of all Jacobi polynomials of odd degree. The
dashed line curve was obtained with this choice of τ
2
.
IV. SECOND ILLUSTRATION: NONLINEAR FEEDBACK
CONTROL
A. System description
Consider with [5] the mechanical system, which consists
of a DC-motor joined to an inverted pendulum through a
torsional spring:
J
m
¨
θ
m
= κ
θ
l
− θ
m
− B
˙
θ
m
+ K
τ
u
J
l
¨
θ
l
= −κ
θ
l
− θ
m
− mgh sin(θ
l
)
y = θ
l
