Universit¨at Hannover

Institute of Mechanics and Computational Mechanics

Inverse numerische Simulation zur

Berechnung statisch ¨aquivalenter Lasten aus

CT-Aufnahmen des menschlichen Femurs

Bekim Berisha

---

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 Inge-

nieuren 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 Vorrech-

nung .

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 Vor-

rechnung

.

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¨afte-

berechnung 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 Proble-

me 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¨uck-

gef¨uhrt. Die Idee besteht darin, aus einer bekannten bzw. mittels CT1-Aufnahmen erfassten

optimalen Dichteverteilung die Lasten zu ermitteln. Dadurch k¨onnen z.B. individuelle Thera-

pien 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 Problem-

stellung 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 gradientenbasier-

ten Verfahren im Vordergrund stehen. Um mit Gradienten-Verfahren die resultierende Optimie-

rungsaufgabe 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.

1Abk¨urzung f¨ur Computer Tomographie

9

---

Kapitel 2

Knochenumbau

In diesem Kapitel werden die wichtigsten Aspekte des Knochenumbauprozesses und der Mo-

dellbildung 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 Elemen-

te 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

der ¨ausseren Gestalt in ¨

Ubereinstimmung mit mathematischen Gesetzten. Die Untersuchung des

Knochenumbauprozesses ist nicht Gegenstand dieser Arbeit. Aus diesem Grund wird hier nicht

auf Einzelheiten eingegangen. Interessierte Leser werden besonders auf die Ver¨offentlichun-

gen 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 Pro-

blemstellungen von der folgenden Aussage ausgegangen: Knochen sind in ihrer ¨ausseren Ge-

stalt 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¨ohren-

knochen. 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

E

n

=

(2.1)

E0

0

wobei E0 ein mittlerer Elastizit¨atsmodul bzw. 0 eine mittlere Dichte zu Beginn des Umbaupro-

zesses 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

bzw.

0 E Ec

(2.2)

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

=

(2.3)

0

aus (2.1) und (2.3) folgt

E = E0n

(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

= nC0

(2.6)

f¨ur die Verzerrungsenergiedichte gilt folgendes

1

1

1

=

U =

T =

T C

(2.7)

2

2

wobei U die Verzerrungsenergie und C der vierstufige Materialtensor ist.

(2.6) in (2.7) eingesetzt liefert

1

=

nT C0

(2.8)

2

Dass der Exponent n = 2 sein muss, kann durch eine konsistente thermodynamische Betrach-

tung 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 Knochen-

substanz ausdr¨uckt, Nackenhorst[21].

= k

- 1

(2.9)

ref

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 ¨uber-

nommen und k = 0 gesetzt. Dadurch resultiert die folgende Evolutionsgleichung

=

- 1

(2.10)

ref

Ist gr¨osser als der Referenzwert ref , so kommt es zum Anbau von Knochenmasse. Ist al-

lerdings < ref , so wird Knochenmasse abgebaut. Die Evolutionsgleichung (2.10) wird mit

Hilfe der Finiten-Differenzen-Methode gel¨ost:

=

- 1

(2.11)

t

ref

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

t

t+t =

- 1 t

(2.12)

ref

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 Elemen-

te 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 im-

mer 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¨alti-

ger Form in der Literatur dargestellt, siehe z.B. (Bathe[8], Braess[9]). Aus diesem Grund wer-

den 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 Semidiskretisie-

rung 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 Feldglei-

chungen sowie die konstitutive Annahme ¨uber das Werkstoffgesetz ber¨ucksichtigt werden.

Statische Feldgleichung

Div + b = 0

(3.1)

Konstitutive Beziehung

= C :

(3.2)

Kinematische Feldgleichung

1

=

Gradu + GradT u

(3.3)

2

Zu diesem Satz von Gleichungen kommen noch die Randbedingungen hinzu. Diese sind die

Spannungsrandbedingungen auf dem Spannungsrand t und die Verschiebungsrandbedingun-

gen auf dem Verschiebungsrand u

u = ¯

u

auf

u

(3.4)

n = t = t

auf

t

(3.5)

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 Gleichge-

wichts erforderlich sind. Durch Multiplikation der Ausgangsgleichung (3.1) mit einer zul¨assi-

gen Testfunktion , die die wesentlichen Randbedingungen erf¨ullt, und Integration ¨uber das

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

(Div + b) dV = 0

(3.6)

14

---

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

t - t dA = 0

(3.7)

t

(3.6) und (3.7) ergeben

(Div + b) dV +

t - t dA = 0

(3.8)

t

Das Ausmultiplizieren und die Anwendung von Rechenregeln1 zur Divergenzbildung eines Fel-

des auf den ersten Integranden liefert

Div () - : Grad + b dV +

t - t dA

(3.9)

t

Die Anwendung des Divergenztheorems2 auf den ersten Term f¨uhrt zu

() n dA -

: Grad dV +

b dV +

t - t dA = 0

(3.10)

t

t

Umformung des ersten und des letzten Terms liefert

(n) dA -

: Grad dV +

b dV +

t dA -

t dA = 0

(3.11)

t

t

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.

: Grad dV -

b dV -

t dA = 0

(3.12)

t

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

1hier 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

2Divergenztheorem:

Div (...) dV =

(...) n dA

t

15

---

G (u, ) =

Grad : C : Gradu dV -

b dV -

t dA = 0

(3.13)

t

Da in der schwachen Form des Gleichgewichts (3.13) nur erste Ableitungen auftreten, muss f¨ur

das Verschiebungsfeld u die C0-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 Randbedingun-

gen 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 Ele-

menten 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 Randbedingun-

gen 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. Ausge-

hend 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 Pa-

rameter¨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 Ziel-

funktion bezeichnet, die wie folgt definiert wird

1

f =

- ^

2

(4.1)

2

Auf Minimierung der Zielfunktion wird im Kapitel Optimierungsverfahren" n¨aher eingegan-

"

gen. Man k¨onnte sich auch eine andere Funktion definieren z.B. die Summe der absoluten Ab-

weichungen. 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 Ziel-

funktion, 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 nu-

merischen L¨osung von Optimierungsaufgaben mittels Gradienten-Verfahren darstellt. Ihre Qua-

lit¨at beeinflusst sehr stark die Rechenzeit. Aus diesem Grund ist eine m¨oglichst effiziente Vari-

ante 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 nichtli-

neare Probleme behandeln.

4.1

Numerische Sensitivit¨atsanalyse

Die bekanntesten und einfachsten Verfahren zur Implementierung der numerischen Sensiti-

vit¨atsanalyse sind das Euler-Vorw¨arts- und das zentrale Differenzen-Verfahren. Wegen der ein-

fachen Implementierung und dem geringen Formulierungsaufwand sind diese Verfahren in der

Praxis sehr beliebt und finden eine breite Anwendung, auch ausserhalb der numerischen Sensiti-

vit¨atsanalyse . Die erw¨ahnten Verfahren k¨onnen durch eine abgebrochene Taylorreihenentwick-

lung 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 Genau-

igkeit 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 nu-

merischen 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

---

Algorithmus 1

initialisiere

for i = 1 : Ns

s(i) = s(i) +

for j = 1 : Nt

FE - Berechnung

durchf ¨

uhren

end

berechne die ¨

Anderung der Dichteverteilung

ds = neu-alt

Auswertung des Gradienten der Zielfunktion

df (i) = ds(alt - ^

)

end

wobei das Inkrement (die St¨orung) der Designvariablen, Ns die Anzahl der Designva-

riablen und Nt die Anzahl der Zeitschritte sind

Man stellt fest, dass f¨ur die Umsetzung der numerischen Sensitivit¨atsanalyse lediglich die Funk-

tionswerte 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¨atz-

liche Anzahl FE-Berechnungen durchgef¨uhrt werden muss (Zeitabh¨angiges Problem: Evoluti-

onsgleichung 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 un-

erw¨unschte Eigenschaften wie Rundungs- und Konditionsfehler auftretten. W¨ahlt man jedoch

die St¨orung zu gross, so treten Ungenauigkeiten aufgrund der abgebrochenen Taylorreihen-

entwicklung auf. Auf die Einzelheiten wird nicht weiter eingegangen und insbesondere auf die

Arbeit von Kimmich [13] verwiesen. Obwohl die erw¨ahnten Schwierigkeiten in der vorgestel-

len 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 kontinuumsmechani-

sche und mathematische Kenntnisse notwendig. Dazu kommen noch die verschiedenen Vorge-

hensweisen wie die semianalytische, diskrete , variationelle usw., die auf unterschiedliche Va-

rianten der analytischen Sensitivit¨asanalyse f¨uhren. Im Rahmen dieser Arbeit wurde die varia-

tionelle 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¨atsaussa-

gen und somit eine ¨

Uberf¨uhrung der Variationen in Ableitungen. Zun¨achst wird die Variation

der schwachen Form des Gleichgewichts im i-ten Zeitschritt betrachtet

Gint

Gint

Gext

i

i

i

sGi =

sui +

si-1 +

s = 0

(4.2)

ui

i-1

s

K

^

i

Si

Fi

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

Gint =

Grad : C() : Gradu dV

(4.3)

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

C()

sGint =

Grad :

s : Gradu dV

(4.4)

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

~

sG int =

T

BT DB

I

I

J uJ dVe

s

(4.5)

e

I

J

e

Se

Die Berechnung der ~

D-Matrix wird exemplarisch f¨ur den ebenen Spannungszustand gezeigt

20

---

1

0

E()

D =

1

0

1 - 2

; und

mit

E() = 2E0

folgt

0

0

1-

2

1

0

~

D

2E

D =

=

0

1

0

1 - 2

0

0

1-

2

Es wird festgestellt, dass nur der Vorfaktor ge¨andert wird und die Struktur erhalten bleibt. ¨

Ahnli-

ches gilt auch f¨ur die Element-Sensitivit¨atsmatrix Se 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 sui bestimmt, aus der dann die Va-

riation der Verzerrungen si berechnet werden kann. Daraus lassen sich die Variation der Ver-

zerrungsenergiedichte si und die Variation des Dichteinkrements s() und schliesslich die

gesuchte si berechnen. (4.2) umgeformt liefert

sui = -K-1(S si-1 - sF)

(4.6)

Daraus wird die Variation der Verzerrung berechnet

si = B sui

(4.7)

Die Variation der Verzerrungsenergiedichte l¨asst sich wie folgt ausrechnen

1

1

si =

iD0sii-1 +

iD0isi-1

(4.8)

0

20

Aus (4.8) kann

1

s() =

tsi

(4.9)

ref

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

si = si-1 + s()

(4.10)

21

---

¨

Ahnlich wird auch bei der Variation des externen Anteils sGext 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

n

sGext =

F

(,s)

s

ii

(4.11)

i=1

(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

---

Schleife ¨uber alle Zeitschritte: i= 1 bis j

Strukturanalyse

Erf¨ullung der Gleichgewichtsbedingungen und L¨osung des LGS

Gi = Gint(, u

(,

i

i, i-1) - G ext

i

i-1, s) = 0

Kiui = Fi

Schleife ¨uber alle Elemente: k = 1 bis n

Berechnung der Verzerrungen

k = Buk

Berechnung der Verzerrungsenergiedichte

1

k =

kD0kk-1

20

Berechne das Inkrement der normierten Dichte

k

k =

- 1 t

ref

Update der normierten Dichte

k = k-1 + k

min,

k < min

k =

max,

k > max

Sensitivit¨atsanalyse

Schleife ¨uber alle Designvariablen

sui = -K-1(S

i

i si-1 - sFi)

innere-Schleife ¨uber alle Elemente k=1 bis n

sk = B suk

1

1

sk =

kD0skk-1 +

kD0ksk-1

0

20

1

s(k) =

tsk

ref

sk = sk-1 + s(k)

0,

k < min

sk =

0,

k > max

23

---

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 Speicher-

bedarf abgesch¨atzt werden. Innerhalb eines Zeitschrittes sind die folgenden Rechenoperationen

durchzuf¨uhren:

· Assemblierung der globalen Steifigkeitsmatrix Ki und L¨osen des reduzierten linearen

Gleichungssystems. Ki 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 voll-

besetzt, 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 sui = -K-1(S

aus der

i

i si-1 - sFi) zu l ¨

osen, wobei K-1

i

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 Re-

chenzeit 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 ge-

sagt 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 Mi-

nimierungsproblemen vorgestellt.

Gegeben sind eine Menge X

n

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 )

xX

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 vor-

gestellt.

5.1

Optimalit¨atsbedingungen

F¨ur die folgenden Betrachtungen wird vorausgesetzt, dass X

n

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) = x3 oder f (x) = -x2 zeigen, dass (5.1) nur eine notwendige, aber keine hin-

reichende 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 Ma-

xima zu unterscheiden, wird jetzt eine notwendige Optimalit¨atsbedingung zweiter Ordnung be-

trachtet.

Die Funktion f sei in einer Umgebung vo x X zweimal stetig differenzierbar. Ist x lokales

Minimum von (PU), dann gilt (5.1) und

xT

2f (x) x 0

x

n

R

(5.2)

d.h.

2f 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 lo-

kales Minimum ist. Das kann man wieder mit einfachen 1-D Beispielen wie f (x) = x3 mit

x = 0 zeigen. Hier gilt zwar (5.1) und (5.2), es liegt aber kein Minimum vor. Um eine hinrei-

chende Optimalit¨atsbedingung f¨ur ein lokales Minimum zu erhalten, muss man die Bedingung

(5.2) etwas versch¨arfen.

Seien X

n

R und f : X R zweimal stetig differenzierbar. Gelten

·

f (x) = 0 und

·

2f (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

n

R angelangt, so sucht man sich eine Richtung

d

n

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 :

n

n

n

R R und x R . Ein Vektor d R heisst Abstiegsrichtung von f in x, wenn es

ein t > 0 gibt mit

f (x + td) < f (x)

t 0, t

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

tung

Lemma 5.2

Seien f :

n

n

n

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 bestim-

men, siehe z.B. Geiger& Kanzow [3], Alt[1] und Polak[2]. In der erw¨ahnten Literatur wird ge-

zeigt, 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 Va-

rianten benutzt wird oder die strenge Wolfe-Powell-Regel, die sich beim CG-Verfahren ein-

geb¨urgert hat, werden nicht weiter betrachtet. Der interessierte Leser sei auf die erw¨ahnte Lite-

ratur verwiesen.

Sei f :

n

R

R stetig differenzierbar, und seien die Zahlen 0, 1 , [, 1) fest

2

vorgegeben. Zu x, d

n

R mit

f (x)T d < 0 bestimme man ein t > 0 mit

f (x + td) f (x) + t f (x)T d

(W P 1)

27

---

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 1 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 Abstiegs-

verfahren ist. Diese Vor¨uberlegungen sollen auch einen besseren ¨

Ubergang zum Quasi-Newton-

Verfahren erm¨oglichen.

Die Funktion f :

n

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

1

qk (x) := f xk +

f xk T x - xk +

x - xk T

2f xk

x - xk

(5.3)

2

zu minimieren versucht, wobei xk

n

R den aktuellen Iterationspunkt bezeichnet. Ist die Hesse-

Matrix

2f xk positiv definit, so ist xk+1 genau dann L¨osung von

min qk (x) ,

(5.4)

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

28

---

wenn xk+1 der Bedingung

qk (x) =

f xk +

2f xk

x - xk = 0

(5.5)

f¨ur einen station¨aren Punkt von qk gen¨ugt. Hieraus folgt f¨ur xk+1

xk+1 = xk -

2f xk -1 f xk

(5.6)

Um die explizite Berechnung der inversen Hesse-Matrix zu vermeiden, bestimmt man zun¨achst

eine Suchrichtung dk

n

R als L ¨

osung des linearen Gleichungssystems

2f xk dk = - f xk

(5.7)

und setzt anschliessend xk+1 := xk + dk. 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 globalisiertes2 Newton-Verfahren. Die Idee dieses globalisierten Ver-

fahrens besteht darin, die L¨osung dk 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

w¨ahlt dk = - f xk , 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 {xk}

n

R

eine Folge mit Grenzwert x. Die Folge heisst linear konvergent, wenn es

ein 0 < L < 1 gibt mit

xk+1 - x L xk - x

f ¨

ur

k 0

2In 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 {xk} heisst superlinear konvergent, wenn folgendes gilt

xk+1 - x

lim

= 0

k

xk - x

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

xk+1 - x c xk - x 2

f ¨

ur

k 0

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 soge-

nannten Quasi-Newton-Verfahren. Da die Berechnung der zweiten partiellen Ableitungen der

Hesse-Matrix

2f xk sehr aufwendig sein kann, versucht man anstelle der exakten Hesse-

Matrix eine hinreichend gute Approximation durch m¨oglichst einfach berechenbare Matrizen

Hk zu erreichen. Das ist auch der wesentliche Unterschied zum Newton-Verfahren. Diese Ap-

proximation 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

PSB3-, DFP4 und die BFGS5-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 xk+1

erh¨alt man die folgende Gleichung

f xk+1 =

f xk +

2f xk

xk+1 - xk + o

xk+1 - xk

(5.8)

Durch elementare Umformungen und gewisse Anforderungen6 an die Approximation Hk+1 der

Hessematrix erh¨alt man die untenstehende Gleichung

Hk+1 xk+1 - xk =

f xk+1 -

f xk

(5.9)

(5.9) wird oft auch als Quasi-Newton-Gleichung oder auch als Sekantengleichung bezeich-

net. Der folgende Algorithmus soll eine Zusammenfassung des globalisierten Quasi-Newton-

Verfahrens mit der direkten BFGS-Aufdatierungsformel und der Wolfe-Powell-Schrittweitenstrategie

3Powell-symmetric-Broyden

4Davidon-Fletcher-Powell

5Broyden-Fletcher-Goldfarb-Shanno

6siehe dazu Geiger & Kanzow S. 130

30

---

geben.

Algorithmus 3

A1:

W¨ahle x0

n

n×n

R , H0 R

symmetrisch und positiv definit,

0, 1 ,

(, 1) ,

0,

und

setze

k = 0

2

A2:

Ist

f xk

:

ST OP

A3:

Bestimme dk aus

dk := -H-1 f xk

k

A4:

Bestimme ein tk > 0, so dass die Wolfe-Powell-Bedingungen

f xk + tkdk f xk + tk f xk T dk

f xk + tkdk T dk f xk T dk

er¨ullt sind.

A5:

Setze

xk+1 := xk + tkdk,

sk := xk+1 - xk,

yk :=

f xk+1 -

f xk

und

yk yk T

Hksk sk T Hk

Hk+1 := Hk +

-

(yk)T sk

(sk)T Hksk

A6:

Setze

k k + 1, und gehe zu A2.

Die in A5 angegebene Formel f¨ur Hk+1 ist die schon erw¨ahnte BFGS-Updateformel der Ap-

proximation der Hesse-Matrix. Auch hier gilt zu erw¨ahnen, dass die globale Konvergenz durch

Gradientenschritte erzwungen werden kann, falls die Winkelbedingung7 nicht erf¨ullt ist. Al-

lerdings stellt sich dann die Frage, wie die neue Matrix Hk+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-

7Die Winkelbedingung besagt, dass der Winkel zwischen dk und - f xk gleichm¨assig von 90 weg be-

schr¨ankt bleibt

31

---

zeugte Folge Hk gegen

2f (x) konvergiert. Dazu muss noch sichergestellt werden, dass stets

die Schrittweite tk = 1 gew¨ahlt wird, sofern diese die Wolfe-Powell-Bedingungen gen¨ugt. Un-

ter 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 beob-

achtet man in der Praxis h¨aufig lokal eine schnelle Konvergenz.

In den folgenden Abschnitten dieses Kapitels werden noch einige Gradienten-Verfahren vor-

gestellt, die, wie sich im laufe der Zeit dieser Arbeit herausgestellt hat, geeigneter f¨ur unse-

re 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 Kom-

binationen 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 Appro-

ximation

1

f xk + d f xk +

f xk T d +

dT

2f xk d = mk (d)

(5.10)

2

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

min mk (d)

dRn

(5.11)

unter

der

N ebenbedingung

d k

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 Minimierungspro-

blem (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

f xk - f xk + dk

k =

(5.12)

f (xk) - mk (dk)

Falls k < 0 ist dk 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

xk+1 = xk + dk 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 allge-

meine 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 :

n

m

R

R eine zweimal stetig differenzierbare Abbildung und m n. Die all-

gemeine Form des Ausgleichsproblems ist

1

min f (x) =

F (x) 2

(5.13)

xRn

2

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

2f (x). F¨ur die Zielfunktion des Ausgleichsproblems (5.13) gilt

f (x) =

F (x)T F (x)

(5.14)

und

m

2f (x) =

F (x)T

F (x) + R (x)

mit

R (x) =

Fi (x)

2Fi (x)

(5.15)

i=1

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 Fi (x) klein sind. Durch diese Annahme wird die Vernachl¨assigung von

R (x) gerechtfertigt, die schliesslich auf die folgende Approximation f¨uhrt:

2f (x)

F (x)T

F (x)

(5.16)

33

---

Wendet man das ged¨ampfte Newton-Verfahren auf das Problem (5.13) an und ersetzt

2f (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

1

min

f xk T d +

dT

F xk T

F xk d

(5.17)

dRn

2

so erh¨alt man das Gauss-Newton-Verfahren. Mit der folgenden Berechnung wird gezeigt, dass

die Gauss-Newton-Richtung dk

n aus der L¨osung des resultierenden linearen Ausgleich-

GN

R

problems bestimmt werden kann.

Ist xk eine N¨aherungsl¨osung des nichtlinearen Ausgleichsproblems (5.13), so gilt f¨ur xk+1 =

xk + dk:

1

1

1

f xk+1 =

F xk+1

2 =

F xk + dk

2

F xk +

F xk dk 2

(5.18)

2

2

2

Es wird festgestellt, dass beim Gauss-Newton-Verfahren nur die ersten Ableitungen auftauchen

und dass ausgehend von einer Startn¨aherung x0 eine Folge linearer Ausgleichsprobleme zu

l¨osen ist. Schiesslich wird noch die Gauss-Newton-Richtung dk

angegeben, die die Anforde-

GN

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

-1

dk

= -

F xk T

F xk

F xk T F xk

(5.19)

GN

Auch das Gauss-Newton-Verfahren konvergiert in der Regel nur lokal, d.h. bei gen¨ugend gu-

ter 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 quadrati-

sche Approximation modifiziert werden.

34

---

1

mk (d) =

F xk +

F xk d 2

(5.20)

2

1

1

=

F xk

2 + dT T F xk F xk + dT T F xk

F xk d

2

2

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

1

min

F xk +

F xk d 2

dRn 2

(5.21)

unter

der

N ebenbedingung

d k

Man sieht, dass das gleiche Problem wie beim Gauss-Newton-Verfahren zu l¨osen ist, aller-

dings 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 Komle-

xit¨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¨osungskandi-

daten) 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 "Individu-

en" (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¨osungskan-

didatenmenge. 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 Ge-

neration.

Der Alorithmus wird ab Schritt zwei (Evaluation) wiederholt, oder nach einem Abbruchkriteri-

um beendet und der beste verf¨ugbare L¨osungskandidat als L¨osung definiert. Als Abbruchkrite-

rium k¨onnte man z.B. die maximale Anzahl der Generationen angegeben.

Im Vergleich zu den Gradienten-Verfahren k¨onnen hier keine aussagekr¨aftige Konvergenzeigen-

schaften 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 Verfah-

ren diskutiert, die eine sehr kleine Konvergenzgeschwindigkeit, daf¨ur aber eine grosse Konver-

genzsicherheit 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 Konver-

genzgeschwindigkeit gew¨ahrleistet. Das kann z.B. erreicht werden, indem zuerst mit dem Ge-

netischen Algorithmus m¨oglichst viele lokale Minima beseitigt werden, um so eine gute Aus-

gangslage 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 Aus-

schnitt 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 Be-

griffen nicht unterschieden.

37

---

Bild 6.1: Funktion mit lokalen Minima

· P1: Man geht davon aus, dass geeignete Startwerte vorhanden sind und versucht mit gra-

dientenbasierten 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 Such-

bereiche 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¨unsch-

te 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 Funkti-

on 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 prak-

tischen Problemen ist man zufrieden, wenn das Residuum um einen gewissen Betrag ver-

kleinert 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 Mo-

dellannahmen. 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

Bild 6.2: Effizienzerh¨ohung durch Kopplung der Verfahren

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 Test-

beispielen 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

Bild 6.3: Fachwerk mit unbekannter Belastung

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

N 2

=

(6.1)

203E0A2

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

N 2

ref =

(6.2)

203E0A2

Gleichung (6.2) wird nach N aufgel¨ost weden

Ni = ±

2ref 03E

i

0A2

(6.3)

und daraus werden die Normalkr¨afte der St¨abe bestimmt. Mit der Formulierung der Gleichge-

wichtsbedingungen am Lastangriffspunkt k¨onnen die einzelnen Lastkomponenten der Kraft F

berechnet werden.

FV = N1 · 0.6

(6.4)

und

FH = -N1 · 0.8 - N2

(6.5)

6.3

2D -und 3D-Beispiele

Im vorigen Abschnitt konnte festgestellt werden, dass das inverse Problem f¨ur eindimensiona-

le Problemstellungen keine grossen Schwierigkeiten bereitet. F¨ur zwei -und dreidimensionale

Systeme wird das nicht mehr der Fall sein. Es wird die Verzerrungsenergiedichte (2.8) betrach-

tet, die f¨ur den ebenen Spannungszustand wie folgt angegeben werden kann

E

1

=

0

2 + 2 +

(1 - µ)2 + 2µ1122

(6.6)

2

11

22

12

0(1 - µ2)

2

Aus (6.6) wird die Schwierigkeit der Problematik klar, da aus der skalaren Gr¨osse die vekto-

rielle 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 Lastangriffs-

punkte 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

Bild 6.4: System mit unbekannter Belastung

Sowohl die kinematischen Randbedingungen als auch der Lastangriffspunkt sind bekannt. Die

optimale Dichteverteilung kann ebenfalls der oberen Abbildung entnommen werden. Da in die-

sem Beispiel nur eine Last angreift und somit zwei Kraftkomponenten vorhanden sein m¨ussen,

kann der Fehler durch unterschiedliche Lastkombinationen berechnet und grafisch dargestellt

werden.

Bild 6.5: Fehler in Abh¨angigkeit der verschiedenen Lastkombinationen

42

---

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.

Bild 6.6: Fehler in Abh¨angigkeit der verschiedenen Lastkombinationen

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 dement-

sprechend (70000, -70000)T einer Druckkraft.

43

---

Bild 6.7: links: Zugkraft, rechts: Druckkraft

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 Dichtever-

lauf 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 Berech-

nung bei der als Startwert (-55000, 100000)T gew¨ahlt wurde, etwas ausf¨uhrlicher.

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

11

5

6

Funktionsauswertungen

57

30

70

Zeit [s]

562

295

699

Residuum

10-11

10-12

10-7

L¨osung F

(-70000, 70000)T

(-70000, 70000)T

(-70000, 70000)T

numerisch

Anzahl Iterationen

11

5

6

Funktionsauswertungen

57

30

70

Zeit [s]

773

406

962

Residuum

10-11

10-12

10-7

L¨osung F

(-70000, 70000)T

(-70000, 70000)T

(-70000, 70000)T

Tabelle 6.1: Vergleich der numerischen Resultate

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 Ab-

schnitt 5.6) . Obwohl auch die Levenberg-Marquardt-Methode f¨ur nichtlineare Ausgleichspro-

bleme 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¨osen-

den 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 Kon-

vergenz gegen die exakte L¨osung spielt die Diskretisierung eine wesentliche Rolle. Die Berech-

nungen 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 Proble-

matik soll der folgende Ausschnitt Abb. 6.6 veranschaulichen.

45

---

Bild 6.8: Ein wichtiger Bereich der Fehlerkurve

Als Startpunkt wurde P1 = (-80000, 15000)T gew¨ahlt. Mit dem Gauss-Newton-Verfahren

wurde nach 931s P2 = (-52192, 21895)T erreicht. Dabei wurde das Residuum von 825.3 in

P1 auf 303.4 in P2 verkleinert. In Abb. 6.8 sieht man sehr sch¨on, dass mit P2 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 P2 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 P3 =

(-78692, 58951)T mit einem Residuum von 135.6. P3 wurde wieder als Startwert f¨ur das

Gauss-Newton-Verfahren benutzt, welches dann das globale Minmum in P4 = (-70000, 70000)T

nach 400s finden konnte.

46

---

6.3.2

Beispiel 2

Das Vorgehen ist analog zum ersten Beispiel

Bild 6.9: System mit unbekannter Belastung

Auch in diesem Fall sind sowohl die kinematischen Randbedingungen als auch der Lastangriffs-

punkt bekannt. Da der Dichteverlauf nicht symmetrisch ist, sollten auch die Kraftkomponenten

unterschiedlich gross sein. Zun¨achst wird wieder die Fehlerkurve betrachtet, die durch unter-

schiedliche Kraftkombinationen ermittelt wurde.

Bild 6.10: Fehler in Abh¨angigkeit der verschiedenen Lastkombinationen

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.

Bild 6.11: links: Zubgbelastung, rechts: Druckbelastung

Die Kraftkomponenten sind wie erwartet unterschiedlich. Auch hier sollen die mit den Optimie-

rungsverfahren erzielten Ergebnisse diskutiert werden. Dazu wird eine Berechnung betrachtet,

bei der als Startpunkt (70000, -20000)T gew¨ahlt wurde.

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

11

6

7

Funktionsauswertungen

69

37

74

Zeit [s]

702

369

734

Residuum

10-8

10-12

10-9

L¨osung F

(70000, -10000)T

(70000, -10000)T

(70000, -10000)T

numerisch

Anzahl Iterationen

11

6

7

Funktionsauswertungen

69

37

74

Zeit [s]

956

508

1010

Residuum

10-8

10-12

10-9

L¨osung F

(70000, -10000)T

(70000, -10000)T

(70000, -10000)T

Tabelle 6.2: Vergleich der numerischen Resultate

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 Startwer-

te 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 konver-

gierten. Aufgrund der unterschiedlichen Suchrichtungen der Verfahren, kann es dazu kommen,

dass durch die relativ grobe Diskretisierung an unterschiedlichen Stellen lokale Minima gene-

riert 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 dis-

kutierten 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

Bild 6.12: System mit unbekannter Belastung

F¨ur dieses Beispiel wurde keine Fehlerkurve durch unterschiedliche Lastkombinationen ermit-

telt, 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 darge-

stellt. Alle drei Verfahren haben den exakten L¨osungsvektor (-100, -300, -300, -300)T ge-

liefert.

50

---

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

34

6

13

Funktionsauswertungen

200

49

122

Zeit [s]

989

243

604

Residuum

10-11

10-20

10-15

numerisch

Anzahl Iterationen

34

6

13

Funktionsauswertungen

200

49

122

Zeit [s]

1646

403

1005

Residuum

10-11

10-20

10-15

Tabelle 6.3: Vergleich der numerischen Resultate

Schon bei dieser Berechnung ist die Konvergenzrate der Levenberg-Marquardt-Methode deut-

lich 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

Bild 6.13: links: Zugbelastung, rechts: Druckbelastung

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

Bild 6.14: 3D-System mit unbekannter Belastung

wobei die Lasten zun¨achst als Platzhalter dargestellt wurden. Erst durch die Berechnungen wer-

den Gr¨osse und Richtung festgelegt. Der Abb. 6.14 kann entnommen werden, dass die Belas-

tung 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 abge-

baut wird (blauer Bereich). Die Berechnungen wurden mit folgenden Startwerten durchgef¨uhrt:

(50, 110, 60, -110, -60, -115)T .

52

---

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

34

8

9

Funktionsauswertungen

308

79

102

Zeit [s]

3737

959

1238

Residuum

10-12

10-20

10-16

numerisch

Anzahl Iterationen

34

8

9

Funktionsauswertungen

308

79

102

Zeit [s]

6999

2308

2352

Residuum

10-12

10-21

10-16

Tabelle 6.4: Vergleich der numerischen Resultate

Alle drei Verfahren lieferten die richtige L¨osung (100, 100, 50, -100, -50, -100)T . Der Un-

terschied 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 geeig-

neter 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 Proble-

men 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 sogen-

nanten "Standardized Femur", [28], diskutiert werden. Als Mass f¨ur die G¨ute der mittels Opti-

mierungsverfahren 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 Dich-

teverteilung dargestellt.

Bild 7.1: Dichteverteilung aus CT-Daten

Die Idee besteht darin, aus dieser Dichteverteilung auf die einwirkenden Lasten zur¨uckzusch-

liessen. Dabei werden Berechnungen mit zwei, drei -und f¨unf Kr¨aften durchgef¨urt. Weiter ste-

hen die ver¨offentlichten Messwerte zur Verf¨ugung, die beim normalen Gehen gemessen wurden,

[25], [26], [27].

54

---

Last [N]

x

y

z

Knotennummer

Messwerte f¨ur zwei Kr¨afte

F1

179.68

111.23

-221.39

1670

F2

-962.62

-466.34

1911.22

3565

Messwerte f¨ur drei Kr¨afte

F1

179.68

111.23

-221.39

1670

F2

-962.62

-466.34

1911.22

3565

F3

62.89

115.23

-114.51

2352

Messwerte f¨ur f¨unf Kr¨afte

F1

179.68

111.23

-221.39

1670

F2

-962.62

-466.34

1911.22

3565

F3

62.89

115.23

-114.51

2352

F4

-11.31

-15.10

90.02

785

F5

25.73

69.58

215.72

1633

Tabelle 7.1: Messwerte f¨ur zwei-, drei- und 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 Modellan-

nahmen 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 wur-

den, anschliessend werden auch Ergebnisse vorgestellt die durch die Kopplung des Genetischen

- Algorithmus und den Gradienten-Verfahren ermittelt wurden.

55

---

Bild 7.2: Femur mit f¨unf Kr¨aften bzw. 15 Kraftkomponenten

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 (F1, F2)T wurden die folgenden Ergebnisse

erzielt.

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

86

11

14

Funktionsauswertungen

92

122

182

Zeit [h]

28

40

22.3

Residuum

2327

2412

2328

Tabelle 7.2: Vergleich der numerischen Resultate: zwei Kr¨afte

56

---

Last [N]

x

y

z

Knotennummer

Quasi-Newton-Verfahren

F1

622.7

579.2

-751.2

1670

F2

-709.92

-438.9

933.9

3565

Gauss-Newton-Verfahren

F1

746.0

658.0

-376.8

1670

F2

-952.3

-575.4

807.1

3565

Levenberg-Marquardt-Verfahren

F1

227.3

664.5

-1066.1

1670

F2

-77.4

-461.7

230.1

3565

Tabelle 7.3: Berechnete Lasten, 6 Designvariablen

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 Mini-

ma 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 Gelenk-

kraft F2 gr¨osser als die Muskelkraft F1, 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.

Bild 7.3: Dichteverteilung: Quasi-Newton-Verfahren, 6 Designvariablen

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.

Bild 7.4: Dichteverteilung: Gauss-Newton-Verfahren, 6 Designvariablen

Auch die Ergebnisse der Gauss-Newton-Methode stimmen sehr gut mit den optimalen Dichte-

werten ¨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 vorange-

gangenen 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 redu-

zieren. 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 Kno-

chens 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 (F1, F2, F3)T wurden die folgenden Ergebnis-

se erzielt.

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

106

11

8

Funktionsauswertungen

164

147

149

Zeit [h]

42.8

25.2

24.7

Residuum

2269

2017

2324

Tabelle 7.4: Vergleich der numerischen Resultate: drei Kr¨afte

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 wich-

tige Erkenntnis aus der Tabelledie 7.5, dass die Gelenkkraft F2 nach wie vor eine Drucklast

ist. Daraus kann man folgern, dass die ermittelten Lasten nicht in Widerspruch zu den physi-

kalischen Gesetzen stehen. Diese Ergebnisse sind aus mathematischer Sicht auch ohne diese

¨

Uberlegungen korrekt und f¨uhren zu den drei gefundenen Minimalstellen.

59

---

Last [N]

x

y

z

Knotennummer

Quasi-Newton-Verfahren

F1

146.8

341.7

-331.0

1670

F2

-898.4

-352.3

1369.7

3565

F3

759.1

260.1

-204.6

2352

Gauss-Newton-Verfahren

F1

112.1

82.9

-70.9

1670

F2

-2044.5

-306.9

1590.3

3565

F3

1771.4

302.0

461.7

2352

Levenberg-Marquardt-Verfahren

F1

133.8

-39.3

237.4

1670

F2

-1471.2

-501.5

1074.1

3565

F3

809.6

581.6

-841.4

2352

Tabelle 7.5: Berechnete Lasten, 9 Designvariablen

Auch in diesem Fall sind die berechneten Kr¨afte unterschiedlich gross und haben unterschied-

liche Vorzeichen. Beim Gauss-Newton-Verfahren ist die Muskelkraft F3 relativ gross. Die Er-

gebnisse, die durch das Quasi-Newton- und das Levenberg-Marquardt-Verfahren erzielt wur-

den 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.

Bild 7.6: Dichteverteilung: Quasi-Newton-Verfahren, 9 Designvariablen

Da die Gesamtkraft auf drei Lasten verteilt wird und somit lokal kleinere Spannungen vorhan-

den sind, ist auch die Dichte in der Umgebung der Muskelkraft F1 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.

Bild 7.7: Dichteverteilung: Gauss-Newton-Verfahren, 9 Designvariablen

Beim Gauss-Newton-Verfahren ist das Residuum kleiner als bei den anderen zwei Verfahren

und die Dichteverteilung ist im Schaftbereich ebenfalls besser. Im proximalen1 Bereich sieht

aber die Dichteverteilung weniger gut aus. Grund daf¨ur k¨onnte die sehr grosse Muskelkraft

F3 sein. Auch die Gelenkkraft F2 ist relativ gross und kann als Ursache f¨ur die lokal hohen

Dichtewerte betrachtet werden.

1Oberst¨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 (F1, F2, F3, F4, F5)T

genommen und die folgenden Resultate erzielt

Sensitivit¨atsanalyse

Optimierungsverfahren

Quasi-Newton

Gauss-Newton

Levenberg-Marquardt

analytisch

Anzahl Iterationen

108

9

13

Funktionsauswertungen

183

211

311

Zeit [h]

65.2

52.3

72.3

Residuum

2278

2286

1949

Tabelle 7.6: Vergleich der numerischen Resultate: f¨unf Kr¨afte

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.

Last [N]

x

y

z

Knotennummer

Quasi-Newton-Verfahren

F1

216.6

304.7

-174.6

1670

F2

-932.1

-415.4

1447.8

3565

F3

390.7

79.5

-132.6

2352

F4

-7.2

-62.9

40.8

785

F5

183.9

348.9

71.8

1633

Gauss-Newton-Verfahren

F1

12.0

-155.5

280.0

1670

F2

-1340.0

-378.3

1649.3

3565

F3

325.2

-691.1

333.5

2352

F4

748.8

-229.8

-106.6

785

F5

198.9

1432.0

-1851.5

1633

Levenberg-Marquardt-Verfahren

F1

68.6

32.9

101.3

1670

F2

-743.1

-252.8

992.8

3565

F3

364.0

-422.4

121.3

2352

F4

-195.1

61.7

638.7

785

F5

298.4

781.1

-1544.7

1633

Tabelle 7.7: Berechnete Lasten, 15 Designvariablen

Es ist weiterhin so, dass die Gelenkkraft F2, 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 F1 und F5 aus. Mit den letzten zwei Verfahren wurde eine zu kleine Muskelkraft

F1 und eine zu grosse Muskelkraft F5 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. Insbe-

sondere die R¨ohrenform des Femurs ist weiterhin deutlich zu sehen. Aufgrund der relativ hohen

Gelenkkraft F2 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 F1 nur sehr schwach vorhanden, daf¨ur sind sie in der Umgebung

64

---

der recht grossen Gelenkkraft F2 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 F2 sind nicht

mehr zu beobachten. Daf¨ur wird aber in der Umgebung der Muskelkraft F1 recht stark Dichte

angebaut. Grund daf¨ur k¨onnte die relativ hohe, negative z-Komponente der Muskelkraft F5 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 ei-

nem 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:

Last [N]

x

y

z

Knotennummer

Genetischer-Algorithmus

F1

-100.3

-56.4

-153.5

1670

F2

-604.4

-380.9

36.5

3565

F3

599.7

7.2

-242.3

2352

F4

96.4

-139.5

-158.2

785

F5

-158.2

-94.36

525.0

1633

Tabelle 7.8: Berechnete Lasten, 15 Designvariablen

Dieser Lastvektor wurde wieder als Startwert f¨ur den Gauss-Newton-Verfahren benutzt, das das

Residuum nach 41h und 174 Funktionsauswertungen auf 1700 reduzierte

Last [N]

x

y

z

Knotennummer

Gauss-Newton-Verfahren

F1

-100.3

-67.2

-159.3

1670

F2

-545.4

-379.1.9

1.0

3565

F3

574.9

-22.9

-326.1

2352

F4

176.2

-218.7

-54.7

785

F5

-50.6

559.6

70.35

1633

Tabelle 7.9: Berechnete Lasten, 15 Designvariablen

Es wird festgstellt, dass das Residuum nicht mehr wie am Anfang der Berechnugen stark redu-

ziert wird. Da die Gelenkkraft F2, 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 Dichte-

verteilung 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 F1 als

auch in der Umgebung der Gelenkkraft F2 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 folgen-

den 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, ¨uberle-

gen 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 Ef-

fizienz 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 Un-

terschied zwischen der numerischen und der analytischen Sensitivit¨atsanalyse deutlich zu ma-

chen, 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 Sen-

sitivit¨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 Konvergezgeschwindig-

keit sehr klein ist. Weiter konnte festgestellt werden, dass mit der Kopplung des Genetischen-

Algorithmus und den Gradienten-Verfahren das Residuum um einen wesentlichen Betrag redu-

ziert 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

Bild 8.1: links: optimale Dichteverteilung, die restlichen drei: berechnet

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 ¨

Ubereinstim-

mung, 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 Zielfunkti-

on. 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 Referenzverzerrungsener-

giedichte ref als Designvariable definiert und in den Optimierungsprozess integriert, k¨onnte

einen wesentlichen Beitrag dazu leisten. Die Ergebnisse der Berechnungen im Kapitel 7 ha-

ben 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 Berech-

nungen f¨uhren. Eine weitere Verbesserungsm¨oglichkeit liegt darin, die wesentlichen Lasten wie

die Gelenkkraft F2 und die Muskelkraft F1 mit einer h¨oheren Gewichtung in die Berechnung

einzubeziehen. Man konnte feststellen, dass auch mit zwei bzw. drei Kr¨aften sehr gute Ergeb-

nisse erzielt werden konnten. Aus diesem Grund ist es naheliegend, die Berechnungen auf die

dominierenden Lasten F1, F2 und F3 (siehe Abb. 7.2) zu beschr¨anken. Um lokal hohe Bean-

spruchnungen zu vermeiden, k¨onnten die Kr¨afte als Fl¨achenlasten definiert werden.

71

---

Zusammengefasst kann gesagt werden, dass mehr Aufwand in den oben angesprochenen Ver-

besserungsm¨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 Opti-

mierungsaufgaben, 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 Op-

timierungsaufgaben, 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 Zim-

mermann, 2. Auflage Springer 2002

[9] D. BRAESS: Finite-Elemente, Theorie, schnelle L¨oser und Anwendungen in der Elasti-

zit¨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, Dis-

sertation, 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 Be-

reich der Mechanik der Universit¨at Hannover, 1998

[15] K. WIECHMANN: Theorie und Numerik zur Berechnung und Optimierung von Struk-

turen 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 Sys-

tems, 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, Techni-

sche 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 enginee-

ring and physics 26, 2004

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

75

---