Modelado de un quadcoptero

En esta entrada vamos a abordar el modelado y el diseño de un sistema de control basado en PIDs para un quadcóptero comercial. Para ello nos basáremos en estos dos artículos:

• Hossein Bolandi, Mohammad Rezaei, Reza Mohsenipour, Hossein Nemati, Seed Majid Smailzadeh. (2013). Attitude Control of a Quadrotor with Optimized PIDController, In Script Intelligent Control and Automation, August 2013, 4, 335-342.
• Voos, H. (2009). Nonlinear Control of a Quadrotor Micro-UAVusing Feedback Linearization, Proceedings of the 2009 IEEE International Conference on Mechatronics, Malaga, Spain, April 14-17, ISBN 978-1-4244-4195-2

En los últimos años el uso de los Vehículos Aéreos No Tripulados (Unmanned Aerial Unmanned ó UAV) ha aumentado considerablemente en diversos ámbitos, desde la cartografía, hasta aplicaciones de agricultura rescate o seguridad. Dentro de los UAV se puede encontrar un sinfín de modelos, de los cuales uno de los más comunes es el quadcóptero.

El quadcóptero es un helicóptero multirotor propulsado por cuatro hélices.  Cabe destacar, que el quadcóptero se clasifica como helicóptero debido a que su sustentación en el aire se consigue mediante cuatro hélices orientadas verticalmente y no por dos alas fijas, como en las aeronaves tradicionales.

Antes de comenzar a modelar el quadcóptero, es necesario estudiar la configuración de sus hélices. En este tipo de UAV se utilizan dos conjuntos de hélices, las dos del primer conjunto (M1 y M3) giran en la dirección de las agujas del reloj, mientras que las otras dos (M2 y M4) giran en el sentido contrario.

Otro aspecto a tener en cuenta son los movimientos que puede realizar el quadcóptero, que en este caso serán, los desplazamientos en x, y, z y la rotación alrededor de estos ejes. La notación que se utilizará para la rotación en los ejes será la de Roll (𝜙)   para la rotación alrededor del eje x, Pitch (𝜃)  para la rotación alrededor del eje y, y Yaw (Ψ) para la rotación alrededor del eje z.

1. Modelado

Para modelar el sistema se tiene que tener en cuenta que el quadcóptero será considerado como un sólido rígido simétrico, cuyo centro de gravedad estará situado en el centro del mismo.

Antes de aplicar el sumatorio de fuerzas en los ejes, cabe destacar que la orientación del dron viene dada por sus ejes 𝜙, 𝜃 y Ψ y que además se puede calcular  con la matriz de rotación R , que se obtiene de realizar un giro de un ángulo 𝜙 alrededor del eje x, otro giro  𝜃 alrededor del eje y, y por último un giro Y sobre z, todo esto respecto del sistema de referencia fijo.

$$R=\left( \begin{matrix} S_{ ϕ }S_{ θ }S_{ ψ }+C_{ θ }C_{ ψ } & \quad S_{ ϕ }S_{ θ }S_{ ψ }-S_{ θ }S_{ ψ } & \quad S_{ θ }C_{ ϕ } \\ C_{ θ }S_{ ψ } & C_{ ϕ }C_{ ψ } & -S_{ ϕ } \\ S_{ ϕ }C_{ θ }S_{ ψ }-S_{ θ }S_{ ψ } & \quad S_{ ϕ }S_{ θ }S_{ ψ }+S_{ θ }S_{ ψ } & \quad C_{ θ }C_{ ψ } \end{matrix} \right) $$ (Ec.1)

Aplicando el sumatorio de fuerzas en los tres ejes se obtiene la ecuación siguiente:

$$\left( \begin{matrix} mẍ \\ mÿ \\ mz̈ \end{matrix} \right) =\left( \begin{matrix} 0 \\ 0 \\ mg \end{matrix} \right) -R\sum { { F }_{ i }\left( \begin{matrix} 0 \\ 0 \\ 1 \end{matrix} \right)  } $$ (Ec.2)

De esta expresión se obtienen las aceleraciones en los ejes x, y, z. Teniendo en cuenta que,  F_i=bw_i², siendo b el coeficiente de empuje de las hélices (Ns/rad), es decir, la fuerza que producen las hélices en función de la velocidad angular del motor, y u₁ la señal de control :

$$\sum { { F }_{ i } } ={ u }_{ 1 }=b(w_{ 1 }^{ 2 }+w_{ 2 }^{ 2 }+w_{ 3 }^{ 2 }+w_{ 4 }^{ 2 })$$(Ec.3)

$$ẍ=-\frac { { u }_{ 1 } }{ m } S_{ θ }C_{ ϕ }$$(Ec.4)

$$\ddot { y } =-\frac { { u }_{ 1 } }{ m } S_{ ϕ }$$(Ec.5)

$$\ddot { z } =-\frac { { u }_{ 1 } }{ m } C_{ ϕ }C_{ θ }+g$$(Ec.6)

Para calcular el valor de las aceleraciones angulares se utilizará la ecuación de momentos en un sólido rígido teniendo en cuenta que d  hace referencia al coeficiente de arrastre (Nms/rad), es decir, el par que crean las hélices en función de la velocidad angular de los motores, e I es la matriz de inercias, que en este caso, al tratarse de un objeto simétrico será una matriz diagonal. Por último L hace referencia a la distancia entre el centro de masa de quadcóptero a los motores:

$$u_{ 2 }=b(ω_{ 1 }^{ 2 }-ω_{ 2 }^{ 2 }-ω_{ 3 }^{ 2 }+ω_{ 4 }^{ 2 })$$(Ec.7)

$$u_3=b(ω_1^2+ω_2^2-ω_3^2-ω_4^2)$$(Ec.8)

$$u_4=d(-ω_1^2+ω_2^2-ω_3^2+ω_4^2)$$(Ec.9)

$$ τ=I\dot { ω } +ω×Iω→\dot { ω } =(τ-ω×Iω){ I }^{ -1 }$$(Ec.10)

$$\left( \begin{matrix} \ddot { ϕ }  \\ \ddot { θ }  \\ \ddot { ψ }  \end{matrix} \right) =\left( \begin{matrix} \frac { I_{ y }-I_{ z } }{ { I }_{ x } } \dot { θ } \dot { ψ } +\frac { L }{ { I }_{ x } }  \\ \frac { I_{ z }-I_{ x } }{ { I }_{ y } } \dot { ϕ } \dot { ψ } +\frac { L }{ { I }_{ y } }  \\ \frac { I_{ x }-I_{ y } }{ { I }_{ z } } \dot { ϕ } \dot { θ } +\frac { 1 }{ { I }_{ z } }  \end{matrix} \right) $$(Ec.11)

De esta forma, las ecuaciones 4, 5, 6 y 11 sintetizan el modelo del quadcóptero. Por un lado, las ecuaciones 4, 5 y 6 corresponden a las magnitudes lineales, y relacionan las velocidades angulares de los motores con las aceleraciones lineales y la señal de control u1. En lo que respecta a las aceleraciones angulares, la ec. 11 relaciona las velocidades angulares del sistema de referencia OXYZ con las aceleraciones y las señales de control u₂, u₃ y u₄, además de con los parámetros geométricos (L) y los momentos de inercia. Por tanto, las ecuaciones anteriormente mencionadas (ecs. 4-11) forman el modelo no lineal del quadcóptero.   Si se introducen las ecuaciones en un modelo de Simulink se consigue lo siguiente: 

Antes de comenzar con el diseño de controladores tipo PID, es necesario obtener un modelo lineal y para ello se linealizan las ecuaciones en un punto de operación. Del trabajo de Bolandi, se selecciona el mismo punto de operación para poder comparar los resultados posteriormente, siendo los valores:

$$ϕ_{ 0 }=θ_{ 0 }=ψ_{ 0 }=0\quad \quad u_{ 20 }=u_{ 30 }=u_{ 40 }=0$$

Una vez definidas dichas variables en el punto de operación, se calcula u₁₀=0, utilizando para ello la ecuación 6 en dicho punto de operación: 

$$u_{ 10 }=\frac { (-\ddot { { z }_{ 0 } } +g)m }{ { C }_{ { θ }_{ 0 } }{ C }_{ { ϕ }_{ 0 } } } $$

Por tanto, partiendo del modelo no lineal y siguiendo el método de aproximación de series de Taylor, se pueden obtener las funciones de transferencia lineales:

$$\frac { \dot { X } (s) }{ θ(s) } =\frac { -g }{ s } $$(Ec.12)

$$\frac { \dot { Y } (s) }{ θ(s) } =\frac { g }{ s } $$(Ec.13)

$$\frac { \dot { Z } (s) }{ { U }_{ 1 }(s) } =\frac { -1/m }{ s } $$(Ec.14)

$$\frac { ϕ(s) }{ { U }_{ 2 }(s) } =\frac { L/I_{ x } }{ { s }^{ 2 } } $$(Ec.15)

$$\frac { θ(s) }{ { U }_{ 3 }(s) } =\frac { L/I_{ y } }{ { s }^{ 2 } } $$(Ec.16)

$$\frac { \dot { ψ } (s) }{ { U }_{ 4 }(s) } =\frac { 1/I_{ z } }{ { s } } $$(Ec.17)

Combinando las ecuaciones lineales 12 y 16, se conseguirá la velocidad en el eje x en función de la entrada u₃ por un lado; mediante la combinación de las ecuaciones  13 y 15 se conseguirá la velocidad en el eje y en función de la señal de control u_2:

$$\frac { \dot { X } (s) }{ { U }_{ 3 }(s) } =\frac { -g }{ s } \frac { L/I_{ y} }{ { { s }^{ 2 } } }$$(Ec.18)

$$\frac { \dot { Y } (s) }{ { U }_{ 2 }(s) } =\frac { g }{ s } \frac { L/I_{ x } }{ { { s }^{ 2 } } }$$ (Ec.19)

De esta forma se consiguen 4 funciones de transferencia, una por cada señal de control (expresiones 14, 17, 18 y 19). De esta forma se puede construir un modelo lineal en Simulink. 

2. Diseño de los controladores

Una vez obtenido el modelo lineal del quadcóptero (ecs 14,17, 18 y 19) en este apartado se va a llevar a  cabo el diseño de varios controladores con el objetivo de asegurar el seguimiento de las consignas de velocidad en los ejes x, y, z así como en la velocidad angular Yaw. Cabe destacar que en el diseño de los controladores se tendrán en cuenta dos casos diferentes. Para el primer caso, se diseñarán dos bucles de control simple, uno para el control de la velocidad en el eje z y otro para la velocidad de giro del ángulo Yaw. En el segundo caso, se diseñaran  dos  controladores  en  cascada,  uno  para  la velocidad en el eje x y la posición del ángulo pitch, y el otro para la velocidad en el eje y y la posición del ángulo roll.

 2.1 Diseño de controladores para el seguimiento de velocidad en el eje z y para la velocidad del ángulo yaw

Para diseñar los bucles de control para la velocidad en el eje z y la velocidad del ángulo yaw se diseñarán dos bucles de control simple. Para ello, en primer lugar, se calculará la ecuación característica de la función de transferencia, teniendo en cuenta el controlador. Para este primer caso se ha decidido diseñar un PI debido a la necesidad de conseguir un error nulo en estado estacionario. Una vez calculada la ecuación característica del sistema se obtendrá la ecuación (Ec.20) teniendo en cuenta que A será el numerador de la función de transferencia la cual se va controlar, en este caso las ecuaciones 14 y 17:

$$0=1+G(s)H(s)=s^2+AK_p+AK_p/T_i$$(Ec.20)

Si se iguala la ec. 20 a la ecuación característica de  segundo grado, es posible  calcular  las constantes del controlador siempre y cuando se fije un valor para  y  ya que los otros parámetros dependen del quadcóptero y son conocidos.

$$s^{ 2 }+AK_{ p }+\frac { A{ K }_{ p } }{ T_{ i } } =s^{ 2 }+2δw_{ n }s+wn^{ 2 }$$

$$K_{ p }=\frac { 2δwn }{ A } \quad \quad y\quad \quad { T }_{ i }=\frac { AK_{ p } }{ wn^{ 2 } } $$

2.2 Diseño de controladores para el seguimiento de velocidad en los ejes x e y 

Para el control de las velocidades en los ejes x e y se diseñarán dos bucles de control en cascada, donde el bucle secundario se encargará del control de la posición de los ángulos pitch y roll y el bucle primario del control de las velocidades lineales en los ejes x e y. Para ello, se debe tener en cuenta que las funciones de transferencia de los bucles secundarios tienen dos  polos en el origen (ecuaciones 15 y 16), por lo que se utilizarán controladores tipo PD en lugar de PI. Por lo tanto, utilizando la misma metodología utilizada para el diseño de los bucles de control simples del apartado anterior, es decir, igualando la ecuación característica de la función de transferencia a la ecuación genérica de segundo grado, el cálculo de las constantes de los PD quedará de la siguiente manera, asumiendo que B serán los numeradores de las funciones de trasferencia 15 y 16. 

$$K_{ p }=\frac { { w }_{ n }^{ 2 } }{ B } \quad \quad y\quad \quad { T }_{ d }=\frac { 2wnδ }{ K_{ p }B } $$

3. Resultados 

Para validar los modelos lineales y no lineales (ecuaciones 14, 17, 18 y 10) se han utilizado los valores del dron comercial que nos proporciona Voos en su artículo. 

Introduciendo como señales de entrada escalones unitarios en las velocidades de los x, y y z así como en la velocidad del ángulo yaw,  se obtiene los siguientes resultados.

Los resultados muestran salidas muy parecidas para ambos modelos, excepto en la velocidad del eje z. En este caso, la respuesta del modelo no lineal presenta un sobre impulso muy acusado debido a un cero dominante. 

Tras este ajuste se puede afirmar que la linealización del modelo lineal es una buena aproximación, ya que habiendo hecho el diseño de los controladores teniendo únicamente en cuenta el modelo lineal, excepto en la velocidad del eje z,  se consiguen buenos resultados en el modelo no lineal.

Una vez ajustados los controladores, se pasará a validar el modelo tomando como referencia los resultados del trabajo de Bolandi, Para ello, se verá  ver la evolución de la posición de  los ángulos roll, pitch y yaw ante una perturbación escalón unitaria.

Para la segunda prueba, se utilizarán referencias sinusoidales de amplitud unitaria y frecuencia 1rad/s,  para las velocidades en x e y, además de una referencia escalón unitario en la velocidad en el eje z, para simular así un movimiento ascendente en espiral. De esta forma se analizará el comportamiento de las velocidades en x y en y a la hora de seguir una trayectoria.

Como se puede apreciar, los bucles de control de las velocidades en x e y responden de una manera muy satisfactoria al seguimiento de las referencias, siendo el valor máximo del error de 0.1, es decir un 10%,  sin tener en cuenta el estado transitorio de la velocidad en y, que cómo se puede ver en la Fig.8 la consigna comienza en 1 y no en 0 como en la velocidad en el eje x.