4. Predicción de tiempos de espera
Como nuestros datos cuentan con más de un mensaje para un mismo vuelo y siguen un patrón de serie temporal, no podemos realizar un random split de los datos. Así pues, hemos decidido separar en entrenamiento y test del siguiente modo:
- Train: 07/11/2024 a 14/01/2025
- Test: 15/01/2025 al 31/01/2025
Esta división evita problemas de contaminación del futuro y escinde los datos en aproximadamente un 80/20%, respectivamente. Por otro lado, en el entrenamiento de los modelos se podrán emplear distintos métodos de validación. Sin embargo, estos métodos de validación también deberán tener en cuenta las series temporales.
Asimismo, para la evaluación de los modelos de forma individual, hemos implementado un dashboard con Plotly que ofrece un informe detallado de las métricas y distribución de los errores. En este se muestra la siguiente información:
- Métricas globales: MAE, RMSE, R² y MAPE en test, con diferencia con respecto a la validación.
- Diagrama de barras con MAE y RMSE por pista / punto de espera
- Mapa de calor del MAE / RMSE por día / hora
- Distribución del error global (histograma y boxplot)
- Diagrama de dispersión del error en función de la variable numérica seleccionada
- Diagrama de dispersión predicciones vs. valores reales
- Boxplots con el error por zonas (pistas / puntos de espera)
- Mapa de calor geográfico con el MAE / RMSE por pista / punto de espera
- Distribución del error por categoría de turbulencia (avión actual y anterior)
- Evolución del MAE y RMSE de forma diaria
Para llevar a cabo un seguimiento de todos los experimentos realizados y poder comparar los resultados, hemos empleado la plataforma MLflow. En concreto, cada vez que hemos realizado un experimento se almacena la siguiente información:
- Datos generales: nombre del modelo, framework, variable respuesta, preprocesado de los datos, conjunto de datos, semilla, estrategia de búsqueda de hiperparámetros
- Hiperparámetros óptimos
- Métricas: tiempo de entrenamiento; MAE y RMSE globales para train, validación y test; MAE y RMSE para test por pistas y puntos de espera; R² y MAPE para test.
- Artefacto del modelo
A continuación explicaremos todos los modelos que hemos entrenado. Los primeros modelos son más simples y sirven como baseline para los siguientes. Evaluaremos todos los modelos en el siguiente apartado.
Esta primera heurística parte de un tiempo base (100s) que se le ajusta sumando (por categoría de turbulencia del avión anterior, tráfico del aeropuerto, tiempo desde el último despegue, horas pico) o restando (días festivos) segundos en función de nuestro conocimiento operativo extraído de la exploración de las variables. Se trata de una heurística rápida y simple que no requiere preprocesamiento adicional de datos.
Esta segunda heurística se define como:
«Un avión esperará el promedio de los últimos tres tiempos de espera de las aeronaves en su misma pista».
A la hora de calcular las predicciones es importante coger aviones con cuyo tiempo de espera no se exceda el timestamp del avión que nos interesa, puesto que estaríamos empleando información futura. De nuevo, no hay preprocesamiento, pero las operaciones de agregación hacen que sea una heurística significativamente lenta.
Para construir este modelo hemos utilizado PySpark ya que la clase LinearRegressor está disponible dentro de la MLlib de Spark. Este modelo servirá como baseline para los modelos de mayor complejidad.
En primer lugar, realizamos el preprocesamiento mediante un pipeline que incluye principalmente la indexación de variables categóricas y normalización de variables numéricas.
Para encontrar la mejor combinación de hiperparámetros empleamos validación cruzada extendida o expanding window, que consiste en entrenar el modelo con un período creciente de datos y validar en el siguiente punto en el tiempo. Para ello dividimos el conjunto de entrenamiento en 5 folds respetando el orden temporal, y después entrenamos y evaluamos con un bucle, almacenando las distintas métricas. Realizamos una validación cruzada "manual" debido a la naturaleza de los datos, asegurándonos de que no se van a producir filtraciones temporales.
| Fold | Train Period | Validation Period |
|---|---|---|
| 1 | 07/11/2024 - 30/11/2024 | 01/12/2024 - 07/12/2024 |
| 2 | 07/11/2024 - 07/12/2024 | 08/12/2024 - 14/12/2024 |
| 3 | 07/11/2024 - 14/12/2024 | 15/12/2024 - 21/12/2024 |
| 4 | 07/11/2024 - 21/12/2024 | 22/12/2024 - 28/12/2024 |
| 5 | 07/11/2024 - 28/12/2024 | 29/12/2024 - 04/01/2025 |
Una vez obtenidos los mejores hiperparámetros, entrenamos el modelo final con todos los datos de entrenamiento disponibles.
Analizando la importancia de las variables, vemos que las de mayor influencia corresponden a la variable categórica que indica la aerolínea, y también dos de las variables que indican si un punto de espera específico está o no ocupado.
Para este modelo hemos usado la librería de scikit-learn, que proporciona una implementación robusta de RandomForestRegressor.
El preprocesamiento se realiza mediante un pipeline, que incluye la transformación de variables categóricas con OneHotEncoder y el manejo de variables numéricas de forma directa. Para encontrar la mejor combinación de hiperparámetros se emplea GridSearchCV con parámetros como el número de árboles y la profundidad máxima.
La validación se realiza con un enfoque de sliding window, que respeta el orden temporal de los datos, entrenando el modelo con una ventana de tamaño fijo y validando en el siguiente punto en el tiempo. En nuestro caso, indicamos que el tamaño de la ventana de entrenamiento es de 15 días y que el paso para mover la ventana también es de 15 días. Es decir, nuestro modelo entrena con un conjunto de datos que abarca 15 días consecutivos y luego valida con el siguiente período de 15 días.
Tras encontrar el mejor modelo realizamos un análisis de la influencia de variables. Se puede observar que las variables más influyentes son las tres variables de nuestro conjunto de datos relacionadas con el tiempo. De forma secundaria lo son las variables meteorológicas y las relacionadas con el tráfico del aeropuerto.
Para este modelo, se emplearon las librerías scikit-learn y pandas para la preparación de los datos. El preprocesamiento se basó en aplicar OneHotEncoder a las variables categóricas, mientras que las variables numéricas se mantuvieron sin transformación (passthrough).
La validación se llevó a cabo mediante TimeSeriesSplit, una estrategia de validación cruzada especialmente diseñada para series temporales, que garantiza el respeto del orden cronológico y evita fugas de información futura.
La optimización de hiperparámetros, como el número de estimadores y la tasa de aprendizaje, se realizó mediante una búsqueda en malla (GridSearchCV).
El modelo de XGBoost se construyó utilizando scikit-learn y pandas para el preprocesamiento de datos. Usamos OneHotEncoder para codificar las variables categóricas y se utilizaron las variables numéricas sin normalizar, dado que los árboles de decisión no son sensibles a la escala de las variables.
Para validar el modelo, utilizamos TimeSeriesSplit, respetando la secuencia temporal de los datos y evitando filtraciones de información futura. La búsqueda de los mejores hiperparámetros se realizó mediante GridSearchCV.
Una vez obtenidos los mejores hiperparámetros, entrenamos el modelo final con los datos de test , que son los mensajes más recientes.
Se han realizado distintos experimentos empleando redes neuronales densas bajo el framework de Tensorflow. Se ha tomado del conjunto de datos de entrenamiento como validación del 1 al 14 de enero (los mensajes más recientes). Asimismo, se han normalizado las variables numéricas y aplicado One-Hot Encoding a las variables categóricas.
Para la optimización de hiperparámetros se ha empleado la búsqueda bayesiana. Además, se ha conseguido una ligera mejora empleando una transformación logarítmica a la variable respuesta junto con MinMaxScaler.
Por otro lado, también se hizo un experimento probando a escoger únicamente un mensaje aleatorio por vuelo, pero los resultados empeoraron.
Asimismo, realizamos un experimento similar al anterior pero construyendo un modelo por pista. Se emplearon el mismo preprocesado y estrategia de búsqueda de hiperparámetros que en el modelo anterior.
Se han realizado distintos experimentos con redes neuronales LSTM, elegidas por ser útiles para tratar con datos de series temporales. Primero, se convierten los datos de entrada de un formato distribuido a un formato NumPy y realiza un preprocesamiento sobre la variable objetivo, aplicando una transformación logarítmica (ya que es una variable sesgada) seguida de un escalado entre 0 y 1. Luego se adapta el formato de las características para que sean compatibles con una red neuronal recurrente.
Para realizar el preprocesado de los datos se ha construido una pipeline con los siguientes pasos: StringIndexer, para aplicar índices numéricos a las variables categóricas; VectorAssembler para crear un vector con todas las características (features_raw) y MinMaxScaler para escalar los datos y facilitar el entrenamiento.
A continuación, construimos un modelo de red basado en capas LSTM, aplicando regularización, normalización por lotes, técnicas de abandono (dropout) y una capa de salida densa. Estas capas ayudan a evitar el sobreajuste. Este modelo se envuelve en un objeto compatible con scikit-learn para permitir su uso dentro de una búsqueda de hiperparámetros, que se hizo inicialmente. Después se realizaron más pruebas y por cuestiones de tiempo de entrenamiento de los modelos se usaron los hiperparámetros que encontramos como óptimos en el primer experimento.
Después, se configura una búsqueda en cuadrícula (GridSearchCV) combinada con validación cruzada específica para series temporales (TimeSeriesSplit), que respeta el orden de los datos para evitar contaminación temporal. Cada entrenamiento puede terminar antes gracias a un criterio de parada temprana basado en la pérdida (EarlyStop).
Finalmente, tras completar la validación cruzada, el mejor modelo se reentrena utilizando todo el conjunto de datos disponible, quedando listo para su evaluación final o para su uso en producción.
También se intentó añadir otra capa LSTM al modelo, pero la métrica era peor al modelo con una sola capa ajustado con hiperparámetros (MAE de 75).
Otro intento fue crear dos modelos, uno para los aviones con un tiempo de despegue relativamente corto (~237 segundos), y otro para los aviones con tiempos de despegue más largos, con los datos de entrenamiento separados en esos dos conjuntos. Para los datos de test se usó el mismo umbral pero aplicado a la columna time_at_holding_point, tomando como referencia el mayor valor para separarlos y saber qué modelo aplicar.
Se desarrolló un modelo predictivo para el tiempo de espera de despegue (takeoff_time) utilizando una arquitectura Transformer, aplicada a secuencias de mensajes de aeronaves enriquecidas con datos meteorológicos. El preprocesamiento, gestionado en gran parte con Apache Spark, fue clave: se aplicó una transformación logarítmica (log1p) a la variable objetivo sesgada, se crearon características cíclicas (seno/coseno) para variables temporales (hora, día), se codificaron las variables categóricas usando StringIndexer y OneHotEncoder (según cardinalidad), se imputaron valores nulos, y las numéricas se escalaron con StandardScaler antes de ensamblar el vector final de características. Estos datos se convirtieron a secuencias NumPy para TensorFlow/Keras.
El modelo Keras incorpora Codificación Posicional y dos bloques Transformer simplificados apilados, cada uno con Auto-Atención Multi-Cabeza y una red Feed-Forward básica con regularización (Dropout, L2). La salida pasa por GlobalAveragePooling1D y capas Dense para la regresión. Para el entrenamiento se empleó el optimizador Adam, la función de pérdida MAE (seleccionada por su robustez ante outliers), y parada temprana (EarlyStopping) monitorizando un conjunto de validación independiente (derivado de datos de prueba separados). Las predicciones se evalúan en segundos (MAE/RMSE) tras invertir las transformaciones logarítmica y de escalado.
Realizamos una separación de los datos de entrenamiento en función del tiempo de despegue (takeoff_time) para entrenar dos modelos de transformer. Para saber qué modelo aplicar a los datos de test, utilizamos el Q3 del umbral de tiempo de despegue del entrenamiento para clasificar las filas de prueba de acuerdo a la columna time_at_holding_point, es decir, el tiempo que lleva la aeronave situada en el punto de espera.
A la hora de entrenar, se entrenaron dos modelos diferentes con los subconjuntos obtenidos, el modelo low (para despegues con tiempos de espera bajos), y el modelo high, para despegues con tiempos de espera altos. El modelo de transformer elegido para ambos fue el mismo que en el apartado anterior por falta de tiempo para entrenar con búsqueda de hiperparámetros, pero sería una posible mejora para el rendimiento del sistema.
Sin embargo, nos hemos dado cuenta de que teníamos contaminación de futuro: agrupábamos todos los datos del mismo vuelo para ser clasificados en un solo modelo. Esto de forma indirecta provocaba que mensajes de aviones que solo llevan escasos segundos esperando pero que realmente tienen tiempos hasta el despegue altos se han estado prediciendo con el modelo high porque mensajes posteriores del mismo vuelo (que en el momento no se deberían saber) tienen tiempos esperando en el punto altos. Por lo tanto, las métricas del modelo eran mejores que las del resto pero no del todo correctas si solo tenemos en cuenta el contexto anterior.
Es por ello que no hemos incluido este modelo en la evaluación. Lo seguiremos desarrollando en los próximos días con el planteamiento correcto para separar los mensajes.
En esta sección vamos a comentar comportamientos que hemos observado en general en los modelos. Para las gráficas tomaremos como referencia los gráficos generados en nuestro dashboard de uno de los experimentos de red densa.
La distribución general de los errores tiende a estar sesgada hacia la izquierda, hacia los valores negativos. Esto lo que nos sugiere es que los modelos tienden a subestimar los tiempos de espera, especialmente los valores grandes. Una de las posibles razones de ello es que parte de los tiempos de espera altos se deban a factores que no hemos captado con nuestras variables.
En las gráficas siguientes podemos observar que las métricas varían en función de la zona del aeropuerto (puntos de espera / pistas). Contrastando esta información con la gráfica que comentamos anteriormente acerca de la ocupación de los puntos de espera, podemos observar que los puntos de espera en los que se obtienen los mejores resultados son precisamente aquellos de los que se contaba con un mayor número de casos. Ejemplos de estos puntos pueden ser K2, Y1, Y2, Y3, Z4, Z2, Z3. Por el contrario, hay otros puntos que en nuestro conjunto de datos son rara vez ocupados, en los cuales las métricas empeoran.
No se observan diferencias significativas en la distribución de errores por categorías de turbulencia.
Todos los modelos han exhibido un comportamiento similar en los mismos días del test, con un pico en el día 27 de enero.
| Modelo | RMSE | MAE | 36L/18R | 36R/18L | 32L/14R | 32R/14L |
|---|---|---|---|---|---|---|
| Heurística I | 128.60 | 89.15 | 91.94 | 74.67 | 112.82 | 79.80 |
| Heurística II | 118.40 | 82.60 | 84.27 | 63.95 | 114.51 | 73.65 |
| Regresión lineal | 104.01 | 73.65 | 77.19 | 58.22 | 96.10 | 64.10 |
| Random Forest | 104.86 | 69.74 | 71.26 | 54.38 | 94.32 | 64.50 |
| AdaBoost | 107.69 | 75.39 | 79.30 | 60.83 | 100.25 | 61.32 |
| XGBoost | 103.75 | 73.47 | 75.42 | 58.62 | 97.04 | 67.46 |
| Red densa | 114.18 | 72.09 | 73.95 | 55.00 | 99.38 | 65.92 |
| Red densa / pista | 148.55 | 82.66 | 77.17 | 78.90 | 107.80 | 68.89 |
| LSTM | 105.18 | 71.46 | 71.34 | 59.75 | 133.70 | 66.94 |
| Transformer | 106.18 | 67.75 | 66.31 | 55.19 | 94.54 | 60.81 |
Podemos observar que la mayoría de modelos presentan un MAE similar, con una mejora sobre las heurísticas de unos 10 segundos. Aún así, el error medio es bastante alto, teniendo en cuenta que la mediana de la variable respuesta son 174 segundos. Esta métrica se puede comprender mejor al analizar la distribución de la variable respuesta, ya que tiene una gran variabilidad y aunque se realiza un escalado y transformaciones para mitigarlo, hay suficientes datos atípicos como para que no esté justificado eliminarlos, pero sí que introducen ruido en las predicciones de los datos más “regulares” del modelo.
Destaca el hecho de que los modelos entrenados por pista presentan métricas mucho peores. Esto puede deberse a que hay pistas que tienen un número mucho menor de datos de entrenamiento y no se captan los patrones generales, lo que dificulta el entrenamiento de estos modelos más específicos, sobre todo para la pista 32L/14R.
Sorprendentemente, no se detecta una gran diferencia entre los modelos basados en árboles y los modelos más complejos basados en redes. Una posibilidad es que las variables tengan una estructura que responde muy bien a la clasificación que realizan los árboles, ya que se pueden separar y agrupar en categorías que tiene sentido que influyan positivamente o no en la variable respuesta (por ejemplo, las condiciones meteorológicas o si la pista está ocupada).
El modelo que comete menor error medio es el Transformer. Esto puede ser debido a que además de usar redes neuronales para capturar dependencias complejas de datos, son capaces de recoger y utilizar información del contexto de los datos (capturando el estado del aeropuerto con más precisión que los otros). Por lo tanto, este es el modelo elegido para pasar a la fase de producción.
Un aspecto a destacar es que el modelo XGBoost presenta el menor número de errores grandes, logrando el mejor RMSE. Además, algunos modelos superan al Transformer en ciertas métricas específicas para algunas pistas. Sin embargo, estas diferencias no son especialmente significativas en comparación con el Transformer, y consideramos que su consistencia general, junto con su mejor capacidad de escalado ante nuevos datos, justifican su elección como modelo principal. No obstante, el margen de mejora sigue siendo muy amplio.
