Frequency Domain Identication, Subspace methods and Periodic Excitation
T. McKelvey
Dept. of ElectricalEngineering
Linkoping University
S-581 83Linkoping, Sweden
tomas@isy.liu.se
Presented at the 2nd Russian-Swedish Control Conference, Sankt Petersburg, Russia, August 1995.
Technical Report LiTH-ISY-R-1765
Abstract
Recent frequency domain identication algorithms based on subspace based techniques are discussed. The algorithms construct a state-space model by means of extraction of the signal subspace from a matrix constructed from frequency data. A singular value decomposition plays a key part in the subspace extraction. The subspace methods are non-iterative methods in contrast to classical iterative parametric optimization techniques. The use of periodic excitation leads to a leakage free discrete Fourier transform of the measured data as well as simple noise reduction possibilities by averaging.
1 Introduction
System identi cation is the technique of tting a linear dynamic model to measured input and output data. Most often the model tted is represented by the fraction of two polynomials (in the forward shift operator
qor the
Z-transform variable
z) and the polynomial coecients constitutes the parameters. The values of the parameters are found by minimizing a criterion function, often represented by the sum of the squares of the prediction errors 4]. This is known as prediction error minimization (PEM). For most model structures the criterion is a complicated function of the parameters and iterative optimization techniques have to be applied. Some drawbacks with this approach can be identi ed: The iterative optimization can lead to convergence problems with local minima and has a high computational complexity. These properties become more pronounced as the model order increases and the system has lightly damped modes, e.g., a exible mechanical structure. A polynomial representation of a high order system is also known to be numerical sensitive with respect to small perturbations in the polynomial coecients.
During the last decade an alternative branch of identi cation methods for multi-input multi- output systems has emerged. These are commonly known as subspace methods since the key step in all methods is the extraction of a signal subspace from a matrix constructed from measured data.
From this signal subspace a state-space model is constructed by utilizing the geometrical properties of the signal subspace. The subspace methods use robust numerical algebra methods, the singular value decomposition (SVD) and QR-factorization. The identi ed state-space model is derived without any iterations except in the SVD calculation. The state-space model is delivered in a state-space basis which is less sensitive to perturbations in the matrix elements compared with the polynomial
This work was supportedby the SwedishResearch Council for Engineering Sciences (TFR),which is gratefully
acknowledged.
1
representation. The lack of iterations and the numerical robustness of the calculations as well as the representation has opened up the possibilities to identify high order systems with lightly damped modes. The most recent time-domain identi cation methods can be found in Van Overschee and De Moor 9] and Verhaegen 11]. A comprehensive overview of the time domain subspace identi cation are is given by Viberg 12].
This paper will highlight recent subspace based identi cation algorithms for use with frequency data.
Preliminaries A discrete time linear system of minimal order
ncan be represented by a state-space model
x
t
+1=
Axt +
But
y
t =
Cxt +
Dut +
vt (1)
where
ut
2 Rm ,
xt
2 Rn ,
yt
2 Rp and
vt
2 Rp are the input, state, output and noise signals, respectively. Since the system is of minimal dimension, the extended observability matrix
O
q =
CT (
CA) T
(
CAq
;1) T
T
2Rqp
n
q>n(2) and the extended controllability matrix
C
r =
B AB Ar
;1B 2Rn
rm
rn(3) both have full rank
n. The transfer function of (1) is given by
G
(
z) =
D+
C(
zI;A)
;1B(4)
2 Subspace Identication
Let us assume that a matrix
Hqr
2 Rpq
mr can be constructed from data such that the following relation holds
H
qr =
Oq
X+
V(5)
where
Oq
Xrepresents the signal part and
Vis a matrix composed of the noise. Consider the SVD
H
qr =
Us
Uo
s 0 0 o
!
V
Ts
V
To
!
(6) where
Us
2Rpq
n and
Vs
2Rmr
n are orthonormal matrices and s
2 Rn
n is a diagonal matrix with the
nlargest singular values. If
kVkF is small compared to the smallest non-zero singular value of
Oq
X, we have approximately
Oq
X Us s
VTs . The observability range space is then accurately given by
Us and we can use it as an estimate of the observability matrix. In the noise free case
V
= 0,
Us =
Oq
Tfor some non-singular matrix
Twhich shows that
Us is the extended observability matrix of some realization of the system (1). Let ( ^
AB^
C^
D^ ) denote the the state-space realization associated with
Us . Then the rst
prows of
Us is equal to ^
C. Denote by
Us the sub-matrix obtained by deleting the
prst rows and denote by
Us the corresponding matrix obtained by deleting the
plast rows. Then in the noise free case
Us
A^ =
Us and a natural estimate of ^
Ais
^
A
=
Uys
Us
:(7)
Given ^
Cand ^
Athe determination of ^
Band ^
Dis straight forward. The transfer function
G(
z) (4) as well as the output
yt is a linear function in ^
Band ^
D. A simple choice is a minimization in the least-square sense
^
BD
^ = arg min BD
Xt
ky
t
;y^ t
k2F
or ^
BD^ = arg min BD
Xt
kG
(
zt )
;G^ (
zt )
k2F (8)
where
G(
zt ) and
yt are the measured quantities and ^
G(
zt ) and ^
yt are the model quantities. Notice that (8) has an analytical solution.
3 Frequency Domain Algorithms
3.1 Identication from Transfer Function Measurements
Assume that
Msamples
Gk of the transfer function
G(
z) is given at equidistant frequencies covering the full unit circle,
z=
ej
2k=M ,
k= 0
:::M;1. Construct the coecients
hk by use of the inverse discrete Fourier transform
h
k = 1
M
M
X;1s
=0 Gs
ej
2sk=M
:Let the Hankel matrix
Hqr be de ned by
H
qr ] st =
hs
+t
;1 s= 1
:::qand
t= 1
:::r:Straight forward calculations 6] reveal that
H
qr =
Oq (
I;AM )
;1Cr +
Vwhere
Vis the Hankel matrix of the noise contribution. A state-space model can now be estimated along the steps presented in the previous section.
3.2 Identication from the Transform of the Input and Output Signals
If, at arbitrary frequencies, the discrete Fourier transform is given for the input and output signals respectively, a matrix
Hqr with the properties (5) can be constructed by the use of some matrix projections 3, 6, 5].
Let us introduce the notation
W
(
!) = 1
ej!
ej
2!
ej
(;1)
!
T
and the lower Toeplitz matrix
; =
0
B
B
B
B
@
D
0
:::0
CB D :::
0
... ... ... ...
CA
;2B CA;3B ::: D 1C
C
C
C
A
:
(9)
and let
Y(
!),
U(
!) and
V(
!) denote the discrete Fourier transform of
yt ,
ut and
vt respectively. By recursive use of the state-space equations (1) we obtain
W
(
!)
Y(
!) =
OX(
!) + ;
W(
!)
U(
!) +
W(
!)
V(
!) (10) where
denotes the Kronecker product.
Assume we have samples of the Fourier transform of the input and output signals at
Mfrequencies
!
k ,
k= 1
::: M. Collect these samples in matrices
Y =
W(
!1)
Y(
!1)
W(
!M )
Y(
!M )
2Cp
M
(11)
U =
W(
!1)
U(
!1)
W(
!M )
U(
!M )
2Cm
M
(12)
V =
W(
!1)
V(
!1)
W(
!M )
V(
!M )
2Cp
M
(13)
X =
X(
!1)
X(
!M )
2Cn
M
:(14)
From (10) we arrive at the matrix equation
Y =
OX + ; U + V
:(15) The rank of the matrix U is a measure of the excitation of the input signal. In the time domain this represents the concept of persistence of excitation 4]. In the single input case U is of full rank whenever
U(
!k ) is non-zero for at least
distinct frequencies. The characterization of the rank condition for the multi-input case is more involved and from a user's point of view it suces to check that U is of full rank for a given set of data.
By use of a matrix projection we will construct a matrix with the desired properties (5). A simple algorithm was proposed in 3]. Analysis in 7] show that this simple algorithm is consistent only under restrictive noise assumptions. Here we will present an algorithm adopted from the time domain technique of instrumental variables 10]. This leads to an algorithm with improved consistency properties.
3.2.1 Matrix Projections
Partition Y , U and V as
Y =
Y p
Y f
!
U =
U p
U f
!
V =
V p
V f
!
:
By following the nomenclature of the time domain case, we call U p and Y p the \past" inputs and outputs, respectively, and conformally U f and Y f , the \future" inputs and outputs. The partition is done such that each sub-matrix will have
qblock rows and consequently 2
q=
block rows. The size requirement of this partition is that
q >nwhere
nis the system order. It is straightforward to show that
Y f =
Oq X f + ; q U f + V f
(16) where X f is given by
X f =
ejq!
1X(
!1)
ejq!
2X(
!2)
::: ejq!
MX(
!M )
:We now have the possibility to remove the future inputs by a projection and then use the past inputs in order to remove the noise 10]. Using (16) we get
Y f
?UHf
U Hp =
Oq X f
?UHf
U Hp + V f
?UHf
U Hp (17)
where
?UHf
denotes the orthogonal projection onto the null-space of U f and is given by
?UHf
=
I;U Hf ( U f U Hf )
;1U f
:(18)
The inverse in (18) exists whenever U f has full rank. If we assume
V(
!k ) to be zero mean independent random variables with uniformly bounded second moments
EV
(
!k )
V(
!k ) H =
R(
!k )
R 8!k
and assume that
U(
!k ) to be uniformly bounded and that U f has full rank, the following relation follows from a standard limit result 2, Theorem 5.1.2]
M lim
!11
M
N f
?UHf
U Hp = 0
w.p. 1
:Hence we can use
H
qr = Y f
?UHf