Programación de un código en MATLAB para la medida de las propiedades mecánicas del tríceps sural mediante la técnica de vibración libre

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Autor: Ignacio López Alonso Tutores: Alberto Barroso Caro Federico París García Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Programación de un código en MATLAB para la medida de las propiedades mecánicas del tríceps sural mediante la técnica de vibración libre. Ignacio López Alonso Trabajo Fin de Grado 1 Ignacio López Alonso Trabajo Fin de Grado 2 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Programación de un código en MATLAB para la medida de las propiedades mecánicas del tríceps sural mediante la técnica de vibración libre. Autor: Ignacio López Alonso Tutores: Alberto Barroso Caro Federico París García Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Ignacio López Alonso Trabajo Fin de Grado 3 (París-García et al., 2015) Ignacio López Alonso Trabajo Fin de Grado 4 Agradecimientos A mis tutores, Alberto Barroso y Federico París García, por su verdadera implicación con el proyecto, el buen trato que me han dado y el ambiente de trabajo favorable que he percibido en todo momento. A mis compañeros y amigos, por su colaboración desinteresada cuando ha sido necesaria y su apoyo. A mi familia, por su constante interés, apoyo y confianza en mí durante este trabajo y en estos cuatro años. A mi pareja Nani, que no ha dudado en apoyarme incondicionalmente y hacerme el camino más llevadero en los momentos más difíciles. Gracias a todos. Ignacio López Alonso Trabajo Fin de Grado 5 Ignacio López Alonso Trabajo Fin de Grado 6 Índice 1. INTRODUCCIÓN ......................................................................................................... 10 1.1 MOTIVACIÓN INICIAL ......................................................................................................... 10 1.2 OBJETIVO ........................................................................................................................ 10 2. DESCRIPCIÓN DE LA TÉCNICA, EL EQUIPO Y EL MÉTODO ............................................ 12 2.1 INTRODUCCIÓN ................................................................................................................ 12 2.2 MODELO DEL TS COMO UN SISTEMA DE 1 GRADO DE LIBERTAD .................................................. 12 2.3 MODELO DE HILL .............................................................................................................. 13 2.4 DESCRIPCIÓN GENERAL DEL MÉTODO Y DE LAS ECUACIONES ....................................................... 14 2.5 RELACIÓN ENTRE PROPIEDADES APARENTES DEL MTC Y LAS DE SUS COMPONENTES ....................... 21 2.6 DESCRIPCIÓN DETALLADA DEL EQUIPO ................................................................................... 23 2.7 DESCRIPCIÓN DEL MÉTODO ................................................................................................. 28 3. TRATAMIENTO DE LOS DATOS EN MATLAB ............................................................... 34 3.1 INTRODUCCIÓN ................................................................................................................ 34 3.2 DESCRIPCIÓN TEÓRICA DEL TRATAMIENTO DE LOS DATOS .......................................................... 34 3.2.1 Primer ajuste por mínimos cuadrados. Ajuste de cada medición ........................... 34 3.2.2 Segundo ajuste por mínimos cuadrados. Obtención de las rigideces del Sóleo y del tendón de Aquiles .............................................................................................................................. 37 3.2.3 Breve descripción del programa Flexor ................................................................... 38 3.3 DESCRIPCIÓN NUMÉRICA DEL TRATAMIENTO DE LOS DATOS. PROGRAMA DE MATLAB ..................... 39 3.3.1 Función de lectura .................................................................................................. 40 3.3.2 Función de filtro ...................................................................................................... 41 3.3.3 Función de primer ajuste por mínimos cuadrados. Ajuste de cada medición ......... 42 3.3.4 Función de segundo ajuste por mínimos cuadrados. Obtención de las rigideces del Sóleo y del tendón de Aquiles ............................................................................................................ 47 4. VERIFICACIÓN Y VISUALIZACIÓN DE RESULTADOS. CONCLUSIONES ........................... 50 4.1 INTRODUCCIÓN ................................................................................................................ 50 4.2 VERIFICACIÓN DE LOS AJUSTES ............................................................................................. 50 4.3 VISUALIZACIÓN DE LOS RESULTADOS ..................................................................................... 52 4.3.1 Ajuste de cada medición ......................................................................................... 52 4.3.2 Resultados del estudio ............................................................................................ 53 5. CONCLUSIONES ......................................................................................................... 56 REFERENCIAS ........................................................................................................................ 58 ANEXOS ................................................................................................................................ 59 ANEXO I. FUNCIÓN DE LECTURA ..................................................................................................... 59 ANEXO II. FUNCIÓN DE FILTRO ....................................................................................................... 61 ANEXO III. FUNCIÓN DE PRIMER AJUSTE POR MÍNIMOS CUADRADOS ...................................................... 62 ANEXO IV. FUNCIÓN DE SEGUNDO AJUSTE POR MÍNIMOS CUADRADOS .................................................. 69 Ignacio López Alonso Trabajo Fin de Grado 7 Ignacio López Alonso Trabajo Fin de Grado 8 Índice de figuras FIGURA 2.1. VARIANTES DEL MODELO DE HILL. (PARÍS-GARCÍA ET AL., 2010) ............................................................ 13 FIGURA 2.2. ESQUEMA GENERAL DEL EQUIPO UTILIZADO. (PARÍS-GARCÍA ET AL., 2010) .............................................. 14 FIGURA 2.3. POSICIÓN QUE EL SUJETO DEBE ADOPTAR DURANTE EL ENSAYO. (PARÍS-GARCÍA ET AL., 2013) ..................... 15 FIGURA 2.4. ESQUEMA SIMPLIFICADO DEL MODELO DEL PIE Y LA ANTEPIERNA, B) EQUILIBRIO DE FUERZAS SOBRE EL PIE, C) POSICIÓN CUALQUIERA DEL PIE DURANTE LA VIBRACIÓN. (PARÍS-GARCÍA ET AL., 2010) ........................................ 16 FIGURA 2.5. A) SISTEMA EQUIVALENTE DE 1 GDL DEL MTC CON VALORES APARENTES DE RIGIDEZ K Y AMORTIGUAMIENTO C. B) EQUILIBRIO DE FUERZAS SOBRE EL SISTEMA. (PARÍS-GARCÍA ET AL., 2010) .................................................... 18 FIGURA 2.6. EQUILIBRIO SOBRE EL SISTEMA DINÁMICO EQUIVALENTE DEL MTC. (PARÍS-GARCÍA ET AL., 2010) ................ 21 FIGURA 2.7. EQUIVALENCIA ENTRE LOS DOS SISTEMAS. A) SISTEMA APARENTE. B) SISTEMA FORMADO POR SÓLEO Y TENDÓN DE AQUILES. (PARÍS-GARCÍA ET AL., 2010) ..................................................................................... 22 FIGURA 2.8. VISIÓN GENERAL DEL EQUIPO DISEÑADO PARA LA PRÁCTICA DE LA TÉCNICA DE VIBRACIÓN LIBRE APLICADA AL TRÍCEPS SURAL. ...................................................................................................................................... 24 FIGURA 2.9. A) CÉLULA DE CARGA OMEGADYNE LCM501. B) SISTEMA DE ESPEJOS CON DIANA QUE PERMITEN APOYAR EL PIE EN EL PUNTO ADECUADO. ......................................................................................................................... 24 FIGURA 2.10. A) MASA CONCENTRADA QUE ACENTÚA LA VIBRACIÓN. B) SISTEMA DE CORREDERA PARA UN AJUSTE LONGITUDINAL DE LA MASA. ..................................................................................................................... 25 FIGURA 2.11. REGULADOR DE ALTURA DE LA SILLA CON GATO. ................................................................................. 26 FIGURA 2.12. REGULADOR DE DISTANCIA ENTRE LA SILLA Y LA CÉLULA DE CARGA. ........................................................ 26 FIGURA 2.13. ASIENTO CON RESPALDO. .............................................................................................................. 27 FIGURA 2.14. VIGA ARTICULADA QUE SOSTIENE LA MASA CONCENTRADA. .................................................................. 27 FIGURA 2.15. POSICIÓN QUE EL SUJETO DEBE ADOPTAR DURANTE EL ENSAYO. (PARÍS-GARCÍA ET AL., 2013) ................... 28 FIGURA 2.16. EJE LONGITUDINAL DEL PIE DE UNA HUELLA EN DINÁMICO. (PARÍS-GARCÍA ET AL., 2010) .......................... 29 FIGURA 2.17. A) PEDÍGRAFO, RODILLO Y BOTE DE TINTA. B) REGLA DE PERTHES. ......................................................... 30 FIGURA 2.18. PUNTOS DE REFERENCIA PARA LA PEDIGRAFÍA. (PARÍS-GARCÍA ET AL., 2010) .......................................... 31 FIGURA 2.19. LÍNEAS DE REFERENCIA EN LA PEDIGRAFÍA PARA DETERMINAR A, B Y C. (PARÍS-GARCÍA ET AL., 2010) ........... 32 FIGURA 3.1. EJEMPLO DE SEÑAL QUE SE OBTIENE DE UNA MEDICIÓN. (PARÍS-GARCÍA ET AL., 2010) ............................... 35 FIGURA 3.2. AMORTIGUAMIENTO DE UNA FUNCIÓN PERIÓDICA SEGÚN LA EXPONENCIAL NEGATIVA. (PARÍS-GARCÍA ET AL., 2010) .................................................................................................................................................. 36 FIGURA 3.3. REPRESENTACIÓN GRÁFICA DE LA RIGIDEZ DEL TRÍCEPS SURAL FRENTE A LA FUERZA QUE LO SOLICITA. (PARÍSGARCÍA ET AL., 2010) ............................................................................................................................. 37 FIGURA 3.4. EJEMPLO DE AJUSTE POR MÍNIMOS CUADRADOS DE LA FUNCIÓN DE RIGIDEZ. (PARÍS-GARCÍA ET AL., 2015) .... 38 FIGURA 3.5. APARIENCIA VISUAL DEL ANTIGUO PROGRAMA FLEXOR. ......................................................................... 39 FIGURA 3.6. EJEMPLO DE APLICACIÓN DE LA FUNCIÓN DE FILTRO. A) ASPECTO DE LA GRÁFICA SIN FILTRAR. B) ASPECTO TRAS APLICAR EL FILTRO. .................................................................................................................................. 42 FIGURA 3.7. IRREGULARIDAD EN UN MÁXIMO DE UNA ONDA QUE PUEDE DAR LUGAR A LA DETECCIÓN DE FALSOS MÁXIMOS. 44 FIGURA 3.8. RANGO TEMPORAL INICIAL DE AJUSTE. ............................................................................................... 45 FIGURA 3.9. REPRESENTACIÓN GRÁFICA DE LA RIGIDEZ DEL TRÍCEPS SURAL FRENTE A LA FUERZA QUE LO SOLICITA. ............. 48 FIGURA 4.1. COMPARATIVA DEL AJUSTE SOBRE UNA MEDICIÓN EN DISTINTAS VENTANAS TEMPORALES. A) BUEN AJUSTE EN VENTANA TEMPORAL COMPRENDIDA ENTRE EL PRIMER Y EL SEGUNDO MÁXIMO. B) MAL AJUSTE ENTRE EL SEGUNDO Y TERCER MÁXIMO. ................................................................................................................................... 52 FIGURA 4.2. FORMATO DE REPRESENTACIÓN GRÁFICA DE LAS MEDICIONES Y SUS AJUSTES. ............................................ 53 FIGURA 4.4. REPRESENTACIÓN GRÁFICA DE LAS PROPIEDADES VISCOSO-ELÁSTICAS DEL SÓLEO Y EL TENDÓN DE AQUILES. .... 54 FIGURA 4.3. TABLA RESUMEN QUE CONTIENE LOS VALORES DE F, K Y R2 PARA CADA MEDICIÓN. ..................................... 54 FIGURA 5.1. EVALUACIÓN DE LA FUNCIÓN ERROR EN FUNCIÓN DE LA FRECUENCIA. ....................................................... 57 Ignacio López Alonso Trabajo Fin de Grado 15 1. Célula de carga: recoge la fuerza de reacción que el pie efectúa sobre el punto de apoyo. 2. Masa concentrada variable ( 𝑀! ) y punto de impacto: masa que acentúa el movimiento longitudinal de la antepierna. El tendón de Aquiles presenta una rigidez constante con la fuerza que lo solicita, pero la rigidez del sóleo es proporcional a la dicha fuerza. Esta es la razón por la que la masa 𝑀! debe ser variable, ya que aplicar distintas cargas es la única manera de poder medir ambas rigideces. Cabe destacar que la masa equivalente que está oscilando, a la que llamaremos 𝑀 , es igual a la suma de la masa concentrada 𝑀! más la masa distribuida de la pierna, y su valor es difícil de medir. La masa concentrada cuenta con un apoyo semicilíndrico bajo la misma que transmite la carga sobre la pierna puntualmente en el plano xy, y de forma distribuida en la dirección z. 3. y 4. Ajustes que permiten regular la altura de la silla y la distancia de la silla a la célula de carga. Estos ajustes son imprescindibles para asegurar una correcta posición del sujeto sobre el que se realiza el ensayo. 5. Silla: la silla cuenta con un respaldo sobre el que el sujeto debe apoyar en todo momento la espalda y los hombros para evitar variaciones de presión sobre la célula que podrían estar provocadas simplemente por movimientos del cuerpo durante el ensayo. Debemos destacar que el desplazamiento vertical que experimenta la antepierna es del orden de milímetros, por lo que cualquier mínima interferencia o perturbación exterior puede afectar a la calidad de la señal recogida o a la precisión del ensayo. El desplazamiento vertical de la antepierna implica una rotación del muslo alrededor de la cadera, la cual se puede suponer del orden de milirradianes. Se podría pensar que el contacto del muslo con la base de la silla durante la rotación podría suponer una de estas perturbaciones. Sin embargo, debido al orden de esta rotación, y a la pequeña profundidad con la que cuenta Figura 2.3. Posición que el sujeto debe adoptar durante el ensayo. (París-García et al., 2013) Ignacio López Alonso Trabajo Fin de Grado 16 el asiento (31 cm.), esto no supone apenas interferencia, por lo que se p suponer que la antepierna oscila libremente. 6. Viga articulada: la masa suspendida sobre la rodilla 𝑀! está sujeta por una viga articulada al bastidor (7) que permite la transmisión del movimiento oscilatorio vertical a la rodilla. 7. Bastidor: es el soporte para el resto de elementos del equipo y presenta una triangularización que aporta suficiente rigidez al conjunto para garantizar la vibración libre de la pierna y una precisa recogida de la señal. La disposición descrita es válida para el ensayo de ambas piernas sin la necesidad de modificar la posición ni la interpretación de los datos. Se desarrollan a continuación las ecuaciones que rigen el ensayo, y que permiten hallar los parámetros viscoso-elásticos del TS. Cuando se introduce un impacto sobre la masa concentrada 𝑀! , la masa equivalente 𝑀 comienza a oscilar. La Figura 2.4 (a) muestra una posición intermedia arbitraria de dicha oscilación en la que se distinguen varios parámetros. Entre ellos, las distancias 𝐿 (longitud de la antepierna), 𝑅 (longitud del antepié desde la articulación del tobillo) y 𝑟 (longitud del retropié desde la articulación del tobillo) permanecen constantes a lo largo de toda la oscilación, mientras que la longitud del TS varía a lo largo de la oscilación. Figura 2.4. Esquema simplificado del modelo del pie y la antepierna, b) equilibrio de fuerzas sobre el pie, c) posición cualquiera del pie durante la vibración. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 17 Se observa que todo el conjunto se sostiene y oscila sobre un punto que se corresponde con el punto de contacto entre el antepié y la célula de carga, y sobre el que se mide la fuerza de reacción en la célula de carga. El pie debe apoyar a la altura de la cabeza del segundo metatarsiano. La masa 𝑀 es, como se ha explicado anteriormente, la suma de la masa concentrada 𝑀! y la masa distribuida de la pierna, y representa la masa total equivalente del sistema dinámico de 1 grado de libertad. El TS se ha modelado como un sistema dinámico de 1 grado de libertad formado por un amortiguador de constante 𝑐 y un resorte elástico de constante 𝑘 en paralelo, de acuerdo con el modelo de Hill (Hill, 1938). Sin embargo, de acuerdo con el estudio realizado por Fukashiro et al. (Fukashiro et al., 2001), se ha realizado una pequeña modificación del modelo original de Hill (Hill, 1938), obviando el elemento elástico en serie, ya que su contribución al movimiento oscilatorio provocado por el impacto es despreciable. Se puede comprobar que, en la posición en la que se sitúa el sujeto a ensayar, manteniendo los ángulos rectos descritos, el Gastrocnemio no trabaja, por lo que los valores aparentes de 𝑐 y 𝑘 del TS se corresponden con los valores de rigidez y viscosidad del Sóleo y del tendón de Aquiles. Aunque en la Figura 2.4 (b) se ha representado el TS inclinado un cierto ángulo respecto a la dirección longitudinal de la antepierna, durante el desarrollo teórico que se describe a continuación se ha considerado que el TS se alarga y se contrae en una dirección paralela a dicha dirección longitudinal. Los desplazamientos 𝑈 y 𝑢 que aparecen en la Figura 2.4 (c) se relacionan entre sí, aplicando la regla de la palanca, como se indica en la Ecuación (2.1): 𝑈=𝑅 𝑟𝑢 (2.1) Se procede a continuación a aplicar al conjunto descrito en la Figura 2.4 (a) equilibrio de momentos alrededor de la articulación del tobillo. En la Figura 2.4 (b) se indican las fuerzas que intervienen en el equilibrio. 𝐹"# es la fuerza de compresión que actúa sobre el conjunto tibia peroné, mientras que 𝑓 es la fuerza de tracción que solicita al TS. Aplicando el equilibrio, se obtiene la siguiente expresión: ,𝐼+𝑚#𝑟$ %0𝑑%𝜃 𝑑𝑡%=𝐹𝑅𝑐𝑜𝑠𝜃−𝑓𝑟𝑐𝑜𝑠𝜃 (2.2) Donde 𝜃 es el ángulo que define la posición instantánea del pie, medido desde la articulación del tobillo; 𝐼 es el momento de inercia del pie respecto de su centro de gravedad; 𝑚# es la masa del pie que interviene en el movimiento; y 𝑟$ es la distancia entre el centro de gravedad del pie y el centro de rotación del conjunto (la articulación del tobillo). Cabe destacar que en la Ignacio López Alonso Trabajo Fin de Grado 18 Ecuación (2.2) no aparecen los términos relacionados con el peso de la masa concentrada 𝑀! ni la fuerza de compresión 𝐹"# que soporta el conjunto tibia-peroné, ya que ambas actúan sobre la articulación del tobillo, por lo que el momento que general alrededor de este punto es nulo. A priori, la determinación de los parámetros que aparecen en el primer miembro de la Ecuación (2.2), 𝐼 , 𝑚# y 𝑟$7 parece compleja, aunque pueden determinarse (Zatsiorsky & Seluyanov, 1982). Sin embargo, puede comprobarse que estos términos son despreciables frente a los del segundo miembro, por lo que dicha ecuación de puede transformar en: 𝐹𝑅 =𝑓𝑟 (2.3) Ecuación que permite calcular la fuerza 𝑓 que solicita al Tríceps Sural en función de 𝐹 (fuerza de reacción que mide la célula de carga), y los brazos de palanca 𝑅 y 𝑟 . El sistema que muestra la Figura 2.4 es equivalente al sistema simplificado que aparece en la Figura 2.5, donde 𝐾 y 𝐶 representan la rigidez y el amortiguamiento aparentes del sistema, y la masa 𝑀 es la suma de la masa concentrada 𝑀! y de la masa distribuida de la pierna. Como ya se ha comentado, esta masa es desconocida, y se determinará en el proceso de ajuste de la respuesta experimental. Aplicando el equilibrio a la masa 𝑀 se obtiene: 𝑀𝑈:=−𝐾𝑈−𝐶𝑈;(=𝐹) (2.4) Figura 2.5. a) Sistema equivalente de 1 gdl del MTC con valores aparentes de rigidez K y amortiguamiento C. b) Equilibrio de fuerzas sobre el sistema. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 19 Que se puede reordenar para obtener: 𝑀𝑈:+𝐶𝑈;+𝐾𝑈 =0 (2.5) Que es la ecuación de un sistema dinámico de un grado de libertad sometido a vibración libre. Para un sistema subamortiguado, la solución adopta la siguiente forma: 𝑈(𝑡)=𝑒&'("(𝐴𝑠𝑒𝑛𝜔)𝑡+𝐵𝑐𝑜𝑠𝜔)𝑡) (2.6) donde: 𝜔)=𝜔D1−𝜉% (2.7) Es habitual emplear la siguiente nomenclatura: 𝛾 =𝜉𝜔 =𝐶 2𝑀𝜔𝜔 = 𝐶 2𝑀 (2.8) luego: 𝜔) %=𝜔%(1−𝜉%)=𝜔%−𝛾% (2.9) Ya podemos, por tanto, obtener una expresión para el valor de 𝐹 . Aplicando (2.4): 𝐹 =(𝑀𝑈:)=𝑒&*"(𝐴+𝑠𝑒𝑛𝜔)𝑡+𝐵+𝑐𝑜𝑠𝜔)𝑡)77 (2.10) Sin embargo, esta no es la fuerza que registrará la célula de carga. El dispositivo registrará una fuerza que llamaremos 𝐹, y que es la suma de la fuerza 𝐹 descrita en la Ecuación (2.10) y la fuerza gravitatoria: 𝐹,=𝐹+𝑀𝑔=𝑒&*"(𝐴+𝑠𝑒𝑛𝜔)𝑡+𝐵+𝑐𝑜𝑠𝜔)𝑡)+𝑀𝑔 (2.11) Cabe destacar que esta fuerza 𝐹, representa siempre una fuerza de compresión sobre el dispositivo, lo que se consigue gracias a la acción de la masa concentrada 𝑀! , garantizando siempre el contacto entre el pie y la célula de carga. Sin embargo, 𝐹 es una fuerza que, medida desde la posición de equilibrio tras aplicar la masa concentrada 𝑀! , tomará valores positivos y negativos en función de la oscilación. Se puede interpretar, por tanto, la fuerza 𝐹, , como la suma Ignacio López Alonso Trabajo Fin de Grado 20 de una función constante de valor el peso de la masa 𝑀 y una oscilación subamortiguada debida a la vibración libre. El valor de 𝐹, se registrará experimentalmente y se ajustará por mínimos cuadrados a la Ecuación (2.11), en la que los parámetros 𝐴+ , 𝐵+ , 𝛾 , 7𝑀 y 𝜔) son desconocidos. Una vez determinados mediante el ajuste los 5 parámetros desconocidos, se pueden obtener los valores de rigidez 𝐾 y amortiguamiento 𝐶 a través de las ecuaciones del movimiento. Así, teniendo en cuenta que: 𝐾 =𝜔%𝑀 (2.12) Se obtiene el valor de 𝐾 : 𝐾 =𝑀(𝜔) %+𝛾%) (2.13) Y de la Ecuación (2.8) se obtiene el valor del coeficiente de amortiguamiento: 𝐶 =2𝛾𝑀 (2.14) Sin embargo, estos valores de rigidez y amortiguamiento no son los correspondientes al TS, sino al sistema equivalente que se ha descrito en la Figura 2.5. Para obtener los valores aparentes del TS 𝑘 y 𝑐 , se debe utilizar la relación (2.3) entre la fuerza que mide la célula de carga 𝐹 y la fuerza que solicita al TS 𝑓 . El equilibrio sobre el MTC representado según el modelo de Hill se muestra en la Figura 2.6, e implica: 𝑓 =−𝑐𝑢; −𝑘𝑢 (2.15) Si aplicamos a (2.4) las relaciones (2.1) y (2.3) se obtiene: J𝑟 𝑅K𝑓 =−𝐶L𝑅 𝑟M%𝑈;−𝐾L𝑅 𝑟M%𝑈 (2.16) Comparando las Ecuaciones (2.15) y (2.16) se deducen los valores aparentes de rigidez y viscosidad del Tríceps Sural en función de los valores del sistema 𝐾 y 𝐶 : Ignacio López Alonso Trabajo Fin de Grado 21 𝑘 =L𝑅 𝑟M%𝐾 (2.17) 𝑐 =L𝑅 𝑟M%𝐶 (2.18) 2.5 Relación entre propiedades aparentes del MTC y las de sus componentes Para obtener los valores de rigidez y viscosidad 𝑘 y 𝑐 del MTC en el apartado anterior, se ha tratado el MTC como un sistema formado por un muelle y un amortiguador en paralelo, fruto de una adaptación del modelo de Hill. Por tanto, se han obtenido las propiedades del Tríceps Sural como conjunto, lo cual puede ser muy útil para determinar la condición física de un sujeto o evaluar el estado del MTC de cara a recuperaciones musculares. Sin embargo, podría ser de gran utilidad conocer los valores de rigidez y viscosidad de los diferentes elementos que forman el Tríceps Sural (Sóleo, Gastrocnemio y Tendón de Aquiles) por separado. Esta es la cuestión que se aborda en este apartado. En lo que respecta a la viscosidad, el tendón de Aquiles no aporta ningún efecto a esta característica. Además, como ya se ha explicado, el Gastrocnemio no trabaja en la posición en la Figura 2.6. Equilibrio sobre el sistema dinámico equivalente del MTC. (París- García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 22 que se ensaya al sujeto. Por tanto, en esta posición, todo el valor de la viscosidad aparente se corresponde con la viscosidad del Sóleo. En lo que respecta a la rigidez, el muelle de rigidez aparente 𝑘 representa la respuesta del Sóleo y del Tendón de Aquiles, ya que, como se ha comentado, el Gastrocnemio no trabaja. Se trata, por tanto, de relacionar este valor de rigidez con los valores individuales del tendón y del Sóleo. Se sabe que el tendón muestra un valor de rigidez constante ( 𝑘" ), que no varía con la fuerza que lo solicita, mientras que el Sóleo se caracteriza por una rigidez ( 𝑘, ) que es función lineal de la carga aplicada 𝑓 : 𝑘,=𝑘-- ⋅𝑓 (2.19) donde 𝑘-- representa la pendiente de la variación de la rigidez del Sóleo con la carga. Para determinar la relación entre 𝑘 , 𝑘" y 𝑘, se han representado en la Figura los dos sistemas equivalentes: el aparente y el formado por el Sóleo y el tendón, cuyos efectos se han puesto en serie. Para ser equivalentes, ambos sistemas deben experimentar el mismo desplazamiento 𝑢 cuando soportan la misma carga 𝑓 . La relación entre 𝑢 y 𝑓 para el caso (a) es: 𝑢 =𝑓 𝑘 (2.20) Figura 2.7. Equivalencia entre los dos sistemas. a) Sistema aparente. b) Sistema formado por Sóleo y tendón de Aquiles. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 23 Mientras que para el caso (b) la relación queda: 𝑢 = 𝑓 𝑘"+𝑓 𝑘, (2.21) Igualando ambas expresiones se obtiene: 𝑓 𝑘=𝑓 𝑘"+𝑓 𝑘, (2.22) Que conduce a: 𝑘 = 𝑘"𝑘, 𝑘"+𝑘, (2.23) Aplicando la relación (2.19): 𝑘 = 𝑘"𝑘--𝑓 𝑘"+𝑘--𝑓 (2.24) Donde 𝑘 y 𝑓 son datos y 𝑘-- y 𝑘" son los parámetros a determinar. Esto se resolverá a través de un ajuste por mínimos cuadrados que se describe en el Capítulo 3. Nótese la diferencia física de los parámetros 𝑘-- y 𝑘" . Mientras que 𝑘" tiene unidades típicas de rigidez (fuerza/longitud, e.g. 𝑘𝑁/𝑚 ), 𝑘-- representa el cambio en la rigidez del Sóleo respecto a la fuerza que lo solicita, por lo que tiene unidades de rigidez/fuerza (i.e. (𝑘𝑁 𝑚) 𝑘𝑁 ⁄⁄ o 𝑚&. ). 2.6 Descripción detallada del equipo En este apartado se realiza una descripción detallada del equipo desarrollado por París García, y se describen algunos aspectos del equipo que motivaron el diseño de un nuevo equipo, ya que dificultan la realización de los ensayos de manera rápida y cómoda. La Figura 2.8 muestra una imagen general del equipo. Ignacio López Alonso Trabajo Fin de Grado 24 A modo de resumen, como se explica en el Apartado 2.4, el dispositivo requiere que el sujeto se halle sentado, manteniendo los ángulos rectos que se indican en la Figura 2.2. La vibración libre de la antepierna la provoca un impacto sobre la rodilla del sujeto, y se acentúa gracias a una carga variable situada sobre la rodilla. El apoyo del pie a la altura de la cabeza del segundo metatarsiano sobre el dispositivo de medida de fuerza (célula de carga) hace que la fuerza 𝑓 que solicita al MTC sea la fuerza 𝐹 que recoge el dispositivo multiplicada por los correspondientes brazos de palanca de antepié y retropié (2.3). El equipo cuenta con una serie de ajustes que hacen que sea adaptable a las medidas antropométricas de diferentes sujetos. Se realiza a continuación una descripción pormenorizada de las diferentes partes que componen el equipo: Figura 2.8. Visión general del equipo diseñado para la práctica de la técnica de vibración libre aplicada al Tríceps Sural. Figura 2.9. a) Célula de carga OmegaDyne LCM501. b) Sistema de espejos con diana que permiten apoyar el pie en el punto adecuado. Ignacio López Alonso Trabajo Fin de Grado 31 2.1) Obtención de los puntos de referencia: 2.1.1. Se impregna de tinta uniformemente la cara interna del elemento elástico. Para ello se verterán varias gotas de tinta sobre el mismo y se extenderá con la ayuda de un rodillo. Se debe evitar un exceso de tinta. 2.1.2. El pedígrafo debe situarse sobre una superficie totalmente lisa y horizontal para el posterior apoyo del pie. 2.1.3. El sujeto debe colocarse en posición anatómica neutra con la espalda apoyada sobre una pared y el pedígrafo también apoyado sobre esa pared. 2.1.4. El apoyo del pie debe realizarse descalzo con toda la planta a la vez, con cuidado de no mover ni el papel ni el pedígrafo. 2.1.5. El pie contrario debe encontrarse descalzo y apoyado en el suelo a la misma altura que el otro pie, con una separación compatible con un apoyo cómodo y uniforme. Se deberá mantener en todo momento el cuerpo erguido, con el rostro mirando al frente, los hombros relajados, brazos estirados y pegados al cuerpo y rodillas en máxima extensión. 2.1.6. Se localizan los relieves anatómicos que se corresponden con los vértices de los maléolos Tibial y Peroneal y, mediante la regla de Perthes, son trasladados a su posición en la pedigrafía (puntos A y B, Figura 2.18). 2.1.7. Localizamos la cabeza del primer metatarsiano mediante la técnica de palpación, y señalamos con el marcador el punto más proximal y distal (puntos C y D). 2.1.8. Localizamos la cabeza del quinto metatarsiano mediante la técnica de palpación, y señalamos con el marcador el punto más proximal y distal (puntos E y F). Figura 2.18. Puntos de referencia para la pedigrafía. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 32 2.1.9. Señalamos con el marcador los espacios interdigitales del primer y segundo metatarsiano en su zona proximal (puntos G y H). 2.1.10. Se debe levantar toda la planta del pie a la vez mientras se sujeta el pedígrafo para que no se mueva. 2.1.11. Se debe evitar en todo momento el contacto de los dedos del explorador con el tapiz de caucho, para evitar marcas innecesarias que afecten a la huella. 2.2) Una vez realizadas las marcas de referencia, se dibujan las líneas de referencia: 2.2.1. Dibujar el eje de los maléolos (1). 2.2.2. Dibujar la línea que une las referencias proximales de las cabezas del primer y del quinto metatarsiano (2). 2.2.3. Dibujar la línea que une las referencias distales de las cabezas del primer y del quinto metatarsiano (3). 2.2.4. Dibujar la línea tangente a la huella en retropié y antepié por la cara interna (4). 2.2.5. Dibujar la línea tangente a la huella en retropié y antepié por la cara externa (5). 2.2.6. Dibujar las líneas (6) y (7) perpendiculares a las líneas (4) y (5) tangentes a la huella en retropié. La intersección de estas dos líneas es el punto 𝑐 . Figura 2.19. Líneas de referencia en la pedigrafía para determinar a, b y c. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 33 2.2.7. Dibujar las líneas (8) y (9) que unen el punto 7𝑐 con los espacios interdigitales del primer y segundo meta (G y H). 2.2.8. Dibujar las diagonales (10 y (11) del cuadrilátero que forman las líneas (2), (3), (8) y (9). 2.2.9. Dibujar la línea (12) que une el punto 𝑐 con la intersección de las diagonales del cuadrilátero. Esta línea es el eje longitudinal del pie. 2.2.10. La intersección de la línea (12) con la línea (1) da lugar al punto de referencia 𝑏 . 2.2.11. La intersección de la línea (12) con la línea (3) da lugar al punto de referencia 𝑎 . Una vez obtenidos los 3 puntos de referencia 𝑎 , 𝑏 y 𝑐 se delimitan los brazos de palanca. El brazo de palanca de retropié 𝑟 se corresponde con la distancia comprendida entre los puntos 𝑏 y 𝑐 . Aunque la obtención del brazo de palanca de antepié 𝑅 es algo más compleja, ya que depende del punto de apoyo 𝑎 y este a su vez depende de la posición con la que se apoye el pie sobre la célula de carga, se va a considerar que este brazo de palanca se corresponde con la distancia comprendida entre los puntos 𝑎 y 𝑏 . Como se ha podido comprobar, las distancias obtenidas dependen en cierta medida de la manera en que se realicen las marcas durante la pedigrafía. Por ello, para cada sujeto se realizará un mínimo de 3 pedigrafías y se tomarán como valores finales los valores medios obtenidos en cada una de las pedigrafías. Ignacio López Alonso Trabajo Fin de Grado 34 3. Tratamiento de los datos en Matlab 3.1 Introducción Una vez se han registrado todas las mediciones, es preciso llevar a cabo un tratamiento de los datos para llegar al objetivo final del ensayo: la obtención de las propiedades viscosoelásticas del tríceps sural. Como se ha explicado en el Capítulo 1, las mediciones y el tratamiento de los datos se realizaban hasta hace poco tiempo en un mismo programa denominado Flexor. Sin embargo, la pérdida de versiones recientes del software ha llevado a la necesidad de diseñar un nuevo programa que se ha desarrollado en Matlab y que se describe en este Capítulo. El nuevo programa se ha diseñado exclusivamente para el tratamiento de los datos, por lo que su medición se debe llevar a cabo a través de algún software externo como el que aporta el fabricante de la célula de carga. En el Apartado 3.2 se describe el desarrollo teórico necesario para llevar a cabo los ajustes por mínimos cuadrados, así como una breve descripción del antiguo programa Flexor, para ofrecer al lector algunas nociones acerca de cómo se trabajaba anteriormente. En el Apartado 3.3 se describen las 4 funciones que componen el programa de Matlab. En el desarrollo del apartado no se incluye explícitamente el código de las funciones, sino se explican las órdenes que lo componen y se da la información necesaria para entender adecuadamente el funcionamiento de las mismas. El código que compone las funciones se ha incluido en forma de anexos al final del documento. 3.2 Descripción teórica del tratamiento de los datos 3.2.1 Primer ajuste por mínimos cuadrados. Ajuste de cada medición En el Capítulo 2 se han desarrollado las ecuaciones que permiten, a través del registro de la fuerza de reacción en la célula de carga, obtener los distintos parámetros que describen las propiedades viscoso-elásticas del TS. La Ecuación (3.1), que se muestra a continuación por comodidad en la lectura, es la ecuación que gobierna la fuerza de reacción que recoge la célula de carga. 𝐹,=𝐹+𝑀𝑔=𝑒&*" ( 𝐴+𝑠𝑒𝑛𝜔)𝑡+𝐵+𝑐𝑜𝑠𝜔)𝑡 ) +𝑀𝑔 (3.1) Esta ecuación contiene 5 parámetros ( 𝛾 , 𝐴+ , 𝐵+ , 𝑀 y 𝜔) ) que, a priori, son desconocidos, y cuya determinación supone el conocimiento de las propiedades viscoso-elásticas del sistema Ignacio López Alonso Trabajo Fin de Grado 35 dinámico equivalente que representa al TS como un conjunto. Teniendo en cuenta que los parámetros 𝛾 , 𝑀 y 𝜔)7 se pueden estimar gráficamente, es relativamente sencillo realizar un ajuste por mínimos cuadrados con el que obtener, cometiendo un error aceptable, unos parámetros que se correspondan con la curva recogida por la célula de carga. Para ello, se genera una función de error cuya minimización proporciona los parámetros buscados. La función se genera de la siguiente manera: - Elección de un rango temporal del registro en el cual se va a realizar el ajuste. Se deben evitar los puntos muy cercanos al impacto, así como aquellos en los que la vibración esté muy amortiguada. - Para cada punto registrado experimentalmente, del que se obtiene el par [ 𝑡,𝐹(𝑡) ], se evalúa el valor de 𝐹(𝑡) en la expresión, que queda en función de las 5 incógnitas ( 𝛾 , 𝐴+ , 𝐵+ , 𝑀 y 𝜔) ). Se denominará al valor obtenido experimentalmente 𝐹/01 , y 𝐹 al evaluado analíticamente en cada 𝑡 =𝑡/01 . - Se calculará el error cuadrático en cada punto como la diferencia entre el valor experimental y el analítico elevada al cuadrado. Nótese que se utiliza el error cuadrático para evitar compensaciones de signo. - La función error se obtiene como la suma de los errores cuadráticos de todos los puntos del rango seleccionado. 𝑒𝑟𝑟𝑜𝑟 =V[𝐹/01,𝑡/010−𝐹(𝐴+,𝐵+,𝛾,𝜔),𝑀)]%777 2 34. (3.2) Dado que esta función de error puede tener varios mínimos relativos, será preciso concretar unas condiciones iniciales en la minimización de esta función para asegurarnos de que se obtiene el mínimo absoluto. Los parámetros 𝐴+ y 𝐵+ , aunque están relacionados con la amplitud de la señal, no tienen una interpretación gráfica clara, por lo que se estimarán los otros 3 parámetros, dejando libres 𝐴+ y 𝐵+ . Figura 3.1. Ejemplo de señal que se obtiene de una medición. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 36 El coeficiente de amortiguamiento 𝛾 es el único responsable del amortiguamiento de la función periódica [𝐴+sin ( 𝜔)𝑡 ) +𝐵+cos7(𝜔)𝑡)] a través del término 𝑒&*" (Figura 3.2). Su valor puede estimarse, por tanto, mediante la variación de 𝐹(𝑡) entre dos máximos consecutivos: 𝛾 ≈lnL𝐹(𝑡) 𝐹(𝑡+𝑇)M 𝑇 (3.3) Donde 𝑇 es el período. Al tratarse de ondas experimentales, dicho valor puede variar entre distintos máximos, por lo que se puede evaluar el parámetro entre cada dos máximos consecutivos y realizar un promedio de los valores obtenidos para obtener un valor más preciso. La frecuencia angular también es un parámetro muy fácil de evaluar de manera gráfica, simplemente a través del período, convirtiendo convenientemente los valores experimentales en valores de frecuencia: 𝜔)=2𝜋 𝑇 (3.4) También se debe evaluar el período entre distintos máximos y hacer el promedio de las frecuencias obtenidas. Por último, la masa que oscila 𝑀 se puede obtener como el valor medio de la carga antes del impacto, que debe ser igual al que se tiene cuando la oscilación se amortigua. Es sencillo estimar este parámetro a simple vista. Por ejemplo, en la gráfica de la Figura 3.1, el valor de 𝑀 es aproximadamente 30 kg. Figura 3.2. Amortiguamiento de una función periódica según la exponencial negativa. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 37 3.2.2 Segundo ajuste por mínimos cuadrados. Obtención de las rigideces del Sóleo y del tendón de Aquiles Como se ha explicado en el Capítulo 2, aunque es interesante conocer las propiedades viscoso-elásticas del Tríceps Sural como conjunto, el objetivo de un ensayo completo es la obtención de las propiedades viscoso-elásticas de cada uno de los elementos que forman el TS por separado. Recordamos que el Tendón de Aquiles presenta una rigidez constante, mientras que la rigidez del Sóleo varía linealmente con la fuerza que lo solicita. Ambos valores de rigidez se relacionan con la rigidez del TS como conjunto a través de: 𝑘 = 𝑘"𝑘--𝑓 𝑘"+𝑘--𝑓 (3.5) Donde 𝑘" es la rigidez constante del tendón y 𝑘-- es la pendiente de rigidez del Sóleo. Si se representa la rigidez global 𝑘 frente a la fuerza que solicita al TS 𝑓 , se obtiene una curva similar a la que muestra la Figura 3.3, donde se aprecia que la rigidez del Sóleo se corresponde con la pendiente de dicha curva en el origen, y la rigidez del tendón se corresponde con la asíntota a la que tiende la curva en el infinito. Esta representación tiene una interpretación física muy clara. Al suponerse el tendón y el Sóleo dos elementos dinámicos en serie, cuando la fuerza que solicita al TS es baja, la rigidez que gobierna la dinámica del MTC es la rigidez del Sóleo, ya que el tendón se comporta como un sólido rígido ante una carga tan pequeña. Sin embargo, a medida que aumenta la fuerza, es la rigidez del tendón la que gobierna el movimiento, ya que la rigidez de Sóleo se hará tan alta que este se comportará como un sólido rígido. Esta es la razón por la que es necesario repetir los ensayos aplicando diferentes cargas concentradas. Por tanto, para cada carga se obtendrá una pareja de resultados ( 𝑓,𝑘 ), obteniéndose una nube de puntos como la que muestra la Figura 3.4. Téngase en cuenta que 𝑓 es el valor promedio de la carga que oscila en cada medición, por lo que se puede calcular como el valor de Figura 3.3. Representación gráfica de la rigidez del Tríceps Sural frente a la fuerza que lo solicita. (París-García et al., 2010) Ignacio López Alonso Trabajo Fin de Grado 38 la carga antes del impacto (valor 𝑀 calculado en el apartado anterior) en sus debidas unidades (se mide en kg, hay que pasarlo a Newton), y 𝑘 se puede calcular como: 𝑘 = L 𝑅 𝑟 M %𝑀(𝜔) %+𝛾%) (3.6) Cabe destacar que será necesario hacer más de una medición para cada carga con el fin de obtener valores promedios. Para obtener los valores individuales de rigidez a partir de la nube de puntos, se puede realizar un ajuste por mínimos cuadrados similar al realizado para obtener los parámetros de cada medición. Para ello, se genera una función error que toma la siguiente forma: 𝑒𝑟𝑟𝑜𝑟 =V[𝑘/01,𝑓/010−𝑘(𝑘" 2 34. ,𝑘--)]% (3.7) La minimización de esta función de error proporcionará los valores de 𝑘" y 𝑘-- que ajustan la curva y representan las rigideces individuales del tendón y el Sóleo respectivamente. 3.2.3 Breve descripción del programa Flexor Flexor es un software de apariencia comercial que permite llevar a cabo todo el proceso de recogida y tratamiento de los datos de manera sencilla, con buen grado de automatización y sin necesidad de que el usuario conozca el fundamento físico del método. El programa, desarrollado en Java, contiene varias ventanas y funciones que se explican a continuación. Figura 3.4. Ejemplo de ajuste por mínimos cuadrados de la función de rigidez. (París-García et al., 2015) Ignacio López Alonso Trabajo Fin de Grado 39 En primer lugar, se distinguen dos ventanas (zona roja) destinadas a la gestión de los archivos. Permiten crear un nuevo proyecto, abrir proyectos y mediciones ya existentes y generar el informe final que proporcionará toda la información necesaria para el conocimiento de las propiedades viscoso-elásticas del tríceps sural del sujeto. En color azul se aprecia otra ventana cuya función es el control del dispositivo de medida. Esta ventana permite activar y desactivar la recogida de datos, y muestra por pantalla en una gráfica la medición de los datos en tiempo real. La zona recuadrada en verde es la más interesante. En ella se distinguen dos ventanas que muestran la señal original recogida y la señal ajustada por mínimos cuadrados. Este ajuste depende de 5 parámetros (ya descritos en el Capítulo 2) cuyo valor aparece en la ventana “Parámetros”. Además, se incluye una casilla extra por cada parámetro que permite modificar los valores calculados por el sistema por si fuera necesario realizar algún ajuste manual. Es posible también modificar el filtro que se aplica sobre la señal. Por último, el programa incluye una ventana denominada “Segundo ajuste” en la que se pueden introducir las medidas de los brazos de palanca 𝑅 y 𝑟 necesarios para relacionar los datos obtenidos con los parámetros que definen las propiedades del MTC. 3.3 Descripción numérica del tratamiento de los datos. Programa de Matlab Para la implementación de los dos ajustes por mínimos cuadrados necesarios en Matlab se han creado 4 funciones que permiten obtener los resultados de manera relativamente automática y sencilla. En este apartado se explica el funcionamiento de las funciones y se proporciona la información necesaria para su entendimiento. En los Anexos I, II, III y IV se proporciona el código de Matlab de las funciones (con las líneas numeradas), en el que aparece Figura 3.5. Apariencia visual del antiguo programa Flexor. Ignacio López Alonso Trabajo Fin de Grado 40 una serie de comentarios que facilitan la comprensión de todos los comandos. Estas 4 funciones son: - Función de lectura: lee todos los ficheros correspondientes a un ensayo completo y almacena los datos en una matriz. Estos datos son pares de valores [𝑡,𝐹 ( 𝑡 ) ] que definen cada onda subamortiguada. - Función de filtro: atenúa las vibraciones de alta frecuencia que factores externos (vibración del equipo y del sujeto principalmente) inducen sobre cada señal original provocando efectos de ruido. - Función de primer ajuste por mínimos cuadrados: con los datos que le proporciona la función de lectura ajusta la curva y obtiene los 5 parámetros correspondientes al primer ajuste para cada medición. - Función de segundo ajuste por mínimos cuadrados: con los datos que le proporciona la función de primer ajuste obtiene la curva de rigidez. Se describe a continuación de manera detallada el contenido y el funcionamiento de cada una de las funciones que componen el programa. 3.3.1 Función de lectura Los ficheros de texto que definen las gráficas de cada una de las mediciones están formados por dos columnas numéricas en las que aparecen los pares [𝑡,𝐹 ( 𝑡 ) ] que las componen. Cuando se abordó la implementación del proceso en Matlab, una de las principales incógnitas era si se podría hacer un programa manejable y si se iba a poder obtener la gráfica final de rigidez sin necesidad de introducir manualmente el nombre de cada uno de los archivos. Para hacer esto posible, se ha diseñado la función de manera que, siendo la entrada un fichero de texto en el que aparezcan los datos de todas las mediciones, una detrás de otra sin separación, el programa separe cada medición y almacene los datos en dos matrices t_glob y F_glob, que son las salidas de la función. Así, la función tiene la forma: [t_glob,F_glob]=lectura(FileName) La unión de los ficheros se puede llevar a cabo en Windows de la siguiente forma: 1. Copiar en el portapapeles el directorio de la carpeta en la que se encuentran los ficheros. En esta carpeta se deben guardar únicamente los ficheros que se quieren unir, que deben tener formato .txt. 2. Presionar la tecla Windows + R. 3. Escribir cmd en la ventana que aparece y presionar la tecla Enter. 4. Escribir cd en la nueva ventana que aparece y pegar el directorio de la carpeta en la que se encuentran los ficheros dejando un espacio de separación (para cambiar el directorio al que se apunta). Presionar la tecla Enter. 5. Escribir el siguiente comando y presionar la tecla Enter: copy *.txt union.txt Ignacio López Alonso Trabajo Fin de Grado 47 El procedimiento se puede entender con mayor claridad consultando el código en el Anexo III. Como se comentaba en el Apartado 3.2.1, las idea a la hora de minimizar la función 𝑒𝑟𝑟𝑜𝑟 es estimar los parámetros 𝛾 , 𝜔) y 𝑀 para introducirlos como condiciones iniciales, y dejar libres los parámetros 𝐴+ y 𝐵+ . Sin embargo, la función fminsearch no permite dejar libre ningún parámetro. Para evitar asignar a 𝐴+ y 𝐵+ valores iniciales aleatorios que puedan dar lugar a algún tipo de error o a una ralentización del programa, se han creado dos funciones de error: - Una primera función en la que las únicas variables son 𝐴+ y 𝐵+ , asignando a 𝛾 , 𝜔) y 𝑀 el valor constante que se ha estimado, con la que se estiman los valores de 𝐴+ y 𝐵+ . - Una segunda función que se corresponde con la función 𝑒𝑟𝑟𝑜𝑟 propiamente dicha, en la que se tienen como variables los 5 parámetros que definen la función. El objetivo de este paso es minimizar la función, en primer lugar, con 𝐴+ y 𝐵+ como únicas variables para estimar sus valores, y en segundo lugar, ahora que ya se tiene una estimación del valor de los 5 parámetros, minimizar la función con todos los parámetros como variables. La razón por la que se hace esta minimización en dos pasos es que la función error puede contar con varios mínimos relativos, por lo que si no se parte de las condiciones iniciales adecuadas, se puede obtener un resultado que diste mucho de ser el mínimo absoluto. Por último, es preciso recordar que el objetivo de la minimización es la obtención de 5 parámetros por medición, de los cuales se obtendrán las propiedades viscoso-elásticas del MTC. Por tanto, estos parámetros se almacenarán en una matriz 5xN (donde N es el número de mediciones) que será la entrada de la función de segundo ajuste que se describe en el Apartado 3.3.4. 3.3.4 Función de segundo ajuste por mínimos cuadrados. Obtención de las rigideces del Sóleo y del tendón de Aquiles El objetivo de esta función es aproximar la nube de puntos que forman los pares de valores [𝑓,𝑘] que se obtienen de cada medición a través de una curva que tiene la siguiente forma: 𝑘 = 𝑘"𝑘--𝑓 𝑘"+𝑘--𝑓 (3.9) Se recuerda que el par [𝑓,𝑘] para cada medición se puede obtener como: 𝑓 =𝑀𝑔 (3.10) 𝑘 =L𝑅 𝑟M%𝑀(𝜔) %+𝛾%) (3.11) Ignacio López Alonso Trabajo Fin de Grado 48 Se deben, por tanto, aportar como entradas de la función de segundo ajuste, la matriz que contiene los 5 parámetros de cada medición, y las longitudes de los brazos de palanca 𝑅 y 𝑟 que se obtienen de las pedigrafías. La función tomará, por tanto, la siguiente forma: [k_t,k_ss]=rigidez(X_glob,R,r) Donde X_glob es la matriz de los parámetros. De esta función se obtendrá una gráfica que toma la forma que muestra la Figura 3.9. Para llevar a cabo este ajuste, se deberá minimizar la función 𝑒𝑟𝑟𝑜𝑟 que describe la Ecuación (3.12), y que se construirá de la misma forma que se ha hecho con la función 𝑒𝑟𝑟𝑜𝑟7 del primer ajuste en el Apartado 3.3.3; es decir, se montará, en primer lugar, a través de un sumatorio en simbólico, y se convertirá en una función tipo “handle” por medio de cadenas de caracteres. 𝑒𝑟𝑟𝑜𝑟 =V[𝑘/01,𝑓/010−𝑘(𝑘" 2 34. ,𝑘--)]% (3.12) Esta función también se minimizará con fminsearch, y en este caso no se necesitará estimar unas condiciones iniciales ya que se trata de una función con únicamente 2 variables, por lo que el ajuste se realizará adecuadamente en la gran mayoría de los casos. La salida de esta función de segundo ajuste son los parámetros 𝑘" y 𝑘-- que describen la rigidez del tendón de Aquiles y del Sóleo respectivamente, y que son el objetivo principal de estos ensayos. Figura 3.9. Representación gráfica de la rigidez del Tríceps Sural frente a la fuerza que lo solicita. Ignacio López Alonso Trabajo Fin de Grado 49 Ignacio López Alonso Trabajo Fin de Grado 50 4. Verificación y visualización de resultados. Conclusiones 4.1 Introducción En el Capítulo 3 se ha descrito con detalle el programa de Matlab que lleva a cabo el ajuste por mínimos cuadrados de cada medición y la obtención de las propiedades viscosoelásticas del tríceps sural, pero no se ha incluido nada acerca de la validez ni de la visualización de los resultados. En este Capítulo se describirán las funciones encargadas de comprobar si el ajuste de mínimos cuadrados para cada onda es válido, y se mostrará y describirá todo el apartado visual que genera el programa. También se incluye un apartado de conclusiones. 4.2 Verificación de los ajustes Cuando se realiza un ajuste por mínimos cuadrados, no es suficiente con observar visualmente las gráficas original y del ajuste superpuestas para comprobar si el ajuste es válido o no, sino que se debe evaluar dicha validez a través de algún parámetro. En este caso, se ha utilizado el conocido como coeficiente de determinación, que se define según la siguiente expresión: 𝑅%=1−𝑆𝑆𝑅 𝑆𝑆𝑇 (4.1) donde: 𝑆𝑆𝑅 =V[𝑦g(𝑥3)−𝑦(𝑥3)]% 3 (4.2) 𝑆𝑆𝑇 =V[𝑦g(𝑥3)−𝑦h(𝑥3)]% 3 (4.3) Y donde 𝑦 !( 𝑥𝑖 ) representa el valor de la función de ajuste en 𝑥3 , 𝑦(𝑥3) representa el valor experimental de la señal en 𝑥3 , y 𝑦 $( 𝑥𝑖 ) se corresponde con el valor medio de todos los puntos que componen la señal experimental. Ignacio López Alonso Trabajo Fin de Grado 51 Así, el parámetro 𝑆𝑆𝑅 se puede definir como la suma de las diferencias cuadráticas entre el punto y la función evaluada en la 𝑥3 de dicho punto, mientras que el parámetro 𝑆𝑆𝑇 se define como la suma de las diferencias cuadráticas entre el valor medio de todos los puntos y la función evaluada en la 𝑥3 de cada punto. Se considera que valores del coeficiente de determinación mayores que 0.9 informan de un ajuste aceptable. Se ha incluido en las funciones de primer y segundo ajuste descritas en el Capítulo 3 esta medida del error para comprobar que los resultados obtenidos son válidos. En el caso de que el ajuste sobre alguna medición no sea válido ( 𝑅%<0.9 ), el programa informará de qué medición se trata y la desechará, por lo que no se tendrá en cuenta para el cálculo de la rigidez en el segundo ajuste. Además, este programa incluye una mejora considerable respecto a Flexor, y es la capacidad de adaptar el rango temporal en el que se realiza el ajuste en función del valor del coeficiente de determinación 𝑅% obtenido. Como se comenta en el Capítulo 3, el programa se ha diseñado para realizar el ajuste por mínimos cuadrados, en principio, en un rango temporal comprendido entre el segundo y el tercer máximo de cada medición. Sin embargo, el programa no siempre consigue realizar el ajuste entre dichos máximos con un valor de 𝑅%>0.9 , pero sí podría ser capaz de realizarlo entre otros dos máximos. Por ello, se han incluido en las funciones varias subrutinas que, si se ha fracasado en el primer intento, tratan de realizar de nuevo el ajuste modificando el rango temporal, ahora comprendido entre el primer y el segundo máximo. Si también se fracasa en este intento, se intentará ajustar entre el primer y el tercer máximo, desechando definitivamente la onda si este intento también fracasa. La diferencia con Flexor es que estaba diseñado para intentar realizar el ajuste en una ventana temporal predeterminada, y si no era capaz de ajustar la onda en dicha ventana, era necesario bien modificar manualmente los parámetros para conseguir un buen ajuste (cosa improbable), o bien repetir la medición. La inclusión de estas subrutinas confiere al nuevo programa gran potencia y adaptabilidad, ya que, aunque el programa se ha configurado conforme a lo descrito en el párrafo anterior, con pequeñas modificaciones en el código se pueden conseguir ajustes en los rangos que se desee, y así hacer incluso una comparativa entre las distintas ventanas para que el ajuste sea óptimo en cada medición. Para ilustrar la posibilidad de obtener buenos y malos ajustes en una misma onda, la Figura 4.1 muestra una medición en la que se consigue un buen ajuste entre el primer y el segundo máximo, y un mal ajuste entre los dos máximos siguientes. Tras procesar todas las mediciones correspondientes a un ensayo completo, el programa sacará por pantalla un mensaje que informe de si alguna medición no se ha ajustado correctamente, o si todos los ajustes son correctos. Para desechar los ajustes que no sean correctos se ha incluido en la función de segundo ajuste un apartado que recalcula los vectores que contienen los valores (𝑓,𝑘) , eliminando los derivados de aquellas mediciones en las que 𝑅%<0.9 . Ignacio López Alonso Trabajo Fin de Grado 52 Téngase en cuenta que en un ensayo completo se realizan unas 4 mediciones con cada peso con el fin de obtener un valor promedio de la rigidez que se obtiene con cada peso, por lo que el hecho de desechar alguna medición no debe afectar al resultado final de manera notable. Aún así, si no se consigue un buen ajuste para una determinada onda, es conveniente repetir dicha medición, ya que lo más probable es que la onda presente irregularidades. 4.3 Visualización de los resultados Aunque ya se haya comprobado numéricamente que los ajustes son correctos, y se hayan obtenido los resultados numéricos correspondientes a las rigideces del Sóleo y del tendón de Aquiles, es bueno proyectar gráficamente todos los resultados para un mayor entendimiento y una visión más global de los mismos. En este Apartado se describe el contenido visual que genera Matlab cuando se ejecutan las funciones. 4.3.1 Ajuste de cada medición La ejecución de la función de primer ajuste por mínimos cuadrados devolverá, además de un mensaje informando de las mediciones incorrectas, una ventana en la que aparecerán representadas las gráficas correspondientes a las mediciones que componen el ensayo con su ajuste superpuesto. En cada gráfica se incluirá un cuadro de texto que informa del valor del coeficiente de determinación 𝑅% . Si el ajuste es aceptable, el texto aparecerá en color verde, y si no es aceptable, aparecerá en color rojo. En la Figura 4.2 se muestra el formato de presentación de estos resultados, en el que únicamente se han incluido 8 gráficas para mayor claridad. Normalmente se generará una ventana con 32 gráficas (8 pesos x 4 repeticiones). Aunque no es lo habitual, aparece en la Figura 4.2 alguna medición incorrecta para dar una idea del aspecto visual en caso de que esto ocurra. Figura 4.1. Comparativa del ajuste sobre una medición en distintas ventanas temporales. a) Buen ajuste en ventana temporal comprendida entre el primer y el segundo máximo. b) Mal ajuste entre el segundo y tercer máximo. Ignacio López Alonso Trabajo Fin de Grado 53 4.3.2 Resultados del estudio La ejecución de la función de segundo ajuste por mínimos cuadrados emitirá una especie de informe sencillo que contendrá la información relevante correspondiente a los resultados del estudio (ensayo completo). Este informe está formado por dos elementos: - Una tabla (Figura 4.3) en la que se recoge el par de valores (𝑓,𝑘) que se obtiene de cada medición, así como el valor del coeficiente de determinación con el que se ha ajustado la Figura 4.2. Formato de representación gráfica de las mediciones y sus ajustes. Ignacio López Alonso Trabajo Fin de Grado 54 medición. Las celdas que contienen el valor del 𝑅% aparecerán en color verde si el ajuste es bueno, y en color rojo si el ajuste no es aceptable. que se incluye en la tabla la información de las mediciones que no se han ajustado correctamente, aunque no se han tenido en cuenta para el cálculo de las propiedades viscoso-elásticas del Sóleo y el tendón de Aquiles. - Una gráfica (Figura 4.4) que muestra las propiedades del Sóleo y del tendón por separado, obtenida del ajuste por mínimos cuadrados de los puntos (𝑓,𝑘) según una curva de la forma que describe la Ecuación (3.9). Figura 4.4. Tabla resumen que contiene los valores de f, k y R2 para cada medición. Figura 4.3. Representación gráfica de las propiedades viscoso-elásticas del Sóleo y el tendón de Aquiles. Ignacio López Alonso Trabajo Fin de Grado 55 Ignacio López Alonso Trabajo Fin de Grado 56 5. Conclusiones El desarrollo de este proyecto ha concluido con la creación de una nueva herramienta que permite, de manera relativamente sencilla, cómoda y autónoma, llevar a cabo el tratamiento de los datos que se obtienen durante la aplicación de la técnica de vibración libre, y visualizarlos de manera bastante intuitiva. Además, el hecho de contar con el código detallado y explicado permite a cualquier usuario que tenga conocimientos básicos de programación entender todo el proceso y modificar el procedimiento de ajuste para adaptar el programa a cualquier situación que lo requiera. Así, gracias a esta nueva herramienta, será posible continuar con la puesta a punto del nuevo equipo diseñado para una aplicación optimizada de la técnica de vibración libre, que fue la motivación inicial de este proyecto. Sin embargo, aunque el programa cumple con su función, se podrían incluir ciertas modificaciones y mejoras que lo harían más robusto, intuitivo y visual. Se detallan a continuación las 3 posibles mejoras que se han considerado más relevantes: - Implementación de alguna subrutina en la función de lectura que sea capaz de recorrer automáticamente todos los ficheros que componen el ensayo y crear el archivo que contiene todos los ficheros, o incluso ir accediendo individualmente a cada fichero de forma automática para ajustar cada medición. Con esto, se evitaría el paso de creación del fichero que contiene todas las mediciones a través de los comandos de Windows descritos en el Apartado 3.3.1, que resulta algo engorrosa. - Creación, tras la ejecución de la función de segundo ajuste por mínimos cuadrados, de un fichero PDF que contenga un informe mejor estructurado y más formal sobre las propiedades del tríceps sural y la calidad de los ajustes. Esto se podría hacer a través de la GUI de Matlab. - Mejora de la robustez de los ajustes. Como se ha comentado a lo largo del documento, la función de error con la que se realizan los ajustes puede contener varios mínimos relativos. Si se representa la función error respecto a la frecuencia (manteniendo el resto de los parámetros constantes), se obtiene una gráfica con la forma de la Figura 5.1, en la que se observa que suele haber dos valores de frecuencia que proporcionan prácticamente el mismo error. Sin embargo, la dependencia de la rigidez del Tríceps Sural con el cuadrado de la frecuencia hace que, aunque se obtenga un error muy parecido con ambos valores de frecuencia, los valores de rigidez puedan ser ligeramente distintos. Aunque gracias a la estimación de la frecuencia que se hace al comienzo de la función, el ajuste se realiza con el valor de frecuencia adecuado en la gran mayoría de los casos, se podría incluir una subrutina que compruebe que el mínimo hallado es el adecuado, y modifique el ajuste en caso negativo. Por último, a nivel personal, el proyecto me ha generado gran interés por el conocimiento de las técnicas utilizadas para la medición de ciertas propiedades biomecánicas, y me siento Ignacio López Alonso Trabajo Fin de Grado 63 238 end 239 end 240 if flag == itercomp 241 val_inic(j)=F(i+paso); 242 t_inic(j)=i; 243 j=j+1; 244 end 245 end 246 247 %cálculo de la media M 248 249 M_0=sum(F(1:t_inic(1)))/t_inic(1); 250 251 %cálculo de los máximos 252 253 paso=1; 254 num_max=3; %número de máximos que se quiere encontrar 255 max=zeros(1,num_max); 256 j=1; 257 i=t_inic(1)+10; 258 while j<(num_max+1) && i<(length(F)-3*paso) 259 if F(i)/F(i+paso)>1 && F(i)/F(i-paso)>1 260 flag1=1; 261 for m=1:20 %evita falsos máximos 262 if F(i)/F(i+m)>1 && F(i)/F(i-m)>1 263 else 264 flag1=0; 265 end 266 end 267 if flag1==1 268 max(j)=F(i); 269 j=j+1; 270 end 271 end 272 i=i+paso; 273 end 274 275 tmax=zeros(1,num_max); 276 for j=1:num_max %calcula el tiempo correspondiente a los máximos 277 for i=1:length(F) 278 if F(i)==max(j) 279 tmax(j)=i; 280 end 281 end 282 end 283 284 %parámetro w_d 285 286 for i=1:num_max-1 287 T_aux(i)=t(tmax(i+1))-t(tmax(i)); 288 end 289 T=sum(T_aux)/length(T_aux); 290 w_d_0=2*pi/T; 291 292 293 294 %parámetro gamma 295 296 for i=1:num_max-1 297 gamma_aux(i)=log(max(i)/max(i+1))/(T_aux(i)); 298 Ignacio López Alonso Trabajo Fin de Grado 64 end 299 gamma_0=sum(gamma_aux)/length(gamma_aux); 300 301 302 %mínimos cuadrados 303 %intenta ajustar la función entre el segundo y tercer máximo 304 305 syms A_F B_F gamma1 Mvar w_d 306 t_r=t(tmax(2)-20:tmax(num_max)+40); %rango temporal algo mayor que 307 F_r=F(tmax(2)-20:tmax(num_max)+40); %entre dos máximos 308 error=0; 309 310 %se construye la función como un sumatorio en simbólico 311 312 for i=1:length(t_r) 313 error=error+(F_r(i)-(exp(-gamma1*(t_r(i)))*(A_F*sin(w_d*(t_r(i)))) 314 +B_F*cos(w_d*(t_r(i)))+Mvar))^2; 315 end 316 317 erroraux=evalc('disp(error)'); %convierte la funcion en vector char 318 errorauxorig=convertCharsToStrings(erroraux); %convierte a string 319 erroraux=replace(errorauxorig,'A_F','x(1)'); %sustituye las variables 320 erroraux=replace(erroraux,'B_F','x(2)'); 321 erroraux=replace(erroraux,'gamma1','x(3)'); 322 erroraux=replace(erroraux,'Mvar','x(4)'); 323 erroraux=replace(erroraux,'w_d','x(5)'); 324 straux="@(x)"; 325 straux2=[straux,erroraux];%mete el @(x) para que al evaluarlo sea 326 %funcion handle 327 erroraux=join(straux2); 328 error=eval(erroraux); %ya tenemos la funcion handle 329 330 %función error_1 para un primer ajuste para estimar A_F y B_F 331 332 erroraux_1=replace(errorauxorig,'A_F','x(1)'); 333 erroraux_1=replace(erroraux_1,'B_F','x(2)'); 334 erroraux_1=replace(erroraux_1,'gamma1','gamma_0'); 335 erroraux_1=replace(erroraux_1,'Mvar','M_0'); 336 erroraux_1=replace(erroraux_1,'w_d','w_d_0'); 337 straux_1="@(x)"; 338 straux2_1=[straux_1,erroraux_1];%mete el @(x) para que al evaluarlo 339 %sea funcion handle 340 erroraux_1=join(straux2_1); 341 error_1=eval(erroraux_1); %ya tenemos la funcion handle 342 343 %Opciones de búsqueda. Se dejan ajustes opcionales como comentario 344 345 %options=optimset('MaxFunEvals',1e10,'MaxIter',1e10,'TolX',1e-10); 346 options=optimset('MaxFunEvals',1e10); 347 348 349 %se estima A_F y B_F 350 351 [X_1]=fminsearch(error_1,[10,10],options); 352 353 %se minimiza la función partiendo de los parámetros estimados 354 355 [X]=fminsearch(error,[X_1(1),X_1(2),gamma_0,M_0,w_d_0],options); 356 357 358 Ignacio López Alonso Trabajo Fin de Grado 65 %cálculo del coeficiente de determinación 359 360 media=sum(F_r)/length(F_r); 361 SER=0; 362 STR=0; 363 for i=1:length(t_r) 364 SER = SER + ((exp(-X(3).*(t_r(i))).*(X(1).*sin(X(5).*(t_r(i))) 365 +X(2).*cos(X(5).*(t_r(i))))+X(4))-F_r(i))^2; 366 end 367 368 for i=1:length(t_r) 369 STR = STR + ((exp(-X(3).*(t_r(i))).*(X(1).*sin(X(5).*(t_r(i))) 370 +X(2).*cos(X(5).*(t_r(i))))+X(4))-media)^2; 371 end 372 373 R2=1-SER/STR; 374 375 376 X_glob(:,count)=X'; 377 378 379 %Si el ajuste no es bueno, se repite el ajuste modificando el rango 380 %temporal. Rango temporal actual: entre primer y segundo máximo 381 382 if R2<0.9 383 384 t_r=t(tmax(1)-20:tmax(num_max-1)+40); %rango temporal algo mayor que 385 F_r=F(tmax(1)-20:tmax(num_max-1)+40); %entre dos máximos 386 error=0; 387 388 %se construye la función como un sumatorio en simbólico 389 390 for i=1:length(t_r) 391 error=error+(F_r(i)-(exp(-gamma1*(t_r(i)))*(A_F*sin(w_d*(t_r(i)))) 392 +B_F*cos(w_d*(t_r(i)))+Mvar))^2; 393 end 394 395 erroraux=evalc('disp(error)'); %convierte la funcion en vector char 396 errorauxorig=convertCharsToStrings(erroraux); %convierte a string 397 erroraux=replace(errorauxorig,'A_F','x(1)'); %sustituye las variables 398 erroraux=replace(erroraux,'B_F','x(2)'); 399 erroraux=replace(erroraux,'gamma1','x(3)'); 400 erroraux=replace(erroraux,'Mvar','x(4)'); 401 erroraux=replace(erroraux,'w_d','x(5)'); 402 straux="@(x)"; 403 straux2=[straux,erroraux];%mete el @(x) para que al evaluarlo sea 404 %funcion handle 405 erroraux=join(straux2); 406 error=eval(erroraux); %ya tenemos la funcion handle 407 408 %función error_1 para un primer ajuste para estimar A_F y B_F 409 410 erroraux_1=replace(errorauxorig,'A_F','x(1)'); 411 erroraux_1=replace(erroraux_1,'B_F','x(2)'); 412 erroraux_1=replace(erroraux_1,'gamma1','gamma_0'); 413 erroraux_1=replace(erroraux_1,'Mvar','M_0'); 414 erroraux_1=replace(erroraux_1,'w_d','w_d_0'); 415 straux_1="@(x)"; 416 straux2_1=[straux_1,erroraux_1];%mete el @(x) para que al evaluarlo 417 %sea funcion handle 418 Ignacio López Alonso Trabajo Fin de Grado 66 erroraux_1=join(straux2_1); 419 error_1=eval(erroraux_1); %ya tenemos la funcion handle 420 421 %Opciones de búsqueda. Se dejan ajustes opcionales como comentario 422 423 %options=optimset('MaxFunEvals',1e10,'MaxIter',1e10,'TolX',1e-10); 424 options=optimset('MaxFunEvals',1e10); 425 426 427 %se estima A_F y B_F 428 429 [X_1]=fminsearch(error_1,[10,10],options); 430 431 %se minimiza la función partiendo de los parámetros estimados 432 433 [X]=fminsearch(error,[X_1(1),X_1(2),gamma_0,M_0,w_d_0],options); 434 435 %cálculo del coeficiente de determinación 436 437 media=sum(F_r)/length(F_r); 438 SER=0; 439 STR=0; 440 for i=1:length(t_r) 441 SER = SER + ((exp(-X(3).*(t_r(i))).*(X(1).*sin(X(5).*(t_r(i))) 442 +X(2).*cos(X(5).*(t_r(i))))+X(4))-F_r(i))^2; 443 end 444 445 for i=1:length(t_r) 446 STR = STR + ((exp(-X(3).*(t_r(i))).*(X(1).*sin(X(5).*(t_r(i))) 447 +X(2).*cos(X(5).*(t_r(i))))+X(4))-media)^2; 448 end 449 450 R2=1-SER/STR; 451 452 453 X_glob(:,count)=X'; 454 end 455 456 %Si el ajuste no es bueno, se repite modificando el rango temporal. 457 %Rango temporal actual: entre primer y tercer máximo 458 459 if R2<0.9 460 461 462 t_r=t(tmax(1)-20:tmax(num_max)+40); %rango temporal algo mayor que 463 F_r=F(tmax(1)-20:tmax(num_max)+40); %entre tres máximos 464 error=0; 465 466 %se construye la función como un sumatorio en simbólico 467 468 for i=1:length(t_r) 469 error=error+(F_r(i)-(exp(-gamma1*(t_r(i)))*(A_F*sin(w_d*(t_r(i)))) 470 +B_F*cos(w_d*(t_r(i)))+Mvar))^2; 471 end 472 erroraux=evalc('disp(error)'); %convierte la funcion en vector char 473 errorauxorig=convertCharsToStrings(erroraux); %convierte a string 474 erroraux=replace(errorauxorig,'A_F','x(1)'); %sustituye las variables 475 erroraux=replace(erroraux,'B_F','x(2)'); 476 erroraux=replace(erroraux,'gamma1','x(3)'); 477 Ignacio López Alonso Trabajo Fin de Grado 67 erroraux=replace(erroraux,'Mvar','x(4)'); 478 erroraux=replace(erroraux,'w_d','x(5)'); 479 straux="@(x)"; 480 straux2=[straux,erroraux];%mete el @(x) para que al evaluarlo sea 481 %funcion handle 482 erroraux=join(straux2); 483 error=eval(erroraux); %ya tenemos la funcion handle 484 485 %función error_1 para un primer ajuste para estimar A_F y B_F 486 487 erroraux_1=replace(errorauxorig,'A_F','x(1)'); 488 erroraux_1=replace(erroraux_1,'B_F','x(2)'); 489 erroraux_1=replace(erroraux_1,'gamma1','gamma_0'); 490 erroraux_1=replace(erroraux_1,'Mvar','M_0'); 491 erroraux_1=replace(erroraux_1,'w_d','w_d_0'); 492 straux_1="@(x)"; 493 straux2_1=[straux_1,erroraux_1];%mete el @(x) para que al evaluarlo 494 sea funcion handle 495 erroraux_1=join(straux2_1); 496 error_1=eval(erroraux_1); %ya tenemos la funcion handle 497 498 %Opciones de búsqueda. Se dejan ajustes opcionales como comentario 499 500 %options=optimset('MaxFunEvals',1e10,'MaxIter',1e10,'TolX',1e-10); 501 options=optimset('MaxFunEvals',1e10); 502 503 504 %se estima A_F y B_F 505 506 [X_1]=fminsearch(error_1,[10,10],options); 507 508 %se minimiza la función partiendo de los parámetros estimados 509 510 [X]=fminsearch(error,[X_1(1),X_1(2),gamma_0,M_0,w_d_0],options); 511 512 %cálculo del coeficiente de determinación 513 514 media=sum(F_r)/length(F_r); 515 SER=0; 516 STR=0; 517 for i=1:length(t_r) 518 SER = SER + ((exp(-X(3).*(t_r(i))).*(X(1).*sin(X(5).*(t_r(i))) 519 +X(2).*cos(X(5).*(t_r(i))))+X(4))-F_r(i))^2; 520 end 521 522 for i=1:length(t_r) 523 STR = STR + ((exp(-X(3).*(t_r(i))).*(X(1).*sin(X(5).*(t_r(i))) 524 +X(2).*cos(X(5).*(t_r(i))))+X(4))-media)^2; 525 end 526 527 R2=1-SER/STR; 528 529 530 X_glob(:,count)=X'; 531 end 532 533 %matriz de coeficientes de determinacion R2 534 535 R2_mat(count)=R2; 536 Ignacio López Alonso Trabajo Fin de Grado 68 537 %representación de cada medición 538 539 nexttile %pasa a la siguiente región de la ventana (nueva gráfica) 540 plot(t,F),shg %gráfica original 541 hold on 542 fplot=exp(- 543 X(3).*(t_r)).*(X(1).*sin(X(5).*(t_r))+X(2).*cos(X(5).*(t_r)))+X(4); 544 %ajuste 545 plot(t_r,fplot),shg 546 xlabel('Tiempo (s)') 547 ylabel('Carga (kg)') 548 549 if R2>0.9 550 tit=sprintf('R^{2}=%.3f',R2); %valor de R2 verde 551 title(tit,'Color','green'); 552 else 553 num_err(err_count)=count; 554 err_count=err_count+1; 555 tit=sprintf('R^{2}=%.3f',R2); %valor de R2 rojo 556 title(tit,'Color','red'); 557 end 558 559 end %end del bucle que recorre todas las mediciones 560 561 %órdenes para sacar el mensaje por pantalla que informa de las 562 %curvas que se han ajustado correctamente y las que no 563 564 if num_err==0 565 fprintf('Todas las curvas se han ajustado correctamente.\n') 566 end 567 if err_count==2 568 fprintf('La curva de la medición %d no se ha ajustado 569 correctamente.\n',num_err(1)) 570 end 571 if err_count>2 572 err_cad_i=''; 573 for i=1:err_count-1 574 if i<err_count-2 575 err_cad_i1=sprintf('%d, ',num_err(i)); 576 end 577 if i==err_count-2 578 err_cad_i1=sprintf('%d y ',num_err(i)); 579 end 580 if i==err_count-1 581 err_cad_i1=sprintf('%d',num_err(i)); 582 end 583 if i>0 && i<err_count 584 err_cad_aux=[err_cad_i,err_cad_i1]; 585 err_cad_i=join(err_cad_aux); 586 end 587 end 588 589 tit_aux='Las curvas de las mediciones '; 590 tit_aux1=' no se han ajustado correctamente.'; 591 tit_aux_tot=[tit_aux,err_cad_i,tit_aux1]; 592 tit_aux_tot1=join(tit_aux_tot); 593 disp(tit_aux_tot1) 594 end 595 596 end 597 Ignacio López Alonso Trabajo Fin de Grado 69 ANEXO IV. Función de segundo ajuste por mínimos cuadrados 598 599 function [k_t,k_ss,f_glob,k_glob]=rigidez(X_glob,R,r,num_err,R2_mat) 600 601 j=1; 602 603 %definición de vectores para cálculos intermedios 604 605 f_glob_plot=(R/r)*9.81.*X_glob(4,:); 606 k_glob_plot=(R/r)^2.*X_glob(4,:).*(X_glob(5,:).^2+X_glob(3,:).^2).*1e- 607 3; 608 609 f_glob_int=f_glob_plot; 610 k_glob_int=k_glob_plot; 611 612 %órdenes que desechan los malos ajustes 613 614 if num_err==0 615 f_glob=(R/r)*9.81.*X_glob(4,:); 616 617 k_glob=(R/r)^2.*X_glob(4,:).*(X_glob(5,:).^2+X_glob(3,:).^2).*1e-3; 618 else 619 for i=length(num_err) 620 f_glob_int(num_err(i))=0; 621 k_glob_int(num_err(i))=0; 622 end 623 for i=i:(length(f_glob_int)-length(num_err)) 624 if f_glob_int(i)==0 625 else 626 f_glob(j)=f_glob_int(i); 627 k_glob(j)=k_glob_int(i); 628 j=j+1; 629 end 630 end 631 end 632 633 634 %--------------------minimos cuadrados-------------------------------- 635 636 syms kt kss 637 error=0; 638 639 %se construye la función como un sumatorio en simbólico 640 641 for i=1:length(f_glob) 642 error=error+(k_glob(i)-(kt*kss*f_glob(i)/(kt+kss*f_glob(i))))^2; 643 end 644 645 erroraux=evalc('disp(error)'); %convierte la funcion en vector char 646 errorauxorig=convertCharsToStrings(erroraux); %covierte en string 647 erroraux=replace(errorauxorig,'kt','x(1)'); %sustituye las variables 648 erroraux=replace(erroraux,'kss','x(2)'); 649 straux="@(x)"; 650 straux2=[straux,erroraux];%mete el @(x) para que al evaluarlo sea 651 %funcion handle 652 erroraux=join(straux2); 653 error=eval(erroraux); %ya tenemos la funcion handle 654 options=optimset('MaxFunEvals',1e10,'MaxIter',1e3,'TolX',1e-7); 655 [X]=fminsearch(error,[1,1],options); %minimiza la función 656 657 Ignacio López Alonso Trabajo Fin de Grado 70 k_t=X(1); %rigidez del tendón 658 k_ss=X(2); %rigidez del sóleo 659 660 661 f_plot=1:0.01:2500; %f no suele ser mayor que 2500 (cambiar si es 662 necesario) 663 k_plot=k_t*k_ss.*f_plot./(k_t+k_ss.*f_plot); %gráfica analítica 664 665 666 rig_soleo=0; 667 i=1; 668 while rig_soleo < 1.1*k_t %gráfica rigidez sóleo (lineal) 669 rig_soleo(i)=k_ss*f_plot(i); 670 i=i+1; 671 end 672 673 %tabla de parámetros 674 675 table_mat=zeros(length(f_glob_plot),3); 676 677 table_mat(1:length(f_glob_plot),1)=f_glob_plot; 678 table_mat(1:length(f_glob_plot),2)=k_glob_plot; 679 table_mat(1:length(f_glob_plot),3)=R2_mat; 680 colnames={'f','k','R2'}; 681 682 fig=uifigure; 683 table1=uitable(fig,'Data',table_mat,'ColumnName',colnames); 684 685 %pinta verde si buena y rojo si mala 686 687 for i=1:length(f_glob_plot) 688 if R2_mat(i)>0.9 689 s = uistyle('BackgroundColor','green'); 690 addStyle(table1,s,'cell',[i,3]); 691 else 692 s = uistyle('BackgroundColor','red'); 693 addStyle(table1,s,'cell',[i,3]); 694 end 695 end 696 697 %dibuja gráfica rigidez 698 699 figure 700 plot(f_glob,k_glob,'.'),shg 701 hold on 702 plot(f_plot,k_plot,'.'),shg 703 hold on 704 plot(f_plot,k_t*ones(1,length(f_plot))),shg 705 hold on 706 plot(f_plot(1:length(rig_soleo)),rig_soleo),shg 707 xlabel('f (N)') 708 ylabel('k (kN/m)') 709 title('Resultados del estudio') 710 711 kt_ann=sprintf('k_{t}(kN/m)=%.4f',k_t); 712 kss_ann=sprintf('k_{ss}(kN/(m/kN))=%.4f',k_ss*1e3); 713 annotation('textbox',[.55 .4 .4 714 .5],'String',kt_ann,'FitBoxToText','on'); 715 annotation('textbox',[.2 .3 .4 716 .5],'String',kss_ann,'FitBoxToText','on'); 717 Ignacio López Alonso Trabajo Fin de Grado 71 718 %mensaje por pantalla 719 720 fprintf('\nSe ha ajustado la curva de rigidez con %d de los %d puntos 721 iniciales.',length(f_glob),length(f_glob_plot)); 722 723 724 end 725