1 Introduction

Transportation systems have a crucial importance for countries because of their social and economic impact. The air transportation in particular is closely linked to the economic development of the areas in which it expands. This is why it is very important for policy makers to ensure a smooth development, even—and especially—in areas where the traffic increase forecasts are the highest. Indeed, the air traffic system will get closer and closer to its actual capacity, especially in Europe and in the USA where the traffic is expected to increase by 50% in the next 20 years [7]. As a consequence, it is important for the air traffic management world (i) to forecast the consequences on the current infrastructures and procedures and (ii) to find the appropriate solutions to cope with this increase. For this reason, large investment programs like SESAR in Europe and NextGen in the USA have been launched.

Apart from airport capacity, one of the important bottlenecks for the increasing traffic flow will be the management of the airspace and sectors, where the controller needs to actively separate flights in order to avoid conflicts. However, solving conflicts in areas of high traffic complexity is already nowadays a demanding task. With the increase in traffic, the cognitive capacities of air traffic controllers will likely reach their limits and drastically increase the number of conflicts or force to cap the capacities of the sectors. As a consequence, navigating through the European sky will become more and more difficult in the future and will require more careful planning capacities for the network manager and for the airlines. In other words, the airspace is becoming a scarce resource, especially in congested situations, like, for example, during major shutdowns of large areas (extreme weather, strikes, volcano eruptions, etc.).

It is thus expected that the airlines will compete fiercely for two of the most important resources: time and space. More specifically, it is foreseen that the allocation of slots at the airports will change and will be structured as a market for companies. On the other end, the airspace will be more densely populated and the airlines will also compete for it. From the point of view of the transportation companies, this increases the effort required to find better route allocation strategies, whose success depends, among other things, on the strategies adopted by the other users.

Motivated by these considerations, in this paper we study a model of the allocation of the flight plans on the airspace from the point of view of the dynamics on a complex network. This point of view is fruitfully used in different fields, like dynamics of epidemiology, information propagation on the Internet, percolation, opinion spreading, systemic risk, etc. [2, 3]. Recently, an increasing attention is being devoted to the network description of transport systems [13], in particular the air transport system [4, 6, 9, 12, 14, 15, 18, 19]; for a recent review, see Zanin and Lillo [20].

The model describes the strategic allocation of flight plans on an idealized airspace, described as a network of interconnected sectors. The sectors are capacity constrained, i.e. a maximum number of flights can be simultaneously present. This implies that companies might not get their optimal flight plan and thus they will fall back to suboptimal solutions, for which they will develop different strategies. By using two different strategies for companies, stylized version of low-cost companies and so-called traditional carriers, we show how different factors explain the satisfaction of different types of company. Some of these factors (the network topology and the pattern of departing times) can be regulated externally by the policy maker, while others (the fraction of airlines of each type) depend on the airline population and on market forces. We then study the evolutionary dynamics of the populations by considering a “reproduction” rate (i.e. the capacity to expand business) of a company, which is based on its past satisfaction, namely how well its past flights were allocated. In other words, the satisfaction plays the role of a fitness function.

Our main findings can be summarized as follows. By using a simplified baseline model, we prove analytically that in the static framework there exists a Nash equilibrium for the fraction of airlines of different types. Extensive numerical simulations on more complicated settings, including empirically calibrated models, confirm the robustness of this finding. We then show that the equilibrium point is distinct from the optimal point for the system and the (static) equilibrium tends to favour the archetype for low-cost companies. When the evolutionary dynamic is considered, we show that the model can be recast into a replicator-type equation. We study the more interesting finite population case, by showing analytically that the dynamics can be described by a linear stochastic differential equation (SDE). Numerical simulations of the model and the analytical study of the SDE show that the noise due to finite populations lead to a further shift in the equilibrium point in the direction of favouring low-cost companies. Finally, we show that the equilibrium is reached more slowly because of finite populations.

The paper is organized as follows. In Sect. 2, we present briefly the model. In Sect. 3, we investigate a simplified version of the model that we are able to investigate analytically. In Sect. 4, we present some results on the static equilibrium of the model, concerning the behaviours of the airlines in different situations, and in Sect. 5 we investigate the population dynamics in an evolutionary environment. Finally, we draw some conclusions in Sect. 6.

2 The Model

In this section, we present our model, which has been introduced in Gurtner et al. [11]. The implementation of the model is open and can be freely downloaded for any non-commercial purpose.Footnote 1 A more detailed version of the model with a tactical part is also availableFootnote 2 [10].

The model describes the strategic allocation of trajectories in the airspace. Mimicking what is done in the European airspace, the model considers airlines submitting their flight plans to the network manager (NM). The NM checks whether accepting the flight plan(s) would lead to a sector capacity violation. If this is not the case, the flight plan is accepted, otherwise it is rejected and the airline submits the second best flight plan (according to its utility or cost function). The process goes on until a flight plan is accepted or a maximal number of rejected flight plans are reached, and in this case the flight is cancelled. The NM keeps track of the allocated flights and checks violations of newly submitted flight plans responding in a determined way to the requests, without making counter-propositions.

Airspace The airspace is modelled as a network of sectors. Each sector has a capacity C, here fixed to 5 for all sectors. Some of these sectors contain airports, and the geometry of a flight plan is a path connecting two airports. Specifically, we use a triangular lattice with 60 nodes. In order to avoid paths having exactly the same duration, which could lead to ties in the optimization, we sample the crossing times between sectors from a log-normal distribution so as to have a 20 min average and a very small variance (inferior to \(10^{-4}\) min). Unless specified otherwise, we fix the number of airports to 5. In the following, we present results in which we drew 10 times the position of the airports randomly, then ran 1000 independent simulations on each of these realizations.

Real airspaces are clearly more complex. Topological properties of the real networks of sectors have been investigated in Gurtner et al. [12]. In “Appendix A”, we present some robustness checks of the model by considering two more realistic set-ups. In the first one, we consider a scale-free network of airports, i.e. not all the airports are equivalent in terms of number of flights/destinations, but hubs and spokes are present. In the second, we use real ECAC data to construct the network of sectors, the sector of capacities, the origin/destination frequencies, and the waves structure. We find that the results are indeed similar to the ones presented thereafter for the stylized model, where a much more controlled setting is used.

Airlines The main agents of the model are the airline operators (AOs) who try to obtain the best trajectories for their flights and the NM who accepts or rejects the flight plans. In the simplified version, we assume that the quality of a trajectory depends on its length (the shorter, the better) and the discrepancy between the desired and actual departing/arrival time (the smaller, the better).Footnote 3 Companies might be different in the relative weight of two components in their cost or utility function. Companies caring more about length are called of type “S” (for shifting) companies, since when their flight plan is rejected by the NM, they prefer to delay the flight, but keeping a short trip length, mainly because of cost—fuel and airspace charges. Companies caring more about departure punctuality are called of type “R” (for rerouting), since when the flight plan is rejected they prefer to depart on time even if they need to use a longer route to destination, mainly to avoid disruption in their network operations—e.g. connecting flights. As a consequence, S companies can be thought as “low cost”, whereas R companies are more like “traditional” ones. Note that each AO has only one flight. Hence, in our model the optimization takes place after the previous, larger strategic allocation of flights where AOs decide or not to operate the route, with which aircraft, etc. For this reason, all the optimizations here are independent from each other for each flight. However, it is important to stress that the capacity constraints of sectors create a dependency between the accepted flight plans.

More quantitatively, for each flight an AO chooses a departing and arrival airport, a desired departing time, \(t_0\), and selects a number \(N_\mathrm{fp}\) of flight plans. The kth flight plan, \(k=1,\ldots ,N_\mathrm{fp}\), is the pair \((t_0^k, \mathbf{p}^k)\), where \(t_0^k\) is the desired time of departure and \(\mathbf{p}^k\) is an ordered set containing the sequence of sectors in the flight plan. The flight plans are selected by an AO according to its cost function. In our model, it has the form

$$\begin{aligned} c(t_0^k,\mathbf{{p}}^k) = \alpha \mathcal{L}(\mathbf{{p}}^k) + \beta (t_0^k-t_0), \end{aligned}$$
(1)

where \(\mathcal{L}( \mathbf{{p}}^k)\) is the length of the path on the network (i.e. the sum of the lengths of the edges followed by the flight). We also assume that flights are only shifted ahead in time (\(t_0^k \ge t_0\)) by an integer multiple of a parameter \(\theta \) which is taken here as 20 min (all durations in this article is in minutes unless specified otherwise). The parameters \(\alpha \) and \(\beta \) define the main characteristics of the company. Given the discussion above, R companies have \(\beta /\alpha \gg 1\), while S companies have \(\beta /\alpha \ll 1\).

Departing waves An important determinant of the allocations is the desired departing time \(t_0\) chosen by the AO. We assume that departing times are drawn from a distribution inside the day characterized by a certain number of peaks or waves. This is indeed typically observed at most airports.

We define first \(T_d=1440\) as the length of the “day” (in minutes), i.e. the time window of departure for all flights. In this time window, we define \(N_\mathrm{p}\) peaks of \(T_0 = 60\) min, by setting a time \(\Delta t\) between the end of the peak and the beginning of the next one (thus, \(N_\mathrm{p} = \lceil T_d/(\Delta t + T_0) \rceil \) ). Then, we define a total number of flights \(N_\mathrm{f}\) and divide them equally between peaks. In the following, we also use the corresponding hourly density \(d = N_\mathrm{f}/24\), i.e. the average number of flights per hour.

Dynamics Given a mixed population of AOs of different types, at each time step, an AO is selected randomly.Footnote 4 The AO chooses the departing and arrival airports and the desired departing time \(t_0\) for its flight, drawing it from the departing time distribution. It then computes the \(N_\mathrm{fp}\) best flight plans for the flight according to its cost function and submits them, one by one in increasing order of cost, to the NM. The NM accepts the first flight plan which does not cross overloaded sectors, i.e. at already maximal capacity. If none of the \(N_\mathrm{fp}\) flight plans is accepted, the flight is rejected and cancelled.

Metrics The metric measuring the satisfaction (or fitness) of a company about a given flight f is

$$\begin{aligned} \mathcal{S}_\mathrm{f}=c_\mathrm{f}^\mathrm{best}/c_\mathrm{f}^\mathrm{accepted}, \end{aligned}$$
(2)

where \(c_\mathrm{f}^\mathrm{best}\) is the cost of the optimal flight plan for the flight f according to the AO cost function (i.e. the first flight plan to be submitted for the flight), and \(c_\mathrm{f}^\mathrm{accepted}\) is the cost of the flight plan eventually accepted for this flight. If no flight plan has been accepted, we set \(\mathcal{S}_\mathrm{f}\) to 0. Note that \(\mathcal{S}_\mathrm{f}\) is always between 0 and 1. The value 1 is obtained when the best flight plan is accepted.

Since the AOs have only one flight, the satisfaction of a flight is also the satisfaction of its company. When several companies are of the same type (same ratio \(\beta /\alpha \)), we make use of the average satisfaction across them. Thus, \(\mathcal {S}^\mathrm{S}\) and \(\mathcal {S}^\mathrm{R}\) are the average satisfactions of S and R companies, respectively. We use also the average satisfaction across all flights as a measure of the global satisfaction of the system:

$$\begin{aligned} \mathcal{S}^\mathrm{TOT} = f_\mathrm{S} \times \mathcal{S}^\mathrm{S} + f_\mathrm{R} \times \mathcal{S}^\mathrm{R}, \end{aligned}$$
(3)

where \(f_i\) and \(\mathcal{S}^{i}\) are the fraction of flights and the average satisfaction of company i, respectively, and \(f_\mathrm{S}+f_\mathrm{R}=1\). Finally, we consider the difference of satisfaction between S companies and R companies to estimate how well they perform with respect to the other type:

$$\begin{aligned} \Delta \mathcal {S} = \mathcal {S}^\mathrm{S} - \mathcal {S}^\mathrm{R}. \end{aligned}$$

The main features of the model have been presented in Gurtner et al. [11] where only the static setting with two airports has been investigated with numerical simulations. Gurtner et al. [11] showed that with a single type of company, there is a (congestion) transition, much like the congestion observed in other transport systems [16], e.g. car traffic, when the number of flights becomes too large. When two extreme types of companies (R and S) are competing for the airspace, Gurtner et al. [11] used numerical simulations to show that there exists a unique fraction of mixing corresponding to a stable Nash equilibrium. The strategies are interacting positively, leading to an absolute maximum in satisfaction for the overall system at a mixing fraction different from 0 and 1.

Compared with Gurtner et al. [11], the main innovations presented in the following are:

  • We present first a baseline version of the model, which can be treated analytically, and we show explicitly the existence of the equilibrium and how it depends on the model parameters.

  • We consider simulations of a more realistic setting with multiple airports.

  • We consider an evolutionary setting where the capability of a type of company of continuing its business depends on its past satisfaction.

3 Analytical Treatment of a Baseline Model

Given the complexity of the model, the results presented in the following sections are obtained through numerical simulations. In order to understand the results more mathematically, in this section we present a simplified baseline model which can be treated analytically. The main result, namely the existence of an equilibrium in the mixed population case, is derived here and later compared with the numerical results of the complete model.

In the baseline model, we take the simplified case of an airspace with only two airports in which the flights can only travel from one to the other. We also assume that the flights travel instantaneously between the two airports, i.e. they are loading all the sectors in between exactly at their time of departure. In this simplified setting, individual flights do not contribute to the satisfaction when they are rejected and contribute by one unit otherwise. Hence, the satisfaction is \(\mathcal {S}= n_\mathrm{a}/N\), where \(n_\mathrm{a}\) is the number of accepted flights and N is the total number of flights. This simplified setup exhibits nevertheless the main features of the full model.

3.1 Pure Populations

Consider first the case where all the companies are of type S, i.e. they shift their flight plans in time. For \(\Delta t=0\) (all waves are consecutive), shifting the flight plan does not yield any improvement, since there are flights departing in the next wave. As a result, the satisfaction is \(n^\mathrm{S}_\mathrm{a} \simeq C N_\mathrm{p} =C \lceil \frac{T_d}{\Delta t + T_0} \rceil \).

When \(\Delta t\) increases and the waves are parting from each other, S companies start to have some opportunity to shift their flight plan if they are rejected during the first waves. They allocate more and more flight plans until they hit their maximum number of flight plans \(N_\mathrm{fp}\). The number of accepted flights is thus

$$\begin{aligned} n^\mathrm{S}_\mathrm{a} \simeq \left\{ \begin{array}{ll}C\lceil \frac{T_d}{\Delta t + T_0} \rceil \left( 1 + \frac{\Delta t}{T_0}\right) &{} \quad \hbox {if} \quad \Delta t < \theta (N_\mathrm{fp}-1) \\ C\lceil \frac{T_d}{\Delta t + T_0} \rceil \left( 1 + \frac{\theta (N_\mathrm{fp}-1)}{T_0}\right) &{} \quad \hbox {if} \quad \Delta t \ge \theta (N_\mathrm{fp}-1),\end{array}\right. \end{aligned}$$
(4)

where we recall that \(\theta \) is the increment by which the flight plans can be shifted. The satisfaction of S companies is initially oscillating with \(\Delta t\), and for \(\Delta t \ge \theta (N_\mathrm{fp}-1)\) it is non increasing.

On the contrary, R companies cannot shift their flight plan and, since the flights are instantaneous, the number of accepted flights is \(n^\mathrm{R}_\mathrm{a} \simeq C N_\mathrm{p} N_\mathrm{fp}^u= C\lceil \frac{T_d}{\Delta t + T_0}\rceil N_\mathrm{fp}^u\), where \(N_\mathrm{fp}^u\) is the number of flight plans among the \(N_\mathrm{fp}\) which have no sector in common. As a consequence, the satisfaction of companies R is monotonically decreasing with \(\Delta t\).

3.2 Mixed Populations

In the mixed populations case, a further complication appears, namely that the arrival of the two types of companies is modelled as a random process and therefore different realizations of the process can lead to different values of the satisfaction. To get some intuition on the result, we consider first the case of the maximum attainable satisfaction and then consider the expected value of the satisfaction over the distribution of airlines arrival.

In the mixed population case, a fraction \(f_\mathrm{S}\) of companies are of type S and a fraction \(1-f_\mathrm{S}\) are of type R. The parameter \(f_\mathrm{S}\) will be called mixing parameter in the following. Let us consider first S companies: on the one hand, since R companies cannot shift in time, S companies are not competing with them when shifting in time and thus the periods between waves are completely available to companies of type S. On the other hand, there is a competition for the periods within waves and thus on average S companies will be able to allocate only \(N_\mathrm{p} C f_\mathrm{S}\) flights within the waves. Hence, the resulting number of flights accepted for S companies is:

$$\begin{aligned} n^\mathrm{S}_\mathrm{a}\simeq \left\{ \begin{array}{ll}C \lceil \frac{T_d}{\Delta t + T_0} \rceil \left( f_\mathrm{S} + \frac{\Delta t}{T_0}\right) &{} \quad \hbox {if} \quad \Delta t < \theta (N_\mathrm{fp}-1) \\ C \lceil \frac{T_d}{\Delta t + T_0}\rceil \left( f_\mathrm{S} + \frac{\theta (N_\mathrm{fp}-1)}{T_0}\right) &{} \quad \hbox {if} \quad \Delta t \ge \theta (N_\mathrm{fp}-1)\end{array}\right. \end{aligned}$$
(5)

Clearly this formula reduces to Eq. 4 when \(f_\mathrm{S}=1\). The behaviour of \(n^\mathrm{S}_\mathrm{a}\) is now a bit different from equation 4, since, apart from the oscillations, it is increasing with \(\Delta t\) when \(f_\mathrm{S}\simeq 0\). For intermediate values of \(f_\mathrm{S}\), the function has now a maximum in \(\Delta t\). This effect will be seen also in the full model (see the right panel of Fig. 2 in Sect. 4). Note that \(n_\mathrm{a}^\mathrm{S}\) increases with \(f_\mathrm{S}\), but slower than the total number of S companies \(N f_\mathrm{S}\), leading to a monotonic decrease in the total satisfaction with \(f_\mathrm{S}\).

Likewise, R companies face competition within each wave and only \(N_\mathrm{p} C (1-f_\mathrm{S})\) of the best flight plans are allocated by them in average. On the contrary, their suboptimal rerouted flight plans are specific to them and thus face no competition from S companies. This is true under the condition that the rerouted flight plans do not cross the best flight plan at any point—because it is also used by S companies. Hence, the number of accepted flights for R companies is:

$$\begin{aligned} n_\mathrm{a}^\mathrm{R} \simeq N_\mathrm{p} C (1-f_\mathrm{S}) + N_\mathrm{p} C (N_\mathrm{fp}^u - 1) = C \lceil \frac{T_d}{\Delta t + T_0} \rceil (N_\mathrm{fp}^u - f_\mathrm{S}). \end{aligned}$$
(6)

On the contrary of \(n_\mathrm{a}^\mathrm{S}\), this quantity decreases monotonously with \(\Delta t\). It also decreases with \(f_\mathrm{S}\), but slower than the total number \(N (1-f_\mathrm{S})\) of R companies. As a consequence, the satisfaction of R companies increases with \(f_\mathrm{S}\), i.e. decreases with its own fraction \(1-f_\mathrm{S}\).

We are now interested in seeing how the relative satisfaction of S and R companies depend on the wave structure. It is clear that there is always a root for \(\Delta \mathcal {S}\) as a function of \(f_\mathrm{S}\). In fact,

$$\begin{aligned} \Delta \mathcal {S}\propto \left\{ \begin{array}{ll} C \lceil \frac{T_d}{\Delta t + T_0} \rceil \left( 1 + \frac{1}{f_\mathrm{S}}\frac{\Delta t}{T_0} - \frac{N_\mathrm{fp}^u - f_\mathrm{S}}{1-f_\mathrm{S}} \right) &{} \quad \hbox {if} \quad \Delta t< \theta (N_\mathrm{fp}-1) \\ C \lceil \frac{T_d}{\Delta t + T_0} \rceil \left( 1 + \frac{1}{f_\mathrm{S}} (N_\mathrm{fp}-1)\frac{\theta }{T_0} - \frac{N_\mathrm{fp}^u - f_\mathrm{S}}{1-f_\mathrm{S}} \right) &{} \quad \hbox {if} \quad \Delta t\ge \theta (N_\mathrm{fp}-1)\end{array}\right. \end{aligned}$$
(7)

For \(\Delta t< \theta (N_\mathrm{fp}-1)\), there is a single value of \(f_\mathrm{S}\) for which \(\Delta \mathcal {S}= 0\), which is given by:

$$\begin{aligned} f_\mathrm{S}^* = \frac{\Delta t/T_0}{\Delta t/T_0 + N_\mathrm{fp}^u -1}. \end{aligned}$$

Since \(N_\mathrm{fp}^u\) is at least equal to 1, and usually larger, the value of \(f_\mathrm{S}^*\) is somewhere between 0 and 1, which means that there is always an equilibrium. Note, however, that for \(\Delta t = 0\) the equilibrium is in fact \(f_\mathrm{S} = 0\). For \(\Delta t \ge \theta (N_\mathrm{fp}-1)\), the value of the root is given by the same expression, replacing \(\Delta t/T_0\) by \((N_\mathrm{fp}-1)/(\theta /T_0)\). In this case, there is always a root in (0, 1) when \(N_\mathrm{fp}^u>1\).

Interestingly, the equilibrium value depends on \(\Delta t\), which means that different wave structures lead to different relative advantages of companies of type S with respect to those of type R. More specifically, the \(f_\mathrm{S}^*\) increases monotonically with \(\Delta t\), until it reaches a plateau for \(\Delta t/T_0 \ge \theta (N_\mathrm{fp}-1)\) where it does not depend on \(\Delta t\) anymore.

3.3 Effect of Randomness on the Arrival of Companies

The above equations give a first good approximation of the satisfaction for each type of company, but it is easy to see that they are incorrect in some cases, in particular for small values of \(\Delta t\). Indeed, the arrival of companies is random, forcing us to reevaluate the simple average behaviour described above. For instance, when \(\Delta t=0\), allocating the first C companies in a wave to S or to R companies is very different. In both cases, C best flight plans are allocated, but the subsequent S companies are blocked because they cannot shift their flight plan, whereas the subsequent R companies can reroute. In other words, when allocating C flights of type S and then C flights of type R, 2C flights are accepted. If C flights of type R are allocated first and then C flights of type S, only C flights are accepted in total. To capture this effect, it is necessary to go into the details of the allocation by computing expected values using probability distributions.

Moreover, another effect plays a role for high values of \(\Delta t\). In this case, only the first time periods after the wave can be reached by all the flights in the wave, but the last time period can only be reached by those which are departing late in the wave. However, many of these flights have already been allocated during the first time periods. The number of time periods which can be reached by all the flights is \(m = \lceil (\theta (N_\mathrm{fp}-1)/T_0)\rceil - 1\). For \(\Delta t/T_0 < m\), all the time periods are fully allocated, but after that the last one is only filled with the remaining flights S which are late enough, which correspond to a fraction \((1-m+\theta (N_\mathrm{fp}-1)/T_0)\) of flights.

Let us denote \(n = N/N_\mathrm{p}\) the deterministic number of flights in each wave. Then, given that \(n_\mathrm{S}\) flights of type S are to be allocated in a wave, and that \(n_\mathrm{S}'\) of them are to be allocated in the first C flights, for a single wave:

$$\begin{aligned} n_\mathrm{a}^\mathrm{S} \simeq \left\{ \begin{array}{ll} n_\mathrm{S}' + C \Delta t/T_0 &{} \quad \hbox {if} \quad \Delta t/T_0 < m \\ n_\mathrm{S}' + C (m-1) + \Gamma &{} \quad \hbox {if} \quad \Delta t/T_0 \ge m\end{array}\right. \end{aligned}$$
(8)

where

$$\begin{aligned} \Gamma =\min (C, (n_\mathrm{S}-(n_\mathrm{S}'+C(m-1))) (1-m+\theta (N_\mathrm{fp}-1)/T_0)) \end{aligned}$$
(9)

and with the extra condition that \(n_\mathrm{a}^\mathrm{S}<n_\mathrm{S}\).

The probability to have \(n_\mathrm{S}\) flights among the n in the wave is described by the binomial distribution with parameter \(f_\mathrm{S}\), and the probability that \(n_\mathrm{S}'\) of the first C flights is of type S is described by the hypergeometric distribution. Therefore, the expected value of the number of accepted flights is

$$\begin{aligned} \mathbb {E}[n_\mathrm{a}^\mathrm{S}] = \sum _{n_\mathrm{S}=0}^n \sum _{n_\mathrm{S}'=0}^C n_\mathrm{a}^\mathrm{S}(n_\mathrm{S}, n_\mathrm{S}') B(n_\mathrm{S}; n, f_\mathrm{S}) H(n_\mathrm{S}'; n, n_\mathrm{S}, C) \end{aligned}$$

and the expected satisfaction of company S is:

$$\begin{aligned} \mathcal {S}^\mathrm{S} = N_\mathrm{p} \times \mathbb {E}[n_\mathrm{a}^\mathrm{S}]/(N f_\mathrm{S}). \end{aligned}$$

Following the same reasoning, the number of accepted R companies given the number \(n_\mathrm{R}\) R companies in a wave and the number \(n_\mathrm{R}'\) of them in the first C allocated ones is given by \(n_\mathrm{a}^\mathrm{R} = n_\mathrm{R}' + (N_\mathrm{fp}^u-1)C\), and its expected value is

$$\begin{aligned} \mathbb {E}[n_\mathrm{a}^\mathrm{R}] = \sum _{n_\mathrm{R}=0}^n \sum _{n_\mathrm{R}'=0}^C n_\mathrm{a}^\mathrm{R}(n_\mathrm{R}, n_\mathrm{R}') B(n_\mathrm{R}; n, 1-f_\mathrm{S}) H(n_\mathrm{R}'; n, n_\mathrm{R}, C), \end{aligned}$$

with the corresponding expected satisfaction:

$$\begin{aligned} \mathcal {S}^\mathrm{R} = N_\mathrm{p} \times \mathbb {E}[n_\mathrm{a}^\mathrm{R}]/(N (1-f_\mathrm{S})). \end{aligned}$$

Figure 1 illustrates the analytical results of the baseline model with their numerical simulations. The left panel shows the difference of satisfaction between the two types of companies as a function of \(f_\mathrm{S}\) for different values of \(\Delta t\), whereas the right panel shows the corresponding total satisfaction. As explained more in detail in the following sections, the roots of the curves in the left panel represent stable equilibrium points for the system, which are distinct from the optimal points of the system—the maxima of the curves in the right panel.

The agreement between the analytical computations and the simulations is very good, showing that the approximations made have a negligible effect. The larger discrepancies are close to \(\Delta t = 0\), where the analytical model predicts a constant difference in satisfaction, whereas the simulations produce a slightly decreasing function. However, both are well below the zero line and as a consequence the equilibrium is located at \(f_\mathrm{S}=0\).

Fig. 1
figure 1

Left: difference of satisfaction between S companies and R companies as a function of the fraction \(f_\mathrm{S}\) of S companies for different values of \(\Delta t\). Right: total satisfaction in the system as a function of \(f_\mathrm{S}\) for different values of \(\Delta t\). The dashed lines are the results of the analytical model, whereas the circles are the results of the numerical simulations (Color figure online)

Full size image

4 Static Equilibrium of the Full Model

We now study the simulations of the full model to investigate the impact of competition between airlines in different environments. Indeed, in a more realistic environment with multiple airports, the analytical model becomes intractable and simulations need to be performed. We study directly the mixed population case. The airspace is now more complex than in the previous section, with five airports, randomly chosen routes between them, non-instantaneous sector crossing times, and bidirectional flight plans between airports.

Note that the topology of the network as well as the number of airports can have non-trivial effects on the results. We show for instance in “Appendix B” that the satisfaction of the system is a function of the number of flights on the airspace and the overlap between the available paths in terms of number of sectors, which is a direct consequence of the topology of the considered airspace.

4.1 Satisfaction of Companies

We first study how the satisfaction of each type of company depends on the wave pattern and on the population composition. For this, we fix the number of flights, \(N_\mathrm{f} = 24 \times d = 480\) and we change the proportion \(f_\mathrm{S}\) of S companies. We also change the structure of the wave pattern, by changing the parameter \(\Delta t\).

Figures 2 and 3 show the satisfaction of the two types of company as a function of \(\Delta t\) and \(f_\mathrm{S}\). The results for R companies are quite intuitive and are consistent with those of the baseline model described in Sect. 3. These companies are better off when they are competing with a large fraction of S companies (Fig. 3 left) and when there are more waves, i.e. when \(\Delta t\) is small (Fig. 2 left). This is expected, since more waves means more “space” for companies when the number of flights is fixed. Moreover, R companies have a stronger dependency on the mixing parameter when the number of waves is small, i.e. \(\Delta t\) is large, because of an increased competition. Note that the different plateaus present on the plot are due to the fact that for these ranges of parameter, the number of waves is constant and they are far from each other. In particular, for \(\Delta t \in [720, 1320]\), there are only two peaks, which come slowly apart as \(\Delta t\) increases. Since they are sufficiently apart, the flights from the previous wave do not interact with the flights in the next one, and the satisfaction does not change with \(\Delta t\). In other words, the first flight from the second wave departs after the last flight from the first wave arrives.

Fig. 2
figure 2

Satisfaction of R (left) and S (right) companies as a function of the time \(\Delta t\) between waves. Different lines refer to different values of \(f_\mathrm{S}\) (Color figure online)

Full size image
Fig. 3
figure 3

Satisfaction of R (left) and S (right) companies as a function of the mixing parameter \(f_\mathrm{S}\). Different lines refer to different values of \(\Delta t\) (Color figure online)

Full size image

The satisfaction of S companies is more complex, since for small values of \(f_\mathrm{S}\), their satisfaction is not monotonous with \(\Delta t\). This behaviour is explained by the following trade-off. On the one hand, the larger the \(\Delta t\), the less waves there are. Hence, companies are competing effectively with a higher number of other companies (which is the reason behind the decreasing curve of R companies in the left panel of Fig. 2). On the other hand, S companies try to delay their flight if the first flight plan is rejected. This means that if the waves are too close to each other, the delayed flight plans will likely conflict with the flights in the next wave. For this reason, their satisfaction increases at the beginning when the waves comes apart and then decreases when the waves are further apart and in smaller number.

Note that the decrease in satisfaction due to concentration within waves is less significant when S companies compete with many R companies. Indeed, in this case, their increasing concentration within a wave is of little importance for them, because they can always shift their flight plan two or three times to get out of the wave and not conflict anymore with R companies. For this reason, their curve is monotonous with \(\Delta t\) for high values of \(f_\mathrm{S}\). More strikingly, their satisfaction is higher for very high \(\Delta t\) than for very small ones if \(f_\mathrm{S} \ll 1\). All these results are consistent with the ones of the baseline model and presented in Sect. 3.

Hence, companies are reacting differently to different wave structure because they are sensitive to different mechanisms. The interplay of the mechanisms leads to interesting patterns that translate in interesting behaviours when framing the model in an evolutionary environment.

4.2 Global Satisfaction

Figure 4 shows the difference of satisfactions \(\Delta \mathcal {S}\) between S and R companies as a function of the mixing parameter \(f_\mathrm{S}\) and for several values of \(\Delta t\). In the left panel, the density of flights is quite small (\(d=20\)), corresponding to the one used in Figs. 2 and 3. In the right panel, we show the result for a much higher density (\(d=80\)), corresponding to a severely congested airspace.

Fig. 4
figure 4

Difference of satisfaction \(\Delta S\) between S companies and R companies versus the mixing parameters for different values of \(\Delta t\). Left: low density of flights, \(d=20\). Right: high density of flights, \(d=80\) (Color figure online)

Full size image

The difference of satisfaction \(\Delta \mathcal {S}\) between the two populations strongly depends on the two parameters. At both densities, the first values of \(\Delta t\) are clearly crippling population S, since in this configuration \(\Delta \mathcal {S}\) is always negative. This is due to the fact that very frequent waves prevent S companies to delay their flight, whereas R companies can find an available path by suitable rerouting. For higher values of \(\Delta t\), the situation becomes more favourable to S companies, since the difference is usually positive. The details of the variations of the difference are quite complex with the two parameters, but it is clear it is always decreasing monotonically with \(f_\mathrm{S}\). The point where it crosses 0 varies with \(\Delta t\), but not wildly (except for small \(\Delta t\)).

It is worth noting that these curves can be considered as fitness curves for two populations competing for the same resources in a given environment. Assuming that the higher fitness affects positively future reproduction rate (i.e. the possibility of continuing and expanding business), we study in Sect. 5 the dynamics of the two populations in an evolutionary framework. Here, we simply recall that the points where the difference of fitness curves vanishes are equilibrium points for the dynamics. The existence of a single root (as in Fig. 4) shows that there is only one equilibrium point (apart from the two absorbing states at \(f_\mathrm{S}=0\) and \(f_\mathrm{S}=1\)). The slope of \(\Delta \mathcal {S}\) at its root measures the stability of the equilibrium. Since the slope is negative, the equilibrium is stable. In other words, when the proportion of S companies is too high, their satisfaction/fitness decreases, thus giving a lower reproduction rate for them, favouring R companies and driving back the system towards the equilibrium.

Another important question is whether the equilibrium point is optimal also for the system. For this reason, we compute also the global satisfaction, Eq. 3, which is the average satisfaction of all the flights. A higher global satisfaction means that globally resources are better allocated, leading to increased profits for airlines and possibly better service for passengers. In the left panel of Fig. 5, we show global satisfaction as a function of the mixing parameter \(f_\mathrm{S}\) for different values of \(\Delta t\). The first conclusion is that the global satisfaction is usually better for \(0<f_\mathrm{S}<1\) than for pure populations. This is expected, because we saw that each population performs better against the other one, as is typical when different populations have different niches and thus their interaction is beneficial for both. The second conclusion is that for all values of \(\Delta t\), there exists a unique maximum and its position varies with \(\Delta t\).

Fig. 5
figure 5

Left: global satisfaction as a function of the mixing parameters for different \(\Delta t\). Right: value of \(f_\mathrm{S}\) maximizing global satisfaction and the equilibrium point as a function of \(\Delta t\). The error bars represent the uncertainty on the positions of the maximum and the point of equilibrium due to the uncertainty on the value of the satisfaction. In particular, they are quite large for the positions of the maxima because \(\mathcal {S}^\mathrm{TOT}\) is quite flat around its maximum for some values of \(\Delta t\), as can be seen in the left panel (Color figure online)

Full size image

On the right panel of Fig. 5, we plot both the value of \(f_\mathrm{S}\) at the global optimum, extracted from the left panel, and the equilibrium point, extracted from the left panel of Fig. 4. Both exhibit similar variations. For small values of \(\Delta t\), both the equilibrium point and the global optimum are at \(f_\mathrm{S}\simeq 0\). When \(\Delta t\) increases, S companies increase their advantage against R companies, because they are not affected by the next wave. Then, both values decrease, stabilizing at a value \(f_\mathrm{S} \gtrsim 0.5\), showing the greater advantage of S companies when the departing pattern is composed by well-separated waves.

More importantly, both curves are clearly distinct for \(\Delta t \gtrsim 100\), even considering error bars. This is an important result, because it shows that the equilibrium mixing condition is not the optimal at the global level. In particular, the evolution of the system towards its equilibrium mixing would tend to favour drastically population S, whereas the global optimum would be reached with a much smaller market share of S companies. This is exactly where policy makers should step in and issue policies driving the system to the optimum.

5 Evolutionary Dynamics

In the previous section, we interpreted the satisfaction of each company as its fitness when competing with the others for the same resources—namely, time and space. Interpreting these fitnesses as the capability of expanding business, i.e. a reproduction rate, it is possible to develop a dynamical evolutionary model for studying the dynamics towards equilibrium and its fluctuations, as well as the role of finite size populations.

We assume that the population size at time \(t+1\) of a company depends on its satisfaction at time t. In order to keep the simulations under reasonable computational time and following what is done in evolutionary biology models [17], we keep the total population fixed. This means that only the mixing parameter \(f_\mathrm{S}\) is changing between time t and \(t+1\). For the reproduction rule, we use an exponential reproduction, i.e. the rate of reproduction of a population is proportional to its fitness and its current population. Combined to the fixed population conditions, this leads to a discretized version of the so-called replicator model [17]:

$$\begin{aligned} f_\mathrm{S}^{t+1} = f_\mathrm{S}^{t} + \Delta \mathcal {S}_t\, f_\mathrm{S}^t\, (1-f_\mathrm{S}^t), \end{aligned}$$
(10)

where \(f_\mathrm{S}^t\) is the mixing parameter at time t and \(\Delta \mathcal {S}_t\) is the difference in satisfaction between S and R companies at time t. In the following, \(\Delta \mathcal {S}_t\) is simply called the fitness function (of S). In the simulations, we also choose to keep the number of companies of each kind to a minimum of 1. This ensures that the equilibria at \(f_\mathrm{S} = 0\) and \(f_\mathrm{S} = 1\) do not act as absorbing barriers (sinks). Indeed, since the populations are finite, a small non-null \(f_\mathrm{S}^t\) could lead to exactly 0 company S, which leads in turn to \(f_\mathrm{S}^{t'} = 0\) for all \(t'>t\). Analogously, the same happens when \(f_\mathrm{S}^t = 1\). All the other parameters (\(\Delta t\), number of airports, airspace structure, etc.) are kept constant throughout the reproduction process, i.e. the environment is stable.

Fig. 6
figure 6

Evolution of the mixing parameters with the generations, averaged over 100 realizations (and only one network realization). The blue lines are the averages, the violet lines are exponential fits, and the error bars are the average standard deviations. The coefficients of determination of the regressions are over 0.98. Left: \(\Delta t = 23\times 60\). Right: \(\Delta t = 0\) (Color figure online)

Full size image

In a finite population case (see for example [1]), the term \(\Delta \mathcal {S}_t\) in Eq. 10 depends on the specific realization of the process and it is therefore a stochastic variable. To study the dynamics towards equilibrium, we use the linearization \(\Delta \mathcal {S}_t\simeq \gamma (f_\mathrm{S}^t-f_\mathrm{S}^*)+\sigma \eta \), where \(f_\mathrm{S}^*\) is the root of \(\Delta \mathcal {S}_t\), \(\sigma \) is a parameter going to zero with population size, and \(\eta \) is a Gaussian iid variable. Under this assumption, in “Appendix C” we show that the dynamics of \(f_\mathrm{S}^t\) can be described by a linear SDE. We find that the dynamical equilibrium point \(f_\mathrm{S}^{eq}>f_\mathrm{S}^*\) and the convergence to equilibrium are exponential in time and slower than in the infinite population case.

We tested these analytical conclusions with numerical simulations. Figure 6 shows the dynamics of \(f_\mathrm{S}^t\) for two distinct values of \(\Delta t\) obtained with numerical simulations. The solid blue lines are averages over 100 runs, and the solid violet lines are the results of an exponential fit with the functional form:

$$\begin{aligned} f_\mathrm{S}^t=f_\mathrm{S}^{eq}+(f_\mathrm{S}^0-f_\mathrm{S}^{eq})e^{-t/\tau }, \end{aligned}$$
(11)

directly inspired by the behaviour of Eq. 12. Both lines are well fitted (\(R^2>0.98\)), and the equilibrium is clearly reached in both cases. On the left, there is only one wave of departure; thus, S companies have an advantage and the point of equilibrium is \(f_\mathrm{S}^{eq} >0.5\). On contrary, the figure on the right shows that when there are no waves (\(\Delta t=0\)), S companies are very disadvantaged and \(f_\mathrm{S}^{eq} \simeq 0\). Both figures are roughly consistent with left panel of Fig. 4, where the root of \(\Delta \mathcal {S}\) is close to 0 when \(\Delta t = 0\) and close to 0.7 when \(\Delta t = 23\times 60\).

In Fig. 7, we plot the position of the equilibrium point—computed by averaging the last 40 generations in each run—as a function of \(\Delta t\) and using the same parameters as in the right panel of Fig. 5. By comparing the curves in the two figures, we note that the one obtained with evolutionary dynamics displays larger values of \(f_\mathrm{S}^{eq}\) than the static one. This is again consistent with the model in “Appendix C” and is an important result, because the noise coming from the fitness function can drive the equilibrium even further from the global optimum than in the deterministic case.

Fig. 7
figure 7

Equilibrium value of the mixing parameter in the evolutionary setting as a function of \(\Delta t\) (Color figure online)

Full size image

We now consider how the system converges to the equilibrium and how external parameters like \(\Delta t\) influence the convergence. It is worth reminding the link between the fluctuations around the equilibrium and the shape of the fitness functions. In the continuous version of the replicator model, larger absolute slopes of the difference of fitnesses at its root translate into a higher stability and faster convergence to the steady state [17]. However, our system does not have a deterministic fitness function, since \(\Delta \mathcal {S}_t\) depends on the specific realization of the model. As shown in “Appendix C”, this additional noise affects both the fluctuations around the equilibrium and the time of convergence \(\tau \).

The magnitude of the fluctuations around the equilibrium, measured as the standard deviation of \(f_\mathrm{S}^t\) after the transient period, is shown in the left panel of Fig. 8. There is a weak trend towards larger fluctuations when \(\Delta t\) increases, but their magnitude reaches a plateau quickly. Note that the standard deviation is far from being negligible, implying that the fluctuations are typically 15% of the value of the equilibrium point. This means that the static analysis performed in Sect. 4 is far from revealing all the features of the model.

Fig. 8
figure 8

Left: Standard deviation of the value of \(f_\mathrm{S}\) when the equilibrium is reached against \(\Delta t\). Right: Typical number of generations before the equilibrium is reached (time to equilibrium) as a function of \(\Delta t\) (Color figure online)

Full size image

The time of convergence to equilibrium \(\tau \) of Eq. 11 is plotted in the right panel Fig. 8. This time scale is quite high for small values of \(\Delta t\)—where the fluctuations are small—but decreases to a small value (around seven or eight generations) when \(\Delta t\) increases—where the fluctuations are high. As already stated, the magnitude of the fluctuations depends on two independent mechanisms, the variance of the fitness function and its slope. This can be proved analytically in the model of “Appendix C”, see Eq. 13. In order to understand which mechanism plays a major role in simulations, we performed an ordinary least square regression of \(1/\tau \), the inverse of the time to equilibrium, on \(\sigma _\mathrm{S}^2\), the variance of the fitness function around the equilibrium point, and \(\gamma \), the slope of the fitness function. The results of the regression are presented in Table 1.

Table 1 Results of the ordinary least square regression of \(1/\tau \) on \(\sigma _\mathrm{S}^2\) and \(\gamma \)
Full size table

The regression is very good, with a coefficient of determination of 0.95. Both variables impact negatively on the inverse of time to equilibrium. This is expected, since a higher variance of the fitness function should increase the time to equilibrium, as well as a higher slope (because the slope is negative). Finally, we can conclude that both mechanisms play an important role, but \(\sigma _\mathrm{S}^2\) explains most of the variations of the times to equilibrium, with a weight six times superior to the weight of the slope \(\gamma \). As a consequence, one cannot simply infer the dynamics from the static considerations made in Sect. 4. Finally, we show in Fig. 9 the results of the regression, plus a plot showing the variation of \(1/\tau \) against the expression found analytically (see Eq. 13 in “Appendix C”). Also in this case, the agreement is overall good, indicating that the SDE model captures quite well the dynamics. The better performance of the regression might be due to the fact that the analytical model is linearized around the equilibrium, whereas the system can in fact start quite far from it.

Fig. 9
figure 9

Inverse of time to equilibrium versus a combination of the standard deviation and the slope of the fitness function. Left: the weights are the results of an ordinary least square (\(R^2 = 0.95\)) regression. The solid red line shows the results of the regression. Right: the weights are coming from analytical arguments (Eq. 13 in “Appendix C”). The solid red line is a linear regression (\(R^2 = 0.67\)) (Color figure online)

Full size image

The conclusion of the regression is that the time to convergence is influenced by the stochastic behaviour of the fitness function as well as its general shape. In physical terms, it means that the air traffic system as idealized by this model can be quite far from the equilibrium, due to (i) inadequate policies (the slope) (ii) the general stochasticity of the system. Hence, policy makers should carefully assess if any changes in policy is likely to have an impact due to the level of randomness of the system.

6 Conclusions

In this paper, we have presented a stylized model of the allocation of flight plans. We used an agent-based model to simulate the behaviour of different air companies and the network manager. In the model, different types of air companies are competing for the best paths on the network of sectors and the best times of departure. Since the sectors are capacity constrained, in some high traffic conditions the companies might be forced to choose suboptimal flight plans, according to their strategies or cost function.

When different types of companies are competing on the same airspace, their relative satisfaction depends highly on the environment—the airspace, the waves of departure—but also the competition—the fractions of different types of companies. In a nutshell, we find that the companies are performing better when they are competing against other types of companies, in a mechanism of “niche” leading to behaviours similar to those of the minority game [5]. This conclusion is quite generic, since we have shown to hold both in a baseline model that can be treated analytically and in a more complex model, which we investigated with numerical simulations.

As a consequence, it is possible to reinterpret the model as an evolutionary game, through the use of the difference of satisfactions as a fitness function which sets the capacity of a type of company to expand its business by having more flights in the future. In this framework, the populations are the types of companies using the “rerouting” or “shifting” strategies. The study of the shape of the fitness function shows the existence of a stable equilibrium point for the mixing parameter for every values of parameters. Interestingly, this equilibrium point is distinct from the point where a global satisfaction is optimal for the system as a whole. This indicates that the system left alone will not converge to the global optimum, but to a different equilibrium point.

In order to study more in details, the real dynamics of the system around the point of equilibrium, we iterated the model with a reproduction rule mimicking the fact that higher satisfaction for an airline may be converted to better possibilities of expanding business. We found that the dynamical point of equilibrium is different from the one derived from the root of the fitness function (static equilibrium). This is a purely dynamical effect which is driving the point of equilibrium even further from the global optimum. Moreover, we found that both the convergence time to the equilibrium and the fluctuations are highly dependent not only on the slope of the fitness function, but also on its variance.

These results have several policy implications. On the one hand, the fact that the root of the fitness function is distinct from the global optimum means that the regulators should step in. Indeed, issuing well-designed regulations could help the system to have a point of equilibrium closer to the global optimum. On the other hand, the dynamical effects could blur the picture. Indeed, long time to convergence and high fluctuations combined with changing business conditions mean in practice that the system is always out of equilibrium. It is thus hard for the regulators to design incentives to drive the system towards the optimum. A more precise set-up and a more detailed calibration would be needed to definitely assert the potential consequences of regulations.

The model presented here is an idealized version of the reality, a simple, yet phenomenologically rich, toy model. It allows to catch some high-level, emergent, phenomena that are inaccessible to more complicated ones due to the large number of parameters. The existence of a point of equilibrium, its behaviour in certain environments, and its dynamics have certainly a scope broader than the present model. Moreover, the model is not really specific to the air traffic. In fact, it could be adapted to other situations, like packets propagation over the physical network of the internet, with minimal effort. As such, it can be viewed as a quite general model of transportation where entities need to send some material over a capacity-constrained network, thus competing for time and space.

Two possible directions for extension of the present work are the following. First, a more detailed modelling could allow to draw some more precise conclusions about the present and future scenario in ATM. A first path has been made in this direction with another version of the model [10], based on navigation points instead of sectors. The model is also coupled to a tactical part, allowing to simulate the conflict resolution of traffic controllers. The code for this model is freely available.Footnote 5 The second direction is towards model calibration. This is in general a challenging problem because data on strategic allocation are owned by companies and hardly available, especially when the details on many airlines are needed. A potential way to overcome this problem is through indirect calibration based on traffic data which contain the original flight plan and the last filled one. Data mining techniques could be useful to infer from these data unobservable parameters of our model and therefore to calibrate it.