Characterizing departure delays of flights in passenger aviation network of United States

Flight delay happens every day in airports all over the world. However, systemic investigation in large scales remains a challenge. We collect primary data of domestic departure records from Bureau of Transportation Statistics of United States, and do empirical statistics with them in form of complementary cumulative distributions functions (CCDFs) and transmission function of the delays. Fourteen main airlines are characterized by two types of CCDFs: shifted power-law and exponentially truncated shifted power-law. By setting up two phenomenological models based on mean-field approximation in temporal regime, we convert effect from other delay factors into a propagation one. Three parameters meaningful in measuring airlines emerge as universal metrics. Moreover, method used here could become a novel approach to revealing practical meanings hidden in temporal big data in wide fields.

Flight delay happens every day in airports all over the world. However, systemic investigation in large scales remains a challenge. We collect primary data of domestic departure records from Bureau of Transportation Statistics of United States, and do empirical statistics with them in form of complementary cumulative distributions functions (CCDFs) and transmission function of the delays. Fourteen main airlines are characterized by two types of CCDFs: shifted power-law and exponentially truncated shifted power-law. By setting up two phenomenological models based on mean-field approximation in temporal regime, we convert effect from other delay factors into a propagation one. Three parameters meaningful in measuring airlines emerge as universal metrics. Moreover, method used here could become a novel approach to revealing practical meanings hidden in temporal big data in wide fields. Flight delays happen every day in the airports all over the world. Passengers in different countries suffer from such ubiquitous events, and economy and traffics of many cities are largely harmed by these events. Delays of individual flights seem to be random at a glance. Actually, they obey certain statistical laws [1] taking a longterm delay records for a large number of flights into account. In general, some delay originating from an upstream flight spreads to downstream flights, which is particularly evident when an aircraft continues its tasks of successive flights. This phenomenon is defined as delay propagation(DP) [2][3][4][5]. To alleviate the effect caused by DP [3,4,7], we need to characterize them quantitatively first, which is still a big difficult task now although empirical statistics [1,3,4,6,[8][9][10] and mechanism modeling [2,5,[11][12][13][14][15] have appeared recently. The key problems remaining challenging us lie on: what kind of statistical law we can learn from the large scale historical records of flight-delay events? And what are possible mechanisms underlying them in terms of propagation and other factors? They motivate our present work from both aspects of empirical statistics and theoretical modeling.
The causes of flight delays have been classified into five categories by the US Bureau of Transportation Statistics (BTS) [16]: (1)Air carrier; (2)extreme weather; (3)national aviation system(NAS)referring heavy traffic volume and ait traffic control; (4)security; (5)late arriving aircraft. Just as cited by Baumgarten et al. [17], the propagation factor(PF)(delay propagation caused by late arrival of last flight of the same aircraft) occupies 35 percent in this classification. In sense of propagation, we can roughly divide them into just two categories: PF and non-propagation factor(NPF).
The statistics on empirical data of flight delays makes up the primary step on the characterization of DP [1,4,6,8,9]. Analytical model [2,5,[11][12][13][14][15] serves as the second step on the investigation of DP. In view of successful contribution from these two steps, we are called to go on the further step, i.e., to make general models to explore generic mechanisms based on airline-specific big data of flight delays.
We start to investigate empirical laws from collecting data of departure delays (DD). Domestic DD data in all the airports over the whole United States are downloaded from BTS of America [16], including flight numbers, tail numbers, scheduled departure and arrival time, and actual departure and arrival time of each flight in the American passenger network. And DD in 2014 of each flight is calculated as the difference between actual departure time and the scheduled one. Then we obtain the statistical results of probability distribution functions (pdf) of DD time for each airline, respectively. Fig.1 in SupplementaryInformation(SI) illustrates the pdfs for 14 main airlines of America. Furthermore, integrating pdf from certain DD to infinity makes up complementary cumulative distribution functions (CCDFs) [26] for every airline. These CCDFs are shown with color filled circles in Fig.1 and Fig.2, respectively. Each of them is set up a model to understand the underlying mechanism.

Model 1
As a rough approximation, we start to describe DD of flights in a unified way, i.e., to account the effect from NPF by converting them into that equivalent to PF. Then, we merely consider the effect of PF in the derivation of CCDFs of airlines. At first, we make two assumptions based on data observation from the spirit of mean field approach. As assumption (1), current DD caused by PF is assumed smaller than last DD of the same aircraft in general cases over one-year statistics, which takes account of the DP and the efforts against it made by the staffs of the management and operation. In practical manipulation, we take DP as follows: for a flight with DD l (in minutes), if last flight operated by the same aircraft is larger than l we assume that the current DD is caused by PF. This is because real delays incurred by PF are not discernable at present [12]. Based on such an assumption, we let the number of flights with DP in the unit delay interval near l be divided by the total number of the flights with the delay larger than l , and define this ratio as the probability of delay transmission q(l ). Then where n B (l ) is the number of delayed flights in the unit interval near l caused by DP, and n(l ) is the total number of delayed flights per delay interval near l , therefore, cumulated number of delayed flights N (l ) = ∞ l n(l )dl . By formula (1) we account the number of propagationcaused delayed flights in the unit interval of (l as the transmission result of all delayed flights with from l to ∞. In the present phenomenological model, the form of the function q(l ) is unknown and needs to be dug out by data-mining later. As a simple result of direct observation, the probability of large delay obviously decrease with l monotonously, so the probability q(l ) of delay transmission per delay interval should be a monotonous decreasing function as shown in Fig.2 of SI, for instance of the real data of airline AA. However, it should be finite in the limit l → 0 . Therefore, it is a reasonable choice to set as its tentative form, with α, β as airline -specific parameters determined by data-mining according to formula (2). In our phenomenological model, the parameters in q(l ) have to be obtained from the data. Taking the American Airlines (AA) as an example, we calculate the value of q(l ) for 237 aircrafts in the whole year of 2014 using formula (2). By taking ∆l = 10 min and select (0,10) (10,20) ..., as statistical interval, then letting the vari-able l be equal to the median of each interval, we plot empirical results in Fig.2 of SI.
As assumption (2), we suppose that propagationcaused fraction of delayed flights per delay interval keeps a constant k , so that we have Comparing formula (1) with (3), we see When the tentative transmission function q(l ) is substituted into it, we obtain Therefore, where α 0 = α/k . And the CCDF of flight DD reads where N 0 is the total number of flights of an airline over a year. Therefore, we obtain shifted power-law (SPL) [22][23][24] to describe the CCDFs of DD of flights. Formula (8) in phenomenological Model 1 is checked for 14 airlines(See Fig.1 of SI). It is valid for AA, MQ, F9, DL, HA and AS airlines with each pair of two parameters β 1 and α 0 well adjusted to fit the primary data from a specific airline, where β 1 means adaptive β and α 0 in formula (8)(see Fig.1). However, it fails to cover the remaining 8 airlines, which demands other types of mechanisms. In the panels of Fig.1, real CCDFs of these 6 airlines are plotted with the filled circles in black, green, blue, wine red, dark yellow and purple, red lines representing analytic results in formula (8) with corresponding shift parameters β 1 = 40, 50, 90, 40, 30 and 90 fit linear parts of empirical data well. While the transmission function q(l ) is checked in its integrated form F (l ) following formula (2) with real data, where β 2 means β in formula (8) adaptive to fitting integrated q(l ) (F (l ) = l 0 q(l )dl ) in Fig.3(also collected in Table (1) of SI). The tailderivations of red lines from filled circles in both Fig.1 and Fig.3 are discussed in SI. In fitting these CCDFs from primary data with the present mechanism of DP, we are actually converting the effect from NPF into an equivalent one from PF, since the CCDFs illustrated by colored symbols contain effect from all kinds of delay factors, while the red lines are based on the all-PF assumption potentially.   (14). Filled circles in black, blue, green, dark yellow, purple, olive, wine red and orange illustrate CCDFs from real data ( Fig.1 of SI) of airlines B6, VX, UA, US, WN, EV, FL and OO, respectively. Red lines illustrated the CCDFs in formula (14). By adapting parameters (λ, β1 and r) in Model 2 with modified mechanism of DP, red lines from formula (14)

Model 2
To understand the CCDF data from the remaining 8 airlines, we need a new mechanism other than what is for those 6 ones. At first, the assumption of a constant k for the propagation-caused fraction of delayed flights per delay interval in model 1 can not be valid for complex cases. Actually, formula (5) would be still valid if k behaves in an infinitesimal in the same order of q(l ) as l → ∞ according to L ′ Hospita rule. Therefore, we can  expect that k has a similar form of q(l ). Tentatively, we let where g and h are both phenomenological constants. Substituting it into formula (5), we obtain then dN (l ) N (l ) = (−gα + −r l + β )dl (11) where, r = α(h − gβ). We have lnN (l ) = −gαl −rln(l +β)+c 0 = lne −gαl +ln(l +β) −r +c 0 (12) where c 0 is a constant. Therefore, where λ = 1/gα, c 1 = e c0 . And CCDF of DD of flights reads where c 2 = c 1 /N 0 . Different from SPL in formula (8), in this mechanism we have a CCDF in the form of exponentially truncated shifted power-law (ETSPL) [25] for the remaining 8 airlines. When we fit empirical results with parameters in analytic functions based on DP, we convert effect by NPF into equivalent effect by DP factor just as in the way of Model 1.
In last model, we assumed that last DD of the same aircraft is greater than current one . In reality, however, the effects of NPFs including cancellations may overlap with PF. At present, we are not able to separate two kinds of factors in the empirical data [12]. Therefore, the empirical statistical results for 8 airlines in Fig.2 definitely contain effect from all kinds of delaying factors, including that of last delay smaller than the current one of the same aircraft. Obviously, assumption (1) in last model leads to a smaller value of n B (l ) than the actual one. Based on this consideration, we modify the definition of n B (l ) in Model 1: the total number of the flights with the delay near l per delay interval, but transmitted by last flight's DD greater than (l − m), where m describes average additional DD affecting current flights. By counting larger interval as effective one affecting the current delay of flights from equivalent DP, we compensate the deficient counting of contribution from all factors. We demand β 2 approaching values of β 1 in Fig.2, where β 2 means β in formula (14) adaptive to fitting F (l ) in Fig.4. Airlinespecific parameters m emerge out as the new metrics see Fig.4).

Aviation meanings of key parameters
Shift parameters β 1 in fitting CCDFs can serve as the feature quantities in minutes characterizing specific airline with the aviation meanings. We define the average flight delay absorption [27,28] time L (in minutes)(See Table (2) of SI) of propagation delay as the average differences between last delay and the current one of the same aircraft for all flights with larger last delays in each airline. From direct statistics on empirical data, we find that the relation for the absorption time L versus β 1 is a negatively correlated linear function. One can see from VX which has the lowest β 1 (0.53) and the smallest number of flights except HA hence an ignorable propagationdelay effect so that its CCDF can decay more exponentially than that of any other airlines. From both Fig.2 and Fig.4, we reckon that the airline with smaller β 1 is much more effective to reduce the delay from last delay to current ones, i.e., it has a larger ability to absorb DP, which leads to the number of flights with small delay account for a larger proportion in all the delayed flights (See Fig.6). For examples behind, from Table 1 in SI, among 6 airlines characterized by Model 1, F9 and AS behave less effective in the absorption of last delays with their larger β 1 compared with others.
In comparison with CCDFs of 6 airlines fit by SPLs, integrated CCDFs of 8 airlines drop down steeper as increased delay l , and they have smaller β 1 , which means that the factor (l + β 1 ) −r behaves roughly as l −r for the cases l > λ because λ > β 1 (see Table 2 in SI) is always satisfied. Actually, since smaller β 1 means larger absorption L(β 1 ), therefore the probability to have propagationcaused delay from larger l decreases with l quickly. While the chance for NPF to act on delay increases. In this category of delay cases, λ acts more importantly than β 1 . The smaller λ an airline has, the smaller interval for it to have a propagation-induced delay with longer delay than that of current one. Among 8 airlines, VX behaves as the typical case since it has both small β 1 and λ. While FL has larger probabilities with the effects of DP than others since it has both larger r as the power in SPL and larger parameter λ emphasize DP dominating tendency. However, such a picture does not mean smaller fraction of flights with DP since the actual distribution for a flight to delay is on average biased to the side with l < λ, even with the fraction as high as 90 percent, which is shown in Fig.6  set to be 1 4 λ for clear look in Fig.6, with the arrows indicating the positions of critical delay λ separating the PF dominating range from that of NPFs.

CONCLUSIONS
In the present work, two new mechanisms are presented to model flight delays based on the big data from BTS of the United States. Two types of CCDFs describing delays are obtained from two models with analytical derivation, they fit empirical data well. Fourteen American airlines can be sorted into two categories. One is more dominated by the propagation of delays, and its CCDF follows SPL. The other is more dominated by PF under a critical delay λ, while it is more dominated by NPF above it, and its CCDF follows ETSPL.
Starting from the mechanism of delay propagation, we convert the effect from NPF and cancellations into it, not only find the key parameter β 1 for the mechanism of propagation, but also find the compensation interval m describing the overlapping effects with NPF and smaller last delays than current ones. Moreover, key parameter λ separating the effect of PF from NPF also emerges out. Three quantities demonstrate their practical meaning in aviation, which can serve as new metrics to measure the quality of management and operation of airlines. In addition, by accounting the current delays as the transmission result of all previous delays incurred by all kinds of factors, we carry the spirit of mean-field approach in temporal regime. The approach that converting other factors into an equivalent propagation one before extracting its quantitative aviation meaning verifies its validity in dealing with temporal big data, which is hopefully applicable to many topics in wide other fields. At the tails of CCDF curves with colored symbols from empirical data, obvious dropping down from linear fitting lines are always observed. From Fig.2 we see that the probability for a flight to be transmitted by a very long delay is very small, it only fluctuates a little bit above zero. So the behavior of actual flights do not follow SPL in the limit of large l . Actually, flights may be canceled when they delay too much. And the longer one delays, the larger probability for it to be canceled. However, formula (8) is derived under a potential prerequisite, i.e., the mechanism of DP is valid for any delay, which is not true, especially for very large delays. In other words, assumption (1) becomes ineffective in such tails.
The effectiveness of SPL shown in formula (8) seems to rely on the validation of the tentative function q(l ). It is checked by the comparison with the results from empirical data. The functions q(l ) are integrated for smooth look in Fig.3 and Fig.4. The filled circles in black, blue, green, dark yellow, purple, olive, wine red and orange represent integrated q(l ) of airlines AA, MQ, F9, DL, HA and AS based on the definition (formula (2)), each with shift parameters β 1 = 35, 40, 90, 30, 30 and 80, respectively. The red lines represent the tentative function in the integrated form of fractional function q(l ) = α l+β2 from formula (3). And the integrated form F (l ) with variable upper bound reads: F (l ) = ∞ l α(ln(l + β 2 ) − lnβ 2 ). They fit linear parts of filled circles representing empirical data of these 6 airlines well. Due to the same reason in fitting curves in Fig.3, F (l ) deviates from increasing red lines and gets saturated when l becomes large enough.
The parameters used in the integrated CCDFs and the function F (l ), are collected in Table 1 which shows obvious inconsistence between corresponding columns β 1 and β 2 for airlines AA, MQ, DL and AS, respectively. n B (l , l + ∆l ) from formula (1) to (3) is defined as the number of current delayed flights per delay caused by last delayed flights and this number is actually not discernable in the practical data, which has been pointed out [? ] before the present work. So, in our simulation, we have to actually substitute it by the number of last flights with larger delay than the current one. This manipulation is carried out in Fig.3 and Fig.4 for plotting filled circles from real data, which is not always true because such substitution results in smaller quantities of n B (l , l + ∆l ) than the existing ones in fact. However, filled circles in Fig.1 and Figf.2 in the main text are plotted by the data directly from the primary records.
And by fitting red lines to them the parameters α 0 and β 1 are obtained. Neither the data symbols nor the fitting line in Fig.1 and Fig.2 in the main text rely on the substitution of the parameters α 1 and β 2 into them from Fig.3 and Fig.4. Actually integrated CCDFs in Fig.1 and Fig.2 have included the effect from all kinds of delay factor calibrated into equivalent DP measure in formula (5). Therefore, such two parameters α 0 and β 1 serve as the baselines for the comparison. Difference between β 2 and β 1 actually reveals the degree of effectiveness of that manipulation in the simulation and the assumption of all-DF, which is also used to understand the earlier deviation of lines F (l ) from the colored data symbols in the variable l + β 2 in the panels of Fig.3 and Fig.4.