Modelado, simulación y control de un sistema de seguimiento de trayectorias para UAVs

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Proyecto Fin de Máster Ingeniería Aeronáutica Modelado, simulación y control de un sistema de seguimiento de trayectorias para UAVs Autor: Antonio José Aguilera Albendín Tutor: Eduardo Fernández Camacho Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Proyecto Fin de Máster Ingeniería Aeronáutica Modelado, simulación y control de un sistema de seguimiento de trayectorias para UAVs Autor: Antonio José Aguilera Albendín Tutor: Eduardo Fernández Camacho Catedrático Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Proyecto Fin de Máster: Modelado, simulación y control de un sistema de seguimiento de trayectorias para UAVs Autor: Antonio José Aguilera Albendín Tutor: Eduardo Fernández Camacho El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos La realización de este Trabajo Fin de Máster supone la finalización de mi etapa universitara y es por eso que quisiera dar las gracias a todos aquellos que han hecho que sea posible. En primer lugar está claro a quiénes debo agradecer todo su esfuerzo y sacrificio: a mis padres, Antonio y M ª Carmen, ya que sin ellos no habría sido posible que hoy tuviera la formación que tengo y es algo que nunca podré terminar de agradecérselo. Tampoco me puedo olvidar de mi familia, en especial de mi abuela Marina y de mi hermano, Juan Carlos. En segundo lugar quiero agradecer a todos aquellos que me han ayudado cuando lo he necesitado: amigos, compañeros de clase, compañeros de piso, profesores. Ellos también han contribuido a que llegue este momento. Finalmente quisiera darle las gracias a mi Tutor, Eduardo, sin el cual la realización de este Proyecto no habría sido factible. Antonio José Aguilera Albendín Sevilla, 2019 I Resumen Los drones o UAVs tienen un gran potencial en áreas muy diversas, ya que pueden desplazarse rápidamente sobre un terreno irregular o accidentado y superar cualquier tipo de obstáculo ofreciendo imágenes a vista de pájaro y otro tipo de información recogida por diferentes sensores. Su utilización ha crecido exponencialemnte durante los último años debido a la gran cantidad de aplicaciones donde pueden ser muy útiles: búsqueda de personas desaparecidas, fotografía, video y cartografía aérea, prevención y control de incendios, agricultura, medio ambiente, construcción e inspecciones, control y análisis de multitudes, exploración de lugares de díficil acceso, etc. Al tratarse de aeronaves autónomas es necesario dotarlas de un sistema de misión que se encargue de implementar la lógica necesaria para poder realizar la misión deseada. En la literatura podemos encontrar numerosas técnicas de control y guiado para un UAV: control mediante bucle cerrado (PID), control predictivo (MPC), neural network guidance, métodos de control y guiado directos, etc, teniendo cada uno sus ventajas e inconvenientes. En este proyecto se ha hecho uso del control PID debido a su buen funcionamiento y amplio uso en la industria. Además, debido a que el objetivo del proyecto no reside en la optimización de trayectorias, el control PID permite el guiado del UAV sin un alto gasto computacional. En el proyecto desarrollado se modela el UAV X8 Skywalker, con configuración de ala volante. El modelo matemático utilizado para realizar la simulación consiste en un método dinámico no lineal que describe el comportamiento del UAV de la manera más realista posible. Se le dotará de un sistema de misión, basado en control mediante bucle cerrado, y diversos sensores útiles para diferentes tipos de aplicaciones. Una vez implementado el sistema que modela el UAV X8 Skywalker se realizarán diversas simulaciones para estudiar el funcionamiento del sistema de misión implementado y su respuesta frente a perturbaciones atmosféricas. Además se mostrará cómo realizar la implementación con el simulador de vuelo FlightGear, software de código libre que permite ver en tiempo real la actuación de la aeronave durante la simulación. III 1 Introducción 1 Introducción 1.1 Alcance del Proyecto El objetivo de este proyecto consiste en desarrollar un modelo matemático que describa el comportamiento dinámico del UAV Skywalker X8 y posteriormente estudiar diferentes misiones que pueda llevar a cabo. En el Capítulo 2 se desarrollan las ecuaciones que modelan el comportamiento de la aeronave, que corresponde a un sistema dinámico no lineal. Para ello haremos uso de los parámetros que modelan la aeronave presentes en [1]. Una vez que tenemos las ecuaciones que modelan la aeronave, se va a proveer al UAV de los sistemas necesarios para realizar la misión que nos interesa: un vuelo horizontal en el que de forma autónoma sea capaz de pasar por una serie de waypoints definidos. Para ello se llevará a cabo la implementación de un sistema autopiloto y de misión en los Capítulos 3 y 4. Además, en el Capítulo 5 se van a modelar diversos sensores que proporcionan datos de vuelo interesantes para poder implementar un Observador basado en el Filtro de Kalman. En los Capítulos 6 y 7 se muestra cómo se ha implementado dicho sistema en el software Matlab-Simulink y los resultados obtenidos durante la simulación, respectivamente. En el Capítulo 8 se va a estudiar cómo afectan las perturbaciones atmosféricas debidas a viento en el comportamiento del UAV. Finalmente, se realizará una implementación en el simulador de vuelo FlightGear y se abordarán las conclusiones que se obtienen a la finalización de la elaboración del proyecto, así como posibles vías de desarrollo futuro. 1.2 UAV de partida La aeronave que se pretende simular en este proyecto se trata del X8 Skywalker. Corresponde a un UAV de tipo ala volante con las siguientes características físicas: •Masa: m= 3.797 kg •Envergadura: b= 2.10 m •Cuerda: c= 0.3571 m •Longitud: L= 0.78 m •Superficie alar: S= 0.75 m2 las cuales han sido obtenidas de [ 1 ]. Al tratarse de un ala volante no dispone de timón de profundidad (rudder), por lo que el control de guiñada y balance se encuentra acoplado, tal y como se desarrollará en el Capítulo 3. Además al tener una masa tan pequeña, su comportamiento se verá condicionado por la intensidad y dirección del viento incidente, desarrollado en el Capítulo 8. 1 Introducción 1 2Capítulo 1. Figura 1.1 X8 Skywalker.. 1.3 Sistemas de referencia Para el desarrollo de las ecuaciones que definen la cinemática y dinámica de la aeronave es necesario hacer uso de diferentes sistemas de referencia debido a la gran cantidad de acciones a las que se encuentra sometida. Los sistemas empleados son: •Sistema de Ejes Tierra (ECEF, Earth Centered, Earth Fixed) : empleado cuando se quiere conocer la posición respecto de la Tierra, ya que es solidario con el movimiento de ella. Los parámetros empleados para conocer la posición en este sistema de ejes son la altitud (h), longitud ( λ ) y latitud ( ϕ ). Figura 1.2 Sistema de coordenadas de Ejes Tierra.. •Sistema de Ejes de Navegación (NED, North, East, Down) : sistema de horizonte local. Se encuentra referenciado respecto a un punto que se puede encontrar en la superficie terrestre o no (aeronave). Está definido por las direcciones Norte (N), Este (E) y vetical hacia abajo (D). •Sistema de Ejes Cuerpo (BFS, Body Fixed System) : sistema de referencia empleado para definir la actitud de la aeronave respecto al sistema de Ejes de Navegación. Caracterizado por los ángulos de Euler de alabeo (φ), cabeceo (θ) y guiñada (ψ). •Sistema de Ejes Viento (WF, Wind Frame) : empleado para orientar el viento respecto de la aeronave. Definido por el ángulo de ataque (α) y ángulo de resbalamiento (β). A partir de las matrices de cosenos directores (DCM) es posible pasar de un sistema de referencia a otro. Multiplicando un vector expresado en un determinado sistema de referencia por la correspondiente matriz de cosenos directores queda expresado en el sistema de referencia que sea de interés. Vb=Rb aVa 1.3 Sistemas de referencia 3 Figura 1.3 Sistema de Ejes de Navegación.. Figura 1.4 Sistema de Ejes Cuerpo.. Figura 1.5 Sistema de Ejes Viento.. Una particularidad de estas matrices es que son ortogonales ( Ra b = Rb aT ), lo cual permite transformar en sentido contrario multiplicando solamente por la matriz traspuesta Va=Ra bVb= (Rb a)TVb(1.1) Es posible pasar de un sistema de referencia a otro a partir de sistemas de referencia intermedios, es decir, la matriz de cosenos directores que permite transformar un vector expresado en un determinado sistema de referencia en otro se puede expresar mediante el producto de matrices de cosenos directores intermedias A→B→C→D 4Capítulo 1. Vd=Rd aVa=Rd cRc bRb aVa(1.2) Las matrices de cosenos directores que interesan para la consecución de este proyecto son aquellas que permiten pasar de un sistema de referencia a otro de los sistemas desarrollados anteriormente •Matriz de rotación de ejes Tierra a ejes de navegación: Rn g=  −sin(ϕ)cos(λ)−sin(ϕ)sin(λ) cos(ϕ) −sin(λ) cos(λ) 0 −cos(ϕ)cos(λ)−cos(ϕ)sin(λ)−sin(ϕ) (1.3) •Matriz de rotación de ejes de navegación a ejes cuerpo: Rb n=  cθcψcθsψ−sθ −cφsψ+sφsθcψcφcψ+sφsθsψsφcθ sφsψ+cφsθcψ−sφcψ+cφsθsψcφcθ (1.4) •Matriz de rotación de ejes viento a ejes cuerpo: Rb w=  cos(β)cos(α)−sin(β)cos(α)−sin(α) sin(β) cos(β) 0 cos(β)sin(α)−sin(β)sin(α) cos(α) (1.5) donde se ha empleado la notación cx≡cos(x),sx≡sin(x). 2 Modelado 2.1 Modelo cinemático y dinámico En este apartado se desarrollarán las ecuaciones que modelan el comportamiento del X8. Primero se realizará un análisis de la cinemática que describe el movimiento de la aeronave y posteriormente se desarrollará la dinámica que produce dicho movimiento. Como se verá a continuación, se trata de un sistema dinámico no lineal, por lo que las variables de estado longitudinal y lateral de la aeronave se encontrarán acopladas, siendo necesario por tanto estudiar de forma conjunta la dinámica longitudinal y lateral. 2.1.1 Variables de estado Las ecuaciones de movimiento para un UAV se describen mediante doce variables de estado. Hay tres estados de posición y tres estados de velocidad asociados con el movimiento de traslación del UAV. De manera similar, hay tres posiciones angulares y tres estados de velocidad angular asociados con el movimiento de rotación. Dichas variables se encuentran recogidas en la tabla 2.1. Tabla 2.1 Variables de estado.. Símbolo Descripción PnPosición inercial del UAV respecto al norte en ejes NED. PePosición inercial del UAV respecto al este en ejes NED. PdPosición inercial del UAV respecto al centro de la tierra en ejes NED. uVelocidad del UAV en dirección longitudinal en ejes Cuerpo. vVelocidad del UAV en dirección lateral en ejes Cuerpo. wVelocidad del UAV en dirección vertical en ejes Cuerpo. φÁngulo de balance (roll). θÁngulo de cabeceo (pitch). ψÁngulo de alabeo (yaw). PVelocidad angular de balance. QVelocidad angular de cabeceo. RVelocidad angular de guiñada. 2.1.2 Cinemática El modelo cinemático desarrollado a continuación corresponde a las ecuaciones que permiten calcular la derivada del vector posición de la aeronave ( ˙pn,˙pe,˙pd) , expresado en ejes de navegación. Dichas derivadas se calculan a partir de las componentes de la velocidad lineal (u,v,w) y su orientación (φ,θ,ψ) , expresada mediante ángulos de Euler.   ˙pn ˙pe ˙pd =  cθcψsφsθcψ−cφsψcφsθcψ+sφsψ cθsψsφsθsψ+cφcψcφsθsψ−sφcψ −sθsφcθcφcθ   u v w (2.1) 5 6Capítulo 2. Modelado Figura 2.1 Variables de estado.. donde se ha empleado la notación cx≡cos(x),sx≡sin(x). Respecto a los ángulos de Euler empleados en (2.1) para expresar los ángulos de alabeo (φ) , cabeceo (θ) y guiñada (ψ) , se pueden obtener sus derivadas (˙ φ, ˙ θ, ˙ ψ) a partir de las componentes de la velocidad angular (p,q,r).   ˙ φ ˙ θ ˙ ψ =  1 sin(φ)tan(θ) cos(φ)tan(θ) 0 cos(φ)−sin(φ) 0 sin(φ)sec(θ) cos(φ)sec(θ)   p q r (2.2) Integrando los resultados obtenidos en (2.1) y (2.2) se obtienen la posición y actitud de la aeronave, respectivamente. 2.1.3 Dinámica Para la obtención del modelo que representa la dinámica de la aeronave se va a considerar un modelo de Tierra plana, el cual es apropiado para vehículos pequeños, y a emplear la segunda Ley de Newton, para la cual es necesario tener un sistema de referencia inercial al que se encuentre asociado el movimiento de la aeronave. Así, se puede obtener la aceleración lineal ( ˙u,˙v, ˙w) y angular ( ˙p, ˙q, ˙r) de la aeronave, expresadas en ejes de navegación.   ˙u ˙v ˙w =  rv −qw pw −ru qu −pv  +1 m  fx fy fz (2.3)   ˙p ˙q ˙r =  Γ1pq −Γ2qr Γ5pr −Γ6(p2−r2) Γ7pq −Γ1qr  +  Γ3l+Γ4n m Iy Γ4l+Γ8n (2.4) 2.2 Fuerzas y momentos 7 donde Γ1=Ixz(Ix−Iy+Iz) Γ Γ2=Iz(Iz−Iy)+I2 xz Γ Γ3=Iz Γ Γ4=Ixz Γ Γ5=Iz−Ix Iy Γ6=Ixz Iy Γ7=Ix(Ix−Iy)+I2 xz Γ Γ8=Ix Γ (2.5) siendo Γ = IxIz−I2 xz 2.2 Fuerzas y momentos En esta sección se van a desarrollar las hipótesis y ecuaciones que dan lugar al modelo de fuerzas y momentos que experimenta la aeronave, que dan lugar al cierre del sistema que modela nuestro UAV. Para modelar la fuerza total que experimenta la aeronave se puede hacer uso del Principio de Superposición y descomponer la fuerza resultante como la suma de tres componentes: fuerza gravitatoria, aerodinámica y propulsiva. X~ F=~ Fgravitatoria +~ Faerodinámica +~ Fpropulsiva (2.6) Para el momento resultante que experimenta la aeronave se vuelve a hacer uso de nuevo del Principio de Superposición, considerando el momento total como el producido por la suma del momento producido por la fuerza aerodinámico y propulsivo. En este caso, la fuerza gravitatoria no produce momento en la aeronave, ya que ésta se considera como una fuerza puntual aplicada en el centro de gravedad de la aeronave. X~ M=~ Maerodinámico +~ Mpropulsivo (2.7) 2.2.1 Modelo gravitatorio Para modelar el gradiente gravitatorio al que se encuentra sometida la aeronave se va a hacer uso del modelo de gravedad Earth Gravitational Model 96 (EGM96), el cual se encuentra incluido en el sistema de coordenadas World Geodetic System 84 (WGS84). Dicho modelo proporciona un valor del gradiente gravitatorio en función de la latitud (ϕ)y altitud de vuelo. Para el modelo matemático empleado por el EGM96 se ha hecho uso del empleado en [3] g(r,ϕ) = −GMe r2"1−3Iya2 2r2(3sin2(ϕ)−1)#+ω2 ercos2(ϕ)(2.8) siendo Gla constante de gravitación universal, Me la masa del esferoide por el que se aproxima la Tierra, a el semieje mayor y ωela velocidad angular de la Tierra. 8Capítulo 2. Modelado Figura 2.2 Esferoide del sistema WGS84.. Considerando que la fuerza gravitatoria se encuentra aplicada en el centro de masas de la aeronave y expresándola en ejes de navegación, se obtiene el siguiente vector de fuerza gravitatoria ~ Fg=  −mgsin(θ) mgcos(θ)sin(φ) mgcos(θ)cos(φ) (2.9) Como ya se comentó anteriormente, debido a que la fuerza gravitatoria se encuentra aplicada en el centro de gravedad de la aeronave, ésta no produce momento. 2.2.2 Modelo aerodinámico Para el desarrollo del modelo aerodinámico de la aeronave se va a considerar que la fuerzas de presión distribuidas a lo largo del perfil del ala van a estar aplicadas como una fuerza puntual en el centro aerodinámico de la aeronave, dando lugar así a la fuerza de sustentación (~ Flift) , resistencia aerodinámica (~ Fdrag) y el momento aerodinámicos (~ Ma). Figura 2.3 Fuerzas y momento aerodinámicos aplicadas en el centro aerodinámico del perfil.. Dichas fuerzas se modelan mediante las siguientes expresiones ~ Flift =1 2ρ~ Va 2SCL(2.10) ~ Fdrag =1 2ρ~ Va 2SCD(2.11) ~m =1 2ρ~ Va 2ScCm(2.12) donde CL , CD y Cm son los coeficientes aerodinámicos, S la superficie alar y c la cuerda media del perfil de 2.2 Fuerzas y momentos 9 la aeronave. Los coeficientes aerodinámicos dependen de las superficies de control de la aeronave: alerones, elevadores y timón de profundidad, que su configuración dependerá del tipo de aeronave estudiada. Superficies de control Antes de desarrollar el modelo que define los coeficientes aerodinámicos de la aeronave en función de las superficies de control, se va a hacer una pequeña introducción sobre diferentes configuraciones de la aeronave. Configuración estándar La configuración estándar de una aeronave se puede ver en la Fig. 2.4, en la cual actúan como superficies de control los alerones (δa), elevadores (δe)y timón de profundidad (δr). La deflexión del alerón se puede considerar como una deflexión compuesta formada por δa=1 2δa,left −δa,right(2.13) Figura 2.4 Variables de control aerodinámicas en configuración estándar.. Para aeronaves pequeñas existen dos configuraciones alternativas: cola en V y ala volante. Configuración con cola en V Este tipo de configuración sustituye el timón de profundidad (rudder) y elevadores (elevators) por los “ruddervators”. La deflexión angular del “ruddervator” derecho se denota como (δrr) , y la deflexión angular del “ruddervator” izquierdo se denota como (δrl) . Deflectarlos en distinto sentido produce el mismo efecto que el timón de profundidad, mientras que deflectarlos en el mismo sentido produce igual efecto que los elevadores. Matemáticamente se puede convertir entre “ruddervators” y timón-elevadores mediante δe δr=1 1 −1 1δrr δrl (2.14) Ala volante En este tipo de configuración (Fig. 2.6) los coeficientes aerodinámicos dependen de los elevones. La deflexión angular del elevón derecho se denota como δer , y la deflexión angular del elevón izquierdo se denota como δel . Deflectar los elevones en distinto sentido tiene el mismo efecto que los alerones (δa) , mientras que deflectar los elevones juntos, tiene el mismo efecto que los elevadores (δe). 16 Capítulo 2. Modelado donde σu , σv y σw son las intensidades de la turbulencia a lo largo de los ejes de la aeronave, Lu , Lv y Lw son longitudes de onda espaciales y Va la velocidad de la aeronave. Los modelos Dryden se implementan típicamente asumiendo una velocidad nominal constante Va0 . Los parámetros para el modelo de ráfaga Dryden se definen en MIL-F-8785C. A continuación, se muestra una tabla extraída de [ 6 ] donde se presentan parámetros apropiados para las distintas condiciones de turbulencia deseadas. Tabla 2.3 Parámetros de ráfaga del modelo Dryden [6]. Tipo de Altitud Lu=LvLwσu=σvσw ráfaga (m) (m) (m) (m/s) (m/s) altitud baja, turbulencia débil 50 200 50 1.06 0.7 altitud baja, turbulencia moderada 50 200 50 2.12 1.4 altitud media, turbulencia débil 600 533 533 1.5 1.5 altitud media, turbulencia moderada 600 533 533 3.0 3.0 Las componentes constantes del viento se giran desde el sistema de referencia inercial al sistema de ejes cuerpo y se agregan a las componentes de ráfaga para producir el viento total en ejes cuerpo. La combinación de términos estables y ráfagas se puede expresar matemáticamente como Vb w=  uw vw ww =Rb v(φ,θ,ψ)  wns wes wds  +  ugs vgs wws  (2.49) A partir de las componentes de la velocidad Vb w y la velocidad respecto a tierra Vb g , se puede calcular las componentes del vector velocidad del aire como Vb a=  ur vr wr =  u−uw v−vw w−ww (2.50) Partiendo de las componentes del vector velocidad del aire Vb a se puede calcular el valor de la velocidad Va , el ángulo de ataque αy ángulo de resbalamiento β Va=qu2 r+v2 r+w2 r(2.51) α= tan−1(wr ur )(2.52) β= sin−1 vr pu2 r+v2 r+w2 r!(2.53) Estas expresiones de Va , α y β se emplean para calcular las fuerzas y momentos aerodinámicos sobre la aeronave. 2.3 Resumen La fuerza y momento total que actúa sobre la aeronave se puede escribir como "fx fy fz#="−mg sin(θ) mg cos(θ)sin(φ) mg cos(θ)cos(φ)# +1 2ρV 2 aS  CX(α)+CXq(α)c 2Vaq+CXδeδe CY0+CYββ+CYp b 2Vap+CYr b 2Var+CYδaδa+CYδrδr CZ(α)+CZq(α)c 2Vaq+Czδeδe  +1 2ρSpropCprop "(kmotorδt)2 −V2 a 0 0# (2.54) 2.4 Trimado UAV 17   l m n =1 2ρV 2 aS     bCl0+Clββ+Clp b 2Vap+Clr b 2Var+Clδaδa+Clδrδr cCm0+Cmaα+Cmq c 2Vaq+Cmδeδe bCn0+Cnββ+Cnp b 2Vap+Cnr b 2Var+Cnδaδa+Cnδrδri      (2.55) donde CX(α) = −CD(α)cos(α) +CL(α)sin(α) CXq(α) = −CDq(α)cos(α) +CLq(α)sin(α) CXδe(α) = −CDδecos(α)+CLδesin(α) CZ(α) = −CD(α)sin(α)+CL(α)cos(α) CZq(α) = −CDqsin(α)−CLqcos(α) CZδe(α) = −CDδesin(α)−CLδecos(α) (2.56) Notas y referencias Las ecuaciones desarrolladas en este capítulo han sido extraídas de [2] 2.4 Trimado UAV En este apartado vamos a calcular las condiciones de trimado para el UAV. Este cálculo es indispensable para la posterior implementación de controladores en el UAV, ya que es necesario obtener las variables de referencia necesarias para poder realizar el control en bucle cerrado, tal y como se desarrolla en el Capítulo 3. Vamos a supooner que el UAV va a realizar un vuelo a 100m de altitud (por simplificación consideraremos que el territorio a sobrevolar se va a encontrar a nivel del mar). Las ecuaciones necesarias para calcular el valor de las variables de control δt , δa , δe van a ser una simplificación de las ecuaciones desarrolladas durante el Capítulo, ya que permiten un cálculo muy sencillo y a la vez muy aproximado sobre el estado de trimado del UAV. Las ecuaciones que nos permiten calcular dicho punto son las correspondientes a vuelo de crucero y para velocidad de 15 m/s, a partir de las cuales se obtienen los siguientes valores: δ∗ t= 0.34875 δ∗ a=−4.1359·10−25 rad ∼0rad δ∗ e=−0.0072586 rad 3 Control y Guiado: Diseño de Autopiloto En términos generales, un piloto automático es un sistema utilizado para guiar una aeronave sin la ayuda de un piloto. Para los UAV, el piloto automático tiene el control completo de la aeronave durante todas las fases del vuelo. Si bien algunas funciones de control pueden residir en la estación de control de tierra, la parte del piloto automático del sistema de control de UAV reside a bordo del UAV. En este capítulo se diseñará un autopiloto para UAVs basado en control mediante bucle cerrado. El objetivo de este diseño es implementar una ley de control que se encargue de calcular en todo momento el valor de las variables de control del sistema ( δt,δa,δe ) para que la aeronave vuele de manera autónoma. El hecho de que la variable de control δr no se haya mencionado es debido a que el UAV estudiado se trata de un ala volante, y por tanto, carece de timón de profundidad. Para la implementación de los controladores se va a hacer uso del Capítulo 6 de [2]. 3.1 Autopiloto Longitudinal El objetivo en el diseño del piloto automático longitudinal será regular la velocidad del aire y la altitud utilizando la señal PWM ( δt ) y elevadores ( δe ) como actuadores. El método utilizado para regular la altitud y la velocidad del aire depende del error de altitud. Los regímenes de vuelo se muestran en la figura 3.1 Figura 3.1 Regímenes de vuelo para el autopiloto longitudinal.. En la zona de despegue, se ordena empuje máximo ( δt= 1 ) y la actitud de cabeceo se regula a un ángulo de cabeceo fijo θc utilizando los elevadores. El objetivo en la zona de ascenso es maximizar la velocidad 19 20 Capítulo 3. Control: Diseño de Autopiloto de ascenso dadas las condiciones atmosféricas actuales. Para maximizar la velocidad de ascenso, se ordena empuje máximo y la velocidad del aire se regula utilizando el ángulo de cabeceo. La zona de descenso es similar a la zona de ascenso, excepto que la señal PWM se ordena a cero. En la zona de retención de altitud, la velocidad del aire se regula ajustando la señal PWM, y la altitud se regula mediante el ángulo de cabeceo. 3.1.1 Controlador de Cabeceo A partir de la Fig. (3.2) se puede obtener la función de transferencia que permite pasar de θcaθ Figura 3.2 Diagrama de bloques del controlador de cabeceo.. Hθ/θc(s) = kpθaθ3 s2+aθ1+kdθaθ3s+aθ2+kpθaθ3(3.1) Si queremos expresar la función de transferencia en la forma canónica H(s) = KθDC ω2 nθ s2+2ζθωnθs+ω2 nθ (3.2) Igualando términos se obtiene ω2 nθ=aθ2+kpθaθ3(3.3a) 2ζθωnθ=aθ1+kdθaθ3(3.3b) Si configuramos la ganancia proporcional para evitar la saturación cuando se experimenta el error máximo de entrada, obtenemos kpθ=δmax e emax θ sign(aθ3)(3.4) donde se toma el signo de aθ3 , ya que aθ3 se basa en Cmδe que suele ser negativo. Para garantizar la estabilidad, kpθ y aθ3 deben ser del mismo signo. A partir de la ecuación (3.3a) se puede calcular el límite del ancho de banda del bucle de cabeceo como ωnθ=saθ2+δmax e emax θ |aθ3|(3.5) y resolviendo la ecuación (3.3b) kdθ=2ζθωnθ−aθ1 aθ3 (3.6) En resumen, conociendo el límite de saturación del actuador δmax e y el máximo error de cabeceo, se puede seleccionar emax θ para determinar la ganancia proporcional kpθ y el ancho de banda del bucle de cabeceo. Seleccionando la relación de amortiguamiento deseada ζθse fija el valor de ganancia derivativa kdθ. 3.1 Autopiloto Longitudinal 21 La ganancia DC viene dada por KθDC =kpθaθ3 aθ2+kpθaθ3(3.7) Por tanto, se realiza el control mediante bucle cerrado con un controlador PD, siendo la respuesta de los elevadores δe=kpθ(θc−θ)−kdθq(3.8) No se ha empleado el término integral debido a que puede limitar considerablemente el ancho de banda del bucle interno. 3.1.2 Controlador de Altitud mediante Cabeceo El controlador de altitud se basa en un controlador PI en el que a partir de la entrada de altitud de referencia hc , se obtiene el ángulo de cabeceo de referencia θc entrante al controlador de cabeceo. El diseño del controlador de altitud se encuentra en la Fig. (3.3). De forma análoga a como se ha realizado en el controlador de cabeceo, Figura 3.3 Diagrama de bloques del controlador de altitud.. obtenemos las siguientes expresiones h(s) =   KθDC Vakphs+kih kph s2+KθDC Vakphs+KθDC Vakih hc(s) + s s2+KθDC Vakphs+KθDC Vakih!dh(s) (3.9) Las ganancias kph y kih deben elegirse de modo que el ancho de banda del bucle de altitud-cabeceo sea menor que el ancho de banda del bucle de actitud-cabeceo. Calculamos la frecuencia natural a partir de la separación del ancho de banda Wh ωnh=ωnθ Wh (3.10) Volviendo a hacer uso de la ecuación (3.2) calculamos los coeficientes ω2 nh=KθDC Vakih(3.11a) 2ζhωnh=KθDC Vakph(3.11b) y despejamos el valor de las ganancias kphykih kih=ω2 nh KθDC Va (3.12) 22 Capítulo 3. Control: Diseño de Autopiloto kph=2ζhωnh KθDC Va (3.13) Por lo tanto, al seleccionar la relación de amortiguamiento deseada ζh y la separación del ancho de banda Whse fija el valor para kphykih. La salida del bucle cerrado de altitud-cabeceo es θc=kph(hc−h)+ kih s(hc−h)(3.14) 3.1.3 Controlador de Velocidad mediante Cabeceo El controlador de velocidad mediante el ángulo de cabeceo se basa en un controlador PI en el que a partir de la entrada de velocidad de referencia Vc a , se obtiene el ángulo de cabeceo de referencia θc entrante al controlador de cabeceo. El diseño de dicho controlador se encuentra en la Fig. (3.4). Figura 3.4 Diagrama de bloques del controlador de velocidad mediante cabeceo.. En el dominio de Laplace Va(s) =     −KθDC gkpV2s+kiV2 kpV2 s2+(aV1−KθDC gkpV2)s−KθDC gkpV2     Vc a(s) + s s2+(aV1−KθDC gkpV2)s−KθDC gkpV2!dV(s) (3.15) Las ganancias kpV2 y kiV2 deben elegirse de modo que el ancho de banda del bucle de altitud-cabeceo sea menor que el ancho de banda del bucle de actitud-cabeceo. Calculamos la frecuencia natural a partir de la separación del ancho de banda WV2 ωnV2=ωnθ WV2 (3.16) Se vuelven a calcular las ganancias kiV2=−ω2 nV2 KθDC g(3.17) kpV2=aV1−2ζV2ωnV2 KθDC g(3.18) 3.1 Autopiloto Longitudinal 23 Seleccionando la relación de amortiguamiento deseada ζV2 y la separación del ancho de banda WV2 se fija el valor para kpV2ykiV2. La salida del bucle cerrado velocidad-cabeceo es θc=kpV2(Vc a−Va)+ kiV2 s(Vc a−Va)(3.19) 3.1.4 Controlador de Velocidad mediante Empuje El controlador de velocidad mediante empuje se basa en un controlador PI en el que a partir de la entrada de velocidad de referencia Vc a , se obtiene el valor de la señal PWM δt . El diseño de dicho controlador se encuentra en la Fig. (3.5). En el dominio de Laplace Figura 3.5 Diagrama de bloques del controlador de velocidad mediante empuje.. Va(s) = aV2(kpVs+kiV ) s2+(aV1+aV2kpV)s+aV2kiV!Vc a(s) + 1 s2+(aV1+aV2kpV)s+aV2kiV!dV(s) (3.20) Si aV1 y aV2 son conocidas, las ganancias kpV y kiV se obtienen mediante las mismas técnicas comentadas anteriormente kiV=ω2 nV aV2 (3.21) kpV=2ζVωnV−aV1 aV2 (3.22) Los parámetros de diseño para este bucle cerrado son la relación de amortiguamiento ζV y frecuencia natural ωnV. La salida del bucle cerrado velocidad-empuje es δt=δ∗ t+kpV(Vc a−Va)+ kiV s(Vc a−Va)(3.23) donde δ∗ tes el valor de la señal PWM en condición de trimado del UAV. 3.1.5 Máquina de Estados de Control de Altitud En este apartado se va a desarrollar el bloque correspondiente a la máquina de estados de control longitudinal. Se trata de un bloque necesario para poder implementar en el autopiloto los regímenes de la Fig. (3.1). En función de la condición de vuelo que se encuentre el UAV, la máquina de estados se encarga de indicar qué controladores deben actuar y cuales no. Podemos ver una representación de dicha máquina de estados en la Fig. (3.6). Se observa que cada régimen de vuelo se encuentra caracterizado por un ley de control: 24 Capítulo 3. Control: Diseño de Autopiloto •Zona de descenso: empuje nulo y control de velocidad mediante cabeceo. •Zona de retención de altitud : control de velocidad mediante empuje y control de altitud mediante cabeceo. •Zona de ascenso: empuje máximo y control de velocidad mediante cabeceo. •Zona de aterrizaje: empuje máximo y ángulo de cabeceo fijo. Figura 3.6 Máquina de estados para control longitudinal.. 3.2 Autopiloto Lateral El objetivo del autopiloto lateral es la de obtener una ley de control que se encargue de la deflexión automática de las superficies de control lateral. Al tratarse de un ala volante, sólo tendremos control mediante los alerones δa . Por esta razón, es necesario implementar una ley de control en cascada que permita obtener la deflexión de los alerones a partir de un ángulo de guiñada de referencia. Para ello es necesario diseñar un controlador de alabeo cuya entra de referencia es la salida del controlador de guiñada, es decir, el control del ángulo de guiñada y alabeo se encuentra acoplado. 3.2.1 Controlador de Alabeo El controlador de alabeo se basa en un controlador PID en el que a partir de la entrada de alabeo de referencia φc , se obtiene la deflexión de los alerones δa necesaria. El diseño de dicho controlador se encuentra en la Fig. (3.7). 3.2 Autopiloto Lateral 25 Figura 3.7 Diagrama de bloques del controlador de alabeo.. En el dominio de Laplace φ= s s3+(aφ1+aφ2kdφ)s2+aφ2kpφs+aφ2kiφ!dφ2 +    aφ2kpφs+kiφ kpφ s3+(aφ1+aφ2kdφ)s2+aφ2kpφs+aφ2kiφ     φc (3.24) De forma análoga a como se hizo en el controlador de cabeceo, obtenemos los parámetros que definen el controlador kpφ=δmax a emax φ sign(aφ2)(3.25) ωnφ=sδmax a emax φ |aφ2|(3.26) kdφ=2ζφωnφ−aφ1 aφ2 (3.27) donde la relación de amortiguamiento ζφes un parámetro de diseño. En cuanto al término integral del controlador, éste puede ser seleccionado empleando el método del lugar de las raíces [2]. La salida del bucle cerrado de alabeo es δa=kpφ(φc−φ)+ kiφ s(φc−φ)−kdφp(3.28) 3.2.2 Controlador de Rumbo Este controlador tiene como función calcular la entrada necesaria al controlador de alabeo φc para obtener el rumbo deseado debido a que el UAV no dispone de timón de profundidad. Se va emplear un controlador PI cuya entrada es el rumbo de referencia χc , cuyo modelo de diagrama de bloques se puede ver en la Fig. (3.8). Volviendo a elegir el valor de la frecuencia natural ωnχ y la relación de amortiguamiento ζχ podemos calcular el valor de las ganancias kpχ y kiχ . Calculamos el valor de la frecuencia natural a partir de la separación del ancho de banda Wχ. ωnχ=ωnφ Wχ (3.29) 32 Capítulo 4. Guiado: Sistema de Misión En la Fig. (4.6) podemos observar un ejemplo del algoritmo programado. 0 100 200 300 400 500 600 700 800 900 0 200 400 600 800 1000 1200 1400 1600 1800 Algoritmo TSP Figura 4.6 Algoritmo TSP.. 5 Implementación de Sensores y Filtro de Kalman Extendido Para llevar a cabo una simulación de mayor interés se va implementar en el modelo un observador basado en el Filtro de Kalman. Para ello es necesario dotar al UAV de una serie de sensores, los cuales se encargan de obtener las medidas necesarias para poder modelar el observador. Al tratarse de un sistema no lineal discreto, se hace uso del Filtro de Kalman Extendido, el cual se explicará en mayor detalle más adelante. 5.1 Sensores a bordo Los sensores introducidos en el UAV corresponden a una serie de 5 sensores básicos de los cuales dispone cualquier aeronave: acelerómetro, giróscopo, altímetro, anemómetro y magnetómetro. Estos sensores proporcionan información necesaria en términos de actitud y sistemas de control de navegación. Se encuentran caracterizados por el tiempo de muestreo ( Ts ), que se ha considerado de 0.2 segundos en cada uno de ellos, y los parámetros de ruido que introducen en la medición. Los parámetros utilizados para modelarlos han sido extraidos del Apéndice H de [2]. 5.1.1 Acelerómetro Este sensor se encarga de medir las aceleraciones que experimenta la aeronave. Uno de los modelos más empleados es el acelerómetro piezoeléctrico, el cual se encarga de convertir el desplazamiento de una masa de prueba en una señal eléctrica. En aeronaves se suelen emplear tres acelerómetros colocados cerca del centro de masa y alineados con los ejes cuerpo de la aeonave. El acelerómetro mide la aceleración del UAV en m/s2 , siendo el modelo matemático empleado el siguiente yaccel,x = ˙u+qw −rv +gsinθ+ηaccel,x (5.1a) yaccel,y = ˙v+ru −pw −gcosθsinφ+ηaccel,y (5.1b) yaccel,z = ˙wpvu−qu −gcosθcosφ+ηaccel,z (5.1c) donde el ruido presente en los sensores es modelado mediante los parámetros ηaccel,x , ηaccel,y y ηaccel,z , que corresponden a ruido gaussiano con covarianzas σ2 accel,x , σ2 accel,y y σ2 accel,z . Para la simulación realiza el acelerómetro que se ha utilizado corresponde al acelerómetro Analog Devices ADXL325. 5.1.2 Giróscopo Un giróscopo es un sensor mecánico basado en el Principio de Coriolis, que mide la velocidad angular de la aeronave con respeco a un sistema de referencia inercial. Es un aparato en el cual una masa que gira velozmente alrededor de su eje de simetría permite mantener de forma constante su orientación respecto a un sistema de ejes de referencia. El rápido movimiento giratorio del rotor de los giróscopos se puede obtener 33 34 Capítulo 5. Sensores y EKF por vacío o por un sistema eléctrico. El modelo matemático del giróscopo es el siguiente: ygyro,x =p+ηgyro,x (5.2a) ygyro,y =q+ηgyro,y (5.2b) ygyro,z =r+ηgyro,z (5.2c) donde el ruido presente en los sensores es modelado mediante los parámetros ηgyro,x , ηgyro,y y ηgyro,z , que corresponden a ruido gaussiano con covarianzas σ2 gyro,x , σ2 gyro,y y σ2 gyro,z . Para la simulación realiza el acelerómetro que se ha utilizado corresponde al acelerómetro Analog Devices ADXRS450. 5.1.3 Altímetro barométrico Un altímetro barométrico es un sensor que se encarga de obtener la altitud de vuelo a partir de la medida de presión atmosférica. Un sensor de presión absoluta se encarga de medir la presión de la atmósfera, utilizando para ello el Modelo de Atmósfera ISA Internacional. El modelo matemático empleado es el siguiente: Pbarométrica [Pa] = P0T0 T0+L0hASL gM RL0(5.3a) yabs pres = (P0−Pbarométrica)+βabs pres +ηabs pres (5.3b) donde el ruido presente en los sensores es modelado mediante el parámetro ηabs pres , que corresponde a ruido gaussiano con covarianza σ2 abs pres. El término βabs pres corresponde a un offset de temperatura. La expresión (5.4) permite calcular el valor de altura en metros sobre un punto de referencia a partir de la salida del sensor de presión absoluta: hbarométrica [m]=0.3048·(1−(Pbarométrica/P0)0.19026)·288.15 0.00198122 (5.4) 5.1.4 Anemómetro El anemómetro es un disposito que se encarga de medir la velocidad del aire a partir de un tubo de pitot con un sensor de presión diferencial. Este sensor de presión diferencial mide la diferencia entre la presión absoluta y presión estática y mediante señales eléctricas es capaz de obtener el valor de la presión diferencia. El modelo matemático empleado es: ydiff pres =ρV 2 a 2+βdiff pres +ηdiff pres (5.5) donde el ruido presente en los sensores es modelado mediante el parámetro ηdiff pres , que corresponde a ruido gaussiano con covarianza σ2 diff pres . El término βdiff pres corresponde a un offset de temperatura. A partir de la salida del sensor de presión diferencial se calcula la velocidad del aire. Para la simulación se ha utilizado al anemómetro Freescale Semiconductor MPXV5004G. 5.1.5 Magnetómetro Un magnetómetro mide la intensidad y dirección del campo magnético en el sistema de referencia ejes cuerpo. Si se conoce el campo magnético de la tierra en ejes fijos, se puede obtener la actitud comparándola con la medida del instrumento. El modelo matemático empleado es el siguiente: φmag =φ+βmag +ηmag (5.6a) θmag =θ+βmag +ηmag (5.6b) 5.2 Sistema GPS 35 ψmag =ψ+βmag +ηmag (5.6c) donde el ruido presente en los sensores es modelado mediante el parámetro ηmag , que corresponde a ruido gaussiano con covarianza σ2 mag . El término βmag corresponde a un error de bias. Para la simulación se ha utilizado el Magnetómetro Honeywell HMR3300. 5.2 Sistema GPS El Sistema de posicionamiento global (GPS) es un sistema de navegación por satélite que proporciona información acerca de la posición de un determinado cuerpo en la superficie terrestre. Está formado por una constelación de 24 satélites que orbitan continuamente alrededor de la Tierra a una altitud de 20180 km. Al medir los tiempos de vuelo de las señales desde un mínimo de cuatro satélites a un receptor en la superficie de la Tierra o cerca de ella, se puede determinar la ubicación del receptor en tres dimensiones. La precisión de una medición de posición GPS se ve afectada por la precisión de las mediciones de pseudodistancia del satélite y por la geometría de los satélites a partir de las cuales se toman las mediciones de pseudodistancia. La precisión de pseudodistancia se ve afectada por errores en el tiempo de medición de vuelo para cada satélite. Dado que las señales de radio electromagnéticas de los satélites viajan a la velocidad de la luz, pequeños errores de tiempo pueden causar errores de posicionamiento significativos. Por ejemplo, un error de tiempo de solo 10 ns puede resultar en un error de posicionamiento de unos 3 m. Para los propósitos de la simulación también hay que tener en cuenta la dinámica del error cometido por el Sistema GPS. Los errores de norte, este y altitud se componen de un sesgo que cambia lentamente junto con ruido aleatorio. Para modelar el comportamiento transitorio del error, seguimos el enfoque de [37] y modelamos el error como un proceso de Gauss-Markov. Los procesos de Gauss-Markov son modelados por v[n+1] = e−kGP STSv[n]+ηGP S [n](5.7) donde donde v[n] es el error que se está modelando, ηGP S[n] es ruido blanco Gaussiano de media cero, 1/kGP S es la constante de tiempo del proceso y TS es el tiempo de muestreo. Figura 5.1 Parámetros del modelo de error de Gauss-Markov.. Un modelo para mediciones GPS que es adecuado para propósitos de simulación es dado por yGP S,n[n] = pn[n]+vn[n](5.8a) yGP S,e[n] = pe[n]+ve[n](5.8b) yGP S,h[n] = −pd[n]+vh[n](5.8c) donde pn , pe es la posición real del UAV en la Tierra y h la altitud sobre el nivel del mar, y n es el índice de muestra. Mediante el efecto Doppler presente en las señales de satélite GPS, la velocidad del receptor se puede calcular a precisiones con desviaciones estándar en el rango de 0.01 a 0.05 mx/s. Hay que tener en cuenta que 36 Capítulo 5. Sensores y EKF la incertidumbre en la medición del rumbo varía con la inversa de la velocidad de avance: para velocidades altas, el error es pequeño y para velocidades bajas el error es grande. Podemos modelar la velocidad del terreno y las medidas del recorrido disponible desde GPS como yGP S,Vg=q(Vacos(ψ)+ wn)2+(Vasin(ψ)+ we)2+ηV(5.9a) yGP S,χ =atan2(Vasin(ψ)+we,Vacos(ψ)+wn)+ηχ(5.9b) donde wn y we es la velocidad del viento en ejes NED y los parámetros ηV y ηχ son ruido gaussiano con varianzas σ2 Vgyσ2 χ. 5.3 Filtro de Kalman Extendido El sistema que hemos modelado corresponde a un sistema no lineal continuo en el tiempo representado por ecuaciones diferenciales ordinarias, mientras que los sensores que hemos desarrollado son sistemas que toman medidas en intervalos de tiempo Ts , es decir, se trata de un sistema discreto, Por tanto, el observador que vamos a implementar estará definido en tiempo discreto. Un sistema no lineal puede escribirse de la forma ˙x=f(x,u,t)+ζ y[n] = h(x[n],u[n])+η[n] donde x representa el estado del sistema, u la entrada del sistema y ζ , η perturbaciones modeladas como ruido blanco. En este apartado explicaremos cómo sabiendo el estado del sistema dado por las ecuaciones y los modelos anteriores, conocidas las entradas del sistema y las medidas (sensores) cada cierto instante de tiempo, podemos obtener una estimación del sistema, denotado por ˆx . Para ello hacemos uso del Filtro de Kalman: conocida u(t) y una estimación inicial ˆx0 , intentar reconstruir las variables de estado conociendo el modelo matemático y las medidas proveniente de los sensores cada cierto intervalo de tiempo. Al tratarse de un sistema no lineal es neceario hacer uso del Filtro de Kalman Extendido (EKF), que permite hacer estimaciones de sistemas discretos no lineales. Supongamos que las ecuaciones de transición de estado y de medición para un sistema no lineal de tiempo discreto tienen términos de proceso y ruido de medición no aditivos con media cero y matrices de covarianza Q y R, respectivamente: x[k+1] = f(x[k],w[k],us[k]) y[k] = h(x[k],v[k],um[k]) w[k]∼(0,Q[k]) v[k]∼(0,R[k]) El algoritmo que se encarga de implementar el EKF es el siguiente: 1. Inicializar el objeto de filtro con los valores iniciales del estado, x[0] , y la matriz de covarianza del error de estimación del estado, P. ˆx[0|−1] = E(x[0]) P[0|−1] = E(x[0]−ˆx[0|−1])(x[0]−ˆx[0| −1])T donde ˆx es el estado estimado y ˆx[ka|kb] denota el estado estimado en el paso de tiempo ka usando medidas en instantes 0,1,...kb . Así, ˆx[0|−1] es la mejor estimación del valor del estado antes de realizar cualquier medición. Este valor es necesario proporcionarlo al comenzar el algoritmo. 2. Para pasos de tiempo k= 0,1,2,3... realizar los siguientes pasos: 5.3 Filtro de Kalman Extendido 37 a) Calcular el jacobiano de la función de medición y actualice la covarianza del error de estimación de estado y estado utilizando los datos medidos, y[k]. C[k] = ∂h ∂xˆx[k|k−1] S[k] = ∂h ∂v ˆx[k|k−1] Calcular estas matrices jacobianas numéricamente a menos que se especifique el jacobiano analítico. K[k] = P[k|k−1]C[k]TC[k]P[k|k−1]C[k]T+S[k]R[k]S[k]T−1 ˆx[k|k] = ˆx[k|k−1]+K[k](y[k]−h(ˆx[k|k−1],0,um[k])) ˆ P[k|k] = P[k|k−1]−K[k]C[k]P[k|k−1] donde Kes la ganancia de Kalman. b) Calcular el jacobiano de la función de transición de estado y la predicción de la covarianza del error de estimación de estado y el estado en el siguiente paso de tiempo. A[k] = ∂f ∂xˆx[k|k] G[k] = ∂f ∂w ˆx[k|k] Calcular estas matrices jacobianas numéricamente a menos que especifique el jacobiano analítico. P[k+1|k] = A[k]P[k|k]A[k]T+G[k]Q[k]G[k]T ˆx[k+ 1|k] = f(ˆx[k|k],0,us[k]) Los pasos del algoritmo descritos anteriormente suponen que tiene términos de ruido no aditivos en las funciones de transición de estado y medición. Si hay términos de ruido aditivo en las funciones, los cambios en el algoritmo son: • Si el ruido del proceso w es aditivo, es decir, la ecuación de transición de estado tiene la forma x[k] = f(x[k−1],us[k−1]+w[k−1]), entonces la matriz jacobiana G[k]es una matriz identidad. • Si el ruido de medición v es aditivo, es decir, la ecuación de medición tiene la forma y[k] = h(x[k],um[k]+v[k]) , entonces la matriz jacobiana S[k]es una matriz identidad. 5.3.1 Filtro de Kalman Extendido para Posición Aplicamos el EKF desarrollado en este apartado para poder realizar una estimación sobre la posición del UAV a partir de las medidas proporcionadas por el Sistema GPS. Las ecuaciones que permiten conocer la derivada del vector posición en ejes NED viene dado por (2.1). Expresando estas ecuaciones en la forma necesaria para aplicar el EKF: pn[k+1] = pn[k]+hu[k]cθ[k]cψ[k]+v[k](sφ[k]sθ[k]cψ[k]−cφ[k]sψ[k])+ w[k](cφ[k]sθ[k]cψ[k]+sφ[k]sψ[k])]dt (5.12a) pe[k+1] = pe[k]+hu[k]cθ[k]sψ[k]+v(sφ[k]sθ[k]sψ[k]+cφ[k]cψ[k])+ w[k](cφ[k]sθ[k]sψ[k]−sφ[k]cψ[k])]dt (5.12b) pd[k+1] = pd[k]+h−u[k]sθ[k]+v[k](sφ[k]cθ[k])+w[k](cφ[k]cθ[k])idt (5.12c) 38 Capítulo 5. Sensores y EKF donde se ha empleado la notación cα≡cos(α) , sα≡sin(α) . El vector de estados corresponde a x= (pn,pe,pd)T y las entradas del sistema u= (u,v,w,φ,θ,ψ)T . Las medidas y[k] del sistema GPS vienen dadas por la expresión (5.8). 5.3.2 Filtro de Kalman Extendido para Actitud Aplicamos el EKF desarrollado en este apartado para poder realizar una estimación sobre la actitud del UAV a partir de las medidas proporcionadas por el magnetómetro. Las ecuaciones que permiten conocer la derivada del vector posición en ejes NED viene dado por (2.2). Expresando estas ecuaciones en la forma necesaria para aplicar el EKF: φ[k+1] = φ[k]+[p[k]+q[k]sin(φ[k])tan(θ[k]) +rcos(φ[k])tan(θ[k])]dt (5.13a) θ[k+1] = θ[k]+[q[k]cos(φ[k])−r[k]sin(φ[k])]dt (5.13b) ψ[k+1] = ψ[k]+[q[k]sin(φ[k])sec(θ[k])+r[k]cos(φ[k])sec(θ[k])]dt (5.13c) donde el vector de estados corresponde a x= (φ,θ,ψ)T y las entradas del sistema u= (p,q,r)T . Las medidas y[k]del magnetómetro vienen dadas por la expresión (5.6). 6 Implementación en Matlab-Simulink En este capítulo se van a realizar las simulaciones correspondientes al UAV modelado en los capítulos anteriores mediante la herramienta Matlab-Simulink. Antes de ir a la presentación de resultados se va a hacer una breve referencia sobre los distintos bloque empleados en Simulink y cómo han sido configurados. 6.1 Bloque de Cinemática y Dinámica Este bloque se encarga del cálculo de la cinemática y dinámica del UAV a partir de las ecuaciones planteadas en el Capítulo 2.1. 6.1.1 Cinemática Para el cálculo de la cinemática del UAV se ha hecho uso de dos bloques: uno que se encarga de hallar las componentes de la velocidad lineal y otro de la velocidad angular. El bloque de cinemática lineal cuenta con 6 variables de entrada: u , v , w , φ , θ , ψ y 3 variables de salida: ˙pn , ˙pe , ˙pd . Para hallar la posición de la aeronave se hace pasar la salida del bloque por un integrador para así obtener la posición del UAV en cada instante de simulación. El bloque de cinemática angular cuenta con 5 entradas: p , q , r , φ , θ y 3 variables de salida: ˙ φ , ˙ θ , ˙ ψ . Haciendo pasar dicha salida por un integrador se obtiene la actitud del UAV. Figura 6.1 Bloque Cinemática lineal en Simulink.. Figura 6.2 Bloque Cinemática angular en Simulink.. 39 40 Capítulo 6. Implementación en Matlab-Simulink 6.1.2 Dinámica Para el cálculo de la cinemática del UAV se ha hecho uso de dos bloques: uno que se encarga de hallar las componentes de la aceleración lineal y otro de la aceleración angular. El bloque de aceleración lineal cuenta con 9 entradas: u , v , w , p , q , r , fx , fy , fz y 3 variables de salida: ˙u , ˙v,˙w. Haciendo pasar dicha salida por un integrador se obtiene la velocidad del UAV. El bloque de aceleración angular cuenta con 6 entradas: p , q , r , l , m , n y 3 variables de salida: ˙u , ˙v , ˙w . Haciendo pasar dicha salida por un integrador se obtiene la velocidad angular del UAV. Figura 6.3 Bloque Dinámica lineal en Simulink.. Figura 6.4 Bloque Dinámica angular en Simulink.. 6.2 Bloque de Fuerzas y Momentos En este apartado se muestran los bloque empleados para el cálculo de las fuerzas y momentos que actúan sobre el UAV. También se hará un breve comentario sobre la obtención de los coeficientes aerodinámicos y modelos gravitatorio y de viento. 6.2.1 Modelo gravitatorio Para el modelo de gravedad se ha hecho uso del bloque WGS84 de Simulink. Se encarga de obtener el valor de la gravedad a partir de la posición del UAV en coordenadas geodésicas. Por lo tanto, es necesario pasar de coordenadas NED a geodésicas antes de entrar al bloque. 6.2 Bloque de Fuerzas y Momentos 41 Figura 6.5 Bloque WGS84 de Simulink.. 6.2.2 Modelo aerodinámico Para el cálculo de los coeficientes aerodinámicos del UAV se hace uso de un bloque con 9 entradas: u , v , w , p,q,r,δa,δe,viento(uv,vv,wv,pv,qv,rv). Los coeficientes de salida son 6: CL,CD,CY,Cl,Cm,Cn. Figura 6.6 Bloque Coeficientes Aerodinámicos en Simulink.. 6.2.3 Viento Este bloque se encarga de modelar el viento presente mediante el modelo de viento continuo de Dryden. El sistema tiene 7 entradas: h , u , v , w , φ , θ , ψ y como salidas las componentes del viento: uv , vv , wv , pv , qv , rvy la velocidad relativa del UAV respecto a la masa de aire circundante: ua,va,wa. Figura 6.7 Bloque de Modelado del Viento en Simulink.. 48 Capítulo 6. Implementación en Matlab-Simulink Figura 6.24 Altímetro barométrico en Simulink.. 6.5.4 Anemómetro Este bloque implementa el anemómetro desarrollado en el Capítulo 5.1.4 y consta de 2 entradas: ρ , Va . La salida del bloque corresponde a la presión dinámica. Figura 6.25 Anemómetro en Simulink.. 6.5.5 Magnetómetro Este bloque implementa el anemómetro desarrollado en el Capítulo 5.1.5 y consta de 3 entradas: φ , θ , ψ . La salida del bloque corresponde a la actitud del UAV: φmag,θmag,ψmag. Figura 6.26 Magnetómetro en Simulink.. 6.6 Sistema GPS Bloque que implementa el Sistema GPS desarrollado en el Capítulo 5.2. A partir de las entradas al bloque se encarga de calcular la posición, velocidad y rumbo del UAV. Para ello se ha introducido una función que se encarga de calcular las componentes del viento en ejes NED: (wn,we,wd)T. 6.7 Filtro de Kalman Extendido 49 Figura 6.27 Sistema GPS en Simulink.. 6.7 Filtro de Kalman Extendido Para la implementación del EKF se ha hecho uso del bloque “Extended Kalman Filter” de Simulink. Este bloque se encarga de aplicar el EKF a partir de una serie de parámetros de entrada: • State Transition: función f(x[k],w[k],us[k]) que modela la dinámica del UAV. Hay que indicar si el ruido es aditivo o no aditivo y proporcionar la matriz de covarianza. •Inizalitation: vector de estado inicial con su respectiva covarianza. • Measurement: función h(x[k],v[k],um[k]) que modela las medidas tomadas por el sensor. Hay que indicar si el ruido es aditivo o no aditivo y proporcionar la matriz de covarianza. •Sample Time: tiempo de muestreo. A partir de estos parámetros el bloque calcula las variables estimadas. Este bloque calcula el jacobiano de manera numérica, aunque se le puede introducir la expresión analítica en el caso de que se disponga de ella, lo cual hace que los resultados proporcionados sean más fiables y con menor tiempo de computación. 6.7.1 EKF Posición EL EKF diseñado para estimar la posición del UAV hace uso de las medidas tomadas por el Sistema GPS. Podemos ver el EKF implementado en la Fig. (6.28) Figura 6.28 EKF de Posición en Simulink.. 50 Capítulo 6. Implementación en Matlab-Simulink 6.7.2 EKF Actitud EL EKF diseñado para estimar la actitud del UAV hace uso de las medidas tomadas por el magnetómetro. Podemos ver el EKF implementado en la Fig. (6.29) Figura 6.29 EKF de Actitud en Simulink.. 7 Resultados 7.1 Misión realizada En este Capítulo se van a mostrar los resultados obtenidos para la misión deseada. Ésta consiste en un vuelo que abarca una superficie de unos 200 hectáreas en las proximidades del Aeropuerto de Sevilla (LEZL (OACI), SVQ (IATA)), en el que el UAV recorre 6025m. Para simplificar la simulación se va considerar sólo la fase correspondiente a vuelo horizontal, obviando el despegue y aterrizaje. Se pretende que el UAV parta de un punto inicial y vuelva a dicho punto pasando por una serie de waypoints definidos. La simulación corresponde a un vuelo con velocidad de trimado del UAV de 15 m/s y viento de intensidad 10 m/s y componente WSW. Además de la misión realizada por el UAV se van a mostrar otros resultados de interés, como la actitud de la aeronave durante toda la simulación, medidas tomadas por los sensores, la salida del EKF de posición y actitud, etc. 7.2 Trayectoria En la Fig. (7.1) se mustra la trayectoria descrita por el UAV, obtenida mediante Matlab-Simulink. En ella se han representado los waypoints por los que se quiere que pase el UAV, el camino óptimo obtenido mediante el algoritmo TSP (color azul) y la trayectoria descrita por el UAV (color naranja). Hay que indicar que se ha decidido utilizar como origen de coordenadas el punto de partida de la aeronave (0,0,100)T por el mero hecho de facilitar el cálculo en Matlab. Posteriormente se calculará la posición geodésica real de los waypoints y puntos que conforman la trayectoria descrita por el UAV. Podemos comprobar cómo el UAV pasa por cada uno de los waypoints definidos y que la trayectoria que describe es prácticamente igual al camino óptimo. Las mayores diferencias las encontramos en aquellas maniobras que requieren virajes más bruscos, debido a que la aeronave está diseñada para realizar maniobras suaves. Aún así podemos afirmar que el modelo de guiado para el UAV cumple las expectativas con los que se comienza el proyecto. Para dotar de mayor realismo a la simulación, en el Anexo A se muestra la misión llevada a cabo por el UAV en Google Earth. Se vuelve a representar los waypoints por los que se quiere que pase el UAV, el camino óptimo obtenido mediante el algoritmo TSP (color azul) y la trayectoria descrita por el UAV (color naranja). También se van a representar los resultados de posición obtenidos mediante la integración de las ecuaciones diferenciales (ecs. 2.1) y las medidas tomadas por el Sistema GPS y la estimación del EKF de posición. 51 52 Capítulo 7. Resultados -200 0 200 400 600 800 1000 Pe(Este)[m] 0 500 1000 1500 2000 2500 Pn(Norte)[m] Misión realizada por el UAV Camino óptimo Trayectoria UAV Figura 7.1 Misión realizada por el UAV en Matlab-Simulink.. Figura 7.2 Posición Norte del UAV.. 7.3 Actitud 53 Figura 7.3 Posición Este del UAV.. Figura 7.4 Altitud del UAV.. A la vista de los resultados, observamos como el EKF implementado para la estimación de la posición del UAV funciona correctamente. La mínima diferencia entre los valores reales y medidos/estimados se debe al error introducido por los sensores. 7.3 Actitud En este apartado se van a mostrar los resultados correspondientes a la actitud del UAV. Para ello vamos a hacer uso del valor de φ , θ , ψ calculados mediante la integración de las ecuaciones diferenciales (ecs. 2.2), las medidas tomadas por el magnetómetro y la estimación del EKF de actitud. Hay que indicar que los resultados obtenidos se esperan que tenga menor exactitud que los obtenidos anteriormente para posición, debido a que la aeronave se encuentra en todo momento cambiando su actitud para poder realizar la misión deseada. 54 Capítulo 7. Resultados Figura 7.5 Alabeo del UAV.. Figura 7.6 Cabeceo del UAV.. A la vista de los resultados, observamos como el EKF implementado para la estimación de la actitud del UAV funciona correctamente. La diferencia entre los valores reales y medidos/estimados se debe al error introducido por los sensores. 7.4 Velocidad 55 Figura 7.7 Guiñada del UAV.. 7.4 Velocidad En este apartado se muestran las 3 componentes de la velocidad del UAV durante la simulación. Comprobamos como la componente longitudinal es mayor que las demás componentes, algo totalmente lógico ya que corresponde a la dirección en la cual se produce la fuerza propulsiva. Ésta se encuentra en torno a 15 m/s , velocidad para la cual se había decidido trimar el UAV. 0 50 100 150 200 250 300 350 Tiempo (s) -10 -5 0 5 10 15 20 Velocidad (m/s) Velocidad UAV Velocidad longitudinal u Velocidad lateral v Velocidad vertical w Velocidad total Va Figura 7.8 Velocidad del UAV.. 7.5 Señal PWM de Empuje En este apartado mostramos la evolución de la Señal PWM encargada de modelar el empuje del UAV. En ella podemos ver cómo durante la fase de retención de altitud, se ha limitado la señal al intervalo [0.3,0.4] , valores dentro de los cuales se encuentra el valor de la Señal PWM para el trimado del UAV. También podemos 56 Capítulo 7. Resultados 0 50 100 150 200 250 300 350 Tiempo (s) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Señal PWM Señal PWM t Señal PWM Señal PWM trimado UAV Figura 7.9 Señal PWM encargada de modelar el empuje.. obsevar cómo cuando el UAV entra en zona de descenso y ascenso toma los valores 0 y 1, respectivamente, tal y como se explicó en el Capítulo 3.1.5. 7.6 Elevones En este apartado mostramos la deflexión de los elevones que actúan como alerones o elevadores del UAV. A partir de la ecuación (2.6) obtenemos la deflexión de los elevones derecho e izquierdo en función de la deflexión de alerones y elevadores obtenida en la simulación. δer =1 2(δe−δa)(7.1a) δel =1 2(δe+δa)(7.1b) 0 50 100 150 200 250 300 350 Tiempo (s) -8 -6 -4 -2 0 2 4 6 8 deflexión elevones (degrees) Elevones UAV Elevón derecho er Elevón izquierdo el Figura 7.10 Deflexión de los elevones del UAV.. 7.7 Coeficientes aerodinámicos 57 7.7 Coeficientes aerodinámicos En este apartado mostramos la evolución del coeficiente de sustentación y resistencia del UAV. Comprobamos como el CLes mucho mayor que el CD 0 50 100 150 200 250 300 350 Tiempo (s) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Coeficientes aerodinámicos UAV CL CD Figura 7.11 Evolución del CLyCDdel UAV.. 7.8 Acelerómetro En este apartado mostramos las medidas tomadas por el acelerómetro. Figura 7.12 Medidas tomadas por el acelerómetro.. Observamos como las mayores variaciones se dan en las componentes longitudinal y vertical, algo totalmente lógico ya que la máquina de control de altitud se encuentra en todo momento controlando la Señal PWM del empuje y la deflexión de los elevadores. El valor negativo de la acelearación vertical se debe a que 64 Capítulo 8. Estudio de Robustez 0 50 100 150 200 250 300 350 Dirección del viento (degrees) 0 5 10 15 20 25 30 m Error cometido en función de la dirección del viento Error máximo Error medio Figura 8.5 Error cometido en función de la dirección del viento.. -200 0 200 400 600 800 1000 Pe(Este)[m] -500 0 500 1000 1500 2000 2500 Pn(Norte)[m] Trayectoria UAV frente a diferentes vientos 0° = Norte 45° = Noreste 90° = Este 135° = Sudeste 180° = Sur 225 ° = Suroeste 270 ° = Oeste 315° = Noroeste WP 8 WP 10 WP 12 WP 13 Figura 8.6 Trayectoria UAV para diferentes vientos.. 8.3 Efecto de la velocidad de vuelo En los apartados anteriores se ha estudiado el efecto del viento en el UAV, comprobándose como la intensidad de este afecta especialmente. Para vientos de intensidad 15 m/s el UAV es capaz de realizar la misión diseñada, pero requiere que el Autopiloto realice deflexiones de las superficies de control de manera muy brusca, algo que no es bueno para la estructura del UAV. Por ello, en este apartado se va estudiar el efecto de incrementar la velocidad de vuelo del UAV para hacer frente a vientos con intensidad considerable. En la tabla 8.3 se muestran los resultados obtenidos para diferentes rachas de viento y velocidad de vuelo del UAV. 8.3 Efecto de la velocidad de vuelo 65 Tabla 8.3 Error cometido frente a diferentes vientos y velocidad de trimado.. Intensidad del Velocidad de Error viento (m/s) trimado (m/s) máximo (m) Error medio (m) Varianza (m2) 5 17 8.6099 3.5131 6.1445 10 17 10.177 5.3421 11.919 15 17 15.276 7.8965 26.54 5 20 9.7684 3.7187 7.8495 10 20 12.739 4.8855 14.793 15 20 19.538 6.8215 38.632 20 20 27.113 7.7207 66.852 10 22 10.064 4.806 11.686 15 22 17.512 5.8693 23.151 20 22 34.96 8.1581 86.942 5 25 10.505 3.7746 10.719 10 25 10.343 4.2825 12.353 15 25 35.865 10.954 124.75 Al aumentar la velocidad de vuelo contrarrestamos los efectos del viento y se reduce considerablemente el tiempo de misión, pero si la ruta definida no es lo suficientemente extensa el error cometido puede aumentar debido a que el UAV no tiene tiempo para poder pasar tan cerca por todos los puntos como cuando se consideraba la misión a velocidades de vuelo más bajas, hasta llegar al extremo de que no es capaz de realizar la misión. Por ejemplo, para un viento de 5 m/s y velocidad de trimado de 25 m/s el UAV se queda estancado en un waypoint describiendo trayectorias circulares alrededor de él. En la Fig. 8.7 se muestra una comparativa de la trayectoria descrita por el UAV para un viento de intensidad 10 m/s y componente WSW para diferentes velocidades de trimado del UAV. -200 0 200 400 600 800 1000 Pe(Este)[m] -500 0 500 1000 1500 2000 2500 Pn(Norte)[m] Trayectoria UAV en función de velocidad de trimado Utrim 15 m/s Utrim 17 m/s Utrim 20 m/s Utrim 22 m/s Utrim 25 m/s Figura 8.7 Trayectoria del UAV en función de la velocidad de trimado.. 9 Implementación con Simulador de Vuelo FlightGear La implementación con el simulador de vuelo FlightGear llevada a cabo durante este Capítulo es un extra opcional, pero que aporta una experiencia visual al objetivo de este Proyecto. FlightGear es un simulador de vuelo multiplataforma y de código libre que es compatible con Matlab-Simulink. A partir de los datos calculados durante la simulación el simulador es capaz de mostrar en tiempo real la misión llevada a cabo por la aeronave. Para poder exportar los datos calculados durante la simulación en Simulink hacemos uso de dos bloque pertenciente al Aeroespace Blockset de Simulink: •FlightGear Preconfigured 6DoF Animation : bloque que permite exportar los valores de posición y actitud de la aeronave al simulador de vuelo de FlightGear. Los valores de posición deben proporcionarse en coordenadas geodésicas, por ejemplo el Sistema WGS84. Para poder configurar el bloque es necesario seleccionar la versión instalada de FlightGear ,introducir la dirección IP del ordenador donde se vaya a realizar la simulación, el puerto de destino (por defecto 5502), y el tiempo de muestro. •Run Generate Script : se encarga de crear el fichero .bat para ser ejecutado por FlightGear. En este bloque es necesario especificar la aeronave con la cual se va realizar la simulación, que en nuestro caso ha sido Malolo1, un ala volante presente en la base de datos de FlightGear con características muy similares. También hay que introducir el aeropuerto en el que se va realizar la simulación (LEZL), la pista de despegue (27), altitud (328 ft) y rumbo (0 o ) iniciales, así como el offset de distancia (NM) y azimut(o). Figura 9.1 Bloques para ejcutar simulación en FlightGear.. A continuación se muestra la configuración introducida en Simulink de ambos bloques. Haciendo click en la ventana ”Generate Script“ creamos el fichero .bat, el cual al abrirlo nos cargará la situación configurada en la Fig. (9.2). 67 68 Capítulo 9. Simulador FlightGear Figura 9.2 Bloque de animación 6DoF de FlightGear.. Figura 9.3 Configuración del Bloque Run Generate Script.. 69 Figura 9.4 Configuración del Bloque Run Generate Script.. Figura 9.5 Configuración del Bloque Run Generate Script.. Así una vez configurados ambos bloques, abrimos el fichero .bat creado y comenzamos la simulación en Simulink. En la Fig. (9.6) se muestra cómo se vería nuestro UAV en FlightGear durante la simulación. 70 Capítulo 9. Simulador FlightGear Figura 9.6 Simulación en FlightGear.. 10 Conclusiones y Desarrollo Futuro El objetivo de este proyecto era la de desarrollar un modelo matemático que describa el comportamiento dinámico del UAV Skywalker X8 y posteriormente estudiar diferentes misiones que pueda llevar a cabo. Una vez concluido podemos hacer una serie de valoraciones acerca de los resultados obtenidos: • Al tratarse de un sistema dinámico no lineal el gasto computacional es elevado. Esto se traduce en la necesidad de emplear un equipo potente y el seleccionar una ruta no demasiado compleja. Particularmente, la simulación se ha realizado en un ordenador provisto con un procesador Intel Core i7-8750H 2.20GHz, memoria RAM de 8GB y tarjeta gráfica NVIDIA GeForce GTX 1060. Para llevar a cabo la simulación en FlightGear fue necesario calcular por una parte los parámetros necesarios para FlightGear y posteriormente introducir en un modelo Simulink distinto dichos parámetros desde el workspace de Matlab para poder realizar la simulación, ya que si se intentaba ejecutar en el mismo archivo la simulación de Simulink y la generación en FligthGear el ordenador proporcioaba resultados incorrectos. • El sistema de misión implementado implica que el UAV pase por todos los waypoints definidos, lo cual puede dar lugar a que en el caso de que dos puntos se encuentren muy juntos tras un viraje el UAV deba salirse de la ruta óptima calculada por el algoritmo TSP. Una posible solución puede ser el diseñar el sistema de misión de manera que realice maniobras circulares, aunque no llegue a pasar exactamente por el punto en custión. • El UAV al tener una masa de tan solo 4kg se ve afectado notablemente por la presencia de viento. A mayor valor de la intensidad del viento más perjudicado se ve el UAV, hasta llegar a un determinado momento en el que no es capaz de hacer frente al viento incidente. Para solucionar este problema es necesario aumentar la velocidad de vuelo para disminuir el efecto de dicho viento, pero tal y como se vio en el Estudio de Robustez en el Capítulo 8, aumentar la velocidad no simpre es posible, ya que puede introducir error en la trayectoria del UAV o que no sea capaz de realizar la misión. Por tanto, se concluye que a la hora de diseñar una misión hay que tener en cuanta una serie de factores muy importantes: previsión meteorológica, extensión de la zona a cubrir, elegir estratégicamente los waypoints de paso y la velocidad de vuelo del UAV. Teniendo en cuenta estas variables el UAV realiza la misión de manera satisfactoria. • La simulación de un sistema tan complejo muestra cómo hay que tener cuidado a la hora de seleccionar todos y cada uno de los parámetros de diseño, ya que pequeños cambios en cualquier parámetro puede dar lugar a que el sistema no se comporte como se espera. En cuanto a posibles líneas de desarrollo se pueden incluir las siguientes: • Implementar un sistema de misión que sea capaz de detectar obstáculos y de realizar rutas tridimensionales. •Sustituir los controladores PID implementados por leyes de control óptimo de trayectorias. •Introducir diferentes modelos de perturbaciones atmósfericas para comprobar resultados. •Comparar los resultados obtenidos en la simulación frente a datos de telemetría reales. 71 Anexos 73 Índice de Figuras 1.1 X8 Skywalker. 2 1.2 Sistema de coordenadas de Ejes Tierra. 2 1.3 Sistema de Ejes de Navegación. 3 1.4 Sistema de Ejes Cuerpo. 3 1.5 Sistema de Ejes Viento. 3 2.1 Variables de estado. 6 2.2 Esferoide del sistema WGS84. 8 2.3 Fuerzas y momento aerodinámicos aplicadas en el centro aerodinámico del perfil. 8 2.4 Variables de control aerodinámicas en configuración estándar. 9 2.5 Variables de control aerodinámicas para cola en V. 10 2.6 Variables de control aerodinámicas del ala volante. 10 2.7 CLfrente a α. 12 2.8 CDfrente a α. 12 3.1 Regímenes de vuelo para el autopiloto longitudinal. 19 3.2 Diagrama de bloques del controlador de cabeceo. 20 3.3 Diagrama de bloques del controlador de altitud. 21 3.4 Diagrama de bloques del controlador de velocidad mediante cabeceo. 22 3.5 Diagrama de bloques del controlador de velocidad mediante empuje. 23 3.6 Máquina de estados para control longitudinal. 24 3.7 Diagrama de bloques del controlador de alabeo. 25 3.8 Diagrama de bloques del controlador de rumbo. 26 4.1 Altitud deseada para diseño de ruta longitudinal. 28 4.2 Algoritmo empleado para el cálculo de la ruta en línea recta. 29 4.3 Criterio de actualización de segmento de ruta.. 29 4.4 Algoritmo de gestión de ruta. 30 4.5 Ejemplo del Problema TSP. 31 4.6 Algoritmo TSP. 32 5.1 Parámetros del modelo de error de Gauss-Markov. 35 6.1 Bloque Cinemática lineal en Simulink. 39 6.2 Bloque Cinemática angular en Simulink. 39 6.3 Bloque Dinámica lineal en Simulink. 40 6.4 Bloque Dinámica angular en Simulink. 40 6.5 Bloque WGS84 de Simulink. 41 6.6 Bloque Coeficientes Aerodinámicos en Simulink. 41 6.7 Bloque de Modelado del Viento en Simulink. 41 6.8 Bloque de Dryden de Simulink. 42 6.9 Viento predominante a 10 m de altitud. 42 81 82 Capítulo A. Índice de Figuras 6.10 Viento predominante a 80 m de altitud. 42 6.11 Bloque de Fuerzas en Simulink. 43 6.12 Bloque de Momentos en Simulink. 43 6.13 Bloque de controladores en Simulink. 44 6.14 Máquina de estados de control longitudinal en Simulink. 44 6.15 Controlador de Cabeceo en Simulink. 44 6.16 Controlador de Altitud en Simulink. 45 6.17 Controlador de Velocidad mediante Cabeceo en Simulink. 45 6.18 Controlador de Velocidad mediante Empuje en Simulink. 45 6.19 Controlador de Alabeo en Simulink. 46 6.20 Controlador de Rumbo en Simulink. 46 6.21 Sistea de misión en Simulink. 46 6.22 Acelerómetro en Simulink. 47 6.23 Giróscopo en Simulink. 47 6.24 Altímetro barométrico en Simulink. 48 6.25 Anemómetro en Simulink. 48 6.26 Magnetómetro en Simulink. 48 6.27 Sistema GPS en Simulink. 49 6.28 EKF de Posición en Simulink. 49 6.29 EKF de Actitud en Simulink. 50 7.1 Misión realizada por el UAV en Matlab-Simulink. 52 7.2 Posición Norte del UAV. 52 7.3 Posición Este del UAV. 53 7.4 Altitud del UAV. 53 7.5 Alabeo del UAV. 54 7.6 Cabeceo del UAV. 54 7.7 Guiñada del UAV. 55 7.8 Velocidad del UAV. 55 7.9 Señal PWM encargada de modelar el empuje. 56 7.10 Deflexión de los elevones del UAV. 56 7.11 Evolución del CLyCDdel UAV. 57 7.12 Medidas tomadas por el acelerómetro. 57 7.13 Medidas tomadas por el giróscopo. 58 7.14 Comparación entre velocidad real y velocidad dada por GPS. 59 7.15 Rumbo del UAV dado por GPS. 59 8.1 Viento registrado en el Aeropuerto de Sevilla en 2018. 61 8.2 Error cometido en función de la intensidad del viento. 62 8.3 Trayectoria UAV para diferentes vientos. 62 8.4 Trayectoria UAV para diferentes vientos. 63 8.5 Error cometido en función de la dirección del viento. 64 8.6 Trayectoria UAV para diferentes vientos. 64 8.7 Trayectoria del UAV en función de la velocidad de trimado. 65 9.1 Bloques para ejcutar simulación en FlightGear. 67 9.2 Bloque de animación 6DoF de FlightGear. 68 9.3 Configuración del Bloque Run Generate Script. 68 9.4 Configuración del Bloque Run Generate Script. 69 9.5 Configuración del Bloque Run Generate Script. 69 9.6 Simulación en FlightGear. 70 A.1 Vuelo a 15 m/s con viento de intensidad 10 m/s y dirección WSW. 76 A.2 Trayectoria del UAV en función de la intensidad del viento. 77 A.3 Trayectoria del UAV en función de la dirección del viento. 78 A.4 Trayectoria del UAV en función de la dirección del viento. 79 Índice de Tablas 2.1 Variables de estado. 5 2.2 Parámetros UAV Skywalker X8. 14 2.3 Parámetros de ráfaga del modelo Dryden [6] 16 8.1 Error cometido frente a diferentes vientos. 62 8.2 Error cometido frente a diferentes vientos. 63 8.3 Error cometido frente a diferentes vientos y velocidad de trimado. 65 83 Bibliografía [1] Gryte, K. High Angle of Attack Landing of an Unmanned Aerial Vehicle.MSc Tesis. Norwegian University of Science and Technology. Trondheim, Norway. 2015. [2] Randal W. Beard, Timothy W. McLain. Small Unmanned Aircraft: Theory and Practice. Princeton University Press. Princeton and Oxford. [3] Rísquez Ruiz, A. Modelado y simulación del Skywalker X8. 2017. [4] R. F. Stengel. Flight Dynamics. Princeton, NJ: Princeton University Press. 2004. [5] T. R. Yechout, S. L. Morris, D. E. Bossert,W. F. Hallgren. Introduction to Aircraft Flight Mechanics. AIAA Education Series, American Institute of Aeronautics and Astronautics. 2003. [6] J. W. Langelaan, N. Alley, and J. Niedhoefer. Wind field estimationfor small unmanned aerial vehicles. AIAA Guidance, Navigation, and Control Conference, Toronto, Canada. 2010. [7] MathWorks. Traveling Salesman Problem: Solver-Based.https:// es.mathworks.com/ help/ optim/ ug/ travelling-salesman-problem.html?lang=en [8] Fernández Camacho, E. Apuntes de la asignatura “Sistemas de control y guiado”. Universidad de Sevilla. Sevilla, España. 2016. [9] Limón Marruedo, D. Apuntes de la asignatura “Sistemas de control y guiado”. Universidad de Sevilla. Sevilla, España. 2016. [10] Gavilán Jiménez, F. Apuntes de la asignatura “Mecánica del Vuelo Avanzada”. Universidad de Sevilla. Sevilla, España. 2017. [11] Bolzern, Paolo ; Scattolini, Riccardo, coaut. Fundamentos de Control Automático. 2009. [12] Página web de FlightGear. https:// www.flightgear.org/ [13] Enlace de descarga del modelo de aeronave para FlightGear. http:// ftp.igh.cnrs.fr/pub/ flightgear/ ftp/ Aircraft-2.0.0/ Malolo1_0.0.zip [14] Enlace de descarga del escenario utilizado en FlightGear. http:// ns334561.ip-5-196-65.eu/ ~fgscenery/ WS2.0/ scenery-2.0.1.html [15] MathWorks. FlightGear Preconfigured 6DoF Animation.https:// es.mathworks.com/ help/ aeroblks/ flightgearpreconfigured6dofanimation.html 85