10 Parameteridentifikation

Es wird zunächst ein allgemeines Verfahren vorgestellt, mit dem sich die Parameter des $PT_2$-Systems ($T_1$, $T_2$, $k_{P12}$), sowie des $PDT_1$-Systems ($T_1$, $k_{P1}$, $T_D$) bestimmen lassen. Danach werden die Ergebnisse der Parameteridentifikation vorgestellt.

10.1 Die allgemeine Fehler-Quadrat-Methode

Gegeben seien die $N$ Messpunkte

$$t_1, y_1; t_2, y_2; \dots; t_N, y_N \,.$$ (29)

Die Parameter $a_1,\dots,a_M$ einer gegebenen Kurvenschar

$$y=f(t;a_1,\dots,a_M)$$ (30)

sollen diesen Daten möglichst gut angepasst werden. Dazu wird folgendes Minimumforderung mit dem Fehlerfunktional $E(f)$ verwendet [17]:

$$E(f)=\sum^{N}_{i=1} (y_i-f(t_i;a_1,\dots,a_M))^2=min!.$$ (31)

10.2 Die Parameteridentifikation des aktiven Verhaltens

10.2.1 Die Eingangs- und Ausgangssignale

Zur Bestimmung der Systemparameter werden zunächst aus den Messdaten der Proben 6-C, 7-C und 9-F das normierte Eingangssignal $x(t)$ und das Ausgangssignal $y(t)$ aus den Messwerten gebildet. Die Normierung wurde für das Eingangssignal so gewählt, dass der größte gemessene Wert mit eins zusammenfällt. Die Normierung des Ausgangssignals wurde mit den Normierungskonstanten $F_0$ aus Tabelle 2 durchgeführt. Da die zeitkontinuierlichen Messgrößen diskret abgetastet werden, liegen $N$ Messpunkte jeweils des Ein- und Ausgangssignals in der Form vor:

$$t_1, x_1; t_2, x_2; \dots; t_N, x_N \quad \hbox{bzw.} \quad t_1, y_1; t_2, y_2; \dots; t_N, y_N$$ (32)

10.2.2 Modellansatz

Um das aktive Verhalten zu modellieren, wählen wir einen Ansatz gemäß Gleichung 21 mit den Parametern $T_1$, $T_2$ und $k_{P12}$. Da die noch zu ermittelnden Parameterschätzwerte für das aktive Verhalten mit denen für das passive Verhalten im späteren Verlauf verglichen werden sollen, ist es sinnvoll, neue Bezeichner einzuführen, um Verwechslungen zu vermeiden. Um deutlich herauszustellen, dass es sich bei bestimmten Größen um Parameterschätzwerte handelt, werden diese zusätzlich mit einem Dach gekennzeichnet. Für die Parameterschätzwerte des aktiven Verhaltens werden wir im folgenden schreiben: $T_1 \to \hat{T}^{akt}_1$, $T_2 \to \hat{T}^{akt}_2$ und $k_{P12} \to \hat{k}^{akt}_P$. Analog zu Gleichung 30 schreiben wir hier:

$$y=f(t;T^{akt}_1,T^{akt}_2,k^{akt}_P)=y^\ast(t) \quad \mbox{mit}\quad T^{akt}_1,T^{akt}_2,k^{akt}_P=const.$$ (33)

Dabei gilt nach Gleichung 21:

$$T^{akt}_1\,T^{akt}_2\,\ddot y^\ast(t)+(T^{akt}_1+T^{akt}_2)\,\dot y^\ast(t)+y^\ast(t)=k^{akt}_P\,x(t)\,.$$ (34)

Um $y^\ast(t)$ zu bestimmen, wird die inhomogene Differenzialgleichung 34 numerisch mit Hilfe des Programmes Mathematica 3.0 gelöst. Mathematica greift dabei auf die Funktion ''LSODE'' aus der ''ODEPACK''-Bibliothek von Alan Hindmarsh [12] zurück, die auf dem BDF-Verfahren von Gear basiert.

Nun kann das parameterabhängige Fehlerfunktional $E(f)=E(T^{akt}_1,T^{akt}_2,k^{akt}_P)$ mit der numerischen Lösung $y^\ast(t)$ nach Gleichung 30 bestimmt werden. Das Minimum von $E$ wird dadurch angenähert, dass ausgehend von einem Starttripel $(T^a_1,T^a_2,k^a_P)$ für alle Kombinationen der Parameterwerte $T^a_1 \pm \Delta_{T1}$, $T^a_2 \pm \Delta_{T2}$ und $k^a_P \pm \Delta_{kp}$ der Fehler $E$ bestimmt wird. Dabei sind die $\Delta_{T1}$, $\Delta_{T2}$ und $\Delta_{kp}$ geeignet zu wählen. Das Tripel mit dem kleinsten Fehler wird als neuer Startwert gewählt. Für den Fall, dass das Starttripel den kleinsten Fehler besitzt, werden die $\Delta_{T1}$, $\Delta_{T2}$ und $\Delta_{kp}$ verkleinert. Dieses Verfahren wird iterativ solange wiederholt, bis die änderung der Parameter $<1\%$ ist, oder der Fehler um weniger als $1\%$ pro Iteration abnimmt.

Etablierte Verfahren zur Parameterbestimmung von linearen Differentialgleichungen, wie sie zum Beispiel in [13] und [14] angegeben sind, wurden erprobt, lieferten aber augenscheinlich im Vergleich zum oben erklärten Verfahren keine gute Fehlerminimierung.

Abbildung 50 zeigt beispielhaft einen Vergleich zwischen gemessener und angepasster Systemantwort mit absolutem Fehler (Probe 9-F). Insgesamt wurde mehrmals für verschiedene Längen von drei verschiedenen Proben (6-C, 7-C und 9-F) eine Parameterbestimmung des aktiven Verhaltens durchgeführt. Die Normierungskonstanten der einzelnen Proben sind in Tabelle 2 angegeben. Abbildung 51 zeigt die bestimmten Parameter $\hat{T}^{akt}_1$ für verschiedene Längen und Proben, Abbildung 52 für $\hat{T}^{akt}_2$ und Abbildung 53 für $\hat{k}^{akt}_P$.

Probe 9-F zeigte schon bei der Messung ein anormales Kontraktionsverhalten. Sie wird in Abbildung 51 und 53 als Messausreißer behandelt. In Tabelle 6 sind die Mittelwerte über verschiedene Längen und Proben der ermittelten Parameter angegeben. Eine Mittelwertbildung der Systemparameter über die verschiedenen Längen scheint bei Probe 6-C und 7-C sinnvoll, da keine deutliche Längenabhängigkeit der Parameter zu erkennen ist.

Abbildung 50: Kontraktionsverlauf-Vergleich zwischen gemessenem und berechnetem Ausgangssignal (Probe 9-F) mit absolutem Fehler. Normierungskonstante $F_0$: siehe Tabelle 2

Abbildung: Links: Parameter $\hat{T}^{akt}_1$ über der Länge aufgetragen für drei verschiedene Proben, Mitte: Auschnittsvergrößerung mit Box-Whisker-Plot, Rechts: Histogramm; Normierungskonstante $l_0$: siehe Tabelle 2

Abbildung: Links: Parameter $\hat{T}^{akt}_2$ über der Länge aufgetragen für drei verschiedene Proben, Rechts: Histogramm; Normierungskonstante $l_0$: siehe Tabelle 2

Abbildung: Parameter $\hat{k}^{akt}_P$ über der Länge aufgetragen für drei verschiedene Proben, Rechts: Histogramm; Normierungskonstante $l_0$: siehe Tabelle 2

Tabelle 6: Parameterschätzwerte des aktiven Verhaltens für Gleichung 34

Parameter	Mittelwert	Standardabweichung
$\hat{T}^{akt}_1$	$0.231\,s$	$0.15\,s$
$\hat{T}^{akt}_2$	$4.49\,s$	$0.85\,s$
$\hat{k}^{akt}_P$	$2.50$	$0.90$

10.3 Die Parameteridentifikation des passiven Verhaltens

Tabelle 7: Parameterschätzwerte des passiven Verhaltens für Gleichung 36

Parameter	Mittelwert	Standardabweichung
$\hat{T}^{pas}_1$	$4.94\,s$	$1.45\,s$
$\hat{k}^{pas}_P$	$0.130$	$0.068$
$\hat{T}^{pas}_D$	$50.8\,s$	$31.\,s$

Es sollen nun die Parameter $T_1$, $k_{P1}$ und $T_D$ der Gleichung 26 an das passive Verhalten angepasst werden. Auch hier werden wie in Abschnitt 10.2.2 neue Bezeichner für die Parameter eingeführt, und insbesondere die ermittelten Parameterschätzwerte auch mit einem Dach versehen: $T_1 \to \hat{T}^{pas}_1$, $k_{P1} \to \hat{k}^{pas}_P$ und $T_D \to \hat{T}^{pas}_D$. Gleichung 30 schreiben wir hier analog zu Gleichung 33:

$$y=f(t;T^{pas}_1,k^{pas}_P,T^{pas}_D) =y^\ast(t) \quad \mbox{mit}\quad T^{pas}_1,k^{pas}_P,T^{pas}_D=const.$$ (35)

Dabei gilt nach Gleichung 26

$$T^{pas}_1\,\dot y^\ast(t)+y^\ast(t)=k^{pas}_P\,[x(t)+T^{pas}_D\, \dot x(t)]\,.$$ (36)

Nun können mit Hilfe des in Abschnitt 10.2.2 vorgestellten Verfahrens die Parameter $T^{pas}_1$, $k^{pas}_P$ und $T^{pas}_D$ der Gleichung 36 an das gemessene passive Verhalten angepasst werden. Abbildung 54 zeigt beispielhaft einen Vergleich zwischen gemessener Systemantwort auf eine sprunghafte Längenänderung und der Systemantwort des angepassten Systems. Die Parameterschätzwerte $\hat{T}^{pas}_1$, $\hat{k}^{pas}_P$ und $\hat{T}^{pas}_D$ wurden für drei verschiedene Proben bei unterschiedlichen Längen jeweils für sprunghafte Verkürzung und Längung der Probe bestimmt. Abbildung 55 zeigt beispielhaft die Längenabhängigkeit des Parameterschätzwertes $\hat{T}^{pas}_1$ für Probe 11-F bei sprunghafter Längung. Auch bei den anderen Parameterschätzwerten und Proben — sowohl für Längung als auch Verkürzung — ist die Abhängigkeit der Parameterschätzwerte von der Probenlänge zu vernachlässigen. Abbildung 56 zeigt die ermittelten Mittelwerte und Standardabweichungen für die verschiedenen Kategorien. Die Gesamtmittelwerte sind in Tabelle 7 angegeben.

10.4 Ergebnisse

Der Versuch, das aktive Verhalten mit einem $PT_2$-System und das passive Verhalten mit einem $PDT_1$-System zu identifizieren, liefert brauchbare Parameterschätzwerte (Tabellen 6 und 7). Das Reduzieren des mechanischen Systemverhaltens des Muskelgewebes auf diese einfachen linearen Systeme stellt natürlich eine Vereinfachung dar. Trotzdem sind diese einfachen Modelle in der Lage, die wesentlichen, gemessenen Systemeigenschaften der Muskulatur angenähert wiederzugeben.

Beim Vergleichen der ermittelten Systemparameter für das aktive und passive Verhalten fällt auf, dass die Zeitkonstanten, die für das Abklingverhalten verantwortlich sind (aktiv: $\hat{T}^{akt}_2=4.49\pm 0.85\,s$, passiv: $\hat{T}^{pas}_1=4.94\pm 1.45\,s$), bei beiden Systemen nahezu übereinstimmen. Dies deutet darauf hin, dass sowohl beim aktiven als auch beim passiven Verhalten teilweise die gleichen Prozesse ablaufen. Ein Modell, das diese Erkenntnis nutzt, wird im folgenden Abschnitt 11 vorgestellt.

Beim Betrachten der Box-Whisker-Plots in Abbildung 56 fällt auf, dass sich die Standardabweichung der einzelnen Proben kaum von der Gesamtstandardabweichung unterscheidet. Hier deutet sich an, dass selbst bei einer größeren Anzahl von Messungen die ermittelten Parameterschätzwerte nicht genauer bestimmt werden können. Dies ist eine für biologische Gewebe typische Eigenschaft und muß nicht von vornherein als Hinweis auf eine unpassende Modellierung verstanden werden.

Abbildung 54: Vergleich von gemessener und theoretischer Systemantwort (Probe 11-F). Normierungskonstante $F_0=3.1\,mN$

Abbildung: Parameter $\hat{T}^{pas}_1$ über der Länge aufgetragen (Probe 11-F, Längung der Probe). Normierungskonstante $l_0=4.9\,mm$

Abbildung: Box-Whisker-Plots der Parameter $\hat{T}^{pas}_1$, $\hat{k}^{pas}_P$ und $\hat{T}^{pas}_D$ für verschiedene Kategorien, oben: Histogramme der gesamten Daten