A Bruyère

L Illoul

F Chinesta

S Clénet

Comparison Between NEM and FEM in 2-D Magnetostatics Using an Error Estimator

Keywords: Error estimation, finite element method, magnetostatics, natural element method

This communication deals with a comparison between two methods of discretization: the well known finite element method and the natural element method that is a meshless method. An error estimator, based on the nonverification of the constitutive law, is used. This estimation has been applied to two examples: a device with permanent magnets and a variable reluctance machine.

I. INTRODUCTION

N OWADAYS, the finite element method (FEM) is well es- tablished and allows solving many problems in applied physics. However, difficulties can appear with the FEM when the system includes moving parts. Actually, a remeshing step can be required to avoid a significant distortion of elements or its overlapping. Recently, novel methods have been proposed to circumvent, or at least alleviate, the main difficulties related to the mesh constraints. The most usual mesh constraints are related to the need of efficient remeshing procedures, which in the 3-D case constitutes a difficulty nowadays under active investigation, as well as to the field projection from the old mesh to the new mesh that induces an inevitable numerical diffusion whose unfavorable impact on the solution after numerous resmeshing steps is widely accepted. It was in this context that the so-called "meshless methods" were introduced one decade ago. Their applications in several domains experienced a spectacular progression in the last decade, after proving its ability for avoiding, or at least for alleviating significantly, the remeshing requirements. The natural element method [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] (from now on referred as NEM) belongs to this new and vast family of meshless methods and its interest compared with the vast majority of meshless techniques concern its interpolation property that allows enforcing essential boundary conditions in an easy way, as in the finite element method. We recall that the imposition of essential boundary conditions is the main issue in the other meshless methods. On the other hand, it has been proved that the solution accuracy does not depend on the quality of the subjacent triangulation ("meshless character"). Thus, we could affirm, crudely speaking, that the NEM is a kind of FEM in which the triangles could be distorted without a significant impact on the solution accuracy.

The NEM, which was firstly applied in the framework of solid and fluid mechanics [START_REF] Martinez | Updated Lagrangian free surface flow simulations with the natural neighbour Galerkin methods[END_REF], [START_REF] Alfaro | Recent advances in the meshless simulation of extrusion and other related forming processes[END_REF], has already been successfully tested in electromagnetic systems for accounting, in a simple manner, the existence of moving parts [START_REF] Illoul | Application of the natural element method to model moving electromagnetic device[END_REF]. However, in the electromagnetic framework the NEM and FEM solution accuracy comparison has never been addressed. This is a question needing a clear response and consequently justifying or disappointing the use of meshless approaches in this field. In this communication, we compare both methods NEM and FEM by using an error estimator based on the non verification of the constitutive law on two examples [START_REF] Synge | The method of the hypercircle in function-space for boundary-value problems[END_REF], [START_REF] Rikabi | An error based approach to complementary formulations of static field solutions[END_REF]. The first example concerns two permanent magnets in opposition and the second example a Variable Reluctance Machine.

II. PROBLEM DESCRIPTION AND ERROR ESTIMATOR

The magnetostatic problem, defined in the domain D, can be written as and

(1) and

(2)

is the magnetic flux density, is the magnetic field, and is the current density that is assumed to be known in the domain . is bounded by the surface where represents the outward normal unity vector. Moreover, and verify the constitutive relationship given by:

(3) with the coercitive magnetic field and the permeability. In the sequel, the problem is assumed to be linear. Let's call "admissible" fields and that verify (1) and (2) respectively. It is obvious that an admissible field couple verifying the constitutive relationship (3) is the exact solution denoted . We consider now an admissible couple and we define a constitutive law error 'e': [START_REF] Illoul | Application of the natural element method to model moving electromagnetic device[END_REF] with

It can be proven [START_REF] Synge | The method of the hypercircle in function-space for boundary-value problems[END_REF] that this error e is equal to the distance between an average field and the exact solution :

with (5) Therefore, the scalar e is a numerical error estimator for the admissible solution. To obtain the admissible field couple, several methods have been already proposed. The more accurate method consists in solving both potential formulations. The vector potential formulation ( and on ) yields an 2) and (3), can be solved by discretizing the associated Galerkin's weak form. For this purpose we need to define the approximation to be considered for both the trial and the test functions. This approximation is usually built from the values that the unknown fields have at nodes distributed on the domain D. A shape function ( is the position vector) is associated to each node . Thus, both potentials can be approximated from a linear combination of such shape functions [START_REF] Yvonnet | The meshless constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF] The only difference between FEM and NEM lies on the expression of the shape functions that in both methods are defined in a different manner. The ones related to the FEM are well known and therefore they do not need additional explanations. On the contrary, as the ones related to the NEM being less known, a brief introduction seems to be suitable. The NEM shape functions have been introduced by Sibson and are calculated using the Voronoï diagram [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. We are illustrating the concepts related to the construction of such shape functions by considering the simple distribution of five nodes depicted in Fig. 1. Some basic constructions, as the Delaunay triangulation and the associated Voronoi diagram, that serve to define the NEM shape functions in 2-D will be illustrated.

A. Voronoï Diagram

We define the Voronoï diagram associated to a set of nodes as the geometric subdivision of the plane into regions , called first order Voronoï cells, each one related to a distinct node , such that [START_REF] Misoum-Benziane | A new fully coupled two-scales modelling for mechanical problems involving microstructure[END_REF] with the coordinates of the node and the distance between node and point . Thus, is the region of the plane that contains the points closest to node than to any other node in . In Fig. 1 the nodes are numbered from "1" to "5" and the Voronoï cell related to node 5 is depicted. The construction of that Voronoï cell is quite simple. We start identifying the natural neighbors [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] of node 5. Now, we consider the segments joining node 5 to each neighbor node, and a straight line, normal to each one of these segments, traced at the central point of those segments as depicted in Fig. 1. These straight lines defined the Voronoï cell. Fig. 2 illustrates the construction of the Voronoï cells related to the other nodes. For additional details the interested reader can refers to [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] and the references therein.

B. Sibson Shape Functions

Once the Voronoï diagram is constructed, we can calculate the value of the NEM shape function , related to the node at the position . For this purpose we consider, in the scenario depicted in Fig. 3, the insertion of a point and the computation of the shape function related to node 5 at that position, that is, . We consider two possible placements of such Voronoï diagram is modified because the consideration of the cell related to the point just introduced that in our example consists of the polygon "abcd". The cell just inserted, the one related to point , intersects the former cells related to the original nodes. In Fig. 3 we emphasize in grey color the intersection between this cell and the one initially related to node 5. Now, the shape function related to node 5 at position is computed from the ratio between the intersected area (the grey area in Fig. 3) and the one related to new cell (polygon "abcd"), that is

" " " " (9)
It is easy to prove [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] that the support of the shape function related to a node is given by the union of the circles passing through this node and any other two neighbor nodes of such that no other node falls within the circle. Fig. 4(a) depicts those circles and the corresponding support for node 5. By computing , by using the procedure just described at different positions , the shape function can be represented. Fig. 4(b) depicts the shape function related to node 5. One of the reasons of the robustness of the Sibson interpolation against the distortion of the background mesh defining the Delaunay triangles and the associated Voronoï diagram, is that even if a triangle becomes very distorted, the support of the shape functions is not restricted to this triangle (as it is the case in the finite element method). Nowadays there is not a mathematical proof of the higher robustness of this kind of interpolation; however, there are numerous examples showing the low sensibility of NEM to the mesh distortion [START_REF] Yvonnet | The meshless constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF]. Several simulations in solid and fluid mechanics [START_REF] Martinez | Updated Lagrangian free surface flow simulations with the natural neighbour Galerkin methods[END_REF], [START_REF] Alfaro | Recent advances in the meshless simulation of extrusion and other related forming processes[END_REF] have been successfully performed using the same cloud of nodes throughout the whole simulation.

C. Nodal Sibson Shape Function Properties

The NEM shape functions are except at the positions of nodes where they are and on the boundary of the Delaunay circles where they are . The shape functions have a compact support delimited approximately by the neighbour nodes, leading to a sparse system. The band width of the resulting linear systems is a little bit higher than the bandwith associated to the FEM. That is precisely one of the reason of the natural element interpolation robustness. The shape functions have other interesting properties (10) where is the position of node and the Kronecker's delta. The so-called delta property (the first one in (10)) enables to enforce essential boundary conditions as in FEM. The other two properties in (10) prove the linear completeness, that is, the Sibson interpolation can represent exactly linear functions. Moreover, the partition of unity, second relation in (10) allows functional enrichments in the context of the partition of unity paradigm.

The computing time related to the construction of the Sibson shape functions is higher than the one associated with the definition of such functions in the finite element framework. However, nowadays, there are fast algorithms coming from the computational geometry community that are robust, fast even in 3-D, and some of them freely distributed. If we take the CPU time as comparison criterion and consider simulations involving massively remeshning operations (that define the natural context of application of meshless methods), finite elements simulations are nowadays faster than natural element method, because the existence of efficient and robust remeshing procedures in 2-D. If one consider convergence and accuracy criteria, the convergence rate for linear FEM and NEM is the same (first order in the energy norm) but as we proved in some of our former works in mechanics [START_REF] Misoum-Benziane | A new fully coupled two-scales modelling for mechanical problems involving microstructure[END_REF] (and we will prove also in magnetostatics later in this paper) the NEM accuracy is higher than the one related to the finite elements. Moreover, in models requiring field transfer after each remeshing process the superiority of the NEM is more relevant as argued previously.

IV. APPLICATIONS

First, two magnets face to face are considered (Fig. 5). Ten meshes with an increasing number of nodes were used. Fig. 5 gives the Voronoï diagram associated to the clouds consisting of 107, 800 and 2200 nodes. For the FEM simulation, linear interpolation on the Delaunay triangles was used.

For each mesh, both potential formulations have been solved using both the NEM and the FEM. Fig. 6 depicts the evolution of energy with the number of nodes. We can notice that the energies calculated with both formulations bound the exact energy. However, the energies computed from NEM seem to be more accurate that the ones coming from FEM calculations. By representing the evolution of the numerical error computed from (6) versus the number of nodes in a log-log scale (see Fig. 7) for both the finite element and the natural element solutions, we notice that both evolutions are linear and that the convergence rate, slope of such evolutions, is the same and approximately equal to one. We can also notice the higher accuracy of the NEM solution.

A more complex model in now considered. It concerns the variable reluctance machine (VRM) depicted in Fig. 8. The modelling of VRM is a very trick problem because of regions with sharp angles where the concentration of the magnetic flux density is very dramatic. Some simple error indicators are available [START_REF] Yvonnet | A simple error indicator for meshfree methods based on natural neighbors[END_REF] allowing the nodal refinement. The relative permeability of the rotor and stator yokes is equal to 1000. Different nodal distributions have been tested (the coarsest one containing 6000 nodes and the finest one more than one million) and both potential formulations solved by using the NEM and the FEM. The errors and have been computed from [START_REF] Synge | The method of the hypercircle in function-space for boundary-value problems[END_REF]. The error ratio evolves from a value close to 2 related to the coarsest nodal distribution, to 2.5 for the finest one.

V. CONCLUSION

In this paper, the accuracy of meshless natural element method has been compared to the FEM one. For the comparison purposes we proposed the use of an error estimator

Fig. 1 .

 1 Fig. 1. Construction of the Voronoï cell related to node 5.

Fig. 2 .

 2 Fig. 2. Construction of the Voronoï diagram.

Fig. 3 .

 3 Fig. 3. Shape function calculation related to the node 5, at 2 different positions x, T being grey area.

Fig. 4 .

 4 Fig. 4. (a) Voronoï diagram (heavy lines), Delaunay triangulation (thin lines) and the associated Delaunay circumcircles. (b) Sibson shape function.

Fig. 5 .

 5 Fig. 5. Studied system: two magnets surrounded by an air box.

Fig. 6 .

 6 Fig. 6. Evolution of the energy with the number of nodes.

Fig. 7 .

 7 Fig. 7. Global relative error estimator.

Fig. 8 .

 8 Fig. 8. Variable Reluctance Machine. The windings that have been supplied for the test is emphasized.

based on the non verification of constitutive equation. For the considered examples, superiority in terms of accuracy of the NEM solutions has been noticed.


Tambet Treimuth
email: treimuth@recherche.enac.fr

Daniel Delahaye
email: delahaye@recherche.enac.fr

Sandra Ulrich Ngueveu
email: ngueveu@laas.fr

A branch-and-price algorithm for Dynamic Sector Configuration

Keywords: Airspace Sectorization, Air Traffic Management, Branch and Price, Column Generation

Air traffic generates workload for the air traffic controllers in charge of the airspace. For a large airspace, a single air traffic controller is not able to manage all this workload and the airspace is divided into sectors, each of them being assigned to a controller. When the traffic demand is decreasing during the night, the sectors are gathered together into groups to reduce the number of controllers in operation. Nowadays, this regrouping is performed empirically by airspace experts. In this paper, we show how the branch-and-price method can be used to compute a balanced grouping of air traffic control sectors to optimally reduce the number of controller teams during daily low flow periods.

Introduction

Airspace is divided into sectors each of them being assigned to a controller team. For instance, the French airspace is divided into 89 sectors with different altitude extensions. Over the course of a day of normal traffic, the control workload fluctuates based on traffic demands between various origindestination pairings. Traffic demand is reduced at night inducing less workload in the airspace and therefore sectors are gathered together into groups, each group constituting a larger sector assigned to a controller. The objective of the algorithm developed is to adapt the airspace according to the dynamic of traffic. In order to reach this objective, we consider that airspace is divided into elementary sectors which have to be merged or split with time. Each sector having a workload, we are searching for groupings with balanced workloads between groups of sectors and this, for each time period. When an aircraft is transferred from one sector to another, there is some information exchanges between the controllers concerned and the pilots in order to ensure a safe transfer of the aircraft between sectors. This information exchange between controllers and pilots induce an additional coordination workload for controller which is linked to the flow cut by the sector borders, related to the number of flights that cross sector borders. In order to reduce the transfer workload, the algorithm will have to minimize the flow cuts as well. Sector design requires also that sectors have to be connected when belonging to the same group. In the following of this paper, the problem of finding dynamic grouping of airspace sectors will be called Dynamic Sector Configuration (DAC) and a sector configuration will be referred to as configuration.

Previous work

DAC has a number of conflicting objectives, related to either the static (coordination, sector control and transition workloads) or the dynamic aspects of the problem. Therefore, multi-objective approaches have been suggested as a method for the problem [START_REF] Brinton | Airspace partitioning using flight clustering and computational geometry[END_REF]. Some 2D DAC optimization approaches have been developed, such as constraint programming [START_REF] Trandac | Optimised Sectorisation of Airspace with Constraints[END_REF] and geometric algorithms [START_REF] Basu | Geometric Algorithms for Optimal Airspace Design and Air Traffic Controller Workload Balancing[END_REF], but these did not consider the multi-objective aspect of the problem. There exists also genetic programming based methods [START_REF] Delahaye | 3D airspace design by evolutionary computation[END_REF][START_REF] Chen | Dynamic airspace sectorization via improved genetic algorithm[END_REF] [START_REF] Xue | Airspace Sector Redesign Based on Voronoi Diagrams[END_REF] improved the Genetic Algorithm efficiency with an iterative algorithm. Initial works were limited to 2D design but recent works include the third dimension [START_REF] Sergeeva | Dynamic Airspace Configurationg Generated by Evolutionary Algorithms[END_REF]. [START_REF] Tang | A multi-objective evolutionary method for Dynamic Airspace Re-sectorization using sectors clipping and similarities[END_REF]. proposed an improved agent-based model (iABM) combined with GA to achieve the optimized sectorization. Authors of [START_REF] Li | A Spectral Clustering Based Algorithm for Dynamic Airspace Configuration[END_REF]) use a spectral clustering based algorithm and [START_REF] Kulkarni | Static sectorization approach to dynamic airspace configuration using approximate dynamic programming[END_REF] has proposed an algorithm based on Dynamic Programming.

Those previous work address partially the full problem and do not fully take into account the stability of the design between transitions. Indeed, air traffic controllers have to adapt to the changes which translate into extra workload. Therefore, transitions between configurations have to maintain some level of stability with time.

The remainder of the paper is organised as follows. The next section presents the DAC problem statement. Then a set partitioning-based mathematical formulation is introduced. The Branch and Price approach is described in the following section. Finally, computational results are presented before the Conclusion.

Problem statement

Let T={1,…,|T|} represent the discretized time horizon. The airspace is modelled with a timedependent weighted graph t G (V, E, W , t T) = ∈ where V is the set of vertices, E is the set of edges and W t is the set of edge weights and vertices weights at time period t T ∈ . Each vertex i V ∈ represents a (elementary) sector i of the airspace and its weight corresponds to the sector workload.

Each edge ( e (i, j) =

) represents the border between two neighbour sectors (i and j). Its weight corresponds to the coordination workload for the planes crossing the border if the two sectors belong to two different controllers. Our goal is to define a valid airspace configuration for each time period t T ∈ . A solution of DAC is therefore a graph partitioning of G for each time period t, where nodes belonging to the same partition at time period t correspond to sectors grouped together and assigned to the same controller at that time period.

We consider three main objectives to minimize: (a) the air traffic controllers' workload, (b) the workload difference between controllers, and (c) the changes between consecutive configurations. Therefore, the static cost of a solution is covered by the first two objectives and measures the quality of the sector grouping of each time period. The dynamic cost of a solution is covered by the third objective and measures the change rate between consecutive configurations. The user can therefore define his unique objective-function as a weighted sum of the three objectives ( a b c) α +β + γ , using γ to increase or decrease the relative importance of stability (the higher γ , the less changes will be favoured among consecutive sectorizations).

Set-partitioning-based mathematical formulation

The mathematical model proposed is based on the definition of the set of feasible configurations C t for each time period t. These sets are of course exponential in size, hence the Dantzig-Wolfe decomposition and column-generation-based solution method. Each configuration i is assigned a cost denoted c i t which reflects its static cost if it is applied at time period t. An edge e E ∈ is said to be frontier for configuration i if and only if its extremities belong to different controllers in the configuration i. An edge is said to be frontier at time period t if it is frontier in the configuration applied at time period t.

We introduce three sets of binary decision variables to model the DAC problem using an extended formulation: X i t , Y e t and Z e t . Variable X i t is equal to 1 if configuration t i C ∈ is applied at time period t and 0 otherwise. Variable Y e t is equal to 1 if edge e E ∈ is a frontier at time t and 0 otherwise.

Variable Z e t is equal to 1 if edge e E ∈ was not frontier at time t-1 but became frontier at time t, or if edge e was frontier at time t-1 and is no longer frontier at time t. This variable is equal to 0 if the edge status (frontier or not) has not changed from t-1 to t.

The resectorization problem can be formulated as follows:

CFB 1 1 (P ) min t T T t t t i i e t t e E i C f cX Z γ = =∈ ∈ = + ∑∑ ∑∑ (1) s.t. 1 . . } , { 1 t t i i C X t T ∈ = ∀∈ ∑ (2) ( ) 0, {1. . } , t t t e i i C e Y X t T e E ∈ + - = ∀∈ ∀ ∈ ∑ (3) 1 0 , { 2 . . } t t t e e e Z Y Y t T -≥ -+ ∀ ∈ (4) 1 0, . {2 . } t t t e e e Z Y Y t T - ≥ - + ∀ ∈ (5) {0,1} , {1.. }, t e Y t T eE ∈ ∀ ∈ ∀∈ (6) {0,1}, {1.. }, t t i X t T iC ∈ ∀ ∈ ∀ ∈ (7) {0,1} , {2.. }, t e Z t T eE ∈ ∀ ∈ ∀∈ (8)
where γ is a predefined dynamic/static weighting factor and C t (e) is the subset of configurations valid for time t that use edge e as a frontier.

The objective-function (1) minimizes the sum of the static and dynamic costs. The static cost measures the quality of the configuration with regards to the traffic at the time it is applied, whereas the dynamic cost measures the stability of configurations over time. The static cost is computed as a weighted sum of the workload difference among controllers and the total amount of coordination workload. Indeed some configurations might be balanced but generate a higher overall workload for coordination. The dynamic cost is proportional to the number of edges that change status between two consecutive time periods. Constraints (2) ensure that one feasible configuration is applied per time period. Constraints (3) link the status of each edge (frontier or not) with the set of configurations (C t (e)) in which the edge is a frontier. Constraints (4) and (5) ensure that the changes of status between consecutive instants are correctly computed for each edge. This model has a polynomial number of binary variables Y e t and Z e t , but an exponential number of binary variables X i t . Let (P') be the problem obtained by replacing constraints ( 7) and ( 8) with ( 9) and (10) has the same optimal solution as (P). Its validity is stated by Theorem 1.

0 ≤ X i t , ∀t∈{1..T},∀i∈C t (9) 0 ≤ Z e t ≤ 1, ∀t∈{2..T},∀e∈E (10)
Theorem 1: An optimal solution of MILP (P') provides an optimal solution for problem (P).

Proof outline: Although variables X i t are continuous, the model ensures that their values are always binary in feasible solutions, equal to 1 if configuration i is being applied at time t, and 0 otherwise. The same needs to be proven for variables Z e t . Proving that equations (4), ( 5), ( 6), (10) and objective-function (1) result into Z e t =0 if Y e t-1 =Y e t and Z e t =1 otherwise, is straightforward. Regarding X i t , let S be a feasible solution of (P'): Y∈[0,1] |C||T| , Z∈[0,1] |E||T| and X∈IR |L||T | . Given an instant t * ∈ T, let us denote:

• C >0 the subset of configurations used at instant t * . In other words C >0 ={i : X i t* > 0, ∀i ∈ C t * }. • E >0 the subset of edges that are frontier in at least one configuration i of C >0 .

• C >0 (e) the subset of configurations from C >0 that use edge e as a frontier. With a reasoning similar to the one applied for the extended formulation of (Ngueveu et al, 2016), it can be proven that C >0 (e) = C >0 , ∀e ∈ E >0 , from which it can be deduced that:

• either |C >0 |=1 and therefore all X i t* are binary

• or |C >0 |>1 but all configurations from C >0 are identical.

There are no identical configurations in C t of (P') therefore X i t are binary and (P') is equivalent to (P).

Branch-and-Price for DAC

Problem (P') has a polynomial number of constraints but an exponential number of continuous variables and some binary variables, therefore column generation can be applied to solve its linear relaxation and a branch-and-price solution method can be applied to solve (P'). Note that (P) could also be solved with a branch-and-price. However, allowing variables X i t to be continuous as it is done in (P') does ensure that the optimal solution can be obtained without needing to branch on the variables that are being generated by the pricing. As a consequence, the pricing procedure proposed remains valid at any node of the exploration tree. This means that the branching procedure can be left up to the mixed integer linear solver. Only the pricing procedure needs to be specifically implemented.

Column generation for DAC

A master problem (P r ) for DAC is a linear relaxation of (P'), where the set of valid configurations C t is replaced with a subset C tr ⊂ C t . Each master problem can therefore be solved with any linear programming solver. Let us consider the dual of (P') (resp. (P r )) denoted (D')(resp. (D r )). Let u t ,∀t ∈ {1..T} be the dual variables associated with constraints (12) and let v e t ,∀t ∈ {1..T},∀e ∈ E be the dual variables associated with constraints (13). Any variable X i t in (P') corresponds to a dual constraint in (D') expressed with equation ( 14).

( ) 14), where F(i) is the set of edges which are frontier in configuration i. Let t j C ∈ be a feasible configuration that does not belong to the restricted subset C tr , then the corresponding variable X j t does not exist in (P r ) and the corresponding dual constraint (14) j does not exist in (D r ). Let S Dr be the optimal solutions of (D r ). Only one of two different situations can happen: either constraint (14) j is violated by S Dr , or it is not. If the constraint is not violated, then adding it would not change the optimum dual solution, which means that adding primal variable X i t would not change the optimal solution of (P r ). There is therefore no need to add configuration j in C tr . If the constraint is violated, then adding it would change the optimal dual solution, meaning that the corresponding variable X i t should be added to the master problem and the configuration j should be added to C tr . Consequently, our pricing procedure searches for missing and violated constraints ( 14), which is equivalent to searching for configurations of negative reduced cost where the reduced cost of a configuration i is computed with expression ( )

t t t e i e F i u v c ∈ + ≤ ∑ (
t t t i e e E i c u v ∈ --∑ .
The major advantage of using an extended formulation for modelling and solving DAC is to let the model decide optimally what is the best selection and sequence of configurations, while focussing on generating for each time period the configurations of negative reduced cost. Note that this reduced cost takes into account both the static and the dynamic impact of the addition of the configuration at time period t, by including the values of c i t and v e t in its computations.

Pricing procedure for DAC: how to find configurations of negative reduced cost

Given a weighted graph H representing the airspace at a specific time period, there exists different graph partitioning algorithms that can be used to find the best configuration for that specific time period. Our pricing method modifies the graph weights before using the partitioning algorithm to ensure that the configurations generated are the ones with the best reduced costs. Recall that the reduced cost of a configuration is given by expression :

( ) t t t i e e E i c u v ∈ --∑ .
For a given time period t, u t is constant and has the same value for all configurations, therefore finding the configuration that has the best reduced cost at time t is equivalent to finding the configuration i with the best ( )

t t i e e F i c v
∈ -∑ value, and then substracting u t . Let α and β be the weights assigned to the controllers' workload difference ( t i δ ) and the edge costs in the computation of the static cost of a configuration, i.e.

( )

t t t i i e e F i c w α δ β ∈ = Δ + ∑ , then ( ) t t i e e F i c v ∈ -∑ = ( ) ( ) t t t i e e e F i e F i v w α δ β ∈ ∈ Δ + - ∑ ∑ = ( ) ( ) t t t e i e e F i v w α δ β β ∈ Δ + - ∑ .
As a consequence, finding the configuration that has the best reduced cost on the graph G with edge costs w e is equivalent to finding the configuration that has the minimum total cost on graph H with modified edge costs w e tr = w e t -v e t /β and then subtracting u t from that total cost. Therefore, it suffices to apply classic graph partitioning on the graph with modified edge costs w e tr and then subtracting u t from the total cost (on the modified graph) of the configuration obtained. If a negative value is obtained then a configuration of negative reduced cost has been found and should be added to the restricted set; otherwise, the current solution is optimal and the column generation stops.

Any airspace sectorization algorithm based on graph partitioning is applicable on the modified graph in the pricing procedure. This is another major advantage of the method proposed, because any additional restriction or regulation could be simply integrated in this pricing. Note that if the pricing procedure is heuristic, the resulting branch-and-price code is also heuristic, yet valid bounds and feasible solutions of good quality can be obtained. If the pricing procedure is an exact method, then the optimality of the final solution is guaranteed for the DAC.

Computational Results

The proposed algorithm was tested on instances from the Reims ATCC. The real radar data was used for those computations. In Reims, the sectorization can be modified every 15 minutes and the planning has to be done for the next 2 hours, therefore T=8 time periods. This computation was done at 6 different times: 8h00, 8h15, … 9h15, resulting into 6 time steps. For the subsequent computational evaluation, we used a Multi-Level heuristic [START_REF] Bichot | Graph partitioning[END_REF] in the pricing procedure, because of it is fast yet efficient. The solutions computed were compared with the configurations used in the real situation in terms of (a) load balance between controllers and (b) smoothness (stability) of the sectorizations. The results showed that the proposed algorithm produced solutions with a better workload balance [Figure 3]. It means that controllers have more similar workload. The solutions were also better from a smoothness viewpoint. This means that consecutive sectorizations were more similar with proposed algorithm than the sectorizations that were used in practice [Figure 4]. The initial master problem was created with |T|=8 initial columns (configurations), one per time period, generated by applying the pricing procedure on the original weighted graph for each time period with all dual variables values set to zero. Then the branch-and-price procedure was launched and new columns (configurations) were generated until no configuration of negative reduced cost could be found. The algorithm was run at each of the 6 times steps. [Table 1] shows the computing times with an implementation in C++ and with the framework SCIP 3.1.0 on an Intel Core i5-4210U CPU and 6 GB of RAM. Each master problem was solved with IBM CPLEX 12.6 whereas each pricing was done with the Multi-Level heuristic.

Conclusion

This paper presented a branch-and-price-based method for solving an airspace dynamic sector configuration problem. The method proposed is generic, adaptive and was fast on real world instances. It allows the user to define his unique objective-function by choosing the relative weight of the static versus the dynamic aspects of the problem: controller workload balance versus stability of the sectorization. It can therefore be used for an automatic planning of sectorization, to decide when to reconfigure the airspace and how often. Future work may include the automation of the parameter setting procedure that fine-tunes the values assigned to the different weighting coefficients in function of the user expressed preferences.

Figure 1

 1 Figure 1 Column generation scheme. Figure 2 Dynamic sector configuration algorithm.

Figure 3 :

 3 Figure 3: Quality of Load Balance.Figure4: Quality of Smoothness.

Figure 4 :

 4 Figure 3: Quality of Load Balance.Figure4: Quality of Smoothness.

Table 1 :

 1 Performance of the algorithm (SCIP, C++) at the different time

	Time-steps	1 (8h00) 2 (8h15) 3 (8h30) 4 (8h45) 5 (9h00) 6 (9h15)
	# configurations generated 55	48	40	35	72	45
	Computing time	175 s	155 s	130 s	119 s	201 s	145 s