IEEE
TRANSACTIONS ON AUTOMATIC CONTROL, VOL.
40,
NO.
3,
MARCH
1995
45
1
A
Generalized Orthonormal Basis
for Linear Dynamical Systems
Peter
S.
C.
Heuberger, Paul
M.
J.
Van
den
Hof,
Member,
IEEE,
and
Okko
H. Bosgra
Abstract-In many areas of signal, system, and control theory,
orthogonal functions play an important role in issues of analysis
and design. In this paper, it is shown that there exist orthogonal
functions that, in a natural way, are generated by stable linear
dynamical systems and that compose an orthonormal basis for
the signal space
e;.
To this end, use is made of balanced re-
alizations of inner transfer functions. The orthogonal functions
can be considered as generalizations
of,
e.g., the pulse functions,
Laguerre functions, and Kautz functions, and give rise to an
alternative series expansion of rational transfer functions. It
is shown how we can exploit these generalized basis functions
to increase the speed of convergence in a series expansion,
i.e., to obtain a good approximation by retaining only a finite
number of expansion coefficients. Consequences for identification
of expansion coefficients are analyzed, and a bound is formulated
on the error that is made when approximating a system by a finite
number of expansion coefficients.
I.
INTRODUCTION
ONSIDER
a linear time-invariant stable discrete-time
C
system
G,
represented by its proper transfer function
G(z)
in the Hilbert space
?fa,
i.e.,
G(z)
is analytic outside
the unit circle,
IzI
2
1.
A
general and common representation
of
G(z)
is in terms of its Laurent expansion around
z
=
CO,
as
k=O
with
{Gk}k=O,l,
...
the sequence of Markov parameters.
In constructing this series expansion we have employed a
set of orthogonal functions:
{zo,
z-l,
z-',
. .
-},
where orthog-
onality is considered in terms of the inner product in
?fa.
In
a generalized form we can write
(1)
as
m
k=O
with
{fk(z)}k=O,J2,
...
a sequence of orthogonal functions.
Manuscript received November 17, 1992; revised December
3,
1993.
recommended by Past Associate Editor, B. Pasik-Duncan.
This
work was
supported in part by Shell Research B.V., The Hague, and the Center
for
Industrial Control Science, The University of Newcastle, Newcastle NSW,
Australia.
P.
S.
C. Heuberger was with the Mechanical Engineering Systems and
Control Group, Delft University of Technology, 2 628 CD Delft, The
Netherlands and is now with the Dutch National Institute of Public Health
and Environmental Protection (RIVM),
P.O.
Box 1, 3720 BA Bilthoven, The
Netherlands.
P.
M.
J.
Van den Hof and
0.
H.
Bosgra are with the Mechanical Engineering
Systems and Control Group, Delft University
of
Technology, Mekelweg 2,
2
628
CD Delft, The Netherlands.
IEEE Log Number 9408269.
There are a number of research areas that deal with the
question of either approximating a given system
G
with a
finite number of coefficients in a series expansion as in (2), or
(approximately) identifying an unknown system in terms of a
finite number of expansion coefficients through
N
e'(.)
=
Lkfk(2).
(3)
k=O
The problem that will be analyzed in this paper is the
Can we construct a sequence of orthogonal basis functions
a) to some extent, the basis can be adapted to a linear
stable system
G
to be described, implying that
G
can
be accurately described by only a small number of
coefficients in the expansion, and
b) the basis allows the construction of an error bound for
the approximation of a linear stable system
G
by a finite
length expansion in the basis
fk,G,
i.e., an upper bound
on
ll~(z)-~~=o~k~k,G(r)ll
in some prechosen norm,
whenever
G
and
G
do not match exactly.
The use of orthogonal functions with the aim of adapting
the system and signal representation to the specific properties
of the systems and signals at hand has a long history. The
classical work of Lee and Wiener during the 1930's on network
synthesis in terms of Laguerre functions [24], [46] is summa-
rized in [25]. Laguerre functions have been used in the 1950's
and 1960's to represent transient signals [45], [7]. During
the past decades, the use of orthogonal functions has been
studied in problems of filter synthesis [22], [30] and for system
identification [23], [32], [31], [6] and approximation [35], [36].
In these approaches to system identification, the input and
output signals are transformed to a (Laguerre) transformed
domain and standard identification techniques are applied to
the signals in this domain. Data reduction has been the main
motivation in these studies. Identification of continuous-time
models with the aid of orthogonal functions is considered
in e.g., [38] and [29]. In recent years, a renewed interest
in Laguerre functions has emerged. The approximation of
(infinite dimensional) systems in terms of Laguerre functions
has been considered in [27], [28], [12], [13], and [15]. In the
identification of coefficients in finite length series expansions,
Laguerre function representations have been considered from
a statistical analysis point of view in [43], [42], and [16].
following.
{fk,G(z)}k=O,
...
with
G
E
?fa,
such that
0018-9286/95$04.00
0
1995 IEEE
1
.-
I
111
I
"
I
452
IEEE
TRANSACTIONS ON AUTOMATIC CONTROL,
VOL.
40,
NO.
3,
MARCH
1995
The use of Laguerre-function-based identification in adaptive
control and controller tuning is studied in
[47]
and
[9].
A
second-order extension to the basic Laguerre functions using
the
so
called Kautz functions
[21]
is subject of discussion in
[41]
and
[44].
In this paper we will expand and generalize the orthogonal
functions as basis functions for dynamical system representa-
tions. Specifically we will generalize the Laguerre functions
and Kautz functions to a situation where a higher degree of
flexibility is present in the choice of basis functions, and
where consequently a smaller error bound as meant in part
b) of the problem can be obtained. Laguerre functions are
specifically appropriate for accurate modeling of systems with
dominant first-order dynamics, whereas Kautz functions are
directed toward systems with dominant second-order resonant
dynamics. The generalized basis functions, introduced in this
paper, will
be
suited also for systems with a wide range of
dominant dynamics, i.e., dominant high frequency and low
frequency behavior.
We will restrict attention to the transfer function space
7&
being equipped with the usual inner product. This choice,
rather than the &-space where orthogonality is abandoned,
is motivated by the fact that our main intended application of
these results is in the area of approximate system identification.
As the main stream of approaches in system identification is
directed toward prediction error methods and the use of least-
squares types of identification criteria,
[26],
the choice of a
two-norm is quite straightforward and natural in this respect.
Note that the two problems a) and b) should
be
treated as a
joint problem. One of the (trivial) solutions to problem a) only
is the use of a Gram-Schmidt orthogonalization procedure on
the impulse response of the system
G
itself
[l].
In that case
the system can be described by a series expansion of only one
single term. In this situation, however, no results are available
for part b) of the problem.
In an identification context, the use of the orthogonal
functions as in
(1)
leads to the so-called finite impulse response
(FIR)-model
[26]
where
~(t)
is the one-step-ahead prediction error, and
{y(t),
u(t)}
are samples of the output, input of the dynamical
system to be identified. The identification of the unknown co-
efficients
{
Gk
(B)}k=O,
...,N
through least squares minimization
of
E(t)
over the time interval is an identification method that
has some favorable properties. First, it is a linear regression
scheme, which leads to a simple analytical solution: second, it
is
of
the type of output-error-method, which has the advantage
that the input/output system
G(z)
can be estimated consistently
whenever the unknown noise disturbance on the output data
is uncorrelated with the input signal
[26].
It is well known, however, that for moderately damped
systems, and/or in situations of high sampling rates, it may take
a large value of
N,
the number coefficients to be estimated, to
capture the essential dynamics of the system
G
into its model.
If we would
be
able to improve the basis functions in such a
way that an accurate description of the model to
be
estimated
can
be
achieved by a small number of coefficients in a series
expansion, then this is beneficial from both aspects of bias and
variance of the model estimate.
For the series expansion in
(1)
with
fk
=
z-',
it is
straightforward to show that a system
G
will have a finite
length series expansion if and only if all system poles are at
z
=
0.
Moreover, in the scalar case the length of the expansion,
i.e., the index of the last nonzero coefficient, equals the total
number of poles at
z
=
0.
As a generalized situation, we can consider Laguerre poly-
nomials
[37]
that are known to generate a sequence of
orthog-
onal functions
[14]
Similar to above, a system
G
will have a finite length series
expansion if and only if all system poles are at
z
=
a,
with
the length of the expansion being equal to the total number
of poles at
z
=
a.
In dealing with the problem
of
finding similar results for any
general stable dynamical system
G(z),
we have considered
the question of whether a linear system in a natural way
gives rise to a set of orthogonal functions. The answer to
this question appears to be affirmative. It will be shown that
every stable system gives rise to a complete set of orthonormal
functions based on input (or output) balanced realizations, or
equivalently based on a singular value decomposition of a
corresponding Hankel matrix. These generalized orthogonal
basis functions will be shown to provide solutions to problems
a) and b).
In Section I11 we will first briefly state the main result of
this paper. Next in Section IV it will be shown how inner
functions generate two sets of orthonormal functions that are
complete in the signal space
l2.
This is the basic ingredient of
the main result. Next an interpretation of these results is given
in terms of balanced state-space representations. After showing
the relations of the new basis functions with existing ones, we
will focus on the dynamics that implicitly are involved in the
inner functions generating the basis. It will be shown that if the
dynamics of a stable system match the dynamics of the inner
function that generates the basis, then the representation of
this system in terms of this basis becomes extremely simple.
Consequences for a related identification and approximation
problem are discussed in Section VIII.
Due to space limitations, a complete statistical analysis of
the related system identification problems that result from
these basis functions can not
be
given in this paper.
A
statistical analysis along similar lines as
[43]
and
[44]
is
presented elsewhere
[39].
The proofs of all results are collected in an appendix.
11.
PRELIMINARIES
We will use the following notation.
Transpose of a matrix.
HEUBERGER
et
al.:
ORTHONORMAL BASIS FOR LINEAR DYNAMICAL SYSTEMS
~
453
Complex conjugate transpose of a matrix.
Set of complex-valued matrices of
dimension
p
x
m.
Real-valued matrix with dimension
p
x
m.
Set of nonnegative integers.
Space of squared summable sequences on
the time interval
Z+.
Space of matrix sequences
{Fk
E
Cmxn}],=0,1,2,
...
such that
tr(F,*Fk)
is finite.
Set of real
p
x
m
matrix functions,
analytic for
integrable on the unit circle.
Set of real rational
p
x
m
matrix
functions, analytic for
IzI
2
1,
that are
squared integrable on the unit circle.
Induced 2-norm or spectral norm of a
constant matrix, i.e., its maximum singular
value.
H,-norm.
Vector-operation on a matrix, stacking its
columns on top of each other.
Kronecker matrix product.
(Block) Hankel matrix related to transfer
function
G
=
Xij(G)
=
Gi+j-l
being the (i,j)-block
element.
ith Euclidian basis vector in
R".
n
x
n
Identity matrix.
2
1,
that are squared
Gkz-k,
defined by
In this paper we will consider discrete-time signals and
systems.
A
linear time-invariant finite-dimensional system will
be represented by its rational transfer function
G
E
RX;'",
with m the number of inputs in
U,
and
p
the number of outputs
in
y.
State-space realizations will
be
considered of the form
~(k
+
1)
=
Ax(k)
+
Bu(k)
y(k)
=
Cz(k)
+
Du(k)
(6)
(7)
with
A
E
CYXn,
B
E
Cnxm,
C
E
CpXn,
and
D
E
Cpxm.
(A,
B,
C,
D)
is an n-dimensional realization of
G
if
G(z)
=
C(zI
-
A)-lB
+
D.
A
realization is stable if all eigenvalues
of
A
lie strictly within the unit circle. If a realization is
stable, the controllability gramian
P
and observability gramian
Q
are defined as the solutions to the Lyapunov equations
APA*
+
BB*
=
P
and
A*QA
+
C*C
=
Q,
respectively.
A
stable realization is called (intemally) balanced if
P
=
Q
=
C,
with
C
=
diag(a1,.
. .
,
an),
01
2
.
.
.
2
on,
a diagonal matrix
with the positive Hankel singular values as diagonal elements.
A
stable realization
is
called input balanced if
P
=
I,
Q
=
C2,
and output balanced if
P
=
C2,
Q
=
I.
A
system
G
E
RX;Xm
is called inner if it satisfies
GT(z-l)G(z)
=
I.
As
G
is analytic outside and on the unit
circle, it has a Laurent series expansion
Er=o
GkZ-'.
111.
THE
MAIN RESULT
We will start the technical part of this paper by giving the
basic result first and then consecutively give the analysis that
provides the ingredients for making the result plausible.
Theorem
3.1:
Let
G
be an
m
x
m inner transfer function
with McMillan degree
n
>
0,
having a Laurent expansion
G(z)
=
zr=oGkz-k
and satisfying
I1Goll2
<
1,
and let
(A, B,
C,
D)
be a balanced realization of
G(z).
Denote
Vk(z)
=
z(z1
-
A)-'BGk(z).
(8)
Then the set of functions
{eTVk(Z)}i=l, ...n;k=O,...oo
consti-
tutes an orthonormal basis of the function space
Xixm.
A
direct consequence of this theorem is the following
corollay
.
0
Corollary
3.2:
Let
G
be an inner function with McMillan
degree
n
as in Theorem
3.1,
with a corresponding sequence
of basis functions
vk(Z).
Then for every proper stable transfer
function
H
E
there exist unique
D,
E
Rpxm,
and
L
{Lk}k=O,l,
...
E
l;'""[O,CO),
such that
a2
H(Z)
=
D,
+
2-l
Lkvk(z).
(9)
of
H(z).
0
k=O
We refer to
D,,
Lk
as the orthogonal expansion coefficients
Note that due to the fact that
&(Z)
is an
n
x
m-matrix of
transfer functions, the dimension of each
Lk
is
p
x
n.
IV. ORTHONORMAL FUNCTIONS GENERATED
BY
INNER TRANSFER
FUNCTIONS
In this section we will show that a square and inner transfer
function gives rise to an infinite set of orthonormal functions.
This derivation is based on the fact that a singular value
decomposition of the Hankel matrix associated to a linear
system induces a set of left (right) singular vectors that are
orthogonal. Considering the left (right) singular vectors as
discrete time functions, they are known to be orthogonal in
&-
sense, thus generating a number of orthogonal functions being
equal to the McMillan degree of the corresponding system.
We will embed an inner function with McMillan degree
n
into a sequence of inner functions with McMillan degree
kn,
for which the left (right) singular vectors of the Hankel matrix
span a space with dimension
kn.
If we let
k
-+
00
the set
of left (right) singular vectors will yield an infinite number of
orthonormal functions, which can be shown to be complete
in
l2.
First we have to recapitulate some properties of inner
transfer functions.
Proposition
4.1:
Let
G(z)
be an inner transfer function
with a Laurent expansion
G(z)
=
Gkz-'.
Then
00
Gc+'+,Gk
=
I
for
i
=
0;
(10)
k=O
=
0
for
i
>
0.
(1 1)
0
The Hankel matrix of an inner transfer function has some
specific properties, reflected in the following two results.
I
I'
1
'I
1
454
IEEE
TRANSACTIONS ON AUTOMATIC CONTROL, VOL.
40,
NO.
3,
MARCH
1995
-Go G1 G2
.-.-
0
Go G1 G2
..'
0
0
Go GI
...
;
0
Go
**a
..
T,
=
.
. .
..
.....
....
-Go
0 0
...
...-
G1
Go
0
... ...
T,
=
G2 G1 Go
...
...
..
G1 Go
.'.
.
.
-.
.
.....
...
..
Proposition 4.2:
Let
G(z)
be an inner function with
McMillan degree
n
>
0.
Then a singular value decomposition
constitute a singular value decomposition of
X(
Gk),
through
,
(12)
(13)
.
(svd) of
X(G)
satisfies
X(G)
=
UoV;
'With slight abuse
of
notation we will use this notation
to
indicate
an
operator
C"
-+
ez(0,
CO).
x(G~)
=
rgrc,.
(16)
The matrix sequence
{U;,
&};=0,1,
...
is unique up to
postmultiplication of each
U;
and
V,
with one and the
same unitary matrix.
Let
G(z)
have a Laurent expansion
G(z)
=
czo
Giz-;,
and consider the block Toeplitz matrices
Tu,
T,
as in (12), (13) then the matrix sequence
{U;,
V,};=o,l
,...
satisfies
V,*
=
Vc-lT,,
(17)
Uk
=
TUuk-1
for
k
=
1,2,.
.
f
.
(18)
0
The theorem shows the construction of orthogonal matrices
rg,r;
that have a nesting structure. The suggested svd of
X(Gk)
incorporates svd's of
X(Gi)
for all
i
<
k.
In
this
way orthogonal matrices and
r;
are constructed with an
increasing rank. Note that the restriction on the structure of
the consecutive svd's is
so
strong that, according to b), given
a singular value decomposition
X(G)
=
UoV<,
the matrix
sequence
{U;
,
,
i
=
1,2,
.
.
.}
is uniquely determined. Note
also that there is a clear duality between the controllability
part
I?;
and the observability part
rg.
To
keep the exposition
and the notation as simple as possible we will further restrict
attention to the controllability part of the problem. Dual results
exist for the observability part.
Proposition
4.6:
Let
G(z)
be an
m
x
m
inner function with
McMillan degree
n
>
0,
and consider any sequence of unitary
matrices
{V,}i=o,l,
.,.
satisfying (17) in Theorem 4.5. Denote
for
k
E
Z+
m
Vk(2)
=
Mk(z)z-;,
with
Mk(2)
E
c"'"
defined by
i=O
v,*
=:
[hfk(o)
Mk(1) Mk(2)
.'.I.
(19)
Then
Vk(.)
=
V0(,Z)Gk((Z).
0
The proposition actually is a z-transform-equivalent of
the result in Theorem 4.5. It shows the construction of the
controllability matrix
r;.
In the next stage we show that this controllability matrix
generates a sequence of orthogonal functions that is complete
in
e;.
Theorem
4.7:
Let
G(z)
be
an
m
x
m
inner function with
McMillan degree
n
>
0,
such that
llGoll2
<
1;
consider a
sequence of unitary matrices
{&}i=o,l,
...
as meant in Theorem
4.5. For each
k
E
Z+
consider the function
:
Z+
4
c",
defined by
[4k(O)
4k(l)
4k(2)
'.'
]
=
v,*.
Then the set of functions
Q(G)
:=
{4k}&
constitutes an
0
orthonormal basis of the signal space
Cz[O,
m).
1
I
I1
1
1
HEUBERGER
et
al.:
ORTHONORMAL
BASIS
FOR LINEAR DYNAMICAL SYSTEMS
Fig.
1.
Sene
basis
Ik(G).
expansion of
a
transfer function in terms of
an
orthonormal
Remark
4.8:
This basis has been derived from the singular
value decomposition of the Hankel matrix X(G).
As
stated in
Proposition 4.2 this svd is unique up to postmultiplication of
UO,
VO
with a unitary matrix. Consequently-within this con-
text-both
VL, Vk(Z),
and the corresponding basis functions
{
$hk}
are unique up to unitary premultiplication.
For use later on we will formalize the class of inner
functions that have the property as mentioned in the previous
theorem.
Definition
4.9:
We define the class of functions:
91
:=
{all
square inner functions G with McMillan degree
>
0
such that
As
a result of the fact that the proposed orthonormal
functions constitute a basis of
lz,
each square inner function
generates an orthonormal basis that provides a unique transfor-
mation of lg-signals to an orthogonal domain. Similarly, when
given such an orthonormal basis, each stable rational function
can
be
expanded in a series expansion of basis functions
vk
(2)
as defined in Proposition 4.6.
Corollary4.10:
Let G
E
El,
and let 9(G)
be
as defined
in Theorem 4.7. Then
a) For every signal
x
E
lz[O,m)
there exists a unique
transform
x
=
{Xk}k=O,l,
...
E
!gx"[o,
m)
such that
llGoll2
<
1).
0
00
x(t)
=
Xk$hk(t).
(20)
k=O
We refer to
XI,
E
cqx"
as the orthogonal expansion
coefficients of
x.
b) For every proper stable transfer function
H(z)
E
X;',
there exist unique
D,
E
Rpxm,
and
L
=
{Lk}k=o,l,
...
E
l;""[O,
m),
such that
00
H(Z)
=
D,
f
Z-l
Lkvk(Z).
(21)
We refer to
D,, Lk
as the orthogonal expansion coeffi-
cients of
H(z).
0
We will refer to the sequence
{Vk(Z)}k=O,l,
...,
as defined
in Proposition 4.6, as the sequence of generating transfer
functions for the orthonormal basis 9(G).
The series expansion as reflected in (21) is schematically
depicted in the diagram in Fig. 1, where
q
reflects the time
To
find appropriate ways to calculate the orthogonal func-
tions, as well as to determine the transformations as meant in
the corollary, we will now first analyze the results presented
so
far in terms of state-space realizations.
k=O
shift,
qu(t)
=
u(t
+
1).
455
V.
BALANCED STATE-SPACE REPRESENTATIONS
To
represent the orthogonal controllability matrix in a state-
space form, we will use a balanced state-space realization of G.
We first present the following, rather straightforward, lemma.
Lemma
5.1:
Let G
be
a square inner transfer function
with minimal realization
(A,
B,
C,
D).
Then the realization
is (internally) balanced if and only if it is both input balanced
Next we examine how the property that a transfer function
is inner, is reflected in a state-space realization of the function.
Proposition
5.2:
Let G be a transfer function with realiza-
tion
(A,
B,
C,
D),
such that
(A,
B)
is a controllable pair, and
the realization is output balanced, i.e.,
A*A+
C*C
=
I.
Then
GT(zP1)G(z)
=
I
if and only if
i)
D*C
+
B*A
=
0,
and
ii)
D*D
+
B*B
=
I.
0
Note that for this proposition there also exists a dual
form, conceming the transfer function
GT
with realization
(A*,
C*,
B*, D*),
that can
be
applied if G is square inner.
The characterization of the inner property in the above
proposition is made for output balanced realizations. Since,
according to Lemma
5.1
output balancedness is implied by
balancedness, it also refers to balanced realizations.
The class of functions
81
can simply be characterized in
terms of a balanced realization.
Proposition
5.3:
Let G
be
an
m
x
m
inner function with
minimal balanced realization
(A, B,
C,
D).
Then G
E
91
if
and only if rank
B
=
m,
or equivalently rank
C
=
m.
0
The following proposition shows that we can use a balanced
realization
of
G to construct a balanced realization for any
power of G.
Proposition
5.4:
Let G be a square inner transfer function
with minimal balanced realization
(A,
B,
C,
D)
having state
dimension
n
>
0.
Then for any
IC
>
1
the realization
and output balanced.
0
(Ak
>
Bk
i
ck,
Dk)
with
is a minimal balanced realization of
Gk
with state dimension
n
.
IC.
0
Examining the realization in the above proposition, reveals
a similar structure
of
observability and controllability matrices,
as has been discussed in the previous section; e.g., taking the
situation
IC
=
2,
it shows that the controllzbility matrix
of
(A2, B2)
contains the controllability matrix of
(A,
B)
as its
first block row.
-1.1
I
''