domain: info math

text:
Céline Demouge
email: celine.demouge@enac.fr

Marcel Mongeau
email: mongeau@recherche.enac.fr

Nicolas Couellan
email: nicolas.couellan@recherche.enac.fr

Daniel Delahaye
email: delahaye@recherche.enac.fr

A time-dependent subgraph-capacity model for multiple shortest paths and application to CO2/contrail-safe aircraft trajectories

Keywords: Subgraph constrained time-dependent shortest path problem, Contrails, Air transportation, Mixed Integer Linear Programming

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A time-dependent subgraph-capacity model for multiple shortest paths and application to CO 2 /contrail-safe aircraft trajectories Transportation, and in particular air transportation, offers an important source of operations research problems, such as shortest path under constraints. Motivated by the question of more eco-responsible transport and in a sustainable development approach, these problems generally focus on the minimization of environmental costs.

In particular, for air transport, driven by a global will of the sector, the issues of green operations are gaining importance. The air transportation environmental impact is not solely due to CO 2 effects. The non-CO 2 effects are other greenhouse effect gases such as nitrogen oxide (NO x ) or other more complex phenomena like condensation trails (named contrails hereafter). Contrails form white trails at the back of the aircraft engines under certain conditions of humidity and temperature. They can disappear very quickly or on the contrary persist and turn into clouds, called cirrus [START_REF] Kärcher | Formation and radiative forcing of contrail cirrus[END_REF]. These can cause a parasol effect, usually during the day, preventing the solar rays from reaching the earth. On the contrary, they can also create a greenhouse effect, by reflecting the rays emitted by the earth. In total, despite the important uncertainties that remain, contrails have a negative effect on the climate [START_REF] Lee | The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018[END_REF].

Avoiding contrails therefore requires aircraft to avoid certain areas that are too humid and/or too cold. However, this should not be done at the expense of flight safety. This study proposes to plan flights several hours in advance, on a strategic scale. At this scale, ensuring safety means satisfying airspace capacity constraints. Indeed, the airspace is divided into sectors under the responsibility of one or several air traffic controllers. Thus, by limiting the number of aircraft per sector per time period, it is understood that the air traffic controllers workload will be acceptable to ensure safety. On a shorter time scale, a few minutes before arriving in an airspace sector for example, the air traffic controllers ensure safety by enforcing vertical and horizontal separations between aircraft.

To satisfy airspace capacity constraints, two approaches are possible: an airlinecentered or a network-centered approach. On the one hand, flight operators propose several flight plans and the network manager1 looks for a valid and best potential combination. On the other hand, the computations are done directly by the network manager, assimilated in the sequel to the air traffic control. This network-centered approach allows one to reach an optimum, to avoid more cases without a solution satisfying the capacity constraints and to evaluate the consequences of politics to avoid the negative impact of contrails. This approach is the one chosen in this paper.

The considered problem is close to the classical operations research issue of finding shortest paths, which includes several types of problems. A first example is the wellknown standard shortest path problems which can be solved by efficient polynomialtime algorithms. There are several extensions of this classical problem, for instance, the time-dependent shortest-path problem with dynamic costs and constraints. Variants involving constraints are NP-hard problems even when dealing with a single vehicle. Other problems deal with several vehicles simultaneously. The problem under study is

Previous related works

This section presents previous related works by focusing first, in Subsection 2.1, on individual route optimization for contrails mitigation. Then, in Subsection 2.2, different works done on assigning routes for several vehicles are presented.

Individual route optimization for contrail mitigation

In recent years, driven by various initiatives, the issue of green aviation has become more prominent in the literature. For instance, finding trajectories with the least possible CO 2 emissions is a topic well represented in the literature in particular by proposing studies calculating the optimal wind trajectories (see for instance [START_REF] Girardet | Wind-optimal path planning: Application to aircraft trajectories[END_REF][START_REF] Ng | Optimizing aircraft trajectories with multiple cruise altitudes in the presence of winds[END_REF][START_REF] Legrand | Robust aircraft optimal trajectory in the presence of wind[END_REF]). In addition, non-CO 2 effects are a particularly important topic. Several methods taking different points of view and with different resolution strategies have been developed.

To solve the problem in the most general case, optimal control methods have been implemented. Since the constraints related to the mechanics of flight are enforced, the computed trajectory is flyable in the free space. The methods chosen differ according to the dimension of the instance addressed, the objective function considered, and the number of aircraft involved. In [START_REF] Ng | Cross-polar aircraft trajectory optimization and the potential climate impact[END_REF], the problem is solved in 2D, thanks to optimal control by minimizing an objective function that takes into account contrail avoidance, fuel and flight time. It has been used on one-aircraft instances but also for trajectories between 12 city pairs. Hartjes et al. [START_REF] Hartjes | Contrail mitigation through 3D aircraft trajectory optimization[END_REF] solve the problem in 3D for a single aircraft, also by minimizing fuel, time and time in contrail areas. Some other papers present methods taking into account time, such as [START_REF] Matthes | A concept for multi-criteria environmental assessment of aircraft trajectories[END_REF].

Other studies rely on metaheuristics to solve the problem. For instance, Yin et al. [START_REF] Yin | Impact on flight trajectory characteristics when avoiding the formation of persistent contrails for transatlantic flights[END_REF] use genetic algorithms to compute transatlantic flight trajectory to mitigate the impact of contrails. Methods based on graphs are also used, like A* in [START_REF] Rosenow | Impact of multi-criteria optimized trajectories on European airline efficiency, safety and airspace demand[END_REF] or Dijkstra's algorithm [START_REF] Campbell | Multiscale path optimization for the reduced environmental impact of air transportation[END_REF].

Finally, other methodologies rely on MILP [START_REF] Campbell | An operational strategy for persistent contrail mitigation[END_REF] or on Mixed Integer Quadratic Programming (MIQP) [START_REF] Campbell | An optimal strategy for persistent contrail avoidance[END_REF] formulations.

Simorgh et al. remark in [START_REF] Simorgh | A comprehensive survey on climate optimal aircraft trajectory planning[END_REF] that Air Traffic Management (ATM) considerations are less often taken into account when multiple aircraft are considered, although the impact on airspace capacity and controller workload is certain. Indeed, a risk is to empty the spaces favorable to contrails by strongly congesting adjacent airspaces. Addressing problems that aim at avoiding such situations is one of the main contributions of this paper.

Assigning routes for several vehicles

Finding a shortest path is a common problem in operations research. The literature presents various way to solve such problems. The graph version of the problem can be addressed by integer linear programming [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF], or other efficient algorithms such as Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF], A* [START_REF] Hart | A formal basis for the heuristic determination of minimum cost paths[END_REF] or Bellman's algorithm [START_REF] Bellman | On a routing problem[END_REF] (dynamic programming). In some cases, optimal control techniques [START_REF] Luenberger | Introduction to Dynamic Systems: Theory, Models, and Applications[END_REF] can be used, especially if the path is to be computed in a continuous space and not on a graph. Variants of the shortest path problem are subject to constraints that typically involves an upper limit on a function of the arcs. For instance, the goal can be to minimize the distance traveled by a vehicle with an upper bound on the travel time. This type of problem is usually expressed for one vehicle, for one path. It is an NP-hard problem for which some efficient methods have been developed [START_REF] Lozano | On an exact method for the constrained shortest path problem[END_REF][START_REF] Pugliese | A survey of resource constrained shortest path problems: Exact solution approaches[END_REF][START_REF] Santos | An improved solution algorithm for the constrained shortest path problem[END_REF][START_REF] Huang | The shortest path problem on a fuzzy time-dependent network[END_REF].

Some other problems compute several shortest paths, i.e., several vehicles are considered via a global criterion to be minimized. This is the case for the traffic assignment problem (TAP) which aims at reaching an equilibrium for the vehicles or for the whole system. The type of chosen equilibrium determines the objective function to be minimized. An example is the Wardrop user and system equilibrium [START_REF] Wardrop | Correspondance. Some theorical aspects of road traffic research[END_REF][START_REF] Wardrop | Road paper. Some theorical aspects of road traffic research[END_REF]. The problem is subject to flow conservation constraints and capacity constraints on arcs and, in the case of system equilibrium, it minimizes the average journey cost. This problem considers cooperation and possibly a decentralized management of the traffic. On the other hand, the user equilibrium is reached when no vehicle can lower its transportation cost through unilateral action. It is for instance used in road network applications, where roads do not have infinite capacity but the number of vehicles that can use each road in a given amount of time is limited.

The previously mentioned problems take into account the arcs of the graph to define routes, and possibly capacities on these routes. However, there are also other types of problems in which the grain is coarser: nodes are defined by geographical sectors, and a total capacity on each of these sectors is imposed.

This special case therefore involves capacities on the vertices of the graph. This type of problem appears in air transportation: it is then named air traffic flow management problem (ATFMP), originally defined in [START_REF] Bertsimas | The air traffic flow management problem with enroute capacities[END_REF]. The objective function of the ATFMP is the total cost of aircraft delays, but variants can be derived by considering other objective functions. For example, the total cost of trajectories in terms of flight time, distance flown or CO 2 emitted can be taken into account. Optimization models addressing this problem typically involves the following decisions to be made for each flight:

• which sector to fly and when (which may induce speed modulations)? • when to take off (by imposing delays with respect to the scheduled departure time)?

Various capacity upper bounds are imposed on:

• the total number of aircraft in each sector at any given time (sector capacity constraints), • the total number of aircraft in each airport at any given time (airport capacity constraints). This is a large-grained problem, but it is also necessary to take into account the more precise spatial scale of the arcs to know where to fly through a sector at a given time. This is done in the variant of the ATFM problem that involves rerouting (ATFMRP), and is solved in [START_REF] Agustin | On air traffic flow management with rerouting. Part I: Deterministic case[END_REF] (in its deterministic version). This problem is very detailed since it completely defines the trajectory followed, and it decides the speed of each aircraft on flown arcs. It takes into account a cost on each arc. This cost can be the flight time, the distance flown, or estimated CO 2 emissions on this arc, for each aircraft. This paper proposes a method to organize traffic, while ensuring safety by satisfying capacity constraints on sectors. Moreover, it is also able to take into account contrails for minimizing the traffic environmental impact.

Mathematical optimization model

This section presents the mathematical optimization models for the subgraph-capacity multiple shortest path problem: Subsection 3.1 presents the static case while Subsection 3.2 addresses the time-dependent case.

Static case

This subsection focuses on the subgraph-capacity multiple shortest path problem in the static case (input data do not evolve with time).

The classical shortest-path problem on a graph involves only one vehicle and is modelled as follows. Let G = (V, A) be a directed weighted graph, where V is the set of vertices, A is the set of arcs, and w : A → R is the weight function. Let s ∈ V be the start of the path, and e ∈ V be the end of the path. Considering that the decisionvariable vector X has a component x u,v for each arc (u, v) ∈ A, where x u,v indicates whether the arc (u, v) is part of the solution path or not, the optimization formulation of the shortest-path problem is:

min X (u,v)∈A w u,v x u,v (1a) s.t.
(u,v)∈A

x u,v - (v,u)∈A x v,u = 0, u ∈ V \{s, e} (1b)
(s,v)∈A x s,v - (v,s)∈A x v,s = 1 (1c) (e,v)∈A x e,v - (v,e)∈A x v,e = -1 (1d) X ∈ {0, 1} |A| , (1e)
where (1b), (1c) and (1c) are the classical network flow conservation constraints (see for instance [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF]) that ensure that the solution is a path from s to e (the incoming flow is equal to the outcoming flow on each vertex other than the source and the end). This model can be adapted in the case of several vehicles with vehicle-specific weight functions, whose (u, v)-components are noted w u,v,i , (u, v) ∈ A, i = 1, 2, . . . , M , where M is the number of vehicles, and for each vehicle i = 1, 2, . . . , M , a start vertex s i ∈ V , and an end vertex e i ∈ V are given. Defining a decision-variable vector X i for each vehicle, i = 1, 2, . . . , M , the optimization model for the multiple shortest-path problem then reads:

min X M i=1 (u,v)∈A w u,v,i x u,v,i (2a) s.t.
(u,v)∈A

x u,v,i - (v,u)∈A x v,u,i = 0, u ∈ V \{s i , e i }, i = 1, 2, . . . , M (2b)
(si,v)∈A x si,v,i - (v,si)∈A x v,si,i = 1, i = 1, 2, . . . , M (2c)
(ei,v)∈A x ei,v,i - (v,ei)∈A x v,ei,i = -1, i = 1, 2, . . . , M (2d)
X i ∈ {0, 1} |A| , i = 1, 2, . . . , M. (2e)
Finally, a subgraph-capacity extension can be defined, provided a partition N k=0 A k = A of the set of arcs is given with corresponding capacities C k , k = 1, 2, . . . , N , where N is the number of arc subsets considered. In the sequel, we shall call sector each of these arc subsets. We define an auxiliary decision-variable vector, Y i , for each vehicle i, i = 1, 2, . . . , M , whose k th component, y k,i , indicates whether vehicle i uses arcs of sector A k , k = 1, 2, . . . , N . The optimization model for the subgraph-capacity multiple shortest-path problem is then:

min X,Y M i=1 (u,v)∈A w u,v,i x u,v,i (3a)
s.t. (u,v)∈A x u,v,i - (v,u)∈A x v,u,i = 0, u ∈ V \{s i , e i }, i = 1, 2, . . . , M (3b)
(si,v)∈A x si,v,i - (v,si)∈A x v,si,i = 1, i = 1, 2, . . . , M (3c)
(ei,v)∈A x ei,v,i - (v,ei)∈A x v,ei,i = -1, i = 1, 2, . . . , M (3d)
M i=1 y k,i ≤ C k , k = 1, 2, . . . , N (3e)
y k,i = 1 if and only if (u,v)∈A k x u,v,i ≥ 1, i = 1, 2, . . . , M, k = 1, 2, . . . , N (3f)
X i ∈ {0, 1} |A| , i = 1, 2, . . . , M (3g)
Y i ∈ {0, 1} N , i = 1, 2, . . . , M. (3h)
Constraints (3b), (3c) and (3d) are the usual path flow conservation constraints for each vehicle, and constraints (3e) and (3f) are the new subgraph capacity constraints. Constraints (3f) enforce the definition of the auxiliary binary variables Y 1 , Y 2 , . . . , Y N . These constraints can be linearized and the y i 's can be relaxed into continuous variables following the proposition 1: Proposition 1. Each of the constraints (3f), i = 1, 2, . . . , M , can be replaced and linearized by:

y k,i ≥ x u,v,i , (u, v) ∈ A k , (4)
y k,i ≤ (u,v)∈A k x u,v,i , (5)
y k,i ∈ [0, 1]. (6)
Proof. Consider a sector k and a vehicle i. One can easily show that 1 -

y k,i = (u,v)∈A k (1 -x u,v,i ).
Applying then the classical Fortet linearization extended to a product of several binary variables (see [START_REF] Liberti | Reformulation techniques in mathematical programming[END_REF]Subsection 3.3.2]), one obtains:

1 -y k,i ≤ 1 -x u,v,i , (u, v) ∈ A k , 1 -y k,i ≥ 1 - (u,v)∈A k x u,v,i , 1 -y k,i ∈ [0, 1].
This yields straightforwardly the desired result.

Time-dependent case

This subsection defines the mathematical formulation of the problem in the case where input data evolve with time. In the general case, we are considering time-dependent costs. Moreover, constraints (3e) and (3f) can also be time-dependent in case where a sector is only occupied by the vehicle during a certain amount of time and it is not occupied when the vehicle is not in the sector yet/anymore. The associated timediscretized model is detailed in Appendix A.

To take into account the release of the capacity of sectors by vehicles as they move, the paths are optimized for a succession of (sliding) time windows. More precisely, the paths on the graph are computed based on sector occupancy during the time interval under consideration. Then, the time is incremented by the sliding-window length ∆t, and the start vertex s i of each vehicle i is updated: it is replaced by the vertex reached by vehicle i in the previous sliding-window optimization. This process is illustrated by Figure 1. When, at the end of a time window, a vehicle is on an arc but not at a node, an artificial node is created. The process is stopped when each vehicle i has reached its final destination vertex e i , i = 1, 2, . . . , M . The entry time in the simulation of each vehicle is not artificially changed to coincide with the beginning of a time window. This time-window approach may lead to suboptimal solutions. It is used as a resolution heuristic for the following reasons. First, the time-dependent shortest path for a single vehicle is already a difficult problem [START_REF] Foschini | On the complexity of time-dependent shortest paths[END_REF], not to mention the case involving several vehicles plus capacity constraints. Moreover, in the context of our air transport application, the time-window methodology is tailored to operational concerns. Indeed, all relevant information is not always known for the subsequent time slots, and this time-sequential approach reduces drastically the uncertainties at each time window, by taking into account updated information. The paths on the graph are computed taking into account the capacity constraints only for a duration corresponding to K times the length of the sliding window, where K is a user-defined input parameter. This is illustrated by Figure 2. More precisely, constraints (3f) are modified as follows:

Initial start

y k,i = 1 if and only if (u,v)∈A k,K ∆t x u,v,i ≥ 1, i = 1, 2, . . . , M, k = 1, 2, . . . , N, (7)
where A k,K ∆t is the subset of arcs from A k that have at least one end reachable by vehicle i within a time less than K ∆t from s i . Fig. 2: Capacity computation along time windows. We consider that a vehicle i consumes capacity of the sectors that are reachable from the start point s i within a time less than K ∆t (here: capacity is consumed only on red arcs). gives more explanations about the required input data.

General description of the application

In the sequel, the French upper airspace will be considered. Aircraft fly above France following a sequence of 3D or 2D points linked by straight-line segment routes. These points are named waypoints, and all the waypoints defined above France are the vertices of the graph. The segment routes defined between the waypoints define the set of arcs of the graph. The graph is thereby a directed graph, since aircraft routes are typically constrained to one direction. The problem of finding a trajectory for a flight therefore boils down to finding a path in a graph. The aim here is to minimize the global impact of all flights on the environment. More precisely, the goal is to minimize the total environmental cost while ensuring that the airspace capacities are not exceeded. The controller point of view is taken into account, and fairness among aircraft and airlines must be kept in mind. For this reason, and in order to avoid suboptimal solutions, a sequential (one aircraft at a time), greedy-like, computation of trajectories cannot be used, even though there are very efficient algorithms to solve single-vehicle shortest-path problems.

The problem is designed to compute only the cruise part of the trajectories. Indeed, departing and arriving an airport is subject to numerous extra operational constraints that leave almost no degree of freedom. Moreover, the altitudes concerned by the contrails are cruising altitudes which are quite high. It is therefore generally not relevant to address low-altitude parts of the trajectories. Then, the application we consider in this study is solved only in the 2D plane, i.e., the altitude is not to be decided by the optimizer. The altitude is generally little modified during the cruise phase for reasons of passenger comfort. It is moreover already optimized by airlines to minimize fuel consumption and reduce engine worn out. Finally, the point of view chosen in this study is that of air traffic control. However, the weight of the aircraft being part of commercial data, it is not known by air traffic control. Without this data, it is not possible to know the ability of the aircraft to climb and to quantify precisely the environmental impact of a climb or a descent.

In the case of ground transportation, applications involving for instance cars on a road network, one can set a limited capacity per arc. In our air transportation application, the capacity limit is on areas (subsets of arcs), called sectors. These sectors represent subdivisions of the upper airspace, and more precisely a partition of the set of arcs. If the trajectories were to be computed sequentially, then some vehicles, the first ones, would be favored over the others. However, in the case of air transport, from the point of view of air traffic control, no airline should be favored over another, so such a sequential approach is not satisfactory. Above all, on the simple point of view of optimization, the greedy-like sequential approach is likely to lead to undesirable, suboptimal solutions.

Remark that satisfying capacity constraints allows one to ensure safety on the strategic time scale (several hours before the flights) that is considered here. It guarantees that the controllers will be able to manage the number of aircraft in their sector, while avoiding possible losses of separation in the shorter term.

Airspace capacities evolve with time and their occupancy evolve as aircraft enter or leave airspaces. Moreover, weather data, and also wind and contrail area data, evolve with time. The sliding-window approach is particularly adapted in this case for several reasons. Airspace capacities are defined on given time slots (in number of aircraft per time slot), and can be adapted from a time slot to another. It is also adapted for contrail area consideration since weather prediction are given for a period of time and are refined as time runs.

It should be noted that sectors can, in some cases, be reconfigured (sector grouping or splitting) to accommodate demand. In this study, the areas are considered static (the airspace configuration does not evolve with time).

Cost computation

As explained before, the goal here is to take into account CO 2 and non-CO 2 effects, and our application focuses on contrails. Bi-objective optimization is therefore a natural point of view to balance the two criteria that may be contradictory (burn more fuel -CO 2 -to avoid contrails -non-CO 2 ). However, the aim is to minimize the global environmental impact, and for this we use a metric that is common in the climatechange literature [START_REF] Simorgh | A comprehensive survey on climate optimal aircraft trajectory planning[END_REF][START_REF] Sridhar | Aircraft trajectory optimization and contrails avoidance in the presence of winds[END_REF][START_REF] Soler | Flight trajectory design in the presence of contrails: Application of a multiphase mixed-integer optimal control approach[END_REF] to balance the two effects.

The cost, w u,v,i of an arc (u, v) for aircraft i is therefore the result of the (weighted) sum of the CO 2 and non-CO 2 costs:

w u,v,i = w CO2 u,v,i + w non-CO2 u,v,i
. As mentioned before, the only non-CO 2 impact considered in this study is the contrails phenomenon:

w non-CO2 u,v,i = w contrails u,v,i
. Since the goal is to minimize the overall impact, it is necessary to quantify the relative impact of contrails versus that of CO 2 . For this, a metric known as Global Warming Potential (GWP) is used. This metric relates the impact of most greenhouse gases to the impact of CO 2 , under the form of a multiplicative factor, considering that the impact of CO 2 corresponds to GWP = 1. This metric depends on a time horizon over which the impact is computed. 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 according to the different considered time horizons, H. In the sequel, the contrail-induced cirrus (CIC) GWP will be noted g H . Depending on the time horizon, H, considered, the contrails have a more or less weight. Notably, they have less impact compared to CO 2 in the long term. In the following, only cirrus clouds induced by contrails will be taken into account, since they are the most impacting effect on the climate. However, other contrails can also easily be considered by our model, by a simple change in the cost function. Then, the cost for an arc (u, v) flown by aircraft i is:

w u,v,i = (1 + λ u,v g H ) w CO2 u,v,i
, where λ u,v ∈ [0, 1] is the proportion of the arc (u, v) that lies is in a persistent contrail area.

Concerning the cost of CO 2 , it is the total quantity of CO 2 emitted by the aircraft as it flies over the arc. The amount of CO 2 emitted per liter of standard jet fuel is constant. This cost is then directly proportional to the fuel consumption. Then,

w CO2 u,v,i = C CO2 f u,v,i t u,v,i
, where C CO2 is the constant quantity of CO 2 emitted by 1 kg of standard jet fuel, f u,v,i is the fuel flow of aircraft i above the arc (u, v), and t u,v,i its flight time. Then, the cost function reads:

w u,v,i = (1 + λ u,v g H ) C CO2 f u,v,i t u,v,i . (8)
In the case studied, the aircraft evolve in the 2D plane and are assumed to not change speed. The fuel flow is then more or less constant and can be approximated by a representative fuel flow (given altitude and speed, representative mass of the cruise phase), f r i . Then, the cost function is rather wirtten as:

w u,v,i = (1 + λ u,v g H ) f r i t u,v,i , (9)
where multiplicative constants have been removed.

As mentioned earlier a sliding-window approach is adapted to the case of such weather-dependent cost functions, since there are uncertainties on weather forecast. In particular, contrails are difficult to predict, and their impact is even more difficult to predict [START_REF] Gierens | How well can persistent contrails be predicted?[END_REF]. Considering short-term path computation mitigates the uncertainties.

Details about how weather data are processed for the cost-function computation can be found in Appendix B.

Fairness consideration

In our application context, because of the air traffic control point of view, ensuring a certain level of fairness between aircraft and therefore between airlines is an issue that may emerge. The lack of fairness can occur on two points. First, satisfying airspace capacity constraints may result in detours for some aircraft and direct routes for others, and therefore some airlines are likely to be disadvantaged. Secondly, the avoidance of contrails has an economic cost for the airlines, and it is necessary to verify that an airline is not systematically penalized.

To address this issue, a penalty for flight time deviations from a reference flight time, noted t ref i for aircraft i, can be added to the objective function. Thus, the objective function becomes:

1 N N i=1 (u,v)∈A w u,v,i x u,v,i + µ max i=1,...,N (u,v)∈A t u,v,i x u,v,i -t ref i , (10)
where µ is a positive user-defined weighting parameter. The objective function then becomes non-linear, but can be linearized with classical reformulation techniques [START_REF] Liberti | Reformulation techniques in mathematical programming[END_REF].

In this new objective function, the average cost is minimized as well as the penalty for the largest detour. The constant µ is used to give the desired relative importance to the worst-case penalty cost. The parameter t ref i is a reference time of flight for flight i. Its value can be, for example, the flight time when no capacity constraint is imputed and no contrail is taken into account if the user wants to consider the two criteria causing loss of fairness. On the one hand, if the user wants to avoid that the capacity constraints lead to loss of fairness, then t ref i can be set to the flight time without capacity constraint but with contrails taken into account. On the other hand, if the user wants to limit the impact of contrails on the loss of fairness, this reference time can be set to the flight time without consideration of the contrails but with consideration of the capacity constraints.

In this paper, it is chosen not to present results with this new objective function for two reasons. First, fairness is mainly taken into account by not computing the trajectories sequentially, in a greedy-like manner, but "all at once" to avoid systematically favoring some flights. Second, the goal is to show the possibility of reducing the impact of contrails at the network level, and to quantify these impacts. The results could therefore be distorted by taking fairness into account. In addition, other techniques can be used to ensure fairness rather than minimizing the impact of the largest detour (see for example [START_REF] Ogryczak | Fair optimization and networks: A survey[END_REF]).

Data

This subsection details the input data of the problem. It explains in particular how the graph is built, and how the sectors are defined. The wind encountered and the areas favorable to contrails are also known data but the process to obtain the related information is detailed in Appendix B.

The present study follows the new principle of Free Route Airspace (FRA) which is applied nowadays to the European upper airspace. The aim of FRA is to remove the previously-established principle of air routes, and to replace it by navigation points, called waypoints, through which aircraft pass freely. A flight plan is therefore a simple sequence of waypoints through which the aircraft flies. This new paradigm allows one to consider an increased number of possible direct routes. As a consequence, the distance flown and thereby the CO 2 emissions, can be decreased.

Rules are still established to fly from a point to another even if there number aims to be decreased. Here, these rules are approximated by the rules established to build the graph G = (V, A) required for the optimization model. The vertex set, V , is the set of considered waypoints. The arcs (set A) connect two waypoints when their inter distance is less than some user-defined threshold distance, D. One could alternatively consider the complete graph, but this would not be coherent with operational practice, not to mention the increase of complexity in regards with the preliminary nature of the present study.

We construct instances based on the French waypoints that are located in the western and southern parts for now, see Figure 3 (the arcs are not displayed as they depend upon the maximum threshold distance, D, chosen).

Our set of arcs is therefore:

A = {(u, v) | u ∈ V, v ∈ V, d u,v ≤ D}, ( 11
)
where V is the set of waypoints, and d u,v is the distance between waypoints u and v.

For the subgraph capacities, the partition of the set A of arcs into N subsets, called sectors, {A k } k=1,2,...,N , is initially defined with respect to the waypoint set, as shown in different colors in Figure 4 for our instances. Then, all arcs with one of its ends in one such waypoint set, noted V k for same k = 1, 2, . . . , N , is considered to be in the (arc-set) sector:

A k = {(u, v), u ∈ V k or v ∈ V k }.
The next section reports computational results.

Fig. 3: The FRA waypoints above France constituting our instance set.

Fig. 4: The vertex subsets V k 's considered above France for building the arc set sectors A k 's for our instances.

Results and sensitivity analysis

This section presents an illustrative instance of the subgraph-capacity multiple shortest path problem together with various results obtained from numerical experiments resulting from a sensitivity analysis of the different parameters involved. An example of results is given in Subsection 5.1. Then, Subsection 5.2 focuses on the impact of wind on the results, while Subsection 5.3 addresses the impact of the time horizon chosen for the GWP computation on the results. Finally, the impact of the imposed airspace capacity, C k , for sector A k , k = 1, 2, . . . , N , is discussed in Subsection 5.4.

Description of the illustrative instance

In the sequel, one instance of the problem is addressed but with various cost functions and various levels (right-hand side) for the capacity constraints. This subsection details the definition of this instance.

The graph is computed thanks to rules detailed in Section 4 (see Figures 3 and4). Twenty aircraft are entering (taking off or entering the French upper airspace) per 30 minutes simultaneously, for three hours. The source-end vertex pair (s i , e i ) of each aircraft i, i = 1, 2, . . . , M , are chosen randomly in the vertex set so that the minimum distance (as the crow flies) is 200 Nautical Miles and the resulting instance is named FRA-200. A Nautical Mile (NM) is the distance unit used in aeronautics, and corresponds to 1.852 km. In order to define the arc set, A, we set the maximum distance between two linked points at D = 75 NM. The airspeed of all aircraft is set to 400 knots (a knot (kt) is the speed unit used in aeronautics and corresponds to 1 NM/h or 1.852 km/h). This speed is chosen in adequacy with the typical airspeed of a standard commercial aircraft, namely the Airbus A320 [START_REF] Airbus | SKYbrary Aviation Safety[END_REF]. Then, to simplify the presentation, the fuel flow is considered the same for all the flights (identical aircraft types, engines and weights). Figure 5 displays the wind encountered, and Figure 6 shows the persistent contrails areas. The length of the sliding window for time-dependence consideration is set to ∆t = 15 minutes. The parameter K for capacity consumption (as defined in Equation ( 7)) is set to K = 1.25. If nothing else is explicitly mentioned, it is considered that:

• the capacity of each sector k is set to C k = 20, k = 1, 2, . . . , N ; • the time horizon chosen for the GWP computation is set to H = 100 years. Table 2 summarizes the features of instance FRA-200. The details about the instance can be found in [START_REF] Demouge | FRA-200 Instance for subgraph-capacity multiple shortest path application to contrail avoidance[END_REF]. Instance FRA-200 is solved using the sliding window process explained in Subsection 3.2. The following results are obtained with the Java API of CPLEX on a computer with an Intel Core i5-10210U, 1.60 Hz, with 8 Go RAM and a Debian Linux OS. The code is available at [START_REF] Demouge | Time dependent contrails avoidance at network scale[END_REF]. Figure 7 shows the solution obtained on instance FRA-200 on maps. The trajectories are represented by groups of 20 aircraft on different maps for reading purposes. The computation time is 137.2 seconds.

Impact of wind

The wind has a certain impact on the results. To quantify this impact, two variants of instance FRA-200 have been solved:

1. without wind and without contrails:

w u,v,i = d u,v , (u, v) ∈ A, i = 1, 2, . . . , M ; 2.
with wind and without contrails:

w u,v,i = t u,v , (u, v) ∈ A, i = 1, 2, . . . , M .
Figure 8 displays the additional flight distance and flight time in the second case (with wind) in comparison to the first case (without wind).

The results show that in general, the flight distance is increased to the benefit of a reduced flight time, which was expected.

Impact of the time horizon used for GWP computation

As explained in Subsection 4.2, the cost factor associated with contrails depends on the chosen time horizon, H. Intuitively, the shorter the time horizon, the more impact the contrails have, and therefore the more beneficial it is to lengthen the trajectories to avoid contrails. To confirm this thought, three variants of instance FRA-200 are solved:

1. no contrail consideration: g H = 0; 2. contrails consideration with H = 100 years: g H = g 100 = 0.63; 3. contrails consideration with H = 20 years: g H = g 20 = 2.2.

Figures 9 and 10 display comparative results for variants 1 and 2, and for variants 1 and 3, respectively.

The flight time is higher because of the avoidance of contrail area. Figures 9 and10 show that the dependence on the time horizon is important since in the case of a short horizon, the optimal trajectories are much longer than in the case of a longer time horizon. If the contrails impact is considered high, the flight time has a lower

w u,v,i = d u,v ).
impact on the objective function, and then is highly increased to avoid as much as possible the areas favorable to contrails. Indeed, Figure 10 shows optimal trajectories with an increase of more than 30% in flight time. Fuel consumption therefore explodes. From an environmental point of view, it is not necessarily relevant to choose such a trajectory. It is also economically costly (increase in fuel and extra cost associated to flight delays). Finally, such an extra consumption has also an impact on safety since a significant part of the extra fuel that every aircraft must carry in order to anticipate a possible diversion to another airport will be used for avoiding contrails. An extra constraint should therefore be added to our optimization model to avoid such undesirable solutions. It is interesting to consider adding other extra constraint so as to satisfy other operational needs. For example, if an airline is planning connections between its flights, it is desirable that the flights involved are not delayed too much so that the connection is guaranteed.

Let t 0,i be the flight time without any contrail consideration. It can be computed thanks to the optimization process with g H = 0. We propose to add to the model

constraints (u,v)∈A t u,v,i x u,v,i ≤ C t 0,i , i = 1, 2, . . . , M, (12)
where t u,v,i is the flight time over the arc (u, v) for the flight i, and where C is a user-defined constant to control the additional flight time allowed. Figure 11 shows an example of results obtained when such constraints are added with C = 1.1, (10% of additional flight time is allowed). Because of the sliding-window process, the extra constraints are not globally satisfied, but they are satisfied locally at each step, and are nearly satisfied at the global scale, as can be observed in Figure 11.

Impact of airspace capacity

A last parameter to be studied is the airspace capacity values C k 's. The same setup concerning aircraft is taken with a GWP computed with, this time, H = 100 years. Several results are impacted by changing the airspace capacity C k , k = 1, 2, . . . , N , especially the obtained optimal objective-function value, and the computation time.

Indeed, in order to satisfy capacity constraints, some aircraft may be forced to fly longer or through contrail zones. The problem can be harder to solve if the capacity constraints are restrictive, and so the computation time increases. Figure 12 shows the evolution of the optimal objective-function value and the computation time obtained for various capacity levels, C k , assumed constant for every sector k, k = 1, 2, . . . , N , and at each time window. The above results show how the right-hand side of the airspace capacity constraints affect the feasibility of the problem and the environmental impact. The higher the capacity of the airspace is, the better the results are. However, beyond a certain value of C k , no change is observed in computation time, nor in objective-function value since the capacity constraints are not saturated.

Conclusion

Motivated by the issue of contrail-avoidance at the network scale from the air traffic control point of view, we introduced a new model for the multiple shortest-path problem that takes into account capacity constraints on subgraphs. In the static case, the number of vehicles going through each considered subgraph is counted and bounded. Dynamic aspects have also been taken into account by establishing an adapted model with time windows over each of which the static model is solved. Taking this timedependent aspect into account is essential because as vehicles move, they free up capacity, and costs can also change over time. The model used does not provide the same quality of solution as if the static model had been discretized in time, but it fits better the operational reality as illustrated by the application addressed in our study. Moreover, the proposal methodology allows one to use a standard solver to solve realistic air traffic instances.

The application addressed in this paper motivated the introduction of a new model for the contrail-avoidance problem when considered at the network scale and from the air traffic control point of view. The time-window strategy is particularly adapted to model time-dependence since it corresponds to operational realities (time between contrail forecasts and uncertainty on their prediction, airspace capacities defined per time period, etc.). The trajectories are computed for all aircraft and not for each individual aircraft to avoid favoring any particular flight. Work has also been done for computing cost of arcs, so that they can be adapted to the Air Traffic Management point of view. This allows one to consider and to weight CO 2 and contrails impact using Global Warming Potential with different time horizons, giving more or less importance to contrail avoidance. Moreover, it allows one to take into account other non-CO 2 effects due to air transport, or to use more elaborated metrics to determine the impact of contrails.

Future tracks of research should focus on the numerous sources of uncertainty to be taken into account when addressing the contrail-avoidance problem. There are several uncertainties to consider when studying contrail-favorable (or persistentcontrail-favorable) areas, as shown by Gierens et al. [START_REF] Gierens | How well can persistent contrails be predicted?[END_REF]. Other source of uncertainty are the wind estimation and the presence of traffic.

Moreover, when considering the air traffic control point of view, one faces another crucial criterion: fairness for which ideas have been evoked in this paper. Future work should adapt the proposed model so as to quantify and enforce fairness between airlines.

Finally, other decision variables could be considered. First, in this study, the airspace configuration is considered static, the sectors are not modified. However, reconfiguration may also be involved in decisions to accommodate the location of contrail areas. Also, the cruise altitude should be envisaged as decision variables since it is an efficient mean of mitigating contrail impact, as shown by Fichter et al. [START_REF] Fichter | The impact of cruise altitude on contrails and related radiative forcing[END_REF].

This leads to several issues including the knowledge of some critical aircraft parameters (such as its weight at any moment) which are not always known from the air traffic control point of view. Moreover, fuel consumption depends on the altitude as well as on the type of aircraft. Models to estimate fuel flows from an air traffic control perspective should be considered. For that, the OpenAP [START_REF] Sun | OpenAP: An open-source aircraft performance model for air transportation studies and simulations[END_REF] database could be used, or other models based on machine learning as it has been done in [START_REF] Jarry | Approach and landing aircraft on-board parameters estimation with LSTM networks[END_REF] for the approach and landing phases. Supplementary information. Not applicable.

B.1 Wind data

The costs defined by ( 9) involves the computation of the flight time over each arc (u, v) ∈ A. To compute these costs, the wind on the arcs, and the distance between u and v, the two ends of the arcs, have to be known.

Let λ u and λ v be the latitude of vertices u and v respectively, and let ϕ u and ϕ v be their longitude. The distance between u and v is given by: Then, the flight time for aircraft i, noted t u,v,i , between points u and v is given by:

d u,v = R
t u,v,i = d u,v GS u,v,i , (B5)
where GS u,v,i is the ground speed of aircraft i on arc (u, v). It can be computed via:

GS u,v,i = V ai + W u,v , (B6)
where V ai is the airspeed of aircraft i (considered constant), and W u,v is the wind encountered on the arc (u, v).

Wind data have been extracted from the website Windy [START_REF] Windyty | Windy API[END_REF] on a square grid of size 0.2 degree above France, as shown on Figure B1. To compute the wind on each node of the graph, a so-called Shepard interpolation [START_REF] Shepard | A two-dimensional interpolation function for irregularly-spaced data[END_REF] was used. More precisely, for each node P located in a 2D-square P 1 P 2 P 3 P 4 of the data grid (see Figure B2), the wind W (P ) at P is calculated from the wind at P k , k = 1, 2, 3, 4, using the distance from P to each of these points (noted respectively d 1 , d 2 , d 3 and d 4 ) as follows:

i=1 W (P i ) * d -p i 4 i=1 d -p i , (B7)
where p > 1 is a user-defined parameter (set to p = 2 in this study). Finally, the wind along an arc (u, v) is simply defined as the average of that at u and at v:

W (u,v) = W (u) + W (v) 2 . ( B8
)
Figure B3 shows data used for the examples presented in the result section (Section 5).

B.2 Contrail data

Contrails are formed in cold and humid areas. They persist and induce cirrus if the air is supersaturated in ice. The computation of persistent contrail areas is performed in two phases:

1. Areas favorable to contrail formation 2. Areas in which contrails will persist (ice supersaturated areas).

In [START_REF] Soler | Flight trajectory design in the presence of contrails: Application of a multiphase mixed-integer optimal control approach[END_REF], contrails areas are computed thanks to the Schmidt-Appleman criterion. This criterion gives a minimum threshold, r min , of relative humidity of the air in liquid water, noted RH w , above which contrails are formed: contrails are assumed to form when RH w ≥ r min , where , EI H2O = 1.25 is the water vapor emission index, C p = 1004 J.kg -1 .K -1 is the heat capacity of the air, P is the ambient pressure (in Pascals), ϵ = 0.6222 is the ratio of the molecular masses of water and dry air, Q = 43 * 10 6 J.kg -1 is the specific heat of combustion, and η = 0.3 is the average propulsion efficiency of a commercial aircraft.

In [START_REF] Soler | Flight trajectory design in the presence of contrails: Application of a multiphase mixed-integer optimal control approach[END_REF], the ice super saturated areas are determined thanks to the following criterion: RH i > 1, where the relative humidity over the ice, noted RH i is computed as follows: and where T is the ambient temperature in Celsius degrees.

RH i =
The relative humidity and temperature are also computed from data extracted from Windy [START_REF] Windyty | Windy API[END_REF] on a 2D-grid, and interpolated via quadratic interpolation. Figure B4 shows the data used for the examples presented in the result section (Section 5), where red areas are persistent-contrail-favorable areas, to be avoided.

  Compute the path from the starting point to the final destination point.

  Define the new starting point.

  Compute the path from the new starting point to the final destination point.

Fig. 1 :

 1 Fig. 1: Sliding window computation of a shortest path. The final path is computed sequentially, in pieces, for each of the time windows.

4

  Application to contrail avoidance and CO 2 minimization Subsection 4.1 shows that contrail mitigation for several aircraft can be seen as an instance of the time-dependent subgraph-capacity multiple shortest path problem modeled in Section 3. Then, the cost function computation is detailed in Subsection 4.2, before focusing on the fairness issue in Subsection 4.3. Finally, Subsection 4.4

Fig. 5 :

 5 Fig. 5: Wind encountered in instance FRA-200.

Fig. 6 :

 6 Fig. 6: Persistent-contrail areas used for instance FRA-200.

Fig. 7 :

 7 Fig. 7: Results obtained on instance FRA-200, grouped by 20 aircraft, sorted by increasing entry time.

Fig. 8 :

 8 Fig. 8: Additional flight time flight distance when wind is considered (flight time minimization) in comparison with results obtained when wind is not taken into account (only distance is minimized: w u,v,i = d u,v ).

Fig. 9 :

 9 Fig. 9: Additional flight time and flight distance when H = 100 years (g H = g 100 = 0.63) when compared with the case without contrails (g H = 0).

Fig. 10 :

 10 Fig. 10: Additional flight time and flight distance when H = 20 years (g H = g 20 = 2.2) when compared with the case without contrails (g H = 0).

  time constraint Mean additional flight time : 2.97%

Fig. 11 :

 11 Fig. 11: Additional flight time and flight distance with extra flight time constraints added when H = 20 years, in comparison with the case without contrails.

Fig. 12 :

 12 Fig. 12: Result comparison in terms of objective-function value and computation time when the airspace capacity, C k , changes.

  c rad in km, (B2) = 60 c degrees in NM, (B3) where R = 6, 371 km is the Earth radius, c rad and c degrees represent the following c values, expressed in radians and in degrees respectively: c = arccos sin(λ u ) sin(λ v ) + cos(λ u ) cos(λ v ) cos(ϕ v -ϕ u ) . (B4)

0. 2 •0. 2 •

 22 Fig. B1: 2D grid used for wind data extraction from Windy [42].

3 Fig. B2 :

 3B2 Fig. B2: Notations for estimating the wind at P via Shepard interpolation.

Fig. B3 :

 B3 Fig. B3: Wind encountered in the example of the result.

  r min = G(T -T c ) + e liq sat (T c ) e liq sat (T ) , (B9) e liq sat (T ) is the saturation vapor pressure over water, and T c is the estimated threshold temperature (in Celsius degrees) for contrail formation at liquid saturation. The later is computed via: T c = -46.46 + 9.43 log(G -0.053) + 0.72 log 2 (G -0.053), (B10) where G = EI H 2 O CpP ϵQ(1-η)

Fig. B4 :

 B4 Fig. B4: Persistent-contrail areas in western France used for our instances.

Table 1 :

 1 Global warming potential for contrail for various time horizons, H.

		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 Features of the illustrative instance, FRA-200.

	Feature	Notation	FRA-200
	Aircraft and trajectory features
	Number of aircraft	M	20 per 30 minutes for 3 hours
	Airspeed of aircraft		400 kts
	Minimum distance between s i and e i		200 NM
	Airspace features
	Number of sectors	N	23
	Maximum distance for arc definition	D	75 NM
	Capacity of each sector	C k	20
	Optimization model features
	Time-window length	∆t	15 minutes
	Time horizon for the GWP computation	H	100 years

  RH w * 6.0612 * exp( 18.102 * T 249.52+T ) 6.1162 * exp( 22.577 * T 273.78+T )

	,	(B11)

National or international entity that aims to ensure the best possible use of the airspace according to its capacities by adapting the traffic. In Europe, the network manager is Eurocontrol. In the United States, the equivalent is the Air Traffic Control System Command Center (ATCSCC).

Acknowledgments. The authors thank DGAC (French civil aviation authority) for prompting and funding this work, and more specifically its DTA and DSNA services. Their inputs and expertise were also essential to the achievement of this work.

Availability of data and materials. Sample data are available at: https://cloud. recherche.enac.fr/index.php/s/i6jxDFM8GnSAgyF.

Code availability. Code is available in the following git repository: http://celine.demouge@git.recherche.enac.fr/time_dependent_contrails_avoidance_ at_network_scale.git.

Author Declarations

Funding. The PhD research work of Céline Demouge is funded by DGAC (the French civil aviation authority).

Conflict of interest/Competing interests. The authors have no conflicts of interest to declare that are relevant to the content of this article.

Ethics approval/declarations. Not applicable.

Appendix A Time-discretized optimization model

The mathematical optimization model obtained in the static case can be discretized to obtain the time-dependent optimization model. This transformation is based on the classical time-dependent shortest path problem [START_REF] Yang | Constraint reformulation and a Lagrangian relaxationbased solution algorithm for a least expected time path problem[END_REF].

Some notations should be defined first:

• T i = {t 0,i , . . . , t f,i }: the set of time slots for vehicle i;

u,v,i : the time necessary for vehicle i, i = 1, 2, . . . , M , to go through arc (u, v) ∈ A at time slot t, t ∈ T i .

The decision variables are for each vehicle i, i = 1, 2, . . . , M :

The model is then:

(u,v)∈A

Constraints (A1b),(A1c) and (A1d) are the flow conservation constraints. Constraints (A1e) and (A1f) enforce consistency of space-time flow conservation. Constraints (A1g) make the link between the time-dependent and the static decision variables. Constraints (A1h) are the time-discretized capacity constraints. Constraints (A1i) and (A1j) define the auxiliary variables y t k,i . The former can easily be linearized, as for the static model using Proposition 1.

Appendix B Data processing for numerical experiments

This section shows how data for computational experiments have been extracted and computed. Subsection B.1 deals with wind data, while Subsection B.2 focuses on contrail data.


Nassia Daouayry
email: nassia.daouayry@airbus.com

Ammar Mechouche
email: ammar.mechouche@airbus.com

Pierre-Loic Maisonneuve
email: pierre-loic.maisonneuve@airbus.com

Vasile-Marian Scuturici
email: marian.scuturici@insa-lyon.fr

Jean-Marc Petit
email: jean-marc.petit@insa-lyon.fr

Data-Centric Helicopter Failure Anticipation: The MGB Oil Pressure Virtual Sensor Case

Keywords: Big Data, Virtual Sensor, Helicopter, HUMS, Predictive maintenance, Machine Learning

This paper presents a virtual sensor for helicopter Main Gear Box (MGB) oil pressure. It is developed for anticipating failures of systems related to MGB lubrication.

The virtual sensor is built using Machine Learning and massive in-service flight data collected from Airbus helicopters flying world-wide. The correlation between oil pressure values and other flight parameters is learnt during stable phases of flights in which the system is in its nominal state.

At each flight, the values continuously estimated by the virtual sensor are compared to the measured ones, and an alert is raised when the difference becomes higher than a statistically predefined threshold.

The virtual sensor was tested using normal and abnormal flights, and the results obtained so far in terms of anomaly detection performance are promising.

I. INTRODUCTION

Failure anticipation, and more generally predictive maintenance, is an important and growing topic in the industry [START_REF] Vianna | Predictive Maintenance Optimization for Aircraft Redundant Systems Subjected to Multiple Wear Profiles[END_REF][2][3][4] [START_REF] Giobergia | Tuninetti: Mining Sensor Data for Predictive Maintenance in the Automotive Industry[END_REF], especially in the aeronautic domain, due to strong requirements in terms of safety and high costs of unscheduled events [START_REF] Alestra | Rare Event Anticipation and Degradation Trending for Aircraft Predictive Maintenance[END_REF][7] [START_REF] Wade | Validating Machine-learned Diagnostic Classifiers in Safety Critical Applications with Imbalanced Populations[END_REF][9] [START_REF] Sundareswara | Spectrum Reconstruction for Oil Cooler On Condition Monitoring[END_REF]. Failure anticipation is a condition-based maintenance concept and consists in anticipating failures -which are unscheduled events -before they happen. In addition to obvious safety enhancement benefits, anticipating failures is also motivated by economic aspects, which may represent a valuable service for customers. It allows them to: 1) prepare and plan for the required maintenance operation instead of suffering it; 2) order the required part in order to have it in stock on time, right before the maintenance operation; and 3) avoid operational interruptions that imply unexpected financial losses for operators, due to mission's cancellation. In summary, increasing helicopters availability at a lower cost is a key element for operators. Moreover, early failure detection of a component may prevent from causing damages to its neighbouring parts, and it can also make this component reparable before its complete failure. For example, failures of some components of the helicopter Main Gear Box (MGB) that are not detected at a very early stage may contaminate other non-replaceable components. This leads to the removal of that MGB which will be definitely lost, at the cost of around 250 K€.

During decades Airbus Helicopters (AH), like its competitors, focused its efforts on anticipating failures mainly for ensuring safety of helicopters [START_REF] Camerini | Fault detection in operating helicopter drivetrain components based on support vector data description[END_REF] [START_REF] Civil | CAA Paper: The application of advanced anomaly detection to tail rotor HUMS data: based on a report prepared for the CAA by GE Aviation Systems Limited[END_REF]. This is based on 'classical' Health and Usage Monitoring System (HUMS) [START_REF] Civil | Report of the CAA Working Group on Health Monitoring[END_REF]. HUMS mostly consists of counters monitoring quantities related to the usage of the machine, and high-frequency vibration measurements acquired from multiple sensors, typically located in the proximity of critical mechanical transmission components and dynamic systems [START_REF] Mechouche | Helicopter Big Data Processing and Predictive Analytics: Feedback & Perspectives[END_REF]. However, from one hand, embedded algorithms for on-board alarming and usage monitoring are rather simple since they are constrained by embedded hardware processing capabilities. They also only analyse "fresh" data generated during the ongoing flight, without any trend monitoring exploiting historical flight data. Consequently, alarms and alerts are raised when monitored sensors, quantities and systems are already faulty (I.e. no longer operational, exceeding a critical threshold limit or reaching an advanced level of wear). They either allow an emergency landing or cancel future flights until the required maintenance is completed. On the other hand, on-ground vibration monitoring is restricted to critical mechanical and dynamic components, whereas plenty of other systems, which are less critical from a safety point of view, are still not monitored. However, it is worth noting that whenever they are gathered to consolidate a global view of a helicopter in operation, they could reveal weak signals very useful to failure anticipation. Thus, there is a need for covering more helicopter systems and improving early detection of failures.

Considering harsh market competition and given the operational costs of helicopters (2/3 of the total cost), helicopters manufacturers are more and more investing beyond safety to generate new revenues. To do so, they develop new types of services to support their customers in reducing their costs. The massive continuous data generated by helicopters during every flight represent an important lever for developing such services. Analysis of these data might help for improving products, offering new services (e.g. to enlarge condition based maintenance to helicopter usage data) or gaining additional business insights.

In this setting, Airbus Helicopters gathers more than 400000 flight hours (FH) data (called in the sequel in-service flight data) for more than 10 years from hundreds of connected helicopters, operated by several customers worldwide and performing a large variety of missions (Figure 1). These time series data are recorded from the beginning until the end of each flight thanks to the Flight Data Continuous Recorder (FDCR) function in HUMS system. Specifically, hundreds of parameters are recorded from the FDCR, coming from each of the vehicle systems connected to the avionics [START_REF] Mechouche | Helicopter Big Data Processing and Predictive Analytics: Feedback & Perspectives[END_REF].

Figure 1: Data collection in the Big Data platform at AH

The overall objective of this paper is to take benefits from a broad variety of collected in-service flight data in order to anticipate and early detect MGB failures, using machine learning techniques. The idea is to build virtual sensors predicting different parameters from other measures as an alternative to physical sensors.

Problem statement: Given huge amount of in-service flight data, the objective is to build a virtual sensor for MGB oil pressure at each timestamp of the flight, in order to raise an alert as early as possible when a significant deviation between values measured by the physical sensor and those provided by the virtual sensor occurs.

The choice of MGB is motivated both for its criticality with respect to safety requirements and for economic reasons, i.e. especially in case of unscheduled and non-anticipated removals. Moreover, MGB turns out to be a highly monitored system with many physical sensors producing continuous data at high frequency, helpful for building a virtual sensor.

A supervised machine learning based anomaly detection approach is devised to build this virtual sensor. It is built thanks to learning on multi-varied time series data representing flight parameters. The parameter to be predicted is the MGB oil pressure, for which a physical sensor does exist in every helicopter MGB. To strengthen this data-centric approach, "domain knowledge" of experts is taken into account to understand the main parameters that should influence the oil pressure of an MGB. Starting from these parameters we have been able to come up with a new relationship, not known from a physical point of view, between the MGB oil pressure and a subset of other indicated flight parameters during a stabilized flight phase 1 . The idea of monitoring MGB oil pressure allows several systems to be covered that are not covered by vibration monitoring, including lubrication system, oil cooling system, oil pump, oil injectors and radiator filter. The point is that any anomaly in the operation of that MGB system should normally be detected by monitoring MGB oil pressure.

The proposed approach presented several challenges. First, correlation and analysis of big amounts of time series data require a very good knowledge of latest big data technologies and an adequate big data platform. At AIRBUS such platforms and tools are available, which facilitates helicopter big data processing and analytics [START_REF] Mechouche | Helicopter Big Data Processing and Predictive Analytics: Feedback & Perspectives[END_REF]. The second challenge consists in the constitution of a curated training dataset. In the aeronautic domain, failures remain very few compared to the number of flights. Hence, anomaly detection proposed here consists in modelling the normal state of the MGB. Fortunately, most data used for learning the model is composed of "healthy flights". Moreover maintenance actions exist since they are tagged by customers in their Maintenance Information Systems (MIS). The third challenge lies in the constitution of a complete training dataset covering all flight conditions. As an example, one of the impacting parameters of the MGB oil pressure is the Outside Air Temperature (OAT). Consequently, if a new flight is operated in a different region / continent with an OAT which has never been seen in the past, then the virtual sensor model is likely to diverge. The retained solution consists in learning on the whole fleet data already collected, and monitoring the drift of the model during the execution phase, in order to re-learn it when new and different data is collected.

To sum up, we have made the following contributions:

-A principled approach to address the problem of failure anticipation of an MGB in time series data using state of the art supervised machine learning techniques (Gradient Boosting).

-Integrating domain knowledge as early as possible into the virtual sensor's construction, to get domain expert as involved and confident as possible in its industrialization process.

-Ongoing deployment of the virtual sensor in the "Flyscan predictive maintenance program" of Airbus Helicopters 2 .

Regarding related work, virtual sensing has been developed in other domains but they do not rely on Machine Learning and big data [START_REF] Kadlec | Soft sensors: Where are we and what are the current and future challenges[END_REF][15] [START_REF] Kabadayi | Virtual sensors: abstracting data from physical sensors[END_REF]. However, similar anomaly detection methods exist, like in [START_REF] Alestra | Rare Event Anticipation and Degradation Trending for Aircraft Predictive Maintenance[END_REF][7] [START_REF] Benedetti | Anomaly detection and predictive maintenance for photovoltaic systems[END_REF], even if they are not presented in the form of virtual sensors. We exploit massive time series data collected from Airbus Helicopters customers, without any reduction as done for example in [START_REF] Camerini | Fault detection in operating helicopter drivetrain components based on support vector data description[END_REF]. 1 In this paper, a stabilized flight phase is one during which the helicopter altitude, the speed and the direction are maintained around specific value. The remainder of the paper is organized as follows: first an overview of helicopter data is provided as well as a general view of the problem and the proposed solution. Then, the proposed method is detailed, including validation results and tests on real 'known' failures. After that, the deployment in production within FlyScan service, which is the Airbus Helicopters service for predictive maintenance proposed to customers, is described. Finally, the obtained results are discussed and some perspectives are presented.

II. A PRINCIPLED APPROACH FOR MGB MAINTENANCE

A. Available data

The data used in this study comes from measures produced by the sensors installed on Airbus's helicopters. They are recorded continuously during each flight and then transferred, processed and analysed according to the different AH services. The FDCR data collected are recorded at the frequencies of 2 Hz. These time series are gathered on a Big Data platform (Hadoop / Spark).

We focused on a recent aircraft type, for which more than 60 000 FH are already collected. For this helicopter fleet, which encompasses more than 30 different aircraft serial numbers, a dataset has been created with flights data related to the MGB system. The size of the data extracted for the experimentation reported here is around ~180 gigabytes. Table 1 sketches how such data are structured and represented. Each parameter (f 0 , …f p ) is measured at the same frequency (2 Hz) and the elapsed time is represented as the number of 500 millisecond steps from the beginning of the flight ("time_step") identified by the "flight id". Additional information is available concerning for example the operating region, the specific helicopter associated to a flightID etc.

f0(1,1) f1(1,1) f2(1,1) f3(1,1) f4(1,1) fp(1,1) 1 2 f0(1,2) f1(1,2) f2(1,2) f3(1,2) f4(1,2) fp(1,2) .. … … … … … … … 2 1 f0(2,1) f1(2,1) f2(2,1) f3(2,1) f4(2,1) fp(2,1) 2 2 f0(2,2) f1(2,2) f2(2,2) f3(2,2) f4(2,2) fp(2,2) 2 3 f0(2,3) f1(2,3) f2(2,3) f3(2,3) f4(2,3) fp(2,3) … … … … … … … … 3 1 f0(3,1) f0(3,1) f2(3,1) f3(3,1) f4(3,1) fp(3,1) … … … … … … … …

B. Data pre-processing

Since the objective is to build a dataset corresponding to data issued from helicopters operating with a known "normal" state of the MGB.

Significant work was required to address the corresponding issues. Two main approaches are used to reduce the variance in the virtual sensor model. From one hand, abnormal flights are identified in the data cleaning phase and, from the other hand, variations of contextual parameters during a flight in the contextualization step.

Data cleaning

Identifying abnormal flight (regarding the MGB system) on the complete historical data is not an easy task. This has been done in collaboration with the MGB system designer. Several workshops have been organized to identify the types of flights that influence the normal operation of the MGB. As a consequence, this process cannot be automated. We started by removing from the dataset flights associated with system failures that were reported in the MIS database. Then, a sub-set of flights preceding each failure was also excluded from the dataset. The number of flights to exclude before each failure tagged in MIS depends on the type of the failure itself.

Data contextualization

As a matter of fact, the MGB system responds in different ways depending on usage and environmental conditions of the flight. Indeed, the systems such as the MGB are followed differently according to the different phases of the flight such as take-offs, hovers, turns, climbs, descents or landings and therefore react differently. In order to reduce the variance of the virtual sensor model and to maximize the chances of finding a representative model of the normal operation of the MGB system, an intelligent filtering of each time series based on expert knowledge has been applied. Thus, only subsequences of each flight have been selected. These subsequences, less influenced by the climatic conditions and the usage of a helicopter, correspond to the level flight phase, where the helicopter flies almost at the same speed, same altitude and same direction (stabilized flight). Other contextualisation strategies (expert-based), more restrictive on the flight behaviour, have been tested in order to obtain the more accurate model possible. However, they were leading to very small stable sequences of the flight which were not often repeatable over the flights history. For these reasons, it was decided to contextualize the training data relying on the flight sequences corresponding to level flight phase (long and repeatable sequences). This work on contextualization allowed analysing each flight according to the similar context. Consequently, we reduced some parameter variability.

Data sampling

Furthermore, to perform the splitting of flight history data, a sampling method has been developed to generate the training and a test dataset. The creation of these datasets is one of the most important steps, because it directly impacts the model's performance and robustness. To provide a representative dataset for all flights performed by AH customers (regardless of the geographical operation area, the weather conditions, the season of the year and the type of mission); a methodology of sampling flight data has been set up also for the creation of the test dataset. In this phase, we sampled the complete flight history performed by all AH customers. Consequently, we have carried out this step with a random sampling based on periodicity, geographical and mission type criteria.

C. Features selection

The features selection of the virtual sensor model has been conducted with the collaboration of experts. Initially, this made it possible to considerably reduce the scope of the analysis and allowed the obvious relationships from the expert's point of view to be verified through the operational data. The initial motivation for this study was to check the physical laws regarding this system and understand relationships between different contextual parameters.

D. Regression with Gradient Boosting

As stated above, the goal is to define a virtual sensor model (regression) from flight data set in normal operation. Then, the algorithm presented here uses this model to detect the moment that the operating system deviates from the normality. The regression was modelled with Gradient Boosting. The Boosting gradient introduced by Friedman [START_REF] Friedman | Greedy Function Approximation A Gradient Boosting Machine[END_REF] is a generalization of Boosting in which a loss function is used in a manner that can be assimilated to a gradient descent. The trees are thus added by reinforcement to the model to minimize the loss function.

The Gradient Boosting regression model has performed well on flight data compared to other usual regression models. An exploratory analysis has been conducted to test the performance of Gradient Boosting in comparison with other regression models. Table 2 represents the results obtained for Random Forest [START_REF] Breiman | Random forests[END_REF], linear regression [START_REF] Dg | Applied Regression Analysis and Multivariable Methods[END_REF]

E. Context influence on model performance

To complete the exploratory analysis, different data samples have been tested in similar ways according to a particular context. The graph shown in the figure below (Figure 2), illustrates the results obtained for different helicopters. Each bar represents the distribution of the R² observed on 5 samples analysed from different training and test datasets. We observe in all the different helicopters that the R² corresponding to the Gradient Boosting (in green) is higher than the other models. Only two exceptions represented by helicopters 6 and 9 do not validate this observation. These two helicopters obtain low R² values for all models. One of the explanations for this anomaly is the change in the geographical area of intervention of these two helicopters. For example, we have observed changes in operating areas from Asia to Northern Europe.

F. General methodology

The methodology implemented is fairly classical but each step has been developed in depth according to business application.

As illustrated in the Figure 4, the first step is dedicated to the creation of a database for the learning stage of the model and a second for the test step. Gradient Boosting is used on a subset of the data to form the learning stage database. The performance indicators are then evaluated on the test dataset. We use the coefficient of determination and mean squared error to evaluate the best regression model. Once the model performance evaluation step is completed, an anomaly score is calculated from the results obtained with the test dataset. In order to distinguish flights showing signs of failure from healthy flights, a measurement way has been set up to evaluate the evolution of the results returned by the prediction of the model (virtual sensor) and the values of the real sensor. This measurement, called here "anomaly score", allows us to obtain the distribution of the prediction errors observed within the learning base. For this purpose, we compute the absolute difference between the value predicted by the model and the measured value. The analysis of this distribution is necessary to introduce a threshold that will indicate a divergence between the real sensor and the virtual sensor.

The following figures (Figure 5 and Figure 6) illustrate results obtained with a Gradient Boosting regression for a normal flight and an abnormal flight. In both examples the oil pressure values are represented in blue, the values predicted by the Gradient Boosting in orange and the anomaly score in green. We can easily observe a difference in the distance between the curves representing the virtual sensor and the real sensor depending on the case of a normal or abnormal flight. We have therefore decided to use the anomaly score as a parameter that will indicate whether a flight is normal or abnormal. Thus, when we are in a normal flight situation, the anomaly score is close to zero and when we move away from this value, the flight has a tendency to deviate from normality.

(blue line), virtual sensor (orange line) and anomaly score (red line) for an abnormal

Then, a threshold on this score allows us to detect the flights coming out of the normal operating criteria. In this thresholding step, the "false alarms" question is the main constraint for monitoring tools in the aerospace industry due to the costs they induce.

III. ANOMALY DETECTION

A. Principle

The anomaly detection based on a predictive model consists in identifying atypical observations compared to a reference model. According to [START_REF] Chandola | Anomaly detection : a survey[END_REF], an anomaly is an observation that does not conform to the expected standard pattern. The threshold set on anomaly scores calculated from the historical data allows for the model to be validated using samples of real cases regarded as normal.

B. Threshold

To set an anomaly detection threshold on the base of the distribution of the anomaly score, the principle of Extreme Values Theory (EVT) was applied. Particularly in recent years, a growing interest has been observed for the application of Extreme Values Theory in the modelling of such events. For a fairly complete presentation of the topic, let us refer to the reference work of Embrechts, Klüppelberg and Mikosch [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF] giving the main theoretical results on the EVT and to Reiss and Thomas [START_REF] Reiss | Statistical Analysis of Extreme Value with Applications to Assurance, Finance, Hydrology and Other Fields[END_REF] who proposed a number of practical examples in finance, insurance and environment domains. In this section we are interested in the anomaly score distribution and more precisely the distribution tail of extreme values. The detection method implemented in this article is based on the notion of extreme quantiles. In order to approximate the lowest frequencies representing rare events and therefore abnormal occurrences, we use the classical techniques of EVT. We have used the Generalized Pareto Distribution (GPD) [22][25] to fit the distribution of the frequencies furthest from the graph below (Figure 7). The generalized Pareto distribution is a useful tool in this context because, in addition to accurately reproducing the shape of the empirical distribution above a threshold, the GPD makes it possible to estimate the high quantiles of the distribution. This means that the data can be extrapolated beyond the available sample.

Figure 7: Distribution of the anomaly score calculated on the test dataset. The x-axis represents the values of the absolute difference between the virtual and real sensors and the y-axis corresponds to their frequency. C. Experimentations

The use of the virtual sensor to detect abnormal events in the operation of the MGB is essentially based on the hypotheses that the training data represent all normal operating conditions of the helicopter.

The graphical representation below on Figure 8 represents the anomaly score calculated from a set of flight points and the threshold defined to detect the atypical points. The time stamp is represented on the x-axis of the graph and the anomaly score values on the y-axis. The anomaly score is displayed as negative values to be in conformity with the usage of Scikit-Learn python library. Figure 8 displays the results associated with the anomaly score calculated for each time stamp of the filtered flight sequence on level flight phase. The detection threshold for an abnormal event is set for the complete fleet and represented in red dotted lines. To initialize the anomaly detection step for all helicopters in this study, the threshold was set at an absolute difference of 1.2 bar pressure between the values of the real sensor and the virtual sensor. This threshold has been also validated by the domain expert. Testing with this threshold level is quite "conservative", i.e. it voluntarily detects events more rarely. On the other hand, the detection model has been evaluated on cases identified as abnormal during the data preparation stage and more precisely, the data cleaning.

Following this data cleaning work, different forms of detections have appeared. For these abnormal cases, we have identified, for example, anomalies with a relatively slow evolution, observed on several consecutive flights before a maintenance operation. We have also identified anomalies impacting only a few flights before the maintenance intervention with a special pattern for the anomaly score.

IV. MODEL VALIDATION

Data sampling work for the creation of test and training datasets via multiple temporal train-test splits and the application of cross-validation made it possible to check the learning capabilities of the Gradient Boosting model.

In addition, to validate the model's detection performance, different cases of anomalies have been selected by domain experts. These are cases of failures that cannot be detected with vibration analysis. We are particularly interested in anomalies produced by: lubrication system, oil pump and radiator filter. Three flights containing anomalies of this kind are represented in the Figure 9 3 . All these cases are correctly identified by our algorithm.

Several flights considered as normal flights are also proposed by domain experts. All cases are identified as normal by our system (Figure 10), except one. After discussion with domain experts, this case seems to have a real anomaly, probably detected by the maintenance team later after the flight. TESTS WITH MIS DATA

Unused flight phases

Level flight phases

The illustration of an anomaly detection case is shown in Figure 11. The three scatterplots correspond respectively to the anomaly scores calculated on flights before the occurrence of a failure, with the failure and after a corrective maintenance operation. The first figure shows the flight performed by a helicopter shortly before the change of the MGB system. This operation was completed as part of AH's scheduled maintenance plan, which involves changing the system after a number of flight hours performed by the helicopter (TBO: Time Between Overhauls). The second figure shows the first flight operated after the same maintenance operation by the same helicopter. The detection of an abnormal event is obvious. most of the anomaly score points from this flight are below the threshold. The third graphical representation displays the flight performed after the new change of the MGB after an on-board alarm has been triggered. This helicopter had been flying for one month after the first maintenance operation and had not triggered an on-board alarm for exceeding the critical damage threshold of the system. Although the severity of this failure allowed the helicopter's missions have been performed without incident, the failure propagation was damaging the MGB system and required the removal of the system only one month after its installation. In this example, we were able to detect the anomaly from the first flight operated by the helicopter.

The proposed detection method can thus prevent a malfunction of the monitored system and can prevent its advanced degradation as an expensive system such as the MGB. In this example, threshold detection is easily interpretable and effective. However, the abnormal flights that are detected by this method involving the virtual sensor may have different patterns. Indeed, some cases of detection are not as obvious as the latter example and require analysis to identify the phenomena explaining this variety of patterns.

The two graphs on Figure 12 show other cases of abnormal flight detection which are less obvious than the previous example. Indeed, in both cases, the points representing the anomaly scores calculated on the flights do not all exceed the limit threshold. However, the point distributions in both graphs are very distinctive. For each of these patterns, the following flights have a similar distribution. In some cases, we were able to identify a maintenance operation that allowed us to observe a return to normal system operation. This is the case for the first pattern described by a very noisy distribution of the scatterplot. The second pattern which is observed to follow a trend does not necessarily allow associating this pattern to a type of failure in the MGB system.

VI. DEPLOYMENT WITHIN AIRBUS HELICOPTERS FLYSCAN SERVICE

Since 2007, HUMS / flight data has been collected from various aircrafts during different missions. Since 2012, all data has been automatically treated by HUMS information systems in order to return validated analysis results back to customers via a web portal. For now, valorised results are related to statistics, vibration indicators, KPI, flight spectrum, etc. In addition, HUMS information systems now aims to add predictive maintenance results to their capability. Airbus Helicopters aims to be the first to propose predictive services to customers in order to anticipate issues impacting helicopter availability.

Predictive services consist in monitoring the trend of the results computed by HUMS / Flight data algorithms deployed in the customer oriented data chain. An alerting based on that monitoring is added in order to allow the support team to validate the alert and eventually advise the concerned customer on the maintenance action to perform (Figure 13).

Figure 13: The integration of the virtual sensor in an AH service A. Virtual Sensor Integration in HUMS / FlyScan Industrial

Process The virtual sensor presented here is a way to achieve the failure anticipation industrial objective in an aeronautical context. The virtual sensor has been deployed recently in FlyScan. With the predictive model presented, it is possible to detect abnormal usage of the helicopter and predict some failures. However, when an alert is raised no decision will be made automatically. For the moment this step is under the responsibility of a dedicated support team, able to understand the alert, troubleshoot the issue and provide advice to customers in terms of maintenance actions. Nevertheless, on-going research projects are exploring methodologies to make decisions more rapid and optimal in the future (Figure 14). FlyScan supports operators who have subscribed to the service. Thus, Airbus Helicopters brings its expertise to customers, by monitoring their helicopter fleets based on collected data. This service aims to maximize the customer's operational efficiency through data analytics, enhancing flight safety and increasing fleet availability (eg. decreasing the number of Aircraft On Ground events per aircraft per year).

FlyScan is part of HCare Technical Support 4 , offering 24/7 technical support service, continuous airworthiness management and global technical assistance. An operational control room is able to display all events detected by predictive algorithms (Figure 15). The alerted support team analyses the event and validates the relevance of the detection. Algorithms results combined with experts' knowledge about the monitored systems enable the right decision to be made on whether or not to trigger a maintenance action at the customer end.

VII. DISCUSSION

The aim of the proposed method is to automatically detect malfunctions in the helicopter MGB, by exploiting operational data and business knowledge provided by experts.

The principle of the virtual sensor, in the context of anomaly detection, is based on the use of parameters recorded during flights to model an already existing sensor. The presented MGB oil sensor model, which is based on a Machine Learning approach, allows simulating the oil pressure values from other parameters monitored during flight. Subsequently, an anomaly score is developed, which consists in measuring the difference between the virtual sensor values and the measured ones during flight. Finally, a threshold has been set to allow the detection of abnormal flights.

The method could be improved with respect to several aspects. First, the learning data base cleaning is based only on MIS data, which could be incomplete since it is not mandatory for customers to record/complete it. Another option could be the use of unsupervised machine learning in order to automatically isolate faulty cases. Moreover, data cleaning was a tedious and long task. In fact, MIS data is textual inputs from customers' technicians for describing their maintenance actions on helicopters. It is often difficult for a non expert to understand all technical terms used, especially when some terms are interchangeable. To use the data we needed the help from an expert in order to identify maintenance actions as well as key terms referring to MGB sub-systems' failures. Thus, use of this data was not optimal and not reusable. Moreover, the process was long due to unavailability of the expert. Actually, Airbus Helicopters team responsible of MIS data is developing a natural language and text processing based tool for better use and more automated analysis of the data. This will allow for a better and easier cleaning of our learning database, and consequently results presented here to be improved.

The learning phase is based on the Gradient Boosting algorithm which is recognized and confirmed as a efficient machine learning algorithm compared to classical supervised algorithms. This step, however, could be improved by exploring the latest progress in deep learning. In fact, the presented learning problem could also be modeled using deep learning, based for example on the Long Short Term Memory (LSTM) architecture.

The flight phase actually considered is the one specified by the expert, namely the level flight phase. It could be interesting to use a data analysis approach to identify the best and acceptable flight phase to consider which may improve the results, and where different contextual parameters have less variation.

The influence of the geographical area of operations on the MGB oil pressure requires considering the changes of missions in the learning process. Actually, the solutions considered to deal with these issues would be to integrate helicopters location parameters such as geographical positions among the input parameters or to trigger the relearning of the model when a change in geographical area is identified.

During tests performed on real failures tagged in MIS, we have observed different patterns when a flight is detected as abnormal. Sometimes dispersion is observed (Figure 12), some other times ascending trends are also observed. Actually these different shapes are not analysed in depth, but would be useful in order to study the possibility of enlarging anticipation intervals based on these shapes.

We have chosen to develop the same detection model for the complete fleet of a helicopter type by integrating the identifier of the aircraft into the input variables. This choice is justifiable because of the hypothesis that the systems are designed identically and therefore are supposed to respond in the same way depending on the use of the helicopter. However, operations may differ from one customer to another and the model may be less effective for less standard helicopter use. This problem could be solved by duplicating the detection model for each helicopter.

Finally, features (contextual parameters) considered in the implemented model were selected according to the expert knowledge. This step could be automated using recent works about features selection that maximize the performance of a Gradient Boosting regression [START_REF] Friedman | Greedy Function Approximation A Gradient Boosting Machine[END_REF]. Further experimentation of automatic features selection would be likely to identify additional parameters in order to model the oil pressure sensor more accurately. Moreover, this would also allow for the generalization of this method to other sensors and other helicopter systems.

VIII. CONCLUSION

This article presents the end-to-end solution for detecting anomalies in industrial preventive maintenance. It addresses the different steps and difficulties in the processing of huge amounts of data for the development of a failure anticipation method. The initial objective of defining the normal operation of a helicopter system was achieved through the cooperation of aeronautical experts. This expert knowledge has been used to exploit HUMS flight data with machine learning methods to model a form of normal operation of the MGB system through the virtual oil pressure sensor. This is a pragmatic solution that allows AH customer support to improve the maintenance service offered to customers.

We are currently working on the drift monitoring of the deployed virtual sensor, in order to trigger its update when new and different data is collected.

  [START_REF] Sa | Primer of Applied Regression and Analysis of Variance[END_REF][START_REF] Mh | Applied Linear Statistical Models[END_REF] and Gradient Boosting on our dataset, using training datasets of different sizes. The two comparison indicators are the coefficient of determination (R²) and the mean squared error (MSE). The ideal values are R 2 = 1 and MSE = 0. For these two indicators, the Gradient Boosting obtains the highest scores. This is reflected both in the accuracy of the estimate with lowest MSE values regardless of the size of the data sample, and in the prediction quality of the model with the R² always the highest.

Figure 2 :

 2 Figure 2: Boxplot represents R² calculated for 10 different aircrafts. (Random Forest: blue; Gradient Boosting: orange; Linear Regression: green)In the same type of analysis, we have compared the Gradient Boosting with the Random Forest as well as the linear regression on different data samples according to the season. We reproduced 5 different samples for winter and summer in the boxplot Figure3. With respect to other models, Gradients Boosting obtains good results on the R² regardless of the season.

Figure 3 :

 3 Figure 3: Boxplot represents MSE calculated for summer and winter. (Summer: blue; winter: orange)

Figure 5 :

 5 Figure 5: Oil pressure sensor (blue line), virtual sensor (orange line) and anomaly score (red line) for a normal flight

Figure 6 :

 6 Figure 6: Oil pressure sensor (blue line), virtual sensor (orange line) and anomaly score (red line) for an abnormal

Figure 4 :

 4 Figure 4: General methodology for the virtual sensor

Figure 8 :

 8 Figure 8: Graphic with the anomaly score and the threshold for a normal flight filtered in level flight phases. The x-axis represents the time-step and the y-axis represents the anomaly score measured in bars (negative value).

Figure 9 :Figure 10 :

 910 Figure 9: Oil pressure sensor (blue line) and virtual sensor (orange line) for three normal flights

Figure 11 :

 11 Figure 11: Anomaly detection starting from the first flight operated after the removal of the MGB system. Anomaly Score is represented by the blue points and the threshold by the horizontal red dashed line.

Figure 12 :

 12 Figure 12: Different abnormal flight pattern. The anomaly score is represented on the y-axis and the time step for each flight on the xaxis.

Figure 14 :

 14 Figure 14: From raw data to autonomous decision making B. Operational HUMS Support RoomFlyScan supports operators who have subscribed to the service. Thus, Airbus Helicopters brings its expertise to customers, by monitoring their helicopter fleets based on collected data. This service aims to maximize the customer's operational efficiency through data analytics, enhancing flight safety and increasing fleet availability (eg. decreasing the number of Aircraft On Ground events per aircraft per year).

Figure 15 :

 15 Figure 15: FlyScan Interface This human step allows a validation of the algorithm output but it is also a means to improve and adjust deployed predictive models. In fact, this allows a labelling of the case as true or false positive. Then, each detected case, rejected or validated by the support team, reinforces the model. Design office

TABLE I .

 I TIME SERIES COMPOSING THE USED DATASET. THE NAME OF

				ATTRIBUTS IS ANONYMISED
	Flight id	time step	f0	f1	f2	f3	…	fp
	1	1

TABLE II .

 II PERFORMANCES OF DIFFERENT REGRESSION MODELS. FOR

EACH TRAINING DATASET SIZE THE EXPERIMENT HAS BEEN REPEATED 10 TIMES WITH DIFFERENT RANDOM DATA RESPECTING OUR CONSTRAINTS AND THE AVERAGE VALUE IS REPORTED IN THE TABLE.

Please note that for confidentiality reasons, the real values of MGB oil pressure are anonymized.

Acknowledgments:

We would like to thank our colleagues Jayant SEN-GUPTA, Olivier REGNIER-COUDERT and Nolwen HUET for their support, and the interesting technical exchanges about this work; done in the framework of the AIRBUS ADVISED project about anomaly detection in time series data.