Céline Demouge

Marcel Mongeau

Nicolas Couellan

Daniel Delahaye

Climate-aware air traffic flow management optimization via column generation

Keywords: Air Traffic Flow Management, Contrails, Column generation, Climate impact, Contrail metrics

Aviation is one of the global warming contributors. Its impact is due to CO 2 and non-CO 2 effects. Trajectory design is one of action levers for minimizing the environmental impact of air transportation. However, it affects the Air Traffic Management and should satisfy airspace constraints, especially airspace capacities. This paper proposes a climate-aware version of the Air Traffic Flow Management (ATFM), focusing on CO 2 and one particular non-CO 2 effect: condensation trails (contrails), although other non-CO 2 effects can be integrated. An ATFM optimization model is proposed, solved by a column generation approach. This problem is solved using different metrics, from simple to more complex and realistic ones. Numerical experiments are conducted both in the lateral case and when the cruise altitude becomes a decision variable. The impact of airspace capacities is also evaluated. The problem instances that are studied are built from realistic open-access data and made publicly available.

Introduction

Significant progress has been made over decades, but the aviation industry is still investing in further improvements in the efficiency of flight operations, particularly about to climate change. When assessing the impact of air transport, it is necessary to take into account both CO 2 and non-CO 2 effects. These non-CO 2 effects are estimated to be responsible for two-thirds of the sector's climate impact [START_REF] Lee | The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018[END_REF]. Among these effects, condensation trails (contrails) can be cited. They are one of the most concerning effects of air transport, but their impact is still the subject of many studies to improve confidence in its prediction [START_REF] Lee | The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018[END_REF].

Based on meteorological and atmospheric parameters, some geographical areas can be considered favorable for persistent contrails at some specific times. These areas should be avoided, but not at the expense of an excessive increase in CO 2 emissions. In addition, flight safety must always be guaranteed. For example, contrail avoidance should not create conflicts between aircraft.

Flows are organized to comply with airspace capacities, and to meet as much as possible the demand. Indeed, after the COVID-19 crisis, the traffic is growing again [START_REF] Airbus | Global market forecast[END_REF]. Air navigation service providers must organize the expected traffic safely. The goal is to minimize delays or cancellations while satisfying airspace capacity constraints and avoiding conflicts between aircraft (loss of vertical and/or horizontal separation). A sustainable development approach should focus on minimizing the environmental impact of air traffic. In that case, the areas that are avoided due to contrails will be less flown, while the surrounding areas are likely to be overloaded (because the detour would be minimized as much as possible). Flight-path computations are to be done well in advance of the flight (from several hours to several weeks before). In this long term perspective, safety is then guaranteed by satisfying airspace capacity constraints. In the shorter term, ensuring horizontal and vertical separations between aircraft is necessary to guarantee safety. Conflicts are not considered in the first case since there are too many sources of uncertainties a long time in advance.

Therefore, this paper proposes the following contributions:

• adaptation of the classical Air Traffic Flow Management problem to environmental, CO 2 and non-CO 2 , issues;

• proposal of a solution to take into account contrails at the level of the complete airspace network, trading off with CO 2 emissions;

• numerical experiments with realistic airspace and traffic data above France;

• numerical experiments considering various contrail metrics;

• software code associated with the proposed method and with the creation of the instances are made available.

The paper is organized as follows. Section 2 reviews the scientific literature on trajectory optimization for contrails mitigation and on air traffic flow management problem. Our mathematical optimization model is introduced in Section 3 together with details on contrail cost computation. Section 4 explains the resolution method. Section 5 presents some case studies and numerical experiments, it also discusses the impact of the different parameters of the problem. Finally, Section 6 details the conclusions of the study and presents some perspectives.

Literature review

This section proposes a literature review on (Subsection 2.1) trajectory optimization for contrails mitigation, and (Subsection 2.2) on the Air Traffic Flow Management problem.

Trajectory optimization for contrail mitigation

Trajectory optimization for contrail mitigation has become an important topic in the literature. Indeed, contrails are becoming a major issue in the computation of the impact of air transport. Different types of modeling have been proposed, with various points of view, constraints, and costs associated with contrails. Different resolution methods are also proposed.

The problem is addressed under the form of an optimal control problem in [START_REF] Sridhar | Aircraft trajectory optimization and contrails avoidance in the presence of winds[END_REF]. It minimizes a weighted sum of total flight time, fuel burned, and flight time in persistent-contrail areas under flight-mechanics constraints. Soler et al. [START_REF] Soler | Flight trajectory design in the presence of contrails: Application of a multiphase mixed-integer optimal control approach[END_REF] propose also an approach with flight-mechanics constraints, considering several possible flight levels, and trade-offs between passenger travel time, fuel, CO 2 emissions, and contrail cost. In both cases, persistent contrail areas are assumed to be supersaturated in ice. The resolution approach proposed in [START_REF] Soler | Flight trajectory design in the presence of contrails: Application of a multiphase mixed-integer optimal control approach[END_REF] is based on multiphase mixed-integer optimal control techniques.

In [START_REF] Campbell | Fuel-optimal trajectory generation for persistent contrail mitigation[END_REF], the problem is modeled as a Mixed-Integer Linear Program (MILP) with a receding horizon framework. The authors' analysis over 20 days shows that 58% of persistent contrails can be avoided with an increase of 0.48% in fuel consumption. Genetic algorithms are used in [START_REF] Yin | Impact on flight trajectory characteristics when avoiding the formation of persistent contrails for transatlantic flights[END_REF]. Based on an air traffic simulation framework, an annual analysis of transatlantic flights is done, showing the seasonal effects of contrails. The authors also show variability in trade-offs between total flight time and flight time in contrail areas. Indeed, when allowing an increase of 2% in total flight time, the reduction of the distance in contrails varies from 20% to more than 80%. Graph-adapted methods are also used, as in [START_REF] Rosenow | Multicriteria-optimized trajectories impacting today's air traffic density, efficiency, and environmental compatibility[END_REF] with the A* algorithm [START_REF] Hart | A formal basis for the heuristic determination of minimum cost paths[END_REF] in a simulation context to consider overall economic and environmental costs.

Interested readers can find in [START_REF] Simorgh | A comprehensive survey on climate optimal aircraft trajectory planning[END_REF] more details about optimizing trajectory taking into account contrails. This survey reveals that Air Traffic Management (ATM) considerations are not covered extensively in the literature. However, [START_REF] Baneshi | Conflict assessment and resolution of climate-optimal aircraft trajectories at network scale[END_REF] shows that avoiding conflicts has a significant impact on the environmental cost. Moreover, avoiding contrail zones means potentially saturating adjacent airspace, and this has a definitive impact on the air traffic controllers workload.

Air Traffic Flow Management problem

The Air Traffic Flow Management (ATFM) problem raised from the operational issues of maximizing the number of flights that can be accepted in a given airspace. The problem of scheduling flights to minimize congestion was first introduced in [START_REF] Amedeo | The flow management problem in air traffic control[END_REF]. Then, [START_REF] Bertsimas | The air traffic flow management problem with enroute capacities[END_REF] stated and modeled the ATFM problem. It is designed to direct air flows from one sector to another, from departure to arrival, with the possibility to act on several flight parameters, while satisfying air sector capacity constraints. The goal is then to minimize total delay incurred by the considered flights. Different approaches have been developed through the years, with different degrees of freedom and various levels of operational realism. Moreover, some studies consider only the passage from one sector to the next one, while other works take a closer look at the routes taken.

In [START_REF] Bertsimas | An integer optimization approach to large-scale air traffic flow management[END_REF], an integer programming approach is proposed, considering aircraft rerouting. It also includes speed modulation options and takes into account flight phases. This model is efficient enough to solve national-size airspace instances, as decisions for rerouting are made efficiently without auxiliary variables. Another way to obtain an efficient model is to use pre-processing, as done in [START_REF] Agustın | On air traffic flow management with rerouting. Part I: Deterministic case[END_REF]. In [START_REF] Weigang | Intelligent computing methods in air traffic flow management[END_REF], the air traffic system is represented by a multi-agent model and the ATFM problem is solved by specific methods.

Some other studies focus on combinatorial optimization techniques, adapted for large-scale practical problems. For instance, in [START_REF] Diao | A sequence model for air traffic flow management rerouting problem[END_REF][START_REF] Balakrishnan | Optimal large-scale air traffic flow management[END_REF][START_REF] Richard | Computing 4D near-optimal trajectories for dynamic air traffic flow management with column generation and branch-and-price[END_REF], Dantzig-Wolfe decomposition and column generation are used to solve the ATFM problem. Such approaches result in efficient resolution methods, taking into account the complete 4D (space and time) trajectories.

A related study is proposed in [START_REF] Hamdan | An environmental air traffic flow management model[END_REF] for CO 2 minimization. A previous study [START_REF] Demouge | A time-dependent subgraph-capacity model for multiple shortest paths and application to CO2/contrail-safe aircraft trajectories[END_REF] introduced a subgraph capacity-constrained multiple shortest path approach considering air-sectors capacities while mitigating contrail impact. In the mentioned paper, the method allows computing two-dimensional CO 2 /contrail-safe trajectories for a set of flights using time window approach for time dependence consideration and a simple metric for contrail consid-eration, the Global Warming Potential. This study introduces an ATFM model adapted to environmental issues, taking into account CO 2 and non-CO 2 effects. Here, we study the benefit of taking into consideration non-CO 2 phenomena, such as contrails in terms of airspace capacity, and then we propose a rearrangement of the air traffic.

Problem formulation

The context of the problem is presented in Subsection 3.1, the mathematical model is introduced in Subsection 3.2, and Subsection 3.3 focuses on environmental cost computation (contrail cost computation and comparison with CO 2 ).

Context

Let us consider an airspace such as that above a country or a set of countries. Aircraft fly in this airspace following a sequence of points (in 2D, or in 3D when considering altitude), linked by straight-line segment routes. These navigation points, or waypoints, define a set of vertices of a graph. The segment routes, defined from one waypoint to another, define the (directed) arcs of this graph (they generally cannot always be flown in both directions). This study takes place in the context of Free Route Airspace (FRA), applied in European countries upper airspace. FRA aims to replace the traditional air route system with one based on freely flying from a waypoint to another. Planning a trajectory then boils down to proposing a sequence of waypoints, similar to the former system, but without airway requirements. This new approach allows aircraft to use more direct routes, reducing thereby the distance flown and CO 2 emissions.

For the moment, rules are still established to regulate flights from waypoints to waypoints, but the aim is to relax such constraints. This study models these rules via the construction of a graph arc set. More precisely, the underlying graph G = (V, A), of our optimization process is built as follows. The vertex set V is the set of waypoints. An arc of the arc set A connects two waypoints if the distance between these two waypoints is less than some user-defined threshold distance, noted D, and greater than another user-defined distance, D. The set of arcs is then:

A = (u, v) | u ∈ V, v ∈ V, D ≤ d u,v ≤ D , (1)
where d u,v denotes the distance between nodes u and v. These minimum and maximum distances can be defined so as to vary according to airspace regions. For example, over areas without airports or above oceanic zones, the allowed distances between two neighboring nodes may be greater, and the graph be less dense.

The airspace is divided into unitary sectors combined into sectors to define areas of responsibility for air traffic controllers. The sectors' configuration is defined for a day or a few hours and have an associated capacity. This capacity limit ensures that air traffic controllers workload is acceptable to guarantee safety. Typically, these sectors can be dynamically reconfigured, by grouping or ungrouping them, but, in this study, for the sake of simplification, they are considered fixed.

The airspace is partitioned into sectors, and each waypoint belongs to one (and only one) sector:

V = s∈S s, ( 2
)
where S is the set of sectors. Moreover, it can be determined at any time in which sector an aircraft is located. At each given time period t ∈ T , each sector s has a capacity C s,t given in terms of number of aircraft. This capacity limit guarantees a certain level of safety at the strategic level. In practice, this limit determines whether air traffic controllers will be able to manage these aircraft in terms of workload. On a shorter timescale, one ensures safety by considering conflicts (loss of sufficient distance separation between aircraft). However, this study focuses on phase which involves several types of uncertainty and conflicts cannot efficiently be taken into account to guarantee safety, as such events are too tactical.

Our purpose is to direct each flight, over a time-window horizon, on the graph, from a starting point, to a destination, minimizing the total environmental impact while satisfying the sectors' capacity constraints at each time period.

Mathematical model

This subsection introduces our mathematical optimization model.

Let us first define the notation for the input data of the problem:

• F: the index set of flights (with corresponding departure and arrival points, departure time, and aircraft type);

• S: the index set of sectors;

• T : the index set of time periods, that are typically defined with a constant time step;

• C s,t : the capacity of sector s at time period t, s ∈ S, t ∈ T ;

• G = (V, A): the graph of waypoints (vertices) and free routes (arcs) on which aircraft can be directed.

We shall use some further notation:

• P f : the index set of the feasible paths for flight f ∈ F;

• P s,t : the index set of the paths passing through sector s ∈ S during time period t ∈ T ;

• ρ f,p : the (environmental) cost for flight f ∈ F to fly path p ∈ P f . Details on environmental cost computation is given in Subsection 3.3.

The problem aims at finding a feasible path for each flight minimizing the total environmental impact, while satisfying the capacity of sectors. Then, it can be written as follows:

minimize Y f ∈F p∈P f ρ f,p y f,p (IMP)
subject to

p∈P f y f,p = 1, f ∈ F, (3a)
f ∈F p∈Ps,t∩P f y f,p ≤ C s,t , s ∈ S, t ∈ T , (3b)
y f,p ∈ {0, 1}, f ∈ F, p ∈ P f , (3c)
where y f,p is equal to 1 if path p is chosen for flight f . Constraint (3a) ensures that one and only one path is chosen for each flight. Constraint (3b) stipulates that the capacity of the sectors must be satisfied. It ensures that, for each time period, the number of aircraft flying through each sector during this time period must not exceed the defined capacity. Indeed, the set P s,t ∩ P f represents the set of paths candidate for flight f crossing the sector s during the time period t. The paths considered can be in 2D or 3D, depending on whether the altitude is taken into account as a decision. The paths can also be time dependent. These different considerations do not affect the proposed model as they simply involve larger sets of possible paths to be considered.

Environmental cost computation

This subsection focuses on the computation of environmental cost of flying a path for a given flight. This cost takes into account CO 2 and non-CO 2 effects. A bi-objective approach could then have been implemented. Here, however, it is desirable to reduce the total environmental impact, and so we have implemented a single objective approach by summing the CO 2 and non-CO 2 effects. Then ρ f,p can be built as

ρ f,p = ρ CO 2 f,p + ρ non-CO 2 f,p
. As only contrails are considered in this study as non-CO 2 effect, ρ f,p = ρ CO 2 f,p +ρ contrails f,p . A path p is then a sequence of N p chosen nodes of the graph: p = (u 1 , u 2 , . . . , u Np ), where u i ∈ V, i = 1, 2, . . . , N p , respectively flown at times (t 1 , . . . , t Np ). Then the cost of this path is the sum of cost of flying this sequence of arcs. The order of arcs is important since it has an impact of the time at which the arc is flown and then on the associated cost. Then, the cost of a path can be written as follows:

ρ f,p = Np-1 i=1 w CO 2 f,(u i ,u i+1 ),(t i ,t i+1 ) + w contrails f,(u i ,u i+1 ),(t i ,t i+1 ) , (4)
where [t i , t i+1 ] is the time interval over which arc (u i , u i+1 ) ∈ A is flown and w CO 2 f,(u i ,u i+1 ),(t i ,t i+1 ) (respectively w contrails f,(u i ,u i+1 ),(t i ,t i+1 ) ) is the cost for flight f to fly arc (u i , u i+1 ) over [t i , t i+1 ] associated to CO 2 (respectively to contrails).

The time interval over which an arc is flown is only determined by the sequence of nodes since aircraft's speeds are not part of the decision variables in this study. Then, we can rewrite:

ρ f,p = Np-1 i=1 w CO 2 f,(u i ,u i+1 ),t i + w contrails f,(u i ,u i+1 ),t i , (5)
since the time interval over an arc can be determined thanks to the arrival time on the first node. In the sequel, all the cost are written as a function of t i according to this remark. First, balance between flight time and flight time in contrails areas is presented before presenting the Global Warming Potential metric and the Algorithmic Climate Change functions. In the sequel, persistent contrail areas are considered to be ice supersaturated areas, determined thanks to weather data.

Flight time balance with flight time in contrail areas

A straightforward approach is to use a weighting sum of the flight time (fuel consumption) and the flight time in persistent contrail areas (named contrail areas in the sequel). Then, for a given flight f and a given path p = (u 1 , . . . , u Np ) ∈ P f , we have:

w CO 2 f,(u i ,u i+1 ),t i = (1 -α) • δt f,(u i ,u i+1 ),t i (6)
w contrails f,(u i ,u i+1 ),t i = α • δt contrails f,(u i ,u i+1 ),t i , (7)
and thus

ρ f,p = (1 -α) • ∆t f,p + α • ∆t contrails f,p , (8)
where δt f,(u i ,u i+1 ),t i (respectively δt contrails f,(u i ,u i+1 ),t i ) is the flight time (respectively flight time in contrail areas) for flight f over arc (u i , u i+1 ) over time interval [t i , t i+1 ], ∆t f,p (respectively ∆t contrails f,p ) is the total flight time (respectively total flight time in contrail areas) for flight f over path p, and α ∈ [0, 1( a user-defined weighting parameter.

Global Warming Potential

Various metrics have been developed to compare the impact of contrails with that of CO 2 . This paper presents two main ones: the Global Warming Potential (GWP) [START_REF] Fuglestvedt | Transport impacts on atmosphere and climate: Metrics[END_REF] and the algorithmic Climate Change Functions [START_REF] Van Manen | Algorithmic climate change functions for the use in eco-efficient flight planning[END_REF].

The GWP metric is a multiplicative factor that compares the impact of various greenhouse gases to CO 2 . This metric depends on a time horizon over which the impact is computed, and it is in general given per kg of greenhouse effect gas. For the particular case of contrails which we consider here, it is given on a per kg of CO 2 emitted in contrail areas basis. More details about GWP can be found in [START_REF] Fuglestvedt | Transport impacts on atmosphere and climate: Metrics[END_REF]. Table 1 gives different values of GWP for contrails for various time horizons, H. In the sequel, the so-called contrail-induced cirrus (CIC) GWP will be noted g H . Depending on the time horizon, H, considered, contrails have a more or less weight. Notably, they have less relative impact compared to CO 2 in the long term. The choice of time horizon in a study generally depends on the lifetime of the gas, political objectives, or comparability purposes (standard horizon). In the following, only contrail induced cirrus clouds will be considered, since they have the largest impact on the climate. However, other contrails can easily be included in our model by a simple adjustment of the cost function. The cost of a path is therefore computed as follows:

w CO 2 f,(u i ,u i+1 ),t i = c f,(u i ,u i+1 ),t i ,p i (9)
w contrails f,(u i ,u i+1 ),t i = λ (u i ,u i+1 ),t i • g H • c f,(u i ,u i+1 ),t i ,p i (10)
ρ f,p = (1 + λ p • g H ) • c f,p , (11)
where c f,(u i ,u i+1 ),t i ,p i (respectively c f,p ) is the quantity of CO 2 emitted by flight f when flying arc (u i , u i+1 ) over time interval [t i , t i+1 ] knowing that the path flown until u i is p i = (u 1 , . . . , u i ) (respectively when flying path p), λ

(u i ,u i+1 ),t i ∈ [0, 1] (respectively λ p ∈ [0, 1]
) is the part of arc (u i , u i+1 ) over time interval [t i , t i+1 ] (respectively path p) in contrail areas. Knowing the path flown until u i is essential for CO 2 computation since it depends on the aircraft mass which depends on the path flown until the considered point.

Remark that the quantity of CO 2 emitted by flight f when flying path p is directly proportional to the fuel consumed by flight f when flying path p.

Then, without loss of generality, c f,p designates in the sequel the quantity of fuel consumed by f on p.

Algorithmic Climate Change Functions

Algorithmic Climate Change Functions (aCCFs) [START_REF] Van Manen | Algorithmic climate change functions for the use in eco-efficient flight planning[END_REF][START_REF] Yamashita | Newly developed aircraft routing options for air traffic simulation in the chemistry-climate model EMAC 2.53: AirTraf 2.0[END_REF] are another possible metric of the impact of contrails relative to that of CO 2 . They are based on Average Temperature Response (ATR), which reflects the temperature change over a time horizon due to a given phenomenon. It is a cumulative measure on the long-term effect on the temperature. In the sequel, the time horizon dependence is omitted for clarity of notations. ACCFs exist for different effects, and this study focuses in particular on CO 2 and contrails. Concerning CO 2 , the value of the associated aCCF is a constant, as it is considered that the CO 2 emissions have the same impact over the world and at any time. It is expressed in kelvin (K) per kg of fuel. Concerning contrails, the value of the associated aCCF depends on the location and the time of contrail creation, and is expressed in K per km flown. The aCCF values are computed on a space-time grid according to weather and radiative data. Indeed, contrail existence depends on weather data, and its impact depends on weather and solar radiative data. Remark that aCCF values associated with contrails can be negative (cooling contrails) or positive (warming contrails), depending on location and time.

When considering aCCF metrics, the cost of path p for flight f becomes:

w CO 2 f,(u i ,u i+1 ),t i = accf CO 2 • c f,(u i ,u i+1 ),p i ,t i (12)
w contrails f,(u i ,u i+1 ),t i = accf contrails ((u i , u i+1 ), t i ) • d u i ,u i+1 (13)
ρ f,p = accf CO 2 • c f,p + Np-1 i=1 accf contrails ((u i , u i+1 ), t i ) • d u i ,u i+1 , ( 14
)
where accf CO 2 is the constant aCCF value for CO 2 , accf contrails ((u i , u i+1 ), t i ) is the aCCF value for contrails on arc (u i , u i+1 ) over time interval [t i , t i+1 ], d u i ,u i+1 is the distance between nodes u i and u i+1 .

Resolution approach

To favor fairness between aircraft, it is not desirable to proceed sequentially computing the best trajectory, aircraft by aircraft, and updating the sectors' capacity as computations are done. The trajectories must rather be simultaneously computed. For this, considering the literature about the ATFM problem [START_REF] Balakrishnan | Optimal large-scale air traffic flow management[END_REF][START_REF] Richard | Computing 4D near-optimal trajectories for dynamic air traffic flow management with column generation and branch-and-price[END_REF], we solve this problem using to the column generation technique, which is efficient for addressing problems involving a large number of variables and parsimonious solutions (numerous variables are put to zero at the optimum) [START_REF] Desaulniers | Column generation[END_REF]. This section presents the general process of column generation, and then it details how it is adapted to the climate-aware application proposed in this paper.

Figure 1 presents the general column generation algorithm. In the sequel, the step numbers are those used in the figure .  For solving the Integer Main Problem (IMP), the linear relaxation, (LMP), is solved. A set of decision variables is selected to form a feasible solution (Step 1). Then, a restricted version of the problem (LMP) is solved by using these variables (Step 2). Using the dual variables associated with the solution found to this problem (RMP), subproblems are solved to find new variables with negative reduced cost, expecting to be improving variables (Steps 3 to 5). Then, (RMP) is solved with these new variables added to the pool of solutions (Step 2). This loop is repeated until no new variable is found (

p∈P f y f,p = 1, f ∈ F, (15a)
f ∈F p∈Ps,t∩P f y f,p ≤ C s,t , s ∈ S, t ∈ T , (15b)
y f,p ≥ 0, p ∈ P f , f ∈ F, (15c)
where y f,p becomes a continuous variable whose upper bound 1 is ensured by constraints (15a).

During the process, at each iteration, one considers the following restricted main problem version (RMP) continuous relaxation of our problem. It is solved with only a subset of each of the path sets P f , P r f (Step 2). The restricted problem is initialized using big-M techniques [START_REF] Desaulniers | Column generation[END_REF], that is to say, it is initialized with high-cost solutions that do not consume sector capacity (Step 1). The rest of the process consists of adding interesting variables to the restricted model, without adding all possible variables, but still having a proof of optimality.

minimize (y

f,p ) f ∈F ,p∈P r f f ∈F p∈P r f ρ f,p y f,p (RMP) subject to p∈P r f y f,p = 1, f ∈ F, (16a)
f ∈F p∈P r s,t ∩P f y f,p ≤ C s,t , s ∈ S, t ∈ T , (16b)
y f,p ≥ 0, p ∈ P f , f ∈ F, (16c)
where P r s,t is the restricted-set analogue of P s,t and P r s,t ⊂ P s,t . New variables are added thanks to the dual values obtained when solving the problem (RMP) (Step 2). To obtain these new variables, the following problem is solved for each considered flight f ∈ F(Step 3):

c f * = minimize X (u,v)∈A t∈T w f,(u,v),t x f u,v,t -µ f - s∈S t∈T π s,t r f s,t
subject to flow conservation constraints

(SP f )
where

• A is the set of arcs of the graph (free routes on which aircraft can be directed);

• w f,(u,v),t is the cost to fly (u, v) ∈ A at time slot t ∈ T (as defined in Equation ( 5));

• µ f is the dual variable associated to constraint (16a), f ∈ F;

• π s,t is the dual variable associated to constraint (16b), s ∈ S, t ∈ T ;

• r f s,t is an auxiliary variable that indicates whether sector s is flown by the considered aircraft during time slot t, s ∈ S, t ∈ T ;

• x f u,v,t is the time-dependent shortest-path variable corresponding to arc (u, v) ∈ A and time slot t ∈ T .

The time at which an arc is flown is only determined by the sequence of nodes since aircraft's speeds are not part of decision variables in this study.

The subproblem (SP f ) is solved following a process similar to that of [START_REF] Balakrishnan | Optimal large-scale air traffic flow management[END_REF]: it boils down to finding a shortest path in a space-time network defined for each flight

G f t = (V f t , A f t ). It is a time expansion of graph G = (V, A)
, where the set of space-time vertices is defined as V f t = {(v, t) | v ∈ V and t ∈ T }, and A t is the set of time-expanded arcs such that a space-time node (u, t u ) is linked to another space-time node (v, t v ) if (u, v) ∈ A and t v -t u is between the minimum and the maximum flight time for the considered aircraft from u to v. Finally, the cost w

′ f,ut,vt of an arc (u t = (u, t u ), v t = (v, t v )) is defined as: w ′ f,(ut,vt) = w f,(u,v),tu - tv-1 t=tu π s (u,v),t , t , (17)
where s (u,v), t is the sector flown by the aircraft considered at time t on the arc (u, v).

A new decision variable for the restricted problem is found for the flight f ∈ F if the shortest path on the time-expanded network has a cost lower than µ f , and the process is done for all the flights. Thanks to powerful techniques such as dynamic programming on time-expanded networks, we shall see that the subproblem can be solved in a reasonable computation time. Since it must be solved for all the flights at each iteration, it can be directly parallelized.

When for some flight f ∈ F a solution is found with c f * < 0, then it is added to the restricted set of solutions in the problem RMP (Steps 4 and 5). Otherwise, no variable is candidate to be added in the main restricted problem, and the process is stopped. The problem (LMP) is then solved as an integer problem (Step 6).

Numerical experiments

This section presents some case studies considering the French airspace. First, Subsection 5.1 describes the data used. Three instances are built to illustrate the outcomes. The results obtained with different contrail metrics are compared in Subsection 3. Then, Subsection 5.3 shows the impact of capacity. Finally, Subsection 5.4 discusses the impact of altitude.

Data description

This subsection provides details about the data used for numerical experiments. First, the aircraft data are detailed before focusing on the traffic and airspace data used to build realistic instances. Then, weather data are detailed and finally, the three illustrative instances used in the sequel are detailed.

Aircraft data

Aircraft data such as airspeed and mass are extracted from the SKYbrary website [START_REF]SKYbrary Aviation Safety[END_REF]. Fuel flow and emissions data are then extracted from the Ope-nAP python library [START_REF] Sun | Openap: An opensource aircraft performance model for air transportation studies and simulations[END_REF]. OpenAP has been used in this study since it is an open-source library with accuracy in line with this study purposes.

Traffic and airspace data

Traffic data are extracted from Eurocontrol R&D data [START_REF]R&D data[END_REF]. A day and a time period are chosen, and only cruising points (from Flight Level (FL) 300 to FL 400) over France are kept. To obtain starting and destination points, we select the vertices of the graph that are closest from the first and the last points of the trajectory. The time of passage over the starting point is considered to be the departure time of the flight concerned. Only long enough trajectories are kept (those with cruising phase longer than 150 Nautical Miles (NM)). The code to process these data and the code associated with this paper are made publicly available at [START_REF] Demouge | Climate aware ATFM[END_REF].

French airspace data are extracted from the French aeronautical information (SIA) website [START_REF]Service de l'Information Aéronautique de la Direction Générale de l'Aviation Civile[END_REF]. Especially, upper airspace waypoints data and unitary sectors are extracted. This is considered as the sector configuration for the experiments. Removing waypoints from some areas that were too dense in points yields the vertex set displayed on Figure 2. The vertex set is then partitioned into sectors as shown by Figure 3. Other airspace instances, for other European countries and considering several countries, are given with the code associated with this paper [START_REF] Demouge | Climate aware ATFM[END_REF]. These instances are not based on real navigation databases but built ran-domly. Points and sectors are drawn randomly over the countries considered. Their number and density can be adjusted through the code provided.

Weather and contrails data

Weather data are extracted from the ERA5 reanalysis dataset [START_REF] Hersbach | ERA5 hourly data on pressure levels from 1940 to present[END_REF]. They are directly used for wind, and processed thanks to the CLIMaCCF library [START_REF] Dietmüller | A python library for computing individual and merged non-CO 2 algorithmic climate change functions[END_REF] to extract contrail areas and aCCFs values.

Then, aCCFs and contrail areas data are processed to evaluate the aCCF above each graph arc, and to evaluate the percentage of each arc that lies in contrail areas. The code to process these data is also given with the code associated with this paper [START_REF] Demouge | Climate aware ATFM[END_REF].

Instances description

Three instances are presented; they are named FRA H1 M1, FRA H2 M2, and FRA H3 M3. The arcs of the underlying graph of these instances are constructed using the distance bounds D = 40 NM, and D = 130 NM and the time discretization is based on a 5-minute time step. The initial mass of each aircraft is set to 80% of its Maximum Take Off Weight (MTOW) as this data was not available. Table 2 details the instances, in particular the period used for weather and traffic data. Instance FRA H1 M1 will be used in Subsection 5.4 for illustrating altitude consideration, whereas instances FRA H2 M2 and FRA H3 M3 will both be used for illustrating the impact of capacity (Subsection 5.3) and metrics (Subsection 5.2).

Comparison of contrail metrics

This subsection compares the results obtained when using different contrail metrics. The last two instances, FRA H2 M2 and FRA H3 M3, are solved at cruising altitude (FL310). Figures 4 and5 display the contrail geographical setup (aCCFs and persistent contrail areas) for both instances. Two capacity scenarios are considered for each instance: an unconstrained-capacity scenario, noted C 0 , and a constrained-capacity scenario, for which the capacity of each sector s, at each time slot t is set to C s,t = 15. The four metrics studied are:

• No contrail consideration: only the CO 2 is minimized;

• GWP with a 20-year time horizon (noted GWP20);

• GWP with a 100-year time horizon (noted GWP20);

• aCCFs. The presented results are obtained using the Java API of CPLEX [START_REF]Ibm ilog cplex optimization studio[END_REF] to solve the main optimization problem and Java coding for the subproblem resolution via dynamic programming. Computations are performed on a personal computer with an Intel Core i5-10210U, 1.60 Hz, with 8 Go RAM, and a Debian Linux OS. The mean computation time for the presented instances is 476 seconds.

Figure 6 shows the results obtained, on CO 2 emissions and flight time in contrails, in total (summed over all the flights) in comparison to the case where no contrails are considered. Figures 7 and8 show the same indicators in the form of distribution on all flights. Comparisons are done with the case of no contrail consideration, with the same capacity scenario. Sensitivity of the results with respect to the sector is studied in the next subsection. One first remarks that the GWP metric, which does not consider contrails with cooling effects, yields a strategy that is different from that one obtained with aCCFs, which do consider cooling effect (as observed for instance FRA H3 M3). Indeed, aCCFs are more complex metrics, that yields different results in some cases, contrail areas can become attractive, or it is less imperative to avoid them. In some cases, GWP induces a high increase in CO 2 emissions. Nevertheless, one must keep in mind that the presented results are always optimal with respect to the metric chosen. Post-processing can be used to avoid undesirable situations, or a constraint can be added to the model in order to impose a maximum acceptable increase in CO 2 emissions.

Impact of capacity

Capacity constraints can reduce the ability of the traffic to avoid contrail areas. In the last figures presented, results are compared for the same capacity scenario. Then, the impact of capacity is not directly studied. This subsection focuses rather on comparing capacity scenarios with the same metric, here aCCFs.

Figure 9 displays the evolution of total aCCFs (CO 2 and contrails) when solving instance FRA H2 M2. Again, these experiments are done in the case where for each sector s ∈ S and each time slot t ∈ T , the capacity has the same value C s,t = C.

Figure 9 shows that when the capacity decreases, the total environmental cost increases. However, the capacity constraints are hard constraints that must be satisfied. It is therefore crucial to carry out the optimization proposed in this paper to mitigate the impact of capacity constraints on the Another observation is that the computation time increases when the capacity is reduced, since the number of iterations increases. In these experiments, the computation time varies from 432 to 521 seconds.

Adding the altitude in decisions

As observed in [START_REF] Williams | Reducing the climate change impacts of aviation by restricting cruise altitudes[END_REF], changing altitude can be an efficient complementary mean of contrail mitigation. This subsection focuses on adding the choice of the cruising altitude as a decision variable for each aircraft in the optimization problem. First, the principle is explained before presenting results on instance FRA H1 M1.

General principle

As mentioned when modeling the problem (see Section 3), the proposed model can be used for 2D or 3D paths, or even in 4D when considering time dependency. Moreover, the resolution approach presented in Section 4 can also be used in both cases, as solely the arc set of the underlying graph is impacted. The sequel describes the instance of graph used in this study.

The choice of cruising altitude can be integrated into the presented problem by introducing artificial departure and ending vertices. For the sake of simplification, the aircraft are assumed to fly by default at the highest possible altitude, and that the only possible decision consists in decreasing the altitude of one or two allowed flight levels compared to the default flight level. Figure 10 illustrates the new arcs that are added to the 2D graph to take into account the possible altitude decisions (here for one aircraft) in the case where only one other flight level is allowed. Remark that the optimization model takes into consideration the fact that the airspeed and fuel flow vary with of the altitude. FL 0 is the default flight level, FL -is the first available flight level below FL 0 , s ′ i (respectively e ′ i ) is an artificial starting (respectively destination) node for including the choice of level and graph between s 0 i and e 0 i (respectively s - i and e - i ) is the 2D graph considering the cost associated with FL 0 (respectively FL -).

v 0 i v - i f 0 i f - i 2D-
The choice of this 3D graph instance is driven by several operational considerations. First, the choice of cruise altitude is driven by fuel consumption concerns, and turbulence or weather hazard avoidance. Moreover, as mentioned before, the cruise altitude is generally chosen as high as possible. It is a sensitive parameter which is difficult to modify by a large amount, and it is nearly impossible to fly at a higher altitude. Then, since the point of view in this study is the one of air traffic control, with no access to critical data such as the current aircraft mass, aircraft performances are not known well enough for the air traffic controller to issue large flight-level change instructions.

Results on instance FRA H1 M1

The sequel of this subsection shows results when considering two possible flight levels below the default flight level, for instance FRA H1 M1, where the default flight level is FL 320. The graph instance is built similarly as the one shown by Figure 10, but with two possible flight levels under the default one. Following common practice [START_REF]True airspeed[END_REF], the airspeed is assumed to decrease by 2% when the altitude decreases by 1000 ft. Figure 11 displays the contrail setup for this numerical experiment. It presents a day with a lot of contrails, and with variability according to the altitude, especially in the south-east of France. The optimization results are illustrated in figure 12 which shows, according to each metric considered, the proportion of flights on each flight level, in the capacity scenario C 1 . These results show that according to different metrics, it can be chosen to fly lower to avoid contrails. This choice can be done because avoiding contrail areas without altitude change can be too CO 2 -consuming. Flying higher is, in general, less CO 2 -consuming. However, if the cost of avoiding contrails areas is not too high, flying higher can be judicious. Moreover, in some cases, it can be optimal to fly within a contrail area at the highest (default) flight level if this contrail area involves a moderate cost and if it is not worth avoiding it, by flying lower or laterally.

Figure 12: Proportion of total simulation flights for each cruise flight level selected as a result of the optimization process, for instance FRA H1 M1 in capacity scenario C 1 . For each metric considered, the proportion of flights at FL300, 310 and 320 is shown in this order, from left to right. This instance presents small differences in contrail distributions between FL 320 and 310 but a significant reduction in south-east France for FL 300. Then, flight level 310 is less chosen, because the reduction of contrails is relatively small, but the fuel consumption is higher at this level than it is at FL 320. Once again, there are differences in the proportion of flights at each level depending on the metric. For instance, using GWP100 as an optimization criterion yields a solution with flight level 310 for one flight. Finally, the chosen instance studied here is capacity constraining. As a consequence, the choice of flying lower can also be made when the lateral avoidance is likely to violate capacity constraints.

Conclusions and future works

After presenting a literature review on trajectory computation for contrail mitigation and the Air Traffic Flow Management problem, this paper proposes to combine both issues into an original mathematical optimization model. We solve it with methods for the ATFM problem adapted to the case of specific environmental costs. In particular, this paper considers objective functions taking into account both CO 2 and non-CO 2 effects, and more specifically contrails. This is crucial since avoiding contrail areas has an impact on the geographic traffic distribution. For example, when favorable areas are avoided in a zone, the corresponding capacity is not used, which yields neighboring areas that are likely to be more congested. Since this paper focuses on the air traffic control point of view, such capacity issues are critical.

The proposed study relies on realistic instances and is applied to traffic instances, extracted from historical flight databases, that are made publicly available. The considered objective functions are based on different metrics for contrails, from simple ones to more complex ones. We showed that results are very sensitive to the choice of the particular metric. However, the presented optimization methodology can be easily adapted to any such metric. Moreover, this study is limited to the phenomenon of contrails for non-CO 2 effects, but other non-CO 2 effects could also be integrated into the proposed model, using an appropriate metric. For example, aCCFs [START_REF] Van Manen | Algorithmic climate change functions for the use in eco-efficient flight planning[END_REF] related to other effects can be considered. Finally, the case where the cruising altitude can be chosen is studied, showing that in some cases, flying lower is worth to minimize the total environmental impact.

The computational experiments conducted in this study is limited to the national geographical scale, but an extension to the European scale could be envisaged. This larger geographical scale means that more flights can be considered in their entirety, reducing their environmental costs.

The range of possible decisions in the proposed optimization problem could also be expanded. The proposed study only considers choosing one cruising altitude per flight. Avoiding contrails by lowering the reference altitude during the cruise phase can be efficient [START_REF] Williams | Reducing the climate change impacts of aviation by restricting cruise altitudes[END_REF], but there are hard constraints on aircraft cruise altitude, so the way to model such constraints will be crucial. In addition, here the configuration of the air sectors is fixed, but it can be optimized to match better the contrail distribution of the day or period under consideration. In particular, sectors can be grouped or ungrouped, as it is done in dense traffic areas, in order to increase capacity in the sectors that are adjacent to contrail areas, and to decrease it in the contrail areas. Finally, given the uncertainties that are associated with weather forecasts and the difficulty of predicting contrails [START_REF] Gierens | How well can persistent contrails be predicted?[END_REF], a promising future track of research is to design an optimization model that takes such uncertainty into account.

Figure 2 :

 2 Figure 2: Points used for numerical experiments.

Figure 3 :

 3 Figure 3: The vertices used for numerical experiments, partitioned into sectors, represented by different colors.

Figure 4 :

 4 Figure 4: Contrail setup in instance FRA H2 M2 at FL 310.

Figure 5 :

 5 Figure 5: Contrail setup in instance FRA H3 M3 at FL 310.

Figure 6 :

 6 Figure 6: Results obtained summing over all flights with the four considered metrics. Comparisons are relative to the case where no contrails are considered (CO 2 minimization) and in the same capacity scenario.

Figure 7 :

 7 Figure 7: Distribution over the flights of the results obtained when solving the FRA H2 M2 instance with the four metrics considered. These results are obtained when comparing to the case where no contrails are considered (CO 2 minimization) and in the same capacity scenario.

Figure 8 :

 8 Figure 8: Distribution over the flights of the results obtained when solving the FRA H3 M3 instance with the four metrics considered. These results are obtained when comparing to the case where no contrails are considered (CO 2 minimization) and in the same capacity scenario.

Figure 9 :

 9 Figure 9: Evolution of total flight time and flight time in contrail areas according to different values of the sector capacity C.

Figure 10 :

 10 Figure10: New arcs for the 3D graph to take into account possible altitude decisions for one aircraft i ∈ F when only one other flight level is allowed. FL 0 is the default flight level, FL -is the first available flight level below FL 0 , s ′ i (respectively e ′ i ) is an artificial starting (respectively destination) node for including the choice of level and graph between s 0 i and e 0 i (respectively s - i and e - i ) is the 2D graph considering the cost associated with FL 0 (respectively FL -).

Figure 11 :

 11 Figure 11: Contrail setup for instance FRA H1 M1, at FL 300, 310, and 320. The top of the figure shows the aCCFs values and the bottom shows the persistent contrail areas for each flight level.

Table 1 :

 1 Global Warming Potential for contrails for various time horizons, H, given on a per kg of CO 2 emitted basis.

		H = 20 years H = 100 Years H = 500 years
	GWP contrail (H) 0.74	0.21	0.064
	GWP CIC (H)	2.2	0.63	0.19

Table 2 :

 2 Instances summary.

	Instance	Date for weather data	Date for traffic data
	name
	FRA H1 M1 January 26, 2019, 14:00 UTC March 5, 2019, 22:00 UTC (1
			hour, 165 flights)
	FRA H2 M2 January 17, 2019, 23:00 UTC March 5, 2019, 15:00 UTC (1
			hour, 474 flights)
	FRA H3 M3 April 25, 2019, 09:00 UTC	March 5, 2019, 12:00 UTC (1
			hour, 518 flights)

Acknowledgements

The authors thank DGAC (French civil aviation authority) for prompting and funding this work, and more specifically its DTA and DSNA services.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.