Daniel Delahaye
email: delahaye@recherche.enac.fr

Stéphane Puechmorel
email: puechmor@recherche.enac.fr

Aircraft Local Wind Estimation From Radar Tracker Data

Keywords: Wind estimation, Speed triangle, Trajectory prediction, Kalman lter, Air traf c management

Accurate wind magnitude and direction estimation is essential for aircraft trajectory prediction. For instance, based on these data, one may compute entry and exit times from a sector or detect potential con ict between aircraft. Since the ight path has to be computed and updated on real time for such applications, wind information has to be available in real time too.The wind data which are currently available through meteorological service broadcast suffer from small measurement rate with respect to location and time. In this paper, a new wind estimation method based on radar track measures is proposed. When on board true air speed measures are available, a linear model is developed for which a Kalman lter is used to produce high quality wind estimate. When only aircraft position measures are available, an observability analysis shows that wind may be estimated only if trajectories have one or two turns depending of the number of aircraft located in a given area. Based on this observability conditions, closed forms of the wind has been developed for the one and two aircraft cases. By this mean, each aircraft can be seen as a wind sensor when it is turning. After performing evaluations in realistic frameworks, our approach is able to estimate the wind vectors accurately.

I. INTRODUCTION

When an air traf c controller observes its traf co nar a d a r screen, he tries to identify convergent aircraft which may be in con ict in a near future, in order to apply maneuvers that will separate them. The problem is to estimate where the aircraft will be located in this near future (5-10 minutes); this process is call trajectory prediction. The trajectory prediction depends mainly on the residual noise after ltering (see [START_REF]Con ict detection and resolution in the presence of prediction error[END_REF]): the weight of the aircraft, the temperature and the wind. The residual noise is integrated with time with a growing covariance matrix indicating that the estimated position is less and less accurate. The weight of the aircraft is relevant in the ight dynamic model but is still a raw data. The engines of aircraft are sensitive to the air temperature and such a data is very useful to model the trust of the aircraft but it is also very dif cult to measure on real time. Finally, the wind in uences strongly the cinematic of the aircraft and limits also the trajectory prediction. Based on the available accuracy, the actual limit of the trajectory prediction is about 20 minutes. It means that after 20 minutes the uncertainty is so big that the estimated position is no more useful for any ATM purposes. Several efforts have been tried to improved the trajectory prediction by estimating or suppressing the wind [START_REF] Mondoloni | Improving trajectory forecasting through adaptive ltering technique[END_REF], [START_REF]Wind prediction accuracy for air traf c management decision support tools[END_REF], [START_REF] Delahaye | Wind eld update using radar track data[END_REF]. The present paper, proposes a new method for estimating the wind around aircraft by the mean of observations of the radar tracks and some down linked data. One goal of this work is to show how it is possible to extract wind information from the radar observations. The paper is organized as follow : the rst part presents the relation between air speed vector, ground speed vector and the wind. The second part gives the observability condition of the wind based on the available measures. The third part presents the different models which have been used for the wind extraction. Finally, the fourth part presents some results and compares the performance of our models.

II. GENERAL RELATIONS BETWEEN SPEED VECTORS

The following notations will be used in the paper. Vectors and matrices are shown with underlined symbols. Speeds

V = v x v y V = v 2 x + v 2 y θ g =arctan v x v y
Where θ g is the The route angle (with reference to north).

The same notation will be used for the air speed (T ,t x ,t y ,θ a ) and the wind (W ,w x ,w y ,θ w ). We have :V = T + W This relation is shown on gure 1.

Air Turning Rate ω a =

with γ x = dvx dt and γ y = dvy dt . It can be shown that the air turning rate (ω a ) and the ground turning rate (ω g ) are related by the following expression :

ω g = T 2 +T * W * cos(θa(t)-θg(t)) T 2 +W 2 +2.T.W cos(θa(t)-θg(t)) ω a (2)

A. Hypothesis

We consider en-route traf c( t r a f c away from airports), with aircraft ying in cruise phase. Aircraft are supposed to y at constant air speed (T = C te ) and turn with constant air turning rate (ω a = C te ).

The average wind is supposed to be constant in the neighborhood of the aircraft.

III. OBSERVABILITY CONDITIONS

A. Mode S Radar

Modee S radars [START_REF] Lefas | Improved tracking with mode-s data-linked velocity measurements[END_REF] are able to access on board parameters and especially airspeed vector. Having such measures the system is fully determined at any time as it is shown in the following equations.

v x = Tsin(θ a )+w x v y = Tcos(θ a )+w y (3)
In this system, there are two equations and two unknowns (w x ,w y ); the other parameters are given by the radar (v x ,v y ,T,θ a ). AK a l m a n lter [START_REF] Grewal | Kalman Filtering Theory and Practice[END_REF], [START_REF] Welch | An introduction to the kalman lter[END_REF] will be used to extract the wind by removing the noise.

B. Classical Radar

If measures come from classical radars, only position measures are available (x, y) and two situations have to be taken into account.

1) First situation: One aircraft: In the case where radar measures come from one aircraft only, wind may be observed only after two asymmetric turns (meaning three straight lines separated by two turns). As a mater of fact, the rst segment brings the following system of equations ;

v 1x = Tsin(θ a1 )+w x v 1y = Tcos(θ a1 )+w y (4)
In this system (4), there are 4 unknowns (T,θ a1 ,w x ,w y )a n d two equations, so two equations are missing.

After the rst turn two new equations are added to this system with one extra unknown (θ a2 ):

v 1x = Tsin(θ a1 )+w x v 1y = Tcos(θ a1 )+w y v 2x = Tsin(θ a2 )+w x v 2y = Tcos(θ a2 )+w y (5)
In this new system (5), there are 5 unknowns (T,θ a1 ,θ a2 ,w x ,w y ) and only 4 equations. It is only after the second turn that the system is fully determined (6 unknowns (T,θ a1 ,θ a2 ,θ a3 ,w x ,w y ) and 6 equations) :

⎧ ⎨ ⎩ v 1x = Tsin(θ a1 )+w x v 1y = Tcos(θ a1 )+w y v 2x = Tsin(θ a2 )+w x v 2y = Tcos(θ a2 )+w y v 3x = Tsin(θ a3 )+w x v 3y = Tcos(θ a3 )+w y (6)
This system has a closed form solution for which the wind is given by :

w x = (v3y-v2y )V 2 1 +(v1y -v3y)V 2 2 +(v2y-v1y )V 2 3 2{v1y (v2x-v3x)+v2y (v3x-v1x)+v3y (v1x-v2x)} (7)
w y = (v2x-v3x)V 2 1 +(v3x-v1x)V 2 2 +(v1x-v2x)V 2 3 2{v1y(v2x-v3x)+v2y(v3x-v1x)+v3y(v1x-v2x)} (8)
with

V 1 = v 2 1x + v 2 1y V 2 = v 2 2x + v 2 2y V 3 = v 2 3x + v 2 3y
For both expressions (w x ,w y ) the denominators must not be equal to zero meaning that turns have to be asymmetric

(V 1 = V 2 = V 3 )
Having expressions for w x and w y it is very easy to extract the other unknowns :

⎧ ⎨ ⎩ t 1x = v 1x -w x t 1y = v 1y -w y t 2x = v 2x -w x t 2y = v 2y -w y t 3x = v 3x -w x t 3y = v 3y -w y
2) Second situation : Two aircraft: When radar measures are available for two aircraft (a and b), only one turn for both trajectories is needed to have enough information for wind estimation. The two rst segments bring the following system of equations with 6 unknowns (T a ,T b ,θ aa1 ,θ a b1 ,w x ,w y )a n d 4 equations : va 1x = T a sin(θ aa1 )+w x va 1y = T a cos(θ aa1 )+w y vb 1x = T b sin(θ a b1 )+w x vb 1y = T b cos(θ a b1 )+w y (9) After the rst turn, the new systems is fully determined with 8 unknowns(T a ,T b ,θ aa1 ,θ a b1 ,θ aa2 ,θ a b2 ,w x ,w y )a n d8 equations :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ va 1x = T a sin(θ aa1 )+w x va 1y = T a cos(θ aa1 )+w y vb 1x = T b sin(θ a b1 )+w x vb 1y = T b cos(θ a b1 )+w y va 2x = T a sin(θ aa2 )+w x va 2y = T a cos(θ aa2 )+w y vb 2x = T b sin(θ a b2 )+w x vb 2y = T b cos(θ a b2 )+w y (10)
The associated closed form of the wind is given by :

w x = (vb1y -vb2y)(Va 2 1 -Va 2 2 )+(va2y-va1y )(Vb 2 1 -Vb 2 2 ) 2{(va1x-va2x)(vb1y-vb2y )-(va1y-va2y )(vb1x-vb2x)} (11) w y = (vb2x-vb1x)(Va 2 1 -Va 2 2 )+(va1x-va2x)(Vb 2 1 -Vb 2 2 ) 2{(va1x-va2x)(vb1y-vb2y )-(va1y-va2y )(vb1x-vb2x)} (12) with Va 1 = va 2 1x + va 2 1y Va 2 = va 2 2x + va 2 2y Vb 1 = vb 2 1x + vb 2 1y Vb 2 = vb 2 2x + vb 2 2y and (V a 1 = V a 2 ;V b 1 = V b 2 )
Like for the one aircraft case, it is very easy to extract the other unknowns :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ta 1x = va 1x -w x ta 1y = va 1y -w y tb 1x = vb 1x -w x tb 1y = vb 1y -w y ta 2x = va 2x -w x ta 2y = va 2y -w y tb 2x = vb 2x -w x tb 2y = vb 2y -w y
Based on those observability conditions, several wind estimation models have been developed which are now in the following section.

IV. MODELS

A. Model with Air Speed and Air Turning Rate (model 1)

The state vector used in the Kalman lter is given by :

X(k)=[x(k) y(k) t x (k) t y (k) w x (k) w y (k)] T (13)
where x(k),y(k) is the position, t x (k),t y (k) t h eT r u eA i r Speed (TAS) and w x (k),w y (k) the wind.

The measure vector consists in the radar position and the true air speed :

Z(k)=[x m (k) y m (k) t xm (k) t ym (k)] T (14)
Having access to the air turning rate (ω a ), it can be included in the prediction matrix. The structure of the system is the following :

X(k +1)=F (k).X(k)+v(k) (15)
where

F (k)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 10 C 1 (ω a ) C 2 (ω a )∆ t 0 01-C 2 (ω a ) C 1 (ω a )0∆ t 00 C 3 (ω a ) C 4 (ω a )0 0 00-C 4 (ω a ) C 3 (ω a )0 0 00 0 0 1 0 00 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (16) Z(k)= ⎡ ⎢ ⎢ ⎣ 100000 010000 001000 000100 ⎤ ⎥ ⎥ ⎦ .X(k)+w(k) (17)
where

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ C 1 (ω a (k)) = sin(ωa(k)∆t) ωa(k) C 2 (ω a (k)) = 1-cos(ωa(k)∆t) ωa(k) C 3 (ω a (k)) = cos(ω a (k)∆ t ) C 4 (ω a (k)) = sin(ω a (k)∆ t ) (18)
The model being exact, the model noise covariance matrix R =0. The measure noise covariance matrix R is given by :

R = ⎡ ⎢ ⎢ ⎣ σ 2 p 000 0 σ 2 p 00 00σ 2 T 0 000σ 2 T ⎤ ⎥ ⎥ ⎦
with σ p = 100 meters and σ T =0.2 kts.

B. Model with Air Speed only (model 2)

For this model, the state vector and the measure vector are the same as in the rst model.

T h ep r e d i c t i o nm a t r i xi sn o wg i v e nb y:

F (k)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 10∆ t 0∆ t 0 0 10∆ t 0∆ t 00 1 0 0 0 00 0 1 0 0 00 0 0 1 0 00 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (19)
The measure matrix H and the associated covariance matrix R are the same as in the rst model. This model is linear but is false for the air speed vector evolution. In order to take into account this model error, the following covariance matrix (Q) is included in the lter : Q is given by :

Q = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
00 0 0 00 00 0 0 00 001 .000 0 0 001 .000 00 0 0 00 00 0 0 00

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

C. Model with Turning Rate only (model 3)

This model is the same as the rst one with the following measure vector ;

Z(k)=[x m (k) y m (k)] T (20)
and the associated measure equation :

Z(k)= 100000 010000 .X(k)+w(k) (21)

D. Models Without Command 1) Model using Kalman Filter (model 4):

Having no access to the turning rate, this parameter may be extracted by a Kalman lter (ω a included in the state vector). If the turning rate is included in the state vector, the evolution of the system is not linear and an EKF [START_REF] Ribeiro | Kalman and extended kalman lters : Concept, derivation and properties[END_REF] or an UKF [START_REF] Julier | A new extension of the kalman lter to nonlinear systems[END_REF] has to be used to manage such state vector.

2) Model using turns: As it has been shown in the previous section (observability conditions), wind may be extracted by observing ground radar track during turns.

The key element of this approach is the turning rate detector based on equation [START_REF] Cohn | The behavior of the forcast error covariances for a kalman lter in two dimensions[END_REF]. When ω g is greater than a given threshold, the aircraft is considered to be in turn. Based on this turn detector, straight line segments are easily identi ed for which ground speed vector averages are computed. Speed vector estimates in straight lines are given by the framework of gure 2. In this gure, the turns are rst detected and a counter is then used to select the right averaging process in order to compute ground speed estimates. When two aircraft are used for wind estimation, the same framework is duplicated and counters select only two averaging blocks in order to produce the following estimates :V a1 ,V a2 ,V b1 ,V b2 .

V. R ESULTS

A. Simulation framework

In order to test and compare our models, a radar tracker simulator has been used to produce the reference trajectory which has been disturbed with a Gaussian noise. The disturbed trajectory is then ltered by a Kalman lters which generates the estimated trajectory which is then compared to the reference trajectory. This framework is summarized on gure 3. For all experiments, a wind of 40 kts has been used with θ w = 240 o ⇒ W =[ -34.64kts, -20kts] T = [-17.82m/s, -10.28m/s] T .

B. Models with Mode-S radar data (model 1,2 and 3)

The rst trajectory used for our experiments is built with 3 straight lines (20 minutes for each) connected with turns as it can be seen on gure 4.This trajectory has been disturbed by a Gaussian noises for which the means are zero and the standard deviations are the following :σ position = 100m and σ TAS =0 .2kts. Those values are given by the performance of the actual radar trackers.

1) Results for models 1 and 2: Models 1 and 2 have been tried on this trajectory for which the residual errors on the wind estimates are given on gures 5 and 6. The gure 5 shows the residual error for the wind strength estimate after the convergence phase which last 2 minutes. Figure 6 shows the same kind of results for the residual wind angle estimates. The models 1 and 2 have the same evolution during the convergence phase. After this convergence phase, the wind strength estimate error stays below 0.2kts for both models 5.

The model 1 has a more accurate behavior at the end of the simulation with an error which stay below 0.1kts. The wind angle residual error stays below 1 degree (in absolute value) for both model 6. As for the wind strength estimation, the second Having now some wind estimates on some points in the airspace (where the aircraft are located), the following step consists in the global wind eld interpolation based on the meteorological model called "Shallow-Water". This model is described in [START_REF] Cohn | The behavior of the forcast error covariances for a kalman lter in two dimensions[END_REF] and the associated discret expansion in given in [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF].

VI. CONCLUSION

This paper has given a new approach for extracting the wind information from the radar tracks. Two approaches has been presented. When True Airspeed measures are available, linear models may be used with regular Kalman lter. In this rst approach, wind estimate are available all along the aircraft trajectory. When only position measures are available, an observability analysis has shown that wind may be estimated only after turns (one turn for two aircraft or two turns for one aircraft). For such approach, a closed form of the wind has been developed. Then, those models have been validated with realistic simulations.

Fig. 1 .

 1 Fig. 1. Speed vectors relations

Fig. 2 .

 2 Fig. 2. This gure presents the framework used to estimate the average speed in straight line between turn.

Fig. 4 .Fig. 5 .

 45 Fig. 4. First trajectory. Each straight lines last 20 minutes. The turning rate of the rst turn is 1 deg and -1 deg for the second turn. The aircraft is considered to y on cruise phase (vertical speed equal zero).

Fig. 6 .Fig. 10 .

 610 Fig. 6. Residual wind angle error for model 1 and 2

Fig. 11 .

 11 Fig. 11. Turning rate estimates of both aircraft

2) Results for models 3: When air speed is removed from measures, model 1 is not able to converge and model 3 (which is models 2 without air speed measures) has to wait the second turn to be able to produce an estimate of the wind (see gure ?? and 8). During the convergence phase, model 3 reduces the wind strength error till 40kts which is the limit for such estimation because the Kalman lter is not able to distinguish in the position measures the part coming from the air speed and the part coming from the wind. After the convergence phase (see gure 7) the lter has to wait the second turn in order to be able to produce a reliable estimate of the wind strength (as it has been previously shown (observability conditions), it is only after the second turn that the wind observability conditions are met). As for the wind strength error, the wind angle error reach zero after the second turn. During the convergence phase the wind angle error reach 160 deg (see gure 8).

C. Models with classical radar data (model 4 and 5)

1) Results for models 4: When on board measures (TAS, heading, turning rate) are not available, the turning rate has to be included in the state vector inducing a non linear prediction equation of the lter. In order to address such non linear evolution, an EKF and an UKF has been developed and tested on the rst trajectory. Both lters are not able to estimate the wind vector even after the second turn. This behavior is due to the residual error of the Taylor expansion used in the EKF or in the UKF. Based on the observation (positions) the lter is not able to nd the part due to the wind and the one due to the error coming from the Taylor expansion.

2) Results for model 5: Based on the turning rate estimate given by equation 1, it is very easy to build a turn detector when |ω g | > threshold. The ground vector speed has been disturbed by a Gaussian noise with zero mean and 0.2kts standard deviation (performance of the actual trackers). This method has been used on the rst trajectory. The ground turning rate estimate is given by gure 9. The two turns are well identi ed with a value of about +/-1.3 degree/sec in turns (it must be noticed that the ground turning rate is different from the air turning rate). The results of this approach are given in the following table :   Wind Strength Error Est and Angle w x = -17.6798m/s 39.65kts 0.35kts w y = -10.1831m/s 240.053deg 0.053deg

Those results has been computed by using equations ( 7), [START_REF] Mondoloni | Improving trajectory forecasting through adaptive ltering technique[END_REF] and the following ground speed estimates coming from the averaging process :


Daniel Delahaye

Stéphane Puechmorel

Stephane Puechmorel

3D AIRSPACE DESIGN BY EVOLUTIONARY COMPUTATION

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

Introduction

Air traffic controllers monitor the traffic and check that aircraft follow their planed trajectories. When two or more aircraft are converging to the same point, the controllers have to change the heading, the speed or the altitude of some aircraft in order to ensure a minimum distance between all aircraft. Everyday about 8000 aircraft are flying in the French airspace and such amount of traffic generates a huge control workload which cannot be managed by a single controller (as any human being, a controller has working limits and may usually be able to manage a maximum of 15 aircraft with 3 conflict points); the airspace is then partitioned into sectors each of them being assigned to a controller team. When airspace is partitioned into sectors, several rules have to be applied in order to insure the safety and an effective traffic ow. The workload in sectors has to be balanced in order to ensure that controllers have roughly the same amount of work to do in all sectors. When an aircraft crosses a sector boundary, controllers in charge of those neighboring sectors have to exchange information with pilot and between each other in order to insure a safe transfer of the ight between sectors. This transfer is called a coordination of an aircraft between sectors. The way airspace is partitioned may induce different amount of coordination and a good sectoring has to minimize this number of coordination because it cost efforts to the controllers and to the pilots too. Those objective have to be optimized by taking into account several operational constraints. The rst constraint ensures that sectors are convex in the airway link directions (An airway is materialized by a list of linear links). By this mean, a pilot will not encounter the same controller twice during his ight. This constraint will then prevent super uous coordinations; this means that an aircraft crossing a sector will encounter 2 and only 2 sector boundaries (route convexity constraint). The second constraint will ensure that con ict point will be located at a minimum distance from the sector boundaries. Aircraft have to follow airways and con ict between aircraft may appear only at airways crossing. A sector boundary has to be at least at a given distance from each airways crossing (safety constraint). When a controller has to solve a con ict, he needs a minimum amount of time to develop a solution. Each controller managing individually his sector, if a sector boundary is too close to a crossing point, he is not able to solve any con icts because he has not enough time between the coordination step (with the previous sector where the aircraft comes from) and the time the aircraft reaches the crossing point. The minimum delay time is xed at 7 minutes and can be converted into a distance once the aircraft speed is known. The third constraint will ensure that aircraft stay enough time in sectors for controllers to be able to manage them(min stay time constraint). As a matter of fact, if an aircraft stays a few seconds in a sector, the controller will not be able to manage such aircraft; then it brings only super uous coordination and the controller is not able to produce any services for this kind of aircraft. We express this constraint by a minimum distance between two boundaries cutting the same airway. Finally, a fourth constraint will limit the possible shape of synthesized sectors (geometrical constraint). Controllers work on a two dimensional 978-1-4244-2208-1/08/$25.00 ©2008 IEEE.

3.B.6-1 radar screen with aircraft tag giving the altitude of traf c. This HMI limitation may induce uncertainty on the sector which contains the plane. As a matter of fact, if the sector border are not perpendicular to the ground in the third dimensions, controller may have problem to know if an aircraft is still in a neighbor sector or already in its own sector. This problem is shown in the Figure 1.

Fig. 1. Vertical Border Constraint

In Figure 1, the second situation (right side) is easier to manage from the controller point of view than in the rst one (left part) because of the shape of the border in the side view. As it is shown on this gure, the sector border is not perpendicular to the ground (xy). The controllers having only the top-view of sector, they have dif culties to identify the point of coordination when the aircraft is exiting S1 and entering S2. On the left of line "a" the aircraft is in sector S1, on the right of line "b" the aircraft is in sector S2 but between "a" and "b" it is dif cult to know when the aircraft will change of sector mainly if it is descending or climbing. The way to solve this problem is to ensure that sector boundaries are perpendicular to the ground in the third dimension as it is shown in the right part of Figure 1. The bottom right shows the new structure of the vertical border which is now perpendicular to the ground. The upper right shows the sectors top view with the "c" line materializing the sector transition. When the aircraft crosses this line, the controllers are sure that it is changing of sector even in case of some uncertainty in the altitude. Furthermore, the shape of sectors in the horizontal dimension has to be polygonal ensuring linear borders which are necessary for controllers to be effective. Then, sectors are cylinders with polygonal section.

Based on the problem objective and constraints de nition one can summarized the framework of this study: one considers a traf c distribution in the airspace which generates a control workload density. The problem consists of partitioning this airspace into K sectors, in order to minimize both objectives: [START_REF] Lai Andw | Anevolutionary algorithm for the rectangular cutting stock problem[END_REF] workload balancing between sectors and (2) the number of coordinations ( ow cut by the sector boundaries). This optimization process has to meet four constraints : [START_REF] Lai Andw | Anevolutionary algorithm for the rectangular cutting stock problem[END_REF] This problem can be viewed as a 3D geometrical cutting problem with typical constraints and the method presented in this paper may be applied to others cutting problem in 3D space. The next part gives a summary of the previous related works. The third part gives a description of the complexity metric used to model the workload. The fourth part proposes a mathematical modeling and identi es the associated complexity of this 3D cutting problem. A fth part gives a brief overview of the Evolutionary Algorithms which has been used to address such problem. The sixth part gives some details about the chromosome coding and the associated operators dedicated to this problem. Finally, the last part presents results on realistic instances of this problem.

Previous Related Works

The problem of cutting optimization can be summarized the following way. Having a limited domain in R n one has to cut this domain into K sectors in order to optimize an objective function. This optimization process has to take into account several problem dependent constraints. This problem may nd many application in the industry mainly for the stock cutting optimization for which a factory produces material which has to be cut into different pieces ordered by the clients. The objective consists in minimizing the waste after the cutting. This problem has been studied in one, two or three dimensions and may be extended to 3.B.6-2 higher dimensionality [START_REF] Lai Andw | Anevolutionary algorithm for the rectangular cutting stock problem[END_REF][START_REF] Cui | Generating optimal multisegment cutting patterns for circular blanks in the manufacturing of electric motors[END_REF][3][START_REF] Kampke | Simulated annealing: Use of a new tool in bin packing[END_REF]. The associated complexity increases with the dimension of the space which has to be cut. For stock cutting problems, the shape and the size of the synthesis sectors are xed and controlled by the demand. From the mathematical point of view, this kind of cutting problem can be view as an packing problem [START_REF] Dyckhoff | A typology of cutting and packing problems[END_REF][START_REF] Suterf | Fastevaluationofvector splinesin three dimensions[END_REF][START_REF] Falkenauer | Ahybrid grouping genetic algorithm for bin packing[END_REF] which is a typical NP_Hard assignment problem. Those techniques cannot be applied to our problem because the size and the shape of our sectors are not known in advance. Our objective is also different as we are not looking at minimizing the waste but we try to nd a full partition of the airspace (no part of the airspace is left out of sectors) in order to balance the weight of the synthesis pieces with a minimum cutting section (coordination minimization). Depending of the dimension of the object that has to be partitioned into sectors, the associated complexity may strongly change.

For the 1D space, the object that has to be cut is materialized by a line with a linear weight density ( (x)) (no punctual weights is considered). The total weight of the object is (where L is the total length of the object).For this simple case, a trivial solution is given by a dichotomy search which rst x X 1 then x X 2 and so on till X K (where K is the total number of sectors). When punctual weights are introduced in the problem (see Figure 2), one are not sure to reach an exact solution.

Figure 2. 1D Case where the Strip is Partitioned into 4 Balanced Sub-Strips

The weight of an individual segment k limited by and is given by:

(1)

where is the position of the punctual weight and the associated weight . We build the following criterion:

(2) and use a classical gradient method to nd the minimum of J. When the coordination minimization is consider the criterion J is modi ed the following way:

(3)

where the second part represents the cutting section minimization. The coordination being a workload for the controller, it means that when the boundary is built, the associated sectors have to be updated with this new weights. This is taken into account in the de nition:

where is a weight factor. The way to solve such problem is the same as in the case where the is no punctual weight. The two dimensions case is much more dif cult to solve because of the size of the state space. When the cutting direction is xed a priori the problem is the same as the 1D case for which the linear weight density is the projection of the surface density (x,y) on the cutting direction. Depending of this direction, one can nd several solutions to the problem. When two cutting directions are possible (the sectors are then rectangles), there is no way to translate this problem into one dimension and the associated complexity increases drastically. Having polygonal shape of sectors, the number of cutting direction may be much higher than two and the associated complexity increases accordingly. This problem has been studied in [START_REF] Delahaye | Genetic algorithms partitioning airspace[END_REF][START_REF] Delahaye | Spatial Evolutionary Modelling (chapter 7)[END_REF] for which the sectors are synthesized by the mean of Voronoi diagram and optimized by genetic algorithms. In this approach, points called class centers are randomly spread in the 2D plan that has to be partitioned and the 3.B.6-3 associated Voronoi diagram is built producing a random sectoring of the plane (see Figure 3).

Figure 3. Random Sectoring by Voronoi Diagram

The Voronoi diagrams produce convex sectors ensuring that any line crossing the sector will encounter two and only two borders. Having this property, such sectoring will always met the route convexity constraint. The other constraint are taken into account by relaxation of the objective function. A discrete model of this approach has been proposed in [START_REF] Delahaye | Airspace sectoring by stochastic optimization[END_REF]. The approach consists in using a mosaic of elementary cells (this maybe the result of a Voronoi diagram synthesis) and in building the associated Delaunay triangulation (see Figure 4) in order to create a connected planar graph. This principle makes the constraint management easier. As a matter of fact, when a mosaic cell is too small or has a con ict point too close to one of its boundaries, it is merge with one of its neighbors in order to create a bigger cell which met the constraint. After such a cell checking, all the cells met the constraint and the method consists in merging cells by the mean of connected component optimization in the Delaunay graph [START_REF] Delahaye | Airspace sectoring by stochastic optimization[END_REF]. This optimization is based on an evolutionary process. This discrete method enable the management of bigger airspace and is faster to reach the optimum. The two previous models may be straightly extended in the third dimension. For the rst one, it is only necessary to used 3D position for the class centers. This modi cation will induce the synthesis on 3D polyhedral sectors. For the second one, one has only to consider Voronoi cells in three dimensions and build the associated Delaunay triangulation inducing a graph with vertex in a 3D space. In both cases, such 3D straight extension will result in sectors which will not met the geometrical constraint. These previous works have been mainly done in the 2D case and may not be extended straightly in the third dimension. The next section introduced the problem modeling which has been used in this new 3D approach.

Complexity Metric

Traf c complexity metrics quantify congestion in airspace more accurately than a simple number of aircraft which is independent of the traf c con guration. This metric captures the level of disorder (or organization) of any set of trajectories in a 4D space (3D for the space and 1D for the time) and represents a good indicator of the underlying control workload. The key idea of this approach is to model the set of aircraft trajectories by a dynamical system which is summarized by the following equation: [START_REF] Dyckhoff | A typology of cutting and packing problems[END_REF] where is the state vector of the system (= [ x(t),y(t),z(t)] T ) and V : C 2 vector eld depending on parameters , describe systems which integral curves may t the observed trajectories. This equation associates a vector speed to a position in the space coordinate and then synthesis a particular vector eld.

Based on the observations of the aircraft (positions, speed vectors and times), the dynamical system has to be adjusted with the minimum error. This tting is done with a Least Square Minimization (LMS) method. For each considered aircraft i, it is supposed that position i =[xi,yi,zi] T and speed vector Vi =[vx i ,vy i ,vz i ] T are given (radar tracker data). An error criterion between the dynamical system model and the observation is computed:

3.B.6-4 (6)
where N is number of aircraft. There are many classical ways of obtaining a class of parametrized vector elds which ful ll the tting requirement. Among them, vector splines allow a control on the smoothness of vector elds, which is important in our case since civil aircraft maneuvers are based on low acceleration guidance laws. Computing traf c complexity for a given traf c situation requires interpolating a vector eld given only samples (positions and speeds of aircraft at a given time).Vector spline interpolation seeks the minimum of a functional of the form: [START_REF] Falkenauer | Ahybrid grouping genetic algorithm for bin packing[END_REF] where is a vector eld de ned on a domain D R n , L is an elliptic differential operator controlling smoothness of the solution and ( i, i)i=1...m are the interpolation data. By introducing the adjoint operator L T , optimal vector eld can be shown to be a linear combination of shifted version of the elementary solution kernel of the differential operator L T L.A special case is the so-called "divcurl" splines with the criterion: [START_REF] Delahaye | Genetic algorithms partitioning airspace[END_REF] with , positive reals controlling the smoothness of the approximation by focusing on constant divergence or constant curl. The metric chosen for complexity computation relies on a measure of sensitivity to initial conditions of the underlying dynamical system called Lyapunov exponents. In order to gure out what Lyapunov exponents are, let consider a point and look at its evolution when transported by the dynamical system. Let be xed (initial point) and let be a point trajectory of the dynamical system given by: [START_REF] Delahaye | Spatial Evolutionary Modelling (chapter 7)[END_REF] Assume now that trajectory is disturbed by a small perturbation we have : [START_REF] Delahaye | Airspace sectoring by stochastic optimization[END_REF] where is the differential of the vector eld at that satis es: [START_REF] Delahaye | Air traffic complexity based on non linear dynamical systems[END_REF] The Lyapunov exponents are closely related to the singular values of the matrix and can be thought as local shear values for the dynamical system. When Lyapunov exponents are high, the trajectory of a point under the action of the dynamical system is very sensitive to initial conditions (or parameters on which the vector eld may depend), so that situation in the future is unpredictable. On the other hand, small values of the Lyapunov exponents mean that the future is highly predictable (expected to be comfortable for a controller). So, the Lyapunov exponent map determines the area where the underlying dynamical system is organized.

It identi es the places where the relative distances between aircraft do not change with time (low real value) and the ones where such distance change a lot (high real value).

Let us now describe the practical procedure for computing 3D complexity maps. First of all, the optimal dynamic div-curl approximation for the observed trajectories is computed, based on the de ning equations. That step requires a linear system solving. The second step computes the second derivatives matrix at each point of the grid for trajectory starting at . This is done by solving the differential equation 11 with a Runge-Kutta integrator. The complexity value at point is then obtained by averaging Lyapunov exponents over the time : [START_REF] Cheng | The Optimal Partitioning of Networks[END_REF] The result of this process is given by a 3D Lyapunov exponents map which has to be cut into sectors. More information about this metric may be found in [START_REF] Delahaye | Air traffic complexity based on non linear dynamical systems[END_REF]. The next section introduces the model which has been used to address this problem.

3.B.6-5

Problem Modeling

The current approach is based on an initial cutting of the airspace into elementary cells. This initial sectoring mayby done by the means of a 2D Voronoi diagram for which the class center are the plane projection of the 3D crossing point between air ows. Based on this initial Voronoi diagram a merging process remove the cells with too close class centers and the too small ones. After this merging step, one has a 2D mosaic on the projection of the 3D airspace (see Figure 5). The right sideof this gure shows the voronoi diagram of the 2D prejection of the 3D cells.

Figure 5. Initial Voronoi Mosaic

A set of altitude layers is then considered and each mosaic cell is then extended in the third dimension by duplication of the 2D cell in each layer (see Figure 6).

Figure 6. 3D Extension

For each cell and for each layer, the associated workload is summarized by a weight located at the position of the cell class center with height given by the average altitude of the given layer. Those weights are symbolized by black circles on the Figure 6. A Delaunay triangulation of this initial mosaic is built in order to create a graph for which the link value represents the ow going from one cell into one of its neighboring cells. This initial Delaunay triangulation is extended to each altitude layer. This information will be used in the ow cut minimization process. Based on this initial mosaic, the following sector building process has been used in our application. First a set of sector centers is de ned: S. Those sector centers are points with coordinates in the 2D mosaic plane. The number of sector centers is K which is the number of sectors that has to be synthesized (|S| = K). Then, K layer intervals are de ned in a way that ensure that each altitude layer is assigned to an interval or more. Such association is called a covering. For instance, if 5 layers are considered (1,2,3,4, 5) a possible covering of this5 layers with 3 intervals is given by :I 1 = [START_REF] Lai Andw | Anevolutionary algorithm for the rectangular cutting stock problem[END_REF][START_REF] Kampke | Simulated annealing: Use of a new tool in bin packing[END_REF]; I 2 = [3,[START_REF] Dyckhoff | A typology of cutting and packing problems[END_REF]; I 3 = [START_REF] Cui | Generating optimal multisegment cutting patterns for circular blanks in the manufacturing of electric motors[END_REF][START_REF] Kampke | Simulated annealing: Use of a new tool in bin packing[END_REF]. As it can be noticed in this example all layers are assigned to one or more intervals. A cell association process to the sector centers is then applied. This process consists in taking each cell and nding the closest sector center for which the layer interval contains the cell layer. Based on the covering property of intervals, we are sure that, after this association process, all the cells are associated to one and only one sector center. In order to build a mathematical model, the following notation have to be introduced. We call C the set of elementary cells from the initial mosaic. A cell c i is characterized by its position where is the layer number of the cell and its weight . The cardinality of this set is given by the number of 2D Voronoi cells , times the number of layers , so . . We introduce K subsets each of them representing a potential sector (K is the number of sectors). A sectoring being a partition of the initial cell set, we have the following properties: The Delaunay triangulation is represented by a graph A for which the link symbolize the neighborhood relation between cells; this means thata link is built when two cells are neighbor. Each link is weighted by the amount of ow transiting from cell i into cell j. Based on those notations, one can formalize the objectives of this problem. The rst objective consists in balancing the weights of sectors. The weight w k of sector k is given by: (14)

The total amount of workload is given by: (15)

Sector weight are balanced when . So the first criterion is modeled by: (16) This criterion may vary from 0, when sectors are full balanced, to (K 1) when all the cells are into one sector and all the others sectors are empty. The denominator induces to be independent of the total weight. The second objective aims at minimizing the ow cut. Based on the Delaunay graph A this second objective may be modeled the following way:

(17) where

. Those two objective are gathered together into a single objective f by the mean of a weighted sum: f = . f1+(1 ). f2. There is not relaxation of constraints in this objective because they are already taken into account by the construction method of sectors. As a matter of fact, the initial Voronoi mosaic ensure that elementary cells met the min stay constraint and the safety constraint. The aggregation of cells to their nearest sector centers ensures that the route convexity constraint is met by sectors. Finally, the use of a covering of layers by intervals ensure that the synthesized sector will have cylinder shapes with polygonal section. Having now a mathematical model of this problem, let have a look to the associated complexity. The problem to be solved can be divided into two separate parts corresponding to the two different goals:

1. Equilibrium of the different sector weights; 2. Minimization of the coordination workload. The second criterion is typically a discrete graph partitioning problem with topological constraints and is then NP-hard [START_REF] Cheng | The Optimal Partitioning of Networks[END_REF]. The rst criterion is a classical classi cation problem with a connectivity constraint. One must nd an optimal grouping among possibilities; where is the second Stirling number: where is the number of cells and K the number of sectors. This problem is NP-hard and stochastic optimization is a good candidate to address it. Moreover, this kind of problem may have several optimal solutions (or near optimal) due to the different possible symmetries in the topological space and one must be able to nd all of them because they have to be re ned by experts (it is impossible to know at this step which one is really the best). This last point makes us reject classical simulated re-annealing [START_REF] Aarts | Simulated Annealing and Boltzmann Machines[END_REF] optimization which updates only one state variable. On the other hand, Evolutionary Algorithms (EAs) maintains and improves a population of numerous state variables according to their tness and will be able to nd several optimal (or near-optimal) solutions. EAs appear to be relevant to solve this sectoring problem.

Evolutionary Algorithms

Evolutionary algorithms use techniques inspired by evolutionary biology such as inheritance, mutation, natural selection, and recombination (or crossover) to nd approximate solutions to optimization problems. An individual, or solution to the problem to be solved, is represented by a list of parameters, called 3.B.6-7 chromosome or genome. Initially several such individuals are randomly generated to form the rst initial population (POP(k) in Figure 7). As can be seen, The rst step consists in the selection of the best individuals from population POP(k). Afterward, recombination operators are applied in order to produce the POP(k+ 1) population.

Figure 7. Genetic Algorithm with Tournament Selection

Then each individual is evaluated, and a value of tness is returned by a tness function. This initial population undergo a selection process which identify the most adapted individual. The one which has been used in our experiments is a deterministic ( ,μ)-tournament selection. This selection begins by randomly selecting individuals from the current population (POP(k) and keep the μ bests(( > μ)). This two steps are repeated until a new intermediate population (POPi) is completed. Following selection, one of the three following operators is applied : nothing, crossover, and mutation. The associated probability of application are respectively (1 pc pm)), pc and pm. Crossover results in two new child chromosomes, which are added to the next generation population. The chromosomes of the parents are mixed during crossover. Mutation is a genetic operator used to maintain genetic diversity from one generation of a population of chromosomes to the next. The purpose of mutation in EAs is to allow the algorithm to avoid local minima by preventing the population of chromosomes from becoming too similar to each other, thus slowing or even stopping evolution. These processes ultimately result in the next generation population of chromosomes (POP(k+1) in Figure 7) that is different from the initial generation. This generational process is repeated until a termination condition has been reached. The next section presents the application of EA to our problem.

Application to Airspace

Based on the previous problem modeling, a data coding has to be developed in order to make run EA on this problem. The state space may be summarized by a set of points in a 2D space with an altitude layer interval. The initial chromosome of our EA has the following structure:

x 1 y 1 x 2 y 2 ... ... x i y i ... ... x K y K L(M 1 )
L(M2) ... ... ...

L(Mi)

... ... ...

L(MK 1)

Each column of this table represents a gene of the chromosome. A gene numbered i consists in position coordinates (x i ,y i ), a layer marker level L( Mi ) and two marker extension limits and . The sector center enable the building of the sectoring for a given layer by the mean of an aggregation process of the Voronoi cell centers to their nearest sector centers (see Figure 8).

Figure 8. Aggregation Process Example

The markers are used to easily build a covering of the altitude levels by several intervals. Consider that K intervals have to be built then K 1 marker have to be generated on the layer domain (on Figure 9, four intervals are built with three markers).

3.B.6-8 The nal interval labeled K, is given by:

(19)
The chromosome is initialized by randomly drawing 2D points in the underlying surface (ground) and makers in the altitude layer domain (those markers are then ordered in the increasing layer number). This process is then repeated for all individuals in the population. An example of chromosome coding is given on Figure 10. This example presents a mosaic with 5 initial Voronoi cells (dash lines) for which 2 sectors are built. The number of layers considered is 10 and the rst sector addresses the layer 0 to 7 and the second one addresses the layers 5 to 10. This has been produced by the marker M 1 located at the layer 4 with = 1 and = 2. The Figure 10 shows one of the3 common layers of the two sectors [START_REF] Dyckhoff | A typology of cutting and packing problems[END_REF][START_REF] Suterf | Fastevaluationofvector splinesin three dimensions[END_REF][START_REF] Falkenauer | Ahybrid grouping genetic algorithm for bin packing[END_REF]. The sector centers are represented by symbol G 1 and G 2 . For the layers 0 to 4 only the rst sector addresses those layers meaning that all the associated cells belong to this sector. On the same way cells of the layers 8 to 10 are associated to the second sector. First, two genes are randomly selected in each parent. The associated sector centers are then connected by a straight line. Sector centers are randomly moved on this line in order to create two children sector centers. The second step of the crossover operator exchanges markers between parents (see Figure 12). The markers of the children come from both parents. The parents are labeled P1 and P2, the children C1 and C2. First the marker of both parent are gathered together and ordered according their layer level. Then, the rst child C1 take one maker over two in the ordered list, and the child C2 takes the remaining ones

Figure 12. Marker Exchange between Parents

The second operator used for our EA is the mutation operator as shown in Figure 13. To mutate a chromosome, one sector center is randomly selected and is moved by adding noises to its geographical coordinates

Figure 13. Mutation Operator

This operator may be decomposed into two classes:

1) strong mutation. This rst operator takes an individual and randomly initializes the position of the sector centers or the makers. The choice between those two options is random. 2) medium mutation. First, the more unbalanced genes are identi ed. This identi cation, is used to put some bias in the drawing of the gene which will undergo the mutation. The extension of the selected marker are then changed by taking into account the sign on the unbalanced weight of their genes. If the sector weight is under the balanced weight, the vertical extension of the associated interval will be increased and it will be decreased if it is over the balanced weight. Those mutations are randomly selected. The tness used by our EA has been built by mixing both objective the following way:

where I is an individual and = 0.8 and = 0.2. This tness has to be maximized and reach its maximum when f1 = f2 = 0 fitness= 100. Having presented the structure of our Evolutionary Algorithm, the next section shows some results on large instance of the problem.

Results

The parameters of the algorithm have been adjusted in order to reach a good performances for realistic instance of the problem (instances with several hundred cells in the initial mosaic). In order to validate the method, one has to produce problems for which exact solutions are known. To reach this goal, one has to create 3.B.6-10 symmetries in the problem which are easy to investigate from the human being perception point of view but which have the same complexity from the computer point of view.

The rst test considers 500 cells which have been randomly thrown on a 10×10 square (see Figure 14). Each cell has been extended in the third dimension on 10 layers. For the rst layer, random weights has been attached to each cell and duplicated on all layers. This mean that the total weight on each layer is the same, meaning that a trivial balanced solution may be reached by cutting the overall cube into horizontal layers. On the rst layer random ows have been introduced between any cells and its neighbor. This distribution of ows has been also duplicated in the third dimension. The following tests have been run on an Intel P4 (1.4Ghz) computer. The Figure 15 shows the evolution of the tness with generations. The solid line represents the fitness of the best individual among the population and the dashed line shows the associated average fitness. The exact solution is reach for the generation 30 for which the tness reach the value 100. The Figure 16 gives the associated evolution of both objective for the best individual. As it can be seen on the gure both objective reach the value zero meaning that the sectors are full balanced without ow cut. A larger instance of this problem with 1000 cells has been tested by our algorithm. For this test, no symmetries have been introduced and random values have been used for all cell weights and for all ows between cells. The solutions produced by the algorithm are very good as it can noticed on Figures 17 and18. In this case the tness may not be able to reach 100 because for any partition there is remaining ow cut. The remaining relative unbalanced weight is 0.45% and the associated relative ow cut is 0.13%. Those tests have been tried on several runs and produced the same kind of results. An HMI has been developed for this application in order to investigate the physical shapes of sectors by air traf c management experts (the Figure14 has been produced by the mean of this interface).

Conclusion

This paper has presented a new airspace cutting method which synthesizes balanced sectors with minimum ow cut. This problem is known to be NP Hard and stochastic optimization approach has been developed in order to address real instance of airspace cutting. A mathematical model has been developed for which two objective functions has been de ned in order to quantify the sector balance and the ow cut objective. The 3D cylinder sector shapes ensure that the constraint are satis ed. A coding process and recombination operators have been dedicated to this problem in order to make run a classical evolutionary algorithm scheme. This algorithm has been tested on realistic instance of the problem for which symmetries has been introduced in order to identify exact solutions. The results show that our EA nd the exact solution. A second test with 1000 cells has been tried for which no symmetry has been included. The algorithm nd very good solution with the same convergence rate as in the symmetrical case showing the robustness of such method.

Figure 1

 1 Figure 1 shows an aircraft crossing two sectors S1 and S2. The upper left, shows the radar top view of sectors S1 and S2. The dash line represents the trajectory of an aircraft going from S1 into S2. The bottom left shows the side view of S1 and S2 with the aircraft descending from S1 into S2.As it is shown on this gure, the sector border is not perpendicular to the ground (xy). The controllers having only the top-view of sector, they have dif culties to identify the point of coordination when the aircraft is exiting S1 and entering S2. On the left of line "a" the aircraft is in sector S1, on the right of line "b" the aircraft is in sector S2 but between "a" and "b" it is dif cult to know when the aircraft will change of sector mainly if it is descending or climbing. The way to solve this problem is to ensure that sector boundaries are perpendicular to the ground in the third dimension as it is shown in the right part of Figure1. The bottom right shows the new structure of the vertical border which is now perpendicular to the ground. The upper right shows the sectors top view with the "c" line materializing the sector transition. When the aircraft crosses this line, the controllers are sure that it is changing of sector even in case of some uncertainty in the altitude.

Figure 4 .

 4 Figure 4. Delaunay Triangulation (The Dash Line Represents the Associated Vonny Diagram)

Figure 9 .

 9 Figure 9. Layer Interval Building Process by the Mean of Layer MarkersThe makers M 1 ,M 2 ,...,M K are ordered according to their layer number :L(M i ) (L(M 1 ) L(M 2 ) ... L(M K )). The rst interval is built by using minimum layer number (0) and L(M 1 )+ built by using minimum layer number (0) and L(M1)+ Ext1 sup . The second interval begins at L(M 1 ) and ends at L(M 2 )+ . More generally the interval numbered has the following structure:

Figure 10 .

 10 Figure 10. Example of Chromosome CodingThis coding process has been completed by operators. Those operators are used to increase the diversity of the population. The crossover operator begin by randomly drawing two parents in the population of individuals after selection. The crossover operates in two steps. The rst step consists in changing the position on one sector center in each parents (see Figure11). In this figure one sector center is selected in each chromosome (one circle and one triangle). Afterward, those sector centers are connected by a straight line and are randomly moved on this line to create two new sectors (triangles).

Figure 11 .

 11 Figure 11. Crossover Operator

Figure 14 .

 14 Figure 14. Initial Distribution of 500 Cells in a 10×10 Square with 3D Extension

Figure

  Figure 15. Fitness Evolution with Generations (500 Cells)

Figure

  Figure 17. Fitness Evolution with Generations (1000cells)

Table 1 . Parameters of the Algorithm

 1

	Number of generations	500
	Population size	500
	Probability of Crossover	0.4
	Probability of Mutation	0.2
	Elitism	Yes
	( , μ)	(5,2)

Email Addresses

Daniel Delahaye: delahaye@recherche.enac.fr

Stephane Puechmorel: puechmor@recherche.enac.fr 27th Digital Avionics Systems Conference October [26][27][28][29][30]2008


Nicolas Barnier
email: barnier@recherche.enac.fr

Cyril Allignol
email: allignol@tls.cena.fr

Deconfliction With Constraint Programming

Keywords: Slot Allocation, Conflict Resolution, Constraint Programming

Current European Air Traffic Control (ATC) system is far exceeded by the demand and the resulting delays are a financial and psychological burden for airlines and passengers. One of their main sources is the hourly capacity constraints, defined on each en-route ATC sector, but poorly representing the workload of controllers. Whereas previous works were mainly focused on optimizing the ground delay slot allocation process performed by the Central Flow Management Unit (CFMU) to meet these constraints, we propose to directly solve all conflicts occurring above a given flight level by ground delaying, while minimizing the maximal delay. We present a Constraint Programming (CP) model of this large scale combinatorial optimization problem and the results obtained by its implementation with the FaCiLe constraint library.

Introduction

In an already saturated European sky, the predicted growth of air traffic volume urges to improve Air Traffic Management (ATM) efficiency, as attested by the ACARE Strategic Agenda 2 [START_REF]Strategic research agenda[END_REF] and the European Single Sky programme SESAR. Current ATM optimization strategies, like reducing the size of sectors or the distance of separation (RVSM, P-RNAV), seem to have reach the structural limits of the system, while the automation of Air Traffic Control (ATC) has known few significant improvements over the last decades [START_REF] Garot | Failures in the automation of air traffic control[END_REF].

The Central Flow Management Unit (CFMU) in Brussels is in charge of reducing these congestion costs by, among other strategic or tactical measures, delaying departure slots for the flights involved in overloaded en-route sectors. The purpose of delay-ing is to respect the en-route capacity constraints provided by each ATC Centre (ATCC) as a number of aircraft per hour, according to their daily schedule. Former studies like [START_REF] Dalichampt | Innovative slot allocation (ISA)[END_REF][START_REF] Barnier | Slot allocation with constraint programming: Models and results[END_REF] aimed at optimizing this process over the greedy algorithm used at the CFMU. However, one of the limitations of this process is that the definition of sectors capacities (hourly rate of aircraft entering the sector) is poorly related to the complexity of the traffic with respect to the controllers workload, as assessed by [START_REF] Gianazza | Selection and evaluation of air traffic complexity metrics[END_REF].

Instead of trying to satisfy en-route capacity constraints, we propose to directly solve the potential conflicts occurring between any two intersecting trajectories with ground delaying. A single delay would be associated with each flight such that all potential conflicts occurring above a given flight level would be avoided. This very fine grain model would of course generates much larger constraints sets than the macroscopic (at the sector level) capacitated ones, but would garantee conflict-free trajectories all along the flight path... provided that aircraft were able to scrupulously follow their predicted route in the four dimensions.

Obviously, the latter hypothesis is far from being met nowadays, even if a lot of studies about ATC automation have demonstrated that the accuracy of Flight Management Systems is a crucial issue for future Air Traffic Flow Management (ATFM) and ATC systems, as advocated by [START_REF] Alliot | Atm: 20 ans d'effort et perspectives[END_REF]. Nevertheless, we believe our approach may reduce air traffic complexity by "deconflicting" it in advance, the remaining conflicts due to deviation from the flight plan or occurring in the lower airspace being taken care of by more standard ATC procedures.

Several optimization paradigms are being evaluated for this purpose, namely meta-heuristics, local search and Constraint Programming (CP). We will focus here on the CP approach as it offers to obtain proved bounds on the maximal delays needed to solve the conflicts, which can be used to draw conclusions on the feasibility of this kind of regula-tions. Moreover, CP is a technology of choice for implementing such preliminary work, as it allows to easily refine the problem by adding new constraints (e.g. connection constraints between flights using the same aircraft) and to experiment with various search strategies without changing the rest of the model.

In the following sections, we first briefly present ATC and ATFM, focusing on ground delaying. Then we describe our model of a conflict-free slot allocation, starting by the details of the conflict constraints generation and search strategy, which lead to the presentation of our first results on a few instances for which optimal solutions were found. We end with planned further works to enhance the approach before concluding. The overall system currently implemented to achieve this goal can be conceptually divided in several layers or filters with decreasing time horizon with respect to the flight date of the traffic concerned: We will focus in the following section on the kind of regulations performed by the CFMU by postponing the takeoff of aircraft.

Context

ATC and ATFM

Air

Ground Delays

As aircraft obviously cannot be paused while airborne whenever the traffic complexity becomes to high to be safely handled by a controller, one of the simplest way to leverage ATC workload is to postpone the takeoff of aircraft1 . This kind of measure is however quite unpopular among airlines, as it can be very costly and may propagate in terms of missed correspondances (see [Uni04]), so the delays should be minimized as much as possible.

Satisfying En-Route Sectors Capacity Constraints

The aim of CFMU regulations is to maintain the number of en-route aircraft entering a given subset of sectors below some bound over given time periods (usually one hour), according to the constraints declared by experts (FMP) in each ATCC for the day of traffic. The CFMU experts first identify the overloaded sectors and responsible flows with the PRE-DICT tool, then compute a slot allocation as ground delays for the involved flights with the CASA tool (cf. [START_REF]Basic CFMU Handbook -General & CFMU Systems[END_REF]).

CASA is able to take into account many operational constraints and updates to optimize its allocation process, but the algorithm used has greedy properties and thus cannot guarantee to find a correct solution (which satisfies all the constraints) or an optimal one. CP technology has been applied with good results to prove and optimize the allocation process with a relaxed model [START_REF] Dalichampt | Innovative slot allocation (ISA)[END_REF] or to smooth the resulting load profiles [START_REF] Barnier | Slot allocation with constraint programming: Models and results[END_REF] with a tighter model.

However, traffic complexity is very hard to define precisely, and sector capacities, expressed as a maximum number of aircraft entering the sector over a given time period, does not take into account many parameters relevant to accurately represent the performance of ATC. Observed actual capacities, as well as merging and splitting subset of sectors, symptomatically present very different profiles than the predicted ones.

To overcome this issue, recent works such as [FPA + 07] use a much more precise and complex workload CP model to dynamically balance the traffic over the sectors of an ATCC in the upper airspace. Other works, like [START_REF] Barnier | Application de la programmation par contraintes à des problèmes de gestion du trafic aérien[END_REF] uses CP technology as well to optimize the ATCC opening schedules to match the predicted traffic more closely, or even attempt to redesign airspace sectorisation with better balancing like [START_REF] Huy | Airspace sectorization by constraint programming[END_REF].

Solving the Conflicts

Within the SESAR programme, the trend currently is the design of conflict-free 4D tubes within crowded airspace (whereas separation could be delegated to aircraft in less dense areas). So instead of only respecting sectors capacities macroscopically, we propose to evaluate the cost of precisely solving all potential conflicts, only with ground delays, while minimizing the worst allocated delay to maintain equity among airlines.

This conflict-free model will of course yield much larger problem instances as all the conflicting trajectories intervals above a given flight level will be taken into account as constraints. The resulting problem is intrinsically disjonctive as for each potential conflict between two flights i and j, either i must precede j or j precede i at each pair of points concerned (see section 3.2).

It is expected that such regulations will result in too high a total amount of delay, as we will optimize the maximum delay only, minimizing the sum of delays being exponentially more search intensive. However, the use of CP technology will provide proved bounds that other optimization techniques like local search or meta-heuristics cannot produce.

Other approaches have been presented to solve conflicts in real-time, automating the task of controllers. Some of the most promising ones are centralized techniques that compute simple horizontal or vertical manoeuvres [START_REF] Granger | Détection et résolution de conflits aériens : modélisations et analyse[END_REF] and small speed adjustements as proposed by the ERASMUS project [START_REF] Archambault | Speed uncertainty and speed regulation in conflict detection and resolution in air traffic control[END_REF]. These solvers, based around a metaheuristic (Genetic Algorithm), can take uncertainties on ground and vertical speed into account and repeatedly compute solutions for a sliding time window.

Conflict-Free Slot Allocation

Conflict Detection

Our input data are provided by the CATS simulator [START_REF] Alliot | CATS: A complete air traffic simulator[END_REF] which takes all filed flight plans concerning the French airspace for a given day of traffic and the relevant airspace configuration (sectors, waypoints...), and outputs the corresponding 4D trajectories. Trajectories are sampled with a 15s time step, which is the largest interval to guarantee that at least two points of the trajectories of facing aircraft at the highest possible speed will be closer than one separation norm, i.e. even the shortest conflicts will be detected. Trajectories are then probed pairwise for potential conflicts, taking the maximal allowed delay into account. The separation norm is thus tested for each pair of points of the two probed trajectories (up to p = 900 points per trajectory for up to n = 9500 flights in O(n 2 p 2 )) as illustrated on figure 2 in the horizontal plane.

Though the maximal allowed delay can be seen as a parameter of the search algorithm only, it also affect the conflict detection. Actually, when the maximal allowed delay is increased, the size of the problem grows as well, as more and more flights tend to be in potential conflict. Ultimately, if a 24h delay would be allowed, the conflict detection could be done in 3D, regardless of time, as any two geometrically conflicting trajectories would generate a constraint. So, whenever a particular instance has been proved inconsistent, it has to be generated again with higher values of the maximal delay, which will capture later potential conflicts on the trajectories pairs.

Operationally, flights originating outside the Eurocontrol countries cannot be delayed and we will thus fix their delay variable to 0 in our constraint model, reducing the number of variables but tightening the constraints as well and offering less opportunity for optimization. Constraints corresponding to conflicts occurring between two such flights will of course be discarded as we cannot delay the flights to solve them. Such remaining conflicting cases would have to be taken care of by other ATC or ATFM techniques that we will not address in this study.

Constraints Model

Conflicts

To compute the constraints of our model, the trajectories are pairwise probed for couples of conflicting points. Given a flight i, we note the input data:

• {p k i } the chronologically ordered sequence of the 3D points of its trajectory;

• t k i the time at which the flight will be at point p k i , should it not be delayed. We define a set D of decision variables:

D = {δ i , ∀i ∈ [1, n]}
of finite domain [min_delay, max_delay] that represent the delay associated with each of the n flights, and we will describe our model using the following auxilliary variables:

• θ k i = t k i + δ i the date at which flight i will be at point p k i if it is delayed by δ i ; • d ij = δ jδ i the difference of the delays of flight j and i.

For any geometrically conflicting points p k i and p l j such that the separation norm is violated (d h being the distance in the horizontal plane and d v in the horizontal plane):

d h (p k i , p l j ) < 5 Nm and d v (p k i , p l j ) < 1000 ft
we must temporally ensure that:

θ k i = θ l j
which can be rewritten with the difference variables d ij :

d ij = t k i -t l j
Starting at the first such point p k i that conflicts with a point of flight j, we take into account the whole continuous segment of trajectory j conflicting with p k i : {p l j , ∀l ∈ [l k , l k+r ]} for some r, and we impose that:

d ij ∈ {t k i -t l j , ∀l ∈ [l k , l k+r ]} d ij ∈ [lb k , ub k ]
with lb k and ub k being respectively the lower and upper bound of the set of t k it l j . If the next point p k+1 i of the trajectory of flight i conflicts with a further segment of flight j, we will obtain another forbidden segment for d ij :

d ij ∈ [lb k+1 , ub k+1 ]
taking part in the same potential conflict. To ensure separation we must then impose:

d ij ∈ [min(lb k , lb k+1 ), max(lb k , lb k+1 )]
as the conflicting segments of flight j overlap.

So if we take into account all the successive points of flight i, starting at p k i , that conflict with overlapping segments of flight j, up to some last point p k+s i , with lb 1 = min{lb k+u , u ∈ [0, s]} and ub 1 = max{ub k+u , u ∈ [0, s]} being the overall lower and upper bounds of the corresponding forbidden intervals for d ij , we can define the first conflict between flights i and j:

d ij ∈ [lb 1 , ub 1 ]
Note that we take as parameters of the problem instance the minimal and maximal allowable delays δ i ∈ [min_delay, max_delay], therefore the domain of d ij = δ jδ i is [-w, w], with w = max_delay -min_delay. We simply discard the conflict whenever ub < -w or lb > w.

A pair of flights may conflicts several disjoint times over their entire trajectories (as illustrated on figure 3), so several such disjoint intervals may be forbidden for the difference of their delays. For two flights i and j conflicting σ times over their entire trajectories:

d ij ∈ [lb 1 , ub 1 ] ∪ • • • ∪ [lb σ , ub σ ]
or, rewritten as a disjonctive constraint over the decision variables:

(-w ≤ δ j -δ i < lb 1 ) ∨ (ub 1 < δ j -δ i < lb 2 ) ∨ • • • ∨ (ub σ-1 < δ j -δ i < lb σ ) ∨ (ub σ < δ j -δ i ≤ w)
Figure 3: Three potential conflicts between to flights. The gray scale corresponds to time along the trajectory, the lighter the later.

provided that lb 1 > -w and ub σ < w, otherwise the first or last part of the disjonction is removed.

Cost

The cost of a solution is then defined as:

cost = max{|δ i |, ∀i ∈ [1, n]}
Note that min_delay = 0 and max_delay > 0 for all the results presented in section 3.4 such that the absolute values were not needed for constraining the cost and such that w = max_delay. Nevertheless, small negative values of the minimal delay could be of interest in an operational context, as aircraft can be ready for takeoff a few minutes before their slot in some cases. We can notice as well that if the cost was defined as | max D -min D|, i.e. the range of the delay values, equivalent solutions would be obtained by translating the domain [min_delay, max_delay], all the constraints being defined over differences between the decision variables δ i .

Other Parameters

The takeoff and landing part of trajectories are truncated around airports within a given radius (usually 10 Nm) as the traffic is considered handled with specific procedures by the Terminal Area control services in these zones.

After the computation of the conflict constraints, the whole instance is scaled down to a more reasonable time step (e.g. 1 min) than the 15 s used during conflict detection, ensuring that the original forbidden intervals are strictly included in the scaled ones.

Moreover, the flight level of the detected conflicts can be filtered, for example to only take into account conflicts occurring within the upper airspace (from FL290 and above). The minimal and maximal altitude of each conflict is recorded during the detection stage and a conflict is discarded if it entirely occurs below or above the specified airspace slice.

We allow as well to filter the time interval during which the conflicts may occur, taking the time bounds of the allowable delay into account. Any conflict strictly occurring outside the given time interval is then discarded.

For conflict constraints with several exclusion intervals, we add an extra parameter min_gap to prevent using very small allowed interval between two disjoints but very close conflicts. Such solutions would be too sensitive to small perturbations of the original schedule to be operationally robust.

Eventually, all the flights that do not have any conflict with any other flight are withdrawn from the instance.

Search Strategy

The constraints of the problem are reminiscent of the disjonctive mutual exclusion constraints modelling scheduling problems. At a coarse grain, we could consider each conflicting area as a machine on which to process two tasks of different lengths (depending on the speed of the aircraft). Several conflicts along a trajectory could be seen as the ordered tasks of a given job, as in the Jobshop Scheduling Problem (JSP).

However, the comparison does not hold much further. First, the time intervals between any two conflict tasks of the same trajectory is fixed, as only one delay variable is associated with each flight (unlike the JSP where all tasks are only related with precedence constraints). Secondly, to consider a potential conflict in 3D only, as the transitive closure of the overlapping conflicting segments, with task lengths proportional to the time spent by the aircraft within the area, is misleading. In this setting, the conflict associated with two catching-up flights on the same route would be the entire trajectory, preventing them from being airborne at the same time! Obviously, our model is much more precise and allows two aircraft on the same route being only separated by 5 Nm. Third, the number of "conflict machines", if not quadratic in the number of "flight jobs" as it could ultimately grow for arbitrary instances, is quite huge anyway as shown on figure 5.

Nevertheless, the branching scheme of our search strategy is inspired by standard scheduling tech-niques, because the essentially disjonctive nature of our problem shares some issues with scheduling ones. Trying to start the search by labelling the delay variables δ i would be highly inefficient, because the constraints are expressed over the differences d ij . Much more filtering is obtained by feeding the propagation of the arithmetic constraints with new domain bounds for the d ij auxilliary variables.

Similarly to some scheduling branching schemes, where tasks performed on the same machine are ordered pairwise (either task A precedes task B or B precedes A), we either add the constraint d ij < lb or d ij > ub in the case of a single conflicting interval. If there are several holes in the domain of d ij , branching is repeated with the bounds of the remaining holes. The variable d ij with highest sparsity, i.e. the smallest ratio between the domain size and the difference of the domain bounds, is chosen first for branching.

When all conflicts are ordered and there is no more holes in the domain of the d ij , we start labelling the decision variables δ i with a standard dom/deg selection heuristic.

After the first solution is found, the branch and bound algorithm then proceeds by dichotomy on the cost domain to find the optimal solution with respect to minimization of the maximal allocated delay.

Results

We have implemented this CP model with the FaCiLe library [START_REF] Barnier | FaCiLe: a Functional Constraint Library[END_REF] and obtained the following results on various day of traffic of years 2006 and 2007, with up to 8 000 flights and 285 000 intersecting pairs of trajectories taken into account for the lowest minimal flight level we could solve. About 10% of the flights are non-European flights, so their delays are fixed to 0.

The resulting constraint graphs typically exhibit only one single large connected component of maximal degree that can be greater than 300, i.e. a single flight may conflict with more than 300 other flights. Large cliques can also be found, as large as involving more than 60 flights, which indicates the presence of very entangled and dense traffic areas. The hardness of the conflict constraints is quite unevenly distributed, with two peaks for very small and very big (w.r.t. the maximal delay) forbidden intervals.

We mainly tune the size of our instances by decreasing the minimal flight level for which the conflicts are taken into account, typically aiming at the upper airspace (above FL290) where most of the cruising traffic occurs. However, we were able to optimally solve one instance with 6 600 flights and 260 000 conflicting pairs, taking all conflicts into account for the whole day, regardless of the flight level. Figure 4 shows the number of flights of our instances as a fonction of the minimal FL. The plots present two parts, one very steep from the maximal FL to FL300, then flights add up at a slower rate in the lower airspace (plots are labelled with their date, e.g. "070123" for the day of traffic on the 23rd of january 2007). The number of conflicting pairs is not quite quadratic with the number of flights, as mentioned in section 3.2 and shown on figure 5, but is quite huge anyhow, almost reaching 300 000 for our largest instance. We were systematically able to obtain optimal solutions within affordable computation times for all the instances where the 4 GB memory of our Core 2 Duo at 2.4 GHz were not exhausted. As shown in figure 6, small instances are solved in a few seconds whereas the biggest ones could take almost half an hour, growing only quadratically with the number of flights. We plan to address larger instances on a computer with more memory. However, the cost of this conflicts-free slot allocation can be quite high for the busiest days (our worst case is 162 min above FL250), but may be more reasonable (around 60 min) for less crowded days. Figure 7 shows that the cost grows steadily for small instances at high minimal FL, but jumps as soon as we add the main flows of traffic around FL350. The optimal cost then seems to be stable for larger instances, triggered only by the flights added around FL350. The corresponding means exhibit of course a smoother behaviour, but as we do not try to optimize the sum of delays, it is a poor indicator for our algorithm (and so is not represented on our graphs).

Note that for september 2008, Eurocontrol reported a mean delay of 25-30 min per delayed flights and a percentage of about 40-55% delayed departures. Some of these graphs show that some allocated delay may be as high as 90 min, but no precise figures are given on the maximal delay, though. For one of the days of traffic (plots labelled "060709s" and "060709d"), we have also tested our model on direct routes, i.e. aircraft fly in straight line in the horizontal plane from origin to destination, disregarding the waypoints of their flight plan (which we call standard routes). Direct routes are the ideal trajectories for airlines, with respect to operational cost, but such a traffic would be hardly controllable for human operators and ATC would have to be fully automated in this context. However, they tend to generate constraint graphs with a much lower tightness, and it is interesting to observe that the bounds obtained with direct routes are much better than the standard ones (70 min vs 162 min above FL250 for the chosen day, with respectively 121 819 and 284 646 conflicting pairs). Flights following standard routes tend to be on closer trajectories, suitable for current ATC procedures to be efficient, but not using airspace to its full capacity.

Further Work

These first results are encouraging for the efficiency of CP technology applied to such large combinatorial optimization problems from ATFM. However, the use of other optimization techniques like Local Search (LS) or meta-heuristics could be needed to further optimize the total amount of delay generated. Several such solutions are being investigated, including an hybrid LS algorithm that uses CP to filter large neighborhoods.

Validation With CATS

We also plan to validate our model and algorithm by taking the generated solutions as input delays for the CATS simulator. If the model is correct, only the non-European conflicting pairs and the conflicts occurring below the chosen flight level should remain. CATS can also be tuned to add uncertainty on the ground and vertical speed of aircraft. The robustness of our solutions could then be checked, e.g. verifying that the min_gap parameter (see section 3.2) has an influence on the remaining number of conflicts in case of perturbations.

Sliding Time Windows

We've only addressed so far the resolution of conflicts within the French airspace alone. However, in a unified European ATC context, all conflicting traffic throughout the Eurocontrol countries should be taken into account. Such instances would comprise up to 30 000 flights per day.

To be able to address larger instances we plan to adapt our algorithm to repeatedly solve slices of the problem on a limited time window T w , then to keep only the earliest part of the solution over a smaller interval λ and to slide the resolution window by λ. Parameters T w and λ must be carefully chosen according to the computation time of the resolution. Similar approaches are used for dynamic conflict resolution in the CATS solver as mentioned in [START_REF] Granger | Optimal resolution of en conflicts[END_REF].

Combined Flight Level Allocation

To be able to address such large instances as aforementioned, while maintaining reasonable maximal delay figures, we also plan to combine our slot allocation algorithm with a prior flight level allocation, possibly using CP technology as described in [START_REF] Barnier | Graph coloring for air traffic flow management[END_REF]. This first step, computed to minimize horizontally conflicting flows by separating them vertically (trying as well to deviate as little as possible from requested FLs), is expected to deconflict the traffic in a substantial amount before time slot allocation. The optimal cost should then remain within much better bounds than with the raw traffic.

Conclusion

We have presented a new ground delaying approach to solve all potential conflicts occurring above a given flight level for a day of traffic in the French airspace. Rather than trying to respect sector capacity constraints, we model each possibly conflicting situations between any two aircraft and impose delays to keep them separated, with the hypothesis that aircraft can precisely follow their planned 4D trajectories.

The resulting problem size is huge, but our CP algorithm is able to reach optimal solutions for all conflicts occurring inside the upper airspace. The resulting maximal delay can be comparable to delays allocated by the CFMU, but for the busiest days, solving all conflicts by ground delaying can be far too costly.

We plan to address larger (European) instances and reduce the costs of solutions with different techniques like repeatedly solving the problem on a sliding time windows and combining our delay algorithm with a prior flight level allocation.

Uncertainties remain a major issue in our approach, and our solutions still have to be checked on the CATS simulator to estimate their behaviour (remaining conflicts) with various degrees of noise in vertical and ground speeds.

Figure 1 :

 1 Figure 1: Vertical and horizontal separation. Another aircraft cannot be inside the cylinder at the same time.

Figure 2 :

 2 Figure 2: Conflicting Points Detection

Figure 4 :

 4 Figure 4: Minimal flight level vs number of flights.

Figure 5 :

 5 Figure 5: Number of flights vs number of conflicts.

Figure 6 :

 6 Figure 6: Number of flights vs computation time (optimality proof) in seconds.

Figure 7 :

 7 Figure 7: Minimal flight level vs optimal cost.

Note that flights might be delayed for other reasons than enroute capacity violation, like bad weather or equipment failures.