Planificación de trayectorias de avión en espacios aéreos estructurados con incertidumbre meteorológica

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 Trabajo Fin de Máster Ingeniería Aeronáutica Planificación de Trayectorias de Avión en Espacios Aéreos Estructurados con Incertidumbre Meteorológica Autor: Manuel Gutiérrez Fernández Tutor: Antonio Franco Espín Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Máster Ingeniería Aeronáutica Planificación de Trayectorias de Avión en Espacios Aéreos Estructurados con Incertidumbre Meteorológica Autor: Manuel Gutiérrez Fernández Tutor: Antonio Franco Espín Profesor Contratado Doctor Dpto. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Trabajo Fin de Máster: Planificación de Trayectorias de Avión en Espacios Aéreos Estructurados con Incertidumbre Meteorológica Autor: Manuel Gutiérrez Fernández Tutor: Antonio Franco Espín 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: A mi familia y amigos A todos los que me han hecho posible llegar hasta aquí Agradecimientos El fin de otra nueva etapa se acerca y con él la realización de este trabajo. Ahora, toca agradecer a todos los profesores el esfuerzo que han hecho para formar no sólo a mí, sino a todo un repertorio de grandes profesionales de la ingeniería. La consecución de este trabajo demuestra de cara al público lo que los ingenieros pueden hacer. Agradecer a la Escuela Superior de Ingeniería de Sevilla el permitirme estudiar en su centro proporcionando todos los medios necesarios para llevar a cabo mi aprendizaje durante estos años. En especial, al Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos por contar con unos profesionales estupendos que contribuyen en gran medida a que esta Escuela tenga el prestigio que se merece. Mencionar a mis padres que mediante su contribución moral y económica, han posibilitado en gran medida que esté aquí en este momento a punto de ser ingeniero superior, escribiendo el trabajo fin de máster. De manera especial, dedicar unas palabras a Antonio Franco Espín, tutor del trabajo, por su implicación constante, su esfuerzo, paciencia, confianza y capacidad para transmitirme los conocimientos desarrollados en el trabajo. También, agradecer a Alfonso Valenzuela Romero por conseguir junto con Antonio dar visibilidad y financiación al proyecto. Por último, no puedo olvidarme de mencionar al Centro de Referencia de Investigación, Desarrollo e Innovación ATM (CRIDA) por la información aportada al trabajo y por el apoyo económico al proyecto. Concretamente, destacar a Nicolás Suárez Tetzlaff y a Christian Eduardo Verdonk Gallego de CRIDA por su seguimiento y realimentación durante estos últimos 6 meses de trabajo. Manuel Gutiérrez Fernández Autor del trabajo Sevilla, 2020 III Índice Resumen V Abstract VII Índice Abreviado IX Notación XIII Índice de Figuras 1 Índice de Tablas 3 1 Introducción 5 2 Estructura del trabajo 7 3 Espacio Aéreo Estructurado 9 3.1 Características del Espacio Aéreo 9 3.1.1 Reducción del Espacio Aéreo 10 3.1.2 Eliminación de Nodos Inconexos 11 3.1.3 Zonas de Espacio Aéreo D-P-R 12 4 Predicción de la Trayectoria 15 4.1 Ecuaciones del movimiento del avión 15 4.2 Introducción de la Incertidumbre Meteorológica 17 4.3 Vuelo de un segmento a curso constante 17 4.4 Cálculo del tiempo de vuelo 18 4.4.1 Cálculo de la Velocidad con respecto de Tierra 18 4.4.2 Método de Heun 21 4.5 Cálculo del pseudotiempo de vuelo 23 5 Planificación de la Trayectoria 25 5.1 Minimización del tiempo/pseudotiempo de vuelo promedio 25 5.2 Minimización del tiempo de vuelo del peor escenario 27 5.2.1 Método heurístico para la resolución del problema min-max. 28 5.3 Minimización del combustible promedio final 29 5.4 Variaciones del coste con el horizonte temporal 32 5.4.1 Introducción explícita del tiempo de vuelo 32 5.4.2 Método iterativo para incluir las variaciones temporales 36 5.5 Limitación del ángulo de viraje 38 XI XII Índice 6 Aplicación práctica. Resultados 41 6.1 Comparación entre distintos criterios de optimización 41 6.2 Minimización del combustible promedio final 43 6.3 Método Heurístico vs Enfoque Combinatorio 45 6.4 Influencia de las variaciones temporales del coste 47 6.5 Efecto de la consideración de regiones D-P-R activas 49 6.6 Influencia de establecer un límite superior al ángulo de viraje 51 6.6.1 Comparación entre el método Dijkstra modificado y el enfoque combinatorio 53 6.7 Estudio de la estabilidad temporal de la solución 54 6.8 Impacto del criterio de optimización en un escenario de tráfico 61 6.9 Influencia de la altitud de vuelo 63 7 Conclusiones. Futuras Mejoras 69 7.1 Conclusiones Generales 69 7.2 Futuras Mejoras 70 Referencias 71 Notación |X|Valor absoluto de X. [X]Vector X. ¯ XValor promedio de la variable X. λLongitud de un punto de la superficie terrestre. φLatitud de un punto de la superficie terrestre. γaCoeficiente de dilatación adiabática del aire. RaConstante de los gases ideales para el aire. χgCurso entre dos puntos. χRumbo entre dos puntos. MNúmero de Mach. AT M Air Traffic Management. RNAV Navegación de Área. SID Standard Instrumental Departure. STAR Standard Terminal Arrival Route. ISA International Standard Atmosphere. EPS Ensemble Prediction System. ECMWF European Centre for Medium-Range Weather Forecast. ILP Integer Linear Programming. XIII XIV Notación INLP Integer Non-Linear Programming. MILP Mixed-Integer Linear Programming. MIQP Mixed-Integer Quadratic Programming. MINLP Mixed-Integer Non-Linear Programming. PSD Positive Semi-Definite Matrix. Índice de Figuras 3.1 Nodos de un espacio aéreo estructurado. Los asteriscos rojos representan los puntos de origen y destino de la trayectoria 10 3.2 Filtrado del espacio aéreo (K=1.1). Los asteriscos rojos representan los puntos de origen y destino de la trayectoria 11 3.3 Filtrado tras eliminar los nodos inconexos 12 3.4 Regiones D-P-R en España. Fuente: Insignia (Enaire) 12 3.5 Región D-P-R. Nodos y conexiones a eliminar (color rojo) 13 4.1 Superficie isobara y de altitud geométrica constante 16 4.2 Interpolación de los vientos en la malla 19 4.3 Cálculo del viento longitudinal y cruzado 20 4.4 Triángulo de vientos 21 4.5 Ejemplo de un segmento de trayectoria 21 4.6 Método de Heun. Fuente: ETSII ([9]) 22 5.1 Costes vs tre f 32 5.2 Grafo de ejemplo 38 6.1 Caminos críticos en función del criterio elegido 42 6.2 Masa Consumida vs Pseudotiempo 44 6.3 Resultados del Método Heurístico y de la Optimización Combinatoria 46 6.4 Caminos óptimos sin dependencia temporal de los costes 47 6.5 Camino óptimo sin dependencia temporal de los costes: tDep = +22 h 48 6.6 Camino óptimo con dependencia temporal de los costes: tDep = +22 h 48 6.7 Espacio Aéreo Estructurado. Regiones D-P-R 49 6.8 Caminos óptimos con regiones D-P-R I 50 6.9 Caminos óptimos con regiones D-P-R II 51 6.10 Limitación del ángulo de viraje 53 6.11 Limitación del ángulo de viraje: Dijkstra vs Combinatoria 54 6.12 Caminos críticos para el 04/01/2019 a las 19:00 55 6.13 Caminos críticos para el 08/01/2019 a las 19:00 56 6.14 Caminos críticos para el 12/01/2019 a las 19:00 57 6.15 Caminos críticos para el 16/01/2019 a las 19:00 57 6.16 Caminos críticos para el 20/01/2019 a las 19:00 58 6.17 Caminos críticos para el 24/01/2019 a las 19:00 59 6.18 Caminos críticos para el 28/01/2019 a las 19:00 60 6.19 Escenarios de tráfico con 1000 vuelos aleatorios 62 1 2Índice de Figuras 6.20 Variación de las constantes ˆ Ayˆ Bcon la altitud 64 6.21 Consumo de combustible vs altitud para ¯ τf=15500 s65 6.22 Caminos críticos en función de la altitud 67 Índice de Tablas 6.1 Parámetros de la aeronave B763 41 6.2 Tiempos y Pseudotiempos en función del criterio elegido 42 6.3 Masa de combustible y tiempo de ejecución del problema cuadrático 43 6.4 Pseudotiempo medio, masa de combustible y tiempo de ejecución del problema lineal 44 6.5 Comparación Heurísitca y Combinatoria 46 6.6 Tiempos de vuelo para tDep = +22 h 49 6.7 Comparación de los tiempos de vuelo al incluir las variaciones temporales del coste 51 6.8 Comparación Dijkstra y Combinatoria ∆χlim =30◦54 6.9 Comparación tiempos de vuelo para el 04/01/2019 a las 19:00 55 6.10 Comparación tiempos de vuelo para el 08/01/2019 a las 19:00 55 6.11 Comparación tiempos de vuelo para el 12/01/2019 a las 19:00 56 6.12 Comparación tiempos de vuelo para el 16/01/2019 a las 19:00 57 6.13 Comparación tiempos de vuelo para el 20/01/2019 a las 19:00 58 6.14 Comparación tiempos de vuelo para el 24/01/2019 a las 19:00 58 6.15 Comparación tiempos de vuelo para el 28/01/2019 a las 19:00 59 6.16 Tiempos de vuelo para 1000 vuelos aleatorios 61 6.17 Tiempos, Pseudotiempos y Consumo de Combustible en función de la altitud 63 3 1 Introducción El sistema ATM (Air Traffic Management) está atravesando un nuevo cambio dirigido a reducir el impacto ambiental de las emisiones de gases contaminantes, principalmente CO 2 y NO x . Para ello, una correcta planificación en la trayectoria de la aeronave que permita minimizar el consumo de combustible, es uno de los pilares fundamentales con los que debe contar el sistema ATM para cumplir con sus objetivos propuestos. Por este motivo, este trabajo se centra en estudiar los problemas que minimizan el consumo de combustible promedio de forma directa, sin utilizar el tiempo de vuelo como variable intermedia. La razón para resolver el problema directamente es que la relación no lineal entre ambas variables puede provocar que el camino óptimo para minimizar el tiempo de vuelo no coincida con el camino obtenido de minimizar el consumo de combustible. Por ello, establecer las diferencias entre ambas funciones objetivos será uno de los asuntos importantes a tratar en este trabajo. Por otro lado, la incorporación de la incertidumbre meteorológica y su posterior estudio desde el punto de vista de sus efectos en la trayectoria es crucial para elegir trayectorias con un menor impacto ambiental. Una tendencia actual en el modelado de la incertidumbre meteorológica es el uso de la predicción meteorológica por conjuntos o EPS (Ensemble Prediction System). Con estos conjuntos, las variables meteorológicas que intervienen en la optimización de la trayectoria no vienen dadas por un único valor, ni siquiera por una distribución de probabilidad que refleja su comportamiento, sino por un conjunto de valores discretos obtenidos mediante una ligera alteración de las condiciones iniciales y/o los parámetros del modelo, o incluso modificando el horizonte temporal de la predicción. De esta forma, se contemplan una serie de evoluciones atmosféricas posibles en función de una combinación de valores entre la meteorología y el resto de parámetros del problema. Además, el hecho de que la meteorología dependa del instante temporal, provocará cambios en las funciones objetivo consideradas. Las implicaciones cualitativas y cuantitativas que dichas variaciones temporales tienen en el camino óptimo de la trayectoria, serán abordadas a lo largo de este trabajo. En el problema de optimización de trayectorias planteado, sólo se estudiará el perfil de vuelo horizontal de una aeronave durante la fase de crucero, al considerar que el consumo en esta fase es mucho más importante que en otras fases como pueden ser las subidas tras despegue o los descensos para aterrizaje. Además las salidas y llegadas suelen estar normalizadas mediante las cartas aeronáuticas SID y STAR, por lo que el margen de elección de la trayectoria en dichas fases está muy limitado por estos procedimientos. Los métodos empleados para la resolución de los problemas de optimización de trayectorias, se basan en la Teoría de Grafos; en ellos, la aeronave recorre una ruta dentro de un espacio aéreo 5 12 Capítulo 3. Espacio Aéreo Estructurado 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N Figura 3.3 Filtrado tras eliminar los nodos inconexos. Aunque no se aprecie en la figura anterior, el número de nodos se ha reducido a 4817, mientras que el de conexiones alcanza un valor de 11868. Esto supone una reducción de aproximadamente un 18.5% en el número de nodos y de un 20% en el número de conexiones. Dicha reducción, se traducirá en una potencial disminución del tiempo de ejecución de los posteriores problemas de optimización. 3.1.3 Zonas de Espacio Aéreo D-P-R También es importante tener en cuenta que existen regiones del espacio aéreo en las cuales la aeronave tiene prohibido o limitado el vuelo. Estas regiones son denominadas como peligrosas, prohibidas o restringidas (D-P-R). Figura 3.4 Regiones D-P-R en España. Fuente: Insignia (Enaire). 3.1 Características del Espacio Aéreo 13 A la hora de realizar la planificación, hay que asegurarse de que ninguna de las conexiones por las que pasa la trayectoria de la aeronave atraviesan las regiones anteriormente citadas. Para ello, es conveniente realizar un filtrado del grafo que elimine tanto los nodos pertenecientes al interior de dichas regiones, como las conexiones entre nodos (los cuales están fuera de la región D-P-R), en la que alguna de ellas atraviesa esas porciones del espacio aéreo. Figura 3.5 Región D-P-R. Nodos y conexiones a eliminar (color rojo). Al realizar el filtrado se han eliminado todas o parte de las conexiones entre algunos de los nodos, sin que estos tengan por qué ser descartados. Esta acción puede provocar que algunos de esos waypoints queden no conectados con el origen. Por este motivo, es adecuado volver a eliminar los nodos inconexos mediante el proceso descrito en el apartado anterior. Para terminar, cabe decir que muchas de las regiones D-P-R, no tienen por qué estar activas permanentemente, ni siquiera durante todo un día completo. Sin embargo, ese cambio dinámico de la estructura del espacio aéreo a lo largo del tiempo conlleva una mayor complejidad y no será tenido en cuenta en este trabajo. Por tanto, se limitará el espacio aéreo a un instante temporal concreto. 4 Predicción de la Trayectoria En este capítulo se desarrolla toda la teoría necesaria para comprender cómo se puede llegar a calcular el tiempo de vuelo y el consumo de combustible de una trayectoria sujeta a incertidumbre meteorológica. En primer lugar, se establecen las ecuaciones que definen el movimiento del avión. Posteriormente, se explica el modo en el que se introduce la incertidumbre en el problema. Por último, a partir de los resultados anteriores se plantean las funciones que calculan el tiempo de vuelo y el consumo de combustible. 4.1 Ecuaciones del movimiento del avión Tal y como se comentó en la introducción, la trayectoria que se va a estudiar en este trabajo desacopla la cinemática y la dinámica del avión en el plano vertical y horizontal; sólo teniendo en cuenta las ecuaciones en el plano horizontal. Por otra parte, la única fase que se considera es la fase de crucero. La explicación de esto es que la aeronave estará volando la mayor parte del tiempo en esta fase. Aunque el resto de fases como subidas tras despegue, descensos para aterrizaje o virajes pueden constituir un porcentaje importante de la trayectoria, el margen de maniobra a la hora de planificar dichos segmentos es muy escaso, ya que la definición de esos segmentos está sujeta a procedimientos normalizados a través de las cartas aeronáuticas (SID y STAR); o bien por las servidumbres aeronáuticas u otra restricción similar. Las ecuaciones del movimiento de un crucero para un segmento de vuelo json las siguientes: dr dt =Vg dm dt =−cTFT FT=D L =mg (4.1) donde FT , D y L son el empuje, la resistencia aerodinámica y la sustentación; m es la masa de la aeronave; g , es la aceleración de la gravedad y cT es el consumo específico de combustible. En principio, se asumirá que dicho consumo específico es constante, y su valor será el correspondiente a la velocidad y a la temperatura de referencia (proporcionada por la ISA a la altitud de crucero). El vuelo de crucero se va a realizar asumiendo como ligaduras de vuelo que el número de Mach y la presión son ctes (M=cte;p=cte) . Cabe comentar que, en presencia de una atmósfera real, un crucero a p=cte no implica necesariamente un vuelo a h=cte , ya que volar a través de una superficie isobara implica cambios en la altitud geométrica. Sin embargo, a efectos prácticos se considerará como la presión proporcionada por la atmósfera ISA a una altitud determinada. 15 16 Capítulo 4. Predicción de la Trayectoria Figura 4.1 Superficie isobara y de altitud geométrica constante. La sustentación y la resistencia pueden ser escritas en función de las ligaduras de vuelo como L=1/2γapM2SCL y D=1/2γapM2SCD , donde γa=1.4 es la relación entre los calores específicos del aire a presión constante y volumen constante, respectivamente; S , es la superficie del ala del avión en cuestión; y CD y CL son los coeficientes de resistencia y de sustentación del avión respectivamente. Por otro lado, el modelo de avión considerado en este trabajo viene definido a través de la base de datos BADA de Eurocontrol. El coeficiente de resistencia puede relacionarse con el coeficiente de sustentación a través de una polar parabólica de coeficientes constantes CD=CD0+CD2C2 L , donde CD0es el coeficiente de resistencia parásita y CD2es el coeficiente de resistencia inducida. Si se sustituyen las expresiones de la resistencia y de la sustentación en la ecuación de la variación de masa se obtiene la siguiente expresión. dm dt =−(A+Bm2)(4.2) Donde A y B son dos constantes que dependen de los parámetros aerodinámicos y propulsivos del problema, y vienen dadas por las expresiones siguientes. A=1 2γapM2SCD0cT B=2cTCD2g2 γapM2S (4.3) Nótese que A y B son constantes porque se está asumiendo que el consumo específico, cT , también lo es. Dividiendo la derivada de la distancia con respecto al tiempo entre la variación de la masa con respecto al tiempo e integrando a lo largo del segmento j, se tiene la siguiente relación. Z(mf)j (mi)j dm A+Bm2=−Z(rf)j (ri)j dr Vg(r)(4.4) En este trabajo se va a considerar que la masa inicial es conocida (mi)1=mi , y que se cumple continuidad de masas para el resto de segmentos (mi)j= (mf)j−1 , 2≤j≤q , siendo q , el número de segmentos de los que consta la trayectoria. Para cada segmento, la integral de la derecha de la ecuación anterior representa el intervalo de tiempo, tj que la aeronave tarda en recorrer dicha conexión. Por tanto, si se resuelve la ecuación para cada segmento el resultado es una expresión explícita para la masa final. arctan"rB A(mf)j#=arctan"rB A(mi)j#−√ABtj(4.5) 4.2 Introducción de la Incertidumbre Meteorológica 17 Sumando todas las conexiones, se puede obtener el tiempo de vuelo total. tf= q ∑ j=1 tj(4.6) Por otro lado, la masa final viene dada por mf= (mf)q . Entonces, la cantidad de fuel consumida durante todo el crucero, mF=mi−mf, tiene el siguiente valor. mF=mi−rA Btan(arctan rB Ami!−√ABtf)(4.7) 4.2 Introducción de la Incertidumbre Meteorológica La novedad que introduce este trabajo es la incorporación de la incertidumbre meteorológica en el problema de planificación de trayectorias, con el objetivo de que la trayectoria que finalmente se decida volar tenga en cuenta las posibles variaciones en tiempos de vuelo y en consumo de combustible debidas a la meteorología presente durante el vuelo. Las variables que se van a introducir son: el campo de temperaturas, T y los valores del viento zonal (viento oeste-este) y del viento meridional (viento sur-norte), wu y wv , respectivamente. Es habitual no considerar valores únicos para las variables meteorológicas sino utilizar predicciones por conjuntos denominadas EPS (Ensemble Prediction System). Cada uno de los miembros del conjunto es denominado ensemble member y es proporcionado por el European Centre for Medium-Range Weather Forecast (ECMWF). El ECMWF proporciona un grupo de predicciones obtenidas mediante ligeras variaciones de las condiciones iniciales de uno o varios modelos. El resultado es un conjunto de predicciones futuras que permiten cuantificar la incertidumbre existente en la naturaleza, ya sea por los posibles errores en la medición o por el desconocimiento de gran parte de la dinámica terrestre, [5] [6]. Cada una de esas predicciones será también denominada con el nombre de escenario con el objetivo de enfatizar que se trata de una predicción entre múltiples predicciones abordadas. A su vez, los valores de vientos y temperaturas serán proporcionados para una determinada malla en latitud y longitud delimitada por las coordenadas φmin,λmin,φmax,λmax ; y con un determinado paso, ∆φ y ∆λ . Las predicciones se cierran con los datos disponibles en un día y en un instante determinado T0 , pero pueden hacer referencia a un horizonte temporal diferente ( +6 horas ,+12 horas,+24 horas ...), lo que se denomina time-lagged. Téngase en cuenta que el ECMWF no es capaz de proporcionarnos la predicción hasta 7 horas después de haberla cerrado, por lo que si se desea conocer la predicción meteorológica de la trayectoria con un día de antelación, habrá que utilizar predicciones con un horizonte temporal de 31 horas. Por lo tanto, los vientos y la temperatura , wu(φ,λ,n,t) , wv(φ,λ,n,t) y T(φ,λ,n,t) , serán función de 4 parámetros: los dos primeros hacen referencia a los valores de latitud y longitud del punto de la malla en cuestión; el tercero indica el escenario que se utiliza y el último parámetro define el instante temporal de la predicción que se emplea en el cálculo. 4.3 Vuelo de un segmento a curso constante La ruta a seguir por la aeronave estará compuesta por una serie de segmentos, los cuales serán recorridos a curso constante. Por este motivo, para llegar al nodo final de cada una de las conexiones habrá que recorrer una distancia en línea recta rj , con un curso χgj . Los valores de la distancia 18 Capítulo 4. Predicción de la Trayectoria recorrida rj y del curso χgj pueden ser calculados empleando las expresiones de la ruta loxodrómica. χgj=arctan    λF−λI lntan(π/4−φI/2) tan(π/4−φF/2)    (4.8) rj       (RE+h)(φF−φI) cos χgj si φF6=φI (RE+h)cosφI|λF−λI|si φF=φI (4.9) con REel radio de Tierra, que se toma constante e igual a 6371 km. El conocimiento de estos parámetros será fundamental para el cálculo posterior de la velocidad con respecto de tierra, Vg, a lo largo de cada uno de los segmentos que contiene la trayectoria. 4.4 Cálculo del tiempo de vuelo En secciones anteriores se ha desarrollado la expresión que permite calcular la masa de combustible consumida durante la trayectoria de crucero de una aeronave. Para obtener dicho valor de consumo es necesario conocer el tiempo de vuelo, tf , además de la masa inicial y otra serie de parámetros aerodinámicos y propulsivos. El tiempo de vuelo es el resultado de calcular la integral derecha de la ecuación (4.4). Z(rf)j (ri)j dr Vg(r)(4.10) Sin embargo, la definición de la integral anterior sólo es válida si el problema es autónomo, es decir, si la velocidad con respecto de tierra depende sólo de la distancia recorrida Vg(r) y no del tiempo de vuelo. La dependencia con el tiempo puede ser posible puesto que la velocidad con respecto de tierra dependerá de la meteorología presente en cada punto de integración, la cual a su vez viene condicionada por el instante temporal en el cual han sido calculados los vientos y temperaturas. Por lo tanto, es conveniente definir el cálculo del tiempo de vuelo para un caso genérico, el cual vendrá dado por la resolución de la siguiente ecuación diferencial. dt dr =1 Vg(r,t)=f(r,t)(4.11) A lo largo de esta sección nos centraremos tanto en el cálculo del integrando Vg(r,t) , a partir de los datos que son conocidos como en la descripción del método elegido para la resolución de la ecuación diferencial. 4.4.1 Cálculo de la Velocidad con respecto de Tierra En este apartado se explicará todo el desarrollo que es necesario seguir para calcular la velocidad con respecto de tierra, Vg(r,t) . El valor de dicha velocidad servirá para obtener el tiempo de vuelo tf, que tarda la aeronave en recorrer el crucero propuesto. 4.4 Cálculo del tiempo de vuelo 19 Para calcular el tiempo de vuelo de cada conexión habrá que dar un número total de pasos Npj , hasta recorrer el segmento completo. Para cada paso de integración, npj , la distancia recorrida según el curso del segmento χgj se llamará rnp,j . Por otro lado, los únicos datos que se conocen a priori son los puntos inicial y final (I,F) de cada conexión y la hora asignada para el despegue tDep . De acuerdo a ello, para calcular la Vg(r,t)habrá que seguir una serie de pasos: 1. Conociendo el punto inicial del segmento (φI,λI) y el curso χgj , se calcula el punto (φint ,λint) al que se llegaría siguiendo una loxodrómica hasta una distancia rnp,j. 2. Con las coordenadas del punto obtenido se realiza una interpolación bilineal de la meteorología (vientos y temperaturas) en el espacio. Para ello, se definen una serie de coeficientes a partir de las coordenadas del punto (φint ,λint) y de sus puntos adyacentes (anterior y posterior), de la malla de predicciones de las que se disponga. a1=(λj+1−λint )(φi+1−φint) (λj+1−λj)(φi+1−φi)(4.12) a2=(λj+1−λint )(φint −φi) (λj+1−λj)(φi+1−φi)(4.13) a3=(λint −λj)(φi+1−φint) (λj+1−λj)(φi+1−φi)(4.14) a4=(λint −λj)(φint −φi) (λj+1−λj)(φi+1−φi)(4.15) Figura 4.2 Interpolación de los vientos en la malla. La temperatura y los vientos interpolados espacialmente se obtienen de multiplicar las constantes por los valores de las esquinas de la malla. wu,esp =a1wui,j+a2wui+1,j+a3wui,j+1+a4wui+1,j+1(4.16) wv,esp =a1wvi,j+a2wvi+1,j+a3wvi,j+1+a4wvi+1,j+1(4.17) Tesp =a1Ti,j+a2Ti+1,j+a3Ti,j+1+a4Ti+1,j+1(4.18) 20 Capítulo 4. Predicción de la Trayectoria 3. A partir de la hora de despeque, tDep , se calculan los vientos y temperaturas que actúan en el instante actual sobre la aeronave. Para obtenerlos, se realiza una interpolación lineal entre las predicciones de vientos y temperaturas interpoladas espacialmente para los instantes anteriores y posteriores al instante inicial de la trayectoria. tint =tDep b1=tk+1−tint tk+1−tk b2=tint −tk tk+1−tk wu,int =b1wu,esp(tk)+b2wu,esp(tk+1) wv,int =b1wv,esp(tk)+b2wv,esp(tk+1) Tint =b1Tesp(tk)+b2Tesp(tk+1) (4.19) 4. Se descompone el viento interpolado en dos componentes: una en la dirección de la trayectoria, wl , y otra de viento cruzado, wx . Esta descomposición es bastante simple, ya que es solamente un giro de ángulo igual al curso del segmento que está volando la aeronave. La transformación queda como: Figura 4.3 Cálculo del viento longitudinal y cruzado. wl=wv,int cos χgj+wu,int sin χgj wx=wv,int sin χgj−wu,int cos χgj (4.20) Por tanto, para poder volar con el curso χgj , orientamos el rumbo de nuestra aeronave, χ , de manera que la composición del vector velocidad aerodinámica, V , y del viento w , coincida en módulo, dirección y sentido con la velocidad con respecto de tierra, Vg , la cual lleva la dirección del curso a seguir. Dicho esto, la velocidad con respecto de tierra puede calcularse mediante el siguiente triángulo. 4.4 Cálculo del tiempo de vuelo 21 Figura 4.4 Triángulo de vientos. Vg=qV2−w2 x+wl=q(γRaT)M2−w2 x+wl(4.21) 4.4.2 Método de Heun Existen numerosos métodos para la resolución de ecuaciones diferenciales. Son muy típicos los métodos que utilizan un paso de integración variable, ya que pueden ser más eficientes de cara a la computación cuando la función que constituye el integrando es continua y derivable. En nuestro caso, la velocidad con respecto de tierra, Vg(r,t) , depende de los vientos y temperaturas, los cuales están definidos para una malla con unas dimensiones concretas. Como es necesario conocer la velocidad en puntos intermedios del segmento, habrá que interpolar la meteorología para obtener su valor en ese punto. Además, cualquier segmento j que pertenezca a la trayectoria, puede estar contenido entre dos valores de latitud y longitud cuya diferencia sea mayor que el paso de la malla |φF−φI|≥∆φ , |λF−λI|≥∆λ . Estos dos hechos provocan que la función f(r,t) sea continua a lo largo del segmento (porque así lo es la interpolación), pero que no sea derivable en los puntos de la malla donde se proporcionan los valores de los vientos y de las temperaturas. Por lo tanto, usar métodos de paso variable puede hacer que sean más lentos a causa de las zonas donde no es derivable la función, ya que la falta de derivabilidad se intenta cubrir con una acumulación de puntos de integración en dicha zona. Como solución a este problema se utilizarán los métodos de paso fijo. La distancia del paso, ∆r , será un parámetro a seleccionar en el problema. Una vez elegida dicha distancia, el número total de pasos de integración para cada conexión vendrá dado por Npj=drj ∆re . El resultado del cociente será redondeado al entero superior, recalculándose la nueva distancia de paso para cada conexión como ∆r=rj Npj . Por tanto, la distancia recorrida hasta el paso de integración, npj será rnp,j=∆r·npj . Fíjese que puesto que el número de pasos total depende de cada conexión, habrá que considerar una distancia de paso, ∆r , lo suficientemente pequeña para que la precisión en media del problema sea aceptable. Figura 4.5 Ejemplo de un segmento de trayectoria. 28 Capítulo 5. Planificación de la Trayectoria Aeq =A0;Ades =      c[1] t−1 c[2] t−1 . . .. . . c[n] t−1       ;(5.10) En este caso, estamos ante un problema de programación lineal entero-mixta (MIPL), ya que la variable x toma valores enteros y la variable y puede tomar cualquier valor real. Además, los algoritmos de programación lineal y de Dijkstra no pueden emplearse en esta ocasión. La razón es que se han incluido nuevas restricciones de desigualdad a través de la matriz Ades , lo que implica que la relajación del problema a variables reales no tenga por qué dar lugar a una solución en variables enteras para x . Este hecho provoca que la resolución del problema tenga que ser abordada mediante la optimización combinatoria. De nuevo y al igual que en el caso anterior, una vez resuelto el problema y obtenido la solución del vector x∗ , se puede reconstruir el camino crítico viendo cuáles de las conexiones permanecen activas, es decir, cuáles de ellas poseen un valor unidad. Análogamente, el problema se puede plantear para el pseudotiempo sin más que intercambiar la matriz de costes temporales ctpor la matriz de pseudotiempos cτ. 5.2.1 Método heurístico para la resolución del problema min-max. Al terminar la sección anterior se llegó a la conclusión de que el problema de determinación del camino crítico suponiendo el peor escenario (worst-case problem) tenía que ser resuelto empleando técnicas de optimización combinatoria. A pesar de ser así, uno puede plantearse si existen métodos alternativos que permitan obtener la solución del problema, o bien una aproximación a la solución, utilizando los algoritmos clásicos de programación lineal o de Dijkstra. El motivo es que la complejidad que incorpora la optimización combinatoria hace que los software que utilizan dichos algoritmos sean mucho menos accesibles y normalmente requieren de licencias de pago para el consumidor. Por tanto, se propone una formulación alternativa, empleando un método de carácter heurístico que consta de los pasos que se describen a continuación: 1. Se calcula el camino que minimiza el tiempo de vuelo final para cada uno de los k escenarios propuestos, siguiendo la siguiente formulación: m´ ın xc[k]T tx s.a.Ax =b x∈{0,1}l (5.11) Se tendrán que resolver n problemas diferentes, en los cuales variará la función objetivo, que incluirá el coste de cada escenario. Para obtener la solución se emplean los algoritmos de Dijkstra o de programación lineal. El resultado de cada uno de los problemas será un camino crítico x[k]∗ . Además, se obtiene también el camino crítico solución de minimizar el coste promedio, el cual se identifica con un escenario ficticio n+1. 2. Suponiendo que la aeronave sigue cada camino x[k]∗ , se extrae el tiempo de vuelo final máximo de entre todos los escenarios. Para ello, se multiplica la matriz de costes por cada uno de los vectores x[k]∗ , obteniéndose un vector de tiempos finales de n componentes (una para cada escenario), de las cuales se selecciona el escenario de mayor tiempo final m´ axcT tx[k]∗ . 5.3 Minimización del combustible promedio final 29 Una vez realizada esta operación para todos los caminos, se dispondrá de un vector de tiempos finales máximos: tf,max = [tf1,tf1, ... ,tfn+1]. 3. El camino solución de nuestro problema, x∗ será aquel camino x[k]∗ , que nos proporcione el menor de los tiempos de vuelo máximos recogidos en el vector de tiempos anterior, t∗ f= m´ ın[tf1,tf1, ... ,tfn+1]. Es importante destacar que el camino que se obtenga a través del método heurístico no tiene por qué ser el camino óptimo para el worst-case problem, ya que en ningún momento se está abordando la resolución del problema directo a través de la formulación propuesta en la sección anterior. Sin embargo, es cierto que si el número de escenarios es suficientemente elevado y la dispersión entre los mismos es pequeña, puede coincidir con la solución óptima (aunque no se tenga ninguna garantía de ello). Como inconveniente, hay que decir que mediante este procedimiento el número de problemas a resolver es mucho más elevado que mediante la optimización combinatoria. A pesar de que cada uno de los problemas individuales es más simple desde el punto de vista computacional, como se verá en el siguiente capítulo, la resolución de todo el problema completo llevará aparejado un mayor tiempo de ejecución. 5.3 Minimización del combustible promedio final En esta sección se aborda la elección de una trayectoria que minimice el combustible final. La importancia de la reducción del combustible es clave para las aerolíneas no sólo porque reducen la carga de emisiones a la atmósfera, sino porque el simple hecho de reducir el combustible consumido, implica una reducción del coste operativo directo de la aeronave al realizar ese vuelo. Uno podría pensar que minimizar el tiempo de llegada al destino final, tiene como consecuencia una reducción del combustible consumido. Obviamente, este hecho es cierto cuando nos planteamos reducir el combustible final para un escenario meteorológico concreto. Sin embargo, el camino óptimo que minimiza la media del combustible total consumido no tiene por qué coincidir con el que minimiza la media del tiempo (o pseudotiempo), debido a la relación no lineal entre ambas variables. Por tanto, no es mala decisión intentar resolver de manera directa los problemas que minimizan el combustible consumido en términos medios. Para obtener una estimación más precisa del combustible consumido se utilizará el pseudotiempo en lugar del tiempo. Como ya se ha comentado, la incorporación del pseudotiempo lleva asociado un modelo más realista de consumo específico que contempla las variaciones de dicho consumo con la temperatura en cada uno de los puntos del espacio recorrido. Por otra parte, la expresión de la masa consumida con respecto al pseudotiempo se había obtenido cuando se demostraron las ecuaciones, resultando en una expresión no lineal como la siguiente: Φ(τf) = mi−sˆ A ˆ Btan   arctan sˆ B ˆ Ami −pˆ Aˆ Bτf   (5.12) Los problemas no lineales pueden ser muy difíciles de resolver. Por ello, se llegará a un compromiso entre recoger parte de la no linealidad del problema y facilitar en la medida de lo posible la resolución del mismo. Esto pasa por realizar una aproximación cuadrática de la función que permita su posterior tratamiento mediante algoritmos de programación cuadrática. Φ(τf)≃Φ(τf,0)+Φ0(τf,0)(τf−τf,0)+ 1 2Φ00(τf,0)(τf−τf,0)2(5.13) 30 Capítulo 5. Planificación de la Trayectoria Las derivadas primera y segunda se calculan a partir de la función principal Φ(τf). Φ0(τf) = ˆ A· 1+tan2   arctan sˆ B ˆ Ami −pˆ Aˆ Bτf    Φ00(τf) = −2pˆ Aˆ B·Φ0(τf)·tan   arctan sˆ B ˆ Ami −pˆ Aˆ Bτf   (5.14) El valor de τf,0 , es el punto donde se realiza la aproximación cuadrática. Se puede tomar como aproximación el resultado del problema de minimización de la media de pseudotiempos. También es necesario conocer una aproximación de la masa inicial del vuelo, mi . Obtener este dato en ocasiones, puede no ser sencillo, ya que no suele ser un dato que las aerolíneas y las compañías aéreas den a conocer por motivos de pura competitividad. A veces, es más común conocer la masa final de la aeronave en el destino mf y realizar un cálculo inverso para poder conocer la masa inicial. Se podría incluso realizar un cálculo previo mediante simulaciones de Monte-Carlo con las predicciones de vientos que se nos proporcionan. Sin embargo, como el cálculo de la masa inicial no es objeto de este trabajo, se asumirá como un valor conocido a priori. Para cada uno de los escenarios la función Φ puede escribirse en función del camino x y de su aproximación x0. Φ(c[k]T τx)≃Φ(c[k]T τx0)+Φ0(c[k]T τx0)(c[k]T τx−c[k]T τx0) +1 2Φ00(c[k]T τx0)(c[k]T τx−c[k]T τx0)2(5.15) Para simplificar la notación y su posterior resolución se introducirán nuevas variables y[k]=c[k]T τx . Además las aproximaciones a la función Φse escribirán con el subíndice 0. Φ(y[k])≃Φ(y[k] 0)+Φ0(y[k] 0)(y[k]−y[k] 0) +1 2Φ00(y[k] 0)(y[k]−y[k] 0)2(5.16) Si se desea obtener el camino crítico que minimice la media de la masa consumida final, se deben de cumplir las restricciones de conectividad del grafo, así como las que definen las nuevas variables intermedias introducidas. m´ ın x,y[k] 1 n n ∑ k=1 mF,k s.a.Ax =b c[k]T τx=y[k]k∈{1,2,3... n} x∈{0,1}l y[k]∈R (5.17) 5.3 Minimización del combustible promedio final 31 Por otro lado, la media de la masa de combustible puede desarrollarse en función de Φ de la siguiente forma: m´ ın 1 n n ∑ k=1 mF,k=m´ ın 1 n n ∑ k=1 Φ(y[k]) = m´ ın(1 n n ∑ k=1Φ[k] 0+Φ0[k] 0(y[k]−y[k] 0) +1 2Φ00[k] 0(y[k]−y[k] 0)2=m´ ın(1 n n ∑ k=1Φ[k] 0−Φ0[k] 0y[k] 0+1 2Φ00[k] 0y[k]2 0 +1 n n ∑ k=1Φ0[k] 0−Φ00[k] 0y[k] 0y[k]+1 n n ∑ k=1 1 2Φ00[k] 0y[k]2) La expresión anterior consta de tres términos: un término constante, el cual puede ser ignorado en la optimización; un término lineal en y[k]y un último término cuadrático en y[k]2. Almacenando las variables de decisión en un vector w=x,y[1],y[2], ... ,y[n]T , el problema puede ser escrito como: m´ ın wcTw+wTQw s.a.Aeqw=beq (5.18) donde los vectores y matrices c , beq , Aeq y Q recogen las restricciones y la función objetivo anterior. c=       0 1 nΦ0[1] 0−Φ00[1] 0y[1] 0 . . . 1 nΦ0[n] 0−Φ00[n] 0y[n] 0        ;beq =       b 0 . . . 0        ; Aeq =         A0 c[1]T τ−1 . . .... c[n]T τ−1         ;Q=          0 0 1 2 Φ00[1] 0 n 0... 1 2 Φ00[n] 0 n          ; (5.19) En vista a la definición de las matrices cabe destacar varios aspectos que son importantes. Por un lado, la introducción de las variables y[k] permite que la matriz Q de tamaño (l+n)×(l+n) , sea en gran medida dispersa. Este hecho facilita tanto la resolución del problema como los potenciales problemas de memoria a la hora de almacenar la matriz. Por otro lado, la formulación obtenida se puede particularizar para el caso en que se asuma una aproximación lineal de la masa de combustible. En ese caso, la única modificación consistiría en eliminar la matriz Q y las derivadas segundas Φ00[k] 0. En esta ocasión nos encontramos ante un problema de programación cuadrática (MIQP), que únicamente puede ser abordado mediante técnicas de optimización combinatoria. Una vez resuelto y obtenido el valor óptimo de w∗ , se podrá calcular simultáneamente tanto el camino a seguir por la aeronave x∗ ; como el consumo de combustible final óptimo m∗ F , evaluando la función Φ en el pseudotiempo final para cada escenario y[k]∗=c[k]T τx∗. 32 Capítulo 5. Planificación de la Trayectoria 5.4 Variaciones del coste con el horizonte temporal Hasta el momento hemos supuesto que los costes han sido obtenidos para un cierto horizonte temporal, sin especificar en ningún momento para qué instante temporal se han calculado. A raíz de esto, se pueden plantear ciertas dudas a la hora de planificar la trayectoria. ¿En qué instante temporal debo seleccionar los costes: a la hora de despegue, a la de aterrizaje o en un valor medio a lo largo de la trayectoria?. ¿Son importantes en términos relativos las variaciones de los costes a lo largo del tiempo?. Para responder a estas preguntas se han calculado los costes para la trayectoria entre Málaga y Arlanda para distintos tiempos de referencia coincidentes con los instantes en los que se nos proporcionan las predicciones (+0 horas,+6 horas,+12 horas,+18 horas,+24 horas,+30 horas). Para los tiempos intermedios se han interpolado linealmente los costes para cada una de las conexiones. Los resultados para dos conexiones aleatorias se muestran a continuación: 0 5 10 15 20 25 30 tref [h] 154 156 158 160 162 164 166 168 170 172 174 ct [s] Coste en t vs tref 0 5 10 15 20 25 30 tref [h] 310 315 320 325 330 335 340 345 350 ct [s] Coste en t vs tref Figura 5.1 Costes vs tre f . A la vista de estas gráficas, se pueden extraer una serie de conclusiones. En primer lugar, se ve que la dispersión entre los escenarios es pequeña, y aumenta conforme aumenta el tiempo de referencia. Esto es esperable, porque a mayor horizonte de predicción la incertidumbre recogida por los EPS debe de aumentar. Por último, la conclusión más importante es que la variación del coste con los escenarios es menor que la propia variación de los costes con el tiempo en un entorno de tiempos de referencia, el cual, abarque la posible duración de la trayectoria anteriormente citada. Las conclusiones anteriores ponen de manifiesto que es necesario intentar incluir la variación temporal de los costes en el problema de optimización. Dichas variaciones pueden provocar cambios en el camino óptimo. El efecto del tiempo, hace que la complejidad del problema crezca enormemente. El motivo es que buscar un camino que minimice unos costes que son función a su vez del tiempo en el que se llega a una conexión, no es nada fácil. Por ello, en las sucesivas subsecciones se desarrollarán alternativas y estrategias para dar respuesta a este problema. 5.4.1 Introducción explícita del tiempo de vuelo La primera de las metodologías que se proponen, consiste en introducir de manera explícita el tiempo en el problema de optimización. Una forma de realizarlo es mediante la introducción de un vector compuesto por nuevas variables de decisión t[k] , que representan los tiempos de llegada a cada nodo i del grafo para el escenario k correspondiente. Por lo tanto, esto supone un incremento en el número de variables de m×n (nodos × escenarios). También se va a definir una nueva matriz 5.4 Variaciones del coste con el horizonte temporal 33 Aorg, de la siguiente forma: Aorg =     001··· 100··· . . .. . .. . .. . . 010···      (5.20) La matriz Aorg es una matriz similar a la matriz A que contiene las restricciones de conectividad del grafo (nodos × conexiones), pero cambiando el valor -1 de los nodos de destino de cada conexión, por el valor 0. De esta manera, la matriz Aorg , sólo tiene en cuenta el origen de los nodos de cada conexión. Esta matriz permite calcular el tiempo en el que se llegará al nodo origen de cada conexión, sin más que multiplicar la matriz traspuesta por el vector de tiempos: AT orgt[k]. Por otro lado el coste ct ahora depende del tiempo. Además, se asumirá que el coste para cada escenario, c[k] t , sólo depende del tiempo de llegada correspondiente a ese escenario c[k] t(t[k]) y no del tiempo para otros escenarios. Con el objetivo de encontrar una función que nos permita modelar los costes, se realizará una aproximación lineal de los mismos, alrededor de un tiempo t0. c[k] t(t[k]) = ct(t0)[k]+dct(t0) dt [k] AT org(t[k]−t0)(5.21) Una vez presentada la expresión es necesario aclarar varios aspectos. En primer lugar, el tiempo t0 , se escogerá como un tiempo previsiblemente cercano al tiempo de llegada a cada nodo en la realidad, por ejemplo como di−O Vre f , es decir, la distancia ortodrómica desde cada nodo al origen entre la velocidad aerodinámica de referencia a la altitud de crucero. Este tiempo será igual para todos los escenarios. Por otro lado, el término dct(t0) dt [k] , es una matriz diagonal de dimensiones l×l (conexiones × conexiones), cuyos elementos de la diagonal representan las derivadas de los costes para cada escenario evaluadas en el tiempo AT orgt0 . El tiempo AT orgt0 , es el tiempo de llegada al origen de cada conexión, según el valor de t0 propuesto. Puesto que conocemos el valor de los costes para varios horizontes temporales concretos (+0 horas,+6 horas,+12 horas,+18 horas,+24 horas,+30 horas), las derivadas se supondrán constantes entre cada intervalo e iguales a la pendiente del tramo correspondiente. Por último, ct(t0)[k] , no es más que el coste para cada conexión y para cada uno de los escenarios, evaluado en el tiempo de llegada AT orgt0 . Nótese que esta expresión es muy útil ya que nos permite conocer el coste de cada conexión sabiendo cuál es el tiempo de llegada a la misma. Modelada la función objetivo, se plantea de forma análoga el problema de minimización del tiempo de vuelo promedio para todos los escenarios. m´ ın x,t[k] 1 n n ∑ k=1 c[k] t(t[k])!T x s.a.Ax =b x∈{0,1}l (5.22) 34 Capítulo 5. Planificación de la Trayectoria Sin embargo, hay que incluir nuevas restricciones para ligar el camino elegido x con los tiempos de llegada a cada nodo t[k]. Esas restricciones son las siguientes: t[k] O=0 diagh−ATt[k]−c[k] t(t[k])ix=0k∈{1,2,3... n}(5.23) La primera de ellas t[k] O=0 , indica que el tiempo para el nodo de partida de la trayectoria debe ser nulo para todos los escenarios. Por otra parte, la segunda, diagh−ATt[k]−c[k] t(t[k])ix=0 , nos dice que si una conexión pertenece al camino recorrido (xj=1), el coste de esa conexión es igual al tiempo que se tarda en recorrer dicha conexión, representado por la operación matricial −ATt[k] . En el caso en el que la conexión no sea recorrida ( xj=0 ), el tiempo de llegada a su origen será indiferente (0=0) . Es bueno observar que esta última restricción es no lineal y eso puede provocar más complejidad a la hora de su resolución. Por tanto, para linealizar la restricción se convierte dicha igualdad en dos restricciones de desigualdad, mediante la introducción de una constante B suficientemente grande. −ATt[k]≤c[k] t(t[k])+B(1−x) −ATt[k]≥c[k] t(t[k])−B(1−x) k∈{1,2,3... n}(5.24) De esta forma se consigue que si una conexión está activa, xj=1 , el término B(1−xj) correspondiente a dicha conexión desaparezca. Con ello, se obliga a que la única posibilidad para que ambas restricciones se cumplan es que −ATt[k]=c[k] t(t[k]) . Aún así, la constante B deberá tener un valor elevado para que en el caso de que xj=0 , ambas restricciones permitan una solución factible para los tiempos de llegada a los nodos que no pertenecen al camino crítico. Por ejemplo, en este trabajo se ha considerado un valor de B=10 ·m´ ax(ct(t0)[k]). Nótese que el número de restricciones de desigualdad a añadir será de 2×l×n (2 × conexiones × escenarios), lo que supone un aumento muy significativo en el número de restricciones del problema. Para disponer de una formulación más cómoda del problema, se trabajará con el incremento de tiempos ∆t[k]=t[k]−t0 en lugar de con t[k] . Si se sustituyen en las restricciones de desigualdad el valor de t[k]y la expresión de los costes c[k] t(t[k]), se obtiene: ∆t[k] O=−t0,O − AT+dct dt [k] AT org!∆t[k]≤ct(t0)[k]+ATt0+B(1−x) − AT+dct dt [k] AT org!∆t[k]≥ct(t0)[k]+ATt0−B(1−x) k∈{1,2,3... n}(5.25) La función objetivo también se puede desarrollar para escribirla en función de las variables de decisión. m´ ın x,t[k] 1 n n ∑ k=1 c[k] t(t[k])!T x=m´ ın x,∆t[k](1 n n ∑ k=1 ct(t0)[k]Tx+∆t[k]TAorg dct dt [k] x!) (5.26) Al ser la función objetivo cuadrática, el problema a resolver es un MIQP. Para facilitar su resolución se agrupan de nuevo las variables de decisión en un vector w=x,∆t[1], ... ,∆t[n]T . Por 5.4 Variaciones del coste con el horizonte temporal 35 tanto, el problema puede ser expresado como: m´ ın wcTw+wTQw s.a.Aeqw=beq Adesw≤bdes (5.27) donde los vectores y matrices se definen por bloques de la siguiente forma: c=     1 n∑n k=1ct(t0)[k] 0 . . . 0      ;beq =     b −t0,O . . . −t0,O      ;Aeq =A0 0Z; bdes =       ct(t0)[1]+B+(ATt0)[1] −ct(t0)[1]+B−(ATt0)[1] . . . ct(t0)[n]+B+(ATt0)[n] −ct(t0)[n]+B−(ATt0)[n]        ;Q=         0 0 Aorg dc dt [1] . . .0 Aorg dc dt [n]         ; Ades =               ˆ B−AT+dc dt [1]AT org0··· 0 ˆ BAT+dc dt [1]AT org0··· 0 . . .. . .. . .. . .. . . . . .. . .. . .. . .. . . ˆ B0 0 ··· −AT+dc dt [n]AT org ˆ B0 0 ··· AT+dc dt [n]AT org               ; (5.28) La matriz ˆ B es una matriz diagonal de dimensión igual al número de conexiones, l , tal que ˆ B=   B ... B    ; mientras que Z tiene dimensiones de n×(m×n) y contiene un valor unidad en cada fila, sólo para el nodo de salida de la trayectoria correspondiente al escenario equivalente a dicha fila. Por ejemplo, si el nodo de salida es el nodo número 1, la matriz será: Z=         1 0 ··· 0 0 0 ··· 0··· ··· 0 0 ··· 0 0 0 ··· 0 1 0 ··· 0··· ··· 0 0 ··· 0 . . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . . . . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . .. . . 0 0 ··· 0 0 0 ··· 0··· ··· 1 0 ··· 0         (5.29) Aunque el problema es cuadrático, puede encontrarse un problema lineal totalmente equivalente. La manera de hacerlo es modificando la función objetivo de tal manera que en lugar de minimizar la media del tiempo de vuelo final, se minimice la media del tiempo de llegada al nodo de destino 36 Capítulo 5. Planificación de la Trayectoria de la trayectoria. m´ ın x,∆t[k] 1 n n ∑ k=1 ∆t[k] D(5.30) Ahora, la única diferencia con el problema anterior es que la matriz Q desaparece y el vector de costes cambia. Suponiendo que el nodo de llegada es el último nodo definido, el vector c se puede escribir como: c=0 0 ··· 1 n··· ··· 0··· 1 nT(5.31) Finalmente, la resolución del problema nos permite conocer no sólo el camino óptimo x∗ , que debe de seguir la aeronave desde el origen O al destino D , sino también el tiempo de llegada t[k] a cada uno de los nodos del camino crítico. Esta formulación que introduce los tiempos de manera explícita es bastante completa desde el punto de vista teórico. Sin embargo, presenta ciertos problemas cuando se desea implementar de forma numérica. Por un lado, el número de variables crece enormemente con respecto al problema que no tiene en cuenta las variaciones temporales del coste. Así, por ejemplo, para m=6031 nodos, l=14836 conexiones y n=50 escenarios, el número de variables es de 316386, mientras que el de restricciones alcanza las 1489681. Esto supone mayores costes computacionales, dando lugar a un tiempo de resolución inasumible (del orden de días). En segundo lugar y siguiendo el hilo de lo anterior, el hecho de incorporar los tiempos en las restricciones y en la función objetivo hace muy difícil su convergencia para garantizar el óptimo global. La consecuencia es que tenemos que conformarnos con soluciones factibles de las cuales no se tiene ninguna garantía de optimalidad. Por último, aunque muchas de las matrices sean de naturaleza dispersa, su elevado tamaño provoca problemas de memoria. Más destacado es el caso de la matriz Q , que al contener las derivadas del coste, deja de ser semidefinida positiva (PSD). Un carácter no PSD implica utilizar métodos numéricos para problemas no convexos, los cuales necesitan de más memoria para su resolución. 5.4.2 Método iterativo para incluir las variaciones temporales En el apartado anterior se ha concentrado todo el esfuerzo en intentar incluir los tiempos de manera explícita en el problema y se ha explicado todas las consecuencias que eso conlleva. Pues bien, llegados a este punto, se investigará una forma más efectiva para incluir el efecto del tiempo en los costes de cada conexión. La clave está en pensar en una buena primera aproximación de los tiempos de llegada a cada conexión que permita calcular el coste actualizado a dicho instante de tiempo. Una vez seleccionada la aproximación y obtenida la trayectoria pertinente, se actualizan los nuevos tiempos de llegada y se vuelve a realizar una nueva iteración. Este proceso se repite tantas veces como sea necesario para que el camino a seguir converja a una solución, x∗ , que se asumirá que es la óptima buscada. Por tanto, el procedimiento que se va a llevar a cabo para determinar x∗ , consta de los siguientes pasos: 1. Partiendo de una hora asignada para el despegue, tDep , se actualizan los costes de cada conexión para dicho instante de tiempo. c[k] t,Dep =ct(tHt)[k]+ct(tHt+1)[k]−ct(tHt)[k] tHt+1−tHt (tDep −tHt)(5.32) Donde los valores de ct(tHt) y ct(tHt+1) hacen referencia a los costes para los horizontes temporales de predicción anterior y posterior a la hora de despegue. 5.4 Variaciones del coste con el horizonte temporal 37 2. Conociendo la distancia recorrida de cada conexión, rj , se calculan los caminos de mínima distancia entre el origen y cada uno de los nodos del grafo, siguiendo las propias restricciones de conectividad impuestas. El resultado será tantos caminos como nodos contenga el grafo, es decir, x(l×m) . Este paso puede resolverse fácilmente con un algoritmo de Dijkstra, tomando como pesos las distancias de cada conexión. 3. Hallamos el tiempo de llegada a cada nodo para cada uno de los escenarios, ta(m×n) siguiendo su camino correspondiente de mínima distancia. ta=cT t,DepxT(5.33) 4. Se actualizan los costes de cada conexión y de cada escenario en función de su tiempo de llegada. tl=AT orgta+tDep c[k] t,1=ct(tHt)[k]+ct(tHt+1)[k]−ct(tHt)[k] tHt+1−tHt (tl−tHt)(5.34) 5. Se resuelve el problema deseado (Ej: minimización del tiempo promedio o worst-case problem) con los costes actualizados y sin variación temporal, obteniéndose el camino óptimo, x∗. m´ ın x 1 n n ∑ k=1 c[k]T t,1x⇒x∗(5.35) 6. A partir del camino crítico anterior, se actualizan los tiempos de llegada de los nodos que contienen dicho camino óptimo. El resto de nodos que no pertenecen a la solución, continúan con su tiempo de llegada original. El resultado es una nueva matriz de tiempos ta,2(m×n). 7. Con la matriz de tiempos ta,2 , se vuelven a actualizar los costes, c[k] t,2 y se calcula el nuevo camino crítico, x∗ 2. 8. Se repite el proceso tantas veces como sea necesario, hasta que x∗≡x∗ 2. Una vez presentado el método iterativo cabe destacar que es esperable no necesitar realizar ninguna iteración para obtener el camino crítico, x∗ , es decir, que en una segunda iteración x∗ 2 , la solución sea la misma a la anterior (x∗≡x∗ 2) . La explicación es que para la gran mayoría de los nodos que pertenecen al camino crítico, dicho camino crítico contiene a su vez el camino de mínima distancia x para cada uno de los nodos. Por lo tanto, cuando se realice la primera actualización de los tiempos del camino, ta , en realidad coincidirán con los tiempos originales, dando lugar a la misma solución óptima tras resolver el problema de optimización. No hay que olvidar que aunque se haya considerado que la presencia de los vientos pueda alterar el camino crítico, en general, volar siguiendo el camino más corto hasta un destino final será una buena aproximación inicial para el camino de mínimo tiempo y con ello, para el tiempo de llegada a cada uno de los nodos. Es más, suponiendo que el camino de mínima distancia para algunos de los nodos no estuviera incluido dentro del camino crítico, las variaciones al actualizar los tiempos de llegada, serían del orden del coste de la conexión. Y teniendo en cuenta que disponemos de predicciones de costes cada 6 horas, las variaciones del coste ∆ct relativas al intervalo de predicción, serían tan pequeñas que es previsible que no alteren el camino crítico. 44 Capítulo 6. Aplicación práctica. Resultados 0 0.5 1 1.5 2 2.5 3 τf×104 0 0.5 1 1.5 2 2.5 3 3.5 mF ×104 mF vs τf Aircraft B763 (mi=153808 kg) real linear aprox quadratic aprox Figura 6.2 Masa Consumida vs Pseudotiempo. A la vista de la expresión 5.12, la función Φ(τ) aparentaba ser altamente no lineal. Sin embargo, se observa que no es así. El comportamiento de la función real (curva azul), es siempre creciente y prácticamente lineal. Por ese mismo motivo, las aproximaciones lineal y cuadrática (curva roja y verde), realizadas en torno a τf,0=¯ τf , coinciden en gran exactitud con la función original. Dado que la función real es siempre creciente y muy lineal, se cumplirá que minimizar el valor medio de Φ(τ) , será equivalente a minimizar el valor medio de su argumento, τ . Por tanto, el camino crítico que minimiza el consumo promedio es exactamente el mismo que el que minimiza el pseudotiempo promedio. m´ ın x,y[k] 1 n n ∑ k=1 mF,k≡m´ ın x 1 n n ∑ k=1 c[k] τ!T x(6.3) Esto supone una gran ventaja esencialmente por dos motivos. En primer lugar, se ha encontrado una variable equivalente τ , cuyo comportamiento es similar al de la masa de combustible consumida mF (siempre que el vuelo se realice a la misma altitud y número de Mach). Eso implica que para conocer como se comporta mF , no es necesario tener un modelo aerodinámico del avión en cuestión, sino simplemente nos basta con el propulsivo. En segundo lugar, los problemas para minimizar el pseudotiempo promedio son problemas mucho más sencillos que los problemas MIQP o MILP originales. Por ello, el problema puede resolverse sin más que aplicar algoritmos clásicos como el de Dijkstra para calcular el camino que minimiza dicho pseudotiempo promedio. Finalmente, se calcula el consumo medio de combustible, siguiendo el camino anterior. Tabla 6.4 Pseudotiempo medio, masa de combustible y tiempo de ejecución del problema lineal. ¯ τf[s]mF[kg]T.Ejecución [s] 15898 18534 1.29 6.3 Método Heurístico vs Enfoque Combinatorio 45 Con este procedimiento se ha conseguido una reducción en el tiempo de ejecución del orden del 70% para obtener la misma solución que en el problema original. 6.3 Método Heurístico vs Enfoque Combinatorio En esta nueva sección, nos centramos en resolver el worst-case problem por dos métodos diferentes. En primer lugar, se obtendrá la solución de dicho problema de forma directa, empleando un enfoque combinatorio. En segundo lugar, se hará uso de un método heurístico mediante el algoritmo de Dijkstra. Finalmente, se realizará una comparación entre las dos propuestas y se analizará si es posible extraer alguna conclusión de los resultados obtenidos. El vuelo y la predicción meteorológica coinciden con los que se han utilizado en apartados anteriores. En este caso, tampoco se incluyen las variaciones temporales en los costes de la función objetivo. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (a) Posibles Caminos. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (b) Heurística. 46 Capítulo 6. Aplicación práctica. Resultados 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (c) Combinatoria. Figura 6.3 Resultados del Método Heurístico y de la Optimización Combinatoria. El camino óptimo que minimiza el peor de los tiempos de entre todos los escenarios, obtenido por el método heurístico, se observa en la figura anterior. En la subfigura a) se ven todos los posibles caminos resultados de aplicar el método heurístico para cada uno de los escenarios definidos. Cabe recordar, que la heurística requería de resolver tantos problemas como escenarios; y posteriormente, seleccionar como óptimo el camino de menor tiempo máximo entre todos los posibles caminos obtenidos. En este caso, el número de escenarios es de n=51 (50 escenarios más 1 escenario ficticio). Sin embargo, el número de caminos posibles es mucho menor a este valor (un total de 7 diferentes). Esto quiere decir que para minimizar el tiempo de vuelo de cada escenario, la solución es la misma para muchos de ellos. La subfigura b) muestra el camino óptimo correspondiente al de menor tiempo de vuelo máximo. Por otra parte, en la subfigura c) se emplea un enfoque combinatorio, es decir, se resuelve el worst-case problem, encontrándose la solución de forma directa. En este caso, el método heurístico encuentra la solución óptima, ya que ambos caminos óptimos coinciden. Sin embargo, de forma general no tendría por qué ser así. Únicamente mediante la resolución directa, utilizando un enfoque combinatorio, se tiene garantía para asegurar que dicho camino obtenido, será el que minimice el peor de los tiempos de entre todos los escenarios. Los tiempos de ejecución de cada uno de los algoritmos se muestran en la siguiente tabla. Tabla 6.5 Comparación Heurísitca y Combinatoria. HEURÍSTICA COMBINATORIA m´ axtf[s]Tiempo Ejecución m´ axtf[s]Tiempo Ejecución 16147 ≃30 s16147 ≃6s 6.4 Influencia de las variaciones temporales del coste 47 El método heurístico requiere de un mayor tiempo de ejecución que el enfoque combinatorio (del orden de 5 veces más). La razón es que aunque el tiempo de resolución de cada uno de los problemas de los que consta la heurística es inferior al del problema combinatorio, el hecho de tener que resolver un total de 50 problemas (uno para cada escenario), hace que el tiempo final sea mucho más elevado. Como conclusión, en vista a los ínfimos tiempos de computación del método combinatorio; y teniendo en cuenta su garantía al proporcionar la solución óptima, es mucho más aconsejable utilizar un enfoque combinatorio antes que un método heurístico. 6.4 Influencia de las variaciones temporales del coste Hasta este momento, todavía no se han introducido las variaciones del coste con el horizonte temporal, ni se ha estudiado la influencia de dichas variaciones en el camino óptimo. En esta sección nos ocuparemos de dicho asunto. En primer lugar, considérese la predicción para el día 11/01/2019 a las 12:00. Las trayectorias óptimas que se obtienen de minimizar el tiempo de vuelo promedio, tomando los costes para las horas de despegue de +18 horas y +24 horas, con respecto al instante de la predicción, son las siguientes: 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (a) tDep = +18 h. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (b) tDep = +24 h. Figura 6.4 Caminos óptimos sin dependencia temporal de los costes. A la vista de las figuras, vemos que los dos caminos críticos son bastante diferentes. Por lo tanto, los costes 6 horas después difieren lo suficiente para provocar cambios significativos en el camino crítico. Nos preguntamos ahora, cuál sería el camino óptimo para un instante intermedio, por ejemplo, para un tDep = +22 horas . Si interpolamos linealmente los costes para dicha hora de despegue, se puede calcular su trayectoria óptima desde el origen al destino. 48 Capítulo 6. Aplicación práctica. Resultados 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N Figura 6.5 Camino óptimo sin dependencia temporal de los costes: tDep = +22 h. La solución óptima coincide con el camino crítico para +18 horas. Sin embargo, a medida que se va recorriendo la trayectoria el tiempo irá avanzando. Por ejemplo, pasadas 2 horas del despegue, el avión estará atravesando Francia. Los costes en dicha región se parecerán más a los que se tienen a +24 horas que a +18 horas, por lo que la trayectoria óptima podría ser la correspondiente a +24 horas u otra algo distinta. Por tanto, existe la duda acerca de la fiabilidad de la solución obtenida sin introducir la variación temporal en la función objetivo. Ahora bien, si aplicamos el método iterativo presentado en el capítulo anterior para introducir dichas variaciones temporales en el problema, el resultado anterior cambia completamente. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N Figura 6.6 Camino óptimo con dependencia temporal de los costes: tDep = +22 h. Al incorporar la variación temporal en los costes, la solución óptima coincide con la que se tiene a +24 horas. Por otra parte, los tiempos de vuelo se observan en la siguiente tabla. En la columna de la izquierda, se muestra el tiempo de vuelo para el caso en el que se decide volar la trayectoria obtenida de minimizar sin dependencia temporal asumiendo variaciones temporales de los costes, 6.5 Efecto de la consideración de regiones D-P-R activas 49 mientras que en la columna de la derecha, se puede ver el tiempo correspondiente a la trayectoria óptima resultante de minimizar teniendo en cuenta las dependencias temporales. Tabla 6.6 Tiempos de vuelo para tDep = +22 h. SIN DEPENDENCIA TEMPORAL CON DEPENDENCIA TEMPORAL Tiempo de Vuelo Tiempo de Vuelo 14561 s14549 s Aunque el camino crítico sufre variaciones apreciables, se observa que el tiempo de vuelo sólo varía un 0.08% por el hecho de incluir la dependencia temporal del coste en la optimización. En cuanto al método iterativo, cabe decir que ha convergido en la primera iteración, por lo que todos los resultados mostrados se corresponden con la iteración 0. Tal y como se dijo cuando se explicó el método iterativo, era esperable obtener la misma solución tras la primera iteración. El motivo era que considerar como primera aproximación que el tiempo de llegada a cada nodo coincide con el camino más corto para llegar al nodo (recorriendo el grafo), es una muy buena aproximación. Y que en el caso de que existan desviaciones con respecto al camino de mínima distancia, esas desviaciones no alterarán de forma significativa los costes, los cuales varían linealmente (según lo propuesto) en el entorno de cada 6 horas, ya que es la ventana temporal en la que se dispone de la información acerca de la meteorología. Como conclusión, decir que el efecto del tiempo es tan importante que sin duda debe de tenerse en cuenta a la hora de elegir la mejor ruta que debe de volar nuestra aeronave. 6.5 Efecto de la consideración de regiones D-P-R activas Una de las cualidades de esta herramienta planificadora consiste en ser capaz de tener en cuenta las restricciones del espacio aéreo. Para ello, como ya se explicó con anterioridad se realiza un filtrado de todos los nodos o conexiones que están dentro o que atraviesan dichas regiones restringidas al vuelo. Para demostrar el funcionamiento de la herramienta se han definido las siguientes áreas D-P-R. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N Figura 6.7 Espacio Aéreo Estructurado. Regiones D-P-R. 50 Capítulo 6. Aplicación práctica. Resultados En este ejemplo se han utilizado los límites aproximados de dos pequeñas regiones: LED47B y LED104, situadas al norte y al sur de Aragón respectivamente. Para tener en cuenta alguna región D-P-R fuera de España, se ha inventado a modo de ejemplo, una gran región en Centro-Europa, simplemente para verificar que la herramienta realiza su cometido. Sin embargo, la herramienta está capacitada para recibir la información ya automatizada acerca de las fronteras de dichas áreas e integrarlas dentro del filtro que elimina los elementos del grafo correspondiente. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (a) Regiones D-P-R Inactivas. tDep = +24 h. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (b) Regiones D-P-R Activas. tDep = +24 h. Figura 6.8 Caminos óptimos con regiones D-P-R I. En la figura de la izquierda se tiene el caso en el que las regiones no están activas. Se observa perfectamente que el camino óptimo cruza dos de las áreas restringidas en verde. En la figura de la derecha, se procede a activar dichas regiones, de tal manera que el camino se redirige para no pasar por dichas áreas prohibidas. En ambos casos, se ha minimizado el tiempo de vuelo promedio para una hora de despegue de 24 horas a partir del instante en el que se toma la predicción. Por otro lado, se puede calcular la trayectoria teniendo en cuenta las variaciones temporales en el coste. El camino óptimo resultante es distinto al caso en el que no se consideran los efectos temporales en la función objetivo. 6.6 Influencia de establecer un límite superior al ángulo de viraje 51 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (a) D-P-R Activas. Sin Variación Temporal. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (b) D-P-R Activas. Con Variación Temporal. Figura 6.9 Caminos óptimos con regiones D-P-R II. En el caso en el que no se tuvieran en cuenta las variaciones temporales, se decidiría por volar la trayectoria óptima de la subfigura a). Pues bien, si fuera así, el tiempo de vuelo medio con los costes reales, sería un 0.61% superior al caso óptimo obtenido en la subfigura b). Si se compara esta variación en los tiempos de vuelo con la variación obtenida en la sección anterior para + 22 horas (0.08%), se observa que la presencia de regiones D-P-R, parece potenciar bastante el efecto de las variaciones temporales del coste. Por lo tanto, se vuelve a dejar constancia en este ejemplo de la importancia de las dependencias temporales a la hora de elegir la ruta adecuada para nuestro trayecto. Tabla 6.7 Comparación de los tiempos de vuelo al incluir las variaciones temporales del coste. TRAYECTORIA A) TRAYECTORIA B) Tiempo de Vuelo Tiempo de Vuelo 15528 s15433 s 6.6 Influencia de establecer un límite superior al ángulo de viraje Es relevante en esta ocasión, fijarnos en ciertos detalles de algunas de las trayectorias que se han estado mostrando anteriormente. Por ejemplo en las trayectorias de la Figura 6.8, se observan ciertas esquinas muy pronunciadas en las que el avión debe de ser capaz de cambiar su rumbo un ángulo muy elevado para alcanzar el curso del segmento que debe seguir. Como ya se explicó en la sección 5.5, se puede limitar dicho ángulo, proporcionando trayectorias más suaves que permitan un mayor confort para el pasajero y para los tripulantes. En la Figura 6.10 a) se obtiene una trayectoria con un cambio de ángulo de curso máximo de 50.89 ° . Si se desea disminuir ese valor se puede limitar el máximo ángulo de viraje a un valor inferior, por ejemplo a 50 ° . El resultado se observa en la Figura 6.10 b), donde ahora el cambio de rumbo máximo pasa a ser de 30.09 ° . En este caso, el camino crítico ha sido modificado por el 52 Capítulo 6. Aplicación práctica. Resultados algoritmo para obtener el valor mínimo del tiempo medio que cumple con la restricción anterior. Si se continúa imponiendo un límite aún más bajo, por ejemplo 30 ° , se tiene el camino óptimo de la Figura 6.10 c). Todas las trayectorias han sido calculadas para una hora de despegue 24 horas posterior a la predicción del 11/01/2019 a las 12:00. Por otro lado, los límites al ángulo de viraje sólo han sido aplicados para nodos situados a una distancia mayor de 500 km del origen y del destino. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (a) Sin limitación de ángulo: ∆χmax =50.89◦. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (b) ∆χlim =50◦;∆χmax =30.09◦. 6.6 Influencia de establecer un límite superior al ángulo de viraje 53 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N (c) ∆χlim =30◦;∆χmax =27.32◦. Figura 6.10 Limitación del ángulo de viraje. Este proceso puede continuar hasta que se alcanza un valor de ∆χmax , a partir del cual, el modelo deja de ser factible. Es decir, no existe ningún camino posible para viajar del origen al destino verificando la restricción impuesta en el ángulo de viraje. En ese momento, se habrá obtenido el mínimo valor de ∆χ que el avión debe ser capaz de virar durante la trayectoria. En este ejemplo, dicho ángulo es el siguiente. m´ ın∆χmax =26.4◦(6.4) 6.6.1 Comparación entre el método Dijkstra modificado y el enfoque combinatorio Como ya comentamos en el capítulo anterior, una de las opciones para limitar el ángulo de viraje podría ser utilizar un algoritmo de Dijkstra modificado. Esa modificación consistiría en eliminar algunos de los nodos vecinos del nodo explorado. Concretamente, se eliminarían aquellos nodos que tienen conexión con el nodo explorado y cuyo ángulo de viraje necesario para alcanzarlos fuera mayor a la restricción impuesta. De esta manera, se conseguiría cumplir con la limitación del ángulo, al igual que se hacía resolviendo el problema directo mediante los algoritmos combinatorios. En primer lugar, se puede empezar imponiendo un ∆χlim =50◦ . Si se realiza esto, se observa que las trayectorias de ambas soluciones coinciden y están representadas en la Figura 6.10 b). Por lo tanto, Dijkstra podría parecer a priori una buena opción para incorporar dicha mejora. Además, hay que tener en cuenta que la implementación del Dijkstra es mucho más sencilla, y su resolución lleva asociada menor tiempo de computación que la optimización combinatoria. Sin embargo, veamos qué ocurre cuando se continúa bajando el límite, por ejemplo a ∆χlim =30◦. 60 Capítulo 6. Aplicación práctica. Resultados 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N 27/01/19 12:00 28/01/19 00:00 28/01/19 12:00 minima distancia Figura 6.18 Caminos críticos para el 28/01/2019 a las 19:00. A la vista de los resultados anteriores, se pueden sacar varias conclusiones. Por un lado, la dispersión definida como tf,max −tf,min disminuye cuanto menor es el horizonte temporal. Esto es lógico, puesto que la incertidumbre se reduce al acercarnos a la hora del despegue. Además, la incertidumbre puede hacer que el problema de tiempo promedio difiera del worst-case problem. Por otro lado, se comprueba que la solución es bastante estable para la mayoría de los casos, salvo para los días 12 y 24 de Enero, donde sí que aparecen cambios en el camino crítico al reducir el horizonte temporal de la predicción. Aún así, dichas variaciones en el camino crítico, no implican diferencias muy elevadas en el tiempo de vuelo final. Por último, cabe decir que la solución para el worst-case problem coincide con la de mínimo tiempo promedio cuando los caminos críticos son estables, pero difiere cuando no lo son (días 12 y 24). La explicación es que la meteorología tan cambiante hace que no sólo el camino óptimo no sea estable, sino que incluso modificaciones en el criterio de optimización dan lugar a variaciones en la trayectoria final. Estos resultados han de ser tenidos en cuenta con cierta cautela. La estabilidad obtenida para 5 de los 7 días analizados no implica necesariamente que la solución vaya a ser estable para el resto de días del ciclo de AIRAC. Para comprobarlo, sería necesario calcular la trayectoria para todos los días de los que consta el ciclo y contemplar los resultados obtenidos. Sin embargo, con este estudio preliminar no se busca obtener garantías de estabilidad, sino averiguar qué ocurre con la tendencia de la solución a lo largo del ciclo de AIRAC. En nuestro caso, puesto que para algo más del 70% de los días analizados la solución es estable, sería bastante esperable encontrar una cierta estabilidad en la solución para el resto de días que quedan sin analizar. A modo de conclusión general, hay que decir que siempre es preferible utilizar la predicción más cercana a la hora de despegue ya que se reduce la incertidumbre en el problema. Sin embargo, con una predicción a 31 horas vista se puede disponer de una muy buena estimación tanto del tiempo de vuelo como del camino crítico que debe de seguir la aeronave para llegar a su destino final. 6.8 Impacto del criterio de optimización en un escenario de tráfico 61 6.8 Impacto del criterio de optimización en un escenario de tráfico En las secciones anteriores se han analizado los resultados en cuanto a tiempos de vuelo para un determinado vuelo con un origen y destino concreto. En la realidad, no sólo se tendrá un vuelo sino un conjunto de vuelos durante un período de tiempo. Puesto que no se posee la información del tráfico correspondiente al ciclo de AIRAC que se está estudiando, se realizará una simulación del tráfico durante el día 12/01/2019. La simulación constará de 1000 vuelos aleatorios cuya hora de despegue estará comprendida entre las 6:00 y las 22:00 del día 12/01/2019. Los vuelos tendrán como origen y destino puntos aleatorios del espacio aéreo separados una distancia ortodrómica también aleatoria pero mayor de 1000 km. En todos los casos, la predicción escogida para definir la incertidumbre meteorológica es la del día 11/01/2019 a las 12:00. Para cada uno de los vuelos se puede calcular su camino óptimo, y con ello el tiempo de vuelo promedio y el tiempo máximo. De forma análoga se resuelve el worst-case problem para cada caso. En la siguiente tabla, se muestra el valor del tiempo promedio acumulado para todos los vuelos, así como el valor del tiempo acumulado para el peor de los escenarios. Tabla 6.16 Tiempos de vuelo para 1000 vuelos aleatorios. hhhhhhhhhhhhhhhh h Tiempo Acumulado Criterio MÍNIMO TIEMPO PROMEDIO WORST-CASE ¯ tf[s]8543857 8543910 m´ axtf[s]8545654 8545648 T. Ejecución [s]2500 7500 (a) Mínimo tiempo promedio. 62 Capítulo 6. Aplicación práctica. Resultados (b) Worst-Case. Figura 6.19 Escenarios de tráfico con 1000 vuelos aleatorios. En base a los números anteriores deben de comentarse dos aspectos que son interesantes. Por un lado, se ve que la diferencia entre el tiempo de vuelo medio acumulado y el tiempo de vuelo acumulado para el peor de los escenarios es del orden de 1700 s. Esta diferencia equivale aproximadamente a unos 29 minutos de pérdida de tiempo en el caso de que la meteorología real coincida con el escenario más desfavorable. En cuanto a las diferencias de tiempos entre ambos criterios, como máximo son del orden de 50 s, lo que sugiere que muchos de los caminos críticos coinciden para los dos casos. Por otro lado, es importante apreciar que la diferencia entre el tiempo del peor escenario para el worst-case problem es de sólo 6 s con respecto al tiempo máximo acumulado obtenido para el problema de minimización del tiempo promedio. Es decir, prácticamente no se contempla ninguna mejora del tiempo máximo al utilizar el criterio de worst-case para un número elevado de vuelos. Tampoco se observan grandes penalizaciones en el tiempo de vuelo promedio cuando se resuelve el worst-case problem. Es más, el tiempo máximo para el worst-case problem podría haber sido incluso mayor que el tiempo máximo resultante de minimizar el tiempo promedio. Por ejemplo, considérese que el camino óptimo para el worst-case en un vuelo A hace que el peor escenario sea el 1, mientras que para otro vuelo B, hace que el peor escenario sea el 2. Por lo tanto, al sumar los tiempos para cada vuelo, el escenario con mayor tiempo podría ser el 1, el 2 o cualquier otro de entre los 50 escenarios disponibles. Sólo en el caso en el que para todos los vuelos el peor escenario sea el mismo, se podría asegurar con garantías que dicho tiempo máximo acumulado sería menor en el worst-case problem que en el de mínimo tiempo promedio. Concluir que aunque el worst-case funcione bien desde el punto de vista local, es decir, reduzca el tiempo de vuelo en el caso de la meteorología más desfavorable, no tiene por qué ser mejor elección que el criterio de mínimo tiempo promedio. Quizás el mejor criterio de optimización posible sería minimizar la suma de todos los tiempos promedios para cada vuelo, o bien, minimizar la suma del tiempo máximo de cada vuelo para el peor escenario. Sin embargo, debido a la complejidad que acarrean dichos problemas, quedan pendientes para una posible mejora futura de la herramienta. Otro asunto importante son los tiempos de ejecución. Se observa que los tiempos de ejecución totales son de 2500 s y 7500 s para los casos de minimización del tiempo de vuelo promedio y de minimización 6.9 Influencia de la altitud de vuelo 63 del peor de los escenarios respectivamente. Estos valores equivalen a tiempos de ejecución por vuelo de 2.5 s y 7.5 s. Por lo tanto, la herramienta permite obtener una solución para un escenario de tráfico en un tiempo razonablemente corto de tiempo. 6.9 Influencia de la altitud de vuelo En esta sección se va a estudiar cómo cambia el tiempo de vuelo y el consumo de combustible cuando varía la altitud. El ECMWF sólo es capaz de proporcionar la meteorología para presiones de 200 hPa y de 250 hPa, correspondientes a altitudes de 11784 m y 10362 m. Para altitudes de crucero intermedias se realizará una interpolación lineal de los vientos y de la temperatura que se usarán posteriormente para el cálculo de la trayectoria. Más concretamente, se van a analizar las altitudes de crucero de 10362, 10500, 11000, 11500 y 11784 m. El instante asignado para el despegue es el día 12/01/2019 a las 19:00, por lo que la predicción meteorológica empleada será la del día anterior (11/01/2019) a las 12:00. El avión para el cual se va a calcular el consumo de combustible es el B763 con una masa inicial estimada de 153808 kg. Los parámetros aerodinámicos y propulsivos de la aeronave son los mismos que se definieron en la sección 6.1. A continuación, se presentan en una tabla los resultados de minimizar el tiempo de vuelo promedio, el pseudotiempo promedio y el consumo de combustible para distintas altitudes. Tabla 6.17 Tiempos, Pseudotiempos y Consumo de Combustible en función de la altitud. Altitud [m]¯ tf[s]¯ τf[s]mF[kg] 10362 15609 15592 19390 10500 15608 15594 19234 11000 15599 15597 18731 11500 15543 15532 18288 11784 15543 15498 18090 Se observa que los tiempos y pseudotiempos no varían mucho con la altitud. Entre la altitud mínima y la altitud máxima analizada se tienen variaciones del orden de 100 s. El tiempo de vuelo disminuye con la altitud mientras que el pseudotiempo parece aumentar hasta llegar a la estratosfera, donde luego empieza a caer. Sin embargo, no puede obtenerse ningún patrón claro de variación ya que intervienen varios efectos. Por un lado, la velocidad aerodinámica cae con la altitud provocando aumentos en los tiempos, pero por otro, los tiempos varían por los vientos presentes en cada capa de altitud. Para entender la variación existente en el consumo de combustible con la altitud, es bueno recordar la expresión 4.28 que se obtuvo para calcular dicho consumo de combustible. En esta expresión, la variación con la altitud viene dada a través de tres parámetros: el pseudotiempo, que cambia con la altitud según lo mostrado en la tabla anterior; pero también a través de las constantes ˆ A y ˆ B , puesto que estas constantes dependen de la presión y del consumo de combustible de referencia, los cuales se ven modificados con la altitud, de acuerdo con las ecuaciones 4.26. 64 Capítulo 6. Aplicación práctica. Resultados 1.04 1.06 1.08 1.1 1.12 1.14 1.16 h [m] 104 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84 0.86 0.88 A [kg/s] A vs h (a) Constante ˆ A. 1.04 1.06 1.08 1.1 1.12 1.14 1.16 h [m] 104 1.8 1.85 1.9 1.95 2 2.05 2.1 2.15 2.2 2.25 2.3 B [1/(kg s)] 10-11 B vs h (b) Constante ˆ B. Figura 6.20 Variación de las constantes ˆ Ayˆ Bcon la altitud. La constante ˆ A disminuye con la altitud, porque así lo hace la presión y el consumo específico de referencia. En cambio ˆ B aumenta porque el consumo de referencia disminuye en menor medida de lo que lo hace la presión, que al estar en el denominador, provoca un aumento de la constante. En el tramo de interés, la tendencia de las dos curvas es aproximadamente lineal. Simplemente por el hecho de modificar la altitud, el consumo de combustible se ve afectado por el cambio en las constantes anteriores. Si representamos el consumo frente a la altitud para un pseudotiempo dado (suponiendo que éste no varía con la altitud), el resultado es el siguiente. 6.9 Influencia de la altitud de vuelo 65 1.04 1.06 1.08 1.1 1.12 1.14 1.16 h [m] 104 1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 mF [kg] 104mF vs h f = 15500 [s] Figura 6.21 Consumo de combustible vs altitud para ¯ τf=15500 s. Se aprecia que el consumo de combustible disminuye claramente con la altitud, experimentando una reducción de 1300 kg desde la capa mínima a la capa de altitud máxima analizada. Esa disminución se produce a pesar de que los pseudotiempos puedan aumentar de una capa de altitud a otra, como se ve en la Tabla 6.17. La razón se ha visto anteriormente, y se debe a la importante variación del consumo con las constantes ˆ A y ˆ B . Las desviaciones en el consumo debidas a cambios del orden de 100 s en el pseudotiempo, sólo suponen incrementos o reducciones de 150 kg por encima y por debajo del valor obtenido a cada altitud. Esto representa en torno a un 0.6% del consumo a cada altitud. Por lo tanto, se concluye que las variaciones del pseudotiempo con la altitud afectan menos al consumo que las propias variaciones provocadas por el cambio en las constantes ˆ Ayˆ B. Veamos qué ocurre con los caminos críticos obtenidos para las altitudes de este estudio. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N min tf 10362 m min τf 10362 m minima distancia (a) h=10362 m. 66 Capítulo 6. Aplicación práctica. Resultados 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N min tf 10500 m min τf 10500 m minima distancia (b) h=10500 m. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N min tf 11000 m min τf 11000 m minima distancia (c) h=11000 m. 6.9 Influencia de la altitud de vuelo 67 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N min tf 11500 m min τf 11500 m minima distancia (d) h=11500 m. 20°W 10°W 0° 10°E 20°E 30°E 36°N 42°N 48°N 54°N 60°N min tf 11784 m min τf 11784 m minima distancia (e) h=11784 m. Figura 6.22 Caminos críticos en función de la altitud. Se observa que el camino crítico coincide para ambos criterios de optimización. Además, la trayectoria óptima sufre cambios con la altitud. Concretamente, se aprecia una variación en el camino a partir de los 11000 m. Esto indica que volar a diferentes altitudes puede suponer no sólo diferencias en los tiempos y en el consumo, sino variaciones en el perfil horizontal del crucero. Finalmente, comentar que el objetivo de este análisis es mostrar al usuario que la herramienta es capaz de proporcionar la trayectoria óptima en función de la altitud de crucero seleccionada para el vuelo. Una posible mejora de la herramienta sería la optimización del perfil vertical de la aeronave, 68 Capítulo 6. Aplicación práctica. Resultados de tal manera que sea capaz de suministrar la altitud óptima a la cual recorrer un determinado segmento, con el objetivo de reducir cualesquiera de las variables estudiadas, por ejemplo, el tiempo de vuelo o el consumo. 7 Conclusiones. Futuras Mejoras En este capítulo se resumen las conclusiones que pueden ser extraídas a raíz de los resultados del capítulo anterior. Posteriormente, se indican las posibles vías de mejora y desarrollos que podría experimentar la herramienta de planificación de trayectorias en un futuro. 7.1 Conclusiones Generales En primer lugar, se recogen algunas de las conclusiones y ventajas más importantes que presenta la herramienta abordada en este trabajo con respecto a otras metodologías existentes. Entre ellas, se pueden destacar las siguientes: • Aunque las diferencias en la función objetivo puedan ser pequeñas, los cambios en el camino crítico pueden llegar a ser importantes. • La metodología propuesta mediante la introducción del pseudotiempo, permite optimizar el consumo de combustible sin necesidad de plantear un problema de programación no lineal. • La influencia de las variaciones temporales de los costes son igual o incluso más importantes que las variaciones de los costes entre los diferentes miembros del ensemble. • La herramienta propuesta permite optimizar las trayectorias teniendo en cuenta unos costes estimados de llegada a cada nodo del espacio aéreo. • En presencia de incertidumbre, el enfoque combinatorio permite resolver ciertos problemas con garantías de optimalidad frente a los algortimos clásicos de búsqueda de caminos críticos como el Dijkstra. • Esta herramienta de planificación permite introducir las regiones prohibidas, peligrosas y restringidas al vuelo (D-P-R), así como la limitación del ángulo de viraje de la aeronave. • El camino crítico estimado en el entorno de 30 horas antes del despegue no parece sufrir grandes variaciones con respecto a la trayectoria estimada pocas horas antes de la salida del vuelo. • Criterios de optimización que funcionan bien a nivel de trayectoria pueden no implicar las mismas ventajas cuando se consideran muchos vuelos dentro de una misma región de espacio aéreo. • La altitud de vuelo provoca no sólo cambios en el tiempo y en el consumo de combustible, sino variaciones en el perfil horizontal óptimo de la trayectoria. • La herramienta permite al usuario seleccionar la altitud y el número de Mach adecuados, de acuerdo a sus necesidades y/o restricciones del crucero. 69