Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

12
Durch den Einsatz von numerischen Modellen für ingenieurtech- nische Problemstellungen, wie z. B. der FEM oder der FDM, kön- nen zunehmend komplexere Berechnungen in immer kürzerer Zeit bewältigt werden. Gleichzeitig ergibt sich jedoch bei dem Einsatz dieser Werkzeuge der Bedarf an Werten für die verschie- denen Modellparameter, von rein konstitutiven Kennwerten bis hin zu geometrischen Angaben, für deren Bestimmung zunehmend inverse Verfahren Anwendung finden. Bei der Nutzung dieser Methoden ist jedoch – insbesondere bei komplizierten Simulatio- nen – mit sehr langen Berechnungszeiten zu rechnen. Gegen- stand dieses Beitrags ist die Vorstellung einer Verfahrensklasse, die eine Abschätzung der Lösung solcher inverser Aufgaben auf der Basis von relativ wenigen Stützstellen ermöglicht. An die Ver- teilung der Stützstellen werden geringste Anforderungen gestellt, so daß diese wahlweise aus vorhergehenden Simulationen oder auch aus alternativen Quellen stammen können. Im Rahmen dieses Beitrags soll ausgehend von einer Einführung in den theoretischen Ansatz eine Strategie zur Beschleunigung der Lösung von inversen Problemstellungen und Optimierungs- aufgaben an einem Beispiel aus dem Gebiet der Geomechanik vorgestellt werden. Effective algorithm for solving inverse problems – geomechani- cal application. When working with numerical models, it is es- sential to determine model parameters which are as realistic as possible. Optimization techniques are being employed more and more frequently for solving this task. However, using these me- thods may lead to very high time costs – in particular, if rather complicated forward calculations are involved. In this paper, we present a class of methods that allows estimat- ing the solution of this kind of optimization problems based on relatively few sampling points. We put very weak constraints on the sampling point distribution; hence, they may be taken from previous forward calculations as well as from alternative sour- ces. Starting from an introduction into the theoretical approach, a strategy for speeding up inverse optimization problems is intro- duced which is illustrated by an example geomechanics. 1 Motivation Zur Gewinnung von Kennwerten für Modellparameter stehen verschiedene Ansätze zur Verfügung. Zum einen können diese durch Messungen aus Feld- und Laborver- suchen gewonnen werden, was jedoch einen erheblichen finanziellen und zeitlichen Aufwand bedeutet. Zum zwei- 470 © Ernst & Sohn Verlag für Architektur und technische Wissenschaften GmbH & Co. KG, Berlin · Bautechnik 83 (2006), Heft 7 ten kann auf ggf. vorhandene „Erfahrungswerte“ und Werte aus der Literatur zurückgegriffen werden, die jedoch oft- mals als zu ungenau anzusehen sind. Als dritte Möglich- keit bietet sich die indirekte Ermittlung durch inverse Ver- fahren an. Diese Optimierungsverfahren versuchen, durch die iterative Anpassung von Modellparametern einer als „Vorwärtsrechnung“ bezeichneten Simulation eine mög- lichst gute Übereinstimmung der Berechnungsergebnisse mit Meßwerten zu erlangen. Der sich über Fachgebiets- grenzen immer weiter verbreitende Einsatz der inversen Methoden zur Bestimmung von Modellparametern spie- gelt sich auch in der verfügbaren Fachliteratur wider: Es finden sich eine Vielzahl von Beiträgen beispielsweise für die Themenbereiche Strukturanalyse (z. B. Matous ˇ et al. 2000 [19]), Strömungsmechanik (z. B. Jeong 2003 [14]), Hydrogeologie (z. B. Carrera et al. 2005 [6]), Maschinen- und Automobilbau (z. B. Fleischer & Broos 2004 [12], Flores Santiago & Bausinger 1998 [13]), Geotechnik (z. B. Schanz et al. 2006 [24], Cui & Sheng 2006 [8]) und viele weitere. Für diese sukzessive Anpassung der Eingangswerte ist jedoch eine zum Teil enorme Anzahl von Vorwärtsrech- nungen nötig, die sich – insbesondere bei rechenaufwen- digen Modellen – in entsprechenden Berechnungszeiten niederschlägt. Zwar kann hier durch Beschleunigung und Parallelisierung der Vorwärtsrechnung eine Minderung des Zeitbedarfs erreicht werden, jedoch sind hierdurch nur Verbesserungen in beschränkten Grenzen möglich. So würden beispielsweise für die Anschaffung und Wartung von Rechentechnik für eine zusätzliche Parallelisierung nicht unerhebliche finanzielle Mittel benötigt. Mit dem Wissen, daß der erforderliche Berechnungs- aufwand hauptsächlich von der Anzahl der aufgerufenen Vorwärtsrechnungen abhängt, lassen sich folgende Ziele festsetzen: Verringerung der Anzahl der benötigten Vorwärts- rechnungen: Durch die Nutzung immer komplizierterer Materialmodelle und komplexerer Geometrien steigt der Zeitbedarf der „einfachen Vorwärtsrechnung“ stark an. Je- der eingesparte Aufruf bedeutet somit einen Zeitgewinn für eine Optimierung. Einbeziehung aller zur Verfügung stehenden Informa- tionen: Die meisten Optimierungsverfahren nutzen aus- schließlich Informationen, die durch sie selbst berechnet wurden. Ergebnisse anderer Verfahren, wie Parametervek- toren und deren Ergebnisse der Vorwärtsrechnung bzw. Effektiver Algorithmus zur Lösung von inversen Auf- gabenstellungen – Anwendung in der Geomechanik Jörg Meier Sebastian Rudolph Tom Schanz Fachthemen DOI: 10.1002/bate.200610040

Transcript of Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Durch den Einsatz von numerischen Modellen für ingenieurtech-nische Problemstellungen, wie z. B. der FEM oder der FDM, kön-nen zunehmend komplexere Berechnungen in immer kürzererZeit bewältigt werden. Gleichzeitig ergibt sich jedoch bei demEinsatz dieser Werkzeuge der Bedarf an Werten für die verschie-denen Modellparameter, von rein konstitutiven Kennwerten bishin zu geometrischen Angaben, für deren Bestimmung zunehmendinverse Verfahren Anwendung finden. Bei der Nutzung dieserMethoden ist jedoch – insbesondere bei komplizierten Simulatio-nen – mit sehr langen Berechnungszeiten zu rechnen. Gegen-stand dieses Beitrags ist die Vorstellung einer Verfahrensklasse,die eine Abschätzung der Lösung solcher inverser Aufgaben aufder Basis von relativ wenigen Stützstellen ermöglicht. An die Ver-teilung der Stützstellen werden geringste Anforderungen gestellt,so daß diese wahlweise aus vorhergehenden Simulationen oderauch aus alternativen Quellen stammen können.Im Rahmen dieses Beitrags soll ausgehend von einer Einführungin den theoretischen Ansatz eine Strategie zur Beschleunigungder Lösung von inversen Problemstellungen und Optimierungs-aufgaben an einem Beispiel aus dem Gebiet der Geomechanikvorgestellt werden.

Effective algorithm for solving inverse problems – geomechani-cal application. When working with numerical models, it is es-sential to determine model parameters which are as realistic aspossible. Optimization techniques are being employed more andmore frequently for solving this task. However, using these me-thods may lead to very high time costs – in particular, if rathercomplicated forward calculations are involved.In this paper, we present a class of methods that allows estimat-ing the solution of this kind of optimization problems based onrelatively few sampling points. We put very weak constraints onthe sampling point distribution; hence, they may be taken fromprevious forward calculations as well as from alternative sour-ces.Starting from an introduction into the theoretical approach, astrategy for speeding up inverse optimization problems is intro-duced which is illustrated by an example geomechanics.

1 Motivation

Zur Gewinnung von Kennwerten für Modellparameterstehen verschiedene Ansätze zur Verfügung. Zum einenkönnen diese durch Messungen aus Feld- und Laborver-suchen gewonnen werden, was jedoch einen erheblichenfinanziellen und zeitlichen Aufwand bedeutet. Zum zwei-

470 © Ernst & Sohn Verlag für Architektur und technische Wissenschaften GmbH & Co. KG, Berlin · Bautechnik 83 (2006), Heft 7

ten kann auf ggf. vorhandene „Erfahrungswerte“ und Werteaus der Literatur zurückgegriffen werden, die jedoch oft-mals als zu ungenau anzusehen sind. Als dritte Möglich-keit bietet sich die indirekte Ermittlung durch inverse Ver-fahren an. Diese Optimierungsverfahren versuchen, durchdie iterative Anpassung von Modellparametern einer als„Vorwärtsrechnung“ bezeichneten Simulation eine mög-lichst gute Übereinstimmung der Berechnungsergebnissemit Meßwerten zu erlangen. Der sich über Fachgebiets-grenzen immer weiter verbreitende Einsatz der inversenMethoden zur Bestimmung von Modellparametern spie-gelt sich auch in der verfügbaren Fachliteratur wider: Esfinden sich eine Vielzahl von Beiträgen beispielsweise fürdie Themenbereiche Strukturanalyse (z. B. Matous et al.2000 [19]), Strömungsmechanik (z. B. Jeong 2003 [14]),Hydrogeologie (z. B. Carrera et al. 2005 [6]), Maschinen-und Automobilbau (z. B. Fleischer & Broos 2004 [12],Flores Santiago & Bausinger 1998 [13]), Geotechnik (z. B.Schanz et al. 2006 [24], Cui & Sheng 2006 [8]) und vieleweitere.

Für diese sukzessive Anpassung der Eingangswerte istjedoch eine zum Teil enorme Anzahl von Vorwärtsrech-nungen nötig, die sich – insbesondere bei rechenaufwen-digen Modellen – in entsprechenden Berechnungszeitenniederschlägt. Zwar kann hier durch Beschleunigung undParallelisierung der Vorwärtsrechnung eine Minderungdes Zeitbedarfs erreicht werden, jedoch sind hierdurchnurVerbesserungen in beschränkten Grenzen möglich. Sowürden beispielsweise für die Anschaffung und Wartungvon Rechentechnik für eine zusätzliche Parallelisierungnicht unerhebliche finanzielle Mittel benötigt.

Mit dem Wissen, daß der erforderliche Berechnungs-aufwand hauptsächlich von der Anzahl der aufgerufenenVorwärtsrechnungen abhängt, lassen sich folgende Zielefestsetzen:– Verringerung der Anzahl der benötigten Vorwärts-rechnungen: Durch die Nutzung immer kompliziertererMaterialmodelle und komplexerer Geometrien steigt derZeitbedarf der „einfachen Vorwärtsrechnung“ stark an. Je-der eingesparte Aufruf bedeutet somit einen Zeitgewinnfür eine Optimierung.– Einbeziehung aller zur Verfügung stehenden Informa-tionen: Die meisten Optimierungsverfahren nutzen aus-schließlich Informationen, die durch sie selbst berechnetwurden. Ergebnisse andererVerfahren, wie Parametervek-toren und deren Ergebnisse der Vorwärtsrechnung bzw.

Effektiver Algorithmus zur Lösung von inversen Auf-gabenstellungen – Anwendung in der Geomechanik

Jörg MeierSebastian RudolphTom Schanz

Fachthemen

DOI: 10.1002/bate.200610040

471Bautechnik 83 (2006), Heft 7

Zielfunktionswerte, werden entweder nicht oder nur rudi-mentär einbezogen. Beispielsweise gibt ein Gradienten-verfahren die Werte aller Parameter für jeden Schritt vor –und läßt keine „fremden“ Vektoren zu. Mit dieser Zielset-zung ist ebenfalls die Forderung nach einer möglichstgroßen Unempfindlichkeit gegenüber der Lage dieser Para-metervektoren verknüpft. So wäre beispielsweise die Be-schränkung auf rein rasterförmig angeordnete Daten sehrunvorteilhaft.– Nutzung der Vorteile bekannter Verfahren: Ziel solltees nicht sein, einen neuen Optimierungsalgorithmus zu er-stellen. Es sollten vielmehr die Vorteile bereits bekannterund etablierter Verfahren genutzt werden. Hierbei bietetsich die Verwendung bereits vorhandener „Schnittstellen“der Optimierungsverfahren an, um auch den Einsatz vonzukünftigen bzw. verbesserten Algorithmen zur Extrem-wertfindung zu erlauben. Weiterhin sind mit diesem An-satz sequentielle, aufeinander aufbauende Läufe verschie-dener Methoden möglich.

2 Grundlagen der Optimierung

Im allgemeinen wird der dem Fachgebiet der angewand-ten Mathematik entstammende Prozeß der iterativen Su-che nach den „günstigsten“ Werten der Parameter einerVorwärtsrechnung („Lösungsvektor“) als Optimierung be-zeichnet. Wie in Bild 1 dargestellt, wird die Güte einesParametersatzes i. d. R mit Hilfe einer Zielfunktion be-stimmt, welche in Abhängigkeit von den Referenzdatenund den Ergebnissen der aktuellen Vorwärtsrechnung einereelle Zahl, den Zielfunktionswert (auch: „Gütewert“), er-rechnet. Einem Parametervektor kann somit eindeutig einZielfunktionswert zugeordnet werden.

Eng verbunden mit der Optimierung ist der Begriffder Sensitivitätsanalyse. Mit Hilfe einer Sensitivitätsana-lyse wird ausgehend von einem vorzugebenden Parameter-vektor die Änderung des Funktionswertes in Relation zuden Modifikationen der Werte der Eingangsparameter ge-setzt. Aufgeschlüsselt nach den Parametern kann so in dernäheren Umgebung des vorgegebenen Parametervektorseine Teilmenge (Subset) von sensitiven Parametern ermit-telt werden. So können mit diesem Werkzeug Parameter-sätze vor einer Optimierung bestimmt und indifferenteoder voneinander abhängige Parameter aus einer Parame-terauswahl entfernt werden. Weiterhin nutzen verschiedeneOptimierungsalgorithmen Techniken der Sensitivitätsana-lyse (siehe Schwarz 2001 [26]).

Durch den für jeden der n Parameter der zu optimie-renden Vorwärtsrechnung vorgegebenen Definitionsbereichwird ein n-dimensionaler Suchbereich W aufgespannt. DieSuche nach dem Optimum kann somit als Extremwert-suche innerhalb dieses Gebietes umschrieben werden. Die

Gesamtheit aller Zielfunktionswerte innerhalb des Such-bereichs wird als Zielfunktionstopologie bzw. Zielfunk-tionshyperfläche bezeichnet. Bild 2 zeigt schematisch dieZielfunktionshyperfläche in Abhängigkeit zweier zu iden-tifizierender Eingangsgrößen (Parameter x1, Parameter x2).

Aus mathematischer Sicht kann von einem lokalenExtremwert gesprochen werden, wenn für den Lösungs-vektor x* in der Umgebung U der Zielfunktion f : W Æ RGl. (1) bzw. (2) gilt. Entspricht der Suchbereich W der Um-gebung U, so handelt es sich bei x* um einen globalen Ex-tremwert. Gleiches gilt, falls der Suchbereich x* enthältund vollständig in U liegt (Boyd & Vandenberghe 2006[4], Connor 1976 [7], Press et al. 1992 [23])

f(x*) £ f(x), "x ŒU für ein lokales Minimum von f in U (1)

f(x*) ≥ f(x), "x ŒU für ein lokales Maximum von f in U (2)

Hierbei kann jedoch eine Untermenge X* von U den glei-chen Funktionswert f(x*) aufweisen. Somit liegt für allef(x) = f(x*) keine eindeutige Lösung vor. Von einem ein-deutigen Extremwert kann nur gesprochen werden, wennder Funktionswert f(x*) ausschließlich für den Lösungs-vektor x* zu finden ist.

Allgemein ist die Existenz eines Optimums oder eineroptimalen Lösungsmenge dann gesichert, wenn W endlichoder f(x) stetig und ein W kompakter topologischer Raumist (Analysis: Satz von Heine-Borel). Von einem endlichenW kann beispielsweise dann ausgegangen werden, wennalle zu bestimmenden Parameter ganzzahlig oder booleschsind. Da für die meisten in der Praxis vorkommenden Op-timierungsaufgaben weder die Endlichkeit von W noch dieStetigkeit von f(x) vorauszusetzen ist, kann nicht von einergrundsätzlichen Existenz eines Optimums ausgegangenwerden. In der Praxis kann jedoch oft von bereits gelöstenAufgabenstellungen auf aktuelle Probleme geschlossen unddieses Wissen durch numerische Experimente abgesichertwerden.

Erschwerend für eine Optimierung kann eine Ziel-funktion durch Rauhigkeiten überlagert werden. Quellenfür solche Störgrößen können beispielsweise numerische

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 1. Schema der Berechnung eines ZielfunktionswertesFig. 1. Calculation scheme for determination an objectivefunction

Bild 2. Schema der Darstellung einer Zielfunktionstopologiefür eine zweidimensionale AufgabenstellungFig. 2. Visualization of objective function topology for atwo-dimensional problem

Ungenauigkeiten der Vorwärtsrechnung, Rundungsfehlerbei der Übergabe der Werte zwischen Simulation und Op-timierungsalgorithmus als auch die Zielfunktion selbst sein.Suchalgorithmen wie auch Zielfunktionen, die ohne einRauschen hervorragende Ergebnisse liefern, können durchdiesen Effekt vollständig unbrauchbar werden. Scheinbarelokale Extremwerte und fehlerhafte Sensitivitäten bzw.Gradienten sind häufige Folgeerscheinungen (Miller 1997[21], Polheim 1999 [22], Bui et al. 2005 [5]).

Das grundlegende Ablaufschema der Optimierungs-algorithmen kann Bild 3 entnommen werden. Nach derVorgabe von Startparametern ruft der Optimierungsalgo-rithmus die Vorwärtsrechnung ein oder mehrere Male aufund extrahiert die relevanten Daten nachfolgend. Durchden Vergleich mit vorgegebenen Soll- bzw. Referenzwer-ten durch die Zielfunktion wird der zugehörige Zielfunk-tionswert errechnet. Auf der Basis dieses Gütewertes wirddurch den Algorithmus ein neuer Satz von Berechnungs-parametern festgelegt und ein weiterer Optimierungszy-klus begonnen. Alternativ kann bei der Erfüllung eines Ab-bruchkriteriums die Schleife verlassen werden.

Die verschiedenen aus der angewandten Mathematikbekannten Optimierungsalgorithmen lassen sich in folgendeGruppen einteilen:– Einfache Verfahren: Zu den einfachen Methoden zäh-len Verfahren, die zur Findung des gesuchten Extremwer-tes die bisher bestimmten Lösungen der Abweichungsfunk-tion nicht einbeziehen. Beispiele sind die Rasterverfahren,das Monte-Carlo-Verfahren und das Latin-Hypercube-Ver-fahren (McKay et al. 1979 [20], Kollig & Keller 2002 [16]).Bei diesen Verfahren wird für eine vorher festgelegte An-zahl von Parameterkombinationen die Abweichung be-stimmt und die „beste“ Kombination als Ergebnis geliefert.Diese Methoden eignen sich nur bei einer geringen Anzahlvon Parametern.– Gradientenbasierende Verfahren: Bei den Gradienten-methoden wird versucht, für die aktuelle Parameterkombi-nation den Änderungsvektor zu finden, der zu einer mög-lichst „guten“, d. h. geringen Abweichung in der näherenUmgebung führt. Diese Methoden zeichnen sich beson-ders bei einer großen Anzahl von Parametern (in der Näheder globalen Lösung) durch eine geringere Anzahl von

472 Bautechnik 83 (2006), Heft 7

Aufrufen aus. Nachteile zeigen diese Verfahren bei „rau-hen“ Zielfunktionen, da sie ggf. in lokalen Senken „hän-gen bleiben“ können.– Populationsbasierende Verfahren: Populationsbasie-rende Verfahren nehmen Anleihen in der Natur. DieseAlgorithmen nutzen eine definierte Anzahl von unabhän-gigen Parametersätzen („Individuen“ bzw. „Partikel“) in-nerhalb eines Optimierungszyklus, die i. d. R. einmal proZyklus Informationen gemäß einem vorgegebenen Regel-satz austauschen können. Beispiele für diese Verfahrens-klasse sind evolutionäre Algorithmen, genetische Verfah-ren und Partikel-Schwarm-Methoden. Vorteile dieser Ver-fahrensklasse sind gewöhnlich gute Parallelisierbarkeitund hohe Robustheit gegenüber Rauhigkeiten und fehlge-schlagenen Vorwärtsrechnungen.

3 Grundlagen der globalen Hyperflächenapproximations-verfahren

Die hier vorgestellten Hyperflächenapproximationsverfah-ren wurden originär von den Autoren entwickelt. Wie be-reits oben dargestellt, ist die Minimierung der Anzahl derAufrufe der Vorwärtsrechnung für eine performante Opti-mierung unerläßlich. Generalisierend kann die Aussagegetroffen werden, daß Optimierungsverfahren i. d. R. Vor-wärtsrechnungen „nur“ für die Bestimmung von weiterenStützpunkten der Zielfunktionshyperfläche verwenden.

Grundidee der Hyperflächenapproximationsverfah-ren ist die Ersetzung eben jener aufwendig zu berechnen-den Zielfunktion durch eine deutlich weniger recheninten-sive Näherungsfunktion. Die Suche nach dem Extremumerfolgt dann, indem die Näherungsfunktion anstelle dereigentlichen Vorwärtsrechnung aufgerufen wird. Mit der soermittelten Extremstelle wird durch einmalige Vorwärts-rechnung eine neue Stützstelle und mit Hilfe dieser eineexaktere Näherungsfunktion bestimmt.

Die Ersetzung der Zielfunktionshyperfläche kannhierbei jedoch nicht – wie in Bild 4a dargestellt – durch li-neare Interpolation erfolgen, da durch diese keine Aus-sagen über evtl. Extremwerte zwischen den Stützstellenmöglich sind. Für die Verfahrensklasse der globalen Hyper-flächenapproximationsverfahren sind somit höherwertigeAnsätze (Bild 4b) nötig.

Aus mathematischer Sicht besteht die Aufgabe alsodarin, die Zielfunktion f : W Æ R, die zu jedem Parameter-vektor x Œ W einen reellen „Gütewert“ f(x) liefert, auf derGrundlage von Stützpunkten, d. h. bereits bekanntenFunktionswerten yi = f(xi) mit i = 1 … k durch eine „Er-satzfunktion“ f : W Æ R zu approximieren. Der hier vorge-stellte Ansatz für die Gewinnung von f besteht darin,

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 3. Grundlegendes Ablaufschema der OptimierungFig. 3. Principal flowchart of optimization procedure

Bild 4. Schema der Approximation von Stützpunkten durcha) Hyperebenen und b) HyperflächenFig. 4. Visualization of approximation based on samplingpoints by a) hyperplanes and b) hypersurfaces

a) b)Parameter Parameter

Abw

eich

ung

Abw

eich

ung

473Bautechnik 83 (2006), Heft 7

zunächst für jede Stützstelle xi eine lokale Approxima-tionsfunktion f i : W Æ R zu bestimmen, welche die Funk-tion f in der Umgebung von xi annähert. Ein Funktions-wert f (x) der globalen Approximationsfunktion wird dannals gewichtetes Mittel der einzelnen Funktionswerte f i(x)bestimmt:

(3)

Hierbei ist wi : W Æ R+ eine Wichtungsfunktion, welche je-dem Vektor des Suchraums eine nichtnegative reelle Zahlzuordnet. Dabei gibt wi(x) an, wie stark die Funktion f i ander Stelle x zur globalen Approximationsfunktion f bei-trägt. Bild 5 visualisiert das zugrunde liegende Prinzip desAnsatzes.

Die notwendigen Kriterien für die zu verwendendenFunktionen f i und wi lassen sich wie folgt spezifizieren:– Übereinstimmung der lokalen Approximationsfunktio-nen mit den zugeordneten Stützstellen: f i(xi) = f(xi). DieseForderung ergibt sich aus der zugrunde liegenden Idee,daß f i in der Umgebung von xi die Zielfunktion annähert.– Mit zunehmender radialer Entfernung vom zugeordne-ten Stützpunkt xi fällt die Wichtungsfunktion wi monoton,d. h., r1 < r2 impliziert wi(xi + r1v) ≥ wi(xi + r2v) für alle po-sitiven reellen Zahlen r1 und r2 sowie beliebige Richtungs-vektoren u Œ W im Suchbereich.– Zur Gewährleistung der Übereinstimmung der globalenApproximationsfunktion mit allen Stützstellen (also f(xi) =f(xi) für alle xi) können zwei alternative Kriterien formu-liert werden:∑ Eine Polstelle der Funktion wi an der zugeordneten Stütz-stelle xi:

Unter dieser Bedingung verschwinden die Beiträge alleranderen lokalen Approximationsfunktionen für x Æ xi. Diein diesem Fall an den Stützstellen entstehenden Lücken inder Funktion f können durch Setzung von f(xi) = f(xi) ver-mieden werden. Überdies bleibt die Funktion stetig, fallsalle f i auf W und alle wi auf W\{xi} stetig sind.∑ Vollständiges Abklingen der Wichtungsfunktionen:wi(xj) = 0 für alle i π j aus 1 … k. Vorteil dieses Kriteriums

lim ( ) .x x

ii

w xÆ

= +•

( ):

( ) ˆ ( )

( )

f x

w x f x

w x

i ii

k

ii

k= =

=

Â

Â1

1

ist, daß die Stützstellen nicht gesondert behandelt werdenmüssen. Problematisch hingegen ist in diesem Fall (unterEinhaltung oben genannter Forderungen) sicherzustellen,daß die Funktion f auf ganz W (zumindest aber auf einemsinnvollen Teilbereich T W wie beispielsweise im Inne-ren des durch die Stützstellen aufgespannten Bereichs) de-finiert bleibt. Dafür muß für alle x Œ T die Summe

größer als Null sein, es dürfen also an keiner Stelle alleWichtungsfunktionen gemeinsam verschwinden.∑ Je nach eingesetztem Optimierungsverfahren muß dieglobale Approximationsfunktion weitere Anforderungenhinsichtlich Stetigkeit bzw. Differenzierbarkeit erfüllen.So kann beispielsweise die rechenaufwendige nume-rische Bestimmung der Gradienten bei einer entspre-chend leistungsfähigen Funktion entfallen und direkt er-folgen.

Bei den hier durchgeführten Untersuchungen wurdenPolynomfunktionen für die lokalen Approximationen unddie in der Visualisierung häufig anzutreffenden Smooth-step-Funktionen (Krüger 2002 [17]) für die Wichtung ver-wendet.

An dieser Stelle wird der Unterschied zu den Response-Surface-Methoden (RSM; Biles 1975 [3], Daugherty &Turnquist 1980 [9], Allen & Yu 2002 [1], Wang 2003 [29]u. v. a.) deutlich: Die RSM versuchen, die Zielfunktiondurch ein vorgegebenes einfaches mathematisches Modellauf der Basis von Stützstellen lokal anzunähern. Die zu-grunde liegenden Modelle werden i. d. R. so gewählt, daßihre Extremwerte direkt aus den Parametern bestimmtwerden können. Zum Einsatz kommen hier oft quadrati-sche oder kubische Ansätze. Durch den iterativen Einsatzdieses lokalen Verfahrens findet – für entsprechend gut-artige Zielfunktionen – eine Annäherung an die gesuchteExtremstelle statt.

Beschränkend wirkt sich sehr oft die Tatsache aus,daß durch die vorgegebenen mathematischen Modelleauch die Anzahl der Freiheitsgrade direkt von der Anzahlder Modellparameter – und nicht von der Anzahl der zurVerfügung stehenden Stützstellen – abhängt. Es wirddemnach eine Mindestanzahl von Stützstellen gefordert.Parametervektoren und deren Abweichungen, die überdiese Forderung hinaus verwendet werden, verbessernnur die Robustheit gegenüber Rauhigkeiten der Zielfunk-tion – liefern aber keine weiterführende Annäherung bzw.komplexere Response Surface. Von der Wahl polynoma-ler mathematischer Modelle höheren als dritten Gradeswird in der Regel abgesehen, da das Verhalten zwischenden Stützstellen nicht mehr zuverlässig kontrolliert wer-den kann beziehungsweise eine direkte Berechnung derExtremstellen aus den Parametern dann nicht mehr mög-lich wäre.

Diese Nachteile der RSM werden durch die Hyper-flächenapproximationsmethoden nicht geteilt. Es bestehtweder die Forderung nach einer festgelegten Anzahl vonStützstellen, noch wird die Komplexität der Approxima-tion durch das mathematische Modell begrenzt. Überdieswird die Zielfunktion global angenähert, was beispiels-weise das Finden mehrerer separater Extrema in ein undderselben Approximationsfunktion ermöglicht.

w xii

k

( )=

Â1

Õ

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 5. Schematische Darstellung des zugrunde liegendenPrinzips der HyperflächenapproximationsverfahrenFig. 5. Basic idea of the hypersurface approximationmethod

w1 w2

f 2

f

f 1

f(x1)

f(x2)Wer

te d

er W

icht

ungs

- un

dA

ppro

xim

atio

nsfu

nktio

nen

Bereich der Dominanz von f 1 Bereich der Dominanz von f 2 x

4 Anwendungsbeispiele4.1 Beispiel 1: Mathematische Testfunktionen

Um die Leistungsfähigkeit und Eigenschaften von Optimie-rungsalgorithmen zu testen, wurden von verschiedenenAutoren Standardtestfunktionen vorgeschlagen, die einezu optimierende Zielfunktion explizit beschreiben. Der Ein-satz dieser Testfunktionen zur Validation und Verifikationder Hyperflächenapproximationsverfahren bietet verschie-dene Vorteile. Zum einen sind die Zielfunktionswerte die-ser Funktionen sehr schnell zu ermitteln. Weiterhin kön-nen ungünstige Einflüsse aus ggf. extern aufzurufendenund auszulesenden Vorwärtsrechnungen ausgeschlossenwerden.

Im Rahmen des Beispiels 1 soll die Anwendung derHyperflächenapproximationsverfahren auf folgende Aus-wahl von Standardtestfunktionen gezeigt werden:– Rastrigins Funktion– Rosenbrocks Funktion– „Moved Axis Parallel Hyper-Ellipsoid“-Funktion

Rastrigins FunktionEine häufig verwendete Testfunktion ist die in Gl. (4) dar-gestellte Rastrigin-Funktion. Diese Funktion ist an De JongsTest-Suite-Funktion 1 angelehnt (De Jong 1975 [10]) undbasiert auf einer quadratischen Funktion, zu der eine Cosi-nus-Modulation addiert wird. Durch diese Konfigurationsind die Minima regelmäßig verteilt. Das globale Minimumbefindet sich bei F(x*) = 0 mit x*i = 0 und i = 1 … n. Wei-terhin kann für diese Funktion die Anzahl der Dimensio-nen bzw. Parameter beliebig festgesetzt werden:

474 Bautechnik 83 (2006), Heft 7

Bild 6 zeigt im oberen Bereich den originalen Verlauf vonRastrigins Funktion für zwei Parameter im Bereich von je-weils [–1 … 1]. Im unteren Bereich ist die Visualisierungder approximierten Zielfunktionstopologie auf der Basiseines 6 ¥ 6-Rasters dargestellt. Das Raster weist als mini-male Zielfunktionswerte die 4 Punkte mit den Koordina-ten (–0,16; –0,16), (0,16; –0,16), (–0,16; 0,16) und (0,16;0,16) aus. Die Approximation hingegen ist in der Lage, dieoriginale Form von Rastrigins Funktion bereits mit diesen36 Stützstellen sehr gut nachzubilden. Nachteilig könnensich auf eine nachgeschaltete Optimierung die Randberei-che der Hyperfläche auswirken. Diese zeigen aufgrundder fehlenden Information außerhalb des Suchbereichsdes zugrunde liegenden Rasters größere Abweichungenvon der „originalen“ Funktion und täuschen nicht vorhan-dene Extremwerte vor. Als Gegenmaßnahme empfiehlt essich, den Suchbereich nachgeschalteter Optimierungenum 5 bis 10 % einzuschränken.

Bild 7a zeigt die Ergebnisse einer Optimierung auf derApproximationshyperfläche mit einem Particle-Swarm-Optimizer (Eberhardt & Kennedy 1995 [11], Kennedy &Eberhardt 1995 [15], Shi & Eberhardt 1998 [27], [28], vanden Berg 2001 [2]) mit 10 Individuen. Nach bereits 4 Be-rechnungszyklen (entspricht 40 Funktionsaufrufen) weistderAlgorithmus als günstigsten Lösungsvektor (5,928E-003;1,292E-002) aus. Ein ähnliches Bild zeigt sich bei der inBild 7b dargestellten Gradientenmethode (modifiziertes

F x n x xi ii

n

( ): cos( )= ◊ + - ◊ ◊ ◊( )=

Â10 10 22

1

p

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 6. Originalverlauf (oben) und zwei Ansichten der Annäherung von Rastrigins Funktion für zwei Parameter auf derBasis eines 6 ¥ 6-Rasters (grün: 6 ¥ 6-Raster; orange: Approximation)Fig. 6. Original shape (above) and two perspective views of the approximation of Rastrigin’s function for two parameters onthe basis of a 6 ¥ 6-raster (green: 6 ¥ 6-raster; orange: approximation)

475Bautechnik 83 (2006), Heft 7

und gedämpftes Newton-Raphson-Verfahren; Press et al.1992 [23]): nach bereits 8 Berechnungsschritten (entspricht40 Funktionsaufrufen) erreicht dieses Verfahren ausgehendvon dem Startvektor (0,4; 0,4) einen Lösungsvektor von(4,250E-004; 4,241E-004).

Rosenbrocks FunktionRosenbrocks Funktion – auch als „Bananenfunktion“ be-kannt – war ursprünglich nur für zwei Dimensionen defi-niert worden, läßt sich aber wie in Gl. (5) dargestellt leichtauf eine beliebige Anzahl von Parametern erweitern (Pol-heim 1999 [22]). Ihr globales Minimum weist die Funktionbei F(x*) = 0 mit x*i = 1 und i = 1 … n auf. Bild 8 zeigt dieFunktion für zwei Parameter im Bereich [–2 … 2].

(5)

Die Approximation von Rosenbrocks Funktion für zweiParameter im Bereich [–2 … 2] auf der Basis eines 6 ¥ 6-Rasters und eines 7 ¥ 7-Rasters ist in Bild 9 dargestellt.Bei dem 6 ¥ 6-Raster wird zufälligerweise das globale Mi-nimum bei (1; 1) getroffen und somit auch durch die Ap-proximation wiedergegeben. Die Approximation gibt je-doch auch sehr präzise das gesamte bogenförmige „Tal“von Rosenbrocks Funktion wieder, so daß der Suchbe-reich auch ohne das Wissen um die Lage des globalenOptimums stark eingeschränkt wird. Die „Ungleich-mäßigkeiten“ des wiedergegebenen Tals sind wiederumauf die geringe Zahl der Stützstellen zurückzuführen. Umden „Glückstreffer“ des 6 ¥ 6-Rasters zu umgehen, sind inBild 9b weiterhin die Ergebnisse für ein 7 ¥ 7-Raster dar-gestellt. Durch dieses Raster wird der beste Lösungsvek-tor mit (0; 0) angegeben. Die Visualisierung der Approxi-mationshyperfläche zeigt ein zu Bild 9a ähnliches Bild:Es wird die bogenförmige Ausprägung des Tals gut wie-dergegeben, und es werden Optima nur innerhalb diesesTals ausgewiesen.

F x x x xi ii

n

i( ): = ◊ -( ) + -( )+=

-

 100 112 2

1

12

Bei der Anwendung von Optimierungsverfahren aufdie so erstellten Approximationshyperflächen werden beimehrfachen Durchläufen mit variierenden Startwerten diedrei in Bild 9 unten erkennbaren Optima gefunden. Da derzeitliche und rechentechnische Aufwand der Optimierungder Ersatzflächen z. B. gegenüber einer aufwendigen FEM-Simulation klein ist, sind auch mehrfache Durchläufe pro-blemlos möglich. Die gefundenen Lösungsvektoren kön-nen wiederum an der zugrunde liegenden Vorwärtsrech-nung geprüft werden und ggf. als neue Startwerte dienen.

Moved Axis Parallel Hyper-Ellipsoid-FunktionDie weniger bekannte Moved Axis Parallel Hyper-Ellipsoid-Funktion eignet sich besonders für Demonstrationszwecke,

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 7. Ergebnisse der Optimierung auf den Approximationshyperflächen von Rastrigins Funktion für zwei Parameter aufder Basis eines 6 ¥ 6-Rasters (links: Particle-Swarm-Optimizer; rechts: Gradientenverfahren)Fig. 7. Results of the optimization based on the approximated hyper surfaces of Rastrigin’s function for two parameters onthe basis of a 6 ¥ 6-raster (left: Particle-Swarm-Optimizer; right: gradient method)

Bild 8. Rosenbrocks Funktion für zwei Parameter im Be-reich [–2 … 2]Fig. 8. Rosenbrock’s function with two parameters in therange [–2 … 2]

476 Bautechnik 83 (2006), Heft 7

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 9. Approximation von Rosenbrocks Funktion für zwei Parameter im Bereich [–2 … 2] auf der Basis eines 6 ¥ 6-Rasters(links) und eines 7 ¥ 7-Rasters (rechts)Fig. 9. Approximation of Rosenbrock’s function for two parameters in the range [–2 … 2] on basis of a 6 ¥ 6-raster (left) anda 7 ¥ 7-raster (right)

da das globale Optimum für jeden Parameter einen ande-ren Wert annimmt (F(x*) = 0 mit x*i = 5i und i = 1 … n). DieGrundlage ist ein zu den Koordinatenachsen parallelesHyperellipsoid (Polheim 1999 [22]). Die Definition derTestfunktion kann Gl. (6) entnommen werden. Bild 10(links) zeigt die Visualisierung der Funktion.

(6)F x i x iii

n

( ): ( )= ◊ -( )=

 52

1

Anhand dieser Testfunktion soll die Leistungsfähigkeit einerAuswahl von Optimierungsverfahren mit den Hyper-flächenapproximationsverfahren für den zweidimensionalenFall verglichen werden. Grundlage dieses Vergleiches sollder geometrische Abstand der besten Lösung der Verfahrenzum globalen Optimum (5; 10) sein. Bild 10 (rechts) zeigtden typischen Verlauf der Abweichung und des Abstandeszum globalen Optimum eines Particle-Swarm-Optimizers(PSO), einer Gradientenmethode (modifiziertes und ge-

477Bautechnik 83 (2006), Heft 7

dämpftes Newton-Raphson-Verfahren (Press et al. 1992[23]) und eines evolutionären Algorithmus (Schilling 2003[25]) über den benötigten Aufrufen der Vorwärtsrechnung.Bei der Interpretation dieser Darstellung ist jedoch zu be-achten, daß der PSO und der evolutionäre Algorithmus je-weils 10 Aufrufe der Vorwärtsrechnung je Schritt bzw. Stepausführen, wogegen das Gradientenverfahren 5 Aufrufebenötigt. Somit hat der PSO mit 90 Aufrufen einen Abstandkleiner 1,0 zum globalen Optimum erreicht. Die Gradien-tenmethode benötigt 50 Aufrufe und der evolutionsbasie-rende Algorithmus 90 Aufrufe. Tabelle 1 kann ein Vergleichdes Berechnungsaufwandes entnommen werden. Weiterhinzeigt sich im Bild 10 (rechts) der Einfluß derZielfunktion: Inder Nähe des Optimums bietet die hierverwendete, einfacheFunktion, eine sehr geringe Auflösung gegenüber einer wei-teren Annäherung an das globale Optimum.

Angesichts der einfachen Zielfunktionstopologie wurdeein 4 ¥ 4-Raster als Basis für die Approximation verwendet.Bei der Suche des Optimums auf der Approximations-hyperfläche weist ein rekursives Rasterverfahren (der Such-bereich wird nach jedem Durchlauf auf die beste Rasterzelleund deren Nachbarn eingegrenzt) nach der dritten Rekur-sionsstufe (4,977; 9,977) und ein PSO mit 10 Individuennach 30 Berechnungsschritten (5,096; 9,972) als günstigstenLösungsvektor aus (dies entspricht einem Abstand von 0,03bzw. 0,1 vom globalen Optimum). Das Hyperflächenappro-ximationsverfahren erreicht somit auf der Basis von nur16 Stützstellen ein besseres Ergebnis, als die „herkömmli-chen“Verfahren mit weitaus mehr Funktionsaufrufen.

Für den Fall der Moved Axis Parallel Hyper-Ellipsoid-Funktion mit drei Parametern wurde als Basis ein 4 ¥ 4 ¥ 4-Raster gewählt, welches die Approximationshyperflächen-verfahren ebenfalls zuverlässig lösen konnten.

Die hier demonstrierte geringe Anzahl von Stützstel-len ist nur bei entsprechend „gutartigen“ Zielfunktionstopo-

logien zielführend. Grundsätzlich gilt, daß mit zunehmen-der Komplexität als auch mit zunehmender Parameterzahlzusätzliche Stützstellen für eine verläßliche Approxima-tion benötigt werden. Gleichzeitig sollten diese Ausgangs-werte möglichst „repräsentativ“ für den zu untersuchendenSuchraum sein. In dem hier vorgestellten Beispiel wurdeaus diesem Grund ein einfaches Rasterverfahren verwen-det, da es den vorgegebenen Bereich gleichmäßig abdeckt.Anwendbar sind aber z. B. auch Latin-Hypercube-Verfah-ren.

Zusammenfassung Beispiel 1Anhand der drei Testfunktionen des Beispiels 1 konnte ge-zeigt werden, daß die Hyperflächenverfahren in der Lagesind, Extremwerte einer Zielfunktionstopologie zu appro-ximieren und somit die Lösung von Optimierungsaufga-

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 10. Moved Axis Parallel Hyper-Ellipsoid-Funktion für zwei Parameter im Bereich [–25 … 25] bzw. [–20 … 30] und derVerlauf verschiedener Optimierungsalgorithmen bei der LösungFig. 10. Moved Axis Parallel Hyper-Ellipsoid function with two parameters in the range [–25 … 25] and [–20 … 30] respec-tively and performance of several optimization algorithms

Tabelle 1. Vergleich der Ergebnisse der Moved Axis ParallelHyper-Ellipsoid-Funktion (Anzahl der Berechnungsschrittebzw. Funktionsaufrufe für das Erreichen eines geometri-schen Abstandes von kleiner 1 zum globalen Optimum)Table 1. Comparison of results for the Moved Axis ParallelHyper-Ellipsoid function (number of calculation steps andfunction calls respectively for obtaining a geometric distanceless than 1 to global optimum)

Verfahren Berechnungs- Funktions-schritte aufrufe

Gradientenverfahren 10 50

evol. Verfahren 9 90

Particle-Swarm-Optimizer 9 90

Hyperflächenapproximations-methode + PSO

~ 16

ben zu beschleunigen. Zusammenfassend lassen sich fol-gende Aussagen ableiten:– Die approximierten Flächen können außerhalb und anden Rändern des durch die Stützstellen abgedeckten Ge-biets Extremwerte aufweisen, die in der zugrunde liegen-den Funktion nicht vorhanden sind. Hinreichend verläß-liche Aussagen über Extrema sind daher nur im Innerendes durch Stützstellen abgedeckten Bereichs möglich. Diesist bei der Wahl der Stützstellen zu beachten.– Durch die Approximation können auf der wahren Ziel-funktionstopologie zusammenhängende Extremwertberei-che als getrennte Gebiete wiedergegeben werden. Da Be-rechnungen auf den angenäherten Flächen vergleichsweiseschnell ausgeführt werden können, sind mehrfache Durch-läufe mit variierenden Anfangsbedingungen günstig, diedann mit höherer Wahrscheinlichkeit die verschiedenenExtremwertbereiche ausweisen.– Die verwendeten Stützstellen sollten den zu untersu-chenden Raum möglichst „repräsentativ“ abdecken undder Parameterzahl bzw. der Funktionstopologie Rechnungtragen. Zur Erzeugung und/oder Ergänzung bereits vor-handener Daten bieten sich z. B. Latin-Hypercube-Verfah-ren und bei geringen Parameterzahlen u. U. auch Raster-methoden an.– Die Überprüfung der Qualität der approximierten Flächekann durch den Vergleich approximierter Wert und Wertder Vorwärtsrechnung an ausgewählten Stellen erfolgen.Hierbei sollten selbstverständlich Vorwärtsrechnungen ge-wählt werden, die nicht zur Generierung der Hyperflächegenutzt wurden. Weiterhin bieten höhergradige Polynom-funktionen i. d. R. bessere Approximationen.

4.2 Beispiel 2: Lage einer Trennfläche in einer Böschung

Als weitere Aufgabenstellung wurde, im Sinn eines numeri-schen Experiments, die Fußentlastung eines Hangs ge-wählt. In Bild 11 sind die geometrischen Parameter des Mo-

478 Bautechnik 83 (2006), Heft 7

dells und die Lage von zwei Inklinometern zur Messung vonVerformungen dargestellt. Bild 12 zeigt die Vernetzung desModells mit den vorgegebenen Randbedingungen mittelsdes FEM-Pakets ABAQUS/Standard. Es wurde vereinfa-chend angenommen, daß ausgewählte Materialkennwertein der oberhalb der Schichtgrenze befindlichen Schicht 1von einem Startwert an der Geländeoberkante linear bis aufeinen Endwert an der Schichtgrenze ansteigen und dann inSchicht 2 konstant bleiben (Bild 11).

Im Rahmen eines ersten „geostatischen“ Berechnungs-schritts wurde der hydrostatische Primärspannungszustandbei einer vollständigen Fixierung des Modells vorgegebenund diese Fixierung in einem weiteren Schritt vollständiggelöst. Nach dem Erreichen des sich darauf einstellendenGleichgewichtszustands wurde die Unterkante des Mo-dells zusätzlich horizontal fixiert, um eine Auflage desMaterials auf einer rauhen Felsoberfläche zu simulieren.Anschließend wurde die Fußschüttung in einem 4. Be-rechnungsschritt entfernt und so die Stützung des Hangsim unteren Bereich aufgehoben.

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 11. Modellaufbau des numerischen Experiments der Fußentlastung eines HangsFig. 11. Model configuration of the numerical experiment: release of the lower part of a slope

(E,j)

Bild 12. Vernetzung des Modells und Randbedingungen fürdie Fußentlastung eines Hangs für den 1. BerechnungsschrittFig. 12. Discretization and boundary conditions for numeri-cal experiment: first calculation step

479Bautechnik 83 (2006), Heft 7

ZielfunktionZur Bewertung der durch eine Vorwärtsrechnung erreich-ten Lösung in bezug auf die Referenzdaten („Messungen“der Inklinometer) wurde in Anlehnung an die Methodeder kleinsten Quadrate die in Gl. (7) dargestellte Zielfunk-tion verwendet. Diese vergleicht die während eines Be-rechnungsschritts an q vorgegebenen Punkten auftreten-den Verschiebungen up,x,calc der Vorwärtsrechnung mitden Referenzwerten up,meas für den Parametervektor x.

(7)

Auf die Verwendung einer komplexeren Zielfunktion(z. B. Berücksichtigung von Wichtungsfunktionen, sieheMalecot et al. 2004 [18]) wurde verzichtet.

Im Rahmen des Beispiels 2 wurden folgende Punkt-mengen für die Bestimmung der Abweichung untersucht:– alle Knotenpunkte entlang beider Inklinometer– Knotenpunkte jeweils am Kopf der beiden Inklinome-ter

Inverse Parameterbestimmung: Materialmodell Mohr-CoulombDie Parameter der einzelnen Modellbereiche für das Ma-terialmodell Mohr-Coulomb können Tabelle 2 entnom-men werden. Bild 13 zeigt die aus der Fußentlastung resul-tierenden Verschiebungen der Referenzsimulation („Meß-daten“).

Das E-Modul des Hangs an der Schichtgrenze, dieTeufenlage der Trennfläche und der innere Reibungswin-kel des Hangs wurden als Parameter für eine Identifizie-rung bzw. Optimierung ausgewählt. Die Zahlenwerte die-ser Parameter sollten auf der Basis einer Referenzsimula-tion zurückgewonnen werden (Tabelle 2). Bild 14 zeigteine Visualisierung der Zielfunktionstopologie auf der Ba-sis eines 8 ¥ 8 ¥ 8-Rasters für einen inneren Reibungswin-kel f von 28°. Die Zielfunktion auf der Basis der Kopf-punkte der Inklinometer (Bild 14b) weist einerseits einenweitaus größeren Wertebereich [5,8686E-02 … 3,2861E-08]als auch eine höhere Rauhigkeit gegenüber der Zielfunk-

F xq

u up x calc p measp

q

( ): , , ,= -( )=

Â1 2

1

tion auf der Basis aller Meßwerte (Bild 14a; Wertebereich[3,9141E-03 … 1,9151E-08]) auf. Das gesuchte Minimummit den Koordinaten (ESchichtgrenze = 7,001E+008 N/m2; t = 12,49 m; f = 30,0°) wird mit einem Particle-Swarm-Op-timizer mit 10 Individuen für den Fall, daß die Verschie-bungsvektoren beider Inklinometer einbezogen werden,

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Bild 13. Visualisierung der resultierenden Verschiebungenfür die Entfernung der Fußschüttung der Referenzsimulationmit dem Materialmodell Mohr-CoulombFig. 13. Visualization of final displacement field: referencesimulation applying Mohr-Coulomb constitutive model

Tabelle 2. Materialparameter der Referenzsimulation (Mate-rialmodell: Mohr-Coulomb)Table 2. Material parameter of reference simulation (consi-tutive model: Mohr-Coulomb)

Parameter Einheit Wert Suchbereich

Hang (Mohr-Coulomb)Dichte r [kg/m3] 2200E-Modul an der GOK EGOK [N/m2] 1,00E+08E-Modul an der Schicht- 5,00E+08 … grenze ESchichtgrenze

[N/m2] 7,00E+081,00E+09

Querdehnungszahl n [–] 0,3innerer Reibungswinkel f [°] 30,0 25 … 35Kohäsion c [N/m2] 1,00E+04Teufenlage der Schicht t [m] 12,50 7,5 … 17,5

Fußschüttung (elastisch)Dichte r [kg/m3] 2200E-Modul E [N/m2] 7,00E+09Querdehnungszahl n [–] 0,3

Bild 14. Visualisierung der Ergebnisse der Optimierung mitdem Materialmodell Mohr-Coulomb für j = 28°; a) alleMeßwerte; b) Ansatzpunkte Inklinometer)Fig. 14. Visualization of optimization results with Mohr-Cou-lomb constitutive model for j = 28°; a) all measured data; b) head-points of Inclinometer)

i. d. R. nach 35 bis 40 Berechnungsschritten sehr gut an-genähert. Liegen für eine Optimierung nur die Vektorenan den beiden Ansatzpunkten vor, benötigt der PSO ca. 50bis 55 Berechnungsschritte. Bei der Nutzung der Datendes Ansatzpunktes nur eines Inklinometers konnte keineverläßliche inverse Bestimmung der gesuchten Modellpa-rameter durchgeführt werden, da aufgrund der Rauhigkeitund damit auftretenden Nebenminima keiner der zur Ver-fügung stehenden Optimierungsalgorithmen reproduzier-bar bei Wiederholung der Berechnung unter Verwendungunterschiedlicher initialer Parametersets erfolgreich dasOptimum fand.

Als Basis für die Anwendung der Hyperflächenapproxi-mationsverfahren für diese dreidimensionale Problemstel-lung wurden durch ein Latin-Hypercube-Verfahren 75 Stütz-stellen erzeugt. Zum Vergleich: Ein einfaches 6 ¥ 6 ¥ 6-Ra-ster impliziert bereits 216 Aufrufe der Vorwärtsrechnung.Ein PSO mit 10 Individuen konnte für die Basispunkt-menge mit den vollständigen Daten eines oder beider In-klinometer auf der Approximation das gesuchte Optimummit (ESchichtgrenze = 6,965E+008 N/m2; t = 10,6 m; f = 29,9°)in nur 20 Berechnungsschritten annähern. Die bereits er-wähnte stärkere Rauhigkeit der Zielfunktion auf der Basisausschließlich der Daten der Inklinometerköpfe erzeugtein der Approximation ein zusätzliches Nebenextremumbei (ESchichtgrenze = 6,600E+008 N/m2; t = 8,9 m; f =28,9°), wobei jedoch auch das gesuchte globlale Optimumin ca. 75 % der Durchläufe des PSO auf der Approxima-tion gefunden wurde.

Zusammenfassung Beispiel 2Mit Hilfe der Berechnungen von Beispiel 2 konnte gezeigtwerden, daß die vorgestellte Verfahrensklasse auch bei Pro-blemstellungen mit erhöhten Rauhigkeiten in der Lage ist,Extremwerte, d. h. Modellparameter zu ermitteln. Insbe-sondere bei dem Einsatz von numerischen Verfahren alsVorwärtsrechnung ist mit einem solchen „Rauschen“ derZielfunktionswerte zu rechnen. Weiterhin konnte demon-striert werden, daß dieses Verfahren auch bei Problemstel-lungen mit mehr als zwei Parametern eingesetzt werdenkann.

5 Zusammenfassung und Ausblick

Mit den Hyperflächenapproximationsverfahren steht eineleistungsfähige Algorithmenklasse zu Beschleunigung vonIdentifikations- und Optimierungsaufgaben zur Verfü-gung. Es konnte gezeigt werden, daß die hier vorgestelltenMethoden in der Lage sind, die Anzahl der benötigten

480 Bautechnik 83 (2006), Heft 7

Vorwärtsrechnungen signifikant zu verringern, was auchfür komplizierte Aufgabenstellungen im Bereich desBauingenieurwesens von großem Interesse ist. Grundideedieses Verfahrens ist die Ersetzung einer aufwendig zu be-rechnenden Zielfunktion durch eine deutlich wenigerrechenintensive Näherungsfunktion. Die Suche nach demExtremum erfolgt dann, indem die Näherungsfunktionanstelle der eigentlichen Vorwärtsrechnung aufgerufenwird.

An dieser Stelle sei angemerkt, daß diese Methodekein „Allheilmittel“ für die Beschleunigung und/oder Lö-sung von Optimierungsproblemen darstellt. Auch für die-sen Algorithmus gilt: die Ergebnisse können nicht besserals die verwendeten Daten bzw. das als Basis dienendeSimulationsverfahren sein. Die Methode ist jedoch geeig-net, für bestimmte Aufgabenstellungen durch eine verfei-nerte Auswertung von bereits zur Verfügung stehendenDaten die Lage von Extremwerten zu approximieren.

Hierbei ist diese Verfahrensklasse nicht auf die Opti-mierung bzw. inverse Parameterbestimmung beschränkt.So sind Anwendungen beispielsweise im Bereich der grafi-schen Aufbereitung von Datenstrukturen und Rekonstruk-tion von fehlenden Informationen denkbar.

Zukünftige Arbeiten sollen speziell die Leistungsfähig-keit des Verfahrens für höherdimensionale Problemstel-lungen untersuchen und verbessern.

Literatur

[1] Allen, T. T., Yu, L.: Low-Cost Response Surface Methodsfrom Simulation Optimization. Quality and Reliability Engi-neering International 18 (2002), S. 5–17.

[2] van den Berg, F.: An Analysis of Particle Swarm Optimi-zers. PhD thesis, University of Pretoria (2001).

[3] Biles, W. E.: A response surface method for experimentaloptimization of multiresponse processes. Industrial and Engi-neering Chemistry: Process Design and Development, 14(1975), S. 152–158.

[4] Bolzon, G., Fedele, R., Maier, G. (2002): Parameter identifi-cation of a cohesive crack model by Kalman filter. Computermethody in applied mechanics and engineering 191, ElsevierScience B. V., S. 2847–2871.

[4] Boyd, S., Vandenberghe, L.: Convex Optimization. Cam-bridge University Press 2006.

[5] Bui, L. T., Essam, D.,Abbass, H. A., Green, D.: Performanceanalysis of evolutionary multi-objective optimization methodsin noisy environments. ALAR Technical Report Series, Uni-versity of New South Wales, Australia 2005.

[6] Carrera, J.,Alcolea,A., Medina,A., Hidalgo, J., Slooten, L. J.:Inverse problem in hydrogeology. Hydrogeological Journal13, Springer-Verlag 2005, S. 206–222.

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik

Tabelle 3. Vergleich des Berechnungsaufwands zur Bestimmung des Optimums (Materialmodell: Mohr-Coulomb)Table 3. Comparison of the efficiency of optimization process (constitutive model: Mohr-Coulomb)

Verfahren Berechnungsschritte Funktionsaufrufe

Basisdaten: alle Verschiebungsvektoren beider InklinometerParticle-Swarm-Optimizer 35 … 40 350 … 400Hyperflächenapproximationsmethode + PSO ~ 75

Basisdaten: Verschiebungsvektoren am Kopf beider InklinometerParticle-Swarm-Optimizer 50 … 55 500 … 550Hyperflächenapproximationsmethode + PSO ~ 75

481Bautechnik 83 (2006), Heft 7

[7] Connor, J. J.: Analysis of Structural Member Systems. Mas-sachusetts Institute of Technology, The Ronald Press Com-pany, New York 1976.

[8] Cui, L., Sheng, D.: Genetic algorithms in probabilistic finiteelement analysis of geotechnical problems. Computers andGeotechnics 32 (2006), S. 555–563.

[9] Daugherty, A. F., Turnquist, M. A.: Simulation optimizationusing response surfaces based on spline approximations,Proceedings of the 1980 Winter Simulation Conference,S. 183–193.

[10] De Jong, K.: An analysis of the behaviour of a class of ge-netic adaptive systems. PhD thesis, University of Michigan(1975).

[11] Eberhardt, R. C., Kennedy, J.: A new optimizer using par-ticle swarm theory. Proceedings of the Sixed InternationalSymposium on Micro Machine and Human Science 1995,S. 39–43, Nagoya, Japan, IEEE Service Center, Piscataway,NJ.

[12] Fleischer, J., Broos, A.: Parameteroptimierung bei Werk-zeugmaschinen – Anwendungsmöglichkeiten und Potentiale.Weimarer Optimierungs- und Stochastiktage 1.0, Weimar(2004).

[13] Flores Santiago, O., Bausinger, R.: Automatische Schweiß-punkt-Optimierung an Karosserien. XXV. FEM – Kongreß,Baden-Baden 1998.

[14] Jeong, S.-J.: Ein Beitrag zur Erzeugung nichtlinearer Ent-wurfsseegänge im numerischen Wellenkanal. Dissertations-schrift 2003, Fakultät für Verkehrs- und Maschinensysteme,TU Berlin.

[15] Kennedy, J., Eberhardt, R. C.: Particle Swarm Optimiza-tion. Proceedings of IEEE International Conference on Neu-ral Networks 1995, Volume IV, S. 1942–1948, Perth, Austra-lia, IEEE Service Center, Piscataway, NY.

[16] Kollig, T., Keller, A.: Efficient Multidimensional Sampling.Eurographics, Volume 21 (2002), Nr. 3.

[17] Krüger, J. H.: Echtzeitsimulation und -darstellung von Wol-ken. Diplomarbeit, Fachbereich Informatik, RWTH Aachen(2002).

[18] Malecot, Y., Flavigny, E., Boulon, M.: Inverse Analysis ofSoil Parameters for Finite Elemsent Simulation of Geotechni-cal Structures: Pressuremeter Test und Excavation Problem.R. B. J. Brinkgreve, H. Schad, H. Schweiger, E. Willand (Hrsg.),Proc. Symp. Geotechnical Innovations, Glückauf Verlag, Es-sen 2004, S. 659–675.

[19] Matous, K., Leps, M., Zeman, J., Sejnoha, M.: Applying ge-netic algorithms to selected topics commonly encountered in

engineering practice. Computer methods in applied mecha-nics and engineering 190 (2000), Elsevier, S. 1629–1650.

[20] McKay, M. D., Conover, W. J., Beckman, R. J.: A Compari-son of Three Methods for Selecting Values of Input Variablesin the Analysis of Output from a Computer Code, Technome-trics, 21 (1979), S. 239–245.

[21] Miller, B.: Noise, Sampling, and Efficient Genetic Algo-rithms. PhD thesis, Department of ComputerScience, Univer-sity of Illinoise at Urbana-Champaign 1997.

[22] Polheim, H.: Evolutionäre Algorithmen-Verfahren, Opera-toren und Hinweise für die Praxis. Springer 1999.

[23] Press,W. H., Flannery, B. P., Teukolsky, S. A.,Vetterling,W. T.:Numerical Recipes in C: The Art of Scientific Computing.Cambridge University Press; 2. Auflage 1992.

[24] Schanz, T., Zimmerer, M., Datcheva, M., Meier, J.: Identifi-cation of constitutive parameters for numerical models via in-verse approach. Felsbau, Vol. 24 (2006), No. 2, S. 11–21

[25] Schilling, S.: Beitrag zur Lösung ingenieurtechnischer Ent-wurfsaufgaben unter Verwendung Evolutionärer Algorithmen.Dissertationsschrift 2003, Fakultät Bauingenieurwesen, Bau-haus-Universität Weimar.

[26] Schwarz, S.: Sensitivitätsanalyse und Optimierung bei nicht-linearem Strukturverhalten. Bericht Nr. 34 (2001), Institut fürBaustatik, Universität Stuttgart.

[27] Shi, Y., Eberhardt, R. C.: Parameter Selection in ParticleSwarm Optimization. Evolutionary Programming VII: Proc.EP98, New York: Springer-Verlag, (1998) S. 591–600.

[28] Shi, Y., Eberhardt, R. C.: A modified particle swarm opti-mizer. Proceedings of the IEEE International Conference onEvolutionary Computation, (1998) S. 69–73. Piscataway, NJ:IEEE Press.

[29] Wang, G. G.: Adaptive Response Surface Method UsingInherited Latin Hypercube Design Points. Transactions of theASME, Journal of Mechanical Design, Vol. 125, (2003)S. 210–220.

Autoren dieses Beitrages:Dipl.-Ing. Jörg Meier, Bauhaus-Universität Weimar, Coudraystraße 11c,99421 Weimar, [email protected] Sebastian Rudolph, Universität Karlsruhe (TH), Institut AIFB, 76128 Karls-ruhe, [email protected], [email protected] Prof. Dr.-Ing. habil. Tom Schanz, Bauhaus-Universität Weimar, ProfessurBodenmechanik, Coudraystraße 11c, 99421 Weimar,[email protected]

J. Meier/S. Rudolph/T. Schanz · Effektiver Algorithmus zur Lösung von inversen Aufgabenstellungen – Anwendung in der Geomechanik