Modelado y control del sistema de suspensión de un autobús

En este Post explicaremos cómo modelar el sistema de suspensión de un autobús, y como aplicar un controlador en lazo cerrado para conseguir una amortiguación lo más suave posible. Para ello, primero abordaremos el modelado del sistema de suspensión y después el diseño del sistema de control. 

Consideraremos el sistema de suspensión automático de la figura:

Cabe destacar que se ha esquematizado en el diagrama de cuerpo libre la suspensión correspondiente a una de las 4 ruedas.

El sistema de suspensión automático debe ser tal que, ante perturbaciones en la carretera, el autobús no presente grades oscilaciones y éstas se disipen rápidamente. Teniendo en cuenta que la deformación del neumático es en general pequeña se puede considerar que la variación de la altura del autobús al suelo es la misma que la de X₁-X_2, de forma que se tomará esta distancia como salida.

Se desea controlar la fuerza que ha de ejercer la suspensión activa (u) para que, ante un cambio brusco en el firme (W) representado por un escalón unitario, la respuesta (X₁-X₂) presente las siguientes características:

   Tiempo de asentamiento < 5s.

          Rebose < 5%

Esto es, cuando el autobús salga de un socavón de 10cm de altura, el cuerpo del autobús oscilará un máximo de 5mm y dejará de oscilar en menos de 5seg.

Suponer para simplificar que la fuerza de la gravedad es despreciable respecto las fuerzas correspondientes a las constantes elásticas de la suspensión.  

Los parámetros del sistema son los siguientes:

M1 (1/4 masa del autobus)= 2500 kg

M2 (masa de la propia suspensión)= 320 kg

K1 (cte elástica de la suspension)= 80000 N/m

b1 (coef de amortiguamiento de la suspensión)= 350 N.s/m

K2(cte elástica del neumático)= 500000 N/m

b2 (coef de amortiguamiento de la rueda)= 15020 N.s/m

U (fuerza aplicada por el control)

Primeramente siempre se ha de observar la respuesta en LA, para ver si el sistema es estable o no,  y si lo es cómo de alejada está su dinámica respecto de las especificaciones deseadas.

1. Hallar la salida del sistema en LA a un escalón unitario

Lo primero que hay que hacer es modelar el sistema, obteniendo la expresión de la salida X₁-X₂ ante las dos entradas U y W.

       Masa1

$$M_{ 1 }(\ddot { { x }_{ 1 } } )+b_{ 1 }((\dot { { x }_{ 1 } } )-(\dot { { x }_{ 2 } } ))+k_{ 1 }(x_{ 1 }-x_{ 2 })=u$$          (Ec.1)

       Masa2

$$M_{ 2 }(\ddot { { x }_{ 2 } } )+b_{ 2 }(\dot { { x }_{ 2 } } -\dot { w } )+k_{ 2 }(x_{ 2 }-w)=b_{ 1 }(\dot { { x }_{ 1 } } -\dot { { x }_{ 2 } } )+k_{ 1 }(x_{ 1 }-x_{ 2 })-u$$       (Ec.2)

Si aplicamos Laplace:

      Masa1

$$M_1 s^2 X_1 (s)+(b_1 s+k_1 )(X_1 (s)-X_2 (s))=U(s) $$                       (Ec.3)

     Masa2

$$M_2 s^2 X_2 (s)+(b_2 s+k_2 )(X_2 (s)-W(s))=(b_1 s+k_1 )(X_1 (s)-X_2 (s))-U(s)$$      (Ec.4)

Para sacar X₁-X₂ en función de U(s) y W(s) primero despejaremos X₁(s)  de la ecuación (Ec.2) dejándolo en función de U(s) y X₁(s)-X₂(s).

$$X_{ 1 }(s)=\frac { U(s)-(b_{ 1 }s+k_{ 1 })(X_{ 1 }(s)-X_{ 2 }(s)) }{ M_{ 1 }s^{ 2 } } $$  (Ec.5)

Para continuar despejaremos X₂(s) de la ecuación (Ec.4) dejándolo en función de W(s), U(s) y X₁(s)-X₂(s).

$${ X }_{ 2 }(s)=\frac { (b_{ 1 }s+k_{ 1 })(X_{ 1 }(s)-X_{ 2 }(s))-U(s)+(b_{ 2 }s+k_{ 2 })W(s) }{ M_{ 2 }s^{ 2 }+b_{ 2 }s+k_{ 2 } } $$(Ec.6)

Por último para conseguir X₁(s)-X₂(s) en función de U(s) y W(s)  restaremos las ecuaciones (Ec.5) y (Ec.6):

$${ { X }_{ 1 }(s)-X }_{ 2 }(s)=\frac { U(s)((M_{ 1 }+M_{ 2 })s^{ 2 }+b_{ 2 }s+k_{ 2 })-W(s)(M_{ 1 }b_{ 2 }s^{ 3 }+M_{ 1 }k_{ 2 }s^{ 2 }) }{ M_{ 1 }M_{ 2 }s^{ 4 }+(M_{ 1 }b_{ 1 }+M_{ 1 }b_{ 2 }+M_{ 2 }b_{ 1 })s^{ 3 }+(M_{ 1 }k_{ 1 }+M_{ 1 }k_{ 2 }+M_{ 2 }k_{ 1 })s^{ 2 }+(b_{ 1 }k_{ 2 }+k_{ 1 }b_{ 2 })s+k_{ 1 }k_{ 2 } } $$ (Ec.7)

A continuación, hallamos la salida del sistema ante entradas escalón unitario. Aplicando el principio de superposición puesto que la T.L. es lineal, la salida de un sistema a dos entradas es igual a la suma de las salidas a cada entrada por separado.

De manera para considerar la entrada de control U(s) únicamente se pone W(s)=0, obteniendo la F.T. correspondiente:

$${ G }_{ 1 }(s)=\frac { { X }_{ 1 }(s)-{ X }_{ 2 }(s) }{ U(s) }$$ 

 $${ G }_{ 1 }(s)=\frac { (M_{ 1 }+M_{ 2 })s^{ 2 }+b_{ 2 }s+k_{ 2 } }{ M_1 M_2 s^4+(M_1 b_1+M_1 b_2+M_2 b_1 ) s^3+(M_1 k_1+M_1 k_2+M_2 k_1 ) s^2+(b_1 k_2+k_1 b_2 )s+k_1 k_2 }$$(Ec.8)

$${ G }_{ 2}(s)=\frac { { X }_{ 1 }(s)-{ X }_{ 2 }(s) }{ W(s) }$$ 

$${ { G }_{ 2 } }(s)=\frac { -(M_{ 1 }b_{ 2 }s^{ 3 }+M_{ 1 }k_{ 2 }s^{ 2 }) }{ M_1 M_2 s^4+(M_1 b_1+M_1 b_2+M_2 b_1 ) s^3+(M_1 k_1+M_1 k_2+M_2 k_1 ) s^2+(b_1 k_2+k_1 b_2 )s+k_1 k_2 }$$ (Ec.9)

Si graficamos las salidas aplicando un entrada escalón unitaria se consiguen los siguientes resultados. 

G1

G2

Cómo se puede apreciar en las figuras anteriores ambas salidas tienen mucha oscilación y error en régimen permanente. Además el sobre impulso en ambas es mayor que 5% y el tiempo de establecimiento más del doble de lo especificado.

La solución por tanto consiste en implementar un esquema de control en LC que permita alcanzar las especificaciones de diseño:

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.

Los polos de los sistemas G1 y G2 sistema pueden hallarse directamente mediante el comando roots o el comando pole siempre que utilicemos Matlab.

Que da el siguiente resultado.

 -23.9758 +35.1869i

 -23.9758 -35.1869i

  -0.1098 + 5.2504i

  -0.1098 - 5.2504i

Como se puede ver hay dos polos muy cerca del eje imaginario, por tanto son sistemas que pueden tener la estabilidad comprometida.

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, valorar si el sistema G1 puede cumplir las especificaciones.

Si utilizamos la expresión de cálculo del coeficiente de amortiguamiento para un sistema de orden 4 conseguimos el siguiente resultado:

$$ς=\frac { ln⁡\left( \frac { R }{ 100 }  \right)  }{ \sqrt { π^{ 2 }+\left[ ln⁡\left( \frac { R }{ 100 }  \right)  \right] ^{ 2 } }  } =0.69$$

Por lo tanto si utilizamos Sisotool y ponemos cómo requisito que el coeficiente de amortiguamiento es 0.69 conseguimos el siguiente gráfico.

Cómo se puede ver en la figura  no podemos conseguir la especificaciones requeridas con un simple control proporcional, que no podemos hacer que el lugar de las raíces queda completamente fuera de la parte blanca del dibujo.

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.

Probar la bondad del ajuste realizado considerando como referencia de altura del autobus r(t) un escalón unitario, de manera que una vez llegado al estado estacionario -esto es, una vez el autobús está en marcha normal a una altura constante- por ejemplo en t=5 pasa un bache de magnitud 0.1 y de 1s de duración.

Sintonizando el controlador con el bloque de autosintonizado de Simulink conseguimos una respuesta con rebose 0% y un tiempo de establecimiento de 0.0025s, siendo las constantes las que s emuestran en la tabla1.

Kp = 18742064.957701

Ki = 173490373.039331

Kd = 449863.199782723

N = 181333.345428777

Si exponemos al sistema a una entrada escalón y un bache de 0.1 de un segundo de duración en el segundo 5 esta es la salida que obtenemos.

Si analizamos con detenimiento la gráfica podemos ver como el sistema apenas tiene rebose, en el peor de los casos, es decir en el momento del bache, presenta un rebose del 0.3% y además se estabiliza en  1s con precisión del 0.05%.