Post on 18-Aug-2019
transcript
Numerische Losung von
(Optimal-)Steuerungsproblemen
fur PDEs:
Beating the Curse of
Dimensionality!?
Peter BennerInstitut fur Mathematik
TU Berlinund
benner@math.tu-berlin.de
http://www.math.tu-berlin.de/~benner
Mathematische Systemtheorie 200310.–12. Februar 2003, Elgersburg
Gliederung
Gliederung
• Steuerungsprobleme fur (parabolische) PDEs
• Linear-quadratische Optimalsteuerung bei PDEs
– Existenz und Struktur der Losung– Numerische Losung von algebraischen Riccati-
gleichungen: ein Quasi-Newtonverfahren– Direkte Feedback-Iteration– Numerische Beispiele
• Modellreduktion
– Balanciertes Abschneiden– Stochastisches Abschneiden– Numerische Beispiele
♦ Peter Benner ♦ Institut fur Mathematik ♦ 1
Steuerung von parabolischen PDEs
Steuerung von parabolischen PDEs
Parabolische partielle DGL in Gebiet Ω ⊂ Rd
(Z.B. Warmeleitungsgleichung, Diffusions-Konvektions-Gleichung)
xt −∇ · (a(ξ)∇x) + d(ξ)∇x+ r(ξ)x = B(ξ)u(t),
ξ ∈ Ω, t > 0
mit Anfangs- und Randbedingungen (∂Ω = Γ1 ∪ Γ2 ∪ Γ3)
x(ξ, t) = B1(ξ)u1(t), ξ ∈ Γ1,
xν(ξ, t) = B2(ξ)u2(t), ξ ∈ Γ2,
γx(ξ, t) + xν(ξ, t) = B3(ξ)u3(t), ξ ∈ Γ3,
x(ξ, 0) = x0(ξ), ξ ∈ Ω,
y = Cx, t ≥ 0.
• B = 0 =⇒ Randsteuerungsproblem
• Bj = 0 ∀j =⇒ Punktsteuerungsproblem
♦ Peter Benner ♦ Institut fur Mathematik ♦ 2
Systeme mit verteilten Parametern
Abstrakte Formulierung
Gegeben: separable Hilbertraume
X – der Zustandsraum,
U – der Steuerungsraum,
Y – der Beobachtungsraum,
und lineare Operatoren
A : dom(A) ⊂ X → X , B : U → X , C : X → Y.
Betrachte Systeme mit verteilten Parametern derForm
x = Ax+Bu, x(0) = x0 ∈ X ,
y = Cx,
also eine lineare Evolutionsgleichung mit Beobach-tungsgleichung.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 3
Systeme mit verteilten Parametern
Beispiel
Schwache Formulierung des Punktsteuerungspro-blems bei parabolischer PDE mit homogenenDirichlet-RBn à teste mit v ∈ H10(Ω)=⇒ abstraktes lineares System:
〈Aw, v〉 := −
∫
Ω
[a∇x · ∇v + (d · ∇x+ rx)v]dξ
Bu := B(ξ)u(t)
Cx := Cx,
mit
X = L2(Ω),
dom(A) = H2(Ω) ∩H10(Ω) ⊂ X ,
U = R
Y = R
♦ Peter Benner ♦ Institut fur Mathematik ♦ 4
Systeme mit verteilten Parametern Beispiele
Einige weitere Beispiele
• Schrodingergleichung
• hyperbolische Gleichungen 1. Ordnung
• Wellengleichung
• instationare Plattengleichung
• Kirchhoff’sche Gleichung
• Euler-Bernoulli Gleichung
[Lasiecka/Triggiani ’91, ’00]
♦ Peter Benner ♦ Institut fur Mathematik ♦ 5
Linear-quadratische Optimalsteuerung
Linear-quadratische Optimalsteuerung
Abstraktes LQR Problem:
Minimiere
J (x0, u) =1
2
∞∫
0
(‖y‖2Y + ‖u‖
2U
)dt,
uber u ∈ L2(0,∞;U) mit dynamischen Nebenbe-dingungen
x = Ax+Bu, x(0) = x0 ∈ X ,
y = Cx,
Abstrakter Rahmen fur
• optimale Steuerung (LQG/H2),
• Vorhersagemodelle, Schatzung und Filterung(Kalman-Filter),
• Risikobewertung,
• Parameteridentifizierung, inverse Probleme.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 6
Linear-quadratische Optimalsteuerung Existenz der Losung
Losung des LQR Problems
Satz [Gibson ’79,. . . ]Voraussetzungen:
• A infinitesimaler Erzeuger einer C0-Halbgruppe.
• B,C linear, beschrankt.
• (A,B) stabilisierbar, d.h. ∃F : X → U linear, beschrankt,
so daß die von A+BF erzeugteC0-Halbgruppe exponentiell
stabil ist.
• (C,A) entdeckbar, d.h. (A∗,C∗) stabilisierbar.
• ∀x0 ∈ X ∃ zulassige Steuerung u (J (x0, u) <∞).
Dann: Optimalsteuerung ist Feedback-Steuerung
u∞(t) = −B∗P∞x(t) = F∞x(t),
wobei P∞ eindeutige, selbstadjungierte Losung derOperator-Riccatigleichung
0 = R(P) := C∗C+A∗P+PA−PBB∗P
ist und folgende Eigenschaften hat:
• P∞ : dom(A)→ dom(A∗) linear, beschrankt,
• P∞ ≥ 0, d.h. P∞ ist positiv semidefinit.
• P∞ ist stabilisierend, d.h. die von A − BB∗P∞ erzeugte
C0-Halbgruppe ist exponentiell stabil.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 7
Linear-quadratische Optimalsteuerung Numerische Losung
Numerische Losung
Galerkin-Ansatz, Ortsdiskretisierung durch FEM
Minimiere
J (x0, u) =1
2
∞∫
0
(yTy + uTu
)dt,
uber u ∈ L2(0,∞;Rm), wobei
Mx = −Kx+Bu, x(0) = x0,
y = Cx,
mit Steifigkeitsmatrix K ∈ Rn×n, Massematrix M ∈Rn×n, B ∈ Rn×m, C ∈ Rp×n.
Optimalsteuerung ist wieder Feedback-Steuerung
u∗(t) = −BTP∗x(t) =: F∗x(t),
P∗ ≥ 0 ist stabilisierende Losung der algebraischenRiccatigleichung (ARG)
0 = R(P ) := CTC +ATP + PA− PBBTP,
mit A := −M−1K, B :=M−1B.
Konvergenz: Banks/Kunisch ’84, Lasiecka/Triggiani ’91 u.v.m.
Numerik: Banks/Ito ’91, Rosen/Wang ’95
♦ Peter Benner ♦ Institut fur Mathematik ♦ 8
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Numerische Losung von ARGs
Ein Zugang: Benutze Zusammenhang mit Hamilto-nischem Eigenwertproblem. [Potter ’66, Laub ’79,...]
Methoden:
• Losung uber (strukturierte, Block-) Schurform:
– QR Algorithmus [Laub ’79]
– SR Algorithmus [Bunse-Gerstner/Mehrmann ’86]
– Multishift-Algorithmus [Ammar/B./Mehrmann ’93]
– Einbettungsalgorithmus [B./Mehrmann/Xu ’97]
• Spektralprojektionsmethoden, z.B.,
– Signumfunktionsmethode
[Roberts ’79, Byers ’87, Gardiner/Laub ’86]
– Diskfunktionsmethode
[Maylshev ’93, Bai/Qian ’94, B./Byers ’95, B. ’97]
• Krylov-Unterraum-Verfahren ⇒ benutze sparseMatrix-Struktur, berechne approximative Losung:
– Block-Arnoldi-Verfahren [Jaimoukha/Kasenally ’94]
– Hamiltonisches Lanczos-Verfahren
[Ferng/Lin/Wang ’95, B./Faßbender ’95]
Approximative Losung i.d.R. nicht stabilisierend!
♦ Peter Benner ♦ Institut fur Mathematik ♦ 9
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
The Curse of Dimensionality
Probleme:
1. Arbeitsspeicher O(n2).
2. Rechenaufwand O(n3).
Beispiel:
Parabolische PDE, Ω = [0, 1]3, uniformes Gitter mith = 0.01 Ã n = 106, also
1. A benotigt 745 GBytes Speicher,
2. Schurzerlegung wurde ca. 20,000 Jahre dauern(auf PC mit Pentium 4, 2 GHz).
A sparse (dunnbesetzt) Ã Reduktion des Speicher-aufwands auf ca. 1 MB, aber sparse Matrix-Strukturwird wahrend der Rechnung zerstort!
♦ Peter Benner ♦ Institut fur Mathematik ♦ 10
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Das Newton-Verfahren fur ARGs
[Kleinman ’68, Mehrmann ’91, Lancaster/Rodman ’95,
B./Byers ’94/’98, B. ’97, Guo/Laub ’99 ]
Anderer Ansatz: Betrachte ARG
0 = R(P ) = CTC +ATP + PA− PBBTP
als nichtlineares Gleichungssystem.
Frechet-Ableitung von R(·) an der Stelle P :
R′
P : Z → (A−BBTP )TZ + Z(A−BBTP )
⇒ Modifiziertes Newton-Kantorovich-Verfahren:
Nj ← −(
R′
Pj
)−1
R(Pj),
Pj+1 ← Pj + tjNj,j = 0, 1, 2, . . .
Berechne Nj als Losung der Lyapunovgleichung
ATj Nj +NjAj = −R(Pj),
wobei Aj := A−BBTPj =: A+BFj.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 11
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Eigenschaften des mod. Newtonverfahrens[B. ’97, B./Byers ’98 ]
• Wahl von tj uber die Losung des zu R(P ) = 0gehorenden Minimierungsproblems
mintf(t) = min
t‖R(P + tN)‖
2F
= mintSpur
(
R(P + tN)TR(P + tN)
)
.
→ Exakte Liniensuche:
– ∃ lokales Minimum bei t∗j ∈ [ 0, 2 ],
– A(t) := A−BBT (Pj+tNj) stabil ∀t ∈ [ 0, 2 ]!
• Mit A stabil (Λ (A) ⊂ C−1) und P0 = 0 erhalteKonvergenzeigenschaften:
– Aj+1 = A(tj) stabil ∀ j.
– limj→∞ ‖R(Pj)‖F = 0, monoton.
– limj→∞ Pj = P∗ ≥ 0 (lokal quadratisch).
• Implementierung fur große Systeme: Benotige ef-fizienten Loser fur große (sparse) Lyapunovglei-chungen.
• ABER: P ∈ Rn×n =⇒ O(n2) Unbekannte!
♦ Peter Benner ♦ Institut fur Mathematik ♦ 12
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Niedrigrang-Approximation
Betrachte das Spektrum der Losung der ARG.
Beispiel:Lineare 1D Warmeleitungsgleichung mit Punktsteuerung,
Ω = [ 0, 1 ], FE-Diskretisierung mit linearen B-Splines,
h = 1/100 (=⇒ n = 101).
0 10 20 30 40 50 60 70 80 90 100 11010
−22
10−20
10−18
10−16
10−14
10−12
10−10
10−8
10−6
10−4
10−2
Spektrum der Lösung Ph für h=0.01
Index k
λ k
eps
Idee:
P = PT ≥ 0 =⇒ P =n∑
k=1
λkxkxTk = ZZT
λk ≈ 0, k > r =⇒ P ≈
r∑
k=1
λkxkxTk = Z(r)(Z(r))T .
♦ Peter Benner ♦ Institut fur Mathematik ♦ 13
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Iteration zur Approximation von Z(r)
[B./Li/Penzl ]
ATj Nj +NjAj = −R(Pj)
⇐⇒
ATj (Pj +Nj)︸ ︷︷ ︸=Pj+1
+(Pj +Nj)︸ ︷︷ ︸=Pj+1
Aj = −CTC − PjBB
TPj︸ ︷︷ ︸
=:−WjWTj
Setze Pj = ZjZTj fur Rang (Zj)¿ n:
ATj
(Zj+1Z
Tj+1
)+(Zj+1Z
Tj+1
)Aj = −WjW
Tj
⇓Modifiziere ADI-Iteration [Wachspress ‘88] zur Losungvon Lyapunovgleichungen: Berechne Zj+1 direkt;Losung der LGS mit Aj wegen m¿ n und
Aj = A+BFj = A − B · (BTZj) · ZTj ,
= sparse − m · ·
effizient mit Sherman-Morrison-Woodbury-Formel:
(A+ BFj)−1= (In − A
−1B(Im + FjA
−1B)
−1Fj)A
−1.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 14
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
ADI-Methode fur Lyapunovgleichungen
[Wachspress ‘88]
Sei A ∈ Rn×n stabil (Λ (A) ⊂ C−), W ∈ Rn×w
(w ¿ n), betrachte Lyapunovgleichung
AQ+QAT = −WW T .
ADI Iteration:
(A+ pjI)Q(j−1)/2 = −WW T −Qj−1(AT − pjI)
(A+ pjI)QjT = −WW T −Q(j−1)/2(A
T − pjI)
mit Parametern pj ∈ C− und pj+1 = pj falls pj 6∈ R.
Mit Q0 = 0 und geeigneter Wahl der pj:
limj→∞
Qj = Q superlinear.
• A selbstadjungiert: Explizite Berechnung der optimalen
Shift-Parameter moglich.
• A nicht-selbstadjungiert: Fur optimale Shift-Parameter lose
rationales min-max-Problem, benotige Λ (A) =⇒ heuri-
stische Wahl der Parameter, basierend auf conv(Λ (A)).
In der Praxis: Wende q Parameter zyklisch an.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 14.1
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Faktorisierte ADI-Iteration[B./Li/Penzl]
Setze Qj = YjYTj , Umformulierung =⇒
V1 ←√
−2Re (p1)(A+ p1I)−1W
Y1 ← V1
FOR j = 2, 3, . . .
Vj ←
√
Re (pj)√
Re (pj−1)
(
I − (pj + pj−1)(A+ pjI)−1)
Vj−1
Yj ←[Yj−1 Vj
]
⇓
Yjmax =[V1 . . . Vjmax
],
mit
Vj = ∈ Cn×w
undYjmaxY
Tjmax≈ Q
♦ Peter Benner ♦ Institut fur Mathematik ♦ 14.2
Linear-quadratische Optimalsteuerung Numerische Losung von ARGs
Newton-ADI fur ARGs
(A+BFj−1)TZjZ
Tj +ZjZ
Tj (A+BFj−1) = −Wj−1W
Tj−1
Faktorisierte Losung mit modifizierter ADI-Iteration
⇓Folge Y
(j)0 , Y
(j)1 , . . . , Y
(j)kmax
vonNiedrigrang-Approximationen an Zj.
Zj := Y(j)kmax
⇓Newton-Verfahren mit faktorisierten Iterierten
Pj = ZjZTj
⇓Faktorisierte Approximation an die Losung der ARG
P∗ ≈ ZjmaxZTjmax
Komplexitat O(nouter · ninner · nsolve), wobei
nouter = # außere Iterationen (Newtonschritte),
ninner = # (max/Durchschnitt) innere Iterationen (ADI-Schritte)
nsolve = Losen eines LGS mit Steifigkeitsmatrix
♦ Peter Benner ♦ Institut fur Mathematik ♦ 15
Linear-quadratische Optimalsteuerung Direkte Feedback-Iteration
Losung des LQR Problems
Erinnere: Losung des LQR Problems uber ARG:
Minimiere J (x0, u) =12
∫∞
0
(yTy + uTu
)dt, wobei
x = Ax+ Bu, x(0) = x0,
y = Cx.
Losung: Feedback-Steuerung
u∗(t) = −BTP∗x(t) =: F∗x(t),
wobei P∗ ≥ 0 stabilisierende Losung der ARG
0 = R(P ) := CTC + A
TP + PA− PBB
TP.
Aber: ARG ist ein Umweg, benotige nur Feedback-Matrix F∗ = −B
TP∗ = −BTZ∗Z
T∗
Besser: Direkte Berechnung der Feedback-Matrix!
• Newton-Iteration:
−Fj = BTZjZTj
j→∞−−−−→ − F∗ = BTZ∗Z
T∗
• Datiere Fj innerhalb der ADI-Iteration auf!
• Benotige nur Arbeitsspeicher fur Fj ∈ Rm×n.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 16
Linear-quadratische Optimalsteuerung Numerische Beispiele
Numerische Beispiele
Optimale Abkuhlung von Stahlprofilen
(Mannesmann/Demag bzw. [Troltzsch/Unger ’99, Penzl ’99, B./Saak ’02])
Mathematisches Modell: Randsteuerung fur linearisierte 2D
Warmeleitungsgleichung.
∂
∂tx =
λ
c · ρ∆x, ξ ∈ Ω
∂
∂nx =
1
λ(uk − x), ξ ∈ Γk, k = 1, . . . , 7,
∂
∂nx = 0, ξ ∈ Γ8.
=⇒ m = 7, p = 6
• FE-Diskretisierung mit ALBERT, li-
neare Elemente.
• Anfangstriangulierung (n = 106)→
• (Adaptive) Verfeinerung durch Bisek-
tion.
• 1/2/3 Schritte globale Gitterverfeine-
rung =⇒ n = 371/1357/5177.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 17
Linear-quadratische Optimalsteuerung Numerische Beispiele
Drei alternative Vorgehensweisen:
1. Lose LQR Problem fur ODE nach Ortsdiskre-tisierung (global verfeinerter Triangulierung) Ãvertikale Linienmethode.
⇒ Feedback-Matrix wird am Anfang berechnetund bleibt im gesamten Zeitverlauf konstant.
2. Adaptive Regelung: Ortsdiskretisierung wie bei 1.,verwende temperaturabhangige Materialgesetze(Ã nichtlin. Warmeleitungsgleichung), Feedback-Matrix wird nach jedem k. Zeitschritt neu berech-net (k = 5, 10, . . .).
3. Adaptive Regelung und adaptives Gitter: AdaptiveGitterverfeinerung/-vergroberung in jedem Zeit-schritt à horizontale Linienmethode, Neuberech-nung der Feedback-Matrix in jedem Zeitschritt.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 18
Linear-quadratische Optimalsteuerung Numerische Beispiele
Vertikale Linienmethode
1.000
0.500
1.000
0.500
1.000
0.500
♦ Peter Benner ♦ Institut fur Mathematik ♦ 19
Linear-quadratische Optimalsteuerung Numerische Beispiele
Adaptive Regelung (k = 10)
0.0000
47540.06
0.500
1.000 0.0000
47540.06
0.500
1.000 0.0000
47540.06
0.500
1.000
0.0000
47540.06
0.500
1.000 0.0000
47540.06
0.500
1.000 0.0000
47540.06
0.500
1.000
0.0000
47540.06
0.500
1.000 0.0000
47540.06
0.500
1.000 0.0000
47540.06
0.500
1.000
♦ Peter Benner ♦ Institut fur Mathematik ♦ 20
Linear-quadratische Optimalsteuerung Numerische Beispiele
Horizontale Linienmethode (k = 10)
0.0000
45009.65
0.500
1.000
0.0000
45009.65
0.500
1.000
0.0000
45009.65
0.500
1.000
0.0000
45009.65
0.500
1.000
0.0000
45009.65
0.500
1.000
0.0000
45009.65
0.500
1.000
0.0000
45009.65
0.500
1.000
45009.65
0.0000
0.500
1.000
45009.65
0.0000
0.500
1.000
♦ Peter Benner ♦ Institut fur Mathematik ♦ 21
Linear-quadratische Optimalsteuerung Fazit/Ausblick
Fazit/Ausblick
• Newton-ADI-Verfahren ist effiziente und zu-verlassige Methode zur Losung großer sparserARGs/LQR Probleme.
• Newton-Verfahren garantiert die Stabilisierungs-eigenschaft der Niedrigrang-Losungen der ARGbzw. der Feedback-Matrix!
• Direkte Berechnung der Feedback-Matrix fur LQRProbleme ist moglich ohne den ARG Umweg.
• Konvergenztheorie =⇒ berechnete Optimallosun-gen des diskreten Problems konvergieren gegenOptimum fur PDE Problem!
• Viele Detailsverbesserungen notig fur effizienteImplementierung!
• Vielfaltige Anwendungen bei (adaptiven)Steuerungs- und Regelungsproblemen fur PDEs.
• Ahnlich furH∞-Optimierung von PDE Problemenmoglich.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 22
Modellreduktion Motivation
Modellreduktion
Ziel: Approximiere System durch System niedrigererZustandsraumdimension.
Wozu?
• Berechnung der optimalen Steuerung bei hoherKomplexitat (3D Probleme) nicht mehr moglich.
• Probleme mit NB: Losung mit QP/NLP notig.
• Echtzeitregelung nur bei niedriger Komplexitatdes Reglers moglich.
– Feedback-Regler: Lineares
System der Ordnung N .
– Moderner Reglerentwurf
(H2-/H∞): N ≥ n.
G
?
-
˙x = Fx+ Ey
u = Hx+Ky
u y
• Wiederholte Simulation mit demselben Modell furviele verschiedene Eingangssignale:
– VLSI Chip Design; Verifikation von Layouts von immer
komplexeren Schaltkreisen in immer kurzer werdenden
Entwicklungszyklen.
– Simulation von gekoppelten PDE Systemen, z.B. bei
Stromungsbeeinflußung.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 23
Modellreduktion Balanciertes Abschneiden
Modellreduktion basierend auf Betrachtung derabsoluten Fehler
Gesucht: reduziertes Modell
˙x(t) = Ax(t) + Bu(t), u(t) ∈ Rm,
y(t) = Cx(t) + Du(t), y(t) ∈ Rp
der Ordnung ` ¿ n mit “kleinem” absoluten Fehler‖G− G‖∞, wobei
‖G‖∞ = ess supω∈R
σmax(G(ıω)).
y(s) = G(s)u(s), y(s) = G(s)u(s) =⇒
‖G− G‖∞ = supu∈L2
‖(G− G)u‖2‖u‖2
= supu∈L2
‖y − y‖2‖u‖2
=⇒ Entferne die Zustandsraumkomponenten, die amwenigsten in den Energietransfer von u nach y invol-viert sind, d.h. die zu den kleinsten HSV gehorendenKomponenten.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 24
Modellreduktion Balanciertes Abschneiden
Balanciertes Abschneiden I
• Steuerbarkeits-, Beobachtbarkeits-Gram’sche Ma-trizen von G losen die Lyapunovgleichungen
AP + PAT +BBT = 0,
ATQ+QA+ CTC = 0.
• (A,B,C,D) ist eine balancierte Realisierung desSystems G, falls
P = Q =
σ1. . .
σn
.
• σ1 ≥ σ2 ≥ . . . ≥ σn > 0 sind die Hankel-Singularwerte (HSV) des Systems; HSV sind Sy-steminvarianten.
• (A,B,C,D) minimal =⇒ Zustandsraumtransfor-mation auf balancierte Realisierung existiert.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 25
Modellreduktion Balanciertes Abschneiden
Balanciertes Abschneiden II
• Zustandsraumtransformation mit T ∈ Rn×n:
TAT−1 =
[A11 A12A21 A22
]
, A11 ∈ R`×`
TB =
[B1B2
]
, CT−1 =[C1 C2
].
Partitioniere
T =
[Tl
Wl
]
, Tl ∈ R`×n
T−1 =[Tr Wr
], Tr ∈ Rn×`.
• Modellreduktion durch Abschneiden:
A = A11 = TlATr
B = B1 = TlB,
C = C1 = CTr.
D = D.
Beachte: limω→∞(G(ıω)− G(ıω)) = D−D = 0.
• Wahl von T , ` so, daß ‖G− G‖∞ “klein”!
♦ Peter Benner ♦ Institut fur Mathematik ♦ 26
Modellreduktion Balanciertes Abschneiden
Balanciertes Abschneiden III
Fur balancierende Transformation, d.h.
G(s) ≡
[TAT−1 TB
CT−1 D
]
mit P = Q =:
[
Σ
Σ2
]
ist das reduzierte Modell
G(s) ≡
[A B
C D
]
≡
[A11 B1C1 D
]
balanciert, minimal, stabil mit Gram’schen Matrizen
P = Q = Σ =
[σ1
. . .
σ`
]
.
=⇒ Berechenbare Fehlerabschatzung
‖G− G‖∞ ≤ 2
n∑
k=`+1
σk.
=⇒ Adaptive Wahl von `!
♦ Peter Benner ♦ Institut fur Mathematik ♦ 27
Modellreduktion Balanciertes Abschneiden
Balanciertes Abschneiden: SR Methode
[Heath/Laub/Paige/Ward ’87, Tombs/Postlethwaite ’87 ]
Gram’sche Matrizen sind positiv semidefinit =⇒
P = STS, Q = RTR.
Numerisch robuster: Verwende S,R statt P,Q:
σ (SRT )2 = λ(PQ), cond(SRT
)=
√
cond (PQ).
Beachte: Mit SVD SRT = UΣV T gilt
S−T (PQ)ST = (SRT )(SRT )T
= (UΣV T )(V ΣUT ) = UΣ2UT
Berechne balancierende Transformation mit SVD:
SRT = [U1 U2]
[
Σ 00 Σ2
] [V T1
V T2
]
Σ = diag(σ21, . . . , σ2` )
⇓Tl = Σ
−1/2V T1 R, Tr = STU1Σ
−1/2.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 28
Modellreduktion Balanciertes Abschneiden
Balanciertes Abschneiden: hochdim. Systeme
Berechnung der SVD: Komplexitat O(n3).
Fur große Systeme: Gram’sche Matrizen haben nied-rigen numerischen Rang =⇒ benutze erneut Fakto-risierungsidee.
Beispiel:Lineare 1D Warmeleitungsgleichung mit Punktsteuerung,
Ω = [ 0, 1 ], FE-Diskretisierung mit linearen B-Splines,
h = 1/1000 (=⇒ n = 1001).
0 200 400 600 800 100010
−25
10−20
10−15
10−10
10−5
100
Index k
λ k, σk
Eigenwerte der Gram’schen, HSV
Λ (P)Λ (Q)Hankel−Singulaerwerte
P ≈ S(s)(S(s))T , S(s) ∈ Rs×n, s ≈ 20
Q ≈ R(s)(R(s))T , R(r) ∈ Rr×n, r ≈ 20
⇒ SVD von S(s)(R(s))T statt SRT ,⇒ Komplexitat O(105) statt O(1010).
♦ Peter Benner ♦ Institut fur Mathematik ♦ 29
Modellreduktion Balanciertes Abschneiden
Berechnung von S(s), R(s)
A dichtbesetzt:Iterative Losung basierend auf Matrix-Signumfunktion (Ã quadratische Konvergenz)
[B./Quintana-Ortı×2 ’97 ]
• Als “Black-Box”-Verfahren: KomplexitatO(n3)=⇒ Parallelisierung [B./Quintana-Ortı×2 ’99–’02 ]
• Daten-sparse Approximation von A mit hier-archischen Matrizen: Komplexitat O(n logk n)
[Hackbusch/Khoromskij/Grasedyck in Arbeit]
• Matrixkompression mit Mulitskalen-/Wavelet-techniken, formatierte Arithmetik???
A sparse:ADI-basierte Iteration
[Penzl ’99, Li/White ’01, B./Li/Penzl ]
Weiterentwicklung, Implementierung, Parallelisie-rung [B./Quintana-Ortı in Arbeit]
♦ Peter Benner ♦ Institut fur Mathematik ♦ 30
Modellreduktion Stochastisches Abschneiden
Modellreduktion basierend auf relativen Fehlern
Berechne reduziertes System, so daß der durch
G(s) = G(s)(I +∆rel).
definierte relative Fehler ‖∆rel‖∞ “klein” wird.
Balanciertes stochastisches Abschneiden (BST):
• Balanciere P nicht gegen Q, sondern gegenLosung W der ARG
0 = CT(DD
T)−1C + (A− B(DD
T)−1C)
TW +
+ W (A− B(DDT)−1C) +WB(DD
T)−1B
TW.
Berechnung mit Newton-ADI-Verfahren (modifi-ziert, da Lyapunovgleichungen nicht mehr semi-definit sind), erhalte faktorisierte Losung der ARG.
• Berechnung der Projektionsmatrizen uber SVDwie bei balanciertem Abschneiden.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 31
Modellreduktion Stochastisches Abschneiden
Vorteile von BST:
• Gleichmaßige Approximation der Transferfunktion uber ge-
samten Frequenzbereich.
• Erhalt minimale Phase (NST ≥ 0)
à erhalt Phasengang (im Gegensatz zu balanciertem
Abschneiden oder ”proper orthogonal decomposition”),
benotigt bei inversen Problemen, Parameteridentifikation.
• Erhalt robuste Stabilitat à Entwurf robuster Regelungen
à Reglerentwurf fur schwach-nichtlineare Systeme mit li-
nearem Modell moglich.
• Berechenbare Fehlerabschatzung (Ã adaptive Wahl des re-
duzierten Modells):
‖∆rel‖∞ ≤
n∏
j=`+1
1 + σj
1− σj
− 1,
wobei σ21, . . . , σ
2n = Λ (PW ), 1 ≥ σj ≥ σj+1 ≥ 0.
Numerische Verfahren:
A dichtbesetzt:
Kombination von Newton- und Sigumfunktionsmethoden,
Parallelisierung [B./Quintana-Ortı ×2 ’00]
A sparse
Kombination von modifizierten Newton- und ADI-Verfahren.
[B. in Arbeit]
♦ Peter Benner ♦ Institut fur Mathematik ♦ 32
Modellreduktion Numerische Ergebnisse
Numerische Ergebnisse
Testbeispiele:
1. Optimale Abkuhlung von Stahlprofilen (Modell von Man-
nesmann/Demag, [Troltzsch/Unger, Penzl 1999])
Wie bei LQR Problem.
2. Warmeleitungsgleichung mit Punktsteuerung
Beeinflussung der Temperaturverteilung in dunnem Draht
mit Warmequelle in der Mitte und gekuhltem Rand =⇒
1D Warmeleitungsgleichung mit homogenen Dirichlet-RBn,
diskretisiert mit FEM, lineare Elemente.
n = Dimension des FE Raums.
m = 1 Warmequelle an einem Punkt angelegt.
p = 1 Temperaturmessung an einem Punkt.
3. Simulation eines katalytischen Reaktors (Modell von ABB)
• FE Diskretisierung eines Randsteuerungsproblems
fur gekoppeltes PDE System (Bilanzgleichungen,
Reaktionsgleichungen, gemischte und Neumann-
Randbedingungen), Linearisierung um Arbeitspunkt.
• Dynamik: Oxidation (o-Xylene nach Anhydrit), Gaspha-
senreaktion im Reaktor.
• Steuerung: Externe Kuhlung des Reaktors.
• n = 1171, m = 6, p = 4.
4. “Zufallige” Systeme mit vorgegebenem McMillan Grad und
vorgegebenem Rang der Gram’schen Matrizen.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 33
Modellreduktion Numerische Ergebnisse
Bode-Diagramm (Frequenzgang, Amplitude)fur balanciertes Abschneiden
Beispiel: Gesteuerte 1D-Warmeleitungsgleichung
Rang (P ) = 32, Rang (Q) = 38(37), ` = 6
10−6
10−4
10−2
100
102
104
106
10−9
10−8
10−7
10−6
10−5
10−4
10−3
n = 500
frequency ω
|| G
(j ω)
− G
r(jω) |
| 2
SLICOT/AB09ADPDGEBTSRerror bound
10−6
10−4
10−2
100
102
104
106
10−9
10−8
10−7
10−6
10−5
10−4
10−3
frequency ω
|| G
(j ω)
− G
r(jω) |
| 2
n = 1000
SLICOT/AB09ADPDGEBTSRerror bound
♦ Peter Benner ♦ Institut fur Mathematik ♦ 34
Modellreduktion Numerische Ergebnisse
Speed-up/Effizienz der parallelen Algorithmen
Beispiel: Zufallig generierte Daten
n = 1000,m = p = 100, n = ` = 50
0 2 4 6 8 10 120
100
200
300
400
500
600
Number of processors
Exec
utio
n tim
e (s
ec.)
SLICOT/AB09AD
Beispiel: Steuerung einer katalytischen Reaktion
n = 1171,m = 6, p = 4, ` = 40
0 2 4 6 8 10 12 140
50
100
150
200
250
300
350
400
450
500
550
600
Number of processors
Exec
utio
n tim
e (s
ec.)
♦ Peter Benner ♦ Institut fur Mathematik ♦ 35
Modellreduktion Numerische Ergebnisse
Gleichmaßige Approximation durch BST
Beispiel: Abkuhlung von Stahlprofilen
• n = 821: Rang (P ) = 165, Rang (Q) = 210, ` = 40
Newton+Signumfunktion
10−2
100
102
104
10−5
10−4
10−3
10−2
10−1
100
Frequency ω (rad./sec.)
Gai
n (m
agni
tude
)
From u1 to y
1
full order systemPDGESTSRAB09HD
Newton+ADI
10−4
10−2
100
102
104
106
10−5
10−4
10−3
10−2
10−1
100
Frequency (ω rad/sec.)
Gai
n
From u3 to y
1
full−order reduced−order
• n = 3113: Rang (P ) = 178, Rang (Q) = 220
♦ Peter Benner ♦ Institut fur Mathematik ♦ 36
Modellreduktion Fazit/Ausblick
Fazit/Ausblick
• Modellreduktion fur diskrete Systeme analog; Be-rechnung der Faktoren der Gram’schen Matri-zen mit faktorisierter Smith-Iteration, sowohl furdichtbesetzte als auch sparse Systeme.
• Implementierung der Methoden fur sparse Syste-me fur SLICOT.
• Aufbau der Softwarebibliothek PLiCMR, Integra-tion in parallele Version von SLICOT.
• Modellreduktion großer Systeme auf PC Clustermit ”remote computing”: E-mail, Web Server.
• Schaltkreissimulation:
– Berechnung passiver reduzierter Systeme (“po-sitive real balancing”).
– Ubertragung auf DAE Systeme.
• Berucksichtigung von PDE Strukturen?
SLICOT = Subroutine Library in Control Theory,
Softwarebibliothek mit ca. 300 Routinen fur Systemanalyse, Systemidentifi-
kation, Reglerentwurf, etc.
♦ Peter Benner ♦ Institut fur Mathematik ♦ 37
Modellreduktion Fazit/Ausblick
Can we beat the curse of
dimensionality?
Ja,
• wenn wir bereit sind, die Struktur der Problemezu analysieren und zu verwenden,
• wenn die “richtige” numerische lineare Algebraverwendet wird:
– sparse Matrix-Algorithmen,– parallele Algorithmen.
Und dafur mussen wir noch nicht mal auf Matlab
verzichten!
♦ Peter Benner ♦ Institut fur Mathematik ♦ 38