Ruixin Wang
email: ruixin.wang@recherche.enac.fr

Cyril Allignol
email: cyril.allignol@enac.fr

Nicolas Barnier
email: nicolas.barnier@enac.fr

Alexandre Gondran
email: alexandre.gondran@enac.fr

Jean-Baptiste Gotteland

Catherine Mancel
email: catherine.mancel@enac.fr

A New Multi-commodity Flow Model to Optimize the Robustness of the Gate Allocation Problem

Keywords: Gate Allocation Problem, Robustness, Multi-commodity Flow Problem, Integer Linear Programming

The Gate Allocation Problem (GAP) is one of the numerous operational problems that all busy airports have to handle and to optimize every day. In this article, we focus on two important objectives for the GAP: the robustness and the taxiing times.We first analyze a basic model which optimizes the robustness of the solution to the GAP and propose a more realistic and compact Multi-commodity Flow Problem model (MFP) to optimize both the robustness and the taxiing times. Based on tests with real data at Paris-Charles-de-Gaulle airport, this new efficient MFP model outperforms the basic one by orders of magnitude.

Introduction

All airports have to handle many challenges regarding traffic increase, delay reduction, traffic predictability and environmental impact. In this context, airports attempt to improve the efficiency of their operations, among which gate allocation is one of the most important [START_REF] Bouras | The airport gate assignment problem: A survey[END_REF]. Given a set of flights, the schedule of their gate occupancy (i.e. their respective arrival and departure times at the gate) and a list of compatible gates, the Gate Allocation Problem (GAP) consists in assigning a gate to each flight such that at most one aircraft occupies a gate at any time, while optimizing various operating costs.

Classic objectives of the GAP include the optimization of passengers walking distance [START_REF] Kim | Airport Control Through Intelligent Gate Assignment[END_REF], the allocation of terminal gates over apron stands or the minimization of the number of towing movements [START_REF] Guépet | Exact and heuristic approaches to the airport stand allocation problem[END_REF]. However, the emergence of Aircraft Collaborative

Literature Review

The Gate Allocation Problem1 (GAP) consists in assigning arriving flights with fixed occupancy periods to available compatible gates while maximizing both conveniences to passengers and operational efficiency of airports [START_REF] Bouras | The airport gate assignment problem: a survey[END_REF]. [START_REF] Steuart | Gate position requirements at metropolitan airports[END_REF] was one of the first to introduce the problem in 1974, but the literature is scarce until the 2000s. However, many variants of the GAP have been studied since, as mentioned in [START_REF] Das | A review on airport gate assignment problems: Single versus multi objective approaches[END_REF].

If there were no compatibility restriction between gates and aircraft, the corresponding decision problem could be modeled as the coloring of an interval graph, which is polynomial [START_REF] Gupta | Efficient algorithms for interval graphs and circular-arc graphs[END_REF] for the minimization of the number of colors (i.e. gates). Airport gates are generally not equivalent resources though, as they are dimensioned to accommodate specific types of aircraft. Therefore, the set of compatible gates for an aircraft is usually a strict subset of all the available gates and the decision version of the allocation problem is rather a list-coloring problem, which is NP-Complete even for interval graphs [START_REF] Biró | Precoloring extension. I. Interval graphs[END_REF].

Moreover, gates may also be endowed with other secondary features (e.g. compatible airlines, domestic or international, terminal or apron, etc.) which should match the characteristics of the flight and the preferences of airlines as much as possible. These preferences can be modeled as costs associated with each possible assignment, and standard GAP objectives often include the minimization of their sum, which is NP-Hard [START_REF] Deng | Study on an improved adaptive PSO algorithm for solving multi-objective gate assignment[END_REF]. From the airport operations perspective, the objective could be to maximize the utilization of the available gates and terminal [START_REF] Steuart | Gate position requirements at metropolitan airports[END_REF][START_REF] Yan | Optimization of multiple objective gate assignments[END_REF][START_REF] Li | Optimized assignment of airport gate configuration based on immune genetic algorithm[END_REF], minimize the flight delay [START_REF] Yan | A heuristic approach for airport gate assignments for stochastic flight delays[END_REF] and maximizing the preferences [START_REF] Dorndorf | Flight gate scheduling: State-of-the-art and recent developments[END_REF] (i.e. certain aircraft should go for particular gates). Other classic objectives can be the passengers walking distance [START_REF] Das | New multi objective models for the gate assignment problem[END_REF][START_REF] Dijk | The recoverable robust stand allocation problem: a GRU airport case study[END_REF][START_REF] Mokhtarimousavi | A non-dominated sorting genetic algorithm approach for optimization of multi-objective airport gate assignment problem[END_REF] (or other connection means like buses), which is similar to the Quadratic Assignment Problem, or the number of towing movements [START_REF] Guépet | Exact and heuristic approaches to the airport stand allocation problem[END_REF].

Even if a flight occupancy is fixed, many real-life factors of uncertainty (e.g. traffic delays, severe weather conditions, equipment failure, etc.), can lead to deviations from the original schedule. If there is not enough buffer time between successive flight occupancies at a given gate, a delay may propagate to other flights, then to other gates, and hinder the operational efficiency of the whole terminal. In order to absorb potential delays and avoid costly disruptions, our study rather focuses on optimizing the robustness of the allocation as proposed by [START_REF] Bolat | Procedures for providing robust gate assignments for arriving aircrafts[END_REF], which minimizes the variance of idle times of gates to balance and spread them over time and resources. Despite its practical importance, research on the robustness of solutions to the GAP is quite limited.

To solve these very diverse variants of the GAP, many classic combinatorial optimization methods were experimented, depending on the linearity of the model, the size of the instances and the requirements on the execution time of the solver. One of the most used tools is Mixed Integer Programming or Integer Linear Programming (ILP) solvers, like Gurobi or CPLEX, to obtain proven optimal solution like in [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF][START_REF] Guépet | Exact and heuristic approaches to the airport stand allocation problem[END_REF]. Several studies present their ILP model as a Multi-commodity Flow Problem (MFP), like [START_REF] Zhang | Optimization for gate re-assignment[END_REF] which models re-assignment while minimizing passengers cost. Another complete combinatorial optimization technique, Constraint Programming, was also experimented in [START_REF] Simonis | Models for global constraint applications[END_REF][START_REF] Wang | Propagation of idle times costs for fixed job scheduling[END_REF] to solve the GAP as a scheduling problem similar to Fixed Job Scheduling (FJS).

As previously mentioned, all considered variants of the GAP are NP-Complete or NP-Hard, so many publications have suggested various heuristic approaches for solving large instances or non-linear models in reasonable time. [START_REF] Yan | A network model for gate assignment[END_REF] developed an algorithm based on the Lagrangian relaxation of the GAP, with subgradient methods, accompanied by a shortest path algorithm and a Lagrangian heuristic to solve an MFP model. Additionally to the ILP model previously mentioned, [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF] also proposed a Genetic Algorithm, which is a population-based meta-heuristic, to provide several alternative solutions in case of schedule deviations exceeding the robustness of the implemented solution. [START_REF] Benlic | Breakout local search for the multi-objective gate allocation problem[END_REF] proposed a multi-objective heuristic approach base on Breakout Local Search, with a particular focus on the perturbation strategy. [START_REF] Deng | Study on an improved adaptive PSO algorithm for solving multi-objective gate assignment[END_REF] proposed an improved Particle Swarm Optimization algorithm to solve a multi-objective (walking distance of passengers, number of flights at parking apron, etc.) optimization model. [START_REF] Zhang | Optimization for gate re-assignment[END_REF] used two heuristic algorithms to minimize the weighted sum of the total flight delays. More research work can be found in a recent survey about the GAP [START_REF] Das | A review on airport gate assignment problems: Single versus multi objective approaches[END_REF].

More generally, there are already many studies on the GAP focusing on different criteria. However, few of them except [START_REF] Bolat | Procedures for providing robust gate assignments for arriving aircrafts[END_REF][START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF] really consider the robustness of the assignment, all the more while minimizing taxiing time simultaneously, which appears to be one of the most important criteria for major international airports like Paris-CDG. Moreover, the resolution process should be efficient enough to be computed in a few minutes, as the GAP does not occur only once (one day before the date of the traffic), but has also to be solved almost in real time when a severe disruption leads to a necessary re-assignment of gates. Using the same decision variables as [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF], which does not identify the GAP as a flow problem, we introduce in this article a much more efficient and compact MFP model to optimize its robustness. For an instance with n flights and m gates, we reformulate the original O(mn 2 ) "continuity" constraints that ensure the feasibility of the assignment, by only O(mn) flow conservation constraints. This tighter model drastically reduces resolution times and allows to optimally solve large-scale instances. Furthermore, we also propose to take taxiing costs, a major concern of large international airports, into account in the objective function, as well as consider subsets of mutually exclusive gates (which occur at Paris-CDG) with additional constraints.

The Gate Allocation Problem

After presenting the operational context of the GAP at Paris-CDG and the importance of optimizing the robustness and taxiing times of solutions, we formalize the combinatorial optimization problem independently from its resolution technique and specific model, thoroughly describing its instances components (except mutually dependent gates which are presented in Section 4.4), its feasibility condition and possible objective costs. Eventually, we define a compatibility predicate to help define compact ILP models in Sections 3 and 4.

Operational Context

The Gate Allocation Problem (GAP) is one of the numerous operational problems that all busy airports have to handle and to optimize every day. It is usually tackled with specific usages and preferences that are difficult to model or list rigorously -some airlines may even manage their own stands area with their own undisclosed strategy. This can lead to a very complex organization, in which the parking stands can have various configurations during each period of the day. Some groups of gates can also be dependent (cf. Section 4.4): only one aircraft can occupy one of them at a time, but these gates cannot be considered as a single one, because they have different entry or exit points, resulting in different taxiing times for the aircraft. Therefore, the GAP is often specific to the considered airport and can be formulated very differently, with various degrees of accuracy.

In the scope of this article, we consider the GAP at Paris-CDG, one of the busiest European international airports. The terminals and their gates (including dependent ones) are well described, and the types of aircraft that can operationally use a given gate can be easily reconstructed from the available actual traffic records. Figures 1 and2 are extracted from the official airport charts and show the location of the gates in the main terminals at Paris-CDG. Figure 1 shows how some gates of different terminals can be dependent at Paris-CDG 1, as they cannot be used at the same time: on this figure, arrows represent the position of the aircraft at each stand and, for example, gates V08 and W02 are mutually exclusive. For this reason, instances of the GAP that involve these terminals cannot be solved separately and have to be merged in a single, larger one, even if the precise terminal of each flight is fixed.

Moreover, Paris-CDG has implemented the Airport Collaborative Decision Making (A-CDM) program during the last decade [15]: this program defines which and how accurate information (arrival and departure times, aircraft delay. . . ) can be shared between the different airport stakeholders [14], in order to help them make more efficient decisions in real time. Not only Paris-CDG, many other European and Chinese airports have also obtained A-CDM certification and are currently in the process of its implementation [15,[START_REF]China airport ACDM implementation route[END_REF]. The A-CDM program also includes a major improvement in the departures management [START_REF] Huet | Local and network impact assessment of airport collaborative decision making (a-cdm)[END_REF]: the delay due to the runways capacities is anticipated and aircraft are preferably delayed at gate, engines off, rather than near the runway after start-up, burning fuel and emitting CO 2 . The benefits are twofold: it results in less taxiing traffic, which decreases congestion and nuisance, and also provides a more accurate departure schedule. However, it has a negative effect on gates occupancy, which can significantly increase, especially during peak hours. For these reasons, the robustness of the gate allocation towards additional gate occupancy due to departures delay becomes more and more important for such airports: significant disruptions can be caused in case of gate conflicts between arriving aircraft and delayed departures.

Another important criterion for the GAP is the resulting taxiing times: especially at big airport like Paris-CDG, taxiing times minimization is one of the main objectives for both airlines and airport managers. However, many other constraints like safety and customs services, airlines organization, passengers information, etc. prevent to be able to change the terminal of a flight for this purpose, as it is generally fixed a long time before the day of traffic. Nevertheless, some big terminals offer several distant gates at Paris-CDG, and have different ways to enter or leave, so that taxiing times can significantly vary from a gate to another in the same terminal and must be minimized as well, as shown in Section 4.3.

Instance Definition

An instance of the basic model of the GAP is defined by:

• F = { f 1 , .
. . , f n }: a set of n flights (or aircraft) to arrive at the airport during the considered planning horizon and to be assigned to gates.

• G = {g 1 , . . . , g m }: a set of m gates available during the considered planning horizon.

• ∀ f i ∈ F :

f s i and f e i : the expected arrival and departing times of flight f i to its assigned gate, with f s i f e i . Note that these arrival and departing times are supposed to be known in advance and not to depend on the gate allocation decisions. Moreover, we suppose that the flights, indexed from 1 to n (noted [1 . . n]), are ordered by increasing arrival time:

∀ f i , f j ∈ F , i < j ⇒ f s i f s j -G i
⊆ G: a subset that specifies the possible gates on which flight f i can be assigned, depending on the type of the aircraft and the characteristics of the gates. We suppose that G i = ∅, otherwise there is no solution.

• ∀g k ∈ G:

g s k and g e k : the opening and closing times of gate g k . Without loss of generality, we can consider that all gates have the same availability [g s . . g e ], i.e. ∀g k ∈ G, g s k = g s and g e k = g e ; to represent a period of unavailability of gate g k between t 1 and t 2 (e.g. for maintenance or while understaffed), we just add an extra flight f i to F with f s i = t 1 , f e i = t 2 and G i = {g k }. We also consider that ∀ f i ∈ F , g s f s i and f e i g e , i.e. all flights are scheduled during the period of availability of the gates [g s . . g e ] (otherwise there is trivially no possible allocation).

-F k ⊆ F : the subset of flights which can be assigned to gate g k , redundantly defined as

F k = { f i ∈ F | g k ∈ G i } to simplify notations in the following sections.
Subsets G i (and F k ) generally have to be built from raw data that specify the type of aircraft for each flight (among a given set of types) and the accepted types for each gate. Other characteristics of the flights and the gates can further restrict the possible gates, e.g. the airline, custom or services requirement. . . In Section 4.4, we refine this model by taking into account mutually exclusive subsets of gates.

Feasibility

An instance of the GAP is feasible if and only if there exists an assignment σ of flights to gates such that each flight is assigned a gate from its possible ones (Equation ( 1)) and flights assigned to the same gate do not overlap (Equation (2)), i.e. ∃σ : F → G such that:

∀ f i ∈ F , σ( f i ) ∈ G i (1)
∀ f i , f j ∈ F s.t. i < j, σ( f i ) = σ( f j ) ⇒ f e i f s j (2)

Cost

As already mentioned, various kinds of cost can be associated to a feasible allocation. The most usual one is the sum of the individual costs c a i,k associated to each assignment σ( f i ) = g k : Definition 1 (Assignment Cost). The assignment cost of a solution to an instance of the GAP is defined by:

cost a = n i=1 c a i,k
However, we focus in this article on the robustness of the allocation w.r.t. possible time deviations of the gate occupancies. As proposed by [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF], the robustness can be optimized by minimizing the variance of the idle times of the gates, such that they are distributed as evenly as possible over the schedule, allowing occupancies to exceed their expected time bounds without incurring gate conflicts. As the overall leeway and the number of idle times are constant (as explained in Section 3), the mean of the idle times is fixed and the minimization of their variance amounts to the minimization of the sum of their squares.

More generally, we can model an idle time as the cost c tr i, j of the transition between two successive flights f i and f j assigned to the same gate, considered to be independent of the gate. The transition cost also takes into account costs c tr 0,i of the transition between the opening of each gate, represented by fictive flight f 0 , and its first flight f i , and c tr . n] because flights are ordered by arrival time and therefore i < succ(i)) of the flight that immediately succeeds f i on its assigned gate, or n + 1 if it is the last flight.

Definition 3 (Idle Time Duration). The idle time duration I i, j between two successive flights f i ∈ { f 0 } ∪ F and f j ∈ F ∪ { f n+1 } assigned to the same gate is defined by:

∀i ∈ [0 . . n], ∀ j ∈ [i + 1 . . n + 1], I i, j =                    f s j -f e i if i 1 ∧ j n (successive flights) f s j -g s if i = 0 ∧ j n (first flight) g e -f e i if i 1 ∧ j = n + 1 (last flight) g e -g s otherwise (i = 0 ∧ j = n + 1, empty gate)
To optimize the robustness of a solution, the transition cost is defined in our model as c tr i, j = I 2 i, j (as explained in Section 3), i.e. the square of the duration of the idle time between two successive flights f i and f j assigned to the same gate (possibly of null duration if succ(i) = j and f e i = f s j ), or between the opening of a gate and its first flight, or between the last flight of a gate and its closing, or between the opening and closing of an unoccupied gate:

Note that even if this abstract model considers uniform gates, i.e. the transition cost is independent from the gate, the ILP model described in the next section can natively represent gate-dependent costs.

Flights and Gates Compatibility

For a given instance, not all the transition costs c tr i, j , ∀i, j ∈ [0 . . n + 1], i < j, need to be defined, but only for the ordered pairs of indices of compatible flights, i.e. that have at least a possible gate in common and do not overlap. We define a compatibility predicate to specify the corresponding compatible triplets of indices (of two possibly successive flights and a given gate) and help the writing of our ILP models in the following sections, which discard unnecessary decision variables:

Definition 4 (Compatibility Predicate). The compatibility predicate γ : [0 . . n] × [1 . . n + 1] × [1 .
. m] → B on two ordered flight indices i and j and a gate index k is defined by:

∀i ∈ [0 . . n], ∀ j ∈ [i + 1 . . n + 1], ∀k ∈ [1 . . m], γ(i, j, k) =                    g k ∈ G i ∩ G j ∧ f e i f s j if i 1 ∧ j n (successive flights) g k ∈ G j if i = 0 ∧ j n (first flight) g k ∈ G i if i 1 ∧ j = n + 1 (last flight) true otherwise (i = 0 ∧ j = n + 1, empty gate)
f i , f j and g k are said to be compatible if γ(i, j, k) = true.

Basic Transition Cost ILP Model

In [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF], the author proposes several mathematical models for the GAP, with different kind of optimization criteria, in order to capture the robustness of gate assignment solutions towards flight schedule disruptions. In particular, he proposes to minimize the variance of the idle times of the gates to evenly distribute the global leeway over resources and time. This criterion allows for short or long pauses that might be necessary in some instances without hindering the optimization of the rest of the schedule An idle time is defined as a time period, possibly of null duration, between two successive flights assigned to the same gate and during which the gate is unoccupied. The idle times occurring between the opening of a gate and its first assigned flight, and between its last flight and its closing, are also taken into account as the gates are available to absorb possible delays (or early arrivals), such that there are exactly n + m idle times in any allocation: one after any of the n flights and one between the opening of any of the m gates and its first flight (or its closing for unoccupied ones). The underlying idea of this minimization criterion is that the more the flights are regularly distributed to gates, the less the corresponding gate allocation solution should be sensitive to schedule disruptions.

Indeed, even if the number and the sum of the idle times is considered to be constant in [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF], because the arriving and departing times of flights are fixed, their respective durations directly depend on the gate assignment. Hence, the more the idle times are balanced w.r.t. their durations, the smaller the sum of their squares (or any other convex function). The author proposes an Integer Linear Programming (ILP) model to solve the GAP with the sum of the squares of the idle times as the minimization criterion, with O(mn 2 ) binary variables x i, j,k stating that flight f j is the successor of flight f i on gate g k , and O(mn 2 ) constraints.

Before introducing our tighter and more compact flow model in Section 4, which only needs O(mn) constraints, we reformulate this ILP model in the following sections, taking into account the heterogeneity of the gates and discarding variables that can statically be set to 0 because of an incompatibility (see Section 2.5), which was not explicitly done in [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF].

Decision Variables

A binary variable is defined to denote the successive assignment of a pair of compatible flights to the same compatible gate, using compatibility predicate γ of Definition 4:

∀i ∈ [1 . . n], ∀ j ∈ [i + 1 . . n], ∀k ∈ [1 . . m], s.t. γ(i, j, k), x i, j,k =
1 if flights f i and f j are successively assigned to gate g k 0 otherwise

(3)

∀ j ∈ [1 . . n], ∀k ∈ [1 . . m], s.t. γ(0, j, k), x 0, j,k = 1 if f j is the first flight assigned to gate g k 0 otherwise (4) ∀i ∈ [1 . . n], ∀k ∈ [1 . . m], s.t. γ(i, n + 1, m), x i,n+1,k = 1 if f i is the last flight assigned to gateg k 0 otherwise (5) ∀k ∈ [1 . . m],
x 0,n+1,k = 1 if there is no flight assigned to gate g k 0 otherwise [START_REF] Bouras | The airport gate assignment problem: A survey[END_REF] If the instance is feasible, there is exactly one x i, j,k variable equal to 1 for each flight f i , so the allocation σ can trivially be deduced from the decision variables:

∀ f i ∈ F , ∃! j ∈ [i + 1 . . n], ∃! k ∈ [1 . . m], s.t. x i, j,k = 1, σ( f i ) = k
or, alternatively, with the linear expression:

σ( f i ) = j∈[i+1..n] k∈[1..m] γ(i, j,k) kx i, j,k
Note that variables x i, j,k are similar to the ones of typical ILP models of the Multicommodity Flow Problem [START_REF] Wang | Multicommodity network flows: A survey, part I: Applications and formulations[END_REF], where they are interpreted as the the amount of a commodity k flowing through an arc from node i to node j in a network. Nevertheless, as shown hereafter, the model proposed in [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF] is not formulated as a flow problem model.

Constraints

Covering constraints: each flight has to be assigned to exactly one compatible gate. In order to guarantee this statement, the following constraints are introduced to ensure that exactly one flight is assigned immediately before a given flight f j (Equation ( 7)), and only one flight is assigned immediately after a given flight f i (Equation ( 8)):

∀ j ∈ [1 . . n], i∈[0.. j-1] k∈[1..m] γ(i, j,k) x i, j,k = 1 (7) ∀i ∈ [1 . . n], j∈[i+1..n+1] k∈[1..m] γ(i, j,k) x i, j,k = 1 (8)
Moreover, Equation ( 9) is added to specify that there are exactly n + m binary variables that should be equal to 1 in a feasible solution:

i∈[0..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k) x i, j,k = n + m (9)
Continuity constraint: the successive assignments have to be feasible. Indeed, if a pair of flights ( f i , f j ) with j n (i.e. f i is not the last flight on its gate) are successively assigned to a gate g k , there must be a flight f j with j ∈ [ j + 1 . . n + 1] such that the pair of flights ( f j , f j ) are also successively assigned to g k , but there cannot be any flight f j such that the pair of flights ( f j , f j ) are successively assigned to a gate g k = g k . This continuity constraint is modeled with Equation ( 10), which states that, for any pair of flights ( f i , f j ) and any gate g k , either x i, j,k = 1, and then no x j, j ,k with j > j and k = k can be equal to 1 ( k ∈[1..m]|k =k, j ∈[ j+1..n+1] x j, j ,k = 0); either x i, j,k = 0, and then at most one x j, j ,k with j > j and k = k equals 1 (which is always true due to Equations ( 7) and ( 8)

). ∀i ∈ [0 . . n -1], ∀ j ∈ [i + 1 . . n], ∀k ∈ [1 . . m], γ(i, j, k), x i, j,k + j ∈[ j+1..n+1] k ∈[1..m]\{k} γ( j, j ,k )
x j, j ,k 1 [START_REF] Das | A review on airport gate assignment problems: Single versus multi objective approaches[END_REF] Note that these continuity constraints are not formulated as a flow balance constraint but as the linearization of the logical constraint:

∀i ∈ [0 . . n -1], ∀ j ∈ [i + 1 . . n], ∀k ∈ [1 . . m], γ(i, j, k), ∀ j ∈ [ j + 1 . . n + 1], ∀k ∈ [1 . . m] \ {k}, γ( j, j , k ),
x i, j,k = 1 ⇒ x j, j ,k = 0 These O(mn 2 ) constraints dominate the number of constraints of this model.

Resource constraints: each gate can be assigned to at most one flight at a time. Due to the previous continuity constraints, it is sufficient to check that each gate should have at most one assigned first flight. Therefore Equation ( 11) is defined by:

∀k ∈ [1 . . m], j∈[1..n+1] γ(0, j,k) x 0, j,k 1 (11)
Domain constraints. The decision variables are binary as mentioned above:

∀i ∈ [0 . . n], ∀ j ∈ [i + 1 . . n + 1], ∀k ∈ [1 . . m], γ(i, j, k), x i, j,k ∈ {0, 1} (12)

Objective

As previously mentioned, the objective function proposed in [START_REF] Bolat | Models and a genetic algorithm for static aircraft-gate assignment problem[END_REF] consists in minimizing the sum of the squares of the idle times, in order to find robust solutions towards operational arrival and departure delays. Therefore, the following minimization criterion is proposed: min

i∈[0..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k) I 2 i, j x i, j,k (13)
As exactly one binary variable x 0, j,k can be equal to 1, ∀k ∈ [1 . . m], i.e. the corresponding flight f j is the first on gate g k , and exactly one x i, j,k can be equal to

1 ∀i ∈ [1 . . n],
because of the constraint of Equation ( 8) (i.e. flight f i has a unique successor f j ), the cost defined by objective of Equation ( 13) corresponds to the transition cost of Definition 1:

cost tr = k∈[1..m] j∈[1..n+1] γ(0, j,k) c tr 0, j x 0, j,k + i∈[1..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k) c tr i, j x i, j,k (14)
Note that with this definition, the gates are not required to be uniform w.r.t. to the transition cost which could depend on k:

cost tr = i∈[0..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k)
c tr i, j,k x i, j,k

Multi-commodity Flow Model

As mentioned in the previous section, the GAP is similar to a Minimum-Cost Flow Problem [START_REF] Wang | Multicommodity network flows: A survey, part I: Applications and formulations[END_REF], more precisely an Unsplittable Multi-commodity Flow Problem (UMFP) [START_REF] Alvelos | Comparing branch-and-price algorithms for the unsplittable multicommodity flow problem[END_REF] with additional constraints, which translates into a tighter and smaller ILP model. We first introduce the corresponding graph flow model and then describe its translation to ILP as an improvement of the basic model constraints. Eventually, we refine this new GAP flow model to take into account other operational costs and constraints: we propose to combine the robustness with the taxiing cost within the objective and to prevent mutually exclusive gates to be used simultaneously (cf. Section 2.1).

Graph Model

The GAP can be seen as a kind of Unsplittable Multi-commodity Flow Problem (UMFP) [START_REF] Alvelos | Comparing branch-and-price algorithms for the unsplittable multicommodity flow problem[END_REF] with unit capacity for each arc and unit demand for each commodity. But our problem is defined on a directed multigraph G = (V, E) with parallel arcs between nodes v i and v j (corresponding to flights f i and f j ) whenever f e i f s j ∧ G i ∩ G j > 1 that restrict the compatible commodities rather than their amount. Arcs are labeled with the corresponding compatible gate index k (represented by a distinct color in the example of Figure 3) s.t. γ(i, j, k) are specified by triplets of V × V × [1 . . m]. More classically, the opening of gates are the sources of G and their closing are the sinks.

Definition 5 (Compatibility (Multi)Graph). The compatibility (multi)graph of a GAP instance G = (V, E, w, d) with weight w : E → R 0 and supply d : V → R is defined by:

• V = V F ∪ V s G ∪ V e G : (15a) V F = {v i , ∀i ∈ [1 . . n]} (15b) V s G = {v k 0 , ∀k ∈ [1 . . m]} (15c) V e G = {v k n+1 , ∀k ∈ [1 . . m]}
V F corresponds to regular flights, V s G to fictive openings of the gates and V e G to their endings, and |V| = 2m + n.

• E = E F ∪ E s G ∪ E e G ∪ E G : E F = {(v i , v j , k), ∀i ∈ [1 . . n -1], ∀ j ∈ [i + 1 . . n], ∀k ∈ [1 . . m] s.t. γ(i, j, k)} (16a) (16b) E s G = {(v k 0 , v j , k), ∀ j ∈ [1 . . n], ∀k ∈ [1 . . m], s.t. g k ∈ G j } (16c) E e G = {(v i , v k n+1 , k), ∀i ∈ [1 . . n], ∀k ∈ [1 . . m], s.t. g k ∈ G i } (16d) E G = {(v k 0 , v k n+1 , k), ∀k ∈ [1 . . m]}
E F corresponds to idle times between two regular successive flights, E s G to idle times between the opening of the gates and their first flights, E e G to the ones between their last flights and the closing of the gates, and E G to empty gates. |E| is equal to the number of compatible pairs of (real) flights plus 2 k∈[1..m] |F k |, the arcs from the opening and closing of the gates to compatible flights, i.e. in O(2mn + n 2 ) in the worst case.

• Each arc has a unit capacity 1 and a weight (or cost) w defined by:

w((v i , v j , k)) = c tr i, j = I 2 i, j
(also valid for start and end vertices v k 0 and v k n+1 ). Note that our model is generic enough to take into account non uniform resources (gates) for which the cost may depend on k as shown in Section 4.3 where taxiing times are added to the idle times.

• Each node has a supply value d defined by: . 4], to actual flights. A coloured arc joins two nodes if they can succeed each other on the corresponding gate, therefore multiple arcs may join the same pair of nodes (as between nodes v 1 and v 2 ). They are labelled with their weight (i.e. the cost of the associated idle time) indicated in hundreds of min 2 . The flow corresponding to the optimal solution is shown with bold arcs.

d(v) =            1 if v ∈ V s G (openings) -1 if v ∈ V e G (closings) 0 otherwise (i.e. v ∈ V F , regular flights) v 1 0 v 2 0 v 3 0 v 1 v 3 v 2 v 4 v 1
As an example, Figure 3 represents the graph model of the following GAP instance:

Example 1 (GAP Instance with 4 Flights and 3 Gates).

F = { f 1 , f 2 , f 3 , f 4 }, G = {g 1 , g 2 , g 3 }
with g s = 06:00, g e = 21:00 (in hh:mm format) and:

i f s i f e i G i
1 06:00 08:00 {g 1 , g 2 } 2 10:30 12:00 {g 1 , g 2 } 3 11:20 14:00 {g 2 , g 3 } 4 18:00 20:00 {g 1 , g 3 } Feasibility. Like a standard flow problem, an instance of the GAP is feasible iff there exists a binary flow φ : E → {0, 1} such that the imbalance between outgoing and incoming flows at each node is equal to the supply:

∀v ∈ V, e∈E + v φ(e) - e∈E - v φ(e) = d(v) (17)
g 1 g 2 g 3 f 1 f 4 f 2 f 3
6 : 0 0 7 : 0 0 8 : 0 0 9 : 0 0 1 0 : 0 0 1 1 : 0 0 1 2 : 0 0 1 3 : 0 0 1 4 : 0 0 1 5 : 0 0 1 6 : 0 0 1 7 : 0 0 1 8 : 0 0 1 9 : 0 0 2 0 : 0 0 2 1 : 0 0 with E + v = {e, ∀e = (u, v, k) ∈ E} the incoming arcs of node v and E - v = {e, ∀e = (v, w, k) ∈ E} its outgoing ones. However, for the GAP, the flow must also satisfy the following additional constraints.

Additional constraints. In general, a solution to a UMFP is not a solution to the GAP, because all internal nodes (corresponding to real flights) must be covered by the flow (Equation ( 18)) and paths must be "monochromatic" (Equation ( 19)):

• A solution to the GAP must cover all flights of F , therefore all nodes of V F must have a unit inflow and outflow:

∀v ∈ V F , e∈E + v φ(e) = e∈E - v φ(e) = 1 (18)
• The same gate must be used to enter and leave a node, therefore the flow imbalance Equation ( 17) must be expressed individually for each gate (or "color"):

∀v i ∈ V F , ∀g k ∈ G k , (v i ,v i ,k)∈E + v φ((v i , v i , k)) = (v i ,v i ,k)∈E - v φ((v i , v i , k)) (19)
The flow of each edge

(v i , v j , k) (or (v k 0 , v j , k), (v i , v k n+1 , k), (v k 0 , v k n+1 , k)
) corresponds to the binary decision variable x i, j,k (or x 0, j,k , x i,n+1,k , x 0,n+1,k ) of the basic ILP model described in Section 3.1, thus the assignment of flights to gate can be deduced accordingly.

Objective. The cost of a flow is defined by: e∈E w(e)φ(e) [START_REF] Kim | Airport Control Through Intelligent Gate Assignment[END_REF] As mentioned in the previous paragraph, this cost corresponds to the one of the basic ILP model described in Section 3.3, and is therefore also equal to the transition cost of Definition 1. Figure 4 depicts the Gantt diagram of an optimal solution to Example 1 with cost: cost tr = 10069 × 100min 2 .

ILP Model

The previously described GAP model as a UMFP with additional constraints can be translated in an ILP model equivalent to the basic model of Section 3, based on the same x i, j,k binary decision variables and cost expression but with tighter and fewer constraints. The UFMP flow imbalance constraints of Equation [START_REF] Gupta | Efficient algorithms for interval graphs and circular-arc graphs[END_REF] ensures that all sources (or gate openings) have a unit outflow as their supply value is 1:

∀k ∈ [1 . . m],
j∈ [1..n+1] γ(0, j,k)

x 0, j,k = 1 [START_REF] Li | Optimized assignment of airport gate configuration based on immune genetic algorithm[END_REF] and that all sinks (gate closings) have a unit inflow (supply value of -1):

∀k ∈ [1 . . m], i∈[0..n] γ(i,n+1,k) x i,n+1,k = 1 ( 22
)
which amounts to O(m) constraints.

The covering constraints of Equation ( 18) for nodes of V F (i.e. regular flights) lead to Equation ( 8) of the basic model, i.e. all internal nodes have a unit outflow:

∀i ∈ [1 . . n], j∈[i+1..n+1] k∈[1..m] γ(i, j,k) x i, j,k = 1 which amounts to O(n) constraints.
More specifically, the constraints of Equation ( 19) ensuring that the path of a commodity is monochromatic, can be directly translated to the flow balance constraint:

∀i ∈ [1 . . n], ∀k ∈ [1 . . m], i ∈[0..i-1] γ(i ,i,k) x i ,i,k = i ∈[i+1..n+1] γ(i,i ,k) x i,i ,k (23)
which amounts to O(mn) constraints.

Other constraints of the basic model that ensure that nodes of V F have a unit inflow (cf. Equation ( 7)), the continuity of path (cf. Equation ( 10)) or the number of variables equal to 1 (cf. Equation ( 9)) would be redundant in our flow model and are therefore omitted. Overall, this model has only O(nm) constraints in the worst case instead of O(mn 2 ) for the basic one.

Cost with Combined Taxiing

Even if the robustness of the schedule is our main operational concern in this study, the fuel consumption during taxiing has become of paramount importance in the last decade. Generally, the taxiing time cost of a flight does not only depend on the assigned gate, but also on the flight itself. As mentioned in Section 5.1, big international airport (like Paris-CDG) usually offer several runways for departure and arrival, thus the taxiing times of two flights assigned on the same gate could differ from each other, as their distance to a given gate is different if they do not use the same runway.

To attempt to minimize the fuel consumption induced by the amount of taxiing times, we introduce a new parameter t i,k that represents the taxiing time cost of aircraft f i if assigned to gate g k . Note that in the results presented in Section 5.5, t i,k is a taxiing time that only depends on the runway used by f i and the gate, but it could also depend on the type of aircraft to represent fuel consumption or CO 2 emission. If flight f i is assigned to g k , there is exactly one x i, j,k variable equal to 1, ∀ j ∈ [i + 1 . . n + 1], and all other x i, j,k with k = k are equal to 0, so the taxiing cost c t i of f i can be expressed as:

∀i ∈ [1 . . n], c t i = j∈[i+1..n+1] k∈[1..m] γ(i, j,k) t i,k x i, j,k (24)
and the global taxiing cost for all flights as:

cost t = i∈[1..n] c t i = i∈[1..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k) t i,k x i, j,k (25)
which can be integrated into the transition (or robustness) cost summation as is apparent from Equation ( 14). Obviously, no taxiing cost corresponds to the openings of the m gates taken into account in the first term of the sum of Equation ( 14), but we can obtain a uniform expression similar to Equation ( 13) by setting t 0,k = 0, ∀k ∈ [1 . . m]:

cost t = i∈[0..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k) t i,k x i, j,k (26)
To combine both costs in the objective, we choose a linear scalarization with a single weight parameter α ∈ [0, 1], and normalize them because the range of robustness (min 2 of gates idle times) is a priori much larger than the one of taxiing times (expressed in min in our study, but which can be chosen as any convenient operational cost, e.g. fuel consumption). To normalize the costs, we first compute their respective lower and upper bounds as follows.

Idle Time Upper Bound. The GAP can be relaxed as a Bin Packing (BP) problem where gates serve as bins of capacity g eg s and flights f i as items of volume f e if s i . The classical items configuration (stacked and divided as needed to completely fill the first bins) that corresponds to the BP lower bound on the number of bins, can be used to compute a simple upper bound for the sum of any idle time convex cost, in particular their square:

ub tr = m δ(g) 2 + (δ(g) -r) 2 (27)
where function δ : F ∪ G → R 0 returns the duration of its parameter, i.e. g eg s for any gate and f e if s i for flight f i , and is extended to take a set parameter and returns the sum of the duration of its elements, e.g.

δ(F ) = f i ∈F δ( f i ), and m = m -δ(G)-δ(F ) δ(g)
the number of empty gates and r = (δ(G)δ(F )) mod δ(g) the time taken by the flights scheduled on the last non-empty gate.

Idle Time Lower Bound. Still relaxing the GAP as a kind of BP problem, the lower bound for the overall sum of the squares of idle times is reached when all idle times have the same size and are evenly distributed along a continuous virtual gate made of all the gates stacked one after another. As mentioned in Section 3, there are exactly n + m idle time periods (possibly of null duration), therefore we can compute the following global lower bound on the idle time cost:

lb tr = (n + m) δ(G) -δ(F ) n + m 2 ( 28
)
Taxiing Time Bounds. As mentioned in Section 5.1, each departure and arrival at Paris-CDG airport is usually assigned to a fixed terminal and a fixed runway. Hence, the taxiing time bounds for a flight f i are the minimal and maximal times taken between any of its compatible gates G i and both its fixed arrival and departure runways. Therefore its upper bound can be obtained as:

ub taxi = i∈[1..n] max g k ∈G i (tt(ar i , g k ) + tt(g k , dr i )) ( 29
)
where ar i and dr i are the arrival and departure runways of flight f i , and tt is a function that returns the taxiing time between the location of its two parameters. The lower bound lb taxi is trivially computed by substituting the min function to the max one in Equation ( 29).

The combined cost that takes into account both robustness and taxiing time, cost α , can then be defined with the normalized terms (e.g. c ub-lb for cost c and its bounds) or equivalently (i.e. with the same optimal solution) by:

∀α ∈ [0, 1], cost α = i∈[0..n] j∈[i+1..n+1] k∈[1..m] γ(i, j,k) αρ tr t i,k + (1 -α)ρ taxi I 2 i, j x i, j,k (30)
with ρ tr = ub trlb tr and ρ taxi = ub taxilb taxi the respective ranges of both costs. The coefficient of the linear scalarization α can then be adjusted to influence the respective importance of the two criteria. Note that this new definition of the objective does not change the rest of our model.

Mutually Exclusive Gates

Depending on the layout of the terminals, some subsets of the gates cannot simultaneously be assigned to more than one flight, as they would otherwise block each other. Therefore, a set H of r subsets of G representing groups of mutually exclusive gates is additionally specified in the data:

H = {H 1 , . . . , H r }, with H h ⊆ G, ∀h ∈ [1 . . r].
We could consider such a subset as a single gate, but gates are not equivalent w.r.t. the taxiing time and their contribution to the cost must be individually accounted for.

To add this constraint in our flow model, we define O(mn) auxiliary variables y i,k that represent the allocation of flight f i to gate

g k ∈ G i : ∀i ∈ [1 . . n], ∀g k ∈ G i , y i,k = j∈[i+1..n+1] γ(i, j,k) x i, j,k (31)
Then for each mutually exclusive gates subset H h , we constrain that all overlapping flights pairs cannot be assigned to two of its gates simultaneously:

∀h ∈ [1 . . r], ∀g k , g l ∈ H h , k < l, ∀ f i ∈ F k , ∀ f j ∈ F l , [ f s i , f e i [ ∩ [ f s j , f e j [ = ∅, y i,k + y j,l 1 (32)
Groups of mutually exclusive gates H h usually comprise only a small number of gates (w.r.t. m), therefore the worst-case O(n 2 h∈[1.

.r] |H h | 2 ) additional constraints does not burden much the model.

Results

The models described in the previous sections have been implemented and run over real traffic data at Paris-CDG. All experiments were carried out on a standard workstation consisting of a 2.0 GHz Intel ® Xeon ® 32-core processor with 48 GB of memory running Debian GNU/Linux 9.6. We used the Gurobi Commercial Optimizer 8.1 [START_REF] Gurobi Optimization | Gurobi optimizer reference manual[END_REF] for the Integer Linear Programming models. Computation times to solve a 1-day instance using flow models for densest terminals are less than 100 s, and about 1 s is needed to optimally solve most other instances. A detailed report is provided in Section 5.3.

In this section, we first describe how we have extracted the traffic demand and the available gates from actual traffic records at Paris-CDG airport in Section 5.1 to generate a set of Gate Allocation Problem instances. Then, in Section 5.3 we compare the computation times of the different methods used to solve these instances. The implementation of DMAN can be very sensitive to the robustness of the allocation, as the slightest delay in a fragile allocation might quickly propagate and disrupt the departure sequence. Moreover, whereas ATC is concerned mainly with the feasibility of the assignment, our approach also focuses on minimizing the taxiing cost, which is directly related to fuel consumption, while preserving enough robustness in the schedule. Hence, Sections 5.4 and 5.5 discuss the benefits of a robust gate allocation in terms of idle time distribution and taxiing time differences between initial and robust traffic.

The results presented in this section refer to the actual allocation implemented by ATC on the day of traffic as "Initial", to the models presented in Section 3 as "JORS", and to the models presented in Section 4 as "Flow". Each model has two variants: R when only robustness is being maximized, and R+T when robustness is being maximized and taxiing time is being minimized.

Data

The data sample used for this study is available at: http://recherche.enac.fr/ ~wangrx/gap It consists in:

• actual traffic demands during a whole heavy month (July, 2017) at Paris-CDG airport; • a set of allowed aircraft types for each gate;

• an estimation of the average taxiing times needed between each runway and each gate (from each runway to each gate for arrivals and from each gate to each runway for departures).

The traffic demand is extracted from the data sample which provides for each flight: the arrival time at the gate, the departure time from the gate, the aircraft type, the runway and the gate used. The corresponding allocation serves as a reference solution named "Initial" in the following sections. Table 1 gives the number of flights (in average per day) and the number of gates by terminal in this data sample.

In some cases (probably due to some gate-to-gate movements that were not recorded), one of the two times was missing. The data were completed as follows: when the arrival (resp. departure) time is missing, the flight is considered at the gate from (resp. until) the beginning (resp. end) of the day if there is no other flight occupying the same gate before (resp. after) it; otherwise, the flight is considered at the gate 30 min before (resp. after) its departure (resp. arrival) time, which did not cause any gate conflict in our data set.

The sets of allowed aircraft types at each gate is also deduced from the actual traffic sample: we consider that the only aircraft types allowed at a given gate are the ones that actually used this gate during the month.

The estimation of taxiing times between gates and runways are computed using the standard routing rules described in the operations manual of the airport, after considering the average taxiing speed on each taxiway portion, as measured on actual radar records.

Bi-objective Parameterization

For the R+T strategies, we have tested different values for α ranging from 0 to 1, in objective function Equation [START_REF] Zhang | Optimization for gate re-assignment[END_REF], with 0 corresponding to optimizing taxiing time only and 1 to optimizing robustness only. Figure 5 shows how each individual cost is impacted by the value of the α parameter, for the 7th July 2017 instance (other instances exhibit a very similar pattern).

In Figure 5a, the value of both robustness cost and total taxiing time is plotted with respect to α. We notice that robustness cost is rapidly increasing for values of α below 0.5. On the other hand, with α over 0.9 (or even 0.8 for many instances), there is not much gain in the taxiing time. Thus, to achieve an interesting compromise between those two criteria, we recommend values of α in the interval [0.5, 0.8]. In the following sections, results are presented for α = 0.7.

Figure 5b provides a Pareto visualization of the same data. In an operational context, this parameter could be adjusted by airport managers to meet their own performance objectives.

Computation Times and Quality of the Solutions

We report in this section the performances of the various models with the criteria described in Sections 3 and 4 to solve the instances at Paris-CDG airport described in Section 5.1. We first show that an optimal solution to small instances of the problem (i.e. at the least busy terminals) can be easily reached by all models. We then discuss the weaknesses of JORS-R and JORS-R+T models on larger instances compared to their Flow-R and Flow-R+T counterparts.

Figure 6 displays the time to optimally solve a 1-day instance (averaged over the entire month of July, 2017) for each terminal and each implemented model. The four models are distinguished by their color in the figure. For the sake of readability, the figure has been separated in two parts, with Figure 6a focusing on small instances (with a density lower than 3 flights per gate) and Figure 6b showing results for larger instances.

For small instances (i.e. terminals A, B, C, I, L, M, R, and TVWXYZ2 ), Figure 6a shows that all models can obtain an optimal solution, and the optimality is always proved in very limited time (less than 1 s with Flow-R and Flow-R+T models). Flow models systematically outperform their JORS counterparts in terms of computation times, sometimes by two orders of magnitude (e.g. on terminals A and I). Interestingly, taking taxiing-time into account with the Flow-R+T model does not add to the computation time compared to the robustness-only Flow-R model.

Large instances (i.e. terminal D, E, F, J, K and Q) are significantly more challenging, as shown in Figure 6b. Considering only robustness as an optimization criterion, both Flow-R and JORS-R models still find solutions, again with an advantage of one order of magnitude in favor of the Flow-R model. For the largest and densest terminal F, the Flow-R model reaches optimality in less than 1 min, where JORS-R needs more than 20 min to obtain the same result. Additionally taking taxiing time into account with the Flow-R+T model does not incur additional costs compared to Flow-R: even the largest instances are optimally solved in less than 1 min. The JORS-R+T model, however, could not optimally solve any single instance within a reasonable time (the computations were stopped after 3 h). Even for the medium-sized terminal K, only 25 % of instances could be optimally solved. For instances at the largest and busiest terminal F, finding a first solution with JORS-R+T model took at least 30 min. If we impose a 5 min time limit for computation, all instances are optimally solved with the Flow-R+T model. As explained before, the JORS-R+T model could not reach (and even less prove) optimality; however, in many instances, it could still find good solutions within the time limit. Table 2 shows the discrepancy to the optimal cost by terminal for the JORS-R+T model, obtained after 5 min of computation. There is no figure for Terminal F as no solution could be found within the time limit. For other terminals, the solutions found are close to optimal, with approximately 1 % to 5 % additional cost.

As the Flow-R and Flow-R+T models are strictly and significantly faster than JORS-R and JORS-R+T respectively, the following sections only present results for the first two, as well as comparisons with the initial gate allocation extracted from the data (see Section 5.1). Also, these computation times are suitable for an allocation process performed the day before operations, or even on the day of operations if some event breaks the solution and a new allocation needs to be devised. In such a case, an additional criterion could be introduced that would help find a new solution as close as possible to the previous one, so that the gate allocation would remain mostly unchanged. This is, however, not covered in the current study.

Robustness of the Gate Allocation

To assess the robustness of a gate allocation, we measure the idle times occurring in the corresponding schedule, as too small idle times are likely to cause issues on the day of operations. Table 3 shows the average idle time: the average for the Initial allocation is 2 h, when the average for both Flow-R and Flow-R+T allocations are 2 h 40 min.

More interestingly, Figure 7 shows the distribution of idle times for the entire airport. Those distributions are plotted for the Initial allocation (in blue) and for both Flow-R (in orange) and Flow-R+T (in green) models. As we can see, the Initial allocation leads to a high number of short idle times: more than 1000 idle times between 0 and 10 min and With the Flow-R+T model, the solutions are still significantly more robust than the Initial allocation, with again a high reduction of short idle times. Overall, such idle times can absorb a large part of the uncertainty that can be observed on flight arrival and departure times, with no re-allocation needed in most cases. The presence of long idle times in all distributions (more than 2 h) comes from the fact that all terminals are aggregated in this figure, even the smallest ones with a very low density.

Table 4 and Figure 8 present similar results, focusing on terminal F only. Even for this particularly busy terminal, both Flow-R and Flow-R+T models manage to provide a solution that significantly reduces occurrences of smaller idle times. Also, there are very few idle times of more than 2 h, as the density (flights-to-doors ratio) is very high at this terminal.

Impact on Taxiing Time

Another main objective of the GAP is the minimization of taxiing times. As mentioned in Section 2, Paris-CDG large terminals offer several gates with various ways to enter or leave the terminal, so that taxiing times can be very different from one gate to another in the same terminal.

Table 5 gives the total amount of saved taxiing time during the month of July 2017, comparing the Initial (actual) gate allocation with the two gate allocation strategies (Flow-R and Flow-R+T). With the Flow-R strategy, the total taxiing time is almost unchanged, which shows that the gain in robustness commented in Section 5.4 does not affect taxiing times. With the Flow-R+T strategy, a slight reduction by 2.67 % of the taxiing time is observed, which represents 180.97 h of saved taxiing time during the month, i.e. 6 h per day in average. As discussed in Section 5.4, this gain in taxiing time does not significantly penalize the robustness of the gate allocation.

Results are also detailed by terminal to point out the ones with the greatest reduction in taxiing time. Terminal R appears to be the most influenced one with a decrease of 14.54 % in taxiing time, but this result is not significant as the number of flights in this terminal is very low. More important terminals, from A to F, reach a decrease between 0.63 % to 4.08 %, which represents 77 h of saved taxiing time during July 2017, contributing to 42 % of the total saved taxiing time.

Merging Terminals

As described in Section 5.3, Flow-R and Flow-R+T models can solve instances of all the terminals in less than 50 s. In this part, we test the performances of our best models on bigger instances by merging the terminals at Paris-CDG 1 and 2 (see Figure 2), where passengers can walk from one terminal to another. These instances also provide an estimation of the gain that can be expected in terms of robustness and taxiing time, were the terminal of each flight not fixed -which would be outside of the scope of current operational constraints at Paris-CDG airport. Three different sets of merged terminals were considered:

• TVWXYZ at Paris-CDG 1;

• ABCD or ABCDEF at Paris-CDG 2. 6 gives the mean computation time needed to find and prove an optimal solution with models Flow-R and Flow-R+T. The resolution time ranges from a few seconds to 45 min, which remains reasonable for a gate allocation process that is usually performed the day before operations. If the allocation process had to be run again to accommodate delays on the day of operations, the instances would be much smaller (by not considering the entire day of traffic), so that computing times would probably be small enough. The most time-consuming instance comes from the merging of all Paris-CDG 2 terminals (i.e. A to F). 7), we can observe a new significant improvement in the robustness of the gate allocation with the two tested criteria, as the number of short idle times (i.e. less than one hour) is further decreased. Table 7 gives the variation of the average taxiing time for the three sets of merged terminals. We can notice the significant taxiing time reduction that can be expected using the Robust+Taxiing criterion at Paris-CDG 2: 13.13 % by merging terminals A, B, C and D, and 14.86 % by adding terminals E and F. Globally, with merged terminals, the R+T criterion applied on the whole airport allows to save 684 h of taxiing time during the month of July 2017, which corresponds to almost 23 h per day. These results demonstrate how beneficial the relaxing of the terminal allocation of some flights could be at Paris-CDG airport, in terms of gate allocation robustness and saved taxiing time.

Conclusion

In this study, we present and test a new efficient multi-commodity flow model for the gate allocation problem at big airports. To the best of our knowledge, this is the first multi-commodity flow model that can handle complex constraints like mutually exclusive gates, while optimizing both the robustness and the taxiing times, which are obviously two main criteria for large airports. The efficiency of the model is tested on actual instances, extracted from traffic records at Paris-CDG airport, during one of its heaviest month (July, 2017). All the instances was successfully solved with the ILP solver.

We compare two versions of the min-cost flow model (Flow-R and Flow-R+T) to the two corresponding versions of a more classical model (JORS-R and JORS-R+T). Instances of small density are optimally solved in less than one minute by all of the four models. For larger instances, we defined a 300 s limit for the time computation, in order to make the approach compatible with a real-time context (e.g. gate re-allocation in case of operational disruption). With this limit, the Flow-R and Flow-R+T models still provide optimal solutions for all the terminals in less than 100 s, while JORS-R and JORS-R+T models not always reach optimality, and sometimes do not even find a solution.

Results show the overall improvement in robustness and taxiing time after gate allocation: the idle times distributions are better balanced, with a significant reduction of the number of small idle times, and the taxiing times are reduced by a total of 180 h during the whole month, which represents 2.67 % of the total taxiing time.

Flow-R and Flow-R+T models can also be applied on larger instances after merging all the terminals of Paris-CDG 2 (i.e. A, B, C, D, E, F) and of Paris-CDG 1 (i.e. T, V, W, X, Y, Z). The results show that the robustness of the gate allocation is still improved while the taxiing time is decreased by 684 h (i.e. 10.1 % of the total taxiing time), moreover our min-cost flow model still finds an optimal solution in a reasonable time for these larger instances.

In further works, we project to test the gate allocation into some more realistic simulations of the airport traffic, considering the conflicts that can appear between aircraft entering or exiting the gates areas. This will provide a more accurate estimation of the gain that can be expected in terms of gates conflicts and taxiing times diminution with our gate allocation model.

Figure 1 :

 1 Figure 1: Paris-CDG 1: terminals T, U, V, W, X, Y, Z.

Figure 2 :

 2 Figure 2: Paris-CDG 2: terminals A, B, C, D, E, F.

  i,n+1 between the last flight of each gate f i and its closing, represented by fictive flight f n+1 : Definition 2 (Transition Cost). The transition cost of a solution to an instance of the GAP is defined by: where fst : [1 . . m] → [1 . . n + 1] is a function such that fst(k) returns the index (in [1 . . n]) of the first flight assigned to gate g k or n + 1 if the gate is unoccupied, and succ : [1 . . n] → [1 . . n + 1] is a function such that succ(i) returns the index (in [1 . . n], or more precisely in [i + 1 .

Figure 3 :

 3 Figure 3: Coloured representation of the compatibility graph of Example 1 where "commodities", i.e. the three gates, are represented by different colors. Fictive nodes v k 0 and v k n+1 are associated respectively to the opening and ending of gate k, and nodes v i , ∀i ∈ [1 .. 4], to actual flights. A coloured arc joins two nodes if they can succeed each other on the corresponding gate, therefore multiple arcs may join the same pair of nodes (as between nodes v 1 and v 2 ). They are labelled with their weight (i.e. the cost of the associated idle time) indicated in hundreds of min2 . The flow corresponding to the optimal solution is shown with bold arcs.

Figure 4 :

 4 Figure 4: Gantt diagram of the optimal solution to Example 1.

  Visualization of the Pareto front. The value of α is indicated next to each dot.

Figure 5 :

 5 Figure 5: Influence of α on bi-objective optimization (7th July, 2017 at Paris-CDG airport).

Figure 6 :

 6 Figure 6: Comparison of computation times to find an optimal solution with JORS and Flow model.

Figure 7 :

 7 Figure 7: Distribution of idle times for all terminals over all tested instances.

Figure 8 :

 8 Figure 8: Distribution of idle times at terminal F over all tested instances.

Figure 9 :

 9 Figure 9: Distribution of idle times for all terminals with ABCDEF and TVWXYZ combinations.

Figure 9

 9 Figure 9 shows the distribution of idle times resulting from the merging of terminals TVWXYZ at Paris-CDG 1 and ABCDEF at Paris-CDG 2. Compared to the distributions

Table 1 :

 1 Number of flights and gates by terminal (per day in average). The density is the ratio of the average number of flights by the number of gates.

	Terminal	Flights Gates Density
	A B C D E F I J K L M Q R TVWXYZ	47 15 14 84 111 185 46 72 53 18 13 50 12 96	17 10 21 18 23 27 51 20 19 7 10 17 12 43	2.8 1.5 0.7 4.6 4.8 6.8 0.9 3.6 2.8 2.6 1.3 3.0 1.0 2.2
	Total	816	295	2.8

Table 2 :

 2 Average gap between the cost of optimal solution and the best solution found with JORS-R+T model within a 5 min time limit.

	Terminal	K	Q	J	D	E	F
	Best solution +1.08 % +1.11 % +4.07 % +3.11 % +4.63 % -

Table 3 :

 3 Average idle time between two consecutive aircraft at the same gate.

	Allocation	Initial Flow-R Flow-R+T
	Average idle time 120 min 157 min	157 min

Table 4 :

 4 Average idle time between two consecutive aircraft at the same gate for terminal F only.

	Allocation	Initial Flow-R Flow-R+T
	Average idle time 65 min	73 min	72 min
	more than 1300 idle times between 10 and 20 min. The Flow-R allocation drastically reduces the number of idle times of less than an hour, thus increasing the robustness of the allocation as expected. With this model, only 150 idle times are less than 10 min.

Table 6 :

 6 Computation time for merged terminals instances with the Flow model.

	Terminals	Flights Gates Time R Time R+T
	ABCD ABCDEF TVWXYZ	167 461 95	59 109 409.6 s 22.7 s 42 2.6 s	6.8 s 2631.9 s 2.94 s
	Table

Table 7 :

 7 Taxiing time comparison between the different gate allocation strategies.

	Terminal	Flights	Taxiing time (and difference w.r.t. Initial) Initial Robust Robust + Taxiing
	ABCD ABCDEF TVWXYZ others	8,615 23,802 5,237 13,243 10 min 06 s 10 min 03 s (-0.49 %) 9 min 40 s (-4.3 %) 9 min 54 s 9 min 50 s (-0.67 %) 8 min 36 s (-13.13 %) 9 min 52 s 9 min 55 s (+0.68 %) 8 min 24 s (-14.86 %) 7 min 12 s 7 min 20 s (+0.18 %) 7 min 04 s (-1.85 %)
	All Total	42,282	9 min 37 s	9 min 39 s (+0.37 %) 8 min 38 s (-10.1 %) +0.33 h -684.04 h
	obtained in Section 5.4 (see Figure

This problem is sometimes called Airport Gate Allocation Problem (AGAP) or Stand Allocation Problem (SAP) in the literature.

Terminals T, V, W, X, Y, and Z must be solved simultaneously as there are mutually exclusive gates belonging to different terminals in this set. Dependent gates information is accessible at http://recherche. enac.fr/ ~wangrx/gap.

Funding

This article is supported by the funds of basic scientific research in central universities, Civil Aviation University of China. [No. 3122021084].


Jean-Luc Faillie
email: jl-faillie@chu-montpellier.fr

Samy Babai

Sabrina Crépin

Virginie Bres

Marie-Laure Laroche

Hervé Le Louet

Pierre Petit

Jean-Louis Montastruc

Dominique Hillaire-Buys

Dominique Hillaire

Pancreatitis associated with the use of GLP-1 analogs and DPP-4 inhibitors: a case/non-case study from the French Pharmacovigilance Database

Keywords: Type 2 diabetes mellitus, glucagon-like peptide 1 analogs, dipeptidyl peptidase 4 inhibitors, pancreatitis, pharmacovigilance

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Since 2005, new pharmacological classes of drugs potentiating the activity of the incretin hormones have been marketed in the USA in second or third-line treatment of type 2 diabetes mellitus: analogs of glucagon-like peptide 1 (GLP-1) and dipeptidyl peptidase-4 (DPP-4) inhibitors. Safety concerns regarding a potential risk for pancreatitis rapidly raised as several cases were reported. In October 2007, the FDA issued a first alert regarding 30 cases of pancreatitis associated with exenatide, the first GLP-1 analog [START_REF]) -8/2008 Update[END_REF]. Postmarketing cases reported in the FDA Adverse Event Reporting System (AERS) database were analyzed by [START_REF] Elashoff | Pancreatitis, pancreatic, and thyroid cancer with glucagon-like peptide-1-based therapies[END_REF] incretin-based drugs were associated with a 10-fold increase of reported pancreatitis [START_REF] Elashoff | Pancreatitis, pancreatic, and thyroid cancer with glucagon-like peptide-1-based therapies[END_REF]. This study has been a large subject of debate regarding the methodology used and the author's conclusions [START_REF] Hawkes | Journal withdraws article after complaints from drug manufacturers[END_REF][START_REF] Hauben | Data mining in pharmacovigilance: the need for a balanced perspective[END_REF]. Another case/non-case study of pancreatitis using the same database for the period 2004-2009, found disproportionality with exenatide (ROR: 1.76 [95% CI: 1.61-1.92]) and with sitagliptin, the first DPP-4 inhibitor (ROR: 1.86 [95% CI: 1.54-2.24]) [START_REF] Raschi | The association of pancreatitis with antidiabetic drug use: gaining insight through the FDA pharmacovigilance database[END_REF]. In France, the first incretin-based drugs were marketed in March 2008. The objective of this study is to investigate, five years after initial marketing, the association between exposure to incretin-based drugs (GLP-1 analogs and DPP-4 inhibitors) and reports of pancreatitis in the French Pharmacovigilance Database.

Methods

The French Pharmacovigilance Database registers all adverse drug reactions (ADR) spontaneously reported by health professionals to the 31 French Regional Pharmacovigilance Centers since 1985. ADRs are coded according to the MedDRA® (Medical Dictionary for Regulatory Activities) classification [START_REF] Hauben | The role of data mining in pharmacovigilance[END_REF]. "Serious" ADRs include reactions resulting in death, life-threatening event, hospitalization, prolongation of existing hospitalization, significant disability, congenital anomaly, birth defect or other significant medical event. We selected serious ADRs reported in type 2 diabetic patients between 2003 and 2013. Analyses were conducted on serious ADRs reported between March 2008 (first marketing of an incretin-based drug in France) and March 2013. Type 2 diabetes was defined by the presence in the report of an antihyperglycemic agent, alone or in association, except when insulin was reported alone.

The case/non-case method was chosen to assess the disproportionality of reports for pancreatitis associated with incretin-based drugs compared to other type 2 diabetic drugs [START_REF] Wilson | Application of data mining techniques in pharmacovigilance[END_REF]. Cases were defined as reports of pancreatitis according to the MedDRA terms corresponding for acute or chronic pancreatitis or elevated pancreatic enzymes considered serious by the reporter. Other reports of serious ADRs were considered non-cases. Exposure was defined by the presence in the report of one or more incretin-based drugs currently marketed in France: GLP-1 analog (exenatide, liraglutide) or DPP-4 inhibitor (sitagliptin, vildagliptin, saxagliptin), whether or not the drug was suspected of causing the reaction. Disproportionality was assessed by calculating reporting odds ratios (ROR) and corresponding 95% confidence intervals for each antihyperglycemic drug and for the pharmacological classes of incretin-based drugs. The ROR is used to compare the rate of exposure to incretin-based drugs in cases of pancreatitis versus that in non cases [START_REF] Bate | Quantitative signal detection using spontaneous ADR reporting[END_REF][START_REF] Moore | Reports of hypoglycaemia associated with the use of ACE inhibitors and other drugs: a case/non-case study in the French pharmacovigilance system database[END_REF]. Logistic regression model was used to adjust for potential confounders including age, gender, history of pancreatitis, other antihyperglycemic drugs (metformin, sulfonylurea derivatives, glinides, thiazolidinediones, acarbose, and insulin) and other drugs associated with pancreatitis (angiotensin converting enzyme inhibitors, angiotensin receptor antagonists, statins, tetracycline, metronidazole, isoniazid, pentavalent antimony, azathioprine, mercaptopurine, 5-aminosalicylic acid drugs, valproic acid, estrogens, codeine, and sulindac) [START_REF] Balani | Drug-induced pancreatitis : incidence, management and prevention[END_REF].

A sensitivity analysis was performed using a stricter definition of pancreatitis when elevated pancreatic enzymes considered serious were excluded. Another model studied disproportionality of incretin-based drugs differentiating whether they were alone or in association with other antihyperglycemic agents. In order to study variation in reporting rate and the presence of a notoriety effect [START_REF] Pariente | Impact of safety alerts on measures of disproportionality in spontaneous reporting databases: the notoriety bias[END_REF], we performed a year by year analysis and compared the RORs before and after March 2011 (when the first study suggesting a risk of pancreatitis in the FDA AERS database was published) [START_REF] Elashoff | Pancreatitis, pancreatic, and thyroid cancer with glucagon-like peptide-1-based therapies[END_REF].

Analyses were conducted using SAS software version 9.2 (SAS Institute, Cary, NC, USA).

Results

During the study period from March 2008 to March 2013, we identified 3,109 serious ADRs involving antihyperglycemic agents. In total, 147 were cases of pancreatitis and 2962 non cases. No alcoholic, lupic or autoimmune pancreatitis were observed. Among the cases, 122 (82.99%) involved incretin-based drugs: 55 (37.41%) GLP-1 analogs, 64 (43.54%) DPP-4 inhibitors and 3 (2.04%) both of them. Since their marketing, reports of pancreatitis involving an incretin-based drug increased every year for both pharmacological classes (figure 1.)

Characteristics of pancreatitis cases are shown in table 1. Mean age was 60.4 ± 11.4 for incretin-based drugs (59.0 ± 10.2 for GLP-1 agonists and 61.5 ± 12.2 for DPP-4 inhibitors) and 65.2 ± 13.8 for other antihyperglycemic drugs. Males represented 58% of pancreatitis cases in association with incretin-based drugs. GLP-1 agonists and DPP-4 inhibitors were the only antihyperglycemic drugs reported in, respectively, 46.6% and 25.4% of the cases. Overall, hospitalization rate was high (89.8%) and fatal outcomes occurred in 2.1% of the reports. BMI or alcohol use were frequently unavailable in the analysed reports. 2). Increased reporting rate was found for any individual drugs of these pharmacological classes. No other antihyperglycemic or concomitant drug showed significant disproportionality. In the model assessing all incretinbased drugs, older age was significantly associated with lower reporting rate: compared to less than 60, adjusted ROR for age 60 to 69, age 70 to 79 and age ≥80 were 0.70 [95%CI: 0.46-1.08], 0.30 [95%CI: 0.17-0.51] and 0.34 [95%CI: 0.18-0.66] respectively. Pancreatitis history and male gender were also associated with a higher reporting rate (adjusted ROR: 4.32 [95%CI: 1.75-10.69] and 1.44 [95%CI: 0.99-2.08] respectively). Sensitivity analysis using a stricter definition of pancreatitis excluded 12 cases of elevated pancreatic enzymes considered serious among the 147 pancreatitis cases. Results were similar to those with the broader definition (adjusted ROR for all incretin-based drugs: 14.86 [95%CI: 9.24-23.92]). The model differentiating whether incretinbased drugs were alone or in association with other antihyperglycemic drugs showed a stronger disproportionality for incretin-based drugs alone (adjusted ROR: 30.30 [95%CI: 13.75-66.80]) than in association (adjusted ROR: 13.56 [95%CI: 8.35-22.02]).

Temporal analysis of cumulative ROR showed that disproportionality was present since the first year of marketing (figure 2). During marketing in France, no evidence for notoriety effect was found according to March 2011 when the first study using the FDA AERS database suggested a risk of pancreatitis [START_REF] Elashoff | Pancreatitis, pancreatic, and thyroid cancer with glucagon-like peptide-1-based therapies[END_REF].

Discussion

Our study demonstrates that reported pancreatitis in France were more than 15-fold increased with incretinbased drugs compared to other antihyperglycemic agents. Significant disproportionality was only found for incretinbased drugs, since their first year of marketing and among these pharmacological classes, a stronger association was found for GLP-1 analogs. Considering that the rate of pancreatitis associated with incretin-based drugs which are effectively reported to pharmacovigilance system could range from 5% to 20% [START_REF] Hazell | Under-reporting of adverse drug reactions : a systematic review[END_REF], the 122 cases reported between 2008 and 2013 in France could hypothetically represent 610 to 2,440 actual cases.

The results presented here are consistent with previous disproportionality studies conducted by Elashoff et al. [START_REF] Elashoff | Pancreatitis, pancreatic, and thyroid cancer with glucagon-like peptide-1-based therapies[END_REF] and Raschi et al. [START_REF] Raschi | The association of pancreatitis with antidiabetic drug use: gaining insight through the FDA pharmacovigilance database[END_REF] using data from the FDA AERS. However, we found higher ROR estimates for the association between incretin-based drugs and reported pancreatitis. This difference may be explained by the fact that we only selected "serious" ADRs which occurred in diabetic patients. Another recent analysis of the FDA AERS monitored reports of "serious" adverse events from July 2011 to July 2012 and found that reports of pancreatitis were 28.5 and 20.8 times more likely with respectively GLP-1 analogs and DPP-4 inhibitors than with control antihyperglycemic drugs (sulfonylureas and metformin) [START_REF] Cohen | Reports of pancreatitis are 20-30 times more likely with GLP-1 drugs, analysis finds[END_REF][START_REF]Perspectives on GLP-1 agents for diabetes[END_REF]. These results are similar to those we found in the present study.

Case/non case analyses are subject to the inherent limitations of data-mining approaches from spontaneous reporting databases. Under-reporting is the main drawback of this type of data [START_REF] Van Der Heijden | On the assessment of adverse drug reactions from spontaneous reporting systems: the influence of under-reporting on odds ratios[END_REF]. Indeed, under-reporting may vary depending on the antihyperglycemic drug involved in the ADR and may be less important for new drugs such as incretin-based drugs. Furthermore, since relationship between pancreatitis and incretin-based drugs was highlighted in previous alerts by the US FDA (before marketing in France) [START_REF]) -8/2008 Update[END_REF] and was mentioned in the corresponding products information [START_REF] Galvus | Summaries of Product Characteristics[END_REF][START_REF] Onglyza | Summaries of Product Characteristics[END_REF][START_REF] Byetta | Summaries of Product Characteristics[END_REF][START_REF] Januvia | Summaries of Product Characteristics[END_REF][START_REF] Onglyza | Summaries of Product Characteristics[END_REF], clinicians could have been more aware of this particular ADR and could have more frequently detected and reported pancreatitis in their patients using incretin-based drugs. Hence, the notoriety effect inherent of studies using pharmacovigilance database cannot be excluded although it is difficult to assess in which extent this could have affected our results.

Incomplete data in spontaneous reports is also a limiting point since it may concern important information such as risk factors, medical or exposure history. Although we performed adjustment for age, gender, pancreatitis history and use of other drugs, confounding bias may still exist since we could not take into account risk factors such as smocking, BMI and alcohol consumption because of too numerous missing data for these variables in the reports. Precisely, adjustment for BMI should have been particularly relevant because obesity is often associated with a higher prescription of GLP-1 analogs which are known to reduce weight. Since obesity is a risk factor for pancreatitis, lack of adjustment for this confounder could have over-estimated corresponding ROR and might be an explanation for the higher ROR observed for GLP-1 analogs compared to DPP-4 inhibitors in this study. Lack of adjustment for alcohol consumption, a major risk factor for pancreatitis, is also an issue. However, since alcohol consumption is unlikely related to drug exposure, it is not believed to have played an important role in this analysis. Moreover, incompleteness of data on exposure history and time delays, which is frequent in pharmacovigilance databases, did not allow us to study the time period from first administration to onset of pancreatitis, or the impact of a newly increase in dosage of incretin-based drugs.

Several studies showed an increased risk of pancreatitis in patients with type 2 diabetes suggesting a role of diabetes itself [START_REF] Gonzalez-Perez | Acute pancreatitis in association with type 2 diabetes and antidiabetic drugs: a population-based cohort study[END_REF][START_REF] Noel | Increased risk of acute pancreatitis and biliary disease observed in patients with type 2 diabetes: a retrospective cohort study[END_REF][START_REF] Girman | Patients with type 2 diabetes mellitus have higher risk for acute pancreatitis compared with those without diabetes[END_REF]. In order to control confounding by this condition, we decided to perform this study only in reports which involved antihyperglycemic agents. However, our study could not take into account the potential effect played by the duration or the severity of diabetes. Thus, given the fact that incretin-based drugs are defined as second or third-line drugs, confounding by indication cannot be excluded. Nevertheless, considering that more severe diabetic patients are often receiving association of antihyperglycemic agents, the stronger disproportionality found for incretinbased drugs when they were reported alone does not support this potential bias.

Using antihyperglycemic drugs (excluding insulin alone) as a proxy for type 2 diabetes could have led to selection bias. For instance, it was not possible to exclude patients treated by metformin for polycystic ovary syndrome. Nonetheless, this indication is rare and probably negligeable. Likewise, exclusion of patients only treated by insulin alone in order to avoid selecting type 1 diabetes patients could have eliminated patients with advanced stage of type 2 diabetes. This could have led to under-estimate ROR for incretin-based drugs since type 2 diabetes patients treated by insulin are suspected to show a lower incidence of pancreatitis [START_REF] Gonzalez-Perez | Acute pancreatitis in association with type 2 diabetes and antidiabetic drugs: a population-based cohort study[END_REF]. However, this selection error might have a minimal impact on our results as it likely concerns a small sample of patients. Another important point regarding population selection is the fact that, in this analysis, we have chosen to focus only on serious ADRs. This approach allowed us to compare pancreatitis reports to other serious events reports and to avoid a masking effect by numerous and irrelevant non-serious ADRs.

It is known that chronic pancreatitis can give rise to diabetes (which can potentially be treated by antihyperglycemic drugs) and can cause several occurrences of acute pancreatitis. Thereby, errors in chronology of events may occur in pharmacovigilance reports, and reverse causation cannot be excluded in the association between antihyperglycemic drugs and pancreatitis. However, this phenomenon is unlileky to happen since all pancreatitis cases analyzed in the present study have been initially validated by physicians.

Despite the fact that the disproportionality approach cannot give quantitative risk assessment for a specific ADR (ROR cannot be interpreted as risk ratio), case/non case studies analyzing "real life" pharmacovigilance data have shown in the past to be useful to detect drug safety signals [START_REF] Montastruc | Benefits and strengths of the disproportionality analysis for identification of adverse drug reactions in a pharmacovigilance database[END_REF][START_REF] Sommet | Was the thrombotic risk of rofecoxib predictable from the French Pharmacovigilance Database before 30 September 2004?[END_REF]. As soon as 2009, animal studies have suggested a link between incretin-based drugs and pancreatic abnormalities including pancreatitis [START_REF] Matveyenko | Beneficial endocrine but adverse exocrine effects of sitagliptin in the human islet amyloid polypeptide transgenic rat model of type 2 diabetes: interactions with metformin[END_REF]. A recent review of the preclinical and postmaketing studies summarised by Butler et al. suggested a plausible pharmacological mechanism which involves pancreatic duct occlusion enhanced by duct proliferation secondary to activation of GLP-1 receptors in the pancreas [START_REF] Butler | A Critical Analysis of the Clinical Use of Incretin-Based Therapies: Are the GLP-1 therapies safe?[END_REF]. Postmortem study in organ donors exposed to incretin-based drugs also showed pancreatic enlargement and precancerous changes [START_REF] Butler | Marked Expansion of Exocrine and Endocrine Pancreas with Incretin Therapy in Humans with increased Exocrine Pancreas Dysplasia and the potential for Glucagon-producing Neuroendocrine Tumors[END_REF]. These studies have raised a large debate [START_REF] Nauck | A critical analysis of the clinical use of incretin-based therapies: The benefits by far outweigh the potential risks[END_REF][START_REF] Kahn | Incretin therapy and islet pathology: a time for caution[END_REF][START_REF] Drucker | Incretin action in the pancreas: Potential promise, possible perils, and pathological pitfalls[END_REF]. In addition, evidence brought by clinical trials were insufficent to support a risk of pancreatitits in meta-analyses of GLP-1 analogs [START_REF] Alves | A meta-analysis of serious adverse events reported with exenatide and liraglutide: acute pancreatitis and cancer[END_REF] or sitagliptin [START_REF] Engel | Safety and tolerability of sitagliptin in type 2 diabetes: pooled analysis of 25 clinical studies[END_REF]. Rencently, two large randomised controlled trials did not found any increased pancreatitis risk with saxagliptin and alogliptin [START_REF] Scirica | Saxagliptin and Cardiovascular Outcomes in Patients with Type 2 Diabetes Mellitus[END_REF][START_REF] White | Alogliptin after Acute Coronary Syndrome in Patients with Type 2 Diabetes[END_REF]. Pharmacoepidemiological studies using large medical database are contradictory. The first studies which focused on this ADR failed to demonstrated significant association [START_REF] Dore | A cohort study of acute pancreatitis in relation to exenatide use[END_REF][START_REF] Dore | Use of a claims-based active drug safety surveillance system to assess the risk of acute pancreatitis with exenatide or sitagliptin compared to metformin or glyburide[END_REF][START_REF] Garg | Acute pancreatitis in type 2 diabetes treated with exenatide or sitagliptin: a retrospective observational pharmacy claims analysis[END_REF][START_REF] Romley | Exenatide therapy and the risk of pancreatitis and pancreatic cancer in a privately insured population[END_REF][START_REF] Wenten | Relative risk of acute pancreatitis in initiators of exenatide twice daily compared with other anti-diabetic medication: a follow-up study[END_REF]. Recently, a case control study using a large US administrative database between February 2005 and December 2008 suggested that current use of GLP-1-based drugs was associated with increased risk of hospitalization for acute pancreatitis (OR 2.24 [95% CI, 1.36-3.68]) [START_REF] Singh | Glucagonlike peptide 1-based therapies and risk of hospitalization for acute pancreatitis in type 2 diabetes mellitus: a population-based matched case-control study[END_REF]. Soon after, a retrospective population based cohort study in a large US claims database did not confirm this relationship [START_REF] Eurich | Comparative safety and effectiveness of sitagliptin in patients with type 2 diabetes: retrospective population based cohort study[END_REF]. Nonetheless, theses observational studies are still subject to substantial methodological drawacks (confounding, information and selection bias, statistical power) which cannot allow to draw firm conclusions.

Conclusions

Despite some limitations, the present study showed that incretin-based drugs were associated with an almost 16-fold increase in reported pancreatitis in France. These results strengthen and extend the pharmacovigilance signal undercovered by post-marketing spontaneous reports. It is all the more urgent to investigate the potential risk of pancreatitis associated with the use of incretin-based drugs as their benefits in terms of cardiovascular morbi-mortality have not been clearly demonstrated.

Table 1 .

 1 Characteristics of pancreatitis cases associated with incretin-based drugs and other antihyperglycemic drugs (French Pharmacovigilance Database, March 2008-March 2013).According to univariate analysis, use of all incretin-based drugs, younger age and pancreatitis history were associated with a higher reporting rate of pancreatitis. Among other antihyperglycemic drugs, thiazolidinediones and sulfonylureas/glinides showed a lower reporting rate (crude RORs: 0.16 [95%CI: 0.04-0.65] and 0.65 [95%CI: 0.46-0.91] respectively), other drugs did not reached statistical significance. After adjustment for age, gender, history of pancreatitis, other antihyperglycemic agents and other drugs associated with pancreatitis, disproportionality was found for all incretin-based drugs (overall adjusted ROR: 15.62 [95%CI: 9.81-24.87]), with a stronger signal for GLP-1 analogs (adjusted ROR: 29.36 [95%CI: 16.02-53.81]) than for DPP-4 inhibitors (adjusted ROR: 12.08 [95%CI: 7.30-20.00]) (table

						Other
		GLP-1 agonists	DPP-4 inhibitors	antihyperglycemic
	Characteristics	(N=58)		(N=67)		drugs (N=25)
		no.	%	no.	%	no.	%
	Age
	<60	32	55.2	29	43.3	9	36.0
	60-69	19	32.8	20	29.9	6	24.0
	70-79	4	6.9	11	16.4	4	16.0
	≥80	2	3.4	5	7.5	6	24.0
	Missing	1	1.7	2	3.0	0	0.0
	BMI
	<20	0	0.0	1	1.5	0	0.0
	20-24	4	6.9	2	3.0	6	24.0
	25-29	6	10.3	11	16.4	3	12.0
	30-34	9	15.5	10	14.9	1	4.0
	≥35	14	24.1	4	6.0	0	0.0
	Missing	25	43.1	39	58.2	15	60.0
	Gender
	Females	24	41.4	28	41.8	13	52.0
	Males	34	58.6	38	56.7	12	48.0
	Missing	0	0.0	1	1.5	0	0.0
	History of pancreatitis	3	5.2	4	6.0	2	8.0
	Alcohol use
	Yes	1	1.7	3	4.5	4	16.0
	No	2	3.4	2	3.0	4	16.0
	Missing	55	94.8	62	92.5	21	84.0
	Regimen
	Monotherapy	27	46.6	17	25.4	17	68.0
	Dual therapy	12	20.7	26	38.8	7	28.0
	Triple therapy	14	24.1	21	31.3	1	4.0
	Quadruple therapy	5	8.6	3	4.5	0	0.0
	Outcome
	Hospitalization	51	87.9	60	89.6	24	96.0
	Death	0	0.0	2	3.0	1	4.0

Table 2 .

 2 Adjusted Reporting Odds Ratio of pancreatitis for antihyperglycemic drugs (French Pharmacovigilance Database, March 2008-March 2013) ROR adjusted for age, gender, history of pancreatitis, use of angiotensin converting enzyme inhibitors, angiotensin receptor antagonists, statins, codeine, morphine and related opioid agonists (other drugs associated with pancreatitis were not included in the model because of lack of events).

	Drug or pharmacological class	Cases of pancreatitis (N=147)	Non-cases (N=2962)	Adjusted ROR a	95% CI	p value
	All incretin-based drugs	122	568	15.62	9.81	24.87	<.0001
	All GLP-1 agonists	58	150	29.36	16.02	53.81	<.0001
	Exenatide	19	52	28.29	12.84	62.34	<.0001
	Liraglutide	39	99	30.36	15.36	60.01	<.0001
	All DPP-4 inhibitors	67	421	12.08	7.30	20.00	<.0001
	Sitagliptin	53	315	12.36	7.29	20.97	<.0001
	Vildagliptin	9	87	7.43	3.14	17.58	<.0001
	Saxagliptin	5	23	15.09	4.32	52.72	<.0001
	Metformin	84	1841	0.86	0.59	1.26	0.438
	Sulfonylureas/Glinides	53	1381	1.24	0.83	1.85	0.295
	Acarbose	2	137	0.41	0.09	1.89	0.255
	Thiazolidinediones	2	235	0.31	0.07	1.32	0.114
	Insulin	14	356	0.84	0.45	1.57	0.587
	a

Acknowledgements

The authors thank the French Network of Pharmacovigilance Centres and Pascal Auriche (ANSM) for access and extraction of data from the French Pharmacovigilance Database. The authors assume responsibility for the content and conclusions of this article.

Funding

This study was not founded.

Conflict of interest

None

Contribution statement

Pascal Auriche extracted the data from the French Pharmacovigilance Database; Jean-Luc Faillie designed the study, analyzed and interpreted the data and drafted the manuscript; Samy Babai, Sabrina Crépin, Virginie Bres, Marie-Laure Laroche, Hervé Le Louet, Pierre Petit, Jean-Louis Montastruc and Dominique Hillaire-Buys reviewed and revised the manuscript; all authors approved the version to be published.