Optimización de dinámica inversa sin plataforma para biomecánica del ciclo de marcha y reacciones con el suelo

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Optimización de dinámica inversa sin plataforma para biomecánica del ciclo de marcha y reacciones con el suelo Autor: Juan José Rodríguez Perea Tutora: Juana María Mayo Núñez Departamento de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Sevilla, 2019 iii Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Optimización de dinámica inversa sin plataforma para biomecánica del ciclo de marcha y reacciones con el suelo Autor: Juan José Rodríguez Perea Tutora: Juana María Mayo Núñez Catedrática de Universidad Departamento de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 v Trabajo de Fin de Grado: Optimización de dinámica inversa sin plataforma para biomecánica del ciclo de marcha y reacciones con el suelo Autor: Juan José Rodríguez Perea Tutora: Juana María Mayo Núñez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal vii A mi familia, amigos y el profesorado de la escuela. ix Agradecimientos n primer lugar, dar la gracias a todo el personal docente de la escuela por todo lo que me han enseñado durante este tiempo. Especialmente agradecer a Juana Mayo su ayuda para guiarme en este trabajo. La paciencia y dedicación recibida en sus explicaciones para cada consulta me permitieron poder seguir avanzando en el estudio. Por último, agradecer a mi familia el apoyo durante estos años de estudio y a los compañeros de clase por la ayuda y apoyo en muchos momentos. Juan José Rodríguez Perea Sevilla, 2019 E xvii ÍNDICE DE TABLAS Tabla 1. Resultados de coeficientes para cada reacción 19 Tabla 2. Resultado coeficiente para doble apoyo 2 en apoyo simple modificado 26 xix ÍNDICE DE FIGURAS Figura 1. Sistema de coordenadas locales (tomado de [1]) 6 Figura 2. Sistema de coordenadas globales (tomado de [1]) 6 Figura 3. Procedimiento paso a paso para calcular reacciones en cuerpo multisegmento cuando el pie izquierdo está libre en [1] 7 Figura 4. Forma de las expresiones (7), gráfica derecha, y (8), gráfica izquierda 10 Figura 5. Resultados para distintos valores de coeficiente 13 Figura 6. Resultados dinámica inversa simplificada para tobillo derecho 14 Figura 7. Resultados dinámica inversa sin simplificar para tobillo derecho 15 Figura 8. Resultados dinámica inversa simplificada para tobillo izquierdo 15 Figura 9. Resultados dinámica inversa sin simplificar para tobillo izquierdo 16 Figura 10. Optimización en pie derecho en globales 17 Figura 11. Optimización en pie izquierdo en globales 17 Figura 12. Optimización en proyección de tobillo derecho en globales 18 Figura 13. Optimización en proyección de tobillo izquierdo en globales 18 Figura 14. Comparación Fx para proyección tobillo derecho e izquierdo 20 Figura 15. Comparación Fy para proyección tobillo derecho e izquierdo 20 Figura 16. Comparación Fz para proyección tobillo derecho e izquierdo 20 Figura 17. Comparación Mx para proyección tobillo derecho e izquierdo 21 Figura 18. Comparación My para proyección tobillo derecho e izquierdo 21 Figura 19. Comparación Mz para proyección tobillo derecho e izquierdo 21 Figura 20. Mz de tobillo derecho (gráfica izquierda) y My de tobillo izquierdo (gráfica derecha) 22 Figura 21. Fy tobillo derecho (gráfica izquierda) y Fy tobillo izquierdo (gráfica derecha) 22 Figura 22. Fz de tobillo derecho (gráfica izquierda) y Fz de tobillo izquierdo (gráfica derecha) 23 Figura 23. Fx de doble apoyo 2 en pie izquierdo 24 Figura 24. Fy de doble apoyo 2 en pie izquierdo 24 Figura 25. Fz de doble apoyo 2 en pie izquierdo 25 Figura 26. Mx de doble apoyo 2 en pie izquierdo 25 Figura 27. My de doble apoyo 2 en pie izquierdo 25 Figura 28. Mz de doble apoyo 2 en pie izquierdo 26 Figura 29. Resultados My para tobillo izquierdo en doble apoyo 2 27 Figura 30. División del pie multisegmento del modelo de Oxford 28 1 1 INTRODUCCIÓN a bioingeniería podría definirse como el uso de la tecnología junto con los conocimientos científicos biológicos con el objetivo de resolver problemas relacionados con la salud y los seres vivos. Una de las ramas de este campo es la denominada biomecánica, la cual estudia el movimiento de los seres vivos, sus causas y resolver los problemas a los que podrían verse sometidos. Existen aplicaciones en este campo desde hace bastante tiempo y mejoras y avances que han ido desarrollándose a través del estudio, como puede ser el ejemplo de la famosa pata de palo y las actuales prótesis ortopédicas. Esta es una disciplina cuyo avance y desarrollo es constante de modo que hoy en día su estudio y aplicación es bastante importante. El análisis de la dinámica del cuerpo es algo muy útil a la hora de estudiar distintos patrones de movimiento y cómo podría afectar esto al cuerpo humano. Más concretamente, un estudio de las reacciones o movimientos a la hora de correr o caminar nos podría permitir saber cómo poder corregir posibles daños debido a unos malos patrones o poder realizar mecanismos, como exoesqueletos y prótesis, que simulen ese movimiento para ayudar a personas con discapacidad física, por ejemplo. Inicialmente, necesitamos saber cómo es el movimiento del sujeto de estudio y la forma habitual para captar el movimiento del cuerpo es un sistema de marcadores con cámaras de infrarrojo y plataformas de fuerzas para recoger las resultantes del pie con el suelo. Sin embargo, este método puede encontrarse con algunos inconvenientes a la hora del uso de las plataformas ya que en algunos casos no es posible disponer de estas placas de fuerza para monitorizar las reacciones. También, aunque se disponga varias plataformas, pueden surgir casos en los que sea complicado su uso, por ejemplo, en personas con discapacidades motrices les podría resultar una complicada tarea coordinarse para pisar en las dos plataformas o también si el sujeto de estudio es un niño, quizás el uso de plataformas no sea lo más adecuado. Por estos motivos surge la necesidad de usar otro método para calcular las reacciones que no use plataforma. Este método, el cual va a ser empleado en este trabajo, es la dinámica inversa basada en la cinemática, con el cual en base a los datos obtenidos de los marcadores y cámaras de infrarrojo se podrán calcular las reacciones del pie sin necesidad de plataforma. Sin embargo, este método sigue contando con algún inconveniente, más concretamente en las fases de doble apoyo del ciclo de marcha. En fases de simple apoyo, cuando sólo existe un pie en contacto con el suelo, la dinámica inversa es capaz de calcular el problema y obtener las reacciones fácilmente, pero en las fases de doble apoyo nos encontramos que para resolver el problema tenemos más incógnitas que ecuaciones. Por lo tanto, necesitamos solventar este inconveniente para poder realizar el cálculo completo del ciclo de marcha con el método de la dinámica inversa basada en cinemática. A raíz de este problema, varios estudios han sido realizados para solucionar este problema, entre ellos nos centraremos en el estudio de Ren [2] y Villeger [3], pero a parte de estos otros muchos modelos han sido obtenidos por diferentes autores. El estudio realizado por Villeger [3] es el objeto de nuestro trabajo, en el que vamos a utilizar este método para resolver la dinámica inversa en el ciclo de marcha, y posteriormente comparar con los resultados obtenidos mediante el método de Ren [2]. En los siguientes apartados se desarrollará más detalladamente en qué consiste cada uno de estos estudios. En conclusión, la dinámica del cuerpo es algo que se ha estudiado bastante en muchas variantes y nosotros nos centraremos en el ciclo de marcha de una persona para calcular las reacciones en cada pie durante la pisada. Utilizaremos la dinámica inversa basada en cinemática para realizar nuestro estudio e intentaremos mejorar y optimizar los estudios anteriormente realizados para calcular las reacciones entre pie y suelo. L Introducción 2 2 1.1 Dinámica inversa La dinámica inversa nos permite conocer las fuerzas y momentos de reacción basándonos en los valores cinemáticos del cuerpo, así como en las fuerzas y momentos de inercia y masas. Pero el objetivo final de todo este modelo es poder calcular las interacciones entre las uniones de los elementos del cuerpo y sobre todo las fuerzas musculares que producen el movimiento. Generalmente para empezar, se procede a la captura del movimiento del sujeto de estudio a través de algún sistema óptico y la medición y obtención de las reacciones con el suelo mediante plataformas de fuerza. Utilizando un sistema de marcadores colocados en el cuerpo y el sistema óptico de cámaras infrarrojas podemos capturar dicho movimiento de las coordenadas que definen al sujeto en un modelo computacional. Con estos datos cinemáticos y de reacciones con el suelo obtenidos, se establece el problema dinámico o dinámica inversa y se calculan las fuerzas de reacción y momentos entre los diferentes elementos del cuerpo. Finalmente, mediante técnicas de optimización, las interacciones y fuerzas musculares son calculados y de este modo se obtiene un método para conocer las fuerzas motoras del movimiento inicialmente recogido, lo cual tiene varias aplicaciones como conocer la carga en los huesos o calcular el coste energético [4]. Cuando realizamos el análisis dinámico inverso podemos obtener las reacciones con el suelo para un movimiento dado cuando sólo un pie está apoyado sin necesidad de usar plataforma de fuerza, pero existe un obstáculo en las fases de doble apoyo ya que se convierte en un problema indeterminado. Generalmente, este problema se ha solucionado con el uso de plataformas de fuerza, como se indicó al principio de este apartado. De esta forma se supera el problema de resolver las incógnitas adicionales de las fases de más de un apoyo. A esto se le denomina dinámica inversa con plataforma, sin embargo, no siempre es posible disponer de placas de fuerzas para resolver este tipo de problemas y lo que se busca es un método virtual que nos permita calcular las reacciones sin la necesidad de plataformas, a lo que se le llama dinámica inversa sin plataformas. Por lo tanto, lo que se intenta buscar es un método que nos permita obtener el valor de las reacciones tanto en simple como en el doble apoyo sin hacer uso de las placas de fuerzas. En el estudio que vamos a realizar en este trabajo analizaremos uno de estos métodos de dinámica inversa sin plataforma para obtener las reacciones con el suelo teniendo como partida los datos cinemáticos y las masas. Un concepto interesante aparece en este escenario, ya que los valores de reacción calculados mediante este método y los obtenidos directamente mediante las placas de fuerzas no serán exactamente iguales. A la diferencia entre ambos se le denomina residuales dinámicos. De esta forma, el método que consiga obtener unos residuales dinámicos menores será el que se ajuste de mejor forma a los resultados experimentales reales. 1.2 Estado del arte Con el objetivo de obtener un método para resolver la dinámica inversa sin el uso de plataformas, diferentes estudios han sido realizados. Por ejemplo, Ren [2] introdujo el concepto de Suposición de Transición Suave. Este método se basa en la suposición de las componentes de fuerzas y momentos de acuerdo con una función que ajusta los valores a los que se obtendrían con las plataformas. De esta forma se supone las componentes del pie que sale de la zona de doble apoyo. Se utilizan dos fórmulas de transición una para la componente de la fuerza anteroposterior y la otra para las cinco restantes. De este modo, este modelo de suposición de la curva permite obtener los datos 3 adicionales que necesitamos en las fases de doble apoyo y así poder calcular la dinámica inversa. Sin embargo, este método está restringido para sujetos sanos sin problemas en los ciclos de marcha y con ciclos de doble apoyo cortos en comparación con el ciclo completo de pisada, cosa que no suele ocurrir ya que en la fase de doble apoyo puede llegar a ser tan grande como el ciclo de pisada [5]. En vez del uso de curvas de transición para calcular las reacciones, otros estudios realizan un modelo de contacto pie-suelo [6, 7, 8]. Se utiliza un modelo de fuerzas puntuales que proporcionan la fuerza de contacto en función de la hendidura entre las dos superficies de contacto. La fuerza de contacto total es dividida entre las componentes normales y tangenciales. Este proceso incluye incógnitas que son variable de diseño y la superficie del pie para contacto entre el pie y suelo está modelada por esferas. Para este caso, la solución consiste en encontrar los parámetros definidos de diseño, de modo que el error entre las componentes de reacción proporcionadas por la dinámica inversa y el modelo de contacto sea mínimo. Siguiendo la línea de los estudios con curvas de transición, el estudio de Villeger [3] propone mejoras al trabajo realizado en Ren [2]. Para empezar, utiliza únicamente una plataforma de fuerza en vez de dos para realizar el estudio de la dinámica inversa. Esto permite ampliar el rango de posibles casos a estudiar, ya que es más difícil disponer de dos plataformas y también al tener sólo una no es necesaria tanta coordinación del sujeto de estudio, ya que no tiene que realizar las pisadas en dos plataformas diferentes. Por este motivo, es posible realizar este estudio por ejemplo en el caso de niños pequeños o en personas con dificultades motrices. Sin embargo, el uso de una sola plataforma hace más compleja la labor de encontrar las fases de doble apoyo. Este problema se resuelve utilizando la velocidad de avance de centro de presiones. Cada vez que este parámetro toma el valor 0, se produce la transición de simple a doble apoyo y de doble a simple [9, 10]. El desarrollo de este estudio de Villeger [3] se centra en optimizar las curvas de transición propuestas en Ren [2] para que se ajusten mejor a los valores reales de la plataforma. En vez de dos funciones diferentes, sólo usa una para todas las componentes y el coeficiente que Ren define en su estudio como 2/3, en este trabajo es incógnita para cada componente y se trata de optimizar para conseguir la máxima aproximación a los resultados reales. Se estudian diversos casos y se realiza una regresión múltiple y al final se concluyen unas ecuaciones, una para cada componente de las recciones entre el pie y el suelo, en función de los valores anteriores a las fases de doble apoyo y teniendo en cuenta la velocidad de marcha para determinar las reacciones en ambos pies. Sin embargo, en nuestro trabajo actual no llegaremos hasta este paso, nos centraremos en el estudio de optimización y mejora del estudio realizado en Ren [2]. 1.3 Objetivos Este trabajo está basado en un estudio previo de dinámica inversa sin plataformas [1]. En dicho trabajo se desarrolla la dinámica inversa de un cuerpo dividido en segmentos con la finalidad de resolver las fuerzas de reacción del pie sin el uso de plataformas. Para resolver los problemas del doble apoyo se basa en el estudio de Ren [2] y todo esto se plasma en unas funciones implementadas en MatLab® para resolver este problema en base a unos datos cinemáticos ya recogidos anteriormente. El objetivo principal de nuestro estudio es por un lado implementar el mismo estudio [1] de una forma más eficiente y por otro lado mejorar el método utilizando el estudio de Villeger [3] para aproximar y ajustar de Introducción 4 4 una forma más exacta las fuerzas de reacción con el suelo calculadas para las fases del doble apoyo durante el ciclo de marcha cuando no se dispone de plataformas dinámicas para medirlas experimentalmente. La primera parte nos permitirá analizar el mismo método realizado de dos formas diferentes con el objetivo de obtener el mismo resultado. En el anterior estudio [1], el cuerpo humano, como objeto de estudio, se definía como un conjunto de segmentos. Sin embargo, nosotros tomaremos el cuerpo como un conjunto completo, por lo que simplificaríamos considerablemente el método en cuestión ya que no tendremos en cuenta la división en diversos segmentos. La segunda parte consistirá en una optimización del problema en la fase de doble apoyo basándonos en Villeger [3]. El anterior trabajo plantea la resolución del doble apoyo en base a Ren [2]. La optimización consistirá en encontrar unos coeficientes más adecuados para las funciones utilizadas en aproximar las curvas de las resultantes con el suelo, minimizando la diferencia entre los valores de las fuerzas de reacción reales, obtenidos de la plataforma, y las fuerzas de reacción que se calculen. La comparación entre los resultados obtenidos mediante los dos estudios y los resultados reales de la plataforma de fuerza nos permitirá conocer si la mejora ha sido realmente efectiva, en cuanto a aproximación a los resultados reales, o si por el contrario los resultados de ambos métodos son similares. 5 2 DESARROLLO En primer lugar, procederemos a realizar una mejora del código de la función implementada con MatLab® para el cálculo de la dinámica inversa, que consistirá en tomar el cuerpo como un sólido único que reacciona con el suelo en lugar de dividir el cuerpo en segmentos e ir calculando las fuerzas y reacciones en cada uno de ellos tal y como se hacía en el estudio anterior. Como sólo disponemos de los datos posición, matriz de masa, etc. de cada uno de los segmentos necesitamos primero hacer el sumario de todos ellos en cada instante de tiempo para poder calcular las fuerzas y momentos de inercia en el cuerpo completo. Cuando consigamos esto, lo único que tendremos que hacer es realizar ecuaciones de equilibrio entre el cuerpo y el suelo para calcular las reacciones. En segundo lugar, trataremos de optimizar los resultados obtenidos en las fases del doble apoyo para conseguir unas curvas más óptimas y de esa forma obtener unos cálculos más aproximados a los valores reales. La optimización consistirá en ajustar de una mejor forma las funciones usadas para suponer los valores de las reacciones en las fases de doble apoyo. Estas funciones dependen de un valor de coeficiente. En Ren [2], este valor es por defecto 2/3 para todas las curvas, sin embargo nosotros intentaremos encontrar el valor más óptimo para cada fuerza y momento, siguiendo el método descrito en Villeger [3], es decir, encontrando el valor que hace mínimo el cuadrado de la diferencia entre el valor real (obtenido de la plataforma) y el valor calculado (obtenido con la función de suposición). 2.1 Sistema de coordenadas A lo largo del trabajo se usarán dos sistemas de coordenadas de referencia, definidos a continuación, para el cálculo de las reacciones. Uno de ello es el sistema de coordenadas locales, teniendo como origen el centro de gravedad de cada pie o del elemento que queramos estudiar. Este sistema, al estar definido en un punto del cuerpo, variará en función del movimiento del cuerpo. Este sistema de coordenada estará definido de la siguiente forma: • El eje x perpendicular al plano frontal del pie, en dirección al sentido de la marcha • El eje y perpendicular al plano sagital del pie • El eje z perpendicular al plano transversal del pie Desarrollo 12 12 Siendo: 𝐺𝑅𝐹𝑖𝑧𝑞𝑢𝑖𝑒𝑟𝑑𝑎: la componente de reacción obtenida en la plataforma izquierda 𝐹2: la componente de reacción del pie izquierdo Hay que destacar que en el apoyo 0, tenemos que poner la componente de reacción del pie derecho en función del pie izquierdo, que es la que estamos suponiendo (ya que es el pie que sale de la fase de doble apoyo), o lo que es lo mismo, ponerlo todo en función del coeficiente de optimización que queremos calcular. Lo mismo ocurre en el apoyo 1, pero en este caso al contrario, poner los parámetros del pie izquierdo en función del derecho. También, en el apoyo 0 no estamos incluyendo los datos de la plataforma izquierda y del pie izquierdo porque básicamente no tenemos datos de la plataforma de fuerza en esa fase para el pie izquierdo. Lo mismo ocurre en el apoyo 2 para los datos del pie derecho. Una opción posible en este caso hubiera sido unificar la optimización del apoyo 0 y 2, ya que se supone que forman parte de la misma fase de apoyo, es decir, cuando entra en contacto el pie derecho y cuando sale del contacto el pie izquierdo. La optimización hubiera resultado algo parecido a la fórmula del apoyo 1. Sin embargo, esto no resultaría del todo válido porque no tenemos datos de la plataforma izquierda en el apoyo 0 ni datos de la derecha en el apoyo 2, y en consecuencia los datos en dicho apoyo no podrían solaparse entre ellos ya que, por ejemplo, reacciones supuestas del pie izquierdo en el apoyo 0 no son iguales a las supuestas del pie izquierdo para el apoyo 2 (pasa lo mismo con las calculadas del pie derecho). Esto traería como consecuencia una optimización bastante deficiente si lo hiciéramos de esta forma. Por otro lado, en el apoyo 2 podríamos haber optimizado las reacciones del pie izquierdo, que es el que estamos calculando y no suponiendo. De esta manera conseguiríamos unos resultados optimizados más ajustados para el pie izquierdo, pero más incorrectos en el pie derecho. Por este motivo, para conseguir una optimización equitativa para los dos pies, hemos decidido unificarlos en el apoyo 1 para hacer la optimización. Pues con lo definido anteriormente, ya podemos proceder a realizar el cálculo de la optimización. Una vez obtenido el coeficiente optimizado, 𝐶𝑜𝑒𝑓, para cada una de las componentes que necesitamos calcular, es hora de emplear la ecuación (8) para obtener los valores de las reacciones. 13 Para hacernos una idea de las formas de las gráficas que se pueden obtener con esta función dependiendo de los distintos valores del coeficiente podemos referirnos a la siguiente figura: Figura 5. Resultados para distintos valores de coeficiente En los siguientes apartados veremos los resultados obtenidos de emplear este método y los compararemos con los estudiados y obtenidos anteriormente por Ren. Resultados y discusión 14 14 3 RESULTADOS Y DISCUSIÓN Una vez que se ha terminado de explicar el objetivo y el desarrollo práctico de nuestro estudio, vamos a ver y comparar los resultados obtenidos de la realización de la simplificación y los que se obtenían anteriormente cuando el cuerpo estaba dividido por partes. De la misma forma los resultados de la optimización también serán analizados con los datos de las gráficas calculadas mediante Ren y con los datos experimentales obtenidos por las plataformas. 3.1 Resultados obtenidos de la simplificación En este apartado, desarrollaremos los resultados obtenidos de la simplificación realizada definiendo el cuerpo como un único sólido en vez de un cuerpo dividido en varios elementos. Como se desarrolló anteriormente, se ha realizado el mismo cálculo de dinámica inversa dividiendo el cuerpo en varios elementos y considerando el cuerpo como un sólido completo, y lo que buscamos es obtener los mismos resultados pero realizados de una manera más simple y comprensible. A continuación en las siguientes gráficas, se pueden observar los resultados obtenidos de las dos formas. En cada gráfica se muestran 6 figuras. Las tres de la parte superior se corresponden con las fuerzas de reacción con el suelo en coordenadas globales, Fx, Fy y Fz ordenadas de izquierda a derecha. Del mismo modo, en la parte inferior se dibujan los momentos de reacción con el suelo, Mx, My y Mz. También hay que comentar que las curvas dibujadas en línea continua se corresponden con los datos reales de plataforma y las curvas en línea discontinua de puntos a los resultados de la dinámica inversa. Figura 6. Resultados dinámica inversa simplificada para tobillo derecho 15 Figura 7. Resultados dinámica inversa sin simplificar para tobillo derecho Puede observarse que los resultados son idénticos para ambos casos, lo cual nos hace confirmar que la simplificación ha funcionado correctamente para el caso del pie derecho. Veamos ahora el caso para el pie contrario. Figura 8. Resultados dinámica inversa simplificada para tobillo izquierdo Resultados y discusión 16 16 Figura 9. Resultados dinámica inversa sin simplificar para tobillo izquierdo Como se puede comprobar, los resultados obtenidos por los dos métodos son equivalentes para ambos casos, por lo que se demuestra que la simplificación ha sido efectiva y el objetivo de eliminar cálculos innecesarios para obtener un método más simple con el mismo resultado se ha cumplido. Este ha sido el primer objetivo que realizar antes de pasar a la optimización de las fases de doble apoyo. Además, esta simplificación también ha permitido que todos los pasos posteriores de optimización se vuelvan mucho más sencillos que si hubiéramos tenido el caso del cuerpo dividido en segmentos. 3.2 Resultados obtenidos de la optimización En primer lugar, vamos a exponer los resultados obtenidos de realizar la dinámica inversa para obtener las reacciones, cogiendo de referencia el centro de cada pie y en coordenadas globales, y optimizando las fases de doble apoyo según los descrito en apartados anteriores. 17 Figura 10. Optimización en pie derecho en globales Figura 11. Optimización en pie izquierdo en globales Además de esto, también hemos calculado las reacciones en la proyección del centro del tobillo. Hemos elegido calcular también en esta referencia porque en anteriores estudios se había elegido este punto de estudio también y porque algunos resultados de las gráficas que se obtienen se ajustan más a los datos de la plataforma en el simple apoyo. Resultados y discusión 18 18 Figura 12. Optimización en proyección de tobillo derecho en globales Figura 13. Optimización en proyección de tobillo izquierdo en globales Como se ha podido observar, los resultados se ajustan bastante a los datos obtenidos por la plataforma. Principalmente, la gráfica de las fuerzas verticales, Fz, es la que más se ajusta a los datos reales. Por otro lado, las gráficas con más desviaciones son las de los momentos, concretamente la gráfica para My. En la siguiente tabla, se pueden ver los coeficientes obtenidos tras hacer la optimización para cada una de las reacciones, tanto de fuerzas como de momentos, para cada una de las fases de doble apoyo: 19 Apoyo Reacción 0 1 2 Fx -1.5546 0.1657 -0.0178 Fy 0.2937 0.4052 0.5063 Fz 0.1961 0.1116 0.0199 Mx 0.7969 0.2645 0.0036 My -26.6417 -0.0319 -7.2786 Mz -3.2147 0.3541 1.0456 Tabla 1. Resultados de coeficientes para cada reacción Como puede observarse, se obtienen resultados bastante dispares, tanto valores negativos como positivos, por ejemplo el valor -26.6417 para la componente de My en el apoyo 0 o el valor de 1.0456 para Mz en el apoyo 2. Sin embargo, con todos estos valores se consiguen resultados más optimizados que utilizando es coeficiente 2/3, como se define en Ren [2], para todas las reacciones. A continuación, se muestra algunos ejemplos de valores de la diferencia entre el valor experimental y en calculado de las reacciones (9) mediante el coeficiente de Ren y los obtenidos de la optimización. Para el caso de Fy en el apoyo 1, utilizando el coeficiente 2/3 se obtenía una diferencia un valor de 51063.94. Por otro lado, utilizando en valor del coeficiente optimizado, 0.4052 para este caso, se obtiene un valor de 14821.93. Para el caso de Mx en el apoyo 0 del pie derecho, utilizando el coeficiente 2/3 se obtenía una diferencia un valor de 3999.09. Por otro lado, utilizando en valor del coeficiente optimizado, 0.7969 para este caso, se obtiene un valor de 1710.79. Se ha reducido bastante la diferencia entre valor experimental y calculado, es decir, hemos ajustado mejor el conjunto de valores de las reacciones calculadas a los valores recogidos por las plataformas de fueras para las fases de doble apoyo. Estos datos numéricos de la diferencia pueden observarse gráficamente en el siguiente apartado, donde se muestra la comparación entre los cálculos usando el método de Ren y el optimizado, teniendo como referencia los valores experimentales de la plataforma. Resultados y discusión 20 20 3.3 Comparación de datos optimizados con no optimizados Ahora que tenemos los datos de los resultados optimizados, para comprobar si la optimización ha sido efectiva o no, vamos a comparar con los datos que no están optimizados. A continuación, se mostrarán las gráficas de comparación, donde se ha cogido como punto de referencia la proyección del tobillo en coordenadas globales. Las componentes del tobillo derecho se representan en las figuras de la izquierda y las de tobillo izquierdo en las figuras derechas. Figura 14. Comparación Fx para proyección tobillo derecho e izquierdo Figura 15. Comparación Fy para proyección tobillo derecho e izquierdo Figura 16. Comparación Fz para proyección tobillo derecho e izquierdo 21 Figura 17. Comparación Mx para proyección tobillo derecho e izquierdo Figura 18. Comparación My para proyección tobillo derecho e izquierdo Figura 19. Comparación Mz para proyección tobillo derecho e izquierdo Como podemos comprobar en los gráficos, las curvas de los datos optimizados (en azul) se ajustan más a datos de la plataforma (líne de puntos) en la mayoría de los casos que los datos no optimizados (en rojo). Hay que comentar que los tiempos en los que empiezan y acaban los doble apoyos son los siguiente: el doble apoyo 0 tiene inicio en el instante 11 y fin en el 31, el doble apoyo 1 tiene su inicio en el instante 68 y finaliza en el 86 y por último el doble apoyo 2 se inicia en el instante 120 y finaliza en el 143. Si nos fijamos más detalladamente en algunos resultados podremos apreciar y analizar de qué forma ha influido la optimización en los resultados. Para los casos por ejemplo de My en el pie derecho o Mz en el pie izquierdo, podemos apreciar como los resultados de la optimización disminuyen considerablemente la diferencia entre valor real y calculado. Conclusión 28 28 ajustan mejor a los valores reales, pero la diferencia no es demasiado significativa. Como nuevos objetivos, se podría estudiar la última parte del trabajo de Villeger [3], donde se define una serie de fórmulas para cada reacción en función de los datos de las reacciones previos al doble apoyo y de la velocidad de marcha. También se podría aplicar este estudio para el pie multisegmento definido en el modelo de Oxford, el cual divide el pie en tres segmentos diferentes, donde habrá más incógnitas que habrá que calcular y suponer algunas reacciones entre los diferentes elementos del pie multisegmento. Figura 30. División del pie multisegmento del modelo de Oxford Con el estudio realizado en este trabajo, podría aplicarse la optimización para obtener los valores de reacción entre los diferentes segmentos del pie multisegmento, donde sí se dispondría de datos de la plataforma y lo que habría que definir son los periodos en los que actúa cada una de las partes en el ciclo de marcha. Esto nos aportaría un estudio más exacto sobre las interacciones entre el pie y el suelo durante el ciclo de marcha y se obtendrían unos resultados más próximos a la realidad. 29 5 CÓDIGO DE MATLAB Función principal de dinámica inversa simplificada function [Fres cdg2 moment2 GRF2 Fres4 Fin Scoef ]=solveEquiWithoutGRF_simplificada (GRF,pos,vel1,acc1,AnimationData,MARKERS,datafile) %MAIN %solveID %solveIDForceEqui %clearvars -except GRF pos vel1 acc1 AnimationData MARKERS datafile MODEL close all clc global MODEL global AMTI off=1; %activate or desactivate Fin %%###################### [position,JointPos] = PlugInGait(MARKERS,datafile); %###################### Fres4 = zeros(size(GRF,1),6*MODEL.numBod); moment2=zeros(size(GRF,1),3); cdg2=zeros(size(GRF,1),3); %Descripción de variables para fuerzas de inercia sumInertialForces = zeros(1,3); Fin = zeros(size(GRF,1),12); Finf = zeros(size(GRF,1),12); Finm = zeros(size(GRF,1),12); Finmf = zeros(size(GRF,1),12); Fres_pie_der=zeros(size(GRF,1),3); Fres_pie_izq=zeros(size(GRF,1),6); %Calculo de Fext Fext =zeros(12,1);%Comprobar si el tamaño debe ser 12 o 6 bodymass=0; for k=1:15 bodymass=bodymass+MODEL.MassMatrix(6*k-5,6*k-5); end for k=1:2 Fext (6*k-3,:) = -bodymass*9.81;%modificar según el significado de MassMatrix end %F gravedad en cada segmento Fext2 =zeros(6*MODEL.numBod,1); for k=1:MODEL.numBod Fext2 (6*k-3,:) = -MODEL.MassMatrix(6*k-5,6*k-5)*9.81; end Fres = zeros(size(GRF,1),12);%Comprobar si el tamaño debe ser 12 o 6 Código de Matlab 30 30 p=zeros(15,4); %pJ=zeros(MODEL.numJoint,3); pC=zeros(15,3); %pJ2=zeros(MODEL.numJoint,3); pd=zeros(size(GRF,1),3); Rleft_rot=zeros(size(GRF,1),1); Rright_rot=zeros(size(GRF,1),1); %finding double support times %[th0, tt0, th, tt, th2, tt2] = DsupportTime(GRF); th0=11; tt0=31; th=68; tt=86; th2=120; tt2=143;%Valores buenos Juana h=0; %TO ENSURE THE ORDER for i=1:size(GRF,1) %plate Forces and moments on mass centre Rleft=Rot_Mat(pos(i,7*7-3:7*7)); Rright=Rot_Mat(pos(i,4:7)); vec3=(AMTI.pos(1:3)-pos(i,7*1-6:7*1-4)); %From the mass centre of the right foot to the centre of the right plate GRF2(i,1:3)=GRF(i,1:3); %GRF2(i,4:6)=GRF(i,4:6); %GRF2(i,4:6)=(Rright'*(GRF(i,4:6)+cross(vec3, GRF(i,1:3)))')'; GRF2(i,4:6)=((GRF(i,4:6)+cross(vec3, GRF(i,1:3)))')'; vec3=(AMTI.pos(4:6)-pos(i,7*7-6:7*7-4)); %From the mass centre of the left foot to the centre of the left plate GRF2(i,7:9)=GRF(i,7:9); %GRF2(i,10:12)=GRF(i,10:12); %GRF2(i,10:12)=(Rleft'*(GRF(i,10:12)+cross(vec3, GRF(i,7:9)))')'; GRF2(i,10:12)=((GRF(i,10:12)+cross(vec3, GRF(i,7:9)))')'; end for i=1:size(GRF,1) %Matrices de rotación para usarlas en la optimización Rleft=Rot_Mat(pos(i,7*7-3:7*7)); Rright=Rot_Mat(pos(i,4:7)); %Cálculo de las Fin (de las F y de los Momentos) for k=1:15 p(k,1:4) = pos(i,7*k-3:7*k); w(k,1:3) = vel1(i,6*k-2:6*k); pC(k,1:3)=pos(i,7*k-6:7*k-4); R = Rot_Mat(p(k,1:4)); % Inertial Forces Finf(i,1:3) =(MODEL.MassMatrix(6*k-5:6*k-3,6*k-5:6*k-3)* acc1(i,6*k- 5:6*k-3)')*off; %Fin son las -F de inercia Fin(i,1:3)=Fin(i,1:3)+Finf(i,1:3);%Fin pie derecho Fin(i,7:9)=Fin(i,1:3);%Fin pie izquierdo (= que pie derecho) % Fuerzas inercias momentos H = ([cross([1 0 0],R'*w(k,:)');cross([0 1 0],R'*w(k,:)');cross([0 0 1],R'*w(k,:)')]) * MODEL.MassMatrix(6*k-2:6*k,6*k-2:6*k) * R' * w(k,:)'; 31 Finm(i,1:3)=(MODEL.MassMatrix(6*k-2:6*k,6*k-2:6*k) * R' * acc1(i,6*k-2:6*k)'+H)*off; Finm(i,4:6)=R*Finm(i,1:3)'; Finm(i,10:12)=R*Finm(i,1:3)'; %Momento que ejerce las f+Fext2 sobre Pie Derecho pd(i,1:3)=pos(i,1:3)-pos(i,7*k-6:7*k-4); if norm(pd(i,1:3))~=0; Finmf(i,1:3)=cross(-pd(i,1:3),(Finf(i,1:3)-Fext2(6*k-5:6*k- 3)')); %Finmf(i,1:3)=(Rright'*(cross(pd(i,1:3),Finf(i,1:3)))')'; else Finmf(i,1:3)=0; end %Momento que ejerce las f+Fext2 sobre Pie Izquierdo pd(i,4:6)=pos(i,7*7-6:7*7-4)-pos(i,7*k-6:7*k-4); if norm(pd(i,4:6))~=0; Finmf(i,4:6)=cross(-pd(i,4:6),(Finf(i,1:3)-Fext2(6*k-5:6*k- 3)')); %Finmf(i,4:6)=(Rleft'*(cross(pd(i,4:6),Finf(i,1:3)))')'; else Finmf(i,4:6)=0; end %GRF2(i,10:12)=(Rleft'*(GRF(i,10:12)+cross(vec3, GRF(i,7:9)))')'; Fin(i,4:6)=Fin(i,4:6)+Finm(i,4:6)+Finmf(i,1:3);%Momentos pie derecho Fin(i,10:12)=Fin(i,10:12)+Finm(i,10:12)+Finmf(i,4:6);%Momentos pies izquierdo end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:size(GRF,1) Rleft=Rot_Mat(pos(i,7*7-3:7*7)); Rright=Rot_Mat(pos(i,4:7)); if (i>=1 && i<th0) || (i>tt && i<th2) %funcion que me devuelva que la pierna derecha está libre %Final forces Fres(i,7:9)=Fin(i,7:9)-Fext(7:9)'; Fres(i,10:12)=(Fin(i,10:12)); elseif i>tt0 && i<th %función que me devuelva que la pierna izq está libre Código de Matlab 32 32 %Final forces in Joint Fres(i,1:3)=Fin(i,1:3)-Fext(1:3)'; Fres(i,4:6)=(Fin(i,4:6)); %########################### Double support 1 ############################# elseif (i>=th && i<=tt) %Suponemos las Fres del pie derecho if h==1 Fx0(1,1:3)=Fin(i,1:3)-Fext(1:3)'; Fx0(1,4:6)=(Fin(i,4:6)); h=2; ds=1; Scoef1=Opt_Scoef_2(Fx0,Fin,Fext,GRF2,th,tt,ds,pos); end Fres(i,1:6)=DSUPPORT_prueba(Fx0,th,tt,i); Fres_prueba1(i,1:6)=DSUPPORT_prueba_Villeger(Fx0,th,tt,i,Scoef1); %Calculamosla Fres de pie izq vec1=pos(i,1:3)-pos(i,7*7-6:7*7-4); Fres(i,7:9)=Fin(i,7:9)-Fext(7:9)'-Fres(i,1:3); %Fres(i,10:12)=(Fin(i,10:12)-Fres_pie_der(i,1:3)); Fres(i,10:12)=(Fin(i,10:12)-Fres(i,4:6)- ((cross(vec1,Fres(i,1:3))))); %prueba optimizada Fres_prueba1(i,7:9)=Fin(i,7:9)-Fext(7:9)'-Fres_prueba1(i,1:3); Fres_prueba1(i,10:12)=(Fin(i,10:12)-Fres_prueba1(i,4:6)- ((cross(vec1,Fres_prueba1(i,1:3))))); %################## double support0 ###################### elseif (i>=th0 && i<=tt0) %Suponemos la Fres en el pie izquierdo if h==0 Fx0(1,1:3)=Fin(i,7:9)-Fext(7:9)'; Fx0(1,4:6)=(Fin(i,10:12)); h=1; ds=0; Scoef0=Opt_Scoef_2(Fx0,Fin,Fext,GRF2,th0,tt0,ds,pos); end 33 Fres(i,7:12)=DSUPPORT_prueba(Fx0,th0,tt0,i); Fres_prueba0(i,7:12)=DSUPPORT_prueba_Villeger(Fx0,th0,tt0,i,Scoef0); %Final forces pie derecho vec0=pos(i,7*7-6:7*7-4)-pos(i,1:3); Fres(i,1:3)=Fin(i,1:3)-Fext(1:3)'-Fres(i,7:9); %i %Fin %Fext(2) %Fres(i,8) %Fres(i,4:6)=(Fin(i,4:6)-Fres_pie_izq(i,1:3)); Fres(i,4:6)=(Fin(i,4:6)-Fres(i,10:12)- ((cross(vec0,Fres(i,7:9))))); %prueba optimizada Fres_prueba0(i,1:3)=Fin(i,1:3)-Fext(1:3)'-Fres_prueba0(i,7:9); Fres_prueba0(i,4:6)=(Fin(i,4:6)-Fres_prueba0(i,10:12)- ((cross(vec0,Fres_prueba0(i,7:9))))); %################## double support 2 ###################### elseif (i>=th2 && i<=tt2) %Suponemos la Fres en pie izq if h==2 Fx0(1,1:3)=Fin(i,7:9)-Fext(7:9)'; Fx0(1,4:6)=(Fin(i,10:12)); h=3; ds=2; Scoef2=Opt_Scoef_2(Fx0,Fin,Fext,GRF2,th2,tt2,ds,pos); end Fres(i,7:12)=DSUPPORT_prueba(Fx0,th2,tt2,i); Fres_prueba2(i,7:12)=DSUPPORT_prueba_Villeger(Fx0,th2,tt2,i,Scoef2); %Final forces pie derecho vec2=pos(i,7*7-6:7*7-4)-pos(i,1:3); Fres(i,1:3)=Fin(i,1:3)-Fext(1:3)'-Fres(i,7:9); Fres(i,4:6)=(Fin(i,4:6)-Fres(i,10:12)- ((cross(vec2,Fres(i,7:9))))); %prueba optimizada Fres_prueba2(i,1:3)=Fin(i,1:3)-Fext(1:3)'-Fres_prueba2(i,7:9); Fres_prueba2(i,4:6)=(Fin(i,4:6)-Fres_prueba2(i,10:12)- ((cross(vec2,Fres_prueba2(i,7:9))))); end %################ fin doble apoyo ###################### Código de Matlab 34 34 end %F de la función optimizada Fres_test=Fres; Fres_test(th0:tt0,:)=Fres_prueba0(th0:tt0,:); Fres_test(th:tt,:)=Fres_prueba1(th:tt,:); Fres_test(th2:tt2,:)=Fres_prueba2(th2:tt2,:); for i=1:size(GRF,1) Rright = Rot_Mat( pos(i,4:7)); Rleft = Rot_Mat(pos(i,46:49)); vecc=AnimationData.MODEL.JOINT.RAJC.childVec; vecci=AnimationData.MODEL.JOINT.LAJC.childVec; vec6d=pos(i,7*1-6:7*1- 4)+(Rright*vecc)'; vec6d(1,3)=0; vec6i=pos(i,7*7-6:7*7- 4)+(Rleft*vecci)'; vec6i(1,3)=0; vec6r=(AMTI.pos(1:3)-vec6d); vec6l=(AMTI.pos(4:6)-vec6i); GRF6(i,4:6)=GRF(i,4:6)+cross(vec6r, GRF(i,1:3)); GRF6(i,10:12)=GRF(i,10:12)+cross(vec6l, GRF(i,7:9)); GRF6(i,1:3)=GRF(i,1:3); GRF6(i,7:9)=GRF(i,7:9); end vecc=-AnimationData.MODEL.JOINT.RAJC.childVec; vecci=-AnimationData.MODEL.JOINT.LAJC.childVec; for ii=1:size(GRF,1) Rderecha = Rot_Mat( pos(ii,4:7)); Rizquierda = Rot_Mat(pos(ii,46:49)); %ponemos a 0 la coordenada z vectord=Rderecha*vecc; vectord(3)=0; vectori=Rizquierda*vecci; vectori(3)=0; %fres6(ii,4:6)=(Rderecha*(Fres(ii,4:6))')'+cross(vectord, Fres(ii,1:3)); %fres6(ii,10:12)=(Rizquierda*(Fres(ii,10:12))')'+cross(vectori, Fres(ii,7:9)); fres6(ii,4:6)=((Fres(ii,4:6))')'+cross(vectord, Fres(ii,1:3)); fres6(ii,10:12)=((Fres(ii,10:12))')'+cross(vectori, Fres(ii,7:9)); fres6_opt(ii,4:6)=((Fres_test(ii,4:6))')'+cross(vectord, Fres_test(ii,1:3)); fres6_opt(ii,10:12)=((Fres_test(ii,10:12))')'+cross(vectori, Fres_test(ii,7:9)); end fres6(:,1:3)=Fres(:,1:3); fres6(:,7:9)=Fres(:,7:9); fres6_opt(:,1:3)=Fres_test(:,1:3); fres6_opt(:,7:9)=Fres_test(:,7:9); 35 %Para quitar la zona del DS al principio en pie izq y al final en pie der %%Fres(100:150,1:6)=0; %Fres(1:50,7:12)=0; Scoef=[Scoef0,Scoef1,Scoef2]; %Scoef=[Scoef0,Scoef2]; %fres2(1:100,1:6)=Fres(1:100,1:6); %fres2(1:100,7:12)=Fres(1:100,7:12); %multiplot (Fres,GRF2) %multiplot (fres6,GRF6) %figure(714) %plot(Fres(:,2)) %hold on %figure(300) %multiplot(Fres_test,GRF2) multiplot(fres6_opt,GRF6) %W=size(Fres_test) %Y=size(GRF2) figure(200) %plot(Fres(:,10),'r') plot(fres6(:,3),'r') hold on %plot(Fres_test(:,10),'b') plot(fres6_opt(:,3),'b') hold on %plot(GRF2(:,5),'m') plot(GRF6(:,3),'m:') legend('reacciones sin optimizar','reacciones optimizadas','datos de plataforma') hold off figure(201) %plot(Fres(:,11),'r') plot(fres6(:,9),'r') hold on %plot(Fres_test(:,11),'b') plot(fres6_opt(:,9),'b') hold on %plot(GRF2(:,11),'m') plot(GRF6(:,9),'m:') legend('reacciones sin optimizar','reacciones optimizadas','datos de plataforma') hold off Código de Matlab 36 36 Función optimizar mediante fminsearch function Scoef= Opt_Scoef_2(Fx0,Fin,Fext,GRF,th,tt,ds,pos) %Scoef=zeros(6,1); x0=2/3; int=tt-th+1; GRFder(1:int,1:6)=GRF(th:tt,1:6); Finder(1:int,1:6)=Fin(th:tt,1:6); GRFizq(1:int,1:6)=GRF(th:tt,7:12); Finizq(1:int,1:6)=Fin(th:tt,7:12); options=optimset('Display','Iter','PlotFcns',@optimplotfval); if ds==0 %for i=1:3 %Scoef(i,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,i),Finder(:,i),Finizq(:,i),Fext(i,1),GRFder(:,i), GRFizq(:,i),th,tt,ds),x0); %Scoef(i+3,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento(Scoef,Fx0(1,i+3),Finder(:,i+3),Finizq(:,i+3),GRFder(: ,i+3),GRFizq(:,i+3),th,tt,ds),x0); %end Scoef(1,1)=fminsearch(@(Scoef1) Optimizar_Scoef(Scoef1,Fx0(1,1),Finder(:,1),Finizq(:,1),Fext(1,1),GRFder(:,1) ,GRFizq(:,1),th,tt,ds),2/3); AA1=Scoef(1,1); [AA,fval0,exitflag0,output0]=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,2),Finder(:,2),Finizq(:,2),Fext(2,1),GRFder(:,2), GRFizq(:,2),th,tt,ds),x0,options); Scoef(3,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,3),Finder(:,3),Finizq(:,3),Fext(3,1),GRFder(:,3), GRFizq(:,3),th,tt,ds),x0); AA2=Scoef(3,1); [XY,fval0,exitflag0,output0]=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0,Finder(:,4:6),Finizq(:,4:6),GRFder(:,4:6) ,GRFizq(:,4:6),th,tt,ds,1,AA1,AA,AA2,pos),x0,options); [XX,fval0,exitflag0,output0]=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0,Finder(:,4:6),Finizq(:,4:6),GRFder(:,4:6) ,GRFizq(:,4:6),th,tt,ds,2,AA1,AA,AA2,pos),x0,options); Scoef(6,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0,Finder(:,4:6),Finizq(:,4:6),GRFder(:,4:6) ,GRFizq(:,4:6),th,tt,ds,3,AA1,AA,AA2,pos),x0); Scoef(4,1)=XY; Scoef(2,1)=AA; Scoef(5,1)=XX; elseif ds==1 %for i=1:3 %Scoef(i,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,i),Finder(:,i),Finizq(:,i),Fext(i,1),GRFder(:,i), GRFizq(:,i),th,tt,ds),x0); %Scoef(i+3,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento(Scoef,Fx0(1,i+3),Finder(:,i+3),Finizq(:,i+3),GRFder(: ,i+3),GRFizq(:,i+3),th,tt,ds),x0); %end Scoef(1,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,1),Finder(:,1),Finizq(:,1),Fext(1,1),GRFder(:,1), GRFizq(:,1),th,tt,ds),x0); AB1=Scoef(1,1); 37 [AB,fval,exitflag,output]=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,2),Finder(:,2),Finizq(:,2),Fext(2,1),GRFder(:,2), GRFizq(:,2),th,tt,ds),x0,options); Scoef(3,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,3),Finder(:,3),Finizq(:,3),Fext(3,1),GRFder(:,3), GRFizq(:,3),th,tt,ds),x0); AB2=Scoef(3,1); Scoef(4,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0,Finder(:,4:6),Finizq(:,4:6),GRFder(:,4:6) ,GRFizq(:,4:6),th,tt,ds,1,AB1,AB,AB2,pos),x0); Scoef(5,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0,Finder(:,4:6),Finizq(:,4:6),GRFder(:,4:6) ,GRFizq(:,4:6),th,tt,ds,2,AB1,AB,AB2,pos),x0); Scoef(6,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0,Finder(:,4:6),Finizq(:,4:6),GRFder(:,4:6) ,GRFizq(:,4:6),th,tt,ds,3,AB1,AB,AB2,pos),x0); Scoef(2,1)=AB; elseif ds==2 %for i=1:3 %Scoef(i,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,i),Finder(:,i),Finizq(:,i),Fext(i,1),GRFder(:,i), GRFizq(:,i),th,tt,ds),x0); %Scoef(i+3,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento(Scoef,Fx0(1,i+3),Finder(:,i+3),Finizq(:,i+3),GRFder(: ,i+3),GRFizq(:,i+3),th,tt,ds),x0); %end Scoef(1,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,1),Finder(:,1),Finizq(:,1),Fext(1,1),GRFder(:,1), GRFizq(:,1),th,tt,ds),x0); [AC,fval,exitflag,output]=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,2),Finder(:,2),Finizq(:,2),Fext(2,1),GRFder(:,2), GRFizq(:,2),th,tt,ds),x0,options); Scoef(3,1)=fminsearch(@(Scoef) Optimizar_Scoef(Scoef,Fx0(1,3),Finder(:,3),Finizq(:,3),Fext(3,1),GRFder(:,3), GRFizq(:,3),th,tt,ds),x0); Scoef(4,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0(1,4),Finder(:,4),Finizq(:,4),GRFder(:,4), GRFizq(:,4),th,tt,ds,0,0,0,0,pos),x0); Scoef(5,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0(1,5),Finder(:,5),Finizq(:,5),GRFder(:,5), GRFizq(:,5),th,tt,ds,0,0,0,0,pos),x0); Scoef(6,1)=fminsearch(@(Scoef) Optimizar_Scoef_Momento_1(Scoef,Fx0(1,6),Finder(:,6),Finizq(:,6),GRFder(:,6), GRFizq(:,6),th,tt,ds,0,0,0,0,pos),x0); Scoef(2,1)=AC; end end