Technische Universit¨at Berlin
Institut f¨ur Mathematik
Model reduction of descriptor systems
Tatjana Stykel
Technical Report 720-01
Preprint-Reihe des Instituts f¨ur Mathematik
Technische Universit¨at Berlin
Report 720-01 May 2015
Abstract
Model reduction is of fundamental importance in many control applications. We consider
model reduction methods for linear continuous-time descriptor systems. The methods
are based on balanced truncation techniques and closely related to the controllability and
observability Gramians and Hankel singular values of descriptor systems. The Gramians
can be computed by solving the generalized Lyapunov equations with special right-hand
sides. The numerical solution of generalized Lyapunov equations is also discussed. A nu-
merical example is given.
Key words: descriptor systems, controllability and observability Gramian, Hankel singular
values, model reduction.
Author address:
Tatjana Stykel
Institut f¨ur Mathematik, MA 4-5
Technische Universit¨at Berlin
Straße des 17. Juni 136
D-10623 Berlin
Germany
Supported by Deutsche Forschungsgemeinschaft, Research Grant ME 790/12-1.
1 Introduction
Consider a linear continuous-time system
E˙x(t) = A x(t) + B u(t), x(0) = x0,
y(t) = C x(t),(1)
where E,A∈Rn,n,B∈Rn,m,C∈Rp,n,x(t)∈Rnis the state vector, u(t)∈Rmis the
control input, y(t)∈Rpis the output and x0∈Rnis the initial value. The number of state
variables nis called the order of system (1). If I=E, then (1) is a standard state space
system. Otherwise, (1) is a descriptor system or generalized state space system. Such systems
arise naturally in many applications such as multibody dynamics [6, 14], electrical circuits
[8, 21], semidiscretization of partial differential equations [43] and they may have very large
order n.
We will assume throughout the paper that the pencil λE−Ais regular, i.e., det(λE−A)6= 0
for some λ∈C. The pencil λE −Ais called c-stable, if it is regular and all finite eigenvalues
of λE −Alie in the open left half-plane.
The model reduction problem consists in an approximation of the descriptor system (1)
by a reduced order system
e
E˙
ex(t) = e
Aex(t) + e
B u(t),ex(0) = ex0,
ey(t) = e
Cex(t),(2)
where e
E,e
A∈Rℓ,ℓ,e
B∈Rℓ,m,e
C∈Rp,ℓ and ℓ≪n. Note that systems (1) and (2) have
the same input u(t). We require for the approximate system (2) to preserve properties of
the original system (1) like regularity and stability. The descriptor system (1) consists of
differential equations that describe the dynamic behavior of the system as well as algebraic
equations characterizing a constraint manifold for the solution. Therefore, it is natural to
require for the reduced order system to have the same algebraic constraints as the original one.
Clearly, it is also desirable that the approximation error is small. Moreover, the computation
of the reduced order system should be numerically stable and efficient.
There exist various model reduction approaches for standard state space systems such as
balanced truncation [29, 31, 34, 35, 40, 41], moment matching approximation [16, 20], singular
perturbation approximation [27, 30] or optimal Hankel norm approximation [17]. Surveys on
system approximation and model reduction can be found in [2, 15]. One of most effective and
well studied model reduction techniques is balanced truncation which is closely related to the
two Lyapunov equations
AP+PAT=−BBT, ATQ+QA=−CTC.
The solutions Pand Qof these equations are called the controllability and observability
Gramians, respectively. The balanced truncation approach consists in transforming the state
space system to a balanced form whose the controllability and observability Gramians become
diagonal and equal together with a truncation of states that are both difficult to reach and
to observe, see [31] for details.
In this paper we generalize controllability and observability Gramians as well as Hankel
singular values for descriptor systems (Section 2). In Section 3 we present an extension of
balanced truncation algorithms [29, 40, 41] to descriptor systems. These algorithms are based
1
on splitting system (1) into its dynamic and algebraic parts and then reducing the order only
for the dynamic part via a standard model reduction method. Section 4 contains a numerical
example.
2 Descriptor systems
Consider the descriptor system (1). If the pencil λE −Ais regular, then it can be reduced to
the Weierstrass canonical form [37], that is, there exist nonsingular matrices Wand Tsuch
that
E=WInf0
0NT, A =WJ0
0In∞T, (3)
where Ikis the identity matrix of order kand Nis nilpotent with index of nilpotency ν. The
number νis the index of the pencil λE −A. Representation (3) defines the decomposition of
Rninto complementary deflating subspaces of dimensions nfand n∞corresponding to the
finite and infinite eigenvalues of the pencil λE −A, respectively. The matrices
Pr=T−1Inf0
0 0 Tand Pl=WInf0
0 0 W−1(4)
are the spectral projections onto the right and left deflating subspaces of the pencil λE −A
corresponding to the finite eigenvalues.
It is well known that if the pencil λE −Ais regular, u(t) is νtimes continuously differen-
tiable and x0is consistent, i.e., it belongs to the set of consistent initial conditions
X0=(x0∈Rn: (I−Pr)x0=
ν−1
X
k=0
F−k−1Bu(k)(0) ),
then the descriptor system (1) has a unique continuously differentiable solution x(t), see
[8, 11], that is given by
x(t) = F(t)Ex0+Zt
0
F(t−τ)Bu(τ)dτ +
ν−1
X
k=0
F−k−1Bu(k)(t).
Here
F(t) = T−1etJ 0
0 0 W−1(5)
is a fundamental solution matrix of the descriptor system (1), and the matrices Fkhave the
form
Fk=T−10 0
0−N−k−1W−1, k =−1,−2,....
Clearly, Fk= 0 for k < −ν.
If the initial condition x0is inconsistent or the input u(t) is not sufficiently smooth (for
example, in most control problems u(t) is only piecewise continuous), then the solution of the
descriptor system (1) may have impulsive modes [10, 11].
The rational matrix-valued function G(s) := C(sE −A)−1Bis called the transfer function
of the descriptor system (1). The transfer function G(s) is called proper if lim
s→∞ G(s)<∞. A
2
quadruple of matrices [ E, A, B, C ] is called a realization of G(s). We will also often denote
a realization of G(s) by sE −AB
C0.
Two realizations [ E, A, B, C ] and [ ˇ
E, ˇ
A, ˇ
B, ˇ
C] are restricted system equivalent if there
exist nonsingular matrices ˇ
Wand ˇ
Tsuch that
sˇ
E−ˇ
Aˇ
B
ˇ
C0=sˇ
WE ˇ
T−ˇ
W A ˇ
Tˇ
WB
Cˇ
T0.
A pair ( ˇ
W, ˇ
T) is called system equivalence transformation. A characteristic quantity of system
(1) is system invariant if it is preserved under a system equivalence transformation. The
transfer function G(s) is system invariant, since
G(s) = C(sE −A)−1B=ˇ
Cˇ
T−1ˇ
T(sˇ
E−ˇ
A)−1ˇ
Wˇ
W−1ˇ
B=ˇ
C(sˇ
E−ˇ
A)−1ˇ
B.
Other important results from the theory of rational functions and realization theory may be
found in [11, 23].
2.1 Controllability and observability
For descriptor systems there are various concepts of controllability and observability, e.g.,
[7, 11, 44].
Definition 2.1 System (1) and the triplet (E, A, B) are called controllable on the reachable
set (R-controllable) if
rank [ λE −A, B] = nfor all finite λ∈C.(6)
System (1) and the triplet (E, A, B) are called controllable at infinity (I-controllable) if
rank [ E, AKE, B] = n, where the columns of KEspan ker E. (7)
System (1) and the triplet (E, A, B) are called strongly controllable (S-controllable) if (6) and
(7) hold.
System (1) and the triplet (E, A, B) are called completely controllable (C-controllable) if (6)
holds and
rank [ E, B] = n. (8)
Observability is a dual property of controllability.
Definition 2.2 System (1) and the triplet (E, A, C) are called observable on the reachable
set (R-observable) if
rank λE −A
C=nfor all finite λ∈C.(9)
System (1) and the triplet (E, A, C) are called observable at infinity (I-observable) if
rank
E
KT
ETA
C
=n, where the columns of KETspan ker ET. (10)
3
System (1) and the triplet (E, A, C) are called strongly observable (S-observable) if (9) and
(10) hold.
System (1) and the triplet (E, A, C) are called completely observable (C-observable) if (9)
holds and
rank E
C=n. (11)
Clearly, conditions (7) and (10) are weaker than (8) and (11), respectively. Equiva-
lent algebraic characterizations of various concepts of controllability and observability for
descriptor systems are presented in [11, 44].
2.2 Controllability and observability Gramians
If we assume that the pencil λE −Ais c-stable, then the integrals
Gpc =Z∞
0
F(t)BBTFT(t)dt (12)
and
Gpo =Z∞
0
FT(t)CTCF(t)dt (13)
exist, where F(t) is as in (5). The matrix Gpc is called the proper controllability Gramian and
the matrix Gpo is called the proper observability Gramian of the continuous-time descriptor
system (1), see [4, 39]. The improper controllability Gramian of system (1) is defined by
Gic =−
−1
X
k=−ν
FkBBTFT
k,
and the improper observability Gramian of system (1) is defined by
Gio =−
−1
X
k=−ν
FT
kCTCFk.
In summary, the controllability Gramian of the descriptor system (1) is given by
Gc=Gpc +Gic
and the observability Gramian of the descriptor system (1) is given by
Go=Gpo +Gio.
If E=I, then Gpc =Gcand Gpo =Goare the usual controllability and observability Gramians
for the standard state space system [45].
The proper controllability and observability Gramians are the unique symmetric, positive
semidefinite solutions of the projected generalized continuous-time Lyapunov equations
EGpcAT+AGpcET=−PlBBTPT
l,
Gpc =PrGpc (14)
4
and ETGpoA+ATGpoE=−PT
rCTCPr,
Gpo =GpoPl,(15)
respectively, where Pland Prare given by (4), see [39]. If λE −Ais in Weierstrass canonical
form (3) and if the matrices
W−1B=B1
B2and CT −1= [ C1, C2]
are partitioned in blocks conformally Eand A, then we can show that
Gpc =T−1G1c0
0 0 T−T,Gpo =W−TG1o0
0 0 W−1,(16)
where G1cand G1osatisfy the standard continuous-time Lyapunov equations
JG1c+G1cJT=−B1BT
1,
JTG1o+G1oJ=−CT
1C1.
The improper controllability and observability Gramians are the unique symmetric, ne-
gative semidefinite solutions of the projected generalized discrete-time Lyapunov equations
AGicAT−EGicET=−(I−Pl)BBT(I−Pl)T,
PrGic = 0
and ATGioA−ETGioE=−(I−Pr)TCTC(I−Pr),
GioPl= 0,
respectively [39]. They can be represented as
Gic =T−10 0
0G2cT−T,Gio =W−T0 0
0G2oW−1,(17)
where G2cand G2osatisfy the standard discrete-time Lyapunov equations
G2c−NG2cNT=−B2BT
2,
G2o−NTG2oN=−CT
2C2.
Unfortunately, we do not know how to express the controllability and observability Gramians
Gcand Gofor the descriptor system (1) via the solutions of single Lyapunov equations.
The controllability and observability Gramians can be used to characterize controllability
and observability properties of system (1).
Theorem 2.3 [4, 39] Consider the descriptor system (1). Assume that λE −Ais c-stable.
1. System (1) is R-controllable if and only if the proper controllability Gramian Gpc is
positive definite on the subspace im PT
r.
2. System (1) is I-controllable if the improper controllability Gramian Gic is negative defi-
nite on the subspace ker PT
r.
5
3. System (1) is C-controllable if and only if the controllability Gramian Gcis positive
definite on im PT
rand negative definite on ker PT
r.
4. System (1) is R-observable if and only if the proper observability Gramian Gpo is positive
definite on the subspace im Pl.
5. System (1) is I-observable if the improper observability Gramian Gio is negative definite
on the subspace ker Pl.
6. System (1) is C-observable if and only if the observability Gramian Gois positive definite
on im Pland negative definite on ker Pl.
Note that the I-controllability (I-observability) of (1) does not imply that the improper
controllability (observability) Gramian is negative definite on ker PT
r(on ker Pl).
Example 2.4 The descriptor system (1) with
E=1 0
0 0 , A =2 0
0 1 , B =1
0, C = ( 1,0 )
is I-controllable and I-observable. We have Gic =Gio = 0 and PT
r=Pl, i.e., neither Gic nor
Gio are negative definite on ker PT
r= ker Pl.
Corollary 2.5 Consider the descriptor system (1). Assume that λE −Ais c-stable.
1. System (1) is R-controllable and R-observable if and only if
rank(Gpc) = rank(Gpo) = rank(GpcETGpoE) = nf.
2. System (1) is I-controllable and I-observable if
rank(Gic) = rank(Gio) = rank(GicATGioA) = n∞.
3. System (1) is C-controllable and C-observable if and only if
rank(Gc) = rank(Go) = rank(GpcETGpoE+GicATGioA) = n.
Proof. The result follows from Theorem 2.3 and representations (3), (16) and (17). ✷
2.3 Hankel singular values
The proper controllability and observability Gramians Gpc and Gpo as well as the improper
controllability and observability Gramians Gic and Gio are not system invariant. Indeed, under
a system equivalence transformation ( ˇ
W , ˇ
T) the proper and improper controllability Grami-
ans Gpc and Gic are transformed to ˇ
Gpc =ˇ
T−1Gpc ˇ
T−Tand ˇ
Gic =ˇ
T−1Gic ˇ
T−T, respectively,
whereas the proper and improper observability Gramians Gpo and Gio are transformed to
ˇ
Gpo =ˇ
W−TGpo ˇ
W−1and ˇ
Gio =ˇ
W−TGio ˇ
W−1, respectively. Then
ˇ
Gpc ˇ
ETˇ
Gpo ˇ
E=ˇ
T−1GpcETGpoEˇ
T ,
ˇ
Gic ˇ
ATˇ
Gio ˇ
A=ˇ
T−1GicATGioAˇ
T .
6
We see from these formulas that the spectra of the matrices GpcETGpoEand GicATGioAare
system invariant. These matrices play the same role for descriptor systems as the product
of the controllability and observability Gramians for standard state space systems [45]. We
have the following result.
Theorem 2.6 Let λE −Abe c-stable. Then the matrices GpcETGpoEand GicATGioAhave
real and non-negative eigenvalues.
Proof. It follows from (12) and (13) that Gpc and ETGpoEare symmetric and positive
semidefinite. In this case there exists a nonsingular matrix ˇ
Tsuch that
ˇ
T−1Gpc ˇ
T−T=
Σ10
Σ2
0
0 0
,ˇ
TTETGpoEˇ
T=
Σ10
0
Σ3
0 0
,
where Σ1, Σ2and Σ3are diagonal matrices with positive diagonal elements [45, p.76]. Then
we get
ˇ
T−1GpcETGpoEˇ
T=Σ2
10
0 0 .
Hence, GpcETGpoEis diagonalizable and it has real and non-negative eigenvalues.
Similarly, we can show that eigenvalues of GicATGioAare real and non-negative. ✷
Definition 2.7 Let nfand n∞be the dimensions of the deflating subspaces of the pencil
λE −Acorresponding to the finite and infinite eigenvalues, respectively. The square roots of
the largest nfeigenvalues of the matrix GpcETGpoEdenoted by ςjare called the proper Hankel
singular values of the c-stable continuous-time descriptor system (1). The square roots of the
largest n∞eigenvalues of the matrix GicATGioAdenoted by ϑjare called the improper Hankel
singular values of system (1).
We will assume that the proper and improper Hankel singular values are ordered decreas-
ingly, i.e.,
ς1≥ς2≥...≥ςnf≥0, ϑ1≥ϑ2≥...≥ϑn∞≥0.
The proper and improper Hankel singular values put together the set of the Hankel singular
values of the continuous-time descriptor system (1). For E=I, the proper Hankel singular
values are the classical Hankel singular values of the standard state space system [17].
Since the proper (improper) controllability and observability Gramians are symmetric and
positive (negative) semidefinite, there exist Cholesky factorizations
Gpc =RRT,Gpo =LTL(18)
and
Gic =−R′(R′)T,Gio =−(L′)TL′,(19)
where the matrices R,L,R′,L′∈Rn,n are the upper triangular Cholesky factors [28]. The
following lemma gives a connection between the proper and improper Hankel singular values
of system (1) and the standard singular values of the matrices LER and L′AR′.
7
Lemma 2.8 Assume that the descriptor system (1) is c-stable. Consider the Cholesky fac-
torizations (18) and (19) of the proper and improper Gramians of (1). Then the proper
Hankel singular values of system (1) are the nflargest singular values of the matrix LER,
and the improper Hankel singular values of system (1) are the n∞largest singular values of
the matrix L′AR′.
Proof. We have
ς2
j=λj(GpcETGpoE) = λj(RRTETLTLE) = λj(RTETLTLER) = σ2
j(LER),
ϑ2
j=λj(GicATGioA) = λj(R′(R′)TAT(L′)TL′A) = λj((R′)TAT(L′)TL′AR′) = σ2
j(L′AR′),
where λj(·) and σj(·) denote, respectively, the eigenvalues and the singular values ordered
decreasingly. ✷
3 Model reduction
In this section we consider the problem of reducing the order of the descriptor system (1).
3.1 Balanced reduction
For a given transfer function G(s), there are many different realizations [11]. Here we are
interesting only in particular realizations that are useful in applications.
Definition 3.1 A realization [ E, A, B, C ] of the transfer function G(s) is called R-minimal
if the triplet (E, A, B) is R-controllable and the triplet (E, A, C) is R-observable.
Definition 3.2 A realization [ E, A, B, C ] of the c-stable transfer function G(s) is called
proper balanced if the proper controllability and observability Gramians Gpc and Gpo are equal
and diagonal.
We will show that for a R-minimal realization [ E, A, B, C ] of the c-stable transfer function
G(s), there exists a system equivalence transformation (WT
b, Tb) such that the realization
[WT
bETb, W T
bATb, W T
bB, CTb] (20)
is proper balanced.
If (E, A, B) is R-controllable and (E, A, C) is R-observable, then by Corollary 2.5 we have
rank(Gpc) = rank(Gpo) = rank(GpcETGpoE) = nf. Consider the Cholesky factorizations (18)
of the proper controllability and observability Gramians. We may assume without loss of
generality that the Cholesky factors R,LT∈Rn,nfhave full column rank. It follows from
Lemma 2.8 and Corollary 2.5 that ςj=σj(LER)>0 for j= 1,...,nfand, hence, the matrix
LER ∈Rnf,nfis nonsingular.
Let
LER =UΣVT(21)
be a singular value decomposition of LER, where Uand Vare orthogonal matrices and
Σ = diag(ς1,...,ςnf) is nonsingular. Consider the matrices
Wb=hLTUΣ−1/2, W∞i, W ′
b=hERV Σ−1/2, W ′
∞i(22)
8
and
Tb=hRV Σ−1/2, T∞i, T ′
b=hETLTUΣ−1/2, T ′
∞i.(23)
Here the columns of matrices W∞and T∞span, respectively, the left and right deflating
subspaces of the pencil λE −Acorresponding to the infinite eigenvalues, and matrices W′
∞
and T′
∞satisfy WT
∞W′
∞= (T′
∞)TT∞=In∞. Clearly, for Prand Plas in (4), we have
I−Pr=T∞(T′
∞)Tand I−Pl=W′
∞WT
∞. Since
(I−Pr)RRT(I−Pr)T= (I−Pr)Gpc(I−Pr)T= 0,
(I−Pl)TLTL(I−Pl) = (I−Pl)TGpo(I−Pl) = 0,
we obtain that
RTT′
∞= 0 and LW ′
∞= 0.(24)
Then
(T′
b)TTb=Σ−1/2UTLERV Σ−1/2Σ−1/2UTLET∞
(T′
∞)TRV Σ−1/2(T′
∞)TT∞=In,
i.e., the matrices Tband T′
bare nonsingular and (T′
b)T=T−1
b. Similarly, we can show that
the matrices Wband W′
bare also nonsingular and (W′
b)T=W−1
b.
Using (18) and (21)-(24), we obtain that the proper controllability and observability
Gramians of the transformed system (20) have the form
T−1
bGpcT−T
b=Σ 0
0 0 =W−1
bGpoW−T
b,
where Σ = diag(ς1,...,ςnf) with the proper Hankel singular values ςj. Thus, (WT
b, Tb) with
Wband Tbas in (22) and (23), respectively, is the balancing transformation and realization
(20) is proper balanced.
Just as for standard state space systems [17, 31], the balancing transformation for de-
scriptor systems is not unique.
Remark 3.3 Note that the pencil λEb−Ab=WT
b(λE −A)Tbis in Weierstrass-like canonical
form. Indeed, from (21)-(23) we have
Eb=Σ−1/2UTLERV Σ−1/2Σ−1/2UTLET∞
WT
∞ERV Σ−1/2WT
∞ET∞=Inf0
0E∞,
Ab=Σ−1/2UTLARV Σ−1/2Σ−1/2UTLAT∞
WT
∞ARV Σ−1/2WT
∞AT∞=A10
0A∞,
where A1= Σ−1/2UTLARV Σ−1/2,E∞=WT
∞ET∞is nilpotent and A∞=WT
∞AT∞is
nonsingular. Clearly, the pencil λEb−Abis regular, c-stable and has the same index as
λE −A.
3.2 Balanced truncation
In the previous subsection we have considered a reduction of an R-minimal realization to
proper balanced form. However, computing the proper balanced realization may be ill-
conditioned as soon as Σ in (21) has small singular values. In addition, if the realization
9
is not R-minimal, then the matrix Σ is singular. In the similar situation for standard state
space systems one performs a model reduction by truncating the state components corre-
sponding to the zero and small Hankel singular values without significant changes of the
system properties, see, e.g., [31, 40]. This procedure is known as projection of dynamics or
balanced truncation. It can be also applied to the descriptor system (1).
The proper controllability and observability Gramians can be used to describe the future
output energy
Ey:= Z∞
0
yT(t)y(t)dt
and the minimal past proper input energy
Eu:= min
u∈Lm
2(R−)Z0
−∞
uT(t)u(t)dt (25)
that is needed to reach from x(−∞) = 0 the state x(0) = x0∈im Pr. Here R−= (−∞,0)
and Lm
2(R−) is the Hilbert space of all square integrable functions f:R−→ Rmsuch that
f(t) = 0 for t≥0.
Theorem 3.4 Consider a descriptor system (1) that is c-stable and R-minimal. Let Gpc and
Gpo be the proper controllability and observability Gramians of (1). If x0∈im Prand u(t) = 0
for t≥0, then
Ey= (x0)TETGpoEx0.
Moreover, for uopt(t) = BTF(−t)G−
pcx0, we have
Euopt = (x0)TG−
pcx0,
where G−
pc is the unique solution of
GpcG−
pcGpc =Gpc,
PT
rG−
pcPr=G−
pc.(26)
Proof. System (1) with x0∈im Prand u(t) = 0 for t≥0 has a unique solution given by
x(t) = F(t)Ex0. Then y(t) = CF(t)Ex0for t≥0 and, hence,
Ey=Z∞
0
yT(t)y(t)dt =Z∞
0
(x0)TETFT(t)CTCF(t)Ex0dt = (x0)TETGpoEx0.
Consider now the minimization problem (25) subject to the constraint for the initial
conditions
x0=Z0
−∞
F(−t)Bu(t)dt. (27)
Let µ∈Rnbe a Lagrange multiplier vector and let
L(u(t), µ) = Z0
−∞
uT(t)u(t)dt +µTx0−Z0
−∞
F(−t)Bu(t)dt
be the Lagrange function. For any variations ∆u(t) and ∆µwe have that
∆L(u(t), µ) = 2 Z0
−∞
uT(t)∆u(t)dt −µTZ0
−∞
F(−t)B∆u(t)dt
+ ∆µTx0−Z0
−∞
F(−t)Bu(t)dt= 0
10
if and only if (27) holds and
uT(t) = 1
2µTF(−t)B=1
2µTPrF(−t)B. (28)
Substitution of (28) in (27) gives
x0=1
2Z0
−∞
F(−t)BBTFT(−t)µ dt =1
2Z∞
0
F(t)BBTFT(t)µ dt =1
2Gpcµ. (29)
Using the representation for Gpc as in (16), where G1cis symmetric and positive definite
(since (E, A, B) is R-controllable), we obtain that (26) has a unique solution G−
pc given by
G−
pc =TTG−1
1c0
0 0 T. (30)
It follows from (29) that 2G−
pcx0=G−
pcGpcµ=PT
rµ. Hence, for the optimal input
uopt(t) = BTFT(−t)G−
pcx0,
we have that
Euopt =Z0
−∞
uT
opt(t)uopt(t)dt =Z0
−∞
(x0)T(G−
pc)TF(−t)BBTFT(−t)G−
pcx0dt
= (x0)T(G−
pc)TZ∞
0
F(t)BBTFT(t)dtG−
pcx0= (x0)TG−
pcx0.✷
Remark 3.5 Using (16) and (30) we obtain the relationships
GpcG−
pc =Pr,G−
pcGpc =PT
r,G−
pcGpcG−
pc =G−
pc
which, together with the first equation in (26), imply that G−
pc is, in general, a (1,2)-pseudo-
inverse of Gpc, see [9]. However, if PT
r=Pr, then G−
pc is the Moore-Penrose inverse [9] of Gpc.
Theorem 3.4 shows that a large input energy Euis required to reach from x(−∞) = 0 the
state x(0) = Prx0which lies in an invariant subspace of the proper controllability Gramian Gpc
corresponding to its small non-zero eigenvalues. Moreover, if x0is contained in an invariant
subspace of the matrix ETGpoEcorresponding to its small non-zero eigenvalues, then the
initial value x(0) = Prx0has a small effect on the output energy Ey. For the proper balanced
system, Gpc and ETGpoEare equal and, hence, they have the same invariant subspaces. In
this case the truncation of the states related to the small Hankel singular values does not
change system properties significantly.
Let [ E, A, B, C ] be a realization (not necessarily R-minimal) of the c-stable transfer func-
tion G(s). Consider the full rank Cholesky factorizations (18), where the matrices R∈Rn,rc,
LT∈Rn,rohave full column rank and rc= rank(Gpc)≤nf,ro= rank(Gpo)≤nf. Let
LER = [ U1, U0]Σ10
0 Σ0[V1, V0]T(31)
11
be an ”economy size” singular value decomposition of LER ∈Rro,rc, where [ U1, U0]∈Rro,r
and [ V1, V0]∈Rrc,r have orthogonal columns,
Σ1= diag(ς1,...,ςℓf) and Σ0= diag(ςℓf+1,...,ςr)
with ς1≥ς2≥... ≥ςℓf> ςℓf+1 ≥... ≥ςr>0 and r= rank(GpcETGpoE)≤min(rc, ro).
Then the reduced order realization can be computed as
"se
E−e
Ae
B
e
C0#=WT
ℓ(sE −A)TℓWT
ℓB
CTℓ0,(32)
where
Wℓ=hLTU1Σ−1/2
1, W∞i∈Rn,ℓ, Tℓ=hRV1Σ−1/2
1, T∞i∈Rn,ℓ (33)
and ℓ=ℓf+n∞. Here W∞and T∞form the bases of the left and right deflating subspaces,
respectively, corresponding to the infinite eigenvalues of λE −A.
Note that computing the reduced order descriptor system can be interpreted as performing
a system equivalence transformation ( ˇ
W , ˇ
T) such that
ˇ
W(sE −A)ˇ
Tˇ
WB
Cˇ
T0=
sEf−AfBf
sE∞−A∞B∞
CfC∞0
,
where the pencil λEf−Afhas the finite eigenvalues only, all eigenvalues of λE∞−A∞are
infinite, and then reducing the order of the subsystem [ Ef, Af, Bf, Cf] with nonsingular Ef.
Clearly, the reduced order system (32) is c-stable and R-minimal.
The described decoupling of system matrices is equivalent to the additive decomposition of
the transfer function as G(s) = Gp(s)+P(s), where Gp(s) = Cf(sEf−Af)−1Bfis the proper
part and P(s) = C∞(sE∞−A∞)−1B∞is the polynomial part of G(s). The transfer function
of the reduced system has the form e
G(s) = e
Gp(s) + P(s), where e
Gp(s) = e
Cf(se
Ef−e
Af)−1e
Bf
is the reduced subsystem. In this case the difference G(s)−e
G(s) = Gp(s)−e
Gp(s) is a
proper rational function, and we have the following upper bound for the H∞-norm of the
error system
kG(s)−e
G(s)kH∞:= sup
ω∈R
kG(iω)−e
G(iω)k ≤ 2(ςℓf+1 +...+ςnf) (34)
that has been derived in [17]. Here k · k denotes the spectral matrix norm.
3.3 Numerical aspects
To reduce the order of the descriptor system (1) we have to compute the full rank Cholesky
factors Rand Lof the proper controllability and observability Gramians that satisfy the
projected generalized Lyapunov equations (14) and (15). We need also the matrices W∞and
T∞whose columns span the left and right deflating subspaces, respectively, corresponding to
the infinite eigenvalues of λE −A.
The classical numerical methods for (generalized) Lyapunov equations are (generalized)
Bartels-Stewart and Hammarling methods [3, 22, 32] based on the preliminary reduction of
the matrix (matrix pencil) to the (generalized) Schur form [18], calculation of the solution of
12
a reduced system and back transformation. To solve numerically the projected generalized
Lyapunov equations (14) and (15) for the full rank Cholesky factors, we can use the genera-
lized Schur-Hammarling method proposed in [38]. Simultaneously, this method produces the
matrices W∞and T∞.
Algorithm 3.1 Solution of the projected generalized Lyapunov equations (14) and (15)
Input: System [E, A, B, C ]such that λE −Ais c-stable.
Output: Full rank Cholesky factors Rand Lof the proper controllability and observability
Gramians Gpc =RRTand Gpo =LTLand the matrices W∞and T∞that form the bases of
the left and right deflating subspaces corresponding to the infinite eigenvalues of λE −A.
Step 1. Use the GUPTRI algorithm [12, 13] to compute orthogonal transformation matrices
Uand Vsuch that
VTEU =EfEu
0E∞and VTAU =AfAu
0A∞,(35)
where Efis upper triangular, nonsingular and E∞is upper triangular with zeros on the
diagonal, Afis upper quasi-triangular and A∞is upper triangular, nonsingular.
Step 2. Use the generalized Schur method [25, 26] or the recursive blocked algorithm [24] to
compute the solution of the generalized Sylvester equation
EfY−ZE∞=−Eu,
AfY−ZA∞=−Au.
Step 3. Compute the matrices
VTB=B1
B∞, CU = [ Cf, C2].
Step 4. Use the generalized Hammarling method [22, 32] to compute the Cholesky factors Rf
and Lfof the solutions Xc=RfRT
fand Xo=LT
fLfof the generalized Lyapunov equations
EfXcAT
f+AfXcET
f=−(B1−ZB∞)(B1−ZB∞)T,
ET
fXoAf+AT
fXoEf=−CT
fCf.
Step 5. If rank(Rf)< nf, then compute the full column rank matrix R1from the QR
decomposition
RT
f=QRRT
1
0.
Otherwise, R1:= Rf.
Step 6. If rank(Lf)< nf, then compute the full row rank matrix L1from the QR decompo-
sition
Lf=QLL1
0.
Otherwise, L1:= Lf.
Step 7. Compute the full rank Cholesky factors
R=UR1
0, L = [ L1,−L1Z]VT.(36)
13
Step 8. Compute the matrices
W∞=V0
In∞and T∞=UY
In∞.(37)
The generalized Schur-Hammarling method costs O(n3) flops and can be used, unfortu-
nately, only for problems of small and medium size. Moreover, it does not take into account
the sparsity or any structure of the system and is not attractive for parallelization. Recently,
iterative methods related to the ADI method and the Smith method have been proposed to
compute the low rank approximation of the solutions of standard large-scale sparse Lyapunov
equations [1, 33]. It is important to extend these methods for projected generalized Lyapunov
equations. This topic is currently under investigations.
The following algorithm is a generalization of the square root balanced truncation method
[29, 40] for the descriptor system (1).
Algorithm 3.2 Generalized Square Root (GSR) method.
Input: A realization [E, A, B, C ]such that λE −Ais c-stable.
Output: A reduced order system [e
E, e
A, e
B, e
C].
Step 1. Use Algorithm 3.1 to compute the full rank factors Rand Lof the proper controlla-
bility and observability Gramians Gpc =RTRand Gpo =LLTas well as the bases W∞and T∞
of the left and right deflating subspaces of λE −Acorresponding to the infinite eigenvalues.
Step 2. Compute the ”economy size” singular value decomposition (31).
Step 3. Compute the matrices Wℓ= [ LTU1Σ−1/2
1, W∞]and Tℓ= [ RV1Σ−1/2
1, T∞].
Step 4. Compute the reduced order system [e
E, e
A, e
B, e
C] = [ WT
ℓETℓ, W T
ℓATℓ, W T
ℓB, CTℓ].
If the original system (1) is highly unbalanced, then the matrices Wℓand Tℓare ill-
conditioned. To avoid accuracy loss in the reduced system, a square root balancing free
method has been proposed for standard state space systems in [41]. This approach can be
generalized for descriptor systems as follows.
Algorithm 3.3 Generalized Square Root Balancing Free (GSRBF) method.
Input: A realization [E, A, B, C ]such that λE −Ais c-stable.
Output: A reduced order system [e
E, e
A, e
B, e
C].
Step 1. Use Algorithm 3.1 to compute the full rank factors Rand Lof the proper controlla-
bility and observability Gramians Gpc =RTRand Gpo =LLTas well as the bases W∞and T∞
of the left and right deflating subspaces of λE −Acorresponding to the infinite eigenvalues.
Step 2. Compute the ”economy size” singular value decomposition (31).
Step 3. Compute the ”economy size” QR decompositions
RV1=QcR0, LTU1=QoL0,
where Qc,Qo∈Rn,ℓfhave orthogonal columns and R0,L0∈Rℓf,ℓfare upper triangular,
nonsingular.
Step 4. Compute the reduced order system [e
E, e
A, e
B, e
C] = [ WT
ℓETℓ, W T
ℓATℓ, W T
ℓB, CTℓ],
where Wℓ= [ Qo, W∞]and Tℓ= [ Qc, T∞].
The GSR and GSRBF methods are mathematically equivalent in the sense that they de-
liver a reduced system with the same transfer function. But the matrices Wℓand Tℓcomputed
by the GSRBF method are often significantly better conditioned than those computed via
the GSR method.
14
Remark 3.6 In fact, we do not need to compute the full rank Cholesky factors Rand Lin
Step 7 and the matrices W∞and T∞in Step 8 of Algorithm 3.1. From (35) and (37) we have
WT
∞ET∞=E∞,WT
∞AT∞=A∞,WT
∞B=B∞and CT∞=CfY+C2=C∞. Moreover, it
follows from (35) and (36) that LER =L1EfR1. Thus, computation of the proper Hankel
singular values in Step 2 of Algorithms 3.2 and 3.3 can be performed working only with the
matrices L1,Efand R1. This reduces the computational cost and the memory requirement.
Note that the singular value decomposition of L1EfR1may be computed without forming
this product explicitly, see [19] for details.
4 Numerical example
Consider the holonomically constrained planar model of a truck [36]. The linearized equation
of motion has the form
˙
p(t) = v(t),
M˙
v(t) = Kp(t) + Dv(t)−GTλ(t) + B2u(t),
0 = Gp(t),
(38)
where p(t)∈R11 is the position vector, v(t)∈R11 is the velocity vector, λ(t)∈Ris the
Lagrange multiplier, Mis the positive definite mass matrix, Kis the stiffness matrix, Dis
the damping matrix, Gis the constraint matrix and B2is the input matrix. System (38)
together with the output equation y(t) = p(t) forms a descriptor system of order n= 23 with
m= 1 input and p= 11 outputs. The dimension of the deflating subspace corresponding to
the finite eigenvalues is nf= 20. This toy-example is presented to illustrate the reliability of
the proposed model reduction methods for descriptor systems.
All of the following results were obtained on an IBM RS 6000 44P Model 270 with relative
machine precision ǫ= 2.22 ×10−16 using the MATLAB mex-functions based on the GUPTRI
routine [12, 13] and the SLICOT library routines [5, 42].
Figure 1(a) shows the proper Hankel singular values ςj,j= 1,...,20. We approximate
system (38) by a model of order ℓ= 5. Note that the Bode plots of the original and reduced
systems are not presented, since they were impossible to distinguish. Figure 1(b) illustrate
how accurate the reduced order model approximate the original one. We display the amplitude
Bode plot of the error system G(iω)−e
G(iω) for a frequency range ω∈[1,103]. Comparison
of this error with the upper bound 2(ς3+...+ς20) = 1.69×10−5shows that the error estimate
(34) is tight.
5 Conclusion
We have generalized the controllability and observability Gramians as well as Hankel singular
values for descriptor systems and studied their important features. Model reduction methods
for descriptor systems have been presented. These methods are based on the balanced trun-
cation technique and deliver reduced order systems that preserve the regularity and stability
properties of the original system. Moreover, for these methods the approximation error bound
is available.
15
2 4 6 8 10 12 14 16 18 20
10−9
10−8
10−7
10−6
10−5
10−4
j
Proper Hankel singular values
ςj
100101102
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1x 10−5
Frequency (rad/sec)
Error system
kG(iω)−
e
G(iω)k
(a) (b)
Figure 1: Proper Hankel singular values of the linearized truck model (a) and Bode plot of
the error system (b)
References
[1] A.C. Antoulas, D.C. Sorensen, and S. Gugercin. A modified low-rank Smith method for
large-scale Lyapunov equations. Technical Report TR01-10, Department of Computa-
tional and Applied Mathematics, Rice University, 6100 Main St.-MS 134, Houston, TX
77005-1892, 2001.
[2] A.C. Antoulas, D.C. Sorensen, and S. Gugercin. A survey of model reduction methods for
large-scale systems. In V. Olshevsky, editor, Structured Matrices in Mathematics, Com-
puter Science and Engineering, Vol. I, Contemporary Mathematics Series, 280. American
Mathematical Society, 2001.
[3] R.H. Bartels and G.W. Stewart. Solution of the equation AX +XB =C.Comm. ACM,
15(9):820–826, 1972.
[4] D.J. Bender. Lyapunov-like equations and reachability/observability Gramians for de-
scriptor systems. IEEE Trans. Automat. Control, 32(4):343–348, 1987.
[5] P. Benner, V. Mehrmann, V. Sima, S. Van Huffel, and A. Varga. SLICOT - A subroutine
library in systems and control theory. Appl. Comput. Control Signals Circuits, 1:499–539,
1999.
[6] K.E. Brenan, S.L. Campbell, and L.R. Petzold. The Numerical Solution of Initial-Value
Problems in Differential-Algebraic Equations. Elsevier, North-Holland, New York, N.Y.,
1989.
[7] A. Bunse-Gerstner, R. Byers, V. Mehrmann, and N.K. Nichols. Feedback design for
regularizing descriptor systems. Linear Algebra Appl., 299:119–151, 1999.
[8] S.L. Campbell. Singular Systems of Differential Equation, I. Pitman, San Francisco,
1980.
16
[9] S.L. Campbell and C.D. Meyer. Generalized Inverses of Linear Transformations. Dover
Publications, New York, 1979.
[10] D. Cobb. Controllability, observability, and duality in singular systems. IEEE Trans.
Automat. Control, 29(12):1076–1082, 1984.
[11] L. Dai. Singular Control Systems. Lecture Notes in Control and Information Sciences,
118. Springer-Verlag, Berlin, Heidelberg, 1989.
[12] J.W. Demmel and B. K˚agstr¨om. The generalized Schur decomposition of an arbitrary
pencil A−λB: Robust software with error bounds and applications. Part I: Theory and
algorithms. ACM Trans. Math. Software, 19(2):160–174, 1993.
[13] J.W. Demmel and B. K˚agstr¨om. The generalized Schur decomposition of an arbitrary
pencil A−λB: Robust software with error bounds and applications. Part II: Software
and applications. ACM Trans. Math. Software, 19(2):175–201, 1993.
[14] E. Eich-Soellner and C. F¨uhrer. Numerical methods in multibody dynamics. B.G. Teub-
ner, Stuttgart, 1998.
[15] L. Fortuna, G. Nunnari, and A. Gallo. Model Order Reduction Techniques with Applica-
tions in Electrical Engineering. Springer-Verlag, London, 1992.
[16] K. Gallivan, E. Grimme, and P. Van Dooren. A rational Lanczos algorithm for model
reduction. Numerical Algorithms, 12(1-2):33–63, 1996.
[17] K. Glover. All optimal Hankel-norm approximations of linear multivariable systems and
their L∞-errors bounds. Internat. J. Control, 39(6):1115–1193, 1984.
[18] G.H. Golub and C.F. Van Loan. Matrix Computations. 3rd ed. The Johns Hopkins
University Press, Baltimore, London, 1996.
[19] G.H. Golub, K. Sølna, and P. Van Dooren. Computing the SVD of a general matrix
product/quotient. SIAM J. Matrix Anal. Appl., 22(1):1–19, 2000.
[20] E.J. Grimme, D.C. Sorensen, and P. Van Dooren. Model reduction of state space systems
via an implicitly restarted Lanczos method. Numerical Algorithms, 12(1-2):1–31, 1996.
[21] M. G¨unther and U. Feldmann. CAD-based electric-circuit modeling in industry. I. Math-
ematical structure and index of network equations. Surveys Math. Indust., 8(2):97–129,
1999.
[22] S.J. Hammarling. Numerical solution of the stable non-negative definite Lyapunov equa-
tion. IMA J. Numer. Anal., 2:303–323, 1982.
[23] V. Ionescu, C. Oar˘a, and M. Weiss. Generalized Riccati Theory and Robust Control: A
Popov Function Approach. John Wiley and Sons, Chichester, UK, 1999.
[24] I. Jonsson and B. K˚agstr¨om. Recursive bloked algorithms for solving triangular systems
– Part I: One-sided and coupled Sylvester-type matrix equations. ACM Trans. Math.
Software, 28(4):392–415, 2002.
17
[25] B. K˚agstr¨om and P. Poromaa. LAPACK-Style algorithms and software for solving the
generalized Sylvester equation and estimating the separation between regular matrix
pairs. ACM Trans. Math. Software, 22(1):78–103, 1996.
[26] B. K˚agstr¨om and L. Westin. Generalized Schur methods with condition estimators for
solving the generalized Sylvester equation. IEEE Trans. Automat. Control, 34:745–751,
1989.
[27] P.V. Kokotovi´c, R.E. O’Malley, and P. Sannuti. Singular perturbations and order reduc-
tion in control theory - an overview. Automatica, 12(2):123–132, 1976.
[28] P. Lancaster and M. Tismenetsky. The Theory of Matrices. Academic Press, Orlando,
FL, 2nd edition, 1985.
[29] A.J. Laub, M.T. Heath, C.C. Paige, and R.C. Ward. Computation of system balancing
transformations and other applications of simultaneous diagonalization algorithms. IEEE
Trans. Automat. Control, AC-32(2):115–122, 1987.
[30] Y. Liu and B.D.O. Anderson. Singular perturbation approximation of balanced systems.
Internat. J. Control, 50:1379–1405, 1989.
[31] B.C. Moore. Principal component analysis in linear systems: controllability, observabil-
ity, and model reduction. IEEE Trans. Automat. Control, AC-26(1):17–32, 1981.
[32] T. Penzl. Numerical solution of generalized Lyapunov equations. Adv. Comput. Math.,
8(1-2):33–48, 1998.
[33] T. Penzl. A cyclic low-rank Smith method for large sparse Lyapunov equations. SIAM
J. Sci. Comput., 21(4):1401–1418, 1999/00.
[34] T. Penzl. Algorithms for model reduction of large dynamical systems. Preprint
SFB393/99-40, Fakult¨at f¨ur Mathematik, Technische Universit¨at Chemnitz, D-09107
Chemnitz, Germany, December 1999. Available from http://www.tu-chemnitz.de/
sfb393/sfb99pr.html.
[35] M.G. Safonov and R.Y. Chiang. A Schur method for balanced-truncation model reduc-
tion. IEEE Trans. Automat. Control, AC-34(7):729–733, 1989.
[36] B. Simeon, F. Grupp, C F¨uhrer, and P. Rentrop. A nonlinear truck model and its
treatment as a multibody system. J. Comput. Appl. Math., 50:523–532, 1994.
[37] G.W. Stewart and J.-G. Sun. Matrix Perturbation Theory. Academic Press, New York,
1990.
[38] T. Stykel. Numerical solution and perturbation theory for generalized Lyapunov equa-
tions. Linear Algebra Appl., 349:155–185, 2002.
[39] T. Stykel. Generalized Lyapunov equations for descriptor systems: stability and in-
ertia theorems. Preprint SFB393/00-38, Fakult¨at f¨ur Mathematik, Technische Uni-
versit¨at Chemnitz, D-09107 Chemnitz, Germany, October 2000. Available from
http://www.tu-chemnitz.de/sfb393/sfb00pr.html.
18
[40] M.S. Tombs and I. Postlethweite. Truncated balanced realization of a stable non-minimal
state-space system. Internat. J. Control, 46(4):1319–1330, 1987.
[41] A. Varga. Efficient minimal realization procedure based on balancing. In A. EL Moudni,
P. Borne, and S.G. Tzafestas, editors, Proc. of IMACS/IFAC Symposium on Modelling
and Control of Technological Systems (Lille, France, May 7-10, 1991), volume 2, pages
42–47, 1991.
[42] A. Varga. Model Reduction Routines for SLICOT. NICONET Report 1999-8, 1999. Avail-
able from ftp://wgs.esat.kuleuven.ac.be/pub/WGS/REPORTS/ nic1999-8.ps.Z.
[43] J. Weickert. Applications of the Theory of Differential-Algebraic Equations to Partial Dif-
ferential Equations of Fluid Dynamics. Ph.D. thesis, Technische Universit¨at Chemnitz,
Chemnitz, 1997.
[44] E.L. Yip and R.F. Sincovec. Solvability, controllability and observability of continuous
descriptor systems. IEEE Trans. Automat. Control, AC-26(3):702–707, 1981.
[45] K. Zhou, J.C. Doyle, and K. Glover. Robust and Optimal Control. Prentice Hall, Upper
Saddle River, NJ, 1996.
19
