at 00/0000

Modellreduktion fu ¨r parametrisierte Systeme durch balanciertes Abschneiden und Interpolation Model Reduction for Parametric Systems Using Balanced Truncation and Interpolation Ulrike Baur, Peter Benner, Technische Universit¨at Chemnitz

Der Beitrag stellt einen neuen Ansatz zur Parameter-erhaltenden Modellreduktion parametrisierter, linearer Systeme vor. Das Verfahren ist eine Kopplung des balancierten Abschneidens mit Interpolation, so dass Fehlerabsch¨atzungen f¨ ur die Qualit¨at des reduzierten Systems ableitbar sind. Durch die Verwendung d¨ unner Gitter zur Diskretisierung des Parameterraumes kann das Verfahren auch auf Systeme mit mehreren Parametern angewendet werden. This paper presents a parameter-preserving approach for model order reduction of parameterized, linear systems. The proposed method is a coupling of balanced truncation with interpolation. Thus, error estimates for the quality of the reduced-order system can be derived. The method can be applied to systems containing several parameters by using sparse grids for the discretization of the parameter space. Schlagw¨ orter: Parametrische Modellreduktion, balanciertes Abschneiden, Interpolation, d¨ unne Gitter. Keywords: Parametric model order reduction, balanced truncation, interpolation, sparse grids.

1 Einleitung Heutzutage hat die numerische Simulation in vielen Anwendungsgebieten das Experiment abgel¨ ost. Dabei wird das zugrundeliegende dynamische System durch mathematische Modellierung in ein System aus gew¨ ohnlichen (ODEs) oder partiellen (PDEs) Differentialgleichungen u uhrt. Die Finite-Elemente-Methode (FEM) ¨ berf¨ erm¨ oglicht das Generieren genauer Modelle, z. B. zur Beschreibung des W¨ armeflusses in einem Bauteil. Die typischerweise große Anzahl auftretender Gleichungen wird als Ordnung oder Dimension des Systems bezeichnet. Sie verhindert den Einsatz der FEM in Simulationen auf Systemebene. Daher gilt es, die Ordnung durch geeignete Verfahren signifikant zu reduzieren, um eine anschließende numerische Simulation effizient durchf¨ uhren zu k¨ onnen. Allerdings haben die reduzierten Systeme, wie sie mit den g¨ angigen Modellreduktionsverfahren berechnet werden, den großen Nachteil, kaum Variationen am urspr¨ unglichen physikalischen

Modell zuzulassen. In vielen Anwendungen, wie z. B. im Designprozess von Mikrosystemen, werden Modifikatio¨ nen aber ben¨otigt, um Anderungen an den Randbedingungen oder an der Geometrie schnell in den Simulationsprozess u ¨bernehmen zu k¨onnen. Daher besteht der Bedarf nach neuen parametrischen Modellreduktionsverfahren, die die Ordnung von parametrisierten Systemen zuverl¨assig reduzieren und zus¨atzlich die Parameter im reduzierten System erhalten. Die Entwicklung von Verfahren dieser Art steht noch am Anfang. Die meisten der bislang entwickelten Ans¨atze basieren auf multivariater Pad´e-Approximation [6, 8, 9, 10, 11, 12, 25, 30]. Sie unterscheiden sich in der Berechnungsart der Momente (implizit/explizit) und in der Auswahl der (gemischten) Momente, die f¨ ur die Berechnung des reduzierten Systems verwendet werden. Die explizite Berechnung der Momente f¨ uhrt hier, wie in vergleichbaren Verfahren f¨ ur nicht-parametrisierte Systeme, zu numerischen Instabilit¨aten. Daher werden in der Praxis meist

c Oldenbourg Verlag at – Automatisierungstechnik 48 (2000) 00

00

1

implizite Ans¨atze verwendet. F¨ ur Systeme, die nur wenige (ein bis drei) Parameter enthalten, liefern diese Verfahren gute Ergebnisse. Allerdings existiert noch keine solide mathematische Basis, d. h. es existieren keine Fehlersch¨atzer (nicht einmal lokal) und es ist unklar, welche und wieviele Momente in die Berechnung eingehen sollen. In dieser Arbeit wird ein anderer Ansatz zur parametrischen Modellreduktion vorgestellt. Er basiert auf einer Kopplung der Methode des balancierten Abschneidens (Balanced Truncation – BT), welche sehr erfolgreich zur Modellreduktion linearer, nicht-parametrisierter Systeme (im folgenden mit deterministisch bezeichnet) eingesetzt wird, mit (st¨ uckweise) polynomialer Interpolation. Durch den systemtheoretischen Hintergrund des balancierten Abschneidens k¨ onnen wir erste Fehlerschranken f¨ ur die Methode herleiten. Desweiteren erlaubt der Gebrauch der D¨ unngittertechnik die Anwendung des Verfahrens auf Systeme mit einer gr¨ oßeren Anzahl von Parametern.

ˆ ˆ ˆ mit (A(p), B(p), C(p)) ∈ Rr×r × Rr×m × Rr×q ersetzt, wobei r ≪ n und ky − yˆk < tol · kuk gelten soll. Ein weit verbreitetes Verfahren zur Modellreduktion linearer, deterministischer Systeme ist das balancierte Abschneiden [23]. Im folgenden wird ein HybridVerfahren aus balanciertem Abschneiden und Interpolation f¨ ur die Ordnungsreduktion parameterabh¨angiger Systeme vorgestellt. Zur besseren Veranschaulichung betrachten zun¨achst den eindimensionalen Parameterraum. Erweiterung des Ansatzes auf h¨oherdimensionale rameterr¨aume wird im darauffolgenden Abschnitt schrieben.

wir Die Pabe-

3 Hybrid-Verfahren zur Modellreduktion parameterabh¨ angiger Systeme 3.1 Der eindimensionale Parameterraum

2 Grundlagen

In diesem Abschnitt beschr¨anken wir uns auf Systeme der Form (1) mit Abh¨angigkeit von nur einem Parameter p ∈ [a, b] und Λ (A(p)) ⊂ C− f¨ ur alle p ∈ [a, b].

Parametrisierte, lineare Systeme k¨ onnen in der folgenden allgemeinen Form geschrieben werden:

Die vorgeschlagene Methode besteht aus den folgenden Rechenschritten:

x(t) ˙ = A(p) x(t) + B(p) u(t), y(t) = C(p)T x(t),

1. W¨ahle Interpolationsst¨ utzstellen p0 , . . . , pk ∈ [a, b]. 2. Auf die k + 1 deterministischen Systeme mit ¨ Ubertragungsfunktionen

(1)

mit (A(p), B(p), C(p)) ∈ Rn×n × Rn×m × Rn×q , wobei die Abh¨angigkeit der Systemmatrizen vom Parametervektor p ∈ Rd auch nichtlinear sein kann. Mit x(t) ∈ Rn wird der Zustandsvektor bezeichnet, u(t) ∈ Rm ist der Eingangs-, y(t) ∈ Rq der Ausgangsvektor. Um die Darstellung zu vereinfachen, ignorieren wir die Parameterabh¨angigkeit des Zustands und des Ausgangs beim Gebrauch der Notation x(t) := x(t, p),

y(t) := y(t, p).

In vielen Anwendungen ist die Zustandsraumdimension n und damit die Anzahl der ODEs in (1) sehr groß und m, q ≪ n. Es wird vorausgesetzt, dass A(p) f¨ ur alle betrachteten Parameterwerte stabil ist, d. h. alle Eigenwerte von A(p) (bezeichnet mit Λ (A(p))) sind in der offenen linken Halbebene C− . Durch den Einsatz von Modellreduktionsverfahren soll die Anzahl der ODEs reduziert werden, bei m¨oglichst guter Abbildung des Eingangs-, Ausgangsverhaltens, ¨ welches durch die Ubertragungsfunktion T

G(s, p) = C(p) (sIn − A(p))

−1

2

j = 0, . . . , k, wende die Methode des balancierten Abschneidens an. Die resultierenden Systeme rj -ter Ord¨ nung mit den dazugeh¨origen Ubertragungsfunktionen ˆj ˆ j (s) = CˆjT (sIrj − Aˆj )−1 B G sind stabil und erf¨ ullen f¨ ur vorgegebene Toleranz ǫ die BT H∞ -Fehlerschranke [13]:   n X ˆ j k∞ ≤ 2  kGj − G σi  < ǫ, (3) i=rj +1

wobei rj minimal gew¨ahlt ist, so dass (3) erf¨ ullt ist. 3. Durch Interpolation, z. B. nach dem Ansatz von Lagrange, erh¨alt man eine Interpolationsfunktion ˆ p) = G(s,

B(p)

ˆ x ˆ A(p) ˆ(t) + B(p) u(t), T ˆ C(p) xˆ(t),

k X

ˆ j (s) lj (p)G

(4)

j=0

beschrieben wird. Zus¨ atzlich sollen alle Parameter im reduzierten System erhalten bleiben. In nachfolgenden numerischen Simulationen wird dann das System (1) durch das reduzierte System x ˆ˙ (t) = yˆ(t) =

Gj (s) := G(s, pj ) = C(pj )T (sIn − (A(pj )))−1 B(pj ),

(2)

mit Lagrange-Polynomen lj (p) =

k Y

i=0,i6=j

p − pi pj − pi

f¨ ur j = 0, . . . , k, die aus den reduzierten Systemmatrizen zusammengesetzt ist.

¨ Dabei kann die Ubertragungsfunktion des parameterabh¨ angigen, reduzierten Systems in (4) wie folgt dargestellt werden: ˆ p) = G(s,

k X

ˆj lj (p)CˆjT (sIrj − Aˆj )−1 B

j=0

=

k X

ˆj Cˆj (p)T (sIrj − Aˆj )−1 B

in Form von Niedrigrangfaktoren approximieren [1, 2, 4, 5, 16, 17, 19, 22, 26, 27, 28]. Durch die Kombination der bekannten Fehlerschranken f¨ ur das balancierte Abschneiden mit den Interpolationsfehlerabsch¨atzungen ist es m¨oglich, eine Absch¨atzung f¨ ur den Simulationsfehler im Parameterintervall [a, b] anzugeben. Unter der Annahme

j=0

ˆ T (sIr − A) ˆ −1 B, ˆ = C(p)

mit r =

k X

G(s, ·) ∈ C k+1 ([a, b] → Cq×m ) rj

erhalten wir folgende Absch¨atzung f¨ ur den Fehler zwi¨ schen Original- und reduzierter Ubertragungsfunktion:

j=0

und   ˆ   C1 (p)T k Y p − p   i . T ˆ T :=   CˆjT , .. C(p)  , Cˆj (p) :=  pj − pi i=0,i6=j Cˆk (p)T 

 Aˆ := 

(sIr0 − Aˆ0 )

−1

..

ˆ0 B  ˆ  ..   , B :=  .  . −1 ˆk − Aˆk ) B 

. (sIrk





ˆ p)k sup kG(s, p) − G(s,

Man hat somit ein reduziertes System zur Verf¨ ugung, mit welchem die parameterabh¨ angige ¨ Ubertragungsfunktion des Originalmodells (1) f¨ ur jeden Parameter im Interpolationsintervall [a, b] approximiert ˆ p) die Interpolatiwerden kann. Zus¨atzlich erf¨ ullt G(s, onsbedingungen ˆ pj ) = G ˆ j (s) G(s, f¨ ur j = 0, . . . , k. Aufgrund der Blockdiagonalgestalt von ˆ p) reduziert sich der Aufwand zur Auswertung in G(s, der numerischen Simulation von O(n3 ) f¨ ur G(s, p) auf ˆ p). O((k + 1) max(rj )3 ) beim Gebrauch von G(s, Zur Ordnungsreduktion der deterministischen Systeme im Schritt 2 des Verfahrens k¨ onnen f¨ ur kleinere Systemdimensionen (bis etwa n = 2000) klassische Varianten des BT Verfahrens [29] eingesetzt werden. Diese finden sich z. B. in der Subroutine Library in Systems and Control Theory SLICOT1 . F¨ ur gr¨oßere Dimensionen, wie sie z. B. bei der Simulation mit Finite-Elemente Modellen in der Mikrosystemtechnik auftreten, bieten sich neuere Entwicklungen an. Diese Implementierungen des BT Verfahrens unterscheiden sich vor allem in der L¨ osung der auftretenden Lyapunovgleichungen, dem aufw¨ andigsten Rechenschritt beim balancierten Abschneiden. F¨ ur große Systeme, deren Systemmatrix A(p) schwachbesetzt (oder daten-sparse) ist, bieten sich iterative Verfahren zur L¨ osung der Lyapunovgleichungen an, welche die L¨ osung http://www.slicot.org

(5)

s∈C+

p∈[a,b]

= sup kG(s, p) − s∈C+

k X

ˆ j (s)k lj (p)G

k X

lj (p)Gj (s)k +

j=0

p∈[a,b]

≤ sup kG(s, p) − s∈C+

Alternativ kann die Parameterabh¨ angigkeit auch in die Diagonalbl¨ocke (sIrj − Aˆj )−1 von Aˆ oder in den Vektor ˆ verschoben werden. B

1

∀s ∈ C+

j=0

p∈[a,b]

sup k s∈C+

p∈[a,b]

k X

ˆ j (s))k lj (p)(Gj (s) − G

j=0

(3)

≤ sup kRk (G, s, p)k + ǫ sup |

k X

lj (p)|.

(6)

p∈[a,b] j=0

s∈C+

p∈[a,b]

Hierbei bezeichnet Rk (G, s, p) = G(s, p)−

k X

lj (p)Gj (s)

j=0

das Restglied mit der Darstellung Rk (G, s, p) =

1 (k + 1)!



Y k ∂ k+1 (p − pi ) G(s, ξ(p)) ∂pk+1 i=0

an einer Zwischenstelle ξ(p) ∈ [minj pj , maxj pj ]. Durch die spezielle Wahl der St¨ utzstellen als Nullstellen des Tschebyscheff-Polynoms kann der Interpolationsfehleranteil in (6) bzgl. p minimiert werden [18]. Alternativ zur polynomialen Interpolation w¨are auch der Gebrauch rationaler Interpolation denkbar. So ist G(s, p) bei polynomialer Abh¨angigkeit der Matrizen vom Parameter auch rational in p, was rationale Interpolation sinnvoll erscheinen l¨asst. Dies wird in der Zukunft n¨aher untersucht.

3.2 Mehrdimensionaler Parameterraum Die Komplexit¨at des f¨ ur den eindimensionalen Parameterraum geschilderten Hybrid-Verfahrens zur Modellreduktion h¨angt neben der Systemgr¨oße stark von der Anzahl der durchgef¨ uhrten Modellreduktionen f¨ ur deterministische Systeme, also von der Anzahl der Interpolationsst¨ utzstellen, ab.

3

Points: 5, Level: 1 Points full grid: 9

T1

Points: 13, Level: 2 Points full grid: 25 1

1

0.8

0.8

0.8

0.8

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0

0.2

0.4

0.6

0.8

1

0

0

0.2

Points: 7, Level: 1 Points full grid: 27

0.4

0.6

0.8

0

1

0.2

0

0.2

Points: 25, Level: 2 Points full grid: 125

0.4

0.6

0.8

1

0

0

1

1

1

0.5

0.5

0.5

0 1

0 1

1 0 0

0.6

0.8

1 0.5

0.5

1

0 1

1 0.5

0.5

0.4

Points: 177, Level: 4 Points full grid: 4913

1

0.5

0.2

Points: 69, Level: 3 Points full grid: 729

0.5

0 1

T3

Points: 29, Level: 4 Points full grid:289

1

0

T2

Points: 29, Level: 3 Points full grid: 81

1

0 0

0.5

1 0.5

0.5

0 0

Bild 2: MATLAB Sparse Grid Interpolation Tschebyscheff-Gauss-Lobatto Gitter.

0 0

Toolbox:

Zur Interpolation ausreichend glatter Funktionen im mehrdimensionalen Parameterraum [0, 1]d ,

Bild 1: Basisfunktionen von T1 , T2 , T3 .

G(s, ·) : [0, 1]d → R, Daher schlagen wir zur mehrdimensionalen Interpolation den Gebrauch von d¨ unnen Gittern [7, 14, 32] f¨ ur die Diskretisierung des Parameterraums vor, die im folgenden kurz eingef¨ uhrt werden. Wir betrachten zur besseren Veranschaulichung zun¨achst den eindimensionalen Raum [0, 1], welchen wir mit einem ¨ aquidistanten (vollen) Gitter der Maschenweite hℓ = 2−ℓ unterteilen. Dazu assoziiert ist der (2ℓ − 1)-dimensionale Raum der stetigen, st¨ uckweise linearen Funktionen Sℓ , die auf dem Rand verschwinden. Dieser Raum kann nun unter Zuhilfenahme der Teilr¨aume T1 , . . . , Tℓ , siehe Abbildung 1, hierarchisch zerlegt werden [31]: Sℓ = T 1 ⊕ · · · ⊕ T ℓ . Dabei ist ein Teilraum Tj derart definiert, dass er alle stetigen, st¨ uckweise linearen Funktionen auf dem Gitter der Maschenweite 2−j enth¨ alt, die an den Gitterpunkten aller gr¨oberen Gitter verschwinden. Es gilt: dim Tj = 2j−1 . 2

F¨ ur eine Funktion G(s, ·) ∈ C ([0, 1]), die auf dem Rand verschwindet (dies dient ausschließlich zur ¨ Ubersichtlichkeit der Darstellung und ist in der Praxis nicht erforderlich), wird der Interpolant GI im Raum Sℓ gesucht: ℓ X gj , mit gj ∈ Tj . GI = j=1

∂ 2d G(s, ·) ∈ C 0 ([0, 1]d ), ∂p21 . . . ∂p2d

werden im folgenden Teilr¨aume Tj mit geringem Anteil am Interpolationsfehler (7) weggelassen. Dies bedeutet, dass der Interpolant GI nicht mehr im (2ℓ − 1)d dimensionalen Raum der st¨ uckweise d-linearen Funktionen Sℓ =

ℓ M

j1 =1

···

ℓ M

Tj ,

j := (j1 , . . . , jd ) ∈ Nd ,

jd =1

gesucht wird, sondern in S˜ℓ =

M

Tj .

|j|1 ≤ℓ+d−1

Damit einhergehend reduziert sich die Anzahl der Gitterpunkte des urspr¨ unglich rechtwinkligen Gitters O(h−d unnes Gitℓ ) und wir erhalten ein sogenanntes d¨ −1 d−1 ter der Dimension O(h−1 ). Der Interpolaℓ (log(hℓ )) tionsfehler ist nahezu unver¨andert beschr¨ankt bei d−1 |G(s, ·) − GI |∞ ≤ O(h2ℓ (log(h−1 ). ℓ ))

Alternativ zu st¨ uckweise linearen Funktionen k¨onnen auch st¨ uckweise polynomiale Funktionen zur Interpolation verwendet werden [3]. Wie bei eindimensionaler Interpolation mit Polynomen bietet sich hier die Wahl von Tschebyscheff-St¨ utzstellen an Stelle von ¨aquidistanten Abst¨anden an. Das daraus resultierende d¨ unne Tschebyscheff-Gauss-Lobatto Gitter ist f¨ ur d = 2, 3 und verschiedene Gittertiefen ℓ = 1, · · · , 4 in Bild 2 dargestellt.

Mit der Halbnorm |G(s, ·)|∞ := max

 2 ∂ G(s, ·) ∂p2

 0 ≤ p ≤ 1 ,

4 Numerische Beispiele

erhalten wir die Schranke f¨ ur den Interpolationsfehler |G(s, ·) − GI |∞ ≤ O(h2ℓ ),

(7)

wobei der Anteil eines Teilraums Tj am Interpolanten abgesch¨atzt werden kann: kgj k∞

4

1 ≤ 4−j |G(s, ·)|∞ . 2

4.1 Eindimensionaler Parameterraum 4.1.1 Konvektions-Diffusionsgleichung Die Konvektions-Diffusionsgleichung in Ω = (0, 1)2 ∂x (t, ξ) = ∆x(t, ξ) + p · ∇x(t, ξ) + b(ξ)u(t), ∂t x(t, ξ) = 0, ξ ∈ ∂Ω,

ξ ∈ Ω,

| G (j ω,p) − Gr (j ω,p) |

−4

10

Bild 4: Anemometer Modell.

−6

10

−8

10 10

8

6

10

6

4

10

4

2

10 2

p

0

10 0

−2

10

ω

ˆ Bild 3: Fehler zwischen G(ω, p) und G(ω, p). wird mit der Finiten-Differenzen-Methode und n = 400 Gitterpunkten approximiert. Der Parameter p ∈ [0, 10] bestimmt den konvektiven Transport in eine Koordinatenrichtung. Das resultierende ODE System x(t) ˙ = (A + pA1 ) x(t) + b u(t) und eine Ausgangsgleichung y(t) = cT x(t) beinhalten Vektoren b, c ∈ Rn mit bT = [1, 0, · · · , 0], cT = [1, 1, · · · , 1]. Bei einer Wahl von 6 TschebyscheffSt¨ utzstellen im Parameterintervall [0, 10] und einer Fehlertoleranz-Vorgabe von ǫ = 10−4 f¨ ur das balancierte Abschneiden werden reduzierte Systeme der Ordnungen rj ∈ {3, 4} f¨ ur j = 1, . . . , 6 berechnet. Der Fehˆ p) nach Interpolation mit ler zwischen G(s, p) und G(s, Lagrange-Polynomen ist in Bild 3 dargestellt. Der Gesamtfehler (5) wurde hier auf einem 50×30 Gitter approˆ ximiert. Das bedeutet, der Fehler |G(ω, p) − G(ω, p)| wurde f¨ ur 50 Frequenzwerte in logarithmischer Skala zwischen 10−2 und 106 und f¨ ur 30 Parameterwerte im Intervall [0, 10] berechnet. Es zeigt sich eine gute Approximation durch das reduzierte System. Der maximale Wert auf dem gew¨ahlten Gitter betr¨ agt 3.3 × 10−5 .

4.1.2 Str¨omungssensor (Anemometer) Eine weit verbreitete Methode zur Messung von Str¨ omungsraten mittels mikrosystemtechnischer Sensoren ist eine Kombination aus einem Heizer und Thermoelementen zu beiden Seiten des Heizers, die von einer Fl¨ ussigkeit oder einem Gas u omt werden2 , ¨ berstr¨ siehe [24] und die dort angegebenen Referenzen. Diese Str¨omung ver¨andert die Temperaturverteilung, indem der symmetrische Diffusionsprozess von einem unsymmetrischen Konvektionsprozess u ¨ berlagert wird. Somit entsteht eine Temperaturdifferenz zwischen den Sensorpositionen, die von der Flussrate abh¨ angt und mit Thermoelementen detektiert werden kann. Dadurch ist es m¨oglich, aus den Daten der Thermoelemente die Flussrate zu bestimmen. Im Designprozess ist hierbei 2

ein kompaktes Modell n¨otig, das die Flussgeschwindigkeit sowie evtl. Materialparameter als Parameter bereitstellt.

http://www.imtek.de/simulation/benchmark

Um ein Modell zu generieren, welches die kompletten physikalischen Effekte beinhaltet, ist es in diesem Beispiel n¨otig, die Navier-Stokes Gleichungen gekoppelt mit der W¨armeleitungsgleichung zu l¨osen. Dies ist allerdings nur dann n¨otig, wenn die eingebrachte W¨arme das Str¨omungsprofil merklich ver¨andert. Da die W¨armemenge bei der vorgeschlagenen Anwendung aber relativ klein ist, kann man hier davon ausgehen, dass das von außen vorgegebene Geschwindigkeitsfeld kaum ver¨andert wird. Somit reduziert sich das mathematische Modell zu linearem Konvektions-DiffusionsW¨armetransport wie in Beispiel 4.1.1. Die Systemmatrizen wurden f¨ ur die zugrundeliegende Geometrie (siehe Bild 4) mit der Finite-Elemente Methode in der Simulationsumgebung ANSYS3 erstellt. Das resultierende System unterscheidet sich vom vorigen Beispiel durch eine zus¨atzliche Massematrix; der Parameter p ∈ [0, 1] beeinflusst weiterhin die Flussgeschwindigkeit in horizontaler Richtung. Es wurde ein feines Gitter mit n = 29008 Freiheitsgraden verwendet. Zur Ordnungsreduktion haben wir 12 TschebyscheffSt¨ utzstellen im Parameterintervall [0, 1] gew¨ahlt und an diesen Stellen das System mit der Methode des balancierten Abschneidens und einer Fehlertoleranz von ǫ = 10−4 reduziert. Dies ergab 12 reduzierte Systeme 12 X rj = 75. Das pamit Ordnungen 2 ≤ rj ≤ 9 und r = i=1

rameterabh¨angige reduzierte System nach der LagrangeInterpolation approximiert das urspr¨ ungliche System gut u ur den ma¨ ber den gesamten Parameterbereich. F¨ ximalen Fehler (5), der wieder auf einem 50 × 30 Gitter approximiert wird (siehe Bild 5), erhalten wir 6.4×10−4.

4.2 Mehrdimensionaler Parameterraum 4.2.1 Konvektions-Diffusionsgleichung mit d = 2 Wir betrachten die Konvektions-Diffusionsgleichung aus Beispiel 4.1.1, einschließlich der gleichen Ortsdiskretisierung mit n = 400 Gitterpunkten. Zus¨atzlich wird der konvektive Transport nun in beide Koordinatenrichtungen variiert, beim Gebrauch der Parameter p1 , p2 ∈ [0, 1] im ODE System x(t) ˙ = (A + p1 A1 + p2 A2 ) x(t) + b u(t). F¨ ur die Diskretisierung des Parameterraums und einer bilinearen Interpolation auf dem resultierenden d¨ unnen 3

http://www.ansys.com

5

| G (j ω,p) − Gr (j ω, p) |

freq 1.00e−02

freq 2.40e−01

−4

−3

10

freq 1.08e+01

−4

x 10

−4

x 10

x 10

1

1

1

0.5

0.5

0.5

−4

10

0 1

1

0.5 p

0 0

2

0 1

0 1

1

0.5 p

0.5 p1

0 0

2

1

0.5 p

0.5 p1

0 0

2

0.5 p1

−6

10

freq 4.89e+02

freq 2.21e+04

−4

freq 1.00e+06

−4

x 10

−4

x 10

x 10

1

1

1

0.5

0.5

0.5

0 1

0 1

−8

10

1 0.8

6

10

0.6

4

10

0.4

2

10 0.2

10 0

p

0.5 p2

0

−2

10

ω

ˆ Bild 5: Fehler zwischen G(ω, p) und G(ω, p).

1 0 0

0.5 p2

0.5 p1

0 1

1

0.5 p2

0.5 p1

0 0

1 0 0

0.5 p1

Bild 7: Fehler f¨ ur verschiedene Frequenzen.

| G (j ω,p) − Gr (j ω,p) | −4

10 −5

4.5 4 3.5 3 2.5

p

2 1

max | G (j ω,p) − Gr (j ω,p) |

x 10

0.8

1 0.6

0.8

−7

10

0.4

0.2

2

−6

10

0.6

0.4

p

−5

10

0.2 0

0

−8

p1

ˆ Bild 6: Fehler zwischen G(ω, p) und G(ω, p). Gitter haben wir die MATLAB Sparse Grid Interpolation Toolbox [20, 21] verwendet. Bei Vorgabe einer maximalen Tiefe von ℓ = 1, demzufolge einer Interpolation u unf D¨ unngitterpunkten, erhalten wir eine gute Ap¨ ber f¨ ¨ proximation der urspr¨ unglichen Ubertragungsfunktion im gesamten Parameterraum mit gesch¨ atztem Interpolationsfehler von 1.8 × 10−4 . Dieser Wert wird von der Toolbox bereitgestellt, gesch¨ atzt durch Ausnutzen der hierarchischen Struktur des Gitters. Die BT Modellreduktion wurde an den D¨ unngitterpunkten mit Vorgabe einer Fehlertoleranz von ǫ = 10−4 durchgef¨ uhrt, was zu Systemen reduzierter Ordnungen rj = 3 f¨ ur j = 1, . . . , 5 ˆ wird als Maximum f¨ uhrte. Der Fehler zwischen G und G der absoluten Fehler zu 30 verschiedenen Frequenzwerten auf einem 32 × 32 Gitter im Parameterraum approximiert, siehe Bild 6. F¨ ur ausgew¨ ahlte Frequenzen zeigt Bild 7 den absoluten Fehler. Der maximale Fehler u ¨ber den gesamten Parameterraum in Abh¨ angigkeit von der Frequenz ist in Bild 8 dargestellt.

4.2.2 Konvektions-Diffusionsgleichung mit d = 3

10

−2

0

10

2

10

4

10

6

10

Frequenz ω

10

Bild 8: Maximaler Fehler im Parameterraum.

und st¨ uckweise polynomialen Interpolationspolynomen [3] durchgef¨ uhrt. Eine vorgegebene absolute Toleranz f¨ ur die Gitterverfeinerung von 10−4 f¨ uhrt zu einer Gittertiefe von ℓ = 2 und 25 Interpolationspunkten. Die vorgegebene Genauigkeit bei der Modellreduktion ǫ = 10−4 ergibt reduzierte Systeme der Gr¨oße rj = 3 f¨ ur j = 1, . . . , 25. Der gesch¨atzte Interpolationsfehler ist 6.6 × 10−4 . Approximationen des Fehlers (5) sind f¨ ur feste Parameterwerte p0 = 1 und p0 = 1.5 (auf einem 32 × 34 Gitter f¨ ur die Parameter p1 und p2 ) in den Bildern 9 und 10 dargestellt. Den approximierten maximalen Fehler u ¨ ber den gesamten Parameterraum in Abh¨angigkeit von der Frequenz zeigt Bild 11. | G (j ω,p) − Gr (j ω,p) | für p0 = 1.0 −5

x 10 2.8 2.7 2.6 2.5 2.4

Nach Hinzunahme eines weiteren (k¨ unstlichen) Parameters p0 ∈ [1, 1.5] erhalten wir ein System mit dreidimensionalem Parameterraum und

2.3 2.2 0 0.2 0.4

x(t) ˙ = (p0 A + p1 A1 + p2 A2 ) x(t) + b u(t).

1 0.8

0.6

0.6

0.8

0.4 1

Die Interpolation im D¨ unngitterraum wird f¨ ur dieses Beispiel mit der Wahl von Tschebyscheff-St¨ utzstellen

6

p

1

0.2 0

p

2

ˆ Bild 9: Fehler zwischen G(ω, p) und G(ω, p).

| G (j ω,p) − Gr (j ω,p) | für p0 = 1.5

−5

x 10 2.6 2.4 2.2 2 1.8 1.6 0

1

0.2 0.8

0.4

0.6

0.6

0.4 0.8

0.2 1

0

p

p

2

1

ˆ Bild 10: Fehler zwischen G(ω, p) und G(ω, p). −4

p

max | G (j ω,p) − Gr (j ω,p) |

10

−5

10

−6

10

−7

10

−8

10

−2

10

0

10

2

10

Frequenz ω

4

10

6

10

Bild 11: Maximaler Fehler im Parameterraum.

5 Ausblick Neben der polynomialen Interpolation soll als weitere M¨ oglichkeit im eindimensionalen Parameterraum die rationale Interpolation untersucht werden. F¨ ur die mehrdimensionale Interpolation muss im folgenden untersucht werden, wie eine explizite Darstellung ¨ der Ubertragungsfunktion bzw. eine Transformation des ˆ p) in ein Zustandsraummodell erhalten Modells“ G(s, ” werden kann. An Stelle des balancierten Abscheidens soll weiterhin der Einsatz anderer Modellreduktionsverfahren im Frequenzraum, wie z. B. der H2 -optimalen Modellreduktion [15], an den Interpolationsst¨ utzstellen getestet werden. Literatur [1] J. Bad´ıa, P. Benner, R. Mayo, and E. Quintana-Ort´ı. Solving large sparse Lyapunov equations on parallel computers. In B. Monien and R. Feldmann, editors, Euro-Par 2002 – Parallel Processing, number 2400 in Lecture Notes in Computer Science, pages 687–690. Springer-Verlag, Berlin, Heidelberg, New York, 2002. [2] J. Bad´ıa, P. Benner, R. Mayo, and E. Quintana-Ort´ı. Parallel algorithms for balanced truncation model reduction of sparse systems. In J. Dongarra, K. Madsen, and J. Wasniewski, editors, Applied Parallel Computing: 7th International Conference, PARA 2004, Lyngby, Denmark, June 20-23, 2004, number 3732 in Lecture Notes

in Computer Science, pages 267–275. Springer-Verlag, Berlin, Heidelberg, New York, 2006. [3] V. Barthelmann, E. Novak, and K. Ritter. High dimensional polynomial interpolation on sparse grids. Adv. Comput. Math., 12(4):273–288, 2000. [4] U. Baur and P. Benner. Gramian-based model reduction for data-sparse systems. SIAM J. Sci. Comput., 31(1):776–798, 2008. [5] P. Benner. Solving large-scale control problems. IEEE Control Systems Magazine, 14(1):44–59, 2004. [6] B. Bond and L. Daniel. Parameterized model order reduction of nonlinear dynamical systems. In Proceedings of the IEEE Conference on Computer-Aided Design, November 2005. [7] H. Bungartz. D¨ unne Gitter und deren Anwendung bei der adaptiven L¨ osung der dreidimensionalen PoissonGleichung. Dissertation, Inst. f. Informatik, TU M¨ unchen, June 1992. [8] L. Daniel, O. Siong, L. Chay, K. Lee, and J. White. A multiparameter moment-matching model-reduction approach for generating geometrically parameterized interconnect performance models. IEEE Trans. Comput.Aided Design Integr. Circuits Syst., 23(5):678–693, May 2004. [9] R. Eid, B. Salimbahrami, B. Lohmann, E. B. Rudnyi, and J. G. Korvink. Parametric order reduction of proportionally damped second-order systems. Sensors and Materials, Tokyo, 19(3):149–164, 2007. [10] O. Farle, V. Hill, P. Ingelstr¨ om, and R. DyczijEdlinger. Ordnungsreduktion linearer zeitinvarianter Finite-Elemente-Modelle mit multivariater polynomieller Parametrierung (Model Order Reduction of Linear Finite Element Models Parameterized by Polynomials in Several Variables). at-Automatisierungstechnik, 54(4):161–169, 2006. [11] L. Feng. Parameter independent model order reduction. Math. Comput. Simulation, 68(3):221–234, 2005. [12] L. Feng and P. Benner. A robust algorithm for parametric model order reduction based on implicit moment matching. PAMM Proc. Appl. Math. Mech., 7:1021501– 1021502, 2008. [13] K. Glover. All optimal Hankel-norm approximations of linear multivariable systems and their L∞ norms. Internat. J. Control, 39:1115–1193, 1984. [14] M. Griebel. A parallelizable and vectorizable multi-level algorithm on sparse grids. In Parallel algorithms for partial differential equations (Kiel, 1990), volume 31 of Notes Numer. Fluid Mech., pages 94–100. Vieweg, Braunschweig, 1991. [15] S. Gugercin, A. C. Antoulas, and C. Beattie. H2 Model Reduction for large-scale dynamical systems. SIAM J. Matrix Anal. Appl., 30(2):609–638, 2008. [16] S. Gugercin and J.-R. Li. Smith-type methods for balanced truncation of large systems. In P. Benner, V. Mehrmann, and D. Sorensen, editors, Dimension Reduction of Large-Scale Systems, volume 45 of Lecture Notes in Computational Science and Engineering, pages 49–82. Springer-Verlag, Berlin/Heidelberg, Germany, 2005. [17] S. Gugercin, D. C. Sorensen, and A. C. Antoulas. A modified low-rank Smith method for large-scale Lyapunov equations. Numer. Algorithms, 32(1):27–55, 2003. [18] M. Hanke-Bourgeois. Grundlagen der numerischen Mathematik und des wissenschaftlichen Rechnens. Teubner BG GmbH, 2006. [19] I. Jaimoukha and E. Kasenally. Krylov subspace methods for solving large Lyapunov equations. SIAM J. Numer. Anal., 31:227–251, 1994. [20] A. Klimke. Sparse Grid Interpolation Toolbox – user’s guide. Technical Report IANS report 2007/017, University of Stuttgart, 2007.

7

[21] A. Klimke and B. Wohlmuth. Algorithm 847: spinterp: Piecewise multilinear hierarchical sparse grid interpolation in MATLAB. ACM Transactions on Mathematical Software, 31(4), 2005. [22] J.-R. Li and J. White. Low rank solution of Lyapunov equations. SIAM J. Matrix Anal. Appl., 24(1):260–280, 2002. [23] B. C. Moore. Principal component analysis in linear systems: controllability, observability, and model reduction. IEEE Trans. Automat. Control, AC-26(1):17–32, 1981. [24] C. Moosmann and A. Greiner. Convective thermal flow problems. In P. Benner, V. Mehrmann, and D. Sorensen, editors, Dimension Reduction of Large-Scale Systems, volume 45 of Lecture Notes in Computational Science and Engineering, pages 341–343. Springer-Verlag, Berlin/Heidelberg, Germany, 2005. [25] C. Moosmann and J. Korvink. Automatic parametric MOR for MEMS design. In Tagungsband GMAFachausschuss 1.30, Modellbildung, Identifizierung und Simulation in der Automatisierungstechnik”, Workshop am Bostalsee, pages 89–99, 2006. [26] T. Penzl. A cyclic low-rank Smith method for large sparse Lyapunov equations. SIAM J. Sci. Comput., 21(4):1401–1418, 1999/00. [27] Y. Saad. Numerical solution of large Lyapunov equations. In M. A. Kaashoek, J. H. van Schuppen, and A. C. M. Ran, editors, Signal Processing, Scattering, Operator Theory and Numerical Methods, pages 503– 511. Birkh¨ auser, 1990. [28] V. Simoncini. A new iterative method for solving largescale Lyapunov matrix equations. SIAM J. Sci. Comput., 29(3):1268–1288, 2007. [29] M. S. Tombs and I. Postlethwaite. Truncated balanced realization of a stable nonminimal state-space system. Internat. J. Control, 46(4):1319–1330, 1987. [30] D. S. Weile, E. Michielssen, E. Grimme, and K. Gallivan. A method for generating rational interpolant reduced order models of two-parameter linear systems. Appl. Math. Lett., 12(5):93–102, 1999. [31] H. Yserentant. On the multilevel splitting of finite element spaces. Numer. Math., 49(4):379–412, 1986. [32] C. Zenger. Sparse grids. In Parallel algorithms for partial differential equations (Kiel, 1990), volume 31 of Notes Numer. Fluid Mech., pages 241–251. Vieweg, Braunschweig, 1991. Manuskripteingang: 09. Januar 2009.

8

1 2

Dr.rer.nat Ulrike Baur ist wissenschaftliche Mitarbeiterin in der Professur f¨ ur Mathematik in Industrie und Technik in der Fakult¨ at f¨ ur Mathematik der Technischen Universit¨ at Chemnitz. Hauptarbeitsgebiete: Numerische Methoden f¨ ur Modellreduktion, hierarchische Matrizen. Adresse: Technische Universit¨ at Chemnitz, Fakult¨ at f¨ ur Mathematik, 09107 Chemnitz, Fax: + 49-(0)371-531-22509, E-Mail: [email protected] Prof. Dr.rer.nat.habil. Peter Benner Inhaber der Professur Mathematik in Industrie und Technik in der Fakult¨ at f¨ ur Mathematik der Technischen Universit¨ at Chemnitz. Hauptarbeitsgebiete: Numerische Methoden f¨ ur Modellreduktion, Regelung und optimale Steuerung dynamischer Systeme sowie Numerische Lineare Algebra und Parallele Algorithmen. Adresse: siehe oben, E-Mail: [email protected]