Parametrización de un sistema de control para la estabilidad lateral de la bicicleta

Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster en Ingeniería Industrial Parametrización de un sistema de control para la estabilidad lateral de la bicicleta Autor: Francisco José Begines Gavira Tutor: Joaquín Ojeda Granja Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 II Trabajo Fin de Máster Máster en Ingeniería Industrial Parametrización de un sistema de control para la estabilidad lateral de la bicicleta Autor: Francisco José Begines Gavira Tutor: Joaquín Ojeda Granja Profesor titular Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 IV Proyecto Fin de Carrera: Parametrización de un sistema de control para la estabilidad lateral de la bicicleta Autor: Francisco José Begines Gavira Tutor: Joaquín Ojeda Granja El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2022 VI El Secretario del Tribunal A mi familia VIII Resumen El objetivo principal de este trabajo será diseñar un controlador PID que logre estabilizar la bicicleta para aquellas velocidades donde esta no lo sea. Para averiguar en qué momentos la bicicleta no es estable se realizará un breve análisis sobre la estabilidad del sistema. Una vez realizado el diseño del controlador se analizará el sentido físico de lso valores de este controlador, así como el comportamiento de la bicicleta para los casos simulados, que consistirán en simular para casos en los que la bicicleta tiene solamente una determinada velocidad inicial y otros casos donde se aplicará un pedaleo de manera que la velocidad sea constante. Para realizar el máximo de iteraciones posibles, se creará un fichero en Matlab que permita realizar un gran número de simulaciones y de esta manera encontrar rápidamente los parámetros idóneos del controlador. Como resultado de estas simulaciones se podrá observar a lo largo de este trabajo se deduce que para controlar la bicicleta para velocidades inferiores a 4,5 m/s y pequeñas perturbaciones, la primera acción a realizar será girar en la dirección de la perturbación, de manera opuesta a la acción que habrá que realizar en el caso en el que la velocidad de avance sea superior a 6.5 m/s, ya sea una pequeña perturbación o una perturbación grande, la acción que tendrá que realizar el ciclista para estabilizar la bicicleta será un giro en dirección opuesta. XVI NOTACIÓN A Matriz q vector AT Matriz traspuesta A-1 Matriz inversa 𝐩󰇗 Derivada primera 𝐩󰇘 Derivada segunda 𝜑 Phi 𝜃 Theta 𝜀 Épsilon 𝜉 Xi 𝛽 Beta 𝛾 Gamma 𝜓 Psi ∂y/ ∂x Derivada parcial de y respecto x 1 1 INTRODUCCIÓN omo bien se sabe, la bicicleta es uno de los vehículos más utilizados por las personas a la hora de desplazarse, ya sea para dar un paseo o para competir en alto nivel. Sin embargo, a su vez, la bicicleta es un sistema muy inestable provocando numerosos accidentes. Por ello se ve motivado este trabajo, en el que se tratará de analizar brevemente la inestabilidad de este vehículo utilizando como base el modelo de Whipple, un modelo simple que reduce la bicicleta a un sistema de cuatro sólidos asociados a unas determinadas restricciones, pero que permitirá analizar la inestabilidad con tan solo dos grados de libertad, esto también es posible gracias al método de partición de coordenadas, ya que dejará el sistema de ecuaciones diferenciales en función de las coordenadas asociadas a los grados de libertad. Una vez obtenido el estudio de la inestabilidad se podrá proceder a lo realmente interesante en este trabajo, el diseño de un controlador PID mediante la herramiento de Matlab, Simulink, en el cual se podrá experimentar con distintos ejemplos realistas en las cuales el sistema está en condiciones de inestabilidad y el controlador actuará estabilizando el sistema. El controlador modera el par aplicado en el manillar y variando este parámetro tendrá como objetivo mantener no solo mantener la bicicleta en pie, sino neutralizar en la medida de lo posible la perturbación inicial. Los casos en los que se va a analizar este comportamiento serán simulaciones en las cuales se aplicará una velocidad inicial y una perturbación en el manillar. Como últimos casos también se analizarán los casos en los que se introduce un par de pedaleo consiguiendo que la velocidad sea constante, ya que en los casos anteriores la velocidad decaía con el avance del tiempo. Con estas simulaciones se podrán obtener la evolución de los parámetros del controlador en función de las distintas velocidades y así de esta manera crear una forma de obtener los parámetros rápidamente. Siendo este el objetivo principal de este trabajo. C Introducción 2 3 2 DESCRIPCIÓN DEL MODELO DE BICICLETA n este apartado se explicará brevemente en qué consiste el modelo de la bicicleta de Whipple, así como las ecuaciones que rigen la dinámica de la bicicleta, tanto el desarrollo del modelo como las restricciones se pueden encontrar con más detalles en la bibliografía [1] y [2]. 2.1 Modelo de Whipple El modelo de Whipple se basa en dividir la bicicleta en un sistema multicuerpo, concrentamente en dividirla en 4 sólidos rígidos que consisten en la rueda delantera, rueda trasera, un conjunto de manillar y cuadro delantero y finalmente el cuadro de la bicicleta junto con el ciclista. Para incorporar al ciclista se añade una masa puntual a la altura del sillín. Esta descomposición se puede apreciar en la Figura 1. Figura 1. Esquema del modelo de Whipple. J. L. Escalona [1] A este modelo hay que asociarle un sistema de coordenadas, para utilizar el menor número de restricciones posible, se utiliza las minimas coordenadas posibles, quedando estas de la siguiente manera:  𝑥𝑐, 𝑦𝑐 definen la posición del punto de contacto (C) de la rueda trasera en el plano X Y  𝜑 es el ángulo de guiñada, es decir, el ángulo que forma el eje 𝑥𝑖1 del sistema intermedio de referencia con el eje X  𝜃 es el ángulo de vuelco, que se forma entre el eje 𝑧𝑖2 del segundo sistema intermedio con el eje 𝑧𝑖1 del primer sistema intermedio. Con estos dos ángulos se puede determinar la orientación del plano que contiene el cuadro de la bicicleta.  𝜓 es el ángulo de rodadura, que forma el eje 𝑧2 del sistema asociado a la rueda trasera con el eje 𝑧𝑖2 del segundo sistema intermedio  𝛽 es el ángulo que forma el eje 𝑧3 del sistema del marco con el eje 𝑧𝑖2 del segundo sistema E Descripción del modelo de bicicleta 4 4 intermedio.  𝛾 es el ángulo de dirección que forma el eje 𝑥4 del sistema de referencia del manillar con el eje 𝑥3 del sistema de referencia del marco.  𝜀, ángulo de rodadura, es el ángulo que forma el eje 𝑧5 del sistema de referencia asociado a la rueda delantera con el eje 𝑧4 del manillar. Una vez definido estos parámetros se define el siguiente vector de coordenadas, al que se va a denominar como q: 𝐪=[ 𝑥𝑐 𝑦𝑐 𝜑 𝜃 𝜓 𝛽 𝛾 𝜀] Estas coordenadas se pueden ver reflejadas en la imagen 2 donde se muestra el sistema de la bicicleta con estas coordenadas. Figura 2. Modelo de la bicicleta [1]. J.L. Escalona 5 5 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta 2.2 Dinámica del sistema 2.2.1 Restricciones de la bicicleta Para recrear el modelo de la bicicleta, además de definir los sólidos que componen el sistema, será necesario definir las restricciones que intervienen en la bicicleta. 2.2.1.1 Restricción de rodadura sin deslizamiento. Con esta restricción se asegura que la bicicleta mantenga en todo momento las ruedas rueden sin deslizar, para conseguir esto, se deberá forzar al sistema a que la velocidad en el punto de contacto de las ruedas con el suelo sea nula. 𝑉𝐶=0 𝑉𝐷=0 (2.1) En la ecuación (2.1) se observa como se define las velocidades del punto C, punto de contacto de la rueda trasera y punto D, punto de contacto de la rueda delantera, que se observan en la Figura 2. Estas ecuaciones proporcionan hasta 4 restricciones ya que cada rueda proporciona dos ecuaciones, en los ejes x e y. 2.2.1.2 Restricción de contacto Esta restricción fuerza a la rueda delantera a mantener en todo momento contacto con el suelo. Para definir esta restricción será necesario incluir una nueva coordenada, 𝜉. Donde 𝜉 es un parámetro, angular en este caso, como se observa en la Figura 2, que marca la orientación del radio que contiene al punto P con respecto al eje 𝑥5. Para forzar esta restricción habrá que obligar al sistema a que la posición de la rueda en el eje z sea nula, así como su tangente. [𝑟𝐷 (𝑝, 𝜉)]𝑧 = 0 [𝑡𝐷 (𝑝, 𝜉)]𝑧 = 0 (2.2) Proporcionando otras dos restricciones. Obteniendo un total de 6 restricciones para 9 coordenadas, por lo que se obtienen un total de 3 grados de libertad que consisten en el ángulo de vuelco (𝜃), ángulo de giro o dirección (𝛾) y el ángulo de rotación de una de las ruedas, para este caso será la rueda trasera (ψ). Estas coordenadas independientes se denominan como coordenadas dinámicas, y aquellas coordenadas de las que su valor depende de las coordenadas independientes se les denominará dependientes o cinemáticas, quedando dividido el vector q de la siguiente manera: 𝐪𝐢𝐧𝐝=[𝜃 𝜓 𝛾] 𝐪𝐜𝐢𝐧=[𝑥𝑐 𝑦𝑐 𝜑 𝜀] 𝐪𝐝𝐞𝐩=[𝛽𝜉] (2.3) Como se observa en (2.3) dentro de las coordenadas dependientes se ha subdividido en cinemáticas y dependientes, esto se debe a que la forma de obtener unas y otras son diferentes tal y como se aprecia en [2]. A las restricciones del sistema habrá que añadir una restricción adicional para poder analizar la estabilidad en función de la velocidad de avance y no de la velocidad de giro de la rueda, para ello se introduce la siguiente restricción: Descripción del modelo de bicicleta 6 6 𝜓−𝑉 𝑅𝑡=0 (2.4) Donde V es la velocidad de avance y R es el radio de la rueda trasera. 2.2.2 Dinámica del sistema Una vez definida las restricciones, habrá que definir las interacciones que sufre con el medio, como la acción del ciclista, gravedad, etc. 2.2.2.1 Fuerza generalizada de gravedad Para modelar las fuerzas de gravedad de cada sólido se utiliza las siguientes expresiones: 𝐅𝐠𝐫𝐚𝐯 𝐢=[0 0 − 𝑚𝑖𝑔]𝑇 (2.5) Siendo m la masa de cada sólido y g la gravedad. Quedando como expresión final: 𝐐𝑔𝑟𝑎𝑣=∑(𝐇𝑖)𝑇 5 𝑖=2 𝐅𝐠𝐫𝐚𝐯 𝐢 (2.6) Siendo H la matriz jacobiana de la velocidad de cada sólido. 2.2.2.2 Fuerza generalizada de acción del ciclista Las acciones que puede realizer el ciclista son el giro, pedaleo y vuelco, que se reflejan de la siguiente manera: 𝐖𝑎𝑐𝑐=𝐰𝟒     ∗𝑇∗ [0 0 𝑀𝑑𝑖𝑟]𝑇+𝐰𝟐     ∗𝑇∗ [0 𝑀𝑝𝑒𝑑 0]𝑇+𝐰𝟑     ∗𝑇∗[𝑀𝑣𝑢𝑒 0 0]𝑇 (2.7) Siendo 𝐰𝒊    ∗ la velocidad angular de cada sólido. 𝐐𝑎𝑐𝑐=[𝐆𝟒     ∗ (3,:)𝑇 𝐆𝟐     ∗ (2,:)𝑇 𝐆𝟑     ∗ (1,:)𝑇]∗[𝑀𝑑𝑖𝑟 𝑀𝑝𝑒𝑑 𝑀𝑣𝑢𝑒]𝑇 (2.8) Siendo 𝐆𝐢 ∗ la matriz jacobiana de las velocidades angulares. 2.2.2.3 Fuerza generalizada de resistencia aerodinámica La Resistencia aerodinámica se modela como la velocidad del sólido más voluminoso al cuadrado multiplicada por el coeficiente aerodinámico. |𝐹𝑎𝑒𝑟𝑜3|=𝑐𝑣∗(𝐯𝟑𝐆)𝑇∗𝐯𝟑𝐆 (2.9) Siendo 𝑐𝑣 el coeficiente aerodinámico, para este caso se ha utilizado un coeficiente de 0.5. 𝐐𝑎𝑒𝑟𝑜=−𝑐𝑣∗|𝐯𝟑𝐆|𝐇𝟑𝐓∗𝐯𝟑𝐆 (2.10) Siendo H3 la matriz jacobiana de la velocidad del sólido 3. Las acciones de resistencia aerodinámica y la de acción del ciclista componen la fuerza generalizada externa. 𝐐𝑒𝑥𝑡=𝐐𝑎𝑐𝑐+𝐐𝑎𝑒𝑟𝑜 (2.11) 7 7 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta Una vez modeladas todas las interacciones que intervienen en el modelo ya se puede resolver la ecuación dinámica de Newton-Euler: [𝐌𝐝𝟎 𝟎 𝐈𝐆][𝐚𝐆 𝛂 ]=[𝐅 𝐌 ]+[ 𝟎 −𝐰∧(𝐈𝐆 ∗𝐰)] (2.12) Siendo 𝐌𝐝 la matriz de masas, 𝐈𝐆 la matriz de inercias y el ultimo término es lo que se denomina como fuerza cuadrática en velocidad o 𝐐𝑣. Si se hacen ciertos cambios en los que la ecuación de Newton-Euler dependa únicamente de las coordenadas independientes, la ecuación general que presenta el sistema es la siguiente: 𝐄𝑇𝐌𝐄𝐩𝑖󰇘=𝐄𝑇(𝐐𝑣+𝐐𝑔𝑟𝑎𝑣+𝐐𝑒𝑥𝑡−𝐌𝐉) (2.13) Siendo M la matriz de masas de la bicicleta, E y J son matriz dependientes de las coordenadas que deja a las coordenadas dependientes en función de las independientes, 𝐩𝐢󰇘, el vector de aceleraciones de las coordenadas independientes y los vectores Q, que representan los vectores de las fuerzas que actúan en el sistema, como son el vector de velocidades cuadráticas, 𝐐𝑣, el vector de la acción de la gravedad, 𝑸𝒈𝒓𝒂𝒗 y el vector de las acciones externas, que para este análisis serán el pedaleo y giro del manillar y la acción de la resistencia aerodinámica. Esta ecuación se puede dejar de forma más simple si se reagrupan los términos de la siguiente manera: 𝐌𝑖(𝐩)𝐩𝑖󰇘=𝐐𝑖(𝐩,𝐩󰇗) Donde 𝐌𝑖=𝐄𝑇𝐌𝐄 y 𝐐𝐢(𝐩,𝐩󰇗)=𝐄𝐓(𝐐𝐯+𝐐𝐠𝐫𝐚𝐯+𝐐𝐞𝐱𝐭−𝐌𝐉) (2.14) Para obtener E y J hay que desglosar el vector 𝐩󰇗 𝐩󰇗=[𝐩𝐢󰇗 𝐩𝐝󰇗]=[𝐈 −𝐃𝐝−𝟏𝐃𝐢]𝐩𝐢󰇗+[𝟎 −𝐃𝐝−𝟏𝐂𝐭]=𝐄𝐩𝐢󰇗+𝐅 (2.15) Derivando esta expresión para obtener 𝐩 󰇘 𝐩󰇘=𝐄𝐩𝐢󰇘+𝐄𝐩𝐢󰇗+𝐅󰇗=𝐄𝐩𝐢󰇘+𝐉 (2.16) Con todo esto ya se podría proceder a realizar los cálculos dinámicos para cualquier bicicleta, para este caso se ha caracterizado a la bicicleta con los siguientes parámetros, con unidades en el Sistema Internacional: Diseño del controlador 14 14 Figura 7. Evolución temporal de los ángulos de dirección y vuelco para 5 m/s y trayectoria respectivamente En la Figura 7 se observa como ante la perturbación inicial de 0.1 rad, el ángulo del manillar sufre una pequeña oscilación antes de alcanzar la estabilidad aproximadamente a los 7 segundos. También se aprecian las similitudes que se encontraban para el caso inestable, en la cual la evolución temporal cambia de signo debido a la inercia del sistema. De hecho, se puede observa que el momento que el giro tiene un mínimo, el vuelco tiene un máximo, demostrando que en el punto de giro mínimo el vuelco es máximo en sentido contrario. Finalmente, en la Figura 7 se puede ver la trayectoria seguida en los 10 segundos de simulación, en este se contempla como después del giro inicial la bicicleta logra trazar una trayectoria recta habiendo neutralizado la perturbación inicial. 15 15 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta 5 OBTENCIÓN Y ANÁLISIS DE RESULTADOS n este capítulo se va a proceder a exponer la metodología que se ha seguido para obtener los resultados que hagan el sistema estable para los dos tramos en los que la bicicleta es inestable, así como el análisis de estos resultados. Además, se incluirá el caso en el que existe un par de pedaleo manteniendo la velocidad constante. 5.1 Metodología para obtener resultados. Normalmente a la hora de obtener los parámetros del controlador PID, se analiza lo que se denomina como función de transferencia del sistema, pero en este caso no se tiene eso, si no que se tiene el modelo matemático, por ello la forma más simple, aunque tediosa, es ir probando combinaciones hasta obtener los resultados que se deseen. Para empezar a saber sobre qué orden de magnitud tienen los parámetros del controlador se aconseja utilizar el método de Zieger Nichols, cuya idea se basa en ir aumentando el valor de Kp hasta que el sistema sobreoscile y reducir este valor hasta apróximadamente obtener el 50% de ese valor ahí. Este método solo se ha utilizado para conocer el orden de magnitud. Antes de empezar a simular hay que decidir que objetivos se tienen que cumplir para aceptar el resultado como óptimo, obviamente el óptimo tiene que haber alcanzado una estabilidad significativa y tener un valor final en torno al 1% del valor inicial de la perturbación, que es de 0.1 rad. Tabla 2. Ejemplo de metodología para obtener los coeficientes óptimos En la Tabla 2 se muestra un ejemplo de como conseguir el mejor resultado para la perturbación de 0.1 rad, en este caso en concreto a 6.5 m/s de velocidad inicial. En esta tabla se pueden apreciar hasta 5 tonalidades distintas. En rojo se muestran las celdas donde la bicicleta no es estable, ya sea porque tiene una sobreoscilación o incluso se cae. En blanco se muestran los valores en los que, aunque la bicicleta no es del todo inestable, los valores en los que se ha estabilizado son insifucientemente buenos. En naranja se puede apreciar que el valor si está dentro del margen que se ha establecido como óptimo, pero realmente no es un candidato válido debido a que si se hace la simulación se puede comprobar fácilmente que el valor final no es constante en torno a 0.70 % por lo que no se puede dar como válido. Por último, en color verde se muestra los valores que si son susceptibles de ser óptimos o al menos dar una buena respuesta ante la perturbación inicial. Con todos los candidatos obtenidos se pueden elegir diferentes factores que hagan una opción mejor que otra, como por ejemplo el tiempo que tarda el sistema en estabilizarse o que el valor máximo de la variable no alcanze un valor elevado. Sin embargo, no hay diferencia apreciable para elegir un candidato u otro en función de estas opciones expuestas, por tanto, se ha decidido optar por elegir como parámetros aquellos que consigan reducir al máximo la perturbación inicial. Este procedimiento como se puede llegar a ver es muy tedioso en caso de que no se conozcan los valores aproximados por los que se encuentra la estabilización óptima, de hecho, como se verá más adelante, para algunos casos resulta muy difícil encontrar candidatos para velocidades que se encuentran por debajo de los 4,5 m/s. Por lo tanto, resulta interesante crear un algoritmo o fichero en Matlab que pueda encontrar los valores KP -4 -3 -2 -1 0 1 2 3 -7 -2,17% 0,18% -0,15% 0,05% -0,64% -1,39% -6 8,21% -0,22% -0,20% -0,06% -0,15% -0,58% -5 0,70% -0,73% -0,50% -0,31% -0,18% -0,24% -4 -1,93% -2,00% -1,65% -1,14% -0,70% -0,47% -3 -4,00% -5,27% -5,07% -4,30% -3,60% -1,40% OBTENCIÓN PARÁMETROS ÓPTIMOS PARA V=6,5 M/S KD E Obtención y análisis de resultados 16 16 óptimos de manera automática sin tener que ir probando de forma aleatoria y de esta manera poder alcanzar un número de simulaciones muy superior al que se podría hacer de forma manual. Por ello se ha creado un fichero que realiza estas simulaciones para velocidades superiores a 6,5 m/s, es decir el segundo tramo inestable. Este algoritmo consiste en un método iterativo donde para una velocidad y perturbación determinada simula todas las combinaciones posibles de los parámetros del controlador, una vez simulada la combinación extrae la evolución de todas las coordenadas para poder analizar la que resulte más interesante. Además lo realmente interesante de este algoritmo es que mediante un condicionante se puede obtener una matriz que contenga todas las combinaciones susceptibles de ser óptimas. Para este filtro se ha impuesto que la respuesta en su instante final sea menor del 1% del valor inicial de la perturbación, así como que la diferencia absoluta entre el valor final y el valor alcanzado al 90% de la simulación sea menor de 0.2, evitando así que la respuesta sea una oscilación que por casualidad su valor final era menor del 1%. Aún así con este método cabe la posibilidad de tener una respuesta que oscile y cuyo periodo coincida en la respuesta final y en el 90%, dando así un resultado erróneo, pero esta probabilidad es tan pequeña para los casos simulados que se ha dado por bueno este algoritmo. 5.2 Resultados para velocidades superiores a 6.5 m/s Como se ha comentado, para este tramo de velocidades se puede obtener los resultados óptimos mediante el fichero de manera automática, se han realizado 6 simulaciones para distintas velocidades entre ellas 6.5, 7, 8, 9, 11 y 12, se ha puesto 12 m/s como límite de velocidad ya que esta velocidad es aproximadamente 43 km/h, practicamente una velocidad media que se puede alcanzar en una etapa de ciclismo de alto nivel. Además hay que recordar que hasta el momento no hay pedaleo en las simulaciones, es decir, se la aplica una velocidad inicial y se simula durante 10 segundos, aunque las restricciones dinámicas que se han impuesto no contemplan el rozamiento pero si la resistencia aerodinámica, y por tanto hay una pequeña caída de velocidad durante estos 10 segundos. De hecho, se puede apreciar en que si se vuelve a la Figura 3, para 6.5 m/s el sistema debería ser estable, pero debido a la pérdida de velocidad a lo largo de la simulación el sistema se vuelve inestable. En la Figura 8 se observa la evolución de la velocidad de la rueda trasera conforme avanza la simulación en rad/s, cabe recordar que el tramo para el cual la bicicleta es estable es de 15 a 22.6 rad/s, debido a que el radio de la rueda trasera se ha tomado como 30 cm. Como se puede apreciar en el tramo final de la Figura 8, la velocidad experimenta una subida, hasta un punto máximo donde la bicicleta se cae. Este aumento de velocidad se debe a que el ciclista empieza en volcarse provocando un incremento de la inercia y de este manera empieza a describir una curva en la dirección en la que el ciclista se tumba, el último instante donde vuelve a decrecer la velocidad, no se debe tener en cuenta, ya que la bicicleta está ya en proceso de caída . Aunque cabe pensar que la caída de velocidad para 6,5 m/s no debería influir mucho ya que esa caída tiene como punto más bajo 5.6 m/s y por lo tanto estaría dentro del rango estable, la continua caída de velocidad no permite a la bicicleta estabilizar adecuadamente provocando la caída en los últimos instantes de la simulación. Figura 8. Evolución temporal de la velocidad de la rueda trasera (rad/s) para una velocidad inicial de 6,5 m/s El objetivo de estas simulaciones será obtener una serie de puntos para las distintas velocidades de manera de 17 17 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta poder obtener una gráfica que pueda representar de forma aproximada la evolución de los parámetros Kp y Kd en función de la velocidad. De hecho como se habrá podido observar, en las velocidades elegidas para simular se ha obviado la velocidad de 10 m/s de forma intencionada, ya que una vez obtenidas las gráficas se va a proceder a testear los parámetros obtenidos para esa velocidad y comprobar si realmente son estables. Como resultado de las simulaciones se obtienen la siguientes gráfica: Figura 9. Evolución de los parámetros en función de la velocidad En la Figura 9 se observa que a medida que la bicicleta tiene mayor velocidad inicial, el sistema necesita un valor cada vez más negativo, es decir, cada vez se necesita más par en el manillar para poder estabilizar la bicicleta En cuanto al término derivativo, como se ya se dijo en el capítulo anterior, tiene como función suavizar la sobreoscilación que puede sufrir a causa de un valor elevado en el término proporcional, por tanto encontrar su valor exacto no es tan importante como el de el término proporcional, de manera que se tiene un margen mayor para cada velocidad, por ejemplo, si según la gráfica para 10 m/s el valor correspondiente para Kd es 4, el sistema podría ser perfectamente estable para valores entre 2 y 6 sin ningún tipo de problema, esto mismo se puede encontrar para el Kp para velocidades muy elevadas, en las cuales también se tiene un cierto margen de maniobra y se sigue estando en valores de estabilidad inferior al 1% de la perturbación inicial, como también se podía apreciar en la Tabla 2. 5.2.1 Resultados para una velocidad inicial de 10 m/s Como se ha comentado, el objetivo de este trabajo era obtener una forma de conocer que parámetros debe aplicar el controlador de la forma más rápida posible, con la Figura 9 ya se pueden obtener estos términos del PID, ahora solo falta verificar que efectivamente son válidos dentro del rango de velocidades que se han analizado, para ello se ha simulado a una velocidad inicial de 10 m/s que según la Figura 9, le pertenece un Kp de -25.65 y un Kd de 4. Para esta simulación se han obtenido los siguientes resultados: Obtención y análisis de resultados 18 18 Figura 10. Evolución temporal de los ángulos de giro y vuelco para 10 m/s y trayectoria respectivamente Como bien se observa en la Figura 10, se ha conseguido controlar el sistema multicuerpo para una velocidad inicial de 10 m/s y una perturbación de 0,1 rad, alcanzando un valor final del 0.18% del valor inicial, es posible que existan valores con un porcentaje un poco mejor, pero este resultado demuestra que se puede obtener unos parámetros que consiguen controlar la bicicleta con un gran resultado de manera fiable a partir de la Figura 9, además al tener un valor mucho menor que el 1 % permite tener un margen de error en los parámetros en los casos en los que el controlador no consiga llegar a los valores necesarios en tiempo real, en esto se hará más hincapié más adelante cuando se trate la simulación en tiempo real. En la evolución del vuelco se obtienen también los resultados equivalentes al caso estable del apartado anterior, en el cual se veía que además de conseguir estabilizar la bicicleta, las oscilaciones que se encuentran en los momentos iniciales coinciden en el mismo instante en los que existen oscilaciones en el ángulo de giro, pero de manera inversa, es decir, en los picos máximos del ángulo de giro se obtiene un valle en el de vuelco y viceversa. Por último, se observa que la trayectoria seguida también se estabiliza en los últimos metros de la simulación, en la que se observa que la distancia que se ha recorrido horinzontalmente es muy inferior, de manera proporcional, a la recorrida a 5 m/s, por lo que se está mejorando incluso el resultado de la autoestabilidad. Figura 11. Momento aplicado en el manillar para 10 m/s En la Figura 11 se puede ver el momento aplicado en la bicicleta, tal y como se observa llama la atención que la evolución de los valores del par aplicado coincide con la evolución del ángulo de giro del manillar de forma proporcional, salvo al inicio de la simulación. Este fenómeno se debe al bloque de saturación, que como ya se explicó, limita al controlador de forma que no le permite superar el valor de 3 Nm para evitar valores que hicieran 19 19 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta al sistema inestable durante simulación, estos valores al inicio de la simulación se pueden deber a la inestabilidad que provoca la integración numérica al inicio de la simulación. 5.2.2 Resultados para perturbaciones mayores Cuando se definió el problema inicial se estableció una perturbación inicial de 0.1 rad o 5.72 grados, pero como se dijo, el abanico de giro de un manillar consiste en ±30 grados aproximadamente, por ello se va a simular, con los mismos parámetros y una velocidad de 10 m/s, para comprobar si el sistema sigue siendo estable o por si el contrario es necesario ajustar los parámetros en función de la perturbación inicial. Figura 12. Evolución temporal de los ángulos de giro y vuelco para 0.5 rad de perturbación y trayectoria respectivamente En la Figura 12 se aprecia como a simple vista el sistema parece haberse estabilizado a la perfección, aunque el valor final es de -1.17% del valor inicial, lo que estaría muy ajustado al límite impuesto al principio de los objetivos, aunque siendo la perturbación más extrema a la que se puede enfrentar el sistema y estando muy cerca del 1%, se puede dar como válido. Como era de esperar, un aumento de la perturbación, provoca un aumento de la respuesta final. Es por ello por lo que a la hora de elegir los parámetros óptimos es preferible aquellos que consigan una respuesta lo más pequeña posible, en vez de conseguir una respuesta en torno a 1% pero que lo consiga de manera rápida, ya que en el caso de que para una perturbación pequeña se tenga un valor cercano al límite podría provocar que para perturbaciones de 0.5 rad, el sistema no sea estable. En cuanto al vuelco, lo único destacable es que tarda más en llegar a la estabilidad que en el caso de 0.1 rad. Por ultimo, evidentemente y tal y como se aprecia en la Figura 12, al tener mayor perturbación, el giro será mayor y por lo tanto la trayectoria describirá una curva más cerrada que en el caso de la perturbación pequeña, aún así se aprecia como en el último tramo describe una trayectoria recta. Obtención y análisis de resultados 20 20 Figura 13. Evolución temporal del par aplicado para 0.5 rad de perturbación Para el par aplicado se puede observar en la Figura 13 como llega más veces al valor límite, debido que al ser mayor la perturbación, mayor es el giro inicial que se tiene que aplicar. 5.2.3 Discusión de los resultados A la vista de los resultados obtenidos anteriormente se pueden concluir los siguientes hechos. A mayor velocidad se necesitará un parámetro proporcional mayor, en términos absolutos, así como para perturbaciones mayores aunque estas se pueden considerar estables dentro de los límites establecidos. La posibilidad de tener varias combinaciones que, aunque no sean las mejores, esten dentro de los límites, permite al controlador a no ser tan preciso y rápido a la hora de ir ajustando sus valores en función de los valores, al igual que permite obtener una gráfica de evolución de parámetros que con escasos puntos tenga un rendimiento más que aceptable. Por último, se puede observar en la evolución de gamma el comportamiento que debería tener el cicilista ante la perturbación, tanto para 0.1 rad como para 0.5 rad, se aprecia en las Figura 10 y Figura 12 como para estabilizar la perturbación lo primero que se debe hacer es girar en sentido contrario a esta perturbación y a partir de ahí estabilizar poco a poco oscilando ligeramente el manillar. 5.3 Resultados para velocidades inferiores a 4.5 m/s Como ya se comentó con anterioridad, para velocidades inferiores a 4,5 m/s los autovalores que se obtuvieron estaban bastante alejados de cero y por lo tanto eran mucho más inestables que los autovalores para velocidades mayores. Por tanto encontrar los parámetros que hagan el sistema estable será una tarea mucho más delicada que para velocidades mayores. Para el caso de velocidades pequeñas no se puede automatizar tan facilmente la obtención de parámetros, ya que en el funcionamiento de Simulink, una caída hace que el programa se quede parado en el momento exacto de la caída de la bicicleta, por lo que habrá que ir probando de manera manual cada valor. Tabla 3. Obtención de parámetros para 4 m/s Tal y como se aprecia en la Tabla 3, encontrar valores que no hagan al sistema inestable no parece una tarea KP -4 -3 -2 -1 0 1 2 3 0 1 2 3 4,00% 4 OBTENCIÓN PARÁMETROS ÓPTIMOS PARA V=4 M/S KD 21 21 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta sencilla, de hecho, la única combinación que no hace al sistema muy inestable no cumple con el límite impuesto para considerar que la bicicleta se ha estabilizado por completo, por lo que no es suficiente con estos parámetros, es por ello por lo que será necesario ampliar con más detalle las distitntas combinaciones. Tabla 4. Obtención de los parámetros óptimos para 4 m/s Como se aprecia en la Tabla 4, si se pueden encontrar valores adecuados para el sistema, esta tabla tiene la misma forma de representar las celdas que la Tabla 1, rojo para las casillas en las que el sistema no es estable, blancas para las celdas donde no se cumplen los requisitos establecidos, naranjas para las celdas donde, aunque su valor final es adecuado el sistema oscila un poco, por lo que no es un gran candidato, verde claro para las celdas que cumplen los requisitos pero no son óptimas y verde oscuro la celda con mejor resultado. Tanto en la Tabla 3 y Tabla 4 se refleja la fragilidad que tiene el sistema para bajas velocidades, esta fragilidad será cada vez mayor conforme se reduzca la velocidad inicial por lo que llegará un punto en el cual encontrar los términos adecuados sea una tarea muy complicada. Para tener una idea de como varían los parámetros en función de la velocidad se han simulado los casos en los que la velocidad inicial es de 3.5, 3.6, 3.7, 3.8, 3.9 y 4 m/s. Hay que recordar que no se está aplicando ningún par de pedaleo por lo que la velocidad final será inferior a la inicial y en casos en los que la velocidad es tan pequeña este hecho cobra bastante relevancia. Para el caso de una velocidad inicial de 4 m/s se han obtenido los siguientes resultados: Figura 14. Evolución temporal de los ángulos de giro y vuelco para 4 m/s y trayectoria respectivamente Para el ángulo de dirección, tal y como se puede apreciar en la Figura 14, se puede apreciar un ligero cambio con respecto a las gráficas de las velocidades mayores de 6.5 m/s, al principio de la simulación se observa como la tendencia del sistema es a seguir subiendo antes de bajar como hacia en los casos anteriores, por lo demás la evolución de la gráfica es bastante similar a las anteriores. En cuanto al vuelco, no existe esa diferencia inicial de la que se hablaba en la Figura 14, pero se observa como la subida que experimenta el ángulo de vuelco al inicio es menos vertical que en los casos de velocidades mayores, esto se puede deber al fenómeno que ocurre en los que el ángulo de giro sube al principio antes de decrecer. KD 2,5 2,6 2,7 2,8 2,9 3 3,1 3,2 -0,5 2,45% 0,58% -0,6 5,80% -1,50% 0,73% 0,67% -1,29% 1,57% -0,7 1,77% 1,50% 0,93% -0,18% -0,8 0,65% 0,20% -0,9 -1,44% -0,94% OBTENCIÓN PARÁMETROS ÓPTIMOS PARA V=4 M/S KP Obtención y análisis de resultados 22 22 Por último, en la trayectoria seguida por la bicicleta lo único interesante es que, para velocidades pequeñas, una pequeña perturbación provoca que se avance más en horizontal que en vertical, aunque el último tramo de la trayectoria sea recto, sería interesante encontrar una combinación de parámetros que lograran avanzar más en vertical que horizontal. 5.3.1 Evolución de parámetros del controlador Como ya se dijo se han hecho varias simulaciones para analizar la evolución de los parámetros en función de la velocidad, obteniendo los siguientes resultados. Figura 15. Evolución de los parámetros para velocidades menores a 4.5 m/s En la Figura 15 se aprecia como la tendencia es que a mayor velocidad el parámetro proporcional debe ser más negativo, pero a diferencia de la Figura 9, donde se apreciaba que el párametro proporcional para velocidades mayores a 6,5 era siempre negativo, aquí se observa justo lo contrario, para velocidades consideradas del tramo 1, el parámetro proporcional tiene que ser positivo. Del parámetro derivativo se obtiene la conclusión equivalente a la del parámetro proporcional, la tendencia es la misma, siendo todos los valores negativos y cuya evolución es ser cada vez mayor. Hay que tener en cuenta que la Figura 15 es de mera información para tener una idea de cómo se comporta los términos del controlador, pero en este caso, no es recomendable para obtener puntos intermedios, ya que como se pudo ver en las Tablas 3 y 4, la sensibilidad de los parámetros es extrema, aunque con estas gráficas se puede tener una idea de qué valores estarían próximos a los coeficientes óptimos. 5.3.2 Discusión de los resultados A la vista de lo anteriormente expuesto se pueden sacar las siguientes conclusiones. En el caso de velocidades menores a 4,5 m/s será necesario girar primero en la dirección de la perturbación para poder encauzar el rumbo adecuadamente, tal y como se apreciaba en la Figura 14, y de ahí el motivo de que el parámetro proporcional sea positivo, en cuanto al parámetro derivativo, ya se ha comentado varias veces que su valor va en concordancia con el parámetro proporcional y de signo contrario a este, y como la tendencia del parámetro proporcional es a bajar, la derivada de este resulta negativa y de ahí su signo. En cuanto al valor de la perturbación, como se ha podido observar, no se ha indicado nada de qué pasaría con estos parámetros para perturbaciones distintas y el motivo no es otro que al tener velocidades tan bajas, sin pedaleo, una perturbación mayor a 0.1 haría el sistema demasiado inestable para poder ser controlado, al menos con los valores que se han obtenido. 23 23 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta 5.4 Resultados para velocidad de pedaleo constante Hasta ahora el pedaleo de la bicicleta era inexistente, por lo que solo se aplicaba una velocidad inicial que iba decreciendo con el tiempo. En este apartado se aplicará un momento de manera que la velocidad sea constante y analizar la evolución de los parámetros del controlador. Para saber qué par hace constante la velocidad para una velocidad inicial se ha utilizado la Figura 16, que se puede encontrar en la bibliografía [7] Figura 16. Par de pedaleo (Nm) en función de la velocidad (m/s) para mantener la velocidad constante V.Hanne [7] Figura 17. Evolución de la velocidad de la rueda trasera para un par de pedaleo. En la Figura 17, se observa como efectivamente la velocidad se mantiene constante aunque al principio hay una breve oscilación, este ejemplo se ha realizado para una velocidad inicial de 7 m/s, que se traduce en 23.3 rad/s debido al diámetro de la rueda. Obteniendose las siguientes gráficas de los parámetros del controlador: Conclusiones 30 30 31 31 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta 7 REFERENCIAS [1] Escalona,José L., Recuero,Antonio. A bicycle model for education in multibody dynamics and real-time interactive simulation. [2] Begines Gavira, Fco José (2019), Análisis de la influencia del ciclista en el comportamiento dinámico de la bicicleta. Trabajo fin de grado. [3] Schwab, A.L., Meijaardb, J.P., & Kooijmana, J.D. (2012). Lateral dynamics of a bicycle with a passive rider model : stability and controllability. [4] Schwab. A. L , Meijaard, J.P (2013). A review on bicycle dynamics and rider control, Vehicle System Dynamics: International Journal of Vehicle Mechanics and Mobility, 51:7,1059-1090 [5] Schwab, Arend , Lange, P.D.L. , Happee, Riender ,Moore, Jason. (2013). Rider control identification in bicycling using lateral force perturbation tests. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. 227. 390- 406. 10.1177/1464419313492317. [6] Dialynas, G. (2020). An experimental approach into the quantification of steering and balance behaviour of bicyclists. https://doi.org/10.4233/uuid:3c6817fd-9d04-4461-9253-f02f0ca78a6a [7] Hanne, Victor. Modelado de las acciones de control de un ciclista sobre la bicicleta (2021). Trabajo Fin de Máster [8] Meijaard J.P, Papadopoulos Jim M, Ruina Andy and Schwab A.L. (2007). Linearized dynamics equations for the balance and steer of a bicycle: a benchmark and reviewProc. R. Soc. A.4631955–1982 [9] Troncoso García, Manuel. Estudio de la dinámica lateral de una bicicleta. Análisis de estabilidad. (2020). Trabajo Fin de Grado Referencias 32 32 33 33 Parametrización de un sistema de control para la estabilidad lateral de la bicicleta 8 ANEXO 8.1 Algoritmo para obtener los parámetros óptimos %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% SCRIPT PARA OBTENER LOS PARÁMETROS %%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global KD KP clear clc clear vectorsol i=1;j=1;k=1;l=1; for KP=-18:1:-15 for KD=2:2:10 sim("testcontrolgamma.slx") %EXTRAE TODAS LAS COORDENADAS PARA CADA SIMULACION coordenada1(:,:)=simout_vector.signals.values(:,1,end); coordenada2(:,:)=simout_vector.signals.values(:,2,end); coordenada3(:,:)=simout_vector.signals.values(:,3,end); coordenada4(:,:)=simout_vector.signals.values(:,4,end); coordenada5(:,:)=simout_vector.signals.values(:,5,end); coordenada6(:,:)=simout_vector.signals.values(:,6,end); coordenada7(:,:)=simout_vector.signals.values(:,7,end); coordenada8(:,:)=simout_vector.signals.values(:,8,end); coordenada9(:,:)=simout_vector.signals.values(:,9,end); coordenada10(:,:)=simout_vector.signals.values(:,10,end); momentos=simout_momentos.signals.values(:,1,end); valor=coordenada7(end)/0.1*100; %OBTIENE EL VALOR FINAL valor90=coordenada7(round(end*.9))/0.1*100; %OBTIENE EL VALOR AL 90% DE LA SIMULACIÓN vectorsol(i,1)=KP; vectorsol(i,2)=KD; vectorsol(i,3)=valor; i=i+1; vectorsol if valor<1 && abs(valor-valor90)<0.2 Kopt(k,1)=KP; %%% GUARDA EN LA PRIMERA COLUMNA KP Kopt(k,2)=KD; %%% GUARDA EN LA SEGUNDA COLUMNA KD Kopt(k,3)=valor; %%% GUARDA EN LA TERCERA COLUMNA EL VALOR FINAL k=k+1; end end end filename= 'Velocidad7autoPEDALEO.xlsx'; % CREA UN EXCEL CON TODAS LAS SIMULACIONES T= table(vectorsol(:,1),vectorsol(:,2),vectorsol(:,3)) ANEXO 34 34 writetable(T,filename,'Sheet',1,'Range','D1') clear T valormin=min(abs(Kopt(:,3))); %ENCUENTRA EL VALOR MINIMO EN Kopt lugar=find(abs(Kopt(:,3))==valormin);%ENCUENTRA EL LUGAR DEL VALOR MINIMO Kpopt=Kopt(lugar,1) %MUESTRA EL VALOR DE KP ÓPTIMO Kdopt=Kopt(lugar,2) %MUESTRA EL VALOR DE KD ÓPTIMO valorfinal=Kopt(lugar,3) %MUESTRA EL VALOR FINAL DE LA SIMULACIÓN