Aplicación del post-procesado estadístico de predicciones meteorológicas por conjuntos a la predicción de trayectoria de aeronaves

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 Grado Grado Ingeniería Aeroespacial Aplicación del post-procesado estadístico de predicciones meteorológicas por conjuntos a la predicción de trayectoria de aeronaves Autor: Manuel Gutiérrez Fernández Tutor: Antonio Franco Espín Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo Fin de Grado Grado Ingeniería Aeroespacial Aplicación del post-procesado estadístico de predicciones meteorológicas por conjuntos a la predicción de trayectoria de aeronaves Autor: Manuel Gutiérrez Fernández Tutor: Antonio Franco Espín Profesor Contratado Doctor Interino Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo Fin de Grado: Aplicación del post-procesado estadístico de predicciones meteorológicas por conjuntos a la predicción de trayectoria de aeronaves 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 del Grado 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í, si no 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 cuatro años. De igual manera, no puedo dejar atrás a todo el personal que contribuye al funcionamiento de este centro y que hace posible 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 graduado, escribiendo el trabajo fin de grado. De manera especial, dar las gracias a Antonio Franco Espín, tutor del trabajo, su dedicación, esfuerzo, paciencia y capacidad para transmitirme los conocimientos desarrollados en el trabajo. También a Damián Rivas Rivas, catedrático del departamento, por su implicación en el tema y por darme la oportunidad de empezar el proyecto. Manuel Gutiérrez Fernández Autor del trabajo Sevilla, 2018 III XÍndice Abreviado 7.3 D. Trayectoria 71 Referencias 75 Í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 Calibración de los modelos 9 3.1 Introducción 9 3.2 Predicciones meteorológicas por conjuntos 10 3.3 Evaluación de modelos de predicción 11 3.3.1 PIT histograms 11 3.3.2 Continuous Ranked Probability Score (CRPS) 15 3.3.3 Absolute Error e Ignorance Score (AE/IGN) 17 3.3.4 Resultados de la aplicación de los Scores 18 3.4 Post-procesado estadístico de modelos de predicción 23 3.4.1 Non-homogeneous regression o ensemble model output statistics (NR/EMOS) 23 3.4.2 Calibración de modelos de una variable 24 3.4.3 Calibración de modelos de varias variables 28 Introducción 28 Ensemble Copula Coupling (ECC) 29 4 Predicción de trayectorias 33 4.1 Introducción 33 4.2 Caracterización de la trayectoria 33 4.3 Ecuaciones de la trayectoria 35 4.4 Resolución Numérica de la trayectoria 37 5 Aplicación. Trayectoria sujeta a incertidumbre 41 5.1 Vuelo de estudio 41 5.2 Conjuntos de calibración 43 XI XII Índice 5.3 Resultados de los vientos tras la calibración 50 5.4 Tiempos de vuelo tras el postprocesado 57 5.5 Conclusiones generales 62 6 Futuras Mejoras. Recomendaciones 65 7 Anexos 67 7.1 B. Cálculo de Scores 68 7.2 C. ECC 69 7.3 D. Trayectoria 71 Referencias 75 Notación F Función de probabilidad acumulada de una variable aleatoria. f Función de densidad de probabilidad de una variable aleatoria. Nµ,σ2 Distribución de probabilidad normal de media µ y de varianza σ2. X Variable aleatoria que se distribuye según una determinada función. PIT Probability Integral Transform. µMedia de una distribución aleatoria. σDesviación típica de una variable aleatoria. σ2Varianza de una variable aleatoria. g Función de densidad de probabilidad de una distribución Gumbel. G(µ,σ) Función de probabilidad acumulada para una Gumbel, de parámetro de localización µ , y de factor de escala σ . E[X]Esperanza matemática de una variable aleatoria X. S2Varianza muestral de una variable aleatoria. ˜ XMedia muestral de una variable aleatoria. γConstante de Euler-Mascheroni. CRPS Continuous Ranked Probability Score. AE Absolute Error. IGN Ignorance Score o Ignorance Error. I{} Función de Heaviside. Φ Función de probabilidad acumulada de una distribución normal de media 0 y desviación típica 1. ϕ Función de densidad de probabilidad de una normal de media 0 y desviación típica 1. CDF Cumulative Distribution Function. PDF Probability Density Function. argmin() Argumento que minimiza (). median(F) Mediana de una variable aleatoria con función de probabilidad acumulada, F. log Logaritmo en base 10. ln Logaritmo neperiano. XIII XIV Notación tα Valor de una variable aleatoria con función de densidad de probabilidad t-student correspondiente a una probabilidad α. NR Non-Homogeneous Regression. EMOS Ensemble Model Output Statistics. |X|Valor absoluto de X. y|X Función de probabilidad y condicionada por la variable aleatoria X. XDistribución de una variable aleatoria. [X]Vector X. ECC Ensemble Copula Coupling. F−1 Inversa de una distribución de probabilidad acumulada. λ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. ξMedia muestral de los CRPS. UT Tiempo universal. hAltura. MNúmero de Mach. Índice de Figuras 3.1 Forma de los PIT según el modelo. Los modelos se han generado siguiendo una normal de forma Nµ,σ2. En este caso, las observaciones se distribuyen como N(0,1). Los modelos en bruto tomados para cada caso son: a) N(0,0.5), b) N(0,2)y c) N(0.5,1)11 3.2 Gráficas Ejemplo 1 12 3.3 Gráficas de Normal y Gumbel con 1000 observaciones para Normal Truth 14 3.4 Gráficas de Normal y Gumbel con 1e6 observaciones para Normal Truth 14 3.5 Gráficas de Normal y Gumbel con 1000 observaciones para Gumbel Truth 14 3.6 Gráficas de Normal y Gumbel con 1e6 observaciones para Gumbel Truth 15 3.7 Definición gráfica del CRPS. Dada una observación, por ejemplo de temperatura, el área sombreada coincide con el valor del indicador CRPS 16 3.8 Cálculo de Scores para la Normal Truth siguiendo distintas distribuciones 19 3.9 Cálculo de Scores para la Gumbel Truth siguiendo distintas distribuciones 19 3.10 Media de los Scores con un intervalo de confianza del 95 %. Caso Normal Truth 21 3.11 Media de los Scores con un intervalo de confianza del 95 %. Caso Gumbel Truth 22 3.12 Ejemplo gráfico de lo descrito en el párrafo anterior. En azul, se muestran las diferentes predicciones; en rojo su valor medio; y en verde un valor posible de una observación del día posterior. Se puede ver claramente la necesidad de calibrar el modelo para que la observación no caiga fuera del rango de las predicciones 23 3.13 La figura muestra un modelo de predicción, en este caso continuo, que hay que ajustar para que se corresponda con las observaciones realizadas 24 3.14 Variables de entrada y salida en la calibración de modelos de una variable 26 3.15 La figura muestra los diferentes modelos de temperaturas de 4 zonas cercanas, así como la temperatura para zonas intermedias. En azul vemos un modelo obtenido suponiendo que la interpolación se ha realizado a través de una minimización del CRPS en cada punto de manera independiente. En rojo, se aprecian modelos más suaves de temperatura que tienen más sentido físico. Los modelos se han exagerado mucho para poder entender el error que existe al aplicar el método descrito en esta sección 28 3.16 Variables de entrada y salida en el ECC 29 3.17 Relaciones de presión y temperatura entre Berlín y Hamburgo. Ilustraciones tomadas de [7] 31 3.18 Mapas de temperaturas sobre Alemania. Ilustración tomada de [8] 32 4.1 Ejemplo de segmento de trayectoria 35 4.2 Vientos y temperaturas para cada punto de la malla 35 4.3 Cálculo del viento longitudinal y cruzado 36 4.4 Triángulo de vientos 36 1 2Índice de Figuras 4.5 Interpolación de los valores del viento en la malla 38 5.1 Representación de la trayectoria en proyección cilíndrica equidistante 42 5.2 Representación de los tiempos de vuelo durante la ida para el primer conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales 48 5.3 Representación de los tiempos de vuelo durante la ida para el segundo conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales 48 5.4 Representación de los tiempos de vuelo durante la vuelta para el primer conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales 49 5.5 Representación de los tiempos de vuelo durante la vuelta para el segundo conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales 49 5.6 Representación del ensemble forecast antes de la calibración, tras realizar la calibración individual, EMOS y una vez aplicado el ECC 51 5.7 Representación del ensemble forecast antes de la calibración, tras realizar la calibración individual, EMOS y una vez aplicado el ECC 52 5.8 Representación de las curvas de nivel para el viento zonal. Primer conjunto de entrenamiento 53 5.9 Representación de las curvas de nivel para el viento meridional. Primer conjunto de entrenamiento 54 5.10 Representación de las curvas de nivel para el viento zonal. Segundo conjunto de entrenamiento 55 5.11 Representación de las curvas de nivel para el viento meridional. Segundo conjunto de entrenamiento 56 5.12 Modelo en bruto de los tiempos de vuelo para el recorrido de ida 57 5.13 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de ida. Primer conjunto de entrenamiento 59 5.14 Modelo en bruto de los tiempos de vuelo para el recorrido de vuelta 59 5.15 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de vuelta. Primer conjunto de entrenamiento 60 5.16 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de ida. Segundo conjunto de entrenamiento 61 5.17 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de vuelta. Segundo conjunto de entrenamiento 62 Índice de Tablas 3.1 Resultados de las medias de los diferentes Scores según las tres distribuciones consideradas. Caso Normal Truth 21 3.2 Resultados de las medias de los diferentes Scores según las tres distribuciones consideradas. Caso Gumbel Truth 22 3.3 Distribuciones según variables. Tabla basada en [8] 27 3.4 Tabla de reordenación de las muestras 30 3.5 Orden global de todos los escenarios 30 5.1 Coordenadas geográficas de los diferentes waypoints 41 5.2 Distancias y cursos entre los diferentes waypoints tanto en ida como en vuelta 42 5.3 Tiempos de vuelos previstos para la ida de los 6 primeros meses. Primer conjunto de entrenamiento 44 5.4 Tiempos de vuelos previstos para la ida de los 6 últimos meses. Primer conjunto de entrenamiento 44 5.5 Tiempos de vuelos correspondientes a cada uno de los días 1 y 15 del 2017 según los vientos observados en cada punto de la trayectoria en el recorrido de ida. Primer conjunto de entrenamiento 45 5.6 Tiempos de vuelos previstos y observados para los primeros 14 días del mes de Diciembre. Recorrido de ida. Segundo conjunto de entrenamiento 45 5.7 Tiempos de vuelos previstos para la vuelta de los 6 primeros meses. Primer conjunto de entrenamiento 46 5.8 Tiempos de vuelos previstos para la vuelta de los 6 últimos meses. Primer conjunto de entrenamiento 46 5.9 Tiempos de vuelos correspondientes a cada uno de los días 1 y 15 del 2017 según los vientos observados en cada punto de la trayectoria en el recorrido de vuelta. Primer conjunto de entrenamiento 47 5.10 Tiempos de vuelos previstos y observados para los primeros 14 días del mes de Diciembre. Recorrido de vuelta. Segundo conjunto de entrenamiento 47 5.11 Resultados de la optimización para el primer conjunto de entrenamiento 50 5.12 Coeficientes de correlación entre ambas ciudades para el primer conjunto de entrenamiento 51 5.13 Resultados de la optimización para el segundo conjunto de entrenamiento 51 5.14 Coeficientes de correlación entre ambas ciudades para el segundo conjunto de entrenamiento 52 3 1 Introducción La predicción de trayectorias resulta muy interesante para toda aquella persona implicada en el mundo de la aeronáutica y la navegación aérea. No sólo para ellos y ellas es esencial el conocimiento de esta disciplina si no que también la meteorología juega un papel muy importante, y por tanto, los meteorólogos son imprescindibles para su desarrollo. El objetivo de este trabajo es dar un paso más en la fusión de la meteorología y la aeronáutica. Como numerosos trabajos anteriores del departamento de Ingeniería Aeroespacial y Mecánica de Fluidos de la Universidad de Sevilla, los vientos son el objeto de estudio de este trabajo. La novedad que incluye este proyecto es que los modelos de viento que proporciona el European Centre for Medium-Range Weather Forecasts (ECMWF) no se consideran tal cual para el cálculo de trayectorias. Lo que se intenta es realizar una serie de modificaciones a esos vientos proporcionados por el ECMWF, de manera que se puedan obtener otros modelos de vientos que predigan mejor la realidad. A partir de ello, mediante una serie de técnicas se resuelven problemas de trayectorias que integran los nuevos modelos de viento calculados, comparándose las características de la trayectoria obtenida, con la trayectoria que se obtendría si se utilizasen los valores en bruto del viento del ECMWF. Este estudio ha sido siempre bastante complejo, ya que se trabaja esencialmente con variables aleatorias. Por ello, su cálculo siempre se basará en probabilidades y estimaciones, siendo la matemática y la estadística dos ramas necesarias para la obtención de los nuevos modelos. El campo de la meteorología es un campo muy abierto, por lo que el trabajo es sólo una primera aproximación, que puede completarse y enriquecerse en un futuro, pero que sirve de base para poder entender las técnicas empleadas para resolver el problema, así como las dificultades que éstas conllevan. Cabe comentar simplemente que muchas de las técnicas recogidas en este trabajo han sido estudiadas antes por numerosos científicos y meteorólogos importantes como Tilmann Gneiting ([2]), Thordis L. Thorarinsdottir ([3]), Michael Scheuerer ([4]) y Roman Schefzik ([8]), entre otros. Las hipótesis realizadas por estos científicos se incluyen en este proyecto, como también algunas nuevas relacionadas con nuestro problema en cuestión. Sin embargo, este trabajo introduce, como novedad, una aplicación de tales técnicas a la predicción de trayectoria de aviones. Además de las técnicas contempladas en el trabajo, existen muchas más en lo referente al postprocesado de modelos de viento. Algunas de ellas, aunque no son tratadas aquí, se recogen en la bibliografía ([4] y [8]). 5 12 Capítulo 3. Calibración de los modelos normal Nµ,σ2, cuyos parámetros son N(25,1). -Se saca un valor observado y, mediante una normal N(µ,9). Este proceso se repite 1000 veces y 1e6 veces. 2. Se define un modelo de predicción X , como una normal N(25,10) . Este modelo se caracteriza por una distribución de probabilidad acumulada F. 3. Se transforman las observaciones mediante la función F, obteniéndose z=F(y). 4. Se lleva a cabo la representación del histograma de las z considerando barras equiespaciadas (en nuestro caso 50). Se obtienen gráficas como las siguientes: 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 Ejemplo 1. Nobs=1000 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 x 104Ejemplo 1. Nobs=1e6 Figura 3.2 Gráficas Ejemplo 1. Se pueden comentar varias cosas. En primer lugar, se ve que con sólo 1000 observaciones el comportamiento es mucho más caótico y no permite sacar ninguna conclusión ya que cada probabilidad tiene una repetición distinta. En cambio, cuando el número de observaciones va aumentando y se va acercando a infinito el comportamiento del modelo se asemeja al de las observaciones. Con un millón, todas las probabilidades se repiten con la misma frecuencia, lo que nos indica que el modelo coincide con el que se han generado las observaciones. Ejemplo 2. En este ejemplo se usará otra función de distribución aparte de la normal, la denominada Gumbel. La función Gumbel es muy utilizada para modelar valores máximos y valores mínimos, por lo que su empleo para predecir terremotos, inundaciones, sequías, etc es habitual. Su función de densidad de probabilidad para modelar el máximo viene dada por g(x;µ,σ) = 1 σexp−x−µ σexp−exp−x−µ σ (3.2) y su función de probabilidad acumulada es igual a G(x;µ,σ) = exp−exp−x−µ σ (3.3) así como para modelar el mínimo simplemente hay que cambiar de signo el argumento de las exponenciales. Antes de proceder con el ejemplo en sí, es necesario definir un método que nos permita, en base a 3.3 Evaluación de modelos de predicción 13 las observaciones, estimar los parámetros de los que dependen sus distribuciones. Existen varios métodos para realizar estimaciones. Los más conocidos son el Método de los Momentos y el Método de Máxima Verosimilitud . En nuestro caso, se utiliza el primero de ellos. Método de los Momentos. Sea (X1,X2,...XN) una muestra aleatoria de una población X con función de probabilidad fθ , con parámetros desconocidos (θ1,θ2,...θk) . Los estimadores de θ se obtienen de resolver el sistema de ecuaciones:      Eθ[X] = 1 N∑N i=1Xi . . .=. . . Eθ[Xk] = 1 N∑N i=1Xk i (3.4) donde Eθ[Xj]se denomina momento de orden j. Simplemente resolviendo el sistema para los dos parámetros de la normal se obtiene, µn=1 N N ∑ i=1 Xi(3.5) σ2 n=s2=1 N−1 N ∑ i=1 (Xi−˜ X)2(3.6) igualmente para la gumbel los parámetros son los siguientes (según [12]) σg=√6 πs(3.7) µg=˜ X−γσg(3.8) siendo γla constante de Euler-Mascheroni γ=0.5772156649015328606. Una vez que se sabe como estimar los parámetros de las distribuciones en cuestión, se está en condiciones suficientes para desarrollar los pasos a seguir en el segundo ejemplo. Esos pasos están detallados a continuación (basados en [3]). 1. La realidad se va a comportar como una normal (Normal Truth), Nµ,σ2 , con µ∼ N(25,1)yσ2=9. 2. Se generan 300 observaciones aleatorias de acuerdo a esa realidad. 3. Se realiza un ajuste de los parámetros de la Normal y de la Gumbel, mediante la estimación basada en el Método de los Momentos, usando las observaciones anteriores. 4. Se simulan otras 1000 o 1e6 observaciones de la realidad. 5. Se transforman las observaciones anteriores mediante la función de probabilidad acumulada de los modelos ajustados anteriormente. Luego se representan los histogramas considerando 50 barras equiespaciadas. 6. Se repiten los pasos del 2 al 5, pero considerando ahora la realidad como una función Gumbel (Gumbel Truth), G(µ,σ), de media µ+γσ ∼N(25,1), y de varianza π2σ2 6=3π2 2. Algunos de los resultados obtenidos se muestran a continuación, tanto para la distribución Normal como para la Gumbel. 14 Capítulo 3. Calibración de los modelos 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 Normal Nob=1000 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60 70 80 Gumbel Nob=1000 Figura 3.3 Gráficas de Normal y Gumbel con 1000 observaciones para Normal Truth. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 x 104Normal Nob=1e6 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7x 104Gumbel Nob=1e6 Figura 3.4 Gráficas de Normal y Gumbel con 1e6 observaciones para Normal Truth. 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 40 Normal Nob=1000 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 35 Gumbel Nob=1000 Figura 3.5 Gráficas de Normal y Gumbel con 1000 observaciones para Gumbel Truth. 3.3 Evaluación de modelos de predicción 15 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 104Normal Nob=1e6 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 x 104Gumbel Nob=1e6 Figura 3.6 Gráficas de Normal y Gumbel con 1e6 observaciones para Gumbel Truth. Cabe comentar varias conclusiones en base a los resultados obtenidos. Por una parte, se observa que cuando la realidad se comporta como una Normal o como una Gumbel, sus respectivas gráficas tienden a un valor uniforme (salvo errores en la estimación del Método de los Momentos). Sin embargo, esto sólo se puede ver cuando el número de observaciones es muy alto, por ejemplo un millón. Las gráficas obtenidas transformando las observaciones por modelos que no coinciden con la realidad, presentan un comportamiento extraño, que nos induce a pensar que no son buenas estimaciones. Es muy importante tener presente, que cuando el número de observaciones es bajo, tanto la Normal como la Gumbel, presentan un comportamiento del que no se puede extraer ningún patrón. Esto nos supone un problema, ya que en la realidad no se dispone de un número alto de observaciones para ajustar los modelos, y a la vista de lo resultados sería bastante difícil diferenciar entre un modelo Gumbel y un modelo Normal. Por tanto, aunque los histogramas PIT permitan ver de forma rápida si un modelo está bien o mal encaminado para representar a una realidad concreta, no es lo más adecuado para realizar buenos ajustes, si no sólo para una primera aproximación. Dicho esto, interesa ahora establecer métricas a la hora de comparar diferentes modelos, así como para comparar los modelos con las observaciones. Surge la necesidad de usar los scores. Los scores no son más que una serie de indicadores definidos con cierto razonamiento lógico que comparan las características de los modelos de predicción con las características de las observaciones, teniendo en cuenta diversos criterios, como pueden ser sus valores, la forma de sus modelos de predicción o cualquier otro, pero siempre traduciendo dicho criterio a un valor cuantitativo de una magnitud definida mediante una fórmula matemática. Aunque existen muchos scores, sólo se mostrarán tres: el Continuous Ranked Probability Score (CRPS), el Absolute Error (AE) y el Ignorance Score (IGN), que se estudiarán en las sucesivas secciones. 3.3.2 Continuous Ranked Probability Score (CRPS) El primero de los scores que definiremos es el CRPS. Su definición (en [5]) es la siguiente. CRPS(F,y) = EF|X−y|−1 2EF|X−X0|(3.9) CRPS(F,y) = Z∞ −∞(F(x)−I{y≤x})2dx (3.10) 16 Capítulo 3. Calibración de los modelos La ecuación (3.9) se utiliza cuando la distribución de probabilidad es discreta, siendo X y X0 variables independientes aleatorias que tienen la misma función de probabilidad acumulada. La ecuación (3.10) hace referencia al caso de una distribución continua y será la que se estudiará en este trabajo. El integrando de la ecuación (3.10) tiene un nombre específico denominado Brier Probability Score. El CRPS trata de medir la diferencia que existe entre la función de probabilidad acumulada del modelo de predicción y una función de probabilidad empírica de salto unidad. La función I{y≤x} es aquella que cuenta o acumula el número de observaciones que son menores o iguales a cada valor considerado de nuestro conjunto de predicción. Su implementación se lleva a cabo mediante la función de Heaviside cuyo valor es unidad si se cumple que y≤x y nulo si no se da esa condición. Figura 3.7 Definición gráfica del CRPS. Dada una observación, por ejemplo de temperatura, el área sombreada coincide con el valor del indicador CRPS. Las ventajas de este indicador es que la función de probabilidad empírica cumple el Teorema de Clivenko-Cantelli. Este teorema dice que si el tamaño de la muestra, n→∞ , entonces la función de probabilidad empírica converge a la función de probabilidad acumulada de nuestro modelo. Por lo tanto, sirve para cuantificar el parecido de nuestro modelo a las observaciones. Cuanto menor sea el valor del score más acertada será la predicción. A pesar de ser un indicador muy utilizado presenta varios problemas a la hora del cálculo numérico. Uno de ellos, es que el intervalo de integración no es acotado (−∞,∞) y hay que definir subintervalos para la integración numérica en los que el peso del integrando no sea significativo. Surge la idea de integrar en probabilidad, pero para ello hay que hacer un cambio de variable en el integrando, el cual tiene ahora el nombre de Quantile Probability Score. CRPS(F,y) = 2Z1 0 (F−1(τ)−y)(I{y≤F−1(τ)}−τ)dτ(3.11) Para trabajar con vientos, presiones o temperaturas, es muy habitual que se modelen dichas variables mediante distribuciones normales. Debido a que la distribución normal es muy utilizada es lógico 3.3 Evaluación de modelos de predicción 17 considerar una expresión analítica que permita calcular el CRPS, sin necesidad de hacer la integral. Esta expresión es la que sigue (en [1]) : crps[Nµ,σ2,y] = σy−µ σ2Φy−µ σ−1+2ϕy−µ σ−1 √π(3.12) siendo Φ y ϕ la CDF y la PDF, respectivamente, de una normal de media 0 y desviación típica 1, evaluadas en y−µ σ. La media del CRPS será por tanto CRPS =1 n n ∑ i=1 crps[Nµ,σ2,yi](3.13) con nel número de observaciones realizadas. De esta manera, la sencillez de cálculo numérico se reduce notablemente dando más rapidez en su implementación. 3.3.3 Absolute Error e Ignorance Score (AE/IGN) Otro score muy sencillo de implementar es el Absolute Error (AE). En primer lugar, se define la siguiente función de Score del Absolute Error (en [3]) como: S(x,y) = |x−y|(3.14) Sin embargo, cuando tenemos una función de probabilidad acumulada basada en un modelo de predicción, hay que elegir con respecto a que x vamos a medir ese error. La teoría dice que para que el score sea consistente debe de calcularse en base al estimador de Bayes definido por, ˆx=argminxEF[S(x,Y)] (3.15) que para el Absolute Error de una función de probabilidad acumulada F , ese valor coincide con su mediana. Por tanto la definición del AE inducida de su función de Score es la que sigue. AE(y) = |median(F)−y|(3.16) Para terminar esta sección, queda explicar el último de los scores que se utilizará en los ejemplos sucesivos. Este es el Ignorance Score (IGN). Su definición es la siguiente (en [3]): IGN(F,y) = −log(f(y)) (3.17) donde f , es la función de densidad de probabilidad correspondiente al modelo de predicción que se está evaluando. La información contenida en el logaritmo, decrece monótonamente con el incremento de la probabilidad de las observaciones, lo que es muy intuitivo. Por tanto con este método, la calidad del modelo se basa en las observaciones. Hay que destacar también que el valor de probabilidad 0 no debe ser asignado ya que su logaritmo diverge. Como en el CRPS, existen expresiones equivalentes para calcular el IGN score (en [1]) que permiten ahorrar tiempo de cálculo en el caso en el que la distribución sea normal. 18 Capítulo 3. Calibración de los modelos ignNµ,σ2,y=1 2ln(2πσ2)+ (y−µ)2 2σ2(3.18) IGN =1 n n ∑ i=1 ignNµ,σ2,yi(3.19) Por último, solamente comentar que el Ignorance Score es muy útil como medida de verosimilitud de las observaciones. Esta se puede ver haciendo el logaritmo de esas probabilidades y multiplicándolos entre sí. De esta manera, aunque no es un mal indicador, sería mejor si la estimación de los parámetros se realiza con el Método de la Máxima Verosimilitud. 3.3.4 Resultados de la aplicación de los Scores Cuando se han definido todos los indicadores que se usarán a la hora de cuantificar el grado de aproximación de un modelo con las observaciones tomadas de la realidad, es la hora de realizar ejemplos que ayuden a entender su aplicabilidad. Antes de proceder a utilizar los scores para la calibración y ajuste de los modelos, se añadirán su cálculo al Ejemplo 2. Los pasos a seguir son los siguientes (basados en [3]): 1. La realidad se va a comportar como una Normal (Normal Truth), Nµ,σ2 , con µ∼ N(25,1)yσ2=9. 2. Se generan 300 observaciones aleatorias de acuerdo a esa realidad. 3. Se realiza un ajuste de los parámetros de la Normal y de la Gumbel, mediante la estimación basada en el Método de los Momentos, usando las observaciones anteriores. 4. Se simulan otras 1000 observaciones de la realidad. 5. Mediante las observaciones simuladas por segunda vez, se calcula el CRPS, el Absolute Error y el Ignorance Score. Los cálculos se realizan según las fórmulas presentadas en la sección anterior, utilizando las funciones de probabilidad y de probabilidad acumulada de los modelos estimados y de la realidad ( F(µn,σ2 n) , F(µg,σg) , F(µ,σ2) , f(µn,σ2 n) , f(µg,σg) y f(µ,σ2) ). 6. Se repiten los pasos del 2 al 5, pero considerando ahora la realidad como una función Gumbel (Gumbel Truth), G(µ,σ), de media µ+γσ ∼N(25,1), y de varianza π2σ2 6=3π2 2. En primer lugar se muestran las curvas de cada una de las observaciones frente al valor de los indicadores. 3.3 Evaluación de modelos de predicción 19 10 15 20 25 30 35 40 0 2 4 6 8 10 12 14 Observaciones Scores Distribución Normal CRPS AE IGN 10 15 20 25 30 35 40 0 2 4 6 8 10 12 Observaciones Scores Distribución Gumbel CRPS AE IGN 10 15 20 25 30 35 40 0 2 4 6 8 10 12 Observaciones Scores Distribución Real CRPS AE IGN Figura 3.8 Cálculo de Scores para la Normal Truth siguiendo distintas distribuciones. 10 15 20 25 30 35 40 0 2 4 6 8 10 12 Observaciones Scores Distribución Normal CRPS AE IGN 10 15 20 25 30 35 40 0 2 4 6 8 10 12 Observaciones Scores Distribución Gumbel CRPS AE IGN 10 15 20 25 30 35 40 0 2 4 6 8 10 12 Observaciones Scores Distribución Real CRPS AE IGN Figura 3.9 Cálculo de Scores para la Gumbel Truth siguiendo distintas distribuciones. 20 Capítulo 3. Calibración de los modelos Se puede ver la diferente forma que presentan las curvas según el tipo de indicador. Sin embargo, no difieren mucho en cuanto a su comportamiento (tendencia, crecimiento, etc). Si que es verdad que se puede notar que en las gráficas cuya distribución del modelo empleado es similar al real (el que se ha tomado para generar las observaciones), tienen un parecido muy similar por no decir casi idéntico a la realidad (Normal-True en la figura 8 y Gumbel-True, en la figura 9). Posteriormente, se verá que eso implica que la media de los indicadores extendida a las observaciones realizadas no estará muy lejos de la distribución real o truth. Una vez vista la representación de las observaciones con respecto al valor que se obtiene de calcular los indicadores, es necesario comparar ahora el valor medio de los Scores según la distribución empleada. Verdaderamente, esta representación será la que nos aporte una medida más útil para poder comparar por ejemplo, entre una distribución Normal y una Gumbel, tendiendo detrás herramientas que justifiquen la elección de una u otra según un mejor ajuste a la realidad. La comparación se realizará siguiendo el siguiente procedimiento: -Primero se calcula la media de todas las observaciones utilizadas según lo dicho anteriormente. -Esa media se obtiene para los tres Scores: CRPS, Absolute Error e Ignorance Score. El cálculo se repite para cada una de las distribuciones y a su vez para cada una de las realidades: Normal y Gumbel, (en total 18 medias). -Para dar cierta seguridad a la media calculamos un intervalo de confianza del 95 % siguiendo la siguiente fórmula ˜x−s √ntα/2,˜x+s √ntα/2(3.20) siendo tα/2 el valor de la variable aleatoria de distribución t-student con probabilidad 1−α , igual a la probabilidad del intervalo que se quiere, en este caso en cuestión 95 % . Obviamente para usar dicho intervalo se supone la hipótesis de normalidad al ser muy elevado el número de observaciones consideradas. -Representamos tanto la media como el intervalo de confianza para cada uno de los indicadores, considerando las distintas distribuciones. Se realiza este paso dos veces: una con la Normal como realidad y otra con la Gumbel. Los resultados de las mismas se muestran a continuación: 3.3 Evaluación de modelos de predicción 21 0 0.5 1 1.5 2 2.5 3 3.5 4 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 2.05 2.1 CRPS Scores Normal Med. Normal Gumbel Med. Gumbel True Med. True 0 0.5 1 1.5 2 2.5 3 3.5 4 2 2.2 2.4 2.6 2.8 3 Absolute Error Scores Normal Med. Normal Gumbel Med. Gumbel True Med. True 0 0.5 1 1.5 2 2.5 3 3.5 4 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1 Ignorance Error Scores Normal Med. Normal Gumbel Med. Gumbel True Med. True Figura 3.10 Media de los Scores con un intervalo de confianza del 95 %. Caso Normal Truth. Tabla 3.1 Resultados de las medias de los diferentes Scores según las tres distribuciones consideradas. Caso Normal Truth. NORMAL TRUTH CRPS AE IGN Normal 1.7432 2.4508 2.5386 Gumbel 1.7612 2.5446 2.6958 True 1.7406 2.4363 2.5390 28 Capítulo 3. Calibración de los modelos entre puntos cercanos del espacio de manera que se pueda aceptar físicamente. Figura 3.15 La figura muestra los diferentes modelos de temperaturas de 4 zonas cercanas, así como la temperatura para zonas intermedias. En azul vemos un modelo obtenido suponiendo que la interpolación se ha realizado a través de una minimización del CRPS en cada punto de manera independiente. En rojo, se aprecian modelos más suaves de temperatura que tienen más sentido físico. Los modelos se han exagerado mucho para poder entender el error que existe al aplicar el método descrito en esta sección. A raíz de lo comentado, se introduce el concepto de la calibración de modelos de varias variables, que nos servirá para dar sentido al problema propuesto anteriormente. Este concepto se tratará con detenimiento en la sección posterior. 3.4.3 Calibración de modelos de varias variables Introducción Como se comentó en la sección anterior, se pueden cometer ciertos errores importantes a la hora de calibrar los modelos simplemente minimizando el CRPS en cada estación de medición de forma independiente, y luego interpolando para zonas intermedias. Para solucionar este problema rápidamente, uno puede pensar en colocar muchas estaciones de medición de variables por todas las zonas posibles y así reducir ese error de interpolación. Sin embargo, eso no es posible ya que un aumento de estaciones de medición supondría un aumento considerable de variables del problema, que repercutiría en un incremento del tiempo de cálculo que haría inabordable el problema. No sólo eso sería el único inconveniente, si no que el hecho de incrementar el número de estaciones en una región supone mucha más inversión, y por tanto, un mayor coste que depende en cierta medida de decisiones políticas, sociales, técnicas o legales en algunos casos. Lo segundo que uno puede proponerse es intentar hacer que los valores en las regiones intermedias cumplan las ecuaciones de Navier-Stokes para los fluidos o la de los campos de temperatura etc. La complejidad que aportaría meter de manera explícita toda la física del problema para reducir la incertidumbre no sólo dificultaría encontrar una solución aceptable, si no que posiblemente, esa solución llegaría más tarde del día en el que nosotros necesitamos conocer dicha predicción. Traduciendo el problema a lenguaje matemático, la covarianza entre las diferentes variables se pierde a la hora de calibrar. Es decir, si la relación entre la predicción de dos puntos es lineal, tras calibrar, deja de serla, pudiendo ser otra cualquiera. Eso es algo que la calibración de modelos de una variable por si sola no puede resolver. Dicho esto, se necesita buscar otro método basado en los anteriores que nos permita reducir esos errores, teniendo en cuenta cómo se ha obtenido el modelo en bruto original. Para ello, la calibración de varias variables puede ser una solución más que fructífera. En el subapartado siguiente se estudiará un método muy utilizado llamado Ensemble Copula Coupling (ECC). La base del planteamiento consiste en integrar todas las muestras de las diferentes estaciones en un optimizador 3.4 Post-procesado estadístico de modelos de predicción 29 cuyo resultado recoja de cierta manera por decirlo así, la covarianza de todas las zonas de medición entre sí, permitiendo obtener un modelo con cierta claridad física. Su desarrollo en profundidad se analiza a continuación. Ensemble Copula Coupling (ECC) Hablar de Ensemble Copula Coupling, es lo mismo que hablar de postprocesado de múltiples variables. Para entender la forma en la que el ECC integra todas las variables independientes en un mismo paquete, se estructurará su descripción en una serie de pasos que servirán de aclaración de cara a una posible programación futura. Se denotará ahora por el subíndice k∈ {1... K} a las distintas variables del problema. Por variable no sólo se entiende a valores de magnitudes en distintos puntos del espacio, si no que también puede referirse a magnitudes diferentes en un mismo punto, como pueden ser el viento zonal, el viento meridional, la temperatura, etc. En primer lugar, el proceso consiste en repetir los pasos del 1 al 4 del método anterior para todas las variables de forma independiente. Se consideran como variables de entrada, la predicción del modelo en bruto, Xi,j,k , las observaciones, Yi,k y los coeficientes ak,bk,ck,dk . Como salida se obtiene el CRPS medio para cada variable, ξk . Se buscan cuáles son los coeficientes a∗ k,b∗ k,c∗ k,d∗ k óptimos que minimizan cada uno de los ξk. Figura 3.16 Variables de entrada y salida en el ECC. Los siguientes pasos a dar se explican de manera detallada a continuación: 5. Se procede a calcular el Multivariate Rank para cada uno de los casos, i y para cada una de las variables, k . Lo primero que hay que hacer es ordenar de menor a mayor los valores de cada uno de los escenarios para los diferentes casos o instancias X1,(1)≤X1,(2)≤... ≤X1,(M)(3.32) . . . (3.33) Xi,(1)≤Xi,(2)≤... ≤Xi,(M)(3.34) . . . (3.35) XN,(1)≤XN,(2)≤... ≤XN,(M)(3.36) donde () denota los valores ordenados. Todo ello, se repite para las Kvariables de estudio. 6. Entre los numerosos Multivariate Rank que existen se elegirá el Multivariate Rank Average por simplicidad. Este dice que la solución es construir una tabla de ordenación, en la cual se establezca para cada uno de los casos, el nuevo orden resultante de cada escenario tras la reordenación. En una de las columnas se calcula además la media de todos los órdenes para cada uno de los casos considerados. 30 Capítulo 3. Calibración de los modelos Tabla 3.4 Tabla de reordenación de las muestras. ESCENARIO/ORDEN i=1i=2··· i=NMEDIA j=12 7 ··· ··· 4.33 j=21 5 ··· ··· 5.2 . . .. . .. . .. . .. . .. . . j=M··· ··· ··· ··· ··· 7. Se ordenan los escenarios de menor a mayor según el valor obtenido en la columna MEDIA. En caso de que dos escenarios incurran en la misma media, se supondrá que es mayor aquel con mayor índice j. El resultado se puede escribir en otra nueva tabla de la siguiente forma: Tabla 3.5 Orden global de todos los escenarios. ESCENARIO NUEVO ORDEN 1 5 2 9 . . . ? M ? 8. En concordancia con el punto 5 de la calibración de una variable, la predicción para cada una de las variables para el día o para el caso siguiente será: ˆ Xj,k . Por tanto, suponiendo que los coeficientes a∗ k,b∗ k,c∗ k,d∗ k se mantienen, se pueden construir las siguientes funciones de probabilidad: ˆ Zk∼Na∗ k+b∗ k˜ ˆ Xk,c∗ k+d∗ kˆ S2 k. para cada una de las Kvariables en cuestión. 9. De cada uno de los modelos calibrados ˆ Zk , se extraen muestras que sean equiprobables. Por ejemplo, siguiendo la siguiente fórmula. ˆ Z(j),k=ˆ Z−1 km M+1(3.37) con m∈{1....M}. Además como las muestras se han extraído de forma equiprobable, los valores obtenidos ya están ordenados de menor a mayor: ˆ Z(1),k≤ˆ Z(2),k≤... ˆ Z(M),k(3.38) para todo k∈{1....K}. 10. Finalmente, sólo hay que entrar para cada muestra, ˆ Z(j),k en la columna NUEVO ORDEN de la Tabla 3.5, e identificarla con el escenario de partida, restableciéndose el orden original. Comentar que la potencia de este método es que nos permite combinar, por una parte el hecho de que se minimiza el CRPS, ajustándose el modelo más a la realidad, y por otra parte que el reordenamiento de la muestra obtenida nos mantiene las relaciones del modelo original entre las distintas variables de cada región. Así, si por ejemplo, la temperatura o el valor del viento en Madrid siempre es superior al de Barcelona para un escenario concreto, después de calibrar, debe de continuar siéndolo para el mismo escenario, ya que si no se rompería con los principios físicos 3.4 Post-procesado estadístico de modelos de predicción 31 del problema. Con sólo la minimización, esa relación se pierde cambiando totalmente los valores de dispersión entre variables. Con este planteamiento, los mapas de vientos, temperaturas o de cualquier otra variable serán más que aceptables para poder usarlos en la práctica. Figura 3.17 Relaciones de presión y temperatura entre Berlín y Hamburgo. Ilustraciones tomadas de [7]. En la figura anterior, se puede ver un claro ejemplo de la pérdida de la relación entre variables cuando se optimiza el CRPS. El caso a) muestra las relaciones originales entre las presiones y temperatura de dos puntos pertenecientes a Berlín y Hamburgo. En el b) se observa como la minimización del CRPS provoca una dispersión de los valores de dichas variables. La solución de esto es emplear el ECC como en el caso c), lo que nos permite mantener las relaciones originales a la misma vez que se consigue un modelo de predicción más realista. 32 Capítulo 3. Calibración de los modelos Figura 3.18 Mapas de temperaturas sobre Alemania. Ilustración tomada de [8]. Otro ejemplo quizás más significativo, es el de los mapas de temperatura de una determinada región. En la Figura 3.18, se observan tres mapas de temperaturas de Alemania. Cada mapa se corresponde con el modelo en bruto, el modelo tras calibrar y el resultado tras aplicar el ECC, respectivamente. Llama claramente la atención el segundo mapa, en el cual, la calibración hace que se pierda la relación entre las temperaturas de cada zona resultando un mapa mucho menos uniforme. 4 Predicción de trayectorias El objetivo del presente capítulo es presentar al lector de una forma sencilla cómo se puede abordar la predicción de una trayectoria sujeta a determinadas restricciones. 4.1 Introducción La predicción de una trayectoria es una de las aplicaciones más inmediata, pero a la vez más importante, para la utilización de los modelos de vientos, los cuales serán usados a lo largo de este capítulo. El estudio realizado de ahora en adelante, es muy útil para dar una primera aproximación de las variaciones que introducen en las actuaciones los vientos, que se sabe que están presentes en la realidad. Mediante consideraciones como la anterior se pueden determinar cuáles son los incrementos o reducciones en el tiempo que se tarda en recorrer la ruta determinada. Este es el objetivo primordial del capítulo ya que obteniendo esos incrementos con respecto al valor nominal se tiene una idea del combustible de más o de menos para volar que se necesita, repercutiendo a su vez en un coste mayor o menor. También desde el punto de vista operativo, conocer con más precisión el tiempo que tarda una aeronave en ir de un sitio a otro permite elaborar una planificación del tráfico mucho más acotada, por lo que puede servir de vía para acelerar la gestión del tráfico aéreo. No sólo se queda la cosa aquí en lo que respecta a aplicaciones. Uno puede pensar en resolver el problema inverso y preguntarse cuál debe ser la mejor ruta para ir de una región a otra con el menor consumo posible o tardando lo mínimo posible. A su vez, los modelos de viento utilizados para la trayectoria han sido predichos minimizando indicadores como el CRPS, descrito en el capitulo anterior. La estructura de este capítulo será la siguiente: por una parte, se describirá cuáles son las hipótesis más importantes que simplifican la trayectoria, y por otra parte, se ilustrará, en base a dichas hipótesis, cuáles son las ecuaciones necesarias para su cálculo y cómo proceder si se quiere implementar el mismo. A su vez se hará hincapié durante todo el capítulo en las ventajas e inconvenientes que tienen las simplificaciones que se hacen, y cuáles son sus limitaciones al respecto. 4.2 Caracterización de la trayectoria Para facilitar el estudio de la trayectoria completa la dividiremos en un número determinado de segmentos. Cada uno de los segmentos tendrá un punto inicial y un punto final, que vendrán dados por su longitud y por su latitud (λ,φ) . Estos puntos se denominan Waypoints y servirán de apoyo para discretizar la trayectoria. Un waypoint es el resultado de la intersección de un paralelo con un meridiano terrestre. En nuestro estudio se supondrá la superficie terrestre como una cuadrícula delimitada por φmin φmax λmin λmin, con un paso de 1°en latitud y en longitud (según [9]). 33 34 Capítulo 4. Predicción de trayectorias A la hora de establecer las trayectorias que unen un punto de origen con un punto de destino, se trabajará con loxodrómicas. Esto es una ecuación de una curva terrestre del tipo φ=φ(λ) . Será muy útil poder conocer la distancia entre cualquiera dos puntos de la loxodrómica, así como el curso que hay que imponer para alcanzar el punto final. Esas dos magnitudes se pueden calcular a partir de la loxodrómica mediante desarrollos matemáticos (los cuales se omiten en este trabajo), obteniéndose: χgj=arctan     λj−λj−1 lntan(π/4−φj−1/2) tan(π/4−φj/2)      (4.1) rj        (RE+h)(φj−φj−1) cos χgj si φj6=φj−1 (RE+h)cosφj−1|λj−λj−1|si φj=φj−1 (4.2) Estas expresiones no pueden usarse tal cual, ya que hay que distinguir varios casos problemáticos (en [10]). En primer lugar, si |λj−λj−1|<πtenemos lo siguiente: -Para λj≥λj−1, si φj<φj−1hay que corregir el curso con χg=χg+π. -Para λj<λj−1 , si φj<φj−1 hay que corregir con χg=χg+π ; en caso contrario corregimos con χg=χg+2π, para tener un curso entre 0 y 360 grados. En cambio, si se da el caso de que |λj−λj−1|>π podemos recorrer la tierra cruzando el meridiano 180º(por el camino corto). Las expresiones modificadas son: -Para λj<0, la expresión del curso es χgj=arctan     λj−λj−1+2π lntan(π/4−φj−1/2) tan(π/4−φj/2)      (4.3) donde si φj<φj−1 , se corrige el curso con χg=χg+π . Cuando el curso obtenido sale negativo se le vuelve a sumar 2π, para obtener uno entre 0 y 360 grados. -Para λj≥0, la expresión del curso es χgj=arctan     λj−λj−1−2π lntan(π/4−φj−1/2) tan(π/4−φj/2)      (4.4) donde si φj<φj−1 , se corrige con χg=χg+π . Cuando el curso obtenido sale negativo, se le vuelve a sumar 2π, para obtener un curso entre 0 y 360 grados. 4.3 Ecuaciones de la trayectoria 35 En vista a las ecuaciones presentadas, se inducen varias hipótesis que se presentan a continuación: -Se asume que la Tierra es esférica de radio constante RE=6371 km. -Los segmentos que forman parte de la trayectoria serán segmentos rectos. Para simular segmentos curvos se tomarán muchos puntos en las zonas donde se quiera realizar una curva más pronunciada. Como consecuencia de lo anterior sólo se considerarán segmentos de crucero, pudiéndose hacer estos a diferentes altitudes h . Ningún segmento de subida, descenso, aterrizaje, despegue o viraje, se tendrá en cuenta, así como la transición entre distintos segmentos de crucero. Se supone que el tiempo y el consumo de combustible en dichos segmentos es mucho menor que en los segmentos de crucero, siendo esto más acertado cuanto mayor alcance tenga la misión. -Se supondrá que no se puede volar justo en los polos donde el curso a tomar siempre tiende a cero, según las expresiones anteriores. 4.3 Ecuaciones de la trayectoria En la figura se puede ver un ejemplo de trayectoria definida por un punto inicial y un punto final. Si se quiere ir de un waypoint a otro de la trayectoria, se tendrá que recorrer una cierta distancia r en línea recta, con un curso χg . La velocidad aerodinámica de la aeronave vendrá dada por V=M a , siendo M , el número de Mach al que se quiere volar y a , la velocidad del sonido a dicha altitud, a=√γaRaT. Figura 4.1 Ejemplo de segmento de trayectoria. En cada punto de la cuadrícula definida en el apartado anterior se consideran tres variables: un valor de temperatura, T , y dos valores de vientos, un viento zonal que viene del este y otro meridional en dirección norte, wuywv, respectivamente. Figura 4.2 Vientos y temperaturas para cada punto de la malla. Quizás uno de los asuntos más importantes del tema de cálculo de trayectorias es integrar los vientos en nuestro modelo de cálculo. La idea es descomponer el vector global de vientos, w , en dos componentes: una en la dirección de la trayectoria, wl, y otra de viento cruzado, wx(positivas 36 Capítulo 4. Predicción de trayectorias en la dirección de la figura). Esa descomposición es bastante simple, ya que es solamente un giro de ángulo el correspondiente al curso de la trayectoria de la aeronave. La transformación queda como Figura 4.3 Cálculo del viento longitudinal y cruzado. wl=wvcos χgj+wusin χgj(4.5) wx=wvsin χgj−wucos χgj(4.6) Por tanto, para poder volar con el curso de la trayectoria requerida, 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, dando lugar al siguiente resultado. Figura 4.4 Triángulo de vientos. Vg=qV2−w2 x+wl=q(γRaT)M2−w2 x+wl(4.7) El número de Mach se supondrá constante en cada segmento, y el resto de variables será función de cada waypoint y del tiempo: T(λ,φ,t)wl(λ,φ,t)wx(λ,φ,t) . En este proyecto, no se considerará la dependencia de estas variables con el tiempo, por lo que se supondrá que el valor de las variables estará tomado en un instante de tiempo fijo. De aquí en adelante, se omitirá la dependencia del mismo en las ecuaciones. Como en las ecuaciones de la Mecánica del Vuelo, las variables que 4.4 Resolución Numérica de la trayectoria 37 aparecen son el tiempo y la distancia, r y t , hay que escribir la longitud y la latitud en función de r : λ(r),φ(r). Para obtener estas expresiones, hay que invertir la loxodrómica de la siguiente forma φ=φj−1+r RE+hcos χgj(4.8) λ=λj−1+tan χjlntan(π/4−φj−1/2) tan(π/4−φ/2)(4.9) teniendo en cuenta que si el curso es χg=90° o χg=270° , estamos recorriendo un paralelo φj=φj−1, por lo que la ecuación de la longitud hay que sustituirla por esta otra: λ=λj−1+rsin χgj (RE+h)cosφj (4.10) Según las simplificaciones anteriores, la ecuación de la trayectoria viene dada según la Mecánica del Vuelo como: dr dt =Vg(4.11) pero como lo que nos interesa es el tiempo de vuelo, no tenemos más que invertir la ecuación anterior para obtenerlo. dt dr =1 Vg (4.12) El tiempo total de la trayectoria será la suma de los tiempos que se tarda en recorrer cada uno de los segmentos. tvuelo = N ∑ j=1 tj(4.13) 4.4 Resolución Numérica de la trayectoria El objetivo de este apartado es dar al lector unas nociones y/o instrucciones necesarias a tener en cuenta a la hora de implementar numéricamente la trayectoria. Los siguientes pasos resumen lo más destacado para realizar su cálculo numérico, una vez vistas las ecuaciones del apartado anterior. En primer lugar se considera que sólo se disponen de datos de un escenario de viento concreto. A continuación, lo que hay que realizar es: 1. Se crea una malla en φ y en λ , sabiendo que el φmin λmin φmax λmax son las coordenadas de los puntos que tengan menor latitud y longitud, y mayor latitud y longitud de la trayectoria respectivamente. El paso de esta malla será 1 ° . Si por simplicidad, se quiere trabajar con números enteros, se puede tomar 1 ° por debajo o por encima de las coordenadas anteriores. El número de puntos será el suficiente para que la malla no sea demasiado gruesa, pudiéndose modificar los límites de la misma si fuese necesario. 2. Se calcula el χgjy el rjpara cada uno de los segmentos según la loxodrómica. 3. En cada punto de la malla se conocen las coordenadas del viento zonal y del meridional, wu wv . Por tanto, como en la ecuación (4.12) la variable independiente es la distancia recorrida, se necesita saber cual es el valor del viento en puntos intermedios de la malla. Ese valor puede obtenerse mediante una interpolación lineal en 2D. Su cálculo puede realizarse si se conocen los valores de las coordenadas del punto a interpolar (expresiones 4.8 4.9 4.10). 44 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre Tabla 5.3 Tiempos de vuelos previstos para la ida de los 6 primeros meses. Primer conjunto de entrenamiento. DÍAS: ENERO-JUNIO/TIEMPOS(s) tvmedio tvmaximo tvminimo 1 ENERO 48219 48747 47775 15 ENERO 46476 46851 46198 1 FEBRERO 49653 50131 49367 15 FEBRERO 49278 49597 49002 1 MARZO 48378 48671 48038 15 MARZO 45857 46331 45301 1 ABRIL 50064 50802 49598 15 ABRIL 49999 50512 49369 1 MAYO 48067 48419 47683 15 MAYO 50413 51764 49208 1 JUNIO 48558 49004 48052 15 JUNIO 45179 45543 44790 Tabla 5.4 Tiempos de vuelos previstos para la ida de los 6 últimos meses. Primer conjunto de entrenamiento. DÍAS: JULIO-DICIEMBRE/TIEMPOS(s) tvmedio tvmaximo tvminimo 1 JULIO 47036 47283 46708 15 JULIO 45673 46284 45222 1 AGOSTO 47639 48045 47150 15 AGOSTO 45612 46573 44925 1 SEPTIEMBRE 46894 47267 46498 15 SEPTIEMBRE 47185 48060 46559 1 OCTUBRE 43307 43652 42958 15 OCTUBRE 45302 46673 43692 1 NOVIEMBRE 49970 51536 48424 15 NOVIEMBRE 47478 48754 46799 1 DICIEMBRE 47483 48259 46638 5.2 Conjuntos de calibración 45 Tabla 5.5 Tiempos de vuelos correspondientes a cada uno de los días 1 y 15 del 2017 según los vientos observados en cada punto de la trayectoria en el recorrido de ida. Primer conjunto de entrenamiento. DÍAS: ENERO-JUNIO/TIEMPOS(s) tvwobs DÍAS: JULIO-DICIEMBRE/ TIEMPOS(s) tvwobs 1 ENERO 48014 1 JULIO 47063 15 ENERO 46318 15 JULIO 45629 1 FEBRERO 49494 1 AGOSTO 47954 15 FEBRERO 49202 15 AGOSTO 45448 1 MARZO 48391 1 SEPTIEMBRE 46962 15 MARZO 45921 15 SEPTIEMBRE 46731 1 ABRIL 50186 1 OCTUBRE 43178 15 ABRIL 50171 15 OCTUBRE 45312 1 MAYO 47838 1 NOVIEMBRE 49473 15 MAYO 49916 15 NOVIEMBRE 47495 1 JUNIO 48657 1 DICIEMBRE 47529 15 JUNIO 45229 Tabla 5.6 Tiempos de vuelos previstos y observados para los primeros 14 días del mes de Diciembre. Recorrido de ida. Segundo conjunto de entrenamiento. DÍAS: DICIEMBRE/TIEMPOS(s) tvmedio tvmaximo tvminimo tvwobs 1 DICIEMBRE 47483 48259 46638 47529 2 DICIEMBRE 48333 49138 47692 48311 3 DICIEMBRE 47373 47744 46897 47387 4 DICIEMBRE 48521 49499 47824 48060 5 DICIEMBRE 47174 47653 46669 46742 6 DICIEMBRE 47652 48481 47062 46884 7 DICIEMBRE 46413 46960 45743 46348 8 DICIEMBRE 46401 46902 46059 46252 9 DICIEMBRE 46384 47072 45705 46339 10 DICIEMBRE 45430 45823 44841 45401 11 DICIEMBRE 48593 49665 47754 47954 12 DICIEMBRE 47480 48013 47055 47478 13 DICIEMBRE 47431 47826 46955 47485 14 DICIEMBRE 49733 50345 49397 49660 En vista a los resultados se puede comprobar que el máximo de todos los tiempos se obtiene para el 15 de Mayo, siendo el tiempo igual a 51764 s=14 h23 min . El tiempo más pequeño se da el día 1 de Octubre, marcando un tiempo de 42958 s=11 h56 min . Es importante remarcar la diferencia en el tiempo de vuelo cuando se tiene un cierto viento en cada punto de la trayectoria. El viento provoca incrementos y reducciones más que notables del tiempo de vuelo que afectan directamente en el consumo de combustible de la aeronave, y por tanto, en el coste del vuelo. En nuestro caso, se obtienen unos incrementos y reducciones relativas con respecto al tiempo de vuelo sin viento que pueden calcularse como: 46 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre ∆t=|tvmax −tv,w=0| tv,w=0 =0.14 ⇒14% (5.2) ∆t=|tvmin −tv,w=0| tv,w=0 =0.049 ⇒4.9% (5.3) Tabla 5.7 Tiempos de vuelos previstos para la vuelta de los 6 primeros meses. Primer conjunto de entrenamiento. DÍAS: ENERO-JUNIO/TIEMPOS(s) tvmedio tvmaximo tvminimo 1 ENERO 42658 43081 42298 15 ENERO 44140 44394 43892 1 FEBRERO 40499 40674 40189 15 FEBRERO 41505 41941 41249 1 MARZO 41757 42051 41483 15 MARZO 42101 42622 41687 1 ABRIL 41706 42158 41416 15 ABRIL 42148 42622 41783 1 MAYO 43637 44201 42933 15 MAYO 40789 41582 40051 1 JUNIO 44546 44993 44145 15 JUNIO 45242 46176 44457 Tabla 5.8 Tiempos de vuelos previstos para la vuelta de los 6 últimos meses. Primer conjunto de entrenamiento. DÍAS: JULIO-DICIEMBRE/TIEMPOS(s) tvmedio tvmaximo tvminimo 1 JULIO 44893 45437 44506 15 JULIO 46747 47778 45186 1 AGOSTO 43967 44403 43572 15 AGOSTO 46278 47031 45727 1 SEPTIEMBRE 45201 45757 44692 15 SEPTIEMBRE 43924 44232 43599 1 OCTUBRE 46667 47204 46188 15 OCTUBRE 45416 46318 44879 1 NOVIEMBRE 42553 42985 42209 15 NOVIEMBRE 41706 42333 40985 1 DICIEMBRE 43268 43729 42537 5.2 Conjuntos de calibración 47 Tabla 5.9 Tiempos de vuelos correspondientes a cada uno de los días 1 y 15 del 2017 según los vientos observados en cada punto de la trayectoria en el recorrido de vuelta. Primer conjunto de entrenamiento. DÍAS: ENERO-JUNIO/TIEMPOS(s) tvwobs DÍAS:JULIO-DICIEMBRE/ TIEMPOS(s) tvwobs 1 ENERO 42698 1 JULIO 45038 15 ENERO 44369 15 JULIO 46957 1 FEBRERO 40448 1 AGOSTO 43620 15 FEBRERO 41515 15 AGOSTO 45960 1 MARZO 41853 1 SEPTIEMBRE 45252 15 MARZO 42084 15 SEPTIEMBRE 43922 1 ABRIL 41637 1 OCTUBRE 46708 15 ABRIL 42189 15 OCTUBRE 45603 1 MAYO 43718 1 NOVIEMBRE 42491 15 MAYO 41049 15 NOVIEMBRE 41771 1 JUNIO 44580 1 DICIEMBRE 43422 15 JUNIO 45326 Tabla 5.10 Tiempos de vuelos previstos y observados para los primeros 14 días del mes de Diciembre. Recorrido de vuelta. Segundo conjunto de entrenamiento. DÍAS: DICIEMBRE/TIEMPOS(s) tvmedio tvmaximo tvminimo tvwobs 1 DICIEMBRE 43268 43729 42537 43423 2 DICIEMBRE 42681 43054 42244 42732 3 DICIEMBRE 43812 44002 43615 43842 4 DICIEMBRE 42163 42415 41900 42167 5 DICIEMBRE 41997 42260 41745 42119 6 DICIEMBRE 44690 45115 44258 45066 7 DICIEMBRE 44628 45282 44096 44608 8 DICIEMBRE 44951 45662 44528 45199 9 DICIEMBRE 45163 45774 44700 44813 10 DICIEMBRE 46483 47010 45908 46513 11 DICIEMBRE 41617 42694 40882 42148 12 DICIEMBRE 42134 42597 41697 42079 13 DICIEMBRE 42831 43129 42374 42860 14 DICIEMBRE 40585 40831 40314 40561 La importancia de implementar el recorrido de vuelta y compararlo con el ida radica en ilustrar con números como afecta la dirección del viento, y no su módulo, en el tiempo total de vuelo. El día 15 de Mayo se predice que el tiempo de vuelo para la vuelta será mínimo e igual a 40051 s=11 h8min . Si nos fijamos en el máximo de los tiempos máximos de vuelo para el recorrido de ida, que es igual a 51764 s=14 h23 min , se puede comprobar que la variación del tiempo según el sentido de la trayectoria puede ser bastante amplia por el hecho del cambio de dirección de los vientos. En términos de porcentajes relativos el tiempo de vuelo puede llegar a ser un 30% mayor en la ida que en la vuelta. 48 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre Para ver más claramente cómo funcionan las predicciones para cada uno de los días, se puede representar en una serie de figuras el valor del tiempo medio, máximo y mínimo de vuelo, así como el tiempo de vuelo correspondiente a las observaciones reales, para cada uno de los conjuntos de entrenamiento. De esta manera, se puede ver gráficamente, donde cae el tiempo real, es decir, si está dentro del rango previsto, muy alejado o no de la media, o incluso fuera de rango (en consonancia con la figura 3.12). 4.75 4.8 4.85 4.9 x 1044.6 4.65 4.7 x 1044.92 4.94 4.96 4.98 5 5.02 x 1044.9 4.92 4.94 4.96 x 104 4.8 4.82 4.84 4.86 4.88 x 1044.5 4.55 4.6 4.65 x 1044.95 5 5.05 5.1 x 1044.9 4.95 5 5.05 5.1 x 104 4.76 4.78 4.8 4.82 4.84 4.86 x 1044.9 5 5.1 5.2 x 1044.8 4.85 4.9 4.95 x 1044.46 4.48 4.5 4.52 4.54 4.56 x 104 4.66 4.68 4.7 4.72 4.74 x 1044.5 4.55 4.6 4.65 x 1044.7 4.75 4.8 4.85 x 1044.45 4.5 4.55 4.6 4.65 4.7 x 104 4.64 4.66 4.68 4.7 4.72 4.74 x 1044.65 4.7 4.75 4.8 4.85 x 1044.28 4.3 4.32 4.34 4.36 4.38 x 1044.3 4.4 4.5 4.6 4.7 x 104 4.8 4.9 5 5.1 5.2 x 1044.65 4.7 4.75 4.8 4.85 4.9 x 1044.65 4.7 4.75 4.8 4.85 x 104 Figura 5.2 Representación de los tiempos de vuelo durante la ida para el primer conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales. 4.65 4.7 4.75 4.8 4.85 x 1044.75 4.8 4.85 4.9 4.95 x 1044.68 4.7 4.72 4.74 4.76 4.78 x 1044.75 4.8 4.85 4.9 4.95 x 104 4.65 4.7 4.75 4.8 x 1044.65 4.7 4.75 4.8 4.85 x 1044.55 4.6 4.65 4.7 x 1044.6 4.65 4.7 x 104 4.55 4.6 4.65 4.7 4.75 x 1044.45 4.5 4.55 4.6 x 1044.75 4.8 4.85 4.9 4.95 5 x 1044.7 4.75 4.8 4.85 x 104 4.65 4.7 4.75 4.8 x 1044.9 4.95 5 5.05 x 104 Figura 5.3 Representación de los tiempos de vuelo durante la ida para el segundo conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales. 5.2 Conjuntos de calibración 49 4.22 4.24 4.26 4.28 4.3 4.32 x 1044.38 4.4 4.42 4.44 x 1044 4.02 4.04 4.06 4.08 x 1044.12 4.14 4.16 4.18 4.2 x 104 4.14 4.16 4.18 4.2 4.22 x 1044.15 4.2 4.25 4.3 x 1044.14 4.16 4.18 4.2 4.22 x 1044.15 4.2 4.25 4.3 x 104 4.25 4.3 4.35 4.4 4.45 x 1044 4.05 4.1 4.15 4.2 x 1044.4 4.45 4.5 x 1044.4 4.45 4.5 4.55 4.6 4.65 x 104 4.45 4.5 4.55 x 1044.5 4.6 4.7 4.8 x 1044.35 4.4 4.45 x 1044.55 4.6 4.65 4.7 4.75 x 104 4.45 4.5 4.55 4.6 x 1044.34 4.36 4.38 4.4 4.42 4.44 x 1044.6 4.65 4.7 4.75 x 1044.45 4.5 4.55 4.6 4.65 x 104 4.22 4.24 4.26 4.28 4.3 x 1044.05 4.1 4.15 4.2 4.25 x 1044.25 4.3 4.35 4.4 x 104 Figura 5.4 Representación de los tiempos de vuelo durante la vuelta para el primer conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales. 4.25 4.3 4.35 4.4 x 1044.22 4.24 4.26 4.28 4.3 4.32 x 1044.36 4.37 4.38 4.39 4.4 4.41 x 1044.18 4.2 4.22 4.24 4.26 x 104 4.16 4.18 4.2 4.22 4.24 x 1044.42 4.44 4.46 4.48 4.5 4.52 x 1044.4 4.45 4.5 4.55 x 1044.45 4.5 4.55 4.6 x 104 4.45 4.5 4.55 4.6 x 1044.55 4.6 4.65 4.7 4.75 x 1044.05 4.1 4.15 4.2 4.25 4.3 x 1044.16 4.18 4.2 4.22 4.24 4.26 x 104 4.22 4.24 4.26 4.28 4.3 4.32 x 1044.02 4.04 4.06 4.08 4.1 x 104 Figura 5.5 Representación de los tiempos de vuelo durante la vuelta para el segundo conjunto de entrenamiento. En azul, los tiempos máximos y mínimos previstos para cada día; en rojo, la media de los tiempos de cada escenario; y en verde, el tiempo real basado en observaciones reales. 50 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre Con la representación de las figuras anteriores se puede ver claramente, las diferentes posiciones de las observaciones en función del día correspondiente. En general, se observa que el valor medio del tiempo de vuelo de las predicciones es bastante cercano al valor observado. En la gran mayoría de los casos la observación cae dentro del intervalo formado por los tiempos máximos y mínimos de los 50 escenarios, lo que induce a pensar que con bastante probabilidad ocurrirá lo mismo para el día 15 de Diciembre, del cual se quiere conocer el tiempo de vuelo. Aunque el histórico de predicciones sea muy acertado, no quiere decir que el día 15 de Diciembre la observación tenga que caer dentro del intervalo. Es más, si se mira con detenimiento, uno puede darse cuenta que para el recorrido de ida del día 6 de Diciembre ( 2 ª fila y 2 ª columna de la figura 5.3), el tiempo correspondiente a la observación es menor que el mínimo tiempo previsto por los 50 escenarios ese día. Si ocurriera esto para nuestro día de estudio, la calibración no sólo sería útil para ganar precisión, si no que sería imprescindible. 5.3 Resultados de los vientos tras la calibración Antes de mostrar los resultados obtenidos es necesario decir que la optimización del CRPS se ha realizado utilizando el comando fmincon de Matlab. Las condiciones iniciales y las restricciones ya fueron aclaradas anteriormente (véase sección 3.4). Falta por mencionar cuál ha sido la tolerancia elegida para la optimización de las variables. Se ha escogido un valor de tolerancia de 10−5 , ya que se ha comprobado que eligiendo un valor mayor no se consigue mucha más precisión, y sí un aumento en el tiempo de cálculo. Dicho esto, se procede a estudiar los resultados obtenidos cuando se utiliza el primer conjun- to de entrenamiento. Como se tienen muchos puntos en la trayectoria, se analizan sólo dos puntos cercanos correspondientes a las ciudades de Berlín y Hamburgo. En cada uno de los puntos se tienen dos variables: viento zonal y viento meridional. Sólo se mostrarán los resultados para el viento zonal y no para el viento meridional al ser análogos. Con ello, el cálculo del CRPS para el modelo en bruto y para el modelo calibrado puede verse a continuación. Tabla 5.11 Resultados de la optimización para el primer conjunto de entrenamiento. wuξ ξoptimo a b c d Berlín 0.6309 0.6116 −0.116 1.0123 0.3416 0.3280 Hamburgo 0.7481 0.7193 −0.1428 1.0288 0 0.9362 De los resultados pueden extraerse varias conclusiones. Por un lado, un valor tan bajo de CRPS induce a pensar que los modelos proporcionados de predicción ya son bastante buenos de por sí, por lo que una optimización sólo ayuda a mejorarlos más bajando su valor, pero no en gran medida. Llama la atención también, que en los dos casos el valor de a es negativo, siendo necesario un sesgo negativo para corregir el modelo. En el caso de Hamburgo la optimización del CRPS se consigue haciendo cero uno de los parámetros, como es, el coeficiente c. 5.3 Resultados de los vientos tras la calibración 51 17 17.5 18 18.5 19 19.5 20 20.5 21 21.5 22 22.5 23 23.5 24 Modelo en bruto Hamburgo Berlín 17 17.5 18 18.5 19 19.5 20 21.5 22 22.5 23 23.5 24 24.5 25 EMOS Modelo Berlín Hamburgo 17 17.5 18 18.5 19 19.5 20 20.5 21.5 22 22.5 23 23.5 24 24.5 ECC Modelo Hamburgo Berlín Figura 5.6 Representación del ensemble forecast antes de la calibración, tras realizar la calibración individual, EMOS y una vez aplicado el ECC . Se puede ver que optimizando a través del EMOS se pierde la tendencia aproximadamente lineal que tenía la nube de puntos del modelo en bruto. Al aplicar el ECC a este modelo se recupera gran parte de esa tendencia motivada por el uso de la reordenación del modelo en bruto usando el Multivariate Rank Average, cómo se explicó en la sección 3.4. Tabla 5.12 Coeficientes de correlación entre ambas ciudades para el primer conjunto de entrenamiento. BRUTO EMOS ECC COEFICIENTES DE CORRELACIÓN 0.7341 −0.0254 0.7619 El coeficiente de correlación del EMOS es prácticamente nulo mientras que al aplicar el ECC la relación entre las variables se recupera, obteniéndose un valor de correlación similar al original. Si se usa ahora el segundo conjunto de entrenamiento (primeros 14 días del mes de Diciembre) los valores del CRPS y de los coeficientes que se obtienen son: Tabla 5.13 Resultados de la optimización para el segundo conjunto de entrenamiento. wuξ ξoptimo a b c d Berlín 0.4814 0.4220 0.1235 1.005 0.5071 0.0957 Hamburgo 0.9807 0.8850 −0.7558 1.0486 0.3759 0.8108 Lo que más sorprende si se comparan estos valores con el primer conjunto de entrenamiento es que el valor del CRPS óptimo para Berlín es menor en el segundo caso que en el primero (0.6116 52 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre vs 0.4220). Sin embargo, para Hamburgo ocurre justamente lo contrario. El CRPS es menor en el primer caso y no en el segundo (0.7193 vs 0.8850). Aun así, no se pueden sacar conclusiones generales en base a estos resultados ya que sólo se muestran dos variables. En la sección siguiente, se calcularán los tiempos de vuelos, que es el objetivo de nuestra aplicación. A partir de ahí, se podrá decidir y valorar sobre cuál modelo merece más la pena usar para ajustar, o si por lo contrario, ninguno de los dos aporta demasiado, visto los valores tan pequeños del CRPS. Comentar también, que al igual que en el primer caso, se obtiene un valor de a negativo para Hamburgo, lo que implica un sesgo negativo al ajustar. Los valores de c y d en general son menores de 1 mientras que los valores de bsuelen alcanzar valores mayores de 1. 17 17.5 18 18.5 19 19.5 20 20.5 21 21.5 22 22.5 23 23.5 24 Modelo en bruto Hamburgo Berlín 16.5 17 17.5 18 18.5 19 19.5 20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 EMOS Modelo Hamburgo Berlín 17 17.5 18 18.5 19 19.5 20 20.5 21 21.5 22 22.5 23 23.5 24 24.5 25 ECC Modelo Hamburgo Berlín Figura 5.7 Representación del ensemble forecast antes de la calibración, tras realizar la calibración individual, EMOS y una vez aplicado el ECC . Ocurre lo mismo que para el primer conjunto de entrenamiento. El ECC permite recuperar algo mejor la relación primitiva existente en el modelo en bruto entre las ciudades de Berlín y Hamburgo, la cual se pierde al calibrar con el EMOS de forma independiente. Tabla 5.14 Coeficientes de correlación entre ambas ciudades para el segundo conjunto de entrenamiento. BRUTO EMOS ECC COEFICIENTES DE CORRELACIÓN 0.7341 −0.0779 0.6309 Al igual que para el primer conjunto de entrenamiento, la tendencia se recupera bastante. Además, se puede destacar que, mientras en el primer conjunto la tendencia resultante era algo más lineal que en el modelo en bruto (0.7341 vs 0.7619), en este caso ocurre lo contrario, siendo ésta algo menos lineal que la original (0.7341 vs 0.6309). En cualquier caso, es evidente la ganancia obtenida mediante el ECC. 5.3 Resultados de los vientos tras la calibración 53 Para ilustrar lo que ocurre al calibrar no sólo en un punto si no en toda la malla de estudio, se presentan en una serie de figuras los valores del viento zonal y meridional para un escenario cualquiera del modelo en bruto, del modelo calibrado mediante el EMOS de forma individual y tras el ECC. −80 −60 −40 −20 0 25 30 35 40 45 50 Longitud Latitud Viento zonal bruto −10 0 10 20 30 40 50 60 70 −80 −60 −40 −20 0 25 30 35 40 45 50 Longitud Latitud Viento zonal EMOS −10 0 10 20 30 40 50 60 70 80 90 −80 −60 −40 −20 0 25 30 35 40 45 50 Longitud Latitud Viento zonal ECC −10 0 10 20 30 40 50 60 70 Figura 5.8 Representación de las curvas de nivel para el viento zonal. Primer conjunto de entrenamiento. 60 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre siguientes: tvmedio =41293 s tvmaximo =41588 s tvminimo =40928 s(5.26) Desv.Tipica S =152.71 s(5.27) tvwobs =41288 s(5.28) |tvmedio −tvwobs |=5s(5.29) |tvmaximo −tvminimo |=660 s(5.30) |tvmaximo −tvwobs |=300 s(5.31) |tvwobs −tvminimo |=360 s(5.32) La mejora relativa que se tiene en el recorrida de vuelta puede calcularse como %Me jora =34 −5 34 ·100 =83% (5.33) la cual es mucho más elevada que para el recorrido de ida. Es curioso observar que a pesar de que el tiempo medio se acerca muchísimo al real, la dispersión de los 50 escenarios aumenta levemente (151.85 s a 152.71 s), haciendo posible que ese acercamiento sea mucho más llamativo que en la ida. Con ello, los tiempos máximos y mínimos se expanden, por así decirlo, aumentando y disminuyendo respectivamente. Gráficamente, se puede ver en la figura siguiente. 4.09 4.1 4.11 4.12 4.13 4.14 4.15 4.16 4.17 x 104 Modelo en bruto vs Modelo Calibrado Tiempos de vuelo (s) Figura 5.15 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de vuelta. Primer conjunto de entrenamiento. Una vez terminado el desarrollo de los tiempos de vuelos provenientes del primer conjunto de entrenamiento se prosigue con los resultados obtenidos utilizando el Ensemble Copula Coupling con el segundo conjunto de entrenamiento. A modo de recordatorio, este segundo conjunto, está formado por las previsiones de vientos para los primeros 14 días del mes de Diciembre junto con los valores reales para esos días. El modelo en bruto a calibrar será el mismo: el correspondiente al día 15 de Diciembre, cuyas características ya se detallaron anteriormente. Para el recorrido de ida el modelo calibrado tiene la siguientes forma y características: tvmedio =48398 s tvmaximo =48856 s tvminimo =47977 s(5.34) Desv.Tipica S =251.24 s(5.35) 5.4 Tiempos de vuelo tras el postprocesado 61 tvwobs =48565 s(5.36) |tvmedio −tvwobs |=167 s(5.37) |tvmaximo −tvminimo |=879 s(5.38) |tvmaximo −tvwobs |=291 s(5.39) |tvwobs −tvminimo |=588 s(5.40) El margen de mejora de este modelo con respecto al modelo en bruto puede calcularse como %Me jora =185 −167 185 ·100 =9.7%.(5.41) Puede apreciarse que es bastante menor que el obtenido con respecto al primer conjunto de entrenamiento 26 % . También lo es la dispersión, S , siendo ésta algo menor para el primer conjunto: 224.84 s en lugar de 251.24 s . Sin embargo, el rango es mayor en el primer caso que en el segundo (903 s vs 879 s). Comparando el modelo calibrado con el modelo en bruto el resultado es lo siguiente. 4.76 4.78 4.8 4.82 4.84 4.86 4.88 4.9 4.92 x 104 Modelo en bruto vs Modelo Calibrado Tiempos de vuelo (s) Figura 5.16 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de ida. Segundo conjunto de entrenamiento. En cuanto a la vuelta, el modelo calibrado con sus características se muestra a continuación: tvmedio =41265 s tvmaximo =41570 s tvminimo =40905 s(5.42) Desv.Tipica S =159.86 s(5.43) tvwobs =41288 s(5.44) |tvmedio −tvwobs |=23 s(5.45) |tvmaximo −tvminimo |=665 s(5.46) |tvmaximo −tvwobs |=282 s(5.47) |tvwobs −tvminimo |=383 s(5.48) Ahora el nuevo margen de mejora del tiempo medio con respecto al valor observado es de %Me jora =34 −23 34 ·100 =32.4%.(5.49) 62 Capítulo 5. Aplicación. Trayectoria sujeta a incertidumbre Comparado con el obtenido en la vuelta para el primer conjunto de entrenamiento, 83 % este nuevo margen es muy inferior, lo que nos hace notar que quizás el primer conjunto sea más adecuado para realizar la calibración. Es más, también la dispersión de los 50 escenarios para el primer conjunto es menor que para el segundo (152.71 s vs 159.86). La diferencia de tiempos es ahora algo mayor pero en menor medida, pasando de 660 s a 665 s. Por último, superponiendo el modelo en bruto junto con el modelo calibrado se llega a 4.09 4.1 4.11 4.12 4.13 4.14 4.15 4.16 4.17 x 104 Modelo en bruto vs Modelo Calibrado Tiempos de vuelo (s) Figura 5.17 Modelo en bruto y calibrado de los tiempos de vuelo para el recorrido de vuelta. Segundo conjunto de entrenamiento. 5.5 Conclusiones generales Con vista a los resultados anteriores se incluye una sección en la que se habla de las conclusiones generales que pueden sacarse tras aplicarse el postprocesado a los modelos en bruto de partida. Se analizan las mismas en base a lo obtenido con su aplicación a la predicción de trayectorias teniendo en cuenta diversos puntos que puedan tener cierta influencia en un futuro. A continuación, se resumen algunas de ellas: -En primer lugar, una vez terminado el desarrollo del Ensemble Copula Coupling es importante discutir de la rentabilidad de su aplicación a la predicción de trayectorias a la hora de ganar precisión en el ajuste del tiempo de vuelo. Si que es verdad que aunque el ECC proporciona bastante precisión, también supone un gasto desde el punto de vista numérico. Lo que hay que mirar en este caso es si la predicción de trayectorias requiere de gran precisión a la hora de su operatividad en la práctica. No cabe olvidar que son muchas las simplificaciones que se asumen en su cálculo. Algunas de ellas son la eliminación de subidas, descensos, aterrizajes y la velocidad constante suponiendo tramos rectos. Por tanto, los valores obtenidos pueden variar en el torno de minutos o incluso de horas si se meten las esperas para aterrizar, o las emergencias. La precisión que nos da el ECC es del orden de varios minutos, por lo que en teoría no sería muy necesario hacerla porque los modelos proporcionados por el ECMWF se ajustan de por si bastante bien a la realidad. -Con respecto a la necesidad de aplicar el ECC comentada en el párrafo anterior hay que añadir ciertas reseñas. Con los datos de partida proporcionados en este proyecto no se puede establecer una valoración más profunda sobre la rentabilidad de la aplicación del método. La justificación es que la media de los 50 escenarios para el día 15 de Diciembre cae muy cerca del valor real, y en ningún caso, fuera del intervalo de los valores máximos y mínimos, que es donde se puede 5.5 Conclusiones generales 63 notar más la utilidad del método. También, el ECC se ha aplicado al viento zonal y al viento meridional y no a los tiempos de vuelo. Los tiempos de vuelo se obtienen de una forma indirecta donde entran en juego otras variables como la velocidad aerodinámica que en muchos casos es de un orden superior a los vientos, por lo que un cambio de precisión en los mismos, no tiene una gran repercusión en el tiempo de vuelo. Quizás en otras aplicaciones de la estadística, meteorología o incluso de la aeronáutica donde los vientos repercutan de manera más directa o dónde su precisión afecte de manera determinante para el cálculo de una magnitud, sea mucho más útil este método. -Por último es muy importante destacar las diferencias existentes al aplicar el método con conjuntos de entrenamiento diferentes. En el primer caso el número de casos, instancias o días N es mayor N=23 , pero están más separados en el tiempo al ser previsiones sobre los días 1 y 15 de cada mes. En el segundo caso, N=14 , pero las previsiones están más juntas en tiempo por ser del mes de Diciembre. Es curioso que visto lo visto con los resultados de la sección anterior en el que el primer conjunto de entrenamiento daba valores mucho más precisos que el segundo, sea en el cual las previsiones están realizadas durante días más espaciados en el tiempo. Nuestro sentido común diría lo contrario. Aunque el número de casos sea más pequeño, al estar éstos más cercanos en tiempo sería más fácil captar la tendencia de las observaciones y por tanto realizar un mejor ajuste del modelo. El caso es que en nuestro ejemplo, un número mayor de N , nos aporta más precisión captando mejor el valor de las observaciones aunque las predicciones estén más equiespaciadas en el tiempo. 6 Futuras Mejoras. Recomendaciones En este capítulo final se abordan líneas que quedan abiertas tras este proyecto, y en las cuales se puede influir en un futuro. Cabe destacar algunas muy interesantes como pueden ser las siguientes: -Probar a realizar la calibración optimizando otros tipos de Scores, por ejemplo, el Ignorance Score. Para ello, se pueden usar otros ajustes como el método de Máxima Verosimilitud, y ver si existen diferencias significativas. Por lo que se vio en el capítulo 3 no es esperable encontrar grandes diferencias siempre y cuando la hipótesis de normalidad de los vientos se siga cumpliendo, cosa que se comprueba que es así en la realidad. -Realizar la optimización de varias variables de manera más realista teniendo en cuenta la restricción real de c+dS2≥0 , para comprobar que no hay ningún problema al igual que señala la literatura. -Cambiar las variables de optimización o añadir otras variables nuevas como pueden ser la presión y la temperatura. Incluso al cambiar las variables, puede probarse a modificar las distribuciones normales por otras que recojan mejor la tendencia real (véase tabla 3.3). -Variar mucho más el número de casos, N , incluso en un horizonte temporal mucho más extenso, para ver si existen grandes diferencias con valores más pequeños del parámetro N. -Incorporar el escenario de alta resolución a la calibración y comprobar si es verdad que éste aporta mucha más precisión que el resto de escenarios. -Emplear otros tipos de Multivariate Rank que no sea el Multivariate Rank Average. El hecho de utilizar otro Rank cambiaría el criterio de reordenación del modelo calibrado, y por tanto, variaría el resultado del ECC. -Utilizar otros métodos de calibración en lugar del Ensemble Copula Coupling. Uno de ellos es conocido en la literatura como el Ensemble Member a Member, contenido en la bibliografía de este proyecto (véase [4]). -En cuanto a la trayectoria, puede probarse a variar el paso de la malla, para tener una idea de la variación que pueda suponer en los tiempos de vuelo cualquier refinamiento de la misma. -Modificar la ruta de destino y origen para probar qué ocurre con otras combinaciones. -Incorporar otros tipos de segmentos a la trayectoria: rectos y curvos con subidas, descensos y 65 66 Capítulo 6. Futuras Mejoras. Recomendaciones aceleraciones, teniendo en cuenta sus respectivas ecuaciones y la presencia del viento en las mismas. -Hacer que en la trayectoria se tengan en cuenta la variación de los parámetros con el tiempo, para distintos horizontes temporales, y comparar los resultados de unos con otros. -Mezclar la optimización de variables realizadas con el ECC junto con la optimización de trayectorias, para encontrar una trayectoria que sea óptima, desde un punto de vista quizás algo más realista. -Por último, se pueden emplear los métodos descritos en este proyecto, para numerosas aplicaciones en la que intervengan los vientos o cualquier variable atmosférica, y que además, se requiera de gran precisión. Todas estas posibles modificaciones hacen que el proyecto sea bastante abierto y con una buena perspectiva de futuro, siempre siendo consciente de la dificultad que entraña realizar algunas de las propuestas. 7 Anexos Se recogen aquí algunos códigos numéricos en Matlab que pueden ser útiles o interesantes. A. Anexo PIT %CÁLCULO DE HISTOGRAMAS PARA UNA NORMAL Y PARA UNA GUMBEL ASUMIENDO QUE LA %REALIDAD SE COMPORTA COMO UNA NORMAL clear all; close all; clc; %GENERAMOS 300 OBSERVACIONES mu_mu=25; var_mu=1; var_Y=9; mu_Y=mu_mu+sqrt(var_mu).*randn(300,1); Y=mu_Y+sqrt(var_Y).*randn(300,1); %AJUSTAMOS LA NORMAL Y LA GUMBEL mu_n=mean(Y) var_n=var(Y,1) %============================================== var_g=(6/pi^2)*var_n; mu_g=mu_n-double(eulergamma)*sqrt(var_g); %TRANSFORMAMOS MEDIANTE 1000 Y 1E6 SIMULACIONES CON LA CDF DEL MODELO %AJUSTADO N=1000000; mu_Y=mu_mu+sqrt(var_mu).*randn(N,1); Y=mu_Y+sqrt(var_Y).*randn(N,1); % Y=mu_n+sqrt(var_n).*randn(N,1); z_n=normcdf(Y,mu_n,sqrt(var_n)); z_g=1-evcdf(-Y,-mu_g,sqrt(var_g)); figure(1) hist(z_n,50) title(’\bf{Normal Nob=1e6}’) figure(2) hist(z_g,50) title(’\bf{Gumbel Nob=1e6}’) 67 68 Capítulo 7. Anexos %CÁLCULO DE HISTOGRAMAS PARA UNA NORMAL Y PARA UNA GUMBEL ASUMIENDO QUE LA %REALIDAD SE COMPORTA COMO UNA GUMBEL clear all; close all; clc; %GENERAMOS 300 OBSERVACIONES var_Y=9; mu_Y=(25-double(eulergamma)*sqrt(var_Y))+sqrt(1).*randn(300,1); Y=-evrnd(-mu_Y,sqrt(var_Y)); %AJUSTAMOS LA NORMAL Y LA GUMBEL mu_n=mean(Y); var_n=var(Y,1); %============================================== var_g=(6/pi^2)*var_n mu_g=mu_n-(double(eulergamma)*sqrt(var_g)) %TRANSFORMAMOS MEDIANTE 1000 Y 1E6 SIMULACIONES CON LA CDF DEL MODELO %AJUSTADO N=1000000; mu_Y=(25-double(eulergamma)*sqrt(var_Y))+sqrt(1).*randn(N,1); Y=-evrnd(-mu_Y,sqrt(var_Y)); % Y=-evrnd(-mu_g,sqrt(var_g),N,1); z_n=normcdf(Y,mu_n,sqrt(var_n)); z_g=1-evcdf(-Y,-mu_g,sqrt(var_g)); % z_g=1-evcdf(-Y,-(25-double(eulergamma)*3),3); figure(1) hist(z_n,50) title(’\bf{Normal Nob=1e6}’) figure(2) hist(z_g,50) title(’\bf{Gumbel Nob=1e6}’) 7.1 B. Cálculo de Scores function [m,CRPS]=calcula_crps_normal(Yob,par_n) Yob=sort(Yob); N=length(Yob); %CALCULO SU CRPS PARA EL CASO DE LA NORMAL for j=1:N CRPS(j)=quadgk(@(tau)myint_normal(tau,Yob(j),par_n),0,1,’AbsTol’,1e-2,’ RelTol’,1e-2); % CRPS1(j)=quadgk(@(tau)myint_normal(tau,Yob(j),par_n),0,1,’AbsTol’,1e -4,’RelTol’,1e-4); end % err=abs((norm(CRPS-CRPS1))/(norm(CRPS))) m=mean(CRPS); % plot(Yob,CRPS) end function int=myint_normal(tau,Yob,par_n) 7.2 C. ECC 69 int=2.*(norminv(tau,par_n(1),sqrt(par_n(2)))-Yob).*((1+sign(norminv(tau, par_n(1),sqrt(par_n(2)))-Yob)).*0.5-tau); end function [m,IGN]=calcula_ignorance_error_normal(Yob,par_n) Yob=sort(Yob); N=length(Yob); %CALCULAMOS EL IGNORANCE SCORE PARA LA NORMAL for i=1:N IGN(i)=-log(normpdf(Yob(i),par_n(1),sqrt(par_n(2)))); end m=mean(IGN); % plot(Yob,IGN) end function [m,AE]=calcula_absolute_error_normal(Yob,mediana_n) Yob=sort(Yob); N=length(Yob); %CALCULAR EL ABSOLUTE ERROR ASOCIADO A UNA NORMAL for i=1:N AE(i)=abs(mediana_n-Yob(i)); end m=mean(AE); % plot(Yob,AE) end 7.2 C. ECC % CALIBRACIÓN DE LOS MODELOS clear all; clc; close all; % ESTABLECEMOS EL NÚMERO DE VARIABLES, NUMERO DE CASOS Y NUMERO DE ESCENARIOS K=5568; %número de variables N=23;%número de casos M=50; %número de escenarios %1o ¯Cargo las observaciones y 2o ¯Me genero el conjunto de training y la muestra a calibrar load genera_entrada.mat % Xc(1,j)=(20+j*0.2)+3.*randn; % Xc(2,j)=(17+j*0.2)+1.*randn; %3o ¯Me cálculo la media de cada una de la variables y la varianza muestral Xb=mean(X,2); Sc=var(X,1,2); 76 Referencias [13] Guidelines on Ensemble Prediction Systems and Forecasting. World Meteorological Organization (WMO), 2012. [14] The ECMWF Ensemble Prediction System. European Centre for Medium-Range Weather Forecasts (ECMWF), 2012. [15] R. Mladek and F. Vitart. (2016, Sep.) Models. [Online]. Available: https://software.ecmwf.int/wiki/ display/TIGGE/Models. [16] Mercedes Cordero Romero. Optimización de trayectorias sujetas a incertidumbre meteorológica. Trabajo Fin de Máster. Sevilla, 2016.