Variation of the Balanced POD Algorithm for Model Reduction of Linear Systems* John R. Singler1 Abstract— We present a variation on an existing model reduction algorithm for linear systems based on balanced proper orthogonal decomposition (POD). In contrast to many computational approaches to balanced truncation, the algorithm variation approximates the reduced order model directly without first transforming the linear system. The algorithm is applicable to large-scale finite dimensional systems and a class of infinite dimensional systems. The algorithm variation is compared to the original balanced POD algorithm on example partial differential equation systems.

I. I NTRODUCTION Balanced truncation is a classical method for model reduction of linear systems of ordinary differential equations [1], [2]. The standard balanced truncation algorithms are only applicable for systems of moderate size, and a large amount of recent research has focused on the development and analysis of algorithms for balanced truncation (and the associated Lyapunov equations) of very high dimensional systems; see, e.g., [3], [4], [5], [6] and the many references therein. Balanced POD was originally introduced by Rowley in [7] as a balanced truncation algorithm for high dimensional linear systems of ordinary differential equations. Rowley’s balanced POD algorithm uses similar ideas to the method of snapshots for standard POD computations [8], however there are now two separate datasets. Balanced POD has now been widely used for linearized fluid flow model reduction and control problems; see, e.g., [9], [10], [11]. In [12], we extended Rowley’s balanced POD algorithm to an infinite dimensional setting that includes a class of partial differential equations (PDEs). In [13], we proved convergence of two data-based variations of the algorithm in [12]. The data-based algorithm variations were designed to be used when matrix approximations of the infinite dimensional operators are not available. This case can occur, e.g., when using an existing simulation code for the computations and access to the source code is difficult or restricted, or when spatial discretization of the PDE system yields matrix approximations of the operators that are not appropriate for existing matrix-based balanced truncation algorithms. In [13], we related the balanced POD modes (or balancing transformation) to the Hankel operator of the system, instead of the product of the system Gramians as in Rowley’s *This work was supported in part by the National Science Foundation under grant DMS-1217122. 1 J. Singler is with the Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO 65409, USA,

[email protected]

original paper. In this work, we combine this line of reasoning with existing results from infinite dimensional balanced truncation theory [14], [15] to produce another variation on Rowley’s balanced POD algorithm. The algorithm variation is transformation-free, i.e., the balanced reduced order model is approximated directly, as opposed to first transforming the system to balanced coordinates and then truncating as is done in the above works. II. T WO M ODEL P ROBLEMS To test the algorithm variation, we consider two example PDE systems from [13]. The first problem is a simple one dimensional hyperbolic PDE system that has a transfer function that can be evaluated exactly for comparison. The second problem is a two dimensional parabolic PDE system whose transfer function must be approximated. A. Model Problem 1 First, consider the following simple one dimensional hyperbolic PDE wt = −a(x)wx + b(x)u(t),

0 < x < 1,

t > 0,

with input u(t), boundary condition w(t, 0) = 0, and output Z y(t) =

1

c(x)w(t, x) dx. 0

We consider a(x) = β − αx with β > α > 0, and assume b(x) and c(x) are square integrable. This problem can be formulated as a linear system w(t) ˙ = Aw(t) + Bu(t),

y(t) = Cw(t) 2

(1)

as follows. Let X be the RHilbert space L (0, 1) with standard 1 inner product (f, g) = 0 f (x)g(x) dx. The operators B : 1 R → X and C : X → R1 are defined by [Bu](x) = b(x)u and Cw = (w, c). The operator A : D(A) ⊂ X → X is given by [Aw](x) = −a(x)wx (x), with domain D(A) = HL1 (0, 1). The space HL1 is the set of functions w ∈ X such that w(0) = 0, w is (weakly) differentiable, and wx ∈ X. We also need the Hilbert adjoint operator A∗ : D(A∗ ) ⊂ X → X, which is given by [A∗ z](x) = (a(x)z(x))x , with 1 1 D(A∗ ) = HR (0, 1). The space HR is defined similarly to 1 HL , except the boundary condition is now z(1) = 0. For a(x) given above, it can be checked that the exact transfer function G(s) = C(sI − A)−1 B is given by Z 1 Z x c(x) (β−αx)s/α b(v) (β−αv)−1−s/α dv dx. G(s) = 0

0

Later, we take β = 0.5, α = 0.4, b(x) = 1−x, and c(x) = x, and then these integrals can be computed exactly. B. Model Problem 2 Next, consider the parabolic PDE given by wt = µ(wxx + wyy ) + c1 (x, y)wx + c2 (x, y)wy + b(x, y)u(t), over the spatial domain Ω = [0, 1] × [0, 1], with Dirichlet boundary conditions on Γ0 (the bottom, right, and top walls):

A. Assumptions and Notation We consider the following general framework. Let X be a real separable Hilbert space with inner product (·, ·) and corresponding norm k · k = (·, ·)1/2 . Assume the operator A : D(A) ⊂ X → X generates an exponentially stable C0 semigroup eAt over X, and the operators B : Rm → X and C : X → Rm are both finite rank and bounded. The latter assumption implies that B and C take the form Bu =

m X

Cx = [ (x, c1 ), . . . , (x, cp ) ]T ,

bj uj ,

(2)

j=1

w(t, x, 0) = 0,

w(t, 1, y) = 0,

w(t, x, 1) = 0,

and a Neumann boundary condition on the left wall: wx (t, 0, y) = 0. The system output is Z η(t) = c(x, y)w(t, x, y) dx dy. Ω

where each bj and cj are in X and u = [ u1 , . . . , um ]T is a vector in Rm (see [16, Theorem 6.1]). Also, for a Hilbert space X, let L2 (0, ∞; X) be the Hilbert space of all functions x such that x(t) ∈ X for all t > 0 with finite norm Z ∞ 1/2 kx(t)k2 dt . kxkL2 (0,∞;X) = 0

B. The Hankel Operator and the Balanced Realization

We assume µ is a positive constant, the convection coefficients c1 (x, y) and c2 (x, y) are bounded, and the functions b(x, y) and c(x, y) are square integrable over Ω. This problem can also be written as a linear system of the form (1). Take the Hilbert space X to be L2 (Ω), the space of square integrable functions defined over Ω, with R standard inner product (f, g) = Ω f (x, y)g(x, y) dx dy. The operators B : R1 → X and C : X → R1 are defined similarly to the first model problem. For the linear operator A, we place the problem in weak form. Let V be the Hilbert space 1

V = {v ∈ H (Ω) : v = 0 on Γ0 }, with inner product (v, w)V = (vx , wx ) + (vy , wy ). Next, multiply the PDE by a test function v ∈ V and integrate by parts to yield (wt , v)X = a(w, v) + (b, v)X u,

The Hankel operator H : L2 (0, ∞; Rm ) → L2 (0, ∞; Rp ) of the linear system (A, B, C) is defined by Z ∞ [Hu](t) = [CBu](t) = CeA(t+s) Bu(s) ds, 0

where the controllability operator C : X → L2 (0, ∞; Rp ) and the observability operator B : L2 (0, ∞; Rm ) → X are defined by Z ∞ At [Cx](t) = Ce x, Bu = eAs Bu(s) ds. 0

In our earlier works [12], [13] we provided alternate forms for these operators. Proposition 1: Under the above assumptions, the operators C and B defined above are also given by [Cx](t) = [ (x, z1 (t)), . . . , (x, zp (t)) ]T , Z ∞X m Bu = uj (s) wj (s) ds. 0

where the bilinear form a : V × V → R is given by

(3) (4)

j=1

∗

a(w, v) = µ(w, v)V + (c1 wx , v)X + (c2 wy , v)X . The operator A : D(A) ⊂ X → X is (roughly) defined by (Aw, v) = −a(w, v),

for all w ∈ D(A) and v ∈ V .

See [13] for details. III. R EPRESENTATIONS OF THE H ANKEL O PERATOR AND THE BALANCED R EALIZATION We derive the transformation-free balanced POD algorithm variation by using an alternate representation of the Hankel operator of the system as well as different representations of the balanced realization.

where zi (t) = eA t ci and wj (t) = eAt bj are in L2 (0, ∞; X) and are the unique solutions of the linear evolution equations z˙i (t) = A∗ zi (t),

zi (0) = ci ,

(5)

w˙ j (t) = Awj (t),

wj (0) = bj ,

(6)

for i = 1, . . . , p and j = 1, . . . , m. Corollary 1: Under the above assumptions, the Hankel operator H : L2 (0, ∞; Rm ) → L2 (0, ∞; Rp ) of the system (A, B, C) is given by Z ∞ [Hu](t) = k(t, s) u(s) ds, (7) 0

where the p × m kernel function k(t, s) has ij entries
kij (t, s) = zi (t), wj (s) ,

∗

and zi (t) = eA t ci and wj (t) = eAt bj are the unique solutions of the linear evolution equations (5) and (6). For the class of systems (A, B, C) considered here, the Hankel operator is known to be trace class (or nuclear), and therefore compact [17, Theorem 4]. Therefore, there exist singular values σ1 ≥ σ2 ≥ · · · ≥ 0 (with repetitions according to multiplicity) and corresponding singular vectors {fk } ⊂ L2 (0, ∞; Rm ) and {gk } ⊂ L2 (0, ∞; Rp ) satisfying

Aϕj = lim+

∗

Hfk = σk gk ,

H gk = σk fk .

h→0

2

The singular vectors are orthonormal with respect to the L inner product, i.e., Z ∞ (fj , fk )L2 (0,∞;Rm ) = fjT (t) fk (t) dt = δjk , 0 Z ∞ (gj , gk )L2 (0,∞;Rp ) = gjT (t) gk (t) dt = δjk . 0

As in [13], we use the Hankel singular values and vectors to define the balancing transformation. If σj is nonzero, define the jth balancing modes ϕj and ψj in X by Z ∞X m −1/2 −1/2 fj,k (t) wk (t) dt, (8) ϕj = σj Bfj = σj 0 −1/2 ∗

ψj = σj

−1/2

k=1 p ∞X

Z

C gj = σj

0

gj,k (t) zk (t) dt,

(9)

k=1

where fj,k and gj,k are the kth components of fj and gj . As in the finite dimensional case, the balancing modes are the eigenvectors of the products of the Gramians, and the balancing modes provide the balancing transformation. Theorem 1 ([13]): Suppose the above assumptions hold, the Hankel singular values are distinct, and the Hilbert space is infinite dimensional. Then the balancing modes satisfy {ϕi } ⊂ D(A), {ψi } ⊂ D(A∗ ), and a balanced realization (Ab , B b , C b ) of the system (A, B, C) is given by Abij = (Aϕj , ψi ) = (ϕj , A∗ ψi ), b Bij = (bj , ψi ),

b Cij = (ϕj , ci ).

(10) (11)

The balanced POD algorithms in [13] are based on this representation of the balanced realization. Once the balancing modes are approximated, the entries of B b and C b are straightforward to construct. If a matrix approximation of A or A∗ is not available, then there are two options. First, for parabolic problems with the A operator derived from a bilinear form a : V × V → R (as in the second model problem), the entries of Ab can be computed as Abij = −a(ϕj , ψi )

(12)

if the bilinear form is available for computations and the computed balancing modes are in the Hilbert space V . Otherwise, the quantity Aϕj can be approximated using   −1/2 Aϕj = −σj B f˙j + Bfj (0) (13) m m Z ∞ X  X −1/2 = −σj wk (t) f˙j,k (t) dt + wk (0) fj,k (0) , 0

k=1

or A∗ ψi can be approximated using a similar formula. The derivatives of the Hankel singular vectors can be approximated by finite differences or other similar approaches, or using an indirect method proposed and analyzed in [13]. Remark 1: The approach of directly approximating Aϕj for a linearized fluid flow system was also taken in [10]. In that work, the authors used the approximation

k=1

eAδ ϕj − ϕj eAh ϕj − ϕj ≈ h δ

for small δ > 0. They approximated eAδ ϕj by simulating the linear PDE for 0 ≤ t ≤ δ with initial data ϕj . This approach will work, however it may require manual tuning of δ to ensure accuracy. Using expression (13) does not require any tuning, and it also requires less computational cost. In this work, we construct the balanced truncation directly using the Hankel singular values and singular vectors. This can be done using the following result. Theorem 2 ([14], [15]): Under the assumptions of Theorem 1, a balanced realization (Ab , B b , C b ) of the system (A, B, C) is given by  σ 1/2 Z ∞ j Abij = giT (t) g˙ j (t) dt, σi 0  σ 1/2 Z ∞ i f˙iT (t) fj (t) dt, (14) = σj 0 1/2

1/2

B b = [σ1 f1 (0), σ2 f2 (0), . . .]T , Cb =

1/2 1/2 [σ1 g1 (0), σ2 g2 (0), . . .],

(15) (16)

and Ab can also be expressed as 1
(σi σj ) 2 σj fiT (0)fj (0) − σi giT (0)gj (0) , i 6= j, 2 2 σi − σj (17) 1 1 (18) Abii = − fiT (0)fi (0) = − giT (0)gi (0). 2 2 The representation (14)-(16) was proved by Curtain and Glover in [14]. A similar formula to the alternate expression for Ab given in (17)-(18) was proved by Glover, Curtain, and Partington in [15, proof of Lemma 4.3] for the output normal realization; their argument is easily modified for the standard balanced realization to prove (17)-(18). This result has also been recently extended to cover a wider class of systems [18]. Below, we use a variation of the balanced POD algorithm to approximate both the Hankel singular values and singular vectors and then use Theorem 2 to approximate the balanced truncation. As far as the author is aware, Theorem 2 has not been used to approximate the balanced truncation of a PDE system except for [19]. (Theorem 2 has been used for theoretical investigations of balanced truncation; see, e.g., [15], [20], [21].) In [19], we used a continuous time eigensystem realization algorithm (ERA) to approximate the Hankel singular values and singular vectors. A numerical example showed the ERA successfully approximated the balanced truncated reduced order model, however balanced

Abij =

POD (using the representation in Theorem 1) was more accurate in certain cases. The goal of this work is to compare balanced POD computations using the different representations of the balanced realization using Theorems 1 and 2. Remark 2: Discrete time ERA [22] uses the singular value decomposition of a Hankel matrix and a representation of the (discrete time) balanced realization that has a similar form to that of Theorem 2. This method has been used to reduce continuous time PDE systems in [23], [24]. In fact, it is shown in [24] that discrete time ERA and discrete time balanced POD are equivalent for finite dimensional systems when all computations are exact. Of course, for complicated systems there will always be errors and the approaches will differ. Furthermore, continuous time balanced POD and ERA will yield different results from their discrete time counterparts—even if the continuous time nature of the system is taken into account in the discrete time algorithms. A thorough comparison of all of these approaches needs to be performed to see if and when one approach is preferable over the others; this work is one step toward this goal.

IV. BALANCED POD We summarize our approach to balanced POD from [13] for two collections of time varying functions {zi (t)} and {wj (t)} in L2 (0, ∞; X). Although these functions can be ar∗ bitrary, we consider the data {zi , wj } given by zi (t) = eA t ci At and wj (t) = e bj , i.e., the unique solutions of the linear evolution equations (5) and (6). In this case, the balanced POD of {zi , wj } consists of the Hankel singular values, Hankel singular vectors, and balancing modes for the system (A, B, C). These quantities can be used to form the balanced realization (Ab , B b , C b ) using either the representation in Theorem 1 derived from the balancing transformation, or the transformation-free representations in Theorem 2. As with continuous time POD, the balanced POD of two datasets can be approximated using many methods. In [13], we extended the method of snapshots and quadrature approaches for standard POD to balanced POD; here, we briefly outline the quadrature approach. The main idea is to approximate the Hankel operator with quadrature—this reduces the balanced POD computations to a matrix singular value decomposition. This is a different approach than was taken in earlier works [7], [12]; this balanced POD variation provides the same approximation of the Hankel singular values and balanced POD modes as in the earlier works, however we now also obtain approximations to the Hankel singular vectors at the quadrature points. For simplicity, we only consider the case of a single function in each dataset; the case of multiple functions is similar and can be treated by “stacking” the data as in the earlier works. w w Nw z Let {αiz , τiz }N i=1 and {αi , τi }i=1 be the weights and nodes of two quadrature rules. Apply the quadrature rules to the equations Hf = σg and H∗ g = σf and evaluate at

the quadrature nodes to obtain the approximate equations Nw X   z
Hf (τi ) ≈ αjw z(τiz ), w(τjw ) f (τjw ) ≈ σg(τiz ), j=1



 H∗ g (τiw ) ≈

Nz X


αjz w(τiw ), z(τjz ) g(τjz ) ≈ σf (τiw ).

j=1

Replace the above approximate equations by equalities and multiply the resulting equations by (αiz )1/2 and (αiw )1/2 , respectively. Define the scaled quantities azi = (αiz )1/2 z(τiz ), uj = (αjw )1/2 f (τjw ),

w 1/2 aw w(τjw ), j = (αj )

vi = (αiz )1/2 g(τiz ),

and let Γ be the Nz × Nw matrix with ij entries Γij = (azi , aw j ). This gives Γu = σv,

ΓT v = σu.

Therefore, the singular value decomposition of Γ gives approximations to the nonzero singular values of H and the corresponding singular vectors evaluated at the quadrature nodes. The balancing modes (8)-(9) are approximated using quadrature on the integrals: − 12

Nw X

− 12

Nz X

ϕk ≈ σk

− 21

αjw f (τjw ) w(τjw ) = σk

j=1

ψk ≈ σk

i=1

Nw X

uk,j aw j , (19)

j=1 − 21

αiz g(τiz ) z(τiz ) = σk

Nz X

vk,i azi .

(20)

i=1

Except for simple cases, we will not have the exact data {zi , wj }, but we will have approximate data {ziN , wjN }. In [13], we proved that if the approximate data converges to the exact data in L2 (0, ∞; X), then the singular values converge and the singular vectors and balancing modes corresponding to distinct singular values converge. V. T HE A LGORITHM VARIATION We now combine balanced POD and the representations of the balanced realization in Theorems 1 and 2 to approximate the balanced truncation (Ar , Br , Cr ), i.e., the rth order truncation of the balanced realization (Ab , B b , C b ). Balanced POD Algorithms for Balanced Truncation: 1) For i = 1, . . . , p, compute approximations ziN (t) and ∗ wjN (t) to the solutions zi (t) = eA t ci and wj (t) = eAt bj of the linear differential equation (5)-(6). 2) Compute approximations {σkN , fkN , gkN } of the Hankel singular values and singular vectors, and, if needed, N compute approximations {ϕN j , ψi } to the balancing modes, e.g., by the balanced POD quadrature approach presented above. 3) Choose r and approximate the balanced truncated model (Ar , Br , Cr ) by either (MA ) using equations (10)-(11) from Theorem 1, and computing Aϕj using (13), or computing A∗ ψi similarly;

(Ma ) using equations (10)-(11) from Theorem 1, and computing Ab using a bilinear form a : V ×V → R as in (12); (H0 ) using equations (17)-(18) and (15)-(16) from Theorem 2; (Hd ) using equations (14)-(16) from Theorem 2. Notes: A a • In step 3, (M ) and (M ) are the balanced POD algorithms using the balancing modes from [13]. The algorithm (Ma ) is only applicable for parabolic problems derived from a bilinear form (as in the second model problem), and one must be able to evaluate the bilinear form. 0 d • In step 3, (H ) and (H ) are the transformation-free variations of the balanced POD algorithm. The variation (H0 ) uses only approximations to the Hankel singular vectors evaluated at t = 0, while (Hd ) also uses approximations to integrals containing time derivatives of the Hankel singular vectors. 0 d • The variations (H ) and (H ) require less computation: specifically, the balancing modes are not computed, and the inner products in the realization (10)-(11) are not computed. • The balancing transformation can be ill-conditioned, even in the finite dimensional case [25]. Since (H0 ) and (Hd ) do not rely on a balancing transformation, they may avoid possible errors due to an ill-conditioned transformation. b d • The integral representation (14) of A in (H ) with i = j can be integrated exactly to give the representation (18) in (H0 ). Therefore, in our numerical experiments with the integral representation, we used the latter formula when i = j. VI. N UMERICAL R ESULTS We report numerical results for the two model problems outlined in Section II-B. We compare transfer functions of the model problems and the balanced POD reduced order models using the H∞ norm, which is the largest singular value of the function on the imaginary axis. Since we consider systems with only one input and one output, the H∞ norm of a transfer function G(s) is the maximum value of |G(iω)| for ω real. For the computations, we approximated all time integrals using the trapezoid rule. Also, we approximated the time derivatives of the singular vectors using two different methods: second order finite differences and an indirect approach from [13] (again with the trapezoid rule to approximate the time integrals). A. Model Problem 1 We begin with the first model problem, the 1D hyperbolic PDE system with transfer function that can be evaluated exactly. To approximate the solution of the PDEs (5) and (6), we used a simple first order accurate discontinuous Galerkin method (with piecewise constant basis functions) for the

spatial discretization and forward Euler for the time discretization. We used an equal spacing ∆x between the spatial nodes, and an equal time spacing ∆t = β∆x for forward Euler. We integrated the discrete equations until tmax + 1/2, where tmax = α−1 ln(β(β − α)) (see [13] for more details). The approximate H∞ norm between the exact transfer function G and the transfer function of the balanced POD reduced order model Gr is shown in Table I for various values of r and N , the number of equally spaced nodes in the discontinuous Galerkin computation. (For the H∞ norm computation, we chose ω in the interval 10−3 ≤ ω ≤ 102 .) We found nearly identical results for all approaches, including using the two approaches to approximating the derivatives of the singular vectors. We observed similar behavior for other cases. TABLE I A PPROXIMATE TRANSFER FUNCTION ERRORS kG − Gr k∞ FOR MODEL PROBLEM 1 WITH VARIOUS VALUES OF r AND N , THE NUMBER OF EQUALLY SPACED NODES . method

r = 3, N = 500

r = 10, N = 1000 1.93 × 10−3

2.00 × 10−2

7.64 × 10−3

1.93 × 10−3

10−2

10−3

1.93 × 10−3

2.00 ×

H0 Hd

2.00 ×

10−2

r = 5, N = 500 10−3

MA

7.64 × 7.64 ×

For this problem, the transformation free algorithm variation using the Hankel singular values and singular vectors and the representation of Theorem 2 works just as well as the approach using the balancing modes and the representation of Theorem 1. B. Model Problem 2 For the second model problem, we find that the different approaches can give quite different results. For our experiments, we take µ = 0.1 and convection coefficients c1 (x, y) = x sin(2πx) sin(πy), c2 (x, y) = y sin(πx) sin(2πy). We also take b(x, y) = 5 if x ≥ 1/2 and b(x, y) = 0 otherwise, and c(x, y) = 5 for all x, y. We use piecewise bilinear finite elements for the spatial discretization and the trapezoid rule for time stepping. We stopped time stepping when the L2 norm of the approximate solution was less than 10−3 . Since we do not have an exact transfer function for this problem, we compare the transfer function of the balanced truncation of the matrix approximating system (using only finite elements for the spatial discretization) and the transfer function of the balanced POD reduced order model. We used 31 equally spaced finite element nodes in each coordinate direction for all of the computations. The approximate H∞ norm errors are shown in Table II for r = 3 and various values of the time step ∆t. (For the H∞ norm computation, we chose ω in the interval 10−4 ≤ ω ≤ 104 .) Most of the algorithm variations give accurate

approximations with similar errors, except the approaches d (MA f ) and (Hf ) which approximate the time derivatives of the Hankel singular vectors using finite differences. The (MA f ) approach fails completely, even though the similar approach (MA i ), where the time derivatives of the Hankel singular vectors are approximated using the indirect method, does give an accurate approximation. It is unclear why the indirect methods are superior to the finite difference approximations—also see (Hdf ) and (Hdi )—to the time derivatives for this example. TABLE II A PPROXIMATE TRANSFER FUNCTION ERRORS FOR MODEL PROBLEM 2 WITH r = 3 AND VARIOUS VALUES OF ∆t. T HE TIME DERIVATIVES WERE APPROXIMATED BY EITHER SECOND ORDER FINITE DIFFERENCES

( SUBSCRIPT f ) OR THE INDIRECT METHOD ( SUBSCRIPT i) OF [13]. method

1 30 −2 10

r = 3, ∆t =

Ma

1.93 ×

MA f

7.48 × 10+2

MA i H0

4.93 ×

10−2

4.93 ×

10−2

Hdf Hdi

1 60 −3 10

r = 3, ∆t = 8.71 ×

9.07 × 100

1 100 −3 10

r = 3, ∆t = 2.88 ×

2.91 × 100

8.73 ×

10−3

2.76 × 10−3

8.68 ×

10−3

2.71 × 10−3

2.61 × 10−1

5.50 × 10−2

2.28 × 10−2

10−2

10−3

2.71 × 10−3

4.93 ×

8.68 ×

For r = 4, most of the errors behave similarly; see Table III. One difference is that the approach (Ma ) using the bilinear form is the most accurate for all values of the time step ∆t. However, as the time step ∆t is refined, all d of the other approaches except (MA f ) and (Hf ) are nearly as accurate. TABLE III A PPROXIMATE TRANSFER FUNCTION ERRORS kG∆t r − Gr k∞ FOR MODEL PROBLEM 2 WITH r = 4 AND VARIOUS VALUES OF ∆t. method

1 30 10−3

r = 4, ∆t =

1 60 10−3

r = 4, ∆t =

1 100 10−3

r = 4, ∆t =

Ma

7.86 ×

MA f

2.16 × 10+1

2.87 × 10+2

6.41 × 100

MA i

6.27 × 10−2

1.32 × 10−2

3.90 × 10−3

H0

6.27 ×

10−2

1.32 ×

10−2

3.85 × 10−3

8.58 ×

10−2

5.50 ×

10−1

4.98 × 10−2

1.32 × 10−2

3.85 × 10−3

Hdf Hdi

6.27 × 10−2

4.43 ×

2.47 ×

VII. C ONCLUSIONS We presented a variation on the balanced POD algorithm for balanced model reduction of a linear system. The variation computes the reduced order model directly without computing a balancing transformation or requiring access to (approximating) system matrices. This saves some computational work and avoids potential errors due to an ill-conditioned balancing transformation. The algorithm variation worked well on two example problems, although we did find that balanced POD with a transformation gave higher accuracy with coarser time steps for one example.

R EFERENCES [1] B. N. Datta, Numerical Methods for Linear Control Systems. San Diego, CA: Elsevier Academic Press, 2004. [2] K. Zhou, J. C. Doyle, and K. Glover, Robust and Optimal Control. Prentice-Hall, 1996. [3] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems. Philadelphia, PA: SIAM, 2005. [4] P. Benner, V. Mehrmann, and D. C. Sorensen, Eds., Dimension Reduction of Large-Scale Systems. Berlin: Springer-Verlag, 2005. [5] V. Druskin, L. Knizhnerman, and V. Simoncini, “Analysis of the rational Krylov subspace and ADI methods for solving the Lyapunov equation,” SIAM J. Numer. Anal., vol. 49, no. 5, pp. 1875–1898, 2011. [Online]. Available: http://dx.doi.org/10.1137/100813257 [6] V. Simoncini, “A new iterative method for solving large-scale Lyapunov matrix equations,” SIAM J. Sci. Comput., vol. 29, no. 3, pp. 1268–1288, 2007. [7] C. W. Rowley, “Model reduction for fluids, using balanced proper orthogonal decomposition,” Internat. J. Bifur. Chaos Appl. Sci. Engrg., vol. 15, no. 3, pp. 997–1013, 2005. [8] L. Sirovich, “Turbulence and the dynamics of coherent structures. I. Coherent structures,” Quart. Appl. Math., vol. 45, no. 3, pp. 561–571, 1987. [9] S. Ahuja and C. Rowley, “Feedback control of unstable steady states of flow past a flat plate using reduced-order estimators,” Journal of Fluid Mechanics, vol. 645, pp. 447–478, 2010. [10] S. Bagheri, L. Brandt, and D. Henningson, “Input–output analysis, model reduction and control of the flat-plate boundary layer,” Journal of Fluid Mechanics, vol. 620, pp. 263–298, 2009. [11] A. Barbagallo, D. Sipp, and P. Schmid, “Closed-loop control of an open cavity flow using reduced-order models,” Journal of Fluid Mechanics, vol. 641, pp. 1–50, 2009. [12] J. R. Singler and B. A. Batten, “A proper orthogonal decomposition approach to approximate balanced truncation of infinite dimensional linear systems,” Int. J. Comput. Math., vol. 86, no. 2, pp. 355–371, 2009. [13] J. R. Singler, “Balanced POD for model reduction of linear PDE systems: Convergence theory,” Numerische Mathematik, vol. 121, no. 1, pp. 127–164, 2012. [14] R. F. Curtain and K. Glover, “Balanced realisations for infinitedimensional systems,” in Operator Theory and Systems (Amsterdam, 1985). Basel: Birkh¨auser, 1986, pp. 87–104. [15] K. Glover, R. F. Curtain, and J. R. Partington, “Realization and approximation of linear infinite-dimensional systems with error bounds,” SIAM J. Control Optim., vol. 26, no. 4, pp. 863–898, 1988. [16] J. Weidmann, Linear Operators in Hilbert Spaces. New York: Springer-Verlag, 1980. [17] R. F. Curtain and A. J. Sasane, “Compactness and nuclearity of the Hankel operator and internal stability of infinite-dimensional state linear systems,” Internat. J. Control, vol. 74, no. 12, pp. 1260–1270, 2001. [18] C. Guiver and M. R. Opmeer, “Model reduction by balanced truncation for systems with nuclear Hankel operators,” submitted. [19] J. R. Singler, “Model reduction of linear PDE systems: a continuous time eigensystem realization algorithm,” in Proceedings of the American Control Conference, 2012, pp. 1424–1429. [20] S. M. Djouadi, “On the optimality of the proper orthogonal decomposition and balanced truncation,” in Proceedings of the 47th IEEE Conference on Decision and Control, 2008, pp. 4221 –4226. [21] ——, “On the connection between balanced proper orthogonal decomposition, balanced truncation, and metric complexity theory for infinite dimensional systems,” in Proceedings of the American Control Conference, 2010, pp. 4911–4916. [22] J.-N. Juang, Applied System Identification. Prentice Hall, 1994. [23] S. M. Djouadi, R. C. Camphouse, and J. H. Myatt, “Empirical reduced-order modeling for boundary feedback flow control,” Journal of Control Science and Engineering, vol. 2008, p. 11, 2008, article ID 154956. [24] Z. Ma, S. Ahuja, and C. Rowley, “Reduced-order models for control of fluids using the eigensystem realization algorithm,” Theoretical and Computational Fluid Dynamics, vol. 25, pp. 233–247, 2011. [25] M. G. Safonov and R. Y. Chiang, “A Schur method for balancedtruncation model reduction,” IEEE Trans. Automat. Control, vol. 34, no. 7, pp. 729–733, 1989.