Nicolas Barnier
email: barnier@recherche.enac.fr

Pascal Brisset
email: brisset@recherche.enac.fr

Solving the Kirkman's Schoolgirl Problem in a Few Seconds

Keywords: Symmetry Breaking, Social Golfer Problem, Resolvable Steiner Systems

The Social Golfer Problem has been extensively used in recent years by the constraint community as an example of highly symmetric problem. It is an excellent problem for benchmarking symmetry breaking mechanisms such as SBDS or SBDD and for demonstrating the importance of the choice of the right model for one problem. We address in this paper a specific instance of the Golfer Problem well known as the Kirkman's Schoolgirl Problem and list a collection of techniques and tricks to find efficiently all its unique solutions. In particular, we propose SBDD+, an generic improvement over SBDD which allows a deep pruning when a symmetry is detected during the search. Our implementation of the presented techniques allows us to improve previous published results by an order of magnitude for CPU time as well as number of backtracks, and to compute the seven unique solutions of the Kirkman's problem in a few seconds.

Introduction

Highly symmetric problems are always challenging for Constraint Programming and breaking, removing, discarding symmetries among solutions has been the subject of much interest among researchers of the CP community in recent years. We focus in this paper on one particular symmetric problem: the Social Golfer Problem [START_REF] Gent | CSPlib: a problem library for constraints[END_REF], also known as resolvable Steiner system in the combinatorial area [START_REF] Jr | Combinatorial Theory[END_REF]. Except for small instances, this problem is open and Constraint Programming gives balanced results: formulation is straightforward but far from being sufficiently efficient to solve all instances.

We are interested in breaking as much symmetry as possible, as well as combining and improving previously proposed techniques to find all solutions of one specific instance of the problem. The first important choice concerns the model. We naturally choose a set model [START_REF] Gervet | Interval propagation to reason about sets: Definition and implementation of a practical language[END_REF] which automatically removes one kind of symmetry. The next step is to statically remove symmetries by adding constraints. Additional redundant constraints may be added to detect failures as soon as possible. The crucial point is then to be able to find an isomorphism relating two solutions quickly; we propose a "lazy" approach which splits the detection into two phases, building and checking, that provides the required efficiency.

One of the key ideas of the paper is to exploit an isomorphism found between two solutions, or partial solutions, to prune as much as possible of the subsequent search tree. We show that this is possible if the structure of the search tree is intrinsically related to the symmetries; in this case an isomorphism which maps a solution to another one similarly maps also some ancestor nodes. Combined with SBDD [START_REF] Fahle | Symmetry breaking[END_REF], we call the technique SBDD+.

Our experiments show that results presented in previous papers [START_REF] Sellmann | Heuristic constraint propagation[END_REF] can be greatly improved both in number of failures and in CPU time. The problem of finding all solutions to the Kirkman's Problem, which might have been considered hard to solve one year ago using Constraint Programming (two hours of CPU mentioned in [START_REF] Fahle | Symmetry breaking[END_REF]) can be solved in few seconds using our approach.

The remainder of the paper is structured as follows: we first define the Social Golfer Problem and model it, giving numerous possible redundant constraints. In the next section, we present our algorithm used to check symmetry and explain extensively our deep pruning technique associated with symmetry finding. Section 4 displays results of our experiments, confirming the approach. We conclude by recalling that many challenges remain.

Model

The classical instance of the Social Golfer Problem is described in the following terms: 32 golfers want to play in 8 groups of 4 each week, in such way that any two golfers play in the same group at most once. How many weeks can they do this for?

The problem may be generalized to w weeks of g groups, each one containing s golfers. This instance will be denoted g-s-w in the sequel of the paper. We note n = g × s the total number of golfers. The most famous and historical instance is the 5-3-7 for which all 7 unique (non-symmetric) solutions were already computed by Kirkman in the early 1850's [START_REF] Kirkman | Note on an unanswered prize question[END_REF]. In the combinatorics area, solutions for s = 3 are known as Kirkman triple systems or resolvable Steiner systems. Such systems have been extensively investigated (see for example [START_REF] Cohen | Kirkman triple systems of order 21 with nontrivial automorphism group[END_REF]).

In the context of constraint programming, different models have been proposed in [START_REF] Smith | Reducing symmetry in a combinatorial design problem[END_REF] to the Golfer Problem. The integer set model1 which automatically removes symmetries inside groups is the one we chose for our experiments. In this model, the variables are the groups themselves and constraints are expressed as operations on sets.

The variables G i,j , with i index of weeks and j index of groups, are sets and their associated domain is a lattice of sets defined by its greatest lower bound (the necessary elements) and its lowest upper bound (the possible elements) [START_REF] Gervet | Interval propagation to reason about sets: Definition and implementation of a practical language[END_REF]. The G i,j 's are subsets of the set of golfers. Each of them contains exactly s elements. All the groups of a week are disjoint and every pair of groups from different weeks share at most one element. All these properties are expressed with the following constraints:

1 ≤ i ≤ w, 1 ≤ j ≤ g G i,j ⊂ {1, 2, ..., n} 1 ≤ i ≤ w, 1 ≤ j ≤ g |G i,j | = s (1) 1 ≤ i ≤ w, 1 ≤ j < j ′ ≤ g G i,j ∩ G i ′ ,j = ∅ (2) 1 ≤ i < i ′ ≤ w, 1 ≤ j, j ′ ≤ g |G i,j ∩ G i ′ ,j ′ | ≤ 1 (3)
The constraints (2) may be basically implemented as w all_disjoint global constraints instead of the wg(g -1)/2 binary disjoint constraints. Note that no global consistency is achieved for this global constraint. The number of constraints of type (3) grows quadratically with the number of groups. It may prevent to solve large instances with this model.

According to [START_REF] Smith | Reducing symmetry in a combinatorial design problem[END_REF] experiments, the naive set model is not the best one. However, we choose it for several reasons. First it is the simplest one and it uses the highest abstraction level. Second, redundant constraints described in the following section are easy to express with this model. Third, this is the model used in [START_REF] Fahle | Symmetry breaking[END_REF][START_REF] Sellmann | Heuristic constraint propagation[END_REF] and it allows us to compare our approach with this previous ones.

Redundant constraints

Several constraints can be added to the original model. While they may help to solve the hardest instances, the induced overhead is sometimes too large for small instances like Kirkman's one.

The fact that a player plays only once per week is not explicit in the original model, but only entailed by the constraints (1) and [START_REF] Cohen | Kirkman triple systems of order 21 with nontrivial automorphism group[END_REF]. The corresponding constraint is written using reified membership constraints:

1 ≤ i ≤ w, 1 ≤ p ≤ n 1≤j≤g (p ∈ G i,j ) = 1 (4)
Warwick Harvey (www.icparc.ic.ac.uk/~wh/golf) proposes to express the fact that the players of a group appear in exactly s groups in other weeks:

1 ≤ i = i ′ ≤ w, 1 ≤ j ≤ g 1≤j ′ ≤g (G i,j ∩ G i ′ ,j ′ = ∅) = s (5)
Taking into account the size of the groups, the global constraint atmost1 proposed by [START_REF] Sadler | Global reasoning on sets[END_REF] may also be set on the list of all groups.

atmost1({G i,j /1 ≤ i ≤ w, 1 ≤ j ≤ g}, s) (6)
where atmost1(S, c) states that sets of S must have cardinal c and must intersect pairwise in atmost one element. The propagation associated with this constraint basically ensures that the possible number of partners of a player p is large enough, i.e. greater or equal to (c -1)N p where N p is the number of occurrences of p. In our case, N p is statically known (equal to w) so the propagation rule of the constraint can be efficiently customized.

Breaking Symmetries Statically

Our first goal is to compute all the unique non-symmetric solutions to the problem. As described in previous papers [START_REF] Smith | Reducing symmetry in a combinatorial design problem[END_REF][START_REF] Fahle | Symmetry breaking[END_REF][START_REF] Sellmann | Heuristic constraint propagation[END_REF], the Social Golfer Problem is highly symmetric:

-Players can be exchanged inside groups (φ P ); -Groups can be exchanged inside weeks (φ G ); -Weeks can be ordered arbitrarily (φ W ); -Players can be renamed among n! permutations (φ X ).

The symmetry inside groups is inherently removed by modelling groups as sets.

The symmetry inside weeks may be handled by ordering the g groups. Because these groups are disjoint, a total order can be achieved by sorting the smallest element of the groups.

1 ≤ i ≤ w, 1 ≤ j ≤ g -1 min G i,j < min G i,j+1
Note that this implies that the first player is in the first group for each week. Following the same idea, weeks can be ordered with the first group as key, which can be easily done with the second smallest element:

1 ≤ i ≤ w -1 min(G i,1 \ {1}) < min(G i+1,1 \ {1})
Symmetries among players are more difficult to handle and only dynamic checks will be able to remove them completely. Statically:

-First week is fixed; -First group of second week is fixed with smallest possible players -"Small" players are put in "small" groups: for every week, p th player is in a smaller group than the p th group

1 ≤ i ≤ w, 1 ≤ p ≤ g Ḡi,p ≤ p
where Ḡi,p is the number of the group of player p in week i, i.e. the Ḡi,p are dual variables defined by

1 ≤ i ≤ w, 1 ≤ p ≤ n Ḡi,p = j iff p ∈ G i,j
-Players together in a same group in the first week are placed in ordered groups in the second week

1 ≤ j ≤ g, p 1 , p 2 ∈ G 1,j , p 1 < p 2 Ḡ2,p1 < Ḡ2,p2
-Groups of the first week are ordered in the second week:

1 ≤ j < j ′ ≤ g Ḡ2,G1,j < lexico Ḡ2,G 1,j ′
where Ḡi,{x1,x2,...} is the tuple ( Ḡi,x1 , Ḡi,x2 , ...) and < lexico stands for the lexicographic order on integer tuples.

Unfortunately, the conjunction of all these constraints does not remove all the symmetries among players. For example, for the 5-2-2 instance, the two following solutions are found: 1 2 3 4 5 6 7 8 9 10 1 3 2 4 5 7 6 9 8 10 1 2 3 4 5 6 7 8 9 10 1 3 2 5 4 6 7 9 8 10 Both solutions satisfy all aforementioned breaking symmetry constraints but the second one is isomorphic to the first one through the functions (φ G is the same for the two rounds):

φ X = {1 → 7, 2 → 8, 3 → 9, 4 → 10, 5 → 1, 6 → 2, 7 → 3, 8 → 4, 9 → 5, 10 → 6} φ G = {1 → 4, 2 → 5, 3 → 1, 4 → 2, 5 → 3} φ W = {1 → 1, 2 → 2}
We notice that even if permutations within weeks and groups may be statically removed by constraints when considered alone, it is still necessary to take them into account when the permutation within players is handled.

Remaining symmetries must be dynamically discarded; we discuss an efficient way to do it in the next section.

Integer Model with a Cardinality Constraint

We give here another model which allows us to solve efficiently the problem. Its originality comes from the use of a global cardinality constraint [START_REF] Régin | Generalized arc consistency for global cardinality constraint[END_REF]. To the best of our knowledge, it is the first time it is proposed.

Decision variables in this model are Ḡi,p , the number of the group of player p in week i.

1 ≤ p ≤ n Ḡi,p ∈ [1..g] 1 ≤ i ≤ w gcc({ Ḡi,p /1 ≤ p ≤ n}, < (1, s), ..., (g, s) >)
where gcc(S, < (v 1 , c 1 ), ..., (v k , c k ), ... >) constrains the number of occurrences of elements of S equal to v k to be equal to c k .

The "not with the same golfer more than once" is straighforward with reified constraints (whose number is increasing quadratically with the number of golfers):

1 ≤ p 1 < p 2 ≤ n 1≤i≤w ( Ḡi,p1 = Ḡi,p2 ) ≤ 1
This model allows for example a short program2 to solve the classic 8-4-9 instance in half a second and 32 backtracks to find the fist solution. However, as explained earlier, it is not the one used for the experiments in this paper.

Handling Symmetries During Search

In this section we present our adaptation to the Social Golfer Problem of a generic symmetry breaking mechanism proposed in [START_REF] Fahle | Symmetry breaking[END_REF].

Generic Techniques for Breaking Symmetries

Symmetry breaking constraints fail to remove statically all symmetries among the players in the Social Golfer problem. Therefore, several solutions have been proposed to prune the search tree taking into account these symmetries dynamically.

Using SBDS [START_REF] Gent | Symmetry breaking during search in contraint programming[END_REF], which needs to list explicitly the symmetries to remove, Barbara Smith in [START_REF] Smith | Reducing symmetry in a combinatorial design problem[END_REF] was able to break most of the symmetries and obtain new results. Later, two generic and similar approaches were proposed in the same time [START_REF] Focacci | Global cut framework for removing symmetries[END_REF][START_REF] Fahle | Symmetry breaking[END_REF]. In the second one, the technique called SBDD (for Symmetry Breaking via Dominance Detection) was applied with success to the Social Golfer Problem and allows us to compute all the non-symmetric solutions of small instances.

In SBDD, states during the search (i.e. nodes of the search tree) are compared to previously explored ones modulo a symmetry mapping function. A (new) state P ′ is dominated by an (old) state P if P ′ is subsumed by φ(P ) where φ is a symmetry mapping function. When searching only for non-symmetric solutions, a state which is dominated by an already explored ones is discarded. Then, it requires to store all explored nodes. However, it can be noticed that if P ′ is dominated by P then it is dominated by the father node of P . It means that when all the sons of a state have been explored, one can remove them from the store and keep only the father. Concretely, in case of depth first search, the store for SBDD can be handled as a stack where states are associated to their depth. When a new state must be added to the store, all the states on top of the stack which have a greater depth may be removed. We will see later that it is worthwhile to store states in a compiled form in order to ease future dominance checks against this state.

One issue of the technique is its efficiency because checking dominance may be very expensive: w(s!) g g! symmetries to check in the Social Golfer Problem. Some restrictions are necessary to get an effective procedure:

storage of explored states may be limited; -checking of dominance may be restricted to some depths.

We propose to specialize the SBDD technique for the Social Golfer Problem, first to be able to check for dominance quickly, second, to better exploit symmetries found.

Filtering Efficiently Symmetric Solutions

It is shown in experiments of [START_REF] Fahle | Symmetry breaking[END_REF] that it is not worth to check dominance for every node during search for golfer solutions. For the 4-4-4 instance, authors conclude that checks every 8-th depth give the best result. Results given in the next section show that lazy dominance checking is effective when solving small instances of the Golfer Problem.

General dominance check for the Golfer Problem as described in [START_REF] Fahle | Symmetry breaking[END_REF] requires to compute a matching in a bipartite-graph, for which the best algorithm is in O(n 5/2 ) [START_REF] Hopcroft | An n 5/2 algorithm for maximum matching in bipartite graphs[END_REF]. However checking that there exists a symmetric mapping function which maps an old solution (a leaf in the search tree) to a new one is significantly easier.

Actually, it can be noticed that a solution to the Golfer Problem is fully specified by the week number of the pairs of players 3 . Precisely, a solution can be described by the following mapping:

(p 1 , p 2 ) → i such that Ḡi,p1 = Ḡi,p2 (7)
where (p 1 , p 2 ) is a pair of golfers. Note however that the mapping is not total for instances where a player does not play with all others.

A check must be done for each possible symmetry. [START_REF] Fahle | Symmetry breaking[END_REF] remarks that the possible symmetries may be easily enumerated looking for a matching from the first week of the first solution (or partial solution) P to any week of the second solution P ′ . The key idea to compute symmetry checking efficiently is to compute it lazily: instead of choosing a complete matching and checking the other weeks afterward, it is worth checking them while the matching is built. Figure 1 illustrates this principle: Suppose a partial isomorphism is built by mapping an exchange of the first two groups of the first week (1 ↔ 4, 2 ↔ 5, 3 ↔ 6). There is enough information in this partial matching to check all pairs among other weeks with index of players less than 6. In this example, one can check that pairs (1, 4) and (2, 5) of the second week would be mapped in the same week of the second solution. Because it is not the case (pair (1, 4) is mapped to (4, 1) which appears in second week, (2, 5) is mapped to (5, 2) which does not appear at all), it is not necessary to try to complete this partial matching and we can "backtrack" to consider another mapping for the first two groups of the first week.

The check for pair mapping may be easily performed at a low cost by a precomputation phase:

-Derive the set of pairs of each week i of the first solution P , sorted according to the greatest element of the pair, denoted C i (P ). Note that this precomputation may be done only once when the solution is found and stored. -Build a table of the week numbers indexed by the pairs of the second solution P ′ (mapping of equation ( 7), noted W P (c)).

The complete procedure is detailed in figure 2. The worst case complexity is the same as a naïve approach but experiments show that this checking algorithm, while only applied on leaves of the search tree, is efficient enough to compute all solutions of the Kirkman's problem in a reasonable time (c.f. section 4).

Pruning Deeply

An efficient symmetry checking procedure applied on leaves allows us to compute all the unique solutions but does not improve the search itself: no subtrees are removed. However, following an idea of [START_REF] Mckay | Practical graph isomorphism[END_REF] used in an algorithm to compute graph isomorphisms 4 , symmetry checking on leaves may be used to prune large subtrees. The idea is illustrated in figure 3. Let P be a first solution (the tree is explored in depth first search, from left to right) and P ′ a second solution proved to be isomorphic to P (P ′ = γ(P )). We note n = P \ P ′ the lowest node in the search tree which is common to paths from the root to P and P ′ , and s (resp. s ′ ) its immediate successor (we suppose that we do only binary branching) leading to P (resp. its immediate successor leading to P ′ ). Under some conditions (which we call "McKay condition" in the sequel), it can be shown that node s ′ is the image of s by the isomorphism γ (more precisely the canonical extension of γ over partial solutions). In this case, the remaining unexplored subtree of s ′ is itself the image of an already explored subtree starting from s. Then it can be pruned because it would lead only to solutions that are images by γ of already found solutions.

In order to be able to apply this idea to a search tree, the structure of the tree (i.e. the labelling procedure) must be compatible with the symmetric mapping function we consider. For the Golfer Problem, all symmetries but the one on players permutation may be removed simply with constraints. So the symmetries discovered at the leaves concern only φ X . Following these remarks, the right choice is to label golfer by golfer to be able to apply φ X extension on a node of the search tree. Note that a full choice for one golfer p amounts to labelling the w Ḡi,p variables.

Unfortunately, this labelling does not ensure the McKay condition if the set of golfers above node n (c.f. figure 4) is not stable through the found isomorphism γ: suppose node n concerns golfer 3 and γ is such that γ(1) = 4; we clearly cannot have s ′ = γ(s) in this case. The problem in this example is that the choice on golfer 1 for P is mapped to the choice on golfer 4 in P ′ , the latter being under the node n. Hence, node n cannot be pruned.

The necessary adaptation of McKay idea for our search tree is to consider two nodes s, s ′ which are descendants of n, leading respectively to P and P ′ , such that the set of choices above s is mapped by γ to the set of choices above s ′ . In this case s ′ can be pruned. Of course, it is better to choose the highest such node in the search tree to prune as much as possible. Figure 4 illustrates the idea; the smallest set of golfers stable for γ which contains 3 (node n) is {1, 2, 3, 4, 5, 6}.

Results given in next section show that this deep pruning is highly effective.

SBDD+

We propose in this section an integration of the idea of deep pruning presented in the last section with the symmetry breaking mechanism SBDD. This mechanism computes isomorphism not only on solutions (leaves of the search tree) but on states of the search tree which can be described with the current domains of the variables. In SBDD+, we will exploit such isomorphism and try to show that it is applicable on ancestors, leading then to a best pruning. Based on SBDD, the approach is generic: no hypothesis are required on the nature of the variables (integer or set variables) or the nature of the isomorphism. However, we present here the application of the method only to the Golfer Problem and we did not try to implement it in a generic way.

Fine symmetry checking on leaves improved with deep pruning allows us to solve our problem in a reasonable time. However, better pruning may be obtained if the method is integrated with SBDD approach. It requires first to be able to apply a refinement of the algorithm given in figure 2 to incomplete solutions (nodes of the search tree), second to call the procedure only at the appropriate times.

Experiments show that most of the symmetries found between two complete solutions involve a mapping from the first week to itself 5 . Moreover, incomplete solutions always get a complete first week. This means that the previous symmetry checking algorithm may be easily modified to be applied to incomplete solution if only this kind of symmetry is considered; building of the symmetry can be kept. However, the checking phase of the algorithm remains the same: a symmetry is found if all pairs of P get an image in P ′ . But it does not mean that the symmetry is necessarily found when it exists, i.e. we possibly do not find all symmetries.

Dominance checking remains expensive and it must not be done too often. The check frequency must of course be related to the structure of the problem. A good compromise for the Golfer Problem is to -Store nodes at every depth of the search tree; -Check dominance for nodes only against stored nodes of smaller depth; -Check dominance only for nodes at depth multiple of s.

The maximum size of the node store may be estimated: the depth of the search tree is the number of golfer gs; at each level, there may be g w (for one golfer, g choices for each week) nodes to store. So the size of the store is bound by gsg w (12 890 625 for Kirkman's Problem). This bound is a bad upper bound due to numerous symmetries removed by constraints (first week fixed, small golfers fixed, ...). For the Kirkman's Problem, with the detailed previous choices, only [START_REF] Régin | Generalized arc consistency for global cardinality constraint[END_REF] nodes are effectively stored in the same time.

SBDD is compatible with the deep pruning mechanism: when a dominance is found, it is usually possible to prune more than the dominated node, just looking for the highest ancestor of the node for which the McKay condition is verified. We call the method SBDD+. It is illustrated in figure 5.

Results

We give in this section results of our experiments with the different techniques described in previous sections. The implementation has been done using FaCiLe [START_REF] Barnier | Facile: a functional constraint library[END_REF], a functional constraint library written in Objective Caml (caml.org). This library provides usual finite domain variables over integers and sets, arithmetic constraints, set constraints [START_REF] Gervet | Interval propagation to reason about sets: Definition and implementation of a practical language[END_REF], global constraints and an extensive mechanism to build search goals. It also includes basic mechanisms to maintain invariant expressions. This library, under LGPL license, and its documentation are available at www.recherche.enac.fr/opti/facile. CPU times are given for a Pentium 700 MHz running Linux 2.4.5.

The set model has been used with the labelling described earlier, golfer by golfer, choosing a group for each week, in increasing order. Additionally, some simple symmetries are dynamically removed during the search: when a golfer is placed in a group which does not contain any other golfers, this choice is fixed 6 . The refinement of the redundant atmost1 (c.f. equation ( 6)) constraint is set in every experiment.

Labeling is slightly improved according to the following remark: the unique solutions of 5-3-7 are extensions of unique solutions of 5-3-2. Let S and S ′ be isomorphic solutions of 5-3-2 such that S ′ = γ(S). If S is extended into a solution P of 5-3-7, then it is clear that γ(P ) is an isomorphic solution of 5-3-7 and also is an extension of 5-3-2. Then, our labeling first computes all unique solutions of 5-3-2 and extend them to 7 weeks. There are only 2 unique solutions for 5-3-2 and it takes 0.1 seconds to compute them.

Table 1 shows the number of created choice points, the number of backtracks, the total number of found solutions, the number of dominance checks and the CPU time for different combinations of methods and tricks presented in this paper to compute all 7 unique solutions to the schoolgirl problem (then 11 found solutions means that isomorphism has been detected only at the leaf of search the tree for 4 of them). First column (Leaves) gives the results for the straightforward search with simple discarding of symmetric solutions at leaves. This time can be compared with the one annouced in [START_REF] Fahle | Symmetry breaking[END_REF] (two hours). Next column (McKay) corresponds to the symmetry detection at leaves with deep pruning. We see that the number of failures and the time are smaller by an order of magnitude from the previous ones. Column SBDD uses our incomplete dominance checking; SBDD+ adds deep pruning. It is the best time achieved with our experiments. It can be compared with the result of [START_REF] Sellmann | Heuristic constraint propagation[END_REF] (400s on a similar CPU).

However, the number of backtracks can still be reduced with redundant constraints: in column "+(4)", the redundant constraint (4) stating that a golfer plays only once per week allows us to further reduce the search but the overhead is too large (with our implementation) and CPU time is not improved. The last column adds redundant constraint (5) which expresses that players of a group are spread among s different groups in other weeks. The overhead is here dissuasive but the number of backtracks is 5 times smaller than what was done in [START_REF] Sellmann | Heuristic constraint propagation[END_REF]. Our combination of tricks to solve the Golfer Problem allowed us to solve open (at least for constraint programming) instances (6-4-6, 7-3-9, 8-3-7, 7-4-6, 6-5-7, ...). Some of these instances, at the time of writing this paper, are no longer open [START_REF] Prestwich | Randomised backtracking for linear pseudo-boolean constraint problems[END_REF] and last updates of Warwick Harvey's web page include all these results (http://www.icparc.ic.ac.uk/~wh/golf/).

We have presented in this paper a combination of techniques which allows us to find efficiently all solutions to the Golfer Problem. The main contribution of the paper is an improvement of the SBDD principle that we call SBDD+. The key idea of the improvement is, while breaking symmetries, to exploit the symmetry function to be able to prune higher in the search tree. Extensive experiments show that this new mechanism can reduce by an order of magnitude the CPU time as well as the number of backtracks on the considered problem.

The deep pruning technique has been applied only to the Golfer Problem but is general. The only restriction is to have a relative compatibility between the structure of the search tree and the considered symmetry mappings in order to be able to prove that a symmetry found between two nodes is also true for ancestors of these two nodes. We believe that the notion of stability through the isomorphism of ancestor nodes, a necessary and sufficient condition for our problem, should be a general property. Further work is needed in this direction.

For the Golfer Problem itself, some instances remain open to constraint programming approaches (even if they are well known by combinatorics, for example 7-3-10 and 7-4-9 are extensively studied). Our model may be improved using the incomplete propagation techniques proposed by [START_REF] Sellmann | Heuristic constraint propagation[END_REF] in order to attack these instances. Note that SBDD+ must be refined and tuned for larger instances to avoid an explosion of node store and an extreme time overhead due to node dominance checking.

In spite of many efforts from the constraint community, the 8-4-10 instance is still open. This challenge is fascinating and it can be considered with the highest priority to show that constraint technology is really suited for this kind of combinatorial problem.

mapped to week 2 FailureFig. 1 .

 21 Fig. 1. Partial symmetry between solutions. Solid arrows show isomorphism building while dashed lines show isomorphism checking.

Fig. 3 .

 3 Fig. 3. Illustration of McKay's idea: The dark subtree can be pruned.

5 PFig. 4 .

 54 Fig. 4. Deep Pruning for the Golfer Problem.

PFig. 5 .

 5 Fig. 5. SBDD+: Deep pruning integrated to SBDD.

  Search for symmetry φX mapping a solution P to P ′

	continue	mapped to the same week in P ′
	return φX	Symmetry found is returned
	return NoSymmetry
	Fig. 2.

Procedure CheckSymmetry(P, P') Compute Ci(P ) for i ∈ 2..w Usually already computed with P Compute W P ′ (c) for all pairs c of P ′ for i ′ ∈ 1..w Map first week of P to week i ′ of P ′ φW [1] ← i for φG ∈ Pg Permute groups within the week for j ∈ 1..g for φP ∈ Ps Permute players within the group Set φX such that P ′ i ′ ,φ G (j) = φX (φP (P1,j)) try for i ∈ 2..w For all other weeks of P for c ∈ Ci(P ) s.t. c ≤ j s For all mapable pairs if c is the first encountered pair of Ci(P ) φW [i] ← W P ′ (φX (c)) Store image week of week i else if W P ′ (φX (c)) <> φW [i] Check if pairs of week i in P are

Table 1 .

 1 Computing the 7 solutions of the Kirkman's Problem

		Leaves McKay SBDD SBDD+ +(4) +(5)
	Choice points	20062206 1845543 107567 29954 18705 18470
	Fails	19491448 1803492 104134 28777 16370 16169
	Solutions	20640	934	11	11	11	11
	Dominance checks			5373	456	456	443
	CPU(s)	5925	484	24	7.8	9.4	36

http://www.icparc.ic.ac.uk/eclipse/examples

This solution for the Golfer Problem is an example provided in the constraint library we use.

Is is the integer model proposed by Barbara Smith in[START_REF] Smith | Reducing symmetry in a combinatorial design problem[END_REF].

The nauty software based on McKay ideas is used by combinatorics people to find resolvable Steiner sytems[START_REF] Cohen | Kirkman triple systems of order 21 with nontrivial automorphism group[END_REF].

For 7 non-symmetric solutions, 927 solutions are found and

of them have this property.

A similar idea is used in graph coloring algorithm: if a "new" unused color is tried for a node, there is no need to consider later other new colors for this node.


Géraud Granger
email: granger@recherche.enac.fr

Steria Eurocontrol

Nicolas Durand
email: durand@recherche.enac.fr

A traffic complexity approach through cluster analysis

Keywords: cluster, graphic sequences, traffic complexity, conflict resolution

The conflict resolution problem is quite simple as long as the number of aircraft involved is small. Many automation projects disregarded the problem of clusters and failed on real traffic tests because they were unable to deal with complex conflicts. The first definition of cluster appeared in the middle of the nineties when theoretical research started on conflict resolution. The n-aircraft conflict resolution problem is highly combinational and cannot be optimally solved using classical mathematical optimisation techniques. The set of admissible solutions is made of many unconnected subsets enclosing different local optima, but the subset enclosing the optimum cannot be found a priori. Using a priority order to solve a n-aircraft conflict is much easier but the solution is not optimal. However it is difficult to determine the best order or even a good order that ensures that a solution exists. In this paper, a theoretical study of the possible structures of clusters is presented. A simulation using French real traffic data compares the structure of clusters with direct and standard routes. The sensitivity of cluster sizes to uncertainties on trajectories forecast is studied.

Introduction

The conflict resolution problem is quite simple as long as the number of aircraft involved is small. Many automation projects disregarded the problem of clusters and failed on real traffic tests because they were unable to deal with complex conflicts. The American AERA project [START_REF] Niedringhaus | A mathematical formulation for planning automated aircraft separation for AERA3[END_REF] did not suggest any algorithm for the MOM 1 level or the AMPF 2 level, which were supposed to deliver two aircraft conflict to the ASF 3 . In Europe, ARC2000 [K + 89, 1 Maneuver Option Manager 2 Airspace Manager Planning functions 3 Automated Separation Function [START_REF] Fron | Arc 2000 : Automatic radar control[END_REF] or FREER [DHN97] used 1-to-n strategies to solve conflicts which were unable to deal with large conflicts. A 1-to-n solver using an on board token allocation strategy was developed by Granger and Durand [START_REF] Durand | Faces: a free flight autonomous and coordinated embarked solver[END_REF], but is not efficient with high traffic densities [START_REF] Granger | Detection et rsolution de conflits ariens : modlisations et analyse[END_REF].

Theoretical approaches using potential fields models [START_REF] Tomlin | Conflict resolution for air traffic management: a case study in multi-agent hybrid systems[END_REF][START_REF] Gosh | Maneuver design for multiple aircraft conflict resolution[END_REF], neural networks [START_REF] Durand | Collision avoidance using neural networks learned by genetic algorithms[END_REF][START_REF] Durand | Neural nets trained by genetic algorithms for collision avoidance[END_REF], linear programming [M 94], or semidefinite programming [START_REF] Frazzoli | Aircraft conflict resolution via semidefinite programming[END_REF][START_REF] Dodin | Rsolution de conflits via la programmation semie-dfinie[END_REF] are limited by the size of the problem considered (never more than 5 aircraft). Furthermore, assumptions made on trajectories forecast are generally not realistic (use of constant speed, no uncertainty. . . ).

A first definition of clusters was given by Durand and Alliot [START_REF] Durand | An optimizing conflict solver for atc[END_REF][START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF] in the middle of the nineties. A global solver using genetic algorithm was presented in the CATS/OPAS ( [ABDM97, DA97]) simulator. It showed its efficiency on conflicts involving up to 30 aircraft.

The size of clusters widely depends on the detection window, the uncertainty on aircraft speed and the separation standards. Dealing with small clusters is essential to ensure the resolution efficiency.

In order to have a better understanding on the cluster problem, we first detail mathematically different possible cluster configurations. The different configurations are then tested with the CATS/OPAS arithmetic simulator on real traffic data in a second part. The size of clusters also depends on the routes followed by the aircraft and on uncertainties. In the third part, the influence of uncertainty on cluster size and the feedback clusters process are discussed.

Cluster analysis 1.Definitions

A cluster is a transitive closing on conflicting pairs of aircraft. For example (see figure 1), if aircraft A and B are in conflict and aircraft B and C are in conflict then A, B, C belong to the same cluster. A cluster cannot be defined only by its size. For example there are 2 sorts of 3 aircraft cluster (see figure 2).

A cluster can be represented by a graph : the nodes are the aircraft and the edges are the conflicts. A cluster or a graph can be characterized by:

• the number of aircraft or nodes of the graph;

• the number of conflicts or edges;

• the diameter of the cluster or the graph (see definition 1.2);

• the graphic sequence of the graph or the cluster (see definition 1.3).

Definition 1.1 -The distance between two nodes of a graph is the minimum number of edges connecting the two nodes.

Definition 1.2 -The diameter of a graph is the maximum distance among every nodes pair.

Theorem 1.1 -The range of the diameter of a connected graph of size

n > 1 is [1, n -1].
This result is self-evident. The following theorem proved by Erdos and Gallai [EG60] can be used to define an algorithm that counts the number of possible graphic sequences.

Theorem 1.4 -

A sequence degree d = (d 1 , d 2 , ..., d n ) is graphic if and only if: ∀r ∈ [1, n -1], r i=1 d i -r (r -1) ≤ n i=r+1 min(r, d i )
Table 1 details for a given cluster size the minimum and maximum number of conflicts, the total number of graphic sequences, its logarithm and the number of conflicts that gives the maximum number of graphic sequences. The growth of the number of graphic sequences is exponential. Even for a five aircraft cluster, there are 19 different possible configurations. This result gives a better understanding of the inability for human controller to handle clusters involving more than 2 or 3 aircraft.

Experimental results

Results presented in this part were computed with the CATS/OPAS traffic simulator [ABDM97] on a French loaded traffic day (May 21 th 1999). The uncertainties used for the trajectory forecast were 5% on the horizontal speed and 15% on the climbing and descending rates. The time window was set to 8 minutes and the detection was done every 2 minutes.

Two scenari were tested: the first one used direct routes from origin to destination and the second one standard routes given by the flight plans. No pre-regulation is done on flight plans.

Table 2 and 3 detail the structure of different clusters with the standard routes scenario (2) and direct routes scenario (3). The possible and observed graphic sequences are given column 2 and 3. The ratio (percentage of observed/possible structures) is given column 4. The number and percentage of clusters of each size are detailed column 5 and 6. For standard routes, up to 5 aircraft, every configuration is represented. For direct routes, there are still 16 out of 19 different five aircraft clusters. It is important to notice that the graphic sequence of a cluster only measures the interactions between aircraft. No information on the cluster geometry is considered (number of leveled, climbing descending aircraft). If such criteria were taken into account, the number of different clusters would probably become huge and no expert system or human controller can handle so many situations. This explains why the CATS/OPAS global solver cannot be used as a decision making tool for big clusters. For such situations, it is most of the time too difficult for a human to understand the solution suggested by the solver.

Figures 2,3 and 4 detail for 3, 4 and 5 aircraft the possible graphic sequences. The (mean/standard deviation) of the number of conflict per aircraft are given below each graphic sequence. The clusters are sorted according to the number of aircraft they contain. At the end of each line, the ratios of observed clusters of the line over the total number of clusters are given.

For the 3 sizes of clusters, the number of conflict is lower with direct routes than with standard routes.

Finally, the mean diameters of clusters are smaller with the standard route scenario than with the direct routes scenario (figure 5). With direct routes, clusters are less dense and more spreaded out than with standard routes. The simulation does not include any regulation on the traffic, but these results show that the existing routes network may not be optimal for conflict resolution.

Sensitivity to uncertainties and resolution

Trajectories forecast uncertainties

When uncertainty is added on the horizontal speed and on the climbing and descending rates, more conflicts are detected and the size of clusters increases. Table 3.1 details the number of maneuvers, the size of the biggest cluster and the remaining conflicts for the different scenari. The global optimisation solver can easily handle clusters of size lower than 30 aircraft which is confirmed by the simulations. These results point up the necessity of a good trajectory forecast in order to limit the size of clusters and build an efficient solver.

Scenario

The feedback clusters effect

When a maneuver is decided for an aircraft, the solver must check that the new trajectory defined does not create a new conflict with other modified or existing trajectories. For example, if cluster A -B -C has been solved by changing the trajectory of aircraft B, and if this new trajectory creates a conflict with aircraft D, then the solver has to consider the new cluster A -B -C -D. Clusters stemmed from this process will be called feedback clusters.

Table 4 gives for the two scenari the total number of clusters, the number of feedback clusters and the ratios feedback clusters/total clusters. The global ratio of feedback clusters/total clusters is close too 28% with standard routes and only 7% with direct routes. The ratio increases with the size of the clusters. As large clusters appear in dense areas the feedback process occurs more often.

The feedback cluster ratio increases with uncertainty: table 5 shows for the direct route scenario, with three hypothesis of uncertainties, the total number of clusters, the number of feedback clusters and the corresponding ratios. Uncertainty has an important influence on the feedback process. The more uncertainty, the more feedback clusters. With direct routes, when increasing the uncertainty on climbing and descending rates to 50% and on horizontal speed to 20% (see table 6), the total number of clusters and the ratios of feedback clusters become close to those observed with standard routes (with different uncertainties). However, for small clusters, the ratios of feedback clusters are lower for direct routes.

Conclusion

The diversity of cluster structures can be modelled by graphic sequences. The number of different clusters increases exponentially with the size of the cluster. Experimental results have shown that, up to 5 aircraft, all the structures were represented. For five aircraft, 19 different structures can be distinguished disregarding aircraft attitudes (climbing, descending, levelled). This diversity explains the structure of the existing air traffic control system: as a human controller cannot recognize and handle so many situations, different filters must guarantee that he will not have to manage more than 2 or 3 aircraft at a time in the same cluster.

The cluster problem must be taken into account when designing an automatic conflict solver. When clusters become too large (more than 20 to 30 aircraft), finding solutions becomes problematical.

Consequently, the key point for any future control system is to limit the cluster sizes of the problem considered. On the French example, it was shown that the routes network has an important influence on the size of clusters. Direct routes seem to be better than standard routes, but no other network has been tested yet. However, the most important point is to improve the trajectory forecast: actually, a good trajectory forecast is essential to limit clusters sizes. With the help of datalink communications, aircraft will be able to give to the control system their future positions for the next ten minutes with a good accuracy. Regarding the current traffic increase, this seems the only way to limit clusters sizes.

There is still much to do on the cluster problem: some cluster structures are probably easier to solve than others but how can they be caracterized ? Some large clusters are easier to split into small ones but how can they be recognized ? Future models will have to take into account the aircraft attitudes (climbing, descending, leveled). The complexity will be increased, but will the information added be useful to design new resolution strategies ?

  Figure 1: 3 aircraft cluster

Figure 2 :

 2 Figure 2: Possible graphic sequences for a 3 aircraft cluster

Figure 4 :Figure 5 :

 45 Figure 4: Possible graphic sequences for a 5 aircraft cluster

Definition 1.3 -A graphic sequence

	nb	nb	nb	nb		nb	nb
	of conf conf	graphic	log	conf	graph
	acft	min max	sequences		opt	sequ max
	2	1	1	1	0	1
	3	2	3	2	0.30	2/3
	4	3	6	6	0.78	3/4
	5	4	10	19	1.28 5/6/7
	6	5	15	68	4.22	8
	7	6	21	236	5.46	11
	8	7	28	863	6.76	16
	9	8	36	3137	8.05	20
	10	9	45	11636	9.36	24
	11	10	55	43306 10.68	29
	12	11	66	162728 12.00	35
	13	12	78	614142 13.33	41
	14	13	91	2330454 14.66	48
	15	14	105	8875656 16.00	55
	16	15	120	33924699 17.34	63
	17	16	136	130038017 18.68	71
	18	17	153	499753560 20.03	79
	19	18	171	1924912505 21.38	89
	20	19	190	7429159770 22.73	98
	21	20	210 28723877046 24.08	108
	d = (d 1 , d 2 , ..., d n ) is
	a sequence of numbers which can be the degree sequence
	of some graph.
	1.2 Theoretical results
	Theorem 1.2 -A graphic sequence d = (d 1 , d 2 , ..., d n )
	is a degree sequence of a connected graph if and only if n i=1 d i ≥ 2 (n -1).
	The proof of this last theorem can be found in [Gra02]
	and uses the following theorem [Die97].
	Theorem 1.3 -A connected n nodes graph is a tree if and
	only if it has exactly n -1 edges.

Table 1 :

 1

Number of possible graphic sequences as a function of the cluster size.

Table 2 :

 2 Standard Routes -5 and 15% of uncertainty -T

	clus	Graphic sequences	clusters
	size	possible	obs	ratio	nb	ratio
	2	1	1	100.00 18119 66.78
	3	2	2	100.00 4952 18.25
	4	6	6	100.00 1883	6.94
	5	19	19 100.00		3.34
	6	68	53	77.94		1.77
	7	236	98	41.53		1.00
	8	863	119 13.79		0.73
	9	3137	92	2.93		0.41
	10	11636	52	0.45	55	0.20
	11	43306	56	0.13	59	0.22
	12	162728	36	0.02	37	0.14
	13	614142	27	0.00	27	0.10
	14	2330454	10	0.00	10	0.04
	15	8875656	8	0.00	8	0.03
	16	33924699	9	0.00	9	0.03
	17	130038017	4	0.00	4	0.01
	18	499753560	3	0.00	3	0.01
	21	28723877046	1	0.00	1	0.00
	clus	Graphic sequences	clusters
	size possible obs	ratio	nb	ratio
	2	1	1	100.00 12118 80.52
	3	2	2	100.00 2049 13.61
	4	6	6	100.00	540	3.59
	5	19	16	84.21	217	1.44
	6	68	24	35.29		0.52
	7	236	13	5.51		0.15
	8	863	14	1.62		0.10
	9	3137	8	0.26	8	0.05
	11	43306	1	0.00	1	0.01
	Table 3: Direct Routes -5 and 15% of uncertainty -T w =
	8mn -δ = 2mn -21 mai 1999

w = 8mn -δ = 2mn -21 mai 1999

Table 4 :

 4 Ratios of feedback with 5 and 15% of uncertainty, for the two scenari

	cluster		Standard			Direct
	size	clusters feedback %age clusters feedback %age
		15922	3237	20.33	9970	527	5.29
		4981	1556	31.24	1855	206	11.11
		2244	860	38.32	586	93	15.87
		1356	535	39.45	218	43	19.72
		753	338	44.89	100	22	22.00
		541	235	43.44	42	14	33.33
		418	203	48.56	24	7	29.17
		311	154	49.52	14	2	14.29
		212	105	49.53	9	3	33.33
		190	97	51.05	5	1	20.00
		158	73	46.20	2	1	50.00
		111	42	37.84	1	1	100.00
		93	56	60.22	1	0	0.00
		85	40	47.06	2	1	50.00
		60	28	46.67	1	0	0.00
		63	32	50.79	1	0	0.00
		61	35	57.38
		57	35	61.40
		51	30	58.82
		34	15	44.12
	..	..	..	..
		5	1	20.00
		2	2	100.00
		1	0	0.00
		1	0	0.00
		1	0	0.00
	total	28152	7984	28.36	12831	921	7.18

Table 5 :

 5 Ratios of feedback for different uncertainties (direct routes)

Table 6 :

 6 Direct routes: large uncertainties

	[DHN97] Vu N Duong, Eric Hoffman, and Jean-Pierre
	Nicolaon. Initial results of investigation into
	autonomous aircraft concept (freer-1). In
	Procceedings of the 1rst USA/Europe Semi-
	nar, 1997.


Jean-Baptiste Gotteland
email: gotteland@recherche.enac.fr

Nicolas Durand
email: durand@recherche.enac.fr

Genetic Algorithms Applied to Airport Ground Traffic Optimization

Due to air traffic growth and especially hubs development, major European airports can easily become bottlenecks in the global air transportation network. Therefore, accurate models of airport traffic prediction become more and more necessary for ground controllers.

In this paper, a ground traffic simulation tool is proposed and applied to Roissy Charles De Gaulle airport. Two global optimization methods, using genetic algorithms at the airport, are developed to minimize taxiing time while respecting aircraft separation and runway capacities.

In order to compare the efficiency of the different methods, simulations are carried out on a one day traffic sample, and ground delay is correlated to the traffic density at the airport.

Introduction

Traffic delay due to airport congestion and ground operations becomes more and more penalizing in the total gateto-gate flight cycle. This phenomenon can be in a large part attributed to recent hubs development, as all departures and arrivals are tending to be scheduled at the same time. Moreover, many ATC 1 and ATM 2 inefficiencies can appear as a result of taxi queueing and take-off time uncertainty. As a consequence, airport ground traffic simulation tools become essential for airport designers and traffic managers. Highly detailed models of airport operations already exist, such as SIMMOD (SIMulation MODel, developed by the Federal Aviation Agency) or TAAM (Total Airspace and Airport Modeler, developed by the Preston Group [START_REF] Preston | TAAM Reference Manual[END_REF]). They can be useful to evaluate qualitatively the relative effects of various airport improvements.

In the context of the A-SMGCS concepts (Advanced Surface Movement Guidance & Control System), all the ground traffic information is supposed to be directly available for the ATC and many tools could be developed to assist ground controllers.

In this paper, a ground traffic simulation tool with a con-1 ATC: Air Traffic Control 2 ATM: Air Traffic Management flict resolution module is introduced in part 2. Parts 3 and 4 describe the different optimization methods used to find the best trajectory and the most adapted holding points for taxiing aircraft. In the last part, these optimization methods are applied on a one day traffic sample at Roissy Charles De Gaulle airport considering the operational airport configurations and speed uncertainties, and the results of the simulations are compared.

Ground traffic simulation

The ground traffic simulator used for this study is associated with a detailed airport model in order to assign a set of realistic paths to each aircraft. Traffic prediction is regularly computed according to a given speed uncertainty. Algorithms estimate the best maneuvers that can be executed by each aircraft to ensure taxiing separations and runway sequencing while optimizing a global criteria such as the total delay. This part briefly sums up the main topics concerning this ground traffic simulator.

Problem modelling

Inputs

The simulator inputs are the airport topological description, the aircraft types with their corresponding weight categories and a one day traffic sample at the airport.

A flight is described by its recorded flight-plan, containing the aircraft type, the runway and the gate used, and the departure or arrival time.

The aircraft type information enables the evaluation of the needed takeoff or landing distance, in order to select the possible runway entry or exit points. Moreover, it is also associated with a wake turbulence category, which determines the time separation between two aircraft (one, two or three minutes) on the runway.

Paths assignment

The airport is modelled as a graph linking its gates, taxiways and runways. Each link is assigned a cost, which is an evaluation of the time spent by an aircraft when it proceeds via this link, added to some particular penalty related to runway areas, gate access or undesirable directions. Classical graph algorithms can be used:

• The Dijkstra algorithm [START_REF] Ahuja | Network Flows, Theory, Algorithms and Applications[END_REF] computes all the best paths and their corresponding minimal cost from a given node to every other ones.

• A k shortest loopless path enumeration algorithm [START_REF] Martins | A new algorithm for ranking loopless paths[END_REF] can make use of the Dijkstra's result to find an acceptable set of different paths for each origindestination (figure 1).

Conflict detection

Separation rules

A first model of aircraft separation rules consisted in the assignment of a maximal capacity to each portion of the airport graph. These capacities were defined as a linear function of the length of the taxiway portion. However, this model had to be refined as it could not ensure realistic aircraft separations near taxiways intersections nor manage dependent taxiways. As a consequence, a discrete model of separation rules is defined as follows:

• The distance between two taxiing aircraft must never be less than 60 meters, except at the gate position.

• On the runway, a time separation of 1, 2 or 3 minutes (depending on the aircraft category) is necessary after a take off to clear next take off or landing from wake turbulence.

• When an aircraft is taking off or landing on a given runway, other aircraft can be taxiing in the same runway area only if they are behind it.

• When an aircraft is following an other one, its speed uncertainty is reduced (as the pilot won't go faster than the first one), so that the two aircraft are considered to be separated.

When the traffic prediction does not respect one of these rules, the two involved aircraft are called a conflicting pair.

Speed uncertainty

Each aircraft trajectory is predicted with a given speed uncertainty, fixed as a constant percentage of the nominal speed (which is function of procedures and turning rate).

Therefore, an aircraft is considered to occupy several possible positions at a given time (as shown figure 2).

Aircraft possible maneuvers

In order to ensure separations, aircraft trajectories can be modified by ground control orders. For each aircraft, a control order is defined by:

• The path that the aircraft must follow, chosen among the set of remaining possible paths for this aircraft;

• Optionally, a position where the aircraft must wait and a time ending this wait.

Simulation steps Prediction ∆ Tw Figure 4: Shifted windows Thus, the length of the segment representing the several possible positions of an aircraft is reduced when and where the aircraft is expected to wait, as the reference is a precise position and a precise end waiting time.

For an arriving aircraft, the waiting position can be in the air (this means that the aircraft is asked to slow down before landing). In this special case and to keep the simulation realistic, the delay cannot exceed λ seconds (usually λ = 30 seconds).

Clusters

In order to lower the complexity of the problem as often as possible, a transitive closure is applied on conflicting aircraft pairs and gives the different clusters of conflicting aircraft [START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF]. The different clusters will be solved independently at first.

If the separated resolution of two clusters creates new conflicting positions between them, the two clusters are merged and the resultant cluster is solved.

Simulation steps

Shifted windows

The simulator works with a prediction time window T w shifted every ∆ minutes.

When the simulated current time is t, flight-plans expected to land or to leave the gate in the time range [t; t+T w ] are activated: a set of appropriate paths is assigned to each of them and the resulting best trajectories are added to the ones of the already taxiing aircraft.

The conflict detection and resolution is performed in the time window and the resulting maneuvers are applied to build the new situation ∆ minutes later. The problem is then reconsidered at the new simulation step t + ∆ (see figure 4).

Global criteria to optimize

The global criteria to minimize is defined by the total taxiing time (including the time spent queueing for runway), increased by the time spent in lengthened trajectories.

With this definition, lengthening trajectory appears to be twice more penalizing than holding position.

The One-to-n method

This resolution method, called 1-to-n, treats a simplified problem. Aircraft are initially sorted and considered one after the other: first considered aircraft have priority on last considered ones.

The optimization problem is therefore reduced to one aircraft: the algorithm must find the best path and the best holding positions for the aircraft, taking into account the already considered aircraft trajectories.

The Branch & Bound algorithm

The 1-to-n method used in [START_REF] Pesic | Aircraft Ground Traffic Optimization using a Genetic algorithm[END_REF] to find the best solution for an aircraft avoiding some other ones had to be refined in order to take into account the limited time window T w and the speed uncertainties. However, this method is still modelled as a shortest path finding problem in a graph:

• A node of the graph represents a position of the aircraft, in a path p i at time t.

• If the position conflicts with another aircraft, the node have no son. Otherwise, it has two sons, corresponding to the two possible choices for the aircraft at each time step: go forward or hold position.

• The root nodes are defined by the current position of the aircraft at current time t 0 on each remaining path p. The initial cost (to reach each root node) is the time length of the path.

• The terminal nodes are the ones describing a non conflicting position of the aircraft at the end of the time window (t = t 0 + T w ).

As the problem is to minimize delay, it is always better (when possible) to go forward than to hold position: thus, a Branch & Bound (instead of an A*) algorithm with a best first search strategy ( [START_REF] Horst | Global Optimization: Deterministic Approaches[END_REF]) can quickly find the best solution for the aircraft (or prove that there is no solution)

Initial classification

As last considered aircraft are extremely penalized (they must avoid all first considered aircraft) the way to sort aircraft, called the classification, is a determining factor. Moreover, the classification must respect some constraints:

• As landing aircraft cannot hold position before exiting runway, their priority level must be higher than all taking off aircraft on the same runway.

• Departures queueing for runway should be sorted according to their position in the queue.

Therefore, a time T ai is computed for each aircraft a i as follows:

... • T ai = t 0 + t r for departures, where t 0 is the current simulation time and t r is the minimal remaining time until the runway.

• T ai = T landδ land for arrivals, where t land is the expected landing time and δ land the landing priority (1 hour for example).

Aircraft are sorted according to increasing values of T a .

Classification reviews

Due to the limited time window for traffic predictions and speed uncertainties influence, the resolution of the current situation (at time t 0 ) can be very different than the previous one (at time t 0 -∆), so that some aircraft can finally be in a position with no solution. A common example is the one of an aircraft leaded in front of another one without having priority on it.

In order to manage these cases, the classification can be reviewed: when the algorithm cannot find any solution for an aircraft a i , the penalizing aircraft a j (j < i) is identified, and the process is restarted after inserting a i before a j in the classification (see figure 5).

Genetic algorithms

Principles

Two methods using classical Genetic Algorithms [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF][START_REF] Michalewicz | Genetic algorithms + Data Structures = Evolution Programs[END_REF] with partially separable functions and their adapted crossover and mutation operators as presented by N. Durand in [START_REF] Durand | Genetic crossover operator for partially separable functions[END_REF] are developed and compared.

• The first method, called GA WAIT, looks for an optimal association of paths and holding positions for aircraft. It was already developed in [START_REF] Pesic | Aircraft Ground Traffic Optimization using a Genetic algorithm[END_REF], and was directly adapted to the new airport structure and traffic model.

• The second method, called GA SORT, looks for an optimal association of paths and priority levels for aircraft. It is combined with the Branch and Bound algorithm described in part 3.

As the two methods implementation is very similar, they are simultaneously described in the next sections.

Data structure

For both methods, an element of the population (a chromosome) must be a set of parameters describing the N aircraft trajectories during the considered time window.

In order to keep the partial separability of the problem, each aircraft trajectory must be associated with its specific parameters, as detailed below.

Encoding for GA WAIT

In this method, a chromosome is defined by 3N variables:

{(n i , t 0i , t 1i )} 1≤i≤N
n i is the number of the path that the aircraft a i must follow (chosen among the set of remaining possible paths for the aircraft) and [t 0i ; t 1i [ the laps of time during witch it must hold position (if t 0i ≥ t 1i , the aircraft a i does not stop).

Encoding for GA SORT

A chromosome is here defined by 2N variables:

{(n i , p i )} 1≤i≤N
n i is the number of the path that the aircraft a i must follow (as for GA WAIT), and p i is its priority level, in the range [1, ..., N ]. The encoding must ensure that

(p i ) (1≤i≤N ) is a deterministic classification: ∀i = j, p i = p j (1)
In this way, the N aircraft trajectories can be developed by running the Branch & Bound algorithm with the classification obtained with (p i ), restricted to one path per aircraft, given by (n i ).

Fitness function

Both methods are implemented with the same fitness function, in harmony with the global criteria to optimize (see part 2.3) constrained by the separation rules (see part 2.2).

The information needed for the evaluation of the N trajectories is stored in a fitness matrix M (see figure 6).

• m ii = d i +2.l i represents the delay penalty of aircraft a i . d i is the time spent in holding position and l i the time spent in lengthened trajectory due to the path n i .

• m ij is the number of conflicts between aircraft a i and a j .

To ensure that a chromosome describing a situation without any conflict is always better than a situation with a conflict, the fitness function F is defined as follows:

d + 2l i i c ij c ij 1 i j N N 1 i j Figure 6: Fitness matrix
If the fitness matrix is diagonal (no conflict),

F = 1 2 + 1 2 + N i=1 m ii (F > 1 2 )
Otherwise:

F = 1 2 + i<j m ij (F < 1 2 )
Moreover, a local fitness f i (to be minimised) can be assigned to aircraft a i :

f i = m ii + K j =i m ij
Where K ≫ m ii (for example: K = 10T w ). These local fitness are used to orient the crossover and mutation operators (see part 4.5).

Initial population

The initial population is generated with some heterogeneous random functions, as the best solutions are more expected to be in some particular subspaces:

• Values of n i for which l i (time lengthening) is low, for both encodings;

• Low waits for GA WAIT: t 0i ≃ t 1i ;

• Priority levels (p i ) for GA SORT near the ones corresponding to the classification by expected runway access time (T ai ), defined in part 3.2.

Crossover and mutation

The crossover operator used for both methods is the one described in [START_REF] Pesic | Aircraft Ground Traffic Optimization using a Genetic algorithm[END_REF] (see figure 7), applying the principles described in [START_REF] Durand | Genetic crossover operator for partially separable functions[END_REF] with the local fitness f i . The aim of this crossover operator is to increase the probability of generating children with a better fitness than their parents:

• when the local fitness f i of an aircraft a i is far greater in one of the two parents, the two children directly inherit of its associated variables (n i , t 0i , t 1i ) or (n i , p i ) (see figure 4.5).

• when the local fitness f i of the two parents are similar, the children inherit of a combination of the two versions of the associated variables.

The mutation operator randomly modifies the parameters of an aircraft having one of the worst local fitness.

In the GA SORT encoding, the new elements of the population must be formated before insertion in the population, in order to restore the classification constraints (1).

Fitness matrix computation

GA WAIT

To compute the fitness matrix (for a chromosome evaluation) in the GA WAIT method, the aircraft trajectories are built according to the 3N variables and a conflict detection between each pair of aircraft is performed.

The conflict detection is the most penalizing operation: each of the N aircraft can have δT w positions at each time step (δ is the speed uncertainty ratio). An elementary conflict detection (between two aircraft at a given time step) is implemented in O(δT w ) operations.

Therefore, the whole conflict detection is achieved in

O(T w * N 2 * (δT w )) = O(δN 2 T w
2 ) operations. As the time to solve the problem would be limited in an operational context, the evaluation process is accelerated by keeping the m ij values of the fitness matrix when a pair (i, j) is unchanged after a crossover or mutation operation.

GA SORT

In the GA SORT method, the Branch & Bound algorithm described in part 3 (limited to a single path per aircraft) is used to develop sequentially all aircraft trajectories, with the classification given by (p i ).

If there is no solution for an aircraft a i , the penalizing aircraft a j (p j < p i ) is identified by the Branch & Bound algorithm, and the m ij (and m ji ) elements of the fitness matrix are increased: in this method, there is no need of any supplementary conflict detection.

For each aircraft, the graph exploration visits O(T w 2 ) nodes. Each node development needs O(N ) elementary conflict detections (O(δT w ) operations for each of them). Thus, the fitness matrix computation is achieved in O(δN 2 T w 3 ) operations. However, this evaluation process cannot be accelerated, as the trajectory of an aircraft depends on all the trajectories of previously handled aircraft in the classification.

Sharing

In order to avoid a premature convergence of the algorithm around local optima, a clusterised sharing process is implemented as detailed in [START_REF] Yin | A fast genetic algorithm with sharing scheme using cluster analysis methods in multimodal function optimization[END_REF]. This sharing process applies the general principles introduced by Goldberg and Richardson [START_REF] Goldberg | Genetic Algorithms with Sharing for Multimodal Function Optimization[END_REF] but can be computed in O(n log(n)) operations instead of O(n 2 ) for classical sharing (n is the population size).

The distance between two chromosomes A, B is defined as follows :

• For GA WAIT encoding:

D(A, B) = N i=1 |d iA -d iB | + |l iA -l iB | 2N T w
• For GA SORT encoding:

D(A, B) = N i=1 |p iA -p iB | 2N + N i=1 |l iA -l iB | 2T w
Where:

• d iA = t 1iA -t 0iA (resp. d iB )
is the time spent holding position by aircraft a i in chromosome A (resp. B).

• l iA (resp. l iB ) is the time spent in lengthened trajectory in chromosome A (resp. B).

In this way, the distance increases with the difference of penalization assigned to aircraft.

Ending criteria

The ending criteria for GA WAIT and GA SORT is defined with two parameters N max , N opt :

• N max = 100 is the maximal number of generations.

If no available solution is found before, the GA will be aborted.

• N opt = 15 is the maximal number of solved generations with the same best fitness. When the fitness of the best element of the population is greater then 1 2 and is stationary during N opt generations, the GA is stopped. The two GA methods are run with a 200 chromosomes population size, a 60% crossover rate, a 15% mutation rate and the stochastic reminder without replacement selection principle, described in [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF].

Simulation results

Simulations

Figure 8 gives the mean generated delay as a function of the number of aircraft involved in the situation. It appears that the 1-to-n method always generates more delay than the two others, and that GA SORT is the most efficient method: solving the situation after sorting aircraft does not seem to be so penalizing, at the condition to get an adapted classification.

Figure 9 gives the mean number of aircraft simultaneously moving for each time period. It appears that the GA SORT method keeps a lower number of moving aircraft during heavy time periods: good resolutions of ground traffic situations allow to decrease delay, but also lead to better situations with less moving aircraft.

In order to observe the GA efficiency, figure 10 gives the mean number of generations required by the two GA methods as a function of the clusters' size and figure 11 gives the mean value of the best fitnesses of each resolution. It can be noticed that the GA SORT method never deals with clusters bigger than 15 aircraft. It generally requires less generations than the GA WAIT method and find a solution with a bigger fitness to the problem. This can be explained by the one wait per aircraft limitation of the GA WAIT method, which appears to be more penalizing than the classification restrictions needed by the GA SORT method.

Conclusions

The simulations carried out in this study show that ground traffic delay at big airports is not only due to runway capacity but is very sensitive to the way that aircraft are managed: at Roissy Charles De Gaulle, one of the busiest European airport, the generated delay was largely reduced in the simulations using genetic algorithms. Even if it is not yet possible to compare the solutions coming from the simulations to the reality (taxiing speed and followed paths can be very different and ground controllers actions are not even recorded), these results are promising : it seams profitable to develop some traffic management tools that could assist ground controllers and provide them some optimized solutions.

Moreover, it can be noticed that the whole airport model was easily upgraded with some new and more detailed data (push-back procedures, undesirable directions, new runways, aircraft wake turbulence categories, speed uncertainties, etc) without changing the algorithm itself. Thus, this kind of simulations can also be useful to evaluate benefits of new airport structures or new traffic procedures.

Further work will concentrate in refining the modelling and the global criteria to optimize, taking into account for example takeoff sequencing needs of approach sectors or priority levels for slotted departures.

Figure 2 :

 2 Figure 1: Airport graph

Figure 3 :

 3 Figure 3: Uncertainty reduction on holding points

  Figure 5: Classification reviews

  Figure 7: Crossover operator

Figure 9 :

 9 Figure 8: Generated delay

Figure 10 :

 10 Figure 10: Mean number of generations

Figure 11 :

 11 Figure 11: Mean value of best fitnesses