5.1 Data

For this research, a database with departures of scheduled commercial flights, domestic and international, of the entire Colombian airport network (58 airports) for the whole year 2018 was used (of the data provided by the Colombian aeronautical authority, only delays in flight departures are available). Such database includes statistical information related to 357,595 flights, with a total of 5,721,520 data entries (each flight data contains 16 associated flight attributes). Each flight attribute is associated with the following relevant information: (a) type of traffic (domestic or international); (b) airline; (c) origin airport (always national); (d) destination airport (national or international); (e) flight number (linked to the airline, route and frequency of said flight); (f) scheduled departure date; (g) scheduled departure time; (h) towing date (push-back); (i) towing time (push-back); (j) delay expressed in fraction of an hour; (k) delay expressed in minutes; (l) flight situation; (m) departure status (according to the standard of the country’s airport authority); (n) delay code (according to the Standard IATA Delay Codes-AHM 730 [72]), in which each code (89 possible) is associated with the reason that generates or causes said delay; (o) reason for the delay; (p) observations associated with the flight and/or its delay (by the airport authority).

Machine learning models require all input and output variables to be numeric. Encoding is a required pre-processing step when working with categorical data, as it is the case of our different flight’s attributes. Then, all the flight attributes or inputs must be converted into number via ordinal encoders functions. In addition, data scaling is a recommended pre-processing step when working with many machine learning algorithms, as it is this case with the ten models’ analysis. Unscaled target variables on regression problems can result in exploding gradients causing the learning process to fail. Also, data scaling is a recommended pre-processing step when working with deep learning neural networks as it is the case of MLP model. Therefore, the target flight attribute data called here time (flight) delay (DEM) must be standardised, that involve rescaling the distribution of data values (of DEM), then the mean of observed values is 0 and the standard deviation is 1. It is sometimes referred to as ‘whitening’. This can be thought of as subtracting the mean value or centreing the data. Like normalisation, standardisation can be useful, and even required in some machine learning algorithms when the input data has values with differing scales. Standardisation assumes the data fit a Gaussian distribution with a well-behaved mean and standard deviation, as it is in this case. But even this standardised Gaussian pre-processing it is recommended to be applied, if this expectation is not met, to get more reliable results.

5.2 Typology of implemented models

ML/DL models mathematically represent some types of universal approximators [Reference Guss and Salakhutdinov73] of the function f, such that it fulfills: Y = f (X). So, in the present work ten ML/DL models are proposed and used to map the Y output (delay time) with a finite and arbitrary (available attributes) set of X inputs (flight attributes, which can be up to ten in its most extensive version of the dataset), although of these ten models dare defined in Table 1, only a single DL-type model is developed and used, the so-called multi-layer perceptron (MLP) [Reference Vang-Mata74].

Table 1. Models used in this research

The MLP model, whose architecture is shown in Fig. 1, has been generated using the keras/tensorflow libraries, which consist of two dense layers (connections of all neurons in the previous layer with all neurons in the posterior layer) of 9 neurons in the first layer and 36 neurons in the second, using the activation function ReLU (rectified linear unit) [Reference Agarap75–Reference Arora and Barak78], so that neurons can learn the non-linear relationships between the output and the input (of each neuron). The ReLU function is defined as: y = max (0, x), in which y represents the neuron output and x the neuron input. The ReLU function has several advantages: its computation is easy and, therefore, training times are short; it also converges quickly and does not present saturations (a problem that does occur when other activation functions such as sigmoid and hyperbolic tangent are used by having a linear output); then ReLU neurons activation do not present the vanishing gradient problem (which prevent changing or updating the model weights because the derivatives of weights sensitivity are nearly zero); they are sparsely activated (sparsity feature), which simplifies the existing connectivity or couplings between all neurons; and, finally, allows the ML/DL model to learn the non-linear relationships between the input x and the output y.

The rest of the models developed in this work, and referred in the Table 1, are the following: the support vector machine (SVR) model, a very robust model in both linear and non-linear predictions because the dimensional space of the features has been increased [Reference Bao, Xiong and Hu79]; the Random Forest (RFR) model, which is an ensemble-type method that improves the results of any of its simplest constituent algorithms by combinations of these [Reference Breiman21, Reference Schonlau and Zou80]; the simple linear regression (LIR) model [Reference Montgomery, Peck and Vining81]; the bagging gradient (BGR) model, which improves predictions by producing multigroup by combinations and repetitions [Reference Rokach82–Reference Friedman84]; the extreme gradient boosting (XGR) model, which from weaker models generates stronger models using the descending gradient as an optimising method [Reference Hanif85, Reference Chen and Guestrin86]; the Huber Regressor (HBR) model, which is also very robust in training the data with outliers (or atypical values) by penalising them through ad hoc cost functions [Reference Wang, Zheng, Zhou and Zhou87, Reference Sun, Zhou and Fan88]; the extra trees (ETR) model, which implements many more decision trees to better fit learning than simple decision tree (DTR) methods [Reference Camana, Ahmed, Garcia and Koo89, Reference Geurts, Ernst and Wehenkel90]; the DTR model, which is the simple model of decision trees, which is reused due to its self-explanatory capacity of the patterns learned by the ML [Reference Izza, Ignatiev and Marques-Silva91]; the k-nearest neighbors (KNN) model, which predicts closest neighbours based on its K parameter [Reference Azadkia92]. The ‘R’ at the end of model names indicates that the model has been implemented for regression analysis. Table 1 lists, at the summary level, the models used.

5.3 Training ML/DL models

Training the models is part of the core ML process. In this work, each of the ten ML/DL models has been trained separately with the dataset under four possible scenarios: the first scenario with the ten full flight attributes (9 inputs + 1 output), the second scenario reducing scenario one to only six flight attributes (5 inputs + 1 output), the third scenario considering only flight attributes that could be known in advance (5 inputs + 1 output), and finally the fourth scenario reducing the third one to only four flight attributes (3 inputs + 1 output) (see Fig. 3). The reduced scenarios are aimed at evaluating whether the initial flight attributes of the dataset can be reduced. Another aspect to highlight from this paper is having carried out the statistical model training by implementing the k-fold and cross-validation techniques, to ensure that there are no biases in the training due to having chosen special subsets from the dataset for training/test, resulting in a more robust evaluation of the performance of ML/DL models.

The training has been implemented by dividing the dataset into batches of 128 samples each, subsequently, the ML/DL model processes the data forward (from the inputs to the output layer), to predict or estimate the model output values. At this point, with the error obtained by measuring the predicted value versus the target value, the training method goes backward (from the output layer to the input layer) to update the model weights (or connections between neurons) of all neurons of all layers of the model, through the backpropagation method [Reference Goodfellow, Bengio and Courville93–Reference Kuo and Chen95]. This method of backward propagation of the error starts obtaining the total error through the appropriate defined cost function. Here, the root-,eam-square-error (RMSE) applied metric results from comparing the values (of the batch) estimated by the model, with the real or target values collected from the known supervised features of the dataset. This core process of ML can be explained as the process distribution of the total error between the different contributions made by all the weights or nodes of the model, which is conducted in turn by the stochastic gradient descent (SGD) method [Reference Liu, Gao and Yin96, Reference Gower, Loizou, Qian, Sailanbayev, Shulgin and Richtarik97]. The SGD method calculates the partial derivatives of the error versus each of the weights, which is a way of measuring the model sensitivity of the error compared to each of the weights of the artificial neural network. The adjective ‘stochastic’ of the SGD method comes from the fact that the gradient is calculated by randomly chosen batches and not for the entire dataset.

The equations corresponding to MLP [Reference Vang-Mata74, Reference Dingari, Reddy and Sumalatha98] are those associated with neuron layers, in which a first layer corresponds to the input (represented by the X tensor), and the last output layer (associated with the Y tensor). Between them, the hidden layers to be implemented are inserted. Each layer contains a variable number of neurons (or nodes) that can be activated with the activation functions, which can be of various types [Reference Hastie, Tibshirani and Friedman99]. The input and output a neuron are shown schematically in Fig. 4.

Figure 3. Layers and neurons type architecture of the developed MLP model. Source: authors.

Figure 4. Neuron connections of MLP networks. Source: Ref. [Reference Kingma and Ba100].

In Fig. 4, x _i represents the inputs of feature i of the input layer, w _ij represents the weight of feature i of the previous layer in connection with neuron j of the posterior layer, $\varSigma$ represents the transfer function (sum of all the inputs coming from the outputs of all the neurons of the previous layer with neuron j corrected with the b _j bias of layer j, which produces the net output (transfer function + bias) Net _j . The chosen activation function $\sigma$ is then applied to each neuron (for example, tanh or sigmoid, although ReLU has been chosen here), which in turn produces the output O _j . From here, the connections of the output of this neuron are established again, together with all the neurons of its layer, with the neurons of the next layer, until reaching the last output layer Y _j (which can be several neurons of output, although in the case of regression prediction it is constituted by a single neuron). This inflow-outflow is called ‘forward flow’ or ‘propagation’ (from the input layer through all the neurons in the intermediate layers, to the final output layer). $\theta$ _j is the threshold value of neuron j, to activate its output function (if it exists).

In this work, the RMSE has been used as the cost function. Equation (1) defines this error with the norm 2 || ||, in which E is the error, and n is the number of samples in the dataset, y (x) is the actual output (provided by the dataset under supervised learning method) and $\hat{y}(x)$ is the output estimated by the model.

(1) \begin{align} E = \sqrt {\frac{1}{n}\;\mathop \sum \limits_x \|{{\left( {y\left( x \right) - \hat y\left( x \right)} \right)\!\|}^2}} \end{align}

In ML, once the cost function is obtained, the fundamental process for the automatic adjustment of the weights of the model is implemented via BP [Reference Goodfellow, Bengio and Courville93]. Such method calculates the gradient of said total error (in this case the RMSE) compared to each of the weights or connections of the neurons of each layer of the model, for which the rule of chain derivation applies, to propagate said error backward or to distribute the error between all the model weights, neuron by neuron, layer by layer. If the method converges, the new adjustment of weights of the model is sought to produce a smaller error when it is iterated with the next batch of input samples X. This is the essence of ML, minimising some loss function (cost or error), adjusting the contributions of the weights of the model neurons connections per each dataset batch and repeating the BP process adjustment to all dataset iterations implemented. Learning that is sufficiently guaranteed when using the SGD method by the theorems of the existence of the null derivation at the points of minimum error. The following group of Equations (2) shows the derivation chain rule in artificial neuron networks (ANN) [Reference Aggarwal101, Reference Hassoun102].

(2a) \begin{align} \frac{{\partial E}}{{\partial {w_{ij}}}} = \frac{{\partial E}}{{\partial {o_j}}}{\rm{\;}}\frac{{\partial {o_j}}}{{\partial {w_{ij}}}} = \frac{{\partial E}}{{\partial {o_j}}}{\rm{\;}}\frac{{\partial {o_j}}}{{\partial ne{t_j}}}{\rm{\;}}\frac{{\partial ne{t_j}}}{{\partial {w_{ij}}}}\end{align}

(2b) \begin{align} \frac{{\partial ne{t_j}}}{{\partial {w_{ij}}}} = \frac{\partial }{{\partial {w_{ij}}}}\left( {\mathop \sum \limits_{k = 1}^n {w_{kj}}{o_k}} \right) = \frac{\partial }{{\partial {w_{ij}}}}{\rm{\;}}{w_{ij}}{o_i} = {\rm{\;}}{o_i}\end{align}

The problems of searching the cost function minimums are solved through different optimisers, generically known as the SGD [Reference Hastie, Tibshirani and Friedman99]. These optimisers are responsible for avoiding the solution from getting trapped in any of the local minimums that the cost function presents when it is not a convex function, for which each of them has implemented different strategies to find the global optimum. In the present work, the adaptive moment estimation (ADAM) optimiser has been used because it is one of the most efficient optimisers in terms of its high convergence speed and robust stability [Reference Kingma and Ba100], which is an improved version of the RMSProp that implements a variable/adaptive speed (or learning step) and a moment or inertia that tries to avoid getting caught in the local minima of the cost function. The following group of Equations (3) describes the update of the weights by the ADAM optimiser [Reference Kingma and Ba100].

(3a) \begin{align} {m_w}^{\left( {t + 1} \right)} \leftarrow {\beta _1}{m_w} + \left( {1 + {\beta _1}} \right){\rm{\;}}{\nabla _{w\;}}{l^{\left( t \right)}}\end{align}

(3b) \begin{align} {v_w}^{\left( {t + 1} \right)} \leftarrow {\beta _2}{v_w}^{\left( t \right)} + \left( {1 - {\beta _2}} \right){\rm{\;}}{({\nabla _{w\;}}{l^{\left( t \right)}})^2}\end{align}

(3c) \begin{align} {\hat m_w} = \frac{{{m_w}^{\left( {t + 1} \right)}}}{{1 + {\beta _1}^{\left( {t + 1} \right)}}}\end{align}

(3d) \begin{align} {\hat v_w} = \frac{{{v_w}^{\left( {t + 1} \right)}}}{{1 + {\beta _2}^{\left( {t + 1} \right)}}}\end{align}

(3e) \begin{align} {w^{\left( {t + 1} \right)}} \leftarrow {w^t} - \eta \frac{{{{\hat m}_w}}}{{\sqrt {{{\hat v}_w}} + \varepsilon }}\end{align}

Table 2. List of flight attributes, and their participation in the different scenarios of the model

Figure 5. RMSE comparison resulting from the learning obtained for each of the ten models, transforming the delay time into a standard distribution (0 mean and 1 for standard deviation), referred to scenario 1 resulting from choosing ten flight attributes, nine inputs (TRF, AEL, ORI, DES, NVU, CDE, MDE, OBS, EAC) and one output (DEM). Source: authors.

Figure 6. Comparative analysis of the ten models’ performance (RMSE), in its original scale of minutes of delay, referred to as scenario 1, choosing ten flight attributes, nine inputs (TRF, AEL, ORI, DES, NVU, CDE, MDE, OBS, EAC) and one output (DEM). Source: authors.

Figure 7. Relative importance coefficients in increasing order, of the nine full flight attributes considered as inputs of scenario 1, when determining the flight delay time (DEM) as output. Source: authors.

Figure 8. Comparative results of the RMSE obtained by the 10 ML/DL models when trained under reduced scenario 2, with only five flight attributes. It can be seen how the model learning performance (in terms of RMSE) is roughly identical to full scenario 1. Source: authors.

Figure 9. RMSE comparative analysis of the 10 models under a more realistic dataset scenario 3 (five inputs: TRF, AEL, ORI, DES, NVU and one output: predicting delay (DEM)), in which the best RFR model obtains an RMSE in predicting flight delay of 33.8 minutes. Source: authors.

Figure 10. Rank of the importance of the five input flight attributes considered in scenario 3 (TRF, DES, ORI, NVU, AEL) in increasing order. Source: authors.

Figure 11. RMSE comparison of models for reduced scenario 4, with three inputs (ORI, NVU, AEL) and DEM as time delay to be predicted, in which the best RFR model obtains the same RMSE of the prediction of flight delay of 33.8 minutes than on scenario 3. Source: authors.

Figure 12. Time distribution of the flight departure delay. Source: authors.

Figure 13. Unimodal distribution (density) function of a flight delay (at departure). Source: authors.

In the Equations (3), ${w^{\left( {t + 1} \right)}}$ represents the weights in iteration t + 1; $\eta $ represents learning speed or learning step; ${\rm{\;}}{\nabla _{w\;}}{l^{\left( t \right)}}$ represents the gradient of the cost function l with respect to the weights; $\epsilon$ is a minimum value to avoid dividing by 0 (in keras 10^-7 is used); ${\beta _1}$ y ${\beta _2}$ are the forgetting factors for the gradient or first moment and the second moment respectively; ${m_w}^{\left( {t + 1} \right)}$ is the moving mean of the weights or first moments in the iteration t + 1; ${v_w}^{\left( {t + 1} \right)}$ is the moving variance (or second moment) of the weights in the iteration t + 1. The variables with hat indicate that they are estimators and the superscripts with the variable t correspond to each of the iterations (or learning cycles). This is the process of fitting the MLP/DL model. On the other hand, in models based on decision trees, the adjustment is made mainly by induction and pruning methods of the decision trees [Reference Aggarwal101, Reference Hassoun102].

The dataset learning process is repeated in this work 150 times (iterations) for the entire dataset and only for the MLP/DL model (considered sufficient to achieve an asymptotic or stable RMSE value), which marks the maximum fine adjustment possible of the model in the learning of the dataset, or until a sufficient adjustment of the error is obtained for the rest of the non-DL models, characterised by belonging to the ensemble or boosting or stacking type, in which they improve the simplest models of decision trees by their appropriate combinations [Reference Rokach82].

To give more robustness to the value of the metrics obtained from the RMSE, k-folds techniques have been implemented [Reference Pedrycz and Chen103, Reference Chang and Lin104]. Specifically, three groups have been applied to divide the dataset between three different ways of grouping the training and test data subsets. The training data, as its name implies, is used to train the model and the test data to evaluate the learning performance of the model. However, the model can never see the test data during the learning process to make an objective evaluation. In addition, with cross-validation techniques this training is repeated a variable number of times. Because target delay flight attribute was standardised as necessary step of data pre-processing, the loss or cost function RMSE of each ML/DL model it is also characterised statistically, in terms of a mean and a standard deviation, resulting all of this in a more robust way performance measured compared to the case when a single configuration of training and test subset is used. It should be emphasised that all these methods are stochastic, due both to the use of the SGD method, and to the own stochastic nature of the dataset, so their results are not deterministic, but rather are probabilistic.

5.4 Statistical examination

Before training, different examinations or statistical controls were applied to have a feature insight. For this purpose, different tests has been implemented to determine if flight attribute distributions are Gaussians (or normal) (Shapiro-Wilk, D’Agostino and Anderson-Darling tests), or whether flight attributes are not correlated with the output (linear Pearson, Spearman Rank, Kendall Rank, and chi-squared tests), or whether the flight attributes are not stationary (ADFuller test), or whether flight attributes have equal distribution as the output (student’s test, paired student, ANOVA, Mann-Whitney U, Kruskal-Wallis H, and Friedman tests). These tests, of statistical characterisation of flight attributes, serve to know amongst other things, whether their learning process will be easy or which data transformation must be required to facilitate it. These techniques are an enhancement of dataset preprocessing to facilitate ML training.

5.5 Determination of the importance of the flight attributes

To rank the importance of flight attributes in determining flight delay time, the permutation importance (PI) method [Reference Altmann, Tolosi, Sander and Lengauer20, Reference Breiman21] has been implemented. This method consists of applying inspection techniques to trainable estimators; each model proposed and configured for this work offers a different response regarding the selection of the most significant attributes of the flights. The PI method consists of randomly permuting the features among themselves and measuring in each case how the output (delay time) varies through the training of the new dataset (permuted or sometimes called corrupted in the scientific literature) with each ML/DL model, which allows establishing the coefficient of importance or the sensitivity of the output (time delay) to each feature j called c _j (see Equation (4)).

(4) \begin{align} {c_j} = s - \;\;\frac{1}{K}\;\mathop \sum \limits_{k = 1}^K {s_{k,j}}\end{align}

In Equation (4) s is the result of training the ML model when it is trained with the D dataset without permuting. Next, each feature j associated with a column of dataset D is randomly permuted with the other features to produce a permuted version of the dataset $\tilde D$ (j, k), in which the subscript k represents the repetitions used to manage the features, K is the maximum number of permutations, and s _{k, j} is the new result of training the ML model for each feature j repeated k times with the permuted (or corrupted) dataset $\tilde D$ (j, k).