Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
SIAM J. CONTROL OPTIM.
c
2011 Society for Industrial and Applied Mathematics
Vol. 49, No. 2, pp. 686–711
LYAPUNOV EQUATIONS, ENERGY FUNCTIONALS, AND MODEL
ORDER REDUCTION OF BILINEAR AND STOCHASTIC SYSTEMS
∗
PETER BENNER
†
AND TOBIAS DAMM
‡
Abstract. We discuss the relation of a certain type of generalized Lyapunov equations to Grami-
ans of stochastic and bilinear systems together with the corresponding energy functionals. While
Gramians and energy functionals of stochastic linear systems show a strong correspondence to the
analogous objects for deterministic linear systems, the relation of Gramians and energy functionals
for bilinear systems is less obvious. We discuss results from the literature for the latter problem and
provide new characterizations of input and output energies of bilinear systems in terms of algebraic
Gramians satisfying generalized Lyapunov equations. In any of the considered cases, the definition
of algebraic Gramians allows us to compute balancing transformations and implies model reduction
methods analogous to balanced truncation for linear deterministic systems. We illustrate the per-
formance of these model reduction methods by showing numerical experiments for different bilinear
systems.
Key words. Lyapunov equations, Gramians, energy functionals, balanced truncation, model
order reduction, bilinear systems, stochastic systems
AMS subject classifications. 93A15, 93A30, 93C10, 93E20, 93B40
DOI. 10.1137/09075041X
1. Introduction. Model order reduction by balanced truncation is a standard
method, which has been introduced by Moore in [40] for linear deterministic control
systems of the form
(1.1) ˙x = Ax + Bu , y = Cx ,
where A ∈ R
n×n
, B ∈ R
n×m
, C ∈ R
p×n
and x(t) ∈ R
n
, y(t) ∈ R
p
, u(t) ∈ R
m
are the state, output, and input of the system, respectively. It preserves stability
and provides guaranteed error estimates. The main obstacle in its realization is the
computation of controllability and observability Gramians as solutions of the dual
Lyapunov equations
(1.2) AP + PA
T
= −BB
T
,A
T
Q + QA = −C
T
C.
Although this requires a higher effort than, e.g., methods, based on Krylov subspace
approximations, there are algorithms which allow balanced truncation for sparse sys-
tems of dimensions O(10
5
) and more; see, e.g., [46, 5, 36, 29].
The appealing features of balanced truncation have motivated similar approaches
for other system classes. In a series of papers, Scherpen and others (see, e.g., [50, 52,
51, 24, 27, 23, 59]) have developed a theory of balancing for nonlinear systems. The
notion of Gramians is replaced by controllability and observability energy functionals.
While on a conceptional base this generalization is quite attractive, often it is hardly
practicable from the computational point of view, since the energy functionals are
∗
Received by the editors February 23, 2009; accepted for publication (in revised form) December
28, 2010; published electronically April 5, 2011.
http://www.siam.org/journals/sicon/49-2/75041.html
†
Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magde-
burg, Germany (benner@mpi-magdeburg.mpg.de).
‡
AG Technomathematik, Fachbereich Mathematik, TU Kaiserslautern, 67663 Kaiserslautern,
Germany (damm@mathematik.uni-kl.de).
686
Downloaded 06/21/12 to 193.175.53.21. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
ENERGY FUNCTIONALS AND MODEL ORDER REDUCTION 687
obtained as solutions of nonlinear Hamilton–Jacobi equations, which are very expen-
sive for large dimensions. Recently in [33], there have been attempts to reduce the
complexity of the optimality equations by POD methods, but nevertheless the scope
of the approach seems to be limited. To overcome this drawback, other generaliza-
tions of Gramians have been considered, especially for bilinear systems in the context
of model order reduction; cf. [1, 25, 26, 63, 64, 11]. These generalized Gramians are
solutions of generalized Lyapunov equations of the forms
AP + PA
T
+
m
j=1
A
j
PA
T
j
= − BB
T
,
A
T
Q + QA +
m
j=1
A
T
j
QA
j
= − C
T
C,
(1.3)
where A, B, C are as in (1.1), (1.2), and A
j
∈ R
n×n
for j =1,...,m.IfA
j
=0for
all j, then the linear matrix equations in (1.3) boil down to (1.2). Therefore, we call
them generalized Lyapunov equations, but they should not be confused with other
types of generalized Lyapunov equations such as
AP E
T
+ EPA
T
= −BB
T
arising in the context of generalized state-space systems [45]. The Gramians defined
by (1.3) have already been considered in [49, 13] to characterize controllability and
observability of bilinear systems
˙x = Ax +
m
j=1
A
j
u
j
x + Bu , y = Cx .(1.4)
A first attempt to give an energy-based interpretation of these algebraic Gramians
apparently was made by Gray and Mesko in [25]. Their results look quite promising
and have been taken up recently, e.g., in [34, 17, 11, 28, 16]. Unfortunately, however,
the characterization of energy functionals given in [25] does not hold in the stated
generality. This issue, together with the nonuniqueness of singular value functions
and balancing of nonlinear system, has been addressed by Gray and Scherpen in [24].
Their analysis is quite subtle and applies to general nonlinear systems. However, it
does not discuss the special role of the algebraic Gramians of bilinear systems from
[49, 13], which is of particular interest to us from the computational point of view.
Moreover, the implications of [24] are not fully accounted for in subsequent papers on
bilinear systems, e.g., [11].
Hence, for the special case of bilinear systems, we try to clarify conditions un-
der which the algebraic Gramians give quantifiable information on reachability and
observability properties of the state vectors. In section 3, we first suggest a new ap-
proach to characterize unreachable and unobservable states (Theorem 3.1) via the
Gramians and then give a simple example to illustrate how an integrability condition
contradicts the characterization of energy functionals in [25]. We also discuss some
patches, which, however, do not give satisfactory error estimates for truncation errors.
Since we mainly aim at practical methods for model order reduction applicable
to large-scale problems, we review solvability conditions for the generalized Lyapunov
equations and, in section 4, provide numerical examples to support the significance
of the generalized Gramians in (1.3) for model order reduction of a bilinear system
(1.4)—at least in special cases.
Downloaded 06/21/12 to 193.175.53.21. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
688 PETER BENNER AND TOBIAS DAMM
On the other hand, it is a well-known fact that generalized Lyapunov equations of
the forms (1.3) are naturally associated to stochastic linear control systems; see, e.g.,
[32, 31, 15]. Therefore it is not surprising that P and Q can be interpreted as Gramians
of stochastic systems and that the method of balanced truncation can immediately
be carried over to this class of systems. Although, of course, work has been done in
this direction, e.g., in [42, 62, 59], to our knowledge this connection between bilinear
and stochastic Gramians has not really been documented in the literature so far; it is
thus another goal of this paper (pursued in section 2) to fill this gap and to open up
the field for further research.
2. Gramians and energy functionals of linear systems. The representation
of input and output energies for deterministic linear control systems as quadratic forms
involving the Gramians is a classical result. Factorizations of the Gramians are used to
compute balanced realizations which can be reduced by truncation. This method has
first been described for time-invariant systems in [40] and for time-varying systems in
[54, 60]. Our adaptation to stochastic systems is quite analogous. To clarify the idea
and the notation as well as for later reference in the discussion of bilinear systems,
we will briefly recapitulate some basic results for time-varying systems.
2.1. Time-varying deterministic linear equations. Let us consider a linear
control system
(2.1) ˙x = A(t)x + B(t)u, y= C(t)x
with coefficient matrices A(t) ∈ R
n×n
, B(t) ∈ R
n×m
,andC(t) ∈ R
p×n
being mea-
surable functions of t.Hereu ∈ R
m
and y ∈ R
p
are called input and output vectors,
while x ∈ R
n
is the state vector. For a given measurable input function u : R → R
m
and an initial vector x
0
∈ R
n
,letx(t, x
0
,u) denote the solution of (2.1) with input u
and x(0,x
0
,u)=x
0
; the corresponding output will be denoted by y(t, x
0
,u). For the
fundamental solution of the homogeneous system ˙x = A(t)x,wewriteΦ(t, τ).
Assuming that the homogeneous system ˙x = A(t)x is exponentially asymptoti-
cally stable, we can define the controllability and observability Gramian by
P =
0
−∞
Φ(0,τ)B(τ)B(τ)
T
Φ(0,τ)
T
dτ ,
Q =
∞
0
Φ(t, 0)
T
C(t)
T
C(t)Φ(t, 0) dt .
(2.2)
Furthermore, for x
0
∈ R
n
, we define the input and output energy functionals as
E
c
(x
0
) = inf
u∈L
2
]−∞,0]
x(−∞,x
0
,u)=0
0
−∞
u(t)
2
dt ,
E
o
(x
0
)=
∞
0
y(t, x
0
, 0)
2
dt .
Note that E
c
(x
0
)=∞ if x
0
cannot be reached from 0 over the time-interval ]−∞, 0].
It is easy to see that this is equivalent to x
0
∈ Im P . The following result is well
known. We present a proof both to motivate similar arguments for other systems and
to discuss some issues of forward and backward solutions (see Remark 2.2), which are
important for the stochastic case. Some details of the argument will also play a role
in the bilinear setup.
Downloaded 06/21/12 to 193.175.53.21. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
ENERGY FUNCTIONALS AND MODEL ORDER REDUCTION 689
Theorem 2.1. Consider the time-varying system (2.1) and the Gr amians P and
Q defined by (2.2).Ifx
0
∈ Im P ,then
E
c
(x
0
)=x
T
0
P
x
0
,
where P
denotes the Moore–Penrose inverse.
For x
0
∈ R
n
we have
E
o
(x
0
)=x
T
0
Qx
0
.
Pro o f.Forfixedx
0
we define u :]−∞, 0] → R
m
by
u(t)=B(t)
T
Φ(0,t)
T
P
x
0
.(2.3)
Then
x(t, x
0
,u)=
t
−∞
Φ(t, τ )B(τ)u(τ) dτ
is well defined by the exponential stability of the homogeneous equation and satisfies
(2.1) as well as the boundary conditions lim
t→∞
x(t, x
0
,u)=0and
x(0,x
0
,u)=
0
−∞
Φ(0,τ)B(τ)u(τ) dτ =
0
−∞
Φ(0,τ)B(τ)B(τ)
T
Φ(0,τ)
T
dτP
x
0
= PP
x
0
= x
0
.
Among all ˜u with x(∞,x
0
, ˜u) = 0 the given control has minimal L
2
-norm. To show
this, let us assume that ˜u = u +ˆu is another solution to the control problem. Then
x
0
=
0
−∞
Φ(0,τ)B(τ)(u(τ)+ˆu(τ)) dτ , whence
0
−∞
Φ(0,τ)B(τ)ˆu(τ) dτ =0.
This implies
0
−∞
u(t)
T
ˆu(t) dt =0,sothat
˜u
2
L
2
= u +ˆu
2
L
2
= u
2
L
2
+ ˆu
2
L
2
≥u
2
L
2
.
Since
u
2
L
2
=
0
−∞
u(t)
2
dt =
0
−∞
x
T
0
P
Φ(0,t)B(t)B(t)
T
Φ(0,t)
T
P
x
0
dt
= x
T
0
P
PP
x
0
= x
T
0
P
x
0
,
the proof of the first assertion is complete.
To prove the second, assume that the system starts in state x
0
and is not con-
trolled. Then the corresponding output is y(t)=C(t)Φ(t, 0)x
0
. The output energy is
the L
2
-norm of y,
E
o
(x
0
)=y
2
L
2
=
∞
0
y(t)
T
y(t) dt = x
T
0
∞
0
Φ(t, 0)
T
C(t)
T
C(t)Φ(t, 0) dt
x
0
= x
T
0
Qx
0
,
which we had to show.
Downloaded 06/21/12 to 193.175.53.21. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
690 PETER BENNER AND TOBIAS DAMM
Remark 2.2.
(i) If a state x
0
minimizes the quadratic form x
T
0
Px
0
, then either it is in ker P
or it maximizes x
T
0
P
x
0
among all x
0
∈ Im P . Hence a state is hard to reach
if x
T
0
Px
0
is small. Similarly, we can say that a state is hard to observe if
x
T
0
Qx
0
is small.
(ii) We will need later for the controllability Gramian to be interpreted as the
observability Gramian of the dual system. Note that Φ(0,t)=Φ(t, 0)
−1
,
whence
d
dt
Φ(0,t)=−Φ(0,t)A(t)Φ(t, 0)Φ(0,t)=−Φ(0,t)A(t) ,
d
dt
Φ(0, −t)
T
= A(−t)
T
Φ(0, −t)
T
;
see, e.g., [55]. Therefore,
P =
∞
0
Φ(0, −τ)B(−τ)B(−τ)
T
Φ(0, −τ)
T
dτ
=
∞
0
˜
Φ(τ,0)
T
B(−τ )B(−τ)
T
˜
Φ(τ,0) dτ ,
where
˜
Φ is the fundamental solution of the equation ˙x = A(−t)
T
x.
(iii) It is customary to define E
c
(x
0
) as the minimal energy needed to steer from
0tox
0
over the interval ] −∞, 0]. Alternatively, one can steer asymptotically
from t
0
to x
0
over an interval [t
0
,t
0
+T ], where t
0
∈ R and T>0 are arbitrary,
and set
E
(t
0
)
c
(x
0
) = inf
u∈L
2
[t
0
,t
0
+T ],T >0
x(t
0
+T,t
0
,u)=x
0
t
0
+T
t
0
u(t)
2
dt .
For time-varying systems in general this yields a different value, E
(t
0
)
c
(x
0
) =
E
c
(x
0
), but in the time-invariant case it is the same (which is well known and
follows also as a special case from our discussion in the next subsection). An
advantage is that we may also consider solutions for positive times, and these
are also defined for stochastic systems.
(iv) For completeness, let us recall that in the time-invariant case, P and Q satisfy
the Lyapunov equations AP + PA
T
= −BB
T
and QA + A
T
Q = −C
T
C.
2.2. Stochastic linear differential equations. Consider a stochastic linear
control system of Itˆo type (see, e.g., [3, 44])
dx = Ax dt +
N
j=1
A
j
xdw
j
+ Bu dt ,
y = Cx .
(2.4)
The w
j
= w
j
(t) are independent zero mean real Wiener processes on a probability
space (Ω, F,μ) with respect to an increasing family (F
t
)
t∈R
+
of σ-algebras F
t
⊂F.
Let L
2
w
(R
+
, R
q
) denote the corresponding space of nonanticipating stochastic pro-
cesses v with values in R
q
and norm
v(·)
2
L
2
w
:= E
∞
0
v(t)
2
dt
< ∞ ,
Downloaded 06/21/12 to 193.175.53.21. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php