5. Problemy symulacji ruchu
5.1. ZAGADNIENIA SYMULACJI DYNAMICZNEJ PROSTEJ I ODWROTNEJ
Pod pojęciem symulacja dynamiczna rozumie się symulację (numeryczną) prowadzoną
z wykorzystaniem modelu dynamiki układu (dynamicznych równań ruchu). Wyró\nia się dwa
podstawowe rodzaje tej symulacji symulacjÄ™ dynamicznÄ… prostÄ… oraz symulacjÄ™ dynamicznÄ…
odwrotną, określane równie\ odpowiednio zadaniem prostym dynamiki (analizy/symulacji
ruchu) i zadaniem odwrotnym dynamiki (wyznaczania sterowania) [16,57,81]. W pierwszym
przypadku znane są siły lub zaburzenia działające na układ oraz stan układu w chwili począt-
kowej, a zadanie polega na całkowaniu równań ruchu w określonym przedziale czasu celem
wyznaczenia przebiegu ruchu układu pod działaniem tych sił. W drugim przypadku znany jest
ruch układu (zmierzony lub zadany w czasie), a wyznaczane są siły wywołujące ten ruch. W
zagadnieniach praktycznych część sił działających na układ jest zwykle zdeterminowana (za-
le\na w ogólności od czasu i stanu ruchu układu), a rozwiązaniem symulacji odwrotnej są
wymagane przebiegi w czasie jedynie (dodatkowych) sił sterujących. Tak rozumiane zadanie
odwrotne dynamiki (symulacji dynamicznej odwrotnej) znajduje od lat szerokie zastosowanie
przy syntezie sterowania robotami [35,81,106,129], obiektami latającymi [16,132] czy ukła-
dami biomechanicznymi [71,84,98,102,124,141,147].
Warto zauwa\yć, \e w literaturze polskiej zdefiniowane wy\ej zadania symulacji dy-
namicznej prostej i odwrotnej rozpoznawane są często inaczej, odpowiednio jako drugie (od-
wrotne) oraz pierwsze (proste, bezpośrednie) problemy/zadania dynamiki [9,85,101,126].
Istnieje zatem zasadnicza rozbie\ność co do tego, które z zadań dynamiki określać prostym, a
które odwrotnym. Aktualnie, w większości zastosowań technicznych (robotyce, biomechani-
ce, dynamice układów wieloczłonowych oraz biomechanice), przewa\ać zaczyna nazewnic-
two odnoszące się do angielskich odpowiedników [57,74,99,121,132,147]: forward dyna-
mics problem ( direct dynamic simulation ), tłumaczone jako zagadnie proste dynamiki
(symulacja dynamiczna prosta, w znaczeniu: do przodu ) oraz inverse dynamics problem
( inverse dynamic simulation ), tłumaczone jako zagadnienie odwrotne dynamiki (symulacja
dynamiczna odwrotna).
1
Symulacja dynamiczna odwrotna jest jednym z podstawowych narzędzi analizy ukła-
dów biomechanicznych [39,49,5071,84,98,102,124,133,147]. Podstawą dla tej analizy są ade-
kwatne dla rozwa\anego przypadku dynamiczne równania ruchu, przykładem których mogą
być równania (4.20), oraz zadane w czasie charakterystyki ruchu qd (t) pochodzące najczę-
ściej z pomiarów kinematycznych.
5.2. POMIARY KINEMATYCZNE I OBRÓBKA DANYCH DLA ZADAC SYMULACJI
DYNAMICZNEJ ODWROTNEJ
Niezbędne dla rozwiązania zadania odwrotnego przebiegi qd (t) uzyskuje się zwykle
na podstawie rejestracji doświadczalnej analizowanych czynności motorycznych, a podsta-
wowym sposobem tej rejestracji jest technika filmowa. Wykorzystuje się w tym celu układy
zsynchronizowanych ze sobą kamer o odpowiednio wysokiej częstotliwości (dla rejestracji
prostych czynności wykonywanych w płaszczyznie strzałkowej, typu wymachu kończyny
górnej lub dolnej, wystarcza często jedna kamera ustawiona prostopadle do płaszczyzny ru-
chu). Film jest zródłem sekwencji zdjęć, z których, po przekształceniu na format binarny,
odczytuje się poło\enia w czasie markerów umiejscowionych w stosownych punktach antro-
pometrycznych. Na rys 5.1 pokazano przykładowe dwie klatki takiego filmu z rejestracji sko-
ku na batucie, pochodzÄ…ce z prac [22,23,37], gdzie uwidocznione sÄ… markery punktowe i li-
niowe, a skok wykonywany jest na tle siatki umo\liwiającej kalibrację odczytów. Techniki te
oraz sposoby zapewnienia wymaganej dokładności pomiarów omawiane są szeroko w litera-
turze fachowej, m.in. [98,112,123,124,141].
Rys. 5.1. Markery na ciele sportowca wykonujÄ…cego skok na batucie
Pomiary kinematyczne dostarczają dyskretnych przebiegów zmian poło\eń markerów
w interwałach czasowych "t wynikających z częstotliwości kamer rejestrujących, przelicza-
nych następnie na (dyskretne) charakterystyki q* (t) . Celem ich przekształcenia na charakte-
d
rystyki ciągłe wymagane w zadaniu symulacji dynamicznej odwrotnej, q* (t) qd (t) , stosuje
d
się zwykle aproksymację/interpolację tych przebiegów odpowiednimi funkcjami analitycz-
nymi. JednÄ… z takich technik, zastosowanÄ… m.in. w pracach [22,23,25,26,37], jest aproksyma-
2
cja z wykorzystaniem wielomianowych funkcji sklejanych trzeciego stopnia (splajnów) [80,
128]. Stosując algorytm opisany w pracy [111] mo\liwe jest te\ sterowanie efektem wygła-
dzania przebiegów qd (t) za pomocą dobieranych przez u\ytkownika, zwykle metodą prób i
błędów, współczynników wagowych w punktach węzłowych rozwiązywane jest zadanie
aproksymacji (wartości wagowe większe od zera), zadanie interpolacji (wartości wagowe
równe zeru) lub zadanie mieszane.
. .
6
.
. .
4 .
. . .
2
.
.
0
0 2 4 6 8 10
\
współrzędna niezale na x
4
6
2
0
0
-2
-6
-4
0 2 4 6 8 10 0 2 4 6 8 10
\ \
współrzędna niezale na x współrzędna niezale na x
2 2 2
Rys. 5.2. Przebiegi y(x) , y (x) oraz y (x) dla jednokrotnego i dwukrotnego wygładzania
Za wykorzystaniem splajnów dla otrzymywania ciągłych (wyra\onych analitycznie)
& &&
charakterystyk qd (t) , a na tej podstawie równie\ qd (t) i qd (t) , przemawia dostępność opro-
gramowania standardowego. Wadą rozwiązania jest jednak jedynie ciągłość otrzymywanych
&&
tą drogą drugich pochodnych qd (t) (przebiegi w postaci łamanej), co odbija się na gładkości
rozwiązań zadania symulacji dynamicznej odwrotnej. Prostym i skutecznym sposobem mini-
malizacji tego efektu okazuje siÄ™ dwukrotne zastosowanie opisanej procedury aproksymacji
za pomocÄ… wielomianowych funkcji sklejanych trzeciego stopnia, wykorzystujÄ…c w drugim
etapie du\ą liczbę punktów węzłowych wziętych z funkcji aproksymującej qd (t) otrzymanej
w pierwszym etapie. Ten prosty zabieg pozwala istotnie (numerycznie) wygładzić przebiegi
&&
qd (t) finalnej funkcji aproksymującej qd (t) . Dla zilustrowania tego problemu posłu\my się
następującym prostym przykładem obliczeniowym. Dla N =11 punktów węzłowych pokaza-
nych na rys. 5.2 jako zaczernione kwadraty, przebieg funkcji y(x) otrzymanej z zastosowa-
3
aproksymowane
y
(
x
)
aproksymowane
y'
(
x
)
aproksymowane
y"
(
x
)
niem opisanej procedury w pierwszym etapie (z u\yciem N =11 punktów wyjściowych) po-
kazany jest jako linia cienka. Linia gruba pokazuje natomiast tÄ™ funkcjÄ™ otrzymanÄ… w drugim
2
etapie z u\yciem N =101 punktów węzłowych wziętych z rozwiązania z pierwszego etapu.
2
Przebiegi te niewiele się ró\nią. Jakościowo podobne są równie\ przebiegi y (x) . Istotną ja-
2 2
kościową ró\nicę widać natomiast dla przebiegów y (x) przebiegi te dla rozwiązań z dru-
giego etapu mają charakter zdecydowanie gładki w porównaniu z łamanym przebiegiem
tych funkcji z pierwszego etapu aproksymacji.
Ciągłe (zapisane analitycznie za pomocą wielomianowych funkcji sklejanych trzeciego
stopnia) przebiegi qd (t) , aproksymujÄ…ce charakterystyki dyskretne q* (t) otrzymywane do-
d
świadczalnie, w swej istocie odpowiadają wielu innym stosowanym technikom wygładzania
danych pomiarowych oraz filtrowania przypadkowych błędów pomiaru i odczytu przy otrzy-
mywaniu q* (t) , omawianym m.in. w [98,112,141]. Analityczny zapis qd (t) pozwala te\ na
d
& &&
bezpośrednie liczenie qd (t) i qd (t) , unikając w ten sposób szacowania tych pochodnych za
pomocą ilorazów ró\nicowych wykorzystujących (wygładzone) dane pomiarowe.
\
Rys. 5.3. Klatki z filmu rejestrującego ruch kończyny górnej obcią onej dodatkowym
\
odwa nikiem trzymanym w dłoni
Jako przykład charakterystyk kinematycznych otrzymanych opisaną metodą (podwój-
nego) aproksymowania przebiegów dyskretnych q* (t) otrzymanych z pomiarów za pomocą
d
splajnów, poka\my przypadek związany z rejestracją ruchu kończyny górnej obcią\onej do-
datkowym odwa\nikiem trzymanym w dłoni. Ze względu na prostotę i płaski charakter ruchu,
dla jego rejestracji za pomocą pojedynczej kamery wystarczyły trzy widoczne na rys. 5.3 mar-
kery (na zdjęciach tych widać równie\ elektrody dla pomiarów EMG, o czym będzie mowa w
dalszej części pracy). Model dynamiczny kończy górnej, analogiczny do opisanego w roz-
działach 4.2.1 i 4.3.1, wymaga rSd (t) = [ xSd (t) ySd (t) ]T (ruch stawu barkowego S) oraz
qd (t) = [ Õ1d (t) Õ2d (t) ]T (wychylenia ramienia i przedramienia od pionu), analogicznie jak
na rys. 4.9. Przebiegi rSd (t) otrzymywane są poprzez aproksymowanie przebiegów dyskret-
*
nych rSd (t) (pomiar poło\eń markera w stawie S). Charakterystyki dyskretne q* (t) wylicza-
d
ne są na podstawie mierzonych poło\eń markerów, a następnie aproksymowane splajnami dla
4
otrzymania qd (t) . Finalny wynik tych zabiegów, z zastosowaniem opisanej podwójnej
aproksymacji, pokazany jest na rys. 5.4. Dane te u\yte były dla obliczeń symulacji dynamicz-
nej odwrotnej, które będą omawiane w następnym rozdziale.
0.4 0.04 -0.04
.
.. . xS
xS [m/s2] xS [m/s] xS [m]
0.2 0 -0.06 ..
xS
0 -0.04 -0.08
xS
t [s]
-0.2 -0.08 -0.1
0 0.4 0.8 1.2 1.6
0.03 0.34
0.04
.
..
yS
yS [m/s2] yS [m/s]
yS [m]
0.015 0.33
0
.
..
yS yS
-0.04 0 0.32
t [s]
-0.08 -0.015 0.31
0 0.4 0.8 1.2 1.6
800 200 60
..
.. .
Õ1
Õ1 [deg/s2] Õ1 [deg/s] Õ1 [deg]
0 100 30
Õ1
-800 0 0
.
Õ1
t [s]
-1600 -100 -30
0 0.4 0.8 1.2 1.6
2000 600 180
.. . ..
Õ2
Õ2 [deg/s2] Õ2 [deg/s] Õ2 [deg]
Õ2
0 300 90
-2000 0 0
.
Õ2
t [s]
-4000 -300 -90
0 0.4 0.8 1.2 1.6
Rys. 5.4. Charakterystyki ruchu wymagane dla symulacji dynamicznej odwrotnej ruchu
kończyny górnej
5.3. WYZNACZENIE WYPADKOWYCH MOMENTÓW SIA MIŚNIOWYCH
W STAWACH
Symulacja dynamiczna odwrotna układów biomechanicznych ograniczona do wyzna-
czania wypadkowych momentów sił mięśniowych w stawach nale\y do zadań relatywnie pro-
stych i jednoznacznych. Dotyczy to w szczególności modeli wyizolowanych kończyn górnych
5
i dolnych, budowanych jako układy o strukturze otwartej zaczepione w przegubie (staw bar-
kowy/biodrowy) poruszajÄ…cym siÄ™ znanym ruchem. Liczba l = lÄ wypadkowych momentów
sił mięśniowych w stawach (przegubach modelu) jest wówczas równa liczbie k stopni swobo-
dy układu, l = k , a dynamiczne równania ruchu układu wyra\one we współrzędnych nieza-
le\nych q = [ q1 L qk ]T przyjmują postać
&& &
M(q)q + d(q,q,t) = f (q) + BÄ Ä , (5.1)
która jest formÄ… równaÅ„ (4.20) dla a" 0 . W równaniach tych u a" Ä = [Ä1 L Ä ]T , a macierz
r k
BÄ jest odwracalnÄ… macierzÄ… kwadratowÄ… o wymiarze k × k . Z zastosowaniem charakterystyk
& &&
ruchu qd (t) , qd (t) i qd (t) z pomiarów kinematycznych, przebiegi w czasie momentów ste-
Ä
rujÄ…cych (t) podczas tego ruchu liczne sÄ… jako
d
&& &
Ä (t) = BÄ-1(M(qd )qd + d(qd ,qd ,t) - f (qd ) ). (5.2)
d
Przykład 5.1. Model kończyny górnej sterowany za pomocą wypadkowych momentów
Å›
sił mię niowych. W odniesieniu do ogólnej postaci B zdefiniowanej w równaniu (4.20) oraz
do treÅ›ci przykÅ‚adów 4.1 i 4.3, Å‚atwo pokazać, \e 2× 2 -wymiarowa macierz BÄ ma postać
1
îÅ‚ -1
Å‚Å‚
BÄ = DTBÄ = , (5.3)
ïÅ‚0 1 śł
ðÅ‚ ûÅ‚
gdzie BÄ i D sÄ… zdefiniowane odpowiednio w równaniach (4.8) i (4.22). Macierz BÄ jest za-
tem zawsze odwracalna, co czyni rozwiÄ…zanie (5.2) oczywistym i jednoznacznym dla M , d
oraz f wygenerowanym zgodnie z opisem zawartym w przykładzie 4.3.
Na tle powy\szego trywialnego przypadku, problem wyznaczania wypadkowych mo-
mentów sił mięśniowych w modelu całego ciała człowieka nie jest ju\ tak oczywisty. Liczba
l = lÄ momentów sterujÄ…cych u a" Ä = [Ä1 L Äl ]T jest w tym przypadku mniejsza od liczby k
współrzędnych uogólnionych q = [q1 L qk ]T , będącej liczbą stopni swobody układu w fazie
lotu (braku kontaktu z podło\em). Dla rozwa\anych w tej pracy modeli płaskich ró\nica ta
wynosi k - l = 3 , co pojęciowo odpowiada trzem niesterowanym w fazie lotu wypadkowym
ruchom: ruchowi środka masy układu (dwa stopnie swobody) oraz wypadkowemu ruchowi
obrotowemu członów względem środka masy. W fazie lotu układ modelujący ciało człowieka
jest więc zawsze globalnie nie w pełni sterowany (równie\, gdy sterowanie modelowane jest z
u\yciem sił mięśniowych). W fazie kontaktu z podło\em na układ nakładane są dodatkowe
więzy, których liczba, dla budowanych w tej pracy płaskich modeli, wynosi k - l = 3 .
Dla uÄ jako parametry sterowania, dynamiczne równania ruchu (4.20) przyjmÄ… postać
T
&& &
M(q)q + d(q,q) = f (q) - Cr (q) + BÄ Ä , (5.4)
r
6
gdzie k × l - wymiarowa staÅ‚a macierz BÄ ma bardzo prostÄ… budowÄ™ (patrz przykÅ‚ad 5.2). Za-
proponować mo\na co najmniej dwa sposoby rozwiązania zadania symulacji dynamicznej
odwrotnej dla tak sformułowanego modelu dynamiki układu. Rozwiązanie najprostsze polega
T
na posÅ‚u\eniu siÄ™ k × k - wymiarowÄ… (odwracalnÄ…) macierzÄ… [ BÄ M - Cr ], co prowadzi do
Ä
(t)
îÅ‚ Å‚Å‚
d
T
&& &
(5.5)
ïÅ‚ (t)śł = [ BÄ M - Cr (qd )]-1(M(qd )qd + d(qd ,qd ) - f (qd )).
ðÅ‚ rd ûÅ‚
& &&
W ruchu opisanym przez qd (t) , qd (t) i qd (t) wyznaczane są tym samym zarówno momenty
Ä
sił mięśniowych (t) zapewniające realizację tego ruchu jak i reakcje modelu z podło\em
d
(t) w tym ruchu. Wzór (5.5) stosować mo\na te\ zarówno wówczas, gdy występuje kon-
rd
takt ciała z podło\em jak i wówczas, gdy kontakt ten nie występuje (patrz rozdział 3.5). W
pierwszym przypadku, obliczane tą drogą przebiegi (t) mo\na ewentualnie skonfrontować
rd
z wartościami mierzonymi na platformie dynamometrycznej. W drugim przypadku, (t)
rd
obliczane z ostatnich trzech równań (5.5) powinny być (z pewną dokładnością) równe zeru,
co odpowiada omawianym wcześniej trzem niesterowanym ruchu modelu jako całości w fazie
lotu. Ewentualne du\e niezerowe wartości (t) z (5.5) w tej fazie ruchu świadczyć mogą o
rd
błędach popełnionych na etapie modelowania matematycznego lub niepoprawności u\ytych
& &&
charakterystyk kinematycznych qd (t) , qd (t) i qd (t) .
Innym sposobem rozwiÄ…zania zadania symulacji dynamicznej odwrotnej dla modelu
ciała człowieka opisanego równaniami dynamicznymi (5.4), zastosowanym m.in. w pracach
[Blax3, Cza], jest technika rzutowania tych równań odpowiednio do l - wymiarowej podprze-
strzeni sterowanej oraz k - l - wymiarowej podprzestrzeni niesterowanej w k - wymiarowej
&
przestrzeni liniowej odpowiadającej q [16-19,75,89]. O ile podprzestrzeń sterowana definio-
wana jest przez kolumny w k × l - wymiarowej macierzy BÄ , podprzestrzeÅ„ niesterowanÄ…
definiuje (k - l)× k - wymiarowa macierz AÄ speÅ‚niajÄ…ca warunek
AÄ BÄ = 0 , (5.6)
czyli macierz uzupeÅ‚nienia ortogonalnego do BÄ . FormuÅ‚a rzutowania równaÅ„ (5.4) do pod-
przestrzeni sterowanej i niesterowanej, analogiczna do formuły (2.52), jest następująca
T
îÅ‚ Å‚Å‚
BÄ M-1
T
&&
(M q + d - f + Cr r - BÄ Ä)= 0 . (5.7)
ïÅ‚ śł
AÄ
ðÅ‚ ûÅ‚
Pierwsze l z tych równań (rzut do podprzestrzeni sterowanej), odpowiadające pierwszym l
równaniom (5.5), prowadzi do
& &&
HÄ (q)uÄ = hÄ (q,q,q, ,t) , (5.8)
r
T
gdzie l × l - wymiarowa macierz HÄ = BÄ M-1BÄ jest z zaÅ‚o\enia odwracalna, l - wektor hÄ
T -1 T
&&
ma postać hÄ = BÄ (q + M (d - f + Cr r )). Ostatnie k - l równaÅ„ (5.7) (rzut do podprze-
7
strzeni niesterowanej), odpowiadające ostatnim k - l równaniom (5.5), prowadzi natomiast
do zale\ności na bezpośrednie wyznaczanie , a mianowicie
r
& && & & &&
(q,q,q,t) = [GÄ (q)]-1AÄ (f (q,q,t) - d(q,q) - M(q)q), (5.9)
r
T
gdzie GÄ = AÄ Cr jest (k - l)× (k - l) - wymiarowÄ… macierzÄ… odwracalnÄ…. Z zastosowaniem
& &&
qd (t) , qd (t) i qd (t) otrzymywane mogą być stąd (t) w analizowanym ruchu, a następnie
rd
z równania (5.8)
& &&
Ä (t) = [HÄ (qd )]-1hÄ (qd ,qd ,qd , ,t) . (5.10)
d rd
Rozwiązanie to jest równowa\ne otrzymywanemu z pierwszych l równań (5.5).
Przykład 5.2. Płaski model skoczka sterowany za pomocą wypadkowych momentów sił
Å›
mię niowych. W odniesieniu do zawartości przykładów 4.2 (rys. 4.5a) i 4.4 , mo\na pokazać,
\e 11×8 - wymiarowa macierz BÄ oraz 3×11-wymiarowa macierz Cr majÄ… postacie:
0 0 0 0 0 0 0 0
îÅ‚ Å‚Å‚
ïÅ‚ śł
0 0 0 0 0 0 0 0
ïÅ‚ śł
ïÅ‚ -1 0 0 0 0 0 0
śł
1
ïÅ‚ śł
0 1 -1 0 0 0 0 0
ïÅ‚ śł
ïÅ‚ śł
0 0 1 0 0 0 0 0
ïÅ‚ śł
BÄ = 0 0 -1 0 0 0 0 , (5.11)
ïÅ‚-1 śł
ïÅ‚ śł
0 0 0 1 -1 0 0 0
ïÅ‚ śł
0 0 0 0 1 -1 0 -1śł
ïÅ‚
ïÅ‚ śł
0 0 0 0 0 1 -1 0
ïÅ‚ śł
ïÅ‚ 0 0 0 0 0 0 1 0 śł
ïÅ‚ śł
0 0 0 0 0 0 0 1
ðÅ‚ ûÅ‚
1 0 l1 cosÕ1 l2 cosÕ2 l3 cosÕ3 0 0 0 0 0 0
îÅ‚ Å‚Å‚
Cr = -ïÅ‚0 -1 - l1 sinÕ1 - l2 sinÕ2 - l3 sinÕ3 0 0 0 0 0 0śł . (5.12)
ïÅ‚ śł
ïÅ‚ śł
0 1 0 0 0 0 0 0ûÅ‚
ðÅ‚0 0 0
T
Pozwala to sformuÅ‚ować 11×11- wymiarowÄ… macierz [ BÄ M - Cr ] , która (jak mo\na spraw-
dzić numerycznie) jest zawsze odwracalna. Schemat (5.5) symulacji dynamicznej odwrotnej
jest tym samy dobrze okreÅ›lony. Macierz AÄ o wymiarze 3×11, speÅ‚niajÄ…cÄ… warunek (5.6),
AÄ BÄ = 0 , mo\na zaproponować jako
1 0 0 0 0 0 0 0 0 0 0
îÅ‚ Å‚Å‚
ïÅ‚0
Ar = 1 0 0 0 0 0 0 0 0 0śł . (5.13)
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚0 0 1 1 1 1 1 1 1 1 1ûÅ‚
8
T
Aatwo stÄ…d sprawdzić, \e 3× 3 - wymiarowa macierz GÄ = AÄ Cr jest zawsze odwracalna,
det(GÄ ) = -1, co gwarantuje dobre uwarunkowanie zale\noÅ›ci (5.9) na wyznaczanie reakcji z
T
podÅ‚o\em. Z zaÅ‚o\enia odwracalna jest te\ 8×8 - wymiarowa macierz HÄ = BÄ M-1BÄ w za-
le\ności (5.10) na wyznaczenie sterowania układem (wypadkowych momentów sił sterują-
cych) w analizowanym ruchu.
5.4. PROBLEM STEROWANIA NADMIAROWEGO ZA POMOC SIA MIÅšNIOWYCH
Modele biomechaniczne sterowane za pomocą sił mięśniowych są najczęściej wyko-
rzystywanym przybli\eniem opisu aparatu ruchu człowieka (lub ewentualnie innych istot \y-
wych), szczególnie kończyn górnych i dolnych. W zadaniach symulacji dynamicznej odwrot-
nej, podstawowym problemem do pokonania jest w takich przypadkach jest nadmiarowość
sterowania ruchem liczba mięśni biorąca udział w realizacji ruchu w poszczególnych sta-
wach zawsze przewy\sza liczbę stopni swobody w danym stawie (niektóre mięśnie oddziałują
te\ na więcej ni\ jeden staw). Dystrybucji wyliczanych (jednoznacznie) wypadkowych mo-
mentów sił mięśniowych w stawach na poszczególne mięśnie dokonuje się zwykle z zastoso-
waniem odpowiednich kryteriów optymalizacyjnych [1,36,50,90,124,134,141,147]. Dla zilu-
strowania tego problemu posłu\my się opisanymi w poprzednim rozdziale modelami kończy-
ny górnej oraz człowieka wykonującego skok z równoległym prowadzeniem kończyn.
PrzyjmujÄ…c naprÄ™\enia mięśni jako parametry sterowania, u a" à = [Ã1 L à ]T , gdzie
l
l = lF jest liczbą mięśni, wprowadzone w równaniu (5.1) dynamiczne równania ruchu koń-
czyny we współrzędnych niezale\nych q = [q1 L qk ]T modyfikują się do postaci
&& & &
M(q)q + d(q,q,t) = f (q,q,t) + BÃ (q) Ã . (5.14)
Wymiar macierzy Bà wynosi w tym wypadku k × l i jest to macierz prostokÄ…tna o wiÄ™kszej
liczbie kolumn ni\ wierszy, l = lF > k . Dla danej chwili czasu t, dla której znane są qd (t) ,
Ã
& &&
qd (t) i qd (t) , bezpośrednie/jednoznaczne wyznaczenie (t) z równania (5.14) jest więc
d
niemo\liwe. Rozwiązanie mo\na otrzymać z zastosowaniem zagadnienia optymalizacyjnego
(optymalizacji statycznej [50,124, 134,141,147]), które formułuje się w postaci
minimalizuj J (Ã)
Å„Å‚
ôÅ‚
Ã
&& &
tak, by BÃ (qd) = M(qd)qd + d(qd ,qd ,t) - f (qd ) , (5.15)
òÅ‚
à à Ã
ôÅ‚
oraz d" d"
ół min max
à Ã
gdzie J (Ã) jest odpowiednio dobranÄ… funkcjÄ… celu, a i sÄ… dopuszczalnymi fizjolo-
min max
gicznie minimalnymi i maksymalnymi wartościami naprę\eń w mięśniach. W ten sposób
Ã
&&
znajdowane sÄ… (t) , które minimalizujÄ… J (Ã) , speÅ‚niajÄ… równość Bà à = M qd + d - f oraz
d
à Ã
mieszczą się w zakresie à min d" d" .
max
9
W literaturze zdefiniowanych zostało wiele funkcji celu, przeglądu których dokonuje
siÄ™ m.in. w pracach [50,134]. Najpowszechniej stosowanym jest kryterium Crowninshielda i
Branda [36], które zapewnia dającą się uzasadnić fizjologicznie dystrybucję sił mięśniowych
poprzez minimalizację sumy podniesionych do potęgi p naprę\eń mięśni
l
p
J = , (5.16)
"Ã j
j=1
gdzie wartość wykładnika p jest zwykle równa dwa lub trzy. Inne mo\liwość to na przykład
l l
p p
J = Aj ) = , (5.17)
"(Ã j "Fj
j =1 j=1
w której wykorzystywane są bezpośrednio siły mięśniowe Fj = à Aj ( Aj jest przekrojem
j
fizjologicznym mięśnia), oraz
l l
p p
J = Aj v ) = , (5.18)
"(Ã j j "N j
j=1 j=1
gdzie vj jest prędkością skracania mięśnia, a minimalizowana jest suma podniesionych do
potęgi p mocy mięśni.
W pracy [148] (patrz równie\ [147]) pokazuje się, \e rozwiązanie zadania optymaliza-
à à Ã
cyjnego (5.15), bez ograniczenia d" d" i z zastosowaniem funkcji celu (5.16) dla
min max
p = 2 , jest równowa\ne (jednoznacznemu) rozwiązaniu analitycznemu
Ã
&& &
(t) = BÃ (qd )[M(qd ) qd + d(qd ,qd ,t) - fg (qd )], (5.19)
d
T T
gdzie BÃ = BÃ (BÃ BÃ )-1 jest l × k - wymiarowÄ… macierzÄ… pseudo-odwrotnÄ… [30] do (prosto-
kÄ…tnej) k × l - wymiarowej macierzy BÃ . WadÄ… tego prostego i numerycznie efektywnego
Ã
rozwiązania jest fakt, \e otrzymywane tą drogą składowe (t) mogą mieć zarówno dodatnie
d
jak i ujemne wartości. Wartości ujemne przeczą zasadzie działania sił mięśniowych jako skra-
cających mięsień. Jak zasugerowano w pracy [147], problem ten mo\na rozwiązać poprzez
Ã
kolejne eliminowanie ujemnych składowych (t) , zerując je lub zastępując dopuszczalnymi
d
wartościami minimalnymi (co pojęciowo odpowiada stanowi braku aktywności danego mię-
śnia), a następnie uruchamiając procedurę obliczeń (5.19) kolejny raz tak długo, a\ rozwiąza-
Ã
nie (t) składać się będzie tylko ze składowych dodatnich (włączając narzucone wartości
d
Ã
minimalne). Podobnie postępować mo\na, gdy rozwiązanie (t) otrzymane z (5.19) zawie-
d
rać będzie składowe przekraczające maksymalne wartości dopuszczalne, zastępując je w tym
Ã
wypadku tymi wartościami maksymalnymi. Tak modyfikowane rozwiązanie (t) otrzymy-
d
à à Ã
wane za pomocą procedury (5.19) spełniać będzie warunek d" d" . Taki sposób roz-
min max
wiązania zadania wyznaczania sił mięśniowych zastosowano m.in. w pracy [26], otrzymując
wyniki to\same z rozwiÄ…zaniami za pomocÄ… wspomnianego zadania optymalizacyjnego.
10
Å›
Przykład 5.3. Wyznaczanie sił mię niowych w kończynie górnej. Rozwa\my jeszcze raz
model kończyny górnej o k = 2 stopniach swobody, sterowany za pomocą lF = 8 sił mię-
śniowych pokazanych na rysunkach 3.3, 3.11b oraz jeszcze raz na rys. 5.5a. W modelu tym
przedramię wraz z dłonią i trzymanym w niej cię\arem (usztywniony nadgarstek), traktuje się
jako jeden człon. Rejestrowany ruch polega na dynamicznym zgięciu i wyprostowaniu koń-
czyny (rys. 5.3), którego charakterystyki kinematyczne pokazane zostały na rys. 5.4. Dyna-
miczne równania ruchu ukÅ‚adu we współrzÄ™dnych niezale\nych q = [Õ1 Õ2]T (rys. 5.5a) majÄ…
postać (5.14), a sposób ich generowania został zilustrowany w przykładach 4.1 i 4.3. U\yte do
obliczeń dane masowo-geometryczne modelu (tab. 5.1) oszacowano dokonując stosownych
pomiarów antropometrycznych [156], natomiast współrzędne miejsc przyczepów mięśni oraz
dane dotyczące promieni działania poszczególnych mięśni w stawach łokciowym (rys. 3.7) i
barkowym przyjęto kierując się wskazówkami zawartymi w literaturze [52,147,149,156].
1 1
a) b)
S S
F8
( )
y t F7 F8
Sd
F7
F6 F6 2
( )
x t
Sd
F5
F3 F3 F6 F3
F4
2
1
3
4
F2
F4 z4
F6
F1
F1
4
F5
F2
F1
E
E F2 3
2
1 z3
F1 1
Y
2
2
X
Rys. 5.5. Siły mięśniowe w kończynie górnej oraz reakcje w stawie łokciowym
m l JC
¾C ·C
człon
kg m kg m2
m m
S1 3.2 0.311 0.184 -0.016 0.035
S2 5.2 0.295 0.256 -0.005 0.022
Tab. 5.1. Podstawowe dane masowo-geometryczne modelu kończyny górnej
Naprę\enia w poszczególnych mięśniach podczas zarejestrowanego ruchu oszacowano
stosujÄ…c procedurÄ™ optymalizacyjnÄ… (5.15) oraz funkcjÄ™ celu (5.15) dla p = 2 , przeliczajÄ…c je
następnie na wartości sił mięśniowych. Wyniki tych obliczeń, ograniczone do mięśni partycy-
pujÄ…cych w realizacji ruchu w stawie Å‚okciowym, pokazane sÄ… na rys. 5.6 odpowiednio dla
zginaczy i prostowników w stawie łokciowym. Wyró\nić mo\na trzy wyrazne fazy ruchu. W
pierwszej, trwającej 0.45 s , cię\ar podciągany jest dynamicznie do góry (patrz równie\ cha-
11
rakterystyki kinematyczne na rys. 5.4) i ruch ten realizowany jest za pomocÄ… zginaczy, ilo-
ściowo przede wszystkim za pomocą mięśni m4 (biceps briachii) oraz m3 (briachialis), a cha-
rakterystyczne obcięcie wartości F4 w środkowej części tej fazy ruchu jest wynikiem nało-
à Ã
\onego ograniczenia d" . W fazie drugiej, trwajÄ…cej od 0.45 s do 0.8 s , ruch przedra-
max
mienia z cię\arem trzymanym w dłoni jest wyhamowywany, co wymaga zaanga\owania pro-
stowników, a przede wszystkim m6 (triceps briachii c.met.lat). W fazie tej zginacze są nieak-
tywne, podobnie jak prostowniki w fazie pierwszej. Sytuacja znów się zmienia w fazie trze-
ciej, gdy ruch opadajÄ…cego przedramienia z ciÄ™\arem jest wyhamowywany za pomocÄ… zgina-
czy, a prostowniki sÄ… nieaktywne.
a) 600
F4
400
F4
F3
200
F3 F2
F1
F2
0
F1 0.4 0.8 1.2 1.6
t [s]
0
b)
600
F6
400
F5
200
0
t [s]
0 0.4 0.8 1.2 1.6
Rys. 5.6. Przebiegi sił mięśniowych w kończynie górnej: a) zginacze w stawie łokciowym,
b) prostowniki w stawie Å‚okciowym
800
2+ 2
3 4
400
4
0
-400
3
-800
t [s]
0 0.4 0.8 1.2 1.6
Rys. 5.7. Przebiegi składowych poziomej i pionowej reakcji w stawie łokciowym
12
i
F
[N]
i
F
[N]
i
[N]
Siły rozwijane w mięśniach generują momenty sił w stawach wymagane dla realizacji
badanego ruchu. Wpływają jednocześnie na wartości sił reakcji w stawach. Reakcje te, ogra-
niczone do stawu łokciowego E (rys. 5.5b) wyliczyć mo\na z zale\ności (4.29), która przyj-
& && 2 T &&
mie postać (q,q,q,Ã,t) = E [f + Bu - M (Dq + Å‚)] (patrz przykÅ‚ad 4.5), gdzie skÅ‚adowe tej
2
zale\ności opisane zostały w przykładach 4.1 i 4.3, natomiast macierz E jest utworzona z
trzeciej i czwartej kolumny macierzy E zdefiniowanej w równaniu (4.33). Uzyskane przebiegi
pokazanych na rys. 5.5b składowych poziomej ( 3 ) i pionowej ( 4 ), a tak\e wypadkowej siły
2
reakcji w stawie Å‚okciowym ( 3 + 2 ) podczas analizowanego ruchu pokazane sÄ… na rys.
4
5.7. Reakcje te osiągają wartości 800 N w fazach pierwszej (dynamiczne zginanie przedra-
mienia) i drugiej (hamowanie przedramienia w górnym poło\eniu). Poziom reakcji w trzeciej
fazie ruchu jest zdecydowanie mniejszy.
Zagadnienie optymalizacyjne (5.15) sformułowane zostało dla modelu kończyny górnej
zaczepionej przegubowo do poruszajÄ…cego siÄ™ znanym ruchem stawu barkowego (model ten
mo\e być bezpośrednio rozszerzony na przypadek wyizolowanej kończyny dolnej zaczepio-
nej do poruszajÄ…cego siÄ™ znanym ruchem stawu biodrowego). RozwiÄ…zanie tego zagadnienia
prowadzi do wyznaczenia sił w poszczególnych mięśniach modelu wymaganych dla realizacji
analizowanego ruchu. Porównując modele dynamiczne (5.1) oraz (5.14), w których sterowa-
nie realizowane jest odpowiednio za pomocÄ… momentów u a" Ä = [Ä1 L Ä ]T oraz naprÄ™\eÅ„ w
k
Ã
poszczególnych mięśniach u a" = [Ã1 L à ]T , wektor uogólnionych siÅ‚ sterujÄ…cych zapisy-
l
wany jest na dwa sposoby
fu = BÄ Ä = BÃ Ã , (5.20)
gdzie macierze BÄ i BÃ majÄ… wymiary odpowiednio k × k oraz k × l ( 2× 2 oraz 2×8 dla
Ä
analizowanego tu modelu kończyny górnej). Zakładając, \e (t) wymagane w danej chwili
d
czasu t dla realizacji analizowanego ruchu wyznaczono zgodnie z procedurÄ… (5.2), zagadnie-
nie optymalizacyjne (5.15) mo\na zapisać w następującej równowa\nej postaci
minimalizuj J (Ã)
Å„Å‚
ôÅ‚
Ã
tak, by BÃ (qd) = BÄ Ä d , (5.21)
òÅ‚
à à Ã
ôÅ‚
oraz d" d"
ół min max
Ã
Uzyskane tak przebiegi (t) będą identyczne z tymi uzyskanymi jako rozwiązanie (5.15). W
d
obu przypadkach dokonywana jest bowiem to\sama dystrybucja wypadkowych momentów
sił mięśniowych w stawach na poszczególne siły mięśniowe.
Idea zmodyfikowanego zagadnienia optymalizacyjnego (5.21) jest szczególnie u\y-
teczna przy wyznaczaniu sił mięśniowych w modelach mięśniowo-szkieletowych całego ciała
człowieka. Zilustrujmy ten problem na przykładzie rozwa\anego wcześniej płaskiego modelu
człowieka wykonującego wyskok w płaszczyznie strzałkowej (rys. 4.5). Dynamiczne równa-
13
nia ruchu modelu w k =11 współrzÄ™dnych uogólnionych q = [ xH yH Õ1 L Õ9 ]T , sterowa-
nego za pomocÄ… lÄ = 8 wypadkowych momentów siÅ‚ mięśniowych we wszystkich stawach
Ä Ä Ä
2 T 2 2 T
u a" Ä = [Ä1 L Ä8]T majÄ… symbolicznÄ… postać (5.4). StosujÄ…c zapis = [ ]T , gdzie
2 2 2
Ä = [Ä1 Ä Ä3]T sÄ… momentami w stawach koÅ„czyn dolnych, a Ä = [Ä Ä5 Ä6 Ä7 Ä8]T gro-
2 4
madzi wypadkowe momenty sił mięśniowych w pozostałych stawach, wektor uogólnionych
sił sterujących zapisać mo\na wariantowo jako
2 2 2 2 2 2
fu = BÄ Ä = BÄ Ä + BÄ Ä , (5.22)
2 2 2 2 2 2
gdzie BÄ = [ BÄ M BÄ ] , a wymiary macierzy BÄ i BÄ sÄ… odpowiednio 11× 3 i 11× 5 . Przyjmu-
jąc model mieszany sterowania za pomocą lF =15 sił mięśniowych w stawach kończyn dol-
Ã
nych (wybierajÄ…c naprÄ™\enia tych mięśni = [Ã1 L Ã15]T jako parametry sterowania) oraz
2 2 2 2
za pomocÄ… lÄ = 5 wypadkowych momentów siÅ‚ mięśniowych Ä w pozostaÅ‚ych stawach (rys.
2 2 T
4.5b), u = [ÃT Ä ]T = [Ã1 L Ã15 Ä L Ä8]T oraz l = lF + lÄ = 20 , dynamiczne równania
4
ruchu modelu modyfikujÄ… siÄ™ do symbolicznej postaci
T
&& &
M(q)q + d(q,q) = f (q) - Cr (q) + B(q)u . (5.23)
r
W równaniach tych, wektor uogólnionych sił sterujących zapisany jest jako
2 2 2 2
fu = B u = BÃ Ã + BÄ Ä , (5.24)
a wymiary macierzy B , BÃ i BÄ 2 sÄ… odpowiednio k × l , k ×lF i k × lÄ 2 (11× 20 , 11×15 i
11× 5 ). PorównujÄ…c równania (5.22) i (5.24) otrzymujemy
2 2
BÄ Ä = BÃ Ã , (5.25)
2
co wyra\a fakt, \e efekt momentów sterujÄ…cych Ä w stawach koÅ„czyn dolnych zastÄ™powany
Ã
jest efektem poszczególnych siÅ‚ mięśniowych (naprÄ™\eÅ„ w mięśniach = [Ã1 L Ã15]T ).
Ä Ä
2
Dysponując (t) , a więc równie\ (t) , otrzymanymi jako rozwiązania równania (5.5) lub
d d
& &&
(5.8) dla ruchu opisanego przez qd (t) , qd (t) i qd (t) , zagadnienie optymalizacyjne pozwala-
Ã
jące oszacować (t) w kończynach dolnych sformułować mo\na następująco
d
minimalizuj J (Ã)
Å„Å‚
ôÅ‚
Ã
2 Ä2 d
tak, by BÃ (qd) = BÄ ; (5.26)
òÅ‚
à à Ã
ôÅ‚
oraz d" d"
ół min max
Procedura ta była z powodzeniem zastosowana w pracach [25,26] przy wyznaczaniu sił mię-
śniowych (i reakcji w stawach) kończyn dolnych podczas wyskoku pionowego. Ze względów
objętościowych wyniki te nie będą tu pokazywane.
W opisanych wy\ej modelach siła j-tego mięśnia wyra\ana była poprzez wartość na-
prę\enia à tego mięśnia (parametr sterowania), Fj = Ajà , gdzie Aj jest przekrojem fizjo-
j j
logicznym mięśnia. Jako parametr sterowania wybrany mo\e być te\ poziom aktywacji mię-
14
śnia a lub poziom jego pobudzenia u , co wymaga zastosowania wzoru (3.4) oraz dodatko-
j j
wo (3.5), gdy parametrem determinującym siłę jest u , a warunek à d" à d" à zastę-
j j min j j max
powany jest wówczas odpowiednio przez 0 d" a d"1 oraz 0 d" u d"1. W obu przypadkach wy-
j j
magane jest zastosowanie daleko bardziej dokładnego (zło\onego) modelu mięśnia oraz iden-
tyfikacji szeregu opisujących pracę mięśnia zale\ności i parametrów (co tylko w zarysie
omówione zostało w rozdziale 3.3).
Omówione zagadnienia optymalizacyjne wyznaczania sił mięśniowych nale\ą do grupy
zadań optymalizacji statycznej [1,3,50102,119,134], bazujących na procedurze symulacji dy-
namicznej odwrotnej dla wyznaczenia zbiorczych momentów sił mięśniowych w stawach i
ewentualnie reakcji z otoczeniem, a następnie dystrybucji tych wyliczonych momentów na
poszczególne mięśnie z zastosowaniem procedur optymalizacyjnych. Wskaznik jakości eks-
tremalizuje się tym samym sukcesywnie w kolejnych chwilach badanej czynności ruchowej.
Procedury te są relatywnie efektywne nawet dla du\ej liczby mięśni w budowanych modelach
biomechanicznych [117-119]. Inną spotykaną praktyką wyznaczania sił mięśniowych podczas
badanej czynności ruchowej jest optymalizacja dynamiczna [3,50]. Wykorzystuje ona proce-
dury symulacji dynamicznej prostej (całkowania równań ruchu), podczas której ruch układu
w określonym czasie wymuszany jest zadanymi przebiegami sterowania, na przykład zmien-
nymi w czasie poziomami aktywacji mięśni u(t) a" a(t) = [a1(t) L al (t)]T ). Rozwiązanie
numeryczne q(t) (odpowiedz układu) jest następnie porównywane z danymi pomiarowymi
qd (t) , a sterowanie a(t) jest iteracyjnie korygowane w celu minimalizacji q(t) - qd (t) z za-
stosowaniem odpowiedniego kryterium J (najczęściej jest to J = q(t) - qd (t) ). W oblicze-
niach narzucane są te\ odpowiednie ograniczenia na wartości sterowań, na przykład 0 d" a d"1.
Pomimo wielu swoich zalet, wadÄ… optymalizacji dynamicznej jest jej relatywnie niska efek-
tywność numeryczna powodowana wielokrotnie powtarzanym całkowaniem równań ruchu w
poszukiwaniu sterowania optymalnego. W wielu wypadkach oba rodzaje optymalizacji dajÄ…
podobne wyniki [3].
5.5. WERYFIKACJA POPRAWNOŚCI WYNIKÓW SYULACJI DYNAMICZNEJ
ODWROTNEJ
Wyniki symulacji dynamicznej odwrotnej mogą być obarczone wieloma niedokładno-
ściami. Istotnym ich zródłem są błędy pomiarów i odczytu charakterystyk kinematycznych
badanych czynności ruchowych, a następnie obróbka numeryczna tych charakterystyk polega-
jąca na aproksymacji charakterystyk dyskretnych charakterystykami ciągłymi (na poziomie
poło\eń, prędkości i przyspieszeń). Najwięcej zastrze\eń budzą jednak uproszczenia stosowa-
ne przy modelowaniu ciała człowieka niedoskonałe modele pracy mięśni oraz często trudne
do jednoznacznego sprecyzowania linie ich działania i miejsca przyczepu, ograniczenie liczby
modelowanych mięśni, uproszczenia przy modelowaniu połączeń w stawach czy pomijanie
15
y
podatności tkanek miękkich. ródłem niedokładności wyników symulacji numerycznych mo-
gą być ponadto stosowane zało\enia przy budowie algorytmów obliczeniowych, czego naj-
lepszym przykładem jest mo\liwy wariantowy wybór kryteriów optymalizacyjnych stosowa-
nych przy dystrybucji wypadkowych momentów sił mięśniowych w stawach na poszczególne
siły mięśniowe. Równie\ dane masowo-geometryczne stosowane w obliczeniach szacowane
są du\ym marginesem niepewności. Wszystko to sprawia, \e do wyników symulacji nume-
rycznych nale\y podchodzić z ograniczonym zaufaniem, traktując je tylko jako zawsze niedo-
skonałe przybli\enie modelowanych zjawisk.
Rys. 5.8. Salto w tył w pozycji łamanej wykonywane z pomocą batutu: przebiegi otrzymane
jako rozwiÄ…zanie zadania symulacji prostej (gruby szary kontur) oraz przebiegi
zmierzone (kontury w postaci cienkiej czarnej linii)
Prostą metodą oceny jakości wyników symulacji dynamicznej odwrotnej jest zastoso-
wanie rozwiÄ…zania ud (t) jako wymuszenia w zadaniu prostym dynamiki, stosujÄ…c ten sam
model dynamiczny i te same wartości początkowe jak w zadaniu odwrotnym, a następnie po-
równanie tak otrzymanych przebiegów ruchu z u\ytymi w zadaniu odwrotnym. Ró\nica mię-
dzy tym rozwiązaniem a charakterystykami otrzymywanymi z pomiarów będzie oczywiście
się pojawiać (i narastać) w czasie symulacji, wywołana głównie niedokładnościami przy po-
miarach kinematycznych i otrzymywaniu na tej podstawie stosowanych w zadaniu odwrot-
& &&
nym qd (t) , qd (t) i qd (t) , a tak\e, w mniejszym stopniu, błędami zaokrągleń przy całkowa-
16
niu numerycznym w zadaniu prostym. Nie powinien to być jednak proces gwałtowny. Szybko
narastające ró\nice między trajektoriami pierwotną (zarejestrowaną) i wtórną (obliczoną)
wskazują zwykle na błędy metodologiczne w modelowaniu i algorytmach obliczeniowych.
Przykład takiego porównania, pochodzący z prac [22,37], pokazany jest na rys. 5.8 dla mode-
lu sportowca wykonującego salto w tył w pozycji łamanej po wybiciu z batutu. Momenty sił
mięśniowych w stawach, wyliczone z zadania symulacji dynamicznej odwrotnej przy u\yciu
zmierzonych charakterystyk kinematycznych, zastosowane zostały jako wymuszenie ruchu
układem w zadaniu symulacji dynamicznej prostej (przy zgodnych warunkach początko-
wych). Efekt symulacji, podany z krokiem 0.08 s , widoczny jest w postaci konturów zazna-
czonych za pomocÄ… cienkich czarnych linii, na tle grubych szarych linii ilustrujÄ…cych zmie-
rzone poło\enia skaczącego. Jak widać, wyniki symulacji prostej są w przybli\eniu zgodne z
danymi wejściowymi do symulacji odwrotnej. Weryfikuje to korzystnie przyjęty model ma-
tematyczno-fizyczny oraz zastosowane algorytmy obliczeniowe przy rozwiÄ…zywaniu zagad-
nienia odwrotnego.
Porównywanie wyników symulacji dynamicznej prostej z zastosowaniem sterowania
otrzymanego z rozwiÄ…zania zadania symulacji dynamicznej odwrotnej, przy wykorzystaniu
tego samego modelu matematycznego układu biomechanicznego, weryfikuje jedynie jakość
procedur otrzymywania i obróbki charakterystyk kinematycznych oraz poprawność u\ytych
algorytmów obliczeniowych. Poprawność modelowania matematycznego mo\e być (w ogra-
niczonym zakresie) oceniona poprzez porównanie wyników symulacji dynamicznej odwrot-
nej otrzymanych przy wykorzystaniu dwu jakościowo ró\nych modeli matematycznych tego
samego modelu fizycznego, na przykład przy wykorzystaniu tradycyjnych współrzędnych
uogólnionych (niezale\nych) oraz współrzędnych naturalnych (zale\nych).
Najlepszym sposobem weryfikacji poprawności metod modelowania i symulacji nume-
rycznych są badania doświadczalne. W biomechanice zakres takich (bezinwazyjnych) badań
jest niestety ograniczony. Do badań takich nale\y pomiar reakcji z podło\em z wykorzysta-
niem platformy dynamometrycznej. Pochodzący z pracy [26] przykład porównania wyliczo-
nych i zmierzonych przebiegów reakcji pionowej wywieranej obunó\ na podło\e przez spor-
towca (o masie m =114.8 kg ) wykonujÄ…cego wyskok pionowy pokazany jest na rysunku 5.9.
Reakcje z podło\em = [ Rx Ry M ]T (rys. 3.14), a więc i reakcja pionowa Ry , liczone
r A
były z zale\ności (5.9), wykorzystując otrzymane na podstawie pomiarów charakterystyki
& &&
kinematyczne qd (t) , qd (t) oraz qd (t) , a wybicie siÄ™ z platformy dynamometrycznej i lÄ…do-
wanie na niej pozwoliło zmierzyć przebieg reakcji pionowej podczas skoku. Widać du\ą ja-
kościową i ilościową zgodność obu przebiegów. W fazie lotu mierzona reakcja pionowa jest
oczywiście równa zeru, niewielkie (bliskie zeru) wartości wyliczanych reakcji pionowych w
tej fazie jest prawdopodobnie wynikiem niedokładności wejściowych danych kinematycz-
nych. Zmierzony silny przyrost reakcji pionowej bezpośrednio po dotknięciu stopami podło\a
(platformy) wią\e się prawdopodobnie nie tylko z zarejestrowaną siłą dynamicznie ale rów-
17
nie\ z efektem dynamicznego wzbudzenia (krótkotrwałych, silnie tłumionych drgań) platfor-
my. Generalnie nale\y stwierdzić, \e przeprowadzony pomiar korzystnie weryfikuje popraw-
ność zbudowanego modelu matematycznego oraz jakość danych kinematycznych.
obliczenia
3
pomiar
2
1
0
t [s]
0 1 2 3
Rys. 5.9. Wyliczone i zmierzone przebiegi reakcji pionowej podczas wyskoku pionowego
Drugim przykładem bezinwazyjnych badań doświadczalnych jest elektromiografia po-
wierzchniowa (sEMG, od angielskiego: surface electromyography) [28,46,54,91,107]. Jest to
podstawowa metoda badawcza dostarczająca informacji o poziomie aktywacji mięśnia, pole-
gająca na rejestracji czynności elektrycznych mięśnia [92]. Mierzony na (naklejonych w od-
powiednich miejscach) elektrodach (rys. 5.3) potencjał polaryzacji, wynoszący w spoczynku
około 80 mV , ulega zmianom w zale\ności od stanu czynnościowego mięśnia. Zakłada się
przy tym, \e uzyskiwany sygnał EMG w trakcie czynności ruchowej przekłada się na poziom
aktywacji mięśnia [15,88,92,100,107]. Sygnały te wykorzystywane są do ustalania czasu ak-
tywacji mięśnia (kiedy rozpoczyna się i kończy pobudzenie mięśnia), ustalania sił rozwija-
nych w mięśniach (poprzez bezpośrednie przeło\enie poziomu sygnału EMG na wartość siły
[Van] lub korekcję metod optymalizacyjnych dystrybucji momentów sił mięśniowych na po-
18
y
R
[kN]
szczególne siły [33,60]), a tak\e określanie tempa zmęczenia mięśnia poprzez analizę widma
częstotliwości sygnału [34].
0.40 0.96
xS [m] yS [m]
0.36 0.95
0.32 0.94
0.28 0.93
t [s] t [s]
0 1 2 3 0 1 2 3
25 180
Õ1 [deg] Õ2 [deg]
0 90
-25 0
-50 -90
t [s]
t [s]
0 1 2 3
0 1 2 3
Rys. 5.10. Charakterystyki ruchu kończyny górnej uzyskane na podstawie pomiarów
1000
1000
m4 (BBC)
m1 (BR)
[µV]
[µV]
500
500
0 0
t [s] t [s]
0 1 2 3 0 1 2 3
Rys. 5.11. Sygnały EMG oraz na tej podstawie oszacowane okresy aktywności mięsni m1 i m4
0.40
0.12
[MPa] Ã4 (BBC)
Ã1 (BR)
[MPa]
0.30
0.08
0.20
0.04
0.10
0.00 0.00
t [s] t [s]
0 1 2 3 0 1 2 3
\
Rys. 5.12. Wyliczone wartości naprę eń w mięśniach m1 i m4 na tle sygnałów EMG
Zilustrujmy podstawowe badanie EMG dla ruchu kończyny górnej obcią\onej dodat-
kowym odwa\nikiem trzymanym w dłoni (rys. 3), wyniki rejestracji którego pokazano na rys.
5.10, gdzie xS (t) i yS (t) oznaczajÄ… zmiany poÅ‚o\eÅ„ stawu barkowego, a Õ1(t) i Õ2(t) odchy-
lenia odpowiednio ramienia i przedramienia od pionu (rys. 4.9). Przebiegi sygnału EMG z
elektrody naklejonej na brzusiec mięśnia m1 (brachioradialis) oraz uśredniony sygnał z elek-
19
trod naklejonych na dwa brzuśce mięśnia m4 (biceps brachii) pokazane są na rys. 5.11 (poka-
zano wartości bezwzględne sygnałów). Grubą linią szarą oszacowano tam te\ czasy aktywacji
mięśni, co zawsze obarczone jest pewną dozą niepewności/niedokładności. Na rys. 5.12 po-
kazano następnie wyliczone przebiegi naprę\eń w mięśniach m1 i m4, z zastosowaniem opi-
sanych wcześniej algorytmów symulacji dynamicznej odwrotnej i kryterium optymalizacyj-
nego (5.16) dla p = 3 , na tle uzyskanych sygnałów EMG. Dopatrzeć się mo\na jakościowej
zgodności obu tych rodzajów przebiegów przede wszystkim w pierwszej fazie ruchu, gdy
cię\ar podciągany jest dynamicznie do góry.
Symulacje komputerowe oraz elektromiografia powierzchniowa sÄ… w zasadzie jedyny-
mi bezinwazyjnymi technikami pozwalającymi oszacować czasy aktywacji i siły rozwijane w
poszczególnych mięśniach w trakcie czynności motorycznych człowieka. Jakość rozwiązań
numerycznych jest silnie uzale\niona od poprawności/szczegółowości stosowanych modeli
biomechanicznych, identyfikacji parametrów tych modeli, precyzji danych pomiarowych oraz
fizjologicznie uzasadnionych metod optymalizacji przy rozwiÄ…zywaniu zadania sterowania
nadmiarowego w stawach. Sygnały EMG dają unikalną szansę doświadczalnej weryfikacji
tych rozwiązań. Ich interpretacja nie zawsze jest jednak łatwa w praktyce, a same pomiary są
bardzo wra\liwe na szereg zaburzeń i sposób prowadzenia badań (niedokładne przyleganie
elektrod, złe miejsce ich przyklejenia, przemieszczanie się elektrod wraz ze skórą, pocenie się
ciała, zmęczenie mięśni, ... [46]). Jakość sygnału EMG zale\y te\ w du\ej mierze od dynami-
ki ruchu [54,92,107] oraz jakości aparatury, na której prowadzone są pomiary [91].
Ć
5.6. WRAśLIWOŚ WYNIKÓW OBLICZEC NA ZMIAN ZAAOśEC MODELOWYCH
I METOD SYMULACJI
Wyniki symulacji numerycznych pozwalające oszacować siły mięśniowe oraz reakcje
w stawach układów biomechanicznych modelujących badane czynności ruchowe silnie zale\ą
od cech (zało\eń upraszczających) budowanych modeli dynamicznych, adekwatności przyję-
tych danych masowo-geometrycznych, dokładności pomiarów kinematycznych i ich obróbki,
a tak\e u\ytych procedur obliczeniowych. Badanie wra\liwości rozwiązań na zmianę tych
parametrów pozwala ocenić margines popełnianych błędów. Wnioski płynące z takiej analizy
mogą te\ mieć istotne znaczenie dla doskonalenia metod modelowania i symulacji poprzez
ocenę, które czynniki mają największy wpływ na ilościowe i jakościowe zmiany uzyskiwa-
nych wyników, wysiłek modelowania i badań eksperymentalno-pomiarowych mo\na skupić
na aspektach najbardziej istotnych. Podejmowana w literaturze analiza wra\liwości tego typu
dotyczy przykładowo zmian przyjmowanych parametrów masowo-antropometrycznych [105,
109], wartości przekrojów fizjologicznych mięśni [20], umiejscowienia punktów zaczepu
mięśni [129], czy parametrów modelu układu mięsień-ścięgno [110]. Wiele uwagi poświęca-
ne jest te\ wpływowi obróbki (wygładzania) danych z pomiarów kinematycznych na wyniki
20
szacowanych sił mięśniowych [124,141]. Poni\ej analiza taka zostanie pokazana dla ró\nych
zało\eń modelowych oraz u\ytych kryteriów optymalizacyjnych w zastosowaniu do modelu
kończyny górnej analizowanym wcześniej w przykładzie 5.3.
600
a)
F4
F6
400
F4
F3
200
F3
F5
0
t [s]
0 0.4 0.8 1.2 1.6
b)
600
F6
F4
400
F4
F3 F5
200
F3
0
-200
t [s]
0 0.4 0.8 1.2 1.6
Rys. 5.13. Oszacowania sił mięśniowych dla realnego (a) i uproszczonego (b) modelu
działania mięśni m5 i m6
Pierwsza ilustracja polega na ocenie wpływu modelu działania sił mięśniowych, który,
jak pisano ju\ o tym w rozdziale 3.3, jest czynnikiem mającym fundamentalny wpływ na wy-
niki (wiarygodność) wyznaczanych sił mięśniowych, a w konsekwencji równie\ na wyzna-
czane siły reakcji w stawach. Problem ten podnoszony jest zresztą powszechnie w literaturze,
m.in. [108,142,153]. W modelu podstawowym linie działania sił mięśniowych zamodelowano
w sposób pokazany na rys. 4.3, tzn. z uwzględnieniem ewentualnego oplatania stawu przez
ścięgno mięśnia i skutkujące niezerowym promieniem działania siły względem stawu nieza-
le\nie od kąta stawowego. Działanie sił mięśniowych przebiega wówczas wzdłu\ odcinka
łączącego efektywne miejsca wprowadzenia tych sił. Wyniki obliczeń wyznaczania sił mię-
śniowych, prezentowane wcześniej w przykładzie 5.3 (rys. 5.6), ograniczono na rys. 5.13a do
czterech mięśni mających decydujący wpływ na realizację ruchu w stawie łokciowym E: zgi-
naczy m3 i m4 oraz prostowników m5 i m6 (rys. 3.3 oraz 5.5). Model uproszczono następ-
nie nie uwzględniając oplatania stawu łokciowego przez prostowniki m5 i m6, tzn. zastąpiono
pokazany na rys. 3.5 realny model działania tych mięśni przez model uproszczony pokazany
21
3
4
5
6
F
,
F
,
F
,
F
[N]
3
4
5
6
F
,
F
,
F
,
F
[N]
na rysunku 3.4. Efekt tego zabiegu na wyniki symulacji numerycznych pokazany jest na rys.
5.13b. Analizowany ruch kończyny górnej wymaga zaanga\owania prostowników m5 i m6
dla dynamicznego hamowania i następnie prostowania kończyny (z obcią\nikiem trzymanym
w dłoni) dla poło\enia zilustrowanego symbolicznie na rysunku 3.4b. Powoduje to gwałtow-
ny wzrost wymaganych sił w tych mięśniach, gdy ich promień działania względem stawu E
maleje do zera. W dalszej fazie, gdy linie działania tych sił przechodzą na prawo od E, reali-
zacja ruchu za pomocą dodatnich (ciągnących) sił mięśniowych jest niemo\liwa, a rozwiąza-
nie numeryczne istnieje tylko przy zało\eniu ujemnych sił mięśniowych (zarówno prostowni-
ków jak i zginaczy). Jest to oczywiście rozwiązanie pozbawione fizycznego sensu.
40 60
F4
F4
30
F4/2 + F4/2
F4/2 + F4/2
40
20
20
10
0
0
t [s]
0 0.4 0.8 1.2 1.6 t [s]
0 0.4 0.8 1.2 1.6
300 600
F4
F4
F4/2 + F4/2
F4/2 + F4/2
200
400
100
200
0
0
t [s]
t [s]
0 0.4 0.8 1.2 1.6
0 0.4 0.8 1.2 1.6
Rys. 5.14. Przebiegi wybranych sił mięśniowych dla modelu podstawowego i modelu,
w którym mięsień m4 podzielono na dwa identyczne aktony
Kolejna ilustracja dotyczy przypadku, gdy, w porównaniu do modelu podstawowego,
mięsień m4 (biceps briachii) został podzielony na dwa identyczne aktony działające identycz-
nie jak mięsień m4 i mające sumaryczny przekrój fizjologiczny taki jak dla m4. Przez ten za-
2
bieg układ o l = lF = 8 mięśniach modyfikowany jest do modelu o lF = 9 mięśniach. Na rys.
5.14 pokazano wyliczone przebiegi sił mięśniowych dla zginaczy w stawie łokciowym dla
obu tych modeli: podstawowego i zmodyfikowanego, oznaczone liniami odpowiednio ciÄ…-
głymi i przerywanymi. Pierwszą obserwacją jest, \e sumy wyznaczonych sił w dwu aktonach
podzielonego mięśnia m4 dla modelu o lF = 9 mięśniach są mniejsze od wyznaczanych sił w
2
mięśniu m4 dla modelu podstawowego o lF = 8 mięśniach. Jest to wynik minimalizacji sumy
kwadratów naprę\eń w mięśniach, zgodnie z kryterium (5.16) dla p = 2 . Dla pierwszego mo-
delu minimalizowane są bowiem dwa (wartościowo) du\e naprę\enia w dwu aktonach zastę-
22
1
2
F
[N]
F
[N]
3
4
F
[N]
F
[N]
pujących m4, w drugim wypadku jedno wartościowo du\e naprę\enie w mięśniu m4. Zmniej-
szony efekt działania podzielonego mięśnia m4 rekompensowany jest następnie przez zwięk-
szone wartości sił w pozostałych zginaczach.
60 90
p = 1 p = 1
p = 2 p = 2
40 60
p = 3 p = 3
20 30
0 0
t [s] t [s]
0 0.4 0.8 1.2 1.6 0 0.4 0.8 1.2 1.6
400 600
p = 1
p = 1
p = 2
p = 2
300
p = 3
400
p = 3
200
200
100
0
0
t [s]
t [s]
0 0.4 0.8 1.2 1.6
0 0.4 0.8 1.2 1.6
250 800
p = 1 p = 1
200
p = 2 p = 2
600
p = 3 p = 3
150
400
100
200
50
0
0
t [s]
t [s]
0 0.4 0.8 1.2 1.6
0 0.4 0.8 1.2 1.6
0 1000
p = 1
800
-200 p = 2
p = 3
600
-400
p = 1
400
p = 2
-600
p = 3
200
-800 0
t [s] t [s]
0 0.4 0.8 1.2 1.6 0 0.4 0.8 1.2 1.6
Rys. 5.15. Przebiegi sił mięśniowych oddziałujących na ruch w stawie łokciowym oraz reakcji
\
w stawie łokciowym liczone z zastosowaniem ró nych kryteriów optymalizacyjnych
Ostatnia ilustracja dotyczy wpływu wyboru kryterium optymalizacyjnego na wyniki
dystrybucji momentów sił mięśniowych w stawach na poszczególne siły mięśniowe. O ile w
modelu podstawowym zastosowano kryterium Crowninshielda i Branda (5.16) dla p = 2 ,
23
1
2
F
[N]
F
[N]
3
4
F
[N]
F
[N]
5
6
F
[N]
F
[N]
3
4
[N]
[N]
zagadnienie optymalizacyjne (5.15) rozwiązano z zastosowaniem tego kryterium równie\ dla
p =1 oraz p = 3 . Wyniki zestawione są na rys. 5.15. Jak widać, wpływ kryterium optymali-
zacyjne na szacowane wartości sił mięśniowych jest jakościowo i ilościowo bardzo istotny.
Dla p =1 zagadnienie optymalizacji jest liniowe, co powoduje, \e dla realizacji wymaganych
momentów sterujących w stawach anga\owane są przede wszystkim mięśnie działające na
du\ym ramieniu względem stawu oraz mające du\e przekroje fizjologiczne (zginacze m3 i m4
oraz prostownik m6 w odpowiednich fazach ruchu). Sytuacja jest odmienna dla p = 3 miÄ™-
śnie anga\owane są bardziej równomiernie, a przypadek p = 2 daje wyniki pośrednie. Jak
pokazano na dwu ostatnich wykresach na rys. 5.15, wpływ wyboru kryterium optymalizacyj-
nego na wartości reakcji w stawach jest natomiast znikomy, pomimo otrzymywanych ró\nych
dystrybucji momentów sterujących w stawach na poszczególne siły mięśniowe.
Podsumowując, wynik dystrybucji momentów sił mięśniowych w stawach na poszcze-
gólne siły mięśniowe jest bardzo wra\liwy na zało\enia stosowanych modeli biomechanicz-
nych oraz rodzaj stosowanego kryterium optymalizacyjnego. O podstawowym znaczeniu sÄ…
zało\enia dotyczące modelu działania sił mięśniowych, uwzględniające zło\onych charakter
oplatania torebek stawowych przez ścięgna mięśni. Bardzo istotną jest równie\ liczba mode-
lowanych mięśni, co powoduje, \e do uproszczeń polegających na zastępowaniu grup mięśni
pojedynczymi aktonami nale\y podchodzić bardzo ostro\nie. Istotne znaczenie będzie mieć
prawdopodobnie równie\ model pracy mięśnia oraz wybór parametru sterowania (bezpośred-
nio siła mięśniowa, naprę\enie mięśnia, poziom aktywacji mięśnia, poziom pobudzenia mię-
śnia). Te ostatnie zało\enia modelowe korespondują z wyborem kryterium optymalizacyjnego
w obliczeniach numerycznych, okazującym się mieć równie\ decydujący wpływ na wyniki.
Wedle doświadczeń autora, wpływ zmian parametrów masowych i antropometrycznych mo-
delu na te wyniki, szacowanych zawsze z pewnym marginesem błędu, jest zdecydowanie
mniejszy i tylko ilościowy. Zmiany zało\eń modelowych i kryterium optymalizacyjnego mają
wpływ równie\ jakościowy na wynik szacowania sił mięśniowych. Uszczegółowienie i ade-
kwatność budowanych modeli biomechanicznych oraz fizjologicznie uzasadniony dobór kry-
terium optymalizacyjnego mają zatem podstawowe znaczenie dla wiarygodności obliczeń
obcią\eń wewnętrznych (sił mięśniowych i reakcji w stawach).
24
Wyszukiwarka
Podobne podstrony:
23 ROZ warunki i tryb postępowania w spr rozbiórek obiekCZ1 roz 1 12matematyka roz odp7 ROZ warunki techniczne baz i stacji paliw [M G ][21 1196 ROZ warunki przy wprowadzaniu ścieków do wód lub do ziroz Vimmunologia molekularna roz 4 5Rozzad z roz31 ROZ rozbiórki obiektów bud metodą wybuchową [M I ][3zjawiska paranormalne i seks roz 1więcej podobnych podstron