Christophe Hurter
email: christophe.hurter@aviation-civile.gouv.fr.

Serge Roux
email: serge.roux@aviation-civile.gouv.fr.

Stéphane Conversy
email: conversy@enac.fr.

Jean-Paul P Imbert
email: jean-paul.imbert@aviation-civile.gouv.fr

 C Hurter

 S Roux

 S Conversy

Visualization of Uncertainty for Air Traffic Controllers

Keywords: Uncertainty visualization, Air Traffic Controller, InfoVis

The goal of Air Traffic Controllers (ATCo) is to maximize both safety and capacity, so as to accept all flights without compromising the life of the passengers or creating delays. Because air traffic is expected to double by 2030, new IT systems must be developed to manage this traffic increase while maintaining or improving safety. ATCo radar screen displays past, current, and future aircraft position that contains uncertainty without displaying it. Uncertainty is mentally integrated by ATCo when assessing the future aircraft position. In this paper, we give some information regarding the radar data uncertainty and we propose new designs to display future aircraft positions more accurately with their uncertainty.

INTRODUCTION

The goal of Air Traffic Controllers (ATCo) is to maximize both safety and capacity, so as to accept all flights without compromising the life of the passengers or creating delays. Because air traffic is expected to double by 2030, new IT systems must be developed to manage this traffic increase while maintaining or improving safety. In order to monitor actual traffic, Air Traffic Controllers (ATCo) use radar screens which display aircraft positions (Figure 1). The radar screen displays the past and current aircraft position with the comet design [START_REF] Hurter | Temporal data visualizations of Air Traffic Controllers[END_REF]. These radar screens display data from the IT system called STR (Système de Traitement Radar). This IT system merges radar stations data to update tracks associated to aircraft. The data merging process is complex and generates uncertainties due to radar and tracking accuracy. In reality, radar visualizations hide these uncertainties and display the aircraft locations as if they were the actual locations. In order to avoid collisions, proximity warnings, and to optimize aircraft flow, ATCo give order to aircraft. For instance, to keep a safe distance between two aircraft, ATCo can ask the pilot to steer 10 degrees left. The radar screen will then display the steering aircraft when the pilot has decided to turn. ATCo have very limited tool to assess the future aircraft position and its uncertainty. They use aircraft current speed (displayed in the first row of its label), its altitude (displayed in the last row of the label). They can also display the speed vector that shows the future aircraft position if it keeps the same speed and direction, but again without its uncertainty. They also can measure distances in Nautical mile (Nm) by selecting points with the mouse pointer on the radar screen. As previously explained, ATCo use the past, current, and future aircraft position that contains uncertainty without displaying it. Even if it is not displayed, ATCo estimate this uncertainty. Therefore we believe that uncertainty visualization (data externalization) is a tool that will help ATCo to improve their efficiency while maintaining their safety level. In this paper, we investigate new radar visualizations that display the data uncertainty through new designs. The remainder of the paper is as follows:

We first detail the uncertainty while computing current aircraft positions. Secondly we detail our new designs that show uncertainty. And finally, we conclude this short paper with a discussion on the usefulness of uncertainty visualization for ATCo.

UNCERTAINTY IN R AD AR P OS I TIO N C OMP U TATI ON

Given the limited space in this document, this section gives simplified information regarding the uncertainty of aircraft position computation. Nevertheless, we attempt to answer this question: what are the uncertainty margins in pixels and millimeters of the displayed aircraft on radar screens. In order to compute the following data, we used the standard radar screen of 1600*2560 pixels, 640mm wide, 400mm high, with the typical zoom ratio that displays one sector and its surroundings which cover the whole screen area. This zoom ratio allows the display of 22 pixels per Nm. Each track (which corresponds to one aircraft) is updated every 4 seconds. Considering the worst case scenario, when an aircraft flight at a maximum speed of 500 Kt (knots, 900km/h), the error will be 0.55 Nm (1 km), i.e. 12 pixels on the radar screen (3 mm). This error corresponds to the distance that one aircraft fly during 4 seconds. ATCo must maintain a minimum safe distance of 5 Nm between aircraft otherwise they have to report the situation to aviation authorities. This mandatory safety distance is 10 times greater than the uncertainty due to the screen refresh rate.

Uncertainty related to the radar tracking system

Information in this section is extracted from internal documentation regarding the evaluation of radar over France [1] [2]. This process happens regularly to monitor the French radar tracking system (STR). For the purpose of clarity, the following uncertainty computations are basic. This section only gives order of magnitude of these uncertainties.

Given that typically 0.1% of displayed positions have an error greater than 0.15 Nm, this error represents 3 pixels on the radar screen.

Aircraft speed uncertainty

Typically 0.1 % of speed plots have an error > 20 Kt. The aircraft speed, which is displayed with the speed vector, shows at its end the future aircraft position in 3 minutes if the aircraft keeps the same heading and speed. This error corresponds to 1 Nm (3 minutes flight duration) on the radar screen, 22 pixels, 5 millimeters. The radar screen does not display the speed uncertainty.

Aircraft heading uncertainty

Typically 0.1% of heading plots have an error of 10 degrees when flying straight, which corresponds to 4 Nm (500Kt speed, 3 min flight duration), 88 pixels, 22 millimeters. When steering, 0.1% of heading plots have an error of 40 degrees which corresponds to 17 Nm (500Kt, 3 min), 376 pixels, 93 mm. This uncertainty is by far the most important one, but most aircraft are flying straight and the current aircraft heading is less important that its actual location.

Summary of radar screen uncertainty

The following table contains the summary of the uncertainty of the radar tracking systems and their impact on the ATCo's radar screen. This table shows that the uncertainties (for 0.1% of displayed positions) impact very little on the displayed aircraft location. Since the size of the symbol that depicts the current aircraft position is 14 pixels wide, it roughly represents the position uncertainty. This uncertainty is important regarding the extrapolated positions, especially when steering. But ATCo use the speed vector to assess minimum crossing distance only when aircraft are flying straight.

PROSPECTIVE DES IGN S FO R ATCO

In this section, we will detail new designs that display radar data with their uncertainties. As previously explained, positions uncertainty, is roughly embedded in the square that represents the current aircraft position. However, the uncertainty of aircraft extrapolated position (speed vector) is not visible. Therefore we propose new designs that display uncertainty of extrapolated positions, one with a curved speed vector, and one with a parrot scale.

Curved speed vector

The curved speed vector shows the user that the aircraft is actually turning. The animation speed gives him an indication of the turn rate.

A speed vector with a small curvature will show a slow turn. Furthermore, with this new design ATCo can see a more realistic future aircraft position compared to the straight speed vector.

Parrot scale design

This design displays future aircraft positions as dots. We connect dots of two aircraft's future positions if they correspond at the same time. We also added the minimum crossing distance (the yellow shape) and its location (the pink line). This design is animated; hence it shows how the crossing will evolve (Figure 3). This parrot design is impaired by the ladder that uses tamed parrot to climb in their cage. While this design is very accurate, it is very dense; it displays a lot of visual features, and is limited to two aircraft. But it enables ATCos to assess and forecast aircraft crossing.

PERSPECTIVES AN D CONCLUSION

As explained in this paper, uncertainty data is not displayed on the radar screen. Some data uncertainty is explicit, like the aircraft location that is embedded in the design of the current aircraft location. Others are mentally integrated by ATCo when assessing the future aircraft position. Radar screen design is a balance between accurate displayed information and cognitive workload to retrieve data. Complex designs with the visualization of uncertainty may spoil the radar screen and thus hinder the proximity detection and the safety level.

Figure 1 :

 1 Figure 1: radar comet 1.1 Uncertainty related to the screen refresh

Figure 2 :

 2 Figure 2: curved speed vector.

Figure 3 :

 3 Figure 3: parrot scale design.

We identified that ATCo need new tools to display future aircraft positions more accurately with their uncertainty. Therefore we have created new designs: curved speed vector, and parrot scale.

In the future, we plan to assess these designs with simulated traffic and ATCo in order to validate their usefulness.


Nour Dougui

Daniel Delahaye
email: delahaye@recerche.enac.fr

Stephane Puechmorel
email: puechmor@recerche.enac.fr

Marcel Mongeau
email: mongeau@math.univ-toulouse.fr

A new method for generating optimal conflict free 4D trajectory

The need for increasing air traffic capacity motivates 4D trajectory planning concept. In order to generate conflictfree 4D trajectories, we introduce a new concept based on light propagation modeling algorithm. This algorithm is a wavefront propagation method that yields a natural solution for the path planning problem specifically in the case of air traffic congestion.

I. INTRODUCTION

The analysis of air traffic growth expects a doubling of the number of flights over the next 20 years. The Air Traffic Management (ATM) will therefore have to absorb this additional burden and increase the airspace capacity, while ensuring at least equivalent standards of safety and interoperability. The European project SESAR was initiated to propose solutions to this problem. It relies on a new concept of air traffic control, known as 4D (3D plus time) trajectory planning, which consists in exploiting the possibilities of the Flight Management System (FMS) to ensure that a given aircraft is at a given position at a given time. For each flight, a reference trajectory, called Reference Business Trajectory (RBT), is requested by the operating airline. During the flight, conflict situations may nevertheless occur, in which two or several aircraft can dangerously approach each other. In this case, it is necessary to modify one or more trajectories to ensure that minimum separation standards (currently 5 Nm horizontally and 1000 ft vertically) are still satisfied. Moreover, it is desirable that proposed new trajectories deviate as little as possible from RBT. In this context, we propose a new algorithm which seeks to ensure sufficient separation between aircraft while producing flyable trajectories.

A. Previous related works

During recent years, several methods have been proposed to find an optimal solution that could solve conflicts in air traffic. The aim of these methods is to find for each aircraft, an optimal 4D trajectory that avoids conflicts with other aircraft, reaches the destination point and optimizes a cost function which depends on the travel duration and on the cost index (a coefficient that takes into account fuel consumption). There are mainly two classes of methods to address this problem: genetic algorithms [START_REF] Durand | Optimisation de trajectoires pour la rsolution de conflits en route[END_REF] and navigation-function based approach [START_REF] Roussos | Control of multiple non-holonomic air vehicles under wind uncertainty using model predictive control and decentralized navigation function[END_REF]. Each one provides only a partial solution to the problem.

The first one, genetic algorithms, consists in generating a new population of aircraft trajectories from a base population using three basic operators: selection, mutation and crossover in order to improve the cost function. This process is iterated until the cost function is no longer improved. The state space is a set of finite maneuvers, which are straight lines, turning points (changing an aircraft heading and then bringing it back on its initial trajectory) and offsets (inducing a lateral shift from the initial trajectory). Those maneuvers are the ones used by air traffic controller. Genetic algorithms generate trajectories with feasible operational maneuvers and with velocities within bounded ranges. They can reach asymptotically optimal solution, but for a given computing time, a feasible (conflictfree) solution is not guarantee.

The second method, based on navigation functions, consists in using an electrostatic modeling of the problem: an electron (which has a negative charge) is subject to an electric field, and is attracted by a positive charge which represents the goal and is pushed away by negative charges which represent obstacles. Thus, the electron is going to move towards the goal and steered by the resultant electric fields. The aircraft (a virtual electron) is represented by a point in 3D space.

If a mathematical function of potential fields can be built to model adequately the destination charge and the distribution of the obstacle charges, then the virtual forces applied on such virtual electron, initially positioned at the departure point, can be computed. This produces a trajectory which connects the departure point with the destination while avoiding obstacles (the other aircraft).

Navigation functions have already demonstrated their effectiveness in motion planning with guaranteed collision avoidance and convergence towards the goal configuration (reach the destination point with the right orientation). However, they do not take into account the constraints imposed by ATM, such as bounded speed, smooth trajectory and time constraints. Besides, they may tolerate large deviations from RBT.

The objective of our approach, based on an optical analogy, is to find for each aircraft a feasible (relevant to ATM constraints) optimal 4D trajectory, avoiding conflicts and which minimizes a criterion based on a local metric.

B. Paper overview

In the next section, we present our method. Numerical results are presented in section III. Then, a conclusion appears in section IV.

II. LIGHT MODELLING ALGORITHM

In order to build our algorithm, a light propagation analogy is used. In the physical framework, light propagates in space under Decarte laws (see Figure 1). These laws are summarized by the following: For a light ray that goes from a region with index n 1 into an other one with index n 2 with an angle i 1 , we have

n 1 sin(i 1 ) = n 2 sin(i 2 ) v = c n2
where, i 2 is the angle of the light ray in region with index n 2 , v its associated velocity and c the velocity of the light in vacuum.

Our algorithm uses such laws in order to build aircraft trajectories using congestion or other aircraft neighborhood as high index areas.

The light modelling algorithm is adjusted from the aircraft point of view. It is assumed that aircraft knows the surrounding aircraft trajectories (trajectories of other aircraft is a given input of our algorithm).

Assume that objective function is an application associating a positive real value to a curve of class C 1 of R 3 . Such value is computed by integrating a local metric along the curve. We can thus represent length, travel time or the cost associated with a trajectory by a suitable choice of local metric. Determining an optimal trajectory will therefore reduce to search a geodesic which is the shortest path between two points on the space regarding the local metric.

In this algorithm, we use the well-known fact that a light ray trajectory is a geodesic when considering the environment refractive index as a local metric. To represent congestion areas and conflicts in air traffic management, we consider the refractive index as a measure of congestion or traffic complexity. We select a barrier index value in the prohibited areas and in the protection volumes surrounding each aircraft. The optimal trajectory will be computed using a technique of ray tracing. The light will be slowed down in congested areas, but despite this, it can pass through. However, it will be completely blocked by aircraft protection volumes, which ensures conflict free-situations. We launch several light rays in various directions from the departure point of the aircraft. The path of the first ray that reaches the arrival point corresponds to an approximation of a geodesic.

We compute the environment index associated to a given congested area using a model based on Lyapunov exponents [START_REF] Puechmorel | Dynamical systems complexity with a view towards air traffic management applications[END_REF].

To generate a trajectory, we use a wavefront propagation algorithm in 3D with a space discretization (the wave propagation is done with a space step ds) from the departure point. We do not propagate the wavefront randomly in all space directions but into directions with highest probability of success. To ensure this, we guide the wavefront by an initial solution obtained by the navigation function method [START_REF] Roussos | Control of multiple non-holonomic air vehicles under wind uncertainty using model predictive control and decentralized navigation function[END_REF]. Consequently, we can guarantee at least one feasible solution.

To avoid a combinatorial explosion, the propagation will be coupled with a branch-and-bound algorithm that interrupts unnecessary shooting rays. In our case, the trajectory obtained by the navigation function method is sampled with halfspheres of radius ds, oriented towards the destination. These half-spheres are, in turn, sampled with an angle step dθ in the horizontal plane and an angle step dϕ in the vertical plane to build an initial search tree for the branch-and-bound algorithm as shown in Figure 2.

A. Branch-and-bound Algorithm

We implement the propagation within a branch-and-bound algorithm [START_REF] Balas | Branch and Bound Methods[END_REF], a classical framework for solving discrete optimization problems. At the beginning, we consider the set of all possible solutions, represented by the root. Procedures to obtain lower and upper bounds for the optimal value of our criterion are applied to the root. If these two bounds are equals, then the optimal solution is found, and the algorithm stops. Otherwise, the solution set is partitioned into two or more subproblems, which become children of the root. The method is then applied recursively on these sub-problems, generating a tree.

The idea behind the building of sub-problems which are relaxations of the original problem, branching process, is to solve them in a reasonable time. If an optimal solution is found for a sub-problem, it is feasible but not necessarily optimal for the original problem. On the other hand, as a feasible solution, it can be used to eliminate partial solutions. Indeed, if the lower bound of a node exceeds the value of an already known solution, then we can say that global optimum solution cannot belong to in the subset of solution represented by this node. It is therefore eliminated. The search goes on until all the nodes are explored or eliminated.

B. Branch and Bound applied to the Light Modeling Algorithm

The initial upper bound used in the Branch and Bound algorithm is the travel time computed by the navigation function method.

In order to build a lower bound for a given search tree node, we first compute a bound, "TimeToDest", for the remaining time to reach the destination. This bound is a weighted sum of two terms. The first one, "integTime", is the summation of the refractive index along the direct route to the destination. The second one, "maxSpeedTime", is the time needed to reach destination in direct route with the maximum speed.

T imeT oDest := α * integT ime + β * maxSpeedT ime. [START_REF] Durand | Optimisation de trajectoires pour la rsolution de conflits en route[END_REF] with weighting parameters α, β such that α + β = 1.

The lower bound is then the summation of TimeToDest and the time needed to reach the node from the origin (TimeToNode). More specifically, the lower bound is given by (see Figure 3): lowerBound := T imeT oN ode + T imeT oDest. Branching, in our context, involves launching rays as straight lines in a spatial half sphere of radius ds oriented towards the arrival point.

Browsing the search tree can be done in different ways. We choose a strategy whose priority is to find quickly a feasible solution (depth-first search or DFS). Here a live node (a node for which children have not yet been generated) with deepest level in the search tree is chosen for exploration. The memory requirement in terms of number of subproblems (stored at any given time) is bounded above by the number of levels in the search tree multiplied by the maximum number of children of any node, which is in our context a manageable number. The drawback of such approach is that nodes which are far from being optimal, may yield large amount of unnecessary bounding computations.

In order to avoid such drawback in our case, DFS is combined with a selection strategy. This consists in selecting the node that has the best lower bound among the nodes at the same level in the search tree (a combination of DFS as the overall principle and best first search, BeFS as a secondary selection criterion).

At each time step dt of the algorithm, the environment index is updated because aircraft change themselves the congestion while moving.

The main steps of the algorithm are as follows:

1.

Compute the navigation function algorithm solution: trajectory T. Set UpperBound := time travel of T.

2.

Descritize T to build a tree, with half spheres having radius ds, an angle step dθ horizontally and an angle step dϕ vertically.

3.

Set TrajSolution := null. While the destination is not reached do: a.

Update the refractive index. b.

While there is still unexplored nodes in the tree do:

• Choose a node N.

• Relaunch rays from node N: For any light ray, if the light beam goes from a region with index n 1 into a region with index n 2 with an angle i 1 , let it continue with a new angle i 2 such that n 1 sin(i 1 ) = n 2 sin(i 2 ) and with a velocity of v = c n2 . c.

Set TrajSolution := Trajectory obtained at b. interrupted by the step time dt. In the following section, we will see the numerical results of a simplified version of the algorithm (2D), first with a static refractive index and then with a dynamic refractive index.

III. NUMERICAL RESULTS

In all our simulations, we work on a 3.2 GHz machine running under Debian Linux operating system, 1024 KB of RAM. The software has been developed in JAVA.

We use a coordinate system that is scaled with separation standards. Thus, we use an (x, y) grid with a standard horizontal separation (5 Nm) unit. The index map used is a square of (15 * 15) standard horizontal separation.

In step 2. of the algorithm, the radius ds of the semicircles is set to a half standard separation distance. And the sampling angle dθ is set to Π 10 . The coefficients in the formula (1) we chose, are: β := 0.1 and α := 0.9.

A. Results in 2D

The algorithm was first tested with a static refractive index function (it does not depend on time) in 2D space in order to highlight the fact that it does find geodesics in simple cases.

Several refractive index functions were tested. Some examples are presented in Figures 4,5 Resulting trajectory with departure point on the bottom right corner and arrival point on the top left corner. The function index is given by: e -((x-a 1 ) 2 +(y-b 1 ) 2 )/k + e -((x-a 2 ) 2 +(y-b 2 ) 2 )/k + e -((x-a 3 ) 2 +(y-b 3 ) 2 )/k + e -((x-a 4 ) 2 +(y-b 4 ) 2 )/k Fig. 5.

Resulting trajectory with departure point on the bottom left corner and arrival point on the top right corner. The function index is given by: e

-((x-a 1 ) 2 +(y-b 1 ) 2 )/k + e -((x-a 2 ) 2 +(y-b 2 ) 2 )/k + e -((x-a 3 ) 2 +(y-b 3 ) 2 )/k + e -((x-a 4 ) 2 +(y-b 4 ) 2 )/k
In Figures 4, 5 and6, trajectories avoid high index area and passe through "valleys" which is the expected behavior. Thus, the aircraft avoids congested areas.

In Figure 7, the trajectory goes through a relatively congested area instead of bypassing it completely through the blue area above. This behavior can be explained by the fact that direct path, although it slows down the aircraft is more advantageous than a long detour with a higher aircraft velocity.

In all cases, the resulting trajectory is a geodesic approximation.

B. Results in 2D + time

Here, we deal with a situation of conflict resolution involving several aircraft. The algorithm controls one aircraft trajectory and we assume that the other aircraft maintain fixed straight line trajectories. The refractive index takes into account two factors. The first one is related to avoidance of the protection zones of other aircraft. The second one drive away the aircraft from congested areas where it could be in conflict.

The index function n is given by the following formula: Consider P aircraft (a i ) i∈{1..P } moving in the space with the associated position vectors ( -→ X i ) i∈{1..P } and velocity vectors ( -→ V i ) i∈{1..P } . For any space point

-→ Y , n( -→ Y ) := n 1 ( -→ Y ) + C( -→ Y ).
The first function n 1 is used to avoid conflict. In order to ensure that the aircraft controlled by the algorithm avoids the other aircraft, we represent them by disks whose radius is the standard distance separation. We set the function n 1 to a very high value N inside these disks and we make it decrease rapidly outside the disk as follows:

For any space point

-→ Y , let || -→ X i - -→ Y || := α.    α ≤ R ⇒ n 1 ( -→ Y ) := N N ≫ 1. α ≥ R ⇒ n 1 ( -→ Y ) := 1 + N -1 1+(α-R) q .
with R the standard distance separation and q is a parameter that determines the speed with which the index decreases outside the separation zone.

The second function C called convergence metric [START_REF] Delahaye | Air traffic complexity: towards intrinsic matrics, 3rd USA/Europe Air Traffic Management R& D Seminar Napoli[END_REF] models the fact that aircraft converge to a point in the space. It is a metric used to mesure congestion. The convergence value at a space point -→ Y is given by the formula:

C( -→ Y ) := i c(i) * e -β * (|| -→ Xi- - → Y ||) ,
where c(i) :=

j = i, r ij < 0 r ij * e -α * (|| -→ Xi- -→ Xj ||)
where

r ij := ( -→ X i - -→ X j ) || -→ X i - -→ X j || * ( -→ V i - -→ V j )
and where α and β are weighting parameters. Some examples of conflict resolutions with 2 to 4 aircraft are shown in Figures 8,9, 10, 11 and 12.

The track of the aircraft controlled by the algorithm appears in blue. This aircraft has its starting point at the bottom right corner of the figure and its arrival point at the top left corner. Its initial trajectory (a straight line) is shown in red. The other aircraft trajectories appear in black. The starting point of each trajectory is represented by a small dot.

In all these cases, the solution is found in less than 30 s of CPU.

In Figure 8, only two aircraft are in conflict. Their trajectories cross at the center of the figure. The controlled aircraft avoids conflict and follows a path that consists in two segments and an arc.

In Figures 9, 10 and 11, three aircraft are involved. The controlled aircraft is initially in conflict with one or the two other aircraft.

And each time, the trajectory generated by the algorithm avoids conflicts while remaining close enough to the original trajectory.

Moreover, these trajectories are composed of a set of segments and arcs.

In Figure 12, four planes are considered. As previously, the resulting trajectory avoids conflicts by using 3 segments and an arc.

Figure 13 represents a situation with four aircraft in conflict. The resulting trajectory successfully avoids conflicts, but the turns used are too steep. This occurs because the metric used for the index (convergence indicator) does not have a big enough predictability horizon. This creates paths that are not quite regular. We will therefore seek for a better congestion metric that avoids the occurence of such irregularities.

IV. CONCLUSION

The trajectory produced with our light-model algorithm avoids conflicts. It is a geodesic approximation that guarantee a speed lower bound, which is critical for an aircraft. Moreover, this trajectory is a sequence of segments and arcs, a trajectory that the FMS can monitor. Indeed, in most cases, it is a smooth curve which is flyable and which apparently does not require a lot of RTA points (Required Time on Arrival. This means that the plane must be in a given position at a given moment). Indeed, RTA points correspond to curvatures changes and in the tested examples, not a lot of curvatures changes appear. In the next step, a better convergence indicator will be used. The curvature of the resulting trajectories will be analyzed in order to check it satisfy aircraft constraints. This approach allows us to know how many RTA points are required to monitor conflict-free trajectories. Current work involves testing a new method for generating geodesic curves on a triangular mesh following ideas from computer graphics [START_REF] Novotni | Computing geodisic distances on triangular meshes[END_REF]. Our aim is to improve our results both in term of geodesic approximation and in term of computational time.

Fig. 1 .

 1 Fig. 1. Light ray deviation under Decarte laws in a region cut into two index areas (n1,n2)

Fig. 2 .

 2 Fig. 2. Initial search tree for the branch-and-bound algorithm

Fig. 3 .

 3 Fig. 3. The lower bound computing.

  , 6 and 7. The index function used is a continuous function, where high values are represented in red and low values in blue. Thus, the congested areas are represented in red and areas that involve little traffic are shown in blue. The resulting solution trajectory is found in less than 5 s of CPU.

Fig. 4 .

 4 Fig. 4.Resulting trajectory with departure point on the bottom right corner and arrival point on the top left corner. The function index is given by: e -((x-a 1 ) 2 +(y-b 1 ) 2 )/k + e -((x-a 2 ) 2 +(y-b 2 ) 2 )/k + e -((x-a 3 ) 2 +(y-b 3 ) 2 )/k + e -((x-a 4 ) 2 +(y-b 4 ) 2 )/k

Fig. 6 .

 6 Fig. 6. Idem Figure 5 with new parameters (a i , b i , ii ∈ {1..4}) for the index function.

Fig. 7 .

 7 Fig. 7. Idem Figure 5 with new parameters (a i , b i , ii ∈ {1..4}) for the index function.

Fig. 8 .

 8 Fig. 8. Conflict resolution with 2 aircraft. Controlled aircraft initial trajectory is in red with a departure point on the bottom right corner and an arrival point on the top left corner. Resulting trajectory is in blue. The other aircraft trajectory starts on the top in the middle and goes to the bottom. Its separation zone is represented by a red disk.

Fig. 9 .

 9 Fig. 9. Conflict resolution with 3 aircraft. Controlled aircraft initial trajectory is in red with a departure point on the bottom right corner and an arrival point on the top left corner. Resulting trajectory is in blue. The first conflicting aircraft trajectory starts on the top in the middle and goes to the bottom. The second conflicting aircraft trajectory starts on the left bottom corner and goes to the right top corner.

Fig. 10 .

 10 Fig. 10. Conflict resolution with 3 aircraft. Controlled aircraft initial trajectory is in red with a departure point on the bottom right corner and an arrival point on the top left corner. Resulting trajectory is in blue. The first conflicting aircraft trajectory starts on the right and goes to the left. The second conflicting aircraft trajectory starts on the right top corner and goes to the left bottom corner.

Fig. 11 .

 11 Fig. 11. Conflict resolution with 3 aircraft. Controlled aircraft initial trajectory is in red with a departure point on the bottom right corner and an arrival point on the top left corner. Resulting trajectory is in blue. The first conflicting aircraft trajectory starts on the left and goes to the left. The second conflicting aircraft trajectory starts on the left bottom corner and goes to the right top corner.

Fig. 12 .

 12 Fig. 12. Conflict resolution with 4 aircraft. Controlled aircraft initial trajectory is in red with a departure point on the bottom right corner and an arrival point on the top left corner. Resulting trajectory is in blue. The first conflicting aircraft trajectory starts on the left and goes to the left. The second conflicting aircraft trajectory starts on the left bottom corner and goes to the right top corner. The third conflicting aircraft trajectory starts from the bottom and goes to the top.

Fig. 13 .

 13 Fig. 13. Conflict resolution with 4 aircraft. Controlled aircraft initial trajectory is in red with a departure point on the bottom right corner and an arrival point on the top left corner. Resulting trajectory is in blue. The first conflicting aircraft trajectory starts on the left bottom corner and goes to the right top corner. The two other conflicting aircraft trajectory starts from the right and goes to the left.


Pascal Brisset

Gautier Hattenberger

Multi-UAV Control with the Paparazzi System

We present in this paper two experiments involving multi-UAV control. In the first one, three UAVs flew in formation flight, controlled by a single ground control station. In the second one, two aircraft flying in different places in the world have been controlled from the same GCS somewhere else. The two experiments have been conducted with Paparazzi, The Free Autopilot. Human interactions with the system were complex during these experiments. Controlling three aircraft in the same time in a safe way requires the help of a carefully studied interface. Controlling a UAV from a distant remote place requires to share the control authority between the local and the distant operators. We believe that the human work for UAVs control and command in this context is more similar to the work of an air-traffic controller than to the work of a pilot. The graphical interface of the Paparazzi system has been developed with this purpose in mind.

Introduction

A digest of the Paparazzi history can be given in a few sentences :

2003 A system is proposed to fly an autonomous MAV (Micro Air Vehicle). 2004 The system is improved to handle several aircraft (i.e. airframe configuration). 2005 The system is augmented to simultaneously control several aircraft. 2006 -2007 Constant evolutions aim to ease the control of multi-UAVs by a single operator.

Thanks to these effort to provide an ergonomic way of controlling UAVs, several results have been achieved:

-In November 2007, operators from the US army have been trained for the Paparazzi system in two days. Then, they were able to play the role of instructors for other operators who were able to control a fully autonomous flight, from takeoff to landing, after a couple of hours of training. -In February 2008, after a one week training by Martin Mueller 1 , a team from the Norwegian university of Bergen were able to fly an autonomous unmanned air vehicle equipped with meteorological sensors over the Spitsbergen island, in the Arctic circle, under very harsh conditions. -Up to five aircraft have been flown in the same time, controlled by a unique operator on a single standard laptop display.

We present in this paper two other experiments with the Paparazzi system where the operators had to interact in a complex way with a multi-UAVs. In the first one, the operator has to ensure safety while competing with a completely automatic flight planner which control a three aircraft formation flight. In the second one, several operators control several aircraft 1000km from each other.

These two experiments are described after a brief overview of the Paparazzi system.

The Paparazzi System

Paparazzi is an open-source autopilot system oriented toward inexpensive autonomous aircraft of all types. The project began in 2003 and has enjoyed constant growth and evolution ever since. The system has been used on dozens of airframes and implemented by several teams around the world. Hundreds of hours of autonomous flight have been successfully achieved with the Paparazzi system.

The Paparazzi system is extensively described in [START_REF] Brisset | PaparaDzIY: do-it-yourself UAV[END_REF][START_REF] Brisset | The Paparazzi solution[END_REF] and cooperatively documented in a wiki 2 .

This section focuses on the distributed architecture, the airborne subsystem, the navigation scripting language and the current available interactions with the navigation kernel.

Distributed Architecture

In order to offer maximum flexibility and openness, the system was designed from start as a distributed one (figure 1). Ground agents use the http protocol for data transfers and the Ivy software bus [START_REF] Buisson | Ivy : Un bus logiciel au service du développement de prototypes de systèmes interactifs[END_REF] for middleware. Ivy provides a basic message subscription mechanism and is available for most languages and architectures.

The system was designed to allow simultaneous operation of several vehicles and multiple ground stations.

The airborne vehicles may share a RF network or use a point to point connection to the ground station.

The "link" agent acts as a gateway between the RF side of the world and Ivy. It can speak different protocols and provides retry and acknowledgment services for hardware lacking those functionalities.

The "server" is a central agent that holds the configuration of the aircrafts (flight plan, settings etc) and other data (maps, topographic). It process the raw aircraft telemetry messages, log them for replay or analysis. It can redispatch synthetic asynchronous message to agents like ground control stations. The ground control station (GCS) allows for configuration, flight plan edition, simulation, flight control and log analysis. It typically runs in a laptop (figure 2).

"Simulators" are just regular ground agents and can be mixed with real vehicles or replay data.

The "Gaia" agent provides simulated data for environmental factors such as wind conditions and GPS coverage. The figure 3 shows a "Tiny" controller board, which at 20g and 70mm*40mm, is the most integrated hardware supported by Paparazzi.

Airborne subsystem

One of the unique features of the system is the use of infrared thermopiles for attitude sensing (figure 4 shows a typical airborne setup). For the airborne code, focus has been placed on security and reliability with the use of formal methods, code generation and critical code segregation [START_REF] Albert | Certification du code embarqué d'un micro-drone[END_REF].

Navigation Scripting

The behavior of the vehicle while in autonomous mode is described using a scripting language. This language defines stages and defines blocks as groups of stages. It provides a complete flow control mechanism (jump, conditional, loops, etc.) as well as expressions evaluation and exceptions.

Parameters used in a flight plan can be computed expressions. In the next example, UAV is asked to perform 5 circles at progressively increasing altitudes for exactly one minute at each altitude: <for var = "i" from = "1" to = "5"> <circle wp = "HOME" radius="75" alt = "ground_alt+50*$i" until = "stage_time>60" /> </for>

We provide an API to access the hardware, so that the control of peripherals such as a pan/tilt camera or landing gears/lights is possible from the flight plan.

Complex functions such as landing procedures can be factorized in libraries and called from the flight plan.

Being a Turing-complete language, it guarantees that any problem can be addressed and offers a very high flexibility. The tradeoff is that its complexity makes it poorly suited to real time interactions as we'll see in the next part.

Formation Flight Experiments

This section focuses on the management of a set of UAVs flying in formation. Formation flight can improve the safety and the efficiency of the UAVs as it allows mutual support in hostile environments and synchronized arrivals on targets. It is therefore particularly suited for military applications.

The goal of these experiments, conduct at the LAAS-CNRS, was to evaluate an approach that manages the formation to achieve mutual support along a pre-planned itinerary [START_REF] Hattenberger | Formation flight: Evaluation of autonomous configuration control algorithms[END_REF]. The configuration of the formation is autonomously selected as threats are encountered, and the UAVs trajectories that achieve the transition between configurations are dynamically planned and adapted. The algorithms constitute an intermediate layer between the mission planning system and the autopilot of each UAV: the planning layer is therefore released from "internal" formation problems.

After introducing the problem of our experiment, the experimental setup will be presented, and then some improvements will be proposed in the HMI and in the autopilot in order to reduce the operator workload.

Problem Statement

A planned mission is defined by a list of waypoints and a set of tactical constraints to satisfy, mainly consisting in inter-UAV distances and protection with respect to threats. Two kinds of threats are considered: early warning radars (EW), whose location is known at the time of mission planning, and missile tracking systems (TF), that are discovered during mission execution. EW can only be jammed by a specific device which is not available on all aircraft, whereas TF threats can be jammed with a device available on all aircraft, that blinds the tracking system in a given direction.

The configuration of the formation is defined by a set of slots corresponding to the relative positions of the aircraft. The formation layer is in charge of (i) the definition of a configuration according to the constraints and the environment, and (ii) the planning of reconfiguration trajectories for safe transitions between two configurations.

The position of the slots hosting the EW jammers is defined in a deterministic way, on the basis of geometrical parameters (the planned formation trajectory and the position of the threats to be jammed), and of an allocation that associate a jammer to a set of threats. To jam TF threats, a potential field approach is used to place the slots and provide mutual support in a reactive way.

Once a configuration is specified, a Branch&Bound algorithm is used to allocate the UAVs to the slots according to the constraints. As safe flight is a prime issue, reconfiguration maneuvers are planned to switch between two configurations, considering the inter-UAV security distance. An algorithm relying on A* is proposed for that purpose.

The Experimental Setup

For those experiments, we had 2 or 3 hobby aircraft (figure. 5) equipped with a Paparazzi board, a Paparazzi GCS, an Ivy bus and an external control architecture where the formation flight controller is implemented. This external system is designed to be embedded and so is fully distributed with local network communication capabilities. As the payload and CPU is very limited on our UAVs, the system is run on a ground laptop with an emulated network.

The control architecture of the UAVs has been developed based on the LAAS architecture [START_REF] Alami | An architecture for autonomy[END_REF]. Each aircraft is a functional layer composed of three modules. A first module is connected to the Paparazzi communication bus as shown figure 6. It gathers information about an aircraft's state and sends back the commands produced by the other modules (altitude, throttle and roll setpoints). The next module is dedicated to the navigation with a trajectory tracking algorithm and waypoint management. The last module is for the formation control (distributed) and configuration management (centralized). The formation module of each aircraft communicates with the others by an independent bus based on YARP [START_REF] Metta | YARP: yet another robot platform[END_REF]. A Mission Control Station (MCS) is used to send the mission data to the appropriate modules.

Achievement of the Experiments

The first steps of our experiments was to evaluate the control of the formation in a given configuration. A typical mission is achieved this way: the UAVs controlled by the Paparazzi GCS are placed on separated waiting circles. When they are ready, the operator sends them on parallel segments as synchronized as possible. In the meantime, the external control is started from the MCS and at the end of the segments the autopilot switch from the Paparazzi flight plan to the external system control.

At the first try, we had two UAVs with safety pilots and one operator with the GCS and MCS on the same laptop. Even with a wide screen, the two GUI cannot be shown on the same desktop, so the operator had to flip between virtual desktops. In those conditions, the workload is quite high and the operator missed to start the gateway that sends commands from the external system to the Ivy bus. As a result, the UAVs go stray to the ground with altitude setpoint to 0. One was taken back by the safety pilot, the other one was too low and land with no harm.

What clearly arise in this situation is that the operation of the aircraft from the GCS should be separated from the conduct of the experiment from the MCS. When the operator realizes on the MCS display that the altitude was wrong, he had to switch to the right desktop and then place the aircraft on their waiting circles. This last actions need at least one click to select the aircraft and a double click to change the current flight plan block. As a result, it took several seconds to place the aircraft in a safe state. The next experiment was conducted by two operators on two separated laptop connected over a network (point H on figure 6). With one operator on GCS, it was safer and easier to manage the UAVs during the initial and final phases. This way, we have been able to perform formation flights with two and three aircraft.

Reducing Operator Workload

This situation is not fully satisfying. It is pretty difficult for the operator of the GCS to perform the synchronized rendezvous during the initial phase for three aircraft. So, we can't have a better ratio than two operators for three UAVs. In order to improve this ratio, the operator's workload must be reduced. This has been done by improving the following points:

the current navigation block can be changed faster; -the UAVs can be automatically synchronized or they can start the mission already in formation flight; -the UAVs can perform basic collision avoidance.

The first point has been easily solved as the flight plan language of the Paparazzi system allows to set shortcuts buttons to go directly to a given block with one click. Moreover, the layout of the GCS can be configured to the needs of the experiments with multiple-UAVs. Figure 7 shows the strips of each aircraft clearly visible. A vertical view has also been added for this mission. Monitoring Fig. 7. On the top, the Paparazzi GCS. On the bottom, the Mission Control Center that allows to set and display virtual threats in the environment. Thanks to the Paparazzi HMI flexibility, safety has been improved by providing an efficient display to the safety operator of the mission.

the altitude of aircraft for air traffic controllers is always a problem. This solution is a good help as it shows all the planes in the same plot. However, from what we have seen, the operator can experience some difficulties to use it as the 2D map mainly catches his attention. Visual warning may be a good solution to make it really helpful.

A strong limit of the external system is that he cannot deal with the initial phase, i.e. placing the UAVs in formation from any initial situation. To overcome this problem, we have add basic formation flight capabilities directly in the Paparazzi system. One aircraft, declared as the leader, follows a standard trajectory (an oval in figure 7, on the right of the GCS). The other two try to fly keeping a given pattern. The main advantage is that the planes are automatically synchronized and in the right configuration. This way, the switch between the MCS and GCS is fairly easy. Once all the services are started on the MCS, the operator just need to click on a button on the leader strip and all UAVs switch their control at the same time. At the end of the mission, they can be placed again in formation or on waiting circles before landing.

The last point is a common problem to all multi-UAVs situation. How can we guarantee that an aircraft will not collide with an other? It is to much workload for an operator to ensure the lateral and vertical separation between UAVs if they fly in a reduced airspace and at the same level. Our solution is to implement a basic collision avoidance system based on the TCAS system used on civil aircraft. This system ensures a minimal vertical separation between two UAVs if they are in autonomous navigation. This way, safety pilots do not need to interrupt the experiment.

All these improvements allows a single operator to deal safely with two or three UAVs, even from two Control Stations. If the number of aircraft were to increase, two operators would be necessary, but the overall ratio between planes and operators would be better. This work is a start to prospect other situations of control shared between several GCS for a large number of UAVs, or for an air traffic control on top of the GCS.

Remote Control over the Internet

An experiment has been conducted in December 2007 with aircraft flying in Toulouse (France), Hildesheim (Germany), both of them controlled from the conference room of the 24C3 conference3 in Berlin (Germany) [START_REF] Mueller | Paparazzi -the free autopilot. build your own uav[END_REF]. Live video from the two aircraft were also pushed to a server providing real time streaming to any client around the world. While a standard DSL link was used in Toulouse, datalink and live video used a GSM connection in Hildesheim.

Architecture

The architecture of the experiment is displayed in figure 8. The different elements connected to the Internet network were :

-Two standard Paparazzi systems were set up in Hildesheim and in Toulouse, both involving

• A fully autonomous UAV connected to the ground station through two wireless connection, one for control and command and another one for video;

• A ground station organized around an Ivy Ethernet network;

• A graphical user interface and its operator. -A graphical user interface located in a conference room in Berlin.

-A server, located somewhere in Germany, responsible for relaying messages between the flight locations and the control interface in Berlin. This relay was convenient 1. to stream video images which were uploaded from Hildesheim and Toulouse; 2. to handle the datalink with the Hildesheim setup because the GSM connection was not symmetric and was not able to receive packets from any host, but only to receive data through an established connection.

When the UAV operators were ready on the flight fields, the operator in Berlin just clicked on the graphical interface to start the motor; the local operators just launched the UAVs. The operator in Berlin had a full control on the two missions, to change the flight altitude, change the navigation pattern ... or switch on and off the airborne video transmitter.

Human Interactions

In this experiment, the authority on UAV control was shared between the local operator and the remote operator. No priority had been defined between the two except that the local operator always had the option to stop (with one mouse click) the Ivy2TCP agent responsible of the connection between the local Paparazzi network and the Internet. The possibility of conflict between the two operators never occurred during the experiment or during the preparation.

Link

Conclusion

We have presented two experiments realized with the Paparazzi system. These experiments involve multi-UAVs control by several operators in the same time.

These results demonstrate that the design of Paparazzi provides the right architecture for distributed agents, vehicles or/and operators. The current Paparazzi graphical user interface, after years of experience and improvement, has been proved to be flexible, configurable enough, to be customized thanks to the operator needs. Hours of flight show that the system must be autonomous enough to give time to the operator to react to unexpected events; the UAV or MAV must autonomously stay in the air as long as possible while audible or visible alarms are shown to the operator.

Next step is to precisely handle authority sharing between different operators.

Fig. 1 .

 1 Fig. 1. The agents of the Paparazzi system. Aircraft, link and gcs agents may appear several times.

Fig. 2 .

 2 Fig. 2. The current ground control station

Fig. 3 .

 3 Fig. 3. The tiny controller board includes a GPS receiver and an antenna

Fig. 4 .

 4 Fig. 4. The Paparazzi airborne system components: (A)utopilot Control Board, (B)attery, (D)atalink Radio-Modem & Antenna, (G)PS Receiver, (I)R Sensors Board, (M)otor & Controller, (R)C Receiver & Antenna, (S)ervos, (P)ayload = Camera & Video Transmitter

Fig. 5 .

 5 Fig. 5. The three autonomous aircraft used at the LAAS-CNRS.

Fig. 6 .

 6 Fig. 6. The ground segment architecture with an external system connected to the Ivy communication bus.

Fig. 8 .

 8 Fig.8. Multi-UAVs control over a worldwide distributed network. Authority is shared between the local operators (in Toulouse and Hildesheim) and the remote operator (in Berlin). Note that the Relay Server allows several remote agents to interact with the UAVs.

Martin Mueller Engineering is a small UAV company which builds, operates and delivers a complete Paparazzi system.

www.recherche.enac.fr/paparazzi

http://events.ccc.de


P Kerlirzin

S Manchon

C Plusquellec

J B Gotteland

BUILDING AND EVALUATING A MINIMAL REGULATION SCHEME

The airspace that we consider is a volume without hole made up of sub-volumes. A capacity, described by a maximum throughput, is calculated in cooperation with ATC controllers several months before the day of operations. It is assigned for each sub-volume. A sub-volume is called a sector or a group of sectors. According to the predicted traffic and the available staff of controllers, sub-volumes number and capacity vary throughout the day and define a scheme for armed control positions, what we call the opening schedule. Depending on the opening schedule and on the predicted traffic sample, it is necessary to apply a set of regulation measures to avoid remaining sectors overloads. A regulation measure is described by a sector or a group of sectors, a time period and a capacity. A flight passes through 7 sectors or groups of sectors in France on average (13 in Europe). Hence, a single regulation measure may protect several sectors or groups of sectors. Knowing that the French control system manages up to 8000 flights per day (27 000 in Europe), it would be very effective to reduce as much as possible the number of treated regulation measures. In this document, we define a working methodology which aims at reducing the size of the regula-tion scheme. Simulations were carried out with SHAMAN experimental platform. Constraints programming was used to solve the problem.

CENA

CENA is in charge of studies related to air traffic management in order to support French and European Air Traffic Control (ATC) in an international cooperation. One domain of division RFM of CENA is to design and to implement mock-ups, prototypes and simulators for airspace and air traffic flow managers (ASM & ATFM respectively). SHAMAN is one of them.

Introduction Overview

Before the ATFM tactical phase, the FMP prepares its ACC's schedule (that we call opening schedule) for controlled sectors or groups of sectors, taking into account a predicted traffic sample and the available staff of controllers. The goal is to make this scheme as capacitive as possible and to define a set of regulation measures to protect every sector of the ACC.

The FMP then transmits the opening schedule and the temporary protection requests for airspace volumes that might be critical, to the CFMU executive unit (CEU). The opening schedule and the requested measures of protections are indissociable. This work is done locally in the ACC and does not take into account the capacities of the other ACCs' sectors. The CEU collects the opening schedule and the associated requests for protection coming from the different FMPs. It must make it consistent for all the ECAC zone. Considering the number of interdependent regulations and some other phenomena consecutive with exploitation of such a mosaic of capacities in the air route network, it is difficult to find a consistent scheme with appropriate protections for airspace volumes. Let us keep in mind that a flight passes through 7 en-route sectors on average in the French airspace and through 13 in Europe. Today, identification of the most penalizing volumes is made by a simple comparison between traffic demand and sectors capacities. This static identification highlights more sectors than necessary, i.e. regulation scheme may contain redundant regulations.

Goal

ATFM tactical process complexity (slot allocation and real time supervision) depends on the number of capacity constraints treated, knowing that they are interdependent. Decreasing the number of treated constraints allows: to simplify the process by suppressing redundant regulations: benefits for ATFM; to reduce total delay: benefits for companies.

The objective of this study is to address and validate a working methodology that allows flow managers to easily work out a minimal regulation scheme using a tool with an integrated take-off-slot allocation module. First, this tool will help the flow managers to identify the airspace volumes that take part in the slot allocation process. Then it will help them to extract only those which generate high delays during short periods or low delays during large periods.

Proposed method

We recommend the following working method to build the minimal regulation scheme :

Step 1 Identification of participating airspace volumes

The flow manager starts a slot allocation that treats every capacity constraint for all airspace volumes. This results in a list of airspace volumes, each of them beeing characterized by the individual delay it generated.

Step 2 Minimal regulation scheme elaboration The flow manager chooses a "participation threshold" i.e. a minimal value of the ratio airspace volume individual delay generated / total generated delay. The airspace volumes and periods during which they must be protected derive from it.

Step 3 Minimal regulation scheme validation

The flow manager starts a slot allocation process taking into consideration the only capacity constraints given by the minimal regulation scheme. Since all the capacity constraints are not treated, the remaining sectors overloads have to be quantified. If these traffic overloads are unacceptable in quality and/or in quantity (the cost function is presented in the "Protocol" of the technical validation part), step 2 must be reiterated with a smaller participation threshold. If not, the regulation scheme can be forwarded to the CEU.

Participation threshold definition

The aim is to find a relevant parameter that allows the flow manager to highlight critical airspace vol-3rd International Seminar on Air Traffic Management R&D CENA ✁ umes. In step 1, the slot allocation process treating all airspace volumes capacity constraints generates a total delay ✂✁ and the list of airspace volumes that generated delays. Let ☎✄ the delay generated by airspace volume ✆ during the regulated period. Among volumes that may need regulation ( ✝✄ ✟✞ ✠ ✡ ), some are characterized by a large delay genera- tion during a short time period and others by a short delay during a long time period.

Let us first define ☛ , the individual participation of ✆ in the total delay ☞✁ generated during the regulated period:

If ✌ ✄
is the total delay generated by ✍ , then ✎ ✑✏ ✓✒ ✕✔ ✒ ✕✖ where ✜ ✂✄ is the number of 30-minute capacity period relative to ✆ during the regulated period. For example, if ✆ is "armed" between 04h00 and 11h00, and then between 17h00 and 20h00, then the number of 30-minute capacity constraints (corresponding to the allocation step) is ✜ ✢✄ ✠ ✣ ✡ .

Moreover, if we decide to solve the problem between 4h00 and 22h00 (peak traffic), therefore

✜ ✤✁ ✠ ✦✥ ★✧ . Let ✚ ✗ the minimal participation delay: ✙ ✚ ✠ ☛ ☞✁ ✜ ✂✄ ✩ ☛ ✪✗ ✂✁ ✜ ✂✄ ✩ ☛ ✪✗ ✂✁ ✜ ✤✁ ✠ ☛ ✪✗ ✙ ✠ ✚ ✗
The minimal regulation scheme will be reduced to the airspace volumes which generate delays above minimal delay ✚ ✗ on a given 30-minute pe- riod.

In practice, the flow manager chooses ☛ ✘✗ ; the system uses ✚ ✗ :

☛ ✫✗ ☞✁ ✜ ✤✁ ✠ ☛ ✫✗ ✙ ✠ ✚ ✗ ☛ ✫✗ is called " participation threshold ".
We show that ✙ ✚ , therefore ✚ ✗ with constant ☛ ✪✗ doesn't increase significantly with ✜ ✬✁ where "capa width" is the width, in minutes, of a section of capacity.

Thus, the flow manager does not need to modify ☛ ✘✗ even if the allocation step is changed.

Simulations Data sample

The data sample we use relates to the French airspace.

Flight plans come from the French initial flight plan data processing system (STIP) archives. Regulation names and delays came from CFMU figures.

We measured traffic overloads (difference between traffic and capacity).

3rd

Slot allocation treating all capacity constraints

There are 1944 30-minute capacity constraints on average (

✂✁ ✚ ☎✄ ✚ ✝✆ ✟✞ ☎✄ ✠ ✣ ✱
). The total delay average is 55295 minutes The "=0" column contains the number of airspace volumes that do not generate delays.

( ✂✁ ✚ ✠✄ ✚ ✡✆ ✟✞ ☎✄ ✠ ☛ ✱ ✡ ✌☞ , ✍ ✏✎ ✒✑ ✠ ✓ ✣ ✔☛
The "0" column contains the number of airspace volumes which participated in strictly more than 0% and in less than 0.5% of the total delay ( -slot allocation with all capacity constraints; -definition of a reduced set of capacity constraints (minimal regulation scheme) with

✚ ✗
values based respectively on different ☛ ✘✗ values;

-slot allocation based on the respect of this reduced set of capacity constraints;

-study of remaining traffic overloads after slot allocation based on the respect of this reduced set of capacity constraints;

The following rates of remaining traffic overloads are considered: between 0% and 10% over the capacity: admissible regarding safety aspect; between 10% and 20% over the capacity: admissible regarding safety aspect only if they are ponctual. 20% over the capacity: inadmissible or non valid declared capacity for airspace volume. To check that a solution based on a minimal regulation scheme provides better protections of airspace volumes even if some traffic overloads remain (which is the case with CFMU), the following traffic overloads can be compared:

-initial traffic overloads (initial demand);

-remaining traffic overloads with allocated traffic after: a slot allocation respecting all capacity constraints applying a noise to the allocated departure hours;

-remaining traffic overloads with allocated traffic after: a slot allocation respecting capacity constraints of a minimal regulation scheme applying a noise to the allocated departure hours;

-remaining traffic overloads after operations.

Noise Adding a noise aims at simulating operational disturbances which affect allocated departure hours (bad meteorological conditions, delay at boarding, uncertainty on rolling time, ...).

A Gaussian noise with a mean delay of 5 minutes and a standard deviation of 16 minutes is applied to allocated departure hours. These parameters were provided by a statistical study based on STIP and CFMU inputs. Standard deviation is consistent with the window allowed by CFMU to take off (window of ✘ ✁ ✱ ✡ ✍ ✏✎ ✒✑ around each departure hour). In fact, departure hours variation is more complex: there is a main traffic peak at ✘ ✍ ✏✎ ✑ , a secondary one at ✡ and a last one at

✁ ✱ ✡ ✍ ✏✎ ✑ .
Choosing a Gaussian variation of departure hours (i.e. adding to an initial entry time of flights in airspace volume a calculated variation) is practical: it provides an entering traffic curve in airspace volume that is close to the observed one. First, for each studied sample, the initial traffic demand is treated with the SHAMAN slot allocation module (all capacity constraints are respected or only minimal regulation scheme constraints are respected according to each simulation). Then 50 pseudo allocated samples are created by applying Gaussian noise. This number of experiments is sufficient according to the standard deviation on number of remaining traffic overloads (the standard deviation is less than ✱ ✡ ✳ of the average).

Results

Simulations are made with different ☛ ✗ values: 0%, 1%, 1.5%, 2%, 2.5%, 3%, 5%. Each point of a given curve corresponds to a number of traffic overloads (before noise) for a given ☛ ✫✗ . For example, for ✚ ✗ corresponding to ☛ ✪✗ ✠ ✱ ✳ (noted ✚ ✗ ✮✭ ✯ ✰☛ ✫✗ ✠ ✲✱ ✘✳ ), only 18% of constraints need to be treated; total delay is reduced to 16%; the number of remaining traffic overloads is not very high: only about 100 flights.

Points on the figures (not on the curves) represent the remaining traffic overloads after slot allocation respecting all capacity constraints, followed by a Gaussian noise ( ✡ ✳ and ✱ ✳ ) applied to calculated departure hours.

The points are called ✡ and ✜ ✡ , ✱ ✡ and ✜ ✱ ✡ , and ✣ ✡ and ✜ ✣ ✡ when they represent respec- tively traffic overloads between ✡ ✳ and ✱ ✡ ✳ , ✱ ✡ ✳ and ✣ ✡ ✳ , and over ✣ ✡ ✳ , over the capacity.

First global quantitative results about con-

straints and delays show that for ✚ ✗ ✭ ✯ ☛ ✫✗ ✠ ✓✱ ✳ , after slot allocation and noise beeing applied, remaining traffic overloads are very small. The mean numbers of remaining traffic overloads are mentioned in the table 1 and in figures 1 and 2. They can be compared to the traffic overloads after slot allocation respecting all capacity constraints and noise beeing applied.

In all cases, traffic overloads are lower than remaining traffic overloads observed in the traffic after operation, especially concerning traffic overloads ✣ ✡ ✳ over the capacity.

Qualitative results about minimal regulation scheme show that the number of regulations of minimal regulation scheme is close to the number published by CFMU.

Morover, the regulations deduced by this technique correspond to the ones published by CFMU.

Case studies

Validation of the reference day

During the pretactical phase (2 days before the day of operations ), the flow manager knows 60% to 90% of the traffic. Since regulation problems are due to around 10% of the traffic, the pretactical preparation of the regulation scheme uses a predicted traffic. Generally, a good approximation of the traffic of day is the traffic of day ☛ which is called the reference day. The simulations we present now aim at verifying that this method for choosing the reference day is valid.

A slot allocation called ✒ ✒ is applied to the traffic of day . The set of capacity constraints that allocation has to respect derives from the regulations of the minimal regulation scheme, with: ✚ ✗ ✮✭ ✯ ✰☛ ✫✗ ✠ ✱ ✘✳ ; slot allocation with the traffic of day and the opening scheme of day .

A slot allocation called ✒ ✆☎ ✞✝ ✒ is applied to the traffic of day . The set of capacity constraints that allocation has to respect due to regulations induced by minimal regulation scheme when:

✚ ✗ ✮✭ ✯ ✰☛ ✫✗ ✠ ✱ ✘✳ ; slot allocation with the traffic of day First results show that more than 70% of regulations of the minimal regulation scheme ✒ ✆☎ ✞✝ ✒ are the same than the minimal regulation scheme ✒ ✒ ones.

( ✚ ✗ ✭ ✯ ☛ ✫✗ ✠ ✱ ✳ ).
The minimal regulation scheme ✒ ✆☎ ✞✝ ✒ offers a good protection of the sectors even if some variations concerning regulation periods are noticed. However, the number of remaining traffic overloads of 20% over the capacity, resulting from slot allocation ✒ ✆☎ ✞✝ ✒ is too big.

Validation of the local preparation of regulation measures

Each FMP prepares locally i.e. independently of others ACCs, its regulation measures. CFMU collects all regulation measures and merges them together trying to take into account the network effect between capacity constraints. This operational process is simulated: for each ACC i.e. for each ACC's opening schedule, a slot allocation is applied. Critical airspace volumes are identified and a local minimal regulation scheme is elaborated for each ACC. Then, local minimal regulation schemes are merged together in a minimal compiled regulation scheme. A slot allocation based on this compiled regulation scheme is performed to identify the regulation measures that generate no delay due to the network effect. Finally, these regulation measures are removed from the compiled minimal regulation scheme. This final scheme is called the compiled minimal regulation scheme. The compiled minimal regulation scheme is compared to the minimal regulation scheme. If they are close enough (number of regulations, names of regulations, periods of regulations...), then the proposed method of regulation measures preparation is applicable. An efficient coordination supposes that locally defined regulation schemes are changed as less as possible because of the consolidation achieved by CFMU. The common part between the two regulation schemes for each ACC is given in the following table: First results show that the percentage is rather high on average. It should not be penalizing that each FMP prepares its own set of local regulation measures if the method we present is applied.

Conclusion

The complexity of the ATFM tactical process (slot allocation and real time supervision) can be lowered, the total delay generated by the ATFM slot allocation decreased, and finally the safety increased (number of potential remaining traffic overloads reduction), by defining a minimal regulation scheme. This requires the use of a fast pretactical ATFM simulator.

The operational start-up of such a method should be easy because no modification of the actual working method is needed (local preparation of regulation measures and synthesis performed by the CEU).

  ) and so on for the other values of ☛ .It appears that on average, 35% of the airspace volumes do not participate to the ATFM process. For (☛ > 1.0%), only 19% of capacity constraints are treated.3rd International Seminar on Air Traffic Management R&D CENA✁Technical validation of minimal regulation scheme definitionPresentation of the simulation Protocol The simulation process goes through the following steps:

Figure 1 Figure 2

 12 Figure 1: Influence of

  of day . The results of these two simulations can be compared. The mean number of remaining traffic overloads resulting of slot allocation ✒ ✆☎ ✞✝ ✒ is compared to the mean number of remaining traffic overloads resulting of slot allocation ✒ ✒ . 3rd International Seminar on Air Traffic Management R&D CENA ✁ RFM Page 7

  ✔☛ ✡ for the 20/06, The average execution time of a slot allocation with 2000 constraints is about 1 min 10 sec. 20 sec are spent to generate constraints and 50 sec to find a solution, on a SUN ULTRA SPARC workstation using about 20 Mo on average. The following table gives, for each day of traffic studied, the number of airspace volumes that need protections according to ☛ , after a slot allocation treating all capacity constraints.

	✍ ✖✕ ✔✗ ✠ ✧ ✙✘ ✡ ✔✚	✥	for the 27/06) after a slot allocation
	that treats all capacity constraints.

Table 1 :

 1 Traffic overloads before and after slot allocation (CFMU, noise)

			after noised noised
		initial	ops	0%	✁ ✄✂
	Overloads	706 575	376	403
	Overloads capa+10%	355 272	118	134
	Overloads capa+20%	165 165	29	35

Table 2 :

 2 Overloads before and after slot allocation

		initial after ✒ ✆☎ ✞✝ ✒ ops.	✒ ✒
	Overloads	706 575	288 111
	Overloads capa+10%	355 272	118	10
	Overloads capa+20%	165 165	51	0
	The characteristics (volume, time period during
	which the regulation is valid, ...) of regulations of the minimal regulation scheme ✒ ✒ ✆☎ ✞✝ ✒ ( ✚ ✗ ✮✭ ✯ ✰☛ ✫✗ ✠ ✱ ✳ ) are compared to the characteristics of regula-
	tions induced by the minimal regulation scheme ✒ ✒

Table 3 :

 3 Regulations ✒ ✆☎ ✞✝

	✒
	✒ ✆☎ ✞✝ ✒ ✂ Aix Bord. Brest PE PW Reims ✁ 77 73 74 66 70 70 ✄ ✆☎ ✞✝ ✠✟ ✡✝ ☞☛ ✍✌ ✎✟ 7 7 16 17 15 11


Sonia Cafieri
email: sonia.cafieri@enac.fr

Aircraft conflict avoidance: a mixed-integer nonlinear optimization approach

Keywords: air traffic management, conflict avoidance, MINLP, modeling, global exact solution, heuristic 1

Detecting and solving aircraft conflicts, which occur when aircraft sharing the same airspace are too close to each other according to their predicted trajectories, is a crucial problem in Air Traffic Management. We focus on mixed-integer optimization models based on speed regulation. We first solve the problem to global optimality by means of an exact solver. The problem being very difficult to solve, we also propose a heuristic procedure where the problem is decomposed and it is locally solved by an exact solver. Computational results show that the proposed approach provides satisfactory results in reasonable time.

Introduction

The air traffic level currently attained in Europe is around tens of thousands of flights per day and it is expected to be multiplied by a factor of two during next 20 years. Air traffic is therefore at the core of the social and economic dynamism of our society. The European project SESAR gives the guidelines to go towards an Air Traffic Management (ATM) characterized by more efficiency and more safety, which should essentially result from a higher level of automation of ATM. The need for automatical tools to integrate human work specially arises in the context of aircraft conflicts detection and resolution.

Aircraft potential conflicts can be solved in different ways. The most commonly exploited is based on the idea of achieving separation changing the trajectory (heading angle) or the flight level of the aircraft involved in the conflict. Another way is based on the idea of separating aircraft by slightly changing their speeds but keeping the predicted trajectories. A speed regulation which occurs in a reasonable small range allows a subliminal control as suggested by the European ERASMUS project [START_REF] Cafieri | A mixed-integer optimization model for air traffic deconfliction[END_REF]. This project showed the advantage of such a control, which is not even perceived by air traffic controllers. Conflict avoidance is expected to be performed while deviating as little as possible from the original aircraft flight plan, minimizing the impact of the separation maneuvers. Various solution strategies have been proposed. A review is provided in [START_REF] Kuchar | A review of conflict detection and resolution modeling methods[END_REF]. Solution algorithms are currently mainly based on evolutionary computation [START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF]. These methods are computationally efficient, but the global optimal solution and even a feasible solution (with no conflicts) is not guaranteed to be achieved in a given time. Recent advances in mixed-integer linear and nonlinear programming open new perspective for modeling and efficiently solving the addressed problem. The first attempt to use mixed-integer optimization is by Pallottino et al. in [START_REF] Pallottino | Conflict resolution problems for air traffic management systems solved with mixed integer programming[END_REF], where, though under very stringent hypothesis, a geometrical construction leads to a mixed-integer linear programming model. More recently, mixed-integer programming has been proposed again for aircraft conflict resolution (see e.g. [START_REF] Martin-Campo | The collision avoidance problem: methods and algorithms[END_REF]).

In this paper, the speed regulation strategy is modeled by mixed-integer nonlinear programming, building on [START_REF] Cafieri | A mixed-integer optimization model for air traffic deconfliction[END_REF]. A deterministic global solution is first proposed, using a generalpurpose solver for MINLP. Then, to deal with the difficulty of the problem, another strategy is also proposed, where the optimality guarantee is forsaken in exchange for the computational efficiency. This solution strategy is based on hybridizing mathematical programming and a heuristic tailored on the problem.

The paper is organized as follows. In Sect. 2 we review mixed-integer nonlinear modeling of the aircraft conflict avoidance problem. In Sect. 3 we present a global exact solution of a few randomly generated instances and we propose a heuristic tailored on the problem to gain computational efficiency. Sect. 4 concludes the paper.

Modeling the aircraft conflict avoidance problem

Aircraft are said to be potentially in conflict when their horizontal or altitude distances are less than given standard separation distances (5NM and 1000 ft1 ). So, assuming the aircraft flying on a horizontal plane, the separation between aircraft i and j at the instant time t is expressed by the following condition:

||x r ij (t)|| ≥ d, (1)
where d is the minimum required separation distance and x r ij (t) is the vector of relative distance between i and j. We assume, as usually done, that speed changes occur instantaneously. We can therefore consider uniform motion laws. Hence, for each t, we have:

x r ij (t)=x rd ij + v r ij t, (2)
where x rd ij is the relative initial position of aircraft i and j and v r ij is their relative speed. Observing that the minimum t is given by t

m = -v r ij x rd ij /(v r ij )
2 the separation condition can be rewritten as follows:

(x rd ij ) 2 - (v r ij x rd ij ) 2 (v r ij ) 2 -d 2 ≥ 0. (3)
Note that condition (3) has to be checked only when the inner product v r ij x rd ij is negative. In this case, indeed, aircraft are converging.

In our model, conflict avoidance is achieved by performing a speed change maneuvre. Aircraft which are in conflict accelerate or decelerate in order to cross their conflict zone at different instant times, solving the conflict. Let A be the set of n aircraft. Decision variables are q k ,k∈ A, expressing the percentage of speed change of each aircraft with respect to its original speed. As prescribed by the ERASMUS project [START_REF] Bonini | How erasmus can support an increase in capacity in[END_REF], we impose bounds on these variables in order to have speed changes ranging between -6% and +3% of the original speeds. In this way, the so-called subliminal control is achieved. We minimize the speed change for each aircraft together with time intervals during which it flies with a modified speed, in order to deviate as less as possible from the original flight plan:

min k∈A q 2 k (t 2k -t 1k ) 2 , ( 4
)
where t 1k and t 2k are decision variables representing starting and ending instant times for aircraft k changing its speed. The order of instant times when aircraft change speed being unknown, 6 possible time configurations have to be considered for each pair of aircraft and, for each of them, 5 time intervals [t s ,t ′ s ]. The constraints of the problem impose aircraft separation (3) for each time configuration and time interval. This needs the introduction of new (integer) variables and constraints. See [START_REF] Cafieri | A mixed-integer optimization model for air traffic deconfliction[END_REF] for details. The described model can be relaxed, for example imposing that aircraft speed changes occur at the instant time t =0and that the new speeds are kept during the trajectories. We consider in the following this relaxed model. First, we have to express aircraft speeds in (3) in terms of their original speed v and speed modification q. We also have to check if t m is greater than 0. Equation for t m gives rise to a constraint, for each pair of aircraft, defining the minimum instant time. To check if t m ≥ 0, a binary variable y ij for each (i,j) is introduced (y ij =1if t m ≥ 0, and 0 otherwise), and constraints are adjoined accordingly. The separation condition is then imposed, for each pair of aircraft, only when t m ≥ 0:

∀i, j ∈ A, i = j, y ij x rd ij (v r ij ) 2 -(v r ij x rd ij ) 2 -(d 2 (v r ij ) 2 ) ≥ 0. (5)
The obtained mathematical programming model has as many (nonlinear) separation constraints as pairs of aircraft.

Solving the conflict avoidance problem

Global exact solution

We use as a testbed n aircraft in 2-dimensional space, placed on a circle of a given radius r, with speed v and a heading angle such that their trajectory is toward the center of the circle (or slightly deviated with respect to such direction). The zone of conflict is around the center of the circle where aircraft are placed, and each aircraft is in conflict with each other. It is easy to see that the number of conflicts is n(n -1)/2, so a large number of conflicts is generated in the same conflict zone. We solve the problem to global optimality using COUENNE [START_REF] Belotti | Branching and bounds tightening techniques for nonconvex MINLP[END_REF], which implements a spatial Branch-and-Bound based on convex relaxations. Results are reported in Table 1a (v=400 NM/h). They show that we are able to obtain global exact solutions up to n =6(i.e. 15 conflicts). However, an exact solution turns to be high memory and time demanding, even for a small number of aircraft, due to the high number of conflicts and the number of variables and constraints largely increasing with n. Hence, we are not able to solve the problem for n>6 even with the relaxed modeling. Objective function values show that aircraft separation is always achieved with very slight speed changes.

A heuristic based on local exact solutions

We then propose a heuristic procedure where we solve at global optimality subproblems involving up to 4 aircraft at a time, based on the observation that a solution can be efficiently computed for problems involving a small number of aircraft. Let a cluster be the transitive closing on conflicting pairs of aircraft (see, e.g., [START_REF] Granger | A traffic complexity approach through cluster analysis[END_REF]). The heuristic is based on the idea of decomposing the problem in subproblems (clusters) and solve the conflict avoidance problem on clusters. Let ncl be the number of clusters. At each step, ncl problems are sequentially solved by using an exact solver (COUENNE). Combining together all the results, in general the conflicts are not all solved because aircraft inside clusters are typically in conflict with aircraft inside other clusters too. After the resolution step on subproblems, the number of remaining conflicts is computed. If it is greater than 0, a new step is performed. To do so, the initial speed (which together with the initial position represents the data of the problem) of aircraft that are still in conflict is re-initialized taking into account the solution obtained at the last step. That is, if the (optimal) solution obtained for cluster i is such that an aircraft in this cluster has been accelerated with respect to its original speed, then its speed is modified by a random slight further increase. If it has been decelerated, then its speed is modified by a random slight further decrease. In this way, the information obtained at the previous step is preserved and the chances to keep aircraft separated inside clusters increase. To update the speeds, a local search is performed testing a number of candidates and choosing the one that minimizes the sum, over all conflicting aircraft, of the maximum violation of the separation constraints for each considered aircraft, divided by the number of remaning conflicts. When only one conflict is to be solved, this search is intensified to increase the chances to solve the problem. Aircraft speeds have to be bounded in the small range [-6%v, +3%v], so when speeds are modified these bounds have to be checked and speeds adjusted to fulfill this requirement. This may eventually lead to change the speed scenario provided by local solutions.

Results are reported in Table 1b. Values are averaged over 10 runs. For all test problems, all conflicts are solved. Comparing with global exact solutions (Table 1a), it appears that decomposing the problem does not significantly affect the quality of the result. In general, increasing n, faster solutions are obtained using a higher number of smaller subproblems. Solutions are obtained in reasonable time on problems involving up to 10 aircraft.

Conclusions

We presented an approach based on mixed-integer nonlinear optimization for the aircraft conflict avoidance problem. We are able to obtain global exact solutions for problems with up to 6 aircraft, while a new heuristic tailored on the problem and based on local exact solutions allow us to obtain good quality results even on problems involving many conflicts at a time.

Table 1a .

 1a Results obtained with COUENNE

	ID	n	r	obj	CPU time
					(sec.)
	pb n2 2 1 ×10 2 0.002531	0.15
	pb n3 3 2 ×10 2 0.001667	1.45
	pb n4 4 2 ×10 2 0.004009	12.87
	pb n5 5 3 ×10 2 0.003033	841.33
	pb n6 6 3 ×10 2 0.006033 51863.37
	n= number of aircraft
	r= radius of the circle (NM)
	obj= objective function value
	ncl= number of aircraft clusters

Table 1b .

 1b Results obtained with the proposed heuristic

	ID	n	r	ncl	obj	CPU time
						(sec.)
	pb n4	4	2 ×10 2	2	0.005151	26.97
	pb n5	5	3 ×10 2	2	0.004729	17.98
	pb n6	6	3 ×10 2	2	0.006402	17.33
	pb n6	6	3 ×10 2	3	0.007438	341.12
	pb n7	7	3 ×10 2	2	0.009215	131.34
	pb n7	7	3 ×10 2	3	0.008144	22.99
	pb n8	8	4 ×10 2	2	0.008220	759.40
	pb n8	8	4 ×10 2	3	0.007551	39.66
	pb n8	8	4 ×10 2	4	0.012034	48.99
	pb n9	9	4 ×10 2	3	0.009238	97.41
	pb n10 10 4 ×10 2	3	0.014047	484.49

1 NM (nautical mile)= 1852 m, 1 ft (feet) = 0.3048 m