Balanciertes Abschneiden in beschränkten ... · MAX-PLANCK-INSTITUT DYNAMIK KOMPLEXER TECHNISCHER...

MAX−PLANCK−INSTITUT

DYNAMIK KOMPLEXER

TECHNISCHER SYSTEME

MAGDEBURG

9. Elgersburg Workshop,2. - 6. Marz 2014

Balanciertes Abschneiden in beschranktenFrequenzintervallen fur hochdimensionale

Systeme

Patrick Kurschnermit

Peter Benner und Jens Saak

Computational Methods in Systems and Control Theory (CSC)Max-Planck-Institut fur Dynamik komplexer Technischer Systeme

Max-Planck-Institut Magdeburg [email protected] P. Kurschner, Frequenzbegrenztes BT 1/20

[email protected]

Einleitung BT in beschrankten Intervallen Matrixlogarithmus Lyapunovgleichungen Numerische Beispiele

Einleitung

LTI System:

Σ :

x = Ax + Bu,

y = Cx .

Annahmen:

A ∈ Rn×n Hurwitz, n groß, dunnbesetzt:

Matrix-Vektorprodukte, sowie Losung von Au = v effizient moglich,Schur-, Eigen-, Singularwertzerlegung nicht.

B ∈ Rn×m mit m n,

C ∈ Rp×n mit p n.


[email protected]


Einleitung Σ : x = Ax + Bu, y = Cx

Balanciertes Abschneiden (BT)

Σ(A,B,C ), heißt balanciert, falls die Losungen P, Q deralgebraischen Lyapunovgleichungen (ALE)

AP+ PAT+ BBT = 0, ATQ +QA + CTC = 0,

P = Q = diag(σ1, . . . , σn) mit σ1 ≥ σ2 ≥ . . . ≥ σn > 0 erfullen.

σ1, . . . , σn sind die Hankel-Singularwerte von Σ.

Balancierte Realisierung durch Zustandsraumtransformation

T : (A,B,C) 7→ (TAT−1,TB,CT−1)

=

([A11 A12

A21 A22

],

[B1

B2

],[

C1 C2

]).

Abschneiden reduziertes Modell: (A, B, C ) = (A11,B1,C1).


[email protected]



Implementierung: SR (square root) Methode

1 Berechne (Cholesky-)Faktoren der Losungen derLyapunovgleichungen,

P = ZZT , Q = YY T .

2 Berechne Singularwertzerlegung (SVD)

Y TZ = [U1, U2 ]

[Σ1

Σ2

] [V T

1

V T2

].

3 DefiniereW := YV1Σ

−1/21 , V := ZU1Σ

−1/21 .

4 Das reduzierte Modell ist: (W TAV ,W TB,CV ).


[email protected]



Implementierung: SR (square root) Methode1 Berechne Niedrigrang-Faktoren der Losungen der

Lyapunovgleichungen,

P ≈ ZZT , Q ≈ YY T .


Y TZ = [U1, U2 ]

[Σ1

Σ2

] [V T

1

V T2

].

3 DefiniereW := YV1Σ

−1/21 , V := ZU1Σ

−1/21 .

4 Das reduzierte Modell ist: (W TAV ,W TB,CV ).


[email protected]


Einleitung

Integraldarstellung der Gramschen bzgl. Frequenzen ω

P =1

2π

∞∫−∞

(iωI − A)−1BBT (iωI − A)−Hdω,

Q =1

2π

∞∫−∞

(iωI − A)−HCTC (iωI − A)−1dω.

Haufig nur Frequenzintervall

Ω = [−ω2, − ω1] ∪ [ω1, ω2], 0 < ω1 < ω2 <∞

interessant / relevant.

Auch moglich: Betrachtung einer einzelnen Frequenz ω ∈ R.(P. Benners Vortrag 2013) [Du/Benner/Yang/Ye ’13]


[email protected]


Einleitung

Einschrankung auf Intervall Ω frequenzbegrenztes balanciertes Abschneiden (BTΩ)

Motivation / ZieleBTΩ erreicht verglichen mit standard BT

- eine hohere Genauigkeit in Ω bei gleicher reduzierterDimension, oder

- eine vergleichbare Genauigkeit in Ω bei niedrigerer reduzierterDimension.

Hauptanliegen dieses VortragesNumerisch durchfuhrbare und effiziente Umsetzung von BTΩ undobigen Zielstellungen.


[email protected]


BT in beschrankten Intervallen

Frequenzbegrenzte Gramsche bzgl. Ω = [−ω2, − ω1] ∪ [ω1, ω2]

PΩ :=1

2π

∫Ω


QΩ :=1

2π

∫Ω

(iωI − A)−HCTC (iωI − A)−1dω

Darstellung durch gewohnliche Gramsche: [Gawronski/Juang ’90]

PΩ := SΩP + PSTΩ ,

QΩ := STΩ Q + QSΩ

mit

SΩ = Re

1π

ω2∫ω1

(iωI − A)−1dω

.


[email protected]




PΩ :=1

2π

∫Ω


QΩ :=1

2π

∫Ω





mit

SΩ = 12π

∫Ω

(iωI − A)−1dω.


[email protected]




PΩ :=1

2π

∫Ω


QΩ :=1

2π

∫Ω





mit

SΩ = Re

1π

ω2∫ω1

(iωI − A)−1dω

.


[email protected]



Alternative Darstellung als Losung von Lyapunovgleichungen

APΩ + PΩAT + SΩBB

T + BBTSTΩ = 0,

ATQΩ + QΩA + STΩ CTC + CTCSΩ = 0

mit SΩ = Re

1

π

ω2∫ω1

(iωI − A)−1dω

= Re

(i

πlog((A + iω1I )

−1(A + iω2I )))

.

[Gawronski/Juang ’90, Petterson ’13]

Dabei ist logM der Hauptzweig des komplexen, matrixwertigen,naturlichen Logarithmus fur M ∈ Cn×n mit Λ(M) ∈ C\R−.z ∈ C\0: log z = log |z |+ i arg z


[email protected]



Implementierung von BTΩ

1 Berechne SΩ = Re(

iπ log

((A + iω1I )

−1(A + iω2I )))

.

2 Berechne Niedrigrang-Losungen PΩ ≈ ZΩZTΩ , QΩ ≈ YΩY

TΩ der

Lyapunovgleichungen


T + BBTSTΩ = 0,

ATQΩ + QΩA + STΩ CTC + CTCSΩ = 0.


Y TΩ ZΩ = [U1, U2 ]

[Σ1

Σ2

] [V T

1

V T2

].

4 Definiere WΩ := YΩV1Σ−1/21 , VΩ := ZΩU1Σ

−1/21 .

5 Das reduzierte Modell ist: (W TΩ AVΩ,W

TΩ B,CVΩ).


[email protected]


Matrixlogarithmus

Lyapunovgleichungen


T + BBTSTΩ = 0,

ATQΩ + QΩA + STΩ CTC + CTCSΩ = 0

mit

BΩ := SΩB, CΩ := CSΩ,

SΩ = Re(

iπ log

((A + iω1I )

−1(A + iω2I ))).

benotigen Matrixfunktion f (A) mit f (z) = log(

z+iω2

z+iω1

)stabile Berechnung z.B. mit “inverse scaling and squaring”(logm in MATLAB®), [Higham ’08]

benotigt u.A. Schurform von A,

Rechenaufwand O(n3), Speicherbedarf O(n2)⇒ fur große n zu teuer.


[email protected]


Matrixlogarithmus

Lyapunovgleichungen

APΩ + PΩAT + BΩB

T + BBTΩ = 0,

ATQΩ + QΩA + CTΩ C + CTCΩ = 0

mit BΩ := SΩB, CΩ := CSΩ,

SΩ = Re(

iπ log

((A + iω1I )

−1(A + iω2I ))).

benotigen “nur“ die Anwendung einer Matrixfunktion: u = f (A)v

erheblich attraktiver fur große n,

zahlreiche Methoden verfugbar,[Z.B.: Van der Vorst ’86, Higham ’08, Knizherman/Simoncini ’10, . . .]

B ∈ Rn×m, C ∈ Rp×n

⇒ m + p n solcher Anwendungen benotigt.


[email protected]


Matrixlogarithmus

1. Moglichkeit: Quadratur

Ausgangsdarstellung von SΩ:

BΩ = SΩB =1

πRe

ω2∫ω1

(iωI − A)−1Bdω

≈ 1

πRe

∑j

βj(iwj I − A)−1B

mit Quadraturgewichten, -punkten βj ∈ (0, 1), wj ∈ [ω1, ω2].

benotigt Losung von LGS (iwj I − A)−1B hohe Genauigkeit via adaptiver Quadratur ⇒ hohe Anzahl von LGS nicht fur alle Systeme effizientMATLAB:

integral(@(t) (1i*t*speye(n)-A)\B,w1,w2,’ArrayValued’,true));.


[email protected]


Matrixlogarithmus

2. Moglichkeit: Projektionsverfahren fur u = f (A)v

Sei v ∈ V = span V ⊂ Cn, dim(V) = k n, V HV = Ik .Approximation: u ≈ VV Hu = V u ∈ V, u := V Hu.Ritz-Galerkin Bedingung:

V u − f (A)VV Hv ⊥ V⇔ u = V H f (A)VV Hv

≈ f (V HAV )V Hv

⇒ u ≈ V f (V HAV )︸︷︷︸∈Ck×k

V Hv

Auch zweiseitige Projektion moglich:

V u − f (A)VW Hv ⊥ W ⊂ Cn

⇒ u ≈ Vf (W HAV )W Hv


[email protected]


Matrixlogarithmus

Verschiedene Moglichkeiten fur V:

Krylovraum

V = spanv ,Av ,A2v , . . . ,Ak−1v

=: Kk(A, v)

Shift-and-Invert Krylovraum

V = span

(A− τ I )−1v , (A− τ I )−2v , . . . , (A− τ I )−kv

Erweiterter Krylovraum: [Druskin/Knizherman ’98]

V = EKk := Kk(A, v) ∪ Kk(A−1, A−1v)


[email protected]


Matrixlogarithmus

Beispiele:

Testsystem Rechenzeiten

Herkunft n,m, p ω1 ω2 logm integral EK

fom,SLICOT1 1006,1,1 102 103 38 sec 1 sec 1.5 sec

fdm 2d,LyaPack2, Ex. 1

6400,4,41 102 3–4 h 2 min 5 sec

103 104 – 1.5 min 50 sec

chip,OWR3 20082,1,5 101 103 – 22 min 2.5 min

ifiss4,Ex. T-CD 2

66049,1,1 10 102 – 7.5 min 1.5 min

1SLICOT: http://slicot.org2Lyapack: http://www.netlib.org/lyapack/3Oberwolfach MOR Collection:

http://portal.uni-freiburg.de/imteksimulation/downloads/benchmark4http://www.cs.umd.edu/~elman/ifiss.html


http://slicot.org

http://www.netlib.org/lyapack/

http://portal.uni-freiburg.de/imteksimulation/downloads/benchmark

http://www.cs.umd.edu/~elman/ifiss.html

[email protected]


Lyapunovgleichungen

Gramsche der Erreichbarkeit

Suchen 0 ≤ P ∈ Rn×n s.d. AP + PAT + BBT = 0

Haben rank(BBT ) = m n P hat niedrigen numerischen Rang.[Penzl ’99, Antoulas/Sorenson/Zhou ’02,

Grasedyck ’04]

rank(P,u) = f n

100 200

100

10−10

u

σ63σ1

≈ u

Normierte Singularwerte

σ(P)

Vorgehen:

Suchen Niedrigrang-Approximation Pk = ZkZTk s.d.,

‖Pk − P‖ < τ 1 mit rank(Pk) = k n


[email protected]


Lyapunovgleichungen

Existenz von Niedrigranglosungen

Fur AX + XAT + W = 0, rank(W ) = m n, ε 1 existiertk = γ(Λ(A), ε)m, Xk = ZkZ

Tk , rank(Zk) = k mit

‖X − Xk‖ ≤ ε.

[Penzl ’00, Grasedyck ’04]

Algorithmen zur Berechnung von Niedrigranglosungen:

Krylovraumverfahren (EKSM, RKSM)[Druskin/Knizherman/Simoncini ’07/’11]

Niedrig-rang (alternating directions implicit) ADI Iteration (LR-ADI)[Penzl ’99, Li/White ’02, Benner/K./Saak ’13]

Ausschlaggebend fur Existenz von Niedrigranglosungen und damit derPerformance der Algorithmen sind m, Λ(A).


[email protected]


Lyapunovgleichungen

Frequenzbegrenzte Gramsche der Erreichbarkeit

Suchen 0 ≤ PΩ ∈ Rn×n s.d. APΩ + PΩAT + BΩB

T + BBTΩ = 0

Habenrank(BBT ) = m

< rank(BΩBT + BBT

Ω ) = 2m n.Erwartung:

rank(PΩ,u) > rank(P,u).Beobachtung:

rank(PΩ,u) < rank(P,u). 100 200

100

10−10

u

Normierte Singularwerte

σ(P)

σ(PΩ)

Anschauliche Erklarung: PΩ = SΩP + PSTΩ .

Spektralradius von SΩ

Sei 0 < ω1 < ω2 <∞, SΩ = Re(

iπ log

((A + iω1I )

−1(A + iω2I )))

.

Dann gilt ρ(SΩ) < 12 .


[email protected]


Lyapunovgleichungen

Frequenzbegrenzte Gramsche der Erreichbarkeit

Suchen 0 ≤ PΩ ∈ Rn×n s.d. APΩ + PΩAT + BΩB

T + BBTΩ = 0

Mit |ω2 − ω1| max |λ(A)| ist ρ(SΩ) 12 moglich

rank(PΩ,u) < rank(P,u),

PΩ theoretisch besser durch Niedrigranglosung approximierbar als P,

theoretisch besseres Konvergenzverhalten von Niedrigranglosern(vgl. mit AP + PAT + BBT = 0).

Aktuelle ForschungLyapunovgleichung AX + XAT + W = 0:Bisher geht nur rank(W ) in die theoretischen Betrachtungen bzgl.der Niedrigrangapproximation von X ein.Genauer Einfluss von W noch nicht ausreichend erforscht.


[email protected]


Numerische Beispiele

Beispiel 1: LyaPack, Example 1: [Penzl ’99]

FDM Diskretisierung des Operators

L(x) := ∆x − 100ξ1dx

dξ1

− 1000ξ2dx

dξ2

.

auf (0, 1)2 fur x = x(ξ1, ξ2), homogene Dirichlet RB.80 Gitterpunkte ⇒ n = 6400, m = p = 4.ω1 = 1, ω2 = 100.

Konvergenzhistorie fur LR-ADI fur PΩ, QΩ:

20 40 60 80 10010−10

10−4

102

tol.45.3 sec

49.5 sec

LR-ADI Iteration

Lya

p.

Res

idu

um P Q

PΩ QΩ


[email protected]




FDM Diskretisierung des Operators

L(x) := ∆x − 100ξ1dx

dξ1

− 1000ξ2dx

dξ2

.

auf (0, 1)2 fur x = x(ξ1, ξ2), homogene Dirichlet RB.80 Gitterpunkte ⇒ n = 6400, m = p = 4.ω1 = 1, ω2 = 100.Konvergenzhistorie fur LR-ADI fur PΩ, QΩ:

20 40 60 80 10010−10

10−4

102

tol.45.3 sec

49.5 sec

LR-ADI Iteration

Lya

p.

Res

idu

um P Q

PΩ QΩ


[email protected]




Reduktion auf Ordnung k = 5 mit/ohne Frequenzbegrenzung.Rechenzeiten ≈50 sec. fur BT & BTΩ.

10−1 101 103 104

100

101

ω

‖H(iω

)‖2

original BT BTΩ

10−1 101 103 10410−6

10−3

100

ω1 ω2

rela

tive

rF

ehle

r


[email protected]




Reduktion auf Ordnung k = 5 mit/ohne Frequenzbegrenzung.Rechenzeiten ≈50 sec fur BT, ≈100 sec BTΩ.

10−1 100 101 102

100

101

ω

‖H(s

)‖2

original BT BTΩ

10−1 100 101 10210−6

10−3

100

ω1 ω2

rela

tive

rF

ehle

r


[email protected]



Beispiel 2: IFISS, Ex. T-CD2:FEM Diskretisierung einer Konvektions-Diffusionsgleichung auf (−1, 1)2.256 Gitterpunkte ⇒ n = 66049, m = p = 1.ω1 = 10, ω2 = 100.


[email protected]



Beispiel 2: IFISS, Ex. T-CD2:Reduktion auf Ordnung k = 15 mit/ohne Frequenzbegrenzung.

10−1 100 103

50

100

ω1 ω2

ω

‖H(s

)‖2

original BT BTΩ

10−1 100 10310−5

10−2

101

101 102

ω

rela

tive

rF

ehle

r


[email protected]


Zusammenfassung / Ausblick

Zusammenfassung:Frequenzbegrenztes balanciertes Abschneiden auch furhochdimensionale Systeme umsetzbar:

Berechnung von f (A)B, f H(A)CT mit f (z) = log(

z+iω2z+iω1

), z.B., via

Quadratur oder Krylovraumverfahren moglich.Numerische Losungen der Lyapunovgleichungen durchNiedrigrangloser(PΩ,QΩ haben niedrigeren numerischen Rang als P,Q).Verglichen mit gewohnlichem BT: hohere Genauigkeit bei gleicherreduzierter Ordnung.

Weiterfuhrende Problemstellungen:BT in begrenzten Zeitintervallen (D) [Gawronski/Juang ’90]

begrenzte Zeit/Frequenzintervalle fur zeit-diskrete Systeme (D)[Horta/Juang/Longman ’93]

Frequenzgewichtetes BT (?) [Z.B. Enns ’84]

PW :=1

2π

∞∫−∞

(iωI − A)−1W (iω)BBTW (iω)H(iωI − A)−Hdω.


[email protected]

Balanciertes Abschneiden in beschränkten ... · MAX-PLANCK-INSTITUT DYNAMIK KOMPLEXER TECHNISCHER...

Documents

Transcript of Balanciertes Abschneiden in beschränkten ... · MAX-PLANCK-INSTITUT DYNAMIK KOMPLEXER TECHNISCHER...