David Gianazza

Kévin Guittet

Evaluation of air traffic complexity metrics using neural networks and sector status

This paper presents an original method to evaluate air traffic complexity metrics. Several complexity indicators, found in the litterature, were implemented and computed, using recorded radar data as input. A principal component analysis (PCA) provides some results on the correlations between these indicators. Neural networks are then used to find a relationship between complexity indicators and the actual sector configurations. Assuming that the decisions to group or split sectors are somewhat related to the controllers workload, this method allows to identify which types of complexity indicators are significantly related to the actual workload.

I. INTRODUCTION

Much research has been conducted over the last decade to help understanding air traffic complexity and controller workload. The inadequacy of the aircraft count to appropriately reflect the traffic complexity has now been acknowledged for a long-time, and complementary indicators such as "traffic mix", "number of potential conflicts" and others, have been (and still are) designed. A linear combination of these variables, often referred to as dynamic density, is likely to better fit traffic complexity than individual indicators. It is used throughout most studies, where the correllation of a set of indicators with a quantifiable variable, assumed to represent the actual traffic complexity, is maximized.

A possible shortcoming of this methodology is that potentially non-linear relations between indicators are missed (see [START_REF] Hilburn | Toward a non-linear approach to modeling air traffic complexity[END_REF] and the concern of Eurocontrol when writing calls for proposals). But, more importantly, the choice of the dependent variable is crucial to determine how well complexity is actually measured. Indeed, physical activity, as used in [START_REF] Laudeman | Dynamic density: An air traffic management metric[END_REF] and [START_REF] Majumdar | The factors affecting airspace capacity in europe: A framework methodology based on cross sectional timeseries analysis using simulated controller workload data[END_REF], miss the important cognitive part of the controller activity. On the other hand, physiological indicators ( [START_REF] Crump | Review of stress in air traffic control: Its measurement and effects[END_REF], [START_REF] Averty | Evaluating a new index of mental workload in real atc situation using psychological measures[END_REF]) seem difficult to exploit and how well they relate to traffic complexity is unclear. Finally, widely used subjective ratings ( [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF], [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]) provide high quality data (as they obviously relate to the kind of complexity investigated), but are often seen as subject to biases (such as the recency effect denounced in [START_REF] Averty | Evaluating a new index of mental workload in real atc situation using psychological measures[END_REF], and the possibility of raters errors in the case of "over-the shoulder workload ratings" [START_REF] Mannings | The relationship between air traffic control events and measures of controller taskload and workload[END_REF]). In all of these cases, data are very expensive to collect, as they require the active participation of controllers. Databases are often small and might exhibit low variability, which may in turn harm the statistical relevance of the results. This phenomenom is acknowledged in [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], where the overfitting of data is clearly a consequence of a lack of observations rather than a misspecification of the neural network. Finally, as these complexity metrics may be used to design computer-assisted control tools or traffic management tools, and to help organizing airspace, it is surprising to notice that the question of the relevance of the complexity measured to the final goal is scarcely discussed. The question really is to understand which complexity is measured and how well it relates to the foreseen application (benchmarking, improvements in airspace organization, design of new tools...).

This paper is motivated by former studies on optimal airspace sector configurations ( [START_REF] Gianazza | Optimization of air traffic control sector configurations using tree search methods and genetic algorithms[END_REF], [START_REF] Gianazza | Optimal combinations of air traffic control sectors using classical and stochastic methods[END_REF]) and intends to improve the criterion used therein to evaluate sector configurations. The basic idea, introduced in [START_REF] Gianazza | Optimisation des flux de trafic aérien[END_REF], is that the decisions to split a sector, mostly taken when the controller is close to overload, are linked to traffic complexity and may therefore provide an acceptable dependent variable. Interestingly enough, collecting data on sector configurations does not require controllers active participation, as current outcomes from control centers can be used, while related flight informations are available from recorded radar tracks. As such, raw data needed in our study are noticeably cheap to collect and might be produced in large quantities. The price to pay is that these data are noisy, as we may not be sure that a sector splitting (resp. merging) decision is directly related to overload (resp. underload). Other factors might distort data, such as training of unexperienced controllers, meteorological hazards, military airspace use... However, we will assume that the impact of these phenomena on the accuracy of the results is limited, particularly because of the kind of complexity we are looking at here. Indeed, this work is conducted in the perspective of future pre-tactical applications (e.g. sector planning) and thus does not ask for as much details as studies of instantaneous workload would (on the opposite, benchmarking of ATC centers would require an even coarser granularity, as indicators are averaged on wide temporal and geographical horizons [START_REF] Flynn | Traffic complexity indicators and sector typology analysis of u.s. and european centres[END_REF], [START_REF]Note de synthèse sur l'indicateur de complexité pru[END_REF]). To investigate the link between complexity indicators and sector configurations, we use neural networks, as nonlinear interactions are suspected.

The paper is organized as follows. Section II briefly describes the indicators used throughout the study, while section III presents the raw data from which the final database is built. A Principal Component Analysis (PCA) is then performed in section IV to restrain the dimensionality of the data. Neural networks are introduced in section V and their results are presented and discussed in section VI. Section VII concludes.

II. AIR TRAFFIC COMPLEXITY INDICATORS

The accuracy of the results of a study related to air traffic complexity is strongly dependent of the diversity and quality of the chosen individual complexity indicators. Many have been suggested to help describe the controllers workload, and it is hardly possible to implement the entire pool. In order to limit the number of variables to be (re)programmed and present indicators that are representative of the dynamic density litterature, we focused on the ones selected by Kopardekar [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF] in its unified complexity metric 1 . These indicators, such as references of studies where they were used and where definitions may be found, are presented in Table I. We also implemented several indicators inspired by studies conducted elsewhere in the SDER (former CENA). Definitions are indicated in appendix 2 . Finally, we also used incoming flows as explanatory variables, as they may be a significant factor in the decision to split (or merge) a sector.

III. INPUT DATA

The indicators are computed every round minute of the day, using recorded radar data, environnement data (sector description), and recorded sector configurations of the five french ATC centers. The sector configurations are recorded every round minute of the day, which explains our choice concerning the frequency at which we compute the indicators. 1 Though we were not always able to find an explicit formula, and thus missed seamingly important indicators like, e.g., "MET_airspace structure". Note that this difficulty to get clear definitions is also reported by Eurocontrol in [START_REF]Cognitive complexity in air traffic control, a litterature review[END_REF]. 2 Further informations and discussions about indicators are to be found in the internal note [START_REF] Guittet | Analyse descriptive des indicateurs de complexité du trafic aérien à partir des données image et courage[END_REF]. Radar data is available in several forms: records made by each center, with one position every twelve seconds, in average, and a global record of the five centers, with one position every three minutes. Several months of global records were available, whereas the centers local records were not readily available, at least for a sufficiently long period of time. So we used the global records (made by the IMAGE system), and interpolated the aircraft positions in order to get one position per minute. As many trajectory changes may occur within three minutes of flight, the computed positions are not highly accurate, and this may introduce a bias in the indicators values. However, this bias is most probably of small importance in our problem: we just want to predict when a sector will be merged into another one, or split in several smaller sectors. We are not considering the instant workload, which may require a very high level of accuracy on the aircraft position, speed, and so on. To be sure that this bias is small, we should compare the computed positions, and maybe also the indicators values, using local centers records, and global records, on small data samples. This is left for future work.

Several months of recorded traffic are available. However, considering the volume of data, it would be tedious to run several experimentations on very large data samples. So, we have restricted our choice, at least for the moment, to one day of traffic (1st june, 2003). Once we have found the most significative complexity indicators, it will be possible to re-train the neural network on larger data samples.

On the chosen day, 103 different sectors were armed. The term "sector" means here either an elementary sector, or a set of elementary sectors merged together, and handled on a single controller's working position. The air traffic complexity indicators were computed for each of these sectors, every minute of the day, together with the sector status (merged, armed, or split). This data was split into two sets : about sixty percent was randomly selected in order to train the neural network, and the rest was used to test the trained network on fresh data.

This single day of traffic already provides a big volume of data, as detailed in table II, with a great diversity of geographic sectors, and with enough data in each class of sector status. Before applying the neural network to complexity indicators and sector statuses, let us first discuss the correlations between the indicators, using a principal component analysis method.

IV. PRINCIPAL COMPONENT ANALYSIS

Including incoming flows, we end up with 27 complexity indicators. Given the neural network greediness in numbers of parameters and the high multicollinearity of the data, we will use principal components rather than individual indicators for our experiments. 6 main components were identified (corresponding to eigenvalues greater than 1), that covered more than 76 % of the variance of the data set. These components are interpreted below. C1: Eigenvalue 12.6, 46.7 % of the variance of the data set. Appart from vprox_1 and vprox_2, all variables are strongly (and positively) correlated with this component, explaining its high associated eigenvalue. This component may be seen as a "size factor", and we follow [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF] on the term "Overall monitoring". This component is strongly representative of the aircraft count. C2: Eigenvalue 2.78, 10 % of the variance of the data set. This component is strongly correlated with avg_vs, σ 2 gs and σ gs /gs (resp. 0.70, 0.69 and 0.68). Accounting for the impact on ground speed of vertical evolution of aircraft, this component may be seen as related to the ground speed variance, and the aircraft vertical evolutions. C3: Eigenvalue 1.96, 7.3 % of the variance of the data set. This component is mainly correlated with incoming flows. C4: Eigenvalue 1.25, 4.6 % of the variance of the data set. Appart from variables directly related to the traffic volume (N , inter_hori), the correlation on this component is high with insen_c and CREED indicators, and thus might be related to converging flows and anticipation of conflicts. C5: Eigenvalue 1.06, 3.9 % of the variance of the data set. This component is strongly correlated with Div and insen_p, and seems therefore mostly linked with divergent flows. C6: Eigenvalue 1.03, 3.8 % of the variance of the data set. This component is strongly correlated with the vertical proximity measures ( [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF]), and could stand for the monitoring of vertical separation (near the minimas). Notice that we extracted only 6 components, thus significatively less than the 12 components (briefly) described in [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF]. This might be explained by the lack, at the point of the project, of indicators related to the sector geometry.

V. NEURAL NETWORKS

A. General presentation

Artificial neural networks are algorithms inspired from the biological neurons and synaptic links. An artificial neural network is a graph, with vertices (neurons, or units) and edges (connections) between vertices. There are many types of such networks, associated to a wide range of applications: pattern recognition (see [START_REF] Bishop | Neural networks for pattern recognition[END_REF] and [START_REF] Ripley | Pattern recognition and neural networks[END_REF]), control theory,... Beyond the similarities with the biological model, an artificial neural network may be viewed as a statistical processor, making probabilistic assumptions about data ( [START_REF] Jordan | Neural Networks[END_REF]). Some train data is used to determine a statistical model of the process which produced this data. Once correctly trained, the neural network uses this model to make predictions on new data.

Neural networks are closely related to the Bayesian probabilities. They may be used for unsupervised learning (density estimation problems), and, mainly, supervised learning problems (regression, classification). Density estimation is not in the scope of our paper, so we will not detail it. The aim of regression is to find a statistical model producing an output y from input variables (let us denote them by x), so that the output y is as close as possible to a target variable, which we shall denote by t. In the case of classification problems, the target variables represent class labels, and the aim is to assign each input vector x to a class.

We will use a specific class of neural networks, referred to as feed-forward networks, or multi-layer perceptrons (when the activation function is logistic). In such networks, the units (neurons) are arranged in fully-connected layers: an input layer, one or several hidden layers, and an output layer. Figure 1 shows an example of such a network. For a network with one hidden layer, the output vector y = (y 1 , ..., y k , ..., y q ) T is expressed as a function of the input vector x = (x 1 , ..., x i , ..., x p ) as follows:

y k = Ψ( q j=1 w jk Φ( p i=1 w ij x i + w 0j ) + w 0k ) (1)
where the w ij and w jk are weights assigned to the connections between the input layer and the hidden layer, and between the hidden layer and the output layer, respectively, and where w 0 j and w 0 k are biases (or threshold values in the activation of a unit). Φ is an activation function, applied to the weighted output of the preceding layer (in that case, the input layer), and Ψ is a function applied, by each output unit, to the weighted sum of the activations of the hidden layer. This expression can be generalized to networks with several hidden layers.

The output errori.e. the difference between the target values t and the output y computed by the network -will depend on the parameters w (weights and biases). The training aims at choosing these parameters, so as to minimize a chosen function of the output error.

In the case of regression, the minimized function is the sum of quadratic errors. For classification problems, it is best to consider a log-likelihood function. The network is then designed with one output unit per class. When the output of such an unit is 1, the input x is assigned to the class corresponding to the unit, and when the output is 0, it is not. Let us consider a problem with C classes. The log-likelihood function minimized during training is the following, known as cross entropy:

E(w) = - N n=1 C k=1 t (n) k ln(y (n) k ) (2)
where t (n) and y (n) are the n th target and output vectors, respectively.

Several optimization methods may be used to minimize E(w), when training the network. Let us cite backpropagation, which consists in successive modifications of the weights assigned to the connections between the layers, starting with the output layer. These modifications take account of the relative importance of each weight in the output error variations. Other local optimizations using the gradient of the error (BFGS, conjugate gradients, for example) are also widely used. A variety of global optimization methods (simulated annealing, evolutionary algorithms) are also proposed in the literature. These global methods usually perform better than the local methods when there are many local minima for the error function, but they are generally much slower.

B. Neural networks applied to our problem

For our problem, we have chosen three-layers feedforward networks, denoted I α H β O γ in the rest of the paper, with α units in the input layer, β units in the hidden layer, and γ units in the output layer. The input variables are normalized, by substracting the mean value and dividing by the standard deviation.

There are many possible choices for the functions Ψ and Φ, depending on the problem being adressed. A common choice for Φ is the logistic function :

Φ(z) = 1 1 + e -z (3)
This is the activation function that was used in our experiments. As we adress a classification problem -assign each input vector (a list of complexity indicators values) to a class representing the sector status (merged, armed, or split) -we have chosen to minimize the cross entropy function. Therefore, the transfer function Ψ applied to the output layer, must be the softmax function:

Ψ(z k ) = e z k C m=1 e zm (4)
The nnet package of the R language was used (see http://www.r-project.org/ for details on the R language and environnment). In this package, developped by Pr B. D. Ripley, a quasi-newton minimization method (BFGS) is used for the network's learning. The parameters of the nnet tool are the range parameter (default 0.5), defining the range into which the initial random weights are chosen, the parameter for weight decay (default 0), and the maximum number of iterations. The training stops either if the fit criterion (the cross entropy in our case) falls below a chosen parameter abstol (default 1.0e-4), or if the improvement of the fit criterion is less than 1 -reltol (the default value for reltol is 1.0e-8).

Several combinations of air traffic complexity indicators, or of principal components, will be tested. The number of input units of the network will be chosen equal to the cardinal of the evaluated set of indicators (or components). The number of hidden units is 15 (this choice is discussed later). The output layer is made of three units, one for each class (merged, armed, or split).

So a target vector t (n) with value (1, 0, 0) means that the considered sector was merged with other sectors when the n th measure of the vector of complexity indicators was made. Armed sectors will be represented by (0, 1, 0). A value of (0, 0, 1) will mean that the sector was previously split in two or more sectors at the time x was measured. Of course, the actual output of the neural network will not be exact values 0 or 1. It will be triples (a, b, c) of floating-point values between 0 and 1, each value being the probability to belong to a class. The input vector x (n) will be assigned to the class of highest probability.

C. Evaluation of the neural network's outputs

A well-known problem, when using neural networks (or other regression methods), is overfitting: with enough parameters and enough training cycles, it is always possible to find a good fit for a given data set. So one may find a perfect fit for a chosen data sample, and then feel disappointed when the trained network makes wrong predictions on fresh data. So, we will systematically proceed as follows: train the network on a randomly chosen data sample (called train), then check the results, first on the same data sample, and second on a fresh data sample (called test), that was not used for the training.

In order to evaluate the outputs of several different models, we have to compare the neural networks predictions to the actual target values. We may use the fit criterion (cross entropy) but it does not reflect the influence of the number of weights (and biases) in the neural network. It is known (see [START_REF] Jordan | Neural Networks[END_REF]) that a network with too few weights may not be able to capture all the variations of the response to the input x, whereas a network with too many weights will more likely be subject to overfitting. In the next sections, we will compare several sets of input variables, of various sizes. Consequently, the number of weights in the network will not remain constant, and this variation will bias the results.

We will therefore use the Akaike information criterion ( [START_REF] Greene | Econometrics analysis[END_REF]): AIC = 2λ -2ln(L), where λ is the number of unadjusted parameters of the model (i.e. the number of weights and biases of the network), and ln(L) is the log-likelihood. This criterion is strongly related to information theory, and more specifically to the Kullback-Leibler distance (K-L) between a candidate model and the "true" model. In our case, the AIC is written as follows:

AIC = 2λ -2 N n=1 C k=1 t (n) k ln(y (n) k ) (5)
One should be aware that AIC is a relative criterion, which can only be used to compare a set of candidate models relative to a same "true" model: as this true model is unknown, the corresponding term in the K-L distance was considered as a constant and dropped, in the AIC. As we would like to compare predictions made on the train test and on the test set, which are of different size, we will divide the AIC by N , the number of data items, and use: AIC avg = AIC N . In addition to the numerical results provided by the Akaike information criterion, we shall also consider the global proportion of correctly classified input vectors, and also the percentage of correct classifications for each class. One must be aware, however, that the rate of correct classifications is not the criterion being maximized by the neural network, so we should remain cautious when comparing the different classification rates. However, these percentages are easily understandable and may allow us to make some interesting statements on the results.

VI. RESULTS

A. Preliminary results, discussion on the parameters

Before making statements about how significant the indicators (or combinations of indicators) are, let us first check if the chosen network is efficient on our problem, and consider how to choose the network and training parameters.

Figure 2 shows the evolution of the cross entropy criterion during the network's training. The input variables were the six main components found by the PCA, and also the sector's volume, so the input layer of the network had 7 units. We have chosen a hidden layer with 15 units. The output layer had 3 units, each one representing a class, as explained in section V. The maximum number of cycles was set to 1500, but the training actually stopped at 1010 cycles, because the algorithm was unable to significantly improve the cross entropy criterion.

Several questions may arise, concerning the parameters choices, and their influence on the results. In this study, we will use the default values for the following nnet parameters: the stop criteria abstol and reltol, and the weight decay parameter, set to decay = 0. The range in which the network's weights and biases are initialized depends on the range of the input variables, that we have normalized. A parametric study with all components plus the sector volume as input and with 15 units in the hidden layer shows that the best results, averaged on two runs for each value of the range parameter, were obtained with a range value of 0.4. So we shall use 0.4 for the nnet range parameter.

It may have been useful to make another parametric study, trying to find the optimum number of hidden units by minimizing the AIC, but this is rather timeconsuming (it should be repeated for each set of input variables). So we have chosen, after a few trials, a hidden layer of 15 units.

B. Model selection

Let us now select the best model among several candidate models. A model is a trained neural network and set of input variables that we expect to provide a good explanation of the sector status (merged, armed, or split). We use the PCA components plus the sector volume as input variables. The sector volume was not analyzed in the PCA, but as we have not implemeted any indicator using the sector geometry, we will use the volume as a proxy for metrics such as "space available around conflicts" or "distance to sector boundary".

An iterative approach is used: we shall first use component C 1 (representative of the number of aircraft) as input to the neural network, then add the volume, and continue with the five other main components, successively added in the order found by the PCA. At last, we will use all the 27 components, and the sector volume, as input. The AIC avg criterion (see section V-C) is used to select the best model.

For each set of input variables, five training runs are made. The reason for this is the following: the training method of nnet is a local optimization method which starts at a randomly chosen point (the initial weights), and which follows the steepest descent of the error function being optimized. This error function may have several local minima, so choosing different initial weight values may lead to different results. Although these local minima are often fairly close, several runs will comfort our results. Figure 3 shows, for train and test data, how the AIC avg criterion evolves when adding components to the set of input variables. The mean values, averaged on the five runs, are presented. This figure should be interpreted as follows: when the AIC significantly decreases when introducing a new variable in the model, this means that this variable improves the prediction of the sector status. When the AIC increases or remains constant, this means that the benefit provided by the additional variables is offset by the complexity it implies on the model (increase in the number of parameters) 3 . The two curves, train and test, are fairly close, except at the last point (V + AllC), which corresponds to the model {V ; C 1 ; ...; C 27 }. For this last point, the AIC criterion is improved on train data, whereas it remains nearly the same on test data, which shows a little overfitting of the neural network on train data. Otherwise, the results on train and test data are quite consistent, so the neural networks are able to generalize efficiently on fresh data.

The AIC criterion is significantly improved when adding the sector volume to the component C 1 (aircraft count). Component C 2 , which is mostly related the speed variance and the aircraft vertical evolutions also improves the criterion, although much less than the sector volume. Small improvements are brought by component C 3 (incoming flows), and also by C 4 , mostly related to the converging flows and the conflict anticipation. The components C 5 (diverging flows) and C 6 (monitoring of vertical proximity) bring no significant improvement. The use of the other components does not improve the prediction, as shown by the last point of the test curve.

It might look surprising that the aircraft proximity (horizontal or vertical) is of little influence on the sector status explanation. But, as was already stated in [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], the aircraft have already been separated (in the other dimension) before the proximity situation occurs, which explains why it is not very useful. The anticipation of future aircraft proximity (component C 4 ), is more significant.

C. Classification rates

Tables III and IV give the proportions of correct classifications made by the neural networks (the best of the five runs), using train data as input, or test data, respectively. The second column of each table shows the global rate of correct classifications. The three last columns detail the results for each class (merged, armed, or split). As previously stated on figure 3, the results on test data are quite consistent with the results on train data. When using only the component C 1 (aircraft count) as input, the neural network already makes more then 72% of correct predictions. However, the rate for the armed class is less than 45%, which is not very good. Once again, we see that the use of the sector volume (denoted V ) greatly improves the prediction of the sector status, with more than 79% of correct classifications, globally, and a rate of about 59% for the armed class. The classification rates when adding C 2 climb to 81% for the global rate, and about 60% for the armed class. Adding C 3 does not increase the rates, although it improved the AIC. C 4 slightly improves the results, whereas C 5 and C 6 provide no improvement. At most, when considering the sector volume and the 4 first components, we have about 82% of correct classifications, and about 63% of correct classifications for the armed class.

The last line of both tables, with all components and the volume, shows the best classification rates, whereas the AIC showed no improvement in the same case. The apparently good results of the last line are certainly due to the higher number of network's parameters than in the previous models. They do not mean that the model with all components is the best one.

In order to verify this assumption, we have tried another network for the {V ; C 1 ; ...; C 4 } model, with about the same number of parameters than the {V ; C 1 ; ...; C 27 } model. This last model was assessed with a network I 28 H 15 O 3 , where there were 483 weights and biases. Let us take a I 5 H 53 O 3 network, with 53 units in the hidden layer, for the {V ; C 1 ; ...; C 4 } model. This network has 480 parameters. Table V shows that the classification results, with this network, are better than the ones obtained with all components and the volume in table III When considering the detailed results, for each class, on the three tables, we see that the merged and split classes have better classification rates than the armed class. This is not a surprise, as the cloud of points representing the measures of the armed class, in the variables space, is located "between" the clouds representing the two other classes. The neural network aims at finding the frontiers between these clouds, and this is more difficult with two frontiers instead of one.

To conclude this section, we have seen that the best model is able to classify correctly about 84% of the input vectors. We may still slightly improve these results by running a parametric study on the number of hidden units, for this model, still minimizing the AIC criterion. But this is rather time-consuming, and our results are already quite good, considering that our data is rather noisy. Let us remind that we have assumed that the decisions to split or to merge sectors had a single cause, that is the workload. This may not always be the case in reality, where there could be other reasons: controllers training, hardware failures, and so on. In a next step of our research, we may try to improve our results by filtering the data -at least the sector split decisions -using an approximation of the workload (the number of aircraft for example).

D. Discussion on the results, comparison with other works

To summarize the previous results, the AIC criterion (figure3) allowed us to select the best model among the ones we have tested. This model uses a subset of only 5 input variables {V ; C 1 ; ...; C 4 }, among the 28 that we have considered.

The proportions of correct classifications (tables III and IV) show that the use of the sector volume significantly improves the network's prediction, by 6.5%, when compared to the sole aircraft count C 1 . There is a relatively small improvement (around 3%), between the sector volume + aircraft count model and the best model. The other components bring no significant improvement.

Other works, like [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF] and [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF] for example, have already stated that the aircraft count model is not very efficient in predicting the controller's workload. In these studies, the improvement brought by other complexity indicators was much higher than what is shown in our results. There may be several explanations to this.

In [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], Chatterji and Sridhar show workload prediction rates (on test data) ranging from less than 16% to 54%, for their medium workload class, and from 0% to 100% for their high workload class. However, their results with test data were not consistent with the ones obtained with train data. They honestly state that they were unable to draw reasonable conclusions from these results, as they had too few measures in the medium and high workload classes.

In [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF], Kopardekar and Magyarits apply a linear regression on subjective complexity ratings. The results were not given in percentages of correct classifications. The R 2 criterion of the regression was used to compare the candidate models. The results presented in [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF] show some significant differences of R 2 values between the aircraft count model and the model based on dynamic density (which is viewed as a linear combination of several indicators). Kopardekar and Magyarits also lacked measures, but for very low and high workload traffic.

Both studies ( [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF], [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF]) use subjective complexity ratings, provided by air traffic controllers who assessed the traffic complexity of several traffic samples, only during periods when the sector was armed 4 . In our study, we use the actual sector status, recorded for several elementary sectors and groups of sectors, assuming that the sector status is related to the controllers workload. So we have only three levels of workload: low, when the sector is merged with other sectors, normal use, when the sector is armed, and too high when the sector has been split into smaller sectors. We have the feeling that the studies [START_REF] Chatterji | Measures for air traffic controller workload prediction[END_REF] and [START_REF] Kopardekar | Measurement and prediction of dynamic density[END_REF] focus on the normal use load interval, and that the use of complexity ratings within this domain magnifies the observed phenomenons.

Aside from these considerations, another, and more straightforward, explanation of the relatively small differences between our candidate models -when comparing the classification rates -resides in the nature of the data used in our study. We simply have a great number of measures for which there is no doubt as to which class they belong to, even when using the worst possible model. For example, when there are 2 aircraft in a given sector, it is most likely merged with other sectors. On the opposite, when there are 60 aircraft in a sector, it has certainly been split into several sectors. We may have exhibited much higher variations in the results, by computing the classification rates using only measures collected around the times when the sector configurations changed.

VII. CONCLUSION

In conclusion, we were able to select the best model, among several candidate models, establishing a functional relationship (equation V-A with the weights of the trained network) between the air traffic complexity indicators and the sector status. Assuming that the decisions to merge or to split a sector are statistically related to the controllers workload, this original method provides an objective way to validate the complexity metrics. Our method also has the advantage, in comparison to other methods, to use widely available data (sector configurations and radar tracks recorded by the ATC centers), noticeably cheap to collect, as the active participation of air traffic controllers is not required.

Neural networks, minimizing the cross entropy function of the output error, showed good results, consistent on train and test data. The Akaike information criterion proved useful in selecting the best model, avoiding the bias due to the different number of parameters in the candidate models. In the iterative approach that was used, the highest improvement, when comparing to the aircraft count model, was brought when introducing the sector volume as a new input variable. Smaller improvements were provided by components C 2 (speed variance and aircraft vertical evolutions), C 3 (incoming flows), and C 4 (converging flows and conflict anticipation).

So far, we have only considered the PCA components in our study. The next step of our research may be to select a subset of individual indicators, issued from the components of the best model, and re-iterate the approach presented in this paper. This would provide a more direct and simple relationship between the indicators and the sector status, by avoiding to compute the components from all the indicators.

The neural network approach used in this study seems appropriate for the granularity we are interested in and the foreseen applications, either strategical (sector design) or pre-tactical (sector planning). We are also fairly confident that decisions to split or merge sectors may allow to assess the instantaneous workload as well, and could therefore be used to improve tactical tools (PRESAGE). To this end, other statistical methods should be investigated to take into account the serial correlation of sector status, looking closely at the sector splitting times. We plan to tackle this issue in a close future, using dynamic discrete choice models.

Finally, another issue that we intend to address, in relation to the complexity indicators, is the prediction of optimal sector configurations. Previous works ( [START_REF] Gianazza | Optimization of air traffic control sector configurations using tree search methods and genetic algorithms[END_REF], [START_REF] Gianazza | Optimal combinations of air traffic control sectors using classical and stochastic methods[END_REF]) proposed several algorithms to compute optimal sector configurations, using sector capacities and incoming flows. The output of the neural network is a triple of probabilities, allowing to decide when a sector should be split, or merged. We may derive a realistic workload indicator -and also threshold values -from these probabilities, which could be used to compute optimal sector configurations. number of aircraft, it is unclear whether they actually are "complexity" or "simplicity" indicators. We thus define a last pair of indicators, insen_c and insen_d as insen_c = Conv 2 Sd+ and insen_d = Div 2 Sd-.

Modified PRU metrics

The work conducted by SDER-RFM for the Performance Research Unit (citeRFM), though initially designed to compare ATC centers on a daily basis, inspired the following indicators :

• inter_hori: number of potential crossings (irrespective of the aircraft direction on their trajectories) with angle greater than 20 degrees.

• avg_vs: this is simply the average vertical speed of controlled aircraft.

Metrics inspired from the CREED project

The work of P. Averty on conflict detection [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] inspired a set of indicators. One of the ideas in [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] is that conflict perception is "plannar". The author thus defines for converging pairs of aircraft the following quantities

• Ed : minimum horizontal distance between aircraft.

• Efl : horizontal distance when the aircraft are vertically separated (after the crossing). • Da : the "anticipation degree", i.e. the distance between the faster aircraft and the intersection of the aircraft trajectories (in the horizontal plan). We replace this variable to a modified Da, DaC, which stands for the greater distance between one of the aircraft and the point where, horizontally, the distance between aircraft is the smallest. For explanations about this substitution, we refer to [START_REF] Guittet | Analyse descriptive des indicateurs de complexité du trafic aérien à partir des données image et courage[END_REF]. Originally, these quantities are defined to describe conflict perception. To translate the idea of [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] in terms of traffic complexity, we assume that a conflict is all the more critic that the expected separation (Ed and Efl) and the anticipation (DaC) are small. We thus set

creed = 1 α Da + (1 -α)(β Ed + (1 -β) Efl) , (18)
where α and β are parameters in [0; 1]6 . Finally, aircraft pairs considered in [START_REF] Averty | Conflit perception by atcs admits doubt but not inconsistency[END_REF] are such that vertical separation occurs prior to separation, as the converse situation is avoided as much as possible by controllers. Accordingly, the complexity associated with these latter pairs is likely to be greater and we distinguished the two kind of conflicts by summing the quantity fintroduced in [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF] on both sets of aircraft, thus creating two distinct indicators, creed_ok ("good pairs") and creed_pb ("bad pairs").

Fig. 1 .

 1 Fig. 1. Example of a feed-forward network with one hidden layer

Fig. 2 .

 2 Fig. 2. Evolution of the cross entropy criterion with the number of training cycles, using the main 6 components and the volume as input (network I 2 H 15 O 3 )

Fig. 3 .

 3 Fig. 3. Values of the AICavg criterion for the different sets of input variables

  • inter_vert: denote by n1, n2 et n3 the numbers of stable/climbing/descending aircraft. The indicator is then defined asinter_vert = (n1n2 + n2n3 + n1n3) (n1 + n2 + n3) .

  and IV.

	Set	Global	Merged	Armed	Split
	{V ; C 1 ; ..; C 4 } train	84.34%	89.74%	67.61%	91.88%
	{V ; C 1 ; ..; C 4 }test	83.84%	89.32%	67.26%	91.09%

TABLE V CORRECT

 V CLASSIFICATION RATES FOR THE {V ; C 1 ; ...; C 4 } MODEL, WITH A I 5 H 53 O 3 NETWORK

We could have used the Schwartz BIC criterion instead of the AIC. As it assigns a heavier penalty to additionnal parameters, it may have been easier to interpret.

This is not explicitly stated in the papers, but we assume so, as they lacked measures for high workload traffic situation.

Reasonable weights were given by P. Averty and M. Tognoni.

As for now, these parameters are set equal to 0.5, but are meant to be adjusted and possibly vary with DaC to reflect the results of ongoing research[START_REF] Averty | Work in progress, presented at an internal SDER seminar[END_REF].

APPENDIX: COMPLEXITY METRICS Delahaye and Puechmorel metrics

To present the geometrical indicators introduced in [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], we need to define several quantities:

• The vector representing the distance between two aircraft is denoted by ---→ XiXj where Xi (resp. Xj ) stands for the location of aircraft i (resp. j).

• The "oblical" distance between two aircraft (i and j) is denoted by

where < ., . > stands for the appropriate scalar product. • We denote by -→ v ij = -→ v j --→ v i the speed difference between two aircraft.

• The derivative of the "oblical" distance between two aircraft is denoted by vij and writes

• We introduce a weighting function f . As suggested in [START_REF] Chatton | Etudes de nouvelles métriques de complexité de la circulation aérienne[END_REF], we used

with α = 0.002, β = 0.01, the d ob being expressed in nautical miles. These indicators are defined pointwise. To get a value on the controlled airspace, they have to be averaged on the different aircraft. In [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], a density indicator is defined as follows

Two indicators are introduced to reflect the variability in headings (track_disorder) and speed (speed_disorder). There are defined as

Indicators Div et Conv respectively describe convergency and divergency of the aircraft in the controlled sector.

Indicators Sd+ and Sdare designed to set a weight on potential conflicts that are difficult to solve. These "sensitivity" indicators are defined by

Note that components of the gradient are weighted so as to reflect the difficulty of the respective manoeuvers 5 . As observed in [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], a situation with a high "sensitivity" is easier to resolve for the air controller than one with a low "sensitivity". As these indicators "increase" with the


Charles-Edmond Bichot

Nicolas Durand Log

A tool to design Functional Airspace Blocks

Keywords: functional airspace blocks, partitioning, metaheuristics

This paper focuses on the selection of a technical tool for the establishment of functional airspace blocks in Europe. This paper shows that the creation of functional airspace blocks is a partitioning problem. Some state-ofthe-art partitioning libraries and two metaheuristics are applied to this problem. Comparisons have been made between these libraries and metaheuristics. Results show that the Fusion Fission metaheuristic performs better. The purpose of this paper is not to give the best partition of Europe into functional airspace blocks. It only presents a study which compares and suggests different tools for the establishment of functional airspace blocks in Europe.

Introduction

In the context of the Single European Sky Regulations of March 2004 [549], the European Commission has issued a mandate to the Eurocontrol agency: "Mandate on Support for Establishment of Functional Airspace Blocks (FABs)". The main idea of the Single European Sky (SES) Regulation about the establishment of FABs is to enforce regional cooperation, with the intention of enabling air traffic control to operate efficiently. The SES Regulation calls upon States to reconfigure upper airspace (upper flight level 285) into FABs.

There is no easy definition of what a functional airspace block is. Indeed, the SES Regulation gives some indications by defining what makes a functional airspace block in [549]. It first defines the airspace block as "an airspace of defined dimensions, in space and time, within which air navigation services are provided". Then it defines a functional airspace block as an airspace block which needs to fulfil three overriding attributes which are central to its functionality:

• designed on the basis of operationals requirements;

• more integrated management of the airspace;

• free from the constraints of national borders.

The SES Regulation recalls that the objectives of the single European sky are:

• to enhance current safety standards and overall efficiency for general air traffic in Europe;

• to optimize capacity;

• to minimize delays.

Furthermore, FABs shall respect seven criteria [551]:

1. be supported by a safety case;

2. enable optimum use of airspace, taking into account air traffic flows;

3. be justified by their overall added value, including optimal use of technical and human resources, on the basis of cost-benefit analysis;

4. ensure a fluent and flexible transfer of responsibility for air traffic control between air traffic services unit;

5. ensure compatibility between the configuration of upper and lower airspace;

6. comply with conditions stemming from regional agreements concluded within the ICAO, and 7. respect regional agreements.

Within the framework of the Eurocontrol mandate, a stakeholder consultation has been made, to support the establishment of FABs. The stakeholder consultation indicates a level of unease over the mechanism for assessing each FAB against the seven criteria presented above [Eur05].

The issue of reconfiguring airspace into FABs is addressed in the airspace part of the SES Regulation, and it is clear form the definition that the design of airspace is a key issue in the development of FABs. However, it is at the same time widely recognized that FABs are much more than an airspace design issue. The establishment of FABs, like any regional cooperation in ATM, is a complex process [START_REF] Hallgren | Restructuring european airspace: functional airspace blocks. Skyway[END_REF].

This article focuses on how to optimize functional airspace blocks design. It presents a powerful tool for searching FABs configurations. This tool is based on a new method, called Fusion Fission, which is described in [START_REF] Bichot | Optimisation par fusion et fission. application au problème du découpage aérien européen[END_REF]. This paper is organized as follows: Section 2 describes which problem this article focuses on and tries to resolve. The functional airspace block optimization process is presented more formerly in section 3. Section 4 presents some state-of-the-art partitioning libraries which are adapted to the FAB problem. The Fusion Fission method is shortly described in section 5. Section 6 shows a comparison between the Fusion Fission tools results and other state-of-the-art softwares. After the conclusion, the Fusion Fission FABs result is drawn in the annex.

Optimizing FABs

The study was done at the LOG1 laboratory on the Functional Airspace Block Optimized Process2 . The FABOP project is a study in the "strategic" level of air traffic control, i.e.. methods are applied long before any tactical control of aircraft. The FABOP project respects the definition of what a functional airspace block is. Thus national borders are not taken into account, when operational requirements and management of the airspace are taken into account.

The FABOP project aim is safety. To increase safety in Air Traffic Control, a solution consists in decreasing the workload of air traffic controllers. The air traffic controller workload depends on three factors: traffic supervision, conflicts resolution and coordination [START_REF] Alliot | Techniques d'optimisation stochastique appliquées aux problèmes du contrôle aérien[END_REF]. The goal of FABOP project is to minimize the coordination load between controllers. To minimize the coordination, an algorithm of Fusion Fission is used. Thus, by minimizing the coordination load, the safety is increased.

Each air traffic controller supervises a limited space, called an air traffic sector. Controllers have qualifications to work only on a set of sectors, which is a functional airspace block. The FABOP project consists in cutting the European airspace into FABs. Currently blocks almost never cross countries borders. A new organization of blocks only based on flows of aircraft and not on borders is studied. Because "it is well known that controllercontroller coordination is easier and more effective inside an air traffic control unit (a block) than between air traffic control units" [START_REF] Hallgren | Restructuring european airspace: functional airspace blocks. Skyway[END_REF], the FABOP project searches to maximize flows of aircraft inside FABs and further to minimize flows of aircraft between FABs.

Functional airspace blocks shall respect seven criteria (see section 1). Some criteria are technical, some others are political. The FABOP project only tries to solve technical criteria. Because the FABOP project aims at enforcing safety, the first criterion is respected. The FABOP project works on European flight plans. These data are converted into flow of aircraft between sectors. Thus the second criterion is respected too. The FABOP project includes the cost-benefits analysis of [Eur06] for the sizes of FABs, which is the third criterion. Because the Fusion Fission algorithm minimizes the coordination, the fourth criterion is respected. The lower airspace is taken into account by the FABOP project, which respect the fifth criterion. And the last two criteria are clearly political.

The study area presented in this paper is limited to the air traffic above eleven states: Germany, France, United Kingdom, Switzerland, Belgium, Netherlands, Austria, Spain, Denmark, Luxembourg and Italy. This area is called the country core area in [START_REF] Bichot | A theoretical approach to defining the european core area[END_REF].

A recommendation of the European organization for the safety of air navigation, Eurocontrol, is to create FABs of around 25 sectors [Eur06]. Let Card opt be this number of sectors, the optimal cardinal of the number of sectors per FAB. A second recommendation is to create FABs with a maximum of 50 sectors. Then, because there are 640 sectors in the study area, approximately 26 blocks have to be defined.

The Air Traffic Control partitioning problem

The SES Regulation idea is to partition the European airspace into FABs. Thus the FABOP project uses partitioning methods to design the FABs.

The airspace composed of FABs, sectors, and flows of aircraft can be viewed as a graph. The graph vertices are air traffic sectors and the edges are flows of aircraft between sectors. The weight of a vertex is the number of aircraft which are going through the corresponding sector in a day. And the weight of an edge is the num-ber of aircraft which are going from one sector of the edge to another. The weight of a vertex can be different than the sum of the weights of its connected edges, because arrival and departure of aircraft. Arrivals and departures could not be considered as a sector, because its size would be too large. Some aircraft are going out of the states of the study area. These aircraft are also counted in sectors weights. Thus, the FABOP problem is a graph partitioning problem.

Therefore, the creation of FABs into the European airspace is a graph partitioning problem. Classically, the graph partitioning problem consists in dividing the vertices of a graph into several parts such that the number of vertices in each part is roughly the same and the number of edges-cut is minimized. Because FABs can have up to 50 sectors, with an average of 25 sectors, the number of vertices in each part can be very different. The ratio of the maximum number of sectors of a part to that of the average number of sectors is called the imbalance of the partition. Then, the vertex's number imbalance is

Imb v = 50 Cardopt = 2. To conclude, the FABOP prob- lem is a relaxed graph partitioning problem into k = 26 parts with vertex's number imbalance Imb v = 2.
The graph partitioning problem can be formerly defined as follows: Given a graph G = (V, E), where V is the set of vertices and E is the set of edges.

∀v ∈ V, w(v) is the weight of the vertex v, W (V ) = u∈V w(v) and ∀e = (u, v) ∈ E, w(u, v) is the weight of the edge e. Find a partition π k of k disconnected sub- sets V 1 , . . . , V k of V such that their union is V , ∀i ∈ {1, . . . , k}, Card(V i ) ≤ Imb v * Card opt , W (V i ) ≤ Imb w * W V k (Card(V i )
is the number of vertices of the set V i ) and the sum of the weights of the edges crossing borders parts is minimized. The sum of the weight of the edges crossing border parts is returned by the cut function:

Cut(π) = i,j,i =j Cut(V i , V j )
where,

Cut(V i , V j ) = u∈Vi,v∈Vj w(u, v)
Because the FABs constraints consists in minimizing aircraft flows between FABs and maximizing aircraft flows in FABs, the N cut objective function is better designed than the Cut objective function. Indeed, the N cut objective function [START_REF] Shi | Normalized cuts and image segmentation[END_REF] minimizes the sum of the weights of the edges between parts while maximizing the sum of the weight of vertices of each part.

N cut(π) = i Cut(V i , V -V i ) W (V i )
Because there is no recommendation from Eurocontrol about the weights (i.e. the number of aircraft ) of the blocks in [Eur06], this paper assumes that the imbalance of the weight must not be larger than the imbalance of the number of vertices. Thus, the imbalance of the weight, Imb w , must be lower than 2.

To summarize, the creation of FABs in Europe is a graph partitioning problem. The number of FABs is roughly 26 regarding the paper study area, and each FAB must have less than 50 sectors.

Presentation of some state-of-theart partitioning libraries

Some state-of-the-art graph partitioning packages have been searched to solve the FABs graph partitioning problem.

The METIS library

The METIS library [START_REF] Karypis | METIS : A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices[END_REF] of George Karypis and Vipin Kumar from the department of computer science and engineering at the university of Minnesota, USA, is used. The METIS library is the most popular graph partitioning library. It is very often used to make some comparisons between graph partitioning problems. This library is based on the Multilevel method. It provides two graph partitioning applications, pmetis and kmetis. Their objective is to partition a graph into k equalsize parts while minimizing the Cut objective function. pmetis uses a recursive bisection, and kmetis uses a multi-way algorithm. This library is not directly designed for the FABs problem, because it minimizes the Cut objective function and not the N cut objective function and moreover because the imbalance of the partition is not an input of the algorithm.

The GRACLUS library

The GRACLUS library of Inderjit Dhillon, Yuqiang Guan and Brian Kullis [START_REF] Dhillon | Kernel k-means, spectral clustering, and normalized cuts[END_REF] from the department of computer sciences at the university of Texas, USA, is based on the METIS library and a kernel k-means algorithm.

The graclus application is an adaptation of the METIS library. It can be set to minimize the N cut objective function. Thus the graclus application is well designed to the FABs problem. After some observations, to obtain results as good as possible, the number of local search step is set to 4, while default is 0.

The JOSTLE library

The serial JOSTLE library of Chris Walshaw [START_REF] Walshaw | The serial JOSTLE library user guide[END_REF] from the university of Greenwich, UK, is another state of the art partitioning package. It is probably, at the time of writing (January 2007), the best graph partitioning software which minimizes the Cut objective function.

The jostle software has a customizable balance tolerance. This tolerance range is from 1. to 1.5. Thus, the JOSTLE library is not perfectly designed for the FABs problem, but is better designed than the METIS library. The lowest N cut value of the graph returned by the jostle software has been searched. There is 51 possible partitions into 26 parts corresponding to the different imbalance parameter values. Results show that the imbalance parameter for which the N cut objective function returns the lowest value is 1.37.

The SCOTCH library

The SCOTCH library of François Pellegrini [START_REF] Pellegrini | Schotch and LibScotch[END_REF] from the LaBRI (Laboratoire Bordelais de Recherche en Informatique) at the university of Bordeaux, France, gives a complex software, but very powerful. This is a high customizable library. The partitioning program is based on a dual bisection mapping algorithm, which uses a Multilevel algorithm as a bisection strategy. Partitions can be refined with different methods: Fiduccia-Mattheyses, Gibbs-Poole-Stockmeyer or greedy graph growing methods. The algorithm searches to decrease the Cut objective function, with a customizable imbalance of the weight of the parts.

It occurs that the defaults parameters give good results. The selection operator was not used in the strategy of the dual recursive bisection. The imbalance ratio of the Fiduccia-Mattheyeses method was changed, which is originally set for a perfectly balanced partition.

Simulated Annealing

Metaheuristics are useful tools to resolve problems for which specific algorithms have not been designed, or are not powerful. One of the earliest but always very powerful metaheuristic, is the Simulated Annealing. The algorithm of Simulated Annealing of this paper is based on the article of David Johnson, Cecilia Aragon, Lyle Mc-Geoch and Catherine Schevon [START_REF] Johnson | Optimization by simulated annealing: an experimental evaluation; part i, graph partitioning[END_REF]. The authors present an application of Simulated Annealing on the unweighted bisection graph partitioning problem. Thus, to resolve the FABs partitioning problem, some details were changed in this algorithm.

Because Multilevel algorithms used in the other packages are extremely fast, less than one second of computation, the computation time of the Simulated Annealing algorithm has to be limited. About 2 minutes of compu-tations are allowed to the algorithm. It is not a to long time, but allows relatively good results.

The Fusion Fission method

Because lots of the state-of-the-art partitioning packages where not well designed for the FABs partitioning problem, a new method was proposed by the LOG (laboratoire d'optimisation globale) laboratory at the DGAC, DSNA/DTI/SDER and ENAC, France. A graph partitioning method which can have large imbalance, and which can be used with the N cut objective function has been created.

The Fusion Fission method is a metaheuristic, and as many metaheuristics, it comes from the real life. The Fusion Fission method comes from the nuclear process. The nuclear process can be viewed as a matter reorganization in an optimization process which tends to create atoms of great internal cohesion. In the nature, the atom with the greatest internal cohesion is the iron atom F e, with 56 nucleons in a range of 2 to 235 for other atoms. We can imagine that a reorganization of the nucleons of the atoms naturally tends to create iron atoms, if the number of nucleons and the sort of nucleons allows it.

The graph partitioning problem has some similarities with the nuclear process. The objective of the graph partitioning problem is to find a low energy organization of the parts of a graph. Thus, it is easy to compare atoms with parts and nucleons with vertices. Then, like in the nature, the process to find an organization of low energy consists in merging or splitting atoms. So, parts of the partition are successively merged or split by the Fusion Fission process.

The Fusion Fission method consists in repeating a perturbation process of a partition till a stop condition is reached. The perturbation process starts with a partition of the graph and successively applies fusion or fission to the partition.

A nuclear reaction can only be created in a very high temperature plasma. Thus the Fusion Fission method includes a "temperature" which is used to channel the process. The temperature is used to stop the perturbation process after a definite number of pass.

In the nature, when a fusion or a fission is done, some nucleons are ejected by the process. These nucleons join another atoms, or make other fissions if they have a very high energy. The number of ejected nucleons is know with a certain probability, and the "rules" of ejected nucleons for the reactions can be written. Then some "rules" are defined to eject randomly in a fusion or in a fission some vertices of its part. The process of fusion consists just in merging two parts. The new part is the contraction of the two former parts minus the random number of vertices "ejected" by the rule. Then, if it is possible, these vertices are merged with a different part with which they are the most connected.

The process of fission is much more complex. To split a part into two new parts, a bisection method must be used. There is a great number of bisection methods. An agglomerative method based on the percolation process was chosen. But before splitting the part, like in the fusion process, some vertices can be "ejected" of the part. Then these vertices can be merged with a partition with which they have big connection, or they can split other parts, in a chained fission reaction.

The Fusion Fission algorithm is more precisely presented in [START_REF] Bichot | Optimisation par fusion et fission. application au problème du découpage aérien européen[END_REF].

Comparisons with state-of-theart partitioning packages

One day of European traffic can be converted to a graph as it is described in section 3. Such a graph can be built for different days of data. But in this paper only results which are obtain for the Friday, 17 th June 2005 were presented. The graph is stored using the CHACO graph input file format describes in [START_REF] Hendrickson | The Chaco User's Guide[END_REF]. This format has been chosen for compatibility reasons. The METIS, JOSTLE, GRACLUS and SCOTCH libraries both used this graph input file format. The aim of this section is to find the graph partitioning package which finds the best European FABs partition. During the graph partitioning process, a method can return a particularly good or bad partition, depend- 3 The connection between a vertex and its part is the sum of the weight of the edges connecting this vertex to vertices of the same part ing on the input configuration of the graph. To avoid this kind of problem, benchmark were made with one hundred permutation of the initial graph, and then, average and range of the results were returned. A permutation only consists in renaming each vertex of the graph. Indeed,the structure of the graph is unchanged.

All results presented in this paper are found with an Intel Pentium 4, 3GHz processor with 1Go RAM computer, which uses a GNU/Linux Debian operating system. Every specific graph partitioning package partitions the graph in less than one second, but metaheuristics are greatly slower with an average of 121 seconds for Fusion Fission and 122 seconds for Simulated Annealing. The computational times of the metaheuristics was limited. They can find better solutions with more time. We have found that a roughly 2 minutes computation time is enough but not too long regarding results found.

The computational results are shown in table 1. Algorithms were sorted by their N cut values. First, it can be noticed that metaheuristics have better results than the other algorithms. But the difference of computation time between metaheuristics and the other algorithms can explain these differences. The Fusion Fission algorithm has the best N cut average. It is just followed by the Simulated Annealing algorithm with an N cut average which is 3.2% worth than the Fusion Fission. Then come the Multilevel methods. The scotch algorithm is the first of them. Its N cut average is 4.3% larger than for the Fusion Fission algorithm. Then comes the graclus algorithm. The three state of the art graph partitioning packages which are not design to solve this relaxed k-way graph partitioning problem come at the end. The jostle algorithm has a better N cut average than the METIS package, because it takes into account algorithm has an imbalance parameter.

The vertex and weight imbalance of the partitions found for each algorithm can be shown in figure 1. The horizontal line represents the threshold of Imb v ≤ 2 or Imb w ≤ 2. The Simulated Annealing vertex average is really above this threshold, and it seems that it has too much weight imbalance above this threshold too. As it was expected, the jostle, pmetis and kmetis algorithms have very low weight imbalance compared to the thresh- As it is explained in subsection 4, this is due to their very strict weight partitioning constraints. However, both of their vertex imbalance have a too big range, with an upper range of approximately 2.5 ≫ 2. Partitions found by Fusion Fission and graclus have good vertex imbalance, both on average and range. Thus, the Fusion Fission algorithm most respects the FABs partitioning problem constraints.

The Fusion Fission algorithm not only found partitions for a fixed number of parts, but for a range of values around this number. Then, when the Fusion Fission algorithm resolves the FABs partitioning problem, it finds partitions with a partition cardinal range: {19, . . . 30}. And, for more than 90% of the graph permutations, a partition cardinal range: {18, . . . 31}. To compare all of the partitions found by the Fusion Fission algorithm with other partitions, the two better Multilevel algorithms, scotch and graclus, computed the one hundred permutations of the air traffic flow graph for a number of parts btween 18 to 31. Results of the computation of the Fusion Fission algorithm with a number of part, 26 and for the 13 computations of the scotch algorithm and for the 13 computations of the graclus algorithm are displayed figure 2. In this figure, each algorithm results are presented as a ratio of the N cut average to that of the N cut average of graclus. The partitions found by the Fusion Fission algorithm are better than those of scotch for a number of parts btween 22 to 28. And the partitions found by the Fusion Fission algorithm outperforms the partitions found by graclus for a number of parts between 21 to 29.

Conclusion

This paper focuses on technicals tools for the establishment of functional airspace blocks in Europe. This paper shows that the creation of functional airspace blocks is a partitioning problem. The state-of-the-art partitioning libraries METIS, GRACLUS, JOSTLE and SCOTCH have been used to solve this problem. Because these libraries where not perfectly designed to solve this problem, due to imbalance constraints, two metaheuristics have been used. The first is a simulated annealing method applied to graph partitioning. The second is a Fusion Fission method applied to graph partitioning. Comparisons have been made between these methods. The result is that the Fusion Fission performs better than the other methods.

The purpose of this paper is not to give the best partition of the Europe into functional airspace blocks. It only presents a study which compares and suggests different tools for the establishment of functional airspace blocks in Europe. Figures that can be shown in the annex never attempts to be the perfect solution for the establishment of functional airspace blocks in Europe. They just show what may be the single European sky with functional airspace blocks. Because some of the seven criteria presented in the introduction are not respected or not fully respected, tools presented in this paper must be enhanced to respect all of this criteria. Annex: one of the results found

The European airspace area presented in this paper is actually partitioned into 55 blocks. 4 is just a solution, it is not the perfect solution and do not attempt to be so. Moreover the criterion used does not take into account all operational requirements and this partition would be operationaly unfeasible. A lot of work remains to be done to apply this algorithm to real world problems.

Nicolas Durand graduated from the École Polytechnique de Paris in 1990 and from the École Nationale de l'Aviation Civile (ENAC) in 1992. He has been a design engineer at the Centre d'Études de la Navigation Aérienne (now DSNA/DTI/SDER) since 1992 and holds a Ph.D. in computer Science (1996). He is currently head of the he Planning, Optimization, and Modelling team of the DSNA/DTI R&D domain.

Figure 1 :

 1 Figure 1: Algorithms average and range results of the vertex and weight imbalance of 100 permutations of the air traffic flow graph

Figure 2 :

 2 Figure 2: Comparisons between Fusion Fission results for k = 26, with GRACLUS and SCOTCH results for k ∈ {18, . . . , 31}. Each algorithm results are presented as ratio of the N cut average to that of the N cut average of graclus.

Figure 3 :

 3 Figure 3: The real European airspace partitioning at flight levels 180 and 280

Figure 4 :

 4 Figure 4: The best European airspace partition found with the Fusion Fission algorithm at flight levels 180 and 280

Table 1 :

 1 Average of the results of 100 permutations of the air traffic flow graph automatically adjusted during the execution of the algorithm by a self-learning function. A rule is just a probability function to eject a certain amount of vertices of the atom. Vertices which have a proportionally low connection 3 with the other vertices of the same part are ejected first.

	These rules are

Table 2 :

 2 The European airspace partition and the best Fusion Fission partition

	Partitions	Number of part N cut w imb v imb
	The 6/17/2005 partition	55	25.10 4.64 3.33
	The best FF partition	26	4.87	1.69 1.64

  Table 2 presents a comparison between the European airspace partition of June, 17 th 2005 and the best partition found by the Fusion Fission algorithm for 26 parts. Because the numbers of parts of the two partitions are very different, the N cut comparison is not relevant. The imbalance of the vertices number of the current European airspace partition into centers of control is really above the Eurocontrol recommendation. The European airspace partition of June, 17 th 2005 at flight level 180 and 280 can be shown figure 3. The best Fusion Fission partition found is presented by figure 4, it takes four hours of computation. The European partition into FABs presented by figure

Laboratoire d'Optimisation Globale, www.recherce.enac.fr

FABOP, www.recherche.enac.fr/opti/FABOP

Author biographies

Charles-Edmond Bichot graduated from the École Nationale de l'Aviation Civile (ENAC) in 2004. He is a Ph.D. student in computer Science at the LOG laboratory of the DSNA/DTI/SDER and of the ENAC.