Modelado y control del péndulo invertido

En este post se abordará el problema clásico de control del péndulo invertido. Dicho problema consiste en un cuerpo pivotante alrededor de una base cuyo centro de masa se encuentra en el extremo contrario a la base, haciendo que el sistema se inestabilice a no ser que se le dote de un sistema de control que se capaz de desplazar la base hacia delante y hacia atrás. Por lo tanto, tendremos dos cuerpos, el carro y el péndulo. 

El sistema de control del carro debe ser tal que el carro se mueva para equilibrar el péndulo y que éste se mantenga en posición vertical a pesar de posibles perturbaciones. Este tipo de sistemas se implementan en diversas aplicaciones, desde un Segway al sistema de control del ángulo de inclinación de cohetes o misiles.

En este caso se considerará el problema bidimensional donde  la entrada es la fuerza de control u(t) que se aplica al carro y las salidas son la  posición angular del péndulo Φ(t) y la posición del carro x(t). Además se considerará que existe una fricción viscosa que se opone al movimiento del carro. 

Se desea diseñar un controlador para que, partiendo de posición vertical, ante un golpe brusco en el carro, restablezca la posición del péndulo de manera que éste vuelva a la posición vertical en menos de 5 segundos y de manera que el péndulo nunca se incline más de 0.05 radianes. Es decir que la respuesta de Φ(t) ante un impulso unitario δ(t) presente las siguientes características:

-          Tiempo de asentamiento < 5s.

-          Rebose ≤ 0.05 rad.

Los parámetros del sistema son los siguientes:

M (masa del carro) = 0.5 kg

m (masa del péndulo) = 0.2 kg

b (coef de fricción del carro) = 0.1 N/m/sec

l (distancia al centro de masa del péndulo) = 0.3 m

I (momento de inercia del péndulo) = 0.006 kg.m^2

Notar que siempre que se trabaja con FT, se están considerando sistemas SISO, por lo que no se han establecido requerimientos para la posición del carro. No obstante es interesante ver el efecto del control del péndulo sobre la posición del carro y valorar cómo se podría mejorar. 

Primeramente siempre se ha de observar la respuesta en LA para ver cómo de alejada está su dinámica respecto de las especificaciones deseadas. 

 1. Hallar la salida del péndulo en LA a un impulso unitario 

Lo primero que hay que hacer es modelar el sistema, obteniendo la FT de cada una de las salidas Φ(s) y X(s) ante la entrada U(s).

                           

Para ello se representan todas las fuerzas sobre el diagrama de cuerpo libre, donde N y P representan las componentes horizontal y vertical, respectivamente, de la fuerza de acción y reacción que el carro ejerce sobre el péndulo y viceversa, que habrá que calcular.

Además y puesto que la T. Laplace solo se aplica a sistemas LTI, y dado que el ángulo Φ siempre se va a mantener con valores muy pequeños, se aproximará el sistema considerando que cos Φ ≈ 1, sen Φ ≈ Φ y  ≈ ≈ 0 al objeto de linealizarlo.

Ecuación de fuerzas en el eje X del carro:                                                   

                                                 $$M\ddot { x } +b\dot { x } +N=u$$                                           (Ec.1)

Ecuaciones de fuerzas y ecuación de momentos del péndulo:

$$\sum { M } \rightarrow Pl\theta +Nl=I\ddot { \theta  } $$                    (Ec.2)

$$\sum { { F }_{ x } } \rightarrow N=m\left( \ddot { x } -lI\ddot { \theta  }  \right) $$          (Ec.3)

$$\sum { { F }_{ y } } \rightarrow P=mg$$                                     (Ec.4)

Por si aplicamos Laplace y combinamos las ecuaciones (Ec.1) y (Ec.3) obtenemos la siguiente expresión:

$$U(s)=X(s)[(M+m)s^{ 2 }+bs]-mlIs^{ 2 }θ(s)$$                              (Ec.5)

A continuación si aplicamos Laplace y combinamos las expresiones (Ec.2), (Ec.3) y (Ec.4) obtendremos la relación entre las dos salidas:

$$lms^2 X(s)=θ(s)[(I+l^2 m)-mgl]  $$                                        (Ec.6)

Si despejamos X(s) de la expresión (Ec.6):

$$X(s)=\frac { (I+l^{ 2 }m)-mgl }{ lms^{ 2 } } θ(s)$$                        (Ec.7)

Si despejamos  de la expresión (Ec.6):

$$θ(s)=\frac { lms^{ 2 } }{ (I+l^{ 2 }m)-mgl } X(s)$$                        (Ec.8)

Por lo tanto si combinamos la ecuación (Ec.6) con la (Ec.5) obtendremos la entrada respecto del ángulo de inclinación del carro:

$$q=(M+m)·(I+ml^{ 2 })-m^{ 2 }l^{ 2 }$$

$$\frac { θ(s) }{ U(s) } =\frac { \frac { ml }{ q } s }{ s^{ 3 }+\frac { b(I+ml^{ 2 }) }{ q } s^{ 2 }-\frac { (M+m)mgl }{ q } s-\frac { bmgl }{ q }  } $$

Por último si combinamos la ecuación (Ec.8) con la (Ec.5) obtendremos la entrada respecto del desplazamiento del carro:

$$\frac { X(s) }{ U(s) } =\frac { \frac { (I+ml^{ 2 })s^{ 2 }-mgl }{ q }  }{ s^{ 4 }+\frac { b(I+ml^{ 2 }) }{ q } s^{ 3 }-\frac { (M+m)mgl }{ q } s^{ 2 }-\frac { bmgl }{ q } s } $$

Si graficamos la función de transferencia del ángulo de inclinación del péndulo obtenemos lo siguiente:

Cómo se ve en la anterior con un simple impulso el ángulo de inclinación del péndulo crece desmesuradamente, de hecho, en 0.2s ya ha dado media vuelta.

La solución consiste en implementar un esquema de control en LC que permita alcanzar las especificaciones de diseño, de manera que se mantenga una referencia de ángulo cero respecto a la vertical para Φ a pesar de la fuerza de perturbación F. A este tipo de controladores que buscan mantener una referencia 0 se les llama a veces reguladores: 

Pero previamente, para tener una idea de la dinámica del sistema es aconsejable también ver la colocación de los polos y los ceros. Hallar dicha posición de polos y ceros y justificar las respuestas en LA obtenidas:

Si obtenemos los polos del sistema conseguimos el siguiente resultado:

$$s_{ 1 }=5.5650\quad s_{ 2 }=-5.6069\quad s_{ 3 }=-0.1428$$

Cómo podríamos haber deducido simplemente fijándonos en la respuesta del sistema, vemos cómo la función de transferencia tiene un polo en el semiplano derecho, es decir, un polo inestable, por lo tanto es lógico que el sistema no llegue a estabilizarse nunca. Además, en la función de transferencia podemos apreciar que el sistema tiene un cero en el origen, cosa que hará que nuestro lugar de las raíces se mueva hacia la derecha, quedándose muy cerca del eje imaginario por lo que comprometerá la estabilidad del sistema. 

Para ver si el sistema podría cumplir las especificaciones de control mediante un control proporcional sencillo se puede echar un vistazo al Lugar de las Raíces.

2. Hallar el lugar de las raíces

Considerando únicamente un control proporcional C(s)=K, ver si el sistema del péndulo P_pend puede cumplir las especificaciones:

Si utilizamos SISOTOOL y graficamos el lugar de las raíces poniéndole cómo condición que el sistema tenga un rebose del 5% vemos que es imposible conseguirlo, debido a que todas las constantes proporcionales que teóricamente van a cumplir el requisito de rebose, tienen un valor que nos desestabiliza el sistema.

3. Diseño de un controlador PID

Se implementa a continuación el esquema de control en LC visto anteriormente, donde se considera para C(s) un controlador PID.

Para hacerlo de una manera más gráfica, se utilizará Simulink, donde se construye el sistema de control en LC:

Hallar las constantes del controlador PID para que se cumplan las especificaciones de diseño requeridas, bien de manera heurística conociendo los efectos de cada una de las tres acciones del controlador, o utilizando el autosintonizado proporcionado por el bloque PID de Simulink.

Como siempre, recordar que a la hora de sintonizar PIDs se debe tener en cuenta que la sintonía PID está pensada para una entrada  y una salida, considerando otras posibles entradas o salidas como perturbaciones y considerándolas nulas a la hora de sintonizar, y suponiendo que el ajuste sea válido en presencia de dichas perturbaciones (haciendo un reajuste fino si es necesario). Esto es así tanto en procedimientos manuales (heurísticos, Z-N, etc.) como en el autosintonizado de Simulink.

Probar la bondad del ajuste realizado considerando como referencia de posición del péndulo Φ = 0 rad  -es decir, que se mantenga en posición vertical- ante una fuerza de perturbación aplicada al carro de F(t) = δ(t).

Sintonizando el controlador con el bloque de autosintonizado de Simulink conseguimos una respuesta con rebose 3.61% y un tiempo de establecimiento de 0.121s, siendo las constantes las que se muestran en la tabla.

Kp	282.976166021455
Ki	365.574479349202
Kd	48.6681477032831
N	25307.295950281

Si exponemos al sistema a una entrada impulso de 0.1s de duración y de magnitud  unitaria conseguimos la siguiente respuesta del péndulo. Cabe mencionar que este impulso será considerado como una perturbación ya que la referencia del sistema serán 0 rad debido a que nosotros no queremos que péndulo se incline.

Si analizamos con detenimiento la gráfica podemos ver cómo cuando el péndulo es sometido a un impulso este se inclina un poco y enseguida vuelve a la verticalidad. Si cuantificamos lo simulado vemos cómo el péndulo sólo se inclina radianes ante una entrada de   y recupera la verticalidad en menos de 2s. Por lo tanto, las especificaciones establecidas al principio se cumplen.

Ver qué pasa con la posición del carro en base al controlador diseñado y razonar por qué.

Para acabar gráficaremos el desplazamiento del carro para el caso anterior, cabe mencionar que a el carro no se le aplicará ningún controlador.

Si analizamos el desplazamiento del carro vemos como cunado recibe el impulso se desplaza hacia delante, pero como rápidamente se desplaza hacia atrás para que el péndulo recupere su verticalidad. A continuación, el carro se vuelve a mover hacia delante para compensar el sobre impulso del péndulo hasta que el péndulo consigue la verticalidad y el carro se para.