Gevorderde Systeemidentificatie

1/243
System Identification
Rik Pintelon
Vrije Universiteit Brussels, dept. ELEC
Version 25 November 2011
2/243
Table of Contents
1. Preliminary Example 3
3. Nonparametric models of LTI systems 29
4. Parametric models of LTI systems 71
5. Improved nonparametric models for LTI systems 113
6. Estimation parametric plant model with known noise model 146
7. Estimation parametric plant model with estimated nonparametric noise model 174
8. Estimation parametric plant and/or noise models 200
9. Guidelines 232
10. Books on system identification 242
3/243
1. Preliminary Example
2.a Measurement of the value of a resistor 4
2.b Least squares estimator of the output error problem 6
2.c Maximum likelihood estimator of the errors-in-variables problem 7
2.d Properties estimators 10
2.e Stochastic convergence estimators 11
2.f Asymptotic variance estimators 12
2.g Asymptotic normality estimators 15
2.h Asymptotic efficiency estimators 17
2.i Robustness estimators 19
2.j Summary properties 21
2.k Extension to arbitrary currents 23
2.l References 28
4/243
2.a Measurement of the value of a resistor
0
5
0 50 100
M
e
a
s
u
r
e
d

r
e
s
i
s
t
a
n
c
e

(
)
0
2
0 50 100
V
o
l
t
a
g
e

(
V
)
0
2
0 50 100
C
u
r
r
e
n
t

(
A
)
k
k k
i k ( ) i
0
n
i
k ( ) + =
u k ( ) u
0
n
u
k ( ) + =
n
i
k ( ) N 0
i
2
, ( )
n
u
k ( ) N 0
u
2
, ( )
u k ( )
i k ( )
R k ( ) u k ( ) i k ( ) =
u k ( )
i k ( )
R
0
1 =
V
A
5/243
2.a Measurement of the value of a resistor (contd)
Observation: the estimate seems not to converge as
Question: what is going wrong?
0
1
2
1 10 100 1000 10000
E
s
t
i
m
a
t
e
d

v
a
l
u
e

R
Number of measurements
R
SA
N ( )
1
N
--- - R k ( )
k 1 =
N
1
N
--- -
u k ( )
i k ( )
---------
k 1 =
N
= =
N
6/243
2.b Least squares estimator of the output error problem
Observation: the estimate converges to a too small value as
Question: what is going wrong?
0
1
2
1 10 100 1000 10000
E
s
t
i
m
a
t
e
d

v
a
l
u
e

R
V
LS
R z , ( ) u k ( ) Ri k ( ) ( )
2
k 1 =
N
=
R
LS
N ( )
arg min
R
V
LS
R z , ( )
u k ( )i k ( )
k 1 =
N
i
2
k ( )
k 1 =
N
---------------------------------
= =
N
7/243
2.c Maximum likelihood estimator of the errors-in-variables problem
Call the vector of the measurements
Gaussian iid random variables probability density function of the measurements
Constrained negative log-likelihood function
s.t.
Maximum likelihood cost function
Eliminate , and
z
z u 1 ( ) u 2 ( ) u N ( ) i 1 ( ) i 2 ( ) i N ( ) , , , , , , , [ ]
T
=
z
f
z
z R u
p
i
p
, , ( ) f
u k ( )
u k ( ) R u
p
i
p
, , ( )f
i k ( )
i k ( ) R u
p
i
p
, , ( )
k 1 =
N
=
1
2
u
i
( )
N
------------------------- -exp
1
2
---
u k ( ) u
p
( )
2
u
2
----------------------------
i k ( ) i
p
( )
2
i
2
-------------------------
+
k 1 =
N
( ) =
f
z
z R u
p
i
p
, , ( ) log constant
1
2
---
u k ( ) u
p
( )
2
u
2
----------------------------
i k ( ) i
p
( )
2
i
2
-------------------------
+
k 1 =
N
+ = u
p
Ri
p
=
V
ML
R u
p
i
p
z , , , ( )
1
2
---
u k ( ) u
p
( )
2
u
2
----------------------------
i k ( ) i
p
( )
2
i
2
-------------------------
+
k 1 =
N
u
p
Ri
p
( ) + =
u
p
i
p
R
ML
N ( )
8/243
2.c Maximum likelihood estimator of the errors-in-variables problem (contd)
Observation: estimate converges to the true value as
R
ML
N ( )
1
N
---- u k ( )
k 1 =
N
1
N
--- - i k ( )
k 1 =
N
-----------------------------
=
0
1
2
1 10 100 1000 10000
E
s
t
i
m
a
t
e
d

v
a
l
u
e

R
N
9/243
2.c Maximum likelihood estimator of the errors-in-variables problem (contd)
Observation: estimate is different for each realisation stochastic process
Questions: - point wise convergence ?
- for almost all realisations?
- ?
- minimal uncertainty?
R
ML
r [ ]
N ( )
1
N
---- u
r [ ]
k ( )
k 1 =
N
1
N
--- - i
r [ ]
k ( )
k 1 =
N
----------------------------------
=
r 1 2 , , =
0.6
0.8
1
1.2
1.4
0 20 40 60 80 100
R
e
s
i
s
t
a
n
c
e

(
)
N
R
ML
r [ ]
N ( )
N
lim R
ML
r [ ]
=
R
ML
r [ ]
R
ML
=
R
ML
R
0
=
10/243
2.d Properties estimators
stochastic convergence
asymptotic variance
asymptotic distribution
asymptotic efficiency
robustness
Question: how to predict these properties?
R
ML
N ( )
1
N
---- u k ( )
k 1 =
N
1
N
--- - i k ( )
k 1 =
N
-----------------------------
=
0
1
2
1 10 100 1000 10000
E
s
t
i
m
a
t
e
d

v
a
l
u
e

R
R
LS
N ( )
u k ( )i k ( )
k 1 =
N
i
2
k ( )
k 1 =
N
---------------------------------
=
R
SA
N ( )
1
N
--- -
u k ( )
i k ( )
---------
k 1 =
N
=
R
SA
N ( )
R
LS
N ( )
R
ML
N ( )
11/243
2.e Stochastic convergence estimators
Basic tool: the law of large numbers
then as
Simple approach
for
no convergence
Least squares
convergence to the wrong value = inconsistent
Maximum likelihood
convergence to the true value = consistent
S N ( )
1
N
--- - x k ( )
k 1 =
N
= S N ( ) S N ( ) { }
1
N
--- - x k ( ) { }
k 1 =
N
= N
R
SA
N ( )
1
N
--- -
u k ( )
i k ( )
---------
k 1 =
N
=
1
N
--- - u k ( ) { }
1
i k ( )
--------

k 1 =
N
= i k ( ) N i
0

i
2
, ( )
R
LS
N ( )
1
N
--- - u k ( )i k ( )
k 1 =
N
1
N
---- i
2
k ( )
k 1 =
N
--------------------------------------
=
1
N
---- u k ( ) { }E i k ( ) { }
k 1 =
N
1
N
--- - i
2
k ( ) { }
k 1 =
N
-----------------------------------------------------------
u
0
i
0
i
0
2

i
2
+
----------------
R
0
1
i
2
i
0
2
+
-----------------------
= =
R
ML
N ( )
1
N
---- u k ( )
k 1 =
N
1
N
--- - i k ( )
k 1 =
N
-----------------------------
=
1
N
--- - u k ( ) { }
k 1 =
N
1
N
---- i k ( ) { }
k 1 =
N
---------------------------------------
u
0
i
0
-----
R
0
= =
12/243
2.f Asymptotic variance estimators
Basic tool: linearisation estimate around the limit value
Least squares
the asymptotic ( ) variance
Notes:
- valid for any voltage and current signal-to-noise ratio
- asymptotic variance is an
R
LS
N ( )
1
N
--- - u k ( )i k ( )
k 1 =
N
1
N
--- - i
2
k ( )
k 1 =
N
--------------------------------------
u
0
i
0
w
1 [ ]
+
i
0
2
w
2 [ ]
+
-------------------------
R
0
1
i
2
i
0
2
+
----------------------- 1
w
1 [ ]
u
0
i
0
---------
w
2 [ ]

i
2
i
0
2

i
2
+
---------------------
+

= =
w
1 [ ]
1
N
--- - u
0
n
i
k ( ) i
0
n
u
k ( ) n
i
k ( )n
u
k ( ) + + ( )
k 1 =
N
= with w
1 [ ]
{ } 0 =
w
2 [ ]
1
N
--- - 2i
0
n
i
k ( ) n
i
2
k ( ) + ( )
k 1 =
N
= with w
2 [ ]
{ }
i
2
=
N
var NR
LS
N ( ) ( )
R
0
2
1
i
2
i
0
2
+ ( )
2
------------------------------var N
w
1 [ ]
u
0
i
0
---------
w
2 [ ]

i
2
i
0
2

i
2
+
---------------------

( )
O N
1
( )
13/243
2.f Asymptotic variance estimators (contd)
Maximum likelihood
the asymptotic ( ) variance
Notes:
-valid for any voltage and current signal-to-noise ratio
- asymptotic variance is an
- expected value does not exist
- variance does not exist
R
ML
N ( )
1
N
---- u k ( )
k 1 =
N
1
N
--- - i k ( )
k 1 =
N
-----------------------------
u
0
w
1 [ ]
+
i
0
w
2 [ ]
+
--------------------- R
0
1
w
1 [ ]
u
0
---------
w
2 [ ]
i
0
---------
+

= =
w
1 [ ]
1
N
--- - n
u
k ( )
k 1 =
N
= with w
1 [ ]
{ } 0 =
w
2 [ ]
1
N
--- - n
i
k ( )
k 1 =
N
= with w
2 [ ]
{ } 0 =
N
var NR
ML
N ( ) ( ) R
0
2
var N
w
1 [ ]
u
0
---------
w
2 [ ]
i
0
---------

( ) R
0
2

i
2
i
0
2
------

u
2
u
0
2
------
+

=
O N
1
( )
R
ML
N ( )
R
ML
N ( )
14/243
2.f Asymptotic variance estimators (contd)
0.001
0.01
0.1
1
0.01 0.1 1
Standard deviation noise
S
t
a
n
d
a
r
d

d
e
v
i
a
t
i
o
n

o
n
R
std R
LS
( )
std R
ML
( )
0.001
0.01
0.1
1
0.01 0.1 1
R
M
S

e
r
r
o
r
rms R
ML
( )
rms R
LS
( )
MSE R
( ) E R
R
0
( )
2
{ } var R
( ) R
{ } R
0
( )
2
+ = =
u
u
0
-----

i
i
0
----
=

u
u
0
-----

i
i
0
----
=
15/243
2.g Asymptotic normality estimators
Basic tool: linearisation estimate around the limit value + the central limit theorem
then as
Least squares and maximum likelihood estimators
Application CLT to , , are asymptotically normally distributed
, are asymptotically normally distributed
S N ( )
1
N
--- - x k ( )
k 1 =
N
=
S N ( ) S N ( ) { }
var S N ( ) ( )
-------------------------------------- AsN 0 1 , ( ) N
R
LS
N ( )
R
0
1
i
2
i
0
2
+
----------------------- 1
w
1 [ ]
u
0
i
0
---------
w
2 [ ]

i
2
i
0
2

i
2
+
---------------------
+

ML
N ( ) R
0
1
w
1 [ ]
u
0
---------
w
2 [ ]
i
0
---------
+

w
1 [ ]
w
2 [ ]
w
1 [ ]
w
2 [ ]
R
LS
N ( ) R
ML
N ( )
16/243
2.g Asymptotic normality estimators (contd)
0
25
0.5 1 1.5 2
0
25
0.5 1 1.5 2
0
25
0.5 1 1.5 2
p
d
f
R
R R R
p
d
f
R
p
d
f
R
N 10 = N 100 =
N 1000 =
R
LS
R
ML
17/243
2.h Asymptotic efficiency estimators
Basic tool: comparison asymptotic uncertainty with the Cramr-Rao lower bound
with
he Cramr-Rao lower bound
with
the corresponding Fisher information matrix

z
f
z
z ( ) f
u k ( )
u k ( ) ( )f
i k ( )
i k ( ) ( )
k 1 =
N
=
1
2
u
i
( )
N
------------------------- -exp
1
2
-- -
u k ( ) Ri
p
( )
2
u
2
------------------------------
i k ( ) i
p
( )
2
i
2
-------------------------
+
k 1 =
N
( ) =
i
p
R , [ ]
T
=
Cov
( ) Fi
1

0
( ) Fi
0
( )

2
f
z
z ( ) log
0
2
-------------------------------

=
Fi
1

0
( )
1
N
----

i
2
u
0
i
2
i
0
2
------
u
0
i
2
i
0
2
------ R
0
2

i
2
i
0
2
------

u
2
u
0
2
------
+

= Fi
1
R
0
( )
R
0
2
N
------

i
2
i
0
2
------

u
2
u
0
2
------
+

=
18/243
2.h Asymptotic efficiency estimators (contd)
Conclusion: ML-estimate is asymptotically efficient
std R
ML
( )
+
CR R
0
( ) Fi
1 2 /
R
0
( ) =
0.001
0.01
0.1
1
0.01 0.1 1
u
u
0
-----

i
i
0
----
=
19/243
2.i Robustness estimators
Basic question: are the (asymptotic) properties of the estimator sensitive to the (noise)
assumptions made to construct the estimator?
Maximum likelihood estimator
- basic assumption: zero mean Gaussian iid noise on current and voltage
measurements
- relax assumption: zero mean non-Gaussian mixing noise on current and
voltage measurements
consequences on the properties:
- consistency, asymptotic normality, convergence rate still valid
- asymptotic efficiency is lost
- expression asymptotic variance still valid for independently distributed noise
Least squares estimator
- basic assumption: zero mean iid noise on voltage measurements, and
current known exactly
- relax assumption:
(i) zero mean non-Gaussian mixing noise on voltage measurements, and
current known exactly still consistent
(ii) also noise on input inconsistent
20/243
2.i Robustness estimators (contd)
Note:
0.6
0.8
1
1.2
1.4
0 100 200 300 400 500
R
e
s
i
s
t
a
n
c
e

(
)
N
R
ML
N ( )
R
LS
N ( )
R
SA
N ( )
u
0
1 V =
i
0
1 A =
n
u
k ( ) 0.5 V 0.5 V , [ ] uniformly
n
i
k ( ) 0.5 A 0.5 A , [ ] uniformly
u k ( ) u
0
n
u
k ( ) + =
i k ( ) i
0
n
i
k ( ) + =
R
SA
N ( ) { } R
0
i
0
2 3
i
---------------
1 3
i
i
0
+
1 3
i
i
0

-----------------------------

log 1.1 < = =
R
LS
N ( )
R
0
1
i
2
i
0
2
+
----------------------- 0.92 =
21/243
2.j Summary properties
Goal asymptotic analysis
- what happens with the estimate if more data are gathered?
- hypothesis: the asymptotic behaviour reflects the finite sample behaviour
- true finite sample behaviour is in general very difficult to establish
(exception: linear least squares)
Consistency
- convergence in stochastic sense to the true value
- does not exclude divergence for some realisations
- guarantees with high probability that the estimate is close to the true value
- consistency does not imply asymptotic unbiasedness
- proof: law of large numbers
Asymptotic unbiasedness
- asymptotically the expected value equals the true value
- does not guarantee that the estimate is close to the true value
- in general very difficult to verify (exception: linear least squares)
Asymptotic normality
- allows to construct uncertainty bounds with a given confidence level
- proof: linearisation around limit value + central limit theorem
22/243
2.j Summary properties (contd)
Asymptotic variance
- measure of the convergence rate
- construction of uncertainty bounds
- proof: linearisation around limit value
Asymptotic efficiency
- minimal uncertainty within the class of asymptotically unbiased estimators
- in practice applied to the class of consistent estimators
- Cramr-Rao lower bound does not always exist (e.g. uniformly distributed
noise), or may be too conservative (minimal uncertainty may be higher than
the CR-bound)
- proof: comparison asymptotic variance with Cramr-Rao lower bound
Robustness
- sensitivity (asymptotic) properties to the noise assumptions made to
construct the estimator
- proof: validity law of large numbers and central limit theorem under relaxed
noise assumptions
23/243
2.k Extension to arbitrary currents
Zero mean Gaussian iid random variables maximum likelihood cost function
Elimination , , , gives
which is a weighted nonlinear least squares problem
- starting values via linear least squares
- iterative Newton-Gauss minimization procedure
u k ( ) u
0
k ( ) n
u
k ( ) + =
i k ( ) i
0
k ( ) n
i
k ( ) + =
n
i
k ( ) N 0
i
2
k ( ) , ( )
n
u
k ( ) N 0
u
2
k ( ) , ( )
u k ( )
i k ( )
R
0
1 =
V
ML
R u
p
i
p
z , , , ( )
1
2
---
u k ( ) u
p
k ( ) ( )
2
u
2
k ( )
----------------------------------
i k ( ) i
p
k ( ) ( )
2
i
2
k ( )
--------------------------------
+
k 1 =
N

k
u
p
k ( ) Ri
p
k ( ) ( )
k 1 =
N
+ =
u
p
k ( ) i
p
k ( )
k
k 1 2 N , , , =
V
ML
R z , ( )
1
2
---
u k ( ) Ri k ( ) ( )
2
u
2
k ( ) R
2
i
2
k ( ) +
--------------------------------------
k 1 =
N
=
24/243
2.k Extension to arbitrary currents (contd)
Consistency
Difficulty: no explicit expression for the minimizer of exists
Strong law of large numbers applied to
Expected value cost function is minimal in consistent
R
V
ML
R z , ( )
V
N
R z , ( ) V
ML
R z , ( ) N =
V
N
R z , ( ) V
N
R ( ) V
N
R z , ( ) { }
1
2N
-------
u
0
k ( ) Ri
0
k ( ) ( )
2
{ }
u
2
k ( ) R
2
i
2
k ( ) +
-------------------------------------------------
k 1 =
N
1
2
---
+ = =
V
N
R ( ) R R
0
=
25/243
Asymptotic variance and asymptotic normality
Mean value theorem on derivative cost function
with ,
Use
(by definition of )
as (consistency)
as (law of large numbers)
then
with
Asymptotic variance
Asymptotic normality
(central limit theorem on )
V
N
' R
z , ( ) V
N
' R
0
z , ( ) V
N
'' R
1
z , ( ) R
R
0
( ) + = R
1
tR
1 t ( )R
0
+ = t 0 1 , [ ]
V
N
' R
z , ( ) 0 = R
R
1
R
0
N
V
N
'' R
1
z , ( ) V
N
'' R
0
( ) N
R
R
0

R
V
N
''
1
R
0
( )V
N
' R
0
z , ( ) =
R
{ } 0 =
var R
( ) var
R
( ) V
N
''
2
R
0
( ) V
N
' R
0
z , ( ) ( )
2
{ } =
N R
R
0
( ) AsN 0 var N
R
( ) , ( ) V
N
' R
0
z , ( )
26/243
Asymptotic efficiency
with
he Fisher information matrix equals

Conclusion
- the ML estimate is in general inefficient
- asymptotically efficient if current is known exactly
z
f
z
z ( ) f
u k ( )
u k ( ) ( ) f
i k ( )
i k ( ) ( )
k 1 =
N
=
1
2 ( )
N
u
k ( )
i
k ( )
k 1 =
N
-------------------------------------------------------exp
1
2
---
u k ( ) Ri
p
k ( ) ( )
2
u
2
k ( )
-------------------------------------
i k ( ) i
p
k ( ) ( )
2
i
2
k ( )
--------------------------------
+
k 1 =
N
( ) =
i
p
1 ( ) i
p
2 ( ) i
p
N ( ) R , , , , [ ]
T
=
Fi
0
( )

2
f
z
z ( ) log
0
2
-----------------------------

= Fi
1
R
0
( )
1
N
----V
N
''
1
R
0
( ) =
V
N
' R
0
z , ( ) ( )
2
{ } V
N
'' R
0
( ) N = var
R
( ) Fi
1
R
0
( ) =
27/243
Notes
- inefficiency ML estimate is not in contradiction with general ML properties because an
infinite number of parameters are estimated
- uncertainty on does not decrease as
- in practice a good approximation of the variance of is given by
where
with
- the variance of the estimate follows as a by product of the Newton-Gauss
minimization procedure
i
p
1 ( ) i
p
2 ( ) i
p
N ( ) , , , N
R
1
N
----V
N
''
1
R
0
( )
R z
0
, ( )
R
0
---------------------

T
R z
0
, ( )
R
0
---------------------

1
R z , ( )
R
------------------

T
R z , ( )
R
------------------

1
=
R z , ( )
1 [ ]
R z , ( )
2 [ ]
R z , ( )
N [ ]
R z , ( )
=
k [ ]
R z , ( )
u k ( ) Ri k ( )
u
2
k ( ) R
2
i
2
k ( ) +
------------------------------------------
=
28/243
2.l References
Billingsley, P. (1995). Probability and Measure. Wiley, New York.
Brillinger, D. R. (1981). Time Series: Data Analysis and Theory. McGraw-Hill, New York.
Caines, P. E. (1988). Linear Stochastic Systems. Wiley, New York.
Chow, Y. S. and H. Teicher (1988). Probability Theory: Independence, Interchangeability,
Martingales (2nd edition). Springer-Verlag, New York.
Fuller, W. A. (1987). Measurement Error Models. Wiley, New York.
Lukacs, E. (1975). Stochastic Convergence. Academic Press, New York.
Pintelon, R. and J. Schoukens (2012). System Identification: a Frequency Domain
Approach. IEEE Press, Piscataway, New Jersey (second edition).
29/243
3. Nonparametric models of LTI systems
3.a Problem statement 30
3.b Periodic excitations - errors-in-variables problem 32
3.c Measurement example: octave bandpass filter 37
3.d Arbitrary excitations - output error problem 38
3.e Simulation example: FRF and variance estimate using arbitrary input 41
3.f Nonlinear systems - introduction 43
3.g Nonlinear systems - noisy output measurements 48
3.h Simulation example: FRF estimate in the presence of nonlinear distortions 52
3.i Nonlinear systems - noisy input/output measurements 54
3.j Measurement example: the open loop gain of an operational amplifier 56
3.k Multivariable systems - periodic excitations 57
3.l MIMO experiment: identification d-axis synchronous machine 61
3.m MIMO experiment: vibrating mechanical structure 65
3.n Multivariable systems - random excitations 67
3.o MIMO simulation example: discrete-time system 68
3.p References 70
30/243
3.a Problem statement
Noisy observations
- periodic signals
and
- arbitrary signals
and
N
g
k ( )
U
g
k ( )
M
U
k ( )
U k ( ) Y k ( )
M
Y
k ( )
N
p
k ( )
G
0
j ( )
N
g
k ( ) generator noise
N
p
k ( ) process noise
M
U
k ( ) input measurement noise
M
Y
k ( )
output measurement noise
Y k ( ) Y
0
k ( ) N
Y
k ( ) + =
U k ( ) U
0
k ( ) N
U
k ( ) + =
Y
0
k ( ) G
0
j
k
( )U
0
k ( ) =
U
0
k ( ) U
g
k ( ) =
N
Y
k ( ) G
0
j
k
( )N
g
k ( ) N
p
k ( ) M
Y
k ( ) + + =
N
U
k ( ) N
g
k ( ) M
U
k ( ) + =
Y
0
k ( ) G
0
j
k
( )U
0
k ( ) =
U
0
k ( ) U
g
k ( ) N
g
k ( ) + =
N
Y
k ( ) N
p
k ( ) M
Y
k ( ) + =
N
U
k ( ) M
U
k ( ) =
31/243
3.a Problem statement (contd)
Nonparametric representation
Frequency response function
,
(Co-)variances circular complex distributed noise
,
and
G
0
j
k
( ) k 1 2 , , =
Y
2
k ( ) var N
Y
k ( ) ( ) N
Y
k ( )
2
{ } = =
U
2
k ( ) var N
U
k ( ) ( ) N
U
k ( )
2
{ } = =
YU
2
k ( ) covar N
Y
k ( ) N
U
k ( ) , ( ) N
Y
k ( )N
U
k ( ) { } = =
k 1 2 , , =
N
Y
2
k ( ) { } 0 =
N
U
2
k ( ) { } 0 =
N
Y
k ( )N
U
k ( ) { } 0 =
32/243
3.b Periodic excitations - errors-in-variables problem
Noisy input-output measurements
Observe consecutive periods of the steady state response
Sample means
,
Sample (co-)variances
, ,
M
m 1 = 2

M
N
T
M N =
data points
N
data points
Y
k ( )
1
M
---- - Y
m [ ]
k ( )
m 1 =
M
= U
k ( )
1
M
----- U
m [ ]
k ( )
m 1 =
M
2
k ( )
1
M M 1 ( )
------------------------ Y
m [ ]
k ( ) Y
k ( )
2
m 1 =
M
2
k ( )
1
M M 1 ( )
------------------------ U
m [ ]
k ( ) U
k ( )
2
m 1 =
M
2
k ( )
1
M M 1 ( )
------------------------ Y
m [ ]
k ( ) Y
k ( ) ( ) U
m [ ]
k ( ) U
k ( ) ( )
m 1 =
M
=
33/243
3.b Periodic excitations - errors-in-variables problem (contd)
Maximum likelihood estimate frequency response function
is independent (over ), circular complex normally distributed
Gaussian probability density function of the measurements
with
, and
Gaussian ML cost function
Eliminate ,
Z
m [ ]
Y
m [ ]
k ( ) U
m [ ]
k ( ) , [ ]
T
= m
Z
f
Z
Z ( ) f
Z
m [ ]
Z
m [ ]
( )
m 1 =
M
=
1
2M
det C
Z
( )
m 1 =
M
--------------------------------------------exp Z
m [ ]
Z
p
( )
H
C
Z
1
Z
m [ ]
Z
p
( )
m 1 =
M
( ) =
Z
p
T
G j
k
( ) , [ ]
T
= Z
p
Y
p
k ( )
U
p
k ( )
= C
Z
Cov Z
m [ ]
( )

Y
2
k ( )
YU
2
k ( )
YU
2
k ( )
U
2
k ( )
= =
V
ML
G j
k
( ) Z
p
Z , , ( ) Z
m [ ]
Z
p
( )
H
C
Z
1
Z
m [ ]
Z
p
( )
m 1 =
M
1 G j
k
( ) , [ ]Z
p
+ =
Z
p

G
ML
j
k
( )
Y
k ( )
U
k ( )
---------- -
M
1
Y
m [ ]
k ( )
m 1 =
M
M
1
U
m [ ]
k ( )
m 1 =
M
-------------------------------------------
= =
34/243
Maximum likelihood estimate frequency response function (contd)
Properties ML estimate FRF
- consistent ( )
- asymptotically circular complex normally distributed
- asymptotically efficient
- robust w.r.t. independence and normality assumptions
- expected value exists bias expression
- variance does not exist
Cramr-Rao lower bound
where
G
ML
j
k
( )
Y
k ( )
U
k ( )
---------- -
1
M
---- - Y
m [ ]
k ( )
m 1 =
M
1
M
----- U
m [ ]
k ( )
m 1 =
M
---------------------------------------
= =
M
var G
j
k
( ) ( ) G
0
j
k
( )
2

Y
2
k ( )
Y
0
k ( )
2
-----------------

U
2
k ( )
U
0
k ( )
2
------------------
2Re

Y
2
k ( )
Y
0
k ( )U
0
k ( )
------------------------- ( ) +

X
2
O
1
M
---- - ( ) =
35/243
Relative bias on FRF estimate
where
if
Asymptotic ( ) variance
with
b k ( )
G
j
k
( ) { }
G
0
j
k
( )
----------------------------
1 exp
U
0
k ( )
2
2
k ( )
------------------
( ) 1

Y
2
k ( )
k ( )
Y
k ( )
--------------------------
U
0
k ( )
U
k ( )
Y
0
k ( )
Y
k ( )
------------------------------

= =
b k ( ) 1
4
10 < min U
0
k ( )
U
k ( ) Y
0
k ( )
Y
k ( ) , ( ) 10 dB >
M
G
ML
j
k
( ) G
0
j
k
( )
G
k ( ) +
G
k ( ) G
0
j
k
( )
1
M
---- -
N
Y
m [ ]
k ( )
Y
0
k ( )
-----------------
m 1 =
M
1
M
---- -
N
U
m [ ]
k ( )
U
0
k ( )
-----------------
m 1 =
M

=
var
G
k ( ) ( ) G
0
j
k
( )
2

Y
2
k ( )
Y
0
k ( )
2
-----------------

U
2
k ( )
U
0
k ( )
2
------------------
2Re

Y
2
k ( )
Y
0
k ( )U
0
k ( )
------------------------- ( ) +

=
36/243
Uncertainty bounds
with
% confidence region = circle with radius
Notes:
- uncertainty bound accurate if
- radius
- decrease uncertainty FRF estimate: large for a given
- increase frequency resolution : small for a given
G
ML
j
k
( ) G
0
j
k
( ) AsN
c
0
G
2
k ( ) , ( )
G
2
k ( ) G
ML
j
k
( )
2

2
k ( )
Y
k ( )
2
---------------

2
k ( )
U
k ( )
2
----------------
2Re

2
k ( )
Y
k ( )U
k ( )
--------------------- ( ) +

=
p
1 p ( ) log
G
k ( )
U
0
k ( )
U
k ( ) 20 dB >
p 0.95 = 3
G
k ( )
M N
T
M N =
f
s
N M N
T
M N =
37/243
3.c Measurement example: octave bandpass filter
- multisine excitation , and
- points per period, periods
-120
-100
-80
-60
-40
-20
0
0 300 600 900 1200 1500 1800
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
-180
-90
0
90
180
0 300 600 900 1200 1500 1800
P
h
a
s
e

(
d
e
g
.
)
Frequency (Hz)
[
[
[
[ [
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[ [
[
[
[
[
[ [
[ [ [
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[ [
[
[
[
[
[
[
[
[ [
[ [
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[ [
[ [
[
[
[
[
[
[ [
[
[
[
[
[
[
[
[
[ [
[
[
[
[
[
[
[ [
[
[
[
[
[
[
[
[
[
[
[
[
[ [
[
[
[
[
[ [
[
[
[
[
[
[
[
[
[ [ [
[
[
[
[
[
[
[
[
[
[
[ [ [
[
[
[
[
[
[ [
[
[
[
[
[
[ [
[
[ [ [
[
[
[
[
[
[ [
[ [
[
[
-100
-80
-60
-40
-20
0
20
0 300 600 900 1200 1500 1800
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
[
[
[
[
[
[ [
[
[
[ [
[
[
[
[
[
[
[ [ [ [
[
[
[
[
[
[ [
[ [
[
[
[
[
[
[
[
[
[
[
[
[ [
[
[
[
[
[
[
[
[
[
[
[ [
[ [
[
[
[
[ [
[
[
[
[
[
[
[ [
[
[ [
[
[
[
[
[
[
[
[ [
[
[
[ [
[
[ [
[
[
[
[
[
[
[
[
[
[
[ [
[
[
[
[
[ [
[
[
[
[ [
[
[
[
[
[
[
[
[ [
[
[ [ [
[
[
[ [
[
[
[
[
[
[
[
[
[
[
[
[
[ [
[
[
[ [
[ [
[
[
[ [
[
[
[
[
[
[ [ [
[
[
[
[
[
[
[
[ [
[
[
[
[ [
[
[
[
[
[
[
-100
-80
-60
-40
-20
0
20
0 300 600 900 1200 1500 1800
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
kf
0
k 1 3 9 11 17 19 737 , , , , , , , = f
0
f
s
2
11
2.38 Hz =
N 2048 = M 8 =
38/243
3.d Arbitrary excitations - output error problem
Noisy output measurements, known input
Divide the input-output records in segments
Leakage reduction (Schoukens et al., 2006)
- as small as possible for a given
- (better than Hann window!)
Sample auto- and cross-power spectra
, ,
M
m 1 = 2

M
N
T
M N =
data points
N
data points
M N
T
Y
D
m [ ]
diff Y
m [ ]
( ) = U
D
m [ ]
diff U
m [ ]
( ) = ,
S
Y
D
k ( )
1
M
---- - Y
D
m [ ]
k ( )
2
m 1 =
M
= S
U
D
k ( )
1
M
----- U
D
m [ ]
k ( )
2
m 1 =
M
=
S
Y
D
U
D
k ( )
1
M
---- - Y
D
m [ ]
k ( )U
D
m [ ]
k ( )
m 1 =
M
=
39/243
3.d Arbitrary excitations - output error problem (contd)
Least squares estimator
Minimise w.r.t.
Sample variance
Asymptotic ( ) properties LS estimate FRF for known fixed input
- consistent
- asymptotically circular complex normally distributed
- asymptotically efficient
- expected value and variance exist
- not robust w.r.t. input measurement noise
V
LS
G j
k 1 2 +
( ) Z , ( ) Y
D
m [ ]
k ( ) G j
k 1 2 +
( )U
D
m [ ]
k ( )
2
m 1 =
M
=
G j
k 1 2 +
( )
G
LS
j
k 1 2 +
( )
S
Y
D
U
D
k ( )
S
U
D
k ( )
--------------------
=
Y
2
k 1 2 + ( )
M
2 M 1 ( )
---------------------- S
Y
D
k ( )
S
Y
D
U
D
k ( )
2
S
U
D
k ( )
-------------------------

=
M
40/243
3.d Arbitrary excitations - output error problem (contd)
Least squares estimator (contd)
Uncertainty bounds for fixed input
with
% confidence region = circle with radius
Notes:
- decrease leakage error: small for a given
- decrease uncertainty FRF estimate: large for a given
- frequency resolution FRF estimate
G
LS
j
k
( ) G
0
j
k
( ) AsN
c
0
G
2
k ( ) , ( )
G
2
k 1 2 + ( )

Y
D
2
k ( ) M
S
U
D
U
D
k ( )
------------------------
2
Y
2
k 1 2 + ( ) M
S
U
D
U
D
k ( )
-----------------------------------------
= =
p
1 p ( ) log
G
k ( )
M N
T
M N =
M N
T
M N =
f
s
N
41/243
3.e Simulation example: FRF and variance estimate using arbitrary input
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
-4
-2
0
2
4
i
n
p
u
t
number of points
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
-6
-4
-2
0
2
4
6
o
u
t
p
u
t
number of points
N 5000 = M , 2 =
42/243
3.e Simulation example: FRF estimate using arbitrary input (contd)
Note: variance estimate averaged over 30 frequencies
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-30
-25
-20
-15
-10
-5
0
5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
-30
-25
-20
-15
-10
-5
0
5
10
15
FRF
Variance output noise
estimate
true value
var(G)
43/243
3.f Nonlinear systems - introduction
Class of excitation signals with the same power spectrum:
- Gaussian noise
- periodic Gaussian noise
- random phase multisine
with
, and independent of
user defined, and
power is independent of
y
0
t ( )
y
S
t ( )
y
0
t ( ) y
R
t ( )
Gaussian noise
Wiener
Systems
u t ( )
G
BLA
j ( ) =
G
0
j ( ) G
B
j ( ) +
u t ( ) F
1 2 /
U
k
e
j 2kf
s
N ( )t
k F = k 0 ,
F
=
F N 2 < F N N
U
k
e
j U
k
{ } 0 =
N
1
u
2
nT
s
( )
n 0 =
N 1
N
44/243
3.f Nonlinear systems - introduction (contd)
Wiener systems:
- steady state response has same period as input (PISPO)
- saturation, discontinuities (e.g. relays) are allowed
- no chaos nor bifurcations
y
0
t ( )
y
S
t ( )
y
0
t ( ) y
R
t ( )
Gaussian noise
Wiener
Systems
u t ( )
G
BLA
j ( ) =
G
0
j ( ) G
B
j ( ) +
45/243
:
- true underlying linear system (if it exists)
:
- bias contribution
- depends on input power spectrum and odd nonlinear distortions only
- smooth function of the frequency
:
- best linear approximation or related linear dynamic system or LTI second
order equivalent
- can be approximated very well by a rational form in
y
S
t ( )
y
0
t ( ) y
R
t ( )
u t ( )
G
BLA
s ( ) =
G
0
s ( ) G
B
s ( ) +
Y
s
k ( )
Y
0
k ( )
U k ( )
G
BLA
j
k
( )
Y
R
k ( )
G
0
s ( )
G
B
s ( )
G
BLA
s ( )
s
46/243
:
- periodic signal
- depends on input power spectrum and the even and odd nonlinearities
- zero mean
- uncorrelated with - but not independent of -
- not normally distributed
y
S
t ( )
y
0
t ( ) y
R
t ( )
u t ( )
G
BLA
s ( ) =
G
0
s ( ) G
B
s ( ) +
Y
s
k ( )
Y
0
k ( )
U k ( )
G
BLA
j
k
( )
Y
R
k ( )
y
s
t ( )
u t ( )
47/243
:
- uncorrelated with but not independent of
-
- is asymptotically uncorrelated over (cumulant mixing over )
- asymptotically normally distributed
- is a smooth function of the frequency
Note: expectations are taken over the different random realisations of the excitation
y
S
t ( )
y
0
t ( ) y
R
t ( )
u t ( )
G
BLA
s ( ) =
G
0
s ( ) G
B
s ( ) +
Y
s
k ( )
Y
0
k ( )
U k ( )
G
BLA
j
k
( )
Y
R
k ( )
Y
s
k ( )
Y
s
k ( )U k ( ) { } 0 = Y
s
k ( ) U k ( )
Y
s
k ( ) { } 0 =
Y
s
k ( ) k k
Y
s
k ( )
2
{ }
48/243
3.g Nonlinear systems - noisy output measurements
:
- independent of
-
- asymptotically normally distributed
- asymptotically independent over
- is a smooth function of the frequency
Question
- distinction between and ?
G
BLA
s ( )
Gaussian noise
y
S
t ( )
y
R
t ( )
n
y
t ( )
y t ( ) y
0
t ( )
u
0
t ( )
N
Y
k ( )
U
0
k ( )
N
Y
k ( ) { } 0 =
k
N
Y
k ( )
2
{ }
y
s
t ( ) n
y
t ( )
49/243
3.g Nonlinear systems - noisy output measurements (contd)
FRF measurement using random phase multisine
FRF measurement using Gaussian noise
periods
transient
1 2
P
U
0
p [ ]
U
0
= Y
s
p [ ]
Y
s
= N
Y
p [ ]
, ,
G
ML
G
BLA
Y
s
U
0
N
Y
U
0
+ + =
var G
ML
( ) var N
Y
p [ ]
( ) P U
0
2
( ) =
segments
1 2
P
U
0
p [ ]
Y
s
p [ ]
N
Y
p [ ]
, ,
G
LS
G
BLA
S
Y
s
U
0
S
U
0
U
0
S
N
Y
U
0
S
U
0
U
0
+ + =
var G
LS
( ) var N
Y
p [ ]
( ) var Y
s
p [ ]
( ) + ( ) PS
U
0
U
0
( ) =
50/243
FRF measurement using random phase multisine (contd)
.
.
.
M

r
e
a
l
i
z
a
t
i
o
n
s
transient
P periods
...
...
...
G
m p , [ ]
Y
m p , [ ]
U
0
m [ ]
-------------
G
BLA
Y
S
m [ ]
U
0
m [ ]
----------
N
Y
m p , [ ]
U
0
m [ ]
--------------
+ + = =
G
1 [ ]
n
2 1 [ ]
,
G
2 [ ]
n
2 2 [ ]
,
G
M [ ]
n
2 M [ ]
,
G
ML

n
2
,
ML
2
G
1 1 , [ ]
G
1 2 , [ ]
G
1 P , [ ]
G
2 1 , [ ]
G
2 2 , [ ]
G
2 P , [ ]
G
2 P , [ ]
G
M 1 , [ ]
G
M 2 , [ ]
G
M P , [ ]
51/243
FRF measurement using random phase multisine (contd)
Sample mean and sample variance of each realisation
Sample mean and sample variance of the mean values
Properties
Note the asymmetric averaging of and
G
m [ ]
1
P
--- G
m p , [ ]
p 1 =
P
n
2 m [ ]
1
P P 1 ( )
--------------------- G
m p , [ ]
G
m [ ]
2
p 1 =
P
n
2
1
M
2
-------
n
2 m [ ]
m 1 =
M
= ,
G
ML
1
M
---- - G
m [ ]
m 1 =
M
ML
2
1
M M 1 ( )
------------------------ G
m [ ]
G
ML
2
m 1 =
M
n
2
{ }
var N
Y
( )
MP U
0
2
---------------------
=

ML
2
{ }
var Y
s
( )
M U
0
2
-----------------
var N
Y
( )
MP U
0
2
---------------------
+ =
Y
s
N
Y
52/243
3.h Simulation example: FRF estimate in the presence of nonlinear distortions
Multisine
- , , period length
-
Gaussian noise
- record length same measurement time
- segment length same averaging of the stochastic NL distortions
- same power spectrum
u t ( )
5tanh
x
5
---

x t ( ) z t ( )
Gaussian noise
1 2
1 2
P
1
2
M
multisine
M
M P N ( )
data points
n
y
t ( )
y t ( )
M 4 = P 2 = N 12500 =
std u t ( ) ( ) 1 =
M P N
N P
53/243
3.h Simulation example: FRF estimate in the presence of nonlinear distortions (contd)
Note:
- random noise is times more averaged for the multisine experiment
0 200 400 600 800 1000 1200 1400 1600 1800 2000
-60
-50
-40
-30
-20
-10
0
10
0 200 400 600 800 1000 1200 1400 1600 1800 2000
-60
-50
-40
-30
-20
-10
0
10
Random phase multisine Gaussian noise
FRF
total variance (noise + NL distortions)
noise variance
P 2 =
54/243
3.i Nonlinear systems - noisy input/output measurements
Noisy input/output measurements
Two cases
1. reference signal is not available
2. reference signal is available
Gaussian noise
G
BLA
s ( )
y
S
t ( )
y
R
t ( )
n
y
t ( )
y t ( ) y
0
t ( )
u
0
t ( )
n
u
t ( )
u t ( )
A s ( )
r t ( )
y t ( ) y
R
t ( ) y
S
t ( ) n
y
t ( ) + + =
u t ( ) u
0
t ( ) n
u
t ( ) + =
r t ( )
r t ( )
55/243
3.i Nonlinear systems - noisy input/output measurements (contd)
Case 1
- different realisations of periods each
- for each realisation: sample mean and sample variance input/output spectra
over the periods FRFs + noise variance (see section 3.b, p. 32-35)
- sample mean and sample variance FRFs over the realisations (see section
3.g, p. 51)
Case 2
- different realisations of periods each
- calculate the input, output, and reference DFT spectra of each period
- projection step
or
- variance analysis of the set of projected input/output DFT spectra
, ,
- calculate uncertainty on FRF (see section 3.b, p. 35)
M P
M P
U
R
m p , [ ] U
m p , [ ]
R
m [ ]
--------------
= Y
R
m p , [ ] Y
m p , [ ]
R
m [ ]
-------------
= , U
R
m p , [ ] U
m p , [ ]
e
j R
m [ ]
---------------
= Y
R
m p , [ ] Y
m p , [ ]
e
j R
m [ ]
---------------
= ,
U
R
2
var N
U
( )
Y
R
2
var N
Y
( ) var Y
s
( ) ,
Y
R
U
R
2
covar N
Y
N
U
, ( )
56/243
3.j Measurement example: the open loop gain of an operational amplifier
10
1
10
2
10
3
10
4
10
5
-80
-40
0
40
80
Frequency (Hz)
A
m
p
l
i
t
u
d
e

(
d
B
)
Open loop gain
10
1
10
2
10
3
10
4
10
5
-80
-40
0
40
80
Frequency (Hz)
A
m
p
l
i
t
u
d
e

(
d
B
)
Open loop gain
10
1
10
2
10
3
10
4
10
5
-100
-60
-20
20
Frequency (Hz)
P
h
a
s
e

(
)
Open loop gain
10
1
10
2
10
3
10
4
10
5
-100
-60
-20
20
Frequency (Hz)
P
h
a
s
e

(
)
Open loop gain
v
-
6.1 mVrms =
v
-
12.3 mVrms =
M 25 P , 5 = =
57/243
3.k Multivariable systems - periodic excitations
Difficulty: from one experiment one cannot calculate the frequency response matrix
Solution: perform experiments
with
:
:
, : -th input/output signal of the -th experiment
regular? FRM
Problem: which experiments guarantee that is regular and well conditioned?
U k ( )
G j
k
( )
Y k ( )
n
u
1
n
y
n
u
n
y
1
n
u
Y k ( ) G j
k
( )U k ( ) =
Y k ( ) n
y
n
u
U k ( ) n
u
n
u
U
p q , [ ]
k ( ) Y
p q , [ ]
k ( ) p q
U k ( )
G j
k
( ) Y k ( )U
1
k ( ) =
U k ( )
58/243
3.k Multivariable systems - periodic excitations (contd)
Guaranteed well conditioned MIMO experiments
1. cyclic shifts of multisines with interleaved frequency grid,
for example,
Notes:
- in theory one experiment is enough
- generators interact perform experiments
- loss in frequency resolution of a factor
- elimination interaction generators by precompensation
n
u
n
u
n
u
3 =
U k ( )
U
1
k ( ) U
3
k ( ) U
2
k ( )
U
2
k ( ) U
1
k ( ) U
3
k ( )
U
3
k ( ) U
2
k ( ) U
1
k ( )
=
U
1
U
2
U
3
n
u
n
u
59/243
Guaranteed well conditioned MIMO experiments (contd)
2. with an orthogonal matrix
,
for example,
Notes:
- maximal frequency resolution
- different power for each input
- nonlinear systems
with , ,
- elimination interaction generators by precompensation
U k ( ) U k ( )W = W
W
p q , [ ]
1
n
u
---------exp j
2 p 1 ( ) q 1 ( )
n
u
----------------------------------------- ( ) = p q , 1 2 n
u
, , , =
n
u
3 =
U k ( )
U k ( )
3
---------- -
1 1 1
1 e
j
2
3
------
e
j
2
3
------
1 e
j
2
3
------
e
j
2
3
------
=
U k ( ) U k ( )diag A
1
k ( ) A
2
k ( ) A
n
u
k ( ) , , , [ ] ( )W =
U k ( ) diag U
1
k ( ) U
2
k ( ) U
n
u
k ( ) , , , [ ] ( )Wdiag 1 e
j
2
k ( )
e
j
n
u
k ( )
, , , [ ] ( ) =
e
j
r
k ( )
{ } 0 = r 2 3 n
u
, , , = k 1 2 F , , , =
60/243
Variance analysis
- over consecutive periods disturbing noise
- over different realisations of full MIMO experiments stochastic NL
distortions + disturbing noise
61/243
3. Estimation parametric plant model with estimated
nonparametric noise model
3.l MIMO experiment: identification d-axis synchronous machine
3x220V 3x220V
H
P
1
4
3
0

A
/
D
-
c
o
n
v
e
r
t
o
r
H
P
1
4
4
5

A
W
-
g
e
n
e
r
a
t
o
r
H
P
1
4
3
0

A
/
D
-
c
o
n
v
e
r
t
o
r
H
P
1
4
3
0

A
/
D
-
c
o
n
v
e
r
t
o
r
H
P
1
4
3
0

A
/
D
-
c
o
n
v
e
r
t
o
r
H
P
1
4
4
5

A
W
-
g
e
n
e
r
a
t
o
r
M
X
I

c
o
n
t
r
o
l
l
e
r
Computer
VXI
measurement
system
Armature Field
Current
controller
Thyristor
rectifier
+
-
Current
controller
+
-
Thyristor
rectifier
i
a
i
f
e
f
e
a
62/243
3.l MIMO experiment: identification d-axis synchronous machine (contd)
-80
-60
-40
-20
0
0
.
0
0
5
0
.
0
1
0
.
11
1
0
1
0
0
3
0
0
I
a
(dBA)
f (Hz)
-80
-60
-40
-20
0
0
.
0
0
5
0
.
0
1
0
.
11
1
0
1
0
0
3
0
0
f (Hz)
I
f
(dBA)
-120
-100
-80
-60
-40
0
.
0
0
5
0
.
0
1
0
.
11
1
0
1
0
0
3
0
0
f (Hz)
E
a
(dBV)
-120
-100
-80
-60
-40
0
.
0
0
5
0
.
0
1
0
.
11
1
0
1
0
0
3
0
0
f (Hz)
E
f
(dBV)
M 8 =
N 65536 =
0.1 Hz 230 Hz , [ ]
periods
63/243
10
2
10
1
10
0
10
1
10
2
120
80
40
0
10
2
10
1
10
0
10
1
10
2
120
80
40
0
10
2
10
1
10
0
10
1
10
2
120
80
40
0
10
2
10
1
10
0
10
1
10
2
120
80
40
0
measured
noise variance
Z
11
dB ( )
Z
12
dB ( )
Z
22
dB ( )
f (Hz)
Z
21
dB ( )
f (Hz)
impedance
64/243
0.01 0.1 1 10 100
0
25
50
75
100
0.01 0.1 1 10 100
0
25
50
75
100
0.01 0.1 1 10 100
0
25
50
75
100
0.01 0.1 1 10 100
0
25
50
75
100
Z
11
( ) arg ( )
Z
12
( ) ( ) arg
Z
22
( ) ( ) arg
f (Hz)
Z
21
( ) ( ) arg
f (Hz)
65/243
3.m MIMO experiment: vibrating mechanical structure
u
1
u
2
r
1
r
2
y
1
y
2
n
u
n
y
2 = =
Shaker Al-plate
r
1
, r
2
: full random orthogonal multisines
66/243
3.m MIMO experiment: vibrating mechanical structure (contd)
0 200 400 600 800 1000
-100
-50
0
50
G11
0 200 400 600 800 1000
-100
-50
0
50
G12
0 200 400 600 800 1000
-100
-50
0
50
G21
0 200 400 600 800 1000
-100
-50
0
50
G22
FRF
total variance
noise variance
M 25 P , 100 = =
67/243
3.n Multivariable systems - random excitations
FRM measurement,
where
Noise variance measurement,
segments
1 2
M
data points
M N
M n
u
j
k 1 2 +
( ) S
Y
D
U
D
k ( )S
U
D
U
D
1
k ( ) =
S
XY
1
M
---- - X
m [ ]
Y
m [ ]H
m 1 =
M
=
M n
u
1 +
C
Y
k 1 2 + ( )
M
2 M n
u
( )
------------------------ S
Y
D
Y
D
k ( ) S
Y
D
U
D
k ( )S
U
D
U
D
1
k ( )S
Y
D
U
D
H
k ( ) ( ) =
Cov vec G
j
k 1 2 +
( ) ( ) ( )
1
M
---- -S
U
D
U
D
T
k ( ) C
Y
D
k ( )
2
M
-----S
U
D
U
D
T
k ( ) C
Y
k 1 2 + ( ) = =
68/243
3.o MIMO simulation example: discrete-time system
0 0.5 1 1.5 2
-40
-20
0
20
0 0.5 1 1.5 2
-20
-10
0
10
20
0 0.5 1 1.5 2
-40
-20
0
20
40
0 0.5 1 1.5 2
-20
0
20
40
0 0.5 1 1.5 2
-40
-20
0
20
40
0 0.5 1 1.5 2
-40
-20
0
20
40
true value
estimate
n
y
3 =
n
u
2 =
M n
u
1 + =
N 5000 =
FRM
69/243
3.o MIMO simulation example: discrete-time system (contd)
true value
estimate
0 1 2
-5
0
5
10
0 1 2
-10
-5
0
5
10
0 1 2
-5
0
5
10
0 1 2
-10
-5
0
5
10
0 1 2
-20
-10
0
10
20
0 1 2
-10
0
10
20
0 1 2
-5
0
5
10
0 1 2
-10
0
10
20
0 1 2
-10
0
10
20
n
y
3 =
n
u
2 =
M n
u
1 + =
N 5000 =
noise power spectrum
averaged over
30 frequencies
70/243
3.p References
Bendat, J. S. and A. G. Piersol (1980). Engineering Applications of Correlations and
Spectral Analysis. Wiley, New York.
Brillinger, D. R. (1981). Time Series: Data Analysis and Theory. McGraw-Hill, New York.
Dobrowiecki, T. and J. Schoukens (2007). Measuring a linear approximation to weakly
nonlinear MIMO systems, Automatica, vol. 43, no. 10, pp. 1737-1751.
Guillaume, P., I. Kollar and R. Pintelon (1996). Statistical Analysis of Nonparametric
Transfer Function Estimates, IEEE Trans. Instrum. Meas., vol. IM-45, no. 2, pp. 594-600.
Guillaume, P., R. Pintelon and J. Schoukens (1996). Accurate Estimation of Multivariable
Frequency Response Functions, Proceedings of the 13th IFAC Triennial World Conference,
San Fransisco (USA), June 30 - July 5, pp. 423-428.
Approach. IEEE Press: Piscataway, New Jersey (second edition).
Pintelon R. and J. Schoukens (2001). Measurement of frequency response functions using
periodic excitations, corrupted by correlated input/output errors, IEEE Trans. Instrum. and
Meas., vol. 50, no. 6, pp. 1753-1760.
Pintelon, R., Y. Rolain and W. Van Moer (2003). Probability density function for frequency
response function measurements using periodic signals, IEEE Trans. Instrum. and Meas.,
vol. 52, no. 1, pp. 61-68.
Schoukens, J., R. Pintelon, T. Dobrowiecki, and Y. Rolain (2005). Identification of linear
systems with nonlinear distortions, Automatica, vol. 41, no. 2, pp. 491-504.
Schoukens, J., Y. Rolain, and R. Pintelon (2006). Analysis of window/leakage effects in
frequency response function measurements. Automatica, vol. 42, no. 1, pp. 27-38.
71/243
4. Parametric models of LTI systems
4.a Band-limited versus zero-order-hold measurement set up 72
4.b Band-limited versus zero-order-hold signal assumption 77
4.c Measurement example: first order system 78
4.d Band-limited versus zero-order-hold noise models 83
4.e Lumped versus distributed LTI systems 91
4.f Lumped systems 93
4.g Relation between input/output DFT spectra - plant model only 95
4.h Simulation example: 4th order discrete-time Butterworth filter 98
4.i Simulation example: 2nd order continuous-time system 102
4.j Measurement example: octave bandpass filter 103
4.k Measurement example: reduction of iron 104
4.l Relation between input/output DFT spectra - plant and noise model 105
4.m Plant and noise model structures 106
4.n Parametrizations LTI systems 108
4.o Multivariable systems 111
4.p References 112
72/243
4.a Band-limited versus zero-order-hold measurement set up
Zero-order-hold (ZOH)
input known exactly output error problem
actuator is part of the model
perfect ZOH set up
- , from DC to absolute calibration
-
Actuator
y
1
t ( )
+
y kT
s
( )
m
y
t ( )
y
AA
t ( )
u t ( )
+
n
p
t ( )
Generator/
u
d
k ( )
ZOH
u
zoh
t ( )
G
y
s
(
)
G
c
s ( )
G
ZOH
z
1
( )
G s ( )
G
act
s ( )
controller
G
act
j ( ) 1 = G
y
j ( ) 1 =
G
ZOH
z
1
( ) 1 z
1
( )Z L
1
G s ( ) s { } { } =
73/243
4.a Band-limited versus zero-order-hold measurement set up (contd)
Band-limited (BL)
errors-in-variables problem
actuator is not modelled
perfect BL set up
- for relative calibration
- for
y
1
t ( ) u
1
t ( )
+
+
+
y kT
s
( )
n
g
t ( )
m
u
t ( ) m
y
t ( )
u
AA
t ( ) y
AA
t ( )
u
g
t ( )
+
n
p
t ( )
Generator
u
d
k ( )
ZOH
u
zoh
t ( )
G
u
s
(
)
G
y
s
(
)
u kT
s
( )
G s ( ) G
act
s ( )
G
u
j ( ) G
y
j ( ) = f f
s
2 <
G
u
j ( ) G
y
j ( ) 0 = = f f
s
2
74/243
Exact modelling: summary
Mixing intersample behaviour and model
- ZOH measurements and CT model
- BL measurements and DT model
Consequences?
- approximation errors increase model complexity
- results not portable to set ups with different intersample behaviour
BL ZOH
Two channel measurement Single channel measurement
Relative calibration Absolute calibration
No flat amplitude/linear phase
required
Flat amplitude/linear phase required
Instrument bandwidth Instrument bandwidth many times
Anti-alias filters required Anti-alias filtering not allowed
f
s
2 f
s
75/243
Mixing intersample behaviour and model (contd)
is a ZOH-signal
and
Although is not a ZOH-signal, there exists a discrete-time model that exactly
relates the input samples to the output samples
Notes:
- also depends on ( intersample behaviour )
- has (in general) a larger order than the original CT system
- identified with cannot be used in a set up with a different
r
zoh
t ( )
L s ( )
G
c
s ( )
u t ( )
y t ( )
r
zoh
t ( )
U z ( )
R z ( )
----------
1 z
1
( )Z{L
1
L s ( ) s { }} =
Y z ( )
R z ( )
--------- -
1 z
1
( )Z{L
1
L s ( )G s ( ) s { }} =
u t ( )
u nT
s
( ) y nT
s
( )
G z
1
( )
Y z ( )
U z ( )
----------
Y z ( ) R z ( )
U z ( ) R z ( )
------------------------
Z{L
1
L s ( )G s ( ) s { }}
Z{L
1
L s ( ) s { }}
---------------------------------------------------
= = =
G z
1
( ) L s ( ) u t ( )
G z
1
( ) G
c
s ( )
G z
1
( ) L s ( ) L s ( )
76/243
Properties models
- BL(cascade) = cascade(BL)
- ZOH(cascade) cascade(ZOH)
-20
-10
0
0 1 2
A
m
p
l
i
t
u
d
e

(
d
B
)
(1)
(2)
(a)
(b)
-10
-5
0
0 0.5
A
m
p
l
i
t
u
d
e

(
d
B
)
-360
-180
0
180
360
0 0.5
P
h
a
s
e

(
)
f/f
s
f/f
s
f/f
s
L j ( )
G j ( )
ZOH G j ( )L j ( ) ( )
ZOH G j ( ) ( )ZOH L j ( ) ( )
77/243
4.b Band-limited versus zero-order-hold signal assumption
-6
0
6
0 10 20
S
i
g
n
a
l
s
Time (s)
0
1
0 1 2 3 4
A
m
p
l
i
t
u
d
e
0
1
0 1 2 3 4

A
m
p
l
i
t
u
d
e
(a)
(b)
(c)
f/f
s
f/f
s
BL
ZOH
BL signal spectrum
ZOH signal spectrum
78/243
4.c Measurement example: first order system
y t ( )
R
C u t ( )
G s ( )
1
1 RCs +
--------------------
=
G
ZOH
z
1
( )
1 e
T
s
RC ( )
( )z
1
1 e
T
s
RC ( )
z
1
----------------------------------------
=
BL setup
ZOH setup
79/243
4.c Measurement example: BL set up, CT model (contd)
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-7
-6
-5
-4
-3
-2
-1
0
0 500 1000 1500 2000
amplitude (dB)
frequency (Hz)
X BL meas.
CT model
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-120
-90
-60
-30
0
0 500 1000 1500 2000
phase ()
frequency (Hz)
X BL meas.
CT model
X
X
XX
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
-0.05
-0.03
-0.01
0.01
0.03
0.05
0 500 1000 1500 2000
amplitude error (dB)
frequency (Hz)
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-0.5
-0.3
-0.1
0.1
0.3
0.5
0 500 1000 1500 2000
phase error ()
frequency (Hz)
G s ( )
1
b
1
s b
0
+
-------------------
=
80/243
4.c Measurement example: ZOH set up, DT model (contd)
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-7
-6
-5
-4
-3
-2
-1
0
0 500 1000 1500 2000
amplitude (dB)
frequency (Hz)
X ZOH meas.
DT model
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-120
-90
-60
-30
0
0 500 1000 1500 2000
phase ()
frequency (Hz)
X ZOH meas.
DT model
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
-0.05
-0.03
-0.01
0.01
0.03
0.05
0 500 1000 1500 2000
frequency (Hz)
X
X
X
X
X
X
X
X
X
XX
X
XX
X
X
X
X
X
X
-1
-0.5
0
0.5
1
0 500 1000 1500 2000
phase error ()
frequency (Hz)
G
ZOH
z
1
( )
a
0
z
1
1 b
0
z
1
+
---------------------
=
81/243
4.c Measurement example: mixing the assumptions (contd)
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-10
-8
-6
-4
-2
0
2
0 500 1000 1500 2000
amplitude (dB)
frequency (Hz)
X BL meas.
DT model
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-120
-90
-60
-30
0
0 500 1000 1500 2000
phase ()
frequency (Hz)
X BL meas.
DT model
G s ( )
1
b
1
s b
0
+
-------------------
=
G
ZOH
z
1
( )
a
0
z
1
1 b
0
z
1
+
---------------------
=
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-10
-8
-6
-4
-2
0
2
0 500 1000 1500 2000
amplitude (dB)
frequency (Hz)
X ZOH meas.
CT model
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-120
-90
-60
-30
0
0 500 1000 1500 2000
phase ()
frequency (Hz)
X ZOH meas.
CT model
82/243
4.c Measurement example: mixing the assumptions (contd)
X
X
XX
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
-0.05
-0.03
-0.01
0.01
0.03
0.05
0 500 1000 1500 2000
frequency (Hz)
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-0.5
-0.3
-0.1
0.1
0.3
0.5
0 500 1000 1500 2000
phase error ()
frequency (Hz)
X
XXX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
-0.05
-0.03
-0.01
0.01
0.03
0.05
0 500 1000 1500 2000
frequency (Hz)
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
-1
-0.5
0
0.5
1
0 500 1000 1500 2000
phase error ()
frequency (Hz)
ZOH setup, CT model
BL setup, DT model
G s ( ) e
s
a
2
s
2
a +
1
s 1 +
b
2
s
2
b +
1
s b
0
+
------------------------------------ =
G
ZOH
z ( ) z

a
2
a
1
z
1
z
2
+ +
b
2
b
1
z
1
b
0
z
2
+ +
------------------------------------------ =
83/243
4.d Band-limited versus zero-order-hold noise models
Problem statement
Noise mostly continuous-time:
- thermal noise
- 1/f noise in semi-conductor devices
-
Goal:
- write noise at output as filtered white noise
Reason:
- white noise source needed to prove asymptotic properties estimators
- nonparametric noise model covariance matrix with Toeplitz structure
input
plant
noise
output
N N
84/243
4.d Band-limited versus zero-order-hold noise models (contd)
First trial
- = white CT-noise
Problem
- either infinite variance or zero power spectrum
H
c
s ( )
e
c
t ( ) n t ( )
S
e
c
f ( )
f
0
S
n
f ( )
f
0
e
c
t ( )
85/243
Assumptions ZOH noise model
- = piecewise constant
- = white DT-noise
- no anti-alias filter
At sampling instances
H
c
s ( )
e
c
t ( )
n t ( )
T
s
n mT
s
( ) v m ( ) =
T
s
e
c
t ( )
e m ( ) e
c
mT
s
( ) =
H z
1
( ) 1 z
1
( )Z{L
1
H
c
s ( ) s { }} =
e t ( )
v t ( )
S
e
f ( )
f
0
S
v
f ( )
f
0
f
s
2 f
s
2
f
s
2
f
s
2
86/243
Assumptions BL noise model
- = band-limited white CT-noise (strm, 1970)
- is Gaussian distributed
- band-limited measurement setup
-
AA f ( )
f
0
f
s
2 f
s
2
1
S
e
c
f ( )
f
0
f
B
f
B
H
c
s ( )
e
c
t ( )
n t ( )
T
s
mT
s
( ) v m ( ) =
t ( )
AA s ( )
T
s
e
c
t ( )
e
c
t ( )
f
B
f
s
2
87/243
Proof
Notes:
1.
2. is independent, normally distributed
H s ( ) H
c
s ( ) =
t ( )
t ( )
S
f ( )
2
f
s
=
f
0
f
s
2
f
s
2
S
f ( )
f
0
f
s
2 f
s
2
T
s
mT
s
( ) v m ( ) =
S
f ( ) H
c
j2f ( )
2
AA j2f ( )
2
S
e
c
f ( ) H j2f ( )
2
S
f ( ) = =
R
( ) F
1
{S
f ( )}
2
sinc f
s
( ) = =
R
nT
s
( ) 0 = e n ( ) nT
s
( ) =
88/243
Relaxed assumption on the driving CT noise source
- = is constant in the band
- for and
where is the weakest decay such that
H
c
s ( )
e
c
t ( )
n t ( )
T
s
mT
s
( ) v m ( ) =
t ( )
AA s ( )
AA f ( )
f
0
f
s
2 f
s
2
1
S
e
c
f ( )
f
0
f
B
f
B
T
s
e
c
t ( )
S
e
c
f ( ) f f
B
<
S
e
c
f ( ) O f
1 + ( )
( ) = f f
B
0 >
O f
1 + ( )
( ) var e
c
t ( ) ( ) S
e
c
f ( ) f d

+
< =
89/243
Generalisation
Diffusion systems (Pintelon et al., 2005)
- plant model often rational form in :
- hence also true for the process noise model:
1/f noise
s G s ( )
H s ( )
H
c
s ( )
1
s
------
=
90/243
Summary
CT- versus DT-noise modelling
driving noise source noise model application
ZOH
(piecewise constant !?)
DT prediction, control
BL
(perfect AA filters !?)
CT physical interpretation
91/243
4.e Lumped versus distributed LTI systems
Lumped system
transfer function
poles
y t ( )
C
L
u t ( )
LC
d
2
y t ( )
dt
2
------------- y t ( ) + u t ( ) =
y 0 ( ) 0
dy t ( )
dt
-----------
t 0 =
, 0 = =
initial conditions:
G s ( )
Y s ( )
U s ( )
----------
1
1 LCs
2
+
----------------------
= =
s j LC =
92/243
4.e Lumped versus distributed LTI systems (contd)
Distributed system
transfer function
with
poles
,
Note
has limited active frequency range rational approximation in given frequency band
y x t , ( )
x l = x 0 =
u t ( )
2
y x t , ( )
t
2
-------------------
E
---
2
y x t , ( )
x
2
------------------- =
y x 0 , ( ) 0 =
y x t , ( )
t
-----------------
t 0 =
0 = ,
y 0 t , ( ) 0 =
y x t , ( )
x
-----------------
x l =
u t ( )
E
--------
= ,
initial conditions:
boundary conditions:
G s ( )
Y l s , ( )
U s ( )
--------------
l
E
---
tanh s ( )
s
-------------------
l
E
---
2
s ( )
2
2k 1 + ( ) 2 ( )
2
+
--------------------------------------------------------
k 0 =
= = =
l
2
E
------- =
s 2k 1 + ( )

2
----- j = k 1 2 , , =
2 s ( )
2
2k 1 + ( ) 2 ( )
2
+ ( )
s j =
93/243
4.f Lumped systems
Continuous-time (CT)
General

Sometimes (diffusion phenomena such as heat or mass transfer)

Discrete-time (DT)

a
n
y
n ( )
t ( )
n 0 =
n
a
b
m
u
m ( )
t ( )
m 0 =
n
b
= G s , ( )
B s , ( )
A s , ( )
----------------
b
r
s
r
r 0 =
n
b
a
r
s
r
r 0 =
n
a
------------------------
= =
a
n
d
n 2 /
y t ( )
dt
n 2 /
----------------- -
n 0 =
n
a
b
m
d
m 2 /
u t ( )
dt
m 2 /
-------------------
m 0 =
n
b
= G s , ( )
B s , ( )
A s , ( )
--------------------
b
m
s
m 2 /
m 0 =
n
b
a
n
s
n 2 /
n 0 =
n
a
--------------------------------
= =
a
n
y t n ( )
n 0 =
n
a
b
m
u t m ( )
m 0 =
n
b
= G z
1
, ( )
B z
1
, ( )
A z
1
, ( )
--------------------
b
r
z
r
r 0 =
n
b
a
r
z
r
r 0 =
n
a
--------------------------
= =
94/243
4.f Lumped systems (contd)
Summary
with
-domain s -domain s
-domain z
Im
Re
Im
Re Re
Im
G , ( )
B , ( )
A , ( )
------------------
b
r
r
r 0 =
n
b
a
r
r
r 0 =
n
a
--------------------------
= =
s continuous-time
s diffusion
z
1
discrete-time
=
95/243
4.g Relation between input/output DFT spectra - plant model only
Input/output DFT spectra
,
with
Periodic signals
if
- steady state response
- integer number of periods are observed
U k ( )
1
N
-------- u tT
s
( )z
k
t
t 0 =
N 1
= Y k ( )
1
N
-------- y tT
s
( )z
k
t
t 0 =
N 1
=
z
k
e
j2k N
=
Y k ( ) G
k
, ( )U k ( ) =
96/243
4.g Relation between input/output DFT spectra - plant model only (contd)
Arbitrary signals
with
,
where
= function of the difference between the initial and final conditions
,
Notes:
- when the initial conditions equal the final conditions
- reduce for by choice
Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( )
k
( ) + + =
G , ( )
B , ( )
A , ( )
------------------
b
r
r
r 0 =
n
b
a
r
r
r 0 =
n
a
--------------------------
= = T
G
, ( )
I
G
, ( )
A , ( )
-------------------
i
g
r
r
r 0 =
n
i
g
a
r
r
r 0 =
n
a
---------------------------
= =
i
g
r
n
i
g
max n
a
n
b
, ( ) 1 =
U k ( ) O N
0
( ) = T
G

k
, ( ) O N
1 2 /
( ) =

k
( )
0 z
1
=
O N
1 2 /
( ) s s , =
=
T
G
, ( )
k
( ) 0 = =

k
( ) s s , = n
i
g
max n
a
n
b
, ( ) 1 >
97/243
4.g Relation between input/output DFT spectra - plant model only (contd)
Summary
with
,
where
and
Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) + =
G , ( )
B , ( )
A , ( )
------------------
b
r
r
r 0 =
n
b
a
r
r
r 0 =
n
a
--------------------------
= = T
G
, ( )
I
G
, ( )
A , ( )
-------------------
i
g
r
r
r 0 =
n
i
g
a
r
r
r 0 =
n
a
---------------------------
= =
z
1
n
i
g
max n
a
n
b
, ( ) 1 =
s s , n
i
g
max n
a
n
b
, ( ) 1 >
=
T
G

k
, ( )
0 for periodic and time-limited signals
O N
1 2 /
( ) arbitrary signals
=
98/243
4.h Simulation example: 4th order discrete-time Butterworth filter
Cut off frequency
0
0.2
0.4
0.6
0.8
1
0 100 200 300
n (samples)
0
0.2
0.4
0.6
0.8
1
0 100 200 300
n (samples)
e
x
c
i
t
a
t
i
o
n

u
(
n
T
s
)
r
e
s
p
o
n
s
e

y
(
n
T
s
)
input output
f
s
4
99/243
4.h Simulation example: 4th order discrete-time Butterworth filter (contd)
Note:
- random behaviour error on FRF due to random behaviour
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[
[[
[
[
[
[[
[
[
[[
[
[
[
[
[
[
[
[[
[
[
[
[[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
-160
-120
-80
-40
0
40
0 0.1 0.2 0.3 0.4 0.5
a
m
p
l
i
t
u
d
e

(
d
B
)
f/fs
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
-180
-90
0
90
180
0 0.1 0.2 0.3 0.4 0.5
p
h
a
s
e

(
)
f/fs
amplitude phase
true value
[
Y k ( )
U k ( )
---------- -
T
G
z
k
1
, ( ) U k ( ) U k ( )
100/243
-350
-300
-250
-200
-150
-100
-50
0
0 0.1 0.2 0.3 0.4 0.5
e
r
r
o
r

(
d
B
)
f/fs
without transient
with transient
difference true value and model
101/243
-0.8
-0.6
-0.4
-0.2
0
0.2
0 50 100 150
t
r
a
n
s
i
e
n
t

r
e
s
p
o
n
s
e
n (samples)
impulse response
T
G
z
1
( )
I
G
z
1
( )
A z
1
( )
---------------
=
102/243
4.i Simulation example: 2nd order continuous-time system
true FRF
error model
without transient
error estimate
n
i
g
10 =
with transient
-120
-90
-60
-30
0
30
60
0 10 20 30 40 50 60 70
a
m
p
l
i
t
u
d
e

(
d
B
)
f (Hz)
103/243
4.j Measurement example: octave bandpass filter
plant model
n
a
6 n
b
, 4 = =
measured FRF
(multisine excit.)
error model with
transient
n
i
g
6 =
error model
error model
without transient
(multisine excit.)
(arbitr. excit.)
(arbitr. excit.)
X
X
X
X
X
X
X
X
X
XX
X
XX
X
X
X
X
X
XX
X
X
XX
XX
X
X
X
X
XX
X
XXX
XX
XX
XX
X
X
XXXXX
X
X
X
X
X
XXXXX
XXX
X
XX
XX
X
X
X
X
X
X
X
X
XXX
X
X
X
XX
X
X
X
X
X
X
X
XX
X
X
X
X
X
XXX
XX
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
X
XX
XX
X
X
X
XX
XX
X
X
XX
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
XX
X
X
XX
X
X
X
X
XXXXX
X
X
X
XX
X
X
X
XX
X
X
X
X
X
XX
X
X
X
XX
X
XX
XX
X
X
X
XX
X
XX
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
XX
X
X
X
X
XX
X
X
X
X
X
X
XXXXXX
X
X
X
XX
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XXX
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XXX
X
X
X
X
X
XX
XX
X
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
XX
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
XX
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XX
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
XXX
XX
X
X
XXX
-120
-100
-80
-60
-40
-20
0
20
0 0.5 1 1.5 2 2.5
a
m
p
l
i
t
u
d
e

(
d
B
)
f (kHz)
104/243
4.k Measurement example: reduction of iron
-30
-20
-10
0
0.01 0.1 1 10 100
P
h
a
s
e

(
)
Frequency (Hz)
-50
0
50
0.01 0.1 1 10 100
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
-50
0
50
0.01 0.1 1 10 100
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
-30
-20
-10
0
0.01 0.1 1 10 100
P
h
a
s
e

(
)
Frequency (Hz)
amplitude
phase
-domain model
without transient term
-domain model
with transient term
model
error
+ sine measurement
standard deviation sine measurement
n
a
n
b
3 = =
n
i
g
5 =
n
a
n
b
3 = =
s
s
105/243
4.l Relation between input/output DFT spectra - plant and noise model
where
, , ,
v t ( )
y t ( )
u t ( )
G , ( )
H , ( )
e t ( )
Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) H
k
, ( )E k ( ) T
H

k
, ( ) + + + =
G , ( )
B , ( )
A , ( )
------------------
= T
G
, ( )
I
G
, ( )
A , ( )
-------------------
= H
k
, ( )
C , ( )
D , ( )
-------------------
= T
H
, ( )
I
H
, ( )
D , ( )
-------------------
=
106/243
4.m Plant and noise model structures
Frequency domain
Model structures
ARX: , ,
ARMAX: , ,
ARMA: ,
BJ:
Notes:
- AR = autoregressive
- MA = moving average
- X = exogenous input
- BJ = Box-Jenkins
- hybrid BJ = CT plant + DT noise model
- AR or MA or ARMA parametric noise modelling = time series analysis or
spectral estimation (DT and CT)
Y k ( ) G
k
, ( )U k ( ) H
k
, ( )E k ( ) T
G

k
, ( ) T
H

k
, ( ) + + + =
C 1 = D A = T
H
0 =
D A = T
H
0 = n
i
g
max n
a
n
b
n
c
, , ( ) 1
G 0 = T
G
0 =
D A
107/243
4.m Plant and noise model structures (contd)
DT domain
with the backward shift operator:
and where
CT domain
with the derivative operator:
and where
y t ( ) G q , ( )u t ( ) H q , ( )e t ( ) + =
q z
1
= qx t ( ) x t 1 ( ) =
G q ( )x t ( ) g r ( )q
r
r 0 =
( )x t ( ) g r ( )x t r ( )
r 0 =
= =
y t ( ) G p , ( )u t ( ) H p , ( )e t ( ) + =
p
d
dt
-----
= px t ( )
dx t ( )
dt
-----------
=
G p ( )x t ( ) g ( )x t ( )d
0
=
108/243
4.n Parametrizations LTI systems
Plant model
Rational form
with
Partial fraction expansion
for
with
G , ( )
B , ( )
A , ( )
------------------
b
r
r
r 0 =
n
b
a
r
r
r 0 =
n
a
--------------------------
= =
T
a
0
a
1
a
n
a
b
0
b
1
b
n
b
[ ] =
G , ( )
L
r

r
---------------
r p =
r 0
p
S
r

r
----------------
r 1 =
q
W w , ( ) + + = s s , =
G z
1
, ( )
L
r
z
1
1
r
z
1
---------------------
r p =
r 0
p
S
r
z
1
1
r
z
1
---------------------
r 1 =
q
W z
1
w , ( ) + + =
q
Re
1
( )Im
1
( )Re
p
( )Im
p
( ) [ =
S
1
S
q
ReL
1
( )ImL
1
( )ReL
p
( )ImL
p
( )w
0
w
n
w
]
109/243
4.n Parametrizations LTI systems (contd)
Plant model (contd)
State space representation for proper transfer functions ( )
with
Pole/zero representation
disadvantage: ill conditioned for multiple poles/zeroes
Systems with time delay
n
b
n
a
G s , ( ) C sI
n
a
A ( )
1
B D + =
G z
1
, ( ) z
1
C I
n
a
z
1
A ( )
1
B D + =
T
vec
T
A ( ) B
T
C D [ ] =
G , ( ) K

r
( )
r 1 =
n
b

r
( )
r 1 =
n
a
------------------------------------ =
G , ( ) e
s
B , ( )
A , ( )
------------------ =
G z
1
, ( ) z
T
s

B z
1
, ( )
A z
1
, ( )
-------------------- =
110/243
4.n Parametrizations LTI systems (contd)
Noise model
Similar as plant model
Note
Starting value problem parametrizations (see chapter 6)
- rational form linear least squares
- state space subspace methods
- other parametrizations via rational form + back transformation
111/243
4.o Multivariable systems
Left matrix fraction (LMF) description
with , and
Right matrix fraction (RMF) description
with , and
Common denominator
with , and
State space
same form as SISO with , , , and
Partial fraction expansion
are residue matrices
Notes:
- advantage LMF:
- common denominator: rank residue matrices cannot be imposed
( too many parameters if residue matrices are not of full rank)
- LMF, RMF, and state space: rank one residue matrices
( with )
G , ( ) A
1
, ( )B , ( ) = A n
y
n
y
B n
y
n
u
G , ( ) B , ( )A
1
, ( ) = A n
u
n
u
B n
y
n
u
G , ( ) B , ( ) A , ( ) = A 1 1 B n
y
n
u
A n n B n n
u
C n
y
n D n
y
n
u
L
r
S
r
, n
y
n
u
Y k ( ) G , ( )U k ( ) = A , ( )Y k ( ) B , ( )U k ( ) =
F n n F
1
adj F ( ) det F ( ) = det adj F ( ) ( ) det F ( ) ( )
n 1
=
112/243
4.p References
Box, G. E. P. and G. M. Jenkins (1976).Time Series Analysis: Forecasting and Control.
Holden-Day, Oakland.
Hannan, E. J. and M. Deistler (1988). Linear Systems. Wiley: New York.
Henrici, P. (1974). Applied and Computational Complex Analysis. Wiley: New York, vol. 1.
Kailath, T. (1980). Linear Systems. Prentice Hall, Englewood Cliffs.
Ljung, L. (1999). System Identification: Theory for the User. Prentice Hall: Upper Saddle
River, New Jersey (second edition).
Oldham, K. B., and J. Spanier (1974). The Fractional Calculus. Academic Press: New York.
Peeters F., R. Pintelon, J. Schoukens and Y. Rolain (2001). Identification of Rotor-Bearing
Systems in the Frequency Domain. Part II: Estimation of Modal Parameters, Mechanical
Systems and Signal Processing, vol. 15, no. 4, pp. 775-788.
Pintelon R. and J. Schoukens (1997). Identification of Continuous-Time Systems Using
Arbitrary Signals, Automatica, vol. 33, no. 5, pp. 991-994.
Pintelon R., J. Schoukens and G. Vandersteen (1997). Frequency Domain System
Identification Using Arbitrary Signals, IEEE Trans. Autom. Contr., vol. AC-42, no. 12, pp.
1717-1720.
Pintelon R., J. Schoukens, L. Pauwels, and E. Van Gheem (2005). Diffusion systems:
stability, modeling and identification, IEEE Trans. Instrum. Meas., vol. 54, no. 5, pp. 2061-
2067.
Sderstrm, T. and P. Stoica (1989). System Identification. Prentice Hall: Englewood Cliffs.
113/243
5. Improved nonparametric models for LTI systems
5.a Arbitrary excitations generalised output error 114
5.b Simulation example: comparison LPM and SA 121
5.c Measurement example: nonlinear electrical circuit 123
5.d Arbitrary excitations errors-in-variables 125
5.e Measurement example: flexural vibrations of a steel beam 127
5.f Periodic excitations steady state 130
5.g Measurement example: longitudinal vibrations of a plexiglass beam 140
5.h Periodic excitations transient regime 141
5.i Measurement example: flexural vibrations of a steel beam 142
5.j Summary local polynomial methods 144
5.k References 145
114/243
5.a Arbitrary excitations generalised output error
Goal: starting from and obtain a nonparametric estimate of
- the frequency response function
- the noise covariance matrix
Difficulties: the presence of
- for estimating
- and for estimating
Property:
- and are smooth functions of the frequency
G ( )
u t ( )
y t ( )
v t ( )
H ( )
e t ( )
Y k ( ) G
k
( )U k ( ) T
k
( ) V k ( ) + + =
V k ( ) H
k
( )E k ( ) =
T
k
( ) T
H

k
( ) T
G

k
( ) + =
U k ( ) Y k ( ),
G ( )
C
V
k ( ) Cov V k ( ) ( ) =
T
k
( ) G
k
( )
T
k
( ) G
k
( )U k ( ) C
V
k ( )
G ( ), H ( ), T ( )
115/243
5.a Arbitrary excitations generalised output error (contd)
Basic idea solution: local polynomial approximation of and
with Taylor series expansion
and
G ( )
u t ( )
y t ( )
v t ( )
H ( )
e t ( )
Y k ( ) G
k
( )U k ( ) T
k
( ) V k ( ) + + =
V k ( ) H
k
( )E k ( ) =
T
k
( ) T
H

k
( ) T
G

k
( ) + =
G ( ) T ( )
Y k r + ( ) G
k r +
( )U k r + ( ) T
k r +
( ) V k r + ( ) + + =
G
k r +
( ) G
k
( ) g
s
k ( )r
s
s 1 =
R
O r N ( )
R 1 + ( )
( ) + + =
T
k r +
( ) T
k
( ) t
s
k ( )r
s
s 1 =
R
N
1 2 /
O r N ( )
R 1 + ( )
( ) + + =
G
k
( ) g
s
k ( )
t
s
k ( ) T
k
( )
r n n 1 + 1 0 1 n 1 n , , , , , , , , =
G
k
( ) g
1
k ( ) g
2
k ( ) g
R
k ( ) T
k
( ) t
1
k ( ) t
2
k ( ) t
R
k ( )
=
116/243
Local polynomial model in the band :
where
and
with
G ( )
u t ( )
y t ( )
v t ( )
H ( )
e t ( )
Y k ( ) G
k
( )U k ( ) T
k
( ) V k ( ) + + =
V k ( ) H
k
( )E k ( ) =
T
k
( ) T
H

k
( ) T
G

k
( ) + =
k n k n + , [ ]
Y
n
K
n
V
n
+ =
X
n X k n ( ) X k n 1 + ( ) X k ( ) X k n + ( )
=
K k r + ( )
K
1
r ( ) U k r + ( )
K
1
r ( )
= K
1
r ( )
1
r
r
R
=
117/243
Local polynomial model
Linear least squares estimate
Residual linear least squares estimate
with
the degrees of freedom of the residuals
G ( )
u t ( )
y t ( )
v t ( )
H ( )
e t ( )
Y k ( ) G
k
( )U k ( ) T
k
( ) V k ( ) + + =
V k ( ) H
k
( )E k ( ) =
T
k
( ) T
H

k
( ) T
G

k
( ) + =
Y
n
K
n
V
n
+ =
Y
n
K
n
H
K
n
K
n
H
( )
1
= G

k
( )
V
n
Y
n

K
n
= C
V
k ( )
1
q
---V
n
V
n
H
=
q 2n 1 R 1 + ( ) n
u
1 + ( ) + =
dof
118/243
Repeat the previous procedure at frequency for all the other frequencies , ,
Consequence
the FRF and noise covariance estimates are correlated over the
frequency
the correlation length is
k k 1 + k 2 +
k 1 + k k n 1 k n k n 1 +

k 1

k n + k n 1 + + k n 2 + +

k
( ) C
V
k ( )
2n
119/243
Bias error local polynomial method (LPM)
Bias error spectral analysis approach
with for win = diff, half sine; and for win = hann
G ( )
u t ( )
y t ( )
v t ( )
H ( )
e t ( )
Y k ( ) G
k
( )U k ( ) T
k
( ) V k ( ) + + =
V k ( ) H
k
( )E k ( ) =
T
k
( ) T
H

k
( ) T
G

k
( ) + =
G

k
( ) { } G
0

k
( ) G
0
R 1 + ( )

k
( )O
intG
n N ( )
R 1 + ( )
( ) O
leakG
n N ( )
R 2 + ( )
( ) + + =
C
V
k ( ) { } C
V
k ( ) C
V
1 ( )
k ( )O
intH
n N ( ) G
0
R 1 + ( )

k
( )O
intG
n N ( )
2 R 1 + ( )
( ) G
0
R 1 + ( )

k
( ) ( )
H
+ + =
G
win
k
( ) { } G
0

k
( )
win
G
0
1 ( )

k
( )O
intG
M N ( )
2
( ) O
leakG
M N ( )
2
( ) + + =
C
V
win
k ( ) { } C
V
k ( )

win
2
----------C
V
2 ( )
k ( )O
intH
M N ( )
2
( )

win
2
----------G
0
1 ( )

k
( )O
intG
M N ( )
2
( ) G
0
1 ( )

k
( ) ( )
H
+ + =
win
1 4 =
win
1 3 =
120/243
Local polynomial approach Spectral analysis method
smaller bias error larger bias error
smaller uncertainty larger uncertainty
larger frequency resolution smaller frequency resolution
larger calculation time smaller calculation timev
correlation length correlation length (rect),
(diff, half sine), or (hanning)
G ( )
u t ( )
y t ( )
v t ( )
H ( )
e t ( )
Y k ( ) G
k
( )U k ( ) T
k
( ) V k ( ) + + =
V k ( ) H
k
( )E k ( ) =
T
k
( ) T
H

k
( ) T
G

k
( ) + =
2n 0
1 2
121/243
5.b Simulation example: comparison LPM and SA
0 0.5 1 1.5 2
100
50
0
50
Frequency (Hz)
F
R
F

(
d
B
)
0 0.5 1 1.5 2
100
50
0
50
Frequency (Hz)
F
R
F

(
d
B
)
0 0.5 1 1.5 2
100
50
0
50
Frequency (Hz)
F
R
F

(
d
B
)
0 0.5 1 1.5 2
100
50
0
50
Frequency (Hz)
F
R
F

(
d
B
)
G
0
G
SA
G
0
LPM
G
0
LPM and SA (hanning window): same dof noise covariance estimates

122/243
5.b MIMO simulation example: comparison LPM and SA (contd)
LPM and SA (hanning window): same dof noise covariance estimates
0 0.5 1 1.5 2
60
40
20
0
20
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
0 0.5 1 1.5 2
40
30
20
10
0
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
0 0.5 1 1.5 2
40
30
20
10
0
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
0 0.5 1 1.5 2
60
40
20
0
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
true value
SA
LPM
123/243
5.c Measurement example: nonlinear electrical circuit
White noise excitation with bandwidth of
samples at :
split in 20 blocks of samples each
spectral analysis and local polynomial estimates on each block of samples
sample means and sample variances over the 20 blocks
y
3
t ( )
u t ( )
y t ( )
L 1 H =
C 9.4 F = R 16 =
1000 V
2
200 Hz
10400 f
s
20 MHz 2
15
600 Hz =
N 520 =
N 520 =
124/243
5.c Measurement example: nonlinear electrical circuit (contd)
0 50 100 150 200 250
80
60
40
20
0
20
Freq. (Hz)
F
R
F

(
d
B
)
0 50 100 150 200 250
120
100
80
60
40
20
Freq. (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
hann, diff
LPM
hann, diff
LPM
BJ
BJ
full lines: FRF
dashed lines: std FRF
full lines: noise var.
dashed lines: std noise var.
125/243
5.d Arbitrary excitations errors-in-variables
Indirect method for measuring the FRF (two possible cases):
Linear plant and (non)linear actuator and controller:
Nonlinear plant and linear actuator and controller:
Solution:
Reference signal needed + modeling from to
n
c
t ( )
Plant
n
p
t ( )
r t ( )
m
y
t ( )
y t ( )
m
u
t ( )
u t ( )
Actuator
n
g
t ( )
Controller
G
0
j
k
( ) S
yr
j
k
( )S
ur
1
j
k
( ) =
G
BLA
j
k
( ) S
yr
j
k
( )S
ur
1
j
k
( ) =
r t ( ) r t ( ) z t ( ) y
T
t ( ) u
T
t ( ) [ ]
T
=
126/243
5.d Arbitrary excitations errors-in-variables (contd)
Modeling from reference to input and output:
with

Note: applicable to both spectral analysis and local polynomial methods
n
c
t ( )
Plant
n
p
t ( )
r t ( )
m
y
t ( )
y t ( )
m
u
t ( )
u t ( )
Actuator
n
g
t ( )
Controller
Z k ( )
G
ry

k
( )
G
ru

k
( )
R k ( )
T
H
Y

k
( )
T
H
U

k
( )
V
Z
k ( ) + + = Z k ( )
Y k ( )
U k ( )
=
G

k
( ) G
ry

k
( )G
ru
1
k
( ) =
127/243
5.e Measurement example: flexural vibrations of a steel beam
Beam parameters: length 61 cm, height 2.47 cm, width 4.93 mm, density 7800 kg/m
3
Excitation in the band [0 Hz, 6 kHz]:
random binary sequence
random phase multisine
samples of the reference, force (input), and acceleration (output) signals at
:
spectral analysis and local polynomial estimates
generalised output error and errors-in-variables estimates
u t ( )
0 x
0
L
y t x , ( )
x
N 50 1024 =
f
s
10 MHz 2
9
19.53 kHz =
128/243
5.e Measurement example: flexural vibrations of a steel beam (contd)
0 2 4 6
-50
0
50
Frequency (kHz)
F
R
F

(
d
B
)
0 2 4 6
-50
0
50
Frequency (kHz)
F
R
F

(
d
B
)
0.92 0.93 0.94 0.95
40
60
80
Frequency (kHz)
F
R
F

(
d
B
)
0.92 0.93 0.94 0.95
40
60
80
Frequency (kHz)
F
R
F

(
d
B
)
OE-LPM EIV-LPM
random binary sequence
129/243
5.e Measurement example: flexural vibrations of a steel beam (contd)
4.12 4.14 4.16 4.18 4.2
-60
-40
-20
0
20
Frequency (kHz)
F
R
F

(
d
B
)
0.16 0.18 0.2
-50
0
50
Frequency (kHz)
F
R
F

(
d
B
)
+
local polynomial approach
multisine
spectral analysis
FRF
std(FRF)
o
local polynomial approach
multisine
spectral analysis
130/243
5.f Periodic excitations steady state
Two possible cases:
Linear plant and (non)linear actuator and controller:
Nonlinear plant and linear actuator and controller:
Required signals:
(for multi-input and/or nonlinear plant only), and
n
c
t ( )
Plant
n
p
t ( )
r t ( )
m
y
t ( )
y t ( )
m
u
t ( )
u t ( )
Actuator
n
g
t ( )
Controller
G
0
j
k
( ) S
yr
j
k
( )S
ur
1
j
k
( ) =
G
BLA
j
k
( ) S
yr
j
k
( )S
ur
1
j
k
( ) =
r t ( ) u t ( ) y t ( )
131/243
5.f Periodic excitations steady state (contd)
with
Noise transients (leakage) errors
Smooth functions of the frequencies
Increase the variability of the FRF measurement
Introduce a correlation between consecutive signal periods
Important for lightly damped systems (e.g. mechanical vibrating structures)
Z k ( ) Z
0
k ( ) H
Z

k
( )E
Z
k ( ) T
H
Z

k
( ) + + = Z k ( )
Y k ( )
U k ( )
=
T
H
Z

k
( ),
132/243
1 2 3
DFT bin 1 2 3
DFT bin
noiseless signal noisy signal
signal noise transient noise
P 2 = DFT spectrum of consecutive periods
Z k ( ) Z
0
k ( ) H
Z

k
( )E
Z
k ( ) T
H
Z

k
( ) + + = H
Z

k
( )E
Z
k ( ) T
H
Z

k
( )
133/243
Basic idea solution: local polynomial expansion at the non-excited DFT lines
with Taylor series expansion for ( is not used!)
and
kP 5 + kP 3 + kP 1 + kP 3 kP 1
kP
DFT bin
kP 2 + kP 4 + kP 2
T
H
Z
( )
Z kP m
i
( ) H
Z

kP m
i
( )E
Z
kP m
i
( ) T
H
Z

kP m
i
( ) + =H
Z

kP m
i
( )E
Z
kP m
i
( ) T
H
Z

kP m
i
( )
i 1 2 n , , , = kP
T
H
Z

kP m
i
( ) T
H
Z

kP
( ) t
r
k ( ) m
i
( )
r
r 1 =
R
1
PN
------------O
m
i
PN
--------

R 1 + ( )
( ) + + =T
H
Z

kP
( ) t
r
k ( )
T
H
Z

k
( ) t
1
k ( ) t
2
k ( ) t
R
k ( )
=
134/243
Local polynomial model in the band
where
and
kP 5 + kP 3 + kP 1 + kP 3 kP 1
kP
DFT bin
kP 2 + kP 4 + kP 2
kP m
i
, i 1 2 n , , , =
Z
n
K
n
V
n
+ =
X
n
X kP m
n
( ) X kP m
1
( ) X kP m
1
+ ( ) X kP m
n
+ ( )
=
K kP m
i
( )
1 m
i
( ) m
i
( )
R
T
=
135/243
Local polynomial model
Linear least squares estimate
Residual linear least squares estimate
with the degrees of freedom of the residuals
kP 5 + kP 3 + kP 1 + kP 3 kP 1
kP
DFT bin
kP 2 + kP 4 + kP 2
Z
n
K
n
V
n
+ =
Z
n
K
n
H
K
n
K
n
H
( )
1
=
T
H
z

kP
( )
V
n
Z
n

K
n
= C
V
k ( )
1
q
---V
n
V
n
H
=
q 2n R 1 + ( ) = dof
136/243
At the excited DFT lines
Sample mean
Noise sample covariance sample mean
where typically giving a 1 dB increase in uncertainty
kP 5 + kP 3 + kP 1 + kP 3 kP 1
kP
DFT bin
kP 2 + kP 4 + kP 2
kP
Z kP ( ) Z
0
kP ( ) V kP ( ) T
H
Z

kP
( ) + + = V kP ( ) T
H
Z

kP
( )
Z
kP ( ) Z kP ( ) T
H
Z

kP
( ) =
C
kP ( ) 1
2
+ ( )C
V
kP ( ) =
0.5 =
137/243
Robust method - linear plant
as many experiments as inputs with orthogonal multisines ( periods)
- transient removal via LPM at non-excited lines
- sample means + noise covariance estimate
FRF estimate + noise covariance
Robust method - nonlinear plant
as many experiments as inputs with full random orthogonal multisines ( periods)
repeat the experiments for independent random phase realisations
- project input-output DFT spectra on reference DFT spectra
- sample mean and total sample covariance (noise + NL distortions)
- FRF estimate + noise and total covariances
P 2
G

k
( ) Y
k ( )U
1
k ( ) =
P 2
M 2
G

k
( ) Y
R
k ( )U
R
1
k ( ) =
138/243
Fast method - (non)linear plant
one experiment with random phase multisines ( periods)
apply LPM to the sample means at the excited frequencies
- local polynomial approximation FRF (note: no transient term!)
- FRF estimate + noise and total covariances
P 2
Z
kP ( )
G
ry

kP
( )
G
ru

kP
( )
R k ( ) V
Z
kP ( ) + =
G

kP
( ) G
ry

kP
( )G
ru
1
kP
( ) =
139/243
Comparision robust and fast methods
Fast Robust
Local polyn. approx. transient + FRF transient only
1 experiment experiments
Measurements periods periods
realisations
Frequency resolution
( = total measurement time)
Note: factor loss in frequency resolution of the fast method w.r.t. arbitrary excitations
n
u
P 2 P 2
M 2
1
P T
-------------
1
n
u
P M T
------------------------------------
T
P
140/243
5.g Measurement example: longitudinal vibrations of a plexiglass beam
0 1 2 3 4
100
50
0
Frequency (kHz)
F
R
F

(
d
B
)
0 1 2 3 4
2
0
2
4
6
Frequency (kHz)
v
a
r
r
e
c
t
/
v
a
r
L
P
M

(
d
B
)
u t ( )
0 x
L 1.983 m =
y t x , ( )
measured FRF total variance noise variance
Force-to-acceleration M 25 P 10 = , =
plexiglass beam
1 cm 2 cm 1.983 m
141/243
5.h Periodic excitations transient regime
Conclusion
The robust and fast methods for periodic excitations are robust to plant transient
(leakage) errors
n
c
t ( )
Plant
n
p
t ( )
r t ( )
m
y
t ( )
y t ( )
m
u
t ( )
u t ( )
Actuator
n
g
t ( )
Controller
Y k ( )
U k ( )
G
ry

k
( )
G
ru

k
( )
R k ( )
T
G
Y

k
( )
T
G
U

k
( )
H
Y

k
( )E
Y
k ( )
H
U

k
( )E
U
k ( )
T
H
Y

k
( )
T
H
U

k
( )
+ + + =
T
G
Y

k
( )
T
G
U

k
( )
T
H
Y

k
( )
T
H
U

k
( )
142/243
5.i Measurement example: flexural vibrations of a steel beam
Experiment (see also p. 127)
- 1 experiment with a random phase multisine excitation in the band [0 Hz, 6 kHz]
- consecutive periods
Processing
1. subtract the last period from the other signal periods
2. calculate the rms value of the residuals over each period
P 50 =
0 10 20 30 40 50
60
40
20
0
Period number
r
m
s

r
e
s
i
d
u
a
l

(
d
B
)
output
input
143/243
5.i Measurement example: flexural vibrations of a steel beam (contd)
0 2 4 6
50
0
50
Frequency (kHz)
F
R
F

(
d
B
)
0 2 4 6
5
0
5
10
15
Frequency (kHz)
v
a
r
r
e
c
t
/
v
a
r
L
P
M

(
d
B
)
0 2 4 6
50
0
50
Frequency (kHz)
F
R
F

(
d
B
)
0 2 4 6
0
10
20
30
40
Frequency (kHz)
v
a
r
r
e
c
t
/
v
a
r
L
P
M

(
d
B
)
last 36 periods
first 36 periods
G
LPM
var(G
LPM
)
var(G
rect
)
var(G
rect
)/var(G
LPM
)
144/243
5.j Summary local polynomial methods
Arbitrary excitations
full frequency resolution (1 experiment)
no distinction between noise and nonlinear distortions
local polynomial approximation of the FRF and the transient (leakage) errors
Periodic excitations - fast method
1/2 full frequency resolution (1 experiment, 2 periods)
separation between noise and nonlinear distortions
local polynomial approximation of the FRF and the transient (leakage) errors
Periodic excitations - robust method
full frequency resolution ( experiments, 2 periods, 2 realisations)
separation between noise and nonlinear distortions
local polynomial approximation of the transient (leakage) errors only
1 4n
u
( ) n
u
145/243
5.k References
Gevers, M., R. Pintelon, and J. Schoukens (2011). The local polynomial method for
nonparametric system identification: improvements and experimentation. IEEE Conf. on
Decision and Contr. and European Contr. Conf., Orlando-Florida (USA), Dec. 12-15.
Pintelon, R., J. Schoukens, G. Vandersteen, and K. Barb (2010a). Estimation of
nonparametric noise and FRF models for multivariable systemsPart I: theory. Mechanical
Systems and Signal Processing, vol. 24, no. 3, pp. 573595.
Pintelon, R., J. Schoukens, G. Vandersteen, and K. Barb (2010b). Estimation of
nonparametric noise and FRF models for multivariable systemsPart II: extensions,
applications. Mechanical Systems and Signal Processing, vol. 24, no. 3, pp. 596616.
Pintelon, R., K. Barb, G. Vandersteen, and J. Schoukens (2011). Improved (non-
)parametric identification of dynamic systems excited by periodic signals. Mechanical
Systems and Signal Processing, vol. 25, no. 7, pp. 2683-2704.
Pintelon, R., G. Vandersteen, J. Schoukens, and Y. Rolain (2011). Improved (non-
)parametric identification of dynamic systems excited by periodic signals - The multivariate
case. Mechanical Systems and Signal Processing, vol. 25, no. 8, pp. 28922922.
Schoukens, J., G. Vandersteen, K. Barb, and R. Pintelon (2009). Nonparametric
preprocessing in system identification: a powerful tool, European Journal of Control, vol. 15,
no. 3-4, pp. 260-274.
Schoukens, J., G. Vandersteen, and R. Pintelon (2010). Frequency response function
measurements using concatenated data records. International Conference on Noise and
Vibration Engineering, 20-22 Sept., Leuven (Belgium), pp. 3423-3430.
146/243
6. Estimation parametric plant model with known noise model
6.a Frequency domain data 147
6.b Maximum likelihood estimator 149
6.c Calculation ML estimates 155
6.d Calculation covariance matrix ML estimates 157
6.e Starting values 158
6.f Simulation example: 2nd order continuous-time system 160
6.g Measurement example: q-axis impedance 3.4 MW synchronous machine 164
6.h Measurement example: flight flutter data 165
6.i Model validation 166
6.j Measurement example: flexible robot arm 168
6.l Model selection 170
6.m Simulation example: SISO-FIR system 172
6.k Multivariable systems 169
6.n References 173
147/243
6.a Frequency domain data
Noisy input/output observations
with known
- periodic signals
- arbitrary signals
N
g
k ( )
U
g
k ( )
M
U
k ( )
U k ( ) Y k ( )
M
Y
k ( )
N
p
k ( )
G
0
j ( )
N
g
k ( ) generator noise
N
p
k ( ) process noise
M
U
k ( ) input measurement noise
M
Y
k ( )
output measurement noise
Y k ( ) Y
0
k ( ) N
Y
k ( ) + =
U k ( ) U
0
k ( ) N
U
k ( ) + =

Y
2
k ( ) var N
Y
k ( ) ( ) N
Y
k ( )
2
{ } = =
U
2
k ( ) var N
U
k ( ) ( ) N
U
k ( )
2
{ } = =
YU
2
k ( ) covar N
Y
k ( ) N
U
k ( ) , ( ) N
Y
k ( )N
U
k ( ) { } = =
U
0
k ( ) U
g
k ( ) =
U
0
k ( ) U
g
k ( ) N
g
k ( ) + =
148/243
6.a Frequency domain data (contd)
Note:
FRF measurements
with known
are a special case of noisy input/output observations
, , , and
G j
k
( ) G
0
j
k
( ) N
G
k ( ) + =
G
2
k ( ) var N
G
k ( ) ( ) N
G
k ( )
2
{ } = =
Y k ( ) G k ( ) = U k ( ) 1 = N
Y
k ( ) N
G
k ( ) = N
U
k ( ) 0 =
149/243
6.b Maximum likelihood estimator
Define
and
Assumption
has zero mean and is independent (over ), circular complex normally
distributed
Gaussian likelihood function
with
and
and where can be singular (e.g. feedback with process noise and no
measurement noise)
Note:
- is the Moore-Penrose pseudo-inverse of
Z k ( ) Y k ( ) U k ( ) , [ ]
T
= N
Z
k ( ) N
Y
k ( ) N
U
k ( ) [ ]
T
=
N
Z
k ( ) k
f
Z
Z Z
p
, ( ) f
Z k ( )
Z k ( ) Z
p
k ( ) , ( )
k 1 =
F
=
1
2F
det C
Z
k ( ) ( )
k 1 =
F
-------------------------------------------------exp Z k ( ) Z
p
k ( ) ( )
H
C
Z
+
k ( ) Z k ( ) Z
p
k ( ) ( )
k 1 =
F
( ) =
Z
p
k ( )
Y
p
k ( )
U
p
k ( )
= C
Z
k ( ) Cov N
Z
k ( ) ( )

Y
2
k ( )
YU
2
k ( )
YU
2
k ( )
U
2
k ( )
= =
C
Z
k ( )
C
Z
+
k ( ) C
Z
k ( )
150/243
6.b Maximum likelihood estimator (contd)
Periodic signals
Gaussian ML cost function
where is parametrized according to one of the plant models in chapter 4
Eliminate ,
where
Arbitrary signals
Same cost function with
V
ML
Z
p
Z , , ( ) Z k ( ) Z
p
k ( ) ( )
H
C
Z
+
Z k ( ) Z
p
k ( ) ( )
k 1 =
F

k
1 G
k
, ( ) , [ ]Z
p
k ( )
k 1 =
F
+ =
G , ( )
Z
p
Z
p
1 ( ) Z
p
2 ( ) Z
p
F ( ) , , , [ ]
T
=
V
ML
Z , ( )
e
k
Z k ( ) , , ( )
2
e
2

k
, ( )
-------------------------------------
k 1 =
F
=
e
k
Z k ( ) , , ( ) Y k ( ) G
k
, ( )U k ( ) =
e
2

k
, ( )
Y
2
k ( ) G
k
, ( )
2
U
2
2ReG
k
, ( )
YU
2
k ( ) ( ) + =
e
k
Z k ( ) , , ( ) Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) =
151/243
Discussion
with
- suppresses frequency bands with poor signal-to-noise ratio
- relative importance input noise w.r.t. output noise (uncorrelated case )
- significance correlation between input/output errors
with
V
ML
Z , ( )
e
k
Z k ( ) , , ( )
2
e
2

k
, ( )
-------------------------------------
k 1 =
F
=
e
k
Z k ( ) , , ( ) Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) =
e
2

k
, ( )
Y
2
k ( ) G
k
, ( )
2
U
2
k ( ) 2ReG
k
, ( )
YU
2
k ( ) ( ) + =
YU
2
0 =
G
k
, ( )
2
U
2
k ( )
Y
2
k ( )
----------------------------------------
k ( )
2Re
YU
2
k ( )G
k
, ( ) ( )
Y
2
k ( ) G
k
, ( )
2
U
2
k ( ) +
-----------------------------------------------------------
=
k ( ) 0 < decrease
e
2

k
, ( )
k ( ) 0 > increase
e
2

k
, ( )
152/243
Discussion (contd)
- periodic signals + open loop + generator noise dominant in a certain frequency band
high weighting in cost function
generator noise does not increase the variability of the estimates
- arbitrary signals: generator noise is a part of the excitation
- neglect transients
estimate depends on power spectrum input only
same results for multisine or pulse excitation
note: in practise the peak value of the excitation is limited (nonlinearities!)
e
2

k
, ( ) G
0

k
( ) G
k
, ( )
2
g
2
V
ML
Z , ( )
e
k
Z k ( ) , , ( )
2
e
2

k
, ( )
-------------------------------------
k 1 =
F
G
k
( ) G
k
, ( )
2
e
2

k
, ( )
----------------------------------------------- U k ( )
2
k 1 =
F
= =
153/243
Properties
- Consistency
expected value cost function is minimal in consistent
- Convergence rate
with
and : dominating stochastic error
: bias error
- Asymptotic normality
with
V
F
( )
1
F
---V
ML
Z , ( )

1
F
---
e
k
Z
0
k ( ) , , ( )
2
e
2

k
, ( )
---------------------------------------
k 1 =
F

1 + = =
V
F
( )
0
=
Z ( )
0

Z ( ) b
Z ( ) + + =
Z ( ) V
F
''
1

0
( )V
F
'
T
0
Z , ( ) =
Z ( ) O
p
F
1 2 /
( ) =
Z ( ) { } 0 =
b
Z ( ) O
P
F
1
( ) =
F
Z ( )
0
( ) AsN 0 Cov F
Z ( ) ( ) , ( )
Cov F
Z ( ) ( ) V
F
''
1

0
( )Q
F

0
( )V
F
''
1

0
( ) =
Q
F

0
( ) F V
F
'
T
0
Z , ( )V
F
'
0
Z , ( ) { } =
154/243
Properties (contd)
- Cramr-Rao lower bound
the ML estimator is inefficient; ineffiency term is, however, small (Pintelon and Mei,
2007)
- Special cases: input noise only, output noise only, totally correlated input/output errors
the ML estimator is asymptotically efficient
- Robustness: mixing non-Gaussian noise consistent, convergence rate, asymptotic
normality remain valid
- Rational transfer function parametrization: estimate independent of the parameter
constraint chosen

- Model errors (unmodelled dynamics, nonlinear distortions): replace everywhere
with independent of the noise level
Cov F
( ) V
F
''
1

0
( )
Q
F

0
( ) V
F
''
0
( ) =
G , ( ) G , ( ) = V
ML
Z , ( ) V
ML
Z , ( ) =
N
Z
N
Z
Z
0
( )
arg min
V
ML
Z , ( ) { } =
Z
0
( )
155/243
6.c Calculation ML estimates
Newton-Gauss
Basic algorithm
where
and
Numerical stable implementation

where
and
Notes:
- continuous-time, scale angular frequencies by median frequency set

- leave all parameters free and impose in each iteration step the parameter
constraint, for example,
- normalize the columns of by their 2-norm
ReJ
i ( )H
J
i ( )
( )
i 1 + ( )
ReJ
i ( )H
i ( )
( ) =
k [ ]
i ( )
e
k

i ( )
Z k ( ) , , ( )
e

k

i ( )
, ( )
-------------------------------------
= J
k r , [ ]
i ( )

k [ ]
r [ ]
i ( )
-----------
=
J
re
i ( )
i 1 + ( )

re
i ( )
=
i 1 + ( )
J
re
i ( )
( )
+
re
i ( )
V
i ( )

i ( )
( )
+
U
i ( )T
re
i ( )
= =
X
re
ReX ( )
ImX ( )
= J
re
i ( )
U
i ( )
i ( )
V
i ( )T
=
med
med
1

2

F
, , , { } = b
r
s
r
b
r
med
r
( ) s
med
( )
r
=
2
2
1 =
J
re
i ( )
156/243
6.c Calculation ML estimates (contd)
Levenberg-Marquardt
Basic algorithm
with
cost function decreases:
cost function increases:
Numerical stable implementation for the constraint
with
J
re
i ( )T
J
re
i ( )

2
I + ( )
i 1 + ( )
J
re
i ( )T
re
i ( )
=
start

max
J
re
( ) 100 =
0.4
10
2
2
1 =
i 1 + ( )
V
i ( )
i ( )
U
i ( )T
re
i ( )
=
J
re
Udiag
1

2

n
1
0 , , , , ( )V
T
=
diag

1
1
2

2
+
------------------

2
2
2

2
+
------------------

n
1
2

2
+
-------------------------- 0 , , , , ( ) =
157/243
6.d Calculation covariance matrix ML estimates
with
Cov
( ) 2ReJ
H
J ( ) ( )
1
J
k r , [ ]

k [ ]
r [ ]
-----------
=
k [ ]
e
k
Z k ( ) , , ( )
e

k
, ( )
--------------------------------
=
158/243
6.e Starting values
Basic idea for any parametrization
Two possibilities
- rational transfer function + linear least squares + back transformation
- state space + subspace algorithm + back transformation
Linear least squares
Properties
- inconsistent
is not minimal in
- overemphasizes high frequency errors
- estimate is constraint dependent since
- a well chosen weighting can be added WLS
- depends on the noise level
V
LS
Z , ( ) e
k
Z k ( ) , , ( )
2
k 1 =
F
A
k
, ( )Y k ( ) B
k
, ( )U k ( ) I
G

k
, ( )
2
k 1 =
F
= =
V
LS
Z , ( ) F { }
1
F
--- e
k
Z
0
k ( ) , , ( )
2
k 1 =
F
1
F
---
e
2

k
, ( )
k 1 =
F
+ =
0
V
LS
Z , ( )
2
V
LS
Z , ( ) =
Z
0
( ) arg min V
LS
Z , ( ) { } =
159/243
6.e Starting values (contd)
Iterative quadratic maximum likelihood (IQML)
Iterative algorithm
started from (W)LS solution
Other methods
- generalized total least squares (consistent + amplification high frequency errors)
- bootstrapped total least squares (consistent + close to ML + not self starting)
- subspace (consistent, frequency weighting not clear)
V
IQML

i 1 + ( )
Z , ( )
e
k

i 1 + ( )
Z k ( ) , , ( )
2
e
2

k

i ( )
, ( )
------------------------------------------------
k 1 =
F
=
160/243
6.f Simulation example: 2nd order continuous-time system
Simulation parameters
-
-
- in the band
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
-10
-5
0
5
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)
f (Hz)
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
-200
-150
-100
-50
0
50
0 0.1 0.2 0.3
P
h
a
s
e

(
)
f (Hz)
true FRF
measured
FRF
G s
0
, ( ) 1 1 s s
2
+ + ( ) =
U
0
k ( ) 1 =
N
U
k ( ) N
Y
k ( ) , N
c
0 0.02 , ( )
F 100 = 0.1 10 , [ ] 2 ( )
2
Hz
161/243
6.f Simulation example: 2nd order continuous-time system (contd)
Observe the influence of the constraint on the LS estimate:
- underbiased
- overbiased
a
0
= 1
-60
-40
-20
0
20
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)
True plant
LS
LS
f (Hz)
-1
-0.5
0
0.5
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)

e
r
r
o
r
f (Hz)
(b
0
= 1)
(a
0
= 1)
Linear least squares
a
0
1 =
b
0
1 =
162/243
NLS-I/O:
NLS-FRF:
Nonlinear least squares
f (Hz)
-60
-40
-20
0
20
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)
True plant
NLS-I/O
NLS-FRF
f (Hz)
-1
-0.5
0
0.5
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)

e
r
r
o
r
Y k ( ) G
k
, ( )U k ( )
2
k 1 =
F
Y k ( )
U k ( )
---------- -
G
k
, ( )
2
k 1 =
F
163/243
Notes:
- LS uses no noise information
- IQML, and ML use noise information
-5
0
5
0 0.1 0.2 0.3
P
h
a
s
e

e
r
r
o
r

(
)
-1
-0.5
0
0.5
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)

e
r
r
o
r
LS
IQML
ML
164/243
6.g Measurement example: q-axis impedance 3.4 MW synchronous machine
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[[[
[[
[[[
[[[
[[[[[
[[
[[[[[
[[[[
[[[
[[[
[[[[
[[[[[
[[[[
[[[
[[
[[[[
[[
[[[[[[[[[
[[[
[
[[[
[[
[[[
[
[[
-40
-30
-20
-10
0
0.01 0.1 1 10 100
A
L
S

(
d
B
)
[
[
[
[
[
[
[[[[
[[
[
[[
[
[[
[[[[[
[
[[[[[
[
[
[
[[[[[
[[
[[
[[[
[[
[[
[
[
[[[
[
[
[
[[
[[
[
[[
[
[
[
[
[
[[[[
[
[
[[
[[
[[
[
[
[
[
[[[
[
[
[
[
[
[[[
[
[
[[[
0
20
40
60
80
0.01 0.1 1 10 100
F
L
S

(
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[[[
[[
[[[
[[[
[[[[[
[[
[[[[[
[[[[
[[[
[[[
[[[[
[[[[[
[[[[
[[[
[[
[[[[
[[
[[[[[[[[[
[[[
[
[[[
[[
[[[
[
[[
-40
-30
-20
-10
0
0.01 0.1 1 10 100
A
I
Q
M
L

(
d
B
)
[
[
[
[
[
[
[[[[
[[
[
[[
[
[[
[[[[[
[
[[[[[
[
[
[
[[[[[
[[
[[
[[[
[[
[[
[
[
[[[
[
[
[
[[
[[
[
[[
[
[
[
[
[
[[[[
[
[
[[
[[
[[
[
[
[
[
[[[
[
[
[
[
[
[[[
[
[
[[[
0
20
40
60
80
0.01 0.1 1 10 100
F
I
Q
M
L

(
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[[[
[[
[[[
[[[
[[[[[
[[
[[[[[
[[[[
[[[
[[[
[[[[
[[[[[
[[[[
[[[
[[
[[[[
[[
[[[[[[[[[
[[[
[
[[[
[[
[[[
[
[[
-40
-30
-20
-10
0
0.01 0.1 1 10 100
A
M
L

(
d
B
)
[
[
[
[
[
[
[[[[
[[
[
[[
[
[[
[[[[[
[
[[[[[
[
[
[
[[[[[
[[
[[
[[[
[[
[[
[
[
[[[
[
[
[
[[
[[
[
[[
[
[
[
[
[
[[[[
[
[
[[
[[
[[
[
[
[
[
[[[
[
[
[
[
[
[[[
[
[
[[[
0
20
40
60
80
0.01 0.1 1 10 100
F
M
L

(
2
2
1 =
b
0
1 =
n
b
4 n
a
, 3 = =
2
2
1 =
b
0
1 =
estimated
model
measured
FRF
165/243
6.h Measurement example: flight flutter data
n
b
11 n
a
, 10 = =
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[
[[
[
[
[
[
[
[[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[
[
[
[[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[[[
[
[
[[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[[[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[[
[
[
[
[
[[
[[
[
[
[
[
[[
[
[
[
[
[
[[
[
[[
[
[
[
[[
[[
[
[
[
[
-20
-15
-10
-5
0
4 5 6 7 8 9 10 11
f (Hz)
A
L
S
(
d
B
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[[
[[
[
[
[
[[[
[
[[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[
[
[[[
[
[
[
[[[
[
[[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
0
50
100
150
4 5 6 7 8 9 10 11
f (Hz)
F
L
S
(
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[
[[
[
[
[
[
[
[[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[
[
[
[
[[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[[[
[
[
[[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[[[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[[
[
[
[
[
[[
[[
[
[
[
[
[[
[
[
[
[
[
[[
[
[[
[
[
[
[[
[[
[
[
[
[
[
-20
-15
-10
-5
0
4 5 6 7 8 9 10 11
f (Hz)
A
I
Q
M
L
(
d
B
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[[
[[
[
[
[
[[[
[
[[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[
[
[[[
[
[
[
[[[
[
[[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
0
50
100
150
4 5 6 7 8 9 10 11
f (Hz)
F
I
Q
M
L
(
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[
[[
[
[
[
[
[
[[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[
[
[
[
[[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[[[
[
[
[[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[[[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[[
[
[
[
[
[[
[[
[
[
[
[
[[
[
[
[
[
[
[[
[
[[
[
[
[
[[
[[
[
[
[
[
[
-20
-15
-10
-5
0
4 5 6 7 8 9 10 11
f (Hz)
A
M
L
(
d
B
)
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[[
[[
[
[
[
[[[
[
[[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[
[
[[[
[
[
[
[[[
[
[[
[
[
[
[
[
[
[
[[
[[
[
[
[
[
[
[[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[[
[
[
[
[
[[
[
[
[[
[
[
[
[
[
[[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
[
0
50
100
150
4 5 6 7 8 9 10 11
f (Hz)
F
M
L
(
)
estimated
measured
FRF
model
166/243
6.i Model validation
Comparison with measured FRF
Compare to the nonparametric LPM estimate (see Chapter 5)
taking into account the uncertainty of the LPM estimate
Analysis cost function
No modelling errors
with
Test
modelling errors
no modelling errors
wrong noise variances?
G
k

, ( )
G
fast

k
( ) G
robust

k
( ) , for periodic signals
G
arb

k
( ) for arbitrary signals
2
k ( )
V
ML

ML
Z ( ) Z , ( ) { } V
noise
var V
ML

ML
Z ( ) Z , ( ) ( ) V
noise
V
noise
F n
2 =
n
number of free parameters =

V
ML

ML
Z ( ) Z , ( ) V
noise
2 V
noise
+ >
V
ML

ML
Z ( ) Z , ( ) V
noise
2 V
noise
V
noise
2 V
noise
+ , [ ]
V
ML

ML
Z ( ) Z , ( ) V
noise
2 V
noise
<
167/243
6.i Model validation (contd)
Whiteness test residuals measured minus modelled FRF
Sample correlation residuals
with
No unmodelled dynamics
- 95% confidence bound
- is asymptotically ( ) a Dirac function
- and noise (co-)variances no nonlinear distortions
- and noise (co-)variances nonlinear distortions
- (co-)variances noise + stochastic nonlinear distortions
R
m ( )
1
F m
------------- -
k

, ( )
k m +

, ( )
k 1 =
F m
=
k

, ( )
G

k
( ) G
k

, ( )

k
( )
------------------------------------------
=
3std R
m ( ) ( ) 3 F m ( ) =
R
m ( ) F
R
m ( ) O F
1 2 /
( ) for m 0 =
R
0 ( ) 1
0 ( ) 1 =
R
0 ( ) 1 >
R
0 ( ) 1 =
168/243
6.j Measurement example: flexible robot arm
Observation: and nonlinear distortions
-60
-40
-20
0
20
0 5 10 15 20 25
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
Model 4/4 Model 6/6
0
20
40
60
-100 -50 0 50 100
A
m
p
l
i
t
u
d
e
Lag number
0
1
2
3
-100 -50 0 50 100
A
m
p
l
i
t
u
d
e
Lag number
-60
-40
-20
0
20
0 5 10 15 20 25
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
(a)
(b)
V
ML
4964.8 =
V
noise
95.5 =
V
ML
220.5 =
V
noise
93.5 =
F 100 =
meas. FRF
estim. model
x
std FRF
residual
meas. - model
sample
correlation
95% conf. bound
50% conf. bound
R
0 ( ) 1 >
G
2
var N
G
( ) =
169/243
6.k Multivariable systems
ML cost function
with
and
Difficulty: calculation Jacobian matrix
Solution: pseudo-Jacobian (see Guillaume and Pintelon, 1996)
V
ML
Z , ( ) e
H
k
Z k ( ) , , ( )C
e
1

k
, ( )e
k
Z k ( ) , , ( )
k 1 =
F
=
e
k
Z k ( ) , , ( ) Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) =
C
e

k
, ( ) C
Y
k ( ) G
k
, ( )C
U
k ( )G
H
k
, ( ) 2hermC
YU
k ( )G
H
k
, ( ) ( ) + =
hermX ( ) X X
H
+ ( ) 2 =
J
k r , [ ]
C
e
1 2 /

k
, ( )e
k
Z k ( ) , , ( )
r [ ]
-----------------------------------------------------------------
=
J
+k r , [ ]
C
e
1 2 /

k
, ( )
e
k
Z k ( ) , , ( )
r [ ]
----------------------------------- + =
1
2
-- -C
e
1 2 /

k
, ( )
C
e

k
, ( )
r [ ]
-------------------------C
e
1

k
, ( )e
k
Z k ( ) , , ( )
170/243
6.l Model selection
Problem statement
true 5th order FIR
standard deviation
estimated 5th order
FIR model
-40
-30
-20
-10
0
0 0.1 0.2 0.3
A
m
p
l
i
t
u
d
e

(
d
B
)
Frequency (Hz)
standard deviation
estimated 50th order
FIR model
system
171/243
6.l Model selection (contd)
As model complexity increases
- cost function decreases AND model variability increases
Optimal model complexity makes a balance between
- model errors AND model variability
Model selection criteria
with ML cost function
and penalty
2F
2F n
------------------

V
F

Z ( ) Z , ( ) 1 p n
F , ( ) + ( )
V
F
Z , ( )
1
F
--- e
H
k
Z k ( ) , , ( )C
e
1

k
, ( )e
k
Z k ( ) , , ( )
k 1 =
F
=
p n
F , ( )
2n
2F n
------------------ AIC
2 n
y
n
u
+ ( )F ( )n
log
2F n
----------------------------------------------- MDL (EIV)

2n
y
F ( )n
log
2F n
------------------------------- MDL (OE)
=
172/243
6.m Simulation example: SISO-FIR system
with
Input/output DFT spectra
, , , ,
Results 1000 runs
model order AIC MDL
17
0 0
18
0 0
19
845 923
20
78 49
21
45 18
22
16 4
23
10 4
24
6 2
mean value std
19.3 0.79 19.1 0.51
G z
1
, ( ) g
r
z
r
r 0 =
R
= g
0
g
1
g
R
, , , [ ]
T
=
U
0
k ( ) 1 = Y
0
k ( ) G
0
z
k
1
( ) =
Y
k ( )
U
k ( ) 0.05 = =
YU
k ( ) 0 = F 40 =
R
173/243
6.n References
De Ridder F., R. Pintelon, J. Schoukens and D.P. Gillikin (2005). Modified AIC and MDL
model selection criteria for short data records, IEEE Trans. Instrum. Meas., vol. 53, no. 5,
pp. 144-150.
Guillaume P. and R. Pintelon (1996). A Gauss-Newton-Like Optimization Algorithm for
Weighted Nonlinear Least-Squares Problems, IEEE Trans. Sign. Proc., vol. SP-44, no. 9,
pp. 2222-2228.
Levy, E. C. (1959). Complex curve fitting. IEEE Trans. Autom.Contr., vol. AC-4, pp. 37-43.
McKelvey, T., H. Akay and L. Ljung (1996). Subspace-based multivariable system
identification from frequency response data, IEEE Trans. Autom. Contr., vol. 41, no 7, pp.
960-979.
Pintelon, R. and M. Hong (2007). Asymptotic uncertainty of transfer function estimates
using nonparametric noise models, IEEE Trans. Instrum. Meas., vol. 56, no. 6, in press.
Pintelon R., and I. Kollr (2005). On the frequency scaling in continuous-time modeling,
IEEE Trans. Instrum. Meas., vol. 53, no. 5, pp. 318-321.
Sanathanan, C. K. and J. Koerner (1963). Transfer function synthesis as a ratio of two
complex polynomials. IEEE Trans. Autom. Contr., vol. AC-8, pp. 56-58.
Schoukens J., Y. Rolain and R. Pintelon (2002). Modified AIC rule for model selection in
combination with prior estimated noise models, Automatica, vol. 38, no. 5, pp. 903-906.
Van Overschee, P. and B. De Moor (1996). Continuous-time frequency domain subspace
system identification, Signal Processing, vol. 52, no 2, pp. 179-194.
174/243
7.a Classical approach: basic idea 175
7.b Classical approach: sample maximum likelihood 176
7.c Measurement example: identification of Youngs modulus 180
7.d Classical approach: problems 183
7.e Local polynomial method: basic idea 184
7.f Local polynomial method: MIMO sample maximum likelihood 185
7.g MIMO experiment: modal analysis Aluminium tooling plate 187
7.h Measurement example: Magnetic Resonance Imaging (MRI) 194
7.i References 199
175/243
7.a Classical approach: basic idea
Extract noise covariance matrix from the experimental data
Ideal situation
- perform independent repeated experiments
, , and
where , are independent (over and ) jointly correlated
circular complex normally distributed disturbances
- calculate sample means and sample (co-)variances
Practise
- consecutive periods of the steady state response to a multisine excitation
- input/output disturbances = filtered white noise
- independence (over and ) and normality are asymptotically valid for
with fixed
M
Y
m [ ]
k ( ) Y
0
k ( ) N
Y
m [ ]
k ( ) + =
U
m [ ]
k ( ) U
0
k ( ) N
U
m [ ]
k ( ) + =
m 1 2 M , , , = k 1 2 F , , , =
N
U
m [ ]
k ( ) N
Y
m [ ]
k ( ) k m
2
k ( )
2
k ( )
2
k ( ) , ,
M
k m
N M
176/243
7.b Classical approach: sample maximum likelihood
Use in ML cost function sample means and sample (co-)variances sample ML
with
Questions
- minimal value ?
- influence on the asymptotic properties of ?
Key property to analyze the asymptotic behaviour of
- for normally distributed errors the sample mean and sample covariance
matrix are independently distributed
V
SML
Z , ( )
e
k
Z
k ( ) , , ( )
2
e
2
k
, ( )
-------------------------------------
k 1 =
F
=
e
k
Z k ( ) , , ( ) Y
k ( ) G
k
, ( )U
k ( ) =
e
2
k
, ( )
2
k ( ) G
k
, ( )
2
2
k ( ) 2ReG
k
, ( )
2
k ( ) ( ) + =
SML
SML
177/243
7.b Classical approach: sample maximum likelihood (contd)
Asymptotic properties
Consistency ( )
minimal in consistent
Convergence rate, systematic and stochastic errors, asymptotic normality ( )
same as ML (see section 6.b on page 153) with
small loss in efficiency (below 12% for )
M 4
V
SML
( ) V
SML
Z , ( ) { }
e
k
Z
k ( ) , , ( )
2
e
2
k
, ( )
-------------------------------------

k 1 =
F
= =
e
k
Z
k ( ) , , ( )
2
{ }
1
e
2
k
, ( )
----------------------

k 1 =
F
=
M 1
M 2
--------------V
ML
( ) =
0
M 7
Cov
SML
Z ( ) ( )
M 2
M 3
--------------Cov
ML
Z ( ) ( )
M 7
178/243
Asymptotic properties (contd)
Calculation uncertainty
with
Analysis cost function
No modeling errors
Cov
SML
Z ( ) ( )
M 1
M 3
------------- - 2ReJ
SML
H
J
SML
( ) [ ]
1
J
SML
( )
k r , [ ]

k [ ]
SML r [ ]
-------------------
=
k [ ]
e
k
Z
k ( ) , , ( )
e

k
, ( )
--------------------------------
=
V
SML

SML
Z ( ) Z , ( ) { }
M 1
M 2
------------- -V
noise
var V
SML

SML
Z ( ) Z , ( ) ( )
M 1 ( )
3
M 2 ( )
2
M 3 ( )
----------------------------------------V
noise
with
V
noise
F n
2 =
n
number of free parameters =
179/243
Asymptotic properties (contd)
Sample correlation residuals (Pintelon et al., 2003)
with
Notes:
- 95% confidence bound with
- has exactly the same properties as in section 6.i on page 167
- compensates for replacing the true variance by the sample variance
R
m ( )

1
m ( )
F m
--------------
, ( )
k m +

, ( )
k 1 =
F m
, ( )
G

k
( ) G
k

, ( )
G

k
( )
------------------------------------------
=
1
m ( )
M 2 ( ) M 1 ( ) m 0 =
M 5 3 ( ) M 11 12 ( ) m 0
2
m ( ) 3 F m ( )

2
m ( )

1
m ( ) 3 M 1 ( )
3 2 /
M 2 ( ) M 3 ( )
1 2 /
( ) m 0 =
1
m ( ) M 1 ( ) M 2 ( ) m 0
=
R
m ( ) R
m ( )
i
m ( )
180/243
7.c Measurement example: identification of Youngs modulus
u t ( )
0 x L 1.983 m =
y t x , ( )
-120
-80
-40
0
40
0 1000 2000 3000 4000
A
m
p
l
i
t
u
d
e

(
d
B
)
frequency (Hz)
measured FRF
total variance
noise variance
Force-to-acceleration
M 25 P 10 = , =
plexiglass beam
1 cm 2 cm 1.983 m
181/243
7.c Measurement example: identification of Youngs modulus (contd)
-120
-80
-40
0
40
0 1000 2000 3000 4000
A
m
p
l
i
t
u
d
e

(
d
B
)
frequency (Hz)
Force-to-acceleration
measured FRF
total variance
difference modelled
and measured FRF
n
a
n
b
34 = = CT-model
V
SML

Z , ( ) 13023 =
V
noise
1602 =
182/243
7.c Measurement example: identification of Youngs modulus (contd)
E
E
E
E
E
E
E
1.5
2
2.5
3
0 1 2 3 4
P
h
a
s
e

(
)
|poles|/(2 ) (kHz)
E
E E
E
E
E
E
194.2
194.4
194.6
194.8
0 1 2 3 4
A
m
p
l
i
t
u
d
e

(
d
B
N
/
m
2
)
|poles|/(2 ) (kHz)
E s
k
( )
Ls
k
k
--------

2
1
k s
k
( )J
L
-----------------------

2
+ =
for
k 1 2 7 , , , =
identified model of order
o measured Youngs modulus
n
a
n
b
2 = =
V
ML

Z , ( ) 6.24
4
10 =
V
noise
4.5 =
s
k
( ) 0.3 ( )
183/243
7.d Classical approach: problems
Drawbacks
periodic steady state only
consecutive periods are correlated due to the noise transients sample means and
sample (co-)variances are not independently distributed
sensitive to noise transient (leakage) errors
sensitive to system transient (leakage) errors
Solution
the local polynomial method (reason: transient suppression)
184/243
7.e Local polynomial method: basic idea
Arbitrary excitations within an generalised output error framework
Step 1: local polynomial estimates (see p. 114)
- FRF:
- noise covariance: with degrees of freedom
Step 2: generalised sample mean and sample covariance (leakage free!)
- sample mean:
- sample covariance sample mean:
where
with defined on p. 116
Note: errors-in-variables and periodic excitations: follow the same lines
G

k
( )
C
V
k ( ) dof
Y
k ( ) G

k
( )U k ( ) =
C
k ( ) q
n
2
2
C
V
k ( ) =
q
n
K
n
H
K
n
K
n
H
( )
1
U k ( )
0
=
K
n
185/243
7.f Local polynomial method: MIMO sample maximum likelihood
SML cost function
with (no transient term!)
Question
- minimal value ?
- influence on the asymptotic properties of ?
Key property to analyze the asymptotic behaviour of
- the generalised sample mean and sample covariance are asymptotically
independently distributed
V
SML
Z , ( ) e
H
k
Z
k ( ) , , ( )C
e
1
k ( )e
k
Z
k ( ) , , ( )
k 1 =
F
=
e
k
Z
k ( ) , , ( ) Y
k ( ) G
k
, ( )U
k ( ) =
C
e

k
, ( ) C
k ( ) G
k
, ( )C
k ( )G
H
k
, ( ) 2hermC
k ( )G
H
k
, ( ) ( ) + =
dof
SML
SML
186/243
7.f Local polynomial method: MIMO sample maximum likelihood (contd)
Asymptotic properties
( )
( )
( )
( )
Notes:
- pseudo jacobian (see p. 169)
- classical approach:
- whiteness test residuals (see p. 167), where , is valid over a
subset of frequencies (see p. 118, correlation length LPM estimates)
Cov
SML
( )
dof dof n
y
( )
dof n
y
1 + ( ) dof n
y
1 ( )
-------------------------------------------------------------------Cov
ML
( ) = dof n
y
8 +
V
SML

SML
Z , ( ) { }
dof
dof n
y
-------------------V
noise
V
noise
n
y
F n
2 = dof n
y
2 +
1
2
var V
SML

SML
Z , ( ) ( ) 3
1
2

1
2
dof
3
dof n
y
( )
2
dof n
y
1 ( )
------------------------------------------------------------V
noise
= dof n
y
4 +
Cov
SML
( )
dof
2
dof n
y
1 + ( ) dof n
y
1 ( )
------------------------------------------------------------------- 2ReJ
SML+
H
J
SML+
( ) ( )
1
dof n
y
8 +
J
SML+
M dof 1 + =
M dof 1 + =
187/243
7.g MIMO experiment: modal analysis Aluminium tooling plate
size Al plate: 30.4 cm x 61.8 cm x 6.7 mm
188/243
7.g MIMO experiment: modal analysis Aluminium tooling plate (contd)
Measurement
two consecutive periods of the transient response to 1 set of uncorrelated
random phase multisines
points per period
Local polynomial method for periodic excitations
fast method with and
correlation length nonparametric estimates is
Parametric model
continuous-time common denominator model
model order selection via MDL and AIC (see p. 171)
n
u
2 =
f
s
2.44 kHz, = N 64 1024 =
R 4 = dof 9 = n 9 =
18
189/243
MDL criterion:
Parsimonious principle:
MDL AIC
6/4 9855 462.5 [29.9, 51.7] 12920 10690
6/5 9828 462.2 [29.9, 51.7] 12980 10690
6/6 527.4 461.6 [29.8, 51.7] 701.9 575.1
7/6 394.4 459.0 [29.7, 51.5] 540.3 434.4
8/6 365.2 456.4 [29.7, 51.4] 514.5 406.0
8/7 363.0 455.8 [29.6, 51.4] 515.0 404.6
n
b
n
a
V
SML V
SML
{ } std V
SML
( )
n
b
n
a
8 6 =
n
b
n
a
7 6 =
190/243
248 250 252 254
100
50
0
Frequency (Hz)
G
[
1
,
1
]

(
d
B
)
248 250 252 254
100
50
0
Frequency (Hz)
G
[
1
,
2
]

(
d
B
)
248 250 252 254
100
50
0
Frequency (Hz)
G
[
2
,
1
]

(
d
B
)
248 250 252 254
100
50
0
Frequency (Hz)
G
[
2
,
2
]

(
d
B
)
G

k
( ) G s
k

SML
, ( ) ,
var G

k
( ) ( ) G

k
( ) G s
k

SML
, ( )
n
b
n
a
8 6 =
191/243
248 250 252 254
200
150
100
50
0
Frequency (Hz)
a
r
g
(
G
[
1
,
1
]
)

(
)
248 250 252 254
400
300
200
100
0
Frequency (Hz)
a
n
g
l
e
(
G
[
1
,
2
]
)

(
)
248 250 252 254
400
300
200
100
0
Frequency (Hz)
a
r
g
(
G
[
2
,
1
]
)

(
)
248 250 252 254
200
150
100
50
0
Frequency (Hz)
a
r
g
(
G
[
2
,
2
]
)

(
)
G s
k

SML
, ( )
G

k
( )
192/243
20 10 0 10 20
0.0
0.5
1.0
1.5
k
R
[
1
,
1
]
(
k
)
20 10 0 10 20
0.0
0.5
1.0
1.5
k
R
[
1
,
2
]
(
k
)
20 10 0 10 20
0.0
0.5
1.0
1.5
k
R
[
2
,
1
]
(
k
)
20 10 0 10 20
0.0
0.5
1.0
1.5
k
R
[
2
,
2
]
(
k
)
*
R
m ( )
50% confidence bound
95% confidence bound
(42% * outside)
(2% * outside)
193/243
Modal parameters for the MDL (8/6) and parsimonious (7/6) models
Difference
Estimates
(Hz)

(Hz)

n
b
n
a
8 6 = n
b
n
a
7 6 =
f
0
std f
0
( )
250.8552 2.0
5
10 250.8550 2.0
5
10 2.0
4
10
std ( )
2.96
4
10 3.2
6
10 2.97
4
10 3.1
6
10 1.4
6
10
f
0
std f
0
( )
251.7844 3.8
5
10 251.7840 3.7
5
10 4.0
4
10
std ( )
9.0
5
10 6.0
6
10 9.2
5
10 5.9
6
10 1.5
6
10
194/243
7.h Measurement example: Magnetic Resonance Imaging (MRI)
MRI scanner:
~ Tesla static magnetic field,
~ MHz oscillating field (pulse) perpendicular to the static field
response measured in two orthogonal directions x and y
complex signal x(t) + jy(t)
195/243
7.h Measurement example: Magnetic Resonance Imaging (MRI) (cont)
0 0.04 0.08 0.12
0
100
200
300
tijd (s)
a
b
s
(
r
e
s
p
o
n
s
i
e
)
absolute value demodulated signal x(t) + jy(t)
(averaged over M = 64 independent measurements)
time (s)
a
b
s
(
r
e
s
p
o
n
s
e
)
196/243
-2000 -1000 0 1000 2000
-40
-20
0
20
40
measured
spectrum
model
residual
meas.-model
noise variance
frequentie (Hz)
A
m
p
l
i
t
u
d
e

(
d
B
)
T z
1
, ( )
b
r
z
r
r 0 =
n 1
a
r
z
r
r 0 =
n
--------------------------
=
a
r
b
r
,
n 9 =
signal model = sum of complex
damped exponentials
NMR spectrum muscle
197/243
-300 -200 -100 0 100 200 300
0
0.4
0.8
1.2
Function of Dependency analysis + : measured - bound 50 % - - bound 95%
Lag (in samples)
A
u
t
o
c
o
r
r
e
l
a
t
i
o
n

r
e
s
i
d
u
a
l
s
autocorrelation
50% uncertainty
bound (fraction
95% uncertainty
bound (fraction
Whiteness test residuals
V
SML

Z , ( ) 584 =
V
noise
502 =
V
noise
1 2 /
22 =
outside = 51.6%)
outside = 5.2%)
198/243
NMR image brains
frequency and damping peak
depend on metabolite
amplitude peak is a measure of
the concentration of the metabolite
concentration of a metabolite (e.g.
creatine) in a color scale as a
function of space NMR image
199/243
7.i References
Mahata K., R. Pintelon and J. Schoukens (2006). On parameter estimation using non-
parametric noise models, IEEE Trans. Autom. Contr., vol. 51, no. 10, pp. 1602-1612 .
Pintelon R. (2002). Frequency domain subspace identification using nonparametric noise
models, Automatica, vol. 38, no. 8, pp. 1295-1311.
Pintelon R., J. Schoukens and Y. Rolain (2003). Uncertainty of transfer function modeling
using prior estimated noise models, Automatica, vol. 39, no. 10, pp. 1721-1733.
Schoukens J., G. Vandersteen, R. Pintelon and P. Guillaume (1997). Frequency Domain
System Identification Using Non-parametric Noise Models Estimated from a Small Number
of Data Sets, Automatica, vol. 33, no. 6, pp. 1073-1086.
200/243
8. Estimation parametric plant and/or noise models
8.a Classical time domain prediction error framework 201
8.b Frequency domain identification in closed loop with known controller 209
8.c Frequency domain identification in closed loop with unknown controller 213
8.d Measurement example: nonlinear electrical circuit 217
8.e Measurement example: Villa Passo bridge 219
8.f Measurement example: flight flutter 221
8.g Multivariable systems in closed loop with known controller 223
8.h Measurement example: Villa Passo bridge - MIMO case 228
8.i Measurement example: flight flutter - MIMO case (contd) 230
8.j References 231
201/243
8.a Classical time domain prediction error framework
Discrete-time plant and noise model
Model
with a monic noise model
Notes:
- and are known exactly and is not observed
- is not observed is monic by appropriate scaling
- and should be stable
v t ( )
y t ( )
u t ( )
G q ( )
H q ( )
e t ( )
y t ( ) G q ( )u t ( ) H q ( )e t ( ) + =
H z
1
( )
1 c
1
z
1
+ +
1 d
1
z
1
+ +
----------------------------------
=
u t ( ) y t ( ) e t ( )
e t ( ) H z
1
( )
H z
1
( ) H
1
z
1
( )
202/243
8.a Classical time domain prediction error framework (contd)
One-step-ahead prediction
One-step-ahead predictor process noise
= best (in least square sense) estimate of based on the knowledge of all past
samples
One-step-ahead predictor output
Prediction error
v t ( ) H q ( )e t ( ) e t ( ) H q ( ) 1 ( )e t ( ) + e t ( ) 1 H
1
q ( ) ( )v t ( ) + = = =
v
t t 1 ( ) 1 H
1
q ( ) ( )v t ( ) =
v t ( )
v t 1 ( ) v t 2 ( ) v ( ) , , ,
y
t t 1 ( ) G q ( )u t ( ) v
t t 1 ( ) + H
1
q ( )G q ( )u t ( ) 1 H
1
q ( ) ( )y t ( ) + = =
t ( ) y t ( ) y
t t 1 ( ) H
1
q ( ) y t ( ) G q ( )u t ( ) ( ) = =
203/243
Prediction error estimator
Prediction error cost function
Notes:
- the effect of the initial conditions is hidden in the filter operations
- the predictor should be stable, viz., and
- according to the particular model structure one distinguishes
ARX: , and linear least squares
ARMAX:
ARMA: time series/spectral analysis
OE:
BJ: and are independently parametrized
- frequency domain equivalent (Parcevals theorem)
V
PE
z , ( )
2
t , ( )
t 0 =
N 1
H
1
q , ( ) y t ( ) G q , ( )u t ( ) ( ) ( )
2
t 0 =
N 1
= =
H
1
H
1
G
C 1 = D A =
D A =
G 0 =
H 1 =
G H
2
t , ( )
t 0 =
N 1
z
k
1
, ( )
2
k 0 =
N 1
=
z
k
1
, ( ) H
1
z
k
1
, ( ) Y k ( ) G z
k
1
, ( )U k ( ) T
G
z
k
, ( ) T
H
z
k
, ( ) ( ) =
204/243
Properties estimated plant model parameters
Consistent
Open loop
OE, BJ: even with wrong noise model structure
ARX, ARMAX: correct noise model structure
Closed loop
OE, BJ, ARX, ARMAX: correct noise model structure
Asymptotically normally distributed
Maximum likelihood for normally distributed noise
Asymptotically efficient for normally distributed noise
Cramr-Rao lower bound: see Ljung (1999)
205/243
Focusing on a particular frequency band - goal
Observation
- the frequency content of a sampled signal extends from DC to Nyquist
Goal
- removal slow trends and/or high frequency disturbances
- removal (harmonic) disturbances that cannot be written as filtered
white noise
- simplified plant and/or noise model
206/243
Focusing on a particular frequency band - time domain
G z
1
( )
u
f
t ( )
F z
1
( ) F z
1
( )
H z
1
( )
u t ( ) y t ( ) G q ( )u t ( ) H q ( )e t ( ) + =
e t ( )
y
f
t ( ) G q ( )u
f
t ( ) F q ( )H q ( )e t ( ) + =
f
t ( ) F q ( )H q ( ) ( )
1
y
f
t ( ) G q ( )u
f
t ( ) ( ) H
1
q ( ) y t ( ) G q ( )u t ( ) ( ) t ( ) = = =
207/243
Focusing on a particular frequency band - time domain
Properties
- input/output relation is preserved
- prefiltering is equivalent to dividing the noise model by the prefilter
Requirement
- noise model should be flexible enough to preserve efficiency (open/closed
loop) and consistency (closed loop)
Conclusion
- classical time domain solution = compromise between suppression
undesired frequency bands and loss in efficiency and consistency
208/243
Focusing on a particular frequency band - frequency domain
Properties
- exact non-causal filtering
- no loss in consistency
- increased variance
0
j
1
s-domain
z-domain
209/243
8.b Frequency domain identification in closed loop with known controller
Closed loop equations
Conditional expectations ( )
U k ( )
G
k
( )
M
0

k
( )
R k ( ) Y k ( )
V k ( ) H
k
( )E k ( ) =

k
, ( )
Y k ( ) G
k
, ( )U k ( )
H
k
, ( )
------------------------------------------------
= prediction error
U k ( )
1
1 G
k
( )M
0

k
( ) +
------------------------------------------R k ( )
M
0

k
( )H
k
( )
1 G
k
( )M
0

k
( ) +
------------------------------------------E k ( ) =
Y k ( )
G
k
( )
1 G
k
( )M
0

k
( ) +
------------------------------------------R k ( )
H
k
( )
1 G
k
( )M
0

k
( ) +
------------------------------------------E k ( ) + =
var E k ( ) ( ) =
Y k ( ) R k ( ) , , { }
G
k
, ( )
1 G
k
, ( )M
0

k
( ) +
------------------------------------------------R k ( ) Y k , ( ) =
var Y k ( ) R k ( ) , , ( ) var
H
k
, ( )
1 G
k
, ( )M
0

k
( ) +
------------------------------------------------E k ( )

S
k
, ( )
2
= =
210/243
8.b Frequency domain identification in closed loop with known controller (contd)
Gaussian ML estimator
Gaussian likelihood function of
Negative Gaussian log-likelihood of all data
Eliminating gives
Notes:
- open loop (Ljung, 1999)
- knowledge controller necessary (McKelvey, 2000) besides store
also
- numerically stable Newton-Gauss scheme of the form
Y k ( )
f
Y k ( )
Y k ( ) R k ( ) , , ( )
1
S
k
, ( )
2
--------------------------------exp
Y k ( ) Y k , ( )
2
S
k
, ( )
2
-------------------------------------

=
log f
Y k ( )
Y k ( ) R k ( ) , , ( ) ( )
k
C
1
S
k
, ( )
2
( ) log
k

k
, ( )
2
------------------------
k
+ + =
V
F
Z , ( )
1
F
---
k
, ( )g
F
( )
2
k
= with g
F
( ) exp
1
F
---
H
k
, ( )
1 G
k
, ( )M
0

k
( ) +
------------------------------------------------

log
k

=
M
0
0 =
u t ( ) y t ( ) ,
r t ( )
J e =
211/243
Improving the finite sample behaviour
with
with prediction error
Model structures
DT-ARX, CT-ARX: , , and nonlinear least squares
DT-ARMAX, CT-ARMAX: , , and
DT-ARMA, CT-ARMA: , and time series/spectral analysis
DT-OE, CT-OE: , and
DT-BJ, CT-BJ: and are independently parametrized
Hybrid-BJ: CT-plant and DT-noise models
V
F
Z , ( )
1
F
---
k
, ( )g
F
( )
2
k
= g
F
( ) exp
1
F
--- S
k
, ( ) ( ) log
k

=

k
, ( )
Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) T
H

k
, ( )
H
k
, ( )
----------------------------------------------------------------------------------------------------------
=
C 1 = D A = T
H
0 =
D A = T
H
0 = n
i
g
max n
a
n
b
n
c
, , ( ) 1
G 0 = T
G
0 =
H 1 = T
H
0 =
G H
212/243
Connection with the classical prediction error framework
with
If
1.
2. covers uniformly the unit circle
then
with the dominant pole of
Notes:
- condition 1 is fulfilled if
monic noise model
at least one sample delay in plant and/or controller
- knowledge controller no longer necessary
- classical prediction error method
V
F
Z , ( )
1
F
--- z
k
1
, ( )g
F
( )
2
k
= g
F
( ) exp
1
F
---
H z
k
1
, ( )
1 G z
k
1
, ( )M
0
z
k
1
( ) +
-------------------------------------------------

log
k

=
S z
1
, ( )
z
lim 1 =
g
F
( ) exp
1
N
--- - S z
k
1
, ( ) ( ) log
k 0 =
N 1

1 O
max
N
N ( ) + = =
max
S z
1
, ( ) ( ) log
g
F
( ) 1 =
213/243
8.c Frequency domain identification in closed loop with unknown controller
Gaussian likelihood function of ( )
U k ( )
G
k
( )
M
k
( )
R k ( ) L
k
( )W k ( ) =
Y k ( )
V k ( ) H
k
( )E k ( ) =
G

k
, ( )
Y k ( ) G
k
, ( )U k ( )
H
k
, ( )
------------------------------------------------
=
M

k
, ( )
U k ( ) M
k
, ( )Y k ( ) +
L
k
, ( )
-------------------------------------------------
=
plant prediction error
controller prediction error
Z k ( ) Y k ( ) U k ( ) , [ ]
T
= var E k ( ) ( ) = , var W k ( ) ( ) =
Z k ( ) , , { } 0 =
C
Z k ( )
Z
H
k ( )Z k ( ) , , { }
1
1 GM +
2
------------------------
GL
2
H
2
+ G L
2
M H
2
G L
2
M H
2
L
2
HM
2
+
= =
f
Z k ( )
Z k ( ) , , ( )
1
2
detC
Z k ( )
------------------------exp Z
H
k ( )C
Z k ( )
1
Z k ( ) ( ) =
214/243
8.c Frequency domain identification in closed loop with unknown controller (contd)
Gaussian ML estimator (contd)
Negative Gaussian log-likelihood function all data
where
Eliminating and gives
Notes:
- symmetric in , and ,
- identification , is coupled with that of , via
- numerically stable Newton-Gauss scheme of the form
(see Pintelon and Schoukens, 2006)
log f
Z k ( )
Z k ( ) , , ( ) ( )
k
C
1
T
k
, ( )
2
( ) log
k

G

k
, ( )
2
---------------------------
k

M

k
, ( )
2
----------------------------
k
+ + + =
T
k
, ( )
H
k
, ( )L
k
, ( )
1 G
k
, ( )M
k
, ( ) +
----------------------------------------------------
=

V
F
Z , ( ) h
F
( )
2
1
F
---
G

k
, ( )
2
k

1
F
---
M

k
, ( )
2
k

=
h
F
( ) exp
1
F
--- T
k
, ( ) ( ) log
k

=
G H M L
G H M L T HL 1 GM + ( ) =
J e =
215/243
with plant and controller prediction errors
V
F
Z , ( ) h
F
( )
2
1
F
---
G

k
, ( )
2
k

1
F
---
M

k
, ( )
2
k

=
h
F
( ) exp
1
F
--- T
k
, ( ) ( ) log
k
G

k
, ( )
Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) T
H

k
, ( )
H
k
, ( )
----------------------------------------------------------------------------------------------------------
=
M

k
, ( )
U k ( ) M
k
, ( )Y k ( ) T
M

k
, ( ) T
L

k
, ( ) + +
L
k
, ( )
------------------------------------------------------------------------------------------------------------
=
216/243
Connection with the time domain joint input-output approach
If
1.
2. covers uniformly the unit circle
then
with the dominant pole of
Notes:
- condition 1 is fulfilled if
monic noise and signal models
at least one sample delay in plant and/or controller
- identification of , is no longer coupled with that of ,
- joint input-output approach
(Ljung, 1999; Sderstrm and Stoica, 1989)
V
F
Z , ( ) h
F
( )
2
1
F
---
G

k
, ( )
2
k

1
F
---
M

k
, ( )
2
k

=
T z
1
, ( )
z
lim 1 =
h
F
( ) exp
1
N
--- - T z
k
1
, ( ) ( ) log
k 0 =
N 1

1 O
max
N
N ( ) + = =
max
T z
1
, ( ) ( ) log
H z
1
, ( )L z
1
, ( )
1 G z
1
, ( )M z
1
, ( ) +
----------------------------------------------------

log =
G H M L
h
F
( ) 1 =
217/243
8.d Measurement example: nonlinear electrical circuit
White noise excitation with bandwidth of
Measurement time 17 s: data samples at
CT-plant model in the band [0 Hz, 200 Hz] with DT- or CT-noise model
+
-
u t ( )
y t ( )
y
3
t ( )
R 16 =
C 9.4 F =
L 1 H =
1000 V
2
G
BLA
s ( )
U k ( )
Y k ( )
V k ( )
200 Hz
N 10400 = f
s
20 MHz 2
15
600 Hz =
218/243
8.d Measurement example: nonlinear electrical circuit (contd)
DT-noise model
n
c
n
d
8 = = n
c
6 = n
d
, 5 =
DT-noise model
CT-noise model
n
c
4 = n
d
, 5 =
CT-plant model
n
a
2 = n
b
, 0 =
0 50 100 150 200
100
50
0
50
Frequency (Hz)
F
R
F

(
d
B
)
G
LPM
var G
LPM
( )
G
LPM
G
0 100 200 300

150
100
50
0
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
0 50 100 150 200
150
100
50
0
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
0 50 100 150 200
150
100
50
0
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
Y G
U T
1 2 /
H
var
1 2 /
H
( )
var G
( )
G
219/243
8.e Measurement example: Villa Passo bridge
Unobserved excitations: traffic + turbulent air flow
Acceleration measurements in test points
Measurement time 14 min: data points per channel at
Model the band [1.18 Hz, 4.14 Hz]: DFT lines
14
13
12
11
10
9
8
7
6
5
4
3
2
1
N 337700 = f
s
400 Hz =
k 999 1000 3499 , , , = F 2501 =
220/243
8.e Measurement example: Villa Passo bridge (contd)
test point 1 test point 6
CT-ARMA 8/8
CT-ARMA 8/8
1 2 /
H
1 2 3 4
60
40
20
0
20
Frequency (Hz)
Y
1
(
k
)

(
d
B
)
1 2 3 4
60
40
20
0
20
Frequency (Hz)
Y
6
(
k
)

(
d
B
)
V
Y T
=
var
1 2 /
H
( )
221/243
8.f Measurement example: flight flutter
Applied excitation: perturbation force at flap right wing
Acceleration measurement at left wing
Unobserved excitation: turbulent air flow
Measurement time 109 s: data points per channel at
Model the band [4.70 Hz, 16.40 Hz]: DFT lines
N 32768 = f
s
300 Hz =
k 513 514 1791 , , , = F 1279 =
222/243
8.f Measurement example: flight flutter (contd)
4 6 8 10 12 14
40
20
0
20
40
Frequency (Hz)
F
R
F

(
d
B
)
4 6 8 10 12 14
60
40
20
0
20
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
CT-ARMAX
CT-Box-J enkins
n
a
6 = n
b
, 8 = n
c
, 10 = n
a
6 = n
b
, 8 = n
c
, 10 = n
d
, 6 =
4 6 8 10 12 14
40
20
0
20
40
Frequency (Hz)
F
R
F

(
d
B
)
4 6 8 10 12 14
60
40
20
0
20
Frequency (Hz)
N
o
i
s
e

v
a
r
.

(
d
B
)
G
LPM
var G
LPM
( )
G
LPM
G
var G
( )
G
Y G
U =

T
1 2 /
H
var
1 2 /
H
( )
223/243
8.g Multivariable systems in closed loop with known controller
with
and
U k ( )
G
k
( )
M
0

k
( )
R k ( ) Y k ( )
V k ( ) H
k
( )E k ( ) =
n
y
outputs
n
u
inputs
log f
Y k ( )
Y k ( ) R k ( ) , , ( ) ( )
k
C
1
det S
k
, ( ) ( )
2
( ) log
k
Flogdet ( ) + + + =
k
, ( )
k
, ( )
k

k
, ( ) H
1
k
, ( ) Y k ( ) G
k
, ( )U k ( ) ( ) =
S
k
, ( ) I
n
y
G
k
, ( )M
0

k
( ) + ( )
1
H
k
, ( ) =
E k ( )E
H
k ( ) { } =
n
y
n
y
224/243
8.g Multivariable systems in closed loop with known controller (contd)
Gaussian ML estimator (contd)
Eliminating gives
and
Notes:
- no problems with calculation
- numerical stable Newton-Gauss scheme of the form
(see Pintelon et al., 2007)
V
F
Z , ( ) det Re
1
F
--- g
F
( )
2

k
, ( )
H
k
, ( )
k
( ) ( ) = with g
F
( ) exp
1
n
y
F
--------- logdet S
k
, ( ) ( )
k
( ) =

k
, ( ) H
1
k
, ( ) Y k ( ) G
k
, ( )U k ( ) ( ) =
S
k
, ( ) I
n
y
G
k
, ( )M
0

k
( ) + ( )
1
H
k
, ( ) =
H
1
, ( )
J e =
225/243
Asymptotic properties ML estimator
Consistency
Open loop
all model structures
Closed loop
all model structures except Hybrid BJ
Asymptotically normally distributed
Asymptotically efficient
Robustness
consistency and asymptotic normality remain valid for non-Gaussian mixing
frequency domain noise
Note:
- consistency plant model parameters always requires the correct noise
model structure, which is not the case for the classical PE method
(see section 8.a on page 204)
226/243
Prediction error
with
- the plant transient term
- the noise transient term
Common denominator parametrization
where

k
, ( ) H
1

k
, ( ) Y k ( ) G
k
, ( )U k ( ) T
G

k
, ( ) T
H

k
, ( ) ( ) =
T
G

k
, ( )
T
H

k
, ( )
G , ( )
B , ( )
A , ( )
------------------
= T
G

k
, ( ) ,
I
G
, ( )
A , ( )
-------------------
=
H , ( )
C , ( )
D , ( )
-------------------
= T
H

k
, ( ) ,
I
H
, ( )
D , ( )
-------------------
=
A: polynomial B: n
y
n
u
polynomial matrix I
G
: n
y
1 vector polynomial , ,
C: n
y
n
y
polynomial matrix D: polynomial I
H
: n
y
1 vector polynomial , ,
227/243
8.g Model selection
Model selection criteria

with ML cost function
and penalty
Notes:
- complex models are more likely to be rejected than in case of known noise
models: compare with (see section 6.l on page 171)
- for very small data sets may perform better for MDL
2F
2F n
------------------

n
y
V
F

Z ( ) Z , ( )e
p n
F , ( )
V
F
Z , ( ) det Re
1
F
--- g
F
( )
2

k
, ( )
H
k
, ( )
k
( ) ( ) =
p n
F , ( )
2 n
1 + ( )
2F n
2
--------------------------- AIC
2n
y
F ( ) n
1 + ( ) log
2F n
2
--------------------------------------------- MDL
=
e
p n
F , ( )
1 p n
F , ( ) + ( )
V
F

Z ( ) Z , ( )e
p n
F , ( )
228/243
8.h Measurement example: Villa Passo bridge - MIMO case
1 6 7
1
6
7
CT-ARMA
H
H
n
c
10 n
d
, 8 = =
V
H
var H
1 2 /
( )
V
Y T
=
229/243
8.i Measurement example: flight flutter - MIMO case
4 5 6 7 8
40
20
0
20
Frequency (Hz)
F
R
F

(
d
B
)
4 5 6 7 8
0
50
100
150
200
Frequency (Hz)
a
n
g
e
l
(
F
R
F
)

(
)
4 5 6 7 8
40
20
0
20
Frequency (Hz)
F
R
F

(
d
B
)
4 5 6 7 8
200
100
0
100
200
Frequency (Hz)
a
n
g
e
l
(
F
R
F
)

(
)
n
a
8 = n
b
, 10 n
c
8 = , =
CT-ARMAX
G
LPM
var G
LPM
( )
G
LPM
G
var G
( )
G
230/243
8.i Measurement example: flight flutter - MIMO case (contd)
n
a
8 = n
b
, 10 n
c
8 = , =
CT-ARMAX
4 5 6 7 8
40
20
0
20
Frequency (Hz)
N
o
i
s
e

C
o
v
.

(
d
B
)
4 5 6 7 8
40
20
0
20
Frequency (Hz)
N
o
i
s
e

C
o
v
.

(
d
B
)
4 5 6 7 8
40
20
0
20
Frequency (Hz)
N
o
i
s
e

C
o
v
.

(
d
B
)
4 5 6 7 8
40
20
0
20
Frequency (Hz)
N
o
i
s
e

C
o
v
.

(
d
B
)
H
V
H
var H
1 2 /
( )
V
Y G
U T
=
231/243
8.j References
Generalized output error
Box, G. E. and G. M. Jenkins (1970). Time Series Analysis: Forecasting and Control.
Holden-Day: Oakland.
Mckelvey, T. (2000). Frequency domain identification. Preprints 12th IFAC Symposium on
System Identification, Santa Barbara, California (USA), June 21-23.
Pintelon R., and J. Schoukens (2006). Box-Jenkins identification revisited - Part I: theory,
Automatica, vol. 42, no. 1, pp. 63-75.
Pintelon R., Y. Rolain, and J. Schoukens (2006). Box-Jenkins identification revisited - Part
II: applications, Automatica, vol. 42, no. 1, pp. 77-84.
Pintelon, R., J. Schoukens, and P. Guillaume (2007). Box-Jenkins identification revisited -
Part III: multivariable systems, Automatica, vol. 43, no. 5, pp. 868-875.
Pintelon, R., J. Schoukens, and P. Guillaume (2006). Continuous-time noise modelling from
sampled data, IEEE Trans. Istrum. Meas., vol. 55, no. 6, pp. 2253-2258.
Errors-in-variables
Pintelon, R., and J. Schoukens (2007). Frequency domain maximum likelihood estimation
of linear dynamic errors-in-variables models, Automatica, vol. 43, no. 4, pp. 621-630.
Sderstrm, T. (2007). Errors-in-variables methods in system identification, Automatica, vol.
43, no. 6, pp. 939-958.
232/243
9. Guidelines
9.a Two basic questions 233
9.b Identification step by step 235
9.c References 241
233/243
9. Guidelines
9.a Two basic questions
Stochastic framework? Output observations only Errors-in-variables
Known input - output disturbed Output error - indep.
by unobserved input(s) plant and noise dyn.
Noise model Parametric Nonparametric
Identification method Prediction error Sample maximum
Frequency domain maximum likelihood
likelihood
234/243
9. Guidelines
9.a Two basic questions (contd)
Application? Digital control Physical interpretation Digital simulation
Analog simulation
Exp. setup Zero-order-hold Band-limited Band-limited
Domain Discrete-time Continuous-time Continuous-time
235/243
9. Guidelines
9.b Identification step by step
Check and selection of the experimental setup
Design of an experiment
Choice noise model
Preprocessing of the data
The identification step
236/243
9. Guidelines
9.b Identification step by step (contd)
Check and selection of the experimental setup
Visit the site of the experimental setup and talk to the operators to learn from their
experience
Check the systematic errors of the complete data acquisition - calibrate the setup
Check the validity of the signal assumptions (ZOH or BL, anti-alias filters)
Pay attention to the synchronisation of the setup
Always save the signal stored in the arbitrary waveform generator together with the
noisy input-output data
237/243
9. Guidelines
Design of an experiment
Choose the excitation (random or periodic) that best fits your specific needs
Select the amplitude range and frequency band of the excitation signal to cover the
frequency band of interest
Check for the presence of nonlinear distortions
Check whether the device is captured in a feedback loop
Keep your application in mind
238/243
9. Guidelines
Choice noise model
Default choice: nonparametric models
Use parametric noise models for the following problems:
- output data only
- known input / output disturbed by unobserved random input(s)
239/243
9. Guidelines
Preprocessing of the data
Removal of trends, drifts, and offsets
- Check for trends by calculating the mean value over the successive periods
- Do not use the DC information during the identification
Dealing with outliers and missing data
- Perform new experiments
- If this is not possible, use simple interpolation methods if
- Last resort: estimate the missing data
Estimate the nonparametric FRF
- Calculate the FRF + noise level + level nonlinear distortions
- Select the frequency band of interest
Check whether the system is time invariant
- Check the presence of skirts in the output DFT spectrum of all samples
f
max
0.1f
s
<
240/243
9. Guidelines
The identification step
Choice of a model class
- , or domain
Model selection/ validation
- Compare the transfer function model to the nonparametric FRF estimate
- Analyse the value of the cost function
- Perform a whiteness test of the residuals
- Detect overmodeling via AIC (prediction) or MDL (physical interpretation)
Impact of transient/leakage errors
- Suppress nonparametrically the transient (leakage) errors in the data
- If the transient time can be neglected, measure the steady state response
- Always add the plant and noise transient terms when estimating parametric
noise models
- Suppress (non)parametrically the transient (leakage) errors in the validation
data
s- z- s-
241/243
9. Guidelines
9.c References
Broersen P. M. T., de Waele S., Bos R. (2004). Application of autoregressive spectral
analysis to missing data problems, IEEE Trans. on Instr. Meas., vol. 53, no. 4, pp. 981-986.
Broersen P. M. T., de Waele S., Bos R. (2004). Autoregressive spectral analysis when
observations are missing, Automatica, vol. 40, no. 9, pp.1495-1504.
Isaksson, A. J. (1993). Identification of ARX-models subject to missing data, IEEE Trans. on
Automatic Control, vol. 38, no 5, pp. 813-819.
Little, R. J. A. and D. B. Rubin (1987). Statistical Analysis with Missing Data. John Wiley,
New York.
Peirlinckx L., P. Guillaume and R. Pintelon (1996). Accurate and Fast Estimation of the
Fourier Coefficients of Periodic Signals Disturbed by Trends, IEEE Trans. Instrum. Meas.,
vol. IM-45, no. 1, pp. 5-11.
Pintelon, R. and J. Schoukens (1999). Identification of Continuous-Time Systems with
Missing Data, IEEE Trans. on Instr. Meas., vol. IM-48, no 3, pp. 736-740.
Pintelon, R. and J. Schoukens (2000). Frequency domain system identification with missing
data. IEEE Trans. on Autom.Contr., vol. AC-45, no 5, pp. 364-369.
242/243
10. Books on system identification
strm, K. J. (1970). Introduction to Stochastic Control Theory. Academic Press: New York.
Box, G. E. and G. M. Jenkins (1970). Time Series Analysis: Forecasting and Control.
Holden-Day: Oakland.
Caines, P. E. (1988). Linear Stochastic Systems. Wiley: New York.
Eykhoff, P. (1974). System Identification, Parameter and State Estimation. Wiley, New York.
Godfrey, K. R., editor (1993). Perturbation Signals for System Identification. Prentice-Hall,
London.
Goodwin, G. C. and R. L. Payne (1977). Dynamic System Identification. Experiment Design
and Data Analysis. Academic Press: New York.
Kollr, I. (1994). Frequency-Domain System Identification Toolbox for Use with Matlab. The
Mathworks, Natick.
Ljung, L. (1995). System Identification Toolbox Users Guide. The MathWorks, Natick.
Ljung, L. and T. Sderstrm (1983). Theory and Practice of Recursive Identification. MIT
Press, Cambridge, MA.
Middleton, R. H. and G. C. Goodwin (1990). Digital Control and Estimation. Prentice-Hall:
London.
243/243
Sorenson, H. W. (1980). Parameter Estimation: Principles and Problems. Marcel Dekker,
New York.
Zarrop, M. B. (1979). Optimal experiment design for dynamic system identification. Series
of Lecture Notes in Control and Information Sciences, vol. 21, Springer-Verlag, Berlin.

Gevorderde Systeemidentificatie

Uploaded by

Copyright:

Available Formats

Gevorderde Systeemidentificatie

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Gevorderde Systeemidentificatie

Uploaded by

Copyright:

Available Formats

1/243

LPM and SA (hanning window): same dof noise covariance estimates

number of free parameters =

----------------------------------------------- MDL (EIV)

------------------------------- MDL (OE)

number of free parameters =

0 100 200 300

You might also like