Axel Parmentier

Frédéric Meunier

F R Éd

Éric Meunier

AIRCRAFT ROUTING AND CREW PAIRING: UPDATED ALGORITHMS AT AIR FRANCE

Keywords: Updated algorithms at Air France Airline management, column generation, cut generation, shortest path algorithm

come   L'archive ouverte pluridisciplinaire

1. Introduction 1.1. Context. Interactions between Operations Research and Air Transport Industry have been successful for at least five decades [5,[START_REF] Etschmaier | Operations research in the management of the airlines[END_REF]. These interactions have taken various forms: yield management, airplane timetabling, ground operations scheduling, air traffic management, etc. Key applications are notably the construction of sequences of flight legs operated by airplanes and crews. As airplane sequences of legs are routes and crew sequences pairings, this construction is called aircraft routing for airplanes, and crew pairing for crews.

The present paper focuses on these two applications and is the fruit of a research partnership with Air France, the main French airline. We aim at providing a reliable and easy to maintain framework that can cope with the specific and challenging industrial context of the company. Air France working rules are more complex than the IATA standards: collective agreements reach hundreds of pages, and two of the most cited references on crew pairing [START_REF] Desaulniers | Crew pairing at air france[END_REF][START_REF] Clarke | The aircraft rotation problem[END_REF] develop ad-hoc approaches to build pairings satisfying the company's rules. While the aircraft routing remains easy, this turns the exact resolution of Air France crew pairing into a challenge.

As crews need time to cross airports if they change airplane, the two problems are linked and the sequential resolution currently in use in the industry is suboptimal. Solving the integrated problem has been identified by academics as a difficult problem. Air France also requested an easy to maintain solution scheme for the integrated problem.

1.2. Literature review. All the versions of the problems considered in this paper are NPhard. We focus here on mathematical programming approaches.

Aircraft routing is considered at Air France as a pure feasibility problem, which contrasts with the recent literature which considers optimization versions. Authors either maximize profit when the fleet is heterogeneous [START_REF] Barnhart | Flight string models for aircraft fleeting and routing[END_REF][START_REF] Desaulniers | Daily aircraft routing and scheduling[END_REF], or minimize delay propagation along sequences of flights [28]. A recent paper introduces tools to deal with richer maintenance constraints [37]. State-of-the-art solution approaches rely on column generation [START_REF] Barnhart | Flight string models for aircraft fleeting and routing[END_REF][START_REF] Desaulniers | Daily aircraft routing and scheduling[END_REF][START_REF] Froyland | The recoverable robust tail assignment problem[END_REF]28], where columns are sequences of flight legs between airports where maintenance checks can be performed. They can solve to optimality large instances of the optimization versions in a few hours. (The solution proposed in [28] is actually for the so-called tail assignment problem, where the airplanes are distinguishable, but it can be adapted to aircraft routing.) Alternative approaches include heuristics [START_REF] Feo | Flight scheduling and maintenance base planning[END_REF] and Lagrangian relaxations [9]. Clarke et al. [9] propose a MIP with few variables. However, their MIP has two exponentialsize families of constraints: the first one is formed by classical subtour elimination constraints; the second one is formed by "minimal violation path" constraints that enforce the maintenance requirements. While the first family can be discarded for the aircraft routing problem we consider, the second one must be kept, and thus cut generation cannot be avoided to solve their MIP. Practically efficient compact integer programming formulations, that is, formulations which require neither column generation, nor cut generation, have recently been proposed by Cacchiani and Salazar-González [7] and Khaled et al. [25]. Compact formulations have the advantage to be more handy and can often be directly implemented in standard MIP solvers. Cacchiani and Salazar-González [7] consider the integrated problem, and propose a compact integer programming formulation for aircraft routing. Since they assume that the airplanes spend alternatively one night in a base and one night outside, their approach does not generalize to the Air France case, where one maintenance must be performed at least every four days, while the maintenance day is not fixed. By the way, due to the different maintenance requirements, the aircraft routing problem considered by Cacchiani and Salazar-González [7] has a polynomial status [22], while Air France problem is NP-hard [41]. Khaled et al. [25] have also recently proposed a compact MIP approach to tail assignment, which can be adapted to aircraft routing. They are able to solve to optimality instances with up to 1,178 legs and 30 airplanes in 3 hours. Their approach could in principle be used to address the version of the aircraft routing problem met at Air France. We discuss later in the paper the advantage of the approach we propose for aircraft routing with respect to theirs.

Crew variable wages and hotel rooms are among airlines first sources of variable costs. As both depend on the sequences of flight legs crews operate, crew pairing is an intensively studied optimization problem; see Gopalakrishnan and Johnson [21] for an extensive review. Regulatory agencies and collective bargaining agreements list numerous working rules that make the crew pairing problem highly non-linear and hence difficult. There is a long tradition of MIP approaches to crew pairing [2, 21] Since the seminal work of Minoux [START_REF] Minoux | Column generation techniques in combinatorial optimization: A new application to crew pairing[END_REF], state-of-theart approaches solve the crew pairing by column generation [START_REF] Anbil | A global approach to crew-pairing optimization[END_REF]3,[START_REF] Chu | Solving large scale crew scheduling problems[END_REF][START_REF] Desaulniers | Crew pairing at air france[END_REF][START_REF] Barnhart | An approximate model and solution approach for the long-haul crew pairing problem[END_REF][START_REF] Klabjan | Airline crew pairing generation in parallel[END_REF][START_REF] Klabjan | Airline crew scheduling with time windows and plane-count constraints[END_REF][START_REF] Clarke | The aircraft rotation problem[END_REF]40,[START_REF] Zeren | A novel column generation strategy for large scale airline crew pairing problems[END_REF]. They consider set partitioning formulations where columns are possible pairings. These methods hide the non-linearity in the pricing subproblem, which can be efficiently solved using resource constrained shortest path approaches [24]. As a large part of the working rules apply to duties, i.e., subsequences of a pairing formed by the flight legs operated on a same day, the subproblem is often split into two parts [START_REF] Desaulniers | Crew pairing at air france[END_REF][START_REF] Zeren | A novel column generation strategy for large scale airline crew pairing problems[END_REF]. The first one builds the set of all non-dominated duties. The second one builds the pairings by solving a path problem in the graph whose vertices are the non-dominated duties, and whose arcs are the pair of duties that can be chained. However, as the number of non-dominated duties is huge, solving the pricing subproblem becomes costly on large instances. When working rules are simple, one can also use compact integer programming approaches [6] where variables indicate if a given connection is used, and set partitioning formulations, where columns are the duties [42]. However, this is generally not the case, and such models are generally turned into initialization heuristics [3].

During the last decade, much attention has been devoted to the integration of aircraft routing and crew pairing. Moving from a sequential to an integrated approach enables to reduce the cost by 5% on average according to Cordeau et al. [11], and 1.6% according to Papadakos [START_REF] Papadakos | Integrated airline scheduling[END_REF]. Solution methods are column generation based heuristics [7,10,[START_REF] Mercier | An integrated aircraft routing, crew scheduling and flight retiming model[END_REF]32,[START_REF] Papadakos | Integrated airline scheduling[END_REF][START_REF] Salazar-González | Approaches to solve the fleet-assignment, aircraft-routing, crew-pairing and crew-rostering problems of a regional carrier[END_REF]39,43]. The heuristics of Cacchiani and Salazar-González [7], Salazar-González [START_REF] Salazar-González | Approaches to solve the fleet-assignment, aircraft-routing, crew-pairing and crew-rostering problems of a regional carrier[END_REF], Weide et al. [43] share many similarities with the ones we propose in this paper. Dunbar et al. consider robust [16] and stochastic [START_REF] Dunbar | An integrated scenario-based approach for robust aircraft routing, crew pairing and re-timing[END_REF] versions of the problem. To the best of our knowledge, the largest instances considered in the literature have 750 legs [43]. 1.3. Contribution and methods. The present paper is the result of a project initiated by Air France to design efficient and easy to maintain solution schemes for aircraft routing, crew pairing, and the integrated problem.

The first author has recently proposed an abstract framework [START_REF] Parmentier | Algorithms for non-linear and stochastic resource constrained shortest paths[END_REF] for computing resource constrained shortest paths. The main contribution of the present paper is the proof that this framework can be used on a concrete problem and considerably improves the size of the instances that can be solved at optimality. Indeed, we apply this framework to the pricing subproblem of a standard column generation approach for the crew pairing problem and solve to optimality from a few minutes to a few hours instances with up to 1, 000 flight legs, which outperforms previous performances on that problem. One key element in the performances of this framework is the use of sets of bounds to discard paths, instead of single bounds: this is useful in a context where any two resources are not necessarily comparable (in Section 4, further details will be given). Even if this idea of sets of bounds is present in the aforementioned paper of the first author, the present paper is the first proof that such a technique is very efficient in practice. We finally emphasize that the framework for shortest path computation does not explain how to model concrete problems like the one met for crew pairing. The modeling we propose is thus also a contribution on its own: more than 70 rules have to be satisfied, and most of them are non-linear. Finally, Desrosiers and Lübbecke [14, p. 16] underline that, in a column generation context, "accelerating the pricing algorithm itself usually leads most significant speeds-up". As all the approaches to the integrated problem use a column generation approach for the crew pairing, we believe that these approaches can be significantly accelerated by using our improved crew pairing pricing subproblem algorithm.

Our second contribution is a simple and compact integer formulation for aircraft routing. Such formulations are desirable in an industrial context since, as mentioned in the literature review, they do not require tricky development and can often be directly implemented in off-the-shelf solvers. In addition, our formulation is very efficient: it enables to solve all Air France industrial instances in at most a few minutes. As mentioned in the literature review, other compact formulations have recently been proposed [7,25]. Even though the formulation of Khaled et al. [25] could in theory be adapted to Air France specific problem, such an adaptation is not straightforward due to the fact that routes are cyclic in Air France aircraft routing problem. Furthermore, on Air France problem, our formulation admits a stronger linear relaxation than theirs (we discuss it in Appendix B), which gives a clear competitive advantage to our approach. We do not claim that our formulation outperforms their one on other versions of the problem such as the one they consider.

Finally, we design a simple cut generating method for solving the integrated problem, which relies on our contributions for aircraft routing and crew pairing. Like the one of Weide et al. [43], our method consists in solving alternatively crew pairing and aircraft routing problems. However, they do not consider the same problem: their version includes a notion of robustness with respect to delay. Experiments show that the method is able to solve to near optimality instances with up to 1, 766 flight legs, which again outperforms previous results on that problem. Due to the specificity of Air France problem, with no aircraft routing costs, our method cannot handle all the problems considered in the literature. However, it is the only one that proves optimality gaps smaller than 0.01% on instances with more than 600 legs.

We emphasize that for the three problems, our solution is easy to use and to maintain by the company. The algorithm for computing the shortest paths is already implemented and can be used as a black-box. The only non-trivial task is the modeling of the rules in the framework, but once a few techniques have been understood (like the ones we use later in the paper, in Section 4.2), even this step is straightforward.

1.4. Organization of the paper. Each of the three problems considered in the paper is addressed in a separate section: the aircraft routing problem is studied in Section 2, the crew pairing problem in Section 3, and the integrated problem in Section 5. Each of these sections gets exactly the same structure. It starts with a subsection describing the problem. A second subsection is then devoted to a modeling of the problem (e.g., the compact integer formulation for the aircraft routing problem). It ends with a subsection explaining the proposed method to solve the problem (e.g., column generation for the crew pairing problem and cut generation for the integrated problem).

Experiments showing the efficiency of the methods proposed in each of these three sections are provided and discussed in Section 6.

Section 4 is a section making a focus on an algorithmic subroutine required by our method for the crew pairing problem. This section is much more technical than the others and can be safely skipped at first reading (and the same holds for Section 6.7 that deals with specific experiments regarding this subroutine). This subroutine is the algorithm solving the pricing subproblem of the column generation. It relies on the shortest path framework of the first author and on bound sets, both described in that section.

The paper ends with a short conclusion (Section 7). All proofs are postponed to Appendix A. The input is formed by a set of airports, a collection L of flight legs, and a number n a of airplanes. Some airports are bases in which maintenance checks can be performed. A flight leg is characterized by departure and arrival airports, as well as departure and arrival times (it is of course assumed that departure time is smaller than arrival time for any flight leg). We consider the flight legs on a weekly horizon: the departure and arrival times are given for a typical week.

Compact integer program for aircraft routing

The purpose of aircraft routing is to determine routes for airplanes so that each flight leg is operated by an airplane each week without using more than n a airplanes. In addition, there are maintenance operations that have to be regularly performed: each airplane has to spend a night in a base at least every ∆ maint days, where ∆ maint is a given parameter, which is equal to 4 at Air France.

Formally, an airplane connection is a pair ( , ) of flight legs which satisfies

• the arrival airport of is the departure airport of • the duration between the departure time of and the arrival time of is bounded from below by a fixed quantity (which can depend on the airport, the time, and the fleet).

We underline that there are connections ( , ) with at the end of the week and at the beginning of the (next) week. A route is a cyclic sequence of distinct flight legs 1 , . . . , k such that any two consecutive flight legs ( , ) = ( i-1 , i ) or ( k , 1 ) is an airplane connection. Routes can last several weeks, but each week, the sequences of flight legs operated by airplanes are the same. In other words, when we consider all airplanes as indiscernible, the solution must have a week periodicity. Figure 1 illustrates two routes, which last respectively one and two weeks. As each flight leg has to be operated each week, routes lasting p weeks require p airplanes. A route satisfies the maintenance requirement if an airplane following this route in a cyclic way (repeating the solution when it reaches the end of the cycle) spends a night in a base at least every ∆ maint days The task consists in partitioning L into routes satisfying the maintenance requirement such that the number of airplanes needed to operate these routes is less than or equal to n a .

• • • • a • • • • c • • • • e • • • • g • • • • i • • • • k • • • • m • • • • o • • • • q • • • • b • • • • d • • • • f • • • • h • • • • j • • • • l • • • • n • • • • p • • • • r Figure 2.
Illustration of the digraph D of Section 2.2 for ∆ maint = 4. Only a few arcs of D (plain arcs) are displayed. The dashed cycle provides a route satisfying the maintenance requirement. The dotted path cannot be completed to a cycle since there is no outgoing arc from ( k , 4): it implies that there is no route of the form b , c , e , g , i , k , . . . that satisfies the maintenance requirement.

2.2. Modeling as an integer program. We first explain how the problem can be modeled as a disjoint cycle problem in a directed graph. This will make the description of the integer program a straightforward task. Define the directed graph D = (V, A) as follows. Its vertex set is L × [∆ maint ]. In other words, each flight leg is duplicated ∆ maint times. Each vertex ( , δ) corresponds to a flight leg ∈ L with the number of days δ ∈ [∆ maint ] since the last night spent in a base. An ordered pair (( , δ), ( , δ )) is in A if ( , ) is an airplane connection and we are in one of the three following situations:

• and are performed during a same day and δ = δ , as illustrated between legs f and h on Figure 2, • and are not performed on the same day, the airport is a base, and δ = 1, as illustrated between legs j and l , • and are not performed on the same day, the airport is not a base, and δ -δ ≥ 0 is the number of days between the arrival of and the departure of , as illustrated between m and o . In other words, an arc corresponds to two flight legs that can be consecutive in a route, with the suitable restrictions on the number of days since the last night spent in a base. A cyclic sequence of legs 1 , . . . , k satisfies the maintenance requirement if and only if there exists δ i for i in {1, . . . , k} such that ( 1 , δ 1 ), . . . , ( k , δ k ) is a cycle in D. Indeed, suppose that a route 1 , . . . , k satisfies the maintenance requirement, and denote by δ i the number of days since the last night spent in a base before i . Then the definition of D ensures that ( 1 , δ 1 ), . . . , ( k , δ k ) is a cycle in D. Conversely, suppose that a route does not satisfies the A 0 Figure 3. A two-week route crosses A 0 twice maintenance requirement. It contains a sequence of legs 1 , . . . , k spending at least ∆ maint successive days out of a base. Let j be the leg in that sequence before the ∆ maint th night. Then any path ( 1 , δ 1 ), . . . , ( j , δ j ) in D necessarily ends in vertex ( j , ∆ maint ), which has no outgoing arc, and there is no cycle in D corresponding to the route.

Figure 2 illustrates such a directed graph D. The route a , c , e , h , i , k , n , o , q satisfies the maintenance requirement and corresponds to the dashed cycle in D. The route b , c , e , g , i , k , n , o , q does not satisfy the maintenance requirement. The dotted path is an attempt to make it a cycle in D, but since there is no arc outgoing from ( k , ∆ maint ), it is not possible.

We choose arbitrarily one instant in the week and we denote by A 0 the set of arcs (( , δ), ( , δ )) "crossing this instant", i.e., such that the time interval between the departure of (included) and the departure of (excluded) contains the instant. Define moreover V to be the set ( , δ)

∈ V : δ ∈ [∆ maint ] .
Proposition 1. Feasible solutions of the aircraft routing problem are in one-to-one correspondence with collections C of vertex disjoint cycles in D such that we have simultaneously (i) for each , exactly one cycle in C has a nonempty intersection with V , and this intersection consists of a single arc. (ii) C has at most n a arcs in A 0 .

Therefore, the aircraft routing problem is equivalent to deciding whether the following integer program has a feasible solution:

a∈δ -(v) x a = a∈δ + (v) x a ∀v ∈ V (AR.1) a∈δ -(V ) x a = 1 ∀ ∈ L (AR.2) a∈A 0 x a ≤ n a (AR.3) x a ∈ {0, 1} ∀a ∈ A. (AR.4)
Equation (AR.1) is the flow equation. Together with (AR.2), it ensures that the solution is composed of vertex disjoint cycles. Equation (AR.2) ensures that (i) is satisfied and Equation (AR.3) ensures that (ii) is satisfied.

Solution method.

The solution we propose is to implement directly the integer program (AR) in any standard MIP solver. Its number of constraints is |L|×∆ maint +|L|+1 and its number of variables is the number of airplane connections times ∆ maint . Program (AR) is therefore of tractable size, and current off-the-shelf solvers can solve industrial instances in a few minutes. See our experiments in Section 6.

3. Column generation approach to crew pairing 3.1. Problem formulation. Roughly speaking, the crew pairing problem is similar to the aircraft routing problem: instead of building sequences of flight legs for the airplanes (the routes), crew pairing requires to build pairings, which are sequence of flight legs operated by the crews. While routes are cyclic sequences, pairings are noncyclic sequences (and they are often quite short). The set of flight legs have to be partitioned into pairings, but the constraints are much more complicated. Before stating formally the crew pairing problem, we introduce some terminology.

A pair of flight legs ( , ) is a connection if it satisfies:

• the arrival airport of is the departure airport of • the departure time of minus the arrival time of is bounded from below by a fixed quantity (which can depend on the airport, the time, and the fleet). This quantity is in general different from the similar one for aircraft routing. If the arrival of and the departure of are on the same day, then it is a day connection. Otherwise, it is a night connection. If the duration of a night connection is smaller than a threshold, then it is a reduced rest. (This term is due to the rest taken by crews performing a night connection.)

A pairing is a sequence of distinct flight legs such that any two consecutive flight legs form a connection. The subsequence of a pairing formed by all flight legs operated during a same day is a duty.

To be feasible, a pairing 1 , . . . , k has to satisfy the following rules: (a) the period between the departure of 1 and the arrival of k spans at most 4 days, (b) 1 starts and k ends in one of the Paris airports, (c) each duty contains at most 4 flight legs. If a duty starts with a leg , and the night connection ( , ) that leads to is a reduced rest, then the number of legs of the duty is at most 3, (d) the total flying duration in a duty does not exceed F (t), where F is a given function and t is the time at which the first leg of the duty departs, as well as more than 70 other rules which encode the IR-OPS regulation of the European Aviation Safety Agency and Air France working rules. A pairing is long if it spans 4 days. A duty is long if it contains more than 3 flight legs, and short otherwise. We denote by P the set of feasible pairings Operating a pairing p in P has a cost c p that corresponds to crew wages and hotel nights. Given a set L of flight legs, solving the crew pairing problem consists in selecting a collection of feasible pairings of minimum total cost so that each leg in L belongs to exactly one of them, and so that the following global constraints are satisfied: the proportion of long pairings in the solution is less than or equal to a quantity α, and the proportion of long duties is at most a quantity β.

Our purpose here is to introduce the main modeling ideas and not to get into details of the intricacies of the IR-OPS and Air France regulations. In the rest of the paper, we therefore present these ideas on a simplified problem with only the four illustrating rules (a), (b), (c), and (d). All the other rules are also taken into account in the numerical results.

3.2. Modeling as an integer program. As pairings must satisfy many non-linear rules such as rule (c), it is difficult to model crew pairing using a compact integer program that has both a good linear relaxation and a tractable size. The literature therefore generally uses a column generation approach where rules complexity are hidden in the set of variables (and dealt with using ad-hoc algorithms in the pricing subproblem). We also use a column generation approach.

The binary variable y p indicates if a pairing p in P belongs to the solution.

(CP)

min p∈P c p y p s.t. p y p = 1 ∀ ∈ L p∈P l y p ≤ α p∈P y p p∈P (1 -β)∆ l (p) -β∆ s (p) y p ≤ 0 y p ∈ {0, 1} ∀p ∈ P,
where p means that the flight leg is present in p, where P l is the set of long pairings, and where ∆ s (p) (resp. ∆ l (p)) is the number of short (resp. long) duties in a pairing p. The first constraint ensures that each leg is covered, the second that the proportion of long pairings is less than or equal to α, and the third that the proportion of long duties is less than or equal to β.

Column generation approach.

We propose an exact method for solving the program (CP). It is based on column generation. We describe the method without assuming special knowledge in column generation. We will in particular be sketchy on the theoretical rationale; more details on that topic can be found for instance in a survey by Lübbecke [START_REF] Marco E Lübbecke | Column generation[END_REF].

Algorithm 1 describes our column generation approach, which maintains a subset of pairings P ⊆ P. The idea of column generation is to solve the master problem, which is the linear relaxation of (CP) on such a subset P , and to check if by chance the optimal solution found on this restricted version is also an optimal solution of the master problem with the full set P. This checking is done exactly as in the classical simplex algorithm: it is the optimal solution for the full problem if all reduced costs are nonnegative. Since the number of elements in P is huge, it is not possible to compute and check all these reduced costs one by one. However, given an element p in P, it is always possible to compute its reduced cost from the value of the dual variables by standard linear programming theory. To find the element p with the smallest reduced cost, an auxiliary optimization problem instantiated by the values of the dual variables is solved: the pricing subproblem. In Algorithm 1, this is done in Step 5. The exact method to solve the pricing subproblem is described in Section 4. Right before Step 10, the linear relaxation of (CP) is fully solved.

Algorithm 1 Column generation algorithm

1: initialize P in such a way that (CP) restricted to P is feasible (e.g., taking all possible pairings of two flight legs makes the job); 2: repeat 3:

solve the linear relaxation of (CP) restricted to P with any standard solver; end if 9: until (c p ≥ 0 for all p ∈ P) 10: solve (CP) restricted to P with any standard solver; 11: denote by c upp its optimal value; 12: add to P all pairings with reduced cost (from the last linear program of Step 3) nonlarger than c upp -c low ; 13: solve (CP) restricted to P with any standard solver; 14: return its optimal solution y * ; Since the total number of possible pairings is finite, Step 3 -which consists in solving the master problem -is repeated only finitely many times, and thus the overall method terminates in finite time. After having performed Step 3 for the last time, c low is a lower bound on the optimal value of (CP).

Step 10 provides a first feasible solution of (CP). The value c upp (Step 11) is thus an upper bound on the optimal value of (CP). At that time, we have thus a lower bound on the optimal value, and a feasible solution.

The purpose of the remaining steps is to "close the gap". The idea consists in generating all pairings p that might be in an optimal solution. These pairings are precisely those whose reduced cost is smaller than or equal to c upp -c low . This is completely formalized by Lemma 2 below. Finding all these pairings in Step 12 is a variant of the pricing subproblem, briefly discussed in Remark 2 of Section 4. At the end, the solution y * is an optimal solution of (CP).

Lemma 2 (Nemhauser and Wolsey [34, Proposition 2.1, p. 389]). Consider an integer program in standard form with variables (z i ) for which the linear relaxation admits a finite optimal value v. Suppose given an upper bound UB on the optimal value of the integer program. Then for every i such that ci > UB -v, the variable z i is equal to 0 in all optimal solutions of the integer program, where ci denotes the reduced cost of the variable z i when the linear relaxation has been solved to optimality. Remark 1. The last steps of Algorithm 1 works only if the gap c upp -c low is small, as otherwise a huge number of pairings might be added at Step 12. We found out numerically that it is the case for all Air France instances, and we thus use this technique. Lemma 2 has been recently used by Cacchiani and Salazar-González [7] on the integrated problem. They underline that when c upp -c low is large, a branch and bound approach is required.

Pricing subproblem

We now introduce a solution scheme for the pricing subproblem (2) min p∈P cp . As pairings can be considered as paths satisfying constraints in the graph whose vertices are the legs and arcs the connections, the pricing subproblem is generally solved as a resource constrained shortest path problem, and we do not depart from this approach. We model it within the Monoid Resource Constrained Shortest Path Problem framework [START_REF] Parmentier | Algorithms for non-linear and stochastic resource constrained shortest paths[END_REF], which we now briefly describe. This framework is rather abstract, but practically, it only requires to implement a few operators on the resource set. This work is the first application of the Monoid Resource Constrained Shortest Path Problem framework to an industrial problem.

4.1. Framework and algorithm. A binary operation ⊕ on a set M is associative if q ⊕ (q ⊕ q ) = (q ⊕ q ) ⊕ q for q, q , and q in M . An element 0 is neutral if 0 ⊕ q = q ⊕ 0 = q for any q in M . A set (M, ⊕) is a monoid if ⊕ is associative and admits a neutral element. A partial order is compatible with ⊕ if the mappings q → q ⊕ q and q → q ⊕ q are non-decreasing according to this order for all q in M . A partially ordered set (M, ) is a lattice if any pair (q, q ) of elements of M admits a greatest lower bound or meet denoted by q ∧ q , and a least upper bound or join denoted by q ∨ q . A set (M, ⊕, ) is a lattice ordered monoid if (M, ⊕) is a monoid, (M, ) is a lattice, and is compatible with ⊕.

Given a digraph D = (V, A), a lattice ordered monoid (M, ⊕, ), elements q a ∈ M for each a ∈ A, origin and destination vertices o and d, and two non-decreasing mappings c : M → R and ρ : M → {0, 1}, the Monoid Resource Constrained Shortest Path Problem seeks an o-d path P of minimum c a∈P q a among those satisfying ρ a∈P q a = 0, where a∈P is always performed in the order of the arcs on the path P (the operation ⊕ is not necessarily commutative). We call such a q a the resource of the arc a. The sum a∈P q a is the resource of a path P , and we denote it by q P . The real number c (q P ) is its cost, and the path P is feasible if ρ (q P ) is equal to 0. We therefore call c and ρ the cost and the infeasibility functions.

We now describe an enumeration algorithm for the Monoid Resource Constrained Shortest Path Problem. It follows the standard labeling scheme [24] for resource constrained shortest paths. The specificity of our algorithm is that it uses, for each v in V a set B v of bounds such that,

(3) for each v-d path Q, there is a b ∈ B v with b q Q .
The lattice ordered monoid framework enables to design procedures to build these sets of bounds; see Section 4.3 for more details. Having defined these bounds, we define key(P ) as ( 4)

key(P ) = min{c(q P ⊕ b) : b ∈ B v , ρ(q P ⊕ b) = 0}
where v is the last vertex of P .

The empty path at a vertex v is the path with no arcs starting and ending at vertex v. By definition of paths resources, its resource is the neutral element of the monoid. A path P dominates a path Q if q P q Q . During the algorithm, a list L of partial paths, an upper bound c U B od on the cost of an optimal solution, and lists (L nd v ) v∈V of non-dominated o-v paths are maintained. Algorithm 2 states our algorithm. We denote by P + a the path composed of a path P followed by an arc a. if v = d, ρ(q P ) = 0, and c(q P ) < c U B od then 9:

c U B od ← c(q P ); 10:

else (extension of P )

11:

for all a ∈ δ + (v) do 12:

Q ← P + a;

13:

w ← last vertex of Q;

14:

if ∃b ∈ B w such that ρ(q Q ⊕ b) = 0 and c(q Q ⊕ b) < c U B od then 15:
if Q is not dominated by any path in L nd w then 16:

L nd w ← L nd w ∪ {Q} and remove from L nd w and L every path dominated by Q;

17:

L ← L ∪ {Q}; The specificity of our approach lies in the use of the bounds of Equation (3) in the algorithm. While it is well known that the use of lower bounds is a key element in the performance of the enumeration algorithms [START_REF] Dumitrescu | Improved preprocessing, labeling and scaling algorithms for the weight-constrained shortest path problem[END_REF], our approach is the first to allow the use of bounds with non-linear constraints such as (c) and (d). Not only these bounds are used to discard more paths, but they are also used to improve the order in which the paths are considered by the algorithm. The two main resource constrained shortest path algorithms in the literature [24] differ by the order in which one they consider paths. The label correcting algorithm is obtained from our one by using c(q P ) as key(P ) and removing the test of Step 14. It is for instance described by Dunbar et al. [16,[START_REF] Dunbar | An integrated scenario-based approach for robust aircraft routing, crew pairing and re-timing[END_REF] in the context of crew pairing. The label setting algorithm considers vertices v in a topological order, and then apply Steps 8 to 21 for each path P in L nd v . Again, the test of Step 14 is removed.

Finally, as is has already been noted, we go further by using sets of bounds rather than singletons (we give additional explanations in Section 4.3).

Remark 2. Step 12 of Algorithm 1 requires to solve the following variant of the Monoid Resource Constrained Shortest Path Problem: generate all the o-d paths P satisfying ρ(q P ) = 0 and c(q P ) ≤ c upp -c low .

Algorithm 2 can be easily adapted to this variant. It suffices to maintain a set S of solutions (initially empty), to replace c U B od by c upp -c low in Steps 8 and 14, to replace Step 9 by S ← S ∪ {P }, and to return S. The set S returned contains all the o-d paths P satisfying ρ(q P ) = 0 and c(q P ) ≤ c upp -c low . Remark 3. The terminology "label setting" and "label correcting" varies in the literature. We stick here to Irnich and Desaulniers [24]. What Dunbar et al. [16,[START_REF] Dunbar | An integrated scenario-based approach for robust aircraft routing, crew pairing and re-timing[END_REF] call a "label setting" algorithm is a label correcting algorithm according to Irnich and Desaulniers [24]. 4.2. Modeling the pricing subproblem. We now explain how to model our pricing subproblem (2) in the Monoid Resource Constrained Shortest Path Problem framework. As already mentioned, we only consider rules (a) to (d) (Section 3.1) to focus on ideas rather than on the full intricacies of the regulation. In other words, the set of feasible pairings is the set of pairings satisfying rules (a) to (d). For simplicity, we also omit in the master problem the long pairings and long duties constraints, and we assume that the cost c p is of the form ( , )∈p c ( , ) . We emphasize that all IROPS and Air France rules, as well as the real costs and the long pairings and long duties constraints, are taken into account in the numerical experiments. The reduced cost is then of the form cp = ( , )∈p c ( , ) + ∈p z , where z is the dual variable associated to the partitioning constraint.

We now model this toy subproblem as a Monoid Resource Constrained Shortest Path Problem. According to Section 4.1, we have to describe the digraph with its origin and destination vertices, the lattice ordered monoid, the resources on the arcs, and the cost and infeasibility functions. We actually solve a shortest path problem for each sequence of four consecutive days in a week, in order to satisfy rules (a) and (b) (Section 3.1). We thus solve seven Monoid Resource Constrained Shortest Path Problem instances per iterations of the pricing subproblem.

4.2.2.

The lattice ordered monoid. The monoid M we use for the resources is of the form

M ρ × R, where M ρ = (Z + × R + ) ∪ (Z + × R + ) 2 ∪ {∞}.
An element (n, f ) ∈ Z + × R + models the resource used over a single day by a pairing: n flight legs and a total flying duration f . An element (n b , f b , n e , f e ) ∈ (Z + × R + ) 2 models the resource used over the first and last days of a pairing lasting more than one day. Quantities n b and f b are the number of flight legs and flying duration on the first day of a pairing, and n e and f e are the number of flight legs and flying duration on the last day of a pairing. The element ∞ is used to capture infeasibility of certain pairings.

Let F m = max t F (t). We define the operator ⊕ on M ρ as follows.

r ⊕ ∞ = ∞ ⊕ r = ∞ for all r ∈ M ρ (n, f ) ⊕ (ñ, f ) = (n + ñ, f + f ) (n, f ) ⊕ (ñ b , f b , ñe , f e ) = (n + ñb , f + f b , ñe , f e ) (n b , f b , n e , f e ) ⊕ (ñ, f ) = (n b , f b , n e + ñ, f e + f ) (n b , f b , n e , f e ) ⊕ (ñ b , f b , ñe , f e ) = ∞ if n e + ñb > 4 or f e + f b > F m , (n b , f b , ñe , f e ) otherwise.
We define on M ρ by (0, 0) q and q ∞ for all r ∈ M ρ

(n, f ) (ñ, f ) if n ≤ ñ and f ≤ f (n b , f b , n e , f e ) (ñ b , f b , ñe , f e ) if n b ≤ ñb , f b ≤ f b , n e ≤ ñe , and
f e ≤ f e ,
and a pair (n, f ) = (0, 0) is not comparable with (n b , f b , n e , f e ).

Lemma 4. (M ρ , ⊕, ) is a lattice ordered monoid.

As (R, +, ≤) is a lattice ordered monoid, the monoid M = M ρ × R is a lattice ordered monoid when endowed with the componentwise sum and order.

4.2.3.

Resources on the arcs. Consider an arc ( , ) of D. If it is a day connection, then we define its resource to be (1, f ( )), c ( , ) + z , where f ( ) is the flying duration of leg , and z is the dual variable of the cover constraint associated to in (CP). If it is a night connection, then we define its resource to be (0, 0, n e , f e ), c ( , ) + z , where n e = 2 if ( , ) is a reduced rest, and 1 otherwise, and f e = f ( ) + F m -F (t), where t is the departure time of . Similarly, each arc (o, ) has resource (0, 0, 1, f e ), z , and each arc ( , d) resource (0, 0, 0, 0), 0 . 4.2.4. Cost and infeasibility functions. Given q = (r, z) ∈ M , we define Proposition 5. The sequence of flight legs p corresponding to an o-d path P is in P (i.e., is a feasible pairing) if and only if ρ(q P ) = 0. In that case, c(q P ) = cp .

ρ ((r, z)) = ρ M ρ (r) and c((r, z)) = z where ρ M ρ is defined on M ρ by ρ M ρ (n, f ) = max 1 (4,∞) (n), 1 (Fm,∞) (f ) , ρ M ρ (n b , f b , n e , f e ) = max 1 (4,∞) (n b ), 1 (Fm,∞) (f b ), 1 (4,∞) (n e ), 1 (Fm,∞) (f e ) , ρ M ρ (∞) = 1,
Q (2, 0) (0, 2) (2, 2) R (1, 1) (1, 1) R R b. q 1 q 2 × q R × q R × q R • • q R ∧ q R ∧ q R q R ∧ q R
Appendix C details how this reduction works on a small pricing subproblem instance and shows a few typical iterations of Algorithm 2.

Bounds on resources.

We give now a simple illustration of why sets of bounds enable to discard more paths than single bounds. Consider the example on Figure 4.a, where we have an o-w path Q, and three w-d paths R, R , and R . Resources, which belong to R 2 endowed with the componentwise sum and order, are indicated on Figure 4.a. The resources of w-d paths are indicated by crosses on Figure 4.b. Consider a situation where ρ (q 1 , q 2 ) is equal to 1 if and only if max(q 1 , q 2 ) > 2. There is no feasible o-d path starting by Q.

Suppose first that we are using single bounds as in the usual approach, i.e., B w contains a unique element b w . Recall that b w is then such that b w q S for every w-d path S. In such a case, b w q R ∧ q R ∧ q R = (1, 1). Hence, ρ(q Q ⊕ b w ) = 0 and the path P is not discarded at Step 14 of Algorithm 2. Suppose instead that we use the two bounds b 1 = q R ∧ q R and b 2 = q R , which are indicated by circles on Figure 4.b. We have then ρ(q Q ⊕ b 1 ) = ρ(q Q ⊕ b 2 ) = 1, and the path Q is discarded at Step 14. This example is very simple, but it is the same mechanism that is in work in the general case.

The first author [START_REF] Parmentier | Algorithms for non-linear and stochastic resource constrained shortest paths[END_REF] introduced a procedure which, given a size κ in input, builds lower bounds sets B w of size κ. Larger sets of bounds B w enable to get larger lower bounds, and hence to discard more paths. However, larger sets of bounds also mean a longer preprocessing is required to compute the bounds. Hence, the parameter κ is chosen to obtain a tradeoff between the quality of the bounds and the time needed to compute them. Regarding the practical choice of κ for the crew pairing pricing subproblem, the following rule of thumbs ensures good results in practice: use κ = 1 if there are fewer than 100 vertices, κ = 50 is there are fewer than 300 vertices, κ = 150 if there are fewer than 1, 500 vertices, and κ = 250 if there are more.

When solving (CP), the preprocessing, which actually consists in building an "extended graph", is done once and for all: the same extended graph is used each time we solve the pricing subproblem. We can thus work with larger sets of bounds than independent resolutions of the Monoid Resource Constrained Shortest Path Problem would have allowed.

Integrated problem

5.1. Problem formulation. If a crew changes airplane during a connection between two flight legs and , its members need time to cross the airport between the arrival of and the departure of . This is not possible if the time between the arrival of and the departure of is too short. A short connection is an ordered pair ( , ) of flight legs that can be operated by a crew only if and are operated by the same airplane. Due to short connections, aircraft routing and crew pairing are linked. Given the collection of all short connections, the integrated problem consists in finding a solution of the aircraft routing problem of Section 2 and a solution of the crew pairing problem of Section 3 such that whenever a short connection is used in a pairing, it is also used in a route of an airplane.

Modeling as an integer program.

The solutions x of (AR) and y of (CP) provide a solution of the integrated problem if an only if (5) p∈Pα y p ≤ a∈Aα x a for every short connection α = ( , ), where we denote by A α (resp. P α ) the sets of arcs (resp. pairings) using the short connection α. For any feasible solution of the aircraft routing problem, there is a solution of the crew pairing problem compatible with it since there is no constraint on the number of crews, but solving the two problems simultaneously allows to spare these additional crews and to reduce the costs, as explained in Section 1.2. The integrated problem aims at performing this task and is thus modeled by the following integer program (Int) min p∈P c p y p s.t. x satisfies constraints of (AR) y satisfies constraints of (CP)

x and y satisfy constraints (5) for all short connections α.

5.3.

A cut generating approach. As we will see in the numerical results, the aircraft routing and crew pairing solution schemes introduced solve most of our industrial instances to optimality in a few hours. It is therefore natural to test the ability of a simple combination of these approaches to tackle with the integrated problem. Instead of solving directly Program (Int), we adopt a cut generating approach using the methods proposed in the previous sections in a rather independent way. Let S(y) denote the set of short connections used in a solution y of (CP). Given a feasible solution y of (CP), if there is no feasible solution x of (AR) satisfying (5), then any solution y such that S(y) ⊆ S(y ) leads to a more constrained (AR), and hence to a similar infeasibility. To avoid such solutions in (CP), we set S = S(y) and add the constraint (6)

p∈P |p ∩ S|y p ≤ |S| -1,
where |p ∩ S| denotes the cardinality of {α ∈ S : p ∈ P α }. It prevents a solution to use all short connections in S but does not restrict otherwise the set of solutions.

We can now describe the algorithm for the integrated problem. The algorithm maintains a set S of short connection cuts. Initially, S is empty. The following steps are repeated.

(i) Solve (CP) with additional constraints (6) for S ∈ S. Let y * be the optimal solution. (ii) Solve (AR) with the additional constraints (5).

• If it is feasible, then stop (we have found the optimal solution of (Int)).

• Otherwise, add S(y * ) to S and go back to (i).

Because of the cuts added along the algorithm, a solution y * is considered at most once. The number of solutions to the crew pairing problem being finite, the cut generation algorithm terminates after a finite number of iterations. The solutions of the last call to (i) and (ii) form an optimal solution to (Int): at each iteration, the only solutions to (CP) that are forbidden by the additional constraints (6) are not feasible for (Int) and at the last iteration, y * is the optimal solution of a relaxation of (Int).

In practice, the algorithm does not converge after thousands of iterations on industrial instances. We therefore replace |S| -1 by γ|S| with γ < 1 in the constraints (6), losing the optimality of the solution returned. During the first iteration, there is no additional constraint (6): the crew pairing problem (CP) is therefore not constrained by the aircraft routing. It is therefore a relaxation of the integrated problem, and its optimal solution provides a lower bound on the optimal solution of the integrated problem. We use this lower bound to evaluate the quality of the solution of the integrated problem returned by the algorithm. Numerical experiments in Section 6 show that γ = 0.9 is a good compromise: we obtain near optimal solutions after a few dozens of iterations.

Experimental results

6.1. Instances. Table 1 describes six industrial instances of Air France. Each instance contains the legs of a fleet on a weekly horizon. The two first columns provide the name of the instance and the number of legs it contains. Columns "Airplane connect." and "Airplanes" respectively give the number of connections that can be done by airplanes, i.e., the number of ordered pair of legs ( , ) that can be operated consecutively in a route, and the number n a of airplanes available. Columns "Crew connect". and "Crew pairings" respectively provide the number of connections that can be taken by crews, and the order of magnitude of the number of pairings in a good solution. Finally, column "Short connect". gives the number of short connections available. These instances are large: for instance, the largest instance considered by Mercier et al. [32] has 707 legs, and the largest instance for the integrated problem in the literature [43] have 750 legs.

The instance A318-9 (resp. A320-fam) contains the legs of the A318 and A319 (resp. A318, A319, A320, and A321) instances, as well as a few extra "fictitious" legs. As Air France crews can operate legs on planes of different subfleets on the same pairing, there is a common crew pairing problem for each of the instance A318-9 and A320-fam. On the contrary, the subfleet of an airplane is fixed: an A318 does not become an A319. Hence, solving aircraft routing for multiple subfleets together consists in solving one separate problem for each subfleet. This is what we do when we solve aircraft routing instances within the integrated problem solution scheme for instances A318-9 and A320-fam. 2 provides the results for aircraft routing. The first column gives the name of the instance. The next two ones give results for the aircraft routing problem on its own. Column "Uncons. CPU time (mm:ss)" the time needed to solve (AR). Column "Optim. CPU time (mm:ss)" gives the time needed to find an optimal solution of the optimization problem that consists in finding the minimum number of airplanes needed to operate the instance: we use the left-hand side of (AR.3) as objective. Note that this problem has not been mentioned previously in the paper. The solution scheme for the integrated problem in Section 5 solves (AR) with additional constraints (5). The two last columns provide numerical results for this constrained version. On all but the last iterations of the integrated problem scheme, aircraft routing is infeasible. Column "Infeas. CPU time" provides the time needed to solve the penultimate iteration, which is infeasible, and column "Feas. CPU time " the last iteration, which is feasible. The typical computing time is a few dozens of seconds on industrial instances. The longest constrained feasible version requires a few minutes. The optimization version in the second column is typically one order of magnitude faster than the one obtained by Khaled et al. [25, Tables 10 and11] on instances with similar number of legs. However, this last statement must be taken with care as the structure of the instances (number of airplanes, number of days of planning, etc.) is very different. This improved performance is likely due to the fact that our relaxation has a better relaxation than their one, as we prove in Proposition 7 in Appendix A and that our formulation has less symmetry. 3. Crew pairing results -Instances are solved to optimality 6.4. Crew Pairing. Table 3 provides the results for crew pairing. All instances are solved to optimality. The first column of Table 3 gives the name of the instance. The next column provides the value of κ determined using the rule of thumb of Section 4.3 and needed by the algorithm building the sets B v . Column "Col. Gen. Iter" provides the number of iterations in the column generation, and column "Pricing time" the percentage of time spent in the pricing subproblem. This pricing time includes the time needed by the computation of the sets B v and by the enumeration algorithm (Algorithm 2). Columns "LP time" and "MIP time" indicate the percentage of the total CPU time spent in Algorithm 1 solving Step 3, and solving Steps 10 and 13. The last column gives the total time needed by the algorithm. On all these instances, the integrality gap does not exceed 0.01%. This explains the fast resolution of Step 13.

Remark 4. One may be tempted to stop the column generation before convergence to exchange quality for speed. Unfortunately, and this is a limit of our method, if we stop the column generation before convergence, the solution found by the MIP solver at Step 10 is poor, and Step 12 is not tractable in practice. 6.5. Integrated problem. Table 4 provides the results for the integrated problem. The constraint strength parameter γ of the end of Section 5 is equal to 0.9, and the bounds sets size κ is equal to 150. The first column provides the instance solved. Columns "Integ. steps" provides the number of steps of the integrated problem algorithm of Section 5 before convergence. Column "CG it. total" provides the total number of column generation iterations realized on the successive integrated problem algorithms steps. Column "(CP) CG time" provides the proportion of the total CPU time spent in the column generation, i.e., solving the pricing subproblem and the linear relaxation of the master problem, and column "(CP) MIP time" the proportion spent solving the integer version of the crew pairing master problem. The column "(AR) time" provides the proportion spent solving aircraft routing integer program (AR). The column "Sho. Con." gives the number of short connections in the final solution. The linear relaxation of the crew pairing master problem (CP) with no short connection constraint is used as the lower bound on the cost of an optimal solution. The gap provided is between the cost of the solution returned and this lower bound. Finally, the last column provides the total CPU time needed by the algorithm. Only instance A320-fam could not be solved, as the algorithm had not converged after one week of computing time.

We emphasize the fact that the solution returned by the approximate algorithm is almost optimal. Practically speaking, the gap obtained is less than or equal to 0.01%. The computation time needed to obtain a near optimal solution of the integrated problem is of the same 5. Sequential versus integrated resolution of aircraft routing and crew pairing order of magnitude than the time needed to obtain a solution of the crew pairing problem in Table 3. Solving aircraft routing and crew pairing sequentially strongly constrains the solution: indeed, when solved in an integrated fashion, around half of the connections in the solution are short connections.

Finally, Table 5 compares the sequential approach, where aircraft routing is solved first, and then crew pairing, to the integrated approach of Section 5. Line "Cost reduction" provides the percentage by which the cost is reduced when using the integrated approach, and line "CPU time ratio" the increase in computing time. On average, using the integrated approach enables to reduce the costs by 0.25%, and computing time is 5.1 times longer. This reduction of cost on our instances is smaller than what is mentioned in the literature: 5% on average according to Cordeau et al. [11], and 1.6% according to Papadakos [START_REF] Papadakos | Integrated airline scheduling[END_REF]. We believe that this comes from the fact that our instances are larger: adding new connections have a stronger impact when few connections are available. The increase in computing time is mainly due to the fact that the crew pairing is longer to solve in the integrated approach due to the addition of short connections. 6.6. Industrial relevance. To be usable in an industrial context, the computing time of the solvers must not exceed eight hours, which represent one night of computing time. Our algorithms enables to solve to near optimality instances of the integrated problem with up to 1, 766 legs within this time constraint. Our solution scheme therefore enables to deal practically with instances larger than those in the literature -the largest instances in the literature [43] have 750 legs. These performances have been made possible by our pricing subproblem algorithm. Further improvements to deal with larger instances cannot be done by only working on the pricing subproblem. Indeed, we can see in Table 3 that most of the CPU time on instance A320-fam is spent in the simplex algorithm. 6.7. Focus on the pricing subproblem. The key element in the performance of our approach is the performance of Algorithm 2. We compare it in this section to the algorithms previously used. As our industrial instances are too large to be solved using these algorithm 6 smaller instances, which we have built by considering only a subset of the legs of the instance A318. Columns "Legs" and "Connections" respectively provide the number of legs and connections in the instances, and column "Pairings" the approximate number of pairings in a solution. Table 7 provides results on the performance of the pricing subproblem scheme on these instances. Its first column gives the instance solved. The next one provides the algorithm used. Parameter κ, introduced in Section 4.3, gives the size of the lower bounds set for algorithms using bounds. The three next columns provide statistics on the resource constrained shortest path (RCSP) algorithms. As mentioned in Section 4.2, seven RCSP instances are solved for each pricing subproblem, one for each period of four consecutive days. The statistics are averaged on all the instances solved along the column generation. Column "RCSP iter av. nb" provides the average number of iterations of the RCSP algorithm, "Cut Dom." provides the proportion of paths cut at Step 15, the remaining being cut at Step 14. Column "RCSP time" provides the average time needed to solve one RCSP instance. Column "Pricing subproblem" provides the proportion of the total computing time spent solving the pricing subproblem, and the last column gives the total computing time of the crew pairing solution scheme. We can see that the use of bounds enables a huge speed-up with respect to the usual algorithms. This speed-up is required to deal with instances with more than 500 legs. Two elements explain this speed-up. First, the condition at Step 14 enables to discard many paths: in Algorithm 2, more than 90% of the paths discarded are discarded at Step 14 and not at Step 15. Second, running Algorithm 2 with c(q P ) instead of the min{c(q P ⊕ b) : b ∈ B v , ρ(q P ⊕ b) = 0} makes it much slower. Hence, min{c(q P ⊕ b) : b ∈ B v , ρ(q P ⊕ b) = 0} seems numerically to be a better evaluation of how P is promising than c(q P ), and enables Algorithm 2 to find good solutions faster than the label correcting algorithm, and also faster than the label setting algorithm that does not use keys at all for choosing the next path to consider.

Conclusion

We have proposed a compact integer program for aircraft routing. Its main strength is its ease of implementation. Numerical results show that it can deal with large industrial instances in at most a few minutes, even when optimization versions are considered. We have used a resource constrained shortest path algorithm recently introduced by the first author for the crew pairing column generation pricing subproblem. This algorithm leverages on the lattice ordered monoid structure of the resource set to build efficient lower bounds. Practically, this enables to solve to optimality very large industrial crew pairing instances. As a side result, we have combined these aircraft routing and crew pairing solution schemes in a cutting plane approach to the integrated problem. The resulting algorithm solves to near optimality large industrial instances of the integrated problem. The nonnegativity of the ci 's and the z * i 's implies the following inequality for every i ∈ N :

z * i ≤ UB -v ci .
In particular, if ci > UB -v, we necessarily have z * i < 1, i.e., z * i = 0.

Proof of Proposition 3. We first prove that the algorithm terminates after a finite number of iterations. By induction on the iterations, we see that there is never two distinct elements in L such that one is a subpath of the other: there is no such two elements when the algorithm starts, and if there is no such two elements at a given iteration, there cannot be such two elements at the next iteration either. In particular, when a path P leaves L at some iteration, it prevents the presence of a subpath of it in L at the current iteration. Thus, because of the update rule, P cannot be added again to L in a subsequent iteration. It implies that a given path P is considered at most once by the algorithm and, as there is a finite number of paths in an acyclic digraph, we get the sought conclusion regarding the termination of the algorithm.

We prove now the part of the statement regarding the cost of c U B od at the end of the algorithm. At any step of the algorithm, c U B od is either equal to +∞ or to the cost of an o-d path P such that ρ(q P ) = 0. Therefore, if there is no feasible solution, then c U B od is never updated, and equal to +∞ at the end of the algorithm. Suppose now that there is a feasible solution, and let P be a feasible o-d path of minimum cost. By definition of P and due to the update rule, we have c(q P ) ≤ c U B od at the end of the algorithm. Suppose for a contradiction that this inequality is strict. Given the update rule of c U B od , this means that neither P nor a feasible o-d path Q dominating P has been considered present in L. Let P be the longest subpath of P , with origin o, such that P , or a path dominating P with the same origin and destination as P , has been present in L. Denote by v the destination of P . Among all o-v paths present in L at some time during the algorithm and that dominates P , pick a path Q that is non-dominated by the others.

The test of Step 15 is necessarily satisfied by Q because, by definition of Q , there is no o-v path in L nd v dominating Q when Q is considered by the algorithm. We now prove that the test of Step 14 is also necessarily satisfied. Indeed, let P be the v-d subpath such that P = P + P , and b a bound in B v such that b q P . (Here, P + P means that we append P to P .) We have q Q ⊕ b q P ⊕ b q P ⊕ q P = q P . We therefore have, when Q is considered, ρ(q Q ⊕ b) ≤ ρ(q P ) = 0, and c(q

Q ⊕ b) ≤ c(q P ) < c U B
od , where the last inequality relies on the fact that c U B od is non-increasing along the algorithm. Hence, the test of Step 14 is satisfied. Therefore, Q satisfies the tests of Steps 14 and 15, and is added to L whatever the combination of these two tests the algorithm uses. Since Q is non-dominated by other o-v paths considered by the algorithm, it is extended in a subsequent iteration. Denote by a the arc in P that follows P and denote by w the head of a. Since q Q ⊕ q a q P ⊕ q a , either Q + a, or another o-w path that dominates P + a is present once in L, which contradicts the length maximality of P .

Proof of Lemma 4. Considering the different cases in the definitions enables to prove that (0, 0) is the neutral element of ⊕, that ⊕ is associative, that is compatible with ⊕, and that (M ρ , ) is a lattice with meet operator

q ∧ ∞ = q (n, f ) ∧ (ñ, f ) = min(n, ñ), min(f, f ) (n b , f b , n e , f e ) ∧ (ñ b , f b , ñe , f e ) = min(n b , ñb ) min(f b , f b ) min(n e , ñe ) min(f b , f b )
and q ∧ q = (0, 0) for any other combinations.

Proof of Proposition 5. Let p be a sequence of flight legs. The definition of D ensures that there is an o-d path P whose vertices correspond to the legs in p if and only if P is a pairing that satisfies rules (a) and (b). In that case, this path P is unique. Let p be such a pairing and P be the corresponding path. Let q P = q ρ P , c P be the resource of P . By definition of the arc resources, c P is the sum of c ( , ) + z for ( , ) in P and hence in p, and we therefore have c P = cp . Hence c(q P ) = cp .

Let n b , f b , n e , and f e be respectively the number of legs and the flying time of the first duty of p, and the number of legs and the flying time of the last duty of p. Let A P be the set of duties of P except the first and the last. We claim that q ρ P is equal to ∞ if and only if there is a duty in A P that does not satisfy both rules (c) and (d), and to (n b , f b , n e , f e ) otherwise. This result is proved by induction on the number of arcs in P . Denoting a the last arc of P , and P the subpath of P obtained by removing a, the induction hypothesis can be applied to P , and the result for P follows by considering the different possible cases for the components in M ρ of q P and q a . The definition of M ρ then ensures that ρ(q P ) = 0 if and only if rules (c) and (d) are satisfied. Hence, ρ(q P ) = 0 if and only if p is a feasible pairing.

Appendix B. Aircraft routing MIP

Khaled et al. [25] propose a compact MIP for the tail assignment problem. Aircraft routing and tail assignment both consist in building the sequences of legs operated by the airplanes of an airline. The main difference between them is that in aircraft routing, airplanes are identical and routes do not need to be assigned to airplanes, while in tail assignment airplane specific costs are taken into account and routes are assigned to airplanes. Aircraft routing and tail assignment problems being fairly similar, MIP formulations for one problem can generally be applied to the other one.

However, the first reason why we do not use Khaled et al.'s formulation on Air France problem is that their formulation does not naturally adapt to this problem. Indeed, an important difference between Khaled et al.'s tail assignment and Air France aircraft routing is that, in Khaled et al.'s problem, routes are not cyclic, while in Air France problem, they are. By cyclic routes, we mean that sequences of legs are built for a typical week, and that airplanes operate the same sequences of legs week after week in a cyclic way. And as we have seen in Section 2.1, one cycle can potentially last several weeks. If Khaled et al.'s compact MIP for tail assignment can easily be adapted to solve a "non-cyclic" version of the aircraft routing problem, this is not the case for the "cyclic" aircraft routing problem considered at Air France that we introduce in Section 2. Indeed, adapting it would require to introduce many new binary variables to encode how the end of a week cycles with the beginning of the next one.

The second reason is computational. As cyclic routes are difficult to take into account in their formulation, we compare the formulations on Air France tail assignment problem. Both Khaled et al.'s formulation and an adapted version of our formulation (AR) have been implemented at Air France, and the adapted version of our formulation is now used in practice by Air France to solve its tail assignment. Our formulation is able to solve Air France medium haul instances (the tail assignment version of instances A318, A319, A320, and A321 of Table 1) to optimality in at most 3 minutes, while the version of Khaled et al. is unable to solve these instances in two hours. On long-haul instances of Air France, which are easier, their formulation had been previously used and was able to find optimal solutions, but each instance took at least 16 minutes, while our formulation solves each instance in at most 20 seconds.

Actually, the better performances of our formulation are easily explained by theoretical considerations on the linear relaxations. Solvers of MIP are based on branch-and-bound, which crucially relies on the quality of the linear relaxation to discard partial solutions. On Air France tail assignment problem, the linear relaxation of our MIP provides bounds that are non smaller than those provided by the Khaled et al.'s MIP and can be strictly larger, even on very simple and natural examples. In the remaining of the appendix, we introduce Air France tail assignment problem, adapt Khaled et al.'s formulation and (AR) to that problem, and prove the result mentioned on linear relaxations.

We now introduce Air France tail assignment problem. The input is formed of a given week, a set of airports, a collection L of flight legs operated between these airports that week, and a set of n a available airplanes. Some airports are bases where maintenance can be performed, and each airplane must still spend a night in a maintenance base at least every ∆ maint days. Airplanes are indexed by j. For each airplane j in [n a ], let k j 0 be the airport where airplane j is at the beginning of the week, and δ j 0 be the number of days since the last maintenance night of j at the beginning of the week. For each airplane j and leg , we have a cost of operating a leg with airplane j. The aim of the tail assignment problem is to build the (non-cyclic) sequence of legs operated by each airplane at minimum cost.

A tail assignment connection is therefore a pair ( , ) of flight legs such that departs from the arrival airport of , and such that the departure time of minus the departure time of is bounded from below by a given quantity. We underline that, if there were airplane connections between a leg at the end of the week and a leg at the beginning of the week in the aircraft routing problem, there is no such tail assignment connection. A tail assignment route r for airplane j is a (non-cyclic) sequence of legs 1 , . . . , k such that i departs from k i 0 and any pair of two consecutive legs ( i , i+1 ) is a tail assignment connection. It satisfies the maintenance requirement if, first, supposing that airplane j follows this route, it spends its first night in a maintenance base after at most ∆ maint -δ j 0 days and then spends a night in a maintenance base at least every ∆ maint days, and second, if k arrives in an airport that is not a base, then the last night of the week spent in a base is at most ∆ maint -1 days before the end of the week. The cost of r is the sum of the costs of operating the legs in r with airplane j.

The task consists in building, for each airplane j, a route satisfying the maintenance requirement in such a way that each leg of is operated by one airplane j, and the sum of the costs of the routes is minimum.

We now generalize our MIP to Air France tail assignment. Let D = (V, A ) be the digraph with vertex set V = L × ∆ maint , and arc set A composed of pairs ( , δ), ( , δ ) such that ( , ) is a tail assignment connection and , δ, , and δ satisfy one of the three conditions defining the arcs of digraph D in Section 2.2. Digraph D is the analogue of digraph D of Section 2.2 where airplane connections are replaced by tail assignment connections. Contrary to digraph D, digraph D is acyclic as there is no connection between the end of the week and the beginning of the week. For each airplane j, let V j and Ãj be copies of V and A . We build a digraph D j as follows. Its vertex set is V j ∪ {s j , t j }, where s j is a source vertex, and t j a sink vertex. Its arc set is denoted by A j and contains Ãj as well as arcs

• s j , ( , δ) such that leg starts from airport k j 0 on day δ -δ j 0 if k j 0 is not a base, • s j , ( , 1) such that leg starts from airport k j 0 if k j 0 is a base, • ( , δ), t j such that leg ends on day d in an airport that is not a base, and δ < ∆ maint + d -8, • ( , δ), t j such that leg ends in a base. The number 8 in the third condition is just the number of days in a week plus 1. The digraph D j is acyclic. For each leg and airplane j, we denote by V j the vertices of V j of the form ( , δ). Given an arc a in A j , we define the cost c a to be equal to 0 if the tail of a is s j and to the cost of operating the leg of the tail vertex of a with airplane j otherwise. Given the definition of the digraph D j , the following proposition is immediate. Proposition 6. A sequence of legs 1 , . . . , k is a tail assignment route r for airplane j if and only if there exists δ 1 , . . . , δ k such that s j , ( 1 , δ 1 ), . . . , ( k , δ k ), t j is an s j -t j path P in D j . In that case, the cost of operating r with j is a∈P c a .

The following integer program therefore enables to model Air France tail assignment. min j∈[n a ] a∈A j c a x a (TA.1) 3)-( 6), ( 8), ( 9), (11), ( 12), (14), and (15) of [25] (KTA.2) d∈{1,...,∆ maint -δ j 0 } y j d ≥ 1 for all j in [n a ] (KTA.3) where F and P are the notations of [25] for the sets of flights legs and available airplanes respectively, and constraint (KTA.3) ensures that airplane j spends its first night in a base after at most ∆ maint -δ j 0 days. The next proposition shows that the linear relaxation of Program (TA) provides bounds that will discard more partial solutions than that of (KTA) in a branch-and-bound.

s.t. a∈δ + (s j ) x a = 1 ∀j ∈ [n a ] (TA.2) a∈δ -(v) x a = a∈δ + (v) x a ∀v ∈ V j , ∀j ∈ [n a ] (TA.3) j∈[n a ] a∈δ -(V j ) x a = 1 ∀ ∈ L (TA.4) x a ∈ {0, 1} ∀a ∈ A j , ∀j ∈ [n a ]. (
Proposition 7. The optimal value of the linear relaxation of (TA) is not smaller than the one of the linear relaxation of (KTA) and there are instances for which it is strictly larger.

Proof. We first prove that any feasible solution of the linear relaxation of Program (TA) can be turned into a feasible solution of the linear relaxation of (KTA). Consider a feasible solution x = (x 1 , . . . , x n a ) of the linear relaxation of Program (TA), where x j is the vector of variables x a with a in A j . As Equation (TA.3) defines the s j -t j flow polyhedron of the acyclic digraph D j , x j can be written as a conic combination P λ P χ P of indicator vectors of s j -t j paths in D j , where the sum is taken over all s j -t j paths and the λ P are non-negative. Equation (TA.2) gives P λ P = 1.

Let R j be the set of tail assignment routes r for airplane j. Given r in R j , we define λ r as the coefficient λ P of the s j -t j path P corresponding to r according to Proposition 6. Denoting c r the cost of operating route r with airplane j, Proposition 6 also implies that (10)

a∈A j c a x a = r∈R j λ r c r .
Besides, as x satisfies Equation (TA.4), we have (11) j r∈R j, λ r = 1, where R j, denotes the subset of R j of routes containing leg .

For each airplane j, route r in R j , leg i, and day d, we set x r ij = 1 if leg i is in route r and 0 otherwise, y r jd = 1 if airplane j operating route r undergoes a maintenance on the night after day d and 0 otherwise, and z r ijd = 1 if i is in route r and airplane j operating r undergoes a maintenance on the night after day d and 0 otherwise. As r satisfies the maintenance requirement, (x r ij , y r jd , z r ijd ) id satisfies the equations of the linear relaxation of (KTA) (restricted to airplane j) except cover constraint (3) of [25]. Furthermore, as x ij is the cost of operating leg i with j, we have i∈F c ij x r ij = c r , where c r is the cost of operating route r with airplane j.

Let (x, y, z) be defined as follows. For each airplane j, leg i, and day d, let x ij = r∈R j λ r x r ij , let y jd = r∈R j λ r y r jd , and let z ijd = r∈R j λ r z r ijd . As λ r ≥ 0 and r∈R j λ r = 1, the components indexed by j of (x, y, z) is a convex combination of the (x r ij , y r jd , z r ijd ) id . Hence, (x, y, z) satisfies the equations of the linear relaxation of (KTA) (restricted to airplane j) except cover constraint (3) of [25]. Equation (11) ensures that (x, y, z) also satisfies cover constraint (3) of [25], and is therefore a solution of the linear relaxation of (KTA). Finally, Equation (10) ensures that i∈F j∈P c ij x ij = j a∈A j c a x a , which concludes the proof that the optimal value of the linear relaxation of (TA) is not smaller than the one of the linear relaxation of (KTA).

Consider now the example with n a = 3 and ∆ maint = 4 on Figure 5. Dashed horizontal lines correspond to airports and arrows to legs between airports. The weekly schedule is composed of six round trips between airport A and airport B, and six round trips between airport A and airport C. There is an outward leg . . ,7}. The table provides, for each airplane j, the airport k j 0 where the airplane starts, and the number of days δ j 0 since the last maintenance on day 1, and the costs of operating the different legs with the airplanes. In the remaining of the proof, we show that the optimal value of the linear relaxation of (TA) on this instance is equal to 160, and exhibit a solution of the linear relaxation of (KTA) with value 157.5. Let x = (x j ) j be a solution of (TA) on that instance. Let V j ≤d be the set vertices of D j composed of s j and vertices ( x a = 3.

Hence, any arc a = (v, v ) with x a > 0 in δ + (V j ≤d ) is such that v is in

V j d b ∪ V j d c ∪ V j d f and v belongs to V j d+1 a ∪ V j d+1 c ∪ V j d+1 f
. Any arc a with x a > 0 is therefore of the form ( , δ)( , δ ) with ( , ) in ) .

Furthermore, as we discussed in the first part of the proof, x j can be written as the conic combination r∈R j λ r χ P (r) , where χ P (r) is the indicator vector the s j -t j path P (r) corresponding to route r, and coefficients λ r are non-negative. Remark that there is no connection ( d f , d+1 c

) in a route that satisfies the maintenance requirement, as such a route would spend ∆ maint +1 days out of a base. Constraint (TA.4) applied to leg d+1 c then ensures the only arcs a such that x a > 0 are of the form ( , δ)( , δ ) with ( , ) in ), ( d c , d+1 f ), ( d f , d+1 a ) . Hence, for each airplane j the only routes r in R j that can satisfy λ r > 0 are r 1 = 1 a , 1 b , 2 c , 3 f , 4 a , 4 b , 5 c , 6 f , r 2 = 1 c , 2 f , 3 a , 3 b , 4 c , 5 f , 6 a , 6 b , and r 3 = 2 a , 2 b , 3 c , 4 f , 5 a , 5 b , 6 c , 7 f . Operating any of these three routes has cost 32 with airplane 1, 48 with airplane 2, and 80 with airplane 3. Hence, x has cost 160. Remark that this is the cost of the optimal integer solution obtained by assigning r j to airplane j.

On the contrary setting and z ijd = 0 for any i, j, d such that z ijd has still not been defined, provides a feasible solution of the linear relaxation of (KTA) with cost 157.5, which concludes the proof.

Appendix C. Example of the Monoid resource constrained shortest path problem

This appendix details the execution of Algorithm 2 on a simple example. Figure 6 provides an example of instance of the crew pairing pricing subproblem and its Monoid Resource Constrained Shortest Path Problem modeling. Only a subpart of the instance is represented. On Figure 6.a, legs are represented as arrows between airports. On this instance, there are three airports A 1 , A 2 , and A 3 , and legs only between A 1 and A 2 , and A 2 and A 3 . The flying durations between A 1 and A 2 and A 2 and A 3 are respectively 6 and 2 hours. There are two reduced rests: ( 1 , 3 ) and ( 5 , 8 ). The maximum duty flying duration F (t) in a duty is taken equal to F m = 9 hours for all t. (0, 0, 2, 2), -6.2 (0, 0, 1, 6), 6.2 (0, 0, 1, 2), 2.1 (0, 0, 1, 2), 4.9 (0, 0, 1, 6), 6.2 (0, 0, 1, 2), 5.4

(1, 2), -2.1 (0,0,0,0),-5.8

3

(2,6,0,0),0.4

4

(1,4,0,0),2.5 5 (0,0,0,0),1.2

6

(1,2,0,0),1.5 Table 8. Sets of bounds B v used (here, singletons) Figure 6.b provides the corresponding digraph D, as well as the resource of each arc. Note that the component in M ρ of the resources of day connections are in Z + × R + and the resources of night connections are in (Z + × R + ) 2 . Furthermore, reduced rests are the only night connections whose resources have a component in M ρ of the form (0, 0, 2, •), the other night connections having a component of the form (0, 0, 1, •). The reduced costs have been chosen arbitrarily. For simplicity, we assume that (0, 0, 0, 0), 0 is a lower bound on the resource of 7 -d paths and on the resource of 8 -d paths.

Algorithm 2 execution. We now provide an example of Algorithm 2 when the bounds in Table 8 are used. We remind the reader that these bounds are computed before the execution of Algorithm 2 in a preprocessing and taken by this latter in input. See Section 4.3 for more

2. 1 .

 1 Problem formulation. Building the sequences of flight legs required for aircraft routing corresponds to solving the following problem.

Figure 1 .

 1 Figure 1. Two routes. A two-week route in plain line. A single-week route in dashed line

4 :

 4 denote by c low its optimal value; 5: find a pairing p of minimum reduced cost cp ; (pricing subproblem) 6: if (c p < 0) then 7: add p to P ; 8:

4.2. 1 .

 1 The digraph with its origin and destination vertices. Let D = (V, A) be the acyclic digraph defined as follows. The vertex set is V = L ∪ {o, d}, where L is the set of legs of four consecutive days, o is a dummy origin vertex, and d is a dummy destination vertex. The arc set A contains an arc (o, ) for all legs starting in Paris on the first day of the period, an arc ( , d) for all legs ending in Paris, and an arc ( , ) for each connection ( , ). With these definitions, pairings starting on the first day of the period and satisfying rules (a) and (b) are in one-to-one correspondence with o-d paths in D.

where 1 I

 1 denotes the indicator function of a set I. With this definition, the feasibility function ρ encodes the satisfaction of rules (c) and (d).

Figure 4 .

 4 Figure 4. a. A digraph, and b. the corresponding bounds on resources.

  TA.5) The MIP proposed by Khaled et al. for tail assignment with maintenance constraints is given by Equations (3-15) of their paper. They use binary variables x, y, z. (We denote their variables with an overline to distinguish them from our variables.) The binary variables x ij indicate if the leg i is operated by the airplane j. The binary variables y jd indicate if a maintenance of airplane j takes place on day d. The binary variables z ijd indicate if a night maintenance of airplane j takes place in the arrival airport of the leg i on day d. Air France tail assignment version is slightly different from the one considered by Khaled et al. For instance, they have a constraint limiting the cumulated flight time of an airplane between two maintenances (modeled by Equation (13) of their MIP). The following formulation adapts the MIP of Khaled et al. to Air France tail assignment problem. min (x,y,z) i∈F j∈P c ij x ij (KTA.1) s.t. Equations (

Figure 5 .

 5 Figure 5. Example used in the proof of Proposition 7.

  d a from A to B and an outward leg from d c from A to C every day d in {1,. . . ,6}, a return leg d b from B to A every day d in {1,. . . ,6}, and a return leg d f from C to A every day d in {2,.

  d • , δ) with d ≤ d. As the flow on the cut δ + (V j ≤d ) has value 1, we have 3 j=1 a∈δ + (V j ≤d ) x a = 3. Besides, (TA.3) and (TA.4) give that 3 j=1 a∈δ + (V j d b

x d a 1 = x d a 2 = x d b 1 = x d b 2 = z d b 1d = z d b 2d = z 6 b 1

 1 d ∈ {1, . . . , 6}, x d f 3 = 1 4for all d ∈ {2, . . . , 7},y 1d = y 2d = 1 4for all d ∈ {1, . . . , 7}, and y 3d = 1 2 for d in {1, . . . , 7},

Figure 6 .

 6 Figure 6. Example of instance of the monoid shortest path problem v B v 1

  Algorithm 2 Enumeration algorithm for the Monoid Resource Constrained Shortest Path Problem 1: input: sets (B v ) v∈V satisfying (3); (see Section 4.3) 2: initialization: c U B od ← +∞, L ← ∅, and L nd v ← ∅ for each v ∈ V ; 3: add the empty path at the origin o to L and L nd o ; 4: while L is not empty do

	5:	P ← a path of minimum key(P ) in L;
	6:	L ← L\{P };
	7:	v ← last vertex of P ;

8:

  Suppose that D is acyclic. Then Algorithm 2 converges after a finite number of iterations, and, at the end of the algorithm, c U B od is equal to the cost of an optimal solution of the Monoid Resource Constrained Shortest Path Problem if such a solution exists, and to +∞ otherwise.

	18:	end if
	19:	end if
	20:	end for
	21:	end if
	22: end while
	23: return c U B od ;
	Proposition 3.

Table 2 .

 2 Aircraft

	Instance	Legs Airplane Airplanes Crew	Crew	Short
			connect.	connect. pairings ( ) connect.
	A318	669	39,564	18	3,742		130	1,230
	A319	957	45,901	41	3,738		240	996
	A320	918	49,647	45	3,813		280	1,103
	A321	778	29,841	25	3,918		165	1,006
	A318-9	1,766		-	(59)	8,070		350	2,226
	A320-fam 3,398		-	(129)	21,563		690	4,398
		Table 1. Air France industrial instances
			(AR) alone		(AR) within (Int)
	Instance Uncons. CPU	Optim. CPU	Infeas. CPU	Feas. CPU
		time (mm:ss)	time (mm:ss)	time (mm:ss)	time (mm:ss)
	A318		00:17		00:58		00:14	01:35
	A319		00:16		01:05		00:22	00:19
	A320		01:02		03:55		00:35	13:28
	A321		00:16		01:03		00:23	00:19

routing results 6.2. Experimental setting. All the numerical experiments are performed on a server with 128 Gb of RAM and 12 cores at 2.4 GHz. CPLEX 12.1.0 is used to solve all linear and integer programs. The algorithms are not parallelized.

6.3. Aircraft Routing. Table

Table 4 .

 4 Numerical results on integrated problem

	Instance Integ.	CG it.	(CP) CG	(CP) MIP	(AR)	Sho.	Gap	Total time
		steps	total	time	time			time	Con.	(hh:mm:ss)
	A318	6	460 95.53%	2.56% 1.91%	323 0.0002%	01:53:47
	A319	4	343 76.99%	13.27% 9.74%	448 0.0013%	00:20:18
	A320	2	240 24.24%	39.38% 36.38%	436 0.0017%	00:38:36
	A321	2	380 96.60%	2.53% 0.88%	413 0.0074%	00:29:18
	A318-9	2	915 97.66%	1.71% 0.63%	790 0.0008%	06:34:31
	A320-fam				Stopped after one week
		Instance		318	319	320	321 318-9
		Cost reduction 0.08% 0.16% 0.33% 0.31% 0.31%
		CPU time ratio 10.6× 3.0× 5.6× 2.0× 5.1×
	Table

Table 6 .

 6 Medium size artificial crew pairing instances we introduce in Table

	Instance Legs Connections Pairings ( )
	CP50	290	1,006	50
	CP70	408	1,705	70
	CP90	516	2,490	90

Table 7 .

 7 Relative performance of pricing subproblems algorithms

	Instance	Algorithm	κ	RCSP iter	Cut	RCSP time	Pricing	Total time
				av. nb.	Dom.	av (mm:ss.ff)	time	(hh:mm:ss)
	CP50	Label setting		-1.020e+04	-	00:00.56 97.55%	00:04:38
	CP50	Label correcting		-1.308e+04	-	00:01.28 97.38%	00:11:37
	CP50	Algorithm 2		1 2.326e+03 6.89%	00:00.03 75.28%	00:00:23
	CP50	Algorithm 2	10 4.914e+02 4.01%	00:00.02 59.87%	00:00:17
	CP50	Algorithm 2	100 2.033e+02 5.03%	00:00.04 77.06%	00:00:33
	CP70	Label setting		-5.644e+04	-	00:11.49 99.52%	05:07:05
	CP70	Label correcting		-7.730e+04	-	00:17.16 99.56%	07:28:22
	CP70	Algorithm 2		1 9.208e+03 7.69%	00:00.24 90.61%	00:04:41
	CP70	Algorithm 2	10 1.994e+03 4.28%	00:00.04 58.48%	00:01:12
	CP70	Algorithm 2	100 8.007e+02 5.77%	00:00.07 77.43%	00:01:43
	CP90	Label setting		-9.779e+04	-	00:40.71 Stopped after 48h
	CP90	Label correcting		-2.007e+05	-	01:42.87 Stopped after 48h
	CP90	Algorithm 2		1 5.000e+04 9.81%	00:05.98 98.86%	02:56:33
	CP90	Algorithm 2	10 9.966e+03 5.88%	00:00.34 81.86%	00:12:36
	CP90	Algorithm 2	100 4.377e+03 5.60%	00:00.25 77.98%	00:10:28
	A318	Label setting		-1.319e+05	-	00:53.01 Stopped after 48h
	A318	Label correcting		-3.802e+05	-	01:36.04 Stopped after 48h
	A318	Algorithm 2		1 7.161e+04 8.99%	00:08.61 97.87%	05:35:42
	A318	Algorithm 2	10 5.472e+04 6.62%	00:05.97 96.02%	05:06:47
	A318	Algorithm 2	100 2.549e+04 3.72%	00:01.65 86.97%	01:32:50

Acknowledgments

We thank the reviewers for their thorough reading and all their useful comments and suggestions that helped us improve the paper. We are grateful to Air France which partially supported the project. We are especially thankful to Alexandre Boissy, who initiated the project, and to Mathieu Sanchez and Mohand Ait Alamara who helped us with the implementation in Air France softwares. We also thank Sourour Elloumi for pointing out the correct reference for Lemma 2.

Appendix A. Proofs

Proof of Proposition 1. The fact that a feasible solution of the aircraft routing problem induces a collection C as in the statement is obvious. Let us prove the other direction, which is almost as easy.

Consider a collection C as in the statement. Each cycle provides a route, possibly of several weeks. We show now that the solution consisting of these routes is feasible. By construction of the graph, the maintenance requirement is satisfied. Moreover, the number of times a cycle intersects A 0 is an upper bound on the number of airplanes required to operate the corresponding route: the number of times it intersects A 0 is equal to the number of weeks this cycle lasts, as illustrated on Figure 3. Thus the number of arcs selected in A 0 by the whole collection is an upper bound on the number of airplanes required to operate the solution. Since this number is at most n a , the solution is feasible.

Proof of Lemma 2. Consider the following integer program in the standard form:

By the theory of the simplex algorithm, the program (7) can be written in the form

where the ci -the reduced costs -are all nonnegative, and where N is the complement of the considered basis (all indices i such that ci > 0 are in N ). Consider now an optimal solution z * of the program (7). We have

-4.8 (0, 0, 1, 6) ( 1 , 3 ) (0, 0, 2, 2) (0, 0, 2, 2) 3 (2, 6, 0, 0) (0, 0, 0, 0) yes ( 1 , 5 ) (0, 0, 1, 6) (0, 0, 1, 6) 5 (0, 0, 0, 0) (0, 0, 0, 0) yes ( 1 , 6 ) (0, 0, 1, 2) (0, 0, 1, 2) 6 (1, 2, 0, 0) (0, 0, 0, 0) yes Let P 1 be an o-1 path with resource (0, 0, 1, 6), 1.0 . Table 9 describes the iterations of Algorithm 2 where P 1 and the paths starting by P 1 are dealt with. Each iteration is separated by an horizontal line. Column P provides the path P considered at Step 5 of Algorithm 2. We assume that L nd 3 , L nd 5 , and L nd 6 are empty when P 1 is considered as path P , and L nd 4 is empty when P 1 , 3 is considered. We also assume that c U B od = +∞ and hence c(q

od during all the iterations detailed. Column key(P ) provides its key defined in Equation [START_REF] Barnhart | Flight string models for aircraft fleeting and routing[END_REF]. As the treatment of reduced costs is standard, to enhance readability, we omit reduced cost in all resources in Table 9 and in the remaining of the discussion. Column q P gives the resource of P , column a provides the arc of Step 11, column q a gives its resource. Path Q of Step 12 is path P followed by a. Then next column gives the resource q Q = q P ⊕q a of Q, and column w provides the destination of Q computed at Step 13. Column b provides the single bound in B w , and the next column provides q Q ⊕ b computed at Step 14. Finally, the last column indicates if Q is added to L at Step 17.

The key of P 1 is equal to -4.8 because, with P = P 1 and B 1 = {b}, we have q P ⊕ b = (0, 0, 1, 6), 1.0 ⊕ (0, 0, 0, 0), -5.8 ) = (0, 0, 0, 0), -4.8 , and c (0, 0, 0, 0), -4.8 = -4.8. The sums (0, 0, 2, 2) ⊕ (2, 6, 0, 0) and (0, 0, 3, 4) ⊕ (1, 4, 0, 0) are equal to (0, 0, 0, 0) because (2, 2) + (2, 6) = (3, 4) + (1, 4) = (4, 8) ≤ (4, F m ), where F m = 9. Since ρ((0, 0, 0, 0)) = 0, the path Q cannot be discarded at the iterations where P = P 1 and a = ( 1 , 3 ), and where P = P 1 , 3 and a = ( 3 , 4 ), and it is added to L. On the contrary (0, 0, 4, 10)⊕(0, 0, 0, 0) = ∞ because (4, 10) (4, F m ). Since ρ(∞) = 1, the path Q is not kept after Step 14 when a = ( 4 , 5 ), and it is not added to L. We have a similar outcome when a = ( 4 , 6 ): in this case, (0, 0, 4, 6) ⊕ (1, 2, 0, 0) = ∞ because (5, 8) (4, F m ). The treatment of Q = P 1 , 3 , 4 , 6 shows the interest of the bounds: although path Q itself satisfies rule (c), the algorithm identifies that any path starting by Q violates rule (c).

Table 10 provides the same informations as Table 9 when we use a set of bounds B 3 = (2, 8, 0, 0), 0.4 , (3, 6, 0, 0), 1.5 instead of the singleton given in Table 8. Fewer iterations are then needed: the "if" condition at Step 14 is not satisfied when P 1 , 3 is considered as path Q, and path P 1 , 3 is never added to L. Even though P 1 , 3 itself satisfies rules (c) P key(P )

( 1 , 5 ) (0, 0, 1, 6) (0, 0, 1, 6) 5 (0, 0, 0, 0) (0, 0, 0, 0) yes ( 1 , 6 ) (0, 0, 1, 2) (0, 0, 1, 2) 6 (1, 2, 0, 0) (0, 0, 0, 0) yes P 1 , 6 4.6 (0, 0, 1, 2) ( 6 , 7 ) (1, 2) (0, 0, 2, 4) 7 etc. P 1 , 5 8.4 (0, 0, 1, 6) ( 5 , 8 ) (0, 0, 2, 6) (0, 0, 2, 6) 8

Table 10. Algorithm execution with B 3 = (2, 8, 0, 0), 0.4 , (3, 6, 0, 0), 1.5 and (d), the algorithm identifies that any o-d path starting by P 1 , 3 does not satisfy at least one of these rules.

Rationale of B 3 . We explain why B 3 is a correct bound set, both in the singleton and non-singleton cases. This explanation can be seen as a rough illustration of the procedure mentioned in Section 4.3 for the bound computation. Any 3 -d path must either start with 3 , 4 , 5 , 8 , or with 3 , 4 , 6 , 7 . Recall that we have assumed that (0, 0, 0, 0), 0 is a lower bound on the resource of 7 -d paths and on the resource of 8 -d paths. Given that q ( 3 , 4 ) ⊕ q ( 4 , 5 ) ⊕ q ( 4 , 8 ) ⊕ (0, 0, 0, 0), 0 = (2, 8, 0, 0), 0.4 and q ( 3 , 4 ) ⊕ q ( 4 , 6 ) ⊕ q ( 6 , 7 ) ⊕ (0, 0, 0, 0), 0 = (3, 6, 0, 0), 1.5 , any path 3 -d path starting by 3 , 4 , 5 has a resource lower bounded by (2, 8, 0, 0), 0.4 , and any path 3 -d path starting by 3 , 4 , 6 , 7 has a resource lower bounded by (3, 6, 0, 0), 1.5 . This explains why (2, 8, 0, 0), 0.4 , (3, 6, 0, 0), 1.5 can be used as set of bounds B 3 , and why (2, 8, 0, 0), 0.4 ∧ (3, 6, 0, 0), 1.5 = (2, 6, 0, 0), 0.4 is a lower bound on the resource of any 3 -d path.

A. Parmentier, Université Paris Est, CERMICS, 77455 Marne-la-Vallée CEDEX, France Email address: axel.parmentier@enpc.fr F. Meunier, Université Paris Est, CERMICS, 77455 Marne-la-Vallée CEDEX, France Email address: frederic.meunier@enpc.fr