Vorwort

Die vorliegende Masterarbeit ist am Institute of Mechanics and Computational Mechanics der Universit¨ at Hannover entstanden.

Ich m¨ ochte mich besonders bei Herrn Prof. Dr.-Ing. Udo Nackenhorst, Frau Dr.-Ing. Karin Wiechmann sowie Herrn Dipl.-Ing. Bastian Ebbecke f¨ ur die intensive Betreuung und f¨ ur die zahlreichen fachlichen Ratschl¨ age bedanken.

Ich hoffe, einen kleinen Einblick in die Zusammenarbeit geben zu k¨ onnen, die zwischen Ingenieuren und Medizinern immer weiter w¨ achst, und mit dieser Arbeit zu weiteren Verbesserun- gen in der Biomechanik beizutragen.

.

Ich versichere, dass ich die vorliegende Arbeit selbst¨ andig und nur unter Benutzung des im Literaturverzeichnis aufgef¨ uhrten Schrifttums erstellt habe.

Hannover, 19 September 2005

Inhaltsverzeichnis

1  Problemstellung und Zusammenfassung  9

2  Knochenumbau  10

2.1  Zusammenhang zwischen Elastizit atsmodul und Dichte  11

2.2  Evolutionsgleichung des isotropen Knochenumbaus  12

3  Finite Elemente in der linearen Elastizit atstheorie  13

3.1  Kontinuumsmechanische Grundgleichungen  14

3.2  Prinzip der virtuellen Verr uckung  14

3.3  FE-Diskretisierung  16

4  Sensitivit atsanalyse  17

4.1  Numerische Sensitivit atsanalyse  18

4.2  Analytische Sensitivit atsanalyse  20

5  Optimierungsverfahren  25

5.1  Optimalit atsbedingungen  25

5.1.1  Notwendige Bedingung erster Ordnung  26

5.1.2  Notwendige Bedingung zweiter Ordnung  26

5.1.3  Hinreichende Bedingungen zweiter Ordnung  26

5.2  Ein allgemeines Abstiegsverfahren  27

5.3  Wolfe-Powell-Schrittweitenstrategie  27

5.4  Globalisiertes Quasi-Newton-Verfahren  28

5.5  Trust-Region-Techniken  32

5.6  Das Gauss-Newton-Verfahren  33

4

5.7 Das Levenberg-Marquardt-Verfahren . . . . . . . . . . . . . . . . . . . . . . . 34 5.8 Ableitungsfreie Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.8.1 Genetische Algorithmen und Evolutionsstrategien . . . . . . . . . . . . 35

6 Kopplung des Genetischen-Algorithmus und den Gradienten-Verfahren 37

6.1 Numerische Testbeispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.2 Ein einfaches 1D-Beispiel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 6.3 2D -und 3D-Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 6.3.1 Beispiel 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.3.2 Beispiel 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6.3.3 Beispiel 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 6.3.4 Beispiel 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.4 Diskussion der numerischen Ergebnisse aus 2D -und 3D-Berechnungen . . . . 53

7 Anwendung am Femur 54

7.1 Berechnungen mit zwei Kr¨ aften . . . . . . . . . . . . . . . . . . . . . . . . . 56 7.2 Berechnungen mit drei Kr¨ aften . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.3 Berechnungen mit f¨ unf Kr¨ aften . . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.4 Berechnungen mit f¨ unf Kr¨ aften und den Genetischen-Algorithmus als Vorrechnung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

8 Diskussion der Ergebnisse 69

9 Zusammenfassung und Ausblick 71

5

Abbildungsverzeichnis

6.1 Funktion mit lokalen Minima . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 6.2 Effizienzerh¨ ohung durch Kopplung der Verfahren . . . . . . . . . . . . . . . . 39 6.3 Fachwerk mit unbekannter Belastung . . . . . . . . . . . . . . . . . . . . . . 40 6.4 System mit unbekannter Belastung . . . . . . . . . . . . . . . . . . . . . . . . 42 6.5 Fehler in Abh¨ angigkeit der verschiedenen Lastkombinationen . . . . . . . . . 42 6.6 Fehler in Abh¨ angigkeit der verschiedenen Lastkombinationen . . . . . . . . . 43 6.7 links: Zugkraft, rechts: Druckkraft . . . . . . . . . . . . . . . . . . . . . . . . 44 6.8 Ein wichtiger Bereich der Fehlerkurve . . . . . . . . . . . . . . . . . . . . . . 46 6.9 System mit unbekannter Belastung . . . . . . . . . . . . . . . . . . . . . . . . 47 6.10 Fehler in Abh¨ angigkeit der verschiedenen Lastkombinationen . . . . . . . . . 47 6.11 links: Zubgbelastung, rechts: Druckbelastung . . . . . . . . . . . . . . . . . . 48 6.12 System mit unbekannter Belastung . . . . . . . . . . . . . . . . . . . . . . . . 50 6.13 links: Zugbelastung, rechts: Druckbelastung . . . . . . . . . . . . . . . . . . . 51 6.14 3D-System mit unbekannter Belastung . . . . . . . . . . . . . . . . . . . . . . 52

7.1 Dichteverteilung aus CT-Daten . . . . . . . . . . . . . . . . . . . . . . . . . . 54 7.2 Femur mit f¨ unf Kr¨ aften bzw. 15 Kraftkomponenten . . . . . . . . . . . . . . . 56 7.3 Dichteverteilung: Quasi-Newton-Verfahren, 6 Designvariablen . . . . . . . . . 57 7.4 Dichteverteilung: Gauss-Newton-Verfahren, 6 Designvariablen . . . . . . . . . 58 7.5 Dichteverteilung: Levenberg-Marquardt-Verfahren, 6 Designvariablen . . . . . 58 7.6 Dichteverteilung: Quasi-Newton-Verfahren, 9 Designvariablen . . . . . . . . . 60 7.7 Dichteverteilung: Gauss-Newton-Verfahren, 9 Designvariablen . . . . . . . . . 61 7.8 Dichteverteilung: Levenberg-Marquardt-Verfahren, 9 Designvariablen . . . . . 62 7.9 Dichteverteilung: Quasi-Newton-Verfahren, 15 Designvariablen . . . . . . . . 64 7.10 Dichteverteilung: Gauss-Newton-Verfahren, 15 Designvariablen . . . . . . . . 64

6

7.11 Dichteverteilung: Levenberg-Marquardt-Verfahren, 15 Designvariablen . . . . 65 7.12 Dichteverteilung: Gauss-Newton-Verfahren, 15 Designvariablen und GA als Vorrechnung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 7.13 Dichteverteilung: Levenberg-Marquardt-Verfahren, 15 Designvariablen und GA als Vorrechnung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

8.1 links: optimale Dichteverteilung, die restlichen drei: berechnet . . . . . . . . . 70

7

Tabellenverzeichnis   

6.1  Berechnungsergebnisse  44

6.2  Berechnungsergebnisse  48

6.3  Berechnungsergebnisse  51

6.4  Berechnungsergebnisse  53

7.1  Gemessene Kr afte beim normalen Gehen  55

7.2  Berechnungsergebnisse f ur zwei Kr afte  56

7.3  Berechnete Lasten, 6 Designvariablen  57

7.4  Berechnungsergebnisse f ur drei Kr afte  59

7.5  Berechnete Lasten, 9 Designvariablen  60

7.6  Berechnungsergebnisse f ur f unf Kr afte  62

7.7  Berechnete Lasten, 15 Designvariablen  63

7.8  Berechnete Lasten, 15 Designvariablen  66

7.9  Berechnete Lasten, 15 Designvariablen  66

8

Kapitel 1

Problemstellung und Zusammenfassung

In der vorliegenden Arbeit wird ein Beitrag zur numerischen Behandlung der inversen Kr¨ afteberechnung f¨ ur Aufgabenstellungen aus der Biomechanik vorgestellt. Die Behandlung inverser Probleme ist insbesondere aus mathematischer Sicht bez¨ uglich der Existenz und Regularit¨ at von L¨ osungen nicht trivial. Auf die theoretischen Aspekte soll hier nicht eingegangen werden. Vielmehr soll das entwickelte algorithmische L¨ osungskonzept zur Behandlung inverser Probleme im Vordergrund stehen.

Diese Arbeit beschr¨ ankt sich auf die Untersuchungen am menschlichen Femur. Dabei wird das komplexe dynamische Problem der Lasten am Femur auf statisch ¨ aquivalente Lasten zur¨ uckgef¨ uhrt. Die Idee besteht darin, aus einer bekannten bzw. mittels CT ¹ -Aufnahmen erfassten optimalen Dichteverteilung die Lasten zu ermitteln. Dadurch k¨ onnen z.B. individuelle Therapien f¨ ur Patienten entwickelt werden. Es wird festgestellt, dass diese Problemstellung um ein vielfaches komplexer als die Problemstellung des Knochenumbaus ist. Beim Knochenumbau wird aus einer vektoriellen Gr¨ osse, n¨ amlich der vorgegebenen Kraft, eine skalare Gr¨ osse also die Dichte ermittelt. Das inverse Problem stellt den umgekehrten Fall dar. Aus dieser Problemstellung resultiert eine nichtlineare Optimierungsaufgabe, die mit Hilfe der Algorithmen der mathematischen Optimierung behandelt wird. Zur L¨ osung der Optimierungsaufgabe werden sowohl ableitungsfreie als auch Gradienten-Verfahren eingesetzt, wobei die gradientenbasierten Verfahren im Vordergrund stehen. Um mit Gradienten-Verfahren die resultierende Optimierungsaufgabe l¨ osen zu k¨ onnen, wird eine Sensitivit¨ atsanalyse implementiert. Im Kapitel 6 wird gezeigt, dass f¨ ur das inverse Problem mehrere gleichwertige L¨ osungen vorhanden sein k¨ onnen.

¹ Abk¨ urzung f¨ ur Computer Tomographie

9

Kapitel 2

Knochenumbau

In diesem Kapitel werden die wichtigsten Aspekte des Knochenumbauprozesses und der Modellbildung zusammengefasst, die f¨ ur das weitere Verst¨ andnis der Arbeit sehr hilfreich und zum Teil notwendig sind. Bekanntlich h¨ angt die Festigkeit eines Bauteils nicht allein von den Materialeigenschaften ab. Zu einem erheblichen Teil wird sie durch die strukturellen Elemente bestimmt. Schon vor mehr als einem Jahrhundert wurde in Wolff[29] das Ph¨ anomen des Knochenumbaus untersucht und beschrieben, das auch als Wolffsches Transformationsgesetz bekannt ist. Die wesentliche Aussage des Wolffschen Gesetzes lautet: Jeder Ver¨ anderung in der Funktion eines Knochens folgt eine definierte Ver¨ anderung in der inneren Architektur und ausseren Gestalt in ¨ Ubereinstimmung mit mathematischen Gesetzten. Die Untersuchung des der ¨

Knochenumbauprozesses ist nicht Gegenstand dieser Arbeit. Aus diesem Grund wird hier nicht auf Einzelheiten eingegangen. Interessierte Leser werden besonders auf die Ver¨ offentlichungen von Nackenhorst [21], Nackenhorst & Schr¨ oder[19], Wolff[29] und Weng[24] verwiesen. Da die tats¨ achlichen Vorg¨ ange des Knochenumbaus sehr komlex sind und die mathematischen Gesetzm¨ assigkeiten bisher noch nicht vollst¨ andig formuliert wurden, wird bei praktischen Problemstellungen von der folgenden Aussage ausgegangen: Knochen sind in ihrer ¨ ausseren Gestalt und innerer Architektur optimal an die mechanischen Belastungen angepasst, das heisst, sie formen mit einem Minimum an Materialeinsatz ein Maximum an mechanischer Festigkeit. Knochen haben - entsprechend ihren einzelnen Funktionen - unterschiedliche Formen aber ein gemeinsames architektonisches Grundprinzip, n¨ amlich das Prinzip ” form follows function“.

Der menschliche Femur, der auch in dieser Arbeit untersucht wurde ist ein sogenannter R¨ ohrenknochen. Er hat einen hohlen Schaft, in welchem sich das Knochenmark befindet.

10

2.1 Zusammenhang zwischen Elastizit¨ atsmodul und Dichte

Der Zusammenhang zwischen den elastischen Eigenschaften und der Dichte des Knochens wird durch die folgende empierische Formel beschrieben

wobei E 0 ein mittlerer Elastizit¨ atsmodul bzw. ρ 0 eine mittlere Dichte zu Beginn des Umbauprozesses sind. Sowohl das Elastizit¨ atsmodul E als auch die Dichte ρ k¨ onnen nur Werte annehmen, die positiv und kleiner oder gleich des kortikalen Knochens sind

0 ≤ ρ ≤ ρ c 0 ≤ E ≤ E c bzw. (2.2)

F¨ ur die numerische Behandlung empfiehlt sich die folgende Normierung der Dichte

aus (2.1) und (2.3) folgt

E = E 0 λ ⁿ (2.4)

F¨ ur die Spannungs-Dehnungs-Beziehung in der linearen Elastizit¨ atstheorie gilt der folgende konstitutive Zusammenhang σ = Cε (2.5)

und mit (2.4) folgt

σ = λ ⁿ C 0 ε (2.6)

f¨ ur die Verzerrungsenergiedichte gilt folgendes

wobei U die Verzerrungsenergie und C der vierstufige Materialtensor ist. (2.6) in (2.7) eingesetzt liefert

Dass der Exponent n = 2 sein muss, kann durch eine konsistente thermodynamische Betrachtung begr¨ undet werden, Nackenhorst & Kristin & Lammering [20].

11

2.2 Evolutionsgleichung des isotropen Knochenumbaus

Basis der Knochenumbausimulation ist eine Evolutionsgleichung, die eine Beziehung zwischen mechanischer Beanspruchung und der daraus resultierenden Zu- bzw. Abnahme der Knochensubstanz ausdr¨ uckt, Nackenhorst[21].

Andere Theorien f¨ uhren auf etwas kompliziertere Evolutionsgleichungen. In dieser Arbeit wird die obere Beziehung verwendet, die auch als Hamburger Theorie bekannt ist, Ebbecke & Nacken-horst [22, 23]. Gleichung (2.9) besagt, dass die zeitliche ¨ Anderung der Dichte linear von der

Verzerrungsenergiedichte abh¨ angt. Da weder die tats¨ achliche Steigung der Geraden, also die Konstante k, noch die Referenzverzerrungsenergiedichte ψ ref bekannt sind, ist man gezwungen weitere Annahmen zu treffen. In der Literatur wurden Berechnungen f¨ ur unterschiedlieche k und ψ ref angegeben. F¨ ur die weiteren Berechnungen werden die ermittelten ψ ref direkt ¨ ubernommen und k = ρ 0 gesetzt. Dadurch resultiert die folgende Evolutionsgleichung

Ist ψ gr¨ osser als der Referenzwert ψ ref , so kommt es zum Anbau von Knochenmasse. Ist allerdings ψ < ψ ref , so wird Knochenmasse abgebaut. Die Evolutionsgleichung (2.10) wird mit Hilfe der Finiten-Differenzen-Methode gel¨ ost:

Das Eulervorw¨ artsverfahren eingesetzt in (2.11) liefert die folgende Beziehung

F¨ ur die Aufdatierung der normierten Dichte folgt

λ ^t+∆t = λ ^t + ∆λ ^t+∆t (2.13)

Die Eigenschaften des Verfahrens werden im Kapitel “Sensitivit¨ atsanalyse“ diskutiert.

12

Kapitel 3

Finite Elemente in der linearen

Elastizit¨ atstheorie

Die Betrachtung von statischen Gleichgewichtszust¨ anden von K¨ orpern im Rahmen der linearen Elastizit¨ atstheorie f¨ uhrt auf partielle Differentialgleichungen. Bei der numerischen Behandlung von partiellen Differetialgleichungen wird heute sehr h¨ aufig die Methode der Finiten Elemente eingesetzt. Wegen ihrer Flexibilit¨ at erlaubt die Methode auch die Behandlung schwieriger Probleme, weil sie - anders als Differenzenverfahren oder Berechnungen mit Finiten Volumen - auf die Variationsformulierung der Differentialgleichung zugeschnitten ist. Im Gegensatz zu gew¨ ohnlichen Differentialgleichungen gibt es bei partiellen Differentialgleichungen nicht immer klassische L¨ osungen und oft nur sogenannte schwache L¨ osungen. Im n¨ achsten Abschnitt werden die kontinuumsmechanischen Grundgleichungen vorgestellt, aus denen anschliessend die schwache Formulierung hergeleitet wird. Die Grundlagen dieses Verfahrens sind in vielf¨ altiger Form in der Literatur dargestellt, siehe z.B. (Bathe[8], Braess[9]). Aus diesem Grund werden in dieser Arbeit nur die Grundz¨ uge des Vorgehens aufgezeigt um insbesondere ein Hinter-grundwissen f¨ ur die weiteren Betrachtungen zu bekommen. Da im Rahmen dieser Arbeit ein instation¨ ares Problem (Evolutionsgleichung der Dichte) behandelt wird, ist es empfehlenswert die Zeitableitungen durch Differenzenverfahren, die Geometrie sowie das Verschiebungsfeld durch finite Elemente zu approximieren. Diese Vorgehensweise ist auch als Semidiskretisierung bekannt.

13

3.1 Kontinuumsmechanische Grundgleichungen

Die statische Feldgleichung ergibt sich durch Aufstellen der Gleichgewichtsbeziehungen am infinitesimal kleinen Volumenelement. Des Weiteren m¨ ussen noch die kinematischen Feldgleichungen sowie die konstitutive Annahme ¨ uber das Werkstoffgesetz ber¨ ucksichtigt werden.

Statische Feldgleichung

Divσ + ρb = 0 (3.1)

Konstitutive Beziehung

σ = C : ε (3.2)

Kinematische Feldgleichung

Zu diesem Satz von Gleichungen kommen noch die Randbedingungen hinzu. Diese sind die Spannungsrandbedingungen auf dem Spannungsrand ∂Ω t und die Verschiebungsrandbedingungen auf dem Verschiebungsrand ∂Ω u

wobei n der Normalenvektor auf dem Rand ist

3.2 Prinzip der virtuellen Verr ¨ uckung

Ausgangspunkt f¨ ur eine FE-Formulierung bildet das Prinzip der virtuellen Verr¨ uckung. Zun¨ achst werden die einzelnen Schritte vorgestellt, die zur Herleitung der schwachen Form des Gleichgewichts erforderlich sind. Durch Multiplikation der Ausgangsgleichung (3.1) mit einer zul¨ assigen Testfunktion η, die die wesentlichen Randbedingungen erf¨ ullt, und Integration ¨ uber das

gesamte Gebiet Ω, erh¨ alt man schliesslich eine ¨ aquivalente Formulierung von (3.1)

Analog geht man auch bei (3.5) vor, und erh¨ alt

(3.6) und (3.7) ergeben

Das Ausmultiplizieren und die Anwendung von Rechenregeln ¹ zur Divergenzbildung eines Feldes auf den ersten Integranden liefert

Die Anwendung des Divergenztheorems ² auf den ersten Term f¨ uhrt zu

∂Ωt

Umformung des ersten und des letzten Terms liefert

∂Ωt

Nach (3.5) gilt auf dem Rand ∂Ω t der Zusammenhang σn = t. Daraus folgt, dass der erste und der letzte Term in (3.11) sich gegenseitig aufheben. Es resultiert die folgende Beziehung, die auch als schwache Form des Gleichgewichts bekannt ist.

Setzt man f¨ ur die Spannung σ die Beziehung (3.2) ein, so erh¨ alt man die sehr oft benutzte Darstellung der schwachen Form des Gleichgewichts

¹ hier wurde die folgende Rechenregel benutzt:

Div (T v) = T : Gradv + DivT v =⇒ DivT v = Div (T v) − T : Gradv

wobei T ein Tensor- und v ein Vektorfeld ist

² Divergenztheorem:

Da in der schwachen Form des Gleichgewichts (3.13) nur erste Ableitungen auftreten, muss f¨ ur das Verschiebungsfeld u die C ⁰ -Stetigkeitsanforderung erf¨ ullt werden. Das gilt nat¨ urlich auch f¨ ur die Testfunktion η.

3.3 FE-Diskretisierung

In der Regel ist es nicht m¨ oglich, Testfunktionen zu finden, die die wesentlichen Randbedingungen im gesamten Gebiet erf¨ ullen. Die Grundidee der Finite Elemente Methode besteht darin, das Gebiet Ω in Teilgebiete Ω e (Elemente) zu unterteilen und Ansatzfunktionen in den Elementen einzuf¨ uhren. Anschliessend findet eine Assemblierung statt, die dann auf ein lineares Gleichungssystem f¨ uhrt. Der Unbekanntenvektor beinhaltet gerade die gesuchten sogenannten Knotenwerte. Die Systemmatrix ist zun¨ achst singul¨ ar. Erst durch Einsetzen der Randbedingungen erh¨ alt man eine positiv definite Matrix. Dadurch wird das reduzierte Gleichungssystem l¨ osbar. Im n¨ achsten Kapitel wird die variationelle Sensitivit¨ atsanalyse etwas ausf¨ uhrlicher be-handelt, die aber analog zur allgemeinen Vorgehensweise der FE-Methode hergeleitet werden kann. Aus diesem Grund wird hier nicht weiter auf Details eingegangen.

16

Kapitel 4

Sensitivit¨ atsanalyse

Im Rahmen der vorliegenden Arbeit nimmt die Sensitivit¨ atsanalyse eine zentrale Rolle ein. Aus diesem Grund wird in diesem Kapitel auf die prinzipiellen M¨ oglichkeiten n¨ aher eingegangen. Insbesondere die Vorgehensweise zur Ermittlung der Sensitivit¨ aten mit Hilfe der numerischen bzw. der analytischen Sensitivit¨ atsanalyse wird ausf¨ uhrlich beschrieben und diskutiert. Ausgehend von einem Parametersatz (Designvariablen) wird bei der Sensitivit¨ atsanalyse ein Teil oder alle Parameter nacheinander variiert, um feststellen zu k¨ onnen, wie die Zielgr¨ osse auf die Parameter¨ anderungen reagiert. Dadurch findet man die Parameter heraus, die die Zielgr¨ osse stark beeinflussen, und diejenigen, auf deren ¨ Anderung die Zielgr¨ osse kaum reagiert. Dies bedeu-

tet, mit Hilfe der Sensitivit¨ atsanalyse kann die Parameterempfindlichkeit einer L¨ osung ermittelt

werden. Aus der Zielgr¨ osse λ die mittels FEM berechnet wird und den gegebenen Daten ˆ λ, die

mittels CT-Aufnahmen bestimmt wurden, wird ein Mass f¨ ur die Abweichung zwischen Ziel-gr¨ osse und den vorgegebenen Daten ˆ λ berechnet. Dieses Mass f (λ) h¨ angt nat¨ urlich implizit

auch von den unbekannten Designvariablen ab und soll minimiert werden. f (λ) wird als Zielfunktion bezeichnet, die wie folgt definiert wird

1

λ − ˆ λ ² f = (4.1)

2

Auf Minimierung der Zielfunktion wird im Kapitel ” Optimierungsverfahren“ n¨ aher eingegangen. Man k¨ onnte sich auch eine andere Funktion definieren z.B. die Summe der absoluten Abweichungen. Die Summe der Abweichungsquadrate (4.1) als Zielfunktion kann wahrschein-lichkeitstheoretisch begr¨ undet werden. Die Methode der kleinsten Quadrate (least square fit) hat aber auch praktische Vorteile. Der wesentliche Vorteil ist die Differenzierbarkeit der Zielfunktion, was eine Erleichterung f¨ ur die numerische Suche des Minimums ist.

17

Da die Implementierung der Sensitivit¨ atsanalyse mit einem zus¨ atzlichen Speicher-, Rechen-und Programmieraufwand verbunden ist, sollen genau diese Kritererien f¨ ur die Effizienz der in den n¨ achsten Abschnitten hergeleiteten Alogrithmen massgebend sein. In der Regel besitzen Dichteoptimierungsprobleme eine geringe Anzahl von Designvariablen, so dass auch grosse, aktuelle Forschungsgegenst¨ ande wie der menschliche Femur mit einem vertretbarem Aufwand behandelt werden k¨ onnen.

Man wird feststellen, dass die Sensitivit¨ atsanalyse eine entscheidende Komponente bei der numerischen L¨ osung von Optimierungsaufgaben mittels Gradienten-Verfahren darstellt. Ihre Qualit¨ at beeinflusst sehr stark die Rechenzeit. Aus diesem Grund ist eine m¨ oglichst effiziente Variante anzustreben. F¨ ur eine kompakte Darstellung der Sensitivit¨ atsanalyse wird auf die Arbeiten von Schwarz[12], Meyer[10] , Barthold[11] und Wiechmann[15] verwiesen, die auch nichtlineare Probleme behandeln.

4.1 Numerische Sensitivit¨ atsanalyse

Die bekanntesten und einfachsten Verfahren zur Implementierung der numerischen Sensitivit¨ atsanalyse sind das Euler-Vorw¨ arts- und das zentrale Differenzen-Verfahren. Wegen der einfachen Implementierung und dem geringen Formulierungsaufwand sind diese Verfahren in der Praxis sehr beliebt und finden eine breite Anwendung, auch ausserhalb der numerischen Sensitivit¨ atsanalyse . Die erw¨ ahnten Verfahren k¨ onnen durch eine abgebrochene Taylorreihenentwicklung hergeleitet werden. Im Rahmen dieser Arbeit wurde die numerische Sensitivit¨ atsanalyse mit Hilfe des Euler-Vorw¨ arts-Verfahrens realisiert. Das Verfahren besitzt zwar nur eine Genauigkeit erster Ordnung, hat aber den Vorteil, dass weniger Rechenoperationen als beim zentralen Differenzen-Verfahren durchgef¨ uhrt werden m¨ ussen. Dadurch wird der Rechenaufwand, der beim inversen Problem sehr gross ist, in Grenzen gehalten. Um den Rechenaufwand der numerischen Sensitivit¨ atsanalyse besser absch¨ atzen zu k¨ onnen, wird zun¨ achst der untenstehende Algorithmus betrachtet, der die wesentlichen Bausteine, die im Rahmen dieser Arbeit ben¨ otigt werden, enth¨ alt.

18

Man stellt fest, dass f¨ ur die Umsetzung der numerischen Sensitivit¨ atsanalyse lediglich die Funktionswerte zu bestimmen sind. Genau aus diesem Grund nimmt der numerische Aufwand bei steigender Anzahl der Designvariablen sehr stark zu, da f¨ ur jede Designvariable eine zus¨ atzliche Anzahl FE-Berechnungen durchgef¨ uhrt werden muss (Zeitabh¨ angiges Problem: Evolutionsgleichung f¨ ur die Dichte), um die “gest¨ orten“ Funktionswerte zu berechnen. Bei zentralen Differenzen w¨ are die Anzahl der Funktionsauswertungen doppelt so gross. Ein wesentlicher Nachteil ist aber auch die Abh¨ angigkeit der Ergebnisse von der gew¨ ahlten St¨ orung ∆. M¨ ochte man z.B. auf der sicheren Seite liegen und die St¨ orung sehr klein w¨ ahlen, so k¨ onnten weitere unerw¨ unschte Eigenschaften wie Rundungs- und Konditionsfehler auftretten. W¨ ahlt man jedoch die St¨ orung ∆ zu gross, so treten Ungenauigkeiten aufgrund der abgebrochenen Taylorreihenentwicklung auf. Auf die Einzelheiten wird nicht weiter eingegangen und insbesondere auf die Arbeit von Kimmich [13] verwiesen. Obwohl die erw¨ ahnten Schwierigkeiten in der vorgestellen Variante auftreten k¨ onnen, liefert die Methode in vielen F¨ allen hinreichend gute Ergebnisse.

19

Diese Methode kann z.B. zur Verifikation der komplexen analytischen Verfahren dienen, oder wenn wenig Designvariablen vorhanden sind. Bei einer grossen Anzahl von Designvariablen ist die Methode nicht zu empfehlen, Darauf wird in den Testbeispielen noch mal eingegangen.

4.2 Analytische Sensitivit¨ atsanalyse

F¨ ur die Formulierung der analytischen Sensitivit¨ atsanalyse sind vertiefte kontinuumsmechanische und mathematische Kenntnisse notwendig. Dazu kommen noch die verschiedenen Vorgehensweisen wie die semianalytische, diskrete , variationelle usw., die auf unterschiedliche Varianten der analytischen Sensitivit¨ asanalyse f¨ uhren. Im Rahmen dieser Arbeit wurde die variationelle Variante formuliert und implementiert. Bei der variationellen Sensitivit¨ asanalyse wird zuerst die Berechnung der Variationen des kontinuierlichen Problems nach den Designvariablen durchgef¨ uhrt. Anschliessend erfolgt die Diskretisierung der kontinuierlichen Sensitivit¨ atsaussagen und somit eine ¨ Uberf¨ uhrung der Variationen in Ableitungen. Zun¨ achst wird die Variation der schwachen Form des Gleichgewichts im i-ten Zeitschritt betrachtet

Der interne Anteil der schwachen Form des Gleichgewichts (3.13) ist

Variation von (4.3) bez¨ uglich Designvariablen f¨ uhrt zu

Diskretisierung von (4.4) liefert die Element-Sensitivit¨ atsmatrix

Die Berechnung der ˜ D-Matrix wird exemplarisch f¨ ur den ebenen Spannungszustand gezeigt

20

Es wird festgestellt, dass nur der Vorfaktor ge¨ andert wird und die Struktur erhalten bleibt. ¨ Ahnliches gilt auch f¨ ur die Element-Sensitivit¨ atsmatrix S e verglichen mit der Element-Steifigkeitsmatrix. Beide Matrizen haben die gleiche Struktur, was eine Erleichterung f¨ ur die Implementierung ist. Jetzt geht es darum, δ s λ zu bestimmen. Dazu geht man wie folgt vor:

Im i-ten Zeitschritt wird die Variation der Verschiebung δ s u i bestimmt, aus der dann die Variation der Verzerrungen δ s ε i berechnet werden kann. Daraus lassen sich die Variation der Verzerrungsenergiedichte δ s ψ i und die Variation des Dichteinkrements δ s (∆λ) und schliesslich die gesuchte δ s λ i berechnen. (4.2) umgeformt liefert

δ s u i = −K ⁻¹ (S δ s λ i−1 − δ s F) (4.6)

Daraus wird die Variation der Verzerrung berechnet

δ s ε i = B δ s u i (4.7)

Die Variation der Verzerrungsenergiedichte l¨ asst sich wie folgt ausrechnen

Aus (4.8) kann

berechnet werden. Schliesslich folgt noch die Aufdatierung der Variation der normierten Dichte

¨ Ahnlich wird auch bei der Variation des externen Anteils δ s G ^ext vorgegangen, wobei hier keine Volumenkr¨ afte ber¨ ucksichtigt werden. Da es sich um diskrete Einzellasten handelt, kann das Integral durch eine Summe ersetzt werden. F¨ ur die Variation des externen Anteils gilt

(4.11) liefert einen Vektor, der nur an der Stelle i eine 1 sonst ¨ uberall Nullen hat. In der n¨ achsten

Abbildung ist der Algorithmus dargestellt, der die hergeleiteten Aussagen zusammenfasst.

22

Nachdem der obige Algorithmus durchlaufen wurde, kann die Zielfunktion (4.1) und deren Ableitung

f = δ s λ(λ − ˆ λ) (4.12)

berechnet werden.

Mit Hilfe des obigen Algorithmus soll nun der wesentliche Rechenaufwand und der Speicherbedarf abgesch¨ atzt werden. Innerhalb eines Zeitschrittes sind die folgenden Rechenoperationen durchzuf¨ uhren:

• Assemblierung der globalen Steifigkeitsmatrix K i und L¨ osen des reduzierten linearen Gleichungssystems. K i hat die Dimension n × n, ist symmetrisch und d¨ unnbesetzt → als Sparse-Matrix speichern

• Berechnung der Verzerrungen durch eine Matrix-Vektor-Multiplikation

• Ermittlung der Verzerrungsenergiedichte f¨ ur jedes Element durch kleine Matrix-Vektor-Multiplikationen

• Assemblierung der globalen Sensitivit¨ atsmatrix. Die Sensitivit¨ atsmatrix ist nahezu vollbesetzt, unsymmetrisch und hat die Dimension n × 1. F¨ ur die Implementierung wurde sie in eine n × n Matrix umgeformt und als Sparse-Matrix abgespeichert. Dadurch wurde die Implementierung einfacher

• F¨ ur jede Designvariable ist δ s u i = −K ⁻¹ i (S i δ s λ i−1 − δ s F i ) zu l¨ osen, wobei K ⁻¹ i aus der Strukturanalyse ¨ ubernommen wird.

• Berechnung der Variation der Verzerrungen durch kleine Matrix-Vektor-Multiplikation

Zusammenfassend kann gesagt werden, dass f¨ ur die Senistivit¨ atsanalyse in gleicher Gr¨ ossen-ordnung Speicherplatz wie bei der Strukturanalyse bereitgestellt werden muss. Auch die Rechenzeit ist etwa gleich gross, da Matrix-Vektor-Produkte durchgef¨ uhrt werden, die in etwa auch bei der Strukturanalyse in ¨ ahnlicher Weise vorkommen.

24

Kapitel 5

Optimierungsverfahren

F¨ ur die diskutierte Problemstellung ist eine unrestringierte Optimierungtsaufgabe, genauer gesagt eine unrestringierte Minimierungsaufgabe zu l¨ osen. In diesem Kapitel wird ein kurzer ¨ Uberblick der implementierten und den von MATLAB bereitgestellten Verfahren zur L¨ osung der Minimierungsaufgabe gegeben. Zun¨ achst werden einige grundlegende Merkmale von Minimierungsproblemen vorgestellt.

Gegeben sind eine Menge X ⊆ R ⁿ und eine Funktion f : X → R.

Gesucht ist ein x ^∗ ∈ X mit der Eigenschaft

f (x ^∗ ) ≤ f (x) ∀ x ∈ X

Im weiteren wird f¨ ur das unrestringierte Minimierungsproblem die folgende Notation benutzt

min f (x) (P U )

x∈X

wobei f (x) die Zielfunktion (4.1) ist, die minimiert werden soll. Der Einfachheit halber wird x in diesem Kapitel das Argument sein. Zun¨ achst werden f¨ ur das unrestringierte Problem (PU) notwendige und hinreichende Bedingungen zur Charakterisierung von lokalen L¨ osungen vorgestellt.

5.1 Optimalit¨ atsbedingungen

F¨ ur die folgenden Betrachtungen wird vorausgesetzt, dass X ⊆ R ⁿ nicht leer und offen ist.

25

5.1.1 Notwendige Bedingung erster Ordnung

Die Funktion f sei in x ^∗ ∈ X stetig differenzierbar. Ist x ^∗ lokales Minimum von (PU), dann gilt

f (x ^∗ ) = 0 (5.1)

Optimierungsverfahren finden in der Regel eine station¨ are Stelle von f . Einfache 1-D Beispiele wie z. B. f (x) = x ³ oder f (x) = −x ² zeigen, dass (5.1) nur eine notwendige, aber keine hinreichende Optimalit¨ atsbedingung ist. Daraus folgt, dass ein station¨ arer Punkt von f keineswegs ein lokales oder globales Minimum sein muss.

5.1.2 Notwendige Bedingung zweiter Ordnung

Die Bedingung (5.1) gilt auch in einem lokalen Maximum von f . Um lokale Minima und Maxima zu unterscheiden, wird jetzt eine notwendige Optimalit¨ atsbedingung zweiter Ordnung betrachtet.

Die Funktion f sei in einer Umgebung vo x ^∗ ∈ X zweimal stetig differenzierbar. Ist x ^∗ lokales Minimum von (PU), dann gilt (5.1) und

x ^{T 2} f (x ^∗ ) x ≥ 0 ∀ x ∈ R ⁿ (5.2)

d.h. ² f ist positiv semidefinit.

5.1.3 Hinreichende Bedingungen zweiter Ordnung

Auch die Bedingungen (5.1) und (5.2) zusammen sind nicht hinreichend daf¨ ur, dass x ^∗ ein lokales Minimum ist. Das kann man wieder mit einfachen 1-D Beispielen wie f (x) = x ³ mit x ^∗ = 0 zeigen. Hier gilt zwar (5.1) und (5.2), es liegt aber kein Minimum vor. Um eine hinreichende Optimalit¨ atsbedingung f¨ ur ein lokales Minimum zu erhalten, muss man die Bedingung (5.2) etwas versch¨ arfen.

Seien X ⊆ R ⁿ und f : X → R zweimal stetig differenzierbar. Gelten

• f (x ^∗ ) = 0 und

• ² f (x ^∗ ) positiv definit

26

so ist x ^∗ ein striktes Minimum von f .

5.2 Ein allgemeines Abstiegsverfahren

Die zentrale Idee des Allgemeinen Abstiegsverfahrens zur L¨ osung von (PU) l¨ asst sich sehr leicht beschreiben: Ist man an einem Punkt x ∈ R ⁿ angelangt, so sucht man sich eine Richtung d ∈ R ⁿ aus, in der es bergab geht. Entlang dieser Richtung geht man dann solange, bis man den Funktionswert von f hinreichend verkleinert hat. Definition Abstiegsrichtung

Seien f : R ⁿ → R und x ∈ R ⁿ . Ein Vektor d ∈ R ⁿ heisst Abstiegsrichtung von f in x, wenn es ein t > 0 gibt mit ∀ t ∈ f (x + td) < f (x) 0, t

Das folgende Lemma liefert ein hinreichendes Kriterium f¨ ur das Vorliegen einer Abstiegsrichtung

Lemma 5.2 Seien f : R ⁿ → R stetig differenzierbar, x ∈ R ⁿ und d ∈ R ⁿ mit f (x) ^T d < 0. Dann ist d eine Abstiegsrichtung von f in x.

5.3 Wolfe-Powell-Schrittweitenstrategie

In der Literatur werden viele M¨ oglichkeiten angegen um geeignete Schrittweiten zu bestimmen, siehe z.B. Geiger& Kanzow [3], Alt[1] und Polak[2]. In der erw¨ ahnten Literatur wird gezeigt, dass f¨ ur das Quasi-Newton-Verfahren, das sp¨ ater behandelt wird, die Wolfe-Powell-Regel zur Bestimmung der Schrittweite effizient eingesetzt werden kann. Die weiteren Schrittweiten-Regeln, wie z.B. die Armijo-Regel, die insbesondere beim Newton-Verfahren und seinen Varianten benutzt wird oder die strenge Wolfe-Powell-Regel, die sich beim CG-Verfahren eingeb¨ urgert hat, werden nicht weiter betrachtet. Der interessierte Leser sei auf die erw¨ ahnte Literatur verwiesen.

0, 1 Sei f : R ⁿ → R stetig differenzierbar, und seien die Zahlen σ ∈ , ρ ∈ [σ, 1) fest

2

vorgegeben. Zu x, d ∈ R ⁿ mit f (x) ^T d < 0 bestimme man ein t > 0 mit

und

f (x + td) ^T d ≥ ρf (x) ^T d (W P 2)

Zur Veranschaulichung wird ϕ := f (x + td) gesetzt.Damit lauten die Wolfe-Powell-Bedingungen (WP1),(WP2):

ϕ (t) ≤ ϕ (0) + σtϕ (0)

und

ϕ (t) ≥ ρϕ (0)

Die Schrittweite ist somit aus einem Bereich zu w¨ ahlen, in welchem einerseits der Graph von ϕ unterhalb der Geraden ¹ durch (0, ϕ (0)) mit der skalierten Steigung ϕ (0) verl¨ auft, und in dem andererseits, grob gesprochen, der Graph von ϕ nicht mehr so steil wie in der N¨ ahe von t = 0 abf¨ allt oder bereits ansteigt. (WP1) garantiert eine gewisse Reduktion von f . Falls (WP1) durch die Wahl von t nicht gesichert ist, halbiert man t solange, bis die Bedingung erf¨ ullt ist. Um un¨ otig kleine Schrittweiten zu vermeiden muss zus¨ atzlich (WP2) erf¨ ullt werden.

5.4 Globalisiertes Quasi-Newton-Verfahren

Zun¨ achst werden die wesentlilchen Eingenschaften des Newton-Verfahrens vorgestellt, da das Newton-Verfahren wegen seiner guten Konvergenzeigenschaften Grundlage f¨ ur viele Abstiegsverfahren ist. Diese Vor¨ uberlegungen sollen auch einen besseren ¨ Ubergang zum Quasi-Newton-Verfahren erm¨ oglichen.

Die Funktion f : R ⁿ → R sei zweimal stetig differenzierbar. Die zentrale Idee des Newton-Verfahrens besteht darin, das unrestringierte Minimierungsproblem

min f (x)

zu l¨ osen, indem man sukzessiv die quadratischen N¨ aherungen

q k (x) := f

zu minimieren versucht, wobei x ^k ∈ R ⁿ den aktuellen Iterationspunkt bezeichnet. Ist die Hesse- Matrix ² f positiv definit, so ist x ^k+1 genau dann L¨ osung von

min q k (x) , (5.4)

¹ In der Literatur wird h¨ aufig auch von Armigo Goldstein Geraden gesprochen

28

wenn x ^k+1 der Bedingung

f¨ ur einen station¨ aren Punkt von q k gen¨ ugt. Hieraus folgt f¨ ur x k+1

Um die explizite Berechnung der inversen Hesse-Matrix zu vermeiden, bestimmt man zun¨ achst eine Suchrichtung d ^k ∈ R ⁿ als L¨ osung des linearen Gleichungssystems

² f d ^k = −−f x k (5.7)

und setzt anschliessend x ^k+1 := x ^k + d ^k . Auf diese Weise erh¨ alt man den gleichen Vektor wie in (5.6). In der Literatur wird (5.7) oft als Newton-Gleichung bezeichnet. Dieses Verfahren ist als lokales Newton-Verfahren bekannt. Das Verfahren hat zwar sehr gute, d.h. superlineare bzw. quadratische Konvergenz in der N¨ ahe der L¨ osung, hat aber den wesentlichen Nachteil, dass die Newton-Gleichung (5.7) nicht immer eine L¨ osung im aktuellen Iterationspunkt besitzt. Auch dann, wo (5.7) eine L¨ osung bestitzt, muss diese nicht unbedingt eine Abstiegsrichtung sein. Um diese unerw¨ unschten Eigenschaften zu vermeiden, wird das Verfahren modifiziert. Es resultiert ein sogenanntes globalisiertes ² Newton-Verfahren. Die Idee dieses globalisierten Verfahrens besteht darin, die L¨ osung d ^k der Newton-Gleichung (5.7) als Suchrichtung zu w¨ ahlen und entlang dieser Richtung eine Schrittweitenbestimmung mittels der weiter oben erw¨ ahnten Armigo-Regel durchzuf¨ uhren. Sollte die L¨ osung der Newton-Gleichung nicht existieren oder ungeeignet f¨ ur die Suchrichtung sein, dann f¨ uhrt man i. a. einen Gradientenschritt, d.h. man x k w¨ ahlt d ^k = −−f , um eine sichere Abstiegsrichtung zu ermitteln. Da die Begriffe wie superlineare bzw. quadratische Konvergenz schon benutzt wurden, wird noch eine allgemeine Definition der Konvergenzgeschwindigkeit der Folgen angegeben:

Sei {x ^k } ⊂ R ⁿ eine Folge mit Grenzwert x ^∗ . Die Folge heisst linear konvergent, wenn es ein 0 < L < 1 gibt mit

x ^k+1 − x ^∗ ≤ Lx ^k − x ^∗ ur k ≥ 0 f ¨

² In der Literatur wird h¨ aufig auch von ged¨ ampftes Newton-Verfahren gesprochen, d.h. die Begriffe globalisiert

und ged¨ ampft werden als Synonyme benutzt. Auch in dieser Arbeit werden beide Varianten benutzt

29

Die Folge {x ^k } heisst superlinear konvergent, wenn folgendes gilt

Die Folge {x ^k } heisst quadratisch konvergent, wenn es ein c > 0 gibt mit

x ^k+1 − x ^∗ ≤ cx ^k − x ∗ ² ur k ≥ 0 f ¨

Mit Hilfe der betrachteten Vor¨ uberlegungen wird ein globalisiertes Quasi-Newton-Verfahren vorgestellt. ¨ Anlich wie bei den Schrittweitenstrategien gibt es auch hier eine Klasse der sogenannten Quasi-Newton-Verfahren. Da die Berechnung der zweiten partiellen Ableitungen der

Hesse-Matrix ² f

Matrix eine hinreichend gute Approximation durch m¨ oglichst einfach berechenbare Matrizen H k zu erreichen. Das ist auch der wesentliche Unterschied zum Newton-Verfahren. Diese Approximation der Hesse-Matrix wird in jeder Iteration aufdatiert, so dass in jedem Schritt die linearen Gleichungssysteme sehr viel leichter zu l¨ osen sind (auch durch Aufdatieren) als beim Newton-Verfahren. Es gibt unterschiedliche Aufdatierungsformeln, wobei die bekanntesten die PSB ³ -, DFP ⁴ und die BFGS ⁵ -Formel sind. In der Praxis hat sich die BFGS-Updateformel am besten bew¨ ahrt. Es wird kurz die Grundidee des Quasi-Newton-Verfahrens skiziert. Durch eine Taylor-Reihenentwicklung des Gradienten der Zielfunktion um den Punkt x ^k+1 erh¨ alt man die folgende Gleichung

x k+1 + ² f x ^k+1 − x k x ^k+1 − x ^k f = f + o (5.8)

Durch elementare Umformungen und gewisse Anforderungen ⁶ an die Approximation H k+1 der Hessematrix erh¨ alt man die untenstehende Gleichung

x ^k+1 − x k = f − f H k+1 (5.9)

(5.9) wird oft auch als Quasi-Newton-Gleichung oder auch als Sekantengleichung bezeichnet. Der folgende Algorithmus soll eine Zusammenfassung des globalisierten Quasi-Newton-Verfahrens mit der direkten BFGS-Aufdatierungsformel und der Wolfe-Powell-Schrittweitenstrategie

³ Powell-symmetric-Broyden

⁴ Davidon-Fletcher-Powell

⁵ Broyden-Fletcher-Goldfarb-Shanno

⁶ siehe dazu Geiger & Kanzow S. 130

30

geben.

Die in A5 angegebene Formel f¨ ur H k+1 ist die schon erw¨ ahnte BFGS-Updateformel der Approximation der Hesse-Matrix. Auch hier gilt zu erw¨ ahnen, dass die globale Konvergenz durch Gradientenschritte erzwungen werden kann, falls die Winkelbedingung ⁷ nicht erf¨ ullt ist. Allerdings stellt sich dann die Frage, wie die neue Matrix H k+1 berechnet werden soll. Auf die verschiedenen M¨ oglichkeiten zur Berechnung der neuen Matrix wird weiter nicht eingegangen und auf den Abschnitt 11.5 von Kosmol[6] verwiesen. Es wird festgestellt, dass das globalisierte Quasi-Newton-Verfahren einen kleineren numerischen Aufwand als das globalisierte Newton-Verfahren hat, aber den Nachteil, dass das Verfahren lokal nur eine superlinerare Konvergenz erreichen kann. Dies allerdings unter der Voraussetzung, dass die durch den Algorithmus er- gleichm¨assig von 90 ^◦ weg be- x k ⁷ Die Winkelbedingung besagt, dass der Winkel zwischen d ^k und −−f

schr¨ ankt bleibt

31

zeugte Folge H k gegen ² f (x ^∗ ) konvergiert. Dazu muss noch sichergestellt werden, dass stets die Schrittweite t k = 1 gew¨ ahlt wird, sofern diese die Wolfe-Powell-Bedingungen gen¨ ugt. Unter gewissen Voraussetzungen kann dies gezeigt werden, siehe dazu Geiger & Kanzow[3] S. 168. Es wird festgestellt, dass viele Voraussetzungen erf¨ ullen werden m¨ ussen. Dennoch beobachtet man in der Praxis h¨ aufig lokal eine schnelle Konvergenz. In den folgenden Abschnitten dieses Kapitels werden noch einige Gradienten-Verfahren vorgestellt, die, wie sich im laufe der Zeit dieser Arbeit herausgestellt hat, geeigneter f¨ ur unsere Aufgabenstellung sind. Allerdings beschr¨ ankt man sich nur auf das wesentliche, da eine ausf¨ uhrliche Erl¨ auterung den Rahmen dieser Arbeit bei weitem sprengen w¨ urde.

5.5 Trust-Region-Techniken

Auch hier gibt es eine Klasse der Trust-Region-Verfahren. Dabei handelt es sich um Kombinationen von Verfahren, wie z.B. das Trust-Region-Newton-Verfahren oder auch das Trust-Region-Quasi-Newton-Verfahren usw. Die Grundidee ist hierbei, dass die quadratische Approximation

nur f¨ ur kleine d hinreichend gut ist. Deswegen schr¨ ankt man die Verwendung von (5.10) in einem Bereich (Trust-Region) ein, in dem man der quadratischen Approximation trauen darf. Die Suchrichtung berechnet man z.B. aus der restringierten Minimierungsaufgabe

Dabei muss noch beachtet werden, dass der Trust-Region-Radius ∆ k geeignet gew¨ ahlt wird. Mit (5.11) wird deutlich, dass in jedem Iterationsschritt ein quadratisches Minimierungsproblem (Teilproblem) mit Nebenbedingungen zu l¨ osen ist. Weiter wird festgestellt, dass beim Trust-Region-Verfahren keine Schrittweitenstrategie ben¨ otigt wird. Als Mass f¨ ur die G¨ ute der

32

Approximation dient der folgende Quotionet

Falls ρ k < 0 ist d ^k zu verwerfen und der Schritt zu wiederholen. Falls ρ k nahe an 1 liegt, stimmt die tats¨ achliche Reduktion (Z¨ ahler) gut mit der vorhergesagten Reduktion (Nenner) ¨ uberein und

x ^k+1 = x ^k + d ^k kann als neuer Punkt akzeptiert werden. Dabei kann der Radius ∆ k in vielen F¨ allen sogar vergr¨ ossert werden.

5.6 Das Gauss-Newton-Verfahren

In den vorigen Abschnitten wurden einige Verfahren vorgestellt, die in vielen F¨ allen f¨ ur allgemeine Optimierungsprobleme ohne Nebenbedingungen eingesetzt werden. Nun werden neue Verfahren vorgestellt, die besonders f¨ ur spezielle nichtlineare Optimierungsaufgaben, n¨ amlich f¨ ur nichtlineare Ausgleichsprobleme geeignet sind. Die zu minimierende Zielfunktion (4.1) hat die Struktur eines nichtlinearen Ausgleichsproblems ohne Nebenbedingungen.

Es sei F : R ⁿ → R ^m eine zweimal stetig differenzierbare Abbildung und m ≥ n. Die allgemeine Form des Ausgleichsproblems ist

In der Regel ist m n, also ein ¨ uberbestimmtes Gleichungssystem. Es ist naheliegend das

Newton-Verfahren zur Minimierung von f anzuwenden. Dazu ben¨ otigt man den Gradient f (x) und die Hesse-Matrix ² f (x). F¨ ur die Zielfunktion des Ausgleichsproblems (5.13) gilt

f (x) = F (x) ^T F (x) (5.14)

und

² f (x) = F (x) ^T F (x) + R (x)

Oft geht man davon aus, dass f¨ ur eine L¨ osung x ^∗ von (5.13) das Residuum klein ist und daher auch die Werte von F i (x) klein sind. Durch diese Annahme wird die Vernachl¨ assigung von R (x) gerechtfertigt, die schliesslich auf die folgende Approximation f¨ uhrt:

Wendet man das ged¨ ampfte Newton-Verfahren auf das Problem (5.13) an und ersetzt ² f (x) durch die obige Approximation F (x) ^T F (x), d.h. ersetzt man das Problem (5.4) bzw. (5.5) beim Newton-Verfahren zur Richtungsbestimmung durch das Problem

so erh¨ alt man das Gauss-Newton-Verfahren. Mit der folgenden Berechnung wird gezeigt, dass die Gauss-Newton-Richtung d k GN ∈ R ⁿ aus der L¨ osung des resultierenden linearen Ausgleichproblems bestimmt werden kann.

Ist x ^k eine N¨ aherungsl¨ osung des nichtlinearen Ausgleichsproblems (5.13), so gilt f¨ ur x ^k+1 = x ^k + d ^k :

x k+1 f

Es wird festgestellt, dass beim Gauss-Newton-Verfahren nur die ersten Ableitungen auftauchen und dass ausgehend von einer Startn¨ aherung x ⁰ eine Folge linearer Ausgleichsprobleme zu l¨ osen ist. Schiesslich wird noch die Gauss-Newton-Richtung d k GN angegeben, die die Anforde-

rungen von Lemma 5.2 erf¨ ullt und somit eine Abstiegsrichtung ist

Auch das Gauss-Newton-Verfahren konvergiert in der Regel nur lokal, d.h. bei gen¨ ugend guter Startn¨ aherung. Konvergenz h¨ angt in der Regel sehr stark von der Gr¨ osse des Residuums F (x ^∗ ) , sowie von der St¨ arke der Nichtlinearit¨ at ab. Lokal superlineare oder quadratische Konvergenz erh¨ alt man nur unter zus¨ atzlichen Voraussetzungen, z.B. wenn F (x ^∗ ) = 0 ist, siehe dazu Alt[1] S.158 und die dort angegebene Literatur. Da wie schon erw¨ aht die Gauss-Newton-Richtung eine Abstiegsrichtung ist, kann analog zu den vorangegangenen Abschnitten eine D¨ ampfungsstrategie eingef¨ uhrt werden.

5.7 Das Levenberg-Marquardt-Verfahren

Die Kombination von Trust-Region-Techniken mit dem Gauss-Newton-Verfahren wird als Levenberg-Marquardt-Verfahren bezeichnet. Es muss lediglich die in Abschnitt 2.5 verwendete quadratische Approximation modifiziert werden.

34

m k (d) =

Zu l¨ osen ist in jedem Schritt also das folgende restringierte Problem

Man sieht, dass das gleiche Problem wie beim Gauss-Newton-Verfahren zu l¨ osen ist, allerdings unter der obigen Nebenbedingung. Minimierungsprobleme mit Nebenbedingungen sind i.a. schwieriger und aufwendiger zu behandeln als solche ohne Nebenbedingungen. MATLAB stellt diverse Implementierungen von Verfahren f¨ ur restringierte Probleme zur Verf¨ ugung.

5.8 Ableitungsfreie Verfahren

In vielen F¨ allen ist die Berechnung der Ableitungen der Zielfunktion sehr aufwendig oder gar nicht vorhanden. Deshalb wurden Verfahren entwickelt, die keine Ableitungen der Zielfunktion ben¨ otigen. Solche Verfahren sind sehr zeitaufwendig, haben aber den Vorteil, dass sie i.d.R. unabh¨ angig vom Startpunkt ein globales Minimum finden, was bei den Gradienten-Verfahren i.a. nicht der Fall ist. Es gibt unterschiedliche ableitungsfreie Verfahren. Im n¨ achsten Abschnitt wird die “biologisch” motivierte Evolutionsstrategie kurz vorgestellt. Auf Grund der Komlexit¨ at und unterschiedlichen Varianten des Verfahrens wird nicht auf Einzelheiten eingegangen, siehe z.B. Schwefel[7]. Auch hierf¨ ur stellt MATLAB Impementierungen von Genetischen Al-gorithmen zur Verf¨ ugung, die auch im Rahmen dieser Arbeit benutzt wurden. Es wird weiterhin das Problem (PU) betrachtet.

5.8.1 Genetische Algorithmen und Evolutionsstrategien

Die Grundidee ist, ¨ ahnlich der biologischen Evolution, eine Anzahl Individuen (L¨ osungskandidaten) zuf¨ allig zu erzeugen und diejenigen auszuw¨ ahlen, die einem bestimmten G¨ utekriterium am besten entsprechen. Deren Eigenschaften werden dann leicht ver¨ andert und miteinander

35

kombiniert, um eine neue Generation zu erzeugen. F¨ ur die Minimierungsaufgabe (4.1) wird wohl der Funktionswert als G¨ utekriterium am besten geeignet sein. Der typische Genetischer Algorithmus umfasst die folgenden Schritte:

• Initialisierung: Erzeugen einer ausreichend grossen Menge unterschiedlicher “Individuen” (L¨ osungskandidaten).

• Evaluation: F¨ ur jeden einzelnen L¨ osungskandidat wird anhand der Zielfunktion ein Wert bestimmt.

• Selection: Zuf¨ allige Auswahl von L¨ osungskandidaten aus der vorhandenen L¨ osungskandidatenmenge. Dabei werden L¨ osungskandidaten mit besseren Zielfunktionswerten mit einer h¨ oheren Wahrscheinlichket ausgew¨ ahlt.

• Rekombination: die Gene verschiedener Idividuen werden gemischt und aus den neuen Parametern eine neue Generation von Individuen erzeugt (Vermehrung)

• Mutation: Zuf¨ allige Ver¨ anderung der Wertekombinationen der Individuen der neuen Generation.

Der Alorithmus wird ab Schritt zwei (Evaluation) wiederholt, oder nach einem Abbruchkriterium beendet und der beste verf¨ ugbare L¨ osungskandidat als L¨ osung definiert. Als Abbruchkriterium k¨ onnte man z.B. die maximale Anzahl der Generationen angegeben. Im Vergleich zu den Gradienten-Verfahren k¨ onnen hier keine aussagekr¨ aftige Konvergenzeigenschaften hergeleitet bzw. bewiesen werden, da hierf¨ ur sehr wenig Konvergenztheorie vorhanden ist.

36

Kapitel 6

Kopplung des Genetischen-Algorithmus

und den Gradienten-Verfahren

Im vorigen Kapitel wurden die wesentlichen Eigenschaften der gradientenbasierten und der gradientenfreien Verfahren diskutiert. Dabei ist man zum Schluss gekommen, dass Gradienten-Verfahren eine hohe Konvergenzgeschwindigkeit und eine relativ kleine Konvergenzsicherheit haben. Es wurde gezeigt, dass Gradienten-Verfahren nur lokale Minima finden k¨ onnen. Ein Spezialfall bei dem das lokale Minimum gleichzeitig das globale Minimum ist, sind konvexe Funktionen, die aber in dieser Arbeit nicht vorkommen. Weiter wurden ableitungsfreie Verfahren diskutiert, die eine sehr kleine Konvergenzgeschwindigkeit, daf¨ ur aber eine grosse Konvergenzsicherheit haben. Es liegt auf der Hand, die beiden Verfahren so zu kombinieren, dass ein neuer Algorithmus entsteht, der sowohl eine grosse Konvergenzsicherheit als auch eine Konvergenzgeschwindigkeit gew¨ ahrleistet. Das kann z.B. erreicht werden, indem zuerst mit dem Genetischen Algorithmus m¨ oglichst viele lokale Minima beseitigt werden, um so eine gute Ausgangslage f¨ ur die Gradienten-Verfahren zu schaffen, d.h. die L¨ osungswerte des Genetischen-Algorithmus werden als Startwerte f¨ ur die Gradienten-Verfahren benutzt. Abbildung 6.1 soll die diskutierten Eigenschaften der Methoden und deren Kopplung, mit dem Ziel, die Effizienz der Algorithmen zu erh¨ ohen, veranschaulichen. Es ist ein interessanter Ausschnitt einer 1D-Funktion dargestellt, die i.a. nicht explizit bekannt ist. In dieser Arbeit ist die zu minimierende Funktion das Residuum. Aus diesem Grund wird zwischen diesen beiden Begriffen nicht unterschieden.

37

• P1: Man geht davon aus, dass geeignete Startwerte vorhanden sind und versucht mit gradientenbasierten Verfahren die Funktion zu minimieren. Das Verfahren bricht ab, wenn ein Extremum gefunden wird. In diesem Fall handelt es sich um ein lokales Minimum in P2.

• P2: Ist P2 erreicht, so f¨ uhren gradientenbasierte Methoden nicht weiter. Es werden Suchbereiche in der Umgebung des lokalen Minimums definiert und man versucht, mit dem Genetischen-Algorithmus geeignete Startwerte f¨ ur die gradientebasierten Methoden zu finden. Dazu muss der Suchbereich hinreichend gross definiert werden, um die gew¨ unschte Reduktion der Funktion zu gew¨ ahrleisten, denn nur eine Verkleinerung der Funktion ist ein Indiz daf¨ ur, dass das lokale Minimum ¨ ubersprungen wurde (siehe Abb. 6.1).

• P3: Sobald ein kleinerer Funktionswert als in P2 durch den Genetischen - Algorithmus erzielt wurde, werden wieder gradientenbasierte Verfahren zur Minimierung der Funktion eingesetzt, bis wieder die Situation wie in P2 auftritt, d.h. bis ein lokales Minimum gefunden wird.

38

• P4: In P4 ist wieder die gleiche Vorgehensweise wie in P2 zu empfehlen. In vielen praktischen Problemen ist man zufrieden, wenn das Residuum um einen gewissen Betrag verkleinert wurde, denn das Residuum auf Null zu minimieren, wird in vielen F¨ allen nicht m¨ oglich sein. Grund daf¨ ur sind z.B. die Diskretisierung und insbesondere die vielen Modellannahmen. F¨ ur analytische Funktionen k¨ onnte in P4 das globale Minimum vorhanden sein.

Das folgende Flussdiagram soll die Kopplung des Genetischen-Algorithmus mit den Gradienten-Verfahren schematisch darstellen

Die Vorg¨ ange k¨ onnen problemabh¨ angig mehrmals wiederholt werden. In vielen F¨ allen ist es besser, wenn man zuerst mit den Gradienten-Verfahren beginnt, denn es kann sein, dass die zeitaufwendigen Berechnungen des Genetischen-Algorithmus gar nicht ben¨ otigt werden.

39

6.1 Numerische Testbeispiele

Um das inverse Problem verst¨ andlicher zu machen, wird zun¨ acht ein einfaches 1D Beispiel, das von Hand gel¨ ost werden kann, vorgestellt. Anschliessend werden 2D -und 3D Testbeispiele vorgestellt, die nicht mehr ohne den Rechner gel¨ ost werden k¨ onnen. Mit diesen einfachen Testbeispielen sollen einerseits die entwickelten Algorithmen auf Effizienz und Robustheit gepr¨ uft werden und andererseits sollen sie zum Verst¨ andnis dieser Arbeit beitragen.

6.2 Ein einfaches 1D-Beispiel

Zun¨ achst wird das folgende System betrachtet

F¨ ur dieses Beispiel sind sowohl die optimale Dichte der jeweiligen St¨ abe, die kinematischen Randbedingungen als auch der Lastangriffspunkt bekannt. (2.8) umgeformt und vereinfacht liefert

Wenn in (2.10) ˙ λ = 0 wird, also keine zeitliche ¨ Anderung der Dichte statt findet, dann ist

ψ = ψ ref . Diese Bedingung muss f¨ ur jeden Stab gefordert werden, was mit (6.1) auf folgende

40

Gleichung f¨ ur ψ ref f¨ uhrt

Gleichung (6.2) wird nach N aufgel¨ ost weden

und daraus werden die Normalkr¨ afte der St¨ abe bestimmt. Mit der Formulierung der Gleichgewichtsbedingungen am Lastangriffspunkt k¨ onnen die einzelnen Lastkomponenten der Kraft F berechnet werden.

F V = N 1 · 0.6 (6.4)

und

F H = −N 1 · 0.8 − N 2 (6.5)

6.3 2D -und 3D-Beispiele

Im vorigen Abschnitt konnte festgestellt werden, dass das inverse Problem f¨ ur eindimensionale Problemstellungen keine grossen Schwierigkeiten bereitet. F¨ ur zwei -und dreidimensionale Systeme wird das nicht mehr der Fall sein. Es wird die Verzerrungsenergiedichte (2.8) betrachtet, die f¨ ur den ebenen Spannungszustand wie folgt angegeben werden kann

Aus (6.6) wird die Schwierigkeit der Problematik klar, da aus der skalaren Gr¨ osse ψ die vektorielle Gr¨ osse u berechnet werden muss, aus der dann die Verzerrung ε bestimmt werden kann. Da f¨ ur die Testbeispiele sowohl die kinematischen Randbedingungen als auch die Lastangriffspunkte bekannt sind, kann mit einer freigew¨ ahlten Belastung

Ku = F (6.7)

gerechnet werden. Beim Femur m¨ ussen allerdings sowohl die kinematischen Randbedingungen als auch die Lastangriffspunkte durch diverse Annahmen festgelegt werden.

41

6.3.1 Beispiel 1

Gegeben ist das folgende System

Sowohl die kinematischen Randbedingungen als auch der Lastangriffspunkt sind bekannt. Die optimale Dichteverteilung kann ebenfalls der oberen Abbildung entnommen werden. Da in diesem Beispiel nur eine Last angreift und somit zwei Kraftkomponenten vorhanden sein m¨ ussen, kann der Fehler durch unterschiedliche Lastkombinationen berechnet und grafisch dargestellt werden.

In Abb. 6.5 sieht es so aus, als ob unendlich viele Extremalstellen vorhanden w¨ aren. Das liegt aber daran, dass durch die relativ grobe Disktetisierung 10 × 10 die hohen “Spitzen“ generiert werden. Es wurde eine zus¨ atzliche Berechnung mit doppelt sovielen Elementen durchgef¨ uhrt, die eine sehr glatte Fehlerkurve ergibt.

W¨ urde man allerdings noch feiner diskretisieren, so w¨ aren keine Spitzen mehr zu sehen. Daf¨ ur w¨ urde aber die Rechenzeit um ein Vielfaches steigen. F¨ ur die Diskretisierung von 10 × 10 Elementen wurden 13 h und f¨ ur die 20 × 20 Elemente 50 h ben¨ otigt. Da in beiden F¨ allen die gleichwertigen globalen Minima sehr deutlich zu sehen sind, werden die n¨ achsten Beispiele nur mit der Diskretisierung von 10 × 10 Elementen durchgef¨ uhrt, um so den Rechenaufwand in Grenzen zu halten. Die beiden globalen Minima deuten darauf hin, dass es zwei gleichwertige L¨ osungen geben muss. Die L¨ osung (−70000, 70000) ^T entspricht einer Zugkraft und dementsprechend (70000, −70000) ^T einer Druckkraft.

43

Diese Bilder best¨ atigen die Aussage aus Kapitel 2, dass Knochen optimale Strukturen sind, die mit minimalem Materialeinsatz maximale Festigkeit erreichen. In Abb. 6.7 ist der Dichteverlauf sehr sch¨ on zu sehen, der einer Fachwerkkonstruktion sehr ¨ ahnlich ist. Mit den im Kapitel

“Optimierungsverfahren“ vorgestellten Methoden wurden diverse Berechnungen durchgef¨ uhrt. Damit die ¨ Ubersicht weiterhin erhalten bleibt, betrachtet man nur die Ergebnisse einer Berechnung bei der als Startwert (−55000, 100000) ^T gew¨ ahlt wurde, etwas ausf¨ uhrlicher.

analytisch

Anzahl Iterationen 11 5 6

Funktionsauswertungen 57 30 70

Zeit [s]

Residuum

L¨ osung F

numerisch

Anzahl Iterationen 11 5 6

Funktionsauswertungen 57 30 70

Der oberen Tabelle kann man entnehmen, dass das Gauss-Newton-Verfahren eine viel gr¨ ossere

44

Konvergenzgeschwingikeit als die beiden anderen Verfahren hat. Das kann einerseits aus dem kleinen (bzw. Null) Residuum resultieren und andererseits kann die schwache Nichtlinearit¨ at zwischen Start- und L¨ osungspunkt zu dieser Konvergenzrate gef¨ uhrt haben (siehe dazu Abschnitt 5.6) . Obwohl auch die Levenberg-Marquardt-Methode f¨ ur nichtlineare Ausgleichsprobleme zugeschitten ist, sind die Unterschiede zu der Quasi-Newton-Methode nicht wie erwartet. Ganz im Gegensatz zu den Erwartungen, ist in diesem Fall das Quasi-Newton-Verfahren besser als das Levenberg-Marquardt-Verfahren. Das kann daran liegen, dass bei kleinen Problemen, wie in diesem Beispiel, beim Levenberg-Marquardt-Verfahren die in jedem Schritt zu l¨ osenden Minimierungsprobleme mit Nebenbedingungen zu mehr Rechenaufwand f¨ uhren. Auch der Unterschied zwischen der numerischen und der analytischen Sensitivit¨ atsanalyse ist schon bei dieser Berechnung mit zwei Designvariablen gut zu sehen. Ein deutlicher Unterschied wird erst bei mehreren Designvariablen zu beobachten sein, dazu die kommenden Beispiele. F¨ ur die Konvergenz gegen die exakte L¨ osung spielt die Diskretisierung eine wesentliche Rolle. Die Berechnungen mit einer 10 × 10 Diskretisierung haben nicht zum Ziel gef¨ uhrt, da f¨ ur jeden beliebigen Startpunkt, in dessen unmittelbaren Umgebung ein lokales Minimium gefunden wurde. F¨ ur ung¨ unstige Startwerte f¨ uhren gradientenbasierte Verfahren nicht zum Ziel. Diese Problematik soll der folgende Ausschnitt Abb. 6.6 veranschaulichen.

45

Als Startpunkt wurde P 1 = (−80000, 15000) ^T gew¨ ahlt. Mit dem Gauss-Newton-Verfahren wurde nach 931s P 2 = (−52192, 21895) ^T erreicht. Dabei wurde das Residuum von 825.3 in P 1 auf 303.4 in P 2 verkleinert. In Abb. 6.8 sieht man sehr sch¨ on, dass mit P 2 nur ein lokales Minimumum gefunden wurde. Um das globale Minimum erreichen zu k¨ onnen, wird, wie am Anfang dieses Kapitels beschrieben, die Kopplung der Verfahren ben¨ otigt. F¨ ur den Genetischen-Algorithmus wird ein Suchbereich definiert der z.B. um ±8% vom P 2 abweichen soll. Weiter werden f¨ ur die Startrechnung 20 “Individuen“ generiert die sich ¨ uber 2 Generationen vermehren

sollen. Mit diesen Parametern sollen mit dem Genetischen-Algorithmus geeignete Startwerte f¨ ur das Gauss-Newton-Verfahren gefunden werden.

Die Berechnung mit dem Genetischen-Algorithmus dauerte 256s und lieferte als L¨ osung P 3 = (−78692, 58951) ^T mit einem Residuum von 135.6. P 3 wurde wieder als Startwert f¨ ur das Gauss-Newton-Verfahren benutzt, welches dann das globale Minmum in P 4 = (−70000, 70000) ^T nach 400s finden konnte.

46

6.3.2 Beispiel 2

Das Vorgehen ist analog zum ersten Beispiel

Auch in diesem Fall sind sowohl die kinematischen Randbedingungen als auch der Lastangriffspunkt bekannt. Da der Dichteverlauf nicht symmetrisch ist, sollten auch die Kraftkomponenten unterschiedlich gross sein. Zun¨ achst wird wieder die Fehlerkurve betrachtet, die durch unterschiedliche Kraftkombinationen ermittelt wurde.

Auch in dieser Grafik sind zwei globale Minima zu sehen. Da der Kurvenverlauf in der Um-

47

gebung der Minima sehr flach und langgezogen ist, wird die Minimierung der Funktion etwas schwieriger als im ersten Beispiel sein.

Die Kraftkomponenten sind wie erwartet unterschiedlich. Auch hier sollen die mit den Optimierungsverfahren erzielten Ergebnisse diskutiert werden. Dazu wird eine Berechnung betrachtet, bei der als Startpunkt (70000, −20000) ^T gew¨ ahlt wurde.

analytisch

Anzahl Iterationen 11 6 7

Funktionsauswertungen 69 37 74

Zeit [s] 702 369 734

Residuum

L¨ osung F

numerisch

Anzahl Iterationen 11 6 7

Funktionsauswertungen 69 37 74

Die Gauss-Newton-Methode ist weiterhin die effizienteste Methode. Der Unterschied zwischen der Quasi-Newton- und der Levenberg-Marquardt-Methode ist zwar kleiner geworden, aber

48

immer noch nicht zufrieden stellend. Bei gr¨ osseren Beispielen sollte ein solcher Vergleich die ¨ Uberlegenheit der Levenberg-Marquardt-Methode deutlich machen. Es gibt aber Startwerte wie z.B. (100000, 100000) ^T bei dem die Levenberg-Marquardt-Methode nicht das globale Minimium finden konnte. das Gleiche gilt auch f¨ ur die Quasi-Newton- und die Gauss-Newton-Methode, die mit dem Startpunkt (80000, −40000) ^T nicht gegen die exakte L¨ osung konvergierten. Aufgrund der unterschiedlichen Suchrichtungen der Verfahren, kann es dazu kommen, dass durch die relativ grobe Diskretisierung an unterschiedlichen Stellen lokale Minima generiert wurden, die zum Abbruch der Berechnungen gef¨ uhrt haben. Da f¨ ur diverse Startpunkte etwa die gleichen Unterschiede zwischen den Verfahren zu beobachten waren, werden die diskutierten Ergebnisse als repr¨ asentativ f¨ ur dieses Beispiel betrachtet.

49

6.3.3 Beispiel 3

Im Gegensatz zu den ersten zwei Beispielen wird hier das System an zwei verschiedenen Stellen belastet

F¨ ur dieses Beispiel wurde keine Fehlerkurve durch unterschiedliche Lastkombinationen ermittelt, da einerseits die Berechnungen f¨ ur mehr als zwei Kraftkomponenten sehr lange dauern und andererseits kann keine aussagekr¨ aftige Darstellung erzeugt werden. In der folgenden Tabelle sind die Ergebnisse der Berechnungen f¨ ur den Startpunkt (−10, −320, −280, −330) ^T dargestellt. Alle drei Verfahren haben den exakten L¨ osungsvektor (−100, −300, −300, −300) ^T geliefert.

50

analytisch

Anzahl Iterationen 34 6 13

Funktionsauswertungen 200 49 122

Zeit [s] 989 243 604

10 ⁻¹¹ 10 ⁻²⁰ 10 ⁻¹⁵ Residuum

numerisch

Anzahl Iterationen 34 6 13

Funktionsauswertungen 200 49 122

Schon bei dieser Berechnung ist die Konvergenzrate der Levenberg-Marquardt-Methode deutlich besser als die der Quasi-Newton-Methode. Die Gauss-Newton-Methode ist weiterhin mit Abstand die bessere Methode. Auch der Unterschied zwischen den Ergebnissen, die mit der numerischen und denen, die mit der analytischen Sensitivit¨ atsanalyse erzielt wurden, ist relativ gross.

In der Abbildung 6.13 sind die durch die Optimierungsverfahren gefundenen Lasten dargestellt

Auch hier sieht man die optimale Dichteverteilung sehr sch¨ on, die aus der unsymmetrischen Belastung resultiert. In Richtung der Last B wird mehr Dichte angebaut, da die dazugeh¨ orige

51

Last gr¨ osser ist.

6.3.4 Beispiel 4

Gegeben ist das folgende 3D-System

wobei die Lasten zun¨ achst als Platzhalter dargestellt wurden. Erst durch die Berechnungen werden Gr¨ osse und Richtung festgelegt. Der Abb. 6.14 kann entnommen werden, dass die Belastung an der Stelle A gr¨ osser als die an der Stelle B sein sollte, da in der Umgebung von A mehr Dichte angebaut wird, als in der Umgebung von B. Weiter wird auch hier klar, dass in den Stellen, wo die Spannungen sehr klein zu erwarten sind, die Dichte nahezu vollst¨ andig abgebaut wird (blauer Bereich). Die Berechnungen wurden mit folgenden Startwerten durchgef¨ uhrt: (50, 110, 60, −110, −60, −115) ^T .

52

analytisch

Anzahl Iterationen 34 8 9

Funktionsauswertungen 308 79 102

Zeit [s] 3737 959 1238

10 ⁻¹² 10 ⁻²⁰ 10 ⁻¹⁶ Residuum

numerisch

Anzahl Iterationen 34 8 9

Funktionsauswertungen 308 79 102

Alle drei Verfahren lieferten die richtige L¨ osung (100, 100, 50, −100, −50, −100) ^T . Der Unterschied zwischen den numerischen und den analytischen Berechnungen wird mit diesem 3D-Beispiel sehr deutlich. Die Rechenzeiten der Guass-Newton- und der Levenberg-Marquard-Methode sind leicht unterschiedlich. Dagegen hat die Quasi-Newton-Methode viel mehr Zeit f¨ ur die Minimierung der Funktion ben¨ otigt.

6.4 Diskussion der numerischen Ergebnisse aus 2D -und 3D-Berechnungen

In einer zusamenfassenden Bewertung k¨ onnen folgende Aussagen getroffen werden: F¨ ur kleine nichtlineare Ausgleichsprobleme, etwa wie die ersten Beispiele dieses Kapitels, sind die drei vorgestellen Gradienten-Verfahren gleichwertig. Es gab Startwerte, die f¨ ur eine Methode geeigneter waren und schnell zum Ziel gef¨ uhrt haben, und es gab auch solche die weniger geeignet waren. F¨ ur gr¨ ossere Probleme wie z.B. das 3D-Beispiel wurden die Vorteile der Gauss-Newton-und Levenberg-Marquardt-Methode f¨ ur die Behandlung nichtlinearer Ausgleichsprobleme sehr deutlich. Auch die Vorteile der analytischen Sensitivit¨ atsanalyse werden bei gr¨ osseren Problemen und insbesondere bei vielen Designvariablen bemerkbar.

53

Kapitel 7

Anwendung am Femur

Mit den gewonnenen Erkenntnissen aus den vorangegangenen Kapiteln, und insbesondere des letzten Kapitels, sollen die im Rahmen dieser Arbeit durchgef¨ uhrten Berechnungen am sogennanten “Standardized Femur“, [28], diskutiert werden. Als Mass f¨ ur die G¨ ute der mittels Optimierungsverfahren erzielten Resultate, soll die durch CT-Aufnahmen erfasste Dichteverteilung dienen. In Abb. 7.1 ist der Femur mit einem Netz von 3080 Elementen und der optimalen Dichteverteilung dargestellt.

Die Idee besteht darin, aus dieser Dichteverteilung auf die einwirkenden Lasten zur¨ uckzuschliessen. Dabei werden Berechnungen mit zwei, drei -und f¨ unf Kr¨ aften durchgef¨ urt. Weiter stehen die ver¨ offentlichten Messwerte zur Verf¨ ugung, die beim normalen Gehen gemessen wurden, [25], [26], [27].

54

Messwerte f¨ ur zwei Kr¨ afte

Messwerte f¨ ur drei Kr¨ afte

Messwerte f¨ ur f¨ unf Kr¨ afte

Der oberen Tabelle kann man entnehmen, dass bei einer Erh¨ ohung der Anzahl der Kr¨ afte, die bestehenden weiterhin benutzt werden. Diese Kr¨ afte werden als Startwerte f¨ ur die Gradienten-Verfahren benutzt. Die 2D- und 3D-Beispiele des letzten Kapitels haben gezeigt, dass es nicht einfach ist, auf die unbekannten Belastungen zur¨ uckzuschliessen, da einerseits viele Modellannahmen vorhanden sind, und andererseits die zu minimierende Funktion mehrere lokale Minima aufweist. In den n¨ achsten Abschnitten werden die Ergebnisse, die mittels des im Kapitel 4 und 6 vorgestellten algorithmischen L¨ osungskonzeptes ermittelt wurden, vorgestellt und diskutiert. Zun¨ achst werden die Ergebnisse betrachtet, die nur mit den Gradienten-Verfahren ermittelt wurden, anschliessend werden auch Ergebnisse vorgestellt die durch die Kopplung des Genetischen - Algorithmus und den Gradienten-Verfahren ermittelt wurden.

55

In Abb. 7.2 sind die gemessenen Werte dargestellt. Auch f¨ ur die weiteren Berechnungen werden die Lastangriffpunke nach wie vor die gleichen bleiben. Die Verschiebung aller Knoten am unteren Ende des Femurs wird in allen drei Richtungen verhindert.

7.1 Berechnungen mit zwei Kr¨ aften

Mit den Gradienten-Verfahren und den Startwert (F 1 , F 2 ) ^T wurden die folgenden Ergebnisse erzielt.

analytisch

Anzahl Iterationen 86 11 14

Funktionsauswertungen 92 122 182

Zeit [h] 28 40 22.3

Die ermittelten Kr¨ afte sind recht unterschiedlich, obwohl das Residuum bei allen Methoden nahezu gleich gross ist. Das deutet darauf hin, dass es sich um drei verschiedene lokale Minima handeln k¨ onnte. Es w¨ are auch denkbar, dass es sich um die gleiche Minimalstelle handelt, die aber sozusagen langgezogen ist bzw. ¨ uber einen grossen Bereich einen sehr flachen Verlauf

haben k¨ onnte. Eine solche Form der Minimalstelle kann auf unterschiedliche Kr¨ afte mit etwa den gleichen Funktionswert f¨ uhren. Bei den ersten zwei Verfahren ist nach wie vor die Gelenkkraft F 2 gr¨ osser als die Muskelkraft F 1 , was auch sein sollte. Im Gegensatzt dazu wird mit der Levenberg-Marquardt-Methode eine kleinere Gelenkkraft ermittelt. Die folgenden Bilder sollen die berechneten Dichtewerte grafisch veranschaulichen.

Obwohl das Residuum relativ gross ist und nur mit zwei Kr¨ aften bzw. sechts Designvariablen

57

gerechnet wurde, ist eine sehr gute ¨ Ubereinstimmung mit der optimalen Dichteverteilung aus Abb. 7.1 festszustellen.

Auch die Ergebnisse der Gauss-Newton-Methode stimmen sehr gut mit den optimalen Dichtewerten ¨ uberein.

Bild 7.5: Dichteverteilung: Levenberg-Marquardt-Verfahren, 6 Designvariablen

Der Dichteverlauf, der in Abb. 7.5 dargestellt ist, ist qualitativ nicht von den zwei vorangegangenen Dichteverteilungen zu unterscheiden, d.h. auch das Levenberg-Marquardt-Verfahren liefert Ergebnisse, die in sehr guter ¨ Ubereinstimmung mit den optimalen Werten sind.

58

Mit den gemessenen Kr¨ aften ist das Residuum 3454 gross. Mit den Optimierungsverfahren wurde es um etwa 1000 reduziert. Das Ziel ist es nun, das Residuum noch weiter zu reduzieren. Der Tabelle 7.2 kann man entnehmen, dass das Levenberg-Marquardt-Verfahren die beste Konvergenzgeschwindigkeit hatte. Erstaunlicherweise ist in diesem Fall mit dem Quasi-Newton-Verfahren sowohl eine schnellere Konvergenz als auch ein kleineres Residuum als beim Gauss-Newton-Verfahren erreicht worden. In allen Abbildungen ist die R¨ ohrenform des Knochens sehr gut zu erkennen. Auf dem Rand ist die Dichte etwas gr¨ osser und im inneren etwas kleiner als bei der optimalen Dichteverteilung in Abb. 7.1. Weiter sieht man in der Umgebung der Lastangriffspunkte die lokal hohen Dichtewerte, die aufgrund der konzentrierenten bzw. Einzellasten zustande kommen.

7.2 Berechnungen mit drei Kr¨ aften

Mit den Gradienten-Verfahren und den Startwert (F 1 , F 2 , F 3 ) ^T wurden die folgenden Ergebnisse erzielt.

analytisch

Anzahl Iterationen 106 11 8

Funktionsauswertungen 164 147 149

Zeit [h] 42.8 25.2 24.7

Es ist offensichtlich, dass die Konvergenzgeschwindigkeit der Quasi-Newton-Methode deutlich kleiner als die der Gauss-Newton- und Levenberg-Marquardt-Methode ist. Das Residuum ist bei der Gauss-Newton-Methode kleiner als bei den anderen zwei Methoden. In Tabelle 7.5 sind die Lasten, die durch die drei Optimierungsverfahren ermittelt wurden, dargestellt. Die wichtige Erkenntnis aus der Tabelledie 7.5, dass die Gelenkkraft F 2 nach wie vor eine Drucklast ist. Daraus kann man folgern, dass die ermittelten Lasten nicht in Widerspruch zu den physikalischen Gesetzen stehen. Diese Ergebnisse sind aus mathematischer Sicht auch ohne diese ¨ Uberlegungen korrekt und f¨ uhren zu den drei gefundenen Minimalstellen.

59

Auch in diesem Fall sind die berechneten Kr¨ afte unterschiedlich gross und haben unterschiedliche Vorzeichen. Beim Gauss-Newton-Verfahren ist die Muskelkraft F 3 relativ gross. Die Ergebnisse, die durch das Quasi-Newton- und das Levenberg-Marquardt-Verfahren erzielt wurden sind plausibel, wobei die Kr¨ afte die durch das Quasi-Newton-Verfahren ermittelt wurden verh¨ altnism¨ assig eine bessere ¨ Ubereinstimmung mit den dynamischen Messwerten liefern.

Da die Gesamtkraft auf drei Lasten verteilt wird und somit lokal kleinere Spannungen vorhanden sind, ist auch die Dichte in der Umgebung der Muskelkraft F 1 nicht mehr so hoch wie bei der Belastung mit nur zwei Kr¨ aften. Die “Wandst¨ arke“ des Schaftes hat im Vergleich zu der

60

Belastung mit nur zwei Kr¨ aften zugenommen.

Beim Gauss-Newton-Verfahren ist das Residuum kleiner als bei den anderen zwei Verfahren und die Dichteverteilung ist im Schaftbereich ebenfalls besser. Im proximalen ¹ Bereich sieht aber die Dichteverteilung weniger gut aus. Grund daf¨ ur k¨ onnte die sehr grosse Muskelkraft F 3 sein. Auch die Gelenkkraft F 2 ist relativ gross und kann als Ursache f¨ ur die lokal hohen Dichtewerte betrachtet werden.

¹ Oberst¨ uck von R¨ ohrenknochen

61

Bild 7.8: Dichteverteilung: Levenberg-Marquardt-Verfahren, 9 Designvariablen

Mit der Levenberg-Marquardt-Methode wurde qualitativ die gleiche Dichteverteilung wie bei der Quasi-Newton-Methode ermittelt. Die R¨ ohrenform des Femurs ist nach wie vor in allen Berechnungen sehr gut zu erkennen.

7.3 Berechnungen mit f ¨ unf Kr¨ aften

F¨ ur die Gradienten-Verfahren wurden als Startwert die gemessenen Kr¨ afte (F 1 , F 2 , F 3 , F 4 , F 5 ) ^T genommen und die folgenden Resultate erzielt

analytisch

Anzahl Iterationen 108 9 13

Funktionsauswertungen 183 211 311

Zeit [h] 65.2 52.3 72.3

Offensichtlich hat die Levenberg-Marquardt-Methode f¨ ur diese Berechnung mehr Zeit ben¨ otigt, daf¨ ur wurde aber das Residuum deutlich st¨ arker reduziert, als bei der Quasi-Newton- und der Gauss-Newton-Methode. H¨ atte man die Berechnung bei der Levenberg-Marquardt-Methode bei einem Residuum von 2250 abgebrochen, so h¨ atte die Berechnung nur 48h gedauert. Ob wirklich

62

alleine das Residuum f¨ ur die G¨ ute der Ergebnisse betrachtet werden sollte, wird wieder andahnd der Dichteverteilung in den n¨ achsten Abbildungen festgestellt . Obwohl das Residum der ersten zwei Methoden nahezu gleich ist, sind die ermittelten Kr¨ afte Tab. 7.7 sehr unterschiedlich. Auch hier k¨ onnte der glieche Fall wie in den vorigen Berechnungen vorliegen, n¨ amlich mehrere lokale Minima oder ein sehr flacher Verlauf der Minimalstelle in der N¨ ahe der L¨ osung. Zun¨ achst werden die mittels der erw¨ ahnten Methoden gewonnenen Ergebnisse tabelarisch aufgelistet und diskutiert.

Es ist weiterhin so, dass die Gelenkkraft F 2 , physikalisch richtig, eine Drucklast ist. Weiter kann festgestellt werden, das die Gelenkkraft bei allen Methoden etwa die gr¨ osste Last ist. Anders sieht es bei F 1 und F 5 aus. Mit den letzten zwei Verfahren wurde eine zu kleine Muskelkraft F 1 und eine zu grosse Muskelkraft F 5 ermittelt. Aus diesem Grund k¨ onnen die Ergebnisse des Quasi-Newton-Verfahrens als repr¨ asentative Kr¨ afte betrachtet werden. In den Abb. 7.9, 7.10 und 7.11 ist der Dichteverlauf grafisch dargestellt.

63

Bild 7.9: Dichteverteilung: Quasi-Newton-Verfahren, 15 Designvariablen

Auch diese Dichteverteilung stimmt sehr gut mit der optimalen Dichteverteilung ¨ uberein. Insbesondere die R¨ ohrenform des Femurs ist weiterhin deutlich zu sehen. Aufgrund der relativ hohen Gelenkkraft F 2 sind auch bei dieser Berechnung die lokal hohen Dichtewerte entstanden.

Bild 7.10: Dichteverteilung: Gauss-Newton-Verfahren, 15 Designvariablen

Bei der Gauss-Newton-Methode ist im Schaftbereich eine bessere ¨ Ubereinstimmung mit der

optimalen Dichteverteilung festzustellen, da im Randbereich weniger Dichte als bei der ersten Methode angebaut wurde. Auch im proximalen Bereich sind die lokal hohen Dichtewerte in der Umgebung der Muskelkraft F 1 nur sehr schwach vorhanden, daf¨ ur sind sie in der Umgebung

64

der recht grossen Gelenkkraft F 2 st¨ arker vertreten. Es sind alo bei etwa gleichem Residuum recht unterschiedliche Dichteverl¨ aufe zu beobachten.

Bild 7.11: Dichteverteilung: Levenberg-Marquardt-Verfahren, 15 Designvariablen

Bei der Levenberg-Marquardt-Methode wird der Unterschied zu den ersten zwei Methoden deutlich gr¨ osser. Die lokal hohen Dichtewerte in der Umgebung der Gelenkkraft F 2 sind nicht mehr zu beobachten. Daf¨ ur wird aber in der Umgebung der Muskelkraft F 1 recht stark Dichte angebaut. Grund daf¨ ur k¨ onnte die relativ hohe, negative z-Komponente der Muskelkraft F 5 sein. Eine geeignetere Erkl¨ arung hierf¨ ur konnte nicht gefunden werden. Das Levenberg-Marquardt-Verfahren hatte das Residuum am st¨ arksten reduziert und weist im proximalen Bereich dennoch eine andere Dichteverteilung auf, als das bei der optimalen Dichtevteilung der Fall ist. Dieses Beispiel zeigt, das auch bei einem kleineren Residuum lokal gr¨ ossere Abweichungen als bei einem gr¨ osseren Residuum auftretten k¨ onnen. Ein ¨ anliches Verhalten wurde beim Gauss-Newton-Verfahren Abb. 7.7 mit 9 Designvariablen beobachtet. Ausser im erw¨ ahnten Bereich, lieferte die Levenberg-Marquardt-Methode einen sehr guten Dichteverlauf, der keine ausgepr¨ agte hohen Dichtewerte aufweist.

65

7.4 Berechnungen mit f ¨ unf Kr¨ aften und den Genetischen-Algorithmus als Vorrechnung

Die folgenden Ergebnisse sind durch die Kopplung des Genetischen-Algorithmus und der Gradienten-Verfahren ermittelt worden. F¨ ur 15 Designvariablen wurde mit dem Gauss-Newton-Verfahren das Residuum von 3603 auf 2284 reduziert und ein lokales Minimum gefunden, Tab. 7.7. F¨ ur den Genetischen Algorithmus wurde ein neuer Suchbereich ±15% um den L¨ osungsvektor des Gauss-Newton-Verfahrens definiert. Der Genetische-Algorithmus lieferte nach 40h Rechenzeit ein Residuum von 1776 und den folgenden Lastvektor:

Dieser Lastvektor wurde wieder als Startwert f¨ ur den Gauss-Newton-Verfahren benutzt, das das Residuum nach 41h und 174 Funktionsauswertungen auf 1700 reduzierte

Es wird festgstellt, dass das Residuum nicht mehr wie am Anfang der Berechnugen stark reduziert wird. Da die Gelenkkraft F 2 , insbesondere die z-Komponente sehr klein ist, wird davon aussgegangen, dass obwohl das Residuum im Vergleich zu den anderen Berechnungen klein ist, lokal hohe Abweichungen vorhanden sein m¨ ussen. F¨ ur diese Resultate wurde die Dichteverteilung in Abb. 7.12 grafisch dargestellt. Mit dem Genetischen-Algorithmus war es m¨ oglich

66

das lokale Minimum zu ¨ uberspringen und weiterzurechnen. Es scheint aber weiterhin so zu

sein, dass je kleiner das Residuum wird, um so ausgepr¨ agter ist die relativ grosse Dichte im proximalen Bereich, Abb. 7.12 (gr¨ unner Bereich).

Bild 7.12: Dichteverteilung: Gauss-Newton-Verfahren, 15 Designvariablen und GA als Vorrech-

nung

Die lokal hohen Dichteabweichungen sind sowohl in der Umgebung der Muskelkraft F 1 als auch in der Umgebung der Gelenkkraft F 2 vollst¨ andig beseitigt. Es wurden diverse andere Berechnungen durchgef¨ uhrt die bis auf ein Residuum von 1615 gef¨ uhrt haben. Es waren aber keine neuen Erscheinungen zu beobachten. Aus diesem Grund wird nur noch der Dichteverlauf, der mit der Kopplung des Levenberg-Marquardt-Verfahrens und des Genetischen-Algorithmus ermittelt wurde, grafisch dargestellt

67

Bild 7.13: Dichteverteilung: Levenberg-Marquardt-Verfahren, 15 Designvariablen und GA als

Vorrechnung

¨ Ahnlich wie bei der Gauss-Newton- Berechnung ist auch mit dem Levenberg-Marquardt-Methode ein sehr ¨ ahnlicher Dichteverlauf ermittelt worden. Nach wie vor, sind die ermittelten Kr¨ afte sehr unterschiedlich. Da die Berechnungen sehr lange dauern und im Rahmen dieser Arbeit die Zeit sehr begrenzt ist, muss man sich mit den ermittelten Ergebnissen zufrieden stellen. Im n¨ achsten Kapitel werden die wichtigsten Erkenntnisse nochmal zusammengefasst.

68

Kapitel 8

Diskussion der Ergebnisse

In einer zusammenfassenden Bewertung der Ergebnisse des vorigen Kapitels k¨ onnen die folgenden Aussagen getroffen werden: Die Vergleiche der drei gradienten Verfahren haben gezeigt, dass sowohl das Gauss-Newton- als auch das Levenberg-Marquardt-Verfahren dem Quasi-Newton-Verfahren bzgl. der Konvergenzgeschwindigkeit f¨ ur die behandelten Problemstellungen, ¨ uberlegen sind. Zwischen den ersten zwei k¨ onnen bzgl. der Vor- und Nachteile keine aussagekr¨ aftigen Bemerkungen gemacht werden, da beide Methoden f¨ ur diese Aufgabe geeignet waren. Die Effizienz der Verfahren wurde von den Startwerten sehr beeinflusst. F¨ ur die gleichen Startwerte wurden unterschiedliche L¨ osungen gefunden, was auf die unterschiedlichen Suchrichtungen der Verfahren zur¨ uckzuf¨ uhren ist. Sehr oft war die Gauss-Newton-Methode schneller als die Levenberg-Marquardt-Methode. Grund daf¨ ur k¨ onnte eine schwache Nichtlinearit¨ at zwischen Start- und Endpunkt der Zielfunktion gewesen sein. Im vorigen Kapitel wurde aufgrund der grossen Rechenzeiten ausschliesslich die analytische Sensitivit¨ atsanalyse benutzt. Um den Unterschied zwischen der numerischen und der analytischen Sensitivit¨ atsanalyse deutlich zu machen, wurde noch eine Berechnung mit der numerischen Sensitivit¨ atsanalyse und dem Gauss-Newton-Verfahren mit neun Designvariablen durchgef¨ uhrt und die Rechenzeiten miteinander verglichen. Man hat festgestellt das die Rechenzeit der Berechnung mit der numerischen Sensitivit¨ atsanalyse etwa das dreifache der analytischen Variante betrug. Diese Berechnung wurde im vorigen Kapitel nicht vorgestellt. Berechnungen, nur mit dem Genetischen-Algorithmus, wurden f¨ ur diese grosse Problemstellung nicht durchgef¨ uhrt, da die Konvergezgeschwindigkeit sehr klein ist. Weiter konnte festgestellt werden, dass mit der Kopplung des Genetischen-Algorithmus und den Gradienten-Verfahren das Residuum um einen wesentlichen Betrag reduziert werden konnte.

69

Die durch den entwickelten Algorithmus erzielten Ergebnisse stimmen mit den gemessenen CT-Daten sehr gut ¨ uberein. Weiter konnte festgestellt werden, dass auch die Berechnungen mit nur zwei Lasten gute Resultate lieferten. Die folgende Abbildung soll einige Vergleiche zusammenstellen

Auch in diesem Vergleich ist die sehr gute ¨ Ubereinstimmung zwischen den berechneten und der mittels CT-Aufnahmen ermittelte optimale Dichteverteilung, gut zu sehen.

70

Kapitel 9

Zusammenfassung und Ausblick

Mit dem entwickelten algorithmischen L¨ osungskonzept war es m¨ oglich, statisch ¨ aquivalente

Lastf¨ alle zu finden, die die Berechnung der Dichteverteilung des Femurs in guter ¨ Ubereinstimmung, mit der aus CT-Aufnahmen erfassten Dichteverteilung erm¨ oglichen. Zur Zeit sind es unterschiedliche berechnete Lastf¨ alle, die auf die gleichen Ergebnisse f¨ uhren k¨ onnen, und so die zielgerichtete Berechnung der statisch ¨ aquivalenten Lasten in Frage stellen.

Eine M¨ oglichkeit dem gesetzten Ziel n¨ aher zu kommen, w¨ are eine Gl¨ attung der Zielfunktion. Dadurch w¨ urden die unz¨ ahlingen lokalen Minima verschwinden und so die Berechnungen mit den Gradienten-Verfahren effizienter machen. Auch eine Reduktion der Modellannahemen k¨ onnte zu einer zielgerichteten Berechnung f¨ uhren. Insbesondere die Referenzverzerrungsenergiedichte ψ ref als Designvariable definiert und in den Optimierungsprozess integriert, k¨ onnte einen wesentlichen Beitrag dazu leisten. Die Ergebnisse der Berechnungen im Kapitel 7 haben gezeigt, dass sogar f¨ ur kleinere Residuen, lokal h¨ ohere Abweichungen auftretten k¨ onnen. Das liegt daran, dass das Residuum eine skalare Gr¨ osse ist, die keine lokalen Aussagen liefern kann. Auch eine andere Definition der Zielfunktion k¨ onnte zu besser kontrollierbaren Berechnungen f¨ uhren. Eine weitere Verbesserungsm¨ oglichkeit liegt darin, die wesentlichen Lasten wie die Gelenkkraft F 2 und die Muskelkraft F 1 mit einer h¨ oheren Gewichtung in die Berechnung einzubeziehen. Man konnte feststellen, dass auch mit zwei bzw. drei Kr¨ aften sehr gute Ergebnisse erzielt werden konnten. Aus diesem Grund ist es naheliegend, die Berechnungen auf die dominierenden Lasten F 1 , F 2 und F 3 (siehe Abb. 7.2) zu beschr¨ anken. Um lokal hohe Beanspruchnungen zu vermeiden, k¨ onnten die Kr¨ afte als Fl¨ achenlasten definiert werden.

71

Zusammengefasst kann gesagt werden, dass mehr Aufwand in den oben angesprochenen Verbesserungsm¨ oglichkeiten ben¨ otigt wird, um zu eindeutigen ¨ aquivalenten statischen Lasten zu

kommen. Weiter k¨ onnen neue Erkenntnisse aus dem Bereich des Knochenumbauprozesses zu einer wesentlich zielgerichteten Behandlung des inversen Problems beitragen.

72

Literaturverzeichnis

[1] W. ALT: Nichtlineare Optimierung, 1. Auflage, Vieweg 2002

[2] E. POLAK: Optimization: algorithms and consistent approximations, Springer 1997

[3] C. GEIGER & Ch. KANZOW: Numerische Verfahren zur L¨ osung unrestringierter Optimierungsaufgaben, Springer 1999

[4] G. STARKE: Numerik nichtlinearer Optimierung, Institut f¨ ur angewandte Mathematik, Universit¨ at Hannover, Vorlesungsunterlagen WS 03/04

[5] F. JARRE & J. STOER: Optimierung, Springer 2003

[6] P. KOSMOL: Methoden zur numerischen Behandlung nichtlinearer Gleichungen und Optimierungsaufgaben, Teubneer-Verlag, 1989

[7] H. -P. SCHWEFEL: Numerical Optimization of Computer Models, John Wiley & Sons,

1981

[8] K.-J. BATHE: Finite-Elemente-Methoden, Aus dem Englischen ¨ ubersetzt von Peter Zimmermann, 2. Auflage Springer 2002

[9] D. BRAESS: Finite-Elemente, Theorie, schnelle L¨ oser und Anwendungen in der Elastizit¨ atstheorie, 3. Auflage Springer 2002

[10] L. MEYER: Formoptimierung in der Strukturdynamik, Forschungs und Seminarberichte aus dem Bereich der Mechanik der Universit¨ at Hannover, 1998

[11] F.-J. BARTHOLD: Theorie und Numerik zur Berechnung und Optimierung von Strukturen aus isotropen, hyperelastischen Materialien, Forschungs und Seminarberichte aus dem Bereich der Mechanik der Universit¨ at Hannover, 1993

73

[12] S. SCHWARZ: Sensitivit¨ atsanalyse und Optimierung bei nichtlinearem Strukturverhalten, Bericht Nr. 34, Institut f¨ ur Baustatik der Universit¨ at Stuttgart , 2001

[13] S. KIMMICH: Strukturoptimierung und Sensibilit¨ atsanalyse mit finiten Elementen, Dissertation, Institut f¨ ur Baustatik der Universit¨ at Stuttgart , 1990

[14] R. MAHNKEN: Theoretische und numerische Aspekte zur Parameteridentifikation und Modellierung bei metallischen Werkstoffen, Forschungs und Seminarberichte aus dem Bereich der Mechanik der Universit¨ at Hannover, 1998

[15] K. WIECHMANN: Theorie und Numerik zur Berechnung und Optimierung von Strukturen mit elastoplastischen Deformationen, Dissertation, Institut f¨ ur Baumechanik und Numerische Mechanik der Universit¨ at Hannover, 2001

[16] R.T. HAFTKA, Z. G ¨ URDAL and M.P. KAMAT: Elements of Structural Optimization, Kluwer Academic Publishers, 1990

[17] E. J. HAUG & K. CHOI & V. KOMKOV: Design Sensitivity Analysis of Structural Systems, Mathematics in Science and Engineering, Volume 177, Academic Press, 1986

[18] A.G. RAMM: Inverse Problems, mathematical and analytical techniques with applications to engineering, Springer, 2005

[19] U. NACKENHORST & U. SCHR ¨ ODER: Numerische Simulation des beanspruchungsa-

daptiven Knochenumbaus endoprothetisch versorger Femura, Bericht aus dem Institut f¨ ur Mechanik, Universit¨ at der Bundeswehr Hamburg , 1996

[20] U. NACKENHORST, N. KRISTIN & R. LAMMERING: Zur konstitutiven Beschreibung des anisotropen beanspruchnungsadaptiven Knochenumbaus, Technische Mechanik, Vol. 20, S. 31-40, 2000

[21] U. NACKENHORST: Numerical simulation of stress simulated bone remodeling, Technische Mechanik, Vol. 17, S. 31-40, 1997

[22] B. EBBECKE and U. NACKENHORST: Numerical Studies on the Biocompatibility of Innovative Hip-Joint-Prosthesis, The Finite Element Method in Biomedical Engineering, Biomechanics and Related Fiels Ulm, 2003

74

[23] B. EBBECKE and U. NACKENHORST: Simulation of Stress adaptive Bone Remodelling - Towards an individual Therapy in Endoprosthetics, PAMM-Proc in Appl. Mathem. and Mechanics, Dresden, 2004

[24] S. WENG: Ein Evolutionsmodell zur mechanischen Analyse biologischer Struktuturen, Mitteilungen aus dem Institut f¨ ur Mechanik, Ruhr Universit¨ at Bochum

[25] G. N. DUDA, M. HELLER and G. BERGMANN: Musculoskeletal loading database: loading conditions of the proximal femur, Theoretical Issues in Ergonomics Science, vol. 6, p. 287-292, 2005

[26] M. O. HELLER, G. BERGMANN, J. P. KASSI, L. CLAES, N. P. HAAS and G. N. DU-DA: Determination of muscle loading at the hip joint for use in pre-clinical testing, Journal of Biomechanics, vol. 38, p. 1155-1163, 2005

[27] M. O. HELLER, G. BERGMANN, G. DEURETZBACHER, L. D ¨ URSELEN, M. POHL,

L. CLAES, N. P. HAAS and G. N. DUDA: Musculo-skeletal loading conditions at the hipduring walking and stair climbing, Journal of Biomechanics, vol. 34, p. 883-893, 2001

[28] T. FULVIA, A. PANCANTI and M. VICECONTI: An improved method for the automatic mapping of computet tomography numbers onto finite element models, Medical engineering and physics 26, 2004

[29] J. WOLFF : Das Gesetz der Transformation der Knochen, Berlin, 1892

75