Analisis dinámico inverso de la marcha humana

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Máster Máster en Ingeniería Industrial Analisis Dinámico Inverso de la Marcha Humana Autor: Eduardo Maza Ortega Tutor: Joaquín Ojeda Granja Dep. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Proyecto Fin de Máster Máster en Ingeniería Industrial Analisis Dinámico Inverso de la Marcha Humana Autor: Eduardo Maza Ortega Tutor: Joaquín Ojeda Granja Profesor Contratado Doctor Interino Dep. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Proyecto Fin de Máster: Analisis Dinámico Inverso de la Marcha Humana Autor: Eduardo Maza Ortega Tutor: Joaquín Ojeda Granja El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos Es te trabajo no ha sido fácil y no hubiera sido posible completarlo sin el apoyo de mi familia y mis amigos. A todos ellos, gracias. I Resumen En el presente trabajo se desarrolla un programa en Matlab para resolver un problema dinámico inverso aplicado a la marcha humana. El objetivo de este problema es obtener el valor de los momentos en cadera, rodilla y tobillo durante la marcha. Para ello se toma como referencia el software de OpenSim pero intentando llevar a cabo un código mas sencillo y que permita al usuario libertad para modificarlo sin necesidad de amplios conocimientos de programación. Destacan el uso de calculo simbólico en Matlab y de un conjunto mínimo de coordenadas de manera que no existan restricciones algebraicas sobre ellas [18]. En los dos primeros capítulos se desarrolla el trasfondo teórico necesario para comprender el problema. En el tercer capítulo se explica el modelo dinámico utilizado. Los capítulos cuarto y quinto suponen el núcleo del trabajo ya que desarrollan el problema cinemático y dinámico inverso presentando los resultados obtenidos en Matlab y su comparación con OpenSim. Por último, en el sexto capitulo se obtienen conclusiones y se comentan aspectos a mejorar. III 2Capítulo 1. Introducción Newton-Euler para resolver el problema dinámico inverso. La principal innovación de este trabajo radica en que presenta una técnica general que permite no tener que derivar manualmente las ecuaciones del movimiento como se había hecho hasta ahora, por ejemplo, en el trabajo de Winter [ 20 ]. En este sentido también hay que destacar el artículo de Riemer et al. [ 15 ], en el que se modela el cuerpo humano como un conjunto de trece segmentos articulados entre sí para calcular los momentos en el plano sagital de las diferentes articulaciones. Para ello se resuelven las ecuaciones de equilibrio de cada segmento de forma sucesiva, empezando por el pie, por lo que es necesario conocer las fuerzas de reacción con el suelo. Una de las principales limitaciones de este modelo es que es bidimensional y por tanto solo permite calcular los momentos de flexión-extensión. En el trabajo de Damsgaard et al. [ 4 ] se hace una revisión del software de simulación AnyBody Modeling System, el cual se utiliza para resolver el problema dinámico inverso de modelos musculoesqueléticos del cuerpo humano. Para ello se hace uso de una formulación general basada en la dinámica de sistemas multicuerpo en la que se usan coordenadas cartesianas y parámetros de Euler como coordenadas generalizadas. Dichas coordenadas están sujetas a ecuaciones de restricción que es necesario resolver antes de plantear las ecuaciones de equilibrio. El artículo de Delp et al. [ 6 ] presenta un modelo gráfico del miembro inferior y se puede considerar el antecesor de OpenSim. En el se modela el miembro inferior como un sistema multicuerpo formado por siete segmentos. Es más tarde, en 2007, cuando Delp et al. [ 5 ] presenta propiamente OpenSim en otro artículo. OpenSim es una plataforma de software para modelar humanos, animales, robots y el ambiente, y simular su interacción y movimiento. OpenSim dispone de una interfaz gráfica de usuario (GUI por sus siglas en ingles) para visualizar modelos y generar y analizar simulaciones. Una captura del aspecto de esta GUI se puede observar en la Figura 1.1. De esta manera, OpenSim constituye una herramienta de fácil manejo para realizar análisis dinámicos inversos de modelos del cuerpo humano. Simbody es el código detrás de OpenSim y surge de la necesidad de la investigación biomédica de tener un software de dinámica multicuerpo orientado específicamente a la simulación de sistemas biomecánicos como el cuerpo humano. Se trata de un motor dinámico pensado para programadores o desarrolladores de aplicaciones que no necesariamente tienen que ser expertos en dinámica [ 19 ]. El código es de acceso libre y se puede encontrar en https://simtk.org/home/simbody. Con este último trabajo se llega a la aparición de este software, que supone el punto de partida del presente trabajo. Figura 1.1 Aspecto de la GUI de OpenSim. Frente a estos trabajos, el presente utiliza un modelo tridimensional del aparato locomotor humano basado en el uso de un conjunto mínimo de coordenadas relativas. Esto permite que estas coordenadas no estén sujetas a ecuaciones de restricción, lo que provoca un ahorro computacional ya que no hay que resolverlas. Por otra parte, al plantear el problema de forma general, es decir, 1.2 Objetivos 3 aplicando las ecuaciones de equilibrio sobre el conjunto de sólidos que conforman el sistema no es necesario conocer el valor de las fuerzas de reacción con el suelo. Los dos códigos de Matlab programados resuelven el problema dinámico inverso emulando a OpenSim pero mediante un código mas sencillo y que permite una mayor personalización por parte del usuario. Aunque el código de OpenSim está publicado y se puede acceder a el para modificarlo, está escrito usando lenguaje C++, que esta muy optimizado para mejorar el rendimiento pero puede resultar difícil para un ingeniero que no este especializado en el. Por ello el presente trabajo se plantea el objetivo de emular a OpenSim pero escribiendo el código en Matlab, un lenguaje de programación mas orientado a la ingeniería y que se imparte y utiliza ampliamente en la ESI. Uno de los retos que se plantean es hacer comparable la velocidad de ejecución del código en Matlab con el código en C++ que utiliza OpenSim y que a priori es mas rápido Esto se debe a que C++ es un lenguaje de programación compilado frente a M (el lenguaje de programación de Matlab), que es un lenguaje interpretado. La diferencia entre ambos es que en el lenguaje compilado el código se traduce a lenguaje maquina en un paso previo, mientras que en el interpretado este paso se va haciendo a medida que se ejecuta el código. La ventaja del lenguaje compilado es que cuando se ejecuta el código ya esta en lenguaje maquina, por lo que la ejecución es mucho mas rápida. La desventaja del uso de lenguaje compilado es que hay que compilar el código cada vez que este cambia. 1.2 Objetivos El objetivo principal de este trabajo es resolver un problema dinámico inverso de la marcha humana. Para ello se generan dos códigos de Matlab que resuelven la cinemática y la dinámica inversa. En el primero de ellos se obtiene el valor de los ángulos en las articulaciones; los resultados obtenidos se usan en el segundo para obtener el valor de los momentos en ellas. El esquema general del trabajo se muestra en la Figura 1.2. Desde el punto de vista computacional, la resolución del problema cinemático inverso usa calculo numérico para resolver un problema de optimización. Por su parte, el problema dinámico inverso consiste en última instancia en la resolución de un sistema de ecuaciones lineal (mas sencillo que una optimización) pero con la dificultad añadida del uso de calculo simbólico para la obtención de las matrices dinámicas. El modelo biomecánico utilizado es el conocido como gait2392 y se presentará en el tercer capítulo. Este modelo se encuentra dentro del paquete de descargas de OpenSim. Como datos de partida se tienen todos los parámetros del modelo y un fichero también presente en el paquete de descargas de OpenSim con los resultados de un experimento de análisis de la marcha humana. CINEMÁTICA INVERSA DINÁMICA INVERSA 𝑞 𝑇 parámetros del modelo coordenadas de los marcadores Figura 1.2 Esquema global. 2 Conceptos Teóricos En este capitulo se explican algunos conceptos importantes de los que se hará uso a lo largo de este trabajo y que conviene tener claros para una buena comprensión del mismo. 2.1 Ejes y planos anatómicos Antes de nada, al hacer referencia continua al movimiento del cuerpo humano y puntos anatómicos de este, es conveniente introducir los ejes y planos anatómicos que establecen una base para poder hablar de la disposición de los miembros del cuerpo. En la Figura 2.1 se observan los ejes y planos anatómicos que definen la geometría del cuerpo, así como los ejes de coordenadas globales (fijos en el espacio) que se van a utilizar a lo largo de todo el trabajo: •Eje X: en la dirección de la marcha. •Eje Y: eje vertical hacia arriba. •Eje Z: producto vectorial de los anteriores. Los tres planos anatómicos principales son: •Plano coronal o frontal: divide al cuerpo en las mitades posterior y anterior. • Plano sagital: divide al cuerpo en las mitades derecha e izquierda. En este plano se dan los movimientos principales relacionados con la marcha. •Plano transversal o axial: divide al cuerpo en las mitades craneal y caudal. Estos planos se cortan entre si dando lugar a los ejes principales. Las direcciones de estos ejes servirán para describir las rotaciones que se producen en las diferentes articulaciones del cuerpo. •Eje anteroposterior: es el corte de los planos sagital y transversal. •Eje craneocaudal: es el corte de los planos sagital y coronal. •Eje latero-lateral: es el corte de los planos coronal y transversal. En la Figura 2.2 se presentan algunos términos que definen las posiciones de los puntos anatómicos respecto a los planos anatómicos Así, a lo largo del eje craneocaudal un punto es proximal frente a uno distal por estar mas cerca del centro del cuerpo, es decir, del plano transversal tal y como esta representado en la figura Figura 2.1. Así, la tibia es un hueso mas distal que el fémur, por ejemplo. Si nos referimos ahora al eje latero-lateral, un punto es medial frente a uno lateral si esta mas cerca del plano sagital. Así, el dedo gordo del pie está en una posición medial frente al pequeño, que esta en una posición lateral. Por ultimo, para la pareja anterior-posterior nos tenemos que fijar en si el 5 6Capítulo 2. Conceptos Teóricos Figura 2.1 Ejes anatómicos. Figura 2.2 Términos anatómicos. punto anatómico esta delante o detrás del plano coronal. Si está delante, como el pecho, la posición es anterior. Por el contrario, si se encuentra detrás, como la espalda, la posición es posterior. Otra pareja de términos que se va a utilizar a lo largo del trabajo es la de dorsal/ plantar. Dorsal hace referencia al empeine del pie frente a plantar que hace referencia a la planta del pie. 2.2 Ecuaciones cinemáticas 7 X Y Z ≡z x y z x ≡ x' y z' y' x' x'' z' y'≡ y'' z'' Figura 2.3 Secuencia de rotaciones. 2.2 Ecuaciones cinemáticas 2.2.1 Movimiento de un sólido En cinemática se define un sistema de coordenadas global, fijo en el espacio, que sirve como referencia. A cada solido se le acopla un sistema de coordenadas que se denomina local. La posición y orientación de dicho sistema local frente al global es lo que permite describir el movimiento del solido. Todo movimiento resulta de una combinación de un movimiento de traslación con uno de rotación. El movimiento de traslación hace referencia a como cambia la posición del origen de coordenadas local respecto al global. El de rotación, por su parte, hace referencia a como varia la orientación de los ejes locales respecto a los globales. Si se piensa en dos dimensiones, describir el movimiento resulta sencillo: la traslación tiene lugar en un plano y la rotación consiste en un giro alrededor de un eje perpendicular al plano. Si bien el movimiento de traslación en tres dimensiones es homologo al caso tridimensional, la definición de la orientación es algo mas complicado. Aunque existen diferentes métodos para esta definición, a lo largo de este trabajo se han utilizado dos de ellos: ángulos de Tait-Bryan y parámetros de Euler. Ángulos de Tait-Bryan Los ángulos de Tait-Bryan describen tres giros sucesivos que permiten pasar de la orientación del sistema global a la del sistema local. El primer giro se realiza a lo largo del eje Z global, como se aprecia en la Figura 2.3 Realizando este giro se obtiene un sistema de coordenadas intermedio x0−y0−z0 . El siguiente giro se produce alrededor del eje x0, como se aprecia en la Figura 2.3 El ultimo giro se produce alrededor del eje y00 , como se aprecia en la Figura 2.3 Este tipo de rotaciones se denominan intrínsecas ya que se producen de forma sucesiva alrededor de los ejes que ya se han rotado previamente. En el modelo utilizado todas las rotaciones descritas con ángulos de Tait-Bryan se realizan de la manera descrita, es decir, alrededor de los ejes z−x0−y00 . En ningún caso las rotaciones se producen alrededor de los ejes sin girar del sistema de coordenadas original. Parámetros de Euler Otro de los métodos utilizados para hallar la orientación de un solido en este trabajo esta basado en los parámetros de Euler. En este caso el paso de un sistema de coordenadas a otro se lleva a cabo mediante la rotación alrededor de un único eje. Para describir esta rotación es necesario conocer la 8Capítulo 2. Conceptos Teóricos dirección de dicho eje: v1,v2,v3 son las coordenadas de un vector unitario en la dirección de dicho eje; por otra parte, el angulo de rotación se denota como α . Con estos valores se construyen los cuatro parámetros de Euler: θ=    θ1 θ2 θ3 θ4    =    cos(α 2) v1sin(α 2) v2sin(α 2) v3sin(α 2)    (2.1) Figura 2.4 Parámetros de Euler. Ecuación cinemática fundamental A partir de los ángulos de Tait-Bryan o de los parámetros de Euler se construye la matriz de orientación del sólido. Para obtener el vector con las coordenadas de un punto cualquiera del solido expresadas en el sistema global hace falta conocer esta matriz de orientación y la posición del origen de coordenadas locales expresado en el sistema de coordenadas global: r=R+Au (2.2) Siendo u el vector de posición del punto expresado en el sistema de coordenadas locales, A la matriz de orientación del solido, R la posición del origen de coordenadas locales expresado en coordenadas globales y r el vector de posición de dicho punto expresado en el sistema de referencia inercial (coordenadas globales). 2.2.2 Movimiento de un conjunto de sólidos Hasta ahora se ha estado centrado en conocer la posición y orientación de un solido en particular, sin embargo, el aparato locomotor humano se trata como un conjunto de sólidos articulados entre sí. Para poder obtener la posición de cada uno de estos sólidos se utilizará un conjunto de coordenadas relativas. Este tipo de coordenadas describen la posición y orientación de cada cuerpo respecto al cuerpo que le precede en la cadena cinemática. La ventaja del uso de estas coordenadas radica en que no se imponen restricciones sobre ellas, lo que permite no tener que resolver ningún tipo de ecuaciones de restricción. 2.3 Ecuaciones dinámicas 9 En el modelo utilizado en este trabajo el origen de la cadena cinemática es la pelvis. La cinemática de este cuerpo se define a partir de la posición del origen del sistema de coordenadas asociado a él y de los tres ángulos de Tait-Bryan que definen su matriz de rotación asociada. A partir de estos datos se puede obtener la posición de cualquier punto que se encuentre en la pelvis: rpelvis =Rpelvis +Apelvisupelvis (2.3) Como se ha dicho, la expresión de Rpelvis y Apelvis se obtienen a partir de la posición del origen del sistema de coordenadas asociado a él y de los tres ángulos de Tait-Bryan, que constituyen las coordenadas generalizadas asociadas a la pelvis. Por tanto, la ecuación (2.3) se escribe como: rpelvis(q) = Rpelvis(q)+Apelvis(q)upelvis (2.4) La matriz de rotación del fémur, el siguiente en la cadena cinemática, se define a partir de la orientación relativa de este solido respecto a la pelvis. Esta orientación se define a partir de la matriz de giro de la cadera, que es la articulación entre ambos sólidos. Por tanto la matriz de rotación del fémur se calcula de la siguiente manera: Afemur =ApelvisAcadera (2.5) El origen del sistema de coordenadas asociado al fémur también se define de manera relativa al de la pelvis, ya que se encuentra en la cadera, la articulación que une ambos sólidos. Con la posición y la matriz de rotación asociada se puede calcular la posición de cualquier punto perteneciente al fémur: rfemur =Rfemur +Afemurufemur (2.6) Esto se va haciendo sucesivamente para cada solido del modelo, desde el más proximal al más distal. Este resulta el principal inconveniente del uso de coordenadas relativas, ya que la posición de los segmentos mas distales se define en función de la posición de todos los segmentos que le anteceden en la cadena cinemática. De esta manera se complican las expresiones de las matrices de rotación de estos segmentos más distales, ya que dependen de un gran numero de coordenadas generalizadas. 2.3 Ecuaciones dinámicas Si en el apartado anterior se explicó la cinemática del modelo, es decir, como se describen los movimientos del mismo, en el presente se explica como se hallan las fuerzas que producen dicho movimiento. La dinámica del modelo se resuelve aplicando las ecuaciones de Newton-Euler al conjunto de sólidos que componen el modelo, como se muestra en la ecuación (2.7)           mi1 ... mi1 ¯ Ii ... ¯ Ii           | {z } ˆ M           ai G . . . ai G ¯ αi . . . ¯ αi           =           Fi . . . Fi ¯ Mi . . . ¯ Mi           | {z } ˆ Q +           0 . . . 0 −¯ ωi∧¯ Ii¯ ωi . . . −¯ ωi∧¯ Ii¯ ωi           | {z } ˆ Qv (2.7) 10 Capítulo 2. Conceptos Teóricos A continuación se explica el desarrollo matemático que permite escribir la ecuación (2.7) en función de las coordenadas generalizadas del modelo tal y como explica en el documento del profesor Escalona [ 7 ]. El primer paso consiste en expresar el vector de aceleraciones lineales y angulares en función de las coordenadas generalizadas, como muestra la ecuación (2.8).           ai G . . . ai G ¯ αi . . . ¯ αi           =           Hi . . . Hi ¯ Gi . . . ¯ Gi           |{z} L ¨ q+           hi . . . hi ¯ gi . . . ¯ gi           |{z} l ˙ q(2.8) Para poder escribir esta ecuación los valores de las coordenadas generalizadas obtenidos en la cinemática inversa se derivan de forma numérica para obtener ˙ q y ¨ q . La construcción de las matrices L y l requiere de los valores de Ri G y de ωi . La posición del centro de gravedad de los segmentos del modelo se extrae del modelo de cinemática inversa. Para las velocidades angulares, en un principio se siguió el documento del profesor Escalona [ 7 ], que presenta una forma general para realizar este cálculo a partir de las matrices de rotación: ˜ ¯ ωi=AiT˙ Ai(2.9) Sin embargo, al usar coordenadas relativas, la expresión de las matrices de rotación de los segmentos distales puede complicarse bastante. Por ejemplo: Atibia =ApelvisAcaderaArodilla (2.10) Es fácil entender la dificultad de tener que derivar una matriz de este tipo. La solución pasa por utilizar esta formulación tan solo para la pelvis y la combinación de velocidades angulares para el restos de segmentos. De esta manera se calcula la velocidad angular de cada segmento a partir de la velocidad angular del segmento anterior: ωi1=ωi i−1+ωi−1 1 (2.11) Donde ωi1 es la velocidad angular absoluta del sólido i, ωi−1 1 la velocidad angular absoluta del sólido i-1 y ωi i−1 es la velocidad angular relativa del sólido i respecto al sólido i-1. Esta velocidad angular es la asociada a la articulación que une los dos segmentos consecutivos. Su expresión es fácil de conseguir debido al uso de coordenadas relativas en el modelo. De esta manera, una vez obtenida la velocidad angular de la pelvis, la expresión de las velocidades angulares asociadas a espalda y cadera es inmediata ya que tienen la misma expresión pero con el conjunto de coordenadas generalizadas asociado a la articulación. A partir de Ri Gyωise obtienen las siguientes matrices: Hi=∂Ri G ∂q, ¯ Gi=∂¯ ωi ∂˙ q(2.12) hi=∂Hi˙ q ∂q, ¯ gi=∂¯ ωi ∂q(2.13) Con estas ecuaciones ya se pueden construir las matrices Lylde la ecuación (2.8) Sustituyendo la ecuación 2.8 en la (2.7) queda ˆ M(L¨ q+l˙ q) = ˆ Q+ˆ Qv(2.14) 2.4 El ciclo de marcha 11 Premultiplicando por LT M¨ q=Q+Qv(2.15) Donde M=LTˆ ML,Q=LTˆ Q,Q=LTˆ Q,Qv=LT(ˆ Qv−ˆ Ml˙ q)(2.16) De la forma en que esta definido el problema no existen fuerzas de reacción como tal ya que no se definen ecuaciones de restricción sobre las coordenadas generalizadas. Por tanto, las únicas fuerzas existentes son las asociadas a las coordenadas generalizadas y las gravitatorias: Q=Qgdl +Qgrav (2.17) Por tanto, la ecuación (2.15) se puede escribir de la siguiente manera: M¨ q=Qgdl +Qgrav +Qv(2.18) Despejando Qgdl en la anterior ecuación se llega a la ecuación (2.19) Qgdl =M(q)¨ q−Qv(q,˙ q)−Qgrav(q)(2.19) Donde: •Qgdl es el vector de fuerzas generalizadas asociadas a las coordenadas generalizadas •M(q)es la matriz de masas del sistema •¨ qes la derivada segunda de las coordenadas generalizadas •Qv(q,˙ q)es el vector de fuerzas cuadráticas en velocidad •Qgrav(q)es el vector de fuerzas gravitatorias En la ecuación (2.19) todos los términos del segundo miembro son conocidos ya que el movimiento es conocido y se trata de obtener Qgdl. Al existir una fuerza generalizadas por cada coordenada generalizada ( ¨q y Qgdl tienen en mismo tamaño) se tienen el mismo numero de ecuaciones que de incógnitas. Una vez obtenido Qgdl no queda mas que calcular la matriz con los vectores de posición de los puntos de aplicación de las fuerzas S para obtener el vector de momentos en las articulaciones T: T=S\Qgdl (2.20) 2.4 El ciclo de marcha La marcha humana se divide en dos fases principales: •Fase de apoyo, en la que el pie de referencia está en contacto con el suelo. Ocupa aproximadamente el 60% del tiempo de un ciclo. 18 Capítulo 3. Modelo Utilizado Figura 3.5 Movimiento del tobillo. para cada articulación se establece un numero de coordenadas generalizadas igual al número de grados de libertad que permite la articulación. Estas coordenadas, por tanto, describen el movimiento relativo entre los sólidos. Las coordenadas constituyen una forma matemática de caracterizar los movimientos de las articulaciones descritos en el apartado anterior. En total se tiene un conjunto de 23 coordenadas generalizadas se denota por q q q q=                                            pelvis_tx pelvis_ty pelvis_tz pelvis_tilt pelvis_list pelvis_rotation hip_f lexion_r hip_adduction_r hip_rotation_r knee_angle_r ankle_angle_r subtalar_angle_r mt p_angle_r hip_f lexion_l hip_adduction_l hip_rotation_l knee_angle_l ankle_angle_l subtalar_angle_l mt p_angle_l lumbar_extension lumbar_bending lumbar_rotation                                            (3.1) • Pelvis: las coordenadas pelvis tx, pelvis ty y pelvis tz describen la traslación en los ejes X, Y y Z globales, respectivamente. Se miden en metros. Las coordenadas pelvis tilt, pelvis list y pelvis rotation son los tres ángulos de Tait-Bryan que permiten pasar del sistema de coordenadas global al sistema de coordenadas local de la pelvis. 3.1 Modelo Original 19 • Espalda: esta articulación se modela como un par esférico, por lo que permite tres rotaciones y las tres coordenadas generalizadas son lumbar extensión, lumbar bending y lumbar rotation, que constituyen los tres ángulos de Tait-Bryan con los que se puede construir la matriz de rotación de esta articulación. Esta matriz de rotación permite pasar del sistema de coordenadas de la pelvis al del torso. • Cadera: se establece un par esférico, que permite tres rotaciones y al que se le asocian las tres coordenadas generalizadas hip flexion, hip adduction y hip rotation, que constituyen los tres ángulos de Tait-Bryan con los que se puede construir la matriz de rotación de esta articulación. Esta matriz de rotación permite pasar del sistema de coordenadas de la pelvis al del fémur. • Rodilla: en ella se produce una rotación entre fémur y tibia, pero el centro de rotación instantáneo no se mantiene fijo durante el movimiento. Se trata de un movimiento plano (en el plano que forman los ejes x y locales de la tibia) del centro de rotación instantáneo que depende del valor de la coordenada knee angle, como se muestra en la Figura 3.6. Más concretamente se utilizan splines que en función del valor de la coordenada knee angle dan los valores del movimiento en el plano. La rotación se produce alrededor del eje z del fémur. Esta articulación permite un grado de libertad, ya que conociendo el valor de la coordenada generalizada asociada (knee angle) se conoce tanto el movimiento del eje instantáneo de rotación como la rotación en sí. -2 -1.5 -1 -0.5 0 0.5 1 1.5 cm -48 -47.5 -47 -46.5 -46 -45.5 cm Origen de la tibia Con spline Sin spline Figura 3.6 Origen de la tibia en el plano sagital. • Tobillo: esta articulación se modela como un par de revolución entre tibia y talón cuya coordenada generalizada asociada es ankle angle. La matriz de rotación de esta articulación se construye a partir de parámetros de Euler donde la dirección del eje de revolución es un parámetro tomado del modelo y el ángulo de rotación es la coordenada generalizada ankle angle. • Articulación subtalar: similar a la articulación del tobillo pero entre talón y calcáneo. La coordenada generalizada asociada es subtalar angle. 20 Capítulo 3. Modelo Utilizado • Articulación metatarso-falángica: similar a la articulación del tobillo pero entre calcáneo y dedos. La coordenada generalizada asociada es mtp angle. Las rotaciones asociadas a las coordenadas generalizadas ankle angle, subtalar angle y mtp angle (metatarsophalangeal angle) son alrededor de unos ejes que especifica el modelo y que se aprecian en la figura 3.7 Figura 3.7 Ejes articulaciones del pie. Conociendo el valor de las 23 coordenadas generalizadas queda definido la posición de cualquier punto del modelo. Ademas estas 23 coordenadas son independientes entre si, no están sujetas a ninguna restricción. Por tanto se tiene un modelo de 23 grados de libertad. Como se ha visto en el apartado de orientación de sólidos en el espacio, la principal desventaja del uso de estas coordenadas es la expresión difícil de los últimos segmentos de la cadena cinemática; sin embargo, estas coordenadas permiten describir de manera univoca la posición de un sistema multicuerpo y son independientes entre sí. Como se vera en el capitulo dedicado a la cinemática inversa, esto permite resolver el problema sin tener que imponer ninguna restricción sobre ellas. 3.1.4 Datos de entrada del modelo Marcadores Se tiene un total de 31 marcadores que se utilizaran en la resolución de la cinemática y la dinámica inversa. •Torso –Sternum –R.Acronium –L.Acronium –TopHead •Pelvis –R.ASIS –L.ASIS 3.1 Modelo Original 21 –V.Sacral •Fémur –Thigh.Upper –Thigh.Front –Thigh.Rear •Tibia –Shank.Upper –Shank.Front –Shank.Rear •Calcáneo –Heel –Midfoot.Sup –Midfoot.Lat –Toe.Lat –Toe.Med –Toe.Tip Figura 3.8 Posición de los marcadores: vista anterior. Conviene aclarar que estos son los marcadores cuya posición se utiliza para resolver la cinemática inversa. La posición de estos marcadores aparece recogida en un fichero que constituye un dato de partida para el modelo. Sin embargo, para la construcción del modelo se entiende que se utilizo un menor número de ellos, como se explica en el siguiente apartado. 22 Capítulo 3. Modelo Utilizado Parámetros cinemáticos Los datos experimentales incluidos en el modelo se tomaron siguiendo las pautas dispuestas en el trabajo de Chand [ 9 ] pero usando un sujeto diferente al que se describe en dicho trabajo. El modelo general de OpenSim corresponde a un sujeto de unos 1.8 metros de altura y 75.16 kilogramos de peso. En el trabajo de Chand se realizaron mediciones sobre 8 sujetos sanos caminando a 4 velocidades diferentes, siguiendo el trabajo de Liu [ 10 ]. Las articulaciones subtalar y metatarso-falángica se bloquearon en ángulos anatómicos neutros. El protocolo de toma de datos fue el mismo definido en el trabajo de Schwartz [ 16 ]. En el se sigue un protocolo Plug-in-Gait usando cámaras Vicon. Este protocolo fue desarrollado en la Universidad de Oxford y se basa en el modelo Newington de marcha [ 11 ]. En el se usan de un numero mínimo de marcadores para reconstruir el movimiento del sujeto. Se trata de un modelo predictivo, lo que quiere decir que se utilizan ecuaciones de regresión para obtener las coordenadas de los centros de articulaciones en el sistema de coordenadas local como función de parámetros antropométricos. La principal ventaja de este protocolo es que es mínimamente invasivo [ 11 ]. Se asume que mediante el uso de este protocolo se obtuvieron los orígenes de los sistemas de coordenadas asociados a cada segmento y las coordenadas de los centros de las articulaciones, que constituyen los parámetros cinemáticos del modelo. Parámetros dinámicos Los parámetros dinámicos del modelo son necesarios para la resolución de la dinámica inversa. Los datos necesarios son las masas y momentos de inercia de cada segmento del modelo. Con estos datos se construye la matriz de masa que aparece en la ecuación dinamica. En este trabajo estos parámetros se han tomado del código de OpenSim. Estos parámetros están basados en datos medios antropométricos obtenidos de cinco sujetos para todos los segmentos excepto para calcáneo y dedos. Para estos últimos los datos se obtienen de representar su volumen a partir de la medida de una zapatilla de tenis [2]. 3.2 Modelos simplificados La primera simplificación realizada consiste en considerar el pie como monosegmento. Aunque, como se ha visto, en el modelo original el pie esta formado por tres cuerpos (talón, calcáneo y toes) unidos por dos articulaciones (subtalar y metatarso-falángica) en las opciones por defecto de OpenSim estas articulaciones se fuerzan a cero, por lo que en realidad el pie se modela como un segmento único. Esta simplificación se entenderá mucho mejor tras haber leído el capitulo dedicado a la cinemática inversa. La segunda simplificación se tomo durante el desarrollo del modulo de dinámica inversa debido a la dificultad del tratamiento del spline de la rodilla. De esta manera, se elimino el desplazamiento relativo entre tibia y fémur convirtiéndose la articulación de la rodilla en un par de rotación sencillo. La posición de este centro de rotación se hizo fija y se tomo como la media de las posiciones para cada uno de los valores del spline. En la Figura 3.6 se representa el movimiento que realiza el origen de la tibia en el plano sagital considerando el spline para todos los valores posibles de angulo de flexión de la rodilla comparado con el eje fijo que se toma al eliminar el spline. 4 Cinemática Inversa En el presente capítulo se explica cómo se ha resuelto el problema cinemático inverso. El objetivo de este problema es calcular el valor de las coordenadas generalizadas a partir de las coordenadas de los marcadores del modelo. La obtención del valor de estas coordenadas se ha explicado en el apartado de medidas experimentales y como se utilizan para obtener los parámetros cinemáticos del modelo se ha explicado en el apartado correspondiente del segundo capitulo. Con estos datos de partida ya se puede plantear el problema cinemático inverso. En este capitulo se explica en primer lugar como se formula, para posteriormente describir como se ha implementado en Matlab y finalmente presentar los resultados obtenidos y su comparación con los de OpenSim. 4.1 Formulación De la manera en que está planteado el problema se tiene una mayor numero de ecuaciones que incógnitas. Las ecuaciones son la posición de 31 marcadores (de valor conocido) en función de las coordenadas generalizadas. Por tanto se tienen 93 ecuaciones (3 por cada marcador) frente a 23 incógnitas (el valor de las coordenadas generalizadas). Se esta por tanto frente a un sistema de ecuaciones sobredeterminado. La solución pasa por tanto por encontrar el valor de las coordenadas generalizadas no que cumplan las 93 ecuaciones si no que produzcan el menor error posible. Tenemos por tanto una optimizacion numérica en la que la función a minimizar es el error cometido y se define de la siguiente manera [5]. min"∑ iεmarcadores ωi||xexp i−xi(q)||2+∑ jεcoordenadas prescritas ωj(qexp j−qj)2#(4.1) Donde: •ωies un peso asociado al marcador i •ωjes un peso asociado a la coordenada generalizada j •xexp ies la posición en globales del marcador i registrada en el laboratorio •xi(q) es la posición en globales del marcador i calculada analíticamente a partir de los valores de las coordenadas generalizadas, como se ha explicado en el capítulo dedicado a los conceptos teóricos. •qexp j es el valor de la coordenada generalizada j que el usuario puede forzar a que tenga un valor determinado (coordenada de valor prescrito) •qjes el valor de la coordenada generalizada j obtenida de manera analítica 23 24 Capítulo 4. Cinemática Inversa Se aprecia que este error tiene una componente asociada a la posición de los marcadores y otra asociada a las coordenadas generalizadas. La componente asociada a la posición de los marcadores es la obvia, ya que se pretende que la posición calculada de los marcadores se parezca lo máximo posible a la obtenida de manera experimental. La contribución del error de cada marcador se controla mediante el peso asociado. En este trabajo se han utilizado los vectores de pesos por defecto de OpenSim. En cuanto al vector de pesos ωi este contiene pesos mas altos para los marcadores mas cercanos al hueso y mas bajos para aquellos que pueden moverse debido al movimiento del tejido blando. Esto se debe a que este movimiento de tejido blando desvirtúa la hipótesis tomada de que cada cuerpo es un solido rígido e introduce errores en la reconstrucción del movimiento muy importantes. De esta manera, para los marcadores R.ASIS y L.ASIS, colocados en el punto mas antero-superior de la cresta iliaca (un punto anatómico fácilmente reconocible y muy cercano al hueso de la pelvis), se tienen unos pesos de 10 frente a los marcadores colocados en el muslo, que tienen un peso asociado unidad ya que estarán sometidos al movimiento de todo la musculatura del muslo. Figura 4.1 Marcadores de la pelvis y del muslo. En la Tabla 4.1 se aprecia el peso ωiasociado a cada marcador. Por su parte, la segunda componente del error tan solo sirve si hay alguna coordenada de valor prescrito. En las opciones por defecto de OpenSim las únicas coordenadas de valor prescrito han sido las asociadas al movimiento del pie, es decir, subtalar_angle y mtp_angle. Ambas se fuerzan a cero con un peso asociado de 1000. De esta manera, aunque originalmente el modelo considera el pie formado por tres segmentos, en la practica se tiene un pie monosegemento. En todo momento se sigue el protocolo de reconstrucción del movimiento conocido como LOM (Local Optimization Method). Este método se centra en minimizar el error cometido al obtener la posición de los marcadores en función de las coordenadas generalizadas sin imponer ninguna restricción cinemática sobre ellas. Este método se contrapone al método AM (Adaptation Method), el cual desprecia el error cometido limitándose a imponer sobre las posiciones de los marcadores las restricciones cinemáticas del modelo. [11] En un principio se utilizo la formula presentada en la ecuación (4.1) para calcular el error. Los resultados obtenidos usando esta formulación se presentan como “modelo original Delp” en la sección de resultados. Posteriormente, se elimino la parte del error asociada al valor de las coordenadas generalizadas y se utilizo la restricción que permite fmincon Aeqx=beq para forzar las 4.1 Formulación 25 Tabla 4.1 Peso asociado a cada marcador. Marcador ωi Sternum 1 R.Acronium 0.5 L.Acronium 0.5 TopHead 0.1 R.ASIS 10 L.ASIS 10 V.Sacral 10 R.Thigh.Upper 1 R.Thigh.Front 1 R.Thigh.Rear 1 R.Shank.Upper 1 R.Shank.Front 1 R.Shank.Rear 1 R.Heel 10 R.Midfoot.Sup 1 R.Midfoot.Lat 1 R.Toe.Lat 1 R.Toe.Med 1 R.Toe.Tip 10 L.Thigh.Upper 1 L.Thigh.Front 1 L.Thigh.Rear 1 L.Shank.Upper 1 L.Shank.Front 1 L.Shank.Rear 1 L.Heel 10 L.Midfoot.Sup 1 L.Midfoot.Lat 1 L.Toe.Lat 1 L.Toe.Med 1 L.Toe.Tip 10 coordenadas a cero. La función objetivo queda como sigue: min"∑ iεmarcadores ωi||xexp i−xi(q)||2#(4.2) Donde: •ωies un peso asociado al marcador i •xexp ies la posición en globales del marcador i registrada en el laboratorio •xi(q)es la posición en globales del marcador i calculada analíticamente Los resultados obtenidos usando esta formulación se presentan como “modelo original con bloqueo” en la sección de resultados. Tanto en estos resultados como los siguientes que se presentan se ha seguido esta formulación, ya que se consigue un menor tiempo de ejecución. 26 Capítulo 4. Cinemática Inversa Una segunda mejora en la ejecución del código se desarrollo cuando se descubre que la parte mas costosa computacionalmente es la que corresponde a la construcción de los splines para hallar la posición del origen de la tibia. Así, se construyó un modelo similar al original pero no considerando el desplazamiento relativo entre fémur y tibia, por lo que la rodilla se modela como un par de rotación. Este modelo produjo una reducción del tiempo de ejecución bastante importante y se presenta en la sección de resultados como “modelo sin splines”. La ultima mejora consistió en eliminar del modelo sin splines las coordenadas generalizadas asociadas a los grados de libertad del pie. De esta manera, los segmentos talón, calcáneo y dedos pasan a combinarse en uno solo denominado pie. Es asumible realizar esta simplificación ya que aunque el modelo original del pie es multisegmento, en la práctica los grados de libertad asociados a dichos segmentos se restringen. Con esta simplificación, al eliminarse los dos últimos segmentos de la cadena cinemática, que generan expresiones mas complicadas, la resolución del problema se vuelve más sencilla. Los resultados cinemáticos de este modelo se presentan en la sección de resultados como “modelo monosegmento sin splines”. 4.2 Modelado en Matlab La optimización numérica se ha llevado a cabo usando la función de Matlab fmincon , un algoritmo que busca el mínimo de una función de varias variables no lineal y restringida [ 1 ]. Este algoritmo plantea el problema de minimización de la siguiente forma. m´ ın xf(x)tal que           c(x)≤0 ceq =0 Ax ≤b Aeqx=beq lb≤x≤ub (4.3) En un primer momento, los únicos parámetros que se utilizaron fueron los limites superior e inferior lb y ub con los mismos valores que toma OpenSim por defecto. Posteriormente se uso la restricción Aeqx=beq. En un primer momento se utilizo lenguaje simbólico para definir la formula del error a minimizar en el modulo. El objetivo era calcular el gradiente del error de forma simbólica para pasarlo a la función fmincon, mejorando de esta manera su rendimiento. Sin embargo, aunque se consiguió hallar de forma simbólica la formula del error, la extensión del gradiente hacia muy poco practico su uso. Los resultados actuales surgen de haber usado calculo numérico en MATLAB, después de haber realizado una optimización del código para mejorar el rendimiento, ya que en un primer momento el modulo empleaba un tiempo de unos 30 minutos, un tiempo prohibitivo. La optimización del código se llevo a cabo usando cuando era posible vectorización, al ser el modo en el que esta pensado Matlab para funcionar y minimizando el numero de operaciones a realizar. Por otra parte se ha utilizado una función recurrente siempre que ha sido posible. Esto consiste en aprovechar los cálculos ya realizados y supone una ventaja cuando se hace uso de coordenadas relativas. En cuanto al algoritmo utilizado, el que mejor resultados ha dado ha sido el sqp (sequential quadratic programming) frente al interior point. Los criterios de convergencia han sido TolFun 1e-4 (como OpenSim) y máximo numero de iteraciones 400 (frente a 1000 de OpenSim). A continuación se describe de manera abreviada la estructura del código: • rot_mat , rot_pie , rot_rodilla : estas funciones sirven para calcular las matrices de rotación empleadas en el modelo. La función rot_mat devuelve la matriz de rotación asociada a pelvis, cadera y espalda dados los tres ángulos de Euler en radianes; la matriz rot_pie devuelve la matriz de rotación asociada a tobillo, articulación subtalar o metatarso-falángica en función de la coordenada generalizada asociada; la función rot_rodilla devuelve la 4.3 Presentación de resultados 27 matriz de rotación asociada a la rodilla dado el valor de la coordenada generalizada knee angle. • calcula_globales : esta función devuelve la posición en coordenadas globales de un punto dados su posición en el sistema de coordenadas local así como el origen y la matriz de orientación asociados a este. • cinematica2392 : esta función devuelve las posiciones de los marcadores del modelo dado el vector de coordenadas generalizadas q. En esta función se llama tanto a las funciones rot_mat , rot_pie y rot_rodilla como a la función calcula_globales . • sq_err : esta función calcula el valor de la función objetivo para un instante de tiempo dados el vector de coordenadas generalizadas y las posiciones de los marcadores medidas de forma experimental. Esta función llama a la función cinematica2392 . • ik_main : es el código principal en el que se resuelve el problema cinemático inverso realizando la optimización descrita al principio del capitulo. En esta función se realiza un bucle en el que se llama a fmincon para minimizar la función sq_err para cada instante de tiempo. Una vez ejecutado este código se obtiene el valor de las coordenadas generalizadas para cada instante de tiempo, que se analizan en el siguiente apartado. 4.3 Presentación de resultados 4.3.1 Modelo original Delp De la Figura 4.2 a la Figura 4.5 se representan los valores obtenidos de las coordenadas generalizadas frente al porcentaje del ciclo de marcha usando esta formulación comparados con los de OpenSim. Se observa que los resultados obtenidos son casi idénticos a los de OpenSim pero con un tiempo de ejecución mucho mayor, 73 segundos frente a tan solo 2. En la Figura 4.3 y la Figura 4.4 se observa que se obtienen valores cercanos a cero para las coordenadas subtalar angle y mtp angle, pero con mucho ruido. El error se ha calculado usando la fórmula propuesta por [ 12 ] para hallar el error residual asociado a la posición de cada marcador: ei=1 nt ∑ j ||xexp ij −xij(q)|| (4.4) Donde xexp ij es la posición en el instante j del marcador i registrada en el laboratorio , xi(q) es la posición en el instante j del marcador i calculada analíticamente y nt el número de instantes de tiempo considerados. Posteriormente se calcula una media ponderada de los errores de cada marcador: error =∑iωiei ∑iωi (4.5) Donde ωies un peso asociado al marcador i Tabla 4.2 Comparación de métodos IK. Modelo Error (cm) Tiempo de ejecución (s) OpenSim 0.7797 2 Modelo original Delp 0.7795 73 34 Capítulo 4. Cinemática Inversa Tabla 4.6 Error asociado a cada marcador. Marcador Modelo original Delp Modelo original con bloqueo Modelo sin splines Modelo sin splines monosegmento Sternum 2.04 2.04 2.19 2.19 R.Acronium 7.75 7.75 7.95 7.95 L.Acronium 6.19 6.19 6.36 6.36 TopHead 3.26 3.26 3.23 3.23 R.ASIS 0.58 0.58 0.65 0.65 L.ASIS 0.72 0.72 0.77 0.77 V.Sacral 0.93 0.93 0.69 0.69 R.Thigh.Upper 1.49 1.49 1.47 1.47 R.Thigh.Front 1.80 1.80 1.85 1.85 R.Thigh.Rear 1.31 1.31 1.36 1.36 R.Shank.Upper 0.80 0.80 0.89 0.89 R.Shank.Front 0.95 0.95 1.14 1.14 R.Shank.Rear 0.63 0.63 0.64 0.64 R.Heel 0.46 0.46 0.56 0.56 R.Midfoot.Sup 0.71 0.71 0.58 0.58 R.Midfoot.Lat 0.67 0.67 0.75 0.75 R.Toe.Lat 0.93 0.93 0.95 0.95 R.Toe.Med 0.80 0.80 0.73 0.73 R.Toe.Tip 0.36 0.36 0.34 0.34 L.Thigh.Upper 1.74 1.75 1.73 1.73 L.Thigh.Front 2.18 2.18 2.16 2.16 L.Thigh.Rear 1.69 1.69 1.67 1.67 L.Shank.Upper 0.93 0.93 0.98 0.98 L.Shank.Front 0.94 0.94 0.98 0.98 L.Shank.Rear 0.78 0.78 0.79 0.79 L.Heel 0.58 0.58 0.62 0.62 L.Midfoot.Sup 1.09 1.09 0.99 0.99 L.Midfoot.Lat 0.62 0.62 0.69 0.69 L.Toe.Lat 0.72 0.72 0.78 0.78 L.Toe.Med 0.74 0.74 0.65 0.65 L.Toe.Tip 0.46 0.46 0.45 0.45 5 Dinámica Inversa En el presente capítulo se explica cómo se ha resuelto el problema dinámico inverso. El objetivo de este problema es calcular el valor de las fuerzas y momentos en las articulaciones a partir del valor de las coordenadas generalizadas obtenidas en la resolución del problema cinemático inverso. En este capitulo se explica en primer lugar como se formula, para posteriormente describir como se ha implementado en Matlab y finalmente presentar los resultados obtenidos y su comparación con los de OpenSim. 5.1 Formulación Para este módulo la formulación se reduce a resolver el conjunto de ecuaciones de Newton-Euler que se han descrito en el capítulo dedicado a los conceptos teóricos. Qgdl =M(q)¨ q−Qv(q,˙ q)−Qgrav(q)(5.1) 5.2 Modelado en Matlab A continuación se presenta un extracto del código en el que se aprecia que la velocidad angular de la pelvis se calcula siguiendo la formulación del profesor Escalona, mientras que las velocidades angulares de espalda y cadera se definen directamente ya que su expresión es análoga a la de la pelvis. Código 5.1 Codigo de generacion de matrices simbolicas. %% Derivadas de las matrices de rotacion dA_pelvis = deriva_matriz_optimizado(A_pelvis,dq); w_pelvis_sk = simplify(A_pelvis' * dA_pelvis); w_pelvis = anti_skew (w_pelvis_sk); % aqui empieza lo interesante w_back = [dq18*cos(q19) - dq17*cos(q18)*sin(q19); dq19 + dq17*sin(q18); dq18*sin(q19) + dq17*cos(q18)*cos(q19)]; % misma expresion que w_pelvis pero con q17,q18,q19 35 36 Capítulo 5. Dinámica Inversa w_hip_r = [dq8*cos(q9) - dq7*cos(q8)*sin(q9); dq9 + dq7*sin(q8); dq8*sin(q9) + dq7*cos(q8)*cos(q9)]; % misma expresion que w_pelvis pero con q7,q8,q9 w_hip_l = [dq13*cos(q14) - dq12*cos(q13)*sin(q14); dq14 + dq12*sin(q13); dq13*sin(q14) + dq12*cos(q13)*cos(q14)]; % misma expresion que w_pelvis pero con q12,q13,q14 El proceso de creación de las matrices dinámicas es el siguiente: en primer lugar uso de calculo simbólico para obtener las expresiones de todas las matrices necesarias. En ese mismo código, uso de la función matlabFunction para convertir todas las expresiones en funciones de Matlab. Posteriormente y en otro código, llamada a dichas funciones para resolver el problema de forma numérica. Código 5.2 Codigo de generacion de matrices simbolicas. q = sym('q',[19 1],'real') ; dq = sym('dq', [19 1],'real'); syms q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 q15 q16 q17 q18 q19 real syms dq1 dq2 dq3 dq4 dq5 dq6 dq7 dq8 dq9 dq10 dq11 dq12 dq13 dq14 dq15 dq16 dq17 dq18 dq19 real G_pelvis = jacobian(w_pelvis,dq); matlabFunction(G_pelvis,'File','G_pelvis_maza','Vars', {q,dq}) Código 5.3 Uso de matrices simbolicas. clear all close all clc load('qopt_sin_spline') for i= 1:119 G_pelvis_modelo = G_pelvis_maza (q(:,i),dq(:,i)); T(:,i) = S\Q_accion; end Conviene aclarar que las coordenadas generalizadas se han filtrado utilizando un filtro pasobajo IIR tomado de OpenSim y traducido a Matlab. La frecuencia de corte escogida han sido 6 Hz, como hace OpenSim de manera predeterminada. Sin este filtro las fuerzas generalizadas tienen mucho ruido ya que son proporcionales a la derivada segunda de las coordenadas generalizadas. El tiempo de ejecución del código que escribe las funciones es bastante elevado, debido a que se hace uso de calculo simbólico. Una de las tareas mas costosas computacionalmente es cuando se ejecuta matlabFunction. En los casos de matrices con una expresión mas compleja es necesario 5.3 Presentación de resultados 37 0.5 1 1.5 2 2.5 t (s) 2.5 3 3.5 4 4.5 5 5.5 6 deg pelvis tilt original filtrado Figura 5.1 Filtrado. 0.5 1 1.5 2 2.5 t (s) -30 -20 -10 0 10 20 30 deg/s derivada de pelvis tilt original filtrado Figura 5.2 Filtrado. desactivar la optimización que Matlab hace por defecto al crear la función. Esto se lleva a cabo utilizando la opción Óptimize’, false. Por su parte, el tiempo de ejecución del código que evalúa numéricamente las matrices y obtiene el valor de los momentos para cada instante de tiempo es de unos 10 segundos. 5.3 Presentación de resultados El valor de las coordenadas generalizadas utilizado como dato de partida ha sido el obtenido con el modelo monosegmento sin splines. De la Figura 5.3 a la Figura 5.6 se presentan los resultados obtenidos y su comparación con los de OpenSim. Se aprecia que se obtienen resultados aproximados a los de OpenSim excepto para los momentos en la pelvis y el momento de flexión del tobillo. En ambos casos se desconoce el motivo por el que no se obtienen unos buenos resultados. En la Tabla 5.1 se presenta el error cuadrático medio cometido para cada fuerza generalizada. 38 Capítulo 5. Dinámica Inversa 0 50 100 % ciclo marcha 450 500 550 600 650 700 750 800 850 900 950 N pelvis ty force calculado opensim 0 50 100 % ciclo marcha -150 -100 -50 0 50 100 150 N pelvis tx force calculado opensim 0 50 100 % ciclo marcha -150 -100 -50 0 50 100 N pelvis tz force calculado opensim 0 50 100 % ciclo marcha -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 Nm pelvis tilt moment calculado opensim 0 50 100 % ciclo marcha -40 -30 -20 -10 0 10 20 30 40 50 Nm pelvis list moment calculado opensim 0 50 100 % ciclo marcha -20 -15 -10 -5 0 5 10 15 20 Nm pelvis rotation moment calculado opensim Figura 5.3 Resultados dinámicos para la pelvis. 0 50 100 % ciclo marcha -80 -60 -40 -20 0 20 40 60 Nm hip r flexion moment calculado opensim 0 50 100 % ciclo marcha -15 -10 -5 0 5 10 15 Nm hip r adduction moment calculado opensim 0 50 100 % ciclo marcha -6 -4 -2 0 2 4 6 8 Nm hip r rotation moment calculado opensim 0 50 100 % ciclo marcha -30 -25 -20 -15 -10 -5 0 5 10 15 20 Nm knee r flexion moment calculado opensim 0 50 100 % ciclo marcha -6 -4 -2 0 2 4 6 8 Nm ankle r flexion moment calculado opensim Figura 5.4 Resultados dinámicos para la pierna derecha. 0 50 100 % ciclo marcha -30 -20 -10 0 10 20 30 40 50 Nm hip l flexion moment calculado opensim 0 50 100 % ciclo marcha -15 -10 -5 0 5 10 15 Nm hip l adduction moment calculado opensim 0 50 100 % ciclo marcha -3 -2 -1 0 1 2 3 4 5 6 Nm hip l rotation moment calculado opensim 0 50 100 % ciclo marcha -25 -20 -15 -10 -5 0 5 10 15 Nm knee l flexion moment calculado opensim 0 50 100 % ciclo marcha -1 -0.5 0 0.5 1 1.5 2 2.5 3 Nm ankle l flexion moment calculado opensim Figura 5.5 Resultados dinámicos para la pierna izquierda. 5.3 Presentación de resultados 39 0 50 100 % ciclo marcha -5 0 5 10 15 20 25 30 Nm lumbar extension moment calculado opensim 0 50 100 % ciclo marcha -50 -40 -30 -20 -10 0 10 20 30 40 Nm lumbar bending moment calculado opensim 0 50 100 % ciclo marcha -6 -4 -2 0 2 4 6 Nm lumbar rotation moment calculado opensim Figura 5.6 Resultados dinámicos para la espalda. Tabla 5.1 Error asociado a cada fuerza generalizada. Fuerza generalizada Error (N o Nm) pelvis_ty_force 63.38 pelvis_tx_force 31.33 pelvis_tz_force 26.82 pelvis_tilt_moment 42.11 pelvis_list_moment 15.97 pelvis_rotation_moment 5.55 hip_r_flexion_moment 14.47 hip_r_adduction_moment 3.08 hip_r_rotation_moment 1.52 knee_r_flexion_moment 5.59 ankle_r_flexion_moment 1.96 hip_l_flexion_moment 10.29 hip_l_adduction_moment 3.46 hip_l_rotation_moment 1.28 knee_l_flexion_moment 5.55 ankle_l_flexion_moment 0.36 lumbar_extension_moment 4.06 lumbar_bending_moment 9.44 lumbar_rotation_moment 1.55 6 Conclusiones La s conclusiones generales que se pueden extraer de la realización de este trabajo es que se obtienen unos resultados muy similares a los obtenidos por OpenSim pero con un tiempo de ejecución mucho mayor. Esto se debe a que OpenSim esta programado en C++, un lenguaje de programación compilado. Matlab mas que un lenguaje de programación es una herramienta de calculo, muy potente, pero que no se puede enfrentar a otros lenguajes de programación. También hay que tener en cuenta que el programa de OpenSim ha sido desarrollado por un equipo con ingenieros de software, lo que da una idea de la importancia que tiene un código bien escrito y optimizado para las simulaciones cinemáticas y dinámicas. En cuanto a los diferentes modelos planteados, se considera que el modelo monosegmento sin splines es adecuado para el objetivo propuesto de emular a OpenSim ya que con el se obtienen valores muy similares a los del programa y además permite llevar a cabo el módulo de dinámica inversa en Matlab. En cuanto a los resultados obtenidos en este último módulo para los momentos en la pelvis y el momento de flexión del tobillo no se obtienen valores similares a los de OpenSim por alguna razón desconocida. Una posible mejora de este trabajo pasaría por obtener los residuales dinámicos para evaluar el error cometido en los resultados dinámicos en lugar del error cuadrático medio. Otra pasaría por abarcar no solo el modelo gait2392, si no muchos más y que al comienzo del programa se le preguntara al usuario que modelo quiere usar y con qué opciones. Esto permitiría que usuarios sin muchos conocimientos de programación pudiesen usar el programa y a la vez, al estar programado en Matlab, un lenguaje mas sencillo que C++, que usuarios con mayores conocimientos de programación pudiesen modificar el código. Otra posible mejora de este trabajo pasaría por hacer uso de Python, un lenguaje de programación orientado a objetos que dispone de librerías muy interesantes de álgebra matricial. Los principales retos a los que se ha tenido que hacer frente en el desarrollo de este trabajo han sido en primer lugar la necesaria comprensión del modelo gait2392 para programar su cinemática. Posteriormente, en el modulo de dinámica inversa Matlab se quedaba sin memoria al hallar las matrices dinámicas utilizando los procedimientos descritos en el documento del profesor JL Escalona. Este problema quedo solventado al calcular las velocidades angulares de manera directa como se explica en el capitulo dedicado a la dinámica inversa. La continuación natural de este trabajo pasa en primer lugar por la programación en Matlab del modulo de Optimización estática. Este modulo consiste en, a partir de los valores de los momentos 41 42 Capítulo 6. Conclusiones en las articulaciones obtener el valor de las fuerzas producidas por cada musculo en cada instante de tiempo. Estas fuerzas se obtienen de una minimización de la suma de las activaciones de los músculos al cuadrado. El hecho de que sea necesaria una optimización se debe a que por cada momento existe mas de un musculo que contribuye a realizarlo, por lo tanto es necesario un criterio para de alguna manera repartir que musculo hace mas o menos fuerza. En este caso el criterio es la minimización de las activaciones musculares, lo cual es una manera de minimizar la potencia muscular necesaria para realizar el movimiento. El otro problema de este modulo se produce al derivar numéricamente el valor de las coordenadas generalizadas para obtener velocidades y aceleraciones, ya que, como es conocido, la derivación hace se produzca ruido en las señales de velocidad y aceleración. Bibliografía [1] fmincon doc page,https:// es.mathworks.com/ help/ optim/ ug/ fmincon.html. [2] Gait 2392 and 2354 models,https:// simtk-confluence.stanford.edu/ display/ OpenSim/ Gait+ 2392+y+2354+Models. [3] B. Bresler and J.P. Frankel, The forces and moments in the leg during level walking, Transactions of the American Society of Mechanical Engineers 72 (1950), 27–36. [4] M. Damsgaard, J. Rasmussen, S.T. Christensen, E. Surma, and M. de Zee, Analysis of musculoskeletal systems in the AnyBody Modeling System, Simulation Modelling Practice and Theory 14 (2006), no. 8, 1100–1111. [5] S.L. Delp, F.C. Anderson, A.S. Arnold, P. Loan, A. Habib, C.T. John, E. Guendelman, and D.G. Thelen, Opensim: Open-source software to create and analyze dynamic simulations of movement, IEEE Transactions on Biomedical Engineering 54 (2007), no. 11, 1940–1950. [6] S.L. Delp, J.P. Loan, M.G. Hoy, F.E. Zajac, E.L. Topp, and J.M. Rosen, An interactive graphicsbased model of the lower extremity to study orthopaedic surgical procedures, IEEE Transactions on Biomedical Engineering 37 (1990), no. 8, 757–767. [7] J.L. Escalona, Dinámica inversa de una bicicleta a partir de datos experimentales, Apuntes de clase. [8] J.R. Gage, M.H. Schwartz, and S.E. Koop, Identification and treatment of gait problems in cerebral palsy, ch. 2.6, 2009. [9] C.T. John, A. Seth, M.H. Schwartz, and S.L. Delp, Contributions of muscles to mediolateral ground reaction force over a range of walking speeds, Journal of Biomechanics 45 (2012), no. 14, 2438–2443. [10] M.Q. Liu, F.C. Anderson, M.H. Schwartz, and S.L. Delp, Muscle contributions to support and progression over a range of walking speeds, Journal of Biomechanics 41 (2008), no. 15, 3243–3252. [11] J. Ojeda, Application of multibody system techniques to human locomotor system, Ph.D. thesis, Escuela Superior de Ingenieros Industriales of the University of Seville, 2012. [12] J. Ojeda, J. Martínez-Reina, and J. Mayo, A method to evaluate human skeletal models using marker residuals and global optimization, Mechanism and Machine Theory 73 (2014), 259– 272, cited By 10. 43